一种流水线式TOF正电子图像快速重建方法与流程

一种流水线式tof正电子图像快速重建方法
技术领域
[0001]
本发明属于正电子发射断层成像领域，特别涉及了一种正电子图像重建方法。


背景技术：

[0002]
正电子发射断层成像(positron emission tomography,pet)，是利用正电子湮灭时所产生的高能γ光子对物体内部状态进行成像刻画的技术。正电子湮灭产生的γ光子对具有稳定的方向和能量特性，不受电场、磁场、温度等外界环境变化的影响，能穿透密度较大的金属及其它材料，这使得pet技术在工业无损检测方面呈现出独特的优势。利用正电子的湮灭现象，正电子湮灭技术已经成为研究金属、半导体、高温超导体、高聚物等材料物质结构(材料表面及内部缺陷)和电子状态(电子动量分布)的无损伤探测分析手段。正电子湮灭无损检测技术由于其穿透力强、可用工质作为探测源及三维成像等特点，一旦解决其在工业无损检测中的动态测量问题，未来很有可能是解决这些检测问题的最佳方案，前景将十分诱人。
[0003]
随着科学技术的不断进步，γ光子探测器能量分辨率已经从过去的15％
--
20％提升到1％，飞行时间(tof)也从不能检测发展到检测不确定度500ps，进一步地还可达到250ps，γ光子图像的重建质量也在不断提高。但是基于目前已有的技术，要得到清晰的正电子图像，仍然需要几分钟甚至几十分钟以上的数据采集时间，这导致正电子技术在工业应用尤其是动态测量工业应用中遇到了极大的瓶颈。在硬件技术不断取得进步的情况下，缩短扫描时间，提高数据采集速度已经有了一定的支撑条件，但是在正电子图像重建的算法上，还需要新的技术和手段来提高图像的质量。
[0004]
飞行时间pet(time-of-flight pet,tof-pet)相比于传统pet能够利用飞行时间信息大致确定湮灭点发生的位置，从而缩小重建像素的范围。由于重建时将对投影数据有贡献的像素限定在了一定的范围内，tof-pet可以显著地提升图像的对比度以及信噪比，从而提高工业件小缺陷的检出率。同时tof-pet可以在保持图像质量的前提下，降低放射性核素的剂量以及减少采集时间。


技术实现要素：

[0005]
为了解决上述背景技术提到的技术问题，本发明提出了一种流水线式tof正电子图像快速重建方法，使电子发射断层图像重建能便捷、高效、准确地进行。
[0006]
为了实现上述技术目的，本发明的技术方案为：
[0007]
一种流水线式tof正电子图像快速重建方法，包括以下步骤：
[0008]
(1)构建传统系统矩阵a，a中的元素a
ij
表示第j个像素发出的光子被第i条响应线探测到的概率，探测器的不确定度为以湮灭确定点为中心的高斯分布，以像素为单位，计算距该高斯分布中心点的不同距离d，构建高斯表g，g中的元素g
d
表示距离d所对应的像素权值高斯分布；
[0009]
(2)设定时间段δt作为一帧采样时间段，通过pet探测器获取光子对符合数据，从
而对第一帧数据进行采集，在该时间段δt中，根据响应线位置和飞行时间的信息，将湮灭点的坐标位置转换成图像中的坐标(x
t
,y
t
)，并通过查表的方式得到传统系统矩阵a的对应行矩阵a
i
，对行矩阵a
i
进行逐点运算，当元素a
ij
为非零时，根据图像分辨率得到该元素所对应的像素点坐标(x
j
,y
j
)，根据该坐标计算距离d，由计算出的d查询高斯表g，得到对应的g
d
，据此计算出行矩阵a
i
中每个元素对应的f
ij
＝a
ij
·
g
d
，将所有f
ij
组成行矩阵f
i
；
[0010]
(3)上一个时间段δt结束时，下一个时间段δt开始，在此δt内，首先将步骤(2)得到的行矩阵f
i
合成得到系统矩阵f
n
，下标n表示第n帧，将上一帧重建的图像矩阵u
n-1
作为输入，通过图像重建算法得到新的图像矩阵u
n
，并将图像矩阵u
n
作为下一次迭代的输入；与此同时，在本时间段内继续进行步骤(2)的数据采集与处理过程；重复步骤(3)，在每一个时间段δt内持续得到能反映对象动态变化过程的正电子图像。
[0011]
进一步地，在步骤(1)中，利用符合探测技术确定第一帧数据内湮灭事件发生的响应线，将采集到的符合数据转换成包含位置和时间的多维度信息数据，从中获取包含光子飞行时间差的tof数据，同时根据响应线的距离和角度信息，得到相应的传统系统矩阵a。
[0012]
进一步地，在步骤(1)中，g
d
的计算方法如下：
[0013][0014]
其中，e为自然常数，δ为探测器不确定度的标准差，δx为探测器的不确定度，c为光速，δi为系统的时间分辨率。
[0015]
进一步地，在步骤(2)中，距离d的计算方法如下：
[0016][0017]
进一步地，在步骤(3)中，采用ml-em重建算法进行数据迭代，给定初始的估计图像，初始再根据估计图像计算估计投影，将估计投影值与实测投影值进行比较，比较后更新估计图像，重复迭代直至满足迭代终止条件；ml-em重建算法的迭代式如下：
[0018][0019]
其中，和分别为第j个元素的活度在第k+1次迭代和第k次迭代的结果，n为湮灭事件的数目，j为像素数目，y
i
为第i条投影线的实测投影。
[0020]
进一步地，在步骤(3)中，计算y
n
＝f
n
·
u
n-1
，其中y
n
为以f
n
为系统矩阵的正弦图向量，根据该式求出新的图像矩阵u
n
各个像素点的活度值，从而得到新的图像矩阵u
n
。
[0021]
采用上述技术方案带来的有益效果：
[0022]
(1)本发明是针对于上述背景技术提出的正电子发射断层图像重建效率低、速度慢等问题，提出了一种流水线式tof正电子图像快速重建方法，实现切片图像的快速重建，进而实现了边采集边处理的图像动态连续重建，对图像的变化可以进行动态跟踪，为解决正电子技术在工业无损检测中的动态测量问题提供了技术支持。
[0023]
(2)本发明当前的计算结果既可以作为图像的输出，也可以作为下一次迭代计算的图像输入值，进而得到下一次计算的结果，并且通过不断迭代实现动态跟踪图像变化，这不仅可以减少计算工作的工作量，还可以保持迭代的优势。
[0024]
(3)本发明利用飞行时间数据实现连续采样，同时可以在低计数条件下保持图像质量的前提下，降低放射性核素的剂量的同时能有效缩短采集时间，为解决正电子技术在工业无损检测中的动态测量问题提供了技术支持。
[0025]
(4)本发明不需要像传统pet一样计算lor经过的所有像素对该投影值的贡献，而是先根据时间信息确定湮灭事件发生的大概位置，然后根据时间分辨率性能计算不确定度的长度，重建时只计算湮灭事件发生在不确定度范围内的概率，并且能够按照以湮灭确定点为中心的高斯分布对lor的计数进行分配，在图像重建时将对投影数据有贡献的像素限定在了一定的范围内，可以显著地提升图像的对比度以及信噪比，从而提高工业件小缺陷的检出率。
[0026]
(5)本发明利用list-mode数据格式逐个存储湮灭事件所含信息的特点进行逐个事例的重建，充分保留了时间信息，存储的时间信息是连续分布的，利用到了每一个湮灭事件的飞行时间信息，使利用到的信息更为全面，使精度有很大提高，在短扫描时间低计数率情况下，并没增加太多的运算时间，但是增加了更多的有用信息，进而图像质量将有很大提高，在低计数条件下相比于传统算法具有显著优势，这种短时、低采样的工作模式尤其适用于工业领域。
附图说明
[0027]
图1是本发明的流程图；
[0028]
图2是本发明中list-mode数据组成示意图；
[0029]
图3是本发明中飞行时间测量原理图；
[0030]
图4是本发明中tof不确定的高斯分布示意图；
[0031]
图5是本发明中平移计算距离所对应权值g
d
示意图；
[0032]
图6是本发明中基于tof的ml-em算法流程图。
具体实施方式
[0033]
以下将结合附图，对本发明的技术方案进行详细说明。
[0034]
本发明提出了一种流水线式tof正电子图像快速重建方法，如图1所示，通过短扫描时间、多帧位地采样来提高图像重建速度，从而通过更短的扫描时间来实现更高效的图像采集，实现边采集边处理的图像动态连续重建，对图像的变化进行动态跟踪。
[0035]
步骤1、构建传统系统矩阵a，a中的元素a
ij
表示第j个像素发出的光子被第i条响应线(lor)探测到的概率，探测器的不确定度为以湮灭确定点为中心的高斯分布，以像素为单位，计算距该高斯分布中心点的不同距离d，构建高斯表g，g中的元素g
d
表示距离d所对应的像素权值高斯分布。
[0036]
在传统未考虑飞行时间(tof)的pet图像重建过程中，系统矩阵描述了探测器对被检测物的探测过程，它连接了图像空间和投影空间，是迭代重建算法的关键，其公式如下：
[0037]
y＝a
·
u
[0038]
其中，y是向量化的正弦图，u是向量化的切片图像，它建立了正弦图和切片图像的关系。重建过程就是通过投影正弦图数据y和系统矩阵a得到图像信息u。
[0039]
然后，通过pet探测环探测获取响应线，利用符合探测技术确定第一帧数据内湮灭事件发生的lor，并统计得到每条lor上的事件数，将从正电子湮灭γ光子探测采集到的符合数据转换成如图2所示包含光子能量、位置以及入射到晶体条时的信息等多维度信息的list-mode数据。根据lor的距离和角度等位置信息，就可获得相应的系统矩阵。从list-mode数据中获取包含光子飞行时间差的tof数据，利用list-mode数据中的时间和位置信息可进行逐个存储湮灭事件所含信息的特点进行逐个事例的重建，这样就可以利用每一个湮灭事件的飞行时间信息。
[0040]
一个湮灭事件产生的两个γ光子入射到晶体的时间可以被探测到，湮灭发生的位置距离响应线lor的中心越远，对应的两个光子入射到晶体的时间差就越大。飞行时间(tof)能够根据两个γ光子飞行到两个探测器的时间差来确定湮灭点发生在lor上的大致位置。根据该飞行时间差理论上可以计算出湮灭发生的精确位置，但由于系统具有一定的时间分辨率，所能确定的湮灭位置并非一个点，而是以该点为中心的一定范围内，称为不确定度长度，不确定度一般符合高斯分布。tof定位湮灭辐射点位置的原理如图3所示，设d为探测器a到扫描中心的距离，d是某次正电子湮灭辐射点到扫描中心的距离。因此探测器a到正电子湮灭辐射点的距离d
a
＝d+d，探测器b到正电子湮灭辐射点的距离d
b
＝d-d，两个γ光子到达两个探测器的时间差为：
[0041][0042]
其中，c是光速。而湮灭位置发生的不确定度为：
[0043][0044]
其中，δi为系统的时间分辨率。该不确定度是以湮灭确定点o为中心的高斯分布，其长度是高斯分布的半高宽，只有在不确定度范围内的像素才对lor有贡献，其权值符合高斯分布，以像素为单位，计算距高斯分布中心点o距离d，对应的高斯权系数g
d
：
[0045][0046]
g
d
表示计算距离d所对应像素的权值高斯分布，是不确定度范围内像素对lor贡献的高斯权值。由高斯分布特点，不确定度范围内的像素距离湮灭中心点越远，其对lor的贡献越小。δ为不确定度的标准差：
[0047][0048]
如图4所示，当传感器对于时间检测存在误差时，tof可以按照不同的权值g
d
(高斯分布)对lor上的计数值进行分配，将湮灭点的位置限定在一定范围内，为此可以构建一个高斯表g，将带有飞行时间的响应数据信息保存在高斯表中，即根据探测器的不确定度可构建高斯表g。
[0049]
步骤2、设定时间段δt作为一帧采样时间段，通过pet探测器获取光子对符合数据，从而对第一帧数据进行采集，在该时间段δt中，根据响应线位置和飞行时间的信息，将
湮灭点的坐标位置转换成图像中的坐标(x
t
,y
t
)，并通过查表的方式得到传统系统矩阵a的对应行矩阵a
i
，对行矩阵a
i
进行逐点运算，当元素a
ij
为非零时，根据图像分辨率得到该元素所对应的像素点坐标(x
j
,y
j
)，根据该坐标计算距离d，由计算出的d查询高斯表g，得到对应的g
d
，据此计算出行矩阵a
i
中每个元素对应的f
ij
＝a
ij
·
g
d
，将所有f
ij
组成行矩阵f
i
。
[0050]
根据坐标(x
j
,y
j
)计算距离d：
[0051][0052]
根据距离和角度等位置信息在传统系统矩阵中进行查询就可得到相应的元素行矩阵，对所有的lor都进行处理后，就可获得每一根响应线对应的行矩阵a
i
。采样是动态进行的，每次采样完成后得到多个响应线，则相应得到多个行矩阵，每个采样时间段δt内的响应线有不同，故得到的系统矩阵也不同，即系统矩阵(用a
n
表示)也是动态变化的，其行矩阵a
i
的元素也会随着向量y的组成发生改变。数据采集是连续不断地进行的，在采集过程中不断地对每个响应线处理，一旦到达设定时间，响应线处理完成，本帧数据的系统矩阵更新完成。
[0053]
步骤3、上一个时间段δt结束时，下一个时间段δt开始，在此δt内，首先将步骤2得到的行矩阵f
i
合成得到系统矩阵f
n
，下标n表示第n帧，将上一帧重建的图像矩阵u
n-1
作为输入，通过图像重建算法得到新的图像矩阵u
n
，并将图像矩阵u
n
作为下一次迭代的输入；与此同时，在本时间段内继续进行步骤2的数据采集与处理过程；重复步骤3，在每一个时间段δt内持续得到能反映对象动态变化过程的正电子图像。
[0054]
步骤3.1，上一个时间段δt结束时，下一个时间段δt开始。即第一帧的数据采集完成后，马上获取第二帧的符合响应数据并对获取到的每一根响应线进行处理。在第二个δt内，首先将步骤2得到的f
i
合成得到带有tof的系统矩阵f1。
[0055]
步骤3.2，利用ml-em算法进行迭代时，先给定初始的估计图像，初始估计图像设为1，再根据估计图像计算估计投影，将估计投影值与实测投影值进行比较，比较后更新估计图像，如果满足迭代终止条件，则终止迭代，否则重复过程。本发明方法在传统ml-em算法的基础上对每个事件的飞行时间信息加以利用，相应的ml-em重建算法公式如下：
[0056][0057]
其中，和分别为第j个元素的活度在第k+1次迭代和第k次迭代的结果，n为湮灭事件的数目，j为像素数目，y
i
为第i条投影线的实测投影。按上式进行迭代重建，其中每次重建迭代次数在总迭代时间不超过采样时间δt，即不影响下一次的采样数据计算情况下，尽可能多的迭代，以使收敛性更佳。
[0058]
步骤3.3，多帧数据引入了图像之间的不同高频细节信息，对于多数据帧的采样方式，当探测器在原始位置进行探测时，得到的数据用正弦图向量表示为y1，如果用ml-em算法进行图像重建，y1和切片图向量u的关系是：
[0059]
y1＝f1·
u0[0060]
其中，f1是y1对应的系统矩阵，初始图像u0设为1。根据响应线的飞行时间(tof)信息，计算距离d，对构造好的高斯表g进行查表，然后与由位置信息查表得到的系统矩阵a1的
行矩阵a
i
按行逐点相乘，可得到矩阵f1。
[0061]
根据上述求出待求图像u1中各个像素点的活度值，从而得到待求图像u1。
[0062]
步骤3.4，在第二帧响应事件采样完成后，在下一个δt内马上对第三帧的符合响应数据进行采样，同时在本时间段内继续步骤2数据处理工作得到a2，其中，计算距离d并查表g所对应的权值g
d
，如图5所示，可通过平移得到，进而获得相应的tof系统矩阵f2。整个过程并不需要完全重算，大大节省了计算时间。
[0063]
由于每个采样时间段内的响应线有不同，所以即tof系统矩阵也是动态变化的，tof系统矩阵f1≠f2，tof系统矩阵的元素也会随着向量y的组成发生改变。
[0064]
然后按步骤3.2的方法对采集到的第二帧数据进行ml-em迭代重建，本次迭代是将上一帧重建的图像矩阵作为输入，通过图像重建算法得到新的图像矩阵，该图像矩阵可以作为图像输出，并作为下一次迭代的图像输入，即对u1图像进行更新获得新的图像u2。
[0065]
这时得到的数据用正弦图向量表示为y2，用ml-em算法，y2和切片图向量u1的关系是：
[0066]
y2＝f2·
u1[0067]
其中，f2是y2对应的系统矩阵，这里u1是上一次迭代得到的图像。根据上式求出待求图像u2中各个像素点的活度值，从而得到待求图像u2。
[0068]
步骤3.5，以此类推，重复步骤3，直到最后一帧数据采集完成，按以上步骤方法对最后一帧数据进行处理，更新得到新系统矩阵f
n
，进行迭代重建，对上一切片图更新得到最终的切片图，图像重建完成。
[0069]
最后一帧数据采集完成后，在本时间段内按步骤2的方法进行数据处理获得a
n
和所对应的g
d
，进而得到tof系统矩阵f
n
，然后将u
n-1
图像矩阵作为输入，按步骤3.2的方法对数据进行ml-em迭代重建，对u
n-1
图像进行更新获得新的图像u
n
，即最终的正电子图像。
[0070]
最后一帧数据用正弦图向量表示为y
n
，用ml-em算法，y
n
和切片图向量u
n-1
的关系是：
[0071]
y
n
＝f
n
·
u
n-1
[0072]
其中，f
n
是y
n
对应的系统矩阵，这里u
n-1
是上一次迭代得到的图像。根据上式求出待求图像u
n
中各个像素点的活度值，从而得到待求图像u
n
，即最终的正电子图像，图像重建完成，并在每一个时间段δt内持续得到能反映对象动态变化过程的正电子图像，如图6所示，m为当前迭代次数，i为终止迭代次数。
[0073]
实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。