相关原理
具体详细的解释可以看论文A Gentle Introduction to Bilateral Filtering
and its Applications,简单的解释可以参考下面的链接:
链接: Bilateral Filters(双边滤波算法)的超简单原理,学不会你打我
对我来说,这个算法本质上也是一个矩阵在图像上移动并进行卷积的过程,这个算法的关键就在于对其计算公式的理解,其公式如下所示:
[图片]
我们来详细解释一下上面两个代码,p代表正在进行滤波的像素点,即矩阵块移动的中心位置,q代表矩阵中的各个像素点的位置。
B
F
[
I
p
]
BF[I_p]
BF[Ip]代表经过滤波后的新的像素值,
W
p
W_p
Wp代表矩阵块移动到位置p时的权重的和,用来进行归一化。
剩下的我们可以拆开来看,
G
σ
s
(
∣
∣
p
−
q
∣
∣
)
G_{\sigma_{s}}(||p-q||)
Gσs(∣∣p−q∣∣)其实等同于进行高斯滤波的卷积核计算,
G
σ
r
(
∣
∣
I
p
−
I
q
∣
∣
)
G_{\sigma_{r}}(||I_p-I_q||)
Gσr(∣∣Ip−Iq∣∣)则是通过高斯函数计算两个像素点之间像素值的相关性,最后将两者作为参数乘上q点的像素值,遍历整个参考框(卷积核范围内的像素值),最后进行加权平均,得到这一个点的新像素值。最后遍历整幅图。
几点简化计算:
1、对于高斯函数函数的计算如下,
[图片]
我们在实际编写滤波算法时,高斯前面的系数实际上时不用计算的,因为当参数固定时,前面的参数也是固定。
公式3可以看作分子和分母同时除以两个高斯函数的系数。
2、计算距离的高斯时,需要计算两个像素点的距离,
(
(
x
0
−
x
1
)
2
+
(
y
0
−
y
1
)
2
)
\sqrt((x_0-x_1)^2+(y_0-y_1)^2)
((x0−x1)2+(y0−y1)2)这里的开根号可以与高斯函数中的平方相抵消,节省计算。
相关代码
import cv2
import os
import math
import numpy as np
def add_noisy_one_img(src, mean, sigma):
# 获取图片的高度和宽度
height, width, channels = src.shape
gauss = np.random.normal(mean, sigma, (height, width, channels))
noisy_img = src + gauss
# noisy_img = noisy_img.astype(np.float32)
# 将加噪图片的像素值缩放放到cv2接收的uint8范围
noisy_img1 = np.clip(noisy_img, a_min=0, a_max=255)
noisy_img1 = noisy_img1.round().astype(np.uint8)
return noisy_img1
# 对图像进行边缘镜像处理
def set_border_sym(image_data, ksize):
h_, w_, c_ = image_data.shape
border = int((ksize - 1) / 2)
result_data = np.zeros((h_ + ksize - 1, w_ + ksize - 1, c_))
for c in range(c_):
for i in range(h_ + 2 * border):
if i < border:
y_ori = border - i
elif i > h_ + border - 1:
y_ori = (h_ + border - 1) - (i - (border + h_ - 1)) - border
else:
y_ori = i - border
for j in range(w_ + 2 * border):
if j < border:
x_ori = border - j
elif j > w_ + border - 1:
x_ori = (w_ + border - 1) - (j - (border + w_ - 1)) - border
else:
x_ori = j - border
result_data[i, j, c] = image_data[y_ori, x_ori, c]
return result_data.astype("uint8")
def bilateral_filtering(image_pad_data, k_size, sigma_c, sigma_s):
# image_pad_data = image_pad_data.astype("float32")
h_, w_, c_ = image_pad_data.shape
result_data = np.zeros_like(image_pad_data)
border = int((k_size - 1) / 2)
# 对图像分图层进行处理
for c in range(c_):
# 输入的是填充后的图像,但索引要从原图像素处开始
for y in range(border, h_ - border):
for x in range(border, w_ - border):
# 获得(x, y)处的像素参考框,然后与k逐像素相乘并求和
sum_w = 0
sum_pixel = 0
for i in range(-border, border + 1):
for j in range(-border, border + 1):
ref_y = y + i
ref_x = x + j
# diff = np.sqrt(i**2 + j**2)
dist2_s = math.exp(-((i**2 + j**2) / (2 * sigma_s**2)))
p_ref = image_pad_data[ref_y, ref_x, c]
p_ori = image_pad_data[y, x, c]
dist2_c = math.exp(-(p_ref-p_ori)**2 / (2 * sigma_c**2))
sum_w += dist2_s * dist2_c
sum_pixel += dist2_s * dist2_c * p_ref
a = sum_pixel / sum_w
result_data[y, x, c] = a
return result_data[border:h_ - border, border:w_ - border, :].astype("uint8")
if __name__ == "__main__":
ori_image = cv2.imread("../../data/image/lena.png")
noisy_image = add_noisy_one_img(ori_image, 0, 15)
# cv2自带的均值滤波 cv2.bilateralFilter(原始图像,核大小)
cv2_result = cv2.bilateralFilter(noisy_image, 5, 30, 150)
sym_border_image = set_border_sym(noisy_image, 5)
result_data_bil = bilateral_filtering(sym_border_image, 5, 30, 150)
cv2.imshow("sym_border_image", sym_border_image)
cv2.imshow("noisy_image", noisy_image)
cv2.imshow("cv2_result", cv2_result)
cv2.imshow("result_data_bil", result_data_bil)
cv2.imwrite("my_bil.png", result_data_bil)
cv2.waitKey(0)
实际结果
上述代码具有边缘填充操作,在使用双边滤波前需进行边缘填充处理,填充的边缘数为滤波框d的
(
(
d
−
1
)
/
2
)
((d-1)/2)
((d−1)/2)。
同时我们注意到有一个被注释的语句:
image_pad_data = image_pad_data.astype(“float32”)
可以试试取消注释结果会如何。
实际结果如下:
实际代码可以参考github