导向滤波(Guided Filter)的解析与实现

时间：2024-03-20 22:08:08

现在从一个最简单的情形来开始我们的讨论。假设有一个原始图像 pp较远的像素则具有更小的权重。

无论是简单平滑，还是高斯平滑，它们都有一个共同的弱点，即它们都属于各向同性滤波。我们都知道，一幅自然的图像可以被看成是有（过渡平缓的，也就是梯度较小）区域和（过渡尖锐的，也就是梯度较大）边缘（也包括图像的纹理、细节等）共同组成的。噪声是影响图像质量的不利因素，我们希望将其滤除。噪声的特点通常是以其为中心的各个方向上梯度都较大而且相差不多。边缘则不同，边缘相比于区域也会出现梯度的越变，但是边缘只有在其法向方向上才会出现较大的梯度，而在切向方向上梯度较小。

因此，对于各向同性滤波（例如简单平滑或高斯平滑）而言，它们对待噪声和边缘信息都采取一直的态度。结果，噪声被磨平的同时，图像中具有重要地位的边缘、纹理和细节也同时被抹平了。这是我们所不希望看到的。研究人员已经提出了很多Edge-perserving的图像降噪（平滑）算法，例如双边滤波、自适应（维纳）平滑滤波（请参见文献【1】）、基于PM方程的各向异性滤波以及基于TV-norm的降噪算法等。本文将考虑在文献【2】中提出的另外一种Edge-perserving的图像滤波（平滑）算法——导向滤波（Guided Filter）。当然，通过阅读文献【2】，我们也知道导向滤波的应用不止有Edge-perserving的图像平滑，还包括图像去雾、图像Matting等等。

欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji，为保证公式、图表得以正确显示，强烈建议你从该地址上查看原版博文。本博客主要关注方向包括：数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。

算法原理

导向滤波之所以叫这个名字，因为在算法框架中，要对pp本身的时候，导向滤波就变成了一个Edge-perserving的滤波器，因此可以用于图像的平滑降噪。

导向滤波滤波的示意图如下所示（图片来自作者原文【2】）。注意这也是导向滤波所依赖的一个重要假设——导向滤波器在导向图像II中，尽管丙处的梯度仍然大于乙处的梯度，进而大于甲处的梯度，但是倍数关系可能会扭曲。例如，丙处的梯度是乙处的1.5倍，而是甲处的6倍（即乙处的梯度是甲处的4倍），这时你就会想象，甲处是区域，而乙处和丙处就变成了边缘。可见非线性关系会使得引导图像对于边缘和区域的指示作用发生错乱。

现在已知的是II），corr为相关，var为方差，cov为协方差。

导向滤波(Guided Filter)的解析与实现

基于MATLAB的算法实现

下面来在MATLAB中具体实现一下导向滤波算法（代码来自Dr. Kaiming He，使用时请尊重原作者权利）。

function q = guidedfilter(I, p, r, eps)

%   - guidance image: I (should be a gray-scale/single channel image)
%   - filtering input image: p (should be a gray-scale/single channel image)
%   - local window radius: r
%   - regularization parameter: eps

[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); 

mean_I = boxfilter(I, r) ./ N;
mean_p = boxfilter(p, r) ./ N;
mean_Ip = boxfilter(I.*p, r) ./ N;
% this is the covariance of (I, p) in each local patch.
cov_Ip = mean_Ip - mean_I .* mean_p; 

mean_II = boxfilter(I.*I, r) ./ N;
var_I = mean_II - mean_I .* mean_I;

a = cov_Ip ./ (var_I + eps); 
b = mean_p - a .* mean_I; 

mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;

q = mean_a .* I + mean_b; 
end
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27

上述代码的实现完全遵照上一小节最后给出的算法流程图。这里需要略作解释的地方是函数boxfilter，它是基于积分图算法实现的Box Filter。关于Box Filter，你也可以参考文献【6】以了解更多。首先来看看它到底做了些什么（因为这个Box Filter 和通常意义上的均值滤波并不完全一样）。

A =

     1     1     1
     1     1     1
     1     1     1

>> B = boxfilter(A, 1);
>> B

B =

     4     6     4
     6     9     6
     4     6     4
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14

从上述代码可以看到，当参数r=1时，此时的滤波器窗口是3×33×3。将这样大小的一个窗口扣在原矩阵中心，刚好可以覆盖所有矩阵，此时求和为9，即把窗口里覆盖到的值全部加和。此外当把窗口中心挪动到左上角的像素时，因为窗口覆盖区域里只有4个数字，所以结果为4。所以你可以看出这里的Box Filter只是做了求和处理，并没有归一化，所以并不是真正的均值滤波（简单平滑）。必须结合后面的一句

mean_I = boxfilter(I, r) ./ N;
  • 1

才算是完成了均值滤波。而这整个过程就相当于文献【6】中介绍的函数imboxfilt。但是你会发现它们二者在执行的时候最终的结果（主要是位于图像四周边缘的数值）会有细微的差异。这是因为MATLAB中的imboxfilt函数在处理位于图像四周边缘的像素时，需要虚拟地为原图像补齐滤波窗口覆盖但是没有值的区域。

下面给出上述boxfilter函数的实现代码。

function imDst = boxfilter(imSrc, r)

%   BOXFILTER   O(1) time box filtering using cumulative sum
%
%   - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
%   - Running time independent of r; 
%   - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);
%   - But much faster.

[hei, wid] = size(imSrc);
imDst = zeros(size(imSrc));

%cumulative sum over Y axis
imCum = cumsum(imSrc, 1);
%difference over Y axis
imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);

%cumulative sum over X axis
imCum = cumsum(imDst, 2);
%difference over Y axis
imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
end
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26

上述代码基于积分图实现，如果你对积分图不是很了解，可以参考文献【7】，这里不再赘述。

下面我们来实验一下上述导引滤波用于edge-perserving的平滑滤波效果。

I = double(imread('cat.bmp')) / 255;
p = I;
r = 4; % try r=2, 4, or 8
eps = 0.2^2; % try eps=0.1^2, 0.2^2, 0.4^2

O = guidedfilter(I, p, r, eps);

subplot(121), imshow(I);
subplot(122), imshow(O);
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9

执行上述代码，结果如下所示。可见效果还是很不错的。但是我们可以来做一下事后分析，看看导向滤波是如果实现edge-perserving的平滑滤波效果的。当I=pI=p考虑两种情况：

  • 情况1：高方差区域，即表示图像II中基本保持固定，此时有σ2k<<ϵ

现在从一个最简单的情形来开始我们的讨论。假设有一个原始图像 pp较远的像素则具有更小的权重。

相关文章