其他分享
首页 > 其他分享> > 计算几何2D

计算几何2D

作者：互联网

零. 三态函数

const double eps = 1e-8;
const double Pi = acos(-1.0);
int sgn(double x) {   
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}

一. 点类及其常见操作

typedef struct Point {
    double x, y;
    Point(double x=0, double y=0): x(x), y(y) {}
    void input() { scanf("%lf%lf", &x, &y); }
    double Len() { return sqrt(x*x + y*y); }  // 模长
    double Len2() { return x*x + y*y; }   // 模长平方
    Point operator + (Point r) { return Point(x+r.x, y+r.y); }  // 向量相加
    Point operator - (Point r) { return Point(x-r.x, y-r.y); }  // 向量相减
    Point operator * (double r) { return Point(x*r, y*r); }     // 向量模改
    Point operator / (double r) { return Point(x/r, y/r); }     // 向量模改
    double operator ^ (Point r) { return x*r.y - r.x*y; }       // 叉积
    double operator * (Point r) { return x*r.x + y*r.y; }       // 点积
    bool operator == (Point r) { return sgn(x-r.x)==0 && sgn(y-r.y)==0; }
    bool operator < (Point r) { return sgn(y-r.y)<0 || (sgn(y-r.y)==0&&sgn(x-r.x)<0); }
    void transXY(double rad) {   // 点绕原点旋转 rad 弧度
        double tx = x, ty = y;
        x = tx*cos(rad) - ty*sin(rad);
        y = tx*sin(rad) + ty*cos(rad);
    }
}Vector;

  1. 点积

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

  2. 叉积

    double Cross(Point a, Point b) { return a.x*b.y - b.x*a.y; }
double Cross(Point pt, Point a, Point b) { return Cross(a-pt, b-pt); }

  3. 向量的旋转（逆时针为正方向）

    Vector Ratate(Vector v, double rad) { return v.transXY(rad), v; }

  4. 向量的垂直向量

    Vector Normal(Vector) { return Vector(-v.y, v.x); }

二. 线类（两点式）

struct Line {
    Point s, e;
    Line(Point s=o, Point e=o): s(s), e(e) {}
    void input() { s.input(), e.input(); }
    pair<int, Point> operator & (Line r) {
        Point res = s;
        if(sgn(Cross(r.e - r.s, e - s)) == 0) {
            if(sgn(Cross(r.e - r.s, e - r.s)) == 0) return make_pair(1, res);
            return make_pair(2, res);
        }
        double t = Cross(s-r.s, r.e-r.s)/Cross(r.e-r.s, e-s);
        res.x += (e.x-s.x) * t;
        res.y += (e.y-s.y) * t;
        return make_pair(0, res);
    }
};

  1. 线段相交

    bool Seg_inter_Seg(Line A, Line B) {
    return
        max(A.s.x, A.e.x) >= min(B.s.x, B.e.x) &&
        max(A.s.y, A.e.y) >= min(B.s.y, B.e.y) &&
        max(B.s.x, B.e.x) >= min(A.s.x, A.e.x) &&
        max(B.s.y, B.e.y) >= min(A.s.y, A.e.y) &&
        sgn(Cross(A.e-A.s, B.e-A.s)) * sgn(Cross(A.e-A.s, B.s-A.s)) <= 0 &&
        sgn(Cross(B.e-B.s, A.e-B.s)) * sgn(Cross(B.e-B.s, A.s-B.s)) <= 0;
}

  2. 点是否在向量上

    bool OnSeg(Point pt, Line L) {
    return sgn(Cross(L.s - pt, L.e - pt)) == 0
        &&  sgn(Dot(L.s - pt, L.e - pt)) <= 0;
}

  3. 判断直线和线段是否相交

    bool Seg_inter_Line(Line seg, Line L) {
    return sgn(Cross(seg.s-L.s, L.e-L.s)) * sgn(Cross(seg.e-L.s, L.e-L.s)) <= 0;
}

  4. 两点间的距离

    double Dist(Point A, Point B) { return (A-B).Len(); }

  5. 点到直线的距离

    double Point_to_Line(Point pt, Line L) {
    return fabs(Cross(pt-L.s, L.s-L.e)/Dist(L.s, L.e));
}

  6. 点到线段的距离

    double Point_to_Seg(Point pt, Line L) {
    if(sgn(Dot(pt-L.s, L.e-L.s)) < 0) return Dist(pt, L.s);
    if(sgn(Dot(pt-L.e, L.s-L.e)) < 0) return Dist(pt, L.e);
    return Point_to_Line(pt, L);
}

  7. 两线段的距离

    double Seg_to_Seg(Line A, Line B) {
    if (Seg_inter_Seg(A, B)) return 0;
    return min(Point_to_Seg(A.e, B), Point_to_Seg(A.s, B));
}

三. 线类（点向式）

struct Plane {
    Point p; // 直线上的一点
    Vector v; // 方向向量
    double ang; // 从 x 轴转向 v 方向的夹角
    Plane(Point p=o, Vector v=Vector(1,0)): p(p), v(v) { ang = atan2(v.y, v.x); }
    void Angle() { ang = atan2(v.y, v.x); }
    bool operator < (Plane r) { return sgn(ang - r.ang) > 0; }
};

  1. 判断点 p p p 在直线 L L L 的左边，即点 p p p 在 L L L 的外边。

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

  2. 判断点是否在直线上

    bool OnLine(Plane L, Point pt) { return sgn(Cross(pt - L.p, L.v)) == 0; }

  3. 判断点是否不在直线的右边

    bool NotOnRight(Plane L, Point pt) { return sgn(Cross(L.v, pt - L.p)) >= 0;  }

  4. 两直线的交点

    Point Cross_point(Plane A, Plane B) {
    Vector u = A.p - B.p;
    double t = Cross(B.v, u)/Cross(A.v, B.v);
    return A.p + A.v * t;
}

  5. 计算到固定点的距离和

    double Dist_and(Point pt, Point* p, int n) {
    double res = 0;
    for(int i = 0; i < n; i ++)
        res += Dist(pt, p[i]);
    return res;
}

四. 圆类

struct Circle{
    Point p; double r;
    Circle(Point p=o, double r=0): p(p), r(r) {}
    void input() { p.input(); scanf("%lf", &r); }
};

  1. 过圆外一点与圆的切线

    Plane Tangent(Point p, Circle c, int clock) {
    Plane res;
    Vector v = c.p - p;
    double d = Dist(c.p, p);
    double rad = asin(c.r/d);
    if(clock > 0) res.v = Ratate(v, rad);
    else res.v = Ratate(v, -rad);
    res.p = p;
    return res;
}

五. 凸包相关

  1. 计算凸包的周长

    double Perimeter_polygon(Point* p, int n) {
    p[n] = p[0];
    double ans = 0;
    for(int i = 0; i < n; i ++)
        ans += Dist(p[i], p[i+1]);
    return ans;
}

  2. 计算凸包的面积

    double Area_polygon(Point* p, int n) {
    p[n] = p[0];
    double ans = 0;
    for(int i = 0; i < n; i ++)
        ans += Cross(p[i], p[i+1]);
    return ans;
}

  3. 判断圆是否在凸多边形内部

    bool Circle_in_polygon(Circle O, Point* p, int n) {
    p[n] = p[0];
    for(int i = 0; i < n; i ++) {
        if(sgn(Point_to_Seg(O.p, Line(p[i], p[i+1])) - O.r) < 0)
            return false;
    }
    return true;
}

  4. 判断是否为凸包，方向不一定是逆时针

    bool Judge_Convex_hull(Point* p, int n){
    int dir = 0;
    p[n] = p[0], p[n+1] = p[1];
    for(int i = 1; i <= n; i ++) {
        int u = sgn(Cross(p[i] - p[i-1], p[i+1] - p[i]));
        if(! dir) dir = u;
        if(u * dir < 0) return false;
    }
    return true;
}

  5. 将凸包变为逆时针顺序

    void Anti_clockwise(Point* p, int n) {
    p[n] = p[0];
    for(int i = 0; i < n - 1; i ++) {
        if(sgn(Cross(p[i], p[i+1], p[i+2])) > 0) return ;
        else if(sgn(Cross(p[i], p[i+1], p[i+2])) < 0) {
            for(int j = 0; j < n/2; j ++) {
                Point pt = p[j];
                p[j] = p[n-j-1], p[n-j-1] = pt;
            }
            return ;
        }
    }
}

  6. 判断点是否在多边形（不限于凸包）内部

    int in_polygon(Point pt, Point* p, int n) {
    for(int i = 0; i < n; i ++)
        if(pt == p[i]) return 3;
    for(int i = 0; i < n; i ++) {
        int j = i + 1 == n ? 0 : i + 1;
        if(OnSeg(pt, Line(p[i], p[j]))) return 2;
    }
    int num = 0;
    for(int i = 0; i < n; i ++) {
        int j = i + 1 == n ? 0 : i + 1;
        int c = sgn(Cross(pt-p[j], p[i]-p[j]));
        int u = sgn(p[j].y - pt.y);
        int v = sgn(p[i].y - pt.y);
        if(c > 0 && u >= 0 && v < 0) num ++;
        if(c < 0 && u < 0 && v >= 0) num --;
    }
    return num != 0;
}

  7. 求凸包（Andrew算法）

    int Convex_hull(Point* p, Point* ch, int n) {
    sort(p, p + n);
    int v = 0;
    for(int i = 0; i < n; i ++) {
        while(v > 1 && sgn(Cross(ch[v-1]-ch[v-2], p[i]-ch[v-2])) <= 0) v --;
        ch[v ++] = p[i];
    }
    int j = v;
    for(int i = n - 2; ~ i; i --) {
        while(v > j && sgn(Cross(ch[v-1]-ch[v-2], p[i]-ch[v-2])) <= 0) v --;
        ch[v ++] = p[i];
    }
    if(n > 1) v --;
    return v;
}

  8. 两个多边形的面积交，两个凸包都是逆时针

    double CPI_Area(Point* a, int na, Point* b, int nb) {
    static Point p[N], tmp[N];
    int tn, sflag, eflag;
    a[na] = a[0], b[nb] = b[0];
    memcpy(p, b, sizeof(Point)*(nb+1));
    /// 整个过程类似半平面交
    /// 枚举一个凸包的所有边对另一个凸包进行切割
    for(int i = 0; i < na && nb > 2; i ++) {
        sflag = sgn(Cross(a[i], a[i+1], p[0]));
        for(int j = tn = 0; j < nb; j ++, sflag = eflag) {
            if(sflag >= 0) tmp[tn ++] = p[j]; /// 当前顶点不在右侧, 直接保留
            eflag = sgn(Cross(a[i], a[i+1], p[j+1]));
            /// 异号说明两直线相交, 计算交点切割边
            if((sflag ^ eflag) == -2)
                tmp[tn ++] = (Line(a[i], a[i+1]) & Line(p[j], p[j+1])).second;
        }
        memcpy(p, tmp, sizeof(Point)*tn);
        nb = tn, p[nb] = p[0];
    }
    if(nb < 3) return 0;
    /// 返回切割后的面积就是面积交
    return Area_polygon(p, nb);
}

  9. 两个多边形的面积并

    double SPI_Area(Point* a, int na, Point* b, int nb) {
    Point t1[4], t2[4];
    double res = 0;int num1, num2;
    a[na] = t1[0] = a[0], b[nb] = t2[0] = b[0];
    /// 以各自的 0 处的点作为基本点, 并保证整个三角形都是都是逆时针
    for(int i = 2; i < na; i ++) {
        t1[1] = a[i-1], t1[2] = a[i];
        num1 = sgn(Cross(t1[0], t1[1], t1[2]));
        if(num1 < 0) swap(t1[1], t1[2]);
        for(int j = 2; j < nb; j ++) {
            t2[1] = b[j-1], t2[2] = b[j];
            num2 = sgn(Cross(t2[0], t2[1], t2[2]));
            if(num2 < 0) swap(t2[1], t2[2]);
            res += CPI_Area(t1, 3, t2, 3) * num1 * num2;
        }
    }
    return fabs(res);
}

六. 半平面交

  1. 算法实现

    int HPI(Plane* L, Point *ch, int n) {
    sort(L, L + n); /// 极角排序
    int h = 0, t = 0, v = 0;
    static Plane q[N]; /// 双端队列
    static Point p[N]; /// 两个相邻半平面的交点
    q[h] = L[0];
    for(int i = 1; i < n; i ++) {
        /// 删除队尾的半平面
        while(h < t && ! OnLeft(L[i], p[t-1])) t --;
        /// 删除队首的半平面
        while(h < t && ! OnLeft(L[i], p[h])) h ++;
        q[++ t] = L[i]; /// 将当前的半平面加入双端队列的尾部
        /// 极角相同的两个半平面保留左边
        if(sgn(Cross(q[t].v, q[t-1].v)) == 0) {
            t --;
            if(OnLeft(q[t], L[i].p)) q[t] = L[i];
        }
        /// 计算队列尾部的半平面交点
        if(h < t) p[t - 1] = Cross_point(q[t-1], q[t]);
    }
    /// 删除队列尾部无用的半平面
    while(h < t && ! OnLeft(q[h], p[t-1])) t --;
    if(t - h <= 1) return 0; /// 空集
    p[t] = Cross_point(q[t], q[h]); /// 计算队列首尾部的交点
    for(int i = h; i <= t; i ++) ch[v ++] = p[i]; /// 复制
    return v; /// 返回凸多边形大小
}

  2. 判断一个点是否有核（单独一个点也算）

    bool HPI_Core(Plane* L, int n){
    sort(L, L + n);
    int h = 0, t = 0;
    static Plane q[N];
    static Point p[N];
    q[0] = L[0];
    for(int i = 1; i < n; i ++) {
        while(h < t && ! OnLeft(L[i], p[t-1])) t --;
        while(h < t && ! OnLeft(L[i], p[h])) h ++;
        q[++ t] = L[i];
        if(sgn(Cross(q[t].v, q[t-1].v)) == 0) {
            t --;
            if(OnLeft(q[t], L[i].p)) q[t] = L[i];
        }
        if(h < t) p[t-1] = Cross_point(q[t], q[t-1]);
    }
    while(h < t && ! OnLeft(q[h], p[t-1])) t --;
    return t - h + 1 > 2;
}

七. 旋转卡壳

  1. 旋转卡壳求凸包直径（宽度）

    double Rotating_cailper(Point* p, int n) {
    double res = 0;
    if(n == 2) return Dist(p[0], p[1]);
    p[n] = p[0];
    for(int i = 0, j = 2; i < n; i ++) {
        while(sgn(Cross(p[i], p[i+1], p[j]) - Cross(p[i], p[i+1], p[j+1])) < 0)
            if(++ j >= n) j -= n;
        /// 计算凸包宽度时取对踵点到直线的距离
        /// res = max(res, Point_to_Line(p[j], Line(p[i], p[i+1])));
        res = max(res, max(Dist(p[j], p[i]), Dist(p[j], p[i+1])));
    }
    return res;
}

  2. 旋转卡壳求两个凸包的最近距离

    double Nearest_dist(Point* p, int n, Point* q, int m) {
    p[n] = p[0], q[m] = q[0];
    double tmp, res = 1e18;
    int u = 0, v = 0;
    for(int i = 0; i < n; i ++)
        if(p[i].y < p[u].y) u = i;
    for(int i = 0; i < m; i ++)
        if(q[i].y > q[v].y) v = i;
    for(int i = 0; i < n; i ++) {
        while(sgn(tmp = Cross(p[u], p[u+1], q[v+1]) - Cross(p[u], p[u+1], q[v])) > 0)
            if(++ v >= m) v -= m;
        if(sgn(tmp) < 0) res = min(res, Point_to_Seg(q[v], Line(p[u], p[u+1])));
        else res = min(res, Seg_to_Seg(Line(p[u], p[u+1]), Line(q[v], q[v+1])));
        if(++ u >= n) u -= n;
    }
    return res;
}

标签：return,Point,int,double,sgn,Cross,2D,计算,几何
来源： https://blog.csdn.net/m0_53373033/article/details/120425562