CC BY
22
Методы математического моделирования лазерного воздействия на глазное дно для оценки терапевтического эффекта
А. С. Широканев А. С. Кибиткина 1, Н.Ю. Ильясова 12, А.А. Дегтярев 1 1 Самарский национальный исследовательский университет имени академика С.П. Королёва, 443086, Россия, г. Самара, Московское шоссе, д. 34, 2 ИСОИ РАН - филиал ФНИЦ «Кристаллография и фотоника» РАН, 443001, Россия, г. Самара, ул. Молодогвардейская, д. 151
Аннотация
Лечение диабетической ретинопатии с использованием лазерной коагуляции сетчатки по зонам глазного дна, выбираемым врачом, обеспечивается на основе наведения лазера, после чего осуществляется обстрел, способствующий образованию коагулятов, которые препятствуют кровоизлияниям. Основная проблема при проведении лазерной коагуляции заключается в выборе зон обстрела так, чтобы охватить все проблемные области и подобрать параметры лазерного воздействия, чтобы излишне не повредить сетчатку. У каждого пациента индивидуальная структура глазного дна. Необходимо учитывать индивидуальное строение патологических и анатомических элементов, чтобы достичь требуемого терапевтического эффекта.
Терапевтический эффект достигается при образовании коагулятов во всех областях, в которых происходит кровоизлияние. Сосудистый слой должен нагреться до достаточной температуры, чтобы образовался коагулят, но при этом не должна быть повреждена сетчатка. Такой эффект можно прогнозировать при помощи математического моделирования лазерного воздействия.
В настоящей работе рассматриваются методы математического моделирования лазерного воздействия, основанные на уравнении теплопроводности. Проводится сравнение методов по вычислительной сложности и устойчивости решения. Выполняется анализ возможности прогнозирования терапевтического эффекта при помощи методов математического моделирования лазерного воздействия.
Ключевые слова: диабетическая ретинопатия, лазерная коагуляция, терапевтический эффект, математическое моделирование, уравнение теплопроводности, начальное и граничные условия.
Цитирование: Широканев, А. С. Методы математического моделирования лазерного воздействия на глазное дно для оценки терапевтического эффекта / А. С. Широканев, А. С. Кибиткина, Н.Ю. Ильясова, А. А. Дегтярев // Компьютерная оптика. - 2020. - Т. 44, № 5. - С. 809-820. - DOI: I0.18287/2412-6179-C0-760.
Citation: Shirokanev AS, Kibitkina AS, Ilyasova NY, Degtyarev AA. Methods of mathematical modeling of fundus laser exposure for therapeutic effect evaluation. Computer Optics 2020; 44(5): 809-820. DOI: 10.18287/2412-6179-C0-760.
Введение
Лазерное излучение применяется во многих сферах деятельности, но особое внимание лазерному излучению уделяют в медицине при лечении различных заболеваний [1]. Одним из современных инструментов оперирования при отёках и кровоизлияниях сосудов является лазерная коагуляция. Лазерная коагуляция очень часто применяется при лечении диабетической ретинопатии - заболевания глазного дна, зачастую возникающего вследствие сахарного диабета.
Сахарный диабет является одним из наиболее распространённых и опасных эндокринных заболеваний, которые поражают кровеносные сосуды сетчатки, потенциально приводя к развитию диабетической ретинопатии. Хотя все части сетчатки подвергаются поражению во время этого заболевания, именно изменения в её центральной части в виде диабетического
макулярного отёка приводят к самой быстрой и необратимой потере зрения [2]. Точная и ранняя диагностика, наряду с адекватным лечением, может предотвратить потерю зрения более чем в 50 % случаев [3 -5]. Лазерная коагуляция является стандартом для лечения диабетической ретинопатии, её эффективность была подтверждена в ходе крупного исследования (БТОЯЕ, 1987) [6].
Воздействие лазером с достаточной мощностью приводит к образованию коагулята, который препятствует кровоизлияниям сосудов. Для достижения нужного терапевтического эффекта необходимо подбирать мощность так, чтобы образовался коагулят, но при этом повреждение сетчатки было минимальным. На эффективность лечения также влияет способ расстановки коагулятов. Сетчатка должна минимально подвергаться лазерному воздействию, то есть количество выстрелов должно быть минимальным.
При лечении диабетической ретинопатии используется излучение с длиной волны 532 нм, соответствующей зеленому свету [7]. При этом во время лечения врач может задавать длительность импульса и мощность лазера. Способ нанесения микроожогов выбирается на основе информации о внутреннем строении глазного дна, в частности, о строении сетчатки. Однако современные методики лечения зачастую основываются на ручном способе нанесения микроожогов.
Наиболее современное оборудование МЛУ1ЬЛ8 позволяет в полуавтоматическом режиме производить выстрелы, однако перед лечением необходимо вручную сформировать предварительный план расположения коагулятов. Ручное формирование предварительного плана характеризуется рядом проблем: врач должен быть достаточно опытным, чтобы локализовать патологические и анатомические элементы на изображении глазного дна. Тщательный выбор плана коагулятов требует очень много времени, поэтому редко применяется на практике.
Предварительный план может быть сформирован автоматически при помощи цифровых методов с возможностью корректировки врачом любого из этапов обработки [8 - 11]. Большинство методов основаны на применении алгоритмов распознавания образов и статистического анализа данных [12 - 13]. Технология, ориентированная на решение данной задачи, позволит существенно улучшить эффективность проводимой лазерной коагуляции. Для формирования предварительного плана коагулятов необходимо знать как минимум 2 параметра: мощность лазера при обстреле заданной области и минимальное расстояние между коагулятами для предотвращения излишнего повреждения сетчатки [14].
Математическое моделирование лазерного воздействия позволяет прогнозировать распределение тепла по истечении времени моделирования. На основе решения задачи математического моделирования могут быть оценены необходимая мощность лазерного излучения, обеспечивающая необходимый терапевтический эффект, расстояние распространения тепла, температура, до которой нагревается каждый слой, в частности сетчатка, и минимальное расстояние между коагулятами [15 - 16]. Терапевтический эффект будет определяться температурой, до которой нагревается сетчатка, и расстоянием, до которого будет распространяться тепло. Если сетчатка не нагревается до критической температуры, но при этом сосудистый слой нагревается до нужной температуры, то терапевтический эффект достигается.
1. Общая постановка задачи математического моделирования лазерного воздействия на сетчатку глаза
Математическое моделирование распространения тепла должно, так или иначе, основываться на моде-
лировании световой энергии, которая преобразовывается в тепловую энергию [17 - 19]. Распространение тепла после лазерного воздействия формируется при преобразовании энергии электромагнитного поля в тепло [18 - 19]. Данное преобразование выводится через уравнение баланса энергии.
Интенсивность лазерного излучения описывается Гауссовой функцией (1).
7 (г е
(1)
где а - радиус пятна; Р - мощность лазера.
Распределение температуры в момент времени, когда лазерное воздействие прекратилось, определяется формулой (2). Импульс длится пренебрежимо малый промежуток времени, вследствие чего воздействие импульса считается мгновенным, а дифракция не учитывается [20 - 21].
ц/( х, у, г) =
-}р(х, у%
е о р/(г)А
С об
+ Тс
(2)
где Тс - температура, сформированная в результате предыдущих выстрелов; р=р (х,у,г) - функция поглощения среды; Соб=Соб (х,у,г) - функция объёмной теплоёмкости среды в зафиксированный момент времени.
В общем виде задачу математического моделирования лазерного воздействия можно сформулировать в виде (3).
Со
дТ
= йх\>(к ■ gгadxyz (Т))
Т1=0 =^(х, У,г ) Т1 г = То
(3)
где Т=Т(х, у, г, 0 - распределение температуры; к=к (х, у, г, Т) - функция температуропроводности среды; Соб=Соб (х, у, г, Т)- функция объёмной теплоёмкости среды, зависящая также от температуры; Г -граница; div - дивергенция векторного поля; gгadxyz -градиент функции по пространственным координатам; То - температура на границах (для первого выстрела температура соответствует температуре ткани).
Важно отметить, что область определения должна быть достаточно большой, чтобы температура в результате лазерного воздействия не достигала ни одной границы области, кроме той, через которую проходит лазер. Лазер проходит через одну единственную границу (несущая граница), которая не может быть с фиксированной температурой. Такая проблема решается одним из двух способов.
Первый способ заключается в том, что задается функция на несущей границе, которая зависит от времени, и вводится специфическая линейная замена, которая сводит задачу к фиксированным (или нулевым)
граничным условиям. Соответствующее уравнение теплопроводности становится неоднородным. Задача будет решаться путем решения однородного уравнения с исходными начальным и граничными условиями и неоднородного уравнения с нулевым начальным и нулевыми граничными условиями. То есть задача так или иначе сведётся к решению задачи (3).
Второй способ заключается в расширении области определения так, чтобы она разделилась на информативную и неинформативную области. В информативной области начальное распределение задается так же по формуле (2). В неинформативной области распределение симметричное. Граница, находящаяся в неинформативной области, имеет фиксированное значение температуры. Распределение тепла в неинформативной области не вызывает никакого интереса: результат моделирования визуализируется только в информативной области. Вследствие своей простоты в настоящей работе применялся второй способ моделирования.
Зависимость функции объёмной теплоёмкости и температуропроводности от температуры приводят к сильной нелинейности поставленной задачи. Однако изменение формы сетчатки можно спрогнозировать по величине нагрева слоев. Поэтому будем рассматривать аналогичную задачу (4), в которой отсутствует зависимость от температуры.
Следует отметить, что для задачи (4) достаточно использовать простую замену Т = Т + Т0, чтобы преобразовать представленные граничные условия в нулевые. Далее рассматривается задача с нулевыми граничными условиями. Результаты моделирования отображают, насколько ткань нагревается в результате лазерного воздействия. При переходе к исходной температуре достаточно прибавить нормальную температуру тканей, которая приблизительно равна 36,5 °С.
дТ
Со6 (х,у,г) — = Лу(к(х,у, г) ■ gradXyZ (Т)) д'
Т| = е-р( х'у 'г> г Р( х, у, г) I (г )А' + |'=0 с
С об (х, у, г)
(4)
Т1Г = Т0
Результат моделирования задачи (4) представляет собой распределение температуры в различные моменты времени. На слое сетчатки температура не должна достигать критической. По мнению врачей, критическая температура на сетчатке приблизительно равна 38 - 40 °С. На сосудистом слое коагулят образуется примерно при 39 °С.
Реализация всех описанных методов и проведение экспериментов проводились с использованием пакета МайаЪ. В пакете имеются необходимые встроенные функции для линейной алгебры, матричных вычислений и распараллеливания некоторых алгоритмов.
Данный пакет также содержит средства для математического моделирования методом конечных элементов, но только для простейших двумерных задач, в которых среда однородная и отсутствуют собственные числа в уравнении. К такой задаче свести поставленную сложную задачу невозможно. В связи с этим метод конечных элементов реализован вручную с использованием пакета МайаЪ.
2. Методы решения задачи математического моделирования лазерного воздействия
2.1. Применение метода разделения переменных для решения поставленной задачи
Решение задачи (4) аналитическими методами затруднено. Метод разделения переменных можно применять для задачи, в которой функции Соб и р можно факторизовать, то есть представить в виде Соб (х, у, г) = С1 (х)С2(у)С3(г) и р (х, у, г) = Р^х^Щг). В этом случае решение задачи (4) сводится к решению однородных уравнений вида (5) и поиску неизвестных коэффициентов.
у''+ к (х) у '+Хс( х) у = 0.
(5)
Если рассматривать фрагмент глаза с однородной средой, например, прямоугольный фрагмент стекловидного тела, то вместо функций Соб, к и р рассматриваются константы. Для такой задачи решение представляется в виде (6).
Т(х, у, г, ') = ^ Лт
'Ят,п,я (х, у, г) ;
(6)
. тл . пп . яп где (х, у, г) = бш—— дхэт—ду эт—дг ;
X У А
Лт
ХУА
Ц | У (х, у, г) ■ дт,п,я (х, у, г)dxdydг ;
0 0 0
т2 п2 я2 Xтт = л2(--1---1--)д ;
т,п,! \х2 у 2 А 2
д
= у/ Соб1 к ;
X, У, А - длины сторон параллелепипеда, ограничивающего область определения задачи.
Решение данной задачи может быть использовано как начальное условие для задачи в слоистой среде. Однако в слоистой среде факторизовать функции Соб=Соб (х,у,г) и р=р (х,у,г) нельзя, поскольку для описания сложной формы сетчатки используются сложные функции. На рис. 1 схематически показаны слоистая структура глазного дна и область лазерного воздействия.
2.2. Применение численных методов
Метод конечных разностей часто используется при решении сложных задач математического моделирования. Для решения задачи (4) наиболее подхо-
т,п,я=1
дящими являются следующие разновидности метода конечных разностей:
1) применение явной разностной схемы;
2) применение неявной разностной схемы;
3) применение метода расщепления для решения однородных задач.
Явная разностная схема является наипростейшей, но требует достаточно большого числа узлов, чтобы схема была устойчивой.
Граничные условия
Рис. 1. Структура сетчатки глаза
Неявная схема обычно является устойчивой при любых значениях шагов дискретизации. Однако вычислительная сложность алгоритма на каждой итерации достаточно высокая.
Схема расщепления позволяет рассматривать одномерные задачи вместо двумерных, благодаря чему можно использовать неявную схему, не решая полноценную систему линейных уравнений, а применяя метод прогонки, что существенно быстрее.
Даже для небольших объёмов данных метод конечных разностей в трёхмерном случае необходимо распараллеливать. Иначе алгоритм будет выполняться очень долго.
Метод конечных элементов представляет собой полуаналитический метод решения задачи математического моделирования [22]. Область определения разбивается на подмножества, в каждом из которых рассматривается непрерывный кусок аппроксимирующей функции итогового решения. Решение описывается суммой слагаемых по базисным функциям, то есть представляется в виде (7).
:=Х С" у»
(7)
где у» - известные базисные функции; С» - неизвестные коэффициенты.
В общем случае для заданного дифференциального оператора £[•], соответствующего задаче Ь[и] = 0, можно записать условие ортогональности в виде (8).
X С»у»
V/
= 0, I = 1, ж.
(8)
Обычно дифференциальный оператор является линейным относительно неизвестных коэффициентов, что позволяет записать условие (8) в виде (9).
N
X С»\Т [у»], V'/ = 0,
I = 1, N.
(9)
Условие, записанное в виде (9), представляет собой систему однородных линейных алгебраических уравнений. Скалярные произведения системы (9) представляют собой двойные интегралы, которым уделяется большое внимание, так как вычисление данных интегралов обладает большим временем выполнения, а вид базисных функций влияет только на качество вычисления упомянутых интегралов.
К системе (9) чаще сводятся задачи, которые предполагают поиск собственных чисел. В таком случае необходимо найти собственные числа, при которых определитель матрицы равен нулю. Таким образом, каждому собственному числу соответствует однородная система линейных уравнений, которая имеет нетривиальное решение.
Перед применением метода конечных элементов задача подвергается разделению переменных, поскольку функции, зависящие одновременно от пространственных координат и от времени, отсутствуют. Решение ищется в виде Т(х,у,2,() = и((^(ху2).
Задачу относительно V можно решать методом конечных элементов, то есть функция V представляется в виде (10).
v(X, у, 2) = X Сиуп (X, у, 2) ,
(10)
где у» (х,у,2) - базисные функции.
Базисные функции в трёхмерном случае удобно представляются при разбиении трёхмерного пространства на тетраэдры. Вид такой функции достаточно сложен. С учетом особенностей поставленной задачи базисная функция может быть факторизована. Однако применение метода конечных элементов в трехмерном пространстве требует чрезмерно высоких вычислительных ресурсов. Описать двумерные базисные функции и одномерные намного проще.
3. Постановка двумерной задачи математического моделирования лазерного воздействия
Рассматривая задачу в сечении слоистой структуры глазного дна, получаем двумерную задачу, в которой уравнение теплопроводности представляется в виде (11). Координата х соответствует смещению лазерного излучения, а координата 2 соответствует глубине.
со6 (х, 2 = А (к (х, 2 )дТ) + д (к (х, 2 (11)
ОТ дх дх 02 02
I
п=1
п=1
В вычислительном плане двумерная задача намного проще, чем трёхмерная. Для решения двумерной задачи можно использовать последовательные алгоритмы. Среда, соответствующая двумерной задаче, описывается снимком (проекцией) слоистой структуры глазного дна. Терапевтический эффект в этом случае можно оценивать на основе достижения нужной температуры сосудов, учитывая при этом температуру на сетчатке, чтобы температура на сетчатке не превышала критическую.
3.1. Применение метода разделения переменных для решения двумерной задачи
Фрагмент стекловидного тела (рис. 2) представляет собой однородную среду. Для данного фрагмента коэффициенты функции соб(х, г) и к(х, г) будут постоянными коэффициентами. Соответствующая двумерная задача при заданном начальном условии у(х, г) и нулевых граничных условиях имеет стандартное решение, которое записывается в виде (12).
Зрительный .
нерв Сетчатка
Рис. 2. Однородный фрагмент стекловидного тела
Т(х, г,') = ^ Лт,пе"
т, п=1 X А
. тл . пл Б1П-х Б1П— г,
X г
(12)
4 г Г ^ ч . тл . пл , , где Лтп = хz ] ] х, z)sln-хx Б1п -z~гdxdг;
0 0
т2 п2
X
= л2 (~~Г + —).
^2 А2'
3.2. Применение численных методов для решения двумерной задачи
Область определения задачи содержит множество малоинформативных зон. Наибольший интерес представляют зоны, близкие к центру лазерного воздействия, а также содержащие неоднородности в слоях глазного дна. Метод конечных элементов позволяет задать сложную структуру элементов в области определения, используя примитивные геометрические фигуры, например треугольники.
На рис. 3 представлена триангуляция Делоне, построенная на проекции слоистой структуры. Серым цветом выделена сетчатка, чёрным цветом - сосудистый слой, а белым - стекловидное тело. Треугольники плотнее скапливаются в области прохождения лазера, на слое сетчатки и сосудистом слое. В области стекловидного тела, где лазер не проходит, треугольники скапливаются более разреженно.
Для построения сетки сгенерируются точки с заданной плотностью вероятностей таким образом, чтобы точки наиболее плотно сосредотачивались на границах между слоями, на сетчатке, сосудистом слое и области лазерного воздействия. То есть вероятность попадания точки в границу между слоями самая высокая, в слой сетчатки или сосудистый слой, а также в зону лазерного воздействия - чуть меньше первой вероятности, в область стекловидного тела - низкая. По сгенерированным точкам строится триангуляция Делоне с использованием встроенных средств среды МаНаЪ.
Рис. 3. Триангуляция, построенная на слоистой структуре глазного дна
Каждому слою соответствует своё значение объёмной теплоёмкости, температуропроводности и коэффициента поглощения. На рис. 4 представлен пример функции объёмной теплоёмкости.
Температура, К*10
г, мм
х10-3 4,15 4,10 4,05 4,00 3,95 3,90 3,85 3,80 3,75 3,70
I
Рис. 4. Вид функции объёмной теплоёмкости
На рис. 5 представлен вид начального условия для двумерной задачи. При воздействии лазером на каждом новом слое происходит скачок температуры. В результате сосудистый слой нагревается больше, чем сетчатка и стекловидное тело.
Температура, К 50 '
45 40 35 30 25 20 15 10 5
- 4 0,2 ' х, мм
Рис. 5. Начальное распределение температуры на глазном дне после лазерного воздействия
Как и в трёхмерном случае, перед применением метода конечных элементов следует разделить переменные на временную координату и пространственные, то есть представить решение в виде ТХ х, 2)=и(()г(х, 2). Так образуются две задачи, которые записываются в виде (13) и (14).
dU ^)
Л
1
Соб
-Ш ^ ) = 0.
±Г к — 1+АГ к —
дх I дх I д21 с2
+ Хг = 0 .
(13)
(14)
Уравнение (13) представляет собой обычное уравнение первого порядка и имеет известное решение в виде U(t) = Ее'Задача (14) с нулевыми граничными условиями может быть решена методом конечных элементов. Решение в таком случае представляется в виде (15).
V(х, 2) = Х С»уп (х, 2).
(15)
Базисная функция описывается полиномом второго порядка, которая удовлетворяет условиям (16).
у(х, 2,- ) =
[1, ( х,, 2,- ) = р [0, (х,-, 2,- ) р'
(16)
где р - центральная точка базисной функции; х,, 2, -узлы триангуляции.
Коэффициенты полинома вычисляются на основе представления базисной функции в виде произведения двух двумерных кусочно-линейных функций У1(х, 2) и у2(х, 2).
Линейная функция представляется в виде (17).
у (х, 2) = |Акх + Вк2 + Ск , (х, 2) 6 Дк
(17)
где Дк - к-й треугольник.
Функция у1(х, 2) в узлах триангуляции соответствует условию (16). Функция у2(х, 2) удовлетворяет условию (18).
у(х,, 2, ) =
[1, ( х,, 2, ) = р [т, (х,-, 2, ) р
(18)
где т - регулирующий параметр; р , х,, 2, - то же, что и для (16).
Оптимальное значение регулирующего параметра равно 2. Легко показать, что функция у(х, 2) = у1(х, 2)у2(х, 2) удовлетворяет условию (16).
Однако у таких базисных функций имеется недостаток: в рёбрах треугольника функция не является гладкой, что может повлиять на итоговый результат. Для смягчения требования к гладкости дифференциальное уравнение обычно сводится к интегральному.
В соответствии с условием (9) записывается система (19).
где
X С» [А»/ + ХВ» ] = 0, / =1, N
п=1
А»/ =Я Т^
(19)
—[ к ду»
дх I дх
д_
д2
.ду» д2
у/dxd2;
В»/ = Ц уnУ/dxd2 .
Как было упомянуто ранее, требование гладкости может быть смягчено, если применить формулу Грина. В таком случае формулы для интегралов будут выглядеть сложнее, однако удастся беспрепятственно использовать базисные функции первого порядка.
Система (19) имеет нетривиальное решение, если определитель матрицы равен нулю. Поиск собственных чисел и неизвестных коэффициентов системы сводится к задаче поиска собственных чисел и собственных векторов. Для этого система (19) переписывается в виде (20).
( в-1 (-А)-ХЕ) / = 0.
(20)
С учетом записи (20) приходим к выводу, что числа X являются собственными числами матрицы -В-1А, а векторы, являющиеся решениями системы (20), -собственными векторами упомянутой матрицы.
Собственный вектор является базисным решением системы (19). Это значит, что результат умножения базисного решения на любое вещественное число тоже является решением системы. Таким образом, для найденного собственного числа Хк решение (15) представляется в виде (21).
г^к (х, 2) = р^/к, у (х, 2)
(21)
где у (х, 2) - вектор значений функций в точке (х, 2).
о
п
»=1
Неизвестные коэффициенты общего решения вычисляются на основе решения системы линейных алгебраических уравнений, формируемых при помощи начального условия. В результате поиска неизвестных коэффициентов итоговое решение представляет собой непрерывную функцию, которая в компактной записи описывается в виде (22).
Т (х, 2, ,) = е"х,Ру(х, 2 ), (22)
где Р - матрица неизвестных коэффициентов; у(х, 2 ) = (х, 2 (х, 2 ))Т ; е= (1'...е~Хк') .
Метод конечных элементов требует очень много времени на вычисление значений матриц системы (19). Численное интегрирование недопустимо, поскольку для достижения нужной точности требуется произвести чрезмерно много операций. Вместо этого все функции-множители, участвующие в интегрировании, интерполированы на треугольной области. Тогда интеграл от произведения сложных функций сводится к интегралу от степенных функций. Погрешность возникает только при интерполяции заданных в узлах функций. Даже такое упрощение не позволяет избежать проблемы высокой вычислительной сложности.
Метод конечных разностей, как простой метод, очень хорошо подходит для поставленной задачи. Единственный его недостаток: необходимость разбивать сетку равномерно. Даже если информативной областью является лишь небольшая зона, все равно необходимо всю область определения дробить мелкой сеткой.
Для явной схемы итерационный процесс в общем случае записывается в виде (23).
и*+1 = Щи* + + + Щи*+1 + . (23)
Ключевой недостаток явной схемы заключается в условии устойчивости. Схема устойчива, если шаг дискретизации по времени достаточно мал по сравнению с шагами дискретизации по пространственным координатам. В таком случае требуется задавать очень большое число узлов по времени.
Неявная схема в общем случае записывается в виде (24).
В\и\+1 + Щи*+1 + Щи*+1 + + = и*. (24)
Матрица, соответствующая системе линейных уравнений, в этом случае будет разреженной, но всё же не позволяющей решать её методом прогонки. Однако неявная схема является более устойчивой, чем явная. Поэтому моделирование может быть проведено не на всём промежутке наблюдения.
Система линейных уравнений, соответствующая (24), не может быть решена методом прогонки, так как рассматривается двумерный случай. Задача мо-
жет быть численно сведена к одномерной при использовании метода расщепления. Метод расщепления позволяет численно решать 2 одномерные задачи вместо одной двумерной. При этом в одномерной задаче при построении неявной разностной схемы матрица является трёхдиагональной, то есть решение может быть найдено методом прогонки.
При таком варианте решения система устойчива, однако появляется дополнительная погрешность из-за преобразования двумерной задачи к одномерным.
4. Экспериментальное исследование задачи математического моделирования лазерного воздействия
Как уже было описано ранее, аналитически поставленная задача решается во фрагменте стекловидного тела. Аналитическое решение может быть использовано для определения начального условия при любом смещении проекции стекловидного тела.
На рис. 6 представлено решение задачи в начальный момент времени. Вычисление решения производилось при использовании 250000 слагаемых. В начальный момент времени температура на конце фрагмента стекловидного тела нагревается на 10 - 15 К.
Рис. 6. Распределение температуры в стекловидном теле в начальный момент времени
Со временем среда остывает, а температура распространяется по боковым направлениям (рис. 7). Форма распределения температуры меняется, в начальной проекции фрагмента температура снижается быстрее. Температура очень быстро падает в стекловидном теле. По истечении 1 секунды температура падает и отличается от температуры стекловидного тела на 0,05 К.
Вообще говоря, на начальной и конечной проекциях рассматриваемого фрагмента граничные условия представляют собой убывающие функции, так как область лазерного воздействия пересекает такие границы. Задача с ненулевыми граничными условиями, описываемыми какими-либо функциями, может быть сведена к задаче с нулевыми граничными условиями. Уравнение теплопроводности новой задачи
будет, скорее всего, уже неоднородным. Для решения неоднородного уравнения решение представляется в виде суммы решения аналогичного однородного уравнения и неоднородного с нулевым начальным условием. Решение неоднородного уравнения обычно основывается на виде решения однородного. Таким образом, ключевая задача заключается в решении именно однородного уравнения с нулевыми граничными условиями.
Температура, К 0,06 " *
0,6 0 °>2 ' х,мм
Рис. 7. Распределение температуры в стекловидном теле после 1 с
Метод конечных элементов, позволяющий вычислять собственные числа, обычно направлен на решение однородного уравнения, чтобы задача поиска собственных чисел функций, соответствующих частным решениям задачи, сводилась к задаче поиска собственных чисел и собственных векторов матрицы.
Для решения двумерной задачи методом конечных элементов были применены базисные функции, сформированные на основе поверхностей второго порядка. При использовании большого количества треугольных элементов алгоритм долго (в рамках двух суток при использовании параллельного алгоритма) вычисляет интегралы, соответствующие условию ортогональности. На рис. 8 представлено распределение температуры в начальный момент времени после применения метода конечных элементов при генерации 10677 треугольников.
Моделирование проводилось при использовании различных форм базисных функций. Упомянутые ранее базисные функции второго порядка позволяют с достаточной точностью аппроксимировать начальное распределение, то есть функцию (2). Однако при вычислении результата моделирования и визуализации решения в начальный момент времени наблюдаются некоторые особенности в собственных числах и собственных функциях, играющих ключевую роль в формировании решения. Часть собственных чисел оказывается комплексными. Остальные значения дают высокие погрешности лишь в некоторых точках частных решений, которые приводят к высокой сред-
ней погрешности общего решения. Визуально заметно, что погрешность уже для начального распределения очень высокая. По истечении 0,025 с. среда почти полностью остывает (рис. 9).
Температура, К 80„ '
г, мм
0,4 °-6 '>•2 х, мм
0 0
Рис. 8. Распределение температуры в начальный момент времени при использовании 10677 конечных элементов
При использовании 4309 треугольников начальное распределение формируется с большей точностью (рис. 10). То есть нельзя утверждать, что увеличение количества треугольников увеличивает точность решения задачи, что говорит о неустойчивости рассмотренного способа моделирования.
Рис. 9. Распределение температуры по истечении 0,025 с при использовании 10677 конечных элементов
Анализ каждого описанного этапа метода конечных элементов показал, что и высокая вычислительная сложность, и неустойчивость возникают при вычислении интегралов, выражающих условия ортогональности (20). Неустойчивость вызывает в основном вид базисной функции.
Для достижения нужной точности использование метода конечных элементов при решении поставленной задачи крайне затруднительно. Основная проблема использования метода конечных элементов -чрезмерно высокая вычислительная сложность. Даже параллельные алгоритмы производят вычисления недопустимо долго.
Температура, К 80
0,6 X, мм
Рис. 10. Распределение температуры в начальный момент времени при использовании 4309 конечных элементов
Наиболее простым и подходящим методом решения поставленной задачи является метод конечных разностей. Как выяснилось, данный метод позволяет достичь нужной точности, при этом вычислительная сложность меньше, чем у метода конечных элементов.
В настоящей работе сравнению подвергался также метод расщепления, который заключается в том, что исходная двумерная задача итерационно решается путем моделирования распределения температуры в следующий момент времени для двух одномерных задач. Одномерные задачи решаются соответственно при помощи неявных разностных схем, в которых системы линейных уравнений решаются методом прогонки.
При сравнении времени работы алгоритмов для одинаковых размерностей сеток явная схема по умолчанию является самой быстрой. В табл. 1 представлены результаты времени работы алгоритмов при использовании метода конечных разностей. В отличие от метода конечных элементов, который может формировать решение более одного дня, метод конечных разностей позволяет формировать решение в рамках 10 минут.
Табл. 1. Результаты времени выполнения алгоритмов,
основанных на методе конечных разностей, при использовании различного количества узлов по времени
Алгоритм Время при различном количестве узлов по времени моделирования, с
1000 10000 100000
Явная схема 0,1 0,6 6,6
Неявная схема 211,5 1746,0 20607,5
Метод расщепления 2,3 22,5 239,2
Явная схема хоть и работает намного быстрее при одинаковом количестве узлов по времени, однако требует задания достаточно большого числа узлов по времени для достижения нужной точности. При количестве узлов 1 000 000 явная схема уже выполняется примерно за 80 секунд. Тем не менее, данный результат лучше простейшей неявной схемы. Однако
метод расщепления оказывается лучше, поскольку использование 10 000 узлов по времени - оптимальное значение для метода расщепления (табл. 1).
Благодаря устойчивости к размерности сетки неявная схема позволяет при меньших размерностях сетки достигать достаточной точности. Если для явной схемы количество узлов по времени должно быть достаточно большим, то неявная схема требует меньшее количество узлов.
Неявная схема в простейшем виде требует выделения большого объема памяти: при моделировании на натурном снимке требуется примерно 950 ГБ. Вследствие этого на персональном устройстве проведение такого эксперимента невозможно.
Вместо неявной схемы исследованию подвергался метод расщепления, который предполагает решение одномерных задач неявными схемами. Благодаря трехдиагональному виду матриц системы линейных алгебраических уравнений решаются методом прогонки, что позволяет существенно сократить вычислительные затраты и сэкономить память. В табл. 2 представлены результаты времени выполнения алгоритмов при использовании оптимальных размерностей сеток, экспериментально подобранных для каждого алгоритма и соответствующих качеству натурных снимков ОКТ.
Табл. 2. Результаты времени выполнения алгоритмов,
основанных на методе конечных разностей, при использовании натурных снимков ОКТ и оптимальных параметров алгоритмов
Алгоритм Время, мин
Явная схема 108,08
Метод расщепления 9,64
Использование явной схемы приводит к тому, что моделирование производится почти 2 часа, когда метод расщепления выполняется всего за 10 минут. В трехмерном случае даже маленькие размерности сеток будут приводить к большому времени выполнения. При решении трехмерных задач необходимо использовать параллельные алгоритмы.
В результате численного моделирования при помощи явной схемы температура сглаживается (рис. 11). С течением времени температура распределяется плавно (рис. 12) и сосредотачивается на сосудистом слое.
На рис. 13 представлен пример численного моделирования распространения температуры на сетчатке глаза. Температура в начальный момент в стекловидном теле высокая, но быстро падает. Сосудистый слой нагревается больше, чем сетчатка, поэтому при правильном выборе мощности сетчатка не повреждается, а в сосудистом слое происходит процесс «кипения», благодаря чему образуется коагулят.
В рассматриваемой задаче формы слоёв не изменяются при воздействии температуры. На практике изменения за короткий промежуток времени незначительны. Иногда визуально изменения наблюдаются
лишь по истечении нескольких дней. Нелинейная задача, учитывающая возможные изменения форм сло-ёв, лучше всего решается методами конечных разностей. Другие способы решения приводят к затруднениям, возникающим из-за нелинейности.
г, мм 0,5 О
Рис. 11. Распределение температуры в момент времени 0,01 с в результате моделирования с использованием явной разностной схемы
Температура, К 20
16 14 12 10 8 6 4 2 О
1,0
0-. 1,5
, мм
г, мм
0,5
О
О
Рис. 12. Распределение температуры в момент времени 0,2 с в результате моделирования с использованием явной разностной схемы
В клинических условиях врачи выявляют примерную зависимость параметров излучения от толщины сетчатки на основе формируемой картины, полученной после лазерного воздействия на глазное дно. Температура при этом не учитывается в связи с невозможностью её измерить. Такие исследования не позволяют учитывать особенности строения сетчатки и других слоев и результат предшествующего лазерного воздействия. Моделирование позволяет более точно анализировать состояние сетчатки и сосудистого слоя после лазерного воздействия. Для выявления новых зависимостей потребуются дополнительные клинические исследования, чтобы сопоставить результаты моделирования и результаты лечения пациентов. Такие зависимости позволят более эффективно проводить лазерную коагуляцию, поскольку врач бу-
дет точнее знать, как планировать лечение диабетической ретинопатии, чтобы минимально повреждать сетчатку и обеспечить максимальное уменьшение патологических областей.
Разница температур, Т—Тс 1,0 I--
I
4,0 3,5 3,0 2,5 2,0 1,5 1,0 0,5 О
1,0 1,5
Расстояние, мм
Рис. 13. Распространение тепла на сетчатке глаза в результате моделирования
Заключение
Рассмотренные методы математического моделирования лазерного воздействия направлены на оценку терапевтического эффекта при лечении диабетической ретинопатии. Благодаря указанным методам оцениваются такие параметры, как расстояние распространения тепла на сосудистом слое, температура, до которой нагревается сетчатка, и температура, до которой нагревается сосудистый слой.
Поставленную задачу аналитически решать затруднительно. Однако можно с достаточной точностью моделировать распространение тепла при помощи простейшей явной разностной схемы. Недостаток данной схемы заключается в том, что необходимо задавать очень малый шаг дискретизации по времени, чтобы схема была устойчивой, что может привести к длительным вычислениям. Однако количество узлов по пространственным координатам не очень большое в связи с невысоким разрешением снимков ОКТ, которые позволяют получить соответствующий аппарат в офтальмологической клинике. Неявная схема лучше подходит при большом количестве узлов по пространственным координатам. Вычислительная сложность решения СЛАУ на каждой итерации может быть компенсирована меньшим количеством итераций по времени по сравнению с явной схемой.
Благодарности
Исследование выполнено при финансовой поддержке РФФИ в рамках научных проектов № 19-3190160, № 19-29-01135 и Министерства науки и высшего образования Российской Федерации в рамках выполнения государственного задания Самарского университета и ФНИЦ «Кристаллография и фотоника» РАН.
Литература
1. Поляков, М.В. Математическое моделирование пространственного распределения радиационного поля в биоткани: определение яркостной температуры для диагностики / М.В. Поляков, А.В. Хоперсков // Вестник Волгоградского государственного университета. - 2016. - Т. 36, № 5. - С. 73-84.
2. Замыцкий, Е.А Анализ интенсивности коагулятов при лазерном лечении диабетического макулярного отека на роботизированной лазерной установке Navilas / Е.А. Замыцкий, А.В. Золотарев, Е.В. Карлова, П. А. Замыцкий // Саратовский научно-медицинский журнал. - 2017. - Т. 13, № 2. - С. 375-378.
3. Ильясова, Н.Ю. Оценивание геометрических признаков пространственной структуры кровеносных сосудов // Компьютерная оптика. - 2014. - Т. 38, № 3. - С. 529-538.
4. Хорин, П.А. Выделение информативных признаков на основе коэффициентов полиномов Цернике при различных патологиях роговицы человеческого глаза / П.А. Хорин, Н.Ю. Ильясова, Р.А. Парингер // Компьютерная оптика. - 2018. - Т. 42, № 1. - С. 159-166.
5. Астахов, Ю.С. Современные подходы к лечению диабетического макулярного отека / Ю.С. Астахов, Ф.Е. Шадричев, М.И. Красавина, Н.Н. Григорьева // Офтальмологические ведомости. - 2009. - Т. 4. - С. 59-69.
6. Kozak, I. Modern retinal laser therapy / I. Kozak, J. Luttrull // Saudi Journal of Ophthalmology. - 2014. - Vol. 29, Issue 2. - P. 137-146.
7. Поляков, М.В. Численное моделирование динамики распространения температуры в биологической ткани / М.В. Поляков; под ред. Д.А. Новиковой, А.А. Ворониной. - Материалы всероссийской школы-конференции молодых ученых. - 2015. - С. 971-978.
8. Ильясова, Н.Ю. Диагностический комплекс анализа изображений сосудов глазного дна / Н. Ю. Ильясова // Биотехносфера. - 2014. - Т. 3, № 33. - С. 20-24.
9. Ильясова, Н.Ю. Методы цифрового анализа сосудистой системы человека. Обзор литературы / Н. Ю. Ильясова // Компьютерная оптика. - 2013. - Т. 37, № 4. - С. 511-535.
10. Ильясова, Н.Ю. Биомеханические характеристики сосудов для цифрового анализа изображений глазного дна / Н. Ю. Ильясова, А. В. Куприянов, Н. А. Гаврилова, Г.А. Шилкин, Н.И. Ланевская // Биомеханика глаза. III семинар: сборник трудов. - 2002. - С. 18-30.
11. Сойфер, В. А. Методы компьютерного анализа диагностических изображений глазного дна / В.А. Сойфер, Н.Ю. Ильясова, А.В. Куприянов, А.Г. Храмов, М.А. Ананьин // Технология живых систем. - 2008. -Т. 5, № 5-6. - C. 61-71.
12. Симчера, В.М. Методы многомерного анализа статистических данных / В.М. Симчера. - М: Финансы и статистика, 2008. - 400 с. - ISBN: 978-5-279-03184-9.
13. Fukunaga, K. Introduction to statistical pattern recognition / K. Fukunaga. - New York, London: Academic Press, 1972. - 369 p.
14. Пушкарева, А.Е. Методы математического моделирования в оптике биоткани : Учебное пособие / А.Е. Пушкарева. - СПб: СПбГУ ИТМО, 2008. - 103 с.
15. Liu, G. Digital focusing of OCT images based on scalar diffraction theory and information entropy / G. Liu, Z. Zhi, R.K. Wang // Biomedical Optics Express. - 2012. - Vol. 3, Issue 11. - P. 2774-2783.
16. Jiang, H. Morphologic features of retina pigment epithelial around fluorescein leakage sites in acute central serous cho-rioretinopathy before and after laser coagulation / H. Jiang, Y. Quanyong, J. Xiaoyan, X. Guoxu // Chinese Journal of Ocular Fundus Diseases. - 2016. - Vol. 32, Issue 3. -P. 266-269.
17. 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.
18. Moës, N. Imposing Dirichlet boundary conditions in the extended finite element method / N. Moes, E. Bechet, M. Tourbier // International Journal for Numerical Methods in Engineering. - 2006. - Vol. 67, Issue 12. - P. 1641-1669.
19. Wolfram, S. The Mathematica book (3rd ed.) / S. Wolfram // Assembly Automation. - 1999. - Vol. 19, Issue 1. -P. 77-77.
20. Самарский, А.А. Схемы повышенного порядка точности для многомерного уравнения теплопроводности / А.А. Самарский // Журнал вычислительной математики и математической физики. - 1963. - Т. 3(5). - С. 812-840.
21. Математический анализ. Числовые и функциональные ряды / Е.М. Рудой. - Новосибирск: НГПУ, 2010. - 95 с.
22. Хватцев, А. А. Дифференциальные уравнения в частных производных / А. А. Хватцев. - Псков: Псковский государственный университет, 2016. - 80 с.
Сведения об авторах
Широканев Александр Сергеевич, 1993 года рождения, аспирант Самарского национального исследовательского университета имени академика С.П. Королёва; ассистент кафедры технической кибернетики, Самарский университет. Сфера научных интересов: интеллектуальный анализ медицинских изображений; цифровая обработка изображений; математическое моделирование; численные методы. E-mail: [email protected] .
Кибиткина Алена Сергеевна, 1997 года рождения. В 2014 году окончила Самарскую государственную академию им. Наяновой, учится в Самарском национальном исследовательском университете им. академика С.П. Королева, специальность «Прикладная математика и информатика». Область научных интересов: обработка графических изображений, программирование, математическое моделирование в физике. E-mail: anacy1222@gmail. com .
Ильясова Наталья Юрьевна, 1966 года рождения. В 1991 году окончила с отличием Самарский государственный аэрокосмический университет имени С.П. Королёва (ныне - Самарский Университет). В 1997 году
защитила диссертацию на соискание степени кандидата технических наук, в 2015 году защитила диссертацию на соискание степени доктора технических наук. В настоящее время работает старшим научным сотрудником в ИСОИ РАН - филиале ФНИЦ «Кристаллография и фотоника» РАН и одновременно доцентом кафедры Технической кибернетики Самарского университета. Круг научных интересов включает цифровую обработку сигналов и изображений, анализ и интерпретацию биомедицинских изображений. Имеет более 100 публикаций, в том числе 35 статей и три монографии (в соавторстве). E-mail: ilyasova@smr. ru .
Дегтярёв Александр Александрович, 1953 года рождения. В 1977 году с отличием окончил Куйбышевский авиационный институт, ныне - Самарский университет, по специальности «Прикладная математика». В 1985 году защитил диссертацию на соискание степени кандидата технических наук. В настоящее время работает доцентом кафедры технической кибернетики Самарского университета. Область научных интересов: численные методы и компьютерное моделирование в физике. E-mail: aadegt@gmail. com .
ГРНТИ: 28.17.19
Поступила в редакцию 25 июня 2020 г. Окончательный вариант - 8 июля 2020 г.
Methods of mathematical modeling of fundus laser exposure for therapeutic
effect evaluation
A.S. Shirokanev ',2, A.S. Kibitkina',2, N. Y. Ilyasova ',2, A.A. Degtyarev1 'Samara National Research University, 443086, Samara, Russia, Moskovskoye Shosse 34, 2IPSIRAS - Branch of the FSRC "Crystallography and Photonics" RAS, 443001, Samara, Russia, Molodogvardeyskaya 151
Abstract
When laser coagulation of eye retina is carried out, the laser beam is directed to target retinal areas selected by an ophthalmologist. The exposure to laser light produces a photocoagulate. When using laser coagulation, the main problem is selecting both the laser exposure areas that would cover all pathological zones and the laser exposure parameters to prevent retina damage. Any patient has an individual fundus structure. The individual structure of pathological and anatomical elements must be taken into account to achieve the desired therapeutic effect.
The formation of coagulates in all hemorrhage-affected areas results in the desired therapeutic effect. The vascular layer must be heated to a sufficient temperature to form a coagulate. Such an effect can be predicted using mathematical modeling of laser exposure.
In this paper, we consider methods of mathematical modeling of laser exposure based on the solution of a heat equation. The methods are compared in terms of their computational complexity and solution stability. An analysis of the possibility of predicting the therapeutic effect using methods of mathematical modeling of laser exposure is carried out.
Keywords: diabetic retinopathy, laser coagulation, therapeutic effect, mathematical modeling, heat equation, initial and boundary conditions.
Citation: Shirokanev AS, Kibitkina AS, Ilyasova NY, Degtyarev AA. Methods of mathematical modeling of fundus laser exposure for therapeutic effect evaluation. Computer Optics 2020; 44(5): 809-820. DOI: I0.18287/2412-6179-C0-760.
Acknowledgements: This work was funded by the Russian Foundation for Basic Research under RFBR projects ##19-31-90160, 19-29-01135 and the Ministry of Science and Higher Education of the Russian Federation within a State assignment to Samara University and the FSRC "Crystallography and Photonics" RAS.
References
[1] Polyakov MV, Hoperskov AV. Mathematical modeling of the spatial distribution of the radiation field in biological tissue: determination of brightness temperature for diagnosis [In Russian]. Bulletin of Volgograd State University 2016; 36(5): 73-84.
[2] Zamytskiy E, et al. Analysis of the coagulates intensity in laser treatment of diabetic macular edema in a Navilas robotic laser system. Saratov Journal of Medical Scientific Research 2017; 13(2): 375-378.
[3] Ilyasova N. Evaluation of geometric features of the spatial structure of blood vessels. Computer Optics 2014; 38(3): 529-538.
[4] Khorin P, Ilyasova N, Paringer R. Informative feature selection based on the Zernike polynomial coefficients for various pathologies of the human eye cornea. Computer Optics 2018; 42(1): 159-166.
[5] Astakhov YS, Shadrichev FE, Krasavina MI, Grigoryeva NN. Modern approaches to the treatment of a diabetic macular edema [In Russian]. Ophthalmologic Sheet 2009; 4: 59-69.
[6] Kozak I, Luttrull J. Modern retinal laser therapy. Saudi Journal of Ophthalmology 2014; 29(2): 137-146.
[7] Polyakov MV. Numerical modeling of the dynamics of temperature distribution in biological tissue [In Russian]. Materials of the All-Russian School-Conference of Young Scientists 2015: 971-978.
[8] Ilyasova NYu. Diagnostic computer complex for vascular fundus image analysis [In Russian]. Biotehnosfera 2014; 3(33): 20-24.
[9] Ilyasova NYu. Methods for digital analysis of human vascular system. Literature review. Computer Optics 2013; 37(4): 511-535.
[10] Ilyasova NYu, Kupriyanov AV, Gavrilova NA, Shilkin GA, Lanevskaya NI. Biomechanical characteristics of blood vessels for digital image analysis fundus [In Russian]. Biomehanika Glaza 2002; 18-30.
[11] Soifer VA, Ilyasova NYu, Kupriyanov AV, Khramov AG, Ananin MA. Methods for computer diagnostics using eye's fundus images [In Russian]. Technologies of the Living Systems 2008; 5(5-6): 61-71.
[12] Simchera VM. Methods of multivariate statistical analysis [In Russian]. Moscow: "Financy I Statistika" Publisher; 2008.
[13] Fukunaga K. Introduction to statistical pattern recognition. New York, London: Academic Press; 1972.
[14] Pushkaryova AE. Mathematical modeling methods in optics of biological tissue [In Russian]. Saint-Petersburg: "SPbGU ITMO" Publisher; 2008.
[15] Liu G. Digital focusing of OCT images based on scalar diffraction theory and information entropy. Biomed Opt Express 2012; 3(11): 2774-2783.
[16] Jiang H, Quanyong Y, Xiaoyan J, Guoxu X. Morphologic features of retina pigment epithelial around fluorescein leakage sites in acute central serous chorioretinopathy be-
fore and after laser coagulation. Chinese Journal of Ocular Fundus Diseases 2016; 32(3): 266-269.
[17] Kistenev Y, Buligin A, Sandykova E, Sim E, Vrazhnov D. Modeling of IR laser radiation propagation in bio-tissues. Proc SPIE 2019; 11208: 112081Q.
[18] Moes N, Bechet E, Tourbier M. Imposing Dirichlet boundary conditions in the extended finite element method. Int J Numer Methods Eng 2006; 67(12): 1641-1669.
[19] Wolfram S. The Mathematica book (3rd ed.). Assem Au-tom 1999; 19(1): 77-77.
[20] Samarskiy AA. High accuracy schemes for the multidimensional heat equation [In Russian]. Journal of Computational Mathematics and Mathematical Physics 1963; 3(5): 812-840.
[21] Rudoi EM. Mathematical analysis. Numerical and functional series. Novosibirsk: "Novosibirsk State Pedagogical University" Publisher; 2010.
[22] Khvatsev AA. Partial differential equations. Pskov: "Pskov State University" Publisher; 2016.
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] .
Alena Sergeevna Kibitkina (b. 1997) graduated from the Samara State Academy named after M.V. Nayanova (2014), studied at the Samara National Research University named after Academician S.P. Korolyov, majoring in Computer Science. Research interests are graphics computer processing, programming, mathematical modeling. E-mail: anacy1222@gmail. com .
Nataly Yurievna Ilyasova (b. 1966), graduated with honors from S.P. Korolyov Samara State Aerospace University (presently, Samara Univercity) (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 Samara Univercity 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] .
Aleksandr Aleksandrovich Degtyaryov (b.1953) graduated with honour (1977) from the S.P. Korolyov Kuibyshev Aviation Institute (presently, Samara Univercity), majoring in Applied Mathematics. He received his PhD in Technical Sciences (1985). At present time he is holding a position of Associate Professor at Samara Univercity Technical Cybernetics sub-department. Research interests include numerical methods and computer simulation in physics. E-mail: [email protected] .
Received June 25, 2020. The final version - July 8, 2020.
