本文根据论文”Closed-form solution of absolute orientation using unit quaternions“对相似变换矩阵 ( s R t 0 1 ) \begin{pmatrix} sR & t \\0 & 1\end{pmatrix} (sR0t1)的求解过程中涉及到的原理和公式进行详细阐述。
概念
图形的相似变换是指由一个图形到另一个图形,在改变的过程中保持形状不变(大小方向和位置可变)的图形,类似于相似三角形的变换形式。相比于单纯的欧式变换相似变换多一个变换尺度信息 s s s,表示图像信息的缩放成度。在论文中提出了通过已知的图像间成功匹配的3对特征点求解图像之间存在的相似变换矩阵。在ORB_SLAM回环检测过程中需要对当前帧和回环备选帧之间的相似变换矩阵进行求解。
功能
用于表示俩个图像信息之间的相似变换,包括旋转矩阵 R R R、平移向量 t t t和图像之间的缩尺度信息 s s s。
步骤
-
已知参与运算的俩个图象上有 n n n对成功匹配的特征点坐标信息 { r l , i } a n d { r r , i } i = 1 ⋯ n . \{r_{l,i}\}\quad and \quad\{r_{r,i}\}\qquad i=1\cdots n. {rl,i}and{rr,i}i=1⋯n.这里的坐标信息表示归一化平面的齐次坐标信息 ( x z , y z , 1 ) \begin{pmatrix}\frac{x}{z},\frac{y}{z},1\end{pmatrix} (zx,zy,1)。
左右俩个图像之间的相似变换有 r r = s R ( r l ) + r 0 r_r=sR(r_l)+r_0 rr=sR(rl)+r0构建重投影误差信息 e i = r r , i − s R ( r l , i ) − r 0 e_i=r_{r,i}-sR(r_{l,i})-r_0 ei=rr,i−sR(rl,i)−r0我们的目标就是最小化该误差 ∑ i = 1 n ∣ ∣ e i ∣ ∣ 2 \sum^n_{i=1}||e_i||^2 i=1∑n∣∣ei∣∣2其中因为旋转矩阵 R R R是单位正交阵,因此 ∣ ∣ R ( r l ) ∣ ∣ 2 = ∣ ∣ r l ∣ ∣ 2 ||R(r_l)||^2=||r_l||^2 ∣∣R(rl)∣∣2=∣∣rl∣∣2接下来进行去质心处理,是已知的俩帧图像中的点坐标信息零均值 r l ˉ = 1 n ∑ i = 1 n r l , i r r ˉ = 1 n ∑ i = 1 n r r , i \bar{r_l}=\frac{1}{n}\sum ^n_{i=1}r_{l,i}\qquad\bar{r_r}=\frac{1}{n}\sum ^n_{i=1}r_{r,i} rlˉ=n1i=1∑nrl,irrˉ=n1i=1∑nrr,i r l , i ′ = r l , i − r l ˉ r r , i ′ = r r , i − r r ˉ r'_{l,i}=r_{l,i}-\bar{r_l}\qquad r'_{r,i}=r_{r,i}-\bar{r_r} rl,i′=rl,i−rlˉrr,i′=rr,i−rrˉ值得注意的是 ∑ i = 1 n r l , i ′ = 0 ∑ i = 1 n r r , i ′ = 0 \sum^n_{i=1}r'_{l,i}=0\qquad \sum^n_{i=1}r'_{r,i}=0 i=1∑nrl,i′=0i=1∑nrr,i′=0于是重投影误差可以写做如下形式 e i = r r , i ′ − s R ( r l , i ′ ) − r 0 ′ e_i=r'_{r,i}-sR(r'_{l,i})-r'_0 ei=rr,i′−sR(rl,i′)−r0′其中 r 0 ′ = r 0 − r r ˉ + s R ( r l ˉ ) r'_0=r_0-\bar{r_r}+sR(\bar{r_l}) r0′=r0−rrˉ+sR(rlˉ)因此误差的平方和形式可写成以下形式 ∑ i = 1 n ∣ ∣ r r , i ′ − s R ( r l , i ′ ) − r 0 ′ ∣ ∣ 2 \sum^n_{i=1}||r'_{r,i}-sR(r'_{l,i})-r'_0||^2 i=1∑n∣∣rr,i′−sR(rl,i′)−r0′∣∣2 即: ∑ i = 1 n ∣ ∣ r r , i ′ − s R ( r l , i ′ ) ∣ ∣ 2 − 2 r 0 ′ ⋅ ∑ i = 1 n [ r r , i ′ − s R ( r l , i ′ ) ] + n ∣ ∣ r 0 ′ ∣ ∣ 2 \sum^n_{i=1}||r'_{r,i}-sR(r'_{l,i})||^2-2r'_0\cdot \sum^n_{i=1}[r'_{r,i}-sR(r'_{l,i})]+n||r'_0||^2 i=1∑n∣∣rr,i′−sR(rl,i′)∣∣2−2r0′⋅i=1∑n[rr,i′−sR(rl,i′)]+n∣∣r0′∣∣2我们的目的是求解上述误差平方和公式的最小值,由于左侧第一项中和 r 0 ′ r'_0 r0′没有关系,为了使得上述公式值最小我们应该有 r 0 ′ = 0 o r r 0 = r r ˉ − s R ( r l ˉ ) r'_0=0\qquad or \qquad r_0=\bar{r_r}-sR(\bar{r_l}) r0′=0orr0=rrˉ−sR(rlˉ)至此我们已知需要求解的平移向量 r 0 r_0 r0可以通过尺度信息和旋转矩阵获取。 -
接下来我们分情况求解尺度信息 s s s,我们接下来的目的是在 r 0 ′ = 0 r'_0=0 r0′=0的情况下求解尺度信息和旋转矩阵是的下式值最小 ∑ i = 1 n ∣ ∣ r r , i ′ − s R ( r l , i ′ ) ∣ ∣ 2 \sum^n_{i=1}||r'_{r,i}-sR(r'_{l,i})||^2 i=1∑n∣∣rr,i′−sR(rl,i′)∣∣2由于 ∣ ∣ R ( r l , i ′ ) ∣ ∣ 2 = ∣ ∣ r l , i ′ ∣ ∣ 2 ||R(r'_{l,i})||^2=||r'_{l,i}||^2 ∣∣R(rl,i′)∣∣2=∣∣rl,i′∣∣2因此 ∑ i = 1 n ∣ ∣ r r , i ′ ∣ ∣ 2 − 2 s ∑ i = 1 n r r , i ′ ⋅ R ( r l , i ′ ) + s 2 ∑ i = 1 n ∣ ∣ r l , i ′ ∣ ∣ 2 \sum^n_{i=1}||r'_{r,i}||^2-2s\sum^n_{i=1}r'_{r,i}\cdot R(r'_{l,i})+s^2\sum^n_{i=1}||r'_{l,i}||^2 i=1∑n∣∣rr,i′∣∣2−2si=1∑nrr,i′⋅R(rl,i′)+s2i=1∑n∣∣rl,i′∣∣2上述公式可以写成 S r − 2 s D + s 2 S l S_r-2sD+s^2S_l Sr−2sD+s2Sl其中 { S r = ∑ i = 1 n ∣ ∣ r r , i ′ ∣ ∣ 2 D = ∑ i = 1 n r r , i ′ ⋅ R ( r l , i ′ ) S l = ∑ i = 1 n ∣ ∣ r l , i ′ ∣ ∣ 2 \begin{cases} S_r= \sum^n_{i=1}||r'_{r,i}||^2 \\ D= \sum^n_{i=1}r'_{r,i}\cdot R(r'_{l,i})\\ S_l = \sum^n_{i=1}||r'_{l,i}||^2 \end{cases} ⎩⎪⎨⎪⎧Sr=∑i=1n∣∣rr,i′∣∣2D=∑i=1nrr,i′⋅R(rl,i′)Sl=∑i=1n∣∣rl,i′∣∣2上述公式可以写成 ( s S l − D / S l ) 2 + ( S r S l − D 2 ) / S l (s\sqrt{S_l}-D/\sqrt{S_l})^2+(S_rS_l-D^2)/S_l (sSl−D/Sl)2+(SrSl−D2)/Sl为了让上述公式尽可能小,我们可以另其中第一项等于0,即另 s = D / S l s=D/S_l s=D/Sl
即 s = ∑ i = 1 n r r , i ′ ⋅ R ( r l , i ′ ) / ∑ i = 1 n ∣ ∣ r l , i ′ ∣ ∣ 2 s=\sum^n_{i=1}r'_{r,i}\cdot R(r'_{l,i}) /\sum^n_{i=1}||r'_{l,i}||^2 s=i=1∑nrr,i′⋅R(rl,i′)/i=1∑n∣∣rl,i′∣∣2至此我们已经获取了求解从左侧坐标系转换到右侧坐标系下是对应的尺度信息 s s s的求解公式。
同理我们获取从右侧坐标系转换到左侧坐标系下时对应的尺度信息 s s s的求解公式。已知 r r = s R ( r l ) + r 0 r_r=sR(r_l)+r_0 rr=sR(rl)+r0进行逆向转换有 r l = s ˉ R ˉ ( r l ) + r 0 ˉ r_l=\bar{s}\bar{R}(r_l)+\bar{r_0} rl=sˉRˉ(rl)+r0ˉ其中
{ s ˉ = 1 / s R ˉ = R − 1 r 0 ˉ = − 1 s R − 1 r 0 \begin{cases} \bar{s}=1/s \\ \bar{R}=R^{-1}\\ \bar{r_0}=-\frac1 sR^{-1}r_0 \end{cases} ⎩⎪⎨⎪⎧sˉ=1/sRˉ=R−1r0ˉ=−s1R−1r0
我们通过相同的求解方式获取了求解从右侧坐标系转换到左侧坐标系下是对应的尺度信息 s ˉ \bar s sˉ的求解公式 s ˉ = D ˉ / S r \bar{s}=\bar{D}/S_r sˉ=Dˉ/Sr。
即 s ˉ = ∑ i = 1 n r l , i ′ ⋅ R ˉ ( r r , i ′ ) / ∑ i = 1 n ∣ ∣ r r , i ′ ∣ ∣ 2 \bar s=\sum^n_{i=1}r'_{l,i}\cdot \bar R(r'_{r,i}) /\sum^n_{i=1}||r'_{r,i}||^2 sˉ=i=1∑nrl,i′⋅Rˉ(rr,i′)/i=1∑n∣∣rr,i′∣∣2我们发现获取的俩个尺度信息求解公式互为倒数,即结果不对称。如果已知左右俩帧图像中的某一侧精度更高我们可以选择使用其中的一个公式解算尺度信息。
如果左俩侧的图像帧精度相差无几,我们可以直接默认俩侧求解的尺度信息对称,于是可以通过 e i = 1 s r r , i ′ − s R ( r l , i ′ ) e_i=\frac{1}{\sqrt{s}}r'_{r,i}-\sqrt sR(r'_{l,i}) ei=s1rr,i′−sR(rl,i′)误差平方和为 1 s ∑ i = 1 n ∣ ∣ r r , i ′ ∣ ∣ 2 − 2 ∑ i = 1 n r r , i ′ ⋅ R ( r l , i ′ ) + s ∑ i = 1 n ∣ ∣ r l , i ′ ∣ ∣ 2 1 s S r − 2 D + s S l \frac{1}{s}\sum^n_{i=1}||r'_{r,i}||^2-2\sum^n_{i=1}r'_{r,i}\cdot R(r'_{l,i})+s\sum^n_{i=1}||r'_{l,i}||^2\\\frac{1}{s}S_r-2D+sS_l s1i=1∑n∣∣rr,i′∣∣2−2i=1∑nrr,i′⋅R(rl,i′)+si=1∑n∣∣rl,i′∣∣2s1Sr−2D+sSl调整 s s s是的上式最小化, s = S r / S l s=S_r/S_l s=Sr/Sl
即 s ˉ = ( ∑ i = 1 n ∣ ∣ r r , i ′ ∣ ∣ 2 / ∑ i = 1 n ∣ ∣ r l , i ′ ∣ ∣ 2 ) 1 2 \bar s=(\sum^n_{i=1}||r'_{r,i}||^2 /\sum^n_{i=1}||r'_{l,i}||^2)^\frac1 2 sˉ=(i=1∑n∣∣rr,i′∣∣2/i=1∑n∣∣rl,i′∣∣2)21值得注意的是,如果我们假设尺度因子的对称时可以直接通过俩帧图像中所给的坐标信息获取 s s s,而无需求解旋转矩阵信息 R R R。
至此,我们获取了求解尺度因子的方法公式。 -
接下来就是求解旋转矩阵信息,我们的目标是使得 S r − 2 s D + s 2 S l \;S_r-2sD+s^2S_l\; Sr−2sD+s2Sl的值尽可能的小,我们已知 S r , D , S l S_r,D,S_l Sr,D,Sl的具体表示形式,由于 ∣ ∣ r r , i ′ ⋅ R ( r l , i ′ ) ∣ ∣ 2 = ∣ ∣ r r , i ⋅ r l , i ∣ ∣ 2 ||r'_{r,i}\cdot R(r'_{l,i})||^2=||r_{r,i}\cdot r_{l,i}||^2 ∣∣rr,i′⋅R(rl,i′)∣∣2=∣∣rr,i⋅rl,i∣∣2所以有 D 2 < = S r ⋅ S l D^2<=S_r\cdot S_l D2<=Sr⋅Sl为了使得 S r − 2 s D + s 2 S l \;S_r-2sD+s^2S_l\; Sr−2sD+s2Sl尽可能小,我们应该是的 D D D 尽可能大,即使得 ∑ i = 1 n r r , i ′ ⋅ R ( r l , i ′ ) \sum^n_{i=1}r'_{r,i}\cdot R(r'_{l,i}) ∑i=1nrr,i′⋅R(rl,i′)尽可能大,通过这一条件来求解旋转矩阵 R R R。
这里开始引入四元素的知识,用四元数求解线性方程来获得最优的旋转信息,有关四元数的知识很常见,且论文中已经给出本文所用到的所有四元数知识,这里不再赘述。根据四元数的知识我们有 ∑ i = 1 n r r , i ′ ⋅ R ( r l , i ′ ) = ∑ i = 1 n ( q ˙ r l , i ′ q ˙ ∗ ) ⋅ r ˙ r , i ′ = ∑ i = 1 n ( q ˙ r l , i ′ ) ⋅ ( r r , i ′ q ˙ ) \sum^n_{i=1}r'_{r,i}\cdot R(r'_{l,i})=\sum^n_{i=1}(\dot{q}r'_{l,i}\dot{q}^*)\cdot \dot{r}'_{r,i}=\sum^n_{i=1}(\dot{q}r'_{l,i})\cdot (r'_{r,i}\dot{q}) i=1∑nrr,i′⋅R(rl,i′)=i=1∑n(q˙rl,i′q˙∗)⋅r˙r,i′=i=1∑n(q˙rl,i′)⋅(rr,i′q˙)已知 r l , i ′ = ( x l , i ′ , y l , i ′ , z l , i ′ ) r'_{l,i}=(x'_{l,i},y'_{l,i},z'_{l,i}) rl,i′=(xl,i′,yl,i′,zl,i′) , r r , i ′ = ( x r , i ′ , y r , i ′ , z r , i ′ ) r'_{r,i}=(x'_{r,i},y'_{r,i},z'_{r,i}) rr,i′=(xr,i′,yr,i′,zr,i′),于是有
q ˙ r ˙ l , i ′ = [ 0 − x l , i ′ − y l , i ′ − z l , i ′ x l , i ′ 0 z l , i ′ − y l , i ′ y l , i ′ − z l , i ′ 0 x l , i ′ z l , i ′ y l , i ′ − x l , i ′ 0 ] q ˙ = R ˉ l , i q ˙ \dot{q}\dot{r}'_{l,i}=\begin{bmatrix} 0 & -x'_{l,i} & -y'_{l,i} & -z'_{l,i}\\ x'_{l,i} & 0 & z'_{l,i} & -y'_{l,i} \\ y'_{l,i} & -z'_{l,i}&0& x'_{l,i} & \\ z'_{l,i}& y'_{l,i}&-x'_{l,i} &0\\ \end{bmatrix}\dot{q}=\bar{\R}_{l,i}\dot{q} q˙r˙l,i′=⎣⎢⎢⎡0xl,i′yl,i′zl,i′−xl,i′0−zl,i′yl,i′−yl,i′zl,i′0−xl,i′−zl,i′−yl,i′xl,i′0⎦⎥⎥⎤q˙=Rˉl,iq˙
q ˙ r ˙ r , i ′ = [ 0 − x r , i ′ − y r , i ′ − z r , i ′ x r , i ′ 0 z r , i ′ − y r , i ′ y r , i ′ − z r , i ′ 0 x r , i ′ z r , i ′ y r , i ′ − x r , i ′ 0 ] q ˙ = R ˉ r , i q ˙ \dot{q}\dot{r}'_{r,i}=\begin{bmatrix} 0 & -x'_{r,i} & -y'_{r,i} & -z'_{r,i}\\ x'_{r,i} & 0 & z'_{r,i} & -y'_{r,i} \\ y'_{r,i} & -z'_{r,i}&0& x'_{r,i} & \\ z'_{r,i}& y'_{r,i}&-x'_{r,i} &0\\ \end{bmatrix}\dot{q}=\bar{\R}_{r,i}\dot{q} q˙r˙r,i′=⎣⎢⎢⎡0xr,i′yr,i′zr,i′−xr,i′0−zr,i′yr,i′−yr,i′zr,i′0−xr,i′−zr,i′−yr,i′xr,i′0⎦⎥⎥⎤q˙=Rˉr,iq˙于是我们可以将上述最大化的目标函数写做 ∑ i = 1 n ( R ˉ l , i q ˙ ) ⋅ ( R ˉ r , i q ˙ ) = ∑ i = 1 n ˙ q T R ˉ l , i R ˉ r , i q ˙ = ˙ q T ( ∑ i = 1 n R ˉ l , i T R ˉ r , i ) q ˙ = q ˙ T ( ∑ i = 1 n N i ) q ˙ = q ˙ T N q ˙ \sum^n_{i=1}(\bar{\R}_{l,i}\dot{q})\cdot(\bar{\R}_{r,i}\dot{q})\\=\sum^n_{i=1}\dot{}q^T\bar{\R}_{l,i}\bar{\R}_{r,i}\dot{q}\\=\dot{}q^T(\sum^n_{i=1}\bar{\R}_{l,i}^T \;\bar{\R}_{r,i})\dot{q}\\=\dot{q}^T(\sum^n_{i=1}N_i)\dot{q}\\=\dot{q}^TN\dot{q} i=1∑n(Rˉl,iq˙)⋅(Rˉr,iq˙)=i=1∑n˙qTRˉl,iRˉr,iq˙=˙qT(i=1∑nRˉl,iTRˉr,i)q˙=q˙T(i=1∑nNi)q˙=q˙TNq˙矩阵 N N N是对称矩阵。
接下来我们引出矩阵 M M M,之后在此基础获得矩阵 N N N的具体表示形式 M = ∑ i = 1 n r l , i ′ r r , i ′ T M = [ S x x S x y S x z S y x S y y S y z S z x S z y S z z ] M=\sum^n_{i=1}r'_{l,i}\;r'^T_{r,i}\\M=\begin{bmatrix} S_{xx} & S_{xy} &S_{xz} \\ S_{yx} & S_{yy} &S_{yz} \\ S_{zx} & S_{zy} &S_{zz} \\ \end{bmatrix} M=i=1∑nrl,i′rr,i′TM=⎣⎡SxxSyxSzxSxySyySzySxzSyzSzz⎦⎤
其中 S x x ∑ i = 1 n x l , i ′ x r , i ′ S x y ∑ i = 1 n x l , i ′ y r , i ′ S_{xx}\sum^n_{i=1}x'_{l,i}x'_{r,i}\qquad \qquad S_{xy}\sum^n_{i=1}x'_{l,i}y'_{r,i} Sxxi=1∑nxl,i′xr,i′Sxyi=1∑nxl,i′yr,i′引出矩阵 N N N
N = [ ( S x x + S y y + S z z ) S y z − S z y S z x − S x z S x y − S y x S y z − S z y ( S x x − S y y − S z z ) S x y + S y x S z x + S x z S z x − S x z S x y + S y x ( − S x x + S y y − S z z ) S y z + S z y S x y − S y x S z x + S x z S y z + S z y ( − S x x − S y y + S z z ) ] N=\begin{bmatrix} (S_{xx}+S_{yy}+S_{zz}) & S_{yz}- S_{zy} & S_{zx}- S_{xz} & S_{xy}- S_{yx} \\ S_{yz}- S_{zy} & (S_{xx}-S_{yy}-S_{zz}) & S_{xy}+ S_{yx} & S_{zx}+S_{xz} \\ S_{zx}-S_{xz} & S_{xy}+ S_{yx} & (-S_{xx}+S_{yy}-S_{zz}) & S_{yz}+S_{zy} \\ S_{xy}-S_{yx} & S_{zx}+S_{xz} & S_{yz}+S_{zy} & (-S_{xx}-S_{yy}+S_{zz}) \\ \end{bmatrix} N=⎣⎢⎢⎡(Sxx+Syy+Szz)Syz−SzySzx−SxzSxy−SyxSyz−Szy(Sxx−Syy−Szz)Sxy+SyxSzx+SxzSzx−SxzSxy+Syx(−Sxx+Syy−Szz)Syz+SzySxy−SyxSzx+SxzSyz+Szy(−Sxx−Syy+Szz)⎦⎥⎥⎤
值得注意的是矩阵 N N N主对角线上的元素和为0。
为了使 q ˙ T N q ˙ \dot{q}^TN\dot{q} q˙TNq˙ 的值最大,我们需要获得矩阵 N N N的最大特征值和对应的特征向量解,对应最大的特征值为 λ m a x \lambda_{max} λmax对应的特征向量为 e m a x e_{max} emax,其中最大特征值对应的特征向量就是我们要求解的旋转矩阵对应的四元数表示形式,我们将问题转换为求解四个线性方程组的问题 [ N − λ m a x I ] e ˙ m a x = 0 [N-\lambda_{max}I]\dot e_{max}=0 [N−λmaxI]e˙max=0论文中介绍的是 F e r r a r i Ferrari Ferrari的求解方法,将问题转换为俩个二次方程进行求解。
至此,我们已经将论文中求解相似变换 S i m 3 Sim3 Sim3的主要部分讲解完毕,再次总结过程如下: -
对已知的左右俩侧图像点信息做去质心等处理;
-
构建重投影误差的平方和,将和平移向量有关的项归0时获得平移向量的求解公式;
-
将剩余的误差平方和项展开,将涉及到尺度因子 s s s的项归0可以获得左右俩种非对称的求解公式,如果我们假设左右尺度因子对称可以直通过已知的匹配点坐标进行求解;
-
最后将剩余的误差平方和部分最小化问题转换为将其中的部分项最大化,引出四元数,并对 N N N矩阵球最大特征值和对应的特征向量,该特征向量就是我们要求解的旋转矩阵的四元数形式;
-
将四元数转换为旋转矩阵,再求解尺度因子,在获得平移向量,于是图像之间的相似变换矩阵 ( s R t 0 1 ) \begin{pmatrix} sR & t \\ 0 & 1 \end{pmatrix} (sR0t1)求解成功。