LightOJ 1130 Intersection between Circle and Rectangle

题目链接 http://lightoj.com/volume_showproblem.php?problem=1130

题意：求所给出的圆与矩形的重叠部分的面积

思路：构造出所求图形，求出矩形与圆的交点，将阴影部分分成若干三角形和扇形，分块计算面积。题目对计算精度要求很高

这个题因为精度问题纠结了好久……最后参考了一下别人的思路。

参考博文 http://blog.csdn.net/pvpishard/article/details/7220876

关于参考博文中提到的奇葩问题，我测试的时候全是WA……应该是LightOJ的测试数据更新了吧。

关于参考博文中提到“最精简的求1/sqrt()”的函数貌似精度只有小数点后3位，应该是不满足本题要求（小数点后6位）。

这个函数我也试着了解了一下，在较低精度要求下真的快不少！

参考链接  http://game.chinaitlab.com/devdoc/723647.html

                  http://www.cnblogs.com/nihgwu/archive/2006/08/04/467597.html

这个题应该还有更好的思路，LightOJ上目前最小MEMORY占用居然只有228……欢迎留言，一起讨论

#include <stdio.h>#include <math.h>const double PI=acos(-1.0);const double STD=1e-9;      //浮点数大小判断标准struct Point{double x,y;void disp (){printf("x=%lf y=%lf\n",x,y);}Point friend Vec (Point a,Point b)   //构造向量{Point res;res.x=a.x-b.x;res.y=a.y-b.y;return res;}}pt[20];struct Line{Point s,e;void set (double sx,double sy,double ex,double ey){s.x=sx;s.y=sy;e.x=ex;e.y=ey;}}line[4];struct Circle{Point c;double r;}cir;int point_num;       //记录需要计算的点的数量bool out[20];        //标记点与圆的位置关系,在圆内和圆上为false,在圆外为truevoid Find (Line l,int mode)  //寻找矩形与圆的交点,结果保存在pt数组中{Point temp;double t,d;switch (mode){case 0:                 //寻找"左至右"边与圆的交点if (fabs(cir.c.y-l.s.y) > cir.r)   //与圆不相交return;d=cir.c.y-l.s.y;t=sqrt(fabs(cir.r*cir.r-d*d));temp.y=l.s.y;temp.x=cir.c.x-t;if (temp.x>l.s.x && temp.x<l.e.x)pt[point_num++]=temp;temp.x=cir.c.x+t;if (temp.x>l.s.x && temp.x<l.e.x)pt[point_num++]=temp;break;case 1:                 //寻找"上至下"边与圆的交点if (fabs(cir.c.x-l.s.x) > cir.r)   //与圆不相交return;d=cir.c.x-l.s.x;t=sqrt(fabs(cir.r*cir.r-d*d));temp.x=l.s.x;temp.y=cir.c.y+t;if (temp.y<l.s.y && temp.y>l.e.y)pt[point_num++]=temp;temp.y=cir.c.y-t;if (temp.y<l.s.y && temp.y>l.e.y)pt[point_num++]=temp;break;case 2:                   //寻找"右至左"边与圆的交点if (fabs(cir.c.y-l.s.y) > cir.r)     //与圆不相交return;d=cir.c.y-l.s.y;t=sqrt(fabs(cir.r*cir.r-d*d));temp.y=l.s.y;temp.x=cir.c.x+t;if (temp.x<l.s.x && temp.x>l.e.x)pt[point_num++]=temp;temp.x=cir.c.x-t;if (temp.x<l.s.x&&temp.x>l.e.x)pt[point_num++]=temp;break;case 3:                //寻找"下至上"边与圆的交点if (fabs(cir.c.x-l.s.x) > cir.r)    //与圆不相交return;d=cir.c.x-l.s.x;t=sqrt(fabs(cir.r*cir.r-d*d));temp.x=l.s.x;if (l.s.y<l.e.y)    //寻找"下至上"边与圆的交点temp.y=cir.c.y-t;if (temp.y>l.s.y && temp.y<l.e.y)pt[point_num++]=temp;temp.y=cir.c.y+t;if (temp.y>l.s.y && temp.y<l.e.y)pt[point_num++]=temp;}}/*            void Find (Line l)   //网上别人代码的修改+注释版{Point temp;double t,d;if (fabs(l.s.x-l.e.x) < STD) //起始横坐标相同,即竖直线{if (fabs(cir.c.x-l.s.x) > cir.r) //起点在圆外return;d=cir.c.x-l.s.x;t=sqrt(fabs(cir.r*cir.r-d*d));temp.x=l.s.x;if (l.s.y<l.e.y) //寻找"下至上"边与圆的交点{temp.y=cir.c.y-t;if (temp.y>l.s.y && temp.y<l.e.y)pt[point_num++]=temp;temp.y=cir.c.y+t;if (temp.y>l.s.y && temp.y<l.e.y)pt[point_num++]=temp;}else  //寻找"上至下"边与圆的交点{temp.y=cir.c.y+t;if (temp.y<l.s.y && temp.y>l.e.y)pt[point_num++]=temp;temp.y=cir.c.y-t;if (temp.y<l.s.y && temp.y>l.e.y)pt[point_num++]=temp;}}else    //起始横坐标不同,即水平线{if (fabs(cir.c.y-l.s.y) > cir.r)     //起点在圆外return;d=cir.c.y-l.s.y;t=sqrt(fabs(cir.r*cir.r-d*d));temp.y=l.s.y;if (l.s.x<l.e.x)          //寻找"左至右"边与圆的交点{temp.x=cir.c.x-t;if (temp.x>l.s.x && temp.x<l.e.x)pt[point_num++]=temp;temp.x=cir.c.x+t;if (temp.x>l.s.x && temp.x<l.e.x)pt[point_num++]=temp;}else       //寻找"右至左"边与圆的交点{temp.x=cir.c.x+t;if (temp.x<l.s.x && temp.x>l.e.x)pt[point_num++]=temp;temp.x=cir.c.x-t;if (temp.x<l.s.x&&temp.x>l.e.x)pt[point_num++]=temp;}}}*/void Input ()         //输入及构造部分{scanf("%lf%lf%lf",&cir.c.x,&cir.c.y,&cir.r);double lx,ly,rx,ry;scanf("%lf%lf%lf%lf",&lx,&ly,&rx,&ry);line[0].set(lx,ry,rx,ry);        //上边-左至右line[1].set(rx,ry,rx,ly);        //右边-上至下line[2].set(rx,ly,lx,ly);        //下边-右至左line[3].set(lx,ly,lx,ry);        //左边-下至上//求出需要考虑的点,即矩形顶点以及矩形与圆的交点point_num=0;pt[point_num++]=line[0].s;Find (line[0],0);pt[point_num++]=line[1].s;Find (line[1],1);pt[point_num++]=line[2].s;Find (line[2],2);pt[point_num++]=line[3].s;Find (line[3],3);//for (int i=0;i<point_num;i++) printf("Point %d ",i),pt[i].disp();}double get_dis (Point a,Point b)        //计算两点间距离{return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}double get_dis2 (Point a,Point b)       //计算两点间距离的平方{return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);}void Judge (){double t;for (int i=0;i<point_num;i++){t=get_dis(pt[i],cir.c);if (t>cir.r+STD)out[i]=true;else out[i]=false;}}double get_cj (Point a,Point b)        //计算叉积{return a.x*b.y-a.y*b.x;}double get_dj(Point a,Point b)        //计算点积{return a.x*b.x+a.y*b.y;}double get_tri (Point a,Point b)     //计算三角形面积{return get_cj (Vec(a,cir.c),Vec(b,cir.c))/2;}double get_sec (Point a,Point b)     //计算扇形面积{double da=get_dis2 (a,cir.c);double db=get_dis2 (b,cir.c);double cj=get_cj (Vec(a,cir.c),Vec(b,cir.c));double dj=get_dj (Vec(a,cir.c),Vec(b,cir.c));double dsin=cj/sqrt(da*db);       //da,db是平方过的,如果先开方再乘会损失精度WA……double alpha=asin(dsin);//asin函数求出的角是0到90度的,需要根据实际情况还原到原来的大小if (dj < -STD)           //当夹角大于90度时if (alpha>0)alpha=PI-alpha;elsealpha=-PI-alpha;return alpha*cir.r*cir.r/2;}double Deal (){if (cir.r < STD)return 0;Judge ();double ans=0;for (int i=0;i<point_num;i++){if (out[i] == false && out[(i+1)%point_num] == false)//如果两个点都不在圆外ans+=get_tri(pt[i],pt[(i+1)%point_num]);else  //若至少有一点在圆外ans+=get_sec(pt[i],pt[(i+1)%point_num]);}return fabs(ans);       //用叉积可能负面积}int main (){int T;scanf("%d",&T);for (int cas=1;cas<=T;cas++){Input ();printf("Case %d: %.8lf\n",cas,Deal());}return 0;}/*24 4 41 1 3 34 4 41 1 7 7OutCase 1: 3.93987658Case 2: 35.75950633*/