制图综合和建筑物数据处理等都涉及到建筑物多边形的化简。制图综合中,由于比例尺的变小,建筑物在小比例尺地图上所占面积变小,这意味着建筑物图形的形状精度也有一定的损失,为了更好地表示原有建筑物的特征(面积、图形和方向),需要对建筑物多边形进行化简。另外,从遥感影像提取的建筑物矢量数据需要经过图形化简等一系列操作才能作为可使用的数据。
首先明确,矢量数据结构中,建筑物多边形是由一些列环(circle)组成,环则由一串首尾相同的点构造。这意味着对建筑物多边形化简其实是对这些环(线状)处理,更具体地说,是对环的各个点进行处理(删除、移位)。
对建筑物多边形的化简到最后得到能较好拟合建筑物的成果多边形,这中间需要一系列操作。首先要进行去除冗余点的操作。
1)冗余点
对于建筑物多边形边界上的邻近三点,若它们所组成两条直线的夹角为平角(或零角),则该三点共线,且位于中间位置那个点的取舍对建筑物多边形的 几 何形状没有影响,因此该三点中位于中间位置那个点是冗余的。
需要注意的是由于数据采集误差、制图不规范等原因,使得地图上建筑物多边形边界上出现冗余点的情况很少,出于方法的合理性考虑,选取5°作为缓冲角度,对冗余点的定义做出以下修正定义:对于建筑物多边形边界上的邻近三点,若它们所组成两条直线的夹角与平角(或零角)相差的角度小于5°时,称该三点中位于中间位置的那个点为冗余点。
2)基本思想
取环上邻近三点作为一个单元,分析夹角:若夹角小于阈值,则删除中间点,加入下一个点构成新的分析单元;若夹角大于阈值,则保留中间点,继续遍历。
3)代码实现
// 去除冗余点,zf,0717void CGeoPolygon::RemoveRedundant(void){ for(int i = 0;i
temp; //用于存放冗余点的下标 bool odinary = true; //设定状态,以确定三点单元的首点 int count = 0; //记录删除的点数,一旦有正常情况恢复为0 int ptsNum = circles[i]->pts.size(); for(int n = 0;n
pts[n]; middle = circles[i]->pts[n+1]; last = circles[i]->pts[n+2]; } else //不正常,存在冗余点的情况 { first = circles[i]->pts[n-count]; //存在冗余点,删除后判断冗余点之后的点是否为冗余点还需要之前三点单元的首点 middle = circles[i]->pts[n+1]; last = circles[i]->pts[n+2]; } double angle = calAngle(first,middle,last); //计算角度 if(angle >= 10&&angle<=170) //角度阈值设为10°(可根据需要调整),当角度大于阈值,情况正常,不需要去除点 { odinary = true; count = 0; //记录删除的点数,一旦有正常情况恢复为0 } else // 存在冗余点 { temp.push_back(n+1); //记录冗余点在circles[i]->pts中下标 odinary = false; //情况不正常 count++; //连续不正常时都要记录 } } vector tempPoints; //定义临时点集 if(tempPoints.size()!=0) vector ().swap(tempPoints); //防止点集不为空 for(int j = 0;j pts[j]); } tempPoints.swap(circles[i]->pts); //交换,得到去除冗余点之后的环circles[i]->pts }} // 计算三点之间的角度,zf,0717,要#include double CGeoPolygon::calAngle(CMyPoint* p1, CMyPoint* p2, CMyPoint* p3){ double x1 = p1->Getx(); double y1 = p1->Gety(); double x2 = p2->Getx(); double y2 = p2->Gety(); double x3 = p3->Getx(); double y3 = p3->Gety(); double angle,p1p2,p2p3,p1p3; p1p2 = sqrt(((x1 - x2) * (x1 - x2) + (y1 -y2) * (y1 -y2))); p2p3 = sqrt(((x3 - x2) * (x3 - x2) + (y3 -y2) * (y3 -y2))); p1p3 = sqrt(((x1 - x3) * (x1 - x3) + (y1 -y3) * (y1 -y3))); angle = acos((p1p2 * p1p2 + p2p3 * p2p3 - p1p3 * p1p3) / (2 * p1p2 * p2p3)) * 180.0 / 3.1415926; if(angle<0) angle = angle*(-1); return angle;}
// 判断某个数是不是在数组中,zf,0717bool CGeoPolygon::isContained(vector temp, int i){ bool answer = false; for(int j = 0;j