基于距离变换的多尺度连通骨架算法
丁颐 刘文予
（华中科技大学 电子与信息工程系 宽带无线与多媒体系统研究中心 武汉 430074）
摘要：传统的基于距离变换的骨架算法不能保证骨架的连通性，需要引入鞍点解决连通问题。
这类算法复杂，且不够准确，会引入伪骨架点，同时鞍点的定义很难推广到三维，限制了传
统算法的发展。本文提出一种新型骨架算法，在图形内根据距离变换的约束，由骨架种子点
开始以单像素宽度逐点生长出各骨架分支，逐点生长保证了连通性。该算法的骨架生长过程
是骨架由粗到精的演变过程，能够方便地实现骨架的多尺度控制。
关键词：距离变换、骨架、多尺度、最大圆
A Hierarchical Connected Skeletonization
Algorithm Based on the Distance
Transform
Ding Yi, Liu Wenyu
Electronics & Information Engineering Department, Huazhong University of Sci. & Tech.
Wuhan, Hubei, 430074 P. R. China.
Abstract: The connectivity property is not guaranteed by the traditional skeletonization algorithm
based on the distance transform, saddle points should be added to solve the connectivity problem.
These methods are complex, inaccurate(pseudo points may be inducted), and not adapted to 3D
case. A novel method has been presented in this paper, the skeleton is obtained by growing from
the skeleton seed with 1 pixel width restricted by the distance transform. The connectivity is
easily assured, and the hierarchical control can be achieved from the growing process.
Key word: Distance transform, skeleton, hierarchical, maximal disks
1 引言
1967 年 Blum 通过烧草模型(Grassfire Model)首先给出了骨架的定义[1]，他假设图形边
界点同时着火，火源向图形内部各个方向等速的燃烧直至熄灭，所有熄灭点的点集构成了
该图形的骨架。Blum 还给出了骨架的另一种等效定义[2]，即骨架是所有最大圆 (Maximal
Disks)的圆心集合，最大圆是完全包含在图形内部的圆，并且不被任何其它包含在图形内部
的圆所包含。
骨架是原始图形的一种压缩表示，与原始图形保持了相同的拓扑结构，并且存在于图
基金项目：国家自然科学基金资助项目()
形的对称轴上，能够同时反映图形的拓扑与形状信息。因此骨架在变形、识别、跟踪、重
建、压缩、路径导航、检索等领域都具有广泛的应用。对于骨架的求法，当前的研究主要
可以分为以下两类，这两类研究是根据前面的两种定义发展而来的：
（1）细化的方法。这种方法就是在模拟烧草的模型，即逐层均匀的剥掉图形的边界，
最后剩下最里层的已经不能再剥掉（否则会影响连通性）的部分就构成了图形的骨架[3,4,5]。
（2）基于距离变换的方法。根据最大圆的定义可以推导出最大圆必是图形的内切圆，
因此考察以某点为圆心作的内切圆，如果该圆不被其它的任何内切圆包含，则此点作的内
切圆是最大圆，此点是骨架点。距离变换记录了图形内所有点到边界的最短距离，利用距
离变换可以生成任意点的内切圆及比较两个内切圆的包含关系[6,7,8,9]。
两种方法在连续域都是完备的，实际应用到离散域中，各有优缺点。细化得到的骨架
在拓扑保持上有较好的表现，但是骨架的位置却不准确。其原因在于离散域里任何点的行
进方向最多只能有八个（二维），远不能模拟火烧过程中任意的行进方向，因而产生误差。
文献[5]引入 snake 模型调整骨架的位置，但是提高了复杂度。基于距离变换的方法在骨架
点的准确度上有明显的优势，但连通性却很难保证。其原因在于离散域内最大圆定义和圆
的包含关系很难把握，传统的最大圆判断准则得到的骨架不保证连通性[6,7]，需要另外引入
一些鞍点(Saddle Points)将不连通的部分连通起来。近年来出现了一些改进的距离变换骨架
算法[8,9]，但是这些方法主要针对于柱状物体，例如医学应用的肠道导航。
针对这两类算法的局限性，本文提出一种新型骨架算法，该算法利用距离变换信息从
骨架种子点以单像素宽度生长出其余的骨架点，该过程保证了骨架的连通性及单像素性，
并且可以根据给定权值控制生长精度，实现骨架的多尺度控制。
2 骨架算法
传统方法的缺点
前面已介绍，传统的基于距离变换的方法考察图形内所有点的内切圆，如果一个内切
圆不被任何其它的内切圆包含，则此内切圆是最大圆。但是如果每个内切圆都要与图形内
所有点的内切圆遍历比较，复杂度太高。传统的方法大都是比较每个内切圆与其邻域点的
内切圆[6,7]，这样做得到的点集是骨架点集的父集，因为某点内切圆不被邻域的内切圆包含
却仍然有可能被图形内其它内切圆包含，这部分伪骨架点应该被去掉，文献[7]对此做了改
进。但是即使上面的问题解决了，骨架点集找全了，结果仍然是不理想的，骨架不能保证
连通性。文献[7]分析了失去连通性的原因，认为该方法在连续域是完备的，但是对应到离
散域，圆的定义（如图 1）造成两圆包含关系的连续域判断准则不能直接挪用，这是失去连
通性的根本原因。图 1中(a)、(b)、(c)、(d)分别代表半径平方值是 1、2、4、5的离散圆，(e)
表示的是一个 H 形的图形及它的骨架，(e)中标出了各像素点的距离变换平方值，图中深色
点是骨架点，因为以它们为圆心的内切圆都不能被其它内切圆完全包含，由图可见该骨架
是不连通的，且不是单像素宽的。
中国科技论文在线
图 1 传统的基于距离变换的骨架
为了避免传统方法的上述缺陷，本文提出一种由骨架种子点向外生长的方法以保证骨
架的连通性和单像素性，并且在生长的过程中利用距离变换信息，保证骨架的准确度。
该算法的步骤如下：
(1) 计算距离变换。
(2) 寻找骨架种子点。该步的目标是找到一个合适的起始骨架点；
(3) 生长骨架。这一步考虑如何从一个已知的种子点经过反复迭代的过程找出其它所
有的骨架点。
下面将对这三个步骤做具体的讨论。
计算距离变换
每个像素或体素的距离变换值是该点到图形边界的最短距离，任何基于距离变换的骨
架算法都要首先计算距离变换。计算距离变换通常有两种方法，一种是基于模板的近似方
法[10]，另一种是精确方法[11]。前者复杂度低，效率高，但是计算的结果有误差，后者计算
结果精确，只是复杂度略高。本文选用文献[11]提出的计算结果精确的距离变换计算方法。
寻找骨架种子点
本文算法能够保证从任何一个骨架点出发生长出全部骨架， 节将给出具体论述，
因此理论上骨架种子点的选择十分自由，任何一个骨架点都可以作为种子点。但是选择不
同的种子点构造多尺度骨架的效果不同，理想效果是由骨架中间部位发散地向外生长，这
样在开始获得较少骨架点的情况下，也能最大程度地表现图形的全局特征。骨架中间部位
点的特征是距离变换值相对较大，因此本文选取图形内距离变换值最大的点作为种子点，
可以证明该点的内切圆一定不会被其它内切圆所包含，所以其必定是骨架点，符合种子点
的基本要求，并且位于骨架的中间部位能够生长出效果满意的多尺度骨架。图２所示为同
一图形选择不同种子点生长骨架的过程，上一行的种子点位于图形的边缘，下一行的种子
点选择的是距离变换值最大的点，可以看出虽然骨架最终结果都很理想，但是下一行更能
反映骨架由粗取精的过程，多尺度性表现得更好。
中国科技论文在线
图２ 种子点的选择
生长骨架
生长骨架的思路来源于骨架的一个应用——图形重建。利用骨架重建图形，就是以每
个骨架点为圆心，该点的距离变换值为半径画圆，所有这些圆总共覆盖的区域就是原始图
形。将这个过程反过来考虑，在骨架点还没有搜索完全时，新的骨架点应该如何生长，即
考虑引进哪些骨架点能够更大程度的覆盖原始图形。生长骨架的过程是一个迭代过程，本
文定义上一轮迭代新产生的骨架点称为下一轮迭代的生长前沿点 f，每轮迭代的任务就是在
各生长前沿点的邻域生长出一个或多个本轮的新骨架点。每轮迭代过程分以下三步进行：
（1）覆盖。覆盖这一步的任务是以生长前沿点为圆心，相应的距离变换值为半径，覆
盖图形。被覆盖的图形部分将被挖去，每轮迭代的操作都会减小图形，直至图形的剩余部
分为零，即找到的骨架点已足以覆盖整个图形，整个迭代过程结束。图３显示了上述的操
作过程，其中（a）为原始图形，（b）为找到骨架种子点后的第一轮迭代，图中表示了由骨
架种子点画圆覆盖后的剩余图形，（c）为第二轮迭代的覆盖结果，此时有３个生长前沿点，
（d）为第六轮迭代的覆盖结果，（e）为第十一轮的覆盖结果，此时有４个生长前沿点，（f）
为第十四轮的覆盖结果，此时剩余图形为零，迭代过程结束，此图表示的就是骨架的最终
结果。
(a) (b) （c） （d） （e） （f）
图３ 覆盖的过程
（2）判断新分支数。每个生长前沿点都对应图形剩余部分的一个连通区域，经过覆盖
后，该连通区域包含的点数将会减少，并且被拆分成了 n个连通区，n为大于等于 1的整数，
这 n 个连通区都是此前沿点应该生长的方向，因此需要产生 n 个分支，继续分别覆盖这 n
个连通区。如图３所示，最初种子点对应的是整个图形的一个完整连通区，当第一轮执行
覆盖后，如（b）所示，连通区分为了三块，于是由（b）中的 f1应该生长出三个新的分支，
如（c）所示，第二轮迭代时就已经有三个生长前沿点代表的三个分支了，分别对应了图形
剩余的三个连通区域。
判断编号为 k 的连通区经覆盖一次后，被拆分成多少个连通区可采用类似于图的深度
遍历的算法，伪码如下所示。其中 image[][]表示剩余图形，若点（i,j）属于连通区 k，则
image[i][j]=k。found[][]对已经搜索到的该编号连通区的点进行标记，若点（i,j）属于连通区
中国科技论文在线
k且已被搜索到，则 found[i][j]=true。search(x,y)函数从点(x,y)出发深度遍历所有与其相连通
的点集。over表示是否已搜索完属于原连通区 k的点。n表示被拆分成的连通区个数，n的
合法输出应是大于等于 1的整数。
int splitnumber()
{
n=0;
over=false;
while(over!=true)
{
over=true;
for(j=0;j<imageheight;j++)
for(i=0;i<imagewidth;i++)
{
if((image[i][j]==k)&&(found[i][j]==false))
{
n++;
over=false;
search(i,j);
continue;
}
}
}
return n;
}
search(int x,int y)
{
found(x,y)=true;
if((image[x-1][y-1]==k)&&(found[x-1][y-1]==false))
search(x-1,y-1);
if((image[x][y-1]==k)&&(found[x][y-1]==k))
search(x,y-1);
if((image[x+1][y-1]==k)&&(found[x+1][y-1]==false))
search(x+1,y-1);
if((image[x-1][y]==k)&&(found[x-1][y]==false))
search(x-1,y);
if((image[x+1][y]==k)&&(found[x+1][y]==false))
search(x+1,y);
if((image[x-1][y+1]==k)&&(found[x-1][y+1]==false))
search(x-1,y+1);
if((image[x][y+1]==k)&&(found[x][y+1]==false))
search(x,y+1);
if((image[x+1][y+1]==k)&&(found[x+1][y+1]==false))
search(x+1,y+1);
}
中国科技论文在线
（3）生长新骨架点。上一步判断了分支数，这一步考虑如何生长分支，即选择前沿点
的哪些邻域点作为本轮新骨架点 f’i(i=1~n)，之后开始新一轮的迭代。新骨架点 f’要求比前
沿点 f距离对应的连通区更近，以便下一轮更多地覆盖连通区，在此基础上 f’的距离变换值
要求在局部上有较大值，以保证骨架位于图形的中央。本算法首先找出 f的 8个邻域点里距
离连通区最近的一个作为 f’的初选结果，然后考察 f的 8邻域里与 f’相邻的两个点，这两个
点加上 f ’共覆盖了 f的 90o的领域范围，由于已经确定了大致的方向，精确的新骨架点位置
不会跳出这个范围，因此计算 f’及相邻两点相对于 f的距离变换值的下降程度，其中最缓慢
的一个显然是距离变换值局部较大值，选这个点作为 f’的最终结果。可以证明初选结果保证
了新骨架点朝着未被覆盖的连通区生长，即明确了大体方向，而最终结果在局部上保证了
所选择的新骨架点是位于图形中央的，满足骨架的中轴性质。如图４所示，表示的是图３
（d）第六轮迭代时右上方的分支的生长情况，此时的前沿点 A坐标为（19，9），此点的距
离变换值是２０，其８邻域的距离变换值在图中标示出来。此点右上方的一块连通的深色
区就是对应此点的剩余图形连通区，这部分剩余图形的平均位置在 B（22，5），在 A 的８
个邻域点中，C是距离 B最近的点，因此 C将作为新骨架点的初选点。接下来考虑 A的８
个邻域中与 C 相邻的两个点 D 和 E，分别计算 C、D、E 三点相对于 A 的距离变换下降梯
度， ( )12020 −=dg , ( ) 21820 −=cg ， ( )11320 −=eg ,其中 gd值最小，
说明 D点相对于 C和 E更接近于图形的对称轴，我们由此选择 D点作为本轮的新骨架点。
下面给出生长新骨架点的伪码描述，其中 uncover表示连通区的平均位置，current表示
当前前沿点，current[]表示 current的 8个邻域点，编号 0～7分别表示左上、上、右上、右、
右下、下、左下和左 8个方位，new表示新骨架点。函数 distance计算两点距离，DT计算
该点的距离变换值。
SelectNewNode()
{
n=0;
dmin=distance(current[0],uncover);
for(i=1;i<8;i++)
{
d[i]=distance(current[i], uncover);
if(d[i]<dmin)
{
图４ 生长新的骨架点
中国科技论文在线
dmin=d[i];
n=i;
}
}//n为距离连通区距离最近的 8邻域点的编号
n[1]=(n+7)%8;
n[2]=n;
n[3]=(n+1)%8;//n[1],n[2],n[3]分别为编号 n及它在 8邻域上的两相邻点的编号
m=1;
decentmin=(DT(current)-DT(current[n[1]]))/distance(current,current[n[1]]);
for(i=2;i<=3;i++)
{
decentmin[i] =(DT(current)-DT(current[n[i]]))/distance(current,current[n[i]]);
if(decentmin[i]<decentmin)
{
decentmin=decentmin[i];
m=i;
}
}
new=current[n[m]];
}
３ 骨架的连通性和单像素性保证
传统算法各骨架点的判断都是独立的，不考虑骨架点之间的位置连通关系，因而不保
证连通性。本文算法由起始种子点经迭代逐点生长骨架，每轮新生长的骨架点都是上一轮
骨架点的邻域，所以连通性一定能够得到保证。
对于单像素性，同样要考虑骨架点之间的关系。本文算法每轮迭代过程中对应每个连
通区生长一个分支点，各分支就是由这样的单点组成的，因此只要不同分支上的点不相邻，
同一分支上的点不跳跃相邻，就能保证单像素性。图 5 说明两种失去单像素性的情况，前
者两分支靠在一起，后者分支自身扭曲，针对这两种情况，本文提出了以下两条限制准则：
（1）每轮迭代中，由同一个前沿点产生的新骨架点不能相邻。如图 5(a)中，f’1 和 f’2
都是由 f 产生的，这种情况下，f’1与 f’2必须被强行间隔开，或者重叠在同一点，(a)的问题
就能解决。
（2）任何新产生的骨架点不能与除产生它的前沿点以外的所有已生成骨架点相邻。如
图 5(b)，s5与 s4、s3、s2都相邻，但只有 s4是 s5的前沿点，因此 s5不能选择图 5(b)中的位置。
(a) (b)
图 5 影响单像素性的两种情况
实验结果表明，骨架的连通性效果理想,同时由于考虑了骨架点之间的位置关系，经过
以上两条准则约束的骨架有效地实现了单像素性。
中国科技论文在线
４ 骨架的多尺度控制和识别信息的提取
本文算法十分有利于多尺度骨架的生成，同时也具有良好的抗边界噪声鲁棒性。
骨架生长过程的结束条件根据是否有剩余的未被覆盖的连通区决定，每个骨架点被选
中时，都会对应一个连通区，而这个连通区范围越大就说明这个骨架点的权重越大，这个
权重可以作为我们实现多尺度控制的标准。我们规定连通区包含的点数必须大于给定权值
k，才可以继续生长骨架，否则该骨架分支结束。权值 k越小，骨架越精细；否则，骨架越
粗糙。为了避免噪声的影响，选取一个合适的 k 值就能达到理想的效果了。图 6 中的图形
大小是 100×100，从左到右 k的取值是 0，10，20，30和 40。
k=0 k=10 k=20 k=30 k=40
图 6 枫叶的多尺度骨架
骨架的多尺度控制可以为图形识别所用，实现图形由粗到精不同程度的匹配及分类。同
时骨架的单像素及连通性，也有利于匹配算法的操作，除此之外，本算法还可提供下列附加
的有利于识别的图形拓扑与形状信息。
（1）骨架点的距离变换值。图形内所有点的距离变换值都预先计算，包括最后生成的
各骨架点。距离变换信息涵盖了图形的形状信息，利用此信息可以实现图形的重建，在图形
匹配过程中可以作为形状匹配的参考。
（2）骨架点的分支数。骨架生长过程中计算了每点的分支数，分支数可以用来判断该
骨架点是分叉点，普通点还是终结点，便于物体的描述与识别，在后续的研究中，我们利用
骨架点的分支数，可以很方便的实现基于多尺度骨架的物体拓扑分类。
5 实现结果及结论
图 7 给出本文方法得到的部分实验结果及与文献[7]方法（未经去除伪点及连接）的比
较。上一行是本文方法的结果，下一行是文献[7]的结果。可以看出本文的方法十分有效的
解决了骨架的单像素及连通问题。
图 7 实验结果
本文算法需要改进之处在于目前还不能适用于带孔的图形，因为带孔图形骨架的分支
中国科技论文在线
数与剩余的连通区域数是不等的，有可能几个分支对应一个剩余连通区，需要对算法作相
应的改进。
本文算法改进了传统的基于距离变换的方法，解决了骨架连通性的问题，且保证了单
像素宽度，实现了骨架的多尺度变换，增加了骨架的附加信息，这些附加信息对于后面的
物体描述和分类非常有用。相对于细化的方法，本文的结果更准确。在算法复杂度上较传
统的距离变换法稍有提高[6,7]，若 n是二维图形边长，传统方法的复杂度是 O(n2)，本文的方
法是 O(n3)，与细化相当。
本方法可以推广到三维，每轮迭代用球代替圆来覆盖图形，得到三维图形的线型而非
面型三维多尺度骨架。
References:
[1] H. Blum. A Transformation for Extracting New Descriptors of Shape. Models for the Perception
of Speech and Visual Form, W. Walthen-Dunn, ed., 1967.
[2] H. Blum. Biological Shape and Visual Science: Part I. J. Theoretical Biology, 1973, 38:205-287.
[3] . Ma, M. Sonka. A Fully Parallel 3D Thinning Algorithm and Its Applications. Computer
Vision and Image Understanding, 1996, 64(3):420-433.
[4] C. Pudney. Distance-Ordered Homotopic Thinning: A Skeletonization Algorithm for 3D Digital
Images. Computer Vision and Image Understanding, 1998, 72(3):404-413.
[5] Che Wu-jun, Yang Xun-nian, Wang Guo-zhao. A Dynamic Approach to Skeletonization. Journal
of Software, 2003, 14(4):818-823(in Chinese).
[6] . Niblack, . Gibbon, and . Capson. Generating Skeletons and Centerlines from the
Distance Transform. CVGIP: Graphical Models and Image Processing, 1992, 54(5):420-437.
[7] Y. Ge, . Fitzpatrick. On the Generation of Skeletons from Discrete Euclidean Distance Maps.
IEEE Trans. on Pattern Analysis and Machine Intelligence, 1996, 18(11):1055-1066.
[8] I. Bitter, . Kaufman, M, Sato. Penalized-Distance Volumetric Skeleton Algorithm. IEEE
Trans. on Visualization and Computer Graphics, 2001, 7(3):195-206.
[9] Y. Zhou, . Efficient Skeletonization of Volumetric Objects. IEEE Trans. on
Visualization and Computer Graphics, 1999, 5(3):196-209.
[10] S. Svensson, G. Borgefors. Digital Distance Transform in 3D Image Using Information from
Neighborhoods up to 5×5×5. Computer Vision and Image Understanding, 2002, 88:24-53.
[11] T. Saito, J. Toriwaki. New Algorithm for Euclidean Distance Transformation of an n-Dimensional
Digitized Picture with Applications. Pattern Recognition, 1994, 27(11):1551-1565.
附中文参考文献：
[5] 车武军，杨勋年，汪国昭。动态骨架算法。软件学报，2003，14（4）：818－823。
中国科技论文在线
