SPOJ 8073 The area of the union of circles（计算几何の圆并）（CIRU）

Description

You are given N circles and expected to calculate the area of the union of the circles !

Input

The first line is one integer n indicates the number of the circles. (1 <= n <= 1000)

Then follows n lines every line has three integers

Xi Yi Ri

indicates the coordinate of the center of the circle, and the radius. (|Xi|. |Yi| <= 1000, Ri <= 1000)

Note that in this problem Ri may be 0 and it just means one point !

Output

The total area that these N circles with 3 digits after decimal point

题目大意：求n个圆覆盖的总面积。

思路：参考http://hi.baidu.com/aekdycoin/item/b8ff6adc73c0e71dd78ed0d6

时间复杂度O(n^2*log(n))

代码（0.04S）：

  1 #include <cstdio>
  2 #include <algorithm>
  3 #include <cstring>
  4 #include <iostream>
  5 #include <cmath>
  6 #include <vector>
  7 using namespace std;
  8 typedef long long LL;
  9 
 10 const double PI = acos(-1.0);
 11 const double EPS = 1e-8;
 12 
 13 inline int sgn(double x) {
 14     return (x > EPS) - (x < -EPS);
 15 }
 16 
 17 struct Point {
 18     double x, y;
 19     Point() {}
 20     Point(double x, double y): x(x), y(y) {}
 21     void read() {
 22         scanf("%lf%lf", &x, &y);
 23     }
 24     double angle() {
 25         return atan2(y, x);
 26     }
 27     Point operator + (const Point &rhs) const {
 28         return Point(x + rhs.x, y + rhs.y);
 29     }
 30     Point operator - (const Point &rhs) const {
 31         return Point(x - rhs.x, y - rhs.y);
 32     }
 33     Point operator * (double t) const {
 34         return Point(x * t, y * t);
 35     }
 36     Point operator / (double t) const {
 37         return Point(x / t, y / t);
 38     }
 39     double length() const {
 40         return sqrt(x * x + y * y);
 41     }
 42     Point unit() const {
 43         double l = length();
 44         return Point(x / l, y / l);
 45     }
 46 };
 47 
 48 double cross(const Point &a, const Point &b) {
 49     return a.x * b.y - a.y * b.x;
 50 }
 51 
 52 double dist(const Point &p1, const Point &p2) {
 53     return (p1 - p2).length();
 54 }
 55 
 56 Point rotate(const Point &p, double angle, const Point &o = Point(0, 0)) {
 57     Point t = p - o;
 58     double x = t.x * cos(angle) - t.y * sin(angle);
 59     double y = t.y * cos(angle) + t.x * sin(angle);
 60     return Point(x, y) + o;
 61 }
 62 
 63 struct Region {
 64     double st, ed;
 65     Region() {}
 66     Region(double st, double ed): st(st), ed(ed) {}
 67     bool operator < (const Region &rhs) const {
 68         if(sgn(st - rhs.st)) return st < rhs.st;
 69         return ed < rhs.ed;
 70     }
 71 };
 72 
 73 struct Circle {
 74     Point c;
 75     double r;
 76     vector<Region> reg;
 77     Circle() {}
 78     Circle(Point c, double r): c(c), r(r) {}
 79     void read() {
 80         c.read();
 81         scanf("%lf", &r);
 82     }
 83     void add(const Region &r) {
 84         reg.push_back(r);
 85     }
 86     bool contain(const Circle &cir) const {
 87         return sgn(dist(cir.c, c) + cir.r - r) <= 0;
 88     }
 89     bool intersect(const Circle &cir) const {
 90         return sgn(dist(cir.c, c) - cir.r - r) < 0;
 91     }
 92 };
 93 
 94 double sqr(double x) {
 95     return x * x;
 96 }
 97 
 98 void intersection(const Circle &cir1, const Circle &cir2, Point &p1, Point &p2) {
 99     double l = dist(cir1.c, cir2.c);
100     double d = (sqr(l) - sqr(cir2.r) + sqr(cir1.r)) / (2 * l);
101     double d2 = sqrt(sqr(cir1.r) - sqr(d));
102     Point mid = cir1.c + (cir2.c - cir1.c).unit() * d;
103     Point v = rotate(cir2.c - cir1.c, PI / 2).unit() * d2;
104     p1 = mid + v, p2 = mid - v;
105 }
106 
107 Point calc(const Circle &cir, double angle) {
108     Point p = Point(cir.c.x + cir.r, cir.c.y);
109     return rotate(p, angle, cir.c);
110 }
111 
112 const int MAXN = 1010;
113 
114 Circle cir[MAXN];
115 bool del[MAXN];
116 int n;
117 
118 double solve() {
119     double ans = 0;
120     for(int i = 0; i < n; ++i) {
121         for(int j = 0; j < n; ++j) if(!del[j]) {
122             if(i == j) continue;
123             if(cir[j].contain(cir[i])) {
124                 del[i] = true;
125                 break;
126             }
127         }
128     }
129     for(int i = 0; i < n; ++i) if(!del[i]) {
130         Circle &mc = cir[i];
131         Point p1, p2;
132         bool flag = false;
133         for(int j = 0; j < n; ++j) if(!del[j]) {
134             if(i == j) continue;
135             if(!mc.intersect(cir[j])) continue;
136             flag = true;
137             intersection(mc, cir[j], p1, p2);
138             double rs = (p2 - mc.c).angle(), rt = (p1 - mc.c).angle();
139             if(sgn(rs) < 0) rs += 2 * PI;
140             if(sgn(rt) < 0) rt += 2 * PI;
141             if(sgn(rs - rt) > 0) mc.add(Region(rs, PI * 2)), mc.add(Region(0, rt));
142             else mc.add(Region(rs, rt));
143         }
144         if(!flag) {
145             ans += PI * sqr(mc.r);
146             continue;
147         }
148         sort(mc.reg.begin(), mc.reg.end());
149         int cnt = 1;
150         for(int j = 1; j < int(mc.reg.size()); ++j) {
151             if(sgn(mc.reg[cnt - 1].ed - mc.reg[j].st) >= 0) {
152                 mc.reg[cnt - 1].ed = max(mc.reg[cnt - 1].ed, mc.reg[j].ed);
153             } else mc.reg[cnt++] = mc.reg[j];
154         }
155         mc.add(Region());
156         mc.reg[cnt] = mc.reg[0];
157         for(int j = 0; j < cnt; ++j) {
158             p1 = calc(mc, mc.reg[j].ed);
159             p2 = calc(mc, mc.reg[j + 1].st);
160             ans += cross(p1, p2) / 2;
161             double angle = mc.reg[j + 1].st - mc.reg[j].ed;
162             if(sgn(angle) < 0) angle += 2 * PI;
163             ans += 0.5 * sqr(mc.r) * (angle - sin(angle));
164         }
165     }
166     return ans;
167 }
168 
169 int main() {
170     scanf("%d", &n);
171     for(int i = 0; i < n; ++i) cir[i].read();
172     printf("%.3f\n", solve() + EPS);
173 }

View Code

以下转自：http://hi.baidu.com/aekdycoin/item/b8ff6adc73c0e71dd78ed0d6

【求圆并的若干种算法,圆并扩展算法】

【问题求解】
给定N 个圆形，求出其并集.

【算法分析】

PS.以下算法基于正方向为逆时针

考虑上图中的蓝色圆,绿色的圆和蓝色的圆交于 A,B 2个交点 ,我们在逆时针系下考虑，那么可以知道对于蓝色的圆，它对应的某个 角度区间被覆盖了

假设区间为 [A, B], 并且角度是按照圆心到交点的向量的极角来定义 (为了方便，我一般都把角度转化到 [0,2pi]区间内) 那么可以知道在这种标识情况下，可能存在以下情况

这种区间的跨度如何解决呢?实际上十分简单，只需要把[A,B] 拆成 [A, 2PI], [0,B]即可,也就是所谓的添加虚拟点

下面介绍一下对于我们当前所求任务的实际运用( 利用上述做法)

首先对于所给的N个圆，我们可以进行去冗杂，表现为:
(1) 去掉被包含(内含 or 内切)的小圆 ()
(2) 去掉相同的圆

枚举一个圆，并对于剩下的圆和它求交点，对于所求的的交点，可以得到一个角度区间 [A,B], 当然区间如果跨越了(例如 [1.5PI, 0.5PI],注意这里是有方向的) 2PI那么需要拆区间

可以知道，最后区间的并集必然是最后所有圆和当前圆的交集的一个边界！

于是我们得到互补区间(所谓互补区间就是[0,2PI]内除去区间并的区间,可能有多个)

假设我们先枚举了橙色的圆，那么得到了许多角度区间，可以知道绿色的和蓝色的角度区间是“未被覆盖的”，对于未被覆盖的

圆弧染色！

而对于其他圆，我们也做同样的步骤，同时把他们未被覆盖的角度区间的圆弧标记为黑色阴影

于是最终的结果就如下图 (染色只考虑圆弧)

通过观察不难发现,圆的并是许多的圆弧+ 许多多边形的面积之和(注意这里为简单多边形,并且面积有正负之别！)

于是我们累加互补区间的圆弧面积到答案，并把互补区间确定的弦的有向面积累加到答案

(例如在上面的例子中，我们在扫描到橙色圆的这一步只需要累加蓝色和绿色的圆弧面积 以及 蓝色和绿色的有向面积，注意这里蓝色和绿色的边必然是最后那个多边形的边！)

这里涉及到一个问题，就是：
圆弧可能大于一半的圆，例如上图中最大的圆，当然如果你推出了公式，那么实际上很容易发现无论圆弧长啥样都能算出正确的答案！
半径为R的圆中，弧度区间为K的圆弧的面积为 : 0.5 * r * r * (K - sin(K))

之后还是有一些疑惑，如果存在"洞"?

美妙之处在于此,由于面积有方向(正负),所以洞会被自然抵消！(可以类比计算简单多边形面积)

枚举圆,并和其他圆求交，求得区间，排序
这里是O(n^2 logn)
至于区间扫描和累加的复杂度则是O(n)
于是复杂度
O(n * (nlogn + n) )

O(n^2 logn)

代码已实现,大约140行,和圆的交的思路十分类似，圆的交的思路请参考:

http://hi.baidu.com/aekdycoin/blog/item/12267a4e9476153bafc3abbd.html

圆并的题目：

http://www.spoj.pl/problems/CIRU/

https://www.spoj.pl/problems/VCIRCLES/

【圆并的数值积分】

// 先贴代码，稍后补上

  1 #include<iostream>
  2 #include<stdio.h>
  3 #include<string.h>
  4 #include<algorithm>
  5 #include<cmath>
  6 using namespace std;
  7 const int maxn = 1005;
  8 typedef double db;
  9 const db EPS = 1e-6;
 10 typedef pair<db, db> PDD;
 11 int x[ maxn ], y[ maxn ], r[ maxn ];
 12 int nx[ maxn ], ny[ maxn ], nr[ maxn ];
 13 int xl[ maxn ], xr[ maxn ];
 14 
 15 int s[ maxn ];
 16 inline bool cmp( int a, int b) {
 17      if( x[ a ] - r [ a ] == x[ b ] - r [ b ] ) return x[ a ] + r[ a ] < x[ b ] + r [ b ];
 18      return x[ a ] - r[ a ] < x[ b ] - r [ b ];
 19 }
 20 inline bool cmp0(int a, int b){return r[ a ] > r [ b ];}
 21 int n;
 22 int L, R;
 23 PDD se[ maxn ];
 24 inline db f( db v){
 25    int sz = 0, i, j ;
 26    db ret = 0.0;
 27    for(i = L; i < R; ++ i){
 28         if( v <= xl[ i ] || v >= xr[ i ] ) continue;
 29         j = s[ i ];
 30         db d = sqrt(r[ j ]- (v - x [ j ]) * (v - x[ j ]));
 31         se[ sz ].first = y[ j ] - d;
 32         se[ sz ].second = y[ j ] +  d;
 33         ++ sz;   
 34    }
 35    sort( se, se + sz);
 36    for(i = 0; i < sz; ++ i){
 37          db nowl , nowr;
 38          nowl = se[ i ].first;
 39          nowr = se[ i ].second;
 40          for( j = i + 1; j < sz; ++ j) if(se[ j ].first > nowr) break;
 41          else nowr = max( nowr, se[ j ].second);
 42          ret += nowr - nowl;
 43          i = j - 1;      
 44    }
 45    return ret;
 46 }
 47 #define fs(x) ((x) < 0 ? (-(x)) : (x))
 48 inline db rsimp( db l,db m, db r, db sl, db sm, db sr,db tot){
 49     db m1 = (l + m) * 0.5, m2 = (m + r) * 0.5;
 50     db s0 = f( m1), s2 = f( m2);
 51     db gl = (sl + sm + s0 + s0 + s0 + s0)*(m-l), gr = (sm + sr + s2 + s2 + s2 + s2)*(r-m);
 52     if( fs(gl + gr - tot) < EPS) return gl + gr;
 53     return rsimp( l, m1, m, sl, s0, sm, gl) + rsimp( m, m2,r, sm, s2, sr, gr);         
 54 }
 55 
 56 bool get(){ 
 57      if(1 != scanf("%d", &n)) return 0;
 58      int i, j = 0, k;
 59      for(i = 0; i < n; ++ i) scanf("%d%d%d", x + i, y + i, r + i), s[ i ] = i;
 60      sort( s, s + n, cmp0);
 61      for(i = 0; i < n; ++ i){
 62            for(k = 0; k < j; ++ k)
 63                  if( (nx [ k ] - x [s[i]]) * (nx [ k ] - x [s[i]])  + (ny [ k ] - y [s[i]]) *(ny [ k ] - y [s[i]])  
 64                      <= (nr[ k ] - r[ s[ i ] ]) * (nr[ k ] - r[ s[ i ] ]) ) break;
 65            if( k == j) {
 66                nx[ j ] = x[ s[ i ] ];
 67                ny[ j ] = y[ s[ i ] ];
 68                nr[ j ] = r[ s[ i ] ];
 69                s[ j ] = j;
 70                j ++;    
 71            }      
 72      }
 73      n = j;
 74      for(i = 0; i < n; ++ i) x[ i ] = nx[ i ], y[ i ] = ny[ i ], r[ i ] = nr[ i ];
 75      return 1;
 76 }
 77   
 78 void work(){
 79      sort( s, s + n, cmp) ;
 80      db lo, hi, ans = 0.0;
 81      int i, j;
 82      for(i = 0; i < n; ++ i) xl[ i ] = x[ s[ i ] ] - r[ s[ i ] ], xr[ i ] = x[ s[ i ] ] + r[ s[ i ] ], r[ s[i] ] *= r[ s[i] ];
 83      for(i = 0; i < n; ++ i){
 84            int ilo, ihi;
 85            ilo = xl[ i ];
 86            ihi = xr[ i ];
 87            for(j = i + 1; j < n; ++ j) {
 88                  if( xl[ j ] > ihi ) break;
 89                  ihi = max( ihi, xr[ j ]);
 90            }
 91            db lo = ilo;
 92            db hi = ihi;      
 93            L = i;
 94            R = j;
 95            db mid = (lo + hi) * 0.5;
 96            db sl = f(lo), sm = f(mid), sr = f(hi);
 97            db tot = sl + sr + sm + sm + sm + sm;
 98            ans += rsimp( lo, mid , hi, sl, sm , sr, tot );
 99            i = j - 1;
100      }
101      printf("%.3f\n", ans / 6.0);
102 }
103   
104 int main(){
105     while( get() ) work();
106     return 0;
107 }