快速双边滤波——Python实现

介绍

课程的最后一个实验是处理雀斑，网上查找了很多方法，最后我选择了快速双边滤波。但是实验又不能直接调用 opencv 的库，因此，我参照了 这个博客，将用 C 写的快速双边滤波改用 Python 重新写了一遍。

快速双边滤波是啥呢，官方滴说，它是一种非线性的滤波方法。它最大的特点就是既使用了颜色的相似度，又利用了空间的距离相似度。也就是说，它在进行滤波的过程中，不光要考虑周围像素值与中点像素值的大小之差，还需要考虑空间上的距离，进而确定该点对中间点的影响因子。

至于为啥要用快速双边滤波呢，，因为不优化的话，实在是太慢了。。。有兴趣可以试一下。

代码地址：实验6

计算公式

刚刚已经讲过原理了，现在看一下它的计算过程。

它主要分为两大部分，一个是计算颜色值的相似度，公式如下：

$r(i,j,k,l) = exp(-\frac{||f(i,j) - f(k,l)||^2}{2\sigma_r^2})$r(i,j,k,l)=exp(−2σr2​∣∣f(i,j)−f(k,l)∣∣2​)

另一个是计算空间距离的相似度，也就是说，离得越近，相似度越高。公式如下：

$d(i,j,k,l) = exp(-\frac{(i-k)^2 + (j-l)^2}{2\sigma_d^2})$d(i,j,k,l)=exp(−2σd2​(i−k)2+(j−l)2​)

这里的 （i, j）代表要处理的像素点的坐标，（k，l）则是其周围一定范围内，可能影响到其值的像素点的坐标。千万不要被一堆字母搞晕了，逻辑还是挺清晰的。

总体公式如下：

$w(i,j,k,l) = d(i,j,k,l) * r(i,j,k,l) = exp(-\frac{(i-k)^2 + (j-l)^2}{2\sigma_d^2} - \frac{||f(i,j) - f(k,l)||^2}{2\sigma_r^2})$w(i,j,k,l)=d(i,j,k,l)∗r(i,j,k,l)=exp(−2σd2​(i−k)2+(j−l)2​−2σr2​∣∣f(i,j)−f(k,l)∣∣2​)

• $f(x,y)$f(x,y) 表示要处理的图像，f(x,y)表示图像在点(x,y)处的像素值
• $(k,l)$(k,l)为模板窗口的中心坐标
• $(i,j)$(i,j)为模板窗口的其他系数的坐标；
• $\sigma_r$σr​为高斯函数的标准差

这个就是矩阵中，像素点的系数。

优化代码

主要的优化部分有三个：

1. 空间距离的系数（$d(i,j,l,k)$d(i,j,l,k)）与像素点的值无关，只与矩阵的大小有关，因此可以首先将它计算出来，需要的时候再去查找。
2. 其次就是灰度值的计算，由于像素点过多，每次滤波，逐个计算的计算量特别大；而像素点的值只是在（0,255）区间内，因此，也可以使用查表的方式进行查找，这样只是（0,$255^2$2552）上面进行查找即可，不需要挨个计算。
3. 将二维的模板转换为一维，降低算法复杂度。（就是更改了一下索引的方式）

具体的 C 语音实现，还需要看文章开始的那篇博客，这里再粘贴一下：

图像处理基础(5)：双边滤波器

代码实现

完整的代码如下：

# -*- coding: UTF-8 -*-
import numpy as np
import cv2
import math

def bilateralFilter(img, radius, sigmaColor, sigmaSpace) : 
    B, G, R = cv2.split(img)
    B_tran, G_tran, R_tran = cv2.split(img)
    img_height = len(B)
    img_width = len(B[0])
    # 计算灰度值模板系数表
    color_coeff = -0.5 / (sigmaColor * sigmaColor)
    weight_color = []       # 存放颜色差值的平方
    for i in range(256) :
        weight_color.append(np.exp(i * i * color_coeff))
    # 计算空间模板
    space_coeff = -0.5 / (sigmaSpace * sigmaSpace)
    weight_space = []     # 存放模板系数
    weight_space_row = [] # 存放模板 x轴 位置
    weight_space_col = [] # 存放模板 y轴 位置
    maxk = 0
    for i in range(-radius, radius+1) :
        for j in range(-radius, radius+1) :
            r_square = i*i + j*j
            r = np.sqrt(r_square)
            weight_space.append(np.exp(r_square * space_coeff))
            weight_space_row.append(i)
            weight_space_col.append(j)
            maxk = maxk + 1
    # 进行滤波
    for row in range(img_height) :
        for col in range(img_width) :
            value = 0
            weight = 0
            for i in range(maxk) :
                m = row + weight_space_row[i]
                n = col + weight_space_col[i]
                if m < 0 or n < 0 or m >= img_height or n >= img_width :
                    val = 0
                else :
                    val = B[m][n]
                w = np.float32(weight_space[i]) * np.float32(weight_color[np.abs(val - B[row][col])])
                value = value + val * w
                weight = weight + w
            B_tran[row][col] = np.uint8(value / weight)
    # 绿色通道
    for row in range(img_height) :
        for col in range(img_width) :
            value = 0
            weight = 0
            for i in range(maxk) :
                m = row + weight_space_row[i]
                n = col + weight_space_col[i]
                if m < 0 or n < 0 or m >= img_height or n >= img_width :
                    val = 0
                else :
                    val = G[m][n]
                w = np.float32(weight_space[i]) * np.float32(weight_color[np.abs(val - G[row][col])])
                value = value + val * w
                weight = weight + w
            G_tran[row][col] = np.uint8(value / weight)
    # 红色通道
    for row in range(img_height) :
        for col in range(img_width) :
            value = 0
            weight = 0
            for i in range(maxk) :
                m = row + weight_space_row[i]
                n = col + weight_space_col[i]
                if m < 0 or n < 0 or m >= img_height or n >= img_width :
                    val = 0
                else :
                    val = R[m][n]
                w = np.float32(weight_space[i]) * np.float32(weight_color[np.abs(val - R[row][col])])
                value = value + val * w
                weight = weight + w
            R_tran[row][col] = np.uint8(value / weight)
    cv2.imshow("beauty_after", cv2.merge([B_tran, G_tran, R_tran]))
    cv2.imwrite("beauty_after.png", cv2.merge([B_tran, G_tran, R_tran]))

img = cv2.imread("beauty1.png")
cv2.imshow("original image", img)

# bilateralFilter(img, 5, 45, 100)
bilateralFilter(img, 3, 30, 80)

img = cv2.imread("beauty_after.png")
bilateralFilter(img, 3, 30, 80)

cv2.waitKey(0)

总结

三个维度的滤波不应该像我这样写，应该一起处理，我分成了三个通道进行处理，代码就显得有点长，但这样十分方便，如果有时间的话，还是不要像我这样写。

需要注意的一点是，这段代码在处理图片的效果上不是很好，相比于 opencv 的库，效果差了很多，尤其是处理大图片的时候，使用时需要注意。