正交小波基函数实现图像引导放射系统中计划CT与锥形束CT图像的形变配准

李登旺 王洪君* 万洪林 尹 勇 孙伟峰（山东大学信息科学与工程学院，济南 5000）

2（山东师范大学物理与电子科学学院，济南 250014）

3（山东省肿瘤医院，济南 250117）

4（中国石油大学信息与控制工程学院，东营 257061）

引言

近年来，医学影像形变配准及其应用的研究正在受到医学界和工程界的高度重视。如何提高形变配准算法的精度和速度，是图像引导放射治疗（image guided radiation therapy）和自适应放疗（adaptive radiation therapy）领域的研究热点。目前，多数配准算法复杂，耗时巨大，有些还存在需要进行手工标记、不能实现自动配准等诸多问题，束缚了它们在临床上的应用［1］。因此，研究能够有效描述形变特点的模型，提高形变配准的精度和速度，以达到临床应用的要求，是当前医学界和工程界一个富有挑战性的课题，具有非常重要的理论意义和临床价值。

例如，随着放射治疗技术的迅速发展，临床上已采用了以三维适形放射治疗（3D-CRT）和调强放疗（IMRT）为代表的精确放疗技术。目前，精确放疗通常的做法是先CT定位，由医生一次性制定放疗计划并依照此计划实施逐次照射。然而，由于在治疗过程中会不可避免的产生器官运动、变形和各种误差（如摆位误差等），这会使肿瘤（靶区）和危及器官偏离射野，从而导致肿瘤区域欠剂量照射和危及器官的过剂量照射。为解决这些问题，科学家将放射治疗机与成像设备结合在一起，在治疗时采集有关的图像信息，确定治疗靶区和重要结构的位置及相应的运动轨迹，并在必要时进行位置和剂量分布的校正，这称为图像引导放射治疗［2］。

近年发展起来的锥形束CT（cone beam CT，CBCT）可以通过多个途径确定和跟踪靶区并引导放疗，大大提高了IGRT的精度，成为目前IGRT开发和应用的热点。IGRT系统可以在直线加速器治疗以前，通过kV级CBCT系统得到病人的CBCT三维重建图像，与制定计划时的Planning CT图像（PCT）进行配准操作，得到CBCT重建图像与计划治疗区的偏移情况，指导摆位信息［3］。由于该系统要求配准速度要快，因此，目前系统中指导摆位的自动配准方法多是基于骨形特征的配准［2，4］，其精度依赖于骨形特征提取的准确性，但是，由于器官运动、变形以及治疗过程等方面因素会使肿瘤相对骨形的位置发生改变，同时，对于盆腔内部等不含骨组织的区域，无法进行骨形的配准。PCT与CBCT的形变配准研究是当前关于自适应放疗的热点问题，国际上多位学者已经进行了理论探索。基于B样条的形变配准方法［5］，基于FEM的形变配准方法［6］，基于Demons的形变配准方法［7］，基于光流学模型［8］的形变配准方法分别被用于PCT与CBCT的形变配准研究，用来指导基于CBCT的图像引导放射系统中的自适应放疗。由于计划CT与CBCT配准时，全局形变和局部形变同时存在，而且对于不同病例的形变差别较大，尽管上述各种方法在恢复形变时分别具有各自的局部优势，但用于自适应放疗恢复形变时，仍存在局限性，在精度、速度和鲁棒性上仍需要改进。

为了改进基于CBCT的图像引导放射系统中的形变配准算法，提出一种基于正交小波变换的计划CT与锥形束CT图像的形变配准方法。我们采用正交小波基函数由粗到精地表示形变域，低分辨率下的小波系数表示全局形变特征，高分辨率下的小波系数表示局部形变特征，这些小波系数由粗糙尺度到精细尺度逐步通过极小化形变能量函数估计得到，形变能量函数由Navier偏微分方程来设计。为了使自适应放疗的形变配准算法具有较高的鲁棒性，形变能量函数的外力部分采用归一化互信息作为约束。仿真实验结果及临床应用均验证了本算法的有效性。

1 Navier偏微分方程

在形变配准中，Navier偏微分方程可以作为描述发生形变的各向同性物质的平衡状态方程［9］。基于薄板样条，B样条等数学径向基函数的配准方法是通过插值和拟合理论来接近形变域，而基于Navier偏微分方程的方法可以通过内力使物质形成各种自由变形，并且通过约束力来保证恢复形变，可以得到形变域的数学物理解。因此Navier偏微分方程是描述非线性形变配准的有效方法。Navier偏微分方程基于连续力学模型，模型中器官或者组织的形变是由内力和外力共同作用而成，当形变过程完成时内力与外力达到平衡。达到平衡状态后的Navier偏微分方程表达为

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

式中，X=（x1，x2，x3）T是发生形变组织的三维空间坐标系，F=（F1，F2，F3）T是作用于形变组织各处的外力。形变常量μ和λ取决于发生形变组织的自身性质。u=（u1，u2，u3）T是配准算法所需恢复参考对象与目标对象之间的形变域。在发生形变组织或器官的边缘，外力近似为零，对应图像的轮廓或者边缘；而对于发生形变的其它部位，均需要外力和内力的平衡。对于给定的形变域表示方法和受力，物体的形变域完全由外力和形变常数来决定。为了利用Navier偏微分方程来解决弹性形变问题，需要适当选择外力和内力模型，并且用合适的方法来表达形变域，这对于能否准确恢复形变域具有决定作用。

平衡方程式（1）的前两项，是造成组织或者器官形变的内力。将式（2）带入式（1），得到如下的内力方程：

式中，λ和μ是Lame系数，它们的取值由弹性形变物质的自身性质决定。在不影响算法精度的前提下，为了简化计算的复杂度，选择λ=0和μ=1。

2 利用正交小波基函数实现形变配准

2.1 正交小波基函数表示形变域

临床上涉及到的病例种类较多，因此形变特征较丰富，时常是全局形变与局部形变共存，同一组织或器官的不同部位大形变区域与小形变区域共存。对形变域进行由粗及精描述，即对全局形变和局部形变分别描述，可以更加准确可靠地恢复形变域。同时，由粗及精的策略可以使配准过程避免陷入局部极值，提高配准的鲁棒性。传统由粗及精策略利用塔式数据结构［10］。小波（Wavelet）变换作为信号的一种多分辨率表示方法，自然适合对形变域由粗及精地描述［11］。

对于形变配准问题，将形变域表示为小波系数的函数，则参考图像、目标图像、形变域有如下数学表达：

式中，（x，y，z）是参考图像的空间坐标系，（x′，y′，z′）是目标图像的空间坐标系。形变域u=（u1，u2，u3）是小波系数c的函数，作为目标参数，它可以由极小化能量函数得到。如前所述，将参考图像和目标图像之间的形变域通过小波变换表示。由于正交小波变换具有能量集中特性，在表示形变特征时具有优势。因此，对于三维医学图像配准，形变域u（X），在三维子空间N=（Nx，Ny，Nz）中，用三维离散正交小波分解表示可得

式（5）中各个参数的含义与普通的三维正交小波分解含义相同［12］。表示形变域的正交小波基函数在三维网格N=（Nx，Ny，Nz）中，以2j为空间间隔进行平移变换。基函数采用的是一维小波基函数的三维张量积形式，并且由小波系数加权。把式（5）中对形变域参数的表示形式带入式（3），可得到Navier偏微分方程内力部分的线性组合形式。这一形式与各向同性的静态Navier偏微分方程的Wavelet-Galerkin离散矩阵成线性正比，运算上可以优化［13］。

2.2 归一化互信息形成外力约束

为了使Navier偏微分方程处于平衡状态，从而恢复形变域，需为式（3）中的平衡方程设计合适的外力。外力的设计需要利用配准数据对的特征。仅仅利用局部特征如曲面引导、对应点特征引导或其它局部相似特征引导等，虽然可以获得较快的配准速度，但由于全局信息的丢失，恢复形变将会受限。因此，为了能够对全局形变和局部形变同时恢复，选取适当的特征来产生外力，成为是否可以准确恢复形变域的关键所在。

在配准过程中，外力的作用是使待配准图像数据对的相似或者相同的区域达到空间上的对应。因此，可以利用参考图像和目标图像的相似性测度来模拟外力。基于体素相似性，特别是互信息的医学图像配准方法精度高，鲁棒性好，易实现全自动配准，并且可以实现多模态图像的准确配准［14］。同时，互信息不仅包含图像局部的信息，而且包含全局的信息，这对于恢复全局形变和局部形变具有重要意义。由于归一化互信息与待配准图像数据对的重叠程度无关，相对于基于互信息的配准，鲁棒性更好，也是近年来医学图像配准时普遍采用的相似性测度准则。基于归一化互信息的相似性测度定义为：

式中，H（X）是图像X的香农熵。因此，我们定义外力与参考图像和目标图像的相似性函数S（u）成比例。为了估计目标函数中的小波系数，需要对模拟平衡状态的能量函数进行极小化，因此，在这里，定义外力是归一化互信息的倒数：

2.3 形变配准算法

至此，得到了用于构成Navier偏微分方程来解决弹性形变问题的内力和外力，其中形变域用正交小波基函数来表示，内力和外力的平衡最终形成了关于形变域u的能量函数。因此，最终的能量函数可以表达为：

式中，w是加权常数，实验中，采用w为常数1。因此，形变域是关于小波系数的函数，即能量函数也是关于小波系数的函数。在配准过程中，代表全局形变特征的的大尺度小波系数首先被估计得到，然后估计对应局部形变特征的小尺度小波系数。为了减少需要被估计的小波参数的数量，根据子带和尺度的不同我们把小波系数进行分类，首先估计大尺度的小波系数，然后顺序估计同一尺度上的其他子带系数，最后估计小尺度上系数。从另一个角度看，两幅图像之间的形变由不同尺度和不同子带的小波系数来描述。我们采用递进策略寻找这些小波系数，先估计大尺度上的小波系数，然后再估计小尺度上的小波系数。图（1）是求解二维形变域时估计小波系数的顺序。

形变域经过小波变换以后，从理论上讲，需要估计的小波系数和待配准目标图像的体素个数是相同的。假设三维的图像大小是N=（Nx，Ny，Nz）=（512，512，512），优化算法一次需要估计的小波系数个数是512×512×512=134 217 128。如果估计全部的小波系数，这样可以精确恢复形变域，但是需要耗费大量的时间，这对于需要实时实现的自适应放疗系统是不能接受的。考虑到大多数的器官和组织产生的形变具有近似光滑性，我们用正交小波基函数表达形变时，只需要用一些低尺度的小波系数即可以较精确地恢复形变域，这对于配准的实时实现具有重要的意义。采用正交小波变换7层分解，仅使用第6和第7尺度上的小波系数来估计形变域，实验表明，配准结果具有较高的精度，而数据运算量却大大减少，对于N=（Nx，Ny，Nz）=（512，512，512）图像大小，仅需要优化算法总共估计8×8×8=512个参数，同时，子带和低尺度上有较多的零系数，可以有效提高配准速度。

图1 二维图示小波系数的估计次序，次序是：1-2-3-4-22-33-44Fig.1 Example for the progressive order in which registration parameters were estimated for 2D，from 1-2-3-4-22-33-44

为了简化在极小化能量函数时的复杂度，我们采用的正交小波基函数应当满足三次正交性质，即这些小波基函数自身相互正交，同时分别与自身的一阶导数和二阶导数相互正交。满足这一性质可以使极小化能量函数的内力部分得到运算的简化［13］。采用接近三次正交的不完全正交三阶样条小波来进行表示形变域。采用Marquardt-Levenberg（M-L）优化算法对小波系数进行估计［13］。详细的配准过程采用如下的步骤：

步骤1：输入参考图像R和目标图像T，去掉图像中的床位投影。

步骤2：把形变域三维小波分解后，把小波系数的初始值置为0。

步骤3：使用M-L优化算法递归的极小化式（8），寻找合适的小波系数，以小波系数的尺度和子带作为循环参数进行递归。具体方法采用步骤4和步骤5。

步骤4：根据小波分解塔式结构，首先估计最粗尺度J上的小波系数，从代表概貌的子带1开始，然后估计子带2，一直到子带8；接下来估计J-1尺度上的小波系数，从子带2开始，然后估计子带3，一直到子带8。

步骤5：重复步骤4，一直到金字塔剩余的各个尺度，估计每一个尺度上子带2到子带8的小波系数，根据所需的精度控制在某个尺度上收敛。

步骤6：根据估计到的小波系数得到参考图像与目标图像之间的形变域。

3 实验及结果分析

实验数据来自山东省肿瘤医院Varian’s CBCT系统，选取20组计划CT（planning CT，PCT）和治疗时现场采集的对应CBCT病例数据进行配准实验。CBCT与普通CT成像区别是一个角度直接形成三维的锥形束投影体数据，一周的照射即可快速形成三维的断层成像，但是这种成像的特点会导致获得的CBCT图像总是含有伪影，同时含有较大的噪声，如图2所示。

图2 CT图像。（a）计划CT；（b）CBCTFig.2 CT images.（a）planning CT；（b）cone beam CT image

为了定量验证算法的有效性，CBCT图像可以认为是由普通CT经过灰度均衡变换（Gray_Trans）并加入乘性斑点噪声（Noise）形成的，即：

式中，ImageCBCT表示变换后形成的模拟CBCT图像，ImageCT表示原始的CT图像，经过式（9）的变换后形成CBCT，与原始的CT形成具有金标准性质的配准对。对配准对中的CBCT图像进行定量的变形，变形方法采用ITK软件包中的B样条矢量变形方法。考虑到实际病例中形变的多样性，采用了3种由小到大的定量变形，如图3所示。

同时，临床治疗时，CBCT与CBCT的配准可以对肿瘤的变化或者病变器官进行跟踪，PCT与PCT的配准可以比较治疗前后的变化或者不同病例之间病变的差别。因此，对PCT和CBCT也进行了同模态的形变配准，同样，对PCT或者CBCT分别进行图3中由小到大3种变形，然后再与自身进行配准。

通过计算已知的形变域与形变配准算法获得的形变域差值进行定量分析配准算法的性能，定义

式中，N是像素个数，Wi是由形变配准算法计算出的第i个形变矢量的峰值，而Ki是第i个已知形变矢量的峰值。由于基于B样条的自由形变配准方法被学者认为在配准医学图像时具有优势［1］，因此，本研究的方法与基于B样条的自由形变配准方法进行了比较。表1是分别使用本算法与基于B样条的自由形变配准算法配准后计算出的形变域差值的结果。为了增加实验结果的可靠性，每个结果均取20组配准数据对的均值。

图3 3种已知的形变域。（a）～（c）从左到右的形变域分别由小到大，分别记为D1，D2，D3；（d）形变域的峰值Fig.3 Three level of deformation fields.From（a）～（c），the deformation is from small to large and we index them using D1，D2，D3；（d）meta-image magnitude for the deformation

表1 基于小波基函数和B样条配准方法得到的形变域差值比较Tab.1 Deformation difference comparison between the wavelet based and B-spline based methods

由表1可以看出，由本方法获得的形变域差总小于基于B样条的自由形变配准方法得到的形变域差，这表明本方法配准时可以获得更高的精度，同时具有稳定的鲁棒性，可以更准确地恢复形变域。同时，研究发现，对于较小的形变1，例如图3中形变1，两种方法均可准确恢复，但在具有较大形变时，例如图3中形变3，基于B样条的自由形变配准方法恢复性能较差，这种情况下，如果首先进行全局的仿射变换，然后再进行基于B样条的配准，会更好的恢复形变，达到和基于小波基函数配准方法接近的性能。所提出的方法，即使在具有较大形变的情况下，也可以准确快速地恢复形变域。

最后，采集了3组计划CT和日常治疗时相应的CBCT图像，包括形变较大的胸部，前列腺图像和形变较小的头颈部图像，并把它们分别用基于本文的配准方法和基于B样条的自由形变配准方法进行了配准。由于这种情况下没有金标准，采用了配准后的棋盘图像来对比结果，并且临床医师和物理师对结果进行了主观评价。

图4是计划CT与CBCT配准实验结果比较图，由图4可以看出（可从计划CT与CBCT棋盘图像衔接的位置观测），基于B样条的自由形变配准方法在恢复形变性能上要比本文方法弱，特别是对于形变较大和较多的部位效果比较明显。由于胸部和前列腺部位的形变相比头颈部位要大，并且全局形变和局部形变特征均存在，本研究提出的方法更能表现出可以由粗到精的恢复全局和局部形变的良好性能。对于头部图像，形变本身较小，两种方法均可以较准确的恢复形变，但是观测棋盘衔接处，本研究提出的方法在恢复头颈部局部形变时同样表现了更好的性能。实验结果经过几位临床医师和物理师的评价和分析，认为本方法在恢复形变上，具有更好的性能，可以指导临床的应用。

在速度上，对于较小的形变时两种方法均可准确快速恢复形变，且平均速度近似相同，对于实验中采用的512×512×64的数据，平均时间均在4min左右。对于较大的形变，基于B样条的配准方法平均时间在8min左右，且在实验数据中发生误配的概率在50％左右，而基于小波的配准方法平均时间在6min左右，且在实验数据中均可以准确配准。因此，相对于基于B样条的自由形变配准方法，本方法具有速度和鲁棒性上的优势。

4 临床应用

计划CT与CBCT的配准最直接的应用就是在IGRT系统放射治疗时，配准结果的形变域参数指导治疗靶区的射线进行自适应调整，使靶区的剂量摄取合理，而对周围正常组织器官达到合理限度的伤害，这需要把算法移植入机器，然后控制机器运动，临床实际应用仍在讨论当中。另一个方面的重要应用是计划CT与CBCT的准确配准可以自动准确提取CBCT的组织器官，代替耗费大量物理师时间的手动分割。由于放射治疗时，计划CT中的器官在治疗之前由物理师进行人工分割，为了在治疗时或者是在分次治疗中计算剂量在靶区和正常器官上的分布，以便于剂量更加合理的分布，可以利用形变配准算法对治疗时或者是分次治疗中的CBCT图像进行器官和靶区的自动分割。

图4 胸部病例，前列腺病例，头颈部病例计划CT与CBCT配准后结果棋盘图像（在每组棋盘图像中，上行是本文提出的方法的配准结果，下行是用基于B样条的自由形变配准方法得到的配准结果）。（a）胸部配准结果；（b）前列腺部位配准结果；（c）头颈部位病例配准结果Fig.4 Checkboards after planning CT and CBCT deformable registration（in each group，the upper row is with wavelet based method，the bottom row is with FFD based method）.（a）registration result for chest images；（b）registration result for prostate images；（c）registration result for head and neck（H＆N）

利用配准算法进行器官自动分割的具体方法是：治疗前，物理师在制定放疗计划时，对CT图像上器官和靶区的轮廓进行勾勒。由于计划 CT与CBCT可以准确配准，这些在计划CT上的器官和靶区根据计划CT与CBCT之间的形变域变化，把计划CT的器官轮廓重新映射到CBCT数据上。对临床20个病例的CBCT图像中肝脏数据进行了自动分割，下图是随机抽取了其中3个病例图像中的8片数据，进行肝脏的分割显示。

如图5所示，经过山东省肿瘤医院多位医师评价，由配准算法自动分割得到的CBCT图像肝脏部分与人工分割得到的CBCT肝脏部分基本相同，经过计算，平均重叠区域在97％以上，可以用于制定计划系统。这主要得益于配准算法可以精确恢复形变，把在计划CT上的人工分割轮廓经过形变域映射到治疗中的CBCT图像上，从而由原先的计划CT肝脏边界轮廓转化为CBCT上的肝脏边界轮廓。其它器官以及靶区的分割类似，这些器官和靶区的自动分割对于在治疗过程中快速合理计算放射剂量的分布具有重要的作用。

图5 由配准算法自动分割CBCT病例图像中肝脏结果（每个病例随机抽取了其中8片的分割结果。（a）～（c）分别为其中3个病例Fig.5 The liver automatic segmentation result using deformable registration algorithm for CBCT images（8 slicers are selected randomly from each case）.（a）～（c）are the results of three patients，respectively

5 结论和讨论

本研究提出了改进的基于CBCT的图像引导放射系统中计划CT与日常治疗中的CBCT图像的配准算法，采用基于正交小波基函数来表述形变域，这些不同尺度的小波基函数分别表示了不同的形变特征，粗尺度的小波系数代表了全局的形变，精细尺度的小波系数代表了局部的形变，不同子带的小波系数也表征特有的形变细节。同时利用Navier偏微分方程来设计形变能量函数，方程的自由膨胀部分作为外力，内力部分利用互信息作为距离函数来平衡方程的外力。表示形变域的小波部分通过优化算法由粗尺度到精细尺度逐步估计得到。实验表明，相对于传统的基于B样条的自由形变配准方法，本文提出的配准方法具有精度和速度上的优势。同时，本方法可用于其它用途和其它种类的医学图像配准。

［1］Holden M.A review of geometric transformations for nonrigid body registration［J］.IEEE Transactions on Medical Imaging，2008，27：111 -128.

［2］Sorcini B，Tilikidis A.Clinical application of image guided radiotherapy，IGRT（on the varian OBI platform）［J］.Cancer Radiothérapie，2006，10：252 -257.

［3］Jean P，Ali H，Chen J，et al.Low-dose megavoltage conebeam CT for radiation therapy［J］.International Journal of Radiation Oncology，Biology，Physics，2005，61：552 -560.

［4］Cheng B，Yang Yong，Li Fang，et al.Performance characteristics and quality assurance aspects of kilovoltage conebeam CT on medical linear accelerator［J］.Medical dosimetry：official journal of the American Association of Medical Dosimetrists，2007，32：80-85.

［5］Paquin D，Levy D，Xing Lei.Multiscale registration of planning CT and daily cone beam CT images for adaptive radiation therapy［J］.Medical Physics，2009，36：4-11.

［6］Brock KK，Dawson LA，Sharpe MB，et al.Feasibility of a novel deformable image registration technique to facilitate classification，targeting，and monitoring of tumor and normal tissue［J］.International Journal of Radiation Oncology，Biology，Physics，2006，64：1245-1254.

［7］Nithiananthan S，Brock KK，Daly MJ，etal.Demons deformable registration for CBCT-guided procedures in the head and neck：convergence and accuracy［J］.Medical Physics，2009，36：4755-4764.

［8］Noe K，Senneville BD，Vindelev U，et al.Acceleration and validation of optical flow based deformable registration for image guided radiotherapy［J］.Acta Oncologica，2008，47：1286 -1293.

［9］Bajcsy R，Kovacic S.Multiresolution elastic matching［J］.Computer Vision，Graphics，and Image Processing，1989，46：1-12.

［10］Pluim JPW，Maintz JBA，Viergever MA.Mutual information matching in multiresolution contexts［J］.Image and Vision Computing，2001，19：45-52.

［11］Mallat SG.A wavelet tour of signal processing［M］.San Diego：Academic Press，1998.

［12］彭玉华.小波变换与工程应用［M］北京：科学出版社，1999.

［13］Gefen S，Tretiak O，Nissanov J.Elastic 3-D alignment of rat brain histological images［J］.IEEE Transactions on Medical Imaging，2003，22：1480-1489.

［14］Pluim JPW，Maintz JBA，Viergever MA.Mutual-informationbased registration of medical images：a survey［J］.IEEE Transactions on Medical Imaging，2003，22：986-1004.