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

  |  

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

【对算法,数学,计算机感兴趣的同学,欢迎关注我哈,阅读更多原创文章】
我的网站:潮汐朝夕的生活实验室
我的公众号:算法题刷刷
我的知乎:潮汐朝夕
我的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