CN103226837A - 一种观察子宫肿瘤放疗总剂量的分布图像的生成方法 - Google Patents

Abstract

本发明涉及X射线辐射的测量，具体涉及一种观察子宫肿瘤放疗总剂量的分布图像的生成方法，该方法由以下步骤组成：先分割出含施源器HDR CT图像的施源器区域，得到二值掩膜图像；利用掩膜图像对HDR CT图像做收缩变形及插值处理后，作为参考图像，以不含施源器的IMRT CT图像为浮动图像进行配准，用该变形场，对IMRT CT图像和IMRT剂量分布图像进行变形；再对变形后的IMRT CT图像的阴道区域进行扩大，作为浮动图像，以施源器区域为空气CT值的HDR CT图像为参考图像进行配准；用该变形场，对阴道区域经扩大的IMRT剂量分布图像进行变形；最后，将该变形结果与HDR剂量分布图像进行叠加并与HDR CT图像进行融合。本方法所生成的图像可准确地评价子宫肿瘤病人受照的总剂量。

Description

一种观察子宫肿瘤放疗总剂量的分布图像的生成方法

技术领域

本发明涉及X射线辐射的测量，具体涉及用于评估子宫肿瘤放疗辐射剂量的图像的生成方法。

背景技术

常规的放射治疗是以第一次获取的病人计算机断层图像（CT，Computed Tomography）作为计划图像，并只以该计划CT图像生成治疗计划，对病人实施治疗。但实际上，在每一次分次治疗中，病人的内部解剖结构都会发生改变，如技师摆位的误差，病人内部脏器的蠕动，肿瘤随着放疗的进行体积会缩小等，如果仍基于首次计划实施放疗，必将带来严重的剂量误差，并有可能使正常组织卷入高剂量区域。自适应放射治疗（ART，Adaptive Radiotherapy）技术可以解决上述问题。ART技术在治疗分次中获取当次治疗的病人CT图像，并基于该图像重新设计计划，对病人实施治疗。而图像变形配准技术（DIR，Image DeformableRegistration）是ART的一项关键技术，它主要用于将计划CT图像上的轮廓线推衍至当前CT图像，避免物理师重新勾画轮廓线的麻烦；另一方面，变形配准可以将各分次的剂量分布变形到同一解剖空间后进行叠加，从而评价病人所受的总剂量照射。

图像变形配准本质上是要寻找一种点对点的空间变换，将一幅图像（浮动图像）的解剖结构变形到另一幅图像（参考图像）上。目前变形配准算法可分为基于物理模型变形配准算法：如基于线性弹性模型的，基于流体力学模型的，基于光流场模型的；也有基于数学函数的变形配准算法：如基于小波变换的，基于薄板样条函数或B样条函数的，基于调和函数的等等。但几乎所有这些算法都基于这样一个假设：两幅图像上点是一一对应的，即一幅图像上的点，总能在另一幅图像上找到。如果在一副图像上出现了另一副图像并不存在的结构，一般常规的变形配准算法无法有效解决这类问题。

在宫颈癌放射治疗中，很多病人需要同时接受高剂量率近距离放射治疗（high-dose-ratebrachytherapy,HDR）以及调强放射治疗（Intensity Modulated Radiation Therapy，IMRT）。在进行HDR内照射的时候，需要先在病人阴道中放置一个施源器，然后将放射源通过施源器输送至阴道内临近肿瘤的区域进行照射，因此，进行HDR内照射放疗时获取的计划CT图像含有施源器。而常规的IMRT属于外照射放疗，病人并不需要放置施源器，因此获取的计划CT图像不包含施源器。

如果要评价该病人受照的总剂量，需要对IMRT和HDR的CT图像进行图像变形配准，得到变形场，然后利用该变形场将剂量矩阵映射到同一空间（例如，将IMRT的剂量矩阵变换到HDR的剂量矩阵空间，或相反），最后再进行剂量的叠加。对于这样的图像配准问题，显然违背了一般变形配准算法的基本假设，因为施源器在两幅图像上并不同时存在。因此，如果直接对IMRT CT和HDR CT图像进行变形配准，必然会产生严重的配准误差，从而导致剂量叠加的不准确。

Christensen G E等2001年公开了一种包含施源器的图像与不包含施源器的图像之间的变形配准方法（Christensen G E,Carlson B,Chao K S,Yin P,Grigsby P W,Nguyen K,Dempsey J F,Lerma F A,Bae K T,Vannier M W and Williamson J F2001Image-baseddose planning of intracavitary brachytherapy:registration of serial-imaging studies usingdeformable anatomic templates Int J Radiat Oncol Biol Phys51227-43），该方法由以下步骤组成：（1）将浮动图像和参考图像中膀胱、直肠、阴道和子宫（包括施源器）三组感兴趣区域分别勾画分割出来，并生成对应的二值掩模图像；（2）在浮动图像和参考图像中手动选择标记点，利用标记点进行初始化配准；（3）在步骤（2）的初始化配准的基础上，利用粘滞性流体力学算法对三组掩模图像进行配准，得到三组变形场；（4）将步骤（3）中得到的配准变形场进行叠加，得到总的变形场，然后对浮动图像进行变形。上述方法尽管也能得到配准结果，且该配准结果可应用于放疗辐射剂量的评估，但是由于该方法是直接用含施源器的图像进行配准，因此如果将其所得到的图像用于放疗辐射剂量评估，显然不准确。

发明内容

本发明所要解决的技术问题是提供一种用于观察子宫肿瘤放疗总剂量的分布图像的生成方法，该方法所生成的图像可准确地评价宫肿瘤病人受照的总剂量。

本发明解决上述问题的技术解决方案是：

一种用于观察子宫肿瘤放疗总剂量的分布图像的生成方法，该方法由以下步骤组成：

（1）读取含施源器的盆腔HDR CT图像、盆腔HDR剂量分布图像、不含施源器的盆腔IMRT CT图像和盆腔IMRT剂量分布图像；

（2）根据图像中每个像素点CT值大小的不同，采用区域生长的分割法分割出含施源器盆腔HDR CT图像中的施源器区域，并赋予所分割出的施源器区域像素值为0，得到施源器区域像素值为0的盆腔CT图像，然后将施源器区域像素值为0的盆腔CT图像中的施源器区域用1表示，其他区域用0表示，得到掩膜图像；

（3）利用变形场将掩膜图像进行收缩变形，该收缩变形的具体步骤如下：

在三维直角坐标系中，令收缩变形力F在x和y方向上的收缩变形力为掩膜图像在x和y方向上的梯度，收缩变形力F在z方向上的收缩变形力为0，然后求解下式（Ⅰ）所示的Navier-Stokes方程，得到收缩的变形场u和收缩变形的掩膜图像，

▽²v+▽(▽·v)+F=0     （Ⅰ）

上式（Ⅰ）中，▽为梯度，▽²为拉普拉斯算子，v为变形场的速度，F为收缩变形力；

（4）用步骤（3）所得的变形场u分别对施源器区域像素值为0的盆腔HDR CT图像进行收缩形变，得到收缩变形的盆腔HDR CT图像；然后搜索出收缩变形的盆腔HDR CT图像中对应于收缩变形的掩膜图像中像素值不等于0的点，并对所述不等于0的点进行线性插值，得到去除施源器后的盆腔HDR CT图像；

（5）以不含施源器的盆腔IMRT CT图像为浮动图像，以去除施源器后的盆腔HDR CT图像为参考图像进行变形配准，同时得到配准变形场r和变形配准后的盆腔IMRT CT图像；

（6）用步骤（5）得到的配准变形场r对盆腔IMRT剂量分布图像进行变形，得到变形配准后的盆腔IMRT剂量分布图像；

（7）分别搜索变形配准后的盆腔IMRT CT图像和变形配准后的盆腔IMRT剂量分布图像中对应于收缩变形的掩膜图像中像素值不等于0的点，并将搜到的点的像素值赋予0，得到变形配准后阴道内部像素值为0的盆腔IMRT CT图像和变形配准后阴道内部像素值为0的剂量分布图像；

（8）以变形配准后阴道内部像素值为0的盆腔IMRT CT图像作为浮动图像，施源器区域像素值为0的盆腔CT图像为参考图像进行变形配准，同时得到配准变形场q和变形配准后阴道区域扩大的盆腔IMRT图像；

（9）用步骤（8）得到的配准变形场q对变形配准后阴道内部剂量值为0的剂量分布图像进行变形，得到变形配准后阴道区域扩大的盆腔IMRT剂量分布图像；

（10）将变形配准后阴道区域扩大的盆腔IMRT剂量分布图像与盆腔HDR剂量分布图像进行叠加，再把叠加结果与含施源器的盆腔HDR CT图像进行图像融合，得到具有人体解剖结构的辐射总剂量的分布图像。

上述方案中，求解下式（Ⅰ）所示的Navier-Stokes方程的方法可以是常规的迭代方法、网格方法或算子分裂方法等，本发明人推荐采用迭代方法，该方法的步骤如下所述：

（A）采用下式（II）进行卷积运算，

ν^(k)＝φ_s*F^(k)     （II）

式（II）中，v为变形场的速度，φ_s为三维高斯核函数，F为收缩变形力，k为迭代循环次数的序号；

（B）采用下式（III）计算变形场的增量

$\frac{{\∂u}^{\left(k\right)}}{\∂t}=v^{\left(k\right)}-(v_x^{\left(k\right)}\frac{\∂u^{(k-1)}}{\∂x}+v_y^{\left(k\right)}\frac{\∂u^{(k-1)}}{\∂y}+v_z^{\left(k\right)}\frac{{\∂u}^{(k-1)}}{\∂z})---\left(\mathrm{III}\right)$

式（III）中，v为变形场的速度，u为变形场，k为迭代循环次数的序号，t为时间；

（C）采用下式（IV）计算本次迭代得到的变形场的增量并将其累加在上一次迭代所得到的变形场上，得到本次迭代的变形场，

$u^{\left(k\right)}=u^{(k-1)}+{\mathrm{\δ}}_t\frac{{\∂u}^{\left(k\right)}}{\∂t}---\left(\mathrm{IV}\right)$

式（IV）中，u、t和k与上式（III）相同；

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

（D）用当次迭代所得的变形场u，对上一次迭代中的掩膜图像进行收缩变形，得到本次迭代的掩膜图像；

（E）重复上述步骤（A）～（D），直至下式（V）中的终止参数ε_u小于或等于1.0×10^-14迭代终止，

${\mathrm{\ϵ}}_u=\mathrm{\Σ}\frac{\∂u^{(k-10)}}{\∂t}/\mathrm{\Σ}\left|u^{(k-10)}\right|-\mathrm{\Σ}\frac{{\∂u}^{\left(k\right)}}{\∂t}/\mathrm{\Σ}\left|u^{\left(k\right)}\right|---\left(V\right)$

式（V）中，u、t和k与上式（III）相同。

上述方案的步骤（5）以及步骤（8）中，所述的变形配准方法可以是常见的医学图像变形配准方法，如，基于线性弹性模型、流体力学模型、光流场模型的变形配准算法、基于小波变换、薄板样条函数、B样条函数或调和函数等变形配准算法，本发明人推荐采用属于光流场模型的变形配准算法中的Demons变形配准算法，该方法是通过迭代方法进行计算的，以步骤（5）为例，该步骤所述的变形配准方法的步骤如下：

1）采用下式（VI）计算配准变形场增量，

${\mathrm{dr}}^{\left(k\right)}=\frac{(m^{(k-1)}-s)\▿s}{{(m^{(k-1)}-s)}^2+{\left|\right|\▿s\left|\right|}^2}+\frac{(m^{(k-1)}-s)\▿m^{(k-1)}}{{(m^{(k-1)}-s)}^2+{\left|\right|\▿m^{(k-1)}\left|\right|}^2}---\left(\mathrm{VI}\right)$

式（VI）中，dr^(k)为配准变形场增量，m^(k-1)为上次迭代变形的浮动图像，s为参考图像，▽m^(k-1)为上次迭代变形的浮动图像的梯度场，▽s为参考图像的梯度场，k为迭代循环次数的序号；

2）按下式（VII）将本次迭代的配准变形场增量与上一次迭代变形场相加，得到当次的变形场r^(k)，

r^(k)＝r^(k-1)+dr^(k)     （VII）

式（VI）中，r为配准变形场，dr为配准变形场增量，k为迭代循环次数的序号；

3）用得到的变形场r^(k)对上次迭代得到的浮动图像进行变形，所得到的变形图像作为下一次迭代的浮动图像；

4）重复上述步骤1）～3），直至下式（VIII）中的终止参数ε_r小于或等于1.0×10^-14迭代终止，得到配准变形场l和配准后的盆腔HDR CT图像，

${\mathrm{\ϵ}}_r=\mathrm{\Σ}\left|{\mathrm{dr}}^{(k-10)}\right|/\mathrm{\Σ}\left|r^{(k-10)}\right|-\mathrm{\Σ}\left|{\mathrm{dr}}^{\left(k\right)}\right|/\mathrm{\Σ}\left|r^{\left(k\right)}\right|---\left(\mathrm{VIII}\right).$

式（VIII）中，r、dr和k与上式（VII）相同。

在步骤（8）中的变形配准方法与上述步骤（5）中的变形配准方法是一样的。

本发明具有以下有益效果：

（1）本发明巧妙地利用图像中施源器区域与其它区域CT值的显著差异，分割出盆腔HDR CT图像中施源器区域，然后对所得图像中施源器对应的区域进行收缩形变等一系列处理，因此，所得到的辐射总剂量的分布图像可准确地反应人体每一体素受照的总剂量及其分布情况，可为医生制订下一步放疗计划提供更准确的科学依据；

（2）根据CT值对图像进行分割，易于实现自动化，既可提高图像处理效率，又可避免手动分割而引入的配准误差；

（3）收缩形变过程中变形力的施加方式符合真实人体的物理形变过程，更有利于提高配准的精度；

（4）摒弃了现有技术需要对人体各组织器官进行标记、勾画分割和预配准的步骤，有利于进一步提高图像的处理速度。

附图说明：

图1：含有施源器的盆腔HDR CT图像。

图2：含有施源器的盆腔HDR剂量分布图像。

图3：不含有施源器的盆腔IMRT CT图像。

图4：不含有施源器的盆腔IMRT CT剂量分布图像。

图5：施源器区域的像素值为0的盆腔HDR CT图像。

图6：掩模图像。

图7：掩模图像收缩的变形力F示意图。

图8：收缩变形后的掩模图像。

图9：去除施源器后的盆腔HDR CT图像。

图10：变形配准后的盆腔IMRT CT图像

图11：变形配准后的盆腔IMRT剂量分布图像

图12：变形配准后阴道内部像素值为0的盆腔IMRT CT图像；

图13：变形配准后阴道内部像素值为0的剂量分布图像

图14：变形配准后阴道区域扩大的盆腔IMRT图像；

图15：变形配准后阴道区域扩大的盆腔IMRT剂量分布图像；

图16：将变形配准后阴道区域扩大的盆腔IMRT剂量分布图像与盆腔HDR剂量分布图像进行叠加后，再与含施源器的盆腔HDR CT图像进行图像融合，得到具有人体解剖结构的辐射总剂量的分布图像；

图17：利用Christensen G E等人的方法，变形配准后的盆腔IMRT剂量分布与不含有施源器的盆腔HDR CT剂量分布进行叠加后，再与含有施源器的盆腔HDR CT图像进行融合后的总剂量分布图像。

具体实施方式：

例1（子宫肿瘤放疗总剂量的分布图像的生成）

本实施例以某宫颈癌患者放射治疗过程中所获得的盆腔HDR CT图像（图1）、盆腔HDR剂量分布图像（图2）、盆腔IMRT CT图像（图3）和盆腔IMRT剂量分布图像（图4）详细描述本发明所述的方法，其具体步骤如下所述：

步骤1：读入如图1所示的分辨率均为256*256*108的含有施源器的盆腔HDR CT图像和如图2所示的盆腔HDR剂量分布图像以及如图3所示的分辨率均为256*256*108的不含有施源器的盆腔IMRT CT图像和如图4所示的盆腔IMRT剂量分布图像。

步骤2：由于施源器是金属的，且内置有X射线辐射源，因此施源器在图1所对应区域的CT值大于200HU，即像素值大于1200，并且像素值最大的点位于施源器内部。基于含施源器盆腔HDR CT图像的上述特征，本例以像素值1200为阈值，以像素值最大的点作为种子点，以像素点的像素值大于等于1200作为生长条件，采用区域生长的分割法，分割出施源器区域。然后，将所分割出的区域像素值赋予为0，即相当于将施源器区域填充为空气，其结果为如图5所示的施源器区域的像素值为0的盆腔HDR CT图像。再将图5中施源器区域用1表示，其他区域用0表示，得到如图6所示的掩膜图像。

步骤3：如图7所示，在三维直角坐标系中，令收缩变形力F在x和y方向上的收缩变形力为掩膜图像在x和y方向上的梯度，收缩变形力F在z方向上的收缩变形力为0。然后迭代求解发明内容中的式（Ⅰ）所示的Navier-Stokes方程，其迭代求解过程如下：

（3.1）按式（Ⅱ）用三维高斯核函数φ_s对力F^(k)进行卷积，得到本次迭代的变形场的速度场ν^(k)。式（Ⅱ）中的

在本实例中它是一个大小为17×17×17的三维矩阵，其运算式中，δ＝4，-8≤x≤8，-8≤y≤8，-8≤z≤8；

（3.2）先将v^(k)和上一循环得到的变形场u^(k-1)代入式（III）求得变形场增量

再将

和u^(k-1)代入式（IV）中求得本次迭代的变形场u^(k)。式（III）和（IV）中，当k＝1的时候，变形场的初始值u⁽⁰⁾＝0；

（3.3）用变形场u^(k)对上一次迭代中的掩膜图像Mask^(k-1)进行收缩变形，得到本次迭代的掩膜图像Mask^(k)；

（3.4）重复上述步骤（3.1）～（3.3），直至式（V）中的终止参数ε_u小于或等于1.0×10^-14迭代终止，得到收缩的形变场u以及如附图8所示的收缩变形后的掩膜图像。

步骤4：用步骤3所得到的变形场u，对施源器区域像素值为0的盆腔HDR CT图像进行收缩形变，得到收缩变形的盆腔HDR CT图像；然后分别搜索出收缩变形的盆腔HDR CT图像中对应于收缩变形的掩膜图像中像素值不等于0的点，并分别对这些点分别进行线性插值，得到如图9所示的去除施源器后的盆腔HDR CT图像。

步骤5：以图3为浮动图像，以图9为参考图像，采用属于光流场模型的变形配准算法中的Demons变形配准算法进行变形配准，所述的Demons变形配准算法是通过迭代方法进行计算的，其步骤如下：

（5．1）由于变形场增量dr^(k)是由上次迭代变形的浮动图像m^(k-1)与参考图像s以及浮动图像的梯度场▽m^(k-1)和参考图像的梯度场▽s共同决定的，本步骤采用式（VI）计算配准变形场增量dr^(k)；

（5．2）按式（VII）将得到的变形场增量dr^(k)与上一次迭代变形场r^(k-1)相加，得到当次的变形场r^(k)；

（5．3）然后用得到的变形场r^(k)对上次迭代得到的浮动图像m^(k-1)进行变形，所得到的变形图像作为下一次迭代的浮动图像m^(k)；

（5．4）重复上述步骤（5．1）～（5．3），直至式（VIII）中的终止参数ε_r小于或等于1.0×10^-14迭代终止，得到配准变形场r和图10所示的变形配准后的盆腔IMRT CT图像。

步骤6：用配准变形场r对图4进行变形，得到如图11所示的变形配准后的盆腔IMRT剂量分布图像。

步骤7：分别搜索图10和图11中对应于收缩变形后的掩膜图像（图8）中像素值不等于0的点，并将图10中相应位置的点的像素值赋予为0，即相当于将阴道内部填充为空气，得到变形后阴道内部像素值为0的盆腔IMRT CT图像（图12）。同时将图11中相应位置的点的像素值赋予为0，得到变形配准后阴道内部像素值为0的剂量分布图像（图13）；

步骤8：以图12作为浮动图像，图5为参考图像进行变形配准，采用属于光流场模型的变形配准算法中的Demons变形配准算法进行变形配准，所述的Demons变形配准算法是通过迭代方法进行计算的，具体计算方法可参照步骤5所述的迭代方法实施，最终得到配准变形场q和图14所示的变形配准后阴道区域扩大的盆腔IMRT图像。

步骤9：用步骤（8）得到的配准变形场q对图13进行变形，得到变形配准后阴道区域扩大的盆腔IMRT剂量分布图像（图15）；

步骤10：先将图15与如图2进行叠加，再把叠加结果用Amide软件与图1进行图像融合，得到如图16所示的具有人体解剖结构的总剂量分布图像。

例2（效果的分析比较）

本发明所提出的方法主要应用于放射治疗中的剂量叠加，从而评价病人在整个放疗期间的实际受照总剂量。在本实施例中，HDR近距离内照射和IMRT外照射的剂量分布分别如附图2和附图4所示。从附图4可以看到，，IMRT的剂量是仅仅是人体组织吸收的剂量分布，在进行剂量叠加之前，必须先把IMRT剂量在阴道的区域进行扩大，然后才能叠加在含有施源器的HDR剂量之上。

为了比较本发明的方法与现有方法的区别，将本发明与Christensen G E（参考文献Christensen G E,Carlson B,Chao K S,Yin P,Grigsby P W,Nguyen K,Dempsey J F,LermaF A,Bae K T,Vannier M W and Williamson J F2001Image-based dose planning ofintracavitary brachytherapy:registration of serial-imaging studies using deformableanatomic templates Int J Radiat Oncol Biol Phys51227-43）等2001年公开的一种包含施源器的图像与不包含施源器的图像之间的变形配准方法进行比较。参照该方法对两幅图像进行配准，步骤如下：（1）以不含有施源器的盆腔IMRT CT图像作为浮动图像，以含有施源器的盆腔HDR CT图像作为参考图像，将浮动图像和参考图像中膀胱、直肠、阴道和子宫（包括施源器）三组感兴趣区域分别勾画分割出来，并生成对应的二值掩模图像；（2）在浮动图像和参考图像中手动选择标记点，利用标记点进行初始化配准；（3）在步骤（2）的初始化配准的基础上，利用粘滞性流体力学算法对三组掩模图像进行配准，得到三组变形场；（4）将步骤（3）中得到的配准变形场进行叠加，得到总的变形场，然后分别对盆腔IMRT CT图像和盆腔IMRT剂量分布图像进行变形，得到变形配准后的盆腔IMRT CT图像以及变形配准后的盆腔IMRT剂量分布图像，然后用Amide软件，将变形配准后的盆腔IMRT剂量分布图像与包含施源器的盆腔HDR剂量分布图像（如图2）进行叠加得到总剂量分布图，再将总剂量分布图与图1融合，得到图17。

将图17与图16进行比较可见，图17中施源器内出现的X射线辐射源的高剂量区域（图17中箭头所指位置）实际上是人体组织的吸收剂量叠加在施源器内部剂量之上，这样，显然不能用于评价子宫肿瘤放疗总剂量及其分布。

Claims (2)

1.一种观察子宫肿瘤放疗总剂量的分布图像的生成方法，该方法由以下步骤组成：

（1）读取含施源器的盆腔HDR CT图像、盆腔HDR剂量分布图像、不含施源器的盆腔IMRT CT图像和盆腔IMRT剂量分布图像；

（2）根据图像中像素点CT值大小的不同，采用区域生长的分割法分割出含施源器盆腔HDR CT图像中的施源器区域，并赋予所分割出的施源器区域像素值为0，得到施源器区域像素值为0的盆腔CT图像，然后将施源器区域像素值为0的盆腔CT图像中的施源器区域用1表示，其他区域用0表示，得到掩膜图像；

（3）利用变形场将掩膜图像进行收缩变形，该收缩变形的具体步骤如下：

在三维直角坐标系中，令收缩变形力F在x和y方向上的收缩变形力为掩膜图像在x和y方向上的梯度，收缩变形力F在z方向上的收缩变形力为0，然后求解下式（Ⅰ）所示的Navier-Stokes方程，得到收缩的变形场u和收缩变形的掩膜图像，

▽²v+▽(▽·v)+F=0     （Ⅰ）

上式（Ⅰ）中，▽为梯度，▽²为拉普拉斯算子，v为变形场的速度，F为收缩变形力；

（4）用步骤（3）所得的变形场u分别对施源器区域像素值为0的盆腔HDR CT图像进行收缩形变，得到收缩变形的盆腔HDR CT图像；然后搜索出收缩变形的盆腔HDR CT图像中对应于收缩变形的掩膜图像中像素值不等于0的点，并对所述不等于0的点进行线性插值，得到去除施源器后的盆腔HDR CT图像；

（5）以不含施源器的盆腔IMRT CT图像为浮动图像，以去除施源器后的盆腔HDR CT图像为参考图像进行变形配准，同时得到配准变形场r和变形配准后的盆腔IMRT CT图像；

（6）用步骤（5）得到的配准变形场r对盆腔IMRT剂量分布图像进行变形，得到变形配准后的盆腔IMRT剂量分布图像；

（7）分别搜索变形配准后的盆腔IMRT CT图像和变形配准后的盆腔IMRT剂量分布图像中对应于收缩变形的掩膜图像中像素值不等于0的点，并将搜到的点的像素值赋予0，得到变形配准后阴道内部像素值为0的IMRT CT图像和变形配准后阴道内部像素值为0的剂量分布图像；

（8）以变形配准后阴道内部像素值为0的盆腔IMRT CT图像作为浮动图像，施源器区域像素值为0的盆腔HDR CT图像为参考图像进行变形配准，同时得到配准变形场q和变形配准后阴道区域扩大的盆腔IMRT图像；

（9）用步骤（8）得到的配准变形场q对变形配准后阴道内部像素值为0的剂量分布图像进行变形，得到变形配准后阴道区域扩大的盆腔IMRT剂量分布图像；

（10）将变形配准后阴道区域扩大的盆腔IMRT剂量分布图像与盆腔HDR剂量分布图像进行叠加，再把叠加结果与含施源器的盆腔HDR CT图像进行图像融合，得到具有人体解剖结构的辐射总剂量的分布图像。

2.根据权利要求1所述的一种用于观察子宫肿瘤放疗总剂量的分布图像的生成方法，其特征在于，所述求解式（Ⅰ）所示的Navier-Stokes方程的方法为迭代的方法，该方法的由以下步骤组成：

（A）采用下式（II）进行卷积运算，

ν^(k)＝φ_s*F^(k)     （II）

式（II）中，v为变形场的速度，φ_s为三维高斯核函数，F为收缩变形力，k为迭代循环次数的序号；

（B）采用下式（III）计算变形场的增量

$\frac{{\∂u}^{\left(k\right)}}{\∂t}=v^{\left(k\right)}-(v_x^{\left(k\right)}\frac{\∂u^{(k-1)}}{\∂x}+v_y^{\left(k\right)}\frac{\∂u^{(k-1)}}{\∂y}+v_z^{\left(k\right)}\frac{{\∂u}^{(k-1)}}{\∂z})---\left(\mathrm{III}\right)$

式（III）中，v为变形场的速度，u为变形场，k为迭代循环次数的序号，t为时间；

（C）采用下式（IV）计算本次迭代得到的变形场的增量并将其累加在上一次迭代所得到的变形场上，得到本次迭代的变形场，

$u^{\left(k\right)}=u^{(k-1)}+{\mathrm{\δ}}_t\frac{{\∂u}^{\left(k\right)}}{\∂t}---\left(\mathrm{IV}\right)$

式（IV）中，u、t和k与上式（III）相同；

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

（D）用当次迭代所得的变形场u，对上一次迭代中的掩膜图像进行收缩变形，得到本次迭代的掩膜图像；

（E）重复步骤上述（A）～（D），直至下式（V）中的终止参数ε_u小于或等于1.0×10^-14迭代终止，

${\mathrm{\ϵ}}_u=\mathrm{\Σ}\frac{\∂u^{(k-10)}}{\∂t}/\mathrm{\Σ}\left|u^{(k-10)}\right|-\mathrm{\Σ}\frac{{\∂u}^{\left(k\right)}}{\∂t}/\mathrm{\Σ}\left|u^{\left(k\right)}\right|---\left(V\right)$

式（V）中，u、t和k与上式（III）相同。

Images