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 ", solve() + EPS); 173 }
以下转自：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 ", ans / 6.0); 102 } 103 104 int main(){ 105 while( get() ) work(); 106 return 0; 107 }
【圆并的扩展算法】
https://www.spoj.pl/problems/CIRUT/
【题目大意】
给出N个不同的圆，求出被覆盖恰好K次的面积，并顺序输出
【思路提示】
参考圆并和圆交的思路， 就是把问题转化为区间来考虑
同时请看下面的图：
上面的图，是对于被覆盖恰好2次的区域的边界连接起来得到的图形，我们可以发现什么呢？
(就不剧透了~~~~)