[计算几何] 2 二维凸包/笨蛋(我)也能看懂的二维凸包算法
二维凸包，这篇博客已经说得够好了，介绍了斜率逼近法、Jarvis算法，Graham算法，还有Andrew算法。我这篇博客只会非常详细的介绍Andrew算法。
数论小白都能看懂的平面凸包详解 - ShineEternal的笔记小屋 - 洛谷博客 (luogu.com.cn)
我相信凭借着我6个粉丝其中5个都是老熟人的传播量，应该不会因为乱贴别人链接导致啥问题（）。也许会有朋友要问了，人家都写的这么好了，你博客的创新点在哪里呢？（谁问你了）我主要解决的是这三个问题：
- 关于三点共线等等的特殊情况，网上的有些代码会被hack掉，我在这里给出一个相对比较靠谱的代码。
- 每个人对Andrew算法的理解可能都有点点不一样，也许我的博客会更适合你。
- 我会尽量探讨，在不同的情况中，板子应该怎么做一些小改动。
那我们开始吧。
Andrew算法
算法思想
要对一个点集求凸包，就是用一个凸多边形把这些点集全部围起来，这些点最多只能在凸多边形上，不能超出这个凸多边形，并且这个凸多边形是所有可以满足上述条件的凸多边形中，面积、周长均是最小的那一个。并且，这样的一个凸多边形是唯一的。我们从感性上就可以感觉到，这个凸多边形的顶点一定是从点集中选取的。
假如黑色的点是原本的点集，点P是形成的凸包中的某一个顶点。P不属于初始，并且凸包把所有的点围了起来。那么，一定存在初始点集中两个点，把它们相连之后，围成的凸多边形仍然是凸包，并且面积比原来的凸包更小。有了这样一个认知，我们要做的事情就变成了：从点集中选出若干点构成凸包。也可以说，这些点肯定是最“外面”的凸包框住的，对吧？
首先，我们把所有点从左到右的排个序，找到点集的最左边的点A和最右边的点B，点A和点B肯定是最“外面”的点，也就是说，肯定是构成凸包的点。我们想象现在有一块柔软的布从上到下垂下来，盖住了这些点，形成了凸包的上边界：
这个上边界很明显，起始于A，终于B。我们先不想一次找完整个凸包，我们先找上边界。
把点按照这个方法排序，可以得到\(P=\{p_1,...,p_7\}\)
bool operator < (Point& B) {
if (!sgn(x - B.x)) return y < B.y;
return x < B.x;
}
初始点集\(U=\{A\}\)，然后，A先试探性的和\(p_1\)相连，然后把\(p_1\)插入点集U。嗯，很显然\(U=\{A,p_1\}\)就是\(A,p_1\)的上边界。再看\(p_2\)，如果把\(p_1p_2\)相连，那\(U=\{A,p_1,p_2\}\)就是这三个点的上边界。
现在压力来到\(p_3\)，把\(p_2p_3\)相连后，情况开始不对了，为了维护美丽的上边界形状，我们删除\(\overrightarrow{p_1p_2},\overrightarrow{Ap_1}\)，然后把\(p_3\)加进来，现在\(U=\{A, p_3\}\)。
重复上述过程，最终我们得到了\(U=\{A, p_3, p_7, B\}\)
纵观整个寻找上边界的过程，其实可以描述为：按照排序从小到大依次遍历节点，当遍历到第\(i\)个节点时，看节点\(i\)是否当前的上边界\(U\)的末尾两点构成的向量\(\overrightarrow{v}\)的右侧，如果是，则将点\(i\)加入\(U\)，否则，从\(U\)中不断弹出最后一个节点，直到满足条件或者弹到没点可弹为止。
维护下边界，就是反着再来一遍：按照排序从大到小依次遍历节点，当遍历到第\(i\)个节点时，……（不要指望着我再告诉你怎么做了！自己推一下）就像这样：
嗯嗯，最后我们可以得到，\(U=\{A, p_3, p_7, B\}; D=\{ B, p_5, p_2, A\}\)
合起来就是凸包：
代码思路
问题已经解决了，现在我们只要把我们的思想用代码表示出来即可。
回过头来看凸包的形成过程，我们为什么能判断\(p_7\)是新的上边界，而把\(p_6\)弹出了呢？通过仔细的观察，我们发现是因为\(p_7\)这个点在向量\(\overrightarrow{p_3p_6}\)的左侧。如果在维护上边界的时候，点在向量左侧，那就说明\(p_7\)在“更上面”，不是吗？维护下边界也是一样：如果当前的点\(p_1\)在下边界末尾两点构成向量\(\overrightarrow{p_5p_2}\)的右侧，就把当前点加入下边界；如果当前的点在向量左侧，就把下边界末尾的点弹出，直到没点可弹或满足点在向量右侧为止。
判断点和向量（其实这里是有方向的直线）的位置关系，如果你已经看过我的[计算几何] 1 二维基础运算与概念 - 跳岩 - 博客园 (cnblogs.com)的话，应该可以已经非常轻松的写出来代码了罢（）
// 判断直线a和点b的关系
// 1: a在直线左边; 0: a在直线上; -1: a在直线右边
int relation(Line a, Point b){
return sgn(cross(a.v, b - a.p));
}
代码实现
[P2742 USACO5.1] 圈奶牛Fencing the Cows /【模板】二维凸包 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)
现在我们就以这道题为例，看看怎么做一个二维凸包。
andrewScan的伪代码如下，真的代码其实就多了一点点边界的判断。注意看，我们给点排序的时间是\(O(n\log n)\)，找上下边界的时间是\(O(n)\)，故总的时间复杂度为\(O(n\log n)\)。
Polygon andrewScan(Polygon p){
sort(p.begin(), p.end()); // 给点排序
Polygon u, d; // u: 上边界; d: 下边界
u.push_back(p[0], p[1]); // 把最左边的两个点push到上边界
d.push_back(p[end], p[end - 1]); // 把最右边的两个点push到下边界
// 找上边界
for (int i = 2; i < p.size(); ++i){
while (u.size() >= 2 && p[i]不在u末尾向量的右边) 弹出u末尾元素;
u.push_back(p[i]);
}
// 找下边界
for (int i = end - 2; i >= 0; --i){
while (d.size() >= 2 && p[i]不在u末尾向量的右边) 弹出d末尾元素;
d.push_back(p[i]);
}
// 现在找到了上下边界, 我们把结果合并起来
// 本人喜欢逆时针存储
result = merge(d, u, anticlockwise);
return result;
}
AC的代码如下：
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ZERO 1e-8
#define xx first
#define yy second
#define RIGHT -1
#define ON 0
#define LEFT 1
int sgn(double x){
if (fabs(x) < ZERO) return 0;
return x > 0 ? 1 : -1;
}
struct Point{
double x, y;
Point (double _x = 0, double _y = 0) : x(_x), y(_y) {}
Point operator + (Point& B) { return Point(x + B.x, y + B.y); }
Point operator - (Point& B) { return Point(x - B.x, y - B.y); }
bool operator < (Point& B) {
if (!sgn(x - B.x)) return y < B.y;
return x < B.x;
}
friend ostream& operator <<(ostream& out, Point& p) { out << "(" << p.x << ", " << p.y << ")"; return out; }
friend void operator >>(istream& in, Point& p) { in >> p.x >> p.y; }
};
typedef Point Vector;
typedef pair<Point, Point> Line;
typedef vector<Point> Polygon;
double cross(Vector v1, Vector v2){
return v1.x * v2.y - v1.y * v2.x;
}
int relation(Line l, Point p){
return sgn(cross(l.yy - l.xx, p - l.xx));
}
double get_dist(Point p1, Point p2){
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
Polygon andrewScan(Polygon p){
Polygon u, d;
// 如果p的点少于3, 可以直接返回了
if (p.size() < 3) return p;
sort(p.begin(), p.end());
u.push_back(p[0]), u.push_back(p[1]);
d.push_back(p.back()), d.push_back(p[p.size() - 2]);
for (int i = 2; i < p.size(); ++i){
int k = u.size();
while (k >= 2 && relation(Line(u[k - 2], u[k - 1]), p[i]) != RIGHT){
k--; u.pop_back();
}
u.push_back(p[i]);
}
for (int i = p.size() - 3; i >= 0; --i){
int k = d.size();
while (k >= 2 && relation(Line(d[k - 2], d[k - 1]), p[i]) != RIGHT){
k--; d.pop_back();
}
d.push_back(p[i]);
}
// 这里就是做了逆时针合并, 记得u.st=A, u.ed=B, d.st=B, d.ed=A
// 所以合并的时候还要注意不要重复存储两遍A和B
reverse(d.begin(), d.end());
for (int i = u.size() - 2; i > 0; --i){
d.push_back(u[i]);
}
return d;
}
int n;
const int N = 1e5 + 5;
int main(void){
cin >> n;
Polygon pg;
Point p;
for (int i = 0; i < n; ++i){
cin >> p.x >> p.y;
pg.push_back(p);
}
Polygon res = andrewScan(pg);
double dist = 0;
for (int i = 0; i < res.size(); ++i){
dist += get_dist(res[i], res[(i + 1) % res.size()]);
}
cout << fixed << setprecision(2) << dist << endl;
return 0;
}
小小的改动
- 注意这段代码：
for (int i = 2; i < p.size(); ++i) {
int k = u.size();
while (k >= 2 && relation(Line(u[k - 2], u[k - 1]), p[i]) != RIGHT) {
k--; u.pop_back();
}
u.push_back(p[i]);
}
我们的代码中，认为如果当前点不在向量的右边，就要把上边界末尾的点弹出。也就是说，如果\(p_1,p_2,p_3\)三点共线，我们最终会把\(p_2\)弹出，就像下面这样：
如果我们遇到需要保留\(p_2\)的情况，那只要把!=RIGHT改为== LEFT就好啦。
for (int i = 2; i < p.size(); ++i) {
int k = u.size();
while (k >= 2 && relation(Line(u[k - 2], u[k - 1]), p[i]) == LEFT) {
k--; u.pop_back();
}
u.push_back(p[i]);
}
- 而且，我们的代码也是经得起一些hack的：
if (p.size() < 3) : 返回的是p;
if (p中所有点都是三点共线): 返回p的左端点和右端点;
好！水完了！有什么表达不清的地方请联系我哦！