配准算法对 PCA 单幅投影肺部重建的影响

范忠银 王建军 朱青松 谢耀钦 邢 磊

(中国科学院深圳先进技术研究院 深圳 518055)

配准算法对 PCA 单幅投影肺部重建的影响

范忠银 王建军 朱青松 谢耀钦 邢 磊

(中国科学院深圳先进技术研究院 深圳 518055)

肺癌放射治疗中，肺部肿瘤位置实成像对于临床意义重大。在一种利用单 X 射线投影进行成像的实时肺部 3D 成像算法中，图像配准过程引入的不准确对于 PCA 模型构建以及重建过程有重大影响。文章分析了光流法、Demons 算法、水平集算法三种配准算法对重建效果的影响，并通过定性以及定量实验分析验证。结果表明，光流法配准在配准结果以及模型构建方面有较好的效果。

实时肺部 3D 成像；主成分分析；弹性配准

1 引 言

在现代肺癌的放射治疗中，治疗过程中对实时肺部肿瘤位置进行精确定位十分重要[1]。肺部肿瘤定位方法可以大致分为：直接肿瘤定位、基于呼吸参照定位以及基于植入标记物定位。现今的直接定位技术存在以下缺点：(1)方法需从前后方向获得图像且肿瘤边界清晰时才能有效运行；(2)需要经过预处理的荧光透视图像；(3)无法精确得到垂直于 X 射线成像平面的第三维肿瘤位置信息。基于呼吸参照物的定位方法不需要额外的放射剂量来进行定位，但是肿瘤位置和参照物内外部关系的微小变化对定位有很大的影响[2,3]。基于植入标记物的定位方法定位精度较高[4]，但这种方法需要侵入人体，且易引起气胸[5]。

要克服这些不足，需构建一个用于直接定位的精确肺部运动模型。Zhang 等[6]提出一种基于 主成分分析(Principal Component Analysis，PCA)的肺部运动模型，只需要较少的特征向量和 PCA 系数就可以精确表示肺部运动。类似的，Li 等[7,8]将序列相位与参考相位进行弹性配准得到偏移矢量场(Deformation Vector Field，DVF)，计算 DVF 本征模式，建立一个类似于Zhang 的 PCA 模型，重建一个 4D-CT。然后利用 PCA 模型对 4D-CT 中的参考相位图像进行变形，将单个的 DRRs(Digitally Reconstructed Radiographs)与 CBCT(Cone Beam Computed Tomography)投影匹配来估计原始 4D CBCT 投影数据集每个时间步长的本征模式系数(即主要的系数)。PCA 肺部运动模型隐式配准且效率高，使其能精确、高效率地得到一个动态肺部运动模型，且只需要少量信息。

要将运动模型和先验图像综合起来，最主要的挑战是得到一个真实的模型，且模型具有合适数量的参数，并可以简单地在任何成像物体上归纳得出。PCA 肺部运动模型基于先验图像(4DCT)之间的弹性配准得到 DVF，这给 PCA 分析引入误差的源头：图像配准过程引入的不准确及近似问题。

在本文中，我们将对比几种主流配准方法在PCA 构建模型以及重建效果的差别，并甄选出最合适的配准算法。

2 基于 PCA 模型的单幅投影肺部重建方法

2.1 PCA 模型

来自于仿真治疗过程的 4D-CT 数据作为先验知识，以此进行三维图像重建。在参考位相和其他位相之间，进行弹性图像配准，产生一系列的 DVF。这些 DVF 能够用 PCA 方法得到的一系列特征向量和系数来表示。通过改变 PCA 系数，能够得到新的 DVF。把新的 DVF 应用于参考图像就能得到新的立体图像。然后优化 PCA系数(运用构建出的模型)使计算出来的投影数据与测量得到的立体图像投影数据一致。将反转的DVF 应用于参考图像就可以得到待重建图像。

在 PCA 肺部运动模型中，DVF(相对于参考图像是一个空间和时间函数)可以通过样本均值向量和一系列具有最大特征值的特征向量的线性组合近似估计。如下：

其中，x(t)是一个参数化的 DVF，即一个时间和空间的函数； 表示不同时间 DVF 的均值；uk表示由 PCA 得到的特征向量，是空间函数；wx(t)表示 PCA 系数，是时间的函数。因为每一个空间像素元在每一个特征向量都有一个自由度。所以，每一个特征向量定义了空间中体素运动的一个方向，两个特征向量就是一个平面，三个特征向量构成肿瘤运动的三维空间。这种特性使得 PCA 运动模型能够准确把握肿瘤运动的轨迹。在本文的方法中，保留具有最大特征值的 3个 PCA 系数。

在得到一个参数化的 PCA 肺部运动模型之后，我们希望能够得到一系列最优化 PCA 系数，使得用新的 DVF 重建的立体图像投影数据与测量得到的 X 射线投影数据能够匹配。f0作为参考图像，f是重建图像，y是测量得到的投影图像，P是用来计算f的投影图像的投影矩阵。代价函数如下：

其中，

其中，U是一个矩阵，它的列是 PCA 的特征向量；w是一个包含了待优化的 PCA 系数的向量；x是参数化的 DVF；为p阶泛数。

对代价函数进行优化：

其中

每次迭代之后，由公式(4)得到更新的 PCA系数，从而更新 DVF，由三线性插值得到重建图像，对应地，必须与插值过程保持一致，从而得到正确的梯度值。

2.2 主流配准方法

在 PCA 模型中，参考相位和其他相位之间配准生成 DVF，这些 DVF 决定了建模的准确度。配准算法对于整个重建有着重大的影响。为研究图像配准方法对 PCA 建模以及重建效果的影响，选取三种经典的配准方法：光流法(Optical Flow Method，OFM)、Demons 法、水平集法(Levelset)进行比较研究。

其中，光流场方法起源于计算机视觉领域，主要用于运动物体检测，由于运动物体的光流场与图像配准的位移场具有一定的相似性，因此光流场的基本原理也可以用于图像配准。主流的光流场方法是基于 Horn 等[9]提出的变分法，该方法基于像素灰度值的局部梯度模型和整体平滑性假设。相对于 Horn 的全局化模型，Lucas 等[10,11]提出了局部化窗口模型，其思想是将每一像素速度矢量的求解放到其邻域中，由邻域内各像素的信息共同约束该点速度矢量。窗口化操作中，每一像素只受其周围邻域像素的影响，可以通过限定小的窗口来减少对图像造成的模糊。本课题研究中采用窗口化的 Lucas 模型进行图像配准。

Demons 算法配准是目前最为流行的配准方法。Demons 模型由 Thirion[12]于 1998 年提出并应用于图像配准领域。将源图像视为一个可变形网络，其形变自由度可以人工设定。图像像素对应于网格的格点，并且都被赋予一定的极性：内点或外点。Demons 根据格点的极性来控制方向，内点只能进入某个区域，外点只能离开某个区域。如果将源图像放入布满 Demons 的目标图像空间内，则网格就会在 Demons 的控制下发生形变，最终所有内点都将进入某一区域内，而外都将被退出到某一区域外。通过合理设置所有格点的极性，可以使得源图像自动发生形变，达到与目标图像相似程度最大。

水平集方法是 Sethian[13]在研究曲线以曲率相关的速度演化时提出来的，用于描述曲线或曲面的演化过程，是处理封闭界面在演化过程中几何拓扑变化的有力工具。其最大的优势在于它的稳定性以及拓扑无关性。但 水平集方法有一个缺点，即计算量太大。

2.3 算法流程

本文选取了 Demons 配准算法、水平集配准算法和 OFM 配准算法作为对比方法，来考察它们在 PCA 构建模型以及重建效果上的差别。

算法过程如下：

(1)数据准备：10 个相位的肺部 CT 三维数据；

(2)选取第一个相位作为参考图像，其他相位作为浮动图像，进行弹性图像配准；

(3)PCA 提取主成分；

(4)利用主成分进行迭代重建。

算法流程如图 1 所示。

3 实验与结果分析

3.1 实验数据

本文研究的图像数据为 4D-CT 数据。CT 扫描结束后，利用 4D 软件对得到的大量不同位置、不同呼吸时相的数据进行处理。以 10% 为时相间隔，把一个呼吸周期分为 10 个呼吸时相。分别命名为 0%(吸气末)、10%(吸气中)、20%、…、50%(呼气末)、…、90%(吸气末)，依次定义为 0 相位至 9 相位。

针对肺部的形变情况，我们截取了 CT 影像序列中肺腔的部分，这样可以减少不必要的数据量，加快算法的速度。本文使用的是连续 0～127层 CT 影像信息，数据大小为 256×256×128，体素分辨率 1 mm×1 mm×1.25 mm。投影大小为 300×400，探测器像素 1 mm×1 mm。

图 1 算法流程图Fig. 1 The flow chart of the proposed algorithm

图 2 是呼吸周期 10 个呼吸时相中间层(第 64层)的图像。在图像中显示为黑色区域的为肺组织，可看出肺轮廓的变化情况。而毫米级的组织运动变化对于放疗计划都将有重大的影响，所以需要采集实时的三维体数据以根据人体组织运动变化制定治疗计划。

3.2 弹性图像配准

选择呼吸周期 0 相位的图像作为参考图像，其他 9 个相位作为待配准图像。

图 2 呼吸周期 10 个呼吸时相中间层(第 64 层)的图像Fig. 2 The 64th slice of each phase in a respiratory circle

图 3 是待配准图像为呼吸周期 50% 时相，用三种配准方法得到的结果。与标准待配准图像对比，我们发现 OFM 配准的结果更加接近待配准图像，在图像边缘部分尤为明显。其中，从上到下依次为从冠状面、矢状面和横断截取中间层的图像。

表 1 为对配准进行定量分析的结果。其中，配准误差采用公式(6)计算。不同配准方法得到的位置偏移不同，Demons 方法对配准偏移估计过大，而 水平集方法对偏移估计较小，OFM 得到了一个较为适中的偏移，且 OFM 配准误差更小。

从图 3 定性的说明以及表 1 定量的计算，都可以看出：OFM 配准的结果较为准确。

图 3 待配准图像为呼吸周期 50% 时相，用三种配准方法得到的结果Fig. 3 The registration deformed image at 50% phase with three methods

表 1 三种配准方法的定量对比Table 1 The quantitative comparison of three methods

3.3 PCA 结果

将配准步骤得到的x、y、z方向上的偏移矢量进行 PCA。

表 2 是三种算法配准偏移矢量进行 PCA 之后x、y、z方向上各个特征向量的贡献率，最后一行为前三个特征向量的累计贡献率。可以看出，Demons 算法和 OFM 算法的结果 PCA 能量较为集中，前三个特征向量累计贡献率也较高，从而对后续的重建更加有利。

3.4 重建结果

图 4 为三种方法结果对比情况。其中，从上到下依次为冠状面、矢状面和横断截取中间层的图像。从图像的轮廓以及灰度级初略查看，三种方法的重建结果差别不大。

采用定量的方法加以考察，其中重建结果的配准误差分别为：Demons 为 0.0145727，水平集为 0.0130284，OFM 为 0.0129043。可以得到，OFM 方法得到的重建图像与标准图像误差最小。

表 2 PCA 累积贡献率Table 2 The accumulative score of the principal component analysis

4 讨 论

本文介绍的单个 X 射线图像重建肺部三维影像算法在临床应用具有广阔的潜力。本项研究中肿瘤定位误差来自两个不同的过程：(1)PCA 肺部运动模型；(2)投影和参考图像之间的 2D/3D配准。一个可靠的弹性配准算法，特别是能够精确地模拟较大位移运动的配准算法，从而提高PCA 肺部运动模型的精确度。因此，我们采用Demons、水平集和 OFM 配准算法来对 PCA 肺部运动模型的精确度和重建效果进行分析。综合以上分析，在三种方法中，光流法配准结果对于PCA 模型构建和重建效果较好，而 水平集算法在配准结果以及重建方面性能都较差。本文选取的配准方法有限，下一步将继续扩大配准算法的范围来对模型构建和重建进行对比。

在治疗过程中，病人可能在分次间发生解剖变化，如果肺部运动模型是基于病人模拟仿真获得的 4D-CT 而建立，那么结果将会很不理想。在这种情况下，病人摆位过程采集的 4D CBCT建立的 PCA 肺部运动模型能克服这个问题。另一个潜在的问题是当肿瘤的位置靠近心脏时，肿瘤的运动不仅会受到呼吸运动影响，同时也会受到无法用 PCA 运动模型表达的心脏运动影响。这是 4D-CT 的一个局限，4D-CT 与呼吸同步而非心脏运动。在研究中，我们假设测量的锥束投影和待重建图像的仿真投影之间的灰度值有线性关系，但是未考虑到不同图像模态直接的差异，用于先验的 4D-CT 和术中重建的 CBCT 的 X 射线能量值和光谱值均不一样。综合这些意见，下一步将采用病人摆位过程采集的 4D CBCT 作为先验图像来改进算法。

[1] Keall PJ, Mageras GS, Balter JM, et al. The management of respiratory motion in radiation oncology report of AAPM Task Group 76 [J]. Medical Physics, 2006, 33(10): 3874-3900.

[2] Hoisak JD, Sixel KE, Tirona R, et al. Correlation of lung tumor motion with external surrogate indicators of respiration [J]. International Journal of Radiation Oncology, Biology, Physics, 2004, 60(4): 1298-1306.

[3] Tsunashima Y, Sakae T, Shioyama Y, et al. Correlation between the respiratory waveform measured using a respiratory sensor and 3D tumor motion in gated radiotherapy [J]. International Journal of Radiation Oncology, Biology, Physics, 2004, 60(3): 951-958.

[4] Harada T, Shirato H, Ogura S, et al. Real-time tumor-tracking radiation therapy for lung carcinoma by the aid of insertion of a gold marker using bronchofiberscopy [J]. Cancer, 2002, 95(8): 1720-1727.

[5] Laurent F, Latrabe V, Vergier B, et al. CT-guided transthoracic needle biopsy of pulmonary nodules smaller than 20 mm: results with an automated 20-gauge coaxial cutting needle [J]. Clinical Radiology, 2000, 55(4): 281-287.

[6] Zhang QH, Pevsner A, Hertanto A, et al. A patientspecific respiratory model of anatomical motion for radiation treatment planning [J]. Medical Physics, 2007, 34(12): 4772-4781.

[7] Li R, Lewis JH, Jia X, et al. On a PCA-based lung motion model [J]. Physics in Medicine and Biology, 2010, 56(18): 6009-6030.

[8] Li R, Jia X, Lewis JH, et al. Real-Time volumetric image reconstruction and 3D tumor localization based on a single X-ray projection image for lung cancer radiotherapy [J]. Medical Physics, 2010, 37(6): 2822-2826.

[9] Horn BKP, Schunck BG. Determining optical flow [J]. Artificial Intelligence, 1981, 17: 185-203.

[10] Lucas B, Kanade T. An iterative image registration technique with an application to stereo vision [C] // Proceedings of the International Joint Conference on Artificial Intelligence, 1981: 674-679.

[11] Lucas BD. Generalized image matching by the method of differences [D]. Pittsburgh: Carnegie Mellon University, 1984.

[12] Thirion JP. Image matching as a diffusion process: an analogy with Maxwell’s Demons [J]. Medical Image Analysis, 1998, 2(3): 243-260.

[13] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulation [J]. Journal of Computational Physics, 1988, 79(1): 12-49.

Deformable Image Registration Effect on PCA-based Reconstruction with Single X-ray Imaging

FAN Zhongyin WANG Jianjun ZHU Qingsong XIE Yaoqin XING Lei

(Shenzhen Institutes of Advanced Technology,Chinese Academy of Sciences,Shenzhen518055,China)

In modern lung cancer radiotherapy, it is important to have a precise knowledge of the real-time lung tumor position during the treatment delivery. For a real-time 3D lung imaging algorithm from a single X-ray projection image, the inaccuracies contributed by the image registration process affects much on the PCA modelling and construction process. We utilize 3 deformable image registration algorithms: Optical Flow method, Demons method and Levelset method to evaluate the effect. By making quantitative analysis and qualitative analysis, we get the conclusion: Optical Flow method works much better in registration and PCA modelling.

real-time 3D tumor localization; principal component analysis; deformable image registration

R 812

A

2014-04-22

： 2015-03-15

广东省引进科研创新团队项目(2011S013)

范忠银，硕士研究生，研究方向为医学图像处理、医学影像系统；王建军，硕士研究生，研究方向为计算机控制与三维可视化；朱青松，硕士，研究方向为计算机视觉、医疗机器人视觉；谢耀钦(通讯作者)，博士，研究员，研究方向为图像引导放射治疗、医学影像处理和分析，E-mail：yq.xie@siat.ac.cn；邢磊，博士，教授，研究方向为肿瘤放射物理。