PET断层重建中动态射线追踪算法的实现

张 斌, 王李栓， 赵书俊

(1.郑州大学 物理工程学院 河南 郑州 450001； 2.郑州市节能监察中心 河南 郑州 450006)

0 引言

正电子发射断层扫描仪(positron emission tomography, PET)利用放射性核素标记化合物作为分子探针，非侵入式地获取生物活体的功能和解剖影像，其功能成像是当前所有成像模式中灵敏度最高的，是获得生物活体的生理生化信息的强有力工具，近年来在诸多科学研究领域得到广泛应用.在临床医学中，已经证明PET在癌症诊断和分期、神经系统疾病评估、心脏病诊疗、放疗计划制定以及化疗监测等方面起到重要作用.在生物医学工程中，已经广泛地利用小动物PET进行基因传送与表达、细胞贩运、信号传输通路中的蛋白质相互作用等方面的研究[1].

在PET扫描过程中，不改变活体对象的生理状态，向其注射放射性示踪剂参与代谢.示踪剂标记物是构成有机体的基本元素的正电子同位素，如11C、13N、15O、18F或其他核素.示踪剂标记物衰变发射正电子，发生正电子湮没反应，产生逆向发射的511 keV的γ光子对.使用符合探测技术探测成对出现的γ光子对，确定符合响应线(line of response, LOR)，通过数据采集得到大量LOR，加以校正和滤波后，进行图像断层重建，产生示踪剂的时间序列浓度分布影像，最后得到活体代谢功能影像.图像断层重建算法是PET设备的关键组成部分，是决定PET设备性能与成像质量的重要环节.目前主要包括迭代重建算法与解析重建算法，两者相比，迭代重建算法能够提供更好的图像重建质量，已成为PET的标准重建方法[2].

迭代重建过程需要借助SRM，频繁在图像域和投影域之间进行变换运算，其执行流程如图1所示.当变换运算沿着一条射线路径进行时，称为“线驱动的前向投影和背投影运算”.前向投影需要计算每一个像素对投影值的贡献.背投影利用之前计算的校正因子，对估算的像素值进行修正.这两个计算过程都要反复用到像素的权重因子，需要占用大量的系统资源，是重建算法中的关键步骤和最耗时的部分.

图1 图像域和投影域之间的变换示意图Fig.1 Diagram of the transformations between the image space and the projection space

SRM由诸多物理因素决定，包括扫描仪几何效应、随机符合、散射、被扫描对象内衰减、正电子行程、非共线性、晶体穿透、晶体内散射等.任何情况下，几何效应都是主要因素.在PET中，巨量LOR线和像素空间中同样数量庞大的像素组合所产生的SRM非常巨大，无论计算、存储和调用都是迭代重建算法中的瓶颈.SRM可以离线计算并存储在磁盘上，在重建过程中调入内存使用.但由于过于庞大，通常不能一次全部装入内存，需要反复从磁盘上读取，成为制约重建速度的一个主要因素.SRM也可以在重建过程中动态计算，即所谓on the fly模式，但计算强度非常大.PET中采集的数据主要有正弦图数据和list-mode数据两种组织模式.list-mode数据是以列表的模式逐一记录每个湮没事件的LOR线空间位置、发生时间、能量以及其他相关信息.正弦图数据是通过对采集到的list-mode数据按照LOR线分类重组而获取.

目前，基于正弦图的迭代断层重建算法在PET中得到广泛应用.基于list-mode数据的重建算法具有许多独特优势，近年来，受到人们的高度关注，处于快速发展和完善过程中[3-5].对于正弦图重建，SRM可以预先计算并存储在硬盘上，也可以采用动态计算模式.而基于list-mode的断层重建，迭代计算是一个事件接一个事件进行，实时计算系统矩阵，可以有效降低数据存取时间和存储负担，提高断层重建的速度.实时计算时，为了提高计算效率，主要考虑几何效应，并选用适当的高斯核，图像质量几乎不受影响[6-7].

本文面向list-mode断层重建，实现了一种敏捷的实时计算系统矩阵的方法，并融合在前向投影过程中.

1 正交距离射线追踪算法原理与实现

该算法以优化Siddon算法为基础，考虑PET探测特性，生成近似高斯形的LOR线，有效建模了探测器响应函数.

1.1 Siddon算法的基本原理

Siddon算法在图像断层重建中得到广泛应用，该算法是找到被LOR线穿越的全部像素并计算LOR线在每个像素中的截断长度，用截断长度代表像素相对于该LOR线的权重因子.它的时间复杂度为O(n)，而直接的射线追踪方法时间复杂度为O(n3).它首先对视场空间像素进行有序化标记，再对LOR线进行参数化表达，然后找到LOR线进入视场空间的入射点和出射点，从而确定射线追踪顺序.这些操作构成Siddon算法的实现基础[8].

2D情况下，用2组分别垂直于x、y轴的平行等间距直线(3D时，用3组分别垂直于x、y、z轴的平行等间距平面)，把视场空间划分成Nx×Ny的像素阵列，实现对视场空间的有序标记，如图2所示.对任一像素v(x,y)，可用(i,j)标记，其中1≤i≤Nx，1≤j≤Ny.

假如P1(x1,y1)、P2(x2,y2)代表2个探测器的中心点，连接2中个心点的射线代表一条LOR线，可参数化表示为：x(α)=x1+α(x2-x1);y(α)=y1+α(y2-y1),其中α∈[0, 1]，表示LOR线上任意一点到起始点的距离与LOR线总长度的比值.参数方程实际上分别代表LOR线在x、y轴上的投影线.

沿着P1到P2，或者P2到P1的方向，分别计算2个投影线与2组平行线的交点，包括入射点和出射点，剔除重复项，并按照升序排列，得到参数化表示的与LOR线相交的一系列像素{α(0),α(1),…,α(n)}.

被LOR线穿越的第m个像素在像素空间中的索引，可通过i(m)=⎣x1+αmid(x2-x1)」和j(m)=⎣y1+αmid(y2-y1)」计算求出.

需要注意的是，在实践中可以根据需要，灵活建立参考坐标系，不会影响最终计算结果，例如，可以将坐标原点设定在视场中心位置.像素空间中像素的排序以坐标系为参照，程序自定义.

1.2 优化的Siddon算法

文[9]通过对Siddon算法的优化，将射线追踪的计算时间减少2/3.在优化算法中，首先求第一个被射线击中的像素的索引坐标，然后通过加、减运算，依次求出射线路径上其余的像素索引坐标，把大量乘法运算简化为加减运算.引入一个新的参量λ替代α，参数化公式转化为：x(λ)=x1+(λ/L)·(x2-x1),y(λ)=y1+(λ/L)·(y2-y1),其中，λ表示当前点到LOR线起点的距离.LOR线在像素中的截断长度l(m)=λ(m)-λ(m-1).

优化算法的主要计算步骤如下：

1)初始化.初始化计算LOR线穿越像素空间时λ的取值范围，即最小值λmin和最大值λmax；初始化计算LOR射线进入图像空间后被x、y方向轴对齐等间距平行线的截断初值λx、λy；初始化计算被射线击中的第一个像素的索引坐标(ix,jy)，其中1≤ix≤Nx，1≤jy≤Ny.

2)找到λx和λy中的最小者，设λξ=min(λx，λy).如果λξ>λmax，则射线追踪结束.否则计算LOR线在当前像素中的横断长度l=λξ-λ.

3)按照单位像素大小，步进增加λξ值，λξ=λξ+L/(ξ2-ξ1),ξ或者是x或者是y.

4)根据上一像素索引计算下一像素索引.然后返回第2)步，循环执行，直到满足结束射线追踪条件.

1.3 正交距离射线追踪算法的实现

用Siddon优化算法生成的LOR射线较细，不能够很好匹配有一定宽度的探测器.因此，具有一定宽度的LOR线(3D时，用响应管(tube of response, TOR))建模探测器之间的响应会更加精确.正交距离射线追踪法以Siddon优化算法为基础，生成近似高斯形的LOR线，可有效建模探测器响应函数，计算公式为aij=1-dij/(cδ),其中dij表示像素几何中心到LOR线的距离，cδ表示归一化因子.δ与PET系统点扩展函数的半高宽相关，c是经验调整系数,aij表示权重因子.

在实时计算系统矩阵时，首先用Siddon算法，找到被LOR线穿越的像素.然后，设置aij的最小阈值，并结合cδ值，计算LOR穿过影像空间时，被LOR线直接穿越的像素及其邻接像素的dij，求出全部关联像素的权重因子.正确设置aij的最小阈值非常关键，需要考虑在计算时间和重建图像品质因数之间进行折中.

可用矢量公式计算各像素的dij值[10].假设LOR线两端晶体条中心点坐标分别为P1(x1,y1)和P2(x2,y2)，求矢量r在v上的投影，也就是像素几何中心P0到LOR线的距离，如图3所示.

生成的LOR线效果如图2所示.图像空间左下角坐标为(bx,by)；dx,dy分别表示x轴平行线和y轴平行线间的距离；Nx,Ny分别表示x轴平行线和y轴平行线的个数；dij表示像素几何中心到LOR线的距离；P1P2确定一条LOR线，阴影部分表示用正交距离法求出的具有一定宽度的LOR线，每个像素的灰度表示其权重.

图2 正交距离射线追踪法原理示意图Fig.2 Schematic of OD-RT method

图3 矢量法计算像素中心到直线距离示意图Fig.3 Schematic of calculating the distance from geometric centre of pixel to LOR using vector method

2 实验装置与材料

Eplus-166是中国科学院高能物理研究所研制的我国第一台拥有自主知识产权的小动物PET.该扫描仪由16个模块组成，呈正十六边型，每个模块由2个block沿着轴向排列组成，每个block包含16×16个晶体条.晶体条的大小为1.9 mm×1.9 mm×10 mm，为了提高闪烁光收集效率，晶体条之间加入了反光材料，使得探测器单元沿轴向和切向断面长度都为2.0 mm.系统总共形成32个探测器环，轴向长度为64 mm，其结构如图4所示[11].

实验采用自制的Derenzo仿体模型，模型内各空心柱直径分别为 1.4,1.6,1.9,2.2,2.5,3.0 mm，空心柱圆心间的距离为其直径的2倍.空心柱内均匀注入18F-FDG溶液，在Eplus-166中扫描，对采集的数据进行归一化、衰减以及散射校后，进行list-mode迭代重建.

3 List-mode重建算法

图像重建采用S-LMEM (ordinary subsetized list-mode EM algorithm)算法，迭代公式为

4 结果分析与结论

本文以2D重建进行算法验证，采用S-MLEM算法，使用正交距离射线追踪方法以on-the-fly模式实时计算系统矩阵.把list-mode数据划分为32个子集，每个子集约50 000个事例.Eplus-166中心视场的点扩展函数半高宽约1.67 mm[11]，分别设cδ=3, 2, 1.5，当cδ=2时，达到最佳效果，能够清晰地区分出1.9 mm热区，结果如图5所示.

图4 Eplus-166结构示意图Fig.4 Schematic for Eplus-166 structure

图5 重建结果图Fig.5 Comparison of image reconstruction result using S-MLEM

从重建结果可以看出，当设置cδ为合适值时，内含的线性核能够更好地近似高斯模糊函数，有效建模探测器响应函数，提高系统矩阵的计算精度，对于分辨率恢复起到决定性作用.实践中发现，aij的下限越小，概率计算精度越高，但计算强度迅速增加.因此，aij并非越小越好，要根据重建目标，合理选择aij的下限值，达到图像重建质量和时间之间的平衡.

无论采用何种算法，PET断层重建的计算强度都非常大.从list-mode算法的原理与实现来看，在一次迭代中某个像素值的更新不受其他像素更新的影响，整个计算过程特别适合采用多线程并行计算，进行实时重建.CUDA技术非常适宜于实现多线程list-mode断层重建[12-14]，在本研究中进行了初步尝试，重建速度明显提高.

本文的研究表明，利用正交距离射线追踪方法以on-the-fly模式实时计算系统响应矩阵，进行PET动态重建是可行的.这对于进一步研究面向list-mode数据的低统计计数率(低放射性活度)情况下的断层重建以及实现包含TOF信息的断层重建都有非常重要的借鉴意义.

参考文献：

[1] Hao Peng , Craig S L. Recent developments in PET instrumentation[J]. Current Pharmaceutical Biotechnology, 2010, 11(6):555-571.

[2] Reader A J, Zaidi H. The promise of new PET image reconstruction[J]. Physica Medica, 2008, 24:49-56.

[3] Barrett H H, White T, Parra L C. List-mode likelihood[J]. J Opt Soc Am A, 1997, 14(11):2914-2923.

[4] Parra L , Barrett H H. List-mode likelihood: EM algorithm and image quality estimation demonstrated on 2-D PET[J]. IEEE Trans Med Imaging, 1998, 17(2):228-235.

[5] Rahmin A, Lenox M, Reader A J, et al. Statistical list-mode image reconstruction for the high resolution research tomograph[J]. Phys Med Biol, 2004, 49:4239-4258.

[6] Aguiar P, Rafecas M, Ortuo J E, et al. Geometrical and Monte Carlo projectors in 3D PET reconstruction[J]. Med Phys, 2010, 37(11):5691-5702.

[7] Colas S. A fast tube of response ray-tracer[J]. Med Phys, 2006, 33(12):4744-4748.

[8] Jacobs F, Sundermann E, Sutter B D, et al. A fast algorithm to calculate the exact radiological path through a pixel or voxel space[J]. Journal of Computing and Information Technology, 1998, 6(1):89-94.

[9] Han Gouping, Liang Zhengrong, You Jiangsheng. A fast ray tracing technique for TCT and ECT studies[C]// IEEE Nuclear Science Symposium Conference Record. Washington, 1999:1515-1518.

[10] Weisstein E W. Point-line distance 2-dimensional[EB/OL].[2011-12-11]. http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html.

[11] 贠明凯．正电子发射断层成像中的三维迭代重建研究[D]．北京：中国科学院高能物理研究所，2009．

[12] 刘琳，何剑锋，王红玲．GPU加速数据挖掘算法的研究[J]．郑州大学学报：理学版，2010，42(2)：31-34．

[13] Pratx G, Chinn G, Habte F, et al. Fully 3-D list-mode OSEM accelerated by graphics processing units[C]// IEEE Nuclear Science Symposium Conference Record. San Diego CA, 2006:2196-2202.

[14] 王君，任永功.基于FP-tree 最大频繁模式超集挖掘算法[J].郑州大学学报：理学版，2011，43(1)：33-36.