您的位置:首页 > 其它

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;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息