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

最新推荐文章于 2024-09-02 22:55:35 发布

鲸~屿

最新推荐文章于 2024-09-02 22:55:35 发布

阅读量487

点赞数 10

文章标签：计算机视觉人工智能

本文链接：https://blog.csdn.net/qq_53034510/article/details/141752991

版权

5.1 外极几何

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

三维场景点X经过4×4的单应性矩阵H变换为HX后，则HX在照相机 $PH^{-1}$ 里得到的图像点和X在照相机 P 里得到的图像点相同。利用照相机方程为 $\lambda x=PX=PH^{-1}HX=\hat{P}\hat{X}$ ，这里的单应性矩阵通常只做刚体变换，即只通过单应矩阵变换了坐标系。将原点和坐标轴与第一个照相机对齐，则：

$P_{1}=K_{1}[I|0]$ 和 $P_{2}=K_{2}[R|t]$

其中 $K_{1}$ 和 $K_{2}$ 是标定矩阵，R 是第二个照相机的旋转矩阵，t是第二个照相机的平移向量。利用这些照相机参数矩阵，我们可以找到点X的投影点 $x_{1}$ 和 $x_{2}$ （分别对应于投影矩阵 $P_{1}$ 和 $P_{2}$ ）的关系。

同一个图像点经过不同的投影矩阵产生的不同投影点必须满足： $x_{2}^{T}Fx_{1}=0$ ，其中F的表达式为 $F=K_{2}^{-T}S_{t}RK_{1}^{-1}$ 。矩阵 $S_{t}$ 为反对称矩阵：

$S_{t}=\begin{bmatrix} 0 & -t_{3} & t_{2}\\ t_{3} & 0 & -t_{1}\\ -t_{2} & t_{1} & 0 \end{bmatrix}$

矩阵F为基础矩阵。基础矩阵可以由两照相机的参数矩阵（相对旋转R和平移t）表示。由于det(F)=0，所以基础矩阵的秩小于等于2。

第二个视图中的图像点 $x_{2}$ ，利用上述公式，我们找到对应第一幅图像的一条直线：

$x_{2}^{T}Fx_{1}=l_{1}^{T}x_{1}=0$

其中 $l_{1}^{T}x_{1}=0$ 是第一幅图像中的一条直线。这条线称为对应于点 $x_{2}$ 的外极线，也就是说，点 $x_{2}$ 在第一幅图像中的对应点一定在这条线上。

5.1.1 一个简单的数据集

我们需要一个带有图像点、三维点和照相机参数矩阵的数据集。我们这里使用一个牛津多视图数据集。

import numpy as np
from PIL import Image
from numpy import genfromtxt
import camera
import matplotlib.pyplot as plt

# 载入一些图像，并确保关闭图像对象
def load_image(file_path):
    img = Image.open(file_path)
    arr = np.array(img)
    img.close()  # 关闭图像对象
    return arr

im1 = load_image('images001.jpg')
im2 = load_image('images002.jpg')

# 载入每个视图的二维点到列表中
points2D = []
for i in range(3):
    points_file = f'00{i + 1}.corners'
    try:
        points2D.append(np.loadtxt(points_file).T)
    except Exception as e:
        print(f"Error loading {points_file}: {e}")

# 载入三维点
try:
    points3D = np.loadtxt('p3d').T
except Exception as e:
    print(f"Error loading 3D points: {e}")

# 载入对应关系
try:
    corr = genfromtxt('nview-corners', dtype='int', missing='*')
except Exception as e:
    print(f"Error loading correspondence data: {e}")

# 载入照相机矩阵到 Camera 对象列表中
P = []
for i in range(3):
    P_file = f'00{i + 1}.P'
    try:
        P.append(camera.Camera(np.loadtxt(P_file)))
    except Exception as e:
        print(f"Error loading camera matrix {P_file}: {e}")

将三维的点投影到一个视图，然后和观测到的图像点比较：

# 将三维点转换成齐次坐标表示，并投影
X = np.vstack((points3D, np.ones(points3D.shape[1])))
x = P[0].project(X)
# 在视图1中绘制点
plt.figure()
plt.imshow(im1)
plt.plot(points2D[0][0], points2D[0][1], '*')
plt.axis('off')
plt.figure()
plt.imshow(im1)
plt.plot(x[0], x[1], 'r.')
plt.axis('off')
plt.show()

运行结果如下

5.1.2 用Matplotlib绘制三维数据

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

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np

# 创建第一个3D图，绘制测试数据
fig1 = plt.figure()
ax1 = fig1.add_subplot(111, projection="3d")

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

# 在三维中绘制点
ax1.plot(X.flatten(), Y.flatten(), Z.flatten(), 'o')
ax1.set_title("Test Data Points")
plt.show()

# 创建第二个3D图，绘制随机生成的点
points3D = [np.random.rand(10), np.random.rand(10), np.random.rand(10)]

fig2 = plt.figure()
ax2 = fig2.add_subplot(111, projection='3d')
ax2.plot(points3D[0], points3D[1], points3D[2], 'k.')
ax2.set_title("Random 3D Points")
plt.show()

5.1.3 计算F：八点法

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

$\begin{bmatrix} x_{2}^{1}x_{1}^{1}& x_{2}^{1}y_{1}^{1} & x_{2}^{1}w_{1}^{1} & ... & w_{2}^{1}w_{1}^{1}\\ x_{2}^{2}x_{1}^{2}& x_{2}^{2}y_{1}^{2} & x_{2}^{2}w_{1}^{2} & ... & w_{2}^{2}w_{1}^{2} \\ .& . & . & & . \\ . & . & . & & . \\ .& . & . & & . \\ x_{2}^{n}x_{1}^{n}& x_{2}^{n}y_{1}^{n} & x_{2}^{n}w_{1}^{n} & ... & w_{2}^{n}w_{1}^{n} \end{bmatrix}\begin{bmatrix} F_{11}\\ F_{12}\\ F_{13}\\ .\\ .\\ F_{33} \end{bmatrix}=Af=0$

其中，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，所以该算法叫做八点法。以下代码是八点法中最小化 $||Af||$ 的函数：

import numpy as np

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 = np.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 = np.linalg.svd(A)
    F = V[-1].reshape(3, 3)
    # 受限F
    # 通过将最后一个奇异值置0，使秩为2
    U, S, V = np.linalg.svd(F)
    S[2] = 0
    F = np.dot(U, np.dot(np.diag(S), V))
    return F

上面的函数忽略了一个重要的步骤：对图像坐标进行归一化，这可能会带来数值问题。我们会在后面加以解决。

5.1.4 外极点和外极线

外极点满足Fe1 =0，因此可以通过计算F的零空间来得到。计算的函数为

def compute_epipole(F):
    """ 从基础矩阵F中计算右极点（可以使用F.T获得左极点）"""
    # 返回F的零空间（Fx = 0）
    U, S, V = np.linalg.svd(F)
    e = V[-1]
    return e / e[2]

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

import numpy as np
import Merton1
from PIL import Image
import matplotlib.pyplot as plt
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 = np.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 = np.linalg.svd(A)
    F = V[-1].reshape(3, 3)
    # 受限F
    # 通过将最后一个奇异值置0，使秩为2
    U, S, V = np.linalg.svd(F)
    S[2] = 0
    F = np.dot(U, np.dot(np.diag(S), V))
    return F

def compute_epipole(F):
    """ 从基础矩阵F中计算右极点（可以使用F.T获得左极点）"""
    # 返回F的零空间（Fx = 0）
    U, S, V = np.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 = np.dot(F,x)
    #外极线参数和值
    t = np.linspace(0,n,100)
    lt = np.array([(line[2]+line[0]*tt)/(-line[1]) for tt in t])
    #仅仅处理位于图像内部的点和线
    ndx = (lt >=0) & (lt < m)
    plt.plot(t[ndx],lt[ndx],linewidth=2)

    if show_epipole:
        if epipole is None:
            epipole = compute_epipole(F)
        plt.plot(epipole[0]/epipole[2],epipole[1]/epipole[2],'r*')

corr = np.genfromtxt('nview-corners',dtype='int')
#在前面两个视图中点的索引
ndx = (corr[:,0]>=0) & (corr[:,1]>=0)

#获得坐标，并将其用齐次坐标表示
x1 = Merton1.points2D[0][:,Merton1.corr[ndx,0]]
x1 = np.vstack( (x1,np.ones(x1.shape[1])) )
x2 = Merton1.points2D[1][:,Merton1.corr[ndx,1]]
x2 = np.vstack( (x2,np.ones(x2.shape[1])) )

#计算F
F = compute_fundamental(x1,x2)
#计算极点
e = compute_epipole(F)
#绘制图像
plt.figure()
plt.imshow(Merton1.im1)
#分别绘制每条线，这样会绘制出很漂亮的颜色
for i in range(5):
    plot_epipolar_line(Merton1.im1,F,x2[:,i],e,False)
plt.axis('off')

plt.figure()
plt.imshow(Merton1.im2)
#分别绘制每个点，这样会绘制出和线同样的颜色
for i in range(5):
    plt.plot(x2[0,i],x2[1,i],'o')
plt.axis('off')
plt.show()

运行结果如图所示：

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

5.2.1 三角剖分

对于两个照相机 $P_{1}$ 和 $P_{2}$ 的视图，三维实物点X的投影点为 $x_{1}$ 和 $x_{2}$ （这里用齐次坐标表示），照相机方程定义了下列关系：

$\begin{bmatrix} P_{1} & -x_{1} & 0 \\ P_{2} & 0 & -x_{2} \end{bmatrix}\begin{bmatrix} X\\ \lambda _{1}\\ \lambda _{2} \end{bmatrix}=0$

我们可以通过SVD算法来得到三维点的最小二乘估值。下面的函数用于计算一个点对的最小二乘三角剖分。

import Merton1
import numpy as np
import matplotlib.pyplot as plt
import sfm
from mpl_toolkits.mplot3d import axes3d
 # 前两个视图中点的索引
ndx = (Merton1.corr[:,0]>=0) & (Merton1.corr[:,1]>=0)
 # 获取坐标，并用齐次坐标表示
x1 = Merton1.points2D[0][:,Merton1.corr[ndx,0]]
x1 = np.vstack( (x1,np.ones(x1.shape[1])) )
x2 = Merton1.points2D[1][:,Merton1.corr[ndx,1]]
x2 = np.vstack( (x2,np.ones(x2.shape[1])) )
Xtrue = Merton1.points3D[:,ndx]
Xtrue = np.vstack( (Xtrue,np.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 = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(Xest[0],Xest[1],Xest[2],'ko')
ax.plot(Xtrue[0],Xtrue[1],Xtrue[2],'r.')
plt.axis('equal')
plt.show()

上面的代码首先利用前两个视图的信息来对图像点进行三角剖分，然后把前三个图像点的齐次坐标输出到控制台，最后绘制出恢复的最接近三维图像点。输出到控制台的信息如下所示：

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

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

每个三维点 $\mathbf{X}_{i}$ （齐次坐标系下）按照 $\lambda_{i} \mathbf{x}_{i}=\mathbf{PX_{i}}$ 投影到图像点 $\mathbf{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$

其中 $\mathbf{p}_{1}$ 、 $\mathbf{p}_{2}$ 和 $\mathbf{p}_{3}$ 是矩阵 $\mathbf{P}$ 的三行。如下所示

$\mathbf{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 = np.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 = np.linalg.svd(M)
    return V[-1,:12].reshape((3,4))

import sfm, camera
import Merton1
import numpy as np
import matplotlib.pyplot as plt
corr = Merton1.corr[:,0] # 视图 1
ndx3D = np.where(corr>=0)[0] # 丢失的数值为-1
ndx2D = corr[ndx3D]
 # 选取可见点，并用齐次坐标表示
x = Merton1.points2D[0][:,ndx2D] # 视图 1
x = np.vstack( (x,np.ones(x.shape[1])) )
X = Merton1.points3D[:,ndx3D]
X = np.vstack( (X,np.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)
# 绘制图像
plt.figure()
plt.imshow(im1)
plt.plot(x[0], x[1], 'bo')
plt.plot(xest[0], xest[1], 'r.')
plt.axis('off')
plt.show()

实验效果如下：

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

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

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

在没有任何照相机内参数知识的情况下，照相机矩阵只能通过射影变换恢复出来。也就是说，如果利用照相机的信息来重建三维点，那么该重建只能由射影变换计算出来（你可以得到整个投影场景中无畸变的重建点）

在无标定的情况下，第二个照相机矩阵可以使用一个（3×3）的射影变换得出。一个简单的方法是：

$P_{2}=[S_{e}F|e]$

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

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

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

$\mathbf{x}_{K}=K^{-1}K[R|\mathbf{t}]\mathbf{X}=[R|\mathbf{t}]\mathbf{X}$

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

$\mathbf{x}_{K_2}^TF\mathbf{x}_{K_1}=0$

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

5.3 多视图重建

让我们来看，如何使用上面的理论从多幅图像中计算出真实的三维重建。由于照相机的运动给我们提供了三维结构，所以这样计算三维重建的方法通常称为SfM（从运动中恢复结构）。

假设照相机已经标定，计算重建可以分为下面4个步骤：

        (1) 检测特征点，然后在两幅图像间匹配；

        (2) 由匹配计算基础矩阵；

        (3) 由基础矩阵计算照相机矩阵；

        (4) 三角剖分这些三维点。

5.3.1 稳健估计基础矩阵

类似于稳健计算单应性矩阵（3.3节），当存在噪声和不正确的匹配时，我们需要估计基础矩阵。和前面的方法一样，我们使用RANSAC方法，只不过这次结合了八点算法。

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 = np.dot(F,x1)
        Fx2 = np.dot(F,x2)
        denom = Fx1[0]**2 + Fx1[1]**2 + Fx2[0]**2 + Fx2[1]**2
        err = ( np.diag(np.dot(x1.T,np.dot(F,x2))) )**2 / denom
        # 返回每个点的误差
        return err

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 = np.mean(x1[:2],axis=1)
    S1 = np.sqrt(2) / np.std(x1[:2])
    T1 = np.array([[S1,0,-S1*mean_1[0]],[0,S1,-S1*mean_1[1]],[0,0,1]])
    x1 = np.dot(T1,x1)
    x2 = x2 / x2[2]
    mean_2 = np.mean(x2[:2],axis=1)
    S2 = np.sqrt(2) / np.std(x2[:2])
    T2 = np.array([[S2,0,-S2*mean_2[0]],[0,S2,-S2*mean_2[1]],[0,0,1]])
    x2 = np.dot(T2,x2)
    #使用归一化的坐标计算F
    F = compute_fundamental(x1,x2)
    #反归一化
    F = np.dot(T1.T,np.dot(F,T2))
    return F/F[2,2]

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

def F_from_ransac(x1,x2,model,maxiter=5000,match_theshold=1e-6):
    """ 使用 RANSAN 方法（ransac.py，来自http://www.scipy.org/Cookbook/RANSAC），
    从点对应中稳健地估计基础矩阵F
    输入：使用齐次坐标表示的点x1，x2（3×n的数组）"""
    import ransac
    data = np.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']

这里，我们返回最佳基础矩阵F，以及正确点的索引，这样可以知道哪些匹配和F 矩阵是一致的。

5.3.2 三维重建示例

import homography
import sfm
import sift2
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import camera

# 标定矩阵
K = np.array([[2394,0,932],[0,2398,628],[0,0,1]])
 # 载入图像，并计算特征
im1 = np.array(Image.open('alcatraz1.jpg'))
sift2.process_image('alcatraz1.jpg','im1.sift')
l1,d1 = sift2.read_features_from_file('im1.sift')
im2 = np.array(Image.open('alcatraz2.jpg'))
sift2.process_image('alcatraz2.jpg','im2.sift')
l2,d2 = sift2.read_features_from_file('im2.sift')
 # 匹配特征
matches = sift2.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 = np.dot(np.linalg.inv(K),x1)
x2n = np.dot(np.linalg.inv(K),x2)
# 使用RANSAC方法估计E
model = sfm.RansacModel()
E,inliers = sfm.F_from_ransac(x1n,x2n,model)
 # 计算照相机矩阵（P2是4个解的列表）
P1 = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0]])
P2 = sfm.compute_P_from_essential(E)

# 选取点在照相机前的解
ind = 0
maxres = 0
infront = 0
for i in range(4):
    # 三角剖分正确点，并计算每个照相机的深度
    X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2[i])
    d1 = np.dot(P1, X)[2]
    d2 = np.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 = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot(-X[0], X[1], X[2], 'k.')
    plt.axis('off')
# 绘制X的投影

# 绘制三维点
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2[ind])
x1p = cam1.project(X)
x2p = cam2.project(X)
# 反K归一化
x1p = np.dot(K, x1p)
x2p = np.dot(K, x2p)
plt.figure()
plt.imshow(im1)
plt.gray()
plt.plot(x1p[0], x1p[1], 'o')
plt.plot(x1[0], x1[1], 'r.')
plt.axis('off')
plt.figure()
plt.imshow(im2)
plt.gray()
plt.plot(x2p[0], x2p[1], 'o')
plt.plot(x2[0], x2[1], 'r.')
plt.axis('off')
plt.show()

5.3.3 多视图的扩展示例

1.多视图

当我们有同一场景的多个视图时，三维重建会变得更准确，包括更多的细节信息。

对于视频序列，我们可以使用时序关系，在连续的帧对中匹配特征。相对方位需要从每个帧对增量地加入下一个帧对，该方法非常有效，同时可以使用跟踪有效地找到对应点。一个问题是误差会在加入的视图中积累。

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

2.光束法平差

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

3.自标定

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

鲸~屿

关注

10
点赞
踩
19

收藏

觉得还不错? 一键收藏
0
评论
计算机视觉第五章多视图几何

如果有一个场景的两个视图以及视图中的对应图像点，那么根据照相机间的空间相对位置关系、照相机的性质以及三维场景点的位置，可以得到对这些图像点的一些几何关系约束。以三维重建为例，恢复出的点的位置和由估计的基础矩阵计算出的照相机矩阵都存在误差。根据在照相机标定矩阵参数上做出的假设，以及可用的图像数据的类型（特征匹配、平行线、平面等），我们有很多不同的自标定算法。上面的代码首先利用前两个视图的信息来对图像点进行三角剖分，然后把前三个图像点的齐次坐标输出到控制台，最后绘制出恢复的最接近三维图像点。
复制链接

扫一扫