三角剖分·圆和多边形的交

三角剖分·圆和多边形的交

POJ 3675 Telescope

​​http://poj.org/problem?id=3675​​

大意：求解圆和多边形的交。

任意一个凸N多边形均可分解成N-2个三角形。因此，这就是讨论分解后的三角形和圆的交的问题。

它有这些情况:

（1）：

（2）：

（3）:

（4）:

(5):

5)又可分为p1在外，p2在里；p1在里，p2在外。

在写代码的过程中遇到了一些无法解决的问题，特别是5.2情况，参考了 ​​​#include #include 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)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&&x1p2.x&&x1p1.y&&y1p2.y&&y1=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&&x1p2.x&&x1p1.y&&y1p2.y&&y1

POJ 2986 A Triangle and a Circle

​​#include #include 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)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&&x1p2.x&&x1p1.y&&y1p2.y&&y1=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&&x1p2.x&&x1p1.y&&y1p2.y&&y1

当我使用模板后,主函数内这样写：

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#include#include#include#includeconst 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<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<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 &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 &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.r0) // 外离 又有两条外公切线 { 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 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 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;}

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。