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

ВОССТАНОВЛЕНИЕ ДЕФОКУСИРОВАННЫХ ИЗОБРАЖЕНИЙ МЕТОДОМ ДВУМЕРНОГО ПРЕОБРАЗОВАНИЯ ФУРЬЕ

Рассматривается задача реконструкции дефокусированных изображений. Задача сводится к решению двумерного интегрального уравнения Фредгольма I рода типа свертки. Для решения прямой задачи (моделирования дефокусированного изображения) используется двумерное преобразование Фурье. Для решения обратной задачи (реконструкции дефокусированного изображения) используется метод преобразования Фурье и метод регуляризации Тихонова. Приведены численные результаты с иллюстрациями.

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

Реконструкция дефокусированных изображений является одной из актуальных задач цифровой обработки изображений [1-5]. Данная задача обычно описывается двумерным интегральным уравнением Фредгольма I рода типа свертки [6-13]. Одним из распространенных способов решения этой некорректной задачи (как прямой, так и обратной) является использование преобразования Фурье (ПФ) с фильтрацией Винера, Тихонова и т.д. (для устойчивости решения). При этом края дефокусированного изображения обычно формируются с использованием так называемых граничных условий [5, с. 108], [14, 15]. Их введение обусловлено тем, что на практике для регистрации изображений применяются системы с матрицами, размер которых совпадает с размером области фокусировки. В данной работе отмечается искусственность граничных условий при формировании краев дефокусированного изображения, а также неадекватность описания физической задачи дефокусировки с помощью ПФ. Разработан ряд программ на МЛТЬЛВ'е для решения прямой и обратной задач. Производится сравнение разработанных алгоритмов и программ с классическими методами (параметрической фильтрации Винера и др.) и соответствующими программами реконструкции ^есоп^'ипг.т, deconvlucy.m, deconvreg.m), а также дана качественная и количественная оценка погрешностей полученных результатов.

Математическое описание задачи получения дефокусированного изображения

Изображение можно описать функцией /(£, п) двух пространственных координат п, представляющей интенсивность или яркость изображения в каждой точке (£, п) плоскости. Математически задача восстановления изображений, искаженных линейной пространственно-инвариантной системой (однородной системой), описывается интегральным уравнением ФредгольмаIрода типа свертки [4-13, 16, 17]:

И РЕГУЛЯРИЗАЦИИ ТИХОНОВА В.В. Шемплинер Научный руководитель - д.т.н., профессор В.С. Сизиков

Введение

| IЦ х-£, у-п) / (£, п) ^ (1ц = g( х, у), -ю< х, у <+<ю , (1)

где g(х, у) - искаженное изображение на выходе системы; /п) - исходное изображение на входе системы; к(х - у - п) - импульсная характеристика системы, или функция рассеяния точки (ФРТ), являющаяся трансляционно-инвариантной (разностной) [4, с. 380].

Рассмотрим физический принцип получения расфокусированного изображения. Считаем, что фотографируемый объект (для простоты полагаемый плоским) и фотопленка фотоаппарата расположены параллельно апертуре линзы фотоаппарата по разные стороны от нее на расстояниях соответственно, / и /2 +5 от линзы (5 - погрешность фокусировки изображения), причем

V /+V I = V/, (2)

где / - фокусное расстояние линзы, /1 > /2. В результате на фотопленке возникнет перевернутое изображение (см. рис. 1).

Рис. 1. Модель расфокусировки изображения

Введем в плоскости объекта прямоугольную систему координат £ 'O 'ц ', на «идеальной» фотопленке, расположенной «в фокусе» (S = 0) - систему координат £ О "ц ", а на реальной фотопленке, расположенной «вне фокуса» (S Ф 0) - систему координат £Оц, а также совпадающую с ней xOy. Обозначим через w'(£', ц') интенсивность, исходящую из некоторой точки A'(£' , ц' ) объекта. Точка A' отобразится на «идеальной» фотопленке также в точку A" с интенсивностью w "(£ '', ц") = w '(£ ', ц' ) и с координатами £ " = -£ '/q, ц " = - ц '¡q, где q = fi / f 2. На реальной же фотопленке A' отобразится в круг радиуса

ab

р=h

с центром в точке A(x,y), где a - радиус апертуры линзы.

Опишем математически задачу дефокусировки. Рассмотрим, помимо круга с центром в точке А(х,у), некоторый другой круг с центром в точке (|,п) (см. рис. 1). Радиусы этих

2

(и других) кругов одинаковы и равны р, а их площади равны $ = пр . В результате некоторая интенсивность и'(! п), соответствующая точке (|, п), будет «размазана» по кругу

22 радиуса р и площади $ = пр с плотностью интенсивности ч,(! п) / пр (постоянной, в

первом приближении, в пределах круга). Интенсивность в точке А(х, у) будет результатом

суммирования (интегрирования) по всем тем кругам, которые накрывают точку А(х, у).

Условие накрытия точки А(х, у) кругом с центром в точке (|, п) и радиуса р есть

х-|)2 + (у-п)2 <р, (4)

в результате интенсивность в точке А( х, у) будет равна

е (х, у) = Ц ^рп) ^ ап. (5)

х-I)2 + (у+п)2 <р Пр Запишем (5) в виде

Я ^^ 0! ап = Е( х, у). (6)

л/сх-I)2 + (у+п)2 <р пр Соотношение (6) есть двумерное интегральное уравнение I рода относительно '№(!, п) [8-10,15]. Однако оно записано не в стандартной форме. Приведем его к стандартной форме. Запишем (6) в виде

ж ж

I | к(х -I, у -пМ!, п)0|0п = Е(х, у), - ж < х, у <ж , (7)

—да —да

где

или

IVпр2, V(x — + (У + n)2 <Р,

k(х — y — n) = Г пр , V(х — + (y + n) <р, (8)

0, иначе

k(X, y) J1/пр2, Vх2 + y 2 <р, (9)

[ 0, иначе.

Соотношение (7) есть двумерное интегральное уравнение Фредгольма I рода типа свертки [6-13]. В нем g(х,y) - интенсивность в плоскости реальной (расположенной «не в фокусе») фотопленки, k (х, y) - ядро интегрального уравнения, причем р определяется согласно (3), где a и f2 известны, a 5 (или р) может быть определено путем подбора. Ядро интегрального уравнения k(х,y) называется функцией рассеяния точки (ФРТ, PSF). Наконец, n) есть искомая интенсивность (которая была бы на снимке при 5 = 0).

На рис. 2 представлен результат применения к исходному изображению (томограмме-фантому) искажающего оператора с ядром k (х, y) (9).

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

Рис. 2. Исходное изображение (а). Результат дефокусировки исходного изображения с использованием искажающей функции в виде диска радиуса 2 (Ь), 5 (с) и 10 (ф пикселей, соответственно

Для моделирования дефокусированного изображения больше подходит методика Гонсалеса и др. [4, 5], когда обрезаются размытые края и решается сложная краевая задача. Если же наблюдательное устройство позволяет зарегистрировать дефокусированный снимок с размытыми краями (большего размера, чем исходное изображение), то более эффективным представляется алгоритм обработки с учетом размытых краев. При этом не нужно решать сложную краевую задачу, так как информация на краю снимка не потеряна.

Реконструкция дефокусированного изображения методом преобразования Фурье и регуляризации Тихонова

Ниже рассмотрена фильтрация методом регуляризации Тихонова, которая в некоторых источниках [4, с. 395] называется методом минимизации сглаживающего функционала со связью, или методом наименьших квадратов со связью [5, с. 187]. В этом методе регуляризация задачи достигается заменой исходной задачи на задачу нахождения экстре-

мума (минимума) некоторого сглаживающего функционала. В качестве такого функционала C[ f ] можно использовать квадрат нормы лапласиана [4, с. 397] M-1N-1/ V2

C[f ] = I I (v2 f(X, y)) (10)

x=0 y=0

с дополнительным ограничением (связью) вида

iig - h f ii2=ihii2. (ii)

2

Норма (10) соответствует норме в пространстве Соболева W2 [16, с. 501]. В отечественной литературе [11, с. 192] чаще используют норму в пространстве L2 . В методе регуляризации Тихонова ставят два условия: условие минимизации невязки типа

iig - Hf ii2 = min (12)

f

и условие минимизации нормы решения типа

iifii2 = min. (13)

f

Рассмотрим одномерное интегральное уравнение Фредгольма I рода типа свертки:

да

Hw = J k(x-£) w(£) d£ = f (x), -да< x <да . (14)

—да

Применительно к нему в методе регуляризации Тихонова решение находится из условия минимума сглаживающего функционала (10):

да да

J [Hw- f (x)]2dx + a J M(o) i Y(o) i2 do = min, (15)

w

-да -да

где

M(o) = i o i2q (16)

- регуляризатор q-го порядка, причем q > 0 - порядок регуляризации, обычно q = 1.

Из условия (15) получается следующее выражение для регуляризованного решения

wa (£)=± да ^(-o)F(o) еdo, (17)

a 2п -да L(o) + a M (o)

где

L(o) =i A,(o) i2 = A,(o) A,(-o) = Re2 A,(o) + Im2 A,(o). (18)

Уравнение (1) как двумерное интегральное уравнение Фредгольма I рода типа свертки может быть решено методом двумерного ПФ (инверсная фильтрация). Решение записывается в виде ОПФ

да да

w(£,п) = —- [ [ Ж(и1,ш2)е—''(а^+^л)dш1dш2, (19)

4П2

—да—да

где ПФ (спектр) решения

W(»1,«2) = GPЩ. (20)

К (ШЬ Ш2)

где С(ш1,Ш2) и K(ш1,Ш2) - преобразования Фурье (спектры) правой части g(x,у) и ядра интегрального уравнения (2.21), равные

а(шьш2) = | | Е(х,у)в1 (Ю1 х+©2у), (21)

—да—да да да

К(©ьШ2) =| | к(х,у)в'4©1 х+©2у)ёхёу. (22)

—да —да

Ядро к(х,у) выражается аналитической формулой (9), поэтому К(ш1,©2)может быть найдено аналитически согласно (16). А С(©1,©2) (а также К(©1,Ш2)) должно быть найдено численно по стандартной программе двумерного ДПФ (обычно в виде БПФ).

Однако задача решения уравнения (1) является некорректной [1-3, 8, 13-16]. Это связано с тем, что функция ^ (х, у) измеряется с погрешностью, и это ведет к сколь угодно большим погрешностям решения п). Поэтому формулы (14)-(16) не годятся для устойчивого решения уравнения (1). Устойчивое решение уравнения (1) методами двумерного ПФ и регуляризации Тихонова имеет вид

дада

^а &П) = -Лу | | Ша (©ь©2)в-г'(Ю1^+Ю2П)^Ш2, (23)

4п —да—да

где

Жа (© 1,©2) = © 1 ч , Ь 2\ , (24)

Д©1, ©2) + а М (©!, ©2) 2

Ц©1,©2) =| К(©1,©2)| = К(©1,©2)К(-©1,-©2),

М (©1, ©2) = (©1 +©2 ) - регуляризатор, а > 0 - параметр регуляризации.

Решение (17)-(18) при правильно выбранных значениях а и 5 (или р) обладает устойчивостью и достаточной разрешающей способностью.

Выбор параметра регуляризации

Важным является вопрос о выборе значения параметра регуляризации а. Разработан ряд способов выбора а в методе регуляризации Тихонова [16, 19]. Выбор а можно осуществлять, например, способом невязки или обобщенным принципом невязки [8-9, 16, 19]. Разработаны также следующие способы выбора параметра регуляризации а: способ квазиоптимального (квазинаилучшего), способ отношения, способ независимых реализаций, способ перекрестной значимости, способ моделирования и др. [16, 20]. Однако для задачи реконструкции изображений, как показала практика, более эффективен способ подбора [10-13, 18]). В данной работе для задачи реконструкции дефокусированных изображений выбор а осуществлялся способом подбора. Согласно ему, для ряда значений а вычисляются решения на (£, у) по вышеизложенным формулам, они выводится на дисплей в графической форме и

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

оценкой следует использовать также количественную оценку среднеквадратического отклонения (СКО) регуляризованного решения wa (£, у) от точного w.

а

ННЙЯ 1

т

□

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

р= 3 (а) и р = 10 (с) методом регуляризации Тихонова с а = 3,5 • 10 5 (Ь, ф

Для количественной оценки погрешности метода реконструкции изображений использовалось относительное среднеквадратическое отклонение (СКО) восстановленного распределения плотности от точного распределения плотности [21]:

а

МЫ 0 2

I I — )

i =1 У =1

отн

1

МЫ 2

i =1} =1

(25)

0

где - значения восстановленной плотности интенсивности; /у - значения точной

плотности; М, N - соответственно, количество строк и столбцов матрицы плотности. При компьютерном моделировании выражения (25) при помощи средств системы программирования МЛТЬЛВ была составлена соответствующая программа КМББО.т. На рис. 3

представлен результат реконструкции дефокусированных изображений (фантомов томограмм, искаженных оператором к (х, у) (9) с р = 3 и р = 10 пикселей) методом преобразования Фурье и регуляризации Тихонова с параметром регуляризации а = 3,5 • 10—5 .

Сравнение методов реконструкции

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

Оценка эффективности того или иного метода производится как субъективно (визуально), так и с помощью расчета относительного СКО. На рис. 4 представлен результат восстановления дефокусированного изображения фотографа (с р = 10 пс) методом фильтрации Винера и методом ПФиРТ. На рис. 5 представлен результат восстановления дефокусированного изображения девушки ( р = 10 пс) методом фильтрации Люси-Ричардсона и методом ПФиРТ. На рис. 6 представлен результат восстановления дефокусированного изображения фотографа ( р = 10 пс) методом «слепой» деконволюции и методом ПФиРТ.

w -опдта1 ¡гладь (а) д - сЫисиэес! ¡гладе, гаа-10 (Ь)

Пегапб!гип1е|!лпаде [\\:\впегI (с) 9есопэ*гис*еп !гпаде [Т\к 1тапиV,:а!рИа_18 -3] Щ

Рис. 4. Результат реконструкции дефокусированного изображения: а) исходное изображение; Ь) дефокусированное изображение, р = 10;

с) метод Винера (аотн = 0,25); с1) метод ПФ и РТ, а = 10—3 (аотн = 0,22)

w - original irriaga (a) g - defocused image, rad=10 (b)

Reco nst mots d i mags [Lucy] (d) ReGüristmciedirnage;Tiktonü^^alpha-i6-S](d)

Рис. 5. Результат реконструкции дефокусированного изображения: a) исходное изображение; b) дефокусированное изображение, р = 10;

с) метод Люси-Ричардсона ( аотн = 0,14 ); d) метод ПФиРТ, а = 10 ( аотн = 0,15 )

Заключение

В результате выполнения данной работы были получены следующие результаты.

1. В работах [4-5, 14-15] при моделировании прямой задачи для расчета интенсивно-стей вблизи краев изображения используется такой прием, как граничные условия (boundary conditions). Однако в случаях, когда этот прием создает резкие края у дефокусированного изображения, возникают помехи (эффект Гиббса и т.д.) на реконструированном изображении. Для устранения таких помех предложено искусственно моделировать размытые края у дефокусированного изображения, но без использования граничного условия типа 'replicate', 'symmetric', 'circular' [5, с. 108], что повышает точность реконструкции.

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

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

Ш - опд!па1 ¡гладе (а) д - с^осиБес) ¡гпаде, гасЫО (Ь)

Рис. 6. Результат реконструкции дефокусированного изображения: a) исходное изображение; b) дефокусированное изображение, р = 10;

_3

с) метод «слепой» деконволюции ( аотн = 0,37 ); d) метод ПФиРТ, а = 10 ( аотн = 0,39 )

3. Рассмотрено два варианта алгоритма решения прямой задачи: с использованием операции накопления (суммирования) в пределах ФРТ и с использованием аппарата преобразования Фурье.

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

5. Разработаны программы в виде собственных m-функций в системе MATLAB.

Литература

1. Бейтс Р., Мак-Доннелл М. Восстановление и реконструкция изображений. - М.: Мир, 1989. - 336 с.

2. Прэтт У. Цифровая обработка изображений. - М.: Мир, 1982. - T. 2. - 792 с.

3. Грузман И.С., Киричук В.С., Косых В.П. и др. Цифровая обработка изображений в информационных системах. - Новосибирск: Изд-во НГТУ, 2000. - 168 с.

4. Гонсалес Р., Вудс Р. Цифровая обработка изображений. - М.: Техносфера, 2006. -

1072 с.

5. Гонсалес Р., Вудс Р., Эддинс С. Цифровая обработка изображений в среде MAT-LAB. - М.: Техносфера, 2006. - 616 с.

6. Арефьева М.В., Сысоев А.Ф. Быстрые регуляризирующие алгоритмы цифрового восстановления изображений // Вычислительные методы и программирование. -

1983. - Вып. 39. - С. 40-55.

7. Василенко Г.И., Тараторин А.М. Восстановление изображений. - М.: Радио и связь, 1986. - 304 с.

8. Тихонов А.Н., Гончарский А.В., Степанов В.В. Обратные задачи обработки фотоизображений // Некорректные задачи естествознания / Под ред. А. Н. Тихонова, А.В. Гончарского. - М.: Изд-во МГУ, 1987. - С. 185-195.

9. Бакушинский А.Б., Гончарский А.В. Некорректные задачи. Численные методы и приложения. - М.: Изд-во МГУ. - 1989. - 199 с.

10. Сизиков В.С., Белов И.А. Реконструкция смазанных и дефокусированных изображений методом регуляризации // Оптический журнал. - 2000. - Т. 67, № 4. - С. 6063.

11. Сизиков В.С. Математические методы обработки результатов измерений. - СПб.: Политехника, 2001. - 240 с.

12. Петров Ю.П., Сизиков В.С. Корректные, некорректные и промежуточные задачи с приложениями: Учебное пособие для вузов. - СПб.: Политехника, 2003. - 261 с.

13. Petrov Yu.P., Sizikov V.S. Well-Posed, Ill-Posed, and Intermediate Problems with Applications. - Leiden-Boston: VSP, 2005. - 234 p.

14. Lee K.P., Nagy J.G., Perrone L. Iterative methods for image restoration: a Matlab object oriented approach, 2002. http://www.matcs.emory.edu

15. Donatelli M. et al. Improved image deblurring with anti-reflective boundary conditions and re-blurring // Inverse Problems. - 2006. - Vol. 22. - P. 2035-2053.

16. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. - Киев: Наук. думка, 1986. - 544 с.

17. Воскобойников Ю.Е., Литасов В.А. Устойчивый алгоритм восстановления изображения при неточно заданной аппаратной функции // Автометрия. - 2006. - Т. 42, № 6. - С. 3-15.

18. Римских М.В., Евсеев В.О., Сизиков В.С. Реконструкция смазанных изображений различными методами // Оптический журнал. - 2007. - Т. 74, № 11. - С. 53-57.

19. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. - М.: Наука, 1990. - 232 с.

20. Воскобойников Ю.Е., Преображенский Н.Г., Седельников А.И. Математическая обработка эксперимента в молекулярной газодинамике. - Новосибирск: Наука,

1984. - 240 с.

21. Пикалов В.В., Непомнящий А.В. Итерационный алгоритм с вэйвлет-фильтрацией в задаче двумерной томографии // Вычислительные методы и программирование. -2003. - Т. 4. - С. 244-253.