[计算几何] 2 二维凸包/笨蛋(我)也能看懂的二维凸包算法

二维凸包，这篇博客已经说得够好了，介绍了斜率逼近法、Jarvis算法，Graham算法，还有Andrew算法。我这篇博客只会非常详细的介绍Andrew算法。

数论小白都能看懂的平面凸包详解 - ShineEternal的笔记小屋 - 洛谷博客 (luogu.com.cn)

我相信凭借着我6个粉丝其中5个都是老熟人的传播量，应该不会因为乱贴别人链接导致啥问题（）。也许会有朋友要问了，人家都写的这么好了，你博客的创新点在哪里呢？（谁问你了）我主要解决的是这三个问题：

• 关于三点共线等等的特殊情况，网上的有些代码会被hack掉，我在这里给出一个相对比较靠谱的代码。
• 每个人对Andrew算法的理解可能都有点点不一样，也许我的博客会更适合你。
• 我会尽量探讨，在不同的情况中，板子应该怎么做一些小改动。

那我们开始吧。

目录

• 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的左端点和右端点;

好！水完了！有什么表达不清的地方请联系我哦！