12 Structure Computation

本章讲述如何在已知基本矩阵$F$和两幅图像中若干对对应点$x\leftrightarrow x^{\prime }$的情况下计算三维空间点$X$的位置。

12.1 Problem statement

我们假设已知摄像机矩阵$P$和$P^{\prime }$，基本矩阵$F$，还有两幅图像中若干对对应点$x\leftrightarrow x^{\prime }$。因为有噪声的存在，图像中的点反投影回去的两条射线不一定相交，$xF{x}^{\prime }$也不一定等于0，所以简单三角化不一定可行。

我们先回忆一下第10章三维重建的知识。我们介绍了好几种不同种类的三维重建，这取决于我们对摄像机矩阵的知晓程度。那么结合本章的三角化，我们希望三角化在不同种类的重建之间能给出同样的结果。我们首先用$\tau$来代表三角化的过程，如果$\tau$能满足下式，那么我们就说三角化在变换$H$下是不变的：
$$\tau (x,x^{\prime },P,P^{\prime })=H^{-1}\tau (x,x^{\prime },P{H}^{-1},P^{\prime }{H}^{-1})$$

为什么需要讨论这个？这是因为我们首先需要确定三维重建的种类，才能决定优化目标的形式。如果我们只知道摄像机矩阵是一个projective matrix，那么我们就不能在三维空间最优化目标函数。因为这样的优化函数在投影变换中不能给出唯一的结果，因为距离和垂直度等概念在projective geometry的背景下无效。所以，本章给出的三角化方法优化的是二维图像上的距离，所以本章的方法在投影变换（projective transformation）中是不变的。

12.2 Linear triangulation methods

对于两幅图像，我们分别有$x=PX,x^{\prime }=P{X}^{\prime }$，我们可以将第一个方程改成$x\times PX=0$，第二幅图也一样。我们继续改写就可以有$AX=0$。

Homogeneous method 找出$A$最小特征值对应的特征向量

Inhomogeneous method 参见4.1.2节，原书P90

讨论
Inhomogeneous method假设点不在无穷远处，不适合projective reconstruction。其实这两个方法都不适合。

Inhomogeneous method适合affine reconstruction。

Homogeneous method不适合affine reconstruction。

12.3 Geometric error cost function

由于图像中有噪声的存在，$x\leftrightarrow x^{\prime }$其实不能满足极线的约束，我们用$\bar{x},\overline{x^{\prime }}$表示没有噪声的点。那么我们可以构建以下优化函数：

$$\begin{gathered} C(x,x^{\prime })=d(x,\hat{x}{)}^2+d(x^{\prime },{\hat{x}}^{\prime }{)}^2 \\ subjectto{\widehat{x^{\prime }}}^{T}F\hat{x}=0 \end{gathered}$$

其中$d$表示两点之间的欧氏距离。这相当于最小化点$X$的重投影误差，该点$X$通过与$F$一致的投影矩阵映射到两个点，如图12.2。

12.4 Sampson approximation (first-order geometric correction)

我们定义$X$与$\hat{X}$之间的差为${\delta }_X$：
$${\delta }_X=-J^T(J{J}^T{)}^{-1}ϵ$$

其中
$$\begin{gathered} ϵ=x^{\prime T}Fx \\ J=\partial ϵ/\partial x=[(F^T{x}^{\prime }{)}_1,(F^T{x}^{\prime }{)}_2,(FX{)}_1,(FX{)}_2] \end{gathered}$$

其中$(F^T{x}^{\prime }{)}_1=f_{11}{x}^{\prime }+f_{21}{y}^{\prime }+f_{31}$，以此类推。
所以我们可以看出该差值其实是基本矩阵方程关于$x$的导数
那么$X$和$\hat{X}$之间的关系可以写成：
$$\hat{X}=X+{\delta }_X$$

我们只需要把${\delta }_X$算出来，然后对计算出的理论点$X$按照上式进行一个纠正就可以了。

12.5 An optimal solution

本节介绍一种可以找到全局最优解的优化函数，并且是非迭代的，我们同时假设噪声服从高斯分布。

12.5.1 Reformulation of the minimization problem

先对问题进行一个梳理。

我们知道第一幅图的极点一定在极线上，第二幅图的极点也满足这个性质。反过来，在极线上的点也满足基本矩阵的约束。那么就能让观测到的点尽可能靠近极线，也就是找观测点到极线的距离，并使其最小。

所以我们就可以构建出以下损失函数
$$d(x,l{)}^2+d(x^{\prime }+l^{\prime }{)}^2$$

我们的策略如下:

1. 将极线方程参数化，所以第一幅图像中的极线方程就可以写为$l(t)$
2. 利用基本矩阵$F$,和$l(t)$来计算第二幅图像中的极线l${}^{\prime }(t)$
3. 将损失函数写成$d(x,l(t){)}^2+d(x^{\prime }+l^{\prime }(t){)}^2$
4. 求解最优的$t$

12.5.2 Details of the minimization

接下来我们讲一下需要注意的一些细节。

首先，两幅图中对应点都不能与极点重合。

并且，我们可以对两幅图都做一个刚体变换，那么$x,x^{\prime }$就可以被放置在原点$(0,0,1)$，那么两幅图的极点分别是$(1,0,f),(1,0,f^{\prime })$。我们知道极点也是要满足$F$的，所以我们有$F(1,0,f{)}^T=(1,0,f^{\prime })F=0$，如此以来我们就可以把基本矩阵表示为一种特殊形式：
$$F=\left[\begin{array}{ccc}f{f}^{\prime }d& -f^{\prime }c& -f^{\prime }d\\ -fb& a& b\\ -fd& c& d\end{array}\right]$$

同时我们也知道极线会通过极点$(1,0,f)$，我们再找一个特殊点，那就是极线与$y$轴的交点 $(0,t,1)$，所以极线就可以写成$(1,0,f)\times (0,t,1)=(tf,1,-t)$，那么该直线到原点的距离就是：
$$d(x,l(t){)}^2=\frac{t^2}{1+(tf{)}^2}$$

紧接着我们找下一个极线：
$$l^{\prime }(t)=F(0,t,1)T=(-f^{\prime }(ct+d),at+b,ct+d{)}^T$$

该极线到原点的距离：
$$d(x^{\prime },l^{\prime }(t){)}^2=\frac{(ct+d{)}^2}{(at+v{)}^2+f^{\prime 2}(ct+d{)}^2}$$

于是我们把 $d(x^{\prime },l^{\prime }(t){)}^2,d(x,l(t){)}^2$ 加在一起，记为 $s(t)$求导数，令导数等于0，就可以了。

一些讨论 $s(t)$是6次多项式，那么它就有6个实根，对应于3个最小值和3个最大值。顺便别忘了检查$x\to \infty$的情况。

下面我们把整个算法流程重复一遍，对应于P318算法12.1。

算法输入：观测到的对应点$x\leftrightarrow x^{\prime }$，基本矩阵$F$

算法输出：寻找一对$\hat{x}\leftrightarrow {\hat{x}}^{\prime }$可以使几何损失函数最小，同时这一对点满足 ${\hat{x}}^{\prime T}F\hat{x}=0$

算法步骤：

1. 定义一对转换矩阵，可以把$x=(x,y,1{)}^T,x^{\prime }=(x^{\prime },y^{\prime },t{)}^T$转换到原点
   $$T=\left[\begin{array}{ccc}1& & -x\\ & 1& -y\\ & & 1\end{array}\right]$$

   $T^{\prime }$的形式与$T$是类似的

2. 将基本矩阵$F$变成$T^{\prime -T}F{T}^{-1}$

3. 计算左极点$e=(e_1,e_2,e_3)$和右极点$e^{\prime }=(e_1^{\prime },e_2^{\prime },e_3^{\prime })$，并且归一化，使得$e_1+e_2=1$

4. 构造两个旋转矩阵，这两个矩阵可以把$e$旋转到$(1,0,e_3)$ $(1,0,e_3^{\prime })$.
   $$R=\left[\begin{array}{ccc}{e}_1& e_2& \\ -e_2& e_1& \\ & & 1\end{array}\right]$$
   $R^{\prime }$与$R$类似

5. 把$F$改成$R^{\prime }F{R}^T$

6. 设置以下等式$f=e_3,f^{\prime }=e_3,a=F_{22},b=F_{23},c=F_{32},d=F_{33}$

7. 将第6步中的等式带入$s(t)$中，求解t

8. 对求得的解进行验证，同时检查$t\to \infty$ 的情况

9. 将$t$带入极线方程，找到$\hat{x}，{\hat{x}}^{\prime }$，极线知道了，观测点$x,x^{\prime }$也知道，求直线上某个点，它要满足到已知点距离最近，由于我们把$x,x^{\prime }$转到了原点，那么问题就转变成了直线上求某一点，它到原点距离最近。书中给出了一个公式，对于一个一般的直线$(\lambda ,\mu ,\nu )$，直线上到原点最近的点是$(-\lambda \nu ,-\mu \nu ,{\lambda }^2+{\mu }^2)$

10. 知道$\hat{x},{\hat{x}}^{\prime }$后，再把他们旋转到原坐标，$\hat{x}=T^{-1}{R}^T\hat{x}$ ${\hat{x}}^{\prime }=T^{-1}{R}^T{\hat{x}}^{\prime }$

11. 可以顺便利用$\hat{x},{\hat{x}}^{\prime }$计算出三维空间点$\hat{X}$（三角化，12.2）

12.5.3 Local minima

$g(t)$有6个自由度，所以它最多有三个最小值。那么如果用迭代的方法去寻找最小值，可能陷在局部最小值里出不来。

12.5.4 Evaluation on real images

本节大概展示了一些实验结果，在P320

12.6 Probability distribution of the estimated 3D point

估计三维点的概率分布。

通过两幅图像估计出来的三维空间点应该是满足一定概率分布的。其准确与否主要取决于从摄像机出发的，两条射线之间的角度。本节就对这个问题进行建模。书中为了简化这个问题，只考虑空间某平面上的点$X=(x,y{)}^T$，其图像上的点分别表示为$x=f(X),x^{\prime }=f^{\prime }(X)$, $f,f^{\prime }$是$2\times 3$的矩阵，而不是$3\times 4$ 如果忘了可以复习一下p175 6.4.2节

我们线考虑第一幅图像上的点$x$，并且我们假设噪声服从均值为0，方差为${\sigma }^2$的高斯分布，那么在已知$X$的条件下$x$的概率分布可以表示为$p(x|X)$，对第二幅图上的点$x^{\prime }$有相同的结论$p(x^{\prime }|X)$。那么当$x,x^{\prime }$已知的时候，我们可以用贝叶斯公式反推$X$的概率分布

$$p(X|x,x^{\prime })=p(x,x^{\prime }|X)p(X)/p(x,x^{\prime })$$

再加上$x,x^{\prime }$独立的假设，上式就可以化成

$$p(X|x,x^{\prime })\sim p(x|X)p(x^{\prime }|X)$$

12.6 线段重建

我们现在要重建空间中的一个线段。它在两幅图像上分别表示为$l,l^{\prime }$。我们可以把$l,l^{\prime }$反投影回去，那么他们在空间中就是两个平面$\pi ,{\pi }^{\prime }$, 这两个平面的交点就是所求直线。我们可以形式化的表示为$\pi =P^{T}l,{\pi }^{\prime }=P^{\prime T}{l}^{\prime }$，那么三维空间中的线就可以用这两个平面来表示 ($L$是一个$2\times 4$的矩阵)
$$L=\left[\begin{array}{c}{l}^{T}P\\ l^{\prime T}{P}^{\prime }\end{array}\right]$$

空间中的点$X$在$L$上，所以$LX=0$

退化的情况

如果这个直线在极平面上，那么上一节的方法就失效了，而且这样直线会和基线相交。在实际情况下，几乎要和基线相交的线也不能用以上方法来重建.

多平面相交的重建

假设有$n$个平面，那么我们就他们像前文$L$一样放在一起，形成一个$n\times 4$的矩阵$A$。对$A$做SDV分解 $A=UD{V}^T$，从$D$中找出两个最大的特征值对应的特征向量，用他们来表示平面，也可以假设空间中直线$L$投影到各个平面，然后计算投影直线和观测直线之间的几何损失函数，用极大似然估计求解。