大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。
Jetbrains全系列IDE使用 1年只要46元 售后保障 童叟无欺
图像拼接首先通过一些手段（标定、SIFT等特征点、其他传感器等）获取图像间的对应关系（2D-2D）。
这个对应关系可以用单应性矩阵 H H H（也称透视变换）描述
[ x ′ y ′ w ′ ] = [ h 00 h 01 h 02 h 10 h 11 h 12 h 20 h 21 h 22 ] [ x y w ] = H [ x y w ] (1) \begin{aligned}\left[\begin{matrix}x’\\y’\\w’\end{matrix}\right]=&\left[\begin{matrix}h_{00}&h_{01}&h_{02}\\h_{10}&h_{11}&h_{12}\\h_{20}&h_{21}&h_{22}\end{matrix}\right]\left[\begin{matrix}x\\y\\w\end{matrix}\right]\\=&H\left[\begin{matrix}x\\y\\w\end{matrix}\right]\end{aligned} \tag{1} ⎣⎡x′y′w′⎦⎤==⎣⎡h00h10h20h01h11h21h02h12h22⎦⎤⎣⎡xyw⎦⎤H⎣⎡xyw⎦⎤(1)
但上述思路的前提是：待拼接的目标处于或者近似处于同一平面。否则就会产生遮挡和角度差异，从而造成拼接图的“虚影”
图1. 由于相机角度差异造成的虚影现象（红圈）
针对上述问题，Autostitch（OpenCV实现的方法）和Photosynth等方法走的是一个全局H + 虚影后处理的路线；
APAP【1】算法则使用多个局部H的方式尝试在拼接时就尽可能精确。
APAP算法思路
APAP算法官方主页提供了matlab的源码
Github上有第三方实现的C++代码
算法以A图为基准，首先计算B图到A图的全局变换 H H H（DLT方法）,
然后将B图等间距划分小网格，逐个计算B图每个小网格与A图的局部变换 H ∗ H_* H∗（Moving DLT方法）,
最后将B图每个网格内的点按对应的 H ∗ H_* H∗做局部变换，与A图融合。
图2. 网格局部变换（摘自APAP论文）
几个问题
SVD在最小二乘中的应用
– 模型 A x = b Ax=b Ax=b的求解
令 A A A尺寸为 m × n m\times n m×n，且 m > n m>n m>n（超定问题）
则优化目标为 m i n ∣ ∣ A x − b ∣ ∣ min||Ax-b|| min∣∣Ax−b∣∣
首先对 A A A做SVD分解，有
A = U Σ V T (2) A=U\Sigma V^T \tag{2} A=UΣVT(2)
其中 U U U和 V V V分别称为 A A A的左、右奇异向量，满足正交性（即 U T = U − 1 U^T=U^{-1} UT=U−1）
Σ \Sigma Σ为对角阵，对角元素是 A A A的奇异值
将式(2)代入目标函数，有
m i n ∣ ∣ A x − b ∣ ∣ → m i n ∣ ∣ U Σ V T x − b ∣ ∣ → m i n ∣ ∣ Σ V T x − U T b ∣ ∣ ( 正 交 性 ) → m i n ∣ ∣ Σ y − b ′ ∣ ∣ ( 令 y = V T x , b ′ = U T b ) (3) \begin{aligned}&min||Ax-b||\\&\rightarrow min||U\Sigma V^Tx-b||\\&\rightarrow min||\Sigma V^Tx-U^Tb|| (正交性) \\&\rightarrow min||\Sigma y-b’|| (令y=V^Tx,b’=U^Tb)\end{aligned} \tag{3} min∣∣Ax−b∣∣→min∣∣UΣVTx−b∣∣→min∣∣ΣVTx−UTb∣∣(正交性)→min∣∣Σy−b′∣∣(令y=VTx,b′=UTb)(3)
注意到 Σ \Sigma Σ是一个对角阵，令 λ i \lambda_i λi为其对角线元素，则有
y i = b i ′ λ i x = V y (4) y_i=\frac{b’_i}{\lambda_i}\\x=Vy \tag{4} yi=λibi′x=Vy(4)
– 模型 A x = 0 Ax=0 Ax=0的求解
令 A A A尺寸为 m × n m\times n m×n，且 m > n m>n m>n； ∣ ∣ x ∣ ∣ = 1 ||x||=1 ∣∣x∣∣=1
则优化目标为 m i n ∣ ∣ A x ∣ ∣ min||Ax|| min∣∣Ax∣∣
首先对优化目标做等价变化：
m i n ∣ ∣ A x ∣ ∣ ⇒ m i n ∣ ∣ x T A T A x ∣ ∣ (5) min||Ax||\Rightarrow min||x^TA^TAx|| \tag{5} min∣∣Ax∣∣⇒min∣∣xTATAx∣∣(5)
对 A A A做SVD分解（式2），考虑到 U V UV UV的正交性，则有
A T A = V T Σ T U T U Σ V T = V Σ T Σ V T = V [ λ 0 2 ⋱ λ n 2 ] V T = λ 0 2 v 0 v 0 T + . . . + λ n 2 v n v n T (6) \begin{aligned}A^TA=&V^T\Sigma^TU^TU\Sigma V^T\\=&V\Sigma^T\Sigma V^T\\=&V\left[\begin{matrix}\lambda_0^2\\&\ddots\\&&\lambda_n^2\end{matrix}\right]V^T\\=&\lambda_0^2v_0v_0^T+…+\lambda_n^2v_nv_n^T\end{aligned} \tag{6} ATA====VTΣTUTUΣVTVΣTΣVTV⎣⎡λ02⋱λn2⎦⎤VTλ02v0v0T+...+λn2vnvnT(6)
其中 λ 0 . . . λ n \lambda_0…\lambda_n λ0...λn由大到小排列
令 x = ∑ k i t i x=\sum{k_it_i} x=∑kiti，其中 t i t_i ti为基向量，结合向量正交性质，则有
x T A T A x = λ 0 2 k 0 2 t 0 T v 0 v 0 T t 0 + . . . + λ n 2 k n 2 t n T v n v n T t n = λ 0 2 k 0 2 + . . . + λ n 2 k n 2 (7) \begin{aligned}x^TA^TAx=&\lambda_0^2k_0^2t_0^Tv_0v_0^Tt_0+…+\lambda_n^2k_n^2t_n^Tv_nv_n^Tt_n\\=&\lambda_0^2k_0^2+…+\lambda_n^2k_n^2\end{aligned} \tag{7} xTATAx==λ02k02t0Tv0v0Tt0+...+λn2kn2tnTvnvnTtnλ02k02+...+λn2kn2(7)
由于 ∣ ∣ x ∣ ∣ = 1 ||x||=1 ∣∣x∣∣=1，则 ∑ k i 2 = 1 \sum{k_i^2}=1 ∑ki2=1
式(7)取最小时有 k n = 1 , k i ( i ≠ n ) = 0 k_n=1,k_i(i\neq n)=0 kn=1,ki(i=n)=0
此时 x = v n x=v_n x=vn,即最小奇异值 λ n \lambda_n λn对应的特征向量
DLT求单应性矩阵H
理论上8个自由度的单应性矩阵 H H H只需要4组数据就可以唯一求解
但实际中考虑到数据的误差，一般希望通过多于4组的数据逼近真实值，此时 H H H的求解属于一个超定问题
将式(1)展开，同时令 w = 1 w=1 w=1，有
h 00 x + h 01 y + h 02 − h 20 x ′ x − h 21 x ′ y − h 22 x ′ = 0 h 10 x + h 11 y + h 12 − h 20 x y ′ − h 21 y ′ y − h 22 y ′ = 0 (8) h_{00}x+h_{01}y+h_{02}-h_{20}x’x-h_{21}x’y-h_{22}x’=0\\h_{10}x+h_{11}y+h_{12}-h_{20}xy’-h_{21}y’y-h_{22}y’=0 \tag{8} h00x+h01y+h02−h20x′x−h21x′y−h22x′=0h10x+h11y+h12−h20xy′−h21y′y−h22y′=0(8)
即
A h = 0 (9) Ah=0 \tag{9} Ah=0(9)
其中
A = [ x y 1 0 0 0 − x ′ x − x ′ y − x ′ 0 0 0 x y 1 − x y ′ − y ′ y − y ′ ] (10) A=\left[\begin{matrix}x&y&1&0&0&0&-x’x&-x’y&-x’\\0&0&0&x&y&1&-xy’&-y’y&-y’\end{matrix}\right] \tag{10} A=[x0y0100x0y01−x′x−xy′−x′y−y′y−x′−y′](10)
h = [ h 00 h 01 . . . h 21 h 22 ] T (11) h=\left[\begin{matrix}h_{00}&h_{01}&…&h_{21}&h_{22}\end{matrix}\right]^T \tag{11} h=[h00h01...h21h22]T(11)
根据(9)式，得优化目标
h ′ = m i n ∣ ∣ A h ∣ ∣ (12) h’=min||Ah|| \tag{12} h′=min∣∣Ah∣∣(12)
通过上一节 A x = 0 Ax=0 Ax=0模型求解的讨论，可知只需对 A A A做SVD分解，取 v n v_n vn即为所求。
Moving DLT求局部单应性矩阵
APAP算法在利用DLT求得全局单应性矩阵 H H H的基础上，对每个小格子应用式(13)求得局部修正矩阵
h ∗ = m i n ∣ ∣ W ∗ A h ∣ ∣ (13) h_*=min||W_*Ah|| \tag{13} h∗=min∣∣W∗Ah∣∣(13)
对比式(12)，局部单应性其实只是对 A A A乘以了不同的权重
算法通过计算每个像素与格子左上角的坐标的欧式距离确定上述权重
w i = m a x ( e − ∣ ∣ x ∗ − x i ∣ ∣ 2 / σ 2 , α ) (14) w_i=max(e^{-||x_*-x_i||^2/\sigma^2},\alpha) \tag{14} wi=max(e−∣∣x∗−xi∣∣2/σ2,α)(14)
注： W ∗ A W_*A W∗A是对式(10)内的值做[0,1]加权，如果只是取0或者1比较好理解。但取中间的数，不知道能不能从数学上严格推导？
实现的小细节
1）找到两张图对应的特征点后，需要将这些点的坐标归一化到均值为0,标准差 2 \sqrt{2} 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 ′ H’ H′，然后通过 T T T恢复到 H H 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.
发布者：全栈程序员-用户IM，转载请注明出处：https://javaforall.cn/190681.html原文链接：https://javaforall.cn
【正版授权，激活自己账号】： Jetbrains全家桶Ide使用，1年售后保障，每天仅需1毛
【官方授权 正版激活】： 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...
