【多维分桶】力扣149:自定义哈希函数&RANSAC算法

字数统计: 2.8k字 | 阅读时长: 13分

2020-05-20

摘要: 力扣 149 题，分桶算法的经典题

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

在文章分桶法中，我们系统梳理了力扣上关于分桶法的题目。今天这个题目是多维分桶的经典题目

按照与固定点连成的线的斜率(long double)对点做哈希，用哈希值分桶
用多层 treemap 实现多属性(横坐标和纵坐标)分桶
STL无序关联容器自定义哈希函数

$1 题目

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

$2 题解

算法1 : 哈希分桶

枚举所有点，当枚举到点 i 时，枚举所有其它的点 j。

用 unordered_map 记录不同的点的个数，这里的相同指的是到点 i 的斜率相同。如果 j1 到 i 的斜率与 j2 到 i 的斜率相同，则 j1, j2 相同，在 unordered_map 中属于同一个键。

斜率的表示如果用 double 或者 long double 会有精度的问题，用 pair<int, int> 比较保险。对两个点求斜率，返回的是一个数对，它是一个向量。

返回的数对应该是没有公约数的，这里用到了辗转相除法。

对斜率做哈希即是对返回的数 pair<int, int> slope 对做哈希，hash(slope) := ((ll)hash<int>()(slope.first) * 10007 + hash<int>()(slope.second)) % MOD;。

自定义哈希函数相关的写法问题，参考 STL无序关联容器自定义哈希函数

每枚举完一个点 i，都可以求出一个出现最多的点，将这个最多的个数如果比当前答案大，就更新答案。

对于后枚举的点 i，如果其余的点 j 中出现最多的点里有在 i 之前的，例如点 k ($k < i$)，则该答案在枚举到 k 时候就更新到答案里了，这里也就不用枚举点 k 了，因此枚举 j 时所有小于 i 的点就不用考虑了。

需要特判的边界条件：

重合的点
两点连线是竖线

代码(C++)

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;
        int ans = 1;
        for(int i = 0; i < n; ++i)
        {
            Point origin(points[i][0], points[i][1]);
            int result = 0;
            int duplicates = 1;
            unordered_map<PPP, int, MyHash, MyCmp> mapping;
            mapping.insert({{origin, origin}, 1});
            for(int j = i + 1; j < n; ++j)
            {
                Point cur(points[j][0], points[j][1]);
                if(cur == origin)
                    ++duplicates;
            }
            for(int j = i + 1; j < n; ++j)
            {
                Point cur(points[j][0], points[j][1]);
                if(cur == origin)
                    continue;
                PPP cur_pair({origin, cur});
                ++mapping[cur_pair];
                result = max(result, mapping[cur_pair]);
            }
            result += duplicates;
            ans = max(ans, result);
        }
        return ans;
    }

private:
    using PII = pair<int, int>;
    class Point
    {
    public:
        Point(int x, int y): x(x), y(y) {}

        PII get_slope(const Point& origin) const
        {
            int dx = x - origin.x, dy = y - origin.y;
            if(dx == 0 && dy == 0) return PII({0, 0});
            if(dx == 0) return PII({0, 1});
            if(dy == 0) return PII({1, 0});
            if(dy < 0)
            {
                dy = -dy;
                dx = -dx;
            }
            int dx_sign = 1;
            if(dx < 0)
            {
                dx = -dx;
                dx_sign = -1;
            }
            int gcd = _gcd(dx, dy);
            dx /= gcd;
            dx *= dx_sign;
            dy /= gcd;
            return PII({dx, dy});
        }

        int get_x() const {return x;}
        int get_y() const {return y;}

        bool operator ==(const Point& point) const
        {
            return point.get_x() == x && point.get_y() == y;
        }

        bool operator !=(const Point& point) const
        {
            return !(*this == point);
        }

    private:
        int x, y;
        int _gcd(int x, int y) const
        {
            return x == 0 ? y : _gcd(y % x, x);
        }
    };

    using PPP = pair<Point, Point>;

    class MyHash
    {
    public:
        int operator()(const PPP& point_pair) const
        {
            PII slope = point_pair.second.get_slope(point_pair.first);
            return ((ll)hash<int>()(slope.first) * 10007 + hash<int>()(slope.second)) % MOD;
        }
        const int MOD = 1e9 + 7;
        using ll = long long;
    };

    class MyCmp
    {
    public:
        bool operator()(const PPP& point_pair1, const PPP& point_pair2) const
        {
            PII slope1 = point_pair1.second.get_slope(point_pair1.first);
            PII slope2 = point_pair2.second.get_slope(point_pair2.first);
            return slope1 == slope2;
        }
    };
};

算法2 : 多维 TreeMap 分桶

class Solution {
public:
    int maxPoints(vector<vector<int>>& points) {
        int n = points.size();
        if(n == 1)
            return 1;
        int max_cnt = 0;
        for(int i = 0; i < n - 1; ++i)
        {
            // 经过 points[i] 的直线
            // mapping[i][j] := 方向向量除以最大公约数后为 (i, j) 的点数, i >= 0
            map<int, map<int, int>> mapping;
            int duplicate = 1;
            for(int j = i + 1; j < n; ++j)
            {
                int a = points[i][0] - points[j][0];
                int b = points[i][1] - points[j][1];
                if(a == 0 && b == 0)
                {
                    ++duplicate;
                    continue;
                }
                if(a < 0)
                {
                    a = -a;
                    b = -b;
                }
                if(a != 0 && b != 0)
                {
                    int gcd_ = gcd<int>(abs(a), abs(b));
                    a /= gcd_;
                    b /= gcd_;
                }
                else if(b != 0)
                    b = INT_MAX / 2;
                else if(a != 0)
                    a = INT_MAX / 2;
                ++mapping[a][b];
            }
            max_cnt = max(max_cnt, duplicate);
            for(auto a: mapping)
                for(auto b: a.second)
                    max_cnt = max(max_cnt, b.second + duplicate);
        }
        return max_cnt;
    }
};

$3 面试题 16.14. 最佳直线

多层 TreeMap

class Solution {
public:
    vector<int> bestLine(vector<vector<int>>& points) {
        int n = points.size();
        int max_cnt = 0;
        int p1 = -1, p2 = -1;
        for(int i = 0; i < n - 1; ++i)
        {
            // 经过 points[i] 的直线
            // mapping[i][j] := 方向向量除以最大公约数后为 (i, j) 的点数, i >= 0
            map<int, map<int, set<int>>> mapping;
            vector<int> duplicate(1, i);
            for(int j = i + 1; j < n; ++j)
            {
                int a = points[i][0] - points[j][0];
                int b = points[i][1] - points[j][1];
                if(a == 0 && b == 0)
                {
                    duplicate.push_back(j);
                    continue;
                }
                if(a < 0)
                {
                    a = -a;
                    b = -b;
                }
                if(a != 0 && b != 0)
                {
                    int gcd_ = gcd<int>(abs(a), abs(b));
                    a /= gcd_;
                    b /= gcd_;
                }
                else if(b != 0)
                    b = INT_MAX / 2;
                else if(a != 0)
                    a = INT_MAX / 2;
                mapping[a][b].insert(j);
            }
            for(auto a: mapping)
                for(auto b: a.second)
                {
                    set<int> &setting = b.second;
                    for(int i: duplicate)
                        setting.insert(i);
                    if((int)setting.size() > max_cnt)
                    {
                        max_cnt = b.second.size();
                        p1 = *setting.begin();
                        p2 = *++(setting.begin());
                    }
                    else if((int)setting.size() == max_cnt)
                    {
                        if(*setting.begin() < p1)
                        {
                            p1 = *setting.begin();
                            p2 = *++(setting.begin());
                        }
                        else if(*setting.begin() == p1 && *++setting.begin() < p2)
                        {
                            p2 = *++(setting.begin());
                        }
                    }
                }
        }
        return {p1, p2};
    }
};

$4 扩展: RANSAC 算法

RANSAC简介

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

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

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

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

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

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

RANSAC 的基本思想和流程

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

具体实现步骤：

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

迭代次数推导

假设内点在数据中的占比为 $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$ 来估算出迭代次数。

代码(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);
}

取 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;
    }
};

潮汐朝夕的生活实验室 \Doge 陪伴一个算法工程师的职业生涯

算法leetcode