2. 中国科学院大学 北京 100049;
3. 上海联影医疗科技有限公司 上海 201807
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Shanghai United Imaging Healthcare Co., Ltd., Shanghai 201807, China
正电子发射断层成像 (Position emission tomography, PET) 在肿瘤学、心血管疾病学、神经系统疾病学和新医药学开发等研究领域有着卓越的性能[1]。磁共振成像 (Magnetic resonance imaging, MRI) 在反映解剖形态和生理功能信息方面也具有无可比拟的优越性，特别是在软组织对比方面，其优势更强且无辐射[2]。一体化PET/MRI扫描在诊断、预测、监控许多重大疾病的发生与转移中发挥着重要的作用[3]。头部PET/MRI扫描对于脑血管疾病、老年性痴呆、癫痫、帕金森、神经退行性疾病的诊断以及神经精神药物研究与脑功能研究等[4]具有重要价值。在一体化PET/MRI同步扫描过程中，由一些疾病（癫痫、帕金森等）引起的头部运动或者由于患者肌肉收缩与舒张带来的运动不可避免，由此带来的运动伪影对于临床上判断病灶的大小与发生位置产生了很大的影响。
为了获取被测物体更精确的运动信息，减少运动伪影对PET图像质量的影响，国内外已经做了很多相关研究。对头部加绑带以抑制运动和麻醉被测患者都可以减少运动伪影，但前者并不能完全抑制患者的运动[5]，而后者会干扰患者神经系统功能[6]。现在较为常见的PET头部运动校正包括自适应校正[7-9]和利用外部设备记录被测物体运动信息的方法[10-11]。可借助的外部设备包括视频监控系统等。由于头部运动是刚性运动，某一点的运动情况即可代表整个头部的运动，用视频监控系统监测头部某一点6个维度的运动信息，以此用于头部运动校正[10]。
在以往的各种手段中，所获取的被测物体的运动信息精确度有限，因此基于帧 (frame) 和基于响应线 (Line of response, LOR) 的两种头部运动校正方法校正效果区别不甚明显[12]。PET/MRI出现以后，相比于之前的运动信息获取手段，由于MRI成像速度快、图像质量高、序列选取灵活，因此基于MRI获取被测模体的运动信息的精确度大为提高，PET的头部运动校正有了长足的发展。
本文旨在验证一体化PET/MRI应用中，由于获取的运动信息精确度高，那么基于LOR的方法则具有更好的校正效果。通过对快速成像的MRI图像进行配准、利用导航序列记录被测物体的运动轨迹等都是PET/MRI中比较重要的快速精确获取运动信息的方法[3, 13]，这为PET/MRI的运动校正奠定了基础和方向，也为临床上选取最适合的提高PET图像质量的方法提供了依据。本文在一体化PET/MRI同步扫描过程中，通过对快速成像的高空间分辨率MRI图像进行配准，获取被扫描模体的高精度运动信息。利用获取到的运动信息对PET原始数据分别进行基于两种方法的头部运动补偿。为比较分析两种补偿方法的效果，拟针对模拟和实验两部分展开研究。使用模拟数据进行校正方法的设计和调试，使用实验数据对校正算法进行验证。本文对两种方法的校正结果做了定性和定量的对比研究。
1 方法 1.1 一体化PET/MRI系统参数我们使用上海联影医疗科技有限公司自主研发的一体化PET/MRI (uPMR790) 设备，该设备结合PET与MRI两种模态的优点，采用自主研发的硅基光电转换探测技术，使用2.76 mm×2.76 mm× 15.50 mm的硅酸钇镥闪烁晶体，7×8的晶体组成一个阵列块，5×14的阵列块组成一个模块。PET探测器共有20个模块，轴向112环，单环700个晶体，内环直径647 mm。实现了32cm的PET成像轴向视野和1.4 mm的空间分辨率，具有飞行时间 (Time of flight, TOF) 功能。
该设备还配备了3 T超导磁体和50 mT·m-1的梯度线圈，最大爬升率200 mT·m-1·ms-1。
1.2 系统模拟使用蒙特卡罗 (Monte Carlo) 方法[14]，利用标准的医学影像模拟软件GATE (GEANT4 Application for Tomographic Emission) 模拟真实PET探测器的结构、物理、模体、放射源以及采集条件和采集时间。GATE模拟能够精确地模拟PET扫描中正电子的产生、湮灭以及产生的光子到达探测器等物理过程[15]。
设置PET扫描时NEMA IQ (National Electrical Manufacturers Association Imaging Quality)[16]标准模体以系统中心为原点，沿系统坐标轴6个方向各运动15 mm，模体在每个位置模拟PET扫描30s，分别存储在6个位置的PET模拟数据，模体如图 1所示。
|
图 1 NEMA IQ模体 Figure 1 NEMA IQ phantom. |
类似地，设置XCAT (the 4D extended cardiac-torso)[17]头部模体沿系统坐标轴6个方向各运动8 mm，并存储6个位置的PET模拟数据。为了更加精确地分析和对比两种校正方法的效果，系统仿真模拟时在XCAT头部加入一个10 mm×10 mm× 10 mm的球形病灶，如图 2所示。
|
图 2 XCAT头部矢状面图像 Figure 2 Sagittal plane of XCAT head. |
利用MRI快速梯度回波序列 (Gradient echo quick 3D) 成像速度快和MRI图像空间分辨率高的特点进行实验，参数为：矩阵大小77×96，像素尺寸3.12 mm×3.13 mm，层数56，层厚3 mm，回波时间1.8ms，重复时间4.5ms。
实验的方法步骤如下：
1) 使用水分较丰富的柚子，将三个PET实验用点源（22Na源）嵌入柚子表皮中并绑定；
2) 将上述绑定了点源的柚子置于PET/MRI扫描床板上，并使用激光灯定位；
3) 设置PET/MRI同步扫描环境，利用快速梯度回波序列进行连续不间断采集25次，获取25帧MRI图像；
4) 在此扫描中保持PET一直采数，采集期间推动柚子发生4次运动；
5) 存储MRI的25次扫描完成时采集到的PET原始数据。
本实验过程中PET一直保持采数，PET/MRI总共同时扫描约280.1s。
实验过程中，柚子发生了4次运动，对应5个空间位置，每个位置都有对应的MRI图像与PET图像。如图 3所示，实验在t0时刻开始扫描，以第一幅MRI图像为参考图像，后续24幅MRI图像与该参考位置配准获取运动信息。模体分别在t1、t2、t3、t4时刻发生了运动，实验在t5时刻结束。分别存储每个位置相对于参考位置的运动矢量，并以t1-t44个时间点为分界点，把PET原始数据分段存储，根据运动补偿方法的不同（基于frame的方法和基于LOR的方法），分别存储为弦图格式sinogramn和表模式listmoden [18]，n=1, 2, 3, 4, 5。
|
图 3 PET/MRI同步扫描 Figure 3 Simultaneous scanning of PET/MRI. |
将采集到的各帧PET数据分别重建后，根据MRI图像配准得到各帧的平移、旋转等运动信息，将各帧图像做相应的平移、旋转反变换，并将变换后的各帧图像相加[7-8]，即获得运动校正后的PET图像。
1.4.2 基于LOR的方法PET采集的原始数据可以采用表模式listmode记录。表模式数据以一维数据表的形式记录了全部事件信息，包括每个γ光子的入射位置、发生时间、能量信息等，直到采集结束。根据一对光子落在探测器环上的位置，PET采集的每一个事件，可以对应一条LOR，因此每一条LOR都有相应的4个信息 (Na、Nb、Ra、Rb)，Na表示a光子所在环上的晶体编号，Ra表示a光子所在探测器环的编号。
如图 4所示，以探测器环所在平面的水平方向为x轴，竖直方向为y轴，从机架指向床板的方向为z轴，建立空间坐标系，探测器可近似为该坐标系内的一个圆柱，坐标原点位于圆柱中心，每一条LOR则对应该空间坐标系内圆柱上的两个点M1 和M2。
|
图 4 基于LOR的校正示意图 Figure 4 Procedure of LOR based motion compensation. |
通过MRI图像配准获取到LOR的运动矢量为α，对M1M2变换到N1N2的位置，变换如下：
| $ \overrightarrow {{N_1}{N_2}} = \alpha \overrightarrow {{M_1}{M_2}} $ | (1) |
由于N1点与N2点不一定落在圆柱壁上，需将N1、N2延长，取其与圆柱的交点
重建基于开源重建软件STIR (Software for tomographic image reconstruction)[19]，采用顺序子集期望最大算法 (Ordered subsets expectation-maximization algorithm, OSEM) 进行图像重建[20]。重建使用9次迭代，5个子集，重建图像尺寸为128×128×223，像素尺寸为2.778 mm×2.778 mm× 1.425 mm。
为了获得精确的PET图像，还需要对PET数据进行各项校正，包括归一化校正、随机校正、衰减校正、散射校正等，与运动校正一样这些校正都基于STIR完成。
1.6 数据分析方法 1.6.1 半高宽PET图像运动校正的效果通常通过比较校正前后病灶点的半高宽 (Full width at half maximum, FWHM) 来分析，如图 5所示。
|
图 5 点扩展函数半高宽定义示意图 Figure 5 FWHM of point spread function. |
由点扩展函数的最大值的一半处的相邻像素间的线性插值即可确定FWHM，最大值用峰值点和其最近邻的两个点的抛物线拟合来确定。
对于XCAT头部模拟数据和绑定了点源的柚子实验数据，本文都分别获取了两种校正方法下的PET重建图像上病灶区域随机5个位置的FWHM来分析校正效果，FWHM越小校正效果越好。
1.6.2 对比度恢复系数对NEMA IQ模体的重建图像通常通过计算对比度恢复系数 (Contrast recovery coefficient, CRC) 来分析图像质量。分析中使用以冷球和热球为中心的横断面图像，在每个热球和冷球体画出相应的感兴趣区域 (Region of interest, ROI)，并在该层的模体本底上画出与冷热球上所画ROI相同尺寸的ROI，接近中间层两侧±1 cm与±2 cm处的其它层上也应该画出ROI。应画不同大小的本底ROI共计60个，共5层，记录每个ROI上本底的平均计数。每个热球j的百分对比度QH, j由式 (2) 计算：
| $ {Q_{{\rm{H}}, j}} = \frac{{{C_{{\rm{H}}, j}}/{C_{{\rm{B}}, j}} - 1}}{{{a_{\rm{H}}}/{a_{\rm{B}}} - 1}} \times 100\% $ | (2) |
式中：CH, j为球体j上ROI内平均计数；CB, j为球体j本底ROI上计数的平均值；aH为热球体内放射性活度浓度；aB为本底放射性活度浓度。对于每个冷球体j的百分对比度QC, j由式 (3) 计算：
| $ {Q_{{\rm{C}}, j}} = \left( {1 - \frac{{{C_{{\rm{C}}, j}}}}{{{C_{{\rm{B}}, j}}}}} \right) \times 100\% $ | (3) |
式中：CC, j为球体j上ROI内的平均计数；CB, j为球体j 60个本底ROI计数的平均值[16]。
分别统计了未运动、未经运动校正和两种运动校正方法后NEMA IQ模体图像的冷、热球体的百分对比度随着重建的迭代次数的变化，即绘成CRC曲线。QH, j和QC, j越高，校正效果越好。
2 结果为了对比相同条件下两种校正方法校正效果的差异，两种校正方法所用PET原始符合事件数保持一致。
2.1 模拟数据对XCAT头部模拟数据进行标准的PET校正、重建，其结果如图 6所示。
|
图 6 XCAT头部校正前后的重建图像 (a) 未运动，(b) 未校正，(c) 基于frame的校正，(d) 基于LOR的校正 Figure 6 Reconstruction images of XCAT head before and after motion compensation. (a) With no motion, (b) Without motion correction, (c) With frame based motion correction, (d) With LOR based motion correction |
分析图 6中 (a)-(d) 图像有病灶区域中心附近随机5个位置的FWHM，结果如图 7所示。
|
图 7 校正前后图像病灶区域的FWHM值 Figure 7 FWHM values around lesion of reconstruction images before and after motion correction. |
对NEMA IQ模体的模拟数据进行标准的PET校正、重建，运动校正前后的重建结果如图 8所示。
|
图 8 NEMA IQ模体校正前后的重建图像 (a) 未运动，(b) 未校正，(c) 基于frame的校正，(d) 基于LOR的校正 Figure 8 Reconstruction images of NEMA IQ phantom before and after motion compensation. (a) With no motion, (b) Without motion correction, (c) With frame based motion correction, (d) With LOR based motion correction |
经计算，两种方法校正前后NEMA IQ图像的CRC曲线如图 9所示。
|
图 9 校正前后NEMA IQ图像冷、热球的对比度恢复系数曲线 (a) 未运动，(b) 未经运动校正，(c) 基于frame的运动校正，(d) 基于LOR的运动校正 Figure 9 CRC curves of NEMA IQ hot and cold spheres before and after motion compensation. (a) With no motion, (b) Without motion correction, (c) With frame based motion correction, (d) With LOR based motion correction |
图 10(a)-(e) 分别为5个位置上MRI横断面图像，图 10(f)-(j) 分别为对应的5个位置上的PET图像。
|
图 10 5个位置上的MRI横断面成像 (a-e) 和PET图像 (f-j) Figure 10 MRI transverse images (a-e) and PET image (f-j) at 5 positions. |
未经运动校正的所有PET数据重建图像如图 11(a)所示，对PET原始数据进行基于frame和基于LOR的运动校正后的重建图像分别如图 11(b)和(c)所示。
|
图 11 钠源运动校正前后的重建图像 (a) 未校正，(b) 基于frame的校正，(c) 基于LOR的校正 Figure 11 Reconstruction images of sodium sources before and after motion correction. (a) Without motion correction, (b) With frame based motion correction, (c) With LOR based motion correction |
分析图 11中 (a)-(c) 图像点源中心区域随机5个位置的FWHM，结果如图 12所示。
|
图 12 校正前后点源中心区域的FWHM值 Figure 12 FWHM value of center of point sources before and after motion correction. |
未经运动校正的PET/MRI融合图像如图 13(a)所示，经过基于frame和基于LOR的方法校正的融合图像分别如图 13(b)和(c)所示。
|
图 13 运动校正前后的PET/MRI融合图像 (a) 未校正，(b) 基于frame的校正，(c) 基于LOR的校正 Figure 13 PET/MRI combined image before and after motion correction. (a) Without motion correction, (b) With frame based motion correction, (c) With LOR based motion correction |
对于XCAT模拟数据，由图 7可以看出，5个位置上，两种运动校正方法下病灶附近区域的FWHM值都比未校正的图像FWHM值要小，可见两种校正方法都有明显的补偿运动伪影的作用；同时5个位置上，基于LOR的校正图像中FWHM值都要比基于frame的校正图像中的小，且更接近未运动的模体重建图像中的FWHM，因而基于LOR的方法要比基于frame的方法校正效果好。
对于NEMA IQ模拟数据，由图 9(a)-(d)对比可以看出，相对于未校正的图像，基于frame和基于LOR的方法中随着迭代次数的增加，冷球和热球的CRC值都明显上升，因此这两种方法都有明显的校正效果；而基于LOR的方法相对于基于frame的方法，冷球和热球的CRC值都上升得更多，且更接近于未运动的模体重建图像的CRC值，因此前者具有更好的运动校正效果。
对于点源实验数据，由图 12可以看出，结果与图 7类似，两种校正方法都有明显的补偿运动伪影的效果，而基于LOR的校正效果更为优异。
由上述模拟和实验结果可以看出，两种方法都具有明显的校正效果，但也各有优缺点。基于frame的运动校正是在重建后的PET图像上的操作，相对于基于LOR的方法更为简单方便；但PET图像重建过程会引入一定的噪声，因而校正效果不如基于LOR的方法。而基于LOR的方法是对重建前的PET原始数据操作，不受重建过程带来的其他噪声的影响，因此在获取同样精度的运动信息的条件下，相比于基于frame的方法校正效果更为明显。另外，基于frame的校正方法对PET图像质量有一定的要求，而高质量的PET图像需要足够多的符合事件。在一定剂量的放射源下，统计量与扫描时间正相关。基于LOR的校正方法的校正过程不需要PET重建，只需要MRI成像，而MRI有一些快速成像的序列，可以获取更多的运动信息，因此，基于LOR的校正方法在MRI的序列选取上也相对更为灵活。
PET/MRI扫描中，根据不同的临床需求，使用的MRI序列不同，MRI成像时间有所不同，图像质量亦有所不同。但是，成像时间越短，MRI的图像质量却相应地下降，给图像配准带来了一定的影响，获取运动信息的精确度下降。因此，在MRI成像时间对获取到运动信息准确性的影响和配准算法对MRI图像质量的要求中间寻求一个均衡点，是PET/MRI头部运动校正未来的重要研究内容之一。
4 结语PET/MRI同步扫描实验以及模拟数据结果都很好地验证了两种方法的校正效果。XCAT头部模体和NEMA IQ模体的模拟结果表明，基于LOR的运动校正结果与没有运动的模体图像更为接近，绑定了点源的柚子的PET/MRI同步扫描实验也进一步证明了基于LOR的方法校正效果更为明显。快速精确成像的MRI图像能够获取到足够准确的运动信息，这都为临床上减少PET/MRI头部图像伪影，提高图像质量提供了研究依据。
| [1] | Strauss L G, Conti P S. The applications of PET in clinical oncology[J]. Journal of Nuclear Medicine, 1991, 32(4): 649–650. |
| [2] | Constantine G, Shan K, Flamm S D, et al. Role of MRI in clinical cardiology[J]. Lancet, 2004, 363(9427): 2162–2171. DOI: 10.1016/S0140-6736(04)16509-4 |
| [3] | Catana C, Benner T, van der Kouwe A, et al. MR-assisted PET motion correction for eurological studies in an integrated MR-PET scanner[J]. Journal of Nuclear Medicine, 2011, 52(1): 154–161. DOI: 10.2967/jnumed.110.079343 |
| [4] | Schlemmer H P, Pichler B J, Schmand M, et al. Simultaneous MRPET imaging of the human brain: feasibility study[J]. International Journal of Medical Radiology, 2008, 248(3): 1028–1035. |
| [5] | Bloomfield P M, Spinks T J, Reed J, et al. The design and implementation of a motion correction scheme for neurological PET[J]. Physics in Medicine & Biology, 2003, 48(8): 959–978. |
| [6] | Lindauer U, Villringer A, Dirnagl U. Characterization of CBF response to somatosensory stimulation: model and influence of anesthetics[J]. American Journal of Physiology, 1993, 264(2): 1223–1228. |
| [7] | Ashburner J A, Jesper L R, Friston K J. Image registration using a symmetric prior-in three dimensions[J]. Human Brain Mapping, 2015, 9(4): 212–225. |
| [8] | Jenkinson M, Smith S. A global optimisation method for robust affine registration of brain images[J]. Medical Image Analysis, 2001, 5(2): 143–156. DOI: 10.1016/S1361-8415(01)00036-6 |
| [9] | Maes F, Collignon A, Vandermeulen D, et al. Multimodality image registration by maximization of mutual information[J]. IEEE Transactions on Medical Imaging, 1997, 16(2): 187–198. DOI: 10.1109/42.563664 |
| [10] | Bloomfield P M, Spinks T J, Reed J, et al. The design and implementation of a motion correction scheme for neurological PET[J]. Physics in Medicine & Biology, 2003, 48(8): 959–978. |
| [11] | Lopresti B J, Russo A, Jones W F, et al. Implementation and performance of an optical motion tracking system for high resolution brain PET imaging[J]. IEEE Transactions on Nuclear Science, 1999, 46(6): 2059–2067. DOI: 10.1109/23.819283 |
| [12] | Wang J, Hu L, Feng T, et al. Improved spatial and temporal resolution of gated myocardial perfusion PET using post reconstruction dual respiratory and cardiac motion compensation[C]. Nuclear Science Symposium and Medical Imaging Conference, IEEE, 2014: 1-4. |
| [13] | Thurfjell L, Lötjönen J, Lundqvist R, et al. Combination of biomarkers: PET[J]. Neurodegenerative Diseases, 2012, 10(1-4): 246–249. DOI: 10.1159/000335381 |
| [14] | Wasserstein R L. Monte Carlo: concepts, algorithms, and applications[M]. New York, USA: Springer-Velag, 1996. |
| [15] | Jan S, Santin G, Strul D, et al. GATE: a simulation toolkit for PET and SPECT[J]. Physics in Medicine & Biology, 2004, 49(19): 4543–4561. |
| [16] | Twitchen D J, Newton M E, Baker J M, et al. Performance measurements of positron emission tomographs[M]. Virginia, US: National Electrical Manufacturers Association, 2001. |
| [17] | Segars W P, Sturgeon G, Mendonca S, et al. 4D XCAT phantom for multimodality imaging research[J]. Medical Physics, 2010, 37(9): 4902–4915. DOI: 10.1118/1.3480985 |
| [18] | Nichols T E, Qi J, Asma E, et al. Spatiotemporal reconstruction of list-mode PET data[J]. IEEE Transactions on Medical Imaging, 2002, 21(4): 396–404. DOI: 10.1109/TMI.2002.1000263 |
| [19] | Thielemans K, Tsoumpas C, Mustafovic S, et al. STIR: software for tomographic image reconstruction release 2[C]. Nuclear Science Symposium Conference Record, IEEE, 2012: 2174-2176. |
| [20] | Boellaard R, Van L A, Lammertsma A A. Experimental and clinical evaluation of iterative reconstruction OSEM in dynamic PET: quantitative characteristics and effects on kinetic modeling[J]. Journal of Nuclear Medicine, 2001, 42(5): 808–817. |