CC BY
58
ОЦЕНИВАНИЕ ГЕОМЕТРИЧЕСКИХ ПАРАМЕТРОВ ВЕТВЕЙ ТРЁХМЕРНЫХ ДРЕВОВИДНЫХ И СЕТЧАТЫХ СТРУКТУР НА ПРИМЕРЕ ИЗОБРАЖЕНИЯ СОСУДИСТОЙ СИСТЕМЫ СЕРДЦА
Н.Ю. Ильясова, А.В. Куприянов, А.В. Устинов, А.Г. Храмов, В.Г. Баранов Институт систем обработки изображений РАН Самарский государственный аэрокосмический университет
Аннотация
На примере сосудистой системы изображения сердца рассматривается класс изображений, содержащих ветви сетчатых и древовидных трёхмерных структур. Описывается один из методов получения 3-х мерных координат точек сканирования сосуда сердца по результатам обработки двух проекций изображения. Предлагается набор основных геометрических характеристик, описывающих объекты двухмерных и трёхмерных структур. Приводятся методы расчета предложенных параметров, а также алгоритмы оценивания признаков сосуда с учётом трёхмерности исходного участка в задаче анализа изображений сердца.
Введение
Целью настоящего исследования является разработка компьютерной системы для динамической визуализации и анализа пространственных объектов, структура которых описывается деревом. К таким структурам относится, в частности, кровеносная система человека - общая и отдельные её фрагменты (кровеносные сосуды в области сердца, кровеносная система глазного дна, и т.п.). Исходными данными для построения пространственной модели кровеносной системы является множество двумерных проекций, полученных, например, при рентгеноскопии. Предлагаемый нами метод пространственной трассировки позволяет восстанавливать пространственную структуру кровеносной системы и оценивать локальные параметры сосудов, необходимые в медицинских приложениях для определения степени патологии. На основе полученного структурного описания может быть легко построена трехмерная модель для динамической визуализации.
В работах [1, 2] мы рассматривали класс изображений, содержащих ветви сетчатых и древовидных двухмерных структур. Предлагали набор геометрических характеристик, описывающих объекты двухмерных структур. Рассматривали методы их расчета. В работе [3] оценивается один из параметров (четкообразность), но используется другой подход для выделения ветвей на изображениях. На-
стоящая работа является развитием методов и алгоритмов, предложенных ранее для плоского случая.
Локальными параметрами, получаемыми в результате трассировки ветви, как в двухмерном так и в трёхмерном случае, являются диаметр и направление ветви в каждой точке. Глобальные параметры (средний диаметр, прямолинейность, чёткообраз-ность, амплитуда колебаний толщины, частота колебаний толщины, извилистость толщины, частота колебаний трассы, амплитуда колебаний трассы, извилистость трассы) вычисляются на основе локальных. Ниже приведём определение данных параметров в двумерных структурах и их отличие для трёхмерного варианта. Введем в рассмотрение также специфический параметр объёмного изображения как «частота вращения», характеризующий такое пространственное свойство ветви, как отклонение трассы от плоскости.
Древовидная структура характеризуется наличием множества узлов, соединенных в виде направленного графа - дерева. Сетчатая структура образует ненаправленный граф на плоскости без пересечения дуг. В данной работе анализируются ветви, которые является общим элементом этих структур. Под ветвью мы понимаем протяженный объект между пересечениями или разветвлениями. Определяющими характеристиками ветви являются трасса (геометрическое место срединных точек), толщина и направление ветви в каждой точке.
щ
W Щ
: >
ж
а) б) в)
Рис. 1. Примеры изображений, характеризуемых наличием ветвей: а) кровеносная система глазного дна человека - древовидная структура;
б) кровеносная система сердца - древовидная структура;
в) сеть трещин в земле вдоль рек - сетчатая структура
Узлы структур и соединяющие их ветви могут располагаться как в пространстве, так и на плоскости.
Ранее мы рассматривали алгоритм трассировки ветвей и вычисления их глобальных параметров в предположении, что поле зрения плоское. В случае глазного дна (сосудистая система) или сетчатых структур на поверхности земли в рамках реальной точности вычислений это вполне приемлемо. Однако такое предположение не имеет места при анализе сосудистой системы объёмных органов (сердце). В этом случае необходимо учитывать трёхмерность наблюдаемого объекта, один из способов - использование проекций объекта на две различные плоскости.
Схема получения трёхмерной трассы
Будем считать, что объект исследования трёхмерной структуры снимается с двух точек зрения на плоскости, расположенные под углом а друг к другу. В каждом из двух проекционных изображений имеется своя система координат, но либо аппаратно, либо операциями сдвига и поворота можно достичь такой конфигурации, что обе координатные системы имеют общее начало координат и одну общую ось, совпадающую с линией пересечения проекционных плоскостей (рис. 2).
На рис. 2 изображён случай, когда совпадающей является ось Оу.
Пусть точка в трёхмерном пространстве имеет (х, у, 2 )
координаты 4 ' ', а на каждой из плоскостей
(х1, у1) и (х2, У2). Пространственные координаты выражаются через координаты на указанных проекциях следующим образом:
ч
у = У1 (= У2) .
2 = х2 (бш а + 008 ас/§а) -
(1)
Рис. 2. Расположение проекционных плоскостей и системы координат (вид сбоку)
Формула (1) выведена для непрерывных значений координат и применима для любых точек проекций. Особенность задачи анализа сосудов заключается в том, что обрабатываются координаты выделенного участка сосуда; которые принимают дискретные значения.
Ниже приведён метод получения трёхмерных координат трассы:
Рис. 3. Параллельная трассировка двух проекций
1. Производится стандартная (плоская) трассировка двух участков сосуда [4]- на первой проекции и соответствующего ему участка на второй проекции трёхмерного изображения (см. рис. 3). При этом 7-координата начальных и конечных точек участка трассы выбирается одинаковой.
2. Выбирается основная проекция и для всех точек трассы участка сосуда последовательно в соответствии с формулой (1) выполняем следующие присвоения:
a. x = Xj,
b. y' = yj, для точки (xj, yj) ищем ближайшую по Y-координате точку (x'2, y'2) трассы второго участка сосуда
c. z' = xJ2 (sinа+ cos actga) - xjctga
Если совпадающей является ось Ox, то выполняются аналогичные действия с использованием
следующей системы:
y = yj
x = xj (= x2) (2)
z = y2 (sina + cos actga) - yjctga
Вычисление признаков участка трёхмерного сосуда
Как указывалось в работе [1] ветвь сетчатой структуры можно описать двумя типами признаков: локальными (первичными) и глобальные, вычисляемые на основе локальных.
Локальные признаки включают в себя трассу (геометрическое место серединных точек), распределение диаметра и направления ветви вдоль трассы и рассчитываются непосредственно по изображению. Глобальные признаки, представленные ниже, характеризуют всю ветвь в целом. Математическая модель ветви древовидной структуры, представленная в [1] определяется функциями:
x = x(t), y = y(t), г = r (t), 0 < t < Lv ,
где x(t), y(t) - дифференцируемые функции, описывающие линию центров - трассу; г(t)- функция толщины ветви (расстояние от трассы до границы ветви, отсчитываемое по перпендикуляру к трассе); t - расстояние от начала трассы, измеренное по трассе; Lv - длина трассы. Данные характеристики однозначно определяют функцию направления трассы в каждой точке <p(t); функцию локальной высоты f (t), определяемую расстоянием от текущей точки трассы до ее проекции на отрезок L , соединяющий начальную и конечные точки трассы;
конфигурацию границ ветви - стенок xj = xj (t),
yj = yj (t), x2 = x2 (t), y2 = y2 (t), 0 < t < Lv .
И в плоском случае, и в объёмном локально сосуд можно считать участком цилиндра. Проекция цилиндра на плоскость, не перпендикулярную его
оси, является фигурой с поперечным размером, равным диаметру цилиндра (кроме концевых участков). Поэтому все признаки ветви, относящихся к толщине - средний диаметр, чёткообразность, амплитуда толщины, частота толщины, извилистость толщины (стенки), а так же параметр прямолинейности имеют одинаковое выражение как для двумерной, так и для трёхмерной структуры дерева изображения. Введём определение приведённых выше признаков ветви.
1. Средний диаметр
Средний диаметр ветви Б вычисляется по формуле:
- 2 N D = 2r = — I г (tn).
Nn=1
(3)
где = пА , N - число отсчётов локального радиуса, измеренных вдоль ветви в достоверных точках, Г -средний радиус ветви.
2. Прямолинейность
Прямолинейность ветви Р характеризует отклонение ветви от прямолинейного направления и определяется как отношение длины Ьу средней линии ветви к длине Ь прямой, соединяющей начальную и конечные точки трассы:
P = Lv = n-1 L
I (n -xn+1 )2 +(n -Уп+1)2 + ((n -zn+1)2]
(1 - xN )2 + (1- yN )2 +((1- zN )2 ¡2
.(4)
3. Чёткообразность
Чёткообразность ветви характеризует неравномерность толщины ветви и определяется как отношение среднеквадратичного отклонения радиуса ветви к его среднему значению:
S = V г2 - г2 / г
(5)
где г - средний радиус ветви, г - средний квадрат радиуса.
4. Амплитуда колебаний толщины Амплитуда колебаний толщины ветви А0 характеризует отклонения стенок ветви от прямой линии и определяется по формуле:
Ao =Л(2(г2 - г2).
(6)
5. Частота колебаний толщины Частота колебаний толщины ветви ю0 характеризует изменение направление стенок ветви на единицу ее длины и определяется по формуле:
N
\
2 R(m2) mj +— агctg-—
П R(mj ) J
R(m) =
,2nnm
I г (tn )e~
(7)
(8)
CO n =
0
N
n=0
m¡ = min(m0, m'), m2 = m¡ +1 = max(m0, m'),
(9)
m0 = argl max R(m) V 1<m<N
где m' - номер наибольшего из двух значений Фурье спектра, соседних с максимальным отсчётом, arg - аргумент вещественной функции f (x):
arg f (x) = x, r (t) - функция толщины ветви. 6. Извилистость толщины Извилистость толщины I0 характеризует скорость изменения функции толщины вдоль трассы аппроксимируемой гармонической функцией с амплитудой A0 и частотой ю0 и вычисляется по формуле:
10 = A0ra0 •
(10)
7. Извилистость трассы (для двухмерной структуры дерева)
Извилистость трассы 11 характеризует скорость изменения функции трассы (траектории ветви) на выделенном участке, которая аппроксимируется гармонической функцией А1 и частотой ю 1, по определению равна 11 = А1ю1 и определяется из уравнения:
P = 2 л/1 + Ii • E(к):
к = I1
(11)
(12)
1 +11
где Р - прямолинейность ветви, Е(к) - полный эллиптический интеграл 2-го рода:
E(к) = |л/1 - к2 sin21dt •
0
(13)
8. Амплитуда колебаний трассы (для двухмерной структуры дерева)
Амплитуда колебаний трассы А1 характеризует степень отклонения хода трассы от прямолинейного и определяется по формуле:
A = -
2f • E (к)
(14)
lnl I¡ +л/1 +1
1 +—-¡=
М1+11
где f - средняя высота ветви:
N
f = N-1 XI f (tn
n-1
(15)
где / ) - функция высоты ветви, 11 - извилистость трассы.
9. Частота колебаний трассы (для двухмерной структуры дерева)
Частота колебания участка ветви (трассы) ю1 характеризует число изменений направления трассы на единицу длины ветви и определяется по формуле:
Ij_
А
(16)
Рассмотрим определение признаков трассы ветви (извилистость, амплитуды колебаний и частота колебаний трассы) для трёхмерной структуры дерева. А также введём ещё один признак, характеризующий отклонение трассы от плоскости. Будем называть его частотой вращения. В качестве пространственной модели ветви рассмотрим участок спирали:
y = R sin rorx z = Rcosro,.x '
где гог - частота вращения; R - расстояние до оси Ох, совмещенной с прямой, соединяющей начальные и конечные точки трассы.
В свою очередь R является гармонической функцией: R = A sin rox, где А - амплитуда трассы; го - частота трассы.
В результате получим параметрическое уравнение трассы:
x = t
y = A sin rot sin ra rt • z = A sin rat cos ra rt
(17)
В системе имеются два параметра: А - амплитуда трассы и ш - частота трассы, которые аналогичны соответствующим параметрам в двумерном случае и один параметр ш r специфичный для трехмерной модели - частота вращения трассы.
Введём следующую функцию расстояния точки трассы от оси Ox:
r( x)=V7(X)+7(X)
Учитывая (17), получаем:
R(x) = A sin ox| .
Согласно этому уравнению среднее значение R = 2 A / п, откуда имеем выражение для амплитуды трассы:
A = nR / 2. (18)
Значение R(x) для каждой точки трассы вычисляется с помощью известной формулы для расстояния от точки до прямой [5]. Выражение (18) даёт приближённое значение параметра амплитуды. В [4] приведена точная формула оценивания параметра в плоском случае, которая при малом значении Аш , принимает вид (18). В силу того, что после перехода от двух проекций к трёхмерным координатам точки трассы не лежат равномерно по какой-либо величине, мы не сможем вывести точную фор-
1
п
мулу. Основываясь на результатах в плоском случае можно показать, что поправка к (18) невелика - в плоском случае при Лю ^ 0 получаем (18), а при Лю ^ да получаем Л = 2R (см. [4]).
Без использования преобразования Фурье в нашем случае невозможно вычислить частоту трассы. Функция R(х) = Л\sin юх| является периодической, но не является синусоидой, поэтому для повышения точности используем следующий подход: выполняем преобразование Фурье от центрированной функции R2 (x), где R2 (x) = 2Л2 (1 - cos 2юх) и получаем оценку частоты трассы:
ю = 1 arg max 3(R (х)),
(19)
Частоту вращения трассы теоретически можно найти аналогично частоте трассы:
ю r = arg max 3(sin ф(x)),
где ф(x) - полярный угол проекции точки трассы на плоскость yOz . Однако при таком подходе неизбежна большая погрешность оценивания в силу малого размера анализируемого участка сосуда и малого его отклонения от плоскости (точка проекции совершит заметно меньше полного оборота вокруг начала координат). В то же время спектральный анализ даёт приемлемые результаты, если имеется не менее одного периода разлагаемой функции. Используем следующий подход: для кривой (17) значение прямолинейности P равно:
P = + Л2 ю 2E(k),
(20)
полный эллиптиче-
,2 Л2(ю2-ю2г) где к2 = —Ь——^ и Е(-)
1+Л2ю2
ский интеграл 2-го рода. (Очевидно, что при ю г = 0 получаем формулу для плоского случая [4]).
Используя найденные значения величин Р , Л и ю, из (20) можно получить оценку частоты вращения трассы:
¡2 -
k = E
Л2
nP
у2-Jl + Л2ю2 ,
(21)
где E (•) - функция, обратная к эллиптическому интегралу 2-го рода.
Заключение
В результате проведённых исследований был разработан ряд методов и алгоритмов трассировки и оценивания геометрических параметров изображений плоских и пространственных древовидных и сетчатых структур.
Основной областью применения разработанных методов оценивания параметров являются биомедицинские задачи, связанные с анализом и измерением характеристик кровеносных систем: сосудов глазного дна (плоский вариант), сосудов сердца (трёхмерная структура), а также ряд геодезических задач (анализ изображений русел рек, оврагов). На основе разработанных методов была построена опытная компьютерная система для измерения геометрических параметров биомедицинских изображений [4]. Внедрение разработанных методов в медицинскую практику расширяет возможности существующих медицинских методик, позволяет автоматизировать диагностику ряда заболеваний, осуществить мониторинг патологических изменений на основе объективных количественных данных.
Другим применением разработанных методов является их использование в мультимедийных и Web-приложениях, ориентированных на визуализацию пространственных древовидных структур, например, кровеносной системы человека. Такие приложения могут быть использованы при обучении студентов-медиков и в справочных системах.
Литература
1. S.L. Branchevsky, A.B. Durasov, N.Yu. Ilyasova, A.V. Ustinov Methods for estimating geometric parameters of retinal vessels using diagnostic images of fundus // Proceedings SPIE, 1998. Vol.3348. P. 316-325.
2. S.L. Brantchevsky, Yu.V. Vasiliev, A.B. Durasov, N.Yu. Ilyasova, A.V. Ustinov, Method for the distinguishing and quantitative evaluation of the elements of pathological patterns in the retina (pathology of microciculation) // Proceedings SPIE, Vol. 2363. P. 236-242.
3. Ching-Wen Computer-aided diagnostic detection system of venous beading in retinal images // Optical Engineering, 2000. Vol.39, N. 5.Р. 1293-1303.
4. Ильясова Н.Ю., Устинов А.В., Баранов В.Г. Экспертная компьютерная система диагностики глазных заболеваний по изображениям глазного дна // Компьютерная оптика, 1999. №19. С. 202209.
5. Бронштейн И. Н., Семендяев К.А. Справочник по
математике, 1956.
п
юг =
1
