CN106340021A - 血管提取方法 - Google Patents

Abstract

一种血管提取方法,包括：获取医学图像，与标准图像配准获取至少一个血管区域；确定所述血管区域的若干个种子点；采用分类器方法获取所述医学图像的血管增强图；在所述血管增强图中，基于水平集算法计算所述种子点之间的中心线；基于所述中心线生长血管，获取血管提取结果。本发明在采用分类器获取的血管增强图像上提取血管，有效改善结果的准确性和鲁棒性。

Description

血管提取方法

【技术领域】

本发明涉及医学图像处理领域，尤其是涉及一种血管提取方法。

【背景技术】

血管成像技术包括CT血管造影(CTA，CT angiography)，磁共振血管造影(Magnetic Resonance Angiography，MRA)，以及MR非造影成像等。血管成像技术是当前诊断血管疾病的重要方法之一,在临床上主要用于诊治多种血管疾病，如动脉瘤、血管狭窄、血管钙化等。

头颈部动脉主要包括颈总动脉(CCA)，颈内动脉(ICA)，颈外动脉(ECA)，椎动脉(VA)，基底动脉(BA)等。颈总动脉分叉为颈内动脉和颈外动脉，如图1a和1b分别为包含椎动脉和颈内动脉的医学图像所示，左右椎动脉在一节一节的椎骨中穿行，最后合并成基底动脉穿过过枕骨，为大脑后部供血。颈内动脉穿过颅骨，给大脑的前部和中部供血；颈外动脉为牙齿和面部神经供血。临床上，分割出的血管可以被用来定量的诊断各种血管疾病，如颈总动脉和颈内动脉易发动脉粥样硬化，椎动脉易狭窄等，这些病变也是大脑中风的主要原因。但是，由于血管与周围骨头紧邻，空间位置错综复杂，同时增强后的血管强度值和骨的强度值存在部分重叠，诸如区域生长等常用算法很容易混淆血管和骨头，造成血管分割失败。

为解决上述技术问题，现有技术中采用的血管提取方法有：

1.基于模型的方法实现穿骨血管的提取。首先将医学图像与已有的血管模型配准，划分不同的血管区域及得到不同血管的起始点，然后通过水平集算法得到血管的中心线。再以中心线为起始，通过单纯网格向外扩张，提取血管。但由于网格计算非常复杂，导致分割速度很慢。

2.基于Graph Cut的算法实现了颈内动脉的提取。该方法需要手动标记血管起始点及血管大致半径，然后计算起始点之间的最短路径。以该路径为中心，标记血管感兴趣区域。最后通过Graph Cut算法，优化感兴趣区域的边界，实现血管提取。该算法计算复杂，血管提取速度慢，且很难实现全自动血管提取。

3.基于Ray-Casting的动态追踪实现了穿颅血管的提取。该算法需要手动确定起始点，之后通过Ray-Casting向四周放射线，获得血管的边界，再通过椭圆拟合实现血管提取。该方法鲁棒性差，不同数据的血管提取需要不同的参数，同时该方法无法提取旋转跨度很大的椎动脉。

因此，有必要改进现有的血管提取方法。

【发明内容】

本发明所要解决的技术问题是提供一种血管提取方法，该方法能够改善现有血管提取方法，特别是改善提取穿颅或紧贴骨骼组织的血管的准确性和鲁棒性。

本发明为解决上述技术问题而采用的技术方案是：提供一种血管提取方法,包括以下步骤：

获取医学图像，与标准图像配准获取至少一个血管区域；

确定所述血管区域的种子点；

采用分类器方法获取所述医学图像的血管增强图；

在所述血管增强图中，基于水平集算法计算所述种子点之间的中心线；

基于所述中心线生长血管，获取血管提取结果。

可选的，所述分类器方法，包括：

采集训练样本；

读取所述训练样本，训练得到弱分类器；

组合所述弱分类器，获取强分类器；

基于所述强分类器，获取所述医学图像的血管增强图。

可选的，所述采集训练样本包括如下步骤：

基于所述标准图像采集若干个训练样本；以及

膨胀所述标准图像，对膨胀前后的标准图像进行异或操作，基于该异或操作结果采集所述训练样本。

可选的，所述训练样本包括若干个正样本和负样本。

可选的，所述获取强分类器包括：

计算所述训练样本的Haar特征，利用AdaBoost算法从所述Haar特征选取有效Haar特征，训练得到若干个弱分类器；

筛选弱分类器并分配对应权重，基于筛选后的分类器及其权重组合构成所述强分类器；

其中，每个有效Haar特征组成单个弱分类器。

可选的，获取所述血管增强图还包括：基于所述强分类器，计算所述医学图像上的血管区域中任一个点的血管增强值，获取所述血管增强图。

可选的，若干个所述强分类器级联构成AdaBoost级联强分类器。

可选的，获取所述中心线包括：

基于血管特征，计算所述血管增强图中各像素点或体素点的相似度值；

以所述相似度作为解，采用快速行进算法构造U场；

基于能量最低原则在所述U场中连接所述种子点，获取所述中心线。

可选的，所述血管特征包括图像灰度指、梯度值、增强值、形状、灰度值方差或前述组合。

可选的，所述中心线采用最速下降法连接所述种子点获取。

本发明对比现有技术有如下的有益效果：本发明采用分类器方法获取所述医学图像的血管增强图，特别是针对血管区域计算血管增强值，以突出血管区域，减少骨头部分对血管提取的影响，提高提取结果的精度和速度；在所述血管增强图中，基于水平集算法计算所述种子点之间的中心线，根据所述血管增强图像的相似度，采用快速行进算法构造U场；基于能量最低原则连接所述种子点，获取所述中心线。本发明提供的血管提取方法相比现有技术可以提高血管提取的准确性和鲁棒性。

【附图说明】

图1a和1b分别为包含椎动脉和颈内动脉的医学图像示意图；

图2为本发明一实施例中的血管提取流程示意图；

图3为本发明一实施例中医学图像中的血管区域示意图；

图4为本发明一实施例中获取血管增强图的流程示意图；

图5a和5b分别为本发明一实施例中椎动脉和颈内动脉的血管增强图的结果示意图；

图6为本发明一实施例中获取中心线的流程示意图；

图7为本发明一实施例中获取U场的结果示意图；

图8为本发明一实施例中提取头颈部血管的结果示意图。

【具体实施方式】

在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施的限制。

其次，本发明利用示意图进行详细描述，在详述本发明实施例时，为便于说明，所述示意图只是实施例，其在此不应限制本发明保护的范围。使本发明的上述目的、特征和优点能够更为明显易懂，下面结合附图和实施例对本发明作进一步的描述。

请参阅图2，本发明的实施例的一种血管提取方法，包括以下步骤：

S10：获取医学图像，与标准图像配准获取至少一个血管区域；本实施例中，所述医学图像可以通过CT血管造影(CTA)，磁共振血管造影(MRA)等血管造影成像技术，以及MR非造影成像获取。所述医学图像可以是三维序列图像，其中包含若干个二维切片图像。所述医学图像包含血管、以及颅骨、枕骨等骨头组织区域。所述血管区域包括静脉窦、脑中动脉、脑前动脉和脑后动脉等血管区域。由于血管与周围骨头紧邻，空间位置错综复杂，同时增强后的血管强度值和骨头的强度值存在部分重叠，从而影响医生对血管的检测及诊断。

本实施例中，所述医学图像与标记有血管图谱的标准图像进行配准，根据映射关系将所述标准图像中的若干个血管区域匹配至所述医学图像，获取至少一个血管区域。所述标准图像是标记正常人体血管系统或者参考病例中对血管区域的信息叠加，例如各血管提取结果的金标准。配准的方法可以基于一定映射关系(例如基于仿射配准关系矩阵和非刚体配准的形变场)，将每个参考病例中血管区域(例如主血管)映射到一个坐标空间，例如，标准图像的坐标空间，并计算出包含有一个主血管区域的概率分布。进一步地，在一些实施例中，可以动态选择适应的配准方法或者空间变换模型。图3所示为配准后的医学图像示意图，分别表示头颈部血管的正面。其中区域A是左颈内动脉，区域B是右颈内动脉，区域C是基底动脉，区域D是右椎动脉，和区域E是左椎动脉。

S20：确定所述血管区域的种子点。所述种子点可以通过手动、半自动或者自动提取于所述血管区域中。本实施例中，可以在前述配准操作获取的血管区域图像，根据血管的特征，例如图像灰度值、梯度值、增强值、形状、灰度值方差、空间形态和位置关系等特征确定若干个种子点。例如，在血管区域内选取初始点，基于水平集算法等追踪方法，获取血管区域中至少两个端点即为种子点。又例如，左、右颈内动脉在头颈部血管图像中呈左右对称状态，可以根据空间构型确定种子点。

S30：采用分类器方法获取所述医学图像的血管增强图。所述分类器方法可以结合图4作详细描述：

S31：采集训练样本。所述训练样本包括正样本和负样本，所述正样本为血管样本点，即血管区域内的像素点或体素点；优选的，所述正样本为血管中心区域上的点，即中心线附近的点，所述中心线可以为位于血管内部沿着血管走向的一条线，也可以指位于血管中心部位或靠近血管中心部位像素点的集合，还可以指到血管边界距离相等或基本相等的像素点的连线。所述负样本为非血管样本点，即背景区域中的像素点或体素点。所述训练样本来源可以在标准图像中选取血管区域内的点为正样本，血管周围区域的点为负样本。还可以对标准图像中的血管区域作膨胀操作，与膨胀结果中选取正、负样本，此处所述膨胀操作是以一定宽度使边界向外部扩张的过程，即以一定宽度将与血管接触的背景点融合之血管区域中。进一步地，对膨胀前后的血管区域作异或操作，即排除膨胀前、后图像中相同的点，提取得到不同的点，即获取血管周围区域，在此区域中提取负样本。所述正样本和负样本的数量可以是任意的。正负样本的数量不一定要保持一致。例如，负样本的数量可以多于正样本。本实施例中选取正样本为m个，负样本为n个。

S32：读取所述训练样本，训练得到弱分类器；训练的目的在于使得该分类器对所有训练样本的分类误差达到最小，即通过一轮分类器的训练找到分类错误率最低，分类效果最好的分类器。本实施例中，计算通过所述正样本和负样本的Haar特征值，基于AdaBoost算法从Haar特征选取有效Haar特征，每个有效Haar特征组成单个弱分类器。所述Haar特征是常用于识别物体的特征描述算子，适用于跟踪结构稳定的物体。目前常用的Haar特征基本可以分为四类：边缘特征、线性特征、点特征和对角线特征，他们反映了图像的灰度变化情况。

本实施例中，通过训练样本获取所述单个弱分类器步骤可以为：对于一幅图像x(例如256×256的二维图像)或者是一个子窗口(完整包括一个血管区域的矩形窗口，例如横断面图像上30×30，此处256和30均是像素点个数)，共计有k个Haar特征，选取其中一个Haar特征t，计算图像x中所有训练样本关于此Haar特征t的值f_t(x)，并选取一个和此特征t相关的阈值θ_t和偏位符p_t，其中所述偏位符p_t用于控制不等式的方向，取值为1或-1，当f_t(x)>θ_t时令p_t为1，当f_t(x)<θ_t时令p_t为-1。弱分类器的计算公式如公式(1)所示：

通过计算每一个训练样本关于该Haar特征t的值f_t(x)，将特征值按照升序进行排列，记为f_t(1)，f_t(2)……f_t(m+n)，从中选取最佳阈值θ_t作为分类的分水岭，用以分类输入的训练样本为正样本或负样本，使得利用该特征t对所有的训练样本的分类误差达到最小。所述阈值θ_t的选取与最小错误分类率有关，最小错误分类率e_t的计算方式如公式(2)所示：

e_t＝min(S_t ⁺+(T_t ^--S_t ^-),S_t ^-+(T_t ⁺-S_t ⁺)) (2)

其中，min表示取最小值函数，T_t ⁺表示所有正样本的权重和，T_t ^-表示所有负样本的权重和，S_t ⁺为在θ_t阈值下所有特征值小于θ_t的正样本权重和，S_t ^-为在θ_t阈值下所有特征值小于θ_t的负样本权重和。初始状态下，所有的样本都被赋予相等的权重值。根据最小错误分类率e_t，确定最佳阈值θ_t,该阈值θ_t对应的特征t即为有效Haar特征，为所需要的弱分类器h(x)。

需要说明的是，伴随着学习过程的进行，根据训练得到的弱分类器，可以判断正、负样本是否被准确分类，如果该像素点被准确分类，则在下一次训练中该像素点的权重将会被降低，否则，若该像素点被分类错误，则该像素点的权重被提高，从而得到一个新的样本分布。即，在每轮的迭代中需对于所有正样本和负样本的权重根据上一轮的训练结果进行归一化，其权重归一化公式为：

其中，T表示迭代次数，ω_T(j)表示第T次迭代中第j个样本的权重。遍历所有特征，利用前述步骤选取的包含一个Haar特征的最佳弱分类器h(x，f_t，p_t，θ_t)，并根据公式(3)计算该弱分类器对m+n个样本的加权错误分辨率，对应每个有效Haar特征的弱分类器群，找出错误率ε_t最小的分类器h_j，选择器对应的特征t。

随后更新下一轮训练样本的权重其中如果样本被正确分类则e_j＝0，否则，如样本被错误分类，则e_j＝1。按照上述AdaBoost迭代算法，通过R轮迭代，将选取得到R个Haar特征，本实施例中所述R选取为50。

S33：组合所述弱分类器，获取强分类器；经过R次训练完成后获取的R个弱分类器及其权重组合在一起，构成强分类器。进一步地，若干个所述强分类器级联构成AdaBoost级联分类器，强分类器H(x)表达式如公式(4)。其中，若H(x)＝1代表该样本分类为正样本，即该点为血管上的点，若H(x)＝0代表该点分类为负样本，该点为非血管区域上的点。

其中为有效Haar特征t对应的弱分类器所赋予的权重，表示该弱分类器的性能。权重越小代表该弱分类器在最终的强分类器中起着较小的作用，权重越大代表该弱分类器在最终的强分类器中起着较大的作用。虽然对于单个弱分类器来说，其分类的准确率不高，但是通过赋予各个弱分类器不同的权重，可以有效的提高分类器的准确率，即。

S34：基于所述强分类器，获取所述医学图像的血管增强图。本实施例中，为提高提取血管的速度，优选计算配准后获取的血管区域的血管增强值，由前述步骤获取的强分类器或者级联强分类器计算血管区域中各像素点或者体素点的血管增强值(即Vessel值)，得到如图5a和5b所示的血管增强图，图5a为针对椎动脉所在区域获得的血管增强图，图5b所示为针对颈内动脉所在区域获得的血管增强图，图5a和图5b中将非血管区域设置为背景区域，针对血管区域计算获取所述血管增强图，有效减少骨骼组织对血管提取的影响，凸显待提取的血管区域，提高血管提取的准确性。本实施例中，血管增强图可以参照公式(5)所示获取：

$Vessel\left(x\right)={\mathrm{\&Sigma;}}_{T=1}^R{\mathrm{\&alpha;}}_t{h}_t\left(x\right)---\left(5\right)$

其中h_t(x)为所述血管区域中某一像素点或者体素点x对应有效Haar特征t的弱分类器，α_t代表表示各个弱分类器对应的权重，Vessel(x)表示点x对应的血管增强值。当某一点的血管增强值大于所设定的阈值时，则可以判断该点为血管内的点可能性较高；当某一点的血管增强值低于阈值时，则判断该点为血管内的点可能性较低。根据判断结果，可以得到血管增强图。所述血管增强值越大，表示该点更像血管点，在所述血管增强图中显示为灰度值大小的不同。

S40：在所述血管增强图中，基于水平集算法计算所述种子点之间的中心线。所述中心线可以指位于血管内部的沿着血管走向的一条线。在一些实施例中，血管中心线可以指位于血管中心部位或靠近血管中心部位像素点的集合。在一些实施例中，血管中心线可以指到血管边界距离相等的像素点的连线。

本实施例中基于水平集方法提取中心线，可以基于前述步骤所确定的种子点，设置起止两个种子点的水平集函数场的初始曲线。这两个水平集函数场会随着时间不断演化。在变化较快的路径上，两个水平集函数场有一个相交点，这个交点即为血管上的点。所述交点出发，利用根据最短距离原则，确定两个起止种子点间的血管中心线。

优选的，可以基于快速行进算法(Fast Marching Method)确定血管中心线。快速行进算法(Fast Marching Method，FMM)通过检查种子点与邻域点内相邻点之间所有可能路径来选择达到能量最低的点，即根据相似度P值选择相邻点之间到达能量最低的点。获取所述中心线的方法可结合附图6作具体描述：

S41：基于血管特征，计算所述血管增强图中各像素点或体素点的相似度值。在前述步骤获取的血管增强图中，以所述种子点为起始点，根据血管特征，例如，根据所述血管增强图中血管区域的图像灰度指、梯度值、增强值、形状等中的一种或多种组合分析种子点n领域内各相邻点与血管的相似度P值，所述领域可是26邻域，也是可以6邻域；所述相似度P值可以根据公式(6)获取：

$P_{i,j,k}=\frac{\left(Grad\right(x)+|Std\left(x\right)-Std\left(x_0\right)\left|\right)}{Vessel\left(x\right)}---\left(6\right)$

其中，Grad(x)为当前点的梯度值，Std(x)为当前点与周围点的方差，Std(x₀)为血管的种子点对应的方差，以及Vessel(x)为根据公式(5)获取的所述血管增强图上对应该点的血管增强值。由于血管内部为均质属性，属于血管内的点之间的梯度值较小，同理血管内的点之间的方差也较小，与血管的种子点对应的方差的差值绝对值|Std(x)-Std(x₀)|也较小，所以对于公式(6)中P值越小的像素点或体素点，则越像血管，即从种子点为起点到达该像素点的能量最低，从图像上显示为围绕血管中心区域的种子点连线。需要说明的是，公式(6)经为本实施例中的一个示例，其中梯度值Grad(x)和方差Std(x)可以根据血管的特征进行替换，可以取得等同效果。

S42：以所述相似度作为解，采用快速行进算法构造U场，所述U场中的U值可以根据公式(7)获取：

$\begin{array}{c}{\left(max\right\{u-U_{i-1,j,k},u-U_{i+1,j,k},0\left\}\right)}^2+{\left(\max \right\{u-U_{i,j-1,k},u-U_{i,j+1,k},0\left\}\right)}^2\\ +{\left(\max \right\{u-U_{i,j,k-1},u-U_{i,j,k+1},0\left\}\right)}^2=P_{i,j,k}^2\end{array}---\left(7\right)$

其中，u为粘性解常数，U值代表所述两个种子点通过点(i，j，k)路径所需的最小能量值，即两个种子点通过点(i，j，k)路径所需能量最低。根据公式5获取的点(i，j，k)的相似度P值，代入公式(7)，得到该点的U值，获取所述血管区域中任意点的U值，得到如图7所示的U场效果图。

S43：基于能量最低原则在所述U场中连接所述种子点，获取所述中心线。本实施例中，优选采用最速下降法在根据所述U场，连接所述种子点，获取所述中心线。

S44：基于所述中心线生长血管，获取血管提取结果。本实施例中，可以基于前述步骤提取的中心线作为初始分割结果，基于血管特征构建血管模型方法，拟合血管结构，不断优化血管边界获取血管分割结果。也可以所述中心线为起始种子点，基于区域生长获取血管分割结果。如图8所示为本实施例中头颈部血管提取结果的示意图。

综上所述，本发明采用分类器方法获取所述医学图像的血管增强图，特别是针对血管区域获取血管增强图，用以突出显示血管区域，减小骨头部分对血管提取的影响，提高提取结果的精度和速度；在所述血管增强图中，基于水平集算法计算所述种子点之间的中心线，根据所述血管增强图像的相似度，以所述相似度作为解，采用快速行进算法构造U场；根据所述U场，基于能量最低原则连接所述种子点，获取所述中心线，有效提高血管中心线提取的精度和鲁棒性。另外，本发明是一种自动化的血管提取方法，工作效率较高。

需要说明的是，通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到本申请的部分或全部可借助软件并结合必需的通用硬件平台来实现。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可包括其上存储有机器可执行指令的一个或多个机器可读介质，这些指令在由诸如计算机、计算机网络或其他电子设备等一个或多个机器执行时可使得该一个或多个机器根据本发明的实施例来执行操作。机器可读介质可包括，但不限于，软盘、光盘、CD-ROM(紧致盘-只读存储器)、磁光盘、ROM(只读存储器)、RAM(随机存取存储器)、EPROM(可擦除可编程只读存储器)、EEPROM(电可擦除可编程只读存储器)、磁卡或光卡、闪存、或适于存储机器可执行指令的其他类型的介质/机器可读介质。

本申请可用于众多通用或专用的计算系统环境或配置中。例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备、多处理器系统、基于微处理器的系统、置顶盒、可编程的消费电子设备、网络PC、小型计算机、大型计算机、包括以上任何系统或设备的分布式计算环境等。

本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。

虽然本发明已以较佳实施例揭示如上，然其并非用以限定本发明，任何本领域技术人员，在不脱离本发明的精神和范围内，当可作些许的修改和完善，因此本发明的保护范围当以权利要求书所界定的为准。

Claims (10)

1.一种血管提取方法,其特征在于，包括：

获取医学图像，与标准图像配准获取至少一个血管区域；

确定所述血管区域的若干个种子点；

采用分类器方法获取所述医学图像的血管增强图；

在所述血管增强图中，基于水平集算法计算所述种子点之间的中心线；

基于所述中心线生长血管，获取血管提取结果。

2.根据权利要求1所述的方法，其特征在于，所述分类器方法，包括：

采集训练样本；

读取所述训练样本，训练得到弱分类器；

组合所述弱分类器，获取强分类器；

基于所述强分类器，获取所述医学图像的血管增强图。

3.根据权利要求1所述的方法，其特征在于，所述采集训练样本包括：

基于所述标准图像采集若干个训练样本；以及

膨胀所述标准图像，对膨胀前后的标准图像进行异或操作，基于该异或操作结果采集所述训练样本。

4.根据权利要求2或3任意一项所述的方法，其特征在于，所述训练样本包括若干个正样本和负样本。

5.根据权利要求2所述的方法，其特征在于，所述获取强分类器包括：

计算所述训练样本的Haar特征，利用AdaBoost算法从所述Haar特征中选取有效Haar特征，训练得到若干个弱分类器；

筛选弱分类器并分配对应权重，基于筛选后的分类器及其权重组合构成所述强分类器；

其中，每个有效Haar特征组成单个弱分类器。

6.根据权利要求2所述的方法，其特征在于，获取所述血管增强图还包括：

基于所述强分类器，计算所述医学图像上的血管区域中任一个点的血管增强值，获取所述血管增强图。

7.根据权利要求2所述的方法，其特征在于，若干个所述强分类器级联构成级联强分类器。

8.根据权利要求1所述的方法，其特征在于，获取所述中心线包括：

基于血管特征，计算所述血管增强图中各像素点或体素点的相似度；

以所述相似度作为解，采用快速行进算法构造U场；

基于能量最低原则在所述U场中连接所述种子点，获取所述中心线。

9.根据权利要求8所述的方法，其特征在于，所述血管特征包括图像灰度指、梯度值、增强值、形状、灰度值方差、位置关系或前述组合。

10.根据权利要求8所述的方法，其特征在于，所述中心线采用最速下降法连接所述种子点获取。