凸包问题的GRAHAM-SCAN解法(附C++代码)

前言

凸包（convex_hull）是一个计算图形学的概念，在二维空间中凸包可以看成一个点集中所有点的最小凸多边形。凸包问题实际上是一个简单的数学问题，只要熟悉原理我们可以很快地编程实现求解，但是网上很多资料都写得比较混乱，推导过程不完整，因此本文将详细介绍凸包问题的推导过程，并附上相应的C++代码。

一、凸包问题

使用GRAHAM-SCAN求解凸包问题的算法流程如下：

上述算法流程出自权威算法类教材，整个流程是绝对正确的，也就是说只要照着上述流程进行代码实现，是可以正确求解凸包问题的。为了让读者更进一步理解SCAN算法，我们将对算法进行更详细的解释。
1、首先找出Q中y坐标最小的点p0，并将它压入S，这是因为p0必然是凸包上的点。
2、按照相对p0的极角从小到大对其他点集中的点进行排序，如果两个点的极角相同，那么我们只保留离p0更远的点，删除另一个点。在代码实现中，极角排序是比较难的部分，推荐读者使用math.h头文件中定义的atan2()函数,可直接对[0,Π]的角度进行排序。在本文代码中，我们使用了更复杂一些的方法，自定义了一种排序方式。
3、我们需要定义一个栈S，用于存储凸包的点。
4、p0进栈
5、p1进栈
6、p2进栈
7、for i=3 to n
8、 前溯过程，当栈顶的前一个元素s1,栈顶元素s2以及当前点pi的叉积小于0时
9、 弹出队顶元素
10、 将pi压入S中
11、返回S
在第8步我们需要计算一个叉积，cross(s1,s2,pi)，若叉积大于0，那么我们不需要回溯，直接将pi压入即可，否则我们需要回溯，不断弹出栈顶元素，并检验新的cross(s1,s2,pi)。假设有点p1(x1,y1)，p2(x2,y2)，p3(x3,y3)，那么我们可以用下图公式计算cross(p1,p2,p3)。

二、凸包实例

上文中我们讲解了凸包问题的SCAN算法，接下来我们通过一个实例展示SCAN算法的求解过程。

如上图所示，我们已经找到了p0点，并按照极角对其进行了排序，顺序为p1、p2、p3,…,p12。
首先我们将p0、p1、p2压入栈中，如图a所示。
在图b中，我们发现p1、p2、p3的叉积小于0，因此弹出p2，而p0、p1、p3的叉积大于0，因此将p3压入栈中。
在图c中，我们发现p1、p3、p4的叉积大于0，因此将p4压入栈中。
在图d中，我们发现p3、p4、p5的叉积小于0，因此弹出p4，而p1、p3、p5的叉积大于0，因此将p5压入栈中。
在图e中，p6进栈。
在图f中，p7进栈。
循环这个过程，直到p12已经被处理，我们就完成了整个求解过程。

三、C++代码

我们使用C++实现了求解凸包问题的SCAN算法，编译环境为VS2019。简单解释一下整个代码的结构，首先我们定义了一个点类c_points，定义了点的x坐标和y坐标，以及它的构造函数，即随机生成点的坐标，作为代码的输入。n是点的数量。最后的输出的凸包所包含的点的坐标。其余的注释代码中基本都有，极角排序函数比较难理解，有问题的读者可私信或留言。

#include <iostream>
#include<algorithm>
#include<stack>
#include<math.h>

using namespace std;

const int n=20;
const double MAXNUM=1e10;


class c_point
{
public:
	double x,y;
	c_point(){x=rand()%100;y=rand()%100;}//构造函数
};

c_point p0;//最低点
c_point points[n];//散点集
stack<c_point>convex_hull;//凸包

//显示栈
void check_stack(stack<c_point>points)
{
	cout<<"==========="<<endl;
	while(!points.empty()){
		cout<<"  ["<<points.top().x<<","<<points.top().y<<"]"<<endl;
		points.pop();
	}
	cout<<"==========="<<endl;
}

//寻找p0
void find_p0()
{
	p0=points[0];
	int ii=0;
	for(int i=1;i<n;i++){
		if(points[i].y<p0.y){
			p0=points[i];
			ii=i;
		}else if(points[i].y==p0.y){
			if(points[i].x<p0.x){
				p0=points[i];
				ii=i;
			}
		}
	}
}

//极角排序
bool cmp(c_point &p1,c_point &p2)
{
	//p0排首位
	if(p1.x==p0.x&&p1.y==p0.y)return true;
	if(p2.x==p0.x&&p2.y==p0.y)return false;

	//计算极角（等于0则赋予一个极大值）
	double angle1=p1.x==p0.x?MAXNUM:(p1.y-p0.y)/(p1.x-p0.x);
	double angle2=p2.x==p0.x?MAXNUM:(p2.y-p0.y)/(p2.x-p0.x);
	//小于0则赋予一个更大的值
	if(angle1<0)angle1+=2*MAXNUM;
	if(angle2<0)angle2+=2*MAXNUM;
	
	//极角排序
	if(angle1<angle2)return true;
	else if(angle1==angle2){
		if(p1.y>p2.y)return true;
		else return false;
	}
	else return false;
}

//叉积
double cross(c_point p1, c_point p2, c_point p3)
{
	return (p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);
}

//搜索凸包
void find_convex_hull()
{
	//p0和p1是凸包中的点
	convex_hull.push(points[0]);
	convex_hull.push(points[1]);

	int i=2;
	//p1,p2为栈顶两个节点
	c_point p1=points[0];
	c_point p2=points[1];
	while(i<n){
		//如果points[i]和points[i-1]在同一个角度，则不再对points[i]进行计算
		if((points[i-1].y-p0.y)*(points[i].x-p0.x)==(points[i-1].x-p0.x)*(points[i].y-p0.y)){
			i++;
			continue;
		}

		//如果叉积大于0，将当前点压入栈
		if (cross(p1, p2, points[i])>=0){
			//假设现在栈中为a,b,c,d,cross(c,d,e)大于等于0
			convex_hull.push(points[i]);//a,b,c,d,e,p1=c,p2=d
			p1=p2;//p1=d
			p2=convex_hull.top();//p2=e
			i++;
		}

		//如果叉积小于0，对栈中节点进行处理
		else{
			while(1){
				//假设现在栈中为a,b,c,d,cross(c,d,e)小于0
				convex_hull.pop();//a,b,c
				convex_hull.pop();//a,b
				p2=p1;//p2=c;
				p1=convex_hull.top();//p1=b
				convex_hull.push(p2);//a,b,c
				//cross(b,c,e)
				if(cross(p1,p2,points[i])>=0){
					convex_hull.push(points[i]);//a,b,c,e
					p1=p2;//p1=c
					p2=convex_hull.top();//p2=e
					i++;
					break;
				}
			}
		}
	}
}

int main()
{
	find_p0();//寻找p0
	sort(points,points+n,cmp);//按极角排序
	find_convex_hull();//搜索凸包
	check_stack(convex_hull);//显示结果

	system("pause");
}

总结

写本文的主要目的是使读者能快速理解并求解凸包问题，网上虽然有各种资料，但是确实都有所缺失，使得读者难以掌握凸包的整个求解过程，本文写得较为详细，并附有完整代码，以便读者学习借鉴。