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

где Pi(расч) – расчетное значение проекции для li для эффективной энергии Еэф, Pi(расч) = m∫μ(l, Eэф)dl ; т – нормировочный коэффи-

li

циент, i =1, ..., n ; Pi(изм) – измеренное значение проекции при li.

Из системы уравнений (3.112) определяют значение коэффициентов аi, которые используют для определения скорректированных рабочих проекций исследуемого объекта.

Технический аспект. В процессе измерения проекций участвует не только детектор, но и электронная система сбора проекционных данных (ССПД), и это вместе представляет измерительный тракт.

Некоторые погрешности измерительного тракта можно уменьшить, предварительно обработав “сырые” измеренные данные с детекторов по соответствующим алгоритмам. К этим погрешностям относятся:

- систематическая погрешность “темного тока” измерительного тракта, возникающего при отсутствии излучения источника; причинами этой погрешности являются: фоновый сигнал на детекторе, выставляемый уровень “подставки” измерительного тракта для обеспечения его помехоустойчивости, систематическая погрешность ССПД;

-случайные и систематические сбои измерительного тракта, вызванные микропробоями детектора и сбоями электронной части ССПД;

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

Для уменьшения систематической погрешности “темнового тока” проводят его вычитание из рабочих “сырых” данных исследуемого объекта для каждого детектора

ni′ = ni − n0i ,

(3.113)

где ni′ – цифровой код рабочих измерений без “темнового тока”

для i-го детектора; ni – цифровой код рабочих измерений с “темновым током” для i-го детектора; n0i – цифровой код “темнового то-

271

ка” для i-го детектора, полученный при отсутствии излучения источника.

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

а) после получения цифровых кодов рабочих измерений на объекте исследования и измерений на воздушном фантоме;

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

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

торам и ракурсам. В этом случае погрешность σ(Рki ) проекции

Рi = ln n0 ((i)) ,

n1 i

где n0(i) – цифровой код i-го детектора при рабочих измерениях исследуемого объекта, будет определяться выражением [29]

σ(k ) =

∂(Pi ) 2

σ2

+

∂(Pi )

2 σ2

+ 2

∂Pi

∂Pi

r

σ

σ

=

,n

n

n

P

n

∂n1

n

∂n1

n

i

∂n0

0

1

∂n0

0 1

0

1

=

1

2

σn2

+

1

2

σn2

−

2

1

1

rn ,n σn

σn

,

(3.114)

n0

0

n1

1

n0

n1

1

0

1

0

где σ(Pk )

– среднеквадратичное отклонение проекции Pi

для корре-

i

лированных значений кодов n0

и n1;

σn

,

σn

– среднеквадратиче-

0

1

ское отклонение кодов, соответственно n0 и n1;

rn ,n

– коэффици-

ент корреляции кодов n0 и n1/

0

1

Значение σ(РH )

для некоррелированных величин n0 и n1

опре-

i

делится из выражения (3.114) при

rn , n = 0 , как

0

1

σ(Н) =

1

2

σ2

+

1

2

σ2 .

(3.115)

Pi

n0

n0

n1

n1

272

Из выражений (3.114) и (3.115) следует, что σ(Pki ) < σ(PHi ) , а для

частого случая при п0 = п1, σn0 = σn1 , r = 1 значение σ(Pki ) будет равным нулю.

Если рассматривать погрешность на этапе получения цифровых кодов, то суммарная погрешность σΣ от измерений п0 и п1 может

быть	(Н)
			2		2
	σ∑	= σn0		+ σn1 ,							(3.116)
для некоррелированных п0 и п1 и
σ(Н)	= σ2	+ σ2	+ 2r			σ	n0	σ	n1	,	(3.117)
∑	n0	n1		n0 ,n1
для коррелированных п0 и п1.
То есть, значение погрешности σ(Pk )					для проекции при корре-
				i

лированных п0 и п1 может быть значительно меньше погрешности при измерениях кодов п0 и п1 для некоррелированных (3.116) и коррелированных (3.117) кодов п0 и п1. Поэтому коррекцию сбоев в измерительной информации об исследуемом объекте целесообразно осуществлять для коррелированных п0 и п1 в проекциях Pi .

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

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

Линейная фильтрация сбоев основана на методах интерполяции, и самым простым методом этой фильтрации может быть коррекция “по среднему”

P =	Pi−1 + Pi+1	,	(3.118)
i	2

где Pi – корректируемое значение проекции для i-го детектора; Pi–1, Pi+1 – соответственно предыдущее и последующее значения проекции.

273

Медианная фильтрация основана на определенном прогнозировании тенденции возрастания или убывания окружающих фильтруемое значение последующих и предыдущих значений проекции. Для фильтруемого значения Pi определяют из окружающих значений Pi+l, где l = −1, 0, +1 (l может быть равно 3, 5, 7, 9…) медиану по определенной сравнительной логике. На рис. 3.37 отражена эта сравнительная логика для l = 3 .

Значение

P

i −1 i i +1

Pi−1 < Pi < Pi+1

Pi медиана

(Pi = Pi )

i −1 i i +1

Pi−1 > Pi < Pi+1 Pi−1 < Pi+1

Pi−1 медиана

(Pi = Pi−1 )

i	−	1 i i +1		i	−	1 i i +1
Pi−1 > Pi < Pi+1				Pi −1 < Pi > Pi +1
Pi−1 > Pi+1				Pi −1 < Pi +1
Pi+1 медиана				Pi +1 медиана
(Pi = Pi+1 )				(Pi = Pi+1 )

Рис. 3.37. Возможные комбинации значений фильтруемой Pi и окружающих ее значений Pi–1 и Pi+1. Сравнительная логика определения медианы, которая заменяет значение Pi

В общем случае, медианной последовательностью P1, ..., Pl, где l = 2к + 1, к = 1, 2, ..., N, является средний по значению член ряда, получающийся при упорядочении последовательности по возрастанию. Так как для широкого класса исследуемых объектов проекции являются кусочно-монотонными функциями (поскольку соседние лучи пересекают объект приблизительно с одинаковой распределенной плотностью, но при медленно изменяющихся размерах объекта), то для устранения одиночного сбоя (выброса) достаточно взять l = 3. Однако в ряде случаев оказывается полезным увеличить l фильтра, так как это приведет к уменьшению динамического диапазона проекционных данных.

274

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

значения проекции Pi среднего значения из окружающих значений
								l=+		nl
проекцию		Pi:		если		(nl +1)Pi − Pi*		> δ , где Pi* = ∑2			Pi+l , то
								l=−		nl
					2
P(0) =	Pi* − Pi		,	где n	– число значений проекции, окружающих
i	nl	−1		l

корректируемое значение проекции Pi; Рi(0) – присваемое значение

проекции вместо Pi при коррекции; δ – величина, определяемая предел сглаживания Pi (выбирается экспериментально).

3.5. Методы восстановления изображения по проекционным данным

3.5.1. Классификация методов восстановления изображения

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

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

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

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

275

рование, и алгоритмы, использующие двумерное преобразование Фурье.

Тот и другой метод предполагают, что известно точное значение проекций P(l,θ) (для параллельной схемы) или P(γ,β) (для веерной схемы сканирования) для всех l (или γ) и θ (или β) и что требуемые интегральные преобразования можно выполнить точно. Однако реально ни то, ни другое условие выполнить невозможно. Измерение реальных проекционных данных проводится с определенной погрешностью; набор, проекций и значений в этих проекциях не для всех l (или γ) и θ (или β), а для конечного числа их дискретных значений. Расчет интегральных преобразований в реконструкторе томографа проводится в дискретной форме по дискретным операторам преобразования, которые действуют на требуемые функции, имеющие конечное число компонент. Невыполнение условий приводит к неустойчивости решения задачи реконструкции изображения, что выражается в виде значительных артефактов на томограмме.

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

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

Метод реконструкции, основанный на разложении искомой функции (в нашем случае функции μ(х, у) ) в ряд, принципиально

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

Алгоритмы с использованием разложения в ряд можно разделить на две группы: итерационные и неитерационные.

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

276

известные величины линейных алгебраических уравнений, свободными числами которых являются проекции. Решаются уравнения итерационными методами, что и дало название данному классу алгоритмов реконструкции. В настоящее время известно несколько итерационных методов реконструкции. Отличаются они, в основном, последовательностью внесения поправок во время итерации. Среди них наиболее известны следующие методы: алгебраический метод реконструкции томограммы (АRТ), метод одновременного итерационного восстановления (SIRT), итерационный метод наименьших квадратов (ILST) и мультипликативный алгоритм алгебраической реконструкции (MART).

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

Выше рассмотренную классификацию методов реконструкции изображений по проекциям можно представить в виде схемы рис. 3.38.

		Алгоритмы реконструкции изображений
	Алгоритмы с использованием			Алгоритмы с использованием
интегральных преобразований					разложения в ряд
Алгоритмы		Алгоритмы
		двумерного		Итерационные		Неитерационные
обратного
		преобразования		алгоритмы		алгоритмы
проецирования
		Фурье
Алгоритмы с фильтрацией сверткой	Алгоритмы с фильтрацией Фурье	Алгоритмы радоновской фильтрацией	Алгебраический метод реконструкций (ART)	Метод одномерного итерационного восстановления (SIRT)	Итерационный метод наименьших квадратов (ILST)	Мультипликативный алгоритмалгебраических реконструкций (MART)
		с

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

277

3.5.2. Особенности аналитического метода восстановления изображения с использованием обратного проецирования с фильтрацией сверткой

Алгоритмы, основанные на методе интегральных преобразований (аналитический метод) состоят из ряда этапов:

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

2)нахождение формул обращения и определение по ней неизвестной функции;

3)адаптация формулы обращения к дискретизированным зашумленным данным.

Рассмотрим первый этап.

Обозначим двумерное распределение физической величины

функцией μ(х, у), вид которой априори не известен. Однако из-

вестно, что в большинстве приложений она ограничена в пространстве, т. е. равна нулю вне некоторой конечной области плоскости, обозначаемой далее через Ω. Можно полагать, что функция μ определена областью Ω, представляющей собой круг радиуса Т с центром в начале координат.

Иногда удобнее записывать функцию μ не в прямоугольных

(х,у), а в полярных координатах (r, φ). В этом случае
μ(х, у) = μ(r cosφ, r sin φ) .	(3.119)

Прямая на плоскости может быть задана двумя параметрами: расстоянием l (со знаком) от начала координат и углом θ относительно оси у (рис. 3.39).

Положение точки P на плоскости определяется ее координатами (х,у) или (r cosφ, r sin φ) , а положение рентгеновского луча – его расстоянием l от начала координат и углом θ.

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

278

T
P(l,θ) = ∫ μ(l cosθ−t sin θ, l sin θ+ t cosθ)dt ,	(3.120)

−T

где пределы интегрирования в общем случае зависят как от параметров l и θ, так и от области Ω.

y

t	T
			T (l)
		θ	P
				l
Ω		r
		l '
Реконструируемый
объект		l	θ	x
		0

− T (l)

Радиус T

Луч

Детектор

− T

Рис. 3.39. Параллельная геометрия сканирования

Функция P(l,θ) называется проекцией, а правая часть выражения

(3.120) – лучевой суммой. Физическая интерпритация выражения (3.120) для рентгеновского излучения показана в гл. 2. Аргументы функции μ(х,у) в выражении (3.120) взяты в координатах (l, t). Связь между неподвижной системой координат (х,у) и подвижной (l, t) задается формулами Эйлера

		х = l cosθ−t sin θ,	(3.121)
		у = l sin θ+t cosθ,

и

279

P(l,θ) = P(−l,θ+ π) .

l = xcosθ+ y sin θ,

(3.122)

t = −xsin θ+ y cosθ.

Когда реконструируемый объект представлен в виде круга радиусом T, будем иметь

T (l ) = (T 2 −l2 )1 2 ,						l		≤T ,	(3.123)
P(l,θ) = 0,		l		>T .					(3.124)
Необходимо отметить, что при		любых l							и θ пары (l,θ) и

(−l, θ + π) задают одну и ту же прямую, поэтому

(3.125)

Следует также отметить, что аргументы функции P(l,θ) отлича-

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

θ1 и θ2 . В общем случае интегралы P(0,θ1 ) и P(0,θ2 ) от функции μ из выражения (3.120) по этим прямым будут различны, но этой

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

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

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

всех l, проходящих через искомую точку.

280