Математические модели адаптивной фильтрации биомедицинских изображений Текст научной статьи по специальности «Медицинские технологии»

2. Омельченко ВМ. Компьютерный анализ биопотенциалов мозга как основа оценки и фармакологической коррекции психопатологических состояний: Дис... д-ра биол. наук.

- Ростов-на-Дону, 1990. - 408 с.

3. Рудкоеский М.В., Омельченко ВМ., Матуа СМ. Дискретный электроэнцефалографиче-ский мониторинг фармакотерапии психоневротических больных с использованием метода многомерного шкалирования // Изд. вузов Сев.-Кавк. регион. Естественные науки.

- Ростов-на-Дону. - 2003. - № 8. - С. 59-67.

Омельченко Виталий Петрович

ГОУ ВПО «Ростовский государственный медицинский университет Росздрава». E-mail: [email protected]. г. Ростов-на-Дону, пер. Нахичеванский, 29.

.: 88632632352.

Omelchenko Vitaliy Petrovich

Rostov State Medical University.

E-mail: [email protected].

29, Nakhichevansky side street, Rostov-on-Don, Russia.

Phone: +78632632352.

УДК 621.396.1.001.24, 681.323:621.391

ЕЛ. Попечителев, И.В. Разин

МАТЕМАТИЧЕСКИЕ МОДЕЛИ АДАПТИВНОЙ ФИЛЬТРАЦИИ БИОМЕДИЦИНСКИХ ИЗОБРАЖЕНИЙ

В работе приведены математические модели адаптивной фильтрации биомедицин -.

Контурное описание; адаптивная фильтрация; гладкая функция; дифференциальный оператор; корреляционная матрица; собственный вектор.

E.P. Popechitelev, I.V. Razin

MATHEMATICAL MODELS OF THE ADAPTIVE FILTRATION OF BIOMEDICAL IMAGES

In this work are presented mathematical models of an adaptive filtration of biomedical images.

Contour description; adaptive filtration; continuously differentiable function; differentiation operator; correlation matrix, eigenvector.

. -

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

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

Задача контурного описания далеко не тривиальна - ее решению посвящены , . Наиболее широкое практическое применение нашли два метода контурного описа-,

14S

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

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

Е(j,а)-кj)]2 = J j[j)]2w[^y+G),ф^]^(y+G)^j0+1), (1)

где [u,Vф(i, j)) - ф”"^ (i, j)cosа + ф^+^г', j)sinа - производная изображения по направлению единичного вектора, u - i cos а + jsin а; ф(1+, ф(/°+1)] - совместная плотность распределения производных изображения ф^^ и ф(°+1) по

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

Emax(i, І) = 2

min 2

( 1------------------------------------------------------------------------2-----------------------------^

2 - -2

НК+[фтГ ФгГ-Nr]2} +4

фГ’фГ1

Средняя энергия приобретает экстремальные значения £шах ( ]), если про-

ш1л

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

1 (i,j) (2), (3) .

а1,2 ( 1 )= агс1ё

„(0+1)

і2

„(1+0)

і2

{[ф(1+0)_ 2 [ф(Г 1

+ 4

Ф0+О,Ф(О*“

У

2

ф(1+0)ф(0+1)

(3)

Таким образом, результатом решения уравнения (1) является поле, каждому элементу которого ставится в соответствие два ортогональных собственных вектора матрицы преобразования. Отметим, что величины векторов, определяемые выражением (2), инвариантны к повороту в силу инвариантности характеристического уравнения к любым преобразованиям системы координат [4]. Важно отметить,

что в (2), (3) используют сглаженные оценки производных Ф(;1+0) и ф-°+1^ по направлениям осей, полученные путем свертки функции изображения Ф с диффе-ренциально-сгл^ивающими операторами, приведенными соответственно на рис. 1,а,Ь. Корень квадратный из Етах (г,у) по определению представляет собой оценку максимального модуля производной изображения по направлению в среднеквадратическом смысле, а а1 (г, у) - оценку направления максимальной скорости изменения для функции яркости изображения.

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

модели перепада яркости составило 0,144°, а максимальное СКО (а^ не превыси-

ло 0,485°

Рис. 1

Ь

а

а

с

Локализация и усреднение оценок энергий производных [фу-] , [ф'у ]

в (2), (3) обеспечивается с помощью операции свертки выражений под чертой сверху с анизотропными операторами, приведенными соответственно на рис. 1,еД Анизотропные свойства проявляются в заметной на рисунках преимущественной вытянутости их весовой функции в направлении одной из осей. Эти операторы синтезированы с помощью гладких функций Лоренца, Гаусса, Моффата и Баттер-ворта и приведены в работе [5]. Специально вводимая для этих операторов анизотропия позволяет «усилить» эффект сглаживания энергий производных в направлении, ортогональном направлению дифференцирования, и обеспечивается заданием единственного параметра выбранных для сглаживания гладких функций.

Так, для сглаживания энергии производной ф(у+0) используется оператор, пред-

ставленный на рис. 1,с, а для сглаживания энергии производной ф-у - оператор

на рис. 1Д Напротив, локализация и усреднение оценки взаимной энергии произ-

оператора. Такой оператор обладает свойством центральной симметрии. Необхо-

,

, ( ) , удовлетворяет условиям определения совместной плотности распределения производных изображения в (1).

Таким образом, за счет дополнительного усреднения энергий выражения (2), (3) -

, -.

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

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

. ,

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

В основе альтернативного подхода к формированию контурного описания изображения лежит задача синтеза дифференциально-сгл^ивающего оператора, весовая функция которого была бы адаптирована сообразно полученной априори оценке локальной ориентации перепада яркости изображения. В качестве таковой предлагается использовать оценку, формируемую согласно выражению (3). При такой постановке задачи весовая функция дифференциально-сглаживающего оператора должна быть пространственно-зависимой функцией, адаптированной в каждой точке изображения в соответствии с оценкой направления максимальной скорости изменения функции яркости изображения в этой точке. В этом случае математическая модель фильтрации в аналоговой форме представима интегралом Фредгольма первого рода [1]:

Процедура адаптации означает, что для каждого элемента изображения с координатами (і, у) необходимо осуществить поворот весовой функции Ь (, п, і, У)

-

угол а1 (і, у). Используя известные соотношения, связывающие координаты ис-

водных

обеспечивается с помощью изотропного сглаживающего

(4)

-^ -^>

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

(k, n, i, j) = h(1) (kcosa1(i, j) + nsinaj(i, j))* (-ksinaj(i, j) + ncosa1(i, j)) ,2(5) h2 (k, n, i, j) = h!^2) (kcosaj (i, j) + nsinaj (i, j))h* (-ksinaj (i, j) + ncosaj (i, j)) ,(6)

где2 hs (k) = exp (-k2/2g0g ), hL (k)=-------------1- , hM (k)= 1

b ’ M V ' / л / 2 \b ’

Gom )

1+(k 2/gOl ) (l+k Vc

hB (k) =--- - функции соответственно Гаусса, Лоренца, Моффата и Бат-

1 + (к/Gob )

\b

З

терворта; hG1 (к) = (- к/o2G) exp (- к2/2o2G), hL1 (к) =

-4кЗ

g4l (i+к V g;

2

j (1) / і \ -4к , (1) /, ч -2к

hM (к) =------------------- , hB (к) =----------------------- — производные первого

Gi2m (l + к7Gi2m ) Gi2b (l + к7Gi2b )

порядка этих функций; a hG2) (k )=[(k2 -g2g j/o^Jexp (- k);

) ( к ) ; hM2) ( к ) =

З2к 12k G2l (l + k / G2L К ,( 2) (7 ) = 24k 4g2M (l + k / G2M К

З 7 ' "M

g2l (l+k V g2l ) g2m (l+k 7G

2M ,

,(2)^,. \ = 8k 2g2B (l + k /G2B )

ЬБ^ (к) =--------- -------3—- - производные второго порядка соответст-

ст2б I1+к 7 ст2б )

вующих гладких функций при Ь=2 [5]. Очевидно, что сомножители а** в зн амена-телях приведенных выражений для производных можно опустить.

Выражения (5), (6) имеют и другую форму записи:

Ь*1 (к,п,г, у) = Ь*1) (ксо8а1(г, у) + дата^г, у))Ь* (ксоШ2(1, у) + дата2(г, ])),

Ь*2 (к, п,г, у) = Ь*2) (ксо8«!(г, у) + датс^(г, у))Ь* (ксо8«2(-, у) + дата2(-, 7)) .

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

Для дискретной формы записи выражения (4) выходной сигнал адаптивного дифференциально-сгл^ивающего фильтра первого порядка с учетом (5) выражается следующим образом:

Ni, j)=

. (7)

Цф(і -к у -п) Ь) (ккка! (і, у)+даіпа1 (і, у))) (ксска2(і, у)+даіпа2(і, у))

к п

В свою очередь выходной сигнал адаптивного дифференциально-сглаживающего фильтра второго порядка с учетом (6) определяется выражением

2 Нижний индекс при Ь заменяется заглавной буквой выбранной гладкой функции О, Ь, М или Б.

152

У2ф(г, 7) = - к, ] - п)Ь*2) (ксгасс^г, Ц) + пяпос^г, 7))Ь*(ксо8а2(г, ]) + пяпо^О', 7)). (8)

к п

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

Ь*, Ь(1), и Ь*2) произведя в них формальную замену переменных.

Рис. 2

Ь

а

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

его среднеквадратической частоты [6]. Тогда выражение Етах(1,Х>/Ет1п(1,1>

определяет максимальное отношение локальных среднеквадратических частот и дает оценку степени анизотропии перепада яркости. Чем больше величина этой оценки, тем более протяженную форму имеет перепад яркости, и тем большую величину параметра «сгл^ивания» О0* мы можем задать, определив ее выражением

а о* (М) = Етах (1,.1)/Ет,п (1, ]) ха1,. Тем самым, изменяя пространственные

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

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Гонсалес Р., Вудс Р. Цифровая обработка изображений: Пер. с англ.; Под ред. ПА. Чо-чиа. - М.: Техносфера, 2005. - 1070 с.

2. Кожевников Н.И., Краснощекова Т.И., Шишкин Н.Е. Ряды и интеграл Фурье. Теория

. . . - .: , 1964. -183 с.

3. Разин КВ., Эмдин B.C. О системе инвариантов энергетического спектра градиентных изображений произвольного порядка применительно к анализу текстуры // Автометрия.

- 2003. - Т. 39, № 4. - С. 93-108.

4. . - . - .: , 1967. - 778 .

5. . .

перепадов яркости изображения // Известия ЮФУ. Технические науки. - 2009. - № 10 (99).

- С. 77-82.

6. . .

полей посредством нулей. // Известия СПбГЭТУ: Сб. тр. «Биотехнические системы в медицине и экологии». - 2006. - Вып. 2. - С. 59-69.

7. . ., . . -

рования контурного сигнала изображения // Известия вузов Радиоэлектроника. - 2007.

- № 1.- С. 24-36.

Попечителев Евгений Парфирович

Санкт-Петербургский государственный электротехнический университет “ЛЭТИ”. E-mail: [email protected].

195257, г. Санкт-Петербург, Северный пр., 65/1, кв. 169.

Тел.: +79219465462.

Разин Игорь Вениаминович

E-mail: [email protected].

197376, г. Санкт-Петербург, ул. Профессора Попова, 5.

Тел.: 88123464487.

Popechitelev Evgenij Parfirovich

Saint Petersburg Electrotechnical University “LETI“.

E-mail: [email protected].

65/1, Severnyj pr., app.169, Saint Petersburg, 195257, Russia.

Phone: +79219465462.

Razin Igor Veniaminovich

E-mail: [email protected].

5, Proffessor Popov street, Saint Petersburg, 197376, Russia.

Phone: +78123464487.

УДК 615.471:612.143

. . , . .

-

ГЕМОДИНАМИЧЕСКИХ ИССЛЕДОВАНИЙ

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

Функциональная проба с дозированной компрессией; объемнометрические изменения; спектрограмма осциллометрического сигнала; биомеханическая фильтрация.