《Python 计算机视觉编程》学习笔记（五）

第 5 章 多视图几何

引言

本章讲解如何处理多个视图，以及如何利用多个视图的几何关系来恢复照相机位置信息和三维结构。通过在不同视点拍摄的图像，我们可以利用特征匹配来计算出三维场景点以及照相机位置。

5.1　外极几何

多视图几何是利用在不同视点所拍摄图像间的关系，来研究照相机之间或者特征之间关系的一门科学。图像的特征通常是兴趣点，本章使用的也是兴趣点特征。多视图几何中最重要的内容是双视图几何。

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

没有关于照相机的先验知识，会出现固有二义性，因为三维场景点 X 经过 4× 4 的单应性矩阵 H 变换为 HX 后，则 HX 在照相机 PH -1 里得到的图像点和 X 在照相机P 里得到的图像点相同。利用照相机方程，可以将上述问题描述为：
$$\lambda \mathbf{x}=\mathbf{PX}=\mathbf{PH}{}^{-1}\mathbf{HX}=\hat{\mathbf{P}}\hat{\mathbf{X}}$$

因此，当我们分析双视图几何关系时，总是可以将照相机间的相对位置关系用单应性矩阵加以简化。这里的单应性矩阵通常只做刚体变换，即只通过单应矩阵变换了坐标系。一个比较好的做法是，将原点和坐标轴与第一个照相机对齐，则：$${\mathbf{P}}_1={\mathbf{K}}_1[\mathbf{I}\mid 0]\text{ 和 }{\mathbf{P}}_2={\mathbf{K}}_2[\mathbf{R}\mid \mathbf{t}]$$

其中 K1 和 K2 是标定矩阵， R是第二个照相机的旋转矩阵， t 是第二个照相机的平移向量。利用这些照相机参数矩阵，我们可以找到点 X 的投影点 x1 和 x2（分别对应于投影矩阵 P1 和 P2）的关系。这样，我们可以从寻找对应的图像出发，恢复照相机参数矩阵。

同一个图像点经过不同的投影矩阵产生的不同投影点必须满足：
$${\mathbf{x}}_2^T\mathbf{F}{\mathbf{x}}_1=0$$

其中：$$\mathbf{F}={\mathbf{K}}_2^{-T}{\mathbf{S}}_{\mathrm{t}}\mathbf{R}{\mathbf{K}}_1^{-1}$$
矩阵$S_t$为反对称矩阵：
$${\mathbf{S}}_{\mathbf{t}}=\left[\begin{array}{rrr}0& -t_3& t_2\\ t_3& 0& -t_1\\ -t_2& t_1& 0\end{array}\right]$$

公式${\mathbf{x}}_2^T\mathbf{F}{\mathbf{x}}_1=0$为外极约束条件。矩阵 F 为基础矩阵。 基础矩阵可以由两照相机的参数矩阵（相对旋转 R 和平移 t）表示。由于 det(F)=0，所以基础矩阵的秩小于等于 2。

上面的公式表明，我们可以借助 F 来恢复出照相机参数，而 F 可以从对应的投影图像点计算出来。在不知道照相机内参数（ K1 和 K2）的情况下，仅能恢复照相机的投影变换矩阵。在知道照相机内参数的情况下，重建是基于度量的。所谓度量重建，即能够在三维重建中正确地表示距离和角度 。

利用上面的理论处理一些图像数据前，我们还需要了解一些几何学知识。

给定图像中的一个点，例如第二个视图中的图像点 x2，利用公式${\mathbf{x}}_2^T\mathbf{F}{\mathbf{x}}_1=0$，我们找到对应第一幅图像的一条直线：$${\mathbf{x}}_2^T\mathbf{F}{\mathbf{x}}_1={\mathbf{l}}_1^T{\mathbf{x}}_1=0$$

其中${\mathbf{I}}_1^T{\mathbf{x}}_1=0$是第一幅图像中的一条直线。这条线称为对应于点 x2 的外极线，也就是说，点 x2 在第一幅图像中的对应点一定在这条线上。所以，基础矩阵可以将对应点的搜索限制在外极线上。外极线都经过一点 e，称为外极点。实际上，外极点是另一个照相机光心对应的图像点。外极点可以在我们看到的图像外，这取决于照相机间的相对方向。因为外极点在所有的外极线上，所以$F_{e_1}=0$。因此，我们可以通过计算 F 的零向量来计算出外极点。同理，另一个外极点可以通过计算$e_2^{T}F=0$得到。外极点和外极线如图 5-1 所示。

加载带有图像点、三维点和照相机参数矩阵的数据集。

from pylab import *
from PCV.geometry import camera
from PIL import Image
import matplotlib.font_manager as fm
import time


#定义字体模板
myfont=fm.FontProperties(fname='C:/Windows/Fonts/simsun.ttc') 


def loadVggdata(imagepath):
    """ 函数功能：加载图像 
        参数说明：imagepath__图像数据
        函数返回：无
    """
    # 载入一些图像
    im1 = array(Image.open(imagepath[0]))
    im2 = array(Image.open(imagepath[1]))
    # 载入每个视图的二维点到列表中
    points2D = [loadtxt(imagepath[4]+'00'+str(i+1)+'.corners').T for i in range(3)]
    # 载入三维点
    points3D = loadtxt(imagepath[2]).T

    # 载入对应
    corr = genfromtxt(imagepath[3])

    # 载入照相机矩阵到 Camera 对象列表中
    P = [camera.Camera(loadtxt(imagepath[4]+'00'+str(i+1)+'.P')) for i in range(3)]

    # 将三维点转换成齐次坐标表示，并投影
    X = vstack( (points3D, ones(points3D.shape[1])))
    x = P[0].project(X)

    # 在视图1中绘制点
    figure()
    subplot(121)
    imshow(im1),plot(points2D[0][0], points2D[0][1],'*'), axis('off'),title('视图 1 与图像点',fontproperties = myfont)
    subplot(122)
    imshow(im1),plot(x[0],x[1],'r.'), axis('off'),title('视图 1 与投影的三维点',fontproperties = myfont)
    savefig('./output/test1.jpg')
    show()

    return


def main():
    imagepath = [r'C:/hqq/document/python/computervision/ch05/MertonCollege1/image/001.jpg',
                r'C:/hqq/document/python/computervision/ch05/MertonCollege1/image/002.jpg',
                r'C:/hqq/document/python/computervision/ch05/MertonCollege1/3d/p3d',
                r'C:/hqq/document/python/computervision/ch05/MertonCollege1/2d/nview-corners',
                r'C:/hqq/document/python/computervision/ch05/MertonCollege1/2d/',
                ]

    loadVggdata(imagepath)

    return

if __name__ =="__main__":

    startTime = time.time()
    main()
    endTime = time.time()
    print("the program running time is :",endTime - startTime)

输出结果如下：

左图为第一个视图以及该视图中的图像点。投影后的点绘制在另一张图上，如图所示。第二幅图比第一幅图多一些点。这些多出的点是从视图 2 和视图 3 重建出来的，而不在视图 1 中。

用Matplotlib绘制三维数据

from mpl_toolkits.mplot3d import axes3d
from pylab import *

def plot3Dimage():
    """ 函数功能：绘制3D图像 
        参数说明：imagepath__图像路径
        函数返回：无
    """

    fig = figure()
    ax = fig.add_subplot(projection='3d')

    # 生成三维样本点
    # get_test_data() 函数在 x, y 空间按照设定的空间间隔参数来产生均匀的采样点。
    X, Y, Z = axes3d.get_test_data(0.25)

    # 在三维中绘制点
    ax.plot(X.flatten(),Y.flatten(),Z.flatten(),'o')
    show()

输出如下:

上图为随机输出的三位点。
画出 Merton 样本数据来观察三维点：

def plot3Dimage(imagepath):
    """ 函数功能：绘制3D图像 
        参数说明：imagepath__图像路径
        函数返回：无
    """

    fig = figure()
    ax = fig.add_subplot(projection='3d')

    # 生成三维样本点
    # get_test_data() 函数在 x, y 空间按照设定的空间间隔参数来产生均匀的采样点。
    # X, Y, Z = axes3d.get_test_data(0.25)

    # 在三维中绘制点
    # ax.plot(X.flatten(),Y.flatten(),Z.flatten(),'o')
    # 载入三维点
    points3D = loadtxt(imagepath[2]).T

    ax.plot(points3D[0],points3D[1],points3D[2],'k.')
    show()

    return

输出的3D图像如上图所示。



上面三幅图像分别为从上面和侧边观测的视图（第一幅）；俯视的视图，展示了建筑墙体和屋顶上的点（第二幅）；侧视图，展示了一面墙的轮廓，以及另一面墙上点的主视图（最后一幅）。

计算F： 八点法

八点法是通过对应点来计算基础矩阵的算法。
外极约束可以写成线性系统的形式：
$$\left[\begin{array}{ccccc}{x}_2^1{x}_1^1& x_2^1{y}_1^1& x_2^1{w}_1^1& \cdots & w_2^1{w}_1^1\\ x_2^2{x}_1^2& x_2^2{y}_1^2& x_2^2{w}_1^2& \cdots & w_2^2{w}_1^2\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_2^n{x}_1^n& x_2^n{y}_1^n& x_2^n{w}_1^n& \cdots & w_2^n{w}_1^n\end{array}\right]\left[\begin{array}{c}{\mathbf{F}}_{11}\\ {\mathbf{F}}_{12}\\ {\mathbf{F}}_{13}\\ ⋮\\ {\mathbf{F}}_{33}\end{array}\right]=\mathbf{Af}=0$$
f 包含 F 的元素，${\mathbf{x}}_1^i=\left[x_1^i,y_1^i,w_1^i\right]\text{ 和 }{\mathbf{x}}_2^i=\left[x_2^i,y_2^i,w_2^i\right]$是一对图像点，共有 n 对对应点。基础矩阵中有 9 个元素，由于尺度是任意的，所以只需要 8 个方程。因为算法中需要 8 个对应点来计算基础矩阵 F，所以该算法叫做八点法。

八点算法的基本步骤：
1）求线性解 由系数矩阵A最小奇异值对应的奇异向量f求的F。
2）奇异性约束 是最小化Frobenius范数‖ F − F ′ ‖中的 F′ 代替F

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

from scipy import *
from numpy 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

我们通常用 SVD 算法来计算最小二乘解。由于上面算法得出的解可能秩不为 2（基础矩阵的秩小于等于 2），所以我们通过将最后一个奇异值置 0 来得到秩最接近 2 的基础矩阵。这是个很有用的技巧。上面的函数忽略了一个重要的步骤：对图像坐标进行归一化，这可能会带来数值问题。
计算基础矩阵步骤如下：
1）sift提取特征
2）RANSAC去除错误点匹配
3）归一化8点算法估计基础矩阵

sift结果如下：

使用ransac去除错误点后：

基础矩阵输出如下：

外极点和外极线

本节开始提到过，外极点满足 Fe1=0，因此可以通过计算 F 的零空间来得到。把下面的函数添加到 sfm.py 中

def compute_epipole(F):
    """从基础矩阵F中计算右极点（可以使用F.T获得左极点）"""

    # 返回F的零空间（Fx=0）
    U,S,V = linalg.svd(F)
    e = V[-1]
    return e/e[2]

如果想获得另一幅图像的外极点（对应左零空间的外极点），只需将 F 转置后输入上述函数即可。

def plotpoint(imagepath):
        # 载入一些图像
    im1 = array(Image.open(imagepath[7]))
    im2 = array(Image.open(imagepath[8]))

    # 载入每个视图的二维点到列表中
    points2D = [loadtxt(imagepath[4]+'/00'+str(i+1)+'.corners').T for i in range(3)]

    # 载入三维点
    points3D = loadtxt(imagepath[2]).T

    # 载入对应
    corr = genfromtxt(imagepath[3],dtype='int',missing_values='*')

    # 载入照相机矩阵到 Camera 对象列表中
    P = [camera.Camera(loadtxt(imagepath[4]+'/00'+str(i+1)+'.P')) for i in range(3)]

    # 将三维点转换成齐次坐标表示，并投影
    X = vstack( (points3D, ones(points3D.shape[1])))
    x = P[0].project(X)
    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):
        sfm.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()

    return

输出图像如下:


上图为5个不同颜色的点和5条不同颜色的线。这些线表明点之间的对应可以在另外一幅图像中找到（线和点之间的颜色编码匹配）两幅图片还是有差别的。

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

给定照相机参数模型，图像点可以通过三角剖分来恢复出这些点的三维位置。基本的算法思想如下。

对于两个照相机 P1 和 P2 的视图，三维实物点 X 的投影点为 x1 和 x2（这里用齐次坐标表示），照相机方程（ 4.1）定义了下列关系：
$$\left[\begin{array}{ccc}{P}_1& -{\mathbf{x}}_1& 0\\ P_2& 0& -{\mathbf{x}}_2\end{array}\right]\left[\begin{array}{l}\mathbf{X}\\ {\lambda }_1\\ {\lambda }_2\end{array}\right]=0$$

由于图像噪声、照相机参数误差和其他系统误差，上面的方程可能没有精确解。我们可以通过 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]

多个点的三角剖分的函数来实现:

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数据集上的三角剖分代码实现:

# 获得前两个视图中点的索引
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])

# 绘制图像
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()

输出图片如下：

上面的代码首先利用前两个视图的信息来对图像点进行三角剖分，然后把前三个图像点的齐次坐标输出到控制台，最后绘制出恢复的最接近三维图像点。输出到控制台的信息如下:

算法估计出的三维图像点和实际图像点很接近。如图所示，估计点和实际点可以很好地匹配

红色点位实际点，黑色点位估计点。可以看到大致上两个颜色的点都有重合。

由三维点计算照相机矩阵

如果已经知道了一些三维点及其图像投影，我们可以使用直接线性变换的方法来计算照相机矩阵 P。本质上，这是三角剖分方法的逆问题，有时我们将其称为照相机反切法。利用该方法恢复照相机矩阵同样也是一个最小二乘问题。我们从照相机方程（ 4.1）可以得出，每个三维点 Xi（齐次坐标系下）按照${\lambda }_i{\mathbf{x}}_i={\mathbf{PX}}_i$投影到图像点${\mathbf{x}}_i=\left[x_i,y_i,1\right]$，相应的点满足下面的关系：
$$\left[\begin{array}{ccccccc}{\mathbf{X}}_1^T& 0& 0& -x_1& 0& 0& \dots \\ 0& {\mathbf{X}}_1^T& 0& -y_1& 0& 0& \dots \\ 0& 0& {\mathbf{X}}_1^T& -1& 0& 0& \dots \\ {\mathbf{X}}_2^T& 0& 0& 0& -x_2& 0& \dots \\ 0& {\mathbf{X}}_2^T& 0& 0& -y_2& 0& \dots \\ 0& 0& {\mathbf{X}}_2^T& 0& -1& 0& \dots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮‘在这里插入代码片‘\end{array}\right]\left[\begin{array}{c}{\mathbf{p}}_1^T\\ {\mathbf{p}}_2^T\\ {\mathbf{p}}_3^T\\ {\lambda }_1\\ {\lambda }_2\\ ⋮\end{array}\right]=0$$

p1、 p2 和 p3 是矩阵 P 的三行。上面的式子可以写得更紧凑，如下所示：$M\mathbf{v}=0$

使用 SVD 分解估计出照相机矩阵实现代码如下：

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

该函数的输入参数为图像点和三维点，构造出上述所示的 M 矩阵。最后一个特征向量的前 12 个元素是照相机矩阵的元素，经过重新排列成矩阵形状后返回。

选出第一个视图中的一些可见点（使用对应列表中缺失的数值），将它们转换为齐次坐标表示，然后估计照相机矩阵：

from pylab import *
from PIL import Image
from mpl_toolkits.mplot3d import axes3d
from PCV.geometry import camera
from PCV.geometry import sfm
from PCV.localdescriptors import sift
# 载入一些图像
im1 = array(Image.open('D:\\123\图像处理\Image Processing\Image Processing\Chapter5\\001.jpg'))
im2 = array(Image.open('D:\\123\图像处理\Image Processing\Image Processing\Chapter5\\002.jpg'))

# 载入每个视图的二维点到列表中
points2D = [loadtxt('D:\\123\图像处理\Image Processing\Image Processing\Chapter5\\00'+str(i+1)+'.corners').T for i in range(3)]

# 载入三维点
points3D = loadtxt('D:\\123\图像处理\Image Processing\Image Processing\Chapter5\p3d').T

# 载入对应
corr = genfromtxt('D:\\123\图像处理\Image Processing\Image Processing\Chapter5\\nview-corners',dtype='int',missing_values='*')

# 载入照相机矩阵到 Camera 对象列表中
P = [camera.Camera(loadtxt('D:\\123\图像处理\Image Processing\Image Processing\Chapter5\\00'+str(i+1)+'.P')) for i in range(3)]

# 将三维点转换成齐次坐标表示，并投影
X = vstack( (points3D, ones(points3D.shape[1])))
x = P[0].project(X)

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

为了检查照相机矩阵的正确性，将它们以归一化的格式（除以最后一个元素）打印到控制台。输出如下所示：

上面是估计出的照相机矩阵，下面是数据集的创建者计算出的照相机矩阵。可以看到，它们的元素几乎完全相同。最后，使用估计出的照相机矩阵投影这些三维点，然后绘制出投影后的结果。真实点用圆圈表示，估计出的照相机投影点用点表示。

可以看到基本重合在了一起，说明估计的较为准确。

由基础矩阵计算照相机矩阵

在两个视图的场景中，照相机矩阵可以由基础矩阵恢复出来。假设第一个照相机矩阵归一化为 P1=[I|0]，现在我们需要计算出第二个照相机矩阵 P2。研究分为两类，未标定的情况和已标定的情况。

1.未标定的情况——投影重建

在没有任何照相机内参数知识的情况下，照相机矩阵只能通过射影变换恢复出来。也就是说，如果利用照相机的信息来重建三维点，那么该重建只能由射影变换计算出来（你可以得到整个投影场景中无畸变的重建点）。在这里，我们不考虑角度和距离。

因此，在无标定的情况下，第二个照相机矩阵可以使用一个（ 3× 3）的射影变换得出。一个简单的方法是：
$${\mathbf{P}}_2=\left[{\mathbf{S}}_{\mathbf{e}}\mathbf{F}\mid \mathbf{e}\right]$$
e 是左极点，满足 eTF=0， Se 是如公式（ 5.2）所示的反对称矩阵。请注意，使用该矩阵做出的三角形剖分很有可能会发生畸变，如倾斜的重建。

2.已标定的情况——度量重建

在已经标定的情况下，重建会保持欧式空间中的一些度量特性（除了全局的尺度参数）。在重建三维场景中，已标定的例子更为有趣。给定标定矩阵 K，我们可以将它的逆 K -1 作用于图像点 xK=K -1x，因此，在新的图像坐标系下，照相机方程变为：
$${\mathbf{x}}_K=K^{-1}K[R\mid \mathbf{t}]\mathbf{X}=[R\mid \mathbf{t}]\mathbf{X}$$

在新的图像坐标系下，点同样满足之前的基础矩阵方程：
$${\mathbf{x}}_{K_2}^T\mathbf{F}{\mathbf{x}}_{K_1}=0$$
在标定归一化的坐标系下，基础矩阵称为本质矩阵。为了区别为标定后的情况，以及归一化了的图像坐标，我们通常将其记为 E，而非 F。从本质矩阵恢复出的照相机矩阵中存在度量关系，但有四个可能解。因为只有一个解产生位于两个照相机前的场景，所以我们可以轻松地从中选出来。

5.3 多视图重建

如何使用上面的理论从多幅图像中计算出真实的三维重建。由于照相机的运动给我们提供了三维结构，所以这样计算三维重建的方法通常称为 SfM（ Structure from Motion，从运动中恢复结构）。

假设照相机已经标定，计算重建可以分为下面 4 个步骤:

1. 检测特征点，然后在两幅图像间匹配；
2. 由匹配计算基础矩阵；
3. 由基础矩阵计算照相机矩阵；
4. 三角剖分这些三维点。

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

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

与单应性矩阵估计相比，我们增加了默认的最大迭代次数，改变了匹配的阈值（之前使用像素，现在使用 Sampson 距离来衡量）。

三维重建示例:

def rebuild3D(imagepath):

    # 标定矩阵
    K = array([[2394,0,932],[0,2398,628],[0,0,1]])
    # 载入图像，并计算特征


    im1 = array(Image.open('C:/hqq/document/python/computervision/ch05/alcatraz1.jpg'))
    sift.process_image('C:/hqq/document/python/computervision/ch05/alcatraz1.jpg', 'im0.sift')
    l1, d1 = sift.read_features_from_file('im0.sift')

    im2 = array(Image.open('C:/hqq/document/python/computervision/ch05/alcatraz2.jpg'))
    sift.process_image('C:/hqq/document/python/computervision/ch05/alcatraz2.jpg', 'im1.sift')
    l2, d2 = sift.read_features_from_file('im1.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]

    # 绘制三维图像
    fig = figure()
    #ax = fig.gca(projection='3d')
    ax = fig.add_subplot(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()
    
    return

输出结果如下：

可以看到，二次投影后的点和原始特征位置不完全匹配，但是相当接近。

多视图的扩展示例

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

对于视频序列，我们可以使用时序关系，在连续的帧对中匹配特征。相对方位需要从每个帧对增量地加入下一个帧对。该方法非常有效，同时可以使用跟踪有效地找到对应点。一个问题是误差会在加入的视图中积累。

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

2.光束法平差
从简单三维重建的例子，我们可以清楚地看出，恢复出的点的位置和由估计的基础矩阵计算出的照相机矩阵都存在误差。在多个视图的计算中，这些误差会进一步累积。因此，多视图重建的最后一步，通常是通过优化三维点的位置和照相机参数来减少二次投影误差。该过程称为光束法平差。

3.自标定
在未标定照相机的情形中，有时可以从图像特征来计算标定矩阵。该过程称为自标定。根据在照相机标定矩阵参数上做出的假设，以及可用的图像数据的类型（特征匹配、平行线、平面等），我们有很多不同的自标定算法。

5.4　立体图像

一个多视图成像的特殊例子是立体视觉（或者立体成像），即使用两台只有水平（向一侧）偏移的照相机观测同一场景。当照相机的位置如上设置，两幅图像具有相同的图像平面，图像的行是垂直对齐的，那么称图像对是经过矫正的。该设置在机器人学中很常见，常被称为立体平台。

通过将图像扭曲到公共的平面上，使外极线位于图像行上，任何立体照相机设置都能得到矫正（我们通常构建立体平台来产生经过矫正的图像对）。假设两幅图像经过了矫正，那么对应点的寻找限制在图像的同一行上。一旦找到
对应点，由于深度是和偏移成正比的，那么深度（ Z 坐标）可以直接由水平偏移来计算，
$$Z=\frac{fb}{x_l-x_r}$$
f 是经过矫正图像的焦距， b 是两个照相机中心之间的距离，$x_l和x_r$ 是左右两幅图像中对应点的 x 坐标。分开照相机中心的距离称为基线。矫正后的立体照相机设置如图 5-9 所示。

立体重建（有时称为致密深度重建）就是恢复深度图（或者相反，视差图），图像中每个像素的深度（或者视差）都需要计算出来。这是计算机视觉中的经典问题，有很多算法可以解决该问题。

5.5 小结

本章主要是介绍了利用多个视图的几何关系来恢复照相机位置信息和三维结构。充分利用照相机方程和双试图的几何关系，将照相机的相对位置关系用单应性矩阵简化。充分利用几何学知识将我们的基础矩阵求解出来。之后便是三位点和照相机矩阵、基础矩阵之间的转换。

忘了，代码和图片，需要自取：
链接：https://pan.baidu.com/s/1wetY1lJcZmbv4gntz7YwlA?pwd=u2uj
提取码：u2uj