Симонов Томографические измерителные информационные системы 2011

геометрии, через q(γ,β) – для веерной геометрии, причем l < Д – радиус окружности, по которой движется источник излучения.

Рис. 3.44. Веерная геометрия сканирования
Из рис. 3.44 следует
l = Дsin γ , θ = β+ γ , P(l,θ) = q(γ,β) ,	(3.169)
r cos(θ−φ) −l = S sin (γ'− γ),	(3.170)
где S – расстояние между точками U и P ; γ' – угол между отрез-

ками 0U и UP .

Соотношение (3.170) вытекает из подобия треугольников 0RL

и LPM , где 0L = l cos(θ− φ) , а РМ = PM = S sin (γ'− γ).
Из	геометрических соотношений следует	равенство
S 2 = (UK )2 + (PK )2 . Учитывая, что (UK )2 = Д+ r sin		(β −φ) 2	, а
(PK )2 = r cos(β −φ) 2 , то
	301

S 2 = r cos	(β −φ) 2		+ Д+ r sin		(β −φ) 2	=
							(3.171)
r2 + Д2 + 2Дr sin (β −φ),
		r cos(β −φ)
γ' = arctg				.			(3.172)
		Д+ r sin (β − φ)

Из формул (3.130), (3.135) и (3.133) для параллельной геометрии луча следует, что

						∞ T 2π
μ(r cosφ,r sin φ) =1 2 ∫							∫ ∫ p(l,θ)×
	(		(		)	−∞−T 0
							)
×exp		i 2πν r cos		θ−φ		−l		W (ν)	ν	dθdl dν.
Перейдя от переменных (l,θ)						к переменным (γ,β) ,

во внимание соотношения (3.169) и (3.170), получим

μ(r cosφ,r sin φ) =

(3.173)

принимая

∞ γm	2π−γ
= (Д 2) ∫ ∫	∫ q(γ,β)cos γexp(i 2πν)S sin (γ'− γ)W (ν)		ν		dβdγdν.
−∞−γm	−γ
	(3.174)

Учитывая, что β изменяется с периодом 2π , пределы интегрирования по β можно изменить от 0 до 2π .

Как и для параллельной геометрии лучей (3.139) и (3.140), операцию реконструкции по формуле (3.174) можно представить в два этапа:

γm
q(γ',β) = ∫ q(γ,β)cos γ (S sin (γ'− γ))dγ ,	(3.175)
−γm
2π
μ(x, y) = (Д 2) ∫ q(γ',β)dβ,	(3.176)
0

где

∞

g (S sin (γ'− γ)) = ∫ ν W (ν) exp(i 2πνS sin (γ'− γ))dν. (3.177)

−∞

302

Рассмотрим выражение (3.177). Ядро g (l ) в соответствии с (3.134) определяется как обратное преобразование Фурье от модуля пространственной частоты ν

∞
g (l ) = ∫ ν exp(i 2πνl ) dν .	(3.178)

−∞

Здесь следует отметить, что выражение (3.178) представляет неотфильтрованное ядро, поэтому выражение (3.178) не совпадает с выражением (3.138).

Для значения аргумента по выражению (3.177) уравнение (3.178) примет вид

∞

g (S sin γ) = ∫

ν

exp(i 2πνS sin γ) dν ,

(3.179)

−∞

где γ = γ'− γ .

Введя новую переменную

S sin γ

ν'

= ν

,

(3.180)

γ

получим

∞

γ

γ

g (S sin γ) = ∫

ν'

exp(i 2πν'γ)

dν' ,

(3.181)

S sin γ

S sin γ

−∞

откуда

g (S sin γ) =

γ2

∞

exp(i 2πν'γ) dν' .

∫

ν'

(3.182)

S

2

sin

2

γ

−∞

Интеграл в выражении (3.182) представляет ядро с аргументом

γ , полученное для параллельных

лучей. Следовательно, ядро для

веерных лучей может быть определено в виде

γ2

g (γ) =

g (γ) ,

(3.183)

S 2 sin2 γ

где g (γ) – ядро, определяемое выражением (3.178). В результате (3.176) для случая веерной геометрии сканирования примет вид

303

μ(x, y) =		1	2π			1	γm	q (γ,β) g (γ'− γ) dγ ,
			∫				dβ ∫		(3.184)
						2
	2		0			S	−γm
где q (γ,β) = q(γ,β)Дcos γ							–	модифицированная	проекция;
g (γ'− γ) = g (γ'− γ)	γ2				– модифицированное ядро.
	sin2 γ

Следовательно, выражение (3.177) представляет собой модифицированное ядро g (γ), которое определяется как для случая па-

раллельных лучей. Тогда для оценки W(ν), определяемого выражениями (3.146) и дискретного представления проекции (3.147), модифицированное ядро для веерной геометрии сканирования определиться, как

g (i,Δγ)			=	3 − 2α				, i = 0 ,
				α2(Δγ)2
g (i,Δγ) = −				α			, i		– четное,	(3.185)
		π2 sin2 (i Δγ)
g (i,Δγ) = −			1−α		, i				– нечетное,
	π2 sin2 (i Δγ)

где Δγ – угол между соседними детекторами, i = –N/2, ..., 0, ..., N/2–1. На практике, чтобы можно было использовать ЭВМ для решения задачи реконструкции, уравнение (3.184) предварительно аппроксимируют. Используя простейший вид аппроксимации, т.е. заменяя интегрирование суммированием, можно представить алго-

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

1. Модификация полученных экспериментальных проекций

q (γi ,βγ ) = q(γi ,βj )Дcos(γ j ),	(3.186)
где i = –N/2, ..., 0, ..., N/2–1 (для N – четно), j =1,..., M ;	q(γi ,βj ) –

исходные экспериментальные проекции.

2. Вычисление ядра свертки по формуле (3.185). Для коэффициента регуляризации α = 0 будем иметь, например, ядро Лакшминараянана−Рамачандрана

304

		1			,
	4(Δγ)2
g (i Δγ) = 0,
					1
	−						,
			2		2	(i Δγ)
		π		sin

i = 0;
i четное;	(3.187)
i нечетное.

3. Свертка проекции с ядром, т. е. вычисление внутреннего интеграла в (3.184)

N 2−1
g (γi ,βj ) = Δγ ∑ q (γi ,βj ) g (γl − γi ) ,	(3.188)

i =−N2

где i = –N/2, ..., 0, ..., N/2–1.

4. Обратное проецирование, т.е. вычисление внешнего инте-

грала (3.184)

	M		g (γ'p ,βj ),
μ(x, y) =	Δβ ∑	1		(3.189)
	2 j =1 S 2
поэтапно:
а) определение γ'p по формуле (3.172); γ'p				– угол для каждой

точки Р восстанавливаемой матрицы изображения Mx ×M y , Δβ – угловой ракурс, Δβ = 2πМ ;

б)	интерполирование свертки с целью определения	g(γ'p , βj )
по дискретным экспериментальным значениям g(γl , βj ) ;
в)	вычисление множителя 1/S2 по формуле (3.171);	g(γ'p , βj ) ,
г)	обратное проецирование свернутой проекции

предварительно умноженной на множитель 1/S2.

3.5.5. Сравнение методов восстановления изображения

Исторически сложилось так, что алгебраические методы реконструкции изображения были применены на первой установке рентгеновского томографа для головы человека, разработанной Г. Хаунсфилдом – будущим луареатом Нобелевской премии. Однако в

305

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

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

Реконструкция изображений по неполному набору проекций.

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

новить μ(х, у) вне этого включения в отсутствии “лучевых сумм”

по пересекающим его прямым.

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

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

306

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

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

3.6. Требования к вычислительным процессам восстановления изображения

Оценим объем вычислений для алгоритма реконструкции методом обратного проецирования с фильтрацией сверткой для веерной схемы сканирования ((3.186)–(3.189)), как одного из применяемых в практике.

Первый этап алгоритма (3.186) требует N ×M количества умножений на заданную величину Дcos γi , которая может быть вы-

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

Второй этап (3.187) вычисления определяет вектор – строку значений ядра свертки размерностью 1× N и может быть определен заранее.

Третий этап – определение свертки модифицированных проекций с ядром. Количество умножений вычисленной на первом этапе модифицированной проекции и вычисленного ядра будет N , т. е.

после третьего этапа требуется провести		(N × M )× N операций
умножения.	γ′		и коэффициента 1/S2
Четвертый этап – определение углов
	p

для матрицы изображения M x ×M y элементов (пикселов) требуют

их вычисление для каждого элемента изображения, т. е. их вычисление по M x ×M y раз. Углы γ′p и коэффициент 1/S2 можно вы-

числить заранее.

Интерполирование свертки по ее дискретным значениям требуется выполнить M x ×M y раз. Обратное проецирование интерполи-

рованной свернутой проекции (3.189) требует выполнения опера-

307

ций умножения М раз, т.е. этот этап требует выполнения M ×(M x ×M y ) операций умножения.

Таким образом, алгоритм реконструкции обратного проецирования с фильтрацией сверткой для веерной геометрии луча требует (N × M ) × N + M ×(M x × M y ) операцийумножения.

Если условно взять следующее соотношение между количеством детекторов, ракурсов и элементов матрицы изображения, а именно: N = M = M x = M y = K , то минимальное количество параметров,

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

Таблица 3.2

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

Параметры	Соотно-		Число параметров
	шения
Независимые элементы
изображения по диа-	K		512		1024	2048
метру объекта исследо-
вания
Элементы (пикселов)	πK 2	4	2,0 105		8,2105	3,3106
томограммы
Отсчеты в проекции			512		1024	2048
(единичных детекто-	K
ров)
Проекции (ракурсы)	K		512		1024	2048
Проекционные (изме-	K 2		2,6 105		1,0 106	4,2 106
рительные) данные
Операции умножения	2K 3		2,6 108		2,2 109	1,7 1010
Количество операций		умножения		требуется		провести

(0,02−2) 1010 . Если томограмму получать в реальном масштабе

времени, т. е. после сканирования объекта, которое длится обычно 3–6 с, должно быть получено реконструируемое изображение, то реконструктор должен обладать скоростью вычислений

308

0,4 (109 −1010) операций в секунду. В настоящее время это воз-

можно достигнуть с помощью специального процессора реконструкции или быстродействующей универсальной ЭВМ.

3.7.Вычислительные средства и программное обеспечение

3.7.1.Характеристики ЭВМ для восстановления

иобработки изображения

Для восстановления изображения в реальном масштабе времени каждая операция алгоритма должна занимать не более t/K2 с, где t – время получения одной проекции. Обычно для реальных томо-

графов третьего поколения t находится в пределах (3 −9) 10−3 с. Это значит, что каждая операция (например, для K = 512) не должна быть более tK 2 < (11,5 −34,6) нс. Такую скорость возможно по-

лучить, если операции вычисления осуществлять через оперативное ЗУ.

Объем оперативной памяти, требуемый для запоминания элементов изображения, определяется этапами свертки и обратного

проецирования и пропорционален 2K3 . Если вычисления вести с высокой точностью, определяемой разрядностью пр, то объем ОЗУ должен быть

V = 2np K3 бит,

(3.190)

и для пр = 32 и K = 512 V = 2,1 Гбайт.

В настоящее время такие характеристики достижимы как при разработке специального процессора, так и для последних версий ЭВМ типа Pentium–IV.

3.7.2. Задачи обработки изображений

Задачи обработки изображений в реконструктивной томографии рассматриваются, как операции над некоторым множеством

квадратных матриц M (Z ) с элементами целых чисел Z , где эле-

309

менты матрицы – соответствующие значения коэффициента линейного поглощения μ(х, у) элементарных объемов сечения ис-

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

В общем виде алгоритмы обработки изображения обеспечивают отображение матрицы M (Z ) в себя либо отображение некото-

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

На видеотерминале элементу матрицы значений μ(х, у) соответствует элемент площади (пиксел) со свечением, пропорциональным значению μ(х, у).

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

матрице μ(х, у), а также линии, контуры, точки и соответствующая

цифробуквенная информация о томограмме.

Рассмотрим некоторые принципиальные задачи обработки изображения.

Выбор диапазона плотностей для визуализации изображения.

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

воспроизвести сразу на видеотерминале невозможно, поскольку градационная способность современных терминалов составляет не более 256 полутоновых ступеней. Для реализации разрешающей способности томографа по плотности при использовании всех полутоновых ступеней градации “серо-белого” предусматривается понятие

“окна” изображения W (Hv) , ширина и положение которого в диа-

пазоне –1024 – +3071 Hv регулируется оператором.

На рис. 3.45 показано прямолинейное “окно” изображения, которе выделяет линейно μ(х, у) в диапазоне W (Hv) в интервале ярко-

сти 0–256, остальные значения μ(х, у)отображаются “нулевой” или “максимальной” яркостью. Значения μ(х, у) в “окнах” могут изменяться по более сложным законам (нелинейно, в виде импульсной

310