poj 3675 多边形与圆的面积交(转) 

  1 /*
  2 代码：(转)http://blog.csdn.net/u010527492/article/details/18190147
  3 题解：(转)http://hi.baidu.com/billdu/item/703ad4e15d819db52f140b0b
  4 */
  5 #include<cstdio>
  6 #include<algorithm>
  7 #include<cmath>
  8 
  9 const double eps = 1e-8;
 10 
 11 using namespace std;
 12 
 13 struct point{
 14     double x,y;
 15     point(){}
 16     point(double x_,double y_):x(x_),y(y_){}
 17 }p[58];
 18 
 19 double r; // 所求圆半径
 20 point origin; // 所求圆圆心
 21 
 22 int n;
 23 double MIN(double x,double y){ return x<y?x:y; }
 24 double MAX(double x,double y){ return x>y?x:y; }
 25 
 26 double cross(point p1,point p2,point p3) {
 27     return (p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);
 28 }
 29 double dot(point p1,point p2) {
 30     return p1.x*p2.x+p1.y*p2.y;
 31 }
 32 double getdis(point p1,point p2) {
 33     return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
 34 }
 35 
 36 point get_intersect(point a, point b) { // 线段与圆的交点
 37     point temp=point(a.x-b.x,a.y-b.y);
 38     point vec=point(temp.y,-temp.x);
 39     point origin2=point(origin.x+vec.x,origin.y+vec.y);
 40     double a1=a.y-b.y;
 41     double b1=b.x-a.x;
 42     double c1=(a.x*b.y-b.x*a.y);
 43 
 44     double a2=origin.y-origin2.y;
 45     double b2=origin2.x-origin.x;
 46     double c2=(origin.x*origin2.y-origin2.x*origin.y);
 47 
 48     double cc=a1*b2-a2*b1;
 49     return point((b1*c2-b2*c1)/cc,(a2*c1-a1*c2)/cc);
 50 };
 51 
 52 int on_line(point p0,point p1,point p2) { // 点p0是否在线段p1-p2上
 53     if(p0.x>MAX(p1.x,p2.x)) return 0;
 54     if(p0.x<MIN(p1.x,p2.x)) return 0;
 55     if(p0.y>MAX(p1.y,p2.y)) return 0;
 56     if(p0.y<MIN(p1.y,p2.y)) return 0;
 57     return 1;
 58 }
 59 double get_area(point a, point b) { // 获取某点为圆心的三角形与该圆的相交面积
 60     double disa = getdis(origin,a);
 61     double disb = getdis(origin,b);
 62 
 63     double angle = acos(dot(a,b)/disa/disb); // 求两向量夹角
 64     point inter = get_intersect(a,b);
 65     double disinter = getdis(inter,origin);
 66     double sgn = cross(origin,a,b) < 0 ? -1 : 1;
 67     double ret = 0;
 68 
 69     if(angle < eps || fabs(cross(origin,a,b)) < eps)
 70         return 0;
 71     else if(disa < r+eps && disb < r+eps) { // 线段在圆内或圆上
 72         ret = fabs(cross(origin,a,b))/2.0;
 73     }
 74     else if(disa < r-eps || disb < r-eps) { // 线段的一端在圆内
 75         if(disb < r-eps) {
 76             swap(a,b);
 77             swap(disa,disb);
 78         }
 79         double disab = getdis(a,b);
 80         point mov = point((b.x-a.x)/disab,(b.y-a.y)/disab);
 81         double len = sqrt(r*r-disinter*disinter);
 82         point interpoint;
 83         interpoint = point(inter.x+mov.x*len, inter.y+mov.y*len);
 84         double angle0 = acos(dot(b,interpoint)/disb/r);
 85         ret = r*r*angle0/2.0;
 86         ret += fabs(cross(origin,interpoint,a))/2.0;
 87     }
 88     else {
 89         int flag = on_line(inter,a,b);
 90         ret = r*r*angle/2.0;
 91         if(flag && disinter<r)
 92         {
 93             double tangle = 2*acos(disinter/r);
 94             ret -= r*r*tangle/2.0;
 95             ret += r*disinter*sin(tangle/2.0);
 96         }
 97     }
 98     return ret*sgn;
 99 }
100 double solve(point center)
101 {
102     double ret = 0;
103     for(int i=1;i<=n;i++) {
104         point a,b;
105         a.x = p[i].x - center.x;
106         a.y = p[i].y - center.y;
107         b.x = p[i+1].x - center.x;
108         b.y = p[i+1].y - center.y;
109         ret += get_area(a,b);
110     }
111     return ret;
112 }
113 int main(void)
114 {
115     origin = point(0,0);
116     while(~scanf("%lf%d",&r,&n))
117     {
118         for(int i=1;i<=n;i++)
119         {
120             scanf("%lf%lf",&p[i].x,&p[i].y);
121             p[i].x++; // +1是为了测试一下代码
122             p[i].y++;
123         }
124         p[n+1]=p[1];
125         point tmp;
126         tmp.x = 1;
127         tmp.y = 1;
128         printf("%.2lf\n",fabs(solve(tmp)));
129     }
130     return 0;
131 }

 

posted @ 2014-05-07 10:10  辛力啤  阅读(223)  评论(0)    收藏  举报