Компьютерная технология восстановления пространственной структуры коронарных сосудов по ангиографическим проекциям Текст научной статьи по специальности «Компьютерные и информационные науки»

ОБРАБОТКА ИЗОБРАЖЕНИЙ, РАСПОЗНАВАНИЕ ОБРАЗОВ

КОМПЬЮТЕРНАЯ ТЕХНОЛОГИЯ ВОССТАНОВЛЕНИЯ ПРОСТРАНСТВЕННОЙ СТРУКТУРЫ КОРОНАРНЫХ СОСУДОВ ПО АНГИОГРАФИЧЕСКИМ ПРОЕКЦИЯМ

Наталья Юрьевна Ильясова1'2 (старший научный сотрудник, e-mail: [email protected]), Николай Львович Казанский1'2(заместитель директора, e-mail: [email protected]), Андрей Олегович Корепанов1 (доцент, e-mail: [email protected]), Александр Викторович Куприянов12 (старший научный сотрудник, e-mail: [email protected]),

Андрей Владимирович Устинов2 (ведущий программист, e-mail: [email protected]), Александр Григорьевич Храмов12 (старший научный сотрудник, e-mail: [email protected]) 1 Самарский государственный аэрокосмический университет имени академика С.П. Королева, 2 Учреждение Российской академии наук Институт систем обработки изображений РАН

Аннотация

Представлена компьютерная технология автоматического восстановления пространственной структуры коронарных артерий по малому числу ангиографических проекций. Технология построения пространственной структуры базируется на пространственной трассировке сосудов на основе одновременного анализа всех доступных проекций. В результате формируется векторное представление пространственных сосудов, которое в дальнейшем используется для визуализа -ции структуры. Подробно представлены этапы восстановления, приводятся результаты работы отдельных этапов. Описывается один из методов получения пространственных координат ска -нирования трассы по результатам обработки двух проекций изображения. Предлагается набор основных геометрических характеристик, описывающих объекты двухмерных и трёхмерных структур. Описаны аппроксимационные методы оценивания диаметра сосудов, методы выделения центральных линий сосудов на изображениях проекций на основе вейвлет преобразования. Представлены экспериментальные исследования на модельных и натурных объектах.

Ключевые слова: коронарные сосуды, восстановление пространственной структуры, геометрические характеристики сосудов, ангиографическая проекция.

Введение

Рентгеновская ангиография используется в кардиологии для диагностирования заболеваний коронарных артерий и определения степени поражения сосудов [1]. Обычно диагностирование производится посредством визуального наблюдения ангио-грамм и носит субъективный характер, т.к. на результаты диагностики существенное влияние оказывает ракурс зарегистрированных проекций. Пространственная модель сосудов сердца служит для визуализации трехмерной геометрической и топологической информации, позволяет увеличить точность измерений и получить численные оценки геометрических характеристик сосудов.

Исходными данными являются последовательности кадров ангиографической съемки сердца в формате БЮОМ, широко распространенном в медицинской практике [2, 3]. На процесс получения проекций накладывается ряд ограничений технического характера, таких как предельные углы поворота камеры и максимальная доза облучения пациента. В связи с этим становится невозможным получить большое количество проекций, что существенно усложняет процесс восстановления. Помимо этого, съемка проекций производится не одновременно и вследствие того, что сердце является динамическим объектом, проекции являются рассогласованными (т.е. представляют собой в разные моменты времени изображения разных объектов). Эти обстоятельства делают невозможным использование методов компьютерной томографии

для восстановления пространственной интенсивности и точное восстановление пространственной структуры в какой-либо определенный момент времени. Следует также отметить еще одно обстоятельство, существенно затрудняющее процесс восстановления. Регистрация проекций производится в расходящемся пучке рентгеновского излучения, вследствие этого результирующее изображение наблюдается в центральной проекции. Восстановление исходного объекта по таким проекциям невозможно без наличия дополнительной информации о взаимном расположении камеры в момент регистрации проекции и тела пациента.

Рис. 1. Задача восстановления пространственной структуры коронарных сосудов

В этих условиях задачу можно сформулировать следующим образом. По имеющимся сильно зашум-лённым центральным проекциями в разные моменты времени некоторого динамического объекта (дерева

сосудов) при отсутствии точной информации о геометрии съемки требуется восстановить пространственную структуру объекта, наиболее близко по некоторому критерию соответствующего данным проекциям (рис.1). Ниже представлена технология автоматического восстановления пространственной структуры коронарных артерий по малому числу ангио-графических проекций. В настоящее время существует множество различных подходов к восстановлению [4, 9-11], которые, по сути, различаются методом выделения двумерной структуры сосудов на изображениях, способом конструирования пространственной структуры из нескольких двумерных структур с известной геометрией, а также степенью автоматиза-

ции. Следует отметить, что здесь представлены основные ключевые моменты методики и опущены некоторые подробности восстановления (такие как анализ стеноз, выбор пространственного направления по сосуду, логический анализ ответвлений, восстановление 2Б дерева, некоторые геометрические преобразования, алгоритм градиентной тесселляции для построения поверхности сосуда на этапе визуализации и т.п.) [12-17]. Цель настоящей работы - продемонстрировать подход к проблеме восстановления пространственной структуры коронарных артерий, позволяющий максимально автоматизировать процесс восстановления пространственной структуры по малому числу проекций.

Рис. 2. Примеры изображений, характеризуемых наличием ветвей: а) кровеносная система глазного дна человека - древовидная структура; б) кровеносная система сердца - древовидная структура; в) сеть трещин в земле вдоль рек - сетчатая структура

пространственной структуры сосуда были введены некоторые понятия, представленные на рисунке 3.

Модель древовидного объекта Рассматриваемые изображения ангиографических проекций принадлежат широкому классу изображений, характеризуемых наличием древовидных структур. К таким изображениям можно отнести изображения сетчатки глаза человека, русла рек, сети дорог, ан-гиографические снимки сосудов кровеносной системы человека и другие биологические структуры (рис .2). Примером пространственных древовидных объектов могут служить ветви деревьев, сосуды кровеносной системы человека. Все описанные объекты имеют схожую структуру, топологию которой легко можно описать бинарным односвязным графом, в котором корневой элемент соответствует началу объекта, листья являются крайними точками объекта, а все промежуточные элементы графа характеризуют точки ветвления. Будем называть узлом часть сосуда, в которой происходит ветвление. Эта область имеет специфику, так как она принадлежит одновременно нескольким ветвям сосуда, являясь их общей частью. Узел характеризуется своей протяженностью в проекциях, а также углом ветвления. Заметим, что помимо узловых точек особыми являются также конечные точки ветвей, или просто концы, и начало (корень) сосуда. Ветвью сосуда называется геометрическое место непрерывной последовательности точек, принадлежащих сосуду и расположенных между двумя соседними особыми точками (узел-узел, узел-конец, корень-узел) [15, 43]. Ветвь характеризуется в каждой точке пространственной толщиной и направлением. В целом ветвь характеризуется длиной. Для описания

ЗО секция

31) сегмент

ЗИузея (предок)

50 узел (первый наследник)

ЗИузея (второй наследник)

Рис. 3. Пространственная модель древовидной структуры

Трасса сосуда - геометрическое место серединных точек ветви. Секция - структура, характеризуемая местоположением одной точки трассы, пространственной толщиной в данной точке, а также локальным направлением сосуда. Сегмент - последовательный набор секций, полностью описывающий структуру ветви и являющийся ее дискретным представлением. Таким образом, сосуд в целом можно описать в виде набора особых точек (концы, узлы, корень), соединенных между собой сегментами, которые, в свою очередь, описываются набором секций.

Исходные данные Входные данные для восстановления пространственной структуры получаются в результате процесса коронарографии - последнего этапа обследования больного со стенокардией перед операцией,

проводимого в современном ангиографическом кабинете (рис.4). Исследование выполняется путем чрескожной пункции бедренной артерии, через которую под рентгеновским экраном в аорту, а затем -в обе (правую и левую) коронарные артерии (последовательно) вводится специальный катетер с посту-

пающим через него контрастным веществом. Снимается кинофильм в нескольких проекциях. Левая коронарная артерия имеет более сложную структуру, что обуславливает более сложное ее диагностирование (рис.5). На практике получают 5-6 проекций левой коронарной артерии и 2-3 - правой.

a)

Рис. 4. Спецоборудование в коронарографии: a) AXIOMMultista; б) AXIOM Artis BC; в) поворот на первичный угол, г) поворот на вторичный угол

С помощью коронарографии определяется тип кровообращения; количество пораженных артерий; степень стенозирования сосудов (рис.6).

о )

Рис. 5. Иллюстрация двух типов артериальных стволов сердца: (а) левой коронарной артерии, (б) правой

Коронарография позволяет врачу увидеть двумерные проекции сосудов сердца пациента, в то время как реально они имеют трехмерную пространственную структуру. Целью врача является определение пораженных сосудов, идентификация стенозы и на основе этих данных вынесение решения о хирургическом вмешательстве. Проблема заключается в том, что зачастую количественное определение характера поражения сосудов затруднительно по имеющимся ангиографическим снимкам, так как сте-нозированный сосуд может выглядеть критическим на одной проекции, в то время как на другой проекции тот же сосуд может не обнаруживать стенозы.

Рис. 6. Сосудистая система сердца: со стенозой (а) и без стенозы (б)

Визуальное наблюдение медицинским экспертом изображений сердечно-сосудистой системы носит субъективный характер и позволяет дать только качественные оценки патологических изменений. Применение методов цифрового компьютерного анализа изображений позволяет получить объективные количественные данные об объекте исследования.

Процесс катетеризации в клиниках производится на специальном оборудовании (в нашем случае, это установка C-ARM Equipment, [5] см. рис.4), и данные о процессе сохраняются в унифицированном формате хранения медицинской информации DICOM (Digital Imaging and Communications in Medicine) [2,3]. Как правило, в этом формате сохраняются фильмы проекций, кардиограммы и все данные о геометрии съемки. Будем считать, что в нашем распоряжении имеются все геометрические характеристики плоскостей проекций, а также фильмы проекций, которые, в свою очередь, состоят из отдельных кадров. Кадры представляют собой ангиографические снимки сосудов сердца в определенные моменты времени, сделанные под некоторым ракурсом.

Установка состоит из источника рентгеновского излучения и приемника (или камеры), соединенных между собой дугообразным держателем. Держатель не жестко закреплен на подвижной опоре. Такая конструкция обеспечивает две степени свободы камеры в пространстве. Кроме этого камера может перемещаться относительно держателя в продольном направлении, обеспечивая масштабирование изображения. На рисунке 4в стрелкой показан поворот камеры на первичный угол а (в дальнейшем будем

его ассоциировать с географической долготой), который обеспечивается поворотом держателя камеры вместе с опорой относительно балки крепления. Вторичный угол поворота р (аналог географической широты) показан на рисунке 4г. Такой поворот осуществляется путем скольжения дугообразного держателя по направляющей опоры. При этом камера и источник излучения перемещаются по дуге окружности. Таким образом, ракурс проекции определяется двумя углами поворота: а и р.

Рис. 7. Примеры про Съемка производится расходящимся пучком рентгеновского излучения, угол расхождения пучка примерно равен 12°. Масштаб изображения определяется положением тела пациента в пространстве между источником и приемником, а также расстоянием от источника излучения до камеры. В процессе съемки одной проекции движения камеры не происходит. Однако может происходить изменение масштаба, а также движение операционного стола вместе с пациентом в гори -зонтальной плоскости. Съемки правой и левой коронарных артерий производятся последовательно: сначала формируются все проекции правого сердца (рис.7) (время съемки одной проекции составляет 3-5 секунд), после чего врач изменяет ракурс съемки и процесс повторяется.

коронарных артерий

В результате коронарографии получаются данные в формате Б1СОМ, который представляет собой набор ранов. Ран - фильм (последовательный набор кадров), снятый под определенным ракурсом. Частота фиксации кадров (фрэймов) составляет 25 кадров в секунду. Обычно в одном ране содержится от 80 до 130 фрэймов, то есть съемка одного ракурса продолжается в среднем 3-5 секунд. Соответственно в пределах одного фильма можно наблюдать 3-4 полных цикла сокращения сердца.

В каждом ране также хранится информация об условиях съемки (первичный и вторичный углы, расстояние от источника излучения до камеры). В каждом отдельном документе формата БЮОМ содержится в среднем 4-6 ранов с изображением левой

коронарной артерии и 2-3 - с изображением правой. Каждый фрейм представляет собой изображение в формате BMP(GreyScale) размером 512x512 пикселей. Формат DICOM позволяет хранить большое количество сопутствующей информации.

Задачи автоматического выбора наиболее информативных ранов и нахождения по выбранным ранам фрэймов, соответствующих изображениям сердца в одинаковой фазе, являются, пожалуй, одними из наиболее сложных задач процесса восстановления пространственной структуры сердца и составляют первую подзадачу общей задачи восстановления трёхмерной структуры коронарных сосудов - синхронизация последовательности кадров и выделение необходимой для восстановления структуры фазы сердца, описанные ниже. Таким образом, входными данными задачи восстановления являются фрэймы, на которых сердце находится в одинаковой фазе, и информация об условиях съемки каждого из кадров. Проекции трехмерного древовидного объекта задаются следующим набором данных: изображение проекции в формате BMP; первичный угол поворота камеры а; вторичный угол поворота камеры Р; базовая (начальная) точка на изображении проекции, автоматическое определение которой представлено ниже.

1. Восстановление пространственной структуры коронарных артерий по малому числу ангиографических проекций

При разработке компьютерной технологии восстановления пространственной структуры коронарных артерий был сделан ряд допущений. Предполагалось, например, что из фильмов, представляющих различные проекции, можно выбрать кадры, при регистрации которых сердце пациента занимало одинаковое пространственное положение. По причине того, что мы не обладаем достаточной информацией о взаимном расположении плоскостей проекций в пространстве, невозможно восстановление исходного объекта по проекциям, полученным при центральном проецировании. В связи с этим, будем считать, что изображения, полученные центральным проецированием, можно достаточно полно компенсировать, имитировав тем самым параллельное проецирование.

Компьютерная технология восстановления пространственной структуры коронарных артерий состоит из 5 основных этапов [45]: (1) временная синхронизация проекций (идентификация фазы сокращения сердца), (2) компенсация искажений, вносимых центральным проецированием; (3) пространственная привязка изображений проекций; (4) идентификация точек ветвления сосудов и концевых точек на изображениях проекций; (5) пространственная трассировка и оценка толщины сосудов. Результатом выполнения этапов является структура, со-держашая векторное представление пространственных сосудов, которая в дальнейшем используется

для визуализации. Далее рассмотрим более подробно каждый из указанных этапов. Этапы решения задачи пространственного восстановления приведены в структурных схемах на рисунках 9 и 10.

Временная синхронизация проекций

Задача временной синхронизации проекций возникает из-за необходимости выбора для каждой проекции одного кадра из последовательности. При этом необходимо выбрать для различных проекций такие кадры, на которых сердце было бы запечатлено в одинаковой фазе сокращения.

Известно, что движение сердца носит периодический характер. Известно также, что скорость сокращения неравномерна и убывает в фазе своего максимального и минимального расширения до нуля. Исходя из очевидных соображений, можно сказать, что следует искать кадры, на которых сердце находится в фазе максимального расширения. В этом случае можно будет говорить о соответствии проекций одному и тому же пространственному объекту (мы пренебрегаем некоторыми апериодическими процессами, вызванными, например, аритмией или дыханием пациента).

Фильм для каждой проекции состоит приблизительно из 100 кадров. Требуется проанализировать кадры из фильма и выбрать те, которые соответствуют фазе максимального расширения сердца.

Временная синхронизация проекций основана на вычислении коэффициента корреляции соседних пикселей на изображении каждого кадра. При этом считаем, что изображение описывается моделью стационарного, в широком смысле случайного поля с изотропной корреляционной функцией. В результате получаем массив коэффициентов корреляции. Фактически оценивание коэффициента корреляции производится по нескольким выделенным областям с последующим усреднением (рис.8). Локальные минимумы указывают на кадры, которые соответствуют фазе максимального расширения сердца. Можно сортировать найденные минимумы по значимости. Наиболее значимые минимумы будут более глубокими и более симметричными, чем другие.

Рис.8. Зависимость внутрикадрового коэффициента корреляции от номера кадра (1- по области, 2 -усреднённый по всем областям)

Компенсация эффекта центрального проецирования

По причине того, что коэффициент преломления рентгеновских лучей очень мал для любого материала, невозможно изготовить коллиматор, который мог бы преобразовать расходящийся пучок лучей в пучок параллельных лучей.

Данные в формате

ШСОМ

Рис. 9. Общая структурная Вследствие этого результирующие изображения наблюдаются в центральной проекции, даже в случае, когда исходный объект плоский. По таким проекциям невозможно восстановление исходного объекта без привлечения дополнительной информации о геометрии съемки, которой не располагаем. В то же время, для восстановления исходного объекта по параллельным проекциям не требуется привлечения дополнительной информации. Предполагаем, что эффекты, проявляющиеся на изображении от центрального проецирования, можно отчасти скомпенсировать, сделав изображение таким, как если бы оно было получено параллель-

восстановления 3Б объекта

ным проецированием. Предлагаем алгоритм ком -пенсации эффекта центральной проекции. Следует отметить, что полностью компенсировать эффект центрального проецирования возможно только для плоского объекта, в противном случае в принципе невозможно найти значения некоторых параметров. Преобразование координат выполняется следующим образом.

Для каждой точки выходного изображения (ираг,ураг) «прототип» (х0,у0) ищется в исходном изображении. Координаты точки (х0,у0) определяются по формулам:

С

Пространственная трассировка сосудов

1

Совмещение проекций, поиск пространственной базовой точки и базового

Формирование наборов сфер возможных напоавлений

Поиск локальных минимумов и выбор одного оптимального

—Оптимальный трехмерный вектор /

Движение в выбранном направлений]

L

X

Новая трехмерная точка

7-

Проецирование

L

Двумерные вектора на тЬЬёти^тттдтения\

Оценивание радиуса сосуда по проекциям

^¿шшт

ШЩЩмШт

т

Идентификация разветвлений и концов сосудов

Формирование пространственной секции

Конструирование дерева сосудов

Рис. 10. Общая схема алгоритма трехмерной трассировки

a4 (upar + b1h ) - a2(Vpar + b2h)

a,(v + b2 h) - a3(u + b,h)

1v par 2 ' 3 v par 1 '

где a = cos p, b =-cosasin p, a2 = sinasin p,

a, p - первичный и

b2 = sin a,

a3 = 0, a4 = cos a .

вторичный углы поворота камеры, соответственно, И - расстояние от объекта до плоскости проекции. Это расстояние может быть приближенно оценено

хо

как расстояние от камеры до центра сердца (в пикселях). Пример работы алгоритма компенсации центральной проекции показан на рисунке 11.

Рис. 11. Компенсация эффекта центральной проекции: а) исходное изображение, б) обработанное

Улучшение качества изображений проекций

Низкое качество изображений проекций связано с тем, что на используемом медицинском оборудовании невозможно получить высококонтрастные ан-гиографические изображения вследствие опасности превысить дозу облучения пациента. Низкое качество получаемых изображений проявляется в низкой контрастности, неравномерности распределения яркости и низком значении отношения сигнал/шум.

Выравнивание яркости проводится для компенсации неравномерности освещения при вводе изображения [6] (рис.12). Для избавления от импульсных шумов использовалась ранговая фильтрация, основанная на быстрых алгоритмах вычисления локальных гистограмм распределения и их характеристик [7, 47].

ispyfaar

'SU

Рис. 12. Результаты работы методов улучшения качества изображения сосудов: (а) исходное изображение, (б) выравнивание яркости, (в) ранговая фильтрация, (г) фильтр BSMGF, (д) адаптивная пороговая обработка, (е) ранговая фильтрация бинарного изображения и удаление крупных шумовых фрагментов

Ранговая фильтрация способствовала также более точному обнаружению узловых точек на изображении. Для вспомогательных методов восстановления пространственной структуры объекта, таких как поиск начальной базовой точки на проекциях, определение базового вектора, использовалась адаптивная пороговая обработка [7, 47], фильтрация шумов на основе колоколообразного гауссовского согласованного фильтра (BSGMF) [8], ранговая фильтрация бинарного препарата и процедура удаления крупных шумовых фрагментов (рис.12).

Выделение сосудов на изображении колоколообразным согласованным фильтром (BSGMF) [8]

Для выделения сосудов используется колоколо-образный гауссов согласованный фильтр (BSGMF -bell-shaped Gaussian matched filter), импульсная характеристика которого имеет следующий вид:

f (х, у) = A [1 ± С ехр (-d2(х, у)/2а2)] ,

где d - расстояние от точки с координатами (х, у) до центральной линии сосуда (нахождение центральной линии описано ниже); а - определяет распространение интенсивности в сосуде; А - интенсивность фона; С - константа, используемая для отделения сосуда и его окрестностей. Выражение «±« представляет два рассматриваемых типа изображений сетчатки сосудов. Плюс означает, что сетчатка светлее фона, минус -фон на изображении светлее сетчатки сосудов. Коло-колообразный согласованный фильтр Гаусса позволяет отделять колоколообразные объекты (профиль интенсивности визуально напоминает колокол) от фона изображения (частным случаем таких колоколообраз-ных объектов являются сосуды). На рисунке 13 представлен результат работы алгоритма выделения сосудов колоколообразным согласованным фильтром.

Выравнивание яркости

Выравнивание яркости изображения проекции производится двумя методами. Первый из методов (скользящее выравнивание с контрастированием) основан на локальном линейном преобразовании яркости с использованием скользящего «окна»:

1 = ((/ - СО/6^ )• "шах , где "шах- максимальное

значение яркости по полю изображения, /Ш^- ми-

2

нимальное значение яркости в «окне», сгц1 - дисперсия яркости в «окне», / и / - исходное и преобразованное значение яркости в центре «окна».

Размер «окна» выбирается минимальным при ус -ловии, что оно не содержит только точки фона при любом его положении на изображении.

Рис. 13. Результат работы алгоритма выделения сосудов колоколообразным согласованным фильтром (ББОМЕ): а) исходное изображение, б) результат обработки: а =2, размер окна 17

Алгоритм блоковой эквализации Требуется выровнять яркость по полю изображения, но при этом не изменять локальный контраст. Это свойство важно для рентгеновских изображений

сосудов сердца (рис.14), на которых сосуды находятся на фоне переменной яркости и профиль яркости самих сосудов асимметричен. Можно достичь такой цели, используя блочное выравнивание яркости. Изо-

бражение разбивается на блоки одинакового размера, и обработка производится независимо в каждом блоке с использованием соотношения: / = /•(/1//0), где

/, / - средние яркости в исходном блоке и результирующем блоке, /, / - исходное и результирующее значения яркости в каждой точке блока.

1V * Гг-^^лт ~ ^ ■ >

I уРШШЖ<

й- {К* г, / -. ■■!* IV 1"

г)\

Рис. 14. Результат выравнивания яркости: а) исходное изображение (512 х 512), б) блочное выравнивание яркости в скользящем окне (размер блока 64 х 64), в) блочное выравнивание яркости без перекрытия (размер блока 64 х 64), г) скользящее выравнивание яркости с контрастированием (размер окна 21 х 21)

Ранговая фильтрация [7, 47]

Ранговая фильтрация была использована для подавления шумов на изображении. В отличие от линейной фильтрации [7, 47] ранговая фильтрация позволяет сохранить мелкодетальные особенности объектов на изображении, что важно для сохранения структуры стенозы на изображениях сосудов. Ранговая фильтрация основана на построении вариаци-

онного ряда значений яркости точек изображения, содержащихся в скользящем окне. Идея заключается в предварительном исключении из вариационного ряда крайних значений (шумовых «хвостов») и сохранении в качестве значения яркости центральной точки «окна» медианы оставшихся в «окне» «полезных» точек. Описанная процедура применяется многократно (рис.15).

Рис. 15. Результат ранговой фильтрации изображений: а) исходное изображение сосудов (256 градаций яркости), б) - г) изображения после фильтрации: б) размер «окна» 17x17, отсечение хвостов по рангу 10, в) размер «окна» 7x7, отсечение хвостов по рангу 24, г) размер «окна» 7x7, отсечение хвостов по рангу 35

Адаптивная пороговая обработка изображений

После ранговой фильтрации и применения BSGM-фильтра производится сегментация изображения сосудов методом пороговой обработки на основе относительной энтропии [28]. Этот метод основан на построении матрицы частот соседних уровней яркостей пикселей полутонового изображения. Результатом является бинарное изображение сосудов.

Матрица частот взаимных уровней яркостей (GLC-матрица - gray-level co-occurrence matrix). Для определения GLC-матрицы введём обозначения: D -поле изображения, содержащее MxN пикселей xm п, (m, п) е D , G - число градаций яркости

f (х х ) =

J i, j У m, п ' m + k, п +1 '

1; x = i, x k

> m, п > m + k, п

0;xm, фi or xm

= j

ф j

i, j = 0,G -1.

Индикаторная функция f показывает, имеют ли пиксели, расположенные на заданном расстоянии

(к,I), заданные уровни яркости (г,у). Нормализованная вЬС-матрица определяется как

PtJ (i, j) =

Ck j (Л j)

HQ, (i, j)

Ck,l (i, j) ^ ft,j (Xm,n, Xm+k,п+ )

(m,n)еD (m+k п+)eD

к,I = 0, ± 1, ± 2,..., г,] = 0, G-1.

Размер матрицы Ркравен G х G. Пусть t - значение порога, которое делит вЬС-матрицу на четыре квадранта А, В, С, Б, как показано на рис.16.

Данные четыре квадранта могут быть в дальнейшем сгруппированы в два класса. Если будем считать, что пиксели, имеющие уровень яркости выше порогового значения, принадлежат сосуду, а пиксели, имеющие уровень яркости ниже порогового значения, принадлежат фону, квадранты А и С соответствуют локальному переходу внутри

фона (В), обозначаемому как ВВ, и сосуду (Б), обозначаемому как ББ, соответственно.

0 1: Ь-1 t

Ь-1

Рис. 16. Деление ОЬС-матрицы

Аналогично, квадранты В и Б представляют переходы через границу фон-сосуд и могут быть обозначены как ВБ и БВ, соответственно. Пусть

00 = {0,..,(}, О1 = {( +1,.., Ь -1}. Тогда четыре квадранта А (ВВ), В (ВБ), С (БВ), и Б (ББ) определяются диапазонами яркостей О0 х О0, О0 х 01, 01 х О0, и

01 х 01, соответственно. Тогда вероятности, ассоциированные с каждым квадрантом, могут быть получены как

( ( ( Ь-1

РА =ххР,, Р=хх Р,,

1=0 , =0 1=0 , =(+1

Рс1 = х¿р,, РБ = x x р,, Р, = РаО',,).

¡=(+1 ,=0 1=г+1,=г+1

Выбор порога

Относительная энтропия используется для измерения информационного расстояния между двумя сигналами. Чем меньше относительная энтропия, тем более близкими являются два сигнала в смысле вероятностного распределения. Пусть два сигнала размерности Ь описываются распределениями вероятностей р = (Р0,...,рЬ-1) и к = (к0,...,кЬ-1) соответственно. Относительная энтропия между р и к (или энтропия к относительно р), обозначаемая 3(р;к), определяется как [28]

Ь -1

3 ( р; к) = х р, 1о§ (Р,1к,).

,=0

Чем меньше относительная энтропия, тем меньше расстояние между Р и к, то есть лучше согласование Р с к. Это определение базируется только на гистограмме уровней яркости и не учитывает пространственной зависимости значений яркости пикселей на изображении. Заметим, что относительная энтропия несимметрична. Считаем, что исходное изображение всегда обозначается как номинальное изображение Р, а полученное после пороговой обработки - как к. Последнее должно согласовываться с первым. Это приводит к идее использования вЬС-матрицы для распространения относительной энтропии первого порядка на адаптивную пороговую обработку второго порядка. Поступая так, полагаем, что ( - это выбранный порог. Назначая 1 всем уровням яркости в О = {(+ 1,..,Ь -1} и 0 всем уровням яркости О = {0,.., (}, получаем бинарное

изображение с соответствующей GLC-матрицей, которая задаётся для каждого квадранта как:

P> р>

q' = A q' =_-B_

(' + 1)2 ' (' + 1)(L - ' -1)'

P' P'

q' =_tC_ q' =_-D_

(' +1)( L - ' -1)' (L - ' -1)2

Таким образом, приходим к методу пороговой обработки, в котором порог находится из соотношений:

'jre = arg[min Н

JRE

Hjre (') = -(Р'Л log qA + PB log q'B + P< log q'c + P'D log q'D ) .

Рисунок 17б показывает результаты применения алгоритма пороговой обработки к изображению, обработанному BSGM-фильтром. Видим, что выделенные кровеносные сосуды были подчёркнуты и отделены от фона.

Устранение шумовых фрагментов

После пороговой обработки получаем бинарное изображение, которое содержит фрагменты «мусора». Под «мусором» понимаем отдельно расположенные фрагменты изображения, имеющие малые размеры и распределённые по всему полю изображения. Такие фрагменты являются препятствием при построении качественного ске -летного изображения на следующем шаге. Предлагаемый алгоритм подавления мусорных элемен -тов состоит из следующих шагов: 1) определение максимальных размеров фрагмента, который необходимо стереть, 2) выделение фрагментов, которые имеют размер меньше, чем определённый максимальный размер, 3) стирание отмеченных фрагментов (рис.17в).

Морфологическое восстановление связности

сосудов на 2D изображениях сердца

Задача морфологического восстановления связности сосудов заключается в выборе оптимальных пар соединяемых точек разрыва [31]. Точки разрыва планируется определять с использованием скелета изображения. Для того, чтобы осуществлялась процедура соединения опреде-лённых пар точек разрывов линий, необходимо выполнение двух условий: 1) расстояние между разрывами должно быть как можно меньше; 2) направления линий должны быть согласованы, т.е. не давать при соединении изгибов (рис.18).

Второе условие требует оценки ориентации линий в окрестностях точек разрыва. Для этого необходимо использовать не само изображение, а его скелет. Однако обычные процедуры скелети -зации бинарного изображения приводят к появлению новых разрывов, что делает задачу ещё более сложной.

А В

С D

Рис. 17. Результаты применения алгоритмов пороговой обработки и стирания мусора: а) изображение, обработанное фильтром ББОМЕ, б) бинарное изображение - результат пороговой обработки,

в) изображение после устранения шумовых фрагментов

ми, сшах - максимальный косинус угла между касательными к линиям в точках разрыва. Тогда для всех 1 = 1,...,N решается следующая задача [31]:

Рис. 18. Пары точек разрыва, подлежащие и не подлежащие соединению

Соединение разрывов. При соединении разрывов для каждого разрыва ищется наиболее подходящий разрыв-кандидат (будем называть его сопряжённым разрывом).

Первый подход к поиску сопряжённых разрывов основан на поиске кандидатов с учётом следующих критериев:

■ расстояние от рассматриваемого разрыва до разрыва-кандидата;

■ разница в наклонах касательных в обеих точках разрыва;

■ угол между касательной в рассматриваемой точке разрыва и прямой, соединяющей разрывы.

Фактически используется третий из перечисленных критериев, первые же два играют роль ограничения. Для оценки касательных используются два метода [31]:

■ метод, основанный на усреднении векторов, соединяющих точку разрыва с элементами скелета в окрестности точки разрыва;

■ метод, основанный на аппроксимации кривой в окрестности точки разрыва кривой второго порядка (для соединения разрывов линий в точках с большой кривизной линии скелета):

А20х2 + 4)1 У2 + Апху + Аюх + 4,1 У + А00 = 0 .

Пусть Б = {Ъ1 - множество радиус-векторов

точек разрыва; Т = {(¡. }д=1 - множество векторов касательных к линиям в точках разрыва, - максимальное расстояние между соединяемыми разрыва-

• (bj - b )|

Mb - b

^ min ,

j=1,...,N, j*i

b -

< dm

\t, • t,

\t,

■ < cm

Второй подход для поиска сопряжённого разрыва основан на минимизации некоторой платёжной функции. Критерием является угол между одной из касательных и прямой, соединяющей точки разрыва. Для всех 1 = 1,...,N решается следующая задача: /(Ь,Ь,,(,,(,) ^ шт , где [(Ь,Ь,,(,,(,)- платёж-

,=1,...,N,,=1 1 ] 1 ]

ная функция. Экспериментально подобрана платёжная функция, которая даёт наилучшие результаты:

f (bi, bj, t,, tj) =

\d, а > amln. = pd + qa + r ß + m <

[0, eise;

+ n

d, ß>ßm

0, eise;

где р,д,г,т,п,аШп,РШш - экспериментально подобранные параметры, ё - расстояние между соединяемыми разрывами, а - косинус угла между касательными к линиям в точках разрыва, в - косинус угла

между касательной в рассматриваемой точке разрыва и прямой, соединяющей точки разрыва. Результат применения этого подхода к тестовому изображению приведён на рисунке 19.

Скелетизация бинарного изображения. Основная идея построения алгоритма скелетизации заключается в выполнении многоэтапной логической фильтрации бинарного изображения с использованием скользящего окна 3х3. Логика этого фильтра меняется на каждом этапе обработки. На первом этапе логический фильтр направляется на выполнение операции утоньшения объекта (удаление краёв объекта). При этом все точки разрыва, соответствующие концам сосудов (в том числе, концам ложных отростков), должны сокращаться.

Рис. 19. Работа процедур соединения разрывов: а) исходное изображение, б) соединение разрывов

Сканирование изображения «окном» выполняется по всему изображению многократно. Итерации осуществляются до достижения стационарного (неизменяющегося после очередной итерации) бинарного изображения (рис. 20). Этот этап совпадает в основном с традиционными алгоритмами скелетизации би -нарных изображений [7, 47].

Удаление отростков. В результате выполнения первого этапа получаем 8-связный «скелет» толщиной в 1 пиксель. На следующем этапе производится удаление ложных «отростков». Ложными отростками на изображении дерева сосудов называем шумовые фрагменты изображения, соответствующие ветвям, имеющим длину меньше некоторой заданной величины й.

Рис. 20. Работа процедуры скелетизации сосудов:

а) изображение с соединёнными разрывами, б) изображение после процедуры утоньшения линий

Для удаления ложных отростков производится сканирование изображения по скелету, начиная с концевых точек. Логический фильтр 3х3 настраивается на стирание точек при условии, что их удале -ние от концевой точки не превышает йшх. и не достигнута точка ветвления. На рисунке 21 приведён пример скелета с отростками и скелета, полученного описанной обработкой при заданной мере длины отростков, равной 7.

Пространственная привязка изображений проекций По имеющимся у нас данным о геометрии съемки проекций можем определить пространственную ориентацию плоскостей проекций.

Рис. 21. Работа процедуры удаления отростков: а) исходное бинарное изображение сосудов, б) скелет с отростками,

в) результат удаления ложных отростков (йшвх = 7 )

Но для того, чтобы произвести восстановление исходного объекта, необходимо согласовать положение плоскостей в пространстве посредством па-

раллельного переноса.

Для пространственной привязки проекций необходимо на всех изображениях отметить точку (базовую точку), которая бы являлась проекцией одной и той же трехмерной точки на соответствующие плоскости. В качестве такой точки была выбрана точка конца катетера. Будем считать, что в пространстве концевая точ-

ка катетера соответствует началу координат. Далее параллельным переносом проекций можно добиться соответствия найденных базовых точек соответствующим проекциям начала координат. Заметим, однако, что идентификация конца катетера на используемых для построения трассы кадрах затруднена, так как после ввода контрастного вещества конец катетера перестает быть видимым.

Предлагается следующая трехэтапная процедура идентификации базовой точки на изображении:

Идентификация точки, где катетер пересекает границу кадра на одном из начальных кадров из последовательности. Эта точка оказывается более подходящей для экстраполяции, чем конец катетера. Можно использовать любой из начальных кадров, кроме первого, потому что его яркостные характеристики существенно отличаются от остальных кадров. Экстраполяция найденной точки на кадры, используемые при построении трасс сосудов. Зная позицию точки пересечения катетера с границей кадра, можно, двигаясь по катетеру, достигнуть его конца. В найденной точке также оценивается направление по катетеру (базовый вектор), которое используется при конструировании начального пространственного направления движения.

Результат работы алгоритма показан на рис.22.

Рис. 22. Определение базовой точки: а) бинарное изображение, б) исходное изображение с отмеченной базовой точкой

Мы определяем базовые точки и базовые вектора на каждой из доступных проекций. Считаем, что базовая точка расположена в начале пространственной системы координат. Эта точка вместе с базовым вектором используется далее для совмещения проекций.

Восстановление трёхмерных координат точки

по её координатам на двух плоских проекциях

Нам даны измеренные координаты (и,М, VМ) и (иМ, vM) проекций точки М на первую и вторую плоскости проекции (координаты на изображениях), требуется получить координаты (хМ, уМ, гМ) точки М. Считаются известными описания координатных систем на ортогональных проекциях, которые определяются началом координат и компонентами направляющих векторов. То есть, даны величины:

01 = (Х0^ Уol, ^Д «1 = (хи ^ Уul, ^лХ ^ =(У^ О

02 = (Х02, Уо2 , «2 = (хи 2 , Уи 2 , ^ 2), ^ = (Xv2, У, 2, ^ 2)

Все координаты задаются в одной системе координат (рис.23). Получение трёхмерных координат в конечном итоге сводится к решению системы четырёх линейных уравнений с тремя неизвестными. Правая часть включает в себя координаты проекций, а матрица системы при неизменных параметрах является константой:

a,, a,2 a,3 и,м - a,4

ь,, Ь,2 Ь,3 XM Ум = v,M Ь,4

a2, a22 a23 U2M - a24

ЬИ Ь22 Ь23 L zM \ V2M - Ь24

M(xM,yM,zM) у

а) X

IM

• \

NX Mi(UiM'V1M) u2

'öX

Xй 1 --—бд

б)

Рис. 23. Координатные системы в пространстве (а) и на проекциях (б)

Выпишем выражения коэффициентов а, и Ъ, через параметры координатных систем.

а11 = ах1 хи 1 + ау1 Уи 1 + zu1, а12 = Ъх1 хи 1 + Ъу1Уи1 + Ъz1 ^ а13 = Сх1 хи1 + СУ1 Уи1 + Zu1,

а14 = - Х01) Хи1 + (- Уо1) Уи1 + (4 - 201) Zu1, Ъ,, = а ,х, + а ,У , + а,z ,,

11 х1 ,1 у1у ,1 z1

Ъ12 = Ъх1 Х,1 + ЪуУ, + К 1 Zv1,

Ъ13 = Сх1 + СУ1 У,1 + Zv1,

Ъ14 = - Х01)Xv1 + - У01)Уv1 + (4 - Z01 )Zv1.

Коэффициенты ау и Ъг} вычисляются по тем же формулам с заменой индекса 1 на 2. Величины ах1 и т. д., в свою очередь, вычисляются по следующим формулам:

А2 + с,2 АВ

ах1 = —-—; Ъх1 =——; сЛ =-Е1 Е1

В4. Ъ = А2 + с,2.

AC.d =-AD.

a, = -

* E

;ь, =-

E, bc

E,

BD

C,A,

E

C,B,

c, =--LA; d, =--

yi E

E

A,2 + B2.. d =- QA.

E,

E,

Е, Е1

Е = А,2 + В,2 + с,2.

Аналогичные формулы имеются для ах2 и др. В этих формулах А,, В,, С,,Б, коэффициенты в уравнении первой плоскости проецирования. Они даются выражениями:

А1 = Уи1 - Уv1 zu1, В1 = Xv1 zu1 - Хи 1 zv1,

С1 = Хи 1^1 - Xv1 Уи1, Б1 = - (А1Х01 + В1У01 + С1 )

и аналогично для второй плоскости.

Выше приведённая система четырёх линейных уравнений с тремя неизвестными является переопределенной. Её решение (хм, ум, гм) находится с использованием метода наименьших квадратов. Идентификация точек ветвления и концов сосудов

Задача обнаружения точек ветвления сосудов и концевых точек сосуда возникает на этапе построения пространственной трассы сосудов. Пространственное ветвление появляется при наличии ветвления сосуда, по крайней мере, на одной проекции. Вследствие этого для обхода всего дерева сосудов необходимо обнаруживать ветвления сосудов на проекциях. Также необходимо знать, в какой момент прекращать трассировку пространственного сосуда, и для этой цели также необходимо обнаруживать двумерные концевые точки. Помимо этого, необходимо давать различие между сосудом, разветвлением сосудов и пересечением сосудов.

Метод основан на анализе локальной круглой области вокруг исследуемой точки [12, 29, 33, 39]. Небольшой сектор выбирается внутри круга. Он поворачивается от 0 до 360 градусов с наложением. Таким образом, обрабатываем круглую область скользящим сектором (рис.24).

выбранного сектора Для каждого положения сектора могут быть вычислены некоторые параметры области, которая лежит в пределах сектора круга. Например, среднее значение интенсивности и дисперсия. Эти статистические функции могут быть названы радиальными, потому что они зависят от позиции лучей. Используя эти параметры, можно определить, в каком направлении исходит сосуд из центральной точки сосуда.

Например, если в качестве радиальной функции используется усреднение, то в направлении сосудов, а также темных пятен будут минимумы (рис.25).

Для увеличения надежности необходимо выполнить лучевое преобразование для некоторого диапазона радиусов.

Для каждого значения радиуса находится усредняющая радиальная функция, и определяются ее локальные минимумы. Эти минимумы накладываются

на исходное изображение в соответствии с направлениями, в которых они появляются. Расстояние от центра равно радиусу, который был использован при радиальном преобразовании.

Средняя интенсивность il 2 3 4 5

пгг

j_I_I_s

0 30 70 165 205 285

Рис. 25. Графический пример радиальной функции (для случая усреднения интенсивностей по сектору)

Таким образом, получаем набор точек на изображении, который формирует цепочки точек.

Те цепочки точек, которые соответствуют сосуду, следуют непосредственно до центра. Если можно обнаружить некоторую цепочку точек от центра до максимального радиуса лучевого преобразования, то данное направление соответствует направлению по сосуду. Результаты работы алгоритма для различных видов сосудов показаны на рисунке 26.

Рис. 26. Последовательность локальных минимумов для различных конфигураций сосудов

Пространственная трассировка сосудов

Необходимо разработать метод построения трассы сосудов в трехмерном пространстве, используя изображения проекций исходного объекта, полученные в различное время при условии, что объект может со временем изменять свою форму в некоторых пространственных пределах. Предполагаются известными параметры плоскостей проекций (т.е. их нормали), которые определяют относительное положение проекций в пространстве. Также известны начальные точки на проекциях (базовые точки) и 3Б вектор, который определяет начальное направление движения в пространстве. Известны также углы поворота проекций вокруг своих нормалей в мировой системе координат. Общая схема алгоритма 3Б трассировки представлена на рисунке 10.

Прежде чем начать трассировку (построение трассы), необходимо совместить проекции в пространстве. Ориентируем проекции соответственно их первичным и вторичным углам поворота камеры, а также делаем сдвиг таким образом, чтобы 2Б базовая точка на каждой проекции совмещалась с соответствующей проекцией 3Б базовой точки (рис .27, 49).

Рис. 27. Иллюстрация процедуры согласования проекций

В основе алгоритма лежит пошаговое движение точки в 3Б пространстве с коррекцией направления движения, которая производится на каждом шаге, на основании анализа пространственной интенсивности [13, 44]. По сути, пространственная трассировка является обобщением принципов двумерной трассировки [15, 16, 34, 43]. Расположив в пространстве изображения проекций, с помощью алгоритма обратного проецирования получаем некоторое пространственное распределение интенсивности.

Для выбора оптимального направления движения производится анализ критериальной функции /хугЯ (рраспределенной на сфере радиусом Я

с центром в точке (х,у,г). Сферу с распределенной на ней критериальной функцией в дальнейшем будем называть сферой возможных направлений, подразумевая то, что на ее основе будет производиться выбор оптимального направления движения на каждом шаге. В простейшем случае эта сфера является срезом пространственной интенсивности.

Для анализа пространственной интенсивности и выбора направления движения используется полусфера возможных направлений (рис.28). Пусть также дана точка р0(х0,уо,1(0) пространства, которая является текущей точкой трассы сосуда. Будем считать заданным

начальный вектор у0 (,уу,у2), соответствующий начальному направлению движения в пространстве.

Рис. 28. Виды полусфер возможных направлений

Построение полусферы возможных направлений происходит следующим образом. Введем новую систему координат О'X'У' 1' в пространстве таким образом, что ось О 1' будет сонаправлена с вектором V, то есть плоскость О'Х'У' будет перпендикулярна ему. Точка О' в координатах мировой системы координат будет совпадать с текущей точкой трассы. В координатах новой системы координат уравнение сферы запишется в виде х2 + у2 + г2 = Я2. При этом нас будет интересовать только верхняя полусфера, поэтому накладывается дополнительное условие г>0.

По аналогии с двумерным случаем [22] набор интенсивности на поверхности полусферы может производиться несколькими способами: набор пространственной интенсивности собственно точек сферы; набор с усреднением по лучам; набор пространственной интенсивности с усреднением по конусу; набор с усреднением по цилиндру.

Выбор методов набора обуславливается качеством изображений проекций, а также методом поиска минимумов на сфере, представленным ниже.

На втором этапе определения оптимального направления движения необходимо найти область минимальной интенсивности, центр которой и даст искомое пространственное направление по сосуду. Разработан подход к поиску минимумов: метод, основанный на квадратичной аппроксимации интен-сивностей. Поиск локальных минимумов здесь основан на аппроксимации функции (обозначим ее /(х, у) ) в пределах скользящей маски квадратичной поверхностью:

Е = ах1 + Ьу2 + сху + йх + еу + / .

Далее определяем точку минимума функции Е(х, у) . Если такая точка существует и лежит достаточно близко к центру маски, то она фиксируется. Стандартный метод нахождения минимума функции

двух переменных состоит в нахождении критиче-

- дЕ 0 дЕ 0

ской точки из условия — = 0, — = 0 и проверки

дх ду

выполнения в ней достаточных условий минимума:

д2Е „ д2Е д2Е д2Е „

-> 0,----> 0 .

дх2 дх2 ду2 дхду

В нашем случае вторые частные производные являются константами:

д2 Е „ д2 Е д2 Е

—— = 2а, —— = 2Ь,-= с .

дх ду дхду

Это позволяет дополнительно сократить объём вычислений.

Так как нельзя требовать близости к началу координат меньше, чем на 1 (в силу дискретизации и наличия шумов), то получаем не отдельные точки минимумов, а группы точек вокруг истинной точки минимума. Если количество точек в группе мало, то эта группа игнорируется - этот минимум ложный, полученный за счёт влияния шумов. В результи-

рующий список точек локальных минимумов заносится по одной точке (центру тяжести) из каждой группы. На примере, изображённом на рисунке 29, имеем 6 точек, соответствующих локальным минимумам, которые обозначены цифрами.

Рис. 29. Точки локальных минимумов

Минимумы критериальной функции на 3Б сфере определяют новое направление движения. В случае, когда на сфере детектируется не один, а несколько минимумов, считается, что текущая точка представляет точку 3Б разветвления сосудов. В этом случае все найденные направления запоминаются и трассировка производится последовательно по всем направлениям.

После обнаружения минимумов новое 3Б направление найдено и в данном направлении производится шаг. Новая 3Б точка и 3Б направление проецируются на все исходные проекции и, таким образом, получаем новые 2Б точки и 2Б направления на проекциях, которые являются начальными данными для следующего шага. Процесс трассировки в заданном направлении прекращается, когда достигаем 2Б конца сосуда на любой из проекций. На каждом шаге 3Б трассировки производится оценка толщины сосуда. Все результаты работы 3Б трассировки, включая оценку толщины сосуда на каждом шаге, сохраняются в специальную структуру, которая называется «древовидной» структурой, или просто «дерево» сосудов.

2. Оценивание геометрических параметров древовидных структур

Локальными параметрами, получаемыми в результате трассировки ветви, как в двумерном, так и в трёхмерном случае, являются диаметр и направление ветви в каждой точке. Глобальные параметры (средний диаметр, прямолинейность, чёткообраз-ность, амплитуда колебаний толщины, частота колебаний толщины, извилистость толщины, частота колебаний трассы, амплитуда колебаний трассы, извилистость трассы) вычисляются на основе локальных [15, 16, 17, 19, 46]. Ниже приведём определение данных параметров в двумерных структурах и их отличие для трёхмерного варианта. Введем в рассмотрение также специфический параметр объёмного изображения, такой как «частота вращения», характеризующий такое пространственное свойство ветви, как отклонение трассы от плоскости. Узлы

структур и соединяющие их ветви могут располагаться как в пространстве, так и на плоскости.

Вычисление признаков участка сосуда Ветвь древовидной структуры можно описать двумя типами признаков: локальные (первичные) и глобальные, вычисляемые на основе локальных.

Локальные признаки включают в себя трассу (геометрическое место серединных точек), распределение диаметра и направления ветви вдоль трассы и рассчитываются непосредственно по изображению. Глобальные признаки, представленные ниже, характеризуют всю ветвь в целом.

& Т,.

Рис. 30. Математическая модель фрагмента сосуда

Математическая модель ветви двумерной древовидной структуры определяется функциями (рис.30): X = х(0, У = У (г), г = г (г ),0 < г < Ц , где х(г), у(г)-дифференцируемые функции, описывающие линию центров - трассу [18, 41]; г(г) - функция радиуса ветви (расстояние от трассы до границы ветви, отсчитываемое по перпендикуляру к трассе); г - расстояние от начала трассы, измеренное по трассе; Ц - длина трассы. Данные характеристики однозначно определяют функцию направления трассы в каждой точке р(г); функцию локальной высоты /г), определяемую расстоянием от текущей точки трассы до ее проекции на отрезок Ь, соединяющий начальную и конечные точки трассы; конфигурацию границ ветви - стенок х, = х,(г), У, = У, (г),

х2 = хг(г), у2 = У2(г), 0 < г < Ь .

И в плоском случае, и в объёмном локально сосуд можно считать участком цилиндра. Проекция цилиндра на плоскость, не перпендикулярную его оси, является фигурой с поперечным размером, равным диаметру цилиндра (кроме концевых участков). Поэтому все признаки ветви, относящиеся к толщине - средний диаметр, чёткообразность, амплитуда толщины, частота толщины, извилистость толщины (стенки), а также параметр прямолинейности - имеют одинаковое выражение как для двумерной, так и для трёхмерной древовидной структуры [14, 42]. Введём определение приведённых выше признаков ветви.

Средний диаметр.

Средний диаметр ветви Б вычисляется по фор- _ 2 Ь _ муле: Б = 2г = — I г (г) А, где г - средний радиус

ь 0

v 0

ветви.

Прямолинейность.

Прямолинейность ветви Р характеризует отклонение ветви от прямолинейного направления и определяется как отношение длины Ьв средней линии ветви к длине Ь прямой, соединяющей начальную и конечные точки трассы:

A = 2 f ■ E(k)/

P =

L

L

L <J(m - X(L„))2 + (y(0) - y(L ))2 '

Чёткообразность.

Чёткообразность ветви 8 характеризует неравномерность толщины ветви и определяется как отношение среднеквадратичного отклонения радиуса

ветви к его среднему значению: £ = Vг 2 - Г2 /г , где

г - средний радиус ветви, г2 - средний квадрат радиуса.

Амплитуда колебаний толщины. Амплитуда колебаний толщины ветви А0 характеризует отклонения стенок ветви от прямой линии

и определяется по формуле: А0 Частота колебаний толщины. Частота колебаний толщины ветви ю0 характеризует частоту изменения направления стенок ветви на единицу ее длины и определяется по формуле:

a0 = argmax

J r (i)e™'dt

Извилистость толщины.

Извилистость толщины I0 характеризует скорость изменения функции толщины вдоль трассы, аппроксимируемой гармонической функцией с амплитудой A0 и частотой a0, и вычисляется по формуле: I0 = A0C0 •

Извилистость трассы (для двухмерной структуры дерева).

Извилистость трассы I1 характеризует скорость изменения функции трассы (траектории ветви) на выделенном участке, которая аппроксимируется гармонической функцией Al и частотой a1, и по определению равна I1 = A1a1 и определяется из

2I уравнения: P = — Jl + lf ■ E(k), где k = __1_, где

P - прямолинейность ветви, E (k) - полный эллиптический интеграл 2-ого рода:

л

E(k) = JV 1- k2 sin21dt.

0

Амплитуда колебаний трассы (для двухмерной структуры дерева).

Амплитуда колебаний трассы A1 характеризует степень отклонения хода трассы от прямолинейного и определяется по формуле:

ln (li )

1 + ;

Iiji+Ií

____1 Lv

где f - средняя высота ветви: f = — J |f (t)|di, fi) -

Lv 0

функция высоты ветви, I1 - извилистость трассы.

Частота колебаний трассы (для двухмерной структуры дерева).

Частота колебания участка ветви (трассы) a1 ха -рактеризует число изменений направления трассы на единицу длины ветви и определяется по формуле: с1 = 11/А1.

Рассмотрим определение признаков трассы ветви (извилистость, амплитуды колебаний и частота колебаний трассы) для трёхмерной структуры дерева. А также введём ещё один признак, характеризую -щий отклонение трассы от плоскости. Будем называть его частотой вращения. В качестве пространственной модели ветви рассмотрим участок спирали:

[y = R sinarx

[z = R cosarx'

где ar - частота вращения; R - расстояние до оси Ох, совмещенной с прямой, соединяющей начальные и конечные точки трассы.

В свою очередь, R является гармонической функцией: R = A sinax, где А - амплитуда трассы; с - частота трассы.

В результате получим параметрическое уравнение трассы:

x = t

< y = A sina isinart.

z = A sin a t cosat

В системе имеются два параметра: А - амплитуда трассы и a - частота трассы, которые аналогичны соответствующим параметрам в двумерном случае, и один параметр ar , специфичный для трехмерной модели - частота вращения трассы.

Введём следующую функцию расстояния точки трассы от оси Ox:

R(x) = Vy2 (x) + z2(x) .

Учитывая параметрическое уравнение трассы, получаем:

R(x) = AIsina^ .

Согласно этому уравнению среднее значение R = 2A / n, откуда имеем выражение для амплитуды трассы: A = nR /2.

Значение R(x) для каждой точки трассы вычисляется с помощью известной формулы для расстояния от точки до прямой. Выражение ам-плитуды трассы даёт приближённое значение па-

2

раметра. Выше приведена точная формула оцени -вания параметра в плоском случае, которая при малом значении Arn принимает представленный выше вид. В силу того, что после перехода от двух проекций к трёхмерным координатам точки трассы не лежат равномерно по какой-либо вели -чине, невозможно вывести точную формулу. Ос -новываясь на результатах в плоском случае, можно показать, что поправка к амплитуде трассы невелика - в плоском случае при Arn ^ 0 получаем представленную выше формулу, а при Arn ^ да получаем A = 2R [16, 40].

Для вычисления частоты трассы используем

преобразование Фурье. Функция R(х) = A|sin®x|

является периодической, но не является синусоидой, поэтому для повышения точности используем следующий подход: выполняем преобразование Фурье

от центрированной функции R 2(х), где R 2( х) = -j A2(1 - cos2®x), и получаем оценку частоты трассы:

ю = -2 arg max 3(R2(х)),

где 3 - дискретное преобразование Фурье.

Частоту вращения трассы теоретически можно найти аналогично частоте трассы:

юг = arg max 3(sinp(х)),

где р( х) - полярный угол проекции точки трассы на плоскость yOz . Однако при таком подходе неизбежна большая погрешность оценивания в силу малого размера анализируемого участка сосуда и малого его отклонения от плоскости (точка проекции совершит заметно меньше полного оборота вокруг начала координат). В то же время спектральный анализ даёт приемлемые результаты, если имеется не менее одного периода разлагаемой функции. Используем следующий подход: для кривой (параметрическое уравнение трассы см. выше) значение прямолинейности Р :

p = -VITA2®7 e( k),

п

k2 =

A2(®2 -ю2) 1 + A2®2

где Е() - полный эллиптический интеграл 2 рода. (Очевидно, что при аг = 0 получаем формулу для плоского случая). Используя найденные значения величин Р , А и а , можно получить оценку частоты вращения трассы:

ю = Лю -

k 2(1 + A ю )

A2

(

где Е-1 (•) - функция, обратная к эллиптическому интегралу 2 рода.

3. Методы оценивания толщины сосудов

Разработано несколько методов оценивания параметров толщины сосудов по их изображениям. К ним относятся методы непосредственного и косвенного измерения параметров [16, 17, 36, 40]; методы, основанные на моделировании; методы, основанные на применении искусственного интеллекта и нейронных сетей [30]; комбинированные методы [20]. В основе аппроксимационных методов [21] оценивания лежит идея параметрической аппроксимации яркостного профиля, выделенного на изображении фрагмента исследуемого сосуда. Для определения параметра толщины Б объекта предлагается три различных модели яркостного профиля.

Модель профиля яркости изображения ветви древовидного объекта. Экспериментальные исследования показали, что функция яркости сосудов обладает оп -ределенными свойствами, которые необходимо учесть в разрабатываемой модели (рис.31). На рисунке 31 че -рез / (п) обозначена функция яркости поперечного сечения сосуда. Приведем примеры возможных видов функции профиля.

k = E-1

пР

2лД + A2ю2

Рис. 31. Яркостные профили изображений древовидных объектов

Проекция эллипса. Профиль сегмента древовидного объекта в общем случае представляет собой эллипс с известными параметрами и проецируется на изображение в виде полуокружности (рис.32а). Такая модель адекватна в случае ангиографической съемки сегментов больших диаметров.

Прямоугольный импульс. Наиболее простая модель. В этом случае профиль сегмента древовидного объекта аппроксимируется функцией прямоугольного импуль-

са (рис.32б). Такая модель адекватна в случае съемки непрозрачных для данного излучения объектов.

г;

Рис. 32. Различные виды функции профиля Вторая производная от функции Гаусса. Эта модель с подчеркнутыми границами и ярко выраженным центром профиля (рис.32в) хорошо подходит для аппроксимации тонких сегментов древовидного объекта в случае рентгеновской ангиографиче-ской съемки объекта.

Разность второй производной от функции Гаусса и функции Гаусса. В этой модели резко подчеркивается контур сегмента и менее - центральная точка (рис.32г).

Рассмотрим более подробно модели представления профиля изображения древовидного объекта, которые в дальнейшем используются в методах оценивания диагностических признаков ветвей древовидных объектов.

Модель 1. Функция Яей:

х-хс А,|х- хс\ < Б/2; Б ) " 0,1 х - х1 > Б/2.

y( x) = A Rect

Модель 2. Функция Гаусса:

(x - x )2

y(x) = A expl -

D

Модель 3. Дробно-рациональная функция (ДРФ):

A -D

У( X) (x - xc У + D Параметрами моделей являются значение толщины D , xc - координата центра сосуда; A - яркость точки в центре сосуда.

Непосредственная аппроксимация по методу наименьших квадратов

В качестве целевой используем функцию:

T

J (D) = {((( x, D) - y(x))2 dx ^ min,

о

где y(x, D) - модельная функция профиля яркости; y(x) - исходный профиль яркости, T - длина отрезка анализа функции яркости. Аппроксимации подвергается срез функции яркости изображения,

взятый перпендикулярно направлению сосуда. Недостаток этого подхода состоит в значительной вычислительной сложности и необходимости предварительного определения точного направления сосуда. Достоинство - в высокой помехоустойчивости и хорошей точности при правильном подборе модели аппроксимации.

Прямоугольная аппроксимация Представим срез функции яркости в виде прямоугольного сигнала:

17 , |х-х I < Б/2;

у(^ хс , /ш1n, "ш^ Б) = 1™

I /ш■n, иначе,

где хс - координата центра сосуда; /шах - яркость

сосуда; /шп - яркость фона, Б - толщина сосуда.

Для оценивания параметров модели хс, /шП, /шах,

Б воспользуемся методом моментов, приравняв теоретические обобщенные моменты модельной функции к эмпирическим:

Т т

| х"у (х)йх = | х"у(х хс , /ш,„ , /ш,х, Б)йх ,

0 0

к = 0,1,2,3 .

Тогда толщина сосуда оценивается следующим об-

1 т

разом: Б = Т • у , где у = — |у(х)йх, у(х)- норми-

Т 0

рованная к диапазону [0;1] функция яркости сечения сосуда.

Результат работы метода прямой аппроксимации исходного профиля моделью Яей представлен на рисунке 33.

Рис. 33. Пример аппроксимации профиля функцией Яес1

Аппроксимация дробно-рациональной функцией

Модель, основанная на использовании дробно-рациональной функции:

y( x, D) = -

D2

х + Б (1)

где Б - значение, пропорциональное искомой толщине. Для оценивания параметров модели 2 также применяем метод моментов, из которого получаем уравнение относительно оцениваемого параметра Б:

2Б Т _

-аг^я-= у .

Т 2Б

Пример результата аппроксимации, полученной таким образом, можно видеть на рисунке 34.

Рис.34. Пример аппроксимации профиля моделью ДРФ Аппроксимация профиля ступенчато-многочленной моделью

Особенностью ступенчато-многочленной модели является то, что многочлен аппроксимирует область самого сосуда, а константы - фоновую область.

Парабола 2-й степени + ступень (симметричный случай без разрывов, рис.35)

В данном случае аппроксимирующая функция имеет вид (начало координат - в центре симметричного профиля):

с,

х > А;

У=

ах2 + с, х < А.

Условие непрерывности:

ах + с = С.

1х=± А

(2)

(3)

В силу симметрии функция чётная, а вершина параболы лежит в точке (0, с) .

Рис.35. График параболы 2-й степени + ступень,

симметричный случай без разрывов Подлежащий минимизации функционал:

- А

3 = I (С - у(х))2 йх+

-Т/2

А Т/2

+ | (ах2 + с - у(х))2йх + | (С - у(х))2йх.

(4)

Используя условие непрерывности (3) в (4), получаем выражение для эквивалентного функционала:

= |[ а( х2 - А2) + С - у (х) ] 2йх,

где С =

1

Т - 2 А

- А 1 /2

I у(х)йх + I у(х) йх

(5)

Численно минимизируя значение функционалов (4) и (5) с учётом условия непрерывности (3), получаем оценку параметра А и находим оценку толщины сосуда Б:

Б = 2 А,

(6)

2

Коэффициент 3 в (6) выбран эмпирически на

основе анализа реальных изображений сосудов.

В первом приближении модель, описанная выше, вполне достаточна, но обладает следующими недостатками. Если центр профиля недостаточно близок к центру сосуда, то на одном из симметризованных профилей в центре будет не один минимум, а два минимума и максимум между ними. Это не очень согласуется с моделью параболы, в то же время отказ от симметризации сильно увеличивает объём вычислений. Оба вида профиля можно учесть, если использовать вместо параболы многочлен четвёртой степени.

Парабола 4-й степени + ступень (симметричный случай без разрывов)

Аппроксимирующая функция (начало координат в центре симметричного профиля) и условие непрерывности:

ГР, |х| > А;

У Ч 4 _ . . (7)

ах4 + Ъх2 + с

= Р.

(8)

В силу симметрии функция чётная, вершина кривой лежит в точке (0, с) . Точка минимума лежит в точке (х0, у0), координаты которой зависят от числа экстремумов кривой. А именно: если -Ъ/2а < 0, то х0 = 0, у0 = с (рис.36); если наоборот

-Ъ/2а > 0, то х0 = Ъ/2а , у0 = с - Ъ2¡4а

(рис.37). Здесь предположили, что а > 0 - ветви направлены вверх.

Рис.36. График параболы 4-й степени + ступень, симметричный случай без разрывов (один экстремум)

Рис.37. График параболы 4-й степени + ступень, симметричный случай без разрывов (три экстремума)

Подлежащий минимизации функционал выражается формулой:

-А

3 = | (Р - у(х))2йх +

| (ах4 + Ъх2 + с - у(х))2 йх + | (Р - у(х))2 йх.

(9)

С учётом того, что параметр Р оценивается из соотношения

А/

Р=

Т - 2 А

(10)

1 ' I У(х)х + I у(х)йх

-Т /2 А

получаем функционал в виде:

А

= ![ а( х4 - А2) + Ъ( х2 - А2) + Р - у(х) ]йх . (11)

-А

Из условия минимизации функционала находятся значения параметров а и Ъ , параметр с определяется из условия непрерывности (8). Затем находятся численно параметр А как минимум функции (11) одной переменной. Оценка толщины сосуда определяется как Б = 2г, где г - наибольший по модулю корень уравнения:

аг4 + Ъг2 + с = У0 + 2(Р - У0). (12)

Выбор уровня функции вычисления толщины

При определении толщины сосуда существуют два плохо формализуемых вопроса. Первый - выбор размера маски, который в ряде методов равен длине профиля. Будем обозначать размер маски (длину профиля) как 2Я. Будем называть Я - радиусом рамки. Второй - в аппроксимационных методах, если аппроксимирующая функция отлична от константы (прямоугольника), необходимо выбрать, какой уровень значения функции соответствует границе сосуда. Сейчас оба эти вопроса решаются экспериментальным или интуитивным способом. Однако здесь приводятся подходы, позволяющие эти вопросы решить конструктивным образом. Рассмотрим симметризованный и обращённый профиль -сосуд светлее фона. Аппроксимирующая функция считается нормированной (максимум равен единице, минимум - нулю). Аппроксимирующая функция при возрастании аргумента стремится к нулю. Направление сосуда считается известным, поэтому профиль является перпендикулярным к трассе.

Далее будет рассмотрено несколько аппроксимирующих функций и для них будут приведены уровни значений функции и соответствующие толщины.

Прямоугольник (Яесг, рис.38)

В силу простоты модели нет вопроса выбора уровня, но на ней можно привести доказательство правила выбора радиуса рамки, и в дальнейшем, действуя по аналогии, доказать выбор уровня.

Аппроксимирующая функция (Я -радиус рамки):

У( х) = ■

|1, |х| < г; |0, г < |х| < Я.

(13)

Для двухградационного изображения наибольшая дисперсия яркости достигается в случае равенства площадей, занимаемых сосудом и фоном.

0

X

Рис.38. График аппроксимирующей функции Яесг Докажем это для одномерного случая и получим правило для выбора радиуса рамки.

-2 Я г — 2 Я г у=~2я Iу(х)йх=Я' у = 2ЯIу2(х)йх=Я'

Бу = У2 -(у )2 = г/Я - г 7 Я2.

(14)

Максимум функции (14) равен одной четверти и достигается при Я = Я0 = 2г . Отсюда вытекает правило выбора радиуса рамки: оптимальный ра -диус равен двум толщинам сосуда. В таком виде правило нельзя применить на практике, так как толщина ещё неизвестна. Реально это правило можно использовать тремя способами: 1) взять пробное значение радиуса рамки, вычислить начальное значение толщины и на его основе получить новое значение радиуса (возможно выполнение нескольких итераций); 2) аналогично предыдущему, но новое значение радиуса используется не для текущей, а для следующей точки трассы; 3) изменяя размер рамки, найти такой, при котором достигается максимум дисперсии яркости на профиле (не аппроксимированной). Оптимальность понимается в следующем смысле: конфигурация с максимумом дисперсии является наиболее устой -чивой к шумам. При модели прямоугольника толщина г равна параметру г, то есть можно записать равенство:

г = Я0/2. (15)

Мы распространяем формулу (15) на произвольную модель. Толщина сосуда равна половине оптимального радиуса рамки при выбранной модели аппроксимации. Уровень функции, строго говоря, приводится только из соображений наглядности. Далее приведём результаты действия этого правила на нескольких моделях аппроксимации.

Степенная функция (рис.39)

Аппроксимирующая функция имеет вид:

П -(|х|/г)п, |х| < г; У(х) = ] УП' > ' 11 ' |0, г < х < Я.

(16)

Параметр п не обязательно целый, но должен быть положительным.

Вычислим дисперсию яркости:

Ву=г—2п^--г! (17)

Я (п +1)(2п +1) Я2 (п +1)2

а2

0 г Д 2

Рис. 39. График аппроксимирующей степенной функции Максимум функции (17) равен п2/(2п +1)2 и достигается при R = R0 = г ■ (2п + 1)/(п +1). Исходя из формулы (15), получаем значение толщины:

t = г-(2и + 1)/(2п + 2). (18)

Отметим, что, как и следовало ожидать, при п ^ ж получим значения для модели прямоугольника. Значение толщины (18) соответствует уровню

функции: у(х = 0 = 1 - ((2п +1)/(2и + 2))".

Приведём значения уровня для частных случаев.

п Уровень

1 14 = 0,25

2 1136 и 0,306

3 169/512 и 0,33

4 3439/10000 и 0,344

ж 1-]Д/ё и0,393

Из этой таблицы следует, что при модели параболы 2-й степени (формула (6)) - точное значение коэффициента равно не 2/3 и 0,667 , а (с учётом обращения профиля) 25/36«0,694. Что касается формулы (12), то здесь ситуация хуже - в модели параболы 4-й степени присутствуют два слагаемых, кроме того, может быть ямка. Если нет ямки, то коэффициент лежит между 0,656 ( " = 4 ) и 0,694, поэтому 2/3 вполне допустим. Если ямка есть, то обосновать коэффициент, исходя из отдельных степеней, трудно.

Дробно-рациональная функция В данном случае аппроксимирующая функция имеет вид:

У( х) =

г2/(х2 + г2), |х| < R,

(19)

0, |х| > Я

Выражение для дисперсии принимает вид:

11 1 2

°У = ™-Г + Т"аг^(и) —-аг^ и , (20)

2(1 + и ) 2и и

где п = Яг .

Максимум (20) достигается при п и 3,5 (Я и 3,5г). В практических расчётах используется оценка толщины сосуда в виде t = 1,75г, которое соответствует значению функции яркости У^) и 0,246.

Функция Гаусса

Аппроксимирующая функция имеет вид:

У( х) = ■

(21)

(22)

*=и • 2П^ /-7п (и /

2 2 2

где и = Я/а, Ф(2) = е-гА.

0

Чтобы найти максимум этой функции, необходимо решить следующее уравнение:

•Л 1

2л/2

й (Ву)

йи п 1

-■ф

(ил/2 )

+—■ — -ф2 (и /-л/л — -Ф(и )е~"2 = 0.

2 иъ и1

Приближённое (с точностью 0,05) значение его корня равно 2,5, то есть оценка толщины сосуда получается в виде t и 2,5а при отсечении по уровню яркости у^) = е.

Парабола 4-й степени (рис.40) Аппроксимирующая функция имеет вид:

У( х) =

[1-(х/а/ ±(х/Ь )2, X < г; [0, г < IX < Я.

(23)

Знак минус соответствует кривой без ямки, а плюс - кривой с ямкой.

а)

о I Л

б)

1 /

0 г Я х

Рис.40. Аппроксимация параболой 4-й степени: а) без локального минимума, б) с локальным минимумом

Величина г - такое значение х, при котором

у(х) = 0 . Рассмотрим сначала вариант со знаком

минус. Величина г равна:

г = а^Уа4 + 4Ь4 - а Вычислим дисперсию.

Ву = А я - в/Я2 ,

где введены обозначения: 2 г

Ал/2.

(24)

(25)

,12 ^ г5

А = г--+ |---I —+ -

" Ь а4 I 5 7 а 4Ь2

2 г7

3Ь2 ^ Ь4 а41 5 7 а4Ь2 9а8 В =(г - г 3/3Ь2 - г75а4 )2. Максимум функции (25) равен А2/4В и достигается при Я = Я0 = 2В/А . Исходя из формулы (15), получаем значение толщины: t = В А.

В случае варианта со знаком плюс изменяются выражения для г, А, В :

= a^\fa

A = r + -

I a + 4b4 + а 2/b/z, 2r M_ _ 2r1 r9

ibЧ ~ьа ~ a J т _ 7 a 4Ь 2 + 9a7'

B = (r + r73b2 _ r75a4 )2.

Аналогичным образом можно для любой модели аппроксимации найти теоретически обоснованное значение толщины и оптимальный размер рамки.

Спектральные методы

В основе этой группы методов [21, 32] лежит преобразование Фурье. Дискретный спектр F(m) конечной последовательности fn) (дискретный профиль яркости) определяется следующим образом:

F (m) = £ f (n)exp(-/—— nm).

n=0 N

Для обеспечения инвариантности к сдвигу профиля относительно центра сосуда будем производить аппроксимацию квадрата модуля спектра, который является вещественной четной функцией. Таким образом, задача параметрической аппроксимации спектра профиля заключается в минимизации целевой функции:

S(D) = £ (F2 (m, D) _ |F(m)|2 V ^ min,

m=0

N-1

2n

где F(m) = £ f (n) exp(_/ —— nm), дискретный спектр

П=0 N

исходного профиля. Для повышения вычислительной эффективности получения оценки применяется метод последовательного приближения, при этом сначала осуществляется определение толщины каким-либо относительно быстрым методом, а затем методом бинарного поиска производится уточнение полученного значения.

Модель Sync. Эта модель представляет собой спектральный аналог модели Rect. Дискретный спектр Фурье для этого случая:

*(1 + D)

F (т, D) = —

sin-

nm

Высокочастотные компоненты спектра отвечают за шумовые составляющие сигнала, поэтому будем проводить аппроксимацию только главного лепестка спектра, ширина которого обратно пропорциональна ширине сосуда.

Этот факт позволяет нам упростить процедуру оценивания толщины до поиска первого локального минимума спектра исходного профиля. Для ускорения работы метода вычисления проводятся в два этапа. На первом этапе оценка принимается равной оценке, полученной применением модели Ке& на исходном профиле, затем уточняем результат, проводя поиск локального минимума спектра исходного профиля. Для найденного таким образом значения частоты О вычисляем оценку толщины по формуле: Б = . Такой подход позволяет сократить время

расчётов в среднем на 30%. На рисунке 41 приведены результаты аппроксимации одного и того же профиля или его спектра рассмотренными выше моделями.

Рис.41. Применение различных аппроксимационных методов: а) модель Rect; б) модель ДРФ; в) модель парабола; г) модель Sync

Данные методы оценивания толщины были исследованы с использованием имитационного моделирования тестовых изображений, для которых параметры были известны априорно. При этом вид тестовых изображений максимально соответствовал реальным изображениям. В качестве тестового объекта использовались прямолинейные и криволинейные участки ветвей древовидных объектов. Примеры результатов оценивания толщины приведены на рисунке 42.

i ГО "Я х £ ш if a a a> I § ¡4 ||| - s ß ill HI If it Ill % SSS II If Ii «Г 5 ii* III If Hal

ш 6 6.59 e.60 6,84 6,14 5,19

ш 14 13.99 14,67 13.16 13,60 13.05

V 5 5.79 5,91 6,13 6,4В 5.53

J а 5 5.19 5.35 4.32 4.96 4,03

20 19.91 20.02 18,05 18,61 18,26

а 10 9.69 9.92 8.67 9.23 8,52

Рис.42. Экспериментальные результаты по оцениванию толщины сосудов

Достоинством описанных методов является их высокая точность и помехоустойчивость. Недостат-

ком является необходимость указания профиля сечения, а также вычислительная сложность.

Косвенный метод оценивания толщины сосуда

Предполагая, что точка центра сосуда найдена и яркость точек сосуда меньше яркости точек фона, имеем косвенный метод измерения толщины сосуда в данной точке [36].

Метод использует круглую рамку диаметром больше, чем диаметр сосуда в 1,5-2 раза (рис.43), и включает следующие этапы:

Рис. 43. Фрагмент сосуда с круглой маской

1. Формирование одномерного массива чисел, пропорциональных яркости точек изображения, которые лежат на окружности круглой рамки.

2. Применение к сформированному массиву сглаживающего фильтра для уменьшения уровня шума.

3. Пороговая обработка сглаженного массива с целью получения бинарного массива, состоящего из нулей и единиц, для выделения отсчётов изображения, принадлежащих сосуду, а не фону, окружающему сосуд. При этом порог вычисляется на основании значений яркости всех точек, лежащих в квадрате, описанном около круглой рамки.

4. Применение к сформированному на шаге 3 массиву цифрового рангового фильтра для устранения точечных локальных минимумов, связанных с импульсным шумом изображения.

5. Формирование упорядоченного по убыванию угловых длин отрезков нулевых значений сглаженного бинарного одномерного массива и выбор не более четырёх наибольших по длине от -резко в нулевых отсчётов (условно «угловых дуг») по правилу: если число «угловых дуг» меньше че -тырёх, то все «дуги» сохраняются, если больше четырёх, то вычисляется суммарная угловая длина «дуг» до тех пор, пока она не превысит 80% от угловой длины всех «дуг» или пока их число не будет равно четырём.

6. Если число «дуг», найденных на шаге 5, равно двум, что соответствует участку сосуда без пересечения, ветвления и окончания внутри рамки (в противном случае понятие толщины становится некорректным), то диаметр сосуда определяется по формуле:

D = 2R sin—sin 2 2

где R- радиус круглой рамки; в - среднеарифметическое угловых размеров двух «дуг», образованных при пересечении сосуда круглой рамкой:

в = (#1 +£)/2; где у- центральный угол, образованный двумя радиусами круглой рамки, направленными на центры двух «дуг» (рис.44).

Если число «дуг», найденных на шаге 5, не равно двум, то диаметр сосуда не вычисляется, а условно принимается равным нулю и данная точка сосуда

4. Выделение центральных линий сосудов с использованием вейвлет анализа и методов дифференциальной геометрии

Выделение центральных линий В работах [24-27, 48] представлены различные подходы к решению задачи выделения центральных линий: подход, основанный на пороговой обработке и анализе связанных областей, подход, основанный на согласованной фильтрации [25], подход с использованием дифференциальных свойств поверхности изображения [24]. Предлагаемый в работе метод совмещает элементы различных подходов, что позволяет ему избежать ряда недостатков, присущих перечисленным работам. Выделение центральных линий производится на основе анализа нечеткого поля направлений [27], определяемого формулой:

ф (а, т) = Ж (а, т)• е-2"в(°,ш).

Здесь для каждой фиксированной точки ш е 22 величину Ж (а, ш) - весовую функцию поля направлений - можно трактовать как характеристическую

Л -2"Цат)

функцию принадлежности элементов ха = е , а > 0. Оценивание нечеткого поля направлений производится с использованием методов вейвлет анализа, при этом «нечеткость» возникает как результат введения в преобразование параметра масштаба. Предлагается два метода оценивания нечеткого поля направлений: (1) метод, основанный на вейвлет преобразовании с неизотропным порождающим вейвлетом и (2) метод, основанный на анализе главных кривизн изображения.

Непрерывное двумерное вейвлет преобразование с неизотропным вейвлетом Рассмотрим изображение как непрерывную двумерную функцию / (х)е Ь2 (Я2, С2 х), х = (х1, х2) .

Семейство базисных функций вейвлет преобразования получается из порождающего вейвлета посредством сдвигов, растяжений и поворотов последнего

¥а.вЬ (Х) = а ^ Г-в(Х -b)), где Г - оператор

поворота. Тогда прямое вейвлет преобразование имеет вид:

[X/](аД Ь)= С-^| /(х)^~(х)с2х .

Выбор порождающего вейвлета Целью применения вейвлет преобразования является последовательное выделение частот, соответствующих различным деталям объекта с учетом их ориентации. В работе [23] рассмотрены направленные

вейвлеты, которые хорошо подходят для решения подобных задач. В работе рассмотрены два вейвлета. Направленный вейвлет Морле Вейвлет Морле в спектральной области имеет вид:

) = Vsexp[ --2•

£• ш2 +

(ш , - ш Х )2 ]

где го0 - параметр смещения вейвлета в спектральной области, е - параметр, отвечающий за угловой раствор конуса, содержащего эффективный носитель вейвлета в частотной области (рис.45а).

а) б)

Рис.45. Примеры вейвлетов в спектральной и пространственной областях: вейвлет Морле (а) и вейвлет Коши (б)

Направленный конический вейвлет, или вейвлет Коши Конический вейвлет имеет вид:

ц (го ) =

(со,еа)((о,ед )exp1 [ю2 ]J,ю е С(Да)

0, ю г С (в, а)

еа =(sin (а), - cos (а)), ев =(-sin (в), cos (в)).

Эффективный носитель в частотной области лежит в пределах конуса С (а, в), в < а (рис. 45б). В

работе предпочтение было отдано вейвлету Коши. Из практических результатов было выбрано а = -в = П 4. Для решения задачи выделения сосудов (в силу симметрии последних) был сконструирован веерный вейвлет - объединение двух вейвле-тов, носители которых в спектральной области расположены в противоположных конусах.

Реализация непрерывного вейвлет преобразования Дискретизация вейвлета производится в спектральной области. Как правило, толщина сосудов на изображениях такова, что эффективный носитель спектра вейвлета, соответствующего наиболее тонким сосудам, лежит в пределах области, имеющей линейные размеры в 4 раза меньше размеров исходного изображения. Таким образом, вейвлет преобразование для дискретного изображения вычисляется по следующим формулам.

Пусть f(x) - исходное изображение имеет размер N х N, x е Z2.

f (ю) = 3 (f (x)) - ДПФ исходного сигнала считается 1 раз.

Фавн (го) - рассчитывается аналитически для каждого значения масштаба а и угла в, имеет размеры N/4 х N4 . Для всех значений а и в рассчитывается

[#/ ] (а, в, го) = f (го )• фв (го),

[WJ](а,в,b) = З-1 ((/](а, в,го)) .

Результирующий вейвлет-образ при фиксированном значении а и в имеет размер N/4 х N/4.

Вейвлет метод оценивания нечеткого поля направлений

Оценивание весовой функции и направлений производится по следующим формулам - трехмерный набор вейвлет коэффициентов:

W(а,m) = max[Wf](а,в,m), трехмерный набор направлений (углов) (рис.46):

D( а, m) = arg max [Wf ](а ,в, m).

0<в<п

Такое представление позволяет сократить количество избыточной информации и является достаточным для анализа центральных линий.

Метод оценивания нечеткого поля направлений,

основанный на анализе кривизны поверхности

В ряде работ [23, 24, 26, 34] для выделения центральных линий сосудов используются дифференциальные свойства поверхности изображения. В работе [23], например, предлагается метод выделения тонких сетчатых структур на изображении, основанный на анализе максимальной главной кривизны изображения. Мы использовали схожий подход к оцениванию нечеткого поля направлений.

Рис.46. Примеры работы вейвлет метода оценки поля направлений для двух значений масштабной переменной а : а), б) - весовые функции поля направлений, в), г) - угловые функции поля направлений

Будем рассматривать функцию изображения как некоторую регулярную поверхность в трехмерном пространстве. Пусть задана область Ос Я2 на плоскости. Рассмотрим некоторую гладкую функцию (изображение) z=f(x). Отображение г : О^ Я3,

определенное как г :(х) ^ (х, f (х)), задает регулярную поверхность. В каждой точке (х, у) определяется минимальная главная кривизна Л(х, у) поверхности и соответствующее ей главное направление £(ху).

С целью одновременного анализа сосудов, имеющих различную толщину, и с целью уменьшения влияния шумов необходимо производить фильтрацию исходного изображения л х) некоторым низкочастотным фильтром, прежде чем анализировать изображение. Фильтр должен принимать различные значения ширины полосы пропускания, чтобы обеспечить анализ сосудов на различных масштабах. В качестве фильтра рассмотрим гауссову двумерную функцию

- -Ч х| 2

g (а, х ) = е 2а- , где а > 0 - масштабная перемен -ная, х = (х1, х2). Результирующее изображение будет: f (а,х) = Л0 х g (а, х). Для вычисления производных первого и второго порядка (в матрицах первой и второй фундаментальных форм) вос -пользуемся формулой:

(а, х ) = / (а, х ).

дх. У ' дх. V '

Рассмотрим двумерные гауссовы вейвлеты:

уУ (а,х) = -а-3 • х. • g(а,х) ,

у/2 (а,х) = а-3 •((а-1 х1) -1) • g(а,х) ,

уУ (а,х) = а~5 • х1 • х2 • g(а,х) , I = 1,2 .

Тогда частные производные функции f (а, х) в

матрицах первой и второй фундаментальных форм можно рассчитать по формулам:

Л (а, х ) = Ъху (а, х), (a, х) = Лх у. (a, х),

Лх2 (a, х) = /0 х у (а, х) , . = 1 2 .

Таким образом, производится оценивание минимальной главной кривизны и соответствующего главного направления для различных масштабов (рис.47). В этом случае оценкой весовой функции нечеткого поля направления является функция минимальной главной кривизны изображения, а оценкой направлений - соответствующие главные направления: Ж (а, т) =Л(а, т), Б( а,т) = £ (а, т).

Метод выделения центральных линий на основе анализа нечеткого поля направлений

Выделение центральных линий производится на основе анализа значений функции Ж (а, т) и направлений Б(а,т) последовательно для каждого значения масштаба.

Рис.47. Примеры работы вейвлет метода оценки поля направлений для двух значений масштабной переменной а : а), б) - весовые функции поля направлений, в), г) -угловые функции поля направлений

Рис.48. Примеры центральных линий на различных масштабах и результат их комбинации: а), г) - для одного значения масштаба а; б), д) - для другого значения масштаба а; в), е) - комбинирование центральных линий; верхний и нижний ряды изображений относятся к двум различным исходным изображениям

При фиксированном значении а для каждой точки т0 производится поиск ближайшего локального максимума ттах вдоль прямой, проходящей через данную точку перпендикулярно направлению Б(а,т0). Если найденный локальный максимум от-

стоит от исходной точки на один или менее пикселей (Цт^ - т01| < 2), то данная точка считается точкой центральной линии. Далее производится фильтрация центральных линий по их протяженности и комбинирование центральных линий, обнаруженных на различных масштабах (рис.48).

Результирующий набор центральных линий используется для восстановления пространственной структуры коронарных сосудов по малому числу наблюдаемых ангиографических проекций.

5. Результаты восстановления 3Б структуры

Представленный подход трёхмерной реконструкции был применён для восстановления пространственной структуры фантомов (искусствен-ный объект, снятый в реальных условиях) с целью проверки работоспособности алгоритмов и точности восстановления, а также для восстановления пространственной структуры реальных сосудов сердца.

Для оценки различных геометрических параметров были использованы искусственные модельные изображения (таблица 1), а также два ви -

да фантомов, представленных в таблицах 2-5. Геометрические характеристики (такие как длины ветвей, углы разветвления, толщина ветвей) фан -томов были измерены заранее. Восстановление фантомных объектов производилось по ангиогра-фическим проекциям, полученным с разным ракурсом. Как правило, использовались 2-3 проек -ции для восстановления.

Фантом представляет собой некоторый пространственный объект, сделанный из непрозрачного для рентгеновского излучения материала. По своей структуре он приближен к сосудам сердца. Фантом имеет два разветвления и, соответственно, три ветви различные по диаметру. Съемка проекций фантома производится аналогично съемке реального сердца.

Таблица 1. Восстановление ространстеенной структуры модельного объекта

Модельный объект Результаты восстановления

о

Углы поворота камеры

а = 0,0° /? = 0,С

Базовая точка х = 86 II

£

Углы поворота камеры

а = 90,0° /3 = 0, 0°

Базовая точка х = 35 г = 140

Проекции восстановленного модельного объекта

По полученным проекциям необходимо восстановить пространственную структуру исходного объекта. Для всех проекций использовался один и тот же пространственный шаг трассировки, равный 4 пикселям. Следует отметить, что в процессе тестирования исследовалась сама возможность восстановления исходного объекта по данным проекциям. Большой акцент также сделан на проверку правильности геометрических преобразований.

Проводилось тестирование алгоритма обнаружения пространственных узлов, выбора опти-мальных направлений движения в узловых точках. Чёрными точками на рисунках полусферы обозначены точки минимума интенсивности. Как

видно, в точке разветвлений было найдено два истинных минимума. Оценка качества восстановле -ния пространственной структуры проводилась на основе оценивания различных геометрических параметров фантомов.

На фантомах первого вида измерялась точность восстановления разветвлений, длин ветвей и толщины ветвей. Тестирование производилось на 120 тройках различных проекций. Параметры исходного фантома приведены в таблице 3. В таблице 4 приведены относительные погрешности восстановления по указанным геометрическим параметрам. Относительные погрешности вычислялись как среднее по всей выборке троек проекций.

Таблица 2. Восстановление фантома по трём проекциям

Проекции фантома

Результаты восстановления

Углы поворота камеры а = 0,4° ß = 0,4° Базовая точка x = 122 y = 210

Углы поворота а =-10,4° ß = 30,5° Базовая точка x = 172 y = 252

Углы поворота камеры а = 17,4 ß = -43,9 Базовая точка x = 148 y = 248

1 узел 2 узел

Пространственная полусфера возможных направлений

Восстановленный фантом

Таблица 3. Геометрические параметры фантома первого вида

Фантом Название ветви Длина Толщина

Абсолютная Отно сительная Абсолютная Отно сительная

AI А2 ES A1 22 0,222 6,3 1

A2 20 0,202 6,3 1

A3 57 0,576 6,3 1

B 50 0,505 3,2 0,508

C 80 0,808 3,2 0,508

Угол AB Угол AC

45° 53°

Следует отметить, что наличие погрешности в основном связано с рассогласованием проекций, вызванных неполной компенсацией центрального проецирования, а также неточной информацией о геометрии съемки. Этим объясняется «укорачивание» тонких длинных ветвей фантома при восстановлении. Так как наша методика предполагает со-

гласование проекций по начальной точке сосуда (точке входа катетера в сосуд), то рассогласование проекций усиливается при удалении от этой точки. Заметим также, что погрешность восстановления углов между ветвями фантома в основном обуславливается неполной компенсацией центрального проецирования.

Таблица. 4. Погрешности восстановления фантома первого вида

Фантом Название ветви Относительная погрешность длин (%) Отно сительная погрешность толщины (%)

/]\А( A1 3,57 1,12

Л2 3,32 1,07

A3 2,07 0,42

В 2,57 1,92

С 3,81 2,01

Относительная погрешность по углу АВ (%) Относительная погрешность по углу АС (%)

1,26 1,83

Таблица 5. Погрешности восстановления фантома второго вида

Фантом

Название ветви Отно сительная погрешность (%)

Л 5,72

В 7,37

С 4,21

Б 5,60

Второй вид фантомов использовался для оценки точности восстановления изогнутых ветвей. Алгоритмы были протестированы на 50 различных тройках проекций. Для каждой изогнутой ветви исходного фантома было рассчитано отношение расстояния между начальной и конечной точками к длине всей ветви. Сравнение результатов производилось по этой величине с учетом того, что восстановленные ветви могут оказаться короче исходных (в относительных величинах). Ошибки восстановления усреднялись по всем восстановленным фантомам для соответствующих ветвей. Результаты приведены в таблице 5. Как и следовало ожидать, в случае изогнутых ветвей возрастает относительная погрешность восстановления длин кривых, что также связано с наличием эффекта центрального проецирования и отсутствием точной информации о геометрии съемки проекций. Технология восстановления была использована для реконструкции пространственной структуры реальных коронарных артерий. В этом случае использовались 3 проекции. Сравнительная оценка геометрических параметров не проводилась в силу отсутствия параметров исходного объекта (сердца пациента). Была проведена визуальная оценка правильности восстановления реальной трёхмерной структуры.

Результаты совмещения проекций и восстановления пространственной структуры представлены на

рисунках 49-52 (первая проекция: углы поворота камеры а = -7,2, в = -30,8, базовая точка х = 170, у = 50; вторая проекция: углы поворота а = 55,7, в = 45,0 , базовая точка х = 226,у = 264 ; третья проекция: углы поворота а = -11,9,в = 1,5 , базовая точка х = 74, у = 69).

Разрабатывая методику восстановления, мы стремились максимально автоматизировать процесс реконструкции пространственного объекта. Все приведенные примеры восстановления были получены полностью автоматически. Существует множество проблем, связанных с восстановлением реальных объектов и обусловленных динамичностью объекта исследования, наличием искажений, возникающих при регистрации проекций, неточными данными о геометрии съемки и т.п.

Все эти факторы приводят к рассогласованию проекций и значительному понижению точности реконструкции. Все это накладывает довольно жесткие условия на выбор исходных проекций.

Экспериментальные исследования показали, что при сильных рассогласованиях проекций минимум функций на полусфере возможных направлений смещается от истинного направления сосуда, что приводит к понижению точности восстановления объекта. Вследствие этого был использован альтернативный подход к выделению центральных линий сосудов, опи-

описанный выше, который показал работоспособность алгоритмов на рассогласованных проекциях.

Заключение Для реализации информационной технологии восстановления пространственной структуры коронарных сосудов по ангиографическим проекциям были решены следующие задачи.

Задача идентификации фазы сердца на изображениях проекций возникает из-за необходимости выбора для каждой проекции одного кадра из последовательности. Требуется выбрать такие кадры, чтобы разные проекции представляли один и тот же пространственный объект. Для решения этой задачи используется метод корреляционного анализа функции яркости изображения проекции.

Исследования показали, что критериальная функция может быть построена на основе коэф-

фициента корреляции соседних отсчетов для каждого кадра.

Задача компенсации эффекта центрального проецирования возникает из-за особенностей распространения рентгеновского излучения. С увеличением угла наклона камеры масштабы становятся различными на различных частях изображения, что в результате приводит к нарушению пропорций объекта. Разработан алгоритм компенсации эффекта центрального проецирования. Аналогичная задача связана с компенсацией влияния дисторсии.

Задача геометрической привязки изображений проекций заключается в нахождении начальной точки и начального направления для выполнения трассировки дерева сосудов. Производится автоматический поиск согласованных изображений катетера на проекциях.

Рис.49. Совмещение проекций: а) - первая проекция с базовой точкой, б) пространственное совмещение двух проекций

Рис.50. Экран компьютерной системы визуализации восстановленной пространственной структуры коронарных сосудов

Рис.51. Исходные изображения проекций (а) и восстановленная по ним пространственная структура (б)

Рис.52. Интерфейс пользователя системы Задача идентификации точек ветвления сосудов решается на основе анализа локальной пространственной интенсивности внутри сферической области вокруг исследуемой точки методом скользящего сектора (лучевого пространственного преобразования). Используются статистические параметры сфе-

восстановления коронарных сосудов рической области для определения, в каком направлении исходит сосуд из центральной точки сосуда.

Задача пространственной трассировки сосудов является основной для восстановления пространственного дерева сосудов. Решение задачи заключается в построении трассы сосудов в трехмерном пространстве

на основе изображения проекций исходного объекта, полученных в различное время при условии, что объект может со временем изменять свою форму в некоторых пространственных пределах. Для компенсации динамических искажений производится согласование проекций на каждом шаге трассировки.

При реализации технологии использовались новые методы анализа и обработки изображений, имеющие фундаментальный характер, которые могут быть использованы для решения широкого класса научно-технических и практических задач, не связанных с медицинской тематикой. К таким методам относятся, например, метод многомерного поля направлений и методы обработки изображений, заданных на сферической поверхности, метод пространственного веерного преобразования [33], метод нечеткого поля направлений [27].

Разработанная компьютерная технология восстановления пространственной структуры коронарных сосудов по данным рентгеновской ангиографии позволит создать коммерческую компьютерную систему (рис. 50-52) для использования во врачебной практике для диагностики и лечения сердечных заболеваний.

Благодарности

Работа выполнена при поддержке российско-американской программы «Фундаментальные исследования и высшее образование» (BRHE), гранта Президента РФ поддержки ведущих научных школ (НШ-3086.2008.9) и программы Президиума РАН «Фундаментальные науки - медицине».

Литература

1. Coatrieux, JL. Computer Vision Approaches for the Three Dimensional Reconstruction of Coronary Arteries: Review and Prospects / JL. Coatrieux, M Garreau, R. Collorec, C. Roux // Critical Reviews in Biomedical Engineering. - 1994.

- 22(1). - P. 1-38.

2. Digital Imaging and Communications in Medicine (DICOM)/ Standart PS 3.1-2008 // National Electrical Manufacturers Association. - Rosslyn, Virginia, USA. - 2008.

3. Pianykh, Oleg S. Digital Imaging and Communications in Medicine (DICOM) / Oleg S. Pianykh // A Practical Introduction and Survival Guide. - Springer; Visit. - 2008.

4. Wang, Y. Vessel extraction in coronary X-ray Angiography / Y. Wang; C. Toumoulin, H.Z. Shu, Z.D. Zhou; J.L. Coatrieux // 27th Annual International Conference of the Engineering in Medicine and Biology Society. - 2005. - P.1584 - 1587.

5. Schrijver, M. АРМ Angiographic Assessment of Coronary Stenoses: A Review of the Techniques / M Schrijver // Archives Of Physiology And Biochemistry. - 2003. - Vol. 111, No. 2. - P. 77-158.

6. Ильясова, Н.Ю. Методы и алгоритмы оценивания геометрических параметров диагностических изображений / Н.Ю. Ильясова. Диссертация на соискание учёной степени кандидата технических наук - Самара.

- 1997. - 155 с.

7. Методы компьютерной обработки изображений (издание второе, исправленное) / под ред. В.А Сойфера - М.: Физматлит, 2003. - 780 c.

8. Yang, Ching-Wen. Computer-Aided Diagnostic Detection System of Venous Beading in Retinal Images / Ching-Wen Yang, Dye-Jyun Ma, Shuenn-Ching Chao, Chuin-Mu Wang,

Chia-Hsin Wen, Chien-Shun Lo, Pau-Choo Chung, and Chein-I Chang // Opt. Eng. - 2000. - 39 (5). - P. 1293-1303.

9. Mourgues, F. 3D+t Modeling of Coronary Artery Tree from Standard Non Simultaneous Angiograms/ F. Mourgues, F. Devernay, G. Malandain and E. Coste-Maniere // Medical Image Computing and Computer-Assisted Intervention Proceedings, Lecture Notes in Computer Science Springer, Utrecht. - 2001. - V. 2208.

10. Messenger, J.C. 3D coronary reconstruction from routine single-plane coronary angiograms: clinical validation and quantitative analysis of the right coronary artery in 100 patients / J.C. Messenger, S.Y. Chen, J.D. Carroll, J.E. Burche-nal, K. Kioussopoulos and B.M. Groves // The International Journal of Cardiac Imaging. - 2000. - N 16 (6). - P. 413-427.

11. Ding, Z Quantification of 3-D coronary arterial motion using clinical biplane cineangiograms / Z. Ding and M.H. Friedman // The International Journal of Cardiac Imaging. - 2000. -N 16 (5). - P. 331-346.

12. Kupriyanov, A.V. Method of the Functional Tag Field in the Analysis of Ocular Fundus Images/ A.V. Kupriyanov, A.M. Malafeev // Pattern Recognition and Image Analysis. -2003. - Vol. 13, No. 2. - P. 287-289.

13. Ilyasova, N.Yu. A Method for Determination of an Optimal Spatial Direction of Vessels in the Problem of Reconstructing the 3D Topology of a Coronary System / N.Yu Ilyasova, A.V. Kupriyanov, A.G. Khramov, A.O. Korepanov // Pattern Recognition and Image Analysis. - 2003. - Vol. 13, No. 2. -P. 287-289.

14. Ilyasova, N.Yu. Estimating the Geometrical Parameters of Branches of 3D Treelike and Netlike Structures Exemplified by the Cardiovascular Image/ N.Yu. Ilyasova, A.V. Ustinov, V.G. Baranov, A.V. Kupriyanov and A.G. Khramov // Optical Memory & Neural Networks. - 2002. - Volume 11, N. 1. - P. 65-72

15. Soifer, V.A Estimating the geometrical parameters of branches of treelike and netlike structures / V.A.Soifer, N.Yu. Ilyasova, V.G. Baranov, A.G. Khramov // Proceedings of the 12th Scandinavian Conference on Image Analysis, Bergen, Norway, June 11-14. - 2001. - P. 49-53.

16. Ilyasova, N.Yu. An Expert Computer System for Diagnosing Eye Diseases from Retina Images/ N.Yu. Ilyasova, A.V.Ustinov, V.G.Baranov // Optical Memory & Neural Networks. - 2000. - V.9, N.2. - P.113 -120.

17. Ilyasova, N.Yu. Measuring the geometric parameters using image processing and diffractive optics methods / N.Yu. Ilyasova, V.A. Soifer, V.V. Kotlyar, S.N. Khonina, A.G. Khramov // Seventh International Symposium on Laser Metrology Applied to Science, Industry, and Everyday Life. Proceedings of SPIE. - 2002. - V. 4900. - P. 996-1006.

18. Ilyasova, N.Yu. Measuring Biomechanical Characteristics of Blood Vessels for Early Diagnostics of Vascular Retinal Pathologies / N.Yu. Ilyasova, A.V. Kupriyanov, MA Ananin, N.A Gavrilova // Proceedings of 7th International Conference Medical Image Computing and Computer Assisted Intervention. - MICCAI 2004, Saint-Malo, France, September, Part II. - 2004. - P. 251-258.

19. Ilyasova, N.Yu. Estimation of Geometrical Diagnostic Features the Optical Disk Images/ N.Yu. Ilyasova, A.V. Kupriyanov, A.V. Ustinov, M.A. Ananin // Pattern Recognition and Image Analysis. - MAIK «Nauka/Interperiodica». - 2005. -V. 15, N. 3. - P. 602-605.

20. Anan'in, M.A. Estimating Directions of Optic Disk Blood Vessels in Retinal Images / M.A. Anan'in, N.Yu. Ilyasova, A.V. Kupriyanov // Pattern Recognition and Image Analysis. - MAIK «Nauka/Interperiodica». - 2007. - V.17, N. 4. - P. 523-526.

21. Ilyasova, N.Yu. Methods for formation of features of treelike structures on fundus images / N.Yu. Ilyasova, V.V. Yatul'chik // Pattern Recognition and Image Analysis. -

MAIK Nauka/Interperiodica distributed exclusively by Springer Science+Business Media, Inc. - 2006. - V. 16, N. 1.

- P. 124-127. - ISSN: 1054-6618.

22. Udupa, J.K 3D Imaging in Medicine. Second ed. / J.K. Udupa, G.T. Herman // CRC Press. - 2000. - 384 p.

23. Armande, N. Thin nets extraction using multi-scale approach / N. Armande, P. Montesinos, and O. Monga // Computer Vision and Image Understanding. - 1999. - V. 73. - P. 248-257.

24. Devernay, F. 3D+t Modeling of Coronary Artery Tree from Standard Non Simultaneous Angiograms / F. Devernay, F. Mourgues, G. Malandain, and E. Coste-Maniere // in Proc. on MICCAI. - 2001. - P. 1320-1322.

25. Poli, R. An algorithm for real-time vessel enhancement and detection / R. Poli and G. Valli // Comp. Methods and Prog. in Biomed - 1997. - V. 52. - P. 1-22.

26. Korepanov, A.O. Extraction of Central Lines in Diagnostic Images of Vessels Using Methods of Wavelet Analysis and Differential Geometry / A.O. Korepanov // Pattern Recognition and Image Analysis. - 2005. - V. 15, N. 4. - P. 636-639.

27. Soifer, V.A. Fuzzy Direction Field Method for Fringe and Tree-like Patterns Analysis / V.A. Soifer, A.G. Khramov, A.O. Korepanov // in Proc. on 17th International Conference on Pattern Recognition. - 2004.

28. Stockman, G. Computer Vision / G. Stockman & L. Shapiro, // Prentice Hall. - 2001.

29. Ilyasova, N.Yu. Geometrical Parameters Estimation of the Retina Images for Blood Vessels Pathology Diagnostics / N.Yu Ilyasova, A.V. Kupriyanov, M.A Ananin // Proceedings of 15th European Signal Processing Conference. -EUSIPCO 2007, Poznan, Poland, September 3-7. - 2007. -P. 1251-1254.

30. Ильясова, Н.Ю Применение искусственных нейронных сетей для оценивания диагностических параметров на биомедицинских изображениях / Н.Ю. Ильясова, Д.Е. Липка, А.В. Куприянов // Компьютерная оптика. -2003. - № 25. - C. 151-154.

31. Ковалев, А. А. Восстановление связности линий на бинарных изображениях древовидных структур / Н.Ю. Ильясова, АА. Ковалев, А.В. Куприянов, А.Г. Храмов // Компьютерная оптика. - 2002. - № 23. - C. 58-61.

32. Сойфер, В. А. Методы компьютерного анализа диагностических изображений глазного дна / В.А. Сойфер, Н.Ю. Ильясова, А.В. Куприянов, АГ. Храмов, М.А Ананьин, // Технологии живых систем. - издательство «Радиотехника». - 2008. - T. 5, № 5-6. - C.61-71.

33. Куприянов, А. В. Определение параметров русла сосудов с использованием трехмерного локального веерного преобразования / А.В. Куприянов, Н.Ю. Ильясова, А.Г. Храмов, А. М. Малафеев, О. А. Титова // Компьютерная оптика. - 2004. - № 25. - C. 154-157.

34. Корепанов, А. О. Метод выделения центральных линий кровеносных сосудов на диагностических изображениях/ Н.Ю. Ильясова, П.М. Чикулаев // Компьютерная оптика.

- 2006. - № 29. - C. 146-151.

35. Ильясова, Н.Ю Математические модели и методы оценивания диагностических параметров древовидных структур: учебное пособие / Н.Ю. Ильясова, А.В. Куприянов, М.А. Ананьин //Изд-во Самар. гос. аэрокосм. ун-та. - 2007. - 104 с.

36. Ильясова, Н.Ю Оценивание геометрических параметров биомедицинских диагностических изображений: учебное пособие / Н.Ю. Ильясова, А.В. Куприянов, А. О. Корепанов // Изд-во Самар. гос. аэрокосм ун-та. - 2007. - 142 с.

37. Сойфер, В.А. Восстановление пространственной структуры древовидных объектов на основе нечеткого поля направлений: учебное пособие / В.А. Сойфер, А.О. Коре-панов, Н.Ю. Ильясова, А.В. Куприянов // Изд-во Самар. гос. аэрокосм. ун -та - 2007. - 64 с.

38. Сойфер, В. А. Методы обработки и анализа данных рентгеновской томографии : учебное пособие / В.А Сойфер, Н.Ю. Ильясова, А.О. Корепанов, А.В. Куприянов // Изд-во Самар. гос. аэрокосм. ун-та. - 2007. - 112 с.

39. Ильясова, Н.Ю. Оценивание диагностических параметров сосудов на изображениях глазного дна в области диска зрительного нерва / Н.Ю. Ильясова, АВ. Куприянов, МА.Ананьин // Компьютерная оптика. - 2006. -№ 29. - С. 141-146.

40. Ильясова, Н.Ю. Экспертная компьютерная система диагностики глазных заболеваний по изображениям глазного дна / Н.Ю. Ильясова, В. Г. Баранов, А.В. Устинов // Компьютерная оптика. - 1999. - № 19. - С. 202-209.

41. Ильясова, Н.Ю. Измерение биомеханических характеристик сосудов для ранней диагностики сосудистой патологии глазного дна / Н.Ю. Ильясова, А.В. Куприянов, М.А. Ананьин // Компьютерная оптика. -2005. - № 27. - С. 165-170.

42. Ильясова, Н.Ю. Оценивание геометрических параметров ветвей трёхмерных древовидных и сетчатых структур на примере изображения сосудистой системы сердца / Н.Ю. Ильясова, А.О. Корепанов, А.В. Куприянов, В.Г. Баранов, А.Г. Храмов// Компьютерная оптика - 2002. - № 23. - С. 48-52.

43. Ильясова, Н.Ю. Анализ структуры сосудистой системы сердца методом трассировки изображений проекций / Н.Ю. Ильясова, А.О. Корепанов, АВ. Куприянов, В.Г. Баранов, А.Г. Храмов // Компьютерная оптика. -2002. - № 23. - С. 53-57.

44. Корепанов, АО. Метод определения оптимального пространственного направления сосудов в задаче восстановления 3D топологии коронарной системы / А.О. Корепанов, Н.Ю. Ильясова, А.В. Куприянов, А.Г. Храмов // Компьютерная оптика. - 2002. - № 24. - С. 152-154.

45. Ильясова, Н.Ю. Метод пространственного восстановления коронарных артерий по малому числу ангиографиче-ских проекций / Н.Ю. Ильясова, А.О. Корепанов, АВ. Куприянов, А.В. Устинов, А.Г. Храмов, А.О. Ковалёв // Компьютерная оптика. - 2004. - № 26. - С. 61-72.

46. Ильясова, Н.Ю. Оценивание геометрических параметров области диска зрительного нерва на изображении глазного дна / Н.Ю. Ильясова, А.В. Куприянов, М.А Ананьин, АМ. Малафеев, А.В. Устинов // Компьютерная оптика. - 2005. - № 28. - С. 136-139.

47. Computer Image Processing, Part I: Basic concepts and theory / edited by Victor A Soifer // VDM Verlag . - 2009. -296 p.

48. Computer Image Processing, Part II: Methods and algorithms / edited by Victor A. Soifer // VDM Verlag . - 2009. - 584 p.

References

1. Coatrieux, JL. Computer Vision Approaches for the Three Dimensional Reconstruction of Coronary Arteries: Review and Prospects / JL. Coatrieux, M. Garreau, R. Collorec, C. Roux // Critical Reviews in Biomedical Engineering. -1994. - 22(1). - P. 1-38.

2. Digital Imaging and Communications in Medicine (DICOM)/ Standart PS 3.1 -2008 // National Electrical Manufacturers Association. - Rosslyn, Virginia, USA. - 2008.

3. Pianykh, Oleg S. Digital Imaging and Communications in Medicine (DICOM) / Oleg S. Pianykh // A Practical Introduction and Survival Guide. - Springer; Visit. - 2008.

4. Wang, Y. Vessel extraction in coronary X-ray Angiography / Y. Wang; C. Toumoulin, H.Z. Shu, Z.D. Zhou; J.L. Coatrieux // 27th Annual International Conference of the Engineering in Medicine and Biology Society. - 2005. - P.1584 - 1587.

5. Schrijver, M. АРМ Angiographic Assessment of Coronary Stenoses: A Review of the Techniques / M. Schrijver //

Archives Of Physiology And Biochemistry. - 2003. - Vol. 111, No. 2. - P. 77-158.

6. Ilyasova, N. Yu. Methods and Algorithms for the Estimation of the Geometrical Parameters of the Diagnostic Images / N. YU. Ilyasova // Dissertation for candidate's degree application - Samara. - 1997. - 155 p. (in Russian).

7. Methods of Computer Image Processing (Secondary Edition)/ edited by V.A. Soifer - Мoscow: Fizmatlit, 2003. - 780 p. (in Russian)

8. Yang, Ching-Wen. Computer-Aided Diagnostic Detection System of Venous Beading in Retinal Images / Ching-Wen Yang, Dye-Jyun Ma, Shuenn-Ching Chao, Chuin-Mu Wang, Chia-Hsin Wen, Chien-Shun Lo, Pau-Choo Chung, and Chein-I Chang // Opt. Eng. - 2000. - 39 (5). - P. 1293-1303.

9. Mourgues, F. 3D+t Modeling of Coronary Artery Tree from Standard Non Simultaneous Angiograms/ F. Mourgues, F. Devernay, G. Malandain and E. Coste-Maniere // Medical Image Computing and Computer-Assisted Intervention Proceedings, Lecture Notes in Computer Science Springer, Utrecht. - 2001. - V. 2208.

10. Messenger, J.C. 3D coronary reconstruction from routine single-plane coronary angiograms: clinical validation and quantitative analysis of the right coronary artery in 100 patients / J.C. Messenger, S.Y. Chen, J.D. Carroll, J.E. Burchenal, K. Kioussopoulos and B.M. Groves // The International Journal of Cardiac Imaging. - 2000. - N 16(6). - P. 413-427.

11. Ding, Z. Quantification of 3-D coronary arterial motion using clinical biplane cineangiograms / Z.Ding and M.H.Friedman // The International Journal of Cardiac Imaging. - 2000. - N 16(5) - P. 331-346.

12. Kupriyanov, A.V. Method of the Functional Tag Field in the Analysis of Ocular Fundus Images/ A.V.Kupriyanov, A.MMalafeev // Pattern Recognition and Image Analysis. -2003. - Vol. 13, No. 2. - P. 287-289.

13. Ilyasova, N.Yu. A Method for Determination of an Optimal Spatial Direction of Vessels in the Problem of Reconstructing the 3D Topology of a Coronary System/ N.YuIlyasova, AV.Kupriyanov, AG.Khramov, AO.Korepanov// Pattern Recognition and Image Analysis. - 2003. - Vol. 13, No.2. - P. 287-289.

14. Ilyasova, N.Yu. Estimating the Geometrical Parameters of Branches of 3D Treelike and Netlike Structures Exemplified by the Cardiovascular Image/ N.Yu.Ilyasova, A.V.Ustinov, V.G.Baranov, A.V.Kupriyanov and A.G.Khramov// Optical Memory & Neural Networks. -2002. - Volume 11, N. 1. - P. 65 -72

15. Soifer, V.A. Estimating the geometrical parameters of branches of treelike and netlike structures/ V.A.Soifer, N.Yu Ilyasova, V.G. Baranov, A.G. Khramov// Proceedings of the 12th Scandinavian Conference on Image Analysis, Bergen, Norway, June 11-14. - 2001. - P.49-53.

16. Ilyasova, N.Yu. An Expert Computer System for Diagnosing Eye Diseases from Retina Images/ N.Yu. Ilyasova, A.V.Ustinov, V.G.Baranov// Optical Memory & Neural Networks. - 2000. - V.9, N.2. - P.113 -120.

17. Ilyasova, N.Yu. Measuring the geometric parameters using image processing and diffractive optics methods/ N.Yu. Ilyasova, V.A.Soifer, V.VKotlyar, S.N.Khonina, AG.Khra-mov // Seventh International Symposium on Laser Metrology Applied to Science, Industry, and Everyday Life. Proceedings of SPIE. - 2002. - V.4900. - P.996-1006.

18. Ilyasova, N.Yu. Measuring Biomechanical Characteristics of Blood Vessels for Early Diagnostics of Vascular Retinal Pathologies/ N.Yu.Ilyasova, A.V.Kupriyanov, M.A.Ananin, N.A.Gavrilova// Proceedings of 7th International Conference Medical Image Computing and Computer Assisted Intervention. - MICCAI 2004, Saint-Malo, France, September, Part II. - 2004. - P. 251-258.

19. Ilyasova, N.Yu. Estimation of Geometrical Diagnostic Features the Optical Disk Images/ N.Yu. Ilyasova, A.V. Kupriyanov, A.V. Ustinov, M.A. Ananin // Pattern Recognition and Image Analysis. - MAIK «Nauka/tnterperiodica». - 2005. -V. 15, N.3. - P. 602-605.

20. Anan'in, M.A. Estimating Directions of Optic Disk Blood Vessels in Retinal Images/ MAAnan'in, N.YuIlyasova , A.V.Kupriyanov// Pattern Recognition and Image Analysis. -MAIK «Nauka/Interperiodica». - 2007. - V.17, N. 4. - P. 523-526.

21. Ilyasova, N.Yu. Methods for formation of features of treelike structures on fundus images/ N.Yu. Ilyasova, V.V.Yatul'chik., // Pattern Recognition and Image Analysis. -MAIK Nauka/Interperiodica distributed exclusively by Springer Science+Business Media, Inc. - 2006. - V.16, N.1. -P. 124-127. - ISSN: 1054-6618.

22. Udupa, J.K. 3D Imaging in Medicine. Second ed. / J.K. Udupa, G.T.Herman// CRC Press. - 2000. - 384 p.

23. Armande, N. Thin nets extraction using multi-scale approach / N. Armande, P. Montesinos, and O. Monga // Computer Vision and Image Understanding. - 1999. - V. 73. - P. 248-257.

24. Devernay, F. 3D+t Modeling of Coronary Artery Tree from Standard Non Simultaneous Angiograms / F. Devernay, F. Mourgues, G. Malandain, and E. Coste-Maniere// in Proc. on MICCAI. - 2001. - P. 1320-1322.

25. Poli, R. An algorithm for real-time vessel enhancement and detection / R. Poli and G. Valli // Comp. Methods and Prog. in Biomed. - 1997. - V. 52. - P. 1-22.

26. Korepanov, A.O. Extraction of Central Lines in Diagnostic Images of Vessels Using Methods of Wavelet Analysis and Differential Geometry/ A.O. Korepanov//, Pattern Recognition and Image Analysis. - 2005. - V. 15, N. 4. - P. 636-639.

27. Soifer, V.A. Fuzzy Direction Field Method for Fringe and Tree-like Patterns Analysis/ V.A. Soifer, A.G. Khramov, A.O. Korepanov// in Proc. on 17th International Conference on Pattern Recognition. - 2004.

28. Stockman, G. Computer Vision/ G.Stockman & L.Shapiro, //Prentice Hall. - 2001.

29. Ilyasova, N.Yu. Geometrical Parameters Estimation of the Retina Images for Blood Vessels Pathology Diagnostics/ N.Yu. Ilyasova, A.V. Kupriyanov, MA Ananin// Proceedings of 15th European Signal Processing Conference. - EUSIPCO 2007, Poznan, Poland, September 3-7. - 2007. - P. 1251-1254.

30. Ilyasova, N.Yu. Application of the artificial neural networks to estimate the diagnostic parameters of the biomedical images / N.Yu. Ilyasova, D.E. Lipka, A.V. Kupriyanov // Computer optics. - 2003. - N 25. - P. 151-154. - (in Russian).

31. Kovalyov, A.A. Fractures elimination on the binary images of the treelike structures/ N.YuIlyasova, A.A.Kovalyov, A.V.Kupriyanov, A.G.Khramov // Computer Optics. - 2002.

- N 23. - P. 58-61. - (in Russian).

32. Soifer, V.A. Methods for Computer Diagnostics using Eye's Fundus Images / V.A. Soifer, N.Yu. Ilyasova, A.V. Kupriyanov, A.G. Khramov, M.A. Ananin // Technologies of the living systems. - Radiotechnika. - 2008. - V. 5, N 5-6. - P. 61-71. - (in Russian).

33. Kupriyanov, A.V. Estimation of the vessel route parameters with local 3D fan-beam transformation / AV.Kupriyanov, N.Yu. Ilyasova, A.G. Khramov, AM Malapheev, O.A. Ti-tova // Computer optics. - 2004. - N 25. - P. 154-157. - (in Russian).

34. Korepanov, A.O. Methods to extract the central lines of the blood vessel upon the diagnostic images / N.Yu. Ilyasova, A.O. Korepanov, P.M. Chikulaev // Computer optics. - 2006.

- N 29. - P. 146-151. - (in Russian).

35. Ilyasova, N.Yu. Mathematical models and method for estimation of the diagnostic parameters of the tree-like structures: methodical guidelines / N.Yu. Ilyasova, A.V. Kupriyanov, M.A. Ananin // SSAU publishing. - 2007. - 104 p. - (in Russian)

36. Ilyasova, N.Yu. Estimation of the geometrical parameters of the biomedical diagnostic images: methodical guidelines / N.Yu. Ilyasova, A.V. Kupriyanov, A.O. Korepanov // SSAU publishing. - 2007. - 142 p. - (in Russian)

37. Soifer, V.A. Reconstruction of the spatial structure of the tree-like objects base on the fuzzy direction field: methodical guidelines / V.A. Soifer, N.Yu. Ilyasova, A.V. Kupriyanov, A.O. Korepanov // SSAU publishing. - 2007. - 64 p. - (in Russian)

38. Soifer, V.A. Methods for the processing and analysis of the x-ray tomography data: methodical guidelines / V. A. Soifer, N.Yu. Ilyasova, A.V. Kupriyanov, A.O. Korepanov // SSAU publishing. - 2007. - 112 p. - (in Russian)

39. Ilyasova, N.Yu. Estimation of the diagnostic parameters of vessels upon the eye fundus images in the optic disk region / N.YuIlyasova, A.V.Kupriyanov, M.A Ananin // Computer Optics. - 2006. - N29. - P. 141-146. - (in Russian).

40. Ilyasova, N.Yu. Expert computer system for the diagnostics of the eye diseases on the retina images / N.Yu. Ilyasova, V.G. Baranov, A.V.Ustinov // Computer Optics. - 1999. - N19. - P. 202-209. - (in Russian).

41. Ilyasova, N.Yu. Measurement of the biomechanical vessels parameter for the diagnostics of the early stages of the retina vascular pathology / N.Yu.Ilyasova, A.V.Kupriyanov, M.A. Ananin // Computer Optics. -2005. - N27. - P. 165-170. - (in Russian).

42. Ilyasova, N.Yu. Estimation of the geometrical parameters of the tree-like and net-like structures branches exemplified the coronary vessels system images / N.Yu. Ilyasova, A.O. Korepanov, A.V. Kupriyanov, V.G. Baranov,

A.G. Khramov // Computer Optics. - 2002. - N23. - P. 48-52. - (in Russian).

43. Ilyasova, N.Yu. Analysis of the coronary vessels structure based on the method of projections tracing / N.Yu. Ilyasova, A.O. Korepanov, A.V. Kupriyanov, V.G. Baranov, A.G. Khramov // Computer Optics. - 2002. - N23. - P. 53-57. - (in Russian).

44. Korepanov, A.O. The method for the estimation of the optimal spatial vessel route direction if the problem of the 3d reconstruction of coronary vessels structure topology / A.O. Korepanov, N.Yu. Ilyasova, A.V. Kupriyanov, A.G. Khramov // Computer Optics. - 2002. - N24. - P. 152-154. - (in Russian).

45. Ilyasova, N.Yu. The method for the spatial reconstruction of the coronary arteries upon the small number of angiographic projections / N.Yu. Ilyasova, A.O. Korepanov, A.V. Kupriyanov, V.G. Baranov, A.V. Ustinov, A.G.Khramov, A.A. Kovalyov // Computer Optics. -

2004. - N26. - P. 61-72. - (in Russian).

46. Ilyasova, N.Yu. Estimation of the geometrical parameters of the optic disk region in the eye fundus images / N.Yu. Ilyasova, A.V. Kupriyanov, M.A. Ananin, A.M. Malapheev, A.V. Ustinov // Computer Optics. -

2005. - N28. - P. 136-139. - (in Russian).

47. Computer Image Processing, Part I: Basic concepts and theory / edited by Victor A. Soifer // VDM Verlag . -2009. - 296 p.

48. Computer Image Processing, Part II: Methods and algorithms / edited by Victor A. Soifer // VDM Verlag . - 2009. -584 p.

COMPUTER TECHNOLOGY FOR THE SPATIAL RECONSTRUCTION OF THE CORONARY VESELS STRUCTURE FROM ANGIGOGRAPHIC PROJECTIONS

Natalya Yurievna Ilyasova1'2 (senior researcher, e-mail: [email protected]), Nikolay Lvovich Kazansky1'(deputy director, e-mail: [email protected]), Andrey Olegovich Korepanov1 (docent, e-mail: [email protected]), Alexander Victorovich Kupriyanov12 (senior researcher, e-mail: [email protected]),

Andrey Vladimirovich Ustinov2 (lead. programmer, e-mail: [email protected]), Alexander Grigorievich Khramov12 (senior researcher, e-mail: [email protected]) 1 S.P. Korolyov Samara State Aerospace University, 2 Image Processing Systems Institute of the RAS

Abstract

The technology for the automatic reconstruction of the 3D structure of the coronary arteries from the small number of the angiographic projections is considered. The methods for construction of the structure are based on the sequential obtaining of the vessels route during the simultaneous analysis of the available projections. As result the vector representation of the spatial vessels is formed, which subsequently used for 3D structure visualization. In this review the stages of the reconstruction are thoroughly discussed and the results of some stages are presented. The method to obtain the spatial coordinates upon two image projections of the vessels route scanning is described. A set of the basic geometric characteristic is presented to describe objects of 2D and 3D structures. The approximation methods for the estimation of the vessels diameter and methods for the vessels central lines extraction, which are based on the wavelet transformations, are presented with full details. The experimental result provided from the test images are discussed, along with the results of the spatial structure reconstruction for the phantom objects and left coronary artery.

Key words: biomedical diagnostic images, coronary vessels, spatial structure reconstruction, geometrical characteristics, approximation methods, wavelet transformation.

Поступила в редакцию 12.05.2009 г