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

来源：互联网发布：指纹机可以恢复数据么编辑：程序博客网时间：2024/05/21 06:30

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

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

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

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

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

[cpp] view plaincopyprint?

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

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

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

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