1. 简介
之前已经介绍过了双边滤波核联合双边滤波，其中双边滤波是一种非线性的保边滤波器，而联合双边滤波相当于将值域高斯核的来源从原始影像替换成另外一副引导图。本文介绍的导向滤波，其与联合双边滤波类似，也需要除原始影像之外另外一副引导图，是一种保边滤波器，当然其也可以用作图像去雾、HDR压缩等。
2. 算法原理
2.1 导向滤波框架
在算法框架中，要对p进行滤波而得到q，还得需要一个引导图像I。此时，滤波器的数学公式为：
注意，这里的Wij(I)就表示由引导图像I来确定加权平均运算中所采用的权值，并且权值与p无关，该滤波器相对于p是线性的。注意，引导图像可以是单独的一幅图像，也可以输入的图像p本身。当引导图像就是p本身的时候，导向滤波就变成了一个Edge-perserving的滤波器。
2.2 重要假设
导向滤波器在引导图像I和滤波输出q之间在一个二维窗口内是一个局部线性模型，假设q是以像素k为中心的窗口wk中的I的线性变换：
ak和bk是该线性函数的系数，窗口半径为r。这样假设是因为我们希望引导图像提供信息来指示哪里是图像边缘，哪里是平滑区域，接下来我们相应的保留边缘，磨平平滑区域，只有当I和q之间是线性关系的这种引导关系才有意义。
2.3 求解单个窗口线性变换系数
现在已知的是I和p，要求的是q。而如果能求得参数ak和bk，显然就能通过I和q之间的线性关系来求出q。另一方面，p是q受到噪声污染而产生的退化图像，假设噪声是n，则有qi=pi−ni。根据无约束图像复原的方法，我们需要最小化q核p之间的差异，也就是说，需要最小化以下函数：
进而我们通过最小二乘回归求解线性系数ak核bk，也就是分别对ak和bk求偏导，并令其导数为零：
将上式展开成如下形式：
通过上图第二个等式先求解bk如下：
进而通过克莱姆法则求解ak：
上式中，为了表示方便，我们使用n表达窗口wk中的像素总数。此外，uk是I中窗口wk的均值，σ2为I中窗口wk的方差，pk¯是待滤波图像p在窗口wk中的均值。
之前的求解过程我们一直忽略了正则化系数epsilon，加上正则化项之后的ak和bk结果如下：
2.4 求解整图各位置线性变换系数
此外，在计算每个窗口的线性系数时，以滑动窗口为例，我们可以发现一个像素会被多个窗口包含，也就是说，每个像素都由多个线性函数所描述。因此，如之前所说，要具体求某一点的输出值时，只需将所有包含该点的线性函数值平均即可，如下 ：
2.5 算法流程
3. 算法分析
3.1 边缘保留特性
这可以直观地解释如下，考虑I = p的情况。很明显，如果ε= 0，那么ak = 1和bk = 0。如果epsilon> 0，我们可以考虑两种情况：
平坦图像块：如果图像I在wk中是常数，ak=0, bk=pk¯,相当于平滑了原始图像，起到平坦区域模糊平滑的作用；
高方差区域(边缘)：如果图像I在wk中发生了很大变化，并且近似P也在wk中有边缘，则ak接近于1，而bk接近于0，此时q=I，则相应的边缘就会保留；
3.2 O(N)时间复杂度
O(N)时间意味着时间复杂度与窗口半径r无关，因此我们可以在应用程序中自由使用任意内核大小。
正常来看，当内核变大时，其计算复杂度会增加。但是，我们不是直接执行卷积，而是根据其定义Σ(5)(6)(8)计算滤波器输出。这些方程中的所有求和都是盒式滤波器。我们应用O(N)时间积分图像技术来计算盒式滤波器的输出。因此，引导滤波器可以在O(N)时间内计算。
这里简介O(N)时间积分图像技术：
积分图像(integral image)是一种快速计算矩形区域之和的数据结构，常利用它对算法进行加速，常见于均值滤波。
以下主要以灰度图像来讲解。一个图像内矩形区域的积分是指这个矩形区域内所有灰度值的和，下图左图为原始图像，有图为积分图像(从左到右，从上到下)：
如果每个像素点都计算其积分，其运算量就太大了，为了减少计算量，经过推导给出了如下计算方法：如图所示：
区域1 : = sum(A);
区域2 : = sum(A + B);
区域3 : = sum(A + C);
区域4 : = sum(A + B + C + D);
所以，如果需要计算D区域中的灰度和，则
sum(D) = 区域4 - 区域2 - 区域3 + 区域1 (都是灰度值)
现在可以通过将总和除以该窗口中的总像素来轻松计算平均值。
4. python实现
import numpy as np
import matplotlib.pyplot as plt
import cv2
class GuidedFilter:
"""
References:
K.He, J.Sun, and X.Tang. Guided Image Filtering. TPAMI'12.
"""
def __init__(self, I, radius, eps):
"""
Parameters
----------
I: NDArray
Guided image or guided feature map
radius: int
Radius of filter
eps: float
Value controlling sharpness
"""
if len(I.shape) == 2:
self._Filter = GrayGuidedFilter(I, radius, eps)
def filter(self, p):
"""
Parameters
----------
p: NDArray
Filtering input which is 2D or 3D with format
HW or HWC
Returns
-------
ret: NDArray
Filtering output whose shape is same with input
"""
p = (1.0 / 255.0) * np.float32(p)
if len(p.shape) == 2:
return self._Filter.filter(p)
class GrayGuidedFilter:
"""
Specific guided filter for gray guided image.
"""
def __init__(self, I, radius, eps):
"""
Parameters
----------
I: NDArray
2D guided image
radius: int
Radius of filter
eps: float
Value controlling sharpness
"""
self.I = (1.0 / 255.0) * np.float32(I)
self.radius = radius * 2 + 1
self.eps = eps
def filter(self, p):
"""
Parameters
----------
p: NDArray
Filtering input of 2D
Returns
-------
q: NDArray
Filtering output of 2D
"""
# step 1
meanI = cv2.blur(self.I, (self.radius, self.radius))
meanp = cv2.blur(p, (self.radius, self.radius))
corrI = cv2.blur(self.I * self.I, (self.radius, self.radius))
corrIp = cv2.blur(self.I * p, (self.radius, self.radius))
# step 2
varI = corrI - meanI * meanI
covIp = corrIp - meanI * meanp
# step 3
a = covIp / (varI + self.eps)
b = meanp - a * meanI
# step 4
meana = cv2.blur(a, (self.radius, self.radius))
meanb = cv2.blur(b, (self.radius, self.radius))
# step 5
q = meana * self.I + meanb
return q
def double2uint8(I, ratio=1.0):
return np.clip(np.round(I * ratio), 0, 255).astype(np.uint8)
if __name__ == "__main__":
image = cv2.imread('lena512.bmp', 0)
plt.figure()
for i, radius in enumerate([2, 4, 8]):
for j, e in enumerate([0.1**2, 0.2**2, 0.4**2]):
GF = GuidedFilter(image, radius, e)
plt.subplot(3, 3, i*3+j+1)
plt.axis('off')
plt.title('radius: %d, epsilon: %.2f' % (radius, e))
plt.imshow(GF.filter(image), cmap='gray')
plt.show()
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
