Opencv 笔记7 凸包算法-Graham扫描法

一、凸包

凸包（convex hull）是一个计算几何的概念，通俗的来说就是在一个二维的平面中有n 个点，我们可通过在外面的点画出一个多边型使得所以的点都在这个多边形里面，如下所示

 二、Graham扫描法

思路：先找到一个最左边同时最下面的点作为起始点，然后从这个点开始按照逆时针逐个找凸包上的点，实际上就是进行极角排序依次放入栈中，同时对栈进行操作，排除点那些不非顶点的点。

步骤：

1. 把所有点放在二维坐标系中，那么y方向最小的一定是一个顶点，如果这样的点有两个以上那么可以选择x方向最小的那个点，这样的话就可以找到第一个顶点了。如图中的p0。
2.  p0 作为原点，建立一个坐标系；

   

3. 以p0为顶点，计算剩下几个点的角度，如p1的角度为∠Ap0P1=α,p2点的角度为∠Ap0P2=α1,根据图像的这些调度的大小进行排序，像p2和p3这样的点角度是一样的∠Ap0P2=∠Ap0P3，我们将计算x方向的大小，将小的点放到前面，如p2-->p3;
4. 我们先找到了p0 p1,将其按顺序入栈，然后找第上面已经排序的后一个未入栈的点（待入栈点）,这一步主要是判断已经入栈的栈顶元素是不是凸包的顶点,需要的工具是已经排序的下一个入栈的点。此时我们可以 连接原点和栈顶的那个点成一条直线L，判断那个待入栈点的位置，有两种可能性     a: 该点在直线L的上方，如： 

5. 

   

       如上图所示，p2 在直线L的上方，p1是凸包上的顶点，那么我们可以将p2入栈，此时栈里面的原始依次是p0、p1、p2,如上图所示。      b：在直线的直线上或者下方那么这个栈顶的元素就不是凸包上的顶点，例如判断待判断点p3：此时栈顶元素是p1。综合a、b两种情况循环判断凸包顶点。

   

 5、重复步骤4，直到最后一个点（补充一下，步骤四中，点不会出现在那条直线的下方，因为在第三步中已经根据角度排序了，我看了好几个博客第四步的分类是有问题的，这是我个人的理解，如有问题请不吝指出，谢谢）

　最后，栈中的元素就是凸包上的点了。

图像如下：

在实际写代码的时候看到了下面的方法，这个方法比直线的更加方便，让人容易理解，前面的三步都是一样的，就是后面的第四步有所区别

补充：浅谈叉积的意义：

a=AC=(C-A)=(0-0,2-0)=(0,2);

b=AB=(B-A)=(1-0,1-0)=(1,1)

同理：若 a(x1,y1), b(x2,y2),则：

记作result=a×b=AC×AB=x1y2-y1x2；

result>0  ----->AC 在AB 的左边

result=0  ----->AC AB 共线，这里就要查方向了

result<0  ----->AC 在AB 的右边

代码实现

Mat src(800, 800, CV_8UC3);
vector<Point> points;
vector<Point> statk_points;
Point  firstPoint = {0,0}; // 找到的第一个点，也就是最下面，最左边的点

//在图像上随机生成点  找到第一个点
void createPoints() 
{

	RNG rng = theRNG();
	int count = (unsigned)rng % 100 + 21;

	for (int i = 0; i < count; ++i)
	{
		Point p;
		p.x = rng.uniform(100, 600);
		p.y = rng.uniform(100, 600);
		points.push_back(p);
	}


	int n = points.size();
	// 寻找所有点中的最左边和左下面的点
	for (int i = 0; i < n; i++)
	{
		if (points[i].y> firstPoint.y||(points[i].y == firstPoint.y&& points[i].x == firstPoint.x))
		{
			firstPoint = points[i];
		}
	}
	printf("%d    %d\n", firstPoint.x, firstPoint.y);


	for (int  i = 0; i < n; i++)
	{
		if (firstPoint == points[i])
		{
			statk_points.push_back(points[i]);
			circle(src, points[i], 5, Scalar(0, 0, 255), -1);
		}
		else {
			circle(src, points[i], 2, Scalar(0, 0, 0), -1);
		}
		
	}

}

// 计算两个点之间的距离
double distance(Point p1, Point p2)
{
	return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}

//计算叉积，小于0说明p1在p2的逆时针方向(右边)，即p0p1的极角大于p0p2的极角
//a×b，公式：若a = （x1, y1），b = (x2, y2)，则
//AC×AB = x1y2 - x2y1=b
double cross_product(Point p0, Point p1, Point p2)
{
	return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y);
}

// com 
bool  com(Point& p1,  Point& p2)
{
	double temp = cross_product(firstPoint, p1, p2);
	if (fabs(temp) < 1e-6)//极角相等按照距离从小到大排序
	{
		return distance(firstPoint, p1) < distance(firstPoint, p2);
	}
	else
	{
		return temp > 0;
	}
}
vector<Point> Graham_Scan()
{
	// 1、生产点 并找出第一个点
	createPoints();
	int top = 2;
	printf("%d\n", points.size());
	// 2、极坐标排序
	sort(points.begin() + 1, points.end(), com);
	statk_points.push_back(points[1]);
	statk_points.push_back(points[2]);
	imshow("生产点-第一个点", src);

	for (int i = 3; i < points.size(); ++i)
	{
		while (top > 0 && cross_product(statk_points[top - 1], points[i], statk_points[top]) >= 0)
		{
			--top;
			statk_points.pop_back();
		}
		statk_points.push_back(points[i]);
		++top;
	}
	return  statk_points;
}


int main()
{
	
	vector<Point> p=Graham_Scan();
	int count = p.size();
	Mat newimgage = src.clone();
	for (int i = 0; i < count-1; ++i)
	{
		line(newimgage, p[i], p[i+1], Scalar(0, 0, 255),1,8,0);
	}
	line(newimgage, p[count - 1], p[0], Scalar(0, 0, 255), 1, 8, 0);

	imshow("凸包", newimgage);
	waitKey(0);
	return 0;

}