youcans 的 OpenCV 学习课—10. 图像复原与重建

本系列面向 Python 小白，从零开始实战解说 OpenCV 项目实战。

图像复原是对图像退化过程建模，并以图像退化的先验知识来恢复退化的图像。

图像增强是一种主观处理，而图像复原是一种客观处理。

本文提供上述各种算法的完整例程和运行结果。

欢迎关注『youcans 的 OpenCV 学习课』系列，持续更新
欢迎关注『youcans 的 OpenCV 例程 200 篇』系列，持续更新中

文章目录

youcans 的 OpenCV 学习课—10. 图像复原与重建
- 1. 图像退化/复原处理模型
- - 1.1 图像退化模型
- 2. 噪声模型
- - 2.1 高斯噪声（Gauss Noise）
  - 例程 9.1：高斯噪声（Gauss Noise）
  - 2.2 瑞利噪声（Rayleigh Noise）
  - 例程 9.2：瑞利噪声（Rayleigh Noise）
  - 2.3 爱尔兰噪声（Ireland Noise）
  - 例程 9.3：伽马噪声（Gamma Noise）
  - 2.4 指数噪声（Exponential Noise）
  - 例程 9.4：指数噪声（Exponent Noise）
  - 2.5 均匀噪声（Uniform Noise）
  - 例程 9.5：均匀噪声（Uniform Noise）
  - 2.6 椒盐噪声（Salt-pepper Noise）
  - 例程 9.6：椒盐噪声（Salt-pepper Noise）
  - 2.7 图像噪声小结
  - 例程 9.7：各种噪声模型的直方图分布
- 在这里插入图片描述
- 3. 仅噪声存在的空间滤波图像复原
- - 3.1 算术平均滤波器（Arithmentic mean filter）
  - 例程 9.8：算术平均滤波器
  - 3.2 几何均值滤波器（Geometric mean filter）
  - 例程 9.9：几何均值滤波器
  - 3.3 谐波平均滤波器（Harmonic mean filter）
  - 例程 9.10：谐波平均滤波器
  - 3.4 反谐波平均滤波器（Invert-harmonic mean filter）
  - 例程 9.11：反谐波平均滤波器
  - 3.5 统计排序滤波器
  - 例程 9.12：统计排序滤波器
  - 3.6 修正阿尔法均值滤波器（Modified alpha-mean filter）
  - 例程 9.13：修正阿尔法均值滤波器
  - 3.7 自适应局部降噪滤波器
  - 例程 9.14：自适应局部降噪滤波器
  - 3.8 自适应中值滤波器（Adaptive median filter）
  - 例程 9.15：自适应中值滤波器
- 4. 频率域滤波图像复原
- - 4.1 陷波滤波器**（Notch Filter）**
  - 例程 9.16：陷波带阻滤波器的传递函数
  - 例程 9.17：陷波带阻滤波器消除周期噪声干扰
  - 4.2 最优陷波滤波器
- 5. 估计退化函数
- - 5.1 观察法估计退化函数
  - 5.2 试验法估计退化函数
  - 5.3 模型法估计退化函数
  - 例程 9.18：运动模糊退化模型
  - 例程 9.19：湍流模糊退化模型
- 6. 退化图像复原
- - 6.1 退化图像的逆滤波（Inverse filter）
  - 例程 9.20：湍流模糊退化图像的逆滤波
  - 6.2 退化图像的维纳滤波（Wiener filter）
  - 例程 9.21: 维纳滤波（Wiener filter）
  - 6.3 约束最小二乘方滤波（Constrained Least Squares Filtering）
  - 例程 9.21: 约束最小二乘方滤波
  - 6.4 几何均值滤波
  - 例程 9.22: 几何均值滤波
- 7. 投影重建图像
- - 7.1 计算机断层成像（Computed tomography）
  - 7.2 投影和雷登变换（Radon transform）
  - 例程 9.23：雷登变换正弦图
  - 7.2 雷登变换反投影重建图像（Radon transform back projection）
  - 例程 9.24：雷登变换反投影重建图像
  - 7.3 滤波反投影重建图像（Filtered back projection）
  - 例程 9.25：滤波反投影重建图像

1. 图像退化/复原处理模型

1.1 图像退化模型

图像退化：

图像退化是指图像在形成、记录、处理和传输过程中，由于成像系统、记录设备、传输介质和处理方法的不完善，导致图像质量下降。

把图像退化表示为退化算子 H\mathcal{H}H，退化算子与一个加性噪声项 η(x,y)\eta (x,y)η(x,y) 共同对输入图像 f(x,y)f(x,y)f(x,y) 进行运算，生成退化图像 g(x,y)g(x,y)g(x,y)。

图像复原：

图像复原是指利用退化过程的先验知识，消除或降低在图像获取、传输过程中造成的图像品质下降，恢复图像的本来面目。具体地，要根据退化的原因，分析引起退化的环境因素，建立相应的数学模型，并沿着使图像降质的逆过程恢复图像。因此，图象复原是一个求逆问题，就是把退化模型化，采用相反的过程处理，以复原出原始图像。

已知退化图像 g(x,y)g(x,y)g(x,y)，关于退化算子 H\mathcal{H}H 与加性噪声项 η(x,y)\eta (x,y)η(x,y) 的概括信息，通过图像复原处理得到输入图像 f(x,y)f(x,y)f(x,y) 的估计 f^(x,y)\hat{f}(x,y)f^(x,y) ，尽可能接近原图像。

一般来说，关于退化算子 H\mathcal{H}H 与加性噪声项 η(x,y)\eta (x,y)η(x,y) 的信息越多，得到的估计 f^(x,y)\hat{f}(x,y)f^(x,y) 就越接近原图像 f(x,y)f(x,y)f(x,y) 。

由于引起退化的因素众多、性质不同，为了描述图像退化过程所建立的数学模型多种多样，而恢复的质量标准也存在差异性，因此图像复原是一个复杂的数学过程，相关的方法和技术也各不相同。

退化模型：

若退化算子 H\mathcal{H}H 是一个线性位置不变算子，则空间域的退化图像为：
g(x,y)=(h⋆f)(x,y)+η(x,y)g(x,y) = (h \star f)(x,y) + \eta(x,y) g(x,y)=(h⋆f)(x,y)+η(x,y)
频率域的等效公式为：
G(u,v)=H(u,v)F(u,v)+N(u,v)G(u,v) = H(u,v) F(u,v) + N(u,v) G(u,v)=H(u,v)F(u,v)+N(u,v)

2. 噪声模型

数字图像中的噪声源主要来自图像获取和传输过程。在获取图像时，光照水平和传感器温度影响图像中的噪声。在传输图像时，传输信道中的干扰对图像产生污染。

白噪声：

白光等比例地包含可见光谱中的所有频率。当噪声的傅里叶谱是常量时，称为白噪声。

2.1 高斯噪声（Gauss Noise）

高斯噪声中空间域和频率域中都很方便进行数学处理，因而得到了广泛的应用。

高斯噪声的概率密度函数为：
p(z)=12πσe−(z−zˉ)2/2σ2p(z) = \frac{1}{\sqrt{2\pi} \sigma} e^{-(z-\bar{z})^2/2\sigma ^2} p(z)=2πσ1e−(z−zˉ)2/2σ2

高斯噪声的均值为 zˉ\bar{z}zˉ，标准差为：σ2\sigma ^2σ2。

例程 9.1：高斯噪声（Gauss Noise）

    # 9.1：高斯噪声 (GaussNoise)img = cv2.imread("../images/Fig0503.tif", 0)  # flags=0 读取为灰度图像# img = np.ones([256, 256]) * 128mu, sigma = 0.0, 20.0noiseGause = np.random.normal(mu, sigma, img.shape)imgGaussNoise = img + noiseGauseimgGaussNoise = np.uint8(cv2.normalize(imgGaussNoise, 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("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()

2.2 瑞利噪声（Rayleigh Noise）

瑞利噪声的概率密度函数为
p(z)={2/b∗(z−a)e−(z−a)2/b,z≥a0,z<ap(z) = \begin{cases} 2/b * (z-a) e^{-(z-a)^2 /b} &, z \ge a\\ 0&, z < a \end{cases} p(z)={2/b∗(z−a)e−(z−a)2/b0,z≥a,z<a

瑞利噪声的均值和标准差为：
zˉ=a+πb/4σ2=b(4−π)/4\bar{z} = a + \sqrt{\pi b/4} \\ \sigma ^2 = b(4-\pi)/4 zˉ=a+πb/4σ2=b(4−π)/4

瑞利噪声概率密度分布到原点的距离及密度的基本形状右偏，常用于倾斜形状直方图的建模。

例程 9.2：瑞利噪声（Rayleigh Noise）

    # # 9.2：瑞利噪声  (RayleighNoise)img = cv2.imread("../images/Fig0503.tif", 0)  # flags=0 读取为灰度图像# img = np.ones([256, 256]) * 128a = 30.0noiseRayleigh = np.random.rayleigh(a, size=img.shape)imgRayleighNoise = img + noiseRayleighimgRayleighNoise = np.uint8(cv2.normalize(imgRayleighNoise, 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("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()

2.3 爱尔兰噪声（Ireland Noise）

爱尔兰噪声的概率密度函数为
p(z)={abzb−1(b−1)!e−az,z≥00,z<0p(z) = \begin{cases} \frac{a^b z^{b-1}}{(b-1)!} e^{-az} &, z \ge 0\\ 0&, z < 0 \end{cases} p(z)={(b−1)!abzb−1e−az0,z≥0,z<0

爱尔兰噪声的均值和标准差为：
zˉ=b/aσ2=b/a2\bar{z} = b/a \\ \sigma ^2 = b/a^2 zˉ=b/aσ2=b/a2

当标准差的分母 a2a^2a2为伽马函数时，称为伽马噪声。

例程 9.3：伽马噪声（Gamma Noise）

    # # 9.3：伽马噪声 (Gamma Noise)img = cv2.imread("../images/Fig0503.tif", 0)  # flags=0 读取为灰度图像# img = np.ones([256, 256]) * 128a, b = 10.0, 2.5noiseGamma = np.random.gamma(shape=b, scale=a, size=img.shape)imgGammaNoise = img + noiseGammaimgGammaNoise = 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()

2.4 指数噪声（Exponential Noise）

指数噪声的概率密度函数为
p(z)={ae−az,z≥00,z<0p(z) = \begin{cases} ae^{-az} &, z \ge 0\\ 0&, z < 0 \end{cases} p(z)={ae−az0,z≥0,z<0

指数噪声的均值和标准差为：
zˉ=1/aσ2=1/a2\bar{z} = 1/a \\ \sigma ^2 = 1/a^2 zˉ=1/aσ2=1/a2

显然，指数噪声是爱尔兰噪声在 b=1 时的特殊情况。

例程 9.4：指数噪声（Exponent Noise）

    # # 9.4：指数噪声 (Exponential noise)img = cv2.imread("../images/Fig0503.tif", 0)  # flags=0 读取为灰度图像# img = np.ones([256, 256]) * 128a = 10.0noiseExponent = np.random.exponential(scale=a, size=img.shape)imgExponentNoise = img + noiseExponentimgExponentNoise = 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()

2.5 均匀噪声（Uniform Noise）

均匀噪声的概率密度函数为
p(z)={1/(b−a),a≤z≤b0,otherp(z) = \begin{cases} 1/(b-a) &, a \le z \le b\\ 0&, other \end{cases} p(z)={1/(b−a)0,a≤z≤b,other

均匀噪声的均值和标准差为：
zˉ=(a+b)/2σ2=(b−a)2/12\bar{z} = (a+b)/2 \\ \sigma ^2 = (b-a)^2/12 zˉ=(a+b)/2σ2=(b−a)2/12

例程 9.5：均匀噪声（Uniform Noise）

    # # 9.5：均匀噪声 (Uniform noise)img = cv2.imread("../images/Fig0503.tif", 0)  # flags=0 读取为灰度图像# img = np.ones([256, 256]) * 128mean, sigma = 10, 100a = 2 * mean - np.sqrt(12 * sigma)  # a = -14.64b = 2 * mean + np.sqrt(12 * sigma)  # b = 54.64noiseUniform = np.random.uniform(a, b, img.shape)imgUniformNoise = img + noiseUniformimgUniformNoise = 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()

2.6 椒盐噪声（Salt-pepper Noise）

椒盐噪声的概率密度函数为
p(z)={Ps,z=2k−1Pp,z=01−(Ps+Pp),z=Vp(z) = \begin{cases} Ps &, z = 2^k -1 \\ Pp &, z = 0 \\ 1-(Ps+Pp)&, z = V \end{cases} p(z)=⎩⎪⎨⎪⎧PsPp1−(Ps+Pp),z=2k−1,z=0,z=V

当 Ps、Pp 都不为 0 时，噪声值是白色的 (2k−1)(2^k-1)(2k−1) 或黑色的 (0)(0)(0)，就像盐粒或胡椒粒那样随机地分布在整个图像上，因此称为椒盐噪声，也称为双极冲击噪声。当 Ps 或 Pp 为 0 时，称为单极冲击噪声。

椒盐噪声的均值和标准差为：
zˉ=(0)Pp+K(1−Ps−Pp)+(2k−1)Psσ2=(0−zˉ)2Pp+(K−zˉ)2(1−Ps−Pp)+(2k−1)2Ps\bar{z} = (0)Pp+K(1-Ps-Pp)+(2^k-1)Ps \\ \sigma ^2 = (0-\bar{z})^2 Pp+(K-\bar{z})^2(1-Ps-Pp)+(2^k-1)^2Ps zˉ=(0)Pp+K(1−Ps−Pp)+(2k−1)Psσ2=(0−zˉ)2Pp+(K−zˉ)2(1−Ps−Pp)+(2k−1)2Ps

像素被白色盐粒、黑色胡椒粒污染的概率 P 称为噪声密度：
P=Ps+PpP = Ps + Pp P=Ps+Pp
例如，Ps=0.05，Pp=0.02，则噪声密度 P=0.07，表示图像中约 5% 的像素被盐粒噪声污染，约 2% 的像素被胡椒粒噪声污染，噪声密度为 7%，即图像中 7% 的像素被椒盐噪声污染。

例程 9.6：椒盐噪声（Salt-pepper Noise）

    # # 9.6：椒盐噪声 (Salt-pepper)img = cv2.imread("../images/Fig0503.tif", 0)  # flags=0 读取为灰度图像ps, pp = 0.05, 0.02mask = np.random.choice((0, 0.5, 1), size=img.shape[:2], p=[pp, (1-ps-pp), ps])imgChoiceNoise = img.copy()imgChoiceNoise[mask==1] = 255imgChoiceNoise[mask==0] = 0plt.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()

2.7 图像噪声小结

上述噪声模型为各种图像噪声的建模提供了有效的工具。例如：

高斯噪声：常用于电子电路及传感器噪声（光照不足和/或高温引起）等因素所导致噪声的建模；

瑞利噪声：常用于距离成像中的噪声建模；

指数噪声和伽马噪声：常用于激光成像中的噪声建模；

冲激噪声：常用于在成像期间的快速瞬变（如开关故障）的噪声建模；

均匀噪声：是对实际情况最基本描述，也是数字图像处理中各种随机数发生器的基础。

例程 9.7：各种噪声模型的直方图分布

    # # 9.7:各种噪声模型的直方图分布img = cv2.imread("../images/Fig0503.tif", 0)  # flags=0 读取为灰度图像mu, sigma = 0.0, 20.0noiseGause = np.random.normal(mu, sigma, img.shape)imgGaussNoise = img + noiseGauseimgGaussNoise = np.uint8(cv2.normalize(imgGaussNoise, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]a = 30.0noiseRayleigh = np.random.rayleigh(a, size=img.shape)imgRayleighNoise = img + noiseRayleighimgRayleighNoise = np.uint8(cv2.normalize(imgRayleighNoise, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]a, b = 10.0, 2.5noiseGamma = np.random.gamma(shape=b, scale=a, size=img.shape)imgGammaNoise = img + noiseGammaimgGammaNoise = np.uint8(cv2.normalize(imgGammaNoise, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]a = 10.0noiseExponent = np.random.exponential(scale=a, size=img.shape)imgExponentNoise = img + noiseExponentimgExponentNoise = np.uint8(cv2.normalize(imgExponentNoise, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]mean, sigma = 10, 100a = 2 * mean - np.sqrt(12 * sigma)  # a = -14.64b = 2 * mean + np.sqrt(12 * sigma)  # b = 54.64noiseUniform = np.random.uniform(a, b, img.shape)imgUniformNoise = img + noiseUniformimgUniformNoise = np.uint8(cv2.normalize(imgUniformNoise, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]ps, pp = 0.05, 0.02mask = np.random.choice((0, 0.5, 1), size=img.shape[:2], p=[pp, (1-ps-pp), ps])imgChoiceNoise = img.copy()imgChoiceNoise[mask==1] = 255imgChoiceNoise[mask==0] = 0plt.figure(figsize=(12, 8))plt.subplot(231), plt.title("Gauss noise")histNP, bins = np.histogram(imgGaussNoise.flatten(), bins=255, range=[0, 255], density=True)plt.bar(bins[:-1], histNP[:])plt.subplot(232), plt.title("Rayleigh noise")histNP, bins = np.histogram(imgRayleighNoise.flatten(), bins=255, range=[0, 255], density=True)plt.bar(bins[:-1], histNP[:])plt.subplot(233), plt.title("Gamma noise")histNP, bins = np.histogram(imgGammaNoise.flatten(), bins=255, range=[0, 255], density=True)plt.bar(bins[:-1], histNP[:])plt.subplot(234), plt.title("Exponential noise")histNP, bins = np.histogram(imgExponentNoise.flatten(), bins=255, range=[0, 255], density=True)plt.bar(bins[:-1], histNP[:])plt.subplot(235), plt.title("Uniform noise")histNP, bins = np.histogram(imgUniformNoise.flatten(), bins=255, range=[0, 255], density=True)plt.bar(bins[:-1], histNP[:])plt.subplot(236), plt.title("Salt-pepper noise")histNP, bins = np.histogram(imgChoiceNoise.flatten(), bins=255, range=[0, 255], density=True)plt.bar(bins[:-1], histNP[:])plt.tight_layout()plt.show()

3. 仅噪声存在的空间滤波图像复原

当一幅图像中唯一存在的退化是噪声时，退化模型简化为：
g(x,y)=f(x,y)+η(x,y)G(u,v)=F(u,v)+N(u,v)g(x,y) = f(x,y) + \eta(x,y) \\ G(u,v) = F(u,v) + N(u,v) g(x,y)=f(x,y)+η(x,y)G(u,v)=F(u,v)+N(u,v)
当仅存在加性随机噪声时，可以采用空间滤波方法来估计原图像 f(x,y)f(x,y)f(x,y)，即对退化图像 g(x,y)g(x,y)g(x,y) 去除噪声。

空间滤波方法在《7. 空间域图像滤波》中进行了详细介绍，本章简要讨论空间滤波的降噪性能。

3.1 算术平均滤波器（Arithmentic mean filter）

算术平均滤波器是最简单的均值滤波器，与空间域滤波中的盒式滤波器相同。

令 SxySxySxy 表示中心在点 (x,y)(x,y)(x,y) 、大小为 m∗nm*nm∗n 的矩形子窗口（邻域）的一组坐标，算术平均滤波器在由 SxySxySxy 定义的区域中，计算被污染图像 g(x,y)g(x,y)g(x,y) 的平均值。

复原的图像 f^\hat{f}f^ 在点 (x,y)(x,y)(x,y) 的值，是使用 SxySxySxy 定义的邻域中的像素算出的算术平均值，即：
f^(x,y)=1mn∑(r,c)∈Sxyg(r,c)\hat{f}(x,y) = \frac{1}{mn} \sum _{(r,c) \in Sxy} g(r,c) f^(x,y)=mn1(r,c)∈Sxy∑g(r,c)
这一操作可以通过一个大小为 m∗nm*nm∗n 的空间滤波器（核）来实现，核的所有系数都是 1/mn1/mn1/mn。

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

例程 9.8：算术平均滤波器

    # 9.8: 算术平均滤波器 (Arithmentic mean filter)# 参见例程 1.70：图像的低通滤波 (盒式滤波器核)img = cv2.imread("../images/Fig0507b.tif", 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()

3.2 几何均值滤波器（Geometric mean filter）

使用几何均值滤波器复原图像，复原图像 f^\hat{f}f^ 在点 (x,y)(x,y)(x,y) 的值是邻域中的像素的几何平均值：
f^(x,y)=[∏(r,c)∈Sxyg(r,c)]1/mn\hat{f}(x,y) = \begin{bmatrix} \prod _{(r,c) \in Sxy} g(r,c) \end{bmatrix} ^{1/mn} f^(x,y)=[∏(r,c)∈Sxyg(r,c)]1/mn
几何均值滤波器实现的平滑与算术平均滤波器相当，但损失的图像细节更少。

例程 9.9：几何均值滤波器

    # 9.9: 几何均值滤波器 (Geometric mean filter)img = cv2.imread("../images/Fig0507b.tif", 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, 3order = 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()

3.3 谐波平均滤波器（Harmonic mean filter）

谐波平均滤波器的复原图像 f^\hat{f}f^ 由下式给出：
f^(x,y)=mn∑(r,c)∈Sxy1/g(r,c)\hat{f}(x,y) = \frac{mn}{\sum _{(r,c) \in Sxy} 1/{g(r,c)}} f^(x,y)=∑(r,c)∈Sxy1/g(r,c)mn
谐波平均滤波器既能处理盐粒噪声（白色噪点），又能处理类似于高斯噪声的其他噪声，但不能处理胡椒噪声（黑色噪点）。

例程 9.10：谐波平均滤波器

    # 9.10: 谐波平均滤波器 (Harmonic mean filter)img = cv2.imread("../images/Fig0507b.tif", 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, 3order = m * nkernalMean = 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-8imgHarMean = 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 / sumTempplt.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()

3.4 反谐波平均滤波器（Invert-harmonic mean filter）

反谐波平均滤波器适用于降低或消除椒盐噪声，复原图像 f^\hat{f}f^ 由下式给出：
f^(x,y)=∑(r,c)∈Sxyg(r,c)Q+1∑(r,c)∈Sxyg(r,c)Q\hat{f}(x,y) = \frac {\sum _{(r,c) \in Sxy} {g(r,c)}^{Q+1}} {\sum _{(r,c) \in Sxy} {g(r,c)}^Q} f^(x,y)=∑(r,c)∈Sxyg(r,c)Q∑(r,c)∈Sxyg(r,c)Q+1
Q 称为滤波器的阶数，Q 取正整数时可以消除胡椒噪声，Q 取负整数时可以消除盐粒噪声，但不能同时消除这两种噪声。使用反谐波平均滤波器必须知道噪声是亮噪声还是暗噪声，据此选择 Q 的正负号。

反谐波平均滤波器当 Q=0Q=0Q=0 时简化为算术平均滤波器；当 Q=−1Q=-1Q=−1 时简化为谐波平均滤波器。

例程 9.11：反谐波平均滤波器

    # 9.11: 反谐波平均滤波器 (Inv-harmonic mean filter)img = cv2.imread("../images/Fig0508a.tif", 0)  # flags=0 读取为灰度图像img_h = img.shape[0]img_w = img.shape[1]m, n = 3, 3order = m * nkernalMean = 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")Q = 1.5  # 反谐波平均滤波器 阶数epsilon = 1e-8imgHarMean = img.copy()imgInvHarMean = img.copy()for i in range(hPad, img_h + hPad):for j in range(wPad, img_w + wPad):# 谐波平均滤波器 (Harmonic mean filter)sumTemp = np.sum(1.0 / (imgPad[i-hPad:i+hPad+1, j-wPad:j+wPad+1] + epsilon))imgHarMean[i-hPad][j-wPad] = order / sumTemp# 反谐波平均滤波器 (Inv-harmonic mean filter)temp = imgPad[i-hPad:i+hPad+1, j-wPad:j+wPad+1] + epsilonimgInvHarMean[i-hPad][j-wPad] = np.sum(np.power(temp, (Q+1))) / np.sum(np.power(temp, Q) + epsilon)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("Harmonic mean filter")plt.imshow(imgHarMean, cmap='gray', vmin=0, vmax=255)plt.subplot(133), plt.axis('off'), plt.title("Invert harmonic mean")plt.imshow(imgInvHarMean, cmap='gray', vmin=0, vmax=255)plt.tight_layout()plt.show()

3.5 统计排序滤波器

统计排序滤波器是空间滤波器，其响应是基于滤波器邻域中的像素值的顺序，排序结果决定了滤波器的输出。

统计排序包括中值滤波器、最大值滤波器、最小值滤波器、中点滤波器和修正阿尔法均值滤波器。

中值滤波器：用预定义的像素邻域中的灰度中值来代替像素的值，与线性平滑滤波器相比能有效地降低某些随机噪声，且模糊度要小得多。

f^(x,y)=median(r,c)∈Sxy{g(r,c)}\hat{f}(x,y) = {median} _{(r,c) \in Sxy} \{g(r,c)\} f^(x,y)=median(r,c)∈Sxy{g(r,c)}

由于需要排序操作，中值滤波消耗的运算时间很长。
OpenCV 提供了 cv.medianBlur 函数实现中值滤波算法，详见《例程 1.73：图像的非线性滤波—中值滤波器》。

最大值滤波器：用预定义的像素邻域中的灰度最大值来代替像素的值，可用于找到图像中的最亮点，或用于消弱与明亮区域相邻的暗色区域，也可以用来降低胡椒噪声。

f^(x,y)=max⁡(r,c)∈Sxy{g(r,c)}\hat{f}(x,y) = \max _{(r,c) \in Sxy} \{g(r,c)\} f^(x,y)=(r,c)∈Sxymax{g(r,c)}

最小值滤波器：用预定义的像素邻域中的灰度最小值来代替像素的值，可用于找到图像中的最暗点，或用于削弱与暗色区域相邻的明亮区域，也可以用来降低盐粒噪声。

f^(x,y)=min⁡(r,c)∈Sxy{g(r,c)}\hat{f}(x,y) = \min _{(r,c) \in Sxy} \{g(r,c)\} f^(x,y)=(r,c)∈Sxymin{g(r,c)}

中点滤波器：用预定义的像素邻域中的灰度的最大值与最小值的均值来代替像素的值，注意中点的取值与中值常常是不同的。中点滤波器是统计排序滤波器与平均滤波器的结合，适合处理随机分布的噪声，例如高斯噪声、均匀噪声。

f^(x,y)=[max⁡(r,c)∈Sxy{g(r,c)}+min⁡(r,c)∈Sxy{g(r,c)}]/2\hat{f}(x,y) = [\max _{(r,c) \in Sxy} \{g(r,c)\} + \min _{(r,c) \in Sxy} \{g(r,c)\}]/2 f^(x,y)=[(r,c)∈Sxymax{g(r,c)}+(r,c)∈Sxymin{g(r,c)}]/2

例程 9.12：统计排序滤波器

    # 9.12: 统计排序滤波器 (Statistical sorting filter)img = cv2.imread("../images/Fig0508a.tif", 0)  # flags=0 读取为灰度图像img_h = img.shape[0]img_w = img.shape[1]m, n = 3, 3kernalMean = 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, 7))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()

程序说明：

需要说明的是，本例程和各种滤波器图像处理结果是为了说明滤波器实现的编程方法和程序运行结果。图中一些滤波器图像处理的效果较差，并不能全面反映该滤波器的性能，只能说明该滤波器不适合处理某些类型的噪声。关于统计滤波器在处理不同噪声的选择和比较，可以参考冈萨雷斯《数字图像处理（第四版）》第五章的相关内容。

3.6 修正阿尔法均值滤波器（Modified alpha-mean filter）

修正阿尔法均值滤波器也属于统计排序滤波器，其思想类似于比赛中去掉最高分和最低分后计算平均分的方法。

令 SxySxySxy 表示中心在点 (x,y)(x,y)(x,y) 、大小为 m∗nm*nm∗n 的矩形子窗口（邻域）的一组坐标，修正阿尔法均值滤波器在由 SxySxySxy 定义的邻域中，删除 d 个最低灰度值和 d 个最高灰度值，计算剩余像素 gR(r,c)g_R(r,c)gR(r,c) 的算术平均值作为输出结果，即：

f^(x,y)=1mn−2d∑(r,c)∈SRgR(r,c)\hat{f}(x,y) = \frac{1}{mn-2d} \sum _{(r,c) \in S_R} g_R(r,c) f^(x,y)=mn−2d1(r,c)∈SR∑gR(r,c)
d 的取值范围是 [0,mn/2−1][0, mn/2-1][0,mn/2−1]。选择 d 的大小对图像处理的效果影响很大，当 d=0d=0d=0 时简化为算术平均滤波器，当 d=mn/2−1d=mn/2-1d=mn/2−1 简化为中值滤波器。d 取其它值时，适合于处理多种混合噪声，如高斯噪声和椒盐噪声。

例程 9.13：修正阿尔法均值滤波器

    # 9.13: 修正阿尔法均值滤波器 (Modified alpha-mean filter)img = cv2.imread("../images/Fig0507b.tif", 0)  # flags=0 读取为灰度图像img_h = img.shape[0]img_w = img.shape[1]m, n = 5, 5kernalMean = 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")imgAlphaFilter0 = np.zeros(img.shape)imgAlphaFilter1 = np.zeros(img.shape)imgAlphaFilter2 = np.zeros(img.shape)for i in range(img_h):for j in range(img_w):# 邻域 m * npad = imgPad[i:i+m, j:j+n]padSort = np.sort(pad.flatten())  # 对邻域像素按灰度值排序d = 1sumAlpha = np.sum(padSort[d:m*n-d-1])  # 删除 d 个最大灰度值, d 个最小灰度值imgAlphaFilter0[i, j] = sumAlpha / (m*n-2*d)  # 对剩余像素进行算术平均d = 2sumAlpha = np.sum(padSort[d:m*n-d-1])imgAlphaFilter1[i, j] = sumAlpha / (m*n-2*d)d = 4sumAlpha = np.sum(padSort[d:m*n-d-1])imgAlphaFilter2[i, j] = sumAlpha / (m*n-2*d)plt.figure(figsize=(9, 7))plt.subplot(221), plt.axis('off'), plt.title("Original")plt.imshow(img, cmap='gray', vmin=0, vmax=255)plt.subplot(222), plt.axis('off'), plt.title("Modified alpha-mean(d=1)")plt.imshow(imgAlphaFilter0, cmap='gray', vmin=0, vmax=255)plt.subplot(223), plt.axis('off'), plt.title("Modified alpha-mean(d=2)")plt.imshow(imgAlphaFilter1, cmap='gray', vmin=0, vmax=255)plt.subplot(224), plt.axis('off'), plt.title("Modified alpha-mean(d=4)")plt.imshow(imgAlphaFilter2, cmap='gray', vmin=0, vmax=255)plt.tight_layout()plt.show()

3.7 自适应局部降噪滤波器

前述滤波器直接应用到图像处理，并未考虑图像本身的特征。自适应滤波器的特性根据 m∗nm*nm∗n 矩形邻域 SxySxySxy 定义的滤波区域内的图像的统计特性变化。通常，自适应滤波器的性能优于前述的滤波器，但滤波器的复杂度和计算量也更大。

均值和方差是随机变量最基础的统计量。在图像处理中，均值是像素邻域的平均灰度，方差是像素邻域的图像对比度。

令 SxySxySxy 表示中心在点 (x,y)(x,y)(x,y) 、大小为 m∗nm*nm∗n 的矩形子窗口（邻域），滤波器在由 SxySxySxy 定义的邻域操作。

噪声图像在点 (x,y)(x,y)(x,y) 的值 g(x,y)g(x,y)g(x,y)，噪声的方差 ση2\sigma^2_{\eta}ση2 由噪声图像估计；SxySxySxy 中像素的局部平均灰度为 zˉSxy\bar{z}_{Sxy}zˉSxy，灰度的局部方差为 σSxy2\sigma^2_{Sxy}σSxy2。

使用自适应局部降噪滤波器复原的图像 f^\hat{f}f^ 在点 (x,y)(x,y)(x,y) 的值，由如下自适应表达式描述：
f^(x,y)=g(x,y)−ση2σSxy2[g(x,y)−zSxyˉ]\hat{f}(x,y) = g(x,y) - \frac{\sigma^2_{\eta}}{\sigma^2_{Sxy}} [g(x,y)-\bar{z_{Sxy}}] f^(x,y)=g(x,y)−σSxy2ση2[g(x,y)−zSxyˉ]

例程 9.14：自适应局部降噪滤波器

    # 9.14: 自适应局部降噪滤波器 (Adaptive local noise reduction filter)img = cv2.imread("../images/Fig0507b.tif", 0)  # flags=0 读取为灰度图像hImg = img.shape[0]wImg = img.shape[1]m, n = 5, 5imgAriMean = cv2.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")# 估计原始图像的噪声方差 sigmaEtamean, stddev = cv2.meanStdDev(img)sigmaEta = stddev ** 2print(sigmaEta)# 自适应局部降噪epsilon = 1e-8imgAdaLocal = np.zeros(img.shape)for i in range(hImg):for j in range(wImg):pad = imgPad[i:i+m, j:j+n]  # 邻域 Sxy, m*ngxy = img[i,j]  # 含噪声图像的像素点zSxy = np.mean(pad)  # 局部平均灰度sigmaSxy = np.var(pad)  # 灰度的局部方差rateSigma = min(sigmaEta / (sigmaSxy + epsilon), 1.0)  # 加性噪声假设：sigmaEta/sigmaSxy < 1imgAdaLocal[i, j] = gxy - rateSigma * (gxy - zSxy)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("Adaptive local filter")plt.imshow(imgAdaLocal, cmap='gray', vmin=0, vmax=255)plt.tight_layout()plt.show()

3.8 自适应中值滤波器（Adaptive median filter）

中值滤波器的窗口尺寸是固定大小不变的，不能同时兼顾去噪和保护图像的细节，在噪声的密度较小时的性能较好，当噪声概率较高时的性能就会劣化。

自适应中值滤波器根据预先设定的条件，在滤波的过程中动态改变滤波器的窗口尺寸大小；进一步地，根据条件判断当前像素是否噪声，由此决定是否用邻域中值替换当前像素。

自适应中值滤波器可以处理较大概率的脉冲噪声，平滑非脉冲噪声，尽可能保护图像细节信息，避免图像边缘的细化或者粗化。

令 SxySxySxy 表示中心在点 (x,y)(x,y)(x,y) 、大小为 m∗nm*nm∗n 的矩形子窗口（邻域），滤波器在由 SxySxySxy 定义的邻域操作。

Zmin、Zmax、ZmedZmin、Zmax、ZmedZmin、Zmax、Zmed 表示 SxySxySxy 中的最小灰度值、最大灰度值和灰度值的中值，ZxyZxyZxy 是点 (x,y)(x,y)(x,y) 的灰度值，SmaxSmaxSmax 是允许的最大窗口尺寸。

自适应中值滤波器分为两个过程：

Step A：
- A1 = Zmed - Zmin
- A2 = Zmed - Zmax
- 如果 A1>0 且 A2<0，则跳转到 Step B；否则，增大窗口尺寸
- 如果增大后的尺寸 ≤ Smax，则重复 A；否则，输出 Zmed
Step B：
- B1 = Zxy - Zmin
- B2 = Zxy - Zmax
- 如果 B1>0 且 B2<0，则输出 Zxy；否则，输出 Zmed

StepA 的目的是确定当前窗口内得到的中值 Zmed 是否是噪声。如果 Zmin<Zmed<Zmax，则中值 Zmed 不是噪声，转到StepB。如果Zmin<Zxy<Zmax，则 Zxy 不是噪声，滤波器输出 Zxy。如果不满足上述条件，则判定 Zxy 是噪声，输出中值 Zmed。

如果在StepA中，Zmed 不满足条件 Zmin<Zmed<Zmax，则可判断得到的中值 Zmed 是噪声。这时要增大滤波器的窗口尺寸，在更大的范围内寻找一个非噪声点的中值。

因此，如果图像中噪声的概率较低，自适应中值滤波器可以使用较小的窗口尺寸，以提高计算速度；反之，如果噪声的概率较高，则需要增大滤波器的窗口尺寸，以改善滤波效果。

例程 9.15：自适应中值滤波器

    # # 9.15: 自适应中值滤波器 (Adaptive median filter)img = cv2.imread("../images/Fig0514a.tif", 0)  # flags=0 读取为灰度图像hImg = img.shape[0]wImg = img.shape[1]smax = 7  # 允许最大窗口尺寸m, n = smax, smaximgAriMean = cv2.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")imgMedianFilter = np.zeros(img.shape)  # 中值滤波器imgAdaMedFilter = np.zeros(img.shape)  # 自适应中值滤波器for i in range(hPad, hPad+hImg):for j in range(wPad, wPad+wImg):# 1. 中值滤波器 (Median filter)ksize = 3k = int(ksize/2)pad = imgPad[i-k:i+k+1, j-k:j+k+1]  # 邻域 Sxy, m*nimgMedianFilter[i-hPad, j-wPad] = np.median(pad)# 2. 自适应中值滤波器 (Adaptive median filter)ksize = 3k = int(ksize/2)pad = imgPad[i-k:i+k+1, j-k:j+k+1]zxy = img[i-hPad][j-wPad]zmin = np.min(pad)zmed = np.median(pad)zmax = np.max(pad)if zmin < zmed < zmax:if zmin < zxy < zmax:imgAdaMedFilter[i-hPad, j-wPad] = zxyelse:imgAdaMedFilter[i-hPad, j-wPad] = zmedelse:while True:ksize = ksize + 2if zmin < zmed < zmax or ksize > smax:breakk = int(ksize / 2)pad = imgPad[i-k:i+k+1, j-k:j+k+1]zmed = np.median(pad)zmin = np.min(pad)zmax = np.max(pad)if zmin < zmed < zmax or ksize > smax:if zmin < zxy < zmax:imgAdaMedFilter[i-hPad, j-wPad] = zxyelse:imgAdaMedFilter[i-hPad, j-wPad] = zmedplt.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("Median filter")plt.imshow(imgMedianFilter, cmap='gray', vmin=0, vmax=255)plt.subplot(133), plt.axis('off'), plt.title("Adaptive median filter")plt.imshow(imgAdaMedFilter, cmap='gray', vmin=0, vmax=255)plt.tight_layout()plt.show()

4. 频率域滤波图像复原

通过频率域滤波可以有效分析并滤除周期噪声，其理论基础是傅里叶变换后周期噪声在对应周期干扰的频率显示为集中突发的能量，因此可以使用选择性滤波器来分离滤除噪声。

4.1 陷波滤波器（Notch Filter）

陷波滤波器阻止或通过预定的频率矩形邻域中的频率，可以很好地复原被周期性噪声干扰的图像。在《5.2 陷波滤波器（Notch Filter）》已进行介绍并给出了例程。

陷波滤波器可以在某一个频率点迅速衰减输入信号，以达到阻碍此频率信号通过的滤波效果的滤波器。

陷波带阻滤波器的传递函数是中心平移到陷波中心的各个高通滤波器的乘积：
HNR(u,v)=∏k=1QHk(u,v)H−k(u,v)H_{NR}(u,v) = \prod_{k=1}^Q H_k(u,v) H_{-k}(u,v) HNR(u,v)=k=1∏QHk(u,v)H−k(u,v)

其中，滤波器的距离计算公式为：
Dk(u,v)=(u−M/2−uk)2+(v−N/2−vk)2D−k(u,v)=(u−M/2+uk)2+(v−N/2+vk)2D_k(u,v) = \sqrt{(u-M/2-u_k)^2 + (v-N/2-v_k)^2} \\ D_{-k}(u,v) = \sqrt{(u-M/2+u_k)^2 + (v-N/2+v_k)^2} Dk(u,v)=(u−M/2−uk)2+(v−N/2−vk)2D−k(u,v)=(u−M/2+uk)2+(v−N/2+vk)2
例如，具有 3个陷波对的 n 阶巴特沃斯陷波带阻滤波器为：
HNR(u,v)=∏k=13[11+[D0k/Dk(u,v)]n][11+[D−k/Dk(u,v)]n]H_{NR}(u,v) = \prod_{k=1}^3 [\frac{1}{1+[D_{0k}/D_k(u,v)]^n}] [\frac{1}{1+[D_{-k}/D_k(u,v)]^n}] HNR(u,v)=k=1∏3[1+[D0k/Dk(u,v)]n1][1+[D−k/Dk(u,v)]n1]

例程 9.16：陷波带阻滤波器的传递函数

    # 9.16: 陷波带阻滤波器的传递函数def ideaBondResistFilter(shape, radius=10, w=5):  # 理想带阻滤波器u, v = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))D = np.sqrt((u - shape[1]//2)**2 + (v - shape[0]//2)**2)D0 = radiushalfW = w/2kernel = np.piecewise(D, [D<=D0+halfW, D<=D0-halfW], [1, 0])kernel = 1 - kernel  # 带阻return kerneldef butterworthNRFilter(img, radius=10, uk=60, vk=80, n=1):  # 巴特沃斯陷波带阻滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))Dkm = np.sqrt((u - M//2 - uk)**2 + (v - N//2 - vk)**2)  # D_+kDkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-kD0 = radiusn2 = n * 2epsilon = 1e-6kernel = (1 / (1 + (D0 / (Dkm+epsilon))**n2)) * (1 / (1 + (D0 / (Dkp+epsilon))**n2))return kerneldef gaussNRFilter(img, radius=10, uk=60, vk=80):  # 高斯陷波带阻滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))Dkm = np.sqrt((u - M//2 - uk)**2 + (v - N//2 - vk)**2)  # D_+kDkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-kD0 = radiuskernel = (1 - np.exp(-(Dkm**2)/(D0**2))) * (1 - np.exp(-(Dkp**2)/(D0**2)))return kerneldef ideaNRFilter(img, radius=10, uk=60, vk=80):  # 高斯陷波带阻滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))Dkm = np.sqrt((u - M//2 - uk)**2 + (v - N//2 - vk)**2)  # D_+kDkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-kD0 = radiusk1 = Dkm.copy()k1[Dkm>D0] = 1k1[Dkm<=D0] = 0k2 = Dkp.copy()k2[Dkp>D0] = 1k2[Dkp<=D0] = 0kernel = k1 * k2return kernel# 理想、高斯、巴特沃斯陷波带阻滤波器传递函数img = np.zeros([128, 128])shape = img.shaperadius = 15INRF = ideaNRFilter(img, radius=radius, uk=20, vk=30)  # (uk,vk) 陷波中心GNRF = gaussNRFilter(img, radius=radius, uk=20, vk=30)BNRF = butterworthNRFilter(img, radius=radius, uk=20, vk=30, n=2)filters = ["INRF", "GNRF", "BNRF"]u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]fig = plt.figure(figsize=(10, 8))for i in range(3):nrFilter = eval(filters[i]).copy()ax1 = fig.add_subplot(3, 3, 3*i+1)ax1.imshow(nrFilter, 'gray')ax1.set_title(filters[i]), ax1.set_xticks([]), ax1.set_yticks([])ax2 = plt.subplot(3,3,3*i+2, projection='3d')ax2.set_title("transfer function")ax2.plot_wireframe(u, v, nrFilter, rstride=2, linewidth=0.5, color='c')ax2.set_xticks([]), ax2.set_yticks([]), ax2.set_zticks([])ax3 = plt.subplot(3,3,3*i+3)profile = nrFilter[:, 30]ax3.plot(profile), ax3.set_title("profile"), ax3.set_xticks([]), ax3.set_yticks([])plt.show()

例程 9.17：陷波带阻滤波器消除周期噪声干扰

    # 9.17: 陷波带阻滤波器消除周期噪声干扰def butterworthNRFilter(img, radius=10, uk=10, vk=10, n=2):  # 巴特沃斯陷波带阻滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))Dm = np.sqrt((u - M//2 - uk)**2 + (v - N//2 - vk)**2)Dp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)D0 = radiusn2 = 2 * nkernel = (1 / (1 + (D0 / (Dm + 1e-6))**n2)) * (1 / (1 + (D0 / (Dp + 1e-6))**n2))return kernel# (1) 读取原始图像img = cv2.imread("../images/Fig0505a.tif", flags=0)  # flags=0 读取为灰度图像imgFloat32 = np.float32(img)  # 将图像转换成 float32rows, cols = img.shape[:2]  # 图片的高度和宽度fig = plt.figure(figsize=(9, 6))plt.subplot(231), plt.title("Original image"), plt.axis('off'), plt.imshow(img, cmap='gray')# (2) 中心化, centralized 2d array f(x,y) * (-1)^(x+y)mask = np.ones(img.shape)mask[1::2, ::2] = -1mask[::2, 1::2] = -1fImage = imgFloat32 * mask  # f(x,y) * (-1)^(x+y)# (3) 快速傅里叶变换rPadded = cv2.getOptimalDFTSize(rows)  # 最优 DFT 扩充尺寸cPadded = cv2.getOptimalDFTSize(cols)  # 用于快速傅里叶变换dftImage = np.zeros((rPadded, cPadded, 2), np.float32)  # 对原始图像进行边缘扩充dftImage[:rows, :cols, 0] = fImage  # 边缘扩充，下侧和右侧补0cv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换dftAmp = cv2.magnitude(dftImage[:,:,0], dftImage[:,:,1])  # 傅里叶变换的幅度谱 (rPad, cPad)dftAmpLog = np.log(1.0 + dftAmp)  # 幅度谱对数变换，以便于显示dftAmpNorm = np.uint8(cv2.normalize(dftAmpLog, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]plt.subplot(232), plt.axis('off'), plt.title("DFT spectrum")plt.imshow(dftAmpNorm, cmap='gray')plt.arrow(445, 370, 25, 30, width=5, length_includes_head=True, shape='full')  # 在图像上加上箭头plt.arrow(550, 490, -25, -30, width=5, length_includes_head=True, shape='full')  # 在图像上加上箭头# (4) 构建陷波带阻滤波器 传递函数BRFilter = butterworthNRFilter(dftImage, radius=15, uk=25, vk=16, n=3)  # 巴特沃斯陷波带阻滤波器, 处理周期噪声plt.subplot(233), plt.axis('off'), plt.title("Butterworth notch resist filter")plt.imshow(BRFilter, cmap='gray')# (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 陷波带阻滤波器dftFilter = np.zeros(dftImage.shape, dftImage.dtype)  # 快速傅里叶变换的尺寸(优化尺寸)for i in range(2):dftFilter[:rPadded, :cPadded, i] = dftImage[:rPadded, :cPadded, i] * BRFilter# 频域滤波傅里叶变换的傅里叶谱nrfDftAmp = cv2.magnitude(dftFilter[:, :, 0], dftFilter[:, :, 1])  # 傅里叶变换的幅度谱nrfDftAmpLog = np.log(1.0 + nrfDftAmp)  # 幅度谱对数变换，以便于显示nrfDftAmpNorm = np.uint8(cv2.normalize(nrfDftAmpLog, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]plt.subplot(234), plt.axis('off'), plt.title("BNRF DFT Spectrum")plt.imshow(nrfDftAmpNorm, cmap='gray')# (6) 对频域滤波傅里叶变换 执行傅里叶逆变换，并只取实部idft = np.zeros(dftAmp.shape, np.float32)  # 快速傅里叶变换的尺寸(优化尺寸)cv2.dft(dftFilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)# (7) 中心化, centralized 2d array g(x,y) * (-1)^(x+y)mask2 = np.ones(dftAmp.shape)mask2[1::2, ::2] = -1mask2[::2, 1::2] = -1idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)plt.subplot(235), plt.axis('off'), plt.title("g(x,y)*(-1)^(x+y)")plt.imshow(idftCen, cmap='gray')# (8) 截取左上角，大小和输入图像相等idftCenClip = np.clip(idftCen, 0, 255)  # 截断函数，将数值限制在 [0,255]imgFiltered = idftCenClip.astype(np.uint8)imgFiltered = imgFiltered[:rows, :cols]plt.subplot(236), plt.axis('off'), plt.title("BNRF filtered image")plt.imshow(imgFiltered, cmap='gray')plt.tight_layout()plt.show()print("image.shape:{}".format(img.shape))print("imgFloat32.shape:{}".format(imgFloat32.shape))print("dftImage.shape:{}".format(dftImage.shape))print("dftAmp.shape:{}".format(dftAmp.shape))print("idft.shape:{}".format(idft.shape))print("dftFilter.shape:{}".format(dftFilter.shape))print("imgFiltered.shape:{}".format(imgFiltered.shape))

程序说明：

傅立叶变换的频谱反映能量分布。通过中心化移频到原点以后，傅里叶变换的频谱图是以原点为中心对称分布的。中心化不仅可以清晰地看出图像频率分布，还可以分离出有周期性规律的干扰信号。

本例程中的图像受到正弦噪声的干扰，从中心化的频谱图可以看出，除原点以外还存在一对对称分布的亮点（箭头指示处），这就干扰噪声产生的。因此，在该亮点位置设计陷波带阻滤波器，可以消除干扰正弦噪声。

4.2 最优陷波滤波器

当存在多个干扰分量时，陷波滤波器的设计比较困难，而且容易过多地滤除图像信息。特别是当干扰分量不是单频率突发能量，而是携带干扰模式信息的宽裙摆，在正常的变换背景中就难以检测和分离。

最优陷波滤波方法在最小化复原估计的局部方差方面是最优的，该方法首先分离干扰模式的各个主要贡献，然后从被污染图像中减去该模式的一个可变加权部分。

5. 估计退化函数

估计图像复原中所用的退化函数，主要有三种方法：观察法、试验法和数学建模方法。

5.1 观察法估计退化函数

为了降低噪声的影响，选择一个信号内容很强的区域（如一个高对比度区域）。

令gs(x,y)g_s(x,y)gs(x,y) 表示观察子图像，f^s(x,y)\hat{f}_s(x,y)f^s(x,y) 表示处理后的子图像，则：
Hs(u,v)=Gs(u,v)f^s(u,v)H_s(u,v) = \frac{G_s(u,v)}{\hat{f}_s(u,v)} Hs(u,v)=f^s(u,v)Gs(u,v)
由该函数的特性，可以根据位置不变假设来推断完整的退化函数 H(u,v)H(u,v)H(u,v)。

5.2 试验法估计退化函数

如果设备与获取退化图像的设备相似，原则上就可能获得退化模型的精确估计——如果能使用原设备当然更好。在各种系统设置下，获取类似于退化图像的图像，然后再该系统设置条件下对一个冲激（小光点）成像，就得到退化的冲激响应。
H(u,v)=G(u,v)AH(u,v) = \frac{G(u,v)}{A} H(u,v)=AG(u,v)
G(u,v)G(u,v)G(u,v) 是观察图像的傅里叶变换，A 是描述冲激强度的常量。

5.3 模型法估计退化函数

分析导致退化的原因，根据基本原理提出退化模型，如湍流导致的模糊、匀速运动导致的模糊，可以基于模型更加准确地估计退化函数。

下面以运动模糊和大气湍流模型为例，采用退化模型对图像的退化建模。

例程 9.18：运动模糊退化模型

运动模糊是相机，物体，背景间相对运动造成的效果，通常由于长时间曝光或场景内的物体快速移动导致，在摄影中可以借助移动镜头追踪移动的物体来避免。

对匀速线性运动模糊建模，假设图像 f(x,y)f(x,y)f(x,y) 做平面运动，运动在 x、y 方向的时变分量分别为： x0(t)=at/Tx_0(t)=at/Tx0(t)=at/T、y0(t)=bt/Ty_0(t)=bt/Ty0(t)=bt/T。记录介质上任意点的总曝光量是瞬时曝光量的积分，可以建立运动模糊退化函数模型：
H(u,v)=Tπ(ua+vb)sin[π(ua+vb)]e−jπ(ua+vb)H(u,v) = \frac{T}{\pi (ua+vb)} sin[\pi(ua+vb)]e^{-j \pi (ua+vb)} H(u,v)=π(ua+vb)Tsin[π(ua+vb)]e−jπ(ua+vb)

    # 9.18: 运动模糊退化图像 (Motion blur degradation)def motionBlur(image, degree=10, angle=45):image = np.array(image)center = (degree/2, degree/2)  # 旋转中心M = cv2.getRotationMatrix2D(center, angle, 1)  # 无损旋转kernel = np.diag(np.ones(degree) / degree)  # 运动模糊内核kernel = cv2.warpAffine(kernel, M, (degree, degree))blurred = cv2.filter2D(image, -1, kernel)  # 图像卷积blurredNorm = np.uint8(cv2.normalize(blurred, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]return blurredNorm# 运动模糊图像img = cv2.imread("../images/Fig0526a.tif", 0)  # flags=0 读取为灰度图像imgBlur1 = motionBlur(img, degree=30, angle=45)imgBlur2 = motionBlur(img, degree=40, angle=45)imgBlur3 = motionBlur(img, degree=60, angle=45)plt.figure(figsize=(9, 6))plt.subplot(131), plt.title("degree=20"), plt.axis('off'), plt.imshow(imgBlur1, 'gray')plt.subplot(132), plt.title("degree=40"), plt.axis('off'), plt.imshow(imgBlur2, 'gray')plt.subplot(133), plt.title("degree=60"), plt.axis('off'), plt.imshow(imgBlur3, 'gray')plt.tight_layout()plt.show()

例程 9.19：湍流模糊退化模型

湍流是自然界中普遍存在的一种复杂的流动现象。物体通过湍流大气成像时，受到湍流效应的影响，出现光强闪烁、光束方向漂移、光束宽度扩展及接收面上相位的起伏，造成图像模糊和抖动，甚至扭曲变形。

Hufnagel and Stanley 根据大气湍流的物理特性提出一种退化模型：

H(u,v)=e−k(u2+v2)5/6H(u,v) = e^{-k (u^2+v^2)^{5/6}} H(u,v)=e−k(u2+v2)5/6

k 是湍流常数，反映湍流的强烈程度。

    # 9.19: 湍流模糊退化模型 (turbulence blur degradation model)def getDegradedImg(image, Huv):  # 根据退化模型生成退化图像rows, cols = image.shape[:2]  # 图片的高度和宽度# (1) 中心化, centralized 2d array f(x,y) * (-1)^(x+y)mask = np.ones((rows, cols))mask[1::2, ::2] = -1mask[::2, 1::2] = -1imageCen = image * mask# (2) 快速傅里叶变换dftImage = np.zeros((rows, cols, 2), np.float32)dftImage[:, :, 0] = imageCencv2.dft(dftImage, dftImage, cv2.DFT_COMPLEX_OUTPUT)  # 快速傅里叶变换 (rows, cols, 2)# (4) 构建 频域滤波器传递函数:Filter = np.zeros((rows, cols, 2), np.float32)  # (rows, cols, 2)Filter[:, :, 0], Filter[:, :, 1] = Huv, Huv# (5) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数dftFilter = dftImage * Filter# (6) 对修正傅里叶变换 进行傅里叶逆变换，并只取实部idft = np.ones((rows, cols), np.float32)  # 快速傅里叶变换的尺寸cv2.dft(dftFilter, idft, cv2.DFT_REAL_OUTPUT + cv2.DFT_INVERSE + cv2.DFT_SCALE)  # 只取实部# (7) 中心化, centralized 2d array g(x,y) * (-1)^(x+y)mask2 = np.ones(dftImage.shape[:2])mask2[1::2, ::2] = -1mask2[::2, 1::2] = -1idftCen = idft * mask2  # g(x,y) * (-1)^(x+y)# (8) 截取左上角，大小和输入图像相等imgDegraded = np.uint8(cv2.normalize(idftCen, None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]# print(image.shape, dftFilter.shape, imgDegraded.shape)return imgDegradeddef turbulenceBlur(img, k=0.001):  # 湍流模糊传递函数# H(u,v) = exp(-k(u^2+v^2)^5/6)M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))radius = (u - M//2)**2 + (v - N//2)**2kernel = np.exp(-k * np.power(radius, 5/6))return kernel# 读取原始图像img = cv2.imread("../images/Fig0525a.tif", 0)  # flags=0 读取为灰度图像# 生成湍流模糊图像HBlur1 = turbulenceBlur(img, k=0.001)  # 湍流模糊传递函数imgBlur1 = getDegradedImg(img, HBlur1)  # 生成湍流模糊图像HBlur2 = turbulenceBlur(img, k=0.0025)imgBlur2 = getDegradedImg(img, HBlur2)plt.figure(figsize=(9, 6))plt.subplot(131), plt.title("origin"), plt.axis('off'), plt.imshow(img, 'gray')plt.subplot(132), plt.title("turbulence blur(k=0.001)"), plt.axis('off'), plt.imshow(imgBlur1, 'gray')plt.subplot(133), plt.title("turbulence blur(k=0.0025)"), plt.axis('off'), plt.imshow(imgBlur2, 'gray')plt.tight_layout()plt.show()

6. 退化图像复原

图像复原是对图像退化的过程进行估计，并补偿退化过程造成的失真，以便获得未经退化的原始图像或原始图像的最优估值，从而改善图像质量的一种方法。

典型的图像复原方法是根据图像退化的先验知识建立退化模型，以退化模型为基础采用滤波等手段进行处理，使复原后的图像符合一定的准则，达到改善图像质量的目的。

因此，图像复原是沿着质量降低的逆过程来重现真实的原始图像，通过去模糊函数而去除图像模糊。

6.1 退化图像的逆滤波（Inverse filter）

图像退化表示为退化算子 H\mathcal{H}H，退化函数可以用观察法、试验法或建模法估计，则通过逆滤波就可以直接实现图像复原。用退化图像的傅里叶变换除以退化函数的傅里叶变换，得到原始图像的傅里叶变换估计：
F^(u,v)=G(u,v)H(u,v)\hat{F}(u,v) = \frac{G(u,v)}{H(u,v)} F^(u,v)=H(u,v)G(u,v)

但是，由于实际上退化图像是退化算子与加性噪声项共同作用的结果，由此得到：
F^(u,v)=F(u,v)+N(u,v)H(u,v)\hat{F}(u,v) = F(u,v) + \frac{N(u,v)}{H(u,v)} F^(u,v)=F(u,v)+H(u,v)N(u,v)

这表明即使获得退化函数 H(u,v)H(u,v)H(u,v) 的估计，由于噪声项是未知的，因此也不能准确地复原原始图像。

进一步地，如果退化函数为 0 或很小，则噪声项的影响将非常严重（信噪比低）。这时，需要将频率限制到原点附近进行分析，可以减少遇到零值的可能性。

例程 9.20：湍流模糊退化图像的逆滤波

如前所述，通过湍流退化模型可以得到退化图像。使用该退化模型进行逆滤波，退化函数与生成退化图像所用的退化函数相反：
H(u,v)=e−k[(u−M/2)2+(v−N/2)2]5/6H(u,v) = e^{-k [(u-M/2)^2+(v-N/2)^2]^{5/6}} H(u,v)=e−k[(u−M/2)2+(v−N/2)2]5/6

但是，直接使用退化模型 H(u,v) 逆滤波的结果（D0=full）很差，用理想低通滤波器对退化模型 H(u,v) 在半径 D0 之外截止后，则视觉效果较好。

    # 9.20: 湍流模糊退化图像的逆滤波def turbulenceBlur(img, k=0.001):  # 湍流模糊传递函数: H(u,v) = exp(-k(u^2+v^2)^5/6)M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))radius = (u - M//2)**2 + (v - N//2)**2kernel = np.exp(-k * np.power(radius, 5/6))return kerneldef getDegradedImg(image, Huv, eps):  # 根据退化模型生成退化图像# (1) 傅里叶变换, 中心化fft = np.fft.fft2(image.astype(np.float32))  # 傅里叶变换fftShift = np.fft.fftshift(fft)  # 将低频分量移动到频域图像中心# (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数fftShiftFilter = fftShift * Huv  # Guv = Fuv * Huv# (3) 对修正傅里叶变换 进行傅里叶逆变换，逆中心化invShift = np.fft.ifftshift(fftShiftFilter)  # 将低频分量逆转换回图像四角imgIfft = np.fft.ifft2(invShift)  # 逆傅里叶变换，返回值是复数数组imgDegraded = np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]return imgDegradeddef ideaLPFilter(img, radius=10):  # 理想低通滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))D = np.sqrt((u - M//2)**2 + (v - N//2)**2)kernel = np.zeros(img.shape[:2], np.float32)kernel[D <= radius] = 1return kerneldef inverseFilter(image, Huv, D0):  # 根据退化模型逆滤波# (1) 傅里叶变换, 中心化fft = np.fft.fft2(image.astype(np.float32))  # 傅里叶变换fftShift = np.fft.fftshift(fft)  # 将低频分量移动到频域图像中心# (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数if D0==0:fftShiftFilter = fftShift / Huv  # Guv = Fuv / Huvelse:lpFilter = ideaLPFilter(image, radius=D0)fftShiftFilter = fftShift / Huv * lpFilter  # Guv = Fuv / Huv# (3) 对修正傅里叶变换 进行傅里叶逆变换，逆中心化invShift = np.fft.ifftshift(fftShiftFilter)  # 将低频分量逆转换回图像四角imgIfft = np.fft.ifft2(invShift)  # 逆傅里叶变换，返回值是复数数组imgRebuild = np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]return imgRebuild# 读取原始图像img = cv2.imread("../images/Fig0525a.tif", 0)  # flags=0 读取为灰度图像# 生成湍流模糊图像HTurb = turbulenceBlur(img, k=0.0025)imgBlur = np.abs(getDegradedImg(img, HTurb, 0.0))print(imgBlur.max(), imgBlur.min())# # 逆滤波imgRebuild = inverseFilter(imgBlur, HTurb, 480)  # Huv 全滤波器imgRebuild1 = inverseFilter(imgBlur, HTurb, D0=40)  # 在半径 D0 之外 Huv 截止imgRebuild2 = inverseFilter(imgBlur, HTurb, D0=70)imgRebuild3 = inverseFilter(imgBlur, HTurb, D0=100)plt.figure(figsize=(9, 7))plt.subplot(231), plt.title("origin"), plt.axis('off'), plt.imshow(img, 'gray')plt.subplot(232), plt.title("turbulence blur"), plt.axis('off'), plt.imshow(imgBlur, 'gray')plt.subplot(233), plt.title("inverse filter(D0=full)"), plt.axis('off'), plt.imshow(imgRebuild, 'gray')plt.subplot(234), plt.title("inverse filter(D0=40)"), plt.axis('off'), plt.imshow(imgRebuild1, 'gray')plt.subplot(235), plt.title("inverse filter(D0=70)"), plt.axis('off'), plt.imshow(imgRebuild2, 'gray')plt.subplot(236), plt.title("inverse filter(D0=100)"), plt.axis('off'), plt.imshow(imgRebuild3, 'gray')plt.tight_layout()plt.show()

6.2 退化图像的维纳滤波（Wiener filter）

逆滤波对加性噪声特别敏感，使得恢复的图像几乎不可用（例程9.20中的退化图像未加入噪声项）。

最小均方误差滤波用来去除含有噪声的模糊图像，其目标是寻找未污染图像的一个估计，使均方误差最小：
e2=E{(f−f^)2}e^2 = E\{ (f - \hat{f})^2 \} e2=E{(f−f^)2}
最小均方差滤波由 Wiener [1942] 首先提出，是最早提出的线性图像复原方法，因此称为维纳滤波。

信噪比 SNR(f)=S(f)/N(f)SNR(f) = S(f)/N(f)SNR(f)=S(f)/N(f)，是信息承载信号功率（未退化的原图像）水平与噪声功率水平的测度。低噪声图像的信噪比较高，高噪声图像的信噪比较低。

将复原后的图像视为“信号”，而复原图像与原图像的差视为“噪声”，则可以将空间域的信噪比定义为：
SNR=∑x=0M−1∑y=0N−1f^(x,y)2/∑x=0M−1∑y=0N−1[f(x,y)−f^(x,y)]2SNR = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} \hat{f}(x,y)^2 / \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} [f(x,y) -\hat{f}(x,y)]^2 SNR=x=0∑M−1y=0∑N−1f^(x,y)2/x=0∑M−1y=0∑N−1[f(x,y)−f^(x,y)]2
最小均方误差滤波的传递函数描述为：
G(f)=H∗(f)S(f)∣H(f)∣2S(f)+N(f)=H∗(f)∣H(f)∣2S(f)+N(f)/S(f)G(f) = \frac{H^* (f) S(f)}{|H(f)|^2 S(f) + N(f)} = \frac{H^* (f)}{|H(f)|^2 S(f) + N(f)/S(f)} G(f)=∣H(f)∣2S(f)+N(f)H∗(f)S(f)=∣H(f)∣2S(f)+N(f)/S(f)H∗(f)
式中 G(f)、H(f)G(f)、H(f)G(f)、H(f) 是 g 和 h 在频率域的傅里叶变换，S(f)、N(f)S(f)、N(f)S(f)、N(f) 是信号 x(t)、噪声 n(t) 的功率谱。

用信噪比将上式进一步表示为：
G(f)=1H(f)∣H(f)∣2∣H(f)∣2+1/SNR(f)G(f) = \frac{1}{H(f)} \frac{|H(f)|^2}{|H(f)|^2 + 1/SNR(f)} G(f)=H(f)1∣H(f)∣2+1/SNR(f)∣H(f)∣2
当处理白噪声时∣N(u,v)∣2|N(u,v)|^2∣N(u,v)∣2 是一个常数，但未退化图像和噪声的功率谱通常未知或不能估计，则可用下式近似：
F^(u,v)=1H(u,v)∣H(u,v)∣2∣H(u,v)∣2+KG(u,v)\hat{F}(u,v)= \frac{1}{H(u,v)} \frac{|H(u,v)|^2}{|H(u,v)|^2 + K} G(u,v) F^(u,v)=H(u,v)1∣H(u,v)∣2+K∣H(u,v)∣2G(u,v)
K 是加到 ∣H(u,v)∣2|H(u,v)|^2∣H(u,v)∣2的所有项上的一个规定常数。

例程 9.21: 维纳滤波（Wiener filter）

    # 9.21: 退化图像的维纳滤波 (Wiener filter)def getMotionDsf(shape, angle, dist):xCenter = (shape[0] - 1) / 2yCenter = (shape[1] - 1) / 2sinVal = np.sin(angle * np.pi / 180)cosVal = np.cos(angle * np.pi / 180)PSF = np.zeros(shape)  # 点扩散函数for i in range(dist):  # 将对应角度上motion_dis个点置成1xOffset = round(sinVal * i)yOffset = round(cosVal * i)PSF[int(xCenter - xOffset), int(yCenter + yOffset)] = 1return PSF / PSF.sum()  # 归一化def makeBlurred(image, PSF, eps):  # 对图片进行运动模糊fftImg = np.fft.fft2(image)  # 进行二维数组的傅里叶变换fftPSF = np.fft.fft2(PSF) + epsfftBlur = np.fft.ifft2(fftImg * fftPSF)fftBlur = np.abs(np.fft.fftshift(fftBlur))return fftBlurdef inverseFilter(image, PSF, eps):  # 逆滤波fftImg = np.fft.fft2(image)fftPSF = np.fft.fft2(PSF) + eps  # 噪声功率，这是已知的，考虑epsilonimgInvFilter = np.fft.ifft2(fftImg / fftPSF)  # 计算F(u,v)的傅里叶反变换imgInvFilter = np.abs(np.fft.fftshift(imgInvFilter))return imgInvFilterdef wienerFilter(input, PSF, eps, K=0.01):  # 维纳滤波，K=0.01fftImg = np.fft.fft2(input)fftPSF = np.fft.fft2(PSF) + epsfftWiener = np.conj(fftPSF) / (np.abs(fftPSF)**2 + K)imgWienerFilter = np.fft.ifft2(fftImg * fftWiener)imgWienerFilter = np.abs(np.fft.fftshift(imgWienerFilter))return imgWienerFilter# 读取原始图像img = cv2.imread("../images/Fig0526a.tif", 0)  # flags=0 读取为灰度图像hImg, wImg = img.shape[:2]# 不含噪声的运动模糊PSF = getMotionDsf((hImg, wImg), 45, 100)  # 运动模糊函数imgBlurred = np.abs(makeBlurred(img, PSF, 1e-6))  # 生成不含噪声的运动模糊图像imgInvFilter = inverseFilter(imgBlurred, PSF, 1e-6)  # 逆滤波imgWienerFilter = wienerFilter(imgBlurred, PSF, 1e-6)  # 维纳滤波# 带有噪声的运动模糊scale = 0.05  # 噪声方差noisy = imgBlurred.std() * np.random.normal(loc=0.0, scale=scale, size=imgBlurred.shape)  # 添加高斯噪声imgBlurNoisy = imgBlurred + noisy  # 带有噪声的运动模糊imgNoisyInv = inverseFilter(imgBlurNoisy, PSF, scale)  # 对添加噪声的模糊图像进行逆滤波imgNoisyWiener = wienerFilter(imgBlurNoisy, PSF, scale)  # 对添加噪声的模糊图像进行维纳滤波plt.figure(figsize=(9, 7))plt.subplot(231), plt.title("blurred image"), plt.axis('off'), plt.imshow(imgBlurred, 'gray')plt.subplot(232), plt.title("inverse filter"), plt.axis('off'), plt.imshow(imgInvFilter, 'gray')plt.subplot(233), plt.title("Wiener filter"), plt.axis('off'), plt.imshow(imgWienerFilter, 'gray')plt.subplot(234), plt.title("blurred image with noisy"), plt.axis('off'), plt.imshow(imgBlurNoisy, 'gray')plt.subplot(235), plt.title("inverse filter"), plt.axis('off'), plt.imshow(imgNoisyInv, 'gray')plt.subplot(236), plt.title("Wiener filter"), plt.axis('off'), plt.imshow(imgNoisyWiener, 'gray')plt.tight_layout()plt.show()

程序说明：

对于不含噪声的运动模糊图像，在已知运动模糊退化模型和参数的前提下，使用逆滤波可以很好地复原退化图像，逆滤波的性能优于维纳滤波。但是，考虑实际退化图像往往含有一定水平的加性噪声，此时即使已知退化模型，逆滤波的后的噪声几乎掩盖了图像内容，而维纳滤波的结果则较好。

6.3 约束最小二乘方滤波（Constrained Least Squares Filtering）

维纳滤波建立在退化函数和信噪比已知的前提下，这在实践中并不容易满足。

约束最小二乘方滤波仅要求噪声方差和均值的知识或估计，这些参数通常可以由一幅给定的退化图像算出，因而具有更为广泛的应用。而且，维纳滤波是以最小化一个统计准则为基础，因此是平均意义上的最优，而约束最小二乘方滤波则对每幅图像都会产生最优估计。

与维纳滤波相比，最小二乘方滤波对于高噪声和中等噪声处理效果更好，对于低噪声二者处理结果基本相同。当手工选择参数以取得更好的视觉效果时，约束最小二乘方滤波的效果有可能比维纳滤波的效果更好。

约束最小二乘方滤波的核心是解决退化函数对噪声的敏感性问题。降低噪声敏感的一种方法是，以平滑度量的最佳复原为基础，如图像的二阶导数（拉普拉斯算子），其数学描述是求一个带约束条件的准则函数的最小值：
minC=∑x=0M−1∑y=0N−1[▽2f(x,y)]2s.t.:∣∣g−Hf^∣∣2=∣∣η∣∣2min \ C = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} [\triangledown ^2 f(x,y)]^2 \\ s.t.: ||g - H \hat{f}||^2 = ||\eta||^2 min C=x=0∑M−1y=0∑N−1[▽2f(x,y)]2s.t.:∣∣g−Hf^∣∣2=∣∣η∣∣2
其中，∣∣a∣∣2||a||^2∣∣a∣∣2是欧几里德范数，f^\hat{f}f^ 是未退化图像的估计，$ \triangledown ^2 $ 是拉普拉斯算子。求准则函数 C 的最小值，便得到最好的平滑效果即退化图像的最佳复原。

该约束最小二乘方问题的解是：
F^(u,v)=[H∗(u,v)∣H(u,v)∣2+γ∣P(u,v)∣2]G(u,v)\hat{F}(u,v) = [\frac{H^*(u,v)}{|H(u,v)|^2 + \gamma |P(u,v)|^2}] \ G(u,v) F^(u,v)=[∣H(u,v)∣2+γ∣P(u,v)∣2H∗(u,v)] G(u,v)
式中，γ\gammaγ 是一个必须满足约束条件的参数，P(u,v)P(u,v)P(u,v) 是函数 p(x,y)p(x,y)p(x,y) 的傅里叶变换：
p(x,y)=[0−10−14−10−10]p(x,y) = \begin{bmatrix} 0 & -1 & 0 \\ -1 & 4 & -1\\ 0 & -1 &0 \end{bmatrix} p(x,y)=⎣⎡0−10−14−10−10⎦⎤
注意函数 P(u,v)P(u,v)P(u,v) 和 H(u,v)H(u,v)H(u,v) 的大小必须相等。如果 H(u,v)H(u,v)H(u,v) 的大小为 M∗NM*NM∗N，则 p(x,y) 必须嵌入 M∗NM*NM∗N 零阵列的中心，并保持偶对称。

例程 9.21: 约束最小二乘方滤波

    # 9.21: 约束最小二乘方滤波def getMotionDsf(shape, angle, dist):xCenter = (shape[0] - 1) / 2yCenter = (shape[1] - 1) / 2sinVal = np.sin(angle * np.pi / 180)cosVal = np.cos(angle * np.pi / 180)PSF = np.zeros(shape)  # 点扩散函数for i in range(dist):  # 将对应角度上motion_dis个点置成1xOffset = round(sinVal * i)yOffset = round(cosVal * i)PSF[int(xCenter - xOffset), int(yCenter + yOffset)] = 1return PSF / PSF.sum()  # 归一化def makeBlurred(image, PSF, eps):  # 对图片进行运动模糊fftImg = np.fft.fft2(image)  # 进行二维数组的傅里叶变换fftPSF = np.fft.fft2(PSF) + epsfftBlur = np.fft.ifft2(fftImg * fftPSF)fftBlur = np.abs(np.fft.fftshift(fftBlur))return fftBlurdef wienerFilter(input, PSF, eps, K=0.01):  # 维纳滤波，K=0.01fftImg = np.fft.fft2(input)fftPSF = np.fft.fft2(PSF) + epsfftWiener = np.conj(fftPSF) / (np.abs(fftPSF)**2 + K)imgWienerFilter = np.fft.ifft2(fftImg * fftWiener)imgWienerFilter = np.abs(np.fft.fftshift(imgWienerFilter))return imgWienerFilterdef getPuv(image):h, w = image.shape[:2]hPad, wPad = h - 3, w - 3pxy = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]])pxyPad = np.pad(pxy, ((hPad//2, hPad - hPad//2), (wPad//2, wPad - wPad//2)), mode='constant')fftPuv = np.fft.fft2(pxyPad)return fftPuvdef leastSquareFilter(image, PSF, eps, gamma=0.01):  # 约束最小二乘方滤波fftImg = np.fft.fft2(image)fftPSF = np.fft.fft2(PSF)conj = fftPSF.conj()fftPuv = getPuv(image)# absConj = np.abs(fftPSF) ** 2Huv = conj / (np.abs(fftPSF)**2 + gamma * (np.abs(fftPuv)**2))ifftImg = np.fft.ifft2(fftImg * Huv)ifftShift = np.abs(np.fft.fftshift(ifftImg))imgLSFilter = np.uint8(cv2.normalize(np.abs(ifftShift), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]return imgLSFilter# # 读取原始图像img = cv2.imread("../images/Fig0526a.tif", 0)  # flags=0 读取为灰度图像hImg, wImg = img.shape[:2]# 带有噪声的运动模糊PSF = getMotionDsf((hImg, wImg), 45, 100)  # 运动模糊函数imgBlurred = np.abs(makeBlurred(img, PSF, 1e-6))  # 生成不含噪声的运动模糊图像scale = 0.01  # 噪声方差noisy = imgBlurred.std() * np.random.normal(loc=0.0, scale=scale, size=imgBlurred.shape)  # 添加高斯噪声imgBlurNoisy = imgBlurred + noisy  # 带有噪声的运动模糊imgWienerFilter = wienerFilter(imgBlurNoisy, PSF, scale, K=0.01)  # 对含有噪声的模糊图像进行维纳滤波imgLSFilter = leastSquareFilter(imgBlurNoisy, PSF, scale, gamma=0.01)  # 约束最小二乘方滤波plt.figure(figsize=(9, 7))plt.subplot(231), plt.title("blurred image (dev=0.01)"), plt.axis('off'), plt.imshow(imgBlurNoisy, 'gray')plt.subplot(232), plt.title("Wiener filter"), plt.axis('off'), plt.imshow(imgWienerFilter, 'gray')plt.subplot(233), plt.title("least square filter"), plt.axis('off'), plt.imshow(imgLSFilter, 'gray')scale = 0.1  # 噪声方差noisy = imgBlurred.std() * np.random.normal(loc=0.0, scale=scale, size=imgBlurred.shape)  # 添加高斯噪声imgBlurNoisy = imgBlurred + noisy  # 带有噪声的运动模糊imgWienerFilter = wienerFilter(imgBlurNoisy, PSF, scale, K=0.01)  # 维纳滤波imgLSFilter = leastSquareFilter(imgBlurNoisy, PSF, scale, gamma=0.1)  # 约束最小二乘方滤波plt.subplot(234), plt.title("blurred image (dev=0.1)"), plt.axis('off'), plt.imshow(imgBlurNoisy, 'gray')plt.subplot(235), plt.title("Wiener filter"), plt.axis('off'), plt.imshow(imgWienerFilter, 'gray')plt.subplot(236), plt.title("least square filter"), plt.axis('off'), plt.imshow(imgLSFilter, 'gray')plt.tight_layout()plt.show()

程序说明：

对于含有不同水平噪声的运动模糊图像，最小二乘方滤波的图像复原效果都较好，特别是对于高噪声和中等噪声的处理效果更好。

对于最小二乘方滤波，可以交互地或循环地调整参数 γ\gammaγ 以取得更好的复原效果，例如参数 γ\gammaγ 可以从较小的取值递增以获得最优估计。

6.4 几何均值滤波

几何均值滤波器是维纳滤波的推广，其传递函数由括号内幂次分别为 α\alphaα 和 1−α1-\alpha1−α 的两个表达式组成：
F^(u,v)=[H∗(u,v)∣H(u,v)∣2]α[H∗(u,v)∣H(u,v)∣2+β[Sη(u,v)/Sf(u,v)]]1−αG(u,v)\hat{F}(u,v) = \Bigg[ \frac{H^*(u,v)}{|H(u,v)|^2} \Bigg] ^{\alpha} \ \Bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \beta [S_{\eta} (u,v) /S_f(u,v)]}\Bigg]^{1-\alpha}\ G(u,v) F^(u,v)=[∣H(u,v)∣2H∗(u,v)]α [∣H(u,v)∣2+β[Sη(u,v)/Sf(u,v)]H∗(u,v)]1−α G(u,v)
式中，$ \alpha$ 和 β\betaβ 是非负的实常数。

当 α=1/2\alpha=1/2α=1/2 时几何均值滤波器是幂次相同的两个量的乘积，这就是几何均值的含义。

当 α=1/2,β=1\alpha=1 / 2, \beta=1α=1/2,β=1 时，称为频谱均衡滤波器。当 α=1\alpha=1α=1 时，简化为逆滤波器；当 α=0\alpha=0α=0 时，简化为带参数的维纳滤波器，并在 β=1\beta=1β=1 时成为标准维纳滤波器。

例程 9.22: 几何均值滤波

    # 9.22: 几何均值滤波器def getMotionDsf(shape, angle, dist):xCenter = (shape[0] - 1) / 2yCenter = (shape[1] - 1) / 2sinVal = np.sin(angle * np.pi / 180)cosVal = np.cos(angle * np.pi / 180)PSF = np.zeros(shape)  # 点扩散函数for i in range(dist):  # 将对应角度上motion_dis个点置成1xOffset = round(sinVal * i)yOffset = round(cosVal * i)PSF[int(xCenter - xOffset), int(yCenter + yOffset)] = 1return PSF / PSF.sum()  # 归一化def makeBlurred(image, PSF, eps):  # 对图片进行运动模糊fftImg = np.fft.fft2(image)  # 进行二维数组的傅里叶变换fftPSF = np.fft.fft2(PSF) + epsfftBlur = np.fft.ifft2(fftImg * fftPSF)fftBlur = np.abs(np.fft.fftshift(fftBlur))return fftBlurdef wienerFilter(input, PSF, eps, K=0.01):  # 维纳滤波，K=0.01fftImg = np.fft.fft2(input)fftPSF = np.fft.fft2(PSF) + epsfftWiener = np.conj(fftPSF) / (np.abs(fftPSF)**2 + K)imgWienerFilter = np.fft.ifft2(fftImg * fftWiener)imgWienerFilter = np.abs(np.fft.fftshift(imgWienerFilter))return imgWienerFilterdef geometricMeanFilter(image, PSF, eps, K=1, alpha=1, beta=1):  # 几何均值滤波器fftImg = np.fft.fft2(image)fftPSF = np.fft.fft2(PSF)conj = fftPSF.conj()squarePSF = (fftPSF * conj).realHuv = np.power(conj / (squarePSF), alpha) * np.power(conj / (squarePSF + beta * K), 1-alpha)ifftImg = np.fft.ifft2(fftImg * Huv)ifftShift = np.abs(np.fft.fftshift(ifftImg))imgGMFilter = np.uint8(cv2.normalize(np.abs(ifftShift), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]return imgGMFilter# # 读取原始图像img = cv2.imread("../images/Fig0526a.tif", 0)  # flags=0 读取为灰度图像hImg, wImg = img.shape[:2]# 带有噪声的运动模糊PSF = getMotionDsf((hImg, wImg), 45, 100)  # 运动模糊函数imgBlurred = np.abs(makeBlurred(img, PSF, 1e-6))  # 生成不含噪声的运动模糊图像scale = 0.01  # 噪声方差noisy = imgBlurred.std() * np.random.normal(loc=0.0, scale=scale, size=imgBlurred.shape)  # 添加高斯噪声imgBlurNoisy = imgBlurred + noisy  # 带有噪声的运动模糊imgWienerFilter = wienerFilter(imgBlurNoisy, PSF, scale, K=0.01)  # 对含有噪声的模糊图像进行维纳滤波imgGMFilter = geometricMeanFilter(imgBlurNoisy, PSF, scale, K=0.01, alpha=0.5, beta=1)  # 约束最小二乘方滤波plt.figure(figsize=(9, 7))plt.subplot(231), plt.title("blurred image (dev=0.01)"), plt.axis('off'), plt.imshow(imgBlurNoisy, 'gray')plt.subplot(232), plt.title("Wiener filter"), plt.axis('off'), plt.imshow(imgWienerFilter, 'gray')plt.subplot(233), plt.title("geometric mean filter"), plt.axis('off'), plt.imshow(imgGMFilter, 'gray')scale = 0.1  # 噪声方差noisy = imgBlurred.std() * np.random.normal(loc=0.0, scale=scale, size=imgBlurred.shape)  # 添加高斯噪声imgBlurNoisy = imgBlurred + noisy  # 带有噪声的运动模糊imgWienerFilter = wienerFilter(imgBlurNoisy, PSF, scale, K=0.01)  # 维纳滤波imgGMFilter = geometricMeanFilter(imgBlurNoisy, PSF, scale, K=0.01, alpha=0, beta=1)  # 约束最小二乘方滤波plt.subplot(234), plt.title("blurred image (dev=0.1)"), plt.axis('off'), plt.imshow(imgBlurNoisy, 'gray')plt.subplot(235), plt.title("Wiener filter"), plt.axis('off'), plt.imshow(imgWienerFilter, 'gray')plt.subplot(236), plt.title("geometric mean filter"), plt.axis('off'), plt.imshow(imgGMFilter, 'gray')plt.tight_layout()plt.show()

程序说明：

7. 投影重建图像

7.1 计算机断层成像（Computed tomography）

图像重建（Image Reconstruction）的基本思想，就是通过探测物体的投影数据，重建物体的实际内部构造。

断层成像就是要获得物体内部的截面图像。利用 X射线、超声波等射线穿透被遮挡物体的透视投影图，可以计算恢复物体的断层图，进而可以利用断层图或直接从二维透视投影图重建物体的形状和内部结构。其原理是射线在穿过不同组织时的吸收率不同，在成像面上得到不同的投射强度，由此反演求得内部分布的图像。

X 射线、CT 技术就是应用断层重建的医学诊断方法。计算机断层扫描（CT）已经发展成为一种非常成功和必不可少的医学诊断工具，被认为是 X 射线发现以来医学影像领域最伟大的发明。

X 射线计算机断层成像的基本原理是：使用 X 射线从多个不同方向和角度对物体进行扫描，通过反投影算法获取物体内部结构的切片，堆叠这些切片就可以得到人体的三维表示。

投影重建还应用于地矿探测，接收不同地层和矿体反射的超声波，按照超声波在媒质的透射率和反射规律，对透射投影图进行分析计算，就可以恢复重建地下的矿体形状。

7.2 投影和雷登变换（Radon transform）

雷登变换是三维重建的数学基础。

一条直线 y=ax+by=ax+by=ax+b 的法线束为：
xcosθ+ysinθ=ρx cos \theta + y sin \theta = \rho xcosθ+ysinθ=ρ
沿该直线对函数 f(x,y)f(x,y)f(x,y) 进行积分，这个积分结果就是雷登变换后的函数 Rf\mathcal{R}fRf 在这条直线上的值：
g(ρ,θ)=∫−∞+∞∫−∞+∞f(x,y)δ(xcosθ+ysinθ−ρ)dxdyg(\rho, \theta) = \int ^{+\infty}_{-\infty} \int ^{+\infty}_{-\infty} f(x,y) \ \delta(x cos \theta + y sin \theta - \rho) dxdy g(ρ,θ)=∫−∞+∞∫−∞+∞f(x,y) δ(xcosθ+ysinθ−ρ)dxdy
雷登变换的离散形式：
g(ρ,θ)=∑x=0M−1∑y=0N−1f(x,y)δ(xcosθ+ysinθ−ρ)g(\rho, \theta) = \sum ^{M-1}_{x=0} \sum^{N-1}_{y=0} f(x,y) \ \delta(x cos \theta + y sin \theta - \rho) g(ρ,θ)=x=0∑M−1y=0∑N−1f(x,y) δ(xcosθ+ysinθ−ρ)

雷登变换将二维空间 xyxyxy 映射到了另一个二维空间 αs\alpha sαs，称为直线空间。

雷登变换常被用于医学影像处理，利用反雷登变换可以进行三维图像重建。

例程 9.23：雷登变换正弦图

    # # 9.23: 雷登变换正弦图from scipy import ndimagedef discreteRadonTransform(image, steps):channels = image.shape[0]res = np.zeros((channels, channels), dtype=np.float32)for s in range(steps):rotation = ndimage.rotate(image, -s * 180/steps, reshape=False).astype(np.float32)res[:, s] = sum(rotation)return res# # 读取原始图像plt.figure(figsize=(9, 7))fileImg = ["Fig0534a.tif", "Fig0534b.tif", "Fig0534c.tif"]  # 原始图像 文件名for i in range(len(fileImg)):img = cv2.imread("../images/"+fileImg[i], flags=0)  # flags=0 读取为灰度图像imgRadon = discreteRadonTransform(img, img.shape[0])  # Radon 变换print(img.shape, imgRadon.shape)plt.subplot(2, len(fileImg), i + 1), plt.axis('off')  # 绘制原始图像plt.title("origin image"), plt.imshow(img, 'gray')plt.subplot(2, len(fileImg), i + len(fileImg) + 1), plt.axis('off')  # 绘制 sinogram 图plt.title("Radon transform"), plt.imshow(imgRadon, 'gray')plt.tight_layout()plt.show()

7.2 雷登变换反投影重建图像（Radon transform back projection）

空间点 X 通过摄像机 P 被作用到图像平面的图像点 m = PX ，获得采集到的图像数据，这种投影关系称为摄像机的正向投影（forward projection），简称投影。

反向投影是针对图像平面的基本几何元素而言的，图像平面点 m 的反投影是指在摄像机 P 的作用下具有像点 m 的所有空间点的集合。

反投影一个点形成部分图像的过程，是将直线 L(ρj,θk)L(\rho_j,\theta_k)L(ρj,θk) 复制到图像上，直线上每点的灰度值是 g(ρj,θk)g(\rho_j,\theta_k)g(ρj,θk)。遍历投影信号中的每个点，得到：
fθ(x,y)=g(ρ,θ)=g(xcosθ+ysinθ)f_{\theta}(x,y) = g(\rho,\theta) = g(x cos \theta + ysin \theta) fθ(x,y)=g(ρ,θ)=g(xcosθ+ysinθ)
对所有反投影图像积分，得到最终的图像：
f(x,y)=∫0πfθ(x,y)dθf(x,y) = \int_0^{\pi}f_{\theta}(x,y)d\theta f(x,y)=∫0πfθ(x,y)dθ
其离散形式为：

f(x,y)=∑θ=0πfθ(x,y)f(x,y) = \sum_{\theta=0}^{\pi}f_{\theta}(x,y) f(x,y)=θ=0∑πfθ(x,y)

反投影的图像有时称为层图，可以理解为投影图像的一个近似。

例程 9.24：雷登变换反投影重建图像

    # 9.24: 雷登变换反投影重建图像from scipy import ndimagedef discreteRadonTransform(image, steps):  # 离散雷登变换channels = image.shape[0]resRadon = np.zeros((channels, channels), dtype=np.float32)for s in range(steps):rotation = ndimage.rotate(image, -s * 180/steps, reshape=False).astype(np.float32)resRadon[:, s] = sum(rotation)return resRadondef inverseRadonTransform(image, steps):  # 雷登变换反投影channels = image.shape[0]res = np.zeros((steps, channels, channels))for s in range(steps):expandDims = np.expand_dims(image[:, s], axis=0)repeat = expandDims.repeat(channels, axis=0)res[s] = ndimage.rotate(repeat, s * 180/steps, reshape=False).astype(np.float32)invRadon = np.sum(res, axis=0)return invRadon# 读取原始图像img1 = cv2.imread("../images/Fig0534a.tif", 0)  # flags=0 读取为灰度图像img2 = cv2.imread("../images/Fig0534c.tif", 0)# 雷登变换imgRadon1 = discreteRadonTransform(img1, img1.shape[0])  # Radon 变换imgRadon2 = discreteRadonTransform(img2, img2.shape[0])# 雷登变换反投影imgInvRadon1 = inverseRadonTransform(imgRadon1, imgRadon1.shape[0])imgInvRadon2 = inverseRadonTransform(imgRadon2, imgRadon2.shape[0])plt.figure(figsize=(9, 7))plt.subplot(231), plt.axis('off'), plt.title("origin image"), plt.imshow(img1, 'gray')  # 绘制原始图像plt.subplot(232), plt.axis('off'), plt.title("Radon transform"), plt.imshow(imgRadon1, 'gray')  # 绘制 sinogram 图plt.subplot(233), plt.axis('off'), plt.title("Inv-Radon transform"), plt.imshow(imgInvRadon1, 'gray')plt.subplot(234), plt.axis('off'), plt.title("origin image"), plt.imshow(img2, 'gray')plt.subplot(235), plt.axis('off'), plt.title("Radon transform"), plt.imshow(imgRadon2, 'gray')plt.subplot(236), plt.axis('off'), plt.title("Inv-Radon transform"), plt.imshow(imgInvRadon2, 'gray')plt.tight_layout()plt.show()

7.3 滤波反投影重建图像（Filtered back projection）

直接由正弦图得到反投影图像，会存在严重的模糊，这是早期 CT 系统所存在的问题。

傅立叶中心切片定理表明，投影的一维傅立叶变换是得到投影区域的二维傅立叶变换的切片。滤波反投影重建算法在反投影前将每一个采集投影角度下的投影进行卷积处理，从而改善点扩散函数引起的形状伪影，有效地改善了重建的图像质量。

以 f(x,y)f(x,y)f(x,y) 表示需要重建的图像，F(u,v)F(u,v)F(u,v) 的二维傅里叶反变换是：
f(x,y)=∫−∞∞∫−∞∞F(u,v)ej2π(ux+vy)dudvf(x,y) = \int^{\infty}_{-\infty} \int^{\infty}_{-\infty} F(u,v) e^{j 2 \pi (ux+vy)} dudv f(x,y)=∫−∞∞∫−∞∞F(u,v)ej2π(ux+vy)dudv
利用傅里叶切片定理，可以得到：
f(x,y)=∫0π[∫−∞∞∣ω∣G(ω,θ)ej2πωρdω]dθf(x,y) = \int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} |\omega|G(\omega,\theta) e^{j 2 \pi \omega \rho} d\omega \big] d\theta f(x,y)=∫0π[∫−∞∞∣ω∣G(ω,θ)ej2πωρdω]dθ
括号 [] 内部是一个一维傅里叶反变换，可以认为这是一个一维滤波器的传递函数。由于 ∣ω∣|\omega|∣ω∣ 是一个不可积的斜坡函数（Slope function），可以通过对斜坡加窗进行限制，典型地如汉明窗（Hamming window）、韩窗（Hann window）。

该式也可以使用空间卷积来实现：
f(x,y)=∫0π[∫−∞∞∣ω∣G(ω,θ)ej2πωρdω]dθ=∫0π[∫−∞∞g(ρ,θ)s(xcosθ+ysinθ−ρ)dρ]dθf(x,y) = \int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} |\omega|G(\omega,\theta) e^{j 2 \pi \omega \rho} d\omega \big] d\theta \\ = \int^{\pi}_{0} \big[ \int^{\infty}_{-\infty} g(\rho,\theta) s(xcos \theta + ysin \theta - \rho) d\rho \big] d\theta \\ f(x,y)=∫0π[∫−∞∞∣ω∣G(ω,θ)ej2πωρdω]dθ=∫0π[∫−∞∞g(ρ,θ)s(xcosθ+ysinθ−ρ)dρ]dθ
这表明，将对应的投影 g(ρ,θ)g(\rho, \theta)g(ρ,θ) 与斜坡滤波器传递函数 s(ρ)s(\rho)s(ρ) 的傅里叶反变换进行卷积，可以得到角度 θ\thetaθ 的各个反投影，整个反投影图像可以通过对所有反投影图像积分得到。

例程 9.25：滤波反投影重建图像

SL 滤波器是使用 Sinc 函数对斜坡滤波器进行截断产生，其离散形式为：
hSL(nδ)=−2π2δ2(4∗n2−1),n=0,±1,±2,...h_{SL}(n \delta) = - \frac{2}{\pi ^2 \delta ^2 (4*n^2 -1)}, \ n=0,\pm 1,\pm 2,... hSL(nδ)=−π2δ2(4∗n2−1)2, n=0,±1,±2,...
利用投影值和 SL 滤波器进行卷积，然后再进行反投影，就可以实现图像重建。

   # 9.25: 滤波反投影重建图像from scipy import ndimagedef discreteRadonTransform(image, steps):  # 离散雷登变换channels = image.shape[0]resRadon = np.zeros((channels, channels), dtype=np.float32)for s in range(steps):rotation = ndimage.rotate(image, -s * 180/steps, reshape=False).astype(np.float32)resRadon[:, s] = sum(rotation)return resRadondef inverseRadonTransform(image, steps):  # 雷登变换反投影重建图像channels = image.shape[0]backproject = np.zeros((steps, channels, channels))for s in range(steps):expandDims = np.expand_dims(image[:, s], axis=0)repeat = expandDims.repeat(channels, axis=0)backproject[s] = ndimage.rotate(repeat, s * 180/steps, reshape=False).astype(np.float32)invRadon = np.sum(backproject, axis=0)return invRadondef SLFilter(N, d):  # SL 滤波器, Sinc 函数对斜坡滤波器进行截断filterSL = np.zeros((N,))for i in range(N):filterSL[i] = - 2 / (np.pi**2 * d**2 * (4 * (i-N/2)**2 - 1))return filterSLdef filterInvRadonTransform(image, steps):  # 滤波反投影重建图像channels = len(image[0])backproject = np.zeros((steps, channels, channels))  # 反投影filter = SLFilter(channels, 1)  # SL 滤波器for s in range(steps):value = image[:, s]  # 投影值valueFiltered = np.convolve(filter, value, "same")  # 投影值和 SL 滤波器进行卷积filterExpandDims = np.expand_dims(valueFiltered, axis=0)filterRepeat = filterExpandDims.repeat(channels, axis=0)backproject[s] = ndimage.rotate(filterRepeat, s * 180/steps, reshape=False).astype(np.float32)filterInvRadon = np.sum(backproject, axis=0)return filterInvRadon# 读取原始图像img1 = cv2.imread("../images/Fig0534a.tif", 0)  # flags=0 读取为灰度图像img2 = cv2.imread("../images/Fig0534c.tif", 0)# 雷登变换imgRadon1 = discreteRadonTransform(img1, img1.shape[0])  # Radon 变换imgRadon2 = discreteRadonTransform(img2, img2.shape[0])# 雷登变换反投影重建图像imgInvRadon1 = inverseRadonTransform(imgRadon1, imgRadon1.shape[0])imgInvRadon2 = inverseRadonTransform(imgRadon2, imgRadon2.shape[0])# 滤波反投影重建图像imgFilterInvRadon1 = filterInvRadonTransform(imgRadon1, imgRadon1.shape[0])imgFilterInvRadon2 = filterInvRadonTransform(imgRadon2, imgRadon2.shape[0])plt.figure(figsize=(9, 7))plt.subplot(231), plt.axis('off'), plt.title("origin image"), plt.imshow(img1, 'gray')  # 绘制原始图像plt.subplot(232), plt.axis('off'), plt.title("inverse Radon trans"), plt.imshow(imgInvRadon1, 'gray')plt.subplot(233), plt.axis('off'), plt.title("filter inv-Radon trans"), plt.imshow(imgFilterInvRadon1, 'gray')plt.subplot(234), plt.axis('off'), plt.title("origin image"), plt.imshow(img2, 'gray')plt.subplot(235), plt.axis('off'), plt.title("inverse Radon trans"), plt.imshow(imgInvRadon2, 'gray')plt.subplot(236), plt.axis('off'), plt.title("filter inv-Radon trans"), plt.imshow(imgFilterInvRadon2, 'gray')plt.tight_layout()plt.show()

注：本文中部分原始图片来自 Rafael C. Gonzalez “Digital Image Processing, 4th.Ed.”，特此致谢。

youcans 的 OpenCV 学习课原创作品
转载必须标注原文链接：https://blog.csdn.net/youcans/article/details/123192478
Copyright 2021 youcans, XUPT
Crated：2022-02-28

欢迎关注『youcans 的 OpenCV 学习课』系列，持续更新
欢迎关注『youcans 的 OpenCV 例程 200 篇』系列，持续更新中
youcans 的 OpenCV 学习课—1. 安装与环境配置
youcans 的 OpenCV 学习课—2. 图像读取与显示
youcans 的 OpenCV 学习课—3. 图像的创建与修改
youcans 的 OpenCV 学习课—4. 图像的叠加与混合
youcans 的 OpenCV 学习课—5. 图像的几何变换
youcans 的 OpenCV 学习课—6. 灰度变换与直方图处理
youcans 的 OpenCV 学习课—7. 空间域图像滤波
youcans 的 OpenCV 学习课—8. 频率域图像滤波（上）
youcans 的 OpenCV 学习课—9. 频率域图像滤波（下）
youcans 的 OpenCV 学习课—10. 图像复原与重建