【图像处理】引导滤波（guided image filtering）——附C++、python实现彩图与灰度图代码

古城微笑少年丶 2021-10-29 23:31 1083阅读 0赞

### 【fishing-pan：[https://blog.csdn.net/u013921430][https_blog.csdn.net_u013921430] 转载请注明出处】 ###

## 前言 ##

引导滤波是何凯明等人在2010年提出的一种滤波方法，该方法与之前博客中介绍的双边滤波（回顾点[这里][Link 1]）都属于边缘保护滤波。引导滤波除了可以用于图像平滑，还可以用于HDR压缩、细节增强、图像去雾、联合上采样等图像处理任务。

本文主要根据论文中的内容介绍引导滤波的数学推断过程，并用C++与Python（因为最近在学Python）分别实现彩色图像与灰度图的引导滤波。

## 引导滤波（guided image filtering） ##

引导滤波的思想用一张引导图像产生权重，从而对输入图像进行处理，这个过程可以表示为公式  ( 1 ) (1) (1) 中的内容。  
 q i = ∑ j W i j ( I ) ⋅ p j                 ( 1 ) q\_\{i\}=\\sum\_\{j\} W\_\{ij\}\\left ( I \\right )\\cdot p\_\{j\}\\, \\, \\, \\, \\, \\, \\,\\: \\: \\: \\: \\: \\: \\: \\: (1) qi=j∑Wij(I)⋅pj(1)  
  公式中  q q q、 I I I、 p p p分表表示输出图像、引导图像和输入图像 ， i i i、 j j j 分别表示图像中像素点的索引。可以看到上方公式中权重  W W W 仅与引导图像  I I I 有关，而在双边滤波中权重  W W W 由输入图像自身决定。下面将根据论文中的内容讲述算法的推导过程。

引导滤波的核心的假设就是输出图像与引导图像在局部上是一个线性的模型，在局部窗口  w k w\_\{k\} wk 中，输出图像与引导图像的线性关系可以表示成式子  ( 2 ) (2) (2) 中的形式，而 ( a k ， b k ) (a\_\{k\}，b\_\{k\}) (ak，bk) 是窗口中的线性系数，在每个窗口中是常数。  
 q i = a k I i + b k , ∀ i ∈ w k                 ( 2 ) q\_\{i\}=a\_\{k\}I\_\{i\}+b\_\{k\},\\forall i\\in w\_\{k\}\\, \\, \\, \\, \\, \\, \\,\\: \\: \\: \\: \\: \\: \\: \\: (2) qi=akIi\+bk,∀i∈wk(2)  
  根据这个模型可以保证输出图像  q q q 受到引导图像  I I I 的约束，因为  ▽ q = a ▽ I \\triangledown q=a\\triangledown I ▽q=a▽I ，即当引导图像在局部有梯度变化时，输出图像对应地也会有梯度变化。  
![在这里插入图片描述][20190818080123999.PNG_pic_center]  
  为了求解式子  ( 2 ) (2) (2) 中的系数  ( a k ， b k ) (a\_\{k\}，b\_\{k\}) (ak，bk) ，假设所需要的系数能够使滤波输入  p p p 与输出  q q q 的差异最小，那么在窗口  w k w\_\{k\} wk 中的代价函数可以表示成式子  ( 3 ) (3) (3) 中的形式。  
 E ( a k , b k ) = ∑ i ∈ w k ( ( a k I i + b k − p i ) 2 + ϵ a k 2 )                 ( 3 ) E(a\_\{k\},b\_\{k\})=\\sum\_\{i\\in w\_\{k\}\}((a\_\{k\}I\_\{i\}+b\_\{k\}-p\_i\{\})^\{2\}+\\epsilon a\_\{k\}^\{2\})\\, \\, \\, \\, \\, \\, \\,\\: \\: \\: \\: \\: \\: \\: \\: (3) E(ak,bk)=i∈wk∑((akIi\+bk−pi)2\+ϵak2)(3)  
  其中  ϵ \\epsilon ϵ 是  L 2 L2 L2 范数正则化系数，他的作用是防止  a k a\_\{k\} ak 过大，（**他的物理意义会在后面讲到**）。该式子将  ( a k ， b k ) (a\_\{k\}，b\_\{k\}) (ak，bk) 的求解变成了线性回归的问题，根据最小二乘法可以很快得出我们想要的解（需要回顾最小二乘法的可以点[这里][Link 2]），在窗口  w k w\_\{k\} wk 中两者的解表示为下面的公式。  
 a k = 1 ∣ w ∣ ∑ i ∈ w k I i p i − μ k p k σ k 2 + ϵ                 ( 4 ) a\_\{k\}=\\frac\{\\frac\{1\}\{\\left | w \\right |\}\\sum\_\{i\\in w\_\{k\}\} I\_\{i\}p\_\{i\}-\\mu\_\{k\}p\_\{k\}\}\{\\sigma \_\{k\}^\{2\}+\\epsilon \}\\, \\, \\, \\, \\, \\, \\,\\: \\: \\: \\: \\: \\: \\: \\: (4) ak=σk2\+ϵ∣w∣1∑i∈wkIipi−μkpk(4)  
 b k = p k ˉ − a k μ k                 ( 5 ) b\_\{k\}=\\bar\{p\_\{k\}\}-a\_\{k\}\\mu \_\{k\}\\, \\, \\, \\, \\, \\, \\,\\: \\: \\: \\: \\: \\: \\: \\: (5) bk=pkˉ−akμk(5)

在公式中， μ k \\mu \_\{k\} μk 与  σ k \\sigma \_\{k\} σk 表示  I i I\_\{i\} Ii 在窗口内的均值、标准差，  ∣ w ∣ \{\\left | w \\right |\} ∣w∣ 表示窗口内像素块的总数， p k ˉ \\bar\{p\_\{k\}\} pkˉ 表示窗口内输入图像  p p p 的均值。而由于窗口是有尺寸的，输出图像中的一个像素点  q i q\_\{i\} qi 可以被不同的窗口中的线性系数求得，且不同窗口得到的输出值不同，因此对这些值求均值，最终得到的  q i q\_\{i\} qi 如式子  ( 6 ) (6) (6) 。  
 q i = a k ˉ I i + b k ˉ                 ( 6 ) q\_\{i\}=\\bar\{a\_\{k\}\}I\_\{i\}+\\bar\{b\_\{k\}\}\\, \\, \\, \\, \\, \\, \\,\\: \\: \\: \\: \\: \\: \\: \\: (6) qi=akˉIi\+bkˉ(6)  
  通过上面的公式，可以看出  ( a k ， b k ) (a\_\{k\}，b\_\{k\}) (ak，bk) 在空间上也是有变化的，所以输出图像  q q q 不是引导图像  I I I 的简单的尺度缩放。这个很好理解，因为  ( a k ， b k ) (a\_\{k\}，b\_\{k\}) (ak，bk) 的值也取决于输入图像  p p p ，但是由于公式  ( 6 ) (6) (6) 中对系数取了平均值，因此 ( a k ， b k ) (a\_\{k\}，b\_\{k\}) (ak，bk) 的梯度变化是明显小于  I I I 的，所以  q q q 的像素值变化主要还是取决于  I I I ，因此依然可以近似地认为 ▽ q ≈ a ▽ I \\triangledown q\\approx a\\triangledown I ▽q≈a▽I 。

经过上面的推导，可以将算法流程描述出来；下图中  r r r 表示窗口的边长， f m e a n ( I , r ) f\_\{mean\}\\left (I,r \\right ) fmean(I,r) 表示在大小为  ( r , r ) (r,r) (r,r) 窗口的内对图像做均值滤波。  
![在这里插入图片描述][20190818090102695.PNG_pic_center]  事实上，上面的公式还可以继续推导，将  ( a k ， b k ) (a\_\{k\}，b\_\{k\}) (ak，bk) 表示成下面的形式。  
 a k = ∑ j A k j ( I ) p j                 ( 7 ) a\_\{k\}=\\sum\_\{j\}A\_\{kj\}(I)p\_\{j\}\\, \\, \\, \\, \\, \\, \\,\\: \\: \\: \\: \\: \\: \\: \\: (7) ak=j∑Akj(I)pj(7)  
 b k = ∑ j B k j ( I ) p j                 ( 8 ) b\_\{k\}=\\sum\_\{j\}B\_\{kj\}(I)p\_\{j\}\\, \\, \\, \\, \\, \\, \\,\\: \\: \\: \\: \\: \\: \\: \\: (8) bk=j∑Bkj(I)pj(8)  
  其中  A k j A\_\{kj\} Akj 与  B k j B\_\{kj\} Bkj 仅与  I I I 有关，从而可以得到形式如公式  ( 1 ) (1) (1) 中的式子。（这个公式比较复杂，我就不写出来了，直接上截图）  
![在这里插入图片描述][20190818092027250.PNG_pic_center]

## 参数分析 ##

我们可以从双边滤波出发来思考引导滤波，在双边滤波中权重  W W W 由空间域与色彩域共同决定，距离越近的像素点的贡献越大，色彩域的贡献与颜色的相似性正相关，其实引导滤波也同样如此。引导滤波中空间域的贡献自然取决于窗口的大小，即由参数  r r r 决定。而标准差则是评判颜色差异性的参数，窗口中标准差越大，说明局部的像素相似性越差。

当输入图像作为自身的引导图时，很显然，当  ϵ = 0 \\epsilon = 0 ϵ=0 ， a = 1 a=1 a=1、 b = 0 b=0 b=0 ，输入输出之间不存在变化。而当  ϵ > 0 \\epsilon > 0 ϵ>0 时，考虑两种极端情况，第一种是引导图像完整平整，即在窗口内  I I I 为固定常数，那么 a = 0 a=0 a=0、 b = p k ˉ b=\\bar\{p\_\{k\}\} b=pkˉ，此时模型退化成为均值滤波。第二种情况是引导图像极其不平整，即在窗口内  I I I 的标准差接近无穷大，那么  a ≈ 1 a\\approx1 a≈1、 b ≈ 0 b\\approx0 b≈0 那么输出与输入几乎无变化。而图像的平整与否由方差评价，方差的大小由参数  ϵ \\epsilon ϵ 评判，当  ϵ \\epsilon ϵ 越大，引导图中被判定为“平整”区域的越多，输出图像约平滑，反之输出图像的平滑效果越差。

再反观作者提出公式  ( 3 ) (3) (3) 的假设，如果仅仅是输入图像与输出图像之间的差异最小，公式中的正则化项便可以省去，此时公式的意义就是输入图像与输出图像完全一致。而作者添加  ϵ a k 2 \\epsilon a\_\{k\}^\{2\} ϵak2 这一项，其实让输入图像与输入图像之前能够存在差别，我们姑且认为这种差别是噪声导致的，差异性的大小由  ϵ \\epsilon ϵ 决定，所以这其实也是证明了当  ϵ \\epsilon ϵ 越大，输出图像越平滑。

## 代码 ##

经过上述的推导与分析后，使用C++、Python分别利用OpenCV库实现了上述过程，其中求图像的均值滤波图采用OpneCV的boxfliter函数。在这里贴出C++中的主要函数(只需要自己编写一个主函数即可运行)，以及完整的Python脚本。完整的代码和测试图像，我放在了[这里][Link 3]。

### C++ ###

//-------------------------------------------------------------
    //作者：不用先生，2019.8.11
    //自实现的导向滤波去噪算法
    //guidedFilter.cpp
    //-------------------------------------------------------------
    
         
    //--------------------------------------------------------------
    //函数名：my_guidedFilter_oneChannel
    //函数功能：用于单通道图像（灰度图）的引导滤波函数；
    //参数：Mat &srcImg：输入图像，为单通道图像；
    //参数：Mat &guideImg：引导图像，为单通道图像，尺寸与输入图像一致；
    //参数：Mat &dstImg：输出图像，尺寸、通道数与输入图像吻合；
    //参数：const int rad：滤波器大小，应该保证为奇数，默认值为9；
    //参数：const double eps ：防止a过大的正则化参数，
    
    bool my_guidedFilter_oneChannel(Mat &srcImg, Mat &guideImg, Mat &dstImg, const int rad = 9, const double eps = 0.01)
    {
        
    	//--------------确保输入参数正确-----------
    	{
        
    		if (!srcImg.data || srcImg.channels() != 1)
    		{
        
    			cout << "输入图像错误，请重新输入图像" << endl;
    			return false;
    		}
    
    		if (!guideImg.data || guideImg.channels() != 1)
    		{
        
    			cout << "输入引导图像错误，请重新输入图像" << endl;
    			return false;
    		}
    
    		if (guideImg.cols != srcImg.cols || guideImg.rows != srcImg.rows)
    		{
        
    			cout << "输入图像与引导图像尺寸不匹配，请重新确认" << endl;
    			return false;
    		}
    		if (dstImg.cols != srcImg.cols || dstImg.rows != srcImg.rows || dstImg.channels() != 1)
    		{
        
    			cout << "参数输出图像错误，请重新确认" << endl;
    			return false;
    		}
    		if (rad % 2 != 1)
    		{
        
    			cout << "参数“rad”应为奇数，请修改" << endl;
    			return false;
    		}
    			
    	}
    
    	//--------------转换数值类型，并归一化-------------
    	srcImg.convertTo(srcImg, CV_32FC1, 1.0 / 255.0);
    	guideImg.convertTo(guideImg, CV_32FC1, 1.0 / 255.0);
    
    	//--------------求引导图像和输入图像的均值图
    	Mat mean_srcImg, mean_guideImg;
    	boxFilter(srcImg, mean_srcImg, CV_32FC1, Size(rad, rad));
    	boxFilter(guideImg, mean_guideImg, CV_32FC1, Size(rad, rad));
    
    	Mat mean_guideImg_square, mean_guideImg_srcImg;
    	boxFilter(guideImg.mul(guideImg), mean_guideImg_square, CV_32FC1, Size(rad, rad));
    	boxFilter(guideImg.mul(srcImg), mean_guideImg_srcImg, CV_32FC1, Size(rad, rad));
    
    	Mat var_guideImg = mean_guideImg_square - mean_guideImg.mul(mean_guideImg);
    	Mat cov_guideImg_srcImg = mean_guideImg_srcImg - mean_guideImg.mul(mean_srcImg);
    
    	Mat aImg = cov_guideImg_srcImg / (var_guideImg + eps);
    	Mat bImg = mean_srcImg - aImg.mul(mean_guideImg);
    
    	Mat mean_aImg, mean_bImg;
    	boxFilter(aImg, mean_aImg, CV_32FC1, Size(rad, rad));
    	boxFilter(bImg, mean_bImg, CV_32FC1, Size(rad, rad));
    
    	dstImg = (mean_aImg.mul(guideImg) + mean_bImg);
    
    	dstImg.convertTo(dstImg, CV_8UC1, 255);
    
    	return true;
    }
    
    //--------------------------------------------------------------
    //函数名：my_guidedFilter_threeChannel
    //函数功能：用于三通道图像（RGB彩色图）的引导滤波函数；
    //参数：Mat &srcImg：输入图像，为三通道图像；
    //参数：Mat &guideImg：引导图像，为三通道图像，尺寸与输入图像一致；
    //参数：Mat &dstImg：输出图像，尺寸、通道数与输入图像吻合；
    //参数：const int rad：滤波器大小，应该保证为奇数，默认值为9；
    //参数：const double eps ：防止a过大的正则化参数，
    
    bool my_guidedFilter_threeChannel(Mat &srcImg, Mat &guideImg, Mat &dstImg, const int rad = 9, const double eps = 0.01)
    {
        
    	//----------------确保输入参数正确-------------
    	{
        
    		if (!srcImg.data || srcImg.channels() != 3)
    		{
        
    			cout << "输入图像错误，请重新输入图像" << endl;
    			return false;
    		}
    
    		if (!guideImg.data || guideImg.channels() != 3)
    		{
        
    			cout << "输入引导图像错误，请重新输入图像" << endl;
    			return false;
    		}
    		if (guideImg.cols != srcImg.cols || guideImg.rows != srcImg.rows)
    		{
        
    			cout << "输入图像与引导图像尺寸不匹配，请重新确认" << endl;
    			return false;
    		}
    		if (rad % 2 != 1)
    		{
        
    			cout << "参数“rad”应为奇数，请修改" << endl;
    			return false;
    		}
    	
    	}
    
    	vector<Mat> src_vec, guide_vec,dst_vec;
    	split(srcImg, src_vec);
    	split(guideImg, guide_vec);
    
    	for (int i = 0; i < 3; i++)
    	{
        
    		Mat tempImg = Mat::zeros(srcImg.rows, srcImg.cols, CV_8UC1);
    		/*Mat tempImg;*/
    		my_guidedFilter_oneChannel(src_vec[i], guide_vec[i], tempImg, rad, eps);
    		dst_vec.push_back(tempImg);
    	}
    
    	merge(dst_vec, dstImg);
    
    	return true;
    }

### Python ###

# -*- coding: utf-8 -*-
    """
    Created on Sat Aug 17 18:46:20 2019
    
    @author: 不用先生
    """
    
    import cv2
    import numpy as np
    
    input_fn='03.jpg'
    
    def my_guidedFilter_oneChannel(srcImg,guidedImg,rad=9,eps=0.01):
        
        srcImg=srcImg/255.0
        guidedImg=guidedImg/255.0
        img_shape=np.shape(srcImg)
        
    #    dstImg=np.zeros(img_shape,dtype=float)
    #    
    #    P_mean=np.zeros(img_shape,dtype=float)
    #    I_mean=np.zeros(img_shape,dtype=float)
    #    I_square_mean=np.zeros(img_shape,dtype=float)
    #    I_mul_P_mean=np.zeros(img_shape,dtype=float)
    #    var_I=np.zeros(img_shape,dtype=float)
    #    cov_I_P=np.zeros(img_shape,dtype=float)
    #    
    #    a=np.zeros(img_shape,dtype=float)
    #    b=np.zeros(img_shape,dtype=float)
    #    a_mean=np.zeros(img_shape,dtype=float)
    #    b_mean=np.zeros(img_shape,dtype=float)
        
        P_mean=cv2.boxFilter(srcImg, -1, (rad, rad), normalize=True) 
        I_mean=cv2.boxFilter(guidedImg,-1, (rad, rad), normalize=True) 
        
        I_square_mean=cv2.boxFilter(np.multiply(guidedImg,guidedImg), -1, (rad, rad), normalize=True) 
        I_mul_P_mean=cv2.boxFilter(np.multiply(srcImg,guidedImg), -1, (rad, rad), normalize=True)
        
        var_I=I_square_mean-np.multiply(I_mean,I_mean)
        cov_I_P=I_mul_P_mean-np.multiply(I_mean,P_mean)
        
        a=cov_I_P/(var_I+eps)
        b=P_mean-np.multiply(a,I_mean)
        
        a_mean=cv2.boxFilter(a, -1, (rad, rad), normalize=True) 
        b_mean=cv2.boxFilter(b, -1, (rad, rad), normalize=True) 
        
        dstImg=np.multiply(a_mean,guidedImg)+b_mean
        
        return dstImg*255.0
        
    
    def my_guidedFilter_threeChannel(srcImg,guidedImg,rad=9,eps=0.01):
        
        img_shape=np.shape(srcImg)
    
        dstImg=np.zeros(img_shape,dtype=float)
    
        for ind in range(0,img_shape[2]):
            dstImg[:,:,ind]=my_guidedFilter_oneChannel(srcImg[:,:,ind],
                  guidedImg[:,:,ind],rad,eps)
        
        dstImg=dstImg.astype(np.uint8)
        
        return dstImg
    
    
    def main():
        img=cv2.imread(input_fn)
        print(np.shape(img))
    
        dstimg=my_guidedFilter_threeChannel(img,img,9,0.01)
        print(np.shape(dstimg))
    #    cv2.imwrite('output.jpg',dstimg)
        cv2.imshow('output',dstimg)
        cv2.waitKey(0)
        
    if __name__ == '__main__':
        main()

## 测试结果 ##

测试中选取了不同的参数  r r r 与  ϵ \\epsilon ϵ 形成对照，从结果中可以看出，当  r r r 与  ϵ \\epsilon ϵ 越大，图像被平滑的程度越大。  
![在这里插入图片描述][20190818101702277.png_pic_center]  
![在这里插入图片描述][20190818102018819.png_pic_center]

## 最后 ##

经典的双边滤波的时间复杂度为  O ( N r 2 ) O (Nr^\{2\}) O(Nr2) ，而双边滤波的时间复杂度为 O ( N ) O (N) O(N)，虽然后来有人对双边滤波进行了加速，但是却牺牲了精度。作者在文中对比了两者的效果，发现在细节保持上引导滤波要优于双边滤波，计算耗时也更短很多。后来作者又对引导滤波进行了加速，即 Fast Guided Filter，时间复杂度更低。

## 参考 ##

> He K, Sun J, Tang X. Guided image filtering\[J\]. IEEE transactions on pattern analysis and machine intelligence, 2012, 35(6): 1397-1409.
> 
> He K, Sun J. Fast guided filter\[J\]. arXiv preprint arXiv:1505.00996, 2015.

[https_blog.csdn.net_u013921430]: https://blog.csdn.net/u013921430
[Link 1]: https://blog.csdn.net/u013921430/article/details/84532068
[20190818080123999.PNG_pic_center]: /images/20211029/7df3b044c00440f6be0935b8e1bd5b0a.png
[Link 2]: https://blog.csdn.net/u013921430/article/details/84344510
[20190818090102695.PNG_pic_center]: /images/20211029/8ae6124b9b26451295abb686ce649ac1.png
[20190818092027250.PNG_pic_center]: /images/20211029/8625622c3a8146f88d30b4e666b773ed.png
[Link 3]: https://download.csdn.net/download/u013921430/11576121
[20190818101702277.png_pic_center]: /images/20211029/270b96e6aa82426d81c93e3b8b2b07ab.png
[20190818102018819.png_pic_center]: /images/20211029/c97738dde5c041e69d9fce0bc1187dd2.png