【计算机视觉】多视几何之基础矩阵估计

阅读量：

一、对极几何

1.1 对极几何的概念

1.2 基本模型

1.3 对极几何约束

二、基础矩阵（F）

三、基础矩阵估计方法

3.1 八点估算法

3.2 ransac随机采样一致法估算

四、基础矩阵估计代码实现

4.1 代码

4.2 代码运行结果：

一、对极几何

1.1 对极几何的概念

对极几何（Epipolar Geometry）是Structure from Motion问题中，在两个相机位置产生的两幅图像的之间存在的一种特殊几何关系**，** 是sfm问题中2D-2D求解两帧间相机姿态的基本模型。

1.2 基本模型

c0、c1为两个相机中心；p为空间中一点；p在c0、c1对应像平面上的投影分别为x0、x1。

极平面（Epipolar Plane） ：两个相机坐标原点C0、C1和物体P三点组成的平面。

对极平面束（epipolar pencil） ：以基线为轴的平面束

极线（Epipolar Lines） ：极平面和两个像平面的交线，I0和I1

极点（Epipoles） :e0为右相机原点在左像平面的投影;e1为左相机原点在右像平面的投影

1.3 对极几何约束

下面等式称为对极约束（Epipolar Constraint）

对极约束的几何意义：x1、t、Rx0三者混合积为0，表示这三个向量共面，即上图中三角形的三边共面。

令E=t×RE=t×R，得到对极约束的新形式：

E称为本质矩阵（Esential Matrix ），由外参数R和t决定。

本质矩阵的几何意义：，即x1在直线上，表示E将x0投影到另一帧图像中的直线l1上。

二、基础矩阵（F）

基本矩阵 就表示了这种从点到直线的射影映射关系。

存在这么一个矩阵F，使得空间中不在两图像平面上的任意点X分别在两图像的投影坐标x,x'满足等式(x')TFx=0，即x'的转置乘以F，再乘以x的结果为0，那么F就是左边图像到右边图像的基本矩阵，从公式上可以看出基本矩阵是有方向的，右图到左图的基本矩阵就是F的转置。F矩阵有如下性质：
1、秩为2；
2、F矩阵是一个7个自由度的33矩阵(33矩阵本身9个自由度，因为相差一个常数因子和行列式值为0两个条件，减掉2个自由度)，相差一个常数因此的意思是：kF(k!=0)也是基本矩阵，也就是说如果F是基本矩阵，那么kF也是基本矩阵，所以基本矩阵不唯一，在相差一个倍数的前提下是唯一的，也就是我们可以固定矩阵中某一个非零元素的值，这样自然少一个自由度。

三、基础矩阵估计方法

3.1 八点估算法

八点法是通过对应点来计算基础矩阵的算法。
外极约束可以写成线性系统的形式：

其中，ff包含FF的元素，Xi1=[xi1,yi1,wi1]X1i=[x1i,y1i,w1i]和Xi2=[xi2,yi2,wi2]X2i=[x2i,y2i,w2i]是一对图像点，共有nn对对应点。
基础矩阵中共有9个元素，由于尺度是任意的，所以只需要8个方程。
我们通常用SVD算法来计算最小二乘解。由于算法得出的解可能秩不为2（基础矩阵的秩小于等于2），所以我们通过将最后一个奇异值置0来得到秩最接近2的基础矩阵。

3.2 ransac随机采样一致法估算

ransac方法是一种鲁棒性算法。该方法使用是满足可行条件下尽可能少的点集来计算初始解，然后检验点集初始解的吻合程度，从而选择最佳吻合点集为内点。该方法的思想是通过迭代随机抽样来找到能够使正确的匹配所占比例最高的最小点集。求解基础矩阵的步骤如下：
1、从m个点对中随机抽取k个子样本，其中每个子样本中含有8个点对。

2、对每个子样本用8点算法计算基础矩阵Fi。

3、对计算出来的Fi，从而计算其中每对匹配点到各自极线距离的平方和dj，统计dj<T匹配点对的估算，T为开始就设定的阈值。

4、选择满足dj<T匹配点最多的点对并根据此基础矩阵剔除误差较大的点，得出内点，并重新估计基础矩阵。

给定NN个数据点组成的集合PP，假设集合中大多数的点都是可以通过一个模型来产生的，且最少通过nn个点可以拟合出模型的参数，则可以通过以下的迭代方式拟合该参数：

从集合PP中随机选择nn个点
使用这nn个点估计出一个模型MM
对PP中剩余的点，计算每个点与模型MM的距离，距离超过阈值则认为是外点；不超过阈值则认为是内点，并记录该模型MM所对应的内点的个数mimi

将上面步骤重复kk次，选择最大mimi所对应的模型MM作为最终结果。

四、基础矩阵估计代码实现

4.1 代码

下面给出的代码是使用ransac算法。

复制代码

 from PIL import Image

    
 from numpy import *
    
 from pylab import *
    
 import numpy as np
    
 import camera
    
 import homography
    
 import sfm
    
 import sift
    
  
    
 im1 = array(Image.open('111.jpg'))
    
 sift.process_image('111.jpg', 'im1.sift')
    
  
    
 im2 = array(Image.open('222.jpg'))
    
 sift.process_image('222.jpg', 'im2.sift')
    
 l1, d1 = sift.read_features_from_file('im1.sift')
    
 l2, d2 = sift.read_features_from_file('im2.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):
    
     """ Robust estimation of a fundamental matrix F from point
    
     correspondences using RANSAC (ransac.py from
    
     http://www.scipy.org/Cookbook/RANSAC).
    
   36.     input: x1, x2 (3*n arrays) points in hom. coordinates. """
    
  
    
     import ransac
    
     data = np.vstack((x1, x2))
    
     d = 10 # 20 is the riginalo
    
     # 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：',F)
    
 P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
    
 P2 = sfm.compute_P_from_fundamental(F)
    
 print ('P2：',P2)
    
  
    
 # triangulate inliers and remove points not in front of both cameras
    
 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('333.jpg'))
    
 sift.process_image('333.jpg', 'im3.sift')
    
 l3, d3 = sift.read_features_from_file('im3.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 ('P2：',P2)
    
 print ('P3：',P3)