CN101599176A - 一种对管状结构内层分割的方法 - Google Patents

Abstract

本发明公开了一种对管状结构内层分割的方法，该方法首先提取管状结构三维序列数据的纵向截平面，并在纵向截平面提取管状结构的内层轮廓线，然后把纵向截平面的轮廓线变换到横向截平面，接着把横向截平面的轮廓点样条插值得到样条轮廓，利用小波变换产生横向截平面的多尺度图像，并把样条轮廓变换到不同尺度作为活动轮廓的初始轮廓，然后进化轮廓，最后结合不同尺度活动轮廓的分割结果，得到管状结构的最终内层轮廓。应用本发明的技术方案减少了初始轮廓的时间。为提高分割精度，该方法在多尺度上使用活动轮廓进化轮廓。

Description

一种对管状结构内层分割的方法

技术领域

本发明涉及计算机图像处理领域，特别涉及一种基于模糊聚类和活动轮廓的结合纵向和横向截平面的管状结构内层分割的方法。

背景技术

图像分割是指将图像中若干感兴趣的区域提取出来，感兴趣的区域通常称为前景，其它区域称为背景。在分割的基础上才可能对目标进行定量分析或者识别，从而使更高层的分析和理解成为可能。

九十年代以来，随着技术的发展，可以获得更高分辨率的图像，基于形变模型的方法也开始大量应用。形变模型是一种寻找图像中物体边界的方法。给定初始轮廓，然后通过图像信息及先验知识定义外力，推动轮廓到达物体的边界。形变模型分为两大类：参数形变模型和几何形变模型。

参数化形变模型最早可以追溯到Kass等人提出的活动轮廓模型。分割过程就是模型在外力和内力的作用下向物体边缘靠近。内力控制曲线的形状，外力推动曲线运动。然而传统的活动轮廓模型对初始轮廓的敏感性及有限的捕捉范围。对此，许多学者提出了多种外力的定义。M.Kass和D.Terzopoulos使用高斯函数对图像进行平滑，虽然扩大了捕捉范围，但是随着σ的增加，可能会使边缘平滑的太厉害，从而找不到真实的边缘。L.D.Cohen等提出了balloon模型，在外力中增加了压力(pressureforce)项，该模型改善了活动轮廓对初始轮廓的敏感性，并且能够跨越伪边缘点，但是需要控制压力的大小及方向。而且该模型无法收敛到深度凹陷部分。L.D.Cohen和I.Cohen提出了距离力可以扩大捕捉范围，但是距离力使模型上的点靠近最近的边界点，这使得轮廓无法收敛到深度凹陷处。

1998年Xu等提出了梯度矢量流来解决不能收敛到深度凹陷处的问题，他运用向量扩散公式将梯度扩散到远离边缘的区域中，于是产生了另外一个力场，被称为梯度矢量流(gradient vector flow，GVF)场。GVF改善了边缘凹陷问题。虽然GVF解决了活动轮廓捕捉范围小的问题，但是如果初始轮廓位置不正确的话，仍然有可能导致活动轮廓收敛到错位的位置。这导致了基于活动轮廓的方法都需要人工设定初始轮廓，然而如果图像的数据量很大，人工设定初始轮廓即费力又费力，重复性又差，因此必须使用自动化的方法进行分割。

发明内容

本发明的目的是提供一种自动的图像分割方法，能够对管状结构内层进行准确、自动分割的方法。

为实现上述目的，本发明的一种对管状结构内层分割的方法，结合纵向和横向截平面的管状结构内层自动分割的步骤包括：

步骤a：对管状结构三维序列数据进行各向异性扩散滤波，用于除去随机噪声；

步骤b：从管状结构三维序列数据中提取纵向截平面，并在纵向截平面上利用模糊聚类提取管状结构内层的轮廓线；

步骤c：把纵向截平面提取的轮廓线变换到横向截平面，每一个横向截平面得到对应的轮廓点；

步骤d：对横向截平面上的轮廓点进行样条插值，得到一条封闭的样条轮廓；

步骤e：利用小波变换产生横向截平面的多尺度图像，用来在多尺度图像中使用活动轮廓，多尺度图像包含三个尺度，即第一尺度、第二尺度和第三尺度；

步骤f：把样条轮廓变换到不同尺度，用来作为多尺度图像中活动轮廓的初始轮廓；

步骤g：计算横向截平面的多尺度图像的梯度矢量流场；

步骤h：在横向截平面的多尺度图像中利用活动轮廓进化轮廓，活动轮廓的外力使用梯度矢量流场；

步骤i：将多尺度图像中活动轮廓进化得到的轮廓通过变换及加权，得到管状结构的最终内层轮廓。

优选的实施例，在纵向截平面提取内层轮廓线的步骤包括：

步骤b1：从管状结构三维序列数据中等角度的方向上提取N，N∈[3，6]个纵向截平面，N为需要处理的纵向截平面的个数。

步骤b2：移除纵向截平面中不包含有用信息的导管区域，得到两个纵向子截平面；

步骤b3：使用a*a，a＝{3，5，7，9}的滑动窗口，滑动窗口中心点遍历纵向子截平面的每一个像素点，并统计窗口内像素灰度的均值和方差，并以统计得到的均值和方差作为特征进行模糊聚类，a为滑动窗口的边长；

步骤b4：利用模糊聚类把每一个纵向子截平面分为两类，一类为内层区域，另一类为管壁及其它区域，从这两类的边界可提取管状结构在纵向截平面的内层轮廓线。

优选的实施例，其中，提取最终内层轮廓的步骤包括：

步骤i1：将不同尺度的活动轮廓进化得到的轮廓变换到第一尺度(即原始图像所在的尺度)，得到不同尺度的内层轮廓；

步骤i2：对不同尺度的内层轮廓分别进行形态学操作，填充轮廓，得到内层区域；

步骤i3：对不同尺度的内层区域进行距离变换，得到内层区域距离图D_i，i＝1，2，3；

步骤i4：对内层区域距离图D_i，i＝1，2，3进行加权，通过以下公式进行加权：

D＝ω₁D₁+ω₂D₂+ω₃D₃，

其中ω₁、ω₂、ω₃为加权因子，加权因子的取值范围为：0≤ω₁≤1，0≤ω₂≤1，0≤ω₃≤1，ω₁+ω₂+ω₃＝1，并设定一个阈值T，min{ω₁，ω₂，ω₃}≤T≤1，最终内层区域可以通过公式D_final＝D＞T得到，大于该阈值的区域即为管状结构的最终内层区域，小于该阈值的区域被舍弃，D_final为最终内层区域；最终内层轮廓即为最终内层区域的边界。

本发明的有益效果：本发明利用模糊聚类在纵向截平面上提取管状结构内层的轮廓线，然后把纵向截平面的轮廓线变换到横向截平面，进行轮廓的初始化。不仅避免了手工初始化管状结构内层的轮廓，而且大大提高了初始化管状结构内层轮廓的速度。活动轮廓算法要尽可能减少噪声的干扰，然而图像噪声的存在，导致活动轮廓不能得到精确的内层轮廓。为了减少噪声对活动轮廓的影响，首先，使用各向异性滤波去除随机噪声，各向异性扩散滤波能在去除随机噪声的同时还能保持边缘信息并能增强边缘；其次，利用小波变换产生多尺度图像，并在不同尺度上利用活动轮廓进化轮廓。最后，把不同尺度上得到的轮廓加权变换到原始图像尺度，得到管状结构内层的最终轮廓。本发明利用了模糊聚类和活动轮廓在多尺度上提取管状结构的内层轮廓，在保持了抗噪性能好的情况下，实现了初始化的自动化，从而实现了整个分割方法的自动化。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的不当限定，在附图中：

图1是本发明的管状结构内层分割方法的实现流程图；

图2是本发明实施例中提取纵向截平面的示意图；

图3a-图3f是本发明实施例中在纵向截平面上提取管状结构内层轮廓线的示意图；

图4是本发明实施例中纵向截平面管状结构内层轮廓点与横向截平面对应关系的示意图；

图5a-图5c是本发明实施例中样条插值的示意图；

图6a-图6c是本发明实施例中多尺度图像中活动轮廓进化示意图；

图7是本发明实施例中提取管状结构最终内层轮廓示意图。

具体实施方式

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

作为一种具体的实现方案，使用C++语言实现了本发明所描述的算法，所有的实现都是在一台酷睿2双核2.2G，2G内存，操作系统为Windows XP的计算机上完成的。图1示出了可用于实施在此公开的方法的流程图。该方法包括：

步骤100：读入三维序列数据；

步骤101：各向异性扩散滤波如下所述：图像中常含有随机噪声，因此图像预处理是图像分析与识别必要的前期工作。为了减少噪声对分割算法的影响，采用各向异性扩散滤波来去除管状结构三维序列数据的随机噪声，同时还能保持边缘信息和增强边缘位置。各向异性扩散滤波(如下所述：P.Perona and J.Malik，“Scale-space and edge detection usinganisotropic diffusion，”IEEE Trans.Pattern Anal.Machine Intell.，vol.12，no.7pp.629-639，July 1990)最早是由Perona和Malik提出的，可以用下面的微分模型表示：

$\frac{\∂I(x,y,t)}{\∂t}=\mathrm{div}[g\left(\right||\▿I|\left|\right)\▿I],$

其中x是图像的横坐标，y是图像的纵坐标，t是时间变量，div是散度算子，

是图像的梯度。扩散系数是关于梯度幅值的非负单调递减函数g(·)，满足 $g\left(0\right)=1,\underset{x\→\∞}{\lim }g\left(x\right)=0,$ 这样可以保证区域内部(梯度小)，扩散系数大；在边缘(梯度大)，扩散系数小，从而起到边缘增强的作用。由于Perona和Malik提出的扩散系数比较大，可能会消除一些不明显而又重要的边缘细节，本发明采用基于Tukey的Biweight Estimator误差模型(如下所述：Black M.J.，Snpiro G..Marimont D.H.，Heeger D.，“Robustanisotrupic diffusiun”.IEEE Truns.On Image Pmc.，Vol.7.No 3，pp.421-432，March 1998)的边缘终止函数为：

作为扩散系数，其中σ是梯度阈值，基于Tukey函数比Perona-Malik函数递减程度更快，当梯度超过σ时，扩散系数为0，即模型不再扩散，从而避免了边缘弱化问题。

步骤102：从管状结构三维序列数据中提取纵向截平面，并在纵向截平面上利用模糊聚类提取管状结构内层的轮廓线如下所述：模糊C均值聚类(FCM)算法是Bezkek(如下所述：J.C.Bezdek，Pattern recognitionwith fuzzy objective function algorithms，Plenum Press，New York；1981)于1981年提出的。其基本思想是：{x_i，i＝1，2，…，N}是N个样本组成的样本集合X，C为预定的类别数目，c_i，i＝1，2，…，C为每个聚类的中心，μ_ij是第i个样本对第j个样本的隶属度函数。目标准则函数可以写为

$J_m=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{C}{\mathrm{\Σ}}}{u}_{\mathrm{ij}}^m{\left|\right|x_i-c_j\left|\right|}^2,$ 其中m＞1是一个可以控制聚类结果的模糊程度的常数。模糊C均值算法是通过对目标准则函数的迭代来获取对数据集的模糊分类，即迭代

$u_{\mathrm{ij}}=\frac{1}{\underset{k=1}{\overset{C}{\mathrm{\Σ}}}{\left[\frac{\left|\right|x_i-c_j\left|\right|}{\left|\right|x_i-c_k\left|\right|}\right]}^{\frac{2}{m-1}}},$ $c_j=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{u}_{\mathrm{ij}}^m\·x_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{u}_{\mathrm{ij}}^m},$ 使目标准则函数收敛到一个局部极小值或鞍点，得到X的一个最优模糊C划分。u_ij是第i个样本对第j个样本的隶属度函数，C为聚类数目，c_j为第j个聚类中心，c_k为第k个聚类中心，x_i为第i个样本，m为控制聚类结构的模糊程度的常数，且m＞1，N为样本的数目。

由于噪声的影响，仅利用像素的灰度信息对图像进行聚类，聚类结果对噪声较为敏感。为此，使用一个a*a，a＝{3，5，7，9}滑动窗口，其中a为滑动窗口的边长，滑动窗口的中心点遍历纵向截平面的每一个像素，并统计窗口内的均值和方差作为模糊聚类的特征x_i，即x_i＝{x_iS，σ_iS}，x_iS为窗口内像素的均值，σ_iS为窗口内像素的方差。以滑动窗口内的均值和方差作为特征进行模糊聚类，不仅利用了像素的灰度信息，还利用了窗口内像素的方差信息，对噪声有一定的抑制能力。

活动轮廓分割图像需要初始轮廓。然而管状结构三维序列数据有上百张乃至上千张的图像，如果每张都手动描绘初始轮廓，这样工作量大，耗时多，而且可重复性差。有些半自动的方法需要初始化第一张图像的轮廓，然后利用上一张图像的分割结果作为下一张图像的初始轮廓，如果相邻两张图像的轮廓变化较大，这样会导致初始化位置不正确，可能出现分割错误。为此，本发明利用模糊聚类在纵向截平面上提取管状结构内层的轮廓线，然后变换到横向截平面初始轮廓。每一张图像的初始轮廓都是独立的，每一张图像的初始轮廓不会受相邻图像初始轮廓的影响，包括以下步骤：

步骤102a：请参考图2提取纵向截平面的示意图，提取三维序列数据的纵向截平面是从管状结构三维序列数据20中等角度的方向上提取N，N∈[3，6]个纵向截平面21，22，N为需要处理的纵向截平面的个数。

步骤102b：移除纵向截平面的导管区域如下所述：由于纵向截平面21，22的图像中包括三个区域：导管区域，内层区域、管壁和其它区域。导管区域不包含任何有用的信息，可以移除导管区域。纵向截平面22移除导管区域后得到左、右两个纵向子截平面30。

步骤102c：请参考图3a-图3f本发明实施例中在纵向截平面22上提取管状结构内层轮廓线的示意图；其中：如图3a示出一个纵向截平面，在纵向截平面上进行模糊聚类如下所述：使用a*a，a＝{3，5，7，9}的滑动窗口，其中a为滑动窗口的边长；如图3b示出移除导管区域后得到的两个纵向子截平面，滑动窗口的中心点遍历纵向子截平面30的每一个像素，并统计窗口内像素的均值和方差。并以统计得到的均值和方差作为特征进行模糊聚类。每个纵向子截平面聚为两类，一类为内层区域，一类为管壁和其它区域。如图3c示出了模糊聚类的结果。由于斑点噪声的影响，虽然使用了均值和方差作为特征进行聚类，但是聚类结果中可能还会出现噪声点。因此使用形态学操作消除噪声点和不规则的边缘，聚类后的图像31经过形态学操作得到纵向截平面管状结构内层的分割结果32，如图3d示出了形态学操作后的图像。

步骤102d：提取内层轮廓线如下所述：利用模糊聚类把每一个纵向子截平面分为两类，一类为内层区域，另一类为管壁及其它区域，从这两类的边界可提取管状结构在纵向截平面的内层轮廓线。续请参考图3d，对于左子截平面321，从右到左扫描图像，当像素从黑色变为白色时，该白色点即为轮廓点，然后跳转到下一行扫描。对于右子截平面322，从左到右进行扫描。两个子截平面都处理完毕后，得到纵向截平面上的内层轮廓线33，如图3e示出提取的管状结构内层轮廓线。轮廓线叠加到原始纵向截平面图像34，如图3f示出管状结构内层轮廓线叠加到原始纵向截平面。

步骤103：把纵向截平面提取的轮廓线变换到横向截平面，每一个横向截平面得到对应的轮廓点；如图4示出本发明实施例中纵向截平面管状结构内层轮廓点和横向截平面对应关系的示意图；将纵向截平面内层轮廓线变换到横向截平面如下所述：为了能在横向截平面使用活动轮廓，需要在横向截平面上初始轮廓，因此，需要把纵向截平面提取的轮廓线变换到横向截平面。纵向截平面提取的轮廓线40与横向截平面42有唯一的对应关系，两个“☆”标志为同一轮廓点在纵向截平面和横向截平面中的位置41。本发明提取了管状结构三维序列数据中等角度的N＝5个纵向截平面，每一个纵向截平面能提取两条内层轮廓线。把五个纵向截平面的轮廓线变换到横向截平面，每一个横向截平面会得到对应的十个轮廓点50。

步骤104：如图5a-图5c示出本发明实施例中样条插值的示意图；其中：图5a示出横向截平面的轮廓点，对横向截平面上的轮廓点进行样条插值，得到一条封闭的样条轮廓；如下所述：纵向截平面的轮廓线变换到横向截平面后，只有十个轮廓点50，初始轮廓需要的是一条封闭的曲线，而不是独立的点。因此，需要对轮廓点进行样条插值。三次样条在灵活性和计算速度之间提供了一个合理的折衷方案，与低次多项式相比，三次样条在模拟任意曲线时显得更灵活，与更高次多项式相比，三次样条只需较少的计算和存储空间。因此，初始轮廓使用Cardinal三次样条(如下所述：Donald Hearn，M Pauline Baker.ComputerGraphics[M].Prentice Hall Press Company.1998)进行插值得到的样条轮廓。Cardinal样条完全由连续4个控制点给出，中间2个点是曲线的端点，第1个点和第4个点用来计算斜率。插值公式可以由以下多项式表示：

P(u)＝P_k-1(su³+2su²-su)+P_k[(2-s)u³+(s-3)u²+1]，

      +P_k+1[(s-2)u³+(3-2s)u²+su]+P_k+2(su³-su²)

其中s＝(1-t)/2，t为张力参数，u为点k和点k+1间的相对偏移量，u∈[0，1]，P(u)为相对偏移量是u时，插值点的位置。P_k，P_k+1为曲线的端点，P_k-1，P_k+2为连续四个点的第1个点和第4个点。轮廓点使用上述公式进行插值得到样条轮廓51，如图5b示出Cardinal三次样条插值得到的样条轮廓。把样条轮廓叠加到原始横向截平面图像52，如图5c示出样条插值的轮廓叠加到原始横向截平面。

步骤105：利用小波变换产生横向截平面的多尺度图像，用来在多尺度图像中使用活动轮廓；小波变换得到多尺度图像如下所述：轮廓进化要尽量减少噪声的干扰，并同时尽量多保留感兴趣的边缘。但是两者往往不能同时兼顾。多尺度分析提供了有效的解决方法。小尺度时图像比较精细，但是噪声也比较严重，大尺度时图像比较粗略，但噪声得到抑制。因此，可以使用多尺度的方法提取管状结构的内层轮廓。使用小波变换生成k层小波系数A_k，k＝1，2。将原始横向截平面图像记为第一尺度图像42，A₁为第二尺度图像61，A₂为第三尺度图像63，则第一尺度图像42、第二尺度图像61、第三尺度图像63可以看成是由精细到粗略的多尺度图像。

步骤106：把样条轮廓变换到不同尺度，用来作为多尺度图像中活动轮廓的初始轮廓；将样条轮廓变换到不同尺度如下所述：步骤104得到了一个样条轮廓51，这个样条轮廓51只能作为第一尺度图像42的初始轮廓，第二尺度图像61、第三尺度图像63的初始轮廓还没有得到。因此，需要把样条轮廓51变换到第二尺度和第三尺度。将样条轮廓51变换到第二尺度，得到第二尺度图像的初始轮廓60，将样条轮廓51变换到第三尺度，得到第三尺度图像的初始轮廓62。

步骤107：计算横向截平面的多尺度图像的梯度矢量流场。梯度矢量流场(GVF Field)(如下所述：XU C Y，P rince J L.Snakes，shapes，andgradient vector flow[J].IEEE T ransactions on Imaging Processing，1995)可以表示为矢量场v(x，y)＝[u(x，y)，v(x，y)]，并且满足以下能量最小化函数：

$\mathrm{\ϵ}=\∫\∫\mathrm{\μ}(u_x^2+u_y^2+v_x^2+v_y^2)+{|\▿f|}^2{|v-\▿f|}^2\mathrm{dxdy},$

f是边缘图像，

是边缘图像的梯度，v为矢量场，u，v分别为矢量场的第一和第二个分量，u_x，u_y，v_x，v_y是u，v分别对x，y求一阶偏导，ε为能量最小化函数，μ为平衡能量最小化函数中前后两项权重的系数，μ随着噪声的增加而增大。依据变分法，梯度矢量流场可以通过求解下列Euler方程得到：

$\mathrm{\μ\Δu}-(u-f_x){|\▿f|}^2=0$

                     ，

$\mathrm{\μ\Δv}-(v-f_y){|\▿f|}^2=0$

其中，f_x，f_y为边缘图像f对x，y的偏导，Δ是拉普拉斯算子。以梯度矢量流场作为外力极大地增加了活动轮廓的捕获范围。每一个尺度的图像都需要计算对应尺度的梯度向量流场。

步骤108：在横向截平面的多尺度图像中利用活动轮廓进化轮廓；如图6a-6c示出了多尺度图像中活动轮廓进化示意图，在多尺度图像中利用活动轮廓进化轮廓如下所述：

活动轮廓是由Kass等人提出的(如下所述：M.Kass，A.Witkin，andD.Terzopoulos，“Snakes：Active contourmodels，”Int.J.Comput.Vis.，vol.1，pp.321-331，1987)，活动轮廓定义为一条参数化曲线X(s)＝(x(s)，y(s))，s∈[0，1]，曲线在图像中移动并满足下面最小化能量函数：

$E={\∫}_0^1\frac{1}{2}[\mathrm{\α}{\left|X^{\′}\left(s\right)\right|}^2+\mathrm{\β}{\left|X^{\′\′}\left(s\right)\right|}^2]+E_{\mathrm{ext}}\left(X\left(s\right)\right)\mathrm{ds},$

其中，E为最小化能量函数，X(s)为参数化曲线，α为加权因子，β为加权因子，E_ext(X(s))为与图像有关的外力，X′(s)为关于参数s的一阶导数，X″(s)为关于参数s的二阶导数，最小化能量函数E必须满足以下Euler方程：

$\mathrm{\α}{X}^{\′\′}\left(s\right)-\mathrm{\β}{X}^{\′\′\′\′}\left(s\right)-\▿E_{\mathrm{ext}}=0,$

其中X″(s)为关于参数s的二阶导数，X″″(s)为关于参数s的四阶导数，

为梯度算子，取 $\▿E_{\mathrm{ext}}=v(x,y),$ v(x，y)为步骤107计算的梯度矢量流场。

参数化曲线的初始位置即为活动轮廓的初始轮廓，步骤106可以得到不同尺度活动轮廓的初始轮廓。通过求解上述的Euler方程，可以得到轮廓进化的结果。

由于噪声的干扰，单纯的应用活动轮廓并不能得到精确的轮廓。因此，本发明在三个不同尺度上使用活动轮廓进化轮廓。三个尺度的轮廓进化可以同时进行，这样可以使用并行程序实现活动轮廓进化，减少了轮廓进化时间。第一尺度图像42和第一尺度的初始轮廓51使用活动轮廓进化得到的轮廓611，如图6a示出第一尺度活动轮廓进化的示意图、第二尺度图像61和第二尺度的初始轮廓60使用活动轮廓进化得到的轮廓621，如图6b示出第二尺度活动轮廓进化的示意图、第三尺度图像63和第三尺度的初始轮廓62使用活动轮廓进化得到的轮廓631，如图6c示出第三尺度活动轮廓进化的示意图。

步骤109：如图7示出了提取管状结构最终内层轮廓示意图，将多尺度图像中活动轮廓进化得到的轮廓通过变换及加权，得到管状结构的最终内层轮廓如下所述：不同尺度的图像有不同的特征，小尺度图像比较精细，大尺度图像比较粗略。因此，本发明结合不同尺度的分割结果，得到最终内层的轮廓。首先，将不同尺度的活动轮廓进化得到的轮廓变换到第一尺度(即原始图像所在的尺度)，由于不同尺度图像的大小不一样，为了结合不同尺度活动轮廓进化得到的轮廓，需要把第二尺度和第三尺度活动轮廓进化得到的轮廓变换到第一尺度上，第一尺度活动轮廓进化得到的轮廓不需要变换，因此第一尺度内层轮廓611即为第一尺度活动轮廓进化得到的轮廓611，第二尺度活动轮廓进化得到的轮廓621变换到第一尺度，得到第二尺度内层轮廓622；第三尺度活动轮廓进化得到的轮廓631变换到第一尺度，得到第三尺度内层轮廓632。

其次，对不同尺度的内层轮廓分别进行形态学操作，填充轮廓，得到内层区域。第一尺度内层轮廓611进行区域填充，得到第一尺度内层区域72；第二尺度内层轮廓622进行区域填充，得到第二尺度内层区域71；第三尺度内层轮廓632进行区域填充，得到第三尺度内层区域70；

然后，对不同尺度的内层区域进行距离变换，得到内层区域距离图D_i，i＝1，2，3，D₁为第一尺度内层区域距离图75，D₂为第二尺度内层区域距离图74，D₃为第三尺度内层区域距离图73。

最后，对内层区域距离图D_i，i＝1，2，3进行加权，通过以下公式进行加权：

D＝ω₁D₁+ω₂D₂+ω₃D₃，

其中ω₁、ω₂、ω₃为加权因子，加权因子的取值范围为：0≤ω₁≤1，0≤ω₂≤1，0≤ω₃≤1，ω₁+ω₂+ω₃＝1，并设定一个阈值T，min{ω₁，ω₂，ω₃}≤T≤1，最终内层区域可以通过公式D_final＝D＞T得到，大于该阈值的区域即为管状结构的最终内层区域，小于该阈值的区域被舍弃，D_final为最终内层区域。通过阈值判断，得到管状结构最终内层区域77。最终管状结构的内层轮廓78即为管状结构最终内层区域77的边界。

本发明采用的是从不同纵向截平面提取内层轮廓线，并把纵向截平面的轮廓线变换到横向截平面来初始化活动轮廓的轮廓，而不是使用上一张图像的轮廓作为下一张图像的初始轮廓，这样避免了初始化错误，提高了分割的精度。由于每张图像都有一个初始轮廓，可以使用并行程序同时对多张图像的轮廓进行进化。本发明还采用了多尺度分割管状结构内层来提高分割精度。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (3)

1.一种对管状结构内层分割的方法，其特征在于，包括步骤：

步骤a：对管状结构三维序列数据进行各向异性扩散滤波，用于除去随机噪声；

步骤b：从管状结构三维序列数据中提取纵向截平面，并在纵向截平面上利用模糊聚类提取管状结构内层的轮廓线；

步骤c：把纵向截平面提取的轮廓线变换到横向截平面，每一个横向截平面得到对应的轮廓点；

步骤d：对横向截平面上的轮廓点进行样条插值，得到一条封闭的样条轮廓；

步骤e：利用小波变换产生横向截平面的多尺度图像，用来在多尺度图像中使用活动轮廓，多尺度图像包含第一尺度、第二尺度和第三尺度；

步骤f：把样条轮廓变换到不同尺度，用来作为多尺度图像中活动轮廓的初始轮廓；

步骤g：计算横向截平面的多尺度图像的梯度矢量流场；

步骤h：在横向截平面的多尺度图像中利用活动轮廓进化轮廓，活动轮廓的外力使用梯度矢量流场；

步骤i：将多尺度图像中活动轮廓进化得到的轮廓通过变换及加权，得到管状结构的最终内层轮廓。

2.按照权利要求1所述方法，其特征在于，在纵向截平面提取内层轮廓线的步骤包括：

步骤b1：从管状结构三维序列数据中等角度的方向上提取N，N∈[3，6]个纵向截平面，N为需要处理的纵向截平面的个数。

步骤b2：移除纵向截平面中不包含有用信息的导管区域，得到两个纵向子截平面；

步骤b3：使用a*a，a＝{3，5，7，9}的滑动窗口，滑动窗口中心点遍历纵向子截平面的每一个像素点，并统计窗口内像素灰度的均值和方差，并以统计得到的均值和方差作为特征进行模糊聚类，a为滑动窗口的边长；

步骤b4：利用模糊聚类把每一个纵向子截平面分为两类，一类为内层区域，另一类为管壁及其它区域，从这两类的边界可提取管状结构在纵向截平面的内层轮廓线。

3.按照权利要求1所述方法，其特征在于，所述提取最终内层轮廓的步骤包括：

步骤i1：将不同尺度的活动轮廓进化得到的轮廓变换到第一尺度(即原始图像所在的尺度)，得到不同尺度的内层轮廓；

步骤i2：对不同尺度的内层轮廓分别进行形态学操作，填充轮廓，得到内层区域；

步骤i3：对不同尺度的内层区域进行距离变换，得到内层区域距离图D_i，i＝1，2，3；

步骤i4：对内层区域距离图D_i，i＝1，2，3进行加权，通过以下公式进行加权：

D＝ω₁D₁+ω₂D₂+ω₃D₃，

其中ω₁、ω₂、ω₃为加权因子，加权因子的取值范围为：0≤ω₁≤1，0≤ω₂≤1，0≤ω₃≤1，ω₁+ω₂+ω₃＝1，并设定一个阈值T，min{ω₁，ω₂，ω₃}≤T≤1，最终内层区域可以通过公式D_final＝D＞T得到，大于该阈值的区域即为管状结构的最终内层区域，小于该阈值的区域被舍弃，D_final为最终内层区域；最终内层轮廓即为最终内层区域的边界。

Images