基于对称性匹配滤波器组和区域生长的三维血管分割方法

摘要

本发明公开了一种基于对称性匹配滤波器组和区域生长的三维血管分割方法，首先设计一个满足三维血管的形状特征和灰度分布特点的管状结构的基本匹配滤波器模板，并用两个正交平面将基本匹配滤波器等分为四个象限；然后设计方向和尺度均可变化的匹配滤波器组，与待增强的医学图像进行空间卷积，确定以得到最大卷积响应值为融合规则，实现多尺度多方向的对称性匹配滤波器组的血管增强，得到增强图像，最后对增强图像使用区域生长提取三维血管。本发明公开的方法能够有效的对血管可视化效果进行增强，尤其是细小的分支末梢血管和靠近病变肿瘤的供血动脉有较好的增强作用，同时对于边缘皮质等杂质的抑制明显，使得血管和背景的对比度显著提升，为最终血管的提取打下坚实基础。

说明书

技术领域

本发明涉及医学图像处理领域的CTA图像的三维血管增强及分割技术，具体涉及一种基于对称性匹配滤波器组和区域生长的医学图像三维血管分割方法。

背景技术

血管疾病已经成为最重要的公共卫生问题之一，面对特定大小和复杂的血管造影图像，从中分割出血管是走向精确可视化、诊断和血管病变定量分析的关键一步。人体不同部位的血管可以提供人体相关组织的大量信息，其形态(主要包括血管的直径、血管的分叉角度、血管的弯曲度)、分布都是诊断相关血管疾病的重要指标。人类血管的准确可视化和精确量化对诊断和治疗血管疾病有着重要的意义，成为许多临床实践的重要前提。而血管的狭窄程度是血管疾病严重程度的重要指标，由于其决定着后续的具体治疗方案。不论是介入手术还是搭桥手术都需要一个性能优异的术中导航系统，其能帮助外科医生对三维血管结构有更好地观察。以上临床需求，对三维血管分割技术的完整性、准确度提出了很大的挑战，目前，多种血管成像技术已被应用于临床实践中，如数字减影血管造影(DSA)、CT血管造影(CTA)、核磁共振血管造影(MRA)等，但是CTA/MRA数据较为普遍，所以基于该两类数据，发明了本方法，并在该两类数据上进行实践应用。

由于CTA/MRA图像中三维血管的主干部分具有较高的灰度值，因此能够利用阈值分割法提取血管的主干部分，然而对于直径较小的细小血管分支末梢以及穿过皮质、病变肿瘤的血管，由于和背景的灰度值相近，使用一般的成像工具很难利用肉眼观察到，并且直接使用全局固定阈值的区域生长算法分割血管容易出现误分割和过分割的情况，影响分割效果。因此本发明提出了一种先增强后分割的方法，首先使用三维匹配滤波器组对原始图像进行增强显示，考虑到CTA/MRA图像边缘的高亮皮质区域也会产生较的响应值，为了抑制这种情况，本文利用了三维血管基于中心线的对称性信息，设计了基于对称性的改进匹配滤波器组与原图进行空间卷积，当滤波器的方向和尺度均最匹配时，输出卷积响应值作为增强结果。然后对血管增强，背景抑制的结果图使用全局阈值区域生长算法提取血管，并融合到原图得到最后的分割结果。

发明内容

发明目的：为了克服现有技术中存在的不足，有效对三维血管尤其是细小分支末梢进行增强显示和分割，本发明提供了基于对称性匹配滤波器组的三维血管增强方法，并结合区域生长算法分割得到完整的血管。

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

一种基于对称性匹配滤波器组和区域生长的三维血管分割方法，包括如下步骤：

步骤1：设计由三维高斯函数和复正弦函数实部乘积组成的符合血管管状结构的基本匹配滤波器模板；

步骤2：构建正交平面将基本匹配滤波器模板等分成四个象限的匹配滤波器模板之和；

步骤3：构造多尺度、多方向的对称性的匹配滤波器组，与待增强的图像进行空间卷积，对于尺度大于设定阈值的匹配滤波器，使用基本匹配滤波器减去四象限卷积两两做减法的绝对值，对于尺度不大于设定阈值的匹配滤波器，使用基本匹配滤波器减去对称象限的卷积响应的绝对值，同一尺度下的多个方向取最大值作为增强图像的卷积响应值，最后融合对称性匹配滤波器组的增强卷积响应，各个体素点选取与其尺度最相符的滤波器的空间卷积结果作为其增强结果值，得到增强图像；

步骤4：利用三维区域生长算法提取增强图像的血管，获取最终的分割结果。

具体的，步骤1中设计的在第j个尺度m_j下第i个方向的基本匹配滤波器模板h_ij(x_j,y_j,z_j)＝g_j(x_j,y_j,z_j)×sr_ij(x_j,y_j,z_j)；

其中(x_j,y_j,z_j)表示三维坐标，

表示尺度为m_j的3D高斯球的第j个点的值，其中σ_x表示g_j(x_j,y_j,z_j)在x轴上的标准差，σ_y表示g_j(x_j,y_j,z_j)在y轴上的标准差，σ_z表示g_j(x_j,y_j,z_j)在z轴上的标准差；

尺度为m_j、母线方向为(cosβ_icosα_i,cosβ_isinα_i，sinβ_i)的圆柱体的第j个点的值，其中F_j表示径向空间频率，β表示单位球面上某个离散点在球坐标系下的仰角，α是单位球面上某个离散点在球坐标系下的方位角。

具体的，所述步骤3中包括：

(3.1)选定m_j尺度(0<j<＝n)的l个方向的滤波器与CTA/MRA图像做空间卷积，根据图像实际血管尺度设定尺度阈值k，当m_j<＝k时，得到原图与对称性匹配滤波器的卷积响应值3DVSMFR_j为：

当m_j>k时，得到原图与对称性匹配滤波器的卷积响应值3DVSMFR_j为：

其中3DMFR_j表示基本匹配滤波器与原图的卷积响应，3DMFR_j^a、3DMFR_j^b、3DMFR_j^c、3DMFR_j^d分别表示四个象限的基本匹配滤波器与原图的卷积响应值，各个体素点选取与其方向最相符的滤波器的空间卷积结果作为该尺度下的最终响应值，最终形成m_j尺度下的响应图像；

(3.2)对m_j尺度下的响应图像，采取对于边缘高亮区域和响应异常点赋0值处理进行异常点的去除；

(3.3)对所有n个尺度的最终响应值进行归一化处理，各个体素点选取与其尺度最相符即卷积响应输出值最大的滤波器的空间卷积结果作为其增强结果值。

具体的，所述步骤4中包括：

(4.1)在增强结果图的血管主干入口处设置生长种子点；

(4.2)设置生长阈值范围；

(4.3)使用三维全局区域生长提取增强结果图的血管部分，并将分割结果作为掩膜获取原图中的血管分割结果。

有益效果：本发明提供了一种基于对称性匹配滤波器组和区域生长的三维血管分割方法，首先设计了由三维高斯函数和复正弦函数实部乘积组成的基本匹配滤波器，其在三维空间中的形状为管状与血管形状相似，使用此匹配滤波器与原图进行空间卷积滤波；为了抑制对边缘区域的错误增强，本发明提出了在基本匹配滤波器中加入对称性信息，将基本匹配滤波器在三维空间中等分成四个象限的匹配滤波器模块之和。考虑血管方向和尺度的多变性，本发明使用多尺度多方向的对称性匹配滤波器组，当尺度较小时使用对称象限的对称性信息，当尺度较大时使用四个象限的对称性信息进行增强，此方法在增强血管区域使其高亮的同时有效抑制了背景的边缘区域的响应，使得血管和背景组织的对比度进一步得到提升，最后使用区域生长算法提取得到完整的血管，结果表明对于小血管末梢的提取效果较好。本发明公开的方法能够有效对三维血管进行增强显示和分割，为以后的血管疾病临床诊断和手术提供基础信息。

附图说明

图1为本发明的方法流程图；

图2为基本匹配滤波器设计示意图；图2(a)-(c)分别为3D高斯球、圆柱体和管状滤波器的示意图，图2(d)-(f)分别为图2(a)-(c)截取x＝0平面对应的元素值分布图；

图3为本发明提出的四等分匹配滤波器的点乘因子示意图；图3(a)-(f)分别为用于不同象限之间做减法运算的点乘因子；

图4为人工合成数据增强效果图；图4(a)为模拟血管和边缘区域的人工合成数据；图4(b)人工合成数据经本方法对称性匹配滤波器组增强后的效果图；

图5为本发明方法用于肾脏CTA图像数据提取血管效果图；图5(a)为包含病变肿瘤的肾脏CTA图像数据，图5(b)为肾脏CTA图像数据经过本发明的对称性匹配滤波器组增强后通过最大密度投影法的三维显示，图5(c)肾脏CTA图像数据经过本方法提取得到肾动脉；

图6为本发明方法用于心脏CTA图像数据提取血管效果图；图6(a)为心脏CTA图像数据，图6(b)为心脏CTA图像经过对称性匹配滤波器组增强后通过最大密度投影法的三维显示，图6(c)为心脏CTA图像数据经过本方法提取得到的心脏冠状动脉。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如图1所示为一种基于对称性匹配滤波器组和区域生长的三维血管分割方法，首先是从人体3D血管的形状特征入手，设计出一个符合血管管状特征的滤波函数，并确定径向空间频率，使在任意尺度的滤波器模板中，形成一个处于中心位置的高亮管状结构；然后使用两个三维空间中的正交切面将此管状结构分割成四个象限部分，依次对每个象限赋予点乘因子1，其余象限赋予点乘因子0，获得每个象限对应的匹配滤波器部分模板。然后是通过调节尺度，生成不同直径的管状结构，同时通过对3D坐标进行旋转，生成母线沿着不同方向的滤波器；利用构造好的滤波器组与读入的CTA/MRA数据(必要时要经过一些预处理)进行空间卷积，对于尺度集合的对称滤波器组融合规则，当尺度较小时，使用中心对称的匹配滤波器部分模板的对称性与基本模板融合卷积，当尺度较大时，使用四个象限的匹配滤波器部分模板与基本模板融合卷积，完成所有匹配滤波器集合与原图的空间卷积并融合，以得到最终血管增强图像，对于增强图像，设定生长种子点和生长阈值规则，使用区域生长算法提取增强图像的三维血管，最后将其与原图进行掩膜处理得到最终的血管分割结果。

下面结合一个实施例对本发明作进一步说明，该实施例包括如下步骤：

步骤一：设计符合血管形状特征的滤波器并将其等分成四象限之和形式。具体包括如下步骤：

(1.1)在右手三维笛卡尔坐标系(x,y,z)上建立边长为m_j的正方体，该正方体的中心与坐标系的原点O重合，该正方体的边与坐标系的轴平行；将该正方体内第j个点P_j的三维坐标记为(x_j,y_j,z_j)，x_j,y_j,z_j均为整数，且

(1.2)基于如下数学表达式，生成尺度为m_j的3D高斯球，第j个点P_j对应的值为g_j(x_j,y_j,z_j)：

其中：σ_x表示g_j(x_j,y_j,z_j)在x轴上的标准差，σ_y表示g_j(x_j,y_j,z_j)在y轴上的标准差，σ_z表示g_j(x_j,y_j,z_j)在z轴上的标准差，且σ_x＝σ_y＝σ_z；

在三维空间笛卡尔坐标系下基于3D高斯函数g_j(x_j,y_j,z_j)生成的高斯球，如图2(a)所示，截取x＝0平面，其对应的元素值分布情况如图2(d)，其他垂直于x轴的平面的元素值分布情况类似于图2(d)；

(1.3)基于如下数学表达式，生成尺度为m_j、母线方向为(cosβ_icosα_i,cosβ_isinα_i，sinβ_i)的圆柱体，第j个点P_j对应的值为sr_ij(x_j,y_j,z_j)：

其中：P_j＝[x_j，y_j，z_j]^T，P_ij＝[x′_ij，y′_ij，z′_ij]^T，P_ij＝R_i×P_j，其中表示单位球面上某个离散点在球坐标系下的仰角，α(-π≤α≤π)是单位球面上某个离散点在球坐标系下的方位角。其中j表示第j个尺度，i表示该尺度下第i个方向，F_j表示径向空间频率，

基于复指调和函数的实部sr_ij(x_j,y_j,z_j)生成的圆柱体，如图2(b)所示，其任意垂直

于x轴的平面的元素值分布情况如图2(e)；

(1.4)基于如下数学表达式，生成尺度为m_j、主轴方向为(cosβ_icosα_i,cosβ_isinα_i，sinβ_i)的匹配滤波器模板h_ij(x_j,y_j,z_j)：

h_ij(x_j,y_j,z_j)＝g_j(x_j,y_j,z_j)×sr_ij(x_j,y_j,z_j)

对h_ij(x_j,y_j,z_j)进行归零处理，得到第j个滤波器小组中第i个滤波器，该滤波器的尺度为m_j，方向为(α_i,β_i)。

基于滤波器函数h_ij(x_j,y_j,z_j)，生成符合血管管状特征的滤波器，如图2(c)，截取x＝0平面，其对应的元素值分布情况如图2(f)，其他垂直于x轴的平面的元素值分布情况类似于图2(f)；

(1.5)基于(1.4)中的管状结构的滤波器，以滤波器在三维空间右手笛卡尔坐标系中的中心点p₀(x₀,y₀,z₀)和滤波器方向(v_x,v_y,v_z),选取空间中的随机一点P₁(x₁,y₁,z₁)，过该点和中轴线得到第一个切面，该切面的法向量N₁＝(n_1x,n_1y,n_1z)：

N₁＝(x₁-x₀,y₁-y₀,z₁-z₀)×(v_x,v_y,v_z)

得到平面A的点法式方程：

n_1x(x-x₀)+n_1y(y-y₀)+n_1z(z-z₀)＝0

由于两个切面相互垂直且都过中轴线，因此第二个切面B的法向量N₂＝(n_2x,n_2y,n_2z)可以由向量N₁和向量V的叉积获得：

N₂＝(n_1x,n_1y,n_1z)×(v_x,v_y,v_z)

得到平面B的方程如下：

n_2x(x-x₀)+n_2y(y-y₀)+n_2z(z-z₀)＝0

由平面A和平面B将滤波器管状结构等分为a、b、c、d四个象限，将步骤1中的基本匹配滤波器h_ij(x_j,y_j,z_j)分解成四个象限的匹配滤波器之和：

h_ij(x_j,y_j,z_j)＝h_{ij_a}(x_j,y_j,z_j)+h_{ij_b}(x_j,y_j,z_j)+h_{ij_c}(x_j,y_j,z_j)+h_{ij_d}(x_j,y_j,z_j)

其中h_{ij_a}(x_j,y_j,z_j)、h_{ij_b}(x_j,y_j,z_j)、h_{ij_c}(x_j,y_j,z_j)、h_{ij_d}(x_j,y_j,z_j)分别表示a、b、c、d四象限的匹配滤波器模板。

步骤二：设计多尺度多方向的对称性匹配滤波器组

(2.1)设计8个尺度的滤波器小组，M＝[3,5,7,9,11,13,15,17]；

(2.2)第j个滤波器小组包括l个滤波器，l个滤波器分别对应l个方向，l个方向按照如下步骤确定：

在单位球体的球面上均匀选取2l个点，鉴于球体的对称性，选取出半个球面上的l个点；以单位球体的球心为原点O，建立右手三维笛卡尔坐标系(x,y,z)，球面上第i个点P_i的三维坐标记为(x_i,y_i,z_i)，1≤i≤l；

将三维坐标(x_i,y_i,z_i)转换为球坐标(α_i,β_i)，α_i表示OP_i在x-y面的投影与x轴正半轴的夹角，β_i表示OP_i与x-y面的夹角，-π＜α_i＜π，(α_i,β_i)决定第i个滤波器的方向，该方向对应的三维坐标旋转矩阵为；

(2.3)对于第j个滤波器小组中第i个基本滤波器以及a、b、c、d四象限中的滤波器部分模板由步骤一得出。

步骤三、确定融合规则，实现多尺度多方向的对称性匹配滤波器组的血管增强，具体步骤如下：

(3.1)对于选定m_j尺度的l个滤波器与CTA/MRA图像做空间卷积，当m_j<＝5时，得到原图与对称性匹配滤波器的卷积响应3DVSMFR_j为：

当5<m_j<＝17时，得到原图与对称性匹配滤波器的卷积响应3DVSMFR_j为：

其中3DMFR_j表示基本匹配滤波器与原图的卷积响应，3DMFR_j^a、3DMFR_j^b、3DMFR_j^c、3DMFR_j^d分别表示四个象限的基本匹配滤波器模板与原图的卷积响应值。其中3DMFR_j^a-3DMFR_j^b、3DMFR_j^a-3DMFR_j^c、3DMFR_j^a-3DMFR_j^d、3DMFR_j^b-3DMFR_j^c、3DMFR_j^b-3DMFR_j^d、3DMFR_j^c-3DMFR_j^d由图3中的点乘因子和3DMFR点乘获得，各个体素点选取与其方向最相符的滤波器的空间卷积结果作为该尺度下的最终响应值，最终形成m_j尺度下的响应图像；

(3.2)对m_j尺度下的响应图像，采取去除边缘高亮区域和响应异常点赋0值处理；

(3.3)然后对所有n个尺度的最终响应值进行归一化处理，各个体素点选取与其尺度最相符的滤波器的空间卷积结果作为其增强结果值，以显示不同方向的大小血管。

基本匹配滤波器的对于管状结构的组织有较好的增强作用，但与此同时，对于平面结构的边缘部分也有较高的增强响应值，这对于后续提取血管的后处理步骤有较大的干扰。本发明提出的对称性匹配滤波器组调到血管增强方法在对血管中心点实行增强显示的同时，能够抑制平面结构边缘的错误响应，使得血管和背景的对比度提升，对后续血管分割效果的提升有明显帮助。图4(a)为模拟血管、半径分别为1和2的人工合成管状圆柱和模拟肾皮质的平面结构，图4(b)为使用对称性匹配滤波器组进行增强的结果图，结果对比显示，圆柱截面中心点增强的同时，对称性匹配滤波器组对于平面边缘的抑制效果较好。

步骤四：对增强图像使用区域生长提取血管，并融合到原图中得到最终的分割结果；具体过程如下：

(4.1)在增强结果图的血管主干入口处设置生长种子点P₀。

(4.2)设置生长阈值范围[T₁，T₂]，其中T₁和T₂的设定原则是在提取到血管的同时不误提取任何皮质、骨组织等背景。

(4.3)使用三维全局区域生长提取增强结果图的血管部分，并将分割结果作为掩膜获取原图中的血管分割结果。

下面利用本发明中的提出的对称性匹配滤波器组，对从外部机器设备上的获取的并经过基本预处理后的CTA图像或MRA图像进行增强处理，然后利用区域生长算法对血管进行分割，以说明本方法的广泛应用性。

1、病变肾脏CTA图像的肾动脉的增强与分割

准确实施肾动脉分支阻断肾部分切除手术提供可靠保障重要的一环是清晰显示肾动脉的结构，这也是为医师对血管造影进行观察进行相关判断和确认所必需的，肾脏动脉的主干尺寸较大，具有较高的灰度值，与背景组织区分明显，能够借助三维影像显示平台较为清楚地显示，而对于直径较小、灰度值较低的肾动脉分支末梢以及在病变肿瘤附近的细小血管则是本发明的增强与分割重点。

具体处理步骤如下：

①截取包含左侧肾脏的感兴趣区域(ROI)区域(如图5(a)所示)，对其进行重采样，保证采样后图像在各方向上的分辨率相同，为0.5mm。

②本发明中的对称性滤波器组进行滤波增强；

③对增强后的图像选择肾动脉主干一点作为生长种子点，设定生长灰度阈值范围为[70,600]，使用区域生长算法得到三维血管的分割结果，与原图进行掩膜处理获得最后的三维血管分割结果。

病变肾脏CTA图像肾动脉增强与分割结果说明：图5(a)和图5(b)显示了在同一窗宽窗位下，基于最大密度投影体绘制方法的原始图像和本文发明方法增强后的图像，可以看到，增强后图像对于肾脏动脉有其是灰度值较低、直径较小的肾动脉分支末梢和肾肿瘤附近的供血动脉有较好的增强作用，同时，对肾皮质边缘的抑制较好，没有让肾皮质区域产生较高的错误响应。图5(c)显示了使用区域生长算法对增强结果图中对血管的分割结果，同样的，对于肾脏动脉末梢和肿瘤供血动脉的提取效果明显。

2、心脏CTA图像中心脏冠状动脉的增强与分割

随着冠心病的发病率的居高不下，对于心脏CTA图像，从中提取冠状动脉树，可以使医生对心血管疾病进行分析，及时发现病变。图6(a)为截取包含心脏冠状动脉的ROI区域的MIP体绘制结果，图6(b)为使用本发明提出的对称性匹配滤波器组增强后的MIP体绘制结果，图6(c)为使用区域生长算法对增强图中冠状动脉的提取结果，可以看到，经过本发明提出的对称性匹配滤波器增强后的心脏CTA图像，冠状动脉和背景组织的对比度明显提高，灰度值差异明显，使用区域生长算法完成了对于冠状动脉树的准确提取。