Беляев Физика ядерной медицины Ч.2 Учебное пособие 2012

2.5. Регуляризация

Инверсия уравнения (2.1) относится к плохо обусловленным задачам и его решение (2.3) является нестабильным по отношению к небольшим возмущениям данных p(s, ). Такая ситуация может

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

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

f (x, y) f (x, y) F{W (vs ) | vs | F1{ p(s, )}}d	(2.5)
0
где f (x, y) – реконструированная оценка f(x,y); W(vs) –	сглажи-

вающая (аподизационная) функция.

В литературе по реконструкции томографических изображений термин "фильтр" закреплен за единственным пилообразным фильром |vs|. В то время как функция сглаживания может иметь произвольный вид, который представляется предпочтительным с точки зрения отношения сигнал/шум или по другим соображениям. Отметим, что если W(vs) ≠1, то уравнение (2.5) является смещенной оценкой f(x,y).

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

данных и некоррелированного случайного шума, pM (s, )

p(s, ) n(s, ) , то соответствующий измеренный спектр

мощности равен PPM (vs , ) Pp (vs , ) P n (vs , ). По причине конечной частотной чувствительности системы детектирования мощность сигнала (PP) будет постепенно затухать, в то время как мощность шума (Pn) будет оставаться практически постоянной, что демонстрируется на рис. 2.6. Эффект пилообразного фильтра, также показанный, будет усиливать высокочастотные компоненты спек-

61

тра мощности, которые станут доминирующими в шуме при высоких частотах. Общепринятая аподизационная функция представляет косинусоидальное аподизационное окно, называемое также окном Ханнинга. Специфическая форма W(vs) будет определять компромисс шум/разрешение в реконструированном изображении и также влиять на структуру шума в изображении. На практике затруднительно строго обосновать выбор функции W(vs), в принципе, ее выбор должен оптимизировать конкретную задачу, например, детектирование болезни.

Рис. 2.6. Иллюстрация использования W(vs,)|vs|, аподизирующей пилообразный фильтр для подавления усиления мощности шума высокой частоты выше граничной частоты vc [3]

3. Аналитическая 3-М реконструкция изображений

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

62

3.1. Алгоритм 3-М обратного проецирования

Пространственное изменение чувствительности сканера возникает как следствие конечной аксиальной протяженности сканера, приводящее к усечению проекций. Для 2-М кольцевого сканера регистрируемая интенсивность от точечного источника остается почти постоянной независимо от положения источника в пределах FOV, что иллюстрируется на рис. 2.7. В то время как для 3-М цилиндрического сканера детектируемая интенсивность точечного источника изменяется внутри FOV, особенно при аксиальном перемещении источника.

Рис. 2.7. Иллюстрация пространственной инвариантности 2-М визуализации и пространственного изменения 3-М визуализации при рентгеновском преобразовании (адаптировано из [3])

Пространственная изменчивость чувствительности усложняет процесс реконструкции. Например, если область определения не измеряется из-за усечения, то точный расчет 2-М фурьепреобразования проекций для алгоритма реконструкции с помощью метода фильтрованного обратного проецирования становится невозможным при простом использовании быстрого преобразования Фурье.

Общепринятым методом восстановления пространственной инвариантности измеряемых 3-М рентгеновских проекций является алгоритм 3-М обратного проецирования [7]. В этом методе неизме-

63

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

сти [8].

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

3.2. Методы перегруппировки

При 3-М рентгеновской трансформации данных, применяя определенную форму усреднения сигналов, можно использовать измеренные данные для оценки массива 2-М поперечных синограмм. Такая процедура называется алгоритмом перегруппировки (или пересортировки). Преимущество алгоритмов перегруппировки заключается в возможности эффективной реконструкции каждой перегруппированной синограммы либо аналитическими, либо итерационными методами. Дополнительно, перегруппировка может существенно уменьшить размер данных. Недостаток методов перегруппировки состоит в появлении пространственно изменяющегося возмущения и/или усиления статистического шума. Хотя методы перегруппировки не относятся, по существу, к процедурам реконструкции, они служат важным дополнением к семейству методов, на которые опирается решение проблемы 3-М реконструкции изображений.

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

64

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

4. Итеративная реконструция изображений

4.1. Основные элементы

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

Все итеративные методы содержат пять основных компонентов. Первый и часто упускаемый из виду компонент заключается в модели изображения. Обычно это дискретизация области изображения на N отдельных пикселей (2-М элементы изображения) или вокселей (3-М элементы изображения). В некоторых моделях предлагаются сферические элементы [10,11] с перекрывающимися границами. Но на практике из-за сложности расчетов они не используются. Как правило, элементы изображения считаются детерминистскими переменными, но в некоторых случаях эти же элементы могут рассматриваться как случайные переменные и при реконструкции изображений может применяться байесовский подход.

Второй основной компонент представляет модель системы, связывающую изображение и данные. Элемент Hi,j модели системы H характеризует систему визуализации и представляет вероятность того, что эмиссия из вокселя j детектируется в проекции i. Поэтому

65

N
pi Hi, j f j ,	(2.6)
j 1
где pi – среднее значение i-проекции; fi	– активность в вокселе j
(рис. 2.8).

Рис. 2.8. Иллюстрация одиночного элемента модели системы Hi,j (адаптировано из

[3])

Большинство клинических методов использует пространственно инвариантные модели системы с функциями отклика, упрощенными для повышения эффективности расчета. В недавних исследованиях показано, что применение более точных, пространственно зависимых моделей системы улучшает разрешение и уменьшает погрешности позиционирования в системах высокого разрешения для небольших животных [12] и при ПЭТ-визуализации всего тела

[13].

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

66

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

M	p pi exp( p )
P(P p | f )	i	i	(2.7)
		pi !
i 1

Хотя модель Пуассона вполне подходит для ПЭТ визуализации с концептуальной точки зрения, однако после введения поправок на случайные события, рассеяние и ослабление излучения данные уже не совсем корректно считать пуассоновскими. Поэтому для повышения точности модели и расчетных соображений были предложены другие модели такие, как аппроксимации Пуассона [14], смещенный Пуассон [15], гауссовские модели.

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

(2.7) представляет функцию правдоподобия объекта f , т.е. задача

состоит в такой оценке объекта f , которая обеспечивает макси-

мальное значение L( ).

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

приближается к истинному изображению ( E( f ) ftrue ). Доказано, что эти оценки обеспечивают минимальную дисперсию и соз-

67

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

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

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

4.2. Алгоритм максимизации ожидания максимального правдоподобия

Метод максимизации математического ожидания максимально-

го правдоподобия (англ. maximum likelihood expectation maximization (MLEM)) или, короче, метод максимального правдоподобия (МП) является численным методом определения оценки максимального правдоподобия. Алгоритм MLEM , начиная с его первого применения в области реконструкции изображений в 1982 г. [16], остается основой для наиболее популярного статистического метода реконструкции изображений и является фундаментом для многих других методов.

Приложение алгоритма MLEM к реконструкции изображений в ПЭТ приводит к простому итерационному уравнению

68

				ˆ (n 1)		p
		f j(n 1)		f j	Hij
						i	,	(2.8)
				Hi j		ˆ (n)
					i	Hik fk
				i		k
ˆ (n 1)	следующая оценка вокселя j, основанная на текущей
где f j
ˆ	(n)	.
оценке f

Рис. 2.9. Диаграмма MLEM алгоритма. Начало с позиции первоначального угады-

ˆ (0)

вания ( f ), алгоритм итеративно выбирает новую оценку изображения, основанную на измеренных проекциях p (адаптировано из [3])

Опишем этот алгоритм более подробно с помощью рис. 2.9. Итерационный процесс начинается с "угадывания" начального

ˆ (0)	, показанного в верхнем левом
приближения к изображению f

углу рисунка и присутствующего как знаменатель в уравнении (2.8). За начальное приближение обычно берется равномерное распределение источников эмиссии по всему объему. Первый шаг (1) заключается в прямом проецировании этого изображения в пространство проекций. Затем (2) эти проекции сравниваются с измеренными проекциями p и формируется мультипликативный корректирующий фактор для каждой проекции. Далее (3) данный фак-

69

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

В процессе работы алгоритма после нескольких первых итераций в создаваемом изображении появляется низкочастотный компонент. При приближении к решению максимального правдоподобия в реконструируемом изображении появляется все больше высокочастотных проявлений, что существенно увеличивает вариацию в реконструкции. Как отмечалось выше, эту вариацию часто уменьшают (ценой увеличивающегося смещения), применяя раннюю остановку алгоритма или постсглаживание реконструкции. Сходимость алгоритма MLEM зависит от изображения, но на практике обычно для достижения решения требуется от 20 до 50 итераций. Так как MLEM требует на каждой итерации одного прямого и одного обратного проецирования, то общее время процессинга значительно больше, чем для метода фильтрованных обратных проекций, однако алгоритм MLEM позволяет получать более точную реконструкцию.

4.3. Алгоритм максимизации ожидания упорядоченных подмножеств

Алгоритм ожидания упорядоченных подмножеств (англ. ordered subsets expectation maximization (OSEM)) был внедрен в 1994 г. [17]

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

	ˆ (n)		p
f j(n 1)	f j	Hij
			i	,	(2.8)
			ˆ (n)
	Hi j i Sb
			Hik fk
	i Sb		k
		70