计算几何模板

二维几何

坐标运算

点和向量共用的坐标Coord运算，使用stl中的容器complex实现（参考地址）。
复数相乘的几何意义为长度相乘，极角相加。

const double EPS=1e-9,PI=acos(-1);typedef complex<double> Coord;Coord(double x=0,double y=0);//构造函数double real(Coord a);//a的实部（复平面的横坐标）,也可写作a.real()double imag(Coord a);//a的虚部（复平面的纵坐标）,也可写作a.imag()double abs(Coord a);//向量a的模长，或是点a到原点的距离double norm(Coord a);//abs的平方，比abs快，但是要注意浮点数精度溢出double arg(Coord a);//a的幅角，与atan2(a.real(),a.imag())等价Coord polar(double r,double t);//极坐标生成方式，r为幅值，t为幅角//运算符重载,<< >>运算符输出括号形式的坐标+、-、*、/（以及对应的赋值运算，但是赋值运算不能写在表达式中，详见参考地址）、<<、>>int sgn(double x)//浮点数判断符号 {    return fabs(x)<EPS?0:           x<0?-1:1;}bool operator!=(Coord a,Coord b)//不等运算符，涉及到浮点数比较要重写{    return sgn(a.real()-b.real())||sgn(a.imag()-b.imag());}bool operator==(Coord a,Coord b){    return !(a!=b);}

点和向量

typedef Coord Point;typedef Coord Vec;double Dot(Vec A,Vec B)//点乘{    return A.real()*B.real()+A.imag()*B.imag();}double Cross(Vec A,Vec B)//叉乘 {    return A.real()*B.imag()-A.imag()*B.real();}Vec Rotate(Vec A,double rad)//向量逆时针旋转rad度(弧度){    return A*polar(1.0,rad);}double Angle(Vec A,Vec B)//向量夹角{    return acos(Dot(A,B)/abs(A)/abs(B));}double Area2(Point A,Point B,Point C)//向量AB、向量AC三角形有向面积的二倍{    return Cross(B-A,C-A);}Vec Normal(Vec A)//计算非零向量A的单位法向量，左转90°，把长度归一{    return Vec(-A.imag(),A.real())/abs(A);}

直线和线段

用直线上的一点p和方向向量v表示一条经过p的直线，直线上的所有点q满足q=p+t*v，其中t是参数；当限制t≥0时，该参数方程表示射线；限制0≤t≤1时，该参数方程表示线段。
此外，如果已知线段端点a1和a2，可以通过Line(a1,a2-a1)来得到对应的参数形式。

struct Line{    Point p;    Vec v;    Line() {}    Line(Point p,Vec v):p(p),v(v) {}    Point point(double t)//线上一点     {        return p+t*v;    }};bool Onleft(Point P,Line L)//p是否在有向直线L左侧，不含线上{    return Cross(L.v,P-L.p)>0;}double DistanceToLine(Point P,Line L)//点到直线距离（有向）{    return Cross(L.v,P-L.p)/abs(L.v);}double DistanceToSegmentS(Point P,Point A,Point B)//点到线段的距离（无向）{    if(A==B)return abs(P-A);    Vec v1=B-A,v2=P-A,v3=P-B;    if(sgn(Dot(v1,v2))<0)return abs(v2);    if(sgn(Dot(v1,v3))>0)return abs(v3);    return fabs(DistanceToLine(P,Line(A,B-A)));}Point GetLineProjection(Point P,Line L)//点在直线上的投影{    return L.point(Dot(L.v,P-L.p)/norm(L.v));}Point GetLineIntersection(Line L1,Line L2)//直线交点,须确保两直线相交。{    return L1.point(Cross(L2.v,L1.p-L2.p)/Cross(L1.v,L2.v));}bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2)//线段相交判定，交点不在一条线段的端点{    double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1),           c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);    return sgn(c1)*sgn(c2)<0&&sgn(c3)*sgn(c4)<0;}bool OnSegment(Point P,Point a1,Point a2)//判断点是否在线段上，不包含端点{    return sgn(Dot(a1-P,a2-P))<0&&!sgn(Cross(a1-P,a2-P));}

多边形和凸包

Morley定理：三角形每个内角的三等分线相交成等边三角形。
欧拉定理：平面图的点数V、边数E和面数F满足V+F-E=2。

typedef vector<Point> Polygon;//表示一个点集、多边形或凸包int isPointInPolygon(Point p,Polygon poly)//点在多边形内的判定，转角法，正值为内部，0为外部，-1在边界上{    int ans=0;    for(int i=0,k,d1,d2,n=poly.size(); i!=n; ++i)    {        if(OnSegment(p,poly[i],poly[(i+1)%n]))return -1;//在边界上        k=sgn(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));        d1=sgn(poly[i].imag()-p.imag());        d2=sgn(poly[(i+1)%n].imag()-p.imag());        if(k>0&&d1<=0&&d2>0)++ans;        if(k<0&&d2<=0&&d1>0)++ans;    }    return ans;}double PolygonArea(Polygon p)//计算多边形的有向面积，凸多边形即为面积{    double area=0;    for(int i=1,n=p.size()-1; i<n; ++i)        area+=Area2(p[0],p[i],p[i+1]);    return area/2;}struct PointCmp//先比横坐标再比纵坐标的比较类{    bool operator()(Point a,Point b)    {        return sgn(a.real()-b.real())?               a.real()<b.real():               a.imag()<b.imag();    }};int ConvexHull(Polygon p,Polygon &sol)//获得凸包；不希望凸包的边上有输入点，把两个<=改成<{    sort(p.begin(),p.end(),PointCmp());//先比横坐标再比纵坐标     for(int i=0; i!=p.size(); ++i)    {        while(sol.size()>1&&            Area2(sol[sol.size()-2],sol[sol.size()-1],p[i])<=0)            sol.pop_back();        sol.push_back(p[i]);    }    for(int i=sol.size()-2,k=sol.size(); i>=0; --i)    {        while(sol.size()>k&&                Area2(sol[sol.size()-2],sol[sol.size()-1],p[i])<=0)            sol.pop_back();        sol.push_back(p[i]);    }    if(p.size()>1)sol.pop_back();    return sol.size();}Polygon CutPolygon(Polygon poly,Point A,Point B)//用有向直线A->B切割多边形poly， 返回“左侧”。 如果退化，可能会返回一个单点或者线段，复杂度O(n2){    Polygon newpoly;    for(int i=0,n=poly.size(); i!=n; ++i)    {        Point C=poly[i],D=poly[(i+1)%n];        if(sgn(Cross(B-A,C-A))>=0)newpoly.push_back(C);        if(!sgn(Cross(B-A, C-D)))        {            Point ip=GetLineIntersection(Line(A,B-A),Line(C,D-C));            if(OnSegment(ip,C,D))newpoly.push_back(ip);        }    }    return newpoly;}struct LineCmp//比较极角的比较类{    bool operator()(Line a,Line b)    {        return arg(a.v)<arg(b.v);    }};Polygon HalfplaneIntersection(vector<Line> L)//半平面交{    sort(L.begin(),L.end(),LineCmp());//按极角排序    int first,last;//双端队列的第一个元素和最后一个元素    vector<Point> p(L.size(),Point()); //p[i]为q[i]和q[i+1]的交点    vector<Line> q(L.size(),Line());//双端队列    q[first=last=0]=L[0]; //队列初始化为只有一个半平面L[0]    for(int i=0,n=L.size(); 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(!sgn(Cross(q[last].v,q[last-1].v)))        {            --last;            if(Onleft(L[i].p,q[last]))                q[last]=L[i];        }        if(first<last)            p[last-1]=GetLineIntersection(q[last-1], q[last]);    }    while(first<last&&!Onleft(p[last-1],q[first]))        --last;//删除无用平面    if(last-first<=1)return Polygon();//空集    p[last]=GetLineIntersection(q[last],q[first]);    return Polygon(p.begin()+first,p.begin()+last+1);//从deque复制到输出中}

圆

struct Circle{    Point c;    double r;    Circle(Point c,double r):c(c),r(r) {}    Point point(double t)//通过圆心角t确定圆上坐标    {        return c+polar(r,t);    }};Line getTangent(Point C,Point P)//圆心C，圆上一点P处切线{    return Line(P,Normal(C-P));}int getLineCircleIntersecion(Line L,Circle C,Polygon &sol)//直线和圆的交点，返回交点数量，sol保存解{    double a=L.v.real(),b=L.p.real()-C.c.real(),           c=L.v.imag(),d=L.p.imag()-C.c.imag();    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(sgn(delta)<0)return 0;//相离    if(sgn(delta)==0)    {        double t=-f/2/e;        sol.push_back(L.point(t));        return 1;    }    double 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,Polygon &sol)//两圆相交，返回值为-1时表示重合，其余同上{    double d=abs(C1.c - C2.c);    if(sgn(d)==0&&sgn(C1.r-C2.r)==0)return -1;//同心    if(C1.r+C2.r<d||d<fabs(C1.r-C2.r))return 0;//相离或内含    double a=arg(C2.c - C1.c),//向量C1C2的极角           da=acos((d*d+C1.r*C1.r-C2.r*C2.r)/(2*d*C1.r));//C1C2到C1P1的角    if(sgn(da)==0)//相切    {        sol.push_back(C1.point(a));        return 1;    }    sol.push_back(C1.point(a-da));    sol.push_back(C1.point(a+da));    return 2;}int getTangents(Point p,Circle C,Polygon &sol)//过点p做圆C的切线，返回切线个数{    Vec u=p-C.c;    double d=abs(u);    if(d<C.r)return 0;//点在圆内    else if(sgn(d-C.r)==0)//点在圆上    {        sol.push_back(p);        return 1;    }    double base=arg(u),ang=acos(C.r/d);    sol.push_back(C.point(base+ang));    sol.push_back(C.point(base-ang));    return 2;}int getTangents(Circle A, Circle B,Polygon &a,Polygon &b)//两圆的公切线，返回切线的个数，-1表示有无数条公切线，a[i]、b[i]表示第i条切线在圆A，圆B上的切点{    int d=abs(A.c-B.c),rdiff =fabs(A.r-B.r),rsum =A.r+B.r;    if(sgn(d)==0&&sgn(rdiff)==0)return -1;//重合     if(d<rdiff)return 0;//内含    double base=arg(B.c-A.c);    if(sgn(d-rdiff)==0)//内切    {        a.push_back(A.point(base));        b.push_back(B.point(base));        return 1;    }    //有外共切线,两条夹住两圆的切线     int cnt=2;    double ang=acos(rdiff/d);    a.push_back(A.point(base+ang));    a.push_back(A.point(base-ang));    b.push_back(B.point(base+ang));    b.push_back(B.point(base+ang));    if(sgn(d-rsum)==0)//外切    {        ++cnt;        a.push_back(A.point(base));        b.push_back(B.point(PI+base));    }    else if(d>rsum)//外离，再加两条     {        cnt+=2;        double ang = acos(rsum/d);        a.push_back(A.point(base+ang));        a.push_back(A.point(base-ang));        b.push_back(B.point(PI+base+ang));        b.push_back(B.point(PI+base-ang));    }    return cnt;}

三维几何

坐标运算

增加了一个add_noise增加扰动函数，可以让坐标发生微小偏移。

struct Coord3{    double x,y,z;    Coord3(double x=0,double y=0,double z=0):x(x),y(y),z(z) {}    Coord3& operator+=(const Coord3 &a){return x+=a.x,y+=a.y,z+=a.z,*this;}    Coord3& operator-=(const Coord3 &a){return x-=a.x,y-=a.y,z-=a.z,*this;}    Coord3& operator*=(const double &d){return x*=d,y*=d,z*=d,*this;}    Coord3& operator/=(const double &d){return x/=d,y/=d,z/=d,*this;}    double randeps(){return (rand()/(double)RAND_MAX-0.5)*EPS;}    void add_noise(){*this+=Coord3(randeps(),randeps(),randeps());}};bool operator!=(const Coord3 &a,const Coord3 &b)//不等运算符，涉及到浮点数比较要重写{    return sgn(a.x-b.x)||sgn(a.y-b.y)||sgn(a.z-b.z);}bool operator==(const Coord3 &a,const Coord3 &b){    return !(a!=b);}Coord3 operator+(Coord3 a,const Coord3 &b){    return a+=b;}Coord3 operator-(Coord3 a,const Coord3 &b){    return a-=b;}Coord3 operator*(Coord3 a,const double &d){    return a*=d;}Coord3 operator/(Coord3 a,const double &d){    return a/=d;}

点和向量

typedef Coord3 Point3;typedef Coord3 Vec3;double Dot(const Vec3& A,const Vec3& B){    return A.x*B.x+A.y*B.y+A.z*B.z;}Vec3 Cross(const Vec3& A, const Vec3& B){    return Vec3(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 abs(const Vec3& A){    return sqrt(Dot(A,A));}double Angle(const Vec3& A,const Vec3& B){    return acos(Dot(A,B)/abs(A)/abs(B));}double Area2(Point3 A,Point3 B, Point3 C){    return abs(Cross(B-A,C-A));}double Volume6(Point3 A, Point3 B, Point3 C, Point3 D)//四面体体积{    return Dot(D-A,Cross(B-A, C-A));}Point3 Centroid(Point3 A, Point3 B, Point3 C, Point3 D)//四面体的重心{    return (A+B+C+D)/4.0;}

线和平面

double DistanceToPlane(Point3 p, Point3 p0, const Vec3& n)// 点p到平面p0-n的有向距离。n必须为单位向量{    return Dot(p-p0,n);}Point3 GetPlaneProjection(Point3 p, Point3 p0, const Vec3& n)// 点p在平面p0-n上的投影。n必须为单位向量{    return p-n*Dot(p-p0,n);}Point3 LinePlaneIntersection(Point3 p1, Point3 p2, Point3 p0, Vec3 n)//直线p1-p2 与平面p0-n的交点，假设交点唯一存在{    Vec3 v = p2-p1;    double t = Dot(n, p0-p1) / Dot(n, p2-p1);//分母为0，直线与平面平行或在平面上    return p1 + v*t; //如果是线段 判断t是否在0~1之间}double DistanceToLine(Point3 P,Point3 A,Point3 B)// 点P到直线AB的距离{    Vec3 v1=B-A,v2=P-A;    return abs(Cross(v1, v2))/abs(v1);}double DistanceToSeg(Point3 P, Point3 A, Point3 B)//点到线段的距离{    if(A==B)return abs(P-A);    Vec3 v1=B-A,v2=P-A,v3=P-B;    if(sgn(Dot(v1,v2))<0)return abs(v2);    if(sgn(Dot(v1,v3))>0)return abs(v3);    return fabs(DistanceToLine(P,A,B));}bool LineDistance3D(Point3 p1, Vec3 u, Point3 p2, Vec3 v, double& s)//求异面直线 p1+s*u与p2+t*v的公垂线对应的s，如果平行|重合，返回0{    double b = Dot(u, u) * Dot(v, v) - Dot(u, v) * Dot(u, v);    if(sgn(b) == 0) return 0;    double a = Dot(u, v) * Dot(v, p1-p2) - Dot(v, v) * Dot(u, p1-p2);    s = a/b;    return 1;}bool SameSide(Point3 p1, Point3 p2, Point3 a, Point3 b)// p1和p2是否在线段a-b的同侧{    return sgn(Dot(Cross(b-a, p1-a), Cross(b-a, p2-a))) >= 0;}bool PointInTri(Point3 PP, Point3 *P)// 点P在三角形P0, P1, P2中{    return SameSide(PP,P[0],P[1],P[2])&&           SameSide(PP,P[1],P[0],P[2])&&           SameSide(PP,P[2],P[0],P[1]);}bool TriSegIntersection(Point3* P,Point3 A,Point3 B,Point3& PP)// 三角形P0P1P2是否和线段AB相交{    Vec3 n = Cross(P[1]-P[0], P[2]-P[0]);    if(sgn(Dot(n, B-A)) == 0) return false; // 线段A-B和平面P0P1P2平行或共面    double t = Dot(n, P[0]-A) / Dot(n, B-A);// 平面A和直线P1-P2有惟一交点    if(sgn(t) < 0 || sgn(t-1) > 0) return false; // 不在线段AB上    PP = A + (B-A)*t; // 交点    return PointInTri(PP, P);}bool TriTriIntersection(Point3* T1, Point3* T2)//空间两三角形是否相交{    Point3 P;    for(int i = 0; i < 3; i++)    {        if(TriSegIntersection(T1, T2[i], T2[(i+1)%3], P)||                TriSegIntersection(T2, T1[i], T1[(i+1)%3], P)) return true;    }    return false;}double SegSegDistance(Point3 a1, Point3 b1, Point3 a2, Point3 b2,Point3 &ans1,Point3 &ans2)//空间两直线上最近点对 返回最近距离 点对保存在ans1 ans2中{    Vec3 v1 = (a1-b1), v2 = (a2-b2);    Vec3 N = Cross(v1, v2);    Vec3 ab = (a1-a2);    double ans = Dot(N, ab) / abs(N);    Point3 p1 = a1, p2 = a2;    Vec3 d1 = b1-a1, d2 = b2-a2;    double t1, t2;    t1 = Dot((Cross(p2-p1, d2)), Cross(d1, d2));    t2 = Dot((Cross(p2-p1, d1)), Cross(d1, d2));    double dd = abs((Cross(d1, d2)));    t1 /= dd*dd;    t2 /= dd*dd;    ans1 = (a1 + (b1-a1)*t1);    ans2 = (a2 + (b2-a2)*t2);    return fabs(ans);}bool InsideWithMinDistance(Point3 PP, Point3 *P, double mindist)// 判断PP是否在三角形P中，并且到三条边的距离都至少为mindist。保证P, A, B, C共面{    return PointInTri(PP,P)&&           DistanceToLine(PP,P[0],P[1])>=mindist||           DistanceToLine(PP,P[1],P[2])>=mindist||           DistanceToLine(PP,P[2],P[0])>=mindist;}

空间多边形和凸包

struct ConvexPolyhedron//凸包问题{    struct Face    {        int v[3];        Face(int a, int b, int c)        {            v[0] = a,v[1] = b,v[2] = c;        }        Vec3 Normal(const vector<Point3>& P) const        {            return Cross(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);        }        bool CanSee(const vector<Point3>& P, int i) const// f是否能看见P[i]        {            return Dot(P[i]-P[v[0]], Normal(P)) > 0;        }    };    vector<Point3> P;    vector<Face> faces;    void read(vector<Point3> P2)    {        P=P2;        for(int i = 0; i < P2.size(); i++)            P2[i].add_noise();        faces = CH3D(P2);    }    vector<Face> CH3D(const vector<Point3>& P)//增量法求三维凸包，P已经过扰动，没有考虑各种特殊情况（如四点共面）    {        int n = P.size();        vector<vector<int> > vis(n);        for(int i = 0; i < n; i++) vis[i].resize(n);        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(int j = 0; j < cur.size(); j++)// 计算每条边的“左面”的可见性            {                Face& f = cur[j];                int 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(int j = 0; j < cur.size(); j++)                for(int k = 0; k < 3; k++)                {                    int a = cur[j].v[k], b = cur[j].v[(k+1)%3];                    if(vis[a][b] != vis[b][a] && vis[a][b]) // (a,b)是分界线，左边对P[i]可见                        next.push_back(Face(a, b, i));                }            cur = next;        }        return cur;    }    Point3 centroid()//三维凸包重心    {        Point3 C = P[0],tot(0,0,0);        double totv = 0;        for(int i = 0; i < faces.size(); i++)        {            Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];            double v = -Volume6(p1, p2, p3, C);            totv += v;            tot = tot + Centroid(p1, p2, p3, C)*v;        }        return tot / totv;    }    double mindist(Point3 C)//凸包重心到表面最近距离    {        double ans = 1e30;//INF        for(int i = 0; i < faces.size(); i++)        {            Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];            ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));        }        return ans;    }};