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

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


# 载入一些图像
im1 = array(Image.open('ch05\\5.1\\images\\001.jpg'))
im2 = array(Image.open('ch05\\5.1\\images\\002.jpg'))

# 载入每个视图的二维点到列表中
points2D = [loadtxt('ch05/5.1/2D/00'+str(i+1)+'.corners').T for i in range(3)]
# 载入三维点
points3D = loadtxt('ch05/5.1/3D/p3d').T
# 载入对应
corr = genfromtxt('ch05/5.1/2D/nview-corners',dtype='int',missing_values='*')
# 载入照相机矩阵到Camera 对象列表中
P = [camera.Camera(loadtxt('ch05/5.1/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()

实验分析：
仔细观察这些点的分布情况可知，图像点和三维点有些点是重合的，但是并不是完全重合。蓝色这些点代表的是Merton1数据集中的点，红色的点代表的是特征映射过来的坐标。如下图所示。

5.1.2 用Matplotlib绘制三维数据

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

根据上面Merton1下载的数据信息，绘制出该图像的三维点图，编写代码如下：

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

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

# -*- coding: utf-8 -*-
from pylab import *
from matplotlib.pyplot import figure, show
import camera
import camera
import numpy as np
from PIL import Image
from pylab import *

# 载入一些图像
im1 = array(Image.open('ch05\\5.1\\images\\001.jpg'))
im2 = array(Image.open('ch05\\5.1\\images\\002.jpg'))

# 载入每个视图的二维点到列表中
points2D = [loadtxt('ch05/5.1/2D/00'+str(i+1)+'.corners').T for i in range(3)]
# 载入三维点
points3D = loadtxt('ch05/5.1/3D/p3d').T
# 载入对应
corr = genfromtxt('ch05/5.1/2D/nview-corners',dtype='int',missing_values='*')
# 载入照相机矩阵到Camera 对象列表中
P = [camera.Camera(loadtxt('ch05/5.1/2D/00'+str(i+1)+'.P')) for i in range(3)]


plt.rcParams['font.sans-serif']=['SimHei']  # 用来正常显示中文标签  
plt.rcParams['axes.unicode_minus']=False  # 用来正常显示负号

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

俯视的视图，展示了建筑墙体和屋顶上的点：

5.1.3 计算F ：八点法

八点法是通过对应点来计算基础矩阵的算法。

八点算法的原理：
由于基础矩阵F定义为：

\(x^{T}F{x}'=0\)

任给两幅图像中的匹配点x与x'，令\(x=(u,v,1)^T,x'=(u',v',1)^T\)，基础矩阵F是一个3×3的秩为2的矩阵，一般记基础矩阵F为：

\(\text{F}=\begin{bmatrix}\text{f}_{11}&\text{f}_{12}&\text{f}_{13}\\\text{f}_{21}&\text{f}_{22}&\text{f}_{23}\\\text{f}_{31}&\text{f}_{32}&\text{f}_{33}\end{bmatrix}\)

有相应方程：

\(\mathrm u\mathrm u'\mathrm f_{11}+\mathrm u\mathrm v'\mathrm f^{12}+\mathrm u\mathrm f^{13}+\mathrm v\mathrm u'\mathrm f^{21}+\mathrm v\mathrm v'\mathrm f^{22}+\mathrm v\mathrm f^{23}+\mathrm u'\mathrm f_{31}+\mathrm v'\mathrm f_{32}+\mathrm f_{33}=0\)

由矩阵乘法可知有：

\(\begin{bmatrix}x_2^1x_1^1&x_2^1y_1^1&x_2^1w_1^1&\cdots&w_2^1w_1^1\\x_2^2x_1^2&x_2^2y_1^2&x_2^2w_1^2&\cdots&w_2^2w_1^2\\\vdots&\vdots&\vdots&\ddots&\vdots\\x_2^nx_1^n&x_2^ny_1^n&x_2^nw_1^n&\cdots&w_2^nw_1^n\end{bmatrix}\begin{bmatrix}F_{11}\\F_{12}\\F_{13}\\\vdots\\F_{33}\end{bmatrix}=Af=0\)

在实际计算中，可以直接用\(A^{T}A\)的分解来求解参数。也可以用非线性优化，通过搜索f使得\(||Af||\)最小化, 同时满足\(||f||=1\)的约束。上述求解后的F不一定能满足秩为2的约束，因此还要在F基础上加以约束。

因为算法中需要8 个对应点来计算基础矩阵F，所以该算法叫做八点法。

新建一个文件sfm.py，写入下面八点法中最小化\(||Af||\)的函数：

from matplotlib.pyplot import plot
from numpy import *
from scipy import linalg


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):
    """在图像中，绘制外极点和外极线 Fx=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*')

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

5.1.4 外极点和外极线

外极点满足\(F_(e_1))\)，因此可以通过计算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 转置后输入上述函数即可。

在之前样本数据集的前两个视图上运行这两个函数：

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

# 载入一些图像
im1 = array(Image.open('ch05\\5.1\\images\\001.jpg'))
im2 = array(Image.open('ch05\\5.1\\images\\002.jpg'))

# 载入每个视图的二维点到列表中
points2D = [loadtxt('ch05/5.1/2D/00'+str(i+1)+'.corners').T for i in range(3)]
# 载入三维点
points3D = loadtxt('ch05/5.1/3D/p3d').T
# 载入对应
corr = genfromtxt('ch05/5.1/2D/nview-corners',dtype='int',missing_values='*')
# 载入照相机矩阵到Camera 对象列表中
P = [camera.Camera(loadtxt('ch05/5.1/2D/00'+str(i+1)+'.P')) for i in range(3)]

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

结果：
选择两幅图像的对应点，然后将它们转换为齐次坐标，这里的对应点是从一个文本文件中读取得到的；而实际上，可以按照sift提取图像特征的方式，然后通过匹配来找到它们。由于缺失的数据在对应列表corr中为-1，所以程序中有可能选取这些点。因此，上面的程序通过数组操作符&只选取了索引大于等于0的点。最后，在第一个视图中画出了前5条外极线，在第二个视图中画出了对应5个匹配点。还借助了plot_epipolar_line()函数，这个函数将x轴的范围作为直线的参数，所以直线超出了图像边界的部分会被截断。如果show_epipole为真，外极点会被画出来（如果输入参数中没有外极点，外极点会在程序中计算获得）。用不同的颜色将点和对应的外极线对应起来。

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

5.2.1 三角剖分

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

数学定义： 假设V是二维实数域上的有限点集，边e是由点集中的点作为端点构成的封闭线段, E为e的集合。那么该点集V的一个三角剖分T=(V,E)是一个平面图G，该平面图满足条件：

　　　　a.除了端点，平面图中的边不包含点集中的任何点。

　　　　b.没有相交边。

　　　　c.平面图中所有的面都是三角面，且所有三角面的合集是散点集V的凸包。

对于两个照相机P1和P2的视图，三维实物点X的投影点为x1和x2（用齐次坐标表示），照相机方程定义了下面的关系：

\(\begin{bmatrix}P_1&-\mathbf{x}_1&0\\P_2&0&-\mathbf{x}_2\end{bmatrix}\begin{bmatrix}\mathbf{X}\\\lambda_1\\\lambda_2\end{bmatrix}=0\)

由于图像噪声、照相机参数误差和其他系统误差，上面的方程可能没有精确解。可以通过SVD算法来得到三维点的最小二乘估值。

下面的函数用于计算一个点对的最小二乘三角剖分，把triangulate_point()和triangulate()函数添加到sfm.py中：

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

第一个函数用于计算一个点对的最小二乘三角剖分，它的最后一个特征向量的前4个值就是齐次坐标系下的对应三维坐标。第二个函数用于实现多个点的三角剖分，这个函数的输入是两个图像点数组，输出为一个三维坐标数组。

利用下面的代码来实现Merton1 数据集上的三角剖分：

from pylab import *
from mpl_toolkits.mplot3d import axes3d
import sfm

import camera
from PIL import Image
from pylab import *


# 载入一些图像
im1 = array(Image.open('ch05\\5.1\\images\\001.jpg'))
im2 = array(Image.open('ch05\\5.1\\images\\002.jpg'))

# 载入每个视图的二维点到列表中
points2D = [loadtxt('ch05/5.1/2D/00'+str(i+1)+'.corners').T for i in range(3)]
# 载入三维点
points3D = loadtxt('ch05/5.1/3D/p3d').T
# 载入对应
corr = genfromtxt('ch05/5.1/2D/nview-corners',dtype='int',missing_values='*')
# 载入照相机矩阵到Camera 对象列表中
P = [camera.Camera(loadtxt('ch05/5.1/2D/00'+str(i+1)+'.P')) for i in range(3)]



# 在前两个视图中点的索引
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('auto')

show()

这个算法估计出的三维图像点和实际图像点很接近。估计点和实际点可以很好地进行匹配。

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

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

从照相机方程可以得出，每个三维点\(X_{i}\)（齐次坐标系下）按照\( \lambda _{i}x_{i}=PX_{i}\)投影到图像点\(x_{i}=[x_{i},y_{i},1]\)，相应的点满足下面的关系：

\(\begin{bmatrix}\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\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\end{bmatrix}\begin{bmatrix}\mathbf{p}_1^T\\\mathbf{p}_2^T\\\mathbf{p}_3^T\\\lambda_1\\\lambda_2\\\vdots\end{bmatrix}=0\)

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

\(\text{Mv=0}\)

我们可以使用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.")
    # 创建用于计算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))

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

import camera
from PIL import Image
from pylab import *
import sfm

# 载入一些图像
im1 = array(Image.open('ch05\\5.1\\images\\001.jpg'))
im2 = array(Image.open('ch05\\5.1\\images\\002.jpg'))

# 载入每个视图的二维点到列表中
points2D = [loadtxt('ch05/5.1/2D/00'+str(i+1)+'.corners').T for i in range(3)]
# 载入三维点
points3D = loadtxt('ch05/5.1/3D/p3d').T
# 载入对应
corr = genfromtxt('ch05/5.1/2D/nview-corners',dtype='int',missing_values='*')
# 载入照相机矩阵到Camera 对象列表中
P = [camera.Camera(loadtxt('ch05/5.1/2D/00'+str(i+1)+'.P')) for i in range(3)]

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

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

真实点用圆圈表示，估计出的照相机投影点用点表示。

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

在两个视图场景中，照相机矩阵可以由基础矩阵恢复出来。假设第一个照相机矩阵归一化为P1=[I/0]，然后需要计算出第二个照相机矩阵P2。分成两类情况，一类是未标定，另一类是已标定。

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

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

\(P_{2}=[S_{e}F|e]\)

其中，e是左极点，满足\(e^{T}F=0\)，\(S_{e}\)是公式所示的反对称矩阵。需要注意的是，使用该矩阵做出的三角形剖分很有可能会发生畸变，如倾斜的重建。

下面是具体的代码，把它添加到sfm.py中：

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

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

在已经标定的情况下，重建会保持欧式空间中的一些度量特性（除了全局的尺度参数）。
定义：给定标定矩阵K，可以将它的逆\(K^{-1}\)作用于图像点\(x_{K}=K^{-1}x\)，因此，在新的图像坐标系下，照相机方程变为：

\(\mathrm x_\mathrm K=\mathrm K^{-1}\mathrm K[\mathrm R|\mathrm t]\mathrm X=[\mathrm R|\mathrm t]\mathrm X\)

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

\(\mathbf{x}_{K_{2}}^{T}F\mathbf{x}_{K_{1}}=0\)

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

从本质矩阵中恢复出的照相机矩阵中存在度量关系，但有四个可能解。因为只有一个解产生位于两个照相机前的场景，所以可以从中选出来。

下面是具体计算4个解的代码，把它添加到sfm.py中：

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

5.3 多视图重建

下面介绍如何使用从多幅图像中计算出真实的三维重建。由于照相机的运动给我们提供了三维结构，所以这样计算三维重建的方法通常称为SFM（从运动中恢复结构）。
假设照相机已经标定，计算重建可以分为下面的4个步骤：

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

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

5.3.1 稳健估计基础矩阵

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

这个需要将下面的代码添加到sfm.py文档中即可：

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

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]

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

其中，compute_fundamental_normalized()函数将图像点归一化为零均值固定方差。

5.3.2 三维重建示例

import homography
from matplotlib.pyplot import *
from numpy import *
import sfm
import sift
import camera
from PIL import Image
from pylab import *

from PCV.geometry import homography, camera, sfm

# 标定矩阵
K = array([[2394,0,932],[0,2398,628],[0,0,1]])
# 载入图像，并计算特征
im1 = array(Image.open('pic\\alcatraz1.jpg'))
sift.process_image('pic\\alcatraz1.jpg','ch05\\5.2\im1.sift')
l1,d1 = sift.read_features_from_file('ch05\\5.2\im1.sift')

im2 = array(Image.open('pic\\alcatraz2.jpg'))
sift.process_image('pic\\alcatraz2.jpg','ch05\\5.2\im2.sift')
l2,d2 = sift.read_features_from_file('ch05\\5.2\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)
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()

5.3.3 多视图的扩展示例

多视图重建有一些步骤和进一步的扩展，下面提供关于一些有关内容的介绍。

多视图

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

        对于视频序列，可是使用时序关系，在连续的帧对中匹配特征。相对方位需要从每个帧对增量地加入下一个帧对。同时可以使用跟踪有效地找到对应点。

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

光束法平差

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

自标定

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

5.4 立体图像

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

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

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

\(Z=\frac{fb}{x_l-x_r}\)

其中，f是经过矫正图像的焦距，b是两个照相机中心之间的距离，\(x_{l}\)和\(x_{r}\)左右两幅图像中对应点的x坐标。分开照相机中心的距离称为基线。

立体重建就是恢复深度图（或者相反，视差图），图像中每个像素的深度都需要计算出来。

腾讯云开发者社区

腾讯云面向开发者汇聚海量精品云计算使用和开发经验，营造开放的云计算技术生态圈。

更多推荐

Elasticsearch复杂数据类型终极指南：从入门到精通

Elasticsearch作为功能强大的搜索引擎，支持多种复杂数据类型，让开发者能够灵活处理各种结构化和非结构化数据。本文将带你全面了解Elasticsearch中的复杂数据类型，从基础概念到实际应用，助你轻松掌握数据建模的核心技巧。## 内部对象：构建层级化数据结构在Elasticsearch中，对象类型（Object）是最基础的复杂数据类型之一，用于表示具有嵌套关系的数据。例如，我们可

腾讯云开发者社区

终极指南：Flink SQL连接器版本管理从混乱到有序的升级之路

Apache Flink作为流处理领域的佼佼者，其SQL连接器的版本管理一直是开发者面临的核心挑战。本文将系统讲解Flink SQL连接器版本管理的最佳实践，帮助你轻松应对版本兼容性问题，实现从混乱到有序的升级之旅。## 连接器版本管理的常见痛点 😫在Flink应用开发中，连接器版本管理常常让开发者头疼不已。不同版本的连接器可能导致各种兼容性问题，例如API变更、功能差异甚至运行时错误。

腾讯云开发者社区

如何快速搭建Neon无服务器PostgreSQL：面向初学者的完整指南

Neon是一款革命性的无服务器PostgreSQL解决方案，它通过分离存储和计算层，实现了自动扩缩容、类代码式数据库分支以及零级扩展能力。本指南将帮助你从零开始搭建Neon开发环境，体验这款创新数据库的强大功能。## 准备工作：环境要求与依赖项在开始搭建Neon环境前，请确保你的系统满足以下要求：- Linux操作系统（推荐Ubuntu 20.04+或Debian 11+）- Git