CC BY
20
Разработка векторного алгоритма по технологии CUDA для трехмерного моделирования процесса лазерной коагуляции сетчатки
А.С. Широканев12, Н.А. Андриянов 3, Н.Ю. Ильясова1,2 1 Самарский национальный исследовательский университет имени академика С.П. Королёва, 443086, Россия, г. Самара, Московское шоссе, д. 34, 2 ИСОИ РАН - филиал ФНИЦ «Кристаллография и фотоника» РАН, 443001, Россия, г. Самара, ул. Молодогвардейская, д. 151, 3 Финансовый университет при Правительстве Российской Федерации, 125167, Россия, г. Москва, Ленинградский пр-т, д. 49
Аннотация
Для лечения диабетической ретинопатии в современной практике применяется лазерная коагуляция. В процессе лазерной операции параметры лазерного воздействия подбираются вручную врачом, что требует от врача достаточного опыта и знаний, чтобы достичь терапевтического эффекта. На основе математического моделирования процесса лазерной коагуляции можно оценить основные параметры без проведения операции. Однако сетчатка имеет достаточно сложную структуру, и при применении даже низкозатратных численных методов для моделирования требуется значительное время для получения результата. В связи с этим разработка эффективных по времени алгоритмов трехмерного моделирования является актуальной задачей, поскольку применение таких алгоритмов позволит обеспечить проведение комплексного исследования в рамках ограниченного времени.
В настоящей работе проводится исследование времени выполнения алгоритмов, реализующих различные вариации применения метода расщепления и метода конечных разностей, адаптированных под поставленную задачу теплопроводности, выявляется наиболее эффективный алгоритм, который далее подвергается векторизации и реализации с использованием технологии CUDA. Исследование проводилось с использованием Intel Core i7-10875H и Nvidia RTX 2080 MAX Q и показало, что аналог векторного алгоритма, ориентированного на решение многомерной задачи теплопроводности, обеспечивает ускорение не более, чем в 1,5 раза, по сравнению с последовательным вариантом. Разработанный векторный алгоритм, ориентированный на применение метода прогонки по всем направлениям трехмерной задачи, существенно снижает временные издержки, затрачиваемые на копирование в память видеокарты, и обеспечивает 40-кратное ускорение по сравнению с последовательным алгоритмом трехмерного моделирования. На основе такого же подхода разработан параллельный алгоритм математического моделирования, который обеспечил 20-кратное ускорение при полной загрузке процессора.
Ключевые слова: диабетическая ретинопатия, лазерная коагуляция, математическое моделирование, уравнение теплопроводности, параллельные алгоритмы, векторные алгоритмы, CUDA.
Цитирование: Широканев, А.С. Разработка векторного алгоритма по технологии CUDA для трехмерного моделирования процесса лазерной коагуляции сетчатки / А.С. Широканев, Н.А. Андриянов, Н.Ю. Ильясова // Компьютерная оптика. - 2021. - Т. 45, № 3. - С. 427-437. - DOI: 10.18287/2412-6179-C0-828.
Citation: Shirokanev AS, Andriyanov NA, Ilyasova NY. Development of vector algorithm using CUDA technology for three-dimensional retinal laser coagulation process modeling. Computer Optics 2021; 45(3): 427-437. DOI: 10.18287/2412-6179-C0-828.
Введение
В настоящее время наблюдается широкое применение лазерного излучения для решения различных прикладных задач: в передовых научных исследованиях, например, при термоядерном синтезе, в промышленности, системах связи и военной сфере. Однако одной из важнейших сфер применения технологий лазерного излучения является медицина [1 - 4]. Одним из современных инструментов оперирования при отёках и кровоизлияниях сосудов является лазерная коагуляция. Лазерная коагуляция является эф-
фективным способом лечения диабетической ретинопатии при ранней диагностике заболевания.
Самым важным фактором, определяющим риск развития ретинопатии, является продолжительность сахарного диабета. Сахарный диабет поражает кровеносные сосуды сетчатки, потенциально приводя к развитию диабетической ретинопатии. Несмотря на то, что поражению подвергаются все части сетчатки, наиболее опасным является поражение её центральной части - макулярной области. Возникновение макуляр-ного отёка сопровождается снижением остроты зре-
ния, искривлением видимых предметов (метамор-фопсии). Длительное развитие диабетической ретинопатии приводит к необратимым последствиям, поэтому очень важно обнаружение и лечение заболевания на ранних стадиях, иначе последствия могут привести к полной потере зрения [5 - 8]. Точная и ранняя диагностика, наряду с адекватным лечением, может предотвратить потерю зрения более чем в 50 % случаев [3, 9].
Лазерная коагуляция является стандартом для лечения диабетической ретинопатии [10], её эффективность была подтверждена в ходе крупного исследования (ЕТОИ^, 1987) [11, 12]. При воздействии на сетчатку лазерное излучение приводит к образованию коагулята, препятствующего кровоизлияниям сосудов. При этом необходимо подбирать оптимальные параметры лазерного воздействия (мощность, длительность импульса), чтобы излишне не повреждать сетчатку.
При лечении диабетической ретинопатии используется излучение с длиной волны 532 нм, соответствующей зеленому свету [13]. Применяемое при лечении оборудование позволяет врачу выбирать параметры лазерного изучения, важнейшими из которых являются длительность импульса и мощность. Микроожоги наносятся в соответствии со способом, выбранным в ходе анализа строения глазного дна, в том числе строения сетчатки. Тем не менее большинство современных методов лечения зачастую основываются на ручном способе нанесения микроожогов [14].
Современное оборудование МЛУ1ЬЛ8 позволяет обстреливать лазером области с целью образования коагулятов в полуавтоматическом режиме, однако для этого перед лечением требуется предварительный план расположения коагулятов, как правило, формируемый врачом непосредственно вручную. Такой способ имеет ряд существенных недостатков, среди которых можно выделить необходимость наличия у врача достаточного опыта, позволяющего корректно локализовать патологические и анатомические элементы на изображении глазного дна.
Для автоматизации формирования такого плана предлагаются цифровые методы [15 - 22], однако в настоящее время не существует алгоритмов полной автоматизации решения указанной задачи [14, 23]. Важно отметить, что в таких методах зачастую используются алгоритмы распознавания образов и статистического анализа данных [15, 18 - 24].
Учитывая вышесказанное, актуальной задачей является реализация алгоритмов для автоматического формирования плана коагулятов. Для этого необходимо выполнить адекватную оценку параметров излучения.
При этом для формирования предварительного плана необходима оптимальная оценка как минимум 2 параметров: мощности лазера при обстреле заданной области и минимального расстояния между коагулятами, обеспечивающего предотвращение излиш-
него повреждения сетчатки [14, 24] и повышения эффективности плана лазерной коагуляции сетчатки.
Оценка указанных параметров может быть выполнена путём математического моделирования лазерного воздействия на глазное дно. Моделирование лазерного воздействия сводится к решению трёхмерной задачи теплопроводности, поскольку интерес представляет изменение распределения температуры с течением времени. Для решения трёхмерной задачи эффективнее использовать метод расщепления [26], который позволит рассматривать устойчивые методы для решения задач меньшей размерности. Решение трёхмерной задачи обеспечивается алгоритмами с высокой вычислительной сложностью, в связи с чем применяют параллельные алгоритмы [27, 28], ориентированные на загрузку процессора (CPU), и векторные алгоритмы [29] по технологии CUDA, ориентированные на загрузку графического устройства (GPU), для ускорения последовательных алгоритмов. Также применяют многопроцессорные системы с распределенной, общей памятью и системы гибридного типа [27 - 29].
В настоящей работе решается задача математического моделирования распространения тепла в трёхмерном пространстве с использованием параллельных и векторных алгоритмов, реализующих метод расщепления и метод конечных разностей. При помощи разрабатываемых алгоритмов предлагается исследовать влияние анатомических и патологических элементов на сетчатке глаза на радиус лазерного воздействия, обеспечивающего терапевтический эффект в зоне поражения, и анализ влияния плана коагуляции на результирующий терапевтический эффект.
1. Трёхмерная задача математического
моделирования процесса лазерной коагуляции на глазном дне
Поскольку решается задача моделирования распространения тепла, то в случае применения лазера необходимо учитывать модель световой энергии, которая преобразуется в тепловую [30 - 31]. При лазерном излучении формируется электромагнитное поле, энергия которого впоследствии преобразуется в тепло [32 - 33].
Интенсивность световой энергии на заданном расстоянии при известных мощности лазера и радиусе пятна может быть посчитана по формуле (1).
P
I (r ) = -P- e-(r /a )2, (1)
na2
где a - радиус пятна, P - мощность лазера.
Интерес представляет распределение температуры после прекращения лазерного воздействия, определяемое формулой (2). На практике длительность импульса лазерного воздействия At задается достаточно малой величиной, и поскольку эта величина достаточно мала, то при анализе можно пренебречь и дифракцией света [29 - 34].
ехр(-|р(х, у, %)а -р/ (г )Д/
х, у, г) =-о----+ Тс
Со!
(2)
где Тс - температура, сформированная в результате предыдущих выстрелов, р = р (х,у, г) - функция поглощения среды, С об = Сов (х, у, г) - функция объёмной теплоёмкости среды в зафиксированный момент времени;
г =у](х-х0)2 + (у-у0)2 + (г-г0)2 -
удаление от точки, куда был произведен выстрел, имеющей координаты (хо, уо, го).
Формула (2) характеризует начальное распределение температуры на глазном дне. В слое эпителия температура максимальная и со временем будет, очевидно, уменьшаться. Однако вследствие перераспределения температуры на соседние слои со слоя эпителия происходит повышение температуры соседних слоев. Таким образом, сетчатка, которая является соседним слоем эпителия, нагревается не только за счет лазерного воздействия, но и за счет перераспределения температуры со слоя эпителия. Поэтому важно учитывать изменение температуры на глазном дне с течением времени, то есть решать задачу математического моделирования. Также необходимо учитывать, на какое расстояние распределяется температура. Таким образом, сформулируем задачу моделирования лазерного воздействия на основе распространения температуры в виде (3).
дТ
Сов — = Лу(к - gгadxyг (Т)) д/
П=о=^(х, у,г) ,
Т Г =Т-
(3)
где Т = Т (х, у, г, /) - распределение температуры, к = к (х, у, г, Т) - функция температуропроводности среды, Соб = Соб (х,у, г, Т) - функция объёмной теплоёмкости среды, зависящая также от температуры, Г - граница воздействия лазера, div - дивергенция векторного поля, gгadxyг - градиент функции по пространственным координатам, То - температура, соответствующая нормальной жизнедеятельности живых тканей.
При использовании выражения (3) необходимо, чтобы область определения теплового поля в пространстве была достаточно большой, обеспечивая при лазерном воздействии распространение тепла лишь до той границы, через которую проходит лазер. Для достижения этого могут быть использованы два способа.
Во-первых, функция Т (х, у, г) должна быть задана на несущей границе, т.е. такой, которая изменяется со временем. Кроме того, необходимо выполнить специальную линейную замену для того, чтобы свести поставленную задачу к фиксированным (или нулевым) граничным условиям. Таким образом, теплопровод-
ность будет описываться неоднородным дифференциальным уравнением. Тогда решение задачи возможно с разложением на решение соответствующего однородного уравнения с исходными граничными и начальными условиями. Далее находится решение неоднородного уравнения, для которого обеспечиваются нулевые граничные и начальные условия. Следовательно, окончательное решение будет находиться из уравнения вида (3).
Тем не менее, в общем случае при поиске такой замены и решении неоднородного уравнения возникает ряд сложностей, поэтому в данной работе было принято решение использовать второй способ. Такой способ базируется на том, что сперва требуется искусственно расширить область определения, а затем на основе полученного расширения разделить её на информативную и неинформативную области. При этом для первой области выполняется условие относительно пренебрежимо малой длительности импульса Д/, что позволяет описать распределение температур в начальный момент времени в соответствии с выражением (2). Определение распространения тепла в неинформативной области не составит труда, поскольку для решения данной задачи достаточно будет использовать симметричное отображение. Однако границе в неинформативной области будет соответствовать уже фиксированное значение температуры. При этом при анализе математической модели можно будет использовать только данные в информативной области, что упрощает задачу.
К сожалению, в силу зависимости объёмной теплоёмкости и теплопроводности от температуры поставленная задача будет нелинейной. Тем не менее возможно оценить изменение формы сетчатки, опираясь на температурные значения её слоёв, т.е. оценивая, какой слой и до какой температуры нагрелся. Рассмотрим аналогичную задаче (3) задачу (4), в которой функции объёмной теплоёмкости Соб и теплопроводности к зависят только от пространственных координат.
дТ
Сов (х, у, г)— = div(k (х, у, г) - gгadXyZ (Т)) д/
Т| = ехр(-Р( х, у, г) г) -р( х, у, г) / (г )Д/ +Т
|/= о ^ ' - о
Сов (х, у, г)
(4)
Т1Г =То
Как видно из уравнения (4), перейти к нулевым граничным условиям можно с помощью линейной замены Т = Т + То. Такой подход позволяет упростить решение задачи и в тексте статьи далее будет решаться задача ускорения математического моделирования именно в случае нулевых граничных условий. На выходе математической модели будет формироваться многомерное поле, отображающее, как сильно нагрелась ткань в той или иной точке пространства после
лазерного воздействия. Поскольку моделируется непосредственно изменение температуры (Т ), вызванное лазером, то определение абсолютной температуры ткани легко получить, прибавив к смоделированному значению изменения температуры нормальную температуру То (~ 36,5 °С).
Кроме того, следует отметить, что при моделировании ситуации (4) учитывается процесс изменения поля температур с течением времени. Важно, чтобы на слое сетчатки температура не превышала критической температуры, которая приблизительно равна 38 - 40 °С. Для образования коагулята на сосудистом слое необходимо нагреть его до 39 °С.
Однако на самом деле сетчатка имеет более сложную структуру и не может рассматриваться как однородная среда (рис. 1). Поэтому решение задачи (6) может быть использовано лишь в качестве начального условия при моделировании распространения тепла на сложной структуре сетчатки. Также на рис. 1 представлены граничные условия и область лазерного воздействия.
Граничные условия
Рис. 1. Структура сетчатки глаза
Следует отметить, что рис. 1 представляет собой один из снимков сетчатки в сечении в плоскости оХ2, формируемых в результате оптической когерентной томографии (ОКТ). Обычно при одной регистрации накапливается 85 подобных снимков по разным координатам У. Используя данные снимки и зная координату У для каждого, можно реконструировать трёхмерную модель.
Следовательно, для решения задачи математического моделирования лазерного воздействия необходимо использовать численные методы. Метод конечных разностей часто используется при решении сложных задач математического моделирования [34, 35].
В частности, схема расщепления позволяет рассматривать одномерные задачи вместо многомерных, благодаря чему можно использовать неявную схему, не решая полноценную систему линейных уравнений, а применяя метод прогонки, что существенно быст-
рее. При этом при любом количестве входных данных желательно распараллеливать используемый метод, поскольку алгоритм для трёхмерного случая может выполняться крайне долго.
Таким образом, перепишем задачу (4) в виде задачи (5), для решения которой необходимо применить метод расщепления
ч дТ д т чТ соб (X, У, 2)-= — (к(X, у, 2) —) +
t дх дх
д дТ д дТ
+—(к (Х ^ 2)—) + —(к (Х У, z)^), (5)
ду ду д2 д2
Т =о =Ц>(X, У, 2), Т |г = То.
Благодаря методу расщепления численное решение исходной задачи представляется в удобной форме, которое может производиться алгоритмами с более низкой вычислительной сложностью. При использовании метода расщепления отрезок по времени подвергается равномерной дискретизации, а исходная задача (5) преобразовывается к итерационным задачам (6) и (7). Координата у отщепляется в первую очередь, поскольку снимки ОКТ располагаются вдоль оси у.
( ) д^ д(к )дЖ)
Соб (X ^ 2)— = — {к(X ^ 2)—X дt ду ду
= Т\,=_к, (6)
Ж \Г= 0.
соб(x, У, 2)— = — {к( X У, 2 ) —) +
д/,, \д¥. +—k(x, ^ 2)—) (7)
д^ д2
V |г= о.
Алгоритм применения метода расщепления следующий: сначала решается задача (6) на отрезке ^к, tk +1], в которой функция Ж в момент tk совпадает с искомой функцией Т в тот же момент tk; затем решается задача (7), в которой функция V в момент ^ является результатом моделирования функции Ж в момент ^ +1. В рамках метода расщепления в момент времени ^ +1 результат моделирования может быть приблизительно представлен функцией V, т.е.
Т| и VI .
Следует отметить, что задача (6) представляет собой, вообще говоря, набор задач, так как присутствует зависимость от всех пространственных координат. То же самое можно сказать про задачу (7). Задачу (6) рекомендуется решать методом конечных разностей с использованием неявной разностной схемы. Задача
(7) может быть решена разными численными методами, в том числе методом расщепления.
На основе уравнения задачи (6) с помощью метода конечных разностей получим дискретное уравнение в виде (8). Переход к дискретным координатам в непрерывном пространстве (х, у, г) и времени / осуществляется путем разбиения этих координат на счетное количество, которые будут составлять дискретную сетку со счётными индексами (, }, к) и 5, т.е. ху, гк, 4).
^ijk "
*ijk
h
vijk _ k wj+1k ~ wijk1 _ kijk
hi
- k j
Kij-1k
Vij-1k
hi
(8)
где к, ку - достаточно маленькие шаги по дискретной сетке времени / и оси у.
Наконец, перепишем полученное выражение (8) в виде (9)
f
vjk
1+h
hi
(1
w
ij-1k
V cijk
ijk
vij+1k
//
kijkht
h2c ''y^ijk
-wj+1k^ _ Aywj-k + Bywj+ + CyWj+1k _ wjk,
(9)
где Ay _ —
^ h2
ijk y
kij-1kht
c h2
ijk y
By_
1 +
kijk + kij-1k
S\
C _ ^^
Cy :
/у
hyCijk
Для упрощения полученной записи коэффициентов Ау, Ву, Су можно ввести следующие замены:
у _ ht D000 _ kijk D0-0 _ kj-1k
ly ~ TT' Dijk —-iDijk —
Cj
hi
ijk
Тогда коэффициенты Ay, By, Cy перепишутся в виде (10)
Ay _ -уyDj-o,
By _1 + уy(Dj0 + Dj-°),
Cy _ -у y А°Г,
(10)
где у коэффициента Б обозначение «-» сверху сообщает, что произведено смещение по соответствующей оси, а обозначение «о» характеризует, что смещения по соответствующей оси не было. Обозначения «-» и «о» выстраиваются в ряд из трех символов в соответствии с осями х, у, г.
Аналогичные уравнения вида (9) с учётом (10) могут быть получены и для координат х, г. Тогда можно записать все коэффициенты в системах (11) и (12).
Ax _ -уxD-™, Bx _1 + уx(Dj00 + D->°) Cx _ -уD
ijk
A; _ -уzDi
jk
Bz _1 + уz (Dj0 + D»-)
с. _ -уD
(11)
(12)
ijk
где у х = к, /кх2, у г =к, /кг2.
В дальнейшем для удобства записи будем рассматривать тензоры Бооо,Б - оо,Бо - 0,Б00 - в виде (13).
k , k
П000 _ ijk П-00 - i-
Dijk — ;Dijk — Ciik Ci
1 jk
П0-0 _ j-1k П00- _ Dijk —-;Dijk —
cijk
(13)
ijk
Рассматривая задачу в сечении слоистой структуры глазного дна, получаем двумерную задачу. Координата x соответствует смещению лазерного излучения, а координата z соответствует глубине.
Таким образом, трёхмерная задача сводится к решению уравнений (6) и (7). При этом двумерная задача также может быть разложена на одномерные задачи либо решена с использованием численных методов.
2. Численные методы решения задачи математического моделирования распространения тепла в структуре глазного дна
Сложные задачи математического моделирования часто решаются с использованием разностных схем [34 - 35]. Например, решение задачи (4) можно получить, используя такие разновидности данного метода, как:
1) применение явной разностной схемы;
2) применение неявной разностной схемы;
3) применение метода расщепления для решения однородных задач.
Самым простым способом является применение явной разностной схемы. Однако для обеспечения устойчивости метода требуется достаточно большое количество узлов. При этом время работы алгоритма возрастает из-за необходимости повторения множества итераций, хотя одна итерация и выполняется достаточно быстро.
Преимуществом неявной схемы является её устойчивость вне зависимости от числа шагов дискретизации. Но основной недостаток неявной схемы -высокая вычислительная сложность алгоритма на каждой итерации.
Наконец, метод расщепления используется именно для представления многомерных задач в виде одномерных задач. Это позволяет применять неявную схему без решения полноценной системы линейных уравнений, что приводит к значительному сокращению вычислительных затрат при применении метода прогонки.
Сравнительные характеристики работы описанных методов показаны в табл. 1 для размерности сетки 100 х 100. Расчеты выполнялись на базе AMD Ryzen 5 2600 Six-Core. При этом сравнение выполнялось для последовательных алгоритмов.
Обратим внимание, что результаты в табл. 1 получены при маленьких размерах сетки (100 х 100). Очевидно, что неявная схема работает значительно медленнее явной схемы и метода расщепления. В связи с
этим далее не будет рассматриваться явная схема. Рассмотрим работу на более реалистичных размерах сетки и количестве узлов.
Табл. 1. Результаты времени выполнения алгоритмов,
основанных на методе конечных разностей, при использовании различного количества узлов по времени
Алгоритм Время при различном количестве узлов, с
100000 10000 1000
Явная схема 6,58 0,59 0,08
Метод расщепления 239,20 22,52 2,33
Неявная схема 20607,55 1746,02 211,55
При этом стоит отметить, что для сеток размерностей 640x800 были получены следующие результаты, представленные в табл. 2. Следует также отметить, что для ускорения явной схемы на CUDA использовалась видеокарта NVIDIA GeForce GTX 1070 Ti. Алгоритм на CUDA для явной схемы был реализован с использованием Parallel Computing Toolbox среды Matlab. Явная схема описывается через матричные операции (покомпонентное умножение, сложение матриц), которые автоматически векторизуются с использованием упомянутого инструмента параллельных вычислений.
Табл. 2. Результаты времени выполнения явной схемы, явной схемы с ускорением и метода расщепления
Алгоритм Время, с
Явная схема, CUDA 405,08
Явная схема 6484,94
Метод расщепления 578,17
Важно, что в табл. 2 рассмотрены значения для оптимального для разных методов количества узлов. Следует отметить, что в данном случае в методе расщепления реализовывался только последовательный алгоритм. Анализируя эти значения, можно сделать 2 вывода. Во-первых, последовательный алгоритм, основанный на расщеплении, приводит к лучшим результатам по скорости, чем явная схема и неявная (тоже последовательные). Действительно, он оказывается быстрее методов конечных разностей явной и неявной схемы. Во-вторых, применение СИБЛ позволяет существенно ускорить явную схему и значительно превзойти по скорости метод расщепления.
Несмотря на эффективность использования СИБЛ-алгоритма, реализующего явную схему, по сравнению с последовательным алгоритмом, реализующим метод расщепления и решение одномерных задач с использованием неявных схем, метод расщепления может быть векторизован, благодаря чему сможет оказаться более эффективным по сравнению с рассмотренными алгоритмами.
3. Параллельные и векторные алгоритмы моделирования распространения тепла на структуре глазного дна
Для решения уравнения вида (6) используется метод прогонки [34]. При этом применение метода расщепления позволяет использовать метод прогонки для решения одномерной задачи в виде (9). Прогонка
должна проводиться по каждому направлению, однако при использовании GPU перераспределение данных перед каждым запуском единой реализации прогонки чревато выполнением затратных операций копирования в память видеокарты и обратно. В связи с этим эффективнее векторизовать сам метод расщепления, который разбивает трёхмерную задачу на наборы одномерных задач.
Следует отметить, что векторный алгоритм может быть реализован как алгоритм с использованием CUDA (векторный алгоритм), так и алгоритм, использующий все ядра процессора (параллельный алгоритм). Векторный алгоритм представляет собой псевдокод, описанный математически. Реализация векторного алгоритма для CUDA описывается функциями ядра на языке С++ и кодом запуска ядра, а параллельный алгоритм с учетом независимости вычислительных задач может быть описан на высокоуровневом языке программирования, например C#.
Для разработки векторного алгоритма, реализующего метод расщепления, рассмотрим в первую очередь отщепление переменной y от переменных (x, z), которое приводит к множеству однотипных одномерных задач, каждая из которых решается относительно оси oY с фиксированными значениями по другим координатам (x0, z0). В таком случае образуется совокупность независимых задач вида (6), которые могут быть распараллелены между ядрами процессора или графическими ядрами видеокарты. В дискретной форме происходит переход от координат (x, y, z) к координатам (i,j, k) и требуется решение однотипных задач вида (9) при фиксированных координатах (i 0, k 0).
Векторный алгоритм сводится к решению независимых систем линейных уравнений, описываемых трёхдиагональными матрицами, то есть к выполнению метода прогонки [36, 37] в отдельных вычислительных задачах.
Трёхдиагональная матрица может быть описана тремя векторами a, b и c, соответствующими диагоналям. Набор трёхдиагональных матриц будем описывать матрицами А, B и C, где i-й столбец соответствует вектору-диагонали (a, b или c соответственной i-й матрицы. Векторы правых частей будут храниться в соответствующих столбцах матрицы F.
Для решения трёхдиагональной системы в методе прогонки используются векторы а и р. Аналогично приведенной логике будем рассматривать матрицы а и р, соответствующие набору трёхдиагональных систем линейных уравнений. Отметим, что далее будут использоваться операции покомпонентного умножения (.*) и покомпонентного деления (./).
На первом шаге векторного алгоритма, реализующего метод прогонки, выполняются операции по формуле (14).
а(1,:) = -С (1,:)./ 5(1,:), ß(1,:) = F (1,:)./ B(1,:).
(14)
Здесь а и р характеризуют матрицы прогоночных коэффициентов. «:» означает, что операция выполняется для всех строк или столбцов (то есть всех индексов по соответствующему измерению).
На втором шаге выполняется векторная версия прямой прогонки по формуле (15).
a(i,:) =
-С (i,:).
ß(i,:) =
A(i, [ F (i,
A(i,
).* a(i -1,:) + B(i,:) ) - A(i,:).* ß(i -1,:)].
(15)
).* a(i -1,:) + B(i,:)
На третьем шаге векторно вычисляется последний коэффициент решения по формуле (16).
x( N ,:) =
[F (N,:) - A(N ,:).*ß(N -1,:)].
B(N,:) + A(N,:).*a(N-1,:) :
(16)
где N - размерность одной матрицы в системе.
На последнем шаге производится векторный вариант обратной прогонки по формуле (17).
x(i,:) = a(i,:).* x(i +1,:) + ß(i,:).
(17)
Данная задача должна рассматриваться по всем трём измерениям. Поэлементные операции в формулах (12 - 17) легко векторизуются, а так как систем линейных уравнений достаточно большое количество, то векторный метод прогонки эффективен при решении набора трёхдиагональных систем.
Векторный алгоритм, реализующий метод прогонки, эффективен при решении одномерной задачи, когда коэффициенты трёхдиагональных матриц не меняются. Однако при реализации метода расщепления необходимо на каждой итерации производить прогонку в трёх разных направлениях. Для этого требуется постоянно копировать в память видеокарты новые матрицы, а также перераспределять данные, что может сказаться на снижении эффективности алгоритма. В связи с этим для уменьшения необходимости излишних копирований в память видеокарты и обратно были предложены 3 алгоритма, реализующих применение неявных схем для решения одномерных задач, соответствующих разным направлениям при использовании метода расщепления.
Идея алгоритмов заключается в использовании тензоров вида (13) для вычисления на каждой итерации значений трёхдиагональных матриц. То есть алгоритм по любому направлению включает вычисление значений трёхдиагональных матриц и проведение на основе вычисленных значений метода прогонки.
Все 3 алгоритма идентичны, поэтому рассмотрим одно из трёх направлений - направление У. Для рассматриваемого направления входными данными являются двумерные индексы текущей задачи (/, к) и одномерный индекс задачи t, который должен образовывать биекцию с двумерным индексом, но при этом допускается использовать любое биективное
отображение. При использовании СИБЛ одномерный индекс t используется для матриц а и р.
Далее в тексте используются следующие обозначения: а и р - матрицы прогоночных коэффициентов; коэффициенты а, Ь, с вычисляются локально в рамках алгоритма; 3 - размерность сетки по оси У.
тЫ {/, к, ^
b = 1+1 у (0 + АУ) c = -1yAZ
a0t =-c / b ß0t = Xi0k / b
forj = 1: J - 2
a = -jyAy b =1 + i y ( + АУ)
c=-ЪА7
d = a ■ а j-11 + b а jt =-c / d ß jt =(Xjk - aß j-1,)/ d
a = -jyAJA k
b = 1 + 1 y ((J-1 k + А0-- k )
Xij-ik =(XU-1 k - aßj-21)/(a ■aj-2, + b)
forj = J-"275
[Xjjk =ajt ■Xij+1 k +ßjt
На основе приведенного псевдокода реализуется функция ядра для GPU на языке С++ для соответствующего направления. В результате были реализованы 3 функции ядра, соответствующие описанным алгоритмам. Функции ядра запускаются поочередно в рамках каждой итерации. При этом функции ядра на каждой итерации обеспечивают изменение результирующего тензора X без необходимости копирования данных в память видеокарты и обратно. После проведения необходимого количества итераций по времени в тензоре X будет храниться результат моделирования в требуемый момент времени.
4. Экспериментальные исследования алгоритмов моделирования процесса лазерной коагуляции
Исследование разработанных алгоритмов проводилось на основе вычисления ускорения - отношения времени выполнения последовательного алгоритма ко времени параллельного или векторного. Для получения наиболее достоверных результатов будем вычислять среднее время одной итерации. Для этого будем замерять время проведения 30 итераций, а затем усреднять по ним. Для запуска параллельного алгоритма использовался процессор Intel Core i7-10875H, для запуска векторного алгоритма на CUDA - видеокарта Nvidia RTX 2080 MAX Q.
В табл. 3 показано среднее время итерации для различных алгоритмов на сетке 150x150x500. В данной таблице последовательный файловый алгоритм представляет собой алгоритм, основанный на использовании внешней памяти для хранения всех тензоров, чтобы была возможность проводить моделирование на крупных трёхмерных сетках, не влезающих в оперативную память. Однако такой алгоритм производит вычисления недопустимо долго. Последовательный алгоритм RAM использует оперативную память для проведения вычислений. Алгоритм, основанный на векторизации прогонки, применяет функцию ядра CUDA, соответствующую векторной прогонке. Векторизация метода расщепления реализована в виде параллельного алгоритма (полной загрузки CPU) и векторного алгоритма (полной загрузки GPU).
На рис. 2 продемонстрирована зависимость ускорения параллельного (штриховая линия) и векторного (сплошная линия) алгоритмов от объёма памяти, занимаемой трёхмерными сетками.
Табл. 3. Результаты времени выполнения для различных алгоритмов
Алгоритм Время выполнения, мс
Последовательный файловый 75G34
Последовательный RAM 24S1
Векторизация прогонки, CUDA 1611
Векторизация метода расщепления, CPU 136
Векторизация метода расщепления, CUDA 67
Анализ результатов, представленных в табл. 3, показывает, что векторизация прогонки не обеспечивает достаточную эффективность: копирование данных в память видеокарты сильно замедляет алгоритм. Однако алгоритм, основанный на векторизации метода расщепления, обеспечивает существенное ускорение, несмотря на дублирование операций вычисления значений трехдиагональных матриц, по сравнению с последовательным алгоритмом. Параллельный алгоритм также обеспечивает большое ускорение за счёт отсутствия операции перераспределения данных под разные направления.
Ускорение
50-
40 30 20 10 0
Г
-•J
Intel Core i7-10875H NVidia RTX 2080 MAX Q
0,57 1,84
3,12
4,4
5,68 6,96 8,24 Размерность, ГБ
Рис. 2. Результаты ускорения в зависимости от объема данных при использовании разработанных алгоритмов
Из представленных результатов можно сделать вывод, что векторный алгоритм стабильно в 2 раза
быстрее параллельного. При низких размерностях сетки, что эквивалентно малым объёмам памяти, не наблюдается ярко выраженной зависимости ускорения для обоих алгоритмов. Однако при переходе к большим размерностям сеток ускорение стабилизируется и не изменяется при увеличении размерности сетки (объёма памяти).
Следует отметить, что параллельный алгоритм обеспечивает ускорение на уровне 20 (рис. 2), что выше теоретического, которое для процессора Intel Core i7-10875H составляет 12. За счет векторизации метода расщепления, а не метода прогонки, наблюдается сокращение издержек на работу с памятью и существенное ускорение алгоритма. Как демонстрируют результаты на рис. 2, использование GPU и CUDA позволяет обеспечивать ускорение на уровне 40 по сравнению с последовательным алгоритмом.
При моделировании лазерного воздействия следует учитывать возможность смещения лазера. Для попадания в нужную точку сетчатки луч лазера наводится под углом, который обеспечит прижигание лазером нужного участка. Моделирование лазерного излучения при наведении лазера в центр сетчатки является самым простым случаем. Для такого случая закон Бугера действует вдоль оси z. На рис. 3 демонстрируется результат моделирования лазерного воздействия на трёхмерной структуре сетчатки. Интенсивность белого цвета характеризует величину температуры в рассматриваемой точке. Для наглядности распределение температуры в трёхмерном случае наложено на перпендикулярные центральные сечения структуры глазного дна (рис. 3).
Рис. 3. Результат моделирования лазерного воздействия на глазное дно в трёхмерном пространстве
При задании угла лазерного воздействия вводится прямая, вдоль которой распространяется лазерное излучение. Алгоритм вычисления начального распределения температуры при задании угла является более сложным, чем алгоритм вычисления начального распределения температуры вдоль оси г. На рис. 4 представлено распределение температуры для трёх разных
углов наклона направления лазерного излучения, которые вычисляются по значениям смещения в плоскости oXZ. На рис. 4 интенсивность белого цвета характеризует величину температуры в соответствующих точках. Для наглядности распределения температур наложены на снимки ОКТ.
Анализ представленных результатов моделирования (рис. 4а- в) позволяет сделать следующие ключевые выводы. Во-первых, независимо от формы сетчатки в целом, температура распространяется примерно одинаково, т.е. температура нагрева незначительно зависит от формы сетчатки, которая различна для рассмотренных смещений. При этом следует отметить, что толщина сетчатки для рассмотренных случаев приблизительно одинаковая. Во-вторых, наблюдается некоторое незначительное различие в рассеивании тепла при разных смещениях. Такая раз-
ница становится важной, когда происходит не одиночная коагуляция, а, например, формируются 2 коагулята рядом. В таком случае возникает необходимость учитывать смещение, поскольку в некоторой области сетчатки суммарное воздействие двух выстрелов может привести к преувеличению критической температуры, что недопустимо.
По представленным рисункам и реализованным моделям можно оценить один из важнейших параметров моделирования: на какое расстояние распространяется тепло. Также при помощи математического моделирования возможно оценить безопасное расстояние между коагулятами. Это, как было отмечено ранее, позволит использовать математическое моделирование как полезный инструмент при определении необходимых параметров излучения и построении плана коагулятов.
Рис. 4. Коагуляция при различном смещении лазерного воздействия относительно оси оХ: 200 мкм (а), 350 мкм (б), 500 мкм (в)
Заключение
В настоящей работе были разработаны и исследованы алгоритмы численного моделирования распространения тепла по результатам лазерного воздействия на глазное дно на основе метода расщепления и метода конечных разностей. Разработанные векторные и параллельные алгоритмы позволили значительно повысить скорость вычислений по сравнению с последовательными алгоритмами. Такое моделирование позволяет оценить основные параметры, необходимые при реальной практике лечения диабетической ретинопатии.
Последовательные алгоритмы реализации численных методов, в частности рассмотренного в статье метода расщепления, приводили к крайне долгому процессу реализации модели даже на современных устройствах вычислительной техники. Задачу ускорения последовательных алгоритмов можно решать при поиске эффективных параллельных алгоритмов, ориентированных на нагрузку всех ядер процессора CPU, или векторных алгоритмов, ориентированных
на работу с графическими устройствами GPU с использованием CUDA-ядер.
Полученные результаты моделирования и оценки ускорения показали, что векторный алгоритм, реализованный на видеокарте Nvidia RTX 2080 MAX Q, обеспечивает приблизительно двукратное ускорение по сравнению с параллельным алгоритмом, реализованным на процессоре Intel Core i7-10875H. При этом при увеличении размерностей сетки ускорение переходит в постоянное значение.
Разработанные алгоритмы моделирования лазерного воздействия для формирования распределения температуры могут быть использованы для повышения эффективности лазерной коагуляции сетчатки благодаря подбору оптимальных параметров лазерного воздействия.
Благодарности
Исследование выполнено при финансовой поддержке РФФИ в рамках научных проектов № 19-3190160, № 19-29-01135 и Министерства науки и высшего образования Российской Федерации в рамках
выполнения государственного задания Самарского
университета и ФНИЦ «Кристаллография и фотоника» РАН.
Литература
1. Гафуров, С.Д. Особенности применения лазеров в медицине / С.Д. Гафуров, Ш.М. Катахонов, М.М. Холмонов // European Science Journal. - 2019. -№ 3(45). - С. 92-95.
2. Коцур, Т.В. Эффективность лазерной коагуляции в ма-куле и микрофотокоагуляции высокой плотности в лечении диабетической макулопатии / Т.В. Коцур,
A.С. Измайлов // Офтальмологические ведомости. -2016. - Т. 9, № 4. - C. 43-45. - DOI: 10.17816/0V9443-45.
3. Замыцкий, Е.А. Лазерное лечение диабетического ма-кулярного отека / Е.А. Замыцкий // Аспирантский вестник Поволжья. - 2015. - № 1-2. - С. 74-80.
4. Kozak, I. Modern retinal laser therapy / I. Kozak, J. Luttrull // Saudi Journal of Ophthalmology. - 2014. - Vol. 29, Issue 2. - P. 137-146.
5. Doga, A.V. Modern diagnostic and treatment aspects of diabetic macular edema / A.V. Doga, G.F. Kachalina, E.K. Pedanova, D.A. Buryakov // Ophthalmology, Diabetes. - 2014. - No. 4. - P. 51-59.
6. Whiting, D.R. IDF diadetes atlas: global estimates of the prevalence of diabetes for 2011 and 2030 / D.R. Whiting [et al.] // Diabetes Res. Clin. Pract. - 2011. - Vol. 94(3). - P. 311-321.
7. Братко, Г.В. К вопросу о ранней диагностике и частоте встречаемости диабетического макулярного отека и формировании групп риска его развития / Г.В. Братко,
B.В. Черных, О.В. Сазонова // Сибирский научный медицинский журнал. - 2015. - Т.35, №1. - С. 33-36.
8. Воробьёва, И.В. Диабетическая ретинопатия у больных сахарным диабетом второго типа. Эпидемиология, современный взгляд на патогенез. Обзор / И.В. Воробьёва, Д.А. Меркушенкова // Офтальмология. - 2012. - Т. 9, № 4. - С. 18-21. - DOI: 10.18008/1816-5095-2012-4-18-21.
9. Амиров, А.Н. Диабетический макулярный отёк: эпидемиология, патогенез, диагностика, клиническая картина, лечение / А.Н. Амиров, Э.А. Абдулаева, Э.Л. Минхузи-на // Казанский медицинский журнал. - 2015. - Т. 96, №1. - С. 70 - 74.
10. Астахов, Ю.С. Современные подходы к лечению диабетического макулярного отека / Ю.С. Астахов, Ф.Е. Шадричев, М.И. Красавина, Н.Н. Григорьева // Офтальмологические ведомости. - 2009. - Т. 2, № 4. -
C. 59-69.
11. Исхакова, А.Г. Результаты клиникоэкономического анализа лечения больных диабетической ретинопатией с макулярным отеком / А.Г. Исхакова // Аспирантский вестник Поволжья. - 2014. - № 1. - С. 96 - 98.
12. Уманец, Н.Н. Интравитреальное введение ранибизумаба как метод лечения больных кистозным диабетическим ма-кулярным отеком / Н.Н. Уманец, З.А. Розанова, М. Альзин // Офтальмологический журнал. - 2013. - № 2. - С. 56-60.
13. Cohen, S.M. Laser energy and dye fluorescence transmission through slood in vitro / S.M. Cohen, J.H. Shen, W.E. Smiddy // American Journal of Ophthalmology. -1995. - Vol. 119, Issue 4. - P. 452-457.
14. Замыцкий, Е.А. Анализ интенсивности коагулятов при лазерном лечении диабетического макулярного отека на роботизированной лазерной установке Navilas / Е.А. Замыцкий, А.В. Золотарев, Е.В. Карлова, П.А. Замыцкий // Саратовский научно-медицинский журнал. - 2017. - Т. 13, № 2. - С. 375-378.
15. Ильясова, Н.Ю. Оценивание геометрических признаков пространственной структуры кровеносных сосудов / Н.Ю. Ильясова // Компьютерная оптика. - 2014. - Т. 38, № 3. - С. 529-538. - DOI: 10.18287/0134-2452-2014-38-3529-538.
16. Хорин, П. А. Выделение информативных признаков на основе коэффициентов полиномов Цернике при различных патологиях роговицы человеческого глаза / П.А. Хорин, Н.Ю. Ильясова, Р.А. Парингер // Компьютерная оптика. - 2018. - Т. 42, № 1. - С. 159-166. - DOI: 10.18287/2412-6179-2018-42-1-159-166.
17. Широканев, А.С. Исследование алгоритмов расстановки коагулятов на изображение глазного дна / А.С. Широканев, Д.В. Кирш, Н.Ю. Ильясова, А.В. Куприянов // Компьютерная оптика. - 2018. -Т. 42, № 4. - С. 712-721. - DOI: 10.18287/2412-61792018-42-4-712-721.
18. Ильясова, Н.Ю. Диагностический комплекс анализа изображений сосудов глазного дна / Н.Ю. Ильясова // Биотехносфера. - 2014. - Т. 3(33). - С. 20-24.
19. Ильясова, Н.Ю. Методы цифрового анализа сосудистой системы человека. Обзор литературы / Н.Ю. Ильясова // Компьютерная оптика. - 2013. - Т. 37, № 4. - С. 511-535. - DOI: 10.18287/0134-2452-2013-37-4511-535.
20. Ильясова, Н.Ю. Измерение биомеханических характеристик сосудов для ранней диагностики сосудистой патологии глазного дна / Н.Ю. Ильясова, А.В. Устинов, А.В. Куприянов, М.А. Ананьин, Н.А. Гаврилова // Компьютерная оптика. - 2005. - № 27. - С. 165-169.
21. Сойфер, В. А. Методы компьютерного анализа диагностических изображений глазного дна / В.А. Сойфер, Н.Ю. Ильясова, А.В. Куприянов, А.Г. Храмов, М.А. Ананьин // Технология живых систем. - 2008. -Т. 5, № 5-6. - C. 61-71.
22. Симчера, В.М. Методы многомерного анализа статистических данных / В.М. Симчера. - М: Финансы и статистика, 2008. - 400 с. - ISBN: 978-5-279-03184-9.
23. Dementyiev, V.E. Use of images augmentation and implementation of doubly stochastic models for improving accuracy of recognition algorithms based on convolutional neural networks / V.E. Dementyiev, N.A. Andriyanov, K.K. Vasilyiev. - In: 2020 Systems of Signal Synchronization, Generating and Processing in Telecommunications (SYNCHROINFO). - 2020. - P. 1-4. - DOI: 10.1109/SYNCHROINFO49631.2020.9166000.
24. Vasiliev, K.K. Using probabilistic statistics to determine the parameters of doubly stochastic models based on autoregression with multiple roots / K.K. Vasiliev, V.E. Dementyiev, N.A. Andriyanov // Journal of Physics: Conference Series. - 2019. - Vol. 1368. - 032019. - DOI: 10.1088/1742-6596/1368/3/032019.
25. Jung, J.J. NAVILAS Laser System Focal Laser Treatment for Diabetic Macular Edema - One Year Results of a Case Series / Jung J.J., Gallego-Pinazo R., Lleo-Perez A., Huz J.I., Barbazetto I.A. // Open Ophthalmology Journal, 2013. - Vol. 7. - P. 48-53.
26. Ковеня, В.М. Алгоритмы расщепления при решении многомерных задач аэрогидродинамики / В.М. Ковеня; под ред. Ю.И. Шокина. - Новосибирск: Изд-во СО РАН, 2014. - 280 c.
27. Чернов, В.М. Параллельная машинная арифметика для рекуррентных систем счисления в неквадратичных полях / В.М. Чернов // Компьютерная оптика. - 2020. -Т. 44, № 2. - С. 274-281. - DOI: 10.18287/2412-6179-CO-666.
28. Якобовский, М.В. Введение в параллельные методы решения задач: Учебное пособие / М.В. Якобовский. - М.: Издательство Московского университета, 2012. - 328 с.
29. Fadeev, D.A. High performance 2D simulations for the problem of optical breakdown / D.A. Fadeev // Computer Optics. - 2016. - Vol. 40(5) - P. 654-658. - DOI: 10.18287/2412-6179-2016-40-5-654-658.
30. Поляков, М.В. Математическое моделирование пространственного распределения радиационного поля в биоткани: определение яркостной температуры для диагностики / М.В. Поляков, А.В. Хоперсков // Вестник Волгоградского государственного университета. - 2016.
- Т. 36, № 5. - С. 73-84.
31. Поляков, М.В. Численное моделирование динамики распространения температуры в биологической ткани / М.В. Поляков; под ред. Д.А. Новиковой, А.А. Ворониной.
- В кн.: Материалы всероссийской школы-конференции молодых ученых. - 2015. - Т. 1. - С. 971-978.
32. Пушкарева, А.Е. Методы математического моделирования в оптике биоткани. Учебное пособие / А.Е. Пушкарева. - СПб: СПбГУ ИТМО, 2008. - 103 с.
33. Kistenev, Y. Modeling of IR laser radiation propagation in bio-tissues / Y. Kistenev, A. Buligin, E. Sandykova, E. Sim, D. Vrazhnov // Proceedings of SPIE. - 2019. - Vol. 11208.
- 112081Q.
34. Самарский, А.А. Схемы повышенного порядка точности для многомерного уравнения теплопроводности / А.А. Самарский // Журнал вычислительной математики и математической физики. - 1963. - Т. 3(5). -С. 812-840.
35. Ануфриев, И.Е. Математические методы моделирования физических процессов. Метод конечных разностей. С решениями типовых задач: Учебное пособие / И.Е. Ануфриев, П.А. Осипов. - СПб: Издательство Санкт-Петербургского государственного политехнического университета, 2014. - 130 с.
36. Федоров, А.А. Сравнение двух методов распараллеливания прогонки на гибридных ЭВМ с графическими ускорителями / А.А. Федоров, А.Н. Быков // Вопросы атомной науки и техники. Серия: математическое моделирование физических процессов. - 2016. -№ 4. - С. 40-50.
37. Широканев, А. С. Разработка векторного алгоритма параметрической идентификации трёхмерных кристаллических решёток на основе оценки расстояний между двумерными слоями / А.С. Широканев, Д.В. Кирш, А.В. Куприянов // Информационные технологии и нанотехнологии (ИТНТ-2017). Сборник трудов по материалам III Международной конференции и молодежной школы. - 2017. - С. 1615-1619.
Сведения об авторах
Широканев Александр Сергеевич, 1993 года рождения, аспирант Самарского национального исследовательского университета имени академика С.П. Королёва; ассистент кафедры технической кибернетики, Самарский университет. Сфера научных интересов: интеллектуальный анализ медицинских изображений; цифровая обработка изображений; математическое моделирование; численные методы. E-mail: alexandrshirokanev@,gmail.com .
Андриянов Никита Андреевич, 1990 года рождения. В 2013 году окончил с отличием Ульяновский государственный технический университет (УлГТУ) по специальности «Инфокоммуникационные технологии и системы связи». В 2017 году защитил диссертацию на соискание учёной степени кандидата технических наук. В настоящее время работает доцентом департамента анализа данных и машинного обучения Финансового университета при Правительстве Российской Федерации. Область научных интересов: обработка изображений, интеллектуальный анализ данных, статистическая обработка сигналов. E-mail: [email protected] .
Ильясова Наталья Юрьевна, 1966 года рождения. В 1991 году окончила с отличием Самарский государственный аэрокосмический университет имени С.П. Королёва (СГАУ). В 1997 году защитила диссертацию на соискание степени кандидата технических наук, в 2015 году защитила диссертацию на соискание степени доктора технических наук. В настоящее время работает старшим научным сотрудником в Учреждении Российской академии наук Институте систем обработки изображений РАН и одновременно доцентом кафедры технической кибернетики СГАУ. Круг научных интересов включает цифровую обработку сигналов и изображений, анализ и интерпретацию биомедицинских изображений. Имеет более 100 публикаций, в том числе 35 статей и три монографии (в соавторстве). E-mail: [email protected]
ГРНТИ: 27.35.45
Поступила в редакцию 2 ноября 2020 г. Окончательный вариант - 9 марта 2021 г.
Development of vector algorithm using CUDA technology for three-dimensional retinal laser coagulation process modeling
A.S. Shirokanev12, NA.Andriyanov3, N.Y. Ilyasova1,2 1 Samara National Research University, 443086, Samara, Russia, Moskovskoye Shosse 34, 2IPSIRAS - Branch of the FSRC "Crystallography and Photonics" RAS,
443001, Samara, Russia, Molodogvardeyskaya 151, 3 Financial University under the Government of the Russian Federation, 125167, Moscow, Russia, Leningradskii prospekt 49
Abstract
For diabetic retinopathy treatment, laser coagulation is used in modern practice. During the laser surgery process, the parameters of laser exposure are selected manually by a doctor, which requires the doctor to have sufficient experience and knowledge to achieve a therapeutic effect. On the basis of mathematical modeling of the laser coagulation process, it is possible to estimate the crucial parameters without performing an operation. However, the retina has a rather complex structure, and when even low-cost numerical methods are used for modeling, it takes a long time to obtain a result. In this regard, the development of time-efficient algorithms for three-dimensional modeling is an urgent task, since the use of such algorithms will provide a comprehensive study within a limited time.
In this paper, we study the execution time of algorithms that implement various variations in the application of the splitting method and the finite difference method, adapted to the set problem of heat conduction. The study reveals the most efficient algorithm, which is then vectorized and implemented using the CUDA technology. The study was carried out using Intel Core i7-10875H and Nvidia RTX 2080 MAX Q and showed that an analog of the vector algorithm, focused on solving a multidimensional heat conduction problem, provides an acceleration of no more than 1.5 times compared to the sequential version. The developed vector-based algorithm, focused on the application of the sweep method in all directions of the three-dimensional problem, significantly reduces the time spent on copying into the memory of the video card and provides a 40-fold acceleration in comparison with the sequential three-dimensional modeling algorithm. On the basis of the same approach, a parallel algorithm of mathematical modeling was developed, which provided a 20-fold acceleration at full processor load.
Keywords: diabetic retinopathy, laser coagulation, mathematical modeling, heat conduction equation, parallel algorithms, vector algorithms, CUDA.
Citation: Shirokanev AS, Andriyanov NA, Ilyasova NY. Development of vector algorithm using CUDA technology for three-dimensional retinal laser coagulation process modeling. Computer Optics 2021; 45(3): 427-437. DOI: 10.18287/2412-6179-C0-828.
Acknowledgements: The reported study was funded by RFBR, project numbers 19-31-90160, 1929-01135, and by the Ministry of Science and Higher Education of the Russian Federation within the State assignment to the Samara University and FSRC "Crystallography and Photonics" RAS.
References
[1] Gafurov SD, Katakhonov ShM, Holmonov MM. Features of the use of lasers in medicine [In Russian]. Eur Sci J 2019; 3(45): 92-95.
[2] Kotsur TV, Izmailov AS. Comparative estimation of laser coagulation efficiency in macular and microphotocoagula-tion of high density in diabetic maculopathy treatment [In Russian]. Ophthalmology Journal 2016; 9(4): 43-45. DOI: 10.17816/0V9443-45.
[3] Zamytsky EA Laser treatment of diabetic macular edema [In Russian]. Postgraduate Bulletin of the Volga Region 2015; 1-2: 74-80.
[4] Kozak I, Luttrull J. Modern retinal laser therapy. Saudi J Ophthalmol 2014; 29(2): 137-146.
[5] Doga AV, Kachalina GF, Pedanova EK, Buryakov DA Modern diagnostic and treatment aspects of diabetic macular edema. Ophthalmology, Diabetes 2014; 4: 51-59.
[6] Whiting DR, Guariguata L, Weil C, Shaw J. IDF diadetes atlas: global estimates of the prevalence of diabetes for 2011 and 2030. Diabetes Res. Clin. Pract 2011; 94(3): 311-321.
[7] Bratko GV, Chernykh VV, Sazonova OV, Kovaleva MV, Sidorova EG, Shishko AP, Mirochnik Lyu. On the early diagnosis and frequency of occurrence of diabetic macular edema and group formation at risk of its development. Siberian Scientific Medical Journal 2015; 35(1): 33-36.
[8] Vorobieva IV, Merkushenkova DA. Diabetic retinopathy in type two diabetic patients: Epidemiology and modern view on the pathogenesis. Review [In Russian]. Ophthalmology in Russia 2012; 9(4): 18-21. DOI: 10.18008/18165095-2012-4-18-21.
[9] Amirov AN, Abdulaeva EA, Minkhuzina EL. diabetic macular edema. Epidemiology, pathogenesis, diagnosis, clinical features, treatment. Kazan Medical Journal 2015; 96(1): 70-74.
[10] Astakhov YuS, Shadrichev FE, Krasavina MI, Grigoiieva NN. Modern approaches to diabetic macular edema treatment [In Russian]. Ophthalmology Journal 2009; 2(4): 59-69.
[11] Iskhakova AG. The results of the clinical economic analysis of the treatment of patients with diabetic retinopathy with macular edema. Postgraduate Bulletin of the Volga Region 2014; 1: 96 - 98.
[12] Umanets NN, Rozanova ZA, Alzein M. Intravitreal ranibi-zumab injection in the treatment of cystoid diabetic macular edema [In Russian]. Oftalmol Zh 2013; 2: 56-60.
[13] Cohen SM, Shen JH, Smiddy WE. Laser energy and dye fluorescence transmission through blood in vitro. Amer J Ophthalmol 1995; 119(4): 452-457.
[14] Zamytsky EA, Zolotarev AV, Karlova EV, Zamytsky PA. Analysis of the coagulates intensity in laser treatment of diabetic macular edema in a Navilas robotic laser system [In Russian]. Saratov Journal of Medical Scientific Research 2017; 13(2): 375-378.
[15] Ilyasova NYu. Estimating the geometric features of a 3D vascular structure. Computer Optics 2014; 38(3): 529-538. DOI: 10.18287/0134-2452-2014-38-3-529-538.
[16] Khorin PA, Ilyasova NYu, Paringer RA. Informative feature selection based on the Zernike polynomial coefficients for various pathologies of the human eye cornea. Computer Optics 2018; 42(1): 159-166. DOI: 10.18287/2412-61792018-42-1-159-166.
[17] Shirokanev AS, Kirsh DV, Ilyasova NYu, Kupriyanov AV. Investigation of algorithms for coagulate arrangement in fundus images. J Computer Optics 2018; 42(4): 712721. DOI: 10.18287/2412-6179-2018-42-4-712-721.
[18] Ilyasova NYu. Diagnostic complex for the analysis of fundus vessels [In Russian]. Biotechnosphere 2014; 3(33): 20-24.
[19] Ilyasova NYu. Methods for digital analysis of the human vascular system. Literature review. Computer Optics 2013; 37(4): 511-535. DOI: 10.18287/0134-2452-2013-37-4-511-535.
[20] Ilyasova NYu, Ustinov AV, Kupriyanov AV, Ananin MA, Gavrilova NA Measurement of biomechanical characteristics of vessels for early diagnosis of vascular pathology of the fundus [In Russian]. Computer Optics, 2005; 27: 165-169.
[21] Soifer VA, Ilyasova NA, Kupriyanov AV, Khramov AG, Ananyin MA. Methods for computer diagnostics using eye's fundus images [In Russian]. Technologies of Living Systems 2008; 5(5-6): 61-71.
[22] Simchera VM. Methods of multivariate analysis of statistical data [In Russian]. Moscow: "Financy i Statistika" Publisher; 2008. ISBN: 978-5-279-03184-9.
[23] Dementyiev VE, Andriyanov NA, Vasilyiev KK. Use of images augmentation and implementation of doubly stochastic models for improving accuracy of recognition algorithms based on convolutional neural networks. In Book: Proceedings of 2020 Systems of Signal Synchronization, Generating and Processing in Telecommunications (SYN-CHROINFO) 2020: 1-4. DOI: 10.1109/SYNCHROINFO49631.2020.9166000.
[24] Vasiliev KK, Dementyiev VE, Andriyanov NA. Using probabilistic statistics to determine the parameters of doubly stochastic models based on autoregression with multiple roots. J Phys Conf Ser 2019; 1368: 032019. DOI: 10.1088/1742-6596/1368/3/032019.
[25] Jung JJ, Gallego-Pinazo R, Lleo-Perez A, Huz JI, Barba-zetto IA. NAVILAS Laser System Focal Laser Treatment for Diabetic Macular Edema - One Year Results of a Case Series. Open Ophthalmology Journal 2013; 7: 48-53.
[26] Kovenya VM. Splitting algorithms for solving multidimensional problems of aerohydrodynamics [In Russian]. Novosibirsk: Publishing house of the SB RAS; 2014.
[27] Chernov VM. Parallel machine arithmetic for recurrent number systems in non-quadratic fields. Computer Optics 2020; 44(2): 274-281. DOI: 10.18287/2412-6179-CO-666.
[28] Yakobovsky MV. Introduction to parallel methods for solving problems: Textbook [In Russian]. Moscow: Publishing house of Moscow University; 2012.
[29] Fadeev DA. High performance 2D simulations for the problem of optical breakdown. Computer Optics 2016; 40(5): 654-658. DOI: 10.18287/2412-6179-2016-40-5654-658.
[30] Polyakov MV, Khoperskov AV. Mathematical modeling of the spatial distribution of the radiation field in biological tissue: determination of the brightness temperature for diagnostics [In Russian]. Bulletin of the Volgograd State University 2016; 36(5): 73-84.
[31] Polyakov MV. Numerical modeling of the dynamics of temperature distribution in biological tissue [In Russian]. In Book: Proceedings of the All-Russian School-Conference of Young Scientists 2015; 1: 971-978.
[32] Pushkareva AE. Methods of mathematical modeling in biotissue optics. Textbook [In Russian]. Saint-Petersburg: "SPbGU ITMO" Publisher; 2008.
[33] Kistenev Y, Buligin A, Sandykova E, Sim E, Vrazhnov D. Modeling of IR laser radiation propagation in bio-tissues. Proc SPIE 2019; 11208: 112081Q.
[34] Samarsky AA. Schemes of the increased order of accuracy for the multidimensional heat conduction equation. J Com-put Math Math Phys 1963; 3(5): 1107-1146.
[35] Anufriev IE, Osipov PA. Mathematical methods for modeling physical processes. Finite Difference Method. With solutions of typical problems: a textbook [In Russian]. Saint-Peterburg: Saint-Petersburg State Polytechnic University Publisher; 2014.
[36] Fedorov AA, Bykov AN. Comparison of two methods for parallelizing a run on hybrid computers with graphics accelerators [In Russian]. Problems of Atomic Science and Technology. Series: Mathematical Modeling of Physical Processes 2016; 4: 40-50.
[37] Shirokanev AS, Kirsh DV, Kupriyanov AV. Development of a vector algorithm for parametric identification of three-dimensional crystal lattices based on estimating the distances between two-dimensional layers [In Russian]. Proc ITNT-2017 Conference 2017; 1: 1615-1619.
Authors' information
Aleksandr Sergeevich Shirokanev, (b. 1993), graduated (2017) with a master's degree in Applied Mathematics and Informatics. At present he is a postgraduate student of Samara University. At present he is a assistant of the Technical Cybernetics department at the Samara University. The area of interests includes digital image processing, mathematical modeling, numerical analysis and intellectual analysis of medical images. E-mail: [email protected] .
Nikita Andreevich Andriyanov (b. 1990) graduated with honor (2013) from Ulyanovsk State Technical University with a master's degree in Infocommunication Technologies and Communication Systems. He received his PhD (2017) in Technical Sciences. At present he is Associate Professor at Data Analysis and Machine Learning department in Financial University under the Government of the Russian Federation. The area of his interests includes image processing, data mining, statistical signal processing. E-mail: nikita-and-nov@mail. ru .
Nataly Yurievna Ilyasova (b. 1966), graduated with honors from S.P. Korolyov Samara State Aerospace University (SSAU) (1991). She received her PhD (1997) and DSc (2015) in Technical Sciences. At present, she is a senior researcher at the Image Processing Systems Institute of the Russian Academy of Sciences, and holding a part-time position of Associate Professor at SSAU's Technical Cybernetics sub-department. The area of interests includes digital signals and image processing, pattern recognition and artificial intelligence, biomedical imaging and analysis. She's list of publications contains more than 100 scientific papers, including 35 articles and 3 monographs published with coauthors. E-mail: [email protected] .
Received November 2, 2020. The final version - March 9, 2021.
