Подавление шумов как задача диффузии

Hello, Хабр.

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

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

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

Представление изображения


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

image

Таким образом каждому пикселю изображения соответствует элемент матрицы. Значения в матрице — числа типа float из промежутка [0,255].

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



Изображение как функция


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

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

image

image

Тут k==1 так как разница между соседними пикселями = 1 пиксель.

Таким образом, производные будут тоже матрицами такого же размера как и исходное изображение:




Имея производные по обеим координатам логично посчитать градиент:
image

И его величину:
image

Величину градиента можно вывести в виде изображения:


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

Связь с процессом диффузии


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

image

Где D — Тензор диффузии, который задает направление и величину диффузии.

Запишем так же закон сохранения массы, который говорит о том, что разница входящего и исходящего потоков равна изменению массы в точке:

image

Если поток не ограничивать, то есть D==1 — масса будет распространяться во всех направлениях равномерно, то есть мы получим эквивалент Гауссовского размытия, а значит потеряем важные детали изображения. Нас интересует метод, с помощью которого поток можно ограничивать.
Одним из самых простых методов является ограничение величины потока, не ограничивая его направления. Про него и будем говорить далее.

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

image

Тогда уравнение сохранения массы принимает вид:

image

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



Таким образом мы получили итерационную схему, которая стабильна при значениях параметра tau<=0.25. Она довольна медленна, так как метод базовый и более успешные алгоритмы базируются на нем.

Результат


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

UPD: Приношу искренние извинения, изображения были перепутаны. Заменил изображение, которое было раньше на то, что получается в результате обработки изображения с таким количеством шума как представлено в статье. Еще раз приношу извинения.



Для сравнения — оригинал и версия с шумом (еще раз):



Литература:
  • P. Perona, J. Malik: Scale-space and edge detection using anisotropic diffusion. IEEE
    Transactions on Pattern Analysis and Machine Intelligence 12(7):629-639, 1990.
  • J. Weickert: Anisotropic Diffusion in Image Processing. Teubner, Stuttgart, 1998. The book is
    not available anymore, but the PhD thesis of Joachim Weickert, which was the basis of his
    book, is available on his website.
  • T. Brox: Lectures on Computer Vision. Freiburg im Breisgau, 2011


P.S. При обнаружении неточностей — буду рад услышать критику и пожелания.
Поделиться публикацией
Комментарии 22
    +3
    Красивая математика, внушительный список литературы, но написано суховато (IMHO). Чувствуется преподавательский стиль учебных заведений;)

    Думается, не помешало бы представить предложенные алгоритмы в псевдо-коде, а еще лучше на чём-нибудь удобочитаемом, Python например. Тем самым вдохнув в тему жизнь…
      0
      присоеденяюсь к пожеланию.
      Sourcecode?
        0
        Спасибо за комментарии :)
        С кодом проблемы нет, надо его только причесать, так как все равно для себя писал маленький отдельный примерчик для создания картинок, которые использовались в статье. Единственное что — библиотека с разными матрицами и производными — под копирайтом.
        Отсюда вопрос: как удобнее читателям:
        — оставить в с++, переписав маленькие частички, чтобы не использовать местную библиотечку
        — переписать во всесильном питоне (должно быть не очень долго)
        — остановиться на псевдокоде, оставив реализацию некоторых мелочей людям? :)

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

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

          Опыта написания статей подобного рода у меня откровенно говоря нет, поэтому благодарен за указания недочетов. :)
        0
        Если уже приводите в сравнение размытие по Гауссу, было бы совсем неплохо увидеть отличия работы этих алгоритмов (особенно на разных картинках)
          0
          Ах да, конечно, даже картинку подготовил, но забыл вставить. Добавлю сегодня вечерком скорее всего, как с работы прийду.
          +1
          Забавно. Особенно удивила степень восстановления маленькой надписи справа снизу. Как я понимаю в существующих коммерческих шумодавах вместо использования диффузии принимается пороговое решение о том, размывать или нет?
          Ради любопытства прогнал картинку через стоящий на компе Neat Image:
          image
            +1
            Если честно, я не очень знаю как действуют коммерческие шумодавы, но точно не так, потому что этот метод реально очень, ОЧЕНЬ медленный. На сколько я знаю, они используют взвешивание по окрестности каждого пикселя и таким образом сохраняют так же текстуры на изображении.
            По поводу надписи — меня это тоже поразило. Если честно, я сейчас не на 100% уверен, что результат получен именно с этой картинки, так как готовил эту лекцию когда-то довольно давно, и может быть брал менее ужасный шум на изображении (хотя в папке лежала только эта фотка). Теперь, после Вашего коментария, захотелось еще раз самому перепроверить. Так что скоро отрапортую.
            +1
            Прошу прощения, когда писал статью, перепутал изображения, чем мог ввести в заблуждение. Исправил на правильный и куда более правдоподобный вариант.
              +1
              Открыл что ли для себя, что физику не зря в институте изучают?
                0
                Да, примерно так оно и было :) Ну правда не ожидал, что ее можно применить в обработке изображений.
                  0
                  Абсолютно no offence, да и вы, видимо, дружелюбно восприняли.

                  Сообщество поборников высшего образования (в пику «самоучкам») с нетерпением ждет вашего появления в своих рядах.

                  Спасибо!
                    +1
                    Да, конечно, никаких обид, Вы что? :) Я если честно сейчас становлюсь фанатом отечественного образования со всеми вытекающими — решениями примеров повышенной сложности у доски без подготовки, шпоры, странные преподы, отношение большинства профессоров «я начальник — ты дурак». Мне все это правда нравится, потому что это закаляет. Сейчас я перебрался в Германию и я в сравнении с местыми студентами очень неплохо шарю, не смотря на то, что я уже ни не помню когда я сдавал экзамен без шпор. Не очень понимаю как, но наше образование хорошо функционирует :)
                0
                www.ipol.im/pub/demo/bcm_non_local_means_denoising/
                10 сек, и все готово. И это еще не самый крутой алгоритм денойзинга…

                image
                  0
                  Да-да, этот алгоритм знаю. Как раз про него где-то в ветке ответил Злодею Баалу, что они используют окрестность пикселя и считают по нему статистики. Алгоритм хорош так же тем, что сохраняет текстуры. Хотя конекретно картинка приведенная Вами выше навевает мне неверие в то, что они делают это честно. Я сам кодил несколько алгоритмов шумоподавления, в том числе и этот, но таких результатов я не видел. Хотя может это у меня руки не оттуда ростут.

                  Как бы там ни было, описанный в статье метод даже не претендует на практическое использование — он красив исключительно с описательной точки зрения.

                  Если есть желание — про алгоритм non local means так же могу статейку накидать. (Правда я не смотрел, может на Хабре уже такая и есть)
                    0
                    Как бы там ни было, описанный в статье метод даже не претендует на практическое использование
                    Но тем не менее, даже такой простой и прозрачный алгоритм дает хорошие результаты.
                      +1
                      Ну, если бы он не давал результатов — то им бы никто и не занимался, наверное. :)
                  0
                  Стало похоже на теплую ламповую фотографию!!!
                  Статье не хватает примера численной реализации данного метода, хотя бы в матлабе.
                    0
                    Да, по поводу этого уже кто-то писал. Код появится в статье на выходных скорее всего. В рабочие дни ни времени ни желания, к сожалению.

                  Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                  Самое читаемое