Санкт-Петербургский государственный университет
Фундаментальные информатика и информационные технологии
Математическое и программное обеспечение вычислительных машин,
комплексов и компьютерных сетей
Монькин Сергей Александрович
Модуль расчета деформации трехмерной модели мягкой ткани моделью
несжимаемого имплантата для компьютерной системы планирования
хирургических операций
Магистерская диссертация
Научный руководитель:
д.ф.-м.н., проф. Терехов А.Н.
Рецензент:
инженер-программист Николаев С.Н.
Санкт-Петербург
2016
SAINT-PETERSBURG STATE UNIVERSITY
Fundamental Computer Science and Information Technologies
Software of Computers, Complexes and Networks
Monkin Sergei
Calculation module of 3d soft tissue model deformation by incompressible implant
model for surgery planning computer system
Master’s Thesis
Scientific supervisor:
Professor Terekhov A.N.
Reviewer:
software engineer Nikolaev S.N.
Saint-Petersburg
2016
Содержание
1Введение.......................................................................................................................................4
2Постановка задачи.......................................................................................................................6
3Обзор.............................................................................................................................................7
3.1Помещение имплантата под мягкие ткани........................................................................7
3.2Подходы к физическому моделированию.........................................................................8
3.3Модель масс и пружин........................................................................................................9
3.4Симуляция..........................................................................................................................11
3.5Обработка столкновений...................................................................................................11
3.6Обзор существующих решений........................................................................................13
4Модель........................................................................................................................................14
4.1Модель масс и пружин......................................................................................................14
4.2Иерархии ограничивающих объемов...............................................................................16
4.3Импульсный метод............................................................................................................18
5Архитектура...............................................................................................................................19
5.1Система планирования хирургических операций...........................................................19
5.2Представление трехмерных объектов..............................................................................20
5.3Архитектура........................................................................................................................20
6Библиотека физических расчетов............................................................................................23
6.1Структура библиотеки.......................................................................................................23
6.2Модель масс и пружин......................................................................................................24
6.3Иерархия ограничивающих объемов...............................................................................24
6.4Обработка столкновений...................................................................................................25
6.5Примеры..............................................................................................................................26
7Модуль помещения имплантата...............................................................................................31
7.1Описание модуля................................................................................................................31
7.2Симуляция..........................................................................................................................32
7.3Примеры..............................................................................................................................32
8Результаты..................................................................................................................................37
9Список литературы....................................................................................................................38
3
1 Введение
При планировании различных операций в медицине широко применяют
компьютерные системы. Например, в пластической и реконструктивной
хирургии при планировании определенных хирургических операций по
снимкам (например, томографическим) строят трехмерную модель части тела
пациента, по полученной модели производят различные замеры, имитируют
операцию, результат используют для принятия решения о проведении самой
операции, демонстрации ее возможного результата пациенту.
Современные компьютеры позволяют имитировать операции с достаточной
точностью непосредственно во время консультации. При планировании
определенных операций необходимо имитировать физическое поведение
тканей тела человека. Представляет определенный интерес физическая
симуляция в реальном времени, например, для визуального наблюдения за
моделируемым поведением или проведения серии моделирований с разными
параметрами.
В реконструктивной хирургии занимаются восстановлением утраченных или
деформированных
частей
организма
человека.
Тканевое
растяжение
применяют для устранения различных дефектов покровных тканей,
например, для восстановления мягких тканей головы и шеи при ожоговых
поражениях, восстановлении утраченного волосяного покрова головы.
Предварительно
наращивают
некоторый
избыток
кожи,
например,
баллонным растяжением. При планировании таких операций представляет
интерес вычисление необходимой площади для пересадки, вычисление
необходимого объема имплантата для растяжения определенной площади,
вычисление возможных вариантов пересадки различных фрагментов тканей,
моделирование возможного конечного результата, внешнего вида.
Целью работы является разработка модуля симуляции поведения мягких
4
тканей и применение этого модуля к симуляции деформации трехмерной
модели
мягкой
ткани
моделью
несжимаемого
имплантата.
Модуль
физических расчетов разрабатывается как часть компьютерной системы
планирования
хирургических
операций
Phoenixcas
3D
Viewer
для
расширения ее функциональности. Разрабатываемым модулем симулируется
результат растяжения (наращивание) кожи для пересадки. В рамки задачи
входит задание параметров имплантата, имитация помещения имплантата
под мягкие ткани для получения определенной площади при растяжении.
5
2 Постановка задачи
Целью работы является разработка модуля, позволяющего симулировать
возможный результат помещения имплантата под мягкие ткани по
трехмерной модели поверхности тела пациента.
Для этого необходимо выполнить следующие задачи:
1. определить архитектуру модуля
2. разработать библиотеку физических расчетов
3. разработать модуль помещения имплантата
В рамках поставленной задачи стоит ограничение на время выполнения
физических расчетов. Время работы модуля не должно превышать 30 секунд
на моделях пациентов системы планирования. Также стоит задача имитации
сдавливания имплантата при помещении.
6
3 Обзор
3.1 Помещение имплантата под мягкие ткани
В [1] приведен обзор подходов пластики дефектов мягких тканей, анализ
литературы, посвященной дермотензии, представлен в [2]. Тканевое
растяжение в восстановительной хирургии применяют для устранения
различных дефектов покровных тканей. Например, для восстановления
мягких тканей головы и шеи при ожоговых поражениях, восстановлении
утраченного волосяного покрова головы. При наращивании ткани не меняют
своих свойств, толщину.
Для проведения операций предварительно наращивают некоторый избыток
кожи,
например,
баллонным
растяжением.
Баллонное
растяжение
заключается в помещении под ткани имплантата (баллончика) определенной
формы и размера и последующем дозированном заполнении его жидкостью в
течение нескольких недель до определенного объема, получая увеличение
площади растягиваемых тканей в несколько раз. Получаемый избыток кожи
используют для пересадки или устранения дефекта рядом с наращенными
тканями. Также устанавливают баллоны различной емкости с разных сторон
от дефекта.
Имплантат представляет собой силиконовый баллончик определенной
формы,
размера
и
объема
с
трубкой
и
клапаном,
заполняемый
физиологическим раствором.
Помещение имплантата под мягкие ткани моделируют различными
способами.
В
[3]
описан
параметрический
подход,
моделирование
происходит последовательностью деформаций исходной области. В [4] для
моделирования
тканевого
растяжения
применяют
метод
конечных
элементов. В [5] описана модификация модели масс и пружин для
моделирования помещения имплантата под мягкие ткани.
7
3.2 Подходы к физическому моделированию
Подходы к моделированию физического поведения различных объектов
приведены в [6]. Для проведения симуляции объекту задают решетку. По
способам задания решетки рассматривают методы Эйлера, в которых
решетка привязана к пространству, и методы Лагранжа, в которых решетка
привязана к объекту. Методы Лагранжа также делят по фиксированной и
свободной топологиям. К решеткам, привязанным к пространству, относится,
например, воксельная модель. К решеткам, привязанным к объекту,
относятся, например, метод конечных элементов, модель масс и пружин.
Распространенными подходами для моделирования физического поведения
различных объектов являются метод конечных элементов (Finite Element
Method (FEM)), модель масс и пружин, модель SPH (Smoothed Particle
Hydrodynamics). Метод конечных элементов применяют для моделирования
поведения объемных объектов, модель масс и пружин для, например,
одежды, волос, модель SPH для жидкостей, газов. Метод конечных
элементов является наиболее точным и сложным в вычислениях. Модель
масс и пружин отличается простотой реализации, меньшей сложностью
вычислений.
В методе конечных элементов объект задают множеством конечных
элементов,
например,
кубов.
Деформация
объекта
описывается
перемещением узловых точек между конечными элементами. Задают
описывающие перемещения полиномы с неизвестными коэффициентами.
Строят систему уравнений, решая ее, получают перемещения узловых точек.
В модели масс и пружин объект задают множеством частиц и соединяющих
их пружин. Частицам задают массы и позиции, пружинам — длину и
коэффициент жесткости. Пружинами описывают взаимодействие между
парами частиц. Симуляция происходит итеративно применением различных
сил к частицам. В модели SPH объект представляют множеством частиц.
8
Частицам задают свойства объекта и области влияния (функции ядра).
Симуляция выполняется итеративно.
Для задания значений параметров модели (например, для модели масс и
пружин — масс частиц, коэффициентов жесткости пружин) выделяют два
подхода — по некоторым известным данным (data-driven) и аналитически [7].
Первый подход заключается в том, чтобы подобрать такие параметры, при
которых симуляция повторяет результат некоторой известной деформации.
Второй подход заключается в определении аналитических выражений для
параметров, используя физическую модель.
3.3 Модель масс и пружин
Объект задают множеством частиц и соединяющих их пружин. Частицам
задают массы и позиции, пружинам — длину и коэффициент жесткости.
Пружинами описывают взаимодействие между парами частиц. Симуляция
происходит итеративно применением различных внешних (external) и
внутренних (internal) сил к частицам. К внутренним относят силы,
воздействующие на частицы при натяжении пружин системы, к внешним
относят,
например,
гравитацию,
трение,
силы,
возникающие
при
столкновениях. Выделяют position based [8] подход к симуляции.
Массы частицам задают, например, константой для каждой частицы или
константой на весь объект. Во втором случае массу всего объекта
распределяют на частицы пропорционально площадям состоящих из этих
частиц полигонов или объемам включающих в себя эти частицы тетраэдров.
Общую массу объекта можно вычислить как произведение плотности на
объем объекта. Если объект представляет собой поверхность, общую массу
можно вычислить, задав объем через произведение площади поверхности на
ее толщину.
9
Коэффициенты жесткости пружин задают, в частности, для решетки
определенной топологии [9] или с использованием генетических алгоритмов
[10].
В
[11]
для
пружин,
составляющих
треугольную
топологию,
коэффициент жесткости задают через отношение площади содержащего
пружину треугольника к длине пружины.
После обновления скоростей частиц умножением скорости на коэффициент
от 0 до 1 вводят затухание (damping).
Задавая определенным образом ограничения (constraints), имитируют
некоторое поведение модели. Например, для имитации несжимаемости
вводят специальные пружины, воздействующие на частицы при изменении
объема объекта. Также для различных задач вводят таким образом
ограничения на сохранения углов между пружинами, расстояний между
парами частиц, площади поверхности, формы объекта.
В [12] описаны ограничения для сохранения площади поверхности через
сохранение площадей треугольников, составляющих поверхность. При
изменении площади треугольника к его вершинам добавляют силы по
направлениям ребер треугольника. Сохранение объема объекта описано
через сохранение объемов тетраэдров, состоящих из треугольников
поверхности и некоторой выбранной точки в пространстве. При изменении
объема тетраэдра к его вершинам добавляют силы по направлениям ребер.
В [13] описано сохранение объема добавлением специальных пружин от
треугольников поверхности к некоторой точке в пространстве. При
изменении
общего
объема
объекта
эти
пружины
добавляют
пропорциональные своим текущим длинам силы вершинам поверхности.
Таким образом объект может значительно менять свою форму, но сохранять
объем.
10
3.4 Симуляция
Методы симуляции физического поведения объектов рассмотрены в [6]. К
ним относятся, например, одношаговые и многошаговые методы, методы
прогноза и коррекции, методы Верле. Методы, использующие текущую
скорость и положение, называют явными (explicit methods). Методы, в
которых для обновления на текущем шаге используют значения в конце
шага, называют неявными (implicit methods). Например, явный, неявный
методы Эйлера первого порядка точности, метод Рунге-Кутты четвертого
порядка.
Для вычисления скоростей и позиций задают некоторый размер шага
симуляции. Чем больше шаг — тем быстрее происходит симуляция, но хуже
стабильность. При маленьком шаге стабильность лучше, но ниже скорость
симуляции. В [14] описывают динамическое изменение размера шага
симуляции. Пока симуляция стабильна, шаг увеличивают. Если на какой-то
итерации
система
перестает
быть
стабильной,
шаг
уменьшают
и
пересчитывают итерацию с меньшим шагом, затем некоторое число итераций
шаг не увеличивают.
Симуляцию проводят некоторое заранее определенное количество шагов или
до достижения завершающего условия, например, минимума суммарного
изменения позиций всех частиц.
3.5 Обработка столкновений
Для
моделирования
физического
взаимодействия
между
объектами
необходимо определять пересечения (столкновения) этих объектов и
разрешать их, т.е. воздействовать на объекты таким образом, чтобы
минимизировать
количество
и
глубину
пересечений.
Для
объектов,
состоящих из вершин и полигонов, обнаружение столкновений —
11
определение пересекающихся полигонов между объектами, разрешение
столкновений — "расталкивание" вершин.
Обнаружение столкновений является вычислительно сложным этапом.
Описание различных подходов к обнаружению столкновений представлено в
[15], [16]. Выделяют две фазы — широкую и узкую. В широкой фазе
происходит отсев непересекающихся пар полигонов. В узкой фазе — точный
поиск столкновений среди оставшихся пар. В узкой фазе также получают
данные для разрешения столкновения, например, глубину пересечения,
направление расталкивания. Используют иерархии ограничивающих объемов
(Bounding
Volume
Hierarchies),
пространственное
разбиение
(Spatial
Partitioning), иерархии BSP деревьев (BSP Tree Hierarchies), алгоритм ЛинаКэнни, алгоритм Гильберта-Джонсона-Керти.
Представление объекта ограничивающими объемами является одним из
распространенных
подходов.
Задают
ограничивающую
геометрию,
например, сферы, выравненные по осям прямоугольные параллелепипеды
(axis-aligned bounding box (AABB)), ориентированные прямоугольные
параллелепипеды
(oriented
bounding
box
(OBB)).
Каждый
объект
представляют множеством ограничивающих объемов. По ограничивающим
объемам строят иерархию — различные деревья, например, АВЛ деревья,
восьмеричные деревья (Octree), k-d деревья. Поиск столкновений в узкой
фазе проводят по тем полигонам объектов, ограничивающие объемы которых
пересекают ограничивающие объемы других объектов.
Разбиение
пространства
воксельной
моделью
также
является
распространенным подходом. Задают некоторый размер вокселя, для
полигонов объектов вычисляют позиции в воксельном пространстве. Удобно
определять только те воксели, в которых находится объект. Таким образом
воксели, в которых находятся полигоны данного объекта, определяют
данный объект с точностью до размера вокселя. Поиск столкновений в узкой
12
фазе проводят только в тех вокселях, в которых находятся полигоны разных
объектов.
Разрешение столкновений рассмотрено в [17]. Среди методов разрешения
столкновений различают методы, основанные на ограничениях (constraint
based), импульсные методы (impulse based) [18], методы штрафных сил
(penalty based).
3.6 Обзор существующих решений
Среди коммерческих продуктов для планирования в пластической хирургии
можно выделить Canfield [19]. В Canfield представлены решения для задач
восстановления трехмерных моделей, визуализации, проведения замеров,
измерения
разницы
объемов,
анализа
состояния
кожи,
симуляции
результатов различных операций на лице и теле.
Среди коммерческих инструментов физического моделирования можно
выделить ANSYS [20], ABAQUS [21]. Эти продукты разработаны для
проектирования различной техники, анализа прочности конструкций,
вязкоупругих материалов, теплопереноса. Симуляция основана на методе
конечных элементов.
Среди инструментов с открытым исходным кодом для моделирования
поведения мягких тел в задачах медицины можно выделить SOFA (Simulation
Open Framework Architecture) [22], FEBio [23]. SOFA написан на C++,
поддерживает различные физические модели (масс и пружин, конечных
элементов, SPH), методы обработки столкновений, схемы интегрирования,
свободен для коммерческого использования. FEBio основан на методе
конечных элементов, свободен для некоммерческого использования.
13
4 Модель
В качестве основы физической модели выбрана модель масс и пружин.
Модель масс и пружин проста в реализации, в поставленной задаче позволяет
достигнуть достаточной скорости физических расчетов. Для проведения
симуляции выбран явный метод Эйлера. Для обнаружения столкновений
выбраны иерархии ограничивающих объемов. Иерархии ограничивающих
объемов
позволяют
ускорить
время
обработки
столкновений.
Для
разрешения столкновений выбран импульсный метод.
4.1 Модель масс и пружин
Модель масс и пружин описана в [6]. Объект задается множеством частиц и
соединяющих их пружин. Частицам задаются массы и позиции, пружинам —
длина и коэффициент жесткости. Симуляция происходит итеративно
применением различных сил к частицам.
Масса mi частицы i через объемы содержащих частицу тетраэдров v i ,
объем всего объекта V и общую массу M задается следующим образом:
mi= M
∑i∈t vt
.
(1)
3V
Через площади содержащих частицу полигонов a i и общую площадь
поверхности A масса частицы задается аналогично:
mi= M
∑i∈t at
.
(2)
3A
Общая масса объекта вычисляется как произведение плотности на объем
объекта. Для объекта, представляющего собой поверхность, общая масса
вычисляется при объеме, заданном как произведение площади поверхности
на ее толщину.
14
Для
пружин,
составляющих
треугольную
топологию,
коэффициент
жесткости k c задается через отношение площади содержащего пружину
треугольника к длине пружины c следующим образом [11]:
( )
( )
E 2 2 A abc
E 2 v a 2+b2−c 2
,
k c=
+
1+v
4 A abc
c2
1−v 2
где E 2=Et , E —
модуль
Юнга, t —
толщина
(3)
ν —
пружины,
коэффициент Пуассона, A abc — площадь треугольника со сторонами
a , b, c , каждая из которых является пружиной. Общий коэффициент
жесткости пружины — сумма коэффициентов жесткости, вычисленных по
треугольникам, в которые входит пружина. Таким же образом коэффициенты
жесткости задаются для остальных пружин.
Для пружины между частицами i, j сила воздействия пружины на частицу i
описывается законом Гука:
f i=k s (|x ij|−l ij )
x ij
,
|x ij|
(4)
где k s — коэффициент жесткости пружины, |x ij|=x j− x i , xi, xj – позиции
частиц i, j, |x ij| — текущая длина пружины,
lij — начальная длина
пружины.
После обновления скоростей частиц умножением скорости на коэффициент
от 0 до 1 вводится затухание (damping).
Для сохранения объема вводятся дополнительные пружины от частиц
поверхности к центру масс объекта [13]. Для частицы поверхности i сила,
действующая для сохранения объема:
f i=−c i k v (Δ l tot +C Δ l i )u i ,
где k v —
коэффициент
Δ l i=li −l 0i , li —
текущая
(5)
1
жесткости, 0<C <1 , Δ l tot =
3
длина
15
дополнительной
ns
∑ ci Δ li
,
i=1
0
пружины, li
—
начальная длина дополнительной пружины, ui —
направление от центра
объекта к частице, значение c i вычисляется через треугольники, в которые
ni
входит частица, как c i =∑
j=1
Sj
, где S j – площадь треугольника j , S AV
S AV
– средняя площадь.
Явный метод Эйлера для обновления скорости и позиции частицы:
fi
h ,
mi
(6)
x i+1=x i +v i h ,
(7)
v i+ 1=v i +
где v i+ 1 , x i+1 — скорость и позиция частицы на новом шаге симуляции,
v i , x i — скорость и позиция частицы на текущем шаге симуляции, f i —
сумма сил, действующих на частицу, mi — масса частицы, h — заданный
шаг симуляции.
4.2 Иерархии ограничивающих объемов
Иерархии ограничивающих объемов (Bounding Volume Hierarchies) подробно
описаны в [15], [24].
Иерархию ограничивающих объемов над множеством объектов определяют
следующим образом. Если число объектов не превышает какого-то заданного
числа, то иерархия является листом, содержащим множество объектов и
ограничивающий объем над ним. Если число объектов больше заданного
числа, то иерархия является узлом, содержащим ограничивающий объем над
множеством
объектов
и
множество
узлов
нижнего
уровня
над
подмножествами множества объектов, каждый из которых является
иерархией.
К ограничивающим объемам относят, например, сферы, выровненные по
16
осям прямоугольные параллелепипеды (Axis Aligned Bounding Box —
AABB),
ориентированные
прямоугольные
параллелепипеды
(Oriented
Bounding Box — OBB), дискретно-ориентированные многогранники (Discrete
Oriented Polytope — k-DOP).
Выделяют три подхода к построению иерархий ограничивающих объемов:
снизу-вверх (bottom-up), сверху-вниз (top-down), вставкой (insertion). При
построении
снизу-вверх
сначала
для
каждого
объекта
строится
ограничивающий объем, каждый объект представляется листом. Затем по
некоторому критерию происходит объединение листов в узлы более
высокого уровня. Такое объединение производится с узлами каждого уровня
до достижения корневого узла. При построении вставкой для каждого
объекта вставка в узел происходит рекурсивно следующим образом. Если
узел является листом, объект помещается в него, иначе среди подузлов узла
выбирается узел согласно некоторому критерию.
При построении сверху-вниз сначала строится ограничивающий объем по
всему множеству объектов — корневой узел. Узел разбивается согласно
некоторому определенному критерию на множество узлов нижнего уровня,
затем в каждом таком узле построение иерархии происходит рекурсивно.
Критерием разбиения может быть, например, разделение пространства на
равные части, минимизация суммы объемов или площадей подузлов,
минимизация объема пересечения объемов подузлов. Исходя из того, чтобы
при проверке пересечения узла лучом вероятность пересечения подузлов
лучом была одинаковой, узел разбивают таким образом, чтобы общая
площадь поверхности в одном подузле совпадала с общей площадью
поверхности в другом подузле (surface area heuristics). Также разбиение узла
на подузлы производят до достижения некоторой определенной глубины
разбиения или до достижения какого-то минимального числа объектов в
узле.
17
Обнаружение столкновений между двумя иерархиями ограничивающих
объемов рекурсивно, начиная с корневых узлов, описывается следующим
образом. Если узлы являются листьями, производится точная проверка
пересечений объектов в листьях; если узлы пересекаются, попарно
проверяются столкновения между подузлами как между узлами.
4.3 Импульсный метод
Разрешение столкновения двух объектов импульсным методом
[18]
заключается в вычислении импульса в точке контакта и применении его к
обоим объектам для изменения их скоростей и направлений движений.
Изменяют линейные и угловые скорости объектов, также вводят трение.
Для изменения линейной скорости двух объектов импульс j в точке
контакта:
j=
−(1+e)(v 2−v 1)⋅n
,
−1
m−1
+m
1
2
(8)
где e — коэффициент восстановления скорости, 0⩽e⩽1 , v 1 , v 2 —
скорости объектов, n — направление удара, m1 , m2 — массы объектов.
Линейные скорости объектов после столкновения:
v 1=v 1+ j
n
,
m1
(9)
v 2=v 2− j
n
.
m2
(10)
18
5 Архитектура
5.1 Система планирования хирургических операций
Система планирования хирургических операций разрабатывается в среде
Microsoft
Visual
Studio
2010.
Архитектура
системы
планирования
хирургических операций разделена на ряд уровней. Среди них можно
выделить уровень модулей, уровень общих компонент и уровень сторонних
библиотек. Структура компонент системы представлена на Рис. 1.
Рис. 1: Структура компонент системы планирования
Модули
первого
пользовательские
уровня
задачи.
(Modules)
Модули
реализуют
разделены
на
определенные
логику
(Logic)
и
графический интерфейс пользователя (UI). Каждый модуль реализует логику
своей работы и пользовательский интерфейс. Логика написана на языке C++,
графический интерфейс на HTML, Javascript. На уровне общих компонент
(Shared Components) реализованы алгоритмы (Algorithms) и структуры
данных (Data), доступные всем модулям. Уровень сторонних библиотек
19
(External Libraries) закрыт от прямого доступа для модулей, работа с ним
производится на уровне общих компонент. В системе реализована, в
частности, визуализация (Visualization), ввод, вывод данных, работа с
потоками.
5.2 Представление трехмерных объектов
Трехмерные модели пациентов для системы планирования строятся с
помощью
пассивной
фотограмметрии.
Эти
модели
являются
поверхностными незамкнутыми (не имеющими объема) полигональными
сетками. Не имеют дублирующихся вершин, дыр в поверхности.
Полигональные сетки задаются списками вершин и треугольников. Вершины
задаются тремя координатами в трехмерном пространстве, треугольники —
индексами вершин. Площадь поверхности вычисляется как сумма площадей
треугольников, составляющих поверхность. Объем замкнутой трехмерной
модели вычисляется как сумма объемов тетраэдров, построенных от
некоторой точки внутри объекта к треугольникам поверхности. Объемы
складываются с знаками плюс и минус в зависимости от направления
нормали треугольника относительно выбранной точки.
5.3 Архитектура
Библиотека физических расчетов разрабатывается на уровне общих
компонент,
модуль
помещения
имплантата
под
мягкие
ткани
разрабатывается на уровне модулей. Библиотека физических расчетов
предоставляет
имплантата
интерфейс
описывается
к
на
физической
уровне
модели,
модуля.
логика
помещения
Выделение
проведения
симуляции модулем в отдельный поток определяется длительностью
выполнения симуляции и возможностью взаимодействия пользователя с
20
системой во время выполнения симуляции, например, для того, чтобы
отменить
симуляцию.
имплантата
Основная
определена
в
логика
главном
работы
потоке,
модуля
помещения
выполнение
симуляции
определено в новом потоке.
Симуляция мягких тканей реализуется моделью масс и пружин, обнаружение
столкновений реализуется через иерархии ограничивающих объемов,
разрешение столкновений — импульсным методом. Структура классов
библиотеки
физических
расчетов
и
модуля
помещения
имплантата
представлена на Рис. 2.
Рис. 2: Структура классов модуля
Модуль
помещения
имплантата
состоит
из
нескольких
классов,
реализующих графический интерфейс для взаимодействия с пользователем и
основную логику работы. На диаграмме представлены классы Application и
ImplantInsertionModule этого модуля. Класс Application реализует основную
логику работы модуля, класс ImplantInsertionModule реализует помещение
21
имплантата взаимодействием с классами библиотеки физических расчетов.
Архитектура остальных классов модуля определена системой планирования,
поэтому здесь они не представлены. Визуализация, ввод, вывод данных здесь
также не представлены. Библиотека физических расчетов представлена
классом модели масс и пружин (MSS), описывающим симуляцию мягких
тканей, классами иерархии ограничивающих объемов (BVH) и методов
обработки столкновений (CollisionSolver), описывающими обнаружение и
обработку столкновений. Класс ImplantInsertionModule содержит в себе
объекты классов MSS, BVH, CollisionSolver для проведения симуляции, класс
CollisionSolver реализует обработку столкновений, используя данные
объектов классов MSS, BVH. Для выполнения симуляции работа класса
ImplantInsertionModule вынесена в другой поток.
22
6 Библиотека физических расчетов
6.1 Структура библиотеки
Структура классов библиотеки физических расчетов представлена на Рис. 3.
Модель масс и пружин реализуется классами MSS, Particle, Spring, иерархия
ограничивающих объемов — классами BVH, Node, Leaf, обработка
столкновений — классом CollisionSolver.
Рис. 3: Структура классов библиотеки физических расчетов
Шаг симуляции состоит из следующих действий:
1. вычисление сил, воздействующих на частицы;
2. вычисление новых скоростей частиц;
3. обработка столкновений;
4. вычисление новых позиций частиц;
5. проверка критерия остановки симуляции.
Библиотека предоставляет интерфейс к методам для выполнения каждого
действия. Методы для выполнения действий 1, 2, 4, 5 определены в классе
MSS, для обновления иерархии ограничивающих объемов — в классе BVH,
для обнаружения и разрешения столкновений — в классе CollisionSolver.
23
6.2 Модель масс и пружин
Класс MSS реализует модель масс и пружин, методы симуляции модели,
методы задания параметров модели, содержит списки частиц (Particle) и
пружин (Spring).
Класс Particle реализует частицы модели масс и пружин, содержит массу и
позицию частицы. Класс Spring реализует пружины модели масс и пружин,
содержит индексы частиц, начальную длину, коэффициент жесткости.
Массы частиц задаются по формулам (1), (2), коэффициенты жесткости
пружин по формуле (3). Возникающие при натяжении пружин силы
вычисляются по формуле (4). Силы для сохранения объема вычисляются по
формуле (5). Симуляция производится явным методом Эйлера (6), (7).
Введено затухание после вычисления скоростей.
Шаг симуляции задан константой. Остановка симуляции производится при
достижении некоторого минимума суммарного перемещения частиц. На
каждом шаге вычисляется сумма перемещений позиций частиц. Если за
последние некоторое количество шагов отношение минимального значения к
максимальному
близко к
единице
(меньше заданной
погрешности),
симуляция считается завершенной.
6.3 Иерархия ограничивающих объемов
Класс BVH реализует иерархию ограничивающих объемов, методы
построения и обновления, содержит списки узлов (Node) и листьев (Leaf).
Класс Node реализует узлы иерархии ограничивающих объемов, содержит
указатели на подузлы, позицию и размеры ограничивающего объема узла.
Класс Leaf реализует листья иерархии ограничивающих объемов, содержит
список полигонов, находящихся в листе.
24
Разделено построение и обновление иерархии. Построение происходит перед
симуляцией, обновление — на каждом шаге симуляции.
Геометрия
ограничивающего
прямоугольный
объема
параллелепипед.
—
Иерархия
выровненный
по
ограничивающих
осям
объемов
строится сверху вниз. При построении каждый узел, начиная с корневого,
содержащего все треугольники объекта, разбивается на два подузла таким
образом, что сумма площадей треугольников в одном узле равняется сумме
площадей треугольников в другом. Разбиение на подузлы происходит, пока
количество треугольников в узле более заданного числа. При количестве
треугольников менее заданного числа дальнейшее разбиение не происходит,
строится лист.
Во
время
Обновление
симуляции
изменяются
ограничивающих
позиции
объемов
узлов
вершин
треугольников.
и листьев
происходит
рекурсивно, начиная с корневого узла. Для узла, не являющегося листом,
рекурсивно обновляются подузлы и если хотя бы для одного подузла
ограничивающий объем был пересчитан, пересчитывается ограничивающий
объем узла. Если узел является листом, проверяется, отличается ли
ограничивающий объем листа от ограничивающего объема треугольников,
входящих в лист. Если да, ограничивающий объем листа пересчитывается.
6.4 Обработка столкновений
Класс
CollisionSolver
реализует методы
обнаружения
и
разрешения
столкновений. Обнаружение столкновений — между парой трехмерных
моделей, представленных двумя иерархиями, между иерархией по одной
модели и воксельным представлением по второй модели. Разрешение
столкновений реализуется импульсным методом.
Обнаружение столкновений между иерархией ограничивающих объемов и
25
воксельным представлением реализуется исходя из того, что одна
трехмерная
модель
имеет
значительно
большее
число
вершин
и
треугольников, чем вторая. Для второй модели достаточно воксельного
представления.
Обнаружение
столкновений
производится
аналогично
обнаружению столкновений между двумя иерархиями. Для двух трехмерных
моделей иерархия построена для первой модели, воксельное представление
для второй. Для каждого вокселя второй модели проверяется пересечение с
узлами иерархии первой модели, начиная с корневого. При достижении листа
иерархии,
производится
точная
проверка
пересечений
между
треугольниками.
Проверка пересечений между парой треугольников производится проверкой
пересечений между ребрами треугольника одной модели с треугольником
второй модели. По ребрам треугольников построены пружины. При
пересечении ребра с треугольником вычисляется импульс по формуле (8),
изменяются скорости частиц треугольников по формулам (9), (10).
6.5 Примеры
Модель сферы продавливает модель плоскости снизу. Центр плоскости
находится в нулевых координатах, сфера смещена вниз на некоторое
расстояние. Во время симуляции сфера с некоторой постоянной скоростью
смещается в нулевые координаты. При достижении нулевых координат
симуляция считается завершенной.
Ниже приведены примеры симуляции с плоскостями разной плотности сетки.
Трехмерная модель плоскости представляет собой квадрат со стороной 1 см.
Выделена область в центре модели плоскости. По области построены модель
масс
и
пружин,
иерархия
ограничивающих
объемов.
Область
продублирована и смещена, между получившимися слоями добавлены
26
пружины. Модель плоскости состоит из 4225 вершин (на картинках —
слева), 16641 вершины (на картинках — посередине), 66049 вершин (на
картинках — справа). При меньшем числе вершин (несколько сотен)
значения площадей и объемов после симуляции значительно расходятся со
значениями при выбранных числах вершин. Работа библиотеки со
значительно большим числом вершин не предполагается. Трехмерная модель
сферы радиуса 0.25 см состоит из 994 вершин, 1984 треугольников. По
модели сферы построена модель масс и пружин.
Размеры, площадь и толщина плоскостей во всех примерах одинаковы.
Общая масса модели плоскости равняется произведению площади на
толщину. Массы частиц вычислены пропорционально площадям содержащих
их треугольников. Константные значения для жесткости пружин одинаковые
во всех примерах.
На Рис. 4 показана сетка плоскости, выделенной области и построенного
второго слоя до проведения симуляции, на Рис. 5 — пружины области.
Рис. 4: Сетка до симуляции
Рис. 5: Пружины до симуляции
27
На Рис. 6 показана сетка выделенной области, сетка построенного второго
слоя и сетка плоскости после симуляции. На Рис. 7 показаны вершины сетки.
Здесь видно, что для плоскостей разной плотности деформация различается
незначительно.
Рис. 6: Сетка после симуляции
Рис. 7: Вершины после симуляции
Ниже представлена таблица (Таблица 1) результатов симуляции для
плоскостей,
показанных
на
картинках
выше.
Во
второй
колонке
представлены результаты для плоскости, показанной на картинках слева, в
третьей результаты для плоскости на картинках посередине, в четвертой —
соответственно для плоскости на картинках справа. Единицы измерения
площадей — см2, объемов — см3, время симуляции — в секундах. Площади и
объемы выделенных областей после симуляции отличаются незначительно.
Площади выделенных областей до симуляции отличаются в том числе из-за
размера треугольников при разной плотности сетки. Время выполнения
симуляции (при одинаковом числе шагов симуляции) с увеличением
28
плотности сетки значительно увеличивается.
Число вершин плоскости
4225
16641
66049
Число треугольников плоскости
8192
32768
131072
Число вершин выделенной области
3922
15090
59154
Число треугольников выделенной области 7504
29512
116984
Число пружин
22176
87872
349632
Площадь области до симуляции, см2
0.458
0.450
0.446
Площадь области после симуляции, см2
0.672
0.689
0.698
Объем области после симуляции, см3
0.062
0.063
0.064
Время симуляции, секунды
11
31
97
Таблица 1. Результаты симуляции плоскости разной плотности сетки
Обнаружение столкновений с помощью иерархии ограничивающих объемов
позволяет производить симуляцию с указанными плотностями сеток в
реальном времени.
Корневой узел иерархии ограничивающих объемов — узел первого уровня.
На Рис. 8 показаны узлы 7-го, 8-го, 9-го уровней иерархии ограничивающих
объемов, построенной по выделенной области. Число узлов 7-го уровня —
64, 8-го — 128, 9-го — 256, общее число узлов — 7503. Узлы иерархии
полностью покрывают собой область. На картинках узлы показаны ребрами
параллелепипедов.
Рис. 8: Узлы 7-го, 8-го, 9-го уровней иерархии ограничивающих объемов
29
Ниже (Рис. 9) показан пример сохранения объема. Модель сферы
продавливает модель куба сверху. Модель куба со стороной 2 см, объема 8
см3 представляет собой замкнутую трехмерную модель, состоит 386 вершин,
768
треугольников.
Модель
сферы
состоит
из
134
вершин,
264
треугольников, объем — 3.9 см3. Построены модели масс и пружин по
трехмерным моделям куба и сферы. Для модели масс и пружин куба
построены пружины для сохранения исходного объема. При деформации куб
сохраняет исходный объем, на картинке посередине и справа объем
соответственно 7.97 и 7.96 см3.
Рис. 9: Сохранение объема
30
7 Модуль помещения имплантата
7.1 Описание модуля
В
модуле
реализован
графический
интерфейс
для
ввода
данных
пользователем и логика обработки этих данных. Исходные данные для
проведения симуляции — трехмерная модель участка тела пациента,
положение и параметры имплантата. Пользователем через графический
интерфейс задаются область для помещения имплантата, его форма, размеры,
глубина помещения. Область определяется вокруг имплантата при выборе
его положения на поверхности исходной трехмерной модели. Затем
пользователем запускается симуляция. По ее завершении отображаются
область после симуляции, разницы площадей и объемов области до и после
симуляции.
Взаимодействие с пользователем в модуле разделено на шаги. На первом
шаге пользователь выбирает модель пациента, на втором задает параметры и
положение имплантата на модели пациента, при переходе на третий шаг
выполняется симуляция, на третьем шаге отображаются результаты
симуляции.
Разница объемов между областями до и после проведения симуляции
вычисляется по полигональной модели, полученной объединением обеих
областей, т.к. выделяемая для симуляции область является незамкнутой.
Класс модуля, реализующий взаимодействие с библиотекой физических
расчетов, содержит две модели масс и пружин, построенные соответственно
по модели области и по модели имплантата, иерархию ограничивающих
объемов для обнаружения столкновений, построенную по полигональной
модели области. В том же классе осуществляются вызовы методов
обновления состояний систем масс и пружин, методов обновления иерархии
ограничивающих объемов, методов обработки столкновений.
31
7.2 Симуляция
При инициализации моделей масс и пружин для вычисления масс частиц
общая масса выделенной области вычисляется через ее площадь, в качестве
толщины задается толщина кожи. Общая масса имплантата вычисляется
через его объем. Плотность задается плотностью воды. Другие константные
значения в формулах подобраны для плотности сетки моделей пациентов
системы планирования.
Перед запуском симуляции строятся два слоя пружин по выделенной
области: один по исходной области, второй — копированием и смещением
первого слоя на глубину помещения имплантата. Между слоями для
сохранения толщины строится множество пружин. Частицам нижнего слоя
добавляются пружины нулевой длины, возвращающие частицы в исходные
положения. Для обработки столкновений между моделью масс и пружин
выделенной области и моделью имплантата иерархия строится по нижнему
слою пружин области.
Частицы слоев на границе заданной области считаются неподвижными,
остальные — подвижными. Перед симуляцией все подвижные частицы слоев
смещаются на высоту имплантата. С началом симуляции частицы слоев
начинают смещаться в исходные положения. Во время симуляции
столкновения обрабатываются между нижним слоем пружин и моделью
имплантата. Также во время симуляции деформируется модель имплантата,
при деформации вычисляются силы для сохранения объема. Симуляция
останавливается через какое-то количество шагов по достижении некоторого
минимума суммарного перемещения частиц.
7.3 Примеры
Помещение модели имплантата разного объема под модель мягких тканей.
32
Модель плоскости представляет собой квадрат со стороной 20 см. Модель
имплантата в основании представляет собой круг диаметра 7 см. Модель
плоскости состоит из 4225 вершин, число вершин выделенной области —
2016, число пружин — 11140. Площадь области до симуляции составляет
92.57 см2. Модель имплантата состоит из 514 вершин. Модель имплантата
объемом 50, 100, 150 см3, на картинках соответственно слева, посередине,
справа.
На Рис. 10, Рис. 11 показаны выделенная область, модель имплантата до
проведения симуляции.
Рис. 10: Выделенная область до симуляции и модель имплантата
Рис. 11: Сетка выделенной области до симуляции и модели имплантата
На Рис. 12 показана выделенная область после симуляции. На Рис. 13
показана сетка выделенной области.
33
Рис. 12: Выделенная область после симуляции
Рис. 13: Сетка выделенной области после симуляции
Ниже представлена таблица (Таблица 2) результатов симуляции помещения
модели имплантата разного объема под модель мягких тканей. Единицы
измерения площадей — см2, объемов — см3, время симуляции — в секундах.
При увеличении объема увеличивается натяжение пружин системы масс и
пружин, стабилизация происходит дольше.
Объем имплантата, см3
50
100
150
Площадь области после симуляции, см2
104.06
124
149.73
Разница площадей, см2
11.49
31.43
57.16
Разница объемов, см3
53.14
104.74
151.79
Число шагов симуляции
2800
2900
3300
Время симуляции, секунды
Таблица 2. Результаты симуляции
23
27
36
34
Ниже приведен пример работы модуля в системе планирования. Имплантат
объема 100 см3 представляет собой круг диаметра 7 см в основании. Площадь
выделенной области — 91 см2, площадь после симуляции — 135 см 2, разница
площадей — 44 см2, полученный объем при растяжении — 109 см3, время
симуляции — 23 секунды. На Рис. 14 показано задание параметров для
проведена симуляции.
Рис. 14: Задание параметров для симуляции
Ниже показана модель после симуляции. На Рис. 15 показана сетка модели,
отмечена выделенная область, на Рис. 16 показана полученная модель, сетка
выделенной области.
35
Рис. 15: Сетка полученной модели после симуляции
Рис. 16: Сетка выделенной области после симуляции
36
8 Результаты
Разработан модуль, позволяющий симулировать возможный результат
помещения имплантата под мягкие ткани по трехмерной модели поверхности
тела пациента.
Достигнуты следующие результаты:
1. определена архитектура модуля
2. разработана библиотека физических расчетов
3. разработан модуль помещения имплантата
37
9 Список литературы
[1] Цыдик И. С. и др. ПЛАСТИКА ДЕФЕКТОВ МЯГКИХ ТКАНЕЙ (ОБЗОР
ЛИТЕРАТУРЫ) //Журнал Гродненского государственного медицинского
университета. – 2006. – №. 3 (15).
[2] Пасичный Д. А. Дермотензия в лечении повреждений покровных тканей
стопы и голени //Международный медицинский журнал. – 2009.
[3] Chen D.T. Modeling for Plastic and Reconstructive Breast Surgery / D.T.
Chen, I.A. Kakadaris, M.J. Miller, R.B. Loftin, C. Patrick // Medical Image
Computing and Computer-Assisted Intervention – MICCAI 2000. Lecture Notes in
Computer Science. – Berlin; Heidelberg: Springer, 2000. – Vol. 1935. – P. 10401050.
[4] Zollner A. M. et al. Growth on demand: reviewing the mechanobiology of
stretched skin //Journal of the mechanical behavior of biomedical materials. –
2013. – Т. 28. – С. 495-509.
[5] Николаев С. Н. Нелинейная система масс-с-пружинками для
моделирования больших деформаций мягких тканей, 2013
[6] Nealen A. et al. Physically based deformable models in computer graphics
//Computer graphics forum. – Blackwell Publishing Ltd, 2006. – Т. 25. – №. 4. –
С. 809-836.
[7] Lloyd B. A., Székely G., Harders M. Identification of spring parameters for
deformable object simulation //Visualization and Computer Graphics, IEEE
Transactions on. – 2007. – Т. 13. – №. 5. – С. 1081-1094.
[8] Müller M. et al. Position based dynamics //Journal of Visual Communication
and Image Representation. – 2007. – Т. 18. – №. 2. – С. 109-118.
[9] Baudet V. et al. Integrating tensile parameters in mass-spring system for
deformable object simulation //Rapport de recherché, Technical Report. – 2009
[10] Bianchi G. et al. Simultaneous topology and stiffness identification for massspring models based on fem reference deformations //Medical Image Computing
38
and Computer-Assisted Intervention–MICCAI 2004. – Springer Berlin Heidelberg,
2004. – С. 293-301.
[11] Gelder A. V. Approximate simulation of elastic membranes by triangulated
spring meshes //Journal of graphics tools. – 1998. – Т. 3. – №. 2. – С. 21-41.
[12] Teschner M. et al. A versatile and robust model for geometrically complex
deformable solids //Computer Graphics International, 2004. Proceedings. – IEEE,
2004. – С. 312-319.
[13] Vassilev I. T., Spanlang B. A mass-spring model for real time deformable
solids, Proceedings of East-West Vision 2002, pp. 149-154, 2002
[14] Baraff D., Witkin A. Large steps in cloth simulation //Proceedings of the 25th
annual conference on Computer graphics and interactive techniques. – ACM, 1998.
– С. 43-54.
[15] Ericson C. Real-time collision detection. – CRC Press, 2004.
[16] Собинов Д. И., Коробицын В. В. Алгоритмы обнаружения
столкновений //Математические структуры и моделирование. – 2010. – №. 21.
– С. 82-95
[17] Faure F., Allard J., Nesme M. Eulerian Contact for Versatile Collision
Processing : дис. – INRIA, 2007.
[18] Mirtich B., Canny J. Impulse-based dynamic simulation. – California :
Computer Science Division (EECS), University of California, 1994.
[19] Canfield, www.canfieldsci.com
[20] ANSYS, www.ansys.com
[21] ABAQUS, www.3ds.com/products-services/simulia/products/abaqus/
[22] SOFA, www.sofa-framework.org
[23] FEBio, www.febio.org
[24] Zachmann G., Langetepe E. Geometric data structures for computer graphics.
– Eurographics Assoc., 2003.
39
Отзывы:
Авторизуйтесь, чтобы оставить отзыв