Аналитическое и статистическое моделирование формирования изображений рассеянного излучения в эмиссионной томографии Текст научной статьи по специальности «Физика»

УДК 621.391 И.Г. Казанцев

ИВМиМГ СО РАН, Новосибирск И.П. Яровенко, И.П. Прохоров ИПМ ДВО РАН, Владивосток

АНАЛИТИЧЕСКОЕ И СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ФОРМИРОВАНИЯ ИЗОБРАЖЕНИЙ РАССЕЯННОГО ИЗЛУЧЕНИЯ В ЭМИССИОННОЙ ТОМОГРАФИИ

Рассмотрена задача моделирования изображений рассеянных фотонов в эмиссионной томографии. Выведенные интегральные уравнения требуют сравнительного статистического эксперимента, которое проведено методами Монте-Карло. Получено хорошее совпадение изображений, генерированных статистически и с помощью формул. Работа поддержана грантом РФФИ 10-07-00131-а.

I.G. Kazantsev

Institute of Computational Mathematics and Mathematical Geophysics SB RAS

6, prospect Akademika Lavrentjeva, Novosibirsk, 630090, Russian Federation

I.P. Yarovenko, I.VProkhorov

Institute of Applied Mathematics FEB RAS

7 Radio Street, Vladivostok, 690041, Russian Federation

ANALYTICAL AND STATISTICAL MODELING OF SCATTER RADIATION IMAGE FORMATION IN EMISSION TOMOGRAPHY

The mathematical model of single scatter image formation in emission tomography is investigated. Integral equations are derived under idealized assumptions and therefore they need physical or statistical verification. Statistical validation of the analytical approach using Monte Carlo simulation techniques is presented. The numerical experiments were performed. Good fit of formula-based and statistically generated profiles of scatter images is observed.

Идеализированная интегральная модель эмиссионной томографии на не рассеянных (прямых) фотонах хорошо известна. Для известного распределения активности (изотопа) f(x,y,z) внутри среды с коэффициентом ослабления ^(x,y,z) и детекторов малых размеров А и B на поверхности просвечиваемого тела, модель формирования данных имеет вид:

В В

РАВ =exp[-j//(.x',/,z')<i/]- j'f(x,y,z)dl (1)

A A

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

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

В практических задачах оценки однократного рассеяния в сканерах, широко используется аппроксимация [1], являющаяся упрощенной версией уравнения переноса [2],[3]. Это приближение оценивает поток однократно рассеянных (по Комптону) пар фотонов, зарегистрированных парой детекторов А и В как интеграл по всему объему У8 возможных событий рассеяния:

_)Л_д°^[1А+1 в^

S1AB = S\dVs (

Ап|AS||BS|

(Тс дО

(2)

где

О D О ODD

!Л =£as8'bs^P[-(J/^+ \M'dl)]\fdl, IB =s'ASsBSexp[-(J/£//+ J//V//)]jfd/ (3)

ASA ASS

Здесь gas и aBS - геометрические поперечные сечения детекторов A и B, f -активность источников, ^ = ^(E,S) - линейный коэффициент ослабления, зависящий от энергии и положения точки рассеяния, sAS и sBS - величины, связанные с эффективностью детекторов A и B, daC /d Q - дифференциальное поперечное сечение рассеяния. Знак штриха означает значение величины после рассеяния. Уравнение (2) симметрично относительно A и B, и прямые фотоны регистрируются в обоих детекторах A и B. Модифицируем это уравнение в предположении, что регистрируется пара фотонов, один из которых попал в детектор A нерассеянным, а детектор B зафиксировал рассеянный фотон и определил его энергию. Тогда можно считать, что выражение f в (3) равно нулю. Для удобства изложения обозначим ц дас

ас 60

£as£'bs ехр[-(\jjdl + J/Л//)] |/б// , (4)

вид

где (х,у,7) декартовы координаты точки рассеяния Б. Тогда (2) принимает

(5)

Это уравнение будет использовано для вывода интегральных уравнений однократного комптоновского рассеяния с определенным углом 0.

S

A

S

A

Аппроксимация (2)-(5) применима к объектам различной формы или объема УБ. Однако для идеальных детекторов уравнение (5) допускает

упрощение благодаря свойствам поверхностей, на которых сосредоточены точки рассеяния с одинаковым углом. Тогда объем У8 описывается точно и интегрирование в (5) приобретает конкретные пределы.

Предположим, что детекторы А и В - круги малого радиуса 5, 5 << 1 (Рис. 1 справа). Детектор А регистрирует прямые, детектор В - рассеянные фотоны. Из геометрии однократного рассеяния можно видеть, что все точки рассеяния Б (с энергией фотона Е', или углом рассеяния 0) лежат на поверхности тела вращения (вокруг оси Т) дуги окружности АБВ радиуса

Я =| АВ11{2%т 6). (6)

Выберем некоторый угол 0 < Т < п/2. Тогда все события, состоящие в регистрации в В рассеянных фотонов с углами, меньшими или равными Т, при условии, что прямой фотон из пары, испускаемой из точки С, зарегистрирован в А, составляют в сферических координатах (\|/,ср,г) множество (тело вращения)

К = {(уу, ф,г)\у/е [0,2л-), <р е [0, Т\ г е [0,1АБ |]}. (7)

Выберем некоторый угол 0 <0 < п/2. Тогда все события, состоящие в рассеянии фотонов с углом, точно равным 0, составляют поверхность тела вращения

Ъв ={(у/,<р,г) I ^ е [0,2*), <р е [О, Г], г = [О, IАБ |]}. (8)

Хорды |А8| и |8В| вычисляются по формулам:

| АБ |=| АВ | $хп(6-(р)!ъхпб , | БВ |=| АВ | зт (р/ъхпб . (9)

Будем искать полное рассеяние с порогом Т, учитывая формулу (5), в виде:

Б^ = 1ш1----1---- \дА [(^ )дВ. (10)

г (лб I А) \ 1 7 47

Для малых 5 геометрические поперечные сечения детекторов А и В для падающих лучей АБ и БВ являются приближенно эллипсами с площадями

(7А5 ~ {п82/4)со8^> , ст5В х (7г32 / А)со$>(6 - (р) . (11)

Тогда (10) принимает вид

~ (р])МГм(х, у, 2). (12)

г т Ап\АБ\2\ВБ \2

Заменим переменные в (12) с декартовых (х,у,г) на криволинейные (\|/,ф,0):

х = Х{ц/,(р,&) =| АБ | бш (рсоъцг,

у = У{ц/,(р,в) =| АБ | $,т<р$,тц/, (13)

г = 2(ц/,(р,9) =| АБ | со$,<р,

Вычислим элементарный объем интегрирования с1УТ = йхйуй2= det(У)Jy/J<7*i0, где I - якобиан преобразования (13) и его определитель может быть вычислен:

= (14)

вт в |АВ|

Тогда интеграл (12) приводится к искомому виду (положим е=1 для краткости)

5 В

сав 1 ь/г, т л'Ь м(¥,(Р,\АБ\)дас ^

^ =——— 1^ |^СО8^СО8(0-^) 1^--------------------------------------• \Я¥,<Р,г)^ . (15)

Аж | АВ | * I I ос дп *

выводится интегральное уравнение рассеяния с углом

Аналогично рассеяния 0:

sin2# вг, cos<pcos(0-(p)2!r ju(y/,<p,\AS\) да1

rtAB _ ^в ~

4л | АВ\‘

■ U(p-j 011

sm(psm(9-(p)

07

дО.

'SAS^BS

-(J .иЯ+J ji’dl) \AS\

■e A s ■ ^f{y/,(p,r)dr. (16)

В данной работе исследуется упрощенная оценка ^3) рассеяния с

о

углом рассеяния е=60 , в виде

5 В

М(.Ч>,<Р,\Л8\)8стс и?

О 2 к

£e(A,B)= ^d<pcos<pcos(9-<р) jёц/‘-

дП

\f(V,<P,r)dr • (I7)

Для проверки формулы (17) методами статистического моделирования была посчитана величина

^ = И2/4)24±дй(Лй), (18)

где 5 - диаметр детектора, Д0 - некоторый диапазон изменения угла рассеяния, который при моделировании траекторий обеспечивает попадание в некоторую фиксированную точку. Вычислительные эксперименты проведены для цилиндрического фантома радиусом 8 см и заполненного водой (рис. 2). Внутрь цилиндра помещены четыре стационарно излучающих сферических источника, заполненные изотопом с единичной активностью: No. 1 в центре и No. 2, 3, 4, смещенные на 4см от центра. Коэффициент ослабления воды ^ принят одинаковым для всех энергий ^=^’=0.096см-1 . Линейки A и B детекторов лежат в плоскости YOZ, и проекция ^0 (A,B) вычисляется для пар точек (A,B) с одинаковой координатой Y. Параметры теста: 0=60°, 5=0.16 см, s = s’=1. При моделировании траекторий отслеживалось до 10 актов рассеяния. Для метода Монте-Карло данные о сечениях взаимодействия брались из таблиц Хабла [4]. Моделировалось 1011 траекторий (время расчета на компьютере Pentium 3.2 GHz около 5 суток). Для генерации случайных чисел использовался датчик, описанный в [5]. Результаты численных экспериментов (рис. 3) показывают хорошее соответствие между проекциями, полученными по формуле (17) и статистическим моделированием (18).

0

0

0

с

0

0

0

Рис. 2

Рис. 3

Алгоритм состоит из следующих шагов:

1. Розыгрыш точки аннигиляции фотона. Разыгрываем точку рождения позитрона внутри источника активности. Ее же считаем точкой аннигиляции.

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

3. Розыгрыш свободного пробега фотона в веществе. Учитывая экспоненциальный закон уменьшения интенсивности потока фотонов с расстоянием, найдем длину свободного пробега в веществе. Если фотон в результате пробега покинул область, проверяется: не попал ли он в детектор А. Если попал, то запоминаем номер детектора и отслеживаем направление -ю, если же фотон не попал в детектор А, то перестаем отслеживать данную траекторию и переходим к пункту 1. В случае, когда фотон остался внутри области, разыгрываем тип взаимодействия фотона с веществом.

4. Розыгрыш взаимодействия фотона с веществом. Для рассматриваемого диапазона энергий, среди всех видов взаимодействия излучения с веществом преобладают фотоэлектрическое поглощение, комптоновское рассеяние и рассеяние по Рэлею. Для определения типа взаимодействия излучения с веществом вычисляются вероятности поглощения фотона, рассеяния по Комптону и рассеяния по Рэлею: Ра, Рс, Рг. Если а < Ра, то фотон поглотился, и такая траектория дальше не отслеживается. Если а > Ра, то фотон рассеется и необходимо определить тип рассеяния. При а < Рг, фотон рассеется по Рэлею, в противном случае фотон рассеется по закону Комптона.

5. Рассеяние фотона по Рэлею. Если фотон рассеялся по Рэлею, то разыгрываем новое направление ю’, используя функцию рассеяния Рэлея.

6. Рассеяние фотона по Комптону. При рассеянии фотона по Комптону изменяется не только направление, но и энергия фотона. Вероятность рассеянному фотону иметь энергию от Е до Е' в зависимости от значения

случайного числа а определяется из уравнения \^^dE = a , где dac /dE -

Е dE

нормированное на единицу дифференциальное сечение комптоновского рассеяния в интервал энергий от E до E+dE, определяемое формулой Кляйна-Нишины [3]. Для ускорения решения этого уравнения использовалась заранее насчитанная двумерная таблица, содержащая решения уравнения с заданной точностью на равномерных сетках по E и а. Найденное значение энергии рассеянного фотона позволяет определить косинус угла рассеяния: у = \+Е0/Е-Е0/Е' , где Е0 - энергия покоя электрона. Далее, разыгрывая равновероятно в диапазоне от 0 до 2п азимутальный угол рассеяния, находим новое направление движения фотона после рассеяния.

7. Возврат и выход. Возврат на шаг 3 к розыгрышу свободного пробега и повторение этого пока либо фотон попадает в детектор или покидает объем, или поглощается.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Watson C.C. “New, faster, image-based scatter correction for 3D PET,” IEEE Transactions on Nuclear Science. -2000. - Vol. 47. -P. 1587-1594.

2. Аниконов Д.С., Ковтанюк А.Е., Прохоров И.В. “Использование уравнения переноса в томографии “, M: Логос, 2000.

3. Фано У, Спенсер Л., Бергер М. “Перенос гамма излучения” М.: Госатомиздат, 1963.

4. Hubbell J.H., Seltzer S.M. “Tables of X-Ray Mass Attenuation Coefficients and Mass Energy-Absorption Coefficients 1 Kev to 20 Mev for Elements Z = 1 to 92 and 48 Additional Substances of Dosimetric Interest,” NISTIR 5632, 1995.

5. L'Ecuyer, P. and Cote, S. “Implementing a Random Number Package with Splitting Facilities.” ACM Transactions on Mathematical Software, 17: 98-111 (1991).

© И.Г. Казанцев, И.П. Яровенко, И.П. Прохоров, 2011