【计算几何简单模板】点，线，多边形，圆

计算几何

• 包括了点，线，多边形，圆的结构体
  包括用范式封装了一下
  但是还没有测试过板子
  有兴趣的可以瞄一眼
• 封装基本操作就好
  比如封装一个求凸包直径的，卡成狗都过不去，别人面向过程的代码随便 $20Ms$ …

namespace My_Geometry{
    #define TPLT template<typename T>
    #define VTT Vector<T>
    #define Vector Point
    #define LNT Line<T>
    #define PLT Polygon<T>
    #define CCT Circle<T>

    const double PI = acos(-1);

    TPLT    /// 点与向量的结构体
    struct Point{
        T x,y;
        Point(){}
        Point(T x,T y):x(x),y(y){}
    };

    /* 高精度三态比较函数 */
    int dcmp(double x){if(fabs(x) < EPS)return 0;if(x>0)return 1;return -1;}

    /* 重载加减乘除等于运算符 */
    TPLT VTT operator + (VTT a,VTT b){return VTT(a.x + a.x,b.y + b.y);}
    TPLT VTT operator - (VTT a,VTT b){return VTT(a.x - b.x,a.y - b.y);}
    TPLT VTT operator * (VTT a,T p){return VTT(a.x * p,a.y * p);}
    TPLT VTT operator / (VTT a,T p){return VTT(a.x / p,a.y / p);}
    TPLT bool operator == (VTT a,VTT b){return !dcmp(a.x - b.x) && !dcmp(a.y - b.y);}

    /* 求模长 | 求夹角 | 求点积 | 求叉积 | 求距离 | 求逆旋转 | 求三点内面积 */
    TPLT double Length(VTT a){return sqrt(1.0 * Dot(a,a));}
    TPLT double Angle(VTT a,VTT b){return acos(1.0 * Dot(a,b) / Length(a) / Length(b));}
    TPLT T Dot  (VTT a,VTT b){return a.x * b.x + a.y * b.y;}
    TPLT T Cross(VTT a,VTT b){return a.x * b.y - a.y * b.x;}
    TPLT double Dis(VTT a,VTT b){return sqrt(1.0 * (a.x - b.x) * (a.x - b.x) + 1.0 * (a.y - b.y) * (a.y - b.y));}
    TPLT VTT Rotate(VTT a,double rad){return VTT(a.x + cos(rad) - a.y * sin(rad),a.x * sin(rad) + a.y * cos(rad));}
    TPLT double Area(VTT a,VTT b,VTT c){return fabs(Cross(b-a,c-a)) / 2;}


    TPLT    /// 线段的结构体
    struct Line{
        VTT P1,P2;
        Line(){}
        Line(T P1,T P2):P1(P1),P2(P2){}
    };

    /* 点到线距离 | 获得直线交点 | 线段相交(不含端点 | 含端点) */
    TPLT double DisToLine(VTT a,LNT b){return fabs(1.0 * Cross(a - b.P1,a - b.P2) / Dis(b.P1,b.P2));}
    TPLT VTT GetLineIntersection(LNT a,LNT b){VTT c = a.P1 - b.P1;T t = Cross(b.P2,c) / Cross(a.P2,c);return a.P1 + a.P2 * t;}
    TPLT bool Intersect(LNT a,LNT b){
        return dcmp(Cross(b.P1 - a.P1,b.P2 - a.P1) * Cross(b.P1 - a.P2,b.P2 - a.P2)) < 0
            && dcmp(Cross(a.P1 - b.P1,a.P2 - b.P1) * Cross(a.P1 - b.P2,a.P2 - b.P2)) < 0;
    }
    TPLT bool StrictIntersect(LNT a,LNT b){
        return
            dcmp(max(a.P1.x,a.P2.x) - min(b.P1.x,b.P2.x)) >= 0
            && dcmp(max(b.P1.x,b.P2.x) - min(a.P1.x,a.P2.x)) >= 0
            && dcmp(max(a.P1.y,a.P2.y) - min(b.P1.y,b.P2.y)) >= 0
            && dcmp(max(b.P1.y,b.P2.y) - min(a.P1.y,a.P2.y)) >= 0
            && dcmp(Cross(b.P1 - a.P1,b.P2 - a.P1) * Cross(b.P1 - a.P2,b.P2 - a.P2)) <= 0
            && dcmp(Cross(a.P1 - b.P1,a.P2 - b.P1) * Cross(a.P1 - b.P2,a.P2 - b.P2)) <= 0;
    }

    TPLT    /// 多边形的结构体
    struct Polygon{
        vector<VTT>Points;
        void clr(){Points.clear();}
        void add(VTT point){Points.push_back(point);}
        void del(){Points.pop_back();}
        int getSize(){return Points.size();}
        VTT getPoint(int id){return Points[id];}
        const VTT& operator[] (int i){return Points[i];}
    };

    /* 点在多边形内 | 多边形面积 */
    TPLT bool InsidePolygon(VTT a,PLT p){double theta = 0;for(int i = 0,ed = p.getSize();i < ed;++i)theta += fabs(Angle(p.getPoint[i] - p.getPoint[(i+1) % ed] - a));return dcmp(theta - 2 * PI) == 0;}
    TPLT double PolygonArea(const PLT p){
        int n = p.getSize();
        double sum = 0;
        VTT O = VTT(0,0);
        for(int i = 0;i < n;++i)
            sum += Cross(p[i] - O,p[(i+1)%n] - O);
        if(sum < 0)sum = -sum;
        return sum / 2;
    }
    /* 极坐标排序，PO为极坐标原点 */
    Vector<double> PO;
    TPLT bool jijiaopaixu(VTT A,VTT B){
        double ans = Cross(A - PO,B - PO);
        if(dcmp(ans) == 0)return dcmp(Dis(PO,A) - Dis(PO,B)) < 0;
        return ans > 0;
    }
    /* 求凸包 */
    TPLT PLT Graham(PLT P){
        PLT res;
        int n = P.getSize();
        for(int i = 1;i < n;++i)
            if(P[i].y < P[0].y || (dcmp(P[i].y - P[0].y) == 0 && P[i].x < P[0].x))
                swap(P[i],P[0]);
        PO = P[0];
        sort(P.Points.begin(),P.Points.end(),jijiaopaixu);
        res.add(P[0]);
        res.add(P[1]);
        int top = 1;
        for(int i = 2;i < n;++i){
            while(top >= 1 && Cross(res[top] - res[top-1],P[i] - res[top-1]) < 0){
                top--;
                res.del();
            }
            res.add(P[i]);
            top++;
        }
        return res;

    }

    TPLT    /// 圆的结构体
    struct Circle{
        VTT C;T R;
        Circle(){}
        Circle(T C,T R):C(C),R(R){}
        VTT getPoint(double rad){return VTT(C.x + cos(rad) * R,C.y + sin(rad) * R);}
    };
    /* 求两个圆的交点，返回值为一个vector , vector[0].R 表示状态
        -1   0      1      2    3
        重合 有交点 同心圆 外离 内含
    */
    TPLT double Angle(VTT a){if(a.y >= 0)return Angle(a,VTT(1,0));return 2 * PI - Angle(a,VTT(1,0));}
    TPLT vector<CCT> GetCC(CCT c1,CCT c2){
        vector<CCT> ans;
        double d = Length(c1.c - c2.c);
        if(dcmp(d) == 0){
            if(dcmp(c1.r - c2.r) == 0){ans.push_back(CCT(VTT(0,0),-1));}  // 重合
            else{ans.push_back(CCT(VTT(0,0),1));}   // 同心圆
        }else{
            if(dcmp(c1.r + c2.r - d) < 0){ans.push_back(CCT(VTT(0,0),2));}  // 外离
            else if(dcmp(fabs(c1.r - c2.r) - d > 0)){ans.push_back(CCT(VTT(0,0),3));}    // 内含
            else{
                ans.push_back(CCT(VTT(0,0),0));   // 有交点
                double rad = Angle(c2.C - c2.C);
                double da = acos((c1.R * c1.R + d * d - c2.R * c2.R) / (2.0 * c1.R * d));
                VTT p1 = c1.getPoint(rad - da),p2 = c1.getPoint(rad + da);
                ans.push_back(p1);
                if(p1 != p2)ans.push_back(p2);
            }
        }
        return ans;
    }
    /* 不共线三点 求圆心(外心) */
    TPLT inline T pow2(const T &x){return x * x;}
    TPLT VTT circumcenter(VTT p1,VTT p2,VTT p3){
        double a = p1.x - p2.x;
        double b = p1.y - p2.y;
        double c = p1.x - p3.x;
        double d = p1.y - p3.y;
        double e = (pow2(p1.x) - pow2(p2.x) +
        pow2(p1.y) - pow2(p2.y)) / 2;
        double f = (pow2(p1.x) - pow2(p3.x) +
        pow2(p1.y) - pow2(p3.y)) / 2;
        return VTT((d * e - b * f) /
        (a * d - b * c),
        (a * f - c * e) /
        (a * d - b * c));
    }
}
using namespace My_Geometry;