凸包算法详解

一些废话：写这篇博客是因为在听算法网课老师讲分治的时候提到过分治法求凸包的问题，当时没怎么听懂也没太在意这个（因为老师只是提了一嘴，没有细讲），结果好巧不巧，第二天力扣的每日一题就是凸包问题，题目在这：力扣

然后就想着好好学一下凸包算法，但因为我比较懒。。。立了计划后给这篇博客起了个头就没管了，直到几乎20天后，力扣的每日一题又又又有一题可以用凸包做，题目：力扣

所以我特意抽出一个下午好好看了一下各种凸包算法，写下了这篇博客。

注：一下算法均围绕力扣题目 587.安装栅栏 来讨论

前置知识：

设平面上右三个点p1, p2, p3，选其中任意两点作为线段的两个端点，利用向量叉乘，我们可以快速的判断另一点使在线段的左边还是右边，设向量 p1p2 = (x1, y1);向量 p1p3 = (x2, y2);

当 p1p2 × p1p3 = x1y2 - y1x2 ,当值为正时，p3在线段p1p2的左边，即向量p1p2绕p1逆时针旋转一个角度(0 ~ π)可以和向量p1p3共线，当值为负时，结论与上述刚好相反，当值为0时，p1、p2、p3三点共线。

还有一个结论（p1p2 × p1p3）/ 2 的值为p1 p2 p3三点围成的三角形的面积。

首先想到的肯定是暴力法，枚举点集中所有点对，对点对所连成的直线，判断其余点是否在该直线的同一侧，若是，则这两个点是凸包上的点。时间复杂度为O(n^3)，但这个算法懂得都懂，一般算法题或面试官是不会考你这么无脑的问题的。

方法一、Graham扫描法

算法时间复杂度：O(n log n)

基本思想：先找到凸包上的一个点，然后从那个点开始按逆时针方向逐个寻找凸包上的点。

算法步骤：（假设平面中一共有 m + 1 个点）

1. 找到二维平面中y坐标最小的点（若有多个这样的点，随机选取一个即可），给该点编号为p0，由常理可得p0一定为凸包上的一点
2. 求出其余m个点到p0的极角的值，按极角从小到大对该m个点进行排序（由于p0是所有点中位置最低的点，所以其余点到p0的极角的范围为[0, π]），若有多个点到p0的极角大小相等，则按到p0的距离升序排序，对排好序后的点编号为<p1, p2, ......, pm>
3. 可知p1是到p0极角最小的一点，所以p1一定是凸包上的点，若p1不是凸包上的点，则一定存在一点到p0的极角更小，同理pm也一定是凸包上的点
4. 创建一个栈S，栈S的含义：若栈顶元素为pi，则该栈中存储了点集<p0, p1, ......, pi>所形成的凸包
5. 将p0，p1， p2三个点压入栈S中，很显然，此时栈中存储的是点集<p0, p1, p2>所形成的凸包的顶点（由于p1是整个点集形成的凸包上的点，所以在遍历中，p1也是所有子集形成的凸包上的点，所以栈S中的元素在整个过程中不可能小于2）
6. 依次遍历余下的点<p3, p4, ......, pm>，设某次遍历到点pj，我们假设此时栈顶的点为x，顶点下面的点为y，则对于向量yx和xpj，利用向量叉乘判断xpj相当于yx是不是往外拐（即向右拐），若向右拐，则证明栈顶元素不是点集<p0, p1, ......, pj>所形成凸包的顶点，弹出栈顶元素，继续迭代比较，xpj对于yx不是向右拐，此时将pj入栈，则栈S中存储了<p0, p1, ......, pj>所形成凸包的顶点
7. 直到遍历完点pm后，此时栈S中存储了<p0, p1, ......, pm>所形成的凸包的所有顶点，即我们要求的答案

图片实例（图片来自于算法导论）： 

关于该算法的证明可以查看算法导论33.3节，我反正是感觉很玄学，就是感觉自己看懂了，但是转眼就会忘记（我学算法的时候比较钻牛角尖，看过一些算法的证明，其实都很玄学，需要很强的逻辑思维能力和数学功底才能看懂），其实这个算法我不建议大家看证明，因为凭感觉这个算法就是对的，hhh，像我这种菜鸡看完证明后要我证明这个算法的正确性我还是不会。

算法代码：

int cross(const std::vector<int>& p, const std::vector<int>& x, const std::vector<int>& y) {//计算向量px × py， 当返回值大于0时点y在向量px左边
	return (x[0] - p[0]) * (y[1] - p[1]) - (x[1] - p[1]) * (y[0] - p[0]);
}

std::vector<std::vector<int>> Graham(std::vector<std::vector<int>>& points) {
	int n = points.size(), p0 = 0;
	if (n <= 3)return points;//当n小于等于3时直接返回
	for (int i = 1; i < n; i++) {//寻找p0
		if (points[i][1] < points[p0][1]) p0 = i;
		else if (points[i][1] == points[p0][1] && points[i][1] < points[p0][1]) p0 = i;
	}
	std::vector<int> t = points[0];
	points[0] = points[p0];
	points[p0] = t;
	t = points[0];
	sort(points.begin() + 1, points.end(), [t](const std::vector<int>& a, const std::vector<int>& b)->bool {
		int k = cross(t, a, b);
		if (k == 0)return (a[0] - t[0]) * (a[0] - t[0]) + (a[1] - t[1]) * (a[1] - t[1]) < (b[0] - t[0]) * (b[0] - t[0]) + (b[1] - t[1]) * (b[1] - t[1]);//当叉积为0时证明t, a, b在同一条直线上，此时与t距离更小的点排在前面
		return k > 0;//当叉积大于0时a到t的极角要小于b到t的极角
		});
	std::vector<std::vector<int>> S;
	S.emplace_back(points[0]);
	S.emplace_back(points[1]);
	S.emplace_back(points[2]);
	int size = 3;//栈的大小
	for (int i = 3; i < n; i++) {
		//对于求凸包，可能凸包的一条边上包含多个顶点，有的题需要把所有顶点都算入凸包，有的题只需考虑两个端点，若只需考虑两个端点，把下边的<换成<=即可
		while (cross(S[size - 2], S[size - 1], points[i]) < 0) {
			size--;
			S.pop_back();
		}
		S.emplace_back(points[i]);
		size++;
	}
	//由于当有顶点和p0pn-1共线时（即p0和排序后的最后一个点），我们的算法会漏掉这些点，所以需要单独考虑
	for (int i = n - 2; i >= 0; i--) {
		if (cross(points[0], points[n - 1], points[i]) == 0 && S.size() < n)S.emplace_back(points[i]);//加上S.size() < n是为了防止那些恶心的所有点共线的测试案例
		else break;
	}
	return S;
}

方法二、分治法

算法时间复杂度：O(n log n)

算法思路：

先在点集中找出横坐标最大和最小的连个点p1和pn（若有多个，则任取两个），这两个点必定是凸包上的点，将这两个点加入凸包，这两个点所连成的直线将点集分成了两部分，直线上的部分和直线下的部分，对这两部分分别递归来求凸包上的点。我们拿直线上的那部分点来做例子，找出直线上部所有点中距离直线最远的点pmax，该点一定是凸包上的点（距离点集中任意两点所构成的直线最远的点一定在凸包上，可以自己画一下图证明），连接p1,pmax，pnpmax，现在上部的点集分为3个部分了：在p1,pn,pmax所围成的三角形内的点，p1，pmax左部的点，pn，pmax右部的点，其中第一部分不可能是凸包上的点，我们对第二部分和第三部分的点递归的讨论，这样我们可以把原问题不断的分解为较小的子问题来解决。

对于p1,pmax下部的点，同上。当然还有很多需要注意的小细节，我会在代码里写出来。

小技巧：上面讲过向量叉乘的值除2后是向量所围成的三角形的面积，所以叉乘的值越大，则点离直线的距离越远

算法代码：

int cross(const std::vector<int>& p, const std::vector<int>& x, const std::vector<int>& y) {//计算向量px × py， 当返回值大于0时点y在向量px左边
	return (x[0] - p[0]) * (y[1] - p[1]) - (x[1] - p[1]) * (y[0] - p[0]);
}

std::vector<std::vector<int>> fenzhi(std::vector<std::vector<int>>& points) {
	int min = 0, max = 0;
	int n = points.size();
	if (n <= 3)return points;
	for (int i = 1; i < n; i++) {//找到横坐标最小和最大的两个点
		if (points[i][0] < points[min][0])min = i;
		else if (points[i][0] > points[max][0])max = i;
	}
	std::vector<std::vector<int>> ans;//存放凸包上的点
	std::vector<int> up_points;//存放min,max所连成直线上方的点
	std::vector<int> low_points;//存放min,max所连成直线下方的点
	std::vector<int> on_line;//存放在min,max所连直线上的点
	ans.emplace_back(points[min]);
	ans.emplace_back(points[max]);
	for (int i = 0; i < n; i++) {
		int k;
		if (i != max && i != min) {
			k = cross(points[min], points[max], points[i]);
			if (k > 0) {//该点在直线上方
				up_points.push_back(i);
			}
			else if(k < 0) {//该点在直线下方
				low_points.push_back(i);
			}
			else {//点在直线上
				on_line.push_back(i);
			}
		}
	}
	if (up_points.empty() || low_points.empty()) {//若直线上部或下部没有点，此时直线上的点也属于凸包上的点
		for (int x : on_line)ans.emplace_back(points[x]);
	}
	if (up_points.size() <= 1) {//当上部的点小于1个时，上部的点为凸包上的点
		if (up_points.size())ans.emplace_back(points[up_points[0]]);
	}
	else _fenzhi(min, max, points, up_points, ans);
	if (low_points.size() <= 1) {
		if (low_points.size())ans.emplace_back(points[low_points[0]]);
	}
	else _fenzhi(min, max, points, low_points, ans);
	return ans;
}


//s1, s2是划分点集的两个点，points是原点集，handle是待存放的点集，ans存放答案
void _fenzhi(int s1, int s2, std::vector<std::vector<int>>& points, std::vector<int> handle, std::vector<std::vector<int>>& ans) {
	int n = handle.size(), max = handle[0], max_value;
	max_value = cross(points[s1], points[s2], points[max]);
	if (max_value < 0)max_value = -max_value;
	for (int i = 1; i < n; i++) {//找到距离最远的点
		int t = cross(points[s1], points[s2], points[handle[i]]);
		if (t < 0)t = -t;
		if (t > max_value) {
			max_value = t;
			max = handle[i];
		}
	}
	ans.emplace_back(points[max]);
	bool positive = false;
	std::vector<int> s1_max;//直线s1,max所需要处理的点
	std::vector<int> s2_max;//直线s2,max所需要处理的点
	std::vector<int> on_s1;
	std::vector<int> on_s2;
	//判断条件里面函数的返回值小于0说明s2在直线s1,max的右边，此时只需考虑s1,max左边的点即可;同样可得还需考虑s2,max右边的点
	if (cross(points[s1], points[max], points[s2]) < 0)positive = true;
	for (int i = 0; i < n; i++) {
		if (handle[i] == max)continue;
		int k1 = cross(points[s1], points[max], points[handle[i]]);
		int k2 = cross(points[s2], points[max], points[handle[i]]);
		if (positive) {
			if (k1 >= 0) {
				if (k1 > 0)s1_max.push_back(handle[i]);
				else on_s1.push_back(handle[i]);
			}
			else if (k2 <= 0) {
				if (k2 < 0)s2_max.push_back(handle[i]);
				else on_s2.push_back(handle[i]);
			}
		}
		else {
			if (k1 <= 0) {
				if (k1 < 0)s1_max.push_back(handle[i]);
				else on_s1.push_back(handle[i]);
			}
			else if (k2 >= 0) {
				if (k2 > 0)s2_max.push_back(handle[i]);
				else on_s2.push_back(handle[i]);
			}
		}
	}
	if (s1_max.size() <= 1) {
		if (s1_max.empty()) {//当s1,max直线外没有点时，与这两个点共线的点也是凸包上的点
			for (int x : on_s1) {
				ans.emplace_back(points[x]);
			}
		}
		else ans.emplace_back(points[s1_max[0]]);
	}
	else {
		_fenzhi(s1, max, points, s1_max, ans);
	}
	if (s2_max.size() <= 1) {
		if (s2_max.empty()) {
			for (int x : on_s2) {
				ans.emplace_back(points[x]);
			}
		}
		else {
			ans.emplace_back(points[s2_max[0]]);
		}
	}
	else {
		_fenzhi(s2, max, points, s2_max, ans);
	}
}

 方法三、Jarvis步进法

算法时间复杂度：O(nH)，n是集合中点的总数，H是凸包上点的个数

Jarvis步进法和Graham扫描法十分相似，如果方法一看懂了的话，Jarvis步进法是很好懂的。

Jarvis步进法和Graham算法一样，都是先找出纵坐标最小的点p0（当然，如果你愿意的话，也可以先找出横坐标最小的点，总之就是先找出一个凸包上的点，只不过取纵坐标最小的点，算法实现比较方便），它们的不同处在于，Graham算法每次取出还未处理的，到p0极角最小的点来判断栈中的点是否符合要求，若不符合要求，则弹出，即压入栈中的点不一定就是凸包上的点，而Javis步进法每一步即确定一个凸包上的点，所以该算法叫步进法。当凸包上的点数H = log n 时，Jarvis步进法在渐进意义上比Graham扫描法更快。

算法思路：

Jarvis算法的思路非常简单，必须先从凸包上的某一点开始（比如点集中最左边的点或最下面的点），此处我们选择最下面的点p0，然后寻找点p1，使得所有点落在p0p1的左边或右边，此处我们让所有点都落在向量p0p1的左边，则p1是使向量p0p1极角最小的那个点，然后找点p2，使其他所有点落在向量p1p2的左边，则p2是使夹角β最小的点（β为p0p1与p1p2的夹角），依此类推，每一步找到的点都是凸包上的点。当我们再次找到点p0使，递推结束。

 算法代码：

int cross(const std::vector<int>& p, const std::vector<int>& x, const std::vector<int>& y) {//计算向量px × py， 当返回值大于0时点y在向量px左边
	return (x[0] - p[0]) * (y[1] - p[1]) - (x[1] - p[1]) * (y[0] - p[0]);
}

std::vector<std::vector<int>> Jarvis(std::vector<std::vector<int>>& points) {
	int min = 0, n = points.size();
	if (n <= 3)return points;
	for (int i = 1; i < n; i++) {//寻找p0
		if ((points[min][1] > points[i][1]) || (points[min][1] == points[i][1] && points[min][0] > points[i][0]))min = i;
	}
	std::vector<int> store;
	std::vector<std::vector<int>> ans;
	ans.emplace_back(points[min]);
	points[min][0] = -1;//由于题目数据是坐标在0到100之间，所以将横坐标改为-1表示该点已被选取到
	while (1) {
		min = -1;
		for (int k = 0; k < n; k++) {//找到凸包中β夹角最小的点（即最右边的点）
			if (points[k][0] < 0)continue;
			if (min == -1) {
				min = k;
				store.emplace_back(min);
				continue;
			}
			int y = cross(ans.back(), points[min], points[k]);
			if (y < 0) {
				store.clear();
				store.emplace_back(k);
				min = k;
			}
			else if (y == 0) {
				store.emplace_back(k);
			}
		}
		if (ans.size() == n || cross(ans.back(), points[min], ans[0]) < 0) break;//当点集中所有的点都被选入凸包中或当下一个使β角最小的点为p0时结束循环
		if (store.size() == 1) {
			ans.emplace_back(points[store[0]]);
			points[store[0]][0] = -1;
		}
		else {
			int max = store[0], l = store.size();
			int max_value = (points[max][0] - ans.back()[0]) * (points[max][0] - ans.back()[0]) + (points[max][1] - ans.back()[1]) * (points[max][1] - ans.back()[1]);
			for (int i = 1; i < l; i++) {
				if (max_value < (points[store[i]][0] - ans.back()[0]) * (points[store[i]][0] - ans.back()[0]) + (points[store[i]][1] - ans.back()[1]) * (points[store[i]][1] - ans.back()[1])) {
					max = store[i];
					max_value = (points[store[i]][0] - ans.back()[0]) * (points[store[i]][0] - ans.back()[0]) + (points[store[i]][1] - ans.back()[1]) * (points[store[i]][1] - ans.back()[1]);
				}
			}
			for (int i = 0; i < l; i++) {
				if (store[i] != max) {
					ans.emplace_back(points[store[i]]);
					points[store[i]][0] = -1;
				}
			}
			ans.emplace_back(points[max]);//将距离最远的点作为最后访问的点
			points[max][0] = -1;
		}
		store.clear();
	}
	return ans;
}

方法四、Andrew算法

Andrew算法和Graham算法也极其相似，甚至让我觉得它们只是同一种算法的不同实现方式而已。Graham算法是按极角的大小进行排序，而Andrew算法是按横坐标从小到大进行排序（若横坐标大小相同，则按纵坐标从小到大进行排序）。排完序后的第一个元素和最后一定是凸包上的点，且分别位于凸包的最左侧和最右侧，则我们可以将凸包分为两个上凸壳和下凸壳，我们假设从最左边的点开始沿逆时针出发，则剩下的过程就和Graham算法大致相同了，先将排序后数组最前面的两个点p0, p1依次压入栈S中，遍历数组中的点，如果p2使向量p1p2相对于向量p0p1向左拐或与之共线，则将p2压入栈中，否则将栈顶元素弹出（当栈中元素个数小于2时直接将该元素压入栈中），将数组遍历一遍后，我们可以展出下凸壳中的所有点，此时我们再反向遍历一遍数组，即可找出上凸壳上的所有点。这个方法中，我们不需要显式地考虑共线的点，因为这些点已经按照 x 坐标排好了序。所以如果有共线的点，它们已经被隐式地按正确顺序考虑了(但是当所有点都共线时我们还是需要单独考虑)。同时，由于最左边的点和最右边的点同时属于上凸壳和下凸壳，所以对这两个点要进行特殊处理。

算法代码：

std::vector<std::vector<int>> Andrew(std::vector<std::vector<int>>& points) {
	int n = points.size();
	if (n <= 3)return points;
	std::vector<std::vector<int>> S;
	int s_size = 2;
	sort(points.begin(), points.end(), [](const std::vector<int>& a, const std::vector<int>& b)->bool {
		if (a[0] != b[0])return a[0] < b[0];
		return a[1] < b[1];
		});
	S.emplace_back(points[0]);
	S.emplace_back(points[1]);
	for (int i = 2; i < n; i++) {
		while (s_size >= 2 && cross(S[s_size - 2], S[s_size - 1], points[i]) < 0) {
			s_size--;
			S.pop_back();
		}
		S.emplace_back(points[i]);
		s_size++;
	}
	//当所有点在同一直线上时，我们再去求下凸包的话就会把中间的点重复入栈，所以这里对这种情况特殊处理
	if (s_size == n)return S;
	for (int i = n - 2; i >= 0; i--) {
		while (s_size >= 2 && cross(S[s_size - 2], S[s_size - 1], points[i]) < 0) {
			s_size--;
			S.pop_back();
		}
		S.emplace_back(points[i]);
		s_size++;
	}
	S.pop_back();//由于会重复压入points[0]
	return S;
}

终于结束了，我记得这篇博客我是4月22号打算写的，一直拖到6月4号才写完。。。。。