PET/CT成像呼吸运动B样条校正

潘李鹏，贺建峰，封 硕，崔 锐，马 磊，相 艳，易三莉，张 俊

昆明理工大学 信息工程与自动化学院，昆明 650500

1 引言

在PET/CT扫描成像过程中，病人的各种运动往往会造成伪影的产生，尤其是呼吸心跳等运动产的伪影，会严重影响PET/CT图像的质量。它将会干扰医生的诊断、制定治疗计划和治疗效果评估。针对呼吸运动校正，目前常用方法之一有呼吸门控方法[1]。它把呼吸运动划分成不同的时间相位（时相），运用特定外部设备跟踪呼吸运动的位置状态，再将不同呼吸周期中具有相同时相的数据叠加后再重建。这样各时相的呼吸伪影便大大减轻。一般门控方法需要额外的硬件（如红外监测装置、电磁装置、弹性传感带等）[2-4]来监测呼吸运动的位置状态。

呼吸门控除了用上面提到的硬件跟踪方法外，也采用基于扫描数据的数据驱动（data-driven）门控方法[5-9]。这种方法不需要硬件跟踪设备，而是基于已采集的数据通过算法估算出具有相同位置状态的时相。但是目前现有的数据驱动方法基本上是在图像空间内进行操作来进行门控校正的，即在已重建的图像序列上，运用不同的算法分析处理图像序列估算出运动时相。由于是在图像序列上估算运动参数，这类方法需要重建大量图像序列，运算量非常大，而且其精度仍需进一步改善[5]。例如，数据驱动方法中常用的质点中心方法（Centre of Mass），根据每幅图像序列不同状态的质心点（例如图像灰度值质心点）差异，估算出呼吸时相进行门控[7，9]。这种方法在呼吸运动速度过快，超出扫描范围或扫描数据采集不足的情况下会影响估算结果的精确性。

利用硬件或数据驱动的门控方法虽然能够降低运动成像伪影，但会导致图像的信噪比大大降低。这是因为它只提取了扫描总数据的10%～20%进行门控处理，而大部分有效数据被白白浪费了[10]。为了提高信噪比，只有利用更多的扫描数据进行门控。这可以将扫描时间加大到正常时间的5～10倍；或者成比例地增加放射性示踪剂的含量，但势必会使患者接受更多的辐射剂量，是不可取的。

为此，研究人员提出用图像配准的方法利用所有扫描采集数据进行图像校正。这种方法是先估算出已门控时相配准参数，然后再根据配准参数将所有已门控时相数据全部配准到参考相位以提高图像校正信噪比[11-14]。例如光流法（optical flow algorithm），在已重建所有图像序列的基础上，运用算法估算出每两幅图像之间每个像素的运动状态来确定其配准运动参数，然后进行图像配准[15-16]。它需要重建图像并对所有序列图像里的像素进行处理，其算法的复杂性和运算量不应忽视[5]。

为此，本文提出一种基于运动特征的呼吸运动校正方法，通过时间帧对应提取与PET图像相对应的CT图像运动特征，来配准校正呼吸运动图像。

2 方法

PET图像因为其固有的缺陷和不足使得图像质量较差噪声也较高，这就造成了PET图像的结合边缘比较模糊。这样的图像中所提取的运动信息的可靠度较低，难以满足临床运用对精度的需求。在PET/CT的一次扫描中就可以同时获取PET和CT的图像，因此本文提出通过对CT图像进行运动特征信息提取，然后使用这些信息对PET图像进行呼吸运动校正。本文所使用的运动特征信息是一系列的控制点，通过这些控制点对PET图像进行B样条变换就可以达到校正目的。

2.1 基于B样条的FFD配准方法

在本文中进行图像配准的目的是通过对不同时间帧状态下不同呼吸相的肺部CT图像进行对比和校正，得出用于校正PET图像所需的运动特征信息。而这种配准的实质是寻找一最优变换T：(x，y，z)→(x′，y′，z′)，此式映射在一个动态图像序列I(x，y，z，t)在时刻t上任意一 点对应到时刻t0图像I(x′，y′，z′，t0) 上点的关系。分析人体肺部呼吸运动特性可知，肺部呼吸运动主要表现形式为非刚性形变运动。

基于B样条[17-18]的配准已被用于心脏图像运动分析和乳房MRI图像配准等应用上，证明其可以被用于人体非刚性组织的运动配准。这里选择基于B样条的FFD（Free form Deformation）方式来对肺部影像进行处理。本文首次提出使用基于B样条的配准来提取CT图像的运动信息，用于PET图像的校正。FFD是一种用于对立体模型进行任意自由变形的一种技术，其可以对任意平面、角度、参数的表面进行变形。这种变换可以是全局性的也可以是局部性的。FFD的基本思想是使用由一系列的控制所组成的网格来控制定义一个物体的形变。

为了定义一个基于FFD的样条，定义图像空间为Ω={(x，y)|0≤x＜X，o≤y＜Y}。Φ表示一个在均匀空间σ内nx×ny网格控制点Φi，j。然后FFD可以被写成熟悉的一维三次B样条张量积：

其中i=(x/nx)-1，j=(y/ny)-1，u=x/nx-(x/n)，v=y/ny-(y/n)，并且Bl表示第l个B样条基函数：

相对于 thin-plate样条[19]或者 elastic-body 样条[20]，B样条拥有更好的局部控制特性，能够有效地计算大量的控制点。

控制点集Φ作为B样条FFD和非刚性变换角度的参数。非刚性变换角度可以通过控制点Φ构成网格的分辨率来表示。当控制点网格间距比较大也就是网格比较稀疏时，可用于全局非刚性变换，当网格间距比较小时也就是网格比较密集时可用于局部变换。同时，网格的间距越小，计算的复杂度就越高。因此在选择网格密度时，需要权衡多方因素，在计算时间合理的前提下找到变换效果最好的网格间距。

为了能够更好地对呼吸运动进行校正，引入塔形多层网络结构校正的方法。首先使用分辨率较低的控制点网格对图像进行整体校正，随后使用分辨率高的网格对图像进行校正：初始数据第一次经过分辨率32×32的控制点网格校正后得到第一次中间图像数据，随后对该数据用分辨率64×64的控制点网格进行校正并得到第二个中间图像数据，最后用分辨率128×128的控制点网格进行最后一次校正，最终得到所需要的校正图像。

2.2 基于CT图像的运动特征提取

首先分析呼吸运动特性，在呼吸运动期间，肺部的运动变化呈现出周期性的特点。也就是说在一个周期内，呼吸运动相位均有差异，而在周期间这一运动过程是不断重复的过程，类似正弦曲线的波形如图1所示。其中横轴t表示呼吸周期，每一个重复的弧形曲线表示一个呼吸周期，曲线上升段表示吸气状态此时肺部膨胀，曲线下降端表示呼吸状态，此时肺部缩小。在每个周期的起始和末尾肺部都达到自然态，而曲线在一个呼吸周期的中点也就是达到波峰的时候肺部体积达到最大，在此处图像面积也最大。

图1 呼吸运动曲线图

为了能更好达到校正的目的，正确有效地提取运动信息至关重要。在运用CT图像提取运动信息对PET图像进行校正之前，要解决PET图像和CT图像之前的时间周期匹配问题。目前CT的成像速度要大大快于PET图像。这就造成了CT图像每一帧的图像要多于PET图像，如图2所示。本文使用时间轴对称采样的方法来解决这一问题。提取CT图像中时间片时刻和PET图像完全相同的图像，这些图像的位置和几何结构与同时刻的PET图像完全相同。事实上，这一过程可以视为对CT图像周期数据的一个稀疏化过程。对这些CT图像进行运动特征提取后，就可用于PET图像呼吸运动校正。

图2 经过周期同步后的CT和PET图像

运动特征提取是找寻PET连续图像中各个不同时间帧的对应图像之间的差异，用于对这些图像做一些补偿处理使得这些图像都能够向基准图像靠拢，这样就可以去除由于运动所造成图像伪影和图像质量降级。通过寻找图像间差异方式实现运动特征提取。

选取基于B样条的FFD配准方法来获取呼吸特性，那么这种呼吸特性的表现形式为描述呼吸周期内各个呼吸相与基准相差异的一组控制点Φi。这一过程可由以下公式表示：

其中MR表示一个B样条配准的过程，fi表示在一个呼吸周期内的某一帧，ft表示用于配准的基准帧。根据公式（2）对周期内每一帧图像均做配准特征提取处理后得：

最终得到Φ={Φ1，Φ2，Φ3，…}为呼吸运动一周期内的运动特性。

为了获取运动特征信息，首先需要选定用于差异特征提取的基准帧图像，也就是基准呼吸相。所有的图像都将把基准帧图像做为参考，做配准来提取运动特征（控制点集Φ）。本文使用NCAT生成的肺部体模[21]，作为研究对象，设置一个呼吸周期有11个呼吸相，也就是说用11个图像序列来模拟一个呼吸过程。

2.3 选取基准图像帧

先确定一个基准帧，在配准过程中这一基准帧可以是动态变化的。用B样条获取特征，当获取到运动信息后通过这些控制点来对图像做变换使其向基准靠拢。如图3所示，本文选取第6帧作为初始基准帧，经过一次配准后图像会作为新的基准帧与相邻的帧来做配准。其中SF表示初始基准帧，在获得第5帧的运动特征信息后，对第5帧图像做B样条变换得到第5帧图像的动态基准帧DSF5，用于和第4帧配准，其他帧也可以以此类推，最终得到所有的运动特征信息（控制点集Φ）。

2.4 校正效果定量分析方法

在运动校正后需要对校正结果进行评估，可以通过视觉观察图像质量评价校正效果。但是这样往往难以发现不同校正方式和参数在校正结果上细微的差异，而且人眼有的时候难以做到准确判断图像质量。本文提出使用RMCR（Respiration Motion Correction Rate）呼吸运动校正率来对校正结果进行定量分析。

图3 动态基准方法

在校正过程中，所有的图像切片都是与基准图像做配准和变换的，如果对校正后图像和基准图像做差就可以得到一个配准后图像与基准图像的面积差DAIT，如公式（3）所示。其中IMST表示基准图像，IMA表示校正后图像。同理，对未校正图像和基准图像做面积差可得配准前图像差DBIT，如公式（4）所示，其中IMB表示配准前图像的面积。

定义RMCR未校正前图像差和校正后图像差的差与未校正前图像差的比值，它描述了在校正的过程中有多少差异被校正，从理想状态来说，校正后的图像应该与基准图像相同或者非常接近那么RMCR的数值就应该较大，证明大量的差异已经被校正。RMCR如公式（5）所示。因为在实际应用中要校正的图像都是多帧组成，而每个时间帧中又包含很多图像切片，因此可以取所有切片的平均RMCR作为对整个校正过程的效果评估，如公式（6）所示，其中M表示一个呼吸周期内的帧数，而N表示一个时间帧内所包含的图像切片数。

3 实验结果

实验校正的具体步骤如图4所示，首先对CT图像进行同步处理。由于CT图像与PET图像存在采集时间上的差异，造成两种图像存在周期性的差异。由于CT采集图像速度快于PET图像，在一个相同时长的周期中CT图像可以分割出的相位多于PET图像。那么根据PET图像的相位时间点来匹配CT图像，去除CT图像多余PET图像中的相位图像以实现周期同步的目的。第二步，对处理后的CT图像进行运动特征提取。在这一过程中，首先对同一周期不同时间相位的CT图像进行和基准帧的配准，可以得到每一个相位的CT图像要变换成基准图像所需的控制点集Φ。第三步，就是对所有的PET图像进行B样条变换，其中每一个相位的PET图像都使用与其相对应的CT图像的运动信息（控制点信息Φ）进行变换。经过这一系列的变换后，可以得到所需的经过运动校正的PET图像。本文设计了两种实验：几何位移形变测试和像素体模测试。

图4 基于运动信息呼吸校正方法流程图

3.1 几何位移形变测试

肺部的运动是一个非常复杂的过程，简单的使用线性位移和规则形变难以模拟肺部的复杂运动。因此，本文设计了一个结合了线性位移和形变位移的运动模型来模拟肺部运动，如图5所示。128像素×128像素的二值化图像中绘制一个50×50的正方形，然后加入膨胀缩小和线性位移信息。图像融合后的一些重要信息都被膨胀信息所遮挡使得图像效果下降。而本文的校正方法使用网格尺寸为32×32的网格进行校正，如图6所示。校正后图像遮挡的问题得到了解决，也没有出现由于位移而造成的运动模糊。但是依然存在图像角落模糊的情况，这一问题可以通过提高控制点网格分辨率的方法改善。

图5 位移形变模型

图6 位移形变模型

3.2 基于NCAT模拟数据的运动校正实验结果

简单的几何图像测试还不能说明本文方法的实际执行情况，因此采用更为接近真实的数据，基于NCAT体模呼吸模型对提出的方法进行验证。NCAT可以被当作人体模型应用于PET/CT扫描仿真平台的GATE软件包中[22]。

NCAT是一个人体全身的模拟，由于其数据源自于真实的人体数据，所以包含较多的细节使得处理过程会很长。而本文的主要研究对象为人体的肺部运功，这里运用NCAT生成一个包含运动信息的呼吸周期图像。在这一周期中一共有11帧图像（其中1～5帧为吸气7～11为呼气），每一帧包含着31个图像切片。每一帧图像的各个切片都包含着相邻帧图像对应切片图像的差异运动信息，就可以运用本文提出的方法（如图4所示）进行呼吸校正。

在获取可以用于模拟呼吸过程的运动体模信息后，就可以将该模型引入到基于GATE的PET仿真环境中，用于模拟PET扫描过程获取包含呼吸运动的PET仿真图像。图7所示为一个只包含人体胸腔的NCAT体模在经过GATE仿真并还原后的图像。

图7 含人体胸腔的NCAT在经过GATE仿真并还原后的图像

图8给出了校正结果，11帧图像的第15个切片在校正前和校正后的比较。可以发现所有的图像都向基准帧第6帧靠拢，这样也就达到了校正的目的。

图8 校正实验结果

从实验结果来看，校正后的图像和未校正图像有着明显的差异，所有图像的形状都接近本文所选定的初始基准第6帧图像（作为基准帧其本身不需要校正；由于校正的时候使用了差值的方法使得在对周围的像素起到了平均化的效果，使得图像出现了一些花纹，这些随着图像分辨率的提高这一问题就会减弱）。

3.3 结果讨论

根据三网格层迭代校正的实验数据，对其进行定量分析，用本文提出的RMCR方法如公式（6），本文的实验数据中M为11，N为31，计算这次实验的平均呼吸运动校正率。由于DAIT和DBIT在之前校正程序中已经作为参数存在，所以可以直接在进行校正的同时获取；可通过简单的迭代程序实现来计算最终整体的RMCR。

经过计算本实验的RMCR值为0.882，也就是说有10帧图像中的每一切片和基准帧图像切片的差异有88.2%被消除。实验图像为128像素×128像素，每个图像的面积约为150～600像素，校正后大约还剩余11.8%的差异也就是17.7～70.8像素的差异（面积大的差异像素多，面积小的差异像素少），放在整体128像素×128像素的图像中这个误差为0.108%～0.43%之间，效果还是好的。

可以发现第5帧图像在校正前和校正之后没有什么差别，这是因为第5帧图像和基准图像十分接近（其DBIT＜10，两幅图像的面积差小于10个像素），在误差允许范围之内，没有必要对其进行校正，对其校正可能会使图像质量下降。

迭代校正次数与图像质量的关系也十分密切，不是处理次数越多效果就越好，需要根据实际图像的不同来确定迭代处理的次数，如果迭代次数过多也会降低图像质量。本文使用的图像数据为128像素×128像素，对其进行3次迭代校正效果较好，如果迭代4次就会出现过度校正现象。

另外，在实验的过程中发现网格分辨率对校正图像质量有很大的影响，对于简单的几何图像，当网格分辨率较低的时候对图像质量影响有限，但是对于复杂的不规则图像做非刚性运动校正时，网格分辨率低会严重影响校正效果。将NCAT图像分别用2×2，4×4，8×8以及32×32，64×64，128×128这几种分辨率进行3次迭代，其效果差异十分明显，在低网格分辨率下的非刚性运动校正出现了畸变，而在高分辨率下效果比较理想。

4 结论

针对PET/CT成像中的呼吸运动影响图像质量问题，本文提出了基于CT图像提取呼吸运动特征的B样条校正方法。通过基于位移形变刚性运功实验和基于NCAT呼吸模型的非刚性运动校正实验，证明本文方法在控制点分辨率大小以及迭代次数合理的情况下，能够明显改善PET/CT呼吸运动图像质量。

[1]许全盛，袁克虹，于丽娟，等.PET/CT图像呼吸运动伪影校正研究进展[J].中国生物医学工程学报，2009，28（4）：573-580.

[2]Schaerer1 J，Fassi A，Riboldi M，et al.Multi-dimensional respiratory motion tracking from markerless optical surface imaging based on deformable mesh registration[J].Physics in Medicine and Biology，2012，57：357-373.

[3]Chang G，Chang T，Pan T，et al.Implementation of an automated respiratory amplitude gating technique for PET/CT：clinical evaluation[J].Journal of Nuclear Medicine，2010，51（1）：16-24.

[4]Chamberland M，Wassenaar R，Spencer B，et al.Performance evaluation of real-time motion tracking using positron emission fiducial markers[J].Medical Physics，2011，38（2）：810-819.

[5]Kesnera A，Kuntner C.A new fast and fully automated software based algorithm for extracting respiratory signal from raw PET data and its comparison to other methods[J].Medical Physics，2010，37（10）：5550-5559.

[6]Kesner A，Bundschuh R.Respiratory gated PET derived in a fully automated manner from raw PET data[J].IEEE Transactions on Nuclear Science，2009，56（3）：677-686.

[7]Buther F，Dawood M，Stegger L，et al.List mode-driven cardiac and respiratory gating in PET[J].The Journal of Nuclear Medicine，2009，50（5）：674-681.

[8]He J，G.Keefe O，Gong S，et al.A novel method for respiratory motion gated with geometrical sensitivity of the scanner in 3D PET[J].IEEE Transactions on Nuclear Science，2008，55（5）：2557-2565.

[9]Schleyer P，Doherty M，Marsden P.Extension of a datadriven gating technique to 3D，whole body PET studies[J].Physics in Medicine and Biology，2011，56：3953-3965.

[10]Nehmeh S，Erdi Y.Respiratory motion in positron emission tomography/computed tomography：a review[J].Seminars in Nuclear Medicine，2008，38（3）：167-176.

[11]Hahn D，Daum V，Hornegger J.Automatic parameter selection for multimodal image registration[J].IEEE Transactions on Medical Imaging，2010，29（5）：1140-1155.

[12]Loeckx D，Slagmolen P，Maes F，et al.Nonrigid image registration using conditional mutual information[J].IEEE Transactions on Medical Imaging，2010，29（1）：19-29.

[13]Vandemeulebrouckea J，Rit S，Kybic J，et al.Spatiotemporal motion estimation for respiratory-correlated imaging of the lungs[J].Medical Physics，2011，38（1）：166-178.

[14]Huang X，Ren J，Guiraudon G，et al.Rapid dynamic image registration of the beating heart for diagnosis and surgical navigation[J].IEEE Transactionson MedicalImaging，2009，28（11）：1802-1814.

[15]Dawood M，Lang N，Jiang X，et al.Respiratory motion correction in 3-D PET data with advanced optical flow algorithms[J].IEEE Transactionson MedicalImaging，2008，27（8）：1164-1175.

[16]Dawood M，Brune C.A continuity equation based optical flow method for cardiac motion correction in 3D PET data[C]//Proceedings of the 5th International Workshop on Medical Imaging and Augmented Reality，2010，6326：88-97.

[17]Duijster A，De Backer S，Scheunders P.Multicomponent image restoration，an experimental study[C]//Kamel M.Proceedings of the 4th International Conference on Image Analysis and Recognition，MontrealCanada，August 22-24，2007，4633：58-68.

[18]Jiang M，Wang G，Skinner M W，et al.Blind deblurring of spiral CT images[J].IEEE Transactions on Medical Imaging，2003，22（7）：837-845.

[19]Khare A.A new method for deblurring and denoising of medical images using complex wavelet transform[C]//Proceedings of IEEE Engineering in Medicine and Biology 27th Annual Conference，2005，7：1897-1900.

[20]Herraiz J L，Samuel S，Vicente E，et al.Optimal and robust PET data sinogram restoration based on the response of thesystem[C]//ProceedingsofIEEE NuclearScience Symposium Conference，2006，6：3404-3407.

[21]Segars W P.Development and application of the new dynamic NURBS-based Cardiac-Torso（NCAT）phantom[D].The University of North Carolina at Chapel Hill，2001.

[22]Jan S，Santim G，Strul D，et al.GATE：a simulation toolkit for PET and SPECT[J].Physics in Medicine and Biology，2004，49（19）：4543-4561.