计算几何学习笔记

SUNCHAOYI · 2024-10-30 00:22:38 · 算法·理论

跟着这个题库完善自己的模板。

三角库函数

sin(x)/cos(x)/tan(x) 三角函数

asin(x)/acos(x)/atan(x) 反三角函数

atan2(x,y) 返回点 (x,y) 的反正切值，以弧度为单位。返回值的范围是 [-\pi,\pi]（扩展到四个象限），与 atan(x) 的具体区别如下：

每种函数关于精度问题有三种写法，以 atan2(x,y) 为例：

double atan2 (double y, double x);
float atan2f (float y, float x);
long double atan2l (long double y, long double x);

acos (-1.0) 计算 \pi 的值。

极角排序的应用

极角

一般以原点为中心，x 正半轴为始边，逆时针转过的角度（\theta \in [0,2\pi]）。
极角排序（第三象限 < 第四象限 < 第一象限 < 第二象限）

利用 atan(x,y) 直接排序（精度较低）

for (int i = 1;i <= n;++i) a[i] = atan2l (y[i],x[i]);
sort (a + 1,a + 1 + n);

利用叉积排序

根据 A,B,O 围成的面积（存在正负）即可判断。若面积为正，则 B 在 A 的正方向（当然，角度小于 \pi）；反之亦然。

bool cmp (Point &x, Point &y) 
{
  Point p({0, 0}); // 极点 
  Point v1 = a - p,v2 = b - p;
  if (v1 ^ v2 == 0) return a.x < b.x;//共线 
  else return v1 ^ v2 > 0;
}

先按照象限区分，然后按照极角排序【需要包含坐标轴】

int Qua (Point a)
{
if (a.x > 0 && a.y >= 0) return 1;
if (a.x <= 0 && a.y > 0) return 2;
if (a.x < 0 && a.y <= 0) return 3;
if (a.x >= 0 && a.y < 0) return 4;
}
bool cmp (Point a,Point b)
{
if (Qua (a) == Qub (b)) return atan2l (a.y,a.x) < atan2l (b.y,b.x);
else return Qua (a) < Qua (b);
}

【例题】 [ICPC2021 Macao R] Laser Trap

本质是三点能不能形成一个包围住原点的三角形。那么就可以先极角排序，然后用双指针，如果两点的极点差值小于 \pi，那么就不会经过原点，可以保留。为了方便书写代码，可以破环为链。同时，这道题需要注意精度问题。核心代码如下：

for (int i = 1;i <= n;++i) a[i + n] = a[i] + 2 * pi;//pi = acosl(-1);
for (int i = 1;i <= n;++i)
{
    while (p <= 2 * n && a[p] - a[i] < pi) ++p;
    ans = min (ans,p - 1 - i);
}

Computational Geometry(2D)

using LD = long double;
const LD pi = acos (-1.0);
const LD eps = 1e-8;
int dcmp (LD x) {return x < -eps ? -1 : (x > eps ? 1 : 0);}
struct Point {LD x,y;Point (LD x = 0,LD y = 0) : x (x),y (y) {}};
struct Circle {Point O;LD r;Circle (Point O = Point (),LD r = 0) : O (O),r (r) {}};
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,LD k) {return Vector (A.x * k,A.y * k);} 
Vector operator / (Vector A,LD k) {return Vector (A.x / k,A.y / k);} 

LD dot (Vector A,Vector B) {return A.x * B.x + A.y * B.y;} 
LD dis (Point A,Point B) {return sqrt ((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));} 
LD cross (Vector A,Vector B) {return A.x * B.y - A.y * B.x;} // A -> B counter-clockwise if cross (A,B) > 0
LD len (Point A)  {return sqrt (A.x * A.x + A.y * A.y);}
LD angle (Vector A,Vector B) {return acos (dot (A,B) / (len (A) * len (B)));}
Vector proj (Vector A,Vector B) {return A * (dot (A,B) / dot (A,A));} //project onto A
Point foot (Point P,Point A,Point B) {Vector AP = P - A,AB = B - A;return A + proj (AB,AP);} //foot
Point reflect (Point P,Point A,Point B) {Point F = foot (P,A,B);return F * 2 - P;} //symmetry point
Point rotate (Point P,LD theta) {return (Point){P.x * cos (theta) - P.y * sin (theta),P.x * sin (theta) + P.y * cos (theta)};}
bool on_line (Point P,Point A,Point B) {return dcmp (cross (P - A,B - A)) == 0;}
bool on_seg (Point P,Point A,Point B) {return on_line (P,A,B) && dcmp (dot (P - A,P - B)) <= 0;} //judge whether on segment AB
LD dis_seg (Point P,Point A,Point B)
{
    if (dcmp (dot (B - A,P - A)) < 0) return dis (P,A);
    if (dcmp (dot (A - B,P - B)) < 0) return dis (P,B);
    return fabs (cross (P - A,P - B)) / dis (A,B);
}
Point inter_line (Point A,Point B,Point C,Point D) {return A + (B - A) * cross (C - A,D - C) / cross (B - A,D - C);}
bool pd_ll_inter (Point A,Point B,Point C,Point D) {return dcmp (cross (B - A,D - C)) != 0;} // line - line
bool pd_ls_inter (Point A,Point B,Point C,Point D) {return on_line (inter_line (A,B,C,D),C,D);} //The intersection of AB(line) and CD (line) is on the CD (seg).
bool pd_ss_inter (Point A,Point B,Point C,Point D) // seg - seg
{
    LD c1 = cross (B - A,C - A),c2 = cross (B - A,D - A);
    LD d1 = cross (D - C,A - C),d2 = cross (D - C,B - C);
    if (dcmp (c1) * dcmp (c2) < 0 && dcmp (d1) * dcmp (d2) < 0) return true;
    if (dcmp(c1) == 0 && on_seg (C,A,B)) return true;
    if (dcmp(c2) == 0 && on_seg (D,A,B)) return true;
    if (dcmp(d1) == 0 && on_seg (A,C,D)) return true;
    if (dcmp(d2) == 0 && on_seg (B,C,D)) return true;
    return false;
}

LD area (vector <Point> P)
{
    int n = P.size ();
    LD res = 0;
    for (int i = 0;i < n;++i) res += cross (P[i],P[(i + 1) % n]);
    return res / 2.0;
}
bool is_convex (vector <Point> P)
{
    int n = P.size ();
    for(int i = 0;i < n - 1;++i) 
        if (dcmp (cross (P[i + 1] - P[i],P[(i + 2) % n] - P[i])) < 0) return false;
    return true;
}
int in_Poly (vector <Point> P,Point A)
{
    int cnt = 0,n = P.size ();
    for (int i = 0;i < n;++i)
    {
        int j = (i + 1) % n;
        if (on_seg (A,P[i],P[j])) return 2;// on the edge
        if (A.y >= min (P[i].y,P[j].y) && A.y < max (P[i].y,P[j].y)) // the intersection is on the right
            cnt += dcmp (((A.y - P[i].y) * (P[j].x - P[i].x) / (P[j].y - P[i].y) + P[i].x) - A.x) > 0;
    }
    return cnt & 1;
}
vector <Point> convex_hull (vector <Point> P) // strict convex hull (<= 0)
{
    int n = P.size ();
    sort (P.begin (),P.end (),[] (Point &x,Point &y) {return x.x == y.x ? x.y < y.y : x.x < y.x;});
     vector <Point> hull;
    hull.resize (2 * n + 1);
    int k = 0;
    for (int i = 0;i < n;++i) 
    {
        while (k >= 2 && dcmp (cross (hull[k - 1] - hull[k - 2],P[i] - hull[k - 2])) <= 0) --k;
        hull[k++] = P[i];
    }
    for (int i = n - 2,t = k;i >= 0;--i) 
    {
        while (k > t && dcmp (cross (hull[k - 1] - hull[k - 2],P[i] - hull[k - 2])) <= 0) --k;
        hull[k++] = P[i];
    }
    hull.resize (k - 1);
    return hull;
}
LD diameter (vector <Point> P)
{
    int n = P.size ();
    if (n <= 1) return 0;
    if (n == 2) return len (P[1] - P[0]);
    LD res = 0;
    for (int i = 0,j = 2;i < n;++i)
    {
        while (dcmp (cross (P[(i + 1) % n] - P[i],P[j] - P[i]) - cross (P[(i + 1) % n] - P[i],P[(j + 1) % n] - P[i])) <= 0) j = (j + 1) % n;
        res = max (res,max (len (P[j] - P[i]),len (P[j] - P[(i + 1) % n])));
    }
    return res;
}

bool in_cir (Circle C,Point P) {return dcmp (len (P - C.O) - C.r) <= 0;}
Point get_cir_p (Circle C,LD theta) {return {C.O.x + C.r * cos (theta),C.O.y + C.r * sin (theta)};}
int pd_lc_inter (Point A,Point B,Circle C) 
{
    double d = dis_seg (C.O,A,B);
    if (dcmp (d - C.r) == 0) return 0; // tangent
    if (dcmp (d - C.r) > 0) return -1; // separation
    return 1; // intersection
}
int pd_cc_inter (Circle A,Circle B) // the number of tagent lines
{
    LD d = len (A.O - B.O);
    if (dcmp (A.r + B.r - d) < 0) return 4; // externally separate
    if (dcmp (A.r + B.r - d) == 0) return 3; // externally tangent
    if (dcmp (fabs (A.r - B.r) - d) == 0) return 1; // internally tangent
    if (dcmp (fabs (A.r - B.r) - d) > 0) return 0; // one circle inside the other
    return 2; // intersection
}
pair <Point,Point> lc_inter (Point A,Point B,Circle C)
{
    Point F = foot (C.O,A,B);LD d = dis (C.O,F); 
    Vector E = (B - A) / dis (A,B);
    Point P1 = F - E * sqrt (C.r * C.r - d * d);
    Point P2 = F + E * sqrt (C.r * C.r - d * d);
    return {P1,P2};
}
pair <Point,Point> cc_inter (Circle A,Circle B)
{
    Vector k = B.O - A.O;
    LD d = len (k);
    LD alpha = atan2 (k.y,k.x),beta = acos ((A.r * A.r + d * d - B.r * B.r) / (2 * A.r * d));
    Point P1 = get_cir_p (A,alpha - beta),P2 = get_cir_p (A,alpha + beta);
    return {P1,P2};
}
pair <Point, Point> tan_cir (Point P,Circle C)
{
    LD d = len (C.O - P),theta = asin (C.r / d);
    Vector E = (C.O - P) / d;
    Vector P1 = P + (rotate (E,theta) * sqrt (d * d - C.r * C.r));
    Vector P2 = P + (rotate (E,-theta) * sqrt (d * d - C.r * C.r));
    return {P1,P2};
}
Circle triangle_incir (Point A,Point B,Point C) 
{
    LD a = dis (B,C),b = dis (A,C),c = dis (A,B);
    Point O = (A * a + B * b + C * c) / (a + b + c);
    return {O,dis_seg (O,A,B)};
}
Circle triangle_circum (Point A,Point B,Point C) 
{
    LD Bx = B.x - A.x,By = B.y - A.y,Cx = C.x - A.x,Cy = C.y - A.y;
    LD D = 2 * (Bx * Cy - By * Cx);
    LD x = (Cy * (Bx * Bx + By * By) - By * (Cx * Cx + Cy * Cy)) / D + A.x;
    LD y = (Bx * (Cx * Cx + Cy * Cy) - Cx * (Bx * Bx + By * By)) / D + A.y;
    Point P (x,y);
    return Circle (P,dis (A,P));
}
vector <pair <Point,Point>> get_tangents (Circle A,Circle B)
{
    vector <pair <Point,Point>> tangents;
    LD d = len (A.O - B.O),dif = A.r - B.r,sum = A.r + B.r;
    if (dcmp (d - fabs (dif)) < 0) return tangents;
    LD base = atan2 (B.O.y - A.O.y,B.O.x - A.O.x);
    if (dcmp (d - fabs (dif)) == 0) 
    {
        tangents.push_back ({get_cir_p (A,base + (A.r < B.r ? pi : 0)),get_cir_p (A,base + (A.r < B.r ? pi : 0))});
        return tangents;
    }
    LD theta = acos (dif / d);
    tangents.push_back ({get_cir_p (A,base + theta),get_cir_p (B,base + theta)});
    tangents.push_back ({get_cir_p (A,base - theta),get_cir_p (B,base - theta)});
    if (dcmp (d - sum) == 0) tangents.push_back ({get_cir_p (A,base),get_cir_p (A,base)});
    if (dcmp (d - sum) > 0)
    {
        theta = acos (sum / d);
        tangents.push_back ({get_cir_p (A,base + theta),get_cir_p (B,base + theta + pi)});
        tangents.push_back ({get_cir_p (A,base - theta),get_cir_p (B,base - theta + pi)});
    }
    return tangents;
}
LD tri_ploy_area (Point A,Point B,Circle C)
{
    Vector OA = A - C.O,OB = B - C.O;
    LD S = cross (OA,OB),sign = dcmp (cross (OA,OB)) > 0 ? 1 : -1;
    bool da = dcmp (len (OA) - C.r) < 0,db = dcmp (len (OB) - C.r) < 0;
    if (dcmp (S) == 0) return 0;
    if (da && db) return S * 0.5; // triangle
    if (!da && !db) 
    {
        if (pd_lc_inter (A,B,C) == 1)// arc + triangle + arc
        {
            auto [P1,P2] = lc_inter (A,B,C);
            Vector OP1 = P1 - C.O,OP2 = P2 - C.O;
            if (dis (A,P1) > dis (A,P2)) swap (P1,P2);
            return cross (OP1,OP2) * 0.5 + sign * 0.5 * C.r * C.r * (angle (OA,OP1) + angle (OB,OP2));
        }
        else return sign * 0.5 * C.r * C.r * angle (OA,OB); // arc
    }
    else // triangle + arc
    {
        auto [P1,P2] = lc_inter (A,B,C);
        if (on_seg (P2,A,B)) swap (P1,P2); 
        Vector OP1 = P1 - C.O;
        if (dcmp (len (OA) - C.r) < 0) return cross (OA,OP1) * 0.5 + sign * 0.5 * C.r * C.r * angle (OP1,OB);
        else return cross (OP1,OB) * 0.5 + sign * 0.5 * C.r * C.r * angle (OP1,OA);
    }
}
LD cc_area (Circle C1,Circle C2)
{
    int op = pd_cc_inter (C1,C2);
    if (op <= 1) return pi * min (C1.r,C2.r) * min (C1.r,C2.r);
    else if (op == 4) return 0;
    else
    {
        LD d = dis (C1.O,C2.O);
        LD alpha = 2 * acos ((C1.r * C1.r - C2.r * C2.r + d * d) / (2 * C1.r * d));
        LD beta = 2 * acos ((C2.r * C2.r - C1.r * C1.r + d * d) / (2 * C2.r * d));
        return 0.5 * (C1.r * C1.r * (alpha - sin (alpha)) + C2.r * C2.r * (beta - sin (beta)));
    }
}

计算几何 学习笔记

三角库函数

极角排序的应用

计算几何学习笔记