基于Mojette变换的CT重建方法

摘要

本发明提供了基于Mojette变换的CT重建方法，属于计算机层析成像技术领域。该方法建立了投影像元与当前待重建像素间的对应关系，根据投影像元的不同重建次序形成当前待重建像素不同的重建路径，分析了不同重建路径对应不同迭代次数和噪声累积量，找到了最少迭代次数的迭代重建算法，解决了迭代重建过程中噪声在重建路径上繁衍到一个很大数值的问题，使得基于最优重建路径的重建算法在小噪声条件下也能重构出效果最好的重建图像。本发明在有限投影集的情况下，提出了最少迭代次数的迭代重建算法，减少了投影更新次数，避免了重建路径上极小的噪声也会在求解的迭代过程中不断被放大，保证少量的噪声条件下也能重构出效果最好的重建图像。

权利要求

1.一种基于Mojette变换的CT重建方法，其特征在于，步骤如下：

S1.设待重建目标图像大小为M×N，实际采样的探测器分辨率为DetRowNum，则基于Mojette变换，依据公式(5)将求解得到所有投影方向Ptotal，其中B(M，N，pi，qi)是投影方向(pi,qi)下的投影像元个数，随着投影方向(pi,qi)变化而变化，且所有满足公式(5)的投影方向整合在一起表达为Ptotal＝{(p1,q1)(p2,q2)...(pi,qi)...(pI,qI),i＝1,2,...,I}，I是总投影个数；

S2.在没有噪声情况下，从步骤S1求解得到的所有投影方向Ptotal中任意选择一些投影构成投影子集，该投影子集满足Katz引理，则能够精确重建断层图像；在噪声情况下，满足Katz引理投影子集不一定能完成重建图像，因为投影数据中的噪声在迭代步骤中不断被放大，使得重建断层遭到严重破坏，所以需要增加更多可利用的投影冗余信息换取待重建图像的质量；因此，在Katz引理中，通过引入N2增加Katz引理的上限，引入N1降低Katz引理的下限，从而获得扩展Katz引理，即公式(6)；当扩展Katz引理成立时，能够获得更多投影参与重建目标图像，进而减少噪声繁衍次数，改进重建图像质量，其中，N1，N2是整数，I*是参与目标图像重建的总投影数，即为最优投影角度集Poptimal，I*≤I；

S3.针对最优投影角度集Poptimal，确定最优重建路径的第一个投影方向：首先遍历所有投影方向Ptotal，计算出每一个投影方向下只穿过一条投影射线的射线个数，将射线个数最多的投影方向作为最优重建路径的第一投影方向；即通过求解目标函数(7)，求得初次迭代中可直接重建的像素点数量最多的初始投影方向(pi,qi)，在初始投影方向(pi,qi)下，第1次迭代步骤中重建出来的像素点总数N1最大；

N1＝{(p1，q1)＝max(|pi|·|qi|)s.t.|pi|·N+|qi|·M+1≤DetRowNum}    (7)S4.确定反投影方案中映射矩阵Jrad(m,n)：映射矩阵建立了初始投影方向(pi,qi)上的投影像元与对应像素之间的联系，找出初始迭代步骤中投影方向下可直接重建的像素点；

类似传统的CBI重建算法，计算初始迭代步骤中投影方向(pi,qi)上每个投影像元对应的一次迭代能够重建的像素个数，映射矩阵将初始迭代步骤中能够重建最多图像像素的投影直接反投影到对应的离散图像像素(m,n)中，其中映射矩阵公式记为Jrad(m,n)＝bin；

S5.针对确定最优重建路径的第一个投影方向后，进一步确定最优重建路径的第二个投影方向：遍历所有投影方向Ptotal，利用步骤S4中得到的离散图像像素(m,n)寻找当前迭代过程中待重建像素，计算出当前迭代步骤中每一个投影方向下只穿过一条投影射线的射线个数，将射线个数最多的投影方向作为最优重建路径的第二个投影方向(p2,q2)，其中待重建像素数目N2为；

其中 表示当前迭代步骤中投影方向(pi,qi)能重建出来的像素数目；

确定第二个投影方向(p2,q2)后，开始第二次反投影重建，将当前迭代步骤中的投影像元直接反投影到对应像素索引坐标中，具体步骤同S4；

S6.返回步骤S5，重复迭代步骤，直到所有像素被重建；对于尺寸大小为M×N待重建目标图像，总像素个数为M×N；在第j次迭代步骤中，被重建得到的像素个数Nj表示为公式(9)；

其中 表示当前迭代过程中(pi，qi)方向的投影能重建出来的像素数目，J为总迭代次数，Etotal为目标图像重建完成后的总迭代次数；当前迭代步骤中，变量值越大，总迭代次数Etotal越小；

在每一次迭代过程中，首先需要遍历所有投影方向Ptotal，找到投影方向(pi,qi)上一次就能重建出来的像素个数，将投影射线上的实际投影直接反投影到对应像素上，重复迭代过程，直到所有像素被重建。

说明书

基于Mojette变换的CT重建方法

技术领域

[0001] 本发明属于计算机层析成像(Computed Tomography，CT)技术领域，它从属于图像处理领域，涉及一种在图像投影变换域稀疏采样及高效复原的方法。

背景技术

[0002] 计算机断层成像技术通过X射线穿过物体产生的衰减信息，能够在不接触、不破坏物体结构的前提下，对物体的内部结构进行成像。因此CT技术在医疗诊断、无损检测等领域有广泛的应用。

[0003] 计算机断层重建技术的投影过程可用Radon变换表示，并利用Radon反变换来重建断层。Radon变换是对物体连续采样，得到无限多角度下的投影数据。而在实际应用中，由于X射线对人体存在较大的辐射和病人运动状态的约束，只能采集到有限的投影数据。目前重建算法中，利用有限角的投影数据重建二维图像，已有很多研究工作，但是为了保证CT图像的质量，仍需采集投影的空间频率满足奈奎斯特准则，才能获得理想的断层图像。

[0004] 因此，如何在保证图像重建质量的前提下，设计新的高效采集模式，感知CT成像系统中的有用信息，减少投影采集的数量，从而减少处理时间和对生物样品的辐射剂量，并从有限采集的投影中高质量重建CT图像，已成为传统CT重建的热点和关键问题。

[0005] 因为基于有限角度的离散Radon逆变换本身是一种病态的逆问题，重建结果总是一个近似值，是不可能得到精确的重建图像。因此，在Radon变换的基础上，Guedon J P等提出了Mojette变换的概念，该变换是离散Radon变换的一种特殊形式，它不同于在离散域近似重建的Radon变换，在满足Katz引理的前提下，当投影数据中不含噪声时，基于Mojette变换的重建算法是一种精确的强健的重建算法，能够重建出与目标图像完全相同的精确的重建结果；当投影数据中包含有噪声时，Mojette逆变换算法对噪声太敏感，因为基于Mojette逆变换的精确重建算法只包含减法运算，逐点重建每个像素，像素的重建顺序是从四个角到图像中间，噪声从图像边沿像素繁衍到图像中间像素，如果完成图像重建的迭代步骤少，累积到图像中心的像素噪声也小；如果完成图像重建的迭代步骤多，累积到图像中心的像素噪声也大。总之，在循环更新投影过程中，噪声在投影数据中依次迭代，最后更新的投影数据受到噪声污染程度最大。所以，在重建过程中，投影迭代次数是非常重要的，迭代次数越少，噪声繁衍次数越少，重建图像的质量越高。本发明提出了最少迭代次数的重建算法，也称为最优重建路径重建算法。

[0006] 然而，为了更好理解本文阐述的重建算法，接下来，需要简述Mojette变换以及它的重要性质。Mojette变换中用一对互质的整数(pi,qi)来表达投影方向，一般有pi∈Z,qi∈Z+，pi代表了图像水平方向上的整数位移，qi代表了图像垂直方向上的整数位移，投影方向(pi,qi)对应的投影角度为θ＝tan-1qi/pi,I是选取的投影角度总个数，如图1所示。那么，该投影方向(pi,qi)上的所有平行投影射线表达为b＝qi·m-pi·n，Mojette变换公式为：

[0007]

[0008] 其中f(m,n)代表待重建的图像切片上索引坐标为(m,n)的像素点的灰度值，bi值代表(m,n)这一像素点经过Mojette变换打到第bi个探测器像元上，

[0009] Mojette变换是一种特殊的Radon变换，Mojette变换与Radon变换相比最大的不同在于，每个投影方向(pi,qi)下所需的探测器像元个数B(M,N,pi,qi)不同，这也意味着在覆盖相同直径范围内的切片时，不同投影角度下的投影射线之间的间距hi不同，即：

[0010]

[0011] 其中，M代表待重建图像在水平方向上的尺寸，而N代表了待重建图像在垂直方向上的尺寸。Mojette变换的重要性质如下所述：

[0012] 1.在Mojette投影域精确重建的条件。

[0013] 不同于在离散域近似重建的Radon变换，Mojette变换在满足Katz引理的前提下，可以在离散域内完成精确重建。针对于Mojette变换可精确重建的条件，1978年，Katz给出一个约束投影方向个数的公式，即著名的Katz引理，该定理指出对于由一些互质整数构成的投影方向(pi,qi)，在重建图像的大小为M×N情况下，将满足公式(3)：

[0014]

[0015] 即可通过I个投影完成精确重建。

[0016] 式(3)的重建条件可进一步化简为：

[0017]

[0018] 下面本节通过一个具体的实例来说明一下Katz引理。例如，对于一个3×3的待重建图像，即M＝3，N＝3，根据Katz引理，当选择的投影方向为(1,1)，(1,0)，(-1,1)时，有满足精确重建条件。若将这个Mojette投影过程看成是建立线性方程组的过程，则Mojette反重建被看作是求解线性方程组的过程。即可通过验证方程组的可解性来验证Katz引理的有效性。

[0019] 将3×3的待重建图像列矢量化为9×1的一维向量{x1,x2,x3,x4,x5,x6,x7,x8,x9}，根据图3(a)中所示的投影关系，列出在投影方向(1,1)下的投影线性方程组，如图4(a)所示；根据图3(b)中所示的投影关系，列出在投影方向(1,0)下的线性方程组，如图4中(b)所示；根据图3(c)中所示的投影关系，列出在投影矢量(-1,1)下的线性方程组，如图4中(c)所示。如果利用图4中的方程组联合解出离散像素{x1,x2,x3,x4,x5,x6,x7,x8,x9}的灰度值，则说明在满足Katz引理的条件下，Mojette正变换可逆。从图4(a),(b),(c)中看出，离散像素x1，x3，x7，x9的灰度值可以直接由投影值解出，且在同一个投影方程组中，不会有重复出现的变量。然后，利用投影方向(1,1)下解出的离散像素x1和x9去更新投影方向(1,0)和(-1,1)下的方程组。即，根据解出的离散像素x1，x3，x7和x9，化简图4，得到图5，且依据图5(a)求出x5的值，依据图5(b)求出离散像素x4，x6的值。再根据所求出的离散像素x1，x3，x5，x6，x7，x9的值又可求出离散像素x2和x8的值，如图6所示。

[0020] 2.分析不同重建路径上噪声累积的过程。

[0021] 基于Mojette逆变换的原理，Mojette逆变换的求解过程可以视为是一种串行的求解模式，即每次迭代求出一批像素点，再根据这些已知像素点的灰度值，来更新各个方向上投影，从而产生出投影值即为单个像素的灰度值，再重复这一步骤，直至所有的待求像素点都与探测器像元上更新完毕后的一个投影值一一对应，迭代结束。

[0022] 在这一求解逆变换的过程中，投影角度选择不同，每次迭代步骤中被重建出来的像素点就不同，且每一像素点所在的迭代次序也会随着前面迭代重建结果的不同而不同，称每一种情况下的重建次序为一条重建路径，如图7所示。

[0023] 在Mojette逆变换算法中，像素点的灰度值是被依次重建出来的，这与滤波反投影的批量重建，不断叠加的重建原理是截然不同的。投影方向、重建路径选择的不同，即像素点的重建次序的不同会直接造成算法的运算量和抗噪性能的不同。

[0024] 图8中给出了一个具体的实例来说明这一现象，在(0,1),(1,1),(-1,1)投影方向下重建3×3的待重建图像，图8中从探测器像元指向像素点的直线代表重建过程，从像素点指向探测器像元的直线代表投影更新过程，图8中的红色小圈里的数字，代表该位置处的像素点在哪一步迭代中被重建出来，如1指在第一次迭代重建步骤中。如图8中(b)所示，首先可以利用投影方向(1,1)下的投影直接求解出索引坐标为(1,1)和(3,3)的像素，然后通过图8(d)可知在(b)中解出的像素点值的基础上，可利用投影方向(-1,1)下的投影求出中心像素；即对索引坐标为(2,2)的中心像素点，若采用投影方向(1,1)和(-1,1)方向下穿过中心点的射线投影来重建，且共计迭代2次即可重建出该离散像素，重建结果为8-(e2+e3)，累加噪声为e2+e3。若选择(0,1)投影方向下穿过中心点的射线投影来重建，则迭代更新步骤如图8中(a)(b)(c)(e)(f)所示，共计迭代4次，重建结果为8-(e1+e2+e1+e2)，噪声误差比第一种重建路径增加了近一倍。

[0025] 通过这个例子可以看出，在Mojette重建过程中，从求解外围点到求解中心点的过程中，迭代次数逐渐增加，噪声也不断累加，每个点的迭代次数k反映了其噪声累加的次数，设在第k次迭代步骤中重建出来的像素点总数为Nk，重建所有像素点所需的迭代次数为K，若重建M×N的图像，则有N1+N2+…+Nk+…+NK＝M·N。若想让噪声累加次数Etotal最少，即使得总迭代次数K越小越好，同时初始迭代步骤中的像素点数N1越大越好，理想的极端情况中，假设所有的点一次迭代就可求出，则Emin＝M·N，此时为误差最小的理想情况；若每次迭代中都只能求解出一个点，Emax＝1+2+3+…+(M·N-1)+(M·N)＝(1+M·N)·(M·N)/2，此时为误差最大的极端不理想情况。因此，最大程度上减少投影误差对重建结果的影响，需要在最大程度上减少重建迭代次数，且尽量避免误差点重复参与运算，通过选择投影角度和合理安排重建点的次序，从而有效的降低噪声的累加次数，从而获得更良好的重建效果。

发明内容

[0026] 该方法提出了噪声最少繁衍次数的迭代重建算法，它基于Mojette变换和依据图像精确重建的Katz引理，然后给出了不同可能的重建路径和噪声在不同重建路径上不断累积过程，解决了迭代重建过程中噪声在重建路径上繁衍到一个很大数值的问题，使得在噪声条件下最优重建路径算法能重构出效果最好的重建图像。

[0027] 本发明的技术方案：

[0028] 基于Mojette变换的CT重建方法，步骤如下：

[0029] S1.设待重建目标图像大小为M×N，实际采样的探测器分辨率为DetRowNum，则基于Mojette变换，依据公式(5)将求解得到所有投影方向Ptotal，

[0030]

[0031] 其中B(M,N,pi,qi)是投影方向(pi,qi)下的投影像元个数，它随着投影方向(pi,qi)变化而变化，且所有满足公式(5)的投影方向整合在一起表达为Ptotal＝{(p1,q1)(p2,q2)...(pi,qi)...(pI,qI),i＝1,2,...,I}，I是总投影个数。

[0032] S2.在没有噪声情况下，从步骤S1求解得到的所有投影方向Ptotal中任意选择一些投影构成投影子集，该投影子集满足Katz引理，则能够精确重建断层图像。在噪声情况下，满足Katz引理投影子集不一定能完成重建图像，因为投影数据中的噪声在迭代步骤中不断被放大，使得重建断层遭到严重破坏，所以需要增加更多可利用的投影冗余信息换取待重建图像的质量；因此，在Katz引理中，通过引入N2增加Katz引理的上限，引入N1降低Katz引理的下限，从而获得扩展Katz引理，即公式(6)。当扩展Katz引理成立时，能够获得更多投影参*与重建目标图像，进而减少噪声繁衍次数，改进重建图像质量，其中，N1，N2是整数，I是参与目标图像重建的总投影数，即为最优投影角度集Poptimal，I*≤I；

[0033]

[0034] S3.针对最优投影角度集Poptimal，确定最优重建路径的第一个投影方向：首先遍历所有投影方向Ptotal，计算出每一个投影方向下只穿过一条投影射线的射线个数，将射线个数最多的投影方向作为最优重建路径的第一投影方向；即通过求解目标函数(7)，求得初次迭代中可直接重建的像素点数量最多的初始投影方向(pi,qi)，在初始投影方向(pi,qi)下，第1次迭代步骤中重建出来的像素点总数N1最大；

[0035] N1＝{(p1,q1)＝max(|pi|·|qi|)s.t.|pi|·N+|qi|·M+1≤DetRowNum}  (7)[0036] S4.确定反投影方案中映射矩阵Jrad(m,n)：映射矩阵建立了初始投影方向(pi,qi)上的投影像元与对应像素之间的联系，找出初始迭代步骤中投影方向下可直接重建的像素点；类似传统的CBI重建算法，计算初始迭代步骤中投影方向(pi,qi)上每个投影像元对应的一次迭代能够重建的像素个数，映射矩阵将初始迭代步骤中能够重建最多图像像素的投影直接反投影到对应的离散图像像素(m,n)中，其中映射矩阵公式记为Jrad(m,n)＝bin；

[0037] S5.针对确定最优重建路径的第一个投影方向后，进一步确定最优重建路径的第二个投影方向：遍历所有投影方向Ptotal，利用步骤S4中得到的离散图像像素(m,n)寻找当前迭代过程中待重建像素，计算出当前迭代步骤中每一个投影方向下只穿过一条投影射线的射线个数，将射线个数最多的投影方向作为最优重建路径的第二个投影方向(p2,q2)，其中待重建像素数目N2为；

[0038]

[0039] 其中 表示当前迭代步骤中投影方向(pi,qi)能重建出来的像素数目；

[0040] 确定第二个投影方向(p2,q2)后，开始第二次反投影重建，将当前迭代步骤中的投影像元直接反投影到对应像素索引坐标中，具体步骤同S4；

[0041] S6.返回步骤S5，重复迭代步骤，直到所有像素被重建；对于尺寸大小为M×N待重建目标图像，总像素个数为M×N；在第j次迭代步骤中，被重建得到的像素个数Nj表示为公式(9)；

[0042]

[0043] 其中 表示当前迭代过程中(pi,qi)方向的投影能重建出来的像素数目，J为总迭代次数，Etotal为目标图像重建完成后的总迭代次数；当前迭代步骤中，变量值越大，总迭代次数Etotal越小；

[0044] 在每一次迭代过程中，首先需要遍历所有投影方向Ptotal，找到投影方向(pi,qi)上一次就能重建出来的像素个数，将投影射线上的实际投影直接反投影到对应像素上，重复迭代过程，直到所有像素被重建。

[0045] 本发明的效果和益处：基于Mojette变换的精确图像重建算法，提出了最少迭代次数的迭代重建算法，避免了极小的噪声也会在求解的迭代过程中不断被放大到一个很大数值的问题，保证噪声条件下也能重构出效果最好的重建图像。

附图说明

[0046] 图1是投影方向(pi,qi)示意图。图中，Mojette变换用一对互质的整数(pi,qi)来表达投影方向，一般有pi∈Z,qi∈Z+，pi代表了图像水平方向上的整数位移，qi代表了图像垂直方向上的整数位移，投影方向(pi,qi)表达的投影角度为θ＝tan-1qi/pi，I是选取的投影角度总数。

[0047] 图2是对3×3的待重建图像进行Mojette正变换的示例图。图中左侧是3×3的待重建图像，右侧是经过Mojette正变换的三个投影方向(1,0),(1,1),(2,1)上的投影。

[0048] 图3是三个投影方向(1,0),(1,1),(-1,1)上的Mojette投影过程。图(a),(b),(c)分别是三个投影方向上的Mojette投影过程，每一个投影方向的投影射线构成了一组线性方程，对应图4。

[0049] 图4是与Mojette投影过程对应的线性方程组。离散像素x1，x3，x7，x9的灰度值，将通过投影射线直接反投影求解得到，更新线性方程组。

[0050] 图5是第一次更新后的投影线性方程组。离散像素x4，x5，x6的灰度值，将通过第一次更新的投影值求解得到，更新线性方程组。

[0051] 图6是第二次更新后的投影线性方程组。离散像素x2，x8的灰度值，将通过第二次更新的投影值求解得到，到此完成所有像素求解。

[0052] 图7是不同的投影角度对应不同的重建路径。图(a)是采用投影(-1,1)，(1,1)，(1,2)得到的重建路径，图(b)是采用投影(2,1)，(-2,1)，(1,2)，(-1,2)得到的重建路径。

[0053] 图8是噪声累加过程示意图。第一条重建路径需要经过图(a)，(b)，(d)，3次迭代步骤即能求解得到像素坐标为(2,2)的像素灰度值；第二条重建路径需要经过图(a)，(b)，(c)，(e)，(f)，5次迭代步骤才能求解得到像素坐标(2,2)的像素灰度值。

[0054] 图9是投影方向(5,4),(3,4)的图例。对于9×9的待重建图像，投影方向(5,4)只穿过两个离散像素，投影方向(3,4)穿过三个离散像素。

[0055] 图10是投影方向(5,4)上仅穿过一个像素的投影射线图例。图中，投影方向(5,4)上一些投影射线只穿过灰色区域的一个像素。

[0056] 图11是三个不同投影方向上一些投影射线只穿过一个像素的图例。图(a)，(b)，(c)，分别表示三个投影方向，图(a)投影方向上仅穿过一个像素的投影射线最多。

[0057] 图12是最优投影方向集P1的重建路径。最优路径上7次投影更新，即可完成图像重建。

[0058] 图13是非最优投影方向集P2的重建路径。非最优重建路径上需要109次投影更新，才能完成图像重建。

[0059] 图14是最优重建路径和非最优重建路径的重建结果。图14(a)是原图；图(b)(c)是二组非最优重建路径下的重建结果；图(d)是最优重建路径在噪声取值范围[-0.15％Mmax,0.15％Mmax]的重建结果；图(e)是最优重建路径在噪声取值范围[-0.1％Mmax,0.1％Mmax]的重建结果；从实验结果可知，最优重建路径下的重建结果远远好于非最优重建路径下的重建结果。

具体实施方式

[0060] 以下结合技术方案和附图详细叙述本发明的具体实施方式。

[0061] 实施例1：

[0062] 步骤1，找到所有投影方向Ptotal，对一个16×16待重建图像，探测器分辨率为DetRowNum＝136，则求解得到所有投影方向Ptotal为，

[0063]

[0064] 步骤2，针对寻找到的所有投影方向Ptotal，找到投影方向(p,q)上只穿过一个像素的投影射线；投影方向通过互质整数(p,q)表示，p表示图像水平方向上的整数位移，q表示-1图像垂直方向上的整数位移，投影方向(p,q)表达的投影角度为θ＝tan q/p，如图9所示。图
9中，对于一个9×9待重建图像，投影方向(5,4)上某一条投影射线穿过两个像素，投影方向(3,4)上某一条投影射线穿过三个像素。图10中，灰色区域的虚线表示在投影方向(5,4)上这些投影射线只穿过一个像素，当投影方向(5,4)作为已知投影方向时，该方向上仅有一个像素的投影射线的实际投影值直接反投影到对应的投影像素坐标中，即灰色区域像素的灰度值只需一次迭代即可求解得到。

[0065] 步骤3，计算出每个投影方向上仅穿过一个像素的射线个数后，确定第一个投影方向(5,4)。知道了一个投影方向上只穿过一个像素的投影射线，接下来要找到所有投影方向Ptotal中哪一个投影方向上只穿过一个像素的投影射线个数最多，即射线个数最多的投影方向就是第一个投影方向，如图11(a)中，投影方向(5,4)上投影射线个数最多，为20＝5×4；图11(b)中，投影方向(4,3)上投影射线个数为12＝4×3；图11(c)中，投影方向(5,1)上投影射线个数为5＝5×1；

[0066] 步骤4，找到第一个投影方向后，确定第二个投影方向(-5,4)。将第一次迭代重建出来的像素作为已知信息，更新投影数据，然后再次遍历所有投影方向Ptotal，找到第二次迭代步骤中投影方向上只穿过一个像素的最多投影射线，记为第二个投影方向(-5,4)，将第二个投影方向(-5,4)上单个投影像元直接反投影到对应的像素坐标。

[0067] 步骤5，重复迭代步骤，求解得到所有像素的灰度值，并找到最优投影方向集合P1＝{(5,4),(-5,4),(5,1),(1,5)}，最优投影方向集合P1仅需7次投影更新就能完成所有像素求解，任意给出一组非最优投影集都将需要更多投影更新才能完成所有像素求解。例如，非最优投影集P2需要109次投影更新。

[0068] 步骤6，给出最优投影方向集P1和非最优投影方向集P2在每一次投影更新中采用的投影方向，投影方向集P1总共需要7次投影更新，(5,4)→(-5,4)→(5,1)→(5,4)→(-5,4)→(1,5)→(5,4)；投影方向集P2总共需要109次投影更新，(5,1)→(-4,1)→(3,1)→(5,1)→(3,1)→(5,1)…→…(5,1)→(3,1)。

[0069] 步骤7，给出最优投影方向集合P1的重建路径和非最优投影方向集合P2的重建路径；最优重建路径如图12，非最优重建路径如图13。为了方便描述，给定一个16×16的待重建图像，如图12，此刻每个离散单元上的灰度值代表着该点被求出时所在的迭代步数。例如，图12中迭代次数为1的栅格点，则表示第1步迭代步骤。图12中，连通线段上迭代次数为1，3，5，6，7的离散像素依次被求解得到。首先，迭代次数为1的离散像素通过第一个投影方向(5,4)一次投影更新求得的，该离散像素上噪声累积一次。同样，迭代次数为2的离散像素通过第二个投影方向(-5,4)一次投影更新求得，该离散像素上噪声累积一次。迭代次数为3的离散像素从第三个投影方向(5,1)投影更新求解得到，因为该投影方向穿过三个离散像素，迭代次数为1和迭代次数为2在第一次迭代和第二次迭代先后被求解得到，因此，通过减去迭代次数为1的离散像素和迭代次数为2的离散像素求解得到迭代次数为3的离散像素，最后迭代次数为3的离散像素噪声累积4次。同样，依次求解得到迭代次数为4，5，6，7的离散像素，并计算出迭代次数为7的离散像素上噪声累积17次。又因为Mojette成像模型为平行束成像，重建路径是对称的，且对称的重建路径上每一离散像素的求解过程完全一样。

[0070] 图13分析了非最优投影集P2的重建路径，非最优重建路径需要109次投影更新才能完成所有像素求解。与最优重建路径相同，非最优重建路径上迭代次数为1的栅格点在投影方向(5,1)上，由第一次投影更新中被求解得到；迭代次数为1和迭代次数为3的点在投影方向(3,1)上，当迭代次数为1的像素灰度值已知，更新投影，求解得迭代次数为3的像素灰度值，该像素上噪声累积了2次，依次求解迭代次数为109的像素灰度值，该像素上的误差累积次数远远高于投影更新次数，同时误差被放大很多倍，求解得到的灰度值近似崩溃。

[0071] 实施例2：

[0072] 实施例1依据技术方案中迭代步骤图了解算法的每一步迭代步骤，实施例2在实施例1的基础上，则给出实际探测器分辨率的模拟仿真重建结果。

[0073] 基于以上的步骤说明，在实际尺寸的模拟仿真实验中，对于分辨率大小为1024的探测器来说，选择的投影方向及重建图像大小需要满足(|p|+|q|)·(M-1)+1<1024，∑iI|pi|≥M，利用这两个不等式计算可知，当最大可重建64×64大小的图像断层时，有|p|+|q|≤16，投影集Ptotal包含162个投影，162个投影的具体值如下，为了便于存储和显示，此处将投影集Ptotal中的p和q分别存储在一维向量p和向量q中，取p中的第一个值p(1)和q中的第一个值q(1)，构成投影方向(p1,q1)，同理取p中的第i个值p(i)和q中的第一个值q(i)，可构成另一对投影方向(pi,qi)。一维向量p中各值为：

[0074] p＝[0,1,1,2,3,3,4,4,5,5,5,5,6,6,7,7,7,7,7,7,8,8,8,8,9,9,9,9,9,10,10,10,11,11,11,11,11,12,13,13,13,14,15,1,1,2,1,3,1,2,3,4,1,5,1,2,3,4,5,6,1,3,5,
7,1,2,4,5,7,1,3,5,1,2,3,4,5,1,1,2,3,1,1,-1,-2,-3,-3,-4,-4,-5,-5,-5,-5,-6,-6,-
7,-7,-7,-7,-7,-7,-8,-8,-8,-8,-9,-9,-9,-9,-9,-10,-10,-10,-11,-11,-11,-11,-11,-
12,-13,-13,-13,-14,-15,-1,-1,-2,-1,-3,-1,-2,-3,-4,-1,-5,-1,-2,-3,-4,-5,-6,-
1,-3,-5,-7,-1,-2,-4,-5,-7,-1,-3,-5,-1,-2,-3,-4,-5,-1,-1,-2,-3,-1,-1]；

[0075] 与一维矢量q中的各值为：

[0076] q＝[1,0,1,1,1,2,1,3,1,2,3,4,1,5,1,2,3,4,5,6,1,3,5,7,1,2,4,5,7,1,3,5,1,2,3,4,5,1,2,3,1,1,1,2,3,3,4,4,5,5,5,5,6,6,7,7,7,7,7,7,8,8,8,8,9,9,9,9,9,10,
10,10,11,11,11,11,11,12,13,13,13,14,15,1,1,1,2,1,3,1,2,3,4,1,5,1,2,3,4,5,6,1,
3,5,7,1,2,4,5,7,1,3,5,1,2,3,4,5,1,2,3,1,1,1,2,3,3,4,4,5,5,5,5,6,6,7,7,7,7,7,
7,8,8,8,8,9,9,9,9,9,10,10,10,11,11,11,11,11,12,13,13,13,14,15]；

[0077] 但在重建过程中，是不需要这些投影全部都参与重建的，需要去除一些冗余，在满足Katz引理的前提下，只需选取其中的一部分即可。根据技术方案步骤，选出的最优投影方向为(11,5)，(15,1)，(-9,7)，(-1,15)，(9,7)，(5,11)，(-13,3)，(-15,1)，(-5,11)，(1,15)。下面就选出的最优路径与任意非最优重建路径下的10个投影方向(8,7),(8,3),(8,5),(8,
1),(7,4),(7,3),(7,1),(4,5),(1,1),(7,2)的重建结果进行对比，图14(a)是原图；图(b)(c)是二组非最优重建路径下的重建结果；图(d)是最优重建路径在噪声取值范围[-0.15％Mmax,0.15％Mmax]的重建结果；图(e)是最优重建路径在噪声取值范围[-0.1％Mmax,0.1％Mmax]的重建结果。通过对比，可以发现，重建路径的不同对于重建结果有着很大的影响，最优重建路径算法避免了极小的噪声也会在求解的迭代过程中不断被放大到一个很大数值的问题，保证噪声条件下也能重构出效果最好的重建图像。