在地图制图和GIS领域，骨架线是用于描述地图目标内部或地图目标群间的拓扑结构特征的主要手段之一。相关应用涉及地图自动综合、地图注记的自动配置、栅格数据的压缩存储等众多方面。提取骨架线的算法大致分为基于栅格的方法和基于矢量的方法两类^[1]，基于栅格的方法主要包括传统的细化算法^[2]和基于数学形态学、地图代数的算法^{[1, 3]}；基于矢量的方法主要有Delaunay三角网法^[4]和Voronoi图法等^[5]。基于栅格的方法能够较准确地获取研究区域的中轴，避免产生过多“毛刺”，但需要事先对图形进行栅格化、二值化等，增加了计算量。与之相比，基于矢量的方法更易于与面向目标的地图数据模型无缝结合，且效率高。由于各种商用GIS软件中均提供了对TIN数据结构的支持，使得基于Delaunay三角网的骨架线的提取和应用十分便利。

目前，制图自动综合中基于Delaunay三角网的骨架线的应用已有很多，包括线状目标的弯曲识别及简化^[6]、双线河流或街道的中轴线提取及降维^{[7, 8]}、河网汇水区域的剖分与综合^[9]、道路网的空间冲突识别与分类^[10]、多边形(群)的简化、合并和移位^{[11, 12]} 等。在艾廷华提出的基于Delaunay 三角网的空间场表达的形式化数据模型(formal triangulation data model,FTDM)中，明确将骨架线提取(骨架化)定义为模型的四种基本操作之一，并给出了概念化定义^[13]。该算法的实现一般有质心法、中线中点法和非约束边中点法几种方式^[8]。以上研究多以单个地图目标内部骨架线的提取为主，缺乏对地图目标群或复杂地图目标内部进行骨架线提取时的特殊情况的考虑。本文结合算法的详细设计和研究，对基于约束Delaunay三角网的骨架线提取算法进行了改进，试图处理一些特殊情况，增强算法的稳健性。

1 算法基本思路

首先，要建立平面约束Delaunay三角网。Delaunay三角网算法已是计算几何领域的成熟的算法，此处不作讨论。

提取骨架线时，按照三角网中三角形之间的邻接关系追踪覆盖地图目标之间空白区域的三角形序列，将这些三角形的中轴点(顶点、重心或虚边上的中点)按序串连起来，从而形成树状或网状的骨架线结构。本文将三角网中覆盖一对地图目标之间空白区域的三角形序列称为通道；通道中提取的中轴线称为骨架线弧段(简称弧段)；通道的第一个三角形称为入口三角形，最后一个三角形称为出口三角形，它们统称为结点三角形。按照三角形包含约束边的数目，将三角网中的三角形分为0、1、2和3四类(见图 1)：0类没有约束边；1类包含1条约束边；2类包含2条约束边；3类包含3条约束边。为叙述方便，将三角形的约束边称为实边，非约束边称为虚边。三角网中相邻的三角形由公共虚边连通，被公共的实边隔离。

追踪过程中，对不同类型的三角形采用不同的处理方式。如图 1所示，一个0类三角形的3条虚边可分别通向3条通道，将其重心作为3条对应弧段的交叉结点分别与3条虚边的中点相连，作为3条弧段的起始线段或终止线段；一个1类三角形的2条虚边可分别连通2个相邻的三角形，该三角形一般处于一条通道的中部，将其两条虚边的中点相连，构成所属弧段中的一条线段；一个2类三角形的唯一虚边仅可与一个相邻的三角形连通，其必为一条通道的入口或出口三角形，将其虚边的中点和与之相对的顶点相连，构成对应弧段的起始线段或终止线段；3类三角形比较特殊，其三条实边将其与相邻三角形完全隔离，处理时仅需将三角形重心分别与三个顶点连接构成三条简单弧段即可。以上分析为一般情况，对于三角网边界上的三角形可能并不成立，这时还需特殊处理。例如，在三角网边界上的1类三角形可能有一条虚边在边界上，此时仅剩另一条虚边与一个相邻的三角形连通(单连通)，因此，也可作为一条通道的入口或出口。此时，将其两虚边的中点相连，构成弧段的起始线段或终止线段。文中将此类三角形称为1′类三角形。

图 1 四类三角形及各自内部的骨架线段 Fig. 1 Four Types of Triangles and the Corresponding Skeleton Segments

图选项

经以上分析可知，0类、2类和1′类三角形均可作为一条通道的入口或出口三角形，为弧段提供端点。不同的是，2类或1′类三角形仅能为一条弧段提供端点，而0类三角形可同时为3条弧段提供端点。

但当三角网中存在首尾相连的双连通1类三角形构成环路时，由于环路上没有可作为入口的0类、2类以及1′类三角形，用常规方法无法提取其骨架线(见图 2中的“环路”示例)。此时，可将环路中的任意一个三角形作为入口，选择其两条虚边中的任意一条追踪相邻的三角形，沿着选定的方向一直追踪，直到回到入口三角形为止。

2 算法设计

首先，给出相关数据结构和概念的形式化定义。

1) 平面三角网图结构。由非空点集合 V={v₁,v₂,…,v_k}、非空边集合E={e₁,e₂,…,e_m}和非空三角形集合T={t₁,t₂,…,t_n}构成三元组图结构G_t(V,E,T)。

2) 骨架线图结构。由非空骨架线结点(简称结点)集合O={o₁,o₂,…,o_p}、非空骨架线弧段(简称弧段)集合A={a₁,a₂,…,a_q}构成二元组图结构G_s(O,A)。

3) 出入口信息。追踪骨架线的过程中，进入或离开一个通道时所穿过的结点三角形中的虚边记为(e,t),其中e为虚边，t为该虚边所属的结点三角形。

本文提取骨架线的过程就是对隐含在G_t中的G_s进行遍历，并记录下所经过的各弧段和结点的过程。主要函数和操作定义如下。

1) 三角形类型函数Type(t)。输入一个三角形t，返回t的类型，结果标识为 0、1、1′、2或3，分别代表 0类三角形、1类三角形、单连通1类三角形、2类三角形或3类三角形(虽然前文中仅有四类三角形，但此处从算法设计角度考虑，将1类中的单连通和双连通三角形作为1和1′两类，因此Type(t)的返回值有5种)。

2) 邻接三角形函数Neighbor( e,t) 。输入一个三角形t及其包含的一条边e ，返回与t相邻且以e为公共边的相邻三角形t′。

3) 通道入口属于0类三角形时的操作P₀^s(e,t,&a,&S)。输入入口信息(e,t)、对应弧段a、入口信息栈S，输出t的重心作为a的起点。本操作将t的重心和e的中点作为a的第1、2点加入弧段中，如果t未被访问过，则将t中另外两条边对应的入口信息加入栈S。

4) 通道出口属于0类三角形时的操作P₀^e(e,t,&a,&S)。输入出口信息(e,t)、对应弧段a、入口信息栈S，输出t的重心作为a的终点。本操作将t的重心和e的中点作为a的倒数第1、2点加入弧段中，如果t未被访问过，则将t中另外两条边对应的入口信息加入栈S，否则在S中查找 (e,t)，如果找到，将其删除。

5) 1类三角形操作P₁(e,t,&a)。输入1类三角形t及其虚边e。本操作将e的中点加入弧段a中。

6) 通道入口属于1′类三角形时的操作P_1′^s(e,t,& a)。输入入口信息(e,t)和对应弧段a，输出t的不连通虚边的中点作为a的起点。本操作将t的不连通虚边的中点和e的中点作为a的第1、2点加入弧段。

7) 通道出口属于1′类三角形时的操作P_1′^e(e,t,&a)。输入出口信息(e,t)和对应弧段a，输出t的不连通虚边的中点作为a的终点。本操作将t的不连通虚边的中点和e的中点作为a的倒数第1、2点加入弧段。

8) 通道入口属于2类三角形时的操作P₂^s(e,t,&a)。输入入口信息(e,t)和对应弧段a，输出t的虚边e相对的顶点作为a的起点。本操作将e相对的顶点和e的中点作为a的第1、2点加入弧段。

9) 通道出口属于2类三角形时的操作P₂^e(e,t,&a)。输入出口信息(e,t)和对应弧段a，输出t的虚边e相对的顶点作为a的终点。本操作将e相对的顶点和e的中点作为a的倒数第1、2点加入弧段。

10) 3类三角形P₃(t,G_s)。分别将t的3个顶点与其重心连接，形成3条弧段加入骨架线结构图G_s中。

算法伪代码如下：

Algorithm Skeleton

输入：三角网图结构 G_t(V,E,T)

输出：骨架线图结构G_s(O,A)

S←，i←0

while S ≠ 或 i < n(n为三角形的总数)

if S = then t←t_i，i++

if t为未被访问过的结点三角形then

a←创建弧段

switch(Type(t))

case 0： e←t的第1条边, o←P₀^s(e,t,&a,&S) case 1′： e←t的连通虚边, o←P_1′^s(e,t,&a)

case 2： e←t的虚边, o←P₂^s(e,t,&a)

end switch

else if Type(t)=3 then

P₃(t,G_s),continue

end if

else then

(e,t)←Pop(S),a←创建弧段,

o←P₀^s(e,t,&a,&S)

end if

将o加入G_s

t←Neighbor(e,t)

whilet≠NULL 且 Type(t)=1

e←t的另一条虚边,P₁(e,t,&a),

t← Neighbor(e,t)

end while

switch(Type(t))：

case 0： o←P₀^e(e,t,&a,&S)

case 1′： o←P_1′^e(e,t,&a)

case 2： o←P₂^e(e,t,&a)

end switch

将a和o加入G_s

end while

while 三角网中还有未追踪的三角形

a←创建弧段

t←任取一未追踪的1类三角形，

e←t的任一条虚边,o←e的中点

将o作为弧段起点加入a

t_s←t，t← Neighbor(e,t)

while t ≠t_s

e←t的另一条虚边,P₁(e,t,&a)

t← Neighbor(e,t)

end while

将a加入G_s，并再次将o加入a end while

3 应用与分析

为了验证该算法的可行性和健壮性，分别以等高线的内插和街区地图目标群邻近分析为例，对算法进行验证与分析。

图 2(a)为采用本文算法提取的等高线簇间的骨架线，算法成功地将图中的一个三角形环路识别，并提取出对应的环形弧段。所得骨架线图中不仅包含了相邻等高线之间的中轴线，而且包含了单条等高线的内部弯曲或封闭区域内的中轴线。将这些内部中轴线删除，然后对剩余的骨架线进行拓扑重建，即可得到内插等高线(见图 2(c))。如果要得到更好的效果，可对内插所得的等高线进行优化，如删除由于0类三角形产生的局部“抖动”，对线划进行光滑，以及相邻线对之间的进一步协调等。

图 2 等高线内插 Fig. 2 Contour Interpolation

图选项

图 3描述了对地图上部分街区的空白区域提取骨架线并构建地图目标之间邻近图的基本过程。图中骨架线近似地对地图目标之间的间隔区域进行等距离剖分。每个目标周围由数条骨架线弧段环绕构成一个闭合区域，该区域表达了对应目标的空间势力范围，总体上形成一种类似Voronoi 图的剖分结构。图中每条弧段对应相邻两个目标之间的一阶邻近关系。基于这种特性，图 3绘出了对应的邻近图，显然两者存在对偶性，都可为进一步的空间分析与地图综合操作提供丰富的空间邻域知识。

图 3 街区骨架线与邻近图 Fig. 3 Skeleton and Proximity Graph of District

图选项

本文改进了一种基于约束Delaunay三角网的地图目标群间骨架线提取算法。在深入研究算法的数据结构和控制流程的基础上，对三角网中的环路和3类三角形的特殊情况进行了特殊处理，使算法更具稳健性。分别以地图中等高线的内插和街区地图目标群的邻近分析为例，对算法进行了验证。下一步将继续研究如何将该算法用于道路网、居民地等地图目标之间的空间冲突的检测与解决。