Источник: http://cgm.computergraphics.ru/content/view/74

Проблема подавления шума на изображениях и видео и различные подходы к ее решению

Авторы: Калинкина Д., Ватолин Д.

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

Введение в проблему подавления шума

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

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

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

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

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

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

Зачем нужно подавлять шум?

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

Также шумоподавление играет важную роль при сжатии видеопоследовательностей и изображений [25, 26]. И в видео, и в изображениях сжатие основано на пространственной корреляции значений пикселей, а в случае видео - еще и временной корреляции (на сходстве последовательно идущих кадров между собой). Одной из основных проблем в алгоритмах сжатия является определение локальной зашумленности данной области изображения, поскольку при сжатии сильный шум может быть принят за детали изображения, и это может, во-первых, привести к увеличению сложности с точки зрения сжатия и, во-вторых, отрицательно повлиять на результирующее качество сжатого изображения. При сжатии видеоизображений наличие шума приводит, кроме того, к увеличению межкадровой разницы, понижая, таким образом, степень сжатия, и отрицательно влияет на точность работы алгоритма компенсации движения [27].

Откуда же берется шум? Источники шума могут быть различными:

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

Соответственно, шумы тоже бывают разных видов. Самые распространенные:

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

Менее распространенные виды шума - шум вида <соль и перец> (как правило, вызывается помехами в электросети) и вертикальные царапины (характерны для старых черно-белых кинолент, возникают при механическом повреждении эмульсии на пленке).

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

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

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

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

Еще одной сложностью является оценка качества подавления шума. Как правило, оно оценивается следующим образом: на незашумленное изображение накладывается искусственный шум, затем полученное изображение фильтруется алгоритмом шумоподавления и сравнивается с исходным с помощью какой-нибудь метрики. Чаще всего для этой цели используют метрику PSNR (peak signal-to-noise ratio - пиковое соотношение сигнал/шум), которая определяется формулой:
,
где xi, yi - i-е пиксели двух сравниваемых изображений, N - общее число пикселей на каждом изображении, а d(xi, yi) - разность между цветами соответствующих пикселей (для оттенков серого это просто разность значений пикселей, а для цветных изображений - евклидово расстояние между пикселями в трехмерном цветовом пространстве). Соответственно, чем ближе отфильтрованное изображение к оригинальному, тем больше значение PSNR, и тем выше считается качество работы алгоритма.

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

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

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

Обзор существующих методов

Можно выделить следующие базовые подходы к пространственному шумоподавлению:

  1. Линейное усреднение пикселей по соседям
  2. Медианная фильтрация
  3. Математическая морфология
  4. Гауссовское размытие
  5. Методы на основе вейвлет-преобразования
  6. Метод главных компонент
  7. Анизотропная диффузия
  8. Фильтры Винера

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

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

Рассмотрим каждый подход подробнее.

Линейное усреднение пикселей

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

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

Эти алгоритмы очень простые, но они не дают хорошего результата.

Интересная модификация этого метода была предложена Де Хааном [4-6]. Он предложил в качестве значения центрального пикселя также брать взвешенную сумму соседних пикселей, только соседей брать не подряд, а через один или два пикселя (см. рисунок 2.1.1б). Утверждается, что при таком подходе удается подавить низкочастотный шум, который заметнее на глаз, чем высокочастотный.

 
Рисунок 2.1.1а.
Обычный вид окна
 Рисунок 2.1.1б.
Вид окна в алгоритме Де Хаана

Аналогично можно применять этот метод во временной области, только усреднение будет производиться уже между соседними кадрами, и окно соответственно будет браться по времени (то есть каждый пиксель будет усредняться по пикселям, расположенным в той же позиции в соседних кадрах). В общем виде такую схему шумоподавления можно выразить следующей формулой:
,
где x - пиксель, а t - номер кадра. Веса могут браться в зависимости от близости значений пикселей и расстояния между кадрами. Также усреднение может проводиться рекурсивно:
,
где - значение, посчитанное для этого пиксела в предыдущем кадре. Для предотвращения возникновения ореолов вокруг движущихся объектов, о которых было сказано выше, во временные фильтры встраивают алгоритмы определения движения. При этом возможно два варианта: простое детектирование движения (пикселы в движущихся блоках просто остаются без изменения, и шум вдоль движущихся объектов не подавляется) или построение скомпенсированного предыдущего и/или следующего кадра (см. [27]) и смешивание текущего с ним. В последнем случае компенсация движения должна быть выполнено качественно, иначе будут артефакты на месте неправильно найденных блоков.

Простой метод, усредняющий пиксели по семи кадрам и учитывающий движение, описан в [7].

Медианная фильтрация

Медианная фильтрация - это стандартный способ подавления импульсного шума.

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

Медиану также можно определить формулой:
,
где W - множество пикселей, среди которых ищется медиана, а fi - значения яркостей этих пикселей.

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

Однако в чистом виде медианный фильтр размывает мелкие детали, величина которых меньше размера окна для поиска медианы, поэтому на практике практически не используется. Пример усовершенствованной медианной фильтрации можно найти в [1].

Математическая морфология

Шумоподавление можно также осуществлять с использованием двух основных морфологических операций (см. [2]): сужения (erosion) и расширения (dilation), а также их комбинаций - закрытия (closing) и раскрытия (opening). Раскрытие (сначала сужение, потом расширение) убирает выступы на границах объектов, а закрытие (сначала расширение, потом сужение) заполняет отверстия внутри и на границах.

В качестве примера можно предложить следующий алгоритм, описанный в [3]. Сначала по исходному изображению I вычисляется новое изображение I', равное полусумме открытия-закрытия и закрытия-открытия исходного изображения. Получим сглаженное изображение, не содержащее шума. Тогда изображение D, равное разнице I и I',будет содержать весь шум и все те детали исходного изображения, размер которых меньше размера структурного элемента, примененного при морфологических операциях. Предполагая, что амплитуда у шума меньше, чем у деталей, обнулим в D все значения, меньшие некоторого порога, и снова сложим с I'.

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

Гауссовское размытие

Гауссовское размытие - это свертка изображения с функцией
,
где параметр s задает степень размытия, а параметр A обеспечивает нормировку. Фактически, это то же усреднение, только пиксель смешивается с окружающими по определенному закону, заданному функцией Гаусса. Матричный фильтр, посчитанный по указанной формуле, называется гауссианом; чем больше его размер, тем сильнее размытие (при фиксированном s). Поскольку данный фильтр сепарабелен, то есть представим в виде
,
то свертку можно производить последовательно по строкам и по столбцам, что приводит к значительному ускорению работы метода при больших размерах фильтра.

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

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

Методы на основе вейвлет-преобразования

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

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

Схема одномерного вейвлет-преобразования

Рисунок 1. Схема одномерного вейвлет-преобразования

Здесь H1 (G1) - высокочастотный фильтр анализа (синтеза), а H2 (G2) - низкочастотный фильтр анализа (синтеза).

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

Схема двумерного вейвлет-преобразования

Рисунок 2. Схема двумерного вейвлет-преобразования

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

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

Схема трехмерного вейвлет-преобразования

Рисунок 3. Схема трехмерного вейвлет-преобразования

Подробнее о ДВП можно прочитать в [16, 22, 23].

Итак, общую схему шумоподавления с помощью ДВП можно описать следующим образом (см. рисунок 4):

Общая схема подавления шума на основе вейвлет-преобразования

Рисунок 4. Общая схема подавления шума на основе вейвлет-преобразования

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

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

Методы шумоподавления с использованием ДВП можно классифицировать по двум параметрам: по типу преобразования и способу подавления коэффициентов.

Рассмотрим вначале различные типы ДВП и аналогичных им преобразований:

1.    Обычное ДВП

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

Самый простой способ получения избыточности - не производить прореживание, то есть использовать следующую схему (см. рисунок 5):

Схема одномерного избыточного вейвлет-преобразования

Рисунок 5. Схема одномерного избыточного вейвлет-преобразования

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

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

2.    Комплексное ДВП (dual-tree complex DWT)

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

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

Какие преимущества у такого ДВП? Во-первых, избыточность стала намного меньше - теперь она составляет 3N2 при любом количестве уровней разложений. Во-вторых, такие вейвлеты имеют лучшую направленность и различают, в отличие от обычных вейвлетов, направления в 45о и 135о, что способствует более качественной обработке диагональных границ деталей на изображении (примеры тестовых изображений можно посмотреть на сайте [16]).

В трехмерном случае действительное ДВП будет состоять уже из четырех параллельных ДВП, число полученных в результате вейвлетов будет равно 28 и направлены они будут в 28-ми разных направлениях, а комплексное ДВП также будет состоять из двух параллельных действительных преобразований. Избыточность будет составлять уже 7N2 для любого числа уровней разложения (против +7N2 на каждом уровне разложения у обычного ДВП без прореживания). Кроме того, комплексные трехмерные вейвлеты разделяют движение по 28 различным направлениям (обычные ДВП - только по 7-ми).

Реализацию на МатЛабе и примеры НЧ/ВЧ-фильтров для комплексного ДВП можно найти в [16-18].

3.    Несепарабельное ДВП (non-separable DWT)

Обычное ДВП является сепарабельным. Это значит, что оно может вычисляться отдельно по строкам и по столбцам, причем и для строк, и для столбцов применяются одинаковые фильтры. Несепарабельные ДВП могут лучше адаптироваться к локальным свойствам изображения. Например, анализируя для каждого пикселя его соседей в некотором окне можно минимизировать полезную энергию изображения в ВЧ-полосе и лучше сохранить детали. Алгоритм подавления шума с использованием несепарабельного ДВП описан в [19].

4.    Направленные вейвлеты (directional wavelets)

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

5.    Упрощенные преобразования

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

Теперь перейдем к рассмотрению методов подавления коэффициентов.

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

1.    Резкое уменьшение по порогу (hard thresholding)

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

2.    Плавное уменьшение по порогу (soft thresholding)

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

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

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

1.    Адаптация по соседним коэффициентам

Этот подход основан на том, что степень подавления каждого вейвлет-коэффициента зависит не только от величины этого коэффициента, но и от взвешенной суммы модулей (при анализе вейвлет-коэффициентов имеет значение именно их амплитуда, а не знак) соседних коэффициентов, расположенных в некоторой пространственной окрестности вокруг этого коэффициента. Зависимость может быть самая различная. Можно, кроме того, учитывая взаимосвязь коэффициентов между уровнями, включить в эту сумму коэффициенты из следующего уровня разложения. Веса также можно задавать различными способами. Как вариант, можно задать бóльшие веса вдоль основного направления границ в полосе: либо взять одно из трех основных направлений (вертикальное, горизонтальное или диагональное) в зависимости от полосы [8], либо для большей точности самостоятельно оценить доминирующее направление и искать веса, например, по формуле:
,
где Yi,j - отдельный вейвлет-коэффициент рассматриваемой полосы, w - вектор, содержащий искомые веса, а ui,j - вектор, состоящий из соседних по отношению к Yi,j вейвлет-коэффициентов [9]. То есть ищутся такие веса, чтобы взвешенная сумма соседей максимально приближала каждый коэффициент в полосе.

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

Следующие два подхода ориентированы на подавление белого гауссовского шума (хотя на практике это методы хорошо работают и для других широкополосных шумов). Пусть исходное изображение содержит аддитивный белый гауссовский шум с нормальным распределением: нулевым математическим ожиданием и дисперсией σ2. Аддитивность означает отсутствие корреляции между сигналом и шумом: NSi = Si+Ni, где NS - обозреваемый зашумленный сигнал, S - исходный чистый сигнал, а N - шум. Тогда для полученных вейвлет-коэффициентов wi будет также верно:

wi = yi + zi, i=1,...,n, где

yi - значение вейвлет-коэффициентов без шума, а zi ~ N(0, σ2) - вклад шума. В случае, если дисперсия шума σ2 неизвестна, она обычно оценивается как медиана среди абсолютных значений вейвлет-коэффициентов ВЧ-полосы (HH), деленная на 0.6745.

2.    Метод, основанный на контексте (Context-based method) [9]

В этом методе шумоподавление осуществляется по формуле плавного уменьшения по порогу, где порог задается следующим образом:

Здесь σw - дисперсия зашумленного сигнала, которая вычисляется таким образом, чтобы обеспечить локальную адаптивность величины порога. Сначала все вейвлет-коэффициенты группируются по своему контексту, то есть для каждого коэффициента вычисляется среднее абсолютное значение коэффициентов в некоторой его окрестности (контекст), и коэффициенты с близкими средними значениями объединяются в одну группу. Иными словами, для каждого коэффициента с контекстом Z определяется некое подмножество коэффициентов, имеющих самые близкие к Z значения контекстов (как меньшие Z, так и большие Z). После этого значение σw вычисляется для каждого вейвлет-коэффициента отдельно как дисперсия определенного для этого коэффициента подмножества. Мощность этого подмножества важно выбрать оптимально, поскольку слишком мало коэффициентов не обеспечат объективной оценки, а слишком много нарушат локальность. Таким образом, мы как бы выделяем похожие области на изображении, которые пространственно могут располагаться далеко друг от друга, и по ним считаем дисперсию.

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

3.    Байесовский подход (Bayesian approach)

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

Гистограмма распределения вейвлет-коэффициентов для изображения без шума

Гистограмма распределения вейвлет-коэффициентов для изображения с шумом

Рисунок 6а. Гистограмма распределения вейвлет-коэффициентов для изображения без шума

Рисунок 6б. Гистограмма распределения вейвлет-коэффициентов для изображения с шумом

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

Дисперсия такого распределения вычисляется по формуле:
.

Для изображений параметр υ обычно Î [0, 1], и его можно просто положить равным 1.

Шум описывается функцией нормального гауссовского распределения N(0, σ2):

Тогда для зашумленных коэффициентов w имеем:

.

Поскольку можно вычислить (это дисперсия исходного зашумленного сигнала), выразим из этой формулы параметр s (положим υ = 1):

.

Таким образом, модель незашумленного исходного сигнала полностью восстановлена.

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

при υ = 1.

Здесь - неполная гамма-функция.

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

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

Тогда для вейвлет-коэффициентов зашумленного изображения функции плотности распределения и могут быть получены как результат свертки нормального распределения N(0, σ2) с функциями плотности распределения и соответственно. Теперь для каждого вейвлет-коэффициента можно посчитать следующую величину:
.

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


где
,
а M - среднее абсолютное значение вейвлет-коэффициентов в некотором окне вокруг вейвлет-коэффициента w.

Метод главных компонент

Метод главных компонент (МГК) позволяет выделить структуру в многомерном массиве данных и применяется в основном для распознавания или для сжатия изображений. В области шумоподавления этот подход является довольно новым и мало исследованным. Работает он лучше всего для изображений с белым гауссовским шумом.

Метод главных компонент заключается в нахождении таких базисных векторов исследуемого многомерного пространства, которые бы наилучшим образом отражали расположение некоторых исходных данных в этом пространстве (характеризовали эти данные). Кратко опишем алгоритм построения таких векторов (математические понятия и подробное описание алгоритма см. в [14]).

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

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

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

Далее, вокруг каждого блока рассматривается некоторое окно, из которого выбираются всевозможные блоки того же размера, что и центральный. Вектора, соответствующие этим блокам вкупе с центральным, образуют исходный набор данных, для которого применяется описанный выше алгоритм МГК. Строится матрица ковариаций, ищутся собственные вектора и собственные значения, и затем производится сортировка векторов по убыванию собственных значений. Шумоподавление осуществляется путем уменьшения длин собственных векторов и/или их обнуления в зависимости от величины соответствующих собственных значений. Пример адаптивного алгоритма подавления величин собственных векторов можно посмотреть в [15].

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

Итак, как же уменьшение длин собственных векторов влияет на подавление шума? Дело в том, то найденные таким способом вектора будут соответствовать основным направлениям узора (границ, деталей) на данной области изображения. Поскольку шум распределен по изображению хаотично, безо всякой системы, ему будут соответствовать собственные вектора с маленькими собственными значениями. Подавляя их, мы будем уменьшать амплитуду значений пикселей в этом направлении, а следовательно, подавлять и шум. Именно поэтому метод главных компонент хорошо справляется лишь с белым шумом: у шума не должно быть никакой структуры, иначе МГК примет его за узор и не подавит.

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

Анизотропная диффузия

Основная идея данного подхода заключается в следующем. Яркость каждого пикселя интерпретируется как значение температуры в данной точке изображения, таким образом, все изображение представляется в виде карты температур. Шумоподавление производится путем выравнивания температур (фактически, интенсивностей пикселей) с помощью моделирования процесса теплопереноса:
,
где c - коэффициент теплопроводности, а - разность яркостей (<температур>) двух пикселей.

Коэффициент теплопроводности c зависит от градиента интенсивности изображения:
.

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

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

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

Фильтры Винера основаны на преобразовании Фурье [24] и применяются, в основном, для подавления шума в аудио, поэтому автором они подробно не исследовались. Посмотреть соответствующий алгоритм можно, например, в [28].

Литература

  1. Szymon Graboskwi, Wojciech Bienieck "A two-pass median-like filter for impulse noise removal in multi-channel images". KOSYR 2003, str. 195-200
  2. Презентация: <Введение в математическую морфологию> (http://graphics.cs.msu.su/courses/cg01b/mat_morph.zip)
  3. Richard Alan Peters II «A New Algorithm for Image Noise Reduction using Mathematical Morphology». IEEE Transactions on Image Processing, Volume 4, Number 3, pp. 554-568, May 1995
  4. G. de Haan, T. G. Kwaaitaal-Spassova and O.A. Ojo «Automatic 2-D and 3-D noise filtering for high-quality television receivers». International Workshop on HDTV'94,4-B-2,1994-10
  5. R. J. Schutten, G. de Haan and A. H. M. van Roermund «Noise filtering for television receivers with reduced memory». Proc. of the Int. Workshop on HDTV and the Evolution of Television, Taipei, Taiwan, Nov. 1995, pp. 6A15--6A22
  6. G. de Haan, T.G. Kwaaitaal-Spassova, M.M. Larragy, O.A. Ojo, and R.J. Schutten «Television noise reduction IC». Digest of the ICCE'96, June 1996, Chicago, pp. 184-185
  7. Håkan Norell, Bengt Oelmann and Youshi Xu «Spatio-Temporal Noise Reduction ASIC for Real-Time Video Processing». In Proc. IEEE Nordic Signal Processing Symposium, Kolmården, Sweden, 13-15 June, 2000
  8. Detlev Marpe, Hans L. Cycon, Gunther Zander, Kai-Uwe Barthel «Context-based Denoising of Images Using Iterative Wavelet Thresholding». Proc. SPIE, Vol. 4671, pp. 907-914, Jan. 2002
  9. S.G. Chang, B. Yu, and M. Vetterli, «Spatially adaptive wavelet thresholding with context modeling for image denoising». IEEE Trans. Image Proc., Sept. 2000.
  10. Aleksandra Pizurica, Vladimir Zlokolica and Wilfried Philips «Noise Reduction in Video Sequences Using Wavelet-Domain and Temporal Filtering». SPIE Conference on Wavelet Applications in Industrial Processing 27-31 October 2003.
  11. Aleksandra Pizurica and Wilfried Philips «Adaptive Probabilistic Wavelet Shrinkage for Image Denoising». IEEE Transactions on Image Processing, July 11, 2003.
  12. Aleksandra Pizurica, Wilfried Philips, Ignace Lemahieu, and Marc Acheroy «A Versatile Wavelet Domain Noise Filtration Technique for Medical Imaging». IEEE Trans Med Imaging. 2003 Mar;22(3):323-31
  13. Aleksandra Pizurica, Vladimir Zlokolica and Wilfried Philips «Combined Wavelet Domain and Temporal Video Denoising». IEEE Conference on Advanced Video and Signal Based Surveillance (AVSS'03), July 2003.
  14. Lindsay I Smith «A tutorial on Principal Components Analysis». 2002 (http://www.cs.otago.ac.nz/cosc453/student_tutorials/principal_components.pdf)
  15. D. Darian Muresan, Thomas W. Parks «Adaptive principal components and image denoising». IEEE ICIP 2003
  16. Интернет-ресурс: http://taco.poly.edu/WaveletSoftware/.
  17. Ivan W. Selesnick and Ke Yong Li «Video denoising using 2D and 3D dual-tree complex wavelet transforms». Appl Signal Image Proc. X (Proc. SPIE 5207), Aug 2003
  18. Julia Neumann Gabriele Steidl «Dual Tree Complex Wavelet Transform in the Frequency Domain and an Application to Signal Classification». Technical Report TR-03-013, Department for Mathematics and Computer Science, University of Mannheim, 2003
  19. Miroslav Vrankić, Damir Seršić «Image Denoising Based on Adaptive Quincunx Wavelets». In Proc. of the 2004 IEEE 6th Workshop on Multimedia Signal Processing, (MMSP 2004), Siena, Italy, September 29 - October 01, 2004, pp. 251-254
  20. Vladan Velisavljevic, Baltasar Beferull-Lozano1, Martin Vetterli, Pier Luigi Dragotti «Discrete Multi-Directional Wavelet Bases». Proc. IEEE International Conference of Image Processing (ICIP2003)
  21. Felix C. A. Fernandes, Rutger L. van Spaendonck, C. Sidney Burrus «A Directional, Shift-Insensitive, Low-Redundancy, Wavelet Transform». IEEE Proceedings on International Conference on Image Processing (ICIP2001)
  22. Daubechies «Ten Lectures on Wavelets», Philadelphia: SIAM, 1992.
  23. S. Mallat «A wavelet tour of signal processing», Academic Press, London, 1998.
  24. Интернет-ресурс: http://astronomy.swin.edu.au/~pbourke/analysis/.
  25. Д. Ватолин, А. Ратушняк, М. Смирнов, В. Юкин «Методы сжатия данных. Устройство архиваторов, сжатие изображений и видео». Диалог-МИФИ, 2002
  26. Generic coding of moving pictures and associated audio information: Video. ISO/IEC 13818-2:2000
  27. G. de Haan «Progress in motion estimation for video format conversion». IEEE Transactions on Consumer Electronics Vol 46 No 3 Aug 2000 pp 449-450
  28. Javier Portilla, Vasily Strelay, Martin J. Wainwrightz, Eero P. Simoncelli «Adaptive Wiener Denoising Using A Gaussian Scale Mixture Model In The Wavelet Domain». Proceedings of the 8th IEEE Int'l Conf on Image Processing Thessaloniki, Greece. 7-10 October 2001

Вверх