基于管状特征增强滤波和脊线跟踪的三维血管中轴线提取方法与流程

本发明属于医学图像特殊组织分析与提取方法领域，涉及一种三维血管中轴线提取方法。
背景技术：
：目前对血管相关的医学图像处理，首先运用成熟的方法对相关进行分割，针对不同数据的特点运用不同的方法，针对血管造影成像进行分析处理，血管造影数据在图像上的血管部分，往往表现具有“高亮”的特征，即其值相对局部范围比较高。针对这一特点，目前对医学图像的三维血管中轴线提取的技术主要分为以下三大类：第一类、形态学处理方法即用形态学的腐蚀方法，对预处理过的血管数据进行三维细化。由Palágyi,Kálmán、Kuba,Attila等人提出最初的三维细化算法，而后许多研究都在此框架上提出自己的改进方法，针对不同的细化对象进行了改进。第二类、距离变换方法对体数据进行距离变换获得距离图，从而提取血管的中轴线。HeTS等人针对虚拟内窥镜的特点提出了ReliablePath的中心路径提取算法、WanM等人和BitterI等人，提出了各自距离变换的算法，应用在医学图像的拓扑细化中。第三类、模型演化法利用轴点检测，或者上述两种方法获得初始的中轴线，再根据现有模型进行迭代，使其能量达到最优，从而得到最合适的中轴线。K.Krissian等分析了圆柱模型的特点提出了响应的血管分割算法并提取出中心线。此外还有大量相关研究在此类方法上有深入地探讨。以上三种方法各有各的特点，形态学的方法能够得到比较连续的中轴线，但它的结果是基于原始图像的分割结果，若是分割结果不精准或者存在一些噪点伪影等干扰，容易生成很多不正确的或多余的轴线；距离变换方法受原图干扰较小且具有一定的旋转不变性，缩放不变性和平移不变性，获得的中轴线位置准确性高，但是中轴线大多不连续，需要进行后续的处理，且在三维情况中距离变换的方法相对比二维的情况更复杂；模型法需要一个初始轴线，所以可能引入误差，且模型的参数及演化的方法都需进一步考量，算法的复杂度高。随着计算机技术和医学成像技术日新月异的巨大发展，医学图像在对辅助医生对病人病情诊断和治疗的发挥其越来越大的作用。在医学图像中存在丰富的信息反映病人的身体情况，然而目前医生只通过自身的经验和软件简单的显示来观察，容易受到图像自身噪声或者医生主观的判断的影响，可能产生不准确诊断结果。计算机医学图像处理技术，作为在计算机和医学之间的一个交叉学科，通过对医学和成像知识和计算机图像处理两个方面对医学图像进行分析，能够有效地利用医学图像反映的信息，把它以更直观更全面的方式展现处理，从而达到辅助医生进行更有效更准确地诊断的目的。由于成像方式的不同，所得到的医学图像数据也有各自的特点，目前研究较多的，一种是核磁共振成像MRI，它的数据分辨率高，图像质量较好，并且根据不同的信号加权可以获得很多类型的成像，便于医生和学者进行分析研究；另一种则是计算机断层扫描成像CT，它的成像速度相比MRI要快很多，而且成本较低，各个组织成分的CT值在CT数据中都有一定的范围，但是，它也有成像质量一般，容易受到干扰产生伪影等一系列的缺点。针对不同成像方式得到的医学图像数据，国内外学者用不同的方法去分析他们，达到的想要的结果。技术实现要素：为了克服已有三维血管中轴线提方式的准确性较差的不足，本发明提供了一种准确性较好的基于管状特征增强滤波和脊线跟踪的三维血管中轴线提取方法。本发明解决其技术问题所采用的技术方案是：一种基于管状特征增强滤波和脊线跟踪的三维血管中轴线提取方法，所述提取方法包括以下步骤：1)血管造影图像的预处理，过程如下：1.1)最大密度投影处理四维脑部血管造影数据；1.2)多尺度管状增强滤波；1.3)误增强部分剔除；2)根据血管方向信息，加上血管中心值是局部最大的特点，采用类似脊线跟踪的方法来提取三维脑补血管中轴线，过程如下：2.1)确定起始跟踪位置：在初始化数据后选择全局最大值做为血管跟踪的起始位置；2.2)消除非轴点：跟踪一根血管后，利用三维区域生长的方法，将该血管所处于的局部区域内的非轴点体素进行消除；2.3)继续跟踪或者停止跟踪：在此找到下一全局最大值做为下一根血管的初始跟踪位置，若全局最大值低于设定的门限值，则表示整个数据已经处理完毕，停止跟踪。进一步，所述步骤2.1)的过程如下：2.1.1)正方向的跟踪：根据当前点位置想pos(x0,y0,z0)和当前方向E1(a,b,c)，构建其法平面方程a*(x-x0)+b*(y-y0)+c*(z-z0)＝0，设定一个局部范围找到当前位置所处的局部平面；2.1.2)判断是否处于当前平面的局部最大值，若不是，则进行局部最大偏移校正；若是，则根据血管方向E1确定下一跟踪位置搜索区域；2.1.3)在下一跟踪位置搜索区域中找到最大值，判断是否是图像边缘或者与当前位置的值相差过大或者方向差异在一定范围内，若符合要求，该位置是跟踪的下一位置nextpos则返回2.1.1)步骤继续跟踪；若不符合，则换另外以位置继续判断；如果所有位置都不符合，停止正方向的跟踪；2.1.4)逆方向跟踪：与正方向相似，只把跟踪的轴向旋转180°。再进一步，所述提取方法还包括以下步骤：3)血管中轴的分叉和连结处理：整个血管连结操作按中轴线的首部、尾部分为以下两端点的两种操作；对于普通连结情况：3.1)首部连结尾部的情况，2个连接点构成新的血管，将2个中轴线合并；3.2)首部连结首部的情况，2个连接点构成新的血管，将2个中轴线合并，且将被连结的中轴线的所有轴点方向进行反置；3.3)尾部连结首部，2个连接点构成新的血管，将2个中轴线合并；3.4)尾部连结尾部，2个连接点构成新的血管，将2个中轴线合并，且将被连结的中轴线的所有轴点方向进行反置；对于分叉处理情况：需补齐相的轴点；当我们找到2个端点的对应关系时，获得当前位置的一个方向，根据方向进行逐步跟踪来补全中轴线。一直遍历之前得到的初始中轴线，直到不会产生新的轴线的连结操作，完成全局血管中轴提取。更进一步，所述步骤1.1)中，从一个时间序列中获得的不同时刻的造影数据，由于动脉血流在造影剂的作用下能够在成像中反应出高值的特性，利用最大密度投MIP进行时间维度的合并，同时分析各个时间点每个体素的值的变化率，剔除由于成像时造成的伪影的造影数据，从而获得完整的一个血管造影数据。所述步骤1.2)中，针对不同半径大小的血管，利用尺度空间理论，不同半径的血管在符合它尺度下的滤波能够取得相应的最大响应值，通过不同尺度的高斯滤波，构造多尺度的管状增强滤波。所述步骤1.3)中，通过对每一切层的梯度幅值图像，进行最大类间分割得到梯度幅值的门限值，并且根据原数据建议的窗宽窗位值设定原数据的一个门限值，将获得的目标区域进行剔除，目标区域即误增强的区域，再最终结果归一化并与原数据点乘，得到待提取中轴的数据。本发明针对人脑进行计算机断层扫描得到的血管造影CT，提取血管中轴线并跟踪血管走向，从而分析整个脑部血管树，便于后续的脑部血管三维重建，血流分析等提供可靠且有效的信息。此外，血管的中轴线能够对血管进一步进行精确分割或者不同模态图像的精确配准提供了一定的拓扑先验，能够是这方面的工作更加有效的进行。本发明的有益效果主要表现在：准确性较好。附图说明图1是在时间序列上同一位置的横断面脑补血管造影原数据，最大密度投影法的运用的示意图。图2是不同时刻同一位置的像素变化情况的示意图。图3是图2中data1、data2、data3、data4的示意图，其中，(a)表示data1，(b)表示data2，(c)表示data3，(d)表示data4。图4是圆柱血管模型中心体素的示意图，E1表示最小幅值特征所对应的方向(血管轴向)，E2和E3是剩下两个特征向量，构成了血管的垂直平面。图5是不同尺度下同一横断面的数据，(a)为MIP的结果，(b)多尺度管状增强滤波后进行尺度融合的结果，可以看到误增强了最外一圈，(c)是最大尺度响应图，(d)是进行误增强剔除后的结果，(e)是不同尺度下管状增强滤波后的结果，(f)是不同尺度下的梯度幅值。图6是血管三维跟踪提取中轴线流程图。图7是血管的两种连结情况，(a),(c)表示普通断裂的连结，(b),(d)表示分叉的连结。具体实施方式下面结合附图对本发明作进一步描述。参照图1～图7，一种基于管状特征增强滤波和脊线跟踪的三维血管中轴线提取方法，所述提取方法包括以下步骤：1)血管造影图像的预处理，过程如下：1.1)最大密度投影处理四维脑部血管造影数据；1.2)多尺度管状增强滤波；1.3)误增强部分剔除；步骤1.1)中，从一个时间序列中获得的不同时刻的造影数据，由于动脉血流在造影剂的作用下能够在成像中反应出高值的特性，利用最大密度投影(MIP)进行时间维度的合并(图1)，同时分析各个时间点每个体素的值的变化率，剔除由于成像时造成的伪影的造影数据(图2)，从而获得完整的一个血管造影数据。从图2的曲线分析图中可以观察到，造影部位，普通部位和伪影部位的值随时间变化的情况，造影的部位(Data2)出现类似抛物线的均匀变化，普通组织部位(Data3)和背景部分(Data4)的值基本不变，伪影部位(Data1)在某一时刻产生突变，而且那些时刻造影剂的部分已经失效，那些部位的CT值和普通部位接近。因此我们可以提出出现伪影的那些时间刻的成像，在进行最大密度投影操作。步骤1.2)中，除去原数据中非血管的部分，增强原数据中的血管部分把图像的每个体素看成一个函数I(x,y,z)，将它进行泰勒二阶展开得到I(χ0+Δχ)≈I(χ0)+ΔχT▽I(χ0)+ΔχTH(χ0)Δχ其中▽I(χ0)代表梯度，H(χ0)代表曲率，反应了梯度的变化，每个体素经过不同维度的二次微分都能够得到6种二阶导，它们构成了体素的Hessian矩阵通过求解该矩阵的特征方程，能够得到3组个特征值和它们相对应的的特征向量，将它们按照幅值从小到大排列，表示为λ1、λ2、λ3(|λ1|＜|λ2|＜|λ3|)，通过分析这3个特征值的大小，可以找到一下规律，如表1：λ1λ2λ3结构LLL噪声或无结构LLH+黑暗的片状结构LLH-明亮的片状结构LH+H+黑暗的管状结构LH-H-明亮的管状结构H+H+H+黑暗的团状结构H-H-H-明亮的团状结构表1在血管造影CT中，血管结构在图上表现出明亮的管状结构，我们可以利用|λ1|＜＜|λ2|＜|λ3|和λ2＜0、λ3＜0这两个条件判读体素是否是明亮的管状结构。根据Krissian等人的研究分析，绝对值最小的特征值所λ1对应的特征向量e1代表了管状结构中轴线的切线方向，而其他两个绝对值较大的特征值λ2、λ3构成了与中轴线相垂直的平面，而且特征值的绝对值越大，代表越接近管的中心。如图4所示。根据以上规律，Frangi等人研究设计了血管增强滤波，公式如下：α和β代表两个项的敏感系数，一般取值0.5，C代表噪声控制系数，取全局最大值的四分之一此外，针对不同半径大小的血管，利用尺度空间理论，不同半径的血管在符合它尺度下的滤波能够取得相应的最大响应值，我们加入了多尺度的概念，通过不同尺度的高斯滤波，构造了多尺度的管状增强滤波经过1.2)的操作，大部分的非血管结构被剔除了，但是在图像边缘由于梯度剧烈变化的原因会造成相应部分的误增强，我们通过对每一切层的梯度幅值图像，进行最大类间分割得到梯度幅值的门限值，并且根据原数据建议的窗宽窗位值设定原数据的一个门限值，将获得的目标区域(即误增强的区域)进行剔除，再最终结果归一化并与原数据点乘，得到待提取中轴的数据。2)根据图4E提供的血管方向信息，加上血管中心值是局部最大的特点，本发明设定了一种类似脊线跟踪的方法来提取三维脑补血管中轴线，主要分为以下几步：2.1)确定起始跟踪位置：在初始化数据后选择全局最大值做为血管跟踪的起始位置；2.1.1)正方向的跟踪：根据当前点位置想pos(x0,y0,z0)和当前方向E1(a,b,c)，构建其法平面方程a*(x-x0)+b*(y-y0)+c*(z-z0)＝0，设定一个局部范围找到当前位置所处的局部平面2.1.2)判断是否处于当前平面的局部最大值，若不是，则进行局部最大偏移校正；若是，则根据血管方向E1确定下一跟踪位置搜索区域2.1.3)在下一跟踪位置搜索区域中找到最大值，判断是否是图像边缘或者与当前位置的值相差过大或者方向差异在一定范围内，若符合要求，该位置是跟踪的下一位置nextpos则返回2.1.1)步骤继续跟踪；若不符合，则换另外以位置继续判断；如果所有位置都不符合，停止正方向的跟踪2.1.4)逆方向跟踪：与正方向相似，只把跟踪的轴向旋转180°；2.2)消除非轴点：跟踪一根血管后，利用三维区域生长的方法，将该血管所处于的局部区域内的非轴点体素进行消除；2.3)继续跟踪或者停止跟踪：在此找到下一全局最大值做为下一根血管的初始跟踪位置，若全局最大值低于设定的门限值，则表示整个数据已经处理完毕，停止跟踪；流程图如图6。进一步，所述提取方法还包括以下步骤：3)血管中轴的分叉和连结处理：步骤2)能够快速获得初始的中轴线，但是在分叉处由于多尺度数据和三维区域生长终止的原因造成了断裂，还有其他一些由于细小血管或者曲率很大的血管造成了跟踪失效导致中轴线断裂，针对上述情况进行分叉和连结处理，如图6所示。整个血管连结操作按中轴线的首部、尾部分为以下两端点的两种操作；对于普通连结情况：3.1)首部连结尾部的情况(a)，2个连接点构成新的血管，将2个中轴线合并；3.2)首部连结首部的情况(c)，2个连接点构成新的血管，将2个中轴线合并，且将被连结的中轴线的所有轴点方向进行反置；3.3)尾部连结首部，2个连接点构成新的血管，将2个中轴线合并；3.4)尾部连结尾部，2个连接点构成新的血管，将2个中轴线合并，且将被连结的中轴线的所有轴点方向进行反置；对于分叉处理情况：无论是首部连结分叉处(另一中轴线的中间部分)还是尾部连结分叉处，都只需补齐相的轴点；当我们找到2个端点的对应关系(即图中五角星点和虚线箭头所示的待连结点)时，我们可以获得当前位置的一个方向(两位置坐标相减)，根据方向进行逐步跟踪来补全中轴线。整个过程操作一直遍历之前得到的初始中轴线，直到不会产生新的轴线的连结操作，全局血管中轴提取就完成了。当前第1页1 2 3