基于多边形视野对称性的旋转双平板PET系统及其成像方法

摘要

本发明公开了一种基于多边形视野对称性的旋转双平板PET系统，可以获取多角度投影数据，解决了平板PET系统无法获得完整角度投影数据的问题。本发明还公开了一种基于该系统的成像方法，包括以下步骤：通过数据采集模块获取完整角度投影数据；划分视野、编码晶体和响应线，确定计算所需的体素和相应线编码；计算多边形视野系统响应矩阵利用系统响应矩阵和投影数据进行迭代重建；进行多面体体素到正方体体素的转化，并进行图像显示。本发明提出了利用多边形视野结构的对称性，打破了正方形视野结构对旋转系统对称性的限制，减小了系统矩阵的计算量，提高了PET系统的重建效率；通过调整视野中扇区的大小，本发明可以广泛的应用于旋转平板PET系统中。

说明书

技术领域

本发明涉及核医学影像成像领域，具体涉及一种基于多边形视野 对称性的旋转双平板PET系统及其成像方法，可应用于小动物在体成 像。

背景技术

PET(Positron Emission Tomography,PET)全称为正电子发射 计算机断层扫描，是核医学领域先进的影像技术。PET技术采用解剖 形态方式进行功能、代谢和受体显像技术，尤其在肿瘤的早期检测方 面占据不可忽视的地位。传统的PET探测器大多采用静态环形模式， 虽然可以获得完整角度的角度的采样，但系统的复杂度高，代价大， 同时环形结构封闭性强，给采样对象造成不必要的精神压力，同时不 便于医生或实验员对采样对象的技术指导；为了降低成本，增加探测 器的灵活度，有研究者提出平板PET探测器，由于无法获得完整角度 的投影数据，同时，存在严重的反应深度响应，对成像效果有所影响， 例如Chien-Min Kao等人在“Evaluation of 3D Image Reconstruction  Methods for a Dual-Head Small-Animal PET Scanner”文章中提出 一种静态的双头平板PET探测器，在降低了系统结构的复杂度，但是 无法获得完整角度的投影数据；

系统响应矩阵(System Response Matrix，SRM)是联系探测器 和示踪剂的桥梁，作为探测器系统成像过程的一种数学描述，其计算 主要有实验方法，解析方法和蒙特卡罗仿真方法，实验方法对实验条 件的要求较高，而且操作复杂，蒙特卡罗仿真方法对噪声的抑制能力 比较低，容易受到噪声的影响，而且仿真时间比较长，解析方法虽然 无法精确的描述示踪剂衰变的物理过程，但计算时间较少，而且方便 实现；系统响应矩阵的计算作为PET重建的一个关键因素，其规模往 往非常大，很难保证对其进行实时处理，其计算、存储和调用效率成 为研究者研究的一个主要方向；为了减少系统响应矩阵的复杂度，有 研究者提出利用探测器系统的对称性，减少计算响应线的条数或者减 少计算体素的数目。清华大学在其申请专利“一种用于双平板PET探 测器的图像重建方法”(专利申请号201310055434.7，公开号CN 103099637A)中公开了基于响应线的静态双平板PET对称性方法， 一定程度上简化了系统响应矩阵的计算，提高了系统响应矩阵的压缩 比，但在仿真时需要将体素进行扩展，增加了仿真过程的复杂度，因 此在一定程度上限制了系统响应矩阵的计算效率，而且这种方法只适 用于静态平板PET，对于旋转系统不再具有适用性；June等 人在“Efficient methodologies for system matrix modelling in  iterative image reconstruction for rotating high-resolution  PET”文章中提出利用系统的轴向对称性和面间对称性提高系统矩阵 的计算效率，鉴于其正方体像素结构，限制了旋转系统对称性的发挥。

发明内容

针对现有技术的不足，本发明旨在提供一种基于多边形视野对称 性的旋转双平板PET系统及其成像方法，通过采用多边形视野结构， 打破正方形视野结构对旋转系统对称性的限制，缩减需要计算的像素 数、提高重建效率；

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

基于多边形视野对称性的旋转双平板PET系统包括数据采集模 块、控制模块、成像模块和传输模块，所述数据采集模块、成像模块 和传输模块均与所述控制模块相连接，所述数据采集模块同时还与所 述成像模块相互连接，其中：

数据采集模块：负责进行数据采集，包括双平板PET探测器，所 述双平板PET探测器由两个相距一定距离的相对摆放的PET探测板组 成，两个探测板绕中心轴旋转；每个PET探测板由一定数目且大小相 等的晶体组成，晶体接收核素衰变过程中产生的光子；

控制模块：由计算机系统控制完成，包括系统控制子模块和数据 处理子模块，所述系统控制子模块负责完成对数据采集模块采样间 隔、采样时刻、旋转速率等方面的控制，所述数据处理子模块负责处 理采集到的数据进行信号增强、滤波、衰减校正；

成像模块：由多核计算机系统完成，主要负责进行核素重建，获 取核素的分布情况，其中包括负责多线程处理数据的系统矩阵处理子 模块以及完成对体素的核素重建以及体素的转化工作的重建子模块；

进一步地，所述系统矩阵处理子模块包括如下四个子模块：

视野划分子模块：用来对成像视野空间进行多面体体素划分，划 分视野边界保持与数据采集模块的双平板PET探测器平行；

晶体编码子模块：用来对数据采集模块旋转过程中晶体所在位置 的编码，将晶体编码为一个虚拟晶体环；

系统矩阵计算子模块：利用类线追踪方式，完成对重建过程中系 统响应矩阵的生成；

对称性处理子模块：用于对视野空间和探测器结构对称性的处 理；

进一步地，所述重建子模块包括体素重建子模块和转化子模块两 个子模块：

体素重建子模块：利用计算所得系统响应矩阵，利用迭代MLEM 算法完成对划分多面体体素的重建工作；

转化子模块：完成由多面体体素到正方体体素的转化，获取核素 分布；

传输模块：由帽式滑环组成，用于数据采集模块与计算机处理系 统的数据传输。

需要说明的是，组成所述数据采集模块的双平板PET探测器的初 始位置保持与水平地面平行，在空间旋转时，所述两个PET探测板之 间的间距保持不变并且始终保持平行。

一种基于上述系统的成像方法包括如下步骤：

步骤1，采集数据

利用所述数据采集模块进行数据采集，采集的过程由所述控制模 块的系统控制子模块对所述数据采集模块的采样间隔、采样时刻和旋 转速率进行控制；

步骤2，数据预处理

将采集到的数据通过所述传输模块传送至所述控制模块的数据 处理子模块，完成对采集到的数据进行信号增强、滤波和衰减校正等 操作；

步骤3，划分视野空间

利用所述视野划分子模块，保持视野空间的多边形结构，并按层、 扇区、环、体素逐步细化的方式，将视野步步细化为截面为三角形的 多面体体素；

步骤4，编码晶体

通过所述晶体编码子模块，对所述数据采集模块旋转过程中晶体 所在位置的编码，将晶体编码为一个虚拟晶体环；将连接晶体的连线 作为响应线，响应线采用(晶体i，晶体j)的方式进行调用，其起 始点设在晶体的中心位置，若晶体j和晶体i不在两条相对边界上， 进行剔除，只对连接相对边界晶体的响应线进行保存处理；

步骤5，通过所述系统矩阵计算子模块，利用类线追踪方式，完 成对重建过程中系统响应矩阵的计算；

步骤5.1，确定响应线穿过扇区的每个扇区环的交点坐标 h_i＝(x_i,y_i,z_i)，i表示环编码；通过线面交点确定此条响应线穿过本环 的体素数，并确定穿过第一体素的长度；

步骤5.2，利用交点坐标h，确定响应线穿过每层的第一个体素 索引为：

index＝t+(h-1)²+(z-1)×H²×S

其中，t为此体素在所在环的位置，h为体素所在环，z为体素所 在层，H为环总数，S为扇区总数；

步骤5.3，确定响应线穿过本环其他体素的长度；通过判断体素 是否在本环，若在，利用点面间的距离公式计算两者的交点坐标，以 步骤5.2的方式获取交点索引值，并计两交点长度，作为系统响应矩 阵元素；否则，结束计算；

步骤6，通过所述对称性处理子模块，利用视野空间和双平板PET 探测器结构对称性，简化系统响应矩阵的计算；

步骤6.1，利用扇区间对称性，确定计算所需体素数和响应线：

只需确定一个扇区的体素，选取一个扇区中体素作为计算体素， 记所述计算体素数为M；针对确定计算体素，利用数据采集模块旋转 过程的对称性，只需确定在起始位置时的晶体构成的响应线，记响应 线条数为N；

步骤6.2,利用扇区内对称性，进一步减少计算响应线数：

对于平行响应线满足穿过体素的上截面的位置相同时，只计算穿 过扇区最长的响应线与相交体素的长度；对于与之平行的响应线，只 计算穿过扇区的体素数目即可。

步骤7，重建体素

在所述体素重建子模块内，利用计算所得系统响应矩阵，通过迭 代算法MLEM重建算法重建出多面体体素的体素值，对划分多面体体 素进行重建；

步骤8，体素转化，获取成像结果

由所述转化子模块完成由多面体体素到正方体体素的转化，获取 体内核素分布，以显示核素成像结果。

需要说明的是，所述步骤3的具体实施方法如下：

顺时针编码扇区，记扇区总数为S，由里及外等高度编码环数， 记环总数为H，每一环划分为不同的三角体素，三角截面为等腰三角 形，正三角体素和倒三角体素交叉排列，组成多边形视野空间；每个 环上的三角体素数与环总数H保持正比例增加，呈2×H-1的关系。

需要说明的是，步骤4中对晶体编码的具体实施方法是，以数据 采集模块起始位置为起点，将停留时刻的晶体进行统一编码，虚拟成 一个晶体环，晶体编码以顺时针方向依次递增，两条相对的边界的晶 体编码相差n×m，其中n为晶体采样次数，m为每条边的晶体数。

需要进一步说明的是，所述步骤6.1利用扇区间对称性，确定计 算所需体素数和响应线，包括具体步骤如下：

步骤6.1.1，确定与第一个扇区体素index₁相对应的第i个扇区的 体素index_i：

由层c，扇区i，环h，体素v等体素索引方式，对于第一个扇区 的体素标记为(c,1,h,v)，则第i个扇区上的体素标记为(c₁,i,h₁,v₁)，若 第i个扇区的体素与第一个扇区相对应，则需满足：

c＝c₁

h＝h₁

v＝v₁

index_i＝v₁+(h₁-1)²+H²×(i-1)+c×C_num

index₁＝v+(h-1)²+c×C_num；

其中，C_num表示每层的体素数；

步骤6.1.2，确定体素index_i与体素index₁对应的响应线编码：

设定穿过体素index₁的响应线l为(晶体i'，晶体j'),计算晶体i' 和晶体j'在探测板晶体环上的所在边晶体位置：

c_i'＝i'/CR_num

c_j'＝(j'+n×m)/CR_num；

其中/表示取余操作，CR_num表示每个扇区的体素总数；

设定穿过体素index_i的响应线l_i为(晶体ii,晶体jj),计算晶体ii与 晶体jj在探测板晶体环上所在边的晶体位置：

c_ii＝(ii+(i-1)×n)/CR_num

c_jj＝(jj+(i-1)×n)/CR_num；

其中/表示取余操作，CR_num表示每个扇区的体素总数；

如果满足条件

c_i'＝c_ii

c_j'＝c_jj；

则系统响应矩阵为：

Q(l_i,index_i)＝Q(l,index₁)；

需要说明的是，步骤8体素转化，获取成像结果的具体实施流程 如下：

步骤8.1，令多面体体素的厚度和正方体体素的边界长度相等， 记正方体截面面积为SF，正方体体素数为FC_num，多面体体素数为 VC_num；

步骤8.2，计算体素的重叠体积RP，确定转化矩阵T，T的行列 数分别为VC_num,FC_num；T(i,j)表示第i个三角体素与第j个正方体体 素的截面的重叠面积和正方体体素截面积的比值，即

步骤8.3，利用转化矩阵T求解出每层上正方体体素的体素值：

FC_i＝T(j,i)×VC_j；

步骤8.4，对每层重复步骤步骤8.3，获取最终的正方体体素划 分体素值。

本发明的有益效果在于：

第一、采用多边形视野结构，打破了正方形视野结构对旋转系统 对称性的限制，充分利用扇区间对称性和扇区内对称性，大幅度的缩 减了所需计算的像素数，以此为基础，系统响应矩阵的计算也大幅度 的降低，因而可以在不减弱PET重建质量的同时，缩短PET重建时间， 提高重建效率；

第二、采用旋转系统结构，在保留平板PET成本低，灵活性强等 优点的同时，可以获取完整角度的投影数据，改善了平板PET重建效 果差等问题；

第三、采用三角体素进行像素划分以及采用解析方法获取系统响 应矩阵，提高了系统响应矩阵的计算效率，加快了图像重建速率；

第四、特殊的视野结构不需要对平板探测器晶体或者体素数进行 扩展，相比于现有技术更方便、快捷。

附图说明

图1为基于多边形视野对称性的旋转双平板PET系统示意图；

图2为基于多边形视野对称性的旋转双平板PET系统成像方法的 实施流程图；

图3为图1中数据采集模块1的位置结构图；

图4为旋转双平板PET系统视野划分图，其中左图为视野扇区编 码分布图，右图为扇区内像素编码分布图；

图5为晶体结构二维编码图；

图6为基于多边形视野对称性的旋转双平板PET系统对称性分析 图；

图7为多边形视野到四边形视野过程转化示意图。

具体实施方式

以下将结合附图对本发明作进一步的描述，需要说明的是，本实 施例以本技术方案为前提，给出了详细的实施方式和实施步骤，但并 不限于本实施例。

如图1所示，一种基于多边形视野对称性的旋转双平板PET系统 包括数据采集模块1、控制模块2、成像模块3和传输模块4，所述 数据采集模块1、成像模块3和传输模块4均与所述控制模块2相连 接，所述数据采集模块1同时还与所述成像模块3相互连接，其中：

数据采集模块1：负责进行数据采集，其中包括双平板PET探测 器101；所述双平板PET探测器101由两个相距一定距离的相对摆放 的PET探测板组成，两个探测板绕中心轴旋转；每个PET探测板由一 定数目且大小相等的晶体组成，晶体接收核素衰变过程中产生的光 子；如图3所示，组成所述数据采集模块1的双平板PET探测器101 的初始位置T201保持与水平地面平行，在空间旋转时，所述两个PET 探测板之间的间距保持不变并且始终保持平行。

控制模块2：由计算机系统控制完成，包括系统控制子模块201 和数据处理子模块202，所述系统控制子模块201负责完成对数据采 集模块1的控制，所述数据处理子模块202负责处理采集到的数据；

成像模块3：由多核计算机系统完成，主要负责进行核素重建， 获取核素的分布情况，包括负责多线程处理数据的系统矩阵处理子模 块31以及完成对体素的核素重建以及体素的转化工作的重建子模块 32；

其中，所述系统矩阵处理子模块31由视野划分子模块311、晶 体编码子模块312、系统矩阵计算子模块313和对称性处理子模块314 构成：

视野划分子模块311用于对成像视野空间进行多面体体素划分， 划分视野边界保持与数据采集模块1的双平板PET探测器101平行；

晶体编码子模块312用于对所述数据采集模块1旋转过程中晶体 所在位置的编码，将晶体编码为一个虚拟晶体环；

系统矩阵计算子模块313负责利用类线追踪方式，完成对重建过 程中系统响应矩阵的生成；

对称性处理子模块314用于对视野空间和双平板PET探测器101 结构对称性的处理；

所述重建子模32块则由体素重建子模块321和转化子模块322 构成，体素重建子模块321利用计算所得系统响应矩阵，通过MLEM 迭代算法完成对多面体体素的重建工作；转化子模块322负责完成由 多面体体素到正方体体素的转化，获取核素分布；

传输模块4：由帽式滑环组成，用于数据采集模块1与控制模块 2的数据传输。

如图2所示，基于上述系统，一种基于多边形视野对称性的旋转 双平板PET系统成像方法包括如下步骤：

步骤1，采集数据

利用所述数据采集模块1进行数据采集，采集的过程由所述控制 模块2的系统控制子模块201对所述数据采集模块1的采样间隔、采 样时刻和旋转速率进行控制；

步骤2，数据预处理

将采集到的数据通过所述传输模块4传送至所述控制模块2的数 据处理子模块202，所述数据处理子模块202对采集到的数据进行信 号增强、滤波和衰减校正的操作；

步骤3，划分视野空间

利用所述视野划分子模块311，保持视野空间的多边形结构，并 按层、扇区、环、体素逐步细化的方式，将视野步步细化为截面为三 角形的多面体体素；

进一步地，所述步骤3的具体实施流程如下：

如图4所示，顺时针编码扇区，扇区总数为S，分别为扇区编码 为S1-S8，由里及外等高度编码环H301，环总数为H，每一环划分为 不同的三角体素，三角截面为等腰三角形，正三角体素V301和倒三 角体素V302交叉排列；每个环上的体素数与环数H保持正比例增加， 呈2×H-1的关系；

步骤4，编码晶体

通过所述晶体编码子模块，对数据采集模块旋转过程中晶体所在 位置的编码，将晶体编码为一个虚拟晶体环；将连接探测器晶体的连 线作为响应线，起始点设在晶体的中心位置，响应线采用(晶体i， 晶体j)的方式进行调用；晶体i(C401)和晶体j(C402)不在两条 相对边界上，故对响应线L401进行剔除，只对连接相对边界晶体的 响应线L402进行保存处理；

其中，将晶体编码为一个虚拟晶体环的具体方法为：

如图5所示，利用所述晶体编码子模块312，以对数据采集模块 初始位置T201为起点，将停留时刻的晶体进行统一编码，虚拟成一 个晶体环403，晶体编码以顺时针方向依次递增，如图5中将晶体按 顺时针编码为1-14、15-28、29-42、43-56、57-70、71-84、85-98、 99-112；两条相对的边界(401、402)的晶体编码相差n×m，其中n 为晶体采样次数，m为每条边的晶体数；

步骤5，通过所述系统矩阵计算子模块313，利用类线追踪方式， 完成对重建过程中系统响应矩阵的计算；所述步骤5具体包括如下步 骤：

步骤5.1，确定响应线穿过扇区的每个扇区环的交点坐标； h_i＝(x_i,y_i,z_i)，i表示环编码，通过线面交点确定此条响应线穿过本体 素环的体素数，并确定穿过第一体素的长度；

步骤5.2，利用交点坐标h，确定响应线穿过每层的第一个体素 索引；

index＝t+(h-1)²+(z-1)×H²×S；

其中，t为此体素在所在环的位置，h为体素所在环，z为体素所 在层；

步骤5.3，确定响应线穿过本环其他体素的长度；通过判断体素 是否在本环，若在，利用点面间的距离公式计算两者的交点坐标，以 步骤5.2的方式获取交点索引值，并计算两交点长度，作为系统响应 矩阵元素；否则，结束计算。

步骤6，如图4、图5、图6所示，通过所述对称性处理子模块 314，利用视野空间F401和双平板PET探测器101结构对称性，简化 系统响应矩阵的计算，具体为：

步骤6.1，利用扇区间对称性，确定计算所需体素数和响应线：

只需确定一个扇区的体素，选取一个扇区中体素S1作为计算体 素，体素数为M；针对确定体素，利用数据采集模块1旋转过程的对 称性，只需确定在起始位置T201时的晶体构成的响应线，响应线条 数为N：

步骤6.1.1，确定与第一个扇区体素V502相对应的第i个扇区 的体素V501的索引号index₁、index_i；由体素索引方式(层c，扇区i， 环h，体素v)，对于第一个扇区的体素标记为(c,1,h,v)，则第i个扇区 上的体素标记为(c₁,i,h₁,v₁)，若第i个扇区的体素与第一个扇区S1相 对应，则需满足

c＝c₁

h＝h₁

v＝v₁

index_i＝v₁+(h₁-1)²+H²×(i-1)+c×C_num

index₁＝v+(h-1)²+c×C_num；

其中，C_num表示每层的体素数；

步骤6.1.2，确定体素V501与体素V502对应的响应线编码

设定穿过体素V502的响应线l为(晶体i'，晶体j'),计算晶体i' 和晶体j'在探测板晶体环上的所在边晶体位置；

c_i'＝i'/CR_num

c_j'＝(j'+n×m)/CR_num；

其中，/表示取余操作，CR_num表示每个扇区的体素总数；

设定穿过体素V501的响应线l_i为(晶体ii,晶体jj),计算晶体ii与 晶体jj在探测板晶体环406上所在边的晶体位置：

c_ii＝(ii+(i-1)×n)/CR_num

c_jj＝(jj+(i-1)×n)/CR_num；

其中，/表示取余操作，CR_num表示每个扇区的体素总数；

c_i'＝c_ii

如果，c_j'＝c_jj,则系统响应矩阵为：

Q(l_i,index_i)＝Q(l,index₁)；

步骤6.2，利用扇区内对称性，进一步减少计算响应线数

对于平行响应线满足穿过体素的上截面的位置相同时，只需计算 穿过扇区最长的响应线L502与相交体素的长度；对于与之平行的响 应线L501，只需计算穿过扇区的体素数目即可；

以二维情况进行说明，如图6所示，L502穿过的像素为 C″＝(6₁₁,6₁₀,6₉,6₈,5₅,5₄,5₃)，其中i_j表示第i个像素环的第j个像素，n 为7，计算穿过此n个像素长度值存储在数据D中，假定 D＝(a,b,c,d,e,f,g)，对于响应线L501穿过的像素为 C＝(6₇,6₆,6₅,6₄,4₁)，n为5，则有以下参数化表示：

h(C(i),L501)＝h(C″(i),L502)；

步骤7，在所述体素重建子模块321内，利用计算所得系统响应 矩阵，通过迭代算法MLEM重建算法重建出多面体体素的体素值，对 划分多面体体素进行重建；

步骤8，由所述转化子模块322完成由多面体体素到正方体体素 的转化，获取核素分布。如图7所示，转化是设定多面体体素与正方 体体素的厚度相同，通过计算两者截面的重叠面积，如M601、M602， 获取转化矩阵元素。具体步骤如下：

步骤8.1，令多面体体素的厚度和正方体体素的边界长度相等， 记正方体截面面积为SF，正方体体素数为FC_num，多面体体素数为 VC_num；

步骤8.2，计算体素的重叠体积RP，确定转化矩阵T，T的行列 数分别为VC_num,FC_num；T(i,j)表示第i个三角体素与第j个正方体体 素的截面的重叠面积和正方体体素截面积的比值，即

步骤8.3，利用转化矩阵T求解出每层上正方体体素的体素值

FC_i＝T(j,i)×VC_j；

步骤8.4，对每层重复步骤步骤8.3，获取最终的正方体体素划 分体素值。

对于本领域的技术人员来说，可以根据以上的技术方案和构思， 作出各种相应的改变和变形，而所有的这些改变和变形都应该包括在 本发明权利要求的保护范围之内。