基础矩阵

1. 基础矩阵简介

在计算机视觉中，基础矩阵F是一个3×3的矩阵，表达了立体像对的像点之间的对应关系。在对极几何中，对于立体像对中的一对同名点，它们的齐次化图像坐标分别为p与 p’，Fp表示一条必定经过p’的直线（极线）。这意味着立体像对的所有同名点对都满足：

F矩阵中蕴含了立体像对的两幅图像在拍摄时相互之间的空间几何关系（外参数）以及相机检校参数（内参数），包括旋转、位移、像主点坐标和焦距。因为F矩阵的秩为2，并且可以自由缩放（尺度化），所以只需7对同名点即可估算出F的值。本质矩阵与基础矩阵之间的关系可由下式表达：

其中K和K’分别为两个相机的内参数矩阵。

2. 基础矩阵推导

在双相机的拍摄场景中建立一个空间直角坐标系，称为世界坐标系（如图1中蓝色坐标系）。物点就是场景中物体表面上的点，比如说点P在世界坐标系中的坐标为 ：


相机的光心从物理上讲就是相机镜头组的光学中心。以光心为原点，主光轴为Z轴建立空间直角坐标系，称为相机坐标系（如图1中绿色和红色坐标系）。像平面在相机坐标系中的方程即为z=f，像点就是在物点在像平面上的投影，这个投影关系是透视投影。

用一句话来概括相机的拍摄模型，就是物点、像点、光心三点一线，此模型称为针孔相机模型。在此模型中，世界坐标系到左右相机坐标系的变换是刚性变换，即只包含旋转和平移，因此我们分别用增广矩阵[R|t]和[R’|t’]表示，其中R和R’是3x3的旋转矩阵，t和t’为平移向量。
令为P的齐次化坐标，那么物点P在左右相机坐标系下的坐标分别为

以一台相机为例，如图2所示，C为相机光心，Z轴为主轴。物点在相机坐标系下的坐标和以相片左下角为原点的像点坐标p有如下关系：


式中（x0,y0,f）为像主点在相机坐标系下的坐标。
设两相机内参数矩阵同为：

那么物点与像点之间的关系为：

又因为：

并令下式得：


由于物点、像点、光心三点一线，那么物点、一对同名点和2个光心这5个点一定处于同一个平面上，我们将这个平面称为𝜋平面。𝜋平面和像平面的交线称为极线l’。显然，左片上的每一个像点p对应于右片上的一条极线l’，且p’一定在l’上。两个相机光心的连线与右片像平面的交点称为极点，用e’表示。

在右片像平面内，极线 l’的方程可以表示为Ax+By+C=0。这个平面直线方程的一般式可以视为：

因此，我们可以用一个三维向量 (A,B,C)来表示极线l’，并且 l’的方程可以简单的由e’坐标向量与 p’坐标向量做向量积得到，即

其中：

令[e’]x表示向量积的矩阵形式，那么再将同名点之间的变换关系代入，得到极线的方程为：

因为p’在 l’上，所以显然有：

令

即得到：

3. 实验代码

3.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('C:/picture/3.jpg'))
sift.process_image('C:/picture/3.jpg', 'img1.sift')

im2 = array(Image.open('C:/picture/4.jpg'))
sift.process_image('C:/picture/4.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,
                                   7, 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[