CN103673923A - 基于数字图像处理的曲线纤维网络结构形貌特征测量方法 - Google Patents

Abstract

一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，包括以下步骤：一、纤维网络的骨骼化，得到纤维网络骨架图案；二、纤维网络骨架图案的拓扑特征分类，识别出纤维网络中的单纤维及纤维交叉点；三、删除短纤维，保留主干纤维的网络骨架图案；四、纤维片段重组，识别纤维网络中的交叉、分叉及重叠关系，并重组纤维网络图；通过上述步骤得到亚细胞尺度的曲线纤维网络结构形貌特征。本发明可用于定量分析亚细胞尺度的曲线纤维网络结构，本发明在完成了对曲线网络中纤维结构的识别和分割后，还能够不依赖人为的干预和判断完成对显微图像中曲线网络结构的定量拓扑分析。

Description

基于数字图像处理的曲线纤维网络结构形貌特征测量方法

技术领域

本发明属于生物实验图像数据分析领域，特别是一种应用于亚细胞尺度的基于数字图像处理的曲线纤维网络结构形貌特征测量方法。

背景技术

大量的生物实验已经积累的大量的实验图片，过去的生物学研究以定性为主，一般通过直接观察或简单的量化分析得出研究结论。越来越多的实验图像的产生让我们进入了一个大数据的时代，挖掘数据背后的定量规律成了我们面临的巨大挑战，以数据为驱动的研究将生物学带往定量时代。在医学图像处理蓬勃发展的当下，生物图像处理的发展显得与生物学领域的发展不够匹配。具有曲线纤维网络结构特点的一类生物图像是典型的生物图像，生物体中有大量的网络结构，组织水平有血管网络；细胞水平有神经细胞网络、骨细胞网络；亚细胞水平有细胞骨架网络、核骨架网络等。

亚细胞水平的纤维网络的纤维尺度已经接近光学显微镜的分辨率极限（abby极限），即几百纳米。更小尺度的网络结构需要利用分子成像技术，如光激活定位显微技术(photoactivated localization microscopy,PALM)和随机光学重构显微技术(stochastic optical reconstruction microscopy,STORM)等进行观测。因此基于受限于abby极限的普通光学显微的网络结构中，亚细胞水平的网络结构的分析是最具有挑战性的。发展一套定量的方法来研究这种曲线纤维网络的拓扑结构，离散和提取每一根纤维，分析其拓扑信息，将带来很多直接观察所不能获得的定量信息，更有效地帮助人们找到更精确的量化规律，甚至不易直接察觉的内在关系。这种基于纤维网络结构的图像分析方法将架起一大类生物结构生理功能关系研究的桥梁。

 亚细胞水平的网络结构显微图像往往都非常模糊，并伴随大量的噪声。其原因是显微结构已经接近abby极限，而且胞内的网络一般具有三维结构，显微镜的点扩散函数不仅对聚焦面内信号的模糊产生巨大影响，还将焦外的信号叠加到聚焦面内产生模糊，此外，微小的显微结构直接导致此类图片中的信噪比较低。因此，如何从低信噪比的模糊图像中准确地提取出曲线纤维网络结构并分析其拓扑结构得到每根纤维的拓扑信息是一个高度挑战性的问题。 

现有技术中未见到通过分析单根纤维的形态及纤维间的具体拓扑关系来研究亚细胞纤维网络的拓扑形貌结构。以往研究纤维网络的方法一般都是基于图像分割阈值或轮廓跟踪技术。基于图像分割阈值的方法是通过设定一个全局灰度阈值来确定纤维的中心线进而研究整个纤维网络结构。针对高分辨率相差显微图像发展出了基于能量泛函并结合最优化阈值及水平集函数分割技术的方法来提取微结构轮廓。此外，一种基于方向可调滤波器原理开发出来的曲线纤维结构探测器被用于寻找微血管网络的中心线。这些基于阈值设定的方法都是把图片视作光学信号的简单集合而没有考虑图案本身的纤维状特点，因此这些方法都无法得到纤维网络结构中单根纤维的拓扑形貌信息。

海赛矩阵（Hessian matrix）被用来估计每个像素点所在处的纤维方向以及轮廓的中心线，以用于神经突触骨架结构的标记和冠状动脉搏动波的跟踪。

另外，一种基于约束优化的方法被用于跟踪严重噪声情况下的曲线纤维结构。基于展开动态轮廓法（Stretching Open Active Contours , SOACs）全自动方法能从显微图像中提取肌动纤维网络。然而，基于追踪算法或SOACs算法的图像处理方法均只能提取单根纤维，而不能获取整个纤维网络的形貌信息。

此外也有很多其他处理曲线纤维图像的方法，比如基于蛇模型的算法获取光滑骨架结构、三维线结构增强过滤器从其他结构中提取线型结构等，然而这些方法仍然无法获得纤维网络结构中单根纤维的形貌特征。

发明内容

本发明所要解决的技术问题是提供一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，能够得到亚细胞尺度的曲线纤维网络结构形貌特征，以用于定量分析亚细胞尺度的曲线纤维网络结构。

为解决上述技术问题，本发明所采用的技术方案是：一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，包括以下步骤：

一、纤维网络的骨骼化，得到纤维网络骨架图案；

二、纤维网络骨架图案的拓扑特征分类，识别出纤维网络中的单纤维及纤维交叉点；

三、删除短纤维，保留主干纤维的网络骨架图案；

四、纤维片段重组，识别纤维网络中的交叉、分叉及重叠关系，并重组纤维网络图；

通过上述步骤得到亚细胞尺度的曲线纤维网络结构形貌特征。

纤维网络的骨骼化即为用一个像素宽度的曲线网络来取代纤维网络。

优选的方案中，将二值化后的纤维网络图案进行一次骨骼化，然后进行一次膨胀处理，继续进行一次骨骼化操作，最终得到纤维网络的骨骼图案。

优选的方案中，步骤二中，依据拓扑数N，即骨架像素周围8邻域中出现骨架像素的个数，将信号像素分成不同的拓扑特征结构类别，即孤立点、曲线端点、曲线内点和曲线交叉点。

优选的方案中，步骤三中，删除仅有1个或0个曲线内点的纤维。

进一步优选的方案中，一、拓扑特征结构分类分析，获得交叉点及端点数据；

二、删除交叉点及端点，得到仅含内部点的离散纤维片段图；

三、删除孤立像素点，恢复交叉点和端点，再次删除孤立像素点，得到消除了具有1个内部点的短纤维的网络骨架图案；

四、拓扑特征结构分类分析，获得交叉点数据；

五、删除交叉点，获得含内部点及端点的离散纤维片段图；

六、删除孤立像素点，恢复交叉点，得到消除了具有0个内部点的短纤维的网络骨架图案；

七、将当前图案与本次迭代初始的网络骨架图案进行比对计算，当残差≠0，则重复上述步骤，当残差＝0，则结束，获得短纤维删除后的网络骨架图案。 

优选的方案中，步骤四分为以下步骤:

一、建骨架片段端点数据库；

二、建立端点间距矩阵；

三、四端点紧邻体连接模式识别及交叉纤维片段重组；

四、三端点偶联体连接模式识别及重合纤维片段重组；

五、三端点孤立体连接模式识别及交叉纤维片段重组；

通过上述步骤重组纤维网络图。

优选的方案中，步骤一之前还有以下预处理步骤：

一、原始图像增强；

二、采用Tubeness过滤器识别纤维结构信息；

三、去除微弱纤维信号。

进一步优选的方案中，图像增强采用图像平方增强、谐波均值滤波、二维中值滤波和图像非线性调整方法中的一种或多种组合。

进一步优选的方案，预处理步骤三中，对经过Tubeness过滤器处理的图案进行联通域识别及性质分析，获得每个联通域的几何尺寸和平均光强信息，每个联通域对应着若干根相互连接的纤维组成的网络或是单根孤立的纤维；

针对图案中每个联通域定义纤维系数，如果纤维系数大于阈值就判断此联通域对应的纤维或纤维网络信号为焦平面内实际存在的纤维，予以保留；

如果纤维系数小于阈值就判断为焦外纤维在焦平面的投影，予以删除。

纤维系数阈值为0.005±0.005。

本发明提供的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，以下简称IFS算法，可用于定量分析亚细胞尺度的曲线纤维网络结构，如由肌动蛋白构成的应力纤维骨架网络、微管蛋白骨架网络、微丝蛋白骨架网络、核骨架网络等，且能够识别纤维间的连接关系，并提取单根纤维的长度、方向、位置等拓扑信息，呈现整个网络中所有纤维不同拓扑性质的分布特征。与同类算法相比，IFS算法的最大优点在于它在完成了对曲线网络中纤维结构的识别和分割后，还能够不依赖人为的干预和判断完成对显微图像中曲线网络结构的定量拓扑分析。其他针对曲线网络结构显微图片的图像处理方法都仅是完成了对曲线网络结构的识别和分割，进一步针对此网络结构的定量分析依然需要借助人为的判断及干预。经过验证，本发明的IFS算法具有很强的抗噪能力，可以被用于分析各类显微图像，如普通明场显微镜、激光共聚焦显微镜Confocal，全内反射纤维镜TIRF等。因此能很大程度地提高获取大量生物实验图片背后隐藏着的信息和规律的能力。

附图说明

下面结合附图和实施例对本发明作进一步说明：

图1为本发明方法的整体流程图。

图2为本发明方法的短纤维删除迭代算法流程图。

图3为本发明方法中初始端点间距矩阵以及端点间距矩阵的数据结构示意图。

图4为本发明方法中纤维间三种连接模式及连接处像素点示意图。

图5为采用本发明方法分析真是细胞应力纤维骨架网络的过程图。

图6为采用本发明方法处理后得到的纤维网络拓扑信息。

图7为本发明中纤维系数的参数示意图。

具体实施方式

1.本发明的流程如图1所示。 

具体为：

2.1增强纤维状结构信息的图像前处理算法 

2. 1.1 图像增强（Enhance image） 

采用光学显微镜获得的细胞显微图像由于Abby极限、背景噪声的存在而往往具有较低的清晰度和较差的分辨率。直接在原始显微图像基础上提取和分析显微结构信息的效果比较差，一般需要先做图像增强处理，使得结构性信号得到增强，离散噪声信号得到抑制。传统图像增强方法通常是依赖于某种截断阈值，如灰度阈值、灰度梯度阈值等，高于此截断阈值的信号被加强，而低于此阈值的信号被视为噪声而消除。本发明方法中，尽量避免使用类似的截断阈值，因为设置任何截断阈值都可能将破坏纤维结构的完整性，即将一根纤维人为地打断成若干段相互分离的子纤维结构。本发明将多种图像增强方法组合，实现加强结构性信号减弱离散噪声信号的效果，本质上是通过一些非线性映射和滤波的方法增加了图像的对比度。具体而言，针对显微图像采用了图像平方增强、谐波均值滤波（对抑制盐噪声效果突出）、二维中值滤波、图像非线性调整等几种现有技术中的图像增强方法进行处理。本发明的图像处理操作都是在Matlab中编程实现。

2.1.2 采用Tubeness过滤器识别（discriminate）纤维结构信息 

Tubeness是开源生物图像处理软件ImageJ的一款用于识别生物图像中的曲线结构并恢复其宽度的图像过滤器插件。此过滤器的图像分割算法是基于海森矩阵（Hessian matrix）特征值分析设计的，有关此过滤器的细节可以参见Yoshinobu Sato的文章。经过图像增强处理的显微图像已经获得了较高的图像对比度，但是图像中除了网状曲线连接结构外仍然存在大量弥散的噪声信号和呈团块状的结构信号，这将对后续识别和分析网状曲线连接结构带来很大的干扰。因此，采用Tubeness过滤器对已经经过图像增强的图像进行处理，从各种其它结构信号中只识别和提取出曲线纤维结构，并保持了对应曲线结构的光强强度信息。

2.1.3 去除微弱纤维信号 

经过Tubeness过滤器处理后得到了光强强度各异的一系列曲线纤维结构，有的相互连接构成局部网络，有的孤立地存在。基于显微成像理论，把那些面积微小、光强微弱的纤维或纤维网络信号视为焦平面之外的纤维或纤维网络在焦平面的投影，而不是焦平面处实际存在的纤维。在程序实现中，需要对经过Tubeness处理的图片进行联通域识别及性质分析，获得每个联通域的几何尺寸和平均光强等信息，每个联通域对应着若干根相互连接的纤维组成的网络或是单根孤立的纤维。我们针对图像中每个联通域定义了纤维系数 ，如果纤维系数F大于阈值F0就判断此联通域对应的纤维或纤维网络信号为焦平面内实际存在的纤维，反之，就判断为焦外纤维在焦平面的投影，予以删除。 

纤维系数 的定义见式（1）。  

(1)中I和A分别为当前联通域的平均光强和面积，Imax和Amax分别表示图片内所有联通域中最大的平均光强值和最大的面积。由定义可知纤维系数F的为0至1范围之间的数值，纤维系数值F越接近1意味着此联通域的光强越强、面积越大；反之，纤维系数值F越接近0意味着此联通域的光强越小、面积越小。本例中，将缺省的纤维系数阈值F0设为0.005，在缺省设置下联通域若同时满足平均光强小于最大平均光强的7%，以及联通域面积小于最大面积的7%两个条件，此联通域对应的纤维或纤维网络将被视为焦外纤维在焦平面的投影而被删除，当然以上两个条件只是缺省设置下联通域被删除的一种充分而非必要条件。

通过基于纤维系数阈值的判断，可以成功地去除一些面积微小、光强微弱的纤维信号，而不会破坏原有纤维网络结构的联通性，即不会将一根纤维或一个局部纤维网络人为地打断成若干段相互分离的子纤维结构。这符合前述的对具有纤维网络结构特点的图像进行处理的要求。

2.2 IFS图像分析处理核心算法 

2.2.1 纤维网络的骨骼化（Skeletonized） 

为了便于分析纤维网络的拓扑特征，将纤维网络用一个像素宽度的曲线网络来取代。此曲线网路必须保持和原纤维网络具有相同的连通性，此曲线网络中的曲线必须保持和其对应的纤维具有相同的长度和取向。骨骼化算法和细化算法都能确保细化图案的同时，不破坏图案的连通性以及不改变图案的尺寸。与细化算法相比，骨骼化算法处理后的骨架图案不容易出现细小的毛刺。因此，本例中采用骨骼化算法处理纤维网络图案获得纤维网络的骨骼图案。在具体的程序实现中，对二值化后的纤维网络图案进行一次骨骼化，紧接着进行了一次膨胀处理，目的是为了将由于数值图像处理误差而人为造成断裂（通常断裂长度不超过2个像素）的纤维重新连接起来，然后继续进行一次骨骼化操作，最终得到纤维网络的骨骼图案。

2.2.2 纤维网络骨架图案的拓扑特征分类算法

通过纤维网络骨架图案的拓扑特征分类，进一步得到纤维网络中单根纤维的长度、取向、位置等信息以及不同纤维之间的分叉、交叉、重叠等关系。为了将纤维网络骨架图案的拓扑特征进行分类，需要采用拓扑数N，拓扑数N被定义为骨架像素周围8邻域中出现骨架像素的个数。因此N是范围在1到8之间的整数。依据拓扑数N的不同就可以将信号像素分成不同的拓扑特征结构类别 (Topological Classification)。在骨骼化图案中，包括三种拓扑特征结构类别，分别是曲线端点 (End arc)、曲线内点 (Inner arc)、曲线交叉点 (Arc junctions)。拓扑数N和拓扑特征结构类别的对应关系如下： 

Class Ⅰ：N=0，孤立点 (Isolated point)；

Class Ⅱ：N=1，曲线端点 (End point of arc)；

Class Ⅲ：N>=2，曲线内点 (Inner point of arc) ；

Class IV：T>=3，曲线交叉点 (Junction point of arc) ；

针对纤维网络的骨架图案，能够计算所有骨架像素的拓扑数，并能按拓扑数大小将对应的骨架像素分成四种拓扑特征结构类型。将具有相同拓扑特征的骨架像素用同一种的颜色标记出来后，就可以直观观察骨架网络的连接关系。 

2.2.3短纤维删除算法（flow chat）

 经过拓扑特征分析，就能观察纤维网络骨架图案连接关系。发明人发现大部分曲线交叉点 (Arc junctions) 的出现是由非常短的分支结构 (仅含有1个inner point) 造成的。通常这些图案的研究者是期望研究那些主要的长纤维结构之间的相互交叉情况，而这些短的分支结构引起的交叉点数量很大，会模糊研究焦点。另外，这些短的分支结构并不全是物理存在的短小纤维，而在很大比例上是数字图像处理算法对纤维形状过于敏感引入的毛刺。因此，需要设计新的图像处理算法来删除这些短分支结构，消除它们对纤维网络骨架结构拓扑特性的影响。在本发明将该算法称为短纤维删除算法。 

短纤维删除算法的核心思路是基于拓扑特征分析结果删除骨架网络的交叉点从而获得了判断和删除离散的短纤维片段的可能。此算法中，需要被删除的短纤维被定义为仅有1个或0个曲线内点的纤维，拥有0个曲线内点的纤维意味着只由一个端点组成。整个短纤维删除算法是一个迭代过程，首先被删除的是初始的网络骨架图案中拥有1个曲线内点的短纤维，然后再删除拥有0个曲线内点的短纤维。每次迭代过程中，在删除短纤维的同时网络骨架图案的拓扑特性也在发生变化，即可能会有新的短纤维被再次识别出来，例如：由多个交叉点组成的结构体的一些突起会被再次识别成新的短纤维而被逐渐地删除。这就是为什么不能一次就删除所有短纤维，而需要一个迭代过程才能完全删除所有潜在可能的短纤维的原因。当某次迭代结束后得到的网络骨架图案与此迭代初始骨架网络完全相同时，认为已经完全删除了骨架网络中的短纤维，迭代过程结束。短纤维删除迭代算法中的更细致的环节请参见图2。 

短纤维删除算法是一个迭代过程，在一次迭代中首先需要分析网络骨架图案的拓扑特性，进而删除纤维结构中的交叉点和端点，只留下由仅由内部点组成的离散纤维片段图，然后依次通过删除图像中的孤立点、恢复交叉点和端点、再次删除图像中的孤立点这三步操作就删除了具有1个内部点的短纤维；随后继续分析网络骨架图案的拓扑特性，删除当前图像中的交叉点后获得由内部点和端点组成的离散纤维片段图，然后依次通过删除图像中的孤立点、恢复交叉点这两部操作就删除了具有0个内部点组成的短纤维。此后需要计算当前网络骨架图案与本次迭代的初始网络骨架图案的残差，当此残差为零时，就认为已经完全删除了骨架网络中的短纤维，迭代过程至此结束。 

2.2.4 纤维片段重组算法（Regroup the segments） 

本发明的目的是获取纤维网络的拓扑特征信息，包括纤维网络中单根纤维的长度、取向、位置等信息。若要获得上述纤维的拓扑特征信息，就需要设计算法来从纤维网络中辨识出每一根独立的纤维，本例中将该算法称为纤维片段重组算法。纤维片段重组算法实现了对不同纤维之间分叉、交叉、重叠等关系的识别，并在此基础上将在连接处被打断的纤维片段按照正确的连接关系重新组合成纤维网络。经过纤维片段重组算法处理后，就能得到纤维网络中所有纤维各自的长度、取向、位置等信息，也能得到不同纤维之间的分叉、交叉、重叠等拓扑学关系。

纤维片段重组算法的基本思路是将纤维网络骨架结构在连接处打断从而形成相互孤立的纤维片段，根据这些纤维片段端点之间的相互距离设定特定的判据，识别出不同的连接模式，根据不同连接模式的固有特点将纤维片段重新组合形成连接关系明确的纤维网络。骨架纤维片段的联通域性质被提取出来并存放在骨架纤维片段联通域性质矩阵中。

在本发明中，定义了纤维片段之间的三种连接模式，分别是代表分叉关系的“三端点孤立体”连接模式，代表交叉关系的“四端点紧邻体”连接模式和代表重叠关系的“三端点偶连体”连接模式。在同一种连接模式下，纤维片段能通过相同的规则确定相互之间的连接关系。 

建骨架片段端点数据库。 

纤维片段重组算法的第一步就是将纤维骨架网络在连接处打断，形成纤维骨架片段，即删除拓扑数（由2.2.2节定义）大于等于3对应的骨架像素。然后将纤维骨架片段图的联通域进行编号，并重新计算每个骨架像素的拓扑数，下一步将第m个纤维端点（拓扑数为1）所处的联通域编号、x坐标、y坐标提取出来，依次存入n×3数组的第m行，n为所有纤维端点的个数。最后按照数组第一列（即联通域编号）进行升序排列，调整数组的行号。本例中将这个包含所有纤维端点所处的联通域编号及坐标信息的矩阵称为骨架片段端点数据库。参见图3。 

建立端点间距矩阵。 

纤维之间的不同连接模式的识别判据是基于纤维片段端点之间的距离，因此需要建立端点间距矩阵。首先建立n╳n的初始端点间距矩阵D ，计算任意两个端点(i, j)之间的间距D(i, j)，由矩阵的定义可知D是对称矩阵，D的数据格式如图3.A所示。端点间距矩阵的建立旨在考察分属不同联通域间的端点之间的距离，因此将对称轴上的间距值D(i, i)以及属于同一个联通域的不同端点间的距离设为无穷大。然后对初始端点间距矩阵D的每一列进行升序排列形成端点间距矩阵DS，DS矩阵可以帮助我们很容易地找到与任意端点m的距离在某设定值之内的若干个端点的信息。图3.B中，DS(i, m)的值为D(t, m), 这表示第i个最靠近端点m的端点的编号为t，它们之间的距离为D(t, m)。由于初始端点间距矩阵对称轴上的间距值D(i, i)以及属于同一个联通域的不同端点间的距离被设为无穷大，因此它们在DS中总是被排到最后几行。

图3.A为初始端点间距矩阵D，其行号和列号和端点序号重合，D(i, j)表示端点i和端点j之间的距离，对称轴上的间距值D(i, i)以及属于同一个联通域的不同端点间的距离设为无穷大。

图3.B为对D的每一列进行升序排列后的矩阵DS的数据结构，称为端点间距矩阵。DS矩阵的列号与端点矩阵序号重合，第m列的值表示纤维骨架片段图中所有的n个端点与端点m的距离值的升序排列。其中DS(i, m)的值为D(t, m), 这表示第i个最靠近端点m的端点的编号为t，它们之间的距离为D(t, m)。 

删除纤维网络骨架连接处的像素后，网络结构被打断成了大量孤立的纤维骨架片段，这些纤维片段在未被离散前具有若干种连接模式。其中最典型的连接模式是如下三种，图4(a) 反应纤维分叉的三端点孤立体连接模式；图4(b) 反应纤维交叉的四端点紧邻体连接模式；图4(c) 反应纤维重叠的三端点偶连体连接模式。图4(a)、(b)、(c)中被删除的连接处像素子集如图4(d)所示。图4(d) 中第一排像素体具有三个端点和一个交叉点，他们代表了图4(a) 中被删除的三端点孤立体和图4(c) 中被删除的三端点偶连体中的两个三端点紧邻体。图4(d) 中第二排像素体均具有4个端点，而两个交叉点之间分别错开0至4个像素，它们都是图4(b) 中被删除的四端点紧邻体可能的形式。图4(d)中的有效宽度W是指连接到对应连接体的纤维片段端点间的最大距离，这个有效宽度被用来设定识别不同的连接模式的判定阈值。图4(d) 中第三排像素体具有4个端点，而两个交叉点之间的公共纤维片段超过2个像素，这是图4(c) 中被删除的三端点偶连体的具体形式。 

四端点紧邻体连接模式识别及交叉纤维片段重组。

纤维网络骨架在连接处被打断后形成了大量孤立的纤维骨架片段，这些纤维片段在未被离散前具有若干种连接模式。其中最典型的连接模式之一是代表纤维交叉的四端点紧邻体连接模式，如图4. b所示，其中四端点紧邻体如图4(d)中后4至8号小图所示。

形成四端点紧邻体的判据为：

1) 四个端点两两之间的距离均小于阈值Gap_4p；

2) 不同的四端点紧邻体相互间不重复；

3) 四个端点分别属于不同联通域。 

在程序实现中，将判定四端点紧邻体的阈值Gap_4p设置为8个像素，这意味这把两根纤维直接交叉，如图4(d) 4号小图所示；以及错开1至4个像素交叉，如图如图4(d) 5~7号小图所示的5种纤维间交叉形式都认为属于由四端点紧邻体介导的纤维交叉模式。阈值Gap_4p被设置为这些不同的连接体的最大的有效宽度，即为8个像素。当两根纤维错开大于4个像素交叉时，两个交叉点之间具有连续的超过两个拓扑数为2的内部点，2.3小节中的短纤维删除算法会保留这些内部点，这种交叉模式会被自动识别成符合后文所述的三端点偶连体介导的纤维重复连接模式。图4 (d)中4至8号小图是四端点紧邻体的典型连接体形式，但不是全部，仍然还有一些连接体形式的变体没有用图呈现出来，但我们设定的判据会甄选出所有的四端点紧邻体连接形式。因此，只要满足条件：端点间距矩阵DS(3, m)小于等于阈值 8个像素, 就能找到与端点m距离小于8个像素的另外三个端点。对满足以上条件的每一个四端点子集，进一步计算每两个端点间的距离，排除有间距大于阈值8个像素的情况。之后，比较不同的四端点子集，删除重复的子集。最后，通过骨架片段端点数据库将四端点子集中的端点序号更新为与之对应的联通域序号，删除联通域序号中有重复的四端点子集。最后剩下的四端点子集就是满足判据的四端点紧邻体了。将由四端点紧邻体介导的纤维交叉中方向相近的纤维片段连接起来，重组成新的纤维，并更新骨架纤维片段联通域性质矩阵。

以图4(b)为例，方向相近的纤维片段S1和S4，以及S2和S3分别被连接起来构成两条新的纤维。可以将两条纤维分别用红色和蓝色表示，交叉连接部分显示为紫色。 

三端点偶联体连接模式识别及重合纤维片段重组。

纤维骨架片段在未被离散前的另一种典型的连接模式是代表纤维重叠的三端点偶联体连接模式，如图4（c）所示，其中三端点偶连体中的三端点紧邻体如图4(d)中1、2和3号小图所示。 

形成三端点紧邻体的判据为：

1) 三个端点两两之间的距离均小于Gap_3p个像素；

2) 不同的三端点紧邻体相互间不重复；

3) 三个端点分别属于不同联通域。

 三端点偶连体由两组三端点紧邻体构成，这两组三端点紧邻体中各有一个端点属于某一条共同的纤维片段，如图4 (c)中的纤维片段S1。形成三端点偶连体的判据为：

1) 两组三端点紧邻体各有一个端点属于同一条纤维片段；

2) 相同的三端点紧邻体不能在不同的三端点偶联通中重复出现；

3) 三端点偶联体中的两个三端点紧邻体都不能是任意一个四端点紧邻体的子集。 

在程序实现中，将判定三端点紧邻体的阈值Gap_3p设置为4个像素，这是因为对应于三端点紧邻体的两种连接体和的有效宽度不超过4个像素。只要满足条件：端点间距矩阵DS(2, m)小于等于阈值4个像素, 就能找到与端点m距离小于四个像素的另外两个端点。对满足以上条件的每一个三端点子集，进一步计算每两个端点间的距离，排除有间距大于阈值4个像素的情况。之后，比较不同的三端点子集，删除重复的子集。最后，通过骨架片段端点数据库将三端点子集中的端点序号更新为与之对应的联通域序号，继而删除联通域序号中有重复的三端点子集。最后剩下的三端点子集就是满足判据的三端点紧邻体了。搜索三端点偶联的第一步是找出成对的三端点紧邻体，他们中分别具有一个端点属于同一条纤维片段。然后遍历搜索所有三端点偶连体中的紧邻体，排除其中有重复的三端点紧邻体出现的情况。最后再次遍历所有三端点偶连体中的紧邻体，若其中的三端点紧邻体属于某四端点紧邻体子集，则将其删除。最后剩下的三端点紧邻体对就是符合判定标准的三端点偶连体了。将由三端点偶连体介导的纤维交叉中除共享纤维片段之外的方向相近的纤维片段连接起来，重组成新的纤维，并更新纤维骨架片段联通域性质矩阵。以图4.(C)为例，S1为共享纤维，方向相近的纤维片段S2和S5，以及S3和S4被分别与共享纤维S1连接起来构成两条新的纤维。这两条纤维可以分别用红色和蓝色表示，重合部分显示为紫色。 

三端点孤立体连接模式识别及分叉纤维片段重组。

纤维骨架片段在未被离散前的第三种典型的连接模式是代表纤维分叉的三端点孤立体连接模式，如图4（a）所示，其中三端点孤立体与三端点紧邻体形状类似，如图4(d)中1、2和3号小图所示。 

三端点孤立体是特殊的三端点紧邻体，与三端点紧邻体相比，它限定了与其最近的纤维片段端点和其中任意一个端点的距离必须大于阈值Gap_3p_Isolate，并且已经构成三端点偶连体的三端点紧邻体不再构成三端点孤立体。

因此，形成三端点孤立体的判据为：

1) 三个端点两两之间的距离均小于4个像素；

2) 此三端点与距离其最近的第4个端点间的距离均大于4个像素；

3) 三个端点分别属于不同联通域；

4) 此三端点不属于三端点偶连体的子集。 

在程序实现中，将判定三端点孤立体的阈值Gap_3p_Isolate也设置为4个像素，这是因为对应于三端点孤立体的两种连接体和的有效宽度不均超过4个像素。只要满足条件：端点间距矩阵DS(2, m)小于等于阈值4个像素, 并且DS(3, m)大于阈值4个像素，就能找到满足判据2)的与端点m距离小于四个像素的另外两个端点。对满足以上条件的每一个三端点子集，进一步计算每两个端点间的距离，排除有间距大于阈值4个像素的情况。之后，比较不同的三端点子集，删除重复的子集。然后，通过骨架片段端点数据库将三端点子集中的端点序号更新为与之对应的联通域序号，进一步删除联通域序号中有重复的三端点子集。最后，遍历已近得到的三端点偶连体中的紧邻体，删除与之相重复的三端点孤立体。剩下的三端点子集就是满足判据的三端点孤立体了。将由三端点孤立体介导的纤维分叉中方向相近的纤维片段连接起来，重组成新的纤维，并更新纤维骨架片段联通域性质矩阵。以图4 (a)为例，方向相近的纤维片段S1和S2被视为一根纤维而连接起来，另一根纤维S3被认为是从主干纤维上分叉出来的枝干纤维。主干纤维用红色表示,枝干纤维用蓝色表示，分叉部位显示为紫色。

实施例1：

分析真实细胞应力纤维网络图像 

将本发明的IFS算法应用于分析真实细胞应力纤维网络图像。成骨细胞(Osteoblast cell line：MC3T3)的肌动纤维(actin filament)骨架网络被XX荧光标记，在激光共聚焦显微镜下观测获得显微图像，见图5(a)。此图像中心为细胞中心，图像四周已近接近细胞边缘，图像中每个像素的尺寸为0.154微米。

用IFS程序分析真实细胞应力纤维骨架网络的过程主要过程包括：

图5(b)为图像预处理；

图5(c)为删除微弱离面纤维信号及骨骼化；

图5(d)为骨架像素拓扑分类；

图5(e)为删除短纤维；

图5(f)为拓扑连接关系识别及纤维网络重组。

经过IFS程序分析，得到分析的真实细胞应力纤维骨架网络中有101条纤维，它们之间形成了76个交叉点，这些交叉将整个骨架网络分成205根离散的纤维片段。构成骨架网络的101根纤维的方向主要集中在45°的方向上，更具体的方向分布如图6(a)所示。这些纤维的长度比较集中在40个像素以下，而最长的纤维也可以达到150个像素以上。纤维长度的分布比较符合Logistic分布的特点，本例中用Logistic分布函数：

；

拟合了纤维长度的分布，如图6(b)所示。也可以得到骨架网络中76个交叉点离图像中心的距离。通过距离的分布图图6(c)，可以发现交叉点很少分布在接近细胞中心的细胞核区域，这个中心区域的半径大约为20个像素；交叉点在远离细胞核区域有较多的分布，且分布相对比较均匀。纤维交叉点距离图像中心距离的分布比较符合LogNormal分布的特点，用Lognormal分布函数：

；

拟合其分布，如图6(c)所示。

在IFS算法中，本例定义了一个阈值参数，即纤维系数阈值F0，若纤维系数F小于此阈值F0，此纤维就会被认为是焦平面以外的纤维在焦平面上的投影而被删除。若纤维系数阈值F0值被设置得越小，就有越少的纤维被删除，反之则有越多的纤维被删除。在1e-4到0.1范围之间，本例选取了7个值作为不同的纤维系数，分别对同一张细胞应力纤维网络骨架图进行了IFS分析, 结果如图7所示。结果显示，当F0值从0.1减小到0.005时，经IFS分析所得的网络结构交叉数和纤维片段数都在不断增加；当F0值从0.005继续减小到1e-4时，经IFS分析所得的交叉数几乎不发生变化，而纤维片段数还在继续增加。可以认为0.005是一个拐点，小于此拐点的F0值设置不会改变网络的连接关系，而只会引入一些孤立的纤维片段。因此，本例就把这个拐点值0.005设置为纤维系数阈值。用来作为判断某条纤维是否属于焦平面以外的纤维在焦平面上的投影的判据。当然0.005的纤维系数阈值可以在±0.005的范围内调整。

发明人认为纤维系数阈值F0只与显微成像系统和观测对象的几何尺寸有关，因此分析同一种纤维网络结构的同批次显微图像只需要设置一个纤维系数阈值。在IFS程序中，F0的缺省值被设置为0.005，若要分析一批新的显微图像，可以按照本小节所述的方法自行律定合适的纤维系数阈值。  

图5 用IFS程序分析真实细胞应力纤维骨架网络的过程。图5(a)经XX染色的成骨细胞应力纤维骨架网络；图5(b)经过IFS方法中的图像预处理后获得的应力纤维网络灰度图；图5(c)经过离面微弱信号去除算法及骨骼化算法后得到的应力纤维网络骨架图；图5(d)应用短纤维删除算法前应力纤维网络骨架的拓扑分类图，可见图中有大量的连接点和短小的毛刺纤维。图5(e) 应用短纤维删除算法后应力纤维骨架的拓扑分类图，相比于图5(d)图用红色表示的连接点及短小的毛刺纤维有明显的减少，只剩下由较长的纤维构成的网络连接图。图5(f) 经过纤维重组算法处理后得到的应力纤维骨架网络图。  

图6 经IFS算法处理后得到的纤维网络拓扑信息。用IFS算法处理图5(a)图所示的应力纤维骨架网络图，得到应力纤维的方向分布图6(a)，应力纤维的长度分布图6(b)，及应力纤维间的连接点距离图片中心点距离的分布图6(c)。由图图6(a)可知，应力纤维主要分布在45°范围内。应力纤维的长度分布的特点比较接近Logistic分布，构成应力纤维骨架网络的纤维以较短的纤维为主(0~40 pixels)，长度大于40 pixel的纤维数量较少，纤维的长度可以达到150 pixel以上。应力纤维间的连接点距离图片中心点距离的分布特点比较接近Lognormal分布，连接点在细胞的中心区域出现概率较低，在远离细胞中心的区域出现概率比较平稳且均比较高。 

图7是IFS算法中纤维系数的参数研究。利用IFS算法处理应力纤维骨架网络图fig.5(a)的过程中，随着纤维系数F从0.1开始被设置得越来越小，连接数、纤维片段数和纤维数都将越来越大；直到F=0.005后，连接数基本保持不变，而纤维片段数和纤维数继续增大，直到趋于某个最大值。 

Claims (11)

1.一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是包括以下步骤：

一、纤维网络的骨骼化，得到纤维网络骨架图案；

二、纤维网络骨架图案的拓扑特征分类，识别出纤维网络中的单纤维及纤维交叉点；

三、删除短纤维，保留主干纤维的网络骨架图案；

四、纤维片段重组，识别纤维网络中的交叉、分叉及重叠关系，并重组纤维网络图；

通过上述步骤得到亚细胞尺度的曲线纤维网络结构形貌特征。

2.根据权利要求1所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是：纤维网络的骨骼化即为用一个像素宽度的曲线网络来取代纤维网络。

3.根据权利要求1或2所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是：将二值化后的纤维网络图案进行一次骨骼化，然后进行一次膨胀处理，继续进行一次骨骼化操作，最终得到纤维网络的骨骼图案。

4.根据权利要求1所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是：步骤二中，依据拓扑数N，即骨架像素周围8邻域中出现骨架像素的个数，将信号像素分成不同的拓扑特征结构类别，即孤立点、曲线端点、曲线内点和曲线交叉点。

5.根据权利要求1所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是：步骤三中，删除仅有1个或0个曲线内点的纤维。

6.根据权利要求5所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是：

一、拓扑特征结构分类分析，获得交叉点及端点数据；

二、删除交叉点及端点，得到仅含内部点的离散纤维片段图；

三、删除孤立像素点，恢复交叉点和端点，再次删除孤立像素点，得到消除了具有1个内部点的短纤维的网络骨架图案；

四、拓扑特征结构分类分析，获得交叉点数据；

五、删除交叉点，获得含内部点及端点的离散纤维片段图；

六、删除孤立像素点，恢复交叉点，得到消除了具有0个内部点的短纤维的网络骨架图案；

七、将当前图案与本次迭代初始的网络骨架图案进行比对计算，当残差≠0，则重复上述步骤，当残差＝0，则结束，获得短纤维删除后的网络骨架图案。

7.根据权利要求1所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是步骤四分为以下步骤:

一、建骨架片段端点数据库；

二、建立端点间距矩阵；

三、四端点紧邻体连接模式识别及交叉纤维片段重组；

四、三端点偶联体连接模式识别及重合纤维片段重组；

五、三端点孤立体连接模式识别及交叉纤维片段重组；

通过上述步骤重组纤维网络图。

8.根据权利要求1所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是在步骤一之前还有以下预处理步骤：

一、原始图像增强；

二、采用Tubeness过滤器识别纤维结构信息；

三、去除微弱纤维信号。

9.根据权利要求8所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是：图像增强采用图像平方增强、谐波均值滤波、二维中值滤波和图像非线性调整方法中的一种或多种组合。

10.根据权利要求8所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是：预处理步骤三中，对经过Tubeness过滤器处理的图案进行联通域识别及性质分析，获得每个联通域的几何尺寸和平均光强信息，每个联通域对应着若干根相互连接的纤维组成的网络或是单根孤立的纤维；

针对图案中每个联通域定义纤维系数，如果纤维系数大于阈值就判断此联通域对应的纤维或纤维网络信号为焦平面内实际存在的纤维，予以保留；

如果纤维系数小于阈值就判断为焦外纤维在焦平面的投影，予以删除。

11.根据权利要求8所述的一种基于数字图像处理的曲线纤维网络结构形貌特征测量方法，其特征是：纤维系数阈值为0.005±0.005。

Images