三角刨分和面积计算

1. POJ2986 A Triangle and a Circle

题意：给定一个三角形，一个圆的圆心和半径，求圆和三角形的面积交

利用三角剖分，计算简单多边形和圆的相交面积

三角剖分的步骤：

1. 多边形上的每条边都与圆心构成三角形
2. 算出每个三角形与圆的相交面积
3. 根据有向面积的正负累加到答案中

计算每个三角形与圆的相交面积，分为5种情况:

1. 线段与圆心共线：返回0
2. 线段完全在圆内：1个三角形的有向面积
3. 线段完全在圆外：1个扇形的有向面积
4. 线段一端在圆内：1个三角形+1个扇形的有向面积
5. 线段是圆的割线：1个三角形+2个扇形的有向面积

注意：

1. 圆心位于坐标的原点
2. da,db 是圆心到线段两端点的距离
3. d 是圆心到线段的最近距离
4. pa,pb 是线段与圆的两个交点

时间：O(n*T)

 1 #include <iostream>
 2 #include <cstring>
 3 #include <algorithm>
 4 #include <cmath>
 5 #define x first
 6 #define y second
 7 using namespace std;
 8 
 9 typedef pair<double,double> Point;
10 const double eps=1e-8;
11 const double PI=acos(-1.0);
12 double R;
13 Point p[4],o; //顶点和圆心
14 
15 Point operator+(Point a,Point b){ //向量+
16   return Point(a.x+b.x,a.y+b.y);
17 }
18 Point operator-(Point a,Point b){ //向量-
19   return Point(a.x-b.x,a.y-b.y);
20 }
21 Point operator*(Point a,double t){ //数乘
22   return Point(a.x*t,a.y*t);
23 }
24 Point operator/(Point a,double t){ //数除
25   return Point(a.x/t,a.y/t);
26 }
27 double operator*(Point a,Point b){ //叉积
28   return a.x*b.y-a.y*b.x;
29 }
30 double operator&(Point a,Point b){ //点积
31   return a.x*b.x+a.y*b.y;
32 }
33 double len(Point a){ //模长
34   return sqrt(a&a);
35 }
36 double dis(Point a,Point b){ //距离
37   return len(b-a);
38 }
39 Point getNode(Point a,Point u,Point b,Point v){ //直线交点
40   double t=(a-b)*v/(v*u);
41   return a+u*t;
42 }
43 Point rotate(Point a,double b){ //逆转角
44   return Point(a.x*cos(b)-a.y*sin(b),a.x*sin(b)+a.y*cos(b));
45 }
46 bool onSegment(Point p,Point a,Point b){ //p在线段ab上
47   return fabs((a-p)*(b-p))<eps && ((a-p)&(b-p))<=0;
48 }
49 Point norm(Point a){ //单位向量
50   return a/len(a);
51 }
52 double getDP2(Point a,Point b,Point& pa,Point& pb){
53   Point e=getNode(a,b-a,o,rotate(b-a,PI/2)); //垂足
54   double d=dis(o,e);
55   if(!onSegment(e,a,b)) d=min(dis(o,a),dis(o,b));
56   if(R<=d) return d; //线段在圆外
57   double len=sqrt(R*R-dis(o,e)*dis(o,e));
58   pa=e+norm(a-b)*len;
59   pb=e+norm(b-a)*len; //pa,pb:线段与圆的两交点
60   return d;           //d:圆心到线段的最近距离
61 }
62 double sector(Point a,Point b){ //扇形面积
63   double angle=acos((a&b)/len(a)/len(b)); //[0,Pi]
64   if(a*b<0) angle=-angle;
65   return R*R*angle/2;
66 }
67 double getArea(Point a,Point b){ //面积的交
68   if(fabs(a*b)<eps) return 0; //共线
69   double da=dis(o,a),db=dis(o,b);
70   if(R>=da && R>=db) return a*b/2; //ab在圆内
71   Point pa,pb;
72   double d=getDP2(a,b,pa,pb); //d:圆心到线段的最近距离
73   if(R<=d) return sector(a,b); //ab在圆外
74   if(R>=da) return a*pb/2+sector(pb,b); //a在圆内
75   if(R>=db) return sector(a,pa)+pa*b/2; //b在圆内
76   return sector(a,pa)+pa*pb/2+sector(pb,b); //ab是割线
77 }
78 int main(){
79   while(scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf",
80   &p[0].x,&p[0].y,&p[1].x,&p[1].y,&p[2].x,&p[2].y,&o.x,&o.y,&R)!=-1){
81     for(int i=0;i<3;i++) p[i].x-=o.x,p[i].y-=o.y; //三角形顶点平移
82     o=Point(0,0); //圆心平移到原点
83     double res=0;
84     for(int i=0;i<3;i++) res+=getArea(p[i],p[(i+1)%3]);
85     printf("%.2lf\n",fabs(res)); //点可能顺时针
86   }
87   return 0;
88 }