CC BY
45
АЛГОРИТМЫ НЕЛИНЕЙНОГО СГЛАЖИВАНИЯ ЯМР-ТОМОГРАММ, ОСНОВАННЫЕ НА НЕИЗОТРОПНОЙ ДИФФУЗИИ
Institute of Measurement Science, Slovak Academy of Sciences, Dubravska cesta 9, 842 19 Bratislava, Slovak Republic
1.Введение
В ЯМР-томографии (томографии на основе ядерного магнитного резонанса), которая представляет одну из наиболее быстро развивающихся областей медицинской интроскопии, значительный интерес исследователей уделяется в последнее время быстрым методам сбора ЯМР-данных. Кроме прочих доводов, мотивирующих это направление, следующие два кажутся быть особенно важными:
1)элиминирование или минимизация влияния артефактов, возникающих в результате движения органов,
2)образование достаточного объема данных, нужных для трехмерной (3D) визуализации изображения в приемлемом для медицины промежутке времени (минуты).
Более того, быстрые методы ЯМР-томографии (напр."3D Flash Sequence") могут увеличить значение томографического отображения на клинической практике, так как они обеспечивают увеличение числа пациентов обследованных на томографической установке, а также доставляют пациенту больший комфорт.
Паралельно с прогрессом достигнутым в ЯМР-томографии развиваются и методы анализа и описания ЯМР изображений, которые улучшают процесс медицинской диагностики и совершенствуют обучение в рамках университетских курсов по невроанатомии и неврологии ([1 ] / [2 ]) • Так как наиболее увлекательным видом представления визуальных познаний, особенно в медицине, является (3D) объемная визуализация, появление быстрых ЯМР методов послужило как бы новым вызовом в этой области. Сегментация и улучшение изображений представляют два основных метода обработки медицинских изображений на нижнем уровне ("low-level processing"), которые необходимы для решения более сложных задач интерпретации этих изображений. Однако, если такие задачи относятся к исходным данным, полученным быстрыми ЯМР методами, то перед решением этих задач приходится заниматься проблемой низкого отношения сигнал/шум (SNR).
В 1990 году Перона и Малик ([3]) предложили новый метод нелинейной фильтрации, который в значительной степени уменьшает нежелательный эффект быстроого сбора ЯМР-данных. Их подход к "off-line" обработке основывается на физическом понятии Чеиэотропной диффузии. Первоначальное применение этого метода сглаживания изображений к построению масштабного пространства ("scale-space") и к детекции границ областей было расширено Геригом (Г 41 ) на сглаживание "dual-echo" ЯМР изображений мозга. метод неизотропной лиффузии, который рассматривается
в статьи Герига, относится исключительно к двухмерным случаям. Для того, чтобы коэффициент диффузии - функция пространственных координат и времени - отражал отношение шума и резкости границ в данном изображении, при его определении вводится некоторый параметр К. Этот параметр [4] вычисляется на основе кумулятивных гистограмм градиентов (их модуля) исходного изображения и не зависит от итеративного процесса неизотропной диффузии. В статье исследованы также возможности использования многомерной природы ЯМР сигналов в задаче сглаживания и последующей сегментации. Ряду вопросов, однако, внимание уделено небыло. Перечислим те из них, которые будут рассмотрены в данной статье:
-зависимость значений параметра К от этапов итерации, -альтернативная схема дискретизации уравнения неизотропной диффузии и цифровых градиентов,
-количественный критерий для остановления процесса итераций ,
-характеристика конечного эффекта сглаживания изображения.
2.Сглаживание с помощью неизотропной диффузии
В общем случае диффузия рассматривается как изотропный процесс, который характеризуется константой для данного объема вещества - коэффициентом диффузии с . Для двухмерного скалярного поля 1(х,у) (в данном случае поля яркостей изображения) процесс диффузии описывается следующим дифференциальным уравнением в частных производных:
- сНу[с-У/(х,у,«)] (1)
где х,у-пространственные координаты, <:-время или параметр эволюции процесса; для дискретного случая t обозначает номер итерации. Сглаживание изображений, основанное на мультимас-штабном ("пш1г:1зса1е") подходе, которое предложено в статье
[3], требует частичного или полного подавления диффузии
в изображении на границах областей. Это требование не что
иное, как введение подходящей функции проводимости с(х,у,Ь) вместо константы коэффициента диффузии. Таким образом уравнение (1) принимает вид
= <Иу[с(д,у,0’ У/(г,у,*)] (2)
Если предположить, что в данном моменте времени (на данном масштабном уровне) t, положение границ, соответвуюших этому моменту, известно, то цель сглаживания при помощи неизотропной диффузии можно сформулировать как предпочтительное усиление сглаживания изображения внутри однородных областей и одновременно предпочтительное подавление сглаживания через граничные пикселы. Эту цель можно было бы идеально достичь, если бы функция с(х,у,Ь) принимала только два значения: значение "1" в пикселах внутри областей, и значение "О" в пикселах на
границах областей.
Естественно, что положение границ заранее неизвестно и функцию с(х,у,Ь) можно строить только как некоторое приближение к идеальной функции со значениями из отрезка <0,1>. Пе-рона и Малик предложили использовать в качестве индикатора границ следующую показательную функцию градиентов яркостей изображения:
с(У/(х,у,*)) = ехр < -
Р7/(х,у,0|
(3)
Параметр К , который входит в экспонент, как мы уже упомянули в введении, играет роль некоторого фактора взвешивающего влияние шума на определение границ областей посредством оператора градиента яркостей изображения. На рис.1 иллюстрируется зависимость формы графика функции проводимости с(х,у,Ь) от значений параметра К , причем конкретное значение можно интерпретировать как введение некоторого порога градиентов для классификации граничных пикселов. Для уточнения этого соображения рассмотрим некоторый фиксированный уровень диффузии, (напр, с = 0.1), которому будет соответствовать разное сглаживание. Тогда для отдельных значений параметра К , которые определенным образом отражают уровень шума в изображении, получаем различные значения градиентов , Ух2 , . . . , которые можно считать порогами для пикселов, принадлежащих границам областей. Чем больше пороговый градиент, тем большие значения градиентов допускаются для сглаживания диффузией.
При выборе метода вычисления параметра К авторы [3], [5] исходят из следующих соображений. Предполагается, что распределения амплитуд градиентов для шумовых точек и точек границ отличаются. Вычисляется кумулятивная гистограмма градиентов во всем изображении (глобальная гистограмма) и принимается предположение о том, что нижних 80 % представляет только энергию шума. Перона и Малик используют ту же "оценку шума", однако в качестве параметра К берут 90 % квантиль глобальной гистограммы градиентов.
Как показали наши предварительные компьютерные эксперименты, первоначальное вычисление параметра К , которое используется как константа во всем последующем процессе диффузии, недостаточно, поскольку отрицательно сказывается на качестве сглаживаемых изображений, получаемых на последующих итерациях. На основе этих результатов мы модифицировали процедуру неизотропной диффузии таким образом, чтобы параметр К вычислялся для каждой итерации заново.
На процесс неизотропной диффузии будет оказывать влияние кроме параметра К также выбор формулы для цифрового градиента VI. Перона и Малик ([3]) утверждают, что этот выбор дает результаты, которые визуально неотличимы. Чтобы убедиться в этом, в следующем параграфе рассмотрим четыре вида определения таких градиентов и проведем тестовые вычисления, необходимые для обработки реальных ЯМР-томограмм.
Рис. 1
3.Некоторые вопросы дискретного вычисления 20 неизотропной диффузии
Опишем коротко стандартный переход от непрерывной формулировки задачи диффузии к ее дискретному варианту. Для дискре-тизации в [3], [4] используется уравнение (2), переписанное в виде
Сначала аппроксимируется частная производная яркости в скобках: для каждой пространственной координаты получается конечная разность, затем проводится конечная аппроксимация внешней производной. В результате получаем формулу:
~~оі'П' ~ дЬ{с(х+ +Д*.У»0“Л*.У.01 -
А т 'І
Ф - • №(Х>У> 0 - 1(х - Д*»У.О] ) +
с(х,у- • [/(х,у,0 - /(х,у - Ду,0] } (5)
Далее производная яркости относительно времени аппроксимируется разностью, соответствующей двум последующим итерациям, причем специально вычислена интегральная константа времени, которую будем обозначать Л ([4]). В результате получаем дискретное уравнение для пиксела (1,3), которое после некоторого преобразования принимает вид:
1\Г = + с5П51 + съ-Ву,! + сЕ-ВЕ1] (6)
обозначают конечные разности яркостей для соответствующих соседних пикселов и центрального пиксела (1,3) :для (1-1,3), ? для (1+1,з), х=Б ; для (1,3-1), х=М ; для
(1,3+1), х=Е . Следует отметить, что значения коэффициентов
диффузии с^, с3, см, с£ как видно из уравнения (5) вычисляются для дуг (на расстоянии 1/2 дх, 1/2 ду между пикселом
(1,3) и соседними пикселами Ы, Б, И, Е .
Для того, чтобы максимально упростить итеративное решение этого дискретного уравнения неизотропной диффузии, в работах
[3] и [4] значения амплитуд градиентов при вычислении функции проводимости (определенной на растре изображения) аппроксимируются несимметричными разностями яркостей в соседних пикселах и в центральном пикселе. Таким образом параметр шума К вычисляется на основе гистограммы разностей, а не на основе градиентов.
С одной стороны можно сказать, что с появлением быстрых процессоров класса "Н13С" вычислительные доводы для всех упрощений, принятых в упомянутых работах, теряют свой весь. С другой стороны требования более точной и надежной сегментации ЯМР изображений, специально мозга, которая в значительной
степени зависит от достаточного отношения сигнал/шум, приводят в данном случае к необходимости устранения упрощающих элементов алгоритмов 20 неизотропной диффузии. Поэтому мы считали целесообразным перейти к более точным вариантам этого алгоритма, провести соответствующие вычислительные эксперименты, и наконец сравнить полученные результаты.
С целью однозначного описания отдельных вариантов алгоритма 20 неизотропной диффузии нам понадобится схема растра изображения, которая приводится на рис.2. Как мы уже упомянули, дальше для упрощения записи индексов отдельных величин определяемых в 4-окресности центрального пиксела (1,3), будем использовать символы: И, Б, У, Е (рис. 2). Расстояние между
пикселами обозначим через дх (вдоль оси х) и ду (вдоль оси у). В данном случае полагаем = 1. Стандартные не-
симметричные разности обозначим следующим образом:
ю -^ - 07 - 4} - *1-1.].
Рис. 2
3.1 Цифровые градиенты для определения функции проводимости
с(х,у,с;
.А..
Первым возможным уточнением вычисления функции проводимости с(х,у,Ь) является переход от разностей ОН и ОУ к симметричным градиентам СС1 или СС2 , которые определены (относительно центрального пиксела) известными формулами:
С&(»,з) = 1/2 • (|/£ - ы\ + |/5 - /д-1) (7)
или
СС2(г\;) = 1/2 • \/(/£ - Ы)2 + (/5 - /л-)2 (8)
Поскольку для вычисления дискретного уравнения (6) необходимы значения функции с() в точках между пикселами (на расстояниях Дх/2 и Ду/2) (на рис.2 они отмечены крестиками), соответствующие значения симметричных градиентов СС-^ или СС2 следует интерполировать.
Б - Отметим, что несмотря на большую чувствительность градиентов относительно шума, на данном этапе обработки, т.е. для вычисления параметра шума К на основе гистограмм градиентов, нас интересует более детальная характеристика изображения. Поэтому, вторая возможность уточнения дискретного алгоритма неизотропной диффузии кажется также подходящей. Она состоит в переходе к симметричным градиентам Робертса (СБ2), которые даны (напр. [6]) следующей формулой:
65,(1,Я = 11 у/2- (|/5 - 1е\ + 11г+и+1 - 1ц|) (9)
В последующем изложении мы не будем приводит подобные формулы для градиентов, определяемых при помощи квадратного корня. На рис.2 используются следующие обозначения Д+ = 13 - 1Е Д" = 1^+1 -1 + 1 - 1д_-; . Таким образом, в этом случае получаем точные значения градиентов Робертса в субпикселах, обозначенных относительно центрального пиксела ИУ?, БЕ, SW, ИЕ. Опять-таки, нужные значения градиентов для точек, расположенных на растояниях Дх/2 и Ду/2 (на рис.2 - крестики) можно вычислить при помощи интерполяции (напр.линейной).
3.2 Альтернативная схема дискретизации уравнения неизотропной диффузии
Вместо того, чтобы последовательно аппроксимировать уравнение диффузии, записанное в виде (4), можно сначала раскрыть скобки и преобразовать уравнение следующим образом:
Теперь можно перейти к аппроксимациям производных конечными разностями, причем на этот раз получаем также разности функции проводимости с(х,у, Ь); &хс() ~ для горизонтального и 0^с() - вертикального направления. Для неизотропной диффузии получаем следующее альтернативное дискретное уравнение (Ь - номер итерации):
/'«(*, ») = /'(*,») + А lD.eQ-D.IQ + 0,с() С„/() + с() ■ Д/()] (11)
где Л -интегральная константа того же вида как прежде. Через ДІ обозначаем цифровой Лапласиан определенный для центрального пиксела (і,і).
В .
Естественно, что подобный выбор вида цифровых градиентов, который мы рассмотрели в частях А, Б, можно сделать и в случае второй дискретной схемы. Если взять, например, центральные градиенты, то для данного растра изображения получаем полный набор точных значений с(х,у,Ь) , т.е. значений коэффициентов диффузии для каждого пиксела (по сравнению со стандартной схемой Перона и Малика не требуется интерполяция значений в промежуточных точках). Используя обозначение из рис.2, уравнение (11) запишем в следующем виде:
/‘+1(»,Л = /*(*,Я+а{[с‘(*>я - с\у] • юн + [с‘(»,я - 4,] • /ж + с*(1,я • д/‘(г,;)}
(12)
37 -
Наконец опишем интересную схему, которая сочетает альтернативную дискретизацию уравнения неизотропной диффузии с вычислением всех нужных конечных величин на субпикселах. Введем сначала (суб)растер (полные точки на рис.2), который повернут относительно исходного растра на Я/4 и шаг которого = ЁУ' new !) относительно шага исходного растра (old 1)
= 1^( 2 Y~2) . Следует отметить, что так как имеет место
Д* = Д / ДУ = Д ' (Рис-2)» значения градиентов Робертса в исходном растре становятся в новом растре значениями центральных градиентов. При помощи линейной интерполяции, из значений яркостей в пикселах исходного растра получаем значения в соответствующих субпикселах: INW , igE , lsw , iNE. Таким образом в соответствующем уравнении, полученном из (12), используем:
а) точно вычисленные конечные разности функции проводимости (CNW , CgE , CgW , CNE), __ __
б) разности яркостей в пикселах и субпикселах (DH, DV
- рис.2) .
Значения величин определенных в субпикселах используются при этом исключительно в процессе промежуточных вычислений, но диффузия пробегает только в пикселах исходного растра.
4.Определение функции для критерия остановления итеративного процесса дискретной неизотропной диффузии
Для сглаживания ЯМР-томограмм при помощи неизотропной диффузии служит итеративный алгоритм. В работах [3], [4] процесс итераций основан на эмпирических выводах, т.е. оптимальное число итераций устанавливается при помощи качественной оценки результатов. Вопрос подходящей численной характеристики окончательного эффекта сглаживания с сохранением границ ("Edge-Preserving Smoothing") также не расматривался.
Как мы уже упомянули в введении, одной из актуальных задач медицинской интроскопии является задача 3D визуализации мозга на основе 2D томограмм, полученных быстрыми ЯМР методами. Окончательной целью такой визуализации является, конечно, разработка программ, которые позволяли бы врачу, без дополнительных знаний и усилий получать 3D изображения мозга (или других органов) на дисплее. Это приводит к требованию - на всех этапах предобработки ЯМР-данных, к которым относится и сглаживание, использовать алгоритмы, управляемые минимумом параметров задаваемых пользователем.
Выше перечисленные моменты послужили доводом к построению подходящей функции - численного критерия для автоматической остановки процесса итераций ([7]). Количественные меры, которые используются для исследований в области компьютерного моделирования сбора и обработки томографических изображений, содержат важный методический элемент, а именно -информацию об идеальных изображениях, т.н. фантомах. Опять-таки, поскольку медицинская практика нуждается в количественной оценке применения алгоритмов к реальным томографическим данным, численные характеристики, используемые в вычислениях, которые должны быть транспарантны для пользователя, нельзя основывать на фантомах. Таким образом нашей целью в этом параграфе является выведение разумной функции в качестве останавливающего критерия, а также построение количественной меры эффекта сглаживания.
Шринивасан ([8]) предложил численную меру, отношение сигнал/шум (SNR), для оценки выполнения алгоритмов реконструкции изображения по проекциям. Она характеризует глобальную точность реконструированного изображения относительно идеального изображения (фантома). SNR определяется следующей формулой:
SNRldB) = 10. log „Щ&кЗ- (13)
2^1=1 z^j=i \l\j Jtj\
где ■ - обозначены уровни яркостей входного фантома,
- уровни яркостей в выходном (реконструированном) изображении. Размер изображений: N х N пикселов . В качестве фантома
(tji} обычно используется фантом Шеппа и Логана ([9]). SNR
- функция явно ведет себя подобно стандартной характеристике отношения сигнала и шума, используемой в теории сигналов. Такое поведение совпадает с основной мотивировкой задачи згла-живания ЯМР изображений мозга при помощи неизотропной диффузии, именно с целью этой операции - повысить значения отношения сигнал/шум. Это обстоятельство оправдывает наш выбор
- взять как раз формулу (13) за основу искомой функции - кри-
терия для остановления итераций . Так как при характеризации результатов сглаживания томограммы после (р+1) итерации (р-1), данными (tji) (идеальным фантомом) мы не располагаем, предлагаем модифицированную функцию, т.н.относительное отношение сигнал/шум ('^snr^ ‘ Определим эту функцию следующей
формулой:
(Р+1 г. л2
р+'115Ын{с1В) = 10.log ,——------------------------* (14)
Е& (р/о -р+1 М
В которой / Р+1^1п уровни яркостей изображений, получен-
ных после р-той и (р+1)-ой итерации. Критерий остановления итераций потом определяется как
»'8С = 11»^ -’’•** < ,кг (15)
Итак, если задаться частным значением порога 1:11г (обычно = 0.05), то процесс итераций будет остановлен при выполнении условия (15). Вместо функции (14) можно использовать ее нормированную версию. Пусть р+1л?ах и определены
следующим образом:
p+lmax = maxij (р+1/о)> v+xmaxdif = max.j \pfij -p+1 /tJ|.
Тогда формула для нормированного относительного коэффициента определяется как
'+‘ЯЛГ5т|(1Ш) = 10.108 (16)
и’=1 \ * таг<Н} )
Если обозначим через { иГ^ ^) уровни яркости входного изображения и через { -1^-;) уровни яркости окончательного изображения, полученного при оптимальном числе итераций сглаживания, и определим величины пах и тахсШ как минимальные значения для всей диффузии, то достигнутый эффект сглаживания, которое сохраняет границы областей, можно характеризовать функцией
5. Результаты вычислительных экспериментов 2D сглаживания
Для сравнения численных характеристик, а также для визуальной оценки результатов применения отдельных алгоритмов, отличающихся выбором цифрового градиента и схемой дискретизации уравнения неизотропной диффузии, в задаче сглаживания ЯМР томограмм головы (мозга), мы провели набор вычислительных экспериментов. В качестве тестовых изображений были избраны три томограммы (сечения) из 3D куба (2563 вокселов), полученные при помощи ЯМР-томографической установки "Magnetom SP (Siemens - 1.5 Т)" в радиодиагностической клинике больницы имени "Derera" в городе Братислава (Словакия). Была использована т.н. "3D Flash" измерительная последовательность со следующими параметрами:
TR = 30 ms, ТЕ = 5 ms, F0V = 256 mm, flip angle = 40
размер изотропного воксела lxlxl mm.
Эксперименты проводились на рабочей станции НР-Аро11о 720. Для автоматической остановки процесса итераций мы использовали критерий - функцию RNsnr (16)
Так как результаты для отдельных срезов лишь минимально отличались друг от друга, мы приводим данные, относящиеся только к одному выбранному случаю. Во всех случаях вычислительных экспериментов, результаты которых описываем, использовались градиенты, определенные через абсолютные значения разностей, так как использование второй версии градиентов, которая требует больше вычислений, не показало никакой разницы.
На рис.З показаны исходное и обработанное изображения (сечения головы на уровне глаз). Рис.За показывает зашумленный оригинал, на рис.36 изображение, сглаженное неизотропной диффузией при помощи стандартного алгоритма Перону и Малика (дискретизация по первоначальной схеме (4), исходный растр, Для вычисления функции проводимости c(x,y,t) и шумового параметра К использованы вместо градиентов разности).
На рис.Зв показано изображение, сглаженное при помощи той же самой схемы диффузии и центральных градиентов для получения значений функции проводимости c(x,y,t). Наконец, рис.Зг показывает результат неизотропного диффузионного сглаживания на основе альтернативной схемы дискретизации уравнения (11) в повернутом растре, включительно субпикселов и с использованием центральных градиентов в этом субрастре (или градиентов Робертса в исходном растре) для вычисления Функции c(x,y,t).
Рис • 3
Для большей наглядности на рис.4 изображены увеличенные одинаковые фрагменты тех же самых изображений как на предыдущем рисунке.
Рис. 4
Оказывается, что все рассмотренные варианты цифровых градиентов и схемы итеративного процесса неизотропной диффузии дают лучшие результаты, чем исходный алгоритм. Отметим, что в изображении полученном при помощи этого алгоритма, особенно очевидны негладкие границы областей, которые уменьшают качество последующей ЗБ визуализации.
Хотя взаимное визуальное сравнение результатов рассмотренных вычилительных альтернатив показывает небольшие отличия, нам кажется, что альтернативная схема и использование субпикселов приводят к большему сглаживанию при удовлетворительном сохранении границ. Эти выводы потверждаются количественными результатами экспериментов, которые приведены в таблицах 1,
2, 3. В этих таблицах в первом столбце представлены значения параметра К (порогового градиента, определенного как 90 %
квантиль гистограммы градиентов). Второй столбец представляет номер итерации (последний номер является одновременно оптимальным числом итераций). Наконец, в последнем столбце приводятся значения относительного отношения сигнал/шум, полученного для каждой итерации сглаживания.
Таб.1 относится к алгоритму Перону и Малика (изображение на рис.36), таб.II к случаю с центральными градиентами (соответствует рис.Зв ), таб.III- соответствует случаю альтернативной схемы и субрастра (рис.Зг).
Таб.1
Таб.II
Таб.III
к t
21.0 1 26.84
17.0 2 33.04
15.0 3 37.07
13.0 4 39.80
К 1 ЛДГ5Л^Я
33.0 1 24.19
27.0 2 31.99
25.0 3 35.92.
23.0 4 38.70
К 1 К-МSN Я
23.0 1 26.59
19.0 2 31.55
17.0 3 34.70
16.0 4 36.67
Для этих трех алгоритмов получены следующие значения функции Т2МК, характеризующей улучшение отношения сигнал/шум
N Т^Л’Я
I. 21.96
II. 20.42
III. 19.32
Отметим, что поскольку предложенная формула на этот раз включает результаты двух последующих итераций сглаживания, так чем ниже значение функции тем получаем большее
улучшение отношения сигнал/шум. Поскольку приведенные значения относятся к результатам сглаживания , которые получены как оптимум - в смысле критерия (15)- имеем единую основу для
всех сравниваемых алгоритмов. Итак, можно сделать вывод, что наибольший эффект сглаживания неизотропной диффузией (при сохранении границ) дает третьий (III) вычислительный алгоритм.
6.ВЫВОДЫ
На основе проведенных компьютерных экспериментов можно сделать следующие выводы:
1) Выбор формулы цифрового градиента определяет значение квантиля гистограммы градиентов, т.е. значение параметра К показательной функции с() проводимости . Следовательно, конкретная форма этой функции влияет на выполнение неизотропной диффузии, т.е. на эффект сглаживания ЯМР-томограмм.
2) Альтернативная схема дискретизации дифференциального уравнения неизотропной диффузии, использующая субрастер изображения, обеспечивает более плотную аппроксимацию непрерывной функции проводимости конечным числом точек, что обеспечивает более гладкое сглаживание при сохранении границ. Более того, дискретное уравнение неизотропной диффузии в этом случае содержит разности яркостей изображения и разности функции проводимости с(), причем обе эти величины полностью определены в точках субрастра и растра. Таким образом, исключение аппроксимации функции проводимости значениями на дугах между соседними пикселами позволяет единым способом рассматривать дискретное уравнение неизотропной диффузии для точек любой окрестности центрального пиксела, в том числе и ранговой, что будет одной из целей последующих исследований.
3) С одной стороны переход от несимметричных разностей, использованных в стандартной схеме, к симметричным градиентам упрощает вычисления, так как не нужны горизонтальные и вертикальные разности, и таким образом компенсируются дополнительные затраты на вычисление самых градиентов и на интерполяцию нужных значений. С другой стороны этот переход делает обработку изображений независимой от их характера (их симметрии).
4) Вычислительные эксперименты подтвердили, что разработанный нами критерий для автоматического остановления процесса итераций работает точно и надежно. Оптимальное число итераций для данного типа ЯМР-томограмм было 4-6, что меньше, чем число 10, которое приводится как оптимум в работах [3],
[4].
Для проведения вычислительных экспериментов с предложенными версиями алгоритмов сглаживания ЯМР-томограмм мозга при помощи неизотропной диффузии были разработаны программные средства, которые будут использоваться в последующих компьютерных экспериментах. В настоящее время продолжается разра-
ботка 30 алгоритмов сглаживания и исследование оптимальных цифровых схем для этого более общего случая. Следует отметить, что вычислительные затраты отдельных алгоритмов (стандартных и альтернативных) отличались друг от друга лишь немного: обработка одного изображения размером 256 х 256 пикселов длилась несколько секунд (<10) . Вопросы вычислительных
требований в случае 30 сглаживания нуждаются, конечно, в более подробном анализе.
Автор статьи выражает свою благодарность коллегам И.Холлендеру и А.Писару за помощь, оказанную при проведении компьютерных экспериментов.
Работа поддержана грантом N.2/1016/94 Словацкого грантового агентства для науки.
7.ЛИТЕРАТУРА
[1] Bomans М. et al. 3-D segmentation of MR images of the head for 3-D display // IEEE Trans, on Medical Imaging, Vol. MI 9, No. 2, pp. 177-183, 1990.
[2] Stanier J. et al. Segmentation schemes for knowledge-based construction of individual atlases from slice-type medical images // Proc. of SPIE, Vol. 1898: Image Processing, pp. 252-262, 1993.
[3] Perona P. , Malik J. Scale-space and edge detection using anisotropic diffusion // IEEE Trans, on Pattern Analysis and Machine Intelligence, Vol. PAMI 12, No. 7, pp. 629-639, 1990.
[4] Gerig G. et al. Nonlinear anisotropic filtering of MRI Data // IEEE Trans, on Medical Imaging, Vol. MI 11, No. 2, pp. 221-232, 1992.
[5] Canny J. A computational approach to edge detection // IEEE Trans, on Pattern Analysis and Machine Intelligence, Vol. PAMI 8, pp. 679-698, 1986.
[6] Gonzalez R. C. , Woods R. E. Digital image processing // Addison-Wesley , Reading, 1992.
[7] Bajla I. et al. Anisotropic filtering of MRI data based upon image gradient histogram. Ed. by Chetverikov D. , Kropatsch W. G. // Proc. of CAIP 93-Computer Analysis of Images and Patterns, Springer-Verlag, Berlin, pp. 90-97, 1993.
[8] Shrinivasan V. et ad. Image reconstruction by a Hopfield neural network // Image and Vision Computing, Vol. 11, No. 5, pp. 278-281, 1993.
[9] Shepp L. A. , Logan B. F. The Fourier reconstruction of a head section // IEEE Trans, on Nucl. Sci. , Vol. NS-21, No. 3, pp. 21-24, 1974.
