一种基于图像梯度矢量流场的血管脊点提取方法及装置

摘要

本申请提供一种基于图像梯度矢量流场的血管脊点提取方法及装置，其中，方法包括：对血管造影图像中的血管目标进行增强；利用梯度矢量流场模型求取增强后图像的梯度矢量流场；根据增强后图像的梯度矢量流场，分别计算增强后图像中的一像素点处矢量与该像素点沿一坐标轴方向的两个邻像素点处矢量夹角的余弦值，若所述余弦值满足阈值条件，则该像素点为该坐标轴方向的血管脊点。本申请能够有效压制背景噪声影响、提升血管造影图像中血管脊点的提取效率、增加血管目标脊点的探测数量及提高血管脊点探测精度，为后续血管中心线的提取和血管建模提供基础。

说明书

技术领域

本申请属于医学图像处理领域，特别涉及一种基于图像梯度矢量流场的血管脊点提取方法及装置。

背景技术

随着CT血管造影、磁共振血管造影(MRA)等医学成像技术的快速发展，利用图像后处理技术获取血管脊点和中心线信息对于描述血管的形态结构、实现血管目标重建等具有非常关键的意义，这一问题的有效解决，能够提高血管疾病诊断、手术计划和手术导航的自动化水平和实施精度，从而提高诊断的效率和手术成功率。

现有技术中，在管状目标(血管)为高信号、背景为低信号、管状目标截面的轮廓呈高斯状分布这一标准条件下，垂直于管状目标径向的局部灰度极大值点被认为管状目标脊点。传统的血管脊点的提取方法是根据局部灰度极大值的定义，通过图像的一阶微分和二阶微分进行判断得到局部脊点。具体的，首先，通过一阶微分判断极值点，对图像中的像素点来说，梯度值为0的点为局部极值点的充分条件，即满足其次，根据二阶微分特性：即沿垂直于管状目标径向方向的Hessian矩阵的特征矢量v_i对应的特征值λ_i为负的点，作为存在血管脊点的必要条件。

上述血管脊点的提取方法对图像背景中存在的噪声非常敏感，通常会错误地把局部噪声亮点识别为极值点，从而很大程度增加了图像中非目标脊点数量。

发明内容

本申请提供一种基于图像梯度矢量流场的血管脊点提取方法及装置，用于解决现有技术中，基于局部极值点定义的脊点提取方法存在易受高强度噪声点影响、提取到的脊点精度不高的问题。

为了解决上述技术问题，本申请的一技术方案为提供一种基于图像梯度矢量流场的血管脊点取值方法，包括：

对血管造影图像中的血管目标进行增强；

利用梯度矢量流场模型求取增强后图像的梯度矢量流场；

根据增强后图像的梯度矢量流场，分别计算增强后图像中的一像素点处矢量与该像素点沿一坐标轴方向的两个邻像素点处矢量夹角的余弦值，若所述余弦值满足阈值条件，则该像素点为该坐标轴方向的血管脊点。

本申请另一技术方案为提供一种基于图像梯度矢量流场的血管脊点提取装置，该装置包括：

增强模块，用于对血管造影图像中的血管目标进行增强；

梯度矢量流场计算模块，用于利用梯度矢量流场模型求取增强后图像的梯度矢量流场；

脊点探测模块，根据增强后图像的梯度矢量流场，分别计算增强后图像中的一像素点处矢量与该像素点沿一坐标轴方向的两个邻像素点处矢量夹角的余弦值，若所述余弦值满足阈值条件，则该像素点为该坐标轴方向的血管脊点。

本申请能够有效压制背景噪声影响、提升血管造影图像中血管脊点的提取效率、增加血管目标脊点的探测数量及提高血管脊点探测精度，为后续血管中心线的提取和血管建模提供基础。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本申请实施例的基于图像梯度矢量流场的血管脊点提取方法的流程图；

图2a为本申请实施例的孤立脊点示意图；

图2b为本申请实施例的孤立脊点团示意图；

图3a为本申请实施例的磁共振三维血管造影图像；

图3b为图3a经过多尺度血管增强函数处理后的图像数据的最大密度投影图像；

图4a为本申请实施例的血管增强后的图像；

图4b为对图4a所示图像利用梯度矢量流场模型求取图像的梯度矢量流场的示意图；

图5a为本申请实施例的血管造影图像的三维仿真图；

图5b为对图5a所示图像进行脊点提取得到的结果图；

图6为本申请实施例的基于图像梯度矢量流场的血管脊点提取装置的结构图。

具体实施方式

为了使本申请的技术特点及效果更加明显，下面结合附图对本申请的技术方案做进一步说明，本申请也可有其他不同的具体实例来加以说明或实施，任何本领域技术人员在权利要求范围内做的等同变换均属于本申请的保护范畴。

如图1所示，图1为本申请实施例的基于图像梯度矢量流场的血管脊点提取方法的流程图。本实施例能够有效压制背景噪声影响、提升血管造影图像中血管脊点的提取效率、增加血管目标脊点的探测数量及提高血管脊点探测精度，为后续血管中心线的提取和血管建模提供基础。

具体的，该方法包括：

步骤101：对血管造影图像中的血管目标进行增强。

步骤102：利用梯度矢量流场模型求取增强后图像的梯度矢量流场。

梯度矢量流场模型是一种全局优化的矢量场，图像的梯度矢量流场中各点(像素点)的矢量指向血管目标的脊点，即血管目标径向的局部灰度极大值点，根据该特征，通过如下步骤103能够精确、稳定、快速的提取到血管造影图像中的脊点。

步骤103：根据增强后图像的梯度矢量流场，分别计算增强后图像中的一像素点处矢量与该像素点沿一坐标轴方向的两个邻像素点处矢量夹角的余弦值，若所述余弦值满足阈值条件，则该像素点为该坐标轴方向的血管脊点。

本申请所述的血管造影图像可以为二维图像，也可以为三维图像，对于二维血管造影图像，其坐标轴方向包括±x、±y方向，对于三维血管造影图像，其坐标轴方向包括±x、±y、±z方向。实施时，重复步骤103，直至遍历增强后图像的所有像素点及所有坐标轴方向，得到全部血管脊点为止。

一实施例中，得到全部血管脊点之后还包括：从得到的全部血管脊点中剔除孤立点。详细的说，从得到的全部血管脊点中剔除孤立点的过程包括：

判断某一脊点A在以其为圆心，半径为d和d+d₀的两个同心圆所组成的圆环范围内是否存在其它脊点，如果不存在，则将该脊点A为圆心，半径为d形成的圆范围内的脊点删除，其中，d、d₀为距离常数。

如图2a、图2b所示，图2a中，圆环中不存在其他脊点，则将孤立脊点A删除；图2b中，圆环中不存在其他脊点，则将脊点A、B、C、D组成的孤立脊点团删除。实施时，可依据视觉上对孤立点的主观把握，结合噪声伪脊点分布、真实脊点分布及脊线特性三者的关联性确定距离常数d和d₀的值，一般取d＝5～8个像素距离，d₀＝3～4个像素距离。

一实施例中，对于三维血管造影图像，上述步骤101可通过如下第一多尺度血管增强函数对三维血管造影图像中的血管目标进行增强：

v₀(s)为第一多尺度血管增强函数；R_A、R_B和S为三个测度函数；R_A用来区分片状和线状结构；R_B用来区分点状结构和线状结构；S用于区分背景像素；α、β和c为阈值，用于控制血管增强算法对R_A、R_B和S的敏感性；λ₁、λ₂和λ₃为Hessian矩阵H的三个特征值，且满足|λ₁|≤|λ₂|≤|λ₃|，D为图像的维度；s为某个像素点。

需要说明的是，Hessian矩阵H为一个三阶偏导数构成的方阵，可采用现有方法计算得到，阈值α和β通常情况下取0.5，阈值c的取值依赖于图像的灰度范围，通常取最大值Hessian矩阵范数的一半。

对于二维血管造影图像，上述步骤101通过如下第二多尺度血管增强函数对二维血管造影图像中的血管目标进行增强：

其中，v₀'(s)为第二多尺度血管增强函数，R_B'和S'为测度函数，R_B'用来区分点状结构和线状结构，S'用于区分背景像素，β和c为阈值，λ₁和λ₂为Hessian矩阵H的二个特征值，且满足|λ₁|≤|λ₂|，D为图像的维度。

血管造影图像经过公式(1)或(2)的多尺度血管增强函数处理后，使得血管目标灰度值得到增强，背景噪声得到抑制。

多尺度血管增强函数通过对血管造影图像进行多尺度高斯滤波及Hessian矩阵特征值计算得到，具体计算过程参考Frangi,A.F.,Niessen,W.J.,Vincken,K.L.,Viergever,M.A.,1998.Multiscale vessel enhancement filtering.In:Proceedings of the International Conference on Medical Image Computing Computer Assisted Intervention.Lect.Notes Comp.Sci.,1496,pp.130–137，本申请不再进行详细说明。

一具体实施例中，如图3a、3b所示，图3a为本申请实施例的磁共振三维血管造影图像，图3b为图3a经过血管增强函数处理后的图像数据的最大密度投影图像。由图3a及图3b可以看出，经过多尺度血管增强函数处理后的图像背景噪声得到了抑制，同时血管中心线上的点具有垂直于血管方向的最大灰度值。

详细的说，梯度矢量流场定义为V(x,y)＝(u(x,y),v(x,y))，其中，u(x,y)，v(x,y)分别为矢量V的两个分量，可通过最小化能量泛函得到，具体计算公式如下：

其中，ε为能量泛函；(x,y)为增强后图像像素点的坐标；u_x、u_y、v_x、v_y为分量u、v分别对于x、y的一阶偏导；f(x,y)为增强后图像的边缘；为f(x,y)的梯度，μ为控制参数。

μ可根据图像的质量(如噪声情况)来设定，一般μ值越小，梯度矢量流场动态范围越小，能探测较细的血管，但越容易受到噪声影响。

利用变分法可知，梯度矢量流场满足如下欧拉方程组：

其中，f_x(x,y)和f_y(x,y)为增强后图像的血管边缘函数f(x,y)在x和y方向的值，u(x,y)、v(x,y)为梯度矢量流场V(x,y)＝(u(x,y),v(x,y))的两个分量，(x,y)为增强后图像像素点的坐标，为laplacian算子，μ为控制参数。

对于二维血管造影图像，上述步骤102中，利用梯度矢量流场模型求取增强后图像的梯度矢量流场进一步包括：

求解上述公式(4)，得到增强后二维血管造影图像的梯度矢量流场两个分量的迭代公式为：

其中，n为迭代次数。

同理，可以得到增强后三维血管造影图像的梯度矢量流场三个分量的迭代公式为：

求取公式(5)的收敛值，得到增强后二维血管造影图像的梯度矢量流场：

V(x,y)＝(u(x,y),v(x,y))。

求取公式(6)的收敛值，得到增强后三维血管造影图像的梯度矢量流场：

V(x,y,z)＝(u(x,y,z),v(x,y,z),w(x,y,z))。

一具体实施例中，如图4a、4b所示，图4a为本申请实施例的血管增强后的图像，图4b为对图4a所示图像利用梯度矢量流场模型求取图像的梯度矢量流场的示意图，由图4b可以得出，矢量的箭头指向血管目标中心线，血管目标中心线处梯度最强，而远离中心线处的梯度场为0。

下面以两个具体实施例说明上述步骤103进行脊点探测的过程：

(1)二维血管造影图像脊点探测

对于二维血管造影图像的x轴方向，判断某一像素点是否满足如下公式，如果满足，则该像素点为x轴方向的脊点：

其中，(i,j)为某一像素点的坐标，(i-1,j)和(i+1,j)为(i,j)在x轴方向的两个邻像素点，V_x^(i,j)、V_x^(i-1,j)、V_x^(i+1,j)分别为(i,j)、(i-1,j)、(i+1,j)像素点处矢量在x轴方向的分量，V_y^(i,j)、V_y^(i-1,j)、V_y^(i+1,j)分别为(i,j)、(i-1,j)、(i+1,j)像素点处矢量在y轴方向的分量，T₀为阈值。

对于二维血管造影图像的y轴方向，判断某一像素点是否满足如下公式，如果满足，则该像素点为y轴方向的脊点：

其中，(i,j)为某一像素点的坐标，(i,j-1)和(i,j+1)为(i,j)在y轴方向的两个邻像素点，V_x^(i,j)、V_x^(i,j-1)、V_x^(i,j+1)分别为(i,j)、(i,j-1)、(i,j+1)像素点处矢量在x轴方向的分量，V_y^(i,j)、V_y^(i-1,j)、V_y^(i,j+1)分别为(i,j)、(i,j-1)、(i,j+1)像素点处矢量在y轴方向的分量，T₀为阈值。

需要说明的是，每一像素点处的矢量为归一化处理后的矢量，该矢量的模为1。阈值可根据提取精度进行选取，通常T₀选取[0,0.5]范围内的值。

(2)三维血管造影图像脊点探测

对于三维血管造影图像的x轴方向，判断某一像素点是否满足如下公式，如果满足，则该像素点为x轴方向的脊点：

其中，(i,j,k)为某一像素点的坐标；(i-1,j,k)和(i+1,j,k)为(i,j,k)在x轴方向的两个邻像素点；V_x^(i,j,^k)、V_x^(i-1,j,k)、V_x^(i+1,j,k)分别为(i,j,k)、(i-1,j,k)、(i+1,j,k)像素点处矢量在x轴方向的分量；V_y^(i,j,k)、V_y^(i-1,j,k)、V_y^(i+1,j,k)分别为(i,j,k)、(i-1,j,k)、(i+1,j,k)像素点处矢量在y轴方向的分量；V_z^(i,j,k)、V_z^(i-1,j,k)、V_z^(i+1,j,k)分别为(i,j,k)、(i-1,j,k)、(i+1,j,k)像素点处矢量在z轴方向的分量；T₀为阈值。

对于三维血管造影图像的y轴方向，判断某一像素点是否满足如下公式，如果满足，则该像素点为y轴方向的脊点：

其中，(i,j,k)为某一像素点的坐标；(i,j-1,k)和(i,j+1,k)为(i,j,k)在y轴方向的两个邻像素点；V_x^(i,j,k)、V_x^(i,j-1,k)、V_x^(i,j+1,k)分别为(i,j,k)、(i,j-1,k)、(i,j+1,k)像素点处矢量在x轴方向的分量；V_y^(i,j,k)、V_y^(i,j-1,k)、V_y^(i,j+1,k)分别为(i,j,k)、(i,j-1,k)、(i,j+1,k)像素点处矢量在y轴方向的分量；V_z^(i,j,k)、V_z^(i,j-1,k)、V_z^(i,j+1,k)分别为(i,j,k)、(i,j-1,k)、(i,j+1,k)像素点处矢量在z轴方向的分量；T₀为阈值。

对于三维血管造影图像的z轴方向，判断某一像素点是否满足如下公式，如果满足，则该像素点为z轴方向的脊点：

其中，(i,j,k)为某一像素点的坐标；(i,j,k-1)和(i,j,k+1)为(i,j,k)在z方向的两个邻像素点；V_x^(i,j,k)、V_x^(i,j,k-1)、V_x^(i,j,k+1)分别为(i,j,k)、(i,j,k-1)、(i,j,k+1)像素点处矢量在x轴方向的分量；V_y^(i,j,k)、V_y^(i,j,k-1)、V_y^(i,j,k+1)分别为(i,j,k)、(i,j,k-1)、(i,j,k+1)像素点处矢量在y轴方向的分量；V_z^(i,j,k)、V_z^(i,j,k-1)、V_z^(i,j,k+1)分别为(i,j,k)、(i,j,k-1)、(i,j,k+1)像素点处矢量在z轴方向的分量；T₀为阈值。

同样，每一像素点处的矢量为归一化处理后的矢量，该矢量的模为1。阈值可根据提取精度进行选取，通常T₀选取[0,0.5]范围内的值。

一具体实施例中，如图5a、5b所示，采用本申请所述的基于图像梯度矢量流场的血管脊点提取方法提取图5a所示三维仿真图中的脊点，脊点的提取结果如图5b所示。

基于同一发明构思，本申请实施例中还提供了一种基于图像梯度矢量流场的血管脊点提取装置，如下面的实施例所述。由于该装置解决问题的原理与血管脊点提取方法相似，因此该装置的实施可以参见血管脊点提取方法的实施，重复之处不再赘述。

如图6所示，图6为本申请实施例的基于图像梯度矢量流场的血管脊点提取装置。该装置可以通过逻辑电路实现运行于智能终端，例如手机、平板电脑等设备中，或者以功能模块的方式由软件实现各部件的功能，运行于所述智能终端上。具体的，该装置包括：增强模块601，用于对血管造影图像中的血管目标进行增强；

梯度矢量流场计算模块602，用于利用梯度矢量流场模型求取增强后图像的梯度矢量流场；

脊点探测模块603，根据增强后图像的梯度矢量流场，分别计算增强后图像中的一像素点处矢量与该像素点沿一坐标轴方向的两个邻像素点处矢量夹角的余弦值，若所述余弦值满足阈值条件，则该像素点为该坐标轴方向的血管脊点。

进一步的，上述脊点提取装置还包括：剔除模块604，用于从得到的全部血管脊点中剔除孤立点。

实施时，所述剔除模块604具体用于判断某一脊点在以其为圆心，半径为d和d+d₀的两个同心圆所组成的圆环范围内是否存在其它脊点，如果不存在，则将该脊点为圆心，半径为d形成的圆范围的脊点删除，其中，d、d₀为距离常数。

一实施例中，所述增强模块601通过如下第一多尺度血管增强函数对三维血管造影图像中的血管进行增强：

v₀(s)为第一多尺度血管增强函数，R_A、R_B和S为测度函数，R_A用来区分片状和线状结构，R_B用来区分点状结构和线状结构，S用于区分背景像素，α、β和c为阈值，λ₁、λ₂和λ₃为Hessian矩阵H的三个特征值，且满足|λ₁|≤|λ₂|≤|λ₃|，D为图像的维度；

通过如下第二多尺度血管增强函数对二维血管造影图像中的血管目标进行增强：

其中，v₀'(s)为第二多尺度血管增强函数，R_B'和S'为测度函数，R_B'用来区分点状结构和线状结构，S'用于区分背景像素，β和c为阈值，λ₁和λ₂为Hessian矩阵H的二个特征值，且满足|λ₁|≤|λ₂|，D为图像的维度。

一实施例中，所述梯度矢量流场计算模块602具体用于：

求解如下表示梯度矢量流场模型的欧拉方程组：

得到增强后二维血管造影图像的梯度矢量流场两个分量相对于时间的迭代公式为：

其中，n为迭代次数，f_x(x,y)和f_y(x,y)为增强后二维血管造影图像的血管边缘函数f(x,y)在x和y方向的值，u(x,y)、v(x,y)为梯度矢量流场V(x,y)＝(u(x,y),v(x,y))的两个分量，(x,y)为增强后二维血管造影图像像素点的坐标，为laplacian算子，μ为控制参数；

同理，增强后的三维血管造影图像的梯度矢量流场三个分量的迭代公式为：

其中，n为迭代次数，f_x(x,y,z)、f_y(x,y,z)和f_z(x,y,z)为增强后三维血管造影图像的血管边缘函数f(x,y,z)在x、y和z方向的值，(x,y,z)为增强后三维血管造影图像像素点的坐标，为laplacian算子，μ为控制参数。

脊点探测模块603实现脊点探测过程可参见上述方法实施例，本申请在此不再赘述。

通过本申请上述实施例所述的技术方案能够有效压制背景噪声影响、提升血管造影图像中血管脊点的提取效率、增加血管目标脊点的探测数量及提高血管脊点探测精度，为后续血管中心线的提取和血管建模提供基础。

本申请实施例还提供一种电子设备，包括处理器及包括计算机可读程序的存储器，所述计算机可读程序在被执行时使所述处理器执行上面实施例所述的基于图像梯度矢量流场的脊点提取方法。

本申请实施例还提供一种计算机可读程序，其中当在电子设备中执行所述程序时，所述程序使得计算机在所述电子设备中执行如上面实施例所述的基于图像梯度矢量流场的脊点提取方法。

本申请实施例还提供一种存储有计算机可读程序的存储介质，其中所述计算机可读程序使得计算机在电子设备中执行上面实施例所述的基于图像梯度矢量流场的脊点提取方法。

应当理解，本申请的各部分可以用硬件、软件、固件或者它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可以用本领域共知的下列技术中的任一项或者他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述仅用于说明本申请的技术方案，任何本领域普通技术人员均可在不违背本申请的精神及范畴下，对上述实施例进行修饰与改变。因此，本申请的权利保护范围应视权利要求范围为准。