本发明涉及核医学pet图像重建领域，尤其涉及一种基于构造虚拟doi及相应系统矩阵的pet图像重建与优化方法。
背景技术：
pet(positronemissiontomography,正电子发射计算机断层扫描)是核医学领域重要的临床检查影像技术。pet的探测原理是将能发射正电子的核素作为特异性示踪剂注入病人体内，示踪剂发射的正电子与人体内的电子湮没产生一对方向相反的γ光子，这对光子被探测装置记录并储存下来，通过图像重建就可以得到示踪剂在病人体内的分布情况，从而实现肿瘤的早期发现和定位。
pet设备中的探测装置一般由若干个闪烁晶体与光子探测器组成的近似环形结构组成，可以探测各种角度发射的光子对。具体的工作原理是：γ光子在闪烁晶体中衰减一段距离，被转换为可见光信号，然后被光子探测器捕获，由光电倍增管转换为电信号记录，探测器每记录下的一次信号被称为一次lor事件(lineofresponse)。对信号传统的处理方式是：忽略探测器的尺寸，将一对探测器(互成180°)捕捉到的lor事件归纳为同一类事件，即光子的发射位置处于这一对探测器的中心连线某点上。这种方法可以大大提高探测器的计数效率，并且计算简单方便，但同时也忽略了光子在闪烁晶体内被探测到时的位置深度，会造成计算光子发射位置时存在很大误差，造成图像的空间分辨率下降。而pet图像的分辨率过低也正是pet技术发展中的一大诟病。如图1所示：传统pet成像设备将某一次由源a所产生的lor事件，归纳为由这对探测器p1p2的中心连线上某点产生的，计算出源处于p1p2连线上某点，而实际上，该次lor事件的源处于a点，故经过图像重建后会造成肿瘤位置的偏差。
为了解决这一问题，可以在物理硬件上重新设计制造出能检测光子在晶体中的作用深度doi(depthofinteraction)的探测器，但物理方法有几个致命的缺点：一是目前探测器的制作工艺困难，太薄难以捕获光子，太厚又无意义；二是会使探测器更复杂导致成本更高；三是很多已有的pet设备不具备升级硬件的条件，不仅需要重新匹配pet系列中的其他模块，还需要重新测试系统性能。
技术实现要素：
本发明所要解决的技术问题是提供一种构造虚拟doi及相应系统矩阵来进行pet图像重建的新方法。首先对晶体虚拟创建出doi，再对每一对doi检测到γ光子的概率进行计算，增大了检测到pet数据(lor事件)的几何范围；并根据虚拟doi构造出相应的系统矩阵用以后续的pet图像重建迭代算法，这个系统矩阵将不同于以往的简便计算方法，考虑了更多的影响图像质量的因素，计算精度将大大提高；最后在云计算和图形处理单元(gpu)平台上，通过图像重建迭代算法对数据进行分布式图像重建计算，降低了计算时间。具体的技术细节如下：
步骤1，虚拟doi的划分，将pet探测器环中的闪烁晶体划分成若干个区域并对应不同的作用深度；
步骤2，计算虚拟doi的探测概率，将虚拟划分出的doi区域进行配对，计算每一对在不同信号源位置被检测到的概率；
步骤3，划分虚拟doi中的虚拟lor事件，将原先一对晶体中的一条lor事件在虚拟doi中划分为多个子事件，并计算每个虚拟子事件的发生概率，将其作为后续重建图像时每条lor子事件的权重系数；
步骤4，对pet重建中的常规系统矩阵进行优化，将步骤2计算获得的虚拟doi对应的探测概率加入到系统矩阵；
步骤5，利用步骤4获得的系统矩阵对图像进行重建，获得优化后的图像。
进一步的，步骤1中根据闪烁晶体的几何排列和尺寸大小进行虚拟doi的划分。
进一步的，步骤2中计算虚拟doi的探测概率的实现方式如下，
通过连接信号源与某对虚拟doi中心的连线，计算覆盖在晶体中的长度，确定检测概率d，对于一个光子，它在晶体内被检测概率d与其在之前相邻晶体中穿透过的距离a和最终某个晶体内部穿越的距离l相关：d＝(1-e-μl)e-μa，其中μ为γ光子在某种晶体中的衰减系数，为常数；
假设源m处发生的一次lor事件，由一对探测器p1和p2探测到，γ光子在探测器p1中先穿透距离a1，最后才在某个晶体内穿越距离l1后被捕获，将a1和l1带入公式可以计算出光子在晶体p1内被检测到的概率d1，同样的，计算出光子在晶体p2内被检测到的概率d2，则光子在这对探测器中被检测到的概率d＝d1d2。
进一步的，步骤4中优化后的系统矩阵的计算公式如下，
pij＝pdet·blu，p，eom
其中pdet.blur代表光子在晶体不同区域检测到的概率，由步骤2获得；pgeom代表光子通过像素点的几何概率，由像素j到光子投影线i的几何距离的远近决定。
进一步的，步骤5中利用图像重建迭代算法osem对图像进行重建。
进一步的，步骤4和步骤5在gpu平台上进行分布运算。
与现有技术相比，本发明的优点和有益效果如下：1)能有效的提高重建后的pet图像的分辨率2)极大降低了制造成本，本发明通过数学方法构造虚拟doi，在提供相当性能的同时，能极大降低设备的制造成本3)本发明所示的方法流程简单方便，适用范围广，可以轻松应用在各种pet设备上，均能提高相应的图像质量。适用范围：构造虚拟doi所适用的pet设备，不局限于设备的结构，可以是全身，局部器官或小动物的设备；不局限于晶体与探测器的几何分布，可以是环形，半环形，平板或多边形分布。
附图说明
图1为传统pet探测器处理lor事件忽略doi产生误差示意图；
图2为本发明实施例中pet探测器的虚拟doi划分示意图；
图3为本发明实施例中光子在闪烁晶体不同位置(doi)被检测到的概率计算示意图；
图4为本发明具体实施的流程图；
具体实施方式
下面结合附图对本发明内容作进一步介绍，本发明实施例的具体实现步骤如下：
步骤1，虚拟doi的划分：将pet探测器环中的闪烁晶体划分成若干个区域并对应不同的作用深度(doi)，以增大探测器探测到的几何范围。
区别于可以检测到实际doi的物理探测器，这些区域是虚拟的。根据探测器的角度和覆盖的几何范围，虚拟doi的个数可以自由调整，可由pet探测器的几何分布与晶体大小决定划分的数量。例如探测器本身所覆盖的几何范围很小，可以划分很少的虚拟doi。而对于边缘部分的晶体(需要更精确的探测更大范围内的γ光子)，可以大量划分虚拟doi。如图2所示，将一对互成180°的探测器p1,p2虚拟划分为p1,1；p1,2；p2,1；p2,2这4个子探测器，对原先探测器仅能计算出源分布在p1p2连线上，通过划分虚拟doi后计算出源分布可以在p1,1p2,1；p1,1p2,2；p1,2p2,1；p1,2p2,2这4条连线上，分辨率较原先有了极大提高。
优化虚拟doi的划分：不同于物理硬件上根据晶体深度方向划分doi，虚拟doi的划分不局限在晶体的深度方向上，可以取决于不同的几何形态，为了模拟所有γ光子可能被检测到的几何位置，虚拟doi可以划分在晶体内的其他方向，例如细分在晶体长度与宽度方向上，如果晶体非长方体，也可根据不同几何形态划分，这些是物理探测器达不到的。
步骤2，计算虚拟doi的探测概率：将虚拟划分出的doi区域进行配对，计算每一对在不同信号源位置被检测到的概率。
以图2所示为例，即将一对探测器p1p2实际所采集到的所有lor事件，分解到p1,1p2,1；p1,1p2,2；p1,2p2,1；p1,2p2,2这4对虚拟探测器上。具体计算方法为：通过连接信号源与某对虚拟doi中心的连线，计算覆盖在晶体中的长度，确定探测概率：对于一个光子，它在晶体内被检测概率d主要与其在之前相邻晶体中穿透过的距离a和最终某个晶体内部穿越的距离l相关：d＝(1-e-μl)e-μa，其中μ为γ光子在某种晶体中的衰减系数，为常数，可查表得知。如图3所示，源m处发生的一次lor事件，由一对探测器p1和p2探测到，γ光子在探测器p1中先穿透距离a1，最后才在某个晶体内穿越距离l1后被捕获，将a1和l1带入公式可以计算出光子在晶体p1内被检测到的概率d1，同样的可以计算出光子在晶体p2内被检测到的概率d2，则光子在这对探测器中被检测到的概率d＝d1d2。
对于一些边缘位置，还可以通过蒙特卡洛数学模拟的方法来获取更精确的探测概率。需要注意的是，探测概率也由一些其他因素决定，如光子的散射，病人器官不同位置对光子的衰减，光子湮灭不对称，正电子偏移范围等。因此可以根据对所成图像的具体需求，将这些因素加入到探测概率的计算，或者可以通过蒙特卡洛计算机模拟的方法获取虚拟doi的探测概率。
步骤3，虚拟doi中的虚拟lor事件：将原先一对晶体中的一条lor事件在虚拟doi中划分为多个子事件，并计算每个虚拟子事件的发生概率，以此作为后续图像重建时的权重系数。
对于一组互成180°的探测器，区别于传统pet图像重建处理一次lor事件，假设将此组探测器的闪烁晶体分别分为m个与n个虚拟doi，则将考虑不同虚拟doi的组合即需要处理m×n个虚拟lor子事件，每个lor子事件发生的概率由计算出的每对虚拟doi的探测概率决定：将所有虚拟doi配对组对应的探测概率求和，做归一化处理，使得其总的探测概率为1，即可求得每一对虚拟doi对应的虚拟lor子事件在当前信号源位置发生的概率。在后续的图像重建中，每计算一条虚拟lor子事件，则将该子事件相应的发生概率作为权重系数计算在内。
以图2所示虚拟doi为例，即将原先实际采集的1条p1p2连线上发生的lor事件，分解为了p1,1p2,1；p1,1p2,2；p1,2p2,1；p1,2p2,2连线上的4条lor子事件，而这4条中每条子事件发生的概率，由步骤2计算出的每对虚拟doi的探测概率决定，例如p1,1p2,1子事件发生的概率为：
步骤4，考虑光子在闪烁晶体不同位置的检测概率对图像重建质量的影响，对常规系统矩阵进行优化；
系统矩阵p(systemmatrix)是数学上图像空间与探测器空间的变换矩阵，图像重建迭代算法中最重要的参数，在图像重建过程中为了获取更准确的图像，需要获取高精度的系统矩阵。其通常由多种因素构成，可以描述成:p＝pdet.senpdet.blurpattnpgeomppositron。其中pdet.sens代表探测器均一化矩阵，pdet.blur代表光子在晶体不同区域检测到的概率，patt用于衰减修正，pgeom代表光子通过像素点的几何概率，ppositron考虑了正电子偏移范围。为了简化计算，在传统pet设备的图像重建过程中通常考虑的是占决定性因素的pgeom，近似认为p＝pgeom。pgeom由像素j到光子投影线i的几何距离dij的远近决定，
其中，σ为常数，dij表示像素j到光子投影线i的几何距离；pgeom的计算方法在pet专业书籍文献中均有详细介绍，为现有技术。
而光子在闪烁晶体不同位置(doi)的检测概率pdet.blur通常被忽略。本发明通过构造虚拟doi，再将pdet.blur考虑进系统矩阵。以图3为例，将光子在每个闪烁晶体对中被检测的概率计算后可以得到：pdet.blur＝d1d2。最终可以得到系统矩阵pij＝pdet.blurpgeom。
步骤5，利用步骤4获得的系统矩阵对图像重建，将系统矩阵pij带入pet图像重建迭代算法osem中求得优化后的图像，。
本发明利用gpu进行分布运算：通过增加了虚拟lor事件的数量(m×n倍)能大大提高pet图像的质量，但也会造成数据计算量的提高，导致计算时间的增加。利用gpu进行分布计算可降低数据处理的时间。由于在gpu上内存读写速度很快，首先将pet源数据(listmode格式，记录每一次lor事件)缓存到gpu的全局内存中。然后分配gpu中的每一个线程并行计算一条lor(子)事件。为每个线程分配共享内存以存储重建中的图像，将二维图像空间初始值存入gpu的共享内存中。每个线程直接计算系统矩阵并根据osem公式对图像计算正向与反向投影，经过若干次迭代计算，最终可以得到重建后的图像。由于共享内存具有很大的带宽，所有的事件可以同时计算，并且图像存储在共享内存中可以提高数据的读写速度。对于nvidiagpu利用cuda或者opencl工具运算，对于amdgpu则利用opencl工具运算。由于gpu可以划分数以千计的线程，因新增加的虚拟doi区域和虚拟lor事件导致增大的运算量可以通过增大线程数量解决。对于数据量更大的系统，还可以通过在云平台上建立gpu集群计算。
具体实施时，本发明可通过计算机流程自动实现。
1.对于一个pet系统，首先获取该系统的配置(如探测器和闪烁晶体的几何排列，晶体类型，尺寸大小等物理参数)，用以建立listmode格式数据所对应的图像空间坐标。
2.根据系统配置划分虚拟doi。为了简化计算，通常默认虚拟doi的划分方法：即沿晶体深度方向划分虚拟区域
3.将采集得到的listmode格式数据导入到计算机gpu的全局内存中，以利用gpu内存读写的高带宽减少数据传输损耗时间。并对listmode格式数据用更小的位数进行重新编码，以节省内存存入更多数据。
4.根据虚拟doi的划分将lor事件划分为若干倍的虚拟子事件。
5.计算某个lor事件中的γ光子在某个晶体中的衰减距离a和反应距离l。
6.计算划分的虚拟doi的检测概率d，并计算相应的虚拟lor子事件的发生概率。
7.遍历每一条lor子事件，计算出传统方法中pgeom对应的系统矩阵。
8.计算虚拟事件在晶体中被检测的概率pdet.blur。最终求得优化后的系统矩阵pij。
9.将计算出系统矩阵代入osem迭代算法重建出图像。
10.整个系统矩阵和图像的计算都是在gpu平台上进行的，计算速度大大提高。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。