CN109498048A - 一种加速正电子图像重建的系统矩阵生成与处理方法 - Google Patents

Abstract

本发明公开了一种加速正电子图像重建的系统矩阵生成与处理方法，属于射线检测成像领域。该方法包括如下步骤：（1）替代原先MLEM算法中系统矩阵的建立方法提出一种新的系统矩阵的建立方法；（2）在MLEM算法的基础上提出一种通过处理系统矩阵加速正电子图像重建过程的方法。本发明提高了正电子图像重建过程的整体效率。

Description

一种加速正电子图像重建的系统矩阵生成与处理方法

技术领域

本发明涉及一种加速正电子图像重建的系统矩阵生成与处理方法，属于射线检测成像领域。

背景技术

PET(Posintron Emission Tomography，正电子发射型计算机断层显像)，通过对测试对象注入放射性示踪剂，通过对示踪剂衰变发射出的正电子产生的湮没反应的检测，可对肿瘤、癌症、生理代谢、神经受体等进行定位、诊断与评估。PET检测主要分为四个步骤：

一、符合线检测：通过对正电子湮没反应产生的成对γ进行检测，获得大量符合响应线(Line of Response,LOR)数据。

二、图像重组：将大量的响应线数据转化为包含投影信息的正弦图(Sinogram)数据。

三、图像重建：将正弦图数据还原为平面切片图。

四、3D图像重建：将平面切片图重建为3D图。

图像重建作为PET检测中最重要的一环，决定了是否能将检测到的抽象的正弦图数据转化为直观的切片图数据。其主要的算法分为解析法与迭代法，解析法以中心切片定理为基础，通过反投影实现图像重建，滤波反投影法(Filter Back-Projection，FBP)为其代表算法，具有操作简单，重建速度快的优点，但是解析法不可避免的会受到投影数据噪声的影响，导致重建出的图像带有较大噪声，降低了图像质量，无法得到令人满意的重建图像；迭代法分为代数迭代和统计迭代，代数迭代法包括代数迭代(ART)算法，和在ART算法上提高了收敛速度的联合代数迭代(SART)算法，统计迭代法中最为经典的是基于泊松分布模型的最大似然期望最大(Maximum Likelihood Expectation Maximization，MLEM)算法，在其基础上通过不断改进还提出来了最大后验(Maximum A Posterior，MAP)算法、基于有序子集期望最大(Ordered Subsets Expectation Maximization)算法。

迭代算法中最为重要的参数称为系统矩阵，系统矩阵充当着正弦图数据通向切片图数据的桥梁的角色，反映了待还原的图像对于响应线的贡献，即正电子在该图像上发生湮没的概率分布，目前最为常用的系统矩阵为解析法系统矩阵，该系统矩阵通过响应线穿过像素的量化长度来表示正电子在此像素点上湮没产生该响应线的概率。对于距像素较远，没有穿过该像素的响应线，正电子在这个像素上湮没产生该条响应线的概率为0，称为零概率数据，这也就导致了在庞大的解析法系统矩阵中记录的大多数数据皆为零投影数据，不仅占用了大量的存储空间，且在进行图像重建时，会对零概率数据进行迭代运算处理，使得图像重建过程缓慢、耗时长、占用资源大，无法满足一些对实时性要求高的检测场合。

发明内容

`针对MLEM算法中系统矩阵数据量庞大且多数数据为零概率数据的情况，本发明提出了一种加速正电子图像重建的系统矩阵生成与处理方法,通过舍去系统矩阵中零概率数据，从而对重建过程进行加速的方法，减少了重建过程中资源的占用，加快了重建速度，提高了正电子图像重建过程的整体效率。

本发明为解决其技术问题采用如下技术方案：

一种加速正电子图像重建的系统矩阵生成与处理方法，包括如下步骤：

(1)替代原先MLEM算法中系统矩阵的建立方法提出一种新的系统矩阵的建立方法；

(2)在MLEM算法的基础上提出一种通过处理系统矩阵加速正电子图像重建过程的方法。

步骤(1)中所述系统矩阵的建立方法步骤如下：

步骤1，所述重建图像尺寸为N×N，探测器总数为X个，以相同的角度差dφ将所有探测器组成环形结构，构成一个圆形探测视野；对重建图像进行量化，定义重建图像的一个像素为边长为1的正方形，以N²个方形像素组成的正方形区域中心为原点建立直角坐标系，同时探测器的圆形探测视野内切于该正方形区域，计算各个像素方格中心在直角坐标上的的坐标(x_i,y_j),i,j＝1,2,3,…,N；

步骤2，将X个探测器以顺时针方向依次编号，把连线过圆形视野圆心的两个探测器划分为一对探测器，一共划分为M对，其中，M＝X/2，即1号探测器与M+1号探测器组成第1对探测器对，2号探测器与M+2号探测器组成第2对探测器对，以此类推；在计算响应线的直线方程时，通过响应线中垂线与x轴的夹角θ_m先确定响应线斜率k_m，根据探测器排布规则得知，每对探测器所连成的直线恰好为对应角度的响应线的中垂线，将以同一对探测器对所连成的直线为中垂线的所有响应线划分为一个响应线组，因此第m对探测器对应的第m组响应线组的每一条响应线的中垂线与x轴的夹角则第m组响应线组的斜率每组响应线组包括了M响应线，根据这M条响应线距离第m号探测器的远近依次对其进行编号，距离第m号探测器最近的响应线为第m组第一条响应线，距离第m号探测器最远的为第m组第条响应线；对于响应线M组的第n条响应线，计算接收到该条响应线的两个探测器与原点连线的夹角φ_n＝(r+1)dφ，其中，dφ为相邻探测之间的角度差，r为接收到该条响应线的两个探测器间相隔的探测器数，则该条响应线方程的斜率其中，N为环形探测器组成的圆形视野的直径，则第m组第n条响应线y_m-n＝k_mx+c_n，其中x为横坐标；根据上述规则，依次计算每组每条响应线方程并将其转化为Ax+By+C＝0的形式；

步骤3，计算各像素中心到每条投影线的距离其中，(x_i,y_j)为各像素中心坐标,A为步骤2中所得到的响应线方程表达式转化后对应于横坐标x_i的系数，B为步骤2中所得到的响应线方程表达式转化后对应于纵坐标y_j的系数，C为步骤2中所得到的响应线方程表达式转化后对应的常量系数；所述像素中心到投影线距离计算顺序为：取第1组第1条响应线，以先行后列的顺序开始计算各个像素中心到该投影线的距离s_1,1,i,j，得到一组1×N²的一维数组，取第1组第2条响应线，以先行后列的顺序计算各个像素中心到该投影线距离s_1,2,i,j，得到第二组1×N²的一维数组，与第一组数组组合成为2×N²的二维数组，以此类推，当遍历第1组响应线组所有响应线后，得到一组M×N²的二维数组，再取第2组第1条响应线，重复以上步骤，当遍历所有响应线后，得到一组M²×N²的二维数组；

步骤4，将步骤3所求二维数组表示为S，其中元素以s_ij重新表示，当像素中心距离响应线线较远时，认定在像素上发生正电子湮没产生这条响应线的概率为0，根据量化规则，将距离大于1的数据视作零概率数据，赋值为0，其余距离小于1的数据视作概率值存储在系统矩阵中，概率值以a_ij表示，则公式表示为：

其中i＝1,2,3,…,M²，j＝1,2,3,…,N²。

步骤(2)的具体过程如下：

步骤A，将正弦图数据从二维矩阵沿列方向拉伸成为一维矩阵，M对探测器得到M×M大小的正弦图数据，通过拉伸将其转化为M²×1的一维矩阵，表示为Y_i，需解算的切片图大小设为N×N，将其沿列方向拉伸为N²×1的一维矩阵，表示为λ_j，得知，二者抽象关系为Y＝a*λ，其中，a为系统矩阵，Y为正弦图数据矩阵，λ为切片图数据矩阵；

步骤B，通过MLEM算法计算λ_j，利用不定点迭代可得MLEM算法迭代公式为：

其中，k为已迭代次数，λ_j ^k+1为迭代k+1次后的一维切片图中的第j个像素值，λ_j ^k为迭代k次后的一维切片图中的第j个像素值，a_in为以i为行信息索引的系统矩阵第ni行的所有有效值，λ_n ^k为迭代k次后的一维切片图中的第n个像素值,i＝1,2,3,…,M²,n,j＝1,2,3,…,N²；

步骤C，根据步骤B中所述的MLEM算法迭代公式进行正电子图像重建，公式拆解为互不影响的两部分，所述互不影响的两部分分为与对于当j为确定值j₁时，代表对λ_j这个一维矩阵中的第j₁个元素进行求解，以j₁为索引，查找所有列下标信息为j₁的数据a_ij1求和，得出部分的结果；

λ_j ^k代表上一代迭代的结果，若当前迭代为初次迭代，即k＝0时，λ_j ^k默认全为1；

对于部分，当j为确定值j₁时，以j₁为索引，查找所有列下标信息为j₁的数据与其相应的行下标信息i，故该部分最终运算形式为:

其中，为一维正弦图数据中的第i₁个数据，为一维正弦图数据中的第i₂个数据，以此类推，为一维正弦图数据中的第i_s个数据；为以(i₁,j₁)为索引出的系统矩阵数据，为以(i₂,j₁)为索引出的系统矩阵数据，为以(i_s,j₁)为索引出的系统矩阵数据，为以i₁为行信息索引出的所有系统矩阵数据，为以i₂为行信息索引出的所有系统矩阵数据，为以i_s为行信息索引出的所有系统矩阵数据；

对于∑_na_inλ_n ^k部分，当i为确定值i₁时，以i₁为索引，查找所有行下标信息为i₁的数据与λ_n相乘后求和，得出∑_na_inλ_n ^k部分的结果；

λ_n ^k代表上一代迭代的结果，若当前迭代为初次迭代，即k＝0时，λ_n ^k默认全为1。

本发明的有益效果如下：

相比于现有的解析法系统矩阵，替代通过响应线穿过像素的长度来表示正电子在该像素上发生湮没的概率值的方法，采用像素中心到响应线的距离表示湮没的概率值来建立系统矩阵，简化了系统矩阵的建立过程，并对系统矩阵中的零概率数据进行舍去，大大地减少了系统矩阵的存储空间，并利用索引表提取系统矩阵有效数据，减少了迭代过程中对系统矩阵数据处理的循环次数，加速了正电子图像的重建过程。

附图说明

图1为响应线穿过像素点示意图。

图2(a)为探测器环形结构示意图；图2(b)为各像素中心坐标示意图。

图3为响应线示意图。

具体实施方案

以下将结合附图，对本发明的技术方案进行详细说明。

图像重建是通过正弦图来重建断层图像形成切片图的过程，系统矩阵作为正弦图重建为切片图的桥梁，记录了各响应线穿过各个像素的量化长度，用以反映正电子在切片图中的每个像素上发生湮没产生不同响应线的概率分布。为了方便说明，假设一个微型正电子探测装置，如图1所示，在该微型正电子探测环中的视野中检测一幅4×4的切片图，该探测装置其中一对探测器AA'接收到符合数据，即AA'这条响应线的数据，由正电子湮没原理可知，只有响应线穿过的像素才有可能产生该条响应线在不考虑散射的情况下，远离该响应线的像素释放的光子能量也不可能属于这条响应线，图1中响应线AA'穿过像素4，像素7，像素10，像素13，像素14，认定正电子有可能在以上像素发生湮没产生响应线AA'，其余像素未被响应线AA'穿过，认定正电子在其余像素上湮没产生响应线线AA'的概率为0。且穿过的像素长度越长，发生湮没的概率就越大，图1中响应线线穿过像素10的长度最长，因此正电子在像素10上湮没产生响应线AA'的概率最大，所述概率值即系统矩阵内元素的值。由于计算响应线穿过像素长度较为复杂，为了方便计算，本发明将像素中心到响应线线的距离进行量化替代上述穿过像素长度的量化来表示概率值，二者成比例关系，都可反映正电子在像素上湮没产生响应线的概率分布情况。下面举例说明系统矩阵的建立方法：

步骤1，解析法系统矩阵的大小由探测器对数与重建图像尺寸决定，如图2(a)所示，6对探测器以相同的角度差将组成环形结构，构成一个圆形探测视野用以重建一幅6×6的切片图，I为探测器组成的圆形探测视野，如图2(b)所示，对重建图像进行量化，定义重建图像的一个像素为边长为1的正方形，以6²个方形像素组成的正方形区域中心为原点建立直角坐标系，同时探测器的圆形探测视野内切与该正方形区域，计算各个像素方格中心在直角坐标上的的坐标(x₁,y₁),(x₁,y₂),…,(x₂,y₁),…,(x₆,y₆)。

步骤2，如图3所示，将探测器以顺时针方向依次编号，将在同一条直径上的探测器划分为一对，即1号探测器与7号探测器组成第1对探测器对，2号探测器与8号探测器组成第二对探测器对，以此类推，一共划分为6对连线过视野圆心的探测器对。以第3对探测器对，即3号探测器与9号探测器组成的探测器对为例，计算与第3对探测器对相对应的第3组响应线组的所有响应线的直线方程，第3对探测器对所在直线为相对应的第3组响应线组的中垂线，则可知第3组响应线组的斜率根据第3组响应线组中包含的所有响应线距离3号探测器的远近对所有的响应线进行编号，距离最近的第1条响应线，最远的为第6条响应线，可知，第3组第1条响应线y_3-1即图3中的I是由2号探测器与4号探测器接收符合数据产生的，2号探测器与4号探测器之间的相隔1个探测器，故这两个探测器与原点连线的夹角则该条响应线方程的截距故第3组第1条响应线方程为对于第3组第2条响应线y_3-2对应图3中的II，该条响应线是由1号探测器与5号探测器接收符合数据产生的，这组探测器之间相隔3个探测器，故1号探测器与5号探测器之间的角度差接收到的响应线方程的截距故第3组第2条响应线方程为根据上述规则，依次计算每组每条响应线并将其转化为Ax+By+C＝0的形式。

步骤3，计算各像素中心到每条投影线的距离其中：(x_i,y_j)为各像素中心坐标，A为步骤2中所得到的响应线方程表达式转化后对应于横坐标x_i的系数，B为步骤2中所得到的响应线方程表达式转化后对应于纵坐标y_j的系数，C为步骤2中所得到的响应线方程表达式转化后对应的常量系数。所述像素中心到投影线距离计算顺序为：取第1组第1条响应线，以先行后列的顺序开始计算各个像素中心到该投影线的距离s_1,1,1,1,s_1,1,1,2,…,s_1,1,6,6，，得到一组1×6²的一维数组，取第1组第2条响应线，以先行后列的顺序计算各个像素中心到该投影线距离s_1,2,1,1,s_1,2,1,2,…,s_1,2,6,6，，得到第二组1×6²的一维数组，与第一组数组组合成为2×6²的二维数组，以此类推，当遍历第1组响应线组所有响应线后，得到一组6×6²的二维数组，再取第2组第1条响应线，重复以上步骤，当遍历所有响应线后，得到一组6²×6²的二维数组。

步骤4，将步骤3所求二维数组表示为S，其中元素以s_ij重新表示，当像素中心距离响应线线较远时，认定在像素上发生正电子湮没产生这条响应线的概率为0，根据量化规则，将距离大于1的数据视作零概率数据，赋值为0；其余距离小于1的数据视作概率值，存储在系统矩阵中以a_ij表示，则公式表示为：

其中：i＝1,2,3,…,36，j＝1,2,3,…,36。

由上述方法得到的大小为4²×4²的系统矩阵，其中零概率数据量占总数据量的62.5％，且随着系统矩阵的不断增大，零概率数据量占总数据量的百分比是在不断增大的，表1表示了不同规模的系统矩阵其零数据量占总数据量的百分比：

表1系统矩阵零数据占比

系统矩阵规模规模	4<sup>2</sup>×4<sup>2</sup>	40<sup>2</sup>×64<sup>2</sup>	92<sup>2</sup>×128<sup>2</sup>
				零数据占比(％)	62.5	97.4	98.7

下面以利用6²×6²系统矩阵将一幅6×6的正弦图转化成一幅6×6切片图为例，说明利用去“零”系统矩阵进行MLEM算法的具体步骤：

步骤1，通过先列后行的方式将6×6的正弦图拉伸为6²×1的一维矩阵，表示为Y_i(i＝1,2,3,…,36)，设置一组6²×1的一维矩阵，表示为λ_j(j＝1,2,3,…,36)，用以存储图像重组后的切片图数据。

步骤2，遍历系统矩阵中每一个元素a_ij,判断a_ij数值，若等于0视为零概率数据，则将其舍去不做记录，若不等于0视为有效概率数据，将其数值与在系统矩阵中的位置信息记录在表3中，以该系统矩阵第一行为例，依次对第一行元素进行判断，舍去零概率数据，记录下有效贡献数据数值与下标位置信息与表2中

表2 6²×6²系统矩阵有效数据及位置信息表

数值a<sub>ij</sub>	0.9019	0.9019	0.9019	0.9019	0.9019	0.9019
							行信息i	1	1	1	1	1	1
列信息j	1	7	13	19	25	31

步骤3，通过MLEM算法计算λ_i，利用不定点迭代可得MLEM算法迭代公式为：

其中：λ_j ^k+1为迭代k+1次后的一维切片图中的第j个像素值，λ_j ^k为迭代k次后的一维切片图中的第j个像素值，a_in为以i为行信息索引的系统矩阵第i行的所有有效值，λ_n ^k为迭代k次后的一维切片图中的第n个像素值。

以求解切片图数据λ₁为例，将j＝1带入公式，得到：对于以j＝1为索引，查找所有列下标信息为j＝1的数据a_i1，以a_i,j表示，i为行下标，j为列下标，本例子中根据j＝1索引出的数据为a_1,1，a_7,1，a_8,1，a_9,1，a_13,1，a_14,1，a_19,1，a_20，1，a_25，1，a_26，1，a_31，1，a_32，1将上述数据求和，得出部分的结果。其中，λ₁ ⁰代表初次迭代的前次迭代，取为1。

同样的，对于部分，以j＝1为索引，查找所有列下标信息为j＝1的数据a_ij1以及其相应的行下标信息i，，以a_i,j表示，i为行下标，j为列下标，本例子中根据j＝1索引出的数据为a_1,1，a_7,1，a_8,1，a_9，1，a_13，1，a_14,1，a_19，1，a_20，1，a_25，1，a_26，1，a_31，1，a_32，1，故该部分最终运算形式为:

其中：Y₁表示上述6²×1一维正弦图数据中的第1个数据，Y₇表示上述6²×1一维正弦图数据中的第7个数据，以此类推，Y₃₂表示上述6²×1一维正弦图数据中的第32个数据，a_1,n为以1为行信息索引的系统矩阵第1行的所有有效值，a_3,n为以3为行信息索引的系统矩阵第3行的所有有效值，以此类推，a_32,n为以32为行信息索引的系统矩阵第32行的所有有效值。

对于∑_na_inλ_n ^k部分，当ⁱ为确定值ⁱ ₁时，以i₁为索引，查找所有行下标信息为i₁的数据与λ_n相乘后求和，得出∑_na_inλ_n ^k部分的结果。

以求解∑_na_1,nλ_n ⁰为例，以i＝1为索引，查找所有行下标信息为i＝1的数据，以a_i,j表示，i为行下标，j为列下标，本例子中根据i＝1索引出的数据为a_1,2，a_1,8，a_1,14，a_1,20，a_1,26，a_1,32，故该求和的最终运算形式为：a_1，2λ₂ ⁰+a_1，8λ₈ ⁰+…+a_1，32λ₃₂ ⁰

其中，λ_n ⁰代表初次迭代的前次迭代，均取为1。

表3为在主频为3.40Ghz计算机下采用本发明处理后的系统矩阵进行MELM重建图像的时间与原MLEM算法重建图像的时间的对比表，可知，采用本发明处理后的去“零”系统矩阵的MLEM算法重建图像时间相比于原MLEM算法有明显的减少。

表3图像重建时间对比表

Claims (5)

1.一种加速正电子图像重建的系统矩阵生成与处理方法，其特征在于，包括如下步骤：

(1)替代原先MLEM算法中系统矩阵的建立方法提出一种新的系统矩阵的建立方法；

(2)在MLEM算法的基础上提出一种通过处理系统矩阵加速正电子图像重建过程的方法。

2.根据权利要求1所述的一种加速正电子图像重建的系统矩阵生成与处理方法，其特征在于，步骤(1)中所述系统矩阵的建立方法步骤如下：

步骤1，所述重建图像尺寸为N×N，探测器总数为X个，以相同的角度差dφ将所有探测器组成环形结构，构成一个圆形探测视野；对重建图像进行量化，定义重建图像的一个像素为边长为1的正方形，以N²个方形像素组成的正方形区域中心为原点建立直角坐标系，同时探测器的圆形探测视野内切于该正方形区域，计算各个像素方格中心在直角坐标上的的坐标(x_i,y_j),i,j＝1,2,3,…,N；

步骤2，将X个探测器以顺时针方向依次编号，把连线过圆形视野圆心的两个探测器划分为一对探测器，一共划分为M对，其中，M＝X/2，即1号探测器与M+1号探测器组成第1对探测器对，2号探测器与M+2号探测器组成第2对探测器对，以此类推；在计算响应线的直线方程时，通过响应线中垂线与x轴的夹角θ_m先确定响应线斜率k_m，根据探测器排布规则得知，每对探测器所连成的直线恰好为对应角度的响应线的中垂线，将以同一对探测器对所连成的直线为中垂线的所有响应线划分为一个响应线组，因此第m对探测器对应的第m组响应线组的每一条响应线的中垂线与x轴的夹角则第m组响应线组的斜率每组响应线组包括了M响应线，根据这M条响应线距离第m号探测器的远近依次对其进行编号，距离第m号探测器最近的响应线为第m组第一条响应线，距离第m号探测器最远的为第m组第条响应线；对于响应线M组的第n条响应线，计算接收到该条响应线的两个探测器与原点连线的夹角φ_n＝(r+1)dφ，其中，dφ为相邻探测之间的角度差，r为接收到该条响应线的两个探测器间相隔的探测器数，则该条响应线方程的斜率其中，N为环形探测器组成的圆形视野的直径，则第m组第n条响应线y_m-n＝k_mx+c_n，其中x为横坐标；根据上述规则，依次计算每组每条响应线方程并将其转化为Ax+By+C＝0的形式；

步骤3，计算各像素中心到每条投影线的距离其中，(x_i,y_j)为各像素中心坐标,A为步骤2中所得到的响应线方程表达式转化后对应于横坐标x_i的系数，B为步骤2中所得到的响应线方程表达式转化后对应于纵坐标y_j的系数，C为步骤2中所得到的响应线方程表达式转化后对应的常量系数；所述像素中心到投影线距离计算顺序为：取第1组第1条响应线，以先行后列的顺序开始计算各个像素中心到该投影线的距离s_1,1,i,j，得到一组1×N²的一维数组，取第1组第2条响应线，以先行后列的顺序计算各个像素中心到该投影线距离s_1,2,i,j，得到第二组1×N²的一维数组，与第一组数组组合成为2×N²的二维数组，以此类推，当遍历第1组响应线组所有响应线后，得到一组M×N²的二维数组，再取第2组第1条响应线，重复以上步骤，当遍历所有响应线后，得到一组M²×N²的二维数组；

步骤4，将步骤3所求二维数组表示为S，其中元素以s_ij重新表示，当像素中心距离响应线线较远时，认定在像素上发生正电子湮没产生这条响应线的概率为0，根据量化规则，将距离大于1的数据视作零概率数据，赋值为0，其余距离小于1的数据视作概率值存储在系统矩阵中，概率值以a_ij表示，则公式表示为：

其中i＝1,2,3,…,M²，j＝1,2,3,…,N²。

3.根据权利要求2所述的一种加速正电子图像重建的系统矩阵生成与处理方法，其特征在于，步骤1中所述重建图像尺寸为4×4。

4.根据权利要求2所述的一种加速正电子图像重建的系统矩阵生成与处理方法，其特征在于，步骤2中所述探测器的对数为6对。

5.根据权利要求1所述的一种加速正电子图像重建的系统矩阵生成与处理方法，其特征在于，步骤(2)的具体过程如下：

步骤A，将正弦图数据从二维矩阵沿列方向拉伸成为一维矩阵，M对探测器得到M×M大小的正弦图数据，通过拉伸将其转化为M²×1的一维矩阵，表示为Y_i，需解算的切片图大小设为N×N，将其沿列方向拉伸为N²×1的一维矩阵，表示为λ_j，得知，二者抽象关系为Y＝a*λ，其中，a为系统矩阵，Y为正弦图数据矩阵，λ为切片图数据矩阵；

步骤B，通过MLEM算法计算λ_j，利用不定点迭代可得MLEM算法迭代公式为：

其中，k为已迭代次数，λ_j ^k+1为迭代k+1次后的一维切片图中的第j个像素值，λ_j ^k为迭代k次后的一维切片图中的第j个像素值，a_in为以i为行信息索引的系统矩阵第ni行的所有有效值，λ_n ^k为迭代k次后的一维切片图中的第n个像素值,i＝1,2,3,…,M²,n,j＝1,2,3,…,N²；

步骤C，根据步骤B中所述的MLEM算法迭代公式进行正电子图像重建，公式拆解为互不影响的两部分，所述互不影响的两部分分为与对于当j为确定值j₁时，代表对λ_j这个一维矩阵中的第j₁个元素进行求解，以j₁为索引，查找所有列下标信息为j₁的数据求和，得出部分的结果；

λ_j ^k代表上一代迭代的结果，若当前迭代为初次迭代，即k＝0时，λ_j ^k默认全为1；

对于部分，当j为确定值j₁时，以j₁为索引，查找所有列下标信息为j₁的数据与其相应的行下标信息i，故该部分最终运算形式为:

其中，为一维正弦图数据中的第i₁个数据，为一维正弦图数据中的第i₂个数据，以此类推，为一维正弦图数据中的第i_s个数据；为以(i₁,j₁)为索引出的系统矩阵数据，为以(i₂,j₁)为索引出的系统矩阵数据，为以(i_s,j₁)为索引出的系统矩阵数据，为以i₁为行信息索引出的所有系统矩阵数据，为以i₂为行信息索引出的所有系统矩阵数据，为以i_s为行信息索引出的所有系统矩阵数据；

对于∑_na_inλ_n ^k部分，当i为确定值i₁时，以i₁为索引，查找所有行下标信息为i₁的数据与λ_n相乘后求和，得出∑_na_inλ_n ^k部分的结果；

λ_n ^k代表上一代迭代的结果，若当前迭代为初次迭代，即k＝0时，λ_n ^k默认全为1。