用Python和NumPy手把手实现八点法：从匹配点到3D坐标的完整流程

用Python和NumPy手把手实现八点法：从匹配点到3D坐标的完整流程

计算机视觉中的三维重建技术正变得越来越普及，从无人机测绘到增强现实应用，Structure from Motion（SfM）都是核心技术之一。本文将带你用Python和NumPy从零实现八点法计算基础矩阵（F矩阵），并完成三角测量获取3D坐标的完整流程。不同于理论讲解，我们聚焦于可运行的代码实现和实际工程中的关键细节，适合有一定线性代数基础但刚接触SfM实践的开发者。

1. 环境准备与数据加载

在开始前，确保你的Python环境已安装以下库：

pip install numpy opencv-python matplotlib

我们将使用NumPy进行矩阵运算，OpenCV提供一些辅助功能，Matplotlib用于可视化。假设你已经通过SIFT、ORB等特征提取算法获得了两张图片的匹配点对，这些点对应该存储为N×2的数组形式：

import numpy as np # 示例数据：8对匹配点（实际应用中建议使用更多点） points1 = np.array([ # 第一张图片中的点 [488.1, 363.8], [381.2, 371.6], [276.3, 398.7], [184.5, 421.9], [502.4, 215.7], [396.8, 223.4], [294.2, 251.6], [205.3, 275.8] ]) points2 = np.array([ # 第二张图片中的对应点 [521.3, 337.5], [413.9, 345.2], [308.7, 372.4], [217.6, 395.7], [536.8, 189.2], [431.5, 196.8], [329.4, 225.3], [241.2, 249.7] ])

注意：实际项目中，建议使用至少15-20对高质量匹配点以提高稳定性。可以使用OpenCV的特征匹配工具获取这些数据。

2. 数据归一化：提升数值稳定性

直接使用像素坐标计算基础矩阵会导致数值不稳定问题。数据归一化是八点法实现中容易被忽视但至关重要的步骤：

def normalize_points(points): """ 将点集归一化为零均值和单位方差 """ mean = np.mean(points, axis=0) std = np.std(points, axis=0) # 防止除零（齐次坐标的尺度项设为1） std = np.where(std < 1e-8, 1.0, std) T = np.array([ [1/std[0], 0, -mean[0]/std[0]], [0, 1/std[1], -mean[1]/std[1]], [0, 0, 1] ]) # 转换为齐次坐标并归一化 homogeneous = np.column_stack((points, np.ones(len(points)))) normalized = (T @ homogeneous.T).T[:, :2] return normalized, T # 归一化两组点 norm_points1, T1 = normalize_points(points1) norm_points2, T2 = normalize_points(points2)

归一化后的点坐标通常在[-1, 1]范围内，这能显著改善后续SVD分解的数值稳定性。记住这两个变换矩阵T1和T2，我们最后需要它们来恢复原始坐标系下的基础矩阵。

3. 构建线性系统求解基础矩阵

基础矩阵F满足对极约束：x₂ᵀFx₁=0。对于每对匹配点，我们可以展开得到一个线性方程：

def build_A_matrix(points1, points2): """ 构建线性系统矩阵A """ A = [] for (x1, y1), (x2, y2) in zip(points1, points2): A.append([x2*x1, x2*y1, x2, y2*x1, y2*y1, y2, x1, y1, 1]) return np.array(A) A = build_A_matrix(norm_points1, norm_points2)

现在我们需要解方程组Af=0。由于F矩阵具有尺度不确定性（乘以任意非零常数仍满足方程），我们通过SVD求最小二乘解：

def compute_fundamental_matrix(A): """ 通过SVD计算基础矩阵 """ U, S, Vt = np.linalg.svd(A) F = Vt[-1].reshape(3, 3) # 强制秩为2约束 U, S, Vt = np.linalg.svd(F) S[2] = 0 # 将第三个奇异值设为0 F = U @ np.diag(S) @ Vt return F F_norm = compute_fundamental_matrix(A)

得到归一化坐标系下的F_norm后，需要转换回原始像素坐标系：

# 反归一化得到最终基础矩阵 F = T2.T @ F_norm @ T1 F = F / F[2, 2] # 将F[2,2]归一化为1

4. 基础矩阵验证与优化

计算得到的基础矩阵质量如何？我们可以通过重投影误差来验证：

def calculate_epipolar_error(F, points1, points2): """ 计算对极几何误差 """ homogeneous1 = np.column_stack((points1, np.ones(len(points1)))) homogeneous2 = np.column_stack((points2, np.ones(len(points2)))) lines2 = F @ homogeneous1.T # 第二幅图中的极线 lines1 = F.T @ homogeneous2.T # 第一幅图中的极线 # 计算点到极线的距离 error = 0 for i in range(len(points1)): x2, y2, _ = homogeneous2[i] a, b, c = lines2[:, i] error += abs(a*x2 + b*y2 + c) / np.sqrt(a**2 + b**2) x1, y1, _ = homogeneous1[i] a, b, c = lines1[:, i] error += abs(a*x1 + b*y1 + c) / np.sqrt(a**2 + b**2) return error / (2 * len(points1)) print(f"平均对极误差: {calculate_epipolar_error(F, points1, points2):.4f} 像素")

如果误差较大（如>1像素），可以考虑：

1. 使用RANSAC去除异常匹配点
2. 增加匹配点数量（16-20对以上）
3. 检查特征匹配质量

5. 从基础矩阵到相机矩阵

有了基础矩阵F，我们可以推导出两个视图间的相对相机姿态。假设第一个相机的投影矩阵P₁=[I|0]，则第二个相机的P₂可以从F推导：

def compute_camera_matrices(F): """ 从F计算相机矩阵P2（假设P1=[I|0]） """ # 计算本质矩阵E（需要相机内参K，这里假设为单位矩阵） E = F # 当K为单位矩阵时，E=F # 对E进行SVD分解 U, S, Vt = np.linalg.svd(E) W = np.array([[0, -1, 0], [1, 0, 0], [0, 0, 1]]) # P2的四种可能解 P2s = [ np.column_stack((U @ W @ Vt, U[:, 2])), np.column_stack((U @ W @ Vt, -U[:, 2])), np.column_stack((U @ W.T @ Vt, U[:, 2])), np.column_stack((U @ W.T @ Vt, -U[:, 2])) ] # 需要选择正确的解（通常通过三角测量验证） return P2s P2_candidates = compute_camera_matrices(F)

实际应用中，我们需要通过三角测量验证哪个P2是正确的。此外，如果知道相机内参K，应该先计算本质矩阵E=K₂ᵀFK₁。

6. 三角测量获取3D坐标

有了两个相机的投影矩阵P₁和P₂，我们可以通过线性三角测量法恢复3D点：

def linear_triangulation(P1, P2, point1, point2): """ 线性三角测量 """ A = np.array([ point1[0] * P1[2] - P1[0], point1[1] * P1[2] - P1[1], point2[0] * P2[2] - P2[0], point2[1] * P2[2] - P2[1] ]) # SVD求解 _, _, Vt = np.linalg.svd(A) X = Vt[-1] X = X / X[3] # 齐次坐标归一化 return X[:3] # 返回3D坐标 # 假设我们已经确定了正确的P2 P1 = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]]) P2 = P2_candidates[0] # 实际中需要通过验证选择 # 对所有匹配点进行三角测量 points3d = [] for pt1, pt2 in zip(points1, points2): X = linear_triangulation(P1, P2, pt1, pt2) points3d.append(X) points3d = np.array(points3d)

7. 结果可视化与验证

最后，我们可以可视化重建的3D点云：

import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D fig = plt.figure(figsize=(10, 7)) ax = fig.add_subplot(111, projection='3d') ax.scatter(points3d[:, 0], points3d[:, 1], points3d[:, 2], c='b', marker='o') ax.set_xlabel('X Axis') ax.set_ylabel('Y Axis') ax.set_zlabel('Z Axis') ax.set_title('Reconstructed 3D Points') plt.show()

为了验证重建质量，可以计算重投影误差：

def project_points(P, points3d): """ 将3D点投影到图像平面 """ homogeneous = np.column_stack((points3d, np.ones(len(points3d)))) projected = (P @ homogeneous.T).T projected = projected[:, :2] / projected[:, [2]] # 齐次坐标归一化 return projected projected1 = project_points(P1, points3d) projected2 = project_points(P2, points3d) error1 = np.mean(np.linalg.norm(points1 - projected1, axis=1)) error2 = np.mean(np.linalg.norm(points2 - projected2, axis=1)) print(f"重投影误差 - 视图1: {error1:.2f} 像素") print(f"重投影误差 - 视图2: {error2:.2f} 像素")

如果重投影误差较大（如>2像素），可能需要检查特征匹配质量或尝试非线性优化方法进一步优化3D点位置。