[模板] 计算几何基础

基本元素

平面几何可以看做是由点、线、向量、多边形、圆等基础元素构成的，下面我们先来了解如何定义这些元素

点和向量

在平面直角坐标系中，点和向量都可以被表示为一个二元组 $(x, y)$ ，我们用两个变量就可以表示点和向量。

struct Dot{
    long double x, y;

    Dot(){}
    Dot(long double _x, long double _y){
        x = _x;
        y = _y;
    }
};

线

我们通常使用一个点 $p$ ，和一个方向向量 $v$ 来表示一条线

$l = p + λ v$

当 $0 \leq λ \leq 1$ 时，这是一条线段
当 $0 \leq λ \leq \infty$ 时，这是一条射线
当 $- \infty \leq λ \leq \infty$ 时，这是一条直线

struct Line{
    Dot p, v;

    Line(){}
    Line(Dot _p, Dot _v){
        p = _p;
        v = _v;
    }
};

圆

圆可以由一个点 $p$ 表示圆心，一个数 $r$ 表示半径来确定

struct Circle{
    Dot p;
    long double r;

    Circle(){}
    Circle(Dot _p, long double _r){
        p = _p;
        r = _r;
    }
    Circle(long double _x, long double _y, long double _r){
        p = Dot(_x, _y);
        r = _r;
    }
};

多边形

多边形可以由各个顶点的坐标来确定

struct Polygon{
    int n;
    Dot p[MAXN];
};

浮点误差

由于计算几何使用了 double 型浮点变量来储存结果，因此浮点误差不可避免，在进行大小比较时容易出现问题。

因此我们需要设立一个 eps 作为误差允许的范围，再进行比较

const long double eps = 1e-9;

int dcmp(long double x){
    if(fabs(x) < eps)
        return 0;
    return x < 0 ? -1 : 1; 
}

调用方法也很简单，当比较 $a$ 与 $b$ 的大小时，只需要调用 dcmp(a - b) 即可，小于，等于，大于对应的返回 $- 1$ ， $0$ 和 $1$

向量的基础运算

向量加法与减法

坐标直接相加或相减即可

Dot operator + (Dot a, Dot b){
    return Dot(a.x + b.x, a.y + b.y);
}

Dot operator - (Dot a, Dot b){
    return Dot(a.x - b.x, a.y - b.y);
}

数乘向量

坐标直接相乘即可

Dot operator * (Dot a, long double b){
    return Dot(a.x * b, a.y * b);
}

Dot operator / (Dot a, long double b){
    return Dot(a.x / b, a.y / b);
}

向量乘法

幅角相加，模长相乘

Dot operator * (Dot a, Dot b){
    return Dot(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

向量点乘

向量 $\vec{a}$ 点乘向量 $\vec{b}$ 相当于 $\vec{a}$ 在 $\vec{b}$ 上的投影长度乘上 $\vec{b}$ 的模长

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

向量模长

向量 $\vec{a}$ 的模长等于 $d o t (\vec{a}, \vec{a})$ 开二次方

long double len(Dot a){
    return sqrt(dot(a, a));
}

long double len2(Dot a){
    return dot(a, a);
}

向量投影长度

计算 $\vec{a}$ 在 $\vec{b}$ 上的投影长度

long double shade(Dot a, Dot b){
    return dot(a, b) / len(b);
}

调整一个向量的长度

将向量 $\vec{a}$ 的模长调整为 $b$ 相当于乘以 $b$ 是 $\vec{a}$ 模长的多少倍

Dot adjust(Dot a, long double b){
    return a * (b / len(a));
}

向量叉乘（叉积）

向量 $\vec{a}$ 叉乘向量 $\vec{b}$ 相当于 $\vec{a}$ 和 $\vec{b}$ 围成的平行四边形的有向面积。

$\vec{a}$ 与左手向的向量叉乘为正，与右手向的向量叉乘为负，与和自己相等的向量叉乘为 $0$

long double cro(Dot a, Dot b){
    return a.x * b.y - a.y * b.x;
}

判断左手向与右手向

利用叉乘的性质， $\vec{b}$ 在 $\vec{a}$ 左手向返回 $0$ ，右手向返回 $1$

bool direct(Dot a, Dot b){
    return dcmp(cro(a, b)) <= 0;
}

垂直向量

求右手向垂直向量（叉积为正，点乘为 $0$ ）

Dot perpendicular(Dot a){
    return Dot(a.y, -a.x);
}

向量夹角

平行四边形的面积可由边 $a$ 乘边 $b$ 乘它们的夹角 $c o s (α)$ 得到，求角度移项即可

long double angle(Dot a, Dot b){
    return acos(dot(a, b) / len(a) / len(b));
}

向量旋转

还记得向量乘法的 幅角相加，模长相乘 吗

要将向量 $\vec{a}$ 旋转 $α$ 角，只需要乘上一个模长为 $1$ 幅角等于 $α$ 的向量

Dot angle(long double a){
    return Dot(cos(a), sin(a));
}

Dot rotate(Dot a, long double b){
    return a * angle(b);
}

判断线段相交

快速排斥试验

性质一 若线段 $A B$ 和线段 $C D$ 相交，则由他们的两个端点构成的矩形一定相交。

如果两条线段不满足这个性质，那么我们可以判断这两条直线不相交，如果满足，则进行跨立实验。

跨立试验

为什么要进行跨立实验？我们来看一种特殊的情况。

这样的线段可以通过快速排斥试验，因此我们还需要再对这种情况进行处理。

性质二 若两条直线相交，一条直线的两个端点一定在另一条直线的两侧。

这个我们可以通过之前的叉积来判断左手向和右手向，非常简单。

bool isIntersect(Line A, Line B){
    Dot a = A.p;
    Dot b = A.p + A.v;

    Dot c = B.p;
    Dot d = B.p + B.v;

    if(!(max(a.x, b.x) >= min(c.x, d.x) && max(c.x, d.x) >= min(a.x, b.x) && max(a.y, b.y) >= min(c.y, d.y) && max(c.y, d.y) >= min(a.y, b.y)))
        return false;

    if(dcmp(cro(A.v, d-a) * cro(A.v, c-a)) <= 0 && dcmp(cro(B.v, a-c) * cro(B.v, b-c)) <= 0)
        return true;

    return false;
}

计算几何基本操作

计算两个点之间的距离

只需要求出两点间向量的模长即可

long double dist(Dot a, Dot b){
    return len(a - b);
}

计算点到直线的距离

计算向量 $\vec{A . v}$ 与点 $A . p$ 到点 $B$ 的向量构成的平行四边形的面积

$B$ 到直线 $A$ 的距离就是这个平行四边形的高，向量 $\vec{A . v}$ 就是这条高的底

long double dist(Line a, Dot b){
    return fabs(cro(b - a.p, a.v) / len(a.v));
}

计算多边形面积

只要将多边形的顶点两两叉积，最后除以 $2$ 即可

对于图中的多边形 $D E F G$ ，就相当于分别计算 $O D E$ ， $O E F$ ， $O F G$ ， $O G D$ 的面积，因为向量叉积求得有向的，因此不在多边形内的面积会刚好被抵消（比如 $O G D$ 的面积抵消了之前三个三角形中计算的多余面积）

又因为叉积求得的面积其实是平行四边形面积，因此要除以 $2$

long double area(Polygon a){
    long double ans = 0;

    a.p[a.n] = a.p[0];

    for(register int i=0; i<a.n; i++){
        ans += cro(a.p[i], a.p[i+1]);
    }

    return fabs(ans / 2);
}

计算两直线交点

要求直线 $A$ 与直线 $B$ 的交点 $C$ ，只需要将点 $B . v$ 沿向量 $\vec{B . v}$ 方向平移一段距离 $x$ 即可

平行向量大小相等，因此可以将向量 $\vec{A . v}$ 平移到以 $B . p$ 为起点

我们构造两个平行四边形，一个是由向量 $\vec{A . v}$ 和向量 $\vec{B . v}$ 构成的粉红色四边形，一个是由向量 $\vec{A . v}$ 和从 $B . p$ 到 $A . p$ 的向量构成的紫色三角形。

这两个平行四边形共底（ $A . v$ ），则他们的面积之比就是高之比，因为两个三角形相似，所以高之比就是要平移的距离 $x$ 与 $\vec{B . v}$ 的模长之比

Dot intersect(Line a, Line b){
    return b.p + b.v * (cro(a.v, a.p - b.p) / cro(a.v, b.v));
}

线的平移

将线 $a$ 向右手向平移 $b$ 的单位长度

只需要构造一个与 $\vec{A . v}$ 垂直的向量，并将这个向量的长度调整至 $b$ ，然后将 $A . p$ 沿这个向量平移到 $B . p$ ，再与向量 $\vec{A . v}$ 重新组成直线即可

Line translate(Line a, long double b){
    return Line(a.p + adjust(perpendicular(a.v), b), a.v);
}

点的对称

将点 $a$ 沿直线 $b$ 作对称变换

先构造经过点 $A$ 且与 $\vec{B . v}$ 垂直的向量（对应着直线 $\vec{A A^{'}}$ ），求出他们的交点，再通过中点坐标公式即可得到结果。

Dot mirror(Line a, Dot b){
    return (intersect(a, Line(b, perpendicular(a.v))) * (long double)2) - b;
}

线的反射

计算直线 $A$ 在直线 $B$ 上反射后的直线

先计算出直线 $A$ 与直线 $B$ 的交点 $O$ ，再构造一条过点 $O$ 且方向向量与 $\vec{B . v}$ 垂直的直线，将 $A . p$ 沿该直线对称得到 $A . p^{'}$ ，过 $O$ 与 $A . p^{'}$ 的直线即为所求

Line reflect(Line a, Line b){
    Dot c = intersect(a, b);
    Dot d = mirror(Line(c, perpendicular(b.v)), a.p);

    return Line(c, d-c);
}

圆的交点

计算圆 $A$ 和圆 $B$ 的两个交点

第一步，先要判断圆 $A$ 和圆 $B$ 是否相交（不能是内含或是外离）
第二步，通过余弦定理计算出角 $C A B$ 的大小，将向量 $\vec{A B . v}$ 旋转，并调整长度为圆 $A$ 的半径长度，可以得到直线 $A C$ 与直线 $A D$

bool circle_intersect(Circle a, Circle b, Dot &c, Dot &d){
    Dot q = a.p - b.p;

    long double length_2 = len2(q);
    long double length = len(q);

    if(dcmp(length - a.r - b.r) > 0 || dcmp(length - fabs(a.r - b.r)) < 0)
        return 0;

    long double aph = acos((a.r * a.r + length_2 - b.r * b.r) / (a.r * length * 2));

    c = a.p + adjust(rotate(q, aph), a.r);
    d = a.p + adjust(rotate(q, -aph), a.r);

    return 1;
}