python导入几何模型

多视图几何

• （一）外极几何

• 1.1 基础矩阵
• 1.2 一个简单的数据集
• 1.3 用 Matplotlib 绘制三维数据
• 1.4 外极点和外极线

• （二）照相机和三维结构的计算

• 2.1 三角剖分
• 2.2 由三维点计算照相机矩阵
• 2.3 由基础矩阵计算照相机矩阵

• （三）三维重建

• 3.1 稳健估计基础矩阵
• 3.2 三维重建

• （四）问题求助

（一）外极几何

如果有一个场景的两个视图以及视图中的对应图像点，那么根据照相机间的空间相对位置关系、照相机的性质以及三维场景点的位置，可以得到对这些图像点的一些几何关系约束。 外极几何是研究两幅图像之间存在的几何。它和场景结构无关，只依赖于摄像机的内外参数。研究这种几何可以用在图像匹配、三维重建方面。

介绍：

1.1 基础矩阵

1.基础矩阵概述：

1）两幅图像之间的约束关系使用代数的方式表示出来即为基本矩阵。

2）基础矩阵F满足：

3）基础矩阵可以用于简化匹配和去除错配特征

2.基础矩阵性质:

• 转置对称性：如果$F$是一对影像$(P,P′)$的基础矩阵（$x′Fx=0$），反过来 $(P$ 的基础矩阵是$FT$。证明很简单，直接对$x′Fx=0$两侧分别转置，得到 $xTFTx′=0$。
• 对极线：对于左影像上任意一点 $x$，其在右影像上的对极线为 $l$
• 对极点：任何对极线都会经过核点，所以有对于左影像上任意一点$x$ ，$e′Tl′=e′T(Fx)=0$，于是有$e′TF=0$。同理有$Fe=0$。
• F具有7自由度：一个33x33 的单应矩阵，具有88个自由度，而 FF 还满足 $detF=0$，所以$F$
• F是相关的：$F$ 将左影像上的一点 $x$ 投影到右影像上一条对极线$l′$，投影本质上是将$x$ 与左极点的连线$l$投影到右影像上的对极线 $l$，所以右影像上的一条对极线 $l$对应的是左影像上的一条对极线$l$，这种点到线的投影不可逆。

2.八点算法估算基础矩阵F

(1) **原理：**任给两幅图像中的匹配点$x$ 与$x’$ ，令 $x=(u,v,1)^T$ ,

基础矩阵F是一个3×3的秩为2的矩阵,一般记基础矩阵F为:

有相应方程：

由矩阵乘法可知有:

　

(2）优点: 线性求解，容易实现，运行速度快 。

缺点：对噪声敏感。

归一算法步骤： 给定n≥8组对应点{$Xi​<−>Xi′​$}，确定基本矩阵$F$使得$x$。

• 归一化：根据$Xi​=TXi​,Xi=T′X′i$,变换图像坐标。其中$T$和$T′$是有平移和缩放组成的归一变化。
• 求解对应匹配的基本矩阵$F′$。
• 求线性解:用由对应点集{Xi<−>x′iXi<−>xi′}确定的系数矩阵AA的最小奇异值的奇异矢量确定$F$
• 奇异性约束:用SVD对$F$进行分解，令其最小奇异值为0,得到$F′$，使得$detF′=0$。
• 解除归一化： 令F=T’TF’T F=T′TF′T。矩阵FF就是数据xi<−>x′ixi<−>xi′对应的基本矩阵。

示例：

# -*- coding: utf-8 -*-
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
from PCV.geometry import camera
from PCV.geometry import homography
from PCV.geometry import sfm
from PCV.localdescriptors import sift
# -*- coding: utf-8 -*-

载入图像，并计算特征， 匹配特征：

im1 = array(Image.open('E:\\test\\2.jpg'))
sift.process_image('E:\\test\\2.jpg', 'E:\\test\\im1.sift')
l1, d1 = sift.read_features_from_file('E:\\test\\im1.sift')

im2 = array(Image.open('E:\\test\\3.jpg'))
sift.process_image('E:\\test\\3.jpg', 'E:\\test\\im2.sift')
l2, d2 = sift.read_features_from_file('E:\\test\\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 = x1.copy()
x2n = x2.copy()
print(len(ndx))

figure(figsize=(16,16))
sift.plot_matches(im1, im2, l1, l2, matches, True)
show()

# Don't use K1, and K2

#def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
    """ Robust estimation of a fundamental matrix F from point
    correspondences using RANSAC (ransac.py from
    http://www.scipy.org/Cookbook/RANSAC).

    input: x1, x2 (3*n arrays) points in hom. coordinates. """

    from PCV.tools import ransac
    data = np.vstack((x1, x2))
    d = 20 # 20 is the original
    # compute F and return with inlier index
    F, ransac_data = ransac.ransac(data.T, model,
                                   8, maxiter, match_threshold, d, return_all=True)
    return F, ransac_data['inliers']

使用 RANSAC 方法估计 E：

model = sfm.RansacModel()
F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-1)
print("F:", F)

P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_fundamental(F)

print("P2:", P2)

# triangulate inliers and remove points not in front of both cameras
X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)

# plot the projection of X
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2)
x1p = cam1.project(X)
x2p = cam2.project(X)

figure(figsize=(16, 16))
imj = sift.appendimages(im1, im2)
imj = vstack((imj, imj))

imshow(imj)

cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1p[0])):
    if (0<= x1p[0][i]<cols1) and (0<= x2p[0][i]<cols1) and (0<=x1p[1][i]<rows1) and (0<=x2p[1][i]<rows1):
        plot([x1p[0][i], x2p[0][i]+cols1],[x1p[1][i], x2p[1][i]],'c')
axis('off')
show()

d1p = d1n[inliers]
d2p = d2n[inliers]

# Read features
im3 = array(Image.open('picture/home/3.jpg'))
sift.process_image('picture/home/3.jpg', 'im3.sift')
l3, d3 = sift.read_features_from_file('im3.sift')

matches13 = sift.match_twosided(d1p, d3)

ndx_13 = matches13.nonzero()[0]
x1_13 = homography.make_homog(x1p[:, ndx_13])
ndx2_13 = [int(matches13[i]) for i in ndx_13]
x3_13 = homography.make_homog(l3[ndx2_13, :2].T)

figure(figsize=(16, 16))
imj = sift.appendimages(im1, im3)
imj = vstack((imj, imj))

imshow(imj)

cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1_13[0])):
    if (0<= x1_13[0][i]<cols1) and (0<= x3_13[0][i]<cols1) and (0<=x1_13[1][i]<rows1) and (0<=x3_13[1][i]<rows1):
        plot([x1_13[0][i], x3_13[0][i]+cols1],[x1_13[1][i], x3_13[1][i]],'c')
axis('off')
show()

P3 = sfm.compute_P(x3_13, X[:, ndx_13])

效果：

sift特征匹配：特征匹配点

基础矩阵：

1.2 一个简单的数据集

想要从书上给的网址 http://www.robots.ox.ac.uk/~vgg/data/data-mview.html 中下载牛津多视图数据集下载Merton2数据的压缩文件。下面的脚本可以加载Merton2的数据： 示例：

# -*- coding: utf-8 -*-

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

im1 = array(Image.open('picture/images (1)/001.jpg'))
im2 = array(Image.open('picture/images (1)/002.jpg'))

points2D = [loadtxt('picture/2D (1)/00'+str(i+1)+'.corners').T for i in range(3)]
points3D = loadtxt('picture/3D (1)/p3d').T
corr = genfromtxt('picture/2D (1)/nview-corners',dtype='int',missing_values='*')
P = [camera.Camera(loadtxt('picture/2D (1)/00'+str(i+1)+'.P')) for i in range(3)]

# make 3D points homogeneous and project
X = vstack((points3D,ones(points3D.shape[1])))
x = P[0].project(X)

# plotting the points in view 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()

处理效果：

分析：

仔细观察这些点的分布情况可知，图像点和三维点有些点是重合的，但是并不是完全重合。蓝色这些点代表的是Merton2数据集中的点，红色的点代表的是特征映射过来的坐标。图像特征点 1、对应不同视图图像点重建后的三维点以及照相机参数矩阵（使用上一章的 Camera 类）。这里使用 loadtxt() 函数读取文本文件到 NumPy 数组中。因为并不是所有的点都可见，或都能够成功匹配到所有的视图，所以对应数据里包含了缺失的数据。加载对应数据时需要考虑这一点。genfromtxt() 函数通过将缺失的数值（在文件中用 * 表示）填充为 -1 来解决这个问题。

1.3 用 Matplotlib 绘制三维数据

为了可视化三维重建结果，我们需要绘制出三维图像。Matplotlib 中的 mplot3d 工具 包可以方便地绘制出三维点、线、等轮廓线、表面以及其他基本图形组件，还可以 通过图像窗口控件实现三维旋转和缩放。 get_test_data() 函数在 x, y 空间按照设定的空间间隔参数来产生均匀的采样点。

# -*- coding: utf-8 -*-
from mpl_toolkits.mplot3d import axes3d

fig = figure()
ax = fig.gca(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()

效果：

初始显示：

俯视的视图：

分析： 通过旋转上面的坐标系，在不同三维角度显示这些三维数据点的分布。

1.4 外极点和外极线

极点满足$Fe1​=0$，因此可以通过计算FF的零空间来得到，将FF转置后输入到函数def compute_epipole(F) 中获得另一幅图像的外极点（对应左零空间的外极点）。

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

# -*- coding: utf-8 -*-
from pylab import *
from mpl_toolkits.mplot3d import axes3d
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
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()

分析：

首先，选择两幅图像的对应点，然后将它们转换为齐次坐标。 由于缺失的数据在对应列表 corr 中为 -1，所以程序中有可能选取这些点。因此，上面的程序通过数组操作符 & 只选取了索引大于等于 0 的点。最后，在第一个视图中画出了前 5 个外极线，在第二个视图中画出了对应匹配点。主要借助 plot() 函数，将 x 轴的范围作为直线的参数，因此直线超出图像边界的部分会被截断。 如果 show_epipole 为真，外极点也会被画出来（如果输入参数中没有外极点，外极点会在程序中计算获得）。在两幅图中，用不同的颜色将点和相应的外极线对应起来。

（二）照相机和三维结构的计算

2.1 三角剖分

数学定义： 假设V是二维实数域上的有限点集，边e是由点集中的点作为端点构成的封闭线段, E为e的集合。那么该点集V的一个三角剖分T=(V,E)是一个平面图G，该平面图满足条件： 　　　　a.除了端点，平面图中的边不包含点集中的任何点。 　　　　b.没有相交边。 　　　　c.平面图中所有的面都是三角面，且所有三角面的合集是散点集V的凸包。

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

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

将triangulate_point()和triangulate()函数添加到sfm.py中：

# -*- coding: utf-8 -*-
execfile('aaaaa.py')
# -*- coding: utf-8 -*-
from pylab import *
from mpl_toolkits.mplot3d import axes3d
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.gca(projection='3d')
ax.plot(Xest[0],Xest[1],Xest[2],'ko')
ax.plot(Xtrue[0],Xtrue[1],Xtrue[2],'r.')
axis('equal')
show()

计算一个点对的最小二乘三角剖分:

def triangulate_point(x1, x2, P1, P2):
    """ Point pair triangulation from
        least squares solution. """

    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]

最后一个特征向量的前 4 个值就是齐次坐标系下的对应三维坐标。可以增加下面的函数来实现多个点的三角剖分：

def triangulate(x1, x2, P1, P2):
    """    Two-view triangulation of points in
        x1,x2 (3*n homog. coordinates). """

    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

效果：

分析：

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

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

从照相机方程可以得出，每个三维点$Xi​$（齐次坐标系下）按照$λi​xi​=PXi​$ 投影到图像点$xi​=[xi​,yi​,1]$，相应的点满足下面的关系：

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

$Mv=0$

SVD分解估计出照相机矩阵。利用上面讲述的矩阵操作,选出第一个视图中的可见点（使用对应列表中缺失的数值），将它们转换成齐次坐标表示，然后估计照相机矩阵：

# -*- coding: utf-8 -*-

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

im1 = array(Image.open('picture/images (1)/001.jpg'))
im2 = array(Image.open('picture/images (1)/002.jpg'))

points2D = [loadtxt('picture/2D (1)/00'+str(i+1)+'.corners').T for i in range(3)]
points3D = loadtxt('picture/3D (1)/p3d').T
corr = genfromtxt('picture/2D (1)/nview-corners',dtype='int',missing_values='*')
P = [camera.Camera(loadtxt('picture/2D (1)/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()

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

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

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

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

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

在没有任何照相机内参数知识的情况下，照相机矩阵只能通过射影变换恢复出来。也就是说，如果利用照相机的信息来重建三维点，那么该重建只能由射影变换计算出来（可以得到整个投影场景中无畸变的重建点）。在这里，不考虑角度和距离。 定义： 在无标定的情况下，第二个照相机矩阵可以使用一个（3×3）的射影变换得到。 一个简单的方法就是：$P2​=[Se​F∣e]$ 其中，$e$是左极点，满足$eTF=0$，$Se$是公式所示的反对称矩阵。需要注意的是，使用该矩阵做出的三角形剖分很有可能会发生畸变，如倾斜的重建。

示例：

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$ 作用于图像点$xK​=K−1x$，因此，在新的图像坐标系下，照相机方程变为：$xK​=K−1K[R∣t]X=[R∣t]X$ 在新的图像坐标系下，点同样满足之前的基础矩阵方程：$xK2T​FxK1​=0$ 在标定归一化的坐标系下，基础矩阵称为本质矩阵。为了区别为标定后的情况，以及归一化了的图像坐标，通常将其记为$E$，而非$F$。 从本质矩阵中恢复出的照相机矩阵中存在度量关系，但有四个可能解。因为只有一个解产生位于两个照相机前的场景，所以可以从中选出来。

计算4个解:

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

实验步骤：

1. 用SIFT算法实现两幅图像的特征点检测，找到对应的匹配点并绘制出来；
2. 使用RANSAC方法估计最佳基础矩阵F，以及正确点的索引，求出不同图像对的基础矩阵；
3. 由基础矩阵计算照相机矩阵；
4. 从照相机矩阵的列表中，对正确点的三维点进行三角剖分，挑选出经过三角剖分后，在两个照相机前均含有最多场景点。

# -*- coding: utf-8 -*-

# -*- coding: utf-8 -*-
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
from PCV.geometry import camera
import homography
from PCV.geometry import sfm
from PCV.localdescriptors import sift
from imp import reload
camera = reload(camera)
homography = reload(homography)
sfm = reload(sfm)
sift = reload(sift)
# Read features
im1 = array(Image.open(' /1.jpg'))
im2 = array(Image.open('  /2.jpg'))
sift.process_image(' /1.jpg','im1.sift')
l0, d0 = sift.read_features_from_file('im1.sift')
sift.process_image(' /2.jpg','im2.sift')
l1, d1 = sift.read_features_from_file('im2.sift')

matches = sift.match_twosided(d1, d2)

ndx = matches.nonzero()[0]
x1 = homography.make_homog(l1[ndx, :2].T)
ndx2 = [int(matches[i]) for i in ndx]
x2 = homography.make_homog(l2[ndx2, :2].T)

d1n = d1[ndx]
d2n = d2[ndx2]
x1n = x1.copy()
x2n = x2.copy()


figure(figsize=(16, 16))
sift.plot_matches(im1, im2, l1, l2, matches, True)
show()

# def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
    """ Robust estimation of a fundamental matrix F from point
    correspondences using RANSAC (ransac.py from
    http://www.scipy.org/Cookbook/RANSAC).
    input: x1, x2 (3*n arrays) points in hom. coordinates. """

    from PCV.tools import ransac
    data = np.vstack((x1, x2))
    d = 10  # 20 is the original
    # compute F and return with inlier index
    F, ransac_data = ransac.ransac(data.T, model,
                                   8, maxiter, match_threshold, d, return_all=True)
    return F, ransac_data['inliers']

# find F through RANSAC
model = sfm.RansacModel()
F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-3)
print (F)

P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_fundamental(F)

print (P2)
print (F)
# triangulate inliers and remove points not in front of both cameras
X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)

# plot the projection of X
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2)
x1p = cam1.project(X)
x2p = cam2.project(X)
figure(figsize=(16, 16))
imj = sift.appendimages(im1, im2)
imj = vstack((imj, imj))

imshow(imj)

cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1p[0])):
    if (0 <= x1p[0][i] < cols1) and (0 <= x2p[0][i] < cols1) and (0 <= x1p[1][i] < rows1) and (0 <= x2p[1][i] < rows1):
        plot([x1p[0][i], x2p[0][i] + cols1], [x1p[1][i], x2p[1][i]], 'c')
axis('off')
show()

d1p = d1n[inliers]
d2p = d2n[inliers]

# Read features
im3 = array(Image.open(' /3.jpg'))
sift.process_image('  /3.jpg','im3.sift')
l3, d3 = sift.read_features_from_file('im3.sift')

matches13 = sift.match_twosided(d1p, d3)

ndx_13 = matches13.nonzero()[0]
x1_13 = homography.make_homog(x1p[:, ndx_13])
ndx2_13 = [int(matches13[i]) for i in ndx_13]
x3_13 = homography.make_homog(l3[ndx2_13, :2].T)

figure(figsize=(16, 16))
imj = sift.appendimages(im1, im3)
imj = vstack((imj, imj))

imshow(imj)

cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1_13[0])):
    if (0 <= x1_13[0][i] < cols1) and (0 <= x3_13[0][i] < cols1) and (0 <= x1_13[1][i] < rows1) and (
            0 <= x3_13[1][i] < rows1):
        plot([x1_13[0][i], x3_13[0][i] + cols1], [x1_13[1][i], x3_13[1][i]], 'c')
axis('off')
show()

P3 = sfm.compute_P(x3_13, X[:, ndx_13])

print (P3)
print (P1)
print (P2)
print (P3)

sift特征匹配结果：

处理过后：

照相机矩阵：

室内：

sift特征匹配：

特征匹配：

照相机矩阵：

分析：

由上面室外的结果看出，第一个照相机位置和第二个照相机位置匹配的结果相比之下，效果比较好；而第一个照相机位置和第三个照相机为止的匹配效果并没那么好，匹配的点比较少。从室内的结果看出，照相机的位置和图像呈现的匹配点有关系，拍摄的角度和清晰度会对图像匹配产生影响。

（三）三维重建

由于照相机的运动给我们提供了三维结构，从多幅图像中计算出真实的三维重建通常称为SFM（从运动中恢复结构）。

计算重建：

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

3.1 稳健估计基础矩阵

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

代码添加到sfm.py文档： compute_fundamental_normalized()函数将图像点归一化为零均值固定方差

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

3.2 三维重建

# -*- coding: utf-8 -*-
from PIL import Image
from numpy import *
from pylab import *
import numpy as np
import homography
import sfm
import sift

# 标定矩阵
K = array([[2394,0,932],[0,2398,628],[0,0,1]])
# 载入图像，并计算特征
im1 = array(Image.open('picture/home/1.jpg'))
im2 = array(Image.open('picture/home/2.jpg'))
sift.process_image('picture/home/1.jpg','im1.sift')
l0, d0 = sift.read_features_from_file('im1.sift')
sift.process_image('picture/home/2.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()

ransac算法算出的模型参数：

分析：

二次投影后的点和原始特征位置不完全匹配， 但是相当接近。可以进一步调整照相机矩阵来提高重建和二次投影的性能。

（四）问题求助

由于本人电脑的一些问题，本章部分实验结果系实验室师姐的实验数据。问题解决后，再更新自己的。 在这两章的代码编写运行在pycharm2020.1上的python3.7遇见。sift算法运行不了，由于版权所限，于是下载了python2.7.18，加载进pycharm中时发现改不了路径，默认环境依然是python37，出现环境交叉。试过一些网上的19年，18年的在pycharm 同时运行python27和36的方法，失败，出现PIL和pylab无法下载安装在python27当中，需要更新pip

当我更新pip时却显示下载成功但是无法更新的情况，这是什么限制呢？

这个路径如何成功修改？

或者怎样在新版本pycharm里面完好运行python27和37？建虚拟环境也没能达到效果…