大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE使用 1年只要46元 售后保障 童叟无欺

图像拼接首先通过一些手段（标定、SIFT等特征点、其他传感器等）获取图像间的对应关系（2D-2D）。

这个对应关系可以用单应性矩阵$H$（也称透视变换）描述
$$\begin{array}{rl}\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ w^{\prime }\end{array}\right]=& \left[\begin{array}{ccc}{h}_{00}& h_{01}& h_{02}\\ h_{10}& h_{11}& h_{12}\\ h_{20}& h_{21}& h_{22}\end{array}\right]\left[\begin{array}{c}x\\ y\\ w\end{array}\right]\\ =& H\left[\begin{array}{c}x\\ y\\ w\end{array}\right]\end{array}\text{\hspace{1em}(1)}$$

但上述思路的前提是：待拼接的目标处于或者近似处于同一平面。否则就会产生遮挡和角度差异，从而造成拼接图的“虚影”

                                                 图1. 由于相机角度差异造成的虚影现象（红圈）

针对上述问题，Autostitch（OpenCV实现的方法）和Photosynth等方法走的是一个全局H + 虚影后处理的路线；
APAP【1】算法则使用多个局部H的方式尝试在拼接时就尽可能精确。

APAP算法思路

APAP算法官方主页提供了matlab的源码
Github上有第三方实现的C++代码

算法以A图为基准，首先计算B图到A图的全局变换$H$（DLT方法）,
然后将B图等间距划分小网格，逐个计算B图每个小网格与A图的局部变换$H_{\ast }$（Moving DLT方法）,
最后将B图每个网格内的点按对应的$H_{\ast }$做局部变换，与A图融合。

                                                        图2. 网格局部变换（摘自APAP论文）

几个问题

SVD在最小二乘中的应用

参考博文《SVD之最小二乘【推导与证明】》

– 模型$Ax=b$的求解
令$A$尺寸为$m\times n$，且$m>n$（超定问题）
则优化目标为$min||Ax-b||$

首先对$A$做SVD分解，有
$$A=U\Sigma V^T\text{\hspace{1em}(2)}$$
其中$U$和$V$分别称为$A$的左、右奇异向量，满足正交性（即$U^T=U^{-1}$）
$\Sigma$为对角阵，对角元素是$A$的奇异值

将式(2)代入目标函数，有
$$\begin{array}{rl}& min||Ax-b||\\ & \to min||U\Sigma V^{T}x-b||\\ & \to min||\Sigma V^{T}x-U^{T}b||(正交性)\\ & \to min||\Sigma y-b^{\prime }||(令y=V^{T}x,b^{\prime }=U^{T}b)\end{array}\text{\hspace{1em}(3)}$$
注意到$\Sigma$是一个对角阵，令${\lambda }_i$为其对角线元素，则有
$$\begin{gathered} y_i=\frac{b_i^{\prime }}{{\lambda }_i} \\ x=Vy\text{\hspace{1em}(4)} \end{gathered}$$

– 模型$Ax=0$的求解
令$A$尺寸为$m\times n$，且$m>n$；$||x||=1$
则优化目标为$min||Ax||$

首先对优化目标做等价变化：
$$min||Ax||\Rightarrow min||x^T{A}^{T}Ax||\text{\hspace{1em}(5)}$$
对$A$做SVD分解（式2），考虑到$UV$的正交性，则有
$$\begin{array}{rl}{A}^{T}A=& V^T{\Sigma }^T{U}^{T}U\Sigma V^T\\ =& V{\Sigma }^T\Sigma V^T\\ =& V\left[\begin{array}{ccc}{\lambda }_0^2& & \\ & \ddots & \\ & & {\lambda }_n^2\end{array}\right]V^T\\ =& {\lambda }_0^2{v}_0{v}_0^T+...+{\lambda }_n^2{v}_n{v}_n^T\end{array}\text{\hspace{1em}(6)}$$
其中${\lambda }_0...{\lambda }_n$由大到小排列

令$x=\sum k_i{t}_i$，其中$t_i$为基向量，结合向量正交性质，则有
$$\begin{array}{rl}{x}^T{A}^{T}Ax=& {\lambda }_0^2{k}_0^2{t}_0^T{v}_0{v}_0^T{t}_0+...+{\lambda }_n^2{k}_n^2{t}_n^T{v}_n{v}_n^T{t}_n\\ =& {\lambda }_0^2{k}_0^2+...+{\lambda }_n^2{k}_n^2\end{array}\text{\hspace{1em}(7)}$$
由于$||x||=1$，则$\sum k_i^2=1$
式(7)取最小时有$k_n=1,k_i(i\ne n)=0$
此时$x=v_n$,即最小奇异值${\lambda }_n$对应的特征向量

DLT求单应性矩阵H

理论上8个自由度的单应性矩阵$H$只需要4组数据就可以唯一求解
但实际中考虑到数据的误差，一般希望通过多于4组的数据逼近真实值，此时$H$的求解属于一个超定问题

将式(1)展开，同时令$w=1$，有
$$\begin{gathered} h_{00}x+h_{01}y+h_{02}-h_{20}{x}^{\prime }x-h_{21}{x}^{\prime }y-h_{22}{x}^{\prime }=0 \\ {h}_{10}x+h_{11}y+h_{12}-h_{20}x{y}^{\prime }-h_{21}{y}^{\prime }y-h_{22}{y}^{\prime }=0\text{\hspace{1em}(8)} \end{gathered}$$
即
$$Ah=0\text{\hspace{1em}(9)}$$
其中
$$A=\left[\begin{array}{ccccccccc}x& y& 1& 0& 0& 0& -x^{\prime }x& -x^{\prime }y& -x^{\prime }\\ 0& 0& 0& x& y& 1& -x{y}^{\prime }& -y^{\prime }y& -y^{\prime }\end{array}\right]\text{\hspace{1em}(10)}$$
$$h={\left[\begin{array}{ccccc}{h}_{00}& h_{01}& ...& h_{21}& h_{22}\end{array}\right]}^T\text{\hspace{1em}(11)}$$

根据(9)式，得优化目标
$$h^{\prime }=min||Ah||\text{\hspace{1em}(12)}$$
通过上一节$Ax=0$模型求解的讨论，可知只需对$A$做SVD分解，取$v_n$即为所求。

Moving DLT求局部单应性矩阵

APAP算法在利用DLT求得全局单应性矩阵$H$的基础上，对每个小格子应用式(13)求得局部修正矩阵
$$h_{\ast }=min||W_{\ast }Ah||\text{\hspace{1em}(13)}$$
对比式(12)，局部单应性其实只是对$A$乘以了不同的权重

算法通过计算每个像素与格子左上角的坐标的欧式距离确定上述权重
$$w_i=max(e^{-||x_{\ast }-x_i|{|}^2/{\sigma }^2},\alpha )\text{\hspace{1em}(14)}$$
注：$W_{\ast }A$是对式(10)内的值做[0,1]加权，如果只是取0或者1比较好理解。但取中间的数，不知道能不能从数学上严格推导？

实现的小细节

1）找到两张图对应的特征点后，需要将这些点的坐标归一化到均值为0,标准差$\sqrt{2}$

  float cx = mat.col(0).mean();
  float cy = mat.col(1).mean();
  mat.array().col(0) -= cx;
  mat.array().col(1) -= cy;
  
  float sqrt_2 = sqrt(2);
  float meandist = (mat.array().col(0)*mat.array().col(0) + mat.array().col(1)*mat.array().col(1)).sqrt().mean();
  float scale = sqrt_2/meandist;
  mat.leftCols<2>().array() *= scale; 

  T << scale, 0, -scale*cx, 0, scale, -scale*cy, 0, 0, 1;

通过转换后的坐标计算$H^{\prime }$，然后通过$T$恢复到$H$
2) 图像拼接的时候首先根据B图warp后的4个顶点确定A图需要做的offX和offY
然后利用4邻域插值反向查找warp后的B图在原始B图中的位置

文献

【1】Zaragoza J, Chin T J, Brown M S, et al. As-projective-as-possible image stitching with moving DLT[C]//Proceedings of the IEEE conference on computer vision and pattern recognition. 2013: 2339-2346.