CC BY
22
8. Терехов, В.И. Тепломассообмен при испарении капель бинарных растворов [Текст] / В.И. Терехов, Н.Е. Шишкин // Тр. 5-й Рос. нац. конф. по тепломассообмену РНКТ-5,— М.: Изд-во МЭИ, 2010,- Т. 4,- С. 302-305.
I Text j / C.T. Crowe // Boca Raton, USA: CRC Press, 2006,- 1156 p.
10. Sazhin, S.S. Particle grouping in oscillating flows [Text| / S.S. Sazhin, T. Shakked, V. Sobolev, D. Katoshevski // European Journal of Mechanics
9. Crowe, C.T. (Ed.) Multiphase flow handbook B/Fluids.- 2008,- Vol. 27,- No 2,- P. 131-149.
УДК 519.245, 535.34, 535.36
A.M. Кривцун, А.Ю. Сетейкин
АНАЛИЗ ПРОЦЕССОВ РАСПРОСТРАНЕНИЯ ОПТИЧЕСКОГО ИЗЛУЧЕНИЯ В БИОЛОГИЧЕСКИХ СРЕДАХ С ИСПОЛЬЗОВАНИЕМ ВЫЧИСЛЕНИЙ НА ГРАФИЧЕСКИХ ПРОЦЕССОРАХ
Моделирование процесса распространения оптического излучения в биологических средах — актуальная задача биомедицинской оптики. Транспорт фотонов в мутных средах, таких как биологические ткани, достаточно хорошо описывается интегрально-дифференциальным уравнением теории переноса излучения [ 1]. Однако аналитическое его решение без введения ряда приближений и упрощений не является тривиальным [2]. Поэтому в оптике биотканей широкое распространение получили численные методы моделирования распространения света [3].
Одним из таких методов является метод Монте-Карло, который широко используется для численного моделирования сложных систем во многих областях науки. К основным его преимуществам над другими методами относятся высокая точность, сравнительная простота реализации, а также возможность работы со средами произвольной пространственной конфигурации. Тем не менее, поскольку метод Монте-Карло стохастический и основывается на многократном повторении испытаний для большого числа частиц, то наиболее явный его недостаток состоит в низкой вычислительной эффективности.
Текущие достижения в области разработки систем на базе многоядерных процессоров сделали возможным преодоление основного недостатка метода [4], позволив тем самым превратить его в удобный инструмент биомедицинской диагностики [5]. В частности, представляется
целесообразным использование методики массивно параллельных вычислений на графических процессорах (ГП). Поскольку расчеты транспорта фотонов в мутных средах осуществляются независимо (частицы в процессе своего распространения не взаимодействуют друг с другом), данный метод хорошо подходит для реализации с использованием массивно параллельных вычислений на ГП.
Метод Монте-Карло
Этот метод широко применяется для численного решения ряда физических и математических задач, так как предлагает гибкий и в то же время достаточно точный подход к моделированию транспорта фотонов в биологических средах [6]. Суть метода состоит в установлении набора общих правил, описывающих поведение фотонов в тканях и определении математических выражений для функций распределения вероятности двух основных характеристик: длины свободного пробега фотонов между двумя последовательными актами взаимодействия, а также направления, в котором фотоны будут продолжать свое движение после акта рассеяния.
Особенно важный этап реализации метода, определяющий его вычислительную эффективность и точность, — это выбор представления моделируемой среды. Наиболее часто биоткань рассматривают в виде многослойной структуры, при этом оптические параметры в пределах каж-
дого слоя считают постоянными. Несмотря на то, что такая модель среды проста и легко реализуема, в действительности биоткани представляют собой геометрически сложные неоднородные структуры. Таким образом, оптически неоднородную среду удобно аппроксимировать разбиением на оптически однородные элементы простой формы (воксели), для каждого из которых определен свой набор оптических параметров: |дй — коэффициент поглощения, — коэффициент рассеяния, g — параметр анизотропии (характеризует угловое распределение рассеяния), п — показатель преломления. К преимуществам такого разбиения следует отнести возможность относительно простого построения точных моделей биологических объектов на основе данных, полученных в результате обработки компьютерных томограмм.
Процесс распространения пакета фотонов в среде состоит в следующем. На этапе инициализации задаются параметры источника излучения: координаты точки входа фотонов в среду и исходное направление их распространения (источникявляется точечным). Длина свободного пробега фотона до следующего акта взаимодействия выбирается в соответствии с выражением [6, 7]:
5 =
-1п'
(1)
И/ = ^0ехр(-^а/),
(2)
где ^ — коэффициент поглощения /-го векселя; Щ) — статистический вес фотона до взаимодействия со средой.
В процессе распространения фотон может пересечь границу раздела двух вокселей с различным набором оптических параметров. В случае, если параметром, претерпевающим скачок, является показатель преломления, то с помощью за-
кона Френеля рассчитывается вероятность фотона пройти через границу раздела двух сред [8]:
2
зт2(а[ -а2) + ^(а1 -а2) вш^а, + а2) tg2(а^ + а2)
,2/
(3)
В случае преломления траектория дальнейшего распространения фотона рассчитывается с помощью закона Снелиуса:
Если изменяется коэффициент рассеяния, то оставшаяся часть длины свободного пробега
нормируется на величину / ^ . Аналогично в случае изменения коэффициента поглощения статистический вес фотона уменьшается в соответствии с коэффициентом поглощения элемента, в котором происходит поглощение. После прохождения фотоном длины свободного пробега следует акт рассеяния: функция распределения вероятности для косинуса угла отклонения (совб) определяется фазовой функцией Хени — Гринштейна [6, 8]:
/?(СО80) =
1
4ТС(1 + £2-2£С089)3/2
(5)
где — случайное число, равномерно распределенное на интервале (0,1]; (д, = |дй + — коэффициент экстинкции.
Фотон перемещается на величину 5, при этом его статистический вес первоначально равный единице, уменьшается пропорционально произведению коэффициента поглощения элемента, в котором происходит взаимодействие, на пройденную длину свободного пробега [7]:
За актом рассеяния следует расчет нового значения шага фотона и описанный процесс повторяется, пока статистический вес частицы не опустится ниже некоторого порогового значения, по достижении которого дальнейшее отслеживание его распространения теряет смысл, поскольку уже не оказывает существенного влияния на конечный результат. Таким образом, после осуществления описанной последовательности действий для каждого из ТУ фотонов, получаем трехмерный массив А(пх, п\\ пг), представляющий собой пространственное распределение поглощенного средой статистического веса фотонов.
Очевидно,что в процессе своего распространения пакеты фотонов не взаимодействуют друг с другом. Это делает возможным одновременное отслеживание траекторий движения множества частиц. Таким образом, задача расчета процесса распространения оптического излучения методом Монте-Карло хорошо подходит для реализации в рамках параллельной архитектуры. Некоторыми авторами были предприняты
попытки такой реализации [4, 5], однако практически в каждой из них предполагается использование дорогостоящего, а следовательно и труднодоступного оборудования. С учетом вышесказанного крайне перспективным представляется применение метода G PG PU (General Purpose Graphics Processing Units — графические процессоры общего назначения), состоящего в использовании графических процессоров (ГП) для неграфических вычислений. Одной из реализаций такого метода является технология CU DA (Compute Unified Device Architecture) [9]. Архитектура позволяет разработчику по своему усмотрению организовывать доступ к набору инструкций ГП и управлять его памятью, таким образом, дает возможность организовывать сложные параллельные вычисления. Скрывая часть сложных технических аспектов программирования ГП, архитектура позволяет значительно упростить, а следовательно и ускорить этап программной реализации математической модели. В рамках методологии программирования CUDA исходный код программы пишется в форме ядер (kernel), которые в большинстве своем схожи с обычными функциями языка программирования СИ. Основное отличие от традиционного исполнения функций центральным процессором (ЦП) состоит в том, что ядра исполняются ГП параллельно. Таким образом, вся задача сначала должна быть разбита на параллельно выполняемые подзадачи, которые затем
будут организованы в блоки и передаваться на исполнение потокам ГП.
Для оценки точности реализации представленного алгоритма в качестве одной из тестовых задач была выбрана задача расчета процесса распространения излучения с длиной волны 633 нм, входящей в диапазон так называемого «терапевтического окна», в области с однородным распределением оптических параметров (|дй = 0,005 мм~', (js= 1,0 мм"1 и g— 0,8), представляющей собой куб со сторонами 40 мм. На рис. 1, а представлено схематичное изображение геометрии моделируемой области. Запуск 106 пакетов фотонов осуществлялся из точки Sc ко-ординатамих= 20 мм, у = 20 mm,z = 0 мм, расположенной на поверхности куба, в направлении оси z перпендикулярно к поверхности расчетной области. Расчетные траектории движения для 1000 пакетов фотонов представлены на рис. 1, б. В среде с таким набором оптических параметров число отдельных актов рассеяния варьируется в диапазоне от 4 до 215 взаимодействий для каждого из пакетов. На рис. 2, а представлено распределение статистического веса фотонов, поглощенных средой (сечение в плоскости xz). Видно, что в этом случае имеет место ярко выраженное рассеяние вперед. На рис. 2, б'представлено распределение поглощенного статистического веса (сечения в плоскостях xz и ху) в среде с таким же, за исключением параметра анизотропии (g— 0,01), набором оптических параметров; имеет место
X, мм
Рис. 1. Схема тестовой расчетной задачи процесса распространения излучения: а — модель однородной оптической среды; 5— источник (стартовая точка) излучения (х = у = 20 мм, I = 0); б— полученные траектории для 1000 пакетов фотонов;
кружки на гранях куба — точки, в которых фотоны покидают расчетную область
Рис. 2. Распределение статистического веса фотонов, поглощенных средой (сечение в плоскости XI, у = 20 мм) для двух различных значений параметра анизотропии g: 0,80 (а) и 0,01 (б); = 0,05 мм4, ц = 1,0 мм"1
практически изотропное рассеяние излучения. Результаты расчетов хорошо согласуются с работами других авторов [6 — 8, 10].
В биомедицинской оптике особый интерес представляет задача расчета распространения излучения в многослойных биологических материалах. В качестве такой среды была выбрана шестислойная модель кожи человека (рис. 3); оптические параметры для каждого из моделируемых слоев приведены в таблице [10]. Расчеты проводились с помощью двух реализаций изложенного метода: однопоточной, для ЦП (систе-
ма с Core i5, частота каждого ядра — 2,27 ГГц), и версии с использованием технологии вычислений на ГП (GeForce GT330, частота ядра — 1,2 ГГц, 48 потоковых процессоров). Версия для ЦП показала производительность, сравнимую с доступными в литературе реализациями [5, 6]; время расчетов составляло 1320 с. В свою очередь версия для ГП продемонстрировала почти 40-кратный, по сравнению с версией для ЦП, прирост вычислительной производительности (время расчетов — 43 с). На рис. 4 представлено распределение поглощенного средой статисти-
Рис. 3. Схема упрощенной модели кожи человека (6 слоев): /— эпидермис; 2— капиллярная дерма; 3, 5— дерма с поверхностным и глубинным сплетением сосудов, соответственно; 4— ретикулярная дерма; 6— гиподерма; 5 — источник излучения (х = у = 5 мм, г = 0)
Рис. 4. Распределение статистического веса фотонов, поглощенных многослойной средой (сечение в плоскости хг)
Параметры модели кожи человека (длина волны X = 633 нм)
Слой Z, мкм я М„) мм ' И,) мм ' g
Эпидермис 96 1,34 0,088 90 0,85
Папиллярная дерма 192 1,39 0,015 19 0,90
Дерма с поверхностным сплетением сосудов 192 1,40 0,013 35 0,95
Ретикулярная дерма 842 1,39 0,010 19 0,76
Дерма с глубинным сплетением сосудов 576 1,40 0,011 25 0,95
Гиподерма 16500 1,44 0,011 13 0,80
Обозначения: z — толщина слоя; п — относительный показатель преломления; |дп — коэффициент поглощения; — коэффициент рассеяния; g — параметр анизотропии.
ческого веса (сечение в плоскости х£). Границы раздела между слоями обозначены пунктиром.
На длине волны источника излучения 633 нм рассеяние превалирует над поглощением; глубина проникновения падающего излучения составляет 2—3 мм. При этом большая часть энергии падающего излучения поглощается приповерхностными слоями: папиллярной дермой, дермой с поверхностным сплетением сосудов и ретикулярной дермой, не проникая на глубину более 1,3 мм. Лишь незначительная часть энергии излучения локализуется гиподермой и дермой с глубинным сплетением сосудов.
Итак, в работе была рассмотрена задача расчета процесса распространения оптического излучения в биологических средах. Моделирование проводилось с использованием метода Монте-Карло, при этом модельная среда аппроксимировалась воксельной кубической сеткой.
Такое разбиение является удобным, поскольку позволяет проводить расчеты для сред с произвольной пространственной конфигурацией, а также сравнительно просто получать модели реальных биологических объектов, например с помощью данных, полученных в результате обработки компьютерных томограмм. Для реализации предложенного подхода была выбрана программная архитектура CUDA, что позволило существенно повысить вычислительную эффективность метода и тем самым преодолеть основной его недостаток без необходимости использования дорогостоящих аппаратных ресурсов. Так, применение в качестве расчетного устройства доступного ГП среднего уровня GeForce GT330, позволило достигнуть 40-кратного прироста производительности вычислений. Полученные результаты моделирования для случаев однородной среды хорошо согласуются с работами других авторов в этой области.
СПИСОК ЛИТЕРАТУРЫ
1. Ishiniaru, A. Wave propagation and scattering in random media [Text] / A. Ishiniaru.— New York: Academic Press, 1978,— 600 c.
2. Hielscher, A.H. Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues [Text] / A.H. Hielscher, R.E. Alcouffe, R.L. Barbour // Physics in Medicine and Biology.— 1998,— Vol. 43,— №'5,- P. 1285-1302.
3. Welch, A.J. Optical-thermal response of laser irradiated tissue |Text| / A.J. Welch, M.J.C. Van Geniert;— New York: Plenum Publishing Corporation, 1995,- 952 p.
4. Page, A.J. Distributed Monte Carlo simulation of light transportation in tissue |Text| / A.J. Page, S. Coyle, T.M. Keane |et al.| // Proc. 20ltl IEEE International Parallel & Distributed Processing Symposi-um.—2006,— P. 254-257.
5. Romero, L.F. Real-time simulation for laser-tissue interaction model [Text] / L.F. Romero, O. Trelles, M.A. Trelles // Parallel Computing: Current & Future Issues of High-End Computing. Proc. Intern. Conf. ParCo- 2005. NIC Series.- № 33,- P. 415-422.
6. Wang, L.H. MCML — Monte Carlo modeling of photon transport in multi-layered tissues [Text] / L.H. Wang, S.L. Jacques, L.Q. Zheng // Computer Methods and Programs in Biomedicine. —1995.— Vol. 47,- P. 131-146.
7. Boas D. Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head [Text] / D. Boas, J. Culver, J. Stott, A. Dunn // Optics Express.— 2002.— Vol. 10,- P. 159-170.
8. Jacques, S.L. Monte Carlo modeling of light transport in tissues [Text| / S.L. Jacques, L.-H. Wang / / Optical-Thermal Response of Laser Irradiated Tissue.— New York: Plenum Publishing Corp. 1995.— P. 73-100.
9. NVIDIA СUDA Programming guide Version 2.3.1 [Электронный pecypc|.— 2009. http:// developer.download.nvidia.com/compute/cuda/ 2_3 j>rod /toolkit/docs/CUDA_C_Program-ming_ Guide.pdf
10. Сетейкин, А.Ю. Моделирование температурных полей с учетом распространения света в биоткани [Текст] / А.Ю. Сетейкин, И.В. Красников, Н.И. Фогель // Известия вузов. Приборостроение,- 2007,- Т. 50,- №9.-'С. 24-28.
