4 курс / Лучевая диагностика / Введение_в_комп_рентг_и_нейтронную_томографию

Рис. 64. Томографическая реконструкция в коническом пучке

Наиболее часто в случае траектории в виде окружности в настоящее время используется алгоритм Фельдкампа [9]. Его сущность заключается в преобразовании проекционных изображений согласно формуле (19):

где D – расстояние между источником и осью z, β – угловое положение источника по отношению к оси z, R(p,ζ,β) – проекционные данные.

Выражения (20) определяют значения ζ и p. Преобразование координат описывается выражениями (21).

Принятые в выражениях (19), (20) и (21) обозначения представлены на рисунке 64.

Преобразование Фельдкампа состоит из 3 этапов:

1. проекционные данные умножаются на весовой коэффициент

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

Алгоритм Фельдкампа можно применять для реконструкции при достаточно малом растворе конуса (<10°). При этом сравнительно точную реконструкцию можно получить, если изменения структуры объекта вдоль оси вращения достаточно медленные.

Точная реконструкция возможна, если удовлетворяется условие достаточности Туй-Смита, т.е. если все плоскости, пересекающие объект, также пересекают траекторию источника хотя бы раз.

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

Исследование практической возможности томографической реконструкции данных, полученных в коническом пучке, с помощью метода максимума функции правдоподобия (ММФП) показало, что за разумное время (несколько часов) с помощью ПК на основе процессора Pentium 4 можно восстанавливать объекты, состоящие не более чем из (60)3–(100)3 элементов. Такая грубая дискретизация объекта приводит к тому, что линейное ослабление в пределах одного элемента уже нельзя считать постоянным.

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

5.6. Реконструкция при использовании обратно рассеянного излучения

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

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

Процедура реконструкции в случае коллимированного и неколлимированного рентгеновского рассеянного излучения описана, например, в работах [93, 94].

Согласно [94], в случае неколлимированного излучения задача реконструкции сводится к решению линейной системы уравнений с неизвестными линейными коэффициентами рассеяния μС(i,j):

где: i,j – номер элемента в реконструируемой плоскости, N0 – исходное количество фотонов в пучке; k – отношение полного коэффициента ослабления к коэффициенту комптоновского рассеяния; α(x,y,θ) – доля зарегистрированных фотонов от всех рассеянных на векторе (x,y,θ) (рис. 11), Nd(x,y,θ) – число зарегистрированных фотонов.

5.7. Реконструкция при использовании фазового контраста

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

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

Контрольные вопросы к главе 5

1.От чего зависит качество томографического изображения?

2.Классификация методов и алгоритмов.

3.Суть метода регуляризации Тихонова (ознакомиться самостоятельно).

4.Формулировка теоремы о проекциях.

5.Достоинства и недостатки метода обратных проекций.

6.Алгоритм «свертки и обратного проецирования».

7.В чем заключается некорректность обратного преобразования Радона?

8.Достоинства и недостатки метода максимального правдоподобия.

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

10.Реконструкция при использовании обратно рассеянного излучения.

11.Методы реконструкции при использовании фазового контраста (ознакомиться самостоятельно).

Глава 6. Процесс измерения

В рассуждениях, приведенных выше, рассматривается связь между удельным поглощением в объекте и линейным интегралом поглощения (проекционными данными) p(h,θ). На самом деле детектор регистрирует интеграл от интенсивности прошедшего излучения Iθk:

Здесь I0(h) обозначает интенсивность падающего излучения, k – номер элемента детектора с координатой центра hk, а – это ширина отдельного элемента, которая считается постоянной. Если пренебрегать конечностью ширины элемента детектора Δ, т.е. считать, что на масштабе величины, входящие в задачу, не меняются, то интеграл (26) равен:

Это позволяет выразить проекционные данные p(hk,θ) через интенсивности излучения:

Если теперь считать, что интенсивность падающего излучения I0(hk)=I0 не зависит от h, то соотношение приобретает вид:

В более общем случае можно просто измерить излучение источника, не ослабленное объектом:

которое, естественно, уже не зависит от угла θ, и выразить проекционные данные как логарифм отношения интенсивностей без объекта и с объектом:

Соотношение (31) является приближенным и обычно используется для того, чтобы из измеренных величин Iθk получить проекционные данные p(hk,θ) для последующей реконструкции объекта. При его выводе не учитывались эффекты, связанные с конечностью ширины элемента детектора . Более того, предполагалось, что падающий пучок монохроматичен, т.е. не учитывалась зависимость поглощения и интенсивности падающего на объект излучения от энергии. Кроме того, падающий пучок считался плоско-параллельным, хотя в реальных измерениях всегда существует конечная расходимость пучка. Все эти эффекты приводят к тому, что вместо формулы (31) следует использовать более сложные выражения.

Для определения поправки к (31), связанной с конечностью ширины элемента детектора в случае одномерного детектора, воспользуемся разложением p(h,θ) в ряд Тэйлора по h в окрестности точки h=hk. Сначала предположим, что I0(h)=I0=const, и рассмот-

Подставляя (32) в (26), получим:

В пределе 2 p' (hk , ) 0 возвращаемся к выражению (31).

обозначениях выражение (33) можно переписать в виде:

Теперь рассмотрим первую поправку по малому параметру α.

Из (35) видно, что при любой зависимости поглощения от координаты вида (32) формула (31) дает заниженное значение для p(hk,θ).

Теперь рассмотрим квадратичную зависимость поглощения от координаты:

78 Глава 6

и предположим, что интенсивность падающего пучка тоже немного меняется в пределах элемента детектора:

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

В применяемом нами методе максимального правдоподобия в качестве исходных данных используются не значения поглощения в центре элемента детектора p(hk,θ), как это принято в других алгоритмах, а средние значения μk по данному элементу:

Подставив в этот интеграл выражение (36), получим:

Таким образом, мы получили соотношение между ηk и μk:

которое легко обобщается на случай реального двумерного детектора:

где: a p(hk , ) и I0 (hk ) .

2

Следует заметить, что интегрирование (40) по элементу детектора можно представить в виде свертки, что в пространстве Фурьеобразов эквивалентно умножению на функцию:

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

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

В томографических измерениях число проекций должно соответствовать числу лучей в проекции. Пусть при томографической съемке в диапазоне (0–180)° число проекций равно M, а число лучей N. Тогда инкремент в Фурье-пространстве между последующими проекциями δ составляет δ = π/M.

Согласно теореме Найквиста, максимальная частота, измеряемая в каждой проекции, равна:

где t – расстояние между соседними лучами.

Максимальная частота ωmax определяет радиус области измерений в частотном пространстве (рис. 65), а число лучей N в координатном пространстве определяет количество измеряемых точек на