CC BY
37
УДК 004.932:004.925
DOI 10.25205/1818-7900-2019-17-3-5-17
Восстановление трехмерной геометрии сосудов по данным компьютерной томографии
В. В. Борисенко \ Н. С. Серова 2, А. М. Чеповский 3
1 Московский государственный университет им. М. В. Ломоносова Москва, Россия
2 Первый МГМУ им. И. М. Сеченова (Сеченовский Университет) Москва, Россия
3 Национальный исследовательский университет «Высшая школа экономики»
Москва, Россия
Аннотация
Рассматриваются алгоритмы построения триангуляции внутренней поверхности коронарных артерий сердца по данным компьютерной томографии. Точное определение геометрии сосудов необходимо для построения гидродинамической модели кровоснабжения сердца и расчета параметров кровотока с помощью уравнений Навье - Стокса. Для восстановления трехмерной геометрии применяется комбинация двух основных методов: трехмерного алгоритма роста области из семени и ячеечного метода, использующего разбиение пространства на тетраэдры. При этом используется тетраэдрическая сеть, предложенная в работах S. Chan, E. Purisima (1998) и V. Skala (2000), в которой тетраэдры строятся в кубической решетке на общих гранях смежных кубов. Эта сеть строится не на всем пространстве, а только в окрестности границы множества вокселей, вычисленного на первом этапе как результат применения трехмерного алгоритма роста области. Ключевые слова
компьютерная томография, триангуляция поверхности, марширующие тетраэдры, алгоритм роста области из семени, сглаживание поверхности Благодарности
Работа выполнена при поддержке РФФИ, грант № 18-29-26012 мк Для цитирования
Борисенко В. В., Серова Н. С., Чеповский А. М. Восстановление трехмерной геометрии сосудов по данным компьютерной томографии // Вестник НГУ. Серия: Информационные технологии. 2019. Т. 17, № 3. С. 5-17. DOI 10.25205/1818-7900-2019-17-3-5-17
Reconstruction of Three-Dimensional Geometry of the Vessels by Computed Tomography Data
V. V. Borisenko , N. S. Serova 2, A. M. Chepovskiy 3
1 Lomonosov Moscow State University Moscow, Russian Federation
2 Sechenov First Moscow State Medical University (Sechenov University) Москва, Россия 3 National Research University "Higher School of Economics" Москва, Россия
Annotation
We consider algorithms of 3D reconstruction for the internal surface of cardiac vessels. The precise reconstruction of vessel geometry is necessary for the creating a hydrodynamic model of blood supply for the heart and computing vari© В. В. Борисенко, H. С. Серова, A. M. Чеповский, 2019
ous parameters of blood flow. To compute a triangulation of blood vessel walls, we use the combination of two methods. At the first stage we apply the 3D seeded region growing algorithm to reconstruct a set of voxels inside vessels. At the second stage we use the isosurface reconstruction algorithm based on the tessellation of 3D space into small tet-rahedral cells. We use the tetrahedral mesh, which was proposed in the works of S. Chan, E. Purisima (1998), and V. Skala (2000). Tetrahedra in this mesh are constructed on common faces of adjacent cubes in a cubic lattice, so it fits well with the voxel model. The mesh is constructed only in the neighborhood of the border of voxel set obtained at the first stage as the result of seeded region growing algorithms. Keywords
computed tomography, isosurface triangulation, tetrahedral mesh, 3D seeded region growing Acknowledgements
The research was supported by the RFBR grant no. №18-29-26012 мк For citation
Borisenko V. V., Serova N. S., Chepovsky A. M. Reconstruction of Three-Dimensional Geometry of the Vessels by Computed Tomography Data. Vestnik NSU. Series: Information Technologies, 2019, vol. 17, no. 3, p. 5-17. (in Russ.) DOI 10.25205/1818-7900-2019-17-3-5-17
Введение
Ишемическая болезнь сердца является одним из самых распространенных заболеваний. Ее основной причиной является нарушение кровоснабжения сердечных мышц, связанное, как правило, с сужением коронарных артерий (стеноз). В настоящее время имеются достаточно эффективные методы ее лечения, главным из которых является установка расширителей сосудов (стентов). Для этого в коронарную артерию вводится катетер. Перед проведением подобной операции необходимо вначале убедиться в ее необходимости. Одним из основных критериев является так называемая величина FFR (fraction flow reserve), показывающая, насколько уменьшается кровоток в месте сужения артерии. FFR определяется как отношение давления крови в точках до и после сужения. В медицине используется так называемый «золотой стандарт» измерения FFR, использующий инвазивный метод - введение катетера с датчиком давления в артерию. Несмотря на то, что в настоящее время эта процедура хорошо отработана и может быть совмещена с установкой стента, желательно все-таки иметь возможность определить FFR каким-либо неинвазивным методом, например, используя данные компьютерной томографии области сердца. Это поможет избежать введения катетера в коронарную артерию в случаях, когда нет необходимости в установке стента.
При исследовании кровеносных сосудов с помощью компьютерной томографии в кровь пациента через вену вводится контрастное вещество, резко увеличивающее степень поглощения рентгеновских лучей кровью. Томографические срезы, полученные компьютерным томографом, представляют собой матрицы размером 512 х 512, содержащие двухбайтовые целые числа со знаком. Их величины обозначаются через HU (housfield units), они отражают степень поглощения рентгеновских лучей тканью пациента в данной точке среза. Величина HU воздуха равна -1 024, воды - 0, крови без контрастного вещества - около 40-50. При введении контрастного вещества в кровь ее HU резко возрастает до нескольких сотен, и кровеносные сосуды становятся хорошо различимыми на фоне окружающих тканей.
Данная работа является частью проекта создания компьютерной модели кровоснабжения сердца для конкретного пациента. Основная идея состоит в том, чтобы определить всевозможные параметры кровотока методами гидродинамики, применяя численные методы решения уравнений течения жидкости - уравнения Навье - Стокса, уравнение состояния и др. Эти численные методы предполагают построение геометрии расчетной области в виде задания ее триангуляции, затем вычисление расчетной сетки и численное решение по этой сетке дифференциальных уравнений в частных производных. Предполагается использовать систему FlowVision [1; 2], успешно применяемую в расчетах, необходимых для кораблей, турбин, летательных аппаратов и т. п. Для ее применения в области медицинской диагностики при изучении системы кровообращения требуется на первом шаге точно воссоздать трехмерную
геометрию кровеносных сосудов, в нашем случае аорты и коронарных артерий. В данной работе рассматривается именно эта задача.
Томографические срезы представляют собой последовательность параллельных сечений тела пациента. Современный томограф позволяет за очень короткое время (около 0,35 с) получить серию срезов с очень высоким разрешением: обычно сканируется область размером 16 см по оси Z (вертикальной оси тела пациента) с шагом 0,25 мм по вертикали, при этом горизонтальные размеры пикселя по осям X и Y равны 0,5 мм. Компьютерный томограф практически мгновенно выдает серию из 640 срезов в области сердца. Эти срезы определяют функцию рентгеновской плотности ткани человеческого тела в трехмерном пространстве d(x,y,z) (в промежутках между срезами она вычисляется с помощью интерполяции). Стандарт DICOM (digital imaging and communication in medicine) определяет следующие направления координатных осей: ось X направлена от правой части тела пациента к левой (Right ^ Left), ось Y - от передней стороны тела к задней (Anterior ^ Posterior), ось Z -от ног к голове (Feet ^ Head).
Задачу трехмерной реконструкции в самом простом случае можно сформулировать как построение триангуляции поверхности, заданной уравнением
d(x,y,z) = t,
где d(x,y, z) - функция плотности, t - пороговое значение (threshold). В зависимости от выбора значения t мы восстанавливаем разные структуры человеческого тела - например, при t = 200 получается поверхность не очень твердых костей, при t = 1 000 - зубов, при t = -200 - опухолей в легких и т. п. При восстановлении внутренней поверхности кровеносных сосудов в случае использования контрастного вещества применяется пороговое значение t в пределах HU от 150 до 400. Отметим также, что практически всегда восстановление такой поверхности выполняется не во всем пространстве, а только в интересуемой области, в медицине для этого используется термин ROI (region of interect) или, реже, VOI (volume of interest). ROI обычно задается либо наборами контуров на каждом срезе, либо наборами битовых маск (чаще всего используются оба метода одновременно).
Типы трехмерных алгоритмов
Существует два основных класса алгоритмов трехмерной реконструкции. Первый класс состоит из так называемых ячеечных методов, среди них наиболее известен метод «марширующих кубов» ([3], 1987). В таких методах пространство разбивается на многогранники маленького размера (в несколько пикселей на томографических срезах), например, на маленькие кубики в методе «марширующие кубы» или тетраэдры в методах типа «марширующие тетраэдры» (методы МТ5, MT6 и др.). (Под словом «тетраэдр» всюду в этой статье подразумевается любая треугольная пирамида, не обязательно правильный тетраэдр.) Значение функции
fx,y,z) = d(x,y,z) - t
вычисляется в вершинах этих многогранников. Рассматриваются ребра каждого многогранника. Если в вершинах ребра функция f принимает значения разных знаков, то при допущении, что в точках ребра функция меняется линейно, вычисляется точка на ребре, в которой функция принимает значение 0. Полученные точки на ребрах соединяются одним или несколькими треугольниками внутри рассматриваемого многогранника. Объединяя полученные таким образом точки и треугольники, получаем триангуляцию всей поверхности (см. далее, рис. 4 для случая тетраэдрических ячеек).
Второй класс алгоритмов состоит из так называемых методов «роста области из семени» (seeded region grow), мы рассматриваем их трехмерные варианты. Алгоритм роста области из семени строит множество S вокселей внутри интересуемой области, которая может задаваться пороговым значением t функции d(x, y, z) в трехмерном пространстве (область состоит
из вокселей, значения функции d в которых не ниже порогового), а также какими-либо другими признаками. В интересуемой области выбирается начальный воксель («семя»), он добавляется к множеству S. Затем программа в цикле добавляет к множеству S воксели, соседние с уже содержащимся в S и удовлетворяющие критерию принадлежности к области. В процессе работы алгоритм использует какую-либо динамическую структуру D, содержащую воксели, которые уже были добавлены к множеству S, но соседи которых еще не были рассмотрены. Алгоритм заканчивается, когда D становится пустой. В качестве D можно использовать очередь, стек, список, но наиболее удобной структурой такого рода является структура Deque (double ended queue - двусторонняя очередь), содержащаяся в стандартной библиотеке языка C++ (std::deque). Вот запись трехмерного алгоритма роста области из семени на псевдокоде:
S: множество вокселей; D: очередь вокселей (deque)
p0: начальный воксель ("семя")
S.добавить(p0)
D.добавить в конец (p0)
цикл пока D не пуста выполнять
| p = D. начало; D.удалить начало
| цикл для каждого соседа n вокселя p выполнить
| | если n не принадлежит множеству S и n удовлетворяет
| | | критерию принадлежности к области, то
| | | S.добавить(n)
| | | D.добавить в кожц(п)
| | конец если
| конец цикла
конец цикла
В качестве соседей вокселя p в трехмерном пространстве берутся либо 6 соседей, одна из координат которых отличается ровно на единицу от соответствующей координаты p, либо 26 соседей, представляющих собой все воксели куба размером 3 х 3 с центром p, отличные от p.
Достоинства и недостатки трехмерного алгоритма роста области из семени
Метод роста области из семени имеет совсем простую программную реализацию. Но главным достоинством этого метода является то, что он вычисляет связную область, содержащую начальный воксель. Он хорошо подходит для построения множества вокселей внутри кровеносных сосудов, поскольку такое множество невозможно задать только на основе порогового значения. Например, костная ткань (плотность которой варьируется в очень широких пределах) может иметь такие же значения плотности HU, как и кровь при использовании контрастного вещества, поэтому ячеечными методами отделить кровеносную систему от костной ткани невозможно. Но алгоритм роста области добавляет к трехмерной модели лишь те воксели, которые связаны с начальным вокселем («семенем»). Если начальный во-ксель находится внутри кровеносного сосуда, то и все добавленные воксели будут находиться внутри кровеносной системы (учитывая, что кровеносная система замкнута) - алгоритм роста не должен «протекать» в пространство вне сосудов. Таким образом, костная ткань автоматически исключается из трехмерной модели.
После построения множества вокселей остается вычислить триангуляцию его поверхности. Эта процедура также несложная - рассматриваются те грани «внешних» вокселей из множества S, к которым не примыкают другие воксели из S. Каждая такая прямоугольная грань разбивается диагональю на 2 треугольника, и получаем требуемую триангуляцию поверхности множества S.
Еще одно достоинство алгоритма роста области состоит в том, что несложно ограничить рост области, заполнив контуры на некоторых срезах и запустив алгоритм роста повторно 1.
На рис. 1 изображены результаты восстановления области сердца, полученные трехмерным алгоритмом роста области. Левая модель получена ростом из начального вокселя, расположенного внутри аорты. В правой модели используется задание области интереса (ROI), которая отделяет часть аорты и коронарные артерии от остальной части кровеносной системы; правая модель представляет собой часть модели, изображенной слева.
Рис. 1. Результаты трехмерного восстановления области сердца. Левая модель получена алгоритмом роста из начального вокселя, расположенного внутри аорты. Правая модель включает часть аорты и коронарных артерий, она получена из левой модели применением того же алгоритма роста при установленных ограничениях
Fig. 1. The results of 3D reconstruction of the heart area. The left model is obtained by the 3D region growing algorithm from the initial voxel inside the aorta. The right model includes a part of aorta and coronal arteries, it obtained from the left model using the same algorithm, but taking into account the restrictions that separates coronal arteries from other parts of the circulatory system
Но у алгоритма роста области из семени есть и недостатки. Первым из них является большее время работы, чем у ячеечных алгоритмов, поскольку рассматриваются все воксели трехмерной модели, которая строится, а их число может быть очень большим. В ячеечных алгоритмах можно выбирать размер шага (т. е. размер многогранников, на которые разбивается пространство), но в алгоритме роста области шаг всегда минимально возможный - он равен размеру одного вокселя. Кроме того, ячеечные алгоритмы добавляют к модели лишь треугольники, расположенные на поверхности модели, а алгоритм роста области проходит по всему ее объему. Однако для наших целей этот недостаток не столь существенен, поскольку в любом случае шаг нужно выбирать минимально возможным, чтобы достичь наилучшего разрешения, так как коронарные артерии могут иметь очень небольшой поперечный размер. К тому же современные компьютеры позволяют достичь очень высокой скорости, и вдобавок отсутствуют существенные ограничения на объем оперативной памяти.
Отметим, что второй недостаток алгоритма роста области, увы, не позволяет напрямую его использовать при построении триангуляции для последующего гидродинамического рас-
1 Подробно этот метод описан в работе: Борисенко В. В., Веселова Т. Н, Терновой С. К, Чеповский А. М. Реконструкция трехмерной геометрии коронарных артерий // Фундаментальная и прикладная математика. 2018. Т. 22 (в печати).
чета. Дело в том, что поверхность построенной модели получается негладкой, она составлена из поверхностей кубиков, представляющих собой воксели полученного множества (рис. 2, 3). Если бы цель состояла только в том, чтобы получить изображение поверхности, то эта проблема была бы не столь существенна - трехмерный алгоритм Фонга с использованием нормалей в вершинах триангуляции позволяет визуально сгладить изображаемую поверхность (см. рис. 2). Но при гидродинамическом расчете поверхность нуждается в реальном, а не только визуальном сглаживании. Алгоритмы сглаживания из класса Subdivision Surface (дробление поверхности), такие как алгоритм Лупа [4] или Кетмулла - Кларка [5], также мало подходят для нашей цели, поскольку они, лишь сглаживая углы, в целом не меняют форму поверхности, к тому же в медицинской диагностике желательно избежать даже минимального изменения формы вычисляемой поверхности.
Рис. 2. Триангуляция поверхности множества вокселей, полученного трехмерным алгоритмом роста области. На левом изображении грани вокселей закрашиваются как плоские многоугольники (flat shading). Правое изображение той же модели получено с помощью алгоритма Фонга (smooth shading): использование нормалей к поверхности, вычисленных как градиент функции плотности, позволяет визуально сгладить триангуляцию без изменения ее формы
Fig. 2. Surface triangulation of the set of voxels, obtained by the 3D region growing algorithm. On the left image, voxel faces are drawn as flat polygons (flat shading). The right image of the same model is obtained using the Phong algorithm (smooth shading): the use of surface normals, calculated as a gradient of the density function, visually smooths the triangulation without changing its shape
Рис. 3. Разница между ячеечным алгоритмом триангуляции и трехмерным алгоритмом роста области. Крупно показан один и тот же фрагмент модели, созданный разными алгоритмами
Fig. 3. The difference between the cell triangulation algorithm and the 3D region growing. The same fragment of two models, created by different algorithms, is shown in large
Отметим, что ячеечные алгоритмы строят более гладкую поверхность, поскольку построенные в них треугольники получаются приблизительно перпендикулярными вектору градиента функции плотности, в отличие от метода роста области, в котором плоскости треугольников параллельны координатным плоскостям. Решение проблемы состоит в последовательном применении двух алгоритмов. Сначала с помощью алгоритма роста области из семени вычисляется множество вокселей, составляющее трехмерную модель интересующего участка кровеносной системы 2. На втором этапе используется ячеечный алгоритм вычисления триангуляции, причем он применяется не ко всему пространству, а только к окрестности поверхности модели, вычисленной на первом этапе.
Самым первым ячеечным алгоритмом построения триангуляции изоповерхности был метод марширующих кубов [3]. Он неплох, если речь идет только об изображении трехмерного объекта. Однако быстро выяснился его серьезный недостаток: в некоторых случаях вершины триангуляции на ребрах куба можно соединить треугольниками разными способами. При несогласованном выборе треугольников внутри смежных кубов в построенной триангуляции сплошного объекта могут образоваться дыры. Всего с точностью до симметрии имеется 14 вариантов расположения вершин триангуляции на ребрах куба, неоднозначность возможна в 5 случаях. Вероятность этих ситуаций мала, и они почти не влияют на визуальное восприятие объекта, но если речь идет о точном воспроизведении топологии, то подобное недопустимо. Этого недостатка лишены ячеечные методы, в которых пространство разбивается на тетраэдры. С точностью до симметрии имеется лишь 2 варианта расположения вершин триангуляции на ребрах тетраэдра. В случае трех точек они соединяются одним треугольником, в случае четырех - двумя (рис. 4), других случаев не бывает.
При использовании тетраэдрических ячеек возникает проблема: как выбрать наилучшую тетраэдрическую сеть? Дело в том, что разбить трехмерное пространство на правильные тетраэдры невозможно (в отличие от плоскости, которая разбивается на правильные треугольники). Это означает, что априори «самой лучшей» тетраэдрической сети не существует. Обычно под качеством сети понимают средний аспект входящих в нее тетраэдров. Аспектом тетраэдра называется отношение радиусов описанной и вписанной сфер. Аспект правильного тетраэдра в точности равен 3, у любого тетраэдра, отличного от правильного, аспект строго больше трех, так что аспект может служить мерой того, насколько тетраэдр отличается от правильного. Также желательными свойствами тетраэдрической сети являются равенство входящих в нее тетраэдров, инвариантность сети относительно сдвигов в трех независимых направлениях, удобство программирования и т. п.
Самым простым методом построения тетраэдрической сети является разбиение пространства сначала на кубы, а затем разбиение каждого куба на тетраэдры. Не добавляя дополнительных вершин, куб можно разбить на тетраэдры 6 способами. В методе МТ5 (см. [6]) куб разбивается на 5 тетраэдров, из которых один правильный с объемом 1/3, он образован скрещивающимися диагоналями граней куба; остальные 4 тетраэдра имеют объем 1/6, каждый
2 См.: БорисенкоВ. В. и др. Реконструкция трехмерной геометрии коронарных артерий.
Тетраэдрические сети
Рис. 4. Построение триангуляции внутри тетраэдрической ячейки
Fig. 4. Construction of triangulation inside a tetrahedral cell
из них образован тремя перпендикулярными ребрами куба, выходящими из одной вершины. Это разбиение изображено слева на рис. 5.
Рис. 5. Разбиение куба на тетраэдры (слева и в середине на 5 и 6 тетраэдров в методах MT5 и MT6 соответственно; справа показан процесс построения тетраэдров в методе Скала) Fig. 5. The various partitions of a cube into tetrahedra. The partitions into 5 and 6 tetrahedra are shown on the left and in the middle. The right shows the process of constructing tetrahedra in the Skala method
В методе MT6 (см. [7]) куб разбивается на 6 равных тетраэдров (см. рис. 5, в центре). Сеть MT6 удобна тем, что она инвариантна относительно сдвигов на вектор любого ребра куба. Однако для нашей цели лучше всего подходит сеть, предложенная в работе S. L. Chan, E. O. Purísima [8]; она стала широко известна также и благодаря статье V. Skala [9], в литературе ее часто называют сетью Скала. Тетраэдры этой сети изображены на рис. 5, справа, они строятся на границах смежных кубов кубической решетки. В сети Скала все тетраэдры равны между собой (это тетраэдры Кокстера, сеть инвариантна при отражениях пространства относительно любых граней тетраэдров). Сеть Скала имеет наилучший аспект 3.162 среди всех известных на данный момент тетраэдрических сетей. Для сравнения: аспект сети MT6 равен 4.182, средний аспект сети MT5 равен 3.878.
Использование тетраэдрической сети Chan - Purísima - Skala для сглаживания ЭБ-модели, построенной алгоритмом роста области
Помимо остальных прекрасных свойств сети Chan - Purisima - Skala (равенство входящих в нее тетраэдров, наилучший аспект среди всех известных сетей, инвариантность относительно сдвигов и отражений), она обладает свойством, крайне ценным для нашей задачи: в ней тетраэдры строятся на гранях смежных кубов кубической решетки. Справа на рис. 5 изображены 4 тетраэдра, построенные на общей грани двух смежных кубов, расположенных горизонтально, аналогично строятся тетраэдры на общей грани соседних кубов, расположенных вертикально (и во всех трех независимых направлениях). Вершины построенных тетраэдров находятся в вершинах исходной кубической решетки, а также в центрах кубов. Если размер ребра куба равен 1, то длины ребер тетраэдров равны 1 (два ребра, совпадающие с ребрами куба) и -\/3/2 « 0,866 (остальные 4 ребра, равные половине большой диагонали куба). Если каждый куб решетки представляет воксель в пространстве, то построенная на этих кубах тетраэдрическая сеть Скала позволяет достичь субпиксельной точности, которая необходима для точного восстановления геометрии кровеносных сосудов.
Ключевой момент состоит в том, что тетраэдрическая сеть строится не во всем пространстве, а только в окрестности границы трехмерной воксельной модели, полученной алгоритмом роста области из семени. Благодаря этому при последующем применении ячеечного метода трехмерного восстановления мы не добавляем нежелательных объектов к модели, а лишь уточняем форму ее поверхности.
Точный алгоритм построения «уточняющей» тетраэдрической сети следующий. Обозначим через M исходную модель, полученную трехмерным алгоритмом роста области из семе-
ни, она состоит из кубических вокселей. Каждый воксель пространства - это куб с 8 вершинами. Какие-то вершины вокселей из M могут лежать на поверхности модели M, какие-то нет. Построим множество вокселей следующим образом: произвольный воксель v трехмерного пространства добавляется к 5, если хотя бы одна из его вершин лежит на поверхности модели М. Таким образом, множество 5 будет содержать как некоторые воксели из модели М, так и некоторые воксели, не входящие в М, но «касающиеся» ее границы.
После определения множества вокселей 5 строится уточняющая тетраэдрическая сеть: для любых двух соседних вокселей из 5, имеющих общую грань, строятся 4 тетраэдра сети (см. рис. 5, справа). Схематично процесс построения сети изображен на рис. 6 в двумерной проекции (трехмерный рисунок получился бы слишком громоздким).
Рис. 6. Построение «уточняющей» тетраэдрической сети для точного вычисления сглаженной границы воксельной модели Рис. 6. Construction of a refining tetrahedral mesh for the exact calculation of the smoothed boundary of the voxel model
На рис. 6 темными квадратами показаны воксели исходной модели (множество М), светлыми - воксели из построенного множества 5, не принадлежащие М, но имеющие вершину на границе М. Каждые 2 соседних вокселя из множества 5 определяют 4 тетраэдра, которые строятся на их общей грани, проекции этих тетраэдров нарисованы тонкими линиями.
После построения уточняющей тетраэдрической сети запускается алгоритм трехмерного восстановления изоповерхности ячеечного типа. Результат применения алгоритма для модели тора показан на рис. 7.
Рис. 7. Уточнение границы для воксельной модели тора, полученной трехмерным алгоритмом роста области (слева). Справа показана модель, полученная вычислением границы с использованием уточняющей сети Chan -Purisima - Skala. Для изображения используется OpenGL и плоское закрашивание треугольников (Flat Shading)
Fig. 7. The refinement of the boundary for the voxel model of the torus obtained by the 3D region growing algorithm (left). On the right is the model obtained by calculating the boundary using the Chan - Purisima - Skala refining mesh. OpenGL and Flat Shading are used for the image
В правой части рис. 7 специально используется плоский режим закрашивания треугольников, чтобы они были лучше видны на изображении. Если применить алгоритм Фонга, использующий нормали в вершинах триангуляции для сглаживания изображения, то результат получается визуально совершенно гладким. На рис. 8 уточненная модель изображена слева как проволочный каркас (Wireframe), справа используется алгоритм Фонга.
Рис. 8. Изображение уточненной модели тора как проволочного каркаса (слева) и с помощью алгоритма Фонга, использующего нормали для сглаживания поверхности (справа)
Fig. 8. Image of the refined model of the torus as a wire frame (left) and by using normals to the surface and Phong algorithm for smooth shading (right)
На рис. 9. показан результат применения уточняющего алгоритма для компьютерной модели кровеносного сосуда. Отметим, что толщина сосуда в этом примере минимальна (3-4 вокселя), а исходная воксельная модель имеет форму, далекую от цилиндрической, однако уточненная триангуляция получается достаточно гладкой.
Рис. 9. Уточнение границы для компьютерной имитации кровеносного сосуда. Слева показана воксельная модель, полученная алгоритмом роста области, в центре и справа - триангуляция ее границы, построенная по уточняющей сети Chan - Purisima - Skala. В правом изображении применен алгоритм Фонга (Smooth Shading)
Fig. 9. Boundary refinement for a computer simulated blood vessel. The voxel model obtained by the region growing algorithm is shown on the left, the triangulation of its boundary built with the Chan - Purisima - Skala refining mesh is shown in the center and on the right. The Phong algorithm (Smooth Shading) is applied to the right image
На рис. 10 и 11 показан результат применения алгоритма к реальным данным. Слева изображена воксельная модель, в центре и справа - триангуляция, полученная уточняющим ал-
горитмом (используется плоское и гладкое закрашивание граней). Рисунок 11 крупно показывает небольшой фрагмент модели, представленной на рис. 10.
Рис. 10. Применение уточняющего алгоритма к реальным данным. Слева изображена воксельная модель, в центре и справа - триангуляция, полученная уточняющим алгоритмом (используется Flat Shading и Smooth Shading) Fig. 10. Application of the refinement algorithm to real data. On the left is a voxel model, in the center and on the right is a triangulation obtained by a refinement algorithm (using Flat Shading and Smooth Shading)
Рис. 11. Крупное изображение фрагмента модели, представленной на рис. 10 Fig. 11. A large image of a fragment of the model shown in Fig. 10
Все алгоритмы, описанные в данной работе, реализованы на языке C++ с использованием библиотеки классов Qt и системы разработки приложений QtCreator. Для получения изображений использовалась библиотека трехмерной графики OpenGL. Был реализован как трехмерный алгоритм роста области из семени, так и ячеечный алгоритм восстановления изопо-
верхности, использующий тетраэдрическую сеть Chan - Purisima - Skala. Последний алгоритм реализован в двух вариантах: классический вариант строит сеть во всем пространстве, в варианте уточнения границы сеть строится только в окрестности границы воксельной модели, как описано выше.
Отметим, что алгоритм уточнения границы воксельной модели работает чрезвычайно быстро, поскольку он проходит не по всему пространству, а только по границе модели, до этого вычисленной алгоритмом роста области. Тетраэдры сети имеют размеры, не превышающие размера вокселя, и строятся очень естественно по вокселям исходной модели. Это позволяет достичь необходимой гладкости и субпиксельной точности результирующей модели. Причем в модель не вносятся никакие искажения формы при сглаживании, неизбежные, к примеру, в методах Лупа или Кетмулла - Кларка [4; 5], поскольку все вычисления используют только исходную функцию плотности объекта в трехмерном пространстве, полученную как результат томографического обследования пациента. Таким образом, комбинация этих двух алгоритмов позволяет точно и эффективно воссоздать сглаженную трехмерную модель сосудов кровеносной системы, необходимую для программ гидродинамических расчетов.
Список литературы / References
1. Аксенов А. А., Гудзовский А. В. Пакет прикладных программ FlowVision // Тр. МФТИ. Серия: Аэрофизика и прикладная математика. М., 1998. С. 45-56.
Aksenov A. A., Gudzovskiy A. V. Paket prikladnikh program FlowVision. In: Trudy MFTI. Ser. Aerofizika i prikladnaya matematika. Moscow, 1998, p. 45-56. (in Russ.)
2. Аксенов А. А. FlowVision: индустриальная вычислительная гидродинамика // Компьютерные исследования и моделирование. 2017. Т. 9, № 1. C. 5-20.
Aksenov A. A. FlowVision: industrialnaya vichislitelnay gidrodinamika. Komputernye issledovaniya i modelirovanie, 2017, vol. 9, no. 1, p. 5-20. (in Russ.)
3. Lorensen W. E., Cline H. E. Marching Cubes: A high resolution 3D surface construction algorithm. Computer Graphics, 1987, vol. 21, no. 4.
4. Loop Ch. Smooth Subdivision Surfaces Based on Triangles. M.S. Mathematics thesis. University of Utah, 1987.
5. Catmull E., Clark J. Recursively generated B-spline surfaces on arbitrary topological meshes. Computer-Aided Design, 1978, vol. 10, iss. 6, p. 350-355.
6. Carneiro B. P., Kaufman A. E. Tetra-Cubes: An algorithm to generae 3D isosurfaces based upon tetrahedra. In: SIGGRAPH, 1996, p. 205-210.
7. Gueziec A. Exploiting Triangulated Surface Extraction using Tetrahedral Decomposition. IEEE Transactions on Visualization and Computer Graphics, 1995, vol. 1, iss. 4, p. 328-342.
8. Chan S. L., Purisima E. O. A new tetrahedral tesselation scheme for isosurface generation. Computers & Graphics, 1998, vol. 22, iss. 1, p. 83-90.
9. Skala V. Precision of Isosurface Extraction from Volume Data and Visualization. In: Conference on Scientific Computing, 2000, p. 368-378.
Материал поступил в редколлегию Received 22.02.2019
Сведения об авторах / Information about the Authors
Борисенко Владимир Витальевич, кандидат физико-математических наук, ведущий научный сотрудник механико-математического факультета Московского государственного университета им. М. В. Ломоносова (Ленинские Горы, 1, Москва, 119991, ГСП-1, Россия)
Vladimir V. Borisenko, PhD (Math.), Leading Research Scientist, Lomonosov Moscow State University (1 Leninskie Gory, Moscow, 119991, Russian Federation)
vladimir_borisen@mail.ru
Серова Наталья Сергеевна, член-корреспондент РАН, доктор медицинских наук, профессор, директор Института электронного медицинского образования, Первый МГМУ им. И. М. Сеченова (Сеченовский Университет) (ул. Трубецкая, д. 8, стр. 2, Москва, 119991, Россия)
Nataliya S. Serova, Corresponding Member of the Russian Academy of Sciences, Dr. Sc. (Med.), Professor, Sechenov First Moscow State Medical University (Sechenov University) (8/2 Tru-betskaya Str., Moscow, 119991, Russian Federation) dr.serova@yandex.ru
Андрей Михайлович Чеповский, доктор технических наук, профессор, профессор НИУ
Высшая Школа Экономики (ул. Мясницкая, 20, Москва,101000, Россия) Andrey M. Chepovskiy, Dr. Sc. (Eng.), Professor, National Research University Higher School of Economics (20 Myasnitskaya Str., Moscow, 101000, Russian Federation)
achepovskiy@hse.ru
