PET时间实时校正及重建方法和系统与流程

本发明涉及医学影像设备技术领域，尤其涉及一种pet时间实时校正及重建方法和系统。

背景技术：

正电子发射断层设备的原理是测量正电子湮没后的γ光子配对，并利用的γ光子对的几何方向和时间特性进行成像。其中γ光子对的时间特性是指，两个方向相反的光子由于飞行时间不同，进入检测器闪烁体晶体的时间产生差异。这个时间差产生的原因是由于晶体材质，晶体形态和前、后端电子学性能的影响，在每一个晶体上产生一定的偏差，从而使符合数据的时间值和真实值之间有所差异。

为了消除这种差异，文献[mewernerandjskarp,tofpetoffsetcalibrationfromclinicaldata.phys.med.biol.58(2013)4031–4046]中使用了一种基于nontof的结果上来统计时间偏差的方法，该方法设定了一个固定的最小时间差，然后使用histgram的方法统计每个晶体上符合事件分布和在nontof图像的时间分布对比处理。该方案存在的以下问题：1.最小精度由最小时间差决定，一旦确定就无法继续提高精度，但由于最小时间差是一个固定值，在晶体时间差异很大时，无法统计，造成严重误差或失效；2.在计算时间校正时需要重新计算scatter，浪费了计算效率。

技术实现要素：

为解决上述问题，本发明的目的在于提供一种pet时间实时校正及重建方法和系统，通过使用经过校正处理后的listmode数据，在不需要重新进行scatter计算或者其他插值计算的情况下，快速计算出晶体延时，同时逐步迭代提高计算精度，保证了结果的稳定性和准确度。

为了实现上述目的，本发明的技术方案如下：

一种pet时间实时校正及重建方法，包括以下步骤：

s1.收集符合事件，形成可读取和使用的listmode数据，并确定随机分数randfra和散射分数scatfra；

s2.对listmode数据进行无飞行时间的重建，获取辐射源的空间分布voxel[x,y,z]，其中xyz为空间三维坐标；

s3.累积计数listmode数据的发生次数，获取晶体的时间直方图分布h1；

累积计数listmode数据发生的校正信息，获取晶体的时间直方图分布h2；

根据辐射源的空间分布的体数据计算得到晶体的时间直方图分布l1；

s4.根据晶体的时间直方图分布，计算晶体时间基准；

s5.对晶体时间基准进行迭代更新，并设置迭代终止条件，若不满足迭代终止条件，则根据更新后的晶体时间基准重新获取晶体的时间直方图分布h1和h2，并进行晶体时间基准的下一次更新；若满足迭代终止条件，则进行图像重建并输出图像。

进一步的，所述步骤s3中的晶体的时间直方图分布h1的获取包括以下步骤：

定义晶体的时间直方图分布h1为一个二维的存储空间h1[c,n]，该空间的大小为c*n，c代表的晶体的序号的最大值，n代表了出现在该晶体上的事件所记录的时间差的分隔数，n个分隔数中每一个分隔所对应的时间宽度为dt；

当符合事件j发生时，与它相关的晶体为c1和c2，时间差为tj，则偏移量tjd＝[tj/dt]，则分别在h1[c1，m-tjd]和h1[c2，m+tjd]上分别记录，记录的方式采用自加，其中m是指的n的1/2，即分隔数的一半。

进一步的，所述步骤s3中的晶体的时间直方图分布h2的获取包括以下步骤：

定义晶体的时间直方图分布h2为一个二维的存储空间h2[c,n]，该空间的大小为c*n，c代表的晶体的序号的最大值，n代表了出现在该晶体上的事件所记录的时间差的分隔数，n个分隔数中每一个分隔所对应的时间宽度为变量dt；

当符合事件j发生时，与它相关的晶体为c1和c2，时间差为tj，则偏移量tjd＝[tj/dt]，则分别在h2[c1,m-tjd]和h2[c2,m+tjd]上分别记录，记录的方式是h2[c1,m-tjd]＝h2[c1,m-tjd]+b^addj，h2[c2,m+tjd]＝h2[c2,m+tjd]+b^addj，b^addj是散射符合事件和随机符合事件的计算值，其中m是指的n的1/2，即分隔数的一半。

进一步的，所述步骤s3中的晶体的时间直方图分布l1的获取包括以下步骤：

定义晶体的时间直方图分布l1为一个二维的存储空间l1[c,n]，该空间的大小为c*n，c代表的晶体的序号的最大值，n代表了出现在该晶体上的事件所记录的时间差的分隔数，n个分隔数中每一个分隔所对应的时间宽度为变量dt；

统计所有从晶体出发的lor线共有lorc条，其中任意一条lor线上，记录lor上任意一个位置的图像像素与lor的中心位置的距离为d，得到时间直方图分布l1的时间差td＝d/vc，偏移量tdoffset＝[d/vc/dt]，其中vc表示光速常量，tdoffset代表了以dt作为单位下的距离；

定义射线追踪法获取到晶体序号pix，偏移量tdoffset时l1的强度为

公式中，lorc是与晶体相关的所有lor，xi(t)，yi(t)，zi(t)是晶体pix的第i条lor上时间差为t的三维坐标，可以通过晶体的几何位置快速推导，σ是系统的时间分辨率，t代表时间变量。

进一步的，所述步骤s4中的晶体时间基准的计算包括以下步骤：

定义l2[c,n]是一个二维的存储空间，则

通过卷积，计算第k个晶体的晶体时间基准tdk，

公式中，h为自变量，优化求解最大值算法argmax使用线搜索(linesearch)实现。

进一步的，所述步骤s5中对晶体时间基准的迭代更新通过下式实现：

td’k＝∈*tdk0+(1-∈)*tdk，

公式中，td’k代表更新后的晶体时间基准，tdk0代表上一次晶体时间基准，∈是经验参数，代表了更新的速度；

根据更新后的td’k，选出晶体c1和c2的晶体时间基准td’c1和td’c2，计算h1和h2中的时间差tj，tj＝tj0+td’c1-td’c2，从而可以获得更新后的h1和h2，tj0为h1和h2中的初始时间差。

进一步的，所述时间宽度dt为变量，且随迭代次数的增加逐次减小，时间宽度dt的最小值为信号采集时间刻度；每迭代一次，即更新一次dt值，并将更新值带入h1、h2和l1中的偏移量的计算式中，进行偏移量的更新。

进一步的，所述时间宽度dt的变化遵循dt[ir]＝dt0*p^ir，其中ir为迭代序号，p为[01]之间的分数，dt0为初始的间隔宽度。

一种pet时间实时校正及重建系统，包括

pet数据采集模块，用于收集符合事件，形成可以读取并使用的listmode数据；

辐射源空间分布获取模块，对listmode数据进行无飞行时间的重建，获得辐射源的空间分布；

符合事件累积模块，累积计数listmode数据的发生次数，得到晶体的时间直方图分布；

校正系数计算模块，累积计数listmode数据的发生的校正信息，得到晶体的时间直方图分布；

lor线追踪累积模块，根据辐射源的空间分布的体数据，计算晶体的时间直方图分布；

晶体时间基准计算模块，根据符合事件累积模块、校正系数计算模块及lor线追踪累积模块中得到的晶体的时间直方图分布，计算晶体时间基准；

迭代更新模块，用于对晶体时间基准进行迭代更新，若更新结果满足迭代终止条件则终止迭代，若更新结果不满足迭代终止条件，则根据更新结果重新计算晶体的时间直方图分布；

重建模块，用于图像重建并输出。

本发明的优点在于：分辨率计算过程联合了图像和计数的多重信息，使分辨率计算结果的可靠性和准确性得到加强，同时逐步迭代提高计算精度，保证了结果的稳定性和准确度，在pet环内的检测器之间的差异大时，本方法有着更好的适应性，能提升重建图像的质量。

附图说明

图1为listmode数据示意图；

图2为实施例中的校正及重建计算流程图；

具体实施方式

以下结合实施例对本发明作进一步详细描述。

对可能出现的参数进行背景介绍响应线：将探测器探测到的γ光子的两个晶体条之间的连线称为响应线(lineofresponse，lor)。

符合时间窗：是为两个γ光子到达探测器的时间差所设的时长。

符合事件：任意两个γ光子在符合时间窗内被探测器探测到的事件。

晶体：是指pet探测器环上检测器前端接受γ光子从而产生前端触发信号的最小几何单元。

真实符合：指湮没辐射产生的两个γ光子，以相反的方向运动并在符合时间窗内被两个晶体分别接收到的符合。

散射符合：指湮没辐射产生的两个γ光子，如果一个到达之前与组织发生散射，但仍在符合时间窗内被探测到的散射符合。

随机符合：是一种假符合，两个γ光子毫无时间和空间的相互关系，但在符合时间窗内被误探测到的符合事件。

listmode数据：是pet机器扫描过程中，记录符合事件的数据，它的特征是：其中的每一个事件至少包括以下信息：位置信息(探测到符合事件的两个γ光子的检测器序号，或者是这两个检测器组成的lor序号)，时间信息(两个光子被探测到的时间差)，校正系数的数据(包括散射和随机事件)，listmode数据示意图如附图1所示。

校正信息：具体描述可见参考文献[wangw,huz,gualtieriee,etal.systematicanddistributedtime-of-flightlistmodepetreconstruction[c]//nuclearsciencesymposiumconferencerecord.ieee,2006]的公式(2)中的b^addj,其中的随机事件的数量和分布分别由符合延迟电路确定，具体操作是从延迟电路采集的事件首先通过乘法因子进行校正，重新绑定到正弦图中，然后使用低通滤波器平滑，参见文献[rokittao,caseym,wienhardkandpictrzyku2000randomcorrectionforpositronemissiontomographyusingsinglescountratesieeenucl.sci.symp.conf.record.337–40]。散射事件的数量和分布由单次散射方法估计计算，参见文献[c.c.watson,“extensionofsinglescattersimulationtoscattercorrectionoftime-of-flightpet”,ieeenss/mic,2005.]

随机分数：定义为随机符合事件数与总符合事件数的比例。

散射分数：定义为随机散射事件数与(散射事件数+真实符合事件数)的比例。信号采集时间刻度：是指符合事件的时间信息的最小记录单位。

pet诊断成像装置，包括：成像区域外部的多个辐射探测器环状结构，用来采集到发出的辐射事件的信息。探测器具有不同的定时延迟。校准体模可以发射时间同时但方向相反的辐射的事件对(符合事件)的辐射源，辐射事件对与辐射探测器的相互作用路径是lor，当检测器环确定后，每一条lor的几何位置是已知的。一旦每个检测器的定时延迟被校准，可以认为符合事件的标记时间是在真实符合发生的位置为中心的，并以特定宽度形成高斯分布。该高斯分布的半高全款宽度就是系统的时间分辨率。

本实施例的pet时间实时校正及重建系统由pet数据采集模块、辐射源空间分布获取模块、符合事件累积模块、校正系数计算模块、lor线追踪累积模块、晶体时间基准计算模块、迭代更新模块及重建模块组成。

pet数据采集模块由外罩、机架及探测器组成，探测器包括晶体、光电放大和电子处理电路。探测器检测从被检体p内部放出的成对湮没γ射线，生成与检测出的成对湮没γ射线的光量相应的脉冲状电信号。后端信号处理器根据来自光电转换器件的电信号生成单事件数据(single)。计数部从重复供给的单事件数据中重复确定容纳在与预先设定的时间范围内的2个单事件有关的事件数据。时间范围被设定为例如6ns～18ns左右。该成对的单事件被推测为由来于从同一成对湮没点产生的成对湮没γ射线。成对的单事件概括地被称为符合事件。符合事件被收集到后，事件相关的时间空间信息进行组合封装，并在磁盘或者其他数据存储介质上保存，形成了重建程序可以读取并使用的listmode数据。同时通过计算总符合、随机符合和估算散射符合的总数，确定出随机分数记为randfra，散射分数记为scatfra。计算随机分数和散射分数的概念和统计公知的，在文献[wangw,huz,gualtieriee,etal.systematicanddistributedtime-of-flightlistmodepetreconstruction[c]//nuclearsciencesymposiumconferencerecord.ieee,2006]中有说明出处。

辐射源空间分布获取模块运行于计算机设备，其中计算机设备包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序。该模块的作用是提取辐射源的实际形态特征，这样可以获得到辐射源的空间分布。它是利用pet诊断成像装置获取到pet符合数据对信号，对信号进行无飞行时间的重建，重建方法例如，有序子集最大似然期望值法osem，最大似然期望值法(mlem)，滤波反投影fbp或者其他方法，从而得到了基本除去噪声和衰减干扰的质量较佳的光子空间分布图。

符合事件累积模块和校正系数计算模块将所有符合事件进行统计，得到每个晶体的时间直方图分布。

lor线追踪累积模块，它的作用是从图像的角度计算晶体的时间直方图分布。

晶体时间基准计算模块，根据符合事件累积模块、校正系数计算模块及lor线追踪累积模块中得到的晶体的时间直方图分布，计算晶体时间基准(即延时)。

迭代更新模块，用于对晶体时间基准进行迭代更新，若更新结果满足迭代终止条件则终止迭代，若更新结果不满足迭代终止条件，则根据更新结果重新计算符合事件累积模块、校正系数计算模块及lor线追踪累积模块中的晶体的时间直方图分布。

重建模块，用于图像重建并输出。

如图2所示，本实施例的pet时间实时校正及重建方法包括以下步骤，

s1.pet数据采集模块收集符合事件，形成了可以读取并使用的listmode数据，并确定随机分数randfra和散射分数scatfra；

s2.辐射源空间分布获取模块利用listmode数据，获取辐射源的空间分布voxel[x,y,z]，其中xyz为空间三维坐标；

s3.符合事件累积模块累积计数listmode数据的发生次数，获取晶体的时间直方图分布h1，具体包括：

定义当第j个符合事件被检测器捕获后，记为事件j，定义晶体的时间直方图分布h1为一个二维的存储空间h1[c,n]，该空间的大小为c*n，c代表的晶体的序号的最大值，也就是检测器上的总的晶体数量，n代表了出现在该晶体上的事件所记录的时间差的分隔数，n个分隔数中每一个分隔所对应的时间宽度为dt；这样，当listmode数据输入后，经过累积运算就会输出h1。

当符合事件j发生时，与它相关的晶体为c1和c2，时间差为tj，则偏移量tjd＝[tj/dt]([]表示对括号里的算术取整操作，取整方式可以是四舍五入，向上或者向下取整)，则分别在h1[c1，m-tjd]和h1[c2，m+tjd]上分别记录，记录的方式采用自加，其中m是指的n的1/2，即分隔数的一半；

校正系数计算模块通过累积计数listmode数据发生的校正信息，获取晶体的时间直方图分布h2，具体包括：

定义晶体的时间直方图分布h2为一个二维的存储空间h2[c,n]，该空间的大小为c*n，c代表的晶体的序号的最大值，n代表了出现在该晶体上的事件所记录的时间差的分隔数，n个分隔数中每一个分隔所对应的时间宽度为变量dt，这样，当listmode数据输入后，经过累积校正信息就会输出h2；

当符合事件j发生时，与它相关的晶体为c1和c2，时间差为tj，则偏移量tjd＝[tj/dt]([]表示对括号里的算术取整操作，取整方式可以是四舍五入，向上或者向下取整)，则分别在h2[c1,m-tjd]和h2[c2,m+tjd]上分别记录，记录的方式是h2[c1,m-tjd]＝h2[c1,m-tjd]+b^addj，h2[c2,m+tjd]＝h2[c2,m+tjd]+b^addj，b^addj是散射符合事件和随机符合事件的计算值，其中m是指的n的1/2，即分隔数的一半；

lor线追踪累积模块根据辐射源的空间分布的体数据计算得到晶体的时间直方图分布l1，具体包括：

定义晶体的时间直方图分布l1为一个二维的存储空间l1[c,n]，该空间的大小为c*n，c代表的晶体的序号的最大值，n代表了出现在该晶体上的事件所记录的时间差的分隔数，n个分隔数中每一个分隔所对应的时间宽度为变量dt；

统计所有从晶体出发的lor线共有lorc条，其中任意一条lor线上，记录lor上任意一个位置的图像像素与lor中心的距离d，得到时间直方图分布l1的时间差td＝d/vc，偏移量tdoffset＝[d/vc/dt]，其中vc表示光速常量，tdoffset代表了以dt作为单位下的距离；

定义射线追踪法获取到晶体序号pix，偏移量tdoffset时l1的强度为

公式中，lorc是与晶体相关的所有lor，xi(t)，yi(t)，zi(t)是晶体pix的第i条lor上时间差为t的三维坐标，可以通过晶体的几何位置快速推导，σ是系统的时间分辨率，t代表时间变量；

s4.晶体时间基准计算模块计算晶体时间基准(即晶体延时)，具体包括：

定义l2[c,n]是一个二维的存储空间，则

其中的sum()是指对数组内所有元素相加，通过卷积，计算第k个晶体的晶体时间基准tdk，

公式中，h为自变量，优化求解最大值算法argmax使用线搜索(linesearch)实现；

s5.迭代更新模块进行迭代更新，具体包括：

设置迭代终止条件；

对晶体时间基准通过td′k＝∈*tdk0+(1-∈)*tdk实现迭代更新，公式中，td’k代表更新后的晶体时间基准，tdk0代表上一次晶体时间基准，∈是经验参数，代表了更新的速度；根据更新后的td’k，选出晶体c1和c2的晶体时间基准td’c1和td’c2，计算晶体c1和c2的时间差tj，tj＝tj0+td’c1-td’c2，tj0为晶体c1和c2中的初始时间差；

为提高计算精确度，本实施例中的时间宽度dt为变量，在第一次迭代过程中，dt一般取很大的值来保证计算范围，在后面的迭代中dt值会逐次减小，最小值是信号采集时间刻度，为方便计算，时间宽度dt的变化遵循dt[ir]＝dt0*p^ir，其中ir为迭代序号，p为[01]之间的分数，dt0为初始的间隔宽度；

每一次更新tj和dt后，返回步骤s3，重新计算h1、h2及l1，使它们动态更新时间计算范围，从而达到更加精确计算的目的；迭代截止条件有不同的设置方法，如规定迭代的次数，或者计算tdk0和tdk之间差异值，差异值足够小时停止迭代；

步骤s6.迭代终止后，利用周知的list-mode-osem重建公式进行图像重建，重建公式如下：可参见文献[mewernerandjskarp,tofpetoffsetcalibrationfromclinicaldata.phys.med.biol.58(2013)4031–4046]。

上述实施例仅用于解释说明本发明的构思，而非对本发明权利保护的限定，凡利用此构思对本发明进行非实质性的改动，均应落入本发明的保护范围。