单应矩阵估计实战：从基础理论到RANSAC算法实现

1次阅读

共计 2251 个字符，预计需要花费 6 分钟才能阅读完成。

在计算机视觉中，单应矩阵（Homography Matrix）扮演着连接两个平面之间对应点的桥梁角色。无论是全景拼接、增强现实还是相机标定，都离不开这个神奇的 3×3 矩阵。然而，实际应用中我们常会遇到匹配点对不准确的问题——噪声和异常值就像捣蛋鬼，让直接计算的单应矩阵变得不可靠。这时候，RANSAC 算法就像一位经验丰富的侦探，能从杂乱的数据中找出真相。

单应矩阵 H 是一个 3×3 的矩阵，表示两个平面之间的投影变换关系。虽然它有 9 个元素，但实际上只有 8 个自由度（DoF），因为通常我们会对其进行归一化（例如令 h33=1）。其数学定义为：

$$
\begin{bmatrix} x’ \ y’ \ 1 \end{bmatrix} = H \begin{bmatrix} x \ y \ 1 \end{bmatrix} = \begin{bmatrix} h_{11} & h_{12} & h_{13} \ h_{21} & h_{22} & h_{23} \ h_{31} & h_{32} & h_{33} \end{bmatrix} \begin{bmatrix} x \ y \ 1 \end{bmatrix}
$$

直接线性变换 (DLT) 是最直观的求解方法，它通过 4 个匹配点对构建方程组来求解 H。但现实很骨感：

特征匹配总是存在误差（哪怕 SIFT/SURF 也不能 100% 准确）
误匹配 (outliers) 可能完全偏离真实变换
噪声会导致 DLT 求得的矩阵严重失真

实验表明，当 outliers 比例超过 15% 时，DLT 的结果就变得不可用了。这就是我们需要 RANSAC 的原因。

随机采样：从匹配点对中随机选取 4 组（求解 H 的最小集）
模型拟合：用这 4 组点计算临时单应矩阵 H_temp
内点判断：用 H_temp 测试所有点对，距离小于阈值的视为内点(inliers)

重复上述过程，最终选择内点最多的 H 作为最优解。

RANSAC 需要运行多少次才能保证成功？这取决于：

数据中内点的比例 w
所需置信度 p（通常取 0.99）
每次采样需要的点数 k（单应矩阵 k =4）

迭代次数 N 的计算公式为：

$$
N = \frac{\log(1-p)}{\log(1-w^k)}
$$

例如，当 w =0.5（内点占 50%）时，要保证 99% 的成功率需要运行 72 次。

import cv2
import numpy as np

# 假设我们已有匹配点对：pts1 和 pts2（均为 N×2 的数组）# 使用 OpenCV 的 RANSAC 实现
H, mask = cv2.findHomography(pts1, pts2, cv2.RANSAC, ransacReprojThreshold=3.0)

# 查看内点数量
inlier_num = np.sum(mask)
print(f"Found {inlier_num} inliers out of {len(pts1)} matches")

# 手动实现 RANSAC（简化版）def ransac_homography(pts1, pts2, threshold=3.0, max_iters=1000):
    best_H = None
    best_inliers = 0

    for _ in range(max_iters):
        # 1. 随机采样 4 对点
        indices = np.random.choice(len(pts1), 4, replace=False)
        sample1 = pts1[indices]
        sample2 = pts2[indices]

        # 2. 计算临时 H（注意 OpenCV 需要转置）H_temp, _ = cv2.findHomography(sample1, sample2, 0)
        if H_temp is None:
            continue

        # 3. 计算所有点的投影误差
        ones = np.ones((len(pts1), 1))
        pts1_h = np.hstack([pts1, ones])
        projected = (H_temp @ pts1_h.T).T
        projected = projected[:, :2] / projected[:, 2:]  # 齐次坐标转笛卡尔

        errors = np.linalg.norm(projected - pts2, axis=1)
        inliers = np.sum(errors < threshold)

        if inliers > best_inliers:
            best_inliers = inliers
            best_H = H_temp

    return best_H