计算几何模板(白皮书)
2013-08-14 14:42
399 查看
const double eps=1e-8;//精度 const int INF=0x3f3f3f3f; const double PI=acos(-1.0); inline int dcmp(const double& x) //判断double等于0或。。。 { if(fabs(x)<eps)return 0; else return x<0?-1:1; } struct Point { double x,y; Point() {} Point(double x,double y):x(x),y(y) {} }; typedef Point Vector; typedef vector<Point> Polygon; inline Vector operator+(const Vector& a,const Vector& b) { return Vector(a.x+b.x,a.y+b.y); //向量+向量=向量 } inline Vector operator-(const Point& a,const Point& b) { return Vector(a.x-b.x,a.y-b.y); //点-点=向量 } inline Vector operator*(const Vector& a,const double& p) { return Vector(a.x*p,a.y*p); //向量*实数=向量 } inline Vector operator/(const Vector& a,const double& p) { return Vector(a.x/p,a.y/p); //向量/实数=向量 } inline 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); } inline bool operator==(const Point&a,const Point&b) { return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0; } inline bool operator!=(const Point&a,const Point&b) { return a==b?false:true; } struct Segment { Point a,b; Segment() {} Segment(Point _a,Point _b) { a=_a,b=_b; } inline bool friend operator<(const Segment& p,const Segment& q) { return p.a<q.a||(p.a==q.a&&p.b<q.b); } inline bool friend operator==(const Segment& p,const Segment& q) { return (p.a==q.a&&p.b==q.b)||(p.a==q.b&&p.b==q.a); } }; struct Circle { Point c; double r; Circle() {} Circle(Point _c, double _r):c(_c),r(_r) {} Point point(double a)const { return Point(c.x+cos(a)*r,c.y+sin(a)*r); } bool friend operator<(const Circle& a,const Circle& b) { return a.r<b.r; } }; struct Line { Point p; Vector v; double ang; Line() {} Line(const Point &_p, const Vector &_v):p(_p),v(_v) { ang = atan2(v.y, v.x); } inline bool operator<(const Line &L)const { return ang < L.ang; } }; inline double Dot(const Vector& a,const Vector& b) { return a.x*b.x+a.y*b.y; //|a|*|b|*cosθ 点积 } inline double Length(const Vector& a) { return sqrt(Dot(a,a)); //|a| 向量长度 } inline double Angle(const Vector& a,const Vector& b) { return acos(Dot(a,b)/Length(a)/Length(b)); //向量夹角θ } inline double Cross(const Vector& a,const Vector& b) { return a.x*b.y-a.y*b.x; //叉积 向量围成的平行四边形的面积 } inline double Area2(const Point& a,const Point& b,Point c) { return Cross(b-a,c-a); //同上 参数为三个点 } inline double DegreeToRadius(const double& deg) { return deg/180*PI; } inline double GetRerotateAngle(const Vector& a,const Vector& b) //向量a顺时针旋转theta度得到向量b的方向 { double tempa=Angle(a,Vector(1,0)); if(a.y<0) tempa=2*PI-tempa; double tempb=Angle(b,Vector(1,0)); if(b.y<0) tempb=2*PI-tempb; if((tempa-tempb)>0) return tempa-tempb; else return tempa-tempb+2*PI; } inline double torad(const double& deg) { return deg/180*PI; //角度化为弧度 } inline Vector Rotate(const Vector& a,const double& rad) //向量逆时针旋转rad弧度 { return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad)); } inline Vector Normal(const Vector& a) //计算单位法线 { double L=Length(a); return Vector(-a.y/L,a.x/L); } inline Point GetLineProjection(const Point& p,const Point& a,const Point& b) //点在直线上的投影 { Vector v=b-a; return a+v*(Dot(v,p-a)/Dot(v,v)); } inline Point GetLineIntersection(Point p,Vector v,Point q,Vector w) //求直线交点 有唯一交点时可用 { Vector u=p-q; double t=Cross(w,u)/Cross(v,w); return p+v*t; } int ConvexHull(Point* p,int n,Point* sol) //计算凸包 { sort(p,p+n); int m=0; for(int i=0; i<n; i++) { while(m>1&&dcmp(Cross(sol[m-1]-sol[m-2],p[i]-sol[m-2]))<=0) m--; sol[m++]=p[i]; } int k=m; for(int i=n-2; i>=0; i--) { while(m>k&&dcmp(Cross(sol[m-1]-sol[m-2],p[i]-sol[m-2]))<=0) m--; sol[m++]=p[i]; } if(n>0) m--; return m; } double Heron(double a,double b,double c) //海伦公式 { double p=(a+b+c)/2; return sqrt(p*(p-a)*(p-b)*(p-c)); } bool SegmentProperIntersection(const Point& a1,const Point& a2,const Point& b1,const Point& b2) //线段规范相交判定 { double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1); double c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1); return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0; } double CutConvex(const int& n,Point* poly,const Point& a,const Point& b, vector<Point> result[3]) //有向直线a b 切割凸多边形 { vector<Point> points; Point p; Point p1=a,p2=b; int cur,pre; result[0].clear(); result[1].clear(); result[2].clear(); if(n==0) return 0; double tempcross; tempcross=Cross(p2-p1,poly[0]-p1); if(dcmp(tempcross)==0) pre=cur=2; else if(tempcross>0) pre=cur=0; else pre=cur=1; for(int i=0; i<n; i++) { tempcross=Cross(p2-p1,poly[(i+1)%n]-p1); if(dcmp(tempcross)==0) cur=2; else if(tempcross>0) cur=0; else cur=1; if(cur==pre) { result[cur].push_back(poly[(i+1)%n]); } else { p1=poly[i]; p2=poly[(i+1)%n]; p=GetLineIntersection(p1,p2-p1,a,b-a); points.push_back(p); result[pre].push_back(p); result[cur].push_back(p); result[cur].push_back(poly[(i+1)%n]); pre=cur; } } sort(points.begin(),points.end()); if(points.size()<2) { return 0; } else { return Length(points.front()-points.back()); } } double DistanceToSegment(Point p,Segment s) //点到线段的距离 { if(s.a==s.b) return Length(p-s.a); Vector v1=s.b-s.a,v2=p-s.a,v3=p-s.b; if(dcmp(Dot(v1,v2))<0) return Length(v2); else if(dcmp(Dot(v1,v3))>0) return Length(v3); else return fabs(Cross(v1,v2))/Length(v1); } inline bool isPointOnSegment(const Point& p,const Segment& s) { return dcmp(Cross(s.a-p,s.b-p))==0&&dcmp(Dot(s.a-p,s.b-p))<0; } int isPointInPolygon(Point p, Point* poly,int n) //点与多边形的位置关系 { int wn=0; for(int i=0; i<n; i++) { Point& p2=poly[(i+1)%n]; if(isPointOnSegment(p,Segment(poly[i],p2))) return -1;//点在边界上 int k=dcmp(Cross(p2-poly[i],p-poly[i])); int d1=dcmp(poly[i].y-p.y); int d2=dcmp(p2.y-p.y); if(k>0&&d1<=0&&d2>0)wn++; if(k<0&&d2<=0&&d1>0)wn--; } if(wn) return 1;//点在内部 else return 0;//点在外部 } double PolygonArea(Point* p,int n) //多边形有向面积 { 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; } int GetLineCircleIntersection(Line L,Circle C,Point& p1,Point& p2) //圆与直线交点 返回交点个数 { double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y-C.c.y; double e = a*a + c*c, f = 2*(a*b+c*d), g = b*b + d*d -C.r*C.r; double delta = f*f - 4*e*g; if(dcmp(delta) < 0) return 0;//相离 if(dcmp(delta) == 0) //相切 { p1=p1=L.p+L.v*(-f/(2*e)); return 1; }//相交 p1=(L.p+L.v*(-f-sqrt(delta))/(2*e)); p2=(L.p+L.v*(-f+sqrt(delta))/(2*e)); return 2; } double rotating_calipers(Point *ch,int n)//旋转卡壳 返回凸包中距离最大的两点的距离的平方 { int q=1; double ans=0; ch =ch[0]; for(int p=0; p<n; p++) { while(Cross(ch[q+1]-ch[p+1],ch[p]-ch[p+1])>Cross(ch[q]-ch[p+1],ch[p]-ch[p+1])) q=(q+1)%n; ans=max(ans,max(Length(ch[p]-ch[q]),Length(ch[p+1]-ch[q+1]))); } return ans; } Polygon CutPolygon(Polygon poly,const Point& a,const Point& b) //用a->b切割多边形 返回左侧 { Polygon newpoly; int n=poly.size(); for(int i=0; i<n; i++) { Point c=poly[i]; Point d=poly[(i+1)%n]; if(dcmp(Cross(b-a,c-a))>=0) newpoly.push_back(c); if(dcmp(Cross(b-a,c-d))!=0) { Point ip=GetLineIntersection(a,b-a,c,d-c); if(isPointOnSegment(ip,Segment(c,d))) newpoly.push_back(ip); } } return newpoly; } int GetCircleCircleIntersection(Circle c1,Circle c2,Point& p1,Point& p2) //求两圆相交 { double d=Length(c1.c-c2.c); if(dcmp(d)==0) { if(dcmp(c1.r-c2.r)==0) return -1;//两圆重合 return 0; } if(dcmp(c1.r+c2.r-d)<0) return 0; if(dcmp(fabs(c1.r-c2.r)-d)>0) return 0; double a=Angle(c2.c-c1.c,Vector(1,0)); double da=acos((c1.r*c1.r+d*d-c2.r*c2.r)/(2*c1.r*d)); p1=c1.point(a-da); p2=c1.point(a+da); if(p1==p2) return 1; return 2; } inline bool isPointOnleft(Point p,Line L) { return dcmp(Cross(L.v,p-L.p))>0; //点在直线左边 线上不算,算的话要>= } int HalfplaneIntersection(Line *L,int n,Point* poly) //半平面交 { sort(L,L+n); int first,last; Point* p=new Point ; Line* q=new Line ; q[first=last=0]=L[0]; for(int i=1; i<n; i++) { while(first<last&&!isPointOnleft(p[last-1],L[i])) last--; while(first<last&&!isPointOnleft(p[first],L[i])) first++; q[++last]=L[i]; if(dcmp(Cross(q[last].v,q[last-1].v))==0) { last--; if(isPointOnleft(L[i].p,q[last])) q[last]=L[i]; } if(first<last) p[last-1]=GetLineIntersection(q[last-1].p,q[last-1].v,q[last].p,q[last].v); } while(first<last&&!isPointOnleft(p[last-1],q[first])) last--; if(last-first<=1) return 0; p[last]=GetLineIntersection(q[last].p,q[last].v,q[first].p,q[first].v); int m=0; for(int i=first; i<=last; i++) poly[m++]=p[i]; return m; } //两点式化为一般式A = b.y-a.y, B = a.x-b.x, C = -a.y*(B)-a.x*(A);,点与圆的切线、圆与圆的共切线、三角形的外接圆和内接圆在白皮书p266
相关文章推荐