数字图像处理——第五章 图像复原与重建

---

前言

对比图像增强来学习图像复原技术的概念，它们两个都是以预先确定的目标来改善图像，但图像增强主要是一个主观过程，而图像复原则大部分是一个客观过程。图像复原试图利用退化现象的某种先验知识来复原被退化的图像。因而，复原技术是面向退化模型的，并且采用相反的过程进行处理，以便恢复出原图像。

5.1 图像退化复原过程的模型

一个线性模型如下图所示：

退化过程被建模为一个退化函数和一个加性噪声项，如果H是一个线性的、位置不变的过程，那么空间域的退化图像可表示为：$$g(x,y)=h(x,y)\star f(x,y)+\eta (x,y)$$
由于空间域中的卷积等于频率域中的乘积，在频率域中等价表示为：$$G(u,v)=H(u,v)F(u,v)+N(u,v)$$

5.2 噪声模型

噪声主要来源于图像获取或者传输过程，如图像获取过程中的环境条件和传感元器件自身的质量，图像传输中传输信道中的干扰。

5.2.1 噪声的空间和频域特性

噪声的频率特性是指傅里叶域中噪声的频率内容，例如，当噪声的傅里叶谱是常量时，噪声通常称为白嗓声。
噪声的空间特性我们一般不讨论（除了空间周期噪声外）。

5.2.2 噪声概率密度函数

一些噪声的概率密度函数图像如下：

1.高斯噪声：$$p(z)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-(z-\overline{z}{)}^2/2{\sigma }^2}$$
$\overline{z}$表示z的平均值，$\sigma$表示z的标准差。

添加高斯噪声代码示例：

import numpy as np
import matplotlib.pyplot as plt
import cv2 as cv
def show(img):
    if img.ndim == 2:
        plt.imshow(img, cmap='gray')
    else:
        img = cv.cvtColor(img, cv.COLOR_BGR2RGB)
        plt.imshow(img)
    plt.show()
img = cv.imread('pic/sea.jpg')
GaussNoise1 = np.random.normal(0, 0.1, img.shape)
GaussNoise2 = np.random.normal(0, 1, img.shape)
img_GaussNoise1 = img + GaussNoise1
img_GaussNoise2 = img + GaussNoise2
img_GaussNoise1 = np.uint8(cv.normalize(img_GaussNoise1, None, 0, 255, cv.NORM_MINMAX))  # 归一化为 [0,255]
img_GaussNoise2 = np.uint8(cv.normalize(img_GaussNoise2, None, 0, 255, cv.NORM_MINMAX))  # 归一化为 [0,255]
show(np.hstack([img,img_GaussNoise1,img_GaussNoise2]))

添加高斯噪声图片示展示：

由图片可以看到，添加高斯噪声时，方差/标准差越大，噪声越强。

2.瑞利噪声：$$p(z)=\{\begin{array}{l}\frac{2}{b}(z-a){\mathrm{e}}^{-(z-a{)}^2/b},z⩾a\\ 0,z<a\end{array}$$
其中，$\begin{array}{c}\overline{z}=a+\sqrt{\pi b/4}，{\sigma }^2=\frac{b(4-\pi )}{4}\end{array}$。
3.爱尔兰（伽马）噪声：$$\begin{array}{rl}p\left(z\right)& =\{\begin{array}{ll}\frac{a^b{z}^{b-1}}{(b-1)!}{\mathbf{e}}^{-az},& z\ge a\\ 0,& z<a\end{array}\end{array}$$
其中，$\begin{array}{c}\overline{z}=\frac{b}{a}，{\sigma }^2=\frac{b}{a^2}\end{array}$。
4.指数噪声：$$p(z)=\{\begin{array}{ll}a{e}^{-az},& z⩾0\\ 0,& z<0\end{array}$$
其中，$\begin{array}{rlrl}& \overline{z}=\frac{1}{a}，& & {\sigma }^2=\frac{1}{a^2}\end{array}$。
5.均匀噪声：$$\begin{array}{r}b\\ p(z)=\{\begin{array}{ll}\frac{1}{b-a},& a\le z\le b,\\ 0，& 其他\end{array}\end{array}$$
其中，$\begin{array}{c}\overline{z}=\frac{a+b}{2}，{\sigma }^2=\frac{{\left(b-a\right)}^2}{12}\end{array}$。
6.脉冲（椒盐）噪声：$$\begin{array}{r}p\left(z\right)=\{\begin{array}{llc}{P}_a,& & z=a\\ P_b,& & z=b\\ 1-P_a-P_b,& & \mathfrak{其他}\end{array}\end{array}$$

添加椒盐噪声代码示例：

def salt_pepper_noise(img, ratio=0.03):
    img = img.copy()
    for _ in range(int(img.size * ratio)):
        i = np.random.randint(0, img.shape[0])
        j = np.random.randint(0, img.shape[1])     
        img[i,j] = np.random.randint(0, 100)
    print(img.min())
    return img

img = cv.imread('pic/sea.jpg', 0)
img_sp = salt_pepper_noise(img)
show(np.hstack([img, img_sp]))

添加椒盐噪声图片示展示：
椒盐噪声就像在图片上撒上了像椒盐一样的黑点或者白点。

5.2.3 周期噪声

周期噪声是本章在空间域上唯一讨论的噪声，它一般是由于在图像获取期间由电力或机电干扰产生的，可通过频率域滤波来显著减少。

5.2.4 噪声参数的估计

噪声参数可以通过检测图像的傅立叶谱来估计，当仅有通过传感器生成的图像可用时，通常可由合理的恒定灰度值的一小部分来估计 PDF参数。

5.3 只存在噪声的复原——空间滤波

在空间域和频率域上用公式可表示为：$$\begin{array}{rlrl}& g\left(x,y\right)& & =f(x,y)+\eta (x,y)\\ & G(u,v)& & =F(u,v)+N(u,v)\\ & & & \end{array}$$

5.3.1 均值滤波器

1.算术均值滤波器：$$\hat{f}(x,y)=\frac{1}{mn}\sum _{(s,t)\in S_n}g(s,t)$$
它模糊了结果，也降低了噪声。
2.几何均值滤波器：$$\hat{f}(x,y)=[\prod _{(s,t)\in S_m}g(s,t){]}^{\frac{1}{min}}$$
相比于算术均值滤波器，它丢失的图像细节更少。
3.谐波均值滤波器：$$\hat{f}(x,y)=\frac{mn}{{\sum }_{(x,t)\le S_g}\frac{1}{g(s,t)}}$$
谐波均值滤波器对于盐粒噪声效果较好，但不适用于胡椒噪声，它善于处理像高斯噪声那样的其他噪声。
4.逆谐波均值滤波器：$$\hat{f}(x,y)=\frac{\sum _{(s,t)\le S_n}g(s,t{)}^{2+1}}{\sum _{(s,t)\le S_n}g(s,t{)}^Q}$$
其中Q为滤波器的阶数，这种滤披器适合减少或在实际中消除椒盐噪声的影响，当Q值为正时，该滤波器消除胡椒噪声；当Q值为负时，该滤波器消除盐粒噪声，但它不能同时消除这两种噪声。注意，当 Q=O 时，逆谐波均值滤波器简化为算术均值滤波器；而当Q= -1 时，则为谐波均值滤破器

总之，算术均值滤波器和几何均值滤波器（尤其是后者）更适合于处理高斯或均匀随机噪声，逆谐波均值滤波器更适合处理脉冲噪声，但它还有一个缺点，即必哪知道噪声是暗噪声还是亮噪声，以便于为Q选择的正确的符号。

5.3.2 统计排序滤波器

1.中值滤波器： f ^ ( x , y ) = median ( s , t ) ∈ S n { g ( s , t ) } \hat{f}(x,y)=\underset{(s,t)\in S_n}{\text{median}}\left\{g(s,t)\right\} f^​(x,y)=(s,t)∈Sn​median​{g(s,t)}
对于某些类型的随机噪声，中值滤波器可提供良好的去噪能力，且比相同尺寸的线性平滑滤波器引起的模糊更少，在存在单极或双极脉冲噪声的情况下，中值滤波器尤其有效。
2.最大值和最小值滤波器：
f ^ ( x , y ) = max ⁡ ( s , t ) ∈ S n { g ( s , t ) } \hat{f}(x,y)=\max\limits_{(s,t)\in S_n}\left\{g(s,t)\right\} f^​(x,y)=(s,t)∈Sn​max​{g(s,t)}
最大值滤波器对于发现图像中最亮点非常有用，可以降低胡椒噪声。
f ^ ( x , y ) = min ⁡ ( s , t ) ∈ S n { g ( s , t ) } \hat{f}(x,y)=\min\limits_{(s,t)\in S_n}\left\{g(s,t)\right\} f^​(x,y)=(s,t)∈Sn​min​{g(s,t)}
最小值滤波器对于发现图像中最暗点非常有用，可以降低胡椒噪声。
3.中点滤波器： f ^ ( x , y ) = 1 2 [ max ⁡ ( s , t ) ∈ S n { g ( s , t ) } + min ⁡ ( s , t ) ≤ S n { g ( s , t ) } ] \hat{f}(x,y)=\frac{1}{2}\left[\max\limits_{(s,t)\in S_n}\left\{g(s,t)\right\}+\min\limits_{(s,t)\le S_n}\left\{g(s,t)\right\}\right] f^​(x,y)=21​[(s,t)∈Sn​max​{g(s,t)}+(s,t)≤Sn​min​{g(s,t)}]
中点滤波器对于随机分布噪声效果最好，如高斯噪声或均匀噪声。
4.修正的的阿尔法均值滤波器：$$\hat{f}(x,y)=\frac{1}{mn-d}\sum _{(s,t)\in S_r}{g}_r(s,t)$$
其中，d的取值范围可为0打破到mn-1，当d=0时，退化为均值滤波器；当d=mn-1时，退化为中值滤波器；当d取其他值时，在包括多种噪声的情况下很有用，例如高斯噪声和椒盐噪声混合的情况下。

统计排序滤波器代码：

 img = cv.imread("pic/rose_gaussiannoise_200x200.jpg", 0)  # flags=0 读取为灰度图像
    show(img)
    img_h = img.shape[0]
    img_w = img.shape[1]

    m, n = 3, 3
    kernalMean = np.ones((m, n), np.float32)  # 生成盒式核

    # 边缘填充
    hPad = int((m-1) / 2)
    wPad = int((n-1) / 2)
    imgPad = np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), mode="edge")

    imgMedianFilter = np.zeros(img.shape)  # 中值滤波器
    imgMaxFilter = np.zeros(img.shape)  # 最大值滤波器
    imgMinFilter = np.zeros(img.shape)  # 最小值滤波器
    imgMiddleFilter = np.zeros(img.shape)  # 中点滤波器
    for i in range(img_h):
        for j in range(img_w):
            # # 1. 中值滤波器 (median filter)
            pad = imgPad[i:i+m, j:j+n]
            imgMedianFilter[i, j] = np.median(pad)

            # # 2. 最大值滤波器 (maximum filter)
            pad = imgPad[i:i+m, j:j+n]
            imgMaxFilter[i, j] = np.max(pad)

            # # 3. 最小值滤波器 (minimum filter)
            pad = imgPad[i:i+m, j:j+n]
            imgMinFilter[i, j] = np.min(pad)

            # # 4. 中点滤波器 (middle filter)
            pad = imgPad[i:i+m, j:j+n]
            imgMiddleFilter[i, j] = int(pad.max()/2 + pad.min()/2)

    plt.figure(figsize=(9, 9))
    plt.subplot(221), plt.axis('off'), plt.title("median filter")
    plt.imshow(imgMedianFilter, cmap='gray', vmin=0, vmax=255)
    plt.subplot(222), plt.axis('off'), plt.title("maximum filter")
    plt.imshow(imgMaxFilter, cmap='gray', vmin=0, vmax=255)
    plt.subplot(223), plt.axis('off'), plt.title("minimum filter")
    plt.imshow(imgMinFilter, cmap='gray', vmin=0, vmax=255)
    plt.subplot(224), plt.axis('off'), plt.title("middle filter")
    plt.imshow(imgMiddleFilter, cmap='gray', vmin=0, vmax=255)
    plt.tight_layout()
    plt.show()

加了高斯噪声的图像：

不同滤波之后的效果展示：
我们可以看到中值滤波器对于某些随机噪声具有较好的降噪能力，最大值滤波器可有效降低图像中胡椒噪声（黑点），最小值滤波器可有效降低图像中盐粒噪声（白点），中点滤波器对于这种高斯噪声具有较好降噪能力。

5.3.3 自适应滤波器

自适应滤波器能适应图像具体特征来进行图像复原，所以它的性能要比之前讨论过的滤波器性能都要好。
1.自适应局部降低噪声滤波器：
滤波器作用于局部区域Sxy。滤波器在该区域中心任意一点(x,y)上的响应基于以下4个量：
(1)g(x, y)，噪声图像在点(x, y)上的值
(2)${\mathbf{\sigma }}_{\eta }^2$,干扰f(x, y)以形成g(x, y)的噪声方差
(3)$m_L$在Sxy上像素点的局部均值
(4)${\mathbf{\sigma }}_L^2$,在Sxy上像素点的局部方差

我们希望滤波器性能如下：
(1)如果${\mathbf{\sigma }}_{\eta }^2$​为零，则滤波器应该简单地返回g(x, y)的值。这无关紧要，在零噪声情况下g(x, y)等于f(x,y)。
2.如果局部方差与${\mathbf{\sigma }}_{\eta }^2$​是高度相关的，则滤波器返回g(x, y)的一个近似值。典型地，高局部方差与边缘相关， 并且应该保护这些边缘。
3.如果两个方差相等，我们]则希望滤波器返回Sxy中像素的算术均值。这种情况发生在局部区域与整个图像 有相同特性的条件下，并且局部噪声将通过简单地求平均来降低。
基于这些假设得到的$\hat{f}(x,y)$自适应表达式可以写成:$$\hat{f}(x,y)=g(x,y)-\frac{{\sigma }_{\eta }^2}{{\sigma }_L^2}[g(x,y)-m_L]$$

自适应局部降低噪声滤波器代码：

img = cv.imread("pic/rose_gaussiannoise_200x200.jpg", 0)  
hImg = img.shape[0]
wImg = img.shape[1]

m, n = 5, 5
imgAriMean = cv.boxFilter(img, -1, (m, n))  # 算术平均滤波

    # 边缘填充
    hPad = int((m-1) / 2)
    wPad = int((n-1) / 2)
    imgPad = np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), mode="edge")

    # 估计原始图像的噪声方差 sigmaEta
    mean, stddev = cv.meanStdDev(img)
    sigmaEta = stddev ** 2
    print(sigmaEta)

    # 自适应局部降噪
    epsilon = 1e-8
    imgAdaLocal = np.zeros(img.shape)
    for i in range(hImg):
        for j in range(wImg):
            pad = imgPad[i:i+m, j:j+n]  # 邻域 Sxy, m*n
            gxy = img[i,j]  # 含噪声图像的像素点
            zSxy = np.mean(pad)  # 局部平均灰度
            sigmaSxy = np.var(pad)  # 灰度的局部方差
            rateSigma = min(sigmaEta / (sigmaSxy + epsilon), 1.0)  # 加性噪声假设：sigmaEta/sigmaSxy < 1
            imgAdaLocal[i, j] = gxy - rateSigma * (gxy - zSxy)

   
    show(np.hstack([img,imgAriMean,imgAdaLocal]))

图片展示：

将第二与第三张图片与第一张对比可以看出来，自适应均值滤波器要比算术平均滤波器效果好的多，滤除了噪声的同时比较好的保留了原图中的细节部分。代码参考来自添加链接描述。

2.自适应中值滤波器：
自适应中值滤波器相比于中值滤波器可以处理更大概率的脉冲噪声，且平滑非脉冲噪声时保留细节。自适应中值滤波器在滤波的过程中，能动态的改变滤波器的窗口尺寸大小，根据预先设定好的条件，改变滤波窗口的尺寸大小，同时还会根据一定的条件判断当前像素是不是噪声，如果是则用邻域中值替换掉当前像素；不是，则不作改变。
自适应中值滤波器算法以两个进程工作，
进程A:
$A_1=Z_{med}-Z_{min}$
$A_2=Z_{med}-Z_{max}$
如果$A_1$>0且$A_2$<0，则转至进程B，否则增大窗口尺寸
如果窗口尺寸$<=S_{max}^{\prime }$，则重复进程A，否则输出$Z_{med}$
进程B：
$B_1=Z_{xy}-Z_{min}$
$B_2=Z_{xy}-Z_{max}$
如果$B_1>0$且$B_2<0$，则输出$Z_{xy}$否则输出$Z_{med}$
它有3个主要目的：去除椒盐(脉冲)噪声，平滑其他非脉冲噪声，并减少诸如物体边界细化或粗化等失真。

5.4 用频率波消除周期噪声

在傅里叶变换中，周期噪声在对应于周期干扰的频率，以集中的能量脉冲形式出现，用频率技术可以有效分析并滤除周期噪声。

5.4.1 带阻滤波器

带阻滤波器在频率域噪声成分的一般位置近似已知的应用中消除噪声。

5.4.2 带通滤波器

带通滤波器执行与带阻滤波器相反的操作，但是通常不会在一幅图像上直接执行带通滤波器，一般用它来屏蔽噪声模式。

5.4.3 陷波滤波器

陷波滤波器阻止（或通过）事先定义中心邻域内的频域，它必须以关于原点对称的形式出现

5.4.4 最佳陷波滤波

当存在几种干扰分量时，采用最佳陷波滤波，能更好的保留图像细节部分，更好复原图像。

5.5 线性、位置不变的退化

退化复原模型中，输入输出关系可表示为：$$\begin{array}{rlr}& & \\ & g\left(x,y\right)& =H\left[f(x,y)\right]+\eta (x,y)\end{array}$$
假设$\eta (x,y)=0$，则$g(x,y)=H[f(x,y)]，则$
1.线性系统：$$H\left[a{f}_1(x,y)+b{f}_2(x,y)\right]=aH\left[f_1(x,y)\right]+bH\left[f_2(x,y)\right]$$
2.加性：$$H\left[f_1(x,y)+f_2(x,y)\right]=H\left[f_1(x,y)\right]+H\left[f_2(x,y)\right]$$
3.均匀性：$$H\left[a{f}_1(x,y)\right]=aH\left[f_1(x,y)\right]$$
4.位置（或空间）不变系统：$$H\left[\bar{f}(x-\alpha ,y-\beta )\right]=g(x-\alpha ,y-\beta )$$
具有加性噪声的线性空间不变退化系统，可在空间域建模为退化（点扩散）函数与一幅图像的卷积，然后再加上噪声。基于卷和定理，在频率域中，同样的过程可表示为图像和退化函数的变换的乘积，然后再加上噪声的变换，可表示为以下公式：
空间域：$$g(x,y)=h(x,y)\star f(x,y)+\eta (x,y)$$
频率域：$$G(u,v)=H(u,v)F(u,v)+N(u,v)$$

5.6 估计退化函数

5.6.1 图像观察估计

从图像本身来收集信息，我们可以观察图像的一个区域，对这一子图像的进行复原处理，根据图像位置不变的假设来估计退化函数。

5.6.2 试验估计

如果可以使用与获取退化图像的设备相似的设备，这样我们可以得到和退化图像相似的图像，从理论上是可以推断出退化函数。

5.6.3 建模估计

建模是从一个基本原理开始来推导出一个数学模型，通过该数学模型推导出退化函数。

5.7 逆滤波

通过以上提到的方法求出退化函数后，由退化函数做一次逆滤波，但是实际上即使知道退化函数，也不能准确地复原未退化的图像（受噪声影响），一般直接逆滤波的性能也是较差的。

5.8 最小均方差（维纳）滤波

最小均方误差（维纳）滤波综合了退化函数和噪声统计特征进行复原处理的方法，对噪声进行了处理，该方法建立在图像和噪声都是随机变量的基础上，目标是找到为污染图像的一个估计，使它们之间均方误差最小，它与逆滤波相比显然效果更好。

5.9 约束最小二乘方滤波

约束最小二乘方滤波不必知道未退化图像和噪声的功率谱，仅要求知道噪声的方差和均值，该方法的核心使H对噪声的敏感性问题，需要使用参数来约束复原，对于高和中噪声的情况更好。

5.10 几何均值滤波

对维纳滤波稍加推广，这种推广就是所谓的几何均值滤波器的形式：$$\begin{array}{c}{\hat{F}}_{\left(u,v\right)}={\left[\frac{H^{\ast }(u,v)}{{\left|H(u,v)\right|}^2}\right]}^{\alpha }{\left[\frac{H^{\ast }(u,v)}{{\left|H(u,v)\right|}^2+\beta \left[\frac{S_{\eta }(u,v)}{S_f(u,v)}\right]}\right]}^{1-\alpha }G(u,v)\end{array}$$
当$\alpha =1$时，该滤波器退化为逆滤波器；
当$\alpha =0$时，该滤波器变为参数维纳滤波器，参数维纳滤波器在$\beta =1/2$时，还原为标准的维纳滤波器；
当$\alpha =1/2$时，该滤波器变为相同次幂的两个量的积，即为几何均值滤波器；

5.11 由投影重建图像

和我们以前在数学中接触到的类似，通过对一个图像或者物体的不同面或者角度的投影来得出该完整的物体或者图像。

5.11.1 计算机断层（CT）原理

X射线计算机断层的目的是使用X射线从许多不同的方向穿过物体而得到该物体内部结构的三维描述。

5.11.2 投影和雷登变换

1.在笛卡尔坐标系中，一条直线可以用斜截式来表示 ：y = a x + b，或者由其法线来描述：$x\cos \theta +y\sin \theta =\rho$，平行射线束的投影可由一组直线建模。
2.雷登变换：$$g(\rho ,\theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\dot{\delta }(x\cos \theta +y\sin \theta -\rho )\text{d}x\text{d}y$$
雷登变换给出了沿着xy平面中任意一条线f(x,y)的投影（线积分）的公式。

5.11.3 傅里叶切片定理

由公式：$$G(\omega ,\theta )={\left[F(u,v)\right]}_{u=\omega \cos \theta ;v=\omega \sin \theta }=F(\omega \cos \theta ,\omega \sin \theta )$$
傅里叶切片定理说明一个投影的傅里叶变换是得到投影区域的二维傅里叶变换的一个切片。

5.11.4 滤波反投影重建图像

直接得到反投影会产生不可接受的模糊效果，我们可以采用平行射线束或者扇形射线束来进行反投影的重建，在反投影之前进行简单的滤波，能够降低模糊的效果。