结构计算
本章介绍 3 D 3D 3D点在两幅视图上的像和这些视图的摄像机矩阵给定时,如何计算该点的位置。因为被测量点 x , x ′ x,x^\prime x,x′都有误差,所以从点反向投影额射线不共面。这意味着不存在准确的满足 x = P X , x ′ = P ′ X x=PX,x^\prime = P^\prime X x=PX,x′=P′X的点 X X X,并且图像点不满足对极几何约束 x ′ T F x = 0 x^{\prime T}F x = 0 x′TFx=0。
三角测量法的一个理想特性是其在重构的某种适当的变换下保持不变。三角测量法的关键思想是估计一个精确满足给定的摄像机几何的 3 D 3D 3D点 X ^ \hat X X^,因此它投影为 x ^ = P X ^ , x ^ ′ = P ′ X ^ \hat x = P\hat X,\hat x^\prime = P^\prime \hat X x^=PX^,x^′=P′X^且目标是由图像测量 x , x ′ x,x^\prime x,x′来估计 X X X
- 线性三角测量法
由于 x = P X , x ′ = P ′ X x=PX,x^\prime = P^\prime X x=PX,x′=P′X,构建 x × ( P X ) = 0 x \times (PX)=0 x×(PX)=0构建 A X = 0 AX=0 AX=0的方程求解。但是这种方法不是射影不变的 - 几何误差代价函数
寻求最小化函数
C ( x , x ′ ) = d ( x , x ^ ) 2 + d ( x ′ , x ^ ′ ) 2 , x ^ ′ F x ^ = 0 \mathcal C(x,x^\prime) = d(x,\hat x)^2 + d(x^\prime,\hat x ^\prime)^2,\hat x^\prime F \hat x = 0 C(x,x′)=d(x,x^)2+d(x′,x^′)2,x^′Fx^=0 - Sampson近似
在这里我们关心的是对测量点进行矫正计算。测量点 X = ( x , y , x ′ , y ′ ) X=(x,y,x^\prime,y^\prime) X=(x,y,x′,y′)的 S a m p s o n Sampson Sampson矫正 δ X \delta_X δX是
δ X = − J T ( J J T ) − 1 ε \delta_X = -J^T(JJ^T)^{-1} \pmb \varepsilon δX=−JT(JJT)−1ε
矫正的点为
X ^ = X + δ X = X − J T ( J J T ) − 1 ε \hat X = X+\delta_X = X-J^T(JJ^T)^{-1} \pmb \varepsilon X^=X+δX=X−JT(JJT)−1ε
前面已经证明,对 x ′ T F x = 0 x^{\prime T}Fx=0 x′TFx=0所定义的族,该误差是 ε = x ′ T F x \pmb \varepsilon =x^{\prime T}Fx ε=x′TFx,其 J a c o b i a n Jacobian Jacobian是
J = ∂ ε / ∂ x = [ ( F T x ′ ) 1 , ( F T x ′ ) 2 , ( F x ) 1 , ( F x ) 2 ] J = \partial \pmb \varepsilon / \partial x = [(F^Tx^\prime)_1,(F^Tx^\prime)_2,(Fx)_1,(Fx)_2] J=∂ε/∂x=[(FTx′)1,(FTx′)2,(Fx)1,(Fx)2]
矫正点的一阶近似为
[ x ^ y ^ x ^ ′ y ^ ′ ] = [ x y x ′ y ′ ] − x ′ T F x ( F T x ′ ) 1 2 2 + ( F T x ′ ) 2 2 2 + ( F x ) 1 2 2 + ( F x ) 2 2 [ ( F T x ′ ) 1 ( F T x ′ ) 2 ( F x ) 1 ( F x ) 2 ] \begin{bmatrix} \hat x \\ \hat y \\ \hat x^\prime \\ \hat y^\prime \end{bmatrix} = \begin{bmatrix} x \\ y \\ x^\prime \\ y^\prime \end{bmatrix} - \frac{x^{\prime T} Fx}{(F^Tx^\prime)_12^2+(F^Tx^\prime)_22^2+(Fx)_12^2+(Fx)_2^2} \begin{bmatrix} (F^Tx^\prime)_1 \\ (F^Tx^\prime)_2 \\ (Fx)_1 \\ (Fx)_2 \end{bmatrix} x^y^x^′y^′ = xyx′y′ −(FTx′)122+(FTx′)222+(Fx)122+(Fx)22x′TFx (FTx′)1(FTx′)2(Fx)1(Fx)2
针对其中最优的几何代价函数,可以通过转换为解一个 6 6 6次多项式的非迭代方法得到
最优解
最小化问题的重新公式化
给定一组测量得到的对应
x
↔
x
′
x \leftrightarrow x^\prime
x↔x′,我们要寻找
x
^
,
x
^
′
\hat x,\hat x^\prime
x^,x^′,它们最小化几何距离平方和并满足对极几何约束
x
^
′
T
F
x
^
=
0
\hat x^{\prime T}F\hat x=0
x^′TFx^=0。任意一对满足对极约束的点必须在两幅图像中对应的对极线上,
x
⊥
,
x
⊥
′
x_\perp,x_\perp^\prime
x⊥,x⊥′分别为在
I
,
I
′
I,I^\prime
I,I′上最接近
x
,
x
′
x,x^\prime
x,x′的点,因此,几何代价函数变为
d
(
x
,
I
)
2
+
d
(
x
′
,
I
′
)
2
d(x,I)^2 + d(x^\prime,I^\prime)^2
d(x,I)2+d(x′,I′)2
最小化的策略如下:
- 用参数 t t t参数化第一幅图像中的对极约束。因此第一幅图像中的对极线可以被记为 I ( t ) I(t) I(t)
- 利用基本矩阵 F F F,计算第二幅图像中对应的对极线 I ′ ( t ) I^\prime(t) I′(t)
- 把距离函数 d ( x , I ( t ) ) 2 + d ( x ′ , I ′ ( t ) ) 2 d(x,I(t))^2+d(x^\prime,I^\prime(t))^2 d(x,I(t))2+d(x′,I′(t))2显示的表示为 t t t的函数
- 求这个函数最小化
最小化细节
假设两个图像点
x
,
x
′
x,x^\prime
x,x′都不是对极点。对每一幅图像都做刚性变换使得两个点
x
,
x
′
x,x^\prime
x,x′都是齐次坐标系的原点
(
0
,
0
,
1
)
T
(0,0,1)^T
(0,0,1)T。此外,将对极点分别置于
x
x
x轴上的点
(
1
,
0
,
f
)
T
(1,0,f)^T
(1,0,f)T和
(
1
,
0
,
f
′
)
T
(1,0,f^\prime)^T
(1,0,f′)T。由于
F
(
1
,
0
,
f
)
T
=
(
1
,
0
,
f
′
)
T
F
=
0
F(1,0,f)^T = (1,0,f^\prime)^T F = 0
F(1,0,f)T=(1,0,f′)TF=0,该基本矩阵有以下特殊形式
F
=
[
f
f
′
d
−
f
′
c
−
f
′
d
−
f
b
a
b
−
f
d
c
d
]
F = \begin{bmatrix} ff^\prime d & -f^\prime c & -f^\prime d \\ -fb & a & b \\\ -fd & c & d \end{bmatrix}
F=
ff′d−fb −fd−f′cac−f′dbd
考虑第一幅图像中过点
(
0
,
t
,
1
)
T
(0,t,1)^T
(0,t,1)T和对极点
(
1
,
0
,
f
)
T
(1,0,f)^T
(1,0,f)T的对极线,记这条极线为
I
(
t
)
I(t)
I(t),这个对极线的向量表示为
(
0
,
t
,
1
)
T
×
(
1
,
0
,
f
)
T
=
(
t
f
,
1
,
−
t
)
T
(0,t,1)^T \times (1,0,f)^T = (tf,1,-t)^T
(0,t,1)T×(1,0,f)T=(tf,1,−t)T给出,所以由原点到该直线的距离的平方是
d
(
x
,
I
(
t
)
)
2
=
t
2
1
+
(
t
f
)
2
d(x,I(t))^2 = \frac{t^2}{1+(tf)^2}
d(x,I(t))2=1+(tf)2t2
利用基本矩阵计算另一幅图像的对极线,得到
I
′
(
t
)
=
F
(
0
,
t
,
1
)
T
=
(
−
f
′
(
c
t
+
d
)
,
a
t
+
b
,
c
t
+
d
)
T
I^\prime(t) = F(0,t,1)^T = (-f^\prime(ct+d),at+b,ct+d)^T
I′(t)=F(0,t,1)T=(−f′(ct+d),at+b,ct+d)T
从原点到该直线的距离的平方为
d
(
x
′
,
I
′
(
t
)
)
2
=
(
c
t
+
d
)
2
(
a
t
+
b
)
2
+
f
′
2
(
c
t
+
d
)
2
d(x^\prime,I^\prime(t))^2 = \frac{(ct+d)^2}{(at+b)^2 +f^{\prime 2}(ct+d)^2}
d(x′,I′(t))2=(at+b)2+f′2(ct+d)2(ct+d)2
因此总的距离的平方是
s
(
t
)
=
t
2
1
+
(
t
f
)
2
+
(
c
t
+
d
)
2
(
a
t
+
b
)
2
+
f
′
2
(
c
t
+
d
)
2
s(t) = \frac{t^2}{1+(tf)^2} + \frac{(ct+d)^2}{(at+b)^2 +f^{\prime 2}(ct+d)^2}
s(t)=1+(tf)2t2+(at+b)2+f′2(ct+d)2(ct+d)2
我们的任务就是求这个式子的最小值,求其导数,
s
(
t
)
s(t)
s(t)的最大值和最小值将在
s
′
(
t
)
=
0
s^\prime(t)=0
s′(t)=0处出现。把
s
′
(
t
)
s^\prime(t)
s′(t)中的两项在一个公分母上合并并且使分子等于
0
0
0,给出条件
g
(
t
)
=
t
(
(
a
t
+
b
)
2
+
f
′
2
(
c
t
+
d
)
2
)
2
−
(
a
d
−
b
c
)
(
1
+
f
2
t
2
)
2
(
a
t
+
b
)
(
c
t
+
d
)
=
0
g(t) = t((at+b)^2 + f^{\prime 2}(ct+d)^2)^2 - (ad-bc)(1+f^2t^2)^2 (at+b)(ct+d) = 0
g(t)=t((at+b)2+f′2(ct+d)2)2−(ad−bc)(1+f2t2)2(at+b)(ct+d)=0
s
(
t
)
s(t)
s(t)的最大值和最小值将在这个多项式的根处出现,这是一个
6
6
6次多项式,它最多有
6
6
6个实根,对应于
s
(
t
)
s(t)
s(t)的
3
3
3个极大值和
3
3
3个最小值
算法流程
目标
已知一组测量得到的点对应 x ↔ x ′ x \leftrightarrow x^\prime x↔x′和基本矩阵 F F F,计算在对极几何约束 x ^ ′ F x ^ = 0 \hat x^\prime F \hat x = 0 x^′Fx^=0下最小化几何误差的矫正对应 x ^ ↔ x ^ ′ \hat x \leftrightarrow \hat x^\prime x^↔x^′
算法
-
定义变换矩阵
T = [ 1 − x 1 − y 1 ] and T ′ = [ 1 − x ′ 1 − y ′ 1 ] T = \begin{bmatrix} 1 & & -x \\ & 1 & -y \\ & &1 \end{bmatrix} \text{ and } T^\prime = \begin{bmatrix} 1 & & -x^\prime \\ & 1 & -y^\prime \\ & &1 \end{bmatrix} T= 11−x−y1 and T′= 11−x′−y′1
这些变换是将 x = ( x , y , 1 ) T \pmb x=(x,y,1)^T x=(x,y,1)T和 x ′ = ( x ′ , y ′ , 1 ) T \pmb x^\prime=(x^\prime,y^\prime,1)^T x′=(x′,y′,1)T变到坐标原点的平移 -
用 T ′ − 1 F T − 1 T^{\prime -1} F T^{-1} T′−1FT−1代替 F F F。新的 F F F对应于平移后的坐标
-
计算满足 e ′ T F = 0 , F e = 0 e^{\prime T}F=0,Fe=0 e′TF=0,Fe=0的右、左对极点 e = ( e 1 , e 2 , e 3 ) T e=(e_1,e_2,e_3)^T e=(e1,e2,e3)T和 e ′ = ( e 1 ′ , e 2 ′ , e 3 ′ ) T e^\prime=(e_1^\prime,e_2^\prime,e_3^\prime)^T e′=(e1′,e2′,e3′)T。归一化 e e e使得 e 1 2 + e 2 2 + e 3 2 = 1 e_1^2+e_2^2+e_3^2=1 e12+e22+e32=1,并对 e ′ e^\prime e′做同样的处理
-
构造矩阵
R = [ e 1 e 2 − e 2 e 1 1 ] and R ′ = [ e 1 ′ e 2 ′ − e 2 ′ e 1 ′ 1 ] R=\begin{bmatrix} e_1 & e_2 &\\ -e_2 & e_1 & \\ && 1 \end{bmatrix} \text{ and } R^\prime = \begin{bmatrix} e_1^\prime & e_2^\prime &\\ -e_2^\prime & e_1^\prime & \\ && 1 \end{bmatrix} R= e1−e2e2e11 and R′= e1′−e2′e2′e1′1
注意 R , R ′ R,R^\prime R,R′是旋转矩阵,并满足 R e = ( 1 , 0 , e 3 ) T , R ′ e ′ = ( 1 , 0 , e 3 ′ ) T Re=(1,0,e_3)^T,R^\prime e^\prime = (1,0,e_3^\prime)^T Re=(1,0,e3)T,R′e′=(1,0,e3′)T -
用 R ′ F R − 1 R^\prime FR^{-1} R′FR−1代替 F F F
-
令 f = e 3 , f ′ = e 3 ′ , a = F 22 , b = F 23 , c = F 32 , d = F 33 f=e_3,f^\prime =e_3^\prime,a=F_{22},b=F_{23},c=F_{32},d=F_{33} f=e3,f′=e3′,a=F22,b=F23,c=F32,d=F33
-
根据上述分析,构造关于 t t t的多项式 g ( t ) g(t) g(t),得到 t t t的 6 6 6个根
-
计算代价函数在 s ( t ) s(t) s(t)在 g ( t ) g(t) g(t)的每个根的实部的值。同时求当 t = ∞ t= \infty t=∞时几何代价函数的渐近线,选择使代价函数取最小值的 t m i n t_{min} tmin
-
计算两条直线 I = ( t f , 1 , − t ) T I=(tf,1,-t)^T I=(tf,1,−t)T和 I ′ ( t ) I^\prime(t) I′(t)在 t m i n t_{min} tmin处的值,并求出这些直线上最接近于原点的点 x ^ , x ^ ′ \hat x,\hat x^\prime x^,x^′,对于一般的直线 ( λ , μ , v ) T (\lambda , \mu,v)^T (λ,μ,v)T,直线上最接近原点的点是 ( − λ v , − μ v , λ 2 + μ 2 ) T (-\lambda v,-\mu v,\lambda^2+\mu^2)^T (−λv,−μv,λ2+μ2)T
-
用 T − 1 R T x ^ T^{-1}R^T\hat x T−1RTx^替代 x ^ \hat x x^并用 T ′ − 1 R ′ T x ^ ′ T^{\prime -1}R^{\prime T}\hat x^\prime T′−1R′Tx^′替代 x ^ ′ \hat x^\prime x^′变回原来的坐标
-
3 3 3维空间点 X ^ \hat X X^由齐次方法线性三角测量法得到