双边滤波的python实现

相关原理
具体详细的解释可以看论文A Gentle Introduction to Bilateral Filtering
and its Applications，简单的解释可以参考下面的链接：
链接: Bilateral Filters（双边滤波算法）的超简单原理，学不会你打我

对我来说，这个算法本质上也是一个矩阵在图像上移动并进行卷积的过程，这个算法的关键就在于对其计算公式的理解，其公式如下所示：
[图片]
我们来详细解释一下上面两个代码，p代表正在进行滤波的像素点，即矩阵块移动的中心位置，q代表矩阵中的各个像素点的位置。
$$BF[I_p]$$代表经过滤波后的新的像素值，$W_p$代表矩阵块移动到位置p时的权重的和，用来进行归一化。
剩下的我们可以拆开来看，$$G_{{\sigma }_s}(||p-q||)$$其实等同于进行高斯滤波的卷积核计算，$$G_{{\sigma }_r}(||I_p-I_q||)$$则是通过高斯函数计算两个像素点之间像素值的相关性，最后将两者作为参数乘上q点的像素值，遍历整个参考框（卷积核范围内的像素值），最后进行加权平均，得到这一个点的新像素值。最后遍历整幅图。
几点简化计算：
1、对于高斯函数函数的计算如下，
[图片]
我们在实际编写滤波算法时，高斯前面的系数实际上时不用计算的，因为当参数固定时，前面的参数也是固定。
公式3可以看作分子和分母同时除以两个高斯函数的系数。
2、计算距离的高斯时，需要计算两个像素点的距离，$$\sqrt{(}(x_0-x_1{)}^2+(y_0-y_1{)}^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)$。
同时我们注意到有一个被注释的语句：

image_pad_data = image_pad_data.astype(“float32”)

可以试试取消注释结果会如何。
实际结果如下：
实际代码可以参考github