Structure Computation

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

12.1 问题概述

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

我们先回忆一下第10章三维重建的知识。我们介绍了好几种不同种类的三维重建(取决于我们到底对摄像机矩阵了解多少)。那么结合本章的三角化，我们希望三角化在不同种类的重建之间能给出同样的结果。所以我们首先用τ\tau来代表三角化的过程，我们说如果τ\tau能满足下式，那么我们就说三角化在变换HH下是不变的：

τ(x,x′,P,P′)=H−1τ(x,x′,PH−1,P′H−1)\tau(x,x^{'},P,P^{'}) = H^{-1}\tau(x,x^{'},P H^{-1},P^{'} H^{-1})

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

12.2 线性三角化方法

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

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

Inhomogeneous method 参见4.1.2节，原书P90

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

Inhomogeneous method适合affine reconstruction Homogeneous method不适合affine reconstruction

12.3 几何损失函数

由于图像中有噪声的存在，x↔x′x \leftrightarrow x^{'} 其实不能满足基本矩阵的性质，我们用xˉ,x′ˉ\bar{x},\bar{x'}表示没有噪声的点。那么我们可以构建以下优化函数

C(x,x′)=d(x,x^)2+d(x′,x^′)2x′^Fx^=0C(x,x') = d(x,\hat{x})^2 + d(x',\hat{x}')^2 \\ \hat{x'}F\hat{x} = 0

Sampson approximation

我们定义XX与X^\hat{X}之间的插差为δX\delta_X

δX=−JT(JJT)−1ϵ\delta_X = -J^T(JJ^T)^{-1} \epsilon

其中

ϵ=x′TFx\epsilon = x'^{T}Fx

J=∂ϵ/∂xJ = \partial \epsilon/ \partial x

所以我们可以看出该差值其实是基本矩阵方程关于xx的导数 那么XX和X^\hat{X}之间的关系可以写成

X^=X+δX\hat{X} = X + \delta_X

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

12. 5 三角化的全局最优解

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

12.5.1 问题梳理

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

所以我们就可以构建出以下损失函数

d(x,l)2+d(x′+l′)2d(x,l)^2 + d(x'+l')^2

我们的策略如下:

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

12.5.2 一些需要注意的细节

首先, 两幅图中对应点都不能与极点重合 第二, 我们可以对两幅图都做一个刚体变换,那么x,x′x,x'就可以被放置在原点(0,0,1)(0,0,1), 那么两幅图的极点分别是(1,0,f),(1,0,f′)(1,0,f),(1,0,f'). 我们知道极点也是要满足FF的,所以我们有F(1,0,f)T=(1,0,f′)F=0F(1,0,f)^T = (1,0,f')F = 0, 如此以来我们就可以把基本矩阵表示为一种特殊形式

F=[ff′d−f′c−f′d−fbab−fdcd]F = \left[ \begin{matrix} ff'd & -f'c & -f'd \\ -fb & a & b\\ -fd & c & d \\ \end{matrix} \right]

同时我们也知道极线会通过极点(1,0,f)(1,0,f),我们再找一个特殊点, 那就是极线与yy轴的交点 (0,t,1)(0,t,1), 所以极线就可以写成(1,0,f)×(0,t,1)=(tf,1,−t)(1,0,f) \times (0,t,1) = (tf,1,-t), 那么该直线到原点的距离就是

d(x,l(t))2=t21+(tf)2d(x,l(t))^2 = \frac{t^2}{1+(tf)^2}

紧接着我们找下一个极线:

l′(t)=F(0,t,1)T=(−f′(ct+d),at+b,ct+d)Tl'(t) = F(0,t,1)T=(-f'(ct+d),at+b,ct+d)^T

该极线到原点的距离

d(x′,l′(t))2=(ct+d)2(at+v)2+f′2(ct+d)2d(x',l'(t))^2 = \frac{(ct+d)^2}{(at+v)^2 +f'^2(ct+d)^2}

于是我们把 d(x′,l′(t))2,d(x,l(t))2d(x',l'(t))^2, d(x,l(t))^2 加在一起,记为 s(t)s(t)求导数, 令导数等于0,就可以了.

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

下面我们把整个算法流程重复一遍, 对应于p318 算法12.1

算法输入: 观测到的对应点x↔x′x \leftrightarrow x', 基本矩阵FF

算法输出:寻找一对x^↔x^′\hat{x} \leftrightarrow \hat{x}'可以使几何损失函数最小. 同时这一对点满足 x^′TFx^=0\hat{x}'^{T}F\hat{x} = 0

算法步骤:

  1. 定义一对转换矩阵,可以把x=(x,y,1)T,x′=(x′,y′,t)Tx=(x,y,1)^{T},x'=(x',y',t)^{T}转换到原点

    T=[1−x1−y1]T=\left[ \begin{matrix} 1 & & -x \\ &1 & -y \\ & & 1\\ \end{matrix} \right]

    T'的形式与T是类似的

  2. 将基本矩阵FF变成T′−TFT−1T'^{-T}FT^{-1}

  3. 计算左极点e=(e1,e2,e3)e=(e_1,e_2,e_3)和右极点e′=(e1′,e2′,e3′)e'=(e'_1,e'_2,e'_3), 并且归一化,使得e1+e2=1e_1+e_2=1

  4. 构造两个旋转矩阵, 这两个矩阵可以把ee旋转到(1,0,e3)(1,0,e_3) (1,0,e3′)(1,0,e'_3).

R=[e1e2−e2e11]R=\left[ \begin{matrix} e_1 &e_2 & \\ -e_2 &e_1 & \\ & & 1\\ \end{matrix} \right]

R′R'与RR类似

  1. 把FF改成R′FRTR'FR^{T}
  2. 设置以下等式f=e3,f′=e3,a=F22,b=F23,c=F32,d=F33f=e_3,f'=e_3,a=F_{22},b=F_{23},c=F_{32},d=F_{33}
  3. 将第6步中的等式带入s(t)s(t)中, 求解t
  4. 对求得的解进行验证,同时检查t→∞t \rightarrow \infty 的情况
  5. 将tt带入极线方程,找到x^,x^′\hat{x},\hat{x}', (极线知道了, 观测点x,x′x,x'也知道,求直线上某个点,它要满足到已知点距离最近, 由于我们把x,x′x,x'转到了原点, 那么问题就转变成了直线上求某一点, 它到原点距离最近. 书中给出了一个公式,对于一个一般的直线 (λ,μ,ν)(\lambda, \mu, \nu), 直线上到原点最近的点是(−λν,−μν,λ2+μ2)(-\lambda \nu, -\mu \nu, \lambda^2+\mu^2)
  6. 知道x^,x^′\hat{x},\hat{x}^{'}后,再把他们旋转到原坐标, x^=T−1RTx^\hat{x} = T^{-1} R^{T} \hat{x} x^′=T−1RTx^′\hat{x}^{'} = T^{-1} R^{T} \hat{x}^{'}
  7. 可以顺便利用x^,x^′\hat{x},\hat{x}^{'}计算出三维空间点X^\hat{X}.(三角化, 12.2章)

12.5.3 局部最优解

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

12.5.4 在真实图像上实验

本节大概展示了一些实验结果,在p320

12.6 估计三维点的概率分布

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

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

p(X∣x,x′)=p(x,x′∣X)p(X)/p(x,x′)p(X|x,x') = p(x,x'|X)p(X) / p(x,x')

再加上x,x′x,x'独立的假设, 上式就可以化成

p(X∣x,x′)∼p(x∣X)p(x′∣X)p(X|x,x') \sim p(x|X)p(x'|X)

12.6 线段重建

我们现在要重建空间中的一个线段. 它在两幅图像上分别表示为l,l′l, l'. 我们可以把l,l′l,l'反投影回去, 那么他们在空间中就是两个平面π,,π′\pi,, \pi', 这两个平面的交点就是所求直线. 我们可以形式化的表示为π=PTl,π′=P′Tl′\pi = P^Tl, \pi' = P'^T l', 那么三维空间中的线就可以用这两个平面来表示 (LL是一个2×42 \times 4的矩阵)

L=[lTPl′TP′]L = \left[ \begin{matrix} l^T P \\ l'^T P' \end{matrix} \right]

空间中的点XX在LL上, 所以LX=0LX=0

退化的情况

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

多平面相交的重建

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