这篇博文是去年发出来的,由于某种原因删除了,现在重新挂出来和大家分享。代码的格式可能有点乱了。
在地理信息系统(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]);
- }
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]);
}
而这个算法里面的求矢量叉积,点之间距离,点与多边形关系判断的函数的代码就不写上了,也比较简单。其最后的测试效果如下图所示:
大家还有什么更好的方法也可以一起讨论,如果写的不好就当是看了一些废话。