求解二維空間內一個簡單多邊形和一個長度為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 }