本章讲解如何处理多个视图，以及如何利用多个视图的几何关系来恢复照相机位置 信息和三维结构。通过在不同视点拍摄的图像，我们可以利用特征匹配来计算出三 维场景点以及照相机位置。本章会介绍一些基本的方法，展示一个三维重建的完整 例子;本章最后将介绍如何由立体图像进行致密深度重建。

1、外极几何

同一个图像点经过不同的投影矩阵产生的不同投影点必须满足:

St为外极约束条件。矩阵 F 为基础矩阵。基础矩阵可以由两照相机的参数矩阵(相对旋转 R 和平移 t)表示。基础矩阵跟本质矩阵的区别是推导采用的坐标系不同，基础矩阵在图像坐标系中推导而出，本质矩阵在相机坐标系中推导而出，所以基础矩阵由两个照相机的内外参决定（如上），而本质矩阵只由外参决定E=StR=t^R。从几何意义上来说，都是因为世界坐标系中任意一点X和基线C1C2共面（（XC1-C1C2）·XC1xC1C2=0）。外极几何示意图如下：

（1）一个简单的数据集

load_vggdata.py获取数据代码实现：

import camera
from pylab import *
from PIL import Image
import numpy as np
# 载入一些图像
im1 = array(Image.open('../data/MertonCollegeI/images/001.jpg')) 
im2 = array(Image.open('../data/MertonCollegeI/images/002.jpg'))
# 载入每个视图的二维点到列表中
points2D = [loadtxt('../data/MertonCollegeI/2D/00'+str(i+1)+'.corners').T for i in range(3)]
# 载入三维点
points3D = loadtxt('../data/MertonCollegeI/3D/p3d').T
# 载入对应
corr = genfromtxt('../data/MertonCollegeI/2D/nview-corners',dtype='int',missing_values='*')
# 载入照相机矩阵到 Camera 对象列表中
P = [camera.Camera(loadtxt('../data/MertonCollegeI/2D/00'+str(i+1)+'.P')) for i in range(3)]

调用load_vggdata.py获取所有的数据，下面是可视化这些数据的代码实现。

# execfile('load_vggdata.py')
# python3 删去了 execfile()，代替方法如下：
with open('load_vggdata.py','r') as f:
    exec(f.read())

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()

牛津 multi-view 数据集中的 Merton1 数据集:视图 1 与图像点(左);视图 1 与投 影的三维点(右)

（2）用Matplotlib绘制三维数据

Matplotlib 中的 mplot3d 工具 包可以方便地绘制出三维点、线、等轮廓线、表面以及其他基本图形组件，还可以 通过图像窗口控件实现三维旋转和缩放。示例代码：

import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import axes3d


# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname="/System/Library/Fonts/PingFang.ttc", size=14)

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

plt.title(u'mplot3d样本点绘制图', fontproperties=font)

plt.show()

# # 绘制三维点
# from mpl_toolkits.mplot3d import axes3d
with open('load_vggdata.py','r') as f:
    exec(f.read())

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

title(u'Merton样本数据绘制图', fontproperties=font)

plt.show()

使用 Matplottib工具包绘制的，mplot3d样本点和Mertoln样本数据绘制图

 

（3）计算F:八点法

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

其中，f 包含 F 的元素，x1和x2 是两组图像点，共有 n 对对应点。基础矩阵中有 9 个元素，由于尺度是任意的，所以只需要 8 个方程。因为算 法中需要 8 个对应点来计算基础矩阵 F，所以该算法叫做八点法。

（4）外极点和外极线

新建一个文件 sfm.py，写入下面八点法中最小化 ||Af||的函数compute_fundamental；外极点满足 F·e1=0，因此可以通过计算 F 的零空间来得到，用compute_epipole函数实现；我们在第一个视图中画出了前 5 个外极线，在第二个视图中画出了对应匹配点，用plot_epipolar_line函数实现。sfm.py代码实现：

from pylab import *

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) 
	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

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


def plot_epipolar_line(im,F,x,epipole=None,show_epipole=True):
	""" 在图像中，绘制外极点和外极线 F×x=0。F 是基础矩阵，x 是另一幅图像中的点 """
	m,n = im.shape[:2] 
	line = dot(F,x)

	# 外极线参数和值
	t = linspace(0,n,100)
	lt = array([(line[2]+line[0]*tt)/(-line[1]) for tt in t])
	# 仅仅处理位于图像内部的点和线 
	ndx = (lt>=0) & (lt<m) 
	plot(t[ndx],lt[ndx],linewidth=2)

	if show_epipole:
		if epipole is None:
			epipole = compute_epipole(F) 
		plot(epipole[0]/epipole[2],epipole[1]/epipole[2],'r*')

外极点和外极线代码实现：

import sfm

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname="/System/Library/Fonts/PingFang.ttc", size=14)

# execfile('load_vggdata.py')
# python3 删去了 execfile()，代替方法如下：
with open('load_vggdata.py','r') as f:
    exec(f.read())
# 在前两个视图中点的索引
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() 

subplot(121)
imshow(im1)
# 分别绘制每个点，这样会绘制出和线同样的颜色 
for i in range(5):
	plot(x2[0,i],x2[1,i],'o') 
title(u'外极点', fontproperties=font)
axis('off')

subplot(122)
imshow(im2)
# 分别绘制每条线，这样会绘制出很漂亮的颜色 
for i in range(5):
	sfm.plot_epipolar_line(im2,F,x2[:,i],e,False) 
title(u'外极线', fontproperties=font)
axis('off')

show()

视图 1 中的外极线示为 Merton1 数据视图 2 中 5 个点对应的外极线。

2、照相机和三维结构的计算

（1）三角剖分

给定照相机参数模型，图像点可以通过三角剖分来恢复出这些点的三维位置。基本的算法思想如下。对于两个照相机 P1 和 P2 的视图，三维实物点 X 的投影点为 x1 和 x2(这里用齐次坐标表示)，照相机方程定义了下列关系:

triangulate_point函数用于计算一个点对的最小二乘三角剖分，triangulate函数来实现多个点的三角剖分，把它们添加到 sfm.py 中。利用下面的代码来实现 Merton1 数据集上的三角剖分：

import sfm

# execfile('load_vggdata.py')
# python3 删去了 execfile()，代替方法如下：
with open('load_vggdata.py','r') as f:
    exec(f.read())
# 前两个视图中点的索引
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')
axis()

show()

如下图所示，使用照相机矩阵和点的对应关系来三角剖分这些数据点。黑色的圆圈表示估计的点， 红色的点表示真实点。

（2）由三维点计算照相机矩阵

每个三维点 Xi(齐次坐标系下)按照 λixi =PXi 投影到图像点 xi=[xi, yi,1]，相应的点满足下面的关系:

其中 p1、p2 和 p3 是矩阵 P 的三行。上面的式子可以写得更紧凑，如下所示: Mv=0

然后，我们可以使用 SVD 分解估计出照相机矩阵。利用上面讲述的矩阵操作，我们可以直接写出相应的代码compute_P函数，将该函数添加到 sfm.py 文件中。

下面，在我们的样本数据集上测试算法的性能。下面的代码会选出第一个视图中的 一些可见点(使用对应列表中缺失的数值)，将它们转换为齐次坐标表示，然后估计照相机矩阵:

import sfm, camera

# execfile('load_vggdata.py')
# python3 删去了 execfile()，代替方法如下：
with open('load_vggdata.py','r') as f:
    exec(f.read())
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()

如下图所示，为了检查照相机矩阵的正确性，将它们以归一化的格式(除以最后一个元素)打印 到控制台。彩图中真实点用圆圈表示，估计出的照相机投影点用点表示。

（3）由基础矩阵计算照相机矩阵

给定标定矩阵 K，我们可以将它的逆 K^(-1) 作用于图像点 ，因此，在新的图像坐标系下，照相机方程变为:

在新的图像坐标系下，点同样满足之前的基础矩阵方程:

在标定归一化的坐标系下，基础矩阵称为本质矩阵。为了区别为标定后的情况，以及归一化了的图像坐标，我们通常将其记为 E，而非 F。

3、多视图重建

下面让我们来看，如何使用上面的理论从多幅图像中计算出真实的三维重建。由于照相机的运动给我们提供了三维结构，所以这样计算三维重建的方法通常称为 SfM (Structure from Motion，从运动中恢复结构)。假设照相机已经标定，计算重建可以分为下面 4 个步骤:

(1) 检测特征点，然后在两幅图像间匹配;

(2) 由匹配计算基础矩阵;
(3) 由基础矩阵计算照相机矩阵;
(4) 三角剖分这些三维点。

我们已经具备了完成上面 4 个步骤所需的所有工具。但是当图像间的点对应包含不 正确的匹配时，我们需要一个稳健的方法来计算基础矩阵。

（1）稳健估计基础矩阵

当存在噪声和不正确的匹配时，我们需要估 计基础矩阵。和前面的方法一样，我们使用 RANSAC 方法，只不过这次结合了八 点算法。注意，八点算法在平面场景中会失效，所以，如果场景点都位于平面上， 就不能使用该算法。

将类RansacModel、归一化的八点算法函数compute_fundamental_normalized以及F_from_ransac函数添加到 sfm.py 文件中。我们返回最佳基础矩阵 F，以及正确点的索引，这样可以知道哪些匹配和 F 矩阵是一致的。与单应性矩阵估计相比，我们增加了默认的最大迭代次数，改变了 匹配的阈值(之前使用像素，现在使用 Sampson 距离来衡量)。

（2）三维重建示例

代码实现：

import homography
import sfm
import sift
from pylab import *
from PIL import Image
  
# 标定矩阵
K = array([[2394,0,932],[0,2398,628],[0,0,1]])
# 载入图像，并计算特征
im1 = array(Image.open('../data/alcatraz1.jpg')) 
sift.process_image('../data/alcatraz1.jpg','im1.sift') 
l1,d1 = sift.read_features_from_file('im1.sift')
im2 = array(Image.open('../data/alcatraz2.jpg')) 
sift.process_image('../data/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)

# 选取点在照相机前的解 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]

# 绘制三维图像
from mpl_toolkits.mplot3d import axes3d

fig = figure()

subplot(223)
ax = fig.gca(projection='3d') 
ax.plot(-X[0],X[1],X[2],'k.') 
axis('off')

subplot(224)
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)

 
subplot(221)
imshow(im1)
gray() 
plot(x1p[0],x1p[1],'o') 
plot(x1[0],x1[1],'r.') 
axis('off')

 
subplot(222)
imshow(im2)
gray() 
plot(x2p[0],x2p[1],'o') 
plot(x2[0],x2[1],'r.') 
axis('off')

show()

运行结果：

 

（3）多视图的扩展示例

多视图重建有一些步骤和进一步的扩展，但是不能在本书中全部介绍。为了方便读者进一步阅读，下面提供关于一些有关内容及其参考文献。

1. 多视图

当我们有同一场景的多个视图时，三维重建会变得更准确，包括更多的细节信息。 因为基础矩阵只和一对视图相关，所以该过程带有很多图像，和之前的处理有些 不同。

对于视频序列，我们可以使用时序关系，在连续的帧对中匹配特征。相对方位需要 从每个帧对增量地加入下一个帧对，(与图 3-12 中我们在全景图例子中加入单应性 矩阵相同)。该方法非常有效，同时可以使用跟踪有效地找到对应点(关于跟踪的更 多知识，参阅 10.4 节)。一个问题是误差会在加入的视图中积累。该问题可以由最 后的优化步骤来解决，参见下文。

对于静止的图像，一个办法是找到一个中央参考视图，然后计算与之有关的所有其 他照相机矩阵。另一个办法是计算一个图像对的照相机矩阵和三维重建，然后增量 地加入新的图像和三维点，参见参考文献 [34]。另外，还有一些方法可以同时由三 个视图计算三维重建和照相机位置(参阅参考文献 [13]);但除此之外，我们还需 要使用增量的方法。

2. 光束法平差
       从图 5-8 简单三维重建的例子，我们可以清楚地看出，恢复出的点的位置和由估计 的基础矩阵计算出的照相机矩阵都存在误差。在多个视图的计算中，这些误差会进 一步累积。因此，多视图重建的最后一步，通常是通过优化三维点的位置和照相机 参数来减少二次投影误差。该过程称为光束法平差。更多内容参阅参考文献 [13] 和 [35]， 或者可以从 http://en.wikipedia.org/wiki/Bundle_adjustment 参阅该方法的概述。

3. 自标定

在未标定照相机的情形中，有时可以从图像特征来计算标定矩阵。该过程称为自标 定。根据在照相机标定矩阵参数上做出的假设，以及可用的图像数据的类型(特征 匹配、平行线、平面等)，我们有很多不同的自标定算法。感兴趣的读者可以参阅参 考文献 [13] 及参考文献 [26] 第 6 章的内容。

作为标定的附注，值得一提的是一个叫做 extract_focal.pl 的有用脚本，它是 SfM 系 统(http://phototour.cs.washington.edu/bundler/)的一部分。对于常见的照相机，该 脚本基于图像 EXIF 数据，使用一个查找表来估计焦距。

4、立体图像

代码实现：

import stereo
import scipy.misc 
from PIL import Image
from pylab import *
#Fig5-10
# im_l = array(Image.open('../data/cones/im0.ppm').convert('L'),'f') 
# im_r = array(Image.open('../data/cones/im1.ppm').convert('L'),'f')
# #Fig5-11
im_l = array(Image.open('../data/cones/im0.ppm').convert('L'),'f') 
im_r = array(Image.open('../data/cones/im1.ppm').convert('L'),'f')

# 开始偏移，并设置步长 
steps = 12
start = 4 

# ncc 的宽度 
wid = 9

res1 = stereo.plane_sweep_ncc(im_l,im_r,start,steps,wid)

 
scipy.misc.imsave('../images/ch05/depth_ncc.png',res1)


res2 = stereo.plane_sweep_gauss(im_l,im_r,start,steps,wid)

 
scipy.misc.imsave('../images/ch05/depth_gauss.png',res2)

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname="/System/Library/Fonts/PingFang.ttc", size=14)
 
figure()

subplot(221)
imshow(im_l)
title(u'im_l', fontproperties=font)
axis('off')

subplot(222)
imshow(im_r)
title(u'im_r', fontproperties=font)
axis('off')

subplot(223)
imshow(res1)
title(u'ncc', fontproperties=font)
axis('off')

subplot(224)
imshow(res2)
title(u'gauss', fontproperties=font)
axis('off')

show()

运行结果：

5、其他

（1）开发环境

         本代码是在mac电脑sublimetxt编辑器python3.7.8下调试出来的，如果你要在windows/linux下编辑器/IDE的其他python版本运行的话，请对代码做相应的调整。

（2）源码和图片

        已经调试过的源码（含harris.py和sift.py文件）和图片详见：

https://github.com/Abonaventure/pcv-book-code.git   或  https://gitlab.com/Abonaventure/pcv-book-code.git

部分图片可能未上传可在https://pan.baidu.com/share/link?shareid=4059473268&uk=3440544532 或者原书目录http://programmingcomputervision.com下载