Преобразование Фурье. Линейная фильтрация в частотной области
Линейная фильтрация изображений может осуществляться как в пространственной, так и в частотной области. При этом считается, что "низким" пространственным частотам соответствует основное содержание изображения - фон и крупноразмерные объекты, а "высоким" пространственным частотам - мелкоразмерные объекты, мелкие детали крупных форм и шумовая компонента.
Традиционно для перехода в область пространственных частот используются методы, основанные на $\textit$. В последние годы все большее применение находят также методы, основанные на $\textit$.
Содержание
Преобразование Фурье.
Преобразование Фурье позволяет представить практически любую функцию или набор данных в виде комбинации таких тригонометрических функций, как синус и косинус, что позволяет выявить периодические компоненты в данных и оценить их вклад в структуру исходных данных или форму функции. Традиционно различаются три основные формы преобразования Фурье: интегральное преобразование Фурье, ряды Фурье и дискретное преобразование Фурье.
Интегральное преобразование Фурье переводит вещественную функцию в пару вещественных функций или одну комплексную функцию в другую.
Вещественную функцию $f(x)$ можно разложить по ортогональной системе тригонометрических функций, то есть представить в виде
где $A(\omega )$ и $B(\omega )$ называются интегральными косинус- и синус-преобразованиями:
Ряд Фурье представляет периодическую функцию $f(x)$, заданную на интервале $[a,b]$, в виде бесконечного ряда по синусам и косинусам. То есть периодической функции $f(x)$ ставится в соответствие бесконечная последовательность коэффициентов Фурье
Дискретное преобразование Фурье переводит конечную последовательность вещественных чисел в конечную последовательность коэффициентов Фурье.
Пусть $\left\ < \right\>, i= 0,\ldots, N-1 $ - последовательность вещественных чисел - например, отсчеты яркости пикселов по строке изображения. Эту последовательность можно представить в виде комбинации конечных сумм вида
Основное отличие между тремя формами преобразования Фурье заключается в том, что если интегральное преобразование Фурье определено по всей области определения функции $f(x)$, то ряд и дискретное преобразование Фурье определены только на дискретном множестве точек, бесконечном для ряда Фурье и конечном для дискретного преобразования.
Как видно из определений преобразования Фурье, наибольший интерес для систем цифровой обработки сигналов представляет дискретное преобразование Фурье. Данные, получаемые с цифровых носителей или источников информации, представляют собой упорядоченные наборы чисел, записанные в виде векторов или матриц.
Обычно принимается, что входные данные для дискретного преобразования представляют собой равномерную выборку с шагом $\Delta $, при этом величина $T=N\Delta $ называется длиной записи, или основным периодом. Основная частота равна $1/T$. Таким образом, в дискретном преобразовании Фурье производится разложение входных данных по частотам, которые являются целым кратным основной частоты. Максимальная частота, определяемая размерностью входных данных, равна $1/2 \Delta $ и называется $\it$. Учет частоты Найквиста имеет важное значение при использовании дискретного преобразования. Если входные данные имеют периодические составляющие с частотами, превышающими частоту Найквиста, то при вычислении дискретного преобразования Фурье произойдет подмена высокочастотных данных более низкой частотой, что может привести к ошибкам при интерпретации результатов дискретного преобразования.
Важным инструментом анализа данных является также $\it$. Мощность сигнала на частоте $\omega $ определяется следующим образом:
Эту величину часто называют $\it$ на частоте $\omega $. Согласно теореме Парсеваля общая энергия входного сигнала равна сумме энергий по всем частотам.
График зависимости мощности от частоты называется энергетическим спектром или спектром мощности. Энергетический спектр позволяет выявлять скрытые периодичности входных данных и оценивать вклад определенных частотных компонент в структуру исходных данных.
Комплексное представление преобразования Фурье.
Кроме тригонометрической формы записи дискретного преобразования Фурье широко используется $\it$. Комплексная форма записи преобразования Фурье широко используется в многомерном анализе и в частности при обработке изображений.
Переход из тригонометрической в комплексную форму осуществляется на основании формулы Эйлера
$$ e^=\cos \omega t+j\sin \omega t, \quad j=\sqrt . $$
Если входная последовательность представляет собой $N$ комплексных чисел, то ее дискретное преобразование Фурье будет иметь вид
а обратное преобразование
Если входная последовательность представляет собой массив вещественных чисел, то для нее существует как комплексное, так и синусно-косинусное дискретное преобразование. Взаимосвязь этих представлений выражается следующим образом:
$$ a_0 =G_0 , \quad G_k =\left( \right)/2, \quad 1\le k\le N/2; $$
остальные $N/2$ значений преобразования являются комплексно сопряженными и не несут дополнительной информации. Поэтому график спектра мощности дискретного преобразования Фурье симметричен относительно $N/2$.
Быстрое преобразование Фурье.
Простейший способ вычисления дискретного преобразования Фурье (ДПФ) - прямое суммирование, оно приводит к $N$ операциям на каждый коэффициент. Всего коэффициентов $N$, так что общая сложность $O\left( \right)$. Такой подход не представляет практического интереса, так как существуют гораздо более эффективные способы вычисления ДПФ, называемые быстрым преобразованием Фурье (БПФ), имеющее сложность $O (N\log N)$. БПФ применяется только к последовательностям, имеющим длину (число элементов), кратную степени 2. Наиболее общий принцип, заложенный в алгоритм БПФ, заключается в разбиении входной последовательности на две последовательности половинной длины. Первая последовательность заполняется данными с четными номерами, а вторая - с нечетными. Это дает возможность вычисления коэффициентов ДПФ через два преобразования размерностью $N/2$.
Для $m < N/2$ тогда можно записать $G_m =G_ \left( m \right)+G_ \left( m \right)\omega _N^m $. Учитывая, что элементы ДПФ с индексом б ольшим, чем $N/2$, являются комплексно сопряженными к элементам с индексами меньшими $N/2$, можно записать $G_ =G_ \left( m \right)-G_ \left( m \right)\omega _N^m $. Таким образом, можно вычислить БПФ длиной $N$, используя два ДПФ длиной $N/2$. Полный алгоритм БПФ заключается в рекурсивном выполнении вышеописанной процедуры, начиная с объединения одиночных элементов в пары, затем в четверки и так до полного охвата исходного массива данных.
Двумерное преобразование Фурье.
Дискретное преобразование Фурье для двумерного массива чисел размера $M\times N$ определяется следующим образом:
а обратное преобразование
В случае обработки изображений компоненты двумерного преобразования Фурье называют $\textit$.
Важным свойством двумерного преобразования Фурье является возможность его вычисления с использованием процедуры одномерного БПФ:
Здесь выражение в квадратных скобках есть одномерное преобразование строки матрицы данных, которое может быть выполнено с одномерным БПФ. Таким образом, для получения двумерного преобразования Фурье нужно сначала вычислить одномерные преобразования строк, записать результаты в исходную матрицу и вычислить одномерные преобразования для столбцов полученной матрицы. При вычислении двумерного преобразования Фурье низкие частоты будут сосредоточены в углах матрицы, что не очень удобно для дальнейшей обработки полученной информации. Для перевода получения представления двумерного преобразования Фурье, в котором низкие частоты сосредоточены в центре матрицы, можно выполнить простую процедуру, заключающуюся в умножении исходных данных на $-1^$.
На рис. 16 показаны исходное изображение и его Фурье-образ.
Полутоновое изображение и его Фурье-образ (изображения получены в системе LabVIEW)
Свертка с использованием преобразования Фурье.
Свертка функций $s(t)$ и $r(t)$ определяется как
$$ s\ast r\cong r\ast s\cong \int\limits_^ r(t-\tau )d\tau . $$
На практике приходится иметь дело с дискретной сверткой, в которой непрерывные функции заменяются наборами значений в узлах равномерной сетки (обычно берется целочисленная сетка):
Здесь $-N$ и $P$ определяют диапазон, за пределами которого $r(t) = 0$.
При вычислении свертки с помощью преобразования Фурье используется свойство преобразования Фурье, согласно которому произведение образов функций в частотной области эквивалентно свертке этих функций во временн ой области.
Для вычисления сверки необходимо преобразовать исходные данные в частотную область, то есть вычислить их преобразование Фурье, перемножить результаты преобразования и выполнить обратное преобразование Фурье, восстановив исходное представление.
Единственная тонкость в работе алгоритма связана с тем, что в случае дискретного преобразования Фурье (в отличие от непрерывного) происходит свертка двух периодических функций, то есть наши наборы значений задают именно периоды этих функций, а не просто значения на каком-то отдельном участке оси. То есть алгоритм считает, что за точкой $x_$ идет не ноль, а точка $x_$, и так далее по кругу. Поэтому, чтобы свертка корректно считалась, необходимо приписать к сигналу достаточно длинную последовательность нулей.
Фильтрация изображений в частотной области.
Линейные методы фильтрации относятся к числу хорошо структурированных методов, для которых разработаны эффективные вычислительные схемы, основанные на быстрых алгоритмах свертки и спектральном анализе. В общем виде линейные алгоритмы фильтрации выполняют преобразование вида
$$ f'(x,y) = \int\int f(\zeta -x, \eta -y)K (\zeta , \eta ) d \zeta d \eta , $$
где $K(\zeta ,\eta )$ - ядро линейного преобразования.
При дискретном представлении сигнала интеграл в данной формуле вырождается во взвешенную сумму отсчетов исходного изображения в пределах некоторой апертуры. При этом выбор ядра $K(\zeta ,\eta )$ в соответствии с тем или иным критерием оптимальности может привести к ряду полезных свойств (гауссовское сглаживание при регуляризации задачи численного дифференцирования изображения и др.).
Наиболее эффективно линейные методы обработки реализуются в частотной области.
Использование Фурье-образа изображения для выполнения операций фильтрации обусловлено прежде всего более высокой производительностью таких операций. Как правило, выполнение прямого и обратного двумерного преобразования Фурье и умножение на коэффициенты Фурье-образа фильтра занимает меньше времени, чем выполнение двумерной свертки исходного изображения.
Алгоритмы фильтрации в частотной области основываются на теореме о свертке. В двумерном случае преобразование свертки выглядит следующим образом:
где $G$ - Фурье-образ результата свертки, $H$ - Фурье-образ фильтра, а $F$ - Фурье-образ исходного изображения. То есть в частотной области двумерная свертка заменяется поэлементным перемножением образов исходного изображения и соответствующего фильтра.
Для выполнения свертки необходимо выполнить следующие действия.
- Умножить элементы исходного изображения на $-1^$, для центрирования Фурье-образа.
- Вычислить Фурье образ $F(u,v)$, используя БПФ.
- Умножить Фурье образ $F(u,v)$ на частотную функцию фильтра $H(u,v)$.
- Вычислить обратное преобразование Фурье.
- Умножить вещественную часть обратного преобразования на $-1^$.
Как правило, фильтры описываются вещественными функциями, в этом случае каждый компонент $H$ умножается на соответствующие элементы действительной и мнимой части Фурье-образа изображения. Если исходная функция $f(x,y)$ и фильтр $H$ не комплексные, то результат свертки $g(x,y)$ также должен быть вещественной функцией. Однако на практике обратное преобразование содержит паразитную мнимую составляющую, которую надо отбросить.
Связь между функцией фильтра в частотной и пространственной области можно определить, используя теорему о свертке
Свертка функции с импульсной функцией может быть представлена следующим образом:
Фурье-преобразование импульсной функции
Пусть $f(x,y) = \delta (x,y)$, тогда свертка
Из этих выражений видно, что функции фильтра в частотной и пространственной областях взаимосвязаны через преобразование Фурье. Для данной функции фильтра в частотной области всегда можно найти соответствующий фильтр в пространственной области, применив обратное преобразование Фурье. То же верно и для обратного случая. Используя данную взаимосвязь, можно определить процедуру синтеза пространственных линейных фильтров.
- Определяем требуемые характеристики (форму) фильтра в частотной области.
- Выполняем обратное преобразование Фурье.
- Полученный фильтр можно использовать как маску для пространственной свертки, при этом размеры маски можно уменьшить по сравнению с размерами исходного фильтра.
$H(u,v)$ имеет вид $$H(u,v) = 1, \quad \mboxD(u,v) < D_0 ,$$ $$H(u,v) = 0, \quad \mboxD(u,v) \ge D_0 ,$$ где $D\left( \right)=\sqrt $ - расстояние от центра частотной плоскости.
После свертки с этим фильтром на результирующем изображении появляются паразитные искажения в виде полутоновых ложных границ.
получается путем инверсии идеального низкочастотного фильтра:
Здесь происходит полное подавление низкочастотных компонент при сохранении высокочастотных. Однако как и в случае идеального низкочастотного фильтра, его применение чревато появлением существенных искажений.
Для синтеза фильтров с минимальными искажениями используются различные подходы. Одним из них является синтез фильтров на основе экспоненты. Такие фильтры привносят минимальные искажения в результирующее изображение и удобны для синтеза в частотной области.
Широко используемым при обработке изображений является семейство фильтров на основании вещественной функции Гаусса.
Чем уже профиль фильтра в частотной области (чем больше $\sigma $), тем он шире в пространственной.
В двумерном случае фильтр гаусса выглядит следующим образом:
гауссовский фильтр имеет вид
Рассмотрим пример фильтрации изображения (рис. 1) в частотной области (рис. 17 - 22). Заметим, что частотная фильтрация изображения может иметь смысл как сглаживания ($\textit$), так и выделения контуров и мелкоразмерных объектов ($\textit$).
Как видно из рис. 17, 19, по мере нарастания "мощности" фильтрации в низкочастотной составляющей изображения все сильнее проявляется эффект "кажущейся расфокусировки" или $\it$ изображения. В то же время в высокочастотную составляющую, где в начале наблюдаются лишь контура объектов, постепенно переходит большая часть информационного содержания изображения (рис. 18, 20 - 22).
Низкочастотная фильтрация с параметрами $(10,10)$ Высокочастотная фильтрация с параметрами $(10,10)$ Низкочастотная фильтрация с параметрами $(50,50)$ Высокочастотная филь\-трация с параметрами $(50,50)$ Высокочастотная фильтрация с параметрами $(100,100)$ Высокочастотная фильтрация с параметрами $(200,200)$ НЧ-фильтрация зашумленного изображения $(10,10)$ ВЧ-фильтрация зашумленного изображения $(10,10)$ НЧ-фильтрация зашумленного изображения $(50,50)$ ВЧ-фильтрация зашумленного изображения $(50,50)$
Рассмотрим теперь поведение высокочастотных и низкочастотных фильтров (рис. 23 - 28) в присутствии аддитивного гауссовского шума на изображении (рис. 7).
Как видно из рис. 23, 25, свойства низкочастотных фильтров по подавлению аддитивной случайной помехи аналогичны свойствам ранее рассмотренных линейных фильтров - при достаточной мощности фильтра помехи подавляются, однако платой за это является сильное размытие контуров и "расфокусировка" всего изображения. Высокочастотная составляющая зашумленного изображения перестает быть информативной, так как помимо контурной и объектовой информации там теперь также полностью присутствует и шумовая компонента (рис. 27, 28).
ВЧ-фильтрация зашумленного изображения $(100,100)$ ВЧ-фильтрация зашумленного изображения $(200,200)$
Применение частотных методов наиболее целесообразно в случае, когда известны статистическая модель шумового процесса или/и оптическая передаточная функция канала передачи изображения. Учесть такие априорные данные удобно, выбрав в качестве восстанавливающего фильтра обобщенный управляемый (параметрами $\sigma$ и $\mu$) фильтр следующего вида:
где $0 < \sigma < 1$, $0 < \mu < 1$ - назначаемые параметры фильтра, $P(w_$, $w_)$ - передаточная функция системы, $Q(w_$, $w_)$ - стабилизатор фильтра, согласованный с энергетическим спектром фона. Выбор параметров $\sigma = 1$, $\mu = 0$ приводит к чисто инверсной фильтрации, $\sigma =\mu = 1$ к \it, что позволяет получить изображение, близкое к истинному в смысле минимума СКО при условии, что спектры плотности мощности изображения и его шумовой компоненты априорно известны. Для дальнейшего улучшения эффекта сглаживания в алгоритм линейной (винеровской) фильтрации вводят адаптацию, основанную на оценке локальных статистик: математического ожидания $M(P)$ и дисперсии $\sigma (P)$. Этот алгоритм эффективно фильтрует засоренные однородные поверхности (области) фона. Однако при попадании в скользящее окно обработки неоднородных участков фона импульсная характеристика фильтра сужается ввиду резкого изменения локальных статистик, и эти неоднородности (контуры, пятна) передаются практически без расфокусировки, свойственной неадаптивным методам линейной фильтрации.
К достоинствам методов линейной фильтрации следует отнести их ясный физический смысл и простоту анализа результатов. Однако при резком ухудшении соотношения сигнал/шум, при возможных вариантах площадного зашумления и наличии высокоамплитудного импульсного шума линейные методы предварительной обработки могут оказаться недостаточными. В этой ситуации значительно более мощными оказываются нелинейные методы.