python计算机视觉第五章----多视图几何

多视图几何

多视图几何主要研究用几何的方法，通过若干个多个视角的二维图像，来恢复三维物体。其中，最重要的内容是双视图几何。

1、外极几何

假设有一场景点X，其经过单应性矩阵H变换后变为$HX$，则$HX$在照相机 $P{H}^{-1}$后得到图像点为$P{H}^{-1}HX=PX$,与图像点X在照相机P里得到的图像点相同。因此，在分析双视图几何关系时总是可以将照相机的相对位置关系用单应性矩阵加以简化。一个比较好的做法是将原点和坐标轴与第一个照相机对齐，$P_1=K_1[I|0],P_2=K_2[R|t]$,利用$K_1,K_2,R$等参数矩阵可以找到点X的投影点之间的关系，继而恢复照相机参数矩阵。具体方法如下：
已知，同一个图像点经过不同的投影矩阵产生的不同投影点必须满足:
$$x_2^{T}F{x}_1=0$$.
其中，矩阵 F 为基础矩阵。基础矩阵可以由两照相机的参数矩阵(相对旋转 R 和平移 t)表示$$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]$$
借助$F$我们就可以恢复出照相机参数。
利用公式$F=k_2^{-T}{S}_{1}R{K}_1^{-1}$可以找到对应第一条图像的一条直线，这条线称为对应于点$x_2$的外极线，点$x_2$在第一幅图像上的对应点一定在这条直线上。外极线都经过一点e，称为外极点。因为外几点在所有的外极线上，所以$F{e}_1=0$,因此可通过计算F的零向量来计算出外极点。

1.1 简单数据集

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

1.2 用matplotlib绘制三维数据

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

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from numpy import loadtxt
points3D = loadtxt('d:/imdata/3D/p3d').T
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')
fig = plt.figure()
ax = fig.gca(projection='3d') 
ax.plot(points3D[0], points3D[1], points3D[2], 'k.')
plt.show()

1.3 八点法

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

其中，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个方程即可求出。算法中只需要8个对应点来计算基础矩阵 F，所以该算法叫做八点法。具体实现代码如下：

 #计算基础矩阵F
 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)
    #受限F
    #最后一个奇异值置 0，使F秩为 2 
    U,S,V = linalg.svd(F)
    S[2] = 0
    F = dot(U,dot(diag(S),V))
    return F

1.4 外极点和外极线

因为外极点满足$F{e}_1=0$,故计算F的零空间即可得到外极点：

#求外极点
def compute_epipole(F):
    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 numpy import *
im1 = array(Image.open('d:/imdata/images/001.jpg'))
im2 = array(Image.open('d:/imdata/images/002.jpg'))
points2D = [loadtxt('d:/imdata/2D/00' + str(i + 1) + '.corners').T for i in range(3)]
corr = genfromtxt('d:/imdata/2D/nview-corners', dtype='int', missing_values='*')
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.computer_fundamental(x1,x2)
e=sfm.computer_epipole(F)
figure()
imshow(im1)
for i in range(5):
    sfm.plot_epipolar_line(im1,Fx2[:,i],e,False)
axis('off')
figure()
imshow(im2)
for i in range(5):
    plot(x2[0,i],x2[1,i],'0')
axis('off')
show()

实验结果如下：

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

2.1 三角剖分

给定照相机参数模型，图像点可以通过三角剖分来恢复出这些点的三维位置。对于两个照相机$P_1$和$P_2$的视图，三维实物X的投影点为$x_1$和$x_2$，照相机方程定义了下列关系:
$$\left[\begin{array}{ccc}{P}_1& -x_1& 0\\ P_2& 0& -x_2\end{array}\right]\left[\begin{array}{c}X\\ {\lambda }_1\\ {\lambda }_2\end{array}\right]=0$$
因为噪声等原因的干扰，上面的方程可能没有i精确解，我们通过SVD来得到三维点的最小二乘估值：

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

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]

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')
axis()
show()

控制台输出结果如下：

上图是我们估计出的三维图像点，下图是创建者提供的图像点，可以发现它们的元素几乎完全相同。

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

每个三维点 Xi(齐次坐标系下)按照 ${\lambda }_i{x}_i=P{X}_i$ 投影到图像点 $x_i=[x_i,y_i,1]$，相应的点满足下面的关系:

其中 p1、p2 和 p3 是矩阵 P 的三行。上面的式子可以写得更紧凑，如下所示:
$$Mv=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))

下面来做一个简单的实验：

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

可以看出我们估计出的照相机矩阵与提供的照相机矩阵几乎完全相同。

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

照相机矩阵可以由基础矩阵恢复出来。

2.3.1 未标定的情况

在没有任何照相机内参数知识的情况下，照相机矩阵只能通过射影变换恢复出来。也就是说，如果利用照相机的信息来重建三维点，那么该重建只能由射影变换计算出来，因此，在无标定的情况下，第二个照相机矩阵可以使用一个（3×3）的射影变换得出。

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
def skew(a):
    """反对称矩阵A，使得对于每个v有a×v=Av""" 
    return array([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])

2.3.1 已标定的情况

在已经标定的情况下，重建会保持欧式空间中的一些度量特性。给定标定矩阵K，我们可以将它的逆$K_{-1}$作用于图像点$x_k=K^{-1}x$，因此，在新的图像坐标系下，照相机方程变为：
$$x_K=K^{-1}K[R\mid t]X=[R\mid t]X$$
在新的图像坐标系下，点同样满足之前的基础矩阵方程。从本质矩阵中恢复出的照相机矩阵中存在度量关系，但有四个可能解。因为只有一个解产生位于两个照相机前的场景，所以可以从中选出来。
具体代码如下：

def compute_P_from_essential(E):
    """从本质矩阵中计算第二个照相机矩阵（假设 P1 = [I 0])
    输出为4个可能的照相机矩阵列表"""
    # 保证E的秩为2
    U, S, V = svd(E)
    if det(dot(U, V)) < 0:
        V = -V
    E = dot(U, dot(diag([1, 1, 0]), V))
    # 创建矩阵（Hartley）
    Z = skew([0, 0, -1])
    W = array([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
    # 返回所有（4个）解
    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

3、多视图重建

由于照相机的运动给我们提供了三维结构，所以这样计算三维重建的方法通常称为 SfM (Structure from Motion，从运动中恢复结构)。假设照相机已经标定，计算重建可以分为下面 4 个步骤:
(1) 检测特征点，然后在两幅图像间匹配;
(2) 由匹配计算基础矩阵;
(3) 由基础矩阵计算照相机矩阵;
(4) 三角剖分这些三维点。
我们已经具备了完成上面 4 个步骤所需的所有工具。但是当图像间的点对应包含不正确的匹配时，我们需要一个稳健的方法来计算基础矩阵。

3.1 稳健估计基础矩阵

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

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.2 多视图的扩展示例

1. 多视图
   当我们有同一场景的多个视图时，三维重建会变得更准确，包括更多的细节信息。 因为基础矩阵只和一对视图相关，所以该过程带有很多图像，和之前的处理有些 不同。
   对于视频序列，我们可以使用时序关系，在连续的帧对中匹配特征。相对方位需要 从每个帧对增量地加入下一个帧对，同时可以使用跟踪有效地找到对应点。
   对于静止的图像，一个办法是找到一个中央参考视图，然后计算与之有关的所有其他照相机矩阵。另一个办法是计算一个图像对的照相机矩阵和三维重建，然后增量地加入新的图像和三维点，参见参考文献 [34]。另外，还有一些方法可以同时由三个视图计算三维重建和照相机位置

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

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