SPOJ - CIRU The area of the union of circles (圆形面积并)
2015-11-17 01:13
477 查看
题意:
给n个圆的圆心与半径,求这些圆的面积并。思路:
两个办法,一个是很优美的几何做法:http://www.cnblogs.com/oyking/p/3424999.html,
一个是用自适应Simpson积分来来求圆并:http://blog.csdn.net/qpswwww/article/details/44201333 。
代码:
//几何方法 #include<bits/stdc++.h> using namespace std; const int inf = 0x3f3f3f3f; const double eps = 1e-8 ; const double pi = acos(-1.0); inline double sqr(double x) { return x * x ; } inline int sgn(double x) { if (fabs(x) < eps) return 0 ; return x > 0? 1 : -1 ; } struct Point{ double x,y ; Point(){} Point(double _x, double _y): x(_x), y(_y){} double norm() { return sqrt(sqr(x) + sqr(y)) ; } void input() { scanf("%lf%lf", &x, &y) ; } friend Point operator + (Point a, Point b) { return Point(a.x + b.x, a.y + b.y) ; } friend Point operator - (Point a, Point b) { return Point(a.x - b.x, a.y - b.y) ; } friend Point operator * (Point a, double b) { return Point(a.x * b, a.y * b) ; } friend Point operator / (Point a, double b) { return Point(a.x / b, a.y / b) ; } friend bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0 ; } friend bool operator < (Point a, Point b) { return sgn(a.x - b.x) < 0 || (sgn(a.x - b.x) == 0 && sgn(a.y - b.y) < 0) ; } Point rotate(Point p, double ang) { Point v = (*this) - p ; Point ret ; ret.x = v.x * cos(ang) - v.y * sin(ang) ; ret.y = v.y * cos(ang) + v.x * sin(ang) ; return ret + p ; } }; double det(Point a, Point b) { return a.x * b.y - b.x * a.y; } int CircleInterCircle(Point c1, double r1, Point c2, double r2, Point *res) { double d = (c1 - c2).norm() ; if (sgn(d) == 0) { if (sgn(r1 - r2) == 0) return -1 ; return 0 ; } if (sgn(r1 + r2 - d) < 0) return 0 ; if (sgn(fabs(r1 - r2) - d) > 0) return 0 ; double ang = atan2((c2 - c1).y, (c2 - c1).x) ; double vang = acos( (sqr(r1) + sqr(d) - sqr(r2)) / (2*r1*d) ) ; res[0] = Point(c1.x+r1, c1.y).rotate(c1, ang + vang) ; res[1] = Point(c1.x+r1, c1.y).rotate(c1, ang - vang) ; if (res[0] == res[1]) return 1 ; return 2 ; } struct Region { double st, ed; Region(){} Region(double _st, double _ed): st(_st), ed(_ed) {} bool operator < (const Region &a)const { return sgn(st - a.st) < 0 || (sgn(st - a.st) == 0 && sgn(ed - a.ed) < 0) ; } }; struct Circle { Point c ; double r ; vector<Region> reg ; Circle(){} Circle(Point _c, double _r): c(_c), r(_r) {} void add(const Region &r) { reg.push_back(r); } double area(double ang = pi) { return ang * sqr(r); } Point makepoint(double ang) { return Point(c.x + r*cos(ang) , c.y + r*sin(ang)); } bool operator < (const Circle &a)const { return sgn(r - a.r) < 0 || (sgn(r - a.r) == 0 && c < a.c) ; } bool operator == (const Circle &a)const { return sgn(r - a.r) == 0 && c == a.c ; } }; double Angle(Point p) { return atan2(p.y, p.x); } double area_CircleUnion(Circle *cir, int n) { bool ok[n+5] ; memset(ok, true, sizeof ok) ; double ans = 0 ; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) if (ok[j]) { if (i == j) continue ; double d = (cir[i].c - cir[j].c).norm() ; if (sgn(d + cir[i].r - cir[j].r) <= 0) { //!!!! ok[i] = false ; break ; } } } for (int i = 0; i < n; i++) if (ok[i]) { Point p[2] ; bool flag = false ; for (int j = 0; j < n; j++) if(ok[j]) { if (i == j) continue ; int k = CircleInterCircle(cir[i].c, cir[i].r, cir[j].c, cir[j].r, p) ; if (k != 2) continue ; flag = true ; double ang1 = Angle(p[1] - cir[i].c), ang2 = Angle(p[0] - cir[i].c) ; if (sgn(ang1) < 0) ang1 += 2*pi; if (sgn(ang2) < 0) ang2 += 2*pi; if (sgn(ang1 - ang2) > 0) cir[i].add(Region(ang1, 2*pi)), cir[i].add(Region(0, ang2)) ; else cir[i].add(Region(ang1, ang2)) ; } if (!flag) { ans += cir[i].area() ; continue ; } sort(cir[i].reg.begin(), cir[i].reg.end()) ; int cnt = 1 ; for (int j = 1; j < int(cir[i].reg.size()); j++) { if (sgn(cir[i].reg[cnt - 1].ed - cir[i].reg[j].st) >= 0) { cir[i].reg[cnt - 1].ed = max(cir[i].reg[cnt - 1].ed, cir[i].reg[j].ed) ; } else { cir[i].reg[cnt++] = cir[i].reg[j] ; } } cir[i].add(Region()) ; cir[i].reg[cnt] = cir[i].reg[0] ; for (int j = 0; j < cnt; j++) { p[0] = cir[i].makepoint(cir[i].reg[j].ed) ; p[1] = cir[i].makepoint(cir[i].reg[j+1].st) ; ans += det(p[0], p[1]) / 2.0 ; double ang = cir[i].reg[j + 1].st - cir[i].reg[j].ed ; if (sgn(ang) < 0) ang += 2*pi ; ans += 0.5 * sqr(cir[i].r) * (ang - sin(ang)) ; } } return ans ; } int n ; Circle cir[1010] ; int main() { while (cin >> n) { for (int i = 0; i < n; i++) { cir[i].c.input() ; scanf("%lf", &cir[i].r) ; } printf("%.3f\n", area_CircleUnion(cir, n) + eps) ; } return 0; }
相关文章推荐
- 初学算法 - 求凸包的Garham's Scan算法的C++实现
- 【Google Code Jam 2009 round2 problem D】Watering Plants (两圆交点求法详解)
- 计算几何模板
- 计算几何小模板
- BZOJ2829信用卡凸包
- HDU 4922 Hello, Your Package! (计算几何+DP)(WA)
- poj 1514&zoj 1185 Metal Cutting(半平面交)
- UVA 10969 Sweet Dream(圆的相交)
- uva 11177 Fighting Against a Polygonal Monster(凸包与圆的面积交)
- POJ1279 && LA2512 Art Gallery(求多边形的核)
- poj 2540 && uva 10084 Hotter Colder(半平面交)
- 【计算几何】POJ 2318 & POJ 2398
- 【计算几何】POJ 2653
- 【计算几何】POJ 1113
- HDU 5128 The E-pang Palace
- POJ 2318 TOYS(叉积+二分or暴力)
- POJ 2398 Toy Storage(叉积+二分)
- POJ 1228 Grandpa's Estate 计算凸包+判断点在线段上
- POJ 1873 The Fortified Forest 计算凸包
- POJ 2007 Scrambled Polygon 极角排序