算法学习笔记(26): 计算几何

合集 - 算法学习笔记(52)

1.算法学习笔记(∞)：杂项2023-08-172.算法学习笔记(1): 欧几里得算法及其扩展2023-01-123.算法学习笔记(2): 欧拉定理与逆元2023-01-134.算法学习笔记(3): 倍增2023-01-125.算法学习笔记(3.1): ST算法2023-10-136.算法学习笔记(4): 并查集及其优化2023-01-127.算法学习笔记(5): 最近公共祖先(LCA)2023-01-128.算法学习笔记(6): 树链剖分2023-01-129.算法学习笔记(7): 二分图2023-01-1310.算法学习笔记(8): 网络流2023-01-1311.算法学习笔记(8.0): 网络流前置知识2023-01-1312.算法学习笔记(8.1): 网络最大流算法 EK, Dinic, ISAP2023-01-1413.算法学习笔记(8.2): 上下界网络流2023-01-1414.算法学习笔记(8.3): 网络最大流 - 模型篇2023-12-1915.算法学习笔记(9): 中国剩余定理(CRT)以及其扩展(EXCRT)2023-01-1516.算法学习笔记(10): BSGS算法及其扩展算法2023-01-1617.算法学习笔记(11): 原根2023-01-1618.算法学习笔记(12): 线性基2023-02-0319.算法学习笔记(13): Manacher算法2023-02-0720.算法学习笔记(14): 字符串哈希2023-02-0621.算法学习笔记(15): Trie（字典树）2023-02-0822.算法学习笔记(16): 组合数学基础2023-02-0823.算法学习笔记(17): 快速傅里叶变换（FFT）2023-02-1024.算法学习笔记(18): 平衡树（一）2023-03-1025.算法学习笔记(19): 树上启发式合并（DSU on tree）2023-03-1826.算法学习笔记(20): AC自动机2023-03-2627.算法学习笔记(21): 平衡树（二）2023-05-1328.算法学习笔记(22): 逆序对与原序列2023-05-3129.算法学习笔记(23): 马尔可夫链中的期望问题2023-06-0130.算法学习笔记(24): 狄利克雷卷积和莫比乌斯反演2023-07-2731.算法学习笔记(25): 矩阵树定理2023-06-16

32.算法学习笔记(26): 计算几何2023-07-22

33.算法学习笔记(27): 后缀结构2023-07-2634.算法学习笔记(28): 筛法2023-07-2635.算法学习笔记(29)：分块2023-09-2036.算法学习笔记(30)：Kruskal 重构树2023-10-1337.算法学习笔记(31): 李超线段树2023-10-2038.算法学习笔记(32): 分治2023-10-2739.算法学习笔记(33): 格路径与计数2023-10-2940.算法学习笔记(34): 矩阵乘法与线段树标记2023-11-0341.算法学习笔记(35): CMD Tree2023-11-0342.算法学习笔记(36): 期望中的停时2023-11-0343.算法学习笔记(37): 点分治，边分治小记2023-11-0744.算法学习笔记(38): 矩阵2023-11-1345.算法学习笔记(39): 2-SAT2023-11-1346.算法学习笔记(40): 具体数学2023-11-2247.算法学习笔记(41): 朴素多项式算法2023-11-2248.算法学习笔记(42): 颜色段均摊2023-11-2349.算法学习笔记(43): 可持久化线段树 - 区间加！2023-11-2550.算法学习笔记(44): 二维问题小计2024-01-3151.算法学习笔记(45): 快速沃尔什变换 FWT2024-02-1952.算法学习笔记(46): 离散余弦变换（DCT）2024-03-16

收起

计算几何

向量

高一知识，略讲。

向量外积

若 $\vec x = (x_1, y_1), \vec y = (x_2, y_2)$，则有 $\vec x \times \vec y = x_1 y_2 - y_1 x_2$。

或者表示为 $|\vec x||\vec y| \sin \theta$，其中 $\theta$ 表示向量间的夹角。

几何意义：两个向量构成的平行四边形的面积（可以为负数）。

性质：

• 若 $\vec x \times \vec y \gt 0$，则 $\vec y$ 在 $\vec x$ 的逆时针方向，反之亦然。

向量旋转

将一个向量旋转 $\alpha$，可以利用矩阵的性质。

对于向量 $\vec x = (x, y)^T$，旋转公式为：

$$(x, y) \begin{pmatrix} \cos \alpha & \sin \alpha \\ -\sin \alpha & \cos \alpha \end{pmatrix}$$

矩阵只是为了方便表示，通过欧拉公式 $e^{ix} = \cos x + i \sin x$ 可以简单的推导。

注意到向量可以表示为长度和角度，那么就是 $k e^{i \alpha}$，旋转则可以表示为乘上一个 $e^{\beta}$，那么：

$$\begin{aligned} k e^{i \alpha} &= k \cos \alpha + k i \sin \alpha \\ k e^{i (\alpha + \beta)} &= k e^{i \alpha} e^{i \beta} \\ &= k (\cos \alpha + i \sin \alpha) (\cos \beta + i \sin \beta) \\ \end{aligned}$$

将 $k \cos \alpha$ 看为 $x$，$k \sin \alpha$ 看为 $y$，那么：

$$\begin{aligned} &= x \cos \beta - y \sin \beta + i x \sin \beta + i y \cos \beta \\ &\to (x \cos \beta - y \sin \beta, x \sin \beta + y \cos \beta) \end{aligned}$$

利用矩阵表示就是上文所示了。

向量的极角

定义为向量与 $x$ 轴正半轴的夹角，一般用 $atan2(y, x)$ 实现。

$atan(y, x)$ 的取值范围为 $[-\pi, \pi]$。

点

利用一个向量 $\vec x$ 表示其坐标。

这个向量等价于原点到目标点的向量。

线

可以用多种方法表示，视情况而定：

• 一般式表达：$Ax＋By + C = 0$。方便表示一条直线，或者无向的线段。

• 两个点表达：$(x_1, y_1) \to (x_2, y_2)$，可以方便表示有向的线段或者射线。

• 点与向量：$(x,y) + k \vec d$，可以方便的表示射线。

补充：

两个点求解一般式，有 $A = y_2 - y_1, B = x_1 - x_2, C = x_2y_1 - x_1 y_2$。

判断两个一般式直线平行，用：$A_1B_2 = A_2B_1$。本质是求斜率。

求一般式的交点，将式子变为：

$$A_1 A_2 x + B_1 A_2 y + C_1 A_2 = 0 \\ A_2 A_1 x + B_2 A_1 y + C_2 A_1 = 0 \\ $$

相减可得：

$$y = (C_1 A_2 - C_2 A_1) / (A_1 B_2 - A_2 B_1)$$

同理可得：

$$x = (C_2 B_1 - C_1 B_2) / (A_1 B_2 - A_2 B_1)$$

需要保证两条直线不平行！

线线交点

首先排除平行与重合的情况，判断是否相交，利用向量外积即可。

判断有无交点：若 $\vec {AC} \otimes \vec {AD}$ 和 $\vec{AC} \otimes \vec {AB}$ 异号，以及 $\vec {BD} \otimes \vec {BA}$ 和 $\vec {BD} \otimes {BC}$ 异号，则有交点。

然后可以利用面积法求交点。

用向量外积求出 $S_{ABCD}$，以及 $S_{ABD}$。

那么 $AO : AC = S_{ABD} : S_{ABCD}$，然后就可以找出 $O$ 的坐标了。

---

点在线上

线上去两点 $S, T$，对于点 $P$，若 $P$ 在线段 $ST$ 上，则 $\vec{SP} \otimes \vec {PT} = 0$。

反之不一定成立，需要再判断坐标范围。

点到线的距离

也就是垂线长度。在直线上任取两个点，利用向量外积求出所构成的三角形的面积，除以底边长度，进而求出垂线长度。

double height(const Point &dot, const Point &st, const Point &ed) {
    const Point vec = ed - st;
    return (vec ^ (dot - st)) / vec.abs();
}

代码中 vec.abs() 表示向量的模长，也就是底边的长度。

点到线的垂线段

还是找到两个点。

利用这两个点，找到在直线方向上的单位向量，旋转 $\frac \pi 2$，乘上垂线长，与原坐标相加即可。

Point foot(const Point &dot, const Point &st, const Point &ed) {
    const Point vec = ed - st;
    double h = (vec ^ (dot - st)) / vec.abs();
    return dot + (Point){vec.y, -vec.x} / vec.abs() * h;
}

点关于直线对称

类似于找垂足的方法。只是再加上的向量的基础上 $\times 2$。

Point flip(const Point &dot, const Point &st, const Point &ed) {
    const Point vec = ed - st;
    double h = (vec ^ (dot - st)) / vec.abs();
    return dot + (Point){vec.y, -vec.x} / vec.abs() * h * 2;
}

可能需要注意方向！

---

多边形

一般按照顺时针或者逆时针的方向一一列出所有的顶点。

判断顺时针或者逆时针

钦定一个起点，编号为 $1$。

枚举 $i$ 利用 $\vec{1i}$ 和 $\vec{1(i+1)}$ 的叉乘，可以算出整个多边形的面积(的两倍)。

但是考虑到叉乘的正负性，如果结果为正，则所给的顺序为逆时针（因为 $\vec{1i}$ 在 $\vec{1(i+1)}$ 的顺时针方向）。

判断凸多边形

判断折线段拐向是否相同：

对于折线 $A \to B \to C$，若 $(B - A) \otimes (C - B) \lt 0$，则为顺时针，反之为逆时针。

判断点在多边形内

有很多方法，这里说两种常用的。

• 通用射线法：从该点做一条射线，如果该射线与多边形的交点数为奇数，则在内部，反之在外部。

• 凸多边形二分法：注意的是，基准点需要作为原点！因为凸多边形可以被划分为若干个不相交的三角形。那么我们只需要二分其右边最左边的边界即可。注意的是需要特判一下在下面的情况：

  // cvx 指的是凸包，c 表示凸包的大小，vec 表示所求的点的坐标
  int checkIn(const Point *cvx, int c, const Point &vec) {
      // 特判
      if ((vec ^ cvx[0]) > 0 || (cvx[c - 1] ^ vec) > 0) return 0;
  
      // + 二分
      int lt = 0;
      for (int w = 1 << (int)log2(c); w; w >>= 1) {
          if (lt + w < c && (cvx[lt + w] ^ vec) > 0)
              lt += w;
      }
  
      // 判断是否在三角形
      return ((cvx[lt + 1] - cvx[lt]) ^ (vec - cvx[lt])) >= 0;
  }
  

---

凸包算法

Graham 算法

首先极角排序，然后单调栈做一遍，好写不易错！

不过需要注意的是，在极角相同的情况下，要按照距离 原点 排序。

这是核心代码：

using namespace std;
const int N = 2e5 + 7;
const double eps = 1e-9;

inline int dcmp(double x, double y) {
	if (x > y + eps) return 1;
	if (y > x + eps) return -1;
	return 0;
}

struct Point {
	double x, y;
	bool operator < (const Point &p) const {
		return dcmp(x, p.x) ? dcmp(x, p.x) == -1 : dcmp(y, p.y) == -1;
	}
	
	Point operator - (const Point &p) const {
		return {x - p.x, y - p.y};
	}
	
	double operator ^ (const Point &p) const {
		return x * p.y - y * p.x;
	}
	
	Point rotate(double theta) const {
		double c = cos(theta), s = sin(theta);
		return {x * c - y * s, x * s + y * c};
	}
	
	double atan() const {
		return atan2(y, x);
	}
} P[N], stk[N];

inline bool cmp(const Point &x, const Point &y) {
	double w[2] = {(x - P[0]).atan(), (y - P[0]).atan()};
	return dcmp(w[0], w[1]) ? w[0] < w[1] : (dcmp(x.x, y.x) ? x.x < y.x : x.y < y.y);
}

inline int slopeCmp(const Point &x, const Point &y, const Point &z) {
	return dcmp((x - z) ^ (y - z), 0);
} 

inline double dis(const Point &x, const Point &y) {
	double dx = x.x - y.x, dy = x.y - y.y;
	return sqrt(dx * dx + dy * dy);
}

void ConvexHull(int n) {
	for (int i = 0; i < idx; ++i) {
		if (P[i] < P[0]) swap(P[i], P[0]);
	}
		
	sort(P + 1, P + idx, cmp);
	
	int top = 0;
	stk[++top] = P[0];
	
	for (int i = 1; i < idx; ++i) {
		while (top >= 2 && slopeCmp(P[i], stk[top], stk[top - 1]) >= 0) --top;
		stk[++top] = P[i]; 
	}
}

例题：

• [USACO5.1] 圈奶牛Fencing the Cows /【模板】二维凸包 - 洛谷

• [SHOI2012] 信用卡凸包 - 洛谷

记录

Andrew 算法

没意思，不好用，虽然蛮快的，但是写着很容易错！

---

旋转卡壳

题面：[USACO03FALL] Beauty Contest G /【模板】旋转卡壳 - 洛谷

经典凸包上应用。用于求各种神秘的最远东西。

偷一张图……QwQ

这就是旋转卡壳的流程。

显然是双指针。但是不一样。这里枚举的不是每一个点，而是每一条边。

为什么？考虑对于点来说，其最远的点不一定是单调的。考虑这样一个图形：

成功的卡掉了点点的双指针。

然而，对于一条线来说，凸包上的点到他的距离（垂线长）是单调的，距离也随之改变。

而距离的改变意味着面积的变化，得出可以利用向量外积！

所以这里的双指针是一个线段，一个点。

lint ans = 0;
for (int l = 0, r = 0; l < csiz; ++l) {
    ans = max(ans, dis(convex[l], convex[l + 1]));
    while (cross(convex[l], convex[l + 1], convex[r])
        < cross(convex[l], convex[l + 1], convex[r + 1])) {
        ++r;
    }

    ans = max(dis(convex[r], convex[l]), ans);
    ans = max(dis(convex[r], convex[l + 1]), ans);
}

---

例题：[HNOI2007] 最小矩形覆盖 - 洛谷

旋转卡壳的高级用法！很适合练手。

这里不仅仅需要求一个方向的最远，而是三个方向。

所以是维护三个指针……QwQ

double ans = 1e18;
int l = 0, m = 1, r = 1, mi;
for (int i = 0; i < n; ++i) {
    // 维护最上端点
    while (height(A[m + 1], A[i], A[i + 1]) > height(A[m], A[i], A[i + 1]) + eps)
        ++m;
    double h = height(A[m], A[i], A[i + 1]);

    const Point vec = A[i + 1] - A[i];
    Point ed = A[i] + (Point){-vec.y, vec.x};
    while (height(A[l + 1], A[i], ed) > height(A[l], A[i], ed) + eps)
        ++l;

    // 保证第一次是往左边跑，其他时候这个东西没有影响
    while (height(A[(l - 1 + n) % n], A[i], ed) > height(A[l], A[i], ed) + eps)
        l = (l - 1 + n) % n;

    double w = height(A[l], A[i], ed);

    ed = A[i] + (Point){vec.y, -vec.x};
    while (height(A[r + 1], A[i], ed) > height(A[r], A[i], ed) + eps)
        ++r;

    w += height(A[r], A[i], ed);

    rct[i] = {i, l, m, r};
    if (w * h + eps < ans) {
        ans = w * h;
        mi = i;
    }    
}

---

闵可夫斯基和

对于两个凸包 $A, B$，其闵可夫斯基和定义为 $\{\vec v| \vec v = \vec a + \vec b , \vec a \in A, \vec b \in B\}$。

如果画一下，可以发现，两个凸包的闵可夫斯基和就是凸包的所有边极角排序一遍。

虽然可以 $O(n \log n)$ 完成，但是显然不够优秀。

考虑求出的两个凸包，一定是按照极角已经有序了的。

所以可以借鉴归并排序的思路搞一遍。

int Minkowski(Point *pa, Point *pb, int n, int m) {
    static Point ta[N], tb[N];
    for (int i = 0; i < n; ++i)
        ta[i] = pa[(i + 1) % n] - pa[i];
    for (int i = 0; i < m; ++i)
        tb[i] = pb[(i + 1) % m] - pb[i];

    int idx = 0, i = 0, j = 0;
    minkow[idx++] = pa[0] + pb[0];
    while (i < n || j < m) {
        if (j == m || (i < n && (ta[i] ^ tb[j]) >= 0))
            minkow[idx++] = ta[i++];
        else
            minkow[idx++] = tb[j++];
    }

    for (int i = 1; i < idx; ++i)
        minkow[i] = minkow[i] + minkow[i - 1];
    return idx;
}

注意一下，这里的求法有一些前提：

• 两个凸包起始点的斜率范围相同！

• unkown error

---

例题：[JSOI2018] 战争 - 洛谷

很好的综合题。

需要使用 Minkowski，以及神秘的二分判断点在凸多边形内。

首先转换模型。

对于题目的要求是判断命题：$\exists \vec b \in A, \vec b + \vec w \in A$。或者表示为：

$$\exists \vec a \in A, \vec b \in B, \vec b + \vec w = \vec a$$

转化来说，成为判断：

$$\vec w \in \{\vec v | \vec v = \vec a - \vec b, \vec a \in A, \vec b \in B \}$$

不难发现，后者其实就是一个 Minkowski！

由于闵可夫斯基和还是一个凸包，所以，可以满足二分法判断是否在凸包内。

然后，就搞定了！

---

半平面交

好抽象，感觉很简单又不太会。

核心还是栈和极角排序。根据 @Bill666 说就是先弹队尾，再弹队首，再弹队尾。

不过在存在平行但是平面不交的情况下会大大的错掉成为 nan，需要特判。大佬的特判比较抽象，考虑什么时候会出现 nan，简单思考下来可能就只有这种情况，那么用 isnan 判为 0 即可。

小Z教化学 需要求凸包的交，也就变成了半平面交。但是存在的问题是如何处理交为一条直线的情况？问了一下 UOJ 也没有回答……只能暂时咕咕。

---

$\downarrow 2024.01.07$
我去看了一圈 SCOI 历年的表现，好像都有计算几何的一道题。
但是 AC < 2？所以需要的并不是正解，而是优雅的暴力？
可能之后需要锻炼一下乱搞的能力。