Python计算机视觉 第5章-多视图几何

5.1　外极几何

在计算机视觉的多视图几何中，外极几何（Epipolar Geometry）描述了两个相机视角之间的几何关系，常用于立体匹配、图像配准和3D重建。以下是主要概念：

1. 外极几何

• 描述两台摄像机之间的几何关系，与场景无关，只取决于摄像机的相对位置和内参。

2. 关键元素

• 极点（Epipole）：每个摄像机的光心在另一摄像机成像平面上的投影。
• 极线（Epipolar Line）：图像中的点在另一幅图像中的潜在匹配点所在的直线。

3. 基本矩阵（Fundamental Matrix）

• 描述两个图像之间点和极线的映射关系，定义了外极几何的基本约束。

4. 本质矩阵（Essential Matrix）

• 基本矩阵在已知相机内参时的特例，结合校准信息描述两相机之间的关系。

5. 应用场景

• 立体匹配：确定同一三维点在两幅图像中的匹配位置。
• 三维重建：使用极线约束来减少匹配搜索空间，重建场景的3D结构。
• 图像校准：利用外极几何推导摄像机的相对位置和旋转。

外极几何是理解多视图关系的核心概念，特别是在多视角协同应用中。

5.1.1　一个简单的数据集

在接下来的几节中，我们需要一个带有图像点、三维点和照相机参数矩阵的数据集。我们这里使用一个牛津多视图数据集；从http://www.robots.ox.ac.uk/~vgg/data/datamview.html 可以下载Merton1 数据的压缩文件。

（点进网址发现已经404了，这边直接把参考资料的代码贴上来）

下面的脚本可以加载Merton1的数据：

import camera

# 载入一些图像
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', missing='*')

# 载入照相机矩阵到 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()

实验结果如下所示：

图5-2 牛津multi-view数据集中的Merton1数据集：视图1与图像点（左）；视图1与投影的三维点（右）

5.1.2　用Matplotlib绘制三维数据

为了可视化三维重建结果，我们需要绘制出三维图像。Matplotlib中的mplot3d工具包可以方便地绘制出三维点、线、等轮廓线、表面以及其他基本图形组件，还可以通过图像窗口控件实现三维旋转和缩放。

可以通过在axes对象中加上projection="3d"关键字实现三维绘图，如下：

from matplotlib.pyplot import figure, show
from mpl_toolkits.mplot3d import axes3d

fig = figure()
ax = fig.add_subplot(111, projection="3d") 

# 生成三维样本点
X, Y, Z = axes3d.get_test_data(0.25)

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

show()

实验图1 绘图结果

get_test_data()函数在x, y空间按照设定的空间间隔参数来产生均匀的采样点。压平这些网格，会产生三列数据点，然后我们可以将其输入plot()函数。这样，我们就可以在立体表面上画出三维点。

现在通过画出Merton样本数据来观察三维点的效果：

# 绘制三维点
from matplotlib.pyplot import figure
from mpl_toolkits.mplot3d import axes3d

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

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

图5-3是三个不同视图中的三维图像点。图像窗口和控制界面外观效果像加上三维旋转工具的标准画图窗口。

图5-3：使用Matplottib工具包绘制的，牛津multi-view数据库中Mertoln1数据集的三维点：从上面和侧边观测的视图（左）；俯视的视图，展示了建筑墙体和屋顶上的点（中）；侧视图，展示了一面墙的轮廓，以及另一面墙上点的主视图（右）

5.1.3　计算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}{F}_{11}\\ F_{12}\\ F_{13}\\ ⋮\\ F_{33}\end{array}\right]=A\mathbf{f}=0$$

因为算法中需要8个对应点来计算基础矩阵F，所以该算法叫做八点法。最小化$||Af||$的函数代码如下：

def compute_fundamental(x1,x2):
    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")
    A = np.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)
    U,S,V = linalg.svd(F)
    S[2] = 0
    F = np.dot(U,np.dot(np.diag(S),V))
    return F

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

5.1.4　外极点和外极线

本节开始提到过，外极点满足$F{e}_1=0$，因此可以通过计算$F$的零空间来得到。把下面的函数添加代码中：

def compute_epipole(F):
    U,S,V = linalg.svd(F)
    e = V[-1]
    return e/e[2]

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

以及以下代码：

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

上面的函数将x轴的范围作为直线的参数，因此直线超出图像边界的部分会被截断。程序结果如图5-4所示：

图5-4：视图1中的外极线示为Merton1数据视图2中5个点对应的外极线。下面一行为这些点附近的特写。可以看到，这些线在图像外左侧位置将相交于一点。这些线表明点之间的对应可以在另外一幅图像中找到（线和点之间的颜色编码匹配）

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

5.2.1　三角剖分

三角剖分（Triangulation） 是计算机视觉和摄影测量中用于从多个视角恢复三维结构的关键技术。它基于两个或多个视角中的点对应关系，通过三角测量的方法来估计三维点的位置。以下是三角剖分的简要概念：

1. 基本定义

• 三角剖分 是从两个或多个图像中已知的二维点对（在不同视角中的点）来恢复三维空间中点的过程。它通过计算这些点在三维空间中的交点来实现。

2. 算法步骤

1. 已知对应点：在不同的视角中，首先需要确定同一三维点在不同图像中的投影位置（对应点）。

2. 相机模型：使用相机的内参和外参（即相机矩阵），描述每个相机的成像模型和它们之间的相对位置。

3. 三角化计算：基于已知的二维点对应关系和相机的内外参数，使用三角测量法计算三维点的坐标。通常，这通过解线性方程组来实现。

3. 应用场景

• 三维重建：从多个视角的图像中重建三维场景或对象。
• 立体视觉：在立体视觉系统中，计算和校准不同视角的图像来恢复场景的三维结构。
• 计算机视觉：用于机器人视觉、增强现实和虚拟现实中，帮助理解和分析三维世界。

4. 数学基础

• 相机矩阵：描述了相机的内外参数，将三维点投影到二维图像平面。
• 三角化方程：通过解两个相机视图的线性方程组来求解三维点的位置。

三角剖分是理解和重建三维结构的基础，广泛应用于计算机视觉、机器人和测量领域。

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

如果已经知道了一些三维点及其图像投影，我们可以使用直接线性变换的方法来计算照相机矩阵$P$。本质上，这是三角剖分方法的逆问题，有时我们将其称为照相机反切法。利用该方法恢复照相机矩阵同样也是一个最小二乘问题。

我们可以使用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个元素是照相机矩阵的元素，经过重新排列成矩阵形状后返回。

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

import sfm, camera
from numpy import where, vstack, ones

# 视图 1
corr = corr[:, 0]
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()

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

实验图2 控制台示例输出结果

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

图5-6：使用估计出的照相机矩阵来计算视图1中的投影点

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

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

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

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

因此，在无标定的情况下，第二个照相机矩阵可以使用一个（3×3）的射影变换得出。一个简单的方法是：
$$P_2=[S_{e}F|e]$$
其中，$e$是左极点，满足$e^{T}F=0$，$S_e$反对称矩阵。请注意，使用该矩阵做出的三角形剖分很有可能会发生畸变，如倾斜的重建。下面是具体的代码：

def compute_P_from_fundamental(F):
    """从基础矩阵中计算第二个照相机矩阵（假设 P1 = [I 0]）"""
    e = compute_epipole(F.T)  # 左极点
    Te = skew(e)
    return vstack((dot(Te, F.T).T, e)).T

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

在已经标定的情况下，重建会保持欧式空间中的一些度量特性（除了全局的尺度参数）。在重建三维场景中，已标定的例子更为有趣。

给定标定矩阵$K$，我们可以将它的逆$K^{-1}$作用于图像点$x_K=K^{-1}x$，因此，在新的图像坐标系下，照相机方程变为：
$$x_K=K^{-1}K[R|t]X=[R|t]X$$

在新的图像坐标系下，点同样满足之前的基础矩阵方程：
$$x_{K_2}^{T}F{x}_{K_1}=0$$

在标定归一化的坐标系下，基础矩阵称为本质矩阵。为了区别为标定后的情况，以及归一化了的图像坐标，我们通常将其记为E，而非F。
从本质矩阵恢复出的照相机矩阵中存在度量关系，但有四个可能解。因为只有一个解产生位于两个照相机前的场景，所以我们可以轻松地从中选出来。

下面是计算这4个解的算法，关键代码如下：

def compute_P_from_essential(E):
    U,S,V = np.svd(E)
    if np.det(np.dot(U,V))<0:
        V = -V
    E = np.dot(U,np.dot(np.diag([1,1,0]),V))
    Z = skew([0,0,-1])
    W = np.array([[0,-1,0],[1,0,0],[0,0,1]])
    P2 = [np.vstack((np.dot(U,np.dot(W,V)).T,U[:,2])).T,
            np.vstack((np.dot(U,np.dot(W,V)).T,-U[:,2])).T,
            np.vstack((np.dot(U,np.dot(W.T,V)).T,U[:,2])).T,
            np.vstack((np.dot(U,np.dot(W.T,V)).T,-U[:,2])).T]
    return P2

5.3　多视图重建

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

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

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

5.3.1　稳健估计基础矩阵

类似于稳健计算单应性矩阵，当存在噪声和不正确的匹配时，我们需要估计基础矩阵。其实现代码为：

class RansacModel(object):
    """用从 http://www.scipy.org/Cookbook/RANSAC 下载的 ransac.py 计算基础矩阵的类"""

    def __init__(self, debug=False):
        self.debug = debug

    def fit(self, data):
        """使用选择的 8 个对应计算基础矩阵"""
        # 转置，并将数据分成两个点集
        data = data.T
        x1 = data[:3, :8]
        x2 = data[3:, :8]

        # 估计基础矩阵，并返回
        F = compute_fundamental_normalized(x1, x2)
        return F

    def get_error(self, data, F):
        """计算所有对应的 x^T F x，并返回每个变换后点的误差"""
        # 转置，并将数据分成两个点集
        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()和get_error()方法。这里采用的错误衡量方法是Sampson 距离。 fit()方法会选择8个点，然后使用归一化的八点算法：

def compute_fundamental_normalized(x1, x2):
    """使用归一化的八点算法，由对应点（x1，x2 3×n 的数组）计算基础矩阵"""
    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))
    return F / F[2, 2]

该函数将图像点归一化为零均值固定方差。

接下来，我们在函数中使用该类。

def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
    """使用 RANSAC 方法（ransac.py，来自 http://www.scipy.org/Cookbook/RANSAC），
    从点对应中稳健地估计基础矩阵 F
    输入：使用齐次坐标表示的点 x1，x2（3×n 的数组）"""
    import ransac
    
    data = vstack((x1, x2))
    
    # 计算 F，并返回正确点索引
    F, ransac_data = ransac.ransac(data.T, model, 8, maxiter, match_threshold, 20, return_all=True)
    
    return F, ransac_data['inliers']

这里，我们返回最佳基础矩阵F，以及正确点的索引，这样可以知道哪些匹配和F矩阵是一致的。与单应性矩阵估计相比，我们增加了默认的最大迭代次数，改变了匹配的阈值（之前使用像素，现在使用Sampson距离来衡量）。

5.4　立体图像

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

通过将图像扭曲到公共的平面上，使外极线位于图像行上，任何立体照相机设置都能得到矫正（我们通常构建立体平台来产生经过矫正的图像对）。

假设两幅图像经过了矫正，那么对应点的寻找限制在图像的同一行上。一旦找到对应点，由于深度是和偏移成正比的，那么深度（Z坐标）可以直接由水平偏移来计算：

$$Z=\frac{fb}{x_l-x_r}$$
其中，$f$是经过矫正图像的焦距，$b$是两个照相机中心之间的距离，$x_l$和$x_r$是左右两幅图像中对应点的x坐标。分开照相机中心的距离称为基线。矫正后的立体照相机设置如图5-9所示。

图5-9：矫正后立体照相机设置的示意图，其中对应点位于两幅图像的同一行

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

计算视差图

在该立体重建算法中，我们将对于每个像素尝试不同的偏移，并按照局部图像周围归一化的互相关值，选择具有最好分数的偏移，然后记录下该最佳偏移。因为每个偏移在某种程度上对应于一个平面，所以该过程有时称为扫平面法。虽然该方法并不是立体重建中最好的方法，但是非常简单，通常会得出令人满意的结果。

关键代码如下：

def plane_sweep_ncc(im_l, im_r, start, steps, wid):
    """使用归一化的互相关计算视差图像"""
    m, n = im_l.shape
    
    # 保存不同求和值的数组
    mean_l = zeros((m, n))
    mean_r = zeros((m, n))
    s = zeros((m, n))
    s_l = zeros((m, n))
    s_r = zeros((m, n))
    
    # 保存深度平面的数组
    dmaps = zeros((m, n, steps))
    
    # 计算图像块的平均值
    filters.uniform_filter(im_l, wid, mean_l)
    filters.uniform_filter(im_r, wid, mean_r)
    
    # 归一化图像
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    
    # 尝试不同的视差
    for displ in range(steps):
        # 将左边图像移动到右边，计算加和
        filters.uniform_filter(roll(norm_l, -displ - start) * norm_r, wid, s)  # 和归一化
        filters.uniform_filter(roll(norm_l, -displ - start) * roll(norm_l, -displ - start), wid, s_l)
        filters.uniform_filter(norm_r * norm_r, wid, s_r)  # 和反归一化
        
        # 保存 NCC 的分数
        dmaps[:, :, displ] = s / sqrt(s_l * s_r)
    
    # 为每个像素选取最佳深度
    return argmax(dmaps, axis=2)

import stereo
from numpy import array
from PIL import Image
import scipy.misc

# 读取图像并转换为灰度
im_l = array(Image.open('scene1.row3.col3.ppm').convert('L'), 'f')
im_r = array(Image.open('scene1.row3.col4.ppm').convert('L'), 'f')

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

# NCC 的宽度
wid = 9

# 计算视差图像
res = stereo.plane_sweep_ncc(im_l, im_r, start, steps, wid)

# 保存视差图像
scipy.misc.imsave('depth.png', res)

def plane_sweep_gauss(im_l, im_r, start, steps, wid):
    """使用带有高斯加权周边的归一化互相关计算视差图像"""
    m, n = im_l.shape
    
    # 保存不同加和的数组
    mean_l = zeros((m, n))
    mean_r = zeros((m, n))
    s = zeros((m, n))
    s_l = zeros((m, n))
    s_r = zeros((m, n))
    
    # 保存深度平面的数组
    dmaps = zeros((m, n, steps))
    
    # 计算平均值
    filters.gaussian_filter(im_l, wid, 0, mean_l)
    filters.gaussian_filter(im_r, wid, 0, mean_r)
    
    # 归一化图像
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    
    # 尝试不同的视差
    for displ in range(steps):
        # 将左边图像移动到右边，计算加和
        filters.gaussian_filter(roll(norm_l, -displ - start) * norm_r, wid, 0, s)  # 和归一化
        filters.gaussian_filter(roll(norm_l, -displ - start) * roll(norm_l, -displ - start), wid, 0, s_l)
        filters.gaussian_filter(norm_r * norm_r, wid, 0, s_r)  # 和反归一化
        
        # 保存 NCC 的分数
        dmaps[:, :, displ] = s / sqrt(s_l * s_r)
    
    # 为每个像素选取最佳深度
    return argmax(dmaps, axis=2)

可以像前面的扫平面函数一样来使用以上函数。图5-10和图5-11所示为这两个扫平面实现操作在一些标准立体基准图像上的结果。

图5-10：使用归一化互相关从一对立体图像中计算视差图

图5-11：使用归一化互相关从一对立体图像中计算视差图