GIS算法基础(七)矢量数据向栅格数据的转换(面转换的边界代数算法实现)
代码已经po上远程仓库,需要的同学可以自取:
https://github.com/XiaoZhong233/GIS_ALG/blob/master/src/scau/gz/zhw/Raster.java
目录
一、边界代数法算法思想
边界代数多边形填充算法是一种基于积分思想的矢量转栅格的转换算法。
这个和格林公式很像,就是根据边界关系求积分,边界代数的思想应该是从这里来的
“设闭区域 由分段光滑的曲线 围成,函数 及 在 上具有一阶连续偏导数,则有
其中 是 的取正向的边界曲线。”
不过我的算法实现过程和高数没什么关系。先来一波高数劝退哈哈。
沿着多边形的边界进行顺时针环绕。当当前环绕方向向上时,则位于该边界左侧(前进方向看为左侧)的具有相同行坐标的所有栅格的值减去a;当当前环绕方向向下时,则位于该边界左侧(前进方向看为右侧)的具有相同行坐标的所有栅格个的值加上a。
当环绕方向为水平方向则不用管。
但在实际操作过程中,如果多边形内部有孔洞,则有可能造成填充值不唯一的情况,不过这没关系,一样可以把多边形的内部和外部区分开,因为外部的值怎么样都为0。内部的值虽然不一定都为a,但是一定不为0。
二、算法步骤
①遍历多边形的每一条边
②判断该边是否为上行方向,若为上行则把前进方向左侧全部-1,若不是则进入步骤③
③判断该边是否为下行方向,若为下行则把前进方向的右侧全部+1。
④重复步骤②,③直至多边形的边遍历完成。
三、算法实现
//矢量化多边形(边际代数法)
public void RasterPolygon2(Polygon polygon) {
List<Line> lines = polygon.getLines();
//RasterLine(polygon, 1);
for(Line line : lines) {
Point start = line.getStart();
Point end = line.getEnd();
int[] startPoint = transformRasterPoint(start, 0, 0);
int[] endPoint = transformRasterPoint(end, 0, 0);
int startRow = Math.max(startPoint[0], endPoint[0]);
int endRow = Math.min(startPoint[0], endPoint[0]);
if(line.isUp()) {
while(startRow>=endRow) {
int column = 0;
while(column<=Math.max(startPoint[1], endPoint[1])) {
Point p = transformVetorPoint(endRow, column, 0, 0);
//保证p在线段的左边
if(Double.doubleToLongBits(p.getX()) < Double.doubleToLongBits(line.getXByY(p.getY()))
&& Double.doubleToLongBits(line.getXByY(p.getY()))>=Double.doubleToLongBits(Math.min(line.getStart().getX(),line.getEnd().getX()))
&& Double.doubleToLongBits(line.getXByY(p.getY()))<=Double.doubleToLongBits(Math.max(line.getStart().getX(), line.getEnd().getX()))) {
addValue(endRow, column, -1);
}
column++;
}
endRow++;
}
}else if(line.isDown()) {
while(startRow>=endRow) {
int column = 0;
while(column<=Math.max(startPoint[1], endPoint[1])) {
Point p = transformVetorPoint(endRow, column, 0, 0);
//保证p在线段的左边
if(Double.doubleToLongBits(p.getX()) < Double.doubleToLongBits(line.getXByY(p.getY()))
&& Double.doubleToLongBits(line.getXByY(p.getY()))>=Double.doubleToLongBits(Math.min(line.getStart().getX(),line.getEnd().getX()))
&& Double.doubleToLongBits(line.getXByY(p.getY()))<=Double.doubleToLongBits(Math.max(line.getStart().getX(), line.getEnd().getX()))) {
addValue(endRow, column, 1);
}
column++;
}
endRow++;
}
}
}
}
四、测试结果
输入数据:
Point d,e,r,t,y,o,q;
o=new Point(10, 22);
d=new Point(12, 25);
e=new Point(20, 25);
r=new Point(22, 18);
t=new Point(13, 18);
y=new Point(16, 10);
Point[] points = new Point[] {o,d,e,r,t,y};
Polygon polygon = new Polygon(points,true);
绘制结果(矢量绘制):
输出结果(控制台输出绘制):
当多边形内部存在孔洞的情况下:
可以发现孔洞处(GIS上的概念对应的应该是岛与洞)的值比外部更大(或者更小,取决于你是顺时针环绕还是逆时针环绕,是上行方向+1还是下行方向+1抑或是别的条件的差异,不过都没关系,不会影响是否在多边形内部关系的判断),这是因为在孔洞的边界多环绕了一次。由这一点我们也可以发现,边界代数法可以很好的保存多边形的拓扑关系,无论这个孔洞是在内部,亦或者是多个多边形邻接,我们可以通过其属性值判断多边形之间的关系。
五、总结
边界代数法与前述面状转换算法(《GIS算法基础(五)》)的不同之处,在于它不是逐个逐个的点判断与边界的关系,而是根据边界的拓扑信息,通过简单的加减代数运算将边界位置信息动态地赋予各栅格点。这就是边界代数法巧妙而优美的地方。实现了矢量格式到栅格格式的高速转换(因为计算简单,量也不大),并且不需要考虑搜索轨迹之间的关系,因此算法简单、可靠性好、各边界弧段只被遍历一次,不需重复计算。