14 Affine Epipolar Geometry
本章主要是在仿射摄像机的情况下重新考虑对极几何,也就是仿射对极几何。
仿射摄像机的优点是它是线性的,所以很多最优化算法可以用线性代数的知识解决。如果是一般的投影摄像机,很多算法就不是线性的了(比如三角化)。
文章目录
- 14 Affine Epipolar Geometry
- 14.1 Affine epipolar geometry
- 14.2 The affine fundamental matrix
- 14.2.1 Derivation
- 14.2.2 Properties
- 14.3 Estimating F A F_{A} FA from image point correspondences
- 14.3.1 The linear algorithm
- 14.3.2 The Gold Standard algorithm
- 14.4 Triangulation
- 14.5 Affine reconstruction
- 14.6 Necker reversal and the bas-relief ambiguity
- 14.7 Computing the motion
14.1 Affine epipolar geometry
仿射摄像机的特点是它的光心在无穷远,所以从三维空间像二维平面投影时,投影射线时互相平行的。在这种情况下,对极几何就可以得到简化。
极线 所有极线都是互相平行的,因为不同点向图像投影的射线是平行的。
极点 因为所有极线互相平行,那么极点在无穷远处。
14.2 The affine fundamental matrix
仿射基本矩阵 F A F_A FA长这样:
[ 0 0 a 0 0 b c d e ] \left[ \begin{matrix} 0 & 0 & a \\ 0 & 0 & b \\ c & d & e \\ \end{matrix} \right] 00c00dabe
因为 a , b , c , d , e a,b,c,d,e a,b,c,d,e都不为0,那么 F A F_A FA的秩就是3,一般情况下, F F F秩是2。
14.2.1 Derivation
接下来讲一下 F A F_A FA的推导过程。
从几何角度推导 F A F_A FA的过程
- 考虑对应点之间的对应关系。因为从图像1到空间平面,再到图像2的所有变换都是基于平面的仿射变换,所以对应点之间的变换也是仿射变换,也就是 x ′ = H A x x'=H_A x x′=HAx
- 构造极线。极线是通过极点和 x ′ x' x′构造的。所以 l ′ = e ′ × H A x = F A x l'=e' \times H_A x =F_A x l′=e′×HAx=FAx,所以 F A = [ e ′ ] × H A F_A=[e']_{\times} H_A FA=[e′]×HA。我们现在考虑仿射矩阵的特殊形式,以及当 e ′ e' e′为无穷大时的偏斜矩阵 [ e ′ ] × [e']_{\times} [e′]×,因此最后一个元素为零,具体展开如下:
F A = [ e ′ ] × H A = [ 0 0 ∗ 0 0 ∗ ∗ ∗ 0 ] [ ∗ ∗ ∗ ∗ ∗ ∗ 0 0 1 ] = [ 0 0 ∗ 0 0 ∗ ∗ ∗ ∗ ] F_A =[e']_{\times} H_A \\ = \left[ \begin{matrix} 0 & 0 & * \\ 0 & 0 & * \\ * & * & 0 \\ \end{matrix} \right] \left[ \begin{matrix} * & * & * \\ * & * & * \\ 0 & 0 & 1 \\ \end{matrix} \right] \\ = \left[ \begin{matrix} 0 & 0 & * \\ 0 & 0 & * \\ * & * & * \\ \end{matrix} \right] FA=[e′]×HA= 00∗00∗∗∗0 ∗∗0∗∗0∗∗1 = 00∗00∗∗∗∗
从代数角度推导 F A F_A FA的过程
F = [ e ′ ] × P ′ P + F=[e']_{\times}P'P^{+} F=[e′]×P′P+,把仿射摄像机的矩阵带入即可。
14.2.2 Properties
F A F_A FA有5个非零元素,所以有4个自由度,是这么算的:2个极点,每个贡献一个自由度,2个平面上的极线互相映射,贡献2个自由度。
因为 F A e = 0 F_A e = 0 FAe=0,所以极点表示为 ( − d , c , 0 ) (-d,c,0) (−d,c,0),所以 e e e在 l ∞ l_{\infty} l∞上。
点 x x x对应的极线 l ′ = F A x = ( a , b , c x + d y + e ) T l'=F_A x=(a,b,cx+dy+e)^T l′=FAx=(a,b,cx+dy+e)T,此式表明所有的极线都是平行的,因为 ( a , b ) (a,b) (a,b)与 ( x , y ) (x,y) (x,y)互相独立
14.3 Estimating F A F_{A} FA from image point correspondences
给定若干对对应点,矩阵是 x ′ T F A x = 0 x^{'T}F_{A}x=0 x′TFAx=0定义的,所以我们可以从对应点中计算 F A F_A FA。
14.3.1 The linear algorithm
因为 F A F_A FA规定了点与点之间的对应关系,所以给出足够过的点,肯定是能把 F A F_A FA计算出来的。比如我们规定 x i = ( x i , y i , 1 ) , x i ′ = ( x i ′ , y i ′ , 1 ) x_i=(x_i,y_i,1),x'_i=(x'_i,y'_i,1) xi=(xi,yi,1),xi′=(xi′,yi′,1)。这样一对对应点可以构造一个线性方程:
a x i ′ + b y i ′ + c x i + d y i + e = 0 ax'_i+by'_i+cx_i+dy_i+e=0 axi′+byi′+cxi+dyi+e=0
写成矩阵形式就是 A f = 0 Af=0 Af=0, A A A是一个 n × 5 n \times 5 n×5的矩阵,最少需要4对点。 F A F_A FA有4个自由度,剩下的1个变量可由其他4个来表示。
奇异性约束
我们知道一般情况的 F F F是奇异矩阵,所以 F A F_A FA也应该是一个奇异性矩阵。而且 F A F_A FA的形式就确定了它的rank不会大于2,所以没有必要将奇异性约束像求普通 F F F那样加入求解的过程,换言之,不用考虑这个约束,直接求解就可以了。
几何解释
求解两个图像上的对应点就是在四维空间中去拟合一个平面( a x i ′ + b y i ′ + c x i + d y i + e = 0 ax'_i+by'_i+cx_i+dy_i+e=0 axi′+byi′+cxi+dyi+e=0 就是四维空间中的平面)。这样做有两个好处:第一,求 F F F就是一个平面拟合过程,容易思考。第二,可以用sampson损失函数,因为它是唯一的一阶近似方法。
14.3.2 The Gold Standard algorithm
知道了求解 F F F的几何解释,我们就用这个几何解释来求解它。求解的过程就是黄金标准算法。我们考虑有噪声的情况,理论点可以表示为 x i ^ , x i ^ ′ \hat{x_i},\hat{x_i}' xi^,xi^′,实际观测点是 x i , x i ′ x_i,x'_i xi,xi′,所谓的黄金标准算法就是优化以下函数:
m i n ∑ i d ( x i , x i ^ ) 2 + d ( x i ′ , x i ^ ′ ) 2 min \sum_{i} d(x_i,\hat{x_i})^2 + d(x'_i,\hat{x_i}')^2 mini∑d(xi,xi^)2+d(xi′,xi^′)2
其中 x i ^ , x i ^ ′ \hat{x_i},\hat{x_i}' xi^,xi^′满足 x i ^ ′ T F A x i ^ = 0 \hat{x_i}^{'T}F_A \hat{x_i}=0 xi^′TFAxi^=0,如果我们考虑几何解释,我们就可以考虑四维空间的一个点 X i = ( x i ′ , y i ′ , x i , y i ) X_i=(x'_i,y'_i,x_i,y_i) Xi=(xi′,yi′,xi,yi),用这个点去拟合一个平面,其参数为 ( a , b , c , d , e ) (a,b,c,d,e) (a,b,c,d,e),那么就是求点到平面的最小距离。
d ⊥ = a x i ′ , b y i ′ + c x i + d y i + e a 2 + b 2 + c 2 + d 2 d_{\perp} = \frac{ax'_i,by'_i+cx_i+dy_i+e}{\sqrt{a^2+b^2+c^2+d^2}} d⊥=a2+b2+c2+d2axi′,byi′+cxi+dyi+e
求解cost function函数的过程如下:
先对 e e e求导数,令其等于0,可以得到:
e = − 1 n ∑ ( N T X i ) = − N T X ˉ e=-\frac{1}{n} \sum(N^T X_i) = -N^T \bar{X} e=−n1∑(NTXi)=−NTXˉ
N = ( a , b , c , d ) N=(a,b,c,d) N=(a,b,c,d)
N ˉ \bar{N} Nˉ就是所有已知 X X X的均值,也就是质心。
把 e e e反带回 d ⊥ d_{\perp} d⊥ 得到:
d ⊥ = 1 ∣ ∣ N ∣ ∣ 2 ∑ i ( N T Δ X i ) d_{\perp} = \frac{1}{||N||^2} \sum_i(N^T \Delta X_i) d⊥=∣∣N∣∣21i∑(NTΔXi)
Δ X i = X i − X ˉ \Delta X_i = X_i - \bar{X} ΔXi=Xi−Xˉ
用 Δ X i \Delta X_i ΔXi的行构造一个矩阵,然后直接SVD分解就可以了。
最后一个需要注意的点:黄金标准算法需要多于4对对应点。那么如果我们只知道4对对应点,该怎么办? 下面的算法来解决。
4对对应点求解 F A F_A FA 的步骤如下:
- 用前三对对应点计算一个仿射变换 H A H_A HA,也就是求解 x i ′ = H A x i x'_i=H_A x_i xi′=HAxi
- 用 H A H_A HA计算 H A x 4 H_A x_4 HAx4,然后 ( H A x 4 × x 4 ′ ) (H_A x_4 \times x'_4) (HAx4×x4′) 得到极线 l ′ l' l′
那么极点就是 e ′ = ( − l 2 ′ , l 1 ′ , 0 ) e'=(-l'_2,l'_1,0) e′=(−l2′,l1′,0) - 所以对任何一个点 x x x,它对应的极线就是 e ′ × ( H A x ) = F A x e' \times (H_A x) = F_A x e′×(HAx)=FAx
14.4 Triangulation
现在假设我们知道一对对应点 ( x , y ) ↔ ( x ′ , y ′ ) (x,y) \leftrightarrow (x',y') (x,y)↔(x′,y′) 和仿射基本矩阵 F A F_A FA,因为已知对应点是含噪声的,我们想要确定不含噪声的点 ( x ^ , y ^ ) ↔ ( x ^ ′ , y ^ ′ ) (\hat{x},\hat{y}) \leftrightarrow (\hat{x}',\hat{y}') (x^,y^)↔(x^′,y^′) 所以我们得到一个带约束的优化:
( x − x ^ ) 2 + ( y − y ^ ) 2 + ( x ′ − x ^ ′ ) 2 + ( y − y ^ ′ ) 2 (x-\hat{x})^2 + (y-\hat{y})^2 + (x'-\hat{x}')^2 + (y-\hat{y}')^2 (x−x^)2+(y−y^)2+(x′−x^′)2+(y−y^′)2
同时:
( x ^ , y ^ , 1 ) F A ( x ^ ′ , y ^ ′ , 1 ) = 0 (\hat{x},\hat{y},1) F_A (\hat{x}',\hat{y}',1) = 0 (x^,y^,1)FA(x^′,y^′,1)=0
那么怎么样求解呢? 除了几何解释的求解法以外,我们还可以用sampson损失函数:
( x ^ ′ y ^ ′ x ^ y ^ ) = ( x ′ y ′ x y ) − a x ′ + b y ′ + c x + d y + e a 2 + b 2 + c 2 + d 2 ( a b c d ) \begin{pmatrix} \hat{x}' \\ \hat{y}' \\ \hat{x} \\ \hat{y} \\ \end{pmatrix}= \begin{pmatrix} x' \\ y' \\ x \\ y \\ \end{pmatrix}- \frac{ax'+by'+cx+dy+e}{a^2+b^2+c^2+d^2} \begin{pmatrix} a \\ b \\ c \\ d \\ \end{pmatrix} x^′y^′x^y^ = x′y′xy −a2+b2+c2+d2ax′+by′+cx+dy+e abcd
14.5 Affine reconstruction
假设我们有多于4对对应点 x i ↔ x i ′ x_i \leftrightarrow x'_i xi↔xi′,如果摄像机是投影摄像机,那么重建的结果就是投影重建。现在如果摄像机是仿射摄像机,那么重建结果就是仿射重建。本节就来说明这个结果。
我们假设空间中有4个不共面的点 X i X_i Xi,我们选择 X 0 X_0 X0作为原点。然后我们构造三个坐标轴,表示为 E i ~ = X i ~ − X 0 ~ \tilde{E_i} = \tilde{X_i} - \tilde{X_0} Ei~=Xi~−X0~。所以对于一个空间中的点 X = ( x , y , z ) X=(x,y,z) X=(x,y,z),它的坐标就可以表示为:
X ~ = X 0 + x E 1 ~ + y E 2 ~ + z E 3 ~ \tilde{X} = X_0 + x \tilde{E_1} + y \tilde{E_2} + z \tilde{E_3} X~=X0+xE1~+yE2~+zE3~
X 0 ~ = ( 0 0 0 ) X 1 ~ = ( 1 0 0 ) X 2 ~ = ( 0 1 0 ) X 3 ~ = ( 0 0 1 ) \tilde{X_0} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ \end{pmatrix} \tilde{X_1} = \begin{pmatrix} 1 \\ 0 \\ 0 \\ \end{pmatrix} \tilde{X_2} = \begin{pmatrix} 0 \\ 1 \\ 0 \\ \end{pmatrix} \tilde{X_3} = \begin{pmatrix} 0 \\ 0 \\ 1 \\ \end{pmatrix} X0~= 000 X1~= 100 X2~= 010 X3~= 001
有了这些公式以后,我们来看看仿射摄像机的投影过程长什么样,该过程可以被表达为:
x ~ = M 2 × 3 X ~ + t ~ \tilde{x} = M_{2 \times 3} \tilde{X} + \tilde{t} x~=M2×3X~+t~
所以说上文的 E i ~ \tilde{E_i} Ei~可以被表示成 e i ~ = M 2 × 3 E i ~ \tilde{e_i}=M_{2 \times 3} \tilde{E_i} ei~=M2×3Ei~,那么对于空间点 X = ( x , y , z ) X=(x,y,z) X=(x,y,z),它在第一幅图像中就可以被表示成:
x e 1 ~ + y e 2 ~ + z e 3 ~ x \tilde{e_1} + y \tilde{e_2} + z \tilde{e_3} xe1~+ye2~+ze3~
对于第二幅图象, e ~ ′ = M 2 × 3 E i ~ \tilde{e}'=M_{2 \times 3} \tilde{E_i} e~′=M2×3Ei~,那么 X X X在第二幅图象中就可以表示成:
x e 1 ~ ′ + y e 2 ~ ′ + z e 3 ~ ′ x \tilde{e_1}' + y \tilde{e_2}' + z \tilde{e_3}' xe1~′+ye2~′+ze3~′
14.6 Necker reversal and the bas-relief ambiguity
本章主要讲述在已标定摄像机的情况下,只用两个图像会产生一系列的歧义.
主要是两种歧义:
-
Necker reversal
主要原因是物体旋转 ρ \rho ρ 和旋转 − ρ -\rho −ρ的镜像,在affine摄像机下会产生同样的投影图像。如果是透视投影,那么每个点都会有不同的深度,所以这种歧义就没有了。 -
The bas-relief ambiguity
主要原因是摄像机进行一个旋转后,从光心出发的光线依旧相交于同一点。这样导致深度和旋转角度是不确定的,参见P357,fig14.9(b),看图容易理解一点。
14.7 Computing the motion
本节主要讲述如果从 F A F_A FA中计算摄像机运动。
我们首先把摄像机的运动表示为 R = R ρ R θ R=R_{\rho}R_{\theta} R=RρRθ。
θ \theta θ是绕视线方向的旋转角度. ρ \rho ρ是绕与图像平面平行的轴的旋转角度, ϕ \phi ϕ是图像x轴的旋转角度.具体参见P358图14.10,图像之间除了旋转还有一个缩放因子 s s s。
具体的推导过程我们省略,在这里直接给出公式:
tan ϕ = b a tan ( ϕ − θ ) = d c s 2 = c 2 + d 2 a 2 + b 2 \tan \phi = \frac{b}{a} \\ \tan (\phi - \theta) = \frac{d}{c} \\ s^2 = \frac{c^2+d^2}{a^2+b^2} tanϕ=abtan(ϕ−θ)=cds2=a2+b2c2+d2
a , b , c , d a,b,c,d a,b,c,d就是黄金标准算法里提到的四维空间的平面参数。