引导滤波推导及实现

最新推荐文章于 2023-02-11 15:35:30 发布

wx_幼儿园的学霸

最新推荐文章于 2023-02-11 15:35:30 发布

阅读量3.3k

点赞数 11

分类专栏： C++ OpenCV 文章标签：图像处理计算机视觉算法

本文链接：https://blog.csdn.net/leonardohaig/article/details/118254849

版权

C++ 同时被 2 个专栏收录

63 篇文章

订阅专栏

OpenCV

44 篇文章

订阅专栏

引导滤波推导及实现

微信公众号：幼儿园的学霸

前言

引导滤波顾名思义，就是有选择（引导）性的滤波，其与我们经常提及的高斯滤波、双边滤波相比，它具有引导性，说具体点就是，它通过输入一副图像（矩阵）作为引导图，这样滤波器就知道什么地方是边缘，以此更好的保护边缘，最终达到在滤波的同时，保持边缘细节。

各向同性滤波对待噪声和边缘信息采取一致的态度，从而导致噪声被磨皮的同时，图像中边缘，纹理和细节也被抹去了

推导

介绍

引导滤波的思想用一张引导图像 $I$ 产生权重 $W$ ，对输入图像 $p$ 进行处理得到输出图像 $q$ ,使得输出图像的梯度和引导图相似，但是灰度(亮度)与输入图像相似，这个过程可以表示为公式(1)中的内容。

Note:梯度和引导图相似，而非纹理等和引导图相似,关注这里是梯度！！！，后面提到的梯度翻转

$q_{i}=\sum_{j} W_{i j}(I) \cdot p_{j} \tag{1}$
公式中 $q, I, p$ 分别表示输出图像、引导图像和输入图像. $i, j$ 分别表示图像中像素点的索引。从公式(1)中可以看到权重 $W$ 仅与引导图像 $I$ 有关.而在双边滤波中，权重由输入图像自身决定。

在引导滤波的定义中，用到了局部线性模型。
该模型认为，某函数上一点与其邻近部分的点成线性关系，一个复杂的函数就可以用很多分段的局部的线性函数来表示，当需要求该函数上某一点的值时，只需计算所有包含该点的线性函数的值并做平均即可。可以暂时用下图简单的理解。

这种模型，在表示非解析函数上，非常有用。

引申: 局部加权线性回归

引导滤波中核心假设就是输出图像与引导图像在局部上是一个线性的模型,在局部窗口 $w_k$ 中，输入图像与引导图像的线性关系可以表示成公式(2)中的形式:
$q_{i}=a_{k} I_{i}+b_{k}, \forall i \in w_{k} \tag{2}$

公式中 $a_k,b_k$ 是窗口中的线性系数，在每个窗口中是常数。可以看到在滤波窗口 $w_k$ 中的所有像素点 $i$ ,其输出 $q_i$ 均是引导图 $I_i$ 的线性变换，这个变换由 $a_k,b_k$ 决定。我们对上式两边求梯度，得到:
$\nabla q=a \nabla I \tag{3}$
可以看到:
1.输出图像 $q$ 受到引导图像 $I$ 的约束，当引导图像在局部有梯度变化时，输出图像对应地也会有梯度变化。
2. $a_k$ 决定了梯度保留能力，如果 $a_k$ 较大，则梯度保留效果好;反之，则平滑效果好。这里也就明白了， $a_k,b_k$ 就是自适应调节因子了，我们希望在边缘处 $a_k$ 大些，其他地方 $a_k$ 小些，那么 $a_k,b_k$ 究竟应该怎样取值，才能实现这样一个自适应调节作用呢.

推导

因为我们是要对原图像进行滤波，所以在满足上面局部线性的基础上，还需要使输出图像尽可能逼近原图像,使得输入和输出之间的差异最小，原论文作者采用的是最小二乘方法，设计了一个代价函数使输出图像 $q$ 与原图像 $p$ 之间的均方误差最小，在窗口 $w_k$ 中的代价函数公式如下:
$E(a_{k}, b_{k})=\sum_{i \in w_{k}}\left(\left(a_{k} I_{i}+b_{k}-p_{i}\right)^{2}+\epsilon a_{k}^{2}\right) \tag{4}$

可以看到，除了均方误差项之外，作者还添加了含ϵ的一项。作者在原文中解释说这是一个正则化参数，用来防止 $a_k$ 过大。我想大部分人应该对这个很熟悉，这不就是深度学习里损失函数的正则项嘛，只不过那个是用来惩罚网络中的权重，防止过拟合罢了。另外我们试想一下，如果没有这一项，在引导图像等于输入图像时，那么最小化上式，一定存在 $a_k=1,b_k=0$ ，因为没有谁比自身更接近自身了嘛，这样带来的结果就是没有任何平滑效果，输出图像就是原图像本身。所以ϵ是一个非常重要的参数，选取不当会影响输出图像，后面会讲如何选取。

添加ϵ的项被称为L2正则项，添加该项后，上面方程为岭回归(Ridge Regression)方程

代价函数公式(4)使得 $a_k,b_k)$ 的求解变成了线性回归问题，现在根据最小二乘法来推导如何求解 $a_k,b_k)$

代价函数分别对 $a_k,b_k$ 求偏导:
$\begin{aligned} &\frac{\delta E}{\delta a_{k}}=\sum_{i \in w_{k}}\left(2\left(a_{k} I_{i}+b_{k}-p_{i}\right) I_{i}+2 \epsilon a_{k}\right) \\ &\frac{\delta E}{\delta b_{k}}=\sum_{i \in w_{k}} 2\left(a_{k} I_{i}+b_{k}-p_{i}\right) \end{aligned} \tag{5}$
设 $N$ 为窗口 $w_k$ 内像素点的数量，根据极小值处偏导数为0:
$\begin{aligned} &\frac{\delta E}{\delta b_{k}} =0 \\ &\Rightarrow a_{k} \sum_{i \in w_{k}} I_{i}-\sum_{i \in w_{k}} p_{i}+\sum_{i \in w_{k}} b_{k}=0 \\ &\Rightarrow N b_{k}=\sum_{i \in w_{k}} p_{i}-a_{k} \sum_{i \in w_{k}} I_{i} \\ &\Rightarrow b_{k} =\frac{\sum_{i \in w_{k}} p_{i}-a_{k} \sum_{i \in w_{k}} I_{i}}{N} = \bar{p}_{k}-a_{k} \mu_{k} \end{aligned} \tag{6}$
其中: $\mu_k$ 为窗口 $w_k$ 范围内引导图 $I$ 的均值; $\bar{p}_k$ 为窗口 $w_k$ 范围内输入图 $p$ 的均值.

接下来求 $a_k$
$\begin{aligned} &\frac{\delta E}{\delta a_{k}} = 0 \\ &\Rightarrow \sum_{i \in w_{k}}\left(a_{k} I_{i}^{2}+b_{k} I_{i}-p_{i} I_{i}+\epsilon a_{k}\right)=0 \\ &\Rightarrow \sum_{i \in w_{k}}\left(a_{k} I_{i}^{2}+ (\bar{p}_{k}-a_{k} \mu_{k}) I_{i}-p_{i} I_{i}+\epsilon a_{k}\right)=0 \\ \Rightarrow a_k &= \frac{\sum_{i \in w_{k}}(p_i I_i) - \sum_{i \in w_{k}}(\bar{p}_{k} I_i) }{ \sum_{i \in w_{k}}(I_{i}{^2}-\mu_k I_i + \epsilon ) } \\ &= \frac{\sum_{i \in w_{k}}(p_i *I_i) - \sum_{i \in w_{k}}(\bar{p}_{k} *I_i) }{ \sum_{i \in w_{k}}I_{i}{^2}- \mu_k*\sum_{i \in w_{k}} I_i + \sum_{i \in w_{k}} \epsilon } \\ &= \frac{\sum_{i \in w_{k}}(p_i *I_i) - (\frac{1}{N}\sum_{i \in w_{k}}p_i) (\sum_{i \in w_{k}}I_i )}{ \sum_{i \in w_{k}}I_{i}{^2}- (\frac{1}{N} \sum_{i \in w_{k}}I_i)*\left( \sum_{i \in w_{k}} I_i \right)+ N \epsilon } \\ &= \frac{\frac{1}{N} \sum_{i \in w_{k}}(p_i *I_i) - (\frac{1}{N}\sum_{i \in w_{k}}p_i) (\frac{1}{N}\sum_{i \in w_{k}}I_i )}{ \frac{1}{N}\sum_{i \in w_{k}}I_{i}{^2}- (\frac{1}{N} \sum_{i \in w_{k}}I_i)*\left( \frac{1}{N} \sum_{i \in w_{k}} I_i \right)+ \epsilon } \end{aligned} \tag{7}$
根据方差、协方差公式，以及方差和期望关系公式:
$\begin{aligned} D(X) &=E\left(X-E(X)\right)^{2}=E\left(X^{2}-2XE(X)+E^2(X)\right) \\ &=E(X^2)-2E(X)E(X)+E^2(X) \\ &=E(X^2)-E^2(X) \end{aligned} \tag{8}$

$\begin{aligned} \operatorname{Cov}(X, Y) &=E[(X-E(X)) \cdot(Y-E(Y))] \\ &=E(X Y)-E[X E(Y)+Y E(X)]+E[E(X) E(Y)] \\ &=E(X Y)-2 E(X) E(Y)+E(X) E(Y) \\ &=E(X Y)-E(X) E(Y) \end{aligned} \tag{9}$
当 $X = Y$ 时，
$\begin{aligned} \operatorname{Cov}(X, X) &=E(X X)-E(X) E(X) \\ &=D(X) \end{aligned} \tag{10}$
根据公式(8)(9)对公式(7)进行整理，可得:
$a_k = \frac{\operatorname{Cov}(p, I)}{D(I)+\epsilon} = \frac{\operatorname{Cov}(p, I)}{\sigma_{k}^{2}+\epsilon} \tag{11}$

得到上述 $a_k,b_k$ 后，可以对每个窗口 $w_k$ 计算一个 $a_k,b_k)$ ，但是，同一个像素会被不同的滤波窗口所包含。如下图所示:

上图中存在黑色和红色两个3x3的滤波窗口，对于两个窗口交叠区域的那个像素q，假设在黑色和红色窗中，求得的 $a_k,b_k$ 分别为 $a_{k1},b_{k1}),((a_{k2},b_{k2}))$ ;
那么在黑色窗中，它的输出可以表示为 $q_1=a_{k1}I_{i1}+b_{k1}$ ,在红色窗中，它的输出可以表示为 $q_2=a_{k2}I_{i2}+b_{k2}$ 。可见，对每个像素，都能计算出多个不同的 $a_k,b_k)$ ，因此作者使用各窗口 $a_k,b_k)$ 计算得到的 $q_i$ 值求平均得到输出的 $q_i$ 值，上述过程描述如下:
$\begin{aligned} q_{i} &=\frac{1}{N} \sum_{k, i \in w_{k}}\left(a_{k} I_{i}+b_{k}\right) \\ &=\bar{a}_{i} I_{i}+\bar{b}_{i} \end{aligned} \tag{12}$

其中, $\bar{a}_{i},\bar{b}_{i}$ 表示窗口 $w_k$ 范围内所有像素计算得到的 $a_k,b_k$ 的均值。（理论上，此处，不做均值处理也是可以的）
需要注意的是，经过公式(12)的变换后，输出图像不再是引导图像的简单的尺度缩放， $\nabla q$ 也不再是 $\nabla I$ 的简单缩放了，因为 $(\bar{a}_{i},\bar{b}_{i})$ 空间可变，因此相当于被一个均值滤波器滤波，是一个平均值，所以输出图像的梯度将会比输入图像的梯度小一些。在这种情况下， $\nabla q \approx \bar{a} \nabla I$ 。

引导滤波的算法流程及实现

计算流程

引导滤波的流程较为简单，按照上面的推导的公式进行计算即可:
0.输入引导图I和待处理图像p
1.计算引导图I的均值与方差、待处理图像p的均值，以及引导图和输入图像的乘积I*p
2.计算 $a_k,b_k)$
$\begin{aligned} &a_k=\left((I*p)_{mean}-I_{\text {mean }}* p_{\text {mean }}\right) /\left(Var(I)+\epsilon\right) \\ &b_k=p_{\text {mean }}-a_k *I_{\text {mean }} \end{aligned}$
3.计算a,b均值 $a_{mean},b_{mean}$
4.根据a,b均值求解导向滤波结果
$q=a_{mean}*I+b_{mean}$
引导滤波流程

快速引导滤波流程

由于引导滤波效率问题，何凯明博士在2015年，对其做了优化，主要是采用了图像金字塔思想，基本原理是将引导图，输入图都进行下采样计算 $a_k,b_k)$ ，然后 $a_k,b_k)$ 进行上采样恢复原始大小，整个算法流程如下：
快速引导滤波流程

自定义实现及效果

引导滤波算法快的关键之处在于计算 $a_k，b_k$ 所需的均值与方差都可以通过盒式滤波去实现，盒式滤波是一种与滤波半径大小无关的算法，其实现方法类似积分图思想，都是一种以空间换取时间的方式。
一般实现:

/*!
 * 自定义引导滤波
 * @param I 引导图CV_8UC1
 * @param p 待滤波图像CV_8UC1
 * @param dst 输出结果
 * @param r 滤波半径
 * @param eps 惩罚因子
 * @author liheng
 */
void guidedFilter(const cv::Mat& _I, const cv::Mat& _p, cv::Mat& dst,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);
    cv::Mat p;
    _p.convertTo(p, CV_64FC1);

    r = 2*r+1;

    //mean_I = boxfilter(I, r) ./ N;
    cv::Mat mean_I;
    cv::boxFilter(I, mean_I, CV_64FC1, cv::Size(r,r));//计算导向均值mean_I
    //mean_p = boxfilter(p, r) ./ N;
    cv::Mat mean_p;
    cv::boxFilter(p, mean_p, CV_64FC1, cv::Size(r,r));//计算输入图像均值mean_p
    //mean_Ip = boxfilter(I.*p, r) ./ N;
    cv::Mat mean_Ip;
    cv::boxFilter(I.mul(p), mean_Ip, CV_64FC1,
                  cv::Size(r, r));//互相关均值mean_Ip
    //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);//协方差Cov(I,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));//自相关均值mean_II
    //var_I = mean_II - mean_I .* mean_I;
    cv::Mat var_I = mean_II - mean_I.mul(mean_I);//方差D(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_b
    cv::Mat mean_b;
    cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(r,r));

    //TODO:感兴趣的可以将mean_a,mean_b 显示出来,观察其结果
    //在图中可以看到，在边缘部分或变化剧烈的部分，a的值接近于1（白色），b的值接近为0（黑色），
    // 而在变化平坦的区域，a的值接近0（黑色），b的值为平坦区域像素的均值。
    //cv::imshow("a",mean_a);
    //cv::imshow("b",mean_b);
    //cv::waitKey(0);


    //q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
    cv::Mat q = mean_a.mul(I) + mean_b;
    q.convertTo(dst,CV_8UC1);
}

快速引导滤波实现:

/*!
 * 自定义实现快速引导滤波
 * @param I_org 引导图CV_8UC1
 * @param p_org 待滤波图像CV_8UC1
 * @param dst 输出结果
 * @param r 滤波半径
 * @param eps 惩罚因子
 * @param s 下采样倍数,为1时,即为上面的普通引导滤波
 * @author liheng
 */
void fastGuidedFilter(const cv::Mat& I_org, const cv::Mat& p_org, cv::Mat& dst,int r, double eps, int s)
{
    /*
    % GUIDEDFILTER   O(N) 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;
    I_org.convertTo(_I, CV_64FC1);

    resize(_I,I,cv::Size(),1.0/s,1.0/s);



    cv::Mat p,_p;
    p_org.convertTo(_p, CV_64FC1);
    //p = _p;
    resize(_p, p, cv::Size(),1.0/s,1.0/s);


    r = (2 * r + 1)/s+1;

    //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));
    cv::Mat rmean_a;
    resize(mean_a, rmean_a, cv::Size(I_org.cols, I_org.rows));

    //mean_b = boxfilter(b, r) ./ N;
    cv::Mat mean_b;
    cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(r, r));
    cv::Mat rmean_b;
    resize(mean_b, rmean_b, cv::Size(I_org.cols, I_org.rows));

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

调用实例及结果:

#include <iostream>
#include <opencv2/ximgproc/edge_filter.hpp> //opencv自带引导滤波实现
#include <opencv2/opencv.hpp>

int main()
{
    cv::Mat srcImage = cv::imread("./guidefilter.bmp", cv::IMREAD_COLOR);
    cv::imshow("src", srcImage);

    //cv::Mat I;
    //cv::cvtColor(srcImage,I,cv::COLOR_BGR2GRAY);

    std::vector<cv::Mat> channels;
    cv::split(srcImage,channels);

    std::vector<cv::Mat> my_guideChannels(3),cv_guideChannels(3);


    for(int i=0;i<3;++i)
    {
        guidedFilter(channels[i],channels[i],my_guideChannels[i],9,30);
        //fastGuidedFilter(channels[i],channels[i],my_guideChannels[i],9,30,3);
        cv::ximgproc::guidedFilter(channels[i],channels[i],cv_guideChannels[i],9,30);
    }


    cv::Mat my_dst,cv_dst;
    cv::merge(my_guideChannels,my_dst);
    cv::imshow("my_dst",my_dst);

    cv::merge(cv_guideChannels,cv_dst);
    cv::imshow("cv_dst",cv_dst);

    //std::cout<<my_dst-cv_dst<<std::endl;
    cv::waitKey(0);


    return 0;
}

程序输入图像及运行结果如下:

说明	图像
输入图像
自定义实现
OpenCV自带实现

从图像上看，自定义实现和OpenCV自带实现效果基本一致，如果将2图像进行减操作，相减结果基本均为0
当然，OpenCV自带实现的算法由于进行一定的加速，其运行效率较高。

引导滤波的应用

应用

1.保边滤波
当引导图为输入图像，即 $I = p$ 时，该算法成为一个保边滤波器，此时可得上述计算 $a_k,b_k)$ 的公式变为:
$\begin{gathered} a_{k}=\frac{\sigma_{k}^{2}}{\sigma_{k}^{2}+\epsilon} \\ b_{k}=\left(1-a_{k}\right) \bar{p}_{k} \end{gathered}$
考虑以下两种情况：

Case 1：平坦区域。如果在某个滤波窗口内，该区域是相对平滑的， $\sigma_{k}^{2}\ll \epsilon$ 时， $a_{k} \approx 0, b_{k} \approx \bar{p}_{k}$ 。相当于对该区域作均值滤波。
Case 2：高方差区域。相反，如果该区域是边缘区域，方差很大 $\sigma_{k}^{2}\gg \epsilon$ 时, $a_{k} \approx 1, b_{k} \approx 0$ 。相当于在区域保持原有梯度。
由上可知， $\epsilon$ 是调整图的模糊程度与边缘检测精度的参数
2.图像去雾
在图像去雾中，导向滤波一般用来细化透射率图像，以原图的灰度图为导向图，以粗投射率图为输入图，能得到非常精细的透射率图像。

优点

引导滤波相对于双边滤波最大的优点在于算法的复杂度与窗口的大小无关，对于处理较为大型的图片时，在效率上有明显的提升
同时，引导滤波可以很好地克服双边滤波中出现的梯度翻转(gradient reversal )的现象，因为假设前提便是存在线性关系，可以保证梯度一致；而对于双边滤波而言，对于梯度变化大的地方，由于周围没有相似的像素，高斯函数的权重不稳定，导致最终梯度出现反转现象