【图像去噪专题】图像降噪：空间域滤波算法

1. 背景介绍

在图像处理中，空间域滤波是通过直接分析图像像素与其邻域像素之间的关系来实现图像降噪的。这种方法通过在图像的二维空间上应用一系列滤波器（filter）或卷积核（kernel）来平滑图像，从而降低噪声影响。

局部卷积通常会设计卷积核，在图像的每个像素的局部邻域实施像素级的数值运算操作。通常这类算法在过滤噪声的同时图像细节也会有一定衰减，体现在图像清晰度降低和细节丢失。

保边滤波是图像处理领域中一类旨在去除图像噪声的同时尽可能保持图像细节和边缘信息的滤波技术，这种滤波方法主要用于在需要消除图像噪声但又不希望模糊重要边缘的情况下，比如医学影像分析、遥感图像处理、计算机视觉以及摄影后期等领域。

在实际应用中，选择何种空间域滤波方法取决于图像噪声的类型和图像本身的特性，以及对降噪效果和保持细节的要求。

---

2. 局部卷积

高斯滤波和均值滤波是去除噪声、平滑数据的重要工具。作为经典的线性滤波器，其输出是输入信号的线性组合。对于图像处理而言，这意味着滤波器对图像上的每一个像素点的值计算时，该点的新值是它周围邻域像素按照某种权重系数加权求和的结果。

2.1 高斯滤波

高斯滤波器就是一种广泛应用的线性滤波器，它使用高斯函数作为权重分配标准，对图像进行平滑处理。高斯滤波器采用高斯分布作为卷积核，根据临近像素到中心像素的距离赋予不同的权重，距离越近的像素影响越大，距离越远的像素影响越小。这一性质可以减轻图像失真，如果离中心像素很远的像素点有过大权重，过度的均匀会导致边缘模糊。

def gaussian_kernel(kernelsize, sigma):
    kernel = np.zeros(shape=(kernelsize, kernelsize), dtype=np.float)
    radius = kernelsize // 2
    for y in range(-radius, radius + 1):  # [-r, r]
        for x in range(-radius, radius + 1):
            val = 1.0 / (2 * np.pi * sigma ** 2) * np.exp(-1.0 / (2 * sigma ** 2) * (x ** 2 + y ** 2))
            kernel[y + radius, x + radius] = val
    kernel = kernel / np.sum(kernel)
    return kernel

def gaussian_blur(img, sigma = 5, kernelsize = 5):
    row, col = img.shape
    nr_img = img.copy()
    bound = (kernelsize - 1) // 2
    kernel = gaussian_kernel(kernelsize, sigma)
    print(kernel)
    for i in range(bound, row - bound):
        for j in range(bound, col - bound):
            nr_img[i, j] = np.sum(img[i-bound:i+bound+1, j-bound:j+bound+1] * kernel)

    return nr_img

kernel size为5时，可以看到高斯滤波边缘损失较严重

2.2 均值滤波

均值滤波是最简单直接的空间域滤波方法，通过计算图像中每个像素点邻域内所有像素的平均值来代替该像素点的值，以此达到平滑图像、去除噪声的目的。但由于它对图像边缘和细节处理过于粗略，可能会导致图像变得模糊。

2.3 中值滤波

中值滤波是一种非线性滤波方法，它将图像中每个像素点邻域内的像素值按大小排序后取中间值替代原像素值。这种方法对椒盐噪声（salt-and-pepper noise）有很好的抑制效果，同时它可以保留边缘信息，减轻边缘模糊的现象。

---

3. 保边滤波

保边滤波技术的关键在于设计出能够在减少噪声的同时保持图像细节特别是边缘信息的滤波算法，从而实现图像去噪的同时提升图像整体的视觉质量和后续处理的准确性。

3.1 双边滤波

双边滤波 (Bilateral Filter) 是保边滤波器的典型代表之一，它结合了图像的空间邻近性和像素间的灰度相似性两个准则来进行滤波。不同于传统的线性滤波器（如高斯滤波器只考虑空间邻近性），双边滤波器对每一个像素不仅考虑其周围像素的距离（空间域），还考虑它们之间的灰度差值（范围域）。通过这种方式，它能够在抑制噪声的同时避免模糊边缘，因为在图像边缘处，相邻像素间灰度变化较大，因此双边滤波器给予这些像素较低的滤波强度。

与高斯滤波相似，p和q代表距离，而Ip和Iq代表像素差异

def bilateral_blur(img, sigma_s = 20, sigma_r = 200, kernelsize = 5):
    row, col = img.shape
    nr_img = img.copy()
    bound = (kernelsize - 1) // 2
    kernel_s = gaussian_kernel(kernelsize, sigma_s)

    for i in range(bound, row - bound):
        for j in range(bound, col - bound):
            diff = np.abs(img[i-bound:i+bound+1, j-bound:j+bound+1] - img[i, j])**2
            kernel_r = np.exp(-diff / (2 * sigma_r ** 2))
            kernel = kernel_r * kernel_s / np.sum(kernel_r * kernel_s)
            nr_img[i, j] = np.sum(img[i-bound:i+bound+1, j-bound:j+bound+1] * kernel)

    return nr_img

可以看到双边滤波在额外添加了灰度差项后，图像边缘保留比高斯滤波更好，平坦区域降噪力度也差不多。

3.2 引导滤波

引导滤波 (Guided Image Filter) 是另一种保边滤波方法，双边滤波根据图像像素邻域的差值分配权重，有很好的保持边缘的效果，但是有时候也会在边缘产生伪影。因此作者思考利用引导图像（通常是原始图像或某个特征图像）的信息来指导目标图像的滤波过程。这种方法能有效保持边缘结构，尤其适合处理含有复杂细节和平滑区域混合的图像，同时也常用于深度图、光流估计等场景下的预处理。

参考论文：Guided Image Filtering

假定引导图I和原始无噪声图像q之间的关系是线性的，在局部窗口我们可以定义输出为：

可以看到在局部滤波窗中当I在边缘位置时，q也必然在边缘位置，输入噪声图像为p，那么：

因此通过最小化上述等式的差值，可以对线性系数进行拟合：

其中ε作为正则化系数可以控制添加项的大小以控制降噪程度，通过最小二乘法求解上式，可以得到系数的解：

其中μ和σ为引导图I在滤波窗中的均值及方差，|w|为在滤波窗中的像素个数

def guideFilter(I, p, winSize=(5,5), eps=1000):
    mean_I = cv2.blur(I, winSize)
    mean_p = cv2.blur(p, winSize)

    mean_II = cv2.blur(I * I, winSize)
    mean_Ip = cv2.blur(I * p, winSize)

    var_I = mean_II - mean_I * mean_I
    cov_Ip = mean_Ip - mean_I * mean_p

    a = cov_Ip / (var_I + eps)
    b = mean_p - a * mean_I

    mean_a = cv2.blur(a, winSize)
    mean_b = cv2.blur(b, winSize)

    q = mean_a * I + mean_b
    return q

如果假设滤波后的图像是引导图像的线性变换，那么导向滤波在一方面可以达到双边滤波的效果有效保持边缘，非迭代计算，但是不会有双边滤波的边缘伪影。另一方面可以使得输出结构化，不至于过于平滑。同时导向滤波的复杂度是O(N)， 而双边滤波在kernel radius为r的时候，复杂度是O(Nr^2)。

一般在降噪任务中，引导图I与输入图 p 为同一个图像。这时导向滤波的效果与双边滤波的效果类似，但是不同于双边滤波的是，导向滤波可以很容易设计一个与滤波半径无关的优化算法。其中窗口半径为平滑半径，参数为平滑项参数，其值越大平滑的越明显。

在一些mask精细化的任务如深度图估计或融合权重计算，当输入图p为一个初始的深度估计图像时，导向图I为对应的原始场景的图像，导向滤波会精细化深度估计图像的边缘准确性。

实际引导滤波在图像平坦区域和边缘区域都拥有相同平滑项参数，边缘附近产生光晕伪影。后续的研究中加权引导滤波、梯度域引导滤波进一步解决了这个问题：

参考论文：Weighted Guided Image Filtering

加权引导滤波（WGIF）引入了边缘感知项，代表像素邻域patch的方差与图像所有patch方差总和的比值：

上式带入GIF中的平滑项参数，代表不同区域根据边缘感知有不同的控制力度：

参考：Gradient Domain Guided Image Filtering（梯度域导向滤波） - ayew - 博客园 (cnblogs.com)

梯度域引导滤波（GGIF）对边缘感知项进行了改进，并加入了线性系数的控制进一步保持边缘：

3.3 变换域滤波

参考论文：Domain Transform for Edge-Aware Image and Video Processing

变换域滤波（DTFilter）相较于双边滤波在性能和保边性上都做了一定改进：首先提出了距离保持域变换的思想将高维的滤波降维到低维空间去做，如通过两个一维卷积代替二维卷积。然后对于变换域中的一维滤波器，设计了三种针对不同任务的卷积方案，在实时性和保边性都有不错的效果。

(a) 原始信号 （b）域变换 （c）信号域变换后 （d）对变换后的信号进行滤波 （e）逆域变换

首先是域变换部分，对于二维RGB图像的滤波可以表示如下：

其中 �=(��,��) 代表空间坐标， �(�)=(��,��,��) 代表像素值，那么 � 可以看成是一个五维的空间滤波器（当�设计为公式2的时候就是双边滤波器），计算成本比较高。如果存在一种图像降维变换可以在低维空间中进行滤波，同时这种变换不会破坏图像的像素间距离以维持保边性，那么计算效率会大幅提高：

这里 � 代表降维的域变换， � 代表新空间维度下的滤波核。通常像素距离（坐标值与像素值距离）越大，参与滤波的权重应该越低（如双边滤波器里的 �2 距离）。如果高维空间到低维空间的域变换后像素距离可以保持不变，那么低维空间滤波器的保边性也不会变：

约束域变换前后，邻域像素的距离恒定

令域变换 ��(�)=�(�,�(�)) ，像素距离用 �1 来表示，域变换表达式可以推导如下：

对于多通道输入，我们只需要把多个通道的像素距离叠加进来即可：

c代表通道数

为了像双边滤波器一样可以控制内核的滤波核形状，我们可以把尺度信息添加到这个域变换上作用在信号中，其实等价于对滤波核的操作：

这里省略部分推导，最终滤波器 � 响应的大小由 �� 、 �� 和梯度 �′ 共同控制，来完成保边强度的控制：

实际上对于高维空间到二维空间的变换，前面讨论的理想域变换并不存在，几乎没有域变换可以满足每个像素位置的距离不变要求。而针对二维的图像信号，文中采用了水平和垂直方向交错进行一维滤波的操作。

接下来是变换域中滤波器的设计，对于一维信号的处理文中针对不同的任务提出了三种不同的处理方案：

Normalized Convolution (NC)

归一化卷积将非均匀采样信号变换到等间距的采样信号并加权平均，在半径 � 的范围计算冲激响应并归一化：

Interpolated Convolution (IC)

插值卷积对离散不均匀信号进行插值，在连续信号上进行卷积：

Recursive Filtering (RF)

递归滤波器对于变换域中，根据 � 的数值控制边缘保留强度：

import cv2
import numpy as np

def NC(I, guide, sigma_s=60, sigma_r=0.4, num_iterations=3):
    h, w, c = I.shape
    if guide.shape != I.shape:
        return I

    #引导图J的多通道ct(u)域变换计算，用来判断保边强度
    J = guide
    dIcdx = np.zeros((h, w, c), dtype=np.float32)
    dIcdy = np.zeros((h, w, c), dtype=np.float32)
    dIdx = np.zeros((h, w), dtype=np.float32)
    dIdy = np.zeros((h, w), dtype=np.float32)

    for chn in range(c):
        dIcdx[:,1:,chn] = J[:,1:,chn] - J[:,:-1,chn]
        dIcdy[1:,:,chn] = J[1:,:,chn] - J[:-1,:,chn]
        dIdx[:,1:] = dIdx[:,1:] + np.abs(dIcdx[:,1:,chn])
        dIdy[1:,:] = dIdy[1:,:] + np.abs(dIcdy[1:,:,chn])
    dHdx = (1 + sigma_s / sigma_r * dIdx)
    dVdy = (1 + sigma_s / sigma_r * dIdy)

    #引导图J的梯度积分(公式11)
    ct_H = np.cumsum(dHdx, 1)
    ct_V = np.cumsum(dVdy, 0)
    ct_V = ct_V.T

    #Normalized Convolution
    sigma_H = sigma_s
    #对输入图I迭代N次1D滤波，垂直水平交替滤波可以消除artifact条纹
    for iters in range(num_iterations):
        print("iters: ", iters)
        #不同尺度核的标准差，随迭代进行衰减
        sigma_H_i = sigma_H * np.sqrt(3) * np.power(2, num_iterations - (iters + 1)) / np.sqrt(np.power(4, num_iterations) - 1)
        box_radius = np.sqrt(3) * sigma_H_i
        print("current box_radius: ", box_radius)
        #行滤波
        I = TransformedDomainBoxFilter_Horizontal(I, ct_H, box_radius)
        I = image_transpose(I)
        #列滤波
        I = TransformedDomainBoxFilter_Horizontal(I, ct_V, box_radius)
        I = image_transpose(I)
        cv2.imwrite("res_sigma_s{}sigma_r{}iters{}.jpg".format(sigma_s,sigma_r,iters), I * 255)
    return I


def TransformedDomainBoxFilter_Horizontal(I, ct_x, box_radius):
    h, w, c = I.shape
    l_pos = ct_x - box_radius
    u_pos = ct_x + box_radius
    l_idx = np.zeros((h, w), dtype=np.int64)
    u_idx = np.zeros((h, w), dtype=np.int64)

    for row in range(1,h):
        ct_x_row = ct_x[row, :]
        l_pos_row = l_pos[row, :]
        u_pos_row = u_pos[row, :]
        local_l_idx = np.zeros((1, w), dtype=np.int64)
        local_u_idx = np.zeros((1, w), dtype=np.int64)
        index_l = np.where(ct_x_row > l_pos_row[1])
        index_u = np.where(ct_x_row > u_pos_row[1])
        local_l_idx[0][1] = index_l[0][0]
        local_u_idx[0][1] = index_u[0][0]
        for col in range(2, w):
            index_l = np.where(ct_x_row[local_l_idx[0][col-1]:] > l_pos_row[col])
            index_u = np.where(ct_x_row[local_u_idx[0][col-1]:] > u_pos_row[col])
            if len(index_l[0]) == 0:
                local_l_idx[0][col] = local_l_idx[0][col - 1]
            else:
                local_l_idx[0][col] = local_l_idx[0][col - 1] + index_l[0][0]
            if len(index_u[0]) == 0:
                local_u_idx[0][col] = local_u_idx[0][col - 1]
            else:
                local_u_idx[0][col] = local_u_idx[0][col - 1] + index_u[0][0]
        l_idx[row,:] = local_l_idx
        u_idx[row,:] = local_u_idx

    SAT = np.zeros((h, w, c), dtype=np.float32)
    F = np.zeros((h, w, c), dtype=np.float32)
    SAT = np.cumsum(I, 1)

    for chn in range(c):
        for row in range(h):
            for col in range(w):
                sat_a = SAT[row][l_idx[row][col]][chn]
                sat_b = SAT[row][u_idx[row][col]][chn]
                F[row][col][chn] = (sat_b - sat_a) / (u_idx[row][col] - l_idx[row][col]+0.0001)
    return F


def image_transpose(img):
    h, w, c = img.shape
    img_T = np.zeros((w, h, c), dtype = np.float64)
    for chn in range(c):
        img_T[:, :, chn] = img[:, :, chn].T
    return img_T


if __name__ == "__main__":
    img = cv2.imread("statue.png")
    img = img.astype(np.float32) / 255
    sigma_s = 60
    sigma_r = 0.4
    method = 1

    if method == 1:
        res = NC(img, img, sigma_s, sigma_r)
    elif method == 2:
        res = IC(img, img, sigma_s, sigma_r)
    else:
        res = RF(img, img, sigma_s, sigma_r)

---

4. Non-local mean（NLM）

非局部均值（Non-Local Means）算法是一种基于自相似性的图像去噪方法，打破了传统局部滤波器只能考虑目标像素及其邻域内像素的局限，而是全局地比较并整合整个图像中的像素点以达到降噪的目的。

NLM算法的核心思想是：噪声在整幅图像中是随机分布的，而图像的纹理和结构信息往往在多个位置有重复表现。因此通过在全局范围内寻找相似性patch，并对其进行加权平均，噪声就会在大量相似像素加权的过程中稀释掉。而真实信号由于在很多地方都有重复出现，所以得以保留和强化。

具体实现流程为：图像中的每个像素与其在整幅图像中所有其他像素的相关性进行比较，如果两像素对应的邻域像素块之间足够相似，则赋予较高的权重。然后用这些相似像素的灰度值加权平均来估计目标像素的真实值，从而降低噪声影响。

参考论文：Non-Local Means Denoising

在每个搜索窗中，会计算每个像素和中心像素的欧几里得距离：

距离越大代表越不相似，给予加权的权重就越小:

def nlm(img, sigma=50, winSize1=3, winSize2=7):
    bound1 = (winSize1 - 1) // 2
    bound2 = (winSize2 - 1) // 2
    row, col = img.shape
    nr_img = img.copy()

    for i in range(bound1, row-bound1):
        for j in range(bound1, col-bound1):
            cur_patch = img[i-bound1:i+bound1+1,j-bound1:j+bound1+1]
            sum_val = 0
            sum_weight = 0
            for m in range(bound2, winSize1-bound2):
                for n in range(bound2, winSize1-bound2):
                    diff = np.sum((cur_patch[m-bound2:m+bound2+1, n-bound2:n+bound2+1] - img[i-bound2:i+bound2+1, j-bound2:j+bound2+1])**2)/(winSize2**2)
                    weight = np.exp(-diff/sigma)
                    sum_val += weight * cur_patch[m, n]
                    sum_weight += weight
            nr_img[i, j] = sum_val / (sum_weight + 0.001)
    return nr_img

NLM原始算法的复杂度还是相当高，导致性能开销较大，针对NLM算法也有相当多性能优化方向的研究：

1. 限制搜索范围或采样搜索，减少比较次数。
2. 简化权重计算公式，或使用如快速最近邻搜索算法（如kd树、BBF树）查找相似像素。
3. 使用积分图像加速以及circle shift。
4. 多尺度处理，从粗略尺度开始，逐步细化至精确尺度。
5. 利用GPU并行计算能力加速处理。