CN104408734B

CN104408734B - 联合图像分割与形变配准技术的自适应目标区域转换方法

Info

Publication number: CN104408734B
Application number: CN201410764128.5A
Authority: CN
Inventors: 李登旺; 石雪; 陈进琥; 李洪升; 尹勇
Original assignee: Shandong Normal University
Current assignee: Shandong Normal University
Priority date: 2014-12-11
Filing date: 2014-12-11
Publication date: 2017-10-27
Anticipated expiration: 2034-12-11
Also published as: CN104408734A

Abstract

本发明公开了一种联合图像分割与形变配准技术的自适应目标区域转换方法，本方法自动分割危及器官及胸腹部肿瘤区域，采用形变配准算法，将CBCT图像与计划CT图像进行多模态配准，通过计算出的形变域将计划CT上的信息自动转移到CBCT上。本发明在计划CT与CBCT图像刚体配准后，为补偿不能恢复的形变范围，采用基于正交小波基函数的形变配准算法，通过计算出的形变域，实现危及器官与肿瘤的轮廓以及放疗计划从计划CT到CBCT的转移，精确定位靶区位置。

Description

联合图像分割与形变配准技术的自适应目标区域转换方法

技术领域

本发明涉及一种联合图像分割与形变配准技术的自适应目标区域转换方法。

背景技术

临床精确放射治疗的目标是在最大可能地提高肿瘤靶区射线吸收剂量的同时尽可能地降低危及器官的射线接收剂量。针对胸腹部肿瘤，目前临床上通常的做法是先CT定位，再按制定的放疗计划实施逐次照射。但在放疗过程中，由于器官运动和变形或各种客观因素使肿瘤和危及器官(主要是肝脏)偏离射野，导致肿瘤欠剂量照射而危及器官却过剂量。

利用近些年发展起来的图像引导放疗(Image guided radiation therapy，IGRT)技术，可对肿瘤靶区和危及器官进行实时监控，根据器官位置的变化调整治疗条件达到使照射野跟踪靶区的目的。临床上主要采用基于CBCT(Cone beam CT)的图像引导放疗技术，但目前该技术仅限用于调整摆位误差，即将计划CT与放疗前获取的KV级CBCT进行基于骨形的刚体配准，再通过计算放疗肿瘤靶区的位置偏移量来指导摆位；同时，目前的放疗计划也仅仅基于静态的计划CT，并未做到真正意义上的“自适应”。此外，计划CT上的危及器官及肿瘤也均由医生手工勾画，费时费力，具有主观性。

发明内容

本发明为了解决上述问题，提出了一种联合图像分割与形变配准技术的自适应目标区域转换方法，本方法自动分割危及器官及胸腹部肿瘤区域，采用形变配准算法，将CBCT图像与计划CT图像进行多模态配准，通过计算出的形变域将计划CT上的信息自动转移到CBCT上。

为了实现上述目的，本发明采用如下技术方案：

一种联合图像分割与形变配准技术的自适应目标区域转换方法，包括以下步骤：

(1)采用基于图谱的自动分割方法，对采集的N组CT图像去除伪影、噪声成分，采用梯度矢量流(Gradient vector flow，GVF)进行胸肋骨自动提取；

(2)在N组图像中随机选取一组作为参考图像，其余作为目标图像，将得到的肝脏训练集和胸肋骨分别进行肝脏图谱和胸肋骨图谱的创建；

(3)输入患者当前计划CT图像，采用GVF进行胸肋骨自动提取，将胸肋骨图谱作为参考图像，患者计划CT图像作为目标图像，以胸肋骨为特征，采用仿射配准算法与胸肋骨图谱进行刚体配准，得到患者计划CT图像中肝脏区域的初始值；

(4)使用基于正交小波基函数的形变配准算法对得到的肝脏区域与肝脏图谱进行形变配准，使肝脏边界更为精准，通过基于图谱的分割方法在患者计划CT上自动分割出肝脏区域；

(5)结合水平集模型与似然函数在计划CT上自动分割胸腹部肿瘤区域；

(6)以胸肋骨为特征，通过仿射变换将最新CBCT图像与计划CT图像进行刚体配准，为补偿刚体配准时不能恢复的形变范围，使用基于正交小波基函数的形变配准算法，计算出CBCT与计划CT的形变域，将计划CT中的危及器官和肿瘤的轮廓以及等剂量线信息自动转移到CBCT上。

所述步骤(2)中，具体方法为：

在N组满足条件的CT中随机选取一组作为参考图像，将其与其余的目标图像采用仿射变换进行配准，所有的图像对齐后，统计空间中每个像素存在的概率；N组图像数据配准后，得到一个概率图谱；采用多次循环迭代，以提高所建图谱的可靠性；

建立概率图谱后，根据具体应用环境选取合适的阈值，按照灰度值随图谱概率增大而增大的准则，把图谱从概率空间转化到灰度空间。

所述步骤(5)中，具体方法为：

采用似然函数与水平集模型相结合的能量函数，在患者计划CT图像上自动获取胸腹部肿瘤区域；似然函数E^p(φ)由感兴趣区域能量函数e_T(x)和背景区域能量函数e_B(x)两部分构成；用似然函数的有限高斯混合模型估计感兴趣肿瘤区域密度分布，似然函数的高斯混合模型估计胸腹部背景区域密度分布，通过以上似然函数的估计来划分图像概率密度；水平集模型用于获得边缘轮廓信息，防止肿瘤边界泄漏。

所述步骤(5)中，肿瘤区域的密度分布采用有限高斯混合模型估计，具体的表达式为：

其中，Ω_T为肿瘤区域，p(X)为肿瘤区域像素的概率，μ_T为高斯函数均值，σ_T为标准差，u₀(x)为患者CT图像；

背景区域密度分布采用高斯混合模型估计，具体的表达式为：

其中，Ω_B为背景区域，p(X)为背景区域像素的概率，和分别为高斯混合模型的均值和标准差，K为背景区域的分类个数，为分类的混合系数；

似然函数的具体表达式为：

其中，e₁(x)＝-e_T(x)，

距离正则化水平集模型的具体表达式为：

其中：μ、λ、α分别是正则化项、长度项、面积项的权重系数，

当似然函数与水平集模型最小时，水平集模型演化曲线收敛到感兴趣区域边界，水平集模型演化表达式为：

其中，div(·)为散度算子。

所述步骤(6)中，具体方法为：

将已导入加速器的计划CT作为目标图像，CBCT图像作为参考图像，采用GVF对CBCT图像进行胸肋骨自动提取，以胸肋骨为特征，利用仿射配准算法与患者计划CT进行刚体配准；为补偿刚体配准时不能恢复的形变范围，使用基于正交小波基函数的形变配准算法，获得形变域，将形变域三维小波分解后，使小波系数的初始值为0；采用非线性最小二乘LM(Levenberg-Marquadt，LM)优化算法，寻找合适的小波系数，根据所需精度控制在某个尺度上收敛；由估计到的小波系数得到参考图像与目标图像之间的形变域；将计划CT中的危及器官和肿瘤的轮廓以及等剂量线信息自动转移到CBCT上，进而实现剂量体积直方图(Dosevolume histograms，DVH)分析。

所述步骤(6)中，在形变配准时，使用Navier偏微分方程描述发生形变的各向同性物质的平衡状态，通过Navier偏微分方程建立形变能量函数，将形变域建模为小波系数的函数，其具体表达式为：

其中，(x,y,z)表示参考图像的空间坐标系，(x′,y′,z′)表示目标图像的空间坐标系，u表示形变域，c表示小波系数，u是c的函数。

所述Navier偏微分方程的外力具体表达式为：

H(X)是图像X的香农熵；

Navier偏微分内力方程具体表达式为：

其中，λ和μ是Lame系数；

形变能量函数具体表达式为：

E(c)＝inter(c)+w*exter(c)

其中，w是加权常数，本发明w为常数1，inter(c)表示内力约束函数，exter(c)表示外力约束函数。

本发明的有益效果为：

(1)在患者计划CT上自动分割危及器官与胸腹部肿瘤，避免手工勾画引入的人为误差，保证剂量使用的准确性；

(2)在计划CT与CBCT图像刚体配准后，为补偿不能恢复的形变范围，采用基于正交小波基函数的形变配准算法，通过计算出的形变域，实现危及器官与肿瘤的轮廓以及放疗计划从计划CT到CBCT图像的自适应转移，精确定位靶区位置。

附图说明

图1为本发明方法流程图；

图2为创建图谱流程图；

图3为基于小波基函数形变配准流程图。

具体实施方式：

下面结合附图与实施例对本发明作进一步说明。

对于胸腹部肿瘤，进行放射治疗时主要的危及器官为肝脏，通过50组CT图像建立正常肝脏图谱，采用基于图谱的分割算法实现患者计划CT中肝脏组织区域的自动分割；同时通过水平集模型与似然函数相结合的方法在计划CT上自动提取胸腹部肿瘤区域，达到自动分割计划CT的目的。之后使用基于正交小波基函数的形变配准方法，将计划CT中各种信息通过计算出的形变域自动转移到CBCT上。

具体步骤如图1所示:

第一步：为实现胸腹部肿瘤患者计划CT中肝脏组织的自动分割，本发明采用基于图谱的自动分割方法，根据患者的年龄、性别、身高、体重、病种、治疗状况、种族和区域等因素对100组增强型计划CT图像进行详细地分类筛选，选取其中50组条件相近，年龄在20-40岁之间且肝脏组织正常的CT图像，这些CT图像中的肝脏区域均被多位物理师和医学专家详细勾画与核对，可确保图谱的准确性。从中随机选择其中一组作为参考图像，其余49组作为目标图像，并去除伪影、噪声成分；采用梯度矢量流GVF进行胸肋骨自动提取。

第二步：对以上50组数据的肝脏训练集和胸肋骨分别进行肝脏图谱和胸肋骨图谱的创建。将随机选取的一组参考图像与其余的49组目标图像采用仿射变换进行配准，所有图像对齐后，统计空间中每个像素存在的概率。为了提高所建图谱的可靠性，本发明采用两次循环进行迭代，以减少由于采用特定初始参考图像而引入的误差。

使用上述同样的方法建立胸肋骨图谱。

建立概率图谱后，根据具体应用环境，选择合适的阈值，按照灰度值随图谱概率增大而增大的准则，将图谱从概率空间转化到灰度空间。选取阈值时要避免选取过高排除连带边缘组织，而选取过低又会包含其他组织的情况。

第三步：输入患者计划CT图像，采用GVF进行胸肋骨自动提取。以胸肋骨为特征，利用仿射配准算法与胸肋骨图谱进行刚体配准，在计划CT图像上得到肝脏区域的初始值。

第四步：为得到计划CT中更为精准的肝脏区域边界，本发明将计划CT作为目标图像，肝脏图谱作为参考图像，使用基于正交小波基函数的形变配准算法，实现患者计划CT肝脏组织与肝脏图谱的形变配准；根据以上几步所述，通过基于图谱的分割方法便可在计划CT上自动获取肝脏位置。基于正交小波基函数的形变配准算法与第七步相同，其具体详情参照步骤(7)。

第五步：构建似然函数与水平集模型相结合的能量函数，用来自动获取计划CT图像中胸腹部肿瘤区域。

本发明将计划CT图像用各个相互独立的区域组成，这些区域分为感兴趣区域Ω_T和背景区域Ω_B，感兴趣区域即计划CT图像中的肿瘤区域。采用的分割模型由似然函数和水平集模型两部分构成，具体表达式如下：

E(φ)＝E^p(φ)+E^e(φ) (1)

其中，E^p(φ)为似然函数，用以估计各像素的概率，进而区分该像素是属于肿瘤区域还是背景区域，似然函数包括感兴趣区域能量函数和背景区域能量函数两部分。E^e(φ)为局部空间边缘能量函数，用以防止分割时边界泄漏。

依据胸腹部CT图像全局密度分布多峰性特点，可用多个高斯分布的叠加拟合实现肿瘤分割。本发明使用有限高斯混合模型估计肿瘤区域的密度分布，具体表达式为：

其中，p(X)为肿瘤区域像素的概率，μ_T为高斯函数均值，σ_T为标准差，u₀(x)为患者CT图像。然后使用高斯混合模型估计背景区域密度分布，具体表达式为：

其中，p(X)为背景区域像素的概率，和分别为高斯混合模型的均值和标准差，K为背景区域的分类个数，为分类的混合系数。

结合水平集方法，用水平集模型的正负表示目标函数和背景函数，目标函数Ω_T＝{x,φ(x)≥0}，背景函数Ω_B＝{x,φ(x)＜0}。通过水平集模型使用阶梯函数H(φ)(Heaviside Function)以区分Ω_T和Ω_B，因此可得出似然函数表达式：

其中，e₁(x)＝-e_T(x)，

接下来本发明采用距离正则化水平集模型以获得肿瘤边缘轮廓信息，水平集模型的具体表达式为：

其中：为加权长度项，用以驱动零水平集曲线向目标边界演化；为加权面积项，用以加快演化速度。为正则化项，用以纠正水平集模型与符号距离函数的误差，保证水平集模型在整个曲线演化过程中最大程度地接近符号距离函数，避免在演化过程中重新初始化。μ、λ、α分别是正则化项、长度项、面积项的权重系数。

通过以上所述即可得出本发明中似然函数与水平集模型相结合的分割模型能量函数，当该能量函数最小时，水平集演化曲线收敛到感兴趣区域边界。水平集模型演化表达式如下：

其中，div(·)为散度算子。

第六步：通过治疗计划系统TPS(Treatment planning system，TPS)在已自动分割危及器官和肿瘤区域的计划CT上制定放疗计划，设计照射野。

第七步：患者开始接受治疗时，将已导入加速器的计划CT作为目标图像，治疗日CBCT图像作为参考图像，采用GVF对CBCT图像进行胸肋骨自动提取，以胸肋骨为特征，通过仿射变换与患者计划CT进行刚体配准，为补偿刚体配准时不能恢复的形变范围，本发明使用与上述第四步自动分割方法中相同的形变配准算法，即基于正交小波基函数的形变配准算法，计算形变域，将计划CT中各种信息自动转移到CBCT图像上，其结果可对离线状态下的放疗计划进行修正或更新，避免人工修订的误差，在CBCT中自动产生危及器官和肿瘤靶区的轮廓线，并自动生成相应等剂量线，进行剂量体积直方图分析，从而有效指导自适应放疗。基于正交小波基函数的形变配准算法的具体流程如图3所示。

计划CT与CBCT刚体配准后，对其进行小波分解，获得形变域。本发明采用接近三次正交的不完全正交三阶样条小波来表示形变域，以减少极小化能量函数过程的复杂度，将形变域三维小波分解后，小波系数的初始值设为0。采用非线性最小二乘LM(Levenberg-Marquadt，LM)优化算法递归的极小化公式，寻找合适的小波系数，以小波系数的尺度和子带作为循环参数进行递归，根据所需精度控制在某个尺度上收敛。最终由估计到的小波系数得出参考图像与目标图像之间的形变域。

小波系数的函数通过形变域建模，参考图像、目标图像、形变域，小波系数之间的数学表达如下：

其中，(x,y,z)表示参考图像的空间坐标系，(x′,y′,z′)表示目标图像的空间坐标系。形变域u＝(u¹,u²,u³)，c表示小波系数，形变域u为小波系数c的函数。作为目标参数，它通过极小化能量函数获得。

对多尺度小波系数进行估计来建立形变能量函数，函数由两部分组成，其中一项是内力形变能量函数，另外一项是外力约束能量函数。通过内力产生各种自由形变，通过外力约束达到平衡状态。在形变配准中，使用Navier偏微分方程描述发生形变的各向同性物质的平衡状态，具体的Navier偏微分方程表达如下：

其中，θ是三维膨胀系数，表达式为：

其中，X＝(x₁,x₂,x₃)^T是形变组织的三维空间坐标系，F＝(F₁,F₂,F₃)^T是作用于形变组织各处的外力，μ和λ均为形变常量，其取值取决于形变组织的自身性质。u＝(u₁,u₂,u₃)^T是配准算法所需恢复目标对象与参考对象之间的形变域。在发生形变的组织边缘，外力近似为零，对应医学图像的器官轮廓边界；而对于发生形变的其它部位，均需选择合适的内力和外力模型，并用适当的方法表示形变域。对于给定形变域的受力和表示方法，物体的形变域完全由形变常数和外力决定。

Navier偏微分方程前两项是造成组织发生形变的内力部分，将(9)代入(8)中前两项，得到的内力方程如下：

其中，λ和μ是Lame系数，其取值由弹性形变物质的自身性质决定。本发明设置λ＝0、μ＝1，这种情况下的医学图像性质满足各向同性物质属性，且不影响算法的精度。

本发明为多模态配准，由于归一化互信息与待配准图像数据对的重叠程度无关，鲁棒性较好，因此采用基于归一化互信息的相似性测度，具体表达式为：

其中，H(X)是图像X的香农熵。通过极小化模拟平衡状态的能量函数来估计目标函数中的小波系数，因此，外力用归一化互信息的倒数表示：

exter(c)＝1/NMI(X,X(u))(12)

由此便得到用于构成Navier偏微分方程来解决弹性形变问题的内力和外力部分，其中形变域用正交小波基函数表示，通过内力和外力的平衡最终得出形变域的能量函数，能量函数具体的表达式为：

E(c)＝inter(c)+w*exter(c) (13)

其中，w是加权常数，本发明w为常数1。因此，得到的形变域及其能量函数均为关于小波系数的函数。

为了减少被估计小波参数的数量，把小波系数根据尺度和子带的不同进行分类，而形变域便由不同尺度和不同子带的小波系数来描述。本发明采用递进策略寻找小波系数，估计大尺度上的小波系数后，再估计小尺度上的小波系数。不同尺度和不同子带小波系数的估计次序为：1-2-3-4-2²-3³-4⁴，具体如下表所示：

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims

1.一种联合图像分割与形变配准技术的自适应目标区域转换方法，其特征是：包括以下步骤：

(1)采用基于图谱的自动分割方法，对采集的N组CT图像去除伪影、噪声成分，实现CT图像中的肝脏组织区域的自动分割，继而采用梯度矢量流方法进行胸肋骨自动提取；

(2)在N组图像中随机选取一组作为参考图像，其余作为目标图像，将得到的肝脏区域和胸肋骨分别进行肝脏图谱和胸肋骨图谱的创建；

(3)输入患者当前计划CT图像，采用梯度矢量流方法进行胸肋骨自动提取，将胸肋骨图谱作为参考图像，患者计划CT图像作为目标图像，以胸肋骨为特征，采用仿射配准算法与胸肋骨图谱进行刚体配准，得到患者计划CT图像中肝脏区域的初始值；

(6)以胸肋骨为特征，通过仿射变换将CBCT图像与计划CT图像进行刚体配准，为补偿刚体配准时不能恢复的形变范围，使用基于正交小波基函数的形变配准算法，计算出CBCT与计划CT的形变域，将计划CT中的危及器官和肿瘤的轮廓以及等剂量线信息自动转移到CBCT上；

所述步骤(6)中，其具体方法为：

将已导入加速器的计划CT作为目标图像，CBCT图像作为参考图像，采用梯度矢量流方法对CBCT图像进行胸肋骨自动提取，以胸肋骨为特征，利用仿射配准算法与患者计划CT进行刚体配准；为补偿刚体配准时不能恢复的形变范围，使用基于正交小波基函数的形变配准算法，获得形变域，将形变域三维小波分解后，使小波系数的初始值为0；采用非线性最小二乘LM优化算法，寻找合适的小波系数，根据所需精度控制在某个尺度上收敛；由估计到的小波系数得到参考图像与目标图像之间的形变域；将计划CT中的危及器官和肿瘤的轮廓以及等剂量线信息自动转移到CBCT上，进而实现剂量体积直方图分析。

2.如权利要求1所述的自适应目标区域转换方法，其特征是：所述步骤(2)中，具体方法为：

在N组满足条件的CT中随机选取一组作为参考图像，将其与目标图像采用仿射变换进行配准，所有的图像对齐后，统计空间每个像素存在的概率；N组图像数据配准后，得到一个概率图谱；采用多次循环迭代，以提高所建图谱的可靠性；

3.如权利要求1所述的自适应目标区域转换方法，其特征是：所述步骤(5)中，其具体方法为：

4.如权利要求1所述的自适应目标区域转换方法，其特征是：所述步骤(6)中，在形变配准中，使用Navier偏微分方程描述发生形变的各向同性物质的平衡状态，通过Navier偏微分方程来建立形变能量函数，将形变域建模为小波系数的函数，其具体表达式为：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <msup> <mi>x</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mi>x</mi> <mo>+</mo> <msup> <mi>u</mi> <mn>1</mn> </msup> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>x</mi> <mo>;</mo> <mi>c</mi> <mo>)</mo> </mtd> </mtr> <mtr> <mtd> <msup> <mi>y</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mi>y</mi> <mo>+</mo> <msup> <mi>u</mi> <mn>2</mn> </msup> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>;</mo> <mi>c</mi> <mo>)</mo> </mtd> </mtr> <mtr> <mtd> <msup> <mi>z</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mi>z</mi> <mo>+</mo> <msup> <mi>u</mi> <mn>3</mn> </msup> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>;</mo> <mi>c</mi> <mo>)</mo> </mtd> </mtr> </mtable> </mfenced>

其中，(x,y,z)表示参考图像的空间坐标系，(x′,y′,z′)表示目标图像的空间坐标系，u表示形变域，c表示小波系数，u是c的函数，u¹、u²、u³分别表示目标图像上x,y,z三个方向上的形变域，它是小波系数c的函数。