李功权，蔡祥云 (长江大学地球科学学院，湖北 武汉 430100)

一种基于约束三角网的道路中心线的提取方法

李功权，蔡祥云 (长江大学地球科学学院，湖北 武汉 430100)

鉴于道路中心线应用的广泛性，研究了基于约束Delaunay三角网的道路中心线的提取算法。以道路边界线作为约束线，采用Delaunay方法构建三角网。通过确定相邻三角形的类型，把获取的节点分为3类，其对应道路网络中的十字路、T型路和环岛路，对其分别进行优化处理，从而形成道路的中心线。在给出详细的算法步骤的同时，并用C#语言实现该算法。实测数据应用分析表明，该算法生成的道路中心线符合原道路多边形的形态，保持了原图形的拓扑特征。

道路中心线；约束Delaunay三角网；道路网络模型

根据道路中心线建立的道路网络模型，充分利用了GIS的网络分析功能，在交通管理和汽车导航中有着广泛的应用。道路中心线的数据一般无法直接得到，在道路规划设计中可能使用的是道路两边的边界线，城市土地管理部门可能有各个地块的边界信息，这就需要从这些数据中提取道路中心线，而不是去实地人工采集。如果把需要求取道路中心线的区域看作一个多边形，不是道路的相关空间信息作为约束信息，道路中心线的求取问题就可以转化为多边形骨架线的求取。所谓骨架线，就是用与原形状连通性和拓扑结构相一致的细曲线作为原对象的一种抽象表示，多边形骨架线的本质是对多边形形状的抽象描述，它反映了多边形的延伸方向和形状特征[1-3]。在GIS中，面状地理空间对象的分析等很多空间操作都需要提取骨架线[4-6]。

1 算法原理

多边形骨架线提取的关键是搜寻多边形内部到边界线上的等距离点集，本质上属于空间邻近分析问题。一般的骨架线的提取算法有：①数学形态学提取骨架线，这种方法本质是矢量化方法；②最大内切圆盘法，最大圆盘完全落于目标图像内，并且至少有2点与目标边界相切。骨架的每一个点都对应于一个最大圆盘的圆心和半径，圆盘的构建特别是小圆盘的构建是该算法的最大难题；③基于Delaunay三角网的多边形骨架线提取算法，Delaunay 三角网是一系列相连但不重叠的三角形的集合，而且这些三角形的外接圆不包含面域中其他任意点，且是Voronoi图的对偶[7]。Delaunay三角剖分可以最大限度地避免狭长三角形的出现，并且可以不管何处开始都能保持三角网络的唯一性，Delaunay三角网是探测空间图形邻近关系的优秀工具。因此，笔者采用Delaunay三角网作为提取骨架线的理论工具。以多边形的边界为约束条件构建三角网，取三角形边线的中点作为骨架线的节点，顺次连接这些节点，得到多边形的骨架线[7]。

1.1约束Delaunay三角网的构建

对于约束Delaunay三角网生成算法，大致可以分为3种：分治算法、逐点插入法和三角网生长法。而逐点插入算法的特点是实现比较简单，占用内存小，因此笔者采用逐点插入法生成无约束的Delaunay三角网，再根据约束边删除多边形外部多余的三角形。具体过程如下：

1)将离散后多边形的顶点，建立一个包含其他数据点的初始多边形，称其为凸包；

图1 约束Delaunay三角行的生成过程

2)在初始多边形中建立初始三角网，对所有初始多边形中数据点循环处理(见图1(a))。

3)插入一个数据点P，在已有三角网中找出包含P的三角形T，把P与T的3个顶点相连，生成

3个新的三角形，用LOP算法优化三角网(见图1(b))。

4)删除不在多边形内部的三角形。判断三角形的一边是否在多边形的内部，如果在其内部保留该边，如果不在则舍弃。具体的实现过程是每次选取一个三角形一边的中点，从该点根据射线法进行判定，最后结果如图1(c)所示。

1.2三角形类型的确定

图2 三角形类型的划分

三角形类型的确定是在Delaunay三角剖分的同时确定三角形的类型。在该过程中，还要标记出新生成的三角形为何种类型，目的是用来识别道路中心线节点的类型。从多边形内部三角形的邻近关系来看，可以分为3种类型的三角形：第Ⅰ类三角形是只有1个邻接三角形；第Ⅱ类三角形是有2个邻接三角形；第Ⅲ类三角形是3条边都有邻接三角形。根据邻近三角形的数目，将三角形分为3类(见图2)：第Ⅰ类三角形是三角网中的边界节点，其中一边的中点作为骨架线的端点；第Ⅱ类三角形是三角网中的桥接三角形，是道路中心线的骨干结构，描述了中心线的延展方向；第Ⅲ类三角形作为中心线分支的交汇处，是向3个方向伸展的出发点。

三角形类型的确定方法主要依据与某个三角形相邻三角形的个数。首先统计分别与三角形的三条边相邻三角形的个数总和，默认的情况下设置与三角形的一条边的邻接三角形的个数为0；其次，根据三角形三边相邻三角形的总和判断三角形的类型，当值为1时就是第Ⅰ类三角形，值为2时就是第Ⅱ类三角形，值为3时就是第Ⅲ类三角形。这样就能判断出三角形的类型。

1.3中心线的提取

首先，判断三角形是否是多余三角形，如果不是就判断它是哪种类型的三角形。如果是第Ⅰ类的三角形，提取桥接边的中点和另外2边中较长的一边的中点；第Ⅱ类的三角形提取2个桥接边的中点；第Ⅲ类三角形则需要提取该三角形的重心和3条桥接边的中点。其次，对于第Ⅰ类和第Ⅱ类的三角形来说提取的2个点便是中心线；对于第Ⅲ类三角形来说，将重心分别和其他的3个点相连便也是中心线。但求出的中点是事先需知道该三角形的哪条边是桥接边，这样便于找出桥接点、端点、分支点，从而确定中心线各节点的类型。这样把Delaunay三角形一个个的单独处理，提取他们各自的这些端点、桥接点、分支点连起来，便可获得最终的结果(见图3)。

图3 中心线的提取

2 算法实现

根据道路多边形提取道路中心线，涉及到2个基本的数据结构。一是三角形的定义，作为三角网的构成基本单元，不仅需要表达自身三角网的特征，还需要考虑三角形之间的关系；二是道路中心线的表达，可以看成是一系列具有特定属性的顶点的集合。

三角形的数据结构可以用一个结构体来定义，具体的数据结构如下：

struct Triangle

{

public long V1Index; //三角形的三个顶点

public long V2Index;

public long V3Index;

public bool edge1; //(v1,v2)是否已经存在

public bool edge2; //(v2,v3)是否已经存在

public bool edge3; //(v1.v3)是否已经存在

public int AdjIndexE1; //edge1的邻近三角形的索引

public int AdjIndexE2; //edge2的邻近三角形的索引

public int AdjIndexE3; //edge3的邻近三角形的索引

public bool bDelete ; //判断多余的Delaunay三角形是否被删

public bool Fkind; //第Ⅰ类三角形

public bool Skind; //第Ⅱ类三角形

public bool Tkind; //第Ⅲ类三角形

}

其中，AdjIndexE1、AdjIndexE2、AdjIndexE3分别用来存储(v1，v2)边、(v2，v3)边、(v1，v3)边的邻近三角形的索引，默认赋值为-1；Fkind、Skind、Tkind 3个bool变量分别用来判断该三角形属于那种类型，Fkind代表第Ⅰ类三角形，Skind代表第Ⅱ类三角形，Tkind代表第Ⅲ类三角形，默认值都为false。

根据相邻三角形的类型可以把形成的中心线节点分为端点、桥接点和分支点，具体的数据结构如下：

struct Vertex

{

public long x; //顶点的x坐标

public long y; //顶点的y坐标

public long ID; //顶点的索引

public int isHullEdge; //凸壳顶点标记

}

在定义了这2个数据结构后，用C#语言根据上述原理实现了道路多边形文件的读取、道路边界点加密、约束三角网的构建、三角形中心点的提取、道路中心线各节点的优化处理5个模块。

3 应用效果评价

如果道路线上的边界点过于稀疏，在生成的三角网中容易出现小内角的狭长三角形。为了解决该问题，可以采用对道路多边形边界点作适当加密的办法，首先对各边界点作3次样条函数拟合，然后把道路的平均宽度作为间距加到边界点中。

为了验证该算法，把同样的数据用ArcGIS 9进行了处理。在ArcGIS中，如果要提取道路的中心线，如果是面要素，要先将面要素转换为线要素(此处的线要素是双线要素)，然后利用制图工具中的提取中心线工具，设定最大值和最小值参数，且必须给定最大值。设计的测试数据不仅考虑了实际数据的特点，而且还设计了多种道路类型。ArcGIS提取的结果如图4所示，虽然在多数区域提取的道路中心线符合实际，但在图4中圆圈和矩形所圈定的区域都存在着明显错误，圆圈所指示的道路中心线明显不存在；对于某些交叉路口(矩形所指示的范围)，和实际的道路网络模型不符。笔者设计的算法提取的道路中心线如图5所示，对于任何类型的道路限制条件都能适用。

图4 ArcGIS提取的道路中心线 图5 笔者设计算法提取的道路中心线

4 结 语

由于形成道路的多边形样式变化多端，道路中心线的提取是一个复杂的过程。 笔者运用约束的Delaunay 三角网法提取道路中心线，算法完全基于矢量数据结构实现，算法建立在多边形的形状分析基础上，高效稳定。从试验结果来看，提取的道路中心线能反映出道路的基本特征。可见笔者设计的算法可以用来求取道路中心线。但该算法可能产生不需要的细长三角形，多边形边界的细微波动可能会导致道路中心线的明显偏移，这个问题值得进一步研究。

[1]Shaked D,Bruckstein A M. The curve axis [J].Computer Vision and Image Understanding, 1996, 63(2):367-369.

[2] LI Z L. Algorithmic Foundation of Multi-scale Spatial Representation [M].CRC Press, 2007：20-23.

[3]Attneav E F. Some informational aspects of visual perception [J].Psychological Review, 1954, 61(3):183-193.

[4] Mcmaster R B. A statistical analysis of mathematical measures for line simplification[J]. The American Cartographer, 1986, 13: 103-116.

[5] Mcmaster R B. Automated line generalization [J]. Cartographica, 1987, 24(2): 74-111.

[6] LI Z L. An examination of algorithms for detection of critical points on digital lines [J]. The Cartographic Journal, 1995, 32(2): 121-125.

[7]Haunert J H,Sester M. Area Collapse and Road Centerlines based on Straight Skeletons [J]. Geoinformatica,2008,12:169-191.

[8] Paul C L. Constrained Delaunay Triangulations [J]. Algorithmica, 1989, 4:97-108.

2012-11-24

中国石油科技创新基金(2010D-5006-0205)

李功权(1971-)，男，博士，副教授，现主要从事GIS、数字油藏方面的教学与研究工作。

TP391.4

A

1673-1409(2013)04-0047-04

[编辑] 洪云飞