【计算机视觉】基础矩阵

阅读量：

本课程材料介绍了外极几何和基础矩阵F的概念及其应用。外极几何描述了两个相机之间的空间关系，包括极平面、极点和对极线等核心概念。基础矩阵F是一个秩为2的齐次矩阵，在描述对极几何关系中起着重要作用。其性质包括与匹配点坐标的关系以及与单应变换之间的联系。通过八点算法可以估计出基本矩阵F，并通过SVD对其进行约束以确保秩为2。此外，材料还提供了使用Python和OpenCV库进行特征匹配、计算基础矩阵以及展示实验结果的代码实现，并讨论了不同场景下的实验结果及潜在问题。

一、外参数几何
二、基本矩阵F
- - 1.基本矩阵原理
- 2.基本矩阵性质
- 3.基本矩阵估计方法
  - - 3.1 八点估算法
      - 3.2 RANSAC随机取样一致法进行估计（略）
第三章基础矩阵的代码实现
- 第一节核心代码模块
  - 第二节实验验证结果
  - 第二节远距离室外场景
    - ** subsection 距离相近情况下的表现
    - ** subsection 距离远近对比分析
    - ** subsection 室内环境测试
- 四、实验总结

一、外极几何

采用两个相机分别从不同角度对同一物体拍照时，在这两张具有部分重叠的景物的照片之间必然存在着某种对应的几何关系被称为外极几何。

外极几何示意图：

在示意图中，
PO₁O₂是极平面，
相机中心分别为O₁和O₂，
空间任意一点P在其图像平面上的投影分别对应于p₁和p₂，
照相机中心连线即基线段OO₁和OO₂与图像平面相交于外极点e₁和e₂。
而对偶几何则刻画了这些元素间的对应关系。

二、基础矩阵 F

1.基础矩阵原理

基于对极几何理论的研究需求，在本研究中引入了基础矩阵F。其极线l'存在于另一幅图像中。与其匹配的点x'必然位于其对应的极线上。我们利用基础矩阵F建立了数学模型来描述这一映射关系。

假设第一幅图到第二幅图之间存在单应性变换H，则 x′=Hx。

又因为l’是表示过对极点e2和图像点x’的直线，可以表示为： l′=e2×x′=[e2]xx′=[e2]xHx。

所以可得 F=[e2]xH。

其中e′=[a1,a2,a3]Te′=[a1,a2,a3]T，则它的反对称阵定义为：

在图像空间中任意一点x，在已获得两个摄像机的射影矩阵P和P’的情况下，通过推导可得该点的反投影射线方程为：X(λ)=P+ x+λC

其中Moore-Penrose广义逆矩阵满足性质PP⁺ = I，并用于表示相机中心位置C。考虑一条经过相机中心的直线l上任意一点X(t)，其坐标可表示为X(t) = C + t*(PP⁺ - C)，其中t∈ℝ。特别地，在t=1时该直线与图像平面相交得到图像质心g₁；而当t→∞时对应的极限位置则位于图像边缘远端g₂处。

可以推得 F=[e′]xP′P+

2.基础矩阵性质

（1）基础矩阵是由秩2、自由度7所组成的齐次矩阵。
（2）当x和x'分别是两张图上的对应点时，则有x'TFx=0成立。
（3）对于点x而言，在其对应的对极线l'上必定满足l'=Fx的关系式。
（4）当e表示第二个摄像机光心在第一幅图像中的极点时，则必然满足Fe=0这一条件。

3.基础矩阵估算方法

3.1 八点估算法

基本矩阵是由该方程定义的： x′TFx=0

其中 $x \leftrightarrow x'$ 代表两张图像中任意一对对应的匹配点。每对匹配点都贡献了一个用于计算 $F$ 矩阵的一次方程。当至少拥有7个对应点（考虑到3×3齐次矩阵中的缩放自由度以及秩2约束条件）时，则可解出未知的 $F$ 矩阵。对于这些点的位置坐标，则通常采用齐次坐标形式表示为 $x = (x,y,1)^T$ 和 $x' = (x',y',1)^T$ 。

又因为F为：

所以可得到方程：

即相应方程式为

给定n组点的集合，我们有如下方程：

当存在明确且非零值的解时，系数矩阵A的秩至多为8。基于齐次性质，在这种情况下，若系数矩阵A的秩达到8，则其解仅相差一个比例因子。可以通过线性计算方法求取结果。

若观测点坐标存在噪声，则矩阵AA的秩可能超过8（即达到9）。这表明我们需要通过求取最小二乘解来恢复原始结构。具体而言，在这种情况下采用奇异值分解（SVD）方法能够有效地解决这一问题：通过计算系数矩阵A的最小奇异值对应的奇异向量即可获得所需的结果——该向量对应于A=UDVT分解后V矩阵的最后一列矢量。这种方法在约束条件下寻求最优解时表现突出，在计算机视觉领域中被广泛应用于估计基本矩阵的过程，并通常被称为8点算法（Eight Point Algorithm）。

计算后的矩阵F可能无法达到秩为2的要求。从而需要在原有矩阵F的基础上施加额外约束。利用奇异值分解方法则是一种可行的方式；令其具有形式 $F = U\Sigma V^T$ 即可满足条件。

因为要秩为2，所以取最后一个元素设置为0，则

最终的解

3.2 RANSAC随机取样一致法估算（略）

三、基础矩阵代码实现

1.相关代码

复制代码

    from PIL import Image
    from numpy import *
    from pylab import *
    from PCV.geometry import camera
    from PCV.geometry import homography
    from PCV.geometry import sfm
    from PCV.localdescriptors import sift
    from PCV.tools import ransac
    import numpy as np
    
    
    # Read features
    im1 = array(Image.open('pic1.jpg'))
    sift.process_image('imim1.jpg', 'img1.sift')
    
    im2 = array(Image.open('pic2.jpg'))
    sift.process_image('imim2.jpg', 'img2.sift')
    
    
    l1, d1 = sift.read_features_from_file('img1.sift')
    l2, d2 = sift.read_features_from_file('img2.sift')
    
    
    matches = sift.match_twosided(d1, d2)
    
    ndx = matches.nonzero()[0]
    x1 = homography.make_homog(l1[ndx, :2].T)
    ndx2 = [int(matches[i]) for i in ndx]
    x2 = homography.make_homog(l2[ndx2, :2].T)
    
    d1n = d1[ndx]
    d2n = d2[ndx2]
    x1n = x1.copy()
    x2n = x2.copy()
    
    
    #figure(figsize=(16,16))
    sift.plot_matches(im1, im2, l1, l2, matches, True)
    show()
    
    
    
    #def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
    def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
    
    data = np.vstack((x1, x2))
    d = 10 # 20 is the original
    # compute F and return with inlier index
    F, ransac_data = ransac.ransac(data.T, model,
                                   8, maxiter, match_threshold, d, return_all=True)
    return F, ransac_data['inliers']
    
    # find F through RANSAC
    model = sfm.RansacModel()
    F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-5)
    print (F)
    
    
    P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
    P2 = sfm.compute_P_from_fundamental(F)
    
    
    X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)
    # plot the projection of X
    cam1 = camera.Camera(P1)
    cam2 = camera.Camera(P2)
    x1p = cam1.project(X)
    x2p = cam2.project(X)
    
    # figure(figsize=(16, 16))
    imj = sift.appendimages(im1, im2)
    imj = vstack((imj, imj))
    
    imshow(imj)
    
    cols1 = im1.shape[1]
    rows1 = im1.shape[0]
    for i in range(len(x1p[0])):
    if (0<= x1p[0][i]<cols1) and (0<= x2p[0][i]<cols1) and (0<=x1p[1][i]<rows1) and (0<=x2p[1][i]<rows1):
        plot([x1p[0][i], x2p[0][i]+cols1],[x1p[1][i], x2p[1][i]],'c')
    axis('off')
    show()
    
    d1p = d1n[inliers]
    d2p = d2n[inliers]
    
    # Read features
    im3 = array(Image.open('pic2.jpg'))
    sift.process_image('pic2.jpg', '3.sift')
    l3, d3 = sift.read_features_from_file('3.sift')
    
    matches13 = sift.match_twosided(d1p, d3)
    
    ndx_13 = matches13.nonzero()[0]
    x1_13 = homography.make_homog(x1p[:, ndx_13])
    ndx2_13 = [int(matches13[i]) for i in ndx_13]
    x3_13 = homography.make_homog(l3[ndx2_13, :2].T)
    
    #figure(figsize=(16, 16))
    imj = sift.appendimages(im1, im3)
    imj = vstack((imj, imj))
    
    imshow(imj)
    
    cols1 = im1.shape[1]
    rows1 = im1.shape[0]
    for i in range(len(x1_13[0])):
    if (0<= x1_13[0][i]<cols1) and (0<= x3_13[0][i]<cols1) and (0<=x1_13[1][i]<rows1) and (0<=x3_13[1][i]<rows1):
        plot([x1_13[0][i], x3_13[0][i]+cols1],[x1_13[1][i], x3_13[1][i]],'c')
    axis('off')
    show()
    
    P3 = sfm.compute_P(x3_13, X[:, ndx_13])
    
    print ('P1',P1)
    print ('P1',P2)
    print ('P1',P3)