跳转至

前置内容

在学习引导滤波之前,最好对高斯滤波和双边滤波有过理解,对于高斯滤波:W_{ij} = \frac{1}{K_i}exp(-\frac{|x_j-x_i|^2}{\sigma^2}),其中W是权重,ij是像素的索引,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;
}