计算机视觉系列教程1-2：单应性矩阵估计

Mr_Winter
1,149 阅读5分钟

「这是我参与2022首次更文挑战的第11天，活动详情查看：2022首次更文挑战」

1 导论

计算机视觉系列教程1-1：透视空间与透视变换中提到，透视空间所有变换都是投影变换的特例，本节进一步研究投影变换矩阵(单应性矩阵)的估计。

透视变换的核心是单应性矩阵HH或单应性向量hh。

H=[h11h12h13h21h22h23h31h32h33]⇔h=[h11h12h13h21h22h23h31h32h33]TH=\left[ \begin{matrix} h_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\\\end{matrix} \right] \Leftrightarrow h=\left[ \begin{matrix} h_{11}& h_{12}& h_{13}& h_{21}& h_{22}& h_{23}& h_{31}& h_{32}& h_{33}\\\end{matrix} \right] ^T

设psrc=[xy1]Tp_{src}=\left[ \begin{matrix} x& y& 1\\\end{matrix} \right] ^T与pdst=[x′y′1]Tp_{dst}=\left[ \begin{matrix} x'& y'& 1\\\end{matrix} \right] ^T是二维透视空间P2\mathbb{P}^2中，一次透视变换前后的对应点，因此其满足

pdst=Hpsrc⟺[x′y′1]=[h11h12h13h21h22h23h31h32h33][xy1]p_{dst}=Hp_{src}\Longleftrightarrow \left[ \begin{array}{c} x'\\ y'\\ 1\\\end{array} \right] =\left[ \begin{matrix} h_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\\\end{matrix} \right] \left[ \begin{array}{c} x\\ y\\ 1\\\end{array} \right]

若将单应性矩阵进行尺度缩放后作用于psrcp_{src}，则

kHpsrc=k[h11h12h13h21h22h23h31h32h33][xy1]=kpdstkHp_{src}=k\left[ \begin{matrix} h_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\\\end{matrix} \right] \left[ \begin{array}{c} x\\ y\\ 1\\\end{array} \right] =kp_{dst}

而透视空间中，kpdstkp_{dst}与pdstp_{dst}实际上对应同一点，因此kHkH与HH相当于同一次透视变换，故 单应性矩阵HH仅有8个自由度，通常通过设置h33=1h_{33}=1或∥h∥22=1\lVert h \rVert _{2}^{2}=1来约束冗余的参数。

下面详细阐述单应性矩阵估计方法。

2 基本直接线性变换(Basic DLT)

将上式改写为齐次形式

[000−x−y−1y′xy′yy′xy1000−x′x−x′y−x′−y′x−y′y−y′x′xx′yx′000][h11h12h13h21h22h23h31h32h33]=[000]\left[ \begin{matrix} 0& 0& 0& -x& -y& -1& y'x& y'y& y'\\ x& y& 1& 0& 0& 0& -x'x& -x'y& -x'\\ -y'x& -y'y& -y'& x'x& x'y& x'& 0& 0& 0\\\end{matrix} \right] \left[ \begin{array}{c} h_{11}\\ h_{12}\\ h_{13}\\ h_{21}\\ h_{22}\\ h_{23}\\ h_{31}\\ h_{32}\\ h_{33}\\\end{array} \right] =\left[ \begin{array}{c} 0\\ 0\\ 0\\\end{array} \right]

其中系数矩阵的秩为2，因此一对变换点仅能确定2个自由度。因此需要无三点共线的四对变换点才能确定单应性矩阵HH。

[000−x1−y1−1y1′x1y1′y1y1′x1y11000−x1′x1−x1′y1−x1′000−x2−y2−1y2′x2y2′y2y2′x2y21000−x2′x2−x2′y2−x2′⋮][h11h12h13h21h22h23h31h32h33]=[000]⇔Ah=0\left[ \begin{matrix} 0& 0& 0& -x_1& -y_1& -1& y_{1}^{'}x_1& y_{1}^{'}y_1& y_{1}^{'}\\ x_1& y_1& 1& 0& 0& 0& -x_{1}^{'}x_1& -x_{1}^{'}y_1& -x_{1}^{'}\\ 0& 0& 0& -x_2& -y_2& -1& y_{2}^{'}x_2& y_{2}^{'}y_2& y_{2}^{'}\\ x_2& y_2& 1& 0& 0& 0& -x_{2}^{'}x_2& -x_{2}^{'}y_2& -x_{2}^{'}\\ & & & & \vdots& & & & \\\end{matrix} \right] \left[ \begin{array}{c} h_{11}\\ h_{12}\\ h_{13}\\ h_{21}\\ h_{22}\\ h_{23}\\ h_{31}\\ h_{32}\\ h_{33}\\\end{array} \right] =\left[ \begin{array}{c} 0\\ 0\\ 0\\\end{array} \right] \Leftrightarrow {Ah=0}

对于形如Ax=bAx=b的非齐次线性方程可以通过伪逆形式计算xx；对于形如Ax=0Ax=0的齐次线性方程的解则对应矩阵AA右奇异向量vpv_p，其中vpv_p对应的奇异值σp≈0\sigma _p\approx 0或不对应奇异值。

3 归一化直接线性变换(Normalized DLT)

基本DLT估计算法的缺陷是：

  • 单应性估计不具有相似不变性。假设在第一次估计下有xdst=Hxsrcx_{dst}=Hx_{src}。现对两张图像分别进行相似变换并重新进行单应性估计，得到(Tdstxdst)=H′(Tsrcxsrc)\left( T_{dst}x_{dst} \right) =H'\left( T_{src}x_{src} \right)，改写为xdst=(Tdst−1H′Tsrc)xsrcx_{dst}=\left( T_{dst}^{-1}H'T_{src} \right) x_{src}，大部分情况下H≠Tdst−1H′TsrcH\ne T_{dst}^{-1}H'T_{src}，这表明基本DLT算法无法抵抗相似变换的干扰。
  • 估计的单应性矩阵容易产生病态条件，鲁棒性差。由于默认透视空间的尺度变换因子w=1w=1，所以齐次坐标下很可能产生分量幅度差异大的情况，例如某特征点 X=[1001011]TX=\left[ \begin{matrix} 100& 101& 1\\\end{matrix} \right] ^T。在这种情况下估计出的单应性矩阵，各个元素数值数量级可能会相差10410^4以上，导致病态条件——特征点的轻微变化都会造成单应性矩阵的剧变。

基于以上两种缺陷，需要将基本DLT算法进行优化,优化的核心就是特征点坐标的归一化。设原图像特征点集合为 ，目标图像特征点集合为 ，则具体的算法为：

  1. 将特征点集合XsrcX_{src}、XdstX_{dst}归一化。使用相似变换矩阵TsrcT_{src}、TdstT_{dst}将特征点集合中心移至原点，且与原点平均距离为2\sqrt{2}。由于默认尺度因子为w=1w=1，所以归一化到2\sqrt{2}可以保持齐次坐标的三个分量有相同的幅度，例如X=[1001001]T⇒Xnormal=[111]TX=\left[ \begin{matrix} 100& 100& 1\\\end{matrix} \right] ^T\Rightarrow X^{normal}=\left[ \begin{matrix} 1& 1& 1\\\end{matrix} \right] ^T。这里给出一种相似变换矩阵的计算方式，设
Xnormal=TX⇔[x~iy~i1]=[stxsty1][xiyi1]X^{normal}=TX\Leftrightarrow \left[ \begin{array}{c} \tilde{x}_i\\ \tilde{y}_i\\ 1\\\end{array} \right] =\left[ \begin{matrix} s& & t_x\\ & s& t_y\\ & & 1\\\end{matrix} \right] \left[ \begin{array}{c} x_i\\ y_i\\ 1\\\end{array} \right]

令

{1N∑i=1Nx~i=1N∑i=1N(sxi+tx)=01N∑i=1Ny~i=1N∑i=1N(syi+ty)=01N∑i=1Nx~i2+y~i2=2\begin{cases} \frac{1}{N}\sum_{i=1}^N{\tilde{x}_i}=\frac{1}{N}\sum_{i=1}^N{\left( sx_i+t_x \right)}=0\\ \frac{1}{N}\sum_{i=1}^N{\tilde{y}_i}=\frac{1}{N}\sum_{i=1}^N{\left( sy_i+t_y \right)}=0\\ \frac{1}{N}\sum_{i=1}^N{\sqrt{\tilde{x}_{i}^{2}+\tilde{y}_{i}^{2}}}=\sqrt{2}\\\end{cases}

解得

{tx=−s1N∑i=1Nxi=−sxˉty=−s1N∑i=1Nyi=−syˉs=21N∑i=1Nx~i2+y~i2=21N∑i=1N(xi−xˉ)2+(yi−yˉ)2\begin{cases} t_x=-s\frac{1}{N}\sum_{i=1}^N{x_i}=-s\bar{x}\\ t_y=-s\frac{1}{N}\sum_{i=1}^N{y_i=-s\bar{y}}\\ s=\frac{\sqrt{2}}{\frac{1}{N}\sum_{i=1}^N{\sqrt{\tilde{x}_{i}^{2}+\tilde{y}_{i}^{2}}}}=\frac{\sqrt{2}}{\frac{1}{N}\sum_{i=1}^N{\sqrt{\left( x_i-\bar{x} \right) ^2+\left( y_i-\bar{y} \right) ^2}}}\\\end{cases}
  1. 运用基本DLT算法由 与 估计单应性矩阵
  2. 解归一化，映射回实际图像

4 鲁棒单应性估计(Robust Homography Estimation)

结合基本DLT算法、归一化DLT算法及RANSAC算法进行单应性矩阵估计，具体流程如下：

  1. 设置迭代次数K=∞K=\infty，内点集Sin=⊘S_{in}=\oslash，模型参数H=H0H=H_0；
  2. 随机从样本数据集SS中选取4对特征点，并通过基本DLT算法确定测试模型HtestH_{test}；
  3. 用HtestH_{test}遍历样本数据集SS，估计误差ε\varepsilon在距离阈值tt内的点加入内点集SinS_{in}。其中阈值t=5.99σt=\sqrt{5.99}\sigma，σ\sigma为估计不确定度，估计误差ε\varepsilon主要有两种度量方式：① 代数误差εi=∥Aihtest∥\varepsilon _i=\left\| A_ih_{test} \right\|，其中
Ai=[000−x−y−1y′xy′yy′−y′x−y′y−y′x′xx′yx′000]A_i=\left[ \begin{matrix} 0& 0& 0& -x& -y& -1& y'x& y'y& y'\\ -y'x& -y'y& -y'& x'x& x'y& x'& 0& 0& 0\\\end{matrix} \right]

② 几何误差(二次投影误差)εi=∥HXsrc,i,Xdst,i∥22+∥Xsrc,i,H−1Xdst,i∥22\varepsilon _i=\left\| HX_{src,i}, X_{dst, i} \right\| _{2}^{2}+\left\| X_{src,i}, H^{-1}X_{dst, i} \right\| _{2}^{2}，可以视作交叉检验。

  1. 若SinS_{in}的大小小于阈值TT，则放弃该模型，重复(2)；若SinS_{in}的大小大于阈值tt，则通过归一化DLT算法或Levenberg Marquardt等迭代优化算法，利用SinS_{in}中的所有点重新估计模型Htest∗H_{test}^{*}；
  2. 计算当前模型Htest∗H_{test}^{*}下的内点率ω=∣Sin∣∣S∣\omega =\frac{|S_{in}|}{|S|}，根据K=ln⁡(1−p)ln⁡(1−ωn)K=\frac{\ln \left( 1-p \right)}{\ln \left( 1-\omega ^n \right)}更新迭代次数；
  3. 至此完成一次迭代，若Htest∗H_{test}^{*}下内点率为最大，则更新H=Htest∗H=H_{test}^{*}，重复(2) ~ (5)直至迭代次数满足要求。

计算机视觉基础教程大纲

章号                                    内容

  0                              色彩空间与数字成像

  1                              计算机几何基础

  2                              图像增强、滤波、金字塔

  3                              图像特征提取

  4                              图像特征描述

  5                              图像特征匹配

  6                              立体视觉

  7                              项目实战

更多内容欢迎来我的AI频道“AI技术社”

avatar
文章
阅读
粉丝