§1.2.10

RANSAC 的内点比例、迭代次数公式 N = log(1-p)/log(1-w^k)？

手写练习

—手写 RANSAC fit Homography 的最小代码

核心概念

RANSAC (RANdom SAmple Consensus, 随机样本一致性) 是一种迭代算法，用于从包含大量离群点 (Outliers) 的数据集中鲁棒地估计数学模型的参数。其核心思想是，一个由真实内点 (Inliers) 构成的样本子集有很大概率可以拟合出正确的模型，而离群点则不会。因此，算法通过反复随机抽取一个最小样本集来拟合模型，然后用该模型检验所有数据点，最终选择那个被最多数据点（一致集）支持的模型作为最佳估计。

原理与推导

RANSAC 的目标是在有限的迭代次数内，以足够高的概率 p 找到一个“纯净”的、仅由内点组成的最小样本集。迭代次数 N 的公式正是为了计算达到此目标所需的尝试次数。

公式定义:

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

参数解释:

$N$ : 所需的迭代总次数。
$p$ : 至少有一次抽样得到“纯净”样本集的期望概率（通常设为 0.99 或 0.999）。这是算法的“成功率”或“置信度”。
$w$ : 数据集中内点 (Inlier) 的比例。例如，如果 100 个点中有 60 个是内点，则 $w = 0.6$ 。
$k$ : 拟合模型所需的最小样本点数量。例如，拟合一条直线需要 $k=2$ 个点，拟合一个单应性矩阵 (Homography) 需要 $k=4$ 对点。

推导过程:

单次抽样成功的概率: 假设我们从数据集中随机抽取一个点，该点是内点的概率为 $w$ 。为了拟合模型，我们需要抽取 $k$ 个点。由于是独立随机抽样，这 $k$ 个点全部都是内点的概率为：
$P(\text{一次抽样成功}) = w \times w \times \dots \times w = w^k$
“一次抽样成功”意味着我们得到了一个“纯净”的样本集，可以用来计算一个好的模型。
单次抽样失败的概率: 相应地，单次抽样失败（即抽到的 $k$ 个点中至少包含一个离群点）的概率为：
$P(\text{一次抽样失败}) = 1 - w^k$
连续 N 次抽样全部失败的概率: 如果我们进行 $N$ 次独立的迭代（抽样），那么这 $N$ 次抽样全部都失败的概率为：
$P(N \text{次全部失败}) = (1 - w^k)^N$
至少有一次成功的概率: 我们所期望的成功事件是“在 $N$ 次迭代中，至少有一次抽样成功”。这个事件是“ $N$ 次全部失败”的对立事件。因此，其概率为：
$p = P(\text{至少一次成功}) = 1 - P(N \text{次全部失败}) = 1 - (1 - w^k)^N$
求解迭代次数 N: 现在，我们从上面的等式中反解出 $N$ ：
$p = 1 - (1 - w^k)^N \\ 1 - p = (1 - w^k)^N \\ \log(1 - p) = \log((1 - w^k)^N) \\ \log(1 - p) = N \cdot \log(1 - w^k) \\ N = \frac{\log(1 - p)}{\log(1 - w^k)}$
推导完成。

直观解释: 这个公式告诉我们，内点比例 $w$ 越高，或者模型越简单（ $k$ 越小），我们找到正确模型所需的尝试次数 $N$ 就越少。反之，如果数据中离群点很多（ $w$ 很小），或者模型很复杂（ $k$ 很大），我们就需要急剧增加迭代次数 $N$ 才能保证以高概率 $p$ 找到一个好的解。

算法复杂度:

时间复杂度: $O(N \times (T_{\text{fit}} + M))$ ，其中 $N$ 是迭代次数， $T_{\text{fit}}$ 是使用 $k$ 个点拟合一次模型的时间， $M$ 是数据点的总数（用于计算一致集大小）。对于单应性矩阵， $T_{\text{fit}}$ 是一个小的常数（解一个小型线性系统），所以复杂度近似为 $O(N \cdot M)$ 。
空间复杂度: $O(M)$ ，主要用于存储原始数据点和内点索引。

代码实现

以下是使用 NumPy 手写 RANSAC 算法来拟合两组点之间的单应性矩阵 (Homography) 的最小实现。

python

 1import numpy as np
 2
 3def compute_homography(src_pts, dst_pts):
 4    """
 5    使用 DLT (Direct Linear Transform) 算法计算单应性矩阵 H。
 6    src_pts 和 dst_pts 均为 (4, 2) 的 NumPy 数组，代表 4 个匹配点对。
 7    """
 8    # 至少需要4对点来构建方程组
 9    assert src_pts.shape[0] >= 4 and dst_pts.shape[0] >= 4
10    
11    A = []
12    for i in range(src_pts.shape[0]):
13        x, y = src_pts[i]
14        u, v = dst_pts[i]
15        # 构建 Ax = 0 中的两行
16        A.append([-x, -y, -1, 0, 0, 0, u*x, u*y, u])
17        A.append([0, 0, 0, -x, -y, -1, v*x, v*y, v])
18    
19    A = np.asarray(A)
20    
21    # 使用 SVD 求解 Ah = 0
22    # h 是 A.T @ A 的最小特征值对应的特征向量
23    # 在 SVD 中，h 对应于 V.T 的最后一行 (或 V 的最后一列)
24    U, S, Vt = np.linalg.svd(A)
25    L = Vt[-1, :] / Vt[-1, -1] # 归一化，使得 h_22 = 1
26    H = L.reshape(3, 3)
27    
28    return H
29
30def fit_homography_ransac(src_pts, dst_pts, threshold, max_iterations=2000):
31    """
32    使用 RANSAC 算法鲁棒地拟合单应性矩阵 H。
33    
34    参数:
35    - src_pts, dst_pts: 匹配的点对，形状为 (N, 2)
36    - threshold: 内点判断阈值（重投影误差）
37    - max_iterations: 最大迭代次数
38    
39    返回:
40    - best_H: 最佳单应性矩阵 (3, 3)
41    - inliers_mask: 布尔数组，标记哪些点是内点
42    """
43    num_points = src_pts.shape[0]
44    k = 4  # 单应性矩阵需要4对点
45
46    best_H = None
47    max_inliers = 0
48    best_inliers_mask = None
49
50    # 将源点转换为齐次坐标 (N, 3)
51    src_pts_homo = np.hstack((src_pts, np.ones((num_points, 1))))
52
53    for i in range(max_iterations):
54        # 1. 随机采样 k 个点
55        # replace=False 确保我们选择的是 k 个不同的点
56        indices = np.random.choice(num_points, k, replace=False)
57        sample_src = src_pts[indices]
58        sample_dst = dst_pts[indices]
59
60        # 2. 拟合模型
61        # 使用采样的点计算一个候选的单应性矩阵
62        try:
63            H_candidate = compute_homography(sample_src, sample_dst)
64        except np.linalg.LinAlgError:
65            # 如果出现奇异矩阵等问题，跳过本次迭代
66            continue
67
68        # 3. 计算一致集
69        # 将所有源点通过候选 H 进行变换
70        transformed_dst_homo = (H_candidate @ src_pts_homo.T).T
71        
72        # 转换为非齐次坐标
73        # 防止除以零
74        transformed_dst_homo[:, 2][np.abs(transformed_dst_homo[:, 2]) < 1e-7] = 1e-7
75        transformed_dst = transformed_dst_homo[:, :2] / transformed_dst_homo[:, 2, np.newaxis]
76
77        # 计算重投影误差 (欧氏距离)
78        errors = np.linalg.norm(dst_pts - transformed_dst, axis=1)
79        
80        # 根据阈值确定内点
81        current_inliers_mask = errors < threshold
82        num_current_inliers = np.sum(current_inliers_mask)
83
84        # 4. 更新最佳模型
85        if num_current_inliers > max_inliers:
86            max_inliers = num_current_inliers
87            best_inliers_mask = current_inliers_mask
88            
89            # (可选但强烈推荐) 使用所有找到的内点重新拟合模型，以获得更精确的结果
90            inlier_src_pts = src_pts[best_inliers_mask]
91            inlier_dst_pts = dst_pts[best_inliers_mask]
92            best_H = compute_homography(inlier_src_pts, inlier_dst_pts)
93
94    print(f"RANSAC 完成: 找到 {max_inliers}/{num_points} 个内点。")
95    return best_H, best_inliers_mask
96
97if __name__ == '__main__':
98    # --- 1. 生成合成数据 ---
99    # 定义一个真实的 Homography
100    H_gt = np.array([
101        [0.9, -0.2, 50],
102        [0.1, 1.1, 20],
103        [0.0001, -0.0002, 1.0]
104    ])
105
106    # 生成内点
107    num_inliers = 80
108    inlier_src = np.random.rand(num_inliers, 2) * 100
109    inlier_src_homo = np.hstack((inlier_src, np.ones((num_inliers, 1))))
110    inlier_dst_homo = (H_gt @ inlier_src_homo.T).T
111    inlier_dst = inlier_dst_homo[:, :2] / inlier_dst_homo[:, 2, np.newaxis]
112    # 添加少量噪声
113    inlier_dst += np.random.randn(num_inliers, 2) * 0.5
114
115    # 生成外点
116    num_outliers = 40
117    outlier_src = np.random.rand(num_outliers, 2) * 100
118    outlier_dst = np.random.rand(num_outliers, 2) * 100 # 外点的目标点是完全随机的
119
120    # 混合数据
121    all_src_pts = np.vstack((inlier_src, outlier_src))
122    all_dst_pts = np.vstack((inlier_dst, outlier_dst))
123    
124    # 打乱顺序
125    shuffle_idx = np.random.permutation(len(all_src_pts))
126    all_src_pts = all_src_pts[shuffle_idx]
127    all_dst_pts = all_dst_pts[shuffle_idx]
128
129    # --- 2. 运行 RANSAC ---
130    # 设置参数
131    reprojection_threshold = 3.0  # 像素误差阈值
132    
133    # 根据理论公式估算迭代次数
134    # 假设我们预估内点比例 w = 0.5 (保守估计), k=4, p=0.99
135    w = 0.5
136    k = 4
137    p = 0.99
138    N_estimated = int(np.ceil(np.log(1 - p) / np.log(1 - w**k)))
139    print(f"理论估算迭代次数 N: {N_estimated}")
140
141    H_ransac, inliers = fit_homography_ransac(all_src_pts, all_dst_pts, reprojection_threshold, max_iterations=N_estimated)
142
143    # --- 3. 验证结果 ---
144    print("\n真实 Homography (GT):\n", H_gt / H_gt[2,2])
145    print("\nRANSAC 估计的 Homography:\n", H_ransac / H_ransac[2,2])
146
147    # 比较 OpenCV 的实现
148    try:
149        import cv2
150        H_cv, _ = cv2.findHomography(all_src_pts, all_dst_pts, cv2.RANSAC, reprojection_threshold)
151        print("\nOpenCV RANSAC 估计的 Homography:\n", H_cv / H_cv[2,2])
152    except ImportError:
153        print("\n未安装 OpenCV，跳过与 OpenCV 的比较。")

工程实践

自适应迭代次数 (Adaptive RANSAC): 在实际应用中，内点比例 $w$ 是未知的。一个常见的优化是动态调整迭代次数 $N$ 。在 RANSAC 循环中，每当找到一个更大的内点集时，就用当前的内点比例 $w_{\text{est}} = \frac{\text{当前内点数}}{\text{总点数}}$ 重新计算 $N_{\text{new}}$ 。如果已经执行的迭代次数超过了 $N_{\text{new}}$ ，就可以提前终止循环。这在内点比例很高时能显著节约计算时间。
阈值 threshold 的选择: 这是 RANSAC 最关键的超参数。它应该基于对数据噪声的理解来设定。对于图像特征点匹配，这个阈值通常是重投影误差，单位是像素。一个经验法则是设为 1 到 3 个像素。如果特征点定位精度高，可以设小一些；反之则设大一些。阈值太小会导致内点被误判为外点，太大则会导致外点被接受，污染模型。
最终模型的精炼: RANSAC 的核心是找到正确的内点集。一旦找到了最佳的内点集，应该使用所有这些内点来重新拟合模型（如代码中所示）。这会利用更多的数据，得到一个比仅用最小样本集拟合的模型更精确、更稳定的结果。
局部优化 (Local Optimization): 在找到一个包含较多内点的候选模型后，可以尝试对该模型进行局部优化。例如，以当前内点集为基础，用非线性最小二乘法（如 Levenberg-Marquardt）微调模型参数，可能会找到一个能容纳更多内点的新模型。USAC (Universal RANSAC) 框架集成了这类思想。
采样策略:
- PROSAC (Progressive Sample Consensus): 如果匹配点对有质量得分（如 SIFT 的匹配分数），可以优先从高质量的匹配点对中采样，而不是纯随机。这可以更快地找到一个好的样本集。
- 退化配置 (Degenerate Configurations): 采样时需要注意避免退化配置。例如，在拟合单应性矩阵时，如果采样的 4 个点中有 3 个或 4 个点共线，就无法唯一确定一个单应性矩阵。好的实现应该包含对这种退化情况的检测和处理。

常见误区与边界情况

误区：RANSAC 返回的模型是用最小样本集计算的。 这是一个常见的初学者错误。RANSAC 的最佳实践是，在迭代结束后，使用找到的最大一致集（所有内点）来重新计算并返回最终模型。
误区：迭代次数 N 是一个固定值。 虽然可以设置一个固定的、足够大的 max_iterations，但更高效的做法是如上所述的自适应调整。将 N 视为一个动态上限而不是固定循环次数。
边界情况：内点比例极低。 当真实内点比例 $w$ 非常低时（例如低于 10%），从公式可知所需的迭代次数 $N$ 会变得非常大，RANSAC 可能会因为达不到最大迭代次数而失败，或者随机找到一个由离群点组成的“伪一致集”。
边界情况：所有点都是内点。 如果数据非常干净，没有离群点，RANSAC 仍然可以工作，但它会做很多不必要的随机采样。在这种情况下，一个简单的最小二乘拟合会更高效。
面试追问:
- 问：如何确定阈值 threshold？ 答：它与特征点的定位精度和图像噪声有关。通常通过经验或实验确定。可以假设重投影误差服从卡方分布，根据置信区间来设定阈值。例如，假设误差的每个维度是标准正态分布，则误差的平方和服从自由度为2的卡方分布，可以取 95% 置信水平对应的值作为阈值。
- 问：RANSAC 和霍夫变换 (Hough Transform) 有什么区别？ 答：两者都用于从含噪声数据中检测模型。主要区别在于：RANSAC 在数据空间采样，然后验证模型；霍夫变换在参数空间进行投票。RANSAC 对高维参数空间更有效，而霍夫变换在参数维度增加时，内存和计算量会急剧增长。
- 问：除了 RANSAC，还有哪些鲁棒估计算法？ 答：LMedS (Least Median of Squares)，它最小化误差的中位数而不是计算内点数，对离群点更鲁棒但计算成本高。还有 M-estimators, PROSAC, USAC, MAGSAC 等改进版本。