二维计算几何基础 学习笔记
观前提醒：烂尾 & 此文年代过于久远，不保证内容的准确性，有锅我是不会修的。
基操
浮点数
正常的计算几何题都是在 \(\R\) 上进行的，所以我们要进行一些 double 运算。普通的运算没什么要说的，重点在于比较两个浮点数的大小。由于精度误差的存在，我们经常需要认为例如 1 == 1.0000001 是 true。通常我们设一个精度范围 \(\epsilon\) 也就是 eps（通常在 \(10^{-4}\sim 10^{-10}\) 之间），然后定义以下符号函数 \(\operatorname{sgn}\)。
constexpr db eps = ?;
int sgn(db x) { return abs(x) <= eps ? 0 : x < 0 ? -1 : 1; }
这个函数可以让我们知道一个浮点数的符号是什么。基于此，我们可以实现各种比较了：
x == y即sgn(x - y) == 0，x != y即sgn(x - y) != 0。x < y即sgn(x - y) < 0，x > y即sgn(x - y) > 0。x <= y即sgn(x - y) <= 0，x >= y即sgn(x - y) >= 0。
顺带定义一个平方函数，省的有时要对一个长式子做平方，如果直接展开写的话长式子要算两遍，否则就要用一个临时变量记录，都称不上体面的写法。
db sq(db x) { return x * x; }
我们可以通过 acos(-1) 获取 \(\pi\) 的值。
constexpr db pi = acos(-1);
二维向量
记录
一般我们认为向量和点是同一回事，毕竟它们在平面直角坐标系中形成双射。直接用两个分量记录它们。
struct point {
db x, y;
point(db x = 0, db y = 0) : x(x), y(y) {}
bool operator==(point v) { return sgn(x - v.x) == 0 && sgn(y - v.y) == 0; }
bool operator!=(point v) { return sgn(x - v.x) || sgn(y - v.y); }
friend ostream &operator<<(ostream &cout, point u) { return cout << mt(u.x, u.y); }
};
线性运算
加法、减法、数乘。加法表示平移，数乘表示放缩。一个向量减去一个点，表示以这个点作为坐标系原点后的坐标。负号为相反向量。
struct point {
point operator+(point v) { return point(x + v.x, y + v.y); }
point operator-(point v) { return point(x - v.x, y - v.y); }
point operator-() { return point(-x, -y); }
point operator*(db v) { return point(x * v, y * v); }
friend point operator*(db u, point v) { return point(u * v.x, u * v.y); }
point operator/(db v) { return point(x / v, y / v); }
};
模长、极角
模长为一点到原点的距离，亦即对应向量的长度。公式 \(|\pmb v|=\sqrt{x^2+y^2}\)。显然，点 \(\pmb u,\pmb v\) 的距离就是 \(|\pmb u-\pmb v|\)。称 \(\dfrac{\pmb v}{|\pmb v|}\) 为 \(\pmb v\) 的单位向量，与 \(\pmb v\) 同向且模长为 \(1\)。
一个角度是 \(\bmod 2\pi\) 的一个剩余类。我们定义两个向量 \(\pmb u,\pmb v\)（都不为 \(\pmb 0\)）的夹角 \(\langle\pmb u,\pmb v\rangle\) 为 \(\pmb u\) 逆时针转多少能与 \(\pmb v\) 共线且同向。显然 \(\langle\pmb u,\pmb v\rangle=-\langle\pmb v,\pmb u\rangle\)。
定义非零向量 \(\pmb v\) 的极角 \(\arg\pmb v=\lang(1,0),\pmb v\rang\)。显然 \(\lang\pmb u,\pmb v\rang=\arg\pmb v-\arg\pmb u\)。如何求极角？容易想到求坐标比值然后利用反三角函数，但是没有任何一个反三角函数的值域是 \([0,2\pi)\) 的，需要我们自己分象限讨论，还要特判分母为 \(0\) 的情况，非常麻烦。幸运的是，C++ 专门提供了一个函数 atan2(y, x) 直接就求出了 \(\arg(x,y)\)。只不过值域是 \((-\pi,\pi]\)，可以通过把 \(<0\) 的加 \(2\pi\) 把值域调整到 \([0,2\pi)\)。
考虑将向量逆时针旋转 \(\theta\) 角。这是个线性变换，标准矩阵可以由 \((1,0)\) 和 \((0,1)\) 转到的向量拼起来得到，是 \(\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}\)。逆时针旋转 \(-\theta\) 等价于以 \(\theta\) 方向上的向量为新的横轴后的坐标。
struct point {
db operator~() { return sqrt(sq(x) + sq(y)); } // 模长
friend db dis(point u, point v) { return ~(u - v); } // 距离
point unit() { return *this / ~*this; } // 单位向量
db operator!() { db res = atan2(y, x); return res < 0 ? res + 2 * pi : res; } // 极角
friend db angle(point u, point v) { db res = !v - !u; return res < 0 ? res + 2 * pi : res; } // 夹角
point operator<<(db theta) { db c = cos(theta), s = sin(theta); return point(c * x - s * y, s * x + c * y); } // 旋转
};
内积 / 点积 / 数量积
跟据线性代数知识，我们知道 \(\pmb u\cdot\pmb v=|\pmb u||\pmb v|\cos\lang\pmb u,\pmb v\rang=x_{\pmb u}x_{\pmb v}+y_{\pmb u}y_{\pmb v}\)，有交换律、线性性，几何意义是 \(\pmb u\) 在 \(\pmb v\) 上的投影长度乘以 \(\pmb v\) 的模长。有了点积，我们就可以避开三角函数和分类讨论，使用小常数高精度做一些判定：
- 垂直判定：\(\pmb u\perp\pmb v\)（即 \(\pmb u\) 或 \(\pmb v\) 为 \(\pmb 0\) 或 \(\lang\pmb u,\pmb v\rang=\pm\dfrac{\pi}2\)）当且仅当 \(\pmb u\cdot\pmb v=0\)。
struct point {
db operator*(point v) { return x * v.x + y * v.y; } // 内积
friend bool perp(point u, point v) { return sgn(u * v) == 0; } // 垂直判定
};
叉积 / 向量积
这玩意怎么说呢……很多 OIer 仿佛都对它有一些误解。
最普适的定义是外积，定义在 \(\R^n\times\R^n\) 中，而取值空间是一个既不是 \(\R\) 也不是 \(\R^n\) 的空间，我也不知道是啥（因为太难了），就记为 \(\mathbb W_n\) 吧。而只有在 \(0,1,3,7\) 维（分别对应实数、复数、四元数、八元数）有 \(\mathbb W_n\) 与 \(\R^n\) 同构，可以定义一个 \(\R^n\times\R^n\to\R^n\) 的叉积。
那么我们说的 \(\R^2\) 上的叉积是啥呢？实际上是 \(\R^3\) 中的叉积的阉割版。\(\R^3\) 中的叉积定义为 \((a,b,c)\times(x,y,z)=\begin{vmatrix}a&x&\hat i\\b&y&\hat j\\c&z&\hat k\end{vmatrix}\)（是个不严格的写法，因为行列式里怎么能同时塞实数和向量呢……只是方便记忆而已）。经过某些证明（此处略去三维叉积的证明），知道对 \(\pmb u,\pmb v\in\R^3\) 有 \(|\pmb u\times\pmb v|=|\pmb u||\pmb v|\sin\theta\)（其中 \(\theta\) 为 \(\pmb u,\pmb v\) 的夹角，选取 \([0,\pi)\) 内的结果，此时 \(\sin\theta\geq 0\) 保证了模长非负；由于在空间中无法定义顺逆时针，所以尖括号运算符不能用），\(\pmb u\times\pmb v\) 的方向与 \(\pmb u,\pmb v\) 垂直，但仅垂直有可能有两种相反的方向，需要使用右手定则确定方向。由于右手定则有反交换律，所以叉积也有反交换律，即 \(\pmb u\times\pmb v=-\pmb v\times \pmb u\)；叉积显然有线性性。我们容易发现 \(|\pmb u\times\pmb v|\) 的式子其实就是 \(\pmb u,\pmb v\) 围成的平行四边形面积，所以叉积的几何意义是「模长为 \(\pmb u,\pmb v\) 围成平行四边形的面积、方向与 \(\pmb u,\pmb v\) 垂直且满足右手定则的向量」。
现在考虑二维向量 \(\pmb u=(a,b),\pmb v=(x,y)\)。考虑加一维 \(0\) 变成三维向量 \((a,b,0),(x,y,0)\)，相当于加了个 \(Z\) 轴，但这两个向量依然安静地躺在 \(XOY\) 平面上。考虑它们的叉积，带入公式发现 \((a,b,0)\times(x,y,0)=(0,0,ay-bx)\)，前两维都是 \(0\)！这样我们就可以定义 \(\R^2\times\R^2\to\R\) 的二维叉积了：\((a,b)\times(x,y)=z_{(a,b,0)\times(x,y,0)}=ay-bx\)。考察其性质，绝对值（也就是对应三维叉积的模长）为 \(|a||b||\sin\lang\pmb u,\pmb v\rang|\)，符号（方向）的话，跟据右手定则，如果 \(\pmb v\) 在 \(\pmb u\) 的逆时针方向则为正（向上），否则为负（向下）。此时我们发现：把 \(\sin\lang\pmb u,\pmb v\rang\) 的绝对值去掉，恰好能表现出 \(\pmb u,\pmb v\) 的顺逆时针关系（因为 \(\sin\) 是奇函数），直接把符号也确定了！于是我们有 \(\pmb u\times\pmb v=|\pmb u||\pmb v|\sin\lang\pmb u,\pmb v\rang\)（其中 \(\pmb u,\pmb v\in\R^2\)），几何意义是 \(\pmb u,\pmb v\) 围成的平行四边形的有向面积。
三维叉积的几何意义我们没有给出证明，但是阉割成二维后的正确性是容易验证的。我们发现 \(ay-bx\) 其实就是 \(\begin{vmatrix}a&x\\b&y\end{vmatrix}\)，那么正确性便不言而喻、一目了然了。于是我们又有：\(\pmb u\times\pmb v=\begin{vmatrix}\pmb u&\pmb v\end{vmatrix}\)。
有了叉积，我们也可以避开三角函数和分类讨论，使用小常数高精度做一些判定：
- 时针方向：\(\pmb v\) 在 \(\pmb u\) 的逆时针方向、共线、顺时针方向分别对应 \(\pmb u\times\pmb v\) 大于、等于、小于 \(0\)（也对应于 \(\langle\pmb u,\pmb v\rang\) 属于 \((0,\pi)\)、\(\{0,\pi\}\)、\((\pi,2\pi)\)）。
- 三点共线判定：\(\pmb u,\pmb v,\pmb w\) 三点共线显然当且仅当 \(\pmb v-\pmb u,\pmb w-\pmb u\) 共线。
struct point {
db operator%(point v) { return x * v.y - y * v.x; } // 叉积
friend int clock(point u, point v) { return sgn(u % v); } // 时针方向（1 表示逆时针，0 表示共线，-1 表示顺时针）
};
极角排序
我们经常需要将一些非零向量按照极角排序（虽然角度是剩余类，但我们依然可以规定其全序关系：取 \([0,2\pi)\) 内的值时的大小关系）。当然可以直接计算极角来排序，但这样有两个显著的缺陷：atan2 常数和精度都很差。常数方面可以预先计算出所有极角，避免进行 \(\mathrm O(n\log n)\) 次 atan2，但仍然改变不了精度差的事实。
注意到叉积可以判断顺逆时针方向，并且只有乘法和加减法，几乎没有精度误差，常数也很小。但是直接判断 \(\pmb v\) 是否在 \(\pmb u\) 的逆时针方向并不能知道它们的极角相对大小。不过注意到如果 \(\pmb u,\pmb v\) 在横轴同侧，那么顺逆时针就可以判断了，\(\pmb v\) 在 \(\pmb u\) 逆时针方向当且仅当 \(\arg\pmb u<\arg\pmb v\)。那么再写 cmp 的时候就可以先计算出 \(\pmb u,\pmb v\) 分别在横轴的哪侧，如果不同测直接判掉，否则用叉积判。对于 \(y=0\) 的向量，认为 \(x>0\) 就在上侧，否则在下侧。
struct point {
friend bool operator<(point u, point v) {
auto up = [&](point u) { return sgn(u.y) > 0 || sgn(u.y) == 0 && sgn(u.x) > 0; };
bool U = up(u), V = up(v);
return U == V ? clock(u, v) > 0 : U;
} // 极角排序
};
有向直线
一条有向直线是在一条直线的基础上规定两个方向中的一个。在 OI 中大部分时候用到的都是有向直线，如果要研究的东西跟方向无关，那直接忽略方向就好了。
对一条有向直线，任取上面的一点当作新坐标系的原点，直线方向为横轴正方向，在横轴上方的点集规定为该有向直线的左侧，在横轴下方则规定为右侧。
记录
通常记录任意一个直线上的点 \(\pmb p\) 与方向向量 \(\pmb v\)（不一定要是单位向量，只要不为 \(\pmb 0\) 且与直线同向即可），那么该直线上的点集可表示为 \(\pmb p+\lambda\pmb v(\lambda\in\R)\)。
struct line {
point p, v;
line() {}
line(point p, point v) : p(p), v(v) {}
};
判断点在直线的哪边
点与直线的位置共有三种：左侧、在直线上、右侧。任取关键点 \(\pmb p\) 和方向向量 \(\pmb v\)，不难发现 \(\pmb u\) 在直线的三种位置分别对应了 \(\pmb u-\pmb p\) 在 \(\pmb v\) 的逆时针方向、共线、顺时针方向，用叉积判一下就好。
struct line {
int side(point u) { return clock(v, u - p); } // 点 u 在直线的哪边（1 表示左边，0 表示在直线上，-1 表示右边）
};
平行、重合、相等判定
两直线平行当且仅当方向向量共线，用叉积判。
两直线重合当且仅当平行且有交点，后者就看一条直线的关键点是否在另一条直线上。
两直线相等当且仅当重合且同向，用方向向量的单位向量相等代替平行判定即可。
struct line {
friend bool para(line l1, line l2) { return clock(l1.v, l2.v) == 0; } // 平行
friend bool coin(line l1, line l2) { return para(l1, l2) && l1.side(l2.p) == 0; } // 重合
friend bool operator==(line l1, line l2) { return l1.v.unit() == l2.v.unit() && l1.side(l2.p) == 0; } // 相等
bool operator!=(line l2) { return !(*this == l2); }
};
求两直线交点
对两不平行（也不重合）的直线 \(l_1,l_2\)，求它们的交点。设 \(l_i\) 的关键点和方向向量分别为 \(\pmb p_i,\pmb v_i\)，设交点为 \(\pmb p_1+\lambda_1\pmb v_1=\pmb p_2+\lambda_2\pmb v_2\)，试求出 \(\lambda_1\)。移项得 \(\lambda_1\pmb v_1-\lambda_2\pmb v_2=\pmb p_2-\pmb p_1\)，这就得到一个二维的满秩矩阵方程。如果要手动求解满秩方阵方程，最理想的方法是克莱姆法则，那么容易得到 \(\lambda_1=\dfrac{\begin{vmatrix}\pmb p_2-\pmb p_1&\pmb v_2\end{vmatrix}}{\begin{vmatrix}\pmb v_1&\pmb v_2\end{vmatrix}}\)。而二维行列式可以写成叉积，即 \(\lambda_1=\dfrac{(\pmb p_2-\pmb p_1)\times\pmb v_2}{\pmb v_1\times\pmb v_2}\)。
struct line {
friend point inter(line l1, line l2) {
point &p1 = l1.p, &v1 = l1.v, &p2 = l2.p, &v2 = l2.v, b = p2 - p1;
db lambda1 = (b % v2) / (v1 % v2);
return p1 + lambda1 * v1;
} // 直线交点
};
点到直线的距离
也很简单……\(\pmb p,\pmb p+\pmb v\) 显然是直线上的两点，记为 \(\pmb p_{1\sim2}\)。则点 \(\pmb u\) 与 \(\pmb p_1,\pmb p_2\) 组成的三角形显然是以 \(\pmb p_1\pmb p_2\) 为底、\(\pmb u\) 到直线的垂线段为高的三角形，显然可以用叉积（的绝对值）求出三角形面积（的两倍），而底边也可以求出，除一下就是高长了。答案为 \(\dfrac{|\pmb v\times(\pmb u-\pmb p)|}{|\pmb v|}\)。
struct line {
db dis(point u) { return abs(v % (u - p)) / ~v; } // 点到直线的距离
};
线段
是不是越来越水了……毕竟有了向量基操，一切都戎胤了起来。
记录
记录两个端点即可，顺序无所谓，因为线段是没有方向的。
附赠一个判断重合（即端点无序对相等）和判断平行（即所在直线平行）。
struct segment {
point a, b;
segment() {}
segment(point a, point b) : a(a), b(b) {}
bool operator==(segment l2) { return a == l2.a && b == l2.b || a == l2.b && b == l2.a; }
bool operator!=(segment l2) { return !(*this == l2); }
friend bool para(segment l1, segment l2) { return para(line(l1.a, l1.b - l1.a), line(l2.a, l2.b - l2.a)); } // 平行
};
判断点是否在线段上
先判断点是否在线段所在直线上。然后判断点是否在以该线段为对角线的边平行于坐标轴的矩形内。不难发现，这两个条件合起来是充要条件。
struct segment {
bool on(point u) {
return line(a, b - a).side(u) == 0 &&
sgn(min(a.x, b.x) - u.x) <= 0 && sgn(u.x - max(a.x, b.x)) <= 0 &&
sgn(min(a.y, b.y) - u.y) <= 0 && sgn(u.y - max(a.y, b.y)) <= 0;
} // 判断点是否在线段上
};
判断两条线段是否相交
这里指的相交不再排斥所在直线重合的情况，只要有至少一个公共点则称为相交。
分成两条线段是否平行讨论。如果所在直线平行而不重合那就寄了。如果所在直线重合，只要判断分别以这两条线段为对角线的边平行于坐标系的矩形的交是否为空即可。
如果不平行，那么相交当且仅当 \(l_1\) 的两端在 \(l_2\) 所在直线的不同侧（在直线上也算一侧），\(l_2\) 对 \(l_1\) 也如此。求四遍 side 即可。
struct segment {
friend bool inter(segment l1, segment l2) {
line L1(l1.a, l1.b - l1.a), L2(l2.b, l2.b - l2.a);
if(para(L1, L2)) {
if(!coin(L1, L2)) return false;
db l = max(min(l1.a.x, l1.b.x), min(l2.a.x, l2.b.x)), r = min(max(l1.a.x, l1.b.x), max(l2.a.x, l2.b.x));
db d = max(min(l1.a.x, l1.b.x), min(l2.a.x, l2.b.x)), u = min(max(l1.a.x, l1.b.x), max(l2.a.x, l2.b.x));
return sgn(l - r) <= 0 && sgn(d - u) <= 0;
} else {
return L1.side(l2.a) != L1.side(l2.b) && L2.side(l1.a) != L2.side(l1.b);
}
} // 判断两条线段是否有交点
};
多边形
多边形是一条封闭的不自交的折线。不一定是凸多边形哦！
记录
按照顺时针或者逆时针开一个 vector 记录所有顶点。设顶点按顺序为 \(\pmb p_{0\sim n-1}\)。
using vp = vector<point>;
struct polygon : vp {
using vp::vp;
};
周长
《直接计算即可，简洁即美德。》
答案是 \(\sum\limits_{i=0}^{n-1}|\pmb p_i-\pmb p_{(i+1)\bmod n}|\)。
struct polygon : vp {
db peri() {
db ans = 0;
REP(i, 0, size() - 1) ans += ~((*this)[i] - (*this)[(i + 1) % size()]);
return ans;
}
};
面积
不难想到，如果多边形是凸的，那么任取多边形内（或边上 / 顶点处）的一点 \(\pmb O\)，答案显然是所有 \(\mathrm S_{\triangle \pmb O\pmb p_i\pmb p_{(i+1)\bmod n}}\) 的和（因为 \(\pmb O\) 到所有顶点的连边形成一个三角剖分），式子为 \(\dfrac 12\left|\sum\limits_{i=0}^{n-1}(\pmb p_i-\pmb O)\times(\pmb p_{(i+1)\bmod n}-\pmb O)\right|\)。令人惊讶的是，无论多边形是不是凸的，也无论 \(\pmb O\) 在不在多边形内，这个式子都是成立的！几何上直观的理解就是，由于叉积算出来的是有向面积，所以外面被多算的面积又被消掉了；严谨的证明似乎需要归纳，挺复杂度，此处就略去吧。不过当 \(\pmb O\) 在凸多边形内时，绝对值符号是可以移到 \(\sum\) 里面的，在外面就寄了。
一般为了图方便，取 \(\pmb O=\pmb 0\)。
struct polygon : vp {
db area() {
db ans = 0;
REP(i, 0, size() - 1) ans += (*this)[i] % (*this)[(i + 1) % size()];
return abs(ans) / 2;
}
};
判断点是否在多边形内
这是一个经典问题。有很多方法，这里介绍一种精度和常数都不错的算法——射线法。
先特判掉 \(\pmb u\) 在多边形的顶点上或边上的情况。
从 \(\pmb u\) 往任意方向引一条射线，如果任意顶点都不在这条射线上，并且不存在边躺在这条射线上，那么射线发出的过程中，每与多边形的边相交一次，便切换一次内外位置。最终到无穷远时，必定在多边形外，所以我们只要看射线与多边形的总交点数的奇偶性即可。
现在我们需要选一个方向向量，使得这条射线不包含顶点或完整的边。在数学意义上，给定的点基本都是有理数甚至是整数，那么只要取一个无理数斜率就必然满足条件。在程序中实现时，可以随便写一个奇奇怪怪的方向向量，例如 \((\pi,\phi)\) 这种东西。
射线与每条边显然最多有一个交点，那就求出各自所在直线交点然后判是否在射线、边上。判是否在边上就是一个点在线段上的判断。那如何判断是否在射线上呢？一个无精度损失的方法是看 \(\pmb u\) 是否在 \(\pmb u+\pmb v\) 与 \(\pmb i\) 之间的线段上（其中 \(\pmb v\) 是方向向量，\(\pmb i\) 是交点），如果在上面则不在射线上，反之在射线上。
时间复杂度 \(\mathrm O(n)\)。不难发现这种方法几乎无精度误差，因为只涉及到叉积和实数比较。
struct polygon : vp {
bool in(point u) {
int cnt = 0;
static const point v(pi, (sqrt(5) + 1) / 2);
line l1(u, v);
REP(i, 0, size() - 1) {
segment l2((*this)[i], (*this)[(i + 1) % size()]);
if(l2.on(u)) return true;
point itr = inter(l1, line(l2.a, l2.b - l2.a));
if(!segment(itr, u + v).on(u) && l2.on(itr)) ++cnt;
}
return cnt & 1;
}
};
