CN105303537A - 一种医学图像三维血管显示增强方法 - Google Patents

Abstract

本发明公开了一种医学图像三维血管显示增强方法，将血管增强看作一个滤波过程，以寻找符合人类血管管状特征的几何结构为主线。该方法首先是基于血管的形状特征，设计出一个类似管状的滤波器；然后是基于血管方向、直径的多样性，构造出一个方向和尺度均可调的滤波器组；最后是将滤波器组与待增强CTA/MRA图像进行空间卷积，经过一系列后处理得到最终增强图像。本发明公开的方法能够有效地改善血管的可视化效果，尤其是血管末梢的微小血管，以及对背景区域(比如相对平坦的面状区域、斑点状区域)与噪声的抑制，为构造一个性能优异的术中导航系统提供了坚实的基础。

Description

一种医学图像三维血管显示增强方法

技术领域

本发明涉及医学图像处理技术，尤其涉及CTA/MRA图像的3D血管显示增强方法。

背景技术

人体不同部位的血管可以提供人体相关组织的大量信息，其形态(主要包括血管的直径、血管的分叉角度、血管的弯曲度)、分布都是诊断相关血管疾病的重要指标。人类血管的准确可视化和精确量化对诊断和治疗血管疾病有着重要的意义，成为许多临床实践的重要前提。而血管的狭窄程度是血管疾病严重程度的重要指标，由于其决定着后续的具体治疗方案。不论是介入手术还是搭桥手术都需要一个性能优异的术中导航系统，其能帮助外科医生对三维血管结构有更好地观察。以上临床需求，对三维血管分割技术的完整性、准确度提出了很大的挑战，鉴于临床获取的医学数据低对比度、背景区域复杂、噪声大这一特点，在对三维血管实施分割技术之前，对其与周围背景区域的对比度进行良好地改善及背景区域和噪声的良好抑制是尤为重要的，这就需要一个鲁棒性好的三维血管增强算法，其性能将直接影响血管分割结果的精度、小血管在可视化中的显示效果。

目前，多种血管成像技术已被应用于临床实践中，如数字减影血管造影(DSA)、CT血管造影(CTA)、核磁共振血管造影(MRA)等，但是CTA/MRA数据较为普遍，所以基于该两类数据，发明了本方法，并在该两类数据上进行实践应用。

血管增强算法主要可以分为两大类：其一是固定某一尺度，构建一个由不同方向微分算子组合而成的滤波器进行滤波，简称为扩散滤波器，常见的有边界增强扩散(EED)、血管增强扩散(VED)，该方法由于尺度是固定的，对包含大尺度范围的血管区域存在问题；其二是基于血管函数来迭代多尺度，并通过分析Hessian矩阵的特征值来确定边缘，常见的有Frangi滤波器、Sato滤波器，该方法的问题是忽略了血管的方向。如果能将二者良好地结合，设计出一个多尺度、多方向的3D滤波器组，并利用该滤波器组对3D的CTA/MRA图像进行空间滤波，当尺度、方向与血管实际宽度、方向最匹配时，此滤波器组的输出最大，这样将成功地设计出一个性能优异的滤波器组，不仅能显示大血管，还能保留血管末梢的小血管结构，更能良好地抑制背景区域，解决了血管与其周围背景区域低对比度这一难题。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种医学图像三维血管显示增强方法，解决现有的CTA/MRA图像中3D血管显示增强效果差这一难题。

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

一种医学图像三维血管显示增强方法，包括如下步骤：

(1)构造多尺度、多方向滤波器组；

(2)设计符合血管管状的滤波器(类似于柱状的椭球体)，并确定滤波器的径向空间频率F；

(3)融合多尺度、多方向滤波器组的滤波结果，获取滤波后的响应图像，显示不同方向的大小血管。

所述步骤(1)中，多尺度、多方向滤波器组的构造方法为：

(11)首先构造n个尺度的滤波器小组，将第j个滤波器小组的尺度记为m_j，将m_j尺度滤波器的高斯标准差记为σ_j，将m_j尺度滤波器的径向空间频率记为F_j，设计径向空间频率记为F_j时，需使F_j的周期T稍微小于尺度m_j；所有n个滤波器小组的尺度为M＝[m₁,m₂,…,m_j,…,m_n]，各尺度下的高斯标准差为Σ＝[σ₁,σ₂,…,σ_j,…,σ_n]，各尺度下的径向空间频率为F＝[F₁,F₂,…,F_j,…,F_n]；

(12)第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个滤波器的方向，该方向对应的三维坐标旋转矩阵为；

$R_i={\left(\left[\begin{array}{ccc}{\mathrm{cos\&alpha;}}_i& -{\mathrm{sin\&alpha;}}_i& 0\\ {\mathrm{sin\&alpha;}}_i& {\mathrm{cos\&alpha;}}_i& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{\mathrm{cos\&beta;}}_i& 0& -{\mathrm{sin\&beta;}}_i\\ 0& 1& 0\\ {\mathrm{sin\&beta;}}_i& 0& {\mathrm{cos\&beta;}}_i\end{array}\right]\right)}^{-1}$

所述步骤(2)中，第j个滤波器小组中第i个滤波器的设计方法为：

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

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

$g_j\left(x_j,y_j,z_j\right)=\left(\frac{1}{{\left(2\mathrm{\&pi;}\right)}^{3/2}{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y{\mathrm{\&sigma;}}_z}\right)\exp \left(-\frac{1}{2}\left(\frac{{\left(x_j\right)}^2}{{\mathrm{\&sigma;}}_x^2}+\frac{{\left(y_j\right)}^2}{{\mathrm{\&sigma;}}_y^2}+\frac{{\left(z_j\right)}^2}{{\mathrm{\&sigma;}}_z^2}\right)\right)$

其中：σ_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＝σ_j；

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

${\mathrm{sr}}_{ij}(x_j,y_j,z_j)=cos\left(2{\mathrm{\&pi;F}}_j\sqrt{\left({\left(y_{ij}^{\&prime;}\right)}^2+{\left(z_{ij}^{\&prime;}\right)}^2\right)}\right)$

其中：P_j＝[x_j,y_j,z_j]^T，P_ij＝[x′_ij,y′_ij,z′_ij]^T，P_ij＝R_i×P_j；

(24)基于如下数学表达式，生成尺度为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)。

所述步骤(3)中，融合多尺度、多方向滤波器组的滤波结果，获取滤波后的响应图像，显示不同方向的大小血管，具体方法为：

(31)选定m_j尺度的l个滤波器与CTA/MRA图像做空间卷积，各个体素点选取与其方向最相符的滤波器的空间卷积结果(最大响应值)作为该尺度下的最终响应值，最终形成m_j尺度下的响应图像；

(32)由于待增强的图像经空间滤波后，其边缘的亮度非常高，对m_j尺度下的响应图像，采取去除边缘高亮区域和响应异常点赋0值处理，以便于进行归一化处理；

(33)由于在不考虑方向的前提下，大尺度滤波器的卷积结果会把小尺度滤波器的卷积结果覆盖掉，所以，要对所有n个尺度的最终响应值进行归一化处理，各个体素点选取与其尺度最相符的滤波器的空间卷积结果作为其增强结果值，以显示不同方向的大小血管。

本发明的滤波函数h_ij(x_j,y_j,z_j)主要由两部分构成，其一是具有平滑作用的高斯因子g_j(x_j,y_j,z_j)，其二是能够调节滤波器形状和方向的调和因子sr_ij(x_j,y_j,z_j)。选择高斯函数作为本方法滤波器的组成部分，原因是：在三维空间中，3D高斯函数在笛卡尔坐标系下将生成一个高斯球，该高斯球不仅具有旋转不变性这一优良性质，而且球体内部的体素点所对应的值呈现中心高亮、向周围逐渐对称降低，直至球面处接近于0的趋势。调和因子选择复指调和函数的实部特征，而没有考虑虚部及幅值特征的原因是，实部特征可看作为一个余弦函数，该函数的形态特征是轴对称，峰值出现在对称轴上，并向两侧对称递减，符合管状这种几何结构特征。

本方法中的径向空间响应频率F_j选取规则是，使余弦函数的周期稍小于滤波器模板的尺度，这样将高斯因子与调和因子(余弦函数)相乘，其滤波器模板中与余弦函数(T是周期)对应的元素是相对较高的正值，滤波器模板中与余弦函数 (T是周期)对应的元素是相对较低的负值，即相当于提高血管区域元素的权重、抑制血管周围区域元素的权重，这与实际中包含血管的CTA/MRA区域是相符的，给予高灰度值血管区域中的体素更高的正权值，相反的，给予稍低灰度值的血管周围区域更低的负权值，以使血管和非血管区域间的对比度更大，达到增强血管区域，抑制血管周围背景区域的目的。

另外，对卷积模板h_ij(x_j,y_j,z_j)进行归零处理的原因是，可以使滤波器在相对平坦区域的响应值更接近于0，达到对背景区域抑制的目的。

在本发明中，考虑到实际中的血管直径是各不相同的，比如说肾脏、心脏、腿部等区域的动脉都是大直径的，静脉的直径相对于动脉来讲，直径相对小些，而叶脉的直径就更小了，这样的话，一种尺度的滤波器是不能满足实际需求的；并且，实际中的血管整体类似一棵树，方向并不是单一的，会在相对较小的区域内，出现分叉，有的地方，甚至会出现螺旋状，比如脑部的血管，这样的话，一种方向的滤波器同样是不能满足实际需求的。因此，本方法考虑了尺度和方向，致力于构建一个多尺度多方向的滤波器组来满足对直径、方向各异的血管增强的目的。

本发明通过旋转(x_j,y_j,z_j)生成旋转坐标(x′_ij,y′_ij,z′_ij)，使得调和因子中母线方向是依赖方向角(α_i,β_i)的，母线方向决定了滤波器的长轴方向，母线也是sr_ij(x_j,y_j,z_j)中余弦函数的对称轴；由于调和因子和余弦函数的对称性，使得滤波器模板中任意一个垂直于长轴方向的平面，调和因子的元素值是以与长主轴的交点为中心，向周围对称递减，模板边界处元素值接近-1。

本发明通过设置不同尺度的滤波器模板，以增强不同直径的血管；基于步骤(2)中滤波器函数形式，应在不同尺度的滤波器模板中内切一个高斯球，但是考虑到在一个正方体中精确地内切一个球体，在实施上是具有一定困难，鉴于生成高斯球的函数实质上是三维标准正态分布的联合概率密度函数，以及真正决定一体素的周围参与其加权和的有效大小的因素是高斯标准差，因此本方法决定寻找联系滤波器尺度与高斯函数标准差之间的正比例因子，只要使足够多的体素落在该尺度下的滤波器模板内并且该模板的边界所对应的高斯因子值足够接近于0即可。

基于步骤(2)构造出的滤波器，其参数比较多，包括滤波器的尺度m_j、高斯标准差σ_j、径向空间频率F_j和方向(α_i,β_i)；发明人通过研究神经生理学相关的知识，发现人类视觉大脑皮层细胞对指定的物体有着对应的半幅值响应空频带宽B；并且径向空间频率F_j和高斯标准差σ_j通常有着如下关系：因此，发明人通过研究各参数间的关系及结合人工模拟数据加以验证，为添加造影剂的血管找出其对应的半幅值响应空频带宽B，成功将参数由四个减少到两个。

在方向给定的前提下，不同尺度滤波器各参数间的对应关系如下：

表1：参数表

经在人工模拟数据上的试验，发现两种方法得到的高斯标准差对最终响应图像，基本没有差别，以及从在肾脏CTA数据上的尝试结果来看，本发明成功地找到，添加造影剂的血管所对应的半幅值响应空频带宽B＝5，进而也从理论和实践的角度上，验证了径向空间频率F的选取规则和滤波器尺度与高斯标准差间正比例因子ratio＝2.5的合理性和最优性。

基于上述分析，我们给出两种参数限定方案以减少参数的个数：

A.所述滤波器的径向空间频率F_j与滤波器的尺度m_j之间存在如下关系所述滤波器的尺度m_j与高斯标准差σ_j之间存在如下关系ratio取值为2.5；

B.所述滤波器的径向空间频率F_j与滤波器的尺度m_j之间存在如下关系所述滤波器的径向空间频率F_j与高斯标准差σ_j之间存在如下关系B为添加造影剂的血管的半幅值响应空频带宽，且B的取值为5。

上述两种参数限定方案看似不同，但是通过在试验数据上的实验，验证了他们之间存在一定的相似或相关性，可以认为达到了一种殊途同归的效果。

有益效果：本发明提供的医学图像三维血管显示增强方法，构造出一组多尺度、多方向的类似管状的滤波器模板，使获取的CTA/MRA图像经该滤波器模板组空间滤波后，不仅能良好地增强血管区域，还能够保留微小的血管结构，更好地抑制背景区域与噪声，使各种血管更加清晰地显示；本方法直接从3D血管的结构特征入手，设计出了一个管状结构的血管滤波器，其通用性比较强，不再局限于人体某部位的血管，几乎适合于各种血管，只需要根据每部分血管的直径和方向复杂度来调用相应的滤波器即可；本方法相比较其它的方法，不再是将3D血管看成2D切片的叠加，其充分考虑了3D空间的特性，使其增强效果更加显著，甚至只通过简单的基于MIP的直接体绘制技术和窗技术就可以达到很好地可视化效果。

附图说明

图1为本发明的实现流程图；

图2为本发明的滤波器设计流程图；

图3(a)为符合血管亮度值分布的人工模拟数据；

图3(b)人工模拟数据经本方法中的滤波器组增强后的效果图；

图4为人工合成数据经本方法中的滤波器组增强后及区域增长后的效果图；

图5为人工合成数据添加不同程度噪声后的示例图；

图6为含有噪声的人工合成数据经本方法中滤波器组增强后的效果图；

图7为右肾脏CTA图像经阈值和本方法中滤波器组滤波后的结果示例；

图8为左肾脏CTA图像经阈值和本方法中滤波器组滤波后的结果示例；

图9为脑部MRA图像经阈值和本方法中滤波器组滤波后的结果示例。

具体实施方式

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

如图1所示为一种医学图像三维血管显示增强方法，首先是从人体3D血管的形状特征入手，设计出一个符合血管管状特征的滤波函数，并确定径向空间频率的选取规则，使在任意尺度的滤波器模板中，形成一个处于中心位置的高亮管状结构；然后是通过调节尺度，生成不同直径的管状结构，同时通过对3D坐标进行旋转，生成母线沿着不同方向的滤波器；最后是通过研究神经生理学相关的知识，找到添加造影剂的血管所对应的半幅值响应空频带宽B，将方向和尺度均可调的滤波器组的参数自由度降到2，即方向和尺度。这样，利用构造好的滤波器组与读入的CTA/MRA数据(必要时要经过一些预处理)进行空间卷积，并将每个滤波器输出的响应图像进行融合，以得到最终血管增强图像。

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

步骤一：设计符合血管形状特征的滤波器

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

(12)与步骤(11)在同一坐标下，基于复指调和函数的实部sr_ij(x_j,y_j,z_j)生成一个圆柱体，如图2(b)所示，其任意垂直于x轴的平面的元素值分布情况如图2(e)；

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

步骤二：构造满足血管直径和方向要求的滤波器组

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

(22)第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个滤波器的方向，该方向对应的三维坐标旋转矩阵为；

$R_i={\left(\left[\begin{array}{ccc}{\mathrm{cos\&alpha;}}_i& -{\mathrm{sin\&alpha;}}_i& 0\\ {\mathrm{sin\&alpha;}}_i& {\mathrm{cos\&alpha;}}_i& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{\mathrm{cos\&beta;}}_i& 0& -{\mathrm{sin\&beta;}}_i\\ 0& 1& 0\\ {\mathrm{sin\&beta;}}_i& 0& {\mathrm{cos\&beta;}}_i\end{array}\right]\right)}^{-1}$

(23)采用步骤一的方法设计出第j个滤波器小组中第i个滤波器。

步骤三、生成符合血管亮度值分布的人工模拟数据，确定空间滤波后各尺度、各方向下响应图像的融合规则；具体过程如下：

(31)人工模拟一组类似于血管的数据，该数据由8个尺度的圆柱体组成，直径依次是3、5、7、9、11、13、15、17，体素值的分布规律与增加造影剂的血管类似，中轴线上的亮度值最高，向两边对称递减，边缘处(即相当于血管壁)的宽度依次是1、2、3、4、5、6、7、8，边缘体素的亮度值相比较血管中心区域明显降低且更加接近背景区域，如图3(a)；

(32)将各尺度下，21个方向的滤波器组，分别对人工模拟数据进行滤波，分析模拟数据的响应值，发现如下规律：针对同一体素点，其响应值，会随着滤波器组尺度的增大而变大，并且图像边缘的响应值非常大，所以，采取去除边缘高亮区域并进行线性归一化的处理；重复本步骤，发现如下规律：不同直径的圆柱体总是对其相应尺度范围内滤波器组的响应值最大，尺度间没有重合现象。利用尺度分别为3、5、7、9、11、13、15、17的滤波器组对人工模拟数据进行增强组，得到结果，如图3(b)；

(33)对某个直径(比如说7)的圆柱体作旋转变换，使其母线方向与滤波器组某个方向相一致，并利用相应尺度(比如说13)的滤波器组进行空间滤波，对增强后的图像提取骨架，发现如下规律：该圆柱体经与其方向相一致的滤波器增强后，其最大响应值在骨架上；

(34)经步骤(32)和(33)后，原始图像经各尺度、各方向的滤波器空间卷积后，其响应图像的融合规则为：同一尺度下，选取最大的；不同尺度下，经归一化处理后，选取最大的。

步骤四、对包含管状、面片状、斑点状及不同程度高斯白噪声的人工合成数据使用本发明中的滤波器组对其增强；具体过程如下：

(41)对包含管状、面片状、斑点状结构的人工合成数据进行增强，原始图像如图4(a)，其中管状结构是由一个亮度值、直径、厚度均下降的3D简单螺旋线表示，斑点状结构是通过加性高斯随机场获取的；

(42)利用本发明中的滤波器组进行增强后的效果图如图4(b)，从图中可看出，本方法中滤波器组对尺度和方向的鲁棒性，及对包含斑点状及面片状结构的强抑制性，虽然，面片状的边缘不能被完全压制，主要原因是该原始图像的边缘亮度值分布比较像管状，所以本方法没有将其完全压制，同时也能看到，有一处出现了断裂，这是由于原始图像中，该处被一个高亮斑点状结构所覆盖，所以本方法对其进行了抑制；

(43)在增强结果的基础上，使用区域增长进行分割的结果图如图4(b)，从图中可以看出，面片状结构的边缘能够完全被去除，但是管状结构仍旧有微小断裂，对于这种断裂，可以通过步骤(44)中的策略解决；

(44)将滤波器组的尺度划分为几个范围，分别将不同尺度范围内滤波器组滤波后的结果进行区域增长，区域增长过程中，将小尺度范围内滤波器增强结果的种子点选在增强结果的小血管处，并且增强范围阈值要选的相对较小，将大尺度范围内滤波器增强结果的种子点选在增强结果的大血管处，并且增强范围阈值要选的相对较大，将分割后的结果进行融合，得到最终的分割结果，如图4(d)；

(45)对原始图像分别添加均值为0，不同方差的高斯噪声，添加方差分别为10、20、30、40、50、60高斯噪声后的原始图像如图5(a)～5(f)，经本发明中的滤波器组滤波后的图像如图6(a)～6(f)，其中前三幅图像使用的尺度是7、9、11、13、15、17，后三幅使用的尺度是9、11、13、15、17，从图8中不难看出，本方法对噪声能够良好地抑制。

步骤五：利用本发明中的滤波器组，对从外部机器设备上的获取的并经过基本预处理后的CTA图像或MRA图像进行增强处理，以说明本方法的广泛应用性，具体步骤如下：

(51)肾脏CTA图像的增强

肾脏作为人体主要的排泄器官，对维持机体内环境的平衡和稳定起到相当重要的作用。一旦肾脏的结构出现异常，将会导致一系列的临床症状，而肾脏中的血管作为运输通道之一，对其结构达到良好地可视化就显得尤为重要，同时，鉴于肾脏中的血管是穿透肾髓质、肾皮质等组织的，对其基于亮度值的增强更是带来了极大的困难，因此，这就更需要一个性能优异的三维血管增强算法。

人体的肾脏主要有两部分组成，分别是左肾脏和右肾脏，二者没有太大的差别，单从血管结构来讲，右肾的静脉相对于左肾短一些，也就是说，左肾的血管末梢更多一些。于是，本实验中，从获取到的CTA图像中，随机地选取一个人的左肾脏和右肾脏数据，具体处理步骤如下：

①双阈值预处理：通过对原始图像数据的分析，发现原始数据的灰度值范围非常大，动态范围达到2400多，而且血管区域的亮度值基本是在900-1500之间，考虑到后面多幅响应图像融合过程需要归一化处理，于是结合专家经验，将阈值分别设置为900和1500，具体如下：

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

③对增强后的图像进行三维可视化：采用的是最大亮度投影(MIP)方式进行可视化，并通过窗技术以达到更好的可视化效果。

人体右肾的处理结果如下：在窗位为0.7125，窗宽为0.79这一条件下，未经阈值处理前的图像如图7(a)，经阈值处理后的图像如图7(b)，经本发明中的方法增强后的图像如图7(c)，从图7中可以看出，本发明中的增强方法在保留血管末梢的细节部分的同时，对背景区域进行了良好地抑制。同样地，通过实践发现，若仅通过窗技术，要达到将背景区域全部压制掉的效果，不论是否经阈值处理，所得图像的血管区域都将会被抑制。

人体左肾的处理结果如下：在窗位为0.6025，窗宽为0.62这一条件下，未经阈值处理前的图像如图8(a)，经阈值处理后的图像如图8(b)，经本发明中的方法增强后的图像如图8(c)，从图中8(a)和8(b)的白色实线方框可以看到，在图8(c)中这种斑点状的背景区域已经全部被抑制，像白色虚线方框中的骨骼区域虽然没有被全部的抑制掉，但是也得到了部分的抑制，没有影响到后面白色箭头所指向血管末梢的可视化,综上，本发明中的增强方法在抑制背景区域的同时，良好地改善了血管与其周围背景区域的对比度及提高了血管末梢的可视化效果，及本滤波器对不同方向血管增强的鲁棒性。

(52)大脑MRA图像的增强

肾脏CTA图像的背景区域相对比较复杂，会包含一些骨骼、髓质组织等相邻解剖结构及获取过程中的噪声，尤其是骨组织，其灰度值与添加造影剂的血管基本在同一灰度范围内，但是其所包含的血管类似于一棵树，其结构相对简单。

接下来，我们来考虑一个背景相对比较简单，但是血管非常丰富的脑部MRA图像，脑血管一般包含拓扑和几何意义上的复杂结构，并被一些伪影所腐蚀。鉴于脑部MRA图像的动态范围并不是很大，并且基本上被血管区域所覆盖，所以本实验中并没有进行阈值预处理。图9中显示了在同一窗宽窗位下，基于MIP体绘制方式的原始图像和经本发明方法增强后的图像，图9(a)和图9(b)展示的是在同一视角下的三维图。从图中不难看出，原始图像的血管比较丰富，出现很多细小的血管，但是这些血管已经对主血管的可视化造成了影响，经过本发明中的增强后，发现对于一些过于细小的血管进行了压制，并提高了主血管的对比度，可以参见图9(a)和图9(b)中的白色箭头，血管的整体脉络更加地清晰可见。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (5)

1.一种医学图像三维血管显示增强方法，其特征在于：包括如下步骤：

(1)构造多尺度、多方向滤波器组；

(2)设计符合血管管状的滤波器，并确定滤波器的径向空间频率F；

(3)融合多尺度、多方向滤波器组的滤波结果，获取滤波后的响应图像，显示不同方向的大小血管。

2.根据权利要求1所述的医学图像三维血管显示增强方法，其特征在于：

所述步骤(1)中，多尺度、多方向滤波器组的构造方法为：

(11)首先构造n个尺度的滤波器小组，将第j个滤波器小组的尺度记为m_j，将m_j尺度滤波器的高斯标准差记为σ_j，将m_j尺度滤波器的径向空间频率记为F_j，所有n个滤波器小组的尺度为M＝[m₁,m₂,…,m_j,…,m_n]，各尺度下的高斯标准差为Σ＝[σ₁,σ₂,…,σ_j,…,σ_n]，各尺度下的径向空间频率为F＝[F₁,F₂,…,F_j,…,F_n]；

(12)第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个滤波器的方向，该方向对应的三维坐标旋转矩阵为；

$R_i={\left(\left[\begin{array}{ccc}{\mathrm{cos\&alpha;}}_i& -{\mathrm{sin\&alpha;}}_i& 0\\ {\mathrm{sin\&alpha;}}_i& {\mathrm{cos\&alpha;}}_i& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{\mathrm{cos\&beta;}}_i& 0& -{\mathrm{sin\&beta;}}_i\\ 0& 1& 0\\ {\mathrm{sin\&beta;}}_i& 0& {\mathrm{cos\&beta;}}_i\end{array}\right]\right)}^{-1}$

所述步骤(2)中，第j个滤波器小组中第i个滤波器的设计方法为：

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

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

$g_j\left(x_j,y_j,z_j\right)=\left(\frac{1}{{\left(2\mathrm{\&pi;}\right)}^{3/2}{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y{\mathrm{\&sigma;}}_z}\right)\exp \left(-\frac{1}{2}\left(\frac{{\left(x_j\right)}^2}{{\mathrm{\&sigma;}}_x^2}+\frac{{\left(y_j\right)}^2}{{\mathrm{\&sigma;}}_y^2}+\frac{{\left(z_j\right)}^2}{{\mathrm{\&sigma;}}_z^2}\right)\right)$

其中：σ_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＝σ_j；

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

${\mathrm{sr}}_{ij}(x_j,y_j,z_j)=cos\left(2{\mathrm{\&pi;F}}_j\sqrt{({\left(y_{ij}^{\&prime;}\right)}^2+{\left(z_{ij}^{\&prime;}\right)}^2)}\right)$

其中：P_j＝[x_j,y_j,z_j]^T，P_ij＝[x′_ij,y′_ij,z′_ij]^T，P_ij＝R_i×P_j；

(24)基于如下数学表达式，生成尺度为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)。

3.根据权利要求2所述的医学图像三维血管显示增强方法，其特征在于：所述滤波器的径向空间频率F_j与滤波器的尺度m_j之间存在如下关系所述滤波器的尺度m_j与高斯标准差σ_j之间存在如下关系ratio取值为2.5。

4.根据权利要求2所述的医学图像三维血管显示增强方法，其特征在于：所述滤波器的径向空间频率F_j与滤波器的尺度m_j之间存在如下关系所述滤波器的径向空间频率F_j与高斯标准差σ_j之间存在如下关系B为添加造影剂的血管的半幅值响应空频带宽，且B的取值为5。

5.根据权利要求2所述的医学图像三维血管显示增强方法，其特征在于：所述步骤(3)中，融合多尺度、多方向滤波器组的滤波结果，获取滤波后的响应图像，显示不同方向的大小血管，具体方法为：

(31)选定m_j尺度的l个滤波器与CTA/MRA图像做空间卷积，各个体素点选取与其方向最相符的滤波器的空间卷积结果作为该尺度下的最终响应值，最终形成m_j尺度下的响应图像；

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

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