BECTH. MOCK. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИВЕРН. 2024. .5» 3. С. 3 16 Lomonosov Computational Mathematics and Cybernetics Journal
УДК 519.67
M. В. Абакумов1
О РЕАЛИЗАЦИИ ГЕМОДИНАМИЧЕСКИХ РАСЧЕТОВ НА ПРОСТРАНСТВЕННЫХ ГРАФАХ
Рассматриваются способы подготовки данных для проведения гемодинамических расчетов в ква-зиодпомерпом приближении па сложных пространственных графах эластичных сосудов. Описывается метод автоматического построения объемных (3D) моделей графов сосудистой системы по данным их каркасных моделей. Обсуждаются основные подходы к проведению расчетов параллельно с визуализацией их результатов па объемных моделях.
Ключевые слова: гемодинамические расчеты. ЗО-модели сосудистой системы, визуализация расчетных данных па графах.
DOI: 10.55959/МSU/0137 0782 15 2024 47 3 3 16
1. Введение. Математические модели движения жидкости но системе эластичных трубок имеют широкую область применения и, в частности, описывают течение крови но сердечнососудистой системе. Существуют различные подходы к построению подобных моделей (ем., например, обзор в [1|). Один из них квазиодномерный, в рамках которого отдельные сосуды или совокупность функционально однородных мелких сосудов представляются прямолинейными ребрами графа. Вершинам графа ставятся в соответствие различные характеристики органов, мышечных тканей или узлов ветвления сосудов. Квазиодномерное приближение позволяет моделировать гемодинамические течения на достаточно разветвленных графах сосудов, в том числе приближенно представляющих замкнутую систему кровообращения человека [1 41.
Помимо широкого круга проблем, связанных с иол учением достоверных данных о физиологических параметрах сосудов, построением адекватных моделей функционирования органов и тканей, совершенствованием численных методов решения возникающих систем уравнений, нс менее существенны вопросы реализации гемодинамических расчетов на сложных пространственных графах. Здесь следует выделить следующие основные задачи. Первая состоит в подготовке данных таких графов, которые могут содержать тысячи ребер и вершин [4|. При этом необходимо задавать физиологически адекватную топологию графа, а также многочисленные расчетные параметры на всех его элементах. Второй задачей является визуализация сеточных функций на графах. Сложность состоит в том, что цветовая индикация значений функций на ребрах, представленных отрезками, в случае сложных трехмерных графов трудно воспринимается визуально. Поэтому необходимо использование более сложных 3D-моделей графов. И, наконец, третья задача состоит в обеспечении проведения расчетов параллельно с интерактивной визуализацией их результатов.
В данной работе описываются методы решения указанных задач, предложенные автором и реализованные им в программном комплексе CVSS (Cardio-Vascular Simulation System), который на протяжении многих лет планомерно разрабатывается [1 61 научным коллективом сотрудников МГУ имени М. В. Ломоносова. Излагаются детали инженерного подхода к “сборке” сложных пространственных графов из простых фрагментов. Приводится метод автоматического построения объемной (3D) модели пространственного графа но данным его каркасной модели. Обсуждаются основные аспекты использования механизма многопоточности [7( для визуализации сеточных функций на объемных моделях в реальном времени расчета.
1 Факультет ВАЖ А1ГУ, проф., д.ф.-м.п., e-mail: vmabkOcs.msu.ru
4
Абакумов М. В.
2. Подготовка каркасных моделей пространственных графов. Геометрия пространственного графа сердечно-сосудистой системы описывается каркасной моделью, содержащей трехмерные координаты вершин и список ребер отройков, соединяющих вершины и представляющих сосуды. Ребрам и вершинам приписаны параметры, необходимые для проведения расчетов, описывающие? свойства сосудов, узлов их ветвления, органов или тканей. Для каждого ребра хранятся значения сеточных функций давления, скорости и концентраций переносимых кровотоком веществ, что позволяет прерывать расчет, редактировать какие-либо параметры модели и продолжать расчет, стартуя с данных сохраненного временного слоя. Имеются также параметры, не используемые в ходе расчетов. В частности, каждой вершине и ребру сопоставлены целочисленный идентификатор, а также дата и время последнего изменения их параметров. Эти значения используются для групповых изменений данных на сходных сосудах и “сборки” модели из фрагментов. К данным модели в целом относится текстовое описание, дата и время ее последнего редактирования, а также два зарезервированных массива для хранения целочисленных и вещественных параметров, относящихся к модели как к единому целому, таких, как расчетное время.
Данные вершин и ребер хранятся в виде динамических однонаправленных списков, что позволяет ограничивать сложность графа лишь доступным объемом оперативной памяти, а также эффективно осуществлять операции вставки и удаления элементов списков. Используется автоматическая нумерация вершин и ребер. Ребра задаются нарой номеров соединяемых ими вершин (первым указывается меньший номер) и упорядочены но возрастанию первых номеров, а в случае их совпадения вторых. Наличие нескольких ребер, задаваемых одинаковой нарой, не допускается. Такая упорядоченность позволяет осуществлять автоматический контроль целостности данных модели. Так, например, при удалении вершины удаляются ребра, задаваемые нарой, содержащей номер удаленной вершины (каскадное удаление). Далее уменьшаются на единицу номера всех вершин, следующих за удаленной, и соответствующие коррективы вносятся в нары, определяющие ребра. Таким образом осуществляется динамическая нумерация, не требующая контроля со стороны пользователя.
В каждый момент времени активны две каркасные модели: основная и добавочная. Каждая модель может быть пустой, т.с. не содержать данных. Добавочная модель как единое целое может произвольно перемещаться в пространстве относительно основной модели, а также масштабироваться. В любой момент можно осуществить слияние двух моделей. При этом данные добавочной модели дополняют основную модель, после чего добавочная модель становится пустой. В ходе слияния отслеживаются совпадающие вершины и ребра. Вершины считаются совпадающими, если расстояние между ними не превосходит задаваемого значения. Ребра считаются совпадающими, если их концы являются совпадающими вершинами основной и добавочной моделей. Вершины или ребра, совпадающие в основной и добавочной моделях, могут иметь различные параметры. Тогда в обновленную основную модель включается набор параметров е более поздним временем последнего изменения.
Использование механизма слияния моделей позволяет реализовать параллельную разработку различных фрагментов модели несколькими пользователями в локальных системах координат е последующей стыковкой этих фрагментов в единое целое, возможность дублирования повторяющихся фрагментов нутом их присоединения к основной модели в различных ракурсах, масштабирование существующей модели (например, для учета роста пациента).
Опишем алгоритм слияния более подробно. На нервом этане список вершин основной модели дополняется вершинами добавочной модели. Если добавляемая вершина не совпадает е вершинами, уже содержащимися в списке, она записывается в конец списка. В противном случае параметры добавляемой вершины замещают параметры существующей, если существующая вершина редактировалась ранее добавляемой. По окончании этого процесса всем вершинам добавочной модели ставятся в соответствие их номера в расширенном списке вершин основной модели. В соответствии с этой нумерацией корректируются начальный и конечный номера ребер добавочной модели. При необходимости эти номера меняются местами, если начальный номер больше конечного. На втором этане список ребер основной модели дополняется скорректированными ребрами
О реализации гемодинамических расчетов па пространственных графах
5
добавочной модели. Напомним, что список ребер является упорядоченным но возрастанию начальных номеров, а при их совпадении но возрастанию конечных. Если добавляемое ребро нс совпадает с ребрами, уже содержащимися в списке, то оно записывается в соответствующее но порядку место списка. В противном случае параметры добавляемого ребра замещают параметры существующего, сели существующее ребро редактировалось ранее добавляемого.
Поясним применение механизма слияния на примере. Пусть в произвольной системе координат создан и сохранен фрагмент, представляющий систему сосудов правой руки. Очистим основную модель, т.е. сделаем со и устой, после чего считаем сохраненный фрагмент в качестве добавочной мод сын. Разместим добавочную модель в пространстве в нужном положении, соответствующем, например, лежащему в плоскости XY телу человека, где ось Y является осью симметрии, т.е. направлена вдоль позвоночника. Используя слияние, перенесем преобразованную добавочную модель в основную, после чего сохраним со. Далее считаем со же в качестве добавочной модели. В этот момент основная и добавочная модели идентичны. Подвергнем добавочную модель преобразованию масштабирования с коэффициентом (—1) относительно оси X. При этом система сосудов добавочной модели будет соответствовать системе сосудов левой руки и автоматически займет нужное положение в пространстве. После слияния получится модель, объединяющая сосуды правой и левой руки. Аналогично можно сформировать модель, объединяющую системы сосудов правой и левой ноги, после чего состыковать со с моделью, содержащей сосуды рук, и так далее.
Нс менее эффективен механизм фильтрации сосудов. Во-первых, фильтрация позволяет “разгрузить” изображение за счет вывода только тех сосудов (ребер графа), которые удовлетворяют заданному критерию отбора. Отбор осуществляется но тину сосуда (артериальный или венозный), идентификатору, а также их произвольным комбинациям. Во-вторых, отобранные ребра текущей основной модели можно сохранить в качестве самостоятельного фрагмента. Это обеспечивает обмен фрагментами между различными моделями, а также возможность их редактирования независимо друг от друга. Впоследствии фрагменты могут вставляться в произвольные модели, причем “устаревшие” параметры совпадающих ребер будут замещены их более актуальными значениями.
Отмстим, что описанные выше возможности конструирования пространственного графа сердечно-сосудистой системы из готовых фрагментов принципиально важны. Это обусловлено тем, что для последующих вычислений для каждой вершины и ребра должны быть заданы десятки различных параметров, а в моделях, хотя бы приближенно представляющих кровеносную систему человека, количество вершин и ребер исчисляются тысячами. Поэтому желательно избежать повторного ввода одинаковых параметров.
Для минимизации ввода данных также используется группировка схожих но параметрам ребер и вершин нутом задания одинаковых идентификаторов. Группировка позволяет осуществлять операции установки одинаковых параметров наборам вершин и ребер. Кроме? того, некоторые параметры всех ребер рассчитываются на основе? известных координат вершин, например, длины ре?бе?р и углы с вертикальной осью, необходимые для учета влияния гравитации на кровообращение.
На реализации учета влияния гравитации следует остановиться отдельно. В ходе расчетов используются углы между направленными (от конца с меньшим номером к концу с большим) отрезками, представляющими сосуды, и вектором ускорения свободного падения. Без ограничения общности считается, что сила тяжести всегда направлена в отрицательном направлении оси Z. Таким образом, указанные углы определяются исключительно трехмерными координатами вершин графа. Для учета изменения положения тела как единого целого достаточно осуществить ортогональное преобразование (поворот) координат всех вершин графа. Для изменения же положения отдельных частей тела (изменение нозы) необходимо осуществлять подобные преобразования е подмножествами вершин.
Для определения таких подмножеств вводятся дополнительные параметры вершин. Первый определяет принадлежность вершины к определенной части тела, например, кисти, предплечью или плечу правой руки. При этом принадлежность считается иерархической. Так, если верши-
6
Абакумов М. В.
на принадлежит кисти, то она автоматически принадлежит предплечью и плечу правой руки. Второй параметр для большинства вершин равен нулю. В противном случае вершина является центром вращения соответствующей части тела. Как правило, для каждой части тела центр вращения один. Однако могут иметься и несколько вершин с одинаковым значением (отличным от нуля) второго параметра. В этом случае задастся составной центр вращения, определяющий положение сустава, координаты которого вычисляются как среднее арифметическое координат указанных вершин. Для изменения положения руки в целом достаточно осуществить поворот всех вершин, принадлежащих плечу, а стало быть и предплечью, и кисти этой руки, относительно центра вращения плеча. Для того, чтобы учесть изгиб руки в локте, достаточно повернуть вершины, принадлежащие предплечью, а значит и кисти, относительно центра предплечья. И, наконец, поворот кисти осуществляется преобразованием принадлежащих ей вершин относительно ее центра вращения. Путем указанных преобразований можно смоделировать произвольное изменение положения тела, после чего пересчитать углы между ребрами и вертикальной осью. Важно отмстить, что указанные преобразования возможно осуществлять непосредственно в ходе проведения расчетов (между временными шагами), что позволяет отслеживать эффекты, связанные с изменением положения тела.
Каркасные модели, используемые для представления сосудистой системы, удобны тем, что для определения их пространственной топологии необходим ввод минимального количества данных. Однако их визуальное восприятие при отображении на мониторе или печати, особенно при большом количестве сосудов, весьма затруднительно. То же можно сказать и об отображении результатов расчетов, например, нутом цветовой индикации сеточных функций на отрезках, соединяющих вершины. Поэтому необходимо иметь возможность наряду с каркасными моделями отображать также объемные (3D) модели. Вариант решения этой задачи нутом использования какого-либо ЗБ-редактора для независимой разработки объемных моделей здесь представляется неприемлемым но ряду причин. Во-первых, создание сложных 3D-моделей даже в современных редакторах занимает существенное время и требует немалых навыков. Во-вторых, неизбежны сложности с установлением соответствия между ребрами каркасной модели и примитивами объемной модели. И, наконец, модифицировать 3D-модель в редакторе после каждого изменения каркасной модели крайне трудозатратно. По этим причинам в комплексе CVSS используется принципиально иной подход. Объемная модель при необходимости генерируется исключительно но данным текущей каркасной модели. Для этого нс требуется ввода дополнительных данных, да и вообще какого-либо участия пользователя.
3. Генерация объемных моделей. Основная идея построения объемной модели графа сосудистой системы состоит в том, что отрезки, представляющие сосуды в каркасной модели (см. и. 2), заменяются гладкими кривыми. Кривые, представляющие соседние сосуды, сопрягаются в общей вершине, т.е. имеют в этой точке общую касательную. В узлах ветвления, где сходятся более двух сосудов, сопрягаются пары кривых, представляющих сосуды с большими сечениями. Набор указанных кривых образует оглаженный граф. Сглаженные графы строятся гго отдельности для артериальных гг венозных подмножеств сосудов. Далее генерируются оболочки каждой кривой, которые представляются дискретными наборамгг узлов, четырехугольными примитивами, задающими ггх связность, гг нормалями к поверхности оболочки в каждом узле.
Заметим, что набор узлов гг примитивов для каждой оболочки должен быть, с одной стороны, достаточным для качественной визуализации, а с другой нс должен быть избыточным, учитывая, что количество сосудов может исчисляться тысячами. Чтобы удовлетворить этим противоречивым критериям, необходим анализ крггвггзньг ребер сглаженного графа.
Оболочки крггвьгх представляют собой аппроксимацию поверхности, .любое сечение которой плоскостью, ортогональной кривой в каждой ее точке, является окружностью с центром в этой точке. Оболочки должны удовлетворять следующим требованиям. Во-первых, средний вдоль кривой радиус ггх круговых сечений должен соответствовать сечению представляемого сосуда. Во-вторых, необходимо обеспечить равенство радиусов сечений оболочек сопрягаемых крггвьгх в вершинах графа. Пргг этом наборы узлов гг нормалей сопрягаемых гранггц оболочек также должны совпадать.
О реализации гемодинамических расчетов па пространственных графах
7
Начнем с задачи построения оболочки кладкой кривой из семейства кривых Безве [8,9]. Параметризация кривой Безье те-го порядка, определяемая опорными точками Pk = (xk,yk, z), к = 0,1, ...,те, имеет вид
П
r(t)= Cktk(1 - t)n-kPk, 0 ф t ф 1,
k=0
k
те!
к!(те — к)!
Здесь r(t) — радиус-вектор точки на кривой. Такая кривая начинается в точке Ро при t = 0, заканчивается в точке Pn при t = 1, касается в этих точках отрезков PoPi и Pn-iPn соответственно, а также целиком лежит в выпуклой оболочке, порождаемой точками Pk-
Для построения объемной модели графа сосудистой системы используются кривые 4-го порядка
r(t) = (1 — t)4Po + 4t(1 — t)3Pi + 6t2(1 — t)2P 2 + 4t3(1 — t)Ps + t4P4.
Приведем также выражения для производных вектор-функции r(t):
r'(t) = 4[ — (1 — t)3Po + (1 — t)2(1 — 4t)Pi + 3t(1 — t)(1 — 2t)P2 + t2(3 — 4t)P3 + t3P4], r"(t) = 12 [(1 — t)2Po + 2(1 — t)(2t — 1)Pi + (6t2 — 6t + 1)P2 + 2t(1 — 2t)P3 + t2P4].
Важной подзадачей является построение разбиения указанной кривой на секции точками Mj = r(Tj), удовлетворяющего двум условиям. Во-первых, плотность точек разбиения должна быть выше на участках с большей кривизной. Во-вторых, общее количество точек должно быть пропорционально интегральной кривизне. Кривизна в точке, определяемой значением параметра t
‘ = фД) х г"
K(t) |r'(t)|3
а величина r = 1/к называется радиусом кривизны, который совпадает с радиусом соприкасающейся окружности в данной точке кривой. Отметим, что для окружности радиуса R кривизна в каждой точке равняется 1/R, а для прямой линии — нулю.
С учетом того, что длина кривой определяется интегралом
1
L = J |r'(t)|dt, о
введем функцию распределения точек разбиения следующего вида:
а(т)= (ко + K(t))|r'(t)|dt, к0 > 0, 0 ф т ф 1.
о
При постоянной кривизне к(к) такая функция задает равномерное распределение точек разбиения но длине кривой. В противном случае предписывает сгущение точек на участках с большей кривизной. Параметр ко задает степень влияния кривизны на отклонение распределения точек от равномерного по длине кривой. С увеличением ко влияние кривизны понижается. Условие ко > 0 гарантирует корректность функции распределения для отрезков прямых линий.
Осталось рассчитать требуемое количество точек разбиения и значения параметра ту определяющие их положение на кривой. Для этого введем среднее значение кривизны к = а(1)/Ь — ко-Отмстим, что указанное выражение даст точные значения постоянной кривизны для окружности и прямой. Задавая в качестве параметров количество отрезков ломаной mc, необходимое для представления окружности, а также минимальное количество точек на кривой mmjn, определим требуемое количество точек разбиения m следующим образом:
m = max
Lk mc c 2n
mmin
Абакумов М. В.
В частности, если интегральная кривизна LK совпадает с таковой для окружности (Lk = 2п), то в качестве m получим значение mc, а для отрезка прямой (к = 0) — значение mmin.
Для нахождения приближенных значений длины кривой L и параметра Tj, определяющих положения точек Mj, введем на отрезке t € [0,1] равномерную сетку ti = i/N, i = 0,
Будем предполагать, что размерность сетки N существенно превышает значение т. Вычислим коэффициенты к = к(С), li = |Р(С)| и интегральные суммы
ai = ai—i +
(к0 + Ki—i )li-1 + (к0 + Ki)li
2N
i = 1, 2 ,...,N, ao = 0,
Li = Li— 1 +
Lo = 0.
li—1 + li 2 N
(Можно использовать квадратурные формулы и более высокого порядка точности.) Тогда L ln, a(1) ~ on- Определяя значение a = on/m, выберем точки Mj следующим образом:
Mj = r(Tj), Tj = ti, где i : ai_i < ja ф ai, j = 1, 2,...,m, Mo = Po.
При таком выборе точек разбиения их позиции приближенно соответствуют функции распределения а(т) и выполняется условие Mm = P4 . Для каждой точки Mj помимо координат будут использоваться значения параметра Tj и касательные век торы Cj = r/(Tj) /1 r/(ту )|.
Следующая подзадача состоит в построении аппроксимации оболочки кривой на основе найденных точек разбиения Mj и векторов Cj. В качестве параметров здесь выступают радиусы круговых сечений начала и конца оболочки Ro и Rm. Радиус кругового сечения, проходящего через точку Mj определим, исходя из значения Tj, как Rj = (1 — Tj)Ro + TjRm. Далее для каждого j равномерно распределим по окружности радиуса Rj, лежащей в плоскости, ортогональной вектору Cj, узлы границы кругового сечения Nji (в качестве локального индекса вновь используем символ i). Для этого определим единичные векторы aj и bj, ортогональные вектору Cj = (cx,cy, cz)j, например, следующим образом:
aj =
__b =____________— a
|ajI, b = |bj|, a
(0, 0, 1), Cz,j =0,
(cz, 0, cx)j , cz,j = 0,
b j = Cj X aj .
Тогда координаты узлов на окружности вычисляются, исходя из равенств
Nji = Mj + Rj nji, nji = cos jiaj + sin jibj, ji = 2ni/ms.
Здесь i = 0,1,...,ms — 1, где параметр ms задает количество узлов на окружности, которое одинаково для всех j и всех ребер графа. Для каждого узла Nji оболочки помимо координат используются значения параметра Tj для привязки к сеточной функции на ребре и единичный вектор нормали оболочки в этом узле nji для отображения объемной модели сосуда.
Для отображения объемной модели сосуда необходимо также задать связность узлов Nji, т.е. индексы узлов, образующих четырехугольные примитивы, из которых состоит аппроксимированная оболочка. Отметим, что указание в качестве таковых четверок вида
{ji, (j + 1)i, (j + 1)(i + 1), j(i + 1)}
может оказаться неприемлемым. Это связано с тем, что не представляется возможным указать такое правило выбора векторов aj, которое гарантирует близость aj и aj+i, если точки Mj и Mj+i на кривой r(t) близки. Поэтому необходим дополнительный анализ узлов оболочки соседних сечений.
Рассмотрим узлы Nji и N(j+i)i соседних сечений. Среди всех i = 0,1,..., ms — 1 найдем номер p для которого значение скалярного произведения (njo, nj+iy) максимально, т.е. найдем узел N(j+i)p, нормаль в котором наиболее близка к нормали в узле Njo- В качестве примитивов выбираются четверки
{ji, (j + 1)(i + p), (j + 1)(i + 1 + p),j(i + 1)} , i = 0,1,...,ms — 1.
О реализации гемодинамических расчетов па пространственных графах
9
Здесь предполагается, что используется циклическая нумерация по индексу i, т.е. после i = ms — 1 следует i = 0. Отметим, что на практике эффективней использовать сплошную нумерацию узлов N, тогда четырехугольные примитивы задаются не восемью, а четырьмя целыми числами.
Итак, для построения оболочек всех сосудов достаточно задать глобальные параметры mc, mmin, ms, а для каждого ребра — определить точки P&, к = 0,1,...,4, и граничные радиусы сечений Ro = R(Po) и Rm = R(P4). Здесь необходим анализ связности ребер графа артериальной либо венозной части кровеносной системы, а также учет заданных постоянных значений сечений соответствующих им сосудов.
Изначально для каждой кривой Безье в качестве Po и P4 устанавливаются начальная и конечная вершины ребра графа, соответствующего сосуду. Прочие опорные точки, определяющие кривую, расставляются равномерно на отрезке PoP4:
Pi = Po + 4 (P4 — Po), i = 1, 2, 3.
В качестве R(Po) и R(P4) (изначально эти значения одинаковы) устанавливаются заданные радиусы сечения сосуда. В дальнейшем изначально установленные значения могут корректироваться.
На следующем этане анализируются все вершины графа. В каждой вершине ищется пара входящих ребер одного тина (артериальные или венозные) с максимальными сечениями. Если установленные на ребрах сечения совпадают, то приоритет отдается ребру с большей длиной. Может оказаться, что такая пара не найдена, если в вершину входит только один сосуд или два сосуда артериальный и венозный. В этом случае никаких действий не требуется. Если пара найдена, то необходимо осуществить сопряжение соответствующей нары кривых Безье в данной вершине.
P''
P4
Пусть для определенности первое из двух найденных ребер заканчивается, а второе начинается в вершине сопряжения. Опорные точки кривых Безье этих ребер будем помечать одним или двумя штрихами соответственно. В рассматриваемом случае вершина сопряжения совпадает с точками P4 и P^. Изменим положение точек P3 и P'( (см. рис. 1) следующим образом:
P' = р' 1
P3 = P4 — 4
P' P' PoP4
р' р''
PoP4
р' P'' PoP4
'' ''
P1 = Po
1
+ 4
'' '' Po P4
p' p''
PoP4
P' P'' PoP4
Тогда точки P3 и P1' окажутся на одной прямой, параллельной вектору PoP)", которая будет общей касательной кривых Безье в вершине сопряжения. Помимо указанной коррекции координат опорных точек заменим значения R(P4) и R(Pq) максимумом этих чисел, что гарантирует совпадение радиусов сечений смежных оболочек в вершине сопряжения.
10
Абакумов М. В.
В результате аналогичных действий во всех вершинах графа будут скорректированы опорные то чки P 1 и/или Рз тех кривых, которые подвергались сопряжению. Прочие кривые останутся отрезками прямых. В любом случае не меняются установленные изначально координаты точек Ро, Р4 и Р2 • Первые две жестко связаны с начальной и конечной вершинами ребра, а положение последней можно рассматривать в качестве свободного параметра кривых. Отмстим, что в качестве базового семейства можно использовать и кривые Безье 3-го порядка, тогда свободного параметра не будет. Однако это нецелесообразно, поскольку за счет выбора, например, Р2 = (Ро + Р4)/2 можно обеспечить меньшее отклонение оболочки от исходного ребра графа.
Рис. 2. Оболочки кривых 4-го порядка (к0 = 0.001)
Поело коррекции опорных точек и граничных радиусов для всех ребер строятся кривые Безье, рассчитываются интегральные суммы, аппроксимирующие функции распределения точек Mj на кривых, и строятся оболочки. При этом касательные векторы к сопрягаемым кривым в общей вершине могут оказаться обратными друг другу. Это может привести к тому, что наборы узлов, представляющие граничные сечения соответствующих оболочек, нс совпадут. Однако такой ситуации легко избежать, если глобальный параметр ms выбрать четным.
Проиллюстрируем результат построения оболочек на примере тех же двух ребер, что приведены на рис. 1, в совокупности с третьим. Придерживаясь прежней нумерации, будем именовать ребра, начиная с правого верхнего угла но часовой стрелке (см. рис. 2), вторым, первым и третьим. В данном случае первое криволинейное ребро (соответствующая ребру кривая) сопрягается как со вторым, так и с третьим, т.е. для него корректируются обе опорные точки Р1 и Р3. Как видно на рис. 2, наибольшей оказывается кривизна первого ребра, а наименьшей второго. В соответствии с этим кривые разбиваются на 15 и 5 секций. В данном примере параметр ко = 0.001, т.е. функция распределения а(т) задает приоритет распределения точек разбиения Mj по кривизне. В результате наиболее подробно представляются участки изгибов оболочек. Заметим также, что постоянные сечения соответствующих сосудов установлены таким образом, что R3 < Ri < R2, где Ri,R2,R3 — радиусы этих сечений. После генерации оболочек формируется объемная модель составного сосуда, радиус сечения которого постепенно убывает от значения R3 до Ri.
На рис. 3 приведены оболочки тех же ребер, однако при их генерации параметр ко функции распределения а(т) был выбран равным 1000. При таком выборе ко функция распределения предписывает равномерное распределение узлов разбиения но длине кривой. Как видно на рис. 3, форма объемной модели в целом нс претерпела существенных изменений. Однако гладкость изгибов оболочек передастся гораздо хуже, чем на рис. 2. Важно отмстить, что представленные на этих двух рисунках объемные модели имеют одинаковую сложность, т.е. равное количество уз-
О реализации гемодинамических расчетов на пространственных графах
11
Рис. 3. Равномерное распределение узлов разбиения (ко = 1000)
лов и четырехугольных примитивов. Приведенные примеры подтверждают, что анализ кривизны позволяет существенно улучшить качество генерируемых объемных моделей без увеличения их сложности.
Рис. 4. Фрагмент сосудистой системы
Примитивы оболочек в комплексе CVSS отображаются средствами библиотеки OpenGL [11]. При этом большая визуальная гладкость изображаемых поверхностей достигается путем использования информации о нормалях пц в узлах оболочек, которые вычисляются на этапе их генерации. Для иллюстрации приведем изображение объемной модели фрагмента сосудистой системы (см. рис. 4). Более темными цветами представлены сосуды артериальной части, а более светлыми венозной.
12
Абакумов М. В.
Для цветовой индикации сеточных функций на объемных моделях используются значения параметров Tj, известные для каждого узла Nj оболочки. Это позволяет сопоставить узлам одного сечения оболочки точку на соответствующем расчетном отрезке и интерполированное значение сеточной функции в этой точке. Диапазону значений функции между со минимумом и максимумом но веем сосудам сопоставляется цвет из определенной палитры. Минимум и максимум рассчитываются отдельно для сосудов артериальной и венозной частей. Цветовые палитры для артериальных и венозных сосудов также различны. Это связано с тем, что характерные значения давления в различных частях кровеносной системы существенно (почти на порядок) отличаются.
В заключение данного пункта отмстим, что для реализации основных операций (скроллинг, масштабирование, изменение ракурса, печать и т.д.) над изображениями каркасных и объемных моделей пространственных графов в комплексе CVSS используется предложенный автором метод “сдвоенного окна” [12].
4. Реализация процесса вычислений. Особенностью комплекса CVSS является относительная программная независимость его интерфейсной части от вычислительной (солвера,), которые написаны на языках Object Pascal [13] и Fortran [14] соответственно. Это позволяет сочетать удобство среды Delphi [13] для программирования визуальных компонентов в ОС Windows и достоинства современных оптимизирующих компиляторов языка Fortran. При этом реализуется лишь односторонняя связь “интерфейсная часть ^ солвер”, что позволяет исключить привязку вычислительной части комплекса к конкретной операционной системе.
Применительно к Windows солвер CVSS компилируется в виде динамически загружаемой библиотеки DLL (Dynamic Link Library), которая содержит набор процедур, а также внутренние структуры данных, размещенные в common-блоках. Для организации управления еолвсром обеспечивается прямой доступ управляющей программы к этим данным, а также к процедурам и функциям верхнего уровня библиотеки. Для этого оказывается достаточным добавить небольшие интерфейсные программные модули как в управляющей программе, так и в еолвсрс, обеспечивающие необходимые связи.
На стороне солвера указывается, какие процедуры и функции являются внешними, способ их вызова и внешние имена. Кроме этого обеспечивается доступ к адресам common-блоков, которые размещаются в оперативной памяти на этане загрузки DLL. На стороне управляющей программы также создастся интерфейсный модуль, в котором описываются соответствующие структуры данных, процедуры и функции.
Далее остановимся на реализации процесса вычислений, который сводится к заполнению данными текущей каркасной модели пространственного графа (ем. н. 2) внутренних структур еол-вера и запуску цикла временных шагов (MakeStep). Несмотря на кажущуюся алгоритмическую простоту указанных действий, имеются определенные нюансы, которые необходимо учитывать. Их изложение проведем на примере классов (объектов) Delphi, однако оказанное ниже актуально и для других сред программирования Windows-приложений.
Когда в качестве реакции на сообщение (событие) в приложении Windows запускается какая-либо процедура, другие сообщения нс обрабатываются вплоть до ее завершения. Если эта процедура продолжительна, то приложение нс реагирует на активность пользователя достаточно длительное время, нс обновляются визуальные компоненты приложения и т.н. (эффект зависания). В Delphi стандартное решение этой проблемы заключается в периодическом вызове в ходе выполнения процедуры метода ProcessMessages класса Application. При вызове этого метода процедура приостанавливается и обрабатываются все сообщения, накопившиеся в очереди сообщений приложения. Однако такое решение неприемлемо, если речь идет о запуске внешней процедуры солвера, которая нс имеет обратной связи с управляющей программой. Важно отмстить, что основная процедура солвера MakeStep может выполняться достаточно длительное время, поскольку реализует итерационный процесс, на каждом шаге которого обращается матрица большой размерности.
Для преодоления указанных трудностей используется механизм многопоточности [7]. Цикл временных шагов солвера запускается в виде отдельного дополнительного потока, в то время как
О реализации гемодинамических расчетов па пространственных графах
13
управляющая цршрамма функционирует в рамках основного потока. Это означает, что в ходе вв1 числений управляющая программа продолжает функционировать фактически независимо, т.е. вес управляющие элементы сохраняют активности. На этане синхронизации, которая необходима после каждого временного шага, дополнительный ноток приостанавливается до тех нор, пока основной ноток нс перейдет в режим ожидания. Этот режим наступает после того, как основным потоком обработаны вес ранее поступившие сообщения, вызванные командами пользователя или дополнительного потока на предыдущем временном шаге. Поело перехода основного потока в режим ожидания дополнительный ноток передаст основному сообщение о завершении очередного шага, после чего, нс дожидаясь его обработки, переходит к вычислениям на следующем шаге. Параллельно в основном потоке обрабатывается указанное сообщение, в качестве реакции на которое может выступать, например, обновление графиков сеточных функций на выбранных сосудах совместно с изменением цветовой индикации давления на ЗБ-модели графа сосудистой системы (см. н.З).
Необходимо отмстить, что но-настоящсму параллельно указанные действия будут осуществляться лишь на компьютерах с многоядерными процессорами, однако и в случае одного ядра указанный механизм взаимодействия потоков позволяет обеспечить интерактивность. На многоядерных системах временные затраты на один шаг' вычислений из-за необходимости синхронизации будут определяться продолжительностью более медленного из двух процессов временного шага солвсра или обработки его результатов. Если же время, затрачиваемое на отображение результатов, сравнительно невелико, то процесс вычислений будет вестись практически безостановочно.
Остановимся на особенностях реализации взаимодействия потоков в Delphi. Следует обратить внимание на то, что объектная модель визуальных компонентов Delphi, которая обеспечивает удобство программирования иод Windows, нс рассчитана на многопоточность. Поэтому параллельная модификация свойств визуальных компонентов несколькими потоками может привести к непредсказуемым последствиям. Соответствующие ошибки в приложениях крайне сложно обнаружить, поскольку они могут проявляться случайным образом в зависимости от конкретной последовательности модификаций свойств асинхронными потоками. Наиболее надежный и эффективный способ избежать подобных проблем состоит в использовании специально предназначенного для реализации дополнительных потоков в Delphi класса TThread.
Нс вдаваясь в излишние детали (см. [13]), опишем основные моменты, которые необходимо учитывать при использовании этого класса применительно к рассматриваемой задаче. В конструкторе дополнительного потока (наследника TThread), обеспечивающего процесс вычислений, указывается, что ноток необходимо запустить сразу же после создания (вис зависимости от реально переданного CreateSuspended False):
constructor TCVSSCalcThread.Create
(CreateSuspended: Boolean); begin
inherited Create(False);
FreeOnTerminate:=True;
Priority:=tpldle; end;
Далее указывается, что ноток автоматически высвобождает занятые им ресурсы после своего завершения. Приоритет потока целесообразно (см. далее) указать минимальнымДрЫ1е).
Для реализации цикла временных шагов солвсра в дополнительном потоке достаточно описать метод Execute:
procedure TCVSSCalcThread.Execute; begin
PErrors~.iErrorCode:=0;
while not Terminated do begin
14
Абакумов М. В.
MakeStep;
Synchronize(UpDateStep); if PErrors~. iErrorCodeoO then Terminate; end; end;
Предполагается, что перед запуском дополнительного потока осуществлены вес необходимые операции инициализации данных еолвера, необходимые для начала либо продолжения расчетов. Перед началом цикла обнуляется внутренняя переменная еолвера (iErrorCode), используемая для возврата кода возможных ошибок. Цикл выполняется до тех пор, пока свойство Terminated нс примет значение True. В теле цикла выполняется один шаг по времени еолвера MakeStep и специальным образом запускается метод обработки его результатов UpDateStep. Если при выполнении временного шага в солвсрс возникла ошибка, вызывается метод Terminate, который фактически устанавливает свойство Terminated равным True. В данном случае это приведет к завершению временного цикла и дополнительного потока (с последующим автоматическим освобождением занятых им ресурсов). Если же метод Terminate вызывается в ходе вычислений из основного потока, например, в результате нажатия управляющей кнопки, то осуществляется “мягкое” завершение дополнительного потока поело окончания очередного шага по времени и обработки его результатов.
Ключевым в реализации дополнительного потока является вызов обработчика временного шага UpDateStep с помощью метода Synchronize класса TThread. Этот метод обеспечивает запуск обработчика в основном потоке и остановку дополнительного потока в точке вызова Synchronize вплоть до завершения кода UpDateStep. Таким образом, в обработчике целесообразно осуществить только тс операции, которые нет возможности выполнять параллельно с вычислениями в процедуре MakeStep:
procedure TCVSSCalcThread.UpDateStep; begin
Application.ProcessMessages; if PErrors~ . iErrorCodeoO then begin CVSSCalculation.SaveStep;
PostMessage(FormClcCtrl.Handle,wm_UpDateStep,0,0); end; end;
Первым в обработчике вызывается метод ProcessMessages класса Application, который уже обсуждался ранее. Это гарантирует, что обработка всех сообщений, поступивших после предыдущих) шага и в ходе выполнения текущих), будет завершена до выполнения следующих действий. Если при выполнении текущих) временного шага возникла ошибка еолвера, никаких дальнейших действий нс требуется. В противном случае вызывается метод SaveStep класса TCVSSCalculation. Этот класс осуществляет обмен данными между управляющей программой и еолвером. Метод SaveStep обеспечивает сохранение сеточных функций, рассчитанных на текущем шаге для всех сосудов, во внутренние структуры данных объекта CVSSCalculation. В дальнейшем именно эти данные используются в процессе визуализации, поскольку они нс подвержены модификациям при вычислениях в солвсрс. Последним действием обработчика является направление сообщения форме, управляющей расчетом, о необходимости отреагировать на полное завершение текущего временного шага. При этом используется функция PostMessage, которая в отличие от SendMessage нс обращается напрямую к оконной функции формы, а лишь заносит сообщение в очередь. В результате выполнение обработчика завершается без ожидания реакции на сообщение и вычисления продолжаются незамедлительно. Обработка сообщения, которая, как правило, заключается в обновлении графиков сеточных функций на сосудах и/или цветовой индикации давления на объемной модели сосудистой системы, осуществляется параллельно с выполнением следующего временного шага еолвера.
О реализации гемодинамических расчетов па пространственных графах
15
На многоядерных процессорах указанный механизм обеспечивает двукратный прирост скорости работы программы, сели процессорное время, необходимое на шаг еолвсра и его визуализацию, приблизительно одинаково. Если это нс так, то, как уже отмечалось, скорость выполнения одного шага определяется более медленным из двух указанных процессов. По этой причине увеличение приоритета дополнительного потока нс приведет к ускорению работы программы. На одноядерных же системах повышение приоритета дополнительного потока приводит к увеличению времени отклика управляющих элементов программы, поскольку обработка соответствующих событий осуществляется в основном потоке. Поэтому целесообразно устанавливать дополнительному потоку минимальный приоритет tpldle.
5. Заключение. Таким образом, в настоящей работе описана методика подготовки сложных каркасных моделей пространственных графов эластичных сосудов, генерации их объемных моделей и реализации процесса гемодинамических расчетов на таких графах параллельно с визуализацией расчетных данных на объемных моделях. Эта методика реализована в программном комплексе CVSS и показала свою эффективность в ходе численного решения многих прикладных задач. Болес детальное описание функционала комплекса и прочих аспектов его разработки может послужить темой последующих публикаций.
СПИСОК ЛИТЕРАТУРЫ
1. Кошелев В.Б., Мухин С.И., Соснин Н.В. Фаворский А.И. Математические модели квазиодиомериой гемодинамики: Методическое пособие. М.: МАКС Пресс, 2010.
2. Борзов А.Г., М у х и и С.И., Соси и и Н.В. Консервативные схемы переноса вещества по системе сосудов, замкнутых через сердце // Дифференц. уравн. 2012. 48. № 7. С. 935 944. (Borzov A.G., М и k h i n S.I., S о s и i и N.V. Conservative schemes of matter transport in a system of vessels closed by the heart // Differential Equations. 2012. 48. N 7. P. 919 928.)
3. Б у и и ч e в а А.Я., M у x и и С.И., Соси и и H.B., X p у л е и к о А.Б. Математическое моделирование квазиодиомериой гемодинамики // ЖВМиМФ. 2015. 55. № 8. С. 1417 1428. (В unicheva A.Y., М и k h i n S.I., S о s n i n N.V., К h rulenk о A.В. Mathematical modeling of qnasi-one-dimensional hemodynamics // Computational Mathematics and Mathematical Physics. 2015. 55. N 8. P. 1381 1392.)
4. Абакумов M.B., Мухин С.И., Мысов а К.М., Покладюк А.Ю., X рул ей ко А.Б. Моделирование динамики веществ в рамках самосогласованной модели большого круга кровообращения // Математическое моделирование. 2023. 35. № 11. С. 3 20.
5. Абак у м о в М.В., Г а в р и л ю к К.В., Е с и к о в а Н.Б. и др. Математическая модель гемодинамики сердечно-сосудистой системы // Дифференц. уравн. 1997. 33. № 7. С. 892 898. (Abakumov M.V., Gavrilyuk K.V., Esikova Х.В. et al. Mathematical model of the hemodynamics of the cardio-vascular system // Differential Equations. 1997. 33. X 7. P. 895 901.)
6. А ш м e т к о в И.В., Б у и и ч е в а А.Я., Л у к ш и и В. А. и др. Математическое моделирование кровообращения на основе программного комплекса CVSS // Компьютерные модели и прогресс медицины. М.: Наука, 2001. С. 194 218.
7. Ncmirovsky М., Tullson D.M. Multithreading architecture // Synthesis Lectures on Computer Architecture. 2013. 8. XI. P. 1 109.
8. Шик и и E.B., Б о p e с к о в А.В. Компьютерная графика. Динамика, реалистические изображения. М.: Диалог МИФИ, 1996.
9. Роджерс Д., Адамс Д. Математические основы машинной графики. М.: Мир, 2001. (R ogers D.F., Adams J.A. Mathematical Elements For Computer Graphics. 2nd ed. X. Y.: McGraw-Hill, 1990.)
10. И о г о p e л о в А.В. Дифференциальная геометрия. М.: Наука, 1974.
11. В у М., Девис Т., Н е й д е р Д., Ш р а й и е р Д. OpenGL. Руководство по программированию. 4-е изд. СПб.: Питер, 2006. (Shreiner D., Woo М., Xeider J., Davis T. OpenGL Programming Guide. 4th ed. Boston: Addison-Wesley, 2004.)
12. Абак у м о в М.В. О некоторых методах визуализации сеточных данных. Препринт ИПМ РАН. № 72. 2004. С'. 1 40.
13. Т е й к с с й р а С., Пачек о К. Delphi 5. Руководство разработчика. М.: Вильямс, 2000. (Pacheco X., Teixeira S. Delphi 5. Developer’s Guide. Indianapolis: SAMS, 1999.)
16
Абакумов М. В.
14. Б а р т о н ь о в О.В. Современный Фортран. М.: Диалог МИФИ. 2005.
Поступила в редакцию 24.01.24 Одобрена после рецензирования 26.02.24 Принята к публикации 26.02.24