详细的引导滤波
在做毕设的过程中，由于需要相对准确地提取一张图像的光照分量（光照分量一般是低频部分，所以平滑滤波可以提取，高斯滤波边缘保持能力很差），便找了许多边缘保持滤波算法，后来发现了何凯明大神提出的引导滤波算法，时间复杂度仅为O(N)，N为图像像素总数，也即与滤波半径无关，极大地提高了运算速度，就学习了一下，顺便记录下自己的思考。主要参考文献是何凯明的《Guided Image Filtering》以及https://www.cnblogs.com/riddick/p/8367591.html。
引导滤波算法利用引导图像与滤波输出图像之间的局部线性关系（这样就能保证输出图像保留引导图像的局部特征了），将计算出来的输出图像与待滤波图像（即原图像）做最小二乘，使输出图像尽可能逼近原图像。如果引导图像选用原图像，那么即可实现在保留原图像局部特征的基础上实现对原图像的滤波平滑。我们来看原论文的公式：
其中ωk是指某一个滤波窗口，在这个滤波窗口内的所有像素点i，输出图像q均是引导图像I的线性变换，这个变化由ak和bk决定。我们对上式两边求梯度，得到：
可以看出，ak决定了梯度保留能力。如果ak较大，则梯度保留效果好；反之，就是平滑效果好了。看到这儿，我们应该知道ak和bk就是自适应调节因子了。我们希望在边缘处ak大些，其他地方ak小些。那ak，bk究竟应该怎样取值，才能实现这样一个自适应调节作用呢？接着往下看：
因为我们是要对原图像进行滤波，所以在满足上面局部线性的基础上，还需要使输出图像尽可能逼近原图像，作者采用的是最小二乘方法，使输出图像q与原图像p之间的均方误差最小。我们来看原论文的公式：
可以看到，除了均方误差之外，作者还添加了含ϵ的一项。作者在原文中解释说这是一个正则化参数，用来防止ak过大。我想大部分人应该对这个很熟悉，这不就是深度学习里损失函数的正则项嘛，只不过那个是用来惩罚网络中的权重，防止过拟合罢了。另外我们试想一下，如果没有ϵ这一项，引导图像选用原图像，那么最小化上式，ak一定为1，bk一定为0，因为没有谁比自身更接近自身了嘛，这样带来的结果就是没有任何平滑效果，输出图像就是原图像本身。所以ϵ是一个非常重要的参数，选取不当会影响输出图像，后面会讲如何选取。
好了，现在最小化上式，即可求出ak，bk了。因为这是一个线性模型的最小二乘求解，它是存在闭式解的，不用搜索求解。我自己当时也是手推过的，就放个截图上来吧。（字比较丑）
总结一下，也就是如下形式：
其中，|ω|为窗口内像素个数，uk，pk分别为窗口内引导图像与原图像的均值，是窗口内引导图的方差。
求出ak，bk后，根据前面局部线性模型就可求出输出图像啦。不过由于同一个像素点会被许多不同滤波窗口所包含，每一个窗口里的ak，bk又不相同，所以作者又对ak，bk取了均值再计算的。如下图所示：
可以看到，ak的分子其实就是I和p的协方差，如果引导图像选择原图像，分子就变成了
这其实就是方差。（
）
这么一来，ak就变成了如下形式：
现在就很明显了，当原图像信息丰富的区域（包括边缘），方差相应会大，根据上式，ak自然也就大了（糖水加糖），原图像的梯度信息就得以保留。原图像平坦的区域，ak相应会小，bk相应增大，当ak减小为0时，bk就为μk了，这样相当于均值滤波了，平滑效果最强。
引导滤波算法快的关键之处在于计算ak，bk所需的均值与方差都可以通过盒式滤波去实现，盒式滤波是一种与滤波半径大小无关的算法，其实现方法类似积分图思想，都是一种以空间换取时间的方式。因为相邻像素之间有许多重复计算的地方，这些地方可以暂时记录下来，后面需要用到时直接读取就行了。可参考这篇文章https://www.cnblogs.com/lwl2015/p/4460711.html
当引导图为三通道图像时，需要对上述模型进行简单修改，其实就是矢量运算而已，因为同一个像素点处会有三个灰度值分别对应三个通道。
这里，
是一个3*1的向量，对应三个通道灰度值与均值。
是单位矩阵，
为3*3的三通道协方差矩阵，其他类似。这样
就是一个3*1的向量，转置后与3*1的
相乘，求出来的bk依然是常数。
原文作者为了体现引导滤波是一种线性滤波器，又将整个模型改写为如下形式（所谓大道至简）：
其中，Wij为下式：
Wij的详细推导，感兴趣的可以去论文里看看。
回到之前一个问题，ϵ应该怎么选取？因为它的取值决定了ak的自适应作用，ϵ很大的时候，不管的大小，ak都接近于0，所以十分重要。个人建议一开始可以选取一个跟最大
同量级的值，这样效果明显。
后面也有人对引导滤波进行改进，对ϵ也进行了自适应调节。在图像信息丰富的区域，使ϵ小一些，让ak取较大值，反之让ϵ大一些。感兴趣的也可以看看。
下图是我自己的一个实现案例，我是利用HSV颜色系统提取出图片的亮度分量再处理的。
最后是我的代码部分：
import cv2 as cv
from matplotlib import pyplot as plt
import numpy as np
R = 16 #滤波半径
r =2 * R +1
s = 1 #快速引导滤波（如果大于0）
def read(path):
"""读入图片"""
src = cv.imread(path)
return src
def To_HSV(image):
"""转化为HSV空间"""
dst = cv.cvtColor(image,cv.COLOR_BGR2HSV) #opencv是BGR模式
H,S,V = cv.split(dst)
cv.imshow('V',V)
return H,S,V
def To_BGR(H,S,V):
dst = cv.merge([H, S, V])
dst = cv.cvtColor(dst, cv.COLOR_HSV2BGR)
cv.imshow('dst', dst)
return dst
def LSE(V):
V_down = cv.resize(V, (V.shape[1] // s, V.shape[0] // s))
V_32 = V_down.astype(np.float32)
V_2 = np.square(V_32)
dst1 = cv.boxFilter(V_2, ddepth=-1, ksize=(r, r), normalize=True,borderType=cv.BORDER_DEFAULT) # 盒式滤波（True为归一化均值滤波），效率高之关键,类似积分图
#print(dst1)
u_k = cv.boxFilter(V_32, ddepth=-1, ksize=(r, r), normalize=True,borderType=cv.BORDER_DEFAULT) #默认为cv.BORDER_REFLECT101，镜像
u_k_2 = np.square(u_k)
sigma_k = dst1 - u_k_2 # D（X）= E（X2）- E（X）2
print('最大方差为：',np.max(sigma_k)) #kesi的取值应与此同量级,以使方差发挥自适应调节作用（ak公式）
kesi = 10000 #重要参数（太小则全部保留梯度信息，没有任何平滑(ak==1)；太大则等价于均值滤波了）
a_k = (dst1 - u_k*u_k)/(sigma_k + kesi)
b_k = u_k - a_k * u_k
#print(a_k)
return a_k,b_k,V_down
def light(a_k,b_k,V_down,V):
a = cv.boxFilter(a_k, ddepth=-1, ksize=(r, r), normalize=True)
b = cv.boxFilter(b_k, ddepth=-1, ksize=(r, r), normalize=True)
I_down = a * V_down + b
#cv.imshow('I_down',I_down)
I = cv.resize(I_down, (V.shape[1], V.shape[0]))
I = np.where(I > 255, 255, I) #可用clip函数
I = np.where(I < 0, 0, I)
I_ = I.astype(np.uint8)
cv.imshow('I', I_)
return I
if __name__ == '__main__':
path = 'samples/8.png' #修改为自己路径
src = read(path)
cv.imshow('src',src)
H,S,V = To_HSV(src)
ak,bk,V_down = LSE(V)
I = light(ak,bk,V_down,V)
cv.waitKey(0)
cv.destroyAllWindows()