ImageMagick 示例 -- 傅里叶变换

ImageMagick 示例：前言与索引

简介

傅里叶变换

• 图像就是波

• 图像中的二维波

  ImageMagick 中的 FFT/IFT

• 使用 FFT

• 仅振幅或仅相位图像
• FFT 频谱图像
• HDRI FFT 图像

• 作为实部/虚部成分的 FFT

  傅里叶变换的性质

• 常量图像的 FFT

• DC 颜色的影响
• 正弦波图像的频谱
• 直接生成 FFT 图像
• 矩形图案图像的频谱
• 平面圆形图案图像的频谱
• 高斯圆形图案图像的频谱

• 网格图案图像的频谱

  实际应用

• 改变图像对比度 - 系数开方

• 模糊图像 - 低通滤波
• 检测图像边缘 - 高通滤波
• 锐化图像 - 高提升滤波

• 去除噪声 - 陷波滤波

  高级应用

  FFT 乘法与除法（低层示例 - 子页面）

简介

图像处理中最难理解的概念之一就是傅里叶变换。原因有两个。首先，它在数学上比较高级；其次，变换得到的图像不像原图，因此很难解释。不过，利用傅里叶变换，可以用新的方式完成一些熟悉的处理，例如增强亮度和对比度、模糊、锐化以及去除噪声。它还可以提供一些普通图像域中做不到的能力，包括对运动模糊、镜头失焦等典型相机畸变进行反卷积（也称去模糊），以及使用归一化互相关进行图像匹配。本页的目标，是尝试解释傅里叶变换的背景和简化数学，并给出使用傅里叶变换可以完成的处理示例。如果你觉得内容太多，可以跳过这部分，直接从 ImageMagick 中的 FFT/IFT 开始，关注性质和示例。感兴趣的读者还可以看看另一篇简单易懂、包含光学类比的讨论：An Intuitive Explanation of Fourier Theory。Vanderbilt University School Of Engineering 的讲义对数学兴趣更强的人也很有帮助：1 & 2 Dimensional Fourier Transforms 和 Frequency Filtering。其他数学参考资料包括 Wikipedia 上的 Fourier Transform、Discrete Fourier Transform、Fast Fourier Transform，以及 Complex Numbers。感谢 Sean Burke 编写最初的演示代码，也感谢 ImageMagick 的创建者将它整合进 ImageMagick。这两件事都很了不起。 许多示例会使用 HDRI 版 ImageMagick，以便保持变换图像的精度。如果你想充分利用这些技术，建议编译一个个人使用的 HDRI 版本。

傅里叶变换

图像通常由一组“像素”数组构成，每个像素由一组值定义：红、绿、蓝，有时还包括透明度。不过这里我们会忽略透明度。因此，红、绿、蓝这三个“通道”各自包含一组“强度”或“灰度”值。这称为“空间域”中的栅格图像。这只是一个更正式的说法，意思是图像由它在每个“位置”或“空间位置”上的“强度值”定义。不过，图像也可以用另一种方式表示，称为图像的“频率域”。在这个域中，每个图像通道都用正弦波来表示。在这样的“频率域”中，每个通道都有“振幅”值，这些值不是按 X、Y“空间”坐标存储，而是按 X、Y“频率”存储。由于这是数字表示，频率是某个“最小”或单位频率的倍数，而像素坐标表示这个单位频率的索引或整数倍。这来自一个基本原理：“任何表现良好的函数都可以表示为正弦波的叠加（组合或求和）”。换句话说，“频率域”表示只是存储并重现“空间域”图像的另一种方式。那么，图像怎样才能表示为“波”呢？

图像就是波

如果我们从 任意 图像中取出一行或一列像素，并把它画成图（这里用脚本「[im_profile](../static/img/scripts/im_profile)」通过 "gnuplot" 生成），你会发现它看起来很像一条波。

  magick holocaust_tn.gif -colorspace gray miff:- |\
    im_profile -s - image_profile.gif

如果这些起伏在间距和振幅上更规则，就会得到更像波形图案的东西，例如……

  magick -size 20x150 gradient: -rotate 90 \
          -function sinusoid 3.5,0,.4   wave.gif
  im_profile -s wave.gif wave_profile.gif

不过，虽然这个规则波形与上面的图像剖面有些相似，它仍然太规则了。如果把更多波相加，就可以得到更接近原图剖面的图案。

  magick -size 1x150 gradient: -rotate 90 \
          -function sinusoid 3.5,0,.25,.25     wave_1.png
  magick -size 1x150 gradient: -rotate 90 \
          -function sinusoid 1.5,-90,.13,.15   wave_2.png
  magick -size 1x150 gradient: -rotate 90 \
          -function sinusoid 0.6,-90,.07,.1    wave_3.png

  magick wave_1.png wave_2.png wave_3.png \
          -evaluate-sequence add added_waves.png

另一个替代示例可参见添加带偏置的渐变。这种“波的叠加”（波相加）已经接近得多，但仍不能精确匹配图像图案。不过，你可以继续这样做，添加更多波并调整它们，使最终合成波越来越接近原始图像的实际剖面。最终，只要加入足够多的波，就可以精确重现图像的原始剖面。这正是数学家 Joseph Fourier 的发现。它的现代表述是：“任何表现良好的函数都可以表示为正弦波的叠加”。换句话说，只要把足够多、频率和振幅恰到好处的正弦波相加，就可以重现任何起伏变化的图案。因此，“频率域”表示只是存储并重现“空间域”图像的另一种方式。“傅里叶变换”就是找出一幅图像由哪些“波”组成的过程，就像上面的例子所做的那样。

图像中的二维波

上面展示了如何用多条正弦波近似图像单行剖面的一个例子。不过图像是二维的，因此在“频率域”中表示图像所用的波也需要是二维的。下面是这样一个二维波的例子。这个波包含多个组成部分。 图像示例

在 ImageMagick 中使用 FFT/IFT

实现说明

ImageMagick 使用 FFTW，离散傅里叶变换库，它要求图像与浮点值（复数）之间相互转换，并最早在 IM 6.5.4-3 中实现。为了让它按人们通常对图像的期待工作，任何非正方形图像，或某个维度为奇数的图像，都会被填充（使用虚拟像素）为以图像最大宽度或高度为边长的正方形。为了让“FFT 原点”能正确居中在图像中心，它还会被强制为偶数（2 的倍数）尺寸。这样做的结果是，应用逆傅里叶变换之后，图像需要裁剪回原始尺寸，以移除填充。由于傅里叶变换由“复数”构成，变换结果无法直接可视化。因此，复数变换会以下面两种形式之一，被分离成两幅成分图像。
复数
实部/虚部

实部和虚部

“复数”通常的数学和数值表示，是一对浮点值，分别由“实部”(a) 和“虚部”(b) 组成。不幸的是，这两个数字可能包含负值，因此无法形成可查看的图像。因此，这种表示不能用于普通版本的 ImageMagick，因为那样会裁剪这些图像（下面会看到由此产生的效果示例）。不过，使用 HDRI 版 ImageMagick 时，你仍然可以生成、使用，甚至保存这种傅里叶变换图像的表示。它们本身作为图像可能没什么用，甚至不能直接查看，但你仍可对其应用许多数学操作。要生成这种表示，需要使用运算符的“加号”形式：「[+fft](https://usage.imagemagick.org/fourier/option_link.cgi#fft)」和「[+ift](https://usage.imagemagick.org/fourier/option_link.cgi#ift)」，下面的作为实部/虚部成分的 FFT 会详细介绍。
复数极坐标
振幅/相位

振幅和相位

“复数”的直接数值表示，对图像工作并不太有用。但如果把这些值绘制到二维平面上，就可以把数值变成由“振幅”(r) 和“相位”(θ) 成分组成的极坐标表示。这种形式在图像处理中非常有用，尤其是振幅成分，它本质上指定了构成图像的所有频率。“振幅”成分只包含正值，并会直接映射为图像值。它没有固定的取值范围，不过除 DC 或零频率颜色外，这些值通常都相当小。因此，振幅图像通常会显得非常暗（几乎是黑色）。通常需要缩放振幅，并对其强度值应用对数变换，才能显现出视觉细节。得到的“对数变换”振幅图像称为图像的“频谱”。不过请记住，应该用于逆变换的是“振幅”图像，而不是“频谱”图像。出现在图像中心“原点”的 DC（“Direct Current”的缩写）或“零频率”颜色，是整幅图像的平均颜色值。另外，由于输入图像不包含“虚部”成分，DC 相位值也总是零相位，从而产生纯灰色。“相位”成分的范围则是 -π 到 +π。它首先被偏置到 0 到 2π 的范围，然后缩放到 0 到 QuantumRange 的实际图像值范围内（由编译时内存质量决定）。因此，零相位会得到纯灰值（对每个通道都合适），而取反的相位会得到纯黑（「0」）值。注意，纯白（「_QuantumRange_」）几乎但并不完全等同。使用普通 FFT 运算符「[+fft](https://usage.imagemagick.org/fourier/option_link.cgi#fft)」和「[+ift](https://usage.imagemagick.org/fourier/option_link.cgi#ift)」可以生成图像的振幅和相位 FFT 表示。下面会先在生成 FFT 图像及其逆变换中介绍。

生成 FFT 图像及其逆变换

（振幅和相位）

现在，我们简单地对 Lena 图像尝试一次傅里叶变换往返。也就是说，先做正变换，然后立即应用逆变换，得到原图。接着比较结果，看看产生的质量水平。

  time magick lena.png -fft -ift lena_roundtrip.png


  echo -n "RMSE = "

  magick compare -metric RMSE lena.png lena_roundtrip.png null:

  echo -n "PAE = "

  magick compare -metric PAE  lena.png lena_roundtrip.png null:

上面的「[compare](basics.html#compare)」程序会返回两幅图像差异程度的度量。在这个例子中，可以看到总体差异很小，大约是 0.22%。至少一个像素上的峰值差异（「PAE」，Peak Absolute Error）也只有约 1%。使用 HDRI 版 ImageMagick 可以改进这一点。（见下方 HDRI FFT）。我们再仔细看看上面这次往返中生成的 FFT 图像。

  magick lena.png -fft    +depth +adjoin lena_fft_%d.png

原图 | |
振幅 |
相位
---|---|---|---

正如 John M. Brayer 谈到傅里叶变换时所说…… 我们通常不显示 PHASE 图像，因为大多数看到它们的人不久后都会屈服于致幻剂，或最终进入一座西藏寺院。 注意，「[-fft](https://usage.imagemagick.org/fourier/option_link.cgi#fft)」运算符生成了两幅图像，第一幅是“振幅”成分（是的，它大多是黑色，中间只有一个彩色点），第二幅几乎像随机噪声的图像则包含“相位”成分。PNG 图像每个文件只能存一幅图像，因此「[+adjoin](https://usage.imagemagick.org/fourier/option_link.cgi#adjoin)」或输出文件名中的「%d」实际上并不是必需的，IM 会自行处理。不过为了完整性，我在上面保留了这些选项，以明确表示我生成的是两个单独的图像文件，而不是一个。更多细节见写出多图像序列。由于生成了两幅图像，振幅图像（第一幅，也就是第零幅）保存为「lena_fft_0.png」，相位图像（第二幅）保存为「lena_fft_1.png」。 | _为了避免保存 FFT 图像导致任何失真，最好完全不要把它们保存到磁盘，而是在处理图像时把它们保留在内存中。

如果必须保存，最好使用 Magick 文件格式「[MIFF](files.html#miff)」，以保留图像的最高质量（位深）。这种格式也可以在一个文件中保存多幅图像。做脚本工作时，也可以使用详细的「[TXT](files.html#txt)」枚举像素格式。

不要 使用「[JPEG](formats.html#jpg)」或「[GIF](formats.html#gif)」图像格式保存它们。

如果你必须为了实际查看而把这些图像保存成文件，例如给 Web 浏览器查看，请使用「[PNG](formats.html#png)」图像格式，并把「[+depth](https://usage.imagemagick.org/fourier/option_link.cgi#depth)」重置为内部默认值（就像这些示例中所做的）。不过它每个文件只能存一幅图像。

「[TIFF](formats.html#tiff)」文件格式也可以使用，虽然对 Web 浏览器来说不太合适，但它确实允许一个文件中包含多幅图像。

_
---|---
把中间图像保存到单个文件中的最佳方法，是使用「[MIFF](files.html#miff)」文件格式……

  magick lena.png -fft  +depth lena_fft.miff

或者，也可以使用「[-write](https://usage.imagemagick.org/fourier/option_link.cgi#write)」把它们保存成完全独立的文件名（见写出图像）……

  magick lena.png -fft  +depth \
          \( -clone 0 -write lena_magnitude.png +delete \) \
          \( -clone 1 -write lena_phase.png +delete \) \
          null:

注意，上面我使用了特殊的「[NULL:](files.html#null)」图像格式来丢弃两幅图像，而它们仍然保留在内存中，可用于进一步处理。最后，我们再次读入这两幅图像，把它变回普通的“空间”图像……

  magick lena_magnitude.png lena_phase.png -ift lena_restored.png

FFT 过程生成的两幅图像都对修改非常敏感，即使很小的变化也可能导致结果严重失真。因此，绝不要用任何可能扭曲这些值的图像格式保存它们，这一点很重要。还要记住，从频率域恢复图像时，两幅图像都是必需的。所以如果计划用它们重建图像，只保存一幅而丢弃另一幅是没有用的。

仅振幅或仅相位图像

最后，我们试着只从振幅成分或只从相位成分重建图像。

  magick lena_fft_0.png  -size 128x128 xc:'gray(50%)' \
                                                  -ift lena_magitude_only.png

  magick -size 128x128 xc:gray1  lena_fft_1.png  -ift lena_phase_only.png

仅振幅 |
仅相位
---|---

从这里可以注意到，真正包含图像大部分位置信息的是相位图像，而振幅实际上保存了很多颜色信息。这并不绝对，因为信息之间有些重叠，但大体如此。“仅振幅”图像总是会有白色角落，因为这里使用的是恒定 50% 相位图像。可以使用随机化相位图像来去掉这些白色斑块。不过请确保中心像素的相位是完美的 50% 灰色，否则整幅图像都会变暗。“仅相位”图像使用了恒定 1% 灰色（几乎纯黑）的振幅图像来转换。即便振幅恒定，它仍会产生一些强度很高的像素斑块，尤其是在边缘附近。你只需要记住，要重建原图，两幅图像都是必需的。

频率频谱图像

你应该已经注意到，振幅图像（第一幅或第零幅图像）看起来几乎全黑。它其实并不是真的全黑，只是对我们的眼睛来说所有值都非常非常小。这样的图像拿来研究并没有太多视觉上的趣味，因此我们用对数变换增强结果，生成一幅“频率频谱”图像。做法是对归一化后的“振幅”图像应用很强的 Evaluate Log Transform。

  magick lena_fft_0.png -auto-level -evaluate log 10000 \
          lena_spectrum.png

现在我们可以在振幅图像的频谱版本中看到细节。你甚至可能在频谱图像中看到一些特定颜色，但通常这些颜色对频谱图像并不重要。更重要的是每个频率的整体强度，以及它们形成的图案。因此，增强后你也可能希望把频谱图像转为灰度。需要使用多强的对数增强取决于图像，所以应不断调整，直到得到足够细节，能清楚看出图像频率频谱的图案。 也可以使用下面这个小 shell 脚本，为特定振幅图像计算一个可用的 对数缩放因子。 |

  scale=`magick lena_fft_0.png -auto-level \
          -format "%[fx:exp(log(mean)/log(0.5))]" info:`
  magick lena_fft_0.png -auto-level \
          -evaluate log $scale    lena_spectrum_auto.png

不过请记住，不能把频谱图像用于逆「[-ift](https://usage.imagemagick.org/fourier/option_link.cgi#ift)」变换，因为它会产生过亮的图像。 |

  magick lena_spectrum.png lena_fft_1.png -ift lena_roundtrip_fail.png

基本上，由于你增强了“振幅”图像，也会以相同方式增强结果图像，从而产生图中那种严重“裁剪”的结果。

HDRI FFT 图像

当我们把傅里叶变换结果映射成图像表示时，会把浮点“复数”值缩放并转换为整数图像值。这自然会产生舍入误差以及其他“Quantum”效应，尤其是在较小的低频振幅中。如果图像处理中的精度很重要，那么要么需要使用较高的位质量（例如 Q32 或 Q64 位版本的 ImageMagick），要么最好使用 HDRI 版 ImageMagick，使数值以浮点数保存。使用 HDRI 版 IM，并采用傅里叶变换的振幅和相位表示时，振幅成分仍然全部为正值，因此仍可按上面所示使用，只是精确得多。相位成分仍会像前面所示那样被偏置和缩放。换句话说，HDRI 中的振幅和相位表示完全相同，只是精度高得多。
例如，这里我使用 HDRI 版 ImageMagick 生成另一次图像“往返”转换。 |

  # HDRI version of IM used
  time magick lena.png -fft -ift lena_roundtrip_hdri.png


  echo -n "RMSE = "

  magick compare -metric RMSE lena.png lena_roundtrip_hdri.png null:

  echo -n "PAE = "

  magick compare -metric PAE  lena.png lena_roundtrip_hdri.png null:

如果把上面的结果与之前的非 HDRI 比较结果进行对比……

作为实部/虚部成分的 FFT

到目前为止，我们只看了傅里叶变换图像的“振幅”和“相位”表示。但如果你编译了 HDRI 版 IM，也可以使用浮点“实部”和“虚部”成分来处理图像。这通过选项「[+fft](https://usage.imagemagick.org/fourier/option_link.cgi#fft)」和「[+ift](https://usage.imagemagick.org/fourier/option_link.cgi#ift)」的“加号”版本完成。例如，这里我使用 HDRI 版 IM 对图像执行一次“往返”FFT，但这次生成实部/虚部图像。 |

  # HDRI version of IM used
  time magick lena.png   +fft +ift   lena_roundtrip_ri.png


  echo -n "RMSE = "

  magick compare -metric RMSE lena.png lena_roundtrip_ri.png null:

  echo -n "PAE = "

  magick compare -metric PAE  lena.png lena_roundtrip_ri.png null:

使用加号形式生成实部/虚部 FFT 图像时，必须使用 HDRI 版本。如果不用，大约一半的值会变成零，结果图像会显得“脏”。例如…… |

  # non-HDRI Q16 version of IM used  -- THIS IS BAD
  magick lena.png   +fft +ift   lena_roundtrip_ri_bad.png

还要记住，不管你生成哪种形式的 FFT 图像，都会影响你想对 FFT 图像应用的所有图像处理操作。它们是非常不同的图像，因此必须以非常不同的方式、用不同的数学操作来处理。和前面一样，要恢复最终图像，必须同时拥有实部和虚部成分图像。例如，如果我们用一幅“黑色”图像替换其中一个成分，就会发生下面的情况。

  # HDRI version of IM used
  magick lena.png +fft -delete 1 \
          -size 128x128 xc:black +ift lena_real_only.png
  magick lena.png +fft -delete 0 \
          -size 128x128 xc:black +ift lena_imaginary_only.png

仅实部 |
仅虚部
---|---

从这里可以看出，实部/虚部两幅 FFT 图像都相当均衡地包含了关于原始图像的重要信息。两者最大的差别是，特殊的 DC 或“平均颜色”没有虚部成分，因此只存在于振幅图像中。你在两幅图像中看到的对角镜像（实际上是 180 度旋转）效果，是由于丢失了另一个成分中包含的“符号”信息造成的。没有另一个成分时，可以认为波相差 180 度，于是产生这种奇怪外观。这种信息损失在两类图像之间是对等的。

傅里叶变换的性质

常量图像的 FFT

我们来演示这些性质中的一些。先简单取一幅常量颜色图像，并得到它的振幅。

  magick -size 128x128 xc:gold constant.png
  magick constant.png -fft +delete constant_magnitude.png

注意，在这个例子中，振幅图像真的几乎是纯黑，只有图像正中心、像素位置 width/2、height/2 处有一个彩色像素。这个像素就是图像的零频率或 DC（“Direct Current”）值，也是唯一不表示正弦波的像素。换句话说，这个值就是 FFT 的常量成分！为了更清楚地看到这个单独像素，我们也把图像的这一区域放大…… |

  magick constant_magnitude.png -gravity center -extent 5x5 \
           -scale 2000% constant_dc_zoom.gif

注意，DC 点的颜色与原始图像相同。实际上，最好记住你看到的是三个值。也就是说，生成的图像实际上是三个独立的快速傅里叶变换：红、绿、蓝三个图像通道各有一个 FFT。FFT 本身并不真正理解颜色，只理解颜色值或“灰度级”。事实上，FFT 变换几乎可以应用于任何色彩空间，因为实际上……它并不在乎！对傅里叶变换来说，图像只是一个数值数组，仅此而已。 | 虽然 DC 值的“相位”并不重要，但它应始终是“零”角度（相位颜色值为 50% 灰）。如果没有设置为 50% 灰，DC 值会具有一个“不真实”的成分，并且其值会受到给定角度的调制。
---|---

DC 颜色的影响

在更典型的非常量图像中，DC 值是图像的平均颜色。它通常就是当你把图像完全模糊、平均，或缩小到单个像素或单种颜色时应得到的颜色。例如，我们从 "Lena" 图像的 FFT 中提取 DC 像素。

  magick lena.png -fft +delete lena_magnitude.png
  magick lena_magnitude.png -gravity center -extent 1x1 \
          -scale 60x60   lena_dc_zoom.gif

可以看到，这幅图像的平均颜色是一种“暗粉色”。另一种理解这个特殊像素的方式是：它表示中心“偏置”水平，其他所有正弦波都围绕它来修改图像颜色。 例如，我们把这个“暗粉色”DC 像素替换成其他颜色，比如更偏橙色的 'tomato'…… |

  magick lena.png -fft \
          \( -clone 0  -draw "fill tomato  color 64,64 point" \) \
          -swap 0 +delete -ift lena_dc_replace.png

真正发生的是，当你改变 FFT 图像中的 DC 值时，会以同样方式改变整幅图像。实际上，DC 值的任何变化（差值）都会被加到（或从中减去）结果图像的每一个像素上。这就像我们真的给原始图像中的每个像素都加上某个常量。因此，重建图像中的最终像素颜色也可能被最大（白）或最小（黑）限制裁剪。因此，这并不是推荐的图像调色方法。它比修改整幅图像中的每个像素更简单，但 FFT 往返会让它总体上成为慢得多的调色技巧。

正弦波图像的频谱

接下来，我们看一幅在图像宽度上有 4 个周期的单一正弦（或余弦）波图像的频谱。

  magick -size 128x129 gradient: -chop 0x1 -rotate 90 -evaluate sine 4 \
          sine4.png
  magick sine4.png -fft +delete \
          -auto-level -evaluate log 100  sine4_spectrum.png

| _上面创建渐变图像的特殊方式，是为了确保得到的正弦波图像能在图像范围内完美平铺。

普通的「[gradient:](canvas.html#gradient)」图像不能完美平铺，因此由它生成的正弦波也不能完美平铺。对这种不完美平铺的图像做 FFT 变换，会在傅里叶变换频谱中得到一组不需要的谐波，而不是单个“点”。

关于这个问题的更多细节，请参见生成完美渐变。

_
---|---
在上面的频谱图像（增强后的振幅图像）中，可以看到它有 3 个点。中心点和前面一样，是平均 DC 值。另外两个点表示傅里叶运算符在图像中找到的完美正弦波。由于图像宽度方向上的频率恰好是 4 个周期，因此两个频率像素距离中心 DC 值正好 4 个像素。为什么是两个像素？这是因为单一正弦波可以用两种完全不同的方式描述（其中一种带有负方向和相位）。这两种描述在数学上都正确，傅里叶变换不会区分它们。如果我们用 16 个周期的正弦波重复这个过程，会再次看到 3 个点，但这些点相距更远。在这个例子中，侧边的点位于中心点左右各 16 个像素处。

  magick -size 128x129 gradient: -chop 0x1 -rotate 90 -evaluate sine 16 \
          -write sine16.png -fft -delete 1 \
          -auto-level -evaluate log 100 sine16_spectrum.png

由此可以看到，完美正弦波会简单地表示为位于适当位置的两个点。这个位置离中心 DC 值的距离决定了正弦波的频率。波长越短，频率越高，因此这些点离 DC 值越远。事实上，用图像尺寸除以频率（点到中心的距离），就可以得到波长（峰值之间的距离）。在上面的例子中：128 像素除以 16 个周期，得到每个“条带”之间 8 像素的波长。这是 FFT 变换最重要的区分特征之一。原始图像上的小特征图案需要短波长，因此需要高频率。这会在频率域中产生大尺度效果。同样，大特征使用较低频率，因此生成小尺度图案，尤其是在接近中心的位置。在傅里叶变换中……

小的变大，大的变小。

这是处理傅里叶变换时最需要记住的方面之一，因为它是从图像中去除噪声（小特征），同时保留图像整体较大结构的关键。 我们通过绘制这三个“频率”的原始振幅（不是对数频谱），来更仔细地看看它们。 |

  magick sine16.png -fft -delete 1  miff:- |\
     im_profile - sine16_magnitude_pf.png

注意，DC 值（图像平均值或偏置）为 1/2，这是预期结果（图像平均值是完美的 50% 灰），但傅里叶变换找到的两个 16 周期正弦波的实际振幅只有最大值的 1/4。原始正弦波的振幅其实是 1/2，但傅里叶变换把这个振幅分成了两份，并把结果分配给两个绘制出的频率波，因此两个成分各自只有 1/4 的振幅。这是傅里叶变换中的正常现象。 FFT 图像中正频率与负频率的这种二重性，解释了为什么所有 FFT 图像频谱（例如左侧再次显示的 Lena 频谱）总是关于中心对称。图像一侧的每个点，都会在图像中心的旋转镜像位置得到一个类似的“点”。同样的事情也会发生在 FFT 图像对的“相位”成分中，只是值还带有 180 度偏移（负相位）。这意味着每幅图像的一半实际上是另一半的副本，但要重新创建原始图像，你仍然需要两幅图像。换句话说，这两幅图像仍包含完全相同数量的信息，一半在一幅图像中，另一半在另一幅图像中。它们合在一起才构成整体。 | _生成期间，FFT 算法只生成图像的左半部分。另一半通过旋转和复制已生成的数据得到。

把频率域图像转换回空间域图像时，算法也只查看图像左半部分。右半部分会被完全忽略，因为它只是副本。

因此，在后面的示例中，当你对 FFT 振幅图像进行“陷波滤波”时，实际上只需要滤掉振幅图像的左半边。你可以通过忽略右半边来省些工作。不过为了清楚起见，我会对两半都做“陷波”。

_
---|---

直接生成 FFT 图像

现在可以利用上面的信息，实际生成一幅正弦波图像。你需要做的只是创建一对黑色图像和 50% 灰色图像，并用合适的振幅和相位添加“点”。例如……

  magick -size 128x128  xc:black \
          -draw 'fill gray(50%)  color 64,64 point' \
          -draw 'fill gray(50%)  color 50,68 point' \
          -draw 'fill gray(25%)  color 78,60 point' \
          generated_magnitude.png
  magick generated_magnitude.png \
          -auto-level -evaluate log 3  generated_spectrum.png
  magick -size 128x128  xc:gray50  generated_phase.png
  magick generated_magnitude.png generated_phase.png \
          -ift  generated_wave.png

_
---|---
FUTURE: Perlin Noise Generator using FFT

垂直线的频谱

显示细线和粗线的 FFT 频谱 演示小特征如何在图像 FFT 中变“大”，大特征如何变“小”。再把它与可看作具有单一谐波的“线”的正弦波联系起来。 旋转这条线

矩形图案图像的频谱

接下来，我们看一下黑色背景中宽 8、高 16 的白色矩形的频谱。

  magick -size 8x16 xc:white -gravity center \
          -gravity center -background black -extent 128x128 rectangle.png
  magick rectangle.png -fft +delete \
          -auto-level -evaluate log 100 rect_spectrum.png

可以看到，结果图像具有非常特定的图案，并包含许多谐波频率。你还可以看到矩形似乎旋转了 90 度。其实并不是旋转了；你看到的是前面提到的同一条规则：大特征变小，小特征变大。因此，矩形较小的尺寸变大了，较大的尺寸变小了。现在，把矩形旋转 45 度。我们会发现频谱也按同一方向旋转了 45 度。

  magick rectangle.png -rotate 45 -gravity center -extent 128x128 \
          -write rect_rot45.png -fft -delete 1 \
          -auto-level -evaluate log 100 rect_rot45_spectrum.png

可以看到频率域中也有同样的旋转。也就是说，某个旋转对象的效果，在其傅里叶变换中也会被旋转。不过，如果现在移动这个矩形……

  magick rectangle.png -rotate 45  -geometry +30+20 -extent 128x128 \
          -write rect_rot45off.png -fft -delete 1 \
          -auto-level -evaluate log 100 rect_rot45off_spectrum.png

频率图案并没有移动。这是因为所有位置信息都包含在相位图像中。频率图案（振幅或其频谱）不会因为移动而改变。这种位置分离，是傅里叶变换的关键特征之一，也正是它如此重要的原因。它允许你在较大的图像中搜索特定图像图案，而不必关心产生该傅里叶频谱图案的对象位于何处。

平面圆形图案图像的频谱

接下来，我们看一幅带白色平面圆形图案的图像的频谱，一种情况直径为 12（半径 6），另一种情况直径为 24（半径 12）。

  magick -size 128x128 xc:black -fill white  \
          -draw "circle 64,64 64,70" -write circle6.png -fft -delete 1 \
          -auto-level -evaluate log 100 circle6_spectrum.png

  magick -size 128x128 xc:black -fill white  \
          -draw "circle 64,64 64,76" -write circle12.png -fft -delete 1 \
          -auto-level -evaluate log 100 circle12_spectrum.png

注意，第一幅图像非常接近我们在上面 jinc 示例中生成的结果。不过它有些破碎。这些伪影来自圆的尺寸太小。由于它以数字方式表示，周长并不是完美圆形。我们再次看到，小细节会在变换后的频率空间中变大。较大圆的变换效果更好，因为它的周长更接近真实圆。因此可以得出结论，平面圆形形状的变换确实是 jinc 函数，并且包含较小直径圆的图像会产生分布更分散、更宽的变换特征。根据傅里叶变换的数学性质，频谱中从中心到第一个暗环中部的距离是 1.22N/d。当圆的直径 d=12 时，距离为 1.22128/12=13。同样，当圆的直径 d=24 时，距离为 1.22*128/24=6.5。

高斯图案图像的频谱

接下来，我们看两幅图像的频谱，它们各自包含一个白色高斯圆形图案，sigma 分别为 8 和 16。

  magick -size 128x128 xc:black -fill white \
          -draw "point 64,64" -gaussian-blur 0x8 -auto-level \
          -write gaus8.png -fft -delete 1 \
          -auto-level -evaluate log 1000 gaus8_spectrum.png

  im_profile -s gaus8.png gaus8_pf.gif
  im_profile -s gaus8_spectrum.png gaus8_spectrum_pf.gif

  magick -size 128x128 xc:black -fill white \
          -draw "point 64,64" -gaussian-blur 0x16 -auto-level \
          -write gaus16.png -fft -delete 1 \
          -auto-level -evaluate log 1000 gaus16_spectrum.png

  im_profile -s gaus16.png gaus16_pf.gif
  im_profile -s gaus16_spectrum.png gaus16_spectrum_pf.gif

除了图案矩形数组产生的噪声之外，结果是：高斯图案会产生几乎相同的高斯频率图案。更重要的是，这个图案外观相当干净。当然，两者大小不同，再次遵循同一条规则：大的变小，小的变大。从数学性质来看，频谱中的 sigma 会正好是 N/(2*sigma)，其中 sigma 来自原始图像。因此，对于尺寸 N=128、sigma=8 的图像，频谱中的 sigma 为 128/16=8。同样，如果图像的 sigma 为 16，那么频谱中的 sigma 为 128/32=4。这就是“大变小，反之亦然”规则的数学关系，了解它会很有用。

网格图案图像的频谱

接下来，我们变换一幅只包含间隔为 16x8 像素的网格线图像。

  magick -size 16x8 xc:white -fill black \
          -draw "line 0,0 15,0" -draw "line 0,0 0,7" \
          -write mpr:tile +delete \
          -size 128x128 tile:mpr:tile \
          -write grid16x8.png -fft -delete 1 \
          -auto-level -evaluate log 100000 grid16x8_spectrum.png

得到的频谱只是一个点阵：网格线间距越近，产生的点间距越远，反之亦然。根据上面的性质，由于网格线间隔为 16x8 像素，这些点的间距应为 N/a=128/16=8 和 M/b=128/8=16，这正是图中测得的结果。这个图案尤其重要，因为它会告诉你傅里叶变换与图像中规则平铺图案之间的关系。这类平铺图案会在傅里叶变换中产生很强的非中心网格图案。这里的关键点是，形状信息位于中心，而平铺信息位于傅里叶变换中心之外的网格状数组中。

更多频谱信息

如果你想进一步了解频谱图像及其性质，下面有一些链接。

• Wikipedia: Fourier Transform
• Fred Weinhaus, Properties of a Fourier Transform
• Wolfram MathWorld: Fourier Transform

实际应用

好了，既然已经讲完基础，使用傅里叶变换有哪些实际应用？可以完成的事情包括：1）增加或降低图像对比度，2）模糊，3）锐化，4）边缘检测，5）去除噪声。

改变图像对比度 - 系数开方

可以通过执行正向傅里叶变换，将振幅图像提升到某个幂次，然后把它与相位一起用于逆傅里叶变换，从而调整图像对比度。要增加对比度，使用略小于 1 的指数；要降低对比度，使用略大于 1 的指数。那么，先用指数 0.9 增加 Lena 图像的对比度，再用指数 1.1 降低对比度。

  magick lena.png -fft \
          \( -clone 0 -evaluate pow 0.9 \) -delete 0 \
          +swap -ift lena_plus_contrast.png

  magick lena.png -fft \
          \( -clone 0 -evaluate pow 1.1 \) -delete 0 \
          +swap -ift lena_minus_contrast.png

不过，对原始图像这样做，效果也与直接对原始图像这样做相同。也就是说，对振幅进行全局修改，与对原始图像做全局修改具有相同效果。

模糊图像 - 低通滤波

傅里叶变换最重要的性质之一，是空间域中的卷积等价于频率域中的简单乘法。在空间域中，可以用小的、方形的简单卷积滤波器（核），通过 -convolve 选项来模糊图像。这称为低通滤波器。最简单的滤波器只是一个等权重的方形数组。也就是说，所有值都是 1，并在应用卷积前通过除以总和来归一化。这等价于局部或邻域平均。另一种低通滤波器是高斯加权的圆形滤波器，由 -gaussian-blur 或 -blur 提供。在频率域中，一种低通模糊滤波器只是一个被黑色包围的恒定强度白色圆。这类似于空间域中的圆形平均卷积滤波器。不过，由于空间域中的卷积等价于频率域中的乘法，我们只需要执行正向傅里叶变换，然后把滤波器与振幅图像相乘，最后执行逆傅里叶变换。注意，小尺寸卷积滤波器会对应频率域中的大圆。乘法通过 -composite 并使用 -compose multiply 设置来完成。那么，我们用两种尺寸的圆形滤波器试一试，一个直径为 40（半径 20），另一个直径为 28（半径 14）。

  magick -size 128x128 xc:black -fill white \
          -draw "circle 64,64 44,64" circle_r20.png
  magick lena.png -fft \
       \( -clone 0 circle_r20.png -compose multiply -composite \) \
       \( +clone -evaluate log 10000 -write lena_blur_r20_spec.png +delete \) \
       -swap 0 +delete -ift lena_blur_r20.png

  magick -size 128x128 xc:black -fill white \
          -draw "circle 64,64 50,64" circle_r14.png
  magick lena.png -fft \
       \( -clone 0 circle_r14.png -compose multiply -composite \) \
       \( +clone -evaluate log 10000 -write lena_blur_r14_spec.png +delete \) \
       -swap 0 +delete -ift lena_blur_r14.png

可以看到，使用较小直径滤波器的图像产生了更强的模糊。还可以注意到结果图像边缘附近有“振铃”或“波纹”效果。这是因为圆的傅里叶变换，如前面所见，是一个 jinc 函数，它会随着从中心向外扩展而产生逐渐减弱的振荡。不过这里 jinc 函数和这些振荡位于空间域，而不是前面演示过的频率域。那么能怎么办？最简单的做法，是使用各种窗函数让圆的边缘渐变。另一种做法，是使用高斯形状这种按定义已经渐变的滤波器。那么我们采用后一种，用两个高斯模糊圆来去除大部分严重的“振铃”效果。

  magick circle_r20.png -blur 0x4 -auto-level gaussian_r20.png
  magick lena.png -fft \
       \( -clone 0 gaussian_r20.png -compose multiply -composite \) \
       \( +clone -evaluate log 10000 -write lena_gblur_r20_spec.png +delete \) \
       -swap 0 +delete -ift lena_gblur_r20.png

  magick circle_r14.png -blur 0x4 -auto-level gaussian_r14.png
  magick lena.png -fft \
       \( -clone 0 gaussian_r14.png -compose multiply -composite \) \
       \( +clone -evaluate log 10000 -write lena_gblur_r14_spec.png +delete \) \
       -swap 0 +delete -ift lena_gblur_r14.png

这当然好得多。理想的低通滤波器并不是去模糊圆，而是实际使用按 sigma 而非 radius 定义的合适高斯曲线。当然，在这个例子中，我们最后还是为了做模糊而做了一次模糊！不过，乘到 FFT 振幅图像上的模糊图案是固定的，实际上可以从预生成缓存中取出。此外，用来相乘的图像不必与原图同样大小，可以使用较小图像。因此，在处理大图像以及大量图像时，上面的做法可以快很多。更重要的是，对于大而强的模糊，频率域图像较小，并且只做一次乘法，而不用对原始图像中的每个像素平均许多像素。对于小尺寸模糊，直接卷积模糊可能更合适。

检测图像边缘 - 高通滤波

在空间域中，从图像中提取边缘的高通滤波器，通常实现为正负权重相加为零的卷积。在频率域中，事情简单得多。这里高通滤波器只是低通滤波器的取反版本。也就是说，低通滤波器明亮的位置，高通滤波器为暗，反之亦然。因此在 ImageMagick 中，我们只需要对低通滤波器图像执行 -negate。那么，先用圆形图像对 Lena 图像应用高通滤波器，再用纯高斯曲线做一次。

  magick circle_r14.png -negate circle_r14i.png
  magick lena.png -fft \
       \( -clone 0 circle_r14i.png -compose multiply -composite \) \
       \( +clone -evaluate log 10000 -write lena_edge_r14_spec.png +delete \) \
       -delete 0 +swap -ift -normalize lena_edge_r14.png

  magick -size 128x128 xc: -draw "point 64,64" -blur 0x14 \
          -auto-level   gaussian_s14i.png
  magick lena.png -fft \
       \( -clone 0 gaussian_s14i.png -compose multiply -composite \) \
       \( +clone -evaluate log 10000 -write lena_edge_s14_spec.png +delete \) \
       -delete 0 +swap -ift -normalize lena_edge_s14.png

仔细观察这两个结果，可以看到简单圆形不如高斯好，因为它有“振铃”伪影，而且不够锐利。

锐化图像 - 高提升滤波

锐化图像最简单的方式，是对它进行高通滤波（不做归一化拉伸），然后与原始图像混合。

  magick lena.png -fft \
       \( -size 128x128 xc: -draw "point 64,64" -blur 0x14 -auto-level \
          -clone 0 -compose multiply -composite \) \
       -delete 0 +swap -ift \
       lena.png -compose blend -set option:compose:args 100x100 -composite \
       lena_sharp14.png

这里在频率域中进行高通滤波，然后把结果变换回空间域，并与原始图像混合，以增强图像边缘。

去除噪声 - 陷波滤波

许多有噪声的图像包含某种图案化噪声。这种噪声在频率域中很容易去除，因为图案会表现为少量点或线组成的图案。回想一下，简单正弦波就是重复图案，在频谱中只显示为 3 个点。为了去除这种噪声，遗憾的是必须手动遮罩（或陷波）掉振幅图像中的这些点或线。做法是先变换到频率域，创建频谱的灰度版本，遮罩点或线，对它做阈值处理，把二值蒙版图像与振幅图像相乘，再变换回空间域。我们在小丑图像上试试，它包含一种对角条纹的类似抖动图案。首先变换小丑图像，创建它的振幅和相位图像。

  magick clown_orig.jpg -fft \
          \( +clone  -write clown_phase.png +delete \) +delete \
          -write clown_magnitude.png  -colorspace gray \
          -auto-level -evaluate log 100000  clown_spectrum.png

原图 | |
频谱 |
相位
---|---|---|---

可以看到，频谱中有四个明亮的星状点，每个象限各一个。这些异常点表示图像中我们想去掉的图案。图像中部的亮点和线不必担心，因为它们表示 DC（图像平均颜色）以及图像边缘的影响，不应被修改。注意，在生成频谱图像时，我强制结果图像为纯灰度图像。这样我现在就可以把图像载入编辑器，并使用任意非灰色颜色（例如红色）遮住那 4 个星状图案所在区域。编辑完成后，可以通过与未编辑版本提取差异图像来提取我涂色的区域。像这样……

  magick clown_spectrum_edited.png clown_spectrum.png \
          -compose difference -composite \
          -threshold 0 -negate clown_spectrum_mask.png

现在只需把蒙版与振幅相乘，并把结果与原始相位图像一起变换回空间域。旁边显示原图以便比较。

  magick clown_magnitude.png clown_spectrum_mask.png \
          -compose multiply -composite \
          clown_phase.png -ift clown_filtered.png

结果非常好。但我们还能做得更好。正如你在前面的示例中看到的，简单“圆形”对 FFT 图像并不特别友好，所以我们稍微模糊一下蒙版…… |

  magick clown_spectrum_mask.png \
          -blur 0x5 -level 50x100%  clown_mask_blurred.png

然后过滤小丑图像，这次在内存中重新生成 FFT 图像。 |

  magick clown_orig.jpg -fft \
          \( -clone 0 clown_mask_blurred.png -compose multiply -composite \) \
          -swap 0 +delete -ift clown_filtered_2.png

结果简直惊人！如果进一步调整蒙版，使它更贴合那些“星形”，也许还能继续改进。 我们甚至可以取原图与结果之间的差异，生成一幅显示噪声被移除区域的图像。 |

  magick clown_orig.jpg clown_filtered_2.png -compose difference \
          -composite -normalize clown_noise.png

再试另一个例子。这次使用 RoboRealm 网站上的 "Twigs" 图像，它包含不规则的水平和垂直条纹图案。和之前一样，我们提取灰度频谱图像。

  magick twigs.jpg -fft +delete -colorspace gray \
          -auto-level -evaluate log 100000 twigs_spectrum.png

在这个例子中，由于图像中的噪声呈水平和垂直方向，它会显示为沿中心线分布、但不在图像真正中心处的粗水平带和垂直带。我们再次用图像编辑器遮住这些部分，这次使用“蓝色”（实际上用什么颜色都无所谓）……

  magick twigs_spectrum_edited.png twigs_spectrum.png \
          -compose difference -composite \
          -threshold 0 -negate twigs_spectrum_mask.png

现在再次把蒙版与 FFT 振幅图像相乘，并重建图像。

  magick twigs.jpg -fft \
          \( -clone 0 twigs_spectrum_mask.png -compose multiply -composite \) \
          -swap 0 +delete  -ift twigs_filtered.png

还可以取原图与结果之间的差异，生成一幅显示噪声被移除区域的图像。

  magick twigs.jpg twigs_filtered.png -compose difference -composite \
          -normalize twigs_noise.png

给蒙版稍微加一点模糊，也许还能进一步改善结果。作为练习，试着从图像中移除那根细线。提示：记住真实图像中一条线的效果，在 FFT 中会旋转 90 度。如果弄错了，你很可能会把树枝移除掉。

高级应用

使用傅里叶变换的其他更高级应用包括：1）对运动模糊和失焦图像进行反卷积（去模糊），以及 2）使用归一化互相关，在大图中找出小图最匹配的位置。FFT 乘法和除法（反卷积）的示例已移到一个子目录，因为它还在等待更正式定义的图像处理运算符。