CN102800087A - 超声颈动脉血管膜的自动分割方法 - Google Patents

Abstract

本发明属于计算机技术与医学图像的交叉领域，具体涉及一种超声图像颈动脉血管横截面方向上的内外轮廓自动分割方法。具体步骤如下：选取当前帧图像；人工定义基准轮廓点并插值形成基准轮廓；将基准轮廓外推形成第一感兴趣区域；自动检测颈动脉轮廓：包括在第一感兴趣区域转化得到的极坐标图像上寻找径向梯度最大值点作为初始外轮廓点、利用混合分布估计第一感兴趣区域内的灰度概率密度分布、水平集演化得到外轮廓、将外轮廓内推形成初始内轮廓、水平集演化得到内轮廓。本发明能对超声颈动脉图像进行有效的分割，减少医生工作量；同时基于本方法分割得到的颈动脉内外轮廓可用于计算颈动脉血管壁体积，为颈动脉粥样硬化病变研究提供准确有效的信息。

Description

超声颈动脉血管膜的自动分割方法

技术领域

本发明属于计算机技术与医学图像的交叉领域，具体涉及到一种超声图像中颈动脉血管横截面方向上的内外膜的自动分割方法。

背景技术

心脑血管疾病是危害人类健康和造成死亡的重要疾病，如心肌梗塞、中风等，其主要病因是颈动脉粥样硬化。与传统的内中膜厚度（Intima-MediaThickness,IMT）度量指标相比，血管壁体积（Vessel Wall Volume,VWV）等度量指标更能准确直观地反映动脉粥样硬化状况，因此成为近年来预测心脑血管疾病风险的重要指标之一。

血管壁体积在超声颈动脉图像上可根据血管腔-内膜边界（Lumen-IntimaBoundary，LIB）和中膜-外膜边界（Media-Adventitia Boundary，MAB）测量得到。在横截面方向上，正常的颈动脉血管腔在超声图像中表现为低回声的均质区域，外膜表现为高回声的非均质带状区域。图1显示了超声成像下的颈动脉横截面及其结构示意图。图1中实心点围成的轮廓是颈动脉内轮廓LIB边界，其内部较为均匀的黑色低回声部分是血管腔（lumen），其外靠近内轮廓边界LIB的一层带状区域是颈动脉内膜；十字形点围成的闭合轮廓线是颈动脉外轮廓边界（MAB），其外部包含高回声的不均匀的区域是颈动脉外膜。斑块为图中实线所围成的区域，其内部灰度分布不均匀，且与周围的灰度特性相近，不易区分。

传统的颈动脉血管膜及斑块的边界主要通过操作者（医生）肉眼观察手工勾勒获得。其主要缺陷在于它很大程度上取决于操作者的主观判断及操作经验，因此在不同操作者对同一个目标边界的判断，甚至同一个操作者在不同的时间点对同一目标边界的判断也会有差异。此外，手工勾勒边界非常耗时，通常经验丰富的操作者勾勒一组超声颈动脉图像的内外膜边界大致需要10分钟，无法满足临床病例分析的需求。

颈动脉血管膜的半自动和自动分割方法是近十年来的研究热点。颈动脉血管膜的分割算法大致分为两类：一是边缘检测与形态学算法结合的算法，如Ali K.Hamou，Mahmoud R.El-Sakka，A novel segmentation techniquefor carotid ultrasound images.ICASSP,2004:p.521-524;Xin Yang,MingyueDing,et al,Common Carotid Artery Lumen Segmentation in B-modeUltrasound Transverse View Images,I.J.Image,Graphics and SignalProcessing,2011:p.15-21等。这种方法首先通过一些简单的预处理对超声颈动脉图像进行对比度增强，使目标和背景之间的差异更加明显，然后采用边缘检测算子提取出目标轮廓，最后用形态学算子填补轮廓线的缺陷或消除多余的轮廓线。此方法缺点在于精度不高，对梯度值不高的边缘检测能力弱，仅适用于目标边界梯度值较大的情况。另一种分割方法是活动轮廓模型算法，如J.D.Gill,et al,Accuracy and variability assessment of asemiautomatic technique for segmentation of the carotid arteries fromthree-dimensional ultrasound images,Medical Physics,2000,27(6):p.1333-1342;Abir Zahalka,Aaron Fenster,An automated segmentation method forthree-dimensional carotid ultrasound images,Physics in Medicine and Biology,2001,46:p.1321-1342;Ali K.Hamou,et al,Carotid Ultrasound SegmentationUsing DP Active Contours,ICIAR,2007:p.961-971;E.Ukwatta,et al,Three-dimensional ultrasound of carotid atherosclerosis:Semi-automatedsegmentation using a level set-based method,Medical Physics,2011,38(5):p.2479-2493.这种算法需要先给出目标轮廓的初始轮廓，使其在定义的能量项的作用下朝着目标轮廓演化。以往此类方法为：操作者首先给出一条边界的初始轮廓，使其在对应于该边界定义的能量项下迭代演化得到最终的分割结果；然后以此分割结果结合一定的先验知识获得第二条边界的初始轮廓，并在对应于第二条边界所定义的能量项下演化，得到最终的轮廓。这种方法需要操作者手工在每一帧图像上给定初始轮廓，或给出较多的初始点以构成初始轮廓，人工干预较多，仍有较大的改进空间。

发明内容

本发明的目的在于提供一种超声图像中颈动脉血管内外膜的自动分割方法，能够有效分割出颈动脉内外膜，为颈动脉粥样硬化的病变观察和药物治疗提供分析参数。

超声颈动脉血管膜的自动分割方法，包括以下步骤：

（1）若当前帧图像为颈动脉三维超声体数据的第一帧图像，则在当前帧图像上选择明显位于血管外轮廓上的像素点集作为基准点，通过基准点间插值形成闭合曲线，将其作为当前帧的基准轮廓；否则，将上一帧得到的备用基准轮廓作为当前帧的基准轮廓；

（2）利用形态学膨胀法将基准轮廓外推形成第一感兴趣区域ROI1；

（3）在第一感兴趣区域ROI1中检测颈动脉血管外轮廓C_MAB：

（3.1）在当前帧图像中提取包含第一感兴趣区域ROI1的矩形窗口图像，将其转化到极坐标图像，在极坐标图像上每隔相同角度在径向方向靠近基准轮廓的像素点中搜索梯度值最大者作为初始轮廓点，将初始轮廓点转化到直角坐标系并连线形成初始颈动脉血管外轮廓

（3.2）利用混合分布估计第一感兴趣区域ROI1的灰度概率密度分布，从而得到混合分布的shape参数K_j和scale参数θ_j，j＝1,…,M，M为混合分布中单项分布的类别数；

（3.3）初始化迭代次数t=1以及符号距离函数

表示第一感兴趣区域ROI1中初始颈动脉血管外轮廓

以外的像素区域，x₁为第一感兴趣区域ROI1中的像素点，||  ||为求欧式距离；

（3.4）计算

${\mathrm{\φ}}_1^t\left(x_1\right)={\mathrm{\φ}}_1^{t-1}\left(x_1\right)-{\mathrm{\ρ}}_1\·\mathrm{\δ}\left({\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right)\·\left[\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ai}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)-\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ae}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)\right]$

$+{\mathrm{\λ}}_1\·\mathrm{\δ}\left({\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right)\·\mathrm{div}\left(\frac{{\▿\mathrm{\φ}}_1^{t-1}\left(x_1\right)}{\left|{\▿\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right|}\right)$

其中，

这里0＜ε＜0.001，0＜ρ₁≤10，0＜λ₁≤20，q_Aj(j)和q_Ae(j)分别为前一次迭代演化轮廓

的内、外区域中第j项单项分布的权重，为参数为K_j,θ_j的第j项单项分布函数，

为像素点x₁的灰度值；

（3.5）计算本次迭代演化外轮廓

$C_1^t=\left\{x_1\right|{\mathrm{\φ}}_1^t\left(x_1\right)=0\}$

（3.6）若本次迭代演化外轮廓与前一次迭代得到的演化外轮廓相同或者t达到迭代次数上限，则本次迭代得到的演化轮廓即为颈动脉血管外轮廓C_MAB，迭代结束，进入步骤（4）；否则，t＝t+1，返回步骤（3.4）；

（4）在当前帧图像中将颈动脉血管外轮廓C_MAB围成的区域作为第二感兴趣区域ROI2；

（5）在第二感兴趣区域ROI2中检测颈动脉血管内轮廓C_LIB：

（5.1）利用形态学腐蚀法将颈动脉血管外轮廓C_MAB内推，将内推得到的轮廓线作为初始颈动脉血管内轮廓

（5.2）初始化迭代次数t′＝1以及符号距离函数

表示第二感兴趣区域ROI2中初始颈动脉血管内轮廓

以外的像素区域，x₂为第二感兴趣区域ROI2中的像素点；

（5.3）计算

${\mathrm{\φ}}_2^{t^{\′}}\left(x_2\right)={\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)-{\mathrm{\ρ}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·\left[\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Li}}\left(j\right)\·G\left(I_{x_2}\right|K_j,{\mathrm{\θ}}_j)-\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Le}}\left(j\right)\·G\left(I_{x_2}\right|K_j,{\mathrm{\θ}}_j)\right]$

$+{\mathrm{\λ}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·\mathrm{div}\left(\frac{{\▿\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)}{{|\▿\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)|}\right)+{\mathrm{\β}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·B_T\left(x_2\right)$

其中，0＜ρ₂≤10，0＜λ₂≤20，0＜β₂≤100，q_Li(j)和q_Le(j)分别为前一次演化轮廓

的内、外区域中第j项单项分布的权重，

为参数为K_j,θ_j的第j项单项分布函数，为像素点x₂的灰度值；

（5.4）计算本次迭代演化内轮廓

$C_2^{t^{\′}}=\left\{x_2\right|{\mathrm{\φ}}_2^{t^{\′}}\left(x_2\right)=0\}$

（5.5）若本次迭代演化内轮廓

与前一次迭代得到的演化内轮廓

相同或者达到迭代次数上限，则本次迭代得到的演化轮廓为颈动脉血管内轮廓C_LIB，迭代结束，进入步骤（6）；否则，t′＝t′+1，返回步骤（5.3）；

（6）跟踪确定下一帧图像的备用基准轮廓：

（6.1）在当前帧图像中，通过形态学膨胀法将颈动脉血管外轮廓C_MAB外推，将外推得到的闭合轮廓线与C_MAB构成的区域作为跟踪区域Ω；

（6.2）在下一帧图像中搜索使得灰度差异和S＝∑_ΩΔI²最小的像素点集(x₄,y₄)，其中ΔI²＝[I(x₃,y₃)-J(x₄,y₄)]²，I(x₃,y₃)为当前帧图像的跟踪区域Ω中像素点(x₃,y₃)的灰度值，J(x₄,y₄)为下一帧图像的像素点(x₄,y₄)的灰度值；

（6.3）令当前帧图像的跟踪区域中的点(x₃,y₃)与下一帧图像的像素点(x₄,y₄)的映射关系为 $\left(\begin{array}{c}{x}_4\\ y_4\end{array}\right)=\left(\begin{array}{c}{\mathrm{\τ}}_x\\ {\mathrm{\τ}}_y\end{array}\right)+\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}{x}_3\\ y_3\end{array}\right),$ 结合像素点(x₃,y₃)与(x₄,y₄)求解映射关系中的参数τ＝(τ_x,τ_y)，从而确定映射关系；

（6.4）根据步骤（6.3）建立的映射关系将当前帧图像的颈动脉血管外轮廓C_MAB映射为下一帧图像的备用基准轮廓。

所述混合分布为混合gamma分布或混合gauss分布或混合rayleigh分布。

本发明的技术效果体现在：

与以往的颈动脉血管内轮廓和外轮廓的分割算法相比，本发明提供的超声图像中颈动脉血管内外轮廓自动分割算法与现有方法有几点区别：1.仅在每组颈动脉体数据的第一帧上采用极少的经验干预，不需在每一帧上都采用手工勾点；2.每组颈动脉体数据的每一帧图像均采用自动检测得到内外膜的初始轮廓线；3.采用了灰度概率密度分布作为水平集中应用到的特征，能够有效克服横截面方向上探头位置造成的伪影和弱边界情况；4.利用当前帧分割得到的结果预测下一帧目标的大概位置，减少了不必要的计算开销。

本发明提供的超声图像中颈动脉血管内外轮廓的计算机自动分割算法和指标测量方法能够达到以下目标：（1）能有效应对超声图像中存在的噪声及弱边界情况；（2）能较为精确地分割出血管内外膜；（3）能较大幅度地减少医生的工作量；（4）基于本方法得到的度量指标能够提供准确有效的信息，有利于医生分析病变程度及治疗效果。

附图说明

图1为一幅颈动脉二维超声图像。

图2为一组颈动脉三维超声图像。

图3为图2颈动脉三维超声图像中提取出的当前帧图像。

图4为本发明血管自动跟踪过程示意图，图4（a）为当前帧颈动脉图像，图4（b）为其下一帧图像。

图5（a）为包含ROI1曲线在内的矩形窗图像，图5（b）为图5（a）经过stick filter得到的图像。

图6为图5（b）经极坐标转换得到的极坐标图。

图7（a）为MAB初始轮廓示意图，图7（b）为包含MAB分割区域和MAB初始轮廓的矩形窗示意图，图7（c）为水平集演化得到的最终的MAB轮廓与手工勾勒的金标准比较的结果示意图，图7（d）为包含LIB分割区域和LIB初始轮廓的矩形窗示意图，图7（e）为水平集演化得到的最终的LIB轮廓（“.-”连接构成）与手工勾勒的金标准（虚线表示）比较的结果示意图。

图8为最大期望算法估计感兴趣区域1（ROI1）内灰度概率密度分布示意图。

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

具体实施方式

下面结合附图和实施例对本发明做进一步的详细描述。

本发明提供的超声图像中颈动脉内外轮廓和斑块的计算机自动分割算法，其实施步骤如下：

（1）确定基准轮廓：

在计算机中载入一个病例的一组颈动脉三维超声体数据，如图2所示。计算机根据体素的大小自动调整体素与体积（mm³）的比例。本发明的分割对象就是颈动脉三维超声体数据的每一帧图像。若当前帧图像为颈动脉三维超声体数据的第一帧图像，则凭经验判断出外轮廓的大致位置，并勾取几个比较明显的位于血管外轮廓上的基准点，再通过插值的方法形成一条闭合曲线作为基准轮廓，如图3中虚线构成的轮廓所示；否则，将上一帧得到的备用基准轮廓作为当前帧的基准轮廓；

（2）利用形态学膨胀法将基准轮廓外推形成第一感兴趣区域ROI1；

将基准轮廓利用形态学膨胀方法向外推移一定距离，此距离一般可取10到20个像素，本实施例中定为10个像素的距离，形成感兴趣区域1（ROI1），见图3中的实心闭合轮廓线；

（3）在第一感兴趣区域ROI1中检测颈动脉血管外轮廓C_MAB：

（3.1）在感兴趣区域1（ROI1）中检测MAB的初始轮廓线

本步骤的技术思路是：在当前帧图像中提取包含第一感兴趣区域ROI1的矩形窗口图像，将其转化到极坐标图像，在极坐标图像上每隔相同角度在径向上靠近基准轮廓的像素点中搜索梯度值最大者作为初始轮廓点。所述间隔角度可任意选取，一般选5到20度；在搜索中可将搜索范围定为与基准轮廓间隔2~10个像素距离的像素点集。具体过程如下：

（3.1.1）计算ROI1区域的几何中心坐标，自动生成包含ROI1在内的矩形窗口图像，如图5（a）所示，图上的轮廓为经过坐标转化后显示在矩形窗口图像中的ROI1的边界；

（3.1.2）作为优化，对矩形窗口图像进行滤波，以消除部分噪声和增强颈动脉血管与周围组织的对比度，如图5（b）。此处可采用滤波的方法有：棒滤波（stick filter）、双边滤波（bilateral filter）、各向异性扩散滤波（SRAD filter）等，本实例采用棒滤波。棒滤波具体实现方法参见S.D.Pathak,V.Chalana,D.R.Haynor,and Y.Kim,“Edge-guided boundary delineation inprostate ultrasound images,”IEEE Trans.Med Imaging 19(12),1211–1219(2000)；

（3.1.3）将矩形窗口图像转化为极坐标图像，将极坐标图像沿径向均匀划分为N个区间，见图6（a），N的取值在18至72之间，在本实施例中区间个数N取24；取图像6（a）中每个区间径向上的中线，计算每条中线上梯度值最大的像素点，标记为“+”，如图6（b）所示。每个标记点的坐标表示为{M_i|r_i,θ_i＝(i-0.5)*2/N,i＝1,2,...N}，将这些点作为极坐标图上外轮廓MAB的初始轮廓点；

（3.1.4）将步骤（3.1.3）中的标记点从极坐标还原到直角坐标系中，平移后得到在原图中的坐标位置，并用细线依次连接起来，形成MAB的初始轮廓

如图7（ａ）所示；

（3.2）利用混合分布估计第一感兴趣区域ROI1的灰度概率密度分布，从而得到混合分布的shape参数K_j和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），C均值聚类（C-means）或模糊聚类（Fuzzy）等，在本实施例中采用了K均值聚类，且聚类类别设为3（聚类类别一般根据目标所包含的组织种类来设定）。为减少计算开销，本实例聚类迭代次数设为50，将样本点按照灰度均值分为3类，K均值聚类方法具体实现参见J.A.Hartigan,et al,A K-Means Clustering Algorithm,Journal ofthe Royal Statistical Society，197928(1),pp:100-108；

（3.2.3）分别在聚类得到的三类样本点上，采用最大似然估计法（Maximum Likelihood）对样本点进行gamma分布估计，得到三组gamma分布的参数(K₁,θ₁),(K₂,θ₂),(K₃,θ₃),K,θ分别对应gamma分布中的shape参数和scale参数。自动计算三类样本点在总样本点中所占比例W_i:

其中N_i表示样本中属于第i个分量的像素个数，i=1,2,3.N_s代表样本点的总个数。最大似然估计法的实现具体参见Aldrich,John,.A.Fisher,The making ofmaximum likelihood,Statistical Science,1997,12(3):162–176；

（3.2.4）将步骤（3.2.3）中的三组参数(W₁,K₁,θ₁),(W₂,K₂,θ₂),(W₃,K₃,θ₃)作为期望最大化（Expectation Maximization）算法的初始值代入到期望最大化算法中，迭代得到最终的三组参数(W₁,K₁,θ₁),(W₂,K₂,θ₂),(W₃,K₃,θ₃)作为混合gamma分布（mixture gamma distribution）的最终参数，这样ROI1内的灰度密度概率分布可近似表示为：

$P=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{W}_i\·G_i\left(I_x\right|K_i,{\mathrm{\θ}}_i)$

其中

$G\left(I_x\right|K_i,{\mathrm{\θ}}_i)=I_x^{K_i-1}\·\exp (-\frac{I_x}{{\mathrm{\θ}}_i})/({\mathrm{\θ}}_i^{K_i}\·(K_i-1)!)$

为gamma分布的函数通式，I_x为灰度级。

期望最大化算法的具体实现参见A.Dempster,N.Laird,and D.Rubin,“Maximum likelihood from incomplete data via the EM algorithm,”J.R.Stat.Soc.(Ser.B),vol.39,pp.1–38,1977。本实例的直方图分布及用期望最大化算法估计得到的gamma混合分布曲线如图8所示，其中三条实心曲线为混合gamma分布（mixture gamma distribution）的3个不同分量，虚线为3个不同分量组合得到的混合gamma分布曲线，表1为期望最大化算法得到的三个gamma分布分量的参数值。

表1  gamma混合分布的参数估计结果

i	Wi	Ki	θi
				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的初始轮廓C₁以外的像素区域，x₁为ROI1图像域中的像素点，‖‖为求欧式距离；

（3.4）计算

为减少计算开销，本方法自动生成一个包含MAB分割区域的矩形窗口图像，具体方法为：将基准轮廓采用形态学腐蚀法内推一段距离形成一条新的轮廓，此距离可根据图像大小恰当选择，本实施例中定位18个像素距离。这条形成的新的轮廓与感兴趣区域1（ROI1）边界（外边界）共同构成MAB的分割区域，水平集演化过程中演化轮廓仅在此分割区域内演化。MAB分割区域如图7（b）中非黑色区域所示，区域内的轮廓为MAB的初始轮廓；

MAB分割的能量泛函定义为：

$E_{{\mathrm{\φ}}_1}={\mathrm{\ρ}}_1\·E_{{\mathrm{pdf}}_1}+{\mathrm{\λ}}_1\·E_{s_1}$

其中，λ₁为

所占的权重，ρ₁为的权重。

为限制曲线长度的能量，保证演化曲线的光滑性，λ₁越大，曲线越光滑,其中

为正则化的Heaviside函数，ε为一个极小的正实数，通常取0＜ε＜0.001，Ω₁为ROI1图像域，

为梯度算子。

$E_{{\mathrm{pdf}}_1}=-{\∫}_{{\mathrm{\Ω}}_1}H\left({\mathrm{\φ}}_1\left(x_1\right)\right)\·\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ai}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)){\mathrm{dx}}_1-{\∫}_{{\mathrm{\Ω}}_1}(1-H\left({\mathrm{\φ}}_1\left(x_1\right)\right))\·\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ae}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)){\mathrm{dx}}_1$ 为水平集区域信息能量项，它利用最大后验概率准则（Maximization aPosteriori）的思想，按照灰度概率密度分布来划分像素点的类别。其中，M为混合gamma分布的分量个数，q_Ai(j),q_Ae(j)分别为演化轮廓内外区域每种gamma分布的权重，

为参数为K_j,θ_j的gamma分布函数,

是像素点x₁的灰度值，取值范围为0到255。

表示整幅图像灰度概率密度分布的似然性，

越小，目标和背景所包括的像素分离得越好。

为了得到使能量

达到最小的解φ₁(x₁)，采用变分法得到Euler-Lagrange方程:

$\frac{{{\∂\mathrm{\φ}}_1}^t\left(x_1\right)}{\∂t}=-\frac{{\∂E}_{{\mathrm{\φ}}_1}}{{{\∂\mathrm{\φ}}_1}^t\left(x_1\right)}$

计算

得到水平集φ₁(x₁)演化的偏微分方程:

$\frac{{\∂\mathrm{\φ}}_1\left(x_1\right)}{\∂t}=-{\mathrm{\ρ}}_1\·\mathrm{\δ}\left({\mathrm{\φ}}_1\left(x_1\right)\right)\·[\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ai}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)-\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ae}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)]$

$+{\mathrm{\λ}}_1\·\mathrm{\δ}\left({\mathrm{\φ}}_1\left(x_1\right)\right)\·\mathrm{div}\left(\frac{{\▿\mathrm{\φ}}_1\left(x_1\right)}{\left|{\▿\mathrm{\φ}}_1\left(x_1\right)\right|}\right)$

上式中，δ(x)为正则化的Dirac函数，定义为：

其中ε为一个很小的正实数，通常取0＜ε＜0.001。

计算φ₁ ^t(x₁)：

${\mathrm{\φ}}_1^t\left(x_1\right)={\mathrm{\φ}}_1^{t-1}\left(x_1\right)-{\mathrm{\ρ}}_1\·\mathrm{\δ}\left({\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right)\·\left[\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ai}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)-\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ae}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)\right]$

$+{\mathrm{\λ}}_1\·\mathrm{\δ}\left({\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right)\·\mathrm{div}\left(\frac{{\▿\mathrm{\φ}}_1^{t-1}\left(x_1\right)}{\left|{\▿\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right|}\right)$

（3.5）计算本次迭代演化外轮廓

$C_1^t=\left\{x_1\right|{\mathrm{\φ}}_1^t\left(x_1\right)=0\}$

（3.6）若本次迭代演化外轮廓

与前一次迭代得到的演化外轮廓

相同或者t达到迭代次数上限，则本次迭代得到的演化轮廓即为颈动脉血管外轮廓C_MAB，迭代结束；否则，t＝t+1，返回步骤（3.4）；

上式中，ρ₁取值范围在0到10之间，λ₁取值范围在0到20之间，在本实例中迭代次数上限为100次，ρ₁=0.2，λ₁=1.2，精确的MAB边界轮廓线C_MAB见图7（c）中“.”构成轮廓所示。

（4）在当前帧图像中将颈动脉血管外轮廓C_MAB围成的区域作为第二感兴趣区域ROI2；

（5）在第二感兴趣区域ROI2中检测颈动脉血管内轮廓C_LIB：

（5.1）利用形态学腐蚀法将颈动脉血管外轮廓C_MAB内推，将内推得到的轮廓线作为初始颈动脉血管内轮廓

将得到的MAB轮廓C_MAB围成的区域作为LIB轮廓分割的感兴趣区域2（ROI2），并利用形态学腐蚀法将MAB轮廓内推一段距离得到的轮廓线作为LIB的初始轮廓如图7（d）中的闭合轮廓所示，本实施例中内推距离取18个像素。

（5.2）初始化迭代次数t′＝1以及符号距离函数

表示第二感兴趣区域ROI2中初始颈动脉血管内轮廓

以外的像素区域，x₂为第二感兴趣区域ROI2中的像素点；

（5.3）计算

定义LIB分割的水平集能量为：

$E_{{\mathrm{\φ}}_2}={\mathrm{\ρ}}_2\·E_{{\mathrm{pqd}}_2}+{\mathrm{\λ}}_2\·E_{s_2}+{\mathrm{\β}}_2\·E_T$

其中

为结合灰度概率密度分布的区域信息能量项，用以估计图像每一点属于目标和背景的概率，为曲率项，用以控制曲线的光滑性，这两项能量项的定义与步骤（3.4）相同；ρ₂，λ₂，β₂分别为三项能量的权重；

E_T为MAB与LIB之间的距离保持项，定义为：

$E_T={\∫}_{{\mathrm{\Ω}}_2}\mathrm{\δ}\left({\mathrm{\φ}}_2\left(x_2\right)\right)\·B_T\left(x_2\right){\mathrm{dx}}_2$

其中Ω₂为ROI2图像域；

D(x₂,y)为LIB演化轮廓上的点x₂与MAB最终轮廓上的点y之间的欧氏距离，d_T为一般情况下MAB和LIB之间相隔的最小距离（经验值）。当点x₂与MAB轮廓的距离小于d_T时，就会受到一个排斥力B_T(x₂)；当点x₂与MAB轮廓的距离大于d_T时，排斥力为0。

为了得到使得能量

最小的解φ₂(x₂)，采用变分法得到Euler-Lagrange方程

$\frac{{{\∂\mathrm{\φ}}_2}^{t^{\′}}\left(x_2\right)}{\∂t}=-\frac{{\∂E}_{{\mathrm{\φ}}_2}}{{{\∂\mathrm{\φ}}_2}^{t^{\′}}\left(x_2\right)}$

计算

得到水平集φ₂(x₂)演化的偏微分方程:

$\frac{{\∂\mathrm{\φ}}_2\left(x_2\right)}{\∂t}=-{\mathrm{\ρ}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2\left(x_2\right)\right)\·\left[\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Li}}\left(j\right)\·G\left(I_{x_2}\right|K_j,{\mathrm{\θ}}_j)-\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Le}}\left(j\right)\·G\left(I_{x_2}\right|K_j,{\mathrm{\θ}}_j)\right]+$

${\mathrm{\λ}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2\left(x_2\right)\right)\·\mathrm{div}\left(\frac{{\▿\mathrm{\φ}}_2\left(x_2\right)}{\left|{\▿\mathrm{\φ}}_2\left(x_2\right)\right|}\right)+{\mathrm{\β}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2\left(x_2\right)\right)\·B_T\left(x_2\right)$

上式中，δ(x)为正则化的Dirac函数，定义为：

其中ε为一个很小的正实数，通常取0＜ε＜0.001。

更新φ₂ ^t′(x₂)：

${\mathrm{\φ}}_2^{t^{\′}}\left(x_2\right)={\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)-{\mathrm{\ρ}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·\left[\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Li}}\left(j\right)\·G\left(I_{x_2}\right|K_j,{\mathrm{\θ}}_j)-\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Le}}\left(j\right)\·G\left(I_{x_2}\right|K_j,{\mathrm{\θ}}_j)\right]$

$+{\mathrm{\λ}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·\mathrm{div}\left(\frac{{\▿\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)}{\left|{\▿\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right|}\right)+{\mathrm{\β}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·B_T\left(x_2\right)$

（5.4）计算本次迭代演化内轮廓

$C_2^{t^{\′}}=\left\{x_2\right|{\mathrm{\φ}}_2^{t^{\′}}\left(x_2\right)=0\}$

（5.5）若本次迭代演化内轮廓

与前一次迭代得到的演化内轮廓

相同或者达到迭代次数上限，则本次迭代得到的演化轮廓为颈动脉血管内轮廓C_LIB，迭代结束；否则，t′＝t′+1，返回步骤（5.3）；

ρ₂的取值范围在0到10之间，λ₂的取值范围在0到20之间，β₂取值范围较大，一般取5以上，d_T取14到20之间。

在本实例中迭代次数上限为8000次，ρ₂=0.15，λ₂=5，β₂=20，最小距离d_T=18个像素，得到最终的LIB边界轮廓线C_LIB，见图7（e）的由“·–”构成的轮廓，虚线构成的轮廓为用于对比的手工勾勒的LIB轮廓。

（6）跟踪确定下一帧图像的备用基准轮廓：

（6.1）确定跟踪区域Ω；将当前帧颈动脉MAB精确轮廓C_MAB通过形态学膨胀外推一定距离（一般取10到30个像素距离，本实施例中为20个像素）与之构成一块区域Ω，见图4（a）中两条轮廓所围成的区域，其中较小轮廓代表分割得到的MAB轮廓，较大轮廓代表MAB轮廓经过外推得到的轮廓。此区域可视为包含颈动脉外膜及其外周组织的区域，一般情况下每一帧颈动脉图像都有这样一个区域，本方法根据这个区域的偏转移动来实现血管的自动跟踪；

（6.2）计算当前帧与下一帧之间的血管跟踪映射关系，具体方法为：

（6.2.1）对区域Ω中任意一点(x₃,y₃)，假设其在下一帧图像上通过映射M:{(x₄,y₄)＝τ+D*(x₃,y₃)}移动到了点(x₄,y₄)，其中τ为平移向量，D为形变矩阵。本方法中认为连续帧的血管发生的形变很小，因此本实施例中，

$D=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$

假设τ＝(τ_x,τ_y)，其中τ_x为点(x₃,y₃)于映射过程中在x方向平移的距离，τ_y为点(x₃,y₃)于映射过程中在y方向平移的距离则对应点之间的映射关系可以表示为：

$\left(\begin{array}{c}{x}_4\\ y_4\end{array}\right)=\left(\begin{array}{c}{\mathrm{\τ}}_x\\ {\mathrm{\τ}}_y\end{array}\right)+\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}{x}_3\\ y_3\end{array}\right)$

并且本方法认为此区域中每个点均满足同样的映射，通过此映射当前帧中的颈动脉外周组织区域投影到了下一帧的新的区域；

（6.2.2）计算当前帧图像的区域Ω中点(x₃,y₃)和其在下一帧图像上的对应点(x₄,y₄)的灰度值I(x₃,y₃)和J(x₄,y₄)之间差异的平方：

ΔI²＝[I(x₃,y₃)-J(x₄,y₄)]²

对这个区域内每个点进行上述操作，将每个点上的灰度差异求和：

S＝∑_ΩΔI²

（6.2.3）分别在下一帧图片的横轴、纵轴方向上一定范围内，搜索到使S达到最小的区域，此区域对应的映射关系即为所求的连续两帧之间的血管跟踪映射关系M，本实施例中搜索范围在横轴纵轴上均取为正负20个像素；

(6.3)将当前帧图片的MAB精确轮廓C_MAB按照此映射关系M投影到下一帧上，构成下一帧颈动脉图像的备用基准轮廓，见图4（b）中的虚线轮廓。

本发明迭代次数上限一般为5000次以上。

分割得到的颈动脉内外轮廓可用于粥样硬化测量指标的计算，以辅助医生进行病理分析和药物治疗效果观察评价，下面对其详细说明：

测量指标的分析：

根据前面步骤精确分割得到的每一帧超声图像中颈动脉的MAB和LIB计算颈动脉粥样硬化的重要指标——血管壁体积（Vessel WallVolume,VWV）：

计算每一帧上分割得到的MAB和LIB轮廓之间的面积A_vessel(k)，（k=1,2，...,N，N为体数据的总帧数）;

计算每两帧图片之间的血管壁体积：

V_vessel(k)＝0.5*(A_vessel(k)+A_vessel(k+1))*d

将每两帧之间的血管壁体积求和，得到最终的VWV；

本实施例中的指标结果见表2，其中k为体数据中帧数的序号，k=1表示离颈动脉分岔处最远的起始帧。

表2测量指标的分析结果

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.超声颈动脉血管膜的自动分割方法，包括以下步骤：

（1）若当前帧图像为颈动脉三维超声体数据的第一帧图像，则在当前帧图像上选择明显位于血管外轮廓上的像素点集作为基准点，通过基准点间插值形成闭合曲线，将其作为当前帧的基准轮廓；否则，将上一帧得到的备用基准轮廓作为当前帧的基准轮廓；

（2）利用形态学膨胀法将基准轮廓外推形成第一感兴趣区域ROI1；

（3）在第一感兴趣区域ROI1中检测颈动脉血管外轮廓C_MAB：

（3.1）在当前帧图像中提取包含第一感兴趣区域ROI1的矩形窗口图像，将其转化到极坐标图像，在极坐标图像上每隔相同角度在径向方向靠近基准轮廓的像素点中搜索梯度值最大者作为初始轮廓点，将初始轮廓点转化到直角坐标系并连线形成初始颈动脉血管外轮廓

（3.2）利用混合分布估计第一感兴趣区域ROI1的灰度概率密度分布，从而得到混合分布的shape参数K_j和scale参数θ_j，j＝1,…,M，M为混合分布中单项分布的类别数；

（3.3）初始化迭代次数t=1以及符号距离函数

表示第一感兴趣区域ROI1中初始颈动脉血管外轮廓

以外的像素区域，x₁为第一感兴趣区域ROI1中的像素点，||  ||为求欧式距离；

（3.4）计算

${\mathrm{\φ}}_1^t\left(x_1\right)={\mathrm{\φ}}_1^{t-1}\left(x_1\right)-{\mathrm{\ρ}}_1\·\mathrm{\δ}\left({\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right)\·[\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ai}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)-\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Ae}}\left(j\right)\·G\left(I_{x_1}\right|K_j,{\mathrm{\θ}}_j)]$

$+{\mathrm{\λ}}_1\·\mathrm{\δ}\left({\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right)\·\mathrm{div}\left(\frac{{\▿\mathrm{\φ}}_1^{t-1}\left(x_1\right)}{\left|{\▿\mathrm{\φ}}_1^{t-1}\left(x_1\right)\right|}\right)$

其中，

这里0＜ε＜0.001，0＜ρ₁≤10，0＜λ₁≤20，q_Ai(j)和q_Ae(j)分别为前一次迭代演化轮廓

的内、外区域中第j项单项分布的权重，

为参数为K_j,θ_j的第_j项单项分布函数，

为像素点x₁的灰度值；

（3.5）计算本次迭代演化外轮廓

$C_1^t=\left\{x_1\right|{\mathrm{\φ}}_1^t\left(x_1\right)=0\}$

（3.6）若本次迭代演化外轮廓

与前一次迭代得到的演化外轮廓

相同或者t达到迭代次数上限，则本次迭代得到的演化轮廓即为颈动脉血管外轮廓C_MAB，迭代结束，进入步骤（4）；否则，t＝t+1，返回步骤（3.4）；

（4）在当前帧图像中将颈动脉血管外轮廓C_MAB围成的区域作为第二感兴趣区域ROI2；

（5）在第二感兴趣区域ROI2中检测颈动脉血管内轮廓C_LIB：

（5.1）利用形态学腐蚀法将颈动脉血管外轮廓C_MAB内推，将内推得到的轮廓线作为初始颈动脉血管内轮廓

（5.2）初始化迭代次数t′＝1以及符号距离函数

表示第二感兴趣区域ROI2中初始颈动脉血管内轮廓

以外的像素区域，x₂为第二感兴趣区域ROI2中的像素点；

（5.3）计算

${\mathrm{\φ}}_2^{t^{\′}}\left(x_2\right)={\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)-{\mathrm{\ρ}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·\left[\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Li}}\left(j\right)\·G\left(I_{x_3}\right|K_j,{\mathrm{\θ}}_j)-\ln (\underset{j}{\overset{M}{\mathrm{\Σ}}}{q}_{\mathrm{Le}}\left(j\right)\·G\left(I_{x_3}\right|K_j,{\mathrm{\θ}}_j)\right]$

$+{\mathrm{\λ}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·\mathrm{div}\left(\frac{{\▿\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)}{\left|{\▿\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right|}\right)+{\mathrm{\β}}_2\·\mathrm{\δ}\left({\mathrm{\φ}}_2^{t^{\′}-1}\left(x_2\right)\right)\·B_T\left(x_2\right)$

其中，0＜ρ₂≤10，0＜λ₂≤20，0＜β₂≤100，q_Li(j)和q_Le(j)分别为前一次演化轮廓

的内、外区域中第j项单项分布的权重，为参数为K_j,θ_j的第j项单项分布函数，

为像素点x₂的灰度值；

（5.4）计算本次迭代演化内轮廓

$C_2^{t^{\′}}=\left\{x_2\right|{\mathrm{\φ}}_2^{t^{\′}}\left(x_2\right)=0\}$

（5.5）若本次迭代演化内轮廓

与前一次迭代得到的演化内轮廓

相同或者达到迭代次数上限，则本次迭代得到的演化轮廓为颈动脉血管内轮廓C_LIB，迭代结束，进入步骤（6）；否则，t′＝t′+1，返回步骤（5.3）；

（6）跟踪确定下一帧图像的备用基准轮廓：

（6.1）在当前帧图像中，通过形态学膨胀法将颈动脉血管外轮廓C_MAB外推，将外推得到的闭合轮廓线与C_MAB构成的区域作为跟踪区域Ω；

（6.2）在下一帧图像中搜索使得灰度差异和S＝∑_ΩΔI²最小的像素点集(x₄,y₄)，其中ΔI²＝[I(x₃,y₃)-J(x₄,y₄)]²，I(x₃,y₃)为当前帧图像的跟踪区域Ω中像素点(x₃,y₃)的灰度值，J(x₄,y₄)为下一帧图像的像素点(x₄,y₄)的灰度值；

（6.3）令当前帧图像的跟踪区域中的点(x₃,y₃)与下一帧图像的像素点(x₄,y₄)的映射关系为 $\left(\begin{array}{c}{x}_4\\ y_4\end{array}\right)=\left(\begin{array}{c}{\mathrm{\τ}}_x\\ {\mathrm{\τ}}_y\end{array}\right)+\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}{x}_3\\ y_3\end{array}\right),$ 结合像素点(x₃,y₃)与(x₄,y₄)求解映射关系中的参数τ＝(τ_x,τ_y)，从而确定映射关系；

（6.4）根据步骤（6.3）建立的映射关系将当前帧图像的颈动脉血管外轮廓C_MAB映射为下一帧图像的备用基准轮廓。

2.根据权利要求1所述的超声颈动脉血管膜的自动分割方法，其特征在于，所述混合分布为混合gamma分布或混合gauss分布或混合rayleigh分布。

Images