计算机视觉编程多视图几何

_Taylor

已于 2023-08-20 12:26:56 修改

阅读量307

点赞数

分类专栏：计算机视觉编程文章标签：计算机视觉人工智能

于 2023-08-20 12:16:42 首次发布

本文链接：https://blog.csdn.net/sketch_2314/article/details/132321540

版权

计算机视觉编程专栏收录该内容

7 篇文章 0 订阅

订阅专栏

第五章多视图几何

前言
5.1 外极几何
5.2 照相机和三维结构的计算
5.3 多视图重建
- 5.3.1 稳健估计基础矩阵
- 5.3.2 三维重建示例

前言

本专栏按《python计算机视觉编程 ——Jan Erik Solem》一书为参考，第五章介绍用多个视图来三维重建，即恢复相机位置和三维结构。主要包含特征匹配、估计基础矩阵和照相机矩阵

5.1 外极几何

多视图几何用于研究和处理从不同视角获取的图像或图像序列之间的几何关系，它涉及如何在多个视角下理解和分析场景，以及如何从多个图像中恢复三维场景的信息。外极几何是多视图几何中的一个关键概念，描述了两个不同摄像机（或视角）之间的关系，通过连接摄像机的光学中心、图像平面上的点和对应点的位置来建立。通过外极几何有助于在多个图像中寻找匹配点、恢复三维场景信息以及进行深度估计等任务

在外极几何中，主要涉及以下概念：

外极线：对于在一个图像中的点，与之对应的外极线是在另一个图像上的一条线，它连接了摄像机光学中心和对应点的投影。这有助于在一个图像中找到与之对应的点，从而减少匹配问题的搜索空间
本质矩阵：它是描述两个摄像机之间关系的矩阵，可以通过已知的相机内参数和匹配点来计算。本质矩阵能够帮助我们在不恢复具体摄像机的姿态情况下，从一组图像中恢复场景的三维信息
基础矩阵：类似于本质矩阵，但考虑了图像平面上的坐标。基础矩阵能够在图像平面上描述两个视图之间的几何关系，它满足 $x_2^TFx_1=0$ ，其中 $F$ 为基础矩阵 $x_2$ 和 $x_1$ 表示一个点在两个视图上的投影点。基础矩阵可以用于简化匹配和去除错配特征

请添加图片描述
如上图，三维点P分别投影在两个视图上形成p1和p2，两个照相机的中心间的基线O1和O2与图像相交于外极点e1和e2，p1和e1的连线以及p2和e2的连线就称为外极线

5.1.1 一个简单的数据集

按照书上的操作，从Multi-view Data上下载官方的数据，此处使用Merton College I的相关数据集，使用下述代码进行数据加载

# 载入图像
img1 = array(Image.open('./images/001.jpg'))
img2 = 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_values='*')
# 载入照相机矩阵到Camera对象列表中
P = [camera.Camera(loadtxt('./2D/00' + str(i + 1) + '.P')) for i in range(3)]

通过上述代码就能够加载两幅图像、三个视图中的所有图像特征点、重建后的三维点和照相机参数矩阵。将上述代码打包文件就能使用execfile()命令执行上述代码操作获取数据

# 在视图1中绘制点
figure()
imshow(img1)
plot(points2D[0][0], points2D[0][1], '*')
axis('off')
figure()
imshow(img1)
plot(x[0], x[1], 'r.')
axis('off')
show()

视图1与图像点	视图1与投影的三维点

绘制图像结果如下，能看到对比于左边未投影的图像，右边图像的三维点明显增多，这些点是由试图2和视图3重建出来的

5.1.2 用Matplotlib绘制三维数据

接下来使用Matplotlib可视化重建结果，通过add_axes(Axes3D(fig)函数实现三维绘图，书上对axes对象加上projection='3d'会产生报错，报错解决方法

from mpl_toolkits.mplot3d import axes3d, Axes3D

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

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

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

show()

请添加图片描述
用类似代码画出之前的样本数据，能得到如下的效果

还可以通过ax.view_init(elev=90, azim=180) # 设置视角，elev上下旋转角度、azim 表示水平旋转角度来观察不同视角的三维图像点

上侧面	俯视图

5.1.3 计算F: 八点法

八点法用于两个相机之间的几何校准和视觉定位，它的具体步骤如下：

收集匹配点对：首先需要从两幅图像中收集至少8对特征点，这些点在两个图像中具有对应关系。这些点可以是从特征检测器（如SIFT、SURF等）得到的关键点，然后通过匹配算法找到两幅图像之间的对应点
构建矩阵方程：对于每个点对 $x_i, x_i')$ ，其中 $x_i$ 是第一张图像中的点， $x_i'$ 是第二张图像中的对应点，你可以构建一个方程 $x_i'^T F x_i = 0$ ，其中 $F$ 是基础矩阵。将所有点对的方程组合成一个矩阵方程
求解线性方程组：将步骤2中得到的矩阵方程整理成一个线性方程组，形式为 $A f = 0$ ，其中 $f$ 是基础矩阵 $F$ 的展开向量（ $F$ 是一个 $\times 3$ 的矩阵， $f$ 是它的展开向量，共有9个分量）。然后使用线性代数技术，如奇异值分解（SVD）来求解该线性方程组，从而得到基础矩阵 $F$

由于八点法的求解可能得到不精确的结果，需要对得到的基础矩阵 $F$ 进行约束，以满足基础矩阵的性质。如对 $F$ 进行奇异值分解，然后将最小的奇异值设为零，从而满足基础矩阵的秩为2的约束。下面是使用归一化的八点算法

def compute_fundamental(x1,x2):
    """    Computes the fundamental matrix from corresponding points 
        (x1,x2 3*n arrays) using the 8 point algorithm.
        Each row in the A matrix below is constructed as
        [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.")
    
    # build matrix for equations
    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] ]
            
    # compute linear least square solution
    U,S,V = linalg.svd(A)
    F = V[-1].reshape(3,3)
        
    # constrain F
    # make rank 2 by zeroing out last singular value
    U,S,V = linalg.svd(F)
    S[2] = 0
    F = dot(U,dot(diag(S),V))
    
    return F/F[2,2]

5.1.4 外极点和外极线

接下来绘制外极点和外极线，之前有外极点满足 $Fe_1=0$ 的约束条件，使用以下代码计算，若要求另一副图像的外极点，将 $F$ 转置即可

def compute_epipole(F):
    """ Computes the (right) epipole from a 
        fundamental matrix F. 
        (Use with F.T for left epipole.) """
    
    # return null space of F (Fx=0)
    U,S,V = linalg.svd(F)
    e = V[-1]
    return e/e[2]

在图像中用特征匹配后使用上述函数绘制外极线的代码如下

from 视觉编程 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(img1)

# 分别绘制每条线
for i in range(5):
    sfm.plot_epipolar_line(img1, F, x2[:, i], e, False)
axis('off')

figure()
imshow(img2)
# 分别绘制每个点
for i in range(5):
    plot(x2[0, i], x2[1, i], 'o')
axis('off')

show()

外极点	外极线

最终得到上图的五个外极点，注意图像右侧三个点中红色的被覆盖了，将不同颜色的点和外极线连接起来得到右图，这些线在图像左侧位置交于一点（图中未显示），说明这些点之间的对应可以在另外一副图像中找到

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

5.2.1 三角剖分

三角剖分用于将凸多边形或复杂的多边形区域分割成一组不重叠的三角形，在已知照相机模型的情况下，可以通过三角剖分可以恢复点的三维位置。对于两个相机p1和p2的视图，三维实物点X的投影为x1和x2，则有 $\begin{bmatrix} P_1 & -x_1 & 0 \\ P_2 & 0 & -x_2 \end{bmatrix}\begin{bmatrix} X \\ \lambda_1 \\ \lambda_2 \end{bmatrix}=0$
使用以下代码对原数据进行三角剖分

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

Xture = points3D[:, ndx]
Xture = vstack((Xture, ones(Xture.shape[1])))

# 检查前三个点
Xest = sfm.triangulate(x1, x2, P[0].P, P[1].P)
print(Xest[:, :3])
print(Xture[:, :3])

# 绘制图像
fig = figure()
ax = fig.add_axes(Axes3D(fig))
ax.plot(Xest[0], Xest[1], Xest[2], 'ko')
ax.plot(Xture[0], Xture[1], Xture[2], 'r.')
axis('equal')

show()

齐次坐标	数据点对比

输出的图像点齐次坐标和三维坐标效果如上，估计出来的三维图像点与实际点很接近，其中黑色的表示估计的点，红色的表示真实的点，通过转换坐标系能看的更加清楚

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

通常直接使用线性变换来计算计算照相机矩阵P，本质上是三角剖分方法的逆问题，也称照相机反切法，已知的三维点和它们在图像中的对应点来计算相机矩阵的步骤简要如下：

构建方程：对于每一个三维点 P（在齐次坐标下表示），可以建立一个方程，将其变换到图像坐标系上的对应点 p。这个方程可以表示为：p = M * P，其中 M 是相机矩阵
解方程：将收集到的多个方程组成一个线性方程组，然后通过最小二乘法或其他优化方法来解这个方程组，从而得到相机矩阵 M
分解矩阵：通常，得到的相机矩阵 M 包含内参数矩阵和外参数矩阵的组合。你可以使用矩阵分解技术，如奇异值分解（SVD）来将其分解为相机的内参数和外参数

下面给出在数据集上测试算法的代码

corr = corr[:,0] # 视图1
ndx3D = where(corr >= 0)[0] # 丢失的数值为1
ndx2D = corr[ndx3D]

# 选取可见点并用齐次坐标表示
x = points2D[0][:, ndx2D]
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(img1)
plot(x[0], x[1], 'bo')
plot(xest[0], xest[1], 'r.')
axis('off')

show()

齐次坐标	估计照相机投影点

左侧为照相机矩阵，右侧为计算出的视图1中的投影点，其中真实的用圆圈表示，估计出的投影点用点表示，可以看到估计的与真实的点都非常接近

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

基础矩阵用于描述两个不同视角下的图像之间的几何关系。假设已知一个归一化的相机矩阵P1，使用它计算第二个相机矩阵时，可以区分为未标定和已标定两种情况，其中未标定指相机的内部参数（例如焦距、主点位置、畸变等）未知或者没有被准确测量的情况

未标定的情况——投影重建
使用一个引射变换计算，基本公式为 $P_2=[S_eF|e]$ ，其中e是左极点，满足 $e^TF=0$ ， $S_e$ 为反对称矩阵，具体代码为

def skew(a):
    """ Skew matrix A such that a x v = Av for any v. """

    return array([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])
    
def compute_P_from_fundamental(F):
    """    Computes the second camera matrix (assuming P1 = [I 0]) 
        from a fundamental matrix. """
        
    e = compute_epipole(F.T) # left epipole
    Te = skew(e)
    return vstack((dot(Te,F.T).T,e)).T

标定的情况——度量重建
在已标定的情况下，重建会保持一些度量特性。给定标定矩阵K，将它的逆矩阵作用于图像点，在新的图像坐标系下，照相机方程为 $x_K=K^-1K[R|t]X=[R|t]X$ 其中点满足基础矩阵 $x_{K2}^TFx_{K1}=0$ 。下面是具体代码

def compute_P_from_essential(E):
    """    Computes the second camera matrix (assuming P1 = [I 0]) 
        from an essential matrix. Output is a list of four 
        possible camera matrices. """
    
    # make sure E is rank 2
    U,S,V = svd(E)
    if det(dot(U,V))<0:
        V = -V
    E = dot(U,dot(diag([1,1,0]),V))    
    
    # create matrices (Hartley p 258)
    Z = skew([0,0,-1])
    W = array([[0,-1,0],[1,0,0],[0,0,1]])
    
    # return all four solutions
    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(Structure from Motion 从运动中恢复结构)，在已标定的情况下主要有以下几个步骤

特征提取与匹配：对于每一对图像，使用特征提取算法（如SIFT、SURF、ORB等）来检测和描述图像中的特征点（如角点、边缘等）。然后，通过特征匹配算法将相邻图像中的特征点进行匹配，从而建立特征点之间的对应关系
计算基础矩阵：使用已知的特征点对应关系，计算基础矩阵。基础矩阵描述了两个相机视角之间的几何关系，它包含了视差信息，可以用来恢复相机之间的运动
由基础矩阵计算照相机矩阵：在已经计算得到基础矩阵的情况下，可以使用相机标定参数来计算照相机矩阵。通过分解基础矩阵，可以得到两个相机的投影矩阵，进而得到三维点的重建结果
三角剖分这些三维点：使用三角剖分方法，将重建得到的三维点从稀疏点云转化为一组连续的三角形，以便于后续的可视化、分析和应用

5.3.1 稳健估计基础矩阵

与前面章节类似，使用RANSAC方法，并结合八点算法，主要代码函数如下

def compute_fundamental_normalized(x1,x2):
    """    Computes the fundamental matrix from corresponding points 
        (x1,x2 3*n arrays) using the normalized 8 point algorithm. """

    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")

    # normalize image coordinates
    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)

    # compute F with the normalized coordinates
    F = compute_fundamental(x1,x2)

    # reverse normalization
    F = dot(T1.T,dot(F,T2))

    return F/F[2,2]

除了上面使用八点算法计算基础矩阵函数外，还有使用RANSAC计算基础矩阵和从点对应中稳健估计基础矩阵的类，具体代码请自行下载查看

5.3.2 三维重建示例

图像重建的过程之前已经提及，这里直接给代码

# 加载图像
imname1 = './images/001.jpg'
imname2= './images/002.jpg'

# 特征提取
sift0 = cv2.SIFT_create()

im1 = cv2.imread(imname1)
l1, d1 = sift0.detectAndCompute(im1, None)

im2 = cv2.imread(imname2)
l2, d2 = sift0.detectAndCompute(im2, None)

# 匹配特征，并计算齐次坐标（复制）
matches = cv2.BFMatcher().knnMatch(d1, d2, k=2)
ndx = matches[0]
# 去除错误匹配
matches = [m for m, n in matches if m.distance < 0.7 * n.distance]
ndx = matches

x1 = make_homog(array([l1[m.queryIdx].pt for m in matches]).T)
ndx2 = [m.trainIdx for m in matches]
x2 = make_homog(array([l2[m.trainIdx].pt for m in matches]).T)

x1n = dot(inv(K), x1)
x2n = dot(inv(K), x2)

model = sfm.RansacModel()
E, inliers = sfm.F_from_ransac(x1n, x2n, model)

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]

# 绘制三维图像
fig = figure()
ax = fig.add_axes(Axes3D(fig))
ax.plot(-X[0],X[1],X[2],'k.')
axis("off")

# 绘制X的投影
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2)

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