计算几何--简单多边形与圆面积交
求解二维空间内一个简单多边形和一个长度为R的圆公共面积。
因为任意简单多边形都可以划分成若干三角形，我们可以把这个简单多边形划分成三角形后，求三角形与圆的面积交，然后在把所有三角形的解合并。
由于可能有凹多边形，我们计算三角形与圆面积交时采用向量叉乘，这样得到的是一个有向面积，刚好可以把凹多边形面积正负抵消掉，最后把总面积取绝对值就行了。
向量叉乘 A x B == 以向量A,B为2邻边，围城平行四边形的有向面积。 A在B顺时针方向值为正，逆时针为负。
AxB==
|A.x , A.y |
|B.x , B.y |
==A.x*B.y-A.y*B.x
计算一个圆与一个三角形的面积交(其中一个三角形顶点是圆心，如上图)，我采用的方法是分4种情况。
1.
另外2个顶点在圆内（上），这个非常好算直接求三角形的有向面积即可。
2.
另外两个顶点有1个再圆内（上）,另外1个再圆外，求得直线与圆一个交点后求一个三角形面积+上一个扇形面积。
3.
2个顶点在圆外，且2个顶点所在边与圆相交，先求圆外2顶点所在直线与圆交点，然后定比分点公式求另外2条直线与圆交点，然后求一个三角形+2个扇形面积即可。
4.
2个顶点都在圆外且2顶点所在边与圆不相交，这个情况求2个交点后算出那个扇形面积就行了。
下面是我写的圆与三角形有向面积交函数，注意三角形其中一个顶点在圆心，如果都不在圆心，可以把这个三角形在划分成3个其中一个顶点在圆心的三角形求解。
1 /************************************** 2 Author : lxgsbqylbk 3 Date : 2012/08/12 4 Function : Direct area of a circle and triangle 5 ***********/ 6 const double eps = 1e-8; //浮点数精度控制 7 8 struct point //点或者向量结构 9 { 10 double x,y; 11 point(double _x=0.0,double _y=0.0) 12 : x(_x),y(_y) {} 13 point operator - (const point & v) 14 { 15 return point(x-v.x,y-v.y); 16 } 17 double sqrx() //向量的模 18 { 19 return sqrt(x*x+y*y); 20 } 21 }; 22 double xmult(point & p1,point & p2,point & p0) //叉乘 23 { 24 return (p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x); 25 } 26 double distancex(point & p1,point & p2) 27 { 28 return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); 29 } 30 point intersection(point u1,point u2,point v1,point v2) //两直线交点 31 { 32 point ret=u1; 33 double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x)) 34 /((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x)); 35 ret.x+=(u2.x-u1.x)*t; 36 ret.y+=(u2.y-u1.y)*t; 37 return ret; 38 } 39 void intersection_line_circle(point c,double r,point l1,point l2,point& p1,point& p2){ 40 point p=c; 41 double t; 42 p.x+=l1.y-l2.y; 43 p.y+=l2.x-l1.x; 44 p=intersection(p,c,l1,l2); 45 t=sqrt(r*r-distancex(p,c)*distancex(p,c))/distancex(l1,l2); 46 p1.x=p.x+(l2.x-l1.x)*t; 47 p1.y=p.y+(l2.y-l1.y)*t; 48 p2.x=p.x-(l2.x-l1.x)*t; 49 p2.y=p.y-(l2.y-l1.y)*t; 50 } 51 point ptoseg(point p,point l1,point l2) //点到线段的最近距离 52 { 53 point t=p; 54 t.x+=l1.y-l2.y,t.y+=l2.x-l1.x; 55 if (xmult(l1,t,p)*xmult(l2,t,p)>eps) 56 return distancex(p,l1)<distancex(p,l2)?l1:l2; 57 return intersection(p,t,l1,l2); 58 } 59 double distp(point & a,point & b) 60 { 61 return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y); 62 } 63 double Direct_Triangle_Circle_Area(point a,point b,point o,double r) 64 { 65 double sign=1.0; 66 a=a-o; 67 b=b-o; 68 o=point(0.0,0.0); 69 if(fabs(xmult(a,b,o))<eps) return 0.0; 70 if(distp(a,o)>distp(b,o)) 71 { 72 swap(a,b); 73 sign=-1.0; 74 } 75 if(distp(a,o)<r*r+eps) 76 { 77 if(distp(b,o)<r*r+eps) return xmult(a,b,o)/2.0*sign; 78 point p1,p2; 79 intersection_line_circle(o,r,a,b,p1,p2); 80 if(distancex(p1,b)>distancex(p2,b)) swap(p1,p2); 81 double ret1=fabs(xmult(a,p1,o)); 82 double ret2=acos( p1*b/p1.sqrx()/b.sqrx() )*r*r; 83 double ret=(ret1+ret2)/2.0; 84 if(xmult(a,b,o)<eps && sign>0.0 || xmult(a,b,o)>eps && sign<0.0) ret=-ret; 85 return ret; 86 } 87 point ins=ptoseg(o,a,b); 88 if(distp(o,ins)>r*r-eps) 89 { 90 double ret=acos( a*b/a.sqrx()/b.sqrx() )*r*r/2.0; 91 if(xmult(a,b,o)<eps && sign>0.0 || xmult(a,b,o)>eps && sign<0.0) ret=-ret; 92 return ret; 93 } 94 point p1,p2; 95 intersection_line_circle(o,r,a,b,p1,p2); 96 double cm=r/(distancex(o,a)-r); 97 point m=point( (o.x+cm*a.x)/(1+cm) , (o.y+cm*a.y)/(1+cm) ); 98 double cn=r/(distancex(o,b)-r); 99 point n=point( (o.x+cn*b.x)/(1+cn) , (o.y+cn*b.y)/(1+cn) ); 100 double ret1 = acos( m*n/m.sqrx()/n.sqrx() )*r*r; 101 double ret2 = acos( p1*p2/p1.sqrx()/p2.sqrx() )*r*r-fabs(xmult(p1,p2,o)); 102 double ret=(ret1-ret2)/2.0; 103 if(xmult(a,b,o)<eps && sign>0.0 || xmult(a,b,o)>eps && sign<0.0) ret=-ret; 104 return ret; 105 }
