Python 在科研中的应用 08：Python 环境下的数字图像降噪

Antoni Buades 提出指标 method noise 对数字图像降噪方法的性能进行了评价和比较。他首先针对几个被广泛使用的降噪算法计算并分析了降噪性能。同时，基于图像中所有像素的非局部平均，提出了全新的数字图像降噪算法 Non Local means Algorithm，并通过实验比较了新算法与常用的平滑滤波方法的性能。

Introduction

课程作业 01，占总成绩30%

在二维平面中构建随机 Vornoi 结构，将交界处壁厚设置为随机数，生成Vornoi结构图像。要求：三维平面的像素画幅、Vornoi的初始种子的数密度、交界处壁厚的随机分布范围可设置，作为输入参数集放置在程序的最开始。输出：二维结构图像，并输出各个交界节点、棱边的结构信息至.txt文件。

数字图像噪声

  数字图像的两个主要限制是模糊及噪声（blur & noise）。模糊是图像采集系统的固有性质，同时数字图像对连续信号离散采样的形式必须遵循 Shannon-Nyquist 采样定律。另一种主要的干扰形式是噪声。数字图像中每个像素点的值 $u(i)$ 都是对局部光强测量的结果，通常通过 CCD 及光学聚焦元件实现。CCD 中每个方形探测单元（captor）都将记录曝光时间内探测区域的入射光子数。在光源强度恒定的情况下，每个探测单元每个曝光周期内采集的光子数的概率分布将遵循中心极限定理，在光强均值周围震荡。另外，如果 CCD 没有经过充分冷却就会接收热光子（heat spurious photons），由此产生的扰动通常称为 obscurity noise。

  数字图像降噪方法的目标是从观测到的噪声图像中还原出原始信号，

$$ v(i)=u(i)+n(i), \tag{1}\label{1} $$

其中 $v(i)$ 为实测图像，$u(i)$ 为原始信号，$n(i)$ 则是噪声信号。评估图像中噪声水平通常采用信噪比（signal noise ratio, SNR）：

$$ SNR = \sigma(u)/\sigma(n), \tag{2}\label{2} $$

其中，$\sigma(n)$为噪声信号标准差，$\sigma(u)$表示真实信号的经验标准差，

$$ \sigma(u) = \sqrt{\frac{1}{|I|}\sum_i(u(i)- \overline{u} )^2}, \tag{3}\label{3} $$

$$ \overline{u}=\frac{1}{|I|}\sum_{i\in I} u(i) \tag{4}\label{4}$$ 为图像的平均灰度值，$|I|$指全图像素数。当噪声模型和参数已知时，噪声的标准差也可以用经验测量法或形式化方法得到。

  迄今为止，图像处理领域已经提出了诸多抑制噪声、还原真实信号的算法。即便它们通常拥有截然不同的数学形式，但都拥有一个共性：平均。这种平均可以在局部进行：高斯滤波器(Gabor 1994)，各向异性滤波(Perona-Malik 1990, Alvarez et al. 1992)，邻域滤波器(Yaroslavsky 1985, Smith et al. 1997)；也可以通过计算variations实现：TV滤波(Rudin-Osher-Fatemi 1992)；或是在频域进行：经验维纳滤波(Yaroslavsky 1985)，小波阈值方法(Coiffman-Donoho 1995)。

Method noise

  不妨令 $u$ 表示实测图像，$D_hu$ 表示降噪方法的输出结果，$h$ 为滤波参数。Antoni Buades 定义 method noise 为降噪前后图像之差：

$$ u-D_hu. \tag{5}\label{5} $$

  完美的降噪算法在应用中不应该改变无噪声图像。因此，当图像具有某种规律性时，method noise 理应很小。对理想的降噪算法，Method noise 必须看起来与随机噪声无异，几乎不包含原始信号的结构。因为即便是质量非常高的实测图像，噪声也是不可避免的，计算 method noise 对评估任何降噪算法都是有意义的，而非传统的“添加噪声，再去除噪声”的把戏。

局部平均算法

高斯滤波 (Gaussian Filtering)

  对数字图像进行各向同性过滤，本质上可以归结为图像与各向同性核的卷积。采用数值呈现出高斯分布的卷积核，既是高斯滤波，是图像处理中最常用的操作之一。通俗的讲，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到（所有的局部平滑滤波方法都是如此）。高斯卷积核：

$$ G_h(x)=\frac{1}{4\pi h^2}e^{-\frac{|x|^2}{4h^2}}. \tag{6}\label{6} $$

Theorem 1 (Babor 1960): 当高斯卷积核的特征尺寸 $h$ 极小时，高斯滤波的 method noise 为：

$$ u-G_h* u=-h^2\Delta u+o(h^2). \tag{7}\label{7} $$

  高斯滤波的 method noise 在图像谐波部分几乎为零，而在边缘、纹理区域非常大。因此，高斯滤波在图像的平坦区域相对表现优秀，但在边缘、纹理区域较为模糊（blurred）。

import cv2
import numpy as np

img_medianBlur=cv2.medianBlur(img01,5)
font=cv2.FONT_HERSHEY_SIMPLEX

#均值滤波
img_Blur=cv2.blur(img01,(5,5))

#高斯滤波
img_GaussianBlur=cv2.GaussianBlur(img01,(7,7),0)

#高斯双边滤波
img_bilateralFilter=cv2.bilateralFilter(img01,40,75,75)

各向异性滤波（Anisotropic Filtering, AF）

  各向异性滤波提出之初旨在解决高斯滤波在边缘及纹理区域的模糊问题。该算法通过只在 $Du(\vec x)$ 正交方向上计算图像 $u$ 在 $\vec x$ 处的卷积。这种想法可以追溯到 Perona & Malik。各向异性滤波算法的定义为：

$$ AF_hu(\vec x)=\int G_h(t)u(\vec x+t\frac{Du(\vec x)^\perp}{|Du(\vec x)|})dt, \tag{8}\label{8} $$

在 $\vec x$ 处，当 $Du(\vec x)\neq0$ 时成立。$(x,y)^\perp=(-y,x)$，且 $G_h$ 代指方差 $h^2$ 的一维高斯函数。假设原始图像 $u$ 在 $\vec x$ 处二次连续可微（twice continuously differentiable ($C^2$)）,将上式二次泰勒展开（second order Taylor expansion）可以推导出：

Theorem 2: 当 $Du(\vec x)\neq0$ 时，各向异性滤波 $AF_h$ 的 method noise 为：

$$ u(\vec x)-AF_hu(\vec x)=-\frac{1}{2}h^2|Du|curv(u)(\vec x)+o(h^2), \tag{9}\label{9} $$

$curv(u)(\vec x)$ 指局部曲率，即经过 $\vec x$ 点的水平直线上曲率半径的逆（signed inverse）。在图像 $u$ 中局部几乎为一条直线的区域，method noise 几乎为零，而在弯曲的边缘或纹理区域较大。因而，各向异性滤波对直边区域得以保持，而平坦、纹理区域图像精度有所退化。

import numpy as np
import warnings

def anisodiff(img,niter=1,kappa=50,gamma=0.1,step=(1.,1.),option=1,ploton=False):
    """
    Anisotropic diffusion.

    Usage:
    imgout = anisodiff(im, niter, kappa, gamma, option)

    Arguments:
            img    - input image
            niter  - number of iterations
            kappa  - conduction coefficient 20-100 ?
            gamma  - max value of .25 for stability
            step   - tuple, the distance between adjacent pixels in (y,x)
            option - 1 Perona Malik diffusion equation No 1
                     2 Perona Malik diffusion equation No 2
            ploton - if True, the image will be plotted on every iteration

    Returns:
            imgout   - diffused image.

    kappa controls conduction as a function of gradient.  If kappa is low
    small intensity gradients are able to block conduction and hence diffusion
    across step edges.  A large value reduces the influence of intensity
    gradients on conduction.

    gamma controls speed of diffusion (you usually want it at a maximum of
    0.25)

    step is used to scale the gradients in case the spacing between adjacent
    pixels differs in the x and y axes

    Diffusion equation 1 favours high contrast edges over low contrast ones.
    Diffusion equation 2 favours wide regions over smaller ones.

    Reference:
    P. Perona and J. Malik.
    Scale-space and edge detection using ansotropic diffusion.
    IEEE Transactions on Pattern Analysis and Machine Intelligence,
    12(7):629-639, July 1990.

    Original MATLAB code by Peter Kovesi
    School of Computer Science & Software Engineering
    The University of Western Australia
    pk @ csse uwa edu au
    <http://www.csse.uwa.edu.au>

    Translated to Python and optimised by Alistair Muldal
    Department of Pharmacology
    University of Oxford
    <alistair.muldal@pharm.ox.ac.uk>

    June 2000  original version.
    March 2002 corrected diffusion eqn No 2.
    July 2012 translated to Python

    April 2019 - Corrected for Python 3.7 - AvW 
    """

    # ...you could always diffuse each color channel independently if you
    # really want
    if img.ndim == 3:
        warnings.warn("Only grayscale images allowed, converting to 2D matrix")
        img = img.mean(2)

    # initialize output array
    img = img.astype('float32')
    imgout = img.copy()

    # initialize some internal variables
    deltaS = np.zeros_like(imgout)
    deltaE = deltaS.copy()
    NS = deltaS.copy()
    EW = deltaS.copy()
    gS = np.ones_like(imgout)
    gE = gS.copy()

    # create the plot figure, if requested
    if ploton:
        import pylab as pl
        from time import sleep

        fig = pl.figure(figsize=(20,5.5),num="Anisotropic diffusion")
        ax1,ax2 = fig.add_subplot(1,2,1),fig.add_subplot(1,2,2)

        ax1.imshow(img,interpolation='nearest')
        ih = ax2.imshow(imgout,interpolation='nearest',animated=True)
        ax1.set_title("Original image")
        ax2.set_title("Iteration 0")

        fig.canvas.draw()

    for ii in range(niter):

        # calculate the diffs
        deltaS[:-1,: ] = np.diff(imgout,axis=0)
        deltaE[: ,:-1] = np.diff(imgout,axis=1)

        # conduction gradients (only need to compute one per dim!)
        if option == 1:
            gS = np.exp(-(deltaS/kappa)**2.)/step[0]
            gE = np.exp(-(deltaE/kappa)**2.)/step[1]
        elif option == 2:
            gS = 1./(1.+(deltaS/kappa)**2.)/step[0]
            gE = 1./(1.+(deltaE/kappa)**2.)/step[1]

        # update matrices
        E = gE*deltaE
        S = gS*deltaS

        # subtract a copy that has been shifted 'North/West' by one
        # pixel. don't as questions. just do it. trust me.
        NS[:] = S
        EW[:] = E
        NS[1:,:] -= S[:-1,:]
        EW[:,1:] -= E[:,:-1]

        # update the image
        imgout += gamma*(NS+EW)

        if ploton:
            iterstring = "Iteration %i" %(ii+1)
            ih.set_data(imgout)
            ax2.set_title(iterstring)
            fig.canvas.draw()
            # sleep(0.01)

    return imgout

def anisodiff3(stack,niter=1,kappa=50,gamma=0.1,step=(1.,1.,1.),option=1,ploton=False):
    """
    3D Anisotropic diffusion.

    Usage:
    stackout = anisodiff(stack, niter, kappa, gamma, option)

    Arguments:
            stack  - input stack
            niter  - number of iterations
            kappa  - conduction coefficient 20-100 ?
            gamma  - max value of .25 for stability
            step   - tuple, the distance between adjacent pixels in (z,y,x)
            option - 1 Perona Malik diffusion equation No 1
                     2 Perona Malik diffusion equation No 2
            ploton - if True, the middle z-plane will be plotted on every
                 iteration

    Returns:
            stackout   - diffused stack.

    kappa controls conduction as a function of gradient.  If kappa is low
    small intensity gradients are able to block conduction and hence diffusion
    across step edges.  A large value reduces the influence of intensity
    gradients on conduction.

    gamma controls speed of diffusion (you usually want it at a maximum of
    0.25)

    step is used to scale the gradients in case the spacing between adjacent
    pixels differs in the x,y and/or z axes

    Diffusion equation 1 favours high contrast edges over low contrast ones.
    Diffusion equation 2 favours wide regions over smaller ones.

    Reference:
    P. Perona and J. Malik.
    Scale-space and edge detection using ansotropic diffusion.
    IEEE Transactions on Pattern Analysis and Machine Intelligence,
    12(7):629-639, July 1990.

    Original MATLAB code by Peter Kovesi
    School of Computer Science & Software Engineering
    The University of Western Australia
    pk @ csse uwa edu au
    <http://www.csse.uwa.edu.au>

    Translated to Python and optimised by Alistair Muldal
    Department of Pharmacology
    University of Oxford
    <alistair.muldal@pharm.ox.ac.uk>

    June 2000  original version.
    March 2002 corrected diffusion eqn No 2.
    July 2012 translated to Python
    """

    # ...you could always diffuse each color channel independently if you
    # really want
    if stack.ndim == 4:
        warnings.warn("Only grayscale stacks allowed, converting to 3D matrix")
        stack = stack.mean(3)

    # initialize output array
    stack = stack.astype('float32')
    stackout = stack.copy()

    # initialize some internal variables
    deltaS = np.zeros_like(stackout)
    deltaE = deltaS.copy()
    deltaD = deltaS.copy()
    NS = deltaS.copy()
    EW = deltaS.copy()
    UD = deltaS.copy()
    gS = np.ones_like(stackout)
    gE = gS.copy()
    gD = gS.copy()

    # create the plot figure, if requested
    if ploton:
        import pylab as pl
        from time import sleep

        showplane = stack.shape[0]//2

        fig = pl.figure(figsize=(20,5.5),num="Anisotropic diffusion")
        ax1,ax2 = fig.add_subplot(1,2,1),fig.add_subplot(1,2,2)

        ax1.imshow(stack[showplane,...].squeeze(),interpolation='nearest')
        ih = ax2.imshow(stackout[showplane,...].squeeze(),interpolation='nearest',animated=True)
        ax1.set_title("Original stack (Z = %i)" %showplane)
        ax2.set_title("Iteration 0")

        fig.canvas.draw()

    for ii in range(niter):

        # calculate the diffs
        deltaD[:-1,: ,:  ] = np.diff(stackout,axis=0)
        deltaS[:  ,:-1,: ] = np.diff(stackout,axis=1)
        deltaE[:  ,: ,:-1] = np.diff(stackout,axis=2)

        # conduction gradients (only need to compute one per dim!)
        if option == 1:
            gD = np.exp(-(deltaD/kappa)**2.)/step[0]
            gS = np.exp(-(deltaS/kappa)**2.)/step[1]
            gE = np.exp(-(deltaE/kappa)**2.)/step[2]
        elif option == 2:
            gD = 1./(1.+(deltaD/kappa)**2.)/step[0]
            gS = 1./(1.+(deltaS/kappa)**2.)/step[1]
            gE = 1./(1.+(deltaE/kappa)**2.)/step[2]

        # update matrices
        D = gD*deltaD
        E = gE*deltaE
        S = gS*deltaS

        # subtract a copy that has been shifted 'Up/North/West' by one
        # pixel. don't as questions. just do it. trust me.
        UD[:] = D
        NS[:] = S
        EW[:] = E
        UD[1:,: ,: ] -= D[:-1,:  ,:  ]
        NS[: ,1:,: ] -= S[:  ,:-1,:  ]
        EW[: ,: ,1:] -= E[:  ,:  ,:-1]

        # update the image
        stackout += gamma*(UD+NS+EW)

        if ploton:
            iterstring = "Iteration %i" %(ii+1)
            ih.set_data(stackout[showplane,...].squeeze())
            ax2.set_title(iterstring)
            fig.canvas.draw()
            # sleep(0.01)

    return stackout

TV滤波（Total Variation Minimization）

  全变分图像去噪算法最早由 Rudin, Osher and Fatemi提出。给定一幅实测图像 $v(\vec x)$，该算法将恢复原始信号 $u(\vec x)$ 的问题转化为下式的最小化问题：

$$ TVF_\lambda(v)={\rm arg},\underset{u}{\rm \mathop {min}},TV(u)+\lambda\int|v(\vec x)-u(\vec x)|^2d\vec x, \tag{10}\label{10} $$

其中 $TV(u)$ 为图像 $u$ 的全变分，$\lambda$ 为给定的拉格朗日乘子（Lagrange multiplier）。上述最小化问题的最小值存在且唯一。参数 $\lambda$ 与噪声的统计信息相关，并控制了滤波程度。

Theorem 3: TV滤波的 method noise 为：

$$ u(\vec x)-TVF_\lambda(u)(\vec x)=-\frac{1}{2\lambda}curv(TVF_\lambda(u))(\vec x). \tag{11}\label{11} $$

  在各向异性的情况下，直边由于曲率小而得以保持。但 $\lambda$ 过小时细节及纹理会被过度平滑。

邻域滤波（Neighborhood Filtering）

  我们称邻域滤波为将临近区域内具有相似灰度值的像素取平均以期恢复原始信号的滤波器。Yaroslavsky(1985)首次提出的方法通过计算空间邻域 $B_\rho(\vec x)$ 内具有相似灰度值像素的平均来恢复信号：

$$ YNF_{h,\rho}u(\vec x)=\frac{1}{C(\vec x)}\int_{B_\rho (\vec x)} u(\vec y)e^{-\frac{|u(\vec y)-u(\vec x)|^2}{h^2}}d\vec y, \tag{12}\label{12} $$

其中，$\vec x\in\Omega$，$C(\vec x)=\int_{B_\rho(\vec x)} e^{-\frac{|u(\vec y)-u(\vec x)|^2}{h^2}}d\vec y$ 为归一化常量，$h$ 为滤波参数。

  较晚提出的邻域滤波算法相较 Yaroslavsky 滤波要更广为人知一些，既 SUSAN 滤波（1995）及双边滤波（1998）。这些滤波算法都引入到参考像素 $\vec x$ 的距离作为权重因子，而非单纯的考虑一个固定范围的邻域，

$$ SNF_{h,\rho}u(\vec x)=\frac{1}{C(\vec x)}\int_{\Omega} u(\vec y)e^{-\frac{|y-x|^2}{\rho^2}}e^{-\frac{|u(\vec y)-u(\vec x)|^2}{h^2}}d\vec y, \tag{13}\label{13} $$

这里 $C(\vec x)=\int_{\Omega}e^{-\frac{|y-x|^2}{\rho^2}}e^{-\frac{|u(\vec y)-u(\vec x)|^2}{h^2}}d\vec y$ 为归一化常量，$\rho$ 为空间滤波参数（spatial filtering parameter）。事实上，$YNF_{h,\rho}$ 与 $SNF_{h,\rho}$ 之间并没有本质区别。如果两个区域的灰度值差异大于 $h$，这些算法都将计算来自于同一区域的像素灰度平均值来恢复参考点的原始信号。因而该算法不会模糊边界区域，这是正是该类算法最核心的用途。

  然而这类算法的问题是，只将单个像素作为参考点，而如若该参考像素恰好被噪声干扰严重，滤波效果将不够鲁棒（robust）。同时，邻域滤波器也会制造人为干扰（artificial shocks），这将会在它的 method noise 中展示出来。

非局部平均算法（Non Local Means Algorithm）

  Antoni Buades 于 2005 年提出非局部平均数字图像降噪算法（Non Local Algorithm）。给定一幅实测图像 $v=\lbrace v(i)|i\in I\rbrace $，像素 $i$ 处的估计值 $NL[v](i)$ 是该图像上所有像素点的加权平均值，

$$ NL(i)=\sum_{j\in I}\omega(i,j)v(j), \tag{14}\label{14} $$

这里的权重系数 $\lbrace\omega (i,j)\rbrace_j$ 取决于像素 $i$ 与 $j$ 之间的相似程度，且始终满足如下标准： $0\leq\omega(i,j)\leq 1$ 且 $\sum_{j} \omega(i,j)=1$（等价于上文介绍的滤波算法中归一化常量）。

  俩个像素 $(i,j)$ 之间的相似程度取决于邻域灰度矩阵 $v(N_i)$ 及 $v(N_j)$（intensity gray level vectors），这里 $N_k$ 指以像素 $k$ 为中心，给定尺寸的方形邻域。这种相似性被定义为加权欧式距离的递减函数，$\parallel v(N_i)-v(N_j)\parallel_{2,a}^2$，其中 $a>0$ 是高斯卷积核的标准差。欧式距离对噪声邻域的应用引入了如下等式：

$$ E\parallel v(N_i)-v(N_j)\parallel_{2,a}^2 = \parallel u(N_i)-u(N_j) \parallel_{2,a}^2+2\sigma^2. \tag{15}\label{15} $$

  这个等式证明了该算法的鲁棒性（robustness）。因为含噪声的实测图像 $v$ 的欧式距离期望恰恰遵循真正的原始信号之间的相似性。

  与 $v(N_i)$ 具有相似灰度邻域的像素在计算平均时的权重因子更大，

$$ \omega(i,j)=\frac{1}{Z(i)}e^{-\cfrac{\parallel v(N_i)-v(N_j) \parallel_{2,a}^2}{h^2}}, \tag{16}\label{16} $$

这里，$Z(i)$ 为归一化常量，

$$ Z(i) = \sum_j e^{-\cfrac{\parallel v(N_i)-v(N_j) \parallel_{2,a}^2}{h^2}}, \tag{17}\label{17} $$

参数 $h$ 控制滤波程度，它直接影响了指数函数的衰减趋势，进而控制欧式距离对权重因子衰减速度的影响。

  NL-means 算法不仅仅考虑单个像素的灰度值，而是考虑该像素整个邻域的几何构型，这正是 NL-means 算法比邻域滤波更鲁棒的原因。图$(1)$ 也说明了这个问题，像素 $q3$ 与 $p$ 具有完全一致的灰度值，而邻域的几何构型完全不同，导致 NL-means 中的权重因子 $\omega(p,q3)$ 几乎为零。

 
  NL-means 算法最终的数学形式：

$$ NL(x)=\frac{1}{C(x)} \int_{\Omega} e^{-\cfrac{(G_a * |v(x+.)-v(y+.)|^2)(0)}{h^2}}v(y)dy, \tag{18}\label{18} $$

$x\in\Omega$，$C(x)=\int_{\Omega} {\rm exp}\lbrack -\frac{(G_a* |u(x+.)-u(z+.)|^2)(0)}{h^2}\rbrack dz$ 为归一化常量，$G_a$ 为高斯核，$h$ 控制过滤程度。

  NL-means 算法的中心思想是：像素 $x$ 处的信息恢复，是由整幅图像内所有邻域与像素 $x$ 邻域相似的点取平均得到的。与局部滤波算法或频域滤波算法相比，NL-means 算法的主要区别在于可以系统地运用整幅图像中所有可能自预测局部结构的信息。

Non-local means 方法代码实现

import numpy as np
import matplotlib.pyplot as plt

from skimage import data, img_as_float
from skimage.restoration import denoise_nl_means, estimate_sigma
from skimage.metrics import peak_signal_noise_ratio
from skimage.util import random_noise


astro = img_as_float(data.astronaut())
astro = astro[30:180, 150:300]

sigma = 0.08
noisy = random_noise(astro, var=sigma**2)

# estimate the noise standard deviation from the noisy image
sigma_est = np.mean(estimate_sigma(noisy, channel_axis=-1))
print(f'estimated noise standard deviation = {sigma_est}')

patch_kw = dict(patch_size=5,      # 5x5 patches
                patch_distance=6,  # 13x13 search area
                channel_axis=-1)

# slow algorithm
denoise = denoise_nl_means(noisy, h=1.15 * sigma_est, fast_mode=False,**patch_kw)

# slow algorithm, sigma provided
denoise2 = denoise_nl_means(noisy, h=0.8 * sigma_est, sigma=sigma_est,fast_mode=False, **patch_kw)

# fast algorithm
denoise_fast = denoise_nl_means(noisy, h=0.8 * sigma_est, fast_mode=True,**patch_kw)

# fast algorithm, sigma provided
denoise2_fast = denoise_nl_means(noisy, h=0.6 * sigma_est, sigma=sigma_est,fast_mode=True, **patch_kw)

fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 6),sharex=True, sharey=True)

ax[0, 0].imshow(noisy)
ax[0, 0].axis('off')
ax[0, 0].set_title('noisy')
ax[0, 1].imshow(denoise)
ax[0, 1].axis('off')
ax[0, 1].set_title('non-local means\n(slow)')
ax[0, 2].imshow(denoise2)
ax[0, 2].axis('off')
ax[0, 2].set_title('non-local means\n(slow, using $\\sigma_{est}$)')
ax[1, 0].imshow(astro)
ax[1, 0].axis('off')
ax[1, 0].set_title('original\n(noise free)')
ax[1, 1].imshow(denoise_fast)
ax[1, 1].axis('off')
ax[1, 1].set_title('non-local means\n(fast)')
ax[1, 2].imshow(denoise2_fast)
ax[1, 2].axis('off')
ax[1, 2].set_title('non-local means\n(fast, using $\\sigma_{est}$)')

fig.tight_layout()

# print PSNR metric for each case
psnr_noisy = peak_signal_noise_ratio(astro, noisy)
psnr = peak_signal_noise_ratio(astro, denoise)
psnr2 = peak_signal_noise_ratio(astro, denoise2)
psnr_fast = peak_signal_noise_ratio(astro, denoise_fast)
psnr2_fast = peak_signal_noise_ratio(astro, denoise2_fast)

print(f'PSNR (noisy) = {psnr_noisy:0.2f}')
print(f'PSNR (slow) = {psnr:0.2f}')
print(f'PSNR (slow, using sigma) = {psnr2:0.2f}')
print(f'PSNR (fast) = {psnr_fast:0.2f}')
print(f'PSNR (fast, using sigma) = {psnr2_fast:0.2f}')

plt.show()

实验及讨论

  本节将针对以下三点比较 Non Local Means 算法与局部平滑滤波的性能： method noise，视觉效果，均方差（mean square error）。这里均方差指修复图像与真实图像的欧几里得差异（Euclidean difference）。

  在实现 NL-means 算法时，我们将相似邻域的搜索范围约束在一个较大的 $S\times S$ pixels 区域内。在所有的实验中均固定该搜索范围为 $21\times 21$ pixels，并指定邻域 $N_i$ 范围为 $7\times 7$ pixels。不妨设输入图像的像素数为 $N^2$，那么该算法的时间复杂度约为 $21^2\times 7^2 \times N^2$。

<img src=”https://s21.ax1x.com/2024/05/07/pkEWwPH.jpg“ width=”70%” alt=”” align=center / title=”图2. NL-means 权重分布。每组中左图中心像素为参考像素点 $x$，每组中右图为权重因子分布，由黑至白对应从0到1。”/>

 
  $7\times 7$ pixels 的邻域范围已经足够大，能够确保算法对噪声是鲁棒的；也足够小，让细节及纹理得以保持。滤波参数 $h$ 被设置为 $10\sigma$，这里 $\sigma$ 为人工添加高斯白噪声的标准差。由于指数权重的快速衰减，距离中心较远像素的权重几乎为零，这发挥了自动剔除远处像素的作用。Non Local means 算法的权重分布见图$(2)$。

<img src=”https://s21.ax1x.com/2024/05/07/pkEW0Gd.jpg“ width = 80% div align=center / title=”图3. 对自然图像的降噪实验。从左至右：含噪声的实测图像（噪声标准差 35），高斯滤波，各向异性滤波，TV滤波，邻域滤波及 NL-means 算法。”>

 
  在前面的段落中，我们明确的计算了各种局部平滑滤波方法的 method noise 理论值。图$(4)$ 的可视化实验证明了前文中公式的正确性。图$(4)$ 对比了数种降噪方法对含噪声lena图像运算的 method noise，噪声为高斯白噪声，标准差 2.5，滤波参数 $h$ 均相同。Method noise 能更好的协助判断降噪算法的性能及局限，因为被移除或改变的纹理、细节将被 method noise 醒目的展示出来。对比图$(4)$ 中数种降噪算法，NL-means 算法的 method noise 几乎无法察觉任何几何纹理。图$(2)$ 可以解释这一现象，因为 NL-means 算法选取的权重因子完全适应图像局部及非局部的几何结构。

<img src=”https://s21.ax1x.com/2024/05/07/pkEWysP.jpg“ width = 55% div align=center / title=”图4. 降噪算法的method noise。从左到右，由上至下：噪声图像（标准差 20），高斯滤波，各向异性滤波，TV滤波，邻域滤波及 NL-means 算法。视觉实验验证了第二节的计算公式。”>

 
  由于算法本身的性质，纹理及周期性结构是 NL-means 算法最适用的情况。因为对任意像素 $i$，纹理图像或周期性图像中将可以找到大量与该像素具有相似邻域的点，如图$(2.e)$。图$(3)$ 展示了局部平滑滤波器及 NL-means 算法对自然纹理的平滑效果。

  自然图像同样含有足够的信息冗余会被 NL-means 恢复。平坦区域内部将呈现出大量相似的几何结构，见图$(2.a)$。直边界、或弯曲边界将被筛选出一条结构相似的像素线，见图$(2. b)，(2.c)$。并且，NL-means 将会在很远的位置寻找到与参考点相似的结构，如图$(2.f)$。图$(5)$ 展示了一次对自然图像的可视化实验，这组结果与图$(4)$ 是相对应的。

<img src=”https://s21.ax1x.com/2024/05/07/pkEWBRA.jpg“ width = 55% div align=center / title=”图5. 对自然图像的降噪实验。从左到右，由上至下：噪声图像（标准差 20），高斯滤波，各向异性滤波，TV滤波，邻域滤波及 NL-means 算法。”>

 
  最后，表$(1)$ 展示了本文介绍的降噪方法的均方差。这种数值测量方法是最客观的，因为它不依赖于任何肉眼视觉上的解释。然而，这个误差在实际问题中是不可计算的（原始信号是未知的），小的均方误差并不能保证高的视觉质量。因此，以上讨论的标准似乎是比较算法性能的必要条件。

<img src=”https://s21.ax1x.com/2024/05/07/pkEWsMt.jpg“ width = 35% div align=center / title=”表1. 均方误差表。均方误差越小，降噪后越接近原始图像。”>

 

 

Code and Data

  Non-local means 方法 C 语言实现：http://www.ipol.im/pub/art/2011/bcm_nlm/

其他滤波器

Top-hat filter

使用Top-hat滤波器去除灰度图像中的小目标：这个例子展示了如何从灰度图像中移除小物体。顶帽变换是一种从给定图像中提取小元素和细节的操作。这里我们使用白色顶帽变换，它被定义为输入图像与其(数学形态学)开口之间的差异。

import matplotlib.pyplot as plt

from skimage import data
from skimage import color, morphology

image = color.rgb2gray(data.hubble_deep_field())[:500, :500]

footprint = morphology.disk(1)
res = morphology.white_tophat(image, footprint)

fig, ax = plt.subplots(ncols=3, figsize=(20, 8))
ax[0].set_title('Original')
ax[0].imshow(image, cmap='gray')
ax[1].set_title('White tophat')
ax[1].imshow(res, cmap='gray')
ax[2].set_title('Complementary')
ax[2].imshow(image - res, cmap='gray')

plt.show()

小波傅里叶变换

小波去噪依赖于图像的小波表示。高斯噪声倾向于用小波域中的小值来表示，并且可以通过将低于给定阈值的系数设置为零(硬阈值)或将所有系数缩小到给定数量的零(软阈值)来去除。
在这个例子中，我们展示了两种不同的小波系数阈值选择方法:BayesShrink和VisuShrink。

import matplotlib.pyplot as plt

from skimage.restoration import denoise_wavelet, estimate_sigma
from skimage import data, img_as_float
from skimage.util import random_noise
from skimage.metrics import peak_signal_noise_ratio


original = img_as_float(data.chelsea()[100:250, 50:300])

sigma = 0.12
noisy = random_noise(original, var=sigma**2)

fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5), sharex=True, sharey=True)

plt.gray()

# Estimate the average noise standard deviation across color channels.
sigma_est = estimate_sigma(noisy, channel_axis=-1, average_sigmas=True)
# Due to clipping in random_noise, the estimate will be a bit smaller than the
# specified sigma.
print(f'Estimated Gaussian noise standard deviation = {sigma_est}')

im_bayes = denoise_wavelet(
    noisy,
    channel_axis=-1,
    convert2ycbcr=True,
    method='BayesShrink',
    mode='soft',
    rescale_sigma=True,
)
im_visushrink = denoise_wavelet(
    noisy,
    channel_axis=-1,
    convert2ycbcr=True,
    method='VisuShrink',
    mode='soft',
    sigma=sigma_est,
    rescale_sigma=True,
)

# VisuShrink is designed to eliminate noise with high probability, but this
# results in a visually over-smooth appearance.  Repeat, specifying a reduction
# in the threshold by factors of 2 and 4.
im_visushrink2 = denoise_wavelet(
    noisy,
    channel_axis=-1,
    convert2ycbcr=True,
    method='VisuShrink',
    mode='soft',
    sigma=sigma_est / 2,
    rescale_sigma=True,
)
im_visushrink4 = denoise_wavelet(
    noisy,
    channel_axis=-1,
    convert2ycbcr=True,
    method='VisuShrink',
    mode='soft',
    sigma=sigma_est / 4,
    rescale_sigma=True,
)

# Compute PSNR as an indication of image quality
psnr_noisy = peak_signal_noise_ratio(original, noisy)
psnr_bayes = peak_signal_noise_ratio(original, im_bayes)
psnr_visushrink = peak_signal_noise_ratio(original, im_visushrink)
psnr_visushrink2 = peak_signal_noise_ratio(original, im_visushrink2)
psnr_visushrink4 = peak_signal_noise_ratio(original, im_visushrink4)

ax[0, 0].imshow(noisy)
ax[0, 0].axis('off')
ax[0, 0].set_title(f'Noisy\nPSNR={psnr_noisy:0.4g}')
ax[0, 1].imshow(im_bayes)
ax[0, 1].axis('off')
ax[0, 1].set_title(f'Wavelet denoising\n(BayesShrink)\nPSNR={psnr_bayes:0.4g}')
ax[0, 2].imshow(im_visushrink)
ax[0, 2].axis('off')
ax[0, 2].set_title(
    'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}$)\n'
    'PSNR=%0.4g' % psnr_visushrink
)
ax[1, 0].imshow(original)
ax[1, 0].axis('off')
ax[1, 0].set_title('Original')
ax[1, 1].imshow(im_visushrink2)
ax[1, 1].axis('off')
ax[1, 1].set_title(
    'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}/2$)\n'
    'PSNR=%0.4g' % psnr_visushrink2
)
ax[1, 2].imshow(im_visushrink4)
ax[1, 2].axis('off')
ax[1, 2].set_title(
    'Wavelet denoising\n(VisuShrink, $\\sigma=\\sigma_{est}/4$)\n'
    'PSNR=%0.4g' % psnr_visushrink4
)
fig.tight_layout()

plt.show()