计算机视觉——对极几何与基础矩阵

一、外极几何

如果有一个场景的两个视图以及视图中的对应图像点，那么根据照相机的空间相对位置关系、照相机的性质以及三维场景点的位置，可以得到对这些图像点的一些几何关系约束。简单说就是用两个相机在不同位置拍摄同一物体，两张照片中的景物有重叠部分，那么理论上这两张照片会存在一定的对应关系即外极几何。其示意图如下所示：

二、基础矩阵

1.概述

为了描述对极几何，引入基础矩阵F。对于一幅图像上的点x（图上p1），在另一幅图像上存在对极线l’，并且在第二幅图像上，与x匹配的点x’（点p2）必然在l’上。为了表示x与l’关系，我们令基础矩阵F定义为： l′=Fx。

假设第一幅图到第二幅图之间存在单应性变换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，则它的反对称阵定义为：
在已知两个摄像机的射影矩阵P和P’时，对于图像上的一点x，可以得到其反投影射线方程为： X(λ)=P+ x+λC
其中P+ 是P的伪逆，即PP+ = I，C为摄像机中心。射线上的两点：P+ x（当λ=0）、C（当λ=∞）在第二个摄像机P’拍摄下，在第二幅视图上的点分别为P′P+ 和P′C，这两点过对极线l’，即 l′=(P′C)×(P′P+ x)=[e′]x(P′P+ )x
可以推得 F=[e′]xP′P+。

2.性质

（1） 基础矩阵是秩为2、自由度为7的齐次矩阵。
（2）若x与x’是两幅图上的对应点，那么x′TFx=0。
（3）l’是对应于x的对极线，l′=Fx。
（4）若e是第二个摄像机光心在第一幅图像上的极点，那么Fe=0 。

3.估算方法

3.1八点估算法

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

其中x↔x′是两幅图像的任意一对匹配点。由于每一组点的匹配提供了计算F系数的一个线性方程，当给定至少7个点（3×3的齐次矩阵减去一个尺度，以及一个秩为2的约束），方程就可以计算出未知的F。我们记点的坐标为x=(x,y,1)T，x′=(x′,y′,1)T
又因为F为：

所以可得到方程：

即相应方程式为

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

如果存在确定（非零）解，则系数矩阵A的秩最多是8。由于F是齐次矩阵，所以如果矩阵A的秩为8，则在差一个尺度因子的情况下解是唯一的。可以直接用线性算法解得。

如果由于点坐标存在噪声则矩阵AA的秩可能大于8（也就是等于9，由于A是n×9的矩阵）。这时候就需要求最小二乘解，这里就可以用SVD来求解，f的解就是系数矩阵A最小奇异值对应的奇异向量，也就是A奇异值分解后A=UDVT中矩阵V的最后一列矢量，这是在解矢量f在约束∥f∥下取∥Af∥最小的解。以上算法是解基本矩阵的基本方法，称为8点算法。

上述求解后的F不一定能满足秩为2的约束，因此还要在F的基础上加以约束。通过SVD分解可以解决，令F=UΣVT，则

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

最终的解

八点算法的优点：
线性求解，容易实现，运行速度快 。

八点算法的缺点：
对噪声敏感。

三、实验

代码
特征匹配：

# coding: utf-8
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
from PCV.geometry import homography, camera, sfm
from PCV.localdescriptors import sift

camera = reload(camera)
homography = reload(homography)
sfm = reload(sfm)
sift = reload(sift)

# 提取特征
im1 = array(Image.open('D:/picture/1.jpg'))
sift.process_image('D:/picture/1.jpg', 'im1.sift')

im2 = array(Image.open('D:/picture/2.jpg'))
sift.process_image('D:/picture/2.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