OpenCV—Python 导向滤波

一、导向滤波原理

导向滤波是使用导向图像作为滤波内容图像，在导向图像上实现局部线性函数表达，实现各种不同的线性变换，输出变形之后的导向滤波图像。根据需要，导向图像可以跟输入图像不同或者一致。

公式及推导

假设 III 是导向图像、ppp 是输入图像、qqq 是导向滤波输出图像；导向滤波是作为局部线性模型描述 导向图像 III 与 输出图像 qqq 之间的关系。

对于任意像素 I\rm II 来说，ωk\omega_kωk窗口下的线性变换可以表述如下：
qi=akIi+bkq_i = a_k I_i + b_kqi=akIi+bk 其中 (ak,bk)(a_k,b_k)(ak,bk)是窗口 ωk\omega_kωk 范围内的参数常量。
为了寻找线性相关性，窗口 ωk\omega_kωk 定义的损失函数为：
E(ak,bk)=∑i∈ωk((akIi+bk−pi)2+ϵak2)E(a_k,b_k) = \sum_{i \in \omega_k} ((a_k I_i + b_k-p_i)^2+\epsilon a_k^2)E(ak,bk)=i∈ωk∑((akIi+bk−pi)2+ϵak2)
其中：
ωk\omega_kωk是对 aka_kak 值过大时侯的正则化补偿。
εak2\varepsilon a_k^2εak2 是抑制 aka_kak 值过大的。
ε\varepsilonε 是调整图的模糊程度与边缘检测精度的参数。
如果导向图 I\rm II 没有边缘信息，输出均值模糊结果；
如果导向图 I\rm II 包含边缘信息，边缘信息则迁移到输出图像中实现边缘保留滤波；
上述的损失函数可以被看成一个线性回归问题，其中两个参数的求解如下：

ak=1ωk∑i∈ωkIiPi−μkpˉkσk2+ϵa_k = \frac{\frac{1}{\omega_k}\sum_{i \in \omega_k} I_iP_i-\mu_k \bar{p}_k}{\sigma_k^2 + \epsilon}ak=σk2+ϵωk1∑i∈ωkIiPi−μkpˉk
bk=pˉk−akμkb_k = \bar{p}_k - a_k\mu_kbk=pˉk−akμk
μk\mu_kμk 与 σk2\sigma_k^2σk2 是导向图在 ωk\omega_kωk 窗口大小均值与方差
∣ω∣|\omega|∣ω∣表示窗口内像素总数
pˉk=1∣ω∣∑i∈ωkpi\bar{p}_k = \frac{1}{|\omega|}\sum_{i \in \omega_k}p_ipˉk=∣ω∣1∑i∈ωkpi 在 ωk\omega_kωk 窗口内输入图像像素均值

使用线性相关参数 (ak,bk)(a_k,b_k)(ak,bk)，滤波输出图像就可以通过qi=akIi+bkq_i = a_k I_i + b_kqi=akIi+bk 线性模型得到。
针对不同的窗口大小我们就会得到不同的 qiq_iqi 值，所以通过它的均值作为最终的输出结果：
qi=1ωk∑i∈ωk(akIi+bk)=aˉiIi+bˉiq_i = \frac{1}{\omega_k}\sum_{i \in \omega_k}(a_k I_i + b_k) = \bar{a}_iI_i+\bar{b}_iqi=ωk1i∈ωk∑(akIi+bk)=aˉiIi+bˉi

最终导向滤波公式为：
aˉ，bˉi\bar{a}_，\bar{b}_iaˉ，bˉi 是所有像素点 iii 上重叠窗口相关因子的均值。

导向滤波算法实现的一般步骤为：

读取导向图像 I\rm II 与输入图像 P\rm PP
积分图计算 I\rm II 的均值与方差、输入图像 P\rm PP的均值、 I\rm II 与 P\rm PP的乘积 IP\rm IPIP
计算线性相关因子 a\rm aa 与 b\rm bb
a=(IP−ImeanPmean)/(IVar+ϵ)a=(IP-I_{mean}P_{mean})/(I_{Var} +\epsilon )a=(IP−ImeanPmean)/(IVar+ϵ)
b=Pmean−aImeanb=P_{mean}-\textrm{a}I_{mean}b=Pmean−aImean
计算a与b的均值
使用均值得到导向滤波结果 Q=amean∗I+bmeanQ = a_{mean}*I+b_{mean}Q=amean∗I+bmean

导向滤波最常用四个功能是：

边缘保留滤波
图像去噪声
图像边缘羽化
图像增强（对比度）

import numpy as np
import cv2def guideFilter(I, p, winSize, eps):mean_I = cv2.blur(I, winSize)      # I的均值平滑mean_p = cv2.blur(p, winSize)      # p的均值平滑mean_II = cv2.blur(I * I, winSize) # I*I的均值平滑mean_Ip = cv2.blur(I * p, winSize) # I*p的均值平滑var_I = mean_II - mean_I * mean_I  # 方差cov_Ip = mean_Ip - mean_I * mean_p # 协方差a = cov_Ip / (var_I + eps)         # 相关因子ab = mean_p - a * mean_I            # 相关因子bmean_a = cv2.blur(a, winSize)      # 对a进行均值平滑mean_b = cv2.blur(b, winSize)      # 对b进行均值平滑q = mean_a * I + mean_breturn qif __name__ == '__main__':eps = 0.01winSize = (5,5)image = cv2.imread(r'./5921.png', cv2.IMREAD_ANYCOLOR)image = cv2.resize(image, None,fx=0.7, fy=0.7, interpolation=cv2.INTER_CUBIC)I = image/255.0        #将图像归一化p =IguideFilter_img = guideFilter(I, p, winSize, eps)# 保存导向滤波结果guideFilter_img  = guideFilter_img  * 255guideFilter_img [guideFilter_img  > 255] = 255guideFilter_img  = np.round(guideFilter_img )guideFilter_img  = guideFilter_img.astype(np.uint8)cv2.imshow("image",image)cv2.imshow("winSize_5", guideFilter_img )cv2.waitKey(0)cv2.destroyAllWindows()

当 I=pI=pI=p 时，导向滤波就变成了边缘保持的滤波操作，此时原来求出的aaa 和 bbb 的表达式就变成了：
ak=σk2σk2+ϵa_k=\frac{σ_k^2}{σ_k^2+ϵ}ak=σk2+ϵσk2

bk=(1−ak)μkb_k=(1−a_k)μ_kbk=(1−ak)μk
考虑两种情况：

情况1：高方差区域，即表示图像 III 在窗口 wkw_kwk 中变化比较大，此时我们有σk2>>ϵσ_k^2>>ϵσk2>>ϵ，于是有ak≈1a_k≈1ak≈1和bk≈0b_k≈0bk≈0。
情况2：平滑区域（方差不大），即图像 III 在窗口 wkw_kwk 中基本保持固定，此时有σk2<<ϵσ_k^2<<ϵσk2<<ϵ，于是有ak≈0a_k≈0ak≈0和bk≈μkb_k≈μ_kbk≈μk 。

也就是说在方差比较大的区域，保持值不变，在平滑区域，使用临近像素平均（也就退化为普通均值滤波）。

快速导向滤波

通过下采样减少像素点，计算 meana\rm mean_ameana & meanb\rm mean_bmeanb 后进行上采样恢复到原有的尺寸大小。
假设缩放比例为s,那么缩小后像素点的个数为N/s2，N/s^2，N/s2，那么时间复杂度变为O(N/s2)O(N/s^2)O(N/s2)

import cv2
import numpy as npdef guideFilter(I, p, winSize, eps, s):# 输入图像的高、宽h, w = I.shape[:2]# 缩小图像size = (int(round(w * s)), int(round(h * s)))small_I = cv2.resize(I, size, interpolation=cv2.INTER_CUBIC)small_p = cv2.resize(I, size, interpolation=cv2.INTER_CUBIC)# 缩小滑动窗口X = winSize[0]small_winSize = (int(round(X * s)), int(round(X * s)))# I的均值平滑 p的均值平滑mean_small_I = cv2.blur(small_I, small_winSize)mean_small_p = cv2.blur(small_p, small_winSize)# I*I和I*p的均值平滑mean_small_II = cv2.blur(small_I * small_I, small_winSize)mean_small_Ip = cv2.blur(small_I * small_p, small_winSize)# 方差、协方差var_small_I = mean_small_II - mean_small_I * mean_small_Icov_small_Ip = mean_small_Ip - mean_small_I * mean_small_psmall_a = cov_small_Ip / (var_small_I + eps)small_b = mean_small_p - small_a * mean_small_I# 对a、b进行均值平滑mean_small_a = cv2.blur(small_a, small_winSize)mean_small_b = cv2.blur(small_b, small_winSize)# 放大size1 = (w, h)mean_a = cv2.resize(mean_small_a, size1, interpolation=cv2.INTER_LINEAR)mean_b = cv2.resize(mean_small_b, size1, interpolation=cv2.INTER_LINEAR)q = mean_a * I + mean_breturn qif __name__ == '__main__':eps = 0.01winSize = (16,16)       #类似卷积核（数字越大，磨皮效果越好）image = cv2.imread(r'./5921.png', cv2.IMREAD_ANYCOLOR)image = cv2.resize(image,None,fx=0.8,fy=0.8,interpolation=cv2.INTER_CUBIC)I = image/255.0       #将图像归一化p =Is = 3 #步长guideFilter_img = guideFilter(I, p, winSize, eps,s)# 保存导向滤波结果guideFilter_img = guideFilter_img  * 255         #(0,1)->(0,255)guideFilter_img[guideFilter_img  > 255] = 255    #防止像素溢出guideFilter_img = np.round(guideFilter_img )guideFilter_img = guideFilter_img.astype(np.uint8)cv2.imshow("image",image)cv2.imshow("winSize_16", guideFilter_img )cv2.waitKey(0)cv2.destroyAllWindows()输出效果图如上：winSize_16

为了写这篇博客，查了许久的资料，可气的是，那些都是只写了一半的代码，如今代码运行通了，与诸位共勉。

鸣谢
https://blog.csdn.net/qq_40755643/article/details/83831071

特别鸣谢
http://blog.sina.com.cn/s/blog_734f70550102wof3.html
https://blog.csdn.net/baimafujinji/article/details/74750283