计算机视觉中的多视图几何＜Part0—基础知识：射影几何、变换和估计＞(3)

3. 估计——2D射影变换

估计是指在某些本质测量的基础上计算某个变换或其他数学量。

• 测量数：考虑自由度个数和约束个数后给出的一个下届；
• 近似解：精确解称为最小配置解，近似解是通过最小化某个代价函数完成；
• 黄金标准算法：通常存在一种最优代价函数，其最优的含义是在一定假设下，使代价函数取最小值的解是最好估计，计算该代价函数最小值的算法称为黄金标准算法；

3.1 直接线性变换(DLT)算法

齐次矢量方程，相差一个尺度因子：
$${\text{x}}_i^{\prime }=H{\text{x}}_i\to {\text{x}}_i^{\prime }\times H{\text{x}}_i=0\to A_i\text{h}=0$$
(1) $A_i\text{h}=0$是未知矢量$\text{h}$的线性方程，矩阵$A_i$的元素是已知点坐标的二次多项式；
(2) $A_i$中只有两行线性独立(包含无穷远点)，$A_i$为$2\times 9$矩阵；
(3) 该方程组对任何齐次坐标成立，非零因子可以通过对范数要求任意选择，$||\text{h}||=1$；

3.1.1 Over-determined solution

对应点多余最小配置，$A\text{h}=0$为超定问题；由于测量存在噪声，可分为零解和近似解两种情况，为避免产生零解($\text{h}=\text{0}$)，附加范数条件$||\text{h}||=1$，等价于最小化$||A\text{h}||\to ||A\text{h}||/||\text{h}||$最小值问题，解是$||A^{T}A||$的最小特征值的(单位)特征矢量，即$A$最小奇异值的单位奇异向量；

3.1.2 Inhomogeneous solution

强加条件$h_j=1$，既然允许相差一个任意因子，可以通过因子选择使得$h_j=1$：
$$A_i\text{h}=0\to M_i\tilde{\text{h}}=\text{b}$$
$\tilde{\text{h}}$是8维矢量，可用高斯消去(线性方程基本解法)或最小二乘法求解，若真正的解$h_j$趋近于零，此方法导致不稳定解；

3.1.3 Degenerate configurations

对于具体的一类变换而言，出现某种配置(共线点或共点线)不能确定唯一解的情形称为退化(无穷解簇或无解)。

3.1.4 Solutions from lines and other entities

单应也可以由直线或二次曲线的对应来计算，约束数必须等于或大于变换的自由度数，混合配置需注意配置可产生的约束数目；

3.2 不同的代价函数

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

3.2.1 Algebraic distance

DLT算法最小化范数$||A\text{h}||$，矢量$\epsilon =A\text{h}$称为残差矢量；矢量${\epsilon }_i$是关联对应点对和单应$H$的代数误差矢量，该矢量的范数是一个标量，称为代数距离：
$$\sum _i{d}_{alg}({\text{x}}_i^{\prime },H{\text{x}}_i{)}^2=\sum _i||{\epsilon }_i|{|}^2=||A\text{h}|{|}^2=||{\epsilon }^2||$$
对于任何两个矢量${\text{x}}_1$和${\text{x}}_2$的代数距离为：$d_{alg}({\text{x}}_1,{\text{x}}_2{)}^2=a_1^2+a_2^2$，其中$\text{a}=(a_1,a_2,a_3{)}^{\top }={\text{x}}_1\times {\text{x}}_2={\epsilon }_i$；其特殊的优越性在于线性(唯一)解以及计算代价小，通常代数距离的解被用作几何或统计代价函数的非线性最小化的起点，非线性最小化给这个解一个最后的“抛光”；

3.2.2 Geometric distance

基于图像上几何距离的测量并最小化图像坐标的测量值与估计值之差：矢量$\text{x}$表示测量的图像坐标，$\hat{\text{x}}$表示该点的估计值，而$\bar{\text{x}}$表示该点的真值。

• 单幅图像误差：误差仅出现在第二幅图像的情形(标定)，最小化转移误差${\sum }_{i}d({\text{x}}_i^{\prime },H{\bar{\text{x}}}_i{)}^2$；
• 对称转移误差：${\sum }_{i}d({\text{x}}_i,H^{-1}{\text{x}}_i^{\prime }{)}^2+{\sum }_{i}d({\text{x}}_i^{\prime },H{\text{x}}_i{)}^2$；
• 对应集合的转移误差是算法要估计的单应使误差取最小值的单应；

3.2.3 Reprojection error – both images

估计每组对应的“校正值”(对测量值校正)：
$$\sum _{i}d({\text{x}}_i,{\hat{\text{x}}}_i{)}^2+d({\text{x}}_i^{\prime },{\hat{\text{x}}}_i^{\prime }{)}^{2}s.t.{\hat{\text{x}}}_i^{\prime }=\hat{H}{\hat{\text{x}}}_i$$
这种估计可对诸如世界平面点的图像对应点${\text{x}}_i\leftrightarrow {\text{x}}_i^{\prime }$测量建模，估计世界平面点${\hat{\text{X}}}^{\prime }$，然后将它重投影到估计$H$上认为是完全匹配的对应${\hat{\text{x}}}_i\leftrightarrow {\hat{\text{x}}}_i^{\prime }$；

3.2.4 Comparison of geometric and algebraic distance

几何距离与代数距离相关，但不相等：
$$d({\text{x}}_i^{\prime },{\hat{\text{x}}}_i^{\prime })=d_{alg}({\text{x}}_i^{\prime },{\hat{\text{x}}}_i^{\prime })/{\hat{w}}_i^{\prime }{w}_i^{\prime }$$
仿射变换下(齐次不变性)，几何距离和代数距离相等，对于仿射变换，最小化几何距离可以用基于代数距离的线性DLT算法；

3.2.5 Geometric interpretation of reprojection error

两平面之间的单应估计可以视为4D空间IR⁴中的“曲面”来拟合点：每对图像点定义测量空间IR⁴的一个点，满足${\text{x}}^{\prime }\times H\text{x}$的点对定义了IR⁴中的代数簇${\mathcal{V}}_H$(IR^N中一个或多个多元多项式公共零点集，是两个超二次曲面的交集)，$H$中的元素确定多项式每一项系数，即定义了一个具体的二次曲面；给定IR⁴中的点${\text{X}}_i=(x_i,y_i,x_i^{\prime },y_i^{\prime }{)}^{\top }$，单应估计的任务变成寻找一个通过点(或几乎通过)${\text{X}}_i$的某个簇${\mathcal{V}}_H$(对应变换$H$)，对每一点${\text{X}}_i$，令${\hat{\text{X}}}_i=({\hat{x}}_i,{\hat{y}}_i,{\hat{x}}_i^{\prime },{\hat{y}}_i^{\prime }{)}^{\top }$为簇${\mathcal{V}}_H$上最靠近${\text{X}}_i$的点；IR⁴中的几何距离等价于两图像中的重投影误差，求${\mathcal{V}}_H$和其上到测量点距离平方和最小的点，等价于求单应和最小化重投影误差函数；

• 在${\mathcal{V}}_H$上离测量点$\text{X}$最近的点$\hat{\text{X}}$与$\text{X}$的连线垂直于过$\hat{\text{X}}$的切平面：
  $$d({\text{x}}_i,{\hat{\text{x}}}_i{)}^2+d({\text{x}}_i^{\prime },{\hat{\text{x}}}_i^{\prime }{)}^2=d_{\perp }({\text{X}}_i,{\mathcal{V}}_H{)}^2$$
  其中$d_{\perp }({\text{X}}_i,{\mathcal{V}}_H)$是点$\text{X}$到簇${\mathcal{V}}_H$的垂直距离。
• 二次曲线类似：
  
  在仿射变换时，该代数簇是两张超平面的交，即是一个2维线性子空间，${\text{x}}^{\prime }=H_A\text{x}$给出$x,x^{\prime },y$和$x,y,y^{\prime }$两个线性约束，每个约束定一个IR⁴中的一个超平面，关系都可以通过最小化点到该代数簇的垂直距离来估计。

3.2.6 Sampson error

几何误差的性质相当复杂，其最小化需要同时估计单应矩阵和满足单应矩阵的对应点对，引入Sampson误差介于几何与代数误差之间；由于代数簇${\mathcal{V}}_H$的非线性本质，点$\text{X}$不能直接估计$\hat{\text{X}}$，必须通过迭代；Sampson误差函数的思想是估计点$\hat{\text{X}}$的一阶近似并假定代价函数在被估计点附近有很好的线性近似；对给定单应$H$，在${\mathcal{V}}_H$上的任何点$\text{X}$满足$A\text{h}=0$；为突出代价函数对$\text{X}$的相关性，记代价函数(二维矢量)${\mathcal{C}}_H(\text{X})=0$，代价函数一阶Taylor逼近如下：
$${\mathcal{C}}_H(\text{X}+{\delta }_{\text{x}})={\mathcal{C}}_H(\text{X})+\frac{\partial {\mathcal{C}}_H}{\partial \text{X}}{\delta }_{\text{x}}=\text{0}$$
$${\delta }_{\text{x}}=\hat{\text{X}}-\text{X}\to {\mathcal{C}}_H(\text{X})+\frac{\partial {\mathcal{C}}_H}{\partial \text{X}}{\delta }_{\text{x}}=\text{0}\to J{\delta }_{\text{x}}=-\epsilon$$
在满足$J{\delta }_{\text{x}}=-\epsilon$条件下使$||{\delta }_{\text{x}}||$取最小值的矢量${\delta }_{\text{x}}$，引入Langrange乘子：
$${\delta }_{\text{x}}^{\top }{\delta }_{\text{x}}-2{\lambda }^{\top }(J{\delta }_{\text{x}}+\epsilon )\to {\delta }_{\text{x}}=J^{\top }\lambda \&\lambda =-(J{J}^{\top }{)}^{-1}\epsilon \to ||{\delta }_{\text{x}}|{|}^2={\epsilon }^{\top }(J{J}^{\top }{)}^{-1}\epsilon$$
范数$||{\delta }_{\text{x}}|{|}^2$是Sampson误差；其中$\epsilon$和$J$都与$H$有关，变参数集仅由$H$的元素组成。
注解：
(1) 2D单应估计中的$\text{X}=(x,y,x^{\prime },y^{\prime }{)}^{\top }$，其中2D测量是$\text{x}=(x,y,1{)}^{\top }$和${\text{x}}^{\prime }=(x^{\prime },y^{\prime },1{)}^{\top }$.
(2) $\epsilon ={\mathcal{C}}_H(\text{X})$是代数误差矢量$A_i\text{h}$.
(3) $J=\partial {\mathcal{C}}_H(\text{X})/\partial (\text{X})$是一个$2\times 4$的矩阵：$J_{11}=\partial (-w_i^{\prime }{\text{x}}_i^{\top }{\text{h}}^2+y_i^{\prime }{\text{x}}_i^{\top }{\text{h}}^3)/\partial x=-w_i^{\prime }{h}_{21}+y_i^{\prime }{h}_{31}$.
(4) Sampson误差与代数误差的相似性，可以解释成Mahalanobis范数.
(5)这里推导Sampson误差针对一组点，若干组点估计单应时误差相加：${\mathcal{D}}_{\perp }={\sum }_i{\epsilon }_i^{\top }(J_i{J}_i^{\top }{)}^{-1}{\epsilon }_i$
线性代价函数：代数误差矢量${\mathcal{C}}_H(\text{X})=A(\text{X})\text{h}$关于$\text{X}$元素是典型多重线性的，由Taylor展开可以给出几何误差的精确一阶近似，表明Sampson误差等同于几何误差；另外，由一组线性方程定义的代数簇是$H$的超平面，求$H$的变成超平面拟合问题。

3.2.7 Another geometric interpretation

• 一个由测量矢量$\text{X}$组成的测量空间IR^N；
• 一个模型，该模型视为IR^N中点的一个子集，测量矢量$\text{X}$满足此模型，该点集将组成IR^N中的某维数流形S，该子流形维数(自由度)等于用来参数化该流形的最少参数数；

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

• 每个测量点$\text{x}$的概率密度函数(pdf)：
  $$Pr(\text{x})=(\frac{1}{2\pi {\sigma }^2})e^{-d(\text{x},\bar{\text{x}}{)}^2/(2{\sigma }^2)}$$
• 单图像误差：给定真实$H$时获得测量 { x i ′ } \{\textbf{x}'_i\} {xi′​}的概率(独立同分布)，该对应集合的对数似然为：
  P r ( { x i ′ } ∣ H ) = ∏ i ( 1 2 π σ 2 ) e − d ( x i ′ , H x ˉ i ) 2 / ( 2 σ 2 ) Pr(\{\textbf{x}'_i\}|H)=\prod_i(\frac{1}{2\pi\sigma^2})e^{-d(\textbf{x}'_i,H\bar{\textbf{x}}_i)^2/(2\sigma^2)} Pr({xi′​}∣H)=i∏​(2πσ21​)e−d(xi′​,Hxˉi​)2/(2σ2)
  log ⁡ P r ( { x i ′ } ∣ H ) = − 1 2 σ 2 ∑ i d ( x i ′ , H x ˉ i ) 2 + c o n s t \log{Pr}(\{\textbf{x}'_i\}|H)=-\frac{1}{2\sigma^2}\sum_id(\textbf{x}'_i,H\bar{\textbf{x}}_i)^2+const logPr({xi′​}∣H)=−2σ21​i∑​d(xi′​,Hxˉi​)2+const
  单应的最大似然(ML)估计$\hat{H}$最大化对数似然，即最小化几何误差函数；
• 双图像误差：
  P r ( { x i , x i ′ } ∣ H , x ˉ i ′ ) = ∏ i ( 1 2 π σ 2 ) e − ( d ( x i , x ˉ i ) 2 + d ( x i ′ , H x ˉ i ) 2 ) / ( 2 σ 2 ) Pr(\{\textbf{x}_i,\textbf{x}'_i\}|H,\bar{\textbf{x}}'_i)=\prod_i(\frac{1}{2\pi\sigma^2})e^{-(d(\textbf{x}_i,\bar{\textbf{x}}_i)^2+d(\textbf{x}'_i,H\bar{\textbf{x}}_i)^2)/(2\sigma^2)} Pr({xi​,xi′​}∣H,xˉi′​)=i∏​(2πσ21​)e−(d(xi​,xˉi​)2+d(xi′​,Hxˉi​)2)/(2σ2)
  ML估计等同于最小化重投影误差函数；
• Mahalanobis距离：在一般高斯分布的情形，可以假定测量矢量$X$满足一个具有协方差矩阵$\Sigma$的高斯分布函数，最大化对数似然则等价于最小化Mahalanobis距离
  $$||X-\bar{X}|{|}_{\Sigma }^2=(X-\bar{X}{)}^{\top }{\Sigma }^{-1}(X-\bar{X})$$
  图像间误差独立：
  $$||X-\bar{X}|{|}_{\Sigma }^2+||X^{\prime }-{\bar{X}}^{\prime }|{|}_{\Sigma }^2$$
  点之间误差独立：
  $$\sum ||{\text{x}}_i-{\bar{\text{x}}}_i|{|}_{{\Sigma }_i}^2+\sum ||{\text{x}}_i^{\prime }-{\bar{\text{x}}}_i^{\prime }|{|}_{\Sigma }^2$$
  这个方程允许用于非各向同性的协方差矩阵，使用两非垂直线相交计算点的位置会出现此情况。

3.4 变换不变性和归一化

评价算法与图像坐标系原点、尺度甚至定向选择之间的关系。

3.4.1 Invariance to image coordinate transformations

$${\tilde{\text{x}}}^{\prime }=T^{\prime }H{T}^{-1}\tilde{\text{x}}$$
(1) 根据公式${\tilde{\text{x}}}_i=T{\text{x}}_i$和${\tilde{\text{x}}}_i^{\prime }=T^{\prime }{\text{x}}_i^{\prime }$变换图像坐标；
(2) 由${\tilde{\text{x}}}_i\leftrightarrow {\tilde{\text{x}}}_i^{\prime }$求变换$\tilde{H}$；
(3) 令$H=T^{\prime -1}\tilde{H}T$；
最小化几何误差的算法关于相似变换不变，DLT算法结果却不是相似变换不变，归一化变换将消除DLT由任意选取图像坐标系原点和尺度所产生的影响。

3.4.2 Non-invariance of the DLT algorithm

令$T^{\prime }$为具有缩放因子$s$的相似变换，$T$为任意的射影变换；此外，假设$H$是任何2D单应并定义$\tilde{H}=T^{\prime }H{T}^{-1}$；那么$||\tilde{A}\tilde{\text{h}}||=s||A\text{h}||$(余因子式：$(M\text{x})\times (M\text{y})=M^{\ast }(\text{x}\times \text{y})$)，然后$||H||=1$不等价于$||\tilde{H}||=1$(约束条件不同），因此变换矩阵的解不同。

3.4.3 Invariance of geometric error

$$d({\tilde{\text{x}}}^{\prime },\tilde{H}\tilde{\text{x}})=d(T^{\prime }{\text{x}}^{\prime },T^{\prime }H{T}^{-1}T\text{x})=d(T^{\prime }{\text{x}}^{\prime },T^{\prime }H\text{x})=d({\text{x}}^{\prime },H\text{x})$$

3.4.4 Normalizing transformations

归一化使代数最小化在一个固定的标准坐标系中进行，因而使DLT算法实际上关于相似变换不变，归一化影响DLT方程组的条件数或准确的说与方程组矩阵A的第一个和倒数第二个奇异值的比率$d_1/d_{n-1}$有关，非各向同性缩放和无穷远附近的点的缩放需特殊考虑。

3.5 迭代最小化方法

• 几何代价函数迭代最小化方法：
  • 代价函数：几何误差、重投影误差、Sampson误差、Mahalanobis距离；
  • 参数化：把要计算的变换(或其他实体)表示成有限数目的参数，不推荐最小参数化方法(代价函数曲面复杂，易陷入局部最小值)，超参数化(优越性)需结合变换具体类型使用；
  • 函数确定：必须确定一个用参数集$P$描述的代价函数；
    • 一个协方差矩阵为$\Sigma$的测量矢量$X\in I{R}^N$
    • 一组参数被表达成一个矢量$P\in I{R}^M$
    • 映射$f:I{R}^M\to I{R}^N$，容许测量集的模型曲面S
    • 最小化代价函数是Mahalanobis距离的平方$||X-f(P)|{|}_{\Sigma }^2$，Sampson近似使重投影误差也能仅对9个参数最小化，LM方法求解m参数的非线性最小化问题，复杂度为$O(m^3)$；
  • 初始化：计算一个适当的初始参数估计(Normalized-DLT & RANSAC)，参数空间密集采样或固定参数初始化；
  • 迭代方法：由初始解开始，在迭代中逐步修正参数以达到最小化代价函数的目的；
    

3.6 鲁棒估计

3.6.1 RANSAC

• 距离阈值：
  
  
• 采样次数：
  
• 一致集大小：
  $$T=(1-\epsilon )n$$
• 自适应采样次数：
  

3.6.2 Robust MLE

最大一致性集合与LM最小初始集冲突，采用鲁棒代价函数在所有数据上最小化：

3.7 单应自动计算

3.8 黄金标准算法

3.9 总结

1. RANSAC用于剔除外点；
2. 几何误差迭代为了减小噪声；
3. DLT的变换不变性：$H_{33}=1$旋转缩放变换下不变但平移变，$H_{31}^2+H_{32}^2=1$相似变换下不变，$H_{31}=H_{32}=0,H_{33}=1$仿射变换不变；
4. 缩放无界点集：
   $$\sum _i{x}_i=\sum _i{y}_i=0;\sum _i{x}_i+\sum _i{y}_i=2;x_i^2+y_i^2+w_i^2=1,\forall i$$
5. 具有非各向同性误差分布的Sampson误差：
   $$||{\delta }_{\text{x}}|{|}_{{\Sigma }_{\text{x}}}^2={\delta }_{\text{x}}^{\top }{\Sigma }_{\text{x}}^{-1}{\delta }_{\text{x}}={\epsilon }^{\top }(J{\Sigma }_{\text{x}}{J}^{\top }{)}^{-1}\epsilon$$
   $${\delta }_{\text{x}}=-{\Sigma }_{\text{x}}{J}^{\top }(J{\Sigma }_{\text{x}}{J}^{\top }{)}^{-1}\epsilon$$
6. 最小化仿射变换的几何误差：
   a. 最佳仿射变换把${\text{x}}_i$的形心映射到${\text{x}}_i^{\prime }$的形心，可通过点的平移使形心移动到原点，变换的平移部分就确定了；
   b. 点$X_i=({\text{x}}_i^{\top },{\text{x}}_i^{\prime \top }{)}^{\top }$在${\mathcal{V}}_H$上满足充要条件是$[H_{2\times 2}|-I_{2\times 2}]X=\text{0}$，因此${\mathcal{V}}_H$是IR⁴中余维度为2的子空间；
   c. 给定行为$X_i^{\top }$的矩阵$M$，$X_i$的最好拟合子空间由对应$M$的两个最大奇异值的特征矢量$V_1$和$V_2$生成，$[H_{2\times 2}|-I][V_1{V}_2]$；