导读
在图像滤波算法中,导向滤波、双边滤波、最小二乘滤波并称三大保边滤波器,他们是各向异性滤波器。相对于常见的均值滤波、高斯滤波等各向同性滤波器,他们最大的特点是在去除噪声的同时,能最大限度保持边缘不被平滑。本文讲解导向滤波及其应用。
总的来讲,导向滤波就是尽可能让输出图像的梯度和导向图相似,同时让输出图像的灰度(亮度)与输入图像相似,以此来保留边缘并且滤除噪声。
原理推导
我们先看下图:
输入图像
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=akIi+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∑(akIi+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∑((akIi+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δ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
由上面推到得出:
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}
akbk=∑i∈wk(Ii2+ϵ)∑i∈wkpiIi−bk∑i∈wkIi=∣w∣∑i∈wkpi−ak∑i∈wkIi=pk−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的均值;pk——窗口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∈wkpiIi−bk∑i∈wkIi
⇒
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+μkIi+ϵ)∑i∈wk(piIi−pkIi)
⇒
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−μkIi−μkIi+μk2+ϵ)∑i∈wk(piIi−pkIi+μkpi−μkpk)
上下两边同除以
∣
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∈wkpiIi−pkμk+μkpk−μkpk
最后:
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∈wkpiIi−μkpk(5)
b
k
=
p
‾
k
+
a
k
μ
k
(6)
\color{red}{b_k = \overline{p}_k+a_k\mu_k} \tag{6}
bk=pk+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=wk1k,i∈wk∑(akIi+bk)=aiIi+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+ϵσk2bk=(1−ak)pk
考虑以下两种情况:- 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≈pk。相当于对该区域作均值滤波。
- 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]