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

Рис. 3.42. График трех функций пространственной частоты ν, произведение которых дает фурье-преобразование свернутых и интерполированных проекций

291

Кривая 2 соответствует функции окна Шеппа–Логана [48], для которой

1

sin (πνΔl )

,

ν

≤1 2

l ;

gˆ (ν) =

πΔl

(3.151)

0,

ν

>1 2

l.

Кривая 3 соответствует выбору

функции окна

для функции

∑ (ν,θ), если в ней доминируют эффекты наложения и шумы при

14 l < ν < 34 l . На рис. 3.42 также показан график и для двух функций J (l ) , соответствующих идеальному фильтру ниж-

них частот (обозначены сплошной линией) и линейной интерполяции (пунктирной линией), при которой

	ˆ	sin(πνΔl)		2
	Ja (ν) =			.	(3.152)
			πνΔl
	Следует отметить, что идеальный фильтр нижних частот при-
	водит к резкому усечению функции Pu			в частотной области, что

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

функции	Pu , что аналогично эффекту Гиббса при усечении ряда
Фурье.
Влияние вида интерполяции		показано на графике функции
и (ν,θ)	для двух вариантов J (l )	и для варианта gˆΣ(ν,θ) , обозна-

ченного цифрой 1. По сравнению с (ν,θ) (выражение (3.148)) график функции и (ν,θ) искажен ошибкой наложения, которая обусловлена усилением исходного наложения в Σ(ν,θ) . На графике, соответствующем случаю линейной интерполяции (пунктирная линия), видно, что и (ν,θ) также содержит ложные высокочастотные компоненты при ν >12 l , обусловленные кусочной линейностью функции pи (l ', θm ) (3.144).

Заметим, что как ошибки наложения, так и ложные высокочастотные компоненты существенно уменьшаются при использовании

292

функции gˆΣ(ν) график которой обозначен цифрой 3, однако при

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

Рассмотрим влияние оценки интеграла от обратной проекции (3.140), вычисленного по формуле трапеций вида (3.142).

Для того чтобы отделить влияние дискретизации по параметру l и выделить дискретизацию по углу θ, будем считать, что функция

p(l ',θ) в точности известна при произвольных значениях l ' . Ис-

пользуя это упрощение, найдем функцию передачи рассеяния точки для алгоритма реконструкции, т. е. рассмотрим реконструкцию “им-

пульсного” объекта (тонкий стержень) с координатами (ξ,η). Объект исследования 0(x, y) и его преобразование Фурье

ˆ	( X ,Y ) описываются функциями
0
						0(х, у) = δ(х−ξ, у −η) ,		(3.153)
					ˆ	( X ,Y ) = exp(−i 2π( X ξ +Y η)).		(3.154)
					0
	Пусть 0(x, у)					– импульсный сигнал с ограниченной шириной
спектра, соответствующей указанному импульсному объекту
0(х, у) =		π νс			(νcosθ,νsin θ)W (ν)exp i 2πν		(xcosθ + уsin θ)
		∫	∫	ˆ					ν		dνdθ,
				0

0 −νс

(3.155)

где νс – частота среза функции окна.

Функция 0(x, у) представляет собой выходной сигнал, полученный за счет процедуры реконструкции по (3.136). Подставив

0(x, у) из (3.154) в (3.155), получим

	π νс
	∫ ∫
0(х, у) =		exp i 2πρνcos(θ− γ) W (ν)		ν	dν dθ, (3.156)

0 −νс

где tg γ = ух−−ηξ , ρ = (х−ξ)2 + ( у −ξ)2 12 .

293

Поменяв порядок интегрирования, увидим, что интеграл по θ не что иное, как интегральное представление функции Бесселя J0 (2πρν) [69]. Таким образом,

νс
0(х, у) = π ∫ J0	(2πρν) W (ν) ν dν.	(3.157)

−νс

€ ( )

Нас интересует функция 0g , , которая получается в случае,

x

у

когда интеграл по θ в (3.155) вычисляется по формуле трапеции, а не аналитически. По определению будем иметь:

	g		M	νс			m
			∑	∫
0		(х, у) = Δθ			exp i 2	π ρνcos(θ		− γ) W (ν)	ν	dν. (3.158)

m =1−νс

Используя свойства производящей функции для функций Бесселя [69], а именно:

∞
еi z cosθ = ∑ im Jm (Z ) ei mθ ,	(3.159)

m =−∞

получим из (3.158)

νc

∫ −∞

0g (x, у) = π ∑

K =−∞

−νc

(−1)M K e−i 2M K γ J2NK (2πρν) W (ν) ν dν.

(3.160)

Из сравнения результата (3.160) с формулой (3.157) следует, что член суммы с K = 0 равен J0 (2πρν) , т. е. 0(х, у) = 0g (х, у) при K = 0 . Таким образом, сумма (3.160) по ненулевым K представляет собой ошибку аппроксимации интеграла по θ для формулы трапеции.

Как следует из асимптотического представления функции Бесселя [69], она имеет первый максимум при Z > m , причем при Z ≤ m − 2 ее значения малы по сравнению с этим максимумом. Необходимо также отметить, что член суммы (3.160) при каждом нену-

294

левом K периодически изменяется по γ с периодом π/MK. В соответствии с этим из (3.160) получим, что равенство 0g (х, у) 0(х, у) будет выполняться для ρ ≤ ρ' , где

	ρ' =		M	.		(3.161)
			πνc
Собственно условие	(3.161)	формально			вытекает из	условия
Z = m при K =1.
Учитывая исходное	предположение об				импульсном	объекте

(сигнале) с ограниченной шириной спектра, приходим к выводу о возможности реконструкции изображения, которое является копией с ограниченной шириной спектра заданного импульсного объекта, для точек, отстоящих от места локализации импульса не далее, чем на расстоянии ρ' . Если рассмотреть поле реконструкции в ви-

де круга диаметром DΩ , как суперпозицию импульсов, то во избежания ошибок численной аппроксимации интеграла от обратной проекции (3.158), необходимо, чтобы ρ' > DΩ . Учитывая, что DΩ = N l , νc =12 l , получим из (3.161) неравенство, связы-

вающее число ракурсов M и число лучей (единичных детекторов и отсчетов) N :

М > πN 2 .

(3.162)

Это неравенство определяет выбор количества ракурсов и единичных детекторов при практическом проектировании томографа.

3.5.3.Отличительные особенности итерационных методов восстановления изображения

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

Восстанавливаемый объект аппроксимируется массивом с равномерным коэффициентом линейного поглощения в ячейке, как показано на рис. 3.43. Коэффициент линейного поглощения i-ой ячейки обозначается μi . Для того чтобы исключить обработку бес-

295

конечного количества данных, обычно ограничивают пространст-
венное разрешение, выбирая конечный размер ячейки d. Тогда об-
щее количество	ячеек в области с	радиусом	R будет равно
K = πN 2 4 , где N = 2R d .
Аналогичным образом можно ограничить количество данных в
проекции. С этой целью каждая проекция разбивается на полосы
(см. рис. 3.43) конечной ширины (в отличие от идеализированных
лучей). Поскольку разрешающая способность проекционных дан-
ных должна быть сравнима с разрешающей способностью изобра-
жения, каждая проекция делится на N лучей с шириной d .
	y
t	2	3	l
R		d
1			x
			d
			N
Рис. 3.43. Аппроксимация исследуемого объекта и проекции в алгебраических
	методах реконструкции:
1 – i-я ячейка размером d ×d ; 2 – заштрихованная область – j-й луч конечной
ширины; 3 – площадь пересечения i-й ячейки с j-м лучом

Предположим, что имеется множество K базисных функций ωi (x, y) линейная комбинация которых дает адекватную аппрок-

296

симацию некоторой функции μ(x, y), т. е. существуют действительные числа μ1,...,μi ,...,μK , такие, что

K
μ(x, y) ≈ μ(x, y) = ∑μi ωi (x, y) .	(3.163)
i =1

Задача состоит в том, чтобы по имеющимся проекциям функции μ(x, y) восстановить саму функцию. Пусть θ – оператор

отображения функции μ(x, y) в проекцию под углом θ. Применим

этот оператор к уравнению (3.163), получим для некоторого значения l проекции

р(l,θ) =	K μ ω	(l ) .	(3.164)
	∑ i θ i

i =1

Выражение (3.164) представляет собой систему линейных алгебраических уравнений, в которой неизвестными являются числа

μi , проекции же базисных функций ωi (l ) для некоторых углов θ и величин l могут быть вычислены.

Проекции искомой функции μ(x, y) обычно определяются из

экспериментальных данных, в которой ширина луча является конечной. В этом случае, используя уравнение (3.164), получим

l2		K		l2	θ	i
∫	p(l,θ) dl =	∑	μ	i ∫				(3.165)
					ω		(l ) dl .
l1		i=1		l1

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

Независимо от того, какое из уравнений (3.164) или (3.165) будет использоваться, мы получим систему линейных алгебраических уравнений, которую нужно разрешить относительно μi . Однако не

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

быть так велико, что прямые методы оказываются непригодными.

297

Во-вторых, число уравнений, которые можно получить из проекционных данных, может быть меньшим, чем число неизвестных μi , и в этом случае необходимо применение оптимизационных ме-

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

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

Существует много способов выбора базисных функций. Рассмотрим базисные функции, которые в основном используются в итерационных алгоритмах восстановления. Пусть восстанавливае-

мая функция μ(x, y) аппроксимируется системой из M квадратных ячеек (рис. 3.43) со средним значением коэффициента линейного поглощения в каждой ячейке μi (i =1,2,..., M ) . Базисные функции

определим следующим образом:
ω	(x, y) =		1 внутри i−й ячейки	(3.166)
i			0 везде, кроме i−й ячейки .

Тогда лучевые суммы (или проекции) для лучей конечной ширины будут представлять собой сумму значений ячеек, пересекаемых лучом, т. е. для j-го луча

р( j,θ) = μ1 ωj,1 + μ2 ωj,2 +... +μK ωj, K , 1 ≤ j ≤ L , (3.167)

где ωj,i – площадь пересечения j-м лучом i-й ячейки, L – полное число лучей во всех проекциях, L должно быть не менее K, т. е. N ×M = L должно быть не мене K .

При определении ωj,i в алгоритмах используются два варианта: весовой, когда интеграл в правой части уравнения (3.165) подсчитывается точно, и безвесовой, в котором интеграл приравнивается единице, если центр i-й ячейки лежит внутри луча, и приравнивается к нулю, если центр i-й ячейки не лежит внутри луча. Второй вариант имеет преимущество в скорости вычисления, однако при этом может вводиться дополнительная аппроксимация.

Следует отметить, что большинство коэффициентов в уравнении (3.167) равны нулю, так как один луч пересекает небольшое число

298

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

Выбор критерия для решения уравнений (3.167) зависит от количества и качества исходных данных в проекциях. Так, в случае недостаточного количества проекций (т.е. количество проекций меньше определенного числа неизвестных) решение μi обычно выбираются с помощью критерия, минимизирующего общее количество отклонений от среднего значения μi . В случае сильно зашум-

ленных проекционных данных искомые величины μi определяются с помощью критерия, сводящего к минимуму среднеквадратическую

L	K	2
погрешность проекции ∑ Pi − ∑ωj,i μi .
	i =1
j =1

При выборе способа численного решения уравнений (3.167) чаще всего используют итерационные алгоритмы. Решение системы уравнений итерационным путем сводится к уточнению значений μi

таким образом, чтобы расчетные проекции соответствовали измеренным проекциям. Для этого выбирается исходное значение μi

ячеек. Затем по исходным μi с помощью уравнений (3.167) вычис-

ляются проекции и сравниваются с измеренными. Предположим, что рассчитанная лучевая проекция меньше измеренной величины. Тогда каждая ячейка, принимающая участие в вычислении значения этого луча проекции, подвергается изменению на величину Δμi ,

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

В общем виде процесс коррекции во время K-й итерации можно описать уравнением

	L
	μiK = μiK −1 + ∑ΔμiK, j ,	(3.168)
	j =1
где μiK −1	– коэффициент линейного поглощения i–й ячейки до K-й
итерации,	μiK – после итерации, ΔμiK, j – поправка, вносимая в i-ю
	299

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

Вопрос о том, когда следует прекратить итерационный процесс, является остается открытым. Обычно число итераций находится в пределах 4–10. В качестве критерия конца итерационного процесса могут быть заданы, например, допустимые пределы погрешности изображения или скорость уменьшения погрешности.

3.5.4. Модифицированный метод восстановления изображения для случая веерных рентгеновских лучей

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

На рис. 3.44 показаны геометрические соотношения для веерного луча. Луч задан двумя параметрами: γ и β. Пусть 0 – начало

координат, а U – положение источника излучения на окружности радиусом Д с центром 0. Тогда β – угол наклона прямой 0U . Пусть

R – точка у основания перпендикуляра, опущенного из начала координат к лучу. Тогда θ – угол между перпендикуляром 0R и осью х, а l – длина перпендикуляра 0R .

Каждый луч, задававшийся в параллельной геометрии парой чисел (l,θ) , теперь рассматривается как элемент множества

расходящихся лучей (γ,β) , где угол β задает положение источника излучения, а γ – угол между рассматриваемым лучом и прямой, соединяющей источник и центр окружности (см. рис. 3.39). Обозначим линейный интеграл P(l,θ), принятый для параллельной

300