Метод Монте-Карло для переноса фотонов - Monte Carlo method for photon transport

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

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

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

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

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

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

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

Радиационная терапия

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

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

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

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

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

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

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

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

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

${displaystyle {egin {align} x & = 0 {ext {Position:}} y & = 0 z & = 0 mu _ {x} & = 0 {ext {Направляющие косинусы:}} mu _ {y} & = 0 mu _ {z} & = 1end {выровнено}}}$

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

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

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

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

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

${displaystyle {egin {выравнивается} x & стрелка влево x + mu _ {x} s y & стрелка влево y + mu _ {y} s z & стрелка влево z + mu _ {z} send {выровнено}}}$

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

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

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

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

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

${displaystyle Wleftarrow W-Delta W,}$

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

${displaystyle cos heta = {egin {case} {frac {1} {2g}} left [1 + g ^ {2} -left ({frac {1-g ^ {2}} {1-g + 2gxi}}) ight) ^ {2} ight] & {ext {if}} geq 0 1-2xi & {ext {if}} g = 0end {case}}}$

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

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

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

${displaystyle {egin {align} mu '_ {x} & = {frac {sin heta (mu _ {x} mu _ {z} cos varphi -mu _ {y} sin varphi)} {sqrt {1-mu _ {z} ^ {2}}}} + mu _ {x} cos heta mu '_ {y} & = {frac {sin heta (mu _ {y} mu _ {z} cos varphi + mu _ {x } sin varphi)} {sqrt {1-mu _ {z} ^ {2}}}} + mu _ {y} cos heta mu '_ {z} & = - {sqrt {1-mu _ {z} ^ {2}}} sin heta cos varphi + mu _ {z} cos heta end {выровнено}}}$

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

${displaystyle {egin {выравнивается} mu _ {z} = 1end {выравнивается}}}$

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

${displaystyle {egin {выровнено} mu '_ {x} & = sin heta cos varphi mu' _ {y} & = sin heta sin varphi mu '_ {z} & = cos heta end {выровнено}}}$

или же

${displaystyle {egin {выравнивается} mu _ {z} = - 1end {выравнивается}}}$

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

${displaystyle {egin {выровнено} mu '_ {x} & = sin heta cos varphi mu' _ {y} & = - sin heta sin varphi mu '_ {z} & = - cos heta end {выровнено} }}$

C-код:

/ *********************** Индикаторная матрица ********************* Косинусы нового направления после рассеяния на угол тета, 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) * ----------------------- ---------------------------------- * Ввод: * muxs, muys, muzs - косинус направления перед столкновением * mutheta , fi - косинус полярного угла и азимутального угла * -------------------------------------- ------------------- * Вывод: * muxd, muyd, muzd - косинус направления после столкновения * ---------------- ----------------------------------------- * / void Indicatrix (двойной мультиплексор, двойной muys, double muzs, double mutheta, double fi, double * muxd, double * muyd, double * muzd) {двойная costheta = mutheta; двойная синтета = sqrt (1.0-costheta * costheta); // sin (theta) double sinfi = sin (fi); двойной cosfi = cos (fi); если (muzs == 1.0) {* muxd = sintheta * cosfi; * muyd = sintheta * sinfi; * музд = костета; } elseif (muzs == -1.0) {* muxd = sintheta * cosfi; * muyd = -sintheta * sinfi; * muzd = -costheta; } else {двойной деном = sqrt (1.0-muzs * muzs); двойной muzcosfi = muzs * cosfi; * muxd = sintheta * (muxs * muzcosfi-muys * sinfi) / denom + muxs * costheta; * muyd = sintheta * (muys * muzcosfi + muxs * sinfi) / деном + muys * costheta; * muzd = -denom * sintheta * cosfi + muzs * costheta; }}

Шаг 4: прекращение действия фотона

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

${displaystyle W = {egin {case} mW & xi leq 1 / m 0 & xi> 1 / исправить {case}}}$

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

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

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

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

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

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

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

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

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

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