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

长平狐 发布于 2013/12/25 17:25
阅读 196
收藏 0

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

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

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

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

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

  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]);
	}

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

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


原文链接:http://blog.csdn.net/zhouxuguang236/article/details/12257763
加载中
返回顶部
顶部