Сглаживание изображений фильтром анизотропной диффузии Перона и Малика

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

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


Крайнее левое изображение — оригинальное, справа от оригинального — фильтрованные с различными параметрами.

Введение


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

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

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

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

Фильтр Перона и Малика


Впервые был описан в статье указанных математиков в 1990 году (ссылка на статью приведена в «источниках» в конце поста). Не будем вдаваться в глубокие теоретические детали, а хотя бы поверхностно рассмотрим как работает фильтр.

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

Теория


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

$$display$$I_{t}(x,y,t)=div(c(x,y,t)\nabla I(x,y,t)) \: \: \: \: \: \: \: \: (1)$$display$$


C начальным условием:

$$display$$I(x,y,0) = I_{0}(x,y)$$display$$


Где:

$inline$I(x,y,t)$inline$ — однопараметрическое семейство изображений; чем больше $inline$t$inline$, тем больше степень размытия исходного изображения;
$inline$I_{0}(x,y)$inline$ — исходное изображение;
$inline$I_{t}$inline$ — производная по $inline$t$inline$;
$inline$div$inline$ — оператор дивергенции;
$inline$\nabla$inline$ — градиент;

С точки зрения теории, изображение — это некоторая непрерывная функция двух переменных, а сама картинка (матрица пикселей) — дискретизация этой функции. Причем 0 соответствует нулевой яркости точки изображения, то есть обозначает черный цвет.

По своей сути фильтр анизотропной диффузии является модификацией фильтра Гаусса.

Если подставить вместо, пока еще не рассмотренной, функции $inline$c(x,y,t)$inline$ единицу, то есть $inline$с(x,y,t)\equiv 1$inline$, то получается уравнение изотропной диффузии:

$$display$$I_{t}(x,y,t)=div(\nabla I(x,y,t))=\frac{\partial^{2}I}{\partial x^{2}}+\frac{\partial^{2}I}{\partial y^{2}}$$display$$


Математики Коендеринк и Гумел показали, что такое семейство размытых изображений по параметру $inline$t$inline$, можно эквивалентно получить как решение уравнения свертки функции изображения с ядром Гаусса (это и есть фильтр Гаусса):

$$display$$I(x,y,t)= I_{0}(x,y)*G(x,y;t)$$display$$


Где:

$inline$*$inline$ — оператор свертки;
$inline$G(x,y;t)={\frac {1}{2\pi t^{2}}}e^{-{\frac {x^{2}+y^{2}}{2t^{2}}}}$inline$ — функция ядра Гаусса c нулевым математическим ожиданием и $inline$t$inline$ среднеквадратичным отклонением;

Очевидно, функция $inline$с(x,y,t)$inline$ играет роль некоторого «регулировщика» сглаживания.

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

$inline$c(x,y,t)=0$inline$ — на границах;
$inline$c(x,y,t)=1$inline$ — внутри областей, иначе говоря, внутри областей должно происходить обычное гауссово размытие;

Перона и Малик использовали градиент функции изображения $inline$\nabla I$inline$ как простую для расчета и достаточно точную аппроксимацию границ областей. Чем больше норма градиента, тем четче граница. Исходя из этого получаем:

$$display$$c(x,y,t)=g(||\nabla I(x,y,t)||)$$display$$


Где $inline$g(.)$inline$ — некоторая функция с областью значений на отрезке $inline$[0;1]$inline$.

Из-за нечеткого определения границ через норму градиента, требуют также, чтобы функция $inline$g$inline$ была монотонной убывающей.

Перона и Малик предложили два варианта функции $inline$g$inline$:

$$display$$g(x)=e^{-(\frac{x}{k})^{2}} \: \: \: \: \: \: \: \: (2)$$display$$


$$display$$g(x)=\frac{1}{1+(\frac{x}{k})^{2}} \: \: \: \: \: \: \: \: (3)$$display$$


Где $inline$k$inline$ — параметр, определяемый либо опытным путем, либо некоторым «измерителем» зашумленности.

Рассмотрим детальнее вторую функцию (формула 3). Для этого построим графики функции $inline$g$inline$ для нескольких разных $inline$k$inline$.



Видно, что чем больше $inline$k$inline$, тем больше значение функции $inline$g$inline$ для любой фиксированной точки.

Например, пусть в некоторой точке изображения, назовем ее $inline$(\tilde{x},\tilde{y})$inline$, мы знаем значение нормы градиента на уровне $inline$\tilde{t}$inline$: $inline$||\nabla I(\tilde{x},\tilde{y}, \tilde{t})||=4$inline$. Тогда $inline$g_{k=2}(||\nabla I(\tilde{x},\tilde{y}, \tilde{t})||)=g_{k=2}(4)=0.2$inline$, в то время как $inline$g_{k=16}(4)\approx 0.94$inline$. Получается в первом случае мы слабо «сгладили» значение в точке, так как по введенному ранее критерию она скорее всего лежит на границе, а во втором — значение функции $inline$g$inline$ практически единица, соответственно точка не считается граничной и будет сглажена обычным гауссовым размытием.

Таким образом, $inline$k$inline$ выступает в роли «барьера» для шумов, и с возрастанием $inline$k$inline$ возрастает «требование» на то, что точка будет считаться граничной.

Дискретизация


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

$$display$$\frac{\partial }{\partial x}I(x,y,t)=\frac{\partial }{\partial x}\left [ c(x,y,t)\frac{\partial }{\partial x}I(x,y,t) \right ]+\frac{\partial }{\partial y}\left [ c(x,y,t)\frac{\partial }{\partial y}I(x,y,t) \right ]$$display$$


Теперь запишем конечно-разностную явную схему:

$inline$ I(x,y,t+\Delta t)-I(x,y,t)=\\ \hspace{55mm}=\frac{1}{(\Delta x)^{2}}\left [ c(x+\frac{\Delta x}{2},y,t)\cdot (I(x+\Delta x,y,t)-I(x,y,t))-\\ \hspace{68mm}-c(x-\frac{\Delta x}{2},y,t)\cdot (I(x,y,t)-I(x-\Delta x,y,t))\right ]+\\ \hspace{55mm}+\frac{1}{(\Delta y)^{2}}\left [ c(x,y+\frac{\Delta y}{2},t)\cdot (I(x,y+\Delta y,t)-I(x,y,t))-\\ \hspace{68mm}-c(x,y-\frac{\Delta y}{2},t)\cdot (I(x,y,t)-I(x,y-\Delta y,t))\right ]$inline$

Записать условие устойчивости для получившейся явной схемы — нетривиальная задача из-за нелинейности уравнения. Но Перона и Малик определили, что при $inline$\Delta x=\Delta y=1$inline$ схема будет устойчива для всех $inline$0\leq \Delta t \leq \frac{1}{4}$inline$. Учитывая этот факт и то, что дискретным представлением функции изображения будет матрица значений пикселей, перепишем основную расчетную схему в матричном виде:

$$display$$I_{i,j}^{t+1}=I_{i,j}^{t}+\Delta t\left [ {C_{N}}_{i,j}^{t}\cdot \bigtriangledown_{N}I_{i,j}^{t}+{C_{S}}_{i,j}^{t}\cdot \bigtriangledown_{S}I_{i,j}^{t}+{C_{E}}_{i,j}^{t}\cdot \bigtriangledown_{E}I_{i,j}^{t}+{C_{W}}_{i,j}^{t}\cdot \bigtriangledown_{W}I_{i,j}^{t} \right ] $$display$$


Где:

$inline$\bigtriangledown_{N}I_{i,j}^{t}=I_{i-1,j}^{t}-I_{i,j}^{t}$inline$
$inline$\bigtriangledown_{S}I_{i,j}^{t}=I_{i+1,j}^{t}-I_{i,j}^{t}$inline$
$inline$\bigtriangledown_{E}I_{i,j}^{t}=I_{i,j+1}^{t}-I_{i,j}^{t}$inline$
$inline$\bigtriangledown_{W}I_{i,j}^{t}=I_{i,j-1}^{t}-I_{i,j}^{t}$inline$

$inline${C_{N}}_{i,j}^{t}=g(||{(\nabla I)}_{i-1/2,j}^{t}||)$inline$
$inline${C_{S}}_{i,j}^{t}=g(||{(\nabla I)}_{i+1/2,j}^{t}||)$inline$
$inline${C_{E}}_{i,j}^{t}=g(||{(\nabla I)}_{i,j+1/2}^{t}||)$inline$
$inline${C_{W}}_{i,j}^{t}=g(||{(\nabla I)}_{i,j-1/2}^{t}||)$inline$

Для расчета коэффициентов $inline$C$inline$ вначале необходимо вычислить норму градиента. Самый простой способ аппроксимировать норму градиента — это заменить ее на длину проекции градиента на соответствующую ось разностной сетки.

$inline${C_{N}}_{i,j}^{t}=g(|\bigtriangledown_{N}I_{i,j}^{t}|)$inline$
$inline${C_{S}}_{i,j}^{t}=g(|\bigtriangledown_{S}I_{i,j}^{t}|)$inline$
$inline${C_{E}}_{i,j}^{t}=g(|\bigtriangledown_{E}I_{i,j}^{t}|)$inline$
$inline${C_{W}}_{i,j}^{t}=g(|\bigtriangledown_{W}I_{i,j}^{t}|)$inline$

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

Окончательно расчетная формула будет выглядеть так:

$$display$$I_{i,j}^{t+1}=I_{i,j}^{t}+\Delta t\left [ g(|N|)\cdot N+g(|S|)\cdot S+g(|E|)\cdot E+g(|W|)\cdot W \right ] \: \: \: \: \: \: \: \: (4) $$display$$


Где:

$inline$N=\bigtriangledown_{N}I_{i,j}^{t}$inline$
$inline$S=\bigtriangledown_{S}I_{i,j}^{t}$inline$
$inline$E=\bigtriangledown_{E}I_{i,j}^{t}$inline$
$inline$W=\bigtriangledown_{W}I_{i,j}^{t}$inline$

Схематично расчетную схему можно изобразить как на рисунке ниже.


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


Чтобы фильтр менял значения яркости в ячейках в рамках максимума и минимума яркости исходного изображения, процесс, описанный основным уравнением, должен быть адиабатическим. Отсюда имеем граничное условие Дирихле вида: $inline$I(D)=0$inline$. То есть просто заполняем границы нулями.


Алгоритм и реализация на языке Fortran


Исходя из расчетной формулы в памяти всегда придется держать как минимум два двумерных массива, первый будет соответствовать оригинальному изображению, второй — размытому. Если провести аналогию с фильтром Гаусса, то для размытия изображения с радиусом $inline$t$inline$ в фильтре Перона и Малика нам нужно выполнить $inline$\frac{t}{\Delta t}$inline$ повторений полного пересчета изображения, каждый раз заменяя изображение с предыдущего слоя размытия на «оригинальное».
Последовательность действий будет примерно такая:

  1. Определяем ширину и высоту изображения (назовем w и h соответственно);
  2. Выделяем память под двумерный массив размерностью (w+2)*(h+2) (назовем I);
  3. Выделяем память под двумерный массив размерностью (w+2)*(h+2) (назовем BlurI);
  4. Заполняем массив нулями I=0;
  5. Считываем изображение и записываем данные пикселей в центр массива I (то есть оставляем слева, справа, сверху, снизу по одному нулевому столбцу или строке);
  6. Повторяем пока level<t (вначале level=0);

    1. Пересчитываем по формуле 4 все ячейки массива I, учитывая смещенную индексацию из-за нулевых полей. Считаем, что $inline$I_{i,j}^{t+1}$inline$ — это BlurI, а $inline$I_{i,j}^{t}$inline$ — I;
    2. Заменяем данные изображение в I на BlurI;
    3. Увеличиваем счетчик, level=level+deltaT, где deltaT — параметр, шаг по времени;

  7. Создаем и сохраняем изображение из данных массива BlurI. Это наше размытое изображение;
  8. Освобождаем память, выходим из программы;

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

Реализация на Fortran


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

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

Самым простым форматом, который к тому же понимают и многие графические редакторы, мне показался PGM P5. Это открытый формат хранения растровых изображений типа bitmap без сжатия в оттенках серого. У него простой ASCII заголовок, а само изображение есть последовательность однобайтных (для оттенков серого от 0 до 255) целых чисел без знака. Подробнее о формате вы можете прочитать непосредственно в самом стандарте, ссылка приведена в разделе источников.

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

Начнем программу с подключения модуля для работы с PGM и объявления всех необходимых переменных.

program blur
  use pgmio

  implicit none
  
  double precision, parameter :: t=10.0d0, deltaT=0.2d0, k=10.0d0
  character*(*), parameter :: input='dif_tomography.pgm', output='output.pgm'

  double precision, allocatable :: u(:,:), nu(:,:)
  double precision :: north, south, east, west, level
  integer :: w, h, offset, n, i, j

end program blur

Параметр t — это уровень размытия, на котором мы хотим прекратить работу алгоритма, deltaT — шаг по времени, k — параметр для пока еще не описанной функции g. Input и output — файл с исходным изображением и выходной файл соответственно.

Теперь считаем размеры входного изображения в переменные w и h.

call pgmsize(input, w, h, offset)

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

Выделим память для массива изображения и массива сглаженного изображения и считаем в первый данные из файла input.

allocate(u(0:w+1,0:h+1))
u=0
allocate(nu(w,h))
call pgmread(input, offset, w, h, u, 0, 0)

Фортран позволяет программисту самому определять индекс первого элемента в массиве, поэтому, учитывая, что вокруг изображения должны быть нули, заранее заполняем матрицу u нулями, задав ей размер от 0 до w+1 по первому измерения, и от 0 до h+1 по второму. Это позволит при дальнейшем пересчете использовать естественную индексацию без смещения.

Процедура pgmread считываем w*h байт, пропуская offset байт (занимаемых заголовком PGM) в массив u. Последние два параметра сообщают процедуре, что отсчет в матрице u начинается с нуля по каждому измерению.

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

Теперь перейдем с самому сглаживанию.

  level = 0 !счетчик уровня сглаживания
  do while (level<t)
    do j=1,h
        do i=1,w
            north = u(i-1,j)-u(i,j)
            south = u(i+1,j)-u(i,j)
            east = u(i,j+1)-u(i,j)
            west = u(i,j-1)-u(i,j)
            nu(i,j) = u(i,j) + deltaT*(g(north)*north+g(south)*south+g(east)*east+g(west)*west)
        end do
    end do
    u(1:w,1:h) = nu(1:w,1:h) !заменяем оригинальное на размытое
    level = level + deltaT
  end do

Тут все очевидно, единственное, что может быть не совсем очевидно для новичков в Fortran — это вложенность циклов, а именно — цикл по j идет раньше цикла по i. Все, опять же, из-за того, что фортран хранит двумерные массивы в памяти как последовательности столбцов. Если помять циклы местами, программа также будет работать, но значительно медленнее — программа будет обходить память не последовательно, а перебегать из разных ячеек «туда-сюда».

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

  deallocate(u)

  call pgmwriteheader(output, w, h) 
  call pgmappendbytes(output, nu, 1, 1) 

  deallocate(nu)

Процедура pgmwriteheader создает файл output и записывает в него заголовок PGM P5. Процедура pgmappendbytes записывает в конец файла output последовательность байт из nu, учитывая, что индексы nu начинаются с 1 по обоим измерениям. Замечу, что pgmappendbytes записывает байты из двумерного массива опять же в порядке столбцов, поэтому, хотя в памяти и находилось транспонированная версия изображения, при записи изображение транспонируется обратно.

Осталось определить функцию g. Например, реализуем функцию по формуле 3.

  contains 
      function g(x) result(v)
        implicit none
        double precision, intent(in) :: x
        double precision :: v
        v = 1/(1+(x/k)**2)
      end function g

Ссылка на полный код программы

Компилировал все через gfortran следующей командой (при условии, что все файлы лежат в одной директории):

gfortran pgmio.f90 blur.f90

Примеры работы фильтра


Для начала сравним размытие фильтра Перона и Малика с гауссовым размытие при
одинаковых t.

Исходное изображение.



Сгладим это изображение фильтром Гаусса и фильтром анизотропной диффузии.



Рассмотрим еще одно изображение.



Сглаживание фильтром Гаусса (t=10).



Сглаживание фильтром Перона и Малика (t=10, k=8).



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

Поэкспериментируем теперь с параметрами. Будем варьировать k при одинаковых t и deltaT (t=10, deltaT=0.2).



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

Посмотрим как поведет себя алгоритм для разных шагов по времени (t=10,k=5).



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



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

Заключение


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

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

Повторно приведу ссылки на исходные коды: модуль для работы с PGM, реализация фильтра.

Источники


  1. Оригинальная статья Перона и Малика. Scale-Space and Edge Detection Using Anisotropic Diffusion
  2. Численный расчет уравнения анизотропной диффузии
  3. Стандарт формата PGM
  4. Учебник по современному фортрану
Поделиться публикацией
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама
Комментарии 24
  • +2
    Судя по заявленным целям, было бы уместнее сравнивать не с Гауссом, а медианой?
    А вообще интересно.
    • +2
      Я бы сказал, что даже не с Гауссом, а с Bilateral.
      • 0
        Спасибо! На самом деле сравнивать так фильтры сложно, каждый выполняет свою задачу. Я привел для сравнения «на глаз» именно сглаживание Гаусса, потому что фильтр Перона-Малика по теории получается модификацией фильтра Гаусса (это на самом деле не совсем так, но когда функция c равняется единице, сглаживание фильтром анизотропной диффузии будет сходиться к гауссову сглаживанию).
      • 0
        Было бы интересно посмотреть, какой эффект дает последовательное применение этих алгоритмов. На взгляд дилетанта, они могут неплохо друг друга дополнять
        • +1
          http://gmic.eu/ тут огромное количество всяких продвинутых фильтров (анизотропное размытие точно есть) и минимально приличный интерфейс. Работает как плагин к GIMP
          • 0
            Спасибо за ссылку! Тут даже не просто фильтр Перона и Малика, а обобщенный для анизотропного размытия. У Перона и Малика ведь функция c — это скаляр, а вообще на месте функции c может стоять какой-то тензор. К тому же, если я правильно понял, плагин предлагает несколько вариантов дискретизации уравнения, что тоже не мало важно.
          • 0
            Это далеко не дилетантский взгляд, действительно — все зависит от задачи и исходных данных, особенно от типа шумов, ведь обычно их природа известна (хотя бы в теории). Вот, например, я сейчас занимаюсь обработкой изображений срезов керна (уже получены с томографа после обратного преобразования Радона) и там много шумов разной природы, даже больше скажу, некоторые «правильные» данные мы также трактуем как шумы, из-за их очень малой значимости для задачи в целом. Так вот, перед сегментацией мы фильтруем изображения кучей разных фильтров с различными параметрами, причем тут еще важен порядок.
            P.S. Параметры мы настраиваем исходя из уже известных данных.
            • 0
              Я к тому, что если будете еще писать на подобные темы — было бы интересно почитать не только про сравнение, но и про совмещение разлиных методов :)
        • +2

          В качестве практической реализации на C++ семейства подобных фильтров, могу порекомендовать библиотеку ITK. В ней наиболее полная реализация этого алгоритма обобщённого для N-мерных данных. Есть также реализация на OpenCL.
          https://itk.org/Doxygen/html/group__ITKAnisotropicSmoothing.html


          А также дополнительный модуль Anisotropic Diffusion LBR:
          https://github.com/InsightSoftwareConsortium/ITKAnisotropicDiffusionLBR
          Публикация: http://insight-journal.org/browse/publication/953

          • 0
            А существует ли обобщение на случай полноцветных изображений? Вроде бы применение к каждому каналу независимо должно приводить к несовсем «корректному» (ожидаемому) результату.
            • 0
              Если честно, не слышал про обобщение. Но если принять, что функция изображения I — вектор-функция,
              например I(x,y)=(r(x,y),g(x,y),b(x,y)), то из-за линейности операции дифференцирования для каждой координатной функции r,g,b будет записано свое уравнение анизотропной диффузии, поэтому, в целом, можно независимо по каналам считать. Правда, надо понимать, что при оценке границ в обычном фильтре мы используем градиент, и это просто вектор. А когда функция I — векторная, то градиент уже будет тензором, соответственно норма градиента будет другая и оценка этой нормы тоже. То есть поменяется расчет функции c. Если же применять в лоб алгоритм для каждого канала отдельно, то результат будет скорее всего «хуже», но все зависит от конкретной задачи.
              • +1

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

              • 0

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

              • +1

                А неявные схемы (чтобы увеличить шаг по времени — т.е. обсчитывать быстрее) пробовали? Или всё равно при увеличении Δt вылезает заметная погрешность, и выигрыша не получается (в смысле, ускорение за счёт доступного увеличения Δt "сжирается" усложнением схемы)?

                • 0
                  Не пробовал, на это есть несколько причин. Во-первых, как вы уже заметили — сильно возрастает погрешность, то есть, в неявной схеме придется решать системы нелинейных уравнений (в зависимости от аппроксимации, в лучшем случае это может быть система многочленов), их тоже надо как-то аппроксимировать, отсюда — появляются дополнительные вычисление, возрастает невязка решения в целом. Во-вторых, решение, приведенное у меня в статье, легко распараллелить на вычислительном кластере, cpu и gpu.
                • 0
                  Не пробовали ли Вы применять данный фильтр для сглаживания изображений текста перед OCR? Если да, то каковы результаты? Стоит ли пробовать вместо билатерального фильтра?
                  • +2

                    Ну вы же понимаете, что все от случая к случаю. Данный фильтр в фотошопе называется surface blur и присутствует там уже лет 8 точно.Можно сразу протестировать.

                    • 0
                      Да, я прекрасно это понимаю :)
                    • +2
                      Пробовал только ради интереса на одной картинке. Конечно, все зависит от конкретного алгоритма. Применять можно. Я взял исходное изображение из статьи Siarshai, прогнал через онлайновый файнридер результат анизотропной диффузии (t=15,k=14, функция g как у меня в реализации), результат билатерального фильтра (gmic online Spatial variance=3, Value variance=32, Iterations=6) и результат автора статьи. Ниже приведены исходное изображение, сглаженные и текст, распознанный файнридером. Ссылки на текст ведут на pastebin, изображения кликабельны.

                      Исходное:


                      Фильтр Перона и Малика (текст):


                      Билатеральный фильтр (текст):


                      Алгоритм предложенный Siarshai (текст):


                      Опять же повторюсь, что данные результаты не показатель качества работы фильтров.
                      • +1
                        Спасибо большое за проведённу работу.
                    • 0
                      "… некоторое размытие еще происходит, но начинают появляться артефакты." — это мы наблюдаем неустойчивость явной схемы. Сделайте шаги по времени раз в 100 помельче.
                      • 0
                        Ну и туда же. Судя по программе, вы аппроксимируете градиенты I в функции «с» или «g» направленными разностями. Пробовали использовать центральные разности?
                        • 0
                          Не совсем понимаю, зачем судить по программе, когда в разделе теории я написал по поводу аппроксимации. Норма градиента аппроксимируется длинной проекции градиента. Как я уже упоминал в статье, Перона и Малик предложили такую аппроксимацию, хотя она и не совсем соответствует исходному уравнению, потому что она достаточно «дешевая» в плане вычислений и использование более точных аппроксимаций качественно не меняет общий результат. Кстати, изначально я как раз через центральные разности делал, разницы практически никакой не было, по крайней мере для алгоритма сегментации.

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

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