导读

在图像滤波算法中，导向滤波、双边滤波、最小二乘滤波并称三大保边滤波器，他们是各向异性滤波器。相对于常见的均值滤波、高斯滤波等各向同性滤波器，他们最大的特点是在去除噪声的同时，能最大限度保持边缘不被平滑。本文讲解导向滤波及其应用。
总的来讲，导向滤波就是尽可能让输出图像的梯度和导向图相似，同时让输出图像的灰度(亮度)与输入图像相似，以此来保留边缘并且滤除噪声。

原理推导

我们先看下图：

输入图像ppp，经过引导图像III, 滤波得到输出图像qqq, 导向滤波算法中有一个重要假设：即在局部窗口wkw_kwk上，导向图III和输出图qqq存在局部线性关系：
qi=akIi+bk,∀i∈wk(1)q_i = a_kI_i + b_k, \forall{i\in{w_k}}\tag{1}qi=akIi+bk,∀i∈wk(1)
同时在窗口wkw_kwk上, 滤波后的图像qqq和输入图像ppp有如下关系：
qi=pi−ni，∀i∈wk(2)q_i = p_i - n_i，\forall{i\in{w_k}}\tag{2}qi=pi−ni，∀i∈wk(2)
这样公式(1)(1)(1)的线性关系保证了如果在每个局部窗口wkw_kwk中，如果导向图III中存在一个边缘，输出图像qqq将保持边缘不变。同时，滤波结果图qqq要尽可能与输入图像ppp相同以此减小滤波带来的信息损失，该算法的最小二乘表示即：
argmin∑i∈wk(qi−pi)2=argmin∑i∈wk(akIi+bk−pi)2argmin\sum_{i\in{w_k}}(q_i-p_i)^2 = argmin\sum_{i\in{w_k}}(a_kI_i+b_k-p_i)^2argmini∈wk∑(qi−pi)2=argmini∈wk∑(akIi+bk−pi)2
这是求解最优值的问题，引入一个正则化参数ϵ\epsilonϵ防止aka_kak过大，得到损失函数：
E(ak,bk)=∑i∈wk((akIi+bk−pi)2+ϵak2)(3)E(a_k,b_k)=\sum_{i\in{w_k}}((a_kI_i+b_k-p_i)^2+\epsilon{a_k^2})\tag{3} E(ak,bk)=i∈wk∑((akIi+bk−pi)2+ϵak2)(3)
运用最小二乘法求解极小值，利用极小值处导数为0，求解过程如下，∣w∣|w|∣w∣为窗口wkw_kwk像素数量：
δEak=∑i∈wk(2(akIi+bk−pi)Ii+2ϵak)=0⇒∑i∈wk(akIi2+bkIi−piIi+ϵak)=0δEbk=∑i∈wk2(akIi+bk−pi)=0⇒ak∑i∈wkIi−∑i∈wkpi+∑i∈wkbk=0⇒∣w∣bk=∑i∈wkpi−ak∑i∈wkIi\begin{aligned} \frac{\delta{E}}{a_k} &= \sum_{i\in{w_k}}(2(a_kI_i+b_k-p_i)I_i+2\epsilon{a_k}) = 0 \\ & \Rightarrow \sum_{i\in{w_k}}(a_kI_i^2+b_kI_i-p_iI_i+\epsilon{a_k}) = 0 \\ \frac{\delta{E}}{b_k} & = \sum_{i\in{w_k}}2(a_kI_i+b_k-p_i)= 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 |w|b_k = \sum_{i\in{w_k}}{p_i}-a_k\sum_{i\in{w_k}}{I_i} \end{aligned} akδEbkδE=i∈wk∑(2(akIi+bk−pi)Ii+2ϵak)=0⇒i∈wk∑(akIi2+bkIi−piIi+ϵak)=0=i∈wk∑2(akIi+bk−pi)=0⇒aki∈wk∑Ii−i∈wk∑pi+i∈wk∑bk=0⇒∣w∣bk=i∈wk∑pi−aki∈wk∑Ii
由上面推到得出：
ak=∑i∈wkpiIi−bk∑i∈wkIi∑i∈wk(Ii2+ϵ)bk=∑i∈wkpi−ak∑i∈wkIi∣w∣=p‾k−akμk(4)\begin{aligned} a_k &= \frac{\sum_{i\in{w_k}}p_iI_i-b_k\sum_{i\in{w_k}}I_i}{\sum_{i\in{w_k}}(I_i^2+\epsilon)} \\ b_k &= \frac{\sum_{i\in{w_k}}{p_i}-a_k\sum_{i\in{w_k}}{I_i}}{|w|} \\ & = \overline{p}_k - a_k\mu_k \tag{4} \end{aligned} akbk=∑i∈wk(Ii2+ϵ)∑i∈wkpiIi−bk∑i∈wkIi=∣w∣∑i∈wkpi−ak∑i∈wkIi=pk−akμk(4)
其中：μk——窗口wk范围内引导图I的均值；p‾k——窗口wk范围内输入图p的均值。\begin{aligned} & \mu_k——窗口w_k范围内引导图I的均值；\\ & \overline{p}_k——窗口w_k范围内输入图p的均值。 \end{aligned}μk——窗口wk范围内引导图I的均值；pk——窗口wk范围内输入图p的均值。
将bkb_kbk代入aka_kak计算可得：
ak=∑i∈wkpiIi−bk∑i∈wkIi∑i∈wk(Ii2+ϵ)a_k = \frac{\sum_{i\in{w_k}}p_iI_i-b_k\sum_{i\in{w_k}}I_i}{\sum_{i\in{w_k}}(I_i^2+\epsilon)} ak=∑i∈wk(Ii2+ϵ)∑i∈wkpiIi−bk∑i∈wkIi
⇒ak=∑i∈wk(piIi−p‾kIi)∑i∈wk(Ii2+μkIi+ϵ)\Rightarrow a_k = \frac{\sum_{i\in{w_k}}(p_iI_i-\overline{p}_kI_i)}{\sum_{i\in{w_k}}(I_i^2+\mu_kI_i+\epsilon)} ⇒ak=∑i∈wk(Ii2+μkIi+ϵ)∑i∈wk(piIi−pkIi)
⇒ak=∑i∈wk(piIi−p‾kIi+μkpi−μkp‾k)∑i∈wk(Ii2−μkIi−μkIi+μk2+ϵ)\Rightarrow a_k = \frac{\sum_{i\in{w_k}}(p_iI_i-\overline{p}_kI_i {\color{red}{+\mu_kp_i-\mu_k\overline{p}_k}})}{\sum_{i\in{w_k}}(I_i^2-\mu_kI_i{\color{red}{-\mu_kI_i+\mu_k^2}}+\epsilon)} ⇒ak=∑i∈wk(Ii2−μkIi−μkIi+μk2+ϵ)∑i∈wk(piIi−pkIi+μkpi−μkpk)
上下两边同除以∣w∣|w|∣w∣。得到：
ak=1∣w∣∑i∈wkpiIi−p‾kμk+μkp‾k−μkp‾k1∣w∣∑i∈wk(Ii−μk)2+ϵa_k = \frac{\frac{1}{|w|}\sum_{i\in{w_k}}p_iI_i-\overline{p}_k\mu_k{\color{red}{+\mu_k\overline{p}_k-\mu_k\overline{p}_k}}}{\frac{1}{|w|}\sum_{i\in{w_k}}(I_i-\mu_k)^2+\epsilon} ak=∣w∣1∑i∈wk(Ii−μk)2+ϵ∣w∣1∑i∈wkpiIi−pkμk+μkpk−μkpk
最后：
ak=1∣w∣∑i∈wkpiIi−μkp‾kσk2+ϵ(5)\color{red}{a_k = \frac{\frac{1}{|w|}\sum_{i\in{w_k}}p_iI_i{ - \mu_k\overline{p}_k}}{\sigma_k^2+\epsilon}} \tag{5}ak=σk2+ϵ∣w∣1∑i∈wkpiIi−μkpk(5)
bk=p‾k+akμk(6)\color{red}{b_k = \overline{p}_k+a_k\mu_k} \tag{6}bk=pk+akμk(6)
得到上述公式后，可以对每个窗口wkw_kwk计算一个(ak,bk)(a_k,b_k)(ak,bk)，但是，每个像素都被包含在多个窗口中，对每个像素，都能计算出多个(ak,bk)(a_k,b_k)(ak,bk)，我么将使用多个(ak,bk)(a_k,b_k)(ak,bk)计算得到的qiq_iqi值求平均得到输出qiq_iqi值，上述过程描述如下:
qi=1wk∑k,i∈wk(akIi+bk)=a‾iIi+b‾i\begin{aligned} q_i &= \frac{1}{w_k}\sum_{k,i\in{w_k}}(a_kI_i+b_k) \\ & = \overline{a}_iI_i+\overline{b}_i \end{aligned} qi=wk1k,i∈wk∑(akIi+bk)=aiIi+bi
其中：a‾i——窗口wk范围内所有像素计算得到的ak的均值；b‾i——窗口wk范围内所有像素计算得到的bk的均值。\begin{aligned} & \overline{a}_i——窗口w_k范围内所有像素计算得到的a_k的均值；\\ & \overline{b}_i——窗口w_k范围内所有像素计算得到的b_k的均值。 \end{aligned}ai——窗口wk范围内所有像素计算得到的ak的均值；bi——窗口wk范围内所有像素计算得到的bk的均值。

导向滤波的应用

保边滤波
当I=pI=pI=p时，该算法成为一个保边滤波器。上述(ak,bk)(a_k,b_k)(ak,bk)计算公式变化为:
ak=σk2σk2+ϵbk=(1−ak)p‾ka_k = \frac{\sigma_k^2}{\sigma_k^2+\epsilon} \\ b_k = (1-a_k)\overline{p}_kak=σk2+ϵσk2bk=(1−ak)pk
考虑以下两种情况：
- Case 1：平坦区域。如果在某个滤波窗口内，该区域是相对平滑的，方差σk2\sigma_k^2σk2将远远小于ϵ\epsilonϵ。从而ak≈0,bk≈p‾ka_k≈0,b_k≈\overline{p}_kak≈0,bk≈pk。相当于对该区域作均值滤波。
- Case 2：高方差区域。相反，如果该区域是边缘区域，方差很大，σk2\sigma_k^2σk2将远远大于ϵ\epsilonϵ。从而ak≈1,bk≈0a_k≈1,b_k≈0ak≈1,bk≈0。相当于在区域保持原有梯度。
  以上可以出:ϵ\epsilonϵ为界定平滑区域和边缘区域的阈值。
图像去雾
在图像去雾中，导向滤波一般用来细化透射率图像，以原图的灰度图为导向图，以粗投射率图为输入图，能得到非常精细的透射率图像。

当然，导向滤波的应用不止以上两种，网上还有图像融合等应用，本人没有去了解。

导向滤波的实现

导向滤波的代码实现较为简单，我们直接贴出计算流程

快速导向滤波的实现

由于导向滤波效率问题，何凯明博士在2015年，对其做了优化，基本原理是将导向图，输入图都进行下采样计算(ak,bk)(a_k,b_k)(ak,bk)，然后对(ak,bk)(a_k,b_k)(ak,bk)进行上采样恢复原始大小，整个算法流程如下：

算法效果

我们演示一下保边滤波效果：

代码

接下来，废话不多说，我们上代码：https://github.com/ZPEthanWen/ImageAlgorithmDraft，为防止github无法访问，我们直接贴上代码：

导向滤波

#include <opencv2/opencv.hpp>
//导向滤波
void GuidedFilter(cv::Mat& srcImage, cv::Mat& guidedImage, cv::Mat& outputImage, int filterSize, double eps)
{try{if (srcImage.empty() || guidedImage.empty() || filterSize <= 0 || eps < 0 ||srcImage.channels() != 1 || guidedImage.channels() != 1){throw "params input error";}cv::Mat srcImageP, srcImageI, meanP, meanI, meanIP, meanII, varII, alfa, beta;srcImage.convertTo(srcImageP, CV_32FC1);guidedImage.convertTo(srcImageI, CV_32FC1);cv::boxFilter(srcImageP, meanP, CV_32FC1, cv::Size(filterSize, filterSize));cv::boxFilter(srcImageI, meanI, CV_32FC1, cv::Size(filterSize, filterSize));cv::boxFilter(srcImageI.mul(srcImageP), meanIP, CV_32FC1, cv::Size(filterSize, filterSize));cv::boxFilter(srcImageI.mul(srcImageI), meanII, CV_32FC1, cv::Size(filterSize, filterSize));varII = meanII - meanI.mul(meanI); alfa = (meanIP - meanI.mul(meanP)) / (varII + eps);beta = meanP - alfa.mul(meanI);cv::boxFilter(alfa, alfa, CV_32FC1, cv::Size(filterSize, filterSize));cv::boxFilter(beta, beta, CV_32FC1, cv::Size(filterSize, filterSize));outputImage = (alfa.mul(srcImageI) + beta);}catch (cv::Exception& e){throw e;}catch (std::exception& e){throw e;}
}

快速导向滤波

#include <opencv2/opencv.hpp>
//快速导向滤波
void FastGuidedFilter(cv::Mat& srcImage, cv::Mat& guidedImage, cv::Mat& outputImage, int filterSize, double eps, int samplingRate)
{try{if (srcImage.empty() || guidedImage.empty() || filterSize <= 0 || eps < 0 ||srcImage.channels() != 1 || guidedImage.channels() != 1 || samplingRate < 1){throw "params input error";}cv::Mat srcImageP, srcImageSubI, srcImageI, meanP, meanI, meanIP, meanII, var, alfa, beta;cv::resize(srcImage, srcImageP, cv::Size(srcImage.cols / samplingRate, srcImage.rows / samplingRate));cv::resize(guidedImage, srcImageSubI, cv::Size(srcImage.cols / samplingRate, srcImage.rows / samplingRate));filterSize = filterSize / samplingRate;srcImageP.convertTo(srcImageP, CV_32FC1);guidedImage.convertTo(srcImageI, CV_32FC1);srcImageSubI.convertTo(srcImageSubI, CV_32FC1);cv::boxFilter(srcImageP, meanP, CV_32FC1, cv::Size(filterSize, filterSize));cv::boxFilter(srcImageSubI, meanI, CV_32FC1, cv::Size(filterSize, filterSize));cv::boxFilter(srcImageSubI.mul(srcImageP), meanIP, CV_32FC1, cv::Size(filterSize, filterSize));cv::boxFilter(srcImageSubI.mul(srcImageSubI), meanII, CV_32FC1, cv::Size(filterSize, filterSize));var = meanII - meanI.mul(meanI);alfa = (meanIP - meanI.mul(meanP)) / (var + eps);beta = meanP - alfa.mul(meanI);cv::boxFilter(alfa, alfa, CV_32FC1, cv::Size(filterSize, filterSize));cv::boxFilter(beta, beta, CV_32FC1, cv::Size(filterSize, filterSize));cv::resize(alfa, alfa, cv::Size(srcImage.cols, srcImage.rows));cv::resize(beta, beta, cv::Size(srcImage.cols, srcImage.rows));outputImage = alfa.mul(srcImageI) + beta;}catch (cv::Exception& e){throw e;}catch (std::exception& e){throw e;}
}

参考

[1] 视觉一只白 .《导向滤波的原理及实现》[DB/OL].
[2] lsflll.《导向滤波(Guided Filter)公式详解》[DB/OL]
[3] SongpingWang.《OpenCV—Python 导向滤波》[DB/OL]

详解——导向滤波(Guided Filter)和快速导向滤波相关推荐

快速引导滤波原理详解（fast guided filter）
符号标记 q : 输出图像 p : 输入图像 I : 引导图 a k , b k : 线性变换模型 r : 滤波半径 s : 下采样因子 q: 输出图像 \\[2ex] p: 输入图像 \\[2ex] ...
三种经典图像滤波方法介绍——双边滤波(Bilateral filter)、导向滤波(Guided Fliter)、滚动导向滤波(RollingGuidedFilter)
文章目录一.前言二.双边滤波(Bilateral filter) 2.1 双边滤波的理论介绍及公式推导 2.2 双边滤波的matlab程序实现三.导向滤波(Guided Fliter) 3.1 ...
cass地籍参数设置快捷命令_南方cass详解+视频教程+插件汇总，小白快速上手！限时领取...
南方cass详解+视频教程+插件汇总,小白快速上手!限时领取对于测量员来说,扎实的基本功尤为重要.测量员需要通过测量为工程建设提供数据和图纸,这样才能在工作中做到游刃有余. 今天给大家整理了南方Ca ...
springboot 详解 (四）redis filter
---------------------------------------------------------------------------------------------------- ...
粒子滤波 particle filter —从贝叶斯滤波到粒子滤波—Part-III（重要性采样序贯重要性采样SIS）
粒子滤波 particle filter -从贝叶斯滤波到粒子滤波-Part-III(重要性采样&序贯重要性采样SIS) 原创不易,路过的各位大佬请点个赞机动目标跟踪/非线性滤波/传感器融合 ...
导向滤波python_导向滤波(Guided Filter)简要介绍
1.介绍提到导向滤波,首先想到的是"何恺明",他的确是真大神,在图像领域,是中国人的骄傲,深度学习流行的时候,也是提出各种新算法,比如ResNets,而最近两年,深度学习的发展已 ...
引导滤波(guided filter)理解和代码实现
最近在学习图片的滤波和去噪的相关知识,查阅了一些资料参考了一些博客,这里做一个整合+理解.参考的博客资料在文末. 引入普通滤波的概念假设输入图像为p,滤波窗口为wk,经过滤波后的输出图像为q,那么q ...
双边滤波(bilateral filter)以及联合双边滤波（joint bilateral filter）
文章目录双边滤波理论公式代码(C++) 数学辅助理解联合双边滤波(joint bilateral filter) 参考链接写在最后双边滤波自用备忘,若侵则删. 理论公式利用二维高斯函数 ...
logback logback.xml常用配置详解（三） filter
<filter>: 过滤器,执行一个过滤器会有返回个枚举值,即DENY,NEUTRAL,ACCEPT其中之一.返回DENY,日志将立即被抛弃不再经过其他过滤器:返回NEUTRAL,有序列表 ...

详解——导向滤波(Guided Filter)和快速导向滤波

文章目录

导读