【模板】多圆面积并

链接：https://vjudge.net/problem/SPOJ-CIRU

辛普森积分：

 

 

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 const int N = 1000+9;
 4 #define eps 1e-6
 5 #define inf 0x3f3f3f3f
 6 struct Cir{
 7     double x,y,r;
 8 }c[N],tmp[N];
 9 pair<double,double> seg[N];
10 bool cmp2(Cir a,Cir b){return a.r > b.r;}
11 bool cmp(Cir a,Cir b){
12     return (fabs(a.x - a.r - b.x + b.r ) < eps) ? a.x + a.r < b.x + b.r : a.x - a.r < b.x - b.r;
13 }
14 bool in_cir(Cir a,Cir b){
15     return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) <= (a.r -b.r) * (a.r - b.r);
16 }
17 int n,st,ed; 
18 void prework(){
19     sort(c+1,c+1+n,cmp2);
20     int cnt = 0;
21     for(int i = 1;i<=n;++i){
22         bool ok = 1;
23         for(int j = 1;j<=cnt;++j){
24             if(in_cir(c[i],tmp[j])){
25                 ok = 0;
26                 break;
27             }
28         }
29         if(ok) tmp[++cnt] = c[i];
30     }
31     n = cnt;
32     memcpy(c,tmp,sizeof(tmp));
33 }
34 double getf(double x){
35     int tot = 0;
36     for(int i = st;i<=ed;++i){
37         if(x < c[i].x + c[i].r && x > c[i].x - c[i].r){
38             double len = sqrt(c[i].r * c[i].r - ( x -c[i].x ) * (x - c[i].x) );
39             seg[++tot] = make_pair(c[i].y-len,c[i].y + len);
40         }
41     }
42     sort(seg+1,seg+1+tot);
43     double ans = 0,segl = -inf ,segr = -inf;
44     for(int i = 1;i<=tot;++i){
45         if(seg[i].first >= segr){
46             ans += segr - segl;
47             segl = seg[i].first;
48             segr = seg[i].second;
49         }
50         else segr = max(segr,seg[i].second);
51     }
52     ans += segr - segl;
53     return ans;
54 }
55 double calc(double len,double fL,double fM,double fR){
56     return (fL + 4*fM + fR) * len / 6;
57 }
58 double Simpson(double L,double M,double R,double fL,double fM,double fR,double sum){
59     double M1 = (L+M)/2 , M2 = (M+R)/2;
60     double fM1 = getf(M1), fM2 = getf(M2);
61     double g1 = calc(M-L,fL,fM1,fM) , g2 = calc(R-M,fM,fM2,fR);
62     if(fabs(sum - g1 - g2) <= eps ) return g1 + g2;
63     return Simpson(L,M1,M,fL,fM1,fM,g1) + Simpson(M,M2,R,fM,fM2,fR,g2);
64 }
65 void solve(){
66     sort(c+1,c+1+n,cmp);
67     double ans = 0;
68     for(int i = 1;i<=n;++i){
69         double L = c[i].x - c[i].r , R = c[i].x + c[i].r;
70         int j;
71         for(j = i+1;j<=n;++j){
72             if(c[j].x - c[j].r > R) break;
73             else R = max(R,c[j].x + c[j].r);
74         }
75         double M = (L+R)/2;
76         st = i,ed = j-1;
77         i = j-1;
78         double fL = getf(L),fM = getf(M),fR = getf(R);
79         ans += Simpson(L,M,R,fL,fM,fR,calc(R-L,fL,fM,fR));
80     }
81     printf("%.3f\n",ans);
82 }
83 int main(){
84     scanf("%d",&n);
85     for(int i = 1;i<=n;++i) scanf("%lf %lf %lf",&c[i].x,&c[i].y,&c[i].r);
86     prework();
87     solve();
88 }
View Code

 

格林公式：

参考：https://www.cnblogs.com/jefflyy/p/9247686.html

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 const int N = 1000+9;
 4 const double PI = acos(-1.0);
 5 struct Point{
 6     double x,y;
 7     bool operator < (const Point& b)const{
 8         return x<b.x || (x==b.x && y<b.y);
 9     }
10     bool operator == (const Point& b)const{
11         return x==b.x && y == b.y;
12     }
13     Point operator - (const Point& b)const{
14         return (Point){x - b.x,y - b.y};
15     }
16     double len(){return sqrt(x*x + y*y);}
17 };
18 struct Cir{
19     Point o;
20     double r;
21     bool operator < (const Cir& b)const{
22         return o<b.o || (o==b.o && r<b.r);
23     }
24     bool operator == (const Cir& b)const{
25         return o == b.o && r == b.r;
26     }
27     double oint(double t1,double t2){
28         return r*(r*(t2 - t1) + o.x * (sin(t2) - sin(t1)) - o.y * (cos(t2) - cos(t1)));
29     }
30 }c[N],tmp[N];
31 pair<double,int> st[N<<1];
32 int n;
33 double calc(int id){
34     int top = 0, cnt = 0;
35     for(int i = 1;i<=n;++i){
36         if(i != id){
37             Point d = c[i].o - c[id].o;
38             double dis = d.len();
39             if(c[id].r + dis <= c[i].r) return 0;
40             if(c[i].r + dis <= c[id].r) continue;
41             if(dis >= c[i].r + c[id].r) continue;
42             double ang = atan2(d.y,d.x);
43             double deta = acos( (c[id].r * c[id].r + dis * dis - c[i].r * c[i].r) / (2 * c[id].r * dis));
44             double l = ang - deta, r = ang + deta;
45             if( l < -PI ) l += 2*PI;
46             if( r >= PI ) r -= 2*PI;
47             if( l > r ) ++cnt; 
48             st[++top] = make_pair(l,1);
49             st[++top] = make_pair(r,-1);
50         }
51     }
52     st[0] = make_pair(-PI,0);
53     st[++top] = make_pair(PI,0);
54     sort(st+1,st+1+top);
55     double res = 0;
56     for(int i = 1;i<=top;cnt += st[i++].second){
57         if(!cnt) res += c[id].oint(st[i-1].first,st[i].first);
58 
59     }
60     return res;
61 }
62 int main(){
63     scanf("%d",&n);
64     for(int i = 1;i<=n;++i) scanf("%lf %lf %lf",&c[i].o.x,&c[i].o.y,&c[i].r);
65     sort(c+1,c+1+n); n = unique(c+1,c+1+n) - c - 1;
66     double ans = 0;
67     for(int i = 1;i<=n;++i) ans += calc(i);
68     ans *= 0.5;
69     printf("%.3lf",ans);
70     return 0;
71 }
View Code

 

posted @ 2019-11-30 14:23  小布鞋  阅读(154)  评论(0)    收藏  举报