基于GIS空间分析的任意多边形骨架线自动提取方法
针对现有任意多边形骨架线提取方法存在的算法设计繁琐、代码实现复杂、实现成本高和周期长的缺点,提出一种基于地理信息系统(geographic information system,GIS)空间分析的任意多边形骨架线自动提取方法。首先,以电子地图为数据源,在ENVI和ArcGIS平台的支持下,通过调用ArcToolbox工具箱中的ENVI分类模型工具提取出空间对象的多边形矢量边界,并对其进行预处理; 其次,综合运用数据处理、空间分析和文件转换工具提取其骨架结点,并对其进行后处理; 再次,采用Python面向对象编程语言结合ArcPy站点包,通过编写脚本程序自动提取其骨架线; 然后,进一步利用ModelBuilder工具,通过构建模型实现骨架线的自动提取; 最后,将该方法分别应用于电子地图道路和建筑物多边形骨架线的提取。实验过程及其结果表明该方法具有一定有效性、实用性和可操作性。
关键词： 任意多边形 ; 骨架线 ; GIS ; 空间分析
引用格式
宋仁波, 朱瑜馨, 丁上珊, 贺巧宁, 王细元, 王月香. 基于GIS空间分析的任意多边形骨架线自动提取方法. 国土资源遥感杂志[J], 2020, 32(1): 51-59 doi:10.6046/gtzyyg.2020.01.08
SONG Renbo, ZHU Yuxin, DING Shangshan, HE Qiaoning, WANG Xiyuan, WANG Yuexiang. An automatic method for extracting skeleton lines from arbitrary polygons based on GIS spatial analysis. REMOTE SENSING FOR LAND & RESOURCES[J], 2020, 32(1): 51-59 doi:10.6046/gtzyyg.2020.01.08
0 引言
骨架线(中轴)是面状要素特征描述的重要指标,它不仅在模式识别、机助制图和计算机视觉的研究中具有重要的意义[1],而且在地理信息系统(geographic information system,GIS)制图中也得到了广泛应用[2-5]。近些年来,多边形骨架的提取一直是计算机辅助设计、计算机软件、测绘和GIS领域研究的热点[6-9],围绕多边形骨架提取算法的设计、程序实现及应用方面的成果相对丰富[10-17]。现有的骨线提取方法可归纳为基于栅格图像的提取方法[18-21]和基于矢量图形的提取方法[22,23]2大类,但这2类方法都需要繁琐的算法设计和复杂的代码实现,且实现成本高、周期长,在推广和应用中都受到成本和技术条件的制约,所以非常有必要研究和开发出一种实用、有效和可靠性强的多边形骨架线提取方法。
现有基于矢量图形的骨架线提取方法主要包括数据预处理、基于约束Delauny三角剖分的骨架线结点生成和骨架线的连接3个过程,上述过程都可利用现有GIS系统的数据处理、空间分析和建模功能实现。ArcGIS系统不仅提供基础的数据处理、分析和制图功能,而且提供脚本编程批处理、空间分析和建模功能,这些优势都为降低开发成本和实现难度提供了充分的技术条件。本文从GIS空间分析的视角,提出一种基于GIS空间分析的复杂多边形骨架线提取方法。采用Python面向对象编程语言结合ArcGIS系统提供的ArcPy站点包,通过调用ArcToolbox工具箱中的多种空间分析和文件转换工具; 同时,利用ModelBuilder可视化建模工具,通过构建模型实现任意多边形骨架线的自动提取。
1 基本原理与实现方法
基于计算几何学的原理,借助GIS的空间分析和建模技术,同时结合计算机编程技术,探讨任意多边形骨架线自动提取方法。
1.1 基本原理
多边形是平面中由有限线段首尾连接起来划出的封闭几何形状。其中,复杂多边形指的是顶点有凹凸性、且含有岛屿的多边形; 与之相对的为简单多边形。
任意多边形骨架线的结点主要处于相邻三角形公共边的中点上,它恰好是任意多边形的约束Delaunay三角剖分获得的内部邻接三角形公共边的中点,其原理如图1所示。因此,只要在任意多边形的约束Delaunay三角剖分的基础上,采用一定的连接策略连接相邻三角形公共边的中点、重心或边界顶点就可提取出其骨架线。
图1 任意多边形骨架线提取原理示意图
Fig.1 Schematic diagram of the extraction of skeleton lines from arbitrary polygons
1.2 实现方法
从GIS空间分析的角度,借助GIS的数据处理、空间分析、建模和制图功能实现任意多边形骨架线的提取。其中,GIS空间分析指的是借助GIS软件,从空间数据中获取有关地理对象的空间位置、分布、形态、形成和演变等信息并进行分析,其基本功能主要包括空间查询与量算、缓冲区分析、叠加分析、路径分析、空间插值和统计分类分析等。ArcGIS系统的ArcToolbox工具箱除可以提供上述各种基本空间分析工具之外,还能够通过3D分析工具箱的不规则三角网(triangulated irregular network,TIN)工具支持任意多边形区域的约束Delaunay三角剖分,都为任意多边形骨架线的提取提供空间分析基础。
同时,ArcGIS系统内嵌的ModelBuilder是一种可视化编程语言,为创建满足用户要求的特定工作流提供有效的途径。其主要特点是将一系列地理处理工具串联在一起,将其中一个工具的输出作为另一个工具的输入,通过构造和执行工作流,使得任意多边形骨架线的自动提取成为现实。此外,通过创建模型并将其共享为工具来提供扩展ArcGIS系统功能的高级方法,甚至还可用于将ArcGIS系统与其他应用程序进行集成,能够增强构建模型的可复用性和扩展性,为相似软件系统的研发提供借鉴。
此外,Python是一种面向对象、解释型、交互式和面向初学者的计算机程序设计语言,当前,它已成为计算机领域的主流开发语言,其具有语法简捷而清晰,易于编写和维护的优点,并具有丰富和强大的类库[24],能够显著提高程序的开发效率。ArcGIS系统提供了Python语言集成开发环境(IDLE)和ArcPy站点包,通过编程能够创建满足用户特定需求的脚本工具,不仅可以有效弥补ArcGIS系统ArcToolbox工具箱功能的不足,而且能够增强ModeBuilder的建模能力。
2 任意多边形骨架线提取
将任意多边形骨架线提取工作分解为技术流程设计、数据获取和预处理、骨架线结点的生成、骨架线的连接和骨架线提取模型的构建等关键环节。
2.1 技术流程设计
多边形骨架线提取技术流程的设计需要综合考虑实现成本和开发周期,在兼顾可操作性、可靠性、稳定性和效率等方面因素之外,还需要考虑用户操作的便捷性和易用性,其操作步骤主要包括: ① 以电子地图为数据源,经过一定预处理提取出空间对象的矢量多边形,并将其保存至GIS数据库; ②利用3D分析工具箱的创建TIN和TIN转三角形工具,同时结合要素折点转点工具生成任意多边形的约束Delaunay三角网格; ③利用数据管理工具箱的折点处分割线、要素转点和空间连接等工具生成骨架线的结点; ④利用数据管理工具箱的删除重复和空间连接工具对生成的骨架线结点进行后处理; ⑤采用Python面向对象编程语言结合ArcPy站点包,通过编程实现骨架线的连接程序,并结合ArcCatalog在ArcToolbox工具箱构建自定义工具自动提取骨架线。
此外,利用ModelBuilder可视化建模工具结合ArcToolbox工具箱中的多种文件转换和空间分析工具,通过构建模型实现骨架线的自动提取,其技术流程如图2所示。
图2 任意多边形骨架线提取技术流程
Fig.2 Technical processes of extracting skeleton lines from arbitrary polygons
2.2 数据获取和预处理
为便于实验数据的获取和预处理,采用电子地图为实验数据,利用高德定制地图功能结合屏幕截图获取道路的电子地图(图3(a)); 同时,采用ArcToolbox工具箱中的ENVI非监督分类模型提取出道路的矢量边界; 然后,利用高级编辑工具(Advanced Editing)的要素概化(Generalization)工具,对其边界进行抽稀处理得到道路边界的矢量多边形,模型处理结果如图3(b)所示。
图3 实验数据的获取及预处理
Fig.3 Acquisition and preprocessing of the experimental data
2.3 骨架线结点的生成
2.3.1 创建约束Delaunay三角网
采用ModelBuilder建模工具结合ArcToolbox工具中的数据处理、文件转换和空间分析工具,通过构建模型自动提取任意多边形的骨架线; 同时,为便于模型的表述,将其分成约束Delaunay三角剖分、骨架线结点生成和骨架线结点后处理3个模型分别进行设计。首先,利用添加数据工具将数据加入道路矢量多边形要素类,并依次加入要素折点线转点、添加字段、计算字段、创建TIN和TIN转三角形工具; 将每一步骤的各个工具进行连接之后,通过运行模型生成建模区域的约束Delaunay三角网格(图4(a)—(b)),设计的模型如图4(c)所示。其中,要素折点转点工具负责提取矢量多边形边界线的结点; 添加字段和计算字段工具负责给结点添加相应的高程字段。
图4 任意多边形的约束Delaunay三角剖分及其模型
Fig.4 Processes of creating constrained Delaunay triangulation of arbitrary polygons and its model
2.3.2 骨架线结点生成及后处理
同样,通过要素转线、在折点处分割线、要素转点和空间连接工具,生成多边形内部骨架线的结点。其中,要素转线负责将每个三角形面转换成线; 在折点处分割线负责将组成三角形的多线段在折点处进行分割,即一个三角形的线可分割成3条线段; 要素转点负责提取线段的中点; 空间连接工具通过搜索多边形内部的结点提取出骨架线的结点文件(图5(a)—(b)),创建的模型如图5(c)所示。这里需要注意的是,为便于后续的脚本编程和程序自动处理,在利用删除重复工具去除重复结点的基础上,采用空间连接工具建立三角面要素类与骨架线结点要素类的属性连接,对骨架线结点进行后处理。设计的模型如图6所示,生成的空间连接三角面要素类的属性表如表1所示。
图5 任意多边形的骨架线结点生成过程及其模型
Fig.5 Processes of generating the skeleton lines of arbitrary polygon and its model
图6 骨架线结点后处理模型
Fig.6 Model for post processing the skeleton line nodes
表1 空间连接三角面要素类属性表①
Tab.1 Attribute table of spatial join triangular feature class
| OID* | Shape* | ID | T1 | J1 |
|---|---|---|---|---|
| 1 | Polygon Z | 1 | 1 | 25 |
| 2 | Polygon Z | 1 | 1 | 41 |
| 3 | Polygon Z | 2 | 2 | 39 |
| 4 | Polygon Z | 2 | 2 | 41 |
| 5 | Polygon Z | 3 | 3 | 35 |
| 6 | Polygon Z | 3 | 3 | 37 |
| … | … | … | … | … |
①: OID*为对象标识字段; Shape*为图形字段; ID为图形标识字段; T1为三角形索引字段; J1为骨架线结点索引字段; Polygon Z为带Z高度的多边形。
2.4 骨架线的连接
2.4.1 骨架线几何模型的构建
多边形骨架线生成的关键在于提取Delaunay三角形内部到边界线上的等距离点集,主要通过三角形之间的空间邻接关系决定骨架线的连接关系。首先,根据前面生成的骨架线结点表和空间连接三角形表计算出每个三角形包含的结点数,并按它们包含的结点个数(1,2和3)将其连接关系分为Ⅰ,Ⅱ和Ⅲ共3种类型; 然后,对每种类型三角面分别生成不同的骨架线: Ⅰ型三角形连接重心O与对边的结点P1; Ⅱ型三角形连接2条边的中点P1和P2; Ⅲ类三角形分别连接重心O与3条边的中点P1,P2和P3。3种三角形类型骨架线连接关系如图7所示。
图7 3种三角形类型的骨架线连接
Fig.7 Three types of triangles for connecting the skeleton lines
2.4.2 骨架线连接脚本程序的实现
利用IDLE集成开发环境编写脚本程序,同时结合ArcToolbox工具箱创建脚本工具,对脚本程序进行调试,其程序实现主要包括4个步骤: ①利用创建要素类管理器(CreateFeatureclass_management)生成骨架线要素类Shapefile文件,并创建插入游标(InsertCursor); ②利用ArcPy站点包查询游标(SearchCursor)的查询功能,读取约束Delaunay三角形要素类,并返回其属性表全部记录集(表2); ③利用SearchCursor的嵌入式查询功能,通过读取空间连接三角面要素类属性表(表1),返回其相邻三角形记录集; ④利用SearchCursor的嵌入式查询功能,通过读取当前及相邻三角形包含的骨架线结点要素类属性表(表3),根据骨架线结点的数量(1,2或3),连接生成骨架线,并将其插入骨架线要素类Shapefile文件。
表2 约束Delaunay三角形要素类属性表
Tab.2 Attribute table of constraint Delaunay triangular feature class
| OID* | Shape * | ID |
|---|---|---|
| 1 | Polygon Z | 1 |
| 2 | Polygon Z | 2 |
| 3 | Polygon Z | 3 |
| 4 | Polygon Z | 4 |
| 5 | Polygon Z | 5 |
| 6 | Polygon Z | 6 |
| … | … | … |
表3 骨架线结点要素类属性表
Tab.3 Attribute table of skeleton node feature class
| OID* | Shape * | ID |
|---|---|---|
| 25 | Point Z① | 25 |
| 41 | Point Z | 41 |
| 39 | Point Z | 39 |
| 41 | Point Z | 41 |
| 35 | Point Z | 35 |
| 37 | Point Z | 37 |
| … | … | … |
①: Point Z为带Z高度的点。
骨架线连接脚本程序的输入、输出、操作步骤和主要代码实现如下。
输入: 约束Delauany三角面要素类,空间连接约束Delauany三角面要素类,骨架线结点要素类
输出: 骨架线的Shapefile文件
#输入、输出变量定义和初始化
fc0=ConversionUtils.gp.GetParameterAsText(0)
fc1=ConversionUtils.gp.GetParameterAsText(1)
fc2=ConversionUtils.gp.GetParameterAsText(2)
outfc=ConversionUtils.gp.GetParameterAsText(3)
s = '\\'
npos = outfc.index(s)
out_workspace = outfc[0: npos]
file_name= outfc[npos: ]
该脚本程序的实现流程如图8所示。
图8 骨架线连接的程序流程
Fig.8 Procedure flow of connecting the skeleton lines
步骤1: 利用创建要素类管理器(CreateFeatureclass_management)生成一个线要素类的Shapefile文件,并利用插入游标(InsertCursor)生成和返回其记录集;
步骤2: 通过查询游标(SearchCursor)访问Delauany三角面要素类并返回全部记录集;
cur0= arcpy.da.SearchCursor(fc0,[″ID″])
步骤3: 通过循环遍历约束Delaunay三角面要素类记录集的每一行记录,同时,利用查询游标(SearchCursor)结合提取的三角面标识字段并作为嵌入式查询的输入变量,访问空间连接Delauany三角面要素类并返回满足条件的三角面集合;
for row0 in cur0:
id1=str(row0[0])
cur1 = arcpy.da.SearchCursor(fc1,[″Target_FID″,″Join_FID″],' Target_FID =\'' + id1 + '\'')
步骤4: 通过循环遍历空间连接约束Delauany三角面要素类记录集的每一行记录,并利用Python数组保存每一行记录Join_FID的值;
步骤5: 根据数组的长度判断结点数及三角面的类型(I,II和 III型),同时,利用查询游标(SearchCursor)结合数组中的Join_FID的值并将其作为嵌入式查询的输入条件变量,访问骨架线结点要素类提取其坐标并分类型连接生成骨架线段。
#提取三角形重心
PntTriArray = arcpy.Array()
for part in row0[1]:
for pnt in part:
PntTriArray.add(arcpy.Point(pnt.X,pnt.Y,0,0,0))
A=PntTriArray[0]
B=PntTriArray[1]
C=PntTriArray[2]
#访问骨架线结点要素类提取其坐标
cur2= arcpy.da.SearchCursor(fc2,[″ID″,″SHAPE
@XY″,″POINT_X″,″POINT_Y″],'ID=\'' + PntId0 + '\'')
cur2=arcpy.da.SearchCursor(fc2,[″ID″,″SHAPE@XY″,″POINT_X″,″POINT_Y″],'ID=\'' + PntId1 + '\'')
cur2=arcpy.da.SearchCursor(fc2,[″ID″,″SHAPE@XY″,″POINT_X″,″POINT_Y″],'ID=\'' + PntId2 + '\'')
for row2 in cur2:
x=row2[2]
y=row2[3]
步骤6: 利用ArcPy数组保存提取的骨架线段,同时,利用插入行函数将所得到的骨架线段插入到骨架线要素类的Shapefile文件,并进行插入事务的提交。
feat1 = cur3.newRow()
feat1.shape = PntArray
cur3.insertRow(feat1)
if cur3:
del cur3
2.5 骨架线提取模型的构建
在脚本程序实现的基础上,利用ArcCatalog在ArcToolbox工具箱中创建骨架线提取脚本工具,通过设定相关的输入和输出参数,并连接之前编写和保存的脚本程序,其操作界面、输入和输出参数如图9(a)所示,提取的骨架线结果如图9(b)—(c)所示。
图9 脚本工具用户界面和骨架线提取结果
Fig.9 User interface of the script tool for extracting the skeleton lines and the extracted results
最后,利用ModelBuilder建模工具将此脚本与前面设计的模型进行连接,骨架线提取模型如图10所示。
图10 骨架线提取模型
Fig.10 Model constructed for extracting skeleton lines
3 应用测试与结果分析
3.1 应用测试
为测试本文方法的可靠性和有效性,选择淮安市的省级道路和南京市某小区的建筑物2种典型地物作为实验对象,对本文方法进行测试和验证。通过高德地图网站获取电子地图,并通过预处理分别获取道路和建筑物的矢量多边形; 同时,为便于操作和测试模型的性能,将之前分段设计的约束Delaunay三角剖分、骨架线结点生成、骨架线结点后处理和骨架线连接脚本程序连接构成一个提取模型,利用ModelBuilder建模工具分别将提取的道路路网和建筑物的矢量多边形模型添加到建模环境,分别运行模型建模并对实验结果进行分析,其骨架线提取效果如图11所示。
图11 应用测试及骨架线提取结果
Fig.11 Application test and the extracted results of skeleton lines
3.2 结果分析
主要从模型的运行时间和骨架线提取的效果2方面进行分析(表4)。
表4 骨架线提取实验分组统计
Tab.4 Experimental extraction of skeleton lines by grouping statistics
| 多边形对象 | 约束Delaunay 三角面数量/个 | 骨架线提取模 型运行时间/s |
|---|---|---|
| 道路 | 475 | 60 |
| 建筑物 | 735 | 80 |
通过表4发现骨架线提取模型的运行时间与约束Delaunay三角剖分获得三角形的数量相关,并且二者运行的时间差异不明显。这可能是由于要素折点转点、要素转点、添加字段、字段计算器和空间连接等操作也消耗了一定的计算时间,导致模型的运行时间差异不是很明显。
在骨架线提取效果方面,通过对道路和建筑物提取的骨架线进行局部放大观察和分析,道路多边形骨架线提取的效果要优于建筑物多边形骨架线提取的效果。一是由于经过预处理生成的建筑物多边形,其内部会存在一定数量碎块,骨架线连接错乱(图12(a)); 二是由于经过简化处理生成的建筑物多边形,其边界结点数据量太少,提取的骨架线过于简化(图12(b))。针对上述2种情况,可利用编辑工具和加密工具删除小于一定面积的碎块多边形结点并对其边界进行加密,再利用模型提取其骨架线,其处理结果如图12(c)—(d)所示。
图12 骨架线提取实验结果对比分析
Fig.12 Comparing analysis for the experimental results of the skeleton lines extraction
4 结论与讨论
近年来,以任意多边形为研究对象已提出多种多边形骨架线提取方法,多局限于基于栅格图像和矢量图形的计算,需要繁琐的算法设计和复杂的代码实现,实现成本高、周期长。与传统方法相比,本文从地理设计的视角,提出的任意多边形骨架线提取方法具有以下优势:
1)所有与骨架线提取相关的数据预处理、编程开发和建模工作都在同一种GIS平台下完成; 由于采用Phyton面向对象编程语言结合ArcPy站点包编程编写脚本程序,更加易于程序的编写和维护,同时提高了开发效率。
2)利用ModelBuilder工具,结合ArcToolbox工具箱中的系统工具和二次开发脚本工具,通过构建模型实现任意多边形骨架线的提取,其操作过程更加符合熟悉GIS用户的操作习惯,具有更强的可操作性和实用性,并且更易于模型的调试、维护和功能扩展。
本文提出基于GIS空间分析的任意多边形骨架线自动提取方法,一方面,通过分析骨架线的几何形态和空间分布规律,为导航系统中道路的自动化提取和路网构建奠定了技术基础; 另一方面,为遥感数字图像处理中的建筑物目标的自动识别、特征提取和建筑物三维重建提供了有效的技术支持。此外,通过GIS建模的方法实现了任意多边形骨架线的自动化,也为开发具有自主知识版权的自动矢量化系统提供了研发技术参考。文中所述方法可对二维空间内任意多边形的骨架线的几何形态和空间分布规律进行探讨,但由于本文方法需要借助GIS系统的数据处理、空间分析和文件转换工具,模型运行会产生大量中间文件,势必会对模型的性能产生一定的影响。对此,笔者认为可采用固态硬盘技术加速中间文件的读写速度,以降低这一影响。此外,本文方法对骨架线形态的优化未做探讨,这将是未来研究的重点。
参考文献略
END
排版：喜马拉雅
审核：晨风小语