Использование графических процессоров для решения задачи распространения света в диффузионной флуоресцентной томографии методом Монте-Карло Текст научной статьи по специальности «Математика»

Информационные технологии Вестник Нижегородского университета им. Н.И. Лобачевского, 2011, № 4 (1), с. 190-195

УДК 004.94+535.36

ИСПОЛЬЗОВАНИЕ ГРАФИЧЕСКИХ ПРОЦЕССОРОВ ДЛЯ РЕШЕНИЯ ЗАДАЧИ РАСПРОСТРАНЕНИЯ СВЕТА В ДИФФУЗИОННОЙ ФЛУОРЕСЦЕНТНОЙ ТОМОГРАФИИ МЕТОДОМ МОНТЕ-КАРЛО

© 2011 г. И.И. Фикс

Институт прикладной физики РАН, Н. Новгород

FiksII@yandex.ru

Поступила в редакцию 05.03.2011

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

Ключевые слова: математическое моделирование, распространение света, СИЭА, параллельные вычисления, метод Монте-Карло.

Введение

Численное моделирование различных задач повышенной сложности обычно требует больших объемов оперативной памяти и скорости вычислений. И если для решения большинства задач требуемый объем оперативной памяти может быть сравнительно легко обеспечен добавлением модулей памяти в существующую систему, то быстродействие одного вычислительного блока (ядра процессора) имеет определенный порог. Таким образом, увеличение производительности вычислений возможно только при увеличении количества ядер процессора и соответствующем распараллеливании алгоритма вычислений. Однако архитектура современных компьютеров не позволяет значительно наращивать число ядер процессора в одной системе.

Другое направление для решения вычислительных задач состоит в использовании в качестве вычислителя графического процессора GPU (graphics processing units). Технология использования графического процессора для общих вычислений GPGPU (general-purpose graphics processing units), которые обычно проводятся на центральном процессоре CPU (central processing units), является достаточно новой - первые «открытые» программные коды стали доступны в начале 2007 года. Технология GPGPU позволяет перейти на новый уровень в производительности вычислительной системы, что обусловлено тем, что, в отличие от CPU, GPU состоит из множества специализированных ядер, изначально создаваемых для расчетов и визуализации трехмерных сцен. Известно, что производительность одного ядра графического

процессора заметно ниже производительности ядра центрального процессора. Однако современные GPU могут объединять сотни таких ядер - тем самым суммарная вычислительная мощность одного GPU (порядка 1000 Гфлопс) оказывается значительно выше СРи (порядка 30 Гфлопс). К сожалению, достижение столь высокой производительности вычислений на GPU возможно лишь в узком классе задач, допускающих распараллеливание.

Отметим, что стоимость наиболее мощных GPU соизмерима со стоимостью высокопроизводительных СРи. Кроме того, специализированные программные средства позволяют реализовывать алгоритмы и общего назначения, пригодные к выполнению на GPU, на языке высокого уровня Си. В результате время вычислений с использованием технологии GPGPU может быть уменьшено в десятки и сотни раз по сравнению с «классическим» вычислением на СРи при минимальных финансовых затратах, включая затраты на переписывание кода под GPU.

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

ближенных аналитических выражений для интенсивности рассеянного излучения, как правило, ограничено определенными размерами области распространения, существенно зависящими от расстояния до источника излучения, и геометрией объекта.

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

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

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

Диффузионная флуоресцентная томография (ДФТ) представляет собой современный оптический метод визуализации внутренней структуры оптически неоднородных сред и объектов с использованием флуоресцентных маркеров [4], которые вводят в исследуемый объект перед его облучением. Эти маркеры поглощают сигнал и переизлучают его в разных участках спектра; на практике облучение исследуемого объекта происходит на длине волны возбуждения флуоресцирующих маркеров, а детектирование сигнала - на длине волны, отвечающей максимальной флуоресценции.

Постановка задачи

Характер распространения излучения в случайно-неоднородной среде, в первую очередь, обусловлен ее оптическими свойствами. Обычно свойства среды описывают следующими ос-

новными параметрами, зависящими в общем случае от частоты (длины волны) излучения [5]:

• коэффициент поглощения ца - величина, обратная расстоянию, на котором поток излучения ослабляется за счет поглощения в е раз;

• коэффициент рассеяния ц^ - величина, обратная расстоянию, на котором поток излучения ослабляется за счет рассеяния в е раз;

• параметр анизотропии среды g, который принимает значения от 0 до 1 (случай g=0 соответствует изотропному рассеянию, я=1 - полному рассеянию вперед);

• показатель преломления поЬ]-.

Для удобства введем коэффициент экстинк-ции (коэффициент взаимодействия или ослабления) ц = ц^ + ц, характеризующий полные потери коллимированного пучка при распространении в среде.

Достаточно строгое математическое описание процесса распространения немодулирован-ного света в рассеивающей среде может быть сделано с помощью стационарной теории переноса излучения (ТПИ). Теория переноса с успехом применяется при решении ряда практических задач из оптики биотканей. Основное стационарное уравнение ТПИ для монохроматического света имеет вид [4-6]:

1 д

+ пУ + ц ((Г)

V дt

К (г, п, t) =

Ц, (Г) 4л

IК (г,п', t)р(п',п, 7)00п, + Q(r, п, t),

(1)

где (г, п, t) - яркость светового поля (лучевая

интенсивность) в точке г в направлении п в момент времени t на длине волны X, V - ско-рость света в среде, Q(r, п, t)- объемная функция источников, dQ.n, - единичный телесный угол в направлении п'. Фазовая функция р(п',п,7) представляет собой функцию плотности вероятности рассеяния в направлении п' фотона, движущегося в направлении п, т.е. характеризует элементарный акт рассеяния. В дальнейшем ограничимся случаем, когда распределение вероятности рассеяния фотона обладает симметрией относительно направления падающей волны, т.е. фазовая функция зависит только от угла 0 между направлениями волновых векторов п и п'. В такой ситуации фазовую функцию можно описать эмпирической формулой Хеньи-Гринстейна [4]:

р(0,7) =

1 - Я 2(7)

(1 + я 2 (7) - 2я 2(7)со8 0)1

Отметим, что здесь параметр анизотропии g равен среднему значению косинуса угла рассеяния.

Исходя из (1), легко видеть, что задачу распространения излучения в ДФТ можно описать двумя уравнениями. Первое из них описывает распространение излучения на частоте длины волны возбуждения 7ех, второе - на частоте длины волны флуоресценции 7ет соответственно. В стационарном случае

[пУ + ц, (7 ех )Кх (г, п) = х

4л х 11ех (г,п',Г)р(п',п,7ех)ЛПп' + 6(г,п),

4л

[пУ + Ц, (7 ^ )Кт (Г, п) = ^1^ X

4л

ц

Х 14т (Г, П', ,)Р(п\ П, 7ет №п' + ^ЛФ(Г)Еех (Г)

(2)

(3)

где Ф(г)- функция концентрации флуорофора, П - квантовый выход флуоресценции, Еех (г) = | Ьех (г, п^Оп - облученность светового

4л

поля. Естественно, что в общем случае систему (2)-(3) необходимо дополнить граничными условиями. Как правило, в большинстве случаев считают границу G частично прозрачной (полагая среду за границей объекта однородной и нерассеивающей с показателем преломления поШ) с коэффициентом отражения R, зависящим от угла падения а, отвечающего направлению волнового вектора п. В этом случае для отраженной от границы интенсивности излучения, считая, что Я не зависит от длины волны, будем иметь:

(г, п)| о = ЯЬ^т (г, п')| о, (4)

где коэффициент отражения вычисляется из закона Френеля [7]:

R =

у cos а — (у^їпа)2

у cos а + д/1—(у^па)2

где у = поШ / поЬі; углы падения а и отражения а ' = л — а лежат в одной плоскости, а есть угол между волновым вектором п и нормалью п а к поверхности а.

Система (2)-(4) полностью описывает задачу распространения излучения в ДФТ: распространение света от источников 0(г,п) на длине волны накачки и его переизлучение на длине волны флуоресценции. Задача состоит в определении функций Ьх(г, п) и Ьет (г, п) или их интегральных характеристик Еех (г) и Ет (г) .

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

Решение методом Монте-Карло задачи ДФТ

Описывать движение фотона будем шестью основными параметрами: тремя координатами -вектором т=(х,уг) и текущим направлением движения - единичным вектором п=( цх, ц , цг)

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

Алгоритм расчета одной реализации траектории движения фотона состоит из нескольких шагов (далее := обозначает операцию присвоения):

1. Запуск фотона

Происходит инициализация начального положения (х0, у0> z0) , направления движения фотона (цх, цу, цг), исходя из функции Q(r,п,,), его энергии Ж и статуса S:

х := х0 цх := цхс>

У := Уо цу :=цу0>

z := zо цг := цг0>

Ж := 1 S := 0.

2. Перемещение фотона

Выполняется прямолинейное перемещение в соответствии с текущим направлением движения (цх, ц , цг) фотона. Длина пробега фотона определяется из закона Буге-Ламберта-Бера [6]:

ц,

х := х + 5 • ц х,

у := у +5 •ц у,

г := г + 5 • ц ъ,

где £ - случайное число, подчиненное равномерному распределению на интервале (0;1).

2

3. Учет текущего положения фотона

3.1. Пересечение с границей объекта

В случае пересечения траектории движения фотона с границей объекта G:

• Вычисляется угол падения аі.

• Вычисляется коэффициент отражения R.

• Если £ < R, то фотон отражается, в противном случае - преломляется.

• Пересчет новых направлений ц х, ц у, ц г.

3.2. Попадание в область флуоресценции Проверяется попадание фотона в область

флуорофора. Если 5=0, то при попадании в область флуорофора 5:=1, и дальнейшие расчеты проводятся для параметров среды y.a,g на

длине волны флуоресценции; в этом случае переходят к шагу 5.

4. Поглощение и рассеяние фотона Изменение энергии фотона на каждом шаге

вычисляется по следующему закону:

W := WW.

Ц,

Угол рассеяния 0 определяется из фазовой функции Хеньи- Гринстейна:

cos 0 := <

1

(

(

1 + g2 -

1-g

г

г

1 - g + 2g\

g * о,

г g

V у

[2£, -1, g = о.

Новые направления движения вычисляются

следующим образом:

цХ := <

sin0(ЦхЦz cos Ф-Цу sinФ) I I ,

Цxcos0 +-------і-----.----- y, ЦІФ 1,

Ц2

Ц У :=

Цz :=і

sin 0 cos Ф, Цz = 1,

sin0(ЦЦ cos Ф - Цx sinФ) І І ,

ЦУ cos 0+--------^—;==-------, Ы Ф 1,

Ц2

sin 0 sin Ф, |цz| = 1,

Ц cos 0-^11 - Ц2 sin 0 cos Ф, |^| Ф 1,

■p-^- cos 9, L I = 1,

JLz\

где ф - случайное число, подчиненное равномерному распределению на интервале (0;2п).

Далее переходят к шагу 2.

5. Переизлучение фотона

В случае, если фотон на длине волны возбуждения попадает в область флуоресценции, то происходит переизлучение фотона. При этом переизлучение зависит от свойств того или иного флуорофора, но зачастую оно является изотропным:

W :="лФ(х, y, z)W, cos 9 := 2Е, -1,

цx := sin 0 cos ф, ц := sin 0 sin ф,

Цz := cos 0,

где ф и £ - случайные числа, подчиненные равномерному распределению на интервалах (0;2п) и (0;1) соответственно. Далее переходят к шагу 2.

Шаги 2-5 повторяют до тех пор, пока фотон не вылетит за границу исследуемого объекта или энергия фотона W станет пренебрежимо малой (порядка 10-8). Далее описанная процедура повторяется. Как правило, на практике обычно число случайных реализаций траекторий фотона составляет 108—1011. После этого полученные данные статистически обрабатываются. В результате такого моделирования могут быть определены различные величины, в том числе яркости световых полей Lx (r, n), Lem (r, n) и их

облученности Eex (r), Eem (r) в каждой точке объекта. Например, для определения яркости в заданной точке необходимо сложить энергии всех фотонов на соответствующей длине волны для заданного направления рассеяния.

3. Компьютерная реализация метода

Алгоритм метода Монте-Карло был написан на языках программирования C++, C++ c использованием технологии CUDA1 и MatLab R2009.

Особенности реализации приведенного выше алгоритма для GPU.

1. Запуск ядра на GPU

Вычислительное ядро запускается параллельно на 1024 блоках по 128 нитей в каждом (итого 131072 нитей). Каждое вычислительное ядро просчитывает 100 траекторий фотонов, таким образом, уменьшается латентность вычислительного процесса. В результате за один вызов просчитывается примерно 1.3-107 различных траекторий фотонов.

2. Генерация псевдослучайных чисел

В связи с тем, что для графического процесс-сора отсутствует библиотечный генератор случайных чисел (ГСЧ), встает необходимость его реализации под GPU. Для генерации псевдослучайных чисел используется алгоритм, основанный на линейном конгруэнтном методе.

3. Использование «быстрых» функций (таких как____sinf(), _logf() и пр.), которые явля-

ются более быстрой, но менее точной реализацией на GPU соответствующих стандартных функций.

1 В качестве среды программирования использовалась среда Visual Studio 2008 Express Edition, поддержка технологии CUDA осуществлялась с использованием CUDA Toolkit 2.3.

Таблица

GPU CPU C++ CPU MatLab

10' траекторий 0.84 354 315

108 траекторий 8.8 3614 3172

4. Результаты численного моделирования

При моделировании определялись величины Ex (r), Eem (r). В качестве исследуемого выбран объект прямоугольной формы с размерами 20 ммх20 ммх15 мм. Оптические показатели среды: на длине волны возбуждения ^ = 3 мм-1, ца = 0.02 мм-1, g = 0.7. На длине волны флуоресценции: = 2 мм-1, ца = 0.01 мм-1, g = 0.7.

Флуорофор: координаты центра — (10,10,10) мм, радиус — 1 мм, п =1. Координаты «точечного» источника фотонов — (10,9,0) мм. Начальное направление движения фотонов — (0,0,1). Считалось, что границы объекта прозрачные — R = 0.

В таблице представлены времена (в секундах), затраченные на моделирование метода Монте-Карло на CUDA (GPU — Nvidia GTX 260, 192 ядра, 576 МГц), Microsoft C++ (CPU — AMD Phenom II x4 920, 2.7 ГГц, однопотоковая реализация), MatLab (CPU — AMD Phenom II x4 920, 2.7 ГГц, векторная реализация).

Как видно из таблицы, использование GPU позволило почти в 350 раз уменьшить время расчета задачи распространения света внутри объекта методом Монте-Карло. При этом реализованная производительность для GPU составила 170 Гфлопс (23% от пиковой производительности 715 Гфлопс), для CPU — 0.5 Гфлопс (7% от пиковой производительности одного ядра CPU 6.75 Гфлопс). Отметим, что повышение эффективности программного кода под CPU потребовало бы кардинального изменения исходного кода (например, использования технологии SSE), т.е. фактически — разработки нового программного обеспечения. Но даже в этом случае при достижении пиковой производительности CPU в 27 Гфлопс (для четырех ядер) выигрыш от использования GPU составил бы 6 раз.

Заключение

Разработан программный комплекс для решения задачи распространения оптического излучения в диффузионной флуоресцентной томографии методом Монте-Карло, выполняемый на графическом процессоре. Использование графического процессора позволило почти в 350 раз уменьшить время расчета по сравнению с аналогичной реализацией на центральном процессоре. Отметим, что такой прирост произ-

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

В дальнейшем применительно к задаче ДФТ планируется использовать GPU с целью расчета распространения света внутри объекта со «сложной» геометрией; проведения имитационного моделирования экспериментов; восстановления оптических свойств среды (коэффициентов поглощения, рассеяния, параметра анизотропии) и т.д.

Работа выполнена при финансовой поддержке грантов РФФИ №№ 10-02-01109, 10-02-00744, Государственных контрактов с Федеральным агентством по науке и инновациям №№ 16.512.11.2140, 02.740.11.0839, грантов Президента Российской Федерации для поддержки молодых ученых №№ МК-1127.2010.2 и 16.120.11.1909-МК, гранта Правительства РФ для государственной поддержки научных исследований, проводимых под руководством ведущих ученых в российских ВУЗах № 11G.34.31.0017, и программы фундаментальных исследований Президиума РАН «Основы фундаментальных исследований нанотехнологий и наноматериалов».

Список литературы

1. Wang L., Jacques S.L. and Zheng L. MCML-Monte-Carlo modelling of light transport in multilayered tissues // Comput. Methods Programs Biomed. 1995. V. 47(2). P. 131-46.

2. Fang Q. and Boas D.A. Monte-Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units // Opt. Express. 2009. V. 17(22). P. 20178-90.

3. Alerstam E., Svensson T. and Andersson-Engels S. Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration // J. Biomed. Opt. 2008. V. 13(6). P. 060504.

4. Тучин В.В. Лазеры и волоконная оптика в биомедицинских исследованиях. Саратов: Изд-во Саратовского ун-та, 1998.

5. Тучин В.В. Оптическая биомедицинская диагностика. Том I. М.: Физматлит, 2007.

6. Исимару А. Распространение и рассеяние волн в случайно-неоднородных средах. М.: Мир, 1981.

7. Борн М., Вольф Э. Основы оптики. М.: Наука, 1973. 720 с.

USING GRAPHICS PROCESSORS TO SOLVE LIGHT PROPAGATION PROBLEM IN DIFFUSE FLUORESCENCE TOMOGRAPHY BY MONTE CARLO METHOD

I.I. Fiks

Monte Carlo based software has been developed to solve light propagation problem in diffuse fluorescence tomography using a graphics processor. Using a graphics processor has made it possible to significantly reduce the calculation time (almost 350 times) as compared to a similar realization on a central processing unit.

Keywords: mathematical simulation, light propagation, CUDA, parallel computing, Monte Carlo method.