引导滤波器（原理及opencv代码）

这篇博客整理自以下博客内容：

https://blog.csdn.net/sinat_36264666/article/details/77990790

https://blog.csdn.net/piaoxuezhong/article/details/78372787

引导滤波定义：

即需要引导图的滤波器，引导图可以是单独的图像或者是输入图像，当引导图为输入图像时，引导滤波就成为一个保持边缘的滤波操作，可以用于图像重建的滤波。

 

引导滤波与其他方向相比优势所在：

像高斯滤波等线性滤波算法所用的核函数相对于待处理的图像是独立无关的，也就意味着，对任意图像都是采用相同的操作。但是，有时候我们是希望在滤波过程中可以加入引导图像中的信息的，例如，在上色处理过程中，结果图像的色度通道需要包含跟给定亮度通道一致的连续边缘。

 

引导滤波的主要特点及应用：

1、引导滤波（导向滤波）是一种图像滤波技术，通过一张引导图I，对初始图像p（输入图像）进行滤波处理，使得最后的输出图像大体上与初始图像P相似，但是纹理部分与引导图I相似。其典型应用有两个：保边图像平滑，抠图。

2、 引导滤波（导向滤波）的目的是，保持双边滤波的优势（有效保持边缘，非迭代计算），而克服双边滤波的缺点（设计一种时间复杂度为 O(1) 的快速滤波器，而且在主要边缘附近没有梯度的变形）。

3、 引导滤波（导向滤波）不仅能实现双边滤波的边缘平滑，而且在检测到边缘附近有很好的表现，可应用在图像增强、HDR压缩、图像抠图及图像去雾等场景。

 

引导滤波的算法原理：

第一步：

假设该引导滤波函数的输出与输入在一个二维窗口内满足线性关系，如下：

其中，q是输出像素的值，I是输入图像的值，i和k是像素索引，a和b是当窗口中心位于k时该线性函数的系数。（当引导图为输入图像时，引导滤波就成为一个保持边缘的滤波操作，即 I= p）

 

对上式两边取梯度，可以得到：

即当输入图I有梯度时，输出 q 也有类似的梯度，现在可以解释为什么引导滤波有边缘保持特性了。

q 即 p去除噪声或者纹理之后的图像：

其中ni表示噪声。

 

第二步：

求出线性函数的系数，也就是线性回归，即希望拟合函数的输出值 q与真实值 p 之间的差距最小，转化为最优化问题，也就是让下式最小：

这里p只能是待滤波图像，并不像I那样可以是其它图像。同时，a之前的系数是用于防止求得的a过大，也是调节滤波器滤波效果的重要参数。

通过最小二乘法，可以得到：

其中， 是图像  在窗口 中的平均值， 是待滤波图像p在窗口  中的均值， 是 在窗口 中的方差， 是窗口 中像素的数量。

 

参数的详解：

1、如果=0，显然a=1, b=0是E(a,b)为最小值的解，从上式可以看出，这时的滤波器没有任何作用，将输入原封不动的输出。

2、如果>0，在像素强度变化小的区域（或单色区域），有a近似于（或等于）0，而b近似于（或等于），即做了一个加权均值滤波；

3、 而在变化大的区域，a近似于1，b近似于0，对图像的滤波效果很弱，有助于保持边缘。

4、 而的作用就是界定什么是变化大，什么是变化小。在窗口大小不变的情况下，随着的增大，滤波效果越明显。

 

第三步：

在计算每个窗口的线性系数时，我们可以发现一个像素会被多个窗口包含，也就是说，每个像素都由多个线性函数所描述。因此，如之前所说，要具体求某一点的输出值时，只需将所有包含该点的线性函数值平均即可，如下：

这里，是所有包含像素i的窗口，k是其中心位置。

 

引导滤波的优势：在滤波效果上，引导滤波和双边滤波差不多，在一些细节上，引导滤波较好。引导滤波最大的优势在于，可以写出时间复杂度与窗口大小无关的算法，因此在使用大窗口处理图片时，其效率更高。

 

整体实现步骤如下：

• 1.利用boxFilter滤波器完成相关系数参数，其中均值包括引导图像均值、原始待滤波图像均值、互相关均值及自相关均值。
• 2.根据均值计算相关系数参数，包括自相关方差var，互相关协方差cov。
• 3.计算窗口线性变换参数系数a、b。
• 4.根据公式计算参数a、b的均值。
• 5.利用参数得到引导滤波输出图像q

 

算法解释：

1、

先计算的分子部分：

• 分子第一项， Ip是在窗口中的和，再除以窗口中像素的个数，刚好就是盒式滤波，因此我们可以将输入的引导图像I和待滤波图像 p 相乘，并对相乘后的图像做box filtering，即得第一项的结果。
• 分子第二项 ，和分别为I 和 p 在窗口中均值，因此分别对I和p进行box filtering，再将box filtering之后的结果相乘即可。

实际上，的分子就是Ip在窗口中的协方差。

 

2、是引导图I在窗口中的方差，方差和期望（均值）之间是有关系的，如下式：

     

• 因此在计算I的方差时，我们可以先计算|w|的均值，再减去I均值的平方即的平方。在方法上，计算的|w|均值和计算Ip的均值是一样的。
• 最后，对计算出来的方差图像，加上常量（每个元素都加），分母就计算完了，自然，在所有窗口中的值也就得到了。

 

3、的计算类似：

• 第一项：为 p 在窗口中均值，因此对 p 进行box filtering。
• 第二项：为I在窗口中均值，因此对进行box filtering，并与相乘。

第一项减去第二项即为的值。

 

4、第二个公式，即我们最后要得到的滤波后图像 q ：

其中，输出值q又与两个均值有关，分别为a和b在窗口中的均值（不是），所以还是box filtering，我们将上一步得到两个图像都进行盒式滤波，得到两个新图：和。然后用乘以引导图像Ii，再加上，即得最终滤波之后的输出图像q。 

结论：当窗口处于平坦区域时，图像的局部方差小，则趋于0，趋于均值，相当于对图像进行了均值滤波；当窗口处于边缘时，局部方差较大，则趋于1，趋于0，滤波器输出相当于原图。这样就可以很好的保护边缘。一般地，图像中间，像素跳跃变化不大，趋于稳定，方差小；图像边缘，像素跳跃变化大，不稳定，方差大。

 

代码解释：

Mat.mul(Mat):

两个矩阵维数必须相同，对应位置元素相乘。

[0, 1, 2, 3;           [ 0, 1, 2, 3;       [0, 1, 4, 9;

0, 1, 2, 3;   .mul         0, 1, 2, 3;   =      0, 1, 4, 9;

 0, 1, 2, 3]            0, 1, 2, 3]         0, 1, 4, 9;]

 

Mat类型转换 ConvertTo：读取一个图片内容后要把图片内容的像素信息转为浮点并把当前的mat作为矩形进行矩阵计算。

代码：

#include <opencv2/core/core.hpp>                    
#include <opencv2/highgui/highgui.hpp>        
#include <opencv2/imgproc/imgproc.hpp>    
#include <iostream>       
using namespace std;  
using namespace cv; 

Mat guidedFilter(Mat &srcMat, Mat &guidedMat, int radius, double eps);//引导滤波器  
  
int main()  
{  
    //------------【0】定义相关变量-------------  
    Mat resultMat;  //最后结果图像  
    vector<Mat> vSrcImage, vResultImage;  
    //------------【1】读取源图像并检查图像是否读取成功------------       
    Mat srcImage = imread("image//left1.png");  
    if (!srcImage.data)  
    {  
        cout << "读取图片错误，请重新输入正确路径！\n";  
        system("pause");  
        return -1;  
    }  
    imshow("【源图像】", srcImage);  
    //-------【2】对源图像进行通道分离，并对每个分通道进行导向滤波操------  
    split(srcImage, vSrcImage);  
    for (int i = 0; i < 3; i++)  
    {  
        Mat tempImage;  
        vSrcImage[i].convertTo(tempImage, CV_64FC1, 1.0 / 255.0);//将分通道转换成浮点型数据  
        Mat cloneImage = tempImage.clone(); //将tempImage复制一份到cloneImage  
        Mat resultImage = guidedFilter(tempImage, cloneImage, 9, 0.01);//对分通道分别进行导向滤波，半径为1、3、5...等奇数  
        vResultImage.push_back(resultImage);//将分通道导向滤波后的结果存放到vResultImage中  
    }  
    //----------【3】将分通道导向滤波后结果合并-----------------------  
    merge(vResultImage, resultMat);  
    imshow("【引导滤波/导向滤波】", resultMat);  
    waitKey(0);  
    return 0;  
}  
  
//导向滤波器  
Mat guidedFilter(Mat &srcMat, Mat &guidedMat, int radius, double eps)  
{  
    //------------【0】转换源图像信息，将输入扩展为64位浮点型，以便以后做乘法------------  
    srcMat.convertTo(srcMat, CV_64FC1);  
    guidedMat.convertTo(guidedMat, CV_64FC1);  
    //--------------【1】各种均值计算----------------------------------  
    Mat mean_p, mean_I, mean_Ip, mean_II;  
    boxFilter(srcMat, mean_p, CV_64FC1, Size(radius, radius));//生成待滤波图像均值mean_p   
    boxFilter(guidedMat, mean_I, CV_64FC1, Size(radius, radius));//生成引导图像均值mean_I     
    boxFilter(srcMat.mul(guidedMat), mean_Ip, CV_64FC1, Size(radius, radius));//生成互相关均值mean_Ip  
    boxFilter(guidedMat.mul(guidedMat), mean_II, CV_64FC1, Size(radius, radius));//生成引导图像自相关均值mean_II  
    //--------------【2】计算相关系数，计算Ip的协方差cov和I的方差var------------------  
    Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);  
    Mat var_I = mean_II - mean_I.mul(mean_I);  
    //---------------【3】计算参数系数a、b-------------------  
    Mat a = cov_Ip / (var_I + eps);  
    Mat b = mean_p - a.mul(mean_I);  
    //--------------【4】计算系数a、b的均值-----------------  
    Mat mean_a, mean_b;  
    boxFilter(a, mean_a, CV_64FC1, Size(radius, radius));  
    boxFilter(b, mean_b, CV_64FC1, Size(radius, radius));  
    //---------------【5】生成输出矩阵------------------  
    Mat dstImage = mean_a.mul(srcMat) + mean_b;  
    return dstImage;  
}