Моделирование распространения излучения в неоднородных средах с использованием вычислений на графических процессорах
УДК 519.245 535.34 535.36
МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ИЗЛУЧЕНИЯ В НЕОДНОРОДНЫХ СРЕДАХ С ИСПОЛЬЗОВАНИЕМ ВЫЧИСЛЕНИЙ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ
© 2012 г.
А. М. Кривцун*; А. Ю. Сетейкин*,**, канд. физ.-мат. наук ** Амурский государственный университет, г. Благовещенск ** Electrophysics Department, Kwangwoon University, Seoul, Korea ** Е-mail: seteikin@mail.ru
В работе рассмотрен ряд вопросов, возникающих при расчете распространения оптического излучения в биологических средах. Получены результаты моделирования для сред с различным набором оптических параметров, включая шестислойную модель кожи человека. Представлен подход для повышения вычислительной эф фективности используемого метода.
Ключевые слова: метод Монте-Карло, воксел, анизотропия, рассеивание, статистический вес.
Коды OCIS: 260.0260, 170.3660.
Поступила в редакцию 11.04.2012.
Введение
Моделирование процесса распространения оптического излучения в биологических средах является актуальной задачей биомедицинской оптики. Транспорт фотонов в таких мутных средах, как биологические ткани, достаточно хорошо описывается интегрально-дифференциальным уравнением теории переноса излучения [1]. Однако его аналитическое решение, без введения ряда приближений и упрощений, не является тривиальным. Поэтому в биомедицинской оптике1 широкое распространение получили численные методы моделирования распространения света [2, 3].
Одним из таких методов является метод Монте-Карло, который широко используется для численного моделирования сложных систем во многих областях науки. К его основным преимуществам относятся высокая точность, сравнительная простота реализации, а также возможность работы со средами произвольной пространственной конфигурации. Однако поскольку метод Монте-Карло является стохастическим и основывается на многократном повторении испытаний для большого
1 Биомедицинская оптика – область науки и техники, основанная на взаимодействии света с биологическими тканями и системами.
числа частиц, то наиболее существенный его недостаток состоит в низкой вычислительной эффективности.2
Текущие достижения в области разработки систем на базе многоядерных процессоров сделали возможным преодоление основного недостатка метода [4], позволив тем самым превратить его в удобный инструмент биомедицинской диагностики [5]. В частности, представляется целесообразным использование методики массивно параллельных вычислений на графических процессорах (ГП). Поскольку расчеты транспорта фотонов в мутных средах являются независимыми (частицы в процессе своего распространения не взаимодействуют друг с другом), данный метод хорошо подходит для реализации с использованием массивно п араллельных вычислений на ГП.
Метод Монте-Карло
Метод Монте-Карло широко применяется для численного решения ряда физических и математических задач, предлагая гибкий и в то же время достаточно точный подход к моделированию транспорта фотонов в био-
2 Вычислительная эффективность оценивается числом итераций, необходимым для получения решения, и затратами машинного времени.
8 “Оптический журнал”, 79, 9, 2012
логических средах [6]. Суть метода состоит в установлении набора общих правил, описывающих поведение фотонов в тканях, и определении математических выражений для функций распределения вероятности двух основных характеристик: длины свободного пробега фотонов между двумя последовательными актами взаимодействия и направления, в котором фотоны будут продолжать свое движение после акта рассеивания.
Особенно важным этапом реализации метода, определяющим его вычислительную эффективность и точность, является выбор представления моделируемой среды. Оптически неоднородную среду удобно аппроксимировать разбиением на оптически однородные элементы простой формы (вокселы3), для каждого из которых определен свой набор оптических параметров: ma – коэффициент поглощения, ms – коэффициент рассеивания, g – параметр анизотропии (характеризует угловое распределение рассеяния) и n – показатель преломления. К преимуществам такого разбиения следует отнести возможность относительно простого построения точных моделей биологических объектов на основе данных, полученных в результате обработки компьютерных томограмм.
Процесс распространения пакета фотонов в среде состоит в следующем. На этапе инициализации задаются параметры источника излучения: координаты точки входа фотонов в среду и исходное направление их распространения (источник является точечным). Длина свободного пробега фотона до следующего акта взаимодействия выбирается в соответствии с выражением [6, 7]
s = -lnx, mt
(1)
где x - случайное число, равномерно распределенное на интервале (0,1];
mt = ma + ms – коэффициент экстинкции. Фотон перемещается на величину s, при этом его статистический вес W, первоначаль-
но равный единице, уменьшается пропорцио-
нально произведению коэффициента поглощения элемента, в котором происходит взаимо действие, на пройденную длину свободного
пробега [7]
3 Воксел (англ. voxel) – это элемент объемного изо бражения, содержащий значение элемента растра в трехмерном пространстве (образовано от слов volumetric и pixel), т. е. это аналог пиксела в трехмерном пространстве.
W = W0 exp(-mai l),
(2)
где mai – коэффициент поглощения i-го воксела, W0 – статистический вес фотона до взаимодействия со средой.
В процессе распространения фотон может пересечь границу раздела двух вокселов с различным набором оптических параметров. Если параметром, претерпевающим скачок, является показатель преломления, то с помощью закона Френеля рассчитывается вероятность фотона пройти через границу раздела двух сред. Если изменяется коэффициент рассеивания, то оставшаяся часть длины свободного пробега нормируется на величину msi–1/mis. Аналогично, в случае изменения коэффициента поглощения, статистический вес фотона уменьшается в соответствии с коэффициентом поглощения элемента, в котором происходит поглощение. После прохождения фотоном длины свободного пробега следует акт рассеивания: функция распределения вероятности для косинуса угла отклонения (cosθ) определяется фазовой функцией Хени-Гринштейна [6, 8] как
p(cos θ)
=
1 4p
(1 +
1- g2 g2 - 2g cos θ)3/2
.
(3)
За актом рассеивания следует расчет нового значения шага фотона и описанный процесс повторяется, пока статистический вес частицы не опустится ниже некоторого порогового значения, по достижении которого дальнейшее отслеживание его распространения теряет смысл, поскольку уже не оказывает существенного влияния на конечный результат. Таким образом, после осуществления описанной последовательности действий для каждого из N фотонов получаем трехмерный массив A(nx, ny, nz), представляющий собой пространственное распределение поглощенного средой статистического веса фотонов.
Очевидно, что в процессе распространения пакеты фотонов не взаимодействуют друг с другом. Это делает возможным одновременное отслеживание траекторий движения множества частиц. Таким образом, задача расчета процесса распространения оптического излучения методом Монте-Карло хорошо подходит для реализации в рамках параллельной архитектуры. Некоторыми авторами были предприняты попытки такой реализации [4, 5], однако практически в каждой из них предполагается использование дорогостоящего, а следовательно, и труднодоступного оборудования.
“Оптический журнал”, 79, 9, 2012
9
С учетом сказанного крайне перспективным представляется применение метода GPGPU (General Purpose Graphics Processing Units – графические процессоры общего назначения), состоящего в использовании ГП для неграфических вычислений.
Одной из реализаций такого метода является технология CUDA (Compute Unified Device Architecture) [9], позволяющая разработчику по своему усмотрению организовывать доступ к набору инструкций ГП и управлять его памятью. Таким образом, она дает возможность организовывать сложные параллельные вычисления. Скрывая часть сложных технических аспектов программирования ГП, технология CUDA позволяет значительно упростить, а следовательно, и ускорить этап программной реализации математической модели.
В рамках методологии программирования CUDA исходный код программы пишется в форме ядер (kernel), которые в большинстве случаев схожи с обычными функциями языка программирования Си. Основное отличие от традиционного исполнения функций центральным процессором (ЦП) состоит в том, что ядра исполняются ГП параллельно. Таким образом, вся задача сначала должна быть разбита на п араллельно выполняемые подзадачи, которые затем организуются в блоки и передаются на исполнение потокам ГП. Следует отметить, что некоторые математические операции являются в большей степени вычислительно затратными: например, операция деления выполняется за 38 тактов, в то время как умножение –
всего лишь за 4. Однако в большинстве случаев операцию деления удается заменить умножением. Например, при многократно повторяющемся расчете шага фотона до следующего акта взаимодействия с помощью выражения (1), представляется целесообразным вместо хранения в памяти значений коэффициента экстинкции µt для каждого элемента среды хранить предварительно рассчитанные значения, обратные этой величине, – 1/mt. Такой подход может быть применен ко многим вычислениям и позволяет повысить производительность р еализации метода.
Анализ результатов
Для оценки точности реализации представленного алгоритма в качестве одной из тестовых задач была выбрана задача расчета процесса распространения излучения с длиной волны 633 нм, входящей в диапазон так называемого “терапевтического окна”, в области с однородным распределением оптических параметров (ma = 0,005 мм–1, ms = 1,0 мм–1 и g = 0,8), представляющей собой куб со сторонами 40 мм. На рис. 1а представлена схема модели однородной среды. Запуск 106 пакетов фотонов осуществлялся в точке с координатами x = 20 мм, y = 20 мм, z = 0 мм, расположенной на поверхности области в направлении оси z перпендикулярно к поверхности расчетной области. На рис. 1б представлено распределение поглощенного средой статистического веса фотонов (срез в плоскости xz).
z, мм z, мм
(а)
Источник излучения (х = 20 мм, y = 20 мм, z = 0 мм)
10
20
30
40
30
x,
мм20
10
10
20 30
y, мм
40
(б)
12 5
11 10
10 15 9
20 8
25 7
30 6 5
35 4
40 5 10 15 20 25 30 35 40
x, мм
Рис. 1. Результаты моделирования. а – схема модели однородной среды, б – распределение поглощенного статистического веса (срез по xz).
10 “Оптический журнал”, 79, 9, 2012
(а)
5 10 15
12 (б)
5 10
10
8 15
600 500 400
z, мм y, мм
20 6 20 300
25 4
30 2 35
0 40
5 10 15 20 25 30 35 40
x, мм
25 200
30
35 100
40 5 10 15 20 25 30 35 40
x, мм
Рис. 2. Распределение поглощенного средой статистического веса фотонов. а – срез в плоскости xz (y = 20 мм), б – срез в плоскости xy (z = 20 мм).
z, мм
Источник излучения (х = 5 мм, y = 5 мм, z = 0 мм)
2 4 6
1 23 4
5
6
8
10
10
8
6
x, мм
4
2
2
4
68
y, мм
10
Рис. 3. Схема упрощенной 6-слойной модели кожи человека. 1 – эпидермис, 2 – папиллярная дерма, 3 – дерма с поверхностным сплетением сосудов, 4 – ретикулярная дерма, 5 – дерма с глубинным сплетением сосудов, 6 – гиподерма.
Видно, что в этом случае имеет место ярко выраженное рассеяние вперед. На рис. 2 представлены распределения поглощенного статистического веса в среде с таким же, за исключением анизотропии (g = 0,01), набором оптических параметров: срезы в плоскостях xz (а) и xy (б). В этом случае имеет место близкое к изотропному рассеяние излучения.
В биомедицинской оптике особый интерес представляет задача расчета распространения излучения в многослойных биологических материалах. В качестве такой среды была выбрана шестислойная модель кожи человека, которая состоит из эпидермиса, папиллярной дермы, дермы с поверхностным сплетением сосудов, ретикулярной дермы, дермы с глубинным сплетением сосудов и гиподермы (рис. 3). Оптические параметры для каждого из моделируемых слоев приведены в таблице [10, 11]. Расчеты проводились с помощью двух реализа-
Параметры модели кожи человека (длина волны l = 633 нм)
Слой
z, мкм
n
ma, мм–1
ms, мм–1
g
Эпидермис
96
1,34
0,088
90
0,85
Папиллярная дерма
192
1,39
0,015
18,7
0,9
Дерма с поверхностным сплетением сосудов
192
1,4 0,013
35
0,95
Ретикулярная дерма
842
1,39
0,01
19
0,76
Дерма с глубинным сплетением сосудов
576
1,4 0,011
25
0,95
Гиподерма
16 500
1,44
0,011
13
0,8
Примечание. z – толщина слоя, n – относительный показатель преломления, ma – коэффициент поглощения, ms – коэффициент рассеивания, g – параметр анизотропии.
“Оптический журнал”, 79, 9, 2012
11
z, мм
1 2 2000
3 4 1500
5
6 1000
7
8 500
9
10 1 2 3 4 5 0
x, мм
Рис. 4. Распределение поглощенного многослойной средой статистического веса фотонов (срез в плоскости xz).
ций изложенного метода – однопоточной, для ЦП (система с Core i5, частота каждого ядра 2,27 ГГц) и версии с использованием технологии вычислений на ГП (GeForce GT330, частота ядра 1,2 ГГц, 48 потоковых процессоров). Версия для ЦП показала сравнимую с доступными в литературе реализациями [5–6] производительность (время расчетов – 1320 с). В свою очередь версия для ГП продемонстрировала почти 40-кратный, по сравнению с версией для ЦП, прирост вычислительной производительности (время расчетов – 43 с). На рис. 4 представлено распределение поглощенного средой статистического веса фотонов (срез в плоскости xz). Границы раздела между слоями обозначены пунктиром. Видно, что значительная часть энергии падающего излучения поглощается
слоями папиллярной дермы, дермы с поверхностным сплетением сосудов и ретикулярной дермой, не проникая на глубину более 1,3 мм. Лишь незначительная часть энергии излучения локализуется гиподермой и дермой с глубинным сплетением сосудов.
Заключение
В работе рассмотрена задача расчета процесса распространения оптического излучения в биологических средах. Моделирование проводилось с использованием метода Монте-Карло, при этом модельная среда аппроксимировалась воксельной кубической сеткой. Такое разбиение является удобным, поскольку позволяет проводить расчеты для сред с произвольной пространственной конфигурацией, а также сравнительно просто получать модели реальных биологических объектов, например, с помощью данных, полученных в результате обработки компьютерных томограмм. Для реализации предложенного подхода была выбрана программная архитектура CUDA, что позволило существенно повысить вычислительную эффективность метода и тем самым преодолеть основной его недостаток без необходимости использования дорогостоящих аппаратных ресурсов. Так, применение в качестве расчетного устройства доступного графического процессора среднего уровня GeForce GT330 позволило достигнуть 40-кратного прироста производительности вычислений. Полученные результаты моделирования для случаев однородной среды хорошо согласуются с данными других работ в этой области [6–10, 12].
Работа выполнена при поддержке исследовательского гранта Kwangwoon University (Корея).
* * * * *
ЛИТЕРАТУРА
1. Ishimaru A. Wave propagation and scattering in random media. N.Y.: Academic Press, 1978. 600 p.
2. Hielscher A.H., Alcouffe R.E., Barbour R.L. Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues // Physics in Medicine and Biology. 1998. V. 43. № 5. P. 1285–1302.
3. Welch A.J., van Gemert M.J.C. Optical-thermal response of laser irradiated tissue. N.Y.: Plenum Publishing Corporation, 1995. 952 p.
12 “Оптический журнал”, 79, 9, 2012
4. Page A.J., Coyle S., Keane T.M., Naughton T.J., Markham C., Ward T. Distributed Monte Carlo simulation of light transportation in tissue // Proc. 20th IEEE International Parallel & Distributed Processing Symposium. 2006. 254 p.
5. Romero L.F., Trelles O., Trelles M.A. Real-time simulation for Laser-Tissue Interaction Model // Parallel Computing: Current & Future Issues of High-End Computing, Proceedings of the International Conference ParCo – 2005. NIC Series. №. 33. P. 415–422.
6. Wang L.H., Jacques S.L. MCML – Monte Carlo modeling of photon transport in multi-layered tissues // Computer Methods and Programs in Biomedicine. 1995. V. 47. P. 131–146.
7. Boas D., Culver J., Stott J., Dunn A. Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head // Optics Express. 2002. № 10. P. 159–170.
8. Jacques S.L., Wang L.-H. Monte Carlo modeling of light transport in tissues // Optical-Thermal Response of Laser Irradiated Tissue. N.Y.: Plenum Publishing Corp. 1995. P. 73–100.
9. Lo W.C.Y., Han T.D., Rose J., Lilge L. GPU-accelerated Monte Carlo simulation for photodynamic therapy treatment planning // Proc. SPIE. 2009. V. 7373. P. 737313.
10. Меглинский И.В. Моделирование спектров отражения оптического излучения от случайно-неоднородных многослойных сильно рассеивающих и поглощающих свет сред методом Монте-Карло // Квант. электрон. 2001. T. 31. № 12. С. 1101–1107.
11. Тучин В.В. Лазеры и волоконная оптика в биомедицинских исследованиях. Саратов: Изд-во Сарат. ун-та, 1998. 348 с.
12. Сетейкин А.Ю., Красников И.В., Фогель Н.И. Моделирование температурных полей с учетом распространения света в биоткани // Изв. вузов. Приборостроение. 2007. Т. 50. № 9. С. 24–28.
“Оптический журнал”, 79, 9, 2012
13
МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ИЗЛУЧЕНИЯ В НЕОДНОРОДНЫХ СРЕДАХ С ИСПОЛЬЗОВАНИЕМ ВЫЧИСЛЕНИЙ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ
© 2012 г.
А. М. Кривцун*; А. Ю. Сетейкин*,**, канд. физ.-мат. наук ** Амурский государственный университет, г. Благовещенск ** Electrophysics Department, Kwangwoon University, Seoul, Korea ** Е-mail: seteikin@mail.ru
В работе рассмотрен ряд вопросов, возникающих при расчете распространения оптического излучения в биологических средах. Получены результаты моделирования для сред с различным набором оптических параметров, включая шестислойную модель кожи человека. Представлен подход для повышения вычислительной эф фективности используемого метода.
Ключевые слова: метод Монте-Карло, воксел, анизотропия, рассеивание, статистический вес.
Коды OCIS: 260.0260, 170.3660.
Поступила в редакцию 11.04.2012.
Введение
Моделирование процесса распространения оптического излучения в биологических средах является актуальной задачей биомедицинской оптики. Транспорт фотонов в таких мутных средах, как биологические ткани, достаточно хорошо описывается интегрально-дифференциальным уравнением теории переноса излучения [1]. Однако его аналитическое решение, без введения ряда приближений и упрощений, не является тривиальным. Поэтому в биомедицинской оптике1 широкое распространение получили численные методы моделирования распространения света [2, 3].
Одним из таких методов является метод Монте-Карло, который широко используется для численного моделирования сложных систем во многих областях науки. К его основным преимуществам относятся высокая точность, сравнительная простота реализации, а также возможность работы со средами произвольной пространственной конфигурации. Однако поскольку метод Монте-Карло является стохастическим и основывается на многократном повторении испытаний для большого
1 Биомедицинская оптика – область науки и техники, основанная на взаимодействии света с биологическими тканями и системами.
числа частиц, то наиболее существенный его недостаток состоит в низкой вычислительной эффективности.2
Текущие достижения в области разработки систем на базе многоядерных процессоров сделали возможным преодоление основного недостатка метода [4], позволив тем самым превратить его в удобный инструмент биомедицинской диагностики [5]. В частности, представляется целесообразным использование методики массивно параллельных вычислений на графических процессорах (ГП). Поскольку расчеты транспорта фотонов в мутных средах являются независимыми (частицы в процессе своего распространения не взаимодействуют друг с другом), данный метод хорошо подходит для реализации с использованием массивно п араллельных вычислений на ГП.
Метод Монте-Карло
Метод Монте-Карло широко применяется для численного решения ряда физических и математических задач, предлагая гибкий и в то же время достаточно точный подход к моделированию транспорта фотонов в био-
2 Вычислительная эффективность оценивается числом итераций, необходимым для получения решения, и затратами машинного времени.
8 “Оптический журнал”, 79, 9, 2012
логических средах [6]. Суть метода состоит в установлении набора общих правил, описывающих поведение фотонов в тканях, и определении математических выражений для функций распределения вероятности двух основных характеристик: длины свободного пробега фотонов между двумя последовательными актами взаимодействия и направления, в котором фотоны будут продолжать свое движение после акта рассеивания.
Особенно важным этапом реализации метода, определяющим его вычислительную эффективность и точность, является выбор представления моделируемой среды. Оптически неоднородную среду удобно аппроксимировать разбиением на оптически однородные элементы простой формы (вокселы3), для каждого из которых определен свой набор оптических параметров: ma – коэффициент поглощения, ms – коэффициент рассеивания, g – параметр анизотропии (характеризует угловое распределение рассеяния) и n – показатель преломления. К преимуществам такого разбиения следует отнести возможность относительно простого построения точных моделей биологических объектов на основе данных, полученных в результате обработки компьютерных томограмм.
Процесс распространения пакета фотонов в среде состоит в следующем. На этапе инициализации задаются параметры источника излучения: координаты точки входа фотонов в среду и исходное направление их распространения (источник является точечным). Длина свободного пробега фотона до следующего акта взаимодействия выбирается в соответствии с выражением [6, 7]
s = -lnx, mt
(1)
где x - случайное число, равномерно распределенное на интервале (0,1];
mt = ma + ms – коэффициент экстинкции. Фотон перемещается на величину s, при этом его статистический вес W, первоначаль-
но равный единице, уменьшается пропорцио-
нально произведению коэффициента поглощения элемента, в котором происходит взаимо действие, на пройденную длину свободного
пробега [7]
3 Воксел (англ. voxel) – это элемент объемного изо бражения, содержащий значение элемента растра в трехмерном пространстве (образовано от слов volumetric и pixel), т. е. это аналог пиксела в трехмерном пространстве.
W = W0 exp(-mai l),
(2)
где mai – коэффициент поглощения i-го воксела, W0 – статистический вес фотона до взаимодействия со средой.
В процессе распространения фотон может пересечь границу раздела двух вокселов с различным набором оптических параметров. Если параметром, претерпевающим скачок, является показатель преломления, то с помощью закона Френеля рассчитывается вероятность фотона пройти через границу раздела двух сред. Если изменяется коэффициент рассеивания, то оставшаяся часть длины свободного пробега нормируется на величину msi–1/mis. Аналогично, в случае изменения коэффициента поглощения, статистический вес фотона уменьшается в соответствии с коэффициентом поглощения элемента, в котором происходит поглощение. После прохождения фотоном длины свободного пробега следует акт рассеивания: функция распределения вероятности для косинуса угла отклонения (cosθ) определяется фазовой функцией Хени-Гринштейна [6, 8] как
p(cos θ)
=
1 4p
(1 +
1- g2 g2 - 2g cos θ)3/2
.
(3)
За актом рассеивания следует расчет нового значения шага фотона и описанный процесс повторяется, пока статистический вес частицы не опустится ниже некоторого порогового значения, по достижении которого дальнейшее отслеживание его распространения теряет смысл, поскольку уже не оказывает существенного влияния на конечный результат. Таким образом, после осуществления описанной последовательности действий для каждого из N фотонов получаем трехмерный массив A(nx, ny, nz), представляющий собой пространственное распределение поглощенного средой статистического веса фотонов.
Очевидно, что в процессе распространения пакеты фотонов не взаимодействуют друг с другом. Это делает возможным одновременное отслеживание траекторий движения множества частиц. Таким образом, задача расчета процесса распространения оптического излучения методом Монте-Карло хорошо подходит для реализации в рамках параллельной архитектуры. Некоторыми авторами были предприняты попытки такой реализации [4, 5], однако практически в каждой из них предполагается использование дорогостоящего, а следовательно, и труднодоступного оборудования.
“Оптический журнал”, 79, 9, 2012
9
С учетом сказанного крайне перспективным представляется применение метода GPGPU (General Purpose Graphics Processing Units – графические процессоры общего назначения), состоящего в использовании ГП для неграфических вычислений.
Одной из реализаций такого метода является технология CUDA (Compute Unified Device Architecture) [9], позволяющая разработчику по своему усмотрению организовывать доступ к набору инструкций ГП и управлять его памятью. Таким образом, она дает возможность организовывать сложные параллельные вычисления. Скрывая часть сложных технических аспектов программирования ГП, технология CUDA позволяет значительно упростить, а следовательно, и ускорить этап программной реализации математической модели.
В рамках методологии программирования CUDA исходный код программы пишется в форме ядер (kernel), которые в большинстве случаев схожи с обычными функциями языка программирования Си. Основное отличие от традиционного исполнения функций центральным процессором (ЦП) состоит в том, что ядра исполняются ГП параллельно. Таким образом, вся задача сначала должна быть разбита на п араллельно выполняемые подзадачи, которые затем организуются в блоки и передаются на исполнение потокам ГП. Следует отметить, что некоторые математические операции являются в большей степени вычислительно затратными: например, операция деления выполняется за 38 тактов, в то время как умножение –
всего лишь за 4. Однако в большинстве случаев операцию деления удается заменить умножением. Например, при многократно повторяющемся расчете шага фотона до следующего акта взаимодействия с помощью выражения (1), представляется целесообразным вместо хранения в памяти значений коэффициента экстинкции µt для каждого элемента среды хранить предварительно рассчитанные значения, обратные этой величине, – 1/mt. Такой подход может быть применен ко многим вычислениям и позволяет повысить производительность р еализации метода.
Анализ результатов
Для оценки точности реализации представленного алгоритма в качестве одной из тестовых задач была выбрана задача расчета процесса распространения излучения с длиной волны 633 нм, входящей в диапазон так называемого “терапевтического окна”, в области с однородным распределением оптических параметров (ma = 0,005 мм–1, ms = 1,0 мм–1 и g = 0,8), представляющей собой куб со сторонами 40 мм. На рис. 1а представлена схема модели однородной среды. Запуск 106 пакетов фотонов осуществлялся в точке с координатами x = 20 мм, y = 20 мм, z = 0 мм, расположенной на поверхности области в направлении оси z перпендикулярно к поверхности расчетной области. На рис. 1б представлено распределение поглощенного средой статистического веса фотонов (срез в плоскости xz).
z, мм z, мм
(а)
Источник излучения (х = 20 мм, y = 20 мм, z = 0 мм)
10
20
30
40
30
x,
мм20
10
10
20 30
y, мм
40
(б)
12 5
11 10
10 15 9
20 8
25 7
30 6 5
35 4
40 5 10 15 20 25 30 35 40
x, мм
Рис. 1. Результаты моделирования. а – схема модели однородной среды, б – распределение поглощенного статистического веса (срез по xz).
10 “Оптический журнал”, 79, 9, 2012
(а)
5 10 15
12 (б)
5 10
10
8 15
600 500 400
z, мм y, мм
20 6 20 300
25 4
30 2 35
0 40
5 10 15 20 25 30 35 40
x, мм
25 200
30
35 100
40 5 10 15 20 25 30 35 40
x, мм
Рис. 2. Распределение поглощенного средой статистического веса фотонов. а – срез в плоскости xz (y = 20 мм), б – срез в плоскости xy (z = 20 мм).
z, мм
Источник излучения (х = 5 мм, y = 5 мм, z = 0 мм)
2 4 6
1 23 4
5
6
8
10
10
8
6
x, мм
4
2
2
4
68
y, мм
10
Рис. 3. Схема упрощенной 6-слойной модели кожи человека. 1 – эпидермис, 2 – папиллярная дерма, 3 – дерма с поверхностным сплетением сосудов, 4 – ретикулярная дерма, 5 – дерма с глубинным сплетением сосудов, 6 – гиподерма.
Видно, что в этом случае имеет место ярко выраженное рассеяние вперед. На рис. 2 представлены распределения поглощенного статистического веса в среде с таким же, за исключением анизотропии (g = 0,01), набором оптических параметров: срезы в плоскостях xz (а) и xy (б). В этом случае имеет место близкое к изотропному рассеяние излучения.
В биомедицинской оптике особый интерес представляет задача расчета распространения излучения в многослойных биологических материалах. В качестве такой среды была выбрана шестислойная модель кожи человека, которая состоит из эпидермиса, папиллярной дермы, дермы с поверхностным сплетением сосудов, ретикулярной дермы, дермы с глубинным сплетением сосудов и гиподермы (рис. 3). Оптические параметры для каждого из моделируемых слоев приведены в таблице [10, 11]. Расчеты проводились с помощью двух реализа-
Параметры модели кожи человека (длина волны l = 633 нм)
Слой
z, мкм
n
ma, мм–1
ms, мм–1
g
Эпидермис
96
1,34
0,088
90
0,85
Папиллярная дерма
192
1,39
0,015
18,7
0,9
Дерма с поверхностным сплетением сосудов
192
1,4 0,013
35
0,95
Ретикулярная дерма
842
1,39
0,01
19
0,76
Дерма с глубинным сплетением сосудов
576
1,4 0,011
25
0,95
Гиподерма
16 500
1,44
0,011
13
0,8
Примечание. z – толщина слоя, n – относительный показатель преломления, ma – коэффициент поглощения, ms – коэффициент рассеивания, g – параметр анизотропии.
“Оптический журнал”, 79, 9, 2012
11
z, мм
1 2 2000
3 4 1500
5
6 1000
7
8 500
9
10 1 2 3 4 5 0
x, мм
Рис. 4. Распределение поглощенного многослойной средой статистического веса фотонов (срез в плоскости xz).
ций изложенного метода – однопоточной, для ЦП (система с Core i5, частота каждого ядра 2,27 ГГц) и версии с использованием технологии вычислений на ГП (GeForce GT330, частота ядра 1,2 ГГц, 48 потоковых процессоров). Версия для ЦП показала сравнимую с доступными в литературе реализациями [5–6] производительность (время расчетов – 1320 с). В свою очередь версия для ГП продемонстрировала почти 40-кратный, по сравнению с версией для ЦП, прирост вычислительной производительности (время расчетов – 43 с). На рис. 4 представлено распределение поглощенного средой статистического веса фотонов (срез в плоскости xz). Границы раздела между слоями обозначены пунктиром. Видно, что значительная часть энергии падающего излучения поглощается
слоями папиллярной дермы, дермы с поверхностным сплетением сосудов и ретикулярной дермой, не проникая на глубину более 1,3 мм. Лишь незначительная часть энергии излучения локализуется гиподермой и дермой с глубинным сплетением сосудов.
Заключение
В работе рассмотрена задача расчета процесса распространения оптического излучения в биологических средах. Моделирование проводилось с использованием метода Монте-Карло, при этом модельная среда аппроксимировалась воксельной кубической сеткой. Такое разбиение является удобным, поскольку позволяет проводить расчеты для сред с произвольной пространственной конфигурацией, а также сравнительно просто получать модели реальных биологических объектов, например, с помощью данных, полученных в результате обработки компьютерных томограмм. Для реализации предложенного подхода была выбрана программная архитектура CUDA, что позволило существенно повысить вычислительную эффективность метода и тем самым преодолеть основной его недостаток без необходимости использования дорогостоящих аппаратных ресурсов. Так, применение в качестве расчетного устройства доступного графического процессора среднего уровня GeForce GT330 позволило достигнуть 40-кратного прироста производительности вычислений. Полученные результаты моделирования для случаев однородной среды хорошо согласуются с данными других работ в этой области [6–10, 12].
Работа выполнена при поддержке исследовательского гранта Kwangwoon University (Корея).
* * * * *
ЛИТЕРАТУРА
1. Ishimaru A. Wave propagation and scattering in random media. N.Y.: Academic Press, 1978. 600 p.
2. Hielscher A.H., Alcouffe R.E., Barbour R.L. Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues // Physics in Medicine and Biology. 1998. V. 43. № 5. P. 1285–1302.
3. Welch A.J., van Gemert M.J.C. Optical-thermal response of laser irradiated tissue. N.Y.: Plenum Publishing Corporation, 1995. 952 p.
12 “Оптический журнал”, 79, 9, 2012
4. Page A.J., Coyle S., Keane T.M., Naughton T.J., Markham C., Ward T. Distributed Monte Carlo simulation of light transportation in tissue // Proc. 20th IEEE International Parallel & Distributed Processing Symposium. 2006. 254 p.
5. Romero L.F., Trelles O., Trelles M.A. Real-time simulation for Laser-Tissue Interaction Model // Parallel Computing: Current & Future Issues of High-End Computing, Proceedings of the International Conference ParCo – 2005. NIC Series. №. 33. P. 415–422.
6. Wang L.H., Jacques S.L. MCML – Monte Carlo modeling of photon transport in multi-layered tissues // Computer Methods and Programs in Biomedicine. 1995. V. 47. P. 131–146.
7. Boas D., Culver J., Stott J., Dunn A. Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head // Optics Express. 2002. № 10. P. 159–170.
8. Jacques S.L., Wang L.-H. Monte Carlo modeling of light transport in tissues // Optical-Thermal Response of Laser Irradiated Tissue. N.Y.: Plenum Publishing Corp. 1995. P. 73–100.
9. Lo W.C.Y., Han T.D., Rose J., Lilge L. GPU-accelerated Monte Carlo simulation for photodynamic therapy treatment planning // Proc. SPIE. 2009. V. 7373. P. 737313.
10. Меглинский И.В. Моделирование спектров отражения оптического излучения от случайно-неоднородных многослойных сильно рассеивающих и поглощающих свет сред методом Монте-Карло // Квант. электрон. 2001. T. 31. № 12. С. 1101–1107.
11. Тучин В.В. Лазеры и волоконная оптика в биомедицинских исследованиях. Саратов: Изд-во Сарат. ун-та, 1998. 348 с.
12. Сетейкин А.Ю., Красников И.В., Фогель Н.И. Моделирование температурных полей с учетом распространения света в биоткани // Изв. вузов. Приборостроение. 2007. Т. 50. № 9. С. 24–28.
“Оптический журнал”, 79, 9, 2012
13