若场景中的特征点都落在同一平面上(比如墙、地面等),则可以通过单应性进行运动估计。这种情况在无人机携带的俯视相机和扫地机器人携带的顶视相机中比较常见。
单应矩阵 H\boldsymbol{H}H 通常描述处于共同平面上的一些点在两张图片之间的变换关系。
设图像 I1I_1I1 和 I2I_2I2 有一对匹配好的特征点 p1\boldsymbol{p}_1p1 和 p2\boldsymbol{p}_2p2(像素坐标系)。特征点落在平面 P\boldsymbol{P}P 上,设这个平面满足方程:
nTP+d=0(1) \boldsymbol{n}^{T} \boldsymbol{P}+d=0 \tag{1} nTP+d=0(1)
其中,n\boldsymbol{n}n 为平面 P\boldsymbol{P}P 的法向量,ddd 为坐标原点到平面的距离。
整理之后,有:
−nTPd=1(2) -\frac{\boldsymbol{n}^{T} \boldsymbol{P}}{d}=1 \tag{2} −dnTP=1(2)
由针孔相机模型有:
p2≃K(RP+t)≃K(RP+t⋅(−nTPd))≃K(R−tnTd)P≃K(R−tnTd)K−1p1(3) \begin{aligned}\boldsymbol{p}_{2} &\simeq\boldsymbol{K}(\boldsymbol{R} \boldsymbol{P}+\boldsymbol{t}) \\&\simeq\boldsymbol{K}\left(\boldsymbol{R} \boldsymbol{P}+\boldsymbol{t} \cdot\left(-\frac{\boldsymbol{n}^{T} \boldsymbol{P}}{d}\right)\right) \\&\simeq\boldsymbol{K}\left(\boldsymbol{R}-\frac{\boldsymbol{t n}^{T}}{d}\right) \boldsymbol{P} \\&\simeq\boldsymbol{K}\left(\boldsymbol{R}-\frac{\boldsymbol{t n}^{T}}{d}\right) \boldsymbol{K}^{-1} \boldsymbol{p}_{1}\end{aligned} \tag{3} p2≃K(RP+t)≃K(RP+t⋅(−dnTP))≃K(R−dtnT)P≃K(R−dtnT)K−1p1(3)
于是,我们得到了一个直接描述图像坐标 p1\boldsymbol{p}_1p1 和 p2\boldsymbol{p}_2p2 之间的变换,把中间的部分记做 H\boldsymbol{H}H,也就是单应矩阵。于是有:
p2≃Hp1(4) \boldsymbol{p}_{2}\simeq\boldsymbol{H} \boldsymbol{p}_{1} \tag{4} p2≃Hp1(4)
它的定义与旋转、平移以及平面的参数有关。与基础矩阵 F\boldsymbol{F}F 类似,单应矩阵 H\boldsymbol{H}H 也是一个 3 × 3 的矩阵,求解时的思路和 F\boldsymbol{F}F 类似,同样地可以先根据匹配点计算 H\boldsymbol{H}H ,然后将它分解以计算旋转和平移。把上式展开,得:
(u2v21)≃(h1h2h3h4h5h6h7h8h9)(u1v11)(5) \left(\begin{array}{c}u_{2} \\v_{2} \\1\end{array}\right)\simeq\left(\begin{array}{lll}h_{1} & h_{2} & h_{3} \\h_{4} & h_{5} & h_{6} \\h_{7} & h_{8} & h_{9}\end{array}\right)\left(\begin{array}{c}u_{1} \\v_{1} \\1\end{array}\right)\tag{5} ⎝⎛u2v21⎠⎞≃⎝⎛h1h4h7h2h5h8h3h6h9⎠⎞⎝⎛u1v11⎠⎞(5)
注意,这里的等号依然是 ≃\simeq≃ 而不是普通的等号,因此 H\boldsymbol{H}H 也可以乘以任何的常数。实际处理过程中可以令 h9=1h_9 = 1h9=1 (当 h9≠0h_9 \neq 0h9=0 时)。 然后根据第3行去掉这个非零因子,于是有:
u2=h1u1+h2v1+h3h7u1+h8v1+h9 u_{2}=\frac{h_{1} u_{1}+h_{2} v_{1}+h_{3}}{h_{7} u_{1}+h_{8} v_{1}+h_{9}} u2=h7u1+h8v1+h9h1u1+h2v1+h3
v2=h4u1+h5v1+h6h7u1+h8v1+h9 v_{2}=\frac{h_{4} u_{1}+h_{5} v_{1}+h_{6}}{h_{7} u_{1}+h_{8} v_{1}+h_{9}} v2=h7u1+h8v1+h9h4u1+h5v1+h6
整理和化简后得:
h1u1+h2v1+h3−h7u1u2−h8v1u2=u2h4u1+h5v1+h6−h7u1v2−h8v1v2=v2(6) \begin{aligned}&h_{1} u_{1}+h_{2} v_{1}+h_{3}-h_{7} u_{1} u_{2}-h_{8} v_{1} u_{2}=u_{2} \\&h_{4} u_{1}+h_{5} v_{1}+h_{6}-h_{7} u_{1} v_{2}-h_{8} v_{1} v_{2}=v_{2}\end{aligned}\tag{6} h1u1+h2v1+h3−h7u1u2−h8v1u2=u2h4u1+h5v1+h6−h7u1v2−h8v1v2=v2(6)
这样一组匹配点对就可以构造出两项约束(事实上有三个约束,但是因为线性相关,只 取前两个),于是自由度为 8 的单应矩阵可以通过 4 对匹配特征点算出(注意:这些特征 点不能有三点共线的情况),即求解以下的线性方程组(当 h9=0h_9 = 0h9=0 时,右侧为零):
(u11v111000−u11u21−v11u21000u11v111−u11v21−v11v21u12v121000−u12u22−v12u22000u12v121−u12v22−v12v22u13v131000−u13u23−v13u23000u13v131−u13v23−v13v23u14v141000−u14u24−v14u24000u14v141−u14v24−v14v24)(h1h2h3h4h5h6h7h8)=(u21v21u22v22u23v23u24v24)(7) \left(\begin{array}{cccccccc}u_{1}^{1} & v_{1}^{1} & 1 & 0 & 0 & 0 & -u_{1}^{1} u_{2}^{1} & -v_{1}^{1} u_{2}^{1} \\0 & 0 & 0 & u_{1}^{1} & v_{1}^{1} & 1 & -u_{1}^{1} v_{2}^{1} & -v_{1}^{1} v_{2}^{1} \\u_{1}^{2} & v_{1}^{2} & 1 & 0 & 0 & 0 & -u_{1}^{2} u_{2}^{2} & -v_{1}^{2} u_{2}^{2} \\0 & 0 & 0 & u_{1}^{2} & v_{1}^{2} & 1 & -u_{1}^{2} v_{2}^{2} & -v_{1}^{2} v_{2}^{2} \\u_{1}^{3} & v_{1}^{3} & 1 & 0 & 0 & 0 & -u_{1}^{3} u_{2}^{3} & -v_{1}^{3} u_{2}^{3} \\0 & 0 & 0 & u_{1}^{3} & v_{1}^{3} & 1 & -u_{1}^{3} v_{2}^{3} & -v_{1}^{3} v_{2}^{3} \\u_{1}^{4} & v_{1}^{4} & 1 & 0 & 0 & 0 & -u_{1}^{4} u_{2}^{4} & -v_{1}^{4} u_{2}^{4} \\0 & 0 & 0 & u_{1}^{4} & v_{1}^{4} & 1 & -u_{1}^{4} v_{2}^{4} & -v_{1}^{4} v_{2}^{4}\end{array}\right)\left(\begin{array}{l}h_{1} \\h_{2} \\h_{3} \\h_{4} \\h_{5} \\h_{6} \\h_{7} \\h_{8}\end{array}\right)=\left(\begin{array}{c}u_{2}^{1} \\v_{2}^{1} \\u_{2}^{2} \\v_{2}^{2} \\u_{2}^{3} \\v_{2}^{3} \\u_{2}^{4} \\v_{2}^{4}\end{array}\right) \tag{7} ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛u110u120u130u140v110v120v130v140101010100u110u120u130u140v110v120v130v1401010101−u11u21−u11v21−u12u22−u12v22−u13u23−u13v23−u14u24−u14v24−v11u21−v11v21−v12u22−v12v22−v13u23−v13v23−v14u24−v14v24⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛h1h2h3h4h5h6h7h8⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛u21v21u22v22u23v23u24v24⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞(7)
这种做法把 H 矩阵看成了向量,通过解该向量的线性方程来恢复 H\boldsymbol{H}H ,又称直接线性变换法(Direct Linear Transform)。与本质矩阵相似,求出单应矩阵以后需要对其进行分 解,才可以得到相应的旋转矩阵 R\boldsymbol{R}R 是和平移向量 t\boldsymbol{t}t。
分解的方法包括数值法与解析法。与本质矩阵的分解类似,单应矩阵的分解同样会返回四组旋转矩阵与平移向量,并且同时可以计算出它们分别对应的场景点所在平面的法向量。如果已知成像的地图点的深 度全为正值(即在相机前方),则又可以排除两组解。最后仅剩两组解,这时需要通过更多的先验信息进行判断。通常我们可以通过假设已知场景平面的法向量来解决,如场景平面与相机平面平行,那么法向量 n\boldsymbol{n}n 的理论值为 1T\boldsymbol{1}^T1T。
单应矩阵的意义
单应性在 SLAM 中具重要意义。当特征点共面,或者相机发生纯旋转的时候,基础矩阵的自由度下降,这就出现了所谓的退化(degenerate)。
现实中的数据总包含一些噪声, 这时候如果我们继续使用八点法求解基础矩阵,基础矩阵多余出来的自由度将会主要由噪 声决定。
为了能够避免退化现象造成的影响,通常我们会同时估计基础矩阵 F\boldsymbol{F}F 和单应矩阵 H\boldsymbol{H}H,选择重投影误差比较小的那个作为最终的运动估计矩阵。