三角剖分·圆和多边形的交
2016-02-15 17:09
281 查看
POJ 3675 Telescope
http://poj.org/problem?id=3675
大意:求解圆和多边形的交。
分析:任意一个凸N多边形均可分解成N-2个三角形。因此,这就是讨论分解后的三角形和圆的交的问题。
它有这些情况:
(1):
![](https://img-blog.csdn.net/20160215170040264?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
(2):
![](https://img-blog.csdn.net/20160215170109764?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
(3):
![](https://img-blog.csdn.net/20160215170141155?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
(4):
![](https://img-blog.csdn.net/20160215170227854?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
(5):
![](https://img-blog.csdn.net/20160215170253980?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
5)又可分为p1在外,p2在里;p1在里,p2在外。
在写代码的过程中遇到了一些无法解决的问题,特别是5.2情况,参考了
http://gzhu-101majia.iteye.com/blog/1128345
不能理解为什么自己的做法不能通过。
求解p1p2直线和圆的交点(假设不与横轴垂直,可参照(5)图):
![](http://latex.codecogs.com/gif.latex?%5C%5C%20k=%5Cfrac%7By_2-y_1%7D%7Bx_2-x_1%7D%5C%5C%20y_3=y_1+k(x_3-x_1)%5C%5C%20x_3%5E2+%5By_1+k(x_3-x_1)%5D%5E2=r%5E2%5C%5C%20(1+k%5E2)x_3%5E2+(2y_1k-2k%5E2x_1)x_3+k%5E2x_1%5E2+y_1%5E2-2y_1kx_1-r%5E2=0%5C%5C%20%5Clongrightarrow%20Ax_3%5E2+Bx_3+C=0)
分解后的三角形和圆的交按有向面积累加:
![](https://img-blog.csdn.net/20160215170454265?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQv/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/Center)
POJ 2986 A Triangle and a Circle
http://poj.org/problem?id=2986
求解三角形和圆的交。圆心是可变的。但是三角形只有一个。
再一次用此方法,但是看似更加简单的这题居然过不了,无语啊,太打击人的积极性了吧!!
WA:
当我使用模板后,主函数内这样写:
输出0. 不能直接给x,y赋值,对于long double数据成员的点%Lf输入居然也输出0。只能利用构造函数来初始化了
主函数内:
代码:
http://poj.org/problem?id=3675
大意:求解圆和多边形的交。
分析:任意一个凸N多边形均可分解成N-2个三角形。因此,这就是讨论分解后的三角形和圆的交的问题。
它有这些情况:
(1):
(2):
(3):
(4):
(5):
5)又可分为p1在外,p2在里;p1在里,p2在外。
在写代码的过程中遇到了一些无法解决的问题,特别是5.2情况,参考了
http://gzhu-101majia.iteye.com/blog/1128345
不能理解为什么自己的做法不能通过。
求解p1p2直线和圆的交点(假设不与横轴垂直,可参照(5)图):
分解后的三角形和圆的交按有向面积累加:
#include <iostream> #include <stdio.h> #include <cmath> using namespace std; const double eps=1e-7,PI=acos(-1.0); struct point { double x,y; }p[55]; double r; double pp_dis(point p1,point p2) //两点距离 { return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); } double xmulti(point p0,point p1,point p2){ return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y); } double angle(point p1,point o,point p2){ double a=pp_dis(o,p1); double b=pp_dis(o,p2); double c=pp_dis(p1,p2); double cos=(a*a+b*b-c*c)/(2*a*b); return fabs(acos(cos)); } double dir(point o,point p1,point p2){ double area=xmulti(o,p1,p2); if(area<-eps) return 1.0; return -1.0; } point intersec(double x,double y) //计算直线与圆的交点 { double k; point temp; if(x!=0) { k=y/x; temp.x=fabs(r)/sqrt(1+k*k); if(x<0) temp.x=-temp.x; temp.y=k*temp.x; } else { temp.x=0; if(y>0) temp.y=r; else temp.y=-r; } return temp; } double get_area(point o,point p1,point p2) //三角剖分 { double f=dir(o,p1,p2); //判断三角形面积加还是减 double s; double a=pp_dis(o,p1); double b=pp_dis(o,p2); double c=pp_dis(p1,p2); double d=fabs(xmulti(o,p1,p2))/c; // 1 if(a<=r && b<=r) { double area=xmulti(o,p1,p2); return fabs(area)/2.0*f; } // 2 else if(a>=r && b>=r && d>=r) { double sita1=angle(p1,o,p2); double s=fabs(sita1*r*r/2.0); //扇形s=θ*r*r/2 return s*f; } // 3 else if(a>=r && b>=r && d<=r && (angle(o,p1,p2)-PI/2>eps || angle(o,p2,p1)-PI/2>eps)) { double sita=angle(p1,o,p2); s=fabs(sita*r*r/2.0); return s*f; } // 4 else if(a>=r && b>=r && d<=r && angle(o,p1,p2)<=PI/2 && angle(o,p2,p1)<=PI/2) { point p3,p4; if(fabs(p1.x-p2.x)>eps){ double k=(p2.y-p1.y)/(p2.x-p1.x); double A=1+k*k; double B=2*p1.y*k-2*k*k*p1.x; double C=k*k*p1.x*p1.x+p1.y*p1.y-2*p1.y*k*p1.x-r*r; double x1=(-B-sqrt(B*B-4*A*C))/(2*A); double x2=(-B+sqrt(B*B-4*A*C))/(2*A); if(fabs(x1-p1.x)<fabs(x2-p1.x)){ p3.x=x1; p4.x=x2; } else { p3.x=x2; p4.x=x1; } p3.y=p1.y+k*(p3.x-p1.x); p4.y=p1.y+k*(p4.x-p1.x); } else { p3.x=p1.x; p4.x=p1.x; double h=sqrt(r*r-d*d); if(p1.y>eps){ p3.y=h; p4.y=-h; } else { p3.y=-h; p4.y=h; } } double ag1=angle(p1,o,p3); double ag2=angle(p4,o,p2); double area=0; area+=r*r*(ag1+ag2)/2; area+=fabs(xmulti(p3,o,p4))/2.0; return area*f; } // 5.1 else if(a>=r && b<=r) { double area=0; point p3; if(fabs(p1.x-p2.x)>eps){ double k=(p2.y-p1.y)/(p2.x-p1.x); double A=1+k*k; double B=-2*p1.x*k*k+2*p1.y*k; double C=p1.y*p1.y+k*k*p1.x*p1.x-2*p1.y*k*p1.x-r*r; double x1=(-B+sqrt(B*B-4*A*C))/(2*A); double x2=(-B-sqrt(B*B-4*A*C))/(2*A); if( (x1>p1.x&&x1<p2.x) || (x1>p2.x&&x1<p1.x)) p3.x=x1; else p3.x=x2; p3.y=p1.y+k*(p3.x-p1.x); } else { p3.x=p1.x; double y1=sqrt(r*r-p3.x*p3.x); double y2=-y1; if( (y1>p1.y&&y1<p2.y) || (y1>p2.y&&y1<p1.y)) p3.y=y1; else p3.y=y2; } area+=fabs(xmulti(p3,o,p2))/2.0; double ag=angle(p3,o,p1); area+=r*r*ag/2; return area*f; } // 5.2 else if(a<=r && b>=r) //a短于半径,b长于半径 { double area=0; point p3; if(fabs(p1.x-p2.x)>eps) { double k=(p2.y-p1.y)/(p2.x-p1.x); double A=1+k*k; double B=-2*p1.x*k*k+2*p1.y*k; double C=p1.y*p1.y+k*k*p1.x*p1.x-2*p1.y*k*p1.x-r*r; double x1=(-B+sqrt(B*B-4*A*C))/(2*A); double x2=(-B-sqrt(B*B-4*A*C))/(2*A); if( (x1>p1.x&&x1<p2.x) || (x1>p2.x&&x1<p1.x)) p3.x=x1; else p3.x=x2; p3.y=p1.y+k*(p3.x-p1.x); } else { p3.x=p1.x; double y1=sqrt(r*r-p3.x*p3.x); double y2=-y1; if( (y1>p1.y&&y1<p2.y) || (y1>p2.y&&y1<p1.y)) p3.y=y1; else p3.y=y2; } //double sita1=angle(p2,o,p3); point pp2=intersec(p2.x,p2.y); double c=pp_dis(pp2,p3); double cos=(r*r+r*r-c*c)/(2*r*r); double sita1=acos(cos); //让我想不通的是为什么求角度非得用这种方法,而上面更加直接的方法 double s1=fabs(sita1*r*r/2.0); // 却错误。 double s3=fabs(p3.x*p1.y-p3.y*p1.x)/2.0; area=s1+s3; return area*f; } else return 0; } int main() { //freopen("cin.txt","r",stdin); while(~scanf("%lf",&r)){ double ans=0; point o; o.x=o.y=0; int n; scanf("%d",&n); for(int i=0;i<n;i++){ scanf("%lf%lf",&p[i].x,&p[i].y); } p =p[0]; ans=0; for(int i=0;i<n;i++) { ans+=get_area(o,p[i],p[i+1]); } printf("%.2lf\n",fabs(ans)); } return 0; }
POJ 2986 A Triangle and a Circle
http://poj.org/problem?id=2986
求解三角形和圆的交。圆心是可变的。但是三角形只有一个。
再一次用此方法,但是看似更加简单的这题居然过不了,无语啊,太打击人的积极性了吧!!
WA:
#include <iostream> #include <stdio.h> #include <cmath> using namespace std; const double eps=1e-7,PI=acos(-1.0); struct point{ double x,y; }; double pp_dis(point p1,point p2) //两点距离 { return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); } double xmulti(point p0,point p1,point p2){ return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y); } double angle(point p1,point o,point p2){ double a=pp_dis(o,p1); double b=pp_dis(o,p2); double c=pp_dis(p1,p2); double cos=(a*a+b*b-c*c)/(2*a*b); return fabs(acos(cos)); } double dir(point o,point p1,point p2){ double area=xmulti(o,p1,p2); if(area<-eps) return 1.0; return -1.0; } point intersec(double r,double x,double y) //计算直线与圆的交点 { double k; point temp; if(x!=0) { k=y/x; temp.x=fabs(r)/sqrt(1+k*k); if(x<0) temp.x=-temp.x; temp.y=k*temp.x; } else { temp.x=0; if(y>0) temp.y=r; else temp.y=-r; } return temp; } double get_area(double r,point o,point p1,point p2) //三角剖分 { double f=dir(o,p1,p2); //判断三角形面积加还是减 double s; double a=pp_dis(o,p1); double b=pp_dis(o,p2); double c=pp_dis(p1,p2); double d=fabs(xmulti(o,p1,p2))/c; // 1 if(a<=r && b<=r) { double area=xmulti(o,p1,p2); return fabs(area)/2.0*f; } // 2 else if(a>=r && b>=r && d>=r) { double sita1=angle(p1,o,p2); double s=fabs(sita1*r*r/2.0); //扇形s=θ*r*r/2 return s*f; } // 3 else if(a>=r && b>=r && d<=r && (angle(o,p1,p2)-PI/2>eps || angle(o,p2,p1)-PI/2>eps)) { double sita=angle(p1,o,p2); s=fabs(sita*r*r/2.0); return s*f; } // 4 else if(a>=r && b>=r && d<=r && angle(o,p1,p2)<=PI/2 && angle(o,p2,p1)<=PI/2) { point p3,p4; if(fabs(p1.x-p2.x)>eps){ double k=(p2.y-p1.y)/(p2.x-p1.x); double A=1+k*k; double B=2*p1.y*k-2*k*k*p1.x; double C=k*k*p1.x*p1.x+p1.y*p1.y-2*p1.y*k*p1.x-r*r; double x1=(-B-sqrt(B*B-4*A*C))/(2*A); double x2=(-B+sqrt(B*B-4*A*C))/(2*A); if(fabs(x1-p1.x)<fabs(x2-p1.x)){ p3.x=x1; p4.x=x2; } else { p3.x=x2; p4.x=x1; } p3.y=p1.y+k*(p3.x-p1.x); p4.y=p1.y+k*(p4.x-p1.x); } else { p3.x=p1.x; p4.x=p1.x; double h=sqrt(r*r-d*d); if(p1.y>eps){ p3.y=h; p4.y=-h; } else { p3.y=-h; p4.y=h; } } double ag1=angle(p1,o,p3); double ag2=angle(p4,o,p2); double area=0; area+=r*r*(ag1+ag2)/2; area+=fabs(xmulti(p3,o,p4))/2.0; return area*f; } // 5.1 else if(a>=r && b<=r) { double area=0; point p3; if(fabs(p1.x-p2.x)>eps){ double k=(p2.y-p1.y)/(p2.x-p1.x); double A=1+k*k; double B=-2*p1.x*k*k+2*p1.y*k; double C=p1.y*p1.y+k*k*p1.x*p1.x-2*p1.y*k*p1.x-r*r; double x1=(-B+sqrt(B*B-4*A*C))/(2*A); double x2=(-B-sqrt(B*B-4*A*C))/(2*A); if( (x1>p1.x&&x1<p2.x) || (x1>p2.x&&x1<p1.x)) p3.x=x1; else p3.x=x2; p3.y=p1.y+k*(p3.x-p1.x); } else { p3.x=p1.x; double y1=sqrt(r*r-p3.x*p3.x); double y2=-y1; if( (y1>p1.y&&y1<p2.y) || (y1>p2.y&&y1<p1.y)) p3.y=y1; else p3.y=y2; } area+=fabs(xmulti(p3,o,p2))/2.0; double ag=angle(p3,o,p1); area+=r*r*ag/2; return area*f; } // 5.2 else if(a<=r && b>=r) //a短于半径,b长于半径 { double area=0; point p3; if(fabs(p1.x-p2.x)>eps) { double k=(p2.y-p1.y)/(p2.x-p1.x); double A=1+k*k; double B=-2*p1.x*k*k+2*p1.y*k; double C=p1.y*p1.y+k*k*p1.x*p1.x-2*p1.y*k*p1.x-r*r; double x1=(-B+sqrt(B*B-4*A*C))/(2*A); double x2=(-B-sqrt(B*B-4*A*C))/(2*A); if( (x1>p1.x&&x1<p2.x) || (x1>p2.x&&x1<p1.x)) p3.x=x1; else p3.x=x2; p3.y=p1.y+k*(p3.x-p1.x); } else { p3.x=p1.x; double y1=sqrt(r*r-p3.x*p3.x); double y2=-y1; if( (y1>p1.y&&y1<p2.y) || (y1>p2.y&&y1<p1.y)) p3.y=y1; else p3.y=y2; } //double sita1=angle(p2,o,p3); point pp2=intersec(r,p2.x,p2.y); double c=pp_dis(pp2,p3); double cos=(r*r+r*r-c*c)/(2*r*r); double sita1=acos(cos); double s1=fabs(sita1*r*r/2.0); double s3=fabs(xmulti(o,p1,p3))/2.0; area=s1+s3; return area*f; } else return 0; } int main() { //freopen("cin.txt","r",stdin); point p[5],o; double r; while(~scanf("%lf",&p[1].x)){ scanf("%lf",&p[1].y); for(int i=2;i<=3;i++){ scanf("%lf %lf",&p[i].x,&p[i].y); } scanf("%lf %lf %lf",&o.x,&o.y,&r); for(int i=1;i<=3;i++){ p[i].x-=o.x; p[i].y-=o.y; } p[4]=p[1]; o.x=0; o.y=0; double ans=0; for(int i=1;i<=3;i++){ ans=ans+get_area(r,o,p[i],p[i+1]); } //printf("%.2f\n",fabs(ans)); printf("%.2f\n",fabs(ans)+eps); } return 0; }
当我使用模板后,主函数内这样写:
Point p[5],o; double r; while(~scanf("%lf",&p[1].x)){ scanf("%lf",&p[1].y); for(int i=2;i<=3;i++){ scanf("%lf %lf",&p[i].x,&p[i].y); } scanf("%lf %lf %lf",&o.x,&o.y,&r); for(int i=1;i<=3;i++){ p[i].x-=o.x; p[i].y-=o.y; } p[4]=p[1]; double ans=0; Circle C=Circle(Zero,r); for(int i=1;i<=3;i++) ans+=common_area(C, p[i], p[i+1]); printf("%.2lf\n",fabs(ans)+eps); } return 0;
输出0. 不能直接给x,y赋值,对于long double数据成员的点%Lf输入居然也输出0。只能利用构造函数来初始化了
主函数内:
Point p[5]; double a[9]; while(~scanf("%lf",&a[0])){ for(int i=1;i<9;i++){ scanf("%lf",&a[i]); } p[0]=Point(a[0]-a[6],a[1]-a[7]); p[1]=Point(a[2]-a[6],a[3]-a[7]); p[2]=Point(a[4]-a[6],a[5]-a[7]); p[3]=p[0]; double ans=0; Circle C=Circle(Zero,a[8]); for(int i=0;i<3;i++) ans+=common_area(C, p[i], p[i+1]); printf("%.2lf\n",fabs(ans)+eps); }
代码:
#include<iostream> #include<cmath> #include<cstdio> #include<algorithm> #include<vector> const double eps=1e-10; const long double PI=acos(-1.0); using namespace std; struct Point{ long double x; long double y; Point(long double x=0,long double y=0):x(x),y(y){} void operator<<(Point &A) {cout<<A.x<<' '<<A.y<<endl;} }; int dcmp(long double x) {return (x>eps)-(x<-eps); } int sgn(long double x) {return (x>eps)-(x<-eps); } typedef Point Vector; Vector operator +(Vector A,Vector B) { return Vector(A.x+B.x,A.y+B.y);} Vector operator -(Vector A,Vector B) { return Vector(A.x-B.x,A.y-B.y); } Vector operator *(Vector A,long double p) { return Vector(A.x*p,A.y*p); } Vector operator /(Vector A,long double p) {return Vector(A.x/p,A.y/p);} ostream &operator<<(ostream & out,Point & P) { out<<P.x<<' '<<P.y<<endl; return out;} // bool operator< (const Point &A,const Point &B) { return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0); } bool operator== ( const Point &A,const Point &B) { return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0;} long double Dot(Vector A,Vector B) {return A.x*B.x+A.y*B.y;} long double Cross(Vector A,Vector B) {return A.x*B.y-B.x*A.y; } long double Length(Vector A) { return sqrt(Dot(A, A));} long double Angle(Vector A,Vector B) {return acos(Dot(A,B)/Length(A)/Length(B));} long double Area2(Point A,Point B,Point C ) {return Cross(B-A, C-A);} Vector Rotate(Vector A,long double rad) { return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));} Vector Normal(Vector A) {long double L=Length(A);return Vector(-A.y/L,A.x/L);} Point GetLineIntersection(Point P,Vector v,Point Q,Vector w) { Vector u=P-Q; long double t=Cross(w, u)/Cross(v,w); return P+v*t; } long double DistanceToLine(Point P,Point A,Point B) { Vector v1=P-A; Vector v2=B-A; return fabs(Cross(v1,v2))/Length(v2); } long double DistanceToSegment(Point P,Point A,Point B) { if(A==B) return Length(P-A); Vector v1=B-A; Vector v2=P-A; Vector v3=P-B; if(dcmp(Dot(v1,v2))==-1) return Length(v2); else if(Dot(v1,v3)>0) return Length(v3); else return DistanceToLine(P, A, B); } Point GetLineProjection(Point P,Point A,Point B) { Vector v=B-A; Vector v1=P-A; long double t=Dot(v,v1)/Dot(v,v); return A+v*t; } bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2) { long double c1=Cross(b1-a1, a2-a1); long double c2=Cross(b2-a1, a2-a1); long double c3=Cross(a1-b1, b2-b1); long double c4=Cross(a2-b1, b2-b1); return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0 ; } bool OnSegment(Point P,Point A,Point B) { return dcmp(Cross(P-A, P-B))==0&&dcmp(Dot(P-A,P-B))<=0; } long double PolygonArea(Point *p,int n) { long double area=0; for(int i=1;i<n-1;i++) { area+=Cross(p[i]-p[0], p[i+1]-p[0]); } return area/2; } Point read_point() { Point P; scanf("%Lf%Lf",&P.x,&P.y); return P; } // ---------------与圆有关的-------- struct Circle { Point c; long double r; Circle(Point c=Point(0,0),long double r=0):c(c),r(r) {} Point point(long double a) { return Point(c.x+r*cos(a),c.y+r*sin(a)); } }; struct Line { Point p; Vector v; Line(Point p=Point(0,0),Vector v=Vector(0,1)):p(p),v(v) {} Point point(long double t) { return Point(p+v*t); } }; int getLineCircleIntersection(Line L,Circle C,long double &t1,long double &t2,vector<Point> &sol) { long double a=L.v.x; long double b=L.p.x-C.c.x; long double c=L.v.y; long double d=L.p.y-C.c.y; long double e=a*a+c*c; long double f=2*(a*b+c*d); long double g=b*b+d*d-C.r*C.r; long double delta=f*f-4*e*g; if(dcmp(delta)<0) return 0; if(dcmp(delta)==0) { t1=t2=-f/(2*e); sol.push_back(L.point(t1)); return 1; } else { t1=(-f-sqrt(delta))/(2*e); t2=(-f+sqrt(delta))/(2*e); sol.push_back(L.point(t1)); sol.push_back(L.point(t2)); return 2; } } // 向量极角公式 long double angle(Vector v) {return atan2(v.y,v.x);} int getCircleCircleIntersection(Circle C1,Circle C2,vector<Point> &sol) { long double d=Length(C1.c-C2.c); if(dcmp(d)==0) { if(dcmp(C1.r-C2.r)==0) return -1; // 重合 else return 0; // 内含 0 个公共点 } if(dcmp(C1.r+C2.r-d)<0) return 0; // 外离 if(dcmp(fabs(C1.r-C2.r)-d)>0) return 0; // 内含 long double a=angle(C2.c-C1.c); long double da=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d)); Point p1=C1.point(a-da); Point p2=C1.point(a+da); sol.push_back(p1); if(p1==p2) return 1; // 相切 else { sol.push_back(p2); return 2; } } // 求点到圆的切线 int getTangents(Point p,Circle C,Vector *v) { Vector u=C.c-p; long double dist=Length(u); if(dcmp(dist-C.r)<0) return 0; else if(dcmp(dist-C.r)==0) { v[0]=Rotate(u,PI/2); return 1; } else { long double ang=asin(C.r/dist); v[0]=Rotate(u,-ang); v[1]=Rotate(u,+ang); return 2; } } // 求两圆公切线 int getTangents(Circle A,Circle B,Point *a,Point *b) { int cnt=0; if(A.r<B.r) { swap(A,B); swap(a, b); // 有时需标记 } long double d=Length(A.c-B.c); long double rdiff=A.r-B.r; long double rsum=A.r+B.r; if(dcmp(d-rdiff)<0) return 0; // 内含 long double base=angle(B.c-A.c); if(dcmp(d)==0&&dcmp(rdiff)==0) return -1 ; // 重合 无穷多条切线 if(dcmp(d-rdiff)==0) // 内切 外公切线 { a[cnt]=A.point(base); b[cnt]=B.point(base); cnt++; return 1; } // 有外公切线的情形 long double ang=acos(rdiff/d); a[cnt]=A.point(base+ang); b[cnt]=B.point(base+ang); cnt++; a[cnt]=A.point(base-ang); b[cnt]=B.point(base-ang); cnt++; if(dcmp(d-rsum)==0) // 外切 有内公切线 { a[cnt]=A.point(base); b[cnt]=B.point(base+PI); cnt++; } else if(dcmp(d-rsum)>0) // 外离 又有两条外公切线 { long double ang_in=acos(rsum/d); a[cnt]=A.point(base+ang_in); b[cnt]=B.point(base+ang_in+PI); cnt++; a[cnt]=A.point(base-ang_in); b[cnt]=B.point(base-ang_in+PI); cnt++; } return cnt; } int n; Point Zero=Point(0,0); long double common_area(Circle C,Point A,Point B) { // if(A==B) return 0; if(A==C.c||B==C.c) return 0; long double OA=Length(A-C.c),OB=Length(B-C.c); long double d=DistanceToLine(Zero, A, B); int sg=sgn(Cross(A,B)); if(sg==0) return 0; long double angle=Angle(A,B); if(dcmp(OA-C.r)<=0&&dcmp(OB-C.r)<=0) { return Cross(A,B)/2; } else if(dcmp(OA-C.r)>=0&&dcmp(OB-C.r)>=0&&dcmp(d-C.r)>=0) { return sg*C.r*C.r*angle/2; } else if (dcmp(OA-C.r)>=0&&dcmp(OB-C.r)>=0&&dcmp(d-C.r)<0) { Point prj=GetLineProjection(Zero, A, B); if(OnSegment(prj, A, B)) { vector<Point> p; Line L=Line(A,B-A); long double t1,t2; getLineCircleIntersection(L,C, t1, t2, p); long double s1=0; s1=C.r*C.r*angle/2; long double s2=0; s2=C.r*C.r*Angle(p[0],p[1])/2; s2-=fabs(Cross(p[0],p[1])/2); s1=s1-s2; return sg*s1; } else { return sg*C.r*C.r*angle/2; } } else { if(dcmp(OB-C.r)<0) { Point temp=A; A=B; B=temp; } long double t1,t2; Line L=Line(A,B-A); vector<Point> inter; getLineCircleIntersection(L, C, t1, t2,inter); Point inter_point; if(OnSegment(inter[0], A, B)) inter_point=inter[0]; else { inter_point=inter[1]; } long double s=fabs(Cross(inter_point, A)/2); s+=C.r*C.r*Angle(inter_point,B)/2; return s*sg; } } int main() { //freopen("cin.txt","r",stdin); Point p[5]; double a[9]; while(~scanf("%lf",&a[0])){ for(int i=1;i<9;i++){ scanf("%lf",&a[i]); } p[0]=Point(a[0]-a[6],a[1]-a[7]); p[1]=Point(a[2]-a[6],a[3]-a[7]); p[2]=Point(a[4]-a[6],a[5]-a[7]); p[3]=p[0]; double ans=0; Circle C=Circle(Zero,a[8]); for(int i=0;i<3;i++) ans+=common_area(C, p[i], p[i+1]); printf("%.2lf\n",fabs(ans)+eps); // C++ AC G++ TLE } return 0; }
相关文章推荐
- My Knowledge Base
- HDU 2255 奔小康赚大钱(二分图最优匹配:模板题)
- 【证明】关于项链问题中旋转操作的置换群的循环节个数、长度的证明【Pólya】【同余】
- 环信聊天列表不刷新数据
- Android color集锦
- 青瓷引擎问答集锦(一)
- hdu 1226 超级密码【BFS】
- Ajax中Get请求与Post请求的区别
- 使用jQuery和CSS3制作数字时钟(CSS3篇)
- spring学习笔记(一)快速搭建一个spring boot应用
- Git版本控制:Git高级教程
- 开启博客之旅
- WebApp开发技巧精华总结
- Git高级教程
- Too many constants, the constant pool for ***_jsp would exceed 65536 entries
- flex添加 播放 gif图片
- 位运算及相关面试题汇总
- USACO 1.3 barn1
- install opencv3 with python2.7 bindings on Mac using Hombrew and git
- moveTaskToBack Android应用退到后台,类似按下Home键