RANSAC算法：应对大比例异常值

摘要: RANSAC 算法思想

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

---

RANSAC 算法是一种随机算法，它用迭代的方式从一组含有离群点的一组观测数据中估计数学模型的参数。本文首先简要介绍这种算法的思路，然后用该算法解决之前用多维分桶解决的 N 个点中最多的共线点数问题，参考文章 【多维分桶】N 个点中最多有多少点共线。

---

RANSAC简介

RANSAC(RAndom SAmple Consensus,随机采样一致) 算法是从一组含有外点(outliers)的数据中正确估计数学模型参数的迭代算法。

外点的含义是数据中的噪声。因此 RANSAC 也称外点检测算法。

RANSAC算法是一种概率算法，它只能在某个概率下产生正确结果，并且这个概率会随着迭代次数的增加而加大。

RANSAC算法来说一个基本的假设: 数据是由内点和外点组成的。

• 内点: 组成模型参数的数据
• 外点: 不适合模型的数据

RANSAC 的另一个假设：在给定一组含有少部分内点的数据，存在一个程序可以估计出符合内点的模型。

RANSAC 的基本思想和流程

RANSAC是通过反复选择数据集去估计出模型，一直迭代到估计出认为比较好的模型。

具体实现步骤

• step1. 选择出可以估计出模型的最小数据集;(对于直线拟合来说就是两个点，对于计算 Homography 矩阵就是 4 个点)
• step2. 使用这个数据集来计算出数据模型；
• step3. 将所有数据带入这个模型，计算出内点的数目；(累加在一定误差范围内的适合当前迭代推出模型的数据)
• step4. 比较当前模型和之前推出的最好的模型的内点的数量，记录最大内点数的模型参数和内点数；
• step5. 重复 step1 - step4，直到迭代结束或者当前模型已经足够好了(内点数目大于一定数量)。

后面我们通过用 RANSAC 算法解决 N 个点中最多的共线点数问题，看一下以上步骤具体是怎么展开的。

迭代次数推导

假设内点在数据中的占比为 $t$。

$$t = \frac{n_{inliers}}{n_{inliers} + n_{outliers}}$$

每次计算模型使用 $N$ 个点的情况下，选取的点至少有一个外点的情况就是 $1 - t^{N}$，也就是在迭代 k 次的情况下，$(1 - t^{N})^{k}$ 就是 k 次迭代模型都至少采样到一个外点去计算模型的概率。

那么能采样到正确的 N 个点取计算出正确模型的概率就是 $P = 1 - (1 - t^{n})^{k}$。

$$k = \frac{\log(1 - P)}{\log(1 - t^{n})}$$

内点的概率 $t$ 通常是一个先验值。然后 $P$ 是我们希望RANSAC得到正确模型的概率。

如果事先不知道 $t$ 的值，可以使用自适应迭代次数的方法。也就是一开始设定一个无穷大的迭代次数，然后每次更新模型参数估计的时候，用当前的内点比值当成 $t$ 来估算出迭代次数。

---

149. 直线上最多的点数

给你一个数组 points ，其中 points[i] = [xi, yi] 表示 X-Y 平面上的一个点。求最多有多少个点在同一条直线上。

提示：

1 <= points.length <= 300
points[i].length == 2
-1e4 <= xi, yi <= 1e4
points 中的所有点 互不相同

示例 1:
输入: [[1,1],[2,2],[3,3]]
输出: 3
解释:

^
|
|        o
|     o
|  o  
+------------->
0  1  2  3  4

示例 2:
输入: [[1,1],[3,2],[5,3],[4,1],[2,3],[1,4]]
输出: 4
解释:

^
|
|  o
|     o        o
|        o
|  o        o
+------------------->
0  1  2  3  4  5  6

算法：RANSAC

按照 RANSAC 算法的思想，处理最多的共线点数问题的算法如下。

首先设置数据和模型之间可接受的差值 eps = 1e-10，设置正确概率 P=0.9999。

迭代之前设置迭代最大次数 iter_num = INT_MAX，每次得到更好的估计会优化 iter_num 的数值。

迭代过程维护的数学模型参数

设置迭代过程维护的数学模型的参数，即直线 y = ax + b 的 a 和 b，初始值均为 0。

int best_a = 0;
int best_b = 0;

迭代过程

在每一轮迭代过程中，按以下步骤处理：

• step1：随机选出两个点
• step2：用选出的两个点求解模型 y = ax + b
• step3：算出内点数目，即与 step2 估计出的 y = ax + b 共线的点数
• step4：判断当前模型是否比之前估的模型好，同时涉及到迭代步数的更新。

在 STL 容器中随机选择

如果用 C++ 实现，涉及到在 STL 容器中随机选择，仿照 <algorithm> 中的写法如下：

template<typename Iter, typename RandomGenerator>
Iter select_randomly(Iter start, Iter end, RandomGenerator *g) {
    std::uniform_int_distribution<> dis(0, std::distance(start, end) - 1);
    std::advance(start, dis(*g));
    return start;
}

template<typename Iter>
Iter select_randomly(Iter start, Iter end) {
    static std::random_device rd;
    static std::mt19937 gen(rd());
    return select_randomly(start, end, &gen);
}

代码(c++)

取 eps = 1e-10, P = 0.9999 时跑的结果比哈希表的算法快很多。

template<typename Iter, typename RandomGenerator>
Iter select_randomly(Iter start, Iter end, RandomGenerator *g) {
    std::uniform_int_distribution<> dis(0, std::distance(start, end) - 1);
    std::advance(start, dis(*g));
    return start;
}

template<typename Iter>
Iter select_randomly(Iter start, Iter end) {
    static std::random_device rd;
    static std::mt19937 gen(rd());
    return select_randomly(start, end, &gen);
}

class Solution {
public:
    int maxPoints(vector<vector<int>>& points) {
        if(points.empty()) return 0;
        int n = points.size();
        if(n == 1 || n == 2) return n;
        // 迭代最大次数，每次得到更好的估计会优化 iter_num 的数值
        int iter_num = INT_MAX;
        // 数据和模型之间可接受的差值
        const long double eps = 1e-10;
        // 最好模型的参数估计和内点数目
        int best_a = 0;
        int best_b = 0;
        int ans = 0;
        // 希望的得到正确模型的概率
        double P = 0.9999;

        for(int i = 0; i < iter_num; ++i)
        {
            // 随机选两个点求解模型
            auto it1 = select_randomly(points.begin(), points.end());
            auto it2 = select_randomly(points.begin(), points.end());
            while(it1 - points.begin() == it2 - points.begin())
                it2 = select_randomly(points.begin(), points.end());
            int x1 = (*it1)[0];
            int y1 = (*it1)[1];
            int x2 = (*it2)[0];
            int y2 = (*it2)[1];

            // y = ax + b 求解出a，b
            long double a, b;
            if(x2 == x1) // 两个点相同，或者垂直
            {
                a = (long double)INT_MAX;
                b = (long double)INT_MAX;
            }
            else if(y2 == y1) // 水平
            {
                a = 0.0;
                b = y2;
            }
            else
            {
                a = (long double)(y2 - y1) / (x2 - x1);
                b = y1 - a * x1;
            }

            // 算出内点数目
            int total_inlier = 0;
            for(int j = 0; j < n; ++j)
            {
                int x = points[j][0], y = points[j][1];
                long double y_estimate = 0.0;
                if(abs(a - (long double)INT_MAX) < eps)
                {
                    if(x == x1)
                        ++total_inlier;
                }
                else if(abs(a - 0.0) < eps)
                {
                    if(y == y1)
                        ++total_inlier;
                }
                else
                {
                    y_estimate = a * x + b;
                    if(abs(y_estimate - y) < eps)
                        ++total_inlier;
                }
            }

            // 判断当前的模型是否比之前估算的模型好
            if(total_inlier > ans)
            {
                best_a = a;
                best_b = b;
                ans = total_inlier;
                iter_num = log(1 - P) / log(1 - pow((double)total_inlier / n, 2));
            }

            if(total_inlier >= n) break;
        }
        return ans;
    }
};

---