As-Projective-As-Possible Image Stitching with Moving DLT
原文链接：As-Projective-As-Possible Image Stitching with Moving DLT
这里主要对文章的 Moving DLT 做介绍，其它关于图像拼接的细节就不描述了，文章的重点就在使用 Moving DLT 来模拟一个复杂映射
摘要：图像拼接的商用化给我们带来了这个问题已经被解决的错误印象。实际上，在某些严厉的成像条件上，这些结果都不能够使人信服，除非在只有视角旋转变化或者场景为平面的条件下。在这些约束下，可以使用二维射影变换或者单应来对齐图片。在一般情况下，这种条件经常被打破，造成错误的对齐或者鬼影。现在这些拼接技术依赖于后处理方法来消除或减轻鬼影。在这篇文章中，提出一种不同的估计技术，名为 “移动直接线性变换 (Moving DLT)”，可以优化射影变换使得其适应输入数据和理想条件的偏差。这生成了一种 尽可能投影的 的图像对齐方法，有效的减少了鬼影并且不损害透视图象拼接的几何写实性。此外，还描述了如何将这种对齐方法使用到多幅图像的大全景图生成中，通过将光束法平差（Bundle Adjustment）同时用于优化以达到准确地对齐。
关键词：图像拼接，图像对齐，射影变换，直接线性变换，移动最小二乘法
介绍
文章的介绍部分主要谈到了当前图像拼接算法模型的不足，使用单应变换来实现图像对齐。但单应变换有着自身的局限性，只适用于视角仅旋转变换和目标场景为平面两种情况。当违反这两个条件后，会造成鬼影并且错误的对齐。为此，提出了一个名为“移动DLT”的估计技术，可以优化全局单应变换带来的误差。其它技术带来的鬼影等不良影响虽然可以通过后期技术优化，但是作者认为在对齐阶段就消除是很有必要的。
尽可能射影的变换 APAP
二维射影变换（DLT）
令 \(\textbf{x} = [x,y]^T,\textbf{x}'= [x',y']^T\)，为两幅具有重叠部分图像 \(I,I'\) 上的对应点，它们之间的关系为： \[ \widetilde {\textbf{x}}' \sim \textbf{H} \widetilde {\textbf{x}} \] 字符上的波浪线代表着齐次坐标归一化。也就是说将 \(\widetilde {\textbf{x}}\) 通过单应矩阵变换后还需要齐次坐标归一化，也就是除第三个分量。其中 \(3\times 3\) 的矩阵 \(\textbf{H}\) 为单应矩阵，上面的等式展开后有： \[ \begin{align} x' = \frac{r_1[x,y,1]^T}{r_3[x,y,1]^T} \\ y' = \frac{r_2[x,y,1]^T}{r_3[x,y,1]^T} \end{align} \] 其中 \(r_j\) 表示矩阵 \(\textbf{H}\) 的第 j 行。齐次坐标需要除以第三个分量使得它不是一个线性变换，直接线性变换（DLT）是一种在含有噪声的匹配点中估计单应矩阵的基础方法。首先，因为 \(\widetilde {\textbf{x}}' \sim \textbf{H} \widetilde {\textbf{x}}\) 有 \(\textbf{0}_{3\times 1} = \widetilde {\textbf{x}}' \times \textbf{H}\widetilde {\textbf{x}}\) ，这个式子可以化为： \[ \begin{align} \textbf{0}_{3\times 1} &= \left[ \begin{matrix} \textbf{0}_{1\times 3} & -\widetilde{\textbf{x}}^T & y'\widetilde {\textbf{x}}^T\\ \widetilde {\textbf{x}}^T & \textbf{0}_{1\times 3} & -x'\widetilde {\textbf{x}}^T\\ -y'\widetilde {\textbf{x}}^T & x'\widetilde {\textbf{x}}^T & \textbf{0}_{1\times 3} \end{matrix} \right] \times \textbf{h} \\ \textbf{h} &= \left[ \begin{matrix} r_1^T \\ r_2^T \\ r_3^T \end{matrix} \right] \end{align} \] \(\textbf{h}\) 就是矩阵 \(\textbf{H}\) 的向量表示，上面 3*9 的矩阵中，只有两行是线性无关的，所以对某一对对应点 \(\{\textbf{x}_i,\textbf{x}_i'\}\) 使用 \(\textbf{a}_i\) 表示矩阵的前两行，对于需要估计的向量 \(\textbf{h}\)，\(\vert\vert \textbf{a}_i \textbf{h}\vert\vert\) 为第 i 个点的代数误差（因为理想情况下 \(a_ih\) 为 0），所以 DLT 的目标就是将所有对应点的代数误差和做到最小，即： \[ \hat{\textbf{h}} = \mathop{argmin}_{\textbf{h}} \sum_{i=1}^{N} \vert\vert \textbf{a}_i \textbf{h}\vert\vert ^2 \qquad \text{s.t.} \vert\vert \textbf{h}\vert\vert = 1 \] 将所有的点累积起来还可以转变为矩阵的形式： \[ \hat{\textbf{h}} = \mathop{argmin}_{\textbf{h}} \vert\vert \textbf{A} \textbf{h}\vert\vert ^2 \qquad \text{s.t.} \vert\vert \textbf{h}\vert\vert = 1 \] 其中 \(\textbf{A} \in \mathbb{R}^{2N\times 9}\)，这就是一个代数最小二乘问题，它的解是最小奇异值对应的右奇异向量。也就是说对于原图像 \(I\) 中的任意一点 \(\textbf{x}_*\)，它映射到对齐图像 \(I'\) 上的点为 \(\textbf{x}'_*\)，理想情况下，有以下关系成立： \[ \widetilde {\textbf{x}}_*' \sim \textbf{H} \widetilde {\textbf{x}}_* \] 前面已经提到了单应的局限性，DLT 算法适合在避免了单应局限性，并存在一些离群点（噪声）的情况下，获得最佳的单应映射。但是在一般成像下，一个全局的单应不足以拟合所有的点，也就是说会有很大一部分点不符合这个单应矩阵表示的映射。文章提供了一种使用单应模拟复杂映射（非纯单应）的方法，也就是下面要提到的移动 DLT。
在做 DLT 之前可以先将数据归一化可以避免很多数值上的错误，归一化可参考：《CV中的多视图几何》——图像之间的单应变换 - 知乎 (zhihu.com)，其中也谈到了归一化后的反归一化。
Moving DLT
Moving DLT 本质上就是使用局部单应拟合映射。此时的点 \(\textbf{x}_*,\textbf{x}_*'\) 满足下式： \[ \widetilde {\textbf{x}}_*' \sim \textbf{H}_* \widetilde {\textbf{x}}_* \] 这里的 \(\textbf{H}_*\) 和点有关，也就是单应的局部表示。而 \(\textbf{H}_*\) 则从下面的加权优化问题中得出： \[ \textbf{h}_* = \mathop{argmin}_{\textbf{h}} \sum_{i=1}^{N}\vert\vert w_*^i \textbf{a}_i \textbf{h}\vert\vert ^2 \qquad \text{s.t.} \vert\vert \textbf{h}\vert\vert = 1 \] 其中的 \(w_*^i\) 是最关键的部分，它给了那些与 \(\textbf{x}_*\) 相近的点更高的权重。可以使用高斯函数来确定： \[ w_*^i = \exp(-\vert\vert \textbf{x}_* - \textbf{x}_i \vert\vert^2/\sigma^2) \] 自然的 \(\sigma\) 表示尺度参数。和 WLS 的权函数类似，只要符合权函数那几个定义即可。这差不多就是 Moving DLT 的定义了。就加了一个权函数。上式可以矩阵化为： \[ \textbf{h}_* = \mathop{argmin}_{\textbf{h}} \vert\vert \textbf{W}_* \textbf{A} \textbf{h}\vert\vert ^2 \qquad \text{s.t.} \vert\vert \textbf{h}\vert\vert = 1 \] 其中的 \(\textbf{W}_*\) 为权重组成的对角矩阵： \[ \textbf{W}_* = \mathop{diag}([w_*^1 w_*^1 w_*^2 w_*^2 \cdots w_*^N w_*^N]) \] 每一个权重都出现两次是因为 \(\textbf{a}_i\) 是一个 2 行的矩阵。矩阵化后就变成了一个加权 SVD（WSVD） 问题，它的解为 \(\textbf{W}_*\textbf{A}\) 最小奇异值的右奇异向量。为了确保有界，对那些远离的点做一个补偿（正则化），当权值低于某个值时，将它设为预定的值，因为有些时候没有匹配点（不再图像的重叠区域内），这时那些有匹配点原理当前点，则会数据不足。所以改进后的权值表达如下： \[ w_*^i = \max(\exp(-\vert\vert \textbf{x}_* - \textbf{x}_i \vert\vert^2/\sigma^2), \gamma) \] 概念上讲 Moving DLT 是 MLS（移动最小二乘）的单应变换版本。
图像拼接的高效计算（复杂度优化）
RANSAC 过滤离群点
数据一般都不会是理想的，总存在噪声或者错误匹配的点（离群值）。文章使用 RANSAC 来移除离群值来减少不必要的数据。RANSAC 还是使用的 DLT，这就使得一些正确匹配点（单应局限性造成）和这些离群值一样会在 RANSAC 中产生误差，但离群点的误差会比单应局限性造成的误差大，所以 RANSAC 还是可以有效的工作。
划分网格
为每一个源图像中的点计算它的单应模型是没有必要的，因为相近的点的单应差别不大，所以讲图像均分为网格可以提高效率并且保证准确度。可以将原图像划分为 \(C_1\times C_2\) 大小的网格，每一个网格的单应由它的中心点 \(\textbf{x}_*\) 代表。这样会使得格子之间的单应变换不连续，因为使用离散的网格中心点导致权值不是连续变化。但作者在实验中认为只要分辨率足够高就不会有明显的不连续感。
WSVD 更新
在计算 Moving DLT 的时候每一次改变的只有权重矩阵，但是权重矩阵中大部分都是 \(\gamma\)（因为里当前点远的点更多），现在考虑是否不用每次都重新计算当前点的 DLT，而是通过上一次的结果来更新。
首先考虑权重矩阵 \(\textbf{W}_\gamma = \gamma \textbf{I}\)，令 \(\textbf{V}\) 为矩阵 \(\textbf{W}_\gamma\textbf{A}\) 右奇异向量，这个矩阵的特征值分解如下： \[ \textbf{A}^T\textbf{W}_\gamma^T \textbf{W}_\gamma\textbf{A} = \textbf{V}\textbf{D}\textbf{V}^T \] 令 \(\tilde {\textbf{W}}\) 为 \(\textbf{W}_\gamma\) 改变第 i 个对角元素得到的矩阵。对它进行特征值分解有： \[ \begin{align} \textbf{A}^T\tilde {\textbf{W}}^T \tilde {\textbf{W}}\textbf{A} &= \textbf{V}\textbf{D}\textbf{V}^T + \rho \textbf{r}_i \textbf{r}_i^T = \textbf{V}(\textbf{D} + \rho \overline{\textbf{r}}_i \overline{\textbf{r}}_i^T)\textbf{V}^T \\ \rho &= (\tilde{w}^2_i/\gamma^2 - 1) \end{align} \] 其中 \(\textbf{r}_i\) 为矩阵 \(\textbf{A}\) 的第 i 行，\(\overline{\textbf{r}}_i = \textbf{V}^T\textbf{r}_i\)。 \(\textbf{D} + \rho \textbf{r}_i \textbf{r}_i^T\) 又可以特征值分解为 \(\tilde{\textbf{C}}\tilde{\textbf{D}}\tilde{\textbf{C}}^T\)，这个可以使用特征方程得到。解特征方程可以使用相关方法加速。最后解为 \(\textbf{V}\tilde{\textbf{C}}\)，这个称为 Rank-One Update，可参考特征值分解（EVD）的rank 1 update - 知乎 (zhihu.com)。
参考图如下：
其中 a 为原图，b 为将原图网格化，c 为将原图通过 Moving DLT 映射到待匹配平面的结果。
同时优化
这里讲的是将 Moving DLT 用到多幅图像（大于三幅）同时拼接用于产生全景图时，产生的相关问题和解决方法。
这里就没有怎么了解了。
主要是拼接顺序的选择以及光束法平差优化（Bundle Adjustment）两个方面的内容。