多视图几何(3)_估计—2D射影变换

估计—2D射影变换

直接线性变换算法DLT

首先讨论已知4组2D到2D的对应点$x_i⟺x_i^{\prime }$，确定$H$的一种简单线性算法。

变换方程由$x_i^{\prime }=Hx$给出，由于等号两边相差一个比例系数，改为$x_i^{\prime }\times H{x}_i=0$

记$H$的$j$行为$h^{jT}$
$$H{x}_i=\left[\begin{array}{c}{h}^{1T}{x}_i\\ h^{2T}{x}_i\\ h^{3T}{x}_i\end{array}\right]$$
记$x_i^{\prime }=(x_i^{\prime },y_i^{\prime },w_i^{\prime })$，则叉积可表示为
$$x_i^{\prime }\times H{x}_i=\left[\begin{array}{c}{y}_i^{\prime }{h}^{3T}{x}_i-w_i^{\prime }{h}^{2T}{x}_i\\ w_i^{\prime }{h}^{1T}{x}_i-x_i^{\prime }{h}^{3T}{x}_i\\ x_i^{\prime }{h}^{2T}{x}_i-y_i^{\prime }{h}^{1T}{x}_i\end{array}\right]$$
因为$h^{jT}{x}_i=x_i^T{h}^j$，有
$$\left[\begin{array}{ccc}{0}^T& -w_i^{\prime }{x}_i^T& y_i^{\prime }{x}_i^T\\ -w_i^{\prime }{x}_i^T& 0^T& -x_i^{\prime }{x}_i^T\\ -y_i^{\prime }{x}_i^T& x_i^{\prime }{x}_i^T& 0^T\end{array}\right]\left[\begin{array}{c}{h}^1\\ h^2\\ h^3\end{array}\right]=A_{i}h=0$$
在$(39)$中，由于方程只有两个是线性独立的，因此可以只使用前两个公式
$$\left[\begin{array}{ccc}{0}^T& -w_i^{\prime }{x}_i^T& y_i^{\prime }{x}_i^T\\ -w_i^{\prime }{x}_i^T& 0^T& -x_i^{\prime }{x}_i^T\end{array}\right]\left[\begin{array}{c}{h}^1\\ h^2\\ h^3\end{array}\right]=A_{i}h=0$$

超定解

当数据存在误差且数据较多时，在通常约束$\parallel h\parallel =1$下最小化范数$\parallel Ah\parallel =1$，即求解$\parallel Ah\parallel /\parallel h\parallel$的最小值，也就是求$A$最小奇异值的单位奇异向量。

至此给出$DLT$算法，给定$n\ge 4$组$2D-2D$点对应 { x i    ⟺    x i ′ } \{\pmb x_i \iff \pmb x_i^\prime\} {xi​⟺xi′​}，确定$2D$单应矩阵使得$x_i^{\prime }=H{x}_i$

1. 根据每组对应计算矩阵$A_i$，只使用前两行
2. 把$n$个$2\times 9$的矩阵$A_i$合成一个$2n\times 9$的矩阵$A$
3. 求$A$的$SVD$。对应于最小特征值的单位特征向量就是解$h$
4. $H$由$h$确定

非齐次解

除去把$h$直接作为齐次向量来解，还有一种办法是转换成非齐次线性方程组，即对向量$h$中的某个元素硬性加上$h_j=1$的条件。如$h_9=1$，则有
$$A_i{h}_{1-8}=\left[\begin{array}{c}-w_i{y}_i^{\prime }\\ w_i{x}_i^{\prime }\end{array}\right]$$
但是如果真正的解是$h_j=0$，则会出现数值计算法上的不稳定解。

退化配置

出现某种配置不能确定唯一解的情形称为退化

不同的代价函数

介绍为确定$H$的超定解而最小化的若干代价函数

代数距离

DLT算法最小化范数$\parallel Ah\parallel$，向量$\epsilon =Ah$称为残差向量。每组对应$x_i⟺x_i^{\prime }$贡献部分误差向量${\epsilon }_i$并积累成总误差向量$\epsilon$。该向量的范数是一个标量，称为代数距离
$$d_{alg}(x_i^{\prime },H{x}_i{)}^2=\parallel {\epsilon }_i{\parallel }^2$$
其中
$$d_{alg}(x_1,x_2{)}^2=a_1^2+a_2^2,a=(x_1\times x_2)$$

几何距离

几何距离是基于图像上几何距离的测量并最小化图像坐标的测量值和估计值之差。用$x$表示测量得到的图像坐标；$\hat{x}$表示点的估计值，$\overline{x}$表示点的真值

• 单幅图像误差

  首先考虑第一幅图像的测量值是非常准确的，误差仅仅出现在第二幅图像。$d(x,y)$表示点$x,y$的欧氏距离
  $$\sum _{i}d(x_i^{\prime },H{\overline{x}}_i{)}^2$$

• 对称转移误差

  误差在两幅图像中同时发生
  $$\sum _{i}d(x_i,H^{-1}{x}_i^{\prime }{)}^2+d(x_i^{\prime },H{x}_i{)}^2$$

重投影误差

寻找一个单应$\hat{H}$和完全匹配的点${\hat{x}}_i$和${\hat{x}}_i^{\prime }$以最小化总的误差函数
$$\sum _{i}d(x_i,{\hat{x}}_i{)}^2+d(x_i^{\prime },{\widehat{x_i}}^{\prime }{)}^2,s.t.{\hat{x}}_i^{\prime }=\hat{H}{\hat{x}}_i,\forall i$$

几何和代数距离比较

假设误差仅出现在第二幅图像，令$x_i^{\prime }=(x_i^{\prime },y_i^{\prime },w_i^{\prime }{)}^T$并定义向量$({\hat{x}}_i^{\prime },{\hat{y}}_i^{\prime },{\hat{w}}_i^{\prime }{)}^T={\hat{x}}_i^{\prime }=H\overline{x}$，$x^{\prime }\times Hx=x^{\prime }\times H\overline{x}=x^{\prime }\times {\hat{x}}^{\prime }$，变为
$$A_{i}h={\epsilon }_i=\left[\begin{array}{c}{y}_i^{\prime }{\hat{w}}_i^{\prime }-w_i^{\prime }{\hat{y}}_i^{\prime }\\ w_i^{\prime }{\hat{x}}_i^{\prime }-x_i^{\prime }{\hat{w}}_i^{\prime }\end{array}\right]$$
得到
$$d(x_i^{\prime },{\hat{x}}_i^{\prime })=d_{alg}(x_i^{\prime },{\hat{x}}_i)/{\hat{w}}_i^{\prime }{w}_i^{\prime }$$
仿射变换下，几何距离和代数距离相等。DLT算法很容易做到强制$H$的最后一行为$(0,0,1)$，因此对于仿射变换，最小化几何距离可以用基于代数距离的线性DLT算法。

重投影误差的几何解释

两平面之间的单应估计被视为$4D$空间$I{R}^4$中用曲面来拟合点，$x^{\prime }\times Hx=0$定义了一个代数族${\mathcal{V}}_H$，它是两个二次超平面的交集。$X=(x,y,x^{\prime },y^{\prime }{)}^T$，寻找最合适的${\mathcal{V}}_H$，使得点到平面的距离最小，$\hat{X}$是曲面上最接近$X$的点
$$\parallel X_i-{\hat{X}}_i{\parallel }_2=d_{\perp }(X_i,{\mathcal{V}}_H)=d(x_i,\widehat{x_i}{)}^2+d(x_i^{\prime },{\widehat{x_i}}^{\prime }{)}^2$$
为了可视化理解，考虑在二次曲线拟合平面上的$n>5$个点，使得基于几何距离的误差最小

Sampson误差

思想是估计点$\hat{X}$的一阶近似，对于给定的单应$H$，任何在${\mathcal{V}}_H$上的点$X=(x,y,x^{\prime },y^{\prime }{)}^T$满足$Ah=0$，记为$C_H(X)=0$。代价函数可以表示为
$$C_H(X+{\delta }_X)=C_H(X)+\frac{\partial C_H}{\partial X}{\delta }_X$$

记${\delta }_X=\hat{X}-X$并希望$\hat{X}$在族${\mathcal{V}}_H$上，使得$C_H(\hat{X})=0$，便得到$C_H(X)+(\partial C_H/\partial X){\delta }_X=0$。记$J{\delta }_X=-\epsilon$，其中$J$是偏导，$\epsilon$是与$X$相关的代价函数$C_H(X)$。问题就变成了求解满足此方程的最小${\delta }_X$。即

求在满足$J{\delta }_X=-\epsilon$条件下使$\parallel {\delta }_X\parallel$取最小值的向量${\delta }_X$

求解方法是使用拉格朗日乘子，问题就变成了最小化${\delta }_X^T{\delta }_X-2\lambda (J{\delta }_X+\epsilon )$

对${\delta }_X$求导并使其等于0，得到
$$2{\delta }_X^T-2{\lambda }^{T}J=0^T$$
对$\lambda$求导则给出
$$J{\delta }_X+\epsilon =0$$
消去${\delta }_X$，得到
$$J{J}^T\lambda =-\epsilon$$
对$\lambda$求解得$\lambda =-(J{J}^T{)}^{-1}\epsilon$，最后得到
$${\delta }_X=-J^T(J{J}^T{)}^{-1}\epsilon$$
范数$\parallel {\delta }_X{\parallel }^2$是Sampson误差

通用几何解释框架

在通用框架中，估计问题包括两部分

• 一个由测量向量$X$组成的测量空间$I{R}^N$
• 一个模型，该模型视为$I{R}^N$中点的一个子集$S$，在此子集内的一个测量向量$X$称为满足此模型。通常，满足此模型的子空间是一个子流形或$I{R}^N$中的代数族

给定$I{R}^N$中的一个测量向量$X$，估计问题就是求一个离$X$最近并满足该模型的向量$\hat{X}$

统计代价函数和最大似然估计

真值是$\overline{x}$，设每个测量点$x$的概率密度函数是
$$\mathrm{Pr}(x)=(\frac{1}{2\pi {\sigma }^2})e^{-d}(x,\overline{x})/(2{\sigma }^2)$$

单图像误差

因为假设每点的误差是独立分布的，受噪声干扰的数据的$PDF$是：
Pr ⁡ ( { x i ′ } ∣ H ) = ∏ i ( 1 2 π σ 2 ) e − d ( x i ′ , H i x ‾ i ) 2 / ( 2 σ 2 ) \operatorname{Pr}\left(\left\{\mathbf{x}_i^{\prime}\right\} \mid \mathrm{H}\right)=\prod_i\left(\frac{1}{2 \pi \sigma^2}\right) e^{-d\left(\mathbf{x}_i^{\prime}, \mathrm{H}_i \overline x_i\right)^2 /\left(2 \sigma^2\right)} Pr({xi′​}∣H)=i∏​(2πσ21​)e−d(xi′​,Hi​xi​)2/(2σ2)
Pr ⁡ ( { x i ′ } ∣ H ) \operatorname{Pr}\left(\left\{\mathbf{x}_i^{\prime}\right\} \mid \mathrm{H}\right) Pr({xi′​}∣H)表示获得真实$H$的情况下 { x i ′ } \{x_i^\prime\} {xi′​}的概率，上述公式的对数似然为
log ⁡ Pr ⁡ ( { x i ′ } ∣ H ) = − 1 2 σ 2 ∑ i d ( x i ′ , H x ‾ i ) 2 +  constant  \log \operatorname{Pr}\left(\left\{\mathbf{x}_i^{\prime}\right\} \mid \mathrm{H}\right)=-\frac{1}{2 \sigma^2} \sum_i d\left(\mathbf{x}_i^{\prime}, \mathrm{H} \overline{\mathbf{x}}_i\right)^2+\text { constant } logPr({xi′​}∣H)=−2σ21​i∑​d(xi′​,Hxi​)2+ constant 
单应的最大似然估计$ML$估计$H$就是最大化对数似然，即最小化
$$\sum _{i}d{\left({\mathbf{x}}_i^{\prime },\mathrm{H}{\overline{\mathbf{x}}}_i\right)}^2$$
$ML$估计等价于最小化几何误差函数

双图像误差

同样已知 { x ‾ i ↔ H x ‾ i = x ‾ i ′ } \left\{\overline{\mathbf{x}}_i \leftrightarrow \mathrm{H} \overline{\mathbf{x}}_i=\overline{\mathbf{x}}_i^{\prime}\right\} {xi​↔Hxi​=xi′​}，数据的PDF是
Pr ⁡ ( { x i , x i ′ } ∣ H , { x ‾ i } ) = ∏ i ( 1 2 π σ 2 ) e − ( d ( x i , x ‾ i ) 2 + d ( x i ′ , H i ) 2 ) / ( 2 σ 2 ) \operatorname{Pr}\left(\left\{\mathbf{x}_i, \mathbf{x}_i^{\prime}\right\} \mid \mathrm{H},\left\{\overline{\mathbf{x}}_i\right\}\right)=\prod_i\left(\frac{1}{2 \pi \sigma^2}\right) e^{-\left(d\left(\mathbf{x}_i, \overline{\mathbf{x}}_i\right)^2+d\left(\mathbf{x}_i^{\prime}, \mathrm{H}_i\right)^2\right) /\left(2 \sigma^2\right)} Pr({xi​,xi′​}∣H,{xi​})=i∏​(2πσ21​)e−(d(xi​,xi​)2+d(xi′​,Hi​)2)/(2σ2)
类似的$ML$估计是
$$\sum _{i}d{\left({\mathbf{x}}_i,{\overline{\mathbf{x}}}_i\right)}^2+d{\left({\mathbf{x}}_i^{\prime },{\mathrm{H}}_i\right)}^2$$

Mahalanobis距离

在一般高斯分布的情况下，假定测量向量$X$满足一个具有协方差矩阵$\Sigma$地高斯分布函数。

最大化对数似然等价于最小化$Mahalanobis$距离
$$\parallel X-\overline{X}{\parallel }_{\Sigma }^2=(X-\overline{X}{)}^T{\Sigma }^{-1}(X-\overline{X})$$
当每幅图像都有误差并假定一幅图像的误差与另一幅图像的误差是独立的时候，代价函数是
$$\parallel X-\overline{X}{\parallel }_{\Sigma }^2+\parallel X^{\prime }-{\overline{X}}^{\prime }{\parallel }_{{\Sigma }^{\prime }}^2$$
最后，假定所有点的误差是独立的，则有
$$\sum \parallel x_i-{\overline{x}}_i{\parallel }_{{\Sigma }_i}^2+\parallel x_i^{\prime }-{\overline{x}}_i^{\prime }{\parallel }_{{\Sigma }_i^{\prime }}^2$$

变换不变性和归一化

假设在算法运行之前，图像坐标由某个相似、仿射甚至射影进行了变换，其结果是否会有实质性改变。

假设一幅图像的坐标$x$被$\tilde{x}=Tx$替代，另一幅图像的坐标$x^{\prime }$被${\tilde{x}}^{\prime }=T^{\prime }{x}^{\prime }$替代，代入等式$x^{\prime }=Hx$，得到${\tilde{x}}^{\prime }=T^{\prime }H{T}^{-1}\tilde{x}$，求出$\tilde{H}=T^{\prime }H{T}^{-1}$。因此，求把$x_i$映射到$x_i^{\prime }$的变换的另一种方法是

• 根据公式$\tilde{x}=Tx,{\tilde{x}}^{\prime }=T^{\prime }{x}^{\prime }$分别对图像进行变换
• 找到${\tilde{x}}_i\leftrightarrow {\tilde{x}}_i^{\prime }$找到对应的$\tilde{H}$
• 令$H=T^{\prime -1}\tilde{H}T$

现在考虑选择不同的变换，是否会对算法的结果有影响。可以证明，

• 最小化几何误差的算法关于相似变换不变
• DLT算法不是相似变换不变的，补救方法是在进行DLT算法之前进行归一化变换。归一化变换将消除由任意选取图像坐标系的原点和尺度所产生的影响，这种组合变换对图像的相似变换不变

求2D单应的归一化的DLT

目标

给定$n\ge 4$组$2D$到$2D$的点对应 { x i ↔ x i ′ } \{x_i \leftrightarrow x_i^\prime\} {xi​↔xi′​}，确定$2D$单应矩阵$H$使得$x_i^{\prime }=H{x}_i$

算法

1. 归一化$x$：计算一个只包含位移和缩放的相似变换$T$，将点$x_i$变到新的点集${\tilde{x}}_i$，使得${\tilde{x}}_i$的形心位于原点，并且它们到原点的平均距离是$\sqrt{2}$
2. 归一化$x^{\prime }$
3. $\mathrm{DLT}$：根据${\tilde{x}}_i\leftrightarrow {\tilde{x}}_i^{\prime }$，求得$\tilde{H}$
4. 解除归一化：令$H=T^{\prime -1}\tilde{H}T$

迭代最小化方法

迭代最小化技术由五步组成：

1. 代价函数：代价函数是最小化的基础
2. 参数化：把要计算的变换表示成有限数目的参数。一般不要求它是最小参数集，事实上超参数化常有优越性（使用最小化参数会使代价函数曲面变得更复杂，程序在局部最小处停止不前的可能性更大）
3. 函数确定：必须确定一个用参数集描述的代价函数
4. 初始化：计算一个适当的初始参数估计。一般由一个线性算法来实现
5. 迭代：由初始解开始，在迭代中逐步修正参数以达到最小化迭代函数的目的

鲁棒性估计

随机抽样一致算法$\text{Random Sample Consensus,RANSAC}$能应对大比例的野值，思路是通过模型去拟合数据（这里是$H$矩阵），随机样本包含足以确定模型的最小数据集（求单应是4组点）。算法流程如下：

目标

一个模型与含有野值的数据集$S$的鲁棒拟合

算法

1. 随机地从$S$中选择$s$个数据点组成的一个样本作为模型的一个示例
2. 确定在模型距离阈值$t$内的数据点集$S_i$，$S_i$称为采样的一致集并定义$S$的内点
3. 如果$S_i$的大小大于某个阈值$T$，用$S_i$的所有点重新估计模拟并结束
4. 如果$S_i$的大小小于$T$，选择一个新的子集并重复上述流程
5. 经过$N$次试验选择最大一致集$S_i$，并用$S_i$的所有点重估计模型

采样次数的计算，$s$是采样大小，$\epsilon$是野值的概率,$s$个点组成的随机样本中至少有一次没有野值的概率为$p$
$$N=\log (1-p)/\log (1-(1-\epsilon {)}^s)$$