Digital Image Processing, 4th ed. Gonzalez & Woods 读书笔记

  三月中旬我的文章刚刚被 Int. J. Plast. 接收，PMI航空泡沫结构性能关系相关的工作也暂时告一段落了。近三个月得空处理了积压已久的历史遗留问题，卸下包袱方便全力攻克新的难点。去年末本团队依靠 Dynamic CT 技术采集了聚氨酯泡沫发泡的动态三维结构演化过程，共采集 15 套动态 CT 数据，每套含 80 组尺寸为 [1024,1024,1024] 的三维图像，共 2400 组。面对如此海量的四维数据集，如若仍然坚持采用传统数据分析方法进行数据分割、量化分析，需要四位经过充分训练的专业人员持续工作九十周整。重复的机械劳动是对高质量人才的浪费，罔论技术能力与科学观点的创新价值不容拖延，传统方法无法接受。因而发展基于机器学习的全自动数据分割方法刻不容缓。

  了解 Gonzalez 所著的 Digital Image Processing, 4th ed. 是在 Prof. Sapiro 的公开课程 Image and Video Processing 中。然而公开课的内容往往浅尝辄止，阅读本书也是为了打开思路，学习计算机视觉领域内本人尚未了解的新颖技术，对知识体系查缺补漏，进而辅助数据分析方法的设计。

Basic Knowledge

  Sampling, Quantization 分别指对图像空间尺度、位阶进行重采样。空间分辨能力（spatial resolution）为图像中可以分辨的最小细节物理尺寸，由于离散采样的数学特性，空间分辨能力必定劣于两倍像素尺寸，此既奈奎斯特采样定理。2016 年东海大学 Prof. Ryuta 提出一种依靠频域空间特征估算模糊图像的卷积核点扩散函数（Point Spread Function, PSF）半高宽的方法[1]，进而直接换算得到图像空间分辨能力。该方法为完全依托图像自身特征的后验方法，无需任何额外实验，非常适合用于辅助CT实验分析、闪烁体对焦等工作，依据已发表论文复现的程序见文末。

随机变量的期望值和矩（Expected Values and Moments）

连续随机变量的均值及方差定义为：

$$\overline{z}=E[z]=\int_{-\infty}^{\infty}zp(z)dz \tag{2-106}\label{2-106}$$

及：

$$\sigma^2=E[(z-\overline{z})^2]=\int_{-\infty}^{\infty}(z-\overline{z})^2p(z)dz = E[z^2]-\overline{z}^2 \tag{2-108}\label{2-108}$$

连续随机变量的n阶中心矩为：

$$\mu_n=E[(z-\overline{z})^n]=\int_{-\infty}^{\infty}(z-\overline{z})^np(z)dz \tag{2-109}\label{2-109}$$

其中$z$为随机变量的取值，亦为积分下的哑变量（a dummy variable of integration）。

对于已观测的离散随机变量，则需通过下式估算样本方差：

$$\widehat{\sigma}^2=\frac{1}{K-1}\sum_{k=1}^{K}(z_k-\overline{z})^2 \tag{2-114}\label{2-114}$$

  这与$Eq.\ref{2-108}$略有不同，是对样本总体方差的无偏估计（unbiased estimate）。无偏第三矩（unbiased third moments）和第四矩的表达式也很重要，因为它们用于定义偏度和峰度（skewness and kurtosis）。种群偏度（poulation skewness）定义为:

$$s = \frac{\mu_3}{\sigma^3} \tag{2-116}\label{2-116}$$

  第三矩是对样本对称性的度量，将其除以标准差立方使其标准化。当使用偏度作为对称性的度量时，我们主要关注的是它是零（严格对称），为正（向左倾斜），还是负值（向右倾斜）。如果数据集代表抽样（sample）而非总体（population），我们则可使用该表达式测算样本偏度。偏度的无偏估计为：

$$\widehat{s}=\frac{\sqrt{K(K-1)}}{K-2}s \tag{2-117}\label{2-117}$$

总体峰度定义为：

$$\kappa=\frac{\mu_4}{\sigma^4} \tag{2-118}\label{2-118}$$

  具有单位方差高斯分布的数据集四阶矩为3。四阶矩除以标准差的四次方以使其归一化。对于任何高斯总体，无论方差如何取值，峰度均为3。从峰度中减三得到总体过剩峰度（poplation excess kurtosis）：

$$\mathcal{K}=\kappa-3=\frac{\mu_4}{\sigma^4}-3 \tag{2-119}\label{2-119}$$

  因此，总体过剩峰度$\mathcal{K}=0$表明测算样本可能来自于高斯总体。如果该集合代表一组样本而非总体，则需使用该表达式获得无偏样本过剩峰度：

$$\widehat{\mathcal{K}}=\frac{K-1}{(K-2)(K-3)}[(K+1)\mathcal{K}+6] \tag{2-120}\label{2-120}$$

贝叶斯规则（Bayes’ Rule）

假设事件$A$已经发生，分支事件$B_k$的概率可以被表示为：

$$P(B_k|A)=\frac{P(A|B_k)P(B_k)}{P(A)}=\frac{P(A|B_k)P(B_k)}{\sum_{i=1}^{n}P(A|B_i)P(B_i)} \tag{2-80}\label{2-80}$$

  最后一个方程的分母遵循总概率定律。这是概率论中一个重要的成果，它告诉我们如何根据实时状况判断当前的条件概率，既如何通过$P(A|B)$判断$P(B|A)$。通常，分区只包含一个集合$B$及其补集$B^c$。则贝叶斯规则可以简写为：

$$P(B|A)=\frac{P(A|B)P(B)}{P(A|B)P(B)+P(A|B^c)P(B^c)} \tag{2-81}\label{2-81}$$

狄拉克函数（Dirac Delta Function）

  这不是一个通常意义上的函数，而是广义函数（generalized function）。然而，文献中往往也称其为脉冲函数（impulse function），狄拉克函数（delta function, Dirac delta function）。其满足如下形式：

$$\delta(t)=\left\{ \begin{aligned} +\infty\ \ & {\rm if\ }t=0\\ 0\ \ \ \ & {\rm if\ }t\neq0 \end{aligned} \right. \tag{4-10}\label{4-10}$$

且在整个实数域的积分为1：

$$\int_{-\infty}^{+\infty}\delta(t)dt=1 \tag{4-11}\label{4-11}$$

由于这样的独特定义，狄拉克函数具有积分下的筛选特性（sifting property）。

$$\int_{-\infty}^{+\infty}f(t)\delta(t)dt=f(0) \tag{4-12}\label{4-12}$$

筛选特性的更一般的表述涉及到位于任意点的脉冲，在这种情况下表示为：

$$\int_{-\infty}^{+\infty}f(t)\delta(t-t_0)dt=f(t_0) \tag{4-13}\label{4-13}$$

由狄拉克函数衍生而来的一个有趣的概念称作脉冲序列$s_{\Delta T}(t)$（impulse train），定义为互相间隔为$\Delta T$的无穷多个狄拉克函数之和：

$$s_{\Delta T}(t)=\sum_{k=-\infty}^{+\infty}\delta(t-k\Delta T) \tag{4-14}\label{4-14}$$

在离散系统中，单位离散脉冲$\delta(x)$（unit discrete impulse）的作用与处理连续变量时的脉冲函数$\delta(t)$是相同的。它被定义为：

$$\delta(x)=\left\{ \begin{aligned} 1\ \ \ \ & {\rm if\ }x=0\\ 0\ \ \ \ & {\rm if\ }x\neq0 \end{aligned} \right. \tag{4-15}\label{4-15}$$

同样在整个定义域内的累加为1：

$$\sum_{x=-\infty}^{+\infty}\delta(x)=1 \tag{4-16}\label{4-16}$$

离散系统下的筛选特性则表现为如下一般形式：

$$\sum_{x=-\infty}^{+\infty}f(x)\delta(x-x_0)=f(x_0) \tag{4-18}\label{4-18}$$

傅里叶变换（Fourier Transform）与傅里叶级数（Fourier Series）

考虑一个自变量为$t$的连续函数$f(t)$，定义域为$\mathbb{R}^1$。若其为以$T$为间隔的周期函数，则该函数可以表示为正弦、余弦函数乘以特定系数的和。这被称做傅里叶级数，有如下形式：

$$f(t)=\sum_{n=-\infty}^{+\infty}C_ne^{j\frac{2\pi n}{T}t} \tag{4-8}\label{4-8}$$

其中系数$C_n$为：

$$C_n=\frac{1}{T}\int_{-T/2}^{T/2}f(t)e^{-\frac{j2\pi nt}{T}}dt,\ for\ \forall n\in\mathbb{Z}^1 \tag{4-9}\label{4-9}$$

考虑一个自变量为$t$的一般连续函数$f(t)$，定义域为$\mathbb{R}^1$，其傅里叶变换记作$\mathcal{F}\{f(t)\}$或$F(\mu)$，有如下形式：

$$F(\mu)=\int_{-\infty}^{+\infty}f(t)e^{-j2\pi\mu t}dt \tag{4-19}\label{4-19}$$

其中$\mu$亦为连续变量。反过来，在已知傅里叶变换$F(\mu)$的情况下可以通过傅里叶逆变换得到空间域的连续函数$f(t)$：

$$f(t)=\int_{-\infty}^{+\infty}F(\mu)e^{j2\pi\mu t}d\mu \tag{4-21}\label{4-21}$$

根据欧拉公式，$Eq.\ref{4-19}$可被改写为：

$$F(\mu)=\int_{-\infty}^{+\infty}f(t)[{\rm cos}(2\pi\mu t)-j{\rm sin}(2\pi\mu t)]dt \tag{4-22}\label{4-22}$$

注意，傅里叶变换$F(\mu)$是$f(t)$乘以各正弦项的展开式，正弦项的频率由$\mu$决定。由于积分后唯一剩余的变量为频率，我们也称傅里叶变换后的空间为频域。一般来说，傅里叶变换包含复数项，因而也常常使用傅里叶变换的幅值$|F(\mu)|$，这被称为傅里叶频谱（Fourier spectrum or frequency spectrum）。

对于发生平移的函数$h(t-\tau)$，$\tau$为常数，其傅里叶变换有如下形式：

$$\begin{aligned} \mathcal{F}\{h(t-\tau)\} & =\int_{-\infty}^{+\infty}h(t-\tau)e^{-j2\pi\mu t}dt\\ & \\ & =e^{-2\pi\mu\tau}\int_{-\infty}^{+\infty}h(t-\tau)e^{-j2\pi\mu (t-\tau)}d(t-\tau)\\ & \\ & =H(\mu)e^{-2\pi\mu\tau} \end{aligned} \tag{4-100}\label{4-100}$$

例如，狄拉克函数的傅里叶变换：

$$\mathcal{F}\{\delta(t-t_0)\}=F(\mu)=\int_{-\infty}^{+\infty}e^{-j2\pi\mu t}\delta(t-t_0)dt=e^{-j2\pi\mu t_0} \tag{4-23}\label{4-23}$$

对于脉冲序列函数$s_{\Delta T}(t)$，其为以$\Delta T$为间隔的周期函数，因而可以通过傅里叶级数展开：

$$s_{\Delta T}(t) = \sum_{n=-\infty}^{+\infty}c_ne^{j\frac{2\pi n}{\Delta T}t}$$

其中：

$$c_n=\frac{1}{\Delta T}\int_{-\Delta T/2}^{\Delta T/2}s_{\Delta T}(t)e^{-j\frac{2\pi n}{\Delta T}t}dt=\frac{1}{\Delta T}$$

因而，脉冲序列函数的傅里叶级数变为：

$$s_{\Delta T}(t)=\frac{1}{\Delta T}\sum_{n=-\infty}^{+\infty}e^{j\frac{2\pi n}{\Delta T}t}$$

考虑到求和过程为线性运算，求和后函数的傅里叶变换与求和公式中各个分量的傅里叶变换之和是一样的。而这些分量的傅里叶变换如下：

$$\mathcal{F}\Big\{e^{j\frac{2\pi n}{\Delta T}t}\Big\}=\delta\Big(\mu-\frac{\mu}{\Delta T}\Big)$$

则脉冲序列函数的傅里叶变换为：

$$S(\mu)=\mathcal{F}\Big\{\frac{1}{\Delta T}\sum_{n=-\infty}^{+\infty}e^{j\frac{2\pi n}{\Delta T}t}\Big\}=\frac{1}{\Delta T}\sum_{n=-\infty}^{+\infty}\delta\Big(\mu-\frac{n}{\Delta T}\Big)$$

这个结果告诉我们周期性脉冲序列函数的傅里叶变换也是一个脉冲序列，其间隔为$1/\Delta T$。

卷积（Convolution）

考虑两个定义域为$\mathbb{R}^1$的连续函数$f(t)$及$h(t)$，其对连续变量$t$的卷积被定义为：

$$(f\otimes h)(t)=\int_{-\infty}^{+\infty}f(\tau)h(t-\tau)d\tau \tag{4-24}\label{4-24}$$

而卷积的傅里叶变换满足如下规律：

$$\begin{aligned} \mathcal{F}\{(f\otimes h)(t)\} & =\int_{-\infty}^{+\infty}\Bigg[\int_{-\infty}^{+\infty}f(\tau)h(t-\tau)d\tau\Bigg]e^{-j2\pi\mu t}dt\\ & \\ & =\int_{-\infty}^{+\infty}f(\tau)\Bigg[\int_{-\infty}^{+\infty}h(t-\tau)e^{-j2\pi\mu t}dt\Bigg]d\tau\\ & \\ & =\int_{-\infty}^{+\infty}f(\tau)[H(\mu)e^{-j2\pi\mu\tau}]d\tau\\ & \\ & =H(\mu)\int_{-\infty}^{+\infty}f(\tau)e^{-j2\pi\mu\tau}d\tau\\ & \\ & =H(\mu)F(\mu) \end{aligned} \tag{4-25}\label{4-25}$$

  值得注意的是，卷积是可交换的，卷积表达式中函数的顺序无关紧要。$Eq.\ref{4-25}$表明两个连续函数在空间域的卷积的傅里叶变换等于两个函数傅里叶变换的乘积。换句话说，$f\otimes h$与$H\cdot F$为一组傅里叶变换对（Fourier transform pair），满足如下规律，注意双向箭头向右表示傅里叶变换，向左表示傅里叶逆变换。

$$(f\otimes h)(t)\Leftrightarrow(H\cdot F)(\mu) \tag{4-26}\label{4-26}$$ $$(f\cdot h)(t)\Leftrightarrow(H\otimes F)(\mu)$$

离散采样与离散傅里叶变换（Discrete Fourier Transform）

  考虑一个自变量为$t$的连续函数$f(t)$，定义域为$\mathbb{R}^1$。我们希望以均匀间隔$\Delta T$对其采样，计算离散采样$\tilde{f}(t)$的一种方法是将目标函数$f(t)$与一个间隔为$\Delta T$的脉冲函数序列$s_{\Delta T}(t)$相乘。

$$\tilde{f}(t) = f(t)s_{\Delta T}(t) = \sum_{n = -\infty}^{+\infty}f(t)\delta(t-n\Delta T) \tag{4-27}\label{4-27}$$

注意$\tilde{f}(t)$为原始连续函数与脉冲函数序列之积，而离散采样在该处真正的取值$f_k$等于局部的积分：

$$f_k=\int_{-\infty}^{+\infty}f(t)\delta(t-k\Delta T)dt=f(k\Delta T),\ for\ \forall k\in\mathbb{Z}^1 \tag{4-28}\label{4-28}$$

令$F(\mu)$表示连续函数$f(t)$的傅里叶变换，采样函数$\tilde{f}(t)$的傅里叶变换$\tilde{F}(\mu)$为：

$$\tilde{F}(\mu)=\mathcal{F}\{\tilde{f}(t)\}=\mathcal\{f(t)s_{\Delta T}(t)\}=(F\otimes S)(\mu) \tag{4-29}\label{4-29}$$

其中，离散脉冲函数序列的傅里叶变换$S(\mu)$为：

$$S(\mu)=\frac{1}{\Delta T}\sum_{n=-\infty}^{+\infty}\Delta(\mu-\frac{n}{\Delta T}) \tag{4-30}\label{4-30}$$

因而可以直接应用卷积的傅里叶变换特性得到：

$$\begin{aligned} \tilde{F}(\mu) & =(F\otimes S)(\mu)=\int_{-\infty}^{+\infty}F(\tau)(\mu-\tau)d\tau\\ & \\ & =\frac{1}{\Delta T}\int_{-\infty}^{+\infty}F(\tau)\sum_{n=-\infty}^{+\infty}\delta\Big(\mu-\tau-\frac{n}{\Delta T}\Big)d\tau\\ & \\ & = \frac{1}{\Delta T}\sum_{n=-\infty}^{+\infty}\int_{-\infty}^{+\infty}F(\tau)\delta\Big(\mu-\tau-\frac{n}{\Delta T}\Big)d\tau\\ & \\ & = \frac{1}{\Delta T}\sum_{n=-\infty}^{+\infty}F\Big(\mu-\frac{n}{\Delta T}\Big) \end{aligned} \tag{4-31}\label{4-31}$$

由此可见，采样函数的傅里叶变换$\tilde{F}(\mu)$为原始连续函数的傅里叶变换$F(\mu)$无限叠加的平移副本求和。在此我们应当注意到，即便$\tilde{f}(t)$为离散采样函数，它的傅里叶变换$\tilde{F}(\mu)$仍然是连续函数。

$Eq.\ref{4-31}$给出了由原始连续函数$f(t)$推导出采样函数傅里叶变换$\tilde{F}(\mu)$的方法，然而却没有给出直接由采样函数$\tilde{f}(t)$计算其自身傅里叶变换的方法，我们直接从$Eq.\ref{4-19}$中傅里叶变换的定义出发：

$$\tilde{F}(\mu)=\int_{-\infty}^{+\infty}\tilde{f}(t)e^{-j2p\mu t}dt \tag{4-39}\label{4-39}$$

代入$Eq.\ref{4-27}$，可以推导出：

$$\begin{aligned} \tilde{F}(\mu) & = \int_{-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty}f(t)\delta(t-n\Delta T)e^{-j2\pi\mu t}dt\\ & \\ & = \sum_{n=-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(t)\delta(t-n\Delta T)e^{-j2\pi\mu t}dt\\ & \\ & = \sum_{n=-\infty}^{+\infty}f_ne^{-j2\pi\mu n\Delta T} \end{aligned} \tag{4-40}\label{4-40}$$

由$Eq.\ref{4-31}$可知，即便$f_n$为离散数列，其傅里叶变换$\tilde{F}(\mu)$仍然是连续的，以$1/\Delta T$为间隔的无限周期函数。因此，我们只需要记录其单个周期即可完全保留该无限周期函数的信息，这是离散傅里叶变换（DFT）的基础。

假使我们获取$\tilde{F}(\mu)$一个周期内M个等间隔的采样，取值范围$\mu=0$到$\mu=1/\Delta T$，这是通过以下频率的采样来实现的：

$\mu=\frac{m}{M\Delta T}, m=0,1,2,…,M-1 \tag{4-41}\label{4-41}$

将上方$\mu$的取值代入$Eq.\ref{4-40}$，并将结果记为$F_m$：

$$F_m=\sum_{n=1}^{M-1}f_ne^{-j2\pi mn/M},\ \ m=0,1,2,...,M-1 \tag{4-42}\label{4-42}$$

这个表达式即为离散傅里叶变换。给定连续函数$f(t)$的M个采样数据集$\{f_m\}$，产生一组复数数据集$\{F_m\}$以表达输入数据集的离散傅里叶变换。相反，假设我们已知$\{F_m\}$，则可通过离散傅里叶反变换（IDFT）恢复样本集$\{f_m\}$。

$$f_n=\frac{1}{M}\sum_{m=0}^{M-1}F_me^{j2\pi mn/M},\ \ n=0,1,2,...,M-1 \tag{4-43}\label{4-43}$$

这两组恒等式表明，对于任意有限值样本集，傅里叶正变换和傅里叶反变换都存在。注意，离散傅里叶变换既不明确依赖于采样间隔，也不依赖于频率间隔，因此，DFT对适用于任意一致获取的有限离散数据集。

在前文的推导过程中，我们均使用$m,n$来表示离散自变量，以便与连续自变量区分开来。然而，为了直观起见，特别是在图像处理的二维空间中，我们通常用符号$x,y$表示图像坐标变量，用符号$u,v$表示频域自变量，并且它们通常为整数。而后，离散傅里叶正变换（forward DFT）及反变换（inverse DFT）写作：

$$F(u)=\sum_{x=0}^{M-1}f(x)e^{-j2\pi ux/M},\ \ u=0,1,2,...,M-1 \tag{4-44}\label{4-44}$$ $$f(x)=\frac{1}{M}\sum_{u=0}^{M-1}F(u)e^{j2\pi ux/M},\ \ x=0,1,2,...,M-1 \tag{4-45}\label{4-45}$$

为了简单起见，我们使用了函数符号而非下标。从现在开始，我们使用$Eq.\ref{4-44},ref{4-45}$表示一维离散傅里叶变换对。对于发生平移的函数，其傅里叶变换有如下形式：

$$f(x-x_0)\Leftrightarrow F(u)e^{-j2\pi ux_0/M}$$

这一公式在基于傅里叶变换方法计算相关性的数字图像相关技术中非常重要。不难看出，无论是傅里叶正变换或是逆变换的输出都是无限周期性函数，周期为$M$：

$$F(u) = F(u+kM) \tag{4-46}\label{4-46}$$ $$f(x) = f(x+kM) \tag{4-47}\label{4-47}$$

其中$k\in\mathbb{R}^1$。

在离散采样的情况下，$Eq.\ref{4-24}$中的卷积表达式转变为：

$$f(x)\otimes h(x) = \sum_{m=0}^{M-1}f(m)\cdot h(x-m),\ \ x=0,1,2,...,M-1 \tag{4-48}\label{4-48}$$

如果有限数据集${f(x)}$由连续函数$f(t)$的$M$个均匀间隔的样本组成，则组成该集合的长度为：

$$T = M\cdot \Delta T \tag{4-49}\label{4-49}$$

而相应的在频域中的均匀间隔为$\Delta u$，

$$\Delta u = \frac{1}{M\Delta T} = \frac{1}{T} \tag{4-50}\label{4-50}$$

离散傅里叶变换后M个分量组成的整个频域范围为：

$$R = M\Delta u = \frac{1}{\Delta T} \tag{4-51}\label{4-51}$$

采样定理（Sampling Theorem）

在上节中我们直观地介绍了抽样的思想。现在我们正式地考虑抽样，并推导了连续函数从其抽样中唯一恢复的条件。

考虑一个自变量为$t$的连续函数$f(t)$，定义域为$\mathbb{R}^1$，其傅里叶变换在有限区间$[-\mu_{max},\mu_{max}]$以外均为零，其称为带限函数（band-limited function）。严格采样函数（critically sampled function）指采样函数的傅里叶变换在平移叠加时恰好没有互相叠加的情况，既$\mu_{max}=1/2\Delta T$。更高的$\Delta T$取值将导致$\tilde{F}(\mu)$中相邻周期峰位互相叠加合并，而更低的$\Delta T$值将令相邻周期峰位更加清晰的分离开来。

若我们可以从无限平移叠加的$\tilde{F}(\mu)$中唯一分离出连续函数的傅里叶变换$F(\mu)$，即可精准的恢复原始信号$f(t)$，即需满足如下条件：

$$1/\Delta T>2\mu_{max} \tag{4-32}\label{4-32}$$

该方程表明，如果以超过函数最高频率的两倍的速率获取一个连续的带限函数，则可以从其抽样中完全恢复原函数。这个非常重要的结果被称为采样定理。带限函数要求$f(t)$的取值必须延伸至无穷，然而实践中不可能满足这样的条件。你很快就会看到，限制函数的定义域会妨碍从示例中完全恢复函数，除非在某些特殊情况下。

混叠（Aliasing）

在数字信号处理领域，混叠是指采样后不同信号变得无法区分的现象，换句话说，即一组信号伪装成了另一组信号。从概念出发，取样与混叠的关系不难理解。离散采样导致的混叠现象的前提是我们只能通过离散采样的方式来数字化描述连续函数。这意味着两个（或更多）完全不同的连续函数有可能在采样点处恰好重合，但我们无从得知非采样点处函数的特征。图$4.9$展示了以相同采样率采集的两个完全不同的正弦函数。

具有上述特征的两个连续函数被称为混叠对（aliased pair）。可以直观地看出，如果对采样进行细化，在采样信号中会越来越多地呈现出两个连续函数之间的差异。

二维傅里叶变换

双连续自变量 $t,z$ 的二维脉冲函数定义为：

$$\delta(t,z)=\left\{ \begin{aligned} 1\ \ \ \ & {\rm if\ }t=z=0\\ 0\ \ \ \ & {\rm otherwise\ } \end{aligned} \right. \tag{4-52}\label{4-52}$$

且其在二维空间的积分为1：

$$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\delta(t,z)dtdz=1 \tag{4-53}\label{4-53}$$

与一维情况下类似，二维脉冲函数仍然具有筛选特性：

$$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(t,z)\delta(t-t_0,z-z_0)dtdz=f(t_0,z_0) \tag{4-55}\label{4-55}$$

不妨令$f(t,z)$为包含两个连续变量$t,z$的二维连续函数，二维连续傅里叶变换由如下表达式给出：

$$F(\mu,\nu) = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}f(t,z)e^{-j2\pi(\mu t+\nu z)}dtdz \tag{4-59}\label{4-59}$$

相应的，二位连续傅里叶逆变换由如下表达式给出：

$$f(t,z) = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}F(\mu,\nu)e^{j2\pi(\mu t+\nu z)}d\mu d\nu \tag{4-60}\label{4-60}$$

其中，$\mu,\nu$为频域空间变量，而$t,z$通常理解为连续实空间变量。

与一维情况类似，二维离散采样也可以使用二维脉冲序列函数构建。二维脉冲序列函数表达式如下：

$$s_{\Delta T\Delta Z}(t,z)=\sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty}\delta(t-m\Delta T,z-n\Delta Z) \tag{4-61}\label{4-61}$$

其中，$\Delta T,\Delta Z$分别为沿$t$轴及$z$轴方向的采样间隔。

若二维函数$f(t,z)$的傅里叶变换在频率空间某个中心对称矩形区域外取值均为0，则我们称其为二维带限（band limited）函数，即：

$$F(\mu,\nu)=0,{\rm\ for}\ \forall\{(\mu,\nu)|\mu\geq\mu_{max}\vee\nu\geq\nu_{max}\} \tag{4-62}\label{4-62}$$

二维采样定理：二维连续带限函数$f(t,z)$可从它的一组离散样本中精准恢复原始连续函数，只需采样率满足：

$$\Delta T<\frac{1}{2\mu_{max}},{\rm\ and\ }\Delta Z<\frac{1}{2\nu_{max}} \tag{4-63}\label{4-63}$$

二维离散傅里叶变换及逆变换表达式如下：

$$F(u,v) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi(ux/M+vy/N)} \tag{4-67}\label{4-67}$$ $$f(x,y) = \frac{1}{MN}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi(ux/M+vy/N)} \tag{4-68}\label{4-68}$$

其中，$f(x,y)$为尺寸$M\times N$的数字图像，同一维情况类似，变量范围约束在$u=0,1,2,…,M-1$，且$v=0,1,2,…,N-1$。

有时你会在某些文献中发现系数$1/MN$出现在离散傅里叶正变换前，而非逆变换。又或是常数$1/\sqrt{MN}$同时出现在正变换和逆变换前，从而产生一个更对称的对。这些傅里叶变换对中的任意一种都是正确的，只要避免发生混用即可。

与一维情况类似，二维傅里叶变换仍然具有平移特性：

$$f(x,y)e^{j2\pi(u_0 x/M+v_0 y/N)}\Leftrightarrow F(u-u_0,v-v_0) \tag{4-71}\label{4-71}$$ $$f(x-x_0,y-y_0)\Leftrightarrow F(u,v)e^{-j2\pi(x_0 u/M+y_0 v/N)} \tag{4-72}\label{4-72}$$

同样的，二维傅里叶变换及其逆变换在$u,v$方向上是无限周期循环的，也就是说：

$$F(u,v)=F(u+k_1 M,v+k_2 N),{\rm\ and\ }f(x,y) = f(x+k_1 M,y+k_2 N),\ (k_1,k_2)\in\mathbb{Z}^1 \tag{4-74}\label{4-74}$$

在二维情况下，如图4.22(c)所示，围绕点$(M/2,N/2)$的区域可以划分为四个象限。在一维情况下，我们期望平移频域函数令$F(0)$取值在$M/2$处，二维情况下同理，我们通过平移频域函数令$F(0,0)$取值在点$(M/2,N/2)$处，即：

$$f(x,y)(-1)^{x+y}\Leftrightarrow F(u-M/2,v-N/2) \tag{4-76}\label{4-76}$$

二维傅里叶变换的特殊性质

泛函分析的一个重要结论是，任何实函数或复函数$w(x,y)$都可以表示为一个偶部和一个奇部的和，每一个都可以为实或复：

$$w(x,y)=w_e(x,y)+w_o(x,y) \tag{4-77}\label{4-77}$$

其中，偶部与奇部分别定义为：

$$w_e(x,y)\triangleq \frac{w(x,y)+w(-x,-y)}{2} \tag{4-78}\label{4-78}$$ $$w_o(x,y)\triangleq \frac{w(x,y)-w(-x,-y)}{2} \tag{4-79}\label{4-79}$$

根据奇偶函数的定义不难看出：

$$w_e(x,y) = w_e(-x,-y) \tag{4-80}\label{4-80}$$ $$w_o(x,y) = -w_o(-x,-y) \tag{4-81}\label{4-81}$$

偶函数是对称的，而奇函数是反对称的。由于DFT和IDFT中的所有指标均为非负整数，所以当我们讨论对称性时，我们指的是相对于中心点的对称性，此时奇偶函数的定义需满足：

$$w_e(x,y) = w_e(M-x,N-y) \tag{4-82}\label{4-82}$$ $$w_o(x,y) =-w_o(M-x,N-y) \tag{4-83}\label{4-83}$$

和通常一样，$x=0,1,2,…M-1, y=0,1,2,…,N-1$，$M,N$分别为二维数组的行列数。

从初等数学分析中我们知道两个偶函数或两个奇函数的乘积是偶函数，而一个偶函数和一个奇函数的乘积是奇函数。此外，一个离散函数可能是奇数的必要条件是它的所有样本之和为零。这些性质可以推导出一条重要的结论：

$$\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}w_e(x,y)w_o(x,y)=0 \tag{4-84}\label{4-84}$$

对于任意离散奇偶函数，均满足上式。也可以说，因$Eq.(\ref{4-84})$的辐角为奇，所以求和的结果为0。函数可以是实数也可以是复数。

利用前面的概念，我们可以建立DFT及其逆的一些重要的对称性质。一个常用的性质是实函数$f(x,y)$的傅里叶变换是共轭对称（Conjugate symmetry）的：

$$F^\ast(u,v) = F(-u,-v) \tag{4-85}\label{4-85}$$

证明如下：

$$\begin{aligned} F^\ast(u,v) & = \Bigg[\sum_{x=0}^{M-1}\sum_{u=0}^{N-1}f(x,y)e^{-j2\pi(ux/M+vy/N)}\Bigg]^\ast\\ & \\ & = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f^\ast(x,y)e^{j2\pi(ux/M+vy/N)}\\ & \\ & = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi([-u]x/M+[-v]y/N)}\\ & \\ & = F(-u,-v) \end{aligned}$$

其中，第二步源自于$(a\cdot b)^\ast = a^\ast\cdot b^\ast$，而第三步则因为$f(x,y)$为实函数。若$f(x,y)$为虚函数，则其傅里叶变换满足共轭反对称，即$F^\ast(-u,-v)=-F(u,v)$。

傅里叶频谱及相位角（Spectrum and Phase Angle）

考虑到二维离散傅里叶变换的结果一般为复函数，因此可以使用极坐标形式描述：

$$F(u,v) = R(u,v)+jI(u,v) = |F(u,v)|e^{j\phi(u,v)} \tag{4-86}\label{4-86}$$

其中，

$$|F(u,v)| = [R^2(u,v)+I^2(u,v)]^{1/2} \tag{4-87}\label{4-87}$$

称作傅里叶变换频谱（Fourier or frequency spectrum），并且，

$$\phi(u,v) = arctan\Bigg[\frac{I(u,v)}{R(u,v)}\Bigg] \tag{4-88}\label{4-88}$$

称作傅里叶变换相位角或相位谱（phase angle or phase spectrum），其中 arctan 务必使用四象限反正切函数（fourquadrant arctangent function）计算，例如 Matlab 中atan2(Imag,Real)。

最后的，能量谱（power spectrum）被定义为：

$$P(u,v) = |F(u,v)|^2 = R^2(u,v)+I^2(u,v) \tag{4-89}\label{4-89}$$

如前所述，$R(u,v)$和$I(u,v)$分别为$F(u,v)$的实部和虚部，所有的计算都是针对离散变量进行的，$u = 0,1,2,…,M-1$，$v=0,1,2,…,N-1$。因此，频率谱、相位谱、能量谱均为尺寸$M\times N$的矩阵。

实函数的傅里叶变换是共轭对称的[见$Eq.\ref{4-85}$]，这意味着谱相对于原点具有偶对称性：

$$|F(u,v)| = |F(-u,-v)| \tag{4-90}\label{4-90}$$

而相位谱相对于原点遵从奇对偶性：

$$\phi(u,v) = -\phi(-u,-v) \tag{4-91}\label{4-91}$$

根据$Eq.\ref{4-67}$，

$$F(0,0) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y) = MN\overline{f}$$

在电气行业中，有时又称 $F(0,0)$ 为傅里叶变换的直流分量（dc component, where “dc” signifies direct current）。

Figure 4.23(a) 展示了一幅矩形图像，而 Figure 4.23(b) 中展示了其频谱，各点取值被缩放至[0,255]范围内，并以图片形式展示了出来。实空间与频率空间坐标原点均位于左上角，如图中所示。正如预期的那样，频谱原点附近区域包含了最大值，因而在图像中显得更亮。然而请注意，频谱的四角均包含类似的极值，原因即是上节讨论过的周期性。为了将频谱对中心，我们只需在对原始图像进行离散傅里叶变换前乘以$(-1)^{x+y}$即可。Figure 4.23(c) 展示了对中心后的频谱，显然对中心后的频谱更适合于可视化的观察判断。由于直流分量主导了频谱峰值，压缩了频谱图像中其它区域强度的动态范围。为了突出这些细节，我们使用了对数变换。Figure 4.23(d) 展示了经公式 $log(1+|F(u,v)|)$ 对数变换后的结果。对数变换对细节的增强作用是明显的，在后续段落中展示的大部分频谱均将采用这种方式进行缩放。

根据$Eq.\ref{4-72}$，频谱对图像的平移变换不敏感（指数项绝对值为1），而敏感于图像的旋转变换，且频谱旋转角度与原始图像旋转角度相同。Figure 4.24 展示了这一特性。Figure 4.24(b) 与 Figure 4.23(d) 中的频谱图像完全一致。然而，他们各自的实空间图像显然不同，则它们的相位谱必然存在差异，如Figure 4.25。Figure 4.25(a,b) 分别为 Figure 4.23(a)与 Figure 4.24(a) 的相位谱。注意相位图像之间非常缺乏相似性，尽管事实上，它们所对应的实空间图像之间的唯一区别是简单的平移。一般来说，相位谱的视觉分析能够获取的直观信息很少。例如，由于其呈 $45^\circ$ 朝向，人们容易误以为Figure 4.25(a) 的相位谱对应于旋转后的实空间图像 Figure 4.24(c)，而不是真正的 Figure 4.23(a)。然而事实上，Figure 4.25(c) 才是真正的旋转后图像所对应的相位谱。

二维离散卷积（2D Discrete Convolution Thorem）

拓展$Eq.\ref{4-48}$至二维情况，可以得到如下二维离散卷积表达式：

$$(f\otimes h)(x,y) = \sum_{m=0}^{M-1}\sum_{n=0}^{N-1}f(m,n)h(x-m,y-n) \tag{4-94}\label{4-94}$$

其中，$x=0,1,2,…,M-1$，$y=0,1,2,…,N-1$。而二维卷积定理则表示为：

$$(f\otimes h)(x,y)\Leftrightarrow (F\cdot H)(u,v) \tag{4-95}\label{4-95}$$ $$(f\cdot h)(x,y)\Leftrightarrow \frac{1}{MN}(F\otimes H)(u,v) \tag{4-96}\label{4-96}$$

正如前文所言，空间卷积是空间滤波的核心，因此$Eq.\ref{4-95}$是构建空间卷积与频域滤波间等价性的重要纽带，是图像频域滤波的理论基石。

强度变换及滤波

  空间域（spatial domain）指图像平面本身，针对空间域的图像处理方法均为对图像平面中各像素的直接处理，这区别于针对频域（frequency domain）的图像处理方法。考虑一个一维函数，假设我们希望消除其在截止值以上的所有频率，同时保留该值以下的所有频率。图$3.38(a)$展示了实现此目标的低通频域滤波函数。术语滤波传递函数（filter transfer function）通常用于表示频域的滤波器函数，类似于将空间域中滤波器函数称作滤波卷积核（filter kernel）。然而，这样的过滤功能是无法通过物理组件实现的，且在数字化实现时还会出现“振铃”现象（ringing）。我们通常使用的所有频域滤波器都是围绕原点对称的，包含正频率和负频率。空间域滤波的等效卷积核是频域中滤波传递函数（乘子）的傅里叶反变换。

  对尺寸为$M\times N$的图像在空间域进行滤波变换，卷积核尺寸为$m\times n$，该操作的时间复杂度为$O(MNmn)$。若该卷积核是可拆分的，则时间复杂度简化至$O(MN(m+n))$。然而通过频域进行完全等价的滤波操作所需的时间复杂度为$O(2MNlog_{2}(MN))$，其中2倍来自于傅里叶变换与傅里叶逆变换两步时间复杂度相同的操作。

Codes and Data

Estimating spatial resolution of real image[1]: [尚未完成]

References

[1]  Mizutani, Ryuta, Rino Saiga, Susumu Takekoshi, Chie Inomoto, Naoya Nakamura, Masanari Itokawa, Makoto Arai et al. A method for estimating spatial resolution of real image in the Fourier domain. Journal of microscopy 261, no.1 (2016): 57-66.