2016-icpc-beijing

 2020.07.07 23:44  2023.09.01 18:14  ACM  计算几何 󰈭 2108字

二分答案然后跑圆的k次交即可。

因为学习完k次交后精疲力竭就不写完整这道题了。

贴一下用自己的码风写的板子。已更新入模版中。

给出n个圆的平面坐标和半径，能够获得所有圆k次交的面积。

 1const int manx = 1e2 + 10;
 2const double eps = 1e-8;
 3const double pi = acos(-1);
 4
 5inline int sgn(double x){return x < -eps? -1: x > eps? 1: 0;}
 6inline double sqr(double x){return x * x;}
 7
 8struct CIRCLE{
 9    double x, y, r, angle;
10    int d;
11    CIRCLE(){d = 1; angle = 0;}
12    CIRCLE(double _x, double _y, double _angle = 0, int _d = 0){
13        x = _x; y = _y; angle = _angle; d = _d; r = 0;
14    }
15    bool operator < (const CIRCLE &b)const{
16        return sgn(r - b.r) == -1;
17    }
18};
19
20inline double dis(CIRCLE a, CIRCLE b){return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));}
21inline double cross(CIRCLE a, CIRCLE b, CIRCLE c){
22    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
23}
24
25//有两个交点返回true，以及两个交点的坐标和方位角，p1按照顺时针，p2按照逆时针
26bool cir_inter_cir(CIRCLE a, CIRCLE b, CIRCLE &p1, CIRCLE &p2){
27    double d = dis(a, b);
28    if(sgn(d - a.r - b.r) >= 0 || sgn(abs(b.r - a.r) - d) >= 0) return false;
29    double cosa = (sqr(a.r) + sqr(d) - sqr(b.r)) / (2 * a.r * d);
30    double sina = sqrt(max(0., 1. - sqr(cosa)));
31    //旋转矩阵 [cosa, -sina; sina, cosa]
32    p1 = p2 = a;
33    p1.x += a.r / d * ((b.x - a.x) * cosa + (b.y - a.y) * -sina);
34    p1.y += a.r / d * ((b.x - a.x) * sina + (b.y - a.y) * cosa);
35    p2.x += a.r / d * ((b.x - a.x) * cosa + (b.y - a.y) * sina);
36    p2.y += a.r / d * ((b.x - a.x) * -sina + (b.y - a.y) * cosa);
37    return true;
38}
39
40bool cmp(const CIRCLE &a, const CIRCLE &b){
41    if(sgn(a.angle - b.angle)) return sgn(a.angle - b.angle) == -1;
42    else return a.d > b.d;
43}
44
45double cal(CIRCLE a, CIRCLE p1, CIRCLE p2){
46    double ans = sqr(a.r) * (p2.angle - p1.angle) - cross(a, p1, p2) + cross(CIRCLE(0, 0), p1, p2);
47    return ans / 2;
48}
49
50CIRCLE dot[maxn << 1];
51double area[maxn];
52void cirunion(CIRCLE cir[], int n){
53    sort(cir, cir + n);
54    //记录每个圆被完全覆盖的次数，初始值显然为1
55    for(int i = 0; i < n; i++) for(int j = i + 1; j < n; j++)
56        if(sgn(dis(cir[i], cir[j]) + cir[i].r - cir[j].r) <= 0) cir[i].d++;
57    CIRCLE p1, p2;
58    for(int i = 0; i < n; i++){
59        //针对每一个圆考虑它与所有其他圆的交
60        int cnt = 0, top = 0;
61        for(int j = 0; j < n; j++){
62            if(i == j) continue;
63            if(!cir_inter_cir(cir[i], cir[j], p1, p2)) continue;
64            p1.angle = atan2(p1.y - cir[i].y, p1.x - cir[i].x);
65            p2.angle = atan2(p2.y - cir[i].y, p2.x - cir[i].x);
66            p1.d = -1; p2.d = 1;
67            dot[top++] = p1; dot[top++] = p2;
68            if(sgn(p2.angle - p1.angle) == 1) cnt++;
69        }
70        
71        //加入起点终点位置
72        dot[top++] = CIRCLE(cir[i].x - cir[i].r, cir[i].y, -pi, -2);
73        dot[top++] = CIRCLE(cir[i].x - cir[i].r, cir[i].y, pi, 2);
74        sort(dot, dot + top, cmp);
75        int now = cir[i].d + cnt;
76        for(int j = 1; j < top; j++){
77            area[now] += cal(cir[i], dot[j - 1], dot[j]);
78            now += dot[j].d;
79        }
80    }
81}
82
83CIRCLE cir[maxn];
84int main(){
85//    Fastin;
86    int t; scanf("%d", &t); while(t--){
87        clr(area, 0);
88        int n; scanf("%d", &n);
89        for(int i = 0; i < n; i++) scanf("%lf %lf %lf", &cir[i].x, &cir[i].y, &cir[i].r);
90        cirunion(cir, n);
91        for(int i = 1; i <= n; i++) printf("%.5f\n", area[i]);
92    }
93    return 0;
94}

~~因为网上基本没有讲解只有代码..oi-wiki甚至好像都没有收录这个算法，所以我学了好久才…~~

下面是简单的个人理解

在主算法cirunion中，先将圆按照半径从小到大排序，接下来顺序遍历这些圆，$O(n^2)$地找出每个圆被另外的圆完全包含的次数，为后面的计算作准备。然后枚举每一个圆，考虑这个圆参与的重叠面积的求解：遍历其他所有圆，记录所有的交点，同时计算出这些交点对于圆心的角度(angle)，给每一个交点一个标记(d)用来之后处理重叠的次数。最后对这些交点排序，按照逆时针顺序遍历交点，更新area。

总体思路如上，但是具体细节我扣了很久。

在计算交点的时候，通过余弦定理获得角度，然后根据旋转的矩阵变换得到两个方向（逆时针和顺时针）的交点。
atan的值域只有(-pi/2, pi/2)，所以采用了atan2(int y, int x)，记录复数的幅角，范围为(-pi, pi]。
我们按照幅角angle排序，那么遍历的时候就是进行逆时针的遍历，从而如果遇到了通过逆时针操作得到的交点应该置d为1，表示此后的这段区域多一个圆重叠，而遇到了顺时针操作得到的交点应该置d为-1，表示此后这个弧就不再相交了，重叠的圆少一个。
在最后的更新的时候，维护当前重叠的圆的个数。我们从-pi的位置开始逆时针转，因为不知道最开始有多少个弧已经经过了-pi这个位置，所以之前在更新的时候要记录cnt值。而cnt值是这样产生的：因为在上述算法中p1是顺时针构造出的交点，p2是逆时针构造出的交点，所以如果不存在(-pi, pi]的限制p1的角度总是大于p2的，但是因为有了这个限制，所以如果这两个角度跨越了-pi，就会出现p2的角度大于p1，因而就可以根据p2的角度是否大于p1来判断这一条弧是否跨越了-pi。但是cnt值只能表示跨越了-pi的弧个数，实际在-pi处不一定有一个交点，言下之意就是起点不一定从-pi开始，那么cnt值难道根据第一个点的角度动态变化？实际上并不需要，我们添加两个冗余交点，坐标均为(x - r, y)，半径无关紧要，角度分别为-pi, pi, d值分别为2, -2。这样就可以保证在cmp比较符下这两个点一定是这个序列的第一个数和最后一个数。那么从这个-pi点开始的话，就可以通过cnt轻而易举的得知相交的那些弧的个数，再加上完全覆盖当前圆的个数，就得到了初值。之后遇到一个交点就根据它的标记更新当前重叠的个数即可。
最后一点是，如何根据遍历相邻的点对来获得重叠面积。这一个问题我想了非常久，因为之前没有计算几何的基础，所以百思不得其解。观察cal函数，前两项sqr(a.r) * (p2.angle - p1.angle) - cross(a, p1, p2)得到了弓形面积（钝角时因为cross是有向面积所以最终结果也是弓形），但是第三项…后来忽然意识到，多边形的面积好像就是通过不断地对一个固定点求有向面积最终求和即可，那么这边也就受到了启发：如果某个区域有若干的圆重叠得到，那么可以预想将那些相邻的交点相连，除去所有的弓形，剩下的就是一个多边形的面积，而这个面积就可以通过对一个定点求有向面积得到。但是这里并不是一次性求出来所有的邻边对应的有向面积，因为该算法遍历所有的圆，并且上述的每一个邻边上的两个点一定属于?