超声颈动脉斑块自动分割方法技术领域
本发明属于计算机技术与医学图像的交叉领域，具体涉及到一种超声图像中颈动脉血管横截面方向上斑块的自动分割方法。
背景技术
根据世界卫生组织统计的数据，每年心血管疾病导致的死亡人数约占全世界死亡总人数的三分之一，颈动脉粥样硬化与心脑血管疾病密切相关。与传统的内中膜厚度（Intima-Media Thickness,IMT）度量指标相比，斑块总体积（Total Plaque Volume,TPV）、斑块面积（Total Plaque Area,TPA）以及斑块成分分析等指标更能准确直观地反映动脉粥样硬化状况，成为近年来预测心脑血管疾病风险的重要指标之一。
在超声颈动脉横截面图像上，颈动脉血管腔表现为低回声的均质区域，外膜表现为高回声的非均质带状区域，血管腔和外膜之间为血管内膜、中膜以及存在的斑块。图1显示了超声成像下的颈动脉横截面及其结构示意图。图1中实心点围成的轮廓是颈动脉血管腔与内膜边界（一般简称为内膜边界），其内部较为均匀的黑色低回声部分是血管腔（lumen）；十字形点围成的闭合轮廓线是颈动脉血管外膜与中膜边界（一般简称为外膜边界），其外部包含高回声的不均匀的区域是颈动脉外膜，内膜边界与外膜边界间区域包括回声较低的健康内膜、中膜以及回声不均匀且相对较高的斑块；图中实线所围成的区域即为斑块，其内部灰度分布不均匀，且与周围的灰度特性相近，不易区分。
传统斑块边界的获得主要是通过操作者（医生）肉眼观察手工勾勒的方法。其主要缺陷在于它很大程度上取决于操作者的主观判断及操作经验， 因此在不同操作者对同一个目标边界的判断，甚至同一个操作者在不同的时间点对同一目标边界的判断也会有差异。此外，手工勾勒边界非常耗时，当一幅超声颈动脉图像上存在多个斑块时，经验丰富的操作者勾勒出斑块边界大致需要2~3分钟，无法满足临床病例分析的需求。
现有的斑块分割方法大多应用在核磁共振（Magnetic Resonance,MR）图像、血管内超声（Intra-vascular Ultrasound,IVUS）图像、超声RF（Radio Frequency）图像和CT血管造影（CT Angiography,CTA）图像中。斑块分割算法主要分为两类，聚类方法（如Mahdi Mazinani et al.,Automatic Segmentation of Soft Plaque by Modeling the Partial Volume Problem in the Coronary Arterv,2010 Fourth International Conference on Digital Society;Francois Destrempes et al.，Segmentation of Plaques in Sequences of Ultrasonic B-Mode Images of Carotid Arteries Based on Motion Estimation and a BayesianModel,IEEE Trans.Biomed.Eng.2011.58(8):p.2202-2211;）和几何活动轮廓模型方法（如Oliver Gloger et al.,A three-stepped coordinated level set segmentation method for identifying atherosclerotic plaques on MR-images,Commum.Numer.Meth.Engng 2009.25:p.615-638;Shawn Lankton et al.,Soft Plaque Detection and Automatic Vessel Segmentation,MICCAI 2009）。聚类算法以相似性为基础，其方法是通过灰度或其他特征，由初始聚类中心开始，寻找具有较高相似度的同一类对象，从而得到最终斑块轮廓。这种算法存在对聚类中心敏感、不适用于混合斑块分割的问题。几何活动轮廓模型方法是将灰度、位置等信息构成使轮廓变化的能量项，从而找到准确的目标轮廓。上述两篇文章，都是在分割得到内外膜轮廓后，将其向内收缩一定距离作为初始轮廓，然后演化得到精确的分割结果，但都只考虑灰度特征对混合斑块的分割，具有局限性，不能获得理想的结果。上述所有分割方法是针对核磁共振图像、CT图像和血管内超声图像，至今还没有任何关于B超颈动脉图像中斑块分割方法的报道。
发明内容
本发明的目的在于提供一种普通超声图像中基于血管膜分割的颈动脉斑块的自动分割方法，能够有效分割出颈动脉斑块，为颈动脉粥样硬化的病变观察和药物治疗提供分析参数。
本发明提供一种超声颈动脉斑块自动分割方法，包括以下步骤：
（1）从颈动脉三维超声体数据提取当前帧图像；
（2）对当前帧图像分割得到内、外膜轮廓；
（3）从当前帧图提取包含内、外膜轮廓的斑块感兴趣区域ROI，感兴趣区域ROI的中心点与内膜轮廓的中心点相同；
（4）在斑块感兴趣区域ROI检测初始斑块外边界C0：
（4.1）以斑块感兴趣区域ROI的中心点为原点，将斑块感兴趣区域ROI从直角坐标系转换到极坐标系；
（4.2）对极坐标系下的斑块感兴趣区域ROI中每列上位于内外膜间的像素点灰度值进行曲线拟合，并在拟合得到的曲线上检测极小值，从极小值对应的像素点中选取离外膜轮廓最近的像素点；
（4.3）将步骤（4.2）选取的离外膜轮廓最近的像素点从极坐标系转换到直角坐标系，并连线形成初始斑块外边界C0；
（5）依据初始斑块外轮廓线C0进行水平集演化得到斑块区域：
（5.1）初始化迭代次数z=1以及符号距离函数
outside(C0)表示斑块感兴趣区域ROI中初始斑块外边界以外的像素区域，x为斑块感兴趣区域ROI的像素点，|| ||为求欧式距离；
（5.2）计算φz(x)：
φ z ( x ) = φ z - 1 ( x ) + δ ( φ z - 1 ( x ) ) · { λ s div ( ▿ φ z - 1 ( x ) | ▿ φ z - 1 ( x ) | ) + λ l ∫ Ω B L ( x , x ^ ) δ ( φ z - 1 ( x ^ ) ) [ ( I ( x ^ ) - u x ) 2 - ( I ( x ^ ) - v x ) 2 ] d x ^ + λ g · 1 N · Σ i = 1 N [ ( u 0 i ( x ) - c + i ) 2 - ( u 0 i ( x ) - c - i ) 2 ] + λ d · B d ( x ) } ]]>
其中，
半径为r1的圆形邻域 为斑块感兴趣区域ROI的像素点，2≤rL≤20；
λs、λl、λg、λd分别为限制曲线长度能量、局部Chan-Vese能量、Gabor滤波器能量、距离限制能量的权重；
I()为像素灰度；
u x = ∫ Ω B L ( x , x ^ ) H ( φ ( x ^ ) ) I ( x ^ ) d x ^ / ∫ Ω B L ( x , x ^ ) H ( φ ( x ^ ) ) d x ^ ; ]]>
v x = ∫ Ω B L ( x , x ^ ) ( 1 - H ( φ ( x ^ ) ) ) I ( x ^ ) d x ^ / ∫ Ω B L ( x , x ^ ) ( 1 - H ( φ ( x ^ ) ) ) d x ^ ; ]]>
Ω为斑块感兴趣区域ROI；
N为Gabor滤波器的个数；
为斑块感兴趣区域ROI经过第i个Gabor滤波器滤波后的输出；
为位于内外膜之间且位于上一次迭代演化轮廓Cz-1内部的像素点经过第i个Gabor滤波得到的响应均值；
为位于内外膜之间且位于上一次迭代演化轮廓Cz-1外部的像素点经过第i个Gabor滤波得到的响应均值；
msd为外膜轮廓CMAB和内膜轮廓CLIB之间的欧式距离，d为内外膜轮廓间的最小距离。
（5.3）计算本次迭代演化轮廓Cz={x|φz(x)=0}；
（5.4）若本次迭代演化轮廓Cz与前一次迭代得到的演化轮廓Cz-1相同或者z达到迭代次数上限，则本次迭代得到的演化轮廓即为最终斑块外边界Cplaque，进入步骤（5.5）；否则，z=z+1，返回步骤（5.2）；
（5.5）最终斑块演化轮廓Cplaque与内膜轮廓CLIB间区域即为斑块区域。
按上述方案，所述步骤（2）具体为：
（2.1）若当前帧图像为颈动脉三维超声体数据的第一帧图像，则在当前帧图像上选择明显位于血管外轮廓上的像素点集作为基准点，通过基准点间插值形成闭合曲线，将其作为当前帧的基准轮廓；否则，将上一帧得到的备用基准轮廓作为当前帧的基准轮廓；
（2.2）利用形态学膨胀法将基准轮廓外推形成第一感兴趣区域ROI1；
（2.3）在第一感兴趣区域ROI1中检测颈动脉血管外轮廓CMAB：
（2.3.1）在当前帧图像中提取包含第一感兴趣区域ROI1的矩形窗口图像，将其转化到极坐标图像，在极坐标图像上每隔相同角度在径向方向靠近基准轮廓的像素点中搜索梯度值最大者作为初始轮廓点，将初始轮廓点转化到直角坐标系并连线形成初始颈动脉血管外轮廓
（2.3.2）利用混合分布估计第一感兴趣区域ROI1的灰度概率密度分布，从而得到混合分布的shape参数Kj和scale参数θj，j=1,…,M，M为混合分布中单项分布的类别数；
（2.3.3）初始化迭代次数t＝1以及符号距离函数
表示第一感兴趣区域ROI1中初始颈动脉血管外轮廓 以外的像素区域，x1为第一感兴趣区域ROI1中的像素点，|| ||为求欧式距离；
（2.3.4）计算
φ 1 t ( x 1 ) = φ 1 t - 1 ( x 1 ) - ρ 1 · δ ( φ 1 t - 1 ( x 1 ) ) · [ ln ( Σ j M q Ai ( j ) · G ( I x 1 | K j , θ j ) - ln ( Σ j M q Ae ( j ) · G ( I x 1 | K j , θ j ) ] + λ 1 · δ ( φ 1 t - 1 ( x 1 ) ) · div ( ▿ φ 1 t - 1 ( x 1 ) | ▿ φ 1 t - 1 ( x 1 ) | ) ]]>
其中， 这里0<ε<0.001，ρ1>0，λ1>0，qAi(j)和qAe(j)分别为前一次迭代演化轮廓 的内、外区域中第j项单项分布的权重， 为参数为Kj,θj的第j项单项分布函数， 为像素点x1的灰度值；
（2.3.5）计算本次迭代演化外轮廓
C 1 t = { x 1 | φ 1 t ( x 1 ) = 0 } ]]>
（2.3.6）若本次迭代演化外轮廓 与前一次迭代得到的演化外轮廓 相同或者t达到迭代次数上限，则本次迭代得到的演化轮廓即为颈动脉血管外轮廓CMAB，迭代结束；否则，t=t+1，返回步骤（2.3.4）；
（2.4）在当前帧图像中将颈动脉血管外轮廓CMAB围成的区域作为第二感兴趣区域ROI2；
（2.5）在第二感兴趣区域ROI2中检测颈动脉血管内轮廓CLIB：
（2.5.1）利用形态学腐蚀法将颈动脉血管外轮廓CMAB内推，将内推得到 的轮廓线作为初始颈动脉血管内轮廓
（2.5.2）初始化迭代次数t′＝1以及符号距离函数
表示第二感兴趣区域ROI2中初始颈动脉血管内轮廓 以外的像素区域，x2为第二感兴趣区域ROI2中的像素点；
（2.5.3）计算
φ 1 t ′ ( x 2 ) = φ 2 t ′ - 1 ( x 2 ) - ρ 2 · δ ( φ 2 t ′ - 1 ( x 2 ) ) · [ ln ( Σ j M q Li ( j ) · G ( I x 2 | K j , θ j ) - ln ( Σ j M q Le ( j ) · G ( I x 2 | K j , θ j ) ] + λ 2 · δ ( φ 2 t ′ - 1 ( x 2 ) ) · div ( ▿ φ 2 t ′ - 1 ( x 2 ) | ▿ φ 2 t ′ - 1 ( x 2 ) | ) + β 2 · δ ( φ 2 t ′ - 1 ( x 2 ) ) · B T ( x 2 ) ]]>
其中，ρ2>0，λ2>0，β2>0，qLi(j)和qLe(j)分别为前一次演化轮廓 的内、外区域中第j项单项分布的权重， 为参数为Kj,θj的第j项单项分布函数， 为像素点x2的灰度值；
（2.5.4）计算本次迭代演化内轮廓
C 2 t ′ = { x 2 | φ 2 t ′ ( x 2 ) = 0 } ]]>
（2.5.5）若本次迭代演化内轮廓 与前一次迭代得到的演化内轮廓 相同或者达到迭代次数上限，则本次迭代得到的演化轮廓为颈动脉血管内轮廓CLIB，迭代结束；否则，t′=t′+1，返回步骤（2.5.3）；
（2.6）跟踪确定下一帧图像的备用基准轮廓：
（2.6.1）在当前帧图像中，通过形态学膨胀法将颈动脉血管外轮廓CMAB外推，将外推得到的闭合轮廓线与CMAB构成的区域作为跟踪区域Ω′；
（2.6.2）在下一帧图像中搜索使得灰度差异和S=∑Ω′ΔI2最小的像素点集(x4,y4)，其中ΔI2=[I(x3,y3)-J(x4,y4)]2，I(x3,y3)为当前帧图像的跟踪区域Ω′ 中像素点(x3,y3)的灰度值，J(x4,y4)为下一帧图像的像素点(x4,y4)的灰度值；
（2.6.3）令当前帧图像的跟踪区域中的点(x3,y3)与下一帧图像的像素点(x4,y4)的映射关系为 x 4 y 4 = τ x τ y + 1 0 0 1 x 3 y 3 , ]]>结合像素点(x3,y3)与(x4,y4)求解映射关系中的参数τ=(τx,τy)从而确定映射关系；
（2.6.4）根据步骤（2.6.3）建立的映射关系将当前帧图像的颈动脉血管外轮廓CMAB映射为下一帧图像的候选基准轮廓。
按所述方案，步骤（2.3.2）中混合分布为混合gamma分布或混合gauss分布或混合rayleigh分布。
本发明提供的超声图像中颈动脉斑块的自动分割算法有以下几个创新点：1.利用血管内、外膜作为先验实现了颈动脉斑块的全自动分割；2.利用斑块和外膜部分比健康的中膜组织灰度值高的特点，在内外膜之间找到极小值点以确定斑块初始外轮廓；3.根据斑块与内外膜相对位置的先验信息设计水平集的距离限制项，保证了检测出的斑块位置的合理性；4.利用Gabor滤波器得到的特征来设计水平集中的全局项，克服了只用灰度值来分割的局限性。
本发明提供的超声图像中颈动脉斑块的计算机自动分割算法能够达到以下目标：（1）能有效应对超声图像中存在的噪声及弱边界情况；（2）能精确分割出血管斑块；（3）能较大幅度地减少医生的工作量；（4）基于本方法得到的度量指标能够提供准确有效的信息，有利于医生分析病变程度及治疗效果。
附图说明
图1为一幅颈动脉二维超声图像。
图2为利用内外膜轮廓得到斑块分割的ROI图像。
图3为用于斑块分割的ROI图像经过坐标转换的极坐标图像。
图4为极坐标图像中初始斑块外轮廓的检测结果示意图。
图5为内外膜围成的mask区域示意图。
图6为实施例1中斑块自动分割得到的精确轮廓结果示意图。
图7为本发明方法整体流程图。
图8为本发明血管分割方法流程图。
图9为颈动脉三维超声体数据中提取的当前帧图像。
图10为本发明血管自动跟踪过程示意图，图10（a）为当前帧劲动脉图像，图10（b）为下一帧图像。
图11（a）为包含ROI1曲线在内的矩形窗图像，图11（b）为图11（a）经过stick filter得到的图像。
图12为图11（b）经极坐标转换得到的极坐标图。
图13（a）为MAB初始轮廓示意图，图13（b）为包含MAB分割区域和MAB初始轮廓的矩形窗示意图，图13（c）为水平集演化得到的最终的MAB轮廓与手工勾勒的金标准比较的结果示意图，图13（d）为包含LIB分割区域和LIB初始轮廓的矩形窗示意图，图13（e）为水平集演化得到的最终的LIB轮廓（“.-”连线）与手工勾勒的金标准（虚线）比较的结果示意图。
图14为最大期望算法估计感兴趣区域1内灰度概率密度分布示意图。
具体实施方式
下面结合附图和实施例对本发明做进一步的详细描述。
本发明提供的超声图像中斑块的斑块自动分割算法，其实施步骤如下：
1 从颈动脉三维超声体数据提取当前帧图像，并读取图像对应的横、纵方向上单位像素的分辨率ρx和ρy，ρx为水平方向1mm对应的像素个数，ρy为垂直方向上1mm对应的像素个数。
2 对当前帧图像分割得到内外膜轮廓。
3 依据内外膜轮廓提取斑块分割的感兴趣区域（ROI），以减少后述分 割步骤的计算开销，如图2所示，其中“+”点为ROI的中心点，“-”线为步骤2中得到的内膜轮廓，“.”线为步骤2中得到的外膜轮廓，白色矩形框内的区域为ROI图像域；
3.1计算内膜轮廓围成区域的几何中心坐标(x0,y0)，见图2中“+”点；
3.2比较外膜的最大、最小横坐标与内膜中心横坐标的差值以及外膜最大、最小纵坐标与内膜中心纵坐标的差值获得最大差值L，L=max(|minx-x0|,|maxx-x0|,|miny-y0|,|maxy-y0|)，其中maxx,minx,maxy,miny分别为血管外膜的最大、最小横坐标和外膜最大、最小纵坐标。裁剪原图像得到矩形ROI区域，其横坐标范围为(x0-L-a,x0+L+a)，纵坐标范围为(y0-L-a,y0+L+a)，中心坐标为(x0,y0)，其中a为ROI比恰好能包含外膜的最小矩形多出的像素余量，一般取1~20个像素。若ROI区域坐标超过图像区域，则以图像区域作为界限，得到的ROI图像如图2所示。
4在ROI图像中自动检测得到初始斑块外轮廓线，其步骤包括：
4.1极坐标转换：
利用步骤3.1中提供的中心点(x0,y0)将ROI图像进行极坐标转换，并将内外膜轮廓点进行相应的极坐标转换，如图3所示，图中上下两条“.”线分别表示坐标转换后的血管内膜和外膜的轮廓；
4.2寻找径向极小值点与滤波：
将4.1中得到的极坐标图像中每列内外膜间的像素值进行曲线拟合，此实例中采用5阶多项式拟合，并在拟合得到的曲线上检测极小值，从极小值对应的像素点中选取离外膜轮廓最近的像素点。图4（a）显示了其中一列上灰度值曲线拟合以及极小值检测的过程，图中“.”为离散像素点，实线是经过拟合后的曲线，空心点为在拟合的曲线上检测的具有灰度极小值且离对应外膜轮廓点最近的点。
按照图4（a）所示的方法在分别每列内外膜间找到相应的极小值坐标， 得到一系列离散的点，见图4（b），图中“.”点为检测的极小值点。作为优化，对这些点的纵坐标进行中值滤波以保持坐标点之间的连续性，滤波后结果见图4（b）中的实线。
4.3将4.2中得到的滤波后线段坐标转换到直角坐标系中，并将这些点围成的闭合曲线作为初始斑块外边界C0。
5利用步骤4检测得到的初始斑块外边界C0进行水平集演化最终得到斑块轮廓，包括以下步骤：
5.1水平集初始化：
利用初始斑块外边界C0将ROI图像域上的水平集函数初始定义为符号距离函数：
outside(C0)表示斑块感兴趣区域ROI中初始斑块外边界以外的像素区域，x为斑块感兴趣区域ROI的像素点，|| ||为求欧式距离；
5.2水平集演化：
能量泛函定义为：Eφ=λs·Es+λl·El+λg·Eg+λd·Ed
其中，λs、λl、λg、λd分别为能量Es、El、Eg、Ed所占的权重。
为限制曲线长度的能量，保证演化曲线的光滑性，λs越大，曲线越平滑。其中， 为正则化的Heaviside函数，0<ε<0.001，Ω为斑块感兴趣区域ROI。
E l = ∫ Ω δ ( φ ( x ) ) ∫ Ω B L ( x , x ^ ) [ H ( φ ( x ^ ) ) ( I ( x ^ ) - u x ) 2 + ( 1 - H ( φ ( x ^ ) ) ) ( I ( x ^ ) - v x ) 2 ] d x ^ dx ]]>为局部Chan-Vese的能量，其作用是使曲线内外在一定邻域内灰度统计特征差异尽 可能大。其中， 为正则化的Dirac函数， 定义一个半径为rL的圆形邻域，I()为像素灰度值。
u x = ∫ Ω B L ( x , x ^ ) H ( φ ( x ^ ) ) I ( x ^ ) d x ^ / ∫ Ω B L ( x , x ^ ) H ( φ ( x ^ ) ) d x ^ ]]>为像素点x邻域范围内位于迭代演化轮廓C内部的像素点的灰度均值。
v x = ∫ Ω B L ( x , x ^ ) ( 1 - H ( φ ( x ^ ) ) ) I ( x ^ ) d x ^ / ∫ Ω B L ( x , x ^ ) ( 1 - H ( φ ( x ^ ) ) ) d x ^ ]]>为像素点x邻域范围内位于迭代演化轮廓C外部的像素点的灰度均值，上述两式中 为ROI区域的像素点。
E g = ∫ Ω δ ( φ ( x ) ) ∫ mask 1 N [ H ( φ ( x ^ ) ) Σ i = 1 N | u 0 i ( x 3 ) - c + i | 2 + ( 1 - H ( φ ( x ^ ) ) ) Σ i = 1 N | u 0 i ( x ) - c - i | 2 ] dx ]]>为Gabor滤波器能量，表示对于外膜轮廓CMAB和内膜轮廓CLIB围成的环状区域内的像素点，曲线内外两侧经过Gabor滤波后特征统计值的差异，CMAB与CLIB围成的mask区域如图5所示。其中，N指Gabor滤波器组中滤波器的个数， （i＝1,2…N）为ROI图像经过Gabor第i个滤波器滤波后的输出。
为mask区域中且位于迭代演化轮廓C内部的像素点经过第i个Gabor滤波得到的响应均值。
为mask区域中且位于迭代演化轮廓C外部的像素点经过第i个Gabor滤波得到的响应均值。
Ed=∫Ωδ(φ(x))·Bd(x)dx为距离限制能量，起到限制曲线与外膜轮廓CMAB的距离的作用。其中，
msd=||CMAB,CLIB||为外膜轮廓CMAB和内膜轮廓CLIB之间的欧式距离，d为内外膜 间的最小距离，根据统计值设为0.4~0.5mm之间，本实例中采用0.4mm，当曲线与外膜距离小于d时该项提供收缩力，而当距离大于msd时，提供与距离成正比的扩张力，F为扩张力与距离的比例因子，一般0<F≤1。
为了得到使得能量Eφ最小的解φ，采用变分法得到Euler-Lagrange方程
∂ φ ( x ) ∂ t = - ∂ E φ ∂ φ ( x ) ]]>
计算 得到水平集φ演化的偏微分方程:
∂ φ ( x ) ∂ t = δ ( φ ( x ) ) · { λ s div ( ▿ φ ( x ) | ▿ φ ( x ) | ) + λ l ∫ Ω B L ( x , x ^ ) δ ( φ ( x ^ ) ) [ ( I ( x ^ ) - u x ) 2 - ( I ( x ^ ) - v x ) 2 ] d x ^ + λ g · 1 N · Σ i = 1 N [ ( u 0 i ( x ) - c + i ) 2 - ( u 0 i ( x ) - c - i ) 2 ] + λ d · B d ( x ) } ]]>其中，第一项为曲率项，控制曲线的光滑性；第二项为局部Chan-Vese项，控制一定邻域内曲线内外的灰度统计值差异；第三项为全局Gabor项，使在内外膜围成的mask区域内，曲线内外的Gabor输出值差异性增大；第四项为距离限制项，保证斑块的外轮廓与外膜的距离控制在合理的范围内。
各能量项的权重取值依不同的图像特点变化，0<λs≤20，0≤λl≤20，0≤λg≤20，0<λd≤10，在本实例中，λs=6,λl=4,λg=2,λd=1；
局部Chan-Vese项的半径2≤rL≤20，本实例中取6个像素；
Gabor滤波器组一般选择4~8个角度，本例中采用 四个角度；
高斯窗函数1≤σ≤20，尺度因子1≤scale≤5，频率 本实例中 σ = 3 π 2 , ]]>scale=3， f = 2 . ]]>
不断迭代更新φ(x)：
φ z ( x ) = φ z - 1 ( x ) + δ ( φ z - 1 ( x ) ) · { λ s div ( ▿ φ z - 1 ( x ) | ▿ φ z - 1 ( x ) | ) + λ l ∫ Ω B L ( x , x ^ ) δ ( φ z - 1 ( x ^ ) ) [ ( I ( x ^ ) - u x ) 2 - ( I ( x ^ ) - v x ) 2 ] d x ^ + λ g · 1 N · Σ i = 1 N [ ( u 0 i ( x ) - c + i ) 2 - ( u 0 i ( x ) - c - i ) 2 ] + λ d · B d ( x ) } ]]>
依据更新得到的φz(x)计算本次迭代的演化轮廓Cz={x|φz(x)=0}；
本次迭代演化轮廓Cz与前一次迭代得到的演化轮廓Cz-1相同或者z达到迭代次数上限，则本次迭代得到的演化轮廓即为最终斑块外边界Cplaque。
本发明斑块外边界演化迭代次数上限一般为5000次以上，本实例中迭代次数上限设为10000次。
5.3取5.2中得到的最终斑块外边界Cplaque与内膜轮廓CLIB之间的区域作为斑块区域，或者对该区域进行形态学开操作得到斑块区域，如图6所示。开操作可采用的结构元素有圆盘结构、矩形结构等，本发明实例中采用圆盘结构，圆盘半径取4。
利用上述步骤对三维超声颈动脉体数据中每帧图像进行分割，得到的斑块区域可以通过体积和面积的测量作为粥样硬化测量指标，从而辅助医生进行病理分析和药物治疗效果观察评价，下面对其详细说明：
测量指标的分析：
根据前面步骤分割得到的每一帧超声图像中颈动脉斑块的轮廓计算以下指标：
TPA：对体数据中任意一帧颈动脉超声图像，统计其上自动分割得到的斑块轮廓Cp(k)（k=1,2，...,Nr,Nr为体数据中颈动脉图像的帧数）包围的区域内像素点的个数Np(k)，根据像素点个数与面积之间的比例关系ρ（ρ=ρx·ρy为每mm2中像素点的个数）计算得到这一帧上的斑块面积，本实施例中ρ=1041，对体数据中每一帧图像重复上述操作，得到每一帧图片上的TPA；
TPV：根据每一帧图片上的斑块面积Ap(k),（k=1,2，...,N，N为体数据的总帧数），求得相邻两帧之间空间上的体积（如第k帧和第k+1帧之间包含的体积Vp(k)，计算公式为：Vp(k)=0.5*(Ap(k)+Ap(k+1))*d，其中d为相邻两帧之间的距离，本实施例中为1mm，将每两帧之间的斑块体积求和，得到最终的TPV；
本实施例中的指标结果见表1，其中k为体数据中帧数的序号，k=1表示离颈动脉分岔处最远的起始帧。
表1测量指标的分析结果
所述步骤2中内外膜轮廓的分割可采用半自动混合水平集方法（具体实施方法参见E.Ukwatta,et al,Three-dimensional ultrasound of carotid atherosclerosis:Semi-automated segmentation using a level set-based method,Medical Physics,2011,38(5):p.2479-2493）、边缘检测及形态学处理算法（具体实施方法参见Xin Yang,Mingyue Ding,et al,Common Carotid Artery Lumen Segmentation in B-mode Ultrasound Transverse View Images,I.J.Image,Graphics and Signal Processing,2011:p.15-21）等，本发明提供一种如下的分割方法：
（1）确定基准轮廓：
若当前帧图像为的第一帧图像，则凭经验判断出外轮廓的大致位置，并勾取几个比较明显的位于血管外轮廓上的基准点，再通过插值的方法形成一条闭合曲线作为基准轮廓，如图9中虚线构成的轮廓所示；否则，将上一帧得到的备用基准轮廓作为当前帧的基准轮廓；
（2）利用形态学膨胀法将基准轮廓外推形成第一感兴趣区域ROI1；
将基准轮廓利用形态学膨胀方法向外推移一定距离，可取10到20个像素，本实施例中为10个像素，形成感兴趣区域1（ROI1），见图9中的实心闭合轮廓线；
（3）在第一感兴趣区域ROI1中检测颈动脉血管外轮廓CMAB：
（3.1）在感兴趣区域1（ROI1）中检测MAB的初始轮廓线
本步骤的技术思路是：在当前帧图像中提取包含第一感兴趣区域ROI1的矩形窗口图像，将其转化到极坐标图像，在极坐标图像上每隔相同角度在径向上靠近基准轮廓的像素点中搜索梯度值最大者作为初始轮廓点。所述间隔角度可任意选取，一般选5到20度；在搜索中可将搜索范围定为与基准轮廓间隔2~10个像素距离的像素点集。具体过程如下：
（3.1.1）计算ROI1区域的几何中心坐标，自动生成包含ROI1在内的矩形窗口图像，如图11（a）所示，图上的轮廓为经过坐标转化后显示在矩形窗口图像中的ROI1的边界；
（3.1.2）作为优化，对矩形窗口图像进行滤波，以消除部分噪声和增强颈动脉血管与周围组织的对比度，如图11（b）。此处可采用滤波的方法有：棒滤波（stick filter）、双边滤波（bilateral filter）、各向异性扩散滤波（SRAD filter）等，本实例采用棒滤波，其具体实现方法参见S.D.Pathak,V.Chalana,D.R.Haynor,and Y.Kim,“Edge-guided boundary delineation in prostate ultrasound images,”IEEE Trans.Med Imaging 2000，19(12),1211–1219；
（3.1.3）将矩形窗口图像转化为极坐标图像，将极坐标图像沿径向均匀划分为N个区间，见图12（a），N的取值在18至72之间，在本实施例中N取24；取图像12（a）中每个区间径向上的中线，计算每条中线上梯度值最大的像素点，标记为“+”，如图12（b）所示。每个标记点的坐标表示为{Mi|ri,θi=(i-0.5)*2/N,i＝1,2,...,N}，将这些点作为极坐标图上外轮廓MAB的初始轮廓点；
（3.1.4）将步骤（3.1.3）中的标记点从极坐标还原到直角坐标系中，平移后得到在原图中的坐标位置，并用细线依次连接起来，形成MAB的初始轮廓 如图13（a）所示；
（3.2）利用混合分布估计第一感兴趣区域ROI1的灰度概率密度分布， 从而得到混合分布的shape参数Kj和scale参数θj，j=1,…,M，M为混合分布中单项分布的类别数；
本发明采用图像的灰度概率密度分布作为水平集模型的特征，估计感兴趣区域1（ROI1）的灰度概率密度分布。本发明将超声图像的灰度概率密度分布近似估计为混合gamma分布（mixture gamma distribution），具体估计步骤如下：
（3.2.1）对感兴趣区域1（ROI1）内的像素点进行随机采样，样本量一般取800到2000个像素，在本实施例中定为2000；
（3.2.2）对采样得到的样本点进行聚类，以便将样本点进行大致分类，聚类方法可采用K均值聚类（K-means），或模糊聚类（Fuzzy）等，在本实施例中采用了K均值聚类，且聚类类别设为3（聚类类别一般根据目标所包含的组织种类来设定），为减少计算开销，本实例聚类迭代次数设为50，将样本点按照灰度均值分为3类，K均值聚类方法具体实现参见J.A.Hartigan,et al,A K-Means Clustering Algorithm,Journal of the Royal Statistical Society,197928(1),pp:100-108；
（3.2.3）分别在聚类得到的三类样本点上，采用最大似然估计法（Maximum Likelihood）对样本点进行gamma分布估计，得到三组gamma分布的参数(K1,θ1),(K2,θ2),(K3,θ3),K,θ分别对应gamma分布中的shape参数和scale参数。自动计算三类样本点在总样本点中所占比例Wi: 其中Ni表示样本中属于第i个分量的像素个数，i＝1,2,3.Ns代表样本点的总个数。最大似然估计法的实现具体参见Aldrich,John,.A.Fisher，The making of maximum likelihood,Statistical Science,1997,12(3):162-176；
（3.2.4）将步骤（3.2.3）中的三组参数(W1,K1,θ1),(W2,K2,θ2),(W3,K3,θ3)作为期望最大化（Expectation Maximization）算法的初始值代入到期望最大化算法中，迭代得到最终的三组参数(W1,K1,θ1),(W2,K2,θ2),(W3,K3,θ3)作为混合 gamma分布（mixture gamma distribution）的最终参数，这样ROI1内的灰度密度概率分布可近似表示为： 其中 G ( I x | K i , θ i ) = I x K i - 1 · exp ( - I x θ i ) / ( θ i K i · ( K i - 1 ) ! ) ]]>为gamma分布的函数通式，Ix为灰度级。
期望最大化算法的具体实现参见A.Dempster,N.Laird,and D.Rubin,“Maximum likelihood from incomplete data via the EM algorithm,”J.R.Stat.Soc.(Ser.B),1977,39:p.1-38。本实例的直方图分布及用期望最大化算法估计得到的gamma混合分布曲线如图14所示，其中三条实心曲线为混合gamma分布（mixture gamma distribution）的3个不同分量，虚线为3个不同分量组合得到的混合gamma分布曲线，表2为期望最大化算法得到的三个gamma分布分量的参数值。
表2gamma混合分布的参数估计结果
i Wi Ki 0i
1 0.7012 4.1467 2.9920
2 0.2177 8.7697 3.4985
3 0.0811 7.5027 9.0800
本发明实例的灰度概率密度分布近似估计为混合gamma分布，除此之外，还可采用混合gauss分布、混合rayleigh分布等混合分布估计方法。
（3.3）初始化迭代次数t=1以及符号距离函数
利用MAB的初始轮廓 将ROI1图像域上的水平集函数初始定义为符号距离函数
表示MAB的初始轮廓， 表示ROI1图像域中MAB的初始轮廓C1以外的像素区域，x1为ROI1 图像域中的像素点，|| ||为求欧式距离；
（3.4）计算
为减少计算开销，本方法自动生成一个包含MAB分割区域的矩形窗口图像，具体方法为：将基准轮廓采用形态学腐蚀法内推一段距离形成一条新的轮廓，此距离可根据图像大小恰当选择，本实施例中定为18个像素。形成的新轮廓与感兴趣区域1（ROI1）边界（外边界）共同构成MAB的分割区域，水平集演化过程中演化轮廓仅在此分割区域内演化。MAB分割区域如图13（b）中非黑色区域所示，区域内的轮廓为MAB的初始轮廓；
MAB分割的能量泛函定义为：
其中，λ1为 所占的权重，ρ1为 的权重。
为限制曲线长度的能量，保证演化曲线的光滑性，λ1越大，曲线越光滑，其中 为正则化的Heaviside函数，ε为一个极小的正实数，通常取0<ε<0.001，Ω1为ROI1图像域， 为梯度算子。
E pd f 1 = - ∫ Ω 1 H ( φ 1 ( x 1 ) ) · ln ( Σ j M q Ai ( j ) · G ( I x 1 | K j , θ j ) ) dx 1 - ∫ Ω 1 ( 1 - H ( φ 1 ( x 1 ) ) ) · ln ( Σ j M q Ae ( j ) · G ( I x 1 | K j , θ j ) ) dx 1 ]]>
为水平集区域信息能量项，它利用最大后验概率准则（Maximization a Posteriori）的思想，按照灰度概率密度分布来划分像素点的类别。其中，M为混合gamma分布的分量个数，qAi(j),qAe(j)分别为演化轮廓内外区域每种gamma分布的权重， 为参数为Kj,θj的gamma分布函数, 是像素点x1的灰度值，取值范围为0到255。 表示整幅图像灰度概率密度分布的似然性， 越小，目标和背景所包括的像素分离得越好。
为了得到使能量 达到最小的解φ1(x1)，采用变分法得到Euler-Lagrange 方程:
∂ φ 1 t ( x 1 ) ∂ t = - ∂ E φ 1 ∂ φ 1 t ( x 1 ) ]]>
计算 得到水平集φ1(x1)演化的偏微分方程:
∂ φ 1 ( x 1 ) ∂ t = - ρ 1 · δ ( φ 1 ( x 1 ) ) · [ ln ( Σ j M q Ai ( j ) · G ( I x 1 | K j , θ j ) - ln ( Σ j M q Ae ( j ) · G ( I x 1 | K j , θ j ) ] + λ 1 · δ ( φ 1 ( x 1 ) ) · div ( ▿ φ 1 ( x 1 ) | ▿ φ 1 ( x 1 ) | ) ]]>
上式中，δ(x)为正则化的Dirac函数，定义为：
其中ε为一个很小的正实数，通常取0<ε<0.001。
计算φ1t(x1)：
φ 1 t ( x 1 ) = φ 1 t - 1 ( x 1 ) - ρ 1 · δ ( φ 1 t - 1 ( x 1 ) ) · [ ln ( Σ j M q Ai ( j ) · G ( I x 1 | K j , θ j ) - ln ( Σ j M q Ae ( j ) · G ( I x 1 | K j , θ j ) ] + λ 1 · δ ( φ 1 t - 1 ( x 1 ) ) · div ( ▿ φ 1 t - 1 ( x 1 ) | ▿ φ 1 t - 1 ( x 1 ) | ) ]]>
（3.5）计算本次迭代演化外轮廓
C 1 t = { x 1 | φ 1 t ( x 1 ) = 0 } ]]>
（3.6）若本次迭代演化外轮廓 与前一次迭代得到的演化外轮廓 相同或者t达到迭代次数上限，则本次迭代得到的演化轮廓即为颈动脉血管外轮廓CMAB，迭代结束；否则，t=t+1，返回步骤（3.4）；
上式中，理论上ρ1>0和λ1>0即可，根据实际经验值可将ρ1取0到2之间，λ1取值范围为1到10之间，在本实例中迭代次数上限为100次，ρ1=0.2，λ1=1.2，精确的MAB边界轮廓线CMAB见图13（c）中“.”构成轮廓所示。
（4）在当前帧图像中将颈动脉血管外轮廓CMAB围成的区域作为第二感兴趣区域ROI2；
（5）在第二感兴趣区域ROI2中检测颈动脉血管内轮廓CLIB：
（5.1）利用形态学腐蚀法将颈动脉血管外轮廓CMAB内推，将内推得到的轮廓线作为初始颈动脉血管内轮廓
将得到的MAB轮廓CMAB围成的区域作为LIB轮廓分割的感兴趣区域2（ROI2），并利用形态学腐蚀法将MAB轮廓内推一段距离得到的轮廓线作为LIB的初始轮廓 如图13（d）中的闭合轮廓所示，本实施例中内推距离取18个像素。
（5.2）初始化迭代次数t′＝1以及符号距离函数
表示第二感兴趣区域ROI2中初始颈动脉血管内轮廓 以外的像素区域，x2为第二感兴趣区域ROI2中的像素点；
（5.3）计算
定义LIB分割的水平集能量为：
其中 为结合灰度概率密度分布的区域信息能量项，用以估计图像每一点属于目标和背景的概率， 为曲率项，用以控制曲线的光滑性，这两项能量项的定义与步骤（3.4）相同；ρ2，λ2，β2分别为三项能量的权重；
ET为MAB与LIB之间的距离保持项，定义为：
其中Ω2为ROI2图像域；
D(x2,y)为LIB演化轮廓上的点x2与MAB最终轮廓上的点y之间的欧氏距 离，dT为一般情况下MAB和LIB之间相隔的最小距离（经验值）。当点x2与MAB轮廓的距离小于dT时，就会受到一个排斥力BT(x2)；当点x2与MAB轮廓的距离大于dT时，排斥力为0。
为了得到使得能量 最小的解φ2(x2)，采用变分法得到Euler-Lagrange方程
∂ φ 2 t ′ ( x 2 ) ∂ t = - ∂ E φ 2 ∂ φ 2 t ′ ( x 2 ) ]]>
计算 得到水平集φ2(x2)演化的偏微分方程:
∂ φ 2 ( x 2 ) ∂ t = - ρ 2 · δ ( φ 2 ( x 2 ) ) · [ ln ( Σ j M q Li ( j ) · G ( I x 2 | K j , θ j ) - ln ( Σ j M q Le ( j ) · G ( I x 2 | K j , θ j ) ] + λ 2 · δ ( φ 2 ( x 2 ) ) · div ( ▿ φ 1 ( x 2 ) | ▿ φ 1 ( x 2 ) | ) + β 2 · δ ( φ 2 ( x 2 ) ) · B T ( x 2 ) ]]>
上式中，δ(x)为正则化的Dirac函数，定义为：
其中ε为一个很小的正实数，通常取0<ε<0.001。
更新φ2t′(x2)：
φ 2 t ′ ( x 2 ) = φ 2 t ′ - 1 ( x 2 ) - ρ 2 · δ ( φ 2 t ′ - 1 ( x 2 ) ) · [ ln ( Σ j M q Li ( j ) · G ( I x 2 | K j , θ j ) - ln ( Σ j M q Le ( j ) · G ( I x 2 | K j , θ j ) ] + λ 2 · δ ( φ 2 t ′ - 1 ( x 2 ) ) · div ( ▿ φ 1 t ′ - 1 ( x 2 ) | ▿ φ 1 t ′ - 1 ( x 2 ) | ) + β 2 · δ ( φ 2 t ′ - 1 ( x 2 ) ) · B T ( x 2 ) ]]>
（5.4）计算本次迭代演化内轮廓
C 2 t ′ = { x 2 | φ 2 t ′ ( x 2 ) = 0 } ]]>
（5.5）若本次迭代演化内轮廓 与前一次迭代得到的演化内轮廓 相同或者达到迭代次数上限，则本次迭代得到的演化轮廓为颈动脉血管内轮廓CLIB，迭代结束；否则，t′=t′+1，返回步骤（5.3）；
ρ2和λ2取值范围与MAB分割中一致，β2取值范围较大，一般取5以上，dT取14到20之间。
在本实例中迭代次数上限为8000次，ρ2=0.15，λ2=5，β2=20，最小距离dT=18个像素，得到最终的LIB边界轮廓线CLIB，见图13（e）的由“·-”构成的轮廓，虚线构成的轮廓为用于对比的手工勾勒的LIB轮廓。
（6）跟踪确定下一帧图像的备用基准轮廓：
（6.1）确定跟踪区域Ω′；将当前帧颈动脉MAB精确轮廓CMAB通过形态学膨胀外推一定距离（一般取10到30个像素距离，本实施例中为20个像素）与之构成一块区域Ω′，见图10（a）中两条轮廓所围成的区域，其中较小轮廓代表分割得到的MAB轮廓，较大轮廓代表MAB轮廓经过外推得到的轮廓。此区域可视为包含颈动脉外膜及其外周组织的区域，一般情况下每一帧颈动脉图像都有这样一个区域，本方法根据这个区域的偏转移动来实现血管的自动跟踪；
（6.2）计算当前帧与下一帧之间的血管跟踪映射关系，具体方法为：
（6.2.1）对区域Ω′中任意一点(x3,y3)，假设其在下一帧图像上通过映射M:{(x4,y4)=τ+D*(x3,y3)}移动到了点(x4,y4)，其中τ为平移向量，D为形变矩阵。本方法中认为连续帧的血管发生的形变很小，因此本实施例中，
D = 1 0 0 1 ]]>
假设τ=(τx,τy)，其中τx为点(x3,y3)于映射过程中在x方向平移的距离，τy为点(x3,y3)于映射过程中在y方向平移的距离则对应点之间的映射关系可以表示为：
x 4 y 4 = τ x τ y + 1 0 0 1 x 3 y 3 ]]>
并且本方法认为此区域中每个点均满足同样的映射，通过此映射当前帧中的颈动脉外周组织区域投影到了下一帧的新的区域；
（6.2.2）计算当前帧图像的区域Ω′中点(x3,y3)和其在下一帧图像上的对应点(x4,y4)的灰度值I(x3,y3)和J(x4,y4)之间差异的平方：
ΔI2=[I(x3,y3)-J(x4,y4)]2
对这个区域内每个点进行上述操作，将每个点上的灰度差异求和：
S=∑Ω′ΔI2
（6.2.3）分别在下一帧图片的横轴、纵轴方向上一定范围内，搜索到使S达到最小的区域，此区域对应的映射关系即为所求的连续两帧之间的血管跟踪映射关系M，本实施例中搜索范围在横轴纵轴上均取为正负20个像素；
(6.3)将当前帧图片的MAB精确轮廓CMAB按照此映射关系M投影到下一帧上，构成下一帧颈动脉图像的备用基准轮廓，见图10（b）中的虚线轮廓。
与以往的颈动脉血管内轮廓和外轮廓的分割算法相比，本发明提供的超声图像中颈动脉血管内外膜自动分割算法与现有方法有几点区别：1.仅在每组颈动脉体数据的第一帧上采用极少的经验干预，不需在每一帧上都采用手工勾点；2.每组颈动脉体数据的每一帧图像均采用自动检测得到内外膜的初始轮廓线；3.采用了灰度概率密度分布作为水平集中应用到的特征，能够有效克服横截面方向上探头位置造成的伪影和弱边界情况；4.利用当前帧分割得到的结果预测下一帧目标的大概位置，减少了不必要的计算开销。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。