CN101551862B - 基于边界距离场的血管中心路径提取方法 - Google Patents

Abstract

本发明公开了一种基于边界距离场的血管中心路径高精度快速提取方法，主要解决在CTA血管造影体数据中血管提取速度慢和中心路径提取精度不高的问题。该方法包括：利用八叉树对体数据进行分块，并记录所有分块信息；利用每一分块的最大值和最小值与双阈值的关系，对血管进行提取和边界距离场计算；用原CTA数据中血管体素的梯度值倒数和拉普拉斯变换值之和作为边界距离初始值；建立基于边界距离场的最大代价生成树，并提取树的主干作为血管的初始中心路径；用重心法修正初始中心路径，使路径上的点都处于血管管腔横截面的重心处，得到最终的血管中心路径。本发明具有血管提取速度快，中心路径提取精度高的优点，适用于各种管腔组织的中心路径提取。

Description

基于边界距离场的血管中心路径提取方法

技术领域

本发明属于生物医学技术领域，特别涉及血管的中心路径提取方法，可用于各种管腔组织的中心路径提取，并应用到虚拟内窥镜中。

背景技术

医学图像中的血管精确提取和描述具有广泛的应用，如梗塞或者狭窄化的血管疾病诊断及不同时段的病人图像配准等；血管的提取方法是计算机血管辅助诊疗系统的关键。虽然目前存在很多种血管提取方法，但因成像模式、应用领域、人机交互以及其他因素的影响而各不相同，目前尚没有一种通用的提取方法能够适用于所有模式的医学图像。可视化技术的发展使得人们可以对所获得的三维数据通过各种重建方法进行显示以获取有用的信息。对于血管疾病诊断而言，利用可视化技术对提取后的血管进行三维显示，不仅能够对血管外表面进行观察，还可以对血管内部进行观察，即血管虚拟内窥镜。要实现血管虚拟内窥镜，必须提取血管的中心路径，即自动漫游路径。中心路径的提取方法一般分为三类：人工指定，拓扑细化和距离变换。

人工指定方法要求用户在体数据切片上人工指定路径中心点，所有的点连起来组成中心路径。这与人工漫游一样，虽然比较费时，但精确度较高，能够更准确的接近感兴趣区域。拓扑细化又称为剥洋葱法，是一个较为传统的方法，主要思想是通过反复剥除物体最外层体素，且每一步只删除不会影响模型拓扑结构的点，直到只剩下单体素宽度的骨架为止，将所获得的骨架作为中心路径。拓扑细化是一个迭代过程，它能够保持模型的整体连接性，但计算量十分庞大，而且不能直接连通到任意给定的两点，即对于不在骨架上的点，需要事先连接到骨架上。距离变换主要基于这样的事实：物体的中轴是物体内部到边界的最大距离点的集合。这种方法通常是先对二值体数据进行距离变换，然后将变换后的距离数据看作一个有向或无向的加权图，最后采用Dijkstra最短路径生成算法得到中心路径。距离变换法的速度较快，但在保持模型的拓扑结构方面不如拓扑细化法。

用于血管漫游的中心路径提取多采用距离变换法，原因主要有两点：一是医学模型通常数据量较大，采用拓扑细化法耗时太多；二是提取中心路径是为了在漫游时指导视点的移动，其主要目的在于更好的观察血管内部结构。目前存在很多基于距离变换的中心路径提取算法。德国的Marburg大学数学与计算机学院的Jan Egger等在CT血管数据中利用距离变换法提取血管的骨架作为虚拟漫游路径，取得了较好的效果，见文献“Egger J，

Z，Groβkopf S，et al.A fast vessel centerline extraction algorithm for cathetersimulation.Twentieth IEEE International Symposium on Computer-Based Medical Systems，Maribor，Slovenia.2007.”。Bitter等和Sato等用边界距离场和Dijkstra最短路径提取算法提取结肠等数据的中心路径，见文献“Sato M，Bitter I，Bender MA，et al.TEASAR：Tree-structure extraction algorithm for accurate and robust skeletons.Proc.of IEEE PacificGraphics 2000，Hong Kong，China.2000：281-289.”和“Bitter I，Sato M，Bender M，et al.CEASAR：A Smooth，Accurate and robust centerline extraction algorithm.Proc.of IEEEVisualization 2000，Salt Lake City，USA.2000：45-52.”。

上述现有的方法主要存在以下三个缺点：

(1)边界距离场变换是针对二值数据而言，因此需要预先提取血管组织。目前的血管提取方法，如匹配滤波器，主动轮廓模型以及水平集等方法虽然能够很好的将血管提取出来，但速度慢是一个很大的应用瓶颈；

(2)以二值数据为目标的边界距离场计算，对于不同位置的血管组织，其初始的边界距离值是相同的，不利于快速提取到位于血管中心的路径点；

(3)根据边界距离场提取的中心路径无法保证中心路径上的点处于血管管腔横截面的重心处，会导致血管的虚拟漫游过程与人类的观察习惯产生偏差。

发明内容

本发明的目的在于克服上述已有技术的不足，提供一种基于边界距离场的血管中心路径高精度快速提取方法，以实现对血管的快速提取和血管中心路径的高精度提取，使血管虚拟漫游过程更符合人类的观察习惯。

实现本发明目的的技术方案，包括如下步骤：

A.对输入的CTA血管造影体数据进行八叉树分解，得到每一分块的信息；

B.根据每个分块的最大值和最小值与双阈值之间的关系，遍历存在血管体素的叶子节点提取血管组织；

C.根据提取的血管组织中所有点的坐标，计算提取后的血管组织边界距离场；

D.在边界距离场的基础上建立最大代价生成树；

E.提取最大代价生成树的主干作为血管组织的初始中心路径；

F.利用重心法修正提取的血管组织初始中心路径，使路径上的点都处于血管管腔横截面的重心位置，得到一条符合人类观察习惯的最终血管中心路径。

上述步骤B所述的根据每个分块的最大值和最小值与双阈值之间的关系，遍历存在血管体素的叶子节点提取血管组织，按如下步骤进行：

B1.从深度为0的八叉树节点开始，判断其最大值和最小值与双阈值之间的关系，如果有至少一个阈值在最大值和最小值范围内，则该节点内存在血管体素，否则忽略该节点；

B2.判断存在血管体素的节点的深度值与设定的八叉树深度值之间的关系，如果节点的深度值与所设定的八叉树的深度值相等，则该节点为叶子节点，并遍历该叶子节点提取血管组织，否则判断该节点的子节点的最大值和最小值与双阈值之间的关系；

B3.如此循环，直到将所有包含血管体素的叶子节点遍历完毕为止，完成血管组织的提取。

上述步骤C所述的计算提取后的血管组织边界距离场，按如下步骤计算：

C1.根据提取的血管组织中所有点的坐标，利用原CTA血管造影体数据中血管体素的梯度值倒数和拉普拉斯变换值之和作为血管体素的边界距离初始值：

f_start(p)＝1/f_G+f_L

其中，p为提取的血管组织中任意一点的坐标，f_G为原CTA血管数据中点p的梯度值，f_L为原CTA血管数据中点p的拉普拉斯变换值；

C2.根据边界距离初始值用如下公式计算点p的边界距离值f_final(p)：

f_final(p)＝min(min(f_start(p₁)+10)，min(f_start(p₂)+14)，min(f_start(p₃)+17))其中，f_start(p₁)为点p的最近邻邻域点p₁的边界距离初始值，f_start(p₂)为点p的面对角邻域点p₂的边界距离初始值，f_start(p₃)为点p的体对角邻域点p₃的边界距离初始值；

C3.利用八叉树遍历所有包含血管体素的叶子节点，完成提取后的血管组织边界距离场计算。

上述步骤F所述的利用重心法对提取的血管组织初始中心路径进行修正，按如下步骤进行：

F1.对于初始中心路径上的每一点，取与该点法向量垂直，并且与该点相距为10个体素的点组成的圆形二维图像Ω；

F2.利用如下重心法公式求得Ω的重心C：

$C=[x_c,y_c]=(\frac{\mathrm{\Σxw}(x,y)}{\mathrm{\Σw}(x,y)},\frac{\mathrm{\Σyw}(x,y)}{\mathrm{\Σw}(x,y)})$

其中，w(x，y)＝α(f_2D(x，y)-m)，α＝-1，m＝max(f_2D(x，y))，f_2D(x，y)为Ω中坐标为(x，y)处的像素值；

F3.连接所有重心点作为最终的血管中心路径。

本发明由于采用了八叉树对CTA血管造影体数据进行分解，并利用原CTA血管数据中的相关信息计算血管组织的边界距离场，以及利用重心法修正血管组织初始中心路径，因此与现有技术相比具有如下优点：

(1)在血管提取和边界距离场计算过程中只遍历包含血管体素的节点，节省了时间；

(2)利用原CTA血管数据中的相关信息计算边界距离场，使处于不同位置的体素边界距离值相差更为明显，提高了中心路径提取的精度；

(3)利用重心法修正后的中心路径更符合人类的观察习惯，有利于血管疾病的诊断。

附图说明

图1为本发明的流程图；

图2为本发明对CTA血管造影体数据进行八叉树分解的示意图；

图3是本发明中三类邻域距离示意图；

图4为本发明中使用的两个CTA血管造影体数据三维重建结果图；

图5为本发明中对头部动脉瘤数据用不同方法提取的中心路径与手动提取的中心路径之间的距离对比示意图；

图6为本发明中对脑血管数据用不同方法提取的中心路径与手动提取的中心路径之间的距离对比示意图；

图7为利用本发明提取的中心路径进行血管虚拟内窥镜的效果图。

具体实施方式

本发明的核心思想是利用八叉树对输入的CTA血管造影体数据进行分解，并利用八叉树中存储的信息加速血管组织的提取和边界距离场的计算；利用原CTA血管数据中的相关信息提高边界距离场计算的精确度；采用重心法对提取的血管初始中心路径进行修正，提高中心路径的精度。

参照图1，本发明具体步骤如下：

步骤1，对输入的CTA血管造影体数据进行八叉树分解，得到每一分块的信息。

(1a)按照如下方式定义八叉树的属性；

Struct Octree

char LevelX，LevelY，LevelZ；//节点X，Y，Z方向上的深度，

unsigned int X_min，X_max，Y_min，Y_max，Z_min，Z_max；//节点的坐标范围，

bool IsLeaf；//是否子节点，

int Min_value，Max_yalue；//节点内的最小最大体数据值，

Octree*Children[8]；//节点的八个子节点；

(1b)根据设定的八叉树深度初始化八叉树，逐层为当前八叉树节点分配八个子节点；八叉树中最上层的节点表示整个体数据场，深度为0，如图2(a)所示，它被剖分成编号为0～7的八个深度为1的子节点，如图2(b)所示；当深度为2时，每个深度为1的节点被剖分为8个字节点，如图2(c)所示；以此类推，直到八叉树节点的成员变量LevelX，LevelY和LevelZ与预先设定的深度值相等时，剖分结束。此时八叉树节点为叶子节点，IsLeaf为True；剖分过程中同时记录每一分块的信息；体数据八叉树分解的树型表示如图2(d)所示。

步骤2，根据每个分块的最大值和最小值与双阈值之间的关系，遍历存在血管体素的叶子节点提取血管组织。

(2a)从深度为0的八叉树节点开始，判断其最大值和最小值与双阈值之间的关系，如果有至少一个阈值在最大值和最小值范围内，则该节点内存在血管体素，否则忽略该节点；

(2b)判断存在血管体素的节点的深度值与设定的八叉树深度值之间的关系，如果节点的深度值与所设定的八叉树的深度值相等，则该节点为叶子节点，并遍历该叶子节点提取血管组织，否则判断该节点的子节点的最大值和最小值与双阈值之间的关系；

(2c)如此循环，直到将所有包含血管体素的叶子节点遍历完毕为止，完成血管组织的提取。

步骤3，计算提取后的血管组织边界距离场。

(3a)根据提取的血管组织中所有点的坐标，利用原CTA血管造影体数据中血管体素的梯度值倒数和拉普拉斯变换值之和作为血管体素的边界距离初始值：

f_start(p)＝1/f_G+f_L

其中，p为提取的血管组织中任意一点的坐标，f_G为原CTA血管数据中点p的梯度值，f_L为原CTA血管数据中点p的拉普拉斯变换值；

(3b)根据边界距离初始值用如下公式计算点p的边界距离值f_final(p)：

f_final(p)＝min(min(f_start(p₁)+10)，min(f_start(p₂)+14)，min(f_start(p₃)+17))其中，f_start(p₁)为点p的最近邻邻域点p₁的边界距离初始值，f_start(p₂)为点p的面对角邻域点p₂的边界距离初始值，f_start(p₃)为点p的体对角邻域点p₃的边界距离初始值；三类邻域距离如图3所示；

(3c)利用八叉树遍历所有包含血管体素的叶子节点，完成血管组织的边界距离场计算。

步骤4，构建提取后的血管组织最大代价生成树。

借鉴Dijkstra动态规划思想，首先将提取的所有血管体素看成树的节点，除根节点S外，每个节点都指向一个父节点，从而形成一棵连接所有血管体素的有向树。在树的生成过程中利用边界距离值，定义A指向B的权值为W_AB＝DFB(B)，即如果B为A的父节点，则A连通到B的代价为DFB(B)。最大代价生成树的建立过程如下：

(4a)从源点S开始，将它所有26-邻域点均指向S，并以S为父节点，然后把它们的边界距离值放入一个排序队列q中；

(4b)从排序队列q中取出边界距离值最大的点M作为当前处理的点，让M的26-邻域中未被处理过的点均指向M，并以M为父节点，则这些连接边的权值都是DFB(M)，再把这些点的边界距离值放入排序队列q中；

(4c)如此循环，直到将所有提取的血管体素都处理完为止。

在建立最大代价生成树过程中，采用快速队列排序方法以提高效率。首先建立一个动态查找表，表的索引范围为0～255；接着在(4a)中把每个边界距离值归到0～255，再把该点的体素序号放入对应的查找表中；在(4b)中始终从索引值最大的不为空的查找表中提取体素序号，直到所有的查找表都为空。

步骤5，提取最大代价生成树的主干作为血管初始中心路径。

最大代价生成树主干的提取一般分为两步：

(5a)如果用户没有指定结束点，则默认源距离值最大的点为路径终点E；

(5b)从终点E开始，根据节点的指向关系直接连通到源点S作为血管初始中心路径。

步骤6，利用重心法对血管初始中心路径进行修正，使路径上的点都处于血管管腔横截面的重心位置，得到一条符合人类观察习惯的最终血管中心路径。

(6a)对于初始中心路径上的每一点，取与该点法向量垂直，并且与该点相距为10个体素的点组成的圆形二维图像Ω；

(6b)利用如下重心法公式求得Ω的重心C：

$C=[x_c,y_c]=(\frac{\mathrm{\Σxw}(x,y)}{\mathrm{\Σw}(x,y)},\frac{\mathrm{\Σyw}(x,y)}{\mathrm{\Σw}(x,y)})$

其中，w(x，y)＝α(f_2D(x，y)-m)，α＝-1，m＝max(f_2D(x，y))，f_2D(x，y)为Ω中坐标为(x，y)处的像素值；

(6c)连接所有重心点作为最终的血管中心路径。

本发明的优点可通过以下实验进一步说明：

1.实验条件

本发明的实验是在德国图宾根大学的两组CTA血管造影体数据上进行的。一组是头部动脉瘤数据，分辨率为512体素×512体素×512体素，实际大小为100×100×100mm³，三个方向上的相邻体素间隔均为0.1953mm，数据量为128MBytes。另一组是脑血管数据，分辨率为416体素×512体素×112体素，实际大小为171mm×210.5mm×46mm，三个方向上的相邻体素间隔均为0.412mm，数据量为22.7MBytes。头部动脉瘤数据的三维重建结果如图4(a)所示，脑血管数据的三维重建结果如图4(b)所示。

实验中提取了四种路径，一是多次手动提取的路径求平均作为一条最优路径，并将该路径作为与其他三条由方法所提取的路径进行比较的基准，称之为基准路径；二是基于边界距离场方法提取的路径，但未使用梯度值及拉普拉斯变换值作为边界距离变换初始值，且未使用重心法修正，称之为二类路径；三是基于边界距离场方法提取的路径，使用梯度值及拉普拉斯变换值作为边界距离变换初始值，但未使用重心法修正，称之为三类路径；四是本发明基于边界距离场方法提取的路径，既使用了梯度值及拉普拉斯变换值作为边界距离变换初始值，也使用了重心法修正，称之为四类路径。实验在这两组体数据中比较用八叉树分解的方法进行血管提取所用的时间与用全局遍历方法进行血管提取所用的时间；比较用八叉树分解的方法进行边界距离场计算所用的时间与用全局遍历方法进行边界距离场计算所用的时间；比较二类路径、三类路径、四类路径与基准路径的距离。

2.实验结果

实验中两种方法在两组数据中所用时间对比如表1和表2所示。两组数据中所提取的二类路径，三类路径，四类路径与基准路径之间的距离对比如图5和图6所示。

表1头部动脉瘤数据八叉树方法与全局遍历方法血管提取及边界距离场计算用时比较

表2脑血管数据八叉树方法与全局遍历方法血管提取及边界距离场计算用时比较

表1中，提取阈值在三种不同情况下，八叉树大小为16体素×16体素×16体素时，血管提取用时与边界距离场计算用时之和分别为1.808秒，1.719秒，1.625秒，比同阈值下全局遍历所用时间之和4.227秒，4.242秒，4.204秒分别要少57.2％，59.4％，61.3％。

表2中，提取阈值在三种不同情况下，八叉树大小为64体素×64体素×64体素时，血管提取用时与边界距离场计算用时之和分别为0.037秒，0.036秒，0.030秒，比同阈值下全局遍历所用时间之和3.969秒，3.938秒，3.906秒分别要少99％，99％，99.2％。

从表1和表2中可以看出，加入八叉树分解后，血管的提取速度和边界距离场的计算速度得到了大幅度的提高，节省了时间消耗。

实验中共提取了两段血管的中心路径，图5中的路径共有85个点构成，图6中的路径共有68个点构成。两图中，黑色线为二类路径与基准路径的距离，蓝色线为三类路径与基准路径的距离，红色线为四类路径与基准路径的距离；图的横坐标为路径上点的序号，纵坐标为该点与基准路径上对应点的距离值，单位为体素；图中*，+，○均表示路径上的点。图5中，黑色线的均值为3.69，即二类路径与基准路径之间的距离均值为3.69个体素；蓝色线的均值为2.85，即三类路径与基准路径之间的距离均值为2.85个体素；红色线的均值为2.21，即四类路径与基准路径之间的距离均值为2.21个体素。图6中，黑色线的均值为3.46，即二类路径与基准路径之间的距离均值为3.46个体素；蓝色线的均值为2.36，即三类路径与基准路径之间的距离均值为2.36个体素；红色线的均值为1.68，即四类路径与基准路径之间的距离均值为1.68个体素。直观上，在两图中，黑色线在蓝色线之上，而蓝色线又在红色线之上，即四类路径最接近基准路径。因此，利用梯度值倒数及拉普拉斯变换值作为边界距离变换初始值，且使用重心法修正从最大代价生成树提取的路径可使最终路径最接近血管空腔的真实中心路径。利用四类路径对头部动脉瘤数据进行血管虚拟内窥镜的效果如图7(a)所示，利用四类路径对脑血管数据进行血管虚拟内窥镜的效果如图7(b)所示。

上述实验结果中，本发明用八叉树方法进行血管提取和边界距离场计算，与目前存在的全局遍历方法进行血管提取和边界距离场计算相比，用时要少；利用原CTA血管造影体数据的相关信息得到的距离场，比目前存在的利用二值血管数据得到的距离场精度要高；利用重心法修正后的路径比直接提取基于边界距离场的最大代价生成树主干得到的中心路径更贴近人类的观察习惯。

Claims (3)

1.一种基于边界距离场的血管中心路径提取方法，包括如下过程：

A.对输入的CTA血管造影体数据进行八叉树分解，得到每一分块的信息；

B.根据每个分块的最大值和最小值与双阈值之间的关系，遍历存在血管体素的叶子节点提取血管组织；

C.根据提取的血管组织中所有点的坐标，计算提取后的血管组织边界距离场：

C1)根据提取的血管组织中所有点的坐标，利用原CTA血管造影体数据中血管体素的梯度值倒数和拉普拉斯变换值之和作为血管体素的边界距离初始值：

f_start(p)＝1/f_G+f_L

其中，p为提取的血管组织中任意一点的坐标，f_G为原CTA血管数据中点p的梯度值，f_L为原CTA血管数据中点p的拉普拉斯变换值；

C2)根据边界距离初始值用如下公式计算点p的边界距离值f_final(p)：

f_final(p)＝min(min(f_start(p₁)+10)，min(f_start(p₂)+14)，min(f_start(p₃)+17))其中，f_start(p₁)为点p的最近邻邻域点p₁的边界距离初始值，f_start(p₂)为点p的面对角邻域点p₂的边界距离初始值，f_start(p₃)为点p的体对角邻域点p₃的边界距离初始值；

C3)利用八叉树遍历所有包含血管体素的叶子节点，完成提取后的血管组织边界距离场计算：

D.在边界距离场的基础上建立最大代价生成树；

E.提取最大代价生成树的主干作为血管组织的初始中心路径；

F.利用如下重心法修正提取的血管组织初始中心路径，使路径上的点都处于血管管腔横截面的重心位置，得到一条符合人类观察习惯的最终血管中心路径：

_F1)对于初始中心路径上的每一点，取与该点法向量垂直，并且与该点相距为10个体素的点组成的圆形二维图像Ω；

F2)利用如下重心法公式求得Ω的重心C：

$C=[x_c,y_c]=(\frac{\mathrm{\Σxw}(x,y)}{\mathrm{\Σw}(x,y)},\frac{\mathrm{\Σyw}(x,y)}{\mathrm{\Σw}(x,y)})$

其中，w(x，y)＝α(f_2D(x，y)-m)，α＝-1，m＝max(f_2D(x，y))，f_2D(x，y)为Ω中坐标为(x，y)处的像素值；

F3)连接所有重心点作为最终的血管中心路径。

2.根据权利要求1所述的方法，其中过程A所述的对输入的CTA血管造影体数据进行八叉树分解，按如下步骤进行：

A1.按照如下方式定义八叉树的属性；

Struct Octree

char LevelX，LevelY，LevelZ；//节点X，Y，Z方向上的深度，

unsigned int X_min，X_max，Y_min，Y_max，Z_min，Z_max；//节点的坐标范围，

bool  IsLeaf；//是否子节点，

int  Min_value，Max_value；//节点内的最小最大体数据值，

Octree*Children[8]；//节点的八个子节点；

A2.根据所设定的八叉树的深度值LevelX，LevelY，LevelZ对输入的CTA血管造影体数据进行逐层分解，同时利用八叉树记录每一分块的信息。

3.根据权利要求1所述的方法，其中过程B所述的根据每个分块的最大值和最小值与双阈值之间的关系，遍历存在血管体素的叶子节点提取血管组织，按如下步骤进行：

B1.从深度为0的八叉树节点开始，判断其最大值和最小值与双阈值之间的关系，如果有至少一个阈值在最大值和最小值范围内，则该节点内存在血管体素，否则忽略该节点；

B2.判断存在血管体素的节点的深度值与设定的八叉树深度值之间的关系，如果节点的深度值与所设定的八叉树的深度值相等，则该节点为叶子节点，并遍历该叶子节点提取血管组织，否则判断该节点的子节点的最大值和最小值与双阈值之间的关系；

B3.如此循环，直到将所有包含血管体素的叶子节点遍历完毕为止，完成血管组织的提取。

Images