基于边界四边形的凸包生成

这篇博文是去年发出来的，由于某种原因删除了，现在重新挂出来和大家分享。代码的格式可能有点乱了。

在地理信息系统(Geograghic Information System ,GIS) 应用中，原始数据经常是一些离散数据,比如雨量分布数据等,由于数据的采集、传输和录入的顺序不同,一般是一些散乱的数据记录,称其为散乱点集 。
     凸包问题是计算几何中一个重要问题,在GIS中，动态计算面积、裁剪区域、生成TIN 及DTM 都需要用到散乱点集的凸包算法。

    上世纪70 年代以来，不少学者提出了点集凸包的计算方法,较为经典的有增量法、格雷厄姆扫描算法。增量法首先取几个点,形成初始凸包，然后不断寻找当前凸包外的新顶点来更新凸包,直到所有的点都在凸包内。该算法的计算复杂度为O( n2 ) 。格雷厄姆扫描法首先找到最小y 坐标点，接着按照其它点和该极值点的连线与x 轴的夹角的角度值排序，通过判断连续3 个点的空间关系，从而得到逆时针排列的凸包顶点,其计算复杂度为O(nlogn)  。

   作者改进了传统的格雷厄姆扫描算法，首先，将x和y方向上的最小和最大值的点，这四个点必为凸包上的点，然后根据点与多边形的算法，判断散乱点集是否在这四个点组成的边界四边形内，如果在四边形外，则将点加入到一个新开辟的数组。然后用格雷厄姆扫描算法对剩下的这些点计算凸包，这样计算的点就比较少了，从而加快了速度。

本人在VS2008下实现了这个算法，其代码如下：

[cpp] view plain copy print ?

1. void Graham_Scan2(GEOPOINT *pointSet,int n,Stack<GEOPOINT> &s)  
2.     {  
3.         if (n<=3)  
4.         {  
5.             return;  
6.         }  
7.         int k = 0,i = 0;  
8.         int a = 0, b = 0,c = 0, d = 0;   //用于存储四个最值点的索引号   
9.   
10.         //首先寻找x和y的最值   
11.         for (i = 1; i < n; i ++)  
12.         {  
13.             if (pointSet[i].x <= pointSet[a].x)//左   
14.             {  
15.                 a = i;  
16.             }  
17.             if (pointSet[i].y <= pointSet[b].y)//下   
18.             {  
19.                 b = i;  
20.             }  
21.             if (pointSet[i].x >= pointSet[c].x)//右   
22.             {  
23.                 c = i;  
24.             }  
25.             if (pointSet[i].y >= pointSet[d].y)//上   
26.             {  
27.                 d = i;  
28.             }  
29.         }  
30.         GEOPOINT* boundPolygon = new GEOPOINT[5];  
31.         boundPolygon[0] = pointSet[a];  
32.         boundPolygon[1] = pointSet[b];  
33.         boundPolygon[2] = pointSet[c];  
34.         boundPolygon[3] = pointSet[d];  
35.         boundPolygon[4] = pointSet[a];  
36.   
37.         vector<GEOPOINT> ptVec;      //存储排除掉之后的点集   
38.         ptVec.clear();  
39.   
40.         //先排除掉这个四边形内的点   
41.         for (i = 0; i < n; i++)  
42.         {  
43.             int flag = PointInPolygon(pointSet[i].x,pointSet[i].y,boundPolygon,5);  
44.   
45.             if (i == a || i == b || i == c || i == d)//左   
46.             {  
47.                 flag = 2;  
48.             }  
49.   
50.             if (0 == flag || 2 == flag)  
51.             {  
52.                 ptVec.push_back(pointSet[i]);  
53.             }  
54.         }  
55.   
56.         delete [] boundPolygon;  
57.   
58.         k = 0;  
59.         for (i = 1; i < ptVec.size(); i ++)  
60.         {  
61.             if (ptVec[i].y <= ptVec[k].y)//下   
62.             {  
63.                 k = i;  
64.             }  
65.         }  
66.   
67.         GEOPOINT tmp;  
68.         tmp = ptVec[0];  
69.         ptVec[0] = ptVec[k];  
70.         ptVec[k] = tmp;  
71.   
72.         //将剩下的n-1个元素按照与0点的极角排序   
73.         for (i = 1; i < ptVec.size() -1; i++)  
74.         {  
75.             k = i;  
76.             for (int j = i+1; j < ptVec.size(); j ++)  
77.             {  
78.                 if (Miltiply(ptVec[j],ptVec[k],ptVec[0])<0  
79.                     || ((Miltiply(ptVec[j],ptVec[k],ptVec[0])==0)  
80.                     && DistanceToPoint(ptVec[j],ptVec[0]) < DistanceToPoint(ptVec[k],ptVec[0])))  
81.                 {  
82.                     k = j;  
83.                 }  
84.                 tmp = ptVec[i];  
85.                 ptVec[i] = ptVec[k];  
86.                 ptVec[k] = tmp;  
87.             }  
88.         }  
89.   
90.   
91.         //前三个点入栈   
92.         s.push(ptVec[0]);  
93.         s.push(ptVec[1]);  
94.         s.push(ptVec[2]);  
95.   
96.         //从第三点开始去除凹点   
97.         for (i = 3;i < ptVec.size(); i ++)  
98.         {  
99.             //当前点，栈中栈顶点和顶点的下面一个点3个点的转折方向是顺时针方向，就要退栈   
100.             while (Miltiply(ptVec[i],s[s.Size()-1],s[s.Size()-2])<=0)  
101.             {  
102.                 s.pop(tmp);      //出栈   
103.             }     
104.             s.push(ptVec[i]);  //入栈   
105.         }  
106.         s.push(s[0]);  
107.     }  

void Graham_Scan2(GEOPOINT *pointSet,int n,Stack<GEOPOINT> &s)
	{
		if (n<=3)
		{
			return;
		}
		int k = 0,i = 0;
		int a = 0, b = 0,c = 0, d = 0;   //用于存储四个最值点的索引号

		//首先寻找x和y的最值
		for (i = 1; i < n; i ++)
		{
			if (pointSet[i].x <= pointSet[a].x)//左
			{
				a = i;
			}
			if (pointSet[i].y <= pointSet[b].y)//下
			{
				b = i;
			}
			if (pointSet[i].x >= pointSet[c].x)//右
			{
				c = i;
			}
			if (pointSet[i].y >= pointSet[d].y)//上
			{
				d = i;
			}
		}
		GEOPOINT* boundPolygon = new GEOPOINT[5];
		boundPolygon[0] = pointSet[a];
		boundPolygon[1] = pointSet[b];
		boundPolygon[2] = pointSet[c];
		boundPolygon[3] = pointSet[d];
		boundPolygon[4] = pointSet[a];

		vector<GEOPOINT> ptVec;	   //存储排除掉之后的点集
		ptVec.clear();

		//先排除掉这个四边形内的点
		for (i = 0; i < n; i++)
		{
			int flag = PointInPolygon(pointSet[i].x,pointSet[i].y,boundPolygon,5);

			if (i == a || i == b || i == c || i == d)//左
			{
				flag = 2;
			}

			if (0 == flag || 2 == flag)
			{
				ptVec.push_back(pointSet[i]);
			}
		}

		delete [] boundPolygon;

		k = 0;
		for (i = 1; i < ptVec.size(); i ++)
		{
			if (ptVec[i].y <= ptVec[k].y)//下
			{
				k = i;
			}
		}

		GEOPOINT tmp;
		tmp = ptVec[0];
		ptVec[0] = ptVec[k];
		ptVec[k] = tmp;

		//将剩下的n-1个元素按照与0点的极角排序
		for (i = 1; i < ptVec.size() -1; i++)
		{
			k = i;
			for (int j = i+1; j < ptVec.size(); j ++)
			{
				if (Miltiply(ptVec[j],ptVec[k],ptVec[0])<0
					|| ((Miltiply(ptVec[j],ptVec[k],ptVec[0])==0)
					&& DistanceToPoint(ptVec[j],ptVec[0]) < DistanceToPoint(ptVec[k],ptVec[0])))
				{
					k = j;
				}
				tmp = ptVec[i];
				ptVec[i] = ptVec[k];
				ptVec[k] = tmp;
			}
		}


		//前三个点入栈
		s.push(ptVec[0]);
		s.push(ptVec[1]);
		s.push(ptVec[2]);

		//从第三点开始去除凹点
		for (i = 3;i < ptVec.size(); i ++)
		{
			//当前点，栈中栈顶点和顶点的下面一个点3个点的转折方向是顺时针方向，就要退栈
			while (Miltiply(ptVec[i],s[s.Size()-1],s[s.Size()-2])<=0)
			{
				s.pop(tmp);		 //出栈
			}   
			s.push(ptVec[i]);  //入栈
		}
		s.push(s[0]);
	}

而这个算法里面的求矢量叉积，点之间距离，点与多边形关系判断的函数的代码就不写上了，也比较简单。其最后的测试效果如下图所示：

大家还有什么更好的方法也可以一起讨论，如果写的不好就当是看了一些废话。