1 引言

目前,多层螺旋电子计算机断层扫描(CT)广泛应用于心脏疾病的分析和手术导航过程。先进的医学影像设备可以提供高分辨率的心脏及周围心血管组织的序列影像^[1-4]。心脏CT影像可辅助评价心脏主动脉及冠状静脉疾病的异常情况,先进的三维显示技术能辅助心血管手术导航等。以上医疗辅助手段均建立在心脏及心血管图像分割结果的基础上,医学领域迫切需要快速且精确的心血管分割方法。尽管在图像分割领域的很多研究人员已提出了多种心血管图像分割方法^[5-9],如基于区域增长法、阈值分割、体素生长算法、深度学习分割、混合算法等,但兼顾分割精度和分割速度的分割算法的提出仍然是一项艰巨的挑战^[10]。基于区域增长的各种改进算法在血管分割领域得到了普遍认可,但这些算法在血管边缘溢出判断及局部感兴趣血管的提取方面仍存在一定局限性。

O'Brien等^[7]使用区域增长法完成了心脏冠状动脉血管的分割,采用低通滤波对心脏CT图像进行预处理,利用冒泡测试算法规划分割后血管结构的中轴线;Schmitt等^[8]以阈值方法为依据,采用区域生长方法获得了血管树的结构;李雪丽等^[9]采用Live Wire与区域生长法相结合的算法对算法进行简化和优化,获得了比较满意的血管图像分割结果。由于医学图像数据量大,以往分割算法采用二维技术对肺部序列CT切片逐张分割,导致计算复杂,分割速度较慢,不适于临床应用。

文献[ 11]提出的自适应体素生长算法完成了基于统一计算设备架构(CUDA)加速的左右肺叶三维分割,明显提高了分割效率,同时文中提到了解决体素泄漏问题的方法是对肺部图像进行预处理,从而消除肺叶的两个泄露通道,解除泄露的预处理方法是以预分割方式,将两片肺叶连通区域彼此分开,或用区域增长算法分割气管,以避免肺叶体素生长过程中肺叶连接处的体素泄露。

虽然基于体素的全自动分割方法在三维分割方面避免了人工交互过程,但该方法不能直接用于心脏图像的分割,原因在于心脏解剖结构复杂,在心肌周围分布着各类动脉和静脉^[5],这些血管与心肌的连通区域用在无交互的全自动分割算法中会导致体素大量泄漏;另外,如果医生无法观察局部任意感兴趣的血管组织,就会导致该算法分割对象不灵活,具有一定的局限性。到目前为止,还没有一种适合于心血管图像的全自动三维分割算法。

体素生长算法是一种可以快速完成三维体分割的技术,该算法被用于许多三维医学组织(如肺部、肝脏、肠道消化系统等)图像的分割。防止体素泄露是应用体素生长算法首先要解决的问题,要求三维模型的连通域尽量少,由于肺叶结构存在的连通区域主要集中在肺叶和气管连接处以及左右肺叶的连通区域^[11],少量的泄露区域在预处理时将其填补便能实现体素生长,但心血管组织在结构上相对肺叶结构更加复杂,心肌的连通区域较多,很难填补所有泄露区,且较难寻找组织边界,不易进行预处理,导致提取局部血管时会出现体素泄漏的情况。综上所述,体素生长算法不适合直接用于心血管的分割^[12],因此,本文提出了一种适合于心脏血管分割的半自动分割方法。

2 体素生长算法的改进

为了解决心血管的分割问题,本研究采用改进的半自动体素生长算法。对算法的改进主要从两个方面进行:1) 对心血管进行预处理,手动设定局部血管的种子点范围,提取局部感兴趣血管的局部AABB包围盒(axis aligned bounding box)对应的三维模型,即感兴趣体积,为心脏三维模型设定空间边界;2) 对局部血管的像素进行体素生长,生长规则为对比相邻体素的灰度相似性,相似的体素完成生长过程,生长终止条件为找到与血管壁空间曲面相交的体素。

以上改进措施具有以下优势:1) 划定局部感兴趣区域血管的范围时可以人为设定体素生长边界,不会出现泄露的情况;2) 提取局部三维模型,代替完整的三维模型,大大减小了计算量,显著提升了算法的分割效率;3) 用血管拓扑形态公式约束体素生长,可以减少因种子点体素选错而导致的越界心肌范围内体素的盲目错误扩充;4) 医生可以在心脏三维模型上划定感兴趣区域种子点,提取任意感兴趣部分的血管,算法灵活,分割对象不受局限。改进后的算法能显著提高局部血管的分割效率,能够直接处理三维模型,体素生长速度快;边界检查采用体素是否与等值面相交的方式约束能提高分割的准确性,兼顾分割速度和准确性,可以满足临床需求。

2.1 体素生长规则

线性八叉树是对三维模型体素进行邻域拓扑的重要手段。体素拓扑遵循26邻域的方式,每次寻找与体素相邻的26个相邻体素,并将每个体素重新分割成体积相等的8个新体素。三维模型被离散成体素后,通过人工交互的方式选择模型中某个或某些体素点作为体素生长的种子点体素,由这些体素寻找周围26邻域中灰度值相似的体素,如果相似则合并入种子点集合中,作为新的种子点继续与周围26邻域的体素进行对比,将相似体素空间不断扩大,最终按照终止条件判断出边缘体素,完成体素空间的生长。

寻找周围的相似体素通过对比相邻两体素的灰度相似度来完成。相似度是通过相似性函数实现的,可通过比较体素灰度均方误差的方式实现两个体素相似度的比较。均方误差公式为

eMSE=1N∑i=1NAi-T(Bi)2,(1)

式中:N为三维空间中的体素数目;A_i为浮动体素空间;T(B_i)为配准后的体素空间;A_i和B_i为相邻两体素在三维空间中相同位置A_i(a_i,b_i,c_i)和B_i(a_i,b_i,c_i)处对应的空间灰度值,对应点的灰度值主要通过读取两体素对应点的灰度直方图来获得。基于均方误差的相似性测度的优点是比较法则相对简单和计算速度快,能够在对应三维空间中进行快速检索,但是要求两体素的模态和对应离散度必须保持一致,以确保两个体素能够配准成功,对体素进行配准预处理可减小体素相似度的误差。表1列出了手动选取起始种子点的26个拓展体素,其中体素相似度排在前5位的分别为:R_a、R_b、R_c、R_d、R_e。从表1可以看出:经过配准预处理后,体素相似度明显变大,尤其是均方误差的相似度明显提高。从表1中3种相似度测度的对比结果可以看出:配准预处理后,均方误差的相似体素查找精确度为17/26,介于相关系数18/26和Matters互信息15/26之间;在运行时间方面,从表3的体素生长时间3.78 s可以看出,选择均方误差完成相似度对比可基本满足精度和时间的要求。

2.2 终止条件

终止条件即体素生长停止的条件,体素生长算法的终止条件是评价分割算法好坏的关键。如何找到血管边缘体素集合是找到终止条件的前提。

2.2.1 等值面的函数表达式

等值面是灰度空间内具有相同灰度值的点集,等值面的三维曲面函数可表示为

{(x,y,z)|f(x,y,z)=a}x=x(u,v),(u,v)∈Δy=y(u,v),(u,v)∈Δz=z(u,v),(u,v)∈Δ,(2)

式中:a为血管壁的灰度值常量;Δ为三维空间的子空间。由于血管壁可用空间曲面表示,可以用(2)式所示的三元参数函数表示等值面曲面,该曲面和体素所在立方体相交的曲线为双曲线^[13],由于计算空间双曲线与立方体体素的每个棱的交点相当耗时,加之体素所在立方体很小,故而在每个面上用直线相似拟合双曲线所求出的交点误差可以被忽略。立方体的边界线用函数可表示为

f(x,y,z)=a0+a1x+a2y+a3z+a4xy+a5yz+a6zx+a7xyz,(3)

式中:a₀、a₁、a₂、a₃、a₄、a₅、a₆、a₇为代入体素立方体顶点坐标后体素立方体各棱函数的系数。

若将等值曲面的灰度值设为a',则通过(4)式可以计算出哪些体素能够与血管壁曲面相交,从而判断出血管的边界体素。(4)式的表达式为

f(x,y,z)=a'f(x,y,z)=a0+a1x+a2y+a3z+a4xy+a5yz+a6zx+a7xyz。(4)

2.2.2 血管边缘体素的查找

可以利用(4)式中解的个数来判断等值面和体素是否相交,当(4)式有解时证明等值面和体素是相交关系,并可以通过解的坐标找到对应体素集合。

三维模型的体素空间中与血管曲面相交的体素是按照一定的逻辑性聚集在一起构成血管壁的体素集合,将这些与血管曲面相交的体素称为相交体素,而没有与血管相交的体素称为非相交体素,相交体素的查找按照以下逻辑进行:

1) 建立一个相似度查找表(SLT),存储每一个相似体素的体素灰度值均方误差相似度。

2) 建立一个标志位(FLAG),用于标志遍历八叉树中哪些像素节点是相交体素,若体素为相交体素,则FLAG=1;若体素为非相交体素,则FLAG=0。

3) 建立两个队列Queue A和Queue B,Queue A存储相交体素的物理位置,Queue B存储非相交体素的物理位置。

4) 线性八叉树通过查找26邻域的方式遍历三维模型的体素,并将体素标记为两个集合,即FLAG=1的相交体素集与FLAG=0的非相交体素集。

5) 在FLAG=1的相交体素集合中通过相似度查找表找到相似体素的相似度,通过相似度对比找到相似的体素,再通过判断相似体素是否在八叉数的相邻子节点中确定这些体素的物理位置是否也是相邻的,体素的物理位置可以在Queue A中找到。FLAG=0的非相交体素集合也通过该方式查找得到,非相交体素的物理位置能够在Queue B中找到。寻找相交体素的算法流程图如图1所示。

图 1. 根据相似性查找表查找边缘体素的流程图

Fig. 1. Flow chart of finding the edge of the voxel according to similarity table

下载图片 查看所有图片

2.2.3 拓扑形态约束

体素在寻找相邻相似体素过程中,在处理非血管区域(如心肌区域)时,由于连通区域过大,导致种子点骤然增加。由于体素生长算法是一种嵌套算法,种子点的突然增加会降低算法速度,影响实时交互性。为了保证体素生长算法在处理黏连组织中依然适用,该算法在传统体素生长算法的基础上增加了拓扑形态约束,以限制种子点空间的立体形态,从而避免出现大体积连通区域。体素在血管内部拓扑的示意图如图2所示。

图 2. 血管中体素的拓扑约束示意图

Fig. 2. Topology constraint diagram of voxels in blood vessels

下载图片 查看所有图片

为种子点设定体素生长范围,体素在每次扩充时向周围增加体素量,该过程中要设置拓扑结构的形态约束,体素拓扑空间变成非血管形状时,证明该体素已经不在血管内部,此时已出现体素泄露的情况,为保证体素仍在血管中扩充,可以设定1个阈值和1个公式来约束体素的生长空间仍然在血管内部,拓扑形态约束为

πrmin2hmin≤b≤πrmax2hmax,(5)

式中:b为血管内体素体积的增加量;r_min为体素外接圆柱的底面半径;r_max为到生长边界体素的最小距离;h_min为体素最小外接圆柱的高;h_max为到生长边界体素的最大距离。

3 感兴趣血管的提取

3.1 血管AABB包围盒提取

图3所示为程序通过约束种子点生长空间坐标获得的局部三维血管簇,医生能够提取心脏三维模型上感兴趣的局部AABB包围感兴趣体积,这种空间的限制能够约束局部血管的边界,从而防止体素生长溢出问题。该粗略的提取结果可能同时包括除感兴趣血管以外的其他心血管组织,不便于医生观察。

3.2 AABB包围盒中感兴趣血管的提取

为了能够提取局部三维血管簇中的某支血管,通过设置种子点,用改进的体素生长算法判断血管内部体素的生长和边缘体素,图4所示为在AABB包围盒上通过鼠标交互方式划定局部感兴趣血管的种子点范围提取的血管模型。由于心脏CT图像生长过程容易出现边缘模糊、毛刺和片状粗糙现象,为了便于医生观察,用线性插值作重采样来提高其清晰度。从图4中的血管图像可以观察到血管三维显示的效果逼真且平滑。

图 3. (a)胸腔三维重建及(b)局部血管簇的提取结果

Fig. 3. (a) Three-dimensional reconstruction of thoracic cavity and (b) extracted local blood vessel cluster

下载图片 查看所有图片

图 4. 心脏局部AABB包围盒中局部血管的提取效果

Fig. 4. Extraction of local blood vessels from cardiac local AABB bounding box

下载图片 查看所有图片

图 5. 算法改进前后的分割效果。(a)改进前;(b)改进后

Fig. 5. Segmentation results before and after the algorithm improving. (a) Before improving; (b) after improving

下载图片 查看所有图片

3.3 性能分析

图5所示为心脏局部血管提取的对比图。将基于传统体素生长算法和改进前的体素生长算法进行性能对比,将传统体素生长算法在MATLAB上调试并实现分割效果,选用分辨率为256 pixel×256 pixel×256 pixel的三维分割精度作为评价标准。从图5中可以明显看出,传统体素生长算法的三维分割结果出现了边缘模糊和过分割的情况,改进的基于体素生长的三维分割算法的准确性明显优于传统体素生长算法,血管立体效果清晰且管壁光滑,无过分割情况。表明本研究算法既可以提高三维血管的显示清晰度,又能够满足血管三维提取的算法实时性,在一定程度上体现了本研究算法的优势。

将改进的体素生长算法与传统体素生长算法进行算法性能评价^[13-14],结果如表2所示,用4组参数对算法性能进行对比,可知:1) 敏感性和特异性分别决定了分割结果过分割和边界模糊的程度,从表2中可以观察到传统算法的敏感度较低,且特异性与改进算法相比没有优势,本研究算法在边缘检测和过分割方面优于传统算法;2) 本研究算法的准确性优于传统体素生长算法,最高准确度达到0.9483,传统算法的最高精确度只有0.8997;3) Kappa系数体现了与标准结果的相似度,本研究算法的分割标准结果由包头医学院第二附属医院的两名专业医师手工分割完成,结果表明本研究算法更加接近标准分割。由表2可以清楚地看到,本研究算法各参数均优于传统体素生长算法,在一定程度上体现了本研究算法的优势。

表3所示为心脏三维模型分辨率均为256 pixel×256 pixel×256 pixel的前提下,全自动体素生长算法和半自动体素生长算法改进前后体素生长时间的评价对比表.从表3中可以看出:全自动体素生长算法只针对模型整体,内部三维模型体素数的生长时间为103.89 s;半自动体素生长算法针对局部模型,取相同局部AABB包围盒提取的心脏三维模型,在经过心脏模型表面和内部体素化后,初始体素数均为3152个,表面体素数均为345个;在初始数据分辨率相同的前提下,未经改进的体素生长算法内部体素生长时间和表面体素生长时间分别为5.59 s和1.46 s;改进的体素生长算法内部体素生长时间和表面体素生长时间分别为3.78 s和1.12 s 。从体素的生长时间可以看出,改进体素生长算法的体素生长时间明显缩短。

4 结论

提出了半自动交互式体素生长血管分割算法,完成了局部心脏感兴趣血管的三维分割。对于庞大的三维数据,与传统体素生长算法相比,所提算法无论在分割效果还是算法性能上都有所提高。与全自动体素生长算法相比,所提算法能够通过设置交互种子点从心脏三维模型中提取出AABB包围盒,包围盒约束了体素集合的空间边界,防止了无范围的体素扩充所带来的体素泄露问题,最终实现了任意局部心脏血管的体素生长。所提算法用体素生长法和寻找边界体素的方式准确分割出局部血管的三维模型,剔除其他不感兴趣以及不便于医生观察的血管组织。与自适应的全自动体素生长算法相比,所提算法虽然存在人工交互设置种子点的过程,但进行体素生长的对象是局部三维模型,而全自动体素生长算法是对整个三维模型进行计算,因此所提算法的体素生长时间明显缩短,实时性有所提升,保证了算法的实时性和精确性,可为医生诊断心脏主动脉及冠状静脉疾病提供可靠的可视化辅助依据。该技术经改进后还可以辅助仿真结肠镜和仿真支气管镜的检查,在医学物理与工程学领域也具有一定的研究意义。