一种颈动脉超声图像内中膜自动分割方法与流程

本发明涉及一种颈动脉超声图像内中膜自动分割方法，属于医学图像处理领域。

背景技术：

心血管疾病是威胁人类健康的首要疾病。动脉粥样硬化会造成动脉血管壁增厚，造成管腔狭窄，因而是心血管疾病发生的主要原因。血管壁从内到外分为外膜层、中膜层和内膜层三个相互紧贴的膜层。内中膜厚度(intima-mediathickness,imt)是指为从管腔-内膜(lumen-intima,li)边界到中膜-外膜(media-adventitia,ma)边界的距离，即内膜和中膜两部分厚度相加。imt的增加是动脉粥样硬化的一个早期临床症状。大量研究证明，颈动脉imt和心脑血管疾病存在显著相关性，并且可以作为未来心血管事件的一个强预测指标。如今，颈动脉内中膜厚度已被看作是初步判别颈动脉粥样硬化程度和心血管病变情况的一个重要的指标。

医学超声成像技术可以对颈动脉进行较清晰的成像，是一种有效的颈动脉粥样硬化普查手段。借助于超声成像技术，可以对颈动脉内中膜进行成像，进而可以分割出内中膜的边界，实现imt的测量。然而，临床采用的由医生手动分割来获取内中膜边界的方法费时、枯燥，并且会存在较大的观察者间的差异，因此迫切需要一种能够自动分割内中膜边界的图像分割方法。

目前已有一些颈动脉内中膜分割算法，能够实现较为准确的内中膜分割。然而，超声图像受固有的斑点噪声影响，信噪比较低，且解剖学细节不足。严重的斑点噪声甚至会掩盖内中膜层，影响了内中膜自动分割方法成功率。本发明针对这一问题，基于蚁群优化等算法，提出了一种具有更优鲁棒性的内中膜自动分割方法。

蚁群优化算法(antcolonyoptimization,aco)是一种受蚁群爬行特征启发而发展起来的一种仿生型的智能优化算法，其特点是通过仿生自适应个体的局部最优性共同确定问题的整体最优解。该算法具有自学习功能和解的强搜索能力，具 有并行化、鲁棒性强、正反馈的优点，现已被引入到图像分割领域当中。对于内中膜分割，可以采用蚁群算法，将ma和li边界的轮廓提取问题转换为最优化问题进行求解。

技术实现要素：

为了替代手动分割步骤，本发明提出了一种超声图像内中膜自动分割方法。本发明根据内中膜结构独特的双层平行界面结构，提出了一种双腿蚁群优化算法，并结合多尺度高斯核相乘方法、边缘检测算子和snake模型，实现了内中膜的自动分割。

本发明的技术方案如下：

一种超声图像内中膜自动分割方法，包括以下步骤：

1)使用多尺度高斯核相乘方法抑制内中膜以外的背景信息，增强内中膜边界，得到一个内中膜边界凸显的边缘图；

2)使用基于otsu阈值法和sobel算子的初始轮廓检测方法，根据步骤2)得到的边缘图，获取li和ma边界的部分轮廓线段；

3)在步骤2)获取的初始轮廓的基础上，设计一种双腿aco算法，将初始轮廓线段连接起来，得到一条完整的轮廓；

4)使用基于snake模型的轮廓优化算法来进一步优化步骤3)得到的轮廓，使最终轮廓更加平滑、连续、逼近真实边缘。

上述步骤1)使用多尺度高斯核相乘方法抑制内中膜以外的背景信息，突显出内中膜的li和ma边界，得到一个边缘图。边缘图定义为两个具有不同尺度的高斯密度核的滤波器与图像卷积结果的乘积，并且只保留值为正的部分：

其中，

是一个小尺度的二维高斯函数，σ1取值1～3。而gσ2(y)是一个大尺 度的一维高斯函数，σ2取值10～20。

步骤2)使用基于otsu阈值法和sobel算子的初始轮廓检测方法，主要分为以下几个步骤：

(1)基于步骤1)得到的边缘图，使用otsu阈值法获取一个二值化图像；

(2)利用水平sobel算子从二值化图像中提取出ma界面的上边缘和下边缘，及li界面的上边缘和下边缘；

(3)删去有冗余的边缘线和有缺损的边缘线；

(4)对于保留下的部分，取上边缘线和下边缘线的中点，作为每对边缘线合并后的最终单根边缘线。

步骤3)采用了一种双腿aco算法，将初始轮廓线段连接起成一条完整的轮廓。其特征为，在内中膜起始点放置了两只蚂蚁，令它们从血管的左端爬向右端。其中一只蚂蚁被固定放置在另一只蚂蚁的上方，然后二者同时爬动，因而可以本看作是一只蚂蚁的两条腿。另外，定义通过初始轮廓检测获取的轮廓线段为蚂蚁的强制路径，当蚂蚁爬行至这部分时，将被强制沿初始轮廓线段爬行，而线段之间的空隙将由aco算法来连接。其特征包括：

(1)转移概率方程

蚂蚁从手动选取的ps(xs,ys)点开始，从左向右爬行，每次固定向右移动1个横坐标。在t时刻，第k只蚂蚁的第l(l＝1,2)只腿从像素(x,i)移动至相邻像素(x+1,j)的过程依据下转移概率方程：

其中τ指信息素，η指像素的密度。α和β分别确定信息素和启发信息的相对影响力。其中，α取值为1～3，β取值为3～5。

代表第k只蚂蚁的第l只腿在下一步爬行中被允许的位置集合，它对应了像素(x,i)的右邻域。必须指出的是，邻域像素不能同时被两条腿共享。因此，如果两条腿的右邻域出现重合情况，重合的像素将会被从中排除。

(2)全局更新法则

当所有蚂蚁完成爬行之后，信息素按照如下公式进行更替：

τ(x,i)(t+1)＝ρτ(x,i)(t)+δτ(x,i)

其中，τ(x,i)(t)是第t次迭代时像素(x,i)处的信息素数量，τ(x,i)(t+1)是下一次迭代时像素(x,i)处的信息素数量。信息素的初始值τ0取值为0.1～1。ρ是衰减常数，用于仿真信息素的挥发，取值为0.5～1。δτ(x,i)是这一次迭代中释放的信息素数量，其计算公式为：

其中q是一个常数，m是蚂蚁的数量。q＝1，m取值为10～50。c(k)是第k只蚂蚁在搜索路径过程中的消耗函数，其定义如下：

其中in(x,y)(0≤in(xi,yi)≤1)是像素(x,i)归一化后的灰度值。d(k)是蚂蚁antk爬行路径的终点距离之前手动选择的终点pe(xe,ye)的距离。a是一个惩罚系数，设定了误差距离d(k)在消耗函数中的权重，取值为1～2。

上述步骤4)使用基于snake模型的轮廓优化算法来进一步轮廓，snake模型包含平滑能量项(smoothingenergy)、边缘能量项(boundaryenergy)和均一能量项(uniformenergy)，并通过最小化如下能量泛函来实现：

其中y1(x)和y2(x)分别代表li和ma交界面的轮廓，参数μ控制平滑能量项的权重，ν控制均一能量项的权重。均一能量项连接着着li和ma两条相互独立的轮廓，使它们保持均一的距离。μ取值为0.1～0.3，ν取值为1～2。

本发明具有以下优点：

本发明首先根据颈动脉超声图像特征，设计了一种多尺度高斯核相乘方法，对图像进行预处理，抑制内中膜以外的背景信息，增强内中膜边界。然后，在经过预处理的图像基础上，设计了一种双腿蚁群优化算法，并和otsu阈值法、sobel边缘检测算子、蚁群优化算法、snake模型等算法结合，实现了颈动脉超声图像内中膜的自动分割。基于临床图像的实验结果表明，本方法取得了准确的分割结果，其误差小于手动分割的观察者间误差；与此同时，本方法在临床数据测试取得了98.7％的成功率，具有良好的鲁棒性，并能够应对被斑点噪声严重污染的图像。

附图说明

图1是本发明分割颈动脉内中膜的流程图。

图2是初始轮廓检测过程中每一步的结果。

图3是蚁群优化算法中，蚂蚁上腿(upperleg)和下腿(lowerleg)的准许位置集合的典型示例。

图4是本发明分割结果示例图。

具体实施方式

以下通过实施例对本发明做进一步说明，以便更好地理解本发明的技术方案。步骤如下：

1.使用多尺度高斯核相乘方法得到一个内中膜边界凸显的边缘图f(x,y)。f(x,y)计算方法为，将两个具有不同尺度的高斯密度核的滤波器与图像卷积结果的相乘，并且只保留值为正的部分：

其中，

是一个小尺度的二维高斯函数，本实施例中σ1取值为1.5。而gσ2(y)是一个大尺度的一维高斯函数，本实施例中σ2取值为15。

2.使用基于otsu阈值法和sobel算子的初始轮廓检测方法，获取li和ma边界的部分初始轮廓线段。如图2所示，包含以下步骤：

1)基于由上一步得到的边缘图(a)，使用otsu阈值法获取一个二值化图像，如图(b)；

2)利用水平sobel算子从二值化图像(b)中提取出ma界面的上边缘和下边缘，及li界面的上边缘和下边缘；如图(c)所示，白色线表示上边缘，灰色线表示下边缘。

3)经过2)步骤后，理想结果应该含有2对(4条)轮廓线，分别对应li和ma界面的上下边缘。然而，使用sobel算子获取的边缘线通常都会含有多余的错误边缘线或缺口。因此，沿纵向看，多于或少于2对边缘线的部分将会被删去，结果如(d)所示。

4)最后，对于保留下的部分，取上边缘线和下边缘线的中点，作为每对边缘线合并后的最终单根边缘线，如(e)所示。

3.在步骤2获取的初始轮廓的基础上，设计一种双腿aco算法，将初始轮廓线段连接起来，得到一条完整的轮廓。其步骤包括：

1)初始化信息素矩阵。

2)将一直蚂蚁放置在手动选择的初始点ps(xs,ys)附近(上腿放置在点上方2个像素处，下腿放置在点下方2个像素处)。

3)蚂蚁从ps(xs,ys)开始向右爬动，每次固定向右移动1个横坐标，爬行的总步数等于xe-xs。

在t时刻，第k只蚂蚁的第l(l＝1,2)只腿从像素(x,i)移动至相邻像素(x+1,j)的过程依据下转移概率方程：

其中τ指信息素，η指像素的密度。α和β分别确定信息素和启发信息的相对影响力。本实施例中取值为α＝1,β＝4。

代表第k只蚂蚁的第l只腿在下一步爬行中被允许的位置集合，它对应了像素(x,i)的右邻域；但是，如果两条腿的右邻域出现重合情况，重合的像素将会被从中排除。这样可以保证蚂蚁的上腿始终在下腿的上方，因而可以保证最终的轮廓不会交叉或重合。图3展示了双腿蚂蚁爬行过程中的一个典型邻域蚂蚁此刻正从横坐标为x的位置向x+1位置爬行，图中浅色方格为蚂蚁上腿(upperleg)和下腿(lowerleg)的准许位置集合，由于深色方格同时存在于上腿和下腿的爬行准许位置集合中，因而被从二者的集合中排除。

4)计算蚂蚁路径上像素点的信息素。

5)重复步骤2-4，直到所有的蚂蚁完成了他们的爬行任务。

6)更新信息素矩阵，并计算信息素的挥发量。

信息素按照如下公式进行更替：

τ(x,i)(t+1)＝ρτ(x,i)(t)+δτ(x,i)

其中，τ(x,i)(t)是第t次迭代时像素(x,i)处的信息素数量，τ(x,i)(t+1)是下一次迭代时像素(x,i)处的信息素数量。本实施例中，信息素的初始值为τ0＝0.5。ρ是衰减常数，用于仿真信息素的会发，本实施例中选取ρ＝0.6。δτ(x,i)是这一次迭代中释放的信息素数量，其计算公式为：

其中q是一个常数，m是蚂蚁的数量。本实施例中选取q＝1，m＝20。c(k)是第k只蚂蚁在搜索路径过程中的消耗函数，其定义如下：

其中in(x,y)(0≤in(xi,yi)≤1)是像素(x,i)归一化后的灰度值。d(k)是蚂蚁antk爬行路径的终点距离之前手动选择的终点pe(xe,ye)的距离。a是一个惩罚系数，设定了误差距离d(k)在消耗函数中的权重，本实施例中取a＝1。

7)重复2-6步骤，进行一定数量的迭代。本实施例中迭代次数设置为10。

8)选取具有最小消耗方程的蚂蚁路径作为最优路径，即为li和ma界面的轮廓。

4.使用基于snake模型的轮廓优化算法来进一步优化轮廓，snake模型通过最小化如下能量泛函来实现：

其中y1(x)和y2(x)分别代表li和ma交界面的轮廓，参数μ控制平滑能量项的权重，ν控制均一能量项的权重。本实施例选取μ＝0.1，ν＝1.4。

图4显示了cgaco在三张不同示例图上的分割结果。图a是一个噪声严重的图像。由于噪声的污染，内中膜已经很不清晰。如图b所示，严重的斑点噪声并没有影响到内中膜的准确分割。图c是一个存在弯曲情况的颈动脉图像，图e是内中膜出现典型的明显增厚的颈动脉图像。图d和图f分别对应它们的分割结果。