一种基于空间中轴的骨架提取算法

数字图像处理之【骨架抽取篇】骨架抽取把一个平面区域简化成图(Graph)是一种重要的结构形状表示法。

利用细化技术得到区域的细化结构是常用的方法。

因此，寻找二值图像的细化结构是图像处理的一个基本问题。

在图像识别或数据压缩时，经常要用到这样的细化结构，例如，在识别字符之前，往往要先对字符作细化处理，求出字符的细化结构。

骨架便是这样的一种细化结构，它是目标的重要拓扑描述，具有非常广泛的应用。

下面首先对数字图像细化概念做简要介绍。

许多数学形态学算法都依赖于击中/击不中变换。

其中数字图像细化，便是一种最常见的使用击中/击不中变换的形态学算法。

对于结构对B=(B1，B2)，利用B细化X定义为即X郦为X与X连的差集。

更一般地，利用结构对序列，，…，迭代地产生输出序列随着迭代的进行，得到的集合也不断细化。

假设输入集合是有限的(即N为有限)，最终将得到一个细化的图像。

结构对的选择仅受结构元素不相交的限制。

事实上，每一个(i=1,2,…,N)都可以是相同的结构对，即在不断重复的迭代细化过程使用同一个结构对。

在实际应用中，通常选择一组结构元素对，迭代过程不断在这些结构对中循环，当一个完整的循环结束时，如果所得结果不再变化，则终止迭代过程。

骨架还可以用中轴表示。

设想在t=0时刻，将目标边界各处同时点燃，火的前沿以匀速向目标内部蔓延，当前沿相交时火焰熄灭，火焰熄灭点的集合就构成了中轴。

图3(a)是这个过程的图示。

另外一种定义骨架的方法使用了最大圆盘概念：目标X的骨架由X内所有最大内切圆盘的圆心组成，如图3(b)、(c)所示。

最大圆盘不是其他任何完全属于X的圆盘的子集，并且至少有两点与目标边界轮廓相切。

骨架的每个点都对应一个相应的最大圆盘和半径r。

最大圆盘定义的骨架与火种方式定义的骨架除在某些特殊情况下端点处存在差异外，绝大多数情况下都是一致的.按照最大圆盘定义骨架的方式，在欧氏二值图像的内部任意给定一点，如果以该点为圆心存在一个最大圆盘，其整个盘体都在图像的内部，且至少有两点与目标边界相切，则该点便是骨架上的点。

中轴变换方法、细化方法和形状分解方法中轴变换方法、细化方法和形状分解方法是数字图像处理中常用的一类方法，用以提取图像中的特征信息。

本文将分别介绍这三种方法，并探讨它们的应用领域和优缺点。

一、中轴变换方法中轴变换（skeletonization）是一种将二值图像转换为其骨架的方法。

所谓骨架，是指将图像中不重要的细节去除后，将图像中的重要特征用线条表示出来。

中轴变换方法通过迭代地操作，将图像中的对象细化为一系列细线。

中轴变换方法的主要应用领域是图像分析和形状识别。

通过获得图像的中轴线，可以进一步分析图像中的形状特征，如曲线的拟合和分割等。

中轴变换方法还可以用于图像的压缩和特征提取，有效地减小图像数据的存储量，同时保留图像的基本特征。

中轴变换方法的优点是能够提取出图像的骨架结构，有助于进一步分析和处理图像。

然而，中轴变换方法在处理复杂图像时可能会存在问题，例如对于细长和叉状物体的处理效果不佳。

二、细化方法细化（thinning）方法是一种将二值图像细化为其最细表示的方法。

所谓最细表示，是指在保持图像中的特征完整性的前提下，将图像中的线条尽可能细化。

细化方法的主要应用是图像的形状分析和对象识别。

通过细化处理，图像中的细线可以更好地表示线条的方向和形状，有助于进一步分析图像中的几何特征和形状信息。

细化方法还可以用于图像的特征提取和匹配，例如指纹识别和虹膜识别等。

细化方法的优点是能够保持图像中的主要特征，并且能够有效地提取线条的方向和形状信息。

然而，细化方法在处理图像中的噪声和连接区域时可能会丢失细节信息，导致图像处理结果不准确。

三、形状分解方法形状分解（shape decomposition）方法是一种将复杂形状分解为简单形状的方法。

通过将图像中的形状划分为几个基本的形状单元，可以更好地理解和描述整个形状。

形状分解方法的主要应用是形状分析和形状识别。

通过形状分解，可以将复杂形状拆分为简单的几何形状，有助于进一步分析形状的几何特征和拓扑关系。

迭代骨架化算法在汉字图像识别中的分析与应用王晓静;吴亚坤;毛红艳;田宝勇【摘要】图像骨架是一种图像对象结构的表示方法,在基于汉字的图像检索中起着举足轻重的作用.研究了几种经典的图像骨架提取算法Zhang-Suen,Rosenfld和Pavlidis.通过程序实现算法分析和比较不同汉字图像骨架的提取效果,并给出针对识别不同的汉字图像应选择的最佳图像识别算法.【期刊名称】《辽宁大学学报（自然科学版）》【年(卷),期】2013(040)003【总页数】6页(P227-232)【关键词】数字图像处理;骨架;Zhang-Suen;Rosenfeld;Pavlidis【作者】王晓静;吴亚坤;毛红艳;田宝勇【作者单位】辽宁大学计算中心,辽宁沈阳110036;辽宁大学计算中心,辽宁沈阳110036;沈阳工程学院信息系,辽宁沈阳110136;辽宁大学计算中心,辽宁沈阳110036【正文语种】中文【中图分类】TP391.410 引言随着计算机的广泛应用和互联网的发展，大量信息以图像的形式出现.在庞大的图像信息量面前，检索成为关键.用关键字的方式检索需要对图像进行人工标注，检索的效率较低而且难以说明图像中所有的信息.在这种情况下，基于内容的图像检索(CBIR)应运而生.它通过颜色、纹理、形状、空间关系等图像特征进行检索.其中基于空间关系的图像检索是利用图像中对象的空间拓扑结构进行检索的.图像中对象的空间拓扑结构表示方式有多种，其中图像骨架是十分有效的一种.图像骨架广泛应用于生物形状描述、模式识别、工业检测、定量金相以及图像压缩编码等领域.利用图像骨架来表示图像的空间关系最早是由Blum[1]提出来的，他使用中轴的概念，用中轴来表示图像的骨架.对于骨架提取算法的一般要求为：连通性，图像骨架连通性必须与原图像保持一致；细化性，图像骨架线条宽度应该为1个像素；中轴性，图像骨架应尽可能是原图的中心线；保持性，图像骨架应尽可能保持原图的细节特征；快速性，算法的速度应尽可能快.在汉字识别领域骨架提取是一个重要研究课题.利用骨架提取技术可以去除图像的冗余信息，骨架提取算法按照实现思路来划分为迭代的骨架提取算法和非迭代的骨架提取算法.本文主要探讨几种经典的迭代骨架提取算法在汉字图像骨架提取中的应用.1 图像骨架提取的几种模型(1)火种传播：在t=0时刻,将图像边界上所有的点同时点燃,火焰以相同速度向图像内部蔓延.当波前相遇时，火焰熄灭，火焰熄灭处所有点的集合就构成了中轴(即骨架SKF，如图1所示).(2)最大圆盘：设D为图像S中的一个内切圆盘，即至少有两点与图像边界相切，如果D不是图像S内部任何其它圆盘的子集，则称为最大圆盘(如图2所示).此时骨架SKM，可定义为图像内部所有最大圆盘圆心的集合.(3)剥皮法：从线的边沿开始，每次剥掉宽度为一个像素的一层，直到最后留下彼此连通的由单个栅格组成的单像素的线条.由于一条线在不同位置可能有不同的宽度，故在剥皮过程中必须注意一个条件，即决不允许剥去会导致曲线不连通的像素.其基本原理是，凡是去掉后不会影响原栅格影像拓扑连通性的像元都应去掉；反之，则应保留.图1 火种传播骨架图2 最大圆盘骨架图3 像素P的8邻域2 几种经典基于剥皮原理的图像骨架提取算法的比较2.1 模板匹配法的原理如图3所示，与像素P相邻的8个像素称为P的8邻域，可以用Pi表示(1≤i≤8).P1,P3,P5,P7称为4-邻域点，P1,P2,P3,P4,P5,P6,P7,P8称为8-邻域点.模板匹配法的原理一般利用一个3*3的模板，在图像中进行搜索，判断像素点是否有与设计的模板相同的情况存在.根据预先设计的与模板相对应的操作(如保留，删除等)进行操作.Rosenfeld算法，Zhang-Suen算法，Pavlidis算法等都是基于模板匹配的骨架提取算法.2.2 Rosenfeld算法Rosenfeld算法[2]是Stefanelli R,Rosenfeld A.在1971年提出的一种并行模板匹配骨架提取算法.该算法定义了最终像素条件，通过标记边缘像素判断它是否满足最终像素条件来决定像素的删除和保留(最终关系像素条件如图4所示，标记为X的点至少有一个为黑点，Y同X，黑点表示前景点，白点表示背景点，灰色的点表示前景点和背景点均可).若像素点满足最终像素条件则保留，否则删除.它的一个处理周期包括4个子步骤，即从上、下、左、右4个方向完成一次边缘像素的删除.重复这个过程，直到在某一次循环中，没有点被删除，则骨架提取结束.图4 最终像素关系Rosenfeld算法骨架提取的过程可描述为：while(仍有点可删除){对于前景点P进行查找；若上边缘点的8邻域与模板2，3，4，5，6，7，8匹配，则继续扫描；若下边缘点的8邻域与模板1，3，4，5，6，7，8匹配，则继续扫描；若左边缘点的8邻域与模板1，2，4，5，6，7，8匹配，则继续扫描；若右边缘点的8邻域与模板1，2，3，5，6，7，8匹配，则继续扫描；否则，标出此点以备删除；}2.3 Zhang-Suen算法Zhang-Suen[3]算法是由Zhang T Y 在1984年提出的一种并行模板匹配骨架提取算法.该算法是一种基于删除的骨架提取算法.该算法定义了删除点的条件(如图5所示)，对于满足删除点条件的进行删除.图5中，其中一个黑色的点，或者两个灰色的点是背景色，则删除该点.算法的一个处理周期分为2个子步骤，第一个步骤按照图5(1)模板进行匹配，第二个步骤按照图5(2)模板进行匹配.重复这个处理过程，直到在某一次循环中，没有点被删除，则骨架提取结束.Zhang-Suen算法骨架提取的过程可描述为：图5 Zhang-Suen细化方法的模板while(仍有点可删除){对于前景点P进行查找；若边缘点的8邻域与模板1匹配，则删除该点；若边缘点的8邻域与模板2匹配，则删除该点；}2.4 Pavldis算法Pavlidis[4]算法是由Pavlidis T 在1982年中提出的一种基于保留的模板匹配骨架提取算法.该算法定义了保留点的条件(如图6所示)，如果满足保留点的条件则保留该点，否则删除.算法处理周期分为四次进行.首先判断右侧边界的像素与保留点条件的模板进行匹配.用同样的操作按顺序对上边界、左边界、下边界进行处理.若在四个子周期中都没有可删除的像素，则骨架提取结束.Pavlidis算法骨架提取的过程可描述为：while(仍有点可删除){对于前景点P进行查找；若右边缘点的8邻域与模板1，2，3，4，5，6匹配，则继续扫描；若上边缘点的8邻域与模板1，2，3，4，5，6匹配，则继续扫描；若左边缘点的8邻域与模板1，2，3，4，5，6匹配，则继续扫描；若下边缘点的8邻域与模板1，2，3，4，5，6匹配，则继续扫描；否则，标出此点以备删除；}图6 保留点的条件(x,y分别至少有一个为“1”)2.5 其它骨架提取算法Hilditch[5]算法的删除条件与Pavlidis算法的保留条件存在一种互补关系，由这两种算法得到的图像骨架类似，在此不再详细讨论.Naccache[6]错误！未找到引用源.算法通过在选定点边界点的条件下，判断其是否为断点、端点或在宽度为2个像素被删除的点.若满足则标记为安全点，否则标记为删除点.扫描结束后，删除标记为删除点的像素点.如此循环，直到没有删除点标记出现.Deutsch[7]算法是一种基于删除的骨架提取算法.它设定了一些删除的条件，如果像素点满足删除条件则删除该像素否则保留.它利用两个子循环交替进行迭代完成骨架提取过程.2.6 骨架提取算法的比较可以看出，这几种骨架提取算法都是通过一系列的条件匹配，不断的进行迭代直到得到最终的骨架.这些算法中，有的通过设定一些条件进行骨架提取，如Hilditch，Deutsch算法，有的则通过一系列的模板进行匹配进行骨架提取，如Zhang-Suen，Rosenfeld.Hilditch，Deutsch，Zhang-Suen是基于删除的骨架提取算法，而Rosenfeld，Pavlidis，Naccache是基于保留的骨架提取算法.Rosenfeld，Pavlidis，Naccache虽然对于保留条件的叫法不一(Rosenfeld称为最终像素条件，Pavlidis称为保留点条件，Naccache称为安全点条件)，但是本质相同.这几种骨架提取算法都是基于剥皮思想的算法，即通过层层剥去边缘像素点，最终得到单像素连通的图像骨架.这些算法迭代的次数不同，Pavlidis和Rosenfeld是由4个子循环组成一个迭代周期，Deutsch和Zhang-Suen则由2个子循环组成一个迭代周期.迭代的过程不同也使得它们的模板的数量和分类不一，Zhang-Suen的模板最少使得它比其它的算法的速度都要快.算法的整体效率受模板数量和迭代次数的影响，当模板数量较少，同时迭代次数较少时算法的整体速度便较快.3 实验结果与分析本文对Hilditch ，Pavlidis，Rosenfeld，Naccache，Deutsch，Zhang-Suen骨架提取算法分别在中国书法，手写文字，印刷文字几个方面(图7)进行了比较(分别见图8，9，10，11).由实验结果可以看出Rosenfeld，FasThin算法在中国书法，手写文字，印刷文字几个方面处理的效果都要好于Hilditch ，Pavlidis,Naccache，Deutsch算法.Hilditch算法和Pavlidis算法得到的图像骨架类似，对手写文字和印刷文字的处理效果较好，但在中国书法的处理中不能像Rosenfeld，Zhang-Suen算法那样精准的反映线条的交叉.Pavlidis算法得到的图像骨架连通性要差于Hilditch算法，如图11书法图像的处理中Pavlidis算法得到的图像骨架出现了明显的断裂. a.手写文字 b.印刷体汉字 c.印刷体英文 d.书法图7 实验图a.Hilditchb.Pavlidisc.Rosenfeldd.Naccachee.Deutschf.Zhang-Suen图8 手写文字a.Hilditchb.Pavlidisc.Rosenfeldd.Naccachee.Deutschf.Zhang-Suen图9 印刷体汉字a.Hilditchb.Pavlidisc.Rosenfeldd.Naccachee.Deutschf.Zhang-Suen图10 印刷体英文a.Hilditchb.Pavlidisc.Rosenfeldd.Naccachee.Deutschf.Zhang-Suen图11 书法由Rosenfeld算法和Zhang-Suen算法得到的骨架不仅线条流畅，没有多余的枝杈，而且能够精准的反映出线的交叉.Zhang-Suen算法能够反映出较小的必要枝杈.Naccache算法处理手写文字和印刷体英文的效果较好，但对于印刷汉字和中国书法处理效果一般，产生了过多的枝杈.Deutsch算法处理手写文字的效果较好，但在处理印刷体汉字和中国书法时效果较差，不仅会产生多余的枝杈，而且会产生一定程度的畸变.4 结论本文介绍一些骨架提取的模型，并对几种经典的基于剥皮思想的骨架提取算法从原理和实验的角度进行了分析和比较.Zhang-Suen算法的执行效率要高于其它算法，Rosenfeld，Zhang-Suen算法得到的图像骨架能够较好的表现出汉字的拓扑特征，而且图像骨架没有多余的枝杈，较为适合应用在汉字图像的骨架提取.参考文献：[1] Blum H.A transformation for extracting new description ofshape.In:Wathen-Dunn W ed.Model for the Perception of Speech and Visual[M].Cambride,Massachusetts:MIT Press,1967:362-380.[2] Stefanelli R，Rosenfeld A.Some parallel thinning algorithms for digital pictures[J]. Journal of ACM，1971(2).[3] Zhang T Y.A fast parallel algorithm for thinning digitalpatterns[J].Communications of ACM,1984,27(3):236-239.[4] Pavlidis T.Algorithms for Graphics and lmage Procesing[M].Rockville：Computer Science Press，1982.[5] Hilditch C J.Linear Skeletons from Square Cupboards[J].Machine Intelligence,1969,4:403.[6] Naccache N J,Shinghal R.SPTA:A proposed algorithm for thinning by binary patterns[J].IEEE Transactions on System,ManCybernetics,1984,14(3):409-418.[7] Deutsch E S.Thinning algorithms on rectangular,hexagonal and triangular arrays[J].Communications of the ACM,1972,15(9):827-837.。

第34卷第12期自动化学报Vol.34,No.12 2008年12月ACTA AUTOMATICA SINICA December,2008一种改进的骨架曲线串行多边形近似算法吕哲1王福利1,2常玉清1,2刘阳1摘要常见骨架提取算法对复杂多变的目标边缘具有较强的敏感性,提取出的骨架曲线结构相对复杂,数据量仍然较大.针对这一问题,提出了一种新的骨架曲线多边形近似算法.该算法结合骨架曲线的特点,在传统串行多边形近似算法的基础上引入了平滑度保持、结构特征保持以及拓扑特征保持等约束条件,既较好地保留了原始骨架的主要拓扑结构特征,又有效地简化了骨架曲线的结构,进一步压缩了数据.仿真研究证明了该方法的有效性.关键词骨架,多边形近似,形状分析中图分类号TP391An Improved Sequential Polygonal Approximation Algorithm on Skeleton Curves LV Zhe1WANG Fu-Li1,2CHANG Yu-Qing1,2LIU Yang1Abstract High sensitivity of most skeletonization algorithms to the variable edge features makes the structure of the extracted skeletons relatively complicated and the corresponding data excessive.To solve this problem,a novel skeleton polygonal approximation algorithm is proposed.Several impactful improvements such as smoothness-preservation, structure-preservation,and topo-preservation have been proposed to characterize the main structure of the original skeleton as well as predigest other trivial structure for compression of redundant information.The new algorithm is thus superior to the conventional one.Simulation tests have verified the effectiveness of the algorithm.Key words Skeleton,polygonal approximation,shape analysis在图像几何形态分析及其相关领域中,骨架具有非常广泛的应用.利用骨架表示原始图像,可以在保持图像重要拓扑结构特征的前提下,减少冗余信息.因此,基于骨架的目标形状表示和识别技术已经成为模式识别领域的重要研究内容,被广泛应用于自动字符识别、指纹识别、工件识别以及医学图像分析等领域[1−2].近年来,针对目标骨架提取问题的研究已经十分深入,新算法及改进算法层出不穷.总的来说,这些算法可以分为两大类型[3]:1)细化算法,这类算法以焚烧草地定义为基础,通过不断地剥离边缘点来获取目标的骨架[4−5];2)中轴变换法,这类算法以最大圆盘定义为基础,通过距离变换或形态学变换搜寻目标中最大圆盘的圆心点,并由此构成目标的骨架[6−7].不同类型的骨架提取算法虽然各有优缺点,但收稿日期2007-10-08收修改稿日期2008-01-14Received October8,2007;in revised form January14,2008国家重点基础研究发展计划(973计划)(2002CB312201),国家自然科学基金(60374003)资助Supported by National Basic Research Program of China(973 Program)(2002CB312201)and National Natural Science Foun-dation of China(60374003)1.东北大学信息科学与工程学院沈阳1100042.东北大学流程工业综合自动化教育部重点实验室沈阳1100041.School of Information Science and Engineering,Northeast-ern University,Shenyang1100042.Key Laboratory of Inte-grated Automation of Process Industry,Ministry of Education, Northeastern University,Shenyang110004DOI:10.3724/SP.J.1004.2008.01467都能够较好地实现对目标骨架特征的提取.然而,在实际应用中,多数骨架提取算法获得的骨架曲线都存在一个共同的问题,那就是它们对目标形状的细节变化比较敏感,自身结构较为复杂,数据难以得到有效压缩.这一不足一方面影响了骨架曲线对目标主要拓扑特征的描述能力,另一方面占用了大量的存储空间和计算时间.因此,有必要对骨架曲线的结构简化问题进行研究.线性近似(Linear approximation,LP),或称为多边形近似(Polygonal approximation,PA),是一种行之有效的离散曲线简化算法[8].该算法通过控制近似偏差,以分段线性线段或圆弧来代替原曲线,从而达到简化曲线结构、压缩曲线数据的目的.目前,多边形近似算法在目标边缘表示等领域已经得到了广泛应用.借鉴多边形近似算法在离散曲线简化问题中的成功应用,结合骨架曲线自身结构特点,本文对传统的多边形近似算法进行了改进,成功地将其应用于骨架曲线结构简化问题,使得简化后的骨架曲线在有效保留目标主要拓扑结构特征的同时,信息量得到进一步压缩,能够更好地满足实际应用的要求,具有较为重要的理论及现实意义.1传统算法简介在图像处理、计算机图形、模式识别及数据压缩等诸多领域中都存在离散曲线结构简化问题.通1468自动化学报34卷常,在此类应用中,人们希望使用某些简单的特征来近似复杂的图形或几何目标,以达到数据压缩的目的.多边形近似是曲线结构简化算法中应用最为广泛的一种.该方法通过控制近似偏差,在偏差允许范围内不断减少数据点(或称为控制点),用分段线性线段或圆弧来代替原始曲线,简单有效地实现了对离散曲线的近似.近似后的分段线性曲线能够携带尽可能多的原始曲线信息,可以很好地描述原始曲线的形状结构特征.多边形近似算法定义如下[9].定义1.令C为任意一条离散曲线上的N(N>1)个点组成的序列,C={P i(x i,y i)|i=1,2,···,N}.其中(x i,y i)为序列中第i个点P i的坐标.定义2.线段S由C中两个不同的点P i和P i+k唯一确定,S=P i P i+k.定义3.线段S相对于原始离散曲线的近似偏差E(S)为E(S)=maxi≤t≤i+k E Pt(S)=maxi≤t≤i+kDist(P t,S)(1)其中,Dist(P t,S)为点P t,t∈(i,i+k)到线段S 的垂直距离的平方.Dist(P t,S)=((x t−x i)(y i+k−y i)−(y t−y i)(x i+k−x i))2(x i−x i+k)2+(y i−y i+k)2(2)定义4.任意曲线的多边形近似曲线P A满足如下条件:1)P A={S1,S2,···,S Card(P A)};2)∀S n,S n+1∈P A,S n的最后一个点为S n+1的第一个点;3)Card(P A)为P A中线段的数目;4)令SP A为全部多边形近似曲线P A的集合.给定一个最大误差阈值,则多边形近似问题可以描述为求取线段序列P A1={S1,S2,···, S Card(P A)},使得∀E(S n)<MaxError,n∈(1,Card(P A1))(3)且对于满足Card(P A2)<Card(P A1)的∀P A2∈SP A有∃E(S m)>MaxError,m∈(1,Card(P A2))(4)常用多边形近似算法可以分为三大类,即串行方法(Sequential approach)、分裂合并法(Split-and-merge approach)及启发式搜索法(Heuristic-search approach)[10].其中,分裂合并法对初始线段的选择有很强的依赖性,而启发式搜索法计算复杂度高,计算效率较低.串行方法虽然难以保证最优近似结果,但对于结构相对简单的骨架曲线仍不失为一种较好的选择.此外,串行方法还具有算法简单、计算速度快等优点.因此,本文采用串行多边形近似算法对骨架曲线进行简化处理.2改进算法骨架曲线通常被用于目标形状的表示和描述,人们更为关心的是其携带的目标的主要拓扑结构信息,而非目标形状的细节变化信息.然而,不论是采用细化算法,还是采用中轴变换方法,骨架曲线的获取在很大程度上都依赖于目标的边缘信息,最终获得的骨架曲线对边缘的细节变化有很强的敏感性.这一方面使得骨架曲线的结构不够平滑、存在噪声,影响其对目标主要结构特征的表示;另一方面使得骨架曲线中含有大量不重要、甚至错误的细节信息,给之后的计算和存储带来困难.基于此,本文提出了一种新的骨架曲线串行多边形近似算法,通过平滑度保持、结构特征保持和拓扑稳定性保持等约束,提高了串行近似算法保留原始骨架主要拓扑结构特征的能力,降低了骨架曲线对边缘细节变化的敏感程度,达到了简化骨架结构、突出主要信息以及进一步压缩数据的目的.改进的骨架曲线串行多边形近似算法包括控制点初始化和近似计算两个主要步骤.2.1控制点初始化对于给定的骨架曲线Skel有Skel={P i(x i,y i)|i=1,2,···,n}(5)其中,n为骨架曲线上骨架点的数目,P i为第i个骨架点,且其坐标为(x i,y i).初始控制点(Dominant point,DP)选择的依据是由这些控制点确定的近似曲线P A与原曲线之间不能存在任何偏差[10].显然,最简单的方法就是将骨架曲线上的所有点均初始化为控制点.但是,这种方法不可避免地包含了大量的冗余信息.相比于上述方法,将骨架曲线上的所有拐点(即曲线走向发生变化的点)初始化为控制点是一种更为有效的方法.它不仅能够保证近似曲线与原曲线之间的零偏差,而且大量减少了冗余信息.因此,本文使用该方法对骨架曲线的控制点进行初始化.如图1(见下页)所示,使用Freeman链码[11]可以方便地实现拐点的确定.对于骨架曲线上的任意点P i,依据其相邻骨架点P i+1在邻域内所处的位置,可以给P i分配一个0到7之间的整数编码c i.那么当任意骨架点P i的链码与其前一个骨架点12期吕哲等:一种改进的骨架曲线串行多边形近似算法1469P i−1的链码不同时,即c i=c i−1(6)则骨架曲线的走向在该点处发生变化,点P i即为曲线的一个拐点.骨架曲线上的所有拐点构成了曲线的初始控制点.图1Freeman链码Fig.1Freeman chain code此外,由于骨架曲线上的端点和连接点对于骨架曲线的拓扑结构的稳定起着重要的作用,因此,无论这些点是否为骨架曲线的拐点,都应将其初始化为控制点.对于端点和连接点的定义,在下文中将详细叙述.2.2近似计算在完成控制点初始化之后,近似算法从骨架曲线的端点开始以串行方式进行曲线近似.对于控制点DP i,为了评价该点与其后第k个控制点DP i+k 之间的连线S n=DP i DP i+k能否用于近似此两点之间的原始骨架曲线,需要计算两点之间各控制点的近似偏差E DPt(S n),t∈(i,i+k),并求出最大近似偏差E(S n).当E(S n)小于预先设定的最大偏差阈值MaxError时,继续判断下一个控制点DP i+k+1;否则,当E(S n)>MaxError时,将上一个控制点DP i+k−1与起始控制点DP i之间的骨架曲线近似为此两点的连线S n=DP i DP i+k−1,并将DP i+k−1设为下一搜索的起始点.但是在实际应用中,若直接根据式(1)和(2)计算近似偏差E DPt(S n)和最大近似偏差E(S n),会导致近似后的骨架曲线无法保持良好的平滑度和稳定的结构特征(其原因将在第3节中详细分析).为此,本文在传统算法的基础上引入了三个约束条件.1)平滑度保持约束.为了保持骨架曲线的平滑程度,避免骨架曲线上的噪声突起给曲线近似带来的影响,改进算法在E DPt(S n)的计算过程中引入了修正项Modify(DP t,S n),用于对原始骨架中的噪声偏差进行约束.E DPt(S n)=Dist(DP t,S n)+Modify(DP t,S n)(7)2)结构特征保持约束.能够准确描述目标的主要结构特征是近似后骨架曲线必须具备的能力.为此,改进算法进一步引入了惩罚项P enalty(S n),对同方向偏差情况下的近似曲线跨度进行约束,避免了近似曲线偏离原始骨架曲线情况的发生.E(S n)=maxi≤t≤i+kE DPt(S n)+P enalty(S n)=maxi≤t≤i+k(Dist(DP t,S n)+Modify(DP t,S n))+P enalty(S n)(8)其中,近似偏差E DPt(S n)是采用带有平滑度保持约束的式(7)得到的,可见改进后的最大偏差E(S n)同时结合了平滑度保持约束和结构特征保持约束.3)拓扑稳定性保持约束.由于骨架曲线由主干和连接于其上的分支构成,所以曲线近似过程还必须能够保证骨架曲线各分支与主干之间的拓扑关系不发生变化.为此,本文还引入了确保近似曲线拓扑不变性的拓扑稳定性保持约束.下面分别对上述三个约束条件进行详细分析.3约束条件分析3.1平滑度保持如图2所示,当目标边缘存在细微变化时,获得的骨架曲线往往会产生相应的突起,从而降低了骨架曲线的平滑程度.在对该骨架曲线进行多边形近似时,为了克服这些突起的影响,提高骨架曲线的平滑程度,往往需要设定相对较大的偏差阈值.但是,较大的阈值又会影响其他部分曲线的近似效果,甚至导致近似后曲线不能有效保留原始曲线的结构特征.图2边缘细微变化对骨架的影响Fig.2Influence of the edge details在研究过程中发现,突起的顶角的角度和组成突起的骨架点数目是这类由边缘的细微变化引起的突起区别于骨架曲线上正常结构变化的两个主要特征.因此,在计算突起处的近似偏差时,若能依据这两个特征进行修正,就能够有效地抑制其对骨架曲线近似过程的影响.文献[12]给出了一种确定曲线上某一拐角的支持区域的方法.借鉴这一方法,能够十分方便地确定组成突起的骨架点的数目.设突起的顶点为骨架点P i,定义连接P i两侧的1470自动化学报34卷骨架点P i−k和P i+k的线段的长度为l i,k=Pi−kP i+k(9)令d i,k为点P i到线段P i−k P i+k的垂直距离.从k=1开始,计算l i,k和d i,k,直到l i,k≥l i,k−1(10)或d i,k l i,k ≥d i,k+1l i,k+1,d i,k>0(11)d i,k l i,k ≤d i,k+1l i,k+1,d i,k<0(12)此时,以P i为顶点的突起所包含的骨架点的集合为A(P i)={P i−k,···,P i,···,P i+k|a or b}(13)集合A(P i)中骨架点的数目为Card(A(P i))=2k+1(14)在确定了组成突起的骨架点的集合的同时,还可以进一步确定出该突起的顶角的角度.令a ik=(x i−x i+k,y i−y i+k)(15)b ik=(x i−x i−k,y i−y i−k)(16)则顶角的角度为angle Pi =arccosa ik·b ik|a ik||b ik|(17)此时,本文利用各控制点DP i的Card(A(DP i))和angle Pi对该点处的近似偏差进行如下修正E DPt(S n)=Dist(DP t,S n)+Modify(DP t,S n)(18)Modify(DP t,S n)=c·1Card(A(DP t))(1+cos(angle DPt))(19)其中,c为负实数,用于保持尺度一致.当控制点DP i位于由边缘细微变化引起的突起上时,Card(A(DP i))和angle Pi均较小,修正作用较强,只需用较小的偏差阈值即能实现对该突起的平滑.当DP i位于骨架曲线上正常的结构变化处时, Card(A(DP i))和angle Pi均较大,修正作用很弱,不会影响骨架曲线正常部分多边形近似的效果.这样,新算法就能够在平滑骨架曲线的同时,有效降低骨架对边缘变化的敏感性给骨架曲线的多边形近似带来的影响.3.2结构特征保持对骨架曲线进行多边形近似处理的目的就是在保留原始骨架曲线主要结构特征的前提下,尽可能地简化骨架曲线的结构,压缩数据.因此,能够更为准确地保留原始骨架曲线的主要结构特征,这对于骨架曲线多边形近似算法是十分重要的.使用传统串行多边形近似算法虽然能够有效简化骨架曲线,但是在某些情况下,其近似结果不能准确地反映原始骨架的主要结构特征.例如,图3中曲线a、曲线b分别为不同偏差阈值下多边形近似的结果,其中MaxError a<MaxError b.可以看出,虽然曲线b较曲线a有更少的控制点,数据量也更小,但是,其主要结构与曲线a的相似性较差,不能很好地反映原始曲线的主要结构特征.其原因主要是传统串行多变形近似算法只能保证控制点到近似后曲线的最大偏差不超过预先设定的幅值,却无法防止多数控制点位于近似曲线的同一侧的情况出现.随着同侧偏差的累积,近似曲线逐渐偏离原始曲线.图3引入结构稳定性约束前后近似结果Fig.3Approximated skeletons with and withoutstructure-preservation为了解决这一问题,本文在计算最大近似偏差E(S n)时引入了一个惩罚项,用于控制近似曲线与原始骨架曲线的偏离程度.E(S n)=maxi≤t≤i+kE DPt(S n)+P enalty(S n)(20) P enalty(S n)=r·i+kt=isgn(DP t)E DPt(S n)(21) sgn(DP t)=1,y t−mx t−n≥0−1,y t−mx t−n<0(22)其中,r为正实数,用于保持尺度一致,m和n 分别为S n=DP i DP i+k的斜率和截距,惩罚项P enalty(S n)为线段S n包含的区域内所有控制点12期吕哲等:一种改进的骨架曲线串行多边形近似算法1471近似偏差的累计值.当多数控制点均位于近似曲线同一侧时,累计值将快速增大,最大近似偏差E(S n)也随之快速增大.在最大偏差阈值的限制下,近似曲线的跨度将得到控制,从而有效防止了近似曲线与原始骨架曲线的偏离程度进一步扩大.由图3可以看出,在引入惩罚项后,所得近似曲线c与曲线a更为相似,能够更好地反映原始曲线的主要结构特征.3.3拓扑稳定性保持在对骨架曲线进行多边形近似处理时,除了要保证近似后的骨架曲线能够准确地反映原始骨架的主要结构特征,还需保证骨架曲线各分支与主干之间的拓扑关系不发生变化.由于骨架曲线中的端点和连接点决定了该骨架曲线的拓扑,因此,只要保证端点和连接点的位置不发生变化,即可有效保持骨架曲线的拓扑稳定性.定义骨架点P i的8邻域内像素沿逆时针方向由0到1及由1到0变换的次数为[13]X R(P i)=8i=1|x i+1−x i|(23)定义骨架点P i的8邻域内的骨架点的数目为X H(P i)=8i=1x i(24)基于此,本文对骨架曲线中的端点EP(Skel)和连接点JP(Skel)分别作如下定义EP(Skel)=P∈Skel{X R(P)≤2and X H(P)≤2}(25) JP(Skel)=P∈Skel{X R(P)≥6and X H(P)≥3}(26)此时,可将拓扑稳定性保持过程表述如下:1)在初始化控制点时,将骨架曲线上的所有端点和连接点初始化为控制点;2)在进行串行近似时,首先将各端点作为起点进行近似搜索;3)当搜索遇到连接点时,无论近似偏差是否超过预先设定的阈值,都强行结束该搜索,转而从另一端点开始进行新的搜索;4)当没有新的未搜索的端点时,将连接点作为搜索起点,沿未搜索过的骨架曲线,进行搜索;5)当搜索遇到连接点时,无论近似偏差是否超过预先设定的阈值,都强行结束该搜索,转而从另一连接点开始进行新的搜索;6)直到对所有骨架曲线均完成近似搜索后过程结束.4实验结果分析及讨论在仿真实验中,本文采用形态学骨架提取算法对二值的Island图像进行骨架提取,然后应用新的骨架曲线多边形近似算法对原始骨架进行简化处理.图4中的圆圈标出了骨架曲线上控制点的位置,这些点是确定该骨架形状所需的最少的数据.由图4(b)可见,简化前骨架曲线的结构仍然较为复杂,数据量较大,不利于表示和存储.而图4(c)、4(d)分别给出了在不同偏差阈值(Threshold,TH)下应用本文提出的骨架曲线多边形近似算法进行简化处理后获得的新骨架曲线.由图可见本文算法在较好地保留原始拓扑结构特征的同时,大大简化了骨架曲线的结构,从而有效地压缩了数据.(a)Island1图像(b)原始骨架曲线(a)Island1(b)Originalskeleton(c)本文算法近似结果(d)本文算法近似结果(T H=5.0)(T H=15.0)(c)Approximated skeleton(d)Approximated skeletonby proposed algorithm by proposed algorithm(T H=5.0)(T H=15.0)图4本文算法在不同偏差阈值下的近似结果Fig.4Approximated skeletons by proposed algorithmwith different thresholds在将多边形近似算法应用于骨架曲线近似问题的过程中,为了解决骨架曲线平滑度保持和结构特征保持的问题,本文分别引入了修正项Modify(DP t,S n)和惩罚项P enalty(S n).图5和图6对改进前后的多边形近似算法在骨架曲线多边形近似问题中的应用效果进行了比较研究.1472自动化学报34卷(a)Island2图像(b)原始骨架曲线(a)Island2(b)Originalskeleton(c)改进前算法近似结果(d)改进后算法近似结果(c)Approximated skeleton by(d)Approximated skeleton byconventional algorithm proposedalgorithm(e)图(c)中虚线框部分放大(f)图(d)中虚线框部分放大(e)Part of (c)(f)Part of (d)图5较小偏差阈值下近似结果比较Fig.5Comparison of the approximated skeletons withlower threshold当预先设定的偏差阈值较小时,传统近似算法虽然能够较为准确地保留原始骨架曲线的拓扑结构特征,但是由于原始骨架曲线中存在由目标边缘细微变化引起的噪声突起(如图5(b)所示),因此,较小的偏差阈值将很难有效地抑制这些不重要的噪声信息(图5(c)、5(e)),从而既影响了近似后曲线的平滑程度,又不能在偏差允许范围内最大程度地压缩数据.为了解决这一问题,本文在计算近似偏差时引入了修正项,使得新算法在相同的偏差阈值下,较好地避免了噪声突起的影响,获得了更为出色的骨架曲线多边形近似结果(如图5(d)、5(f)所示).当预先设定的偏差阈值较大时,人们往往希望在不破坏原始骨架曲线拓扑特征的前提下,尽可能地简化骨架曲线的结构.但是,在实际应用过程中,由于传统多边形近似算法以最大近似偏差作为衡量标准,因此难以防止近似曲线偏于原始骨架一侧的情况出现,使得近似曲线在一些情况下不能正确地反映原始骨架曲线的主要结构特征(图6(a)、6(c)).为此,本文引入了基于累计近似偏差的惩罚项,较好地避免了传统多边形近似算法存在的问题.如图6(b)、6(d)所示,在相同偏差阈值下,改进算法的近似结果能够更好地保留原始骨架曲线的结构特征.(a)改进前算法近似结果(b)改进后算法近似结果(a)Approximated skeleton by(b)Approximated skeleton byconventional algorithm proposedalgorithm(c)图(a)中虚线框部分放大(d)图(b)中虚线框部分放大(c)Part of (a)(d)Part of (b)图6较大偏差阈值下近似结果比较Fig.6Comparison of the approximated skeletons withhigher threshold为了更客观地评价算法性能,本文对改进前后多边形近似算法的骨架曲线近似效果进行了量化分析.在实际应用中,常见的评价多边形近似算法性能的两个指标分别为压缩比(Compression ratio,CR)和误差平方积分(Integral square error,ISE)，其相应的表达式为CR =nnDP (27)ISE =e i (28)其中,CR 为原始曲线骨架点数目n 与近似后曲线控制点数目nDP 的比值;而ISE 为原始骨架曲线所有控制点与近似后骨架曲线间偏差的累计量.但是,这两个指标往往是相互矛盾的.较大的CR 通常伴随着较大的ISE .Sarkar [14]将两个指标结合起来,构成一个更具一般性的指标,即性能指数(Figure of merit,FOM),其表达式为F OM =CR ISE(29)12期吕哲等:一种改进的骨架曲线串行多边形近似算法1473在相同偏差阈值的情况下,如表1所示,当阈值较小时,改进算法在不明显提高累计偏差的同时,获得了更高的压缩比.这里计算的偏差是相对于含有噪声的原始骨架的,从这一层面上说,抑制噪声突起的过程实际上提高了近似骨架曲线相对于真实骨架曲线的精度;当阈值较大时,改进算法在保持较高压缩比的同时,有效地降低了累计偏差,更好地保持了近似曲线的结构稳定性.相应地,改进方法的综合性能指标F OM也优于传统串行多边形近似算法.表1不同偏差阈值下近似结果量化对比Table1Comparison of approximating results withdifferent thresholds方法T H CR ISE F OM57.5122.49000.0612改进前1510.0415.32270.024158.6126.74180.0679改进后159.2179.66120.0514进一步地,文献[15]提出了一种相对量化指标,通过与最优近似结果的比较来衡量串行近似算法的性能.该量化指标基于两个参数:F idelity(偏差度量参数)和Efficiency(压缩比度量参数).F idelity=E optE appr×100(30)Efficiency=N optN appr×100(31)其中,E appr和N appr为待测试串行近似算法的误差平方和ISE及近似后曲线控制点数目nDP.E opt 为相同近似后曲线控制点数目下,最优算法获得的误差平方和,而N opt为相同误差平方和下,最优算法获得的近似曲线的控制点数目.将两个参数结合可以获得一个更为一般的指标Merit=F idelity×Efficiency(32)本文对广泛应用的基于动态规划的Preze and Vidal最优多边形近似算法[16]进行了拓扑稳定性改造,将其应用于骨架曲线的多边形近似计算,并以此为参考来衡量改进前后的串行多边形近似算法.在相同近似后控制点数目的情况下,如表2和表3所示,约束条件的引入使得改进的骨架曲线串行多边形近似算法在基本不增加算法计算代价的同时,获得了更为接近最优近似曲线的近似结果.这也进一步验证了改进方法在近似性能方面的优势.表2改进前后方法的相对评分Table2Relative scores of the conventional algorithmand the improved one方法T H nDP ISE F idelity Efficiency Merit2.036105.3220.741.729.4改进前12.026221.8332.146.238.51.136102.4721.344.430.8改进后15.026179.6639.761.349.3表3计算用时比较(P43.0,512M)Table3Comparison of computation times(P43.0,512M)方法(nDP=26)改进前改进后Perez and Vidal 计算用时(ms)18.419.11071.2通过文中分析和相关仿真研究,我们知道基于寻优的多边形近似算法虽然能够获得较高的近似精度,但其过高的计算用时是多数实际过程不能允许的;而传统串行多边形近似算法虽然计算效率很高,但其容易受到噪声突起和同侧偏差累积的影响,在很多情况下不能准确描述目标的主要拓扑结构特征.相比之下,本文提出的改进的骨架曲线串行多边形近似算法一方面克服了传统方法存在的不足,其相对评分达到或接近文献[15]中列出的多种基于寻优或学习的曲线近似算法的评分,能够较好地满足平滑骨架曲线、保持目标主要拓扑结构特征的要求;另一方面继承了串行方法计算速度快的优点,在算法效率方面明显优于其他曲线近似算法.5结论骨架是目标形状分析、模式识别等领域中应用最为广泛的一种特征表示方式.但是,由于骨架的提取过程对于多变的目标边缘有较强敏感性,所提取的骨架结构相对复杂,给骨架的应用造成了困难.为此,本文提出了一种改进的串行多边形近似算法,通过平滑度保持、结构保持和拓扑保持等约束条件,成功地将传统多边形近似算法应用于骨架曲线的结构简化问题.相关仿真实验表明,应用本文算法得到的简化骨架曲线不仅能够很好地描述原始骨架曲线的主要拓扑结构特征,而且获得了较高的压缩比,有效地减少了骨架曲线中的冗余信息,具有很高的应用价值.References1Cornea N D,Min P.Curve-skeleton properties,applications, and algorithms.IEEE Transactions on Visualization and Computer Graphics,2007,13(3):530−5482Van E M,Macrini D,Telea A,Sminchisescu C,Dickinson S。

MATLAB中的骨架提取代码
在MATLAB中，骨架提取（也称为中轴变换或骨架化）通常用于二值图像，以提取对象的中心线或形状的核心。

MATLAB的Image Processing Toolbox提供了一些函数，如bwmorph和imthin，可以帮助进行骨架提取。

以下是一个简单的示例代码，展示如何使用MATLAB进行骨架提取：
matlab
% 读取二值图像
binaryImage = imread('your_binary_image.png');
% 使用bwmorph函数进行骨架提取
% 'skel', Inf 表示进行无限次骨架提取，直到图像不再变化
skeletonImage = bwmorph(binaryImage, 'skel', Inf);
% 显示原图和骨架图
figure;
subplot(1, 2, 1);
imshow(binaryImage);
title('Original Binary Image');
subplot(1, 2, 2);
imshow(skeletonImage);
title('Skeleton Image');
在这个示例中，your_binary_image.png应替换为您要处理的实际二值图像文件的名称。

这段代码将读取二值图像，然后使用bwmorph函数进行骨架提取，并显示原始图像和提取的骨架图像。

请注意，骨架提取的效果可能因输入图像的质量和特性而异。

您可能需要调整参数或尝试其他方法来获得最佳结果。

matlab骨架提取MATLAB是一种强大而又广泛使用的数学软件，它有能力进行图像处理并提取图像的特征。

其中，骨架提取是一种常见的图像处理技术，可以将图像中的细节信息提取出来并进行分析。

下面介绍MATLAB骨架提取的基本流程及其实现方法。

1.图像预处理图像预处理是骨架提取的重要步骤，它可以使图像中的噪声、模糊和阴影等因素减轻，从而提高骨架提取的效果。

常用的预处理方法有去噪、二值化、边缘检测等。

其中，二值化是将图像转换为黑白的二值图像，边缘检测是将图像中的边缘提取出来。

2.形态学处理形态学处理是指对图像进行膨胀、腐蚀等操作，以改变图像中像素的形状和结构。

在骨架提取中，形态学处理可以设置不同的结构元素，以使得骨架提取结果更加准确。

常用的形态学处理方法有膨胀、腐蚀、开运算、闭运算等。

3.骨架提取骨架提取是指从二值图像中提取出其主干骨架的过程。

常用的骨架提取算法有细化算法、倒数距变换算法、中轴变换算法等。

其中，细化算法是最常用的骨架提取算法之一，它是基于迭代的二值形态学算法，可以得到相对准确的骨架结果。

4.骨架后处理骨架后处理是指对提取得到的骨架进行处理，以得到更加完整的形态信息。

它可以用于去除骨架中的毛刺、孤立点等干扰噪声，也可以用于将骨架与原图像进行匹配，以便于进一步分析和处理骨架。

综上所述，MATLAB骨架提取是一种广泛用于图像处理领域的技术。

它可以帮助用户提取出图像中的主干骨架，并在其基础上进行各种形态分析和处理。

尽管骨架提取存在一些局限性，如难以处理部分凸起物体的骨架、骨架分支过于复杂等问题，但它仍然是一种有效的图像处理技术，可以用于医学图像分析、机器视觉、遥感等领域。