CN102163326A - 血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法 - Google Patents

Abstract

本发明属于计算机技术与医学图像的交叉领域，具体涉及一种血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法。具体步骤如下：载入原始血管超声图像；选择包含血管远侧壁颈动脉内中膜的ROI；人工定义或计算机自动检测初始轮廓线：包括计算单向边缘图、二值化与分段、Hough变换检测两条直线段PL_i和SL_i作为LII和MAI在该子图像中的部分、校正衔接，采用双Snake模型演化得到内中膜边界轮廓线，最后基于分割得到的两条边界计算厚度均匀度。本发明能有效应付图像中的噪声和边界不连续，达到精确分割内中膜边界的目标，减少医生工作量；同时基于本方法分析得到的厚度均匀度参数能为动脉粥样硬化的早期分析提供更多信息。

Description

血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法

技术领域

本发明属于计算机技术与医学图像的交叉领域，具体涉及一种血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法。

背景技术

缺血性心脑血管病是危害人类健康的重大疾病，它的主要病理基础是动脉粥样硬化。颈动脉血管内中膜厚度(Intima Media Thickness，IMT)反映动脉粥样硬化状况，它与心肌梗塞和脑卒中的发生率高度相关，因此，近二十年来其一直是预测心脑血管疾病风险的重要指标之一。

颈动脉血管内中膜厚度是通过在血管超声图像上测量获得。正常血管内中膜在声像图中表现为高回声和低回声两层带状区域，它处于两条边界即血管管腔与内膜边界(Lumen IntimaInterface，LII)和中膜与外膜边界(Media Adventitia Interface，MAI)之间，如图1所示。传统的血管超声图像中内中膜厚度的分析工作主要依靠操作者(检查医生)的肉眼观察和直接测量，即操作者选择自认为病变最严重的部位对内中膜厚度作人为的测量和估计，且对于多点内中膜平均厚度的测量还需要在多点进行。医生在这些测量活动中对于内中膜边界的判断直接取决于其所受的训练、对于边界的主观理解、经验等，因此不同医生判别的内中膜边界不同，甚至同一名医生在不同时间得到的结果也不完全一致，此外，其还有耗时长的缺点。在已公开的专利中，美国专利公开号为20030199762和专利公开号为20050119555的两个专利提出了自动检测边界的方法，其对梯度极大点使用曲线拟合估计边界，这种原始的方法忽视了超声图像中的噪声和边界不连续等众多伪像，只能应用于边界很清晰的情况。国内还没有这方面的专利公开。

内中膜的自动分割在近年来成为研究热点。现有的自动分割方法集中在用户选择的兴趣区(ROI)中自动检测这两条边界。研究的算法主要可分为2类：动态规划算法(如郭翌，汪源源，王威琪，超声图像颈动脉内中膜提取方法.声学技术，2007.26(005)：p.1034-1035.；Wendelhag，I.，et al.， A new automated computerized analyzing system simplifiesreadings and reduces the variability in ult

发明内容

本发明所要解决的技术问题是提供一种血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，该方法可自动检测血管超声图像中颈动脉血管内中膜边界，并能计算得到大范围颈动脉血管内中膜厚度的变化情况，为早期动脉粥样硬化病变提供分析参数。

为解决上述技术问题，本发明采用的技术方案为：一种血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于，包括以下步骤：

(1)在计算机中人工载入一幅采集的原始血管超声图像，见图1，所述的计算机可为一般的PC机；

(2)人工选择包含血管远侧壁颈动脉内中膜的一个矩形感兴趣区域(ROI)，见图1中的白色方框区域。血管远侧壁是指图像中远离探头的血管壁，通常在血管超声图像中管腔下方。在超声图像中因为声波向后散射，近侧壁不能形成清晰的内中膜结构，故内中膜厚度的测量通常选择在远侧壁进行；

(3)在矩形感兴趣区域(ROI)中检测演化得到颈动脉血管内中膜的LII(血管管腔与内膜边界)和MAI(中膜与外膜边界)两条边界的精确轮廓线：

(3.1)在矩形感兴趣区域(ROI)中检测得到颈动脉血管内中膜的LII和MAI两条边界的初始轮廓线；

(3.2)将步骤(3.1)检测得到的LII和MAI两条边界的初始轮廓线代入下述双Snake模型推导得到的差分方程进行反复迭代直至稳定，得到LII的精确轮廓线Γ₁(q)＝(x₁(q)，y₁(q))和MAI的精确轮廓线Γ₂(q)＝(x₂(q)，y₂(q))，所述的差分方程的推导过程如下：

定义双Snake模型的能量泛函为：

E_d-snake＝E_int+E_ext                  (1)

$E_{\mathrm{int}}=\frac{1}{2}\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\∫}_0^1{w}_{i1}{\left|\frac{\∂{\mathrm{\Γ}}_i}{\∂q}\right|}^2+w_{i2}{\left|\frac{{\∂}^2{\mathrm{\Γ}}_i}{\∂q^2}\right|}^2\mathrm{dq}---\left(2\right)$

$E_{\mathrm{ext}}=\frac{1}{2}\mathrm{\α}{\∫}_0^t{\left|\frac{\∂({\mathrm{\Γ}}_1-{\mathrm{\Γ}}_2)}{\∂q}\right|}^2\mathrm{dq}+(1-\mathrm{\α})\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\∫}_0^1-f\left({\mathrm{\Γ}}_i\right)\mathrm{dq}---\left(3\right)$

其中，Γ₁(q)＝(x₁(q)，y₁(q))表示LII的轮廓线，Γ₂(q)＝(x₂(q)，y₂(q))表示MAI的轮廓线，q表示曲线的参数，q∈[0，1]，q为0表示轮廓线的左端点，q为1表示轮廓线的右端点，q为中间值则表示轮廓线中间的一点。w_i1，i＝1，2分别是LII轮廓线的张力系数和MAI轮廓线的张力系数，w_i1，i＝1，2∈[0，1]，w_i2，i＝1，2分别是LII轮廓线的硬度系数和MAI轮廓线的硬度系数，w_i2，i＝1，2∈[0，1]，w_i2，i＝1，2越大，则对应的轮廓线越平滑，α调节轮廓线靠近图像边缘和两条轮廓线平行的权重，α∈[0，1]，α越小，轮廓线越靠近图像边缘，α越大，两条轮廓线越平行；E_int表示轮廓线的内部能量，E_int越小，则轮廓线越光滑。E_ext表示轮廓线的外部能量，E_ext越小，则轮廓线越靠近图像边缘并且轮廓线间的厚度越均匀，E表示总能量。

表示二条轮廓线的距离矢量的变化率，在双Snake模型E_ext中包含的能量项

使得二条轮廓线间的距离保持不变，这样使得LII和MAI两条轮廓线的演化不再是独立的，而是相互作用，互为调整，从而能够弥补超声图像中常有的不连续的边界。

使E_d-snake达到极小，根据梯度下降流可以得到：

$\frac{\∂x_i}{\∂t}=w_{i1}\frac{{\∂}^2{x}_i}{\∂q^2}-w_{i2}\frac{{\∂}^4{x}_i}{\∂q^4}+\mathrm{\α}(\frac{{\∂}^2{x}_i}{{\∂q}^2}-\frac{{\∂}^2{x}_{3-i}}{\∂q^2})+(1-\mathrm{\α})\frac{\∂f}{\∂x}(x_i,y_i)---\left(4\right)$

$\frac{\∂y_i}{\∂t}=w_{i1}\frac{{\∂}^2{y}_i}{\∂q^2}-w_{i2}\frac{{\∂}^4{y}_i}{\∂q^4}+\mathrm{\α}(\frac{{\∂}^2{y}_i}{{\∂q}^2}-\frac{{\∂}^2{y}_{3-i}}{\∂q^2})+(1-\mathrm{\α})\frac{\∂f}{\∂y}(x_i,y_i)---\left(5\right)$

其中，i＝1，2，对(4)和(5)中的4个偏微分方程在空间上使用中心差分，在时间上使用向前差分，得到差分方程；

(4)厚度均匀度分析计算：根据步骤(3)分割得到的颈动脉血管内中膜LII和MAI两条边界的精确轮廓线，计算厚度均匀度。

按上述方案，所述步骤(3.1)中的初始轮廓线检测方式为人工定义两条初始轮廓线，第一条靠近LII，第二条靠近MAI。人工定义的初始轮廓线应在正确的边界附近，以保证双Snake模型可使初始轮廓线收敛到正确的位置而精确地分割内中膜的两条边界。

按上述方案，所述的步骤(3.1)中的初始轮廓线检测方式为计算机自动检测，其步骤包括：

①计算单向边缘图；

②取T∈(0，1)，将单向边缘图二值化得到二值图像，然后将其沿x轴方向分为若干等宽且无重叠的子图像，图像中像素的列对应坐标系统的x坐标，见图2；

③对每个子图像上使用Hough变换检测PL_i和SL_i两条直线段作为LII和MAI在该子图像中的部分；

④找到检测到的直线段与其背景子图像匹配度最大的子图像，以它为参考，向左向右依次校正衔接错误的直线段，在进行下一次校正之前，采用前一次校正的直线段作为参考，校正邻接未经校正处理的子图像；

⑤将各子图像中的边界直线段连接得到颈动脉血管内中膜LII和MAI两条边界的初始轮廓线。

按上述方案，所述的步骤①的单向边缘图由下述单向边缘图定义计算得到：

单向边缘图定义为

$f(x,y)=\max \{0,\frac{\∂(G_{\mathrm{\σ}}*I(x,y))}{\∂y}\}---\left(6\right)$

其中，I(x，y)表示矩形感兴趣区域(ROI)图像。如图2所示：图像的左上角为原点，图像中像素的列对应坐标系统的x坐标，像素的行对应坐标系统的y坐标，

表示y方向的导数，

是一个2D高斯函数，σ是标准差，σ∈[0.5，2](以像素为单位，像素大小为0.066mm)，σ的取值与像素大小成反比例的变化，如像素大小扩大一倍，该参数的取值应缩小一半，像素大小为超声图像中每一个像素对应的实际正方形区域的边长，如像素大小为0.066mm表示每一像素对应实际的一个正方形区域中该区域的边长为0.066mm。σ较小，则单向边缘图越能精确反映正确的边界，但噪声也越多，σ较大，则噪声越少，但单向边缘图也越不能反映精确的边界。

按上述方案，所述的步骤②为：

a.将单向边缘图f(x，y)规范化到[0，1]区间内，取T＝0.3进行二值化，得到二值图像；

b.设选取的矩形感兴趣区域图像(ROI)的像素宽度为N，l_max是每段图像的可能划分的最大宽度(以像素为单位，像素大小为0.066mm)，把二值图像分为N_s个等宽无重叠的子图像{S_i|i＝1...N_s}，

l_max∈[50，150]，l_max的选取使每个分段子图像内的边界近似为直线，l_max的取值与像素大小成反比例的变化，如像素大小扩大一倍，该参数的取值应缩小一半。

按上述方案，所述的步骤③为：

a.使用Hough变换建立二维累计矩阵(Accumulator Array)；

b.从累计矩阵中找到最大值，根据其位置确定子图像中它所对应的一条边界的直线段PL_i，假设PL_i对应某一条颈动脉血管内中膜边界LII或MAI；

c.根据颈动脉血管内中膜的两条边界在分段子图像内近似平行，且两条边界的距离通常在一定的范围内，即PL_i和另一条边界的直线段SL_i之间的夹角θ应满足约束θ≤θ_max；且PL_i和SL_i之间的距离d应满足d_min≤d≤d_max，限制另一条边界所对应的局部最大值的搜索范围，在该范围内搜索最大值并根据其位置找到另一条边界的直线段，见图3，其中d_min为人工定义的原始血管超声图像中血管颈动脉内中膜厚度(IMT)的最小值，单位为像素，d_max为人工定义的原始血管超声图像中血管颈动脉内中膜厚度(IMT)的最大值，单位为像素，θ_max表示分段子图像中PL_i和SL_i的夹角的最大值，单位为度，d_min和d_max的取值与像素大小成反比例的变化。

按上述方案，所述步骤④为：

a.首先计算所有子图像{S_i|i＝1...N_s}中PL_i和SL_i主、次直线对与子图像的匹配度，找出所有子图像{S_i|i＝1...N_s}中匹配度最大的子图像S_j，匹配度计算见公式(7)：

{S_i|i＝1...N_s}中主、次直线对与子图像的匹配度V_i定义为

$V_i=\frac{1}{\left|{\mathrm{PL}}_i\right|}{\∫}_0^{\left|{\mathrm{PL}}_i\right|}{S}_i\left({\mathrm{PL}}_i\left(s\right)\right)\mathrm{ds}+\frac{1}{\left|{\mathrm{SL}}_i\right|}{\∫}_0^{\left|{\mathrm{SL}}_i\right|}{S}_i\left({\mathrm{SL}}_i\left(s\right)\right)\mathrm{ds}---\left(7\right)$

其中，|PL_i|和|SL_i|分别是PL_i，SL_i的长度。PL_i(s)，SL_i(s)是直线段PL_i，SL_i上的一点，S_i(PL_i(s))，S_i(SL_i(s))为该点处S_i子图像的取值；

b.从S_j出发，向左向右依次判断S_j-1和S_j+1子图像中的SL_j-1和SL_j+1是否与S_j衔接正确，设图4中S_i和S_i+1衔接处的4个端点a，b，c，d的纵坐标分别为y_a，y_b，y_c，y_d，将它们排序得到y′_a≤y′_b≤y′_c≤y′_d，如满足下式(8)，则判定S_i+1中的SL_i+1衔接出错。图4显示SL_i+1出错的3种情况：与PL_i+1靠得太近，太远或在错误的一边。

y′_c-y′_b＜max{y′_b-y′_a，y′_c-y′_b，y′_d-y′_c}                     (8)

如衔接错误，则加以校正，依此类推，直至最左和最右的子图像S₁和

在进行下一次衔接校正之前，将前一次校正的直线段作为参考，以校正邻接的未经校正处理的子图像，所述的校正方式为：假设分段子图像S_i中的PL_i和SL_i以及S_i+1中的PL_i+1对应真实边界，根据S_i中PL_i和SL_i的距离D_i将SL_i+1移至PL_i+1上方或下方使SL_i+1与PL_i+1距离为D_i的位置得到校正的SL_i+1，SL_i+1在PL_i+1上方或下方取决于PL_i+1离S_i中上方的直线近还是下方的直线近，如PL_i+1与下方的直线近，则SL_i+1在PL_i+1上方，否则SL_i+1在PL_i+1下方。

按上述方案，由于在血管超声图像中LII和MAI为近似水平方向的曲线，同一x坐标对应的只有唯一的LII和MAI的y坐标，为避免重新参数化，减少计算量，将所述步骤(3.2)中的偏微分方程可简化成新的方程(9)：

$\frac{{\∂y}_i}{\∂t}=w_{i1}\frac{{\∂}^2{y}_i}{\∂x^2}-w_{i2}\frac{{\∂}^4{y}_i}{\∂x^4}+\mathrm{\α}(\frac{{\∂}^2{y}_i}{\∂x^2}-\frac{{\∂}^2{y}_{3-i}}{\∂x^2})+(1-\mathrm{\α})\frac{\∂f}{\∂y}(x,y_i)---\left(9\right)$

其中，i＝1，2，y₁(x)表示LII轮廓线，y₂(x)表示MAI轮廓线，α∈[0，1]，α调整轮廓线靠近图像边缘和两条轮廓线平行的权重。

按上述方案，所述的步骤(4)为：在步骤(3)得到的颈动脉血管内中膜LII和MAI两条边界的轮廓线的每个横坐标位置插值得到点对(相对应的两点)，计算该点对的垂直距离d，同时计算该点对的平均值，得其平均点，然后进行拟合得到光滑曲线，估计该光滑曲线在点对处的局部斜率k，结合单点垂直距离d根据三角形学计算单点厚度为

根据单点厚度统计一段血管中内中膜各点厚度的均值、标准差、变异系数、百分构成比和百分位数。

与血管内中膜的传统分割算法相比，本发明提供的血管超声图像中颈动脉内中膜的自动分割和厚度均匀度分析方法，有两点区别：1.使用Hough变换检测初始轮廓线，考虑了所有边缘点，故不易受图像噪声的影响；2.引入基于形状约束的双Snake模型分割颈动脉血管内中膜的两条边界。该模型不再独立地使用Snake模型来变形演化检测单条内中膜边界，而是同时考虑颈动脉血管内中膜的两条边界，让两条轮廓线在自身张力和扭曲力以及外在图像的边缘力和相互间牵引力的作用下相互参考与修正，因此能够最大程度上克服超声图像中常有的边界不连续对测量的影响，提供精确的内中膜边界分割结果，从而为进一步计算分析内中膜参数提供方便。

本发明的有益效果：本发明提供的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法能有效地应付超声图像中的噪声和边界不连续，达到精确分割内中膜边界的目标，减少医生的工作量；基于本发明方法分析得到的内中膜厚度，平均值、标准差、变异系数、百分构成比和百分位数等反映厚度均匀度变化的参数能为动脉粥样硬化的早期分析提供更多的信息。

附图说明

图1为一幅内中膜轮廓线清晰的原始血管超声图像，图中白色ROI方框为包含血管远侧壁颈动脉内中膜的矩形感兴趣区域(ROI)。图2图1中ROI图像的坐标系统，其中原点在图像左上角，图像中像素的列对应坐标系统的x坐标，行对应坐标系统的y坐标。

图3为分段子图像中血管内中膜边界的检测。(a)中实线和虚线分别为对应的二值分段图像的PL和SL，PL位于LII上，SL位于MAI上。(b)显示(a)的检测过程。其中，中间位置的黑点为累计矩阵中最大的元素，其对应(a)中的PL。根据该元素的位置和内中膜的形状信息，设置距离约束和角度约束，在这些约束构成的矩形框内，找到下方的另一个最大元素即另一个黑点，其对应(a)中的SL。

图4为相邻两分段直线衔接出错的3种情况。(a)SL_i+1与PL_i+1距离太近；(b)SL_i+1与PL_i+1距离太远；(c)SL_i+1相对PL_i+1处在错误的一边。

图5为实施例1中计算机自动分割颈动脉血管内中膜的流程图。

图6显示实施例1的原始图像经过本发明的自动分割算法处理的全过程。(a)为原始的血管超声图像；(b)为计算得到的单向边缘图；(c)显示Hough变换检测的PL和SL；(d)显示校正后的PL和SL；(e)显示采样得到颈动脉血管内中膜的初始轮廓线；(f)显示双Snake模型演化后颈动脉血管内中膜LII和MAI边界的精确轮廓线。

图7显示实施例2中采用本发明的自动分割方法和人工描绘的内中膜边界的结果比较。(a)为一幅原始血管超声图像的ROI图像；(b)为(a)对应的采用本发明的自动分割方法得到的颈动脉血管内中膜边界的轮廓线；(c)为(a)对应的人工描绘的颈动脉血管内中膜边界的轮廓线。

图8显示实施例3中采用本发明的自动分割方法和人工描绘的内中膜边界的结果比较。(a)为一幅原始血管超声图像的ROI图像；(b)为(a)对应的采用本发明的自动分割方法得到的颈动脉血管内中膜边界的轮廓线；(c)为(a)对应的人工描绘的颈动脉血管内中膜边界的轮廓线。

图9显示实施例4中采用本发明的自动分割方法和人工描绘的内中膜边界的结果比较。(a)为一幅原始血管超声图像的ROI图像；(b)为(a)对应的采用本发明的自动分割方法得到的颈动脉血管内中膜边界的轮廓线；(c)为(a)对应的人工描绘的颈动脉血管内中膜边界的轮廓线。

具体实施方式

下面结合实施例对本发明提供的血管超声图像颈动脉血管内中膜的计算机自动分割算法和厚度均匀度分析方法作进一步详细的说明。

实施例1

一种血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于，包括以下步骤，如图2所示：

1、在与超声图像仪相连的计算机中载入一幅采集的原始血管超声图像，若血管超声图像像素大小不同时，计算机根据图像的像素大小(从DICOM格式中自动获取)自动调整具体的参数；

2、人工选择包含血管远侧壁内中膜的一个矩形感兴趣区域(ROI)，见图6中的图像(a)；

3、在矩形感兴趣区域中检测演化得到内中膜的LII(血管管腔与内膜边界)和MAI(中膜与外膜边界)两条边界的精确轮廓线：

3.1在矩形感兴趣区域中自动检测得到血管颈动脉内中膜LII和MAI两条边界的初始轮廓线，其步骤为；

①计算单向边缘图：

由于LII和MAI边界的上方灰度比下方小，本发明使用梯度方向朝下的边缘(下方灰度比上方大)，而不考虑梯度方向朝上的边缘。单向边缘图定义为

$f(x,y)=\max \{0,\frac{\∂(G_{\mathrm{\σ}}*I(x,y))}{\∂y}\}---\left(10\right)$

其中，I(x，y)表示ROI中的图像，如图2所示；坐标系统以图像左上角为原点，图像中像素的列对应坐标系统的x坐标，行对应坐标系统的y坐标，

表示y方向的导数，

G_σ是一个2D高斯函数，σ是标准差，σ＝1(以像素单位，像素大小为0.066mm)；

②二值化与分段：

在血管超声图像中，内中膜边界LII和MAI呈直线或分段直线，本发明对ROI图像二值化，然后分成若干等宽的子图像，以供后面的Hough变换处理，具体步骤为：

a.将单向边缘图f(x，y)规范化到[0，1]区间内，取T＝0.3进行二值化，得到二值图像；

b.设选取的ROI图像的像素宽度为N，l_max是每段图像的最大宽度(以像素为单位)，本实施例中l_max＝100(以像素为单位，像素大小为0.066mm)，把二值图像分为N_s个等宽无重叠的子图像{S_i|i＝1...N_s}，l_max的选取使得每个分段子图像内的边界近似为直线；

③Hough变换：

a.使用Hough变换建立二维累计矩阵(Accumulator Array)，具体参见Duda，R.O.and P.E.Hart，Use the Hough transformation to detect lines and curves in pictures.Comm.Assoc.Comput.，1972.15：p.11-15；

b.在累计矩阵中搜索最大值，该点对应着S_i子图像中一条直线段，此处称为主直线PL_i，假设PL_i对应某一条内中膜边界即LII或MAI；

c.定义本实施例中d_min＝5像素(像素大小为0.066mm)，d_max＝25像素(像素大小为0.066mm)，θ_max＝1°，根据该内中膜边界的位置，即内中膜的二条边界在分段子图像内近似平行，内中膜的两条边界间的距离通常在一定的范围内，则PL_i和SL_i之间的夹角θ应满足约束θ≤θ_max；且PL_i和SL_i之间的距离d应满足约束条件d_min≤d≤d_max，在累计矩阵中找到另一个边界对应的直线段，称其为次直线SL_i。由于Hough变换的二维累计矩阵中元素的位置的纵坐标和横坐标分别对应图像平面中的某一直线与原点的距离和方向(与坐标轴的角度)，所以上面的约束条件等价于在累计矩阵中最大值处(对应PL_i)的上方和下方的二个矩形区域内搜索最大值来定位SL_i。

④校正：

首先计算所有子图像中主、次直线对与子图像的匹配度，匹配度计算见公式(11)，找出所有子图像中匹配度最大的子图像S_j，从S_j出发，向左向右依次判断S_j-1和S_j+1子图像中的SL_j-1和SL_j+1是否与S_j衔接正确，如衔接错误，则加以校正，依此类推，直至最左和最右的子图像S₁和在进行下一次衔接校正之前，将前一次校正的直线段作为参考，以校正邻接的未经校正处理的子图像；

S_i中主、次直线对与子图像的匹配度V_i定义为

$V_i=\frac{1}{\left|{\mathrm{PL}}_i\right|}{\∫}_0^{\left|{\mathrm{PL}}_i\right|}{S}_i\left({\mathrm{PL}}_i\left(s\right)\right)\mathrm{ds}+\frac{1}{\left|{\mathrm{SL}}_i\right|}{\∫}_0^{\left|{\mathrm{SL}}_i\right|}{S}_i\left({\mathrm{SL}}_i\left(s\right)\right)\mathrm{ds}---\left(11\right)$

其中，|PL_i|和|SL_i|分别是PL_i，SL_i的长度。PL_i(s)，SL_i(s)是直线段PL_i，SL_i上的一点，S_i(PL_i(s))，S_i(SL_i(s))为该点处S_i子图像的取值。

假设分段子图像S_i中的PL_i和SL_i以及S_i+1中的PL_i+1对应真实边界，判断S_i+1中的SL_i+1是否衔接正确，只需考察S_i和S_i+1的衔接处4个端点的纵坐标关系即可，设S_i和S_i+1衔接处的4个端点a，b，c，d的纵坐标分别为y_a，y_b，y_c，y_d，将它们排序得到y′_a≤y′_b≤y′_c≤y′_d，如满足下式(12)，则判定S_i+1中的SL_i+1衔接出错。SL_i+1出错时，存在3种情况：与PL_i+1靠得太近，太远或在错误的一边。

y′_c-y′_b＜max{y′_b-y′_a，y′_c-y′_b，y′_d-y′_c}                      (12)

若S_i+1中的SL_i+1衔接无误，则以PL_i+1和SL_i+1处理下一子图像S_i+2，否则，根据S_i中PL_i和SL_i的距离D_i将SL_i+1移至PL_i+1上方或下方使SL_i+1与PL_i+1距离为D_i的位置得到校正的SL_i+1。SL_i+1在PL_i+1上方或下方取决于PL_i+1离S_i中上方的直线近还是下方的直线近。如PL_i+1与下方的直线近，则SL_i+1在PL_i+1上方，否则SL_i+1在PL_i+1下方。图6(d)显示(c)校正的结果。

⑤采样：将校正后的{PL_i|i＝1...N_s}和{SL_i|i＝1...N_s}首尾相连，按一定的采样间隔重新采样，建立初始轮廓线，其中偏上一条为LII的初始轮廓线，偏下一条为MAI的初始轮廓线。

3.2双Snake模型演化：

定义双Snake(d-snake)模型的能量泛函为

E_d-snake＝E_int+E_ext             (13)

上述公式中，Γ₁是关于q的函数Γ₁(q)，Γ₂是关于q的函数Γ₂(q)。q∈[0，1]表示轮廓线的参数，q为0表示轮廓线的左端点，q为1表示轮廓线的右端点，q为中间值则表示轮廓线中间的一点。Γ₁(q)＝(x₁(q)，y₁(q))表示LII的轮廓线，Γ₂(q)＝(x₂(q)，y₂(q))表示MAI的轮廓线。其中：w_i1，i＝1，2分别是LII轮廓线的张力系数和MAI轮廓线的张力系数，w_i2，i＝1，2分别是LII轮廓线的硬度系数和MAI轮廓线的硬度系数，w₁₁＝0.1，w₁₂＝0，w₂₁＝0.1，w₂₂＝0.1，α调节轮廓线靠近图像边缘和两条轮廓线平行的权重，图像边缘即图6(b)的单向边缘图的亮区域，α＝0.2；

使E_d-snake达到极小，根据梯度下降流可以得到：

$\frac{\∂x_i}{\∂t}=w_{i1}\frac{{\∂}^2{x}_i}{\∂q^2}-w_{i2}\frac{{\∂}^4{x}_i}{\∂q^4}+\mathrm{\α}(\frac{{\∂}^2{x}_i}{{\∂q}^2}-\frac{{\∂}^2{x}_{3-i}}{\∂q^2})+(1-\mathrm{\α})\frac{\∂f}{\∂x}(x_i,y_i)---\left(16\right)$

$\frac{\∂y_i}{\∂t}=w_{i1}\frac{{\∂}^2{y}_i}{\∂q^2}-w_{i2}\frac{{\∂}^4{y}_i}{\∂q^4}+\mathrm{\α}(\frac{{\∂}^2{y}_i}{{\∂q}^2}-\frac{{\∂}^2{y}_{3-i}}{\∂q^2})+(1-\mathrm{\α})\frac{\∂f}{\∂y}(x_i,y_i)---\left(17\right)$

其中，i＝1，2。由于在血管超声图像中LII和MAI为近似水平方向的曲线，同一x坐标对应的只有唯一的LII和MAI的y坐标，为避免重新参数化，减少计算量，将(16)和(17)的4个偏微分方程简化成(18)中的2个方程：

其中，i＝1，2，y₁(x)表示LII轮廓线，y₂(x)表示MAI轮廓线。在(18)中的3项中，第一项是光滑力，控制轮廓线的平滑，第二项是均匀力，控制LII和MAI两条轮廓线的厚度一致，第三项是边界力，控制轮廓线靠近图像边缘。本实施例对光滑力使用隐式差分格式，对均匀力和边界力使用显式差分格式，得到差分方程(19)：

$Y_i^{n+1}={(A_i+I_J)}^{-1}[Y_i^n+\mathrm{\αB}(Y_i^n-Y_{3-i}^n)+(1-\mathrm{\α})F_i^n]---\left(19\right)$

其中，表示LII(i＝1)或MAI(i＝2)轮廓线的离散化的J个点在n时间的纵坐标组成的向量(

表示初始轮廓线)，A_i是一个对角线元素分别为w_i2，-w_i1-4w_i2，2w_i1+6w_i2，-w_i1-4w_i2，w_i2的五对角矩阵(根据隐式差分格式得到)，w₁₁＝0.1，w₁₂＝0，w₂₁＝0.1，w₂₂＝0.1。B是一个对角线元素分别为1，-2，1的三对角矩阵(根据显式差分格式得到)，I_J是秩为J的单位矩阵，

表示轮廓线的离散化的J个点在n时间在力场f_y(x，y)中的取值组成的向量。总体上的这种半隐式差分格式避免了全隐式的计算复杂性，同时也改善了全显式的稳定性。

的绝对值被规范化到[0，1]区间内，以防止边界力产生不稳定。α＝0.2。在轮廓线的离散化中本实施例在x轴上每隔10个像素取一个点，在离散化点的边界上使用临近4个点的线性外插处理边界条件，保持边界处倾斜的内中膜轮廓；

将步骤3.1得到的颈动脉血管内中膜的两条边界的初始轮廓线代入差分方程进行反复迭代200次，得到精确的边界轮廓线LII轮廓线y₁(x)和MAI轮廓线y₂(x)。

4、厚度均匀度分析计算：

根据步骤3精确分割得到的颈动脉血管内中膜的LII和MAI两条边界轮廓线，在内中膜边界的两条轮廓线的每个横坐标位置插值得到点对(相对应的两点)，计算该点对的垂直距离d，同时计算该点对的平均值，得其平均点，然后进行拟合得到光滑曲线，估计该光滑曲线在点对处的局部斜率k，结合单点垂直距离d根据三角形学计算得：单点厚度为

再根据单点厚度统计一段血管中内中膜各点厚度的均值、标准差、变异系数、百分构成比和百分位数，结果见表1。

实施例2

根据实施例1的方法，对图7中的(a)图像进行自动分割得到颈动脉血管内中膜LII和MAI两条边界的精确轮廓线，见图(b)。图(c)为人工绘制的颈动脉血管内中膜LII和MAI两条边界的轮廓线。根据精确分割得到的颈动脉血管内中膜的LII和MAI两条边界轮廓线，统计一段血管中内中膜各点厚度的均值、标准差、变异系数、百分构成比和百分位数，结果见表1。该结果表明：本发明对于健康的内中膜(厚度较小)能分割的很好。

实施例3

根据实施例1的方法，对图8中的(a)图像进行自动分割得到颈动脉血管内中膜LII和MAI两条边界的精确轮廓线，见图(b)。图(c)为人工绘制的颈动脉血管内中膜LII和MAI两条边界的轮廓线。根据精确分割得到的颈动脉血管内中膜的LII和MAI两条边界轮廓线，统计一段血管中内中膜各点厚度的均值、标准差、变异系数、百分构成比和百分位数，结果见表1。该结果表明：本发明对于内中膜增厚但还未发生病变(厚度小于1mm)的情况也能正确分割。

实施例4

根据实施例1的方法，对图9中的(a)图像进行自动分割得到颈动脉血管内中膜LII和MAI两条边界的精确轮廓线，见图(b)。图(c)为人工绘制的颈动脉血管内中膜LII和MAI两条边界的轮廓线。根据精确分割得到的颈动脉血管内中膜的LII和MAI两条边界轮廓线，统计一段血管中内中膜各点厚度的均值、标准差、变异系数、百分构成比和百分位数，结果见表1。该结果表明：本发明对于内中膜增厚以致发生病变(厚度大于1mm)的情况也能正确分割。

表1

上表中：百分位数是第75百分位数，均值、标准差、百分位数和区间分布的区间单位为mm，变异系数是％。

Claims (9)

1.一种血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于，包括以下步骤：

(1)在计算机中人工载入一幅采集的原始血管超声图像；

(2)人工选择包含血管远侧壁颈动脉内中膜的一个矩形感兴趣区域；

(3)在矩形感兴趣区域中检测演化得到颈动脉血管内中膜的LII和MAI两条边界的精确轮廓线：

(3.1)在矩形感兴趣区域中检测得到颈动脉血管内中膜的LII和MAI两条边界的初始轮廓线；

(3.2)将步骤(3.1)检测得到的LII和MAI两条边界的初始轮廓线代入下述双Snake模型推导得到的差分方程进行反复迭代直至稳定，得到LII的精确轮廓线Γ₁(q)＝(x₁(q)，y₁(q))和MAI的精确轮廓线Γ₂(q)＝(x₂(q)，y₂(q))，所述的差分方程的推导过程如下：

定义双Snake模型的能量泛函为：

E_d-snake＝E_int+E_ext                      (1)

$E_{\mathrm{int}}=\frac{1}{2}\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\∫}_0^1{w}_{i1}{\left|\frac{\∂{\mathrm{\Γ}}_i}{\∂q}\right|}^2+w_{i2}{\left|\frac{{\∂}^2{\mathrm{\Γ}}_i}{\∂q^2}\right|}^2\mathrm{dq}---\left(2\right)$

$E_{\mathrm{ext}}=\frac{1}{2}\mathrm{\α}{\∫}_0^t{\left|\frac{\∂({\mathrm{\Γ}}_1-{\mathrm{\Γ}}_2)}{\∂q}\right|}^2\mathrm{dq}+(1-\mathrm{\α})\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\∫}_0^1-f\left({\mathrm{\Γ}}_i\right)\mathrm{dq}---\left(3\right)$

其中，Γ₁(q)＝(x₁(q)，y₁(q))表示LII的轮廓线，Γ₂(q)＝(x₂(q)，y₂(q))表示MAI的轮廓线，q表示曲线的参数，q∈[0，1]，q为0表示轮廓线的左端点，q为1表示轮廓线的右端点，q为中间值则表示轮廓线中间的一点，w_i1，i＝1，2分别是LII轮廓线的张力系数和MAI轮廓线的张力系数，w_i1，i＝1，2∈[0，1]，w_i2，i＝1，2分别是LII轮廓线的硬度系数和MAI轮廓线的硬度系数，w_i2，i＝1，2∈[0，1]，α∈[0，1]，使E_d-snake达到极小，根据梯度下降流可以得到：

$\frac{\∂x_i}{\∂t}=w_{i1}\frac{{\∂}^2{x}_i}{\∂q^2}-w_{i2}\frac{{\∂}^4{x}_i}{\∂q^4}+\mathrm{\α}(\frac{{\∂}^2{x}_i}{{\∂q}^2}-\frac{{\∂}^2{x}_{3-i}}{\∂q^2})+(1-\mathrm{\α})\frac{\∂f}{\∂x}(x_i,y_i)---\left(4\right)$

$\frac{\∂y_i}{\∂t}=w_{i1}\frac{{\∂}^2{y}_i}{\∂q^2}-w_{i2}\frac{{\∂}^4{y}_i}{\∂q^4}+\mathrm{\α}(\frac{{\∂}^2{y}_i}{{\∂q}^2}-\frac{{\∂}^2{y}_{3-i}}{\∂q^2})+(1-\mathrm{\α})\frac{\∂f}{\∂y}(x_i,y_i)---\left(5\right)$

其中，i＝1，2，对(4)和(5)中的4个偏微分方程在空间上使用中心差分，在时间上使用向前差分，得到差分方程；

(4)厚度均匀度分析计算：根据步骤(3)分割得到的颈动脉血管内中膜LII和MAI两条边界的精确轮廓线，计算厚度均匀度。

2.根据权利要求1所述的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于：所述步骤(3.1)中的初始轮廓线检测方式为人工定义两条初始轮廓线，第一条靠近LII，第二条靠近MAI。

3.根据权利要求1所述的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于：所述的步骤(3.1)中的初始轮廓线检测方式为计算机自动检测，其步骤包括：

①计算单向边缘图；

②取T∈(0，1)，将单向边缘图二值化得到二值图像，然后将其沿x轴方向分为若干等宽且无重叠的子图像，图像中像素的列对应坐标系统的x坐标；

③对每个子图像上使用Hough变换检测PL_i和SL_i两条直线段作为LII和MAI在该子图像中的部分；

④找到检测到的直线段与其背景子图像匹配度最大的子图像，以它为参考，向左向右依次校正衔接错误的直线段，在进行下一次校正之前，采用前一次校正的直线段作为参考，校正邻接未经校正处理的子图像；

⑤将各子图像中的边界直线段连接得到颈动脉血管内中膜LII和MAI两条边界的初始轮廓线。

4.根据权利要求3所述的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于：所述的步骤①的单向边缘图由下述单向边缘图定义计算得到：

单向边缘图定义为

$f(x,y)=\max \{0,\frac{\∂(G_{\mathrm{\σ}}*I(x,y))}{\∂y}\}---\left(6\right)$

其中，I(x，y)表示矩形感兴趣区域(ROI)图像，图像的左上角为原点，图像中像素的列对应坐标系统的x坐标，像素的行对应坐标系统的y坐标，

表示y方向的导数，

是一个2D高斯函数，σ是标准差，σ∈[0.5，2](以像素为单位)，σ的取值与像素大小成反比例的变化。

5.根据权利要求3所述的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于：所述的步骤②为：

a.将单向边缘图f(x，y)规范化到[0，1]区间内，取T＝0.3进行二值化，得到二值图像；

b.设选取的矩形感兴趣区域图像(ROI)的像素宽度为N，l_max是每段图像的可能划分的最大宽度(以像素为单位)，把二值图像分为N_s个等宽无重叠的子图像{S_i|i＝1...N_s}，

l_max∈[50，150]，l_max的选取使每个分段子图像内的边界近似为直线，l_max的取值与像素大小成反比例的变化。

6.根据权利要求3所述的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于：所述的步骤③为：

a.使用Hough变换建立二维累计矩阵；

b.从累计矩阵中找到最大值，根据其位置确定子图像中它所对应的一条边界的直线段PL_i，假设PL_i对应某一条颈动脉血管内中膜边界LII或MAI；

c.根据颈动脉血管内中膜的两条边界在分段子图像内近似平行，且两条边界的距离通常在一定的范围内，即PL_i和另一条边界的直线段SL_i之间的夹角θ应满足约束θ≤θ_max；且PL_i和SL_i之间的距离d应满足d_min≤d≤d_max，限制另一条边界所对应的局部最大值的搜索范围，在该范围内搜索最大值并根据其位置找到另一条边界的直线段，其中d_min为人工定义的原始血管超声图像中血管颈动脉内中膜厚度(IMT)的最小值，单位为像素，d_max为人工定义的原始血管超声图像中血管颈动脉内中膜厚度(IMT)的最大值，单位为像素，θ_max表示分段子图像中PL_i和SL_i的夹角的最大值，单位为度，d_min和d_max的取值与像素大小成反比例的变化。

7.根据权利要求3所述的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于：所述步骤④为：

a.首先计算所有子图像{S_i|i＝1...N_s}中PL_i和SL_i主、次直线对与子图像的匹配度，找出所有子图像{S_i|i＝1...N_s}中匹配度最大的子图像S_j，匹配度计算见公式(7)：

{S_i|i＝1...N_s}中主、次直线对与子图像的匹配度V_i定义为

$V_i=\frac{1}{\left|{\mathrm{PL}}_i\right|}{\∫}_0^{\left|{\mathrm{PL}}_i\right|}{S}_i\left({\mathrm{PL}}_i\left(s\right)\right)\mathrm{ds}+\frac{1}{\left|{\mathrm{SL}}_i\right|}{\∫}_0^{\left|{\mathrm{SL}}_i\right|}{S}_i\left({\mathrm{SL}}_i\left(s\right)\right)\mathrm{ds}---\left(7\right)$

其中，|PL_i|和|SL_i|分别是PL_i，SL_i的长度。PL_i(s)，SL_i(s)是直线段PL_i，SL_i上的一点，S_i(PL_i(s))，S_i(SL_i(s))为该点处S_i子图像的取值；

b.从S_j出发，向左向右依次判断S_j-1和S_j+1子图像中的SL_j-1和SL_j+1是否与S_j衔接正确，如衔接错误，则加以校正，依此类推，直至最左和最右的子图像S₁和在进行下一次衔接校正之前，将前一次校正的直线段作为参考，以校正邻接的未经校正处理的子图像，所述的校正方式为：假设分段子图像S_i中的PL_i和SL_i以及S_i+1中的PL_i+1对应真实边界，根据S_i中PL_i和SL_i的距离D_i将SL_i+1移至PS_i+1上方或下方使SL_i+1与PL_i+1距离为D_i的位置得到校正的SL_i+1，SL_i+1在PL_i+1上方或下方取决于PL_i+1离S_i中上方的直线近还是下方的直线近，如SL_i+1与下方的直线近，则SL_i+1在PL_i+1上方，否则SL_i+1在PL_i+1下方。

8.根据权利要求1所述的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于：将所述步骤(3.2)中的偏微分方程可简化成新的方程(8)：

$\frac{{\∂y}_i}{\∂t}=w_{i1}\frac{{\∂}^2{y}_i}{\∂x^2}-w_{i2}\frac{{\∂}^4{y}_i}{\∂x^4}+\mathrm{\α}(\frac{{\∂}^2{y}_i}{\∂x^2}-\frac{{\∂}^2{y}_{3-i}}{\∂x^2})+(1-\mathrm{\α})\frac{\∂f}{\∂y}(x,y_i)---\left(9\right)$

其中，i＝1，2，y₁(x)表示LII轮廓线，y₂(x)表示MAI轮廓线，α∈[0，1]。

9.根据权利要求1所述的血管超声图像中颈动脉血管内中膜的计算机自动分割和厚度均匀度分析方法，其特征在于：所述的步骤(4)为：在步骤(3)得到的颈动脉血管内中膜LII和MAI两条边界的轮廓线的每个横坐标位置插值得到点对，计算该点对的垂直距离d，同时计算该点对的平均值，得其平均点，然后进行拟合得到光滑曲线，估计该光滑曲线在点对处的局部斜率k，结合单点垂直距离d根据三角形学计算单点厚度为根据单点厚度统计一段血管中内中膜各点厚度的均值、标准差、变异系数、百分构成比和百分位数。

Images