使用基于水平序的 Graham 扫描算法扫描凸包

凸包是什么？你可以想象一面墙上有许多钉子（平面上的点），我们用一根绷紧的橡皮绳把这些钉子包围起来，这个绷紧的橡皮绳就是这个平面上一个凸包。凸包在计算几何中有很多用途，在此不再赘述。这里主要介绍由基于水平序的 Graham 扫描算法扫描凸包的过程。

向量的向量积（叉积）

对于向量 $\vec{a}=(x_1,y_1),\vec{b}=(x_2,y_2)$ ，它们的叉积 $\vec{a}\times\vec{b}=\det\begin{bmatrix}x_1&x_2\\y_1&y_2\end{bmatrix}=x_1y_2-x_2y_1$ 叉积不具有交换律，从行列式推出的式子也很容易看出（毕竟是减法嘛）。

我们需要叉积这个东西干什么呢？用来判断两个向量的方向关系。如果 $\vec{a}\times\vec{b}>0$ ，说明 $\vec{a}$ 在 $\vec{b}$ 的顺时针方向，反之则在逆时针方向。如果结果等于 0，说明这两个向量共线（方向相同或相反）。

向量的其他运算法则数学必修4（人教 A 版）已经学过，在此不再赘述。

用代码来表示向量运算：重载运算符大法！

cpp

struct vec {
	double x, y;
	double operator ^(const vec &rhs) const {	// 向量的叉积
		return x * rhs.y - y * rhs.x;
	}
	double operator *(const vec &rhs) const {	// 向量的点积
		return x * rhs.y + y * rhs.x;
	}
	vec operator +(const vec &rhs) const {		// 向量相加
		return (vec){x + rhs.x, y + rhs.y};
	}
	vec operator -(const vec &rhs) const {		// 向量相减
		return (vec){x - rhs.x, y - rhs.y};
	}
	bool operator <(const vec &rhs) const {		// 为了排序而写的重载运算符，与向量运算无关
		return y == rhs.y ? x < rhs.x : y < rhs.y;
	}
};

Graham 扫描算法

使用栈储存点

我们使用一个栈来储存这些在凸包上的点。这个栈并不是传统意义上的栈，因为它除了需要访问栈顶的元素外，还需要访问挨着栈顶的那个元素，即栈顶下面那个元素。（想知道为什么？往下看）所以特别要注意判断栈内要保证至少有两个元素才能进行这样的访问。

手写栈代码

cpp

struct stck {
	int s[maxN], idx;
	void push(int x) {s[++idx] = x;}// 入栈操作
	void pop() {s[idx--] = 0;}		// 弹出操作
	int top() {return s[idx];}		// 访问栈顶元素
	int top2() {return s[idx - 1];}	// 访问栈顶下面那个元素
	bool empty() {return idx <= 1;}	// 判断是否为空（保证里面至少有两个元素）
};

向量叉积计算

我们要对当前枚举的点与栈顶两个点进行向量计算。假设当前枚举点为 $N(x,y)$ ，栈顶两个点为 $Q(x_1,y_1),P(x_2,y_2)$ ，则需要计算 $\vec{NQ}$ 与向量 $\vec{NP}$ 的叉积。

如果要扫描像下图所示的一个凸包，通过肉眼观察发现 $Q$ 点一定是要出栈的， $N$ 点一定是要入栈的。这意味着扫描像这样的凸包，当向量 $\vec{NQ}$ 在向量 $\vec{NP}$ 顺时针方向（即叉积大于 0）时，栈顶的点需要弹出，并把当前点入栈。

当然有时候可能弹出后栈顶的点依然不满足在凸包上的要求，那就继续执行弹出操作直到符合要求再入栈。下图就是这样一种情况。

反过来想，当向量 $\vec{NQ}$ 在向量 $\vec{NP}$ 逆时针方向（即叉积小于 0）时，直接把当前点入栈即可。

算法过程

先将输入的点排个序（向量的坐标表示也可以当点用），按照纵坐标从小到大排序，如果总坐标相同则横坐标小的点排在前面。（正如重载运算符所写）

然后让我们从最左下方的点开始进行扫描。

以这组数据作为例子。 $A(0,5),B(1,8),C(3,4),D(5,0),E(6,2),F(6,6),G(8,3),H(8,7)$

表示在平面上是这样的

我们按照排好顺序的点从纵坐标小到大枚举每一个点（如果栈为空则入栈），执行上述的向量计算过程并进行栈操作。最后栈内的点即是凸包上的点。

对样例这样操作下来只能得到扫描一半的凸包，就像这样

所以还需要从头反过来再扫一遍，原先叉积小于 0 的点是要入栈的，而现在变成了出栈；而原先叉积大于 0 的点现在需要入栈。这样找出的就是一个完整的凸包。

模拟

栈为空，点 $D,E$ 入栈。此时栈内元素： $D,E$

枚举下一个点 $G$ ，进行向量计算，发现 $\vec{GE}\times\vec{GD}=3$ 大于 0，弹出栈顶元素。

再次取出栈顶和栈顶下面的点进行向量计算，符合要求，故点 $G$ 入栈。此时栈内元素： $D,G$

继续枚举下一个点 $C$ ，执行计算，符合要求，点 $C$ 入栈。

对所有的点执行类似操作，直到所有点枚举完。

此时栈内元素： $D,G,H,B$ 清空栈（或直接申请一个新的栈）再次更换计算方式重新枚举。

点 $D,E$ 入栈。

枚举到点 $C$ ，发现 $\vec{CG}\times\vec{CE}=-1$ 小于 0，弹出栈顶元素。

发现 $\vec{CE}\times\vec{CD}=-8$ 小于 0，继续弹出栈顶元素。

将点 $C$ 入栈。此时栈内元素： $D,C$

依然执行类似操作，直到所有点枚举完。

这样一个凸包就扫描出来了。

完整代码

这份代码是用来计算凸包周长的，可以通过洛谷的模板题。

cpp

#include <bits/stdc++.h>
using namespace std;
const int maxN = 2e5 + 10;
struct vec {
	double x, y;
	double operator ^(const vec &rhs) const {
		return x * rhs.y - y * rhs.x;
	}
	double operator *(const vec &rhs) const {
		return x * rhs.y + y * rhs.x;
	}
	vec operator +(const vec &rhs) const {
		return (vec){x + rhs.x, y + rhs.y};
	}
	vec operator -(const vec &rhs) const {
		return (vec){x - rhs.x, y - rhs.y};
	}
	bool operator <(const vec &rhs) const {
		return y == rhs.y ? x < rhs.x : y < rhs.y;
	}
} point[maxN];
double cross(vec a, vec b, vec c) {
	return (b - a) ^ (c - a);
}
double dis(vec a, vec b) {
	vec t = a - b;
	return sqrt(t.x * t.x + t.y * t.y);
}
struct stck {
	int s[maxN], idx;
	void push(int x) {s[++idx] = x;}
	void pop() {s[idx--] = 0;}
	int top() {return s[idx];}
	int top2() {return s[idx - 1];}
	bool empty() {return idx <= 1;}
} sta1, sta2;
double ans;
int n;
int main() {
	scanf("%d", &n);
	for (int i = 1; i <= n; i++)
		scanf("%lf%lf", &point[i].x, &point[i].y);
	sort(point + 1, point + n + 1);
	for (int i = 1; i <= n; i++) {
		while (!sta1.empty() && cross(point[i], point[sta1.top()], point[sta1.top2()]) >= 0)
			sta1.pop();
		sta1.push(i);
	}
	for (int i = 1; i < sta1.idx; i++) ans += dis(point[sta1.s[i]],point[sta1.s[i + 1]]);
	for (int i = 1; i <= n; i++) {
		while (!sta2.empty() && cross(point[i], point[sta2.top()], point[sta2.top2()]) <= 0)
			sta2.pop();
		sta2.push(i);
	}
	for (int i = 1; i < sta2.idx; i++) ans += dis(point[sta2.s[i]],point[sta2.s[i + 1]]);
	printf("%.2lf\n", ans);
	return 0;
}

补充

关于向量的叉积计算有很多种写法。

cpp

double cross(vec a, vec b, vec c) {
	return (b - a) ^ (c - a);
}

这段代码的返回值可以写成 (a - b) ^ (a - c) 或 (a - c) ^ (b - c) 等等。感兴趣的可以试着推导它们的原理，在此不再说明。

本文用到的几何工具为 GeoGebra。