前言

在这篇文章里，其实有关图像平移并不是最重要的，最重要的是其中的一些思想。无论你是不是需要用到图像平移模型，都可以看一下这篇文章，相信你会有所收获的。

傅里叶变换

从傅里叶级数到傅里叶变换

在高数中，我们学过傅里叶级数，傅里叶级数向我们说明了一个很重要的事情，那就是一个周期函数可以分解成无数个正弦和余弦函数的和，并且，这个正弦和余弦函数的频率都是原周期函数周期的整数倍，也就是下面这个公式：

$$
f(x)=a_0+\sum_{n=1}^{\infty}\left(a_n \cos \left(\frac{2 \pi n}{T} x\right)+b_n \sin \left(\frac{2 \pi n}{T} x\right)\right)
$$

其中的常数$a_0$对应着零频率。

但是，在现实生活中，很多的信号往往是没有周期性的，那么能不能对非周期函数进行展开呢？这也就是对上面的级数进行推广，得到了傅里叶变换。因此，傅里叶级数面对的是周期信号，而傅里叶变换面对的是非周期信号。

傅里叶变换的表达式是可以从傅里叶级数的理论推广而来的，具体的方法可以去看网上的很多讲解，我这里简单阐述一下我的理解。既然我们是想将傅里叶级数推广到非周期信号，那么能不能把非周期信号也当成周期信号呢？这时可以将非周期的信号的周期看成$\infty$，即$T\rightarrow \infty$。我们知道，傅里叶级数中的正弦和余弦函数的频率是离散的，而当周期趋于$\infty$时，我们可以很容易想到，这个频率将会变成连续的，这就是其中的一个推广：从离散到连续，反映在计算式上，就是从级数到积分。但是这个过程中，如果直接将上式写成连续的形式，那么将几乎无法求解。因此，能不能将上面这个式子中的 cos 和 sin 表示成一个形式，也就是将系数 a 和 b 变成一个，这就引入了傅里叶级数的复数表达。但是请注意，复数表达和上面的级数形式并无大的区别。当你将级数表达为复数形式后，然后就可以很自然的将其中的频率变成$d\omega$，就得到了傅里叶变换的表达式：

$$
F(\omega)=\frac{1}{2 \pi} \int_{-\infty}^{\infty} f(x) e^{-i \omega x} d x
$$

这也可以看作是一种推广，也就是从实数到复数。掌握了这两个关键点，基本上就清楚了怎么从傅里叶级数到傅里叶变换了。然后也可以得到逆傅里叶变换：

$$
f(t)=\operatorname{IFT}(F(w))=\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(w) e^{j w t} d w
$$

离散傅里叶变换

现实中，常常需要处理离散信号，比如在通信领域，传输信号时，需要将连续信号进行采样编码，变成数字信号，这样可以大幅减少噪声的影响。还有在本节内容中的图像数据，也是由一个个的像素点组成，构成了一个离散的数据。那么如何从连续傅里叶变换 (CFT) 到离散傅里叶变换 (DFT) 呢

思想其实也很简单，从连续信号到离散信号，是通过采样得到的，那么这个变换，也是利用采样。这个过程在这里也不详细论述了，大概过程就是，先对时域采样，CFT 的表达式的右边由积分变成了离散求和，这时时域是离散的，频域当然也需要离散，那么再对频域采样，就得到了最后的结果。当然，这个过程中还有很多细节，比如需要满足采样定理，还有一些加窗处理等等。

经过这样一番处理，对频域进行采样后，得到的这一些$\omega_k$就不能成为频率了，但是为了和 CFT 对应，称之为数字频率。DFT 的表达式如下：

$$
F(w)=\sum_{t=0}^{N-1} f(t) e^{-j 2 \pi \frac{w t}{N}}, w=0,1, \cdots, N-1
$$

这个是一维的 DFT，为了在图像上使用，我们可以将其推广到二维 DFT，表达式如下：

$$
\begin{aligned}
F(u, v)=\sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) e^{-j 2 \pi\left(\frac{u x}{M}+\frac{vy}{N}\right)},
u=0,1, \cdots, M-1 ; v=0,1, \cdots, N-1
\end{aligned}
$$

DFT 的矩阵表达

再看这个表达式：

$$
F(w)=\sum_{t=0}^{N-1} f(t) e^{-j 2 \pi \frac{w t}{N}}, w=0,1, \cdots, N-1
$$

对于每一个数字频率，都是有 N 项求和，将$e^{-j 2 \pi \frac{w t}{N}}$成为基函数，那么对于每一个频率$\omega$，都对应着 N 个基函数，

$$
\begin{array}{ll}
w=0 & 1,1, \cdots, 1 \\
w=1 & 1, e^{-j 2 \pi \frac{1}{N}}, \cdots, e^{-j 2 \pi \frac{N-1}{N}} \\
\vdots \\ \omega=N-1 & 1, e^{-j 2 \pi \frac{N-1}{N}}, \cdots, e^{-j 2 \pi \frac{(N-1)(N-1)}{N}}
\end{array}
$$

可以发现，每一行的相邻元素的相位差是恒定的。这是一个很重要的特征，请记住这个特征，后面会考的。

我们可以将其写作一个$N\times N$的矩阵，然后写成矩阵表达：

$$
\left[\begin{array}{c}
F(0) \\
F(1) \\
\vdots \\
F(N-1)
\end{array}\right]=\left[\begin{array}{cccc}
1 & 1 & \cdots & 1 \\
1 & e^{-j 2 \pi \frac{1}{N}} & \cdots & e^{-j 2 \pi \frac{N-1}{N}} \\
\vdots & \vdots & \ddots & \vdots \\
1 & e^{-j 2 \pi \frac{N-1}{N}} & \cdots & e^{-j 2 \pi \frac{(N-1)(N-1)}{N}}
\end{array}\right]\left[\begin{array}{c}
f(0) \\
f(1) \\
\vdots \\
f(N-1)
\end{array}\right]
$$

这就是一维 DFT 的矩阵表达形式。

同样的，我们也可以写出二维 DFT 的矩阵表达：

$$
\begin{aligned}
& {\left[\begin{array}{cccc}
F(0,0) & F(0,1) & \cdots & F(0, N-1) \\
F(1,0) & F(1,1) & \cdots & F(1, N-1) \\
\vdots & \vdots & \ddots & \vdots \\
F(M-1,0) & F(M-1,1) & \cdots & F(M-1, N-1)
\end{array}\right]=\left[\begin{array}{cccc}
1 & 1 & \cdots & 1 \\
1 & e^{-j 2 \pi \frac{1}{M}} & \cdots & e^{-j 2 \pi \frac{M-1}{M}} \\
\vdots & \vdots & \ddots & \vdots \\
1 & e^{-j 2 \pi \frac{M-1}{M}} & \cdots & e^{-j 2 \pi \frac{(M-1)(M-1)}{M}}
\end{array}\right]} \\
& {\left[\begin{array}{cccc}
f(0,0) & f(0,1) & \cdots & f(0, N-1) \\
f(1,0) & f(1,1) & \cdots & f(1, N-1) \\
\vdots & \vdots & \ddots & \vdots \\
f(M-1,0) & f(M-1,1) & \cdots & f(M-1, N-1)
\end{array}\right]\left[\begin{array}{cccc}
1 & 1 & \cdots & 1 \\
1 & e^{-j 2 \pi \frac{1}{N}} & \cdots & e^{-j 2 \pi \frac{N-1}{N}} \\
\vdots & \vdots & \ddots & \vdots \\
1 & e^{-j 2 \pi \frac{N-1}{N}} & \cdots & e^{-j 2 \pi \frac{(N-1)(N-1)}{N}}
\end{array}\right]}
\end{aligned}
$$

可以发现，二维 DFT 分别是左乘和右乘了相应的 DFT 矩阵，更高维可能就要用到张量了。

DFT 求解图像平移问题

那么为什么要引入傅里叶变换以及频域呢？其目的其实很简单，那就是利于我们解决问题，那么体现在哪里呢？其中的一个体现便是 FT 很好的几个性质，一个是时域平移等于频域相移，另一个是时域卷积等于频域乘积，其他的性质我了解的就不多了。与时域不同，频域是一个虚构的空间，在很多领域，都有这种思想，设计一个虚构的空间，然后将问题放到这个空间中进行研究，可以大幅降低难度，并且能够帮助我们更好的理解问题。比如在半导体领域，引入能带模型，那么面对电子被激发时，我们可以从晶格模型出发，一个电子获得能量，挣脱了共价键的舒服，变成了自由电子，也可以从能带模型出发，电子获得足够的能量，从价带跃迁到了导带。

那么对于图像平移来说，就对应着时域的平移。FT 的时移性质：

$$
\begin{aligned}
f(t)=g(t-\tau) \rightarrow
F(\omega)=G(\omega) e^{-j 2 \pi \frac{\omega \tau}{N}}
\end{aligned}
$$

可以发现，时域移动了$\tau$，频率相位移动了$2\pi \omega \tau /N$。同样的，对于二维 DFT：

$$
\begin{aligned}
f(m, n)=g\left(m-m_0, n-n_0\right) \rightarrow
F(u, v)=G(u, v) e^{-j 2 \pi\left(\frac{u m_0}{M}+\frac{v n_0}{N}\right)}
\end{aligned}
$$

那么，对于图像平移问题，似乎问题就变得很简单了。我们可以找到或者构造一个表达式，将这个相位移动的量提取出来，就可以解决问题了。这个量就可以是归一化互功率谱 (NCPS)。

$$
S(u, v)=\frac{F(u, v) G^*(u, v)}{\left|F(u, v) G^*(u, v)\right|}
$$

我们将上面时移的表达式代入进去，可以求得时移前和时移后的 NCPS 表达式为：

$$
S(u, v)=e^{-j 2 \pi\left(\frac{u m_0}{M}+\frac{v n_0}{N}\right)}
$$

这就将相位的移动量提取了出来，我们可以对其进行 IDFT：

$$
\begin{aligned}
s(m, n)=I D F T[S(u, v)]=\delta\left(m-m_0, n-n_0\right)
\end{aligned}
$$

我们将 IDFT 的结果，找到最大值所在的位置，就对应着平移量。OK，问题似乎是解决了。

但是，这种方法求出的结果一定是一个整数，也就是图像平移的整数量。那么能不能实现任意精度亚像元平移呢？比如说平移 3.5 个像元。对于任意精度来说，这种方法似乎就行不通了，我们需要寻求新的方法。

循环移位矩阵

图像数据本质上就是一个矩阵数据，那么图像的平移，就正好对应着矩阵的平移，那么这就对应着移位操作。循环移位矩阵的形式如下：

$$
\mathbf{Q}=\left[\begin{array}{ccccc}
0 & 1 & 0 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
\vdots & \ddots & \ddots & \ddots & 0 \\
0 & \ddots & \ddots & \ddots & 1 \\
1 & 0 & \cdots & \cdots & 0
\end{array}\right]
$$

那么，对于一个矩阵$\mathbf{A}$，如果我们需要将其向右平移两个像元，向下平移两个像元，就可以写成$(Q^T)^2AQ^2$.

这样的话，如果我们想实现小数级的平移，就涉及到一个问题，这个矩阵的任意次方是不是实数？如果是的话，也就意味着可以平移任意次。

对于这个问题，可以从线性变换的角度来看待。我们知道特征值对应着不同的线性变换（这也是特征值反映了矩阵的特征的原因之一）。假如这个矩阵有一个实特征值，如果这个特征值大于 0，那么就对应着伸缩变换，如果这个特征值小于 0，就对应着反射变换，也就是转了 180 度，如果这个矩阵有复特征值，就对应着旋转变换。理解了这个，这个问题就变得很简单了，对矩阵进行任意次方，就可以看作是线性变换了任意次。伸缩变换在变换了任意次之后，依然是伸缩变换，因此，这个可以任意次方。反射变换在变换了任意次之后，并不一定是反射变换，这是一个不连续的变换，因此，反射变换并不能开任意次方。旋转变换是一个连续变换，因此可以开任意次方。

好了，那么结论就已经出来了，如果矩阵包含负的实特征值，那么就不一定能开任意次方（注意这里是不一定，因为如果一个矩阵包含两个负特征值，那么反射两次之后，宏观上表现出了伸缩变换，那么可以开任意次方），其他情况下就可以开任意次方。

这里，我们再来看循环移位矩阵，首先来看它的实特征值。

实特征值

$$
\begin{aligned}
\mathbf{u} & =\left[\begin{array}{lllll}
u_1 & u_2 & \cdots & u_{n-1} & u_n
\end{array}\right]^{\mathrm{T}} \\
\mathbf{Q u} & =\left[\begin{array}{lllll}
u_2 & u_3 & \cdots & u_n & u_1
\end{array}\right]^{\mathrm{T}}
\end{aligned}
$$

如果$\lambda=1$，可以得到：

$$
\mathbf{Q}\mathbf{u}=\mathbf{u}
$$

可以得到

$$
u_1=u_2=\cdots=u_n
$$

那么特征向量变为$[1,1,\cdots,1]$。可以发现，$\lambda=1$是特征值。

如果$\lambda=-1$，可以很容易得到，如果矩阵是偶数阶，可以得到特征向量$[1,-1,\cdots,1,-1]$。如果是奇数阶，则 -1 不是特征值。

因此，我们可以得到结论，如果循环移位矩阵的阶数为偶数，包含特征值 -1，显然不能开任意次方，而如果是奇数阶，则可以开任意次方。因此，我们一般要求图像数据是奇数阶，如果不是，则可以考虑去掉一行或者一列，变成奇数阶。

因此，如果图像数据是奇数阶，我们便可以利用循环移位矩阵，来对图像进行任意精度的平移。

复特征值

观察这个式子：

$$
u_2=\lambda u_1,u_3=\lambda u_2,u_1=\lambda u_n
$$

可以发现，这些相邻$u_{k}$的相位差是恒定的，都是相差了$\lambda$的相位。并且，可以发现，从$u_2$到$u_n$再到$u_1$，相位不断叠加，但最后回到了起点，这就相当于一个圆，即这些相位不断叠加，回到起点时，正好经历了$2\pi$或是它的整数倍。我们可以将$\lambda$的相位写成这个形式：$2\pi \omega /n$，其中$\omega$为整数。到这里，有没有发现这个有点似曾相识！

没错，这个正好是 DFT 矩阵中每一行的相位差。因此，我们可以得到，循环移位矩阵的特征向量矩阵就是 DFT 中的矩阵。

图像平移优化模型

通过上面的分析，我们知道了如果图像数据是奇数阶，就可以利用循环移位矩阵对其进行任意精度的平移。因此，对于求解平移量的问题，我们可以构造如下优化模型：

$$
\min _{x, y} f(x, y)=\min _{x, y}\left|\mathbf{A}-\left(\mathbf{Q}_1^{\mathrm{T}}\right)^x \mathbf{B} \mathbf{Q}_2^y\right|_{\mathrm{F}}
$$

其中，$x,y$就是我们需要求解的平移量。那么应该如何进行求解呢？我们对里面的式子进行一步步的化简。首先我们对循环移位矩阵进行展开：

$$
\mathbf{Q}=\mathbf{U D U}^{\mathrm{H}}
$$

接下来的推导过程如下图：

然后，我们就得到了最后的解：

$$
\mathbf{s}=\left(\mathbf{X}^T \mathbf{X}\right)^{-1} \mathbf{X}^T \mathbf{c}
$$