Восстановление изображений

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



и т.п. Устранение или ослабление искажений с целью повышения резкости относится к задаче восстановления изображений.

Наиболее общая схема формирования изображения представлена на рис. 4.1,



Рис.4.1. Схема формирования изображения



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

.                                          

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

4.1. Модели изображений и их линейных искажений

4.1.1. Формирование изображений

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

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



С учетом вышеизложенного наблюдаемое нерезкое изображение  можно представить как выход линейной системы, показанной на рис. 4.2,



Рис.4.2. Линейная модель формирования изображения

 а  математическая модель процесса его формирования  имеет вид:

,

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



         (4.1)

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

В выражении (4.1) учтено, что изображения, встречающиеся в практических задачах, имеют конечные размеры. Это означает, что яркость изображения полагается равной нулю всюду, кроме некоторой конечной области, которую будем называть кадром и обозначать через . Верхний индекс в кадре  соответствует символу изображения или ФРТ, для которого определен этот кадр. Например, кадр изображения  будем обозначать через . Мы будем рассматривать только прямоугольные кадры, стороны которых параллельны координатным осям. Под размерами кадра будем понимать совокупность длин вертикальной  и горизонтальной  сторон кадра . Часто оказывается удобным совместить центр кадра  изображения с началом координат на плоскости  изображения, тогда

                     (4.2)

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



Рис.4.3. Относительные размеры изображения и ФРТ

Кадр  функции рассеяния точки представляет собой минимальный прямоугольник, содержащий все точки с координатами , в которых величина абсолютного значения ФРТ заметно отлична от нуля, например, когда  . Следует отметить, что кадр  ФРТ не всегда симметричен относительно начала координат. В частности, симметрией ФРТ не обладают каузальные линейные системы. Действие ФРТ сводится к тому, что каждая точка исходного изображения  «размазывается» в некоторую область,  ограниченную кадром . Кадр  исходного изображения может быть построен путем перемещения  и представляет собой область всех точек, охватываемых кадром   при его перемещении по кадру  наблюдаемого изображения.


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

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

,        (4.3)

где

.      (4.4)

Аргументы с индексом 1 обозначают номер строки, а с индексом 2 - номер столбца. В дискретном случае размеры кадра (число отсчетов)  и  определяются отношением длин соответственно вертикальной  и горизонтальной сторон кадра аналогового изображения к величине интервала дискретизации .

Операция свертки, которая имеется в формулах (4.1) и (4.4), эквивалентна произведению в частотной области. Это позволяет выполнить быструю имитацию линейных искажений с помощью ДПФ, заменив обычную свертку циклической (смотри главу 3). Как правило, размеры  кадра ФРТ много меньше размеров кадра исходного изображения, поэтому перед преобразованием массив  должен быть дополнен нулями. Кроме того, полагается, что изображения и ФРТ являются периодически продолженными, которые так же как и в главе 3 обозначаются волнистой линией.

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

,                              (4.5)



где  и  - пространственные частоты. Размеры кадра изображения  , полученного после обратного ДПФ от  , равны размерам кадра  исходного изображения. Для завершения процедуры имитации необходимо «обрезать» края изображения   до размеров кадра  и добавить аддитивную помеху . Здесь и далее результаты преобразования Фурье от функций в пространственной области будут обозначаться соответствующими прописными буквами.

Рассмотрим импульсные и частотные характеристики формирующих систем при смазе и расфокусировке.

4.1.2. Размытие вследствие движения (смаз)

Смаз изображения возникает при взаимном движении камеры и объекта относительно друг друга во время экспозиции. Наблюдаемое изображение окажется как бы результатом наложения со смещением множества исходных изображений. Мы рассмотрим только  тот случай, когда камера перемещается с постоянной горизонтальной скоростью относительно снимаемого объекта. ФРТ  и передаточная функция такой системы определяются выражениями:

                            

,

где длина смаза  - равна произведению скорости движения камеры на время экспозиции.  Соответственно в дискретном случае ФРТ  смаза равна

                (4.6)

где размеры кадра  и . Здесь квадратные скобки обозначают операцию округления до целого. Взаимное расположение кадров изображений и ФРТ при смазе показано на рис. 4.4.

Рис.4.4. Взаимное расположение  изображения и ФРТ при смазе	Рис.4.5. Изображние модуля частотной характеристики искажающей системы

Дополняя ФРТ (4.6) нулями до размеров кадра исходного изображения и применяя двумерное ДПФ, получим частотную характеристику искажающей системы:

       (4.7)

Изображение модуля  приведено на рис.4.5 при   и размерах исходного изображения элементов.

На рис.4.7 приведен искаженный вариант исходного изображения «Сатурн» (рис.4.6). Горизонтальный смаз составляет 15 элементов. Исходное изображение содержит  элементов, а искаженное -  элементов.

Рис.4.6.Исходное изображение “Сатурн”	Рис.4.7. Смазанное изображение “Сатурн”

4.1.3. Расфокусировка

Четкость изображения характеризуется воспроизведением мелких деталей и определяется разрешающей способностью формирующей системы. Разрешающая способность, например, оптической системы численно выражается количеством пар черно-белых линий на 1 мм изображения, которое формируется объективом системы. Если плоскость формируемого изображения находится в фокусе объектива, то пучок лучей,  исходящий от точки на объекте, сходится в точку на изображении. При расфокусировке точка воспроизводится в виде некоторого пятна (кружка  размытия), и две близко расположенные точки на исходном изображении сливаются в одну на наблюдаемом. Величина кружка размытия зависит от фокусного расстояния объектива, а также  от расстояний от объектива до объекта и до плоскости формируемого изображения [4.1]. Дискретное изображение будет четким (сфокусированным), если диаметр кружка размытия не превышает шага дискретизации  наблюдаемого изображения. В противном случае линейные искажения становятся заметными.

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

                            (4.8)

Из (4.8) следует, что размеры кадра . Взяв двумерное преобразование Фурье от  (4.8), получим передаточную функцию оптической системы

,                (4.9)

где  - функция Бесселя первого порядка.

В дискретном случае ФРТ (4.8) имеет вид:

                     (4.10)

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

Рис.4.8. ФРТ  тонкой линзы	Рис.4.9. Изображение модуля частотной характеристики тонкой линзы

<


Земную атмосферу также можно рассматривать как оптическую систему.  В качестве приближенной модели ФРТ такой системы используется двумерный гауссовский импульс

,                                (4.11)

который в дискретном случае имеет вид

,                                (4.12)

где -  нормирующий коэффициент,  - коэффициент пространственной нерезкости. Передаточная функция, соответствующая ФРТ (4.11), определяется выражением

.                  (4.13)

Очевидно, что точки, для которых выполняется условие (4.2), образуют круг радиусом

 .                                               (4.14)

Следовательно, чем больше , тем меньше расфокусировка наблюдаемого изображения. ФРТ для земной атмосферы и соответствующая ей передаточная функция при  приведены на рис.4.10 и 4.11. Радиус кружка размытия примерно равен 15. Размеры пятна ФРТ на рис. 4.10 визуально кажутся меньше чем размеры пятна для тонкой линзы (рис. 4.8), т.к. гауссовский импульс является быстро убывающей функцией.

Рис.4.10. ФРТ атмосферы Земли	Рис.4.11. Изображение модуля частотной характеристики атмосферы Земли

Соотношения между различными кадрами  изображений при расфокусировке  соответствуют тем, что приведены на рис. 4.3, поскольку кадр ФРТ   симметричен   относительно   центра  координат.  На  рис.4.12  приведен

Рис.4.12. Дефокусированное изображение “Сатурн”

искаженный вариант изображения «Сатурн» (рис.4.6). Свертка исходного изображения производилась с гауссовским импульсом при . Искаженное изображение содержит   элементов.

Таким образом,

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

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


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

2. Яркость на краях кадра искаженного изображения зависит от яркости объектов, расположенных вне кадра, за счет свертки исходного изображения с ФРТ. При восстановлении изображений из-за неполной информации о сигнале вне кадра возникают краевые эффекты. Влияние краевых эффектов на качество восстановления в ряде случаев оказывается даже более существенным, чем зашумленность изображения.

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

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

4.2. Алгебраические методы восстановления изображений



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



Удобно представить соотношения (4.3) и (4.4) в матричной форме, используя лексикографическое упорядочивание. Для этого двумерный массив наблюдаемого изображения  размера  преобразуем в  вектор-столбец  размера . Преобразование осуществляется разверткой массива  по строкам. Аналогичным образом преобразуются в  вектор-столбцы ,  и  искаженное при отсутствии шума изображение , исходное изображение  и  шум . Размеры векторов  ,  и  равны . Связь между лексикографически упорядоченными изображениями определяется соотношением

,                                        (4.15)

где искаженное изображение

.                                            (4.16)

Символ  обозначает прямоугольную матрицу размером , с помощью которой вектор  исходного изображения  преобразуется в искаженное изображение . Матрица  имеет блочную структуру [4.2], элементы которой представляют собой  отсчеты ФРТ. Задачи восстановления изображений алгебраическими методами при наличии и отсутствии шумов наблюдения имеют качественные различия.

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

 .                                                  (4.17)

Если бы  была квадратной матрицей и существовала бы обратная матрица , то, очевидно, что решение системы имело бы вид

.                                                 (4.18)

Однако матричное уравнение (4.16) представляет собой недоопределенную

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


Среди всех возможных решений недоопределенной разрешимой системы (4.16) в качестве оценки  выбирается решение, минимизирующее норму ошибки восстановления

,                       (4.19)

где  - символ транспонирования;  - вектор ошибки восстановления. Критерий (4.19) называется критерием наименьших квадратов. Доказано [4.3], что норма ошибки будет минимальной, если оценка

,

где  - обобщенная обратная матрица. В общем случае норма ошибки не равна нулю.

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

Для  искаженных изображений, наблюдаемых в присутствии шумов, к элементам вектора-столбца  добавляются отсчеты вектора-столбца . Это делает систему уравнений, как правило, неразрешимой. Неразрешимость системы означает, что не существует оценки исходного изображения, при которой она перейдет в тождество. Можно найти лишь приближенное решение  неразрешимой системы, которое определяется из условия минимума нормы ошибки  [4.4, 4.5]

 .                  (4.21)

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

Таким образом, в обоих рассмотренных случаях обобщенное обращение матриц дает оптимальное решение, удовлетворяющее критериям  наименьших квадратов (4.19) или (4.21).


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

Основным недостатком алгебраических алгоритмов восстановления изображений является необходимость выполнения трудоемких операций обращения, умножения и транспонирования матриц огромных размеров. Напомним, что размер матрицы  равен произведению числа отсчетов исходного и наблюдаемого изображений. Кроме того, обращение матриц больших размеров представляет собой трудную задачу. Эта задача значительно упрощается, если искаженное изображение формируется из исходного путем циклической свертки с ФРТ. К сожалению, для реальных задач восстановления изображений это условие не выполняется. Альтернативой алгебраическим методам являются методы линейной фильтрации изображений.

4.3. Методы восстановления изображений на основе

 пространственной фильтрации

 

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

При циклической свертке модель (4.4) формирования искаженного изображения определяется соотношением



      (4.22)

где  - кадр изображения, одинаковый для всех изображений и ФРТ, входящих в (4.22). Размеры кадра равны периоду повторения изображений и ФРТ.

Применяя к (4.22) ДПФ, получим

        (4.23)

 

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

                              (4.24)

исходного изображения . В пространственно-частотной области спектр оценки с учетом (4.24) можно записать как

         (4.25)

4.3.1. Инверсный фильтр

Простейшим способом восстановления четкости изображения является обработка наблюдаемого изображения в пространственно-частотной области  инверсным фильтром [4.5]. Передаточная функция инверсного восстанавливающего фильтра определяется соотношением

.                                     (4.26)

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

                (4.27)

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

На рис. 4.13. и 4.14 приведены результаты восстановления изображений «Часы» и «Сатурн» инверсным фильтром: а) исходные изображения размером  элементов; б) дефокусированные изображения, полученные в результате свертки с гауссовским импульсом при  с последующим «обрезанием» краев до размеров  элементов; в) изображения, восстановленные инверсным фильтром.


Восстановить изображение «Часы» инверсным фильтром не удается из- за краевых эффектов. Практически идеальное восстановление изображения «Сатурн»  объясняется тем, что объекты наблюдаются на фоне постоянной яркости и расположены в центре кадра. В этом случае  изображения  и , полученные в результате обычной и циклической свертки с ФРТ, равны друг другу. Отметим, что при этих условиях алгебраический метод  также позволяет точно восстановить изображение. Однако при инверсной фильтрации процедура обращения матриц заменяется на более простую процедуру перемножения массивов в частотной области.

На рис. 4.15 и 4.16 приведены сечения типичных частотных характеристик ФРТ  и соответствующих им инверсных фильтров, из которых следует,  что модуль передаточной функции формирующей системы, как правило, стремится к нулю   на   высоких    частотах.   Кроме    того,   нули   в передаточной функции имеются в рабочей полосе частот при расфокусировке камеры (4.10) и смазе (4.6). В этом случае инверсный фильтр является сингулярным, т.к. модуль его передаточной функции становится бесконечно большим на пространственных частотах, соответствующих нулевым значениям модуля передаточной функции  искажающей системы. Причем наличие даже относительно слабого шума приводит к появлению интенсивных шумовых составляющих на выходе инверсного фильтра, полностью разрушающих изображение. Этот факт иллюстрируется рис.4.17.  К дефокусированному изображению «Сатурн» (рис. 4.14.б) был добавлен аддитивный дельта-коррелированный шум. Из восстановленного изображения видно, что даже при пренебрежимо малом уровне шума (отношение сигнал/шум ) метод инверсной фильтрации дает очень плохой результат.

а)	б)
в)
Рис.4.13. Результаты восстанвления изображения “Часы”

а)	б)
в)
Рис.4.14. Результаты восстанвления изображения “Сатурн”

	Рис.4.15. Частотные характеристики искажающей системы с цилиндрической ФРТ  и инверсного фильтра	Рис.4.16. Частотные характеристики искажающей системы с гауссовской ФРТ  и инверсного фильтра
Рис.4.17. Результат восстановления изображения “Сатурн” при

<


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

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

4.3.2. Фильтр Винера

Инверсная фильтрация обладает низкой помехоустойчивостью, потому что этот метод не учитывает зашумленность наблюдаемого изображения. Значительно менее подвержен влиянию помех и сингулярностей, обусловленных нулями передаточной функции искажающей системы, фильтр Винера (смотри главу 3), т.к. при его синтезе наряду с видом ФРТ используется информация о спектральных плотностях мощности изображения и шума.  При этом полагается, что изображение является реализацией случайного двумерного поля. Частотная характеристика восстанавливающего фильтра Винера, полученная для периодически продолженных изображений, с учетом (2.34) имеет вид [4.6]

 , (4.28)

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

Преобразуем передаточную функцию фильтра Винера  (4.28) следующим образом:

                   (4.29)



Анализируя соотношения (4.28) и (4.29), можно отметить следующее:

1. При отсутствии шума фильтр Винера переходит в инверсный фильтр. Следовательно, в области низких частот, где, как правило, отношение сигнал/шум велико, передаточные функции инверсного и винеровского фильтров практически совпадают.

2. При уменьшении спектральной плотности мощности исходного изображения передаточная функция фильтра Винера стремится к нулю. Для изображений это характерно на верхних частотах.

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

На рис. 4.18 приведены  одномерные сечения типичных передаточных функций винеровских фильтров (сплошная линия). Здесь же для сравнения приведены сечения передаточных функций инверсных фильтров (4.15) и (4.16), которые обозначены штриховой линией.

Рис.4.18. Частотный характеристики фильтра Винера при цилиндрической и гауссовской ФРТ

Рассмотрим результаты моделирования винеровского алгоритма восстановления. На рис. 4.19.а и 4.21.а приведены результаты искажения изображений «Сатурн» и «Часы» сверткой с гауссовской ФРТ  (  ) с последующим «обрезанием» краев и добавлением аддитивного дельта-коррелированного шума (). На рис. 4.20.а и 4.22.б приведены изображения, полученные в результате смаза () изображений «Сатурн» (рис. 4.6) и «Часы» (рис. 4.22.а) () также с последующим «обрезанием» краев и добавлением аддитивного дельта-коррелированного шума ().

а)	б)
Рис.4.19. Восстановление дефокусированного изображения “Сатурн” при

а)	б)
Рис.4.20. Восстановление смазанного изображения “Сатурн” при

Размеры всех наблюдаемых и восстановленных изображений равны  элементов. Результаты восстановления винеровским фильтром изображения «Сатурн» (рис. 4.19.б и рис.4.20.б) свидетельствуют о том, что фильтр Винера значительно лучше подавляет шумы.


Осциллирующая помеха на результатах восстановления изображения «Часы» (рис. 4.21.б и рис.4.22.в) вызвана краевыми эффектами. Очевидно, что ее уровень существенно меньше, чем при инверсной фильтрации (см. рис.4.13.в). Однако винеровский фильтр лишь частично компенсирует краевые эффекты, которые делают качество восстановления неудовлетворительным.

а)	б)
Рис.4.21. Восстановление дефокусированного изображения “Часы” при
а)	б)
в)
Рис.4.22. Восстановление смазанного изображения “Часы” при

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

4.3.3. Компенсация краевых эффектов при

восстановлении линейно-искаженных изображений

          На восстановленных изображениях, приведенных на рис.4.13.в, 4.21.б и 4.22.в, присутствует осциллирующая помеха большой интенсивности, которая возникает из-за того, что инверсный фильтр и фильтр Винера были синтезированы без учета ограниченных размеров наблюдаемых изображений. Вследствие того, что искаженное изображение записывается в кадре конечного размера, в усеченном изображении происходит потеря информации, содержащейся в исходном изображении вблизи границ. Поэтому при коррекции линейных искажений усеченного изображения возникают ложные детали в виде ряби или полос, интенсивность которых особенно велика при цилиндрической форме ФРТ и равномерном смазе.

          К сожалению, решить уравнение Винера-Хопфа для сигналов и изображений, наблюдаемых на ограниченном интервале, не удается. Поэтому отсутствуют оптимальные пространственно-инвариантные фильтры, учитывающие краевые эффекты. Для компенсации краевых эффектов используются различные эвристические алгоритмы.


Некоторые из них будут рассмотрены в данном подразделе.

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

При дефокусировке функция окна является разделимой относительно пространственных координат: . При вертикальном или горизонтальном смазе используется одномерная функция окна, на которую умножаются соответственно столбцы или строки наблюдаемого изображения. Известно много одномерных функций окна , которые могут быть использованы при восстановлении изображений, например, окна Бартлетта , Кайзера, Блэкмана и др [4.7].

          Хорошие результаты дает функция окна [4.8]

,       (4.30)

форма которой определяется двумя независимыми параметрами   и . Параметр  влияет на размеры окна, а   - на скорость спада краев окна к нулю.

           На рис 4.23 и 4.24 приведены результаты восстановления изображения «Часы» при горизонтальном смазе, где а - результаты умножения строк искаженного изображения, приведенного на рис.4.22.б, на окно Кайзера и окно (4.30); б - результаты восстановления фильтром Винера. Параметры окон подбирались, исходя из визуального качества восстанавливаемых изображений.

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

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


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

а)	б)
Рис. 4.23. Восстановление с использованием окна Кайзера
а)	б)
Рис. 4.24. Восстановление с использованием окна (4.30)

                 (4.31)

корреляционной функции неограниченного изображения [4.9]. В этом случае спектральная плотность мощности усеченного изображения равна свертке спектральной плотности мощности неограниченного изображения и спектральной плотности окна (4.31). Подставляя соответствующие спектральные плотности мощности в уравнение Винера-Хинчина  и решая его, получим коэффициент передачи фильтра для усеченного изображения [4.10]

 ,                  (4.32)

где   - спектральная плотность окна (4.31). Следует подчеркнуть, что импульсная характеристика фильтра (4.32) не сводится к произведению импульсной характеристики фильтра Винера и  регуляризирующего двумерного треугольного окна (4.31).

На рис. 4.25 приведен результат восстановления изображения «Часы» фильтром (4.32), откуда следует, что фильтр (4.32) практически полностью компенсирует краевые эффекты. Это позволяет отказаться от предварительной обработки. Качество восстановления изображения в центре и на краях почти одинаковое. Параметры фильтра (4.32) полностью определяются исходными данными и не требуют дополнительной подстройки. При использовании быстрого преобразования Фурье для обработки изображений объем вычислений при реализации фильтра (4.32) такой же, как и для фильтра Винера (4.28).

Рис.4.25. Результат восстановления с компенсацией краевых эффектов

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



          Простейшей является процедура одномерной экстраполяции. Яркость изображения вдоль строк и столбцов за пределами кадра наблюдаемого изображения  задается в виде полинома

.

Коэффициенты  определяются исходя из требований, предъявляемых к свойствам функции яркости. Например, на границе кадра  функция яркости должна равняться нулю, она должна быть неотрицательной, максимальное значение экстраполирующей функции не должно превышать максимального значения наблюдаемого изображения и т.п. Метод экстраполяции иллюстрируется рис.4.25, где а - экстраполированное изображение; б - результат восстановления. Размер наблюдаемого изображения «Часы» (см. рис. 4.22.б)  равен  элементов, экстраполированного - элементов. В качестве экстраполирующей функции использовался полином первой степени.

а)	б)
Рис. 4.26. Восстановление с применением экстраполяции

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

          Улучшить качество восстановления можно, используя одновременно экстраполяцию наблюдаемого изображения  и фильтра  (4.32). Результат такой комбинированной процедуры приведен на рис 4.27. Рассмотренные методы восстановления являются линейными. Их широкое использование обусловлено достаточно простыми методами синтеза и анализа линейных систем, а также высокой вычислительной эффективностью. Однако эти методы не являются оптимальными и не всегда обеспечивают эффективную компенсацию искажений. Линейная обработка является лишь приближением к оптимальной обработке, т.к. статистические характеристики подавляющего большинства изображений являются негауссовскими.


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



Рис.4.27. Восстановление с использованием экстраполяции и

компенсации краевых эффектов

4.4. Итерационные методы восстановления изображений

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

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

.                           (4.33)

Представим передаточную функцию инверсного фильтра   в виде геометрической прогрессии:

.                     (4.34)

Подставляя (4.34) в (4.33), получим

          (4.35)

Соотношение (4.35) позволяет представить процедуру нахождения оценки  в виде последовательных приближений:

                                                                                (4.36)









где каждое последующее приближение вычисляется по предыдущему. Взяв преобразование Фурье от соотношений (4.36), получим итерационную процедуру Ван Циттера [4.11]:



(4.37)



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


В итерационном алгоритме (4.37) нахождение обратного оператора заменяется на многократное вычисление свертки.

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

 .                                           (4.38)

Условие (4.38) выполняется для гауссовской ФРТ.  При цилиндрической ФРТ и равномерном смазе соотношение (4.33) заменяют на эквивалентное соотношение

.

Тогда итерационный алгоритм  (4.37) имеет вид [4.6]



(4.39)



где  и  - импульсные характеристики фильтров с передаточными функциями   и   соответственно. Свертка в (4.37) и (4.39) может быть выполнена с помощь БПФ в предположении, что изображения и импульсные характеристики являются периодически продолженными.

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



Рис. 4.28. Частотные характеристики итерационного фильтра на разных шагах



На рис 4.29  приведены результаты восстановления изображения «Часы», где а и б - повторно приведенные исходное (рис.4.22.а) и искаженное в результате смаза (рис 4.22.б) изображения; в - восстановленное изображение итерационным алгоритмом (4.37) ( число итераций ); г - результат восстановления по экстраполированному наблюдаемому изображению 4.26.а ( число итераций ). В качестве критерия остановки итеративного процесса использовался критерий минимума нормированной среднеквадратической ошибки оценивания:

,(4.40)

а)	б )
в)	г)
Рис. 4.29. Восстановление изображения “Часы” итерационным алгоритмом

где кадр  размером элементов расположен в центре кадра  наблюдаемого изображения. Для рис 4.29.в  , а для рис.4.29.г . Таким образом, итерационный алгоритм практически полностью компенсирует краевые эффекты. Кроме того, применение экстраполяции позволяет снизить среднеквадратическую ошибку оценивания почти в 2 раза.

          Наряду с описанными выше свойствами итерационные алгоритмы могут быть легко преобразованы в нелинейные путем введения нелинейных ограничений для восстанавливаемого изображения [4.6, 4.11]. Ограничения формулируются на основе априорных данных о форме или структуре объектов на исходном изображении. К априорным данным относятся такие свойства изображения, как неотрицательность яркости, ее верхний и нижний пределы, минимальная мощность сигнала, ограниченная пространственная и спектральная протяженность  и. т.п.

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

Итерационный алгоритм, например (4.39), с ограничением имеет вид



(4.41)

,

где  оператор ограничения.

Например, если используется оператор ограничения на неотрицательность  , то из  (4.41) следует, что в тех областях, где яркость оценки  меньше нуля, изменение оценки не происходит.


Нелинейный итерационный алгоритм (4.41) будет сходится, если сходится линейный алгоритм (4.39)  и оператор  является нерасширяющимся оператором. Для пространства сигналов  это означает, что действие оператора на изображение не должно приводить к увеличению его энергии. Очевидно, что к нерасширяющимся операторам относится оператор ограничения на неотрицательность, а также оператор ограничения диапазона, который определяется следующим соотношением

 

                      (4.42)

Для большинства цифровых изображений диапазон изменения яркости равен . Использование нелинейного алгоритма с ограничением диапазона для восстановления изображения «Часы» (рис. 4.29.б) обеспечивает уменьшение среднеквадратической ошибки до 4%. Особенно эффективен этот алгоритм при восстановлении изображений  с распределением яркости, близким к бинарному. На рис.4.30 приведены результаты восстановления изображения «Текст», где а - исходное изображение размером  эл.;  б - часть исходного изображения, попадающая в кадра  размером   эл.; в - наблюдаемое изображение размером  эл., полученное в результате   смаза исходного изображения (, ) ; г - экстраполированное  изображение  размером      эл.; д   и   е   - изображения,  восстановленные по экстраполированному изображению итерационным линейным алгоритмом () и итерационным алгоритмом с ограничением диапазона яркости (). Ошибка вычислялась по кадру   размером  эл. На рис.4.30 размеры изображений увеличены в полтора раза.

а)	б)
в)	г)
д)	е)
Рис. 4.30.Восстановление изображения “Текст” нелинейным итерационным алгоритмом

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


Содержание раздела

---

---