引言

本章将学习如何处理多个视图，以及如何利用多个视图的几何关系来恢复照相机位置信息及三维结构，通过在不同视点拍摄的图像，可以利用特征匹配计算出三维场景点及相机位置，本章将展示三维重建的完整例子。

5.1 外极几何

多视图几何是利用在不同视点所拍摄图像间的关系，研究照相机之间或特征之间的关系。图像特征通常是兴趣点，本章使用的也是兴趣点特征，我们可以应用于双视图几何。

可由一个场景的两个视图及对应图像点，根据照相机的空间相对位置、照相机的性质、三维场景点的位置，得到一些几何约束，我们通过外极几何来描述几何关系

三维场景点X经过单应性矩阵H变换后，在照相机$P{H}^{-1}$里得到的图像点和X在照相机P里得到的图像点相同，可以描述为：
$$\lambda x=PX=P{H}^{-1}HX=\hat{P}\hat{X}$$

所以，分析双视图几何关系时，可以将相对位置用单应性矩阵加以简化，即通过单应性矩阵变化了坐标系，一个做法是，将原点和坐标轴与第一个照相机对齐（平移向量为0）：$$P_1=K_1[I|0],P_2=K_2[R|t]$$
利用这些参数矩阵，可以找到X的投影点$x_1,x_2$,分别对应于投影矩阵$P_1,P_2$。我们可以从寻找对应的图像出发，恢复照相机参数矩阵，也是本章的一大目的。

同一个图像点经过不同的投影矩阵产生的不同投影点必须满足：
$$x_2^{T}F{x}_1=0$$ 其中：$$F=K_2^{-t}{S}_{t}R{K}_1^{-1}$$
$S_t为反对称矩阵$

$$\left[\begin{array}{ccc}0& -t_3& t_2\\ & & \\ t_3& 0& -t_1\\ & & \\ -t_2& t_1& 0\end{array}\right]$$
由平移t组成。F称为基础矩阵，由于反对称矩阵行列式为0，故F的秩小于等于2。

我们可以借助F恢复照相机参数，F可以从对应的投影图像点计算（x），K未知情况下，可以恢复出投影变换矩阵§，K已知，可以在三维重建中正确表示距离和角度。

由$x_2^{T}F{x}_1=l_1^T{x}_1=0$,找到第一幅图像的一条直线$l_1^T{x}_1=0$，第二个点在第一幅图像中的对应点一定在这条线（$x_2$的外极线）上，两条外极线都经过一个外极点e，是另一个照相机光心对应的图像点，$F{e}_1=e_2^{t}F=0$。

5.1.1 简单的数据集

我们需要一个带有图像点、三维点和照相机参数矩阵的数据集。

import camera
from PIL import Image
from numpy import *
from pylab import *
​
# 载入一些图像
im1 = array(Image.open('images/001.jpg'))
im2 = array(Image.open('images/002.jpg'))
# 载入每个视图的二维点到列表中
points2D = [loadtxt('2D/00'+str(i+1)+'.corners').T for i in range(3)]
# 载入三维点
points3D = loadtxt('3D/p3d').T
# 载入对应
corr = genfromtxt('2D/nview-corners',dtype='int')
# 载入照相机矩阵到 Camera 对象列表中
P = [camera.Camera(loadtxt('2D/00'+str(i+1)+'.P')) for i in range(3)]
​
# 将三维点转换成齐次坐标表示，并投影
X = vstack( (points3D,ones(points3D.shape[1])) )
x = P[0].project(X)
# 在视图 1 中绘制点
figure()
imshow(im1)
plot(points2D[0][0],points2D[0][1],'*')
axis('off')
figure()
imshow(im1)
plot(x[0],x[1],'r.')
axis('off')
show()

points2D有三个视图，我们选择第一幅图的坐标，上面的图片显示了矩阵切片的过程；X经过三维点转换为齐次坐标（vstact(,ones(…shape[1]))，再投影成x（二维点，为齐次坐标），x[0],x[1]是x、y坐标点，x[2]均为1组成。

因为三维点文件points3D的缘故，第二幅图会比第一幅图多一些点（视图2、3重建而来）。

5.1.2 用Matplotlib绘制三维数据

使用mpplot3d工作绘制三维点，用于可视化三维重建结果。

from mpl_toolkits.mplot3d import axes3d

fig = figure()
ax = fig.gca(projection="3d")
# 生成三维样本点
X,Y,Z = axes3d.get_test_data(0.25)
# 在三维中绘制点
ax.plot(X.flatten(),Y.flatten(),Z.flatten(),'yo')


show()

注：0.25是空间间隔参数，产生均匀的采样点。

给出Merton建筑的样本数据来观察三维点的结果。

# 绘制三维点
from mpl_toolkits.mplot3d import axes3

fig = figure()
ax = fig.gca(projection='3d')
ax.plot(points3D[0],points3D[1],points3D[2],'k.')

5.1.3 计算F：八点法

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

f包含F的元素，$x_1^i=[x_1^i,y_1^i,w_1^i]与x_2^i=[x_2^i,y_2^i,w_2^i]$是一对图像点，n对。由于是任意尺度，基础矩阵中有9个元素，所以需要8个对应点来计算基础矩阵F。

八点法中最小化||Af||的函数：

def compute_fundamental(x1,x2):
 """ 使用归一化的八点算法，从对应点（x1，x2 3×n 的数组）中计算基础矩阵
 每行由如下构成：
 [x'*x，x'*y' x', y'*x, y'*y, y', x, y, 1]"""
 n = x1.shape[1]
 
 if x2.shape[1] != n:
    raise ValueError("Number of points don't match.")
 # 创建方程对应的矩阵
 A = zeros((n,9))
 for i in range(n):
    A[i] = [x1[0,i]*x2[0,i], x1[0,i]*x2[1,i], x1[0,i]*x2[2,i],
 x1[1,i]*x2[0,i], x1[1,i]*x2[1,i], x1[1,i]*x2[2,i],
 x1[2,i]*x2[0,i], x1[2,i]*x2[1,i], x1[2,i]*x2[2,i] ]
 # 计算线性最小二乘解
 U,S,V = linalg.svd(A)
 # V的最后一行有9个元素
 F = V[-1].reshape(3,3)
 # 受限 F
 # 通过将最后一个奇异值置 0，使秩为 2
 U,S,V = linalg.svd(F)
 S[2] = 0
 F = dot(U,dot(diag(S),V))
 return F

由于基础矩阵秩小于等于2，我们通过将最后一个奇异值置0来得到最接近与2的基础矩阵。

5.1.4 外极点与外极线

外极点们组$F{e}_1=0$，因此可以通过计算F的零空间来得到。

def compute_epipole(F):
 """ 从基础矩阵 F 中计算右极点（可以使用 F.T 获得左极点）"""
 # 返回 F 的零空间（Fx=0）
 U,S,V = linalg.svd(F)
 e = V[-1]
 # 归一化
 return e/e[2]

若想获得另一幅图像的外极点，只需将F转置后输入上述函数即可。

在Merton数据集的前两个视图上运行这两个函数：

import sfm

# 在前两个视图中点的索引
# 第一/二个视图中>0的部分取交集
ndx = (corr[:,0]>=0) & (corr[:,1]>=0)
# 获得坐标，并将其用齐次坐标表示
x1 = points2D[0][:,corr[ndx,0]]
x1 = vstack( (x1,ones(x1.shape[1])) )
x2 = points2D[1][:,corr[ndx,1]]
x2 = vstack( (x2,ones(x2.shape[1])) )
# 计算 F
F = sfm.compute_fundamental(x1,x2)
# 计算极点
e = sfm.compute_epipole(F)
# 绘制图像
figure()
imshow(im1)
# 分别绘制每条线，这样会绘制出很漂亮的颜色
for i in range(5):
    plot_epipolar_line(im1,F,x2[:,i],e,False)
axis('off')
figure()
imshow(im2)
# 分别绘制每个点，这样会绘制出和线同样的颜色
for i in range(5):
    plot(x2[0,i],x2[1,i],'o')
axis('off')
show()

第一个视图画出了前5个外极线，第二个视图中画出了对应匹配点，可以看到，这些线在图片外左侧位置将相交于一点。外极线上一定存在着另一个图像的对应点。

5.2 照相机和三维结构的计算

这一节，将简单介绍计算照相机参数和三维结构的工具。

5.2.1 三角剖分

给定照相机参数模型，图像可通过三角剖分来恢复出这些点的三维位置。照相机方程关系定义如下：

通过SVD算法来得到三维点的最小二乘估值。

def triangulate_point(x1,x2,P1,P2):
 """ 使用最小二乘解，绘制点对的三角剖分 """
 M = zeros((6,6))
 M[:3,:4] = P1
 M[3:,:4] = P2
 M[:3,4] = -x1
 M[3:,5] = -x2
 U,S,V = linalg.svd(M)
 X = V[-1,:4]
 return X / X[3]

最小二乘解的前4个值就是齐次坐标系下的三维坐标

增加下列函数实现多个点的三角剖分

def triangulate(x1,x2,P1,P2):
 """ x1 和 x2（3×n 的齐次坐标表示）中点的二视图三角剖分 """
 n = x1.shape[1]
 if x2.shape[1] != n:

raise ValueError("Number of points don't match.")
 X = [ triangulate_point(x1[:,i],x2[:,i],P1,P2) for i in range(n)]
 return array(X).T

利用下面的代码来实现 Merton1 数据集上的三角剖分

import sfm
# 前两个视图中点的索引
ndx = (corr[:,0]>=0) & (corr[:,1]>=0)
# 获取坐标，并用齐次坐标表示
x1 = points2D[0][:,corr[ndx,0]]
x1 = vstack( (x1,ones(x1.shape[1])) )
x2 = points2D[1][:,corr[ndx,1]]
x2 = vstack( (x2,ones(x2.shape[1])) )
Xtrue = points3D[:,ndx]
Xtrue = vstack( (Xtrue,ones(Xtrue.shape[1])) )
# 检查前三个点
Xest = sfm.triangulate(x1,x2,P[0].P,P[1].P)
print Xest[:,:3]
print Xtrue[:,:3]
# 绘制图像
from mpl_toolkits.mplot3d import axes3d
fig = figure()
ax = fig.gca(projection='3d')
ax.plot(Xest[0],Xest[1],Xest[2],'ko')
ax.plot(Xtrue[0],Xtrue[1],Xtrue[2],'r.')
axis('equal')
show()

实际点与估计点



三维图像点与实际图像点很接近。

5.2.2 由三维点计算照相机矩阵

本质上，这是三角剖分的逆问题，有时我们将其称为照相机反切法，恢复照相机矩阵同样是一个最小二乘问题。
按照$${\lambda }_i{x}_i=P{X}_i$$
投影到图像点$x_i=[x_i,y_i,1]$,相应点满足下列关系：

def compute_P(x,X): 
    n = x.shape[1]
    if X.shape[1] != n:
        raise ValueError("Number of points don't match.")
   # 创建用于计算 DLT 解的矩阵
    M = zeros((3*n,12+n))
    for i in range(n):
        M[3*i,0:4] = X[:,i]
        M[3*i+1,4:8] = X[:,i]
        M[3*i+2,8:12] = X[:,i]
        M[3*i:3*i+3,i+12] = -x[:,i]
        U,S,V = linalg.svd(M)
    return V[-1,:12].reshape((3,4))

最后一个特征向量的前12个元素是照相机矩阵的元素。

下面的代码会选出第一个视图中的一些可见点，将它们转换维齐次坐标表示，然后估计照相机矩阵。

import sfm, camera
corr = corr[:,0] # 视图 1
ndx3D = where(corr>=0)[0] # 丢失的数值为 -1
ndx2D = corr[ndx3D]
# 选取可见点，并用齐次坐标表示
x = points2D[0][:,ndx2D] # 视图 1
x = vstack( (x,ones(x.shape[1])) )
X = points3D[:,ndx3D]
X = vstack( (X,ones(X.shape[1])) )
# 估计 P
Pest = camera.Camera(sfm.compute_P(x,X))
# 比较！
print (Pest.P / Pest.P[2,3]) 
print (P[0].P / P[0].P[2,3])
# 投影
xest = Pest.project(X)
# 绘制图像
figure()
imshow(im1)
plot(x[0],x[1],'bo')
plot(xest[0],xest[1],'r.')
axis('off')
show()

真实点用圆圈表示，照相机矩阵投影点用圆表示，结果基本相同。

5.3 多视图重建

假设照相机已经标定，计算重建可以分为下面 4 个步骤：
(1) 检测特征点，然后在两幅图像间匹配；
(2) 由匹配计算基础矩阵；
(3) 由基础矩阵计算照相机矩阵；
(4) 三角剖分这些三维点。

已完成好4个步骤的代码，但当图像间的点对应包含不正确的匹配时，需要一个稳健的方法来计算基础矩阵。

5.3.1 稳健估计基础矩阵

使用Ransac方法，结合八点算法，当场景点位于平面上时，就不能使用该算法。

RansacModel类

class RansacModel(object):
    def __init__(self,debug = False):
        self.debug = debug
        
    def fit(self,data):
        data=data.T
        # 数据分成两个点集,x1分为前三行，x2分为三行之后
        x1 = data[:3,:8]
        x2 = data[3:,:8]
        
        F = compute_fundamental_normalized(x1,x2)
        return F

    def get_error(self,data,F):
        data = data.T
        x1 = data[:3]
        x2 = data[3:]
        # 将 Sampson 距离用作误差度量
        Fx1 = dot(F,x1)
        Fx2 = dot(F,x2)
        denom = Fx1[0]**2 + Fx1[1]**2 + Fx2[0]**2 + Fx2[1]**2
        err = ( diag(dot(x1.T,dot(F,x2))) )**2 / denom
        # 返回每个点的误差
        return err

fit()会选择8个点，然后使用归一化的八点算法。

def compute_fundamental_normalized(x1,x2):
    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")
 # 归一化图像坐标
    x1 = x1 / x1[2]
    mean_1 = mean(x1[:2],axis=1)
    S1 = sqrt(2) / std(x1[:2])
    T1 = array([[S1,0,-S1*mean_1[0]],[0,S1,-S1*mean_1[1]],[0,0,1]])
    x1 = dot(T1,x1)
    x2 = x2 / x2[2]
    mean_2 = mean(x2[:2],axis=1)
    S2 = sqrt(2) / std(x2[:2])
    T2 = array([[S2,0,-S2*mean_2[0]],[0,S2,-S2*mean_2[1]],[0,0,1]])
    x2 = dot(T2,x2)
 # 使用归一化的坐标计算F
    F = compute_fundamental(x1,x2)
 # 反归一化
    F = dot(T1.T,dot(F,T2))
 # F[2,2]是最后一个坐标
    return F/F[2,2]

将图像归一化为零均值固定方差。

在F中使用Ransac单稳健估计。

def F_from_ransac(x1,x2,model,maxiter=5000,match_theshold=1e-6):
    import ransac
    data = vstack((x1,x2))
    # 计算F，并返回正确点索引
    F,ransac_data = ransac.ransac(data.T,model,8,maxiter,match_theshold,20,
    return_all=True)
    return F, ransac_data['inliers']

5.3.2 三维重建示例

本节中，我们使用已标定矩阵照相机拍摄的两幅图像，观察重建三维场景。我们将代码分成若干块，首先提取、匹配特征，然后估计基础矩阵和照相机矩阵。

import homography
import sfm
import sift

# 标定矩阵
K = array([[2394,0,932],[0,2398,628],[0,0,1]])
# 载入图像，并计算特征
im1 = array(Image.open('alcatraz1.jpg'))
sift.process_image('alcatraz1.jpg','im1.sift')
l1,d1 = sift.read_features_from_file('im1.sift')
im2 = array(Image.open('alcatraz2.jpg'))
sift.process_image('alcatraz2.jpg','im2.sift')
l2,d2 = sift.read_features_from_file('im2.sift')
# 匹配特征
matches = sift.match_twosided(d1,d2)
ndx = matches.nonzero()[0]
# 使用齐次坐标表示，并使用 inv(K) 归一化
x1 = homography.make_homog(l1[ndx,:2].T)
ndx2 = [int(matches[i]) for i in ndx]
x2 = homography.make_homog(l2[ndx2,:2].T)
x1n = dot(inv(K),x1)
x2n = dot(inv(K),x2)
# 使用 RANSAC 方法估计 E
model = sfm.RansacModel()
E,inliers = sfm.F_from_ransac(x1n,x2n,model)
# 计算照相机矩阵（P2 是 4 个解的列表）
P1 = array([[1,0,0,0],[0,1,0,0],[0,0,1,0]])
P2 = sfm.compute_P_from_essential(E)

P2的结果是照相机矩阵的四个可能解（K已标定的情况下度量重建）

从照相机的列表中，挑选出经过三角剖分后，含有在两个照相机前最多场景点的照相机矩阵。

# 选取点在照相机前的解
ind = 0
maxres = 0
for i in range(4):
 # 三角剖分正确点，并计算每个照相机的深度
 X = sfm.triangulate(x1n[:,inliers],x2n[:,inliers],P1,P2[i])
 # 第三个数值表示深度
 d1 = dot(P1,X)[2]
 d2 = dot(P2[i],X)[2]
 if sum(d1>0)+sum(d2>0) > maxres:
 maxres = sum(d1>0)+sum(d2>0)
 # 正确的三维点对应的下标
 ind = i
 infront = (d1>0) & (d2>0)
 # 三角剖分正确点，并移除不在所有照相机前面的点
 X = sfm.triangulate(x1n[:,inliers],x2n[:,inliers],P1,P2[ind])
 X = X[:,infront]

循环遍历P2[i]这四个解，对于正确的三维点进行三角剖分，用深度来排除在照相机前面的点。

接下来绘制三维重建（需要将第一个坐标轴取相反数）：

# 绘制三维图像
from mpl_toolkits.mplot3d import axes3d
fig = figure()
ax = fig.gca(projection='3d')
ax.plot(-X[0],X[1],X[2],'k.')
axis('off')

然后，每个视图中绘制二次投影。

# 绘制 X 的投影
import camera
# 绘制三维点
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2[ind])
x1p = cam1.project(X)
x2p = cam2.project(X)
# 反 K 归一化

x1p = dot(K,x1p)
x2p = dot(K,x2p)
figure()
imshow(im1)
gray()
plot(x1p[0],x1p[1],'o')
plot(x1[0],x1[1],'r.')
axis('off')
figure()
imshow(im2)
gray()
plot(x2p[0],x2p[1],'o')
plot(x2[0],x2[1],'r.')
axis('off')
show()

红点是特征点，蓝点是二次投影点，可以看出，不完全匹配，但是相当接近，三维重建的结果图很显然是错误的，可能与Python版本有关系。