CC BY
13
УДК 519.652
ПРЯМАЯ ЗАДАЧА ФОРМИРОВАНИЯ МНОГОСПЕКТРАЛЬНЫХ ИЗОБРАЖЕНИЙ В ЭМИССИОННОЙ ТОМОГРАФИИ И ОДИН ЧАСТНЫЙ СЛУЧАЙ ЕЕ ОБРАЩЕНИЯ
Иван Гаврилович Казанцев
Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, Новосибирск, пр. Академика Лаврентьева, 6, доктор физико-математических наук, доцент, старший научный сотрудник лаборатории обработки изображений, тел. (383)330-73-32, e-mail: [email protected]
В работе показано, что прямая задача в многоспектральной эмиссионной томографии может быть упрощена в предположении постоянного коэффициента ослабления среды. В этом случае она сводится к послойному лучевому преобразованию, к которому применима известная формула обращения.
Ключевые слова: многоспектральные изображения, томография.
A FORWARD PROBLEM OF MULTISPECTRAL IMAGES FORMATION IN EMISSION TOMOGRAPHY AND PARTIAL CASE OF ITS INVERSION
Ivan G. Kazantsev
Institute of Computational Mathematics and Mathematical Geophysics SB RAS, 630090, Russia, Novosibirsk, 6 pr. Akademika. Lavrentjeva, Ph. D., Senior scientific researcher of Image Processing Laboratory, tel. (383)330-73-32, e-mail: [email protected]
In this work it is shown that a forward problem of formation of multispectral images in emission tomography can be simplified provided the media has a constant linear attenuation coefficient. In this case the problem is reduced to the slice by slice X-ray transform and the known inversion formula can be applied.
Key words: Multispectral images, tomography.
В работе рассматривается задача восстановления изображений по проекциям, регистрируемым детекторами высокого спектрального разрешения. Идеализированная интегральная модель позитронной эмиссионной томографии (ПЭТ) на не рассеянных (первичных) фотонах хорошо известна. Для известного распределения активности изотопа f(x,y,z) с носителем D(f) внутри среды с линейным коэффициентом ослабления ^(x,y,z) и носителем D(^) и детекторов А и B малых размеров на поверхности просвечиваемого тела, модель формирования данных имеет вид:
в в
р(А,В) = Qxpf- z\)dl']\f(x,y, z)dl С)
А А
Здесь dl, dl' - элементы длины на линии AB. Задача состоит в реконструкции f при известной среде ц по данным p(A,B), регистрируемым большим множеством пар детекторов (A,B). Данные детекторов (1), после коррекции на
ослабление, интерпретируются как интегралы от f по отрезкам АВ распространения первичных фотонов с энергией Е=511 кэВ и задача (1) сводится к классической томографии. На рис. 1 (слева) изображена модель эмиссионной томографии на паре прямых фотонов (и,у), разлетающихся из точки аннигиляции С в результате столкновения позитрона, испускаемого изотопом активности £ с одним из свободных электронов среды
Рис. 1. Модели позитронной эмиссионной томографии:
Слева: геометрическая модель на первичных фотонах (u,v); справа: модель комптоновского рассеяния фотона v в точке S под углом 0.
Однако в детекторы попадают и фотоны с меньшими энергиями E'<E, претерпевшие комптоновское рассеяние (рис. 1 справа), т.е. некогерентное отклонение на угол 0. Угол рассеяния 0 может быть вычислен из соотношения:
= 511 (feK) (2)
2 — cos О
С улучшением спектрального разрешения детекторов появляется возможность настраивания ПЭТ-сканеров на регистрацию фотонов, рассеянных под определенным углом 0, с последующим использованием информации в этих данных в дополнение к первичным фотонам. В работе [1] с соавторами получена (и исследована в сравнении с проекциями, генерированными статистически методом Монте-Карло) формула вычисления потока фотонов, рассеянных с углом 0, регистрируемых детектором B при условии, что в A попали первичные фотоны:
в
С
&(Л,В)= ¡¿<рсоБ<рсоз(6>-<р) 1) с-е А * ■ \/(у/,<р,г)Ф (3)
о 0 сгс 30 0
Здесь Б - точка, где происходит комптоновское рассеяние, ёаС /ё О -дифференциальное поперечное сечение рассеяния. Знак штриха означает значение величины после рассеяния. Сферические координаты (у,ф,г) обозначены на Рис. 2 слева.
В статье рассматривается сведение интегрального преобразования (3) к более простой модели, известной как преобразование Радона (лучевое) с послойной сверткой. Ядро свертки зависит от положения - координаты Ъ слоя среды, подвергающегося свертке (рис. 2 справа). Нами доказано ранее [2], что при наличии полных данных, это преобразование имеет обратное, которое
можно применять в виде конструктивного алгоритма реконструкции. Рассмотрим преобразование Радона с послойной сверткой подробнее.
£fa'W,z)*k(x-x',y-y',z) , Р^М(х,у)
УА
4 \ Z1 1 2 Y
. ■ ■ ■ ■ 5» I
Л''' 1 ;...............\1..................
Рис. 2. Иллюстрации математических обозначений:
Слева: Геометрическая модель комптоновского рассеяния, лежащая в основе уравнения (3). Справа: Сечения объекта - версии £ повернутой с полярными параметрами а и в, с координатами Ъ\ и Ъ2 . Все сечения подвергаются свертке и суммируются параллельно оси Ъ.
На рис. 2 (справа) иллюстрируется определение лучевого преобразования Радона с послойной сверткой в общем виде:
Pl,M>s2)= -sx,s2-s2,z)ds\ds2dz
(4)
Здесь f*01'^ - повернутая версия f с углами поворота аир. Можно показать, что в случае /u(E,S) = const в рассматриваемой нами прямой задаче спектральной ПЭТ уравнение (3) может быть преобразовано к виду (4), где ядро свертки принимает вид:
Jul z(zcos0 + psm0) J
\AB\ju
4л-1 AB | p2(z2+p1)
2 \3/2
e
z-+p-
(z+yotan(6'/2))
(5)
Ядро является радиально-симметричным и зависит от ъ - координаты слоя, подвергающегося свертке. Удобно представить уравнение (4) в частотной области:
(6)
Алгоритм послойного обратного проецирования с фильтром обращения, также зависящим от Ъ^ состоит из шагов:
Шаг 1. Фильтрация проекции Р^ р(кх,к2,2) с помощью регуляризации обратного фильтра Щ1 в частотной области для каждого слоя с позицией Я.
-1
»6»
а,(5х
(H2e(kl,k2,z) + A\k\2)
}
(7)
Шаг 2. Вычисление интегрального изображения в виде сум-
мирования обратных проекций Ъва р по всем направлениям св(х? у.(«, (3) е .V2
с0 (х, у, г) = \\Ъваф (х, у, г)дад/3 (8)
Шаг 3. Интегральное изображение подвергается фильтрации типичным для томографии фильтром | к | в частотной области
Г{к^к2,к3)} = С(?(к1,к2,к3)х\ к\ (9)
Нижеследующая теорема [2] утверждает, что реконструкция по полному набору спектральных проекций не зависит от угла рассеяния 0.
Теорема: Интегральные изображения (в области преобразования
Фурье) асимптотически тождественны на высоких частотах к , и тождественны интегральному изображению С<0=<у>(к), С0 (к) полученному по не рассеянным проекционным данным.
Для проверки полученных формул проведены вычислительные эксперименты - с полным углом обзора и с неполными данными и по углу и по пространству. Для экспериментов с полным углом обзора выбран фантом размером 256 . Он состоит из 7 больших и 11 малых сфер (кис. 3 слева), служащих носителем функции активности f=1,помещенных в воду с постоянным коэффициентом ослабления ц=0.096см-1.
Рис. 3. Вычислительный эксперимент с полными данными:
"5
Слева: Цифровой фантом размером 256 . Справа: три сечения фантома (верхний ряд); те же сечения реконструкции по 1000 проекциям.
Сгенерированы 1000 проекций в случайно выбранных направлениях по формуле (3) для угла сканирования 0=300. Восстановление произведено по алгоритму (7)-(9) с параметром регуляризации X = 0.01. На рис. 3 справа - три сечения фантома, верхний ряд - тестовый объект, нижний - реконструкция.
Вычисления по неполным данным иллюстрируются на рис. 4. Рассматривается система из двух линейных панелей с детекторами A и B. Каждая линейка содержит 72 детектора, настроенные на регистрацию фотонов обладающих
о
энергией соответствующей углу рассеяния 0 = 30 . Данные сгенерированы по формуле (3), ослабление постоянно ц=0.096см-1 . Применение алгоритма обратного проецирования с послойной регуляризирующей фильтрацией дает реконструкции с множеством артефактов вследствие существенной неполноты данных.
А
Рис. 4. Вычислительный эксперимент с неполными данными:
Слева: двумерная геометрия ПЭТ сканера с двумя линейками детекторов A и B. Центр: тестовое изображение 72x72 с областями различных значений активности f; Справа: результат реконструкции по 50 итерациям алгебраического метода ART-TV с 20 итерациями метода полной вариации.
Поэтому применен хорошо адаптируемый к различным геометриям просвечивания и теоретически точный алгебраический итерационный метод ART-TV [3] с параметром релаксации 0.03, в котором на каждой итерации осуществляется минимизация полной вариации градиента реконструируемого изображения. Эксперименты позволяют сделать вывод, что многоспектральная информация, содержащаяся в данных рассеяния, может быть успешно использована наравне с традиционным монохроматическим излучением.
Работа частично поддержана грантом РФФИ 16-07-00066.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Казанцев И. Г., Яровенко И. П., Прохоров И. В. Моделирование процесса измерения комптоновского рассеяния в позитронной эмиссионной томографии // Вычислительные технологии. - 2011. - Т. 16, № 6. -С. 27-37.
2. Kazantsev I. G., Klukowska J., Herman G. T., Cernetic L. Fully three-dimensional defocus-gradient backprojection in cryoelectron microscopy // Ultramicroscopy. - 2010. - V. 100. -P.1128-1142.
3. Sidky E.Y., Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization // Physics in Medicine and Biology. - 2008. - V. 53. -P.4777-4807.
© И. Г. Казанцев, 2016
