Sphere Online Judge (SPOJ) - Problem CIRU
【求圆并的若干种算法,圆并扩展算法】_AekdyCoin的空间_百度空间
参考AekdyCoin的圆并算法解释，根据理解写出的代码。圆并这么多题中，最基础一题。
操作如下：
（1）对一个圆，获得所有与其他圆的交点作为时间戳。
a.如果这个圆不被其他任何圆覆盖，这个圆直接保留。
b.如果Case中有两个完全重合的圆，可以保留标号小（或大）的圆，也只保留这一个。
（2）每经过一个时间戳就对计数器加减，如果计数器从0变1，加上这段圆弧和对应的三角形。
（3）所有这样的圆独立计算，最后求出的面积累加即可。
代码如下(1y)：
1 #include <cmath> 2 #include <cstdio> 3 #include <cstring> 4 #include <iomanip> 5 #include <iostream> 6 #include <algorithm> 7 8 using namespace std; 9 10 const double PI = acos(-1.0); 11 const double EPS = 1e-8; 12 inline int sgn(double x) { return (x > EPS) - (x < -EPS);} 13 template<class T> T sqr(T x) { return x * x;} 14 typedef pair<double, double> Point; 15 #define x first 16 #define y second 17 Point operator + (Point a, Point b) { return Point(a.x + b.x, a.y + b.y);} 18 Point operator - (Point a, Point b) { return Point(a.x - b.x, a.y - b.y);} 19 Point operator * (Point a, double p) { return Point(a.x * p, a.y * p);} 20 Point operator / (Point a, double p) { return Point(a.x / p, a.y / p);} 21 22 inline double cross(Point a, Point b) { return a.x * b.y - a.y * b.x;} 23 inline double dot(Point a, Point b) { return a.x * b.x + a.y * b.y;} 24 inline double angle(Point a) { return atan2(a.y, a.x);} 25 inline double veclen(Point a) { return sqrt(dot(a, a));} 26 27 struct Circle { 28 Point c; 29 double r; 30 Circle() {} 31 Circle(Point c, double r) : c(c), r(r) {} 32 Point point(double p) { return Point(c.x + r * cos(p), c.y + r * sin(p));} 33 void get() { scanf("%lf%lf%lf", &c.x, &c.y, &r);} 34 } ; 35 36 bool ccint(Circle a, Circle b, double &a1, double &a2) { // ips for Circle-a (anti-clockwise p1->p2) 37 double d = veclen(a.c - b.c); 38 double r1 = a.r, r2 = b.r; 39 if (sgn(d - r1 - r2) >= 0) return 0; 40 if (sgn(fabs(r1 - r2) - d) >= 0) return 0; 41 double da = acos((sqr(d) + sqr(r1) - sqr(r2)) / (2 * d * r1)); 42 double ang = angle(b.c - a.c); 43 a1 = ang - da, a2 = ang + da; 44 return 1; 45 } 46 47 const int N = 1111; 48 Circle cir[N]; 49 50 inline double gettri(Point o, Point a, Point b) { return cross(a - o, b - o) / 2;} 51 inline double getarc(Circle o, double a, double b) { return (b - a) * sqr(o.r) / 2 - gettri(o.c, o.point(a), o.point(b));} 52 53 typedef pair<double, int> Event; 54 Event ev[N << 2]; 55 56 bool inside(int a, int n) { 57 double d, r1, r2; 58 for (int i = 0; i < n; i++) { 59 if (i == a) continue; 60 d = veclen(cir[i].c - cir[a].c); 61 r1 = cir[i].r, r2 = cir[a].r; 62 if (sgn(r1 - r2) == 0 && sgn(fabs(r1 - r2) - d) == 0) { 63 if (i > a) return 1; 64 continue; 65 } 66 if (sgn(r1 - r2) >= 0 && sgn(fabs(r1 - r2) - d) >= 0) return 1; 67 } 68 return 0; 69 } 70 71 double getarea(int a, int n) { 72 if (inside(a, n)) return 0; 73 Circle &c = cir[a]; 74 double ret = 0; 75 double a1, a2; 76 Point p1, p2; 77 int tt = 0; 78 for (int i = 0; i < n; i++) { 79 if (i == a) continue; 80 if (ccint(c, cir[i], a1, a2)) { 81 if (a2 > PI) ev[tt++] = Event(a1, 1), ev[tt++] = Event(PI, -1), ev[tt++] = Event(-PI, 1), ev[tt++] = Event(a2 - 2 * PI, -1); 82 else if (a1 < -PI) ev[tt++] = Event(a1 + 2 * PI, 1), ev[tt++] = Event(PI, -1), ev[tt++] = Event(-PI, 1), ev[tt++] = Event(a2, -1); 83 else ev[tt++] = Event(a1, 1), ev[tt++] = Event(a2, -1); 84 } 85 } 86 if (tt == 0) return PI * sqr(c.r); 87 sort(ev, ev + tt); 88 int cnt = 0; 89 double last = -PI; 90 for (int i = 0; i < tt; i++) { 91 double cur = ev[i].x; 92 if (ev[i].y == 1 && cnt == 0) ret += getarc(c, last, cur) + cross(c.point(last), c.point(cur)) / 2; 93 cnt += ev[i].y; 94 last = cur; 95 } 96 ret += getarc(c, last, PI) + cross(c.point(last), c.point(PI)) / 2; 97 return ret; 98 } 99 100 double work(int n) { 101 double ret = 0; 102 for (int i = 0; i < n; i++) ret += getarea(i, n); 103 return ret; 104 } 105 106 int main() { 107 int n; 108 while (~scanf("%d", &n)) { 109 for (int i = 0; i < n; i++) cir[i].get(); 110 printf("%.3f ", work(n)); 111 } 112 return 0; 113 }
——written by Lyon