文章目录

• 导读
• 原理推导
• 导向滤波的应用
• 导向滤波的实现
• 快速导向滤波的实现
• 算法效果
• 代码
• 参考

导读

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

原理推导

我们先看下图：

输入图像 p p p，经过引导图像 I I I, 滤波得到输出图像 q q q, 导向滤波算法中有一个重要假设：即在局部窗口 w k w_k wk​上，导向图 I I I和输出图 q q q存在局部线性关系：
q i = a k I i + b k , ∀ i ∈ w k (1) q_i = a_kI_i + b_k, \forall{i\in{w_k}}\tag{1} qi​=ak​Ii​+bk​,∀i∈wk​(1)
同时在窗口 w k w_k wk​上, 滤波后的图像 q q q和输入图像 p p p有如下关系：
q i = p i − n i ， ∀ i ∈ w k (2) q_i = p_i – n_i，\forall{i\in{w_k}}\tag{2} qi​=pi​−ni​，∀i∈wk​(2)
这样公式 ( 1 ) (1) (1)的线性关系保证了如果在每个局部窗口 w k w_k wk​中，如果导向图 I I I中存在一个边缘，输出图像 q q q将保持边缘不变。同时，滤波结果图 q q q要尽可能与输入图像 p p p相同以此减小滤波带来的信息损失，该算法的最小二乘表示即：
a r g m i n ∑ i ∈ w k ( q i − p i ) 2 = a r g m i n ∑ i ∈ w k ( a k I i + b k − p i ) 2 argmin\sum_{i\in{w_k}}(q_i-p_i)^2 = argmin\sum_{i\in{w_k}}(a_kI_i+b_k-p_i)^2 argmini∈wk​∑​(qi​−pi​)2=argmini∈wk​∑​(ak​Ii​+bk​−pi​)2
这是求解最优值的问题，引入一个正则化参数 ϵ \epsilon ϵ防止 a k a_k ak​过大，得到损失函数：
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_kI_i+b_k-p_i)^2+\epsilon{a_k^2})\tag{3} E(ak​,bk​)=i∈wk​∑​((ak​Ii​+bk​−pi​)2+ϵak2​)(3)
运用最小二乘法求解极小值，利用极小值处导数为0，求解过程如下， ∣ w ∣ |w| ∣w∣为窗口 w k w_k wk​像素数量：
δ E a k = ∑ i ∈ w k ( 2 ( a k I i + b k − p i ) I i + 2 ϵ a k ) = 0 ⇒ ∑ i ∈ w k ( a k I i 2 + b k I i − p i I i + ϵ a k ) = 0 δ E b k = ∑ i ∈ w k 2 ( a k I i + b k − p i ) = 0 ⇒ a k ∑ i ∈ w k I i − ∑ i ∈ w k p i + ∑ i ∈ w k b k = 0 ⇒ ∣ w ∣ b k = ∑ i ∈ w k p i − a k ∑ i ∈ w k I i \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​δE​bk​δE​​=i∈wk​∑​(2(ak​Ii​+bk​−pi​)Ii​+2ϵak​)=0⇒i∈wk​∑​(ak​Ii2​+bk​Ii​−pi​Ii​+ϵak​)=0=i∈wk​∑​2(ak​Ii​+bk​−pi​)=0⇒ak​i∈wk​∑​Ii​−i∈wk​∑​pi​+i∈wk​∑​bk​=0⇒∣w∣bk​=i∈wk​∑​pi​−ak​i∈wk​∑​Ii​​
由上面推到得出：
a k = ∑ i ∈ w k p i I i − b k ∑ i ∈ w k I i ∑ i ∈ w k ( I i 2 + ϵ ) b k = ∑ i ∈ w k p i − a k ∑ i ∈ w k I i ∣ w ∣ = p ‾ k − a k μ 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} ak​bk​​=∑i∈wk​​(Ii2​+ϵ)∑i∈wk​​pi​Ii​−bk​∑i∈wk​​Ii​​=∣w∣∑i∈wk​​pi​−ak​∑i∈wk​​Ii​​=p​k​−ak​μk​​(4)
其中： μ k — — 窗 口 w k 范 围 内 引 导 图 I 的 均 值 ； p ‾ k — — 窗 口 w k 范 围 内 输 入 图 p 的 均 值 。 \begin{aligned} & \mu_k——窗口w_k范围内引导图I的均值；\\ & \overline{p}_k——窗口w_k范围内输入图p的均值。 \end{aligned} ​μk​——窗口wk​范围内引导图I的均值；p​k​——窗口wk​范围内输入图p的均值。​
将 b k b_k bk​代入 a k a_k ak​计算可得：
a k = ∑ i ∈ w k p i I i − b k ∑ i ∈ w k I i ∑ i ∈ w k ( I i 2 + ϵ ) 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∈wk​​pi​Ii​−bk​∑i∈wk​​Ii​​
⇒ a k = ∑ i ∈ w k ( p i I i − p ‾ k I i ) ∑ i ∈ w k ( I i 2 + μ k I i + ϵ ) \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​+μk​Ii​+ϵ)∑i∈wk​​(pi​Ii​−p​k​Ii​)​
⇒ a k = ∑ i ∈ w k ( p i I i − p ‾ k I i + μ k p i − μ k p ‾ k ) ∑ i ∈ w k ( I i 2 − μ k I i − μ k I i + μ k 2 + ϵ ) \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​−μk​Ii​−μk​Ii​+μk2​+ϵ)∑i∈wk​​(pi​Ii​−p​k​Ii​+μk​pi​−μk​p​k​)​
上下两边同除以 ∣ w ∣ |w| ∣w∣。得到：
a k = 1 ∣ w ∣ ∑ i ∈ w k p i I i − p ‾ k μ k + μ k p ‾ k − μ k p ‾ k 1 ∣ w ∣ ∑ i ∈ w k ( I i − μ 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∈wk​​pi​Ii​−p​k​μk​+μk​p​k​−μk​p​k​​
最后：
a k = 1 ∣ w ∣ ∑ i ∈ w k p i I i − μ k p ‾ k σ k 2 + ϵ (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∈wk​​pi​Ii​−μk​p​k​​(5)
b k = p ‾ k + a k μ k (6) \color{red}{b_k = \overline{p}_k+a_k\mu_k} \tag{6} bk​=p​k​+ak​μk​(6)
得到上述公式后，可以对每个窗口 w k w_k wk​计算一个 ( a k , b k ) (a_k,b_k) (ak​,bk​)，但是，每个像素都被包含在多个窗口中，对每个像素，都能计算出多个 ( a k , b k ) (a_k,b_k) (ak​,bk​)，我么将使用多个 ( a k , b k ) (a_k,b_k) (ak​,bk​)计算得到的 q i q_i qi​值求平均得到输出 q i q_i qi​值，上述过程描述如下:
q i = 1 w k ∑ k , i ∈ w k ( a k I i + b k ) = a ‾ i I i + 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​​=wk​1​k,i∈wk​∑​(ak​Ii​+bk​)=ai​Ii​+bi​​
其中： a ‾ i — — 窗 口 w k 范 围 内 所 有 像 素 计 算 得 到 的 a k 的 均 值 ； b ‾ i — — 窗 口 w k 范 围 内 所 有 像 素 计 算 得 到 的 b k 的 均 值 。 \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 = p I=p I=p时，该算法成为一个保边滤波器。上述 ( a k , b k ) (a_k,b_k) (ak​,bk​)计算公式变化为:
  a k = σ k 2 σ k 2 + ϵ b k = ( 1 − a k ) p ‾ k a_k = \frac{\sigma_k^2}{\sigma_k^2+\epsilon} \\ b_k = (1-a_k)\overline{p}_k ak​=σk2​+ϵσk2​​bk​=(1−ak​)p​k​
  考虑以下两种情况：
  • Case 1：平坦区域。如果在某个滤波窗口内，该区域是相对平滑的，方差 σ k 2 \sigma_k^2 σk2​将远远小于 ϵ \epsilon ϵ。从而 a k ≈ 0 , b k ≈ p ‾ k a_k≈0,b_k≈\overline{p}_k ak​≈0,bk​≈p​k​。相当于对该区域作均值滤波。
  • Case 2：高方差区域。相反，如果该区域是边缘区域，方差很大， σ k 2 \sigma_k^2 σk2​将远远大于 ϵ \epsilon ϵ。从而 a k ≈ 1 , b k ≈ 0 a_k≈1,b_k≈0 ak​≈1,bk​≈0。相当于在区域保持原有梯度。
    以上可以出: ϵ \epsilon ϵ为界定平滑区域和边缘区域的阈值。
• 图像去雾
  在图像去雾中，导向滤波一般用来细化透射率图像，以原图的灰度图为导向图，以粗投射率图为输入图，能得到非常精细的透射率图像。

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

导向滤波的实现

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

快速导向滤波的实现

由于导向滤波效率问题，何凯明博士在2015年，对其做了优化，基本原理是将导向图，输入图都进行下采样计算 ( a k , b k ) (a_k,b_k) (ak​,bk​)，然后对 ( a k , b k ) (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]