数字图像处理(4版)——第 5 章——图像恢复和重建(上)(Rafael C.Gonzalez&Richard E. Woods)

与图像增强一样，恢复技术的主要目标是在某种预定义的意义上改善图像。尽管存在一些重叠邻域，但图像增强在很大程度上是一个主观过程，而图像恢复在很大程度上是一个客观过程。恢复试图利用对退化现象的先验知识来恢复已退化的图像。因此，恢复技术致力于对退化进行建模，并应用逆过程来恢复原图像。在本章中，我们将讨论适用于各种恢复情况的线性、空间不变的恢复模型。我们还将讨论从投影重建图像的基本技术及其在计算机分片扫描 (CT) 中的应用，CT 是图像处理最重要的商业应用之一，尤其是在医疗保健邻域。

5.1 一个图像退化/恢复过程模型(A model of the image degradation/restoration process)

在本章中，我们将图像退化建模为一个算子 ℋ ，它与一个加法噪声项一起对一幅输入图像 f (x ,y) 进行计算，生成一幅退化图像 g(x ,y) (见图 5.1 )。已知 g(x ,y)，某个关于 ℋ 的知识，和某个关于加法的噪声项 η(x ,y) ，则恢复的目标是获得原图像的估算值 $\hat{f}(x,y)$ 。

我们将在第 5.5 节中证明，如果 ℋ 是线性和位置不变的算子，则退化图像在空间域中表示为

(5-1) $g(x ,y) = (h {\bigstar} f )(x ,y) + \eta(x ,y)$

其中，h(x ,y) 是退化图像的空间表示。与在第 3 和第 4 章中一样，“★” 表示卷积(折积)。根据卷积定理可知，公式 (5-1) 在频率中的等价形式为

(5-2) $G(u ,v) = H(u ,v)F(u ,v) + N(u ,v)$

其中大写字母表示的项是公式 (5-1) 中相应项的 Fourier 变换。这两个公式是本章大部分恢复内容的基础。

在接下来的三节中，我们只讨论噪声引起的图像退化。从 5.5 节开始，我们将讨论在同时存在 ℋ 和 η 的情况下进行图像恢复的几种方法。

-------------------------------图 5.1：一个图像退化/恢复过程模型--------------

5.2 噪声模型(Noise models)

数字图像的主要噪声源出现在图像采集和/或传输过程中。图像传感器的性能受图像采集过程中各种环境因素以及传感元件本身质量的影响。例如，在使用 CCD 相机采集图像时，光照水平和传感器温度是影响最终图像噪声量的主要因素。图像在传输过程中损坏主要是由于传输信道中的干扰。例如，使用无线网络传输的图像可能会受到雷电或其他大气干扰的破坏。

5.2.1 噪声的空间和频率特性(Spatial and frequency properties of noise)

与我们的讨论相关的是定义噪声空间特性的参数，以及噪声是否与图像相关。频率特性是指噪声在Fourier(频域)中的频率成分，这在第4章中详细讨论过。例如，当噪声的Fourier谱为常数时，该噪声被称为白噪声。这个术语源自白光的物理特性，白光以相等的比例包含可见光谱中的所有频率。

除空间周期性噪声外，本章假设噪声与空间坐标无关，并且与图像本身不相关(即像素值与噪声成分值之间没有相关性)。虽然这些假设在某些应用中至少部分无效(例如，量子限制成像，例如X射线和核医学成像)，但处理空间相关和相关噪声的复杂性超出了我们的讨论范围。

5.2.2 一些重要的噪声概率密度函数(Some important noise probability density functions)

在接下来的讨论中，我们将关注图 5.1 中模型噪声分量中强度值的统计行为。这些可以视为随机变量，并以概率密度函数 (PDF) 为特征，正如前面简要提到的那样。图 5.1 中模型的噪声分量是一个与输入图像大小相同的图像 h(x, y)。我们通过生成一个数组来创建用于模拟的噪声图像，该数组的强度值是具有指定概率密度函数的随机数。这种方法适用于所有即将讨论的 PDF，但椒盐(或盐椒)噪声(salt-and-pepper noise)除外，其应用方式有所不同。以下是图像处理应用中最常见的噪声 PDF。

5.2.2.1 Gauss噪声 (Gaussian Noise)

由于Gauss噪声模型在空间和频域上都易于数学处理，因此在实际应用中非常常见。事实上，这种易处理性非常便捷，以至于Gauss模型经常被用于那些充其量只能勉强适用的场合。

一个Gauss随机变量 z 的 PDF 定义为下述我们熟悉的表达式：

(5-3) $\displaystyle p(z)= \frac{1}{\sigma\sqrt{2\pi}}e^{\displaystyle \frac{(z-\overline{z})^{2}}{2{\sigma}^{2}}} (-\infty < z <\infty)$

其中，z 表示密度， $\overline{z}$ 表 z 的均值(mean或average)，而 σ 是其准的偏差(或离差)(deviation)。图 5.2 (a) 表示这个函数的一段图像。我们知道，对于一个 Gauss随机变量，z 的值在范围 $\overline{z} \pm {\sigma}$ 的范围内的概率接近 0.68 ；z 的值在范围 $\overline{z} \pm 2{\sigma}$ 的范围内的概率接近 0.95 。

---------------------------图 5.2：一些重要的概率密度函数----------------------

5.2.2.2 Rayleigh噪声 (Rayleigh Noise)

Rayleigh 噪声的 PDF 定义为

(5-4) $\displaystyle p(z)=\left \{ \begin{array}{l} \displaystyle \frac{2}{b}(z-a)e^{\displaystyle -(z-a)^{2}/b}\hspace{0.2cm}(z \ge a) \\ \\ 0 \hspace{0.2cm}(z<a) \end{array} \right .$

这时，这个随机变量 z 的均值和方差(variance)(译注：样本值与均值之差的平方之和的均值)分别为

(5-5) $\overline{z}=a+\sqrt{{\pi}b/4}$

和

(5-6) $\displaystyle \sigma^{2}=\frac{b(4-\pi)}{4}$

图 5.2(b) 展示了Rayleigh密度图。请注意，图 5.2(b) 偏离原点的位置，以及密度图的基本形状向右倾斜。Rayleigh密度图对于建模偏斜(skewed)直方图的形状非常有用。

5.2.2.3 Erlang(Gamma)噪声 (Rayleigh Noise)

Erlang的PDF 定义为

(5-7) $\displaystyle p(z)=\left \{ \begin{array}{l} \displaystyle \frac{a^{b}z^{b-1}}{(b-1)!}e^{ \displaystyle -az}\hspace{0.2cm}(z \ge 0) \\ \\ 0 \hspace{0.2cm}(z<a) \end{array} \right .$

其中，参数使得 a > b ，b是正整数，且 “！” 号表示阶乘(factorial)。z的均值和方差分别为

(5-8) $\displaystyle \overline{z}=\frac{b}{a}$

和

(5-9) $\displaystyle \sigma^{2}=\frac{b}{a^{2}}$

图 5.2(c) 显示了该密度的图形。虽然公式 (5-9) 通常被称为 γ (gamma) 密度，但严格来说，只有当分母是γ 函数 (b) 时，这才是正确的。当分母如图所示时，该密度更适合称为Erlang密度。

5.2.2.4 指数噪声 (Exponential Noise)

指数噪声的 PDF 定义为

(5-10) $\displaystyle p(z)=\left \{ \begin{array}{l} \displaystyle ae^{ \displaystyle -az}\hspace{0.2cm}(z \ge 0) \\ \\ 0 \hspace{0.2cm}(z<a) \end{array} \right .$

其中，a > 0 。z的均值和方差分别为

(5-11) $\displaystyle \overline{z}=\frac{1}{a}$

和

(5-12) $\displaystyle \sigma^{2}=\frac{1}{a^{2}}$

请注意，此 PDF 是 Erlang 之PDF 的一个特例，其中 b = 1。图 5.2(d) 显示了指数密度函数的图。

5.2.2.5 均匀噪声 (Uniform Noise)

均匀噪声的PDF 定义为

(5-13) $\displaystyle p(z)=\left \{ \begin{array}{l} \displaystyle \frac{1}{b-a} \hspace{0.2cm}(a \le z \le b ) \\ \\ 0 \hspace{0.2cm}(\rm otherwise) \end{array} \right .$

z 的均值和方差分别为

(5-14) $\displaystyle \overline{z}=\frac{a+b}{2}$

和

(5-15) $\displaystyle \sigma^{2}=\frac{(b-a)^{2}}{12}$

图 5.2 (e) 显示了均匀密度的图。

5.2.2.6 椒盐噪声 (Salt-and-Pepper Noise)

若 k 表示一幅数字图像中用于表示强度值的位的数量，则那幅图像可能强度值的范围是

$[0 , 2^{k} - 1]$ ( 例如，8位图像的范围为 [0, 256] ) 。椒盐噪声的 PDF 是

(5-16) $\displaystyle p(z)=\left \{ \begin{array}{l} \displaystyle P_{s} \hspace{0.2cm}(z = 2^{k}-1) \\ \\ P_{p} \hspace{0.2cm}(z=0) \\ \\ 1-(P_{s}+P_{p}) \hspace{0.2cm}(z=V) \end{array} \right .$

其中，V 是一个介于 $0 < V < 2^{k} - 1$ 范围内的一个整数值。

令 η( x , y ) 表示一幅椒盐噪声图像，其强度值满足公式 (5-16)。已知一幅与 η( x , y ) 相同大小的图像 f ( x , y ) ，我们用椒盐噪声破坏它，将 f 中所有在 h 中出现 0 的位置赋值为 0 。类似地，我们将 $2^{k}-1$ 的值赋给 f 中所有在 h 中出现该值的位置。最后，我们保持 f 中所有在 h 中出现 V 的位置不变。

若 $P_{s}$ 和 $P_{p}$ 均不为零，且特别是若它们相等，则满足公式 (5-16) 的噪声将是白色 $(2^{k}-1)$ 或黑色 (0)，并且会像随机分布在图像上的盐粒和胡椒颗粒(pepper granules)一样；因此这种噪声类型因此而得名。文献中还会用到其他名称，例如双极脉冲噪声(bipolar impulse noise) (如果 $P_{s}$ 或 $P_{p}$ 之一为 0，则为单极)，数据丢失噪声(data-drop-out noise)和尖峰噪声(spike noise)。我们在文中交替使用脉冲噪声和椒盐噪声这两个术语。

一个像素被椒噪声或盐噪声干扰的概率是 $P = P_{s} + P_{p}$ 。常用术语 P 表示噪声密度(noise density)。例如，若 $P_{s} = 0.02$ 而 $P_{p} = 0.01$ ，则 P = 0.03 ，我们称图像中接近 2% 的像素被盐噪声干扰，接近 1% 的像素被椒噪声干扰，且噪声密度是 3% ，意味着接近 3% 的图像像素被椒盐噪声干扰。

正如我们所见，尽管椒盐噪声由每一个概率指下，而不由均值和方差指定，但出于完备性，我们在此包括后者。椒盐噪声的均值为

(5-17) $\overline{z} = (0)P_{p} + K(1 - P_{s} - P_{p} ) + (2^{k} - 1) P_{s}$

而方差为

(5-18) $\displaystyle {\sigma}^{2} = (0-\overline{z})^{2}P_{p} + (K-\overline{z})^{2}(1 - P_{s} - P_{p} ) + (2^{k}-1)^{2} P_{s}$

其中，我们在两个方程中都明确包含了 0 这个值，以表明椒噪声的值假定为零。

总体而言，上述概率密度函数 (PDF) 为建模实际中遇到的各种噪声破坏情况提供了实用工具。例如，由于电子电路噪声以及光照不足和/或高温引起的传感器噪声等因素，图像中会出现Gauss噪声。Rayleigh密度有助于表征距离成像中的噪声现象。指数密度和γ密度在激光成像中得到应用。脉冲噪声常见于成像过程中发生快速瞬变(例如开关故障)的情况。均匀密度可能是对实际情况描述能力最差的。然而，均匀密度作为众多在模拟中广泛使用的随机数生成器的基础，非常有用(Gonzalez ，Woods 和 Eddins [2009])。

例子 5.1 ：噪声图像及其直方图。

图 5.3 展示了一个用于说明刚才讨论的噪声模型的测试图案。这是一个合适的图案，因为它由简单、恒定的区域组成，这些区域仅以三个增量跨越从黑色到接近白色的灰阶。这有助于对添加到图像中的各种噪声成分的特征进行视觉分析。

-------------------图5.3：用于说明图 5.2 中的 PDF 特征的测试模式。------------

图 5.4 展示了添加图 5.2 中六种噪声后的测试图。每幅图像下方是直接根据该图像计算得出的直方图。每种情况下，噪声参数的选择都应确保测试图上对应三个强度级别的直方图开始融合。这使得噪声清晰可见，同时又不会遮挡底层图像的基本结构。

将图 5.4 中的直方图与图 5.2 中的概率密度函数 (PDF) 进行比较，我们发现它们之间存在密切的对应关系。椒盐噪声示例的直方图中没有包含 V 的特定峰值，因为正如你所记得的那样，V 仅在创建噪声图像时使用，以使原图像中的值保持不变。当然，除了椒盐噪声峰值之外，图像中其他强度的值也都有峰值。除了整体强度略有不同外，图 5.4 中的前五幅图像很难在视觉上区分，即使它们的直方图明显不同。图 5.4(i) 中椒盐噪声图像的外观是唯一能够直观地表明导致质量下降的噪声类型的图像。

-----------------------------图5.4：图 5.3 中，在图像中添加Gauss、Rayleigh和Erlanga噪声后得到的图像和直方图。------------------------------------

-----------------------------图 5.4 (续)：图 5.3 中，对图像添加了指数噪声、均匀噪声和椒盐噪声后的图像和直方图。在椒盐直方图中，原点(零强度)和刻度远端的峰值略有偏移，因此不会与页面背景融合。------------------------------------

5.2.2.7 周期噪声 (Periodic-Noise)

图像中的周期性噪声通常源于图像采集过程中的电气或机电干扰。这是本章中我们将要讨论的唯一一类空间相关噪声。正如我们将在 5.4 节中讨论的那样，周期性噪声可以通过频域滤波显著降低。例如，考虑图 5.5(a) 中的图像。该图像被加性(空间)正弦噪声破坏。纯正弦波的Fourier变换是一对位于正弦波共轭频率上的共轭脉冲(见表 4.4)(注：请注意不要将频域中的术语“脉冲”与前面讨论过的脉冲噪声中在空间域中使用的相同术语混淆)。因此，如果正弦波在空间域中的振幅足够强，我们期望在图像的频谱中看到图像中每一个正弦波对应的一对脉冲。如图 5.5(b) 所示，情况确实如此。消除或降低频域中的这些脉冲将消除或降低空间域中的正弦噪声。我们将在第 5.4 节中更详细地讨论这个以及其他周期性噪声的例子。

-------------------------------图 5.5：(a)被加性正弦噪声损坏的图像。(b)频谱显示由正弦波引起的两个共轭脉冲。(原图由 NASA 提供。)--------------------------

5.2.2.8 噪声参数估算 (Estimating Noise Parameters)

周期性噪声的参数通常通过检查Fourier谱来估算。周期性噪声往往会产生频率尖峰，即使通过目视分析也常常能够检测到。另一种方法是尝试直接从图像中推断噪声成分的周期性，但这仅适用于较为简单的情况。在噪声尖峰异常明显或已知干扰频率成分大致位置的情况下，可以进行自动分析(参见第 5.4 节)。

噪声概率密度函数 (PDF) 的参数可能部分来自传感器规格，但通常需要针对特定的成像方案进行估算。如果成像系统可用，研究系统噪声特性的一个简单方法是捕获一组“平面”图像。例如，对于光学传感器，这就像对均匀照明的纯灰色板进行成像一样简单。生成的图像通常可以很好地指示系统噪声。

当仅有传感器已生成的图像可用时，通常可以从背景强度相对恒定的小块图像中估算概率密度函数 (PDF) 的参数。例如，图 5.6 所示的垂直条带是从图 5.4 中的Gauss、Rayleigh和均匀图像中裁剪出来的。图中所示的直方图是使用这些小条带的图像数据计算得出的。图 5.4 中与图 5.6 中的直方图相对应的直方图位于图 5.4(d)、(e) 和 (k) 中三个直方图的中间。我们发现，这些直方图的形状与图 5.6 中对应直方图的形状非常接近。它们的高度由于缩放而有所不同，但形状明显相似。

-----------------------图 5.6：使用图 5.4 中 (a) Gauss图像、(b) Rayleigh图像和 (c) 均匀噪声图像中的小条带(如图所示)计算出的直方图。---------------------------

图像条带数据最简单的用途是计算强度等级的均值和方差。考虑一个图像条带(子图像)，记为 S，令 $p_{S}(z_{i})$ (i = 0, 1, 2,…, L – 1)表示 S 中像素强度的概率估算值(归一化直方图值)，其中 L 是整幅图像中可能强度的数量(例如，对于 8 位图像，L 为 256)。如公式 (2-69) 和公式 (2-70) 所示，我们估算 S 中像素值的均值和方差，如下所示：

(5-19) $\displaystyle \overline{z}= \sum_{i=0}^{L-1} z_{i} p_{S}(z_{i} )$

和

(5-20) $\displaystyle \sigma^{2} = \sum_{i=0}^{L-1}(z_{i}-\overline{z})^{2} p_{S}(z_{i} )$

直方图的形状标识了最接近的概率密度函数(PDF)匹配。如果形状近似于Gauss分布，那么我们只需要均值和方差，因为Gauss概率密度函数(PDF)完全由这两个参数指定。对于前面讨论过的其他形状，我们使用均值和方差来求解参数 a 和 b 。脉冲噪声的处理方式有所不同，因为需要估算的是白色和黑色像素出现的实际概率。获得此估算值需要黑色和白色像素都可见，因此图像中需要一个中等灰色、相对恒定的区域，以便能够计算出有意义的噪声直方图。对应于黑色和白色像素的峰值高度是公式 (5-16) 中 $P_{a}$ 和 $P_{b}$ 的估算值。

5.3 仅存在噪声的情况下的图像恢复——空间滤波(Restoration in the presence of noise only—spatial filtering)

当一幅图像的质量下降仅因加性噪声时，公式 (5-1) 和 (5-2) 就成了

(5-21) $g(x, y) = f (x, y) + \eta(x, y)$

和

(5-22) $G(u, v) = F(u, v) + N(u, v)$

噪声项通常是未知的，因此从 g(x, y)[ G(u, v)] 中减去这些噪声项来获得 f (x, y) [F(u, v)] 通常不可行。对于周期性噪声，有时可以根据 G(u, v) 的谱来估算 N(u, v)，如 5.2 节所述。在这种情况下，可以从 G(u, v) 中减去 N(u, v) 来获得原图像的估算值，但这种知识只是例外，而不是常态。

当仅存在加性随机噪声时，空间滤波是估算 f (x, y)(即去噪图像 g(x, y))的首选方法。空间滤波已在第 3 章详细讨论过。除了特定滤波器执行的计算性质不同外，下文所有滤波器的实现机制与第 3.4 至 3.7 节中讨论的完全相同。

5.3.1 均值滤波器(Mean Filters)

在本节中，我们简要讨论第 3.5 节中介绍的空间滤波器的降噪能力，并开发几个其他滤波器，它们的性能在许多情况下优于该节中讨论的滤波器。

5.3.1.1 算术均值滤波器 (Arithmetic Mean Filter)

算术均值滤波器是均值滤波器中最简单的一种(算术均值滤波器与我们在第三章讨论的箱式滤波器相同)。设 $S_{xy}$ 表示以点(x, y)为中心的m × n大小的矩形子图像窗口(邻域)中的坐标集。算术均值滤波器计算损坏图像 g (x, y) 在 $S_{xy}$ 所定义区域中的均值。恢复图像 $\hat{f}$ 在点(x, y)处的值是使用 $S_{xy}$ 所定义区域中的像素计算的算术均值。换言之，

(5-23) $\displaystyle \hat{f}(x, y) = \frac{1}{mn}\sum_{(r,c)\in S_{xy}}g(r,c)$

其中，如在公式 (2-43)中一样，r和c 分别是包含于邻域 $S_{xy}$ 中的像素的行坐标和列坐标。此操作可以使用大小为 m × n 的空间核来实现，其中所有系数的值均为 1/mn 。均值滤波器可以平滑图像中的局部变化，并通过模糊处理来降低噪声。

5.3.1.2 几何均值滤波器 (Geometric Mean Filter)

使用几何平均滤波器恢复的图像由表达式

(5-24) $\displaystyle \hat{f}(x, y) = \left [ \prod_{\displaystyle (r,c){\in}S_{xy}}g(r,c) \right ]^{ \displaystyle \frac{1}{mn}}$

给出，其中 Π 表示乘法。在这里，每一个恢复的像素都由子图像区域中所有像素自乘 1/(mn) 次方。如下面的示例 5.2 所示，几何均值滤波器的平滑效果与算术平均滤波器相当，但在此过程中丢失的图像细节往往较少。

5.3.1.3 调和均值滤波器 (Harmonic Mean Filter)

调和均值滤波运算由表达式

(5-25) $\displaystyle \hat{f}(x, y) = \frac{mn}{\displaystyle \sum_{(r,c){\in}S_{xy}}\frac{1}{g(r,c)}}$

给出。调和均值滤波器对盐噪声效果良好，但对椒噪声无效。它对其他类型的噪声(例如Gauss噪声)也表现出色。

5.3.1.4 反调和均值滤波器 (Contraharmonic Mean Filter)

反调和均值滤波器基于表达式

(5-26) $\displaystyle \hat{f}(x, y) = \frac{\displaystyle \sum_{(r,c){\in}S_{xy}}(g(r,c))^{Q+1}}{\displaystyle \sum_{(r,c){\in}S_{xy}}(g(r,c))^{Q}}$

恢复图像，其中，Q 称为滤波器的阶。该滤波器非常适合降低或几乎消除椒盐噪声的影响。当 Q 值为正时，该滤波器可消除椒盐噪声。当 Q 值为负时，该滤波器可消除盐噪声。它不能同时实现两者。请注意，如果 Q = 0，则反调和滤波器简化为算术均值滤波器；如果 Q = -1，则简化为谐波均值滤波器。

例子 5.2 ：使用空间均值滤波器进行图像降噪。

图 5.7(a) 显示了一张 8 位的电路板 X 射线图像；图 5.7(b) 显示了同一图像，但图像中含有均值为零、方差为 400 的加性Gauss噪声。对于此类图像而言，噪声水平相当高。图 5.7(c) 和 (d) 分别显示了使用大小为 3×3 的算术平均滤波器和大小相同的几何均值滤波器对噪声图像进行滤波的结果。虽然这两个滤波器都能有效地衰减噪声的影响，但几何均值滤波器对图像的模糊程度不如算术滤波器。例如，图 5.7(d) 中图像顶部的连接器指状物比 (c) 中更清晰。图像的其他部分也是如此。

------------------------------------图 5.7：(a)电路板的 X 射线图像。(b) 受加性高斯噪声影响的图像。(c) 使用大小为 3 × 3 的算术均值滤波器滤波的结果。(d) 使用相同大小的几何均值滤波器滤波的结果。(原图由 Joseph E. Pascente 先生，Lixi 公司提供)------------------------------------------------------------

图 5.8(a) 显示了相同的电路图像，但现在受到概率为 0.1 的胡椒噪声破坏。类似地，图 5.8(b) 显示了以相同概率受到盐噪声破坏的图像。图 5.8(c) 显示了使用 Q = 1.5 的反调和均值滤波器对图 5.8(a) 进行滤波的结果，图 5.8(d) 显示了使用 Q = −1.5 的反调和均值滤波器对图 5.8(b) 进行滤波的结果。这两个滤波器都很好地降低了噪声的影响。正阶滤波器在清洁背景方面做得更好，但代价是稍微稀释和模糊了暗区。负阶滤波器则相反。

------------------------------------图 5.8：(a)被概率为 0.1 的椒噪声干扰的图像。(b) 被概率为 0.1 的盐噪声干扰的图像。(c) 使用 3 × 3 反调和均值滤波器 Q = 1.5 对 (a) 进行滤波的结果(d) 使用 Q = −1.5 滤波(b) 的结果-----------------------

一般来说，算术均值滤波器和几何均值滤波器(尤其是后者)非常适合Gauss噪声或均匀噪声等随机噪声。反调和滤波器非常适合脉冲噪声，但它的缺点是必须知道噪声是暗噪声还是亮噪声才能为 Q 选择合适的符号。选择错误的 Q 符号可能会造成灾难性的后果，如图 5.9 所示。以下章节中讨论的一些滤波器可以弥补这一缺陷。

---------------------------------图 5.9：在反调和滤波的过程中选择了错误的符号的结果。(a) 用大小为 3 × 3 和 Q = −1.5 滤波图5.8(a)的结果。(b) 用Q = 1.5 滤波图5.8(b)的结果--------------------------------------------------------------

5.3.2 顺序统计滤波器 (Order-Statistic Filters)

我们在第 3.6 节中介绍了顺序统计滤波器。现在我们扩展该节的讨论，并介绍一些额外的顺序统计滤波器。如第 3.6 节所述，顺序统计滤波器是一种空间滤波器，其响应基于对滤波器所包含的邻域内像素值进行排序(排序)。排序结果决定了滤波器的响应。

5.3.2.1 中值滤波器 (Median Filter)

图像处理中最著名的顺序统计滤波器是中值滤波器，顾名思义，它用像素预定义邻域中强度级别的中值替换像素的值：

(5-27) $\displaystyle \hat{f}(x,y)=\large \mathop{ \text{median}}_{\displaystyle (r,c){\in}S_{xy}}\{ g(r,c) \}$

其中，如前一样， $S_{xy}$ 是中心位于( x , y )子图像(邻域)。位于 ( x , y ) 处的像素值包含于中值的计算中。中值滤波器非常流行，因为对于某些类型的随机噪声，它们具有出色的降噪能力，并且比类似大小的线性平滑滤波器的模糊程度要小得多。中值滤波器在双极和单极脉冲噪声中都特别有效，如下例 5.3 所示。中值的计算和该滤波器的实现将在第 3.6 节中讨论。

5.3.2.2 最大值和最小值滤波器 (Max and Min Filter)

虽然中值滤波器是迄今为止图像处理中最常用的顺序统计滤波器，但它绝不是唯一的。中值代表一组有序数字的第 50 个百分位数，但你会从基础统计学中回想起，排序本身就有很多其他的可能性。例如，使用第 100 个百分位数可以得到所谓的最大滤波器，其公式如下

(5-28) $\displaystyle \hat{f}(x,y)=\large \mathop{ \text{max}}_{\displaystyle (r,c){\in}S_{xy}}\{ g(r,c) \}$

此滤波器可用于查找图像中最亮点或去除与亮区相邻的暗区。此外，由于椒噪声值非常低，因此，此滤波器会通过子图像区域 $S_{xy}$ 中的最大值选择过程来降低椒噪声。

第 0 个百分位数滤波器是最小过滤器：

(5-29) $\displaystyle \hat{f}(x,y)=\large \mathop{ \text{min}}_{\displaystyle (r,c){\in}S_{xy}}\{ g(r,c) \}$

此滤波器可用于查找图像中最暗的点，或腐蚀与暗区相邻的亮区。此外，它还能通过最小值运算降低盐噪声。

5.3.2.3 中点滤波器 (Midpoint Filter)

中点滤波器计算滤波器所涵盖区域中最大值和最小值之间的中点：

(5-30) $\displaystyle \hat{f}(x,y)=\frac{1}{2}\left [\large \mathop{ \text{max}}_{\displaystyle (r,c){\in}S_{xy}}\{ g(r,c) \}+\large \mathop{ \text{min}}_{\displaystyle (r,c){\in}S_{xy}}\{ g(r,c) \} \right ]$

请注意，此滤波器结合了顺序统计和平均算法。它最适合随机分布的噪声，例如Gauss噪声或均匀噪声。

5.3.2.4 α 修剪滤波器 (Alpha-Trimmed Mean Filter)

假设我们删除邻域 $S_{xy}$ 中 g(r,c) 的 d/2 个最低强度值和 d/2 个最高强度值。令 $g_{R}(r,c)$ 表示 $S_{xy}$ 中剩余的 mn - d 个像素。由这些剩余像素求平均值而形成的滤波器称为 α修剪均值滤波器。该滤波器的形式为

(5-31) $\displaystyle \hat{f}(x, y) = \frac{1}{mn-d}\sum_{(r,c){\in}S_{xy}}g_{R}(r,c)$

其中 d 的值范围为 0 到 mn - 1。当 d = 0 时，α修剪滤波器简化为前面讨论过的算术平均值滤波器。如果 d = mn - 1，滤波器则变为中值滤波器。对于 d 的其他值，α修剪滤波器在涉及多种噪声的情况中很有用，例如椒盐噪声和Gauss噪声的组合。

例子 5.3：使用顺序统计滤波器进行图像去噪。

图 5.10(a) 显示了受椒盐噪声污染的电路板图像，噪声概率为 $P_{s} = P_{p} = 0.1$ 。图 5.10(b) 显示了使用大小为 3 × 3 的滤波器进行中值滤波的结果。与图 5.10(a) 相比，效果有了显著改善，但仍然可以看到一些噪声点。对图 5.10(b) 中的图像进行第二次中值滤波器滤波后，大部分噪声点被去除，只留下几个几乎看不见的噪声点。这些噪声点经过第三次滤波后被去除。这些结果很好地体现了中值滤波在处理脉冲状加性噪声方面的强大功能。谨记，反复进行中值滤波器滤波会使图像模糊，因此最好尽可能减少滤波次数。

---------------------------------图 5.10：(a) 受椒盐噪声破坏的图像，概率为 $P_{s} = P_{p} = 0.1$ (b) 使用大小为3 × 3 的中值滤波器处理一次的结果。(c) 使用此滤波器处理 (b) 的结果。(d) 使用相同滤波器处理 (c) 的结果。-------------------------------

图 5.11(a) 展示了将最大值滤波器应用于图 5.8(a) 中的椒噪声图像的结果。该滤波器在去除椒噪声方面做得不错，但我们注意到，它也去除了(设置为亮强度级别)暗物体边界上的一些暗像素。图 5.11(b) 展示了将最小值滤波器应用于图 5.8(b) 中的图像的结果。在这种情况下，最小值滤波器在去除噪声方面比最大值滤波器做得更好，但它去除了亮物体边界周围的一些白点。这使得亮物体更小，而一些暗物体(例如图像顶部的连接器手指)更大，因为这些物体周围的白点被设置为暗强度级别。

--------------------图 5.11：使用大小为 3 × 3 的最大滤波器对图 5.8(a) 进行滤波。(b) 使用大小相同的最小滤波器对图 5.8(b) 进行滤波的结果。 ---------------------

接下来展示的是 α修剪滤波器。图 5.12(a) 显示的电路板图像这次被方差为 800、均值为零的加性均匀噪声损坏。这是一个严重的噪声损坏，而进一步添加 $P_{s} = P_{p} = 0.1$ 的椒盐噪声会使情况更加恶化，如图 5.12(b) 所示。此图像中的高噪声水平需要使用更大的滤波器。图 5.12( c ) 至 ( f ) 分别显示了使用算术平均值、几何平均值、中值和 α修剪平均值(d = 6)滤波器(大小为 5 × 5)获得的结果。正如预期的那样，算术平均值和几何平均值滤波器(尤其是后者)由于存在脉冲噪声而表现不佳。中值滤波器和 α修剪滤波器的表现要好得多，其中 α 修剪滤波器的降噪效果略好。例如，在图 5.12(f) 中，左上角第四个连接点在 α 修剪后的结果中略微平滑一些。这并不意外，因为当 d 值较大时，α 修剪后的滤波器的性能接近中值滤波器，但仍保留了一些平滑能力。

--------------------图 5.12：(a) 图像被加性均匀噪声损坏。(b) 图像还被加性椒盐噪声损坏。(c)-( f ) 图像 (b) 经5 × 5 滤波：(c) 算术均值滤波器；(d) 几何均值滤波器；(e) 中值滤波器；( f ) α 修剪均值滤波器，d = 6。 --------------------------

5.3.3 自适应滤波器 (Adaptive Filters)

一旦选定，以上讨论的滤波器就会应用于图像，而无需考虑图像特征在不同点之间的变化。在本节中，我们将研究两个自适应滤波器，它们的行为会根据滤波区域内图像的统计特征而变化，该滤波区域由 m × n 矩形邻域 $S_{xy}$ 定义。正如下文讨论所示，自适应滤波器的性能优于以上讨论的滤波器。提升滤波能力的代价是滤波器复杂度的增加。请记住，我们仍然处理的是退化图像等于原始图像加上噪声的情况。目前尚未考虑其他类型的退化(degradations)。

5.3.3.1 自适应、局部降澡滤波器 (Adaptive, Local Noise Reduction Filter)

随机变量最简单的统计度量是其均值(mean)和方差(variance)。这些参数是自适应滤波器的合理基础，因为它们与图像的外观密切相关。均值衡量的是计算均值的区域内的平均强度，而方差衡量的是该区域内的图像对比度。

我们的滤波器是基于一个以 (x ,y) 为中心的邻域 $S_{xy}$ 的运算。滤波器在 (x ,y) 点处响应基于下述几个参量(quantities)：噪声图像在 (x ,y) 处的值 g(x ,y)；澡声之方差 $\sigma_{\eta}^{2}$ ； $S_{xy}$ 中的局部平均像素强度 $\overline{z}_{S_{xy}}$ ；以及 $S_{xy}$ 中像素强度的局部方差 $\sigma_{S_{xy}}^{2}$ 。

我们希望滤波器具有以下行为：

(1) 若 $\sigma_{\eta}^{2}$ 为零，滤波器应当仅返回 (x ,y) 处的像素值。这是关系，在零噪声的情况下 g 就应当等于 f 在 (x ,y) 处的值。

(2) 若局部方差 $\sigma_{S_{xy}}^{2}$ 相对于 $\sigma_{\eta}^{2}$ 较高，则滤波器应当返回一个近似于g 在 (x ,y) 处的像素值。一个高局部方差常与边相关，而这些边的值应当保留。

(3) 如果两个方差相等，我们希望滤波器返回 $S_{xy}$ 中像素的算术平均值。当局部区域与整体图像具有相同的属性时，就会发生这种情况，并且需要通过均值化来实现局部降噪。

基于这些假设的自适应表达式可以记为

(5-32) $\displaystyle \hat{f}(x,y)=g(x,y)-\frac{\sigma_{\eta}^{2}}{\sigma_{S_{xy}}^{2}}\left [ g(x,y)-\overline{z}_{S_{xy}} \right ]$

需要知道的一个唯一的先验参数是干扰图像 f (x ,y ) 的噪声方差 $\sigma_{\eta}^{2}$ 。这是一个常量，其可使用公式 (3-26) 根据样本噪声图像估算。其它参数可根据公式 (3-27) 和 (3-28) 根据邻域 $S_{xy}$ 中像素进行计算。

公式 (5-32) 假设这两个方差之比不超过 1 ，这意味着 $\displaystyle {\sigma}_{\eta}^{2} < {\sigma}_{S_{xy}} ^{2}$ 。我们模型中的噪声是加性的且位置独立的，因此这是一个合理的假设，因为 $S_{xy}$ 是 g (x ,y) 的一个子集。然而，我们缺乏对 ${\sigma}_{\eta}^{2}$ 的确切了解。因此，在实践中违背这个条件是可能的。为此，应该在公式 (5-32) 的实现中构建一个测试，以便当条件 $\displaystyle {\sigma}_{\eta}^{2} > {\sigma}_{S_{xy}} ^{2}$ 满足时，将比率设置为 1 。这使得滤波器是非线性的。然而，由于缺乏对图像噪声方差的了解，这种方法可能会避免出现无意义的结果(例如，取决于 $\overline{z}_{S_{xy}}$ 值的负强度水平)。另一种方法是允许出现负值，然后在最后重新调整强度值。结果将导致图像动态范围的损失。

例子 5.4：使用自适应局部降噪滤波进行图像降噪。

图 5.13(a) 显示了电路板图像，这次图像受到了均值为零、方差为 1000 的加性Gauss噪声的干扰。噪声干扰程度相当严重，但它却是一个理想的测试平台，可以用来比较滤波器的相对性能。图 5.13(b) 是使用大小为 7 × 7 的算术平均值滤波器处理噪声图像的结果。噪声被平滑了，但代价是图像明显模糊。图 5.13(c) 也存在类似的情况，它显示了使用大小同样为 7 × 7 的几何平均值滤波器处理噪声图像的结果。这两幅滤波图像之间的差异类似于我们在例 5.2 中讨论的差异，只是模糊程度不同。

图 5.13(d) 显示了使用公式 (5-32) 中 $\sigma_{\eta}^{2}=1000$ 的自适应滤波器的结果。与之前的两个滤波器相比，此结果的改进非常显著。在整体降噪方面，自适应滤波器取得了与算术平均滤波器和几何平均滤波器类似的效果。然而，经过自适应滤波器滤波的图像更加清晰。例如，图 5.13(d) 中图像顶部的连接器指状部分明显更加清晰。其他特征，例如孔洞和图像左下方暗色部分的八个垫脚(eight legs)，在图 5.13(d) 中更加清晰。这些结果是自适应滤波器可以实现的典型结果。如前所述，性能提升的代价是滤波器复杂度的增加。

--------------------------------图 5.13：(a)图像被零均值、方差为 1000 的加性高斯噪声损坏。(b)算术平均值滤波的结果。(c)几何平均值滤波的结果。(d)自适应降噪滤波的结果。所有使用的滤波器尺寸均为 7 × 7。------------------------------------

上述结果使用的 $\sigma_{\eta}^{2}$ 值与噪声方差完全匹配。如果此值未知，且使用的估算值过低，算法将返回与原始图像非常相似的图像，因为校正值会小于应有的数值。过高的估算值会导致方差比被限制在 1.0，并且算法将比正常情况下更频繁地从图像中减去平均值。如果允许使用负值，并在最后重新缩放图像，则会导致动态范围损失，如前所述。

5.3.3.2 自适应中值滤波器 (Adaptive, Local Noise Reduction Filter)

如果椒盐噪声的空间密度较低(根据经验， $P_{s}$ 和 $P_{p}$ 小于 0.2)，公式 (5-27) 中的中值滤波器性能良好。我们将在以下讨论中表明，自适应中值滤波可以处理概率高于这些值的噪声。自适应中值滤波器的另一个优点是，它力求在平滑非脉冲噪声的同时保留细节，而“传统”中值滤波器则无法做到这一点。与前面讨论的所有滤波器一样，自适应中值滤波器也适用于矩形邻域 $S_{xy}$ 。然而，与这些滤波器不同的是，自适应中值滤波器在滤波过程中会根据稍后列出的某些条件改变(增加) $S_{xy}$ 的大小。请记住，滤波器的输出是一个单一值，用于替换 (x, y) 处像素的值，该点是区域 $S_{xy}$ 在已知时间的中心点。

我们使用下述记法：

$z_{\rm min} = S_{xy}$ 中的最小强度值

$z_{\rm max} = S_{xy}$ 中的最大强度值

$z_{\rm med} = S_{xy}$ 中的强度中值

$z_{\rm xy} =$ 坐标 (x, y) 处的强度值

$S_{\rm{max}} = S_{xy}$ 的最大允许大小

自适应中值滤波算法在每一个点 (x, y) 使用两个处理等级，分别使用 A 和 B 表示：

等级 A ：若 $z_{\rm{min}} < z_{\rm{med}} < z_{\rm{max}}$ ，转到等级 B

否则，递增 $S_{xy}$ 的大小

若 $S_{\rm{xy}} \le S_{\rm{max}}$ ，重复等级 A

否则，输出 $z_{\rm{med}}$

等级B ：若 $z_{\rm{min}} < z_{\rm{med}} < z_{\rm{max}}$ ，输出 $z_{xy}$

否则，输出 $z_{\rm{med}}$ 。

其中， $S_{xy}$ 和 $S_{\rm{max}}$ 是奇数，比 1 大的正值。在等级A的最后一步的一个选择是输出 $z_{xy}$ 而不是 $z_{\rm{med}}$ 。这样可以稍微减少模糊程度，但可能无法检测到嵌入在恒定背景中且与椒噪声值相同的盐噪声。

该算法有三个主要目标：消除椒盐噪声(脉冲噪声)、平滑其他非脉冲噪声以及减少失真，例如物体边界过度变细或变粗。从统计学角度来看，算法将 $z_{\rm{min}}$ 和 $z_{\rm{max}}$ 值视为区域 $S_{xy}$ 中的“脉冲状”噪声成分，即使它们不是图像中可能的最低和最高像素值。

在心中记住这些观察结果之后，我们看到，等级 A 的目的是确定是否中值滤波器输出 $z_{\rm{med}}$ 是否为一个脉冲(盐或椒)。若条件 $z_{\rm{min}} < z_{\rm{med}} < z_{\rm{max}}$ 成立，则由于前述段落提到的原因， $z_{\rm{med}}$ 不会是一个脉冲。在这种情况下，我们转到等级 B ，并测试是否邻域中心点其自身是一个脉冲(回顾一下，(x, y) 是正在处理的点的位置，而 $z_{xy}$ 是其强度值)。若条件 $z_{\rm{min}} < z_{xy} < z_{\rm{max}}$ 成立，则 $z_{xy}$ 处的像素不会是一个脉冲的强度，原因同 $z_{\rm{med}}$ 不会是一个脉冲的原因相同。在这种情况下，算法输出不变的像素值 $z_{xy}$ 。通过保持这些“中间等级”的点，滤波后的图像失真减少。若条件 $z_{\rm{min}} < z_{xy} < z_{\rm{max}}$ 不成立，则有 $z_{xy}=z_{\rm{min}}$ 或 $z_{xy}=z_{\rm{max}}$ 。在任一种情况下，这个像素值都是一个极值，而这个算法的输出是中值 $z_{\rm{med}}$ ，而根据等级 A 我们可知这个值不是一个噪声脉冲。最后一步是标准中值滤波器的作用。问题在于，标准中值滤波器会用相应邻域的中值替换图像中的每一个点。这会导致不必要的细节损失。

继续解释，假设 A 层次确实发现了一个脉冲(即，它未能通过导致其分支到 B 层次的测试)。然后，算法会增加邻域的大小并重复 A 层次。此循环持续进行，直到算法找到一个非脉冲的中值(并分支到 B 层)，或者达到最大邻域大小。若达到最大大小，算法返回 $z_{\rm{med}}$ 的值。注意，不保证这个值不是一个脉冲。噪声概率 $P_{a}$ 和/或 $P_{b}$ 越少，或者允许 $S_{\rm{max}}$ 的值越大，则这一层级提前退出的可能性就越小。这是有道理的。随着噪声脉冲密度的增加，我们需要一个更大的窗口来“清理”噪声尖峰，这是理所当然的。

每当算法输出一个值时，邻域中心 $S_{xy}$ 就会移动到图像中的下一个位置。然后，算法重新初始化并应用于邻域所包含的新区域内的像素。如问题 3.37 所示，中值可以从一个位置迭代更新到下一个位置，从而减少计算负载。

例子 5.5：使用自适应中值滤波进行图像去噪。

图 5.14(a) 显示了受椒盐噪声污染的电路板图像，其概率为 $P_{s} = P_{p} = 0.25$ ，其是图 5.10(a) 中所用噪声水平的 2.5 倍。此处的噪声水平足够高，以至于遮挡了图像中的大部分细节。为了进行比较，首先使用 7 × 7 中值滤波器对图像进行滤波，这是在本例中去降低多数可见脉冲噪声痕迹所需的最小滤波器。图 5.14(b) 显示了结果。虽然噪声被有效去除，但滤波器导致图像细节严重丢失。例如，图像顶部的一些指状连接器出现扭曲或断裂。其他图像细节也同样失真。

图 5.14(c) 显示了使用自适应中值滤波器 ( $S_{\rm{max}}=7$ ) 的结果。降噪性能与中值滤波器相似。然而，自适应滤波器在保留清晰度和细节方面表现得更好。指状连接器的失真程度有所降低，而一些其他特征(这些特征被中值滤波器遮挡或扭曲到无法识别)在图 5.14(c) 中显得更加清晰。两个值得注意的例子是遍布电路板的白色小通孔，以及图像左下象限中带有八条腿的暗色元件。

考虑到图 5.14(a) 中的高噪声水平，自适应算法表现相当出色。 $S_{xy}$ 的最大允许尺寸选择取决于具体应用，但可以通过首先试验不同尺寸的标准中值滤波器来估算合理的初始值。这将为自适应算法的性能预期建立一个直观的基准。

----------图5.14：(a)受椒盐噪声破坏的图像，概率为 $P_{s} = P_{p} = 0.25$ 。(b) 使用 7 × 7 中值滤波器进行滤波的结果。(c)自适应中值滤波的结果， $S_{\rm{max}}=7$ 。------------

5.4 使用频域滤波实现周期性噪声降噪(Periodic noise reduction using frequency domain filtering)

使用频域技术可以非常有效地分析和滤除周期性噪声。其基本思想是，周期性噪声在Fourier变换中表现为能量的集中爆发，出现在与周期性干扰频率相对应的位置。方法是使用选择性滤波器(参见第 4.10 节)来隔离噪声。第 4.10 节详细讨论了三种类型的选择性滤波器(带阻滤波器、带通滤波器和陷波滤波器)。这些滤波器在第 4 章中的使用方法与它们在图像恢复中的使用方法并无区别。在修复受周期性干扰损坏的图像时，首选工具是陷波滤波器(a notch filter)。在接下来的讨论中，我们将扩展第 4.10 节中介绍的陷波滤波方法，并开发一种更强大的最优陷波滤波方法。

5.4.1 深入陷波滤波器(More on notch filtering)

如第 4.10 节所述，陷波抑制滤波器的转移函数由高通滤波器转移函数的乘积构成，这些高通滤波器的转移函数中心已平移至陷波的中心。陷波滤波器转移函数的一般形式为

(5-33) $\displaystyle H_{NR}(u ,v) = \prod_{k=1}^{Q}H_{K} (u ,v)H_{-K} (u ,v)$

其中， $H_{K}(u,v)$ 和 $H_{-K}(u,v)$ 是高通滤波器转移函数，其中心分别位于 $(u_{k},v_{k})$ 和 $(u_{-k},v_{-k})$ 。(注：谨记，频域转移函数关于频率矩形的中心对称，因此陷波指定为对称对。另外，回想一下 4.10 节中提到的，在使用陷波滤波器时，我们使用未填充的图像，以简化陷波位置的指定。) 这些中心关于频率矩形的中心 ( floor(M/2) ,floor(N/2)) 而指定，其中，照例，M 和 N 分别表示输入图像的行数和列数。因此，这个滤波器转移函数的距离计算公式为

(5-34) $\displaystyle D_{k}(u ,v) = \left [(u-M/2-u_{k})^{2}+(v-N/2-v_{k})^{2} \right ]^{1/2}$

和

(5-35) $\displaystyle D_{-k}(u ,v) = [(u-M/2+u_{k})^{2}+(v-N/2+v_k)^{2} ]^{1/2}$

例如，以下是具有三个陷波对的 n 阶Butterworth陷波抑制滤波器的转移函数：

(5-36) $\displaystyle H_{NR}(u ,v) = \prod_{k=1}^{3} \left [\frac{1}{(1+[D_{0k}/D_{k} (u ,v)]^{n} )}\right ] \left [\frac{1}{1+[D_{0k}/D_{-k} (u ,v)]^{n} } \right ]$

由于陷波被指定为对称对，因此常数 $D_{0k}$ 对于每对都是相同的。然而，这个常数在不同对之间可能有所不同。其他陷波抑制滤波器函数的构造方式相同，具体取决于所选的高通滤波器函数。如第 4.10 节所述，陷波带通滤波器的转移函数是通过陷波抑制函数使用以下表达式获得的

(5-37) $H_{NP}(u ,v) = 1 - H_{NR}(u ,v)$

其中 $H_{NP}(u ,v)$ 是对应于转移函数为 $H_{NR}(u ,v)$ 的陷波抑制滤波器的陷波带通滤波器的转移函数。图 5.15 分别展示了理想、Gauss 和 Butterworth 陷波抑制滤波器(每个陷波对都为一个)的转移函数透视图。正如我们在第 4 章中讨论过的那样，我们再次看到，Butterworth转移函数的形状代表了理想函数的尖锐度与 Gauss 转移函数的宽阔平滑形状之间的过渡。

正如我们在以下示例的第二部分中所示，我们并不局限于刚才讨论的陷波滤波器转移函数形式。我们可以构建任意形状的陷波滤波器，只要它们是零相移函数，如第 4.7 节中所定义。

-------------------------------图 5.15：(a)理想、(b)Gauss和(c) Butterworth陷波抑制滤波器转移函数透视图。---------------------------------------

例子 5.6：使用陷波滤波对图像进行降噪(减少干扰)。

图 5.16(a) 与图 2.45(a) 相同，我们在 2.6 节中用图 2.45(a) 介绍了频域滤波的概念。现在，我们更详细地研究这幅图像的去噪过程，该图像被一个二维加性正弦波干扰。从表 4.4 可知，纯正弦波的Fourier变换是一对复数共轭脉冲，因此我们预期其频谱在正弦波的频率处会出现一对亮点。如图 5.16(b) 所示，情况确实如此。由于我们可以精确地确定这些脉冲的位置，因此消除它们是一项简单的任务，只需使用一个陷波滤波器，其陷波与脉冲的位置重合即可。

图 5.16(c) 显示了理想的陷波抑制滤波器转移函数，该函数由一个由 1 组成的数组(白色区域)和两个由 0 组成的小圆形区域(黑色区域)组成。图 5.16(d) 显示了使用该转移函数对噪声图像进行滤波的结果。正弦噪声几乎被消除，许多之前被干扰遮挡的细节在滤波后的图像中清晰可见(例如，细小的基准标记以及地形和岩层中的精细细节)。正如我们在示例 4.25 中所示，获取干涉图样的图像非常简单。我们只需将抑制滤波器从 1 中减去，即可将其转换为通过滤波器，然后用它对输入图像进行滤波。图 5.17 显示了结果。

----------------------------------图 5.16：(a) 受正弦干扰损坏的图像。(b) 频谱显示干扰引起的能量爆发。(为了便于显示，能量爆发被放大。)(c) 用于消除能量爆发的陷波滤波器(圆圈半径为 2 像素)。(细边框不属于数据。)(d) 陷波抑制滤波的结果。(原图由 NASA 提供。)------------------------------------------------------------

----------图 5.17：使用陷波滤波器从图 5.16(a) 的 DFT 中提取的正弦波形。--------

图 5.18(a) 显示的图像与图 4.50(a) 相同，但覆盖面积更大(干涉图样相同)。我们在第四章讨论该图像的低通滤波时指出，有更好的方法可以降低扫描线的影响。接下来的陷波滤波方法可以显著减少扫描线，而不会引入模糊。除非出于我们在 4.9 节中讨论过的原因而需要模糊，否则陷波滤波通常会产生更好的效果。

仅通过观察图 5.18(a) 中噪声模式的近乎水平的线条，我们预计其在频域中的贡献会集中在 DFT 的垂直轴上。然而，噪声并不占主导地位，不足以在此轴上形成清晰的模式，正如图 5.18(b) 所示的频谱所示。在这种情况下，应遵循的方法是使用一个沿垂直轴延伸的窄矩形陷波滤波器函数，从而消除沿该轴的所有干扰分量。我们不在原点附近进行滤波，以避免消除直流项和低频，正如你从第四章所知，这些噪声是造成平滑区域之间强度差异的原因。图 5.18(c) 显示了我们使用的滤波器转移函数，图 5.18(d) 显示了滤波结果。大多数细扫描线被消除或显著衰减。为了获得噪声模式的图像，我们像之前一样，将带阻滤波器转换为带通滤波器，然后用它对输入图像进行滤波。图 5.19 显示了结果。

----------------------图 5.18：(a) 佛罗里达州和墨西哥湾的卫星图像。(注意水平传感器扫描线。)(b) (a) 的频谱。(c) 陷波抑制滤波器转移函数。(细黑边框不是数据的一部分。)(d) 滤波后的图像。(原始图像由 NOAA 提供。)---------------------------

--------------图 5.19：通过陷波滤波从图 5.18(a) 中提取的噪声模式。-------------

5.4.2 最优陷波滤波器(Optimum notch filtering)

在迄今为止给出的陷波滤波示例中，干涉图样在频域中很容易识别和表征，从而可以很容易启发式地定义陷波滤波器的转移函数。

当存在多个干扰分量时，启发式的滤波器转移函数规范并不总是可行的，因为它们可能会在滤波过程中去除过多的图像信息(当图像独特且/或获取成本高昂时，这是一个非常不理想的特性)。此外，干扰分量通常不是单频突发，而是往往具有宽阔的裙边状分布，这些裙边状携带着有关干涉图样的信息。这些裙边状并不总是容易从正常的变换背景中检测到。在实践中，减少这些劣化影响的替代滤波方法非常有用。接下来讨论的方法是最优的，因为它最小化了恢复估算 $\hat{f}(x,y)$ 的局部方差。

该过程首先分离干涉图样的主要贡献，然后从损坏的图像中减去该图样的一个可变加权部分。虽然我们是在特定应用的背景下开发该过程的，但其基本方法具有通用性，可以应用于其他存在多重周期性干扰问题的修复任务。

我们首先提取干涉图样的主频分量。与之前一样，我们通过在每个尖峰的位置放置一个陷波滤波器转移函数 $H_{NP}(u ,v)$ 来实现。如果该滤波器的构造使其仅通过与干涉图样相关的分量，则干涉噪声图样的 Fourier 变换由以下表达式给出：

(5-38) $N(u ,v) = H_{NP}(u ,v)G(u ,v)$

其中，照例，G(u ,v ) 是干扰图像的 DFT 。

指定 $H_{NP}(u ,v)$ 需要对什么是干扰尖峰或什么不是进行大量的判断。因此，陷波滤波器通常是通过在显示器上观察 G(u ,v ) 的频谱来交互式构建的。选择特定的滤波函数后，可以使用熟悉的表达式

(5-39) ${\eta}(x ,y) = \frak{F}^{-1}\{ H_{NP}(u ,v)G(u ,v)\}$

来获得空间域中相应的噪声模式。因为假设损坏的图像是由未损坏的图像 f (x ,y) 和干扰 h (x, y) 相加而成的，如果后者完全已知，那么从 g (x, y) 中减去干扰模式以获得 f (x ,y) 将是一件简单的事情。当然，问题在于这种滤波过程通常只能得到真实噪声模式的近似值。为了最小化 h (x, y) 估算中不存在的不完整成分的影响，可以通过从 g (x, y) 中减去 h (x, y) 的加权部分来获得 f (x ,y) 的估算值：

(5-40) $\hat{f}(x ,y) = g(x, y) - w(x ,y){\eta}(x ,y)$

其中，如前一样， $\hat{f}(x ,y)$ 是 f (x ,y) 的估算值，而 w (x ,y) 为待确定值。这个函数称为加权函数(weighting function)或调制函数 (modulation function)，该过程的目标是选择 w (x ,y)，以便使结果按某种有意义的方式最优。一种方法是选择 w (x ,y)，使得 $\hat{f}(x ,y)$ 的方差在每个点 (x ,y) 的指定邻域内最小化。

考虑(奇数)大小为 m×n 中心位于 (x ,y) 的一个邻域 $S_{xy}$ 。 $\hat{f}(x,y)$ 的位于点(x ,y) 的“局部”方差可以通过按 $S_{xy}$ 采样进行如下估算：

(5-41) $\displaystyle \sigma^{2}(x,y)=\frac{1}{mn}\sum_{(r,c){\in}S_{xy}}\left [ \hat{f}(r,c)-\overline{\hat{f}}\right ]^{2}$

其中， $\overline{\hat{f}}$ 是邻域 $S_{xy}$ 内 $\hat{f}$ 的均值；即，

(5-42) $\displaystyle \overline{\hat{f}}=\frac{1}{mn}\sum_{(r,c){\in}S_{xy}}\hat{f}(r,c)$

可以通过考虑部分邻域或用 0 填充边框来处理图像边缘上或边缘附近的点。

将公式 (5-40) 代入公式 (5-41) ，我们得到

(5-43) $\displaystyle \sigma^{2}(x,y)=\frac{1}{mn}\sum_{(r,c){\in}S_{xy}} \bigg \{ \bigg [ g(r,c)-w(r,c)\eta(r,c) \bigg ]-\bigg [ \overline{g}-\overline{w\eta} \bigg ] \bigg\}^{2}$

其中， $\overline{g}$ 和 $\overline{w\eta}$ 分别表示邻域 $S_{xy}$ 内 g 和乘积 wη 的均值。

若我们假设 w 是 $S_{xy}$ 内我们可以在邻域的中心用 w 的值替代 w(r, c) 的值的近似常量：

(5-44) $w(r, c) = w(x, y)$

因为 w(x, y) 假设为是 $S_{xy}$ 内的常量，从而可推断出 $\overline{w}=w(x,y)$ ，即

(5-45) $\overline{w\eta} = w(x, y)\overline{\eta}$

位于 $S_{xy}$ 中，其中， $\overline{\eta}$ 是邻域内 η 的均值。用这些近似值，公式 (5-43) 成为

(5-46) $\displaystyle \sigma^{2}(x,y)=\frac{1}{mn}\sum_{(r,c){\in}S_{xy}} \bigg \{ \bigg [ g(r,c)-w(r,c)\eta(r,c) \bigg ]-\bigg [ \overline{g}-w(x, y)\overline{\eta} \bigg ] \bigg\}^{2}$

为了最小化 $\sigma^{2}(x,y)$ 关于w(x, y) 的值，我们对 w(x, y) 解方程

(5-47) $\displaystyle \frac{{\partial }\sigma^{2}(x,y)}{{\partial}w(x,y)}=0$ 。其解是(见问题 5.17):

(5-48) $\displaystyle w(x, y) = \frac{\overline{g\eta}-\overline{g}{\;}\overline{\eta}}{\overline{\eta^{2}}-\overline{\eta}^{2}}$

为了获得恢复图像在点 (x, y) 处的值，我们使用此方程计算 w(x, y)，然后将其代入公式 (5-40)。为了获得完整的恢复图像，我们对噪声图像 g 中的每一个点都执行此过程。

例子 5.7：使用最佳陷波滤波进行去噪(消除干扰)。

图 5.20(a) 展示了水手 6 号航天器拍摄的火星地形数字图像。该图像被半周期干涉图样破坏，这种图样比我们迄今为止研究的干涉图样要复杂得多(也微妙得多)。如图 5.20(b) 所示，该图像的Fourier频谱中存在许多由干涉引起的“星状”能量爆发。不出所料，这些成分比我们之前见过的成分更难检测到。图 5.21 再次显示了该频谱，但没有进行中心化处理。这幅图可以更清晰地显示干涉成分，因为更突出的直流项和低频“隐藏”在频谱的左上角。

图 5.22(a) 显示了经验丰富的图像分析师判断与干扰相关的光谱成分。对这些成分应用陷波滤波器，并利用公式 (5-39) 得到空间噪声模式 η(x, y)，如图 5.22(b) 所示。请注意，该模式与图 5.20(a) 中噪声结构的相似性。

--------------------------图 5.20：(a)水手6号拍摄的火星地形图像。(b)傅里叶频谱显示周期性干扰。(图片由 NASA 提供)-------------------------------------------

-------图 5.21：图 5.20(a) 中图像的非中心 Fourier 频谱。(图片由 NASA 提供)----

----------------------图 5.22：(a)N (u,v ) 的 Fourier 频谱和(b)相应的空间噪声干涉图样 η(x, y)。(图片由 NASA 提供。)---------------------------------------

最后，图 5.23 显示了使用公式 (5-40) 获得的修复图像，其中包含刚才讨论的干涉图样。函数 w(x, y) 是使用前几段中解释的程序计算的。正如你所见，图 5.20(a) 中的噪声图像中的周期性干扰几乎被消除了。

---------------------图 5.23：修复后的图像。(图片由 NASA 提供)---------------

5.5 线性位置不变退化(Linear, position-invariant degradations)

图 5.1 中在恢复阶段前的输入输出关系表述为

(5-49) $g(x , y) = \mathcal{H}( f(x , y) ) + {\eta}( x , y )$

这时，我们假设 η ( x , y ) = 0，因此 g (x , y) = ℋ( f (x , y) ) 。基于 2.6 节中的讨论，若

(5-50) $\mathcal{H}[ af_{1}(x , y) + bf_{2}(x , y) ] = a\mathcal{H}[f_{1}(x , y)] + b\mathcal{H}[ f_{2}(x , y))]$

(其中，a 和 b 是标量，而 $f_{1}(x,y)$ 和 $f_{2}(x,y)$ 是任意两幅输入图像)，则 ℋ 是线性的。

若 a = b = 1 ，则公式 (5-50) 成为

(5-51) $\mathcal{H}[ f_{1}(x , y) + f_{2}(x , y) ] = \mathcal{H}[f_{1}(x , y)] + \mathcal{H}[ f_{2}(x , y))]$

这称为可加性(additivity)。这个属性表明，若 ℋ 是一个线性算子，则对两个输入之和的响应等于两个响应之和。

若 $f_{2}(x,y)=0$ ，则公式 (5-50) 成为

(5-52) $\mathcal{H}[ af_{1}(x , y)] = a \mathcal{H}[f_{1}(x , y)]$

这称为齐性(homogeneity)。它表示对任何输入的常数倍数的响应等于对该输入乘以相同常数的响应。因此，线性算子既具有可加性，又具有齐性。

对于具有输入输出关系 g (x , y) = ℋ[ f (x , y)] 的一个算子，若对于任意 f (x , y) 和任意两个标量 α 和 β ，都有

(5-53) $\mathcal{H}[f (x - \alpha , y - \beta)] = g(x - \alpha , y - \beta )$

则称此算子是位置(或空间)不变的。这个定义表明，图像中任何一点的响应仅取决于该点的输入值，而不取决于其位置。

利用二维连续脉冲的筛选性质[参见公式 (4-55)]，我们可以将 f (x , y) 写为

(5-54) $\displaystyle f (x , y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\alpha ,\beta) \delta(x-\alpha,y-\beta)d{\alpha}d\beta$

再次假设 η ( x , y ) = 0，将这个公式代入公式 (5-49)得到

(5-55) $\displaystyle g(x , y) = \mathcal{H}[ f (x , y)] = \mathcal{H} \left [ \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f (\alpha ,\beta)\delta(x-\alpha,y-\beta)d{\alpha}d\beta \right ]$

若 ℋ 是线性算子且我们将可加性扩展到积分，则

(5-56) $\displaystyle g(x , y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\mathcal{H} \left [f (\alpha ,\beta)\delta(x-\alpha,y-\beta) \right ]d{\alpha}d\beta$

因为 f (α , β ) 独立于 x 和 y ，并利用齐性可推断出

(5-57) $\displaystyle g(x , y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f (\alpha ,\beta)\mathcal{H} \left [\delta(x-\alpha,y-\beta) \right ]d{\alpha}d\beta$

项

(5-58) $h(x , \alpha , y , \beta ) = \mathcal{H}[\delta( x - \alpha , y - \beta )]$

称为 ℋ 的脉冲响应。换言之，若公式 (5-49) 中的 η( x , y ) = 0 ，则 h(x , α , y , β ) 是一个脉冲在坐标 ( x , y ) 处的 ℋ 的响应。在光学中，脉冲会变成一个光点，h(x , α , y , β ) 通常称为点扩散函数 (PSF(point spread function))。这个名称基于这样一个事实：所有物理光学系统都会在一定程度上模糊(扩散)光点，而模糊的程度取决于光学元件的质量。

将公式(5.58)代入公式(5.57)就得到表达式

(5-59) $\displaystyle g(x , y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f (\alpha ,\beta)h(x ,\alpha ,y ,\beta )d{\alpha}d\beta$

这称为第一类叠加(或 Fredholm)积分。这个表达式是线性系统理论的核心结论。它指出，如果已知 ℋ 对脉冲的响应，则可以使用公式 (5-59) 计算出对任意输入 f (a ,b) 的响应。换句话说，线性系统 ℋ 完全由其脉冲响应表征。

若 ℋ 是位置不变的，则从公式 (5-53)可推断出

(5-60) $\displaystyle \mathcal{H}[ \delta( x -\alpha , y - \beta )] = h( x - \alpha , y - \beta )$

在这种情况下，公式 (5-59) 归结为

(5-61) $\displaystyle g(x , y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\alpha ,\beta) h(x - \alpha ,y - \beta )d{\alpha}d\beta$

该表达式是公式 (4-24) 中引入的针对单变量的卷积积分，并在问题 4.19 中扩展到二维。公式 (5-61) 告诉我们，对于任何输入，线性、位置不变系统的输出都是通过对输入和系统的脉冲响应进行卷积得到的。

在存在加性噪声的情况下，线性退化模型[公式(5-59)]的表达式变为

(5-62) $\displaystyle g(x , y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\alpha ,\beta) h(x ,\alpha ,y ,\beta )d{\alpha}d\beta + \eta( x , y )$

若 ℋ 是位置不变的，则这个公式成为

(5-63) $\displaystyle g(x , y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(\alpha ,\beta) h(x -\alpha ,y -\beta )d{\alpha}d\beta + \eta( x , y )$

噪声项 η ( x , y ) 的值是随机的，并且假设与位置无关。使用第 3 章和第 4 章中介绍的卷积符号，我们可以将公式 (5-63) 写为

(5-64) $g(x , y) = (h \bigstar f)(x , y) + \eta( x , y )$

或者，利用卷积定理，我们将频域中的等效结果写为

(5-65) $G(u ,v) = H(u ,v)F(u ,v) + N(u ,v)$

这两个表达式与公式 (5-1) 和 (5-2) 一致。请记住，对于离散量，所有乘积都是逐元素之积，如第 2.6 节所定义。

综上所述，前面的讨论表明，一个具有加性噪声的线性、空间不变的退化系统可以在空域中建模为图像与系统退化(点扩散)函数的卷积，然后再添加噪声。基于卷积定理，同样的过程可以在频域中表示为图像变换与退化函数的乘积，然后再添加噪声变换。在频域中，我们使用 FFT 算法。然而，与第 4 章不同，在本章讨论的任何频域恢复滤波器的实现中，我们都不采用图像填充。原因是，在恢复工作中，我们通常只能访问退化图像。为了使填充有效，必须在图像退化之前对其进行处理，而这在实践中显然无法满足。如果我们能够访问原图像，那么恢复就毫无意义了。

许多类型的退化可以用线性、位置不变的过程来近似。这种方法的优势在于，线性系统理论的大量工具可以用来解决图像恢复问题。非线性和位置相关技术虽然更通用(通常也更准确)，但也引入了一些难题，这些难题通常没有已知的解决方案，或者很难通过计算解决。本章重点介绍线性、空间不变的恢复技术。由于退化建模为卷积的结果，而恢复则试图找到反向应用该过程的滤波器，因此“图像反卷积”一词经常用来表示线性图像恢复。同样，恢复过程中使用的滤波器通常被称为反卷积滤波器。

5.6 估计退化函数(Estimating the degradation function)

用于图像恢复的退化函数估计方法主要有三种：(1) 观察法；(2) 实验法；(3) 数学建模法。以下章节将分别讨论这些方法。使用通过上述方法估计出的退化函数来恢复图像的过程有时称为盲反卷积(blind deconvolution)，以强调真正的退化函数鲜少完全知晓。

5.6.1 通过图像观察估计(Estimation by image observation)

假设我们得到一张退化图像，但对退化函数 ℋ 一无所知。基于图像由线性、位置不变过程退化的假设，一种估计 ℋ 的方法是从图像本身收集信息。例如，如果图像模糊，我们可以查看图像中一小块矩形区域，其中包含样本结构，例如物体的一部分和背景。为了降低噪声的影响，我们会寻找信号内容较强的区域(例如，高对比度区域)。下一步是处理子图像，以获得尽可能清晰的结果。

令 $g_{s}(x ,y)$ 表示观察的子图像，并令处理过的子图像为 $\hat{f}_{s}$ (这实际上是我们对该区域原图像的估计)。然后，假设由于我们选择了强信号区域，噪声的影响可以忽略不计，根据公式 (5-65) 可知

(5-66) $\displaystyle H_{s}(x,y)=\frac{G_{s}(u,v)}{\hat{F}_{s}(u,v)}$

然后，我们根据该函数的特征，基于位置不变性假设，推导出完整的退化函数 H(u ,v)。

例如，假设 $H_{s}(u ,v)$ 的径向图近似于 Gauss 曲线。我们可以利用该信息构建一个更大尺度但基本形状相同的函数 H(u ,v)。之后，我们将 H(u ,v) 用于以下章节将要讨论的其中一种恢复方法中。显然，这是一个费力的过程，仅在非常特殊的情况下使用，例如修复具有历史价值的老照片。

5.6.2 通过实验估计(Estimation by experimentation)

如果有与用于获取退化图像的设备类似的设备，原则上可以准确估计退化程度。可以使用各种系统设置获取与退化图像类似的图像，直到其退化程度尽可能接近我们想要恢复的图像。然后，我们的想法是，通过使用相同的系统设置对脉冲(小光点)进行成像，来获得退化的脉冲响应。如第 5.5 节所述，线性、空间不变系统完全由其脉冲响应表征。

脉冲可以用一个明亮的光点来模拟，该光点尽可能明亮，以将噪声的影响降至可忽略不计的程度。然后，考虑到脉冲的 Fourier 变换是一个常数，由公式 (5-65) 可知

(5-67) $\displaystyle H(u ,v) = \frac{G(u ,v)}{A}$

其中，如前一样，G(u ,v) 是受观察图像的 Fourier 变换，而 A 是描述脉冲强度的常量。图 5.24 展示了一个例子。

---------------------------------图 5.24：通过脉冲特性评估衰减。(a) 光脉冲(放大显示)。(b) 成像(衰减)脉冲。----------------------------------------------------

5.6.3 通过建模估计(Estimation by modeling)

退化模型因其对图像恢复问题的深刻理解而被应用多年。在某些情况下，该模型甚至可以将导致退化的环境条件考虑在内。例如，Hufnagel 和 Stanley [1964] 提出的退化模型基于大气湍(tuān)流(turbulence)的物理特性。该模型具有以下常见形式：

(5-68) $\displaystyle H(u ,v) = e^{-k(u^{2}+v^{2})^{5/6}}$

其中 k 是一个常数，取决于湍流的性质。除了指数中的 5 次方之外，该公式的形式与 4.8 节中讨论的 Gauss 低通滤波器转移函数相同。事实上，Gauss低通滤波器有时用于模拟轻微、均匀的模糊。图 5.25 展示了使用公式 (5-68) 模拟图像模糊的示例，其中 k = 0.0025(严重湍流)、k = 0.001(轻微湍流)和 k = 0.00025(低湍流)。我们将在本章后面使用各种方法恢复这些图像。

-------------------------------图 5.25：湍流建模。(a)无可见湍流。(b)强烈湍流，k = 0.0025。(c)轻微湍流，k = 0.001。(d)低湍流，k = 0.00025。所有图像尺寸均为 480 × 480 像素。(原图像由 NASA 提供。)------------------------------------

建模中常用的另一种方法是从基本原理出发，推导出一个数学模型。我们通过这样的方式来说明这一点，即通过详细讨论在图像采集过程中，由于图像和传感器之间发生匀速直线运动而导致图像模糊的情况来说明这一过程。假设图像 f ( x , y ) 进行平面运动， $x_{0}(t)$ 和 $y_{0}(t)$ 分别是 x 和 y 方向运动的时变分量。我们通过将瞬时曝光量与成像系统快门打开的时间间隔进行积分，可以得到记录介质(例如胶片或数字存储器)上任意一点的总曝光量。

假设快门的开启和关闭是瞬间发生的，并且光学成像过程完美无缺，我们可以分离出图像运动的影响。然后，设曝光时长为 T，则可以推导出

(5-69) $\displaystyle g(x , y) = \int_{0}^{T}f \bigg [ x-x_{0} (t),y-y_{0} (t) \bigg ] dt$

其中，g (x , y ) 是模糊图像。这个表达式的连续 Fourier 变换是

(5-70) $\displaystyle G( u , v) =\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}g(x ,y) e^{-j2{\pi}(ux+vy)} dxdy$

将公式 (5-69) 代入公式 (5-70) 得到

(5-71) $\displaystyle G( u , v) =\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\bigg [ \int_{0}^{T}f \left [ x-x_{0} (t),y-y_{0} (t) \right ] dt \bigg ] e^{-j2{\pi}(ux+vy)} dxdy$

在表达式中调换积分次序得到

(5-72) $\displaystyle G( u , v) =\int_{0}^{T}\bigg [\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} f \left [ x-x_{0} (t),y-y_{0} (t) \right ]e^{-j2{\pi}(ux+vy)} dxdy \bigg ]dt$

在外部括号内的这一项是位移函数 $f ( x - x_{0}(t), y - y_{0}(t))$ 的 Fourier 变换，则使用表 4.4 中的条目 3 得出表达式

(5-73) $\begin{array}{rl} G( u , v) &\displaystyle = \int_{0}^{T}F(u ,v)e^{-j2{\pi}[ux_{0} (t)+vy_{0} (t)]}dt \\ \\ &=F( u , v)\int_{0}^{T}e^{-j2{\pi}[ux_{0} (t)+vy_{0} (t)]}dt \end{array}$

通过定义

(5-74) $\displaystyle H(u,v)=\int_{0}^{T}e^{-j2{\pi}[ux_{0} (t)+vy_{0} (t)]}dt$

我们可以按熟悉的形式表达(5-73)

(5-75) $G( u , v) = H( u , v)F( u , v)$

若运动变量 $x_{0}(t)$ 和 $y_{0}(t)$ 是已知的，转移函数 H( u , v) 可以直接从公式 (5-74) 而获得。例如，假设所讨论的图像仅在 x 方向上以速率 $x_{0}(t)= at/T$ 进行均匀线性运动 ( 即 $y_{0}(t) = 0$ ) 。当 t = T 时，图像总共位移了距离 a 。当 $y_{0}(t) = 0$ 时，公式 (5-74)变成

(5-76) $\begin{array}{rl} H(u,v)&\displaystyle=\int_{0}^{T}e^{-j2{\pi}ux_{0}(t)}dt \\ \\ &\displaystyle=\int_{0}^{T}e^{-j2{\pi}uat/T}dt \\ \\ &\displaystyle=\frac{T}{{\pi}ua}\sin{({\pi}ua)}e^{-j{\pi}ua} \end{array}$

若我们也允许 y 分量变动，已知运动由 $y_{0}(t) = bt/T$ 描述，则退化函数成为

(5-77) $\displaystyle H( u , v) = \frac{T}{{\pi}(ua+vb)}\sin[{\pi}(ua+vb)]e^{-j{\pi}(ua+vb)}$

为了产生大小为 M × N 滤波器转移函数，我们对这个公式按 u = 0, 1, 2, …, M - 1 和 v = 0 , 1 ,2 ,… ，N - 1 进行采样。

例子 5.8：由于运动导致的图像模糊。

图 5.26(b) 是通过对图 5.26(a) 中的图像进行 Fourier 变换，乘以公式 (5-77) 中的 H( u , v)，然后进行逆变换而得到的图像。此图像大小为 688 × 688 像素，在公式 (5-77) 中，我们使用 a = b = 0.1 和 T = 1。正如我们将在 5.8 和 5.9 节中讨论的那样，从模糊图像中恢复原图像面临着一些有趣的挑战，尤其是在退化图像中存在噪声的情况下。正如 5.5 节末尾所述，我们执行所有 DFT 计算时均不使用填充。