CN103996194A - 一种基于超声颈动脉图像的内中膜自动分割方法 - Google Patents

Abstract

本发明公开了一种基于超声颈动脉图像的内中膜自动分割方法，包括：获取超声颈动脉纵截面图像，对该超声颈动脉纵截面图像进行灰度归一化处理，以得到归一化超声颈动脉纵截面图像，对该归一化超声颈动脉纵截面图像进行处理，提取去噪后的归一化超声颈动脉纵截面图像的特征图，计算特征图中每个像素点与其纵向相邻像素点的特征值梯度，将特征值梯度分布最显著的横轴作为感兴趣区域的中间横轴，从中间横轴分别向上和向下扩展，以得到感兴趣区域，获取感兴趣区域R的灰度梯度图，在估计的内中膜厚度的取值范围内依次取每一个整数值作为可能的内中膜厚度，根据该可能的内中膜厚度设定特征提取算子。本发明具有实施方便、计算迅速、鲁棒性强的优点。

Description

一种基于超声颈动脉图像的内中膜自动分割方法

技术领域

本发明属于计算机技术与医学图像处理交叉领域，更具体地，涉及一种基于超声颈动脉图像的内中膜自动分割方法。

背景技术

缺血性心脑血管疾病严重危害人类健康，颈动脉粥样硬化是此类疾病的主要致病因素，通过对颈动脉粥样硬化的预防、早期检测及治疗，是降低患缺血性心脑血管疾病发病风险的有效手段。已有大量实验证明，颈动脉内中膜增厚是发生颈动脉粥样硬化的前期临床表现，可通过对颈动脉内中膜厚度及其变化情况的预测来判断是否存在颈动脉粥样硬化风险。其中，内中膜厚度(Intima-Media Thickness，简称IMT)被定义为颈动脉血管壁远端内腔-内膜界面(Lumen Intima Interface,简称LII)与中膜-外膜界面(Media Adventitia Interface,简称MAI)之间的距离。医学超声具有无创、无辐射、廉价、实时、方便等特点，可获得清晰的颈动脉图像，故已广泛应用于颈动脉内中膜分析中。

传统测量颈动脉内中膜的方法，是由有经验的医生在从超声机上获取的颈动脉纵截面图像的LII和MAI上标记数组相对应的点对，并用求平均值的方式计算内中膜的平均厚度。这种方法存在明显的缺陷：取点过程完全由临床医生操作，点对位置的选择及标记均会出现组间差异和组内差异，即不同的医生标记的结果会存在差异，同一医生在不同时间标记的结果也不完全相同。计算过程中，因只选择少数点对，故必然会遗漏点对与点对之间的内中膜的信息。

针对传统方法的缺点，近年来研究者们提出了许多测量颈动脉内中膜厚度的方法，包括有监督方法和无监督方法等。因有监督方法存在需要人为初始化，且结果受初始化影响的弊端，故无监督方法是该领域研究的重点。无监督的方法均根据颈动脉超声图像的特点进行设计，通过不同途径将内中膜的边缘位置与颈动脉图像各区域间的差异性以及区域内的一致性相结合，主要包括基于梯度的方法、基于像素聚类的方法、基于函数值最小化的方法、基于活动轮廓模型的方法、基于数学形态学的方法等等。但由于超声图像往往分辨率得不到保证，图像存在较严重的斑点噪声和伪影，导致许多图像分割方法并不能直接应用，而以上各种方法即使针对超声图像进行优化，也会存在各种不同的问题：相对简单的基于图像梯度、灰度等特点的方法鲁棒性不强，实际中往往不能使用，必须与其他方法相结合；基于活动轮廓模型的方法是近年来研究的热点，此类方法的实现较其他方法复杂，而如何减少初始轮廓对最终结果的影响则是研究者们所要考虑的另一个问题；传统的基于动态规划法、A*算法的方法具有较强的理论基础，但受参数影响大问题限制了其扩展到不同类型图像上的应用。因此，如何减少此类方法在实施过程中对参数的依赖是研究的核心问题，且此类方法亦会受到初始化的影响。

综上所述，由于在临床应用中对方法的计算速度和对结果的精度均有要求，故对内中膜的自动分割仍存在较大困难。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于超声颈动脉图像的内中膜自动分割方法，其仅需将待分割的颈动脉纵截面超声图像作为输入，即可输出已分割出内中膜部分的结果图像，以供进一步的参数采集及临床医生诊断治疗，且本发明具有实施方便、计算迅速、鲁棒性强的优点。

为实现上述目的，按照本发明的一个方面，提供了一种基于超声颈动脉图像的内中膜自动分割方法，包括以下步骤：

(1)获取超声颈动脉纵截面图像，对该超声颈动脉纵截面图像进行灰度归一化处理，以得到归一化超声颈动脉纵截面图像，并利用图像去噪算法对该归一化超声颈动脉纵截面图像进行处理；

(2)提取去噪后的归一化超声颈动脉纵截面图像的特征图，计算特征图中每个像素点与其纵向相邻像素点的特征值梯度，将特征值梯度分布最显著的横轴作为感兴趣区域的中间横轴，从中间横轴分别向上和向下扩展，以得到感兴趣区域；

(3)获取感兴趣区域R的灰度梯度图，在估计的内中膜厚度的取值范围内依次取每一个整数值作为可能的内中膜厚度，根据该可能的内中膜厚度设定特征提取算子，并用该算子计算灰度梯度图中每个像素点的特征值，以得到与该可能的内中膜厚度对应的灰度梯度图的特征图，对获得的每个特征图的所有特征值进行求和处理，和最大的特征图所对应的可能的内中膜厚度作为内中膜厚度参考值；

(4)根据内中膜厚度参考值构建能量函数，并利用能量函数采用基于函数值最小化算法获得内中膜的内腔-内膜界面；

(5)根据内中膜厚度参考值构建能量函数，并利用能量函数采用基于函数值最小化算法获得内中膜的中膜-外膜界面。

优选地，步骤(1)中采用的是SRAD算法。

优选地，步骤(2)具体包括以下子步骤：

(2-1)计算去噪后的归一化超声颈动脉纵截面图像中每一个像素点的即向上光滑性游程，以形成去噪后的归一化超声颈动脉纵截面图像的向上光滑性游程图，其中向上光滑性游程的计算公式如下：

$D(x,y)={(I(x,y)-I(x,y+1))}^{n_1}+n_2\×I(x,y),(x,y)\∈R_0$

$\left\{\begin{array}{c}M(x,y,L)=\underset{n=0}{\overset{L}{\mathrm{\Σ}}}D(x,y+n)\\ M(x,y,L-1)\≤K\\ M(x,y,L)>K\\ P(x,y)=L\end{array}\right.,(x,y)\∈R_0$

其中D(x,y)为去噪后的归一化超声颈动脉纵截面图像中每一个像素点的向上光滑性度量，x和y分别为像素点的横坐标和纵坐标，R₀为去噪后的归一化超声颈动脉纵截面图像，I(x,y)为去噪后的归一化超声颈动脉纵截面图像中每个像素点的灰度值，n₁为梯度惩罚系数，n₂为灰度惩罚系数，P(x,y)为去噪后的归一化超声颈动脉纵截面图像中每一个像素点的向上光滑性游程，K为光滑性阈值，L和M均为中间变量；

(2-2)计算步骤(2-1)获取的向上光滑性游程图中每个向上光滑性游程的梯度，以获得向上光滑性游程图的梯度图G₁，具体采用以下公式：

$G_1(x,y)=\left\{\begin{array}{c}P(x,y)-P(x,y-1),P(x,y)\&GreaterEqual;P(x,y-1)\\ 0,P(x,y)<P(x,y-1)\end{array}\right.,(x,y)\&Element;R_0$

(2-3)对向上光滑性游程图的梯度图G₁进行横向求和，并对横向求和结果进行邻域相加处理，以得到数列S₂，具体为：

$S_1\left(y\right)=\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{G}_1(i,y),y\&Element;[1,Y]$

$S_2\left(y\right)=\left[\begin{array}{ccccccc}1& 2& 4& 6& 4& 2& 1\end{array}\right]\left[\begin{array}{c}{S}_1(y-3)\\ S_1(y-2)\\ S_1(y-1)\\ S_1\left(y\right)\\ S_1(y+1)\\ S_1(y+2)\\ S_1(y+3)\end{array}\right],y\&Element;[1,Y]$

其中S₁为所横向求和得到的数列，S₂为经过邻域求和处理得到的数列，X和Y分别为去噪后的归一化超声颈动脉纵截面图像的宽和高；

(2-4)将数列S₂中最大项的序号所对应的去噪后的归一化超声颈动脉纵截面图像中的行号作为感兴趣区域的中间横轴，从该中间横轴分别向上和向下扩展，以得到感兴趣区域R。

优选地，向上扩展和向下扩展的高度可相同，也可不同，其取值是由去噪后的归一化超声颈动脉纵截面图像的尺寸、向上光滑性游程图中向上光滑性游程的最大值决定。

优选地，步骤(3)具体包括以下子步骤：

(3-1)获取感兴趣区域R的梯度图G₂，具体为：

$G_2^{\&prime;}(x,y)=\left[\begin{array}{ccc}-1& 0.8& 0.2\end{array}\right]\left[\begin{array}{c}I(x,y)\\ I(x,y-1)\\ I(x,y-2)\end{array}\right],(x,y)\&Element;R$

$G_2(x,y)=\left\{\begin{array}{c}{G}_2^{\&prime;}(x,y),G_2^{\&prime;}(x,y)\&GreaterEqual;0\\ 0,G_2^{\&prime;}(x,y)<0\end{array}\right.,(x,y)\&Element;R$

其中G₂′为中间变量；

(3-2)根据感兴趣区域R的大小估计内中膜厚度的取值范围[D₁,D₂]；具体而言，该范围被设置为一个像素点至感兴趣区域R宽度的一半；

(3-3)依次取每一个整数值作为可能的内中膜厚度，根据该可能的内中膜厚度设定特征提取算子，并用该算子计算灰度梯度图中每个像素点的特征值，以得到与该可能的内中膜厚度对应的灰度梯度图的特征图F_d，具体为：

F_d(x,y)＝F_d1(x,y)×F_d2(x,y),(x,y)∈R

其中，

$F_{d1}(x,y)=\left\{\begin{array}{c}{G}_2(x,y+d)-G_2(x,y),G_2(x,y+d)\&GreaterEqual;G_2(x,y)\\ 0,G_2(x,y+d)<G_2(x,y)\end{array}\right.,(x,y)\&Element;R$

$F_{d2}(x,y)=\left\{\begin{array}{c}{G}_2(x,y-d)-G_2(x,y),G_2(x,y-d)\&GreaterEqual;G_2(x,y)\\ 0,G_2(x,y-d)<G_2(x,y)\end{array}\right.,(x,y)\&Element;R$

d为估计的内中膜厚度取值范围[D₁,D₂]中的一个整数值，F_d1和F_d2为中间变量；

(3-4)对获得的每个特征图的所有特征值进行求和处理，和S₃最大的特征图所对应的可能的内中膜厚度d*为内中膜厚度参考值，其中S₃具体为：

$S_3\left(d\right)=\underset{(x,y)\&Element;R}{\mathrm{\&Sigma;}}{F}_d(x,y),d\&Element;[D_1,D_2].$

优选地，步骤(4)具体包括以下子步骤：

(4-1)对感兴趣区域进行对比度增强操作，以得到优化图像I₁，具体为：

$I_1(x,y)=1-\cos (I(x,y)\&times;\frac{\mathrm{\&pi;}}{2}),(x,y)\&Element;R$

(4-2)根据内中膜厚度参考值d^*获取感兴趣区域中像素点下方梯度的特征图G_down，具体为：

$G_{\mathrm{down}}(x,y)=\underset{i=d^{*}}{\overset{{3d}^{*}}{\mathrm{\&Sigma;}}}{G}_2(x,y-i),(x,y)\&Element;R$

(4-3)根据上述步骤(4-1)和(4-2)获得的结果构建能量函数E₁，具体为：

E₁＝w₁₁C₁₁+w₁₂C₁₂+w₁₃C₁₃+w₁₄C₁₄

其中，w₁₁、w₁₂、w₁₃、w₁₄为权重，C₁₁为轮廓线性度量函数，C₁₂为灰度项，C₁₃为灰度梯度项，C₁₄为轮廓下方区域梯度项；

(4-4)使用函数最小化方法使能量函数E₁最小化，即得到内中膜的内腔-内膜界面LII。

优选地，步骤(5)具体包括以下子步骤：

(5-1)获取感兴趣区域中像素点与内中膜的内腔-内膜界面LII位置关系的特征图G_up，具体为：

$G_{\mathrm{up}}(x,y)={1000}^{\frac{|y+{2\&times;d}^{*}-\mathrm{LII}|\left(x\right)}{d^{*}}-1},(x,y)\&Element;R$

(5-2)根据上述步骤(5-1)获得的结果构建能量函数E₂，具体为：

E₂＝w₂₁C₂₁+w₂₂C₂₂+w₂₃C₂₃+w₂₄C₂₄

其中，w₂₁、w₂₂、w₂₃、w₂₄为权重，C₂₁为轮廓线性度项，C₂₂为灰度项，C₂₃为灰度梯度项，C₂₄为轮廓与LLI的位置相关项；

(5-3)使用函数最小化方法使能量函数E₂最小化，即得到内中膜的中膜-外膜界面MAI。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

1、本发明在使用中由于不需要使用者在图像中作任何标记，故提供了一种自动分割超声颈动脉纵截面图像内中膜方法；

2、本发明提供的方法步骤简单，计算速度快，鲁棒性好；

3、本发明实际使用时内中膜测量的结果受图像类型变化影响较小。

附图说明

图1是本发明基于超声颈动脉图像的内中膜自动分割方法的流程图。

图2示出超声颈动脉纵截面原始图像。

图3是图像的向上光滑性游程图。

图4是图像向上光滑性游程图的梯度图。

图5示出感兴趣区域。

图6示出超声颈动脉纵截面图像内中膜分割结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示，本发明基于超声颈动脉图像的内中膜自动分割方法包括以下步骤：

(1)获取超声颈动脉纵截面图像，对该超声颈动脉纵截面图像进行灰度归一化处理，以得到归一化超声颈动脉纵截面图像，并利用图像去噪算法对该归一化超声颈动脉纵截面图像进行处理；具体而言，本发明中采用的是斑点噪声各向异性扩散(Speckle Reducing Anisotropic Diffusion,简称SRAD)算法；

(2)提取去噪后的归一化超声颈动脉纵截面图像的特征图，计算特征图中每个像素点与其纵向相邻像素点的特征值梯度，将特征值梯度分布最显著的横轴作为感兴趣区域的中间横轴，从中间横轴分别向上和向下扩展，以得到感兴趣区域；本步骤具体包括以下子步骤：

(2-1)计算去噪后的归一化超声颈动脉纵截面图像中每一个像素点的特征值(即向上光滑性游程)，以形成去噪后的归一化超声颈动脉纵截面图像的特征图(即向上光滑性游程图)，其中向上光滑性游程的计算公式如下：

$D(x,y)={(I(x,y)-I(x,y+1))}^{n_1}+n_2\&times;I(x,y),(x,y)\&Element;R_0$

$\left\{\begin{array}{c}M(x,y,L)=\underset{n=0}{\overset{L}{\mathrm{\&Sigma;}}}D(x,y+n)\\ M(x,y,L-1)\&le;K\\ M(x,y,L)>K\\ P(x,y)=L\end{array}\right.,(x,y)\&Element;R_0;$

其中D(x,y)为去噪后的归一化超声颈动脉纵截面图像中每一个像素点的向上光滑性度量，x和y分别为像素点的横坐标和纵坐标，R₀为去噪后的归一化超声颈动脉纵截面图像，I(x,y)为去噪后的归一化超声颈动脉纵截面图像中每个像素点的灰度值，n₁＝4为梯度惩罚系数，n₂＝0.5为灰度惩罚系数，P(x,y)为去噪后的归一化超声颈动脉纵截面图像中每一个像素点的向上光滑性游程，K为光滑性阈值，其取值与去噪后的归一化超声颈动脉纵截面图像的质量有关，质量越差，K值越大，反之则越小，L和M均为中间变量；

(2-2)计算步骤(2-1)获取的向上光滑性游程图中每个向上光滑性游程的梯度，以获得向上光滑性游程图的梯度图G₁，具体采用以下公式：

$G_1(x,y)=\left\{\begin{array}{c}P(x,y)-P(x,y-1),P(x,y)\&GreaterEqual;P(x,y-1)\\ 0,P(x,y)<P(x,y-1)\end{array}\right.,(x,y)\&Element;R_0$

(2-3)对向上光滑性游程图的梯度图G₁进行横向求和，并对横向求和结果进行邻域相加处理，以得到数列S₂，具体为：

$S_1\left(y\right)=\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{G}_1(i,y),y\&Element;[1,Y]$

$S_2\left(y\right)=\left[\begin{array}{ccccccc}1& 2& 4& 6& 4& 2& 1\end{array}\right]\left[\begin{array}{c}{S}_1(y-3)\\ S_1(y-2)\\ S_1(y-1)\\ S_1\left(y\right)\\ S_1(y+1)\\ S_1(y+2)\\ S_1(y+3)\end{array}\right],y\&Element;[1,Y]$

其中S₁为所横向求和得到的数列，S₂为经过邻域求和处理得到的数列，X和Y分别为去噪后的归一化超声颈动脉纵截面图像的宽和高；

(2-4)将数列S₂中最大项的序号所对应的去噪后的归一化超声颈动脉纵截面图像中的行号作为感兴趣区域的中间横轴，从该中间横轴分别向上和向下扩展，以得到感兴趣区域R；具体而言，向上扩展和向下扩展的高度可相同，也可不同，其取值是由去噪后的归一化超声颈动脉纵截面图像的尺寸、向上光滑性游程图中向上光滑性游程的最大值决定；

(3)获取感兴趣区域R的灰度梯度图，在估计的内中膜厚度的取值范围内依次取每一个整数值作为可能的内中膜厚度，根据该可能的内中膜厚度设定特征提取算子，并用该算子计算灰度梯度图中每个像素点的特征值，以得到与该可能的内中膜厚度对应的灰度梯度图的特征图，对获得的每个特征图的所有特征值进行求和处理，和最大的特征图所对应的可能的内中膜厚度作为内中膜厚度参考值；本步骤具体包括以下子步骤：

(3-1)获取感兴趣区域R的梯度图G₂，具体为：

$G_2^{\&prime;}(x,y)=\left[\begin{array}{ccc}-1& 0.8& 0.2\end{array}\right]\left[\begin{array}{c}I(x,y)\\ I(x,y-1)\\ I(x,y-2)\end{array}\right],(x,y)\&Element;R$

$G_2(x,y)=\left\{\begin{array}{c}{G}_2^{\&prime;}(x,y),G_2^{\&prime;}(x,y)\&GreaterEqual;0\\ 0,G_2^{\&prime;}(x,y)<0\end{array}\right.,(x,y)\&Element;R$

其中G₂′为中间变量；

(3-2)根据感兴趣区域R的大小估计内中膜厚度的取值范围[D₁,D₂]；具体而言，该范围被设置为一个像素点至感兴趣区域R宽度的一半；

(3-3)依次取每一个整数值作为可能的内中膜厚度，根据该可能的内中膜厚度设定特征提取算子，并用该算子计算灰度梯度图中每个像素点的特征值，以得到与该可能的内中膜厚度对应的灰度梯度图的特征图F_d，具体为：

F_d(x,y)＝F_d1(x,y)×F_d2(x,y),(x,y)∈R

其中，

$F_{d1}(x,y)=\left\{\begin{array}{c}{G}_2(x,y+d)-G_2(x,y),G_2(x,y+d)\&GreaterEqual;G_2(x,y)\\ 0,G_2(x,y+d)<G_2(x,y)\end{array}\right.,(x,y)\&Element;R$

$F_{d2}(x,y)=\left\{\begin{array}{c}{G}_2(x,y-d)-G_2(x,y),G_2(x,y-d)\&GreaterEqual;G_2(x,y)\\ 0,G_2(x,y-d)<G_2(x,y)\end{array}\right.,(x,y)\&Element;R$

d为估计的内中膜厚度取值范围[D₁,D₂]中的一个整数值，F_d1和F_d2为中间变量；

(3-4)对获得的每个特征图的所有特征值进行求和处理，和S₃最大的特征图所对应的可能的内中膜厚度d*为内中膜厚度参考值，其中S₃具体为：

$S_3\left(d\right)=\underset{(x,y)\&Element;R}{\mathrm{\&Sigma;}}{F}_d(x,y),d\&Element;[D_1,D_2].$

(4)根据内中膜厚度参考值构建能量函数，并利用能量函数采用基于函数值最小化算法获得内中膜的内腔-内膜界面；本步骤具体包括以下子步骤：

(4-1)对感兴趣区域进行对比度增强操作，以得到优化图像I₁，具体为：

$I_1(x,y)=1-\cos (I(x,y)\&times;\frac{\mathrm{\&pi;}}{2}),(x,y)\&Element;R$

(4-2)根据内中膜厚度参考值d*获取感兴趣区域中像素点下方梯度的特征图G_down，具体为：

$G_{\mathrm{down}}(x,y)=\underset{i=d^{*}}{\overset{{3d}^{*}}{\mathrm{\&Sigma;}}}{G}_2(x,y-i),(x,y)\&Element;R$

(4-3)根据上述步骤(4-1)和(4-2)获得的结果构建能量函数E₁，具体为：

E₁＝w₁₁C₁₁+w₁₂C₁₂+w₁₃C₁₃+w₁₄C₁₄

其中，w₁₁、w₁₂、w₁₃、w₁₄为权重，w₁₁和w₁₂取正值，w₁₃和w₁₄取负值，C₁₁为轮廓线性度量函数，C₁₂为灰度项，C₁₃为灰度梯度项，C₁₄为轮廓下方区域梯度项；C₁₁与y₁,...,y_X构成的轮廓的线性程度有关，当轮廓线性程度越高时，其值越小；C₁₂与y₁,...,y_X构成的轮廓所在位置的像素点灰度有关，当像素点灰度值越小时，其值越小；C₁₃与y₁,...,y_X构成的轮廓所处位置的像素点灰度梯度有关，当像素点灰度梯度值越大时，其值越小，优选地，此处的灰度梯度可定义为像素点与其下方一定邻域内像素点灰度的差值；C₁₄与y₁,...,y_X构成的轮廓所处位置下方一定区域内的像素点灰度梯度有关，当此区域内像素点灰度梯度值越大时，其值越小，其中，邻域的定义应主要以步骤(3)中求得的内中膜厚度参考值作为参考，表示轮廓下方的该区域内存在中膜-外膜界面的可能性较大，优选地，C₁₄的取值与该区域内各像素点灰度梯度总和成正相关。

具体而言，在本实施方式中，

$E_1(y_1,...,y_X)=w_{11}\&times;C(y_1,...,y_X)+w_{12}\&times;\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{I}_1(i,y_i)+w_{13}\&times;\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{G}_2(i,y_i)+w_{14}\&times;\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{G}_{\mathrm{down}}(i,y_i)$

其中C为y₁,...,y_X构成的轮廓的线性度量函数；

(4-4)使用函数最小化方法使能量函数E₁最小化，即得到内中膜的内腔-内膜界面LII；具体而言，本发明使用的是动态规划法。

(5)根据内中膜厚度参考值构建能量函数，并利用能量函数采用基于函数值最小化算法获得内中膜的中膜-外膜界面；本步骤具体包括以下子步骤：

(5-1)获取感兴趣区域中像素点与内中膜的内腔-内膜界面LII位置关系的特征图G_up，具体为：

$G_{\mathrm{up}}(x,y)={1000}^{\frac{|y+{2\&times;d}^{*}-\mathrm{LII}|\left(x\right)}{d^{*}}-1},(x,y)\&Element;R$

(5-2)根据上述步骤(5-1)获得的结果构建能量函数E₂，具体为：

E₂＝w₂₁C₂₁+w₂₂C₂₂+w₂₃C₂₃+w₂₄C₂₄

其中，w₂₁、w₂₂、w₂₃、w₂₄为权重，w₂₁和w₂₄取正值，w₂₂和w₂₃取负值,C₂₁为轮廓线性度项，C₂₂为灰度项，C₂₃为灰度梯度项，C₂₄为轮廓与LLI的位置相关项；C₂₁与y₁,...,y_X构成的轮廓的线性程度有关，当轮廓线性程度越高时，其值越小；C₂₂与y₁,...,y_X构成的轮廓所处位置的像素点灰度有关，当像素点灰度值越大时，其值越小；C₂₃与y₁,...,y_X构成的轮廓所处位置的像素点灰度梯度有关，当像素点灰度梯度值越大时，其值越小，优选地，此处的灰度梯度可定义为像素点与其下方一定邻域内像素点灰度的差值；C₂₄与y₁,...,y_X构成的轮廓与LLI的相对位置有关，当轮廓关于LLI的位置合理时，其值较小，否则，其值较大。优选地，将内中膜厚度参考值作为参考，将LLI向下平移一定距离作为参考轮廓，将参考轮廓周围的一定邻域设为合理位置。

具体而言，在本实施方式中，

$E_2(y_1,...,y_X)=w_{21}\&times;C(y_1,...,y_X)+w_{22}\&times;\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{I}_1(i,y_i)+w_{23}\&times;\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{G}_2(i,y_i)+w_{24}\&times;\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{G}_{\mathrm{up}}(i,y_i)$

(5-3)使用函数最小化方法使能量函数E₂最小化，即得到内中膜的中膜-外膜界面MAI；具体而言，本发明使用的是动态规划法。

实例：

利用本发明所述方法进行了大量实现，现将其中四个实例的实验结果展示如下，其中，实验中涉及到的变量均如表一所示：

表一

此四个实例的原始图像选自四幅不同类型的颈动脉纵截面超声图像(图2)。

基于上述四幅图像，分别使用以上所描述的步骤进行处理。其中，步骤(2)中可得到图像的向上光滑性游程图(图3)、向上光滑性游程图的梯度图(图4)、感兴趣区域(图5)。

最终得到四幅图像的分割结果如图6所示。

四个实例在图2～6中分别编号为(a)、(b)、(c)、(d)。

由实验结果可以看出，本发明涉及的分割方法分割精度高，且具有良好鲁棒性。

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

Claims (7)

1.一种基于超声颈动脉图像的内中膜自动分割方法，其特征在于，包括以下步骤：

(1)获取超声颈动脉纵截面图像，对该超声颈动脉纵截面图像进行灰度归一化处理，以得到归一化超声颈动脉纵截面图像，并利用图像去噪算法对该归一化超声颈动脉纵截面图像进行处理；

(2)提取去噪后的归一化超声颈动脉纵截面图像的特征图，计算特征图中每个像素点与其纵向相邻像素点的特征值梯度，将特征值梯度分布最显著的横轴作为感兴趣区域的中间横轴，从中间横轴分别向上和向下扩展，以得到感兴趣区域；

(3)获取感兴趣区域R的灰度梯度图，在估计的内中膜厚度的取值范围内依次取每一个整数值作为可能的内中膜厚度，根据该可能的内中膜厚度设定特征提取算子，并用该算子计算灰度梯度图中每个像素点的特征值，以得到与该可能的内中膜厚度对应的灰度梯度图的特征图，对获得的每个特征图的所有特征值进行求和处理，和最大的特征图所对应的可能的内中膜厚度作为内中膜厚度参考值；

(4)根据内中膜厚度参考值构建能量函数，并利用能量函数采用基于函数值最小化算法获得内中膜的内腔-内膜界面；

(5)根据内中膜厚度参考值构建能量函数，并利用能量函数采用基于函数值最小化算法获得内中膜的中膜-外膜界面。

2.根据权利要求1所述的内中膜自动分割方法，其特征在于，步骤(1)中采用的是SRAD算法。

3.根据权利要求1所述的内中膜自动分割方法，其特征在于，步骤(2)具体包括以下子步骤：

(2-1)计算去噪后的归一化超声颈动脉纵截面图像中每一个像素点的即向上光滑性游程，以形成去噪后的归一化超声颈动脉纵截面图像的向上光滑性游程图，其中向上光滑性游程的计算公式如下：

$D(x,y)={(I(x,y)-I(x,y+1))}^{n_1}+n_2\&times;I(x,y),(x,y)\&Element;R_0$

$\left\{\begin{array}{c}M(x,y,L)=\underset{n=0}{\overset{L}{\mathrm{\&Sigma;}}}D(x,y+n)\\ M(x,y,L-1)\&le;K\\ M(x,y,L)>K\\ P(x,y)=L\end{array}\right.,(x,y)\&Element;R_0$

其中D(x,y)为去噪后的归一化超声颈动脉纵截面图像中每一个像素点的向上光滑性度量，x和y分别为像素点的横坐标和纵坐标，R₀为去噪后的归一化超声颈动脉纵截面图像，I(x,y)为去噪后的归一化超声颈动脉纵截面图像中每个像素点的灰度值，n₁为梯度惩罚系数，n₂为灰度惩罚系数，P(x,y)为去噪后的归一化超声颈动脉纵截面图像中每一个像素点的向上光滑性游程，K为光滑性阈值，L和M均为中间变量；

(2-2)计算步骤(2-1)获取的向上光滑性游程图中每个向上光滑性游程的梯度，以获得向上光滑性游程图的梯度图G₁，具体采用以下公式：

$G_1(x,y)=\left\{\begin{array}{c}P(x,y)-P(x,y-1),P(x,y)\&GreaterEqual;P(x,y-1)\\ 0,P(x,y)<P(x,y-1)\end{array}\right.,(x,y)\&Element;R_0$

(2-3)对向上光滑性游程图的梯度图G₁进行横向求和，并对横向求和结果进行邻域相加处理，以得到数列S₂，具体为：

$S_1\left(y\right)=\underset{i=1}{\overset{X}{\mathrm{\&Sigma;}}}{G}_1(i,y),y\&Element;[1,Y]$

$S_2\left(y\right)=\left[\begin{array}{ccccccc}1& 2& 4& 6& 4& 2& 1\end{array}\right]\left[\begin{array}{c}{S}_1(y-3)\\ S_1(y-2)\\ S_1(y-1)\\ S_1\left(y\right)\\ S_1(y+1)\\ S_1(y+2)\\ S_1(y+3)\end{array}\right],y\&Element;[1,Y]$

其中S₁为所横向求和得到的数列，S₂为经过邻域求和处理得到的数列，X和Y分别为去噪后的归一化超声颈动脉纵截面图像的宽和高；

(2-4)将数列S₂中最大项的序号所对应的去噪后的归一化超声颈动脉纵截面图像中的行号作为感兴趣区域的中间横轴，从该中间横轴分别向上和向下扩展，以得到感兴趣区域R。

4.根据权利要求3所述的内中膜自动分割方法，其特征在于，向上扩展和向下扩展的高度可相同，也可不同，其取值是由去噪后的归一化超声颈动脉纵截面图像的尺寸、向上光滑性游程图中向上光滑性游程的最大值决定。

5.根据权利要求3所述的内中膜自动分割方法，其特征在于，步骤(3)具体包括以下子步骤：

(3-1)获取感兴趣区域R的梯度图G₂，具体为：

$G_2^{\&prime;}(x,y)=\left[\begin{array}{ccc}-1& 0.8& 0.2\end{array}\right]\left[\begin{array}{c}I(x,y)\\ I(x,y-1)\\ I(x,y-2)\end{array}\right],(x,y)\&Element;R$

$G_2(x,y)=\left\{\begin{array}{c}{G}_2^{\&prime;}(x,y),G_2^{\&prime;}(x,y)\&GreaterEqual;0\\ 0,G_2^{\&prime;}(x,y)<0\end{array}\right.,(x,y)\&Element;R$

其中G₂′为中间变量；

(3-2)根据感兴趣区域R的大小估计内中膜厚度的取值范围[D₁,D₂]；具体而言，该范围被设置为一个像素点至感兴趣区域R宽度的一半；

(3-3)依次取每一个整数值作为可能的内中膜厚度，根据该可能的内中膜厚度设定特征提取算子，并用该算子计算灰度梯度图中每个像素点的特征值，以得到与该可能的内中膜厚度对应的灰度梯度图的特征图F_d，具体为：

F_d(x,y)＝F_d1(x,y)×F_d2(x,y),(x,y)∈R

其中，

$F_{d1}(x,y)=\left\{\begin{array}{c}{G}_2(x,y+d)-G_2(x,y),G_2(x,y+d)\&GreaterEqual;G_2(x,y)\\ 0,G_2(x,y+d)<G_2(x,y)\end{array}\right.,(x,y)\&Element;R$

$F_{d2}(x,y)=\left\{\begin{array}{c}{G}_2(x,y-d)-G_2(x,y),G_2(x,y-d)\&GreaterEqual;G_2(x,y)\\ 0,G_2(x,y-d)<G_2(x,y)\end{array}\right.,(x,y)\&Element;R$

d为估计的内中膜厚度取值范围[D₁,D₂]中的一个整数值，F_d1和F_d2为中间变量；

(3-4)对获得的每个特征图的所有特征值进行求和处理，和S₃最大的特征图所对应的可能的内中膜厚度d*为内中膜厚度参考值，其中S₃具体为：

$S_3\left(d\right)=\underset{(x,y)\&Element;R}{\mathrm{\&Sigma;}}{F}_d(x,y),d\&Element;[D_1,D_2].$

6.根据权利要求5所述的内中膜自动分割方法，其特征在于，步骤(4)具体包括以下子步骤：

(4-1)对感兴趣区域进行对比度增强操作，以得到优化图像I₁，具体为：

$I_1(x,y)=1-\cos (I(x,y)\&times;\frac{\mathrm{\&pi;}}{2}),(x,y)\&Element;R$

(4-2)根据内中膜厚度参考值d^*获取感兴趣区域中像素点下方梯度的特征图G_down，具体为：

$G_{\mathrm{down}}(x,y)=\underset{i=d^{*}}{\overset{{3d}^{*}}{\mathrm{\&Sigma;}}}{G}_2(x,y-i),(x,y)\&Element;R$

(4-3)根据上述步骤(4-1)和(4-2)获得的结果构建能量函数E₁，具体为：

E₁＝w₁₁C₁₁+w₁₂C₁₂+w₁₃C₁₃+w₁₄C₁₄

其中，w₁₁、w₁₂、w₁₃、w₁₄为权重，C₁₁为轮廓线性度量函数，C₁₂为灰度项，C₁₃为灰度梯度项，C₁₄为轮廓下方区域梯度项；

(4-4)使用函数最小化方法使能量函数E₁最小化，即得到内中膜的内腔-内膜界面LII。

7.根据权利要求6所述的内中膜自动分割方法，其特征在于，步骤(5)具体包括以下子步骤：

(5-1)获取感兴趣区域中像素点与内中膜的内腔-内膜界面LII位置关系的特征图G_up，具体为：

$G_{\mathrm{up}}(x,y)={1000}^{\frac{|y+{2\&times;d}^{*}-\mathrm{LII}|\left(x\right)}{d^{*}}-1},(x,y)\&Element;R$

(5-2)根据上述步骤(5-1)获得的结果构建能量函数E₂，具体为：

E₂＝w₂₁C₂₁+w₂₂C₂₂+w₂₃C₂₃+w₂₄C₂₄

其中，w₂₁、w₂₂、w₂₃、w₂₄为权重，C₂₁为轮廓线性度项，C₂₂为灰度项，C₂₃为灰度梯度项，C₂₄为轮廓与LLI的位置相关项；

(5-3)使用函数最小化方法使能量函数E₂最小化，即得到内中膜的中膜-外膜界面MAI。