MVG读书笔记——单应矩阵估计这件小事（二）

还是以两幅图像进行单应矩阵求解为例，上面讲到使用DLT算法一对对应点之间可以构成一个方程组$Ah=0$，其中A为一个$2 \times 9$的矩阵。由此只需要4个点就可以求解出H矩阵。但是在实际的应用中，还有一些问题需要解决。

超定方程

实际情况下，我们得到的两幅图像的对应点往往多于4个点，从而可以得到矩阵$A_{2n \times 9}$，而由于噪声的存在，矩阵的行向量之间并不是线性相关的。也就是说没有一组解h满足方程组。此时把$Ah=0$称为超定方程。

在此情况下，我们无法得到一个精确的解，只能得到一个近似解。为使得这个近似解尽量准确，我们需要建立一个评判标准，称为损失函数$l$。通过最小化损失函数，我们就可以解出一个近似解。

代数误差

很自然的一个想法是去最小化$||Ah||$，其中$||h||=1$，这其实和求解最小化$||Ah||/||h||$是等价的。

这个问题的求解可以通过求解$A^TA$的最小特征值对应的特征向量得到，也可以通过进行SVD分解$A=U\Sigma V^T$找到最小奇异值对应的V中的奇异向量得到。具体证明在此不多说。

至此我们在有噪声情况下得到了多于4点匹配情况的解。使用这种损失函数称为代数距离。向量$\epsilon = Ah$称为残差向量。

假设一对对应点$x_i,x_i'$，它们通过A贡献点残差向量为$\epsilon_i$，称为代数误差向量。它的范数称为代数距离，即

$$d_{alg}(x'_i,Hx_i)^2=||\epsilon_i||^2=\left\Vert \begin{bmatrix}0^T&-w'_i\textbf x^T_i&y_i'\textbf x^T_i \\\ w'_i\textbf x_i^T&0^T&-x'_i\textbf x_i^T \end{bmatrix}\textbf h\right\Vert^2$$

更一般的，对任意两个向量$x_1,x_2$我们可以写作

$d_{alg}(x_1,x_2)^2=a_1^2+a_2^2$, 其中$a=(a_1,a_2,a_3)^T=x_1\times x_2$

给定一系列的对应点，$\epsilon=Ah$是代数误差向量，可以看到

$$||\epsilon||^2=\sum_i d_{alg}(x'_i,Hx_i)^2=\sum_i||\epsilon_i||^2=||Ah||^2$$

代数距离的优点是形式简单，易于计算，可以看到最小化代数误差基本可以看成DLT算法的一个延伸。它的缺点是没有几何学和统计学的意义，某些情况下不能得到最好的效果。因此它可以用作初值求解。

几何误差

对图像中观测到的一点我们记为x,它的实际坐标我们记为$\bar x$,通过各种方法估计到的该点的位置我们记为$\widehat{x}$。

单幅图像中有噪声时的误差

假设原图像测量很精确，即$x=\bar x$，则此时的几何误差就是转换后图像中的对应点$x'$的测量值与它的理论值Hx的欧氏距离。将两点$x,y$之间的欧氏距离记作$d(x,y)$。则所有匹配点对点误差为

$$\sum_i d(x'_i,H\bar x_i)^2$$

对称转移误差

由于两幅图像中的测量点$x,x'$都有误差，假设变换为H，它的逆变换为$H^{-1}$。则此时的几何误差就是

$$\sum_i d(x_i,H^{-1}x_i')^2+d(x'_i,Hx_i)^2$$
第一项为第一幅图中的转移误差，第二项为第二幅图中的转移误差。显然估计出的单应矩阵 $\widehat H$为使得误差最小时 $H$的值

重投影误差

上面可以看到，无论是将x投影到x’还是将x’投影回x得到的投影点均不与观测值重合。由此我们希望通过寻找一对点$\widehat x,\widehat x'$，以及一个单应矩阵 $\widehat H$，使得两点 $\widehat x,\widehat x'$能够完美的相互投影。

为找到这几个量我们需要最小化的误差函数为

$$\sum_id(x_i,\widehat x_i)^2+d(x'_i,\widehat x'_i)^2,其中\widehat x_i'=\widehat H\widehat x_i$$

重投影误差的几何诠释

点对$\textbf x_i,\textbf x'$的非齐次坐标可以构成$IR^4$中的一点$(x_i,y_i,x'_i,y_i')$。对一个给定的H，对应的两点满足$\textbf x\times (H\textbf x')=0$，从而定义了一个$IR^4$上的代数簇$\mathcal V_H$。（代数簇就是若干多元多项式方程定义的公共零点集。）

$\mathcal V_H$是两个二次超曲面的求交。(因为$\textbf x'\times (H\textbf x)=0$的每一行都是关于坐标$x,x',y,y'$的二次多项式，H定义了多项式的系数，故每一行代表了一个二次超平面；又因为有一行可以被其他两行线性表出，故只有两个）。

给定一个$IR^4$上的点$X_i=(x_i,y_i,x_i',y_i')^T$。估计一个单应矩阵就是找到一个通过点的$\mathcal V_H$。$\mathcal V_H$与H对应。令$\widehat X_i=(\widehat x_i,\widehat y_i,\widehat x_i',\widehat y_i')^T$为$\mathcal V_H$上离$X_i$最近点一个点，有

$$||X_i-\widehat X_i||=d(x_i,\widehat x_i)^2+d(x_i',\widehat x_i')^2$$
这恰好就是重投影误差的公式。找到 $\mathcal V_H$和其上的 $\widehat X_i$等价于找到估计的单应矩阵和估计的匹配点对。

进一步的，$\mathcal V_H$上距离$X$最近点就是使得直线$X\widehat X$垂直于$\mathcal V_H$的点。即

$$d(x_i,\widehat x_i)^2+d(x'_i,\widehat x'_i)^2=d_\perp(X_i,\mathcal V_H)^2$$

总结一下，重投影误差即在$IR^4$估计一个代数簇$\mathcal V_H$使得它与所有的点的距离最小。

Sampson误差

接着上一节，投影误差虽然精确，但是很复杂。我们需要同时估计单应矩阵和对应点，换言之我们需要同时估计$IR^4$上的$\mathcal V_H$和$X_i$，由于前者的非线性，后者的计算估计往往需要需要通过迭代来进行估计，当测量点$X_i$很多时这项任务将变得很困难。因此我们想到假设损失函数在所估计的点附近是线性的，从而近似的估计$X_i$。近似后的误差称为Sampson误差。

$\mathcal V_H$上的点满足$Ah=0$记$\mathcal C_H(X)=0$。做泰勒展开有

$$\mathcal C_H(X+\delta_X)=\mathcal C_H(X)+\frac{\partial \mathcal C_H}{\partial X}\delta_x$$

假设$\delta_x=\widehat X-X$,$\widehat X$在$\mathcal V_H$上等价于$\mathcal C_H(X)+\frac{\partial \mathcal C_H}{\partial X}\delta_x=0$。令$\epsilon=\mathcal C_H(X)$（这也是上一节对代数误差的定义）。有$J\delta_x=-\epsilon$。J为$\mathcal C_H(X)$的Jacobian矩阵。此时我们需要解决的问题就是找到满足该方程的最小的$\delta_x$，用规范的语言表述如下

找到向量$\delta_x$使得$||\delta_x||$取最小值，其中$J\delta_x=-\epsilon$

对该问题我们使用拉格朗日乘子法进行求解，需要找到$\delta_x^T\delta_x-2\lambda(J\delta_x+\epsilon)$的最小值。对$\delta_x,\lambda$分别求偏导解出

$$\delta_x = -J^T(JJ^T)^{-1}\epsilon$$
$$||\delta_x||^2=\delta_x^T\delta_x=\epsilon^T(JJ^T)^{-1}\epsilon$$

此时我们已经消除了$\hat X$的影响。重投影误差可以表示为：

$$\mathcal D_\perp=\sum_i\epsilon_i^T(J_iJ_i^T)^{-1}\epsilon_i$$
不管是J还是ϵ都只与H中的元素相关，而与 $\hat X$无关，重投影误差的优化可以很容易的转化成一个只与H中元素有关的最小二乘问题。通过迭代进行优化。

极大似然估计

极大似然估计在只有一张图像有噪声时等价于转移误差，在两张图像都有噪声时等价于重投影误差。给出了以上两种误差的概率学的诠释。具体推导过程不多赘述。