[模板] 计算几何基础

基本元素

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

点和向量

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

struct Dot{
    long double x, y;

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

线

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

$$l = p + \lambda v$$

  • 当 $0 \leq \lambda \leq 1$ 时,这是一条线段
  • 当 $0 \leq \lambda \leq \infty$ 时,这是一条射线
  • 当 $-\infty \leq \lambda \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$ 的模长等于 $dot(\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$ 乘 它们的夹角 $cos( \alpha)$ 得到,求角度移项即可

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

向量旋转

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

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

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

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

判断线段相交

快速排斥试验

性质一 若线段 $AB$ 和 线段 $CD$ 相交,则由他们的两个端点构成的矩形一定相交。



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



跨立试验

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



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

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

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

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);
} 

计算点到直线的距离



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

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

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

计算多边形面积



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

对于图中的多边形 $DEFG$,就相当于分别计算 $ODE$,$OEF$,$OFG$,$OGD$ 的面积,因为向量叉积求得有向的,因此不在多边形内的面积会刚好被抵消(比如 $OGD$ 的面积抵消了之前三个三角形中计算的多余面积)

又因为叉积求得的面积其实是平行四边形面积,因此要除以 $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$ 沿向量 $\overrightarrow{B.v}$ 方向平移一段距离 $x$ 即可

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

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

这两个平行四边形共底( $A.v$ ),则他们的面积之比就是高之比,因为两个三角形相似,所以高之比就是要平移的距离 $x$ 与 $\overrightarrow{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$ 的单位长度



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

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

点的对称

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



先构造经过点 $A$ 且与 $\overrightarrow{B.v}$ 垂直的向量(对应着直线 $\overrightarrow{AA’}$),求出他们的交点,再通过中点坐标公式即可得到结果。

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

线的反射

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



先计算出直线 $A$ 与 直线 $B$ 的交点 $O$,再构造一条过点 $O$ 且方向向量与 $\overrightarrow{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$ 是否相交(不能是内含或是外离)
第二步,通过余弦定理计算出角 $CAB$ 的大小,将向量 $\overrightarrow{AB.v}$ 旋转,并调整长度为圆 $A$ 的半径长度,可以得到直线 $AC$ 与 直线 $AD$

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;
} 

后记

以上是计算几何的一些最基本的操作,其他更高级的算法都是基于这些操作得到。

其他计算几何相关算法:

  • Pick 定理
  • 扫描线
  • 二维凸包
  • 旋转卡壳
  • 半平面交
  • KD-Tree

参考资料