对极几何以及基础矩阵原理求解（一）

一、P在不同坐标系中的表示

空间点P本身在世界坐标系中的表示：$$PW=⎛⎝⎜XWYWZW⎞⎠⎟PW=(XWYWZW)P_W=\left(\begin{array}{c}{X}_W\\ Y_W\\ Z_W\end{array}\right)$$x′TFx=0

F是左边图像到右边图像的基本矩阵，从公式上可以看出基本矩阵是有方向的，右图到左图的基本矩阵就是F的转置。F矩阵有如下性质：

1. 秩为2
2. F矩阵是一个7个自由度的3 * 3矩阵(3 * 3矩阵本身9个自由度，因为相差一个常数因子和行列式值为0两个条件，减掉2个自由度)，相差一个常数因子的意思是：kF (其中k不为0) 也是基本矩阵，也就是说如果F是基本矩阵，那么kF也是基本矩阵，所以基本矩阵不唯一，在相差一个倍数的前提下是唯一的，也就是我们可以固定矩阵中某一个非零元素的值，这样自然少一个自由度。

基础矩阵是对极几何的代数表达方式
基础矩阵描述了图像中任意对应点 x↔x’ 之间的约束关系

假设空间点X的位置未知，我们只知道X在左边图上的投影x的坐标位置，也知道基本矩阵F，但是X一定在射线Cx上，具体哪一点是没法知道的，也就是X可能是Cx上的任意一点(也就是轨迹的意思)，那么X在右图上的投影肯定也是一条直线l’（因为X在一条轨迹上）。即知道空间中某点在一个二维坐标位置以及两幅图之间的基本矩阵可以知道该点在另一个二维坐标的所在直线（还不能知道具体坐标）。所以说基本矩阵提供了三维点到二维点之间的约束条件。

基础矩阵之歌的视频链接，感兴趣的可以去看一下
https://v.youku.com/v_show/id_XMzEzOTc3ODQ0.html

三、基础矩阵的求解

这里用到了很多包你要提前导入

可以用8点算法估计基础矩阵
但是由于矩阵各留的数据尺度差异太大，最小二乘法得到的结果精度一般很低，所以先把x与x’的坐标进行归一化处理，将图像坐标变换到合适的范围内，相当于无量纲化？？
最后还原得到的基础矩阵就可以了。

目标

计算图像的特征匹配，并估计基础矩阵
选择三张图像，挑选其中一对图像，计算三维点和照相机矩阵，然后利用这些对应的三维点使用后方交汇法，计算第三张图像的照相机矩阵。

代码

以下是代码：
1.提取两幅图像的sift特征点（全部）
2.通过ransac随机一致性采样得到基础矩阵

from PIL import Image
from numpy import *
from pylab import *
import numpy as np

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

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

# Read features
im1 = array(Image.open('c.jpg'))
sift.process_image('c.jpg', 'im1.sift')

im2 = array(Image.open('d.jpg'))
sift.process_image('d.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).

    input: x1, x2 (3*n arrays) points in hom. coordinates. """

    #import ransac
    from PCV.tools import ransac

    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)

print P2
print F

# 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('c.jpg'))
sift.process_image('c.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 P3

print P1
print P2

~~

实验结果及其分析

第一组：中山纪念馆

实验结果：
报错：ValueError: did not meet fit acceptance criteria
是在ransac随机一致性采样中出错，没有采集到大于阈值的有效匹配点，我们将阈值条件放宽，放大很多倍数10^6倍还是不可以采集到有效数量的特征点，于是我换了一个角度拍摄中山纪念馆。在原来的两张图像a、c视角中间择中又拍摄一张图片b，将a与b进行匹配。
在原来的阈值下可以检测到有效的特征点了。

F是图1和图2之间的基础矩阵，P1、P2、P3分别是图像1、2、3的照相机矩阵
F=
[[ 6.98535699e-08 9.70288247e-06 -4.33762855e-03]
[-9.75128318e-06 2.96469330e-07 5.22746243e-03]
[ 4.39126050e-03 -7.70323885e-03 1.00000000e+00]]
P1=
[[1 0 0 0]
[0 1 0 0]
[0 0 1 0]]
P2=
[[-1.97759462e+00 2.38327221e+00 4.55921550e+02 7.79982073e+02]
[ 3.38327258e+00 -4.07733673e+00 -7.79977681e+02 4.55913847e+02]
[ 7.53622717e-03 4.67698579e-03 -8.01042467e+00 1.00000000e+00]]
P3=
[[-1.55915515e-01 1.10107243e-01 -2.06737508e-04 1.99347012e-01]
[ 6.68776368e-02 -2.30452052e-01 -3.72988540e-04 4.07818636e-02]
[ 8.38956147e-05 -3.75456712e-04 -4.74038220e-07 1.07336726e-04]]

第二组：室内
F是图1和图2之间的基础矩阵，P1、P2、P3分别是图像1、2、3的照相机矩阵
F=
[-3.53166718e-06 2.95908056e-05 -6.93899238e-03]
[-3.21103903e-05 -3.07027105e-06 1.49675657e-02]
[ 1.06807984e-02 -1.98205552e-02 1.00000000e+00]]
P1=
[[1 0 0 0]
[0 1 0 0]
[0 0 1 0]]
P2=
[[-1.77665613e+00 3.83222725e+00 2.56055054e+02 6.96387013e+02]
[ 4.83222065e+00 -1.04232505e+01 -6.96376332e+02 2.56035234e+02]
[ 2.15108840e-02 6.08329440e-03 -1.65374379e+01 1.00000000e+00]]
P3=
[[ 7.13172109e-04 2.29100260e-03 -9.08334155e-01 6.08818386e-02]
[ 5.53626143e-04 4.29519151e-04 -4.13444588e-01 1.66261332e-02]
[ 1.97076034e-06 2.22758095e-06 -1.65818292e-03 6.20972210e-05]]

                                </div>