数字图像处理(4版)——第 4 章——频域滤波(下)(Rafael C.Gonzalez&Richard E. Woods)

我们首先在公式 (4-67) 中观察到，F(u,v) 的每一项都包含 f(x, y) 的所有值，这些值均由指数项的值修正。因此，除了一些简单情况，通常不可能将图像的特定分量与其变换直接关联。然而，可以对Fourier变换的频率分量与图像的空间特征之间的关系做出一些一般性的表述。例如，由于频率与空间变化率直接相关，因此直观地将Fourier变换中的频率与图像中强度变化的模式联系起来并不困难。我们在第 4.6 节中证明，变化最慢的频率分量 (u = v = 0) 与图像的平均强度成正比。当我们远离变换的原点时，低频对应于图像中缓慢变化的强度分量。例如，在房间的图像中，这些低频可能对应于墙壁和地板上平滑的强度变化。随着我们离原点越来越远，更高的频率开始对应图像中越来越快的强度变化。这些是物体的边缘以及图像中其他以强度突变为特征的组成部分。

频域滤波技术基于修改Fourier变换以实现特定目标，然后计算逆DFT返回空域，如2.6节所述。根据公式(4-87)，我们可以获得的变换的两个分量是变换幅度(频谱)和相位角。我们在4.6节中了解到，对相位分量进行目视分析通常意义不大。然而，频谱可以为生成频谱的图像的总体强度特性提供一些有用的指导。例如，考虑图4.28(a)，它是一张集成电路的扫描电子显微镜图像，放大了约2500倍。

除了器件本身有趣的结构之外，我们还注意到该图像的两个主要特征：大约 ±45°的强边缘，以及两个由热致失效引起的白色氧化物突起。图4.28(b)中的 Fourier 谱显示了沿±45°方向的显著分量，这些分量与刚才提到的边缘相对应。仔细观察图4.28(b)中的垂直轴，我们可以看到变换的垂直分量偏离轴心，略微偏左。该分量是由氧化物突起的边缘引起的。请注意，频率分量相对于垂直轴的角度与长白色元件的倾斜度(相对于图像的水平轴)相对应。还要注意垂直频率分量中的零点，它们对应于氧化物突起的狭窄垂直跨度。

这些是我们通常在频域和空域之间建立的典型关联类型。正如我们将在本章后面展示的那样，即使是这些类型的粗略关联，结合前面提到的图像中频率成分和强度变化率之间的关系，也能得出一些非常有用的结果。我们将在4.8节中展示修改图4.28(a)变换中不同频率范围的效果。

--------------------------------图 4.28：(a) 受损集成电路的 SEM 图像。(b) (a) 的Fourier谱。(原始图像由加拿大Ontario省Hamilton市McMaster大学Brockhouse材料研究所的 J. M. Hudak 博士提供。)----------------------------------------------

4.7.2 频域滤波的基础知识 (Frequency domain filtering fundamentals)

频域滤波包括修改图像的Fourier变换，然后计算其逆变换以获得处理结果的空域表示。因此，给定一个(填充的)大小为P × Q 个像素的数字图像 f (x, y)，我们感兴趣的基本滤波公式的形式如下：𝔉

(4-104) $g (x, y) = {\rm{Real}} \{\frak{F}^{-1} [H( u ,v)F( u ,v)]\}$

其中， $\frak{F}^{-1}$ 是 IDFT ，F( u ,v) 是输入图像 f (x, y) 的 DFT ，H( u ,v) 是一个滤波转移函数(我们通常称为一个滤波器或滤波函数)，而 g (x, y) 是滤波(输出)图像。函数 F，H 和 g 是大小为 P × Q 的数组，与填充后的输入图像大小相同。乘积 H( u ,v)F( u ,v) 通过元素乘法生成，定义见 2.6 节。滤波器转移函数修改输入图像的变换，得到处理后的输出图像 g(x, y)。使用关于其中心对称的函数可以大大简化 H( u ,v) 的指定任务，这要求 F( u ,v) 也是关于其中心对称的。如 4.6 节所述，这可以通过在计算输入图像的变换之前将其乘以 $(-1)^{x+y}$ 来实现。(注：一些二维DFT的软件实现(例如MATLAB)不进行中心化变换。这意味着滤波函数必须按照与非中心变换相同的数据格式排列(即原点位于左上角)。最终结果是滤波器转移函数更难生成和显示。我们在讨论中使用中心化是为了帮助可视化，这对于清晰地理解滤波概念至关重要。只要保持一致性，任何一种方法都可以在实践中使用。)

现在我们可以详细讨论滤波了。我们可以构造的最简单的滤波转移函数之一是函数 H( u ,v)，它在(中心)变换的中心处为 0，其他位置处为 1。当我们计算 H( u ,v)F( u ,v) 的乘积时，该滤波器会拒绝直流项，并“通过”(即保持不变) F( u ,v) 的所有其他项。根据表 4.3 中的属性 7，我们知道直流项决定了图像的平均强度，因此将其设置为零会使输出图像的平均强度降低为零。图 4.29 显示了使用公式 (4-104) 进行此操作的结果。正如预期的那样，图像变得更暗了。平均值为零意味着存在负强度。因此，尽管图 4.29 说明了这一原理，但它并不是原图像的真实表示，因为所有负强度都被显示器截断了(设置为 0)。

--------------------------------图 4.29：使用滤波器转移函数对图 4.28(a) 中的图像进行滤波的结果，该函数将中心Fourier变换中的直流项 F(P/2 ,Q/2) 设置为 0，同时保持所有其他变换项不变。---------------------------------

如前所述，变换中的低频与图像中缓慢变化的强度成分相关，例如房间的墙壁或室外场景中万里无云的天空。在另一方面，高频是由强度的急剧转变引起的，例如边缘和噪声。因此，我们预期，一个衰减高频、允许低频的函数 F( u ,v)(如前所述，称为低通滤波器)会使图像模糊，而具有相反特性的滤波器(称为高通滤波器)会增强清晰的细节，但会导致图像对比度降低。图 4.30 说明了这些效果。例如，该图的第一列显示了低通滤波器的转移函数和相应的滤波图像。第二列显示了高通滤波器的类似结果。请注意图 4.30(e) 和图 4.29 之间的相似性。原因是所示的高通滤波器函数消除了直流项，从而产生了与图 4.29 相同的基本效果。如第三列所示，在滤波器中添加一个小常数，不会显著影响锐化效果，但它确实阻止了直流项的消除，从而保留了调性(tonality)。

-------------------图 4.30：上行：(a) 低通滤波器、(b) 高通滤波器和 (c) 偏移高通滤波器的频域滤波器转移函数。下行：使用公式 (4-104) 获得的相应滤波图像。(c) 中的偏移量为 a = 0.85，H( u ,v) 的高度为 1。将 (f) 与图 4.28(a) 进行比较。---------

公式 (4-104) 涉及频域中两个函数的乘积，根据卷积定理，这意味着空域中发生了卷积。从 4.6 节的讨论中我们知道，如果所讨论的函数未填充，则可以预期会出现环绕误差。图 4.31 显示了当我们应用不带填充的公式 (4-104) 时会发生什么。图 4.31(a) 显示了一幅简单图像，图 4.31(b) 是用图 4.30(a) 所示形式的Gauss低通滤波器对图像进行低通滤波的结果。正如预期的那样，图像变得模糊。然而，模糊并不均匀；顶部的白色边缘模糊，但两侧并不模糊。在应用公式 (4-104) 之前，根据公式 (4-98) 和 (4-99) 对输入图像进行零填充，得到了图 4.31(c) 中的滤波图像。这个结果正如预期的那样，由于零填充，边框呈现出了均匀的深色(有关此效果的解释，请参见图 3.33)。

-------------图4.31：(a)一张简单图像。(b)使用无填充的Gauss低通滤波器进行模糊处理的结果。(c)使用零填充的低通滤波的结果。比较(b)和(c)中的垂直边缘。----------

图 4.32 说明了图 4.31(b) 和 (c) 之间差异的原因。图 4.32(a) 中的虚线区域对应于图 4.31(a) 中的图像。图像的其他副本是由于使用 DFT 时隐含的图像(及其变换)的周期性所致，如第 4.6 节所述。想象一下，将模糊滤波器的空间表示(即相应的空间核)与该图像进行卷积。当核位于虚线图像的顶部中心时，它将包含图像的一部分以及其正上方周期图像底部的一部分。当滤波器下方存在暗区和亮区时，结果将是中灰色的模糊输出。然而，当核位于图像的右上方中心时，它将仅包含图像中的亮区及其右侧区域。由于常数值的平均值与该值相同，因此滤波在此区域不会产生任何效果，结果如图 4.31(b) 所示。用 0 填充图像会在周期序列的每个图像周围创建均匀的边框，如图 4.32(b) 所示。将模糊核与图 4.32(b) 中填充后的“马赛克”进行卷积，会在图 4.31(c) 中得到正确的结果。从此例中可以看出，在滤波之前未填充图像可能会导致意外的结果。

-------------------图 4.32：(a)未进行图像填充的图像周期性。(b)用 0(黑色)填充后的周期性。中间的虚线区域对应于图 4.31(a) 中的图像。使用 DFT 时，周期性是固有的。(为了清晰起见，两幅图中的细白线是叠加的；它们不是数据的一部分。)--------------

到目前为止，讨论主要集中在对输入图像进行填充。然而，公式 (4-104) 还涉及一个滤波器转移函数，该函数可以在空域或频域中指定。填充是在空域进行的，这引出了一个重要的问题：空间填充与直接在频域中指定的滤波器函数之间的关系。

合理的结论是，处理频域转移函数填充的方法是：构造与未填充图像大小相同的函数，计算该函数的离散Fourier变换(IDFT)以获得相应的空间表示，在空域中填充该表示，然后计算其离散Fourier变换 (DFT) 返回频域。图 4.33 中的一维示例说明了这种方法的缺陷。

------------------------图 4.33：(a) 在(中心)频域中指定的滤波器转移函数。(b) 通过计算 (a) 的离散Fourier变换 (IDFT) 获得的空间表示(滤波器核)。(c) 将 (b) 填充至其长度的两倍后的结果(注意不连续性)。(d) 通过计算 (c) 的离散Fourier变换 (DFT) 获得的频域中对应的滤波器。注意 (c) 中不连续性引起的振铃效应(ringing)。图的 (b) 部分位于 (a) 下方，(d) 位于 (c) 下方。----------------------------------------

图 4.33(a) 显示了频域中一维理想低通滤波器的转移函数。该函数为实函数，且具有偶对称性，因此根据表 4.1 中的属性 8，我们知道它的离散Fourier变换 (IDFT) 也为实函数且对称。图 4.33(b) 显示了将转移函数元素乘以 $(-1)^{u}$ 并计算其离散Fourier变换 (IDFT) 的结果，从而得到相应的空间滤波器核。结果如图 4.33(b) 所示。从图中可以明显看出，该空间函数的极值不为零。对函数进行零填充会产生两个不连续点，如图 4.33(c) 所示。回到频域，我们计算空间填充函数的正向离散Fourier变换 (DFT)。如图 4.33(d) 所示，填充函数中的不连续点导致其频域对应函数出现振铃效应。

前述结果表明，我们无法通过填充频域转移函数的空间域表示来避免环绕误差。我们的目标是在频域中使用指定的滤波器形状，而无需担心截断问题。另一种方法是先填充图像，然后直接在频域中创建所需的滤波器转移函数，该函数的大小与填充后的图像相同(谨记，使用DFT时，图像和滤波器转移函数的大小必须相同)。当然，这会导致环绕误差，因为滤波器转移函数没有使用填充，但填充图像提供的分离可以显著减轻这种误差，而且它比振铃效应更可取。平滑的转移函数(例如图4.30中的转移函数)的问题更小。具体来说，本章将采用的方法是将图像填充到P × Q 个像素的大小，并直接在频域中构建相同维度的滤波器转移函数。如前所述，P和Q由公式(4-100)和(4-101)给出。

我们通过分析滤波图像的相位角来结束本节。我们可以用实部和虚部来表示DFT：F(u ,v) = R(u ,v) + jI(u ,v)。这样，公式(4-104)就变为

(4-105) $g (x, y) = \frak{F}^{-1} [H( u ,v)F( u ,v) + j H( u ,v)I(u ,v) ]$

相位角计算为复数虚部与实部之比的反正切[参见公式 (4-88)]。由于 H( u ,v) 同时乘以 R 和 I，因此当该比率成立时，两者会相互抵消。对实部和虚部影响相同，因此对相位角没有影响的滤波器贴切地称为零相移滤波器。本章仅讨论这类滤波器。

图 4.26 生动地展示了相位角在确定图像空间结构方面的重要性。因此，即使相位角发生微小变化，也会对滤波输出产生显著(且通常不理想的)影响，这也不足为奇。图 4.34(b) 和 (c) 展示了改变图 4.34(a) 中离散Fourier变换 (DFT) 的相位角阵列后的效果(两种情况下 |F(u ,v)| 项均未改变)。图 4.34(b) 是通过将公式 (4-86) 中的相位角 𝜙(u ,v) 乘以 −1 并计算离散Fourier变换 (IDFT) 得到的。最终结果是图像中每一个像素绕两个坐标轴的反射。图 4.34(c) 是通过将相位项乘以 0.25 并计算离散Fourier变换 (IDFT) 得到的。即使是尺度变化，图像也几乎无法辨认。这两个结果说明了使用不改变相位角的频域滤波器的优势。

----------------------------------图 4.34：(a)原始图像。(b)在公式(4-86)中将相位角数组乘以 −1 并计算 IDFT 后得到的图像。(c)将相位角乘以 0.25 并计算 IDFT 后的结果。(b)和(c)中使用的变换 F (u , v) 的幅度相同。-------------------------

4.7.3 频域滤波步骤总结 (Summary of steps for filtering in the frequency domain)

频域滤波的过程可以概括如下：

(1) 给定大小为 M × N 的输入图像 f (x, y)，使用公式 (4-102) 和 (4-103) 得出填充大小 P 和 Q；即 P = 2M 和 Q = 2N。

(2) 使用零填充、镜像填充或复制填充，形成大小为 P × Q 的填充图像 $f_{p}(x,y)$ (有关填充方法的比较，请参见图 3.39)。(注：有时，我们在进行“快速”实验以了解滤波器性能时，或者在尝试确定空间特征与其对频域分量的影响之间的定量关系时，会省略填充，特别是在带内滤波和陷波滤波中，如后面第 4.10 节和第 5 章所述。)

(3) 将 $f_{p}(x,y)$ 乘以 $(-1)^{x+y}$ 以使 Fourier 变换位于 P × Q 频率矩形的中心。

(4) 计算步骤 (3) 中图像的 DFT ，得到其结果 F(u ,v)。

(5) 构造一个实数对称滤波器转移函数 H(u ,v)，大小为 P × Q ，中心位于 (P/2, Q /2)。

(6) 使用元素乘法形成乘积 G(u ,v) = H(u ,v)F(u ,v)；即，G(i ,k) = H(i ,k)F(i ,k) 其中 i = 0, 1 , 2 ,… , M - 1 和 k = 0 , 1, 2 ,… , N – 1 。

(7) 通过计算 G(u ,v) 的 IDFT 获得滤波图像(大小为 P × Q )：

$g_{p} (x, y) = ({\rm{real}}[\frak{F}^{-1}\{ G(u ,v)\}])(-1)^{x+y}$

(8) 通过从 $g_{p} (x, y)$ 的左上象限提取 M × N 区域，获得与输入图像大小相同的最终滤波结果 g(x, y)。

我们将在本章的后续章节中讨论滤波器转移函数的构造(步骤 (5))。理论上，步骤 (7) 中的 IDFT 应该是实数，因为 f (x, y) 是实数，而 H (u ,v) 是实数且对称。然而，由于计算误差，IDFT 中出现寄生复数项的情况并不少见。取结果的实部可以解决这个问题。乘以 $(-1)^{x+y}$ 以抵消步骤 (3) 中乘以该因子的运算。

图 4.35 展示了使用零填充的上述步骤。图例解释了每一幅图像的来源。如果放大，图 4.35(c) 将显示图像中交错的黑点，这是因为 $f_{p}$ 与 $(-1)^{x+y}$ 相乘得到的负强度值被显示器截断为 0。请注意图 4.35(h) 中使用零填充获得的低通滤波图像的典型暗边框。

---------------------------------图 4.35：(a) 一幅 M × N 图像 f 。(b) 大小为 P × Q 的填充图像 $f_{p}$ 。 (c) 用 $(-1)^{x+y}$ 乘以 $f_{p}$ 后的结果。 (d) F 的频谱。(e) 大小为 P × Q 的中心对称 Gauss 低通滤波器转移函数 H 。( f ) 乘积 HF 的频谱。( g ) 用 $(-1)^{x+y}$ 乘以 HF 的 IDFT 的实部后的图像 $g_{p}$ ( h ) 通过提取 $g_{p}$ 的前M行N列得到最终结果g 。--------------------------------------------------------

4.7.4 空间域和频域滤波之间的对应关系 (Correspondence between filtering in the spatial and frequency domains)

如前所述，空间域滤波和频域滤波之间的联系在于卷积定理。在本节前面，我们将频域滤波定义为滤波器转移函数 H(u ,v) 与 F(u, v)(输入图像的Fourier变换)的元素乘积。给定 H(u ,v)，假设我们想在空域中求得它的等价核。设 f (x, y) = d(x, y)，则根据表 4.4 可知 F(u, v) = 1 。则根据公式 (4-140) ，滤波输出是 $\frak{F}\{ H(u,v)\}$ 。该表达式是频域滤波器转移函数的逆变换，其对应的是空域中的核。反过来，根据类似的分析和卷积定理，给定一个空间滤波器核，我们可以通过对该核进行正向Fourier 变换来获得其频域表示。因此，这两个滤波器构成一个Fourier变换对：

(4-106) $h( x ,y) \Longleftrightarrow H(u ,v)$

其中，h( x ,y) 是空间核。由于该核可以从频域滤波器对脉冲的响应中获得，因此 h( x ,y) 有时称为 H(u ,v) 的脉冲响应。此外，由于 (4-106) 式的离散实现中的所有量都是有限的，因此此类滤波器称为有限脉冲响应(FIR)滤波器。本书仅讨论此类线性空间滤波器。

我们在3.4节讨论了空间卷积，并在公式(3-35)中实现了它，其中涉及到不同大小的函数的卷积。当我们使用DFT计算卷积定理中使用的变换时，这意味着我们正在对相同大小的周期函数进行卷积，如图4.27所示。因此，正如前面所解释的那样，公式(4-94)称为循环卷积(circular convolution)。

当计算速度、成本和大小是重要参数时，使用公式 (3-35) 进行空间卷积滤波非常适合使用硬件和/或固件的小核，如第 4.1 节所述。然而，当使用通用机器时，使用快速Fourier 变换 (FFT) 算法计算 DFT 的频域方法可以比使用空间卷积快数百倍，具体取决于所用核的大小，如图 4.2 所示。我们将在第 4.11 节讨论 FFT 及其计算优势。

滤波概念在频域中更直观，滤波器设计也通常更容易。利用这两个域的属性的一种方法是，在频域中指定一个滤波器，计算其离散Fourier变换 (IDFT)，然后使用所得的全尺寸空间核的属性作为构建较小核的指导。接下来将对此进行演示(请记住，Fourier变换及其逆变换是线性过程(参见问题 4.24)，因此讨论仅限于线性滤波)。在示例 4.15 中，我们演示了逆过程，即给定一个空间核，并得到其全尺寸频域表示。这种方法对于分析小空间核在频域中的行为非常有用。

频域滤波器可以作为指定我们在第3章中讨论过的一些小核的系数的指南。基于Gauss函数的滤波器尤其令人感兴趣，因为如表4.4所示，Gauss函数的正向和逆Fourier变换都是实Gauss函数。为了说明其基本原理，我们仅讨论一维空间。二维Gauss转移函数将在本章后面讨论。

令 H(u) 表示一维频域 Gauss 转移函数

(4-107) $\displaystyle H(u) = Ae^{-u^{2}/(2{\sigma}^{2})}$

其中，σ 是 Gauss 曲线的标准偏差。空间域中的核可通过取 H(u) 逆 DFT 而求得(见问题 4.48)：

(4-108) $h(x) = (\sqrt{2\pi}){\sigma}Ae^{-2{\pi}^{2} {\sigma}^{2} x^{2} }$

这两个公式之所以重要，有两个原因：(1) 它们是 Fourier 变换对，两个分量都是Gauss实数。这有利于分析，因为我们不必关心复数。此外，Gauss曲线直观易用。(2) 这两个函数具有互逆性。当 H(u) 具有较宽的轮廓( s 值较大 )时，h(x) 具有较窄的轮廓，反之亦然。事实上，当 s 趋向于无穷大时，H(u) 趋向于常数函数，而 h(x) 趋向于脉冲函数，这意味着在两个域中均不进行滤波。

图 4.36(a) 和 (b) 分别展示了Gauss低通滤波器的频域转移函数和空域对应函数的曲线图。假设我们想以图 4.36(b) 中 h(x) 的形状为指南，指定空域中小核函数的系数。图 4.36(b) 中函数的关键特性是其所有值均为正。因此，我们得出结论，可以使用系数均为正的核函数来实现空域低通滤波(就像我们在3.5节中所做的那样)。作为参考，图 4.36(b) 还展示了该节中讨论过的两个核函数。请注意，正如上一段所述，Gauss函数的宽度之间存在互易关系(reciprocal relationship)。频域函数越窄，低频衰减越厉害，从而导致图像模糊加剧。在空间域中，这意味着必须使用更大的核来增加模糊效果，如示例 3.11 所示。

----------------------------图 4.36：(a) 频域中的一维Gauss低通转移函数。(b) 空域中对应的核。(c)频域中的Gauss高通转移函数。(d)对应的核。图中所示的小型二维核是我们在第 3 章中使用的核。----------------------------------------------

正如你在 3.7 节中所知，我们可以用一个低通滤波器构造一个高通滤波器，方法是用一个常数减去一个低通函数。使用Gauss函数，我们可以通过使用所谓的Gauss差分(它涉及两个低通函数)来更好地控制滤波器函数的形状。在频域中，这变成

(4-109) $\displaystyle H(u) = Ae^{-u^{2}/(2{\sigma}_1^{2} )} - Be^{-u^{2}/(2{\sigma}_{2}^{2})}$

其中，A ≥ B 且 ${\sigma}_{1} > {\sigma}_{2}$ 。空间域中相应的函数是

(4-110) $\displaystyle h(x) = ({\sqrt{2\pi}}){\sigma}_{1}Ae^{-2{\pi}^{2} {\sigma}_{1}^{2} x^2} - ({\sqrt{2\pi}}){\sigma}_{2}Be^{-2{\pi}^{2} {\sigma}_{2}^{2} x^2}$

图 4.36 (c) 和 (d) 展示了这两个公式的图形。我们再次注意到宽度上的互易性，但这里最重要的特征是 h(x) 的中心项为正，两侧为负。图 4.36(d) 中所示的小核(我们在第三章中用它来进行锐化)“捕捉”了这一特性，从而说明了如何将频域滤波的知识用作选择空间核系数的基础。

尽管我们付出了巨大的努力才到达这里，但请放心，如果没有我们刚刚建立的基础，就不可能真正理解频域滤波。实际上，频域可以看作是一个“实验室”，我们可以在其中利用频率内容和图像外观之间的对应关系。正如本章后面将多次演示的那样，一些在空域中直接难以公式化的任务在频域中变得几乎微不足道。一旦我们通过频域实验选择了特定的滤波器转移函数，我们就可以选择使用FFT直接在该域中实现该滤波器，或者我们可以对该转移函数进行IDFT以获得等效的空域函数。如图4.36所示，一种方法是指定一个小的空间核，试图在空域中捕捉完整滤波器函数的“本质”。另一种更正式的方法是使用基于数学或统计标准的近似值来设计二维数字滤波器，正如我们在第 3.7 节中讨论的那样。

例子 4.15 从一个空间核获得一个频域转移函数

在本例中，我们从一个空间核开始，并展示如何在频域中生成其对应的滤波器转移函数。然后，我们比较了使用频域和空间技术获得的滤波结果。当希望将给定核与一个或多个“完整”滤波器候选在频域中的性能进行比较，或者希望深入了解核在空域中的性能时，这种分析非常有用。为简单起见，我们使用图 3.50(e) 中的 3 × 3 垂直 Sobel 核。图 4.37(a) 显示了我们希望滤波的 600 × 600 像素图像 f (x, y)，图 4.37(b) 显示了其频谱。

-------------------图 4.37：(a)建筑物图像；(b)其Fourier光谱。-----------------

图 4.38(a) 显示了 Sobel 核 h(x, y)(透视图解释如下)。由于输入图像的大小为 600 × 600 像素，而核的大小为 3 × 3，我们根据公式 (4-100) 和公式 (4-101)，用零填充 f 和 h 至 602 × 602 像素，以避免频域中的环绕误差。乍一看，Sobel 核似乎表现出奇对称性。然而，它的第一个元素并不像公式(4-81) 所要求的那样为 0。要将核转换为满足公式(4-83)的最小尺寸，我们必须在其中添加由 0 组成的前导行和前导列，这会将其变成一个大小为 4 × 4 的数组。我们可以将此数组嵌入到一个更大的全零数组中，如果该更大的数组是偶数维(就像 4 × 4 核一样)，并且它们的中心重合，则仍然可以保持其奇对称性，如示例 4.10 中所述。上述注释是滤波器生成的一个重要方面。如果我们在形成 $h_{p}(x,y)$ 时保留相对于填充数组的奇对称性，则根据表 4.1 中的性质 9 可知 H(u ,v) 将为纯虚数。正如我们在本示例末尾所示，这将产生与使用原始核 h(x ,y) 对图像进行空间滤波相同的结果。如果不保留对称性，结果将不再相同。

生成 H(u ,v) 的步骤如下：(1) 将 $h_{p}(x,y)$ 乘以 $(-1)^{x+y}$ 以使频域滤波器居中；(2) 对 (1) 式的结果进行正向离散Fourier变换 (DFT) 以生成 H(u ,v)；(3) 将 H(u ,v) 的实部设为 0，以解释寄生(parasitic)实部(我们知道 H 必须是纯虚数，因为 $h_{p}$ 为实数且为奇数)；(4) 将结果乘以 $(-1)^{x+y}$ 。最后一步反转了 H(u ,v) 乘以 $(-1)^{x+y}$ 的运算，当 h(x, y) 被手动置于 $h_{p}(x+y)$ 的中心时，此运算是隐式的。图 4.38(a) 显示了 H(u ,v) 的透视图，图 4.38(b) 显示了H (u ,v) 的图像。请注意，该图像关于其中心具有反对称性，这是由于 H(u ,v) 为奇数造成的。函数 H(u ,v) 可用作任何其他频域滤波器转移函数。图 4.38(c) 是使用刚刚获得的滤波器转移函数在频域中对图 4.37(a) 中的图像进行滤波的结果，采用了前面概述的逐步滤波程序。正如导数滤波器所预期的那样，边缘得到了增强，所有恒定强度区域都降为零(灰色调是由于显示缩放造成的)。图 4.38(d) 显示了使用 Sobel 核 h(x, y) 在空域中对同一图像进行滤波的结果，采用了第 3.6 节中讨论的程序。结果完全相同。

------------------图 4.38：(a) 空间核函数及其对应的频域滤波器转移函数的透视图。(b) 转移函数的图像。(c) 使用 (b) 中的转移函数对图 4.37(a) 进行频域滤波的结果。(d) 使用 (a) 中的核函数对同一幅图像进行空域滤波的结果。结果相同。--------------

4.8 使用低通频域滤波器进行图像平滑(Image smoothing using lowpass frequency domain filters)

本章的其余部分将讨论频域中的各种滤波技术，首先从低通滤波器开始。图像中的边缘和其他急剧的强度过渡(例如噪声)对其Fourier变换的高频成分有显著的影响。因此，在频域中，平滑(模糊)是通过高频衰减实现的，也就是通过低通滤波。在本节中，我们将考虑三种类型的低通滤波器：理想滤波器、Butterworth滤波器和Gauss滤波器。这三种滤波器涵盖了从非常急剧(理想)到非常平滑(Gauss)的滤波范围。Butterworth滤波器的形状由一个称为滤波器阶数(filter order)的参数控制。当该参数值较大时，Butterworth滤波器接近于理想滤波器。当该参数值较小时，Butterworth滤波器更像Gauss滤波器。因此，Butterworth 滤波器提供了两个“极值”之间的过渡。本节中的所有滤波均遵循上一节概述的程序，因此所有滤波器转移函数 H(u ,v) 的大小均理解为 P × Q ；也就是说，离散频率变量在 u = 0, 1, 2,…, P - 1 和 v = 0, 1, 2,…, Q - 1 范围内，其中 P 和 Q 分别是公式 (4-100) 和 (4-101) 给出的填充大小。

4.8.1 理想低通滤波器 (Ideal lowpass filters)

一个二维低通滤波器，如果能够使以原点为半径的圆内的所有频率无衰减地通过，并“截止”该圆外的所有频率，则称其为理想低通滤波器 (ILPF)；它由转移函数定义

(4-111) $H(u,v)=\left \{ \begin{array}{rl} 1 \hspace{0.2cm}\text{if}(D(u,v) \le D_{0}) \\ \\ 0 \hspace{0.2cm}\text{if}(D(u,v) > D_{0}) \end{array} \right .$

其中， $D_{0}$ 是一个正常量，且 D(u ,v) 是频率中的一个点 (u ,v)与 P × Q 频率矩形的中心之间的距离；即

(4-112) $D(u ,v) = [(u-P/2)^{2}+(v-Q/2)^{2} ]^{1/2}$

其中，如前一样，P 和 Q 是基于公式 (4-102) 和 (4-103) 的填充大小。图 4.39(a) 显示了转移函数 H (u ,v) 的透视图，图 4.39(b) 则以图像形式显示了该函数。如第 4.3 节所述，“理想”一词表示半径为 $D_{0}$ 的圆上或圆内的所有频率均会无衰减地通过，而圆外的所有频率则会完全衰减(被滤除)。理想低通滤波器的转移函数关于原点呈径向对称。这意味着它完全由径向横截面定义，如图 4.39(c) 所示。将横截面旋转 360° 即可获得该滤波器的二维表示。

----------------------------------图 4.39：(a) 理想低通滤波器转移函数的透视图。(b) 以图像形式显示的函数。(c)径向横截面。------------------------------------

对于 ILPF 的截面，H(u ,v) = 1 和 H(u ,v) = 0 之间的过渡点称为截止频率(cutoff frequency)。在图 4.39 中，截止频率为 $D_{0}$ 。ILPF 的尖锐截止频率无法用电子元件实现，尽管它们可以在计算机中模拟(但受限于像素间距离的限制，最快的过渡速度是有限的)。

本章通过研究低通滤波器在相同截止频率下的行为来比较它们。建立标准截止频率轨迹的一种方法是使用圆来包围指定数量的总图像功率(image power) $P_{T}$ 。 $P_{T}$ 可以通过将每一个点 (u ,v) 处的填充图像的功率谱分量相加来获得，其中 u = 0, 1, 2,…, P - 1 和 v = 0, 1, 2,…, Q - 1；即

(4-113) $\displaystyle P_{T} = \sum_{u=0}^{P-1}\sum_{v=0}^{Q-1}P(u,v)$

其中 P(u ,v) 由公式 (4-89) 给出。如果DFT 已居中，则以频率矩形中心为原点，半径为 $D_{0}$ 的圆包围 α 比例的功率，其中

(4-114) $\displaystyle \alpha = 100 \bigg [\sum_{u}\sum_{v}P(u,v)/P_{T} \bigg ]$

求和针对的是位于圆内或圆边界上的 (u ,v) 的值。

图 4.40 (a) 和 (b) 展示了测试图案图像及其频谱。叠加在频谱上的圆半径分别为 10、30、60、160 和 460 像素，分别包含图像标题中列出的总功率百分比。频谱迅速衰减，接近 87% 的总功率被半径为 10 的较小圆所包含。这一点的重要性将在以下示例中体现。

------------------------------图 4.40：(a) 688 × 688 像素大小的测试图案；(b)其频谱。由于填充，频谱是图像大小的两倍，但为了适应显示，仅显示一半大小。相对于全尺寸频谱，圆的半径分别为 10、30、60、160 和 460 像素。这些半径分别包含填充后图像功率的 86.9%、92.8%、95.1%、97.6% 和 99.4%。------------------------------

例子 4.16：使用低通滤波器在频域中对图像进行平滑处理。

图 4.41 展示了应用截止频率与图 4.40(b) 所示半径对应的 ILPF 的结果。图 4.41(b) 在所有实际应用中都毫无用处，除非模糊的目的是消除图像中的所有细节，除了代表最大物体的“斑点”。该图像的严重模糊清楚地表明，图像中大部分清晰的细节信息都包含在滤波器去除的 13% 功率中。随着滤波器半径的增加，去除的功率越来越少，从而导致模糊程度降低。请注意，图 4.41(c) 至 (e) 中的图像包含明显的“振铃效应”，随着去除的高频内容量减少，其纹理会变得更细腻。即使在仅去除了总功率 2% 的图像中 [图 4.41(e)]，振铃现象仍然可见。正如我们之前多次提到的那样，这种振铃行为是理想滤波器的一个特性。最后，图 4.41(f) 中 α = 99.4% 的结果显示出非常轻微的模糊和几乎难以察觉的振铃效应，但总体而言，该图像接近原始图像。这表明，在被 ILPF 滤除的频谱功率上方 0.6% 的部分中，几乎没有包含边缘信息。

--------------------------图 4.41：(a)原始图像尺寸为 688 × 688 像素。(b)-( f )使用截止频率分别为半径 10、30、60、160 和 460 的 ILPF 进行滤波的结果，如图 4.40(b)所示。这些滤波器消除的功率分别为总功率的 13.1%、7.2%、4.9%、2.4% 和 0.6%。我们使用图像填充来避免零填充带来的黑边问题，如图 4.31(c) 所示。----------------

ILPF 的模糊和振铃特性可以用卷积定理来解释。图 4.42(a) 展示了一个频域 ILPF 传递函数的图像，其半径为 15，大小为 1000 × 1000 像素。图 4.42(b) 是通过对 (a) 进行离散Fourier变换 (IDFT) 得到的 ILPF 的空间表示 h( x ,y)(注意振铃)。图 4.42(c) 展示了通过 (b) 中心的一条直线的强度分布。该分布类似于一个sinc 函数。(注：虽然该轮廓类似于正弦函数，但ILPF的变换实际上是一个Bessel函数，其推导超出了本文的讨论范围。需要记住的重点是，滤波器函数在频域中的“宽度”与波瓣在空间函数中的宽度“扩展”之间仍然保持反比关系。) 空间域滤波是通过将图 4.42(b) 中的函数与图像进行卷积来实现的。想象一下，图像中的每一个像素都是一个离散脉冲，其强度与该位置图像的强度成正比。将这个类正弦函数与一个脉冲进行卷积，会将函数复制(即，将原点移动)到脉冲的位置。也就是说，卷积会复制图 4.42(b) 中以图像中每一个像素位置为中心的函数。该空间函数的中心叶是造成模糊的主要原因，而外部较小的叶则主要造成振铃效应。由于空间函数的“扩散”与 H(u,v) 的半径成反比，因此 $D_{0}$ 越大(即通过的频率越多)，空间函数就越接近脉冲函数，在极限情况下，该脉冲函数与图像进行卷积时不会引起模糊。当 $D_{0}$ 变小时，情况正好相反。这种互易行为现在对你来说应该很熟悉了。在接下来的两节中，我们将展示，在实现模糊的同时，几乎没有或完全没有振铃效应，这是低通滤波的一个重要目标。

----------------------------4.42: (a) 频域 ILPF 传递函数。(b)相应的空域核函数。

(c)通过(b)图中心的水平线的强度分布。-----------------------------------------

4.8.2 Gauss低通滤波器 (Gauss lowpass filters)

Gauss 低通滤波器(GLPF)转移函数为

(4-115) $\displaystyle H(u ,v) = e^{-D^{2} (u,v)/(2{\sigma}^{2}) }$

其中，如公式 (4-112)中所述，D(u ,v) 是从 P × Q 频率矩形的中心到矩形内任意点 (u ,v) 的距离。与我们之前的Gauss函数表达式不同，这里我们不使用乘法常数，

以便与本节及后续章节讨论的滤波器保持一致，这些滤波器的最高值为 1 。与之前一样，σ 是关于中心的扩散度的度量。并令 $\sigma=D_{0}$ ，我们可以用与本节中其他函数相同的符号来表示Gauss转移函数：

(4-116) $\displaystyle H(u ,v) = e^{-D^{2} (u,v)/(2D_{0}^2)}$

其中， $D_{0}$ 是截止频率。当 $D(u ,v) = D_{0}$ ，其 GLPF 转移函数从其最大值 1.0 下降到 0.607 。从表 4.4 可知，频域 Gauss 函数的逆Fourier变换也服从Gauss函数。这意味着，通过计算公式 (4-115) 或 (4-116) 的 IDFT 得到的空间Gauss滤波器核不会产生振铃效应。如表 4.4 中的性质 13 所示，前面解释的 ILPF 的逆关系也适用于 GLPF。频域中较窄的Gauss传递函数意味着空域中较宽的核函数，反之亦然。图 4.43 展示了 GLPF 转移函数的透视图、图像显示和径向横截面。

------------------------------------图 4.43：(a) GLPF 传递函数的透视图。(b)以图像显示的函数。(c) 不同 $D_{0}$ 值的径向横截面。-----------------------------------

例子 4.17 使用 Gauss 低通滤波器在频域中对图像进行平滑。

图 4.44 显示了将公式 (4-116) 中的 GLPF 应用于图 4.44(a) 的结果，其中 $D_{0}$ 等于图 4.40(b) 中的五个半径。与使用 ILPF(图 4.41)获得的结果相比，我们注意到模糊效果随着截止频率的增加而平滑过渡。GLPF 的平滑度略低于 ILPF。关键区别在于使用 GLPF 时可以确保不会出现振铃效应。这在实践中是一个重要的考虑因素，尤其是在任何类型的伪影都不可接受的情况下，例如医学成像。如果需要更好地控制截止频率附近低频和高频之间的过渡，接下来讨论的Butterworth低通滤波器是更合适的选择。这种对滤波器轮廓的额外控制的代价是可能出现振铃效应，正如你稍后将看到的那样。

---------------------------------图 4.44：(a)原始图像，尺寸为 688 × 688 像素。(b)–(f) 使用截止频率为半径的 GLPF 进行滤波的结果，如图 4.40 所示。与图 4.41 进行比较。我们使用镜像填充来避免零填充带来的黑边。-----------------------------

4.8.3 Butterworth低通滤波器 (Butterworth lowpass filters)

一个其截止频率距离频率矩形的中心为 $D_{0}$ 的 n 阶Butterworth低通滤波器(BLPF)的转移函数是

(4-117) $H(u ,v) = \displaystyle \frac{1}{1+\left [D(u,v)/D_0 \right ]^{2n}}$

其中 D(u ,v) 由公式 (4-112) 给出。图 4.45 展示了 BLPF 函数的透视图、图像显示和径向横截面。比较图 4.39，4.43 和 4.45 中的横截面图，我们可以看出，可以控制 BLPF 函数，使其在使用较高 n 值时接近 ILPF 的特性，并在 n 值较低时接近 GLPF 的特性，同时提供从低频到高频的平滑过渡。因此，我们可以使用 BLPF 函数来接近 ILPF 函数的锐度，同时显著减少振铃。

---------------------------------图 4.45：(a) Butterworth低通滤波器转移函数的透视图。(b)以图像显示的函数。(c) 1 至 4 阶 BLPF 的径向截面。--------------------

例子 4.18：使用Butterworth低通滤波器对图像进行平滑处理。

图 4.46 (b)-( f ) 显示了将公式 (4-117) 中的 BLPF 应用于图 4.46(a) 的结果，其中截止频率等于图 4.40(b) 中的五个半径，n = 2.25。模糊程度介于使用 ILPF 和 GLPF 的结果之间。例如，将图 4.46(b) 与图 4.41(b) 和 4.44(b) 进行比较。BLPF 的模糊程度小于 ILPF，但大于 GLPF 。

----------------------图 4.46：(a)原始图像，尺寸为 688 × 688 像素。(b)–( f )使用截止频率为图 4.40 所示半径且 n = 2.25 的 BLPF 进行滤波的结果。与图 4.41 和 4.44 进行比较。我们使用镜像填充来避免零填充特有的黑边。------------------------

从一阶 BLPF 获得的空间域核没有振铃效应。通常，在二阶或三阶滤波器中，振铃效应难以察觉，但在更高阶的滤波器中可能会变得明显。图 4.47 比较了不同阶数 BLPF 对应的空间表示(即空间核)(所有情况下均使用 5 的截止频率)。图中还显示了沿穿过每个空间核中心的水平扫描线的强度分布。对应于一阶 BLPF 的核[见图 4.47(a)]既没有振铃效应也没有负值。对应于二阶 BLPF 的核确实出现了轻微的振铃效应和较小的负值，但这些负值肯定不如 ILPF 明显。如其余图像所示，对于高阶滤波器，振铃效应变得明显。 20 阶 BLPF 的空间核表现出与 ILPF 类似的振铃特性(在极限情况下，两个滤波器相同)。2 至 3 阶 BLPF 在有效的低通滤波和可接受的空间域振铃之间取得了良好的平衡。表 4.5 总结了本节讨论的低通滤波器转移函数。

-------------------------------图 4.47：(a)–(d) 对应于尺寸为 1000 × 1000 像素,截止频率为 5, 阶数分别为 1,2,5 和 20 的 BLPF 传递函数的空间表示(即空间核)。(e)–(h)通过滤波函数中心的相应强度分布。------------------------------------

表 4.5：低通滤波器转移函数。 $D_{0}$ 截止频率，n 是 Butterworth 滤波器的阶。

4.8.4 低通滤波器的其它例子 (Additional examples of lowpass filtering)

在下面的讨论中，我们将展示低通滤波在频域中的几种实际应用。第一个例子来自机器感知领域，应用于字符识别；第二个例子来自印刷出版行业；第三个例子与处理卫星和航拍图像有关。使用3.5节中讨论的低通空间滤波技术可以获得类似的结果。为了保持一致性，我们在所有示例中都使用了GLPF，但使用BLPF也可以获得类似的结果。请记住，图像会被填充到两倍大小以进行滤波，如公式(4-102)和(4-103)所示，滤波器转移函数必须与填充后的图像大小匹配。以下示例中使用的 $D_{0}$ 值反映了这种两倍的滤波器大小。

图 4.48 显示了一段低分辨率文本样本。例如，在传真、复印件和历史记录中，我们经常会遇到类似的文本。这个样本没有其他问题，例如污迹、折痕和撕裂。图 4.48(a) 中的放大部分显示，由于分辨率不足，该文档中的字符形状扭曲，并且许多字符都已损坏。虽然人类可以轻松地通过视觉填补这些空白，但机器识别系统在读取损坏的字符时却遇到了很大困难。解决此问题的一种方法是通过模糊输入图像来弥补其中的细小空白。图 4.48(b) 显示了使用的Gauss低通滤波器，通过这个简单的过程可以很好地“修复”字符。通常，在完成上述“修复”之后，还会进行其他处理，例如阈值处理和细化，以获得更清晰的字符。我们将在第 9 章讨论细化，并在第 10 章讨论阈值化。

-----------------------图 4.48：(a)低分辨率示例文本(请注意放大视图中的破损字符)。(b)使用 GLPF 滤波的结果，显示破损字符中的间隙已被连接。-----------------------

低通滤波是印刷出版行业的重要组成部分，它用于许多预处理功能，包括3.6节中讨论的反锐化掩模。“美化(cosmetic)”处理是低通滤波在印刷前的另一种用途。图4.49展示了低通滤波的一种应用，它使清晰的原图呈现出更平滑、更柔和的效果。对于人脸，通常的目标是降低细小皮肤线条和小瑕疵的锐度。图4.49(b)和(c)中的放大部分清晰地显示了拍摄对象眼周细小皮肤线条的显著减少。事实上，平滑后的图像看起来非常柔和，赏心悦目。

----------------------------图 4.49：(a)原始图像尺寸为 785 × 732。(b) 使用 $D_{0}=150$ 的 GLPF 进行滤波的结果。(c) 使用 $D_{0}=130$ 的 GLPF 进行滤波的结果。请注意 (b) 和 (c) 中放大部分中细纹的减少。---------------------------------------

图 4.50 展示了在同一幅图像上应用低通滤波的两种方法，但目的截然不同。图 4.50(a) 是一幅 808 × 754 的甚高分辨率辐射计 (VHRR) 图像片段，显示了墨西哥湾(暗)和佛罗里达(亮)的部分区域(请注意水平传感器扫描线)。水体之间的边界是由环路电流引起的。这幅图像说明了遥感图像中，传感器倾向于沿着场景扫描方向产生明显的扫描线。 (参见示例 4.24 中 a b c 图 4.50 (a) 808 × 754 卫星图像，显示明显的水平扫描线。(b) 使用 $D_{0}=50$ 的 GLPF 滤波结果。(c) 使用 $D_{0}=20$ 的 GLPF 滤波结果。(原图由 NOAA 提供。) 低通滤波是一种粗略(但简单)的降低这些线影响的方法，如图 4.50(b) 所示(我们将在 4.10 和 5.4 节中讨论更有效的方法)。此图像是使用 $D_{0}=50$ 的 GLFP 获得的。平滑图像中扫描线影响的降低可以简化宏观特征的检测，例如洋流之间的界面边界。

图 4.50(c) 显示了 $D_{0} = 20$ 时，采用明显更激进的Gauss低通滤波的结果。此处，目标是尽可能地模糊细节，同时保留较大的特征。例如，这种滤波可以作为图像分析系统的预处理阶段的一部分，用于在图像库中搜索特征。此类特征的一个例子可以是给定大小的湖泊，例如佛罗里达州东南部地区的奥基乔比湖，如图 4.50(c) 所示，它是一个近乎圆形的暗区，周围环绕着一个较亮的区域。低通滤波有助于通过平均化比目标特征更小的特征来简化分析。

-------------------------------图 4.50：(a) 808 × 754 卫星图像，显示明显的水平扫描线。(b) 使用 $D_{0}=50$ 的 GLPF 进行滤波的结果。(c) 使用 $D_{0} = 20$ 的 GLPF 进行滤波的结果。(原始图像由 NOAA 提供。)---------------------------------------

4.9 使用高通频域滤波器锐化图像(Image sharpening using highpass filters)

我们在上一节中表明，可以通过衰减Fourier变换的高频分量来平滑图像。由于边缘和其他强度的突变与高频分量相关，因此可以通过高通滤波在频域中实现图像锐化，高通滤波可以衰减低频分量，而不会干扰Fourier变换中的高频分量。与第 4.8 节一样，我们只考虑径向对称的零相移滤波器。本节中的所有滤波均基于第 4.7 节中概述的程序，因此假设所有图像都被填充为 P × Q 大小[参见公式 (4-102)和 (4-103)]，并且滤波器转移函数 H(u ,v) 被理解为大小为 P × Q 的中心离散函数。

4.9.1 基于低通滤波器的理想、Gauss、和Butterworth 高通滤波器(Ideal, gaussian, and butterworth highpass filters from lowpass filters)

与空间域中的核的情况一样(参见第 3.7 节)，从 1 中减去低通滤波器的转移函数，可得到频域中相应的高通滤波器的转移函数：

(4-118) $H_{HP}( u ,v ) = 1 - H_{LP}( u ,v )$

其中， $H_{LP}$ 是低通滤波器转移函数。因此，根据公式 (4-111) 可推出理想高通滤波器的转移函数(IHPF)为:

(4-119) $H(u,v)=\left \{ \begin{array}{rl} 0 \hspace{0.2cm}\text{if}(D(u,v) \le D_{0}) \\ \\ 1 \hspace{0.2cm}\text{if}(D(u,v) > D_{0}) \end{array} \right .$

其中，照例，D(u , v) 是到 P × Q 频域矩形中心的距离，类似地，基于公式 (4-116) 可推导出 Gauss 高通滤波器(GHPF)的转移函数为

(4-120) $\displaystyle H(u,v) = 1 - e^{-D^{2} (u,v)/(2D_0^{2})}$

而基于公式 (4-117) Butterworth 高通滤波器(BHPF)的转移函数为

(4-121) $H(u ,v) = \displaystyle \frac{1}{1+[D_{0}/D(u,v)]^{2n}}$

图 4.51 展示了上述转移函数的三维图、图像表示和径向截面。与之前一样，我们看到图中第三行的 BHPF 转移函数表示 IHPF 转移函数的锐度与 GHPF 转移函数的宽广平滑度之间的过渡。

---------------------------图4.51: 顶行：IHPF 转移函数的透视图、图像和径向横截面。中间行和底行：GHPF 和 BHPF 转移函数的相同序列。(为了清晰起见，添加了细图像边框。它们不属于数据。)-----------------------------------------------------

根据公式 (4-118) 可推出，频域中高通滤波器转移函数对应的空间核为

(4-122) $\begin{array}{rl} h_{HP}( x ,y )&=\frak{F}^{-1}[H_{HP}( u ,v )] \\ \\ &=\frak{F}^{-1}[1-H_{LP}( u ,v )] \\ \\ &=\delta(x,y)-h_{LP}(x,y) \end{array}$

其中，我们利用了频域中 1 的离散Fourier变换 (IDFT) 在空域中为单位脉冲这一事实(见表 4.4)。该公式正是 3.7 节讨论的基础，在该节中，我们展示了如何通过从单位脉冲中减去低通核来构造高通核。

图 4.52 展示了以此方式构建的高通空间核，使用公式 (4-122)，并采用 ILPF、GLPF 和 BLPF 转移函数(本图中使用的 M、N 和 $D_{0}$ 值与图 4.42 中的相同，BLPF 的阶数为 2）。图 4.52(a) 展示了使用公式 (4-122) 得到的理想高通核，图 4.52(b) 是通过核中心的水平强度分布图。分布图的中心元素是一个单位脉冲，在图 4.52(a) 的中心可见为一个亮点。请注意，此高通核具有与图 4.42(b) 中对应的低通核相同的振铃特性。您很快就会看到，振铃和以前一样令人讨厌，但这次是在使用理想高通滤波器锐化的图像中。图 4.52 中的其他图像和轮廓分别针对Gauss核和Butterworth核。从图 4.51 可知，频域中的 GHPF 转移函数往往比大小和截止频率相当的Butterworth函数具有更宽的“裙边”。因此，我们预期Butterworth空间核比同类 Gauss 核“更宽”，图 4.52 中的图像及其轮廓证实了这一事实。表 4.6 总结了前几段讨论的三个高通滤波器转移函数。

--------------------------------图 4.52：(a)–(c)：由 IHPF、GHPF 和 BHPF 频域转移函数获得的理想、Gauss和Butterworth高通空间核(细图像边界不包含在数据中。)(d)–(f)：通过核中心的水平强度分布。--------------------------------------

表 4.6：高通滤波器转移函数。 $D_{0}$ 是截止频率，n 是Butterworth转移函数的阶数。

例子 4.19:字符测试模式的高通滤波。

图 4.53 的第一行显示了使用 IHPF、GHPF 和 BHPF 传递函数 ( $D_{0}=60$ [参见图 4.37(b)] 和巴特沃斯滤波器的 n = 2）对图 4.37(a) 中的测试图案进行滤波的结果。我们从第三章知道，高通滤波会产生带有负值的图像。图 4.53 中的图像未缩放，因此负值被显示器截断为 0(黑色)。高通滤波的主要目的是锐化。此外，由于此处使用的高通滤波器将直流项设置为零，因此图像基本上没有色调，正如前面结合图 4.30 所解释的那样。

本例的主要目标是比较三个高通滤波器的性能。如图 4.53(a) 所示，理想的高通滤波器会产生由振铃效应引起的严重失真。例如，大字母“a”笔画内的斑点就是振铃效应造成的。相比之下，图 4.53(b) 和 (c) 都没有出现此类失真。参考图 4.37(b)，滤波器去除或衰减了大约 95% 的图像能量。众所周知，去除图像的低频部分会显著减少其灰度内容，主要留下边缘和其他尖锐的过渡部分，如图 4.53 所示。你在图的第一行看到的细节仅包含在图像能量的前 5% 中。

第二行是在 $D_{0}=160$ 时获得的，结果更有趣。这些图像的剩余能量约为第一行图像能量的 2.5%，即一半。然而，精细细节的差异却令人瞩目。例如，看看大字母“a”的边界现在有多清晰，尤其是在Gauss和Butterworth结果中。所有其他细节，甚至最小的物体，也都呈现出同样的效果。当边缘和边界检测很重要时，这种结果视为是可以接受的。

---------------------------------图 4.53：第一行：图 4.40(a) 中的图像，经 IHPF、GHPF 和 BHPF 转移函数滤波，所有情况下 $D_{0}=60$ (BHPF 中 n = 2）。第二行：相同序列，但 $D_{0}=160$ 。-----------------------------------------------------------

图 4.54 显示了图 4.53 第二行的图像，并使用公式 (2-31) 和 (2-32) 进行了缩放，以显示正负强度的完整强度范围。图 4.54(a) 中的振铃效应表明了理想高通滤波器的不足。相比之下，请注意另外两幅图像的背景平滑度和边缘的清晰度。

---------------------------------图 4.54：图 4.53 第二行的图像使用公式 (2-31) 和 (2-32) 进行缩放，以显示正值和负值。----------------------------------------

例子 4.20: 使用高通滤波和阈值来增强图像。

图 4.55(a) 是一张 962 × 1026 的指纹图像，其中明显存在污迹(一个典型问题)。

自动指纹识别的一个关键步骤是增强指纹脊线并减少污迹。在本例中，我们使用高通滤波来增强脊线并减少污迹的影响。脊线的增强是通过其边界以高频为特征来实现的，而高通滤波器不会改变这些高频成分。另一方面，滤波器会降低低频成分，这些低频成分对应于图像中缓慢变化的强度，例如背景和污迹。因此，增强是通过降低除高频特征(在本例中我们感兴趣的特征)之外的所有特征的影响来实现的。

图 4.55(b) 是使用截止频率为 50 的 4 阶Butterworth高通滤波器的结果。四阶滤波器提供从低频到高频的尖锐(但平滑)过渡，其滤波特性介于理想滤波器和Gauss滤波器之间。所选的截止频率约为图像长边的 5%。这样做的目的是使 $D_{0}$ 接近原点，以便低频被衰减但不会完全消除(直流项除外，直流项设置为 0)，这样脊线和背景之间的色调差异就不会完全消失。将 $D_{0}$ 的值设置为图像长边的 5% 到 10% 之间是一个良好的开端。选择较大的 $D_{0}$ 值会突出精细细节，以至于会影响脊线的清晰度。正如预期的那样，高通滤波后的图像具有负值，在显示屏上显示为黑色。

在高通滤波图像中，一种突出尖锐特征的简单方法是设置阈值，将所有负值设为黑色 (0)，其余值设为白色 (1)。图 4.55(c) 展示了此操作的结果。请注意，脊线变得清晰，污迹的影响也显著减弱了。事实上，图 4.55(a) 中图像右上部分几乎不可见的脊线在图 4.55(c) 中得到了很好的增强。自动算法会发现，在这幅图像上跟踪脊线比在原始图像上更容易。

--------------------------------图 4.55：(a)模糊的指纹。(b)高通滤波结果(a)。(c)阈值处理结果(b)。(原图由美国国家标准与技术研究院提供。)----------------------

4.9.2 频域中的Laplace算子(The laplacian in the frequency domain)

在 3.6 节中，我们使用Laplace算子在空域中进行图像锐化。在本节中，我们将重新讨论Laplace算子，并证明它使用频域技术可以得到等效的结果。可以证明(参见问题 4.52)，Laplace算子可以在频域中使用滤波器转移函数

(4-123) $H(u, v) = -4{\pi}^{2}(u^{2} + v^{2})$

来实现，或者，对于频率矩形的中心，使用转移函数

(4-124) $\begin{array}{rl} H(u, v) &= -4{\pi}^2[(u-P/2)^{2} + (v-Q/2)^{2}] \\ \\ &= -4{\pi}^{2} D^{2}(u, v) \end{array}$

来实现，其中，是 D(u, v) 是公式 (4-112) 中定义的距离函数。使用这个转移函数，一幅图像 f (x,y) 的 Laplace 算子可按相似的方式获得：

(4-125) $\nabla^{2}f(x,y)=\frak{F}^{-1}[H(u,v)F(u,v)]$

其中，F(u, v ) 是 f (x,y ) 的 DFT 。与公式 (3-54) 一样，增强是通过公式

(4-126) $g (x,y) = f (x,y) + c\nabla^{2} f (x,y)$

实现的，在这里，c = -1 ，因为 H(u, v) 是负的。在第 3 章中，f (x,y) 和 $\nabla^{2} f (x,y)$ 都具有可比值。然而，使用公式 (4-125) 计算 $\nabla^{2} f (x,y)$ 引入了 DFT 的缩放因子，此缩放因子这可能比 f 的最大值大几个数量级。因此，f 与其 Laplace 算子之间的差值必须被放到可比的范围内。处理这个问题最简单的方法是将 f (x,y) 的值归一化到 [0, 1] 范围内(在计算其 DFT 之前)，然后用其最大值去除以 $\nabla^{2} f (x,y)$ ，这样就可以将其近似放到 [-1, 1] 范围内。(谨记，Laplace算子有负值。) 然后就可以使用公式 (4-126)。

我们可以将公式 (4-126) 直接写成频域形式

(4-127) $\begin{array}{rl} g(x,y)&\displaystyle=\frak{F}^{-1}\{ F(u, v) - H(u, v)F(u, v)\} \\ \\ &\displaystyle=\frak{F}^{-1}\{ F(u, v)(1-H(u, v))\} \\ \\ &\displaystyle=\frak{F}^{-1}\{ [1 + 4{\pi}^{2} D^{2}(u, v)]F(u, v) \} \end{array}$

虽然这个结果很优雅，但它也存在刚才提到的缩放问题，而且归一化因子的计算也比较困难。因此，在频域中，公式 (4-126) 是首选实现，其中 $\nabla^{2} f (x,y)$ 使用公式 (4-125) 计算，并使用上一段中提到的方法进行缩放。

例子 4.21：使用 Laplace 算子在频域中锐化图像。

图 4.56(a) 与图 3.46(a) 相同，图 4.56(b) 显示了使用公式 (4-126) 的结果，其中Laplace算子是使用公式 (4-125) 在频域中计算的。缩放操作按照公式 (4-126) 中的说明进行。通过比较图 4.56(b) 和 3.46(d)，我们可以看出频域结果更佳。图 4.56(b) 中的图像更加清晰，并且显示了 3.46(d) 中几乎不可见的细节，而 3.46(d) 是使用图 3.45(b) 中的Laplace核获得的，中心为 -8。在频域中取得的显著改进并不意外。空间Laplace核包含非常小的邻域，而公式(4-125) 和 (4-126) 涵盖了整个图像。

--------------------------------图 4.56：(a)原始模糊图像。(b)使用频域Laplace算子增强的图像。与图 3.46(d) 比较。(原始图像由 NASA 提供。)--------------------

4.9.3 反锐化掩模、高频增强滤波和高频增强滤波(Unsharp masking, high-boost filtering, and highfrequency-emphasis filtering)

在本节中，我们讨论第 3.6 节中介绍的反锐化掩模和高提升滤波图像锐化技术的频域公式。使用频域方法，公式 (3-55) 中定义的掩模可表示为

(4-128) $g_\text{mask}(x,y) = f (x,y) - f_{LP}(x,y)$

且

(4-129) $f_{LP}(x,y) = \frak{F}^{-1}[H_{LP}( u ,v )F(u ,v)]$

其中， $H_{LP}( u ,v )$ 是一个低通滤波器转移函数，而 F(u ,v) 是 f (x,y) 的 DFT 。在此， $f_{LP}(x,y)$ 是一个类似于公式 (3-55) 中的 $\overline{f}(x,y)$ 那样的平滑图像。则如公式 (3-56) 所述，

(4-130) $g(x,y) = f (x,y) + kg_\text{mask}(x,y)$

该表达式定义了当 k = 1 时的非锐化掩模和当 k > 1 时的高提升滤波。利用前面的结果，我们可以将公式 (4-130) 完全表示为涉及低通滤波器的频域计算：

(4-131) $g(x,y) = \frak{F}^{-1}\{[1 + k(1 - H_{LP}( u ,v ))]F(u ,v)\}$

我们可以使用公式(4-118)以高通滤波器的形式来表达这个结果：

(4-132) $g(x,y) = \frak{F}^{-1}\{[1 + kH_{HP}( u ,v )]F(u ,v)\}$

方括号内的表达式称为高频增强滤波器转移函数(high-frequencyemphasis

filter transfer function)。如前所述，高通滤波器将直流项设置为零，从而将滤波图像的平均强度降低到0。高频增强滤波器不存在这个问题，因为它在高通滤波器传递函数中添加了1 。常数k控制影响最终结果的高频比例。高频增强滤波的一个更通用的公式是

(4-133) $g(x,y) = \frak{F}^{-1}\{[k_{1} + k_{2} H_{HP}( u ,v )]F(u ,v)\}$

其中， $k_{1} \ge 0$ 偏移传递函数的值，以免将直流项归零[见图 4.30(c)]， $k_{2} > 0$ 控制高频的贡献。

例子 4.22 使用高频增强滤波器增强图像。

图 4.57(a) 显示了一幅 503 × 720 像素的胸部 X 光图像，其强度范围较窄。本例的目标是使用高频增强滤波来增强图像。X 光无法像光学透镜那样聚焦，因此生成的图像通常会略微模糊。由于这幅特定图像的强度偏向灰阶的暗端，我们也借此机会举例说明如何使用空域处理来补充频域滤波。

在医学图像处理中，诸如振铃效应之类的图像伪影是不可接受的，因此我们使用Gauss高通滤波器传递函数。由于Gauss高通滤波器传递函数的空间表示也是Gauss的，因此我们知道振铃效应不会成为问题。所选的 $D_{0}$ 值应提供足够的滤波来锐化边界，同时又不会过度锐化微小细节(例如噪声)。我们使用了 $D_{0}=70$ ，大约是长图像尺寸的 10%，但其他类似的值也可以。图 4.57(b) 是对原始图像(按图 4.54 中的图像缩放)进行高通滤波的结果。正如预期的那样，图像特征相当平淡，但重要的边界(例如肋骨的边缘)清晰可见。图 4.57(c) 显示了高频强调滤波的优势，其中我们使用了公式(4-133)，其中 $k_{1}=0.5$ , $k_{2}=0.75$ 。虽然图像仍然较暗，但灰度色调已得到恢复，并且特征更加清晰。

正如我们在第 3.3 节中讨论的那样，灰度级分布在较窄范围内的图像是直方图均衡化的理想选择。如图 4.57(d) 所示，这确实是一种进一步增强图像的合适方法。请注意骨骼结构和其他细节的清晰度，而这些细节在其他三幅图像中根本无法看到。最终的增强图像略有噪点，但这是 X 射线图像在灰度扩展时的典型特征。高频增强和直方图均衡化相结合的效果优于单独使用其中任何一种方法的效果。

-----------------图 4.57：(a) 胸部 X 光片。(b) 使用 GHPF 函数滤波的结果。(c) 使用相同 GHPF 函数进行高频增强滤波的结果。(d) 对 (c) 图像进行直方图均衡化的结果。(原始图像由密歇根大学医学院解剖科学系 Thomas R. Gest 博士提供。)----------

4.9.4 同态滤波(Homomorphic filtering)

2.3 节中介绍的照度-反射(illumination-reflectance)模型可用于开发一种频域程序，通过同时进行强度范围压缩和对比度增强来改善图像的外观。根据该节的讨论，图像 f (x, y) 可以表示为其光照 i(x, y)分量和反射率分量 r(x, y)分量之积：

(4-134) $f (x, y) =i(x, y) r(x, y)$

该公式不能直接用于计算照明和反射的频率分量，因为乘积的 Fourier 变换不是变换的乘积：

(4-135) $\displaystyle \frak{F}[f (x, y)] \leq \frak{F}[i(x, y)] \frak{F}[ r(x, y)]$

然而，假设我们定义

(4-136) $\begin{array}{rl} z (x, y) &= \ln f (x, y) \\ \\ &=\ln i(x, y) + \ln r(x, y) \end{array}$

则

(4-137) $\frak{F}[z (x, y)] = \frak{F}[\ln i(x, y)] + \frak{F}[\ln r(x, y)]$

或

(4-138) $Z(u ,v)=F_{i}(u ,v) + F_{r}(u ,v)$

其中， $F_{i}(u ,v)$ 和 $F_{r}(u ,v)$ 分别是 $\ln i(x, y)$ 和 $\ln r(x, y)$ 的 Fourier 变换。我们可以使用滤波器转移函数 H(u ,v) 来滤波 Z(u ,v)，使得

(4-139) $\begin{array}{rl} S(u ,v) &= H(u ,v) Z(u ,v) \\ \\ &=H(u ,v) F_{i}(u ,v) + H(u ,v)F_{r}(u ,v) \end{array}$

则空间域的滤波图像是

(4-140) $s (x, y) = \frak{F}^{-1}[H(u ,v) F_{i}(u ,v)] + \frak{F}^{-1}[H(u ,v) F_{r}(u ,v)]$

通过定义

(4-141) $i^{'} (x, y) = \frak{F}^{-1}[H(u ,v) F_{i}(u ,v)]$

和

(4-142) $r^{' }(x, y) = \frak{F}^{-1}[H(u ,v) F_{r}(u ,v)]$

们可以将(4-140)表示为

(4-143) $s (x, y) = i^{'} (x, y) + r^{'} (x, y)$

最后，因为 z(x, y) 是通过对输入图像取自然对数形成的，我们通过对滤波结果取指数来逆转该过程以形成输出图像：

(4-144) $\begin{array}{rl} g(x,y)&\displaystyle=e^{s(x,y)}\\ \\ &\displaystyle=e^{i^{'}(x,y)}e^{r^{'}(x,y)}\\ \\ &\displaystyle=i_{0}(x,y)r_{0}(x,y) \end{array}$

其中，

(4-145) $i_{0}(x, y) = e^{i^{'} (x,y)}$

和

(4-146) $r_{0}(x, y) = e^{r^{'} (x,y)}$

是输出(处理后)图像的照度和反射分量。

图 4.58 是对刚刚推导的滤波方法的总结。该方法基于一类称为同态系统的特例。在这个特定的应用中，该方法的关键在于将照度分量和反射分量分离，如公式 (4-138) 所示。同态滤波器转移函数 H(u ,v) 可以分别对这些分量进行运算，如公式 (4-139) 所示。

-------------------------------图 4.58：同态滤波步骤总结。----------------------------------

图像的照度分量通常具有缓慢的空间变化特征，而反射分量则倾向于发生突变，尤其是在不同物体的交界处。这些特性导致我们将图像对数Fourier变换的低频与光照关联起来，将高频与反射关联起来。虽然这些关联只是粗略的近似，但它们可以在图像滤波中得到充分利用，如示例 4.23 所示。

使用同态滤波器可以对照度和反射分量进行很好的控制。这种控制需要指定一个滤波器传递函数 H(u ,v)，该函数以不同的、可控的方式影响Fourier变换的低频和高频分量。图 4.59 显示了该函数的横截面。如果选择参数 $\gamma_{L}$ 和 $\gamma_{H}$ ，使得 ${\gamma}_{L} < 1$ 且 ${\gamma}_{H} \ge 1$ ，则图 4.59 中的滤波器函数将衰减低频(照度)的贡献，并放大高频(反射)的贡献。最终结果是同时压缩动态范围并增强对比度。

图 4.59 中的函数形状可以用高通滤波器转移函数来近似。例如，使用稍微修改过的 GHPF 函数形式，可以得到同态函数

(4-147) $H(u ,v) = (\gamma_{H} - \gamma_{L})[1 - e^{-cD^{2} (u,v)/D_{0}^{2} }] + \gamma_{L}$

其中，D(u ,v) 定义于公式 (4-112)，而常数 c 控制函数在 $\gamma_{L}$ 和 $\gamma_{H}$ 之间转换时斜率的锐度。该滤波器转移函数类似于上一节讨论的高频增强函数。

-----------------------图 4.59：同态滤波器转移函数的径向截面。------------------------------

例子 4.23：同态滤波。

图 4.60(a) 显示了 1162 × 746 像素的全身 PET(正电子发射断层扫描)扫描图像。图像略显模糊，许多低强度特征被占据显示屏动态范围的“热点”的高强度信号所掩盖。(这些热点是由脑肿瘤和肺肿瘤引起的。)图 4.60(b) 是通过对图 4.60(a) 进行同态滤波获得的，使用公式 (4-147) 中的滤波器转移函数，其中 ${\gamma}_{L} = 0.4 , {\gamma}_{H} = 3.0, c = 5 , D_{0} = 20$ 。该函数的径向横截面与图 4.59 非常相似，但斜率更陡峭，低频和高频之间的过渡更接近原点。

--------------------------图 4.60：(a)全身 PET 扫描。(b)使用同态滤波增强的图像。(原图由 CTI Pet Systems 的 Michael E. Casey 博士提供。)--------------------------

4.10 选择性滤波(Selective filtering)

前两节讨论的滤波器作用于整个频率矩形。有些应用需要处理特定的频带或频率矩形的小区域。第一类滤波器称为带通滤波器。如果频带中的频率被滤除，则该滤波器称为带阻滤波器；同样，如果允许这些频率通过，则该滤波器称为带通滤波器。第二类滤波器称为陷波滤波器(notch filters)。这些滤波器进一步被定义为陷波抑制滤波器或陷波通滤波器，具体取决于陷波区域中的频率是被抑制还是通过。

4.10.1 带阻和带通滤波器(Bandreject and bandpass filters)

正如你在 3.7 节中了解到的那样，频域中的带通和带阻滤波器转移函数可以通过组合低通和高通滤波器转移函数来构建，而高通滤波器转移函数也可以从低通函数推导出来(参见图 3.52)。换言之，低通滤波器转移函数是构建高通、带阻和带通滤波器函数的基础。此外，带通滤波器转移函数可以通过带阻函数获得，其方式与通过低通转移函数获得高通相同：

(4-148) $H_{BP}(u , v) =1 - H_{BR}(u , v)$

图 4.61(a) 展示了如何构建理想带阻滤波器(IBRF)的转移函数。它由截止频率不同的 ILPF 和 IHPF 函数组成。处理带通函数时，我们感兴趣的参数是带宽 W 和带中心 $C_{0}$ 。

通过观察图 4.61(a)，可以轻松得到 IBRF 函数的公式，如表 4.7 最左侧所示。带通转移函数的关键要求是：(1) 函数值必须在 [0, 1] 范围内；(2) 函数值在距函数原点(中心) $C_{0}$ 处必须为零；(3) 我们必须能够指定 W 的值。显然，刚刚开发的 IBRF 函数满足这些要求。

将低通和高通转移函数相加，形成Gauss和Butterworth带阻函数，这会带来一些困难。例如，图 4.61(b) 显示了由具有不同截止点的低通和高通Gauss函数之和构成的带通函数。两个问题显而易见：我们无法直接控制 W，并且 H(u ,v) 的值在 $C_{0}$ 处不为 0 。我们可以对函数进行偏移和缩放，使其值落在 [0, 1] 范围内，但无法找到低通和高通 Gauss 函数相交点的解析解，而要求解 $C_{0}$ 的截止点，就需要该交点。唯一的选择是反复试验或数值方法。

幸运的是，除了添加低通和高通转移函数外，还有一种替代方法是修改Gauss和Butterworth高通传递函数的表达式，使其满足前面提到的三个要求。我们举例说明Gauss函数的计算过程。在本例中，我们首先将公式 (4-120) 中 H(u ,v) = 0 处的点从 D(u ,v) = 0 改为 $D(u ,v) = C_{0}$ :

(4-149) $H(u ,v) = 1 - \displaystyle e^{\displaystyle -\left [\frac{(D(u,v)-C_{0})^{2}}{W^{2} }\right ] }$

该函数的图(图 4.61(c))表明，在 $C_{0}$ 以下，该函数表现为低通 Gauss 函数；在 $C_{0}$ 处，该函数始终为 0；而对于高于 $C_{0}$ 的值，该函数表现为高通Gauss函数。参数 W 与标准差成正比，因此控制着波段的“宽度”。唯一的问题是，该函数在原点处并不总是为 1。对公式 (4-149) 进行简单修改即可消除这一缺陷：

(4-150) $H(u ,v) = 1 - \displaystyle e^{\displaystyle -\left [ \frac{D^{2}(u,v)-C_{0}^{2}}{D(u,v)W} \right ]^{2}}$

现在，当 D(u ,v) = 0 时，指数为无穷大，这使得指数项趋于零，并且 H(u ,v) = 1，正如所期望的那样。在公式 (4-149) 的这个修改中，基本的Gauss形状得以保留，并且满足了前面提到的三个要求。图 4.61(d) 显示了公式 (4-150) 的图形。类似的分析可以得到Butterworth带阻滤波器转移函数的形式，如表 4.7 所示。

---------------图 4.61：径向截面。(a) 理想带阻滤波器转移函数。(b) 由Gauss低通和高通滤波器函数之和构成的带阻传递函数。(最小值不为 0，且与 $C_{0}$ 不对齐。) (c)公式 (4-149) 的径向图。(最小值是 0，且与 $C_{0}$ 正确对齐，但原点处的值不为 1。) (d)公式 (4-150) 的径向图；此Gauss形状的图满足带阻滤波器转移函数的所有要求。------

表 4.7 ：带阻滤波器的转移函数。 $C_{0}$ 为频带中心，W 为频带宽度，D(u,v) 为传递函数中心到频率矩形中点 (u,v) 的距离。

图 4.62 展示了刚刚讨论过的滤波器转移函数的透视图。乍一看，Gauss函数和Butterworth函数似乎大致相同，但与之前一样，Butterworth函数的行为介于理想函数和Gauss函数之间。如图 4.63 所示，通过将三个滤波器函数视为图像更容易看出这一点。提高Butterworth函数的阶数会使其更接近理想的带阻转移函数。

-----------------------------图 4.62：表 4.7 中 (a) 理想带阻滤波器 (b) 改进Gauss带阻滤波器和 (c) 改进Butterworth带阻滤波器(1 阶)转移函数的透视图。所有转移函数元素大小均为 512 × 512，其中 $C_{0} = 128$ ，W = 60。-----------------------------

-----------------------------图 4.63：图 4.62 中的 (a) 理想带通转移函数、(b) Gauss带通转移函数和(c)Butterworth带通转移函数，以图像形式显示。(细边框线不属于图像数据。)---------------------------------------------------------------

4.10.2 陷波滤波器(Notch filters)

陷波滤波器是最常用的选择性滤波器。陷波滤波器抑制(或允许)频率矩形预定义邻域内的频率。零相移滤波器必须围绕原点(频率矩形的中心)对称，因此，中心位于 $(u_{0},v_{0})$ 的陷波滤波器转移函数在 $(-u_{0},-v_{0})$ 位置处必定有一个对应的陷波。陷波抑制滤波器的转移函数由高通滤波器转移函数的乘积构成，这些高通滤波器的转移函数的中心已平移到陷波的中心。其一般形式为：

(4-151) $\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})$ 。这些中心针对频率矩形中心 (M/2 , N/2)指定，其中，照例，M 和 N 分别是输入图像的行和列。因此，每一个滤波器转移函数的距离计算公式分别为

(4-152) $D_{k}(u,v) = [(u-M/2-u_{k} )^{2}+(v-N/2-v_{k} )^2 ]^{1/2}$

和

(4-153) $D_{-k}(u,v) = [(u-M/2+u_{k} )^{2}+(v-N/2+v_{k} )^{2} ]^{1/2}$

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

(4-154) $\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_{k}(u,v)$ 和 $D_{-k}(u,v)$ 分别由公式 (4-152) 和 (4-152) 给出。常量 $D_{0k}$ 对于每一个陷波对是相同的，但对于不同的陷波对可以是不同的。其他陷波抑制滤波器函数的构造方式相同，具体取决于所选的高通滤波器函数。与前面讨论的滤波器一样，陷波带通滤波器的转移函数可以通过陷波抑制函数使用以下表达式获得：

(4-155) $H_{NP}(u,v) = 1 - H_{NR}(u,v)$

正如接下来的两个示例所示，陷波滤波的主要应用之一是选择性地修改DFT的局部区域。通常，这种处理是以交互方式进行的，直接与未使用填充的DFT进行交互。与实际DFT交互操作(而不是必须将填充后的频率值“转换”为实际的频率值)的优势通常大于滤波过程中不使用填充可能导致的任何环绕误差。如有必要，在获得可接受的解决方案后，可以通过调整所有滤波器参数来补偿填充后的DFT大小，从而生成使用填充的最终结果。以下两个示例均未使用填充。要了解DFT值如何随填充而变化，请参见问题4.42。

例子 4.24：使用陷波滤波从数字化印刷媒体图像中消除moiré条纹(pattern)。

图 4.64(a) 是图 4.21 中使用的扫描报纸图像，其中显示出明显的moiré条纹；图 4.64(b) 是其频谱。纯正弦波(一个周期函数)的Fourier变换是一对共轭对称脉冲(见表 4.4)。图 4.64(b) 中对称的“脉冲状”突发脉冲是由于moiré条纹的近似周期性造成的。我们可以利用陷波滤波来衰减这些突发脉冲。

图 4.64(c) 显示了将图 4.64(a) 的 DFT 与Butterworth陷波抑制转移函数相乘的结果，其中 $D_{0}=9$ ，n = 4，适用于所有陷波对(陷波的中心与图中黑色圆形区域的中心重合）。半径值的选取(通过目视检查光谱)以完全覆盖能量爆发，n 值的选取以产生具有急剧过渡的陷波。陷波中心的位置是根据光谱交互确定的。图 4.64(d) 显示了使用此滤波器转移函数获得的结果，使用了第 4.7 节中概述的滤波程序。考虑到原始图像的低分辨率和退化程度，这种改进是显著的。

-------------------图 4.64：(a)采样的报纸图像，显示出moiré条纹。(b) 频谱。(c) Fourier变换乘以Butterworth陷波抑制滤波器的转移函数。(d)滤波后的图像。---------

例子 4.25：使用陷波滤波消除周期性干扰。

图 4.65(a) 展示了环绕土星的部分光环。这幅图像是由Cassini号——第一艘进入土星轨道的航天器——拍摄的。图像中可见的近正弦波图案，是由在数字化图像之前叠加在相机视频信号上的交流信号引起的。这是一个意料之外的问题，损坏了部分任务图像。幸运的是，这种干扰很容易通过后处理进行校正。一种方法是使用陷波滤波。

图 4.65(b) 展示了 DFT 频谱。仔细分析纵轴，可以发现原点附近有一系列小的能量突发(burst)，这些能量突发对应于近正弦波的干扰。一种简单的方法是使用一个窄陷波矩形滤波器，从最低频率的突发开始，一直延伸到纵轴的其余部分。图 4.65(c) 展示了此类滤波器的转移函数(白色代表 1，黑色代表 0)。图 4.65(d) 显示了使用此滤波器处理损坏图像的结果。此结果与原始图像相比有了显著的改善。

----------------------------------图 4.65：(a)土星环图像，显示近乎周期性的干涉。(b)光谱。(原点附近垂直轴上的能量爆发对应于干涉图样。) (c)垂直陷波抑制滤波器的转和多函数。(d) 滤波结果。 ((c)中的细黑边框不是数据的一部分。)(原图由NASA/JPL 的 Robert A.West 博士提供。)-------------------------------------------------

为了获得仅包含干涉图案的图像，我们使用陷波传递函数(通过从 1 中减去陷波抑制函数获得)分离了纵轴上的频率[参见图 4.66(a)]。然后，如图 4.66(b) 所示，滤波后图像的 IDFT 就是空间干涉图案。

-----------------------------图 4.66：(a)用于隔离图 4.65(a) 中 DFT 的垂直轴的陷波滤波器函数。(b)通过计算 (a) 中的 IDFT 获得的空间模式。-------------------

4.11 快速Fourier变换(The fast Fourier transform)

到目前为止，我们主要关注频域滤波的理论概念和示例。现在应该清楚的是，图像处理领域的计算要求并不低。因此，对简化和加速Fourier变换计算的方法有一个基本的了解至关重要。本节将讨论这些问题。

4.11.1 二维DFT的可分离性(Separability of the 2-D DFT)

如表 4.3 中所示，二维DFT 可分离成一维变换。我们可以将 (4-67) 记为

(4-156) $\begin{array}{rl} F(u, v) &\displaystyle= \sum_{x=0}^{M-1}e^{-j2{\pi}ux/M} \sum_{y=0}^{N-1}f(x,y)e^{-j2{\pi}vy/N} \\ \\ &\displaystyle=\sum_{x=0}^{M-1}F(x,v)e^{-j2{\pi}ux/M} \end{array}$

其中，

(4-157) $\displaystyle F(x, v) = \sum_{y=0}^{N-1}f(x,y)e^{-j2{\pi}vy/N}$

对于 x 的一个值，并且 v = 0, 1, 2,…, N - 1，我们看到 F(x, v) 是 f (x, y) 中某一行的一维离散傅里叶变换 (DFT)。通过在公式 (4-157) 中将 x 从 0 变为 M - 1，我们计算 f (x, y) 所有行的一组一维离散Fourier变换 (DFT)。类似地，公式 (4-156) 中的计算是 F(x, v) 各列的一维变换。因此，我们得出结论，f (x, y) 的二维离散Fourier变换 (DFT) 可以通过计算 f (x, y) 每一行的一维变换，然后沿结果的每一列计算一维变换来获得。这是一个重要的简化，因为我们一次只需处理一个变量。类似的推导也适用于使用一维离散 Fourier 变换 (IDFT) 计算二维离散Fourier变换 (IDFT)。然而，正如我们在下一节中所示，我们可以使用设计用于计算前向 DFT 的算法来计算 IDFT，因此所有二维Fourier变换计算都简化为多次执行设计用于计算一维 DFT 的一维算法。

4.11.2 使用 DFT 算法计算 IDFT (Computing the IDFT using a DFT algorithm)

对公式 (4-68) 两边取复共轭，并将结果乘以 MN ，可得

(4-158) $\displaystyle MNf^{*}(x, y) = \sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F^{*} (u,v)e^{-j2{\pi}(ux/M+vy/N)}$

但我们将结果的右侧形式视为 $F^{*}(u,v)$ 的 DFT 。因此，公式 (4-158) 表明，若我们将 $F^{*}(u,v)$ 代入分配给计算2维正向Fourier变换的算法，则结果将为 $MNf^{*}(x, y)$ 。取复共轭并将结果除以 MN 就得到 f (x, y) ，其是 F(u, v) 的逆。

基于一维变换的连续转移(如上一节所述)的二维正向DFT算法计算二维逆变换经常引起混淆，这涉及到复共轭和乘以常数，而这两者都不是一维算法中完成的。需要注意的关键概念是，我们只需将 $F^{*}(u,v)$ 输入我们所有的任意正向算法。结果将为 $MNf^{*}(x, y)$ 。要得到 f (x, y)，我们只需取其复共轭，然后除以常数 MN 。当然，当 f (x, y) 为实数时(通常如此)， $f^{*}(x, y) = f (x, y)$ 。

4.11.3 快速 Fourier 变换 (The fast Fourier transform)( FFT)

如果我们必须直接实现方程 (4-67) 和 (4-68)，那么在频域中进行工作是不切实际的。这些方程的强力实现需要大约 $MN^{2}$ 次的乘法和加法。对于中等大小的图像(例如 2048 × 2048 像素)，这意味着仅进行一次二维 DFT 运算就需要进行大约 17 万亿次乘法和加法运算(不包括指数运算，这些指数运算可以计算一次并存储在查找表中)。如果没有快速Fourier变换 (FFT) 的发现，其将计算量减少到 $MN\log_{2}MN$ 次乘法和加法的量级，可以肯定地说，本章介绍的内容几乎没有实际价值。FFT 带来的计算量减少确实令人印象深刻。例如，计算一幅 2048 × 2048 图像的二维 FFT 运算需要大约 9200 万次乘法和加法运算，这比上面提到的 1 万亿次计算量级大幅减少。

尽管FFT是信号处理文献中广泛讨论的一个主题，但该主题在我们的工作中至关重要，如果我们不提供介绍来解释FFT的工作原理，本章将是不完整的。我们为实现此目标而选择的算法是所谓的逐次加倍法，它是促成整个行业诞生的原始算法。该特定算法假设样本数是2的整数幂，但这并不是其他方法的普遍要求(Brigham[1988])。从上一节我们知道，二维DFT可以通过逐次进行一维变换来实现，因此我们只需关注单个变量的FFT。

在FFT的推导中，通常将公式(4-44)表示为

(4-159) $\displaystyle F(u) = \sum_{x=0}^{M-1}f(x)W_{M}^{ux} (u = 0 , 1 ,2 ,..., M - 1)$

其中，

(4-160) $W_{M}^{ux} = e^{-j2{\pi}/M}$

且 M 假设为具有形式

(4-161) $M = 2^{p}$

其中， p 是正整数。则可推导出 M 可表示为

(4-162) M = 2K

其中，K 也为正整数。将公式 (4-162) 代入公式 (4-159)得到：

(4-163) $\begin{array}{rl} F(u) &\displaystyle= \sum_{x=0}^{2K-1}f(x)W_{2K}^{ux} \\ \\ &\displaystyle=\sum_{x=0}^{K-1}f(2x)W_{2K}^{u(2x)} + \sum_{x=0}^{K-1}f(2x+1)W_{2K}^{u(2x+1)} \end{array}$

然而，可以使用公式 (4-160) 证明 $W_{2K}^{2ux} = W_{K}^{ux}$ ，因此，公式 (4-163) 可以记为

(4-164) $\displaystyle F(u) = \sum_{x=0}^{K-1}f(2x)W_{K}^{ux} + \sum_{x=0}^{K-1}f(2x+1)W_{K}^{ux} W_{2K}^{u}$

定义

(4-165) $\displaystyle F_\text{even}(u) = \sum_{x=0}^{K-1}f(2x)W_{K}^{ux}(u = 0 , 1 ,2 ,... , K - 1)$

和

(4-166) $\displaystyle F_\text{odd}(u) = \sum_{x=0}^{K-1}f(2x+1)W_{K}^{ux} (u = 0 , 1 ,2,... , K - 1)$

则 (4-164) 简化为

(4-167) $F(u) = F_\text{even}(u) + F_\text{odd}(u)W_{K}^{ux}$

此外，由于 $W_{M}^{u+K} = W_{K}^{u}$ 和 $W_{2K}^{u+K} = - W_{2K}^{u}$ ，则可推导出

(4-168) $F(u + K ) = F_\text{even}(u) - F_\text{odd}(u)W_{2K}^{u}$

对公式 (4-165) 至 (4-168) 的分析揭示了这些表达式的一些重要(且令人惊讶)的属性。M 点 DFT 可以通过将原表达式分为两部分来计算，如公式 (4-167) 和 (4-168) 所示。计算 F(u) 的前半部分需要计算公式 (4-165) 和 (4-166) 中给出的两个 (M /2) 点变换。然后将 $F_\text{even}(u)$ 和 $F_\text{odd}(u)$ 的结果代入公式 (4-167) 中，即可得到当 u = 0 , 1 ,2 ,… , (M/2 – 1) 时的 F(u)。另一半则直接由公式 (4-168) 得出，无需进行额外的变换计算。

研究上述过程的计算含义十分有趣。令 𝔪(p) 和 𝔞(p) 分别表示实现该方法所需的复数乘法和加法的次数。与前面一样，样本数为 2p，其中 p 为正整数。首先假设 p = 1，即样本数为 2。二点变换需要计算 F(0)；然后根据公式 (4-168) 可知 F(1)。要获得 F(0)，需要计算 $F_\text{even}(0)$ 和 $F_\text{odd}(0)$ 。在这种情况下，K = 1，公式 (4-165) 和 (4-166) 是单点变换。然而，由于单个样本点的 DFT 就是样本本身，因此无需乘法或加法即可获得 $F_\text{even}(0)$ 和 $F_\text{odd}(0)$ 。将 $F_\text{even}(0)$ 与 $W_{2}^{0}$ 相乘，然后进行一次加法，即可根据公式 (4-167) 得到 F(0)。然后，F(1) 根据公式 (4-168) 再进行一次加法运算(减法视为加法相同)。由于 $F_\text{odd}(0)W_{2}^{0}$ 已经计算完毕，因此两点变换所需的运算总数为 𝔪(1) = 1 次乘法运算和 𝔞(1) = 2 次加法运算。

p 的下一个允许值为 2。根据前面的推导，四点变换可以分为两部分。F(u) 的前半部分需要计算两个两点变换，如公式 (4-165) 和公式 (4-166) 中 K = 2 所示。两点变换需要 𝔪(1) 次乘法和 𝔞(1) 次加法。因此，计算这两个公式总共需要 2 𝔪(1) 次乘法和 2 𝔞(1) 次加法。为了根据公式 (4-167) 得到 F(0) 和 F(1)，还需要进行两次乘法和两次加法。由于已经针对 u = {0,1} 计算了 $F_\text{odd}(0) W_{2K}^{u}$ ，因此再进行两次加法得到 F(2) 和 F(3)。因此，总和为 𝔪(2) = 2 𝔪(1) + 2 和 𝔞(2) = 2 𝔞(1) + 4 。

当 p 等于 3 时，需要进行两次四点变换来计算 $F_\text{even}(u)$ 和 $F_\text{odd}(u)$ 。这两次变换需要 2 𝔪(2) 次乘法和 2 𝔞(2) 次加法。再进行四次乘法和八次加法即可完成完整的变换。因此，总和为 𝔪(3) = 2 𝔪(2) + 4 次乘法和 𝔞(3) = 2 𝔞(2) + 8 次加法。

对任意正整数 p 继续进行此论证，可得到实现 FFT 所需乘法和加法次数的递归表达式：

(4-169) $\frak{m}(p) = 2\frak{m}(p - 1 ) + 2^{p-1} (p \ge 1)$

和

(4-170) $\frak{a}(p) = 2\frak{a} (p - 1 ) + 2^{p} ( p \ge 1 )$

其中，𝔪(0) = 0 和 𝔞(0) = 0 ，因为单点变换不要求任乘法或加法。

刚刚开发的方法称为逐次加倍 FFT 算法，因为它基于从两个单点变换计算两点变换，从两个两点变换计算四点变换等等(对于任何等于 2 的整数幂的 M )，以一公式的证明留作习题(见问题 4.63)

(4-171) $\displaystyle \frak{m}(p) = \frac{1}{2}M\log_{2}M$

和

(4-172) $\frak{a}(n) = M\log_{2}M$

其中， $M = 2^{p}$ 。

FFT 相对于直接实现 1维 DFT 的计算优势定义为

(4-173) $\displaystyle C(M) = \frac{M^{2}}{M\log_{2}M}= \frac{M}{\log_{2}M}$

其中， $M^{2}$ 是“强力”实现一维DFT所需的运算次数。因为假设 $M=2^{p}$ ，因此我们可以根据 p 来书写公式 (4-173) ：

(4-174) $C(p) = \displaystyle \frac{2^{p}}{p}$

该函数的图(图 4.67)显示，计算优势随 p 的增加而迅速增加。例如，当 p = 15(32,768 个点)时，FFT 的计算优势几乎是 DFT 的强力实现的 2,200 倍。因此，我们预计在同一台机器上，FFT 的计算速度几乎是 DFT 的 2,200 倍。正如您在 4.1 节中了解到的那样，FFT 也比空间滤波具有显著的计算优势，这两种方法的交叉点在于相对较小的核。