其他分享
首页 > 其他分享> > 计算几何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