第五章 图像复原与重建

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

图像退化

​ 图像退化的描述(图像退化及其过程描述)如下:
​ 图像的退化是指图像在形成、传输和记录过程中，由于成像系统、传输介质和设备的不完善，使图像的质量下降(变坏)。其典型表现为：模糊、失真、有噪声。产生原因：成像系统像差、传感器拍摄姿态和扫描非线性、成像设备与物体运动的相对运动、大气湍流、成像和处理过程中引入的噪声等。

图像复原

​ 图像复原就是尽可能恢复退化图像的本来面目，它是沿图像退化的逆过程进行处理，也就是如果我们知道图像是经历了什么样的过程导致退化，就可以按其逆过程来复原图像。因此，图像复原过程流程如下:
​ 找退化原因→建立退化模型→反向推演→恢复图像
​ 典型的图像复原是根据图像退化的先验知识，建立退化现象的数学模型，再根据模型进行反向的推演运算，以恢复原来的景物图像。因此，图像复原的关键是知道图像退化的过程，即图像退化模型。并据此采用相反的过程求得原始图像。针对不同的退化问题，图像复原的方法主要有：代数方法恢复、运动模糊恢复、逆滤波恢复、维纳滤波恢复、功率谱均衡恢复、约束最小平方恢复、最大后验恢复、最大熵恢复、几何失真恢复等。

图像复原与图像增强

​ 图像增强不考虑图像是如何退化的，而是主观上试图采用各种技术来增强图像的视觉效果。

​ 因此，图像增强可以不顾增强后的图像是否失真，只要达到想要的目视效果就可以。而图像复原就完全不同，需知道图像退化的机制和过程等先验知识，客观上找出一种相应的逆处理方法，从而得到复原的图像。如果图像已退化，应先作复原处理，再作增强处理。
​ 二者的目的都是为了改善图像的质量。

图像退化的数学模型

​ 输入图像f(x,y)经过某个退化系统后输出的是一幅退化的图像。 为了讨论方便，把噪声引起的退化即噪声对图像的影响一般作为加性噪声考虑。原始图像f(x,y)经过一个退化算子或退化系统H(x,y)的作用，再和噪声n(x,y)进行叠加，形成退化后的图像g(x,y)。

​ 下图表示退化过程的输入和输出之间的关系，其中H(x,y)概括了退化系统的物理过程，它就是我们要寻找的退化数学模型。

​

​ 这里的H就是退化系统。
​ 通常，我们假设图像经过的退化系统是线性非时变系统，线性非时变系统具有如下四个基本性质:
​ H的四个性质(假设噪声n (x，y) =0) :

​ (1)一致性(齐次性):如果 $f_2(x,y)=0$, 则:
$$H[k_1(f_1(x,y)]=k_{1}H[f_1(x,y)]$$
​ 线性系统对常数与任意输入的乘积的响应等于常数与输入的响应的乘积
​ (2)相加性(叠加性) :如果k1=k2=1 ,则:
$$H[f_1(x,y)+f_2(x,y)=H[f_1(x,y)]+H[f_2(x,y)]$$
​ 两个图像和的退化结果等于它们分别退化结果之和(说明线性系统对两个输入图像之和的响应等于它对两个输入图像响应的和)。
​ (3)线性(齐次叠加性) :如果令k1和k2为常数,$ f_1(x,y)$和$f_2(x,y)$为两幅输入图像,则:
$$H[k_1{f}_1(x,y)+k_2{f}_2(xy)]=k_{1}H[f_1(x,y)]+k_{2}H[f_2(x,y)]$$
​ 两个图像的加权和的退化结果等于它们分别退化结果的加权和。
​ (4)位置(空间)不变性:如果对任意f(x,y)以及a和b,有:
$$H[f(x-a,y-b)]=g(x-a,y-b)$$
​ 原始图像偏移多少,响应的退化图像也偏移多少,说明线性系统在图像任意位置的响应只与在该位置的输入值有关而与位置本身无关

​ 根据这些特点，输入信号与其经过线性非时变系统的输出信号之间的关系，以及傅里叶变换的性质，我们得到如下时域及频域的关系表达式:
$$\begin{gathered} g(x,y)=H[f(x,y)]+n(x,y) \\ G(u、v)=H(u,v)F(u,v)+N(u,v) \\ g(x,y)=f(x,y{)}^{\ast }h(x,y)+n(x,y) \end{gathered}$$

​ 即:在时域上分析时，原始图像经过退化系统后得到的退化图像g(x,y)等于 原始输入图像f(x,y) 与系统冲激响应h(x,y)的卷积再加上噪声信号。.
​ 在频(率)域上分析时，退化图像的傅里叶变换G(u,v)等于原始图像的傅里叶变换F(u,v) 与退化系统的频率响应H(u,v)相乘，再加上噪声信号的傅里叶变换N(u,v)

为什么采用线性位移不变系统模型来描述图像退化过程

​ 1)由于许多种退化都可以用线性位移不变模型来近似，这样线性系统中的许多数学工具如线性代数，能用于求解图像复原问题，从而使运算方法简捷和快速。
​ 2)当退化不太严重时，一般用线性位移不变系统模型来复原图像，在很多应用中有较好的复原结果，且计算大为简化。
​ 3)实际上，尽管非线性和位移可变的情况能更加准确而普遍地反映图像复原问题的本质，但在数学上求解困难。只有在要求很精确的情况下才用位移可变的模型去求解，其求解也常以位移不变的解法为基础加以修改而成。
​ 换句话说，就是利用了:线性系统的计算简洁快速，同时它也是非线性和位移可变系统的基础，非线性位移可变系统的求解可以转化为线性移不变求解。

5.2 噪声模型x

​ 数字图像中，噪声主要来源于图像的获取和/或传输过程。成像传感器的性能受各种因素的影响，如图像获取过程中的环境条件和传感元器件自身的质量。

5.2.1 噪声的空间和频率特性

​ 与我们的讨论相关的是定义噪声空间特性的参数，以及噪声是否与图像相关。频率特性是指傅里叶域中噪声的频率内容(即相对于电磁波谱的频率)。例如，当噪声的傅里叶谱是常量时，噪声通常称为白噪声。这个术语是从白光的物理特性派生出来的，它以相等的比例包含可见光谱中的几乎所有频率。以相同比例包含所有频率的函数的傅里叶谱是一个常量。
​ 除了空间周期噪声(见5.2.3节)之外，在本章中我们假设噪声独立于空间坐标，并且噪声与图像本身不相关(即像素值与噪声分量的值之间不相关)。虽然这些假设至少在某些应用中(例如X射线和核医学成像的有限量子成像就是一个很好的例子)是无效的，但处理空间相关和相关性噪声的复杂性超出了我们讨论的范围。

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

​ 概率密度函数（PDF）

高斯噪声

​ 在空间域和频率域中，由于高斯噪声在数学上的易处理性，故实践中常用这种噪声(也称为正态噪声)模型。事实上，这种易处理性非常方便，以至于高斯模型常常应用于在一定程度上导致最好结果的场合。
​ 高斯随机变量z的PDF由下式给出:
$$p(z)=\frac{1}{\sqrt{2\Pi }\delta }{e}^{-(z-\overline{z}{)}^2/2{\delta }^2}$$
​ 其中，z表示灰度值，$\overline{z}$表示z的均(平均)值，$\delta$表示z的标准差。标准差的平方${\delta }^2$称为z的方差。

实验：采用高斯噪声对图像处理

import random

import numpy as np
import matplotlib.pyplot as plt
import cv2 as cv

img = cv.imread('C:/Users/Yuao/Pictures/learningTest/2.png',0)

mu, sigma = 0.0, 20.0
noiseGause = np.random.normal(mu, sigma, img.shape)
imgGaussNoise = img + noiseGause
imgGaussNoise = np.uint8(cv.normalize(imgGaussNoise, None, 0, 255, cv.NORM_MINMAX))  # 归一化为 [0,255]

plt.figure(figsize=(9, 3))
plt.subplot(131), plt.title("Origin"), plt.axis('off')
plt.imshow(img, 'gray', vmin=0, vmax=255)
plt.subplot(132), plt.title("GaussNoise"), plt.axis('off')
plt.imshow(imgGaussNoise, 'gray')
plt.subplot(133), plt.title("Gray Hist")
histNP, bins = np.histogram(imgGaussNoise.flatten(), bins=255, range=[0, 255], density=True)
plt.bar(bins[:-1], histNP[:])
plt.tight_layout()
plt.show()

实验结果：

​ 可以看出在原图上加入高斯噪声，图片变得模糊不清，且其灰色直方图呈现正态分布。

瑞利噪声

​ 瑞利噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{l}\frac{2}{b}(z-b)e^{-(z-a{)}^2/b},z\ge a\\ 0,z<a\end{array}$$
​ 概率密度的均值和方差由
$$\overline{z}=a+\sqrt{\Pi b/4}$$
​ 和
$${\delta }^2=\frac{b(4-\Pi )}{4}$$
​ 给出。瑞利密度对于近似歪斜的直方图十分适用。

实验: 采用瑞利噪声对图像处理

import random

import numpy as np
import matplotlib.pyplot as plt
import cv2 as cv

# 瑞利噪声  (RayleighNoise)
img = cv.imread('C:/Users/Yuao/Pictures/learningTest/2.png', 0)  # flags=0 读取为灰度图像
# img = np.ones([256, 256]) * 128

a = 30.0
noiseRayleigh = np.random.rayleigh(a, size=img.shape)
imgRayleighNoise = img + noiseRayleigh
imgRayleighNoise = np.uint8(cv.normalize(imgRayleighNoise, None, 0, 255, cv.NORM_MINMAX))  # 归一化为 [0,255]

plt.figure(figsize=(9, 3))
plt.subplot(131), plt.title("Origin"), plt.axis('off')
plt.imshow(img, 'gray', vmin=0, vmax=255)
plt.subplot(132), plt.title("RayleighNoise"), plt.axis('off')
plt.imshow(imgRayleighNoise, 'gray')
plt.subplot(133), plt.title("Gray Hist")
histNP, bins = np.histogram(imgRayleighNoise.flatten(), bins=255, range=[0, 255], density=True)
plt.bar(bins[:-1], histNP[:])
plt.tight_layout()
plt.show()

实验结果：

​ 可以看出加入瑞利噪声后，原图像出现模糊，且其直方图距原点的位移和密度的基本形状向右变了。

爱尔兰(伽马)噪声

​ 爱尔兰噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{l}\frac{a^b{z}^{b-1}}{(b-1)!}{e}^{-az},z\ge a\\ 0,z<a\end{array}$$
​ 其中，参数a>0，b为正整数，并且"!"表示阶乘。其概率密度的均值和方差由
$$\overline{z}=\frac{b}{a}$$
和
$${\zeta }^2=\frac{b}{a^2}$$

实验：伽马噪声对图像的影响

import numpy as np
import matplotlib.pyplot as plt
import cv2
# # 9.3：伽马噪声 (Gamma Noise)
img = cv2.imread('C:/Users/Yuao/Pictures/learningTest/2.png', 0)  # flags=0 读取为灰度图像
# img = np.ones([256, 256]) * 128

a, b = 10.0, 2.5
noiseGamma = np.random.gamma(shape=b, scale=a, size=img.shape)
imgGammaNoise = img + noiseGamma
imgGammaNoise = np.uint8(cv2.normalize(imgGammaNoise, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]

plt.figure(figsize=(9, 3))
plt.subplot(131), plt.title("Origin"), plt.axis('off')
plt.imshow(img, 'gray', vmin=0, vmax=255)
plt.subplot(132), plt.title("Gamma noise"), plt.axis('off')
plt.imshow(imgGammaNoise, 'gray')
plt.subplot(133), plt.title("Gray hist")
histNP, bins = np.histogram(imgGammaNoise.flatten(), bins=255, range=[0, 255], density=True)
plt.bar(bins[:-1], histNP[:])
plt.tight_layout()
plt.show()

实验结果

指数噪声

​ 指数噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{l}a{e}^{-az},z\ge 0\\ 0,z<0\end{array}$$
其中，a>0。该概率密度函数的均值和方差是
$$\overline{z}=\frac{1}{a}$$
和
$${\delta }^2=\frac{1}{a^2}$$
这个PDF是当b=1时爱尔兰PDF的特殊情况。

实验：指数噪声对图像的影响

import numpy as np
import matplotlib.pyplot as plt
import cv2
# 指数噪声 (Exponential noise)
img = cv2.imread('C:/Users/Yuao/Pictures/learningTest/2.png', 0) # flags=0 读取为灰度图像
# img = np.ones([256, 256]) * 128

a = 10.0
noiseExponent = np.random.exponential(scale=a, size=img.shape)
imgExponentNoise = img + noiseExponent
imgExponentNoise = np.uint8(cv2.normalize(imgExponentNoise, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]

plt.figure(figsize=(9, 3))
plt.subplot(131), plt.title("Origin"), plt.axis('off')
plt.imshow(img, 'gray', vmin=0, vmax=255)
plt.subplot(132), plt.title("Exponential noise"), plt.axis('off')
plt.imshow(imgExponentNoise, 'gray')
plt.subplot(133), plt.title("Gray hist")
histNP, bins = np.histogram(imgExponentNoise.flatten(), bins=255, range=[0, 255], density=True)
plt.bar(bins[:-1], histNP[:])
plt.tight_layout()
plt.show()

实验结果：

均匀噪声

​ 均匀噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{l}\frac{1}{b-a},a\le z\le b\\ 0,其他\end{array}$$
该密度函数的均值由下式给出：
$$\overline{z}=\frac{a+b}{2}$$
它的方差由下式给出：
$${\delta }^2=\frac{(b-a{)}^2}{12}$$

实验：均匀噪声对图像的影响

import numpy as np
import matplotlib.pyplot as plt
import cv2

# 均匀噪声 (Uniform noise)
img = cv2.imread('C:/Users/Yuao/Pictures/learningTest/2.png', 0) # flags=0 读取为灰度图像
# img = np.ones([256, 256]) * 128

mean, sigma = 10, 100
a = 2 * mean - np.sqrt(12 * sigma)  # a = -14.64
b = 2 * mean + np.sqrt(12 * sigma)  # b = 54.64
noiseUniform = np.random.uniform(a, b, img.shape)
imgUniformNoise = img + noiseUniform
imgUniformNoise = np.uint8(cv2.normalize(imgUniformNoise, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]

plt.figure(figsize=(9, 3))
plt.subplot(131), plt.title("Origin"), plt.axis('off')
plt.imshow(img, 'gray', vmin=0, vmax=255)
plt.subplot(132), plt.title("Uniform noise"), plt.axis('off')
plt.imshow(imgUniformNoise, 'gray')
plt.subplot(133), plt.title("Gray hist")
histNP, bins = np.histogram(imgUniformNoise.flatten(), bins=255, range=[0, 255], density=True)
plt.bar(bins[:-1], histNP[:])
plt.tight_layout()
plt.show()

实验结果

脉冲（椒盐）噪声

​ （双极）脉冲噪声的PDF由下式给出：
$$p(z)=\{\begin{array}{l}{P}_a,z=a\\ P_b,z=b\\ 1-P_a-P_b,其他\end{array}$$
​ 如果b>a,则灰度级b在图像中将显示为一个亮点；反之，灰度级a在图像中将显示为一个暗点。若$P_a$或$P_b$为零，则脉冲噪声称为单极脉冲。如果$P_a$和$P_b$两者均不可能为零，尤其是它们近似相等时，则脉冲噪声值将类似于在图像上随机分布的胡椒和盐粉微粒。由于这个原因，双极脉冲噪声也称为椒盐噪声。这种类型的噪声也可以使用散粒噪声和尖峰噪声来称呼。

实验：椒盐噪声对图像的影响

import numpy as np
import matplotlib.pyplot as plt
import cv2

# 椒盐噪声 (Salt-pepper)
img = cv2.imread('C:/Users/Yuao/Pictures/learningTest/2.png', 0) # flags=0 读取为灰度图像
ps, pp = 0.05, 0.02
mask = np.random.choice((0, 0.5, 1), size=img.shape[:2], p=[pp, (1 - ps - pp), ps])
imgChoiceNoise = img.copy()
imgChoiceNoise[mask == 1] = 255
imgChoiceNoise[mask == 0] = 0

plt.figure(figsize=(9, 3))
plt.subplot(131), plt.title("Origin"), plt.axis('off')
plt.imshow(img, 'gray', vmin=0, vmax=255)
plt.subplot(132), plt.title("Choice noise"), plt.axis('off')
plt.imshow(imgChoiceNoise, 'gray')
plt.subplot(133), plt.title("Gray hist")
histNP, bins = np.histogram(imgChoiceNoise.flatten(), bins=255, range=[0, 255], density=True)
plt.bar(bins[:-1], histNP[:])
plt.tight_layout()
plt.show()

实验结果：

5.2.3 周期噪声

​ 一幅图像中的周期噪声是在图像获取期间由电力或机电干扰产生的。这是在本章中讨论的唯一的一种空间相关噪声。正如5.4节所讨论的那样，周期噪声可通过频率域滤波来显著地减少。

5.2.4 噪声参数的估计

​ 当仅有通过传感器生成的图像可用时，通常可由合理的恒定灰度值的一-小部分来估计PDF的参数。

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

​ 当一幅图像中唯一存在的退化是噪声时，图像退化函数可以简化为
$$g(x,y)=f(x,y)+\eta (x,y)$$
和
$$G(u,v)=F(u,v)+N(u,v)$$
​ 噪声项是未知的，故从g(x, y)或G(u, v)中减去它们不是个现实的选择。在周期噪声的情况下,由G(u,v)的谱来估计N(u, v)通常是可能的。然而，这种类型的知识通常只是例外而不是规律。
​ 当仅存在加性噪声的情况下，可以选择空间滤波方法。

5.3.1均值滤波器

算术均值滤波器

​ 这是最简单的均值滤波器。令$S_{xy}$表示中心在点(x,y)处、大小为mxn的矩形子图像窗口(领域)的一组坐标。算术均值滤波器在$S_{xy}$定义的区域中计算被污染图像g(x,y)的平均值。在点（x,y）处复原图像$\hat{f}$的值，就是简单地使用$S_{x}y$定义地区域中地像素计算出的算术均值，即
$$\hat{f}(x,y)=\frac{1}{mn}\sum _{(s,t)\in S_{xy}}g(s,t)$$
​ 这个操作可以使用大小为mxn的一个空间滤波器来实现，其所有的系数均为其值的1/mn。均值滤波平滑- -幅图像中的局部变化，虽然模糊了结果，但降低了噪声。

实验：使用算术均值滤波器对图像复原

import numpy as np
import matplotlib.pyplot as plt
import cv2

# 算术平均滤波器 (Arithmentic mean filter)
img = cv2.imread('C:/Users/Yuao/Pictures/learningTest/5.jpg', 0) # flags=0 读取为灰度图像
kSize = (3, 3)
kernalMean = np.ones(kSize, np.float32) / (kSize[0] * kSize[1])  # 生成归一化盒式核
imgConv1 = cv2.filter2D(img, -1, kernalMean)  # cv2.filter2D 方法

kSize = (5, 5)
imgConv3 = cv2.boxFilter(img, -1, kSize)  # cv2.boxFilter 方法 (默认normalize=True)

plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title("Original")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132), plt.axis('off'), plt.title("filter2D(kSize=[3,3])")
plt.imshow(imgConv1, cmap='gray', vmin=0, vmax=255)
plt.subplot(133), plt.axis('off'), plt.title("boxFilter(kSize=[5,5])")
plt.imshow(imgConv3, cmap='gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()

实验结果

​ 可以看出均值滤波平滑图像中的局部变化，可以降低图像中的噪声，但会模糊图像。

几何均值滤波器

​ 使用几何均值滤波器复原的一幅图像由如下表达式给出：
$$\hat{f}(x,y)={\left[\prod _{(s,t)\in S_{xy}}g(s,t)\right]}^{\frac{1}{mn}}$$
​ 其中，每个复原的像素由子图像窗口中像素的乘积的1/mn 次幂给出。几何均值滤波器实现的平滑可与算术均值滤波器相比，但这种处理中丢失的图像细节更少。

实验：算术均值滤波器与几何均值滤波器对比

import numpy as np
import matplotlib.pyplot as plt
import cv2

# 几何均值滤波器 (Geometric mean filter)
img = cv2.imread('C:/Users/Yuao/Pictures/learningTest/5.jpg', 0) # flags=0 读取为灰度图像
img_h = img.shape[0]
img_w = img.shape[1]

# 算术平均滤波 (Arithmentic mean filter)
kSize = (3, 3)
kernalMean = np.ones(kSize, np.float32) / (kSize[0] * kSize[1])  # 生成归一化盒式核
imgAriMean = cv2.filter2D(img, -1, kernalMean)

# 几何均值滤波器 (Geometric mean filter)
m, n = 3, 3
order = 1 / (m * n)
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")

imgGeoMean = img.copy()
for i in range(hPad, img_h + hPad):
    for j in range(wPad, img_w + wPad):
        prod = np.prod(imgPad[i - hPad:i + hPad + 1, j - wPad:j + wPad + 1] * 1.0)
        imgGeoMean[i - hPad][j - wPad] = np.power(prod, order)

plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title("Original")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132), plt.axis('off'), plt.title("Arithmentic mean filter")
plt.imshow(imgAriMean, cmap='gray', vmin=0, vmax=255)
plt.subplot(133), plt.axis('off'), plt.title("Geometric mean filter")
plt.imshow(imgGeoMean, cmap='gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()

实验结果

​ 可以看出几何均值滤波器实现的平滑与算术平均滤波器相当，但损失的图像细节更少。

谐波均值滤波器

​ 谐波均值滤波器操作由如下表达式给出：
$$\hat{f}(x,y)=\frac{mn}{{\sum }_{(s,t)\in S_{xy}}\frac{1}{g(s,t)}}$$
谐波均值滤波器对于盐粒噪声效果较好，但不适用于胡椒噪声。它善于处理像高斯噪声那样的其他噪声。

实验：谐波均值滤波器对图像的复原

import numpy as np
import matplotlib.pyplot as plt
import cv2

# 谐波平均滤波器 (Harmonic mean filter)
img = cv2.imread('C:/Users/Yuao/Pictures/learningTest/5.jpg', 0) # flags=0 读取为灰度图像
img_h = img.shape[0]
img_w = img.shape[1]

# 算术平均滤波 (Arithmentic mean filter)
kSize = (3, 3)
kernalMean = np.ones(kSize, np.float32) / (kSize[0] * kSize[1])  # 生成归一化盒式核
imgAriMean = cv2.filter2D(img, -1, kernalMean)

# 谐波平均滤波器 (Harmonic mean filter)
m, n = 3, 3
order = m * n
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")

epsilon = 1e-8
imgHarMean = img.copy()
for i in range(hPad, img_h + hPad):
    for j in range(wPad, img_w + wPad):
        sumTemp = np.sum(1.0 / (imgPad[i - hPad:i + hPad + 1, j - wPad:j + wPad + 1] + epsilon))
        imgHarMean[i - hPad][j - wPad] = order / sumTemp

plt.figure(figsize=(9, 6))
plt.subplot(131), plt.axis('off'), plt.title("Original")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132), plt.axis('off'), plt.title("Arithmentic mean filter")
plt.imshow(imgAriMean, cmap='gray', vmin=0, vmax=255)
plt.subplot(133), plt.axis('off'), plt.title("Harmonic mean filter")
plt.imshow(imgHarMean, cmap='gray', vmin=0, vmax=255)
plt.tight_layout()
plt.show()

实验结果：

​ 不难发现谐波平均滤波器既能处理盐粒噪声（白色噪点），又能处理类似于高斯噪声的其他噪声，但不能处理胡椒噪声（黑色噪点）。

逆谐波均值滤波器

​ 逆谐波均值滤波器基于如下表达式产生一幅复原的图像：
$$\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时，则为谐波均值滤波器。

5.3.2统计排序滤波器

​ 统计排序滤波器是空间域滤波器，空间滤波器的响应基于由该滤波器包围的图像区域中的像素值的顺序(排序)。排序结果决定滤波器的响应。

中值滤波器

​ 最著名的统计排序滤波器是中值滤波器，它使用一个像素领域的灰度级中值来替代该像素的值，即
f ^ ( x , y ) = m e d i a n ( s , t ) ∈ S x y { g ( s , t ) } \widehat{f}(x,y)=median_{(s,t)\in S_{xy}}\{g(s,t)\} f ​(x,y)=median(s,t)∈Sxy​​{g(s,t)}
在(x, y)处的像素值是计算的中值。中值滤波器的应用非常普遍，因为对于某些类型的随机噪声，它们可提供良好的去噪能力，且比相同尺寸的线性平滑滤波器引起的模糊更少。在存在单极或双极脉冲噪声的情况下，中值滤波器尤其有效。

最大值和最小值滤波器

​ 尽管中值滤波器是目前为止图像处理中最常用的一种统计排序滤波器，但它绝对不是唯一的一种。中值相当于顺序排序的数组中间的那个数，我们还可以使用序列中的最后一个数值，称为最大值滤波器，由下式给出
f ^ ( x , y ) = m a x ( x , t ) ∈ S x y { g ( s , t ) } \widehat{f}(x,y)=max_{(x,t)\in S_{xy}\{g(s,t)\}} f ​(x,y)=max(x,t)∈Sxy​{g(s,t)}​
这种滤波器对于发现图像中的最亮点非常有用。同样，因为胡椒噪声的值非常低，作为子图像区域$S_{xy}$中这种最大值选择过程的结果，可以用这种滤波器降低它。

​ 这种起始值的滤波器称为最小值滤波器，它由下式给出
f ^ ( x , y ) = m i n ( x , t ) ∈ S x y { g ( s , t ) } \widehat{f}(x,y)=min_{(x,t)\in S_{xy}\{g(s,t)\}} f ​(x,y)=min(x,t)∈Sxy​{g(s,t)}​
这种滤波器对于发现图像中的最暗点非常有用。同样，作为最小值操作的结果，它可以降低盐粒噪声。

中点滤波器

​ 中点滤波器简单地计算滤波器包围区域中最大值和最小值之间的中点，即
f ^ ( x , y ) = 1 2 [ m a x ( s , t ) ∈ S x y { g ( s , t } + m i n ( s , t ) ∈ S x y { g ( s , t } ] \widehat{f}(x,y)=\frac{1}{2}\left[max_{(s,t)\in S_{xy}\{g(s,t\}}+min_{(s,t)\in S_{xy}\{g(s,t\}} \right] f ​(x,y)=21​[max(s,t)∈Sxy​{g(s,t}​+min(s,t)∈Sxy​{g(s,t}​]
注意，这种滤波器结合了统计排序和求平均。它对于随机分布噪声工作得最好，如高斯噪声或均匀噪声。

修正的阿尔法均值滤波器

​ 假设在领域$S_{xy}$内去掉g(s,t)最低灰度值的d/2和最高灰度值的d/2。令$g_r(s,t)$代表剩下的mn-d个像素。由这些剩余像素的平均值形成的滤波器称为修正的阿尔法均值滤波器：
$$\hat{f}(x,y)=\frac{1}{mn-d}\sum _{(s,t)\in S_{xy}}{g}_r(s,t)$$
其中，d的取值范围可为0到mn-1。当d=0时，修正的阿尔法均值滤波器退化为算术均值滤波器。如果选择d=mn- 1,则修正的阿尔法均值滤波器将退化为中值滤波器。当d取其他值时,修正的阿尔法均值滤波器在包括多种噪声的情况下很有用，例如高斯噪声和椒盐噪声混合的情况下。

5.3.3 自适应滤波器

​ 自适应滤波器的性能要优于迄今为止讨论过的所有滤波器的性能。改善滤波能力的代价是滤波器的复杂度提高了。

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

​ 随机变量最简单的统计度量是其均值和方差。作为自适应滤波器的基础这些是合理的参数，因为它们是与图像表观紧密相关的量。均值给出了在其上计算均值的区域中的平均灰度的度量，而方差则给出了该区域的对比度的度量。

​ 滤波器作用于局部区域$S_{xy}$滤波器在该区域中心任意一点(x, y)上的响应基于以下4个量: (a) g(x, y), 带噪图像在点(x, y) 上的值; (b) ${\delta }_{\eta }^2$，污染f(x, y)以形成8(x, y)的噪声的方差; ©$m_L$，$S_{xy}$中像素的局部均值; (d) ${\delta }_L^2$，

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

自适应中值滤波器

​ 中值滤波器，只要脉冲噪声的空间密度不大，性能就会很好。自适应中值滤波器可以处理具有更大概率的脉冲噪声。自适应中值滤波器的另一个优点是平滑非脉冲噪声时试图保留细节，这是传统中值滤波器所做不到的。

​ 考虑如下符号：
$$\begin{gathered} Z_{min}=S_{xy}中的最小灰度值 \\ {Z}_{max}=S_{xy}中的最大灰度值 \\ {Z}_{med}=S_{xy}中的灰度值的中值 \\ {Z}_{xy}=坐标(x,y)处的灰度值 \\ {S}_{max}=S_{xy}允许的最大尺寸‘ \end{gathered}$$
​ 自适应中值滤算法以两个进程工作，表示为进程A和进程B，如下所示：

进程A：
$$\begin{gathered} A_1=Z_{med}-Z_{min} \\ {A}_2=Z_{med}-Z_{max} \\ 如果A_1>0且A_2<0，则转到进程B \\ 否则增大窗口尺寸 \\ 如果窗口尺寸\le S_{max}，则重复进程A \\ 否则输出Z_{med} \end{gathered}$$
进程B:
$$\begin{gathered} B_1=Z_{xy}-Z_{min} \\ {B}_2=Z_{xy}-Z_{max} \\ 如果B_1>0且B_2<0，则输出Z_{xy}否则输出Z_{med} \end{gathered}$$

​ 理解该算法机理的关键在于，要记住它有3个主要目的:去除椒盐(脉冲)噪声，平滑其他非脉冲噪声，并减少诸如物体边界细化或粗化等失真。值$Z_{min}$ 和$Z_{max}$在算法统计上认为是类脉冲噪声分量，即使它们在图像中并不是最低和最高的可能像素值。

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

​ 用频率域技术可以有效地分析并滤除周期噪声。其基本概念是在傅里叶变换中，周期噪声在对应于周期干扰的频率处，以集中的能量脉冲形式出现。其方法是用一个选择性滤波器分离出噪声。为消除基本的周期噪声，将使用三种类型的选择性滤波器(带阻、带通和陷波滤波器)。

5.4.1带阻滤波器

​ 带阻滤波器的主要应用之一是在频率域噪声分量的一般位置近似已知的应用中消除噪声。一个典型的例子就是一幅被加性周期噪声污染的图像，该噪声可被近似为二维正弦函数。

5.4.2带通滤波器

​ 带通滤波器执行与带阻滤波器相反的操作。使用传递函数为$H_{BR}(u,v)$的带阻滤波器得到带通滤波器的传递函数$H_{BP}(u,v)$:
$$H_{BP}(u,v)=1-H_{BR}(u,v)$$

5.4.3陷波滤波

​ 陷波滤波器阻止(或通过)事先定义的中心频率的邻域内的频率。

5.5线性、 位置不变的退化

​ 具有加性噪声的线性空间不变退化系统，可在空间域建模为退化(点扩散)函数与一幅图像的卷积，然后再加上噪声。基于卷积定理，在频率域中，同样的过程可表示为图像和退化函数的变换的乘积，然后再加上噪声的变换。还应记住，在实现离散傅里叶变换时，需要对函数进行填充。
​ 许多类型的退化可近似为线性、位置不变的过程。这种方法的优点是可以使用许多线性系统理论的工具来解决图像复原问题。与位置有关的非线性技术虽然更普遍(通常会更精确)，但它们会带来没有已知解的困难，或者解决计算问题时非常困难。本章将集中在线性.空间不变的复原技术上。由于退化被建模为卷积的结果，并且图像复原试图找到应用相反过程的滤波器，所以术语图像去卷积通常用于表示线性图像复原。同样，用于复原处理的滤波器常常称为去卷积滤波器。

5.6估计退化函数

​ 在图像复原中，主要有3 种用于估计退化函数的方法:(1)观察法，(2)试验法，(3)数学建模法。使用以某种方式估计的退化函数来复原一幅图像的过程，有时称为盲目去卷积，因为真正的退化函数很少能完全知晓。

5.6.1图像观察估计

​ 假设给我们一幅退化图像，而没有关于退化函数H的任何知识。基于图像被线性、位置不变的过程退化的假设，估计H的一种方法就是从图像本身来收集信息。

5.6.2试验估计

​ 如果可以使用与获取退化图像的设备相似的装置，从理论上讲，得到-一个准确的退化估计是可能的。与退化图像类似的图像可以通过各种系统设置得到，直到这些图像退化到尽可能接近我们希望复原的程度。之后，概念是使用相同的系统对一个冲激(小亮点)成像,得到退化的冲激响应。

5.6.3建模估计

​ 由于退化建模能解决图像复原问题，因此多年来一直被 人们使用。在某些情况下，模型甚至可以考虑引起退化的环境条件。建模的另一个主要方法是从基本原理开始推导一个数学模型。

5.7逆滤波

​ 逆滤波是研究由退化函数H退化的图像复原的最初手段。退化函数已给出，或通过一些方法获得。最简单的复原方法是直接做逆滤波，在这里，我们用退化函数除退化图像的傅里叶变换G(u, v)来计算原始图像傅里叶变换的估计$\hat{F}(x,y)$,即
$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}$$
​ 即使知道退化函数，也不能准确地复原未退化的图像[ F(u, v)的傅里叶反变换],因为N(u, v)未知。还有更糟的情况，如果退化函数是零或是非常小的值，则N(u, v)/H(u, v)之比很容易支配估计值$\hat{F}(u,v)$。
$$\hat{F}(u,v)=F(u,v)+\frac{N(u,v)}{H(u,v)}$$
​ 所以一般直接逆滤波的性能是较差的。

5.8 最小均方误差(维纳)滤波

​ 逆滤波方法并没有清楚地说明怎样处理噪声，最小均方误差（维纳）滤波综合了退化函数和噪声统计特征进行复原处理的方法。该方法建立在图像和噪声都是随机变量的基础上，目标是找到未污染图像f的一个估计$\hat{f}$，使它们之间的误差最小。这种误差度量由下式给出：
$$e^2=E\{(f-\hat{f}{)}^2\}$$
​ 其中，E{·}是参数的期望值。这里假设噪声和图像不相关，其中一个或另一个有零均值，且估计中的灰度级是退化图像中灰度级的线性函数。基于这些条件，上式中误差函数的最小值在频率域中由如下表达式给出:
$$\begin{array}{r}\begin{array}{cc}\hat{F}(u,v)& =\left[\frac{H^{\ast }(u,v)S_f(u,v)}{S_f(u,v)|H(u,v){|}^2+S_{\eta }(u,v)}\right]G(u,v)\\ & =\left[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+S_{\eta }(u,v)/S_f(u,v)}\right]G(u,v)\\ & =\left[\frac{1}{H(u,v)}\frac{|H(u,v){|}^2}{|H(u,v){|}^2+S_{\eta }(u,v)/S_f(u,v)}\right]G(u,v)\end{array}\end{array}$$
注：
$$\begin{array}{rl}& H(u,v)=退化函数\\ & H^{\ast }(u,v)=H(u,v)的复共轭\\ & |H(u,v){|}^2=H^{\ast }(u,v)H(u,v)\\ & S_{\eta }(u,v)=|N(u,v){|}^2=噪声的功率谱\\ & S_f(u,v)=|F(u,v){|}^2=未退化图像的功率谱\\ & 如果噪声为零，则噪声功率谱消失，并且维纳滤波简化为逆滤波\end{array}$$
​ 这里，我们用了这样一个事实：一个复数量与其共轭的乘积等于该复数量幅度的平方。这个结果就是众所周知的维纳滤波，该概念由N. Wiener[1942]首次提出。由方括号中的项组成的滤波器通常也称为最小均方误差滤波器或最小二乘误差滤波器。注意上式的第一、二行，维纳滤波器没有逆滤波中退化函数为零的问题，除非对于相
同的u值和v值，整个分母都是零。

5.9 约束最小二乘方滤波

​ 维纳滤波存在一些困难：未退化图像和噪声的功率谱必须是已知的。而约束最小二乘方滤波仅要求噪声方差和均值的知识。且维纳滤波建立在最小化统计准则的基础上，因此在平均意义上它是最优的，而约束最小二乘方滤波，对于其应用到的每一幅图像都能产生最优结果。

5.10 几何均值滤波

​ 对维纳滤波稍加推广，这种推广就是所谓的几何均值滤波器的形式：
$$\hat{F}(u,v)={\left[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2}\right]}^{\alpha }{\left[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+\beta [\frac{S_{\eta }(u,v)}{S_f(u,v)}]}\right]}^{1-\alpha }G(u,v)$$
其中，$\alpha$和$\beta$是正的实常数。几何均值滤波器由两个括号内的幂次分别为$\alpha$和$1-\alpha$的表达式组成。

​ 当$\alpha$=1时，该滤波器退化为逆滤波器;当$\alpha$= 0时，该滤波器变为所谓的参数维纳滤波器，参数维纳滤波器在β = 1时还原为标准的维纳滤波器。如果$\alpha$ = 1/2，则滤波器变成相同幂次的两个量的积，这是几何均值的定义，这样就给出了这种滤波器的命名。当β= 1时，随着c减小到1/2 以下，滤波器的性能越来越接近逆滤波器。类似地，当$\alpha$增大到1/2 以上时，滤波器更接近维纳滤波器。当$\alpha$= 1/2且β=1时，该滤波器通常也称为谱均衡滤波器。

5.11 由投影重建图像

​ 这一节，我们研究从一系列投影重建一幅图像的问题，集中讨论X射线计算机断层（CT）。这是最早且应用最广泛的CT类型，也是当前数字图像处理在医学中的主要应用之一。

5.11.1 引言

​ 最初接触由投影重建图像这块内容的时候是在应用在车牌识别的特征提取，通过车牌在垂直投影下的特征足够其进行不同字符的识别。

5.11.2 计算机断层(CT)原理

​ x射线计算机断层的目的是使用x射线从许多不同的方向穿过物体而得到该物体内部结构的三维描述。传统的胸部X射线透视是把物体放在对X射线敏感的平板对面，并用圆锥形X射线束照射该个体得到的。X射线平板产生一幅图像， 该图像上一个点的亮度与X射线通过该物体后照射到该点上的X射线能量成正比。这幅图像就是在前节中我们讨论过的投影的二维图像。我们可以对整幅图像进行反投影，并创建一个三维物体。 通过在许多角度重复该过程，并把反投影相加，就可产生胸腔结构的三维再现。计算机断层通过身体产生的切片试图得到相同的信息(或者其局部)。一个三维描述可以用堆积这些切片得到。CT的实现要经济得多，因为得到高分辨率切片所需要的检测器数量要比产生相同分辨率的-个完整二维投影所需要的检测器数量少得多。计算负担和X射线的剂量同样要降低，使得一维投影CT成为一种更实际的方法。

​ 第七代(G7)扫描器(也称为多切片CT扫描器)将要问世，这种扫描器使用“厚”的扇形射线束与平行检测器族相配合同时收集体CT数据,即三维横截“厚片层”，而不是每个X射线脉冲产生单一的横截面图像。除了有效地增加了细节之外，这种方法的优点是它使用的X射线管更经济，从而降低了成本并降低了剂量。

5.11.3 投影和雷登变换

​ 在笛卡尔坐标系中，一条直线可以用斜截式 ：$y=ax+b$ 来表示，在投影这方面，我们需要用到它的法线表达式：$xcos\theta +ysin\theta =\rho$

​ 平行射线束的投影可由一组直线建模，投影信号中的任意一点由沿直线的$xcos\theta k+ysin\theta k=\rho j$给出（k，j都是下标，k表示$\theta$每个不同的角度，j表示每个$\rho $不同的接收面），工作在连续变量的情况下，线求和变为变限积分：
$$g(\rho ,\theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\delta (xcos\theta +ysin\theta -\rho )dxdy$$
​ 这个就是雷登变换，雷登变换是投影重建的基石。

5.11.4 傅里叶切片定理

$$
G(ω,θ)=∫_{−∞}
^∞
​
g(ρ,θ)e ^{−j2πωρ}

dρ=F(ωcosθ,ωsinθ)
$$

​ 上式就是众所周知的傅里叶切片定理(或投影切片定理)。它说明一个投影的傅里叶变换是得到投影区域的二维傅里叶变换的一个切片。