图像变形配准技术(DIR, Deformable ImageRegistration)是自适应放射治疗(ART, Adaptive Radiotherapy)的一项关键技术^[1-3], 它可以将计划轮廓线形变至当前CT图像, 以重新优化生成新的治疗计划; 另一方面, 变形配准可以将各分次的剂量分布变形到同一解剖空间后进行叠加, 从而评价病人所受的总剂量照射。对于部分宫颈癌放疗病人, 需要同时接受高剂量率近距离放射治疗(High-Dose-Rate Brachytherapy, HDR)以及调强放射治疗(IntensityModulated Radiation Therapy, IMRT)。在进行HDR内照射的时候, 需要先在病人阴道中放置一个施源器, 放射源通过施源器输送至阴道内临近肿瘤的区域进行照射。因此, 进行HDR内照射放疗时获取的计划CT图像中含有施源器。而接受常规IMRT外照射放疗时并不需要放置施源器, 因此, 在IMRT CT图像上并不包含施源器。如果要评价该病人在整个放射治疗期间所接收的总受照剂量, 需要对IMRT和HDR这两幅CT图像进行图像变形配准, 得到变形场, 然后利用该变形场将剂量分布矩阵映射到同一空间(例如, 将IMRT的剂量分布变换到HDR的剂量分布空间, 或相反), 最后再进行剂量的叠加。因为施源器在两幅图像上并不同时存在, 这样的图像配准问题显然违背了一般变形配准算法的基本假设:两幅图像上点是一一对应的, 即一幅图像上的点, 总能在另一幅图像上找到。如果一副图像上出现了另一副图像并不存在的结构, 一般常规的变形配准算法无法有效解决这类问题。因此, 使用一般算法对IMRT CT和HDR CT图像进行变形配准, 必然会产生严重的配准误差, 从而导致剂量叠加的不准确。

目前关于这方面的研究并不多, 经笔者文献检索, 仅发现Christensen G E等^[4]提出了一种关于含施源器与不含施源器的CT图像之间的变形配准方法, 该方法先将膀胱、直肠、阴道和子宫(包括施源器)等感兴趣区域分别勾画分割出来, 然后分别进行配准。由于该配准方法仍然是直接利用含施源器的图像进行配准, 并且忽略了整个图像的拓扑结构连续性, 虽能在某些感兴趣器官区域能得到较理想的配准及剂量叠加结果, 但并不适用于病人受照总剂量的评价。

本研究提出一种新的联合图像分割、收缩变形以及配准的方法, 用于实现IMRT-HDR CT图像的变形配准。我们利用同一宫颈癌病人的IMRT CT和HDRCT图像和相应的剂量分布对算法进行了验证, 结果表明本文方法能有效的消除施源器对变形配准的影响, 实现子宫肿瘤病人IMRT CT及HDR CT图像的变形配准及受照总剂量的精确累加与评价。

1 材料与方法

为了消除HDR CT图像中施源器对于图像配准的影响, 本文提出以下联合图像分割、收缩变形以及图像变形配准策略, 主要包括以下步骤:首先分割出HDR CT图像(CT_HDR)中的施源器区域并将其填充成空气CT值得到, 并得到施源器掩膜图像(Mask); 然后, 对掩膜图像进行收缩变形得到Mask^deflated, 将收缩变形场作用于, 得到; 以为参考图像, IMRT CT图像(CT_IMRT)为浮动图像进行变形配准, 并根据收缩后掩膜Mask^deflated中的非零值在变形后的IMRT CT图像中的阴道区域平滑构建一个小的空气隙, 得到; 最后以为浮动图像, 为参考图像, 得到阴道区域经过扩大的并且解剖结构和CT_HDR一致的图像dCT_IMRT。对于剂量分布图像采用与以上步骤一样操作, 完成IMRT剂量和HDR剂量的累加。算法流程图如下图 1所示。

1.1 施源器自动分割

我们以圆柱形施源器内部的金属点作为种子点, 采用基于区域生长的方法分割出CT_HDR图像中的施源器区域, 将其CT值赋值为空气的CT值, 得到并得到以及施源器区域的掩膜图像。

1.2 施源器区域收缩形变

通过计算对施源器掩膜图像进行收缩变形的变形场, 将该收缩变形场收缩施源器区域, 得到去除施源器的HDR CT图像。

在三维直角坐标系中, 令收缩变形力F在x和y方向上的收缩变形力为掩膜图像Mask在x和y方向上的梯度, 收缩变形力F在z方向上的收缩变形力为0, 然后求解Navier-Stokes方程^[5-7]:

得到收缩的变形场u和收缩变形的掩膜图像Mask^deflated, 其中▽为梯度, ▽²为拉普拉斯算子, v为变形场的速度, F为收缩变形力。本文采用如下迭代方法求解Navier-Stokes方程:

(A) 利用三维高斯核函数Ø_s对力F进行卷积运算, 得到速度场v^(k)=Ø_s·F^(k);

(B) 计算变形场的增量:

(C) 更新变形场:

其中δ_u为每一次迭代允许的变形场的增量的最大值, 其取值范围为0＜δ_u＜1);

(D) 用当次迭代所得的变形场u^(k), 对k-1次迭代中的掩膜图像进行收缩变形, 得到变形更新后的掩膜图像;

(E) 重复步骤上述(A)~(D), 直至以下终止参数:

小于或等于1.0×10^-14迭代终止, 得到收缩变形场u和收缩后的掩膜Mask^deflated;

(F) 将u对进行收缩变形:

得到去除施源器后的HDR CT图像。

1.3 基于Demons的变形配准

本文算法需进行两次变形配准:第一次变形配准是初步将CT_IMRT图像上的解剖结构变形至去除施源器后的HDR CT ()图像域; 第二次变形配准主要是把CT_IMRT图像上的阴道区域扩充成施源器的形状, 以实现CT_IMRT与CT_HDR图像的匹配。我们采用Demons算法^[8-9], 对两幅图像进行变形配准, 配准变形场增量:

其中, k为迭代次数, dr^(k)为配准变形场增量, m^(k-1)为k-1次迭代变形的浮动图像, s为参考图像, ▽ m^(k-1)为k-1次迭代变形的浮动图像的梯度场, ▽s为参考图像的梯度场。为了提高配准速度, 我们采用基于金字塔的配准策略, 即配准从最低分辨率图像开始, 当次得到的变形场作为下一级配准的初始变形场, 直至完成最高分辨率级别图像的配准。

1.4 空气隙的构建

将CT_IMRT图像上变形至图像域后, 由于变形后的CT_IMRT图像上的阴道区域依然是闭合的, 需要再通过变形配准, 将其扩充成施源器的形状以进行剂量累加。我们通过搜索收缩后掩膜Mask^deflated中的不等于0的点, 并将变形后的CT_IMRT图像上相应点的CT值设为空气, 这相当于在IMRT图像中的阴道区域构建了一个小的空气隙, 以方便后续与进行变形配准, 对阴道区域进行扩充。

2 结果

我们利用5例宫颈癌患者放疗过程中接受高剂量率后装治疗(使用圆柱形施源器)所获得的CT_HDR图像及相应的HDR剂量分布, 以及接受IMRT治疗时获得的计划CT_IMRT图像及相应的IMRT剂量分布对本文算法进行了验证, 其中1例的实验结果如图 2所示。首先对CT_IMRT和CT_HDR图像和剂量分布图像进行刚性配准和重采样, 使所有输入图像的大小一致, 为256× 256×108。我们以CT值最大的点(施源器内导管的金属)作为种子点, 以像素点的CT值≥ 200 HU作为生长条件, 采用区域生长的分割法, 分割出施源器区域。然后, 将所分割出的区域赋予空气的CT值-1 000 HU, 得到 (图 2(c))。施源器的掩膜图像如图 2(d)所示, 1表示施源器, 其他区域为0。三维高斯核函数, 它是一个大小为17×17× 17的三维矩阵, 其中取经验值4。收缩变形场的示意图如图 2(e)所示, 收缩变形后的掩膜图像Mask^deflated如图 2(f)所示。利用得到的收缩变形场, 对图 2(c)进行收缩变形后的图像如图 2(g)所示。将CT_IMRT变形至图像域, 变形后的IMRT图像及其剂量分布如图 2(h)所示。对图 2(h)中的阴道区域构建空气隙后的图像及其剂量分布如图 2(i)所示, 对阴道区域变形扩充成施源器形状后的图像dCT_IMRT及其剂量分布如图 2(j)所示。如果利用传统的变形配准算法, 如Demons算法, 直接对CT_IMRT和CT_HDR图像进行配准, 变形配准后的IMRT剂量累加至HDR剂量后的分布如图 2(k)所示, 可以看到, IMRT上组织的剂量被叠加在HDR的施源器区域剂量上(如图 2(k)箭头所示), 出现了严重的剂量误差。而本文算法能准确地完成变形后IMRT剂量与HDR上相应组织剂量的累加(如图 2(l)箭头所示)。在对其余4例病人数据的实验中, 我们得到了类似的实验结果。本文算法是利用CUDA语言在图形处理器GPU中编程实现, 运行时间约为25 s, 其中施源器区域收缩变形约需时5 s, 两次Demons变形配准约需时20 s。

3 讨论

本研究针对子宫肿瘤自适应放疗中总剂量的累加和评价问题, 提出了一种新的联合图像分割、收缩变形以及配准的方法, 实现IMRT-HDR CT图像变形配准及剂量的累加。该方法首先分割出施源器区域, 然后通过求解Navier-Stokes方程对其进行收缩变形, 以消除施源器对图像变形配准的影响。并结合收缩变形的施源器掩膜图像, 利用基于Demons的图像变形配准方法, 将IMRT CT图像及其剂量变形至HDR CT图像域, 完成IMRT-HDR CT图像变形配准及剂量的累加。我们利用临床宫颈癌病人的IMRT CT和HDR CT图像及相应的剂量分布对算法进行了验证, 结果表明本文方法能有效实现子宫肿瘤病人受照总剂量的精确累加与评价。