数几何

射线旋转扫描线算法(Angular-Sweep)

字数统计: 1.7k字   |   阅读时长: 7分
,

摘要: 扫描线算法：射线旋转扫描

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

各位好，扫描线算法是计算几何中非常重要的一个算法，并且可以用于解决很多区间类的问题。扫描线算法中的扫描方式分为直线平移扫描和射线旋转扫描两种。在文章 扫描线算法(Line-Sweep) 中我们详细介绍了直线平移算法的原理、代码模板，以及求 n 个矩形的并集的面积的代码模板。

本文我们以力扣第 1453 题来学习射线旋转扫描(Angular Sweep)算法。

题目

Alice 向一面非常大的墙上掷出 n 支飞镖。给你一个数组 darts ，其中 darts[i] = [xi, yi] 表示 Alice 掷出的第 i 支飞镖落在墙上的位置。

Bob 知道墙上所有 n 支飞镖的位置。他想要往墙上放置一个半径为 r 的圆形靶。使 Alice 掷出的飞镖尽可能多地落在靶上。

给你整数 r ，请返回能够落在 任意 半径为 r 的圆形靶内或靶上的最大飞镖数。

提示：

1
2
3
4
5
1 <= darts.length <= 100
darts[i].length == 2
-1e4 <= xi, yi <= 1e4
darts 中的元素互不相同
1 <= r <= 5000

示例 1 ：
输入：darts = [[-2,0],[2,0],[0,2],[0,-2]], r = 2
输出：4
解释：如果圆形靶的圆心为 (0,0) ，半径为 2 ，所有的飞镖都落在靶上，此时落在靶上的飞镖数最大，值为 4 。

示例 2 ：
输入：darts = [[-3,0],[3,0],[2,6],[5,4],[0,9],[7,8]], r = 5
输出：5
解释：如果圆形靶的圆心为 (0,4) ，半径为 5 ，则除了 (7,8) 之外的飞镖都落在靶上，此时落在靶上的飞镖数最大，值为 5 。

题解

算法1: 枚举圆心

枚举点对，特判需要跳过的情况，然后以两个点 P1, P2 和半径求出在 $\vec{P1P2}$ 左侧的圆心：

  • 特判1: 重合
  • 特判2: 距离大于半径的二倍

然后枚举除了选中的两个点外的所有点，计算与圆心的距离，判定是否在圆上或圆内，并计数。

每个点对计算完成后用所得的计数更新答案

代码 (C++)

$O(N^{3})$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
// 返回点 P 以点 O 为圆心逆时针旋转 alpha 后所在的位置
Vector2 rotate(Vector2 O, double alpha, Vector2 P)
{
    Vector2 new_P;
    Vector2 OP = P - O;
    new_P.x = OP.x * cos(alpha) - OP.y * sin(alpha) + O.x;
    new_P.y = OP.y * cos(alpha) + OP.x * sin(alpha) + O.y;
    return new_P;
}

Vector2 get_center(Vector2 P1, Vector2 P2, double r)
{
    Vector2 P1P2 = (P2 - P1);
    Vector2 P_mid = P1 + P1P2 * 0.5;
    double d = sqrt((r * r) - (P1P2.norm() * P1P2.norm() * 0.25));
    Vector2 O = P_mid + rotate(Vector2(0, 0), PI * 0.5, P1P2.normalize()) * d;
    return O;
}

class Solution {
public:
    int numPoints(vector<vector<int>>& points, int r) {
        vector<Vector2> p;
        int n = points.size();
        for(int i = 0; i < n; ++i)
            p.emplace_back(points[i][0], points[i][1]);
        int ans = 1;
        for(int i = 0; i < n; ++i)
        {
            Vector2 P1 = p[i];
            for(int j = 0; j < n; ++j)
            {
                if(i == j) continue;
                Vector2 P2 = p[j];
                if((P2 - P1).norm() < EPS)
                    continue;
                if((P2 - P1).norm() > r * 2.0)
                    continue;
                Vector2 O = get_center(P1, P2, r);
                int cnt = 0;
                for(int k = 0; k < n; ++k)
                {
                    if(k == i || k == j)
                        continue;
                    Vector2 P = p[k];
                    if((P - O).norm() - r < EPS)
                        ++cnt;
                }
                ans = max(ans, cnt + 2);
            }
        }
        return ans;
    }
};

算法2: Angular Sweep

从点集中选出一个点 P，做一个经过点 P 的，半径为 R 的圆周，然后将圆周围绕点 P 旋转。

例如图中，初始圆周的圆心为 C1，逆时针旋转 $\theta$ 之后到达 C2。

AngularSweepTheta

现在的问题转化为：是如何在旋转圆周的时候，得到圆周中的点的个数，记为 cnt。在圆周旋转的过程中，cnt 仅在某个点进入圆周或者出圆周的时机发生变化。

仅此可以枚举点集中的另一个点 Q，当 Q 进入圆周时，cnt 加 1，当 Q 离开圆周时，cnt 减 1。

圆周旋转的过程中，记旋转的角度为 $\theta$，圆心为 C，现在问题转化为，给定另一点 Q，Q 进入圆周时的 $\theta$，记为 $\alpha$，以及离开圆周的 $\theta$，记为 $\beta$。

具体求法：

对于点 Q，记 A 为 PQ 与 x 轴的夹角，B 为 PC 和 PQ

其中 d 为 PQ 的距离。

有了 A 和 B 之后，$\alpha$ 和 $\beta$ 如下

Q 进入圆周

Q 离开圆周

射线旋转扫描算法

1
2
3
4
5
6
预处理所有点对的距离 d[i][j]
对点 P，枚举所有 Q，若 PQ 距离大于 2R 则跳过，否则计算以下信息
    对 Q 计算的信息: 先计算 A, B, 进而计算 alpha, beta，想当与 Q 对应一个区间 [alpha, beta]
    计算 Q 的信息后，将区间的左右端点(即随着射线的旋转点 Q 的进入和离开)作为事件
    对所有事件排序
    遍历事件维护信息

事件的结构

1
2
3
4
5
6
7
8
9
10
11
struct Event
{
    double theta; // 角度 [0, 2PI)
    int ori; // 0 为进入，1 为离开
    bool operator<(const Event& e) const
    {
        if(fabs(theta - e.theta) < EPS)
            return ori < e.ori;
        return theta < e.theta;
    }
};

代码 (C++)

$O(N^{2}\log N)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
// 扫描线事件
struct Event
{
    double theta; // 角度 [0, 2PI)
    int ori; // 0 为进入，1 为离开
    Event(double theta, int ori):theta(theta),ori(ori){}
    bool operator<(const Event& e) const
    {
        if(fabs(theta - e.theta) < EPS)
            return ori < e.ori;
        return theta < e.theta;
    }
};

class Solution {
public:
    int numPoints(vector<vector<int>>& points, int r) {
        int n = points.size();
        vector<vector<double>> d(n, vector<double>(n));
        vector<Vector2> p(n);
        for(int i = 0; i < n; ++i)
        {
            p[i].x = points[i][0];
            p[i].y = points[i][1];
        }
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < n; ++j)
                d[i][j] = (p[i] - p[j]).norm();
        int ans = 0;
        for(int i = 0; i < n; ++i)
        {
            const Vector2& P = p[i];
            vector<Event> events;
            for(int j = 0; j < n; ++j)
            {
                if(i == j) continue;
                if(d[i][j] > r * 2)
                    continue;
                const Vector2& Q = p[j];
                // double A = (Q - P).polar();
                double A = atan((P.y - Q.y) / (P.x - Q.x));
                double B = acos(d[i][j] / (r * 2));
                double alpha = A - B;
                double beta = A + B;
                events.emplace_back(alpha, 0);
                events.emplace_back(beta, 1);
            }
            sort(events.begin(), events.end());
            int m = events.size();
            int cnt = 0;
            int max_cnt = 0;
            for(int j = 0; j < m; ++j)
            {
                if(events[j].ori == 0)
                    ++cnt;
                else
                    --cnt;
                max_cnt = max(max_cnt, cnt);
            }
            ans = max(ans, max_cnt + 1);
        }
        return ans;
    }
};
Share