计算几何模板

in 知识点 with 0 comment

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#include<algorithm>
using namespace std;
const int max_point=1e+5 +5;
const double Max=0x7f7f7f7f;
const double pi=acos(-1.0); 
const double eps=1e-10;

int dcmp(double x)
{
    if(fabs(x)<eps) return 0;
    else return x>0?1:-1;
}

double MAX(double a,double b)
{
    return a>b?a:b;
}

double rand01() { return rand() / (double)RAND_MAX; }
double randeps() { return (rand01()-0.5)*eps; }

struct Point
{
    double x;
    double y;
    Point(double x=0,double y=0):x(x),y(y) {}
};

typedef Point Vector;
typedef vector<Point> Polygon;
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,double p) { return Vector(a.x*p,a.y*p); }
Vector operator / (Vector a,double p) { return Vector(a.x/p,a.y/p); }
bool operator < (const Vector& a,const Vector& b) 
{ 
    return a.x<b.x || (a.x==b.x && a.y<b.y); 
}
bool operator ==(const Vector& a,const Vector& b) 
{ 
    return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0; 
}

struct Point3
{
    double x,y,z;
    Point3(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}
};

typedef Point3 Vector3;
Vector3 operator + (Vector3 a,Vector3 b) { return Vector3(a.x+b.x,a.y+b.y,a.z+b.z); }
Vector3 operator - (Vector3 a,Vector3 b) { return Vector3(a.x-b.x,a.y-b.y,a.z-b.z); }
Vector3 operator * (Vector3 a,double p) { return Vector3(a.x*p,a.y*p,a.z*p); }
Vector3 operator / (Vector3 a,double p) { return Vector3(a.x/p,a.y/p,a.z/p); }
bool operator < (const Vector3& a,const Vector3 &b) 
{ 
    return a.x<b.x || (a.x==b.x && a.y<b.y) || (a.x==b.x && a.y==b.y && a.z<b.z); 
}
bool operator ==(const Vector3& a,const Vector3& b) 
{ 
    return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0 && dcmp(a.z-b.z)==0; 
}

struct Line
{
    Point p;
    Vector v; 
    double ang;
    Line() {} 
    Line(Point p,Vector v):p(p),v(v) { ang=atan2(v.y,v.x); }
    Point point(double t) { return p+v*t; }
    bool operator < (const Line& rhs) const
    {
        return ang<rhs.ang;
    }
};

struct Plane
{
    /*点法式:n为法向量*/ 
    Point3 p;
    Vector3 n;
    double D;
    Plane() {}
    Plane(Point3 p,Vector3 n):p(p),n(n) { D=-(p.x*n.x+p.y*n.y+p.z*n.z); }
    bool isOnPlane(Point3 p0)
    {
        return dcmp(n.x*p0.x+n.y*p0.y+n.z*p0.z+D)==0; 
    }
};

struct Circle
{
    Point c;
    double r;
    Circle(Point c,double r):c(c),r(r) {}
    Point point(double rad) { return Point(c.x+cos(rad)*r,c.y+sin(rad)*r); }
}; 

struct Globe
{
    Point3 c;
    double r;
    double x,y,z;
    Globe(Point3 c,double r):c(c),r(r) {}
    Globe(double r):r(r) {}
    void GetCoord(double lat,double lng)
    {
        lat=lat/180*pi;
        lng=lng/180*pi;
        x=r*cos(lat)*cos(lng);
        y=r*cos(lat)*sin(lng);
        z=r*sin(lat);
    }
};

double Dot(Vector a,Vector b) { return a.x*b.x+a.y*b.y; }
double Dot(Vector3 a,Vector3 b) { return a.x*b.x+a.y*b.y+a.z*b.z; }

double Cross(Vector a,Vector b) { return a.x*b.y-b.x*a.y; }
Vector3 Cross(Vector3 a,Vector3 b) { return Vector3(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x); }

double Length(Vector a) { return sqrt(Dot(a,a)); }
double Length(Vector3 a) { return sqrt(Dot(a,a)); }

double Angle(Vector a) { return atan2(a.y,a.x); }

double Angle(Vector a,Vector b) { return acos(Dot(a,b)/Length(a)/Length(b)); }
double Angle(Vector3 a,Vector3 b) { return acos(Dot(a,b)/Length(a)/Length(b)); }

double Area2(Point a,Point b,Point c) { return Cross(b-a,c-a); }
double Area2(Point3 a,Point3 b,Point3 c) { return Length(Cross(b-a,c-a)); }
double Volume6(Point3 a,Point3 b,Point3 c,Point3 d) { return Dot(d-a,Cross(b-a,c-a)); }

bool OnLeft(Line L,Point p) { return Cross(L.v,p-L.p) > 0; }

Vector Normal(Vector a) { int l=Length(a);return Vector(-a.y/l,a.x/l); }

Vector Rotate(Vector a,double rad) 
{ 
    return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad)); 
}

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;
    /*double t=Cross(v,u)/Cross(v,w);
    return Q/=w*t;*/
} 

double DistanceToLine(Point P,Point a,Point b)
{
    Vector v1=b-a;
    Vector v2=P-a;
    return fabs(Cross(v1,v2)/Length(v1));
} 

double DistanceToLine(Point3 p,Point3 a,Point3 b)
{
    Vector3 v1=b-a;
    Vector3 v2=p-a;
    return Length(Cross(v1,v2))/Length(v1);
} 

double DistanceToSegment(Point P,Point a,Point b)
{
    Vector v1=b-a;
    Vector v2=P-a;
    Vector v3=P-b;
    if(a==b) return Length(v2);
    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));
}

double DistanceToSegment(Point3 p,Point3 a,Point3 b)
{
    Vector3 v1=b-a;
    Vector3 v2=p-a;
    Vector3 v3=p-b;
    if(a==b) return Length(v2);
    if(dcmp(Dot(v1,v2))<0) return Length(v2);
    else if(dcmp(Dot(v1,v3))>0) return Length(v3);
    else return Length(Cross(v1,v2))/Length(v1);
}

double DistanceToPlane(const Point3& p,const Point3& p0,const Vector3& n)
{
    /*n是单位法向量*/
    return fabs(Dot(p-p0,n));
} 

Point GetLineProjection(Point P,Point a,Point b)
{
    Vector v1=b-a;
    Vector v2=P-a;
    return a+v1*(Dot(v1,v2)/Dot(v1,v1));
} 

Point3 GetLineProjection(Point3 p,Point3 a,Point3 b)
{
    Vector3 v1=b-a;
    Vector3 v2=p-a;
    return a+v1*(Dot(v1,v2)/Dot(v1,v1));
}

Point3 GetPlaneProjection(const Point3& p,const Point3& p0,const Vector3& n)
{   
    return  p-n*Dot(p-p0,n);
} 

Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Vector3 n)
{
    Vector3 v=p2-p1;
    double t=Dot(n,p0-p1)/Dot(n,p2-p1);
    return p1+v*t;
}

bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2)
{
    double c1=Cross(a2-a1,b1-a1);
    double c2=Cross(a2-a1,b2-a1);
    double d1=Cross(b2-b1,a1-b1);
    double d2=Cross(b2-b1,a2-b1);
    return dcmp(c1)*dcmp(c2)<0 && dcmp(d1)*dcmp(d2)<0;
}

bool isOnSegment(Point P,Point a,Point b)
{
    return dcmp(Cross(a-P,b-P))==0 && dcmp(Dot(a-P,b-P))<0;
} 

bool PointOnTri(Point3 p,Point3 p0,Point3 p1,Point3 p2)
{
    Vector3 n=Cross(p2-p0,p1-p0);
    Plane Tri(p0,n);
    if(Tri.isOnPlane(p))
    {
        double area1=Area2(p,p0,p1);
        double area2=Area2(p,p1,p2);
        double area3=Area2(p,p2,p0);
        double area=Area2(p0,p1,p2);
        return dcmp(area1+area2+area3-area)==0;
    }
    return 0;
}

bool TriSegIntersection(Point3 p0,Point3 p1,Point3 p2,Point3 a,Point3 b,Point3& p)
{
    Vector3 n=Cross(p1-p0,p2-p1);
    if(dcmp(Dot(n,b-a))==0) return false;
    else
    {
        double t=Dot(n,p0-a)/Dot(n,b-a);
        if(dcmp(t)<0 || dcmp(t-1)>0) return false;
        p=a+(b-a)*t;
        return PointOnTri(p,p0,p1,p2);
    }
} 

bool TriTriIntersection(Point3* T1,Point3* T2)
{
    Point3 p;
    for(int i=0;i<3;++i)
    {
        if(TriSegIntersection(T1[0],T1[1],T1[2],T2[i],T2[(i+1)%3],p)) return true;
        if(TriSegIntersection(T2[0],T2[1],T2[2],T1[i],T1[(i+1)%3],p)) return true;
    }
    return false;
}

bool LineDistance3D(Point3 p1,Vector3 u,Point3 p2,Vector3 v,double& s)
{
    double b=Dot(u,u)*Dot(v,v)-Dot(u,v)*Dot(u,v);
    if(dcmp(b)==0) return false;                            //平行或重合
    double a=Dot(u,v)*Dot(v,p1-p2)-Dot(v,v)*Dot(u,p1-p2);
    s=a/b;
    return false; 
} 

double DistanceSegToSeg(Point3 a,Point3 b,Point3 c,Point3 d)
{
    double s,t,dist[4];
    double ans=Max;
    if(LineDistance3D(a,b-a,c,d-c,s))
        if(LineDistance3D(c,d-c,a,b-a,t))
        {
            Point3 P=a+(b-a)*s;
            Point3 Q=c+(d-c)*t;
            return Length(Q-P);
        }
    dist[0]=DistanceToSegment(a,c,d);
    dist[1]=DistanceToSegment(b,c,d);
    dist[2]=DistanceToSegment(c,a,b);
    dist[3]=DistanceToSegment(d,a,b);
    for(int i=0;i<4;++i)
      ans=(ans<dist[i]?ans:dist[i]);
    return ans;
} 

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,double& t1,double& t2,vector<Point>& sol)
{
    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(delta < 0) return 0;
    else if(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;
    }
}

int GetCircleCircleIntersection(Circle c1,Circle c2,vector<Point>& sol)
{
    double d=Length(c1.c-c2.c);
    if(dcmp(d==0))
    {
      if(dcmp(c1.r-c2.r)==0) return -1;
      else 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);
    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=c2.point(a-da);
    sol.push_back(p1);
    if(p1==p2) return 1;
    sol.push_back(p2);
    return 2;
} 

int GetTangents(Point p,Circle c,vector<Vector>& sol)
{
    Vector u=c.c-p;
    double d=Length(u);
    if(d<c.r) return 0;
    else if(dcmp(d-c.r)==0)
    {
        sol.push_back(Vector(-u.y,u.x));
        return 1;
    }
    else
    {
        double ang=asin(c.r/d);
        sol.push_back(Rotate(u,-ang));
        sol.push_back(Rotate(u,+ang));
        return 2;
    }
}

int GetTangents(Circle c1,Circle c2,Point* a,Point* b)
{
    int cnt=0;
    if(c1.r<c2.r) { swap(c1,c2),swap(a,b); }
    double d=Length(c2.c-c1.c);
    double rdif=c1.r-c2.r;
    double rsum=c1.r+c2.r;
    if(dcmp(d-rdif)<0) return 0;
    
    double base=Angle(c2.c-c1.c);
    if(d==0 && c1.r==c2.r) return -1;
    
    if(dcmp(c1.r-c2.r-d)==0)
    {
        a[cnt]=c1.point(base);
        b[cnt]=a[cnt];
        cnt++;
        return 1;
    } 
    
    double ang=acos((c1.r-c2.r)/d);
    a[cnt]=c1.point(base+ang);
    b[cnt++]=c2.point(base+ang);
    a[cnt]=c1.point(base-ang);
    b[cnt++]=c2.point(base-ang);
    
    if(dcmp(d-rsum)==0)
    {
        a[cnt]=c1.point(base);
        b[cnt]=a[cnt];
        cnt++;
    }
    else if(dcmp(d-rsum)>0)
    {
        double ang=acos((c1.r+c2.r)/d);
        a[cnt]=c1.point(base+ang);
        b[cnt++]=c2.point(pi+base+ang);
        a[cnt]=c1.point(base-ang);
        b[cnt++]=c2.point(pi+base-ang);
    }
    return cnt;
}

int isPointInPolygon(Point p,Polygon poly)
{
    int wn=0;
    int n=poly.size();
    for(int i=0;i<n;++i)
    {
        if(isOnSegment(p,poly[i],poly[(i+1)%n])) return -1;     //边界上 
        int k=dcmp(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));
        int d1=dcmp(poly[i].y-p.y);
        int d2=dcmp(poly[(i+1)%n].y-p.y);
        if(k>0 && d1<=0 && d2>0) wn++;                          //逆时针穿过 
        if(k<0 && d2<=0 && d1>0) wn--;                          //顺时针穿过 
    }
    if(wn != 0) return 1;                                       //内部 
    return 0;                                                   //外部 
}

int  ConvexHull(Point* p,int n,Point* ch)
{
    sort(p,p+n);
    n=unique(p,p+n)-p;
    int cnt=0;
    for(int i=0;i<n;++i)
    {
        while(cnt>1 && Cross(ch[cnt-1]-ch[cnt-2],p[i]-ch[cnt-2]) <= 0) cnt--;
        ch[cnt++]=p[i];
    }
    int k=cnt;
    for(int i=n-2;i>=0;--i)
    {
        while(cnt>k && Cross(ch[cnt-1]-ch[cnt-2],p[i]-ch[cnt-2]) <= 0) cnt--;
        ch[cnt++]=p[i];
    }
    if(n>1) cnt--;
    return cnt;
} 

double diameter2(Point* p,int n,Point* ch)
{
    n=ConvexHull(p,n,ch);
    if(n==1) return 0;
    if(n==2) return Length(ch[1]-ch[0]);
    ch[n]=ch[0];
    double ans=0;
    for(int u=0,v=1;u<n;++u)
    {
        for( ; ; )
        {
            double diff=Cross(ch[u+1]-ch[u],ch[v+1]-ch[v]);
            if(diff<=0)
            {
                ans=MAX(ans,Length(ch[u]-ch[v]));
                if(diff==0) ans=MAX(ans,Length(ch[u]-ch[v+1]));
                break;
            }
            v=(v+1)%n;
        }
    }
    return ans;
} 

Polygon CutPolygon(Polygon poly,Point a,Point 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(isOnSegment(ip,c,d))
              newpoly.push_back(ip);
        }
    }
    return newpoly;
}

int HalfplaneIntersection(Line* l,int n,Point* poly)                
{
    sort(l,l+n);                                                //将直线按极角排序 
    int first=0,last=0;
    Point* p=new Point[n];                                      //p[i]为q[i]与q[i-1]的交点 
    Line* q=new Line[n];
    q[first]=l[0];
    for(int i=1;i<n;++i)
    {
        while(first < last && !OnLeft(l[i],p[last-1])) last--;
        while(first < last && !OnLeft(l[i],p[first])) first++;
        q[++last]=l[i];
        if(dcmp(Cross(q[last].v,q[last-1].v)) == 0)
        {
            last--;
            if(OnLeft(q[last],l[i].p)) 
                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 && !OnLeft(q[first],p[last-1])) last--;
    if(last-first <= 1) return 0;
    p[last]=GetLineIntersection(q[first].p,q[first].v,q[last].p,q[last].v);
    int m=0;
    for(int i=first;i<=last;++i) poly[m++]=p[i];
    return m;
}

double SphericalDistance(double r,double lat1,double lng1,double lat2,double lng2)
{
    Globe g(r);
    g.GetCoord(lat1,lng1);
    Point3 x1=Point3(g.x,g.y,g.z);
    g.GetCoord(lat2,lng2);
    Point3 x2=Point3(g.x,g.y,g.z);
    return 2*g.r*asin(Length(x2-x1)/(2*g.r));
}

struct Face
{
    int v[3];
    Vector3 Normal(Point3* P) const
    {
        return Cross(P[v[1]]-P[v[0]],P[v[2]]-P[v[0]]);
    }
    bool CanSee(Point3* P,int i)
    {
        return Dot(P[i]-P[v[0]],Normal(P))>0?true:false; 
    }
}; 

Point3 add_noise(Point3 p)
{
    return Point3(p.x+randeps(),p.y+randeps(),p.z+randeps());
}

vector<Face> ConvexHall3D(Point3* P,int n)
{
    bool vis[max_point][max_point];
    memset(vis,0,sizeof(vis));
    vector<Face> cur;
    cur.push_back((Face){{0,1,2}});
    cur.push_back((Face){{2,1,0}});
    for(int i=3;i<n;++i)
    {
        vector<Face> next;
        for(unsigned int j=0;j<cur.size();++j)
        {
            Face& f=cur[j];
            bool res=f.CanSee(P,i);
            if(!res) next.push_back(f);
            for(int k=0;k<3;++k)
                vis[f.v[k]][f.v[(k+1)%3]]=res;
        }
        for(unsigned int j=0;j<cur.size();++j)
            for(int k=0;k<3;++k)
            {
                int a=cur[j].v[k];
                int b=cur[j].v[(k+1)%3];
                if(vis[a][b] != vis[b][a] && vis[a][b])
                    next.push_back((Face){{a,b,i}});
            }
        cur=next; 
    }
    return cur;
}

int main()
{
    /*注释:
     1.Dot() 点积,Cross() 叉积; 
     2.Length() 向量的模;
     3.Angle() 两向量夹角; 
     4.Area2() 计算三点所围三角形面积的二倍; 
     5.Valume6() 计算空间四点围成空间四边形体积的六倍; 
     6.OnLeft() 判断点p是否在直线l的左侧;
     7.Rotate() 向量的逆时针旋转 [rad是弧度制]; 
     8.Normal() 用来计算向量的单位法线,将长度归一化; 
     9.GetLineIntersection() 参数方程下的直线交点公式;
    10.DistanceToLine() 点到直线的距离; 
    11.DistanceToSegment() 点到线段的距离; 
    12.DistanceToPlane() 点到平面(点法式)的距离; 
    13.GetLineProjection() 点在直线上的投影; 
    14.GetPlaneIntersection() 点在平面上的投影; 
    15.LinePlaneIntersection() 直线和平面的交点; 
    16.SegmentProperIntersection() 线段规范相交判定; 
    17.isOnSegment() 判定点是否在线段上; 
    18.PointOnTri() 判定点是否在空间三点围成的三角形上;
    19.TriSegIntersection() 判定三角形是否与空间内一线段相交(不考虑在同一平面的情况); 
    20.TriTriIntersection() 判定两个空间三角形是否相交; 
    21.LineDistance3D() 求异面直线p1+su与p2+tv的公垂线对应的s,若两线段平行或重合返回false;
    22.DistanceSegToSeg() 空间内线段到线段的距离; 
    23.PolygonArea() 计算多边形面积; 
    24.GetLineCircleIntersection() 计算直线与圆的交点,返回交点个数;
    25.GetCircleCircleIntersection() 计算圆与圆的交点,返回交点个数; 
    26.GetTangents() 计算点p与圆c的切点,返回切点个数; 
    27.GetTangents() 计算两圆的公切线切点,返回切线个数; 
    28.isPointInPolygon() 判断点是否在多边形内; 
    29.ConvexHull() 求点集的凸包,返回凸包顶点数;
    30.diameter2() 求凸包的直径; 
    31.CutPolygon() 用有向直线a->b切割多边形poly;
    32.HalfplaneIntersection 计算半平面交; 
    33.SphericalDistance() 计算球面距; 
    34.add_noise() 添加扰动,为防止三维凸包计算中多个点在同一面内的情况添加扰动; 
    35.ConvexHull3D() 增量法计算(经过扰动后且去重的)点集的三维凸包,返回凸包上的点的集合; 
    */ 
}

//用复数运算模拟向量运算--效率不高 
//real()实部  imag()虚部  conj()共轭复数 .
//cos(rad)+sin(rad)*i=sqrt(cos(rad)^2+sin(rad)^2)*e^(i*theta)
#include<complex>
using namespace std;
typedef complex<double> Point;
typedef Point Vector;
double Dot(Vector a,Vector b) { return real(conj(a)*b); }
double Cross(Vector a,Vector b) { return imag(conj(a)*b); }
Vector Rotate(Vector a,double rad) { return a*exp(Point(0,rad)); } 
Vector Rotate(Vector a,double rad) { return a*Point(cos(rad),sin(rad)); } 
int main() {}

扫描二维码,在手机上阅读!
Responses