简单平面几何模板(持续更新)

 1. 向量基本操作

struct Point { double x, y;
    Point(){}
    Point(double a, double b): x(a), y(b){}
};

typedef Point Vector;

Vector operator + (Vector a, Vector b) { return Vector( a.x + b.x, a.y + b.y ); }
Vector operator - (Vector a, Vector b) { return Vector( a.x - b.x, a.y - b.y ); }
Vector operator * (Vector a, double p) { return Vector( a.x * p, a.y * p); }
Vector operator / (Vector a, double p) { return Vector( a.x / p, a.y / p); }

int dcmp(double x) { if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; } // 精准判断正负号

// 点积,模长,两向量夹角, 向量角度,叉积,三角形面积,极角,逆时针旋转
double Dot(Vector a, Vector b) { return a.x * b.x + a.y * b.y; }
double Length(Vector a) { return sqrt( Dot(a, a) ); }
double Angle(Vector a, Vector b) { return acos(Dot(a, b) / Length(a) / Length(b)); }
double Angle(Vector a) { return atan2(a.y, a.x);}
double Cross(Vector a, Vector b) { return a.x * b.y - a.y * b.x; } // 大于0时,b在a的逆时针方向
double Area(Point A, Point B, Point C) { return Cross(B - A, C - A) / 2; } // 三角形面积叉乘公式
double PolarAngle(Vector a) { return atan2(a.y, a.x); }
Vector Rotate (Vector a, double rad) { return Vector( a.x*cos(rad) - a.y*sin(rad), a.x*sin(rad) + a.y*cos(rad) ); }

bool operator < (Vector a, Vector b) {
    if(dcmp(PolarAngle(a) - PolarAngle(b)) == 0)
        return dcmp(a.x - b.x) == 0 ? dcmp(a.y - b.y) < 0 : dcmp(a.x - b.x) < 0;
    return dcmp(PolarAngle(a) - PolarAngle(b)) < 0;
} //优先级:极角 > x > y
bool operator == (Vector a, Vector b) { return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0; }


struct Line{ // 有向直线
    Point A, B;
    Line(){};
    Line(Point P, Point Q) : A(P), B(Q) {}
};

关于旋转的解释(r为模长)

2.判断线段相交 hdu 1086

// 判断线段AB是否跨立MN
bool straddle (Point A, Point B, Point M, Point N) {
    Vector p0 = N - M; // 向量 MN
    Vector p1 = A - M; // 向量 MA
    Vector p2 = B - M; // 向量 MB
    return Cross(p1, p0) * Cross(p2, p0) <= 0;
}

// 判断 线段AB为对角线构成的矩形 和 MN线段为对角线构成的矩形 是否存在重叠
bool overlap (Point A, Point B, Point M, Point N){
    Point P = { max(min(A.x, B.x), min(M.x, N.x)), max(min(A.y, B.y), min(M.y, N.y))};
    Point Q = { min(max(A.x, B.x), max(M.x, N.x)), min(max(A.y, B.y), max(M.y, N.y))};
    return P.x <= Q.x && P.y <= Q.y;
}

// 判断线段AB和MN是否相交
bool intersect (Point A, Point B, Point M, Point N) {
    return straddle(A, B, M, N) && straddle(M, N, A, B) && overlap(A, B, M, N);
}

3. 直线交点 POJ1269

// 方法1,先确定直线 A+v*t1 和 C+w*t2 有唯一交点, 及Cross(v, w) != 0
Point GLI(Point P, Vector v, Point Q, Vector w){ // 两直线的交点
    Vector u = P - Q;
    double t1 = Cross(w, u) / Cross(v, w);
    return P + v * t1;
}

//方法2 两直线的交点 前提Cross(v, w) != 0
Point GLI(Line l1, Line l2){
    Vector v = l1.B - l1.A;
    Vector w = l2.B - l2.A;
    Vector u = l1.A - l2.A;
    double t1 = Cross(w, u) / Cross(v, w);
    return l1.A + v * t1;
}

说明(该问题困扰了我很久,请教的理学院的某大佬):

4.点到直线的距离和点到线段的距离

double DTL(Point P, Point A, Point B){ // 点到直线的距离
    Vector v1 = B - A, v2 = P - A;
    return fabs(Cross(v1, v2) / Length(v1));
}

double DTS(Point P, Point A, Point B){  // 点到线段的距离
    if(A == B) return Length(P-A);
    Vector v1 = B - A, v2 = P - A, v3 = P - B;
    if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
    else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
    else return fabs(Cross(v1, v2)) / Length(v1);
}

5.多边形有向面积

double PolygonArea(Point *P, int n){ // 多边形有向面积
    double area = 0;
    for(int i=1; i<n-1; i++) area += Cross(P[i] - P[0], P[i+1] - P[0]);
    return area / 2;
}

6.凸包Andrew算法, POJ 1113

bool com(Point A, Point B){ return (fabs(A.x - B.x) < eps && A.y < B.y) || A.x < B.x; }

int ConvexHull(Point *P, int n, Point *A){ // P为点集合, A为凸包点(不包括边上的点)
    sort(P, P + n, com);
    int m = 0;
    for(int i=0; i<n; i++) {
        while(m > 1 && Cross(A[m-1] - A[m-2], P[i] - A[m-2]) <= 0) m --;
        A[m++] = P[i];
    }
    int k = m;
    for(int i=n-2; i>=0; i--){
        while(m > k && Cross(A[m-1] - A[m-2], P[i] - A[m-2]) <= 0) m --;
        A[m++] = P[i];
    }
    if(n > 1) m -- ;
    return m;
}

7.点在多边形内的判定

方法1:
int isPIP(Point p, Polygon poly) {
    int wn = 0;
    int n = poly.size();
    for(int i = 0; i < n; i++) {
        if(OnSegment(p, poly[i], poly[(i+1)%n])) return -1; //在边界上
        int k = dcmp(Cross(poly[(i+1)%n]-poly[i], p-poly[i]));
        int d1 = dcmp(poly[i].y - p.y);
        int d2 = dcmp(poly[(i+1)%n].y - p.y);
        if(k > 0 && d1 <= 0 && d2 > 0) wn++;
        if(k < 0 && d2 <= 0 && d1 > 0) wn--;
    }
    if(wn != 0) return 1;       //内部
    return 0;                   //外部
}



double DTS(Point P, Point A, Point B){  // 点到线段的距离
    if(A == B) return Length(P-A);
    Vector v1 = B - A, v2 = P - A, v3 = P - B;
    if(dcmp(Dot(v1, v2)) < 0) return Length(v2);
    else if(dcmp(Dot(v1, v3)) > 0) return Length(v3);
    else return fabs(Cross(v1, v2)) / Length(v1);
}

方法2:
int isPIP(Point p, Point *poly, int n) {
    int wn = 0;
    for(int i = 0; i < n; i++) {
        if(dcmp(DTS(p, poly[i], poly[(i+1)%n])) == 0) return -1; //在边界上
        int k = dcmp(Cross(poly[(i+1)%n]-poly[i], p-poly[i]));
        int d1 = dcmp(poly[i].y - p.y);
        int d2 = dcmp(poly[(i+1)%n].y - p.y);
        if(k > 0 && d1 <= 0 && d2 > 0) wn++;
        if(k < 0 && d2 <= 0 && d1 > 0) wn--;
    }
    if(wn != 0) return 1;       //内部
    return 0;                   //外部
}

8. 圆:求线与圆交点个数及交点

struct Point { double x, y;
    Point(){}
    Point(double a, double b): x(a), y(b){}
};

struct Line {
    Point A, B;
    Line(){};
    Line(Point P, Point Q) : A(P), B(Q) {}
};

struct Circle{
    Point c; double r;
    Point point(double a){ return Point(c.x + cos(a) * r, c.y + sin(a) * r); }
};


// 求线与圆交点个数及交点
// 1. (a*t + b)^2 + (c*t + d)^2 = r^2
// 2. e*t^2 + f*t + g = 0
// 每次交点单独计算需清空sol
int getLCI(Line L, Circle C, double& t1, double& t2, vector<Point>& sol) {
    double a = L.B.x - L.A.x, b = L.A.x - C.c.x, c = L.B.y - L.A.y, d = L.A.y - C.c.y;
    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 (dcmp(delta) < 0) return 0; //相离
    if (dcmp(delta) == 0) {        //相切
        t1 = t2 = -f / (2 * e);
        sol.push_back(C.point(t1));
        return 1;
    }
    t1 = (-f - sqrt(delta)) / (2 * e); sol.push_back(C.point(t1));
    t2 = (-f + sqrt(delta)) / (2 * e); sol.push_back(C.point(t2));
    return 2; // 相交
}

9.圆:求圆和圆的交点个数和交点

// 求圆和圆的交点个数和交点
int getCCI(Circle C1, Circle C2, vector<Point>& sol) {
    double d = Length(C1.c - C2.c);
    if (dcmp(d) == 0) {
        if (dcmp(C1.r - C2.r == 0)) return -1;    //两圆完全重合
        return 0;                                //同心圆,半径不一样
    }
    if (dcmp(C1.r + C2.r - d) < 0) return 0;
    if (dcmp(fabs(C1.r - C2.r) - d) > 0) return 0;

    double a = Angle(C2.c - C1.c);               //向量C1C2的极角
    double da = acos((C1.r * C1.r + d * d - C2.r * C2.r) / (2 * C1.r * d));
    //C1C2到C1P1的角
    Point p1 = C1.point(a - da), p2 = C1.point(a + da);
    sol.push_back(p1);
    if (p1 == p2) return 1;
    sol.push_back(p2);
    return 2;
}

10.圆:过定点做圆的切线

//过定点做圆的切线
//过点p做圆C的切线,返回切线个数。v[i]表示第i条切线
int getTangents(Point p, Circle C, Vector* v) {
    Vector u = C.c - p;
    double dist = Length(u);
    if(dist < C.r) return 0;
    else if(dcmp(dist - C.r) == 0) {
        v[0] = Rotate(u, PI/2);
        return 1;
    } else {
        double ang = asin(C.r / dist);
        v[0] = Rotate(u, -ang);
        v[1] = Rotate(u, +ang);
        return 2;
    }
}

11.圆:两圆的公切线

//两圆的公切线
//返回切线的个数,-1表示有无数条公切线。
//a[i], b[i] 表示第i条切线在圆A,圆B上的切点
int getTangents(Circle A, Circle B, Point *a, Point *b) {
    int cnt = 0;
    if(A.r < B.r) {
        swap(A, B); swap(a, b);
    }
    int d2 = (A.c.x - B.c.x)*(A.c.x - B.c.x) + (A.c.y - B.c.y)*(A.c.y - B.c.y);
    int rdiff = A.r - B.r;
    int rsum = A.r + B.r;
    if(d2 < rdiff*rdiff) return 0;   //内含
    double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
    if(d2 == 0 && A.r == B.r) return -1;   //无限多条切线
    if(d2 == rdiff*rdiff) {         //内切一条切线
        a[cnt] = A.point(base);
        b[cnt] = B.point(base);
        cnt++;
        return 1;
    }
    //有外共切线
    double ang = acos((A.r-B.r) / sqrt(d2));
    a[cnt] = A.point(base+ang); b[cnt] = B.point(base+ang); cnt++;
    a[cnt] = A.point(base-ang); b[cnt] = B.point(base-ang); cnt++;
    if(d2 == rsum*rsum) {  //一条公切线
        a[cnt] = A.point(base);
        b[cnt] = B.point(PI+base);
        cnt++;
    } else if(d2 > rsum*rsum) {   //两条公切线
        double ang = acos((A.r + B.r) / sqrt(d2));
        a[cnt] = A.point(base+ang); b[cnt] = B.point(PI+base+ang); cnt++;
        a[cnt] = A.point(base-ang); b[cnt] = B.point(PI+base-ang); cnt++;
    }
    return cnt;
}

12.将多边形转为逆时针顺序,如果有向面积为正数则为逆时针

// 多边形有向面积
double PolygonArea(Point *P, int n){
    double area = 0;
    for(int i=1; i<n-1; i++) area += Cross(P[i] - P[0], P[i+1] - P[0]);
    return area / 2;
}

// 将多边形转为逆时针顺序,如果有向面积为正数则为逆时针
void Reverse(Point *pol, int n) {
    double ans = PolygonArea(pol, n);
    if(ans < 0) for(int i=0; i<n/2; i++) swap(pol[i], pol[n-i-1]);
}

13.模板整理

typedef pair<int,int> pii;
const double pi = 4 * atan(1);

inline int dcmp (double x) { if (fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; }
inline double getDistance (double x, double y) { return sqrt(x * x + y * y); }
inline double torad(double deg) { return deg / 180 * pi; }

struct Point {
    double x, y;
    Point (double x = 0, double y = 0): x(x), y(y) {}
    void read () { scanf("%lf%lf", &x, &y); }
    void write () { printf("%lf %lf", x, y); }

    bool operator == (const Point& u) const { return dcmp(x - u.x) == 0 && dcmp(y - u.y) == 0; }
    bool operator != (const Point& u) const { return !(*this == u); }
    bool operator < (const Point& u) const { return x < u.x || (x == u.x && y < u.y); }
    bool operator > (const Point& u) const { return u < *this; }
    bool operator <= (const Point& u) const { return *this < u || *this == u; }
    bool operator >= (const Point& u) const { return *this > u || *this == u; }
    Point operator + (const Point& u) { return Point(x + u.x, y + u.y); }
    Point operator - (const Point& u) { return Point(x - u.x, y - u.y); }
    Point operator * (const double u) { return Point(x * u, y * u); }
    Point operator / (const double u) { return Point(x / u, y / u); }
    double operator * (const Point& u) { return x*u.y - y*u.x; }
};
typedef Point Vector;
typedef vector<Point> Polygon;

struct Line {
    double a, b, c;
    Line (double a = 0, double b = 0, double c = 0): a(a), b(b), c(c) {}
};

struct DirLine {
    Point p;
    Vector v;
    double ang;
    DirLine () {}
    DirLine (Point p, Vector v): p(p), v(v) { ang = atan2(v.y, v.x); }
    bool operator < (const DirLine& u) const { return ang < u.ang; }
};


namespace Punctual {
    double getDistance (Point a, Point b) { double x=a.x-b.x, y=a.y-b.y; return sqrt(x*x + y*y); }
};


namespace Vectorial {
    /* 点积: 两向量长度的乘积再乘上它们夹角的余弦, 夹角大于90度时点积为负 */
    double getDot (Vector a, Vector b) { return a.x * b.x + a.y * b.y; }

    /* 叉积: 叉积等于两向量组成的三角形有向面积的两倍, cross(v, w) = -cross(w, v) */
    double getCross (Vector a, Vector b) { return a.x * b.y - a.y * b.x; }

    double getLength (Vector a) { return sqrt(getDot(a, a)); }
    double getPLength (Vector a) { return getDot(a, a); }
    double getAngle (Vector u) { return atan2(u.y, u.x); }
    double getAngle (Vector a, Vector b) { return acos(getDot(a, b) / getLength(a) / getLength(b)); }
    Vector rotate (Vector a, double rad) { return Vector(a.x*cos(rad)-a.y*sin(rad), a.x*sin(rad)+a.y*cos(rad)); }
    /* 单位法线 */
    Vector getNormal (Vector a) { double l = getLength(a); return Vector(-a.y/l, a.x/l); }
};

namespace ComplexVector {
    typedef complex<double> Point;
    typedef Point Vector;

    double getDot(Vector a, Vector b) { return real(conj(a)*b); }
    double getCross(Vector a, Vector b) { return imag(conj(a)*b); }
    Vector rotate(Vector a, double rad) { return a*exp(Point(0, rad)); }
};

namespace Linear {
    using namespace Vectorial;

    Line getLine (double x1, double y1, double x2, double y2) { return Line(y2-y1, x1-x2, y1*x2-x1*y2); }
    Line getLine (double a, double b, Point u) { return Line(a, -b, u.y * b - u.x * a); }

    bool getIntersection (Line p, Line q, Point& o) {
        if (fabs(p.a * q.b - q.a * p.b) < eps)
            return false;
        o.x = (q.c * p.b - p.c * q.b) / (p.a * q.b - q.a * p.b);
        o.y = (q.c * p.a - p.c * q.a) / (p.b * q.a - q.b * p.a);
        return true;
    }

    /* 直线pv和直线qw的交点 */
    bool getIntersection (Point p, Vector v, Point q, Vector w, Point& o) {
        if (dcmp(getCross(v, w)) == 0) return false;
        Vector u = p - q;
        double k = getCross(w, u) / getCross(v, w);
        o = p + v * k;
        return true;
    }

    /* 点p到直线ab的距离 */
    double getDistanceToLine (Point p, Point a, Point b) { return fabs(getCross(b-a, p-a) / getLength(b-a)); }
    double getDistanceToSegment (Point p, Point a, Point b) {
        if (a == b) return getLength(p-a);
        Vector v1 = b - a, v2 = p - a, v3 = p - b;
        if (dcmp(getDot(v1, v2)) < 0) return getLength(v2);
        else if (dcmp(getDot(v1, v3)) > 0) return getLength(v3);
        else return fabs(getCross(v1, v2) / getLength(v1));
    }

    /* 点p在直线ab上的投影 */
    Point getPointToLine (Point p, Point a, Point b) { Vector v = b-a; return a+v*(getDot(v, p-a) / getDot(v,v)); }

    /* 判断线段是否存在交点 */
    bool haveIntersection (Point a1, Point a2, Point b1, Point b2) {
        double c1=getCross(a2-a1, b1-a1), c2=getCross(a2-a1, b2-a1), c3=getCross(b2-b1, a1-b1), c4=getCross(b2-b1,a2-b1);
        return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4) < 0;
    }

    /* 判断点是否在线段上 */
    bool onSegment (Point p, Point a, Point b) { return dcmp(getCross(a-p, b-p)) == 0 && dcmp(getDot(a-p, b-p)) < 0; }
    bool onLeftOrLine(DirLine l, Point p) { return dcmp(l.v * (p - l.p)) >= 0; }
}

namespace Triangular {
    using namespace Vectorial;

    double getAngle (double a, double b, double c) { return acos((a*a+b*b-c*c) / (2*a*b)); }
    double getArea (double a, double b, double c) { double s =(a+b+c)/2; return sqrt(s*(s-a)*(s-b)*(s-c)); }
    double getArea (double a, double h) { return a * h / 2; }
    double getArea (Point a, Point b, Point c) { return fabs(getCross(b - a, c - a)) / 2; }
    double getDirArea (Point a, Point b, Point c) { return getCross(b - a, c - a) / 2; }
};

namespace Polygonal {
    using namespace Vectorial;
    using namespace Linear;

    double getArea (Point* p, int n) {
        double ret = 0;
        for (int i = 1; i < n-1; i++)
            ret += getCross(p[i]-p[0], p[i+1]-p[0]);
        return fabs(ret)/2;
    }

    /* 凸包 */
    int getConvexHull (Point* p, int n, Point* ch) {
        sort(p, p + n);
        int m = 0;
        for (int i = 0; i < n; i++) {
            /* 可共线 */
            //while (m > 1 && dcmp(getCross(ch[m-1]-ch[m-2], inp[i]-ch[m-1])) < 0) m--;
            while (m > 1 && dcmp(getCross(ch[m-1]-ch[m-2], p[i]-ch[m-1])) <= 0) m--;
            ch[m++] = p[i];
        }
        int k = m;
        for (int i = n-2; i >= 0; i--) {
            /* 可共线 */
            //while (m > k && dcmp(getCross(ch[m-1]-ch[m-2], inp[i]-ch[m-2])) < 0) m--;
            while (m > k && dcmp(getCross(ch[m-1]-ch[m-2], p[i]-ch[m-2])) <= 0) m--;
            ch[m++] = p[i];
        }
        if (n > 1) m--;
        return m;
    }

   int isPointInPolygon (Point o, Point* p, int n) {
        int wn = 0;
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            if (onSegment(o, p[i], p[j])) return -1; // 边界上
            int k = dcmp(getCross(p[j] - p[i], o-p[i]));
            int d1 = dcmp(p[i].y - o.y);
            int d2 = dcmp(p[j].y - o.y);
            if (k > 0 && d1 <= 0 && d2 > 0) wn++;
            if (k < 0 && d2 <= 0 && d1 > 0) wn--;
        }
        return wn ? 1 : 0; // 1在内,0在外
    }


    /* 旋转卡壳 */
    void rotatingCalipers(Point *p, int n, vector<pii>& sol) {
        sol.clear();
        int j = 1; p[n] = p[0];
        for(int i = 0; i < n; i++) {
            while (getCross(p[j+1]-p[i+1], p[i]-p[i+1]) > getCross(p[j]-p[i+1], p[i]-p[i+1]))
                j = (j+1) % n;
            sol.push_back(make_pair(i, j));
            sol.push_back(make_pair(i + 1, j + 1));
        }
    }

    /* 计算半平面相交可以用增量法,o(n^2),初始设置4条无穷大的半平面 */
    /* 用有向直线A->B切割多边形u,返回左侧。可能退化成单点或线段 */
    Polygon CutPolygon (Polygon u, Point a, Point b) {
        Polygon ret;
        int n = u.size();
        for (int i = 0; i < n; i++) {
            Point c = u[i], d = u[(i+1)%n];
            if (dcmp((b-a)*(c-a)) >= 0) ret.push_back(c);
            if (dcmp((b-a)*(c-d)) != 0) {
                Point t;
                getIntersection(a, b-a, c, d-c, t);
                if (onSegment(t, c, d))
                    ret.push_back(t);
            }
        }
        return ret;
    }

    bool hpicom(DirLine l1, DirLine l2) {
        if(l1.ang == l2.ang) return onLeftOrLine(l2, l1.p);
        return l1.ang < l2.ang;
    }

    Point * inp = new Point[maxn];
    DirLine * que = new DirLine[maxn];

    /* 半平面相交, 保证半平面交有限*/
    int halfPlaneIntersection(DirLine* li, int n, Point* poly) {
        sort(li, li + n, hpicom);
        int first, last;
        que[first= last=0] = li[0];
        for (int i = 1; i < n; i++) {
            if(que[last].ang == li[i].ang) continue;
            while (first < last && !onLeftOrLine(li[i], inp[last - 1])) last--;
            while (first < last && !onLeftOrLine(li[i], inp[first])) first++;
            que[++last] = li[i];
            if (first < last)
                getIntersection(que[last - 1].p, que[last - 1].v, que[last].p, que[last].v, inp[last - 1]);
        }

        while (first < last && !onLeftOrLine(que[first], inp[last - 1])) last--;
        if (last - first <= 1) return 0;
        getIntersection(que[last].p, que[last].v, que[first].p, que[first].v, inp[last]);

        int m = 0;
        for (int i = first; i <= last; i++) poly[m++] = inp[i];
        return m;
    }
};

using namespace Vectorial;
using namespace Polygonal;

14.将有向直线左移d距离

//向左平移d距离
Line Translation(Line l, double d) {
    Vector AB = l.B - l.A;
    Vector AP = Rotate(AB, PI / 2);
    Point C = l.A + AP * d / Length(AP);
    Point D = l.B + AP * d / Length(AP);
    return Line(C, D);
}

 

  • 3
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值