4 курс / Лучевая диагностика / ТОМОГРАФИЧЕСКИЕ_ИЗМЕРИТЕЛЬНЫЕ_ИНФОРМАЦИОННЫЕ_СИСТЕМЫ

Еще один общий подход ослабления артефактов из-за конечного энергетического спектра излучения источника основан на использовании излучения с различными длинами волн.

Если измерения выполняются на нескольких длинах волн в разных узких энергетических полосах, что можно выполнить подбором соответствующих фильтров, то можно положить, что в преде-

лах каждой полосы μ(x, y) постоянно. Взяв логарифм от результа-

та каждого измерения и суммируя далее величины этих логарифмов, получим искомый линейный интеграл от усредненного

коэффициента поглощения:
ln	J0	+ ln	J0	+... = ∫μ1 (x, y)dl + ∫μ2 (x, y)dl .
	J1x
			J2x		L	L

Для двух энергий последнее выражение запишется:

2

∫

μ1 (x, y) + μ2

(x, y) dl

1

ln

J0

=

L

=

∫

μ

ср

(x, y)dl

2

J

2

J

1x

2x

L

или

J0

∫

μср (x, y)dl ,

ln

=

(J

J

2x

)1/2

1x

L

где J1x , J2x

– измеренная интенсивность выходного потока излу-

чения, соответственно, для двух энергий; μ1 (x, y) , μ2 (x, y) – ко-

эффициент линейного поглощения объекта исследования (однородного водяного фактора) соответственно для двух энергий; μср –

средний коэффициент линейного поглощения.

Для двух энергий этот метод называется методом «дуальных энергий». Он является достаточно «грубым» и требует применения селективных детекторов энергии, фильтрующих элементов для создания узких энергетических полос излучения либо селективных источников излучения по энергии. Все это накладывает большое ограничение для практического использования этого метода.

131

Более эффективный подход основан на использовании знаний о физическом механизме поглощения излучения.

Для обеспечения высокой точности реконструкции зависящей от энергии коэффициент ослабления можно разбить, как это показано в разделе 2.1, на функции, характеризующие фотоэлектрическое поглощение и комптоновское рассеяние

μ(E ) = aФ fФ (E) + aК fК (E) ,

(2.64)

где fФ (E ), fК (E) – зависящие от энергии функции вклада соответ-

ственно, фотоэлектрического и комптоновского поглощения; аф, ак

– независящие от энергии коэффициенты, характерные для объекта исследования (вещества).

Таким образом, выполняя два измерения для различных спектров источника излучения, как это следует из (2.56), мы использу-

ем два спектра N01 (E ) и N02 (E) . Отсюда получаем два уравне-

ния относительно двух неизвестных переменных, которые являются линейными интегралами от следующих коэффициентов:

AФ = ∫aФ dL ,	AК = ∫aK dL .	(2.65)
L	L

Решая нелинейные интегральные уравнения относительно АФ и Ак с помощью того или иного численного метода, можно по результатам измерений получить решение в виде ряда. После расчета значений АФ и Ак при всех углах (ракурсах) и положениях детектора можно получить свободную от артефактов реконструкцию, поскольку эти коэффициенты не зависят от энергии.

Помимо устранения спектральных артефактов этот подход обеспечивает так же раздельную реконструкцию распределений коэффициентов AФ (x, y) , AК (x, y) , что может быть использовано

для идентификации различных веществ. (Данный подход рассматривается достаточно подробно в гл. 3.)

Разработан также итерационный подход [2] уменьшения спектральных артефактов, в котором используются заданные пределы изменения величины поглощения в тканях человеческого тела. Здесь начальное приближение учитывает нелинейные спектраль-

132

ные ошибки в проекциях. Каждое значение интеграла реконструируется в предположении ограниченного диапазона изменения плотности вещества. Затем данные измерения проекций корректируются и реконструкция повторяется.

Другим подходом в уменьшении спектральных атефактов является нахождение по экспериментальным данным однородного водяного фантома корректирующего полинома проекционных данных [65]. Коррекция осуществляется путем подбора коэффициентов полинома соответствующей степени, в зависимости от требуемой точности, таким образом, чтобы «ликвидировать» провал μ в

середине фантома (см. рис. 2.9). (Этот подход, как и предыдущий, также подробно будут рассмотрены с точки зрения его практического применения в гл. 3.)

Нелинейности, обусловленные конечной шириной пучка из-

лучения. Помимо отсутствия спектральных сдвигов выражение (2.54) подразумевает также использование пучка излучения бесконечно малой ширины в направлении Z . Однако на практике размер пучка конечен, поскольку он должен передавать детектору конечную энергию. Для простоты будем полагать, что используется коллимированный пучок с распределением интенсивности излучения в

поперечном сечении S (x, z) , с помощью которого при угле θ = 0 получается проекция в направлении оси Y:

d t
P(l ) = ∫∫S (x, z) exp −∫μ(x, y; z)dy		dxdz ,	(2.66)
0 0	L
где d – ширина входного окна детектора; t		– ширина пучка в на-
правлении Z ; S (x, z)	принимает значения от 0 до 1.

Выражение (2.66) можно разложить в такой же степенной ряд, как и (2.59), следовательно, здесь будут иметь место такие же искажения, но зависящие не от E , а от Z . Эти искажения будут наибольшими, если поперечное сечение пучка излучения включает большие изменения линейного интеграла.

Предположим, например, что некоторой части поперечного сечения пучка соответствует костная ткань большой плотности. В

133

идеальном случае принимаемый сигнал должен соответствовать линейному среднему значению интегралов для частей поперечного сечения пучка, содержащих и не содержащих костную ткань. Однако вместо этого принимаемый сигнал будет характеризовать среднее значение экспонент от этих интегралов. Следовательно, область, не содержащая костной ткани, будет доминировать, что приведет к искажениям.

Для коррекции нелинейности этого типа в томографии сделано достаточно мало, если не считать применения очень узкого луча по ширине в направлении Z. Для уменьшения нелинейности этого типа можно использовать решетку очень малых детекторов в направлении Z. То есть, линейка блока детекторов должна содержать несколько детекторных слоев малой ширины в направлении Z.

Элементарная область Sm (x, y) , покрываемая каждым таким

малым детектором, будет достаточно небольшой, что обеспечит пренебрежимо малые изменения линейного интеграла. Следовательно, интенсивность излучения, измеряемая каждым малым детектором, будет равна

d tm
Jxm = ∫∫	Sm (x, z) exp −∫μ(
0 0		L

= J0m exp −∫μ(x, y, z)L

x, y, z)dy dx dz =

dy , (2.67)

где tm – ширина малого детектора в направлении Z;

d tm

J0m = ∫∫ Sm (x, z)dx dz – эффективная интенсивность излучения

0 0

источника, регистрируемая отдельным малым детектором. Функцию, описывающую источник J0m , можно вынести за знак

интеграла (2.67), если положить, что поглощение слабо меняется внутри области, покрываемой каждым парциальным пучком, приходящимся на малый детектор.

134

Искомый интеграл от усреднения линейного коэффициента поглощения будет в данном случае определяться выражением

		1	M		J0m
	∫μср (x, y, z)dy =		∑ ln				,	(2.68)
	L	M m=1			Jxm

где М – количество малых детекторов в направлении Z, охватываемых ширину пучка.

Рассмотрим далее класс искажений, связанных с недостаточностью проекционных данных.

Реконструкция основана на использовании полного набора линейных интегралов или проекционных данных. Однако реальная томографическая система в той или иной степени отличается от такой идеальной модели. Очевидным примером в томографии является использование дискретизированных данных как по положению лучей в пучке излучения за счет дискретных ячеек (единичных детекторов) в линейке детекторов, так и по положению пучка в угловом пространстве (пространстве ракурсов). Для дискретизированных данных возможно наличие эффекта наложения спектральных пространственных составляющих сигнала на детекторе.

Для удовлетворения условий теоремы отсчетов Котельникова [12] и устранения эффектов наложения дискретизация должна проводится с пространственной частотой, вдвое превышающей пространственную частоту наивысшей частотной составляющей сиг-

нала. Оценка непрерывной проекционной функции P(l,θ) для случая параллельной геометрии сканирования будет иметь вид

P	M	N
	(l,θ) = ∑∑P(l,θ)δ(l − m l )δ(θ− n Δθ) ,		(2.69)
	m=1n=1
где M , N –	соответственно количество дискретных отсчетов ра-
курсов и детекторов;		l, Δθ – интервалы дискретизации, соответ-

ственно вдоль проекции и по ракурсу.

В некоторых одномерных рентгеновских системах для устранения эффектов наложения сигнал перед дискретизацией подвергается низкочастотной фильтрации. В реконструктивной томографии сделать это затруднительно, поскольку приходится иметь дело с

135

двумерными объектами: дискретизация проекционных данных производится по углу (ракурсу) и по проекции при каждом угле.

Эффект наложения по ракурсам данных возникает в том случае, когда угловая дискретизация недостаточна для покрытия угловых гармоник объекта. Единственное видимое решение здесь – использовать большее число ракурсов.

Эффект наложения проекционных данных возникает главным образом из-за использования решеток (блоков) детекторов, состоящих из множества единичных детекторов.

В томографах первого поколения информация о проекции получалась посредством сканирования объекта одним детектором:

P	D
	(x) = ∫ S (l − х) P(l )dl = S (x)* P(x) ,	(2.70)
	0

где P(x) – получаемая при сканировании оценка проекционного

сигнала; S(l) – аппаратная функция детектора; D – диаметр исследуемого объекта; звездочка * означает операцию свертки.

Следовательно, преобразование проекции объекта подвергается низкочастотной фильтрации за счет умножения на фурье-образ ап-

паратной функции детектора S (ν) . Последнее выражение в фурьеобласти можно записать, основываясь на фурье-представлении свертки, как P(ν) = S (ν) P(ν) , где P(ν), S (ν), P (ν) – фурье-

образы соответствующих функций выражения (2.70).

Если для устранения наложения необходима дополнительная низкочастотная фильтрация, то можно использовать дополнитель-

ный низкочастотный фильтр W (ν) .

Поскольку любую требуемую фильтрацию для томографа с одним детектором можно осуществить до дискретизации, тем самым можно устранить ошибки, обусловленные наложением низкочастотных составляющих измеряемого с детектора сигнала.

В томографах третьего поколения, с веерной геометрией пучка излучения, используют решетку единичных детекторов, с помощью которой сразу получают все данные о полной проекции. В этом

случае оценка проекции P(l ) будет равна

136

P(l) = (S (l)* P(l )) Ш	l		,	(2.71)
d

где, как и раньше, S (l ) – аппаратная функция единичного детектора; d – расстояние между детекторами; Ш(x) – гребенка дель- та-функций, определяемая выражением

N

Ш(x) = ∑δ(x − n) ,

n=1

где N – количество единичных детекторов.

Взяв преобразование Фурье от выборочной оценки (2.71), полу-

чим
P(ν) = d (P(ν)S (ν))* Ш(d ν) ,	(2.72)

т. е. решетку частично перекрывающихся спектров.

Здесь в отличие от случая единичного детектора, как это происходит в томографе первого поколения, отсутствует реальная возможность низкочастотной фильтрации, служащей для устранения перекрытия или наложения.

Минимальное перекрытие для максимально широкой аппаратной функции S(l) соответствует детектору шириной d ; эти величины связаны соотношением

l
S (l ) = rect			,	(2.73)
d

где прямоугольная функция rect(.) равна 1 на интервале ±12 и 0 вне его.

Результирующее преобразование S (ν) равно d sin (d ν) и имеет первый нуль на частоте ν = ±1d , а в области высоких частот падает как 1π ν . Поскольку этот спектр повторяется с периодом 1d , результирующие спектры значительно перекрываются. На рис. 2.18 показан спектр оценки P(l ) для объекта исследований в виде широкой «белой» полосы.

137

P (ν)

ν

1/d 2/d

Рис. 2.18. Спектр сигнала на выходе решетки из единичных детекторов

Как уже отмечалось выше, низкочастотный фильтр, который устранил бы имеющееся наложение, в данном случае отсутствует.

Интересный подход к решению этой проблемы основан на применении «смещенных» решеток детекторов, используемых в сканерах с вращающимся веерным пучком [3]. Решетка детекторов сконструирована здесь таким образом, что после поворота на 180° все детекторы сдвигаются относительно объекта на половину ширины детектора. Тем самым, эффективный интервал дискретизации уменьшается вдвое, т. е. становится равным d2 . Составляющие

результирующего спектра повторяются теперь с периодом 2d , а составляющие ±1d устраняются. Таким образом, расстояние ме-

жду центральными частями отдельных спектральных составляющих увеличивается вдвое, что приводит к значительному уменьшению эффектов перекрытия и наложения.

Отсутствие данных о проекциях. В ряде случаев набор дис-

кретизованных данных о проекциях может оказаться неполным. Например, наличие сильно поглощающего объекта может привести к потере части проекций. Точно так же в ряде случаев может оказаться, что объект можно «просветить» с ограниченного малого числа ракурсов. В подобных случаях приходится выполнять реконструкцию по неполным данным. Восстановление изображения по неполным данным требует знания априорной информации о симметрии плотностного и пространственного расположения локаль-

138

ных частей исследуемого объекта для доопределения необходимого набора проекционных данных.

Вопросы малоракурсной томографии и алгоритмов восстановления по неполным проекционным данным более подробно рассматриваются в гл. 6.

2.3.Математические проблемы получения томографического изображения

2.3.1. Постановка задачи

Проблема восстановления структуры многомерного объекта по совокупности его проекций возникла давно [3]. И в настоящее время имеется много методов ее решения. Из всего имеющегося многообразия методов для восстановления изображения можно выделить методы, которые получили наибольшее распространение и сыграли важную роль в развитии реконструктивной компьютерной томографии. Эти методы можно разделить на два основных класса: аналитические и итерационные.

Аналитические методы основаны на точных математических решениях уравнений восстановления изображения. В основе большинства из них используется аппарат преобразования Радона и преобразования Фурье. В связи с этим их можно разделить на две основные группы: обратная проекция с фильтрацией и двумерное восстановление Фурье. При этом методы обратной проекции с фильтрацией имеют три варианта, отличающиеся различными способами фильтрации: радоновская фильтрация, фильтрация Фурье и фильтрация сверткой.

Итерационные методы восстановления изображения используют аппроксимацию восстанавливаемого объекта массивом ячеек равномерной плотности (равномерного μ ), представляющих собой

неизвестные величины линейных алгебраических уравнений, свободными членами которых являются проекции. Решаются системы уравнений итерационными методами, что и дало название данному классу методов восстановления.

139

С математической точки зрения основную задачу реконструктивной томографии, а именно: восстановление изображения по проекциям, можно рассматривать просто, как восстановление некоторой функции в двумерном R2 или трехмерном R3 пространстве по известным значениям интегралов вдоль определенных линий или плоскостей.

В проблеме восстановления функции можно выделить три аспекта. Во-первых, необходимо найти, в какой мере исследуемый объект описывается имеющимися исходными проекционными данными. Во-вторых, необходимо знать, насколько процесс восстановления устойчив относительно погрешностей этих исходных данных. И, в-третьих, необходимо разработать сам алгоритм восстановления изображения.

Первые два аспекта рассматриваются в этой главе, последний – в гл. 3.

2.3.2. Преобразование Радона и формулы обращения

В общем виде преобразование Радона формулируется следующим образом. Пусть в n-мерном пространстве Rn – задана бесконечно дифференцируемая функция f (x) = f (x1,..., xn ) . Так как для практических целей можно рассматривать функцию f (x) , опреде-

ленную в ограниченной области D пространства	Rn и полагать
f (x) = 0 для x D , то в этом случае функция f (x)	автоматически

будет быстро убывать вместе со всеми своими производными. Будем

считать, что пространство Rn ориентировано, и элемент объема в этом пространстве задается дифференциальной формулой

dx = dx1...dxn . Интегрируя функцию f (x) вдоль произвольных гиперплоскостей (ξ, x) = l , где (ξ, x) = ξ1 x1 + ξ2 x2 +... + ξn xn , получим значения интеграла P(ξ,l ) как функцию этих гиперплоскостей. Элемент объема гиперплоскости (ξ, x) = l зададим с помощью дифференциальной формы ω, которая определяется из соотношения

140