快速导向滤波

快速导向滤波

导向滤波、双边滤波、最小二乘滤波并称三大保边滤波器。

原理推导

快速引导滤波可以在滤波的同时，保持边缘细节，效果与双边滤波类似，但算法效率更优。

输入图像p，引导图像I，输出图像q
引导滤波用到局部线性模型，等式两边取梯度，可知输出图q与引导图I具有类似的梯度特征，即边缘特性接近。

输出图像与输入图像关系

求出线性函数的系数，也就是线性回归，即希望拟合函数的输出值与真实值p之间的差距最小。

得到：


其中：μ_k为I在窗口w中的平均值，σ_k^2为在窗口w中的方差，|w|为窗口w中像素的数量，p ̅_k为p在窗口w中的均值。

伪代码

其中第5步方差为平方的均值减去均值的平方

快速导向滤波

通过降采样求出a和b，在升采样回原尺寸，并用来对 I 线性变换。这样时间可以由 O(N) 降低到 O(N/s^2) ，其中s是采样率。

代码

Mat guidedfilter(Mat &img, Mat &dst, int r, double esp)
{
    int row = img.rows;
    int col = img.cols;
    img.convertTo(img, CV_32FC1);
    dst.convertTo(dst, CV_32FC1);
    Mat boxResult, mean_I, mean_P, result, mean_IP, cov_IP, mean_II, var_I;
    Mat a, b, mean_a, mean_b;
    boxFilter(Mat::ones(Size(col, row), img.type()), boxResult, CV_32FC1, Size(r, r));//计算均值N
    boxFilter(img, mean_I, CV_32FC1, Size(r, r));
    mean_I = mean_I / boxResult;//计算导向均值mean_I 
    boxFilter(dst, mean_P, CV_32FC1, Size(r, r));
    mean_P = mean_P / boxResult;//计算原始均值mean_P 
    boxFilter(img.mul(dst), mean_IP, CV_32FC1, Size(r, r));
    mean_IP = mean_IP / boxResult;//计算互相关均值mean_IP corrIP
    cov_IP = mean_IP - mean_I.mul(mean_P);
    boxFilter(img.mul(img), mean_II, CV_32FC1, Size(r, r));
    mean_II = mean_II / boxResult;//计算自相关均值mean_II   corrI
    var_I = mean_II - mean_I.mul(mean_I);//计算方差var_I 
    a = cov_IP / (var_I + esp);//计算相关系数 
    b = mean_P - a.mul(mean_I);
    boxFilter(a, mean_a, CV_32FC1, Size(r, r));
    boxFilter(b, mean_b, CV_32FC1, Size(r, r));
    mean_a = mean_a / boxResult;//计算mean_a 
    mean_b = mean_b / boxResult;//计算mean_b
    result = mean_a.mul(img) + mean_b;//计算P
    return result;
}

优化快速导向滤波

```cpp
cv::Mat fastGuidedFilter(cv::Mat I_org, cv::Mat p_org, int r, double eps, int s)
{

    cv::Mat I, _I;
    _I = I_org.clone();
    resize(_I, I, Size(), 1.0 / s, 1.0 / s, 1);//下采样I

    cv::Mat p, _p;

    _p = p_org.clone();
    resize(_p, p, Size(), 1.0 / s, 1.0 / s, 1); //下采样P

    //[hei, wid] = size(I);    
    int hei = I.rows;
    int wid = I.cols;

    r = (2 * r + 1) / s + 1;//因为opencv自带的boxFilter（）中的Size,比如9x9,我们说半径为4   

                            //mean_I = boxfilter(I, r) ./ N;    
    cv::Mat mean_I;
    cv::boxFilter(I, mean_I, CV_32FC1, cv::Size(r, r)); //求I 的均值

    //mean_p = boxfilter(p, r) ./ N;    
    cv::Mat mean_p;
    cv::boxFilter(p, mean_p, CV_32FC1, cv::Size(r, r)); //求P 的均值

    //mean_Ip = boxfilter(I.*p, r) ./ N;    
    cv::Mat mean_Ip;                                     //求IP的均值
    cv::boxFilter(I.mul(p), mean_Ip, CV_32FC1, 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_32FC1, cv::Size(r, r));

    //var_I = mean_II - mean_I .* mean_I;    
    cv::Mat var_I = mean_II - mean_I.mul(mean_I);    //求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_32FC1, cv::Size(r, r));
    Mat rmean_a;
    resize(mean_a, rmean_a, Size(I_org.cols, I_org.rows), 1);//上采样

    //mean_b = boxfilter(b, r) ./ N;    
    cv::Mat mean_b;
    cv::boxFilter(b, mean_b, CV_32FC1, cv::Size(r, r));
    Mat rmean_b;
    resize(mean_b, rmean_b, Size(I_org.cols, I_org.rows), 1);//上采样

    //q = mean_a .* I + mean_b; % Eqn. (8) in the paper;    
    cv::Mat q = rmean_a.mul(_I) + rmean_b;

    return q;
}

运行

using namespace std;
using namespace cv;
Mat guidedfilter(Mat &img, Mat &dst, int r, double esp);
Mat fastGuidedFilter(cv::Mat I_org, cv::Mat p_org, int r, double eps, int s);
int main()
{
    clock_t start_t, end_t;
    Mat image, dst_image, temp_image, P, GQ, FGQ,MQ;
    vector<Mat>src, dst(3), dst2(3);
    image = imread("1.png");
    split(image, src);
    start_t = clock();
    for (int i = 0; i < 3; i++)
    {
        src[i].convertTo(temp_image, CV_32FC1, 1.0 / 255.0);
        P = temp_image.clone();
        dst_image = guidedfilter(temp_image, P, 10, 0.01);
        //dst_image = fastGuidedFilter(temp_image, P, 10, 0.01,8);
        dst_image.convertTo(dst_image, CV_8UC1, 255.0); 
        dst[i] = dst_image;
    }
    end_t = clock();
    double total_t = (double)(end_t - start_t);
    cout << "guidefilter_time: " << total_t << endl;
    start_t = clock();
    for (int i = 0; i < 3; i++)
    {
        src[i].convertTo(temp_image, CV_32FC1, 1.0 / 255.0);
        P = temp_image.clone();
        //dst_image = guidedfilter(temp_image, P, 10, 0.01);
        dst_image = fastGuidedFilter(temp_image, P, 10, 0.01, 8);
        dst_image.convertTo(dst_image, CV_8UC1, 255.0);
        dst2[i] = dst_image;
    }
    end_t = clock();
    double total_t2 = (double)(end_t - start_t);
    cout << "fastguidefilter_time: " << total_t2 << endl;
    cv::blur(image, MQ, Size(10, 10));
    merge(dst, GQ);
    merge(dst2, FGQ);
    imshow("src image", image);
    imshow("blur image", MQ);
    imshow("guidefilter image", GQ);
    imshow("fastguidefilter image", FGQ);
    waitKey(0);
    system("pause");
    return 0;

}

结果