前置内容¶
在学习引导滤波之前，最好对高斯滤波和双边滤波有过理解，对于高斯滤波:W_{ij} = \frac{1}{K_i}exp(-\frac{|x_j-x_i|^2}{\sigma^2})，其中W是权重，i和j是像素的索引，K是归一化的常量。公式可以看出，权重只和像素之间的空间距离有关系存在关系，滤波效果和图像内容无关。
所以为了增加对不同的图像有不同的滤波效果，引入了双边滤波，双边滤波的形式为:W_{i,j}=\frac{1}{K_i}exp(-\frac{|x_j-x_i|^2}{\sigma_s^2})exp(-\frac{|I_j-I_i|^2}{\sigma_r^2})，其中I是像素的强度值，所以双边滤波在强度差距大的地方(边缘)，权重会减小，滤波效应会减少。也就是说，双边滤波在像素强度变化不大的区域，有类似于高斯滤波的效果，而在图像边缘等像素强度梯度较大的地方可以保持梯度。
引导滤波¶
这个方法是何凯明提出的，在引导滤波中，首先利用了局部线性模型。这个模型认为某函数上一点与其近邻部分的点成线性关系，一个复杂的函数就可以用很多局部的线性函数来表示，当需要求该函数上某一点的值时，只需要计算所有包含该点的线性函数的值并取平均值即可。这种模型，在表示非解析函数上，非常有用。
同理，我们可以认为图像是一个二维函数，并且假设该函数的输出与输入在一个二维窗口内满足线性关系，如下：q_i=a_kI_i+b_k, \forall i \in w_k其中，q是输出像素的值，I是输入图像的值，i和k是像素索引，a和b是当窗口中心位于k时该线性函数的系数。
其实，输入图像不一定是待滤波的图像本身，也可以是其他图像即引导图像，这也是为何称为引导滤波的原因。对上式两边取梯度，可以得到\partial q = a \partial I即当输入图像I有梯度时，输出q也有类似的梯度，现在可以解释为什么引导滤波有边缘保持特性了。
下一步是求出线性函数的系数，也就是线性回归，即希望拟合函数的输出值与真实值p之间的差距最小，也就是让下式最小E(a_k, b_k)=\sum_{i\in w_k}{((a_kI_i + b_k - p_i)^2 + \varepsilon a_k^2)}，这里p只能是待滤波图像，并不像I那样可以是其他图像。
同时，a之前的系数用于防止求得的a过大，也是调节滤波器滤波效果的重要参数(相当于L2正则化的权重惩罚)。接下来利用最小二乘法的原理令\frac{\partial E}{\partial a_k}=0和\frac{\partial E}{b_k}=0得到2个二元一次方程，联立求解得到a_k = \frac{\frac{1}{w}\sum_{i\in w_k}I_ip_i-u_k\bar p_k}{\sigma_k^2+\varepsilon}，b_k=\bar p_k - a_ku_k，其中u_k是I在窗口w_k的平均值，\sigma_k^2是I在窗口w_k的方差，|w|是窗口w_k中的像素个数，\bar p_k是是待滤波图像p在窗口w_k中的均值。在计算每个窗口的线性系数时，我们可以发现一个像素会被多个窗口包含，也就是说，每个像素都由多个线性函数所描述。
因此，如之前所说，要具体求某一点的输出值时，只需将所有包含该点的线性函数值平均即可，如下:q_i=\frac{1}{|w|\sum_{k:i\in w_k}(a_kI_i+b_k)}=\bar a_iI_i + \bar b_i这里，w_k是所有包含像素i的窗口，k是其中心位置。
当把引导滤波用作边缘保持滤波器时，往往有 I = p ，如果\varepsilon =0，显然a=1, b=0是E(a,b)为最小值的解，(这里请参考协方差和方差的概念：https://baike.baidu.com/item/%E5%8D%8F%E6%96%B9%E5%B7%AE/2185936?fr=aladdin )从上式可以看出，这时的滤波器没有任何作用，将输入原封不动的输出。如果\varepsilon > 0，在像素强度变化小的区域（或单色区域），有a近似于（或等于）0，而b近似于（或等于）\bar p_k，即做了一个加权均值滤波；而在变化大的区域，a近似于1，b近似于0，对图像的滤波效果很弱，有助于保持边缘。而\varepsilon的作用就是界定什么是变化大，什么是变化小。在窗口大小不变的情况下，随着e的增大，滤波效果越明显。
在滤波效果上，引导滤波和双边滤波差不多，然后在一些细节上，引导滤波较好(在PS的磨皮美白中，经过亲生实践，效果更好)。引导滤波最大的优势在于，可以写出时间复杂度与窗口大小无关的算法，因此在使用大窗口处理图片时，其效率更高。
伪代码¶
C++代码实现¶
在别人博客上找了个代码，自己实现的不知道为什么滤波后只剩下图像的轮廓了，分析下原因，写出自己的代码再放上来。 以下代码来自：https://blog.csdn.net/wangyaninglm/article/details/44838545
Mat getimage(Mat &a)
{
int hei =a.rows;
int wid = a.cols;
Mat I(hei, wid, CV_64FC1);
//convert image depth to CV_64F
a.convertTo(I, CV_64FC1,1.0/255.0);
//normalize the pixel to 0~1
/*
for( int i = 0; i< hei; i++){
double *p = I.ptr<double>(i);
for( int j = 0; j< wid; j++){
p[j] = p[j]/255.0;
}
}
*/
return I;
}
Mat cumsum(Mat &imSrc, int rc)
{
if(!imSrc.data)
{
cout << "no data input!\n" << endl;
}
int hei = imSrc.rows;
int wid = imSrc.cols;
Mat imCum = imSrc.clone();
if( rc == 1)
{
for( int i =1;i < hei; i++)
{
for( int j = 0; j< wid; j++)
{
imCum.at<double>(i,j) += imCum.at<double>(i-1,j);
}
}
}
if( rc == 2)
{
for( int i =0;i < hei; i++)
{
for( int j = 1; j< wid; j++)
{
imCum.at<double>(i,j) += imCum.at<double>(i,j-1);
}
}
}
return imCum;
}
Mat boxfilter(Mat &imSrc, int r)
{
int hei = imSrc.rows;
int wid = imSrc.cols;
Mat imDst = Mat::zeros( hei, wid, CV_64FC1);
//imCum = cumsum(imSrc, 1);
Mat imCum = cumsum(imSrc,1);
//imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
for( int i = 0; i<r+1; i++)
{
for( int j=0; j<wid; j++ )
{
imDst.at<double>(i,j) = imCum.at<double>(i+r,j);
}
}
//imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
for( int i =r+1; i<hei-r;i++)
{
for( int j = 0; j<wid;j++)
{
imDst.at<double>(i,j) = imCum.at<double>(i+r,j)-imCum.at<double>(i-r-1,j);
}
}
//imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);
for( int i = hei-r; i< hei; i++)
{
for( int j = 0; j< wid; j++)
{
imDst.at<double>(i,j) = imCum.at<double>(hei-1,j)-imCum.at<double>(i-r-1,j);
}
}
imCum = cumsum(imDst, 2);
//imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
for( int i = 0; i<hei; i++)
{
for( int j=0; j<r+1; j++ )
{
imDst.at<double>(i,j) = imCum.at<double>(i,j+r);
}
}
//imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
for( int i =0 ; i<hei;i++)
{
for( int j = r+1; j<wid-r ;j++ )
{
imDst.at<double>(i,j) = imCum.at<double>(i,j+r)-imCum.at<double>(i,j-r-1);
}
}
//imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
for( int i = 0; i< hei; i++)
{
for( int j = wid-r; j<wid; j++)
{
imDst.at<double>(i,j) = imCum.at<double>(i,wid-1)-imCum.at<double>(i,j-r-1);
}
}
return imDst;
}
Mat guidedfilter( Mat &I, Mat &p, int r, double eps )
{
int hei = I.rows;
int wid = I.cols;
//N = boxfilter(ones(hei, wid), r);
Mat one = Mat::ones(hei, wid, CV_64FC1);
Mat N = boxfilter(one, r);
//mean_I = boxfilter(I, r) ./ N;
Mat mean_I(hei, wid, CV_64FC1);
divide(boxfilter(I, r), N, mean_I);
//mean_p = boxfilter(p, r) ./ N;
Mat mean_p(hei, wid, CV_64FC1);
divide(boxfilter(p, r), N, mean_p);
//mean_Ip = boxfilter(I.*p, r) ./ N;
Mat mul_Ip(hei, wid, CV_64FC1);
Mat mean_Ip(hei, wid, CV_64FC1);
multiply(I,p,mul_Ip);
divide(boxfilter(mul_Ip, r), N, mean_Ip);
//cov_Ip = mean_Ip - mean_I .* mean_p
//this is the covariance of (I, p) in each local patch.
Mat mul_mean_Ip(hei, wid, CV_64FC1);
Mat cov_Ip(hei, wid, CV_64FC1);
multiply(mean_I, mean_p, mul_mean_Ip);
subtract(mean_Ip, mul_mean_Ip, cov_Ip);
//mean_II = boxfilter(I.*I, r) ./ N;
Mat mul_II(hei, wid, CV_64FC1);
Mat mean_II(hei, wid, CV_64FC1);
multiply(I, I, mul_II);
divide(boxfilter(mul_II, r), N, mean_II);
//var_I = mean_II - mean_I .* mean_I;
Mat mul_mean_II(hei, wid, CV_64FC1);
Mat var_I(hei, wid, CV_64FC1);
multiply(mean_I, mean_I, mul_mean_II);
subtract(mean_II, mul_mean_II, var_I);
//a = cov_Ip ./ (var_I + eps);
Mat a(hei, wid, CV_64FC1);
for( int i = 0; i< hei; i++){
double *p = var_I.ptr<double>(i);
for( int j = 0; j< wid; j++){
p[j] = p[j] +eps;
}
}
divide(cov_Ip, var_I, a);
//b = mean_p - a .* mean_I;
Mat a_mean_I(hei ,wid, CV_64FC1);
Mat b(hei ,wid, CV_64FC1);
multiply(a, mean_I, a_mean_I);
subtract(mean_p, a_mean_I, b);
//mean_a = boxfilter(a, r) ./ N;
Mat mean_a(hei, wid, CV_64FC1);
divide(boxfilter(a, r), N, mean_a);
//mean_b = boxfilter(b, r) ./ N;
Mat mean_b(hei, wid, CV_64FC1);
divide(boxfilter(b, r), N, mean_b);
//q = mean_a .* I + mean_b;
Mat mean_a_I(hei, wid, CV_64FC1);
Mat q(hei, wid, CV_64FC1);
multiply(mean_a, I, mean_a_I);
add(mean_a_I, mean_b, q);
return q;
}
/*****************
http://research.microsoft.com/en-us/um/people/kahe/eccv10/
推酷上的一篇文章：
http://www.tuicool.com/articles/Mv2iiu
************************/
cv::Mat guidedFilter2(cv::Mat I, cv::Mat p, int r, double eps)
{
/*
% GUIDEDFILTER O(1) time implementation of guided filter.
%
% - 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
*/
cv::Mat _I;
I.convertTo(_I, CV_64FC1);
I = _I;
cv::Mat _p;
p.convertTo(_p, CV_64FC1);
p = _p;
//[hei, wid] = size(I);
int hei = I.rows;
int wid = I.cols;
//N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
cv::Mat N;
cv::boxFilter(cv::Mat::ones(hei, wid, I.type()), N, CV_64FC1, cv::Size(r, r));
//mean_I = boxfilter(I, r) ./ N;
cv::Mat mean_I;
cv::boxFilter(I, mean_I, CV_64FC1, cv::Size(r, r));
//mean_p = boxfilter(p, r) ./ N;
cv::Mat mean_p;
cv::boxFilter(p, mean_p, CV_64FC1, cv::Size(r, r));
//mean_Ip = boxfilter(I.*p, r) ./ N;
cv::Mat mean_Ip;
cv::boxFilter(I.mul(p), mean_Ip, CV_64FC1, cv::Size(r, r));
//cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.
cv::Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
//mean_II = boxfilter(I.*I, r) ./ N;
cv::Mat mean_II;
cv::boxFilter(I.mul(I), mean_II, CV_64FC1, cv::Size(r, r));
//var_I = mean_II - mean_I .* mean_I;
cv::Mat var_I = mean_II - mean_I.mul(mean_I);
//a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
cv::Mat a = cov_Ip/(var_I + eps);
//b = mean_p - a .* mean_I; % Eqn. (6) in the paper;
cv::Mat b = mean_p - a.mul(mean_I);
//mean_a = boxfilter(a, r) ./ N;
cv::Mat mean_a;
cv::boxFilter(a, mean_a, CV_64FC1, cv::Size(r, r));
mean_a = mean_a/N;
//mean_b = boxfilter(b, r) ./ N;
cv::Mat mean_b;
cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(r, r));
mean_b = mean_b/N;
//q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
cv::Mat q = mean_a.mul(I) + mean_b;
return q;
}
本文总阅读量次