Гонсалес Р., Вудс Р. Цифровая обработка изображений (3-е изд., 2012) (1246138), страница 94
Текст из файла (страница 94)
Н. Тихонова. Такой подход называется тихоновским методом регуляризации, и рассматриваемый ниже метод является его весьма частным случаем. В связи со сказанным мы далее будем называть этот метод фильтрации,основанный на минимизации сглаживающего функционала со связью, фильтрациейпо Тихонову или тихоновской регуляризацией. — Прим. перев.5.9. Фильтрация методом минимизации сглаживающего функционала со связью421⎡ 0 −1 0⎤(5.9-5)p( x, y ) = ⎢−1 4 −1⎥ ,⎥⎢⎢⎣ 0 −1 0⎥⎦т. е.
той функции, с помощью которой в разделе 3.6.2 был определен операторЛапласа18. Как уже отмечалось выше, важно помнить о том, что функция p(x,y),равно как и все остальные функции в пространственной области, должны бытьправильно дополнены нулями перед вычислением их Фурье-преобразованиядля использования в (5.9-4). Вопрос о дополнении нулями подробно рассматривался в разделе 4.6.6. Обратим внимание, что при обращении параметра регуляризации γ в нуль выражение (5.9-4) сводится к инверсной фильтрации.Пример 5.14.
Сравнение винеровской фильтрации и фильтрации по Тихонову.■ На рис. 5.30 представлены результаты восстановления изображенийна рис. 5.29(а), (г) и (ж) с помощью фильтрации по Тихонову. Значения параметра γ подбирались вручную так, чтобы обеспечить наилучшее визуальное качество восстановления. Тот же метод выбора параметров был использован дляполучения результатов фильтрации по Винеру, которые были представленына рис. 5.29(в), (е), (и). Сравнивая между собой фильтрации по Тихонову и по Винеру, можно отметить, что первая дала несколько лучшие результаты в случаях большого и среднего по величине шума, при том что обе дали практическиодинаковые результаты в случае низкого шума.
То, что фильтрация по Тихоновуоказалась лучше фильтрации по Винеру в случае подбора значений параметроввручную, не является неожиданностью. Дело в том, что параметр регуляризации γв (5.9-4) является числом, в то время как параметр K в (5.8-6) является по сути значением постоянной функции, используемой для приближения отношения двухнеизвестных функций в частотной области, которое редко является постоянным.Поэтому естественно ожидать, что результат, основанный на подборе параметра γ«вручную», даст более точную оценку неискаженного изображения.■а б вРис.
5.30. Результаты фильтрации по Тихонову. Сравните (а), (б) и (в) с результатами винеровской фильтрации на рис. 5.29(в), (е) и (и) соответственно18Можно, в принципе, использовать и функцию H(u,v) = –(u2 + v2), которая в соответствии с разделом 4.9.4 определяет «частотный» оператор Лапласа. — Прим. перев.422Глава 5. Восстановление и реконструкция изображенийКак показывает приведенный пример, значения параметра регуляризации γможно перебирать в интерактивном режиме до тех пор, пока приемлемый результат не будет получен.
Однако если мы хотим получить оптимальное решение, то значение параметра должно быть выбрано так, чтобы выполнялось условие связи (5.9-3). Ниже приведена итеративная процедура такого выбора.Определим вектор невязки r следующим образом:ˆ(5.9-6)r = g − Hf .ˆПоскольку решение Fˆ(u, v) (и соответствующий вектор f ), определяемоепо (5.9-4), есть функция параметра γ, вектор невязки r также зависит от этогопараметра. Можно показать (см. [Hunt, 1973]), что функционал невязкиφ( γ ) = rT r = r2(5.9-7)является монотонной функцией параметра γ. Мы хотели бы так выбрать параметр γ, чтобы выполнялось условиеη −B ≤ S ≤ η +B ,(5.9-8)где коэффициент a задает приемлемую точность выполнения условия связи.В случае a = 0 имеет место ||r||2 = ||η||2, и ввиду (5.9-6) условие (5.9-3) выполняетсяточно.Поскольку функционал невязки φ(γ) является монотонной функцией параметра регуляризации, нахождение искомого значения γ не представляет трудности.
Один из алгоритмов состоит в следующем.1. Задать начальное значение γ.2. Вычислить ||r||2.3. Если условие (5.9-8) выполняется, то цель достигнута. В противном случае увеличить значение γ, если ||r||2 < ||η||2 – a, или уменьшить его, если||r||2 > ||η||2 + a, и повторить шаг 2, используя новую оценку Fˆ(u, v), пересчитанную по (5.9-4) с новым значением γ.Для увеличения скорости сходимости можно использовать более продвинутыеметоды, например метод касательных Ньютона.Для того чтобы использовать данный алгоритм, необходимы величины ||r||2и ||η||2. Для вычисления ||r||2 заметим, что из (5.9-6) следует равенствоR (u,v ) = G (u,v ) − H (u,v )F (u,v ) ,(5.9-9)и функция r(x,y) может быть получена вычислением обратного Фурьепреобразования от функции R(u,v).
ДалееM −1 N −1r = ∑ ∑ r 2 ( x, y ).2(5.9-10)x =0 y =0Вычисление ||η||2 приводит к интересному результату. Прежде всего рассмотримдисперсию шума полного изображения, которую мы оцениваем методом выборочного среднего, рассмотренным в разделе 3.3.4:σ 2η =1MNM −1 N −1∑ ∑ (η( x, y ) − m )x =0 y =02η,(5.9-11)5.9. Фильтрация методом минимизации сглаживающего функционала со связью423гдеmη =1MNM −1 N −1∑ ∑ η( x, y )(5.9-12)x =0 y =0есть выборочное среднее значение шума. Сравнивая два последних выраженияс выражением для ||η||2, которое имеет тот же вид, что и (5.9-10), мы видим, чтоη = ./ (σ 2η + Nη2 ) .(5.9-13)Это очень важный результат, который показывает, что для реализации оптимального алгоритма восстановления нам достаточно знать лишь среднее значение и дисперсию шума.
Эти величины нетрудно оценить (см. раздел 5.2.4)в предположении, что значения шума и изображения не коррелированы.Это предположение является основным для всех рассматриваемых в настоящейглаве методов.Пример 5.15. Итеративная оптимальная фильтрация по Тихонову.■ На рис. 5.31(а) представлен результат восстановления изображенияна рис. 5.25(б), полученный с помощью оптимального алгоритма фильтрациина основе только что рассмотренной итеративной процедуры. Начальное значение γ было равно 10 –5, величина изменения значений γ в итеративной процедуребыла равна 10 –6, а значение коэффициента a равнялось 0,25.
Параметры шумав алгоритме восстановления выбирались такими же, как и при генерации изображения на рис. 5.25(б) (шум имел нулевое среднее и дисперсию 10 –5). Результатвосстановления столь же хорош, что и результат на рис. 5.28(в), для получениякоторого была использована винеровская фильтрация, а параметр K подбирался вручную так, чтобы обеспечить наилучший в визуальном смысле результат.Изображение на рис. 5.31(б) демонстрирует, что может случиться при использовании неправильных оценок для параметров шума. В данном случае считалось,а бРис. 5.31.(а) Результат восстановления изображения на рис.
5.25(б) с помощьюитеративной оптимальной фильтрации по Тихонову, использующейправильные значения параметров шума. (б) Аналогичный результат,полученный при неправильных значениях параметров шума424Глава 5. Восстановление и реконструкция изображенийчто дисперсия шума равна 10 –2, а среднее значение по-прежнему равно нулю.Результат восстановления расфокусирован заметно сильнее.■Как было сказано в начале этого раздела, важно понимать, что оптимальноев смысле минимизации сглаживающего функционала со связью восстановление не обязательно является «лучшим» в смысле визуального качества. В зависимости от природы и величины искажений и шума разные факторы влияютна вид окончательного результата, получаемого с помощью итеративного алгоритма построения оптимальной оценки. Вообще фильтры для восстановления,которые строятся в автоматическом режиме, дают результаты хуже, чем фильтры, параметры которых настраиваются вручную под конкретное изображение.Это относится, в частности, к рассмотренному фильтру, который полностью задается одним числовым параметром.5.10.
Ñðåäíåãåîìåòðè÷åñêèé ôèëüòðРассмотренный в разделе 5.8 винеровский фильтр можно несколько обобщить.Это обобщение приводит к так называемому среднегеометрическому фильтру:1−α⎛⎞α ⎜⎟∗∗ˆ⎛ H (u,v ) ⎞ ⎜H (u,v )⎟F (u,v ) = ⎜2⎟⎛ S η (u,v ) ⎞ ⎟⎝ H (u,v ) ⎠ ⎜2⎜ H (u,v ) + β ⎜⎟⎟⎝ S f (u,v ) ⎠ ⎠⎝G (u,v ) ,(5.10-1)где α и β — положительные вещественные константы. Среднегеометрическийфильтр представляет собой произведение двух выражений в скобках, которыесоответственно возводятся в степени α и 1 – α.При α = 1 этот фильтр сводится к инверсному фильтру. При α = 0 фильтрпревращается в так называемый параметрический винеровский фильтр, которыйсводится к обычному винеровскому фильтру при β = 1.
При α = 1/2 фильтр представляет собой произведение двух величин, возведенных в одинаковую степень1/2, т. е. является средним геометрическим этих величин, откуда и происходитназвание фильтра. При α = 1/2 и β = 1 этот фильтр известен как фильтр эквализации спектра. При β = 1 по мере увеличения α больше 1/2 работа фильтравсе больше напоминает работу инверсного фильтра. Аналогично при убыванииα меньше 1/2 фильтр ведет себя как винеровский фильтр.
Выражение (5.10-1)весьма полезно, поскольку с его помощью в рамках одного выражения в действительности удается описать целое семейство фильтров.5.11. Ðåêîíñòðóêöèÿ èçîáðàæåíèÿ ïî ïðîåêöèÿìВ предыдущих разделах настоящей главы изучались методы восстановленияизображения по его искаженному варианту. В данном разделе будет рассмотрена проблема реконструкции изображения по множеству его проекций; основ-5.11. Реконструкция изображения по проекциям425ное внимание будет уделено рентгеновской компьютерной томографии (КТ).Это самый первый и до сих пор наиболее широко используемый тип КТ, являющийся одним из основных приложений цифровой обработки изображенийв медицине.Как отмечалось в главе 1, наравне с термином «компьютерная томография» (КТ)используется термин «компьютерная аксиальная томография» (КАТ).5.11.1.
ВведениеЗадача реконструкции в принципе проста и может быть качественно объясненанепосредственным, интуитивно понятным способом. Для начала рассмотримрис. 5.32(а), состоящий из единственного объекта на однородном фоне. Чтобывнести в последующее объяснение физический смысл, предположим, что этоизображение представляет собой сечение некой трехмерной области тела человека. Также допустим, что фон на изображении соответствует мягким тканям,а круглый объект — опухоли, тоже однородной, но с более высоким коэффициентом поглощения.Далее предположим, что мы посылаем тонкий плоский пучок рентгеновскихлучей слева направо (вдоль плоскости изображения), как показано на рис. 5.32(а),и будем считать, что объект поглощает больше энергии лучей, чем фон, как этообычно и бывает.
Поставив с противоположной стороны области линейку детекторов рентгеновских лучей, получим показанный сигнал (профиль поглощения),амплитуда которого (яркость) пропорциональна поглощению19. Можно рассматривать каждую точку полученного сигнала как сумму значений коэффициентапоглощения вдоль того луча из пучка, который пространственно соответствуетименно этой точке (такую сумму часто называют суммой (интегралом) вдоль луча20). Положение дел таково, что полученный одномерный сигнал поглощения —это вся информация, которую мы в данный момент имеем об объекте.По одной проекции у нас нет способа определить, имеем ли мы дело толькос одним объектом или со многими, расположенными вдоль направления лучей; но мы начнем процесс реконструкции с построения изображения, базирующегося именно на такой информации. Подход состоит в проецированииодномерного сигнала назад вдоль направления распространения пучка лучей,как показано на рис.