SfM实战:从八点法到三角测量,构建三维视觉的几何基石
1. 对极几何理解三维视觉的钥匙想象你站在街边用手机拍下两栋相邻的建筑。当你移动几步再拍一张时这两张照片之间隐藏着丰富的三维信息——这就是对极几何要解决的问题。作为SfM运动恢复结构的核心数学工具它建立了不同视角下图像点的对应关系。对极几何中最关键的数学表达是基础矩阵F。这个3×3的矩阵看似简单却蕴含了两个相机之间的全部几何关系。用生活场景类比假设你闭上一只眼睛用手指对准远处的路灯然后换另一只眼睛观察——你会发现手指似乎跳到了路灯旁边。这种视差现象正是基础矩阵描述的几何约束。在实际操作中我们通过匹配两幅图像中的特征点如SIFT或ORB特征来计算F。每对匹配点满足方程x₂ᵀFx₁0其中x₁和x₂是齐次坐标。展开这个方程会得到一个线性约束当收集到足够多的匹配点时就能解出F矩阵的具体数值。提示齐次坐标在三维视觉中至关重要它将二维像素点(x,y)表示为(x,y,1)使得投影变换可以用统一的矩阵乘法表示。2. 八点法从理论到代码实现2.1 线性方程组的构建八点法之所以得名是因为至少需要8对匹配点才能求解基础矩阵。实际操作中我通常会收集20-30对匹配点来提高鲁棒性。每对点产生一个方程最终形成形如Af0的线性系统# 构建系数矩阵A的Python示例 def build_A_matrix(points1, points2): A [] for (x1, y1), (x2, y2) in zip(points1, points2): row [x2*x1, x2*y1, x2, y2*x1, y2*y1, y2, x1, y1, 1] A.append(row) return np.array(A)这个矩阵的构造有个容易踩坑的地方坐标顺序。我在第一次实现时就曾把x1和x2的顺序弄反导致计算结果完全错误。正确的顺序应该保证方程[x₂x₁ x₂y₁ x₂ y₂x₁ y₂y₁ y₂ x₁ y₁ 1]·f0。2.2 SVD求解与秩2约束得到A矩阵后我们通过奇异值分解(SVD)来求解U, S, Vt np.linalg.svd(A) F Vt[-1].reshape(3, 3) # 取最小奇异值对应的向量但这样得到的F矩阵通常不满足秩为2的条件。解决方法是对F再做一次SVD将第三个奇异值置零U, S, Vt np.linalg.svd(F) S[2] 0 # 强制秩为2 F U np.diag(S) Vt这个步骤非常关键。我曾在无人机视觉定位项目中忽略这一点导致后续的相机姿态估计产生严重漂移。秩2约束保证了基础矩阵确实描述了两个相机之间的本质几何关系。3. 标准化八点法提升精度的秘密武器3.1 数据标准化的必要性原始八点法对噪声非常敏感。通过实验发现当图像坐标范围在[0,1000]像素时直接计算的F矩阵误差可能达到10%以上。这是因为数值计算中大数值会主导小数值。标准化操作包括两个步骤将点集平移至原点减去均值缩放使平均距离为√2除以标准差def normalize_points(points): mean np.mean(points[:, :2], axis0) std np.std(points[:, :2], axis0) T np.array([ [1/std[0], 0, -mean[0]/std[0]], [0, 1/std[1], -mean[1]/std[1]], [0, 0, 1] ]) return T points.T, T3.2 完整的标准化流程实际实现时需要特别注意对两幅图像的点分别进行标准化计算标准化后的F矩阵反标准化得到最终的F# 标准化八点法完整实现 def fundamental_matrix_8point(points1, points2): # 标准化 points1_norm, T1 normalize_points(points1) points2_norm, T2 normalize_points(points2) # 构建A矩阵并求解 A build_A_matrix(points1_norm.T, points2_norm.T) U, S, Vt np.linalg.svd(A) F Vt[-1].reshape(3, 3) # 秩2约束 U, S, Vt np.linalg.svd(F) S[2] 0 F_norm U np.diag(S) Vt # 反标准化 F T2.T F_norm T1 return F / F[2,2] # 归一化在AR眼镜开发项目中采用标准化八点法将重投影误差降低了约60%。这个改进对于保证虚拟物体稳定停留在真实场景中至关重要。4. 三角测量从二维到三维的跨越4.1 相机矩阵与投影几何有了基础矩阵和相机内参我们可以恢复相机的外参矩阵[R|t]。假设第一个相机位于世界坐标系原点其投影矩阵为P₁K[I|0]第二个相机的投影矩阵则为P₂K[R|t]。在三维重建中我们需要求解的是三维点X的坐标使得 P₁X x₁ P₂X x₂这个看似简单的方程组实际上包含深度信息的丢失——我们无法直接从单幅图像确定物体的距离。4.2 线性三角测量实现最直接的解法是构建线性方程组def triangulate(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] ]) _, _, Vt np.linalg.svd(A) X Vt[-1] return X / X[3] # 齐次坐标转非齐次但这种方法对噪声敏感。在实际的SLAM系统中我通常会结合多帧观测和优化算法来提高精度。4.3 非线性优化提升精度更鲁棒的方法是先用线性解法获得初始值再用重投影误差作为代价函数进行优化from scipy.optimize import least_squares def reprojection_error(X, P1, P2, x1, x2): proj1 P1 X proj2 P2 X error np.concatenate([ (proj1[:2]/proj1[2] - x1[:2]), (proj2[:2]/proj2[2] - x2[:2]) ]) return error # 初始估计 X_initial triangulate(P1, P2, x1, x2) # 非线性优化 res least_squares(reprojection_error, X_initial, args(P1, P2, x1, x2)) X_optimized res.x在工业零件三维检测系统中这种优化方法将测量精度从±2mm提升到了±0.5mm大幅提高了质检的可靠性。5. 实战中的挑战与解决方案5.1 异常值处理真实场景中特征匹配总会包含错误匹配。RANSAC算法是解决这个问题的利器def ransac_fundamental(points1, points2, n_iter1000, threshold0.01): best_inliers [] best_F None for _ in range(n_iter): # 随机选择8对点 indices np.random.choice(len(points1), 8, replaceFalse) sample1 points1[indices] sample2 points2[indices] # 计算F矩阵 F fundamental_matrix_8point(sample1, sample2) # 计算误差 errors compute_sampson_distance(F, points1, points2) inliers np.where(errors threshold)[0] if len(inliers) len(best_inliers): best_inliers inliers best_F F # 用所有内点重新计算F return fundamental_matrix_8point(points1[best_inliers], points2[best_inliers])5.2 尺度不确定性三角测量面临的一个根本问题是尺度不确定性。在视觉里程计中我通常通过以下方法解决使用IMU提供初始尺度估计利用场景中已知尺寸的物体在SLAM系统中通过闭环检测校正尺度漂移5.3 退化配置当所有三维点共面或相机纯旋转时八点法和三角测量会失效。这时需要检测平面场景通过单应性矩阵结合其他传感器数据使用基于平面假设的算法替代在开发室内导航系统时我们就遇到过墙面纹理单一导致的退化问题。最终通过结合深度相机数据解决了这一挑战。