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

一、图像退化/复原过程的模型

    成像过程中的”退化”，是指由于成像系统各种因素的影响，使得图像质量降低。
    上图中f(x,y)表示一幅输入图像，g(x,y)是f(x,y)产生的一幅退化图像，H表示退化函数，$\eta (x,y)$表示外加噪声。
    如果系统H是一个线性、位置不变性的过 程，退化图像可以表示为：
    $g(x,y)=h(x,y)★h(x,y)+\eta (x,y)$
    $\Rightarrow G(u,v)=H(u,v)F(u,v)+N(u,v)$
（空间域上的卷积等同于频率域上的乘积）其中h(x,y)表示退化函数的空间表示。

二、噪声模型

    数字图像的噪声主要来源于图像的获取和传输过程。图像获取的数字化过程，如图像传感器的 质量和环境条件；图像传输过程中传输信道的噪声干扰，如 通过无线网络传输的图像会受到光或其它 大气因素的干扰。

2.1一些重要的噪声概率密度函数

2.1.1高斯噪声

    高斯噪声的概率密度函数(PDF)由下式给出：
$$p(z)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-(z-\overline{z}{)}^2/2{\sigma }^2}$$
z的值有70%落在$[\mu -\sigma ,\mu +\sigma ]$范围内，有95%落在$[\mu -2\sigma ,\mu +2\sigma ]$范围内。

2.1.2瑞利噪声

    瑞利噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{cc}\frac{2}{b}(z-a)e^{-(z-a{)}^2/b}& z\ge a\\ 0& z<a\end{array}$$
$$\mu =a+\sqrt{\pi b/4},{\sigma }^2=\frac{b(4-\pi )}{4}$$
    距离原点的位移是a，函数曲线向右变形。

2.1.3爱尔兰（伽马）噪声

    伽马噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{cc}\frac{a^b{z}^{b-1}}{(b-1)!}{e}^{-az}& z\ge a\\ 0& z<a\end{array}$$
$$\mu =\frac{b}{a},{\sigma }^2=\frac{b}{a^2}$$

2.1.4指数噪声

    指数噪声的PDF由下式给出（a>0）：
$$p(z)=\{\begin{array}{cc}a{e}^{-az}& z\ge 0\\ 0& z<0\end{array}$$
$$\mu =\frac{1}{a},{\sigma }^2=\frac{1}{a^2}$$
    指数分布的PDF是当b=1时爱尔兰分布的特殊情况。

2.1.5均匀噪声

    均匀分布噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{cc}\frac{1}{b-a}& a\le z\le b\\ 0& 其他\end{array}$$
$$\mu =\frac{a+b}{2},{\sigma }^2=\frac{(b-a{)}^2}{12}$$

2.1.6脉冲（椒盐）噪声

    脉冲噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{cc}{P}_a& z=a\\ P_b& z=b\\ 1-P_a-P_b& 其他\end{array}$$
    如果$P_a$或$P_b$为零，则脉冲噪声称为单极脉冲；如果$P_a$或$P_b$均不为零，则脉冲噪声称为双极脉冲噪声或椒盐噪声。
    脉冲噪声可以为正，也可为负。标定以后，脉冲噪声总是数字化为最大值（纯黑或纯白）。通常，负脉冲以黑点（胡椒点）出现，正脉冲以白点（盐点）出现。

2.1.7周期噪声

    周期噪声是在图像获取中从电力或机电干扰中产生，周期噪声可以通过频率域滤波显著减少。

2.2几种噪声的概率密度函数曲线及应用

2.2.1函数曲线

    一些重要的概率密度函数曲线如下图所示：

    下图为用于噪声模型的测试图，由简单、恒定的区域组成，仅仅有3个灰度级的变化。

    对样本噪声图像添加高斯、瑞利、伽马、指数、均匀和椒盐噪声后的图像与直方图如下图所示。

    结论：
    1、上述噪声图像的直方图和它们的概率密度函数曲线对应相似。
    2、添加各种噪声后的图像并没有显著不同，但它们的直方图具有明显的区别。

2.2.2应用

    1、高斯噪声源于电子电路噪声和由低照明度或高温带来的传感器噪声。
    2、瑞利噪声对分布在图像范围内特征化噪声有用。
    3、伽马分布和指数分布用于激光成像噪声。
    4、均匀密度分布作为模拟随机数产生器的基础。
    5、脉冲噪声用于成像中的短暂停留中，如错误的开关操作。

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

    当唯一退化是噪声时：
$$g(x,y)=f(x,y)+\eta (x,y)$$
$$\Rightarrow G(u,v)=F(u,v)+N(u,v)$$
    噪声项未知，不能从g(x,y)或G(u,v)减去噪声，可以选择空间滤波方法进行图像复原。

3.1均值滤波器

3.1.1算术均值滤波器

$$\hat{f}(x,y)=\frac{1}{mn}\sum _{(s,t)\in S_{xy}}g(s,t)$$
    $S_{xy}$表示中心在(x,y)，尺寸为m×n的矩形窗口，平滑了一幅图像的局部变化，在模糊了结果的同时减少了噪声。

3.1.2几何均值滤波器

$$\hat{f}(x,y)=[\prod _{(s,t)\in S_{xy}}g(s,t){]}^{\frac{1}{mn}}$$
    几何均值滤波器所达到的平滑度可以与算术 均值滤波器相比，但几何均值滤波器在滤波过程中，与算术均值滤波器相比，会丢失更少的图像细节——相对锐化。

3.1.3谐波均值滤波器

$$\hat{f}(x,y)=\frac{mn}{{\sum }_{(s,t)\in S_{xy}}\frac{1}{g(s,t)}}$$
    谐波均值滤波器对于“盐”噪声效果好，但不适用于“胡椒”噪声，善于处理高斯噪声。

3.1.4逆谐波均值滤波器

$$\hat{f}(x,y)=\frac{{\sum }_{(s,t)\in S_{xy}}g(s,t{)}^{Q+1}}{{\sum }_{(s,t)\in S_{xy}}g(s,t{)}^Q}$$
    Q称为滤波器的阶数。当Q为正数时，用于消除“胡 椒”噪声；当Q为负数时，用于消除“盐”噪声，但不能同时消除“椒盐”噪声。
    当Q=0，逆谐波均值滤波器转变为算术均值滤波器；当Q=-1，逆谐波均值滤波器转变为谐波均值滤波器。

    均值滤波器总结：算术均值滤波器和几何均值滤波器适合于处理高斯或均匀等随机噪声。谐波均值滤波器适合于处理脉冲噪声，缺点是必须事先知道噪声是暗噪声还是亮噪声，以便于选择合适的Q符号。

3.2统计排序滤波器

3.2.1中值滤波器

$$\hat{f}(x,y)=media{n}_{(s,t)\in S_{xy}}g(s,t)$$
    在相同尺寸下，比起均值滤波器引起的模糊少，对单极或双极脉冲噪声非常有效。

3.2.2最大值和最小值滤波器

最大值滤波器：
$$\hat{f}(x,y)=ma{x}_{(s,t)\in S_{xy}}g(s,t)$$
    用于发现图像中的最亮点，可以有效过滤“胡椒”噪声（因为“胡椒”噪声 是非常低的值）。
最小值滤波器：
$$\hat{f}(x,y)=mi{n}_{(s,t)\in S_{xy}}g(s,t)$$
    用于发现图像中的最暗点，可以有效过滤“盐”噪声（因为“盐”噪声是非 常高的值）。

3.2.3中点滤波器

$$\hat{f}(x,y)=\frac{1}{2}[ma{x}_{(s,t)\in S_{xy}}g(s,t)+mi{n}_{(s,t)\in S_{xy}}g(s,t)]$$
    结合了顺序统计和求平均，对于高斯和均匀随机分布这类噪声有最好的效果。

3.2.4修正的阿尔法均值滤波器

$$\hat{f}(x,y)=\frac{1}{mn-d}\sum _{(s,t)\in S_{(}xy)}{g}_r(s,t)$$
    在$S_{xy}$邻域内去掉$g(s,t)$最高灰度值的d/2和最低灰度值的d/2。$g_r(s,t)$代表剩余的mn-d个像素，当d=0，退变为算术均值滤波器；当d=(mn-1)/2，退变为中值滤波器；当d取其它值时，适用于包括多种噪声的情况下，例如高斯噪声和椒盐噪声混合的情况。

3.3自适应滤波器

3.3.1自适应局部降低噪声滤波器

    滤波器响应基于以下4个量：
1、g(x,y)，表示噪声图像在点（x,y）上的值；
2、${\sigma }_n^2$，g(x,y)的噪声方差；
3、$m_L$，在$S_{xy}$上像素点的局部均值;
4、${\sigma }_L^2$，在$S_{xy}$上像素点的局部方差。
    滤波器的预期性能如下：
1、如果${\sigma }_n=0$，滤波器返回g(x,y)的值。因为在g(x,y)下零噪声的情况等同于f(x,y)。
2、如果局部方差 ${\sigma }_L^2$与${\sigma }_n^2$高相关，滤波器返回一个g(x,y)的近似值。
3、如果${\sigma }_L^2={\sigma }_n^2$，滤波器返回区域$S_{xy}$上像素的算术均值。这样局部噪声用求平均$m_L$来降低。
    基于这些假设得到的$\hat{f}(x,y)$的自适应表达式可以写成：
$$\hat{f}(x,y)=g(x,y)-\frac{{\sigma }_{\eta }^2}{{\sigma }_L^2}[g(x,y)-m_L]$$
    唯一需要知道或估计的未知量是噪声方差${\sigma }_n^2$，其它参数可以从$S_{xy}$中的像素计算出来。

    通过对比上图，从总体噪声减少情况来看，自适应滤波器所达到的效果与算术和几何均值滤波器的相似，然而使用自适应滤波后的图像更清晰，更尖锐。

3.3.2自适应中值滤波器

    传统中值滤波器只能处理空间密度不大的冲激噪声（$p_a$，$p_b$<0.2），而自适应中值滤波器可以处理具有更大概率的冲激噪声，可以在平滑非冲激噪声时保存细节，而传统中值滤波器无法做到。
    考虑如下符号：
    1、$z_{min}=S_{xy}$中灰度级的最小值
    2、$z_{max}=S_{xy}$中灰度级的最大值
    3、$z_{med}=S_{xy}$中灰度级的中值
    4、$z_{xy}$＝在坐标（x,y）上的灰度级
    5、$S_{max}=S_{xy}$允许的最大尺寸
进程A：
    $A_1=z_{med}-z_{min}$
    $A_2=z_{med}-z_{max}$
    如果$A_1＞0$且$A_2＜0$，则转到进程B，否则增大窗口尺寸
    如果窗口尺寸≤$S_{max}$，则重复进程A
    否则输出$z_{med}$
进程B：
    $B_1=z_{xy}-z_{min}$
    $B_1=z_{xy}-z_{max}$
    $B_1＞0$且$B_2＜0$，则输出$z_{xy}$
    否则输出$z_{med}$

    算法主要目的：除去“椒盐”噪声（冲激噪声）；平滑其它非冲激噪声；减少物体边界细化或粗化等失真。

    通过对比上图，可以看出自适应滤波器在保持清晰度和细节方面做得更好。

3.3.3Python实现

# 自适应中值滤波
# count 为最大窗口数，original 为原图
def adaptiveMedianDeNoise(count, original):
    # 初始窗口大小
    startWindow = 3
    # 卷积范围
    c = int(count/2)
    rows, cols = original.shape
    newI = np.zeros(original.shape)
    for i in range(c, rows - c):
        for j in range(c, cols - c):
            k = int(startWindow / 2)
            median = np.median(original[i - k:i + k + 1, j - k:j + k + 1])
            mi = np.min(original[i - k:i + k + 1, j - k:j + k + 1])
            ma = np.max(original[i - k:i + k + 1, j - k:j + k + 1])
            if mi < median < ma:
                if mi < original[i, j] < ma:
                    newI[i, j] = original[i, j]
                else:
                    newI[i, j] = median

            else:
                while True:
                    startWindow = startWindow + 2
                    k = int(startWindow / 2)
                    median = np.median(original[i - k:i + k + 1, j - k:j + k + 1])
                    mi = np.min(original[i - k:i + k + 1, j - k:j + k + 1])
                    ma = np.max(original[i - k:i + k + 1, j - k:j + k + 1])

                    if mi < median < ma or startWindow > count:
                        break

                if mi < median < ma or startWindow > count:
                    if mi < original[i, j] < ma:
                        newI[i, j] = original[i, j]
                    else:
                        newI[i, j] = median

    return newI


def medianDeNoise(original):
    rows, cols = original.shape
    ImageDenoise = np.zeros(original.shape)
    for i in range(3, rows - 3):
        for j in range(3, cols - 3):
            ImageDenoise[i, j] = np.median(original[i - 3:i + 4, j - 3:j + 4])
    return ImageDenoise

def show(f, s, a, b, c):
    plt.subplot(a, b, c)
    plt.imshow(f, "gray")
    plt.axis('on')
    plt.title(s)

def main():
    original = plt.imread("lena.tiff", 0)
    rows, cols = original.shape

    original_noise = pepperNoise(100000, original)
    adapMedianDeNoise = adaptiveMedianDeNoise(7, original_noise)
    mediDeNoise = medianDeNoise(original_noise)
    plt.figure()
    show(original, "original", 2, 2, 1)
    show(original_noise, "original_noise", 2, 2, 2)
    show(adapMedianDeNoise, "adaptiveMedianDeNoise", 2, 2, 3)
    show(mediDeNoise, "medianDeNoise", 2, 2, 4)
    plt.show()

结果：

（参考博客链接：https://blog.csdn.net/ab136681/article/details/104264811/）

四、用频率域滤波消除周期噪声

4.1带阻滤波器

    阻止一定频率范围内的信号通过而允许其它频率范围内的信号通过，消除或衰减傅里叶变换原点处的频段。
    滤波器的透视图如下图所示：

4.2带通滤波器

    允许一定频率范围内的信 号通过而阻止其它频率范围内的信号通过。
$$H_{BP}(u,v)=1-H_{BR}(u,v)$$
    $H_{BP}(u,v)$表示带通滤波器，$H_{BR}(u,v)$表示相应的带阻滤波器。

4.3陷波滤波器

    阻止或通过事先定义的中心频率邻域内的频率。由于傅里叶变换是对称的，陷波滤波器必须以关于原点对称的形式出现。如果陷波滤波器位于原点处，则以它本身形式出现。
陷波滤波器的三维图如下图所示：

传递函数由下式给出：
$$H_{NP}(u,v)=1-H_{NR}(u,v)$$
    式中，$H_{NP}(U,V)$是陷波滤波器的传递函数，这个陷波滤波器与传递函数为$H_{NR}(U,V)$的陷波滤波器相对应。

五、线性、位置不变的退化

$$g(x,y)=H[f(x,y)]+\eta (x+y)$$
假设$\eta (x+y)=0$,则$g(x,y)=H[f(x,y)]$,如果
$$H[a{f}_1(x,y)+f_2(x,y)]=aH[f_1(x,y)+bH[f_2(x,y)$$
则系统H是一个线性系统，其中a，b是标量，$f_1(x,y)$和$f_2(x,y)$是两幅输入图像。
    具有加性噪声的线性空间不变退化系统，可在空间域建模为退化函数与一幅图像的卷积，然后再加上噪声。基于卷积定理，在频率域中，同样的过程可表示为图像和退化函数的变换的乘积，然后再加上噪声的变换。
    这种方法的优点是可以使用许多线性系统理论的工具来解决图像复原问题。与位置有关的非线性技术虽然更普遍，但它们会带来没有已知解的问题，或解决计算问题时非常困难。

六、总结

    图像复原需要针对图像的退化原因进行补偿，利用图像退化的逆过程去恢复原始图像，使复原后的图像尽可能的接近原图像。