基于系统矩阵Moore-Penrose伪逆的CT和PET迭代重建算法
【专利说明】基于系统矩阵1^1〇〇「6-?611「〇86伪逆的01'和?￡1'迭代重建算法
[0001 ] 1.技术领域:本发明涉及CT (PET)成像技术领域，特别是CT( PET)迭代重建技术领 域，具体是指一种基于投影矩阵Moore-Penrose伪逆(广义逆)的迭代CT(PET)重建技术。 2.
【背景技术】
[0002] X射线计算断层成像(Computed Tomography，CT)因其无损、精确、快速等优点广泛 应用于医疗、科研和工业等领域，成为临床诊断、手术引导、科学研究和质量监控中不可或 缺的工具。高质量和完备的投影数据是重建高质量图像的关键。随着人们对X线辐射剂量越 来越多的关注，低剂量CT(或称绿色CT)成为新的、更迫切的追求目标。为将辐射剂量水平降 到足够低，严格限制扫描密度。稀疏扫描成为降低辐射剂量的主要方式。有些场合受被扫描 物体形状、探测器尺寸及现场扫描空间等条件的限制，如集成电路板成像，长管状物成像和 乳房成像，无法获得满足要求(如平行光(parallel-beam)扫描下180°角范围的扫描）的投 影数据，此时为不完全投影数据CT重建。当投影数据是不完全时，CT重建一般将成为严重的 病态问题，难以精确重建断层图像。近年来，结合先验信息或正则约束（如最小化总变分， total variation minimization)的重建思路成为新的解决途径，其实现通常采用迭代重 建技术。迭代重建技术具有抗噪性好，伪影抑制能力强，方便引入已有的先验信息等优点。 众多学者进行了大量艰苦的研究工作，提出了各种改进算法。但该类算法均存在着一定的、 甚至严重的缺陷，譬如巨大的计算量，重建速度慢、稳定性差等，其实际应用前景甚至理论 基础仍然存在许多争论，使其远未成为CT重建的主流算法。因此，可以认为重建技术在提高 不完全扫描数据重建质量的潜力还远远没有被充分发挥出来，实用性方面存在的问题还远 没解决。本发明专利将以新的角度和方法解决该类算法，也为解决其它类似问题提供新的 途径。 3.
【发明内容】
[0003] 3.1基本算法
[0004] X射线CT的投影过程可以用线性代数模型为：
[0005] X = PI (1)
[0006] 其中I为矢量表示的待重建图像，P为投影矩阵，X是矢量表示的投影数据。由于平 行光(parallel-beam)、扇形光(fan-beam)和锥束光（cone-beam)X射线扫描的投影过程均 可以由公式（1)表示。由于正电子放射断层造影术(positron emission tomography，PET) 的射线是直线传播，其射线发射投影过程也可以由公式(1)表示。
[0007] 根据数学理论，重建图像为
[0008] I = P+X (2)
[0009] 其中P+表示P的Moore-Penrose伪逆(广义逆）。矩阵论中相关理论表明该重建图像 是线性代数模型的最小二乘法最小范数解，是仅根据投影数据重建时，能够得到的最好图 像。
[0010] 实际中，由于其非常大，投影矩阵P是难以获得、存储和使用的，更不要说它的 Moore-Penrose伪逆。如当图像有1024* 1024像素，扫描角度数为360,投影矩阵包括 10243万*360个实数。当采用双精度存储时，大约4.3TB(terabyte)，大于一般工作站存储容 量。因此必须避免直接计算和使用P和P+。应用以下定理可以避免直接计算P+。
[0011 ] 定理:设P e RN ·M矣〇，根据矩阵论有关定理，又设To = PTY〇PT，Y〇 e RN ·M，满足
[0012] p(R〇)=p(Pr(p)-PTo)<1 (3)
[0013]其中p(RQ)是R〇的谱半径(矩阵特征值的模的最大值），PR(P)是P的正交投影矩阵，P T 表示P的转置矩阵，则当时，序列
[0014] Tk+i = Tk+T〇-ToPTk (4)
[0015] 收敛到p+，Bpj^T； =户+。残差| |pR(P)-PTk+i| | < | |PR(p)-PTk| 11 |Pr(p)-pt〇| |<| pR(P)-PTk| I，其中N · 11表示任何乘性范数。
[0016] To为迭代计算P+的初始估计，其选择应满足公式(3)。一种简单选择方法是
[0017] Το = βΡ* (5)
[0018]其中β是正常数满足
[0019]
(6)
[0020]其中λΚΡΡ"是ppt的最大非零特征值。
[0021] 该定理表明通过矩阵的乘、加和减等运算可以近似计算（逼近）矩阵的Moore- Penrose逆。但是，由于投影矩阵P可以大到无法获得或存储，因此，仍然不能直接使用公式 ⑷。
[0022] 现在对公式(4)两边同乘以投影数据X，得到
[0023] Tk+iX = TkX+ToX-ToPTkX (7)
[0024] 其中TkX可以看作厶，=厶，即第k次迭代得到的重建图像;rQX = f。，第0次迭 代得到的重津图像，即初始的重津图像。公式(7)可以表示为
[0025]
(8)
[0026] 公式(8)表示中间重建结果^经过投影得到投影数据对其再用原来重建技术 和参数进行重建得到最后应用公式(8)计算可以得到新的重建结果f A+1。式(8)中的To 表示某种已知的、便于实现的非迭代重建过程。为了满足公式(3)，根据公式(6)，可选择滤 波反投影(Filtered Back Projection，FBP)重建过程作为重建过程(矩阵），记为B(也可 以选择为没有滤波的反投影(Back Projection，BP)重建过程），设α>〇为一个常数，使p (PR(P)-aPB) <1，即满足公式(3)，则迭代重建技术可表示为
[0027] ........
(9)
[0028] 其中&和表示第0步、第k步和第k+Ι步的图像估计值;B表示某种已知的、便 于实现的非迭代重建过程;a>〇为取决于P和B的常数。
[0029] 可以证明，当迭代足够次时，重建图像趋近于在仅有投影数据条件下能够得到的 最佳图像。
[0030] 3.2算法扩展
[0031] 当需要附加约束条件时，所述方法投影过程的线性代数模型可以表示为：
[0032]
(10)
[0033] 其中I、P和X与公式（1 )相同，F表示附加的约束条件，如正则约束，变分 (variation)矩阵（算子）等。根据数学理论，该方程的解近似基于总变分（total variation)最小化的重建结果，其中β表示正则化约束的权重。在此条件下，公式（10)表示