【计算几何】圆的面积并

总之，圆的面积并是一个非常漂亮的算法，虽然说好像没有很强的扩展性以及实用性，但确实训练自己计算几何代码能力的好题目。

一下皆蒯自栗师《圆的并》解题报告:

    试想，如果人来做此题，而不使用计算机计算，那么，人会采取什么样的方法呢？

首先要做一些预处理，如果一个圆完全被另一个圆包围，那么这一个圆可以删除。删除后，如果一个圆是孤立的圆，不与其它任何圆相交，就可以把这个圆的面积现算出来，也将它删去。剩下的工作就是求交点了。

对于一个圆来说，其他的某些圆覆盖了他的圆弧上的某一段，即若干个区间。那么所有的圆覆盖的部分也是若干个区间。如果有k个区间被覆盖了，也就会有k个区间没有被覆盖。这2k个区间被2k个点分开。

 

我们用若干有向弦把这些部分分开，有向是指：沿圆的逆时针方向，

 

对所有的圆都这样处理了以后，只看这些连线，有什么发现？不难发现任意一个分割点一定会有两条连线！一条连线连进来，一条连出去。并且也只会有两条连线。这个试着画就可以了，例如如果3个圆经过通过同一个点，那么其中一定有一个圆，这个点的两边的弧都被覆盖了，即这个点不是这个圆的分割点！

如果任意一个分割点都会有两条连线，那么，这些连线之间形成了若干多边形。这些多边形就是我们要求的多边形的面积。但是，这些多边形中有的面积是负的，这个只需要看这个多边形的连线是顺时针的还是逆时针的，顺的为正，逆的为负。

 

 

代码相当的短：

program ex3;type real=extended;const e=1e-15;var cr,cx,cy,l,r:array[0..300] of real; cover:array[0..300] of boolean; x1,x2,y1,y2,t1,dx,dy,ans,tmp,dist:real; n,i,j,k,t,tt,ss:longint;procedure swap(var a,b:real);begin tmp:=a;a:=b;b:=tmp end;function eq(a,b:real):boolean;begin eq:=abs(a-b)<e end;function dis(a,b,c,d:real):real;begin dis:=sqrt(sqr(a-c)+sqr(b-d)) end;function ji(x,y:real):real;begin if x<0 then ji:=arctan(y/x)+pi else if y<0 then ji:=arctan(y/x)+pi*2 else ji:=arctan(y/x);end;procedure sort(ll,rr:longint);var i,j:longint;begin i:=ll;j:=rr;t1:=r[(i+j)>>1]; repeat while r[i]<t1-e do inc(i); while r[j]>t1+e do dec(j); if i<=j then begin swap(l[i],l[j]);swap(r[i],r[j]); inc(i);dec(j); end; until i>j; if ll<j then sort(ll,j); if i<rr then sort(i,rr);end;begin assign(input,'circle.in');reset(input); assign(output,'circle.out');rewrite(output); readln(n); for i:=1 to n do readln(cx[i],cy[i],cr[i]); for i:=1 to n do begin t:=0; for j:=1 to n do if (i<>j)and not cover[j] then begin dist:=dis(cx[i],cy[i],cx[j],cy[j]); if dist+cr[i]<cr[j]+e then cover[i]:=true else if dist+cr[j]<cr[i]+e then continue else if dist<cr[i]+cr[j]-e then begin inc(t); r[t]:=ji(cx[j]-cx[i]+e,cy[j]-cy[i]); dx:=(sqr(cr[i])-sqr(cr[j])+sqr(dist))/2/dist; dy:=sqrt(sqr(cr[i])-dx*dx); l[t]:=r[t]+ji(dx+e,dy); r[t]:=r[t]-ji(dx+e,dy); if r[t]<0 then begin r[t]:=r[t]+2*pi; l[t]:=l[t]+2*pi; end; end; end; if (t=0)and not cover[i] then ans:=ans+sqr(cr[i])*pi; if (t=0)or cover[i] then continue; sort(1,t);tt:=0;ss:=1; for j:=1 to t do if (tt=0)or(r[j]>l[tt]+e) then begin inc(tt);l[tt]:=l[j];r[tt]:=r[j]; end else if l[tt]<l[j] then l[tt]:=l[j]; while (ss<=tt)and(l[tt]-2*pi>r[ss]-e) do begin if l[tt]-2*pi<l[ss] then l[tt]:=l[ss]+2*pi; inc(ss); end; for j:=ss to tt do begin if j<>tt then k:=j+1 else k:=ss; dx:=r[k]-l[j]; if dx<0 then dx:=dx+2*pi; ans:=ans+sqr(cr[i])/2*(dx-sin(dx)); x1:=cx[i]+cr[i]*cos(l[j]); y1:=cy[i]+cr[i]*sin(l[j]); x2:=cx[i]+cr[i]*cos(r[k]); y2:=cy[i]+cr[i]*sin(r[k]); ans:=ans+(x1*y2-x2*y1)/2; end; end; writeln(ans:0:6); close(input);close(output);end.