Фильтрование. Задача фильтрации изображений

Зашумление изображения. Модели шумов

Выше мы уже писали о том, что под задачей "фильтрации изображений" в широком смысле иногда понимают любые процедуры обработки изображений, при которых на вход процедуры подается (одно) растровое изображение, и на выходе также формируется растровое изображение. Такие процедуры типа (один растровый вход, один растровый выход) называют\linebreak $\it{фильтрами}$.

Однако чаще под "фильтрацией" в более узком смысле понимают так называемую $\textit{помеховую фильтрацию}$, или фильтрацию изображений от "шума". При этом неявно предполагается, что первоначально где-то существовало некое "исходное" идеально чистое (незашумленное) изображение, из которого затем путем $\it{зашумления}$ (определенного вида искажения), было получено то реальное изображение, которое мы наблюдаем. Задача помеховой фильтрации, таким образом, сводится к тому, чтобы путем некоторой обработки наблюдаемого реального изображения как можно лучше "очистить его от шума", то есть получить изображение, наиболее близкое по своим характеристикам к исходному "незашумленному" изображению.

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

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

Наиболее общей моделью независимого зашумления пикселов является $\it{шум замещения}$. Пусть дано исходное ("незашумленное") полутоновое изображение Im$$, каждый пиксел которого может принимать значения в диапазоне $$. Общая модель шума замещения предполагает, что после зашумления каждый пиксел изображения, имевший ранее значение яркости $i$, либо с некоторой известной вероятностью $p(i)$ это значение сохранит, либо данное значение яркости будет случайным образом замещено с вероятностью $q(i,j)$ некоторым другим значением яркости $j$ из того же конечного дискретного диапазона $$. Как видно, для описания такой общей модели случайного замещения нам потребуется задать таблицу $\it{переходных вероятностей}$ размера $I_{\textrm{max}}^{2}$, что составляет весьма значительное количество в случае обычного $8$-битового полутонового изображения (размер таблицы - $256\times 256$ элементов). Такое описание явно является некомпактным и поэтому редко используется на практике для полутоновых изображений. В то же время, для бинарных изображений, в которых $I_{\rm {max}} = 2$, такое описание является наиболее удобным, простым и естественным. Чуть ниже мы еще рассмотрим модель шума замещения на бинарных изображениях - так называемую модель шума "соль и перец".

Для полутоновых изображений, как правило, рассматривают другую, более частную модель зашумления - $\it{аддитивный шум}$, которая предполагает, что зашумленное изображение порождается по закону

$$ {\rm Im}^{\prime} = {\rm Im} + R(x,y), $$

где ${\rm Im}^{\prime} $ - пиксел зашумленного изображения, ${\rm Im} $ - пиксел исходного изображения, а $R(x,y)$ - случайная $\textit{аддитивная шумовая компонента}$. Кроме того, в большинстве приложений зависимость шума от координат пиксела считается несущественной. Наконец, исходя из известного в статистике $\it{закона больших чисел}$, закон распределения аддитивной шумовой компоненты предпочитают описывать удобным параметрическим семейством $\it{нормальных}$ ($\it{гауссовских}$) распределений с нулевым средним. Таким образом, $\textit{гауссовский аддитивный шум}$ описывается выражением

$$ {\rm Im}^{\prime} = {\rm Im} + N(0,\sigma), $$

где $N(a,\sigma)$ - нормальное распределение, $a$ - математическое ожидание нормально распределенного сигнала, $\sigma $ - средний квадрат отклонения (СКО) нормально распределенной величины. Именно такая модель зашумления чаще всего рассматривается в задачах фильтрации полутоновых изображений.

На рис. 2 - 8 показаны примеры искусственного зашумления исходного полутонового изображения лейкоцитов (рис. 1) специально сгенерированным аддитивным гауссовским шумом с различными значениями СКО. Как видно, чем больше параметр зашумления $\sigma $, тем более искаженным выглядит изображение. При больших значениях $\sigma $ (рис. 7, 8) даже человеческий глаз уже с трудом различает общие очертания крупноразмерных объектов изображения (в данном случае - лейкоцитов), более мелкие и менее контрастные объекты становятся практически неразличимы.

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

Оконная фильтрация изображений в пространственной области

Исходя из задачи восстановления исходных значений яркости незашумленного изображения, а также из того, что шумовая компонента каждого пиксела является заранее не известной случайной величиной, следует, что для решения данной задачи необходимо использовать ту или иную процедуру $\textit{статистического оценивания}$. Это может быть Байесовское оценивание, оценивание по методу наибольшего правдоподобия или любой другой метод, известный из курса математической статистики. Однако все эти методы требуют использовать для оценки искомой величины не одно единственное измерение (ведь оно также может быть зашумлено), а большую или меньшую $\textit{статистическую выборку, }$всегда включающую несколько отсчетов, характеризующих данную величину. В связи с этим и основная идея помеховой фильтрации изображений заключается в том, что для оценки исходного значения каждого пиксела изображения используется не только значение самого данного пиксела (как в ранее рассмотренных градационных преобразованиях), но и значения еще нескольких близких к нему пикселов, попадающих в так называемое "$\it{окно}$" или $\it{апертуру}$ фильтра. При этом "близость" пикселов к оцениваемому понимается в буквальном геометрическом смысле.

Наиболее простыми для вычислительной реализации являются традиционно используемые $\textit{прямоугольные окна}$ (апертуры) фильтрации, определяемые простым условием типа "$\textit{все пикселы данного окна отстоят от тестируемого центрального пиксела на более чем на WinX/2 по горизонатали и WinY/2 по вертикали}$", где WinX и WinY - горизонтальный и вертикальный размер окна фильтрации соответственно. Возможны и другие, более сложные способы формирования окон фильтрации - круглой, треугольной или любой другой произвольной формы.

Типовая процедура оконной фильтрации предполагает, что окно фильтрации последовательно движется по $\textit{входному изображению}$ (например, алгоритм может обходить изображение "в порядке чтения": сверху вниз по строкам, слева направо в каждой строке), при этом в каждом положении окна происходит анализ всех пикселов, принадлежащих в данный момент окну, и на основе такого анализа центральному пикселу окна на $\textit{выходном изображении}$ присваивается то или иное финальное значение. Сформированное таким образом выходное изображение также называется $\textit{результатом фильтрации}$.

Процедуры оконной фильтрации могут различаться:

1. размером и формой окна (апертуры);
2. типом собираемых в окне локальных статистик;
3. способом принятия решения на основе собранных статистик.

В любом случае, речь идет об использовании для оценивания значения центрального пиксела апертуры информации о значениях его соседей по изображению. В статистическом смысле это означает, что мы неявно опираемся на предположение о том, что на исходном незашумленном изображении значения яркостей всех этих соседних пикселов были одинаковыми или очень близкими, и наблюдаемые различия в их яркостях на зашумленном изображении определяются только присутствием шумовой компоненты, которую и необходимо исключить. Между тем, как мы уже видели, исследуя профили изображения, содержательное изображение вовсе не представляет собой одну сплошную "плоскость". В тех областях, которые визуально кажутся нам областями одинаковой или медленно меняющейся яркости, значения соседних пикселов действительно различаются незначительно. В то же время, на границах таких областей наблюдаются порой весьма резкие перепады яркости - разница значений составляет от десятков до сотен градаций интенсивности даже между непосредственно соседствующими пикселами. Таким образом, мы видим, что на границах однородных областей оконные фильтры не могут работать эффективно, напротив, здесь они с большой вероятностью будут ошибаться, что визуально приведет к эффекту $\textit{искажения формы контуров}$. Более того, если на исходном изображении присутствуют контрастные объекты (области), размер которых существенно меньше размера окна фильтрации, фильтр может просто "не заметить" такой объект, отфильтровать его как шум, что приведет к $\textit{исчезновению мелкоразмерных объектов}$ на результирующем выходном изображении.

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

Качество фильтрации

Таким образом, конструируя и исследуя оконные процедуры фильтрации изображений, мы всегда должны оценивать наблюдаемое $\textit{качество фильтрации}$ по двум следующим основным позициям:

1. способность фильтра удалять (отфильтровывать) с изображения шум;
2. способность фильтра сохранять на изображении мелкоразмерные детали и форму контуров.

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

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

Фильтрация с использованием функции imfilter

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

I = imread("coins.png"); h = ones(5,5) / 25 h = 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 I2 = imfilter(I,h); imshow(I); % Исходное изображение figure, imshow(I2); % Изображение после фильтрации

Типы данных

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

Таким образом, при использовании функции imfilter необходимо учитывать формат представления данных изображения. Рассмотрим пример, где в результате применения функции imfilter получены отрицательные значения. Исходные данные были представлены в формате double.

A = magic(5) A = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 h = [-1 0 1] h = -1 0 1 imfilter(A,h) ans = 24 -16 -16 14 -8 5 -16 9 9 -14 6 9 14 9 -20 12 9 9 -16 -21 18 14 -16 -16 -2

Отметим, что результат представлен также и отрицательными значениями. Если матрицу A представить не в формате double, а в формате uint8, то получим такой результат:

A = uint8(magic(5)); imfilter(A,h) ans = 24 0 0 14 0 5 0 9 9 0 6 9 14 9 0 12 9 9 0 0 18 14 0 0 0

После того, как исходные данные были представлены в формате uint8, то результирующие данные также представляются в формате uint8, а отрицательные значения будут обрезаны до 0. Таким образом, исходные данные перед применением функции imfilter можно преобразовать и в другие форматы, например, signed integer, single или double.

Опции корреляции и свертки

Функция imfilter может выполнять фильтрацию с использованием корреляции или свертки. По умолчанию используется корреляция, поскольку функции формирования фильтров формируют корреляционные ядра.

Однако, если необходимо выполнить фильтрацию с использованием свертки, нужно указать опцию "conv" в качестве необязательного аргумента в функции imfilter. Рассмотрим пример.

A = magic(5); h = [-1 0 1] imfilter(A,h) % фильтрация с использованием корреляции ans = 24 -16 -16 14 -8 5 -16 9 9 -14 6 9 14 9 -20 12 9 9 -16 -21 18 14 -16 -16 -2 imfilter(A,h,"conv") % фильтрация с использованием свертки ans = -24 16 16 -14 8 -5 16 -9 -9 14 -6 -9 -14 -9 20 -12 -9 -9 16 21 -18 -14 16 16 2

Опция дополнения края изображения

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

Выход значений маски фильтра за границу изображения

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

Заполнение недостающих пикселей нулевыми значениями

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

I = imread("eight.tif"); h = ones(5,5) / 25; I2 = imfilter(I,h); imshow(I), title("Исходное изображение"); figure, imshow(I2), title("Изображение после фильтрации");

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

Заполнение маски крайними пикселями изображения

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

I3 = imfilter(I,h,"replicate"); figure, imshow(I3);

Функция imfilter поддерживает также другие способы заполнения недостающих пикселей при обработке крайних элементов изображения. Для этого есть такие опции как "circular" и "symmetric". Для более детальной информации см. описание функции imfilter.

Многомерная фильтрация

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

Рассмотрим пример обработки цветного изображения, где реализована фильтрация каждого цветового слоя:

Система MATLAB располагает несколькими двумерными и многомерными функциями фильтрации. Функция filter2 выполняет двумерную линейную фильтрацию, conv2 выполняет свертку двумерных сигналов, а convn – свертку многомерных сигналов. При работе с каждой из этих функций необходимо, чтобы исходные данные были представлены в формате double. Результат также будет представлен в формате double. Эти функции всегда используют заполнение нулями недостающих пикселей при работе на краю изображения и не поддерживают других способов заполнения.

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

Использование предопределенных типов фильтров

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

I = imread("moon.tif"); h = fspecial("unsharp") h = -0.1667 -0.6667 -0.1667 -0.6667 4.3333 -0.6667 -0.1667 -0.6667 -0.1667 I2 = imfilter(I,h); imshow(I), title("Исходное изображение ") figure, imshow(I2), title("Обработанное изображение")

Частотные методы преобразований

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

Функция приложения ftrans2 предназначена для реализации частотных методов обработки изображений. Эта функция формирует маску линейного двумерного фильтра, используя метод преобразования частот для трансформации одномерного фильтра с конечной импульсной характеристикой. (Для более детальной информации см. Jae S. Lim, Two-Dimensional Signal and Image Processing, 1990.

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

B = remez(10,,); h = ftrans2(b); = freqz(b,1,64,"whole"); colormap(jet(64)) plot(w/pi-1,fftshift(abs(H))) figure, freqz2(h,)

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

Пример частотного метода

В частотных методах маска линейного фильтра формируется на основе желаемой амплитудно-частотной характеристики.

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

Hd = zeros(11,11); Hd(4:8,4:8) = 1; = freqspace(11,"meshgrid"); mesh(f1,f2,Hd), axis([-1 1 -1 1 0 1.2]), colormap(jet(64)) h = fsamp2(Hd); figure, freqz2(h,), axis([-1 1 -1 1 0 1.2])

Локальные (оконные) методы

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

Рассмотрим это на примере.

Hd = zeros(11,11); Hd(4:8,4:8) = 1; = freqspace(11,"meshgrid"); mesh(f1,f2,Hd), axis([-1 1 -1 1 0 1.2]), colormap(jet(64)) h = fwind1(Hd,hamming(11)); figure, freqz2(h,), axis([-1 1 -1 1 0 1.2])

Желательный двумерный частотный отклик (слева) и реальный двумерный частотный отклик (справа)

Создание желательной амплитудно-частотной характеристики

Функции проектирования фильтров fsamp2, fwind2 и fwind2 создают фильтры на основе матрицы значений желаемой амплитудно-частотной характеристики. Соответствующую желаемую амплитудно-частотную характеристику можно создавать с помощью функции freqspace.

Рассмотрим пример создания кругообразного низкочастотного фильтра со срезом в 0.5. Для этого используют следующий код:

Freqspace(25,"meshgrid"); Hd = zeros(25,25); d = sqrt(f1.^2 + f2.^2)

Идеальный кругообразный низкочастотный отклик

Вычисление частотного отклика фильтра

Функция freqz2 вычисляет частотный отклик двумерного фильтра. Если для фильтра с конечной импульсной характеристикой принять, что

H =;

то рассмотрим вычисление и визуализацию 64×64 точек частотного отклика h.

Freqz2(h)

Частотный отклик двумерного фильтра

Для получения частотного отклика H и вектора частот f1 и f2, используем следующий код

Freqz2(h);

Пример.

Рассмотрим пример размытия изображения. Реализуем эту процедуру с помощью двумерного фильтра низких частот (ФНЧ) с конечной импульсной характеристикой (КИХ), который сформирован из одномерного ФНЧ методом преобразования частот.

Некоторые виды преобразований в обработке изображений

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

1. Преобразование Фурье;
2. Дискретные косинусные преобразования
3. Преобразования Радона
4. Проекционные преобразования

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