基于角平分线的无孔复杂多边形中轴提取算法

黄斌 刘军 范启东 唐嘉胤

摘要：中轴可以有效表达多边形拓扑结构信息，同時可以去除冗余信息，并且具备对局部几何变形不敏感等特点，吸引了学术界的广泛重视。但现有的中轴提取方法存在计算资源消耗大，运行速度慢，数据不稳定等问题。本文在现有草火模型及最大圆盘模型的基础上，提出了一种基于角平分线的无孔复杂多边形中轴提取算法。该算法利用多边形凹凸点的偏置处理来获得中间分支点，然后通过“手挤”方式获得中轴，并利用边界的连续性保持中轴的连续性。实验结果表明该算法能满足中轴的提取并能保持中轴的连续性，而且减少了自交检测的复杂程度。

Abstract The central axis has the characteristics of effectively expressing the topological structure and geometric information of polygons， removing redundant information， and not sensitive to local deformation， which has attracted extensive attention of academia. There are many problems in the existing methods， such as large consumption of computing resources， slow running speed， and data instability and so on. In this paper， based on the existing grass fire model and the largest disk model， an algorithm of extracting the central axis of complex polygon without hole based on angle bisector is proposed. The algorithm uses the offset processing of polygon concave and convex points to obtain the middle branch points， then uses the way of "hand extruding" the boundary to obtain the central axis， and uses the continuous topological structure of the boundary to maintain the continuity of the central axis. The experimental results show that the algorithm can meet the requirements of axis extraction and maintain the continuity of the axis， and reduce the complexity of self-intersection detection.

1. 引言

随着CAD/CAM/CAE的应用日趋广泛及计算机图形学的不断发展，利用多边形对数字模型进行表达的方式在实际应用中越来越多。但该表达方式会隐藏模型的关键信息，不利于拓展模型的应用。另外，数字模型通常蕴含巨大的数据量，对模型上所有数据进行操作，如重建运动场景，分析结构受力，规划加工路径等，对于大多数硬件设备都是一个巨大的挑战[1-3]。对于数量级巨大的多边形模型，如果要以常规硬件实现大数据处理，就需要将模型的结构信息浓缩，并从中提取核心数据，然后以现有硬件进行有限数量级别操作，最后将操作数据逆向重建，实现数量巨大的多边形模型处理。

结构信息浓缩概念在日常生活中并不少见，例如人类肢体及机械结构的运动简图。人类运动可以看成是忽略神经系统和肌肉细胞控制作用的骨骼运动;而机械结构也可以将复杂的零件特征简化为抽象线框。因此，结构信息浓缩的概念也可以用于将多边形模型浓缩为抽象的骨架线框。从多边形中抽取出的骨架线框，学术界将其定义为多边形中轴。多边形中轴由于包含着丰富的几何机构信息，可以直接应用于科学工程领域，包括地理信息系统、人脸识别、图像处理、计算机视觉和机械结构分析等[4]。

多边形中轴的提取模型主要有两种：草火模型和最大圆盘模型。草火模型是假设图形的边界同时点火，火源向图形内部各个方向等速燃烧直至熄灭，熄灭点形成的点集为多边形的中轴[5];而最大圆盘模型是指平面几何图形中的中轴是内切于几何图形边界至少两个或两个以上等距离的点的轨迹[4]。

多边形中轴具备有效表达多边形的拓扑结构和几何信息，同时可以去除冗余信息，并且对局部变形不敏感等特点，吸引了学术界的广泛重视。而中轴提取更是众多研究学者关注的焦点，合理快速的中轴提取算法一直是学术界及工业界追求的目标。本文将在草火模型和最大圆盘模型的基础上提出一种基于角平分线的改良中轴提取算法。

2. 相关工作

目前对多边形中轴的提取方法都是基于草火模型和最大圆盘模型，常用的中轴提取方法包括有形态学细化法[6]、拓扑逻辑瘦身法（topological thinning）[7]、距离变换法（distance transform/DT）[8]、草火法（simulations of grassfire propagation）[5]、泰森多边形法（Voronoi diagram）等[9]。而从计算精度来看，众多的中轴提取方法大致可以分为精确算法和逼近算法两类[12]。

精确中轴提取算法包括形态学细化法、拓扑逻辑瘦身法、距离变换法等。形体学细分法属于图像处理范畴，通过膨胀、腐蚀、开闭运算、击中或击不中变换、细化、骨架及重构对象等一系列运算获得最终的中轴图形[1]。形态细化法最常用的就是最大圆盘法[6]，通过多边形内所有最大圆盘的圆心提取多边形的中轴，但存在着连通性不佳的缺点。拓扑逻辑瘦身法[7]是在形态学细化法的基础上发展出来的中轴提取方法，通过逐步去掉简单点（不影响拓扑性质和形状的点）而获得中轴，但其仅能保证在二维空间的正确性。距离变换法是通过寻找距离梯度脊线来形成骨架，该方法得到的骨架位置精确，但难以保证骨架的连通性。但精确中轴提取算法多采用全局遍历法，存在计算资源消耗大，运行速度慢等问题;另外一些精确中轴提取算法将图形进行分割成简单多边形，然后提取中轴，最后将各部分中轴融合成为最终的总中轴，则在中轴融合时出现准确性不佳及计算复杂等问题。

针对精确中轴提取算法存在的问题，研究人员相继开发出效率更高的逼近算法。逼近算法主要包括扩展草火法[13]、泰森多边形法[9]、Delaunay三角剖分法[10]和栅格法等。扩展草火法采用类似草火法的统一方法来定义二维中轴全局形状度量（EDF）及二维形状的中心（EMA），并通过厚度腐蚀方法得到图形的中轴。该方法在修剪中轴、对齐形状和形状描述方面的具有一定的实用性。泰森多边形法[9]将复杂多边形界线分解成点集，然后通过这些点集构成泰森多边形图（Voronoi diagram），泰森多边形图在多边形内的公共边即是多边形的中轴，但该方法依赖于规模或密度。Delaunay三角剖分法[10]利用Delaunay三角剖分的空球特性（平面上是圆），对边界点集细化，以提取多边形的中轴。栅格法[11]则通过最近边缘点集距离的均值变换，结合新的中轴点判定规则，利用种子点生长判别法提取曲边多边形的中轴。但基本上使用逼近中轴提取算法均存在数据提取不稳定的缺点，迄今为止还没得到很好的解决。

3. 理论算法

针对上述精确及逼近中轴提取算法的不足，本文利用结合角平分线的概念及多边形边界拓扑连续性，实现无孔复杂多边形的中轴提取。无孔復杂多边形的中轴提取的理论主要分为以下四个基本算法：多边形凹凸判断算法，凸多边形偏置算法，非凸多边形 偏置算法，及中轴重构算法。

3.1 多边形凹凸判断算法

要从复杂的封闭多边形进行中轴提取，首先要定义一个实体方向的概念。

多边形实体方向定义（如图1）：

多边形的边界线段是有向线段或多段线。

有向线段或多段线的前进方向的左侧为实体，则该线段或多段线的前进方向为多边形的实体方向。

要确定多边形的凹凸性。在各种多边形中轴定义模型中，多边形的轮廓实体方向的凸点存在中轴，凹点不存在中轴;但凹凸点也会影响中轴的位置，因此需要对多边形的凹凸性进行判断。

多边形顶点凹凸性判断方法主要有角度法、左右点法、矢量面积法、向量积法、射线法和极点顺序法等[14]，这些算法基本上是等效的。本文采用的是向量积法来判断多边形顶点的凹凸性。如图1所示，边界ABC（C）的左侧为实体，顶点C位于矢量线段AB的左侧，则计算向量AB和向量BC的外积，外积公式采用：

向量BD的单位向量与z轴的正方向相同，这时多边形的顶点属于凸顶点。而对于顶点C而言，其位置位于矢量线段AB的右侧，向量BD的单位向量与z轴的正方向相反，这时多边形的顶点属于凹顶点。

由于多边形的轮廓上有多个顶点，因此必须采用历遍法对各个顶点进行凹凸性判断，并最终确定多边形的凹凸性：若全部顶点均为凸顶点，那么，该多边形为凸多边形;反之，若多边形上存在一个或以上的凹顶点，该多边形为非凸多边形。

3.2凸多边形偏置算法

草火模型和最大圆盘模型均采用其描述方式对多边形的中轴进行定：草火模型采用的是由外至内的距离延伸方法，而最大圆盘模型采用的是相反的由内至外的距离延伸方法。通过比较两种定义方法，可以发现其最终的目标都是在寻找距离相等的点或线段。而在平面几何学中，角平分线的性质就刚好满足这一要求，如图2。林福严等[15]利用这角平分线的性质提出了一种平面凸多边形中轴提取方法，但该方法的判据并未考虑多边形的特殊情况。而这里提出的凸多边形偏置算法是在林福严等提出的方法上进行改良，将多边形的特殊情况也进行考虑。

凸多变形偏置算法描述如下：

步骤1：寻找凸多边形的最小厚度。

（1） 按凸多边形实体方向获取各顶点的坐标，然后计算每一个顶点的角平分线的单位向量，求出两顶点间线段与两顶点法向量所构成的三角形的高（图3）。三角形的高的计算公式如下：

（2） 比较凸多边形上所有轮廓线段对应三角形的高，寻找最小值，作为多边形的最小厚度。凸多边形顶点角平分线与顶点简单的多边形轮廓线段构成了三角形，三角形对应高的最小值选择方式如下：

（3） 求出最小厚度所在的三角形上由两多边形顶点角平分线构成的顶点坐标，作为中轴线的中间分支点。计算方式如下（以图2为例，A为第一个中间分支点）：

步骤2：寻找第一个无自交的凸多边形环。

（1） 按凸多边形实体方向对每个顶点角平分线进行截点计算。计算方式如下：

（2） 按凸多边形顶点的拓扑关系对截点进行连接，构成了一个新的凸多边形（图4）。由于最小厚度所在三角形对应的凸多边形轮廓线段已退化成一点，如图7的V1、V2两点退化成P1点，新凸多边形的边数为N-1。

图4利用截点的拓扑关系形成新的凸多边形。

步骤3：找出所有中间分支点。 重复步骤1和2，找出所有无自交的凸多边形环及所有的中间分支点，直到多边形上所有顶点退化到一点。这时，中间分支点到其对应的多边形轮廓线段的距离为多边形在该轮廓线段上的最大厚度。

3.3非凸多边形偏置算法

凸多边形仅是多边形中的特例，而实际上更多存在的是非凸多边形。非凸多边形的偏置算法长期以来都是以一个难以解决的问题，而问题关键是在于如何处理非凸多边形轮廓在偏置后出现的自交现象。自交现象产生的原因是由于偏置距离设置与轮廓线拓扑结构的相互影响而产生的，如图5。正常情况下，轮廓ABCD在偏置d1距离后会生成偏置曲线ABCD，如图5（a），但若偏置距离过大时，轮廓ABCD的偏置曲线ABCD在保持原有拓扑结构的情况下会出现AB和CD两段曲线自交，如图5（b）这就是自交产生的原因。自交现象的出现会增加中间分支点位置的不确定性。对于这种情况，通常的处理情况是在偏置后对所有线段进行自交检测，但会随着多边形变数的增多，算法的执行效率会降低。针对这种情况，在凸多边形偏置算法的基础上，本文通过整合临界偏置处理，如图5（c），提出了一种基于非凸多边形最小厚度的角平分线偏置算法，其核心在于找到偏置处理的临界点。

非凸多变形偏置算法描述如下：

步骤1：寻找非凸多边形的最小厚度。

（1） 按凸多边形偏置算法的步骤1计算出两顶点间线段与两顶点法向量所构成的三角形的高。但不同的地方在非凸多边形上计算出来三角形的高是带正负方向的，因此，这里规定只有正向的三角形高才是有效数值。

（2） 多边形凹凸判断算法找出非凸多边形轮廓线上的凹顶点，并计算凹顶点的角平分线正方向（从顶点交出发指向多边形内部方向）与距离最近的轮廓边界线段的交点。具体方法如下：

寻找跟凹顶点的角平分线正方向相交的轮廓边界线段。假设角平分线与轮廓边界线段相交，轮廓边界线为直线线段，轮廓线段的两端点必在角平分线的两侧。图6所示为角平分线与轮廓边界线段相交的判断，BQ为多边形凹顶点角CBA的正向角平分线，与轮廓边界线段EF交于P点，端点E和F分别位于BQ左侧和右侧。然后引入相交判断系数，并进行计算：

为E到BQ的有向距离，

为F到BQ的有向距离，

为的单位向量

为的单位向量

为相交判断系数。当等于0时，角平分线与轮廓边界线段相交。

求出交点P的位置。借助有向距离的模及边界线段端点坐标，可以快速计算出交点坐标。计算方式如下：

为角平分线与轮廓线段的交点坐标

为角平分线与轮廓线段的交点坐标

计算所有与同一凹顶点的角平分线相交边界线段的距离，找出最短距离，确定最近的边界线段。

计算凹顶点所在角平分线与相交边界线段等高的交点，如图7。计算方法如下：

比较所有凹顶点角平分线及最近轮廓边界线的等高距离，寻找最小值，

和分别为第个凹顶点角平分线及最近轮廓边界线的等高距离及最后一个凹顶点角平分线及最近轮廓边界线的等高距离。

（3） 比较凹顶点等高距离及轮廓边界线对应三角形的高，寻找最小值作为非凸多边形的最小厚度.

步骤2：寻找第一个无自交的非凸多边形环。

（1） 利用式（4）计算非凸多边形实体方向的每个顶点角平分线上的截点计算。

（2） 按凸多边形顶点的拓扑关系对截点进行连接，构成了新的多边形。非凸多边形在这里会出现两种情况，第一种情况，轮廓边界线对应三角形的高为最小值，最小厚度所在三角形对应的凸多边形轮廓线段会退化成一点，新的多边形只有一个且边数为N-1，与凸多边形偏置的情况一样。

第二种情况，凹顶点等高距离为最小值，凹顶点角平分线正方向上与最近轮廓边界线段的偏置线段会出现分支交点（图8），此时，新的多边形会出现两个或以上的多边形。由于分支交点P5既属于角平分线V5P5上的点，又属于轮廓线段V1V7的偏置线段P1P7上的点，但实际上偏置线段P1P7上并无P5这点，需要通过插值的方式将P5放到P1P7上以保持多边形边界的合法性。分支交点也是非凸多边形中的中间分支点。

步骤3：找出所有中间分支点。对于分割出来的新多边形，则对其的多边形凹凸性进行判断，然后按照其属于凸多边形或非凸多边形的偏置算法进行处理。如果新生成的多边形是凸多边形，则按凸多边形偏置方法找出所有无自交的凸多边形环以及所有的中间分支点，直到多边形上所有顶点退化到一點，如果新生成的多边形是非凸多边形，则按步骤1和2 进行操作，找出所有无自交的非凸多边形环以及所有的中间分支点。偏置操作直至该非凸多边形中内再无多边形环可提取。

3.4中轴重构算法

中轴重构算法可分为对两类：凸多边形中轴重构算法及非凸多边形中轴重构算法。而重构算法的关键在于保持中轴拓扑结构的连续性。本文提出一种“手挤（Hand Squeeze/HS）”式中轴重构算法，利用轮廓边界拓扑结构的连续性来保持凸多边形及非凸多边形中轴的连续性。

3.4.1 凸多边形中轴重构算法

凸多边形中轴重构算法如图9， 假设用手握着凸多边形，然后均匀用力握住凸多边形。由于手形的生理特征，多边形较长的方向会与四指弯卷方向呈锐角方向。另外假设凸多边形的硬度为零，轮廓边线可以随意向垂直方向挤压，并在挤压过程中，轮廓边线允许轴向收缩。基于这两点假设，凸多边形可以向内部收缩直到轮廓边线间没有收缩“空间”，这样就成为该凸多边形的中轴。但此时，轮廓边线间还是连续的。通过这种算法，可以快速得到连续中轴。而这种算法的定义为：

定义1：在顶点偏置过程中保持每个顶点间的拓扑关系，不作任何添加和修改。

定义2：在到达最小偏置距离及出现中间分支点时，允许出现重合点，但仍旧保持轮廓的拓扑连续关系。

定义3：从中间分支点继续偏置时，假设在“手挤”过程中，物体已经到达“背靠背”的状态其无法继续挤压，从而形成一条直线段。此时，分支点不再移动，其他点仍然向线段的垂直方向“挤压”。这种情况下，需要在重合点间插入中间分支点并固定，其它点仍沿着顶点的角平分线正方向移动，如图10。

定义4：继续挤压直到所有轮廓边界线段间没有“空间”，得到中轴，如图11。

3.4.2 非凸多边形中轴重构算法

非凸多边形中轴重构算法如图12，与凸多边形的“手挤”中轴重构算法类似。但不同的是凸多边形是单手“手挤”压缩，而非凸多边形是多手在两凹点间构成的半凸多边形进行“手挤”压缩。非凸多边形的“手挤”压缩过程的假设跟凸多边形一样，通过均匀施加压力，挤去非凸多边形上的冗余空间，得到非凸多边形的中轴。非凸多边形的“手挤”压缩中轴提取的定义如下：

定义1：在顶点偏置过程中保持每个顶点间的拓扑关系，不作任何添加和修改，同凸多边形中轴重构算法定义1。

定义2：在到达最小偏置距离及出现中间分支点时，允许出现重合点，但仍旧保持轮廓的拓扑连续关系，同凸多边形中轴重构算法定义2。

定义3：从中间分支点继续偏置时，到达“背靠背”的状态，且形成一条直线段。此时，分支点不再移动，其他点仍然向线段的垂直方向“挤压”。这种情况下，需要分两种情况。第一种情况，如果是凸顶点形成的中间分支点，处理方式同凸多边形中轴重构算法步骤3，即在重合点间插入中间分支点并固定，其它点仍沿着顶点的角平分线移动。第二种情况，如果是凹顶点形成的中间分支点，除了将中间分支点固定以外，还需要插入额外控制点，并沿着先出现的角平分线正方向移动（图13）。

定义4：继续挤压直到所有轮廓边界线段间没有“空间”，得到中轴，如凸多边形中轴重构算法定义4。

4. 实验验证

对于上述理论，需要通过实验进行验证，以下将分别对凸多边形和非凸多边形进行实验。

4.1 凸多边形

为确保凸多边形偏置算法的有效性，实验样本为随机一个凸多边形。在凸多边形的实验样本上，按逆时针顺序作各顶点的角平分线，计算出第一个无自交的凸多边形环的最小偏置距离，利用式（4）求出的中间分支点及第一个无自交凸多边形环上的控制点，连接成新的凸多边形（图14）。然后按3.2凸多边形偏置算法的步骤3按计算得出的各最小偏置距离进行偏置，直至所有凸多边形顶点退化到一点。最后中轴提取方法提取凸多边形的中轴（图15）。

4.2 非凸多边形

同样，为确保非凸多边形偏置算法的有效性，实验样本为随机一个非凸多边形。在非凸多边形的实验样本上，按逆时针顺序作各顶点的角平分线，计算出第一个无自交的非凸多边形环的最小偏置距离，利用式（4）求出的中间分支点及第一个无自交非凸多边形环上的控制点，连接成新的非凸多边形（图16）。然后按3.3凸多边形偏置算法的步骤对所有无自交凸多边形环进行偏置操作直至无法再提取多边形环为止。最后中轴提取方法提取凸多边形的中轴（图17）。

5.结果分析

从凸多边形和非凸多边形的偏置处理及中轴提取实验结果来看，基于角平分线的无孔复杂多边形中轴提取算法可以基本满足多边形的偏置及中轴提取的目的，同时可以减少自交检测的复杂度。

6.结论

本文在现有草火模型及最大圆盘模型的基础上，提出了一种基于角平分线的无孔复杂多边形中轴提取算法。该算法利用多边形凹凸点的偏置处理来获得中间分支点，然后利用“手挤”边界的方式获得中轴，并利用边界的连续拓扑结构保持中轴的连续性。实验结果表明该算法能满足中轴提取的要求并能保持中轴的连续性，而且减少了自交检测的复杂程度。但目前该算法局限于无孔的复杂多边形结构，尚未进行有孔多边形验证。下一步将进行有孔多边形验证及编写程序验证。

参考文献

[1] 黎茂. 中轴变换研究[D]. 电子科技大学， 2006

[2] 徐春. 一种基于手的力反馈虚拟现实系统的研究[J]. 毕业生， 2003.

[3] 巴文兰， 曹利新. 基于中轴变换的刀具优化选择与刀具路径规划[J]. 大连理工大学学报， 2013（01）：63-68.

[4] Smogavec G， Zalik B. A fast algorithm for constructing approximate medial axis of polygons， using Steiner points [J]. Advances in Engineering Software， 2012， 52：p.1-9.

[5] Blum，H. A transformation for extracting new descriptors of shape [J]. Models for the Preception of Speech & Visual Form， 1967， 19：362-380.

[6] 吕哲， 王福利， 常玉清， et al. 改进的形态学骨架提取算法 [J]. 计算机工程， 2009， 035（019）：23-25.

[7] Arcelli， Carlo & Sanniti di Baja， Gabriella. （1985）. A Width- Independent Fast Thinning Algorithm. Pattern Analysis and Machine Intelligence， IEEE Transactions on.4. 463-474.

[8] 徐超， 肖潇， 骆燕， 等. 基于距离变换的新型骨架提取方法 [J]. 仪器仪表学报， 2012（12）：213-218.

[9] Dey T K， Zhao W. Approximate medial axis as a Voronoi subcomplex[J]. Computer Aided Design， 2004， 36（2）：195-202.

[10] Yong， L.K. （2009）. 3D medial axis approximation through Delaunay triangulation. Journal of Science and Technology in the Tropics. 5. 133-139.

[11] 潘鵬， 贺三维， 吴艳兰， 等. 曲边多边形中轴提取的新方法 [J]. 测绘学报， 2012， 041（002）：278-283，290.

[12] 王新生， 谢凯， 姜友华， 等. 复杂多边形中轴构建方法 [J]. 武汉大学学报·信息科学版，2014，39：2， 2014， 39（2）：181-185.

[13] Liu L ， Chambers E W ， Letscher D ， et al. Extended grassfire transform on medial axes of 2D shapes[J]. Computer Aided Design， 2011， 43（11）：1496-1505.

[14] SONG Xiaomei， CHENG Changxiu， ZHOU Chenghu. An Analysis and Investigation of Algorithms for Identifying Convexity-Concavity of a Simple Polygon [J]. Remote Sensing for Land & Resources， 2011， 19（3）：25-31.

[15] 林福严， 郑奇志， 胡于进， 周济. 平面形体的中线提取及其应用[J].计算机应用与软件 2000