现在从一个最简单的情形来开始我们的讨论。假设有一个原始图像
具体而言,在简单平滑中,图像
中第
个像素是由图像
中以第
个像素为中心的一个窗口
中的所有像素取平均而得来的,即
其中,
,
是窗口
中的像素数目。也就是说,在简单平滑中,以像素
为中心的一个窗口
中的像素具有等同的权值。但是在高斯平滑中,权值
将服从二维的高斯分布,结果导致离像素
更接近的像素将具有更高的权重,反之离像素
较远的像素则具有更小的权重。
无论是简单平滑,还是高斯平滑,它们都有一个共同的弱点,即它们都属于各向同性滤波。我们都知道,一幅自然的图像可以被看成是有(过渡平缓的,也就是梯度较小)区域和(过渡尖锐的,也就是梯度较大)边缘(也包括图像的纹理、细节等)共同组成的。噪声是影响图像质量的不利因素,我们希望将其滤除。噪声的特点通常是以其为中心的各个方向上梯度都较大而且相差不多。边缘则不同,边缘相比于区域也会出现梯度的越变,但是边缘只有在其法向方向上才会出现较大的梯度,而在切向方向上梯度较小。
因此,对于各向同性滤波(例如简单平滑或高斯平滑)而言,它们对待噪声和边缘信息都采取一直的态度。结果,噪声被磨平的同时,图像中具有重要地位的边缘、纹理和细节也同时被抹平了。这是我们所不希望看到的。研究人员已经提出了很多 Edge-perserving 的图像降噪(平滑)算法,例如双边滤波、自适应(维纳)平滑滤波(请参见文献【1】)、基于 PM 方程的各向异性滤波以及基于 TV-norm 的降噪算法等。本文将考虑在文献【2】中提出的另外一种 Edge-perserving 的图像滤波(平滑)算法——导向滤波(Guided Filter)。当然,通过阅读文献【2】,我们也知道导向滤波的应用不止有 Edge-perserving 的图像平滑,还包括图像去雾、图像 Matting 等等。
欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji,为保证公式、图表得以正确显示,强烈建议你从该地址上查看原版博文。本博客主要关注方向包括:数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。
导向滤波之所以叫这个名字,因为在算法框架中,要对
注意,这里的
就表示由引导图像
来确定加权平均运算中所采用的权值。注意,引导图像可以是单独的一幅图像,也可以输入的图像
本身。当引导图像就是
本身的时候,导向滤波就变成了一个 Edge-perserving 的滤波器,因此可以用于图像的平滑降噪。
导向滤波滤波的示意图如下所示(图片来自作者原文【2】)。注意这也是导向滤波所依赖的一个重要假设——导向滤波器在导向图像
你可以想象一种情况来理解这种假设。
中的甲处和乙处都是区域,而丙处是边缘。那甲处的梯度和乙处的梯度应该不会相差太大,例如乙处的梯度是甲处梯度的 1.5 倍。但是由于丙处是边缘,所以丙处的梯度会比较大,例如可能是甲处的 6 倍,也就是乙处的 4 倍。如果
和
之间满足线性关系,那么甲乙丙的梯度大小和倍数关系就都不会被扭曲。否则,如果二者之间的关系是非线性的,那么可能的结果是在
中,尽管丙处的梯度仍然大于乙处的梯度,进而大于甲处的梯度,但是倍数关系可能会扭曲。例如,丙处的梯度是乙处的 1.5 倍,而是甲处的 6 倍(即乙处的梯度是甲处的 4 倍),这时你就会想象,甲处是区域,而乙处和丙处就变成了边缘。可见非线性关系会使得引导图像对于边缘和区域的指示作用发生错乱。
现在已知的是
于是便得到了一个最小二乘问题。但是我们也知道普通最小二乘有时候引起一些麻烦,因此要适当地引入惩罚项,即采用正则化的手段。原作者在处理这里时所采用的方法借鉴了从普通线性回归改进到岭回归时所采用的方法(如果你对正则化、岭回归、或者最小二乘法还不够清楚,那么你可以参考文献【5】),即求解下面这个最优化(最小化)目标所对应的参数
和
。
求解上述最优化问题(跟最小二乘法的推导过程一致),便会得到:
其中,
是
中窗口
中的平均值,
是
中窗口
中的方差,
是窗口
中像素的数量,
是待滤波图像
在窗口
中的均值,即
。
此外,在计算每个窗口的线性系数时,我们可以发现一个像素会被多个窗口包含,也就是说,每个像素都由多个线性函数所描述。因此,如之前所说,要具体求某一点的输出值时,只需将所有包含该点的线性函数值平均即可,如下
其中,
,
。
下面给出导向滤波算法的流程,
为一个窗口半径为
的均值滤波器(对应的窗口大小为
),corr 为相关,var 为方差,cov 为协方差。
下面来在 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
上述代码的实现完全遵照上一小节最后给出的算法流程图。这里需要略作解释的地方是函数 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
从上述代码可以看到,当参数 r=1 时,此时的滤波器窗口是
- mean_I = boxfilter(I, r). / N;
才算是完成了均值滤波。而这整个过程就相当于文献【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
上述代码基于积分图实现,如果你对积分图不是很了解,可以参考文献【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);
执行上述代码,结果如下所示。可见效果还是很不错的。但是我们可以来做一下事后分析,看看导向滤波是如果实现 edge-perserving 的平滑滤波效果的。当
考虑两种情况:
也就是说在方差比较大的区域,保持值不变,在平滑区域,使用临近像素平均(也就退化为普通均值滤波)。(这个思想跟文献【1】里设计的自适应降噪滤波器有异曲同工之妙,当然也不完全相同!)
上面的给出的是对灰度图像进行导向滤波的代码。你可能会好奇彩色图像该如何使用导向滤波。一个比较直接的方法就是将引导滤波分别应用到 RGB 三个颜色通道中,然后在组合成结果图像。更多细节可以参考作者原文。不仅如此,在本文开始我们也谈到,导向滤波不仅可以应用与图像平滑,还可以用于图像去雾、图像 Matting 等多个领域,限于篇幅我们无法一一详述,有兴趣的读者可以参考作者原文以了解导向滤波的其他应用。
最后,需要说明的是在新版的 MATLAB 中(R2014 及以后),已经内置了用于导向滤波的函数 imguidedfilter(而这个函数实现的其实是 Fast Guided Filter),也就是说在实际开发中我们已经不再需要编写上面那样的代码,而是只要简单调用 MATLAB 的内置函数就可以了。这个函数的使用方法可以参考文献【3】,本文不再赘述。
【1】 自适应图像降噪滤波器的设计与实现
【2】 Kaiming He, Jian Sun, Xiaoou Tang, Guided Image Filtering. IEEE Transactions on Pattern Analysis and Machine Intelligence, Volume 35, Issue 6, pp. 1397-1409, June 2013
【3】 MATLAB 中的导向滤波函数 imguidedfilter
【4】 约束复原与维纳滤波
【5】 R 语言实战——机器学习与数据分析(P270-P274)
【6】 BoxFilter 的解析
【7】 机器视觉中的图像积分图及其实现
来源: http://blog.csdn.net/baimafujinji/article/details/74750283