CN107480445A - 基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法 - Google Patents

Abstract

本发明公开了一种基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法，包括下述步骤：对模板CT的图像序列进行离散操作，生成对应的三维网格；统计训练集中空间连续的网格点在三维的肿瘤区域内联合出现的频数，得到对应的互关联规则；读入新病人的原发肿瘤数据，运用距离变换，得到各个体素点到原发肿瘤区域的空间距离；设定原发肿瘤区域的概率为1，每个体素点的侵犯概率与邻近体素的概率以及上一个观察状态的概率有关，从原发肿瘤区域开始不断迭代向外计算；最终计算得到所有体素点的被肿瘤侵犯的概率。本发明能辅助临床放射治疗计划的制定，实现不同发射剂量临床放疗靶区范围的自动确定，达到临床应用的勾画标准，提高医生的工作效率。

Description

基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法

技术领域

本发明涉及预测肿瘤生长，模拟肿瘤扩散研究领域，特别涉及一种基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法。

背景技术

放射治疗是目前治疗肿瘤的主流治疗方法，适用与多种肿瘤的治疗，，而放疗靶区范围的精确确定是放射治疗的成功的关键因素。精确的放疗靶区能有效减少肿瘤的复发概率，同时减少放疗带来的副作用。具体治疗过程中，医生将根据不同组织被肿瘤侵犯概率的不同，确定不同的放疗靶区和放射剂量。目前，放疗靶区的确定主要依赖医生临床经验，需要医生手动勾画。另外，ABAS软件提供了的靶区自动勾画系统，该系统原理是通过向原发灶肿瘤区域的三维方向外扩一定距离，得到目标放疗靶区，医生在该结果的基础上，根据临床经验，结合ICRU(InternationalCommissionRadiologicalUnits,83report,2010)给出的大概勾画标准，对结果进行修改。

目前手动勾画的靶区，由于客观原因(勾画时只有水平面的视角参照，缺乏上下层面的视角参考)，主观原因(不同放射医生经验不同造成的主观差异)的存在，造成不同医生之间的靶区的勾画的结果存在很大差异。对后续的治疗方案的确定，治疗效果的评估带来了负面的影响。另一方面，对于一个病人的放疗靶区勾画，通常需要2-3小时，需要很大的人力和时间成本，而ABAS的勾画系统过于简单，单一的距离扩展不适用复杂的结构和组织的情况。

发明内容

本发明的主要目的在于克服现有的技术的缺点和不足，提供了一种基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法。本发明已经成功地应用在鼻咽癌的靶区预测中，模拟了鼻咽癌的生长过程，得到了不同组织和结构的肿瘤侵犯概率，基于计算得到的概率分布图，可以精准地确定不同剂量的放疗靶区。该发明对后续治疗方案，放射剂量参数评估，晚期损伤的相关性研究提供一个一致的评价平台。

为了达到上述目的，本发明采用以下技术方案：

本发明公开了一种基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法，具体包括下述步骤：

S1、对模板CT的图像序列进行数据离散操作，网格化三维的CT图像得到对应的三维网格，读入训练集中医生勾画好的肿瘤轮廓文件，生成肿瘤原发灶区域的二值图像；

S2、对读入的原发灶肿瘤数据做预处理操作，挖掘邻近肿瘤之间的关联关系，统计相邻的体素联合出现的频数，存储统计的结果，建立互关联规则的数据库；

S3、读入新病人的原发肿瘤数据，得到对应的二值图像，按照设定好的转移次序和设定好的相邻体素状态转移公式，从GTV开始向外不断迭代计算相邻体素之间的状态转移概率，最终得到全局的肿瘤概率侵犯概率。

作为优选的技术方案，步骤S1中，还包括下述步骤：

S11、将模板CT的图像序列进行网格化，将CT图像序列划分为width ×height×zLen的三维网格，其中，width和height分别是一个CT切片横向和纵向上的格点数目，zLen为模板CT的层数，网格化得到的网格点即为体素点；

S12、所述二值图像为：载入医生勾画好的肿瘤轮廓文件，判断体素位置是否在肿瘤区域内部；对应肿瘤区域的图像体素点置1，不含肿瘤区域的体素格点置0；所述肿瘤区域将运用判断像素点是否在多边形内的方法，分开非肿瘤区域与肿瘤区域；具体的InPolygon的计算公式为：

IN＝InPolygon(x,y,xv,yv)

式中IN为输出的二值图像，其尺寸与原始输入CT图像I(x,y)相同，(x,y)为像素点的坐标，(xv,yv)为GTV轮廓上的点坐标。

作为优选的技术方案，步骤S2中，还包括下述步骤：

S21、预处理：去掉二值图像中肿瘤表面层面不稳定的数据，所述的不稳定数据是指：肿瘤最外层的数据，由于不同病人拍片时间的不同，导致了肿瘤表面数据存在误差；采用腐蚀算法进行预处理，具体公式为：

其中B(x)是腐蚀结构元素，对空间E中的每一点，结构元素B(x)平移之后， X是B包含于E的元素集合；

S22、挖掘邻近肿瘤区域的关联关系，从训练集中统计所有肿瘤中相邻体素组合之间的频数，存储结果并建立数据库。

具体包括下述步骤：

作为优选的技术方案，步骤S22具体为：

S221、给空间中所有三维网格点按照位置顺序编号，具体公式为 M_xyz＝z×height×width+x×height+y；其中，(x,y,z)是对应着网格点三维坐标； z为网格点所在的CT层面，(x,y)为格点在该CT切片中的坐标位置；

S222、统计空间上所有的连续的体素的组合，在N例训练集数据集中，在GTV内部，体素概率同时为1，且空间距离小的体素点被视为连续，病例N 的个数为2000，而三项集和四项集的频数被统计用于之后的模型计算；具体来说，对于一个三项集和四项集的组合顺序，按照格点的编号的升序排列(i＜j＜k)，当该组合的频数大于1时，该组合和对应的频数按照如下形式保存在数据集中，其中set₃通过查找索引(M_i,M_j,M_k)可以得到该索引对应的频数set₄通过查找索引(M_i,M_j,M_k,M_l)可以得到该索引对应的频数

作为优选的技术方案，步骤S3中，还包括下述步骤：

S31、设定初始状态的肿瘤侵犯概率，位于GTV内部的体素概率设置为1，位于GTV外部的体素概率设置为0，此时为初始观察状态的概率值，即t＝0，具体来说，点坐标为(x,y,z)的体素点M_xyz在第t观察时刻的概率记为

S32、运用距离变换计算所有体素到GTV的欧几里得距离，并四舍五入取整；得到各个点的距离索引；在GTV内部的体素点距离索引为0，点坐标为(x,y,z) 的体素点M_xyz记为d(M_xyz)；

S33、在已知t-1观察时刻的概率下，计算t观察时刻的各个体素点的肿瘤侵犯概率值；按照距离索引的顺序从小到大且从距离索引为1开始，依次计算肿瘤侵犯概率，的值不仅和其周围的体素转移概率有关，同时和该体素上一个观察时刻的概率相关，具体计算公式为：

其中(M_i,M_j,M_k)∈σ(M_xyz)，而σ(M_xyz)是三元组的集合，定义为：

σ(M_xyz)＝{(M_i,M_j,M_k)|

(M_i,M_j,M_k)∈set₃and||M_i,j,k-M_xyz||²＜d₁andd(M_i,j,k)＜d(M_xyz)} 即对于三元组中的各个元素，到M_xyz的距离均在一个常数范围之内，且各个元素的距离序列都小于当前体素M_xyz的距离序列；同时，(i,j,k)为升序排列，并且存在于之前构建的三项集数据库中；

S34、累计计算得到k次时刻的肿瘤侵犯概率之后，加权平均得到最后的全局概率分布图，将k次不同时刻的概率图累加之后除以k，得到最后的概率分布图。

本发明与现有技术相比，具有如下优点和有益效果：

1、目前对于平均一个病例的放疗靶区，人工勾画需要耗时4-6个小时，基于而本系统确定的放疗靶区平均一个病例的消耗时间在5分钟之内，大大减少了医生的工作量和工作效率。

2、现有的ABAS系统靶区勾画系统的工作原理是在硬性扩展用户指定的数值范围作为预防照射靶区。本发明基于大量数据，能精确地预测出肿瘤的侵犯概率，自动确定放疗靶区。预测结果更加精确和合理。

3、本发明基于大量经临床验证治疗有效的训练集，创新性地提出使用互关联规则来挖掘不同位置肿瘤之间的依赖关系，并使用马尔科夫模型模拟肿瘤生长的过程。本方法把原发灶肿瘤在局部位置的扩散行为看作是有互关联的，通过统计找出这一规律，统计子区域中联合出现的频繁项集，并存入数据库，用于新病人的肿瘤侵犯概率计算。

附图说明

图1是本发明方法隐马尔可夫模型的结构示意图；

图2是本发明方法在不同观察状态下的体素概率分布图；

图3是本发明方法的流程图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

本发明中的相关术语如下所示：

肿瘤区(grosstumorvolume，GTV)：指临床可见、可通过诊疗检查手段(包括CT和MRI)证实的肿瘤范围。

原发病灶临床靶区(clinicaltargetvolume，CTV)：除包含GTV外，还包括显微镜下可见的、亚临床灶以及肿瘤可能侵犯的范围。

互关联算法：是一种挖掘关联规则的频繁项集算法，主要包括候选集生成和情节的向下封闭检测两个阶段。

计算机断层扫描(ComputedTomography，CT)：利用精确准直的X线束、γ射线、超声波等，与灵敏度极高的探测器一同围绕人体的某一部位作一个接一个的断面扫描，具有扫描时间快，图像清晰等特点，可用于多种疾病的检查。

像素(pixel)：图像元素的简称。像素是构成数码影像的基本单元，是二维空间的最小单位。

体素(voxel)：三维图像体积元素的简称。数字数据于三维空间分割上的最小单位，用于三维成像、医学影像等领域。概念类似于二维空间的最小单位——像素。

隐马尔可夫模型(Hiddenmarkovmodel，HMM)：统计模型，用来描述一个含有隐含未知参数的马尔可夫过程。一般来说，包含一个隐藏状态和一个输出状态，如图1所示其中P_i是P点在ith观察状态下的侵犯概率，Y_i＝{0,1}，是P 点的输出状态。

如图3所示，本发明提出一种基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法，具体包括下述步骤：

S1、对模板CT的图像序列进行数据离散操作，网格化三维的CT图像得到对应的三维网格，读入训练集中医生勾画好的肿瘤轮廓文件，生成肿瘤原发灶区域的二值图。

S11、所述的网格化具体过程为：把不同层面的模板CT序列图像划分成长方体格点，即体素点，长宽高的物理距离分别为2*0.95756325mm× 2*0.95756325mm×3mm。得到格点数目为width×height×zLen的三维网格。 width和height的值，由训练集中N例病人中肿瘤所侵犯的区域决定，zLen为模板CT横断面的层数。

S12、所述二值图像为：判断网格点是否在肿瘤轮廓(多边形)内部，把三维网格中在肿瘤区域内部的网格点(体素)值设为1，肿瘤区域外部的网格点值设为0；

所述肿瘤区域将运用判断体素点是否在多边形内的方法(InPolygon)，分开非肿瘤区域与肿瘤区域。具体的InPolygon的计算公式为：

IN＝InPolygon(x,y,z,xv,yv,zv)

式中IN为输出的二值图像，其尺寸与原始输入横截面CT图像I(x,y,z)相同， (x,y,z)为像素,(xv,yv,zv)为GTV轮廓的坐标。

S2、对读入的原发灶肿瘤数据做预处理操作，为了挖掘邻近肿瘤之间的互关联关系，本发明统计了相邻的体素联合出现的频数，存储统计的结果，建立互关联规则的数据库；

S21、预处理:去掉二值图像中肿瘤表面层面不稳定的数据，所述的不稳定数据是指：肿瘤最外层的数据，由于不同病人拍片时间的不同，导致了肿瘤表面数据存在误差。本方法采取了腐蚀算法，具体公式为：

其中B是腐蚀结构元素，本发明选取了3×3的基础结构，通过该方法，肿瘤最外层不稳定数据被移除。

S22、从训练集中统计肿瘤中相邻体素之间的关联关系，建立数据库，具体实施过程为：

S221给空间中所有三维网格点按照位置顺序编号，具体公式为 M_xyz＝z×height×width+x×height+y；其中，(x,y,z)是对应着网格点三维坐标。_z为网格点所在的CT层面，(x,y)为格点在该CT切片中的坐标位置。

S222统计空间上所有的连续的体素的组合，在N例训练集数据集中，同时值为1(在GTV内部)的次数，在本例具体实施过程中，空间距离小的体素点被视为连续，病例N的个数为2000，而三项集和四项集的频数被统计用于之后的模型计算。具体来说，对于一个三项集和四项集的组合顺序，按照格点的编号的升序排列(i＜j＜k)，当该组合的频数大于1时，该组合和对应的频数按照如下形式保存在数据集中，其中set₃通过查找索引(M_i,M_j,M_k)可以得到该索引对应的频数set₄通过查找索引(M_i,M_j,M_k,M_l)可以得到该索引对应的频数

S3、读入新病人的原发肿瘤数据，得到对应的二值图像，按照设定好的转移次序和设定好的相邻体素状态转移公式，从GTV开始向外不断迭代计算相邻体素之间的状态转移概率，最终得到全局的肿瘤概率侵犯概率。

S31、设定初始状态的肿瘤侵犯概率，位于GTV内部的体素概率设置为1，位于GTV外部的体素点设置为0，此时为初始观察状态的概率值，即t＝0，具体来说，点坐标为(x,y,z)的体素点M_xyz,在第t个观察时刻的概率记为

S32、运用距离变换计算所有体素到GTV的欧几里得距离，并四舍五入取整。得到各个点的距离索引。在GTV内部的体素点距离索引为0，点坐标为(x,y,z) 的体素点M_xyz记为d(M_xyz)。

S33、在已知t-1观察时刻的概率下，计算t观察时刻的各个体素点的肿瘤侵犯概率值，具体计算过程如下：

按照距离索引的顺序从小到大(从距离索引为1开始)依次计算肿瘤侵犯概率，的值不仅和其周围的体素转移概率有关，同时和该体素上一个观察时刻的概率相关，具体计算公式为：

其中(M_i,M_j,M_k)∈σ(M_xyz)，而σ(M_xyz)是三元组的集合，定义为：

σ(M_xyz)＝{(M_i,M_j,M_k)|

(M_i,M_j,M_k)∈set₃and||M_i,j,k-M_xyz||²＜d₁andd(M_i,j,k)＜d(M_xyz)} 即对于三元组中的各个元素，到M_xyz的空间距离均在一个常数范围之内，且各个元素的距离序列都小于当前体素M_xyz的距离序列。同时，(i,j,k)升序排列，并且存在于之前构建的三项集数据库中。

条件概率值从之前构建的数据集获得，公式如下：

该值为对应的三项集和四项集的频数比例。在所有符合要求的三元组中，概率最大的组合被选取，另外，当前概率值还受到上一个观察状态的影响，综合两个因素，得到体素点M_xyz在t时刻肿瘤侵犯概率

S34、累计计算得到k次时刻的肿瘤侵犯概率，得到不同观察时刻的肿瘤侵犯概率分布图(如图2所示)，当k达到预设定的上限时，停止计算，归一化得到最后的全局概率分布图，在本例具体的实施过程中，取k的具体值为8，将k 次不同时刻的概率图累加之后除以k，得到最后的概率分布图。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (5)

1.一种基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法，其特征在于，具体包括下述步骤：

S1、对模板CT的图像序列进行数据离散操作，网格化三维的CT图像得到对应的三维网格，读入训练集中医生勾画好的肿瘤轮廓文件，生成肿瘤原发灶区域的二值图像；

S2、对读入的原发灶肿瘤数据做预处理操作，挖掘邻近肿瘤之间的关联关系，统计相邻的体素联合出现的频数，存储统计的结果，建立互关联规则的数据库；

S3、读入新病人的原发肿瘤数据，得到对应的二值图像，按照设定好的转移次序和设定好的相邻体素状态转移公式，从GTV开始向外不断迭代计算相邻体素之间的状态转移概率，最终得到全局的肿瘤概率侵犯概率。

2.根据权利要求1所述的基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法，其特征在于，步骤S1中，还包括下述步骤：

S11、将模板CT的图像序列进行网格化，将CT图像序列划分为width×height×zLen的三维网格，其中，width和height分别是一个CT切片横向和纵向上的格点数目，zLen为模板CT的层数，网格化得到的网格点即为体素点；

S12、所述二值图像为：载入医生勾画好的肿瘤轮廓文件，判断体素位置是否在肿瘤区域内部；对应肿瘤区域的图像体素点置1，不含肿瘤区域的体素格点置0；所述肿瘤区域将运用判断像素点是否在多边形内的方法，分开非肿瘤区域与肿瘤区域；具体的InPolygon的计算公式为：

IN＝InPolygon(x,y,xv,yv)

式中IN为输出的二值图像，其尺寸与原始输入CT图像I(x,y)相同，(x,y)为像素点的坐标，(xv,yv)为GTV轮廓上的点坐标。

3.根据权利要求1所述的基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法，其特征在于，步骤S2中，还包括下述步骤：

S21、预处理：去掉二值图像中肿瘤表面层面不稳定的数据，所述的不稳定数据是指：肿瘤最外层的数据，由于不同病人拍片时间的不同，导致了肿瘤表面数据存在误差；采用腐蚀算法进行预处理，具体公式为：

其中B(x)是腐蚀结构元素，对空间E中的每一点，结构元素B(x)平移之后，X是B包含于E的元素集合；

S22、挖掘邻近肿瘤区域的关联关系，从训练集中统计所有肿瘤中相邻体素组合之间的频数，存储结果并建立数据库。

4.根据权利要求3所述的基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法，其特征在于，步骤S22具体为：

S221、给空间中所有三维网格点按照位置顺序编号，具体公式为M_xyz＝z×height×width+x×height+y；其中，(x,y,z)是对应着网格点三维坐标；z为网格点所在的CT层面，(x,y)为格点在该CT切片中的坐标位置；

S222、统计空间上所有的连续的体素的组合，在N例训练集数据集中，在GTV内部，体素概率同时为1，且空间距离小的体素点被视为连续，病例N的个数为2000，而三项集和四项集的频数被统计用于之后的模型计算；具体来说，对于一个三项集和四项集的组合顺序，按照格点的编号的升序排列(i＜j＜k)，当该组合的频数大于1时，该组合和对应的频数按照如下形式保存在数据集中，其中set₃通过查找索引(M_i,M_j,M_k)可以得到该索引对应的频数set₄通过查找索引(M_i,M_j,M_k,M_l)可以得到该索引对应的频数

5.根据权利要求1所述的基于隐马尔可夫模型的肿瘤临床靶区侵犯概率计算方法，其特征在于，步骤S3中，还包括下述步骤：

S31、设定初始状态的肿瘤侵犯概率，位于GTV内部的体素概率设置为1，位于GTV外部的体素概率设置为0，此时为初始观察状态的概率值，即t＝0，具体来说，点坐标为(x,y,z)的体素点M_xyz在第t观察时刻的概率记为

S32、运用距离变换计算所有体素到GTV的欧几里得距离，并四舍五入取整；得到各个点的距离索引；在GTV内部的体素点距离索引为0，点坐标为(x,y,z)的体素点M_xyz记为d(M_xyz)；

S33、在已知t-1观察时刻的概率下，计算t观察时刻的各个体素点的肿瘤侵犯概率值；按照距离索引的顺序从小到大且从距离索引为1开始，依次计算肿瘤侵犯概率，的值不仅和其周围的体素转移概率有关，同时和该体素上一个观察时刻的概率相关，具体计算公式为：

其中(M_i,M_j,M_k)∈σ(M_xyz)，而σ(M_xyz)是三元组的集合，定义为：

σ(M_xyz)＝{(M_i,M_j,M_k)|

(M_i,M_j,M_k)∈set₃ and||M_i,j,k-M_xyz||²＜d₁ and d(M_i,j,k)＜d(M_xyz)}

即对于三元组中的各个元素，到M_xyz的距离均在一个常数范围之内，且各个元素的距离序列都小于当前体素M_xyz的距离序列；同时，(i,j,k)为升序排列，并且存在于之前构建的三项集数据库中。

S34、累计计算得到k次时刻的肿瘤侵犯概率之后，加权平均得到最后的全局概率分布图，将k次不同时刻的概率图累加之后除以k，得到最后的概率分布图。