基于图割的肺4D-CT肿瘤自动分割方法

1.基于图割的肺4D-CT肿瘤自动分割方法，其特征在于，该分割方法包括以下步骤：
(1)读取肺部4D-CT图像，该图像由多幅不同相位的肺部3D-CT图像组成；
(2)获取初始相位3D-CT图像上肿瘤的几何中心点，即初始相位的目标种子点，所述的初始相位3D-CT图像是指多幅不同相位的肺部3D-CT图像中第一幅3D-CT图像；
(3)根据步骤(2)获取的肿瘤的几何中心点以及肿瘤的大小，得到初始相位目标块；
(4)在4D-CT图像除去该初始相位3D-CT图像后所剩下的其余各相位的3D-CT图像上，分别选取与步骤(3)获得的初始相位目标块最相似的目标块；
(5)计算各相邻相位目标块之间的运动位移；
(6)将步骤(2)得到的初始相位的目标种子点与步骤(5)得到的运动位移相加，获得其他相位的目标种子点的位置；
(7)在步骤(3)和步骤(4)得到的各相位目标块上，利用步骤(2)和步骤(6)得到的各相位的目标种子点，获得肿瘤分割结果。
2.根据权利要求1所述的基于图割的肺4D-CT肿瘤自动分割方法，其特征在于：所述步骤(3)中的初始相位目标块是N×N×N大小的立方体，立方体的中心是肿瘤的几何中心点，N的大小由肿瘤大小决定。
3.根据权利要求2所述的基于图割的肺4D-CT肿瘤自动分割方法，其特征在于：所述步骤(4)利用完全搜索块匹配算法选取与步骤(3)获得的初始相位目标块最相似的目标块，具体包括：
(4.1)采用完全搜索块匹配算法中常用的平均绝对误差匹配准则，定义如下：
式中，(i,j,k)为位移矢量，fp(x,y,z)和fp-1(x+i,y+j,z+k)分别为在相位p上(x,y,z)处的灰度值和在相位p-1上(x+i,y+j,z+k)处的灰度值，N×N×N为块的大小，n为像素的编号；
(4.2)在整幅图像的范围内，逐个像素计算MAD值，从所有点中找出MAD最小值的点，该点对应的块即为我们要找的匹配块。
4.根据权利要求3所述的基于图割的肺4D-CT肿瘤自动分割方法，其特征在于：所述步骤(5)中的各相邻相位目标块之间的运动位移是指步骤(4.2)中MAD最小值的点处对应的位移矢量。
5.根据权利要求4所述的基于图割的肺4D-CT肿瘤自动分割方法，其特征在于：所述步骤(7)采用结合星形先验的图割算法进行分割肿瘤，具体包括：
(7.1)将步骤(2)和步骤(6)得到的各相位的目标种子点记作O，目标块的顶点记作背景种子点B；
(7.2)将图像映射成网络图，图像中的像素点对应网络图中的节点，并添加两个额外的顶点，一个为源点s，一个为汇点t，图像中相邻像素点p和q之间的关系对应网络图中节点之间的边；
(7.3)给网络图中每一条边赋予权值，设定权值的原则为：像素间的差异越小则权值越大，像素间差异越大则权值越小，P为所有像素点的集合，具体如下：
若 且 则p和q之间的边权值记为B{p,q}；若 则p和s
之间的边权值记为λ·Rp("bkg")；若p∈O，则p和s之间的边权值记为K；若p∈B，则p和s之间的边权值记为0；若 则p和t之间的边权值记为λ·Rp("obj")；若p∈O，则p和t之间的边权值记为0；若p∈B，则p和t之间的边权值记为K，λ为非负系数；
其中：
Rp("obj")＝-lnPr(Ip|O)     式(3)
Rp("bkg")＝-lnPr(Ip|B)    式(4)
Neighbor表示邻域系统，Ip表示像素p的灰度值，Pr(Ip|O)表示像素p属于目标的概率，Pr(Ip|B)表示像素p属于背景的概率，σ为缩放系数，dist(p,q)表示像素p和q的距离；
(7.4)定义能量函数：
E(A)＝λ·R(A)+B(A)     式(6)
其中：
Rp表示像素p的区域项，Ap表示赋予像素p的标号；
(7.5)在式(6)中加入星形先验，能量函数变为如下形式：
E(A)＝λ·R(A)+B(A)+S(A)    式(9)
其中：
Aq表示赋予像素q的标号；
(7.6)求解得到E(A)的最小值，对应于原始图像的最终分割结果。

基于图割的肺4D-CT肿瘤自动分割方法\n技术领域\n[0001] 本发明涉及医学图像处理技术领域，具体是指一种基于图割的肺4D-CT肿瘤自动分割方法。\n背景技术\n[0002] 4D-CT在传统3D-CT的基础上加入了时间轴，动态采集患者呼吸时的CT图像，而后通过图像重建和重新排序得到患者在不同呼吸相位的3D-CT图像。4D-CT不仅能够明显消除呼吸运动伪影，而且能够真实准确地反映肺部随呼吸运动的变化规律，为实现个体化精确放疗奠定了有力基础。\n[0003] 精确放疗的关键在于靶区的准确定位，靶区的精度会直接影响最终的放疗结果，确定靶区可以采用图像分割技术。然而，目前肺4D-CT肿瘤分割存在着两大主要问题：其一，\n4D-CT生成的图像数量巨大，通常有1000到2000张，甚至更多，如果仅靠放疗科医生手工分割每幅图像中的肿瘤来获得靶区，显然费时费力，也不现实；其二，肺部肿瘤经常与周围肺壁、膈肌、血管等正常器官组织发生黏连，若采用区域生长、边缘提取等自动分割技术，误分割的概率很大，直接影响靶区精度。因此，研究出一种自动化程度高，分割准确，鲁棒性强的肺肿瘤分割算法十分必要。\n发明内容\n[0004] 本发明的目的在于提供一种基于图割的肺4D-CT肿瘤自动分割方法，该方法不仅自动化程度高，而且能够确保分割准确。\n[0005] 本发明的目的可通过下述技术措施来实现：基于图割的肺4D-CT肿瘤自动分割方法，包括以下步骤：\n[0006] (1)读取肺部4D-CT图像，该图像由多幅不同相位的肺部3D-CT图像组成；\n[0007] (2)获取初始相位3D-CT图像上肿瘤的几何中心点，即初始相位的目标种子点，所述的初始相位3D-CT图像是指多幅不同相位的肺部3D-CT图像中第一幅3D-CT图像；\n[0008] (3)根据步骤(2)获取的肿瘤的几何中心点以及肿瘤的大小，得到初始相位目标块；\n[0009] (4)在4D-CT图像除去该初始相位3D-CT图像后所剩下的其余各相位的3D-CT图像上，分别选取与步骤(3)获得的初始相位目标块最相似的目标块；\n[0010] (5)计算各相邻相位目标块之间的运动位移；\n[0011] (6)将步骤(2)得到的初始相位的目标种子点与步骤(5)得到的运动位移相加，获得其他相位的目标种子点的位置；\n[0012] (7)在步骤(3)和步骤(4)得到的各相位目标块上，利用步骤(2)和步骤(6)得到的各相位的目标种子点，获得肿瘤分割结果。\n[0013] 本发明中，所述步骤(3)中的初始相位目标块是N×N×N大小的立方体，立方体的中心是肿瘤的几何中心点，N的大小由肿瘤大小决定。\n[0014] 本发明中，所述步骤(4)利用完全搜索块匹配算法选取与步骤(3)获得的初始相位目标块最相似的目标块，具体包括：\n[0015] (4.1)采用完全搜索块匹配算法中常用的平均绝对误差匹配准则，定义如下：\n[0016]    式(1)\n[0017] 式中，(i,j,k)为位移矢量，fp(x,y,z)和fp-1(x+i,y+j,z+k)分别为当前相位(x,y,z)在和上一相位在(x+i,y+j,z+k)的灰度值，N×N×N为块的大小；\n[0018] (4.2)在整幅图像的范围内，逐个像素计算MAD值，从所有点中找出MAD最小值的点，该点对应的块即为我们要找的匹配块。\n[0019] 本发明中，所述步骤(5)中的各相邻相位目标块之间的运动位移是指步骤(4.2)中MAD最小值的点处对应的运动矢量。\n[0020] 本发明中，所述步骤(7)采用结合星形先验的图割算法进行分割肿瘤，具体包括：\n[0021] (7.1)将步骤(2)和步骤(6)得到的各相位的目标种子点记作O，目标块的顶点记作背景种子点B；\n[0022] (7.2)将图像映射成网络图，图像中的像素点对应图中的节点，并添加两个额外的顶点，一个为源点S，一个为汇点T，图像中相邻像素点p和q之间的关系对应图中节点之间的边。\n[0023] (7.3)给图中每一条边赋予权值，设定权值的原则为：像素间的差异越小则权值越大，像素间差异越大则权值越小，具体如下：\n[0024] 若 且 则p和q之间的边权值记为B{p,q}。若p∈P, 则\np和s之间的边权值记为λ·Rp("bkg")；若p∈O，则p和s之间的边权值记为K；若p∈B，则p和s之间的边权值记为0。若p∈P, 则p和t之间的边权值记为λ·Rp("obj")；若p∈O，则p和t之间的边权值记为0；若p∈B，则p和t之间的边权值记为K。\n[0025] 其中：\n[0026]    式(2)\n[0027] Rp("obj")＝-lnPr(Ip|O)   式(3)\n[0028] Rp("bkg")＝-lnPr(Ip|B)   式(4)\n[0029]    式(5)\n[0030] (7.4)定义能量函数：\n[0031] E(A)＝λ·R(A)+B(A)   式(6)\n[0032] 其中：    式(7)\n[0033]    式(8)\n[0034] (7.5)在式(6)中加入星形先验，能量函数变为如下形式：\n[0035] E(A)＝λ·R(A)+B(A)+S(A)   式(9)\n[0036] 其中：    式(10)\n[0037]    式(11)\n[0038] (7.6)求解得到E(A)的最小值，对应于原始图像的最终分割结果。\n[0039] 与现有技术相比，本发明具有如下有益效果：\n[0040] (1)本发明的分割是在目标块上进行，而非整幅三维图像，从而去掉了尽可能多的干扰信息，背景种子点设置为目标块顶点和加入的星形先验，给分割提供了特定的约束条件，提高了分割精度；\n[0041] (2)本发明提出的方法只需医生在初始相位上选取目标种子点，即可获得同一呼吸周期内所有相位的肿瘤分割结果，提高了算法的自动化程度。\n附图说明\n[0042] 下面结合附图和具体实施方式对本发明做进一步的详细说明。\n[0043] 图1是本发明基于图割的肺4D-CT肿瘤自动分割方法的流程图；\n[0044] 图2是本发明实施例1中选取的0相位的肿瘤几何中心点以及初始目标块，初始目标块为右侧所显示的方形框，肿瘤几何中心点为方形框中的中心黑点，以冠状面为例；\n[0045] 图3是本发明实施例1中利用完全搜索块匹配算法依次得到从1至9其他9个相位的目标块以及目标种子点，以冠状面为例；\n[0046] 图4是本发明实施例1中相位4传统图割算法和本发明分割结果的对比图；从上至下分别为：冠状面；横断面；矢状面；从左至右分别为：待分割原始图像；传统图割算法分割结果；本发明分割结果；\n[0047] 图5是本发明实施例2中利用完全搜索块匹配算法依次得到从1至9其他9个相位的目标块以及目标种子点，以冠状面为例；\n[0048] 图6是本发明实施例2中相位2传统图割算法和本发明分割结果的对比图，从上至下分别为：冠状面；横断面；矢状面。从左至右分别为：待分割原始图像；传统图割算法分割结果；本发明分割结果。\n具体实施方式\n[0049] 实施例1\n[0050] 本发明基于图割的肺4D-CT肿瘤自动分割方法的具体流程如图1所示，具体步骤如下：\n[0051] (1)读取肺部4D-CT图像，该图像由10个不同相位的肺部3D-CT图像组成，各相位3D图像大小为320×224×92像素；\n[0052] (2)获取初始相位3D-CT图像上肿瘤的中心点，即初始相位的目标种子点；\n[0053] (3)得到初始目标块；\n[0054] (4)选取与初始相位目标块最相似的其他相位目标块；\n[0055] (5)估计出各相邻相位目标块之间的运动位移；\n[0056] (6)以步骤(2)得到的初始相位目标种子点和步骤(5)得到的运动位移为基础，获得其他相位的目标种子点；\n[0057] (7)在步骤(3)和步骤(4)得到的各相位目标块上，利用步骤(2)和步骤(6)得到的各相位目标种子点，获得肿瘤分割结果。\n[0058] 本发明中，所述步骤(2)中的肿瘤中心点为肿瘤的几何中心点。\n[0059] 本发明中，所述步骤(3)中的初始目标块是N×N×N大小的立方体，立方体的中心是肿瘤中心点，N的大小由医生决定，此例中N＝45。\n[0060] 本发明中，所述步骤(4)具体利用完全搜索块匹配算法选取。\n[0061] 上述步骤(4)具体包括：\n[0062] (4.1)采用完全搜索块匹配算法中常用的平均绝对误差匹配准则，定义如下：\n[0063]    式(1)\n[0064] 式中，(i,j,k)为位移矢量，fp(x,y,z)和fp-1(x+i,y+j,z+k)分别为当前相位(x,y,z)在和上一相位在(x+i,y+j,z+k)的灰度值，N×N×N为块的大小；\n[0065] (4.2)在整幅图像的范围内，逐个像素计算MAD值，从所有点中找出MAD最小值的点，该点对应的块即为我们要找的匹配块。\n[0066] 本发明中，所述步骤(5)中的运动位移是(4.2)中MAD最小值的点处对应的运动矢量。\n[0067] 本发明中，所述步骤(7)具体采用结合星形先验的图割算法分割肿瘤。\n[0068] 上述步骤(7)具体包括：\n[0069] (7.1)将步骤(2)和步骤(6)得到的各相位目标种子点记作O，目标块的顶点记作背景种子点B；\n[0070] (7.2)将图像映射成网络图，图像中的像素点对应图中的节点，并添加两个额外的顶点，一个为源点S，一个为汇点T，图像中相邻像素点p和q之间的关系对应图中节点之间的边。\n[0071] (7.3)给图中每一条边赋予权值，设定权值的原则为：像素间的差异越小则权值越大，像素间差异越大则权值越小，具体见表1。\n[0072] 表1：像素间的差异与权值之间的关系对照表\n[0073]\n[0074] 其中：\n[0075]    式(2)\n[0076] Rp("obj")＝-lnPr(Ip|O)   式(3)\n[0077] Rp("bkg")＝-lnPr(Ip|B)   式(4)\n[0078]    式(5)\n[0079] (7.4)定义能量函数：\n[0080] E(A)＝λ·R(A)+B(A)   式(6)\n[0081] 其中：    式(7)\n[0082]    式(8)\n[0083] (7.5)在式(6)中加入星形先验，能量函数变为如下形式：\n[0084] E(A)＝λ·R(A)+B(A)+S(A)   式(9)\n[0085] 其中：    式(10)\n[0086]    式(11)\n[0087] (7.6)求解得到E(A)的最小值，对应于原始图像的最终分割结果。\n[0088] 分割结果如图所示，对比利用原始图割算法分割出的结果图和利用本发明方法分割出的结果图。本发明的方法能够准确干净地分割出肿瘤。\n[0089] 除了视觉效果以外，本实例也通过量化指标客观评价本发明的有效性。以医生手工分割的结果作为评价本文实验结果的标准。分别采用Hausdorff距离(Hausdorff Distance，HD)和Dice相似系数(Dice Similarity Coefficient，DSC)去量化评价分割结果。\n[0090] Hausdorff距离定义如下：给定两个有限集合A＝{a1,a2,…,ap}和B＝{b1,b2,…,bq}，则A，B之间的Hausdorff距离定义为：\n[0091] H(A,B)＝max(h(A,B),h(B,A))   式(12)\n[0092] 其中,\n[0093] h(A,B)＝max min||ai-bj||   式(13)\n[0094] h(B,A)＝max min||bi-aj||   式(14)\n[0095] 式中||·||表示点集A，B之间的距离范数。函数h(A,B)和h(B,A)分别称为前向和后向Hausdorff距离。如果h(A,B)＝d，则表示A中所有点到B中点的距离不超过d，也就是说A中点都在B中点的距离为d的范围之内。Hausdorff距离表征了两个点集之间的不相似程度，因此Hausdorff距离越小则代表分割越精确。\n[0096] 利用式(12)计算原始图割算法分割结果和利用本发明方法分割结果的Hausdorff距离。可见，本发明的方法较原始图割算法，Hausdorff距离减小。\n[0097] Dice相似系数用来度量两个集合之间的相似性，定义如下：\n[0098]    (15)\n[0099] 由上式可知，Dice相似系数介于0和1之间，值越大，说明分割结果与标准越接近。\n[0100] 利用式(15)计算原始图割算法分割结果和利用本发明方法分割结果的Dice相似系数。从表2可见，本发明的方法较原始图割算法，Dice相似系数提高。\n[0101] 表2：实施例1中各相位传统图割算法和本发明分割结果的Hausdorff距离(HD)和Dice相似系数(DSC)对比表\n[0102]\n[0103] 实施例2\n[0104] 结合另一组具有10个相位的肺4D-CT数据描述本发明方法的处理过程，其具体步骤如下：\n[0105] (1)读取肺部4D-CT图像，该图像由10个不同相位的肺部3D-CT图像组成，各相位3D图像大小为320×224×86像素；\n[0106] (2)获取初始相位3D-CT图像上肿瘤的中心点，即初始相位的目标种子点；\n[0107] (3)得到初始目标块；\n[0108] (4)选取与初始相位目标块最相似的其他相位目标块；\n[0109] (5)估计出各相邻相位目标块之间的运动位移；\n[0110] (6)以步骤(2)得到的初始相位目标种子点和步骤(5)得到的运动位移为基础，获得其他相位的目标种子点；\n[0111] (7)在步骤(3)和步骤(4)得到的各相位目标块上，利用步骤(2)和步骤(6)得到的各相位目标种子点，获得肿瘤分割结果。\n[0112] 本发明中，所述步骤(3)中的初始目标块是N×N×N大小的立方体，立方体的中心是肿瘤中心点，N的大小由为肿瘤的几何中心点以及肿瘤的大小来决定，此例中N＝75。\n[0113] 本实例中，步骤(2)所述获取初始相位目标种子点，步骤(4)所述选取相似块，步骤(5)所述估计运动位移及步骤(7)分割肿瘤，其具体过程均与实例1相同，不再赘述。\n[0114] 分割结果如图所示，对比利用原始图割算法分割出的结果图和利用本发明方法分割出的结果图。本发明的方法能够准确干净地分割出肿瘤。\n[0115] 除了视觉效果以外，本实例也通过量化指标客观评价本发明的有效性。以医生手工分割的结果作为评价本文实验结果的标准。利用式(15)计算利用本发明方法分割结果的Dice相似系数。从表3可见与现有技术相比，本发明不仅能够提高自动化程度，还能够准确地分割出肿瘤。\n[0116] 表3：实施例2中本发明各相位的Dice相似系数(DSC)表\n[0117]\n相位 DSC\n相位0 0.8434\n相位1 0.8285\n相位2 0.8343\n相位3 0.8294\n相位4 0.8253\n相位5 0.8329\n相位6 0.8409\n相位7 0.8215\n相位8 0.8174\n相位9 0.8261\n平均值±标准差 0.8300±0.008\n[0118] 本发明的实施方式不限于此，可以根据实际需要进行修改，以适应不同的实际需求，因此，在本发明上述基本技术思想前提下，按照本领域的普通技术知识和惯用手段对本发明内容所做出其它多种形式的修改、替换或变更，均落在本发明权利保护范围之内。