颈动脉超声图像内中膜厚度测量方法及系统与流程

本发明涉及颈动脉超声图像处理领域，尤其是颈动脉超声图像的内中膜厚度测量算法。具体讲,涉及颈动脉超声图像内中膜厚度测量方法及系统。

背景技术：

超声成像中颈动脉血管壁的内中膜厚度(Intima Media Thickness,IMT)可以作为评估心血管疾病的早期病变程度的重要指标，对突发心肌梗塞、中风的诊断和预防有很重要的价值。传统上，颈动脉超声图像的内中膜厚度通常由人工手动测量获得，测量者手工描绘图像中的管腔-内膜边界(Lumen-Intima Interface,LII)和中膜-外膜边界(Media-Adventitia Interface,MAI)，然后通过计算得出两边界间的距离获得IMT。然而这种测量方式工作量大，非常耗时并且最终获得的结果受测试者所受训练、个人经历及其主观判断的影响。手动测量的缺陷促进了计算机辅助的IMT测量方法的发展。本发明提出了一种结合水平集分割和马尔可夫随机场(Markov Random Field,MRF)模型的颈动脉超声图像内中膜厚度测量算法。

基于几何变形模型的水平集方法由Osher和Sethian于1982年提出。该方法是利用偏微分方程作为数值分析放过与技术手段被广泛运用于轮廓面或轮廓线的运动跟踪。该方法将低维闭合曲线演化的问题转化为高维空间中谁破解函数全面演化的隐含方式来求解，因此可以适应于拓扑结构变化的处理。水平集方法的思想中最重要的就是Hamilton-Jacobi方程，它为运动的隐式曲面进行基于时间的方程的数值求解。该方法的基本原理是，将曲线或曲面作为零水平集嵌入到高一维的水平集函数中，通过一个更高维的函数来表达低维的曲线或曲面的演化过程。零水平集被嵌入到水平集函数中，只要确定零水平集位置，就可以确定运动曲线或曲面演化的结果。下面简要分析曲线模型的演化过程。平面闭合曲线C(s,t)＝(m(s,t),n(s,t)):0≤s≤1(s为弧长参数，t表示时间，m、n是平面曲线的坐标)被水平集方法隐含地表达为三维连续函数曲面φ(m,n,t):的一个具有相同函数值的等值曲线，令函数φ为符号距离函数，可将曲线C的变化归结为由于函数发生了某种相应的变化所致。因此，可将随时间变化的闭合曲线表示为水平集函数对时间的变化：

上式把对曲线的演化转换为对水平集函数的分析。

可用如下的偏微分方程来表示曲线C沿法线方向的演化过程：

上式被称为水平集方程(Level Set Equation)，函数V被称为速度函数(Speed Function)，表示曲线上个点的演化速度，N是曲线的法线方向。对式(2)求取全微分，得：

式中是φ的梯度，与曲线的法向方向相同。假设Ω为整幅图像，水平集函数位于曲线C以内的部分为负，以外的部分为正，则水平集的内向单位法向量为

将式(2)和式(4)代入式(3)中，整理可得

这样就用水平集方程表示了曲线演化过程，即曲线演化方程的欧氏表达，它是一种Hamilton-Jacobi类型的偏微分方程。给定几何活动轮廓模型的偏微分方程，以及初始的水平集函数φ(m,n):φ(C)＝0，式(5)就可以保证水平集函数φ(m,n,t)随时间的演化满足{φ(C(s,t),t)＝0}条件，即φ(m,n,t)的零水平集始终是轮廓线C(s,t)。

技术实现要素：

为克服现有技术的不足，本发明旨在提出一种可靠、自动颈动脉超声图像的测量方法，使用计算机辅助测量的方式避免手动测量方式存在的缺陷。本发明采用的技术方案是，颈动脉超声图像内中膜厚度测量方法，步骤如下：

1)图像裁剪，裁剪移除无关信息，只保留血管部分；

2)自动提取感兴趣区域ROI(Region of Interest)；

3)图像滤波，滤除回波反射带来的斑块噪声；

4)水平集分割，采用水平集分割算法追踪血管壁的管腔-内膜边界LII(Lumen-Intima Interface,)；

5)初始标记场，根据颈动脉血管壁各层近似平行的特性，在ROI中，将血管壁内膜边界向下平移一定像素距离，作为中膜-外膜边界中膜-外膜边界MAI(Media-Adventitia Interface,)的初始状态，然后依据这两条边界将ROI分割为管腔、内中膜、外膜三个部分，并将获得的分割图像作为Markov模型的初始标号场；

6)再分割，引入Markov模型，根据初始标号场和灰度场对图像再次分割图像；

7)后处理，对分割图像后处理，移除错分小区域，获得最终分割图像和LII、MAI，并测量内中膜厚度IMT(Intima Media Thickness)。

步骤2)具体步骤是：根据超声图像灰度的统计规律，采用最大类间差法自适应获取灰度阈值，二值化超声图像；然后，对二值化超声图像形态学闭运算，再填补图像中的孔洞，利用连通域面积较大以及质心纵坐标较大的原则，选择期望的白色区域，其与黑色区域相交的上边界认为是LII的近似位置；最后以LII的近似位置为标准，向上取40、向下取20像素尺度，分别作为ROI的上、下边界，获得两个分别对应原始图像ROI。

步骤3)图像滤波是基于像素间距离相似性和灰度相似性的滤波，滤波输出为

其中g(r)是滤波后图像，是像素ξ和像素r的相似性函数，为归一化函数，σ₁是定义域方差，σ₂是值域方差；

步骤4)进一步及细化为：

Step1：初始化，假设t_i时刻的水平集函数φ(m,n,t_i)为已知，m、n是平面曲线的坐标；

Step2：求解水平集方程：获得t_i+1时刻的水平集函数φ(m,n,t_i+1)在整个图像区域的值，此时记为Γ_i+1的演化曲线就是φ(m,n,t_i+1)的零等值曲线，即Γ_i+1＝{φ(m,n,t_i+1)＝0}，此时的φ(m,n,t_i+1)不再是符号距离函数；

Step3：重新初始化，使用φ(m,n,t_i+1)代替下式中的初始水平集函数φ₀，迭代求解至稳定解，仍记为φ(m,n,t_i+1)：

其中φ_t为t时刻的水平集函数，φ(m,n)为初始的水平集函数；

Step4：结合φ(m,n,t_i+1)的值，求解物理量的控制方程，得到t_i+1时刻的物理量的值；

Step5：重复步骤，进入下一时刻的计算，直至停止。

步骤6)建立待分割图像的Markov模型，对于Markov模型的求解，流程如下：

Step1:输入待分割的原始图像；

Step2:对图像进行初始化，得到图像的初始标记场；

Step3:在上一步标号场的基础上对图像中每个像素遍历求使之取到满足迭代条件的标记，更新原有的标号场；

对Step3进行迭代直到收敛条件。

步骤6)进一步细化为：

图像分割是基于像素的特征属性和区域属性给每一个像素分配标号的过程。在马尔可夫随机场模型(Markov Random Field,MRF)中，用两个随机场来描述待分割的图像，一个是标号场X，常称为隐随机场，用先验分布P(X)描述标号场的局部相关性；另一个是灰度场或特征场Y，以标号场为条件，用分布函数P(Y|X)描述观测数据或特征向量的分布，其中标号场X是马尔可夫场，满足Markov特性，其先验分布表示为

其中A＝{1,…a}是图像中像素的集合，a是图中像素总数，x_i∈{1,…,q}是像素i的类别标号，q是图像中像素类别的总数，N_i为像素i的邻域，根据Hammersley-Clifford等价定理，所有像素的联合概率分布为

其中Z是归一化常数，U(x)是先验能量

c是一个基团，即包含若干位置的集合，x_i、x_j是像素点i、j的类别标号，极端情况下认为邻域结构的所有子集都是基团，C表示基团的集合，V_c(x)是定义在基团c上的势函数，它仅与像素的邻域像素有关，式中V₁(x_i)对应于邻域相关函数，势团模型选择二阶的Potts模型进行估计，即仅考虑二阶邻域系统中两点间的相互作用，则用伪似然近似逼近的先验概率表达为：

Potts模型：

其中β是控制邻域间作用强度的函数，N_i、n_i是像素i的邻域像素；

另一个随机场Y指待分割图像，它是可观测的随机过程，即P(Y)是固定的，则根据贝叶斯定理：

P(X|Y)∝P(X)P(Y|X) (13)

将模型转化为对似然函数P(Y|X)和先验概率P(X)的求解，这里P(X|Y)是图像的后验概率，其中条件分布函数P(Y|X)，通常用高斯分布描述：

其中为相同标记为x_i的区域像素灰度的方差，为像素i邻域灰度的均值；

引入Markov模型，根据公式(11)、公式(14)求解模型。

选取Potts模型作为公式(10)中的势团表达函数，采用条件迭代算法ICM(Iterative Conditional Mode)求解Markov模型。

颈动脉超声图像内中膜厚度测量系统，由超声波装置及计算机组成，超声波装置测得的图像信息发送给计算机处理，计算机设置有如下模块：

1)图像裁剪模块，裁剪移除无关信息，只保留血管部分；

2)自动提取感兴趣区域ROI(Region of Interest)模块；

3)图像滤波模块，滤除回波反射带来的斑块噪声；

4)水平集分割模块，采用水平集分割算法追踪血管壁的管腔-内膜边界LII(Lumen-Intima Interface,)；

5)初始标记场模块，根据颈动脉血管壁各层近似平行的特性，在ROI中，将血管壁内膜边界向下平移一定像素距离，作为中膜-外膜边界中膜-外膜边界MAI(Media-Adventitia Interface,)的初始状态，然后依据这两条边界将ROI分割为管腔、内中膜、外膜三个部分，并将获得的分割图像作为Markov模型的初始标号场；

6)再分割模块，引入Markov模型，根据初始标号场和灰度场对图像再次分割图像；

7)后处理模块，对分割图像后处理，移除错分小区域，获得最终分割图像和LII、MAI，并测量内中膜厚度IMT(Intima Media Thickness)。

本发明的特点及有益效果是：

水平集算法的本质是用高一维的曲面的零等高线作为目标边界去分割低一维的目标。如果从高维度看低纬度，低纬度的拓扑变化在高纬度中的表现为曲面的形态变化，不会造成曲面的拓扑结构变化。所以水平集算法具备相当强的低纬度的拓扑可变性。本发明引入水平集方法作为ROI的初始分割算法，方法稳定，有助于可靠结果的获得。Markov模型可以有效的刻画图像中邻域像素将的空间信息，能够使计算机的测量结果更加准确。本发明有力的支持了临床中颈动脉超声图像内中膜厚度的测量，为IMT计算机辅助测量技术的进一步优化发展提供了参考，对专家手动测量的方式是很好的补充。

附图说明：

图1算法流程图；

图2初始颈动脉超声图像；

图3裁剪后图像；

图4 ROI；

图5滤波后图像；

图6初始LII边界；

图7初始LII、MAI边界；

图8初始标号场；

图9迭代后标号场；

图10最终分割图像；

图11最终LII、MAI边界。

具体实施方式

本发明的目的就是为了提出一种可靠、自动颈动脉超声图像的测量算法，使用计算机辅助测量的方式避免手动测量方式存在的缺陷。基于此目的，针对目前的国内外研究现状，在本发明中采用水平集算法初步分割测量图像，并引入Markov模型，将图像分割问题转化为每一个像素分配标号的过程。

本发明结合超声颈动脉图像的特点，在所给图像库的测试过程中，该发明能有效的完成分割颈动脉图像并测量IMT的工作，具有较好的理论和使用价值。

为了实现上述目的，本发明采用如下的技术方案：

1)图像裁剪。初始超声图像的周围分布有病人和超声仪器的相关信息，这些信息的存在会影响后续的图像处理步骤，因此需要裁剪移除这些部分，只保留血管部分。

2)感兴趣区域(Region of Interest,ROI)提取。自动提取ROI可以避免手动分割ROI的繁琐，实现自动测量IMT的目的。

3)图像滤波。超声图像以实时性、可重复性、非侵入性及成本低的优点，成为颈动脉检查的首选成像方式。然而，回波反射带来的斑块噪声不仅会减低图像的成像质量，更增加了计算机处理的难度，通过对图像进行滤波，降低噪声的影响，获得更加可信的结果。

4)水平集分割。采用水平集分割算法追踪血管壁的管腔-内膜边界((Lumen-Intima Interface,LII))。

5)初始标号场。根据颈动脉血管壁各层近似平行的特性，在ROI中，将血管壁内膜边界向下平移一定像素距离，作为中膜-外膜边界中膜-外膜边界(Media-Adventitia Interface,MAI)的初始状态。然后依据这两条边界将ROI分割为管腔、内中膜、外膜三个部分，并将获得的分割图像作为Markov模型的初始标号场。

6)再分割。引入Markov模型，根据初始标号场和灰度场对图像再次分割图像。

7)后处理。对分割图像后处理，移除错分小区域，获得最终分割图像和LII、MAI，并测量IMT。

下面结合附图与实施例对本发明作进一步说明。

1)图像裁剪。在保证远端血管壁存在的条件下，截取超声图像中下方320×300像素大小的区域。

2)感兴趣区域(Region of Interest,ROI)提取。

自动选取阈值并对其二值化。理想情况下，超声图像中血管的管腔和血管壁的灰度分别呈现出暗、亮的特点，因此根据图像灰度的统计规律，采用最大类间差法自适应获取灰度阈值，二值化图像。

寻找远端血管壁的近似位置。超声图像中血管位置的不同、斑块噪声的存在，使二值化的图像可能存在多个白色区域。为了寻找与远端血管壁对应的白色区域，对二值图像形态学闭运算，填补图像中的孔洞，并根据远端血管壁位于超声图像下方的特点，利用连通域面积较大以及质心纵坐标较大的原则，选择期望的白色区域，其与黑色区域相交的上边界可认为是LII的近似位置。

提取ROI。正常颈动脉血管的IMT约在0.5mm到1mm的范围内，对应超声图像中约8-16个像素点。考虑到血管可能存在弯曲和病变的情况，我们以LII的近似位置为标准，向上取40、向下取20像素尺度，分别作为ROI的上、下边界。在这一步骤中可以获得两个分别对应原始图像ROI。

3)图像滤波。

本发明选取双边滤波算法分别对两个ROI图像滤波。双边滤波算法是一种非线性滤波方法，计算量小，它在医学图像的去噪处理中有较好的表现，并且在去除噪声的同时能够较好地保持图像的边缘信息。该算法同时考虑了图像像素间的距离相似性和灰度值相似性，并且可以通过高斯分布来刻画。对图像f(r)，基于像素间距离相似性和灰度相似性的滤波输出为：

其中g(r)是滤波后图像，是像素ξ和r的相似性函数，为归一化函数。σ₁是定义域方差，σ₂是值域方差；

4)水平集分割。

Step1：初始化，假设t_i时刻的水平集函数φ(m,n,t_i)为已知，m、n是平面曲线的坐标；

Step2：求解水平集方程：获得t_i+1时刻的水平集函数φ(m,n,t_i+1)在整个图像区域的值，此时记为Γ_i+1的演化曲线就是φ(m,n,t_i+1)的零等值曲线，即Γ_i+1＝{φ(m,n,t_i+1)＝0}，此时的φ(m,n,t_i+1)不再是符号距离函数；

Step3：重新初始化，使用φ(m,n,t_i+1)代替下式中的初始水平集函数φ₀，迭代求解至稳定解，仍记为φ(m,n,t_i+1)：

其中φ_ti+1为t_i+1时刻的水平集函数。

Step4：结合φ(m,n,t_i+1)的值，求解物理量的控制方程，得到t_i+1时刻的物理量的值；

Step5：重复步骤，进入下一时刻的计算，直至停止。

5)初始标号场。

建立初始标号场。根据ROI图像垂直方向边缘图，将初始LII向下平移一定距离(水平方向上平均梯度值最大位置)，作为MAI的初始标记。ROI图像被LII和MAI划分为三部分，由上至下依次是：管腔(标记为1)、内中膜(标记为2)和外膜(标记为3)。

6)再分割。

图像分割是基于像素的特征属性和区域属性给每一个像素分配标号的过程。在马尔可夫随机场模型(Markov Random Field,MRF)中，常用两个随机场来描述待分割的图像，一个是标号场X，常称为隐随机场，用先验分布P(X)描述标号场的局部相关性。另一个是灰度场或特征场Y，常以标号场为条件，用分布函数P(Y|X)描述观测数据或特征向量的分布。其中标号场X是马尔可夫场，满足Markov特性，其先验分布可表示为

其中A＝{1,…a}是图像中像素的集合，a是图中像素总数，x_i∈{1,…,q}是像素i的类别标号，q是图像中像素类别的总数，N_i为像素i的邻域。根据Hammersley-Clifford等价定理，所有像素的联合概率分布为

其中Z是归一化常数，U(x)是先验能量

c是一个基团，即包含若干位置的集合，x_i、x_j是像素点i、j的类别标号，极端情况下可以认为邻域结构的所有子集都是基团，C表示基团的集合。V_c(x)是定义在基团c上的势函数，它仅与像素的邻域像素有关，例式(10)中V₁(x_i)对应于邻域相关函数，势团模型选择二阶的Potts模型进行估计，即仅考虑二阶邻域系统中两点间的相互作用。则用伪似然近似逼近的先验概率表达为：

Potts模型：

其中β是控制邻域间作用强度的函数，N_i、n_i是像素i的邻域像素。

另一个随机场Y指待分割图像，它是可观测的随机过程，即P(Y)是固定的。则根据贝叶斯定理：

P(X|Y)∝P(X)P(Y|X) (13)

将模型转化为对似然函数P(Y|X)和先验概率P(X)的求解。这里P(X|Y)是图像的后验概率，其中条件分布函数P(Y|X)，通常用高斯分布描述：

其中为相同标记(标记为x_i)的区域像素灰度的方差，为像素i邻域灰度的均值。

引入Markov模型，根据公式(11)、公式(14)求解模型。选取Potts模型作为公式(10)中的势团表达函数。本发明中采用条件迭代算法(Iterative Conditional Mode,ICM)求解Markov模型。ICM是一种确定性算法，能够获得局部最优解。算法在每一步的迭代过程中，都要求达到最大化以便求解过程中能过快速地收敛达到局部最优解。ICM的计算量大大减少。

7)后处理。

我们采用最大连通域准则进行处理，移除离散小区域，得到最终分割图像。用Sobel算子提取分割图像中第二部分内中膜(标记为2)的边界，获取LII和MAI的最终轮廓。并计算边界间的平均距离，作为计算机辅助的颈动脉IMT距离。