Размещено на http://www.allbest.ru/

Министерство образования и науки Российской Федерации

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

"Башкирский государственный университет" (БашГУ)

Факультет математики и информационных технологий

Кафедра информационных технологий

Специальность 010501 - Прикладная математика и информатика

Дипломная работа

Цифровая обработка сигналов

Научный руководитель:

д. ф.-м. н., доцент

В.А. Делев

УФА - 2013



Оглавление

Введение

1. Линейные системы

1.1 Дискретные и непрерывные сигналы

1.1.1 Теорема Котельникова

1.1.2 Наложение спектров (алиасинг)

1.2 Свертка

2. Дискретное преобразование Фурье

2.1 Преобразования Фурье

2.2 Комплексное ДПФ

2.3 Двумерное ДПФ

3. Восстановление расфокусированных и смазанных изображений

3.1 Модель процесса искажения

3.1.1 Искажающие функции

3.2 Mодель шума

3.3 Теорема о свертке

3.4 Инверсная фильтрация

3.5 Существующие подходы для деконволюции

3.6 Сравнение изображений на практике

Заключение

Список литературы



Введение

сигнал преобразование фурье искажение

Цифровамя обрабомтка сигнамлов (ЦОС, DSP -- англ. digital signal processing) -- преобразование сигналов, представленных в цифровой форме.

Обработка сигналов во временной области широко используется в современной электронной осциллографии и в цифровых осциллографах. Для представления сигналов в частотной области используются цифровые анализаторы спектра. Для изучения математических аспектов обработки сигналов используются пакеты расширения (чаще всего под именем Signal Processing) систем компьютерной математики MATLAB, Mathcad, Mathematica, Maple и др.

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

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

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

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

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



1. Линейные системы

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

Будем пока для определенности рассматривать одномерные сигналы, зависящие от времени, и обозначать их x(t). Почти весь материал допускает обобщение и на многомерный случай.

Система - это некоторое преобразование сигнала. Система переводит входной сигнал x(t) в выходной сигнал y(t). Будем это обозначать так:

x(t)>y(t).

Обычно все рассматриваемые системы инвариантны к сдвигу, т.е. если x(t)>y(t), то x(t+T)>y(t+T). Это означает, что форма выходного сигнала зависит только от входного сигнала, а не зависит от времени начала подачи входного сигнала. Далее мы будем рассматривать только такие системы.

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

Линейная система - это система, в которой выполняется следующее свойство линейности: если x₁(t)>y₁(t) и x₁(t)>y₁(t), то б•x₁(t)+в•x₂(t)>б•y₁(t)+в•y₂(t). Здесь операции над сигналами следует понимать как операции над функциями от аргумента t.

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

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

Свойства линейных систем:

1. Постоянный (константный) сигнал переводится любой линейной системой в постоянный сигнал.

2. При прохождении через линейную систему синусоида остается синусоидой. Могут измениться лишь ее амплитуда и фаза (сдвиг во времени).

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

Что означает "прохождение синусоиды через линейную систему"? Это значит, что синусоида подается на вход системы бесконечно долго, т.е. от t =?? до t=+?. Если же синусоиду начали подавать лишь в некоторый конкретный момент времени (а до этого подавалось что-то другое, например, - 0), то после начала подачи синусоиды на вход мы можем получить синусоиду на выходе не сразу. Выходной сигнал постепенно начнет приобретать синусоидальную форму. Скорость "стремления к верной синусоиде" на выходе зависит от конкретной линейной системы.

1.1 Дискретные и непрерывные сигналы

1.1.1 Теорема Котельникова

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

Процесс замера величины сигнала через равные промежутки времени называется равномерной (по времени) дискретизацией. Многие устройства для ввода данных в компьютер осуществляют дискретизацию. Например, звуковая карта дискретизирует сигнал с микрофона, сканер дискретизирует сигнал, поступающий с фотоэлемента. В результате дискретизации непрерывный (аналоговый) сигнал переводится в последовательность чисел. Устройство, выполняющее этот процесс, называется аналогово-цифровым преобразователем (АЦП, analogue-to-digital converter, ADC). Частота, с которой АЦП производит замеры аналогового сигнала и выдает его цифровые значения, называется частотой дискретизации.

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

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

Говорят, что сигнал имеет ограниченный спектр, если выше определенной частоты все коэффициенты спектра равны нулю. В этом случае говорят, что спектр сигнала лежит ниже частоты F (ограничен частотой F), где F - частота, выше которой спектр равен нулю.

Теорема Котельникова (Найквиста, Шеннона): если сигнал таков, что его спектр ограничен частотой F, то после дискретизации сигнала с частотой не менее 2F можно восстановить исходный непрерывный сигнал по полученному цифровому сигналу абсолютно точно. Для этого нужно проинтерполировать цифровой сигнал "между отсчетами" специального вида функциями (sinc-функциями).

На практике эта теорема имеет огромное значение. Например, известно, что большинство звуковых сигналов можно с некоторой степенью точности считать сигналами с ограниченным спектром. Их спектр, в основном, лежит ниже 20 кГц. Это значит, что при дискретизации с частотой не менее 40 кГц мы можем потом более-менее точно восстановить исходный аналоговый звуковой сигнал по его цифровым отсчетам. Абсолютной точности достичь не удастся, так как в природе не бывает сигналов с идеально ограниченным спектром.

Устройство, которое интерполирует дискретный сигнал до непрерывного, называется цифро-аналоговым преобразователем (ЦАП, digital-to-analogue converter, DAC). Эти устройства применяются, например, в проигрывателях компакт-дисков для восстановления звука по цифровому звуковому сигналу, записанному на компакт-диск. Частота дискретизации звукового сигнала при записи на компакт-диск составляет 44100 Гц. Таким образом, говорят, что ЦАП на CD-плеере работает на частоте 44100 Гц.

1.1.2 Наложение спектров (алиасинг)

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

Предположим, что мы попытались оцифровать музыку, спектр которой ограничен частотой 20 кГц, но при записи какой-то электроприбор (например, дисплей) сгенерировал сильную помеху с ультразвуковой частотой 39 кГц, которая проникла в аналоговый звуковой сигнал. Мы производим оцифровку с частотой 44.1 кГц. При этом мы предполагаем, что звук, лежащий ниже частоты будет записан правильно (по теореме Котельникова). Но так как помеха лежит выше частоты 22.05 кГц, то возникнет алиасинг, и помеха "отразится" в нижнюю часть спектра, на частоту около 5 кГц. Если мы теперь попробуем пропустить полученный цифровой сигнал через ЦАП и прослушать результат, то мы услышим на фоне музыки помеху на частоте 5кГц. Помеха "переместилась" из неслышимой ультразвуковой области в слышимую область.

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

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

Существуют устройства, называемые фильтрами, которые позволяют изменять спектр сигнала. Например, фильтры низких частот (НЧ-фильтры, low-pass filters) пропускают без изменения все частоты ниже заданной, и удаляют из сиг-нала все частоты выше заданной. Эта граничная частота называется частотой среза (cutoff frequency) фильтра. Одно из важных применений НЧ-фильтров заключается в искусственном ограничении спектра сигнала перед оцифровкой. В этом случае фильтры называются анти-алиасинговыми, т.к. они предотвращают возникновение алиасинга при оцифровке сигнала. Частота среза анти-алиасинговых фильтров устанавливается равной половине частоты дискретизации.

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

В реальные АЦП почти всегда встраивается анти-алиасинговый фильтр. Обычно эффект от искусственного ограничения спектра вполне приемлем, в то время как алиасинг - недопустим. Однако не всегда искусственное ограничение спектра так благотворно влияет на записываемый сигнал. Например, при оцифровке музыки на низкой частоте дискретизации 11 кГц приходится отфильтровывать из спектра музыки все частоты выше 5.5 кГц. В результате этого музыка теряет в качестве (хотя обычно такие потери лучше, чем алиасинг). При оцифровке изображений необходимо аккуратно проектировать анти-алиасинговый фильтр, чтобы изменение спектра изображения не повлекло видимых артефактов (таких как пульсации вблизи резких границ).



1.2 Свертка

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

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

(формула свертки).

Рассмотренная операция получения результирующего сигнала по исходному называется сверткой (convolution). Итак, любая линейная система осуществляет свертку входного сигнала со своей импульсной характеристикой. Это записывается так:

y[n]=x[n]*h[n].

Функция h[n] называется ядром свертки (kernel) или импульсной характеристикой линейной системы.

Обычно все сигналы, обрабатываемые на компьютере, имеют конечную продолжительность (т.е. отличны от нуля лишь на конечном отрезке). Рассмотрим, что происходит с сигналом конечной продолжительности, когда его сворачивают с конечным ядром свертки. Пусть сигнал x[n] отличен от нуля только на отрезке от 0 до N-1 включительно ("имеет длину N"). Пусть ядро свертки h[n] отлично от нуля на отрезке от -m₁ до m₂ включительно, состоящем из M точек (M=m₁+m₂+1). Тогда при подстановке этих сигналов в формулу свертки, мы получим сигнал y[n], который отличен от нуля на отрезке от -m₁ до N-1+m₂ включительно. Таким образом, длина результирующего сигнала равна N+M-1, т.е. сумме длин исходного сигнала и ядра свертки минус один.

Итак, операция свертки расширяет сигнал на M-1 точку, где M - длина ядра свертки.



2. Дискретное преобразование Фурье

2.1 Преобразования Фурье

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

Преобразование Фурье (Fourier transform) - это разложение функций на синусоиды (далее косинусные функции мы тоже называем синусоидами, т.к. они отличаются от "настоящих" синусоид только фазой). Существует несколько видов преобразования Фурье.

1. Непериодический непрерывный сигнал можно разложить в интеграл Фурье.

2. Периодический непрерывный сигнал можно разложить в бесконечный ряд Фурье.

3. Непериодический дискретный сигнал можно разложить в интеграл Фурье.

4. Периодический дискретный сигнал можно разложить в конечный ряд Фурье.

Компьютер способен работать только с ограниченным объемом данных, следовательно, реально он способен вычислять только последний вид преобразования Фурье. Рассмотрим его подробнее.



2.2 Комплексное ДПФ

До сих пор мы рассматривали ДПФ от действительных сигналов. Обобщим теперь ДПФ на случай комплексных сигналов. Пусть x[n], n=0,…,N-1 - исходный комплексный сигнал, состоящий из N комплексных чисел. Обозначим X[k], k=0,…N-1 - его комплексный спектр, также состоящий из N комплексных чисел. Тогда справедливы следующие формулы прямого и обратного преобразований Фурье:

,

.

Если по этим формулам разложить в спектр действительный сигнал, то первые N/2+1 комплексных коэффициентов спектра будут совпадать со спектром "обычного" действительного ДПФ, представленным в "комплексном" виде, а остальные коэффициенты будут их симметричным отражением относительно половины частоты дискретизации. Для косинусных коэффициентов отражение четное, а для синусных - нечетное.

2.3 Двумерное ДПФ

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

,

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

Здесь N₁xN₂ - размер исходного сигнала, он же - размер спектра. k₁ и k₂ - это номера базисных функций (номера коэффициентов двумерного ДПФ, при которых эти функции находятся). Поскольку размер спектра равен размеру исходного сигнала, то k₁ = 0,…,N₁-1; k₂ = 0,…,N₂-1.

n₁ и n₂ - переменные-аргументы базисных функций. Поскольку область определения базисных функций совпадает с областью определения сигнала, то n₁ = 0,…,N₁-1; n₂ = 0,…,N₂-1.

Двумерное ДПФ (в комплексной форме) определяется следующими формулами (здесь x[n₁,n₂] - исходный сигнал, а X[k₁,k₂] - его спектр):

,

.

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

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

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

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



3. Восстановление расфокусированных и смазанных изображений

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

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

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

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



x₁ | x₂ | x₃ | x₄… - Исходное изображение.

После искажения значение каждого пикселя суммируется со значением левого, т.е. x'_i = x_i + x_i-1. Надо еще поделить на 2, но опустим этот шаг для простоты. В результате имеем размытое изображение со значениями пикселей:

x₁ + x₀ | x₂ + x₁ | x₃ + x₂ | x₄ + x₃… - Размытое изображение.

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

x₁ + x₀ | x₂ -- x₀ | x₃ + x₀ | x₄ -- x₀… - Восстановленное изображение.

В итоге, вместо размытого изображения получили исходное, к пикселям которого добавлена неизвестная константа x₀ с чередующимся знаком. Это уже намного лучше - эту константу можно подобрать визуально, можно предположить, что она примерно равна значению x₁, можно автоматически подобрать с таким критерием, чтобы значения соседних пикселей "скакали" как можно меньше и т.д. Но все меняется, когда мы добавляем шум (который всегда есть в реальных изображениях). При описанной схеме на каждом шаге будет накапливаться вклад шума в общую составляющую, что в итоге может дать совершенно неприемлемый результат. Но, как мы убедились, восстановление вполне реально даже таким примитивным способом.



3.1 Модель процесса искажения

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

f(x, y) - исходное неискаженное изображение.

h(x, y) - искажающая функция.

n(x, y) - аддитивный шум.

g(x, y) - результат искажения, т.е. то, что мы наблюдаем в результате (смазанное или расфокусированное изображение).

Сформулируем модель процесса искажения следующим образом:

g(x, y) = h(x, y) * f(x, y) + n(x, y) . (1)

Задача восстановления искаженного изображения заключается в нахождении наилучшего приближения f'(x,y) исходного изображения. Рассмотрим каждую составляющую более подробно. С f(x, y) и g(x, y) все достаточно понятно. Разберем функцию h(x, y): в процессе искажения каждый пиксель исходного изображения превращается в пятно для случая расфокусировки и в отрезок для случая простого смазывания. Либо же можно сказать наоборот, что каждый пиксель искаженного изображения "собирается" из пикселей некоторой окрестности исходного изображения. Все это друг на друга накладывается и в результате мы получаем искаженное изображение. То, по какой функции размазывается или собирается один пиксель и называется функцией искажения.

Другие синонимы - PSF (Point spread function, т.е. функция распределения точки), ядро искажающего оператора, kernel и другие. Размерность этой функции меньше размерности самого изображения. К примеру, в начальном рассмотрении примера, размерность функции была 2, т.к. каждый пиксель складывался из двух.

3.1.1 Искажающие функции

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

PSF в случае размытия по Гауссу функцией fspecial('gaussian', 30, 8).

PSF в случае смазывания фунцией fspecial('motion', 40, 45).

Операция применения искажающей функции к другой функции (в данном случае к изображению) называется сверткой (convolution), т.е. некоторая область исходного изображения сворачивается в один пиксель искаженного изображения. Обозначается через оператор "*" (не путать с умножением). Математически для изображения f с размерами MxN и искажающей функции h c размерами mxn это записывается так:

g(x,y) = h(x,y)*f(x,y) = h(i,j)f(x+i,y+j). (2)

Где , . Операция, обратная свертке, называется деконволюцией (deconvolution), решение такой задачи нетривиально.

3.2 Mодель шума

Осталось рассмотреть последнее слагаемое, отвечающее за шум, n(x, y) в формуле (1). Причины шума в цифровых сенсорах могут быть разными, но основные - это тепловые колебания и темновые токи. На величину шума также влияет ряд факторов, таких как значение ISO, тип матрицы, размер пикселя, температура, электромагнитные наводки и пр. В большинстве случаев шум является Гауссовым (который задается двумя параметрами - средним и дисперсией), а также аддитивным, не коррелирует с изображением и не зависит от координат пикселя.

3.3 Теорема о свертке

Вернемся теперь к первоначальной постановке задачи восстановления - необходимо каким-то образом обратить свертку, при этом, не забывая про шум. Из формулы (2) видно, что получить f(x, y) из g(x, y) не просто - если решать, то получится огромная система уравнений. Но можно воспользоваться преобразованием Фурье. Теорема о свертке гласит, что операция свертки в пространственной области эквивалентна обычному умножению в частотной области (причем умножение поэлементное, а не матричное). Соответственно, операция обратная свертке эквивалентна делению в частотной области, т.е. это можно записать как:

h(x,y)*f(x,y) H(u,v)F(u,v). (3)

Где H(u, v), F(u, v) - Фурье-образы соответствующих функций. Значит процесс искажения из формулы (1) можно переписать в частотной области как:

G(u,v) = H(u,v)F(u,v)+N(u,v). (4).

3.4 Инверсная фильтрация

Поделим равенство (4) на H(u, v) и получаем следующую оценку исходного изображения:

. (5)

Это называется инверсной фильтрацией, но на практике практически никогда не работает. Посмотрим на последнее слагаемое в формуле (5) - если функция H(u,v) принимает значения близкие к нулю или нулевые, то вклад этого слагаемого будет доминирующим. Это практически всегда встречается в реальных примерах - для объяснения этого вспомним, как выглядит спектр после преобразования Фурье.



Берем исходное изображение, преобразуем его в полутоновое и используя Matlab, получаем спектр:



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



noise_var = 0.0000001 noise_var = 0.000005

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

3.5 Существующие подходы для деконволюции

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

. (6)

Этот результат был получен Винером в 1942 году. Функцией S здесь обозначаются энергетические спектры шума и исходного изображения. Поскольку, эти величины редко бывают известны, то дробь S_n / S_f заменяют на некоторую константу K, которую можно приблизительно охарактеризовать как соотношение сигнал-шум.

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



. (7)

Где y - параметр регуляризации, а P(u, v) - Фурье-преобразование оператора Лапласа (матрицы 3 * 3).

Еще один интересный подход предложили независимо Ричардосн [1972] и Люси [1974]. Метод так и называется "метод Люси-Ричардсона". Его отличительная особенность в том, что он является нелинейным, в отличие от первых трех - что может дать лучший результат. Вторая особенность - метод является итерационным, соответственно возникают трудности с критерием останова итераций. Основная идея состоит в использовании метода максимального правдоподобия, для которого предполагается, что изображение подчиняется распределению Пуассона. Формулы для вычисления без использования преобразования Фурье - все делается в пространственной области:

. (8)

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

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

3.6 Сравнение изображений на практике

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



Результаты:

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

Фильтр Люси-Ричардсона

Регуляризация по Тихонову



Слепая деконволюция

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

Запускаем скрипт:



И получаем следующий результат:

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

Теперь рассмотрим пример с реальным смазыванием:



Скрипт примерно тот же, только тип PSF поменяем на "motion":

Результат:



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



Заключение

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

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

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

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

Так же в работе вводятся основные понятия ЦОС: линейные системы, свертка, дискретное преобразование Фурье.



Список литературы

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

2. Глинченко А.С. Цифровая обработка сигналов. - Красноярск: Изд-во КГТУ, 2001 г. - 199 с.

3. Апальков И.В., Хрящев В.В. Удаление шума из изображений на основе нелинейных алгоритмов с использованием ранговой статистики. - Ярославский государственный университет, 2007. - 105 с.

4. Дьяконов В., Абраменкова И. MATLAB. Обработка сигналов и изображений. Специальный справочник. - СПб.: Питер, 2002, 608 с.

5. Сергиенко А.Б. Цифровая обработка сигналов. - СПб.: Питер, 2007 г. - 751с.

6. Оппенгейм А., Шафер Р. Цифровая обработка сигналов. - М.: Техносфера, 2007 г. - 856 с.

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

Размещено на Allbest.ru