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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

${\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 — это расстояние, которое пакет фотонов проходит между местами взаимодействия. Существует множество методов выбора размера шага. Ниже приведена базовая форма выбора размера шага фотона (полученная с использованием метода обратного распределения и закона Бера-Ламберта ), которую мы используем для нашей однородной модели:

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

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

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

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

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

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

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

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

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

$W\leftarrow W-\Delta W\,$

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

$\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 и Обычно предполагается, что $2\pi$. Исходя из этого предположения, можно установить:

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

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

${\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}}$

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

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

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

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

или

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

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

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

C-код:

/*********************** Indicatrix *********************
*New direction cosines after scattering by angle theta, fi.
* 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)
*---------------------------------------------------------
*Input:
* muxs,muys,muzs - direction cosine before collision
* mutheta, fi - cosine of polar angle and the azimuthal angle
*---------------------------------------------------------
*Output:
*  muxd,muyd,muzd - direction cosine after collision
*---------------------------------------------------------
*/
void Indicatrix (double muxs, double muys, double muzs, double mutheta, double fi, double *muxd, double *muyd, double *muzd)
{
 double costheta = mutheta;
 double sintheta = sqrt(1.0-costheta*costheta); // sin(theta)
 double sinfi = sin(fi);
 double cosfi = cos(fi);
 if (muzs == 1.0) {
   *muxd = sintheta*cosfi;
   *muyd = sintheta*sinfi;
   *muzd = costheta;
 } elseif (muzs == -1.0) {
   *muxd = sintheta*cosfi;
   *muyd = -sintheta*sinfi;
   *muzd = -costheta;
 } else {
   double denom = sqrt(1.0-muzs*muzs);
   double muzcosfi = muzs*cosfi;
   *muxd = sintheta*(muxs*muzcosfi-muys*sinfi)/denom + muxs*costheta;
   *muyd = sintheta*(muys*muzcosfi+muxs*sinfi)/denom + muys*costheta;
   *muzd = -denom*sintheta*cosfi + muzs*costheta;
 }
}

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

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

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

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

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

Этот базовый подход сам по себе можно распараллелить, используя несколько связанных вместе графических процессоров. Одним из примеров является «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 .