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

计算几何

向量

高一知识，略讲。

向量外积

若 \(\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]; 
	}
}

例题：

记录

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？所以需要的并不是正解，而是优雅的暴力？
可能之后需要锻炼一下乱搞的能力。

posted @ 2023-07-22 09:06  jeefy  阅读(64)  评论(0)  编辑  收藏  举报