CN102402796B

CN102402796B - 肝脏血管系统的三维结构化描述方法

Info

Publication number: CN102402796B
Application number: CN 201110330817
Authority: CN
Inventors: 房斌; 王鹏; 王翊; 皮净锐
Original assignee: Chongqing University
Current assignee: SHENZHEN YORKTAL DMIT CO Ltd
Priority date: 2011-10-26
Filing date: 2011-10-26
Publication date: 2013-09-18
Anticipated expiration: 2031-10-26
Also published as: CN102402796A

Abstract

本发明公开了一种肝脏血管系统的三维结构化描述方法，属于医学肝脏划分领域。该方法首先对肝脏图像和肝脏血管图像分别进行分割，再对肝脏血管三维图像进行结构分析，将肝脏血管骨架化处理，对骨架化后的体素点进行标记，并且将标记后的体素点存储至多叉树的数据结构中，并且根据体素点在多叉树中的主从关系区分主支血管和分支血管，最后根据肝脏血管系统中血管之间的关系将肝脏体划分成多个肝段。通过本发明，实现了对肝脏的精确划分，从而为肝脏的精准切除提供指导。

Description

肝脏血管系统的三维结构化描述方法

技术领域

本发明涉及肝脏血管系统，尤其是一种肝脏血管系统的三维结构化描述方法。

背景技术

近年来随着科技的发展，采用计算机图形学与数据可视化技术，利用计算机图像处理和分析方法，通过三维重建可以帮助医生将二维影像数据构造成立体的三维可视化模型，清楚显示模型内管道系统复杂的走行、分布及空间结构位置关系，同时实现影像数据的量化分析，提供管道的直径、长度及所支配实质脏器的区域等信息。在三维可视化模型的基础上，进行虚拟切割，建立计算机辅助手术规划平台，使手术可行性的术前评估更加精准，术前准备更加充分，手术方案的制定更加客观、合理，因此手术的创伤更小，出血更少，预后更佳。目前医学影像的三维可视化研究和数字化手术规划平台研发成为国际研究热点。

在肝胆胰外科领域，国内第三军医大学董家鸿教授从2000年以来率先进行了肝胆胰等脏器三维重建的探索；2003年以来，南方医科大学方驰华教授、福建医科大学刘景丰教授等也相继开展了这方面的研究。目前已初步建立了肝胆胰等脏器的三维可视化模型，并在肝脏模型上进行了虚拟切割探讨。但是在上述模型中，仍然存在不同程度信息丢失，肝内管道在肝段水平的显示不佳，同时缺乏数据的量化分析功能，因此严重制约了其临床应用价值。国际上，德国和日本开发的计算机辅助肝脏切除手术规划系统已经进入临床，初步研究提示这个系统有助于提高肝脏切除手术的精准性和安全性。在胆道和胰腺方面，尚未见数字化手术规划系统进入临床的报道。

虽然临床研究初步提示数字化手术规划系统有助于提高肝脏切除手术的精准性和安全性，但是由于数字化手术规划系统在肝脏切除手术中应用时间尚短，因此其在肝脏切除手术中的价值与地位仍然需要进一步检验，尤其是需要按照循证医学的原则进行系统研究。目前关于数字化肝切除手术规划的临床研究只是个案报道和病例回顾，因此本项目中设计了前瞻性的对照研究，这将有助于深入认识数字化手术规划平台的临床应用价值。

关于胆道和胰腺的数字化手术规划平台虽然尚未进入临床，但是鉴于理论上已经成为可能，尽快研发相关软件，使之早日进入临床就成为当务之急。

目前国内外的数字化手术规划平台建设主要是基于对脏器和病灶的解剖学精确评估。然而一个合理的手术规划，不仅取决于解剖学评估，还需要结合病理学、生理学、功能影像学等进行综合考量。“精准外科”理念已经得到国内外同行的广泛赞誉和认同，本项目在数字化手术规划平台构建中融入了预留肝脏功能数字化评估、以及预留结构和功能重建等内容，将“精准外科”理念与数字化手术规划系统融合，形成数字化精准外科手术规划系统，制定符合精准外科理念的个体化手术规划。

鉴于三维可视化模型和数字化手术规划平台对于提高手术精准性和安全性的巨大优势，其在外科临床的应用前景非常广阔。同时，三维可视化模型和数字化手术规划平台还是未来进一步开发虚拟手术系统的基础。利用虚拟现实技术，让外科医生在计算机建立的虚拟环境中进行解剖观察、手术设计和手术操作，获得实际手术中的手感，这将对提高外科整体手术水平具有革命性影响。

发明内容

本发明的目的是提供一种肝脏血管系统的三维结构化描述方法，在对肝脏血管三维图像进行结构分析时，根据血管供血区域划分肝段，实现了对肝脏的精确划分，指导肝脏的精准切除。

为了实现上述目的，本发明提供了一种肝脏血管系统的三维结构化描述方法，其特征在于包括以下步骤：

A、对肝脏体图像进行分割，获得肝脏三维图像的分割结果；

B、对肝脏血管图像进行分割、处理，获得合成的肝脏血管三维图像：

B01、采用基于区域增长的阈值分割方法对肝脏血管图像进行分割，提取肝脏血管的信息；

B02、采用形态学操作和基于空间的中值滤波操作去除噪音，获得粗分割图像集；

B03、采用基于空间连通域的血管树追踪方法对所述粗分割图像集进行三维血管树追踪，获得连通的三维血管树图像集；

B04、采用三维形态学膨胀腐蚀操作对所述连通的三维血管树图像集的表面进行平滑处理，获得合成的肝脏血管三维图像；

其中所述步骤B02由以下步骤组成：

B02-1、经阈值分割后的肝脏血管图像通过形态学操作消除噪音、去除缺陷，并且通过二维中值滤波处理，消除孤立的噪声点，从而获得图像集合R；

B02-2、所述图像集合R中以三张图像作为一个图像单元R_i，将所述图像集合R表示为R=（R₁，R₂，…，R_n），如果所述图像集合R中图像数量不满足3的倍数，则将剩余的图像通过二维中值滤波处理输出，构成的图像单元R_i则执行如下处理；

B02-3、初始化nCount=0,neiborNum=15，其中nCount用于表示像素点的邻域中像素值不为0的像素点的个数，neiborNum用于表示像素点的邻域中像素值不为0的像素点的个数临界值，i、n均为≥1的整数；

B02-4、将所述图像单元R_i中像素值不为0的像素点均投影至一张图像上，计算投影后像素点坐标为（x，y）的像素值大小，如果所述像素点的周围存在26邻域像素点则以该像素点为中心统计其周围26邻域像素点的像素值是否为0，所述26邻域像素点中每存在一个像素值不为0的像素点，nCount加1；

B02-5、如果该像素点的nCount>neiborNum，则将该像素点的像素值置为255，否则置为0；

B02-6、重复执行步骤B02-4～B02-5，直至完成图像集合R中所有图像单元的过滤，获得粗分割图像集；

C、对所述肝脏血管三维图像进行结构分析：

C01、采用三维骨架化方法细化处理所述肝脏血管三维图像；

C02、分别采用端点、曲线点和分叉点标记肝脏血管三维图像的体素点，从而获得标记后的体素点，其中所述端点只具有一组邻接点，所述曲线点具有两组邻接点，所述分叉点具有3组以上邻接点；

C03、采用多叉树的数据结构来存储所述标记后的体素点，一个体素点的坐标存储在多叉树的一个结点中，且将多叉树中的父节点和子节点采用直线连接，从而形成肝脏血管的抽象树结构；

C04、在所述肝脏血管的抽象树结构中将两个连接的分叉点视为一个连通域，采用连通域标记方法对各连通域进行标记，赋予唯一的数值，即将连通域中体素点均标记为相同的唯一的数值，并且采用迭代方法将与各连通域中体素点邻接的体素点均标记为相同的数值，从而区分肝脏血管系统的主支血管和分支血管，获得标记后的肝脏血管系统；

D、根据肝脏血管系统中的主支血管和分支血管标记以及医学需要，将肝脏体划分为多个肝段。

所述步骤A由以下步骤组成：

A01、对肝脏体图像进行预处理，获得原始图像：将所述肝脏体图像转换为Bmp形式，并且去除噪声；

A02、采用四叉树分割方法对所述原始图像进行初始分割，获得多个小区域；

A03、以所述小区域v_i作为Graph-Based图像分割方法中带权图G的节点，所述带权图G=（V，E），其中V=（v₁，v₂，…，v_n）表示所述小区域的集合，由n个小区域组成，E表示所述小区域之间的边集，由m条边组成，n、m均为大于1的整数，i为≥1的整数并且边集E中每条边对应于一个权值，该权值用于表示所述小区域之间的灰度相似度；

A04、根据边集E中各边对应权值的大小，将边集E进行排序E=(e₁,Λ,e_p,Λ,e_m)，其中1<p<m为整数；

A05、设定起始的分割结果为S⁰=（v₁，v₂，…，v_n），即将集合V中各小区域均作为起始分割的连通分量；

A06、如果节点v_i、v_j分别属于第p-1次分割结果S^p-1中的不同连通分量，则两节点之间的边为e_p=(v_i,v_j)，边e_p对应的权值为w(e_p)，分别求取两连通分量v_i、v_j的内部差异Int(v_i)、Int(v_j)，如果权值w(e_p)大于所述内部差异Int(v_i)、Int(v_j)中之一则不合并节点v_i、v_j，否则合并节点v_i、v_j，从而获得第p次分割结果S^p；

A07、循环执行步骤A06，直至完成m次分割，获得肝脏三维图像最终的分割结果。

所述步骤B03由以下步骤组成：

B03-1、将所述粗分割图像集中各粗分割图像的每一个不为零的体素点的坐标均保存至索引Index中，并且将所有像素点的像素值均置为1，其中索引Index_j对应于所述粗分割图像集中第j张图像，下标j为大于0的整数；

B03-2、遍历所述索引Index_j，将所述索引Index_j中所有像素点的坐标均推入栈Stack中；

B03-3、弹出所述栈Stack中的各体素点，搜索各体素点的26邻接点集，如果在像素点的26邻域接点集中存在像素值为1的像素点，则将该像素点的像素值标记为L_j=L_j-1+1，并且将该像素点的坐标推入栈Stack中，其中L_j-1表示遍历索引Index_j-1时，栈Stack中各像素点的26邻接点集中存在像素值为1的像素点时，该像素点对应的像素值；

B03-4、重复执行步骤B03-3，直至所述栈Stack为空；

B03-5、重复执行步骤B03-2～B03-4，直至索引Index遍历完成；

B03-6、根据像素值的大小统计像素点的个数，将像素点个数最多的连通域作为所述连通的三维血管树图像集。

所述步骤C01中采用12次子迭代方法细化所述肝脏血管三维图像。

所述步骤C02由以下步骤组成：

C02-1、由counter来表示所述肝脏血管三维图像中体素点的邻接点组数，初始化counter=0；

C02-2、当体素点p的6-邻接点集N₆(p)中除体素点p外还存在黑点x时，判断体素点p的18-邻接点集N₁₈(p)与6-邻接点集N₆(p)的补集或者体素点p的26-邻接点集N₂₆(p)与18-邻接点集N₁₈(p)的补集中是否具有在黑点x的

范围内的体素点，如果判定两个补集仅其中一个具有在黑点x的

范围内的体素点，则counter加1；

当体素点p的18-邻接点集N₁₈(p)与6-邻接点集N₆(p)的补集中存在黑点x时，判断体素点p的26-邻接点集N₂₆(p)与18-邻接点集N₁₈(p)的补集中是否具有与黑点x相邻的体素点，如果判定不具有与黑点x相邻的体素点则counter加1；

当体素点p的26-邻接点集N₂₆(p)与18-邻接点集N₁₈(p)的补集中具有黑点，则将counter增加的数值等于体素点p的26-邻接点集N₂₆(p)与18-邻接点集N₁₈(p)的补集中黑点的个数；

C02-3、根据体素点p的邻接点组数counter来标记体素点p：如果counter=1则将体素点p标记为端点，如果counter=2则将体素点p标记为曲线点，如果counter>2则将体素点p标记为分叉点。

所述步骤D由以下步骤组成：

D01、设定所述标记后的肝脏血管系统为肝脏血管树数据集Vas，肝脏体掩模数据集为Mask，并且设定所述肝脏血管树数据集Vas中像素值大于0的像素点的个数为vasSum，所述肝脏体掩模数据集Mask中像素值大于0的像素点的个数为maskSum，其中所述肝脏体掩模数据集Mask即为肝脏体的所有数据，所述肝脏血管树数据集Vas为所述肝脏体掩模数据集Mask的子集；

D02、遍历所述肝脏体掩模数据集Mask与所述肝脏血管树数据集Vas的补集中的各点P_i，并且设定坐标为（x_i，y_i，z_i）的各点P_i均存在一个窗体空间Box（P_i），Box（P_i）定义为（x_i-scale：x_i+scale，y_i-scale：y_i+scale，z_i-scale：z_i+scale），即点P_i的窗体空间为以点P_i为中心，边长为2*scale的立方体，其中scale表示窗体尺寸，i为大于0的整数；

D03、判断所述点P_i的窗体空间BoX(P_i)与所述肝脏血管树数据集Vas是否存在交集，如果不存在交集则直接跳过，如果存在交集则计算所述点P_i与所述交集中各段血管的距离，找出所述点P_i与所述交集中血管的最小距离，将该血管记为L_M，并且将所述点P_i划分至所述肝脏血管树数据集Vas的第M段血管L_M中，其中下标M为整数;

D04、计算所述肝脏血管树数据集Vas中像素值大于0的像素点的个数vasSum，并且计算vasSum与所述maskSum的比值ratio，即

D05、如果ratio<1，则重复执行D02～D04，否则结束操作，获得划分好的肝段。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1、在对肝脏血管三维图像进行结构分析时，根据血管供血区域划分肝段：首先采用三维骨架化方法细化肝脏血管三维图像，并且对肝脏血管三维图像的体素点进行标记，存储至树型数据结构中，以主从关系（即子节点和父节点的关系）来表示肝脏血管之间的关系，区分了肝脏血管系统的主支血管和分支血管，从而依据肝脏器内门静脉树或者肝静脉树的空间分布及其分支管径，对于其供血或者引流区域进行划分，并基于这种划分，将肝脏分段，计算其体积，医生可选择多个或者单个区域对肝段形态、体积及附近血管结构进行分析，也可选择任意分支静脉血管，测定该肝脏静脉分支血管的供血区域，根据肝静脉分支的供血区域，实现对肝脏的精确划分，指导肝脏的精准切除；

2、传统的二维分割易造成同一分割内容在邻近层面分割出来的图像变化较大或发生大的偏移二维分割方法而无法满足分割结果在三维空间的一致性要求，本发明在肝脏体图像的分割过程中，采用四叉树分割方法与Graph-Based图像分割方法相结合的方式，提高了分割的效率；

3、本发明在肝脏血管系统的分割中采用基于区域增长的阈值分割方法，不仅可以提取肝脏血管足够的血管细节信息，保证了分割的准确性，而且减少了过多噪音的引入，另外在肝脏血管系统分割完成后本发明还进行了两次消除噪音的处理，且进行平滑处理，保证了肝脏血管三维图像的质量；

4、在肝脏血管系统分割完成后的第一次消除噪声处理中，采用形态学操作并且在原有二维中值滤波的基础上提出了三维中值滤波，使得滤消除噪声的效果更好；

5、在肝脏血管系统的骨架化处理中采用12次迭代方法，细化效果更好；

6、在肝脏体的划分过程中在基于最小空间距离的肝段划分方法的基础上提出加窗的最小空间距离的肝段划分方法，使得划分获得的肝段更加符合实际的情况，更加具有实用价值。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是本发明的流程图；

图2是肝脏体的CT图像；

图3是肝脏体图像的分割结果；

图4是选取不同阈值提取出的血管树的效果图；

图5是肝脏血管树分割最佳阈值的估算曲线图；

图6是粗分割图像和经血管树追踪方法处理后的效果图；

图7是光滑处理前后血管树图像的对比图；

图8是骨架化处理中烧草模型和最大圆盘模型的示意图；

图9是骨架定义的示意图；

图10是三维空间各点之间的邻接关系示意图；

图11是6次子迭代方法中U方向的删除模板示意图；

图12是在MATLAB中N-，E-，S-，W-和D-方向的示意图；

图13是坐标经逆时针90°旋转后位置的示意图；

图14是12次子迭代方法的细化结果图；

图15是12次子迭代方法中12个删除方向的示意图；

图16是12次子迭代方法中方向旋转示意图；

图17是12次子迭代方法中删除模板的示意图；

图18是6次子迭代方法与12次子迭代方法的细化结果比较图；

图19是三维骨架化方法细化后的效果图；

图20是端点、曲线点和分叉点的表示图；

图21是三维骨架化方法细化、分叉点标记以及形成抽象树型结构的效果图；

图22是标记后的肝脏血管图像结果；

图23是肝脏血管的抽象树结构；

图24是迭代标记后的肝脏血管系统示意图；

图25是基于形态学膨胀的肝段划分方法和基于最小空间距离的肝段划分方法的效果图；

图26是经过加窗的最小空间距离的肝段划分方法处理后的效果图。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书（包括任何附加权利要求、摘要和附图）中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。

如图1所示，该肝脏血管系统的三维结构化描述方法由以下步骤组成：

步骤一、对肝脏体图像进行分割，获得肝脏三维图像的分割结果。

医学图像分割技术是医学图像处理和分析中的关键技术。所谓图像分割是指将图像中具有特殊涵义的不同区域区分开来；这些区域是互相不交叉的，每一个区域都满足特定区域的一致性。

定义：令集合R代表整个图像区域，对R的分割可以看作将R分成N个满足以下几个条件的子区域R₁,R₂,...R_K，其中下标i、j、k均为整数；

①

Y_{i = 1}^{K} R_{i} = R;

②对i，1，2，…，K，R_i是连通区域;

③对所有的i和j，i≠j，有R_i I R_j=φ;

④对i=1，2,…,K，有P(R_i)=TRUE;

⑤对i≠j，有P(R_i Y R_j)=FALSE。

上述条件①指出分割所得的所有子区域的并集应能包括图像中的所有像素，即是说分割时应将图像中的每一个像素都分进某一个子区域;条件②要求相同子区域内的像素应是连通的，应该总是根据一些分割的准则对图像进行分割的;条件③表明不同子区域是互不重叠的，也即是说一个像素不能同时被分到两个区域;条件④指出在分割后属于同一个区域中的像素应该具有某些相同的特性;条件⑤指出在分割后属于不同区域中的像素应该具有一些不同的特性。

通常在医学领域，由CT图像来显示肝脏体的状况，其中肝脏图像与周围的腹壁肌肉、隔肌、心脏等器官和组织相邻且缺乏良好的灰度对比，因而肝脏不同组织之间难以找到清晰的边界，肝脏图像分割难度大。Seong-Jae等人提出了利用肝脏的形态学及自身的一些特点，对肝脏图像进行了自动分割，但是该方法仅仅适用于组织没有粘连的理想情况。在2009年出版的“IEEE Transactions onInformation Technology in Biomedicine(IEEE生物医学信息技术汇刊)”中由GangChen等人提出了改进的水平集方法分割MR图像中的肝脏。此外，研究人员还提出了动态轮廓分割方法，并且将动态轮廓法运用到肝及其他器官的分割中，取得了一些成果，但仍然不应用到肝与其他组织严重粘连的情况，将肝从周围连接的组织中分割出来。

近来，研究人员又提出了一种GTaph-Based图论图像分割方法，其特点不仅在于按视觉重要性进行分割，即强调低变区域细节和忽略高变区域细节，使得分割结果较好，而且分割效率最高。该Graph-Based图论图像分割方法的主要思想如下:

给定带权图(即带权值图像)G=(V，E)，其中V表示该带权图的像素点集合，且也可以表示为节点的集合，各节点均为集合V的子集，E表示节点集V中两像素点之间的边集。令v_i，v_j∈V，分别表示像素点集合中的两个像素点，则边(v_i，v_j)对应于一个权值w((v_i，v_j))，其中w((v_i，v_j))用于表示两节点之间的相似度(如灰度、颜色等局部特征)。令Graph-Based图论图像分割方法的分割结果为S，则

在上述公式中，C_i、C_j分别为V的子集，即节点，对应于一个连通分量。

在图像分割中虽然采用不同的分割方法获得的分割效果有所不同，但是分割的基本要求是分割所得区域内部像素点之间差异较小，不同区域之间像素点的差异较大，即分割结果S中各区域内部两像素点的边对应的权值相对较小，分属于不同区域的两像素点的边对应的权值相对较大。

传统的Graph-Based图论图像分割方法中令

定义子集C的内部差异Int(C)为连通分量最小生成树MST(C，E)中的最大权值，即

Int(C)=max_e∈MST(C,E)w(e)

令C_i、定义子集（即节点）C_i、C_j之间的差异Dif(C_i,C_j)为所有连接两连通分量的边的最小权值，即首先求取分属于连通分量C_i、C_j的像素点之间边的权值，再取求得权值中的最小权值作为子集C_i、C_j之间的差异Dif(C_i,C_j)，公式表示为：

Dif (C_{i}, C_{j}) = \min_{v_{i} &Element; C_{1}, v_{j} &Element; C_{2}, (v_{i}, v_{j}) &Element; E} w ((v_{i}, v_{j}))

如果子集C_i、C_j之间没有边相连，则Dif(C_i,C_j)=∞；如果Dif(C_i,C_j)大于C_i和C_j中任何一个内部差异，即Int(C_i)、Int(C_j)，则认为子集C_i、C_j之间差异较大，将对应的图像区域视为不同区域，否则将对应的图像区域合并。

为了控制判决门限，即Dif(C_i,C_j)与Int(C_i)、Int(C_j)的大小关系，Graph-Based图论图像分割方法引入了一个阈值函数τ，由此定义的分割判决为：

D (C_{i}, C_{j}) = \{\begin{matrix} 1, ifDif (C_{i}, C_{j}) > MInt (C_{i}, C_{j}) \\ 0, others \end{matrix}\}

MInt(C_i,C_j)=min(Int(C_i)+τ(C_i),Int(C_j)+τ(C_j))

其中min为取最小值的运算符，当D为1表示C_i、C_j不合并，D为0表示C_i、C_j需要合并。阈值函数τ控制的是当Dif(C_i,C_j)与Int(C_i)、Int(C_j)相差多大时，才认为区域间存在明显的差异。阈值函数τ引入的原因在于，当C中的元素数目较小时，Int(C)并不能很好的放映数据局部特征，因此定义τ为关于连通分量大小的函数，τ(C)=k/|C|，其中|C|表示连通分量的大小，即连通分量的元素个数，K为常数。

传统的Graph-Based图论图像分割方法，由于它是把图像的每个像素点作为带权图的节点，这将导致相似矩阵的规模成倍增长，使得分割速度减慢。因此在本发明中提出了一种结合四叉树分割和Graph-Based分割的方法，四叉树的原则是，首先我们检查每个区域，看它是否满足区域同质性标准，如果不满足则把此区域按照水平和垂直方向进行分割。其他没有必要细分的将会重复这个过程，直至各个邻域的灰度一致性达到给定的阈值。整个分割过程就是，首先对预处理好的图像进行四叉树分割，得到一些小区域，然后把这些小区域作为图的节点，区域间的灰度相似度作为权值，进行分割。

具体地，步骤一由以下步骤组成：

A01、对肝脏体图像进行预处理，获得原始图像：通常获得的肝脏体图像为CT图像，其格式形式为dicom形式，将肝脏体图像转换为Bmp形式，并且去除噪声；

A02、采用四叉树分割方法对所述原始图像进行初始分割，获得多个小区域，由于四叉树分割方法为成熟技术，在此不对其予以累述；

A03、以该小区域v_i作为Graph-Based图像分割方法中带权图G的节点，带权图G=（V，E），其中V=（v₁，v₂，…，v_n）表示所述小区域的集合，由n个小区域组成，E表示所述小区域之间的边集，由m条边组成，n、m均为大于1的整数，i为≥1的整数并且边集E中每条边对应于一个权值，该权值用于表示所述小区域之间的灰度相似度；

A06、如果节点v_i、v_j分别属于第p-1次分割结果S^p-1中的不同连通分量，则两节点之间的边为e_p=(v_i,v_j)，边e_p对应的权值为w(e_p)，分别求取两连通分量v_i、v_j的内部差异Int(v_i)、Int(v_j)，内部差异的求取与上述求取方法相同；

如果权值w(e_p)大于两连通分量的其中一个内部差异（即大于内部差异Int(v_i)、Int(v_j)中之一）则不合并节点v_i、v_j（即不合并小区域v_i、v_j），否则合并节点v_i、v_j（即合并小区域v_i、v_j），从而获得第p次分割结果S^p；

肝脏的实际分割结果如图2～3所示，图2为肝脏体的CT图像，图3是肝脏体图像的分割结果。

步骤二、对肝脏血管图像进行分割，获得肝脏血管三维图像的分割结果，其由以下步骤组成：

第一步、采用基于区域增长的阈值分割方法对肝脏血管系统进行分割，提取肝脏血管的信息。

在血管图像分割中阈值选择非常重要，在分割血管图时若阈值较大，噪声则较少，却会损失较多细节；若阈值较小，得到的血管细节较丰富却也会引入过多的噪声，如图4所示。因而，既能提取足够血管细节信息，保证其准确性，又能保持较少的噪声的分割阈值可以被称为“最佳阈值”。

本发明采用基于区域增长的阈值分割方法，主要步骤为：

交互选取最靠近肝脏的血管入口作为种子点；

从该种子点开始生长，种子点的亮度设为θ_s，生长规则是迭代地将当前体素26邻域点中亮度大于或等于θ_s的体素点保存至表L(θ_s)；

以表L(θ_s)作为新的种子点集，它们的所有邻域点中亮度大于或等于θ_s-1,保存至表L(θ_s-1)中；

不断递减亮度阈值，直到血管的亮度阈值递减到θ_s-end，表L(θ_s-end)中已有体素点超出血管系统为止，如图5所示。

第二步、采用形态学操作和基于空间的中值滤波操作去除噪音，获得粗分割图像集。

对于数学形态学操作，其基础思想是利用若干给定的简单形态，通过几个基本的运算以提取与这些形态有关的信息。它是以整体的概念对图像进行处理，因此形态性强，几何概念明确。特别对二值表示，由于它以逻辑运算为基本形式，适合于硬件实现，效率极高，因而进来收到广泛的重视。

数学形态学的基本操作是用于一个预定的形状，即所谓的结构元，作为探针，来测试被分析图像的空间特性。我们知道，通过收缩和扩张一个二值图像，来消除噪音或去掉不必要的缺陷，是一个早期的且仍被经常使用的方法。而根据结构元的形状进行收缩和扩大的形态学的基本操作，分别是腐蚀(erosion)和膨胀(dilation)。

数学形态学是基于集合论发展起来的，它把二值图像看成是二值平面上的点的集合X，指定一个二值图为结构元，它同样是一个集合，记为B,则用B膨胀X,记为XOB，可定义为集合X和B的并：

X &CirclePlus; B = Y B_{x}, x &Element; X

这里B_x={b+x|b∈B}是结构元B对于点x的平移。腐蚀是膨胀的对偶，用B腐蚀X可表示成集合的包含变换：

XΘB={y|B_y∈X}

这里

是B对原点的对称变换。腐蚀也可以表示成集合的差：

XΘB=∩X_–b，b∈B。

对于基于空间的中值滤波操作，其是一种基于排序统计理论的能有效抑制噪声的非线性信号处理技术。中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的像素值接近的真实值，从而消除孤立的噪声点。中值滤波的方法是取某种结构的二维滑动模板，将该二维滑动模板内像素点按照像素值的大小进行排序，生成单调上升（或下降）的为二维数据序列。二维中值滤波输出为g（x,y）=med{f(x-k,y-l),(k,l∈W)}，其中f(x,y)，g(x,y)分别为原始图像和处理后图像，med是取中值的运算符号，W为二维模板，通常为2*2，3*3区域，也可以是不同的的形状，如线状，圆形，十字形，圆环形等。

然而，一般的中值滤波均基于二维图像的处理，本发明针对肝脏血管三维图像，设计了三维过滤器。

具体地，第二步由以下步骤组成：

B02-2、该图像集合R中以三张图像作为一个图像单元R_i，将该图像集合R表示为R=（R₁，R₂，…，R_n），如果图像集合R中图像数量不满足3的倍数，则将剩余的图像通过二维中值滤波处理输出，构成的图像单元R_i则执行如下处理；

B02-4、将图像单元R_i中像素值不为0的像素点均投影至一张图像上，计算投影后像素点坐标为（x，y）的像素值大小，如果该像素点的周围存在26个邻域像素点则以该像素点为中心统计其周围26个邻域像素点的像素值是否为0，该26个邻域像素点中每存在一个像素值不为0的像素点，nCount加1；

B02-6、重复执行步骤B02-4～B02-5，直至完成图像集合R中所有图像单元的过滤，获得粗分割图像集。

第三步、采用基于空间连通域的血管树追踪方法对所述粗分割图像集进行三维血管树追踪，获得连通的三维血管树图像集。

从三维的角度来看，肝脏各血管系统之间是独立的，而一个血管系统自身是连通的。这种连通性是指同一血管系统在断层图像的上下层之间存在着联系。

集合空间连通性定义：

对一个点P，记P的空间6邻域中的上下邻域{U,D}为N₂(P)。

对两个集合A与B，若B中至少存在一点x，使得N₂(x)∈B成立，则认为A与B空间连通。

粗分割血管图像集经过形态学膨胀腐蚀运算和空间滤波后去除了随机噪声。血管图的每个8连通域便是一个集合，因此血管树的追踪问题可转化为在空间搜索与血管树区域在空间连通的其他区域的问题。

具体地，第三步由以下步骤组成：

B03-4、重复执行步骤B03-3，直至所述栈Stack为空；

B03-5、重复执行步骤B03-2～B03-4，直至索引Index遍历完成；

第四步、采用三维形态学膨胀腐蚀操作对该连通的三维血管树图像集的表面进行平滑处理，获得处理后的肝脏血管三维图像。

在步骤B03中获得的连通的三维血管树图像集中仍然存在毛刺、表面不光滑的现象。本发明采用三维形态学膨胀腐蚀操作进行血管树的平滑，构造一个球体对连通的三维血管树图像集中所有图像进行一次一次滚动，由于球体的大小直接影响着光滑后的效果，通过具体的实验和对比，发现直径为11时效果最好。

采用三维形态学膨胀腐蚀操作进行平滑处理的步骤为：利用matlab工具箱，首先定义一个11*11*11的立方体，然后对该立方体按照球体的相关特性进行赋值，得到一个直径为11的球体，最后对连通的三维血管树图像集进行膨胀和腐蚀操作。

综上所述，在肝脏血管图像的分割实验中，利用西门子64排螺旋CT扫描机采集到门静脉和静脉两期胸部CT图像，首先采用步骤B01中基于区域增长的阈值分割方法进行粗分割，根据大量的实验得出，初始阈值选择为120，然后在初始阈值的周围选取像素值为100的像素点，根据增长规则，最后在像素值为100的像素点周围选择像素值为80的点，将选中的像素点的像素值置为255，将剩余的像素点的像素值置为0；并且采用步骤B02中形态学操作和基于空间的中值滤波操作去除噪声，获得粗分割图像集，效果图如图6的左侧图所示。由于粗分割图像集中仍然存在大量的噪声，本发明采用基于空间连通域的血管树追踪方法对粗分割图像集进行处理，获得连通的三维血管树图像集，效果图如图6的右侧图所示。

虽然图6的右侧图中三维血管树图像集已经连通，但是血管树表面仍然存在较多毛刺，为了消除此等毛刺，采用三维形态学膨胀腐蚀操作进行光滑处理，获得肝脏血管三维图像，如图7所示。从图7中可以看出，血管树表面的毛刺明显减少，同时表面存在小的枝桠，很明显噪音也得到了一定的处理和消除。

步骤三、对步骤二中获得的肝脏血管三维图像进行结构分析。

对于需要为病患做肝脏外科手术的医生而言，能够观察到直观清晰的三维空间中肝脏及其内血管和肿瘤的相邻关系则为制定准确的手术方案提供了有价值的信息。但在临床上遇到的很现实的问题是，当因切除肿瘤的需要而需切断肝内某一主支血管时，这一血管及其下一层级的血管所在的肝脏组织就成为了缺血或淤血的肝脏组织，这时就不得不将该部分的肝脏组织一并切除。所以需要预先确定被切除肝脏组织的范围，并在进行规则性的肝组织切除时，做到在三维空间结构上预先对肝段进行准确的划分，在手术中这些要求是必须要达到的，而仅仅有肝脏的三维图像是很难满足这方面的需求的。由此，在虚拟肝脏手术规划系统中，基于肝脏的三维图像，需要首先提取骨架，然后对肝内血管进行分级，并对上述问题进行一些研究和探讨。

具体地，步骤三由以下步骤组成：

第一步、采用三维骨架化方法细化处理该肝脏血管三维图像。

骨架化的基本原理为：在形态学中对象的骨架，是一种经过降维的物体形态的描述方式，不但把对象的轮廓和区域信息进行了组合，反映出对象重要的视觉上的线索；而且易于将骨架的线形连通结构转化为树或图的抽象形式，方便了对象的特征匹配。骨架用与原始对象形状的连通性和其拓扑结构相一致的细曲线来表示对象形状。骨架最初称为“中轴”，这个定义来源于“烧草”模型，中轴点即是对象边界上的火焰向内同时燃烧得到的相遇点。如图8的左侧图所示，火焰从对象边界上的两点开始向内部推进，燃过的痕迹随着时间流逝形成等距的同心圆，骨架点则是两圆或火焰前沿的交汇处。骨架的另一种定义为最大圆盘模型，该最大圆盘模型是一种更为直观、更为普遍的定义，是指完全包含在对象内部并且与对象边界至少相切于两点的圆，所有的最大圆盘的圆心的集合即为骨架，如图8的右侧图所示。

骨架的数学定义如下所述，设一个图形集合记为D，其骨架记为S(D)，则S(D)是D中所有最大圆的圆心的集合。S(D)所包含的最大圆其半径函数R是一个连续的实函数，函数值则等于每一个骨架点所对应的最大圆的半径长度。如图9所示，对象的骨架就是其所包含的所有最大圆的圆心的集合，最大圆与边界至少有两个切点。

为了保证血管骨架化结果的连续性和准确性，骨架化过程中需要考虑到下述四个较为关键的条件：①必需均衡对称的删除图像对象表层的像素或体素，才能保证最后剩下来的骨架位于管道的中心线；②在对像不断变细的过程中，原始形状的拓扑属性需要保持不变；③有时管道表面是一个不光滑的曲面，当曲面存在一个较小的突起部分形成了边界噪声时，会形成一条错误的骨架分支，即“无关”骨架，也称为伪骨架，因此需要能够对对象的边界噪声不甚敏感，减少伪骨架的形成；④骨架应该是单像素的宽度，尽可能的减少算法的迭代次数。

目前，计算骨架的方法主要集中在了两个方向：第一类是细化的方法，在保证拓扑结构没有改变的约束下，迭代的删除表层的像素或体素，得到最后剩下来的骨架。通过此类方法得到的骨架能够保证连通性，在字符的识别上有重要的应用。第二类是基于距离变换的方法，首先对对象进行距离变换的计算，计算时使用不同的距离标准，这样就计算出了不同的距离分布，然后寻找距离梯度发生突变的点作为骨架线。

在三维拓扑中，令x=(x₁,x₂,x₃)和y=(y₁,y₂,y₃)是三维空间Z³中有整数坐标的两个点。这两个点之间的欧几里得距离是

如果||x-y||≤1，则x和y为6-邻接，如果

则x和y为18-邻接，如果则x和y为26-邻接，采用N_j(p)表示点p的j-邻接点集，j=6,18,26，如图10所示。如果点x∈X且x∈N_j(p)，则称点p为对非空点集X的j-邻接。在二值三维图像中，黑点赋值为1，白点赋值为0。如果一个黑点6-邻接两个白点，则称这个黑点为边界点，对于一个边界点p，如果图10中标记为U的点为白点，则称点p为U-边界点。采用相同的方法还可以用同样方式定义N-，E-，S-，W-和D-边界点。

三维细化算法的目标是提取中心线并把对象收缩为体素链表。一个体素链表包含三种点：端点，曲线点和交叉点。如果一个黑点p仅m-邻接一个黑点，则称它为端点；如果一个黑点p有两个m-邻接的黑点，则称点p为曲线点；如果一个黑点p有超过两个的m-邻接的黑点，则称点p为交叉点。

目前，在骨架化过程中常见的细化方法有6次子迭代方法和12次子迭代方法。6次子迭代方法基于非端点的布尔特征，每个体素的新的赋值依赖于它的3×3×3邻接点。删除点意味着令黑点变为白点而白点保持不变，而删除点的情况是由每次子迭代所使用的匹配模版来描述的。如果一个黑点的3×3×3邻接点至少和一个匹配模版相匹配，则这个黑点需要被删除。

指定给方向U的模版给定了U-边界点的删除条件，模版如图11所示。在模版中用三种元素描述，黑点为“1”，白点为“0”，黑点或白点皆可则用“·”。在标记为“x”的点中至少有一个点为黑点。

模版M1可用如下布尔表达式表示，M1为TRUE则说明匹配成功。其中n(i,j,k)为0或1，i,j,k可取1,2,3。

对于N-，E-，S-，W-和D-边界点，需要先进行旋转。在利用MATLAB实现时，需要先进行方向上的统一。在MATLAB中，令x方向对应East，y方向对应North，z方向对应Up，如图12所示。从方向E向方向U旋转，则可以看成是对E-U平面上的点做逆时针90°旋转。图13中，圆括号中的坐标表示对应坐标的体素经逆时针90°旋转后的位置。

6次子迭代细化方法表述如下：

其中，T代表一次子迭代，并且每次子迭代是并行执行的，这意味着所有满足本次子迭代删除条件的黑点会在同一时刻删除。图14为6次子迭代细化结果的一个例子，左图为合成的三维图像，内容为大写字母“L”，右图为细化结果。

12次子迭代方法与6次子迭代方法相同，其检查3×3×3邻近的每个边界点。迭代步骤分为连续若干个子迭代，只有某一种类型的边界点可以在每个子迭代删除，因此，每个子迭代使用不同的删除规则。该12次子迭代方法被称为定向的或边界有序的。每个子迭代并行执行（即所有满足当前子迭代删除条件的边界点被同时删除）。由于在三维图像中有六个主方向，所以提出了6次子迭代定向细化算法。Srihari等人报告了早期的有序三维细化算法，使用12额外方向。这个算法也使用这12额外方向。主要的六方向和额外的12方向如图15所示。

由12个子迭代组成的定向细化方法可以简述为如下程序：

该方法采用定向策略，每步迭代包括12个连续并行削减操作，参照图13(b)所示的12方向，并提出了一个有序的删除方向列表：

<US,NE,WD,ES,UW,ND,SW,UN,ED,NW,UE,SD>。

这个顺序显示了一种对称性：有序列表可3个项目一组分为4组，每组包含所有六个主方向。因此，细化后的对象都在几何正确的位置上（即原对象的“中央”）。

该12次子迭代方法使用了曲线终点和曲面终点的如下特征。

定义1.在图片中如果集合N₂₆(p)/{p}包含只有一个黑点，则黑点p是一个曲线终点。

定义2.在图片中如果集合N₆(p)包含至少一组相对的白点，黑点p是一个曲面终点。（注意到每个曲线终点是曲面终点。）

子迭代中的一个可删除点是由一组3×3×3的匹配模板给定。黑点可删除如果一套模板中至少一个模板匹配于它。模板通常由三种元素描述，“●”（黑），“○”（白），和“·”（“不关心”），其中“不关心”在给定图片中匹配黑点或白点皆可，这样做是为了减少在模版中使用更多的符号。

指定给删除方向US的第一子迭代可以删除U-或S-边界点；和删除方向NE关联的第二子迭代删除N-或E-边界点，以此类推。

模板T_US（含14个模板）如图17所示，其属于曲线细化算法的第一子迭代，与删除方向US相对应。其它方向上的点需要先做旋转，才能和模板T_US、进行匹配。例如，在处理方向NE时，需要先将立方网格绕z轴顺时针旋转90°，即从棱NE位置转到棱ES位置，再绕y轴逆时针旋转90°，即从棱ES位置旋转到棱US位置，如图16所示，另外也可以从不同的轨迹对立体网格进行旋转，细化结果仍保持不变，原因在于立体网格的点的相对位置保持不变，不影响模板的匹配。

模板T_US指定给曲线细化算法的第一子迭代，对应于方向US。记号：每个标有“●”的位置匹配一个黑点；每个标有“○”的位置匹配一个白点；每个“·”（“不关心”）匹配或是黑色或是白色的点；标有“x”的位置至少有一个匹配黑点；标有“v”的位置至少有一个匹配白点；标有“w”的位置至少有一个匹配白点；两个标明“z”的位置匹配不同的点（其中一个匹配黑点，另一个匹配白点）。

模版T1可用如下布尔表达式表示，T1为TRUE则说明匹配成功。其中n(i,j,k)为0或1，i,j,k可取1,2,3。

图18为12次子迭代细化结果的一个例子，图(a)为合成的三维图像，内容为大写字母“A”，图(b)为12次子迭代算法细化结果，图(c)为6次子迭代算法细化结果。通过对比发现，由于12次子迭代算法将原来的6个方向细分为12个方向，在原来的垂直旋转的基础之上添加了水平旋转，能够得到更为细致的结果，因此12次子迭代算法对于有较复杂分叉点的对象能够更好的进行细化。在本发明中对肝脏血管三维图像的骨架化过程中采用12次子迭代方法进行细化。

综上所述，以肝脏血管三维图像中的肝脏门静脉切片为实验对象，效果图如图19所示，其中（a）表示肝脏门静脉的合成图像，（b）表示肝脏门静脉的合成图像的细化结果，（c）表示肝脏门静脉的合成图像与其细化结果的融合效果，（d）为融合效果的局部图。

第二步、分别采用端点、曲线点和分叉点标记肝脏血管三维图像的体素点，从而获得标记后的体素点，其中端点只具有一组邻接点，曲线点具有两组邻接点，分叉点具有3组以上邻接点。如图20所示，点p(2,2,2)周围共有4个邻接点，分别为点(2,2,1)，点(3,2,1)，点(3,1,1)和点(1,3,3)，此种情况下，将点p视为一个曲线点而非分叉点。

具体地，该第二步由以下步骤组成：

范围内的体素点，则counter加1；

C02-3、根据体素点p的邻接点组数counter来标记体素点p：如果counter=1则将体素点p标记为端点，如果counter=2则将体素点p标记为曲线点，如果counter>2则将体素点p标记为分叉点，标记后的肝脏血管图像结果如图22所示。

第三步、采用多叉树的数据结构来存储所述标记后的体素点，一个体素点的坐标存储在多叉树的一个结点中，且将多叉树中的父节点和子节点采用直线连接，从而形成肝脏血管的抽象树结构，如图23所示。

从采用三维骨架化方法细化处理、血管骨架的标记至形成抽象的树型结构，效果图如图21所示，其中（a）为血管树状结构，（b）为三维骨架化方法细化后的效果图，（c）为分叉点标记后的效果图，（d）为形成抽象树型结构的效果图。

第四步、在该肝脏血管的抽象树结构中将两个连接的分叉点视为一个连通域，采用连通域标记方法对各连通域进行标记，赋予唯一的数值，即将连通域中体素点均标记为相同的唯一的数值，并且采用迭代方法将与各连通域中体素点邻接的体素点均标记为相同的数值，从而区分肝脏血管系统的主支血管和分支血管，获得迭代标记后的肝脏血管系统，如图24所示。

步骤四、根据肝脏血管系统中的主支血管和分支血管标记以及医学需要，将肝脏体划分为多个肝段。

肝段内一个血管分支供养支配着一个相应的肝区域，将该肝区域称为肝亚段。肝段划分常用的方法有形态学膨胀和最小距离两种方法。

该基于形态学膨胀的肝段划分方法在数学上，膨胀定义为集合运算。A被B膨胀，记为

A &CirclePlus; B = {z | {(\hat{B})}_{z} IA &NotEqual; φ}

其中φ为空集，B为结构元素。但是，采用血管膨胀的方式，由于受到结构元素B的影响，造成了分界部分不准。

该基于最小空间距离的肝段划分方法采用L_i,i=1,2,3,...,n来表示第i段血管，V_i表示属于L_i血管段的所有点的集合，对于肝体中任意一点P到血管段L_i的最小距离可以表示为：

d_{i} (P) = \min_{v &Element; V_{i}} (dst (p, v))

此时，若d_M(P)=min(d_i(P)),i=1,2,L,n，则将点P(x,y,z)划分到L_M。

基于形态学膨胀的肝段划分方法和基于最小空间距离的肝段划分方法的效果图如图25所示，图25左侧两图表示基于形态学膨胀的肝段划分方法的效果图，右侧两图表示基于最小空间距离的肝段划分方法的效果图，由此可见，对断层图进行形态学膨胀会造成边界错误。

因此，本发明讨论基于最小空间距离的肝段划分方法，在原有的基础上提出了一种加窗的最小空间距离的肝段划分方法，具体地，由以下步骤组成：

D03、判断所述点P_i的窗体空间Box（P_i）与所述肝脏血管树数据集Vas是否存在交集，如果不存在交集则直接跳过，如果存在交集则计算所述点P_i与所述交集中各段血管的距离，找出所述点P_i与所述交集中血管的最小距离，将该血管记为L_M，并且将所述点P_i划分至所述肝脏血管树数据集Vas的第M段血管L_M中，其中下标M为整数；

经过加窗的最小空间距离的肝段划分方法处理后的效果图如图26所示，其中(a)是提取出的一组门静脉图，其划分根据的是医生的判断；(b)是由(a)通过3D空间距离方法一次迭代的结果；(c)为(a)经过同样方法两次迭代的结果；(d)是迭代完成获得的结果。

如图26所示，肝段之间的边界面形态为呈起伏状的曲面，也比较符合实验情况，并且根据临床统计的实际情况和肝外科专家的判断，本发明中的肝段划分方法划分出的肝段更加准确。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims

1.一种肝脏血管系统的三维结构化描述方法，其特征在于包括以下步骤：

A、对肝脏体图像进行分割，获得肝脏三维图像的分割结果；

其中所述步骤B02由以下步骤组成：

B02-2、所述图像集合R中以三张图像作为一个图像单元R_i，将所述图像集合R表示为R=（R₁，R₂，...，R_n），如果所述图像集合R中图像数量不满足3的倍数，则将剩余的图像通过二维中值滤波处理输出，构成的图像单元R_i则执行如下处理；

C、对所述肝脏血管三维图像进行结构分析：

C01、采用三维骨架化方法细化处理所述肝脏血管三维图像；

2.根据权利要求1所述的肝脏血管系统的三维结构化描述方法，其特征在于：所述步骤A由以下步骤组成：

A03、以小区域v_i作为Graph-Based图像分割方法中带权图G的节点，所述带权图G=（V，E），其中V=（v₁，v₂，...，v_n）表示所述小区域的集合，由n个小区域组成，E表示所述小区域之间的边集，由m条边组成，n、m均为大于1的整数，i为≥1的整数并且边集E中每条边对应于一个权值，该权值用于表示所述小区域之间的灰度相似度；

A04、根据边集E中各边对应权值的大小，将边集E进行排序E=（e₁,...,e_p,...e_m），其中1<p<m为整数；

A05、设定起始的分割结果为S⁰=（v₁，v₂，...，v_n），即将集合V中各小区域均作为起始分割的连通分量；

3.根据权利要求1所述的肝脏血管系统的三维结构化描述方法，其特征在于：所述步骤B03由以下步骤组成：

B03-4、重复执行步骤B03-3，直至所述栈Stack为空；

B03-5、重复执行步骤B03-2～B03-4，直至索引Index遍历完成；

4.根据权利要求1所述的肝脏血管系统的三维结构化描述方法，其特征在于：所述步骤C01中采用12次子迭代方法细化所述肝脏血管三维图像。

5.根据权利要求1所述的肝脏血管系统的三维结构化描述方法，其特征在于：所述步骤C02由以下步骤组成：

范围内的体素点，则counter加1；

6.根据权利要求1所述的肝脏血管系统的三维结构化描述方法，其特征在于：所述步骤D由以下步骤组成：