2. 中国科学院 地理科学与资源研究所，北京 100101
2. Institute of Geographical Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101,China
中轴(medial axis，MA)又称骨架，是对空间形体结构的一种描述，表达区域相连的几何形体的基本特征。中轴不仅应用于医学图像分析[1]、海洋测绘、数学形态学[2]等领域，还被运用于GIS中进行空间分析，如地形分析、制图综合[3]、DEM建模[4]等。
中轴变换方法首先是由文献[5]提出，此后有不少学者对中轴问题进行深入、广泛地研究。当前中轴提取方法可以分为两类：一类是基于连续几何模型(拓扑形状分析)的方法[6, 7]，目前广泛应用的有基于Delaunay三角剖分的外心法[8]、重心法[9]和内心法[10]，这些方法基于矢量模型，所得到的中轴从本质上并不符合中轴的概念[11],同时在处理复杂矢量图形上具有较大难度，不能处理含岛或洞的多边形及边中含有由自曲线的多边形，不是通用的中轴提取方法；另一类是基于离散图像(像素点)的传统方法，主要有细化迭代的方法[12]、基于距离变换的方法两种[13],细化迭代方法所得到的中轴具有良好的拓扑不变性，但是中轴位置不准确[14, 15]，基于距离变换的方法，易于实现，中轴位置准确，这种基于离散域的算法可以处理复杂图形，但是中轴容易受到边界噪声的干扰，容易出现毛刺并且中轴的连通性难以得到保证[16]。
针对上述方法的不足，本文提出一种基于改进距离变换的种子点生长判别法，用于抑制边界噪声，提取连通、准确和光滑的曲边多边形中轴。
1 曲边多边形中轴的定义简单多边形的边均为直线段，但是在实际应用中，多边形的边中不仅包括直线段，还含圆弧、自由曲线(图 1(a)、图 1(b))，内部还可能含不规则形状的“洞”、“岛”，甚至出现自由曲线、圆弧、“洞”、“岛”等任意组合(图 1(c))的复杂情况。
|
| 图 1 曲边多边形 Fig. 1 Curved-polygons |
设E为平面连通(包括单连通和多连通)区域的轮廓边界，边界E由n个顶点{P1,P2,…,Pn}和联结顶点之间的n条线段{P1P2,P2P3,…,Pn－1Pn}组成，这些线段可以是直线段，也可以是圆弧或自由曲线，则称E为曲边多边形。将组成曲边多边形E的顶点和线段称为基本元素，记由顶点元素和线段元素组成的基本元素集为{e1,e2,…,e2n}。由此定义易知：曲边多边形不仅包括一般简单多边形，还涵盖岛洞多边形及边中含有曲线段的特殊简单多边形。
1.2 曲边多边形中轴不同学者对多边形中轴给出了不同定义。文献[17]给出的定义是：到两个或两个以上的基本元素(顶点和边)距离相等的点的轨迹为多边形中轴[17, 18]。文献[11]提出：多边形P的中轴为E内的点集，该点集中的点与∂P不同边或多边形边的延长线中两个或两个以上的点距离相等[18]。
虽然这两个定义应用较普遍，但存在以下问题：
(1) 当P为凹多边形时，在凹顶点处，根据文献[11]的定义计算出的中轴并不是多边形的中轴，而是相对于凹顶点处多边形边延长线的中轴。
(2) 以上中轴的定义是针对一般意义上的简单多边形，在处理曲边多边形时，具有明显的局限性，例如文献[19]中提及的特殊曲边多边形，如图 2所示，根本难以判断出该图形有哪几条边。
|
| 图 2 特殊曲边多边形 Fig. 2 Special curved-polygon |
为克服以上问题，本文提出曲边多边形中轴的定义，引入最近边缘点和中轴点的概念。
最近边缘点：设g为区域D内的任一点，点h为轮廓边界E上的任一点，记点h到点g的距离为d(h,g)，记点g到轮廓边界E的距离(即点g到轮廓边界E上所有点距离的最小值)为d(g,E)，若d(h,g)=d(g,E)，则称点h为点g的最近边缘点。
中轴点：若点g在轮廓边界E上具有两个或两个以上不相同的最近边缘点h1,h2,…,hn(n≥2)，则称点g为中轴点。
则曲边多边形的中轴为区域内所有中轴点的集合。其数学表达形式如下
式中，MA(E)为多边形E的中轴；d(h,g)为点h到点g的距离；d(g,E)为点g到轮廓边界E的距离；h1,h2,…,hn为轮廓边界E上不同的点；n为大于等于2的自然数，{·}表示所有中轴点的集合。
本文给出的多边形中轴定义适用于凸多边形、凹多边形、岛洞多边形、含有曲线段的简单多边形以及如图 2(可将其视为由三条自由曲线围成的曲边多边形)所示的特殊图形，在广泛的自然图形情形下也具有普适性。 2 栅格欧氏距离变换
由于欧氏距离值需要开平方，如果采用距离平方值代替距离值参与运算，这将减少计算量并极大提高距离值精度[20]。将距离平方值记为de2，每个栅格单元八邻域栅格单元的de2值按图 3所示依次记为de21,de22,…,de28。
|
| 图 3 de 2八邻域示意图 Fig. 3 Eight directions of de2 neighborhoods |
栅格欧氏距离变换的具体步骤如下：
(1) 赋所有空间点为足够大的正数M1，并赋所有实体点为零;
(2) 顺序访问，按下式改写各点平方值
de2(i,j)=min(de21(i,j),de22(i,j),de23(i,j),de24(i,j),de2(i,j))
(3) 逆序访问并改写各点平方值
de2(i,j)=min(de25(i,j),de26(i,j),de27(i,j),de28(i,j),de2(i,j))
(4) 改写各点距离平方值为距离值
c(i,j)=int{[de2(i,j)]1/2+0.5}
提取多边形中轴时需要进行内距变换，所谓内距变换是指将图形的区域作为空间点集，将空间区域作为实体点集，计算并记录图形点集中任一点到非图形点集最小距离值的变换[20]。图 4是不同曲边多边形栅格欧氏内距变换结果图。
|
| 图 4 曲边多边形栅格欧氏内距变换结果图 Fig. 4 Inner DT results of curved-polygons |
为抑制图像边界噪声的干扰，本文提出最近边缘点集距离均值变换基本步骤如下：
(1) 栅格欧氏距离变换，对曲边多边形进行栅格欧氏距离变换，获取初始距离值d(g,E)和相应的初始距离源点SrcPt(g)；
(2) 求取最近边缘点集R(SrcPt(g))，对区域D内任一点g，以已求得的距离源点SrcPt(g)为圆心，设定以r为半径的圆形区域，求取圆形区域与曲边多边形边界点集E的交集作为最近边缘点集R(SrcPt(g))；
(3) 计算距离均值d(g,E)，计算出最近边缘点集R(SrcPt(g))内任意点到g点距离值，从中去除一个最大值和一个最小值，将剩余的距离值取算术平均值得到距离均值d(g,E)；
(4) 计算新的最近边缘点SrcPt(g)，从最近边缘点集R(SrcPt(g))中挑选出到g点距离最为接近于距离均值d(g,E)的点作为新的最近边缘点SrcPt(g)；
(5) 更新距离值和距离源点，分别用距离均值d(g,E)和最近边缘点SrcPt(g)替代g点原有的距离值d(g,E)和最近边缘点SrcPt(g)。
其中，步骤(2)中的半径r值是点g的距离值d(g,E)线性函数。g点距离值越小，表示g点越靠近多边形边界，其距离值受边界噪声干扰越严重，反之亦然。经过反复试验，检验到半径r值与g点的距离值平方值de2(g,E)的9组对应数据(因为任一点g的距离值都以其平方值的形式表示)，并计算出g点的距离值d(g,E)。通过这9组试验值建立通过距离de2(g,E)计算半径r值的回归模型，如表 1所示。
| 序号 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| de2(g,E) | 16 | 35 | 100 | 150 | 240 | 300 | 450 | 570 | 752 |
| d(g,E) | 4 | 5.91 | 10 | 12.24 | 15.49 | 17.32 | 21.21 | 23.87 | 27.42 |
| r | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
对以上数据进行回归分析，得出半径r与de2(g,E)的乘幂回归方程(如式(2))，相关系数R2=0.994 2
得出半径r与d(g,E)的线性回归方程(如式(3))，相关系数R2=0.996 3
因此，半径r值可由式(2)计算并取整获得。但当d(g,E)足够大时，半径r值不宜大于10，否则所求得的最近边缘点集中将会包含过多的噪声点。
对原始多边形进行传统栅格欧氏内距变换后，由于栅格数据本身的格式及数据质量等问题，由变换得出的距离值受区域边界噪声影响较大，变换结果反而不准确。图 5中的图 5(a)和图 5(b)分别为无噪声的多边形和加入噪声的多边形，对图 5(a)进行传统的欧氏内距变换，结果如图 5(c)所示。对图 5(b)分别进行传统的欧氏内距变换和最近边缘点集距离均值变换，其结果分别如图 5(d)和图 5(e)所示。对比图 5(c)和图 5(d)，前者距离波为平行线族，后者距离波随噪声分布呈“波浪”传播，由此可知，传统内距变换方法在处理边缘上存在噪声的图形时具有局限性。而图 5(e)中结果与图 5(c)中无噪声多边形的内距变换结果基本一致，表明最近边缘点集均值变换能较好抑制边缘噪声的影响。
|
| 图 5 传统内距变换与距离均值变换的对比图 Fig. 5 Comparison between traditional DT and MDT |
经过两次距离变换后，对于区域D内任意栅格单元g，与其八邻域栅格单元{g1,g2,g3,g4,g5,g6,g7,g8}有如下关系：
(1) g与g1,g2,…,g8均已计算出到轮廓边界的最小距离值及距离源点(最近边缘点)；
(2) g与g1,g2,…,g8具有相近的空间位置，当栅格数据分辨率足够高，栅格单元尺寸足够小时，可以忽略g与其邻域栅格单元的空间距离，g与其邻域单元可以被认为是同一点；
(3) 同理，g与g1,g2,…,g8的最小距离值在数值上十分相近，最小距离值大小上的差异可以忽略不计，可以近似地认为g与g1,g2,…,g8最小距离值相等；
(4) 若栅格单元g是多边形的中轴点，根据多边形中轴的连通性，则在g的八邻域栅格单元g1,g2,…,g8中，必定至少存在一个中轴点。
因此，对于区域D内的任意点g，要在边界上再寻找一个不同于其已有距离源点的最近边缘点，只需要考虑从其八邻域单元的距离源点中判别筛选即可，筛选条件为
式中，SrcPt(gi)和SrcPt(g)分别为点gi和点g的距离源点，即最近边缘点，gi∈{g1,g2,g3,g4,g5,g6,g7,g8}。
换言之，对于区域D内的任一点g，若gi∈g1,g2,g3,g4,g5,g6,g7,g8，使得SrcPt(gi)≠SrcPt(g)，则点g即为多边形中轴点。同时，如果已经确定多边形内某一点g是多边形的中轴点，则根据多边形中轴的连通性，以g为中心，在其8邻域单元中必定找到其他中轴点，从而保证中轴提取的连续性。
4.2 实现步骤距离峰值点是指距离值不小于其8邻域栅格单元距离值的栅格单元点。
根据上述思想，新方法的基本步骤为：首先对曲边多边形区域进行精密欧氏内距变换和最近边缘点集距离均值变换，再在区域内搜寻一个或几个距离峰值点作为初始的中轴点，并以找到初始中轴点为种子点，依据上述中轴点判别方法搜索出新的中轴点，然后以新搜索到的中轴点为种子点，反复这一过程直到不能搜索到新的中轴点为止，所有搜索到的中轴点即构成曲边多边形的中轴。具体实现流程如图 6所示。
|
| 图 6 提取曲边多边形中轴的具体流程 Fig. 6 Detailed flow chart of curved-polygon MA |
上述流程通过提取距离峰值点作为初始中轴种子点，再由种子点生长出曲边多边形的其他中轴，保证所得中轴的连通性。
5 实例验证与对比分析将本文提出的新方法应用于含圆弧(图 1(a))、自由曲线(图 1(b))、洞岛、自由曲线(图 1(c))图形的中轴提取，结果如图 7所示。此外，对简单多边形、岛洞自由曲线界多边形、自然图形也进行验证，图 8、图 9、图 10分别展示了其原始图形、提取过程和提取结果[14, 15, 19]。同时也对道路网和河网进行中轴提取，结果分别如图 11、图 12所示。分析以上提取结果可以得出，新方法提取的中轴连通性良好，能较好地保证图形的拓扑特征；该方法具有普适性，不仅对于简单多边形(图 8)能很好地完成中轴提取，而且对于岛洞多边形及含有由自曲线段的图形(图 9)，甚至自然图形(图 10、图 11、图 12)也均能较好地提取出中轴。
|
| 图 7 各类特殊多边形中轴提取结果 Fig. 7 MA extraction result of three curved-polygons |
|
| 图 11 道路网中轴提取结果 Fig. 11 MA extraction result of road |
|
| 图 12 河网中轴提取结果 Fig. 12 MA extraction result of river |
在抑制噪声方面，将传统距离变换方法与新方法分别用于对已加入噪声图形(图 5(b))进行中轴提取，结果分别如图 13(a)、图 13(b)所示。在图 13(a)中，传统距离变换方法中轴提取结果受区域边界噪声影响十分严重，所得中轴含有较多错误；而图 13(b)中显示中轴提取结果良好，基本不受边缘噪声干扰，所得中轴能很好地反映多边形的拓扑和形状特征。
|
| 图 13 新方法和传统距离变换方法中轴提取结果对比(加入边缘噪声) Fig. 13 MA extraction results comparison between new method and traditional method (apply to polygons with noisy points) |
利用传统距离变换方法提取多边形中轴，原理是依据多边形的中轴是到多边形不同边距离相等的点的集合。基本步骤为首先标记多边形的不同边，再进行内距变换并计算出区域内的点到轮廓边界的最小距离值和源点坐标，然后将源点属于同一条边的点赋予相同的属性，将源点属于不同边的点赋予不同属性，最后将多边形区域内不同属性块的边界记为中轴。
传统方法存在以下几点问题。首先，简单地认为中轴是到多边形不同边距离相等的点的集合不够合理，在处理图 9(a)时具有明显的局限性；其次，需要标记多边形的不同边增加了算法复杂度；最后，在有噪声干扰时，采用传统距离变换所获得的距离值不准确，导致所得中轴不准确。而本文提出的新方法采用了更完善的曲边多边形中轴定义，扩大了中轴提取的适用对象范围，通过制定的判别规则自动筛选中轴点，无需确定多边形的边从而降低算法复杂度，同时，采用最近边缘点集距离均值变换有效地抑制噪声的传播和影响，保证所获取中轴的准确性。新方法在适用范围、复杂度和效果上均优于传统方法。
此外，采用新方法所提取的多边形中轴，其含义不仅仅是中轴本身，各个中轴点还包含着3项重要信息：自身坐标值、距其最近边缘点距离值以及最近边缘点坐标值。利用这一特征，能准确、迅速地完成多边形的重建和恢复，这是一般中轴提取方法所难以做到的，同时也为图形数据压/解缩和加/解密提供了一种新的技术途径。
6 结 论研究提出基于最近边缘点集距离均值变换和种子点生长判别法的多边形中轴提取新方法。经过多组不同类型图形的实例验证，该方法不仅适用于简单多边形，还适用于复杂多边形、岛洞多边形、含曲边的多边形，并对自然图形也具有普适性；该方法能够较好地抑制图形边界噪声的干扰，所提取出的中轴具有良好的光滑性、准确性和连通性，且能够很好地反映曲边多边形拓扑和形状特征。在适用范围上和提取效果上，该方法均优于传统的中轴提取方法。此外，由于所提取的中轴包含丰富的信息，该方法还可用于图形的快速重建和恢复以及图形数据的压/解缩和加/解密。
| [1] | LAM L, LEE S W, SUEN C Y. Thinning Metrologies:A Comprehensive Survey[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence,1992,14(9):869-885. |
| [2] | DAKOWICZ M, GOLD C. Extracting Meaningful Slopes from Terrain Contours[C]//ICCS’02 Proceedings of the International Conference on Computational Science-Part Ⅲ. London : Springer-Verlag,2002:144-153. |
| [3] | ZHAO Fulai, HAO Xiangyang. A Research on Contour Plotting Based on Bone Lines[J]. Acta Geodaetica et Cartographica Sinica, 1999,28(4):345-349.(赵夫来,郝向阳.根据地性线绘制等高线的研究[J].测绘学报,1999,28(4):345-349.) |
| [4] | LI Deren, CHEN Xiaoyong. Automatical Generation of Triangulated Irregular Networks for DTM by Mathematical Morphology[J]. Acta Geodaetica et Cartographica Sinica, 1990, 19(3):162-172.(李德仁,陈晓勇.用数学形态学变换自动生成DTM三角形格网的方法[J].测绘学报, 1990, 19(3):162-172.) |
| [5] | BLUM H. New Algorithm for Medial Axis Transform of Plane Domain[J]. Graphical Models and Image Processing, 1967,59(6):463-483. |
| [6] | LEE D T. Medial Axis Transformation of a Planar Shape[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1982, 4(4):363-369. |
| [7] | CHOI H I, CHOI S W, MOON H P, et al. New Algorithm for Medial Axis Transform of Plane Domain[J]. Graphical Models and Image Processing, 1997, 59(6):463-483. |
| [8] | SHAKED D, BRUCKSTEIN A M. The Curve Axis[J]. Computer Vision and Image Understanding, 1996, 63(2):367-379. |
| [9] | CHEN Tao, AI Tinghua. Automatic Extraction of Skeleton and Center of Area Feature[J]. Geomatics and Information Science of Wuhan University, 2004,29(5):443-446.(陈涛,艾廷华.多边形骨架线与形心自动搜寻算法研究[J].武汉大学学报:信息科学版,2004, 29(5):443-446.) |
| [10] | ZHOU Peide, ZHOU Zhongping. An Algorithm for Determining the Medial Axis of an Arbitrary Simple Polygon[J]. Journal of Beijing Institute of Technology, 2000,20(6):708-711.(周培德,周忠平.确定任意多边形中轴的算法[J].北京理工大学学报, 2000, 20(6):708-711.) |
| [11] | ZHOU Peide. Computational Geometry[M]. Beijing: Tsinghua University Press, 2001:212-235.(周培德.计算几何[M].北京:清华大学出版社,2001:212-235.) |
| [12] | LAM L, SUEN C Y. An Evaluation of Parallel Thinning Algorithms for Character Recognition[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1995, 17(9): 914-919. |
| [13] | CHOI W P, LAM K M, SIU W C. Extraction of the Euclidean Skeleton Based on a Connectivity Criterion[J]. Pattern Recognition,2003,36(3):721-729. |
| [14] | CHE Wujun, YANG Xunnian, WANG Guozhao. A Dynamic Approach to Skeletonization[J].Journal of Software, 2003,14(4):818-823.(车武军,杨勋年,汪国昭.动态骨架算法[J].软件学报,2003,14(4):818-823.) |
| [15] | DING Yi, LIU Wenyu, ZHENG Yuhua. Hierarchical Connected Skeletonization Algorithm Based on Distance Transform[J]. Journal Infrared Millimeter and Waves, 2005,24(4):281-285.(丁颐,刘文予,郑宇化.基于距离变换的多尺度连通骨架算法[J].红外与毫米波学报.2005,24(4):281-285.) |
| [16] | IVANOV D, KUZMIN E, BURTSEV S. An Efficient Integer-based Skeletonization Algorithm[J]. Computers and Graphics, 2000, 24(1):41-51. |
| [17] | CHIN F, SNOEYINK J, WANG CA. Finding the Medial Axis of a Simple Polygon in Linear Time[J]. Discrete and Computational Geometry, 1999, 21(3):405-420. |
| [18] | SHAO Chunli. Problems of Medial Axis of a Simple Polygon in GIS and Its Algorithm[D]. Wuhan: Wuhan University,2004.(邵春丽. GIS中多边形中轴问题和算法研究[D].武汉:武汉大学,2004.) |
| [19] | DEGEN W L F. Exploiting Curvatures to Compute the Medial Axis for Domains with Smooth Boundary[J]. Computer Aided Geometric Design, 2004, 21: 641-660. |
| [20] | HU Peng, YOU Lian, YANG Chuangyong, et al. Map Algebra[M]. Wuhan: Wuhan University Press, 2002.(胡鹏,游涟,杨传勇,等.地图代数[M].武汉:武汉大学出版社, 2002.) |