CC BY
5УДК 004.932
Реализация иерархической модели данных в системе компьютерного планирования хирургических операций в ортопедии
© Гончаренко В.Г., Архипов В.И., Тузиков A.B.
Объединенный институт проблем информатики HAH Беларуси ул. Сурганова, 6, г. Минск, 220012, Беларусь
e-mail: vasUy@mpen.bas-Bet.by, arkhipau@gmaU.com, tuzikov@newman.bas-net.by
Abstract. Hierarchical model describing transformation results of large images is discussed. A procedure of binary image resection by a polygon ("scissor") selected by a user is described. All connected components, completely crossed by the polygon, are resected onto several parts.
Введение
В 1972 году в связи с появлением компьютерной томографии появилась возможность создавать трехмерные изображения внутренних органов и тканей. Эти изображения строятся в виде последовательности двухмерных слоев, составляющих трехмерное объемное изображение. Каждый слой томографического изображения независим от других, при этом ярость каждой точки на томографическом изображении характеризует плотность соответствующей анатомической области, что позволяет по яркости точки с определенной достоверностью классифицировать тип ткани (костная, мышечная, жировая и т.д.). Все это позволяет использовать трехмерные томографические изображения для построения объемных моделей внутренних органов и применять их для планирования хирургических операций при помощи компьютера [1, 2].
Планирование операций выполняется предварительно для вычисления параметров основных процедур хирургической операции в целях подготовки к реальной операции. Объектами операции являются органы и ткани пациента, имеющие выраженную патологию и требующие хирургического вмешательства. Обычно последовательность действий во время хирургической операции определяется хирургом самостоятельно на основании изучения состояния анатомических объектов по результатам предоперационного планирования [3].
Компьютерное планирование выполняется при помощи специальных программных средств путем воспроизведения на компьютере всей последовательности действий, составляющих хирургическую операцию. Компьютерное планирование операции можно назвать виртуальной операцией. Виртуальная операция не несет никакого риска для пациента и, в то же время, позволяет легко определить действия и их параметры, необходимые для проведения реальной операции. Для определения геометрических параметров костей и суставов анализируются цифровые изображения.
Процесс компьютерного планирования хирургических операций в ортопедии обычно состоит из 8 этапов:
1) сканирование пациента с помощью компьютерного томографа;
2) предварительная обработка томографического изображения и построение иерархической модели;
3) сегментация и классификация объектов планирования;
4) оценка визуальной информации;
5) измерение моделей костей;
6) оценка результатов измерений и принятие решений относительно дальнейших действий;
7) виртуальная остеотомия;
8) виртуальный остеосинтез.
Подробное описание этих этапов приводится в работе [3]. Здесь рассматривается понятие иерархической модели, а также задача разрезания объектов иерархической модели в процессе виртуальной остеотомии.
1. Концепция иерархической модели
В целях организации работы с различными типами преобразований отдельных частей изображений была спроектирована так называемая «иерархическая модель». В связи с тем, что различные части изображения могут потребовать различных методов их преобразования, предложено использование одновременно нескольких различных способов выделения объектов на изображении.
Иерархическая модель представляет собой дерево, вершинами которого являются объекты, получившиеся в результате преобразований различных участков исходного изображения. В корне дерева находится исходное изображение. В дереве выделяется один из объектов, который считается текущим - с ним пользователь работает в данный момент. Так, например, сегментация изображения всегда применяется к части текущего объекта, определяемой областью интереса в виде прямоугольного параллелепипеда. Объекты иерархической модели с точки зрения программной реализации представляют собой объекты класса С++ КРОШегагсЫсаЮЪ^есЬ. Это абстрактный класс, имеющий ряд чистых виртуальных функций, таких как GetObjectVolume и Такая концепция позволяет реализовать два типа операций: хорошо параметризуемые и плохо параметризуемые. Объекты, являющиеся результатом хорошо параметризуемой операции, можно реконструировать за приемлемое время, используя в качестве параметров процедуры реконструкции некоторые атрибуты, требующие мало памяти по сравнению с самим объектом. Таким образом, объекты, получившиеся в результате хорошо параметризуемой операции, не хранят в явном виде значения вокселов соответствующего им участка. Вместо этого они хранят атрибуты, которые позволяют строить массив вокселов, используя данные их родителей и исходного изображения. Так, например, пороговую сегментацию можно трактовать как хорошо параметризуемую операцию, атрибутами которой являются пороговые уровни. Значения вокселов (бинарного изображения) можно построить из исходного
изображения, используя атрибуты (значения пороговых уровней). Такая схема позволяет освободить память, выполняя построение массива вокселов только тогда, когда он необходим, после чего он удаляется из оперативной памяти. Объекты, являющиеся результатами плохо параметризуемой операции, используются, когда построение массива вокселов объекта нельзя произвести за приемлемое время, или когда атрибуты занимают больше памяти, чем сам результат операции, В качестве примера плохо параметризуемой операции может выступать ручная сегментация.
Каждый объект иерархической модели имеет свои размеры, а также координаты и матрицу ориентации относительно своего объекта-родителя. Это позволяет быстро выполнять перенос и вращения всей ветви дерева иерархической модели.
Объектами иерархической модели могут быть не только результаты сегментации изображения, но и результаты промежуточных преобразований (например, фильтрации), Предложенная концепция хорошо и плохо параметризуемых операций позволяет учесть множество различных преобразований участков изображения.
После загрузки исходного изображения оно помещается в корень иерархической модели. Как только пользователь выполнит сегментацию какой-либо части исходного изображения, отсегментированное изображение помещается в иерархическую модель как потомок корневого элемента (т.е. исходного изображения). Обычно после сегментации выполняется классификация отсегментированных объектов - они помещаются на третий уровень иерархической модели (рисунок 1), Количество классифицированных объектов соответствует количеству уникальных объектов, имеющихся на классифицируемом участке текущего объекта. Пользователь в любое время может выбрать текущий объект и те объекты, которые необходимо отображать (текущий объект всегда является отображаемым объектом). Если отображается какой-либо объект, то в это же время его потомки и предки отображаться не могут. Все действия по выбору объектов пользователь выполняет при помощи менеджера объектов - специального диалогового окна, где отображена древовидная структура текущей иерархической модели.
Исходное изображение
Сегментированный объект
Классифицированный объект
Классифицированный объект
Классифицированный объект
Сегментированный объект
Классифицированный объект
Рис. 1. Примерная структура обычной иерархической модели
2, Выделение связных компонент и разрезание объектов
Одной из важных плохо параметризуемых операций над бинарными объектами иерархической модели в системе компьютерного планирования выступает разрезание объектов. Эта операция имеет большое значение при планировании остеотомии [2].
В процессе выполнения виртуальной операции появляется необходимость разрезать текущий объект на две или несколько частей замкнутой ломаной линией р = (Б|.....Б/,). образующей /.'-угольник, все вершины которого лежат в одной
плоскости Щ/>). Эта ломаная образует активную часть плоскости разрезания, называемую «резцом» (рисунок 2). Бинарный объект состоит из некоторого множества связных компонент. При использовании резца разрезаются только те связные компоненты, которые полностью пересекаются активной частью плоскости разреза. Сама процедура разрезания использует операцию выделения связных компонент и является, по сути, её модификацией.
Рис. 2. Разрезание объекта инструментом «резец», в котором активная
часть плоскости разрезания задаётся шестиугольником (>',.....). и
разрезаются только объекты, попадающие в активную часть плоскости разреза: а) вид объектов до разрезания; б) вид объектов после разрезания
Рассмотрим, как определяются связи между соседними вокселами трехмерного изображения. Если инструмент «резец» не задан, либо он не пересекается с текущим объектом, то все вокселы, принадлежащие объекту и являющиеся соседями в используемой метрике, имеют связи между собой (рисунок 3(а)), Как только задается «резец», пересекающий текущий объект, вычисляются точки пересечения активной части плоскости резца р и связей между вокселами, образующими объект (см. рисунок 3(в)), «Резец» разрезает только ребра решетки, пересекаемые активной частью плоскости разрезания р.
Опишем более подробно процесс построения связей элементов изображения при использовании инструмента «резец». Пусть задано бинарное изображение В, на котором объекты имеют значение 1, а фон - значение 0. Обозначим через п вектор
® ® ® ® ® ® ® ® ф ® © ® ф ®
® (Ь<р ® ® ф ® ® ® ф ®
® (Мр ® ®
® (Ь<Ъ ® ® ® ® ® ® ® ® ®
(а) (б)
® ® ® ® ® ® ® ® ф ® © ® ф ®
® (Ь<р ® ® © ® ® ® ф ®
® ф ®х®,.ср4 ®
® сЬ<Ь ® (м ® ® ® ® ® ® ® ®
(в) (г)
Рис. 3. Пример разрыва связей пикселов связной компоненты инструментом «резец» с использованием четырехсвязной метрики; а) построение связей между пикселами при выделении связных компонент; б) связные компоненты изображения, помеченные соответствующими номерами; в) разрыв связей между пикселами (пунктирная линия); г) связные компоненты, получившиеся при разрезании объекта,
нормали плоскости разрезания. Введем функцию Л<>(р) положения воксела р относительно плоскости разрезания Щр):
5п(р) = п • (Эх - р), (1)
где «•» обозначает скалярное произведение двух векторов.
Если воксел р принадлежит объекту и у него имеется сосед с|. принадлежащий тому же объекту, то можно определить, пересекает ли плоскость Щ//) связь между р и с|. Связь пересекается, если значение функции Л<> для вокселов р п с| имеет разный знак:
< 0. (2)
Связь между вокселом р и его соседом с| пересекается «резцом» р и разрывается, если она пересекается плоскостью Щр), и точка 1Р)Ч пересечения прямой, проходящей через р п (1- п плоскости разрезания Щ/>) находится внутри многоугольника р.
0 0 0 0 0 0 0
0 1 0 2 0 3 0
0 1 1 0 0 3 0
1 1 1 0 0 3 0
0 1 1 0 3 3 0
0 1 1 0 3 3 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 1 0 2 0 3 0
0 1 1 0 0 3 0
4 4 1 0 0 3 0
0 4 4 0 3 3 0
0 4 4 0 3 3 0
0 0 0 0 0 0 0
Координаты точки 1Р,, вычисляются по формуле:
Ip,q = P+ q-P -7-ч- 3
n • (q - p)
В случае, если знаменатель в (3) равен 0 (прямая, соединяющая р и q, параллельна плоскости разрезания Щ/>) ). связь между р и q считается неразорванной.
Перейдем к двухмерным координатам и обозначим через // = (S'i,..., S'/,) двухмерный многоугольник, координаты вершин которого определяются следующим образом:
с>х _ с ;
% 2 (Л\
с 'У — с ; V /
° i — j 7
где i изменяется от 1 до k, a i.j базисные векторы, определяющиеся следующим образом:
ж ~~ ||s2-si| j = i X п.
Обозначим через I' = (1р.,г 1р.,,) точку 1Р<1 в базисе (i.j):
p,q' рл-1
{ = ^р.Ч ' J-
Tlx т
ap,q р>ч
Рассмотрим задачу определения принадлежности точки многоугольнику. Существует несколько различных способов решения этой задачи [4]. В данном случае, когда количество сторон многоугольника невелико, реализовано следующее решение. Вначале определяется некоторый базовый луч, исходящий из точки 1'р . Затем определяется количество ребер многоугольника // пересекаемых этим лучом. Если количество пересечений нечетное, то точка () лежит внутри многоугольника // и, следовательно, связь между р и с| разрывается. Иначе связь между р и с| сохраняется.
В случае, если многоугольник // является выпуклым, имеется более простое решение данной задачи. Вообще говоря, в большинстве случаев резец вполне можно сделать выпуклым многоугольником, не теряя при этом универсальности его применения и выигрывая в удобстве его использования. Достаточно удобно пользоваться резцом, представляющим собой выпуклый 4-угольник, Точка 1'р () лежит внутри выпуклого многоугольника //. если точка находится всегда с одной стороны по отношению ко всем сторонам многоугольника при обходе их в одном направлении. Другими словами, точка () лежит внутри выпуклого многоугольника //. если величина //,,.,,(/. / + I) имеет одинаковый знак для всех г = 1.2.....к:
г + 1) = (1р)Ч — в/ХЭ^ — ) — (1р)Ч — )(БД1 — Б/), (5)
здесь [лрл(к, к + 1) = ЦрЛ(к, 1).
При реализации данной схемы вычислений необходимо учесть, что для каждого помещаемого в очередь воксела р целесообразно заранее вычислять и сохранять значение 5п(р), Его можно использовать позже для проверки, пересекает ли плоскость Щ/>) связь между выбранным из очереди вокселом р' и его соседом с|. Вычисленное значение Л<>(с|) также желательно сохранять вместе с помещаемым в очередь вокселом с|.
Для случая регулярной метрики необходимо принимать во внимание расположение соседних вокселов. Например, при использовании шестисвязной метрики для трехмерного изображения каждый воксел содержит по 6 соседей (кроме вокселов, расположенных на границе изображения). Расстояние между текущим вокселом и его соседями равно 1, Причем, положение каждого соседа отличается от положения текущего воксела только одной координатой. Соответственно, выражения для Л<)((|) и 1р,, можно вычислять для каждого соседа в зависимости от его положения (таблица 1), Здесь р<. р« р- - координаты вектора р в трехмерном пространстве. Аналогично обозначаются координаты вектора п.
Таблица 1. Оптимизация вычислений в случае шестисвязной метрики на трехмерной регулярной решетке
Координаты соседа q воксела р Значение <5n(q) Значение 1р<1
(Рг - 1.р".р-) <5п(р) + п* (рж + <^п(р)/пж,ру,рг)
(Р'г + 1. р". р") <5п(р) - п* (рж + ^п(р)/пж,ру,рг)
(р'г. р" — 1, р') 8п( Р) + к" (px,pv + 5п(р)/п2',рг)
(Р'г. р" + 1 - Г> ) <fo(p) - п^ (рх,ру + <5п(р)/пг',рг)
(р'г. р". р' — 1) 5п(р) +nz (р'.р^.р- +5п(р)/пг)
(Р'г. р". pz + 1) 5п(р) - п- (р'.р^.р- +5п(р)/пг)
Рассмотрим подробнее алгоритм выделения связных компонент и разрезания объектов изображения (алгоритм 1), Он работает для изображения с любой связностью. Пусть каждый воксел имеет флаг состояния просмотра (TEMP - еще не просматривался, DONE - уже просматривался и включен в состав связной компоненты), Вначале инициализируется номер текущей связной компоненты с •<= 1, Затем алгоритм переходит в стадию последовательного просмотра массива бинарных вокселов. Как только обнаружится непросмотренный (имеющий флаг TEMP) воксел р, принадлежащий объекту (/v//(p) = 1), он помещается в хвост очереди при помощи функции Enqueue(р), присоединяется к текущей связной компоненте ConnectedComponent(р) <= с и ему присваивается флаг DONE. После этого номер текущей связной компоненты увеличивается на 1 и алгоритм переходит в стадию наращивания связной компоненты, которой принадлежит только что просмотренный воксел р. При этом, из головы очереди снимается один воксел р' <= DequeueQ. Затем рассматриваются все его соседи, с которыми он имеет связи. Если воксел q имеет
флаг TEMP, принадлежит объекту и находится с обратной стороны от плоскости разрезания Щр) относительно воксела р', то выполняется проверка на то, попадает ли точка пересечения связи вокселов р' и q с плоскостью Щ//) внутрь многоугольника р. Если связь между вокселами р' и q пересекается «резцом», то она удаляется. Иначе воксел q добавляется в хвост очереди командой Enqueue^q), присоединяется к связной компоненте С onnectedC omponent(р') и ему присваивается флаг DONE. Такая проверка выполняется для всех соседей воксела р', с которыми у него есть связи. После того, как очередь опустеет, алгоритм опять переходит в стадию последовательного просмотра вокселов и переходит к следующему за р вокселу. Процесс прекращается как только будут просмотрены все вокселы объекта, В результате этой процедуры создется изображение, состоящее из вокселов со значениями индексов соответствующих им связных компонент (рисунок 3(6)), Каждая связная компонента с одинаковым значением вокселов будет образовывать отдельный объект.
Алгоритм 1 Выделение связных компонент и разрезание бинарных объектов
На входе: Бинарное изображение В - массив вокселов р G В с множеством значений val(p) G {0,1}; множество R связей п / (p. q) С />' между вокселами p. q С В На выходе: Бинарное изображение ConnectedComponent, состоящее из индексов связных компонент, получившихся после разрезания 1: for all р G В do 2: flag(p)<=TEMP-3: с 4= 1;
{Просмотр вокселов} 4: for all р G В do
5: if flag(р) = TEMP and val(p) = 1 then
6: Enqueue( p); 7: С onnectedC omponent(p) -<= c; 8: flag(p) 4= DONE; 9: c<=c + l;
{Наращивание связной компоненты} 10: while Очередь не пуста do 11: р' -<= DequeueQ;
12: for all сосед q воксела p' do
13: if flag(q) = TEMP and /v//(q) = 1 then
14: connection <= TRUE;
{Проверка на пересечение плоскостью Щр) связи р' с q} 15: if <fo(p')<fo(q) < 0 then
16: if знак //,,'.,,(/. i + 1) одинаков для всех i = 1,..., k then
17: connection <= FALSE;
18: if connection then
19: Enqueue^ q);
20: С onnectedC omponent( q) <= С onnectedC omponent(p');
21: flag(q) DONE;
На рисунке 4 показан пример разрезания бедренной кости.
(а) (б)
Рис.. 4. Пример работы алгоритма разрезания: а) объект до разрезания: б) объект после разрезания.
Заключение
В статье разработана методика иерархического моделирования, позволяющая эффективно сохранять результаты различных преобразований частей изображения. Все преобразования делятся на два тина: хорошо и плохо параметризуемые. Результаты хорошо параметризуемых операций упаковываются таким образом, что в явном виде хранятся только их параметры, а результат строится «на ходу». Эффективность такой методики основывается на том, что параметры хорошо параметризуемой операции занимают значительно меньший объем памяти, чем ее результаты. Иерархическая модель представляет собой дерево, элементами которого выступают результаты преобразований частей одного изображения, В корне этого дерева находится исходное изображение. Каждое последующее преобразование части изображения создает ветвь дерева. Это позволяет связывать все операции, выполняемые над данными одного изображения.
Разработан алгоритм выделения связных компонент и последующего разрезания бинарных объектов специальным инструментом «резец». Этот инструмент представляет собой многоугольник, все точки которого лежат в одной плоскости, называемый активной частью плоскости разрезания. Разрезаются только те связные компоненты
объекта, которые пересекаются активной частью плоскости разрезания. Предложенный инструмент «резец» позволяет выполнять разрезание бинарных объектов любой формы, включая отдельные кости таза,
список литературы
1. Handels Н., Ehrhardt J., Plötz W., Pöppl S.J. Three-dimensional planning and simulation of hip operations and computer-assisted construction of endoprostheses in bone tumor surgery. Computer Aided Surgery, vol. 6, Nu. 2, p. 65-76, 2007.
2. Krivonos 0. Hesser J. Maenner R. Keppler P. Gebhard F. Kinzl L. Sakalouski A.A. Sakalouski A.M. Precise computer aided correction of bone deformities. Proceedings of the World Congress on Medical Physics and Biomedical Engineering WC2003, Sydney, Australia, August 2003.
3. Hancharenka V.G., Tuzikov A.V. Arkhipau V.l., Kryvanos A. Preoperative planning of pelvic and lower limbs surgery by CT image processing. Proceedings of 8th International Conference on Pattern Recognition and Image Analysis PRIA'2007, 7-13 Oct., Yoshkar-Ola, Russia, vol. 1, p. 270-274, 2007.
4. Haines E. Point in Polygon Strategies. Graphics Gems IV, ed. Paul Heckbert, Academic Press, p. 24-46, 1994.
Статья поступила в редакцию 19.04-2008
