基于边界四边形的凸包生成
这篇博文是去年发出来的,由于某种原因删除了,现在重新挂出来和大家分享。代码的格式可能有点乱了。
在地理信息系统(Geograghic Information System ,GIS) 应用中,原始数据经常是一些离散数据,比如雨量分布数据等,由于数据的采集、传输和录入的顺序不同,一般是一些散乱的数据记录,称其为散乱点集 。
凸包问题是计算几何中一个重要问题,在GIS中,动态计算面积、裁剪区域、生成TIN 及DTM 都需要用到散乱点集的凸包算法。
上世纪70 年代以来,不少学者提出了点集凸包的计算方法,较为经典的有增量法、格雷厄姆扫描算法。增量法首先取几个点,形成初始凸包,然后不断寻找当前凸包外的新顶点来更新凸包,直到所有的点都在凸包内。该算法的计算复杂度为O( n2 ) 。格雷厄姆扫描法首先找到最小y 坐标点,接着按照其它点和该极值点的连线与x 轴的夹角的角度值排序,通过判断连续3 个点的空间关系,从而得到逆时针排列的凸包顶点,其计算复杂度为O(nlogn) 。
作者改进了传统的格雷厄姆扫描算法,首先,将x和y方向上的最小和最大值的点,这四个点必为凸包上的点,然后根据点与多边形的算法,判断散乱点集是否在这四个点组成的边界四边形内,如果在四边形外,则将点加入到一个新开辟的数组。然后用格雷厄姆扫描算法对剩下的这些点计算凸包,这样计算的点就比较少了,从而加快了速度。
本人在VS2008下实现了这个算法,其代码如下:
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]);}而这个算法里面的求矢量叉积,点之间距离,点与多边形关系判断的函数的代码就不写上了,也比较简单。其最后的测试效果如下图所示:
大家还有什么更好的方法也可以一起讨论,如果写的不好就当是看了一些废话。