计算几何板子

Posted 2021-02-06 h-riven

struct Point {
    double x, y;
    Point() {}
    Point(double a, double b) {
        x = a, y = b;
    }
    bool operator == (Point a) {
        return sgn(x-a.x) == 0 && sgn(y-a.y) == 0;
    }
    bool operator < (Point a) {     //最左下角
        return sgn(y-a.y) == 0 ? sgn(x-a.x)<0 : y<a.y;
    }
    Point operator - (Point a) {
        return Point(x-a.x, y-a.y);     //向量减
    }
    Point operator + (Point a) {
        return Point(x+a.x, y+a.y);     //向量加
    }
    double operator * (Point a) {   //点乘(内积)
        return x*a.x + y*a.y;
    }
    double operator ^ (Point a) {   //叉乘(外积)
        return x*a.y - y*a.x;
    }
    Point operator * (double k) {
        return Point(x*k, y*k);
    }
    Point operator / (double k) {
        return Point(x/k, y/k);
    }
    Point trunc(double r) {
        //化为长度为 r 的向量
        double l = len();
        if(!sgn(l)) return *this;
        r /= l;
        return Point(x*r, y*r);
    }
    Point rotleft() {
        //原向量逆时针旋转90度
        return Point(-y, x);
    }
    Point rotright() {
        //原向量顺时针旋转90度
        return Point(y, -x);
    }
    Point rotate(Point p, double rad) {
        //原向量绕着 p 点逆时针旋转 rad 度
        Point v = (*this) - p;
        double c = cos(rad), s = sin(rad);
        return Point(p.x + v.x*c - v.y*s, p.y + v.x*s + v.y*c);
    }
    double len() {          //向量长度
        return hypot(x, y);
    }
    double len2() {         //向量长度平方
        return x*x + y*y;
    }
    double distance(Point p) {  //两点之间距离
        return hypot(x-p.x, y-p.y);
    }
    double rad(Point a, Point b) {
        //三个点:p, a, b
        //计算 pa 与 pb 所成的夹角
        Point p = *this;
        return fabs(atan2(fabs((a-p) ^ (b-p)), (a-p) * (b-p)));
    }
};

struct Line {
    Point s, e;
    Line() {};
    Line(Point _s, Point _e) {
        // 两点确定直线
        s = _s, e = _e;
    }
    Line(Point p, double angle) {
        //根据点 p 和倾斜角 angle 确定直线, 0 <= angle <= pi
        s = p;
        if(sgn(angle - pi/2) == 0) {
            e = (s + Point(0, 1));
        } else {
            e = (s + Point(1, tan(angle)));
        }
    }
    Line(double a, double b, double c) {
        //根据 ax+by+c = 0 确定直线
        if(sgn(a) == 0) {
            s = Point(0, -c/b);
            e = Point(1, -c/b);
        } else if(sgn(b) == 0) {
            s = Point(-c/a, 0);
            e = Point(-c/a, 1);
        } else {
            s = Point(0, -c/b);
            e = Point(1, (-c-a)/b);
        }
    }
    void adjust() {
        if(e < s)   swap(s, e);
    }
    bool operator == (Line v) {
        return (s == v.s) && (e == v.e);
    }
    double length() {
        return s.distance(e);
    }
    double angle() {
        // 直线倾斜角
        double k = atan2(e.y-s.y, e.x-s.x);
        if(sgn(k) < 0)  k += pi;
        if(sgn(k-pi) == 0)  k -= pi;
        return k;
    }
    int relation(Point p) {
        /*
        判断直线和点关系
        this line, p point
        1   在左侧
        2   在右侧
        3   在直线上
        */
        int c = sgn((p-s) ^ (e-s));
        if(c < 0)   return 1;
        else if(c > 0)  return 2;
        else    return 3;
    }
    bool pointonseg(Point p) {
        /*
        判断点和线段关系
        this seg, p point
        0   不在线段上
        1   在线段上
        */
        return sgn((p-s) ^ (e-s)) == 0 && sgn((p-s) * (p-e)) <= 0;
    }
    bool parallel(Line v) {
        //判断两向量平行
        return sgn((e-s) ^ (v.e-v.s)) == 0;
    }
    int segcrossseg(Line v) {
        /*
        两线段关系
        2   规范相交
        1   不规范相交
        0   不相交
        */
        int d1 = sgn((e-s) ^ (v.s-s));
        int d2 = sgn((e-s) ^ (v.e-s));
        int d3 = sgn((v.e-v.s) ^ (s-v.s));
        int d4 = sgn((v.e-v.s) ^ (e-v.s));
        if((d1 ^ d2) == -2 && (d3 ^ d4) == -2)
            return 2;
        return (d1 == 0 && sgn((v.s-s) * (v.s-e)) <= 0) ||
               (d2 == 0 && sgn((v.e-s) * (v.e-e)) <= 0) ||
               (d3 == 0 && sgn((s-v.s) * (s-v.e)) <= 0) ||
               (d4 == 0 && sgn((e-v.s) * (e-v.e)) <= 0);
    }
    int linecrossseg(Line v) {
        /*
        判断直线和线段关系
        this line, v seg
        2   规范相交
        1   非规范相交
        0   不相交
        */
        int d1 = sgn((e-s) ^ (v.s-s));
        int d2 = sgn((e-s) ^ (v.e-s));
        if((d1 ^ d2) == -2) return 2;
        else    return (d1 == 0 || d2 == 0);
    }
    int linecrossline(Line v) {
        /*
        两直线关系
        0   平行
        1   重合
        2   相交
        */
        if((*this).parallel(v))
            return v.relation(s) == 3;
        return 2;
    }
    Point crosspoint(Line v) {
        //求两直线的交点
        //要求两直线不平行，不重合
        double a1 = (v.e-v.s) ^ (s-v.s);
        double a2 = (v.e-v.s) ^ (e-v.s);
        return Point((s.x*a2-e.x*a1), (s.y*a2-e.y*a1))/(a2-a1);
    }
    double dispointtoline(Point p) {
        //求点到直线的距离
        return fabs((p-s)^(e-s))/length();
    }
    double dispointtoseg(Point p) {
        if(sgn((p-s)*(e-s))<0 || sgn((p-e)*(s-e)) < 0)
            return min(p.distance(s), p.distance(e));
        return dispointtoline(p);
    }
    double dissegtoseg(Line v) {
        /*
        返回线段到线段的距离
        前提是两线段不相交，相交为0
        */
        return min(min(dispointtoseg(v.s), dispointtoseg(v.e)),
                   min(v.dispointtoseg(s), v.dispointtoseg(e)));
    }
    Point lineprog(Point p) {
        //返回点 p 在直线上的投影点
        return s + (((e-s)*((e-s)*(p-s))) / ((e-s).len2()));
    }
    Point symmetrypoint(Point p) {
        //返回点 p 关于直线p的对称点
        Point q = lineprog(p);
        return Point(2*q.x-p.x, 2*q.y-p.y);
    }
};

struct Circle {
    Point p;    //圆心
    double r;   //半径
    Circle() {}
    Circle(Point _p, double _r) {
        p = _p, r = _r;
    }
    Circle(double x, double y, double _r) {
        p = Point(x, y), r = _r;
    }
    Circle(Point a, Point b, Point c) {
        //求三角形的外接圆
        Line u = Line((a+b)/2, ((a+b)/2) + ((b-a).rotleft()));
        Line v = Line((b+c)/2, ((b+c)/2) + ((c-b).rotleft()));
        p = u.crosspoint(v);
        r = p.distance(a);
    }
    Circle(Point a, Point b, Point c, bool t) {
        //求三角形的内接圆
        //bool t用来区别三角形外接圆的函数
        Line u, v;
        double m = atan2(b.y-a.y, b.x-a.x);
        double n = atan2(c.y-a.y, c.x-a.x);
        u.s = a;
        u.e = u.s + Point(cos((n+m)/2), sin((n+m)/2));
        v.s = b;
        m = atan2(a.y-b.y, a.x-b.x);
        n = atan2(c.y-b.y, c.x-b.x);
        v.e = v.s + Point(cos((n+m)/2), sin((n+m)/2));
        p = u.crosspoint(v);
        r = Line(a, b).dispointtoseg(p);
    }
    bool operator == (Circle v) {
        return p==v.p && sgn(r-v.r)==0;
    }
    bool operator < (Circle v) {
        return (p<v.p) || (p==v.p && sgn(r-v.r)<0);
    }
    double area() {
        //面积
        return pi*r*r;
    }
    double cirperimeter() {
        //周长
        return 2*pi*r;
    }
    int relation(Point b) {
        /*
        点和圆的关系
        0   圆外
        1   圆上
        2   圆内
        */
        double dst = b.distance(p);
        if(sgn(dst-r) < 0)  return 2;
        else if(sgn(dst-r) == 0)    return 1;
        else return 0;
    }
    int relationseg(Line v) {
        /*
        线段和圆的关系
        0   圆外
        1   圆上
        2   圆内
        */
        double dst = v.dispointtoseg(p);
        if(sgn(dst-r) < 0)  return 2;
        else if(sgn(dst-r) == 0)    return 1;
        else return 0;
    }
    int relationline(Line v) {
        /*
        直线和圆的关系
        0   圆外
        1   圆上
        2   圆内
        */
        double dst = v.dispointtoline(p);
        if(sgn(dst-r) < 0)  return 2;
        else if(sgn(dst-r) == 0)    return 1;
        else return 0;
    }
    int relationcircle(Circle v) {
        /*
        两圆的关系
        5   相离
        4   外切
        3   相交
        2   内切
        1   内含
        */
        double d = p.distance(v.p);
        if(sgn(d-r-v.r) > 0)    return 5;
        if(sgn(d-r-v.r) == 0)   return 4;
        double l = fabs(r-v.r);
        if(sgn(d-r-v.r) < 0 && sgn(d-l) > 0)    return 3;
        if(sgn(d-l) == 0)       return 2;
        if(sgn(d-l) < 0)        return 1;
    }
    int pointcrosscircle(Circle v, Point &p1, Point &p2) {
        //求两个圆的交点，返回交点数
        int rel = relationcircle(v);
        if(rel == 1 || rel == 5)    return 0;
        double d = p.distance(v.p);
        double l = (d*d + r*r - v.r*v.r) / (2*d);
        double h = sqrt(r*r-l*l);
        Point tmp = p + (v.p-p).trunc(l);
        p1 = tmp + ((v.p-p).rotleft().trunc(h));
        p2 = tmp + ((v.p-p).rotright().trunc(h));
        if(rel == 2 || rel == 4)
            return 1;
        return 2;
    }
    int pointcrossline(Line v, Point &p1, Point &p2) {
        //求直线和圆的交点，返回交点数
        if(!(*this).relationline(v))    return 0;
        Point a = v.lineprog(p);
        double d = v.dispointtoline(p);
        d = sqrt(r*r-d*d);
        if(sgn(d) == 0) {
            p1 = a;
            p2 = a;
            return 1;
        }
        p1 = a + (v.e-v.s).trunc(d);
        p2 = a - (v.e-v.s).trunc(d);
        return 2;
    }
    int getcircle(Point a, Point b, double r1, Circle &c1, Circle &c2) {
        //过a, b两点, 半径为 r1 的圆, 返回圆数
        Circle x(a, r1), y(b, r1);
        int t = x.pointcrosscircle(y, c1.p, c2.p);
        if(!t)  return 0;
        c1.r = c2.r = r1;
        return t;
    }
    int getcircle(Line u, Point q, double r1, Circle &c1, Circle &c2) {
        //与直线 u 相切, 过点 q, 半径为 r1 的圆,返回圆数
        double dis = u.dispointtoline(q);
        if(sgn(dis-2*r1) > 0)   return 0;
        if(sgn(dis) == 0) {
            c1.p = q + ((u.e-u.s).rotleft().trunc(r1));
            c2.p = q + ((u.e-u.s).rotright().trunc(r1));
            c1.r = c2.r = r1;
            return 2;
        }
        Line u1 = Line((u.s + (u.e-u.s).rotleft().trunc(r1)), (u.e + (u.e-u.s).rotleft().trunc(r1)));
        Line u2 = Line((u.s + (u.e-u.s).rotright().trunc(r1)), (u.e + (u.e-u.s).rotright().trunc(r1)));
        Circle cc = Circle(q, r1);
        Point p1, p2;
        if(!cc.pointcrossline(u1, p1, p2))
            cc.pointcrossline(u2, p1, p2);
        c1 = Circle(p1, r1);
        if(p1 == p2) {
            c2 = c1;
            return 1;
        }
        c2 = Circle(p2, r1);
        return 2;
    }
    int getcircle(Line u, Line v, double r1, Circle &c1, Circle &c2, Circle &c3, Circle &c4) {
        //同时与直线u, v相切, 半径为 r1 的圆
        if(u.parallel(v))return 0;//两直线平行
        Line u1 = Line(u.s + (u.e-u.s).rotleft().trunc(r1),u.e + (u.e-u.s).rotleft().trunc(r1));
        Line u2 = Line(u.s + (u.e-u.s).rotright().trunc(r1),u.e + (u.e-u.s).rotright().trunc(r1));
        Line v1 = Line(v.s + (v.e-v.s).rotleft().trunc(r1),v.e + (v.e-v.s).rotleft().trunc(r1));
        Line v2 = Line(v.s + (v.e-v.s).rotright().trunc(r1),v.e + (v.e-v.s).rotright().trunc(r1));
        c1.r = c2.r = c3.r = c4.r = r1;
        c1.p = u1.crosspoint(v1);
        c2.p = u1.crosspoint(v2);
        c3.p = u2.crosspoint(v1);
        c4.p = u2.crosspoint(v2);
        return 4;
    }
    int getcircle(Circle cx, Circle cy, double r1, Circle &c1, Circle &c2) {
        //同时与不相交圆cx, cy相切，半径为 r1 的圆
        Circle x(cx.p, r1+cx.r), y(cy.p, r1+cy.r);
        int t = x.pointcrosscircle(y, c1.p, c2.p);
        if(!t)  return 0;
        c1.r = c2.r = r1;
        return t;
    }
    int tangentline(Point q, Line &u, Line &v) {
        //过一点做圆的切线,求切线
        int x = relation(q);
        if(x == 2)return 0;
        if(x == 1) {
            u = Line(q, q + (q-p).rotleft());
            v = u;
            return 1;
        }
        double d = p.distance(q);
        double l = r*r/d;
        double h = sqrt(r*r-l*l);
        u = Line(q, p + ((q-p).trunc(l) + (q-p).rotleft().trunc(h)));
        v = Line(q, p + ((q-p).trunc(l) + (q-p).rotright().trunc(h)));
        return 2;
    }
    double areacircle(Circle v) {
        //两相交圆的面积
        int rel = relationcircle(v);
        if(rel >= 4)    return 0.0;
        if(rel <= 2)    return min(area(), v.area());
        double d = p.distance(v.p);
        double hf = (r+v.r+d)/2.0;
        double ss = 2*sqrt(hf*(hf-r)*(hf-v.r)*(hf-d));
        double a1 = acos((r*r+d*d-v.r*v.r) / (2.0*r*d));
        a1 = a1*r*r;
        double a2 = acos((v.r*v.r+d*d-r*r) / (2.0*v.r*d));
        a2 = a2*v.r*v.r;
        return a1+a2-ss;
    }
    double areatriangle(Point a, Point b) {
        //求圆和三角形 pab 的相交面积
        if(sgn((p-a)^(p-b)) == 0)   return 0.0;
        Point q[5];
        int len = 0;
        q[len++] = a;
        Line l(a,b);
        Point p1,p2;
        if(pointcrossline(l,q[1],q[2]) == 2) {
            if(sgn((a-q[1])*(b-q[1])) < 0)  q[len++] = q[1];
            if(sgn((a-q[2])*(b-q[2])) < 0)  q[len++] = q[2];
        }
        q[len++] = b;
        if(len == 4 && sgn((q[0]-q[1])*(q[2]-q[1])) > 0)
            swap(q[1], q[2]);
        double res = 0;
        for(int i = 0; i < len-1; i++) {
            if(relation(q[i]) == 0 || relation(q[i+1]) == 0) {
                double arg = p.rad(q[i], q[i+1]);
                res += r*r*arg/2.0;
            } else {
                res += fabs((q[i]-p)^(q[i+1]-p)) / 2.0;
            }
        }
        return res;
    }
};

struct Polygon {
    int n;
    Point p[maxn];
    Line l[maxn];
    Polygon() {
        n = 0;
    }
    void add(Point q) {
        p[++n] = q;
    }
    void getline() {
        for(int i=1; i<=n; i++) {
            l[i] = Line(p[i], p[i%n+1]);
        }
    }
    struct cmp {
        //极角排序的比较函数
        Point p;
        cmp(Point _p) {
            p = _p;
        }
        bool operator()(Point _a, Point _b) const {
            Point a = _a, b = _b;
            int d = sgn((a-p) ^ (b-p));
            if(d == 0) {
                return sgn(a.distance(p) - b.distance(p)) < 0;
            } else {
                return d > 0;
            }
        }
    };
    void norm() {
        //极角排序
        int id = 1;
        for(int i=2; i<=n; i++) {
            if(p[i] < p[id])
                id = i;
        }
        swap(p[id], p[1]);
        sort(p+1, p+1+n, cmp(p[1]));
    }
    void Graham(Polygon &convex) {  //引用
        //求解凸包
        //convex凸包
        norm();
        mes(convex.p, 0);
        int &top = convex.n = 0;    //引用
        if(n == 1) {
            convex.p[++top] = p[1];
        } else if(n == 2) {
            convex.p[++top] = p[1];
            convex.p[++top] = p[2];
            if(convex.p[1] == convex.p[2]) top--;
        } else {
            convex.p[++top] = p[1];
            convex.p[++top] = p[2];
            for(int i=3; i<=n; i++) {
                while(top>1 && sgn((convex.p[top]-convex.p[top-1]) ^ (p[i]-convex.p[top-1])) <= 0)
                    top--;
                convex.p[++top] = p[i];
            }
            if(top == 2 && convex.p[1] == convex.p[2])
                top--;
        }
    }
    bool isconvex() {
        //判断多边形是不是凸的
        //1是凸的 0不是凸的
        bool s[3];
        mes(s, false);
        for(int i=1; i<=n; i++) {
            int j = i%n + 1;
            int k = j%n + 1;
            s[sgn((p[j]-p[i]) ^ (p[k]-p[i])) + 1] = true;
            if(s[0] && s[2])    return false;
        }
        return true;
    }
    int relationpoint(Point q) {
        /*
        点和多边形关系
        3   点上
        2   边上
        1   内部
        0   外部
        */
        for(int i=1; i<=n; i++) {
            if(p[i] == q) return 3;
        }
        getline();
        for(int i=1; i<=n; i++) {
            if(l[i].pointonseg(q))  return 2;
        }
        int cnt = 0;
        for(int i=1; i<=n; i++) {
            int j = i%n + 1;
            int k = sgn((q-p[j]) ^ (p[i]-p[j]));
            int u = sgn(p[i].y - q.y);
            int v = sgn(p[j].y - q.y);
            if(k > 0 && u < 0 && v >= 0)    cnt++;
            if(k < 0 && v < 0 && u >= 0)    cnt--;
        }
        return cnt != 0;
    }
    bool getdir() {
        //得到方向
        //1表示逆时针，0表示顺时针
        double sum = 0;
        for(int i=1; i<=n; i++) {
            sum += p[i] ^ p[i%n+1];
        }
        if(sgn(sum) > 0)    return 1;
        return 0;
    }
    Point getbarycentre() {
        //得到重心
        Point ans(0, 0);
        double area = 0;
        for(int i=2; i<=n-1; i++) {
            double tmp = (p[i]-p[1]) ^ (p[i+1]-p[0]);
            if(sgn(tmp) == 0)   continue;
            area += tmp;
            ans.x = (p[1].x + p[i].x + p[i+1].x) / 3 * tmp;
            ans.y = (p[0].y + p[i].y + p[i+1].y) / 3 * tmp;
        }
        if(sgn(area))   ans = ans / area;
        return ans;
    }
    void convexcut(Line u, Polygon &po) {
        //直线切割凸多边形左侧,注意直线方向
        int &top = po.n = 0;    //引用
        for(int i=1; i<=n; i++) {
            int d1 = sgn((u.e-u.s) ^ (p[i]-u.s));
            int d2 = sgn((u.e-u.s) ^ (p[i%n+1]-u.s));
            if(d1 >= 0) po.p[++top] = p[i];
            if(d1*d2<0) po.p[++top] = u.crosspoint(Line(p[i], p[i%n+1]));
        }
    }
    double getperimeter() {
        //多边形的周长
        double ans = 0.0;
        if(n == 1) {
            return ans;
        } else if(n == 2) {
            ans += p[1].distance(p[2]);
        } else {
            for(int i=1; i<=n; i++) {
                ans += p[i].distance(p[i%n+1]);
            }
        }
        return ans;
    }
    double getarea() {
        //多边形面积
        double sum = 0;
        for(int i=1; i<=n; i++) {
            sum += (p[i] ^ p[i%n+1]);
        }
        return fabs(sum)/2.0;
    }
    double areacircle(Circle c) {
        //多边形和圆相交的面积
        double ans = 0;
        for(int i=1; i<=n; i++) {
            int j = i%n+1;
            if(sgn((p[j]-c.p) ^ (p[i]-c.p)) >= 0)
                ans += c.areatriangle(p[i], p[j]);
            else
                ans -= c.areatriangle(p[i], p[j]);
        }
        return fabs(ans);
    }
    int relationcircle(Circle c) {
        /*
        多边形和圆的关系
        2   圆完全在多边形内
        1   圆在多边形内，但是碰到了多边形边界
        0   其他
        */
        getline();
        int x = 2;
        if(relationpoint(c.p) != 1) return 0;
        for(int i=1; i<=n; i++) {
            if(c.relationseg(l[i]) == 2)    return 0;
            if(c.relationseg(l[i]) == 1)    x = 1;
        }
        return x;
    }
};