Метод Монте-Карло для транспорта фотонов

Моделирование распространения фотонов с помощью методов Монте-Карло — это гибкий, но строгий подход к моделированию переноса фотонов. В методе локальные правила переноса фотонов выражаются в виде распределений вероятностей, которые описывают размер шага движения фотона между местами взаимодействия фотона с веществом и углы отклонения траектории фотона при возникновении события рассеяния. Это эквивалентно аналитическому моделированию переноса фотонов с помощью уравнения переноса излучения (RTE), которое описывает движение фотонов с помощью дифференциального уравнения. Однако решения RTE в закрытой форме часто невозможны; для некоторых геометрий диффузионное приближение можно использовать для упрощения RTE, хотя это, в свою очередь, вносит множество неточностей, особенно вблизи источников и границ. Напротив, моделирование методом Монте-Карло можно сделать сколь угодно точным, увеличив количество отслеживаемых фотонов. Например, посмотрите фильм, где моделирование методом Монте-Карло карандашного луча , падающего на полубесконечную среду, моделирует как начальный поток баллистических фотонов, так и более позднее диффузное распространение.

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

Моделирование методом Монте-Карло падающего карандашного луча на полубесконечную рассеивающую среду.

Биомедицинские применения методов Монте-Карло

Биомедицинская визуализация

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

• Фотоакустическая томография . При ПАТ диффузный лазерный свет поглощается, что приводит к локальному повышению температуры. Это локальное изменение температуры, в свою очередь, генерирует ультразвуковые волны посредством термоупругого расширения, которые обнаруживаются ультразвуковым преобразователем. На практике варьируются различные параметры установки (например, длина волны света, числовая апертура датчика), и в результате моделирование Монте-Карло является ценным инструментом для прогнозирования реакции ткани до применения экспериментальных методов.
• Диффузная оптическая томография DOT — это метод визуализации, в котором используется ряд источников и детекторов ближнего инфракрасного диапазона для измерения оптических свойств биологических тканей. Можно измерить различные контрасты, включая поглощение окси- и дезокси-гемоглобина (для функциональной нейровизуализации или обнаружения рака) и концентрацию флуоресцентных зондов. Чтобы восстановить изображение, необходимо знать, каким образом свет прошел от данного источника к данному детектору и как измерение зависит от распределения и изменений оптических свойств (так называемая прямая модель). Из-за сильно рассеивающей природы биологической ткани такие пути сложны, а функции чувствительности размыты. Прямая модель часто создается с использованием методов Монте-Карло.

Лучевая терапия

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

Фотодинамическая терапия

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

Реализация транспорта фотонов в рассеивающей среде

Представлена ​​модель фотонного метода Монте-Карло в однородной бесконечной среде. Однако модель легко расширить для многослойных носителей. Для неоднородной среды необходимо учитывать границы. Кроме того, для полубесконечной среды (в которой фотоны считаются потерянными, если они выходят за верхнюю границу) необходимо принять особое внимание. Для получения дополнительной информации перейдите по ссылкам внизу страницы. Мы будем решать задачу, используя бесконечно малый точечный источник (аналитически представленный как дельта-функция Дирака в пространстве и времени). Реакция на произвольную геометрию источника может быть построена с использованием метода функций Грина (или свертки , если существует достаточная пространственная симметрия). Необходимыми параметрами являются коэффициент поглощения , коэффициент рассеяния и фазовая функция рассеяния. (Если рассматриваются границы, необходимо также указать показатель преломления для каждой среды.) Реакции с временным разрешением находятся путем отслеживания общего прошедшего времени полета фотона с использованием длины оптического пути . Реакции на источники с произвольными временными профилями затем можно смоделировать посредством свертки во времени.

В нашей упрощенной модели мы используем следующий метод уменьшения дисперсии, чтобы сократить время вычислений. Вместо распространения фотонов по отдельности мы создаем пакет фотонов с определенным весом (обычно инициализируемым как единица). Когда фотон взаимодействует в мутной среде, он отдает вес за счет поглощения, а оставшийся вес рассеивается в другие части среды. По пути может регистрироваться любое количество переменных, в зависимости от интересов конкретного приложения. Каждый фотонный пакет будет неоднократно подвергаться следующим пронумерованным шагам, пока не будет завершен, отражен или передан. Процесс изображен на схеме справа. Любое количество фотонных пакетов может быть запущено и смоделировано до тех пор, пока полученные смоделированные измерения не будут иметь желаемое соотношение сигнал/шум. Обратите внимание: поскольку моделирование Монте-Карло представляет собой статистический процесс, включающий случайные числа, мы будем использовать переменную ξ как псевдослучайное число для многих расчетов.

Схема моделирования потока фотонов в бесконечной рассеивающей и поглощающей среде с помощью моделирования Монте-Карло.

Шаг 1: Запуск фотонного пакета

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

${\displaystyle {\begin{aligned}x&=0\\{\text{Position: }}y&=0\\z&=0\\\\\mu _{x}&=0\\{\text{Direction cosines: }}\mu _{y}&=0\\\mu _{z}&=1\end{aligned}}}$

Шаг 2: Выбор размера шага и движение пакета фотонов

Размер шага s — это расстояние, которое пакет фотонов проходит между местами взаимодействия. Существует множество методов выбора размера шага. Ниже приведена базовая форма выбора размера шага фотона (полученная с использованием метода обратного распределения и закона Бера-Ламберта ), которую мы используем для нашей однородной модели:

${\displaystyle s=-{\frac {\ln \xi }{\mu _{t}}}}$

где – случайное число, – полный коэффициент взаимодействия (т. е. сумма коэффициентов поглощения и рассеяния). ${\displaystyle \xi }$ ${\displaystyle {\mu _{t}}}$

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

${\displaystyle {\begin{aligned}x&\leftarrow x+\mu _{x}s\\y&\leftarrow y+\mu _{y}s\\z&\leftarrow z+\mu _{z}s\end{aligned}}}$

Шаг 3: Поглощение и рассеяние

Часть веса фотона поглощается в каждом месте взаимодействия. Эта доля веса определяется следующим образом:

${\displaystyle \Delta W={\frac {\mu _{a}}{\mu _{t}}}W}$

где коэффициент поглощения. ${\displaystyle {\mu _{a}}}$

Массовую долю затем можно записать в виде массива, если распределение поглощения представляет интерес для конкретного исследования. Затем вес фотонного пакета необходимо обновить следующим образом:

${\displaystyle W\leftarrow W-\Delta W\,}$

После поглощения пакет фотонов рассеивается. Средневзвешенное значение косинуса угла рассеяния фотонов известно как анизотропия рассеяния ( g ), которая имеет значение от -1 до 1. Если оптическая анизотропия равна 0, это обычно указывает на то, что рассеяние изотропно. Если g приближается к значению 1, это указывает на то, что рассеяние происходит преимущественно в прямом направлении. Чтобы определить новое направление фотонного пакета (и, следовательно, косинусы направления фотонов), нам необходимо знать фазовую функцию рассеяния. Часто используется фазовая функция Хеньи-Гринштейна. Затем угол рассеяния θ определяется по следующей формуле.

${\displaystyle \cos \theta ={\begin{cases}{\frac {1}{2g}}\left[1+g^{2}-\left({\frac {1-g^{2}}{1-g+2g\xi }}\right)^{2}\right]&{\text{ if }}g\neq 0\\1-2\xi &{\text{ if }}g=0\end{cases}}}$

Обычно предполагается, что полярный угол φ равномерно распределен между 0 и . Исходя из этого предположения, можно установить: ${\displaystyle 2\pi }$

${\displaystyle \varphi =2\pi \xi {\frac {}{}}}$

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

${\displaystyle {\begin{aligned}\mu '_{x}&={\frac {\sin \theta (\mu _{x}\mu _{z}\cos \varphi -\mu _{y}\sin \varphi )}{\sqrt {1-\mu _{z}^{2}}}}+\mu _{x}\cos \theta \\\mu '_{y}&={\frac {\sin \theta (\mu _{y}\mu _{z}\cos \varphi +\mu _{x}\sin \varphi )}{\sqrt {1-\mu _{z}^{2}}}}+\mu _{y}\cos \theta \\\mu '_{z}&=-{\sqrt {1-\mu _{z}^{2}}}\sin \theta \cos \varphi +\mu _{z}\cos \theta \\\end{aligned}}}$

Для особого случая

${\displaystyle {\begin{aligned}\mu _{z}=1\end{aligned}}}$

использовать

${\displaystyle {\begin{aligned}\mu '_{x}&=\sin \theta \cos \varphi \\\mu '_{y}&=\sin \theta \sin \varphi \\\mu '_{z}&=\cos \theta \\\end{aligned}}}$

или

${\displaystyle {\begin{aligned}\mu _{z}=-1\end{aligned}}}$

использовать

${\displaystyle {\begin{aligned}\mu '_{x}&=\sin \theta \cos \varphi \\\mu '_{y}&=-\sin \theta \sin \varphi \\\mu '_{z}&=-\cos \theta \\\end{aligned}}}$

C-код:

/*********************** Индикатриса **********************Новые направляющие косинусы после рассеяния на угол тета, фи.* mux new=(sin(theta)*(mux*muz*cos(fi)-muy*sin(fi)))/sqrt(1-muz^2)+mux*cos(theta)* muy new=(sin(theta)*(muy*muz*cos(fi)+mux*sin(fi)))/sqrt(1-muz^2)+muy*cos(theta)* muz new= - sqrt(1-muz^2)*sin(theta)*cos(fi)+muz*cos(theta)*-------------------------------------------------- --------*Вход:* muxs,muys,muzs — направляющий косинус до столкновения* mutheta, fi - косинус полярного угла и азимутального угла*-------------------------------------------------- --------*Выход:* muxd,muyd,muzd - косинус направления после столкновения*-------------------------------------------------- --------*/void Indicatrix (двойные мультиплексоры, двойные муи, двойные музы, двойные мутета, двойные фи, двойные *muxd, двойные *muyd, двойные *muzd){ двойная косета = мутета; двойная синтета = sqrt(1.0-костета*костета); // грех(тета) двойной грех = грех (фи); двойной косфи = соз (фи); если (музс == 1.0) { *muxd = sintheta*cosfi; *muyd = sintheta*sinfi; *muzd = костета; } elseif (muzs == -1.0) { *muxd = sintheta*cosfi; *muyd = -sintheta*sinfi; *muzd = -костета; } еще { двойной номинал = sqrt(1.0-muzs*muzs); двойной muzcosfi = muzs*cosfi; *muxd = sintheta*(muxs*muzcosfi-muys*sinfi)/denom + muxs*costtheta; *muyd = sintheta*(muys*muzcosfi+muxs*sinfi)/denom + muys*costtheta; *muzd = -denom*sintheta*cosfi + muzs*costtheta; }}

Шаг 4: Фотонное прекращение

Если пакет фотонов подвергся множеству взаимодействий, для большинства приложений вес, оставшийся в пакете, не имеет большого значения. В результате необходимо определить средство терминации фотонных пакетов достаточно малого веса. Простой метод будет использовать порог, и если вес фотонного пакета ниже порога, пакет считается мертвым. Вышеупомянутый метод ограничен, поскольку он не экономит энергию. Чтобы сохранить общую энергию постоянной, часто используется метод русской рулетки для фотонов ниже определенного весового порога. Этот метод использует константу рулетки m , чтобы определить, выживет ли фотон. У пакета фотонов есть один шанс выжить в m , и в этом случае ему будет присвоен новый вес mW , где W — начальный вес (этот новый вес в среднем сохраняет энергию). Во всех остальных случаях вес фотона устанавливается равным 0, и фотон уничтожается. Это выражается математически ниже:

${\displaystyle W={\begin{cases}mW&\xi \leq 1/m\\0&\xi >1/m\end{cases}}}$

Графические процессоры (GPU) и быстрое моделирование переноса фотонов методом Монте-Карло

Моделирование миграции фотонов в мутных средах методом Монте-Карло представляет собой задачу с высокой степенью распараллеливания, в которой большое количество фотонов распространяется независимо, но по одинаковым правилам и разным последовательностям случайных чисел. Параллельная природа этого особого типа моделирования Монте-Карло делает его очень подходящим для выполнения на графическом процессоре (GPU). Выпуск программируемых графических процессоров положил начало такому развитию, и с 2008 года появилось несколько сообщений об использовании графических процессоров для высокоскоростного моделирования миграции фотонов методом Монте-Карло. ^{[1] [2] [3] [4]}

Этот базовый подход сам по себе можно распараллелить, используя несколько связанных вместе графических процессоров. Одним из примеров является «GPU Cluster MCML», который можно загрузить с веб-сайта авторов (Моделирование переноса света по методу Монте-Карло в многослойных мутных средах на основе кластеров GPU): http://bmp.hust.edu.cn/ GPU_Cluster/GPU_Cluster_MCML.HTM

Смотрите также

• Уравнение переноса излучения и теория диффузии для транспорта фотонов в биологической ткани
• Метод Монте-Карло
• Свертка для оптических откликов широкого луча в рассеивающих средах
• Методы Монте-Карло для электронного транспорта

Ссылки на другие ресурсы Монте-Карло

• Лаборатория оптической визуализации Вашингтонского университета в Сент-Луисе (MCML)
• Медицинский лазерный центр Орегона
• Исследование миграции фотонов методом Монте-Карло в Лундском университете, Швеция. Ускорение моделирования Монте-Карло с помощью графического процессора и масштабируемый метод Монте-Карло. Открытый исходный код для загрузки.
• Облачный метод Монте-Карло для переноса света в мутной рассеивающей среде. Инструмент можно бесплатно использовать в исследовательской и некоммерческой деятельности.
• Транспорт света в тканях как пример моделирования Монте-Карло (с исходным кодом C++).

Рекомендации

• Ван, Л.Х.; У Синь-И (2007). Биомедицинская оптика: принципы и визуализация . Уайли.
• Л.-Х. Ван; СЛ Жак; Л.-К. Чжэн (1995). «MCML — моделирование переноса света в многослойных тканях по методу Монте-Карло». Компьютерные методы и программы в биомедицине . 47 (2): 131–146. дои : 10.1016/0169-2607(95)01640-F. ПМИД  7587160.
• Л.-Х. Ван; СЛ Жак; Л.-К. Чжэн (1997). «Conv — свертка для реакции на пучок фотонов конечного диаметра, падающий на многослойные ткани» (PDF) . Компьютерные методы и программы в биомедицине . 54 (3): 141–150. дои : 10.1016/S0169-2607(97)00021-7. ПМИД  9421660.
• СЛ Жак; Л.-Х. Ван (1995). «Моделирование транспорта света в тканях методом Монте-Карло» (PDF) . В Эй Джей Уэлче; MJC ван Гемерт (ред.). Оптический тепловой отклик ткани, облученной лазером . Нью-Йорк: Пленум Пресс. стр. 73–100.
• Л.-Х. Ван; С.Л. Жак (1994). «Оптимизированные радиальные и угловые положения при моделировании Монте-Карло» (PDF) . Медицинская физика . 21 (7): 1081–1083. Бибкод : 1994MedPh..21.1081W. дои : 10.1118/1.597351. ПМИД  7968840.

Встроенные ссылки

1. Э. Алерстам; Т. Свенссон; С. Андерссон-Энгельс (2008). «Параллельные вычисления с графическими процессорами для высокоскоростного моделирования миграции фотонов по методу Монте-Карло» (PDF) . Дж. Биомед. Опц . 13 (6): 060504. Бибкод : 2008JBO....13f0504A. дои : 10.1117/1.3041496 . ПМИД  19123645.
2. К. Фанг; Д.А. Боас (2009). «Моделирование миграции фотонов в мутной трехмерной среде методом Монте-Карло с ускорением с помощью графических процессоров». Опция Выражать . 17 (22): 20178–20190. Бибкод : 2009OExpr..1720178F. дои : 10.1364/oe.17.020178. ПМЦ 2863034 . ПМИД  19997242. 
3. Н. Рен; Дж. Лян; С. Цюй; Дж. Ли; Б. Лу; Дж. Тиан (2010). «Моделирование Монте-Карло на основе графического процессора для распространения света в сложных гетерогенных тканях». Опция Выражать . 18 (7): 6811–6823. Бибкод : 2010OExpr..18.6811R. дои : 10.1364/oe.18.006811 . ПМИД  20389700.
4. Александр Доронин; Игорь Меглинский (2011). «Онлайн-объектно-ориентированный вычислительный инструмент Монте-Карло для нужд биомедицинской оптики». Биомед. Опция Выражать . 2 (9): 2461–2469. дои : 10.1364/бнэ.2.002461. ПМК 3184856 . ПМИД  21991540.