Python计算机视觉编程——第五章 多视图几何

1 外极几何

多视图几何是利用在不同视点所拍摄图像间的关系，来研究照相机之间或者特征之间关系的一门科学。最重要的内容是双视图几何。如果有一个场景的两个视图以及视图中的对应图像点，那么根据照相机间的空间相对位置关系、照相机的性质以及三维场景点的位置，可以得到对这些图像点的一些几何关系约束。我们通过外极几何来描述这些几何关系。

同一图像点经过不同的投影矩阵产生不同投影点必须满足：
$${\mathbf{x}}_2^{T}F{\mathbf{x}}_1=0$$
其中：
$$F=K_2^{-T}{S}_{t}R{K}_1^{-1}$$
矩阵$S_{\mathrm{t}}$为反对称矩阵：
$${\mathbf{S}}_1=\left[\begin{array}{ccc}0& -t_3& t_2\\ t_3& 0& -t_1\\ -t_2& t_1& 0\end{array}\right]$$
矩阵F为基础矩阵。基础矩阵可以由两照相机的参数矩阵。
可以借助F来恢复出照相机参数，F可以从对应的投影图像点计算出来。在不知道照相机内参数的情况下，仅能恢复照相机的投影变换矩阵。在知道照相机内参数的情况下，重建是基于度量的。
基础矩阵可以将对应点的搜索限制在外极线上。外极线经过点e，称为外极点。

1.1 简单数据集

从相应网址上下载一个牛津多视图数据集，加载Merton1数据并且可视化数据，将三维的点投影到一个视图。

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

im1 = array(Image.open('001.jpg'))
im2 = array(Image.open('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_values='*')
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)


figure()
imshow(im1)
plot(points2D[0][0],points2D[0][1],'*')
axis('off')

figure()
imshow(im1)
plot(x[0],x[1],'r.')
axis('off')
show()

处理结果如图：

1.2 用Matplpotlib绘制三维数据

Matplotlib中的mplot3d工具包可以方便地绘制出三维点，线，等轮廓线，表面以及其他基本图形组件，还可以通过图像窗口控件实现三维旋转和缩放。
通过再axes对象中加上projection="3d"关键字实现三维绘图，代码如下：

import camera
import numpy as np
from PIL import Image
from pylab import *
from mpl_toolkits.mplot3d import axes3d

im1 = array(Image.open('001.jpg'))
im2 = array(Image.open('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_values='*')
P = [camera.Camera(loadtxt('2D/00'+str(i+1)+'.P')) for i in range(3)]

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

X,Y,Z = axes3d.get_test_data(0.25)

ax.plot(X.flatten(),Y.flatten(),Z.flatten(),'o')
show()

get_test_data()函数再x，y空间按照设定的空间间隔参数来产生均匀的采样点。

得到结果
注意：出现错误，FigureBase.gca() got an unexpected keyword argument ‘projection’。解决方法：使用 figure() 函数来创建一个新的figure，然后使用 add_subplot() 或者 add_axes() 函数来添加一个Axes。

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$$
其中，$f$包含 $F$的元素，$x_1^i=\left[x_1^i,y_1^i,w_1^i\right]$和${\mathbf{x}}_2^i=\left[x_2^i,y_2^i,w_2^i\right]$是一对图像点，共有$n$ 对对应点。基础矩阵中有 9 个元素，由于尺度是任意的，所以只需要 8 个方程。因为算法中需要 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 = 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 = dot(U, dot(diag(S), V))
    return F

通常使用SVD算法来计算最小二乘解。由于上面算法得出的解可能秩不为2，所以通过将最后一个奇异值置0来得到秩最接近2的基础矩阵。

1.4 外极点和外极线

外极点满足$F{e}_1=0$，可以通过计算F的零空间来得到。将函数添加到sfm.py中。

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

在之前样本数据集运行代码如下：

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

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

在第一个视图中画出了前五个外极线，在第二个视图中画出来对应匹配点。

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

2.1 三角剖分

给定照相机参数模型，图像点可以通过三角剖分来恢复出这些点的三维位置。思想如下：
对于两个照相机$P_1$和$P_2$的视图，三维实物点$\mathbf{X}$的投影点为${\mathbf{x}}_1$和${\mathbf{x}}_2$ (这里用齐次坐
标表示),照相机方程定义了下列关系：
$$\left[\begin{array}{ccc}{P}_1& -{\mathbf{x}}_1& 0\\ P_2& 0& -{\mathbf{x}}_2\end{array}\right]\left[\begin{array}{c}\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):
    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])

fig = figure()
ax = fig.add_subplot(projection = '3d')
ax.plot(Xest[0],Xest[1],Xest[2],'ko')
ax.plot(Xtrue[0],Xtrue[1],Xtrue[2],'r.')
axis('equal')
show()

首先利用前两个视图的信息来对图像点进行三角剖分，将前三个图像点的齐次坐标输出到控制台，绘制出恢复的最接近三维图像点。
打印和绘制信息如下：

估计点和实际点可以很好匹配。

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

当知道一些三维点及其图像投影，可以使用直接线性变换来计算照相机矩阵P。称为照相机反切法。利用该方法恢复照相机矩阵同样也是最小二乘问题。每个三维点${\mathbf{X}}_i$(齐次坐标系下)按照${\lambda }_i{\mathbf{x}}_i=\mathbf{P}{\mathbf{X}}_i$
投影到图像点${\mathbf{x}}_i=[x_i,y_i,1]$,相应的点满足下面的关系：
$$\left[\begin{array}{ccccccc}{\mathbf{X}}_1^T& 0& 0& -x_1& 0& 0& \cdots \\ 0& {\mathbf{X}}_1^T& 0& -y_1& 0& 0& \cdots \\ 0& 0& {\mathbf{X}}_1^T& -1& 0& 0& \cdots \\ {\mathbf{X}}_2^T& 0& 0& 0& -x_2& 0& \cdots \\ 0& {\mathbf{X}}_2^T& 0& 0& -y_2& 0& \cdots \\ 0& 0& {\mathbf{X}}_2^T& 0& -1& 0& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\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$$
其中${\mathbf{p}}_1,{\mathbf{p}}_2$和${\mathbf{p}}_3$是矩阵$\mathbf{P}$的三行。
可以使用SVD分解估计出照相机矩阵。

将下面函数添加到sfm.py文件中(输入参数为图像点和三维点)：

def compute_P(x, X):
    n = x.shape[1]
    if X.shape[1] != n:
        raise ValueError("Number of points don't match.")

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

在样本数据集上测试算法性能，选出第一个视图中的一些可见点，将它们转换为齐次坐标表示，然后估计照相机矩阵：

import sfm
import camera
import numpy as np
from PIL import Image
from pylab import *
from mpl_toolkits.mplot3d import axes3d

im1 = array(Image.open('001.jpg'))
im2 = array(Image.open('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_values='*')
P = [camera.Camera(loadtxt('2D/00'+str(i+1)+'.P')) for i in range(3)]


corr = corr[:,0] 
ndx3D = where(corr>=0)[0] 
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])) )
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.3 由基础矩阵计算照相机矩阵

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

在没有任何照相机内参数知识的情况下，照相机矩阵只能通过射影变换恢复出来，在无标定的情况下，第二个照相机矩阵可以使用一个3*3的射影变换得出。
$$\begin{gathered} P_2=[S_{e}F|\mathbf{e}] \\ {\text{其中,e 是左极点,满足 e}}^{T}F=0,S_e是如公式所示的反对称矩阵。 \end{gathered}$$

代码如下：

def skew(a):
    return array([[0, -a[2], a[1]], [a[2], 0, -a[0]], [-a[1], a[0], 0]])


def compute_P_from_fundamental(F):
    e = compute_epipole(F.T) 
    Te = skew(e)
    return vstack((dot(Te, F.T).T, e)).T

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

已经标定情况下，重建会保持欧式空间中的一些度量特性。
给定标定矩阵$K$,我们可以将它的逆$K^{-1}$作用于图像点${\mathbf{x}}_{\kappa }=K^{-1}\mathbf{x}$,因此，在新的图像
坐标系下，照相机方程变为：
$${\mathbf{x}}_K=K^{-1}K[R|\mathbf{t}]\mathbf{X}=[R|\mathbf{t}]\mathbf{X}$$
在新的图像坐标系下，点同一满足基础矩阵方程：${\mathbf{x}}_{K_2}^{T}F{\mathbf{x}}_{K_1}=0$
从本质矩阵恢复出的照相机矩阵中存在度量关系，但有四个可能解。
下面是计算这4个解的算法：

def compute_P_from_essential(E):

    U, S, V = svd(E)
    if det(dot(U, V)) < 0:
        V = -V
    E = dot(U, dot(diag([1, 1, 0]), V))

  
    Z = skew([0, 0, -1])
    W = array([[0, -1, 0], [1, 0, 0], [0, 0, 1]])

    P2 = [vstack((dot(U, dot(W, V)).T, U[:, 2])).T,
          vstack((dot(U, dot(W, V)).T, -U[:, 2])).T,
          vstack((dot(U, dot(W.T, V)).T, U[:, 2])).T,
          vstack((dot(U, dot(W.T, V)).T, -U[:, 2])).T]

    return P2

该函数确保本质矩阵秩为2。

3 多视图重建

假设照相机已经标定，步骤如下:
(1)检测特征点，然后在两幅图像间匹配$,
(2) 由匹配计算基础矩阵，
(3) 由基础矩阵计算照相机矩阵，
(4)三角剖分这些三维点。
当图像间的点对应包含不正确匹配时，需要一个稳健的方法来计算基础矩阵。

3.1 稳健估计基础矩阵

使用RANSAC方法，这次结合了八点算法（八点算法在平面场景中会失效）。
首先添加类到sfm.py文件中

class RansacModel(object):
   
    def __init__(self, debug=False):
        self.debug = debug

    def fit(self, data):

        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):
        """ Compute x^T F x for all correspondences,
            return error for each transformed point. """

        data = data.T
        x1 = data[:3]
        x2 = data[3:]

        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 = compute_fundamental(x1, x2)

    F = dot(T1.T, dot(F, T2))

    return F / F[2, 2]

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

将下列函数添加到sfm.py文件中：

def F_from_ransac(x1, x2, model, maxiter=5000, match_theshold=1e-6):
    import ransac

    data = vstack((x1, x2))

    # compute F and return with inlier index
    F, ransac_data = ransac.ransac(data.T, model, 8, maxiter, match_theshold, 20, return_all=True)
    return F, ransac_data['inliers']

3.2 三维重建

代码如下：

import sfm
import camera
import numpy as np
from PIL import Image
from pylab import *
import sift
import homography
from mpl_toolkits.mplot3d import axes3d

im1 = array(Image.open('001.jpg'))
im2 = array(Image.open('002.jpg'))

sift.process_image('001.jpg', 'im1.sift')
l0, d0 = sift.read_features_from_file('im1.sift')
sift.process_image('002.jpg', 'im2.sift')
l1, d1 = 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()
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()

3.3 多视图的扩展示例

1 多视图

当由同一场景的多个视图时，三维重建会变得更准确，包含更多细节信息。
对于视频序列，可以使用时序关系，在连续的帧对中匹配特征。相对方位需要从每个帧对增量地加入下一个帧对。
对于静止的图像，一个办法是找到一个中央参考视图，然后计算与之有关的所有其他照相机矩阵，另一个办法是计算一个图像对的照相机矩阵和三维重建，然后增量地加入新的图像和三维点。

2 光束法平差

多视图重建的最后一步，通常是通过优化三维点的位置和照相机参数来减少二次投影误差。

3 自标定

在未标定照相机的情形中，有时可以从图像特征来计算标定矩阵。该过程称为自标定。

4 立体图像

一个多视图成像的特殊例子是立体视觉(或者立体成像),即使用两台只有水平(向一侧)偏移的照相机观测同一场景。当照相机的位置如上设置，两幅图像具有相同的图像平面，图像的行是垂直对齐的，那么称图像对是经过矫正的。该设置在机器人学中很常见，常被称为立体平台。
通过将图像扭曲到公共的平面上，使外极线位于图像行上，任何立体照相机设置都能得到矫正(我们通常构建立体平台来产生经过矫正的图像对)。

立体重建 (有时称为致密深度重建)就是恢复深度图 (或者相反，视差图),图像中每个像素的深度 (或者视差)都需要计算出来。

• 计算视差图
  在该立体重建算法中，我们将对于每个像素尝试不同的偏移，并按照局部图像周围归一化的互相关值，选择具有最好分数的偏移，然后记录下该最佳偏移。因为每个偏移在某种程度上对应于一个平面，所有该过程又是称为扫平面法。下面是扫平面法代码，该函数返回每个像素的最佳视差。

def plane_sweep_ncc(im_l,im_r,start,steps,wid):
    m,n = im_l.shape
    mean_l = np.zeros((m,n))
    mean_r = np.zeros((m,n))
    s = np.zeros((m,n))
    s_l = np.zeros((m,n))
    s_r = np.zeros((m,n))
    dmaps = np.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(np.roll(norm_l,-displ-start)*norm_r,wid,s)
        filters.uniform_filter(np.roll(norm_l,-displ-start)*np.roll(norm_l,-displ-start),wid,s_l)
        filters.uniform_filter(norm_r*norm_r,wid,s_r)
        dmaps[:,:,displ] = s/np.sqrt(s_l*s_r)
    return np.argmax(dmaps,axis=2)