Фильтрация изображений

Фильтрация изображений

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

Оптимальная линейная фильтрация. Уравнение Винера-Хопфа.

Пусть наблюдаемое изображение y есть некоторая функция от идеального сигнала x и помехи n. В точке с координатами (i, j): yij = f(xij, nij), i=0¸I-1, j=0¸J-1.

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

Если обе координаты некоторых точек окрестности выше (i, j), то фильтрация некаузальная. Если по одной из координат выполняет принцип каузальности, а по второй нет, то это полукаузальная фильтрация.

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

При линейной фильтрации выходной эффект определяется линейной комбинацией входных данных с весовыми коэффициентами a:

, где S – множество точек образующих окрестность.

Совокупность весовых коэффициентов a представляет собой двумерную импульсную характеристику. Если область S конечна, то импульсная характеристика имеет конечную длину и фильтр называется КИХ-фильтром. Если импульсная характеристика бесконечна, то БИХ-фильтром.

Легко заметить, что от координат рабочей точки (i, j) импульсная характеристика не зависит. Такая процедура обработки, не зависящая от координат, называется однородной.

В качестве критерия оптимальности для оценки фильтрации выбирается критерий минимума среднего квадрата ошибок:

.

Если вычислить производную этого выражения по коэффициенту a(k,l) и приравнять ее 0, то получим

(1)

Мат. ожидания, входящие в это выражение есть корреляционная и автокорреляционная функции:

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

Как изменяется средняя яркость изображения при фильтрации?

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

Рекуррентная каузальная фильтрация

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

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

где Dx - дисперсия последовательности, F – коэффициент ее одношаговой корреляции, параметр a - это ширина спектра последовательности.

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

(2)

Импульсная характеристика этого фильтра имеет экспоненциальный вид

Параметр 0<A<1 это коэффициент усиления фильтра Калмана.

Первое слагаемое в (2) определяет вклад в оценку сигнала на текущем i-ом шаге фильтрации, вносимый оценкой предыдущего шага, и называется одношаговым прогнозом. Второе слагаемое учитывает влияние текущего наблюдения yi и называется новой информацией. Коэффициент усиления А определяет чувствительность фильтра к этой новой информации. При высоком уровне шума параметр А приближается к 0. Это означает удлинение импульсной характеристики и сужение полосы пропускаемых фильтром частот. При снижении шума, наоборот А стремится к 1.

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

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

(3)

Для определения параметра А, воспользуемся уравнением Винера-Хопфа

Выражение в круглых скобках является ошибкой фильтрации

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

Подставив выражение ИХ (3) в выражение для отклика фильтра x*, получим

(4)

Рекуррентный характер (4) обеспечивает сокращение выполняемых операций. На каждом шаге требуется 3 операции умножения и 3 сложения. Структура алгоритма универсальна и не зависит от отношения сигнал/шум.

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

Некаузальная двумерная фильтрация

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

В уравнении Винера-Хопфа заменим конечную область S на бесконечную

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

Применим к обеим частям двумерное ДПФ, частотная характеристика .

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

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

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

Для винеровской фильтрации требуется знать спектральную плотность мощности изображения. Можно предварительно измерить характеристики по реальному изображению и ввести в ЭВМ в виде таблиц. А можно выбрать математическую модель изображения и по реальному изображению измерить параметры модели. В частности, часто используется модель гауссовского двумерного поля с корреляционной функцией Bx(i,j)= Dx F|i|+|j|, где Dx – дисперсия, а F- одношаговая корреляция.

Пространственная реставрация изображений

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

Общая схема формирования изображения: исходное, неизвестное изображение u(k, l) проходит через некоторый известный оператор искажений W (омега), в результате получается наблюдаемое изображение s(x, y)= W{ u(k, l) }. Задача реставрации найти оценку u - .

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

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

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

s(x,y) = z(x,y)+ n(x,y), где n(x,y) – аддитивный двумерный шум.

Изображение z(x,y) есть линейное искажение исходного изображения при отсутствии шума и определяется интегралом свертки

h(x,y) – двумерная импульсная характеристика или функция рассеяния точки искажающей системы.

Конечная область Q(.) определяет размер изображения и называется кадром. За пределами кадра яркость равна 0. Размеры кадра Lx, Ly. Кадр Qh функции рассеяния точки представляет собой минимальный прямоугольник, содержащий все точки с координатами (х, у), в которых величина абсолютного значения ФРТ заметно отлична от нуля.

Следует отметить, что кадр ФРТ не всегда симметричен относительно начала координат. В частности, симметрией ФРТ не обладают каузальные линейные системы. Действие ФРТ сводится к тому, что каждая точка исходного изображения u(k,l) «размазывается» в некоторую область, ограниченную кадром Qh.

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

Совместим центр кадра наблюдаемого изображения s с началом координат, тогда

.

Двумерные функции с непрерывными аргументами x,y, k,l заменим на массивы отсчетов с расстоянием D между узлами.

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

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

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

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

Импульсная характеристика и передаточная функция такой системы

где a - длина смаза равна произведению скорости движения камеры на время экспозиции.

В дискретном случае ФРТ смаза будет

, [] – целая часть числа.

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

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

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

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

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

3. Если передаточные функции искажающей системы имеют нули, то это приводит к полной утрате данных об исходном изображении на соответствующих частотах.

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

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

Линейное искажение исходного изображения при отсутствии шума в дискретном случае имеет вид

,

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

Преобразуем исходное изображение u размером Lu= (Lui +1) x (Luj +1) в вектор столбец u размером Lu x 1, т. е. каждый элемент столбца это целая строка. Аналогично преобразуются, искаженное при отсутствии шума изображение z, наблюдаемое изображение s и шум n. Их размеры совпадают Lz=Ls=Ln = (Lzi +1) x (Lzj +1). Тогда

s=z+n, z=hu, размеры матрицы h Lz x Lu. (5)

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

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

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

Если увеличить число отсчетов нерезкого изображения, то система будет переопределенной. Также переопределенную систему можно получить если выразить некоторые компоненты вектора идеального изображения через априорные данные. Например, известно, что объект расположен на черном фоне, тогда за пределами контура компоненты u можно считать нулевыми.

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

Норма ошибки будет минимальна, если

- обобщенная обратная матрица (6)

В общем случае норма ошибки не равна нулю.

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

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

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

Методы восстановления на основе пространственной фильтрации

Эти методы восстановления изображений реализуются с помощью ДПФ в частотной области. При этом обычная свертка заменяется циклической как в модели формирования искаженного изображения, так и в процедуре восстановления методом пространственной фильтрации. Все изображения полагаются периодически продолженными и имеют одинаковые размеры кадра L = Li х Lj. Размеры кадра равны периоду повторения изображений и ФРТ.

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

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

Применяя к этому выражению ДПФ, получаем

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

В пространственно-частотной области спектр оценки можно записать как

.

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

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

.

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

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

Частотные характеристики искажающей системы

с гауссовской ФРТ, инверсного фильтра и фильтра Винера

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

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

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

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

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

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

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

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

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

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

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

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

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

Другой путь – учесть ограниченные размеры изображения на этапе синтеза фильтра Винера. Корреляционная функция усеченного изображения может быть получена путем умножения корреляционной функции неограниченного изображения на следующее окно:

(Пси)

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

, Y - спектральная плотность окна.

Такой фильтр хорошо компенсирует краевые эффекты и не требует дополнительной подстройки.

Если требуется расширить границы изображения, чтобы извлечь больший объем информации, применяется процедура экстраполяции. Двумерную функцию яркости s(i, j) продолжают за границы кадра, так чтобы она была гладкой, а на границе экстраполированного кадра QЭ равнялась 0. Одномерная экстраполяция: яркость задается в виде полинома

f(i)=a0+a1i+a2i2+..., f(j)=a0+a1j+a2j2+...

Можно одновременно применять экстраполяцию и окно для усечения изображения.

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

Итерационные методы восстановления

Наиболее простой и широко распространенный итерационный метод Ван Циттера. Представим передаточную функцию инверсного фильтра в виде геометрической прогрессии:

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

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

Взяв обратное преобразование Фурье получим итерационную процедуру

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

При использовании итерационных алгоритмов встают 2 вопроса: сходится ли алгоритм и к какому решению он сходится. Алгоритм Ван Циттера сходится к оценке изображения при инверсной фильтрации если передаточная функция искажающей системы 0 < H(w1, w2) < 1.

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

Тогда итерационный алгоритм имеет вид

,

где h1 и h2 – импульсные характеристики фильтров с передаточными функциями H*(w1, w2) и |H(w1, w2)|2.

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

,

где Qe - кадр, расположенный в центре наблюдаемого изображения.

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

Если Y{} – оператор ограничения, то алгоритм с ограничением имеет вид

.

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

Оператор ограничения на неотрицательность:

Оператор ограничения диапазона: