基于投影矩阵广义逆的CT重建迭代算法
王会;徐亚楠;武王将;石宏理;杨智;罗述谦
【摘 要】Objective In the CT examination, the limited-angle projection and sparse projection methods can be used to reduce the X-ray dose.However, these methods would result in the deficiency of the projection data and bringing certain difficulty for image reconstruction.For the purpose of overcoming the challenges and obtaining a higher quality reconstruction image, a new CT reconstruction method is proposed in this paper.Methods In the proposed method, the reconstruction image was expressed as the multiplication of projection matrix and its the Moose-Penrose pseudo-inverse.The first-order recursive formula was used to calculate the pseudo-inverse.However, the projection matrix and its the Moose-Penrose pseudo-inverse were replaced by projection and filtered back projection in the iterative process because of the difficulty of calculation.Then, the data of the parallel beam projection, the limited-angle projection and the sparse projection were reconstructed by different methods, and the mean square error (MSE), universal quality index (UQI) and mutual inform
ation (MI) were used to compare the results of the reconstruction images.Results The results of the MSE, UQI and MI by using the proposed method were better, and the time was less.Conclusions The proposed method can improve the quality of the reconstruction image efficiently without any prior structure information or an artificial assumption on the underlying image.Neither does it need the expression of projection process.The implementation of the proposed method is simple.%目的 在CT检查时,有限角度投影和稀疏矩阵投影能够减少X射线的剂量,然而这会导致投影数据不足,给图像重建带来一定的困难.为了克服这一难题得到较好的重建图像,本文提出一种基于计算投影矩阵广义逆的CT迭代重建算法. 方法 该算法在计算过程中,将重建图像表示为投影矩阵以及其广义逆的乘积.首先使用一阶迭代计算广义逆矩阵,但是由于投影矩阵和其广义逆矩阵都比较大,在迭代过程中以投影和滤波反投影代替.然后通过不同的算法分别对平行束投影、有限角度投影、稀疏矩阵投影的数据进行重建,并对重建结果的均方差、通用图像质量指标以及图像互信息进行比较. 结果 本文提出的方法重建出图像的均方差、通用图像质量指标和图像互信息更优,而且重建时间较短. 结论 该方法能够在没有未知图像先验结构信息和伪影猜想的情况下有效地提高重建图像的质量,而且该算法不需要计算投影过程,重建过程简单易行.
【期刊名称】《北京生物医学工程》
【年(卷),期】2017(036)003
【总页数】6页(P251-256)
【关键词】CT图像重建;迭代算法;广义逆矩阵;稀疏矩阵;有限角度
【作 者】王会;徐亚楠;武王将;石宏理;杨智;罗述谦
【作者单位】首都医科大学生物医学工程学院 北京 100069;首都医科大学生物医学工程学院 北京 100069;首都医科大学生物医学工程学院 北京 100069;首都医科大学生物医学工程学院 北京 100069;首都医科大学生物医学工程学院 北京 100069;首都医科大学生物医学工程学院 北京 100069
【正文语种】中 文
【中图分类】R318.04
计算机断层成像(computed tomography,CT)以非侵入的方式获得人体内部不层叠的解剖结构的信息，然后通过对这些信息进行重建，得到相应的CT图像，为医生的诊断提供依据[1-3]。该图像的重建效果取决于三个因素：重建算法、投影数据的完整性及其质量。当投影数据完整且无噪声干扰时，可以使用解析算法——滤波反投影(filtered back projection,FBP)算法快速地重建出原始图像，例如平行束投影，当对物体进行全方位0°～ 180°扫描时，可以重建出较高质量的图像。但是全方位扫描所需要的X射线的剂量较大，X射线虽然可以杀死病变组织，但同时对于正常的组织也有损害，根据国际辐射防护委员会(International Commission on Radiological Protection,IRCP)研究认证，做一次CT全身扫描检查，会使受检者辐射致癌的危险性增加约8%，因此要尽可能地减少临床上X线的剂量，而这必然会使得投影数据不全，从而导致重建图像中存在条纹状伪影，图像细节信息模糊。此外，在工业中，由于一些实际情况的约束，导致扫描角度的范围变小，在这种有限角度的情况下重建图像会在部分方向上存在严重伪影[4-6]。因此，希望改进的迭代算法能够较好实现不完全投影数据的重建，并尽量抑制伪影的产生。
常用的CT重建算法可以分成两类：解析重建(analytical reconstruction，AR)算法和迭代算法(iterative reconstruction,IR)。AR算法，例如FBP算法以及其改进算法锥形束重建算法，是应
用较为广泛的解析算法，简单且计算速度快，但是对采集数据的完整性要求较高，对于不完全投影数据，用该方法进行图像重建得到的结果比较差。而IR算法，虽然运算量大，花费时间较长，但是对于不完全投影和带有噪声的投影数据而言，迭代算法依然被认为是一种较好的重建方法。常用的迭代算法一般分为两种：代数重建算法(algebraic reconstruction technique,ART)和统计迭代重建算法(statistical iterative reconstruction,SIR)[7-8]。ART算法以及它的改进算法，联合代数重建算法(simultaneous algorithm reconstruction technique,SART)和同时迭代重建算法(simultaneous iterative reconstruction technique,SIRT)会受到不完整投影以及噪声的影响[9]。如果投影数据不完整或受到噪声干扰，它们的收敛速度就会变得缓慢。相对于ART算法，SIR算法的计算更加灵活，尤其在最大后验概率、压缩感知等的模式中，具有更好的效果[10-13]。但是，SIR算法高度依赖于未知图像的先验信息，例如，正则化约束条件和获取图像的系统模型，这些可能会导致图像质量下降。除此之外，迭代重建算法的终止条件并没有一个明确的标准，这意味着为了得到最佳的重建图像，需要不断地调整终止条件来得到最好的重建结果。
本文提出了一种基于投影矩阵广义逆的迭代重建算法。该算法不需要计算投影矩阵以及图像的先验结构信息，仅有投影数据就可以重建出原图像，它将重建图像表示为投影矩阵的广义
逆矩阵与投影值的乘积，其中，投影矩阵的广义逆确保重建图像是投影方程组的最小范数最小二乘法解，理论上讲，这是一种只通过投影数据就可以获得最佳重建结果的方法。由于投影矩阵及其广义逆矩阵比较大，难以获取，在迭代过程中以投影和(滤波)反投影代替矩阵的直接运算，投影矩阵及其广义逆并不会在计算过程中出现，这样计算时间显著减少。此外，由于所有的投影过程都可以表示成线性方程的形式，因此，该方法可以适用于平行投影、扇形束投影以及锥形束投影等所有的投影方式。
一幅大小为N×N的二维图像的投影过程可以表示为：
式中：I表示大小为NN×1的图像矢量；T表示LM×NN的投影系数矩阵，M表示投影角度，L表示每个投影角度上的最大投影数；P是大小为LM×1的投影矢量。
则其重建图像可以用式(2)表示：
式中表示重建图像；大小为NN×1；T +表示投影系数矩阵T的广义逆。由于矩阵T +比较大，在计算过程中会大大增加重建时间，本文中使用其他相关计算进行代替，避免了T +的直接计算，节省了运算时间。
根据下述原理，使用一阶迭代的方法计算T +。
假设大小为LM×NN的投影系数矩阵T的广义逆T +的初始估计值是X0，如果残差R0=PR(T)-TX0满足ρ(R0)<1，其中，ρ(R0)表示R0的谱半径，PR(T)表示T的正交矩阵，则序列[X0,X1,…，Xk,Xk+1,…]可用式(3)表示[14]：
Xk+1=Xk+X0-X0TXk,k=0,1,…(3)
当k→∞时，该式收敛于T+，则残值的对应序列满足：
对于任何相乘矩阵的范数，都满足Rk=PR(T)-TXk。
为了方便，T+初始值的近似值X0，简单地认为它等于：
式中：T t是T的转置；β是一个实标量。β满足公式:
式中：λ1(TTt)为TT t的最大非零特征值。
由于T和T+比较大，为了避免计算它们，将式(3)两边同时乘上投影值P，得到：
Xk+1P=XkP+X0P-X0TXkP,k=0,1，…(7)
Xk+1P和 XkP可以分别表示第k+1次和第k次的重建图像次和；X0P表示初始图像；X0TXkP表示对投影，然后再进行重建，本文中使用FBP算法重建图像的β倍来替代X0TXkP。当NN<LM时，β=1；当NN>LM时，β<2-NN/LM，式(7)又可以简化为：
根据以上原理，可以将该算法分为以下步骤：
(1) 初始化β和终止条件ε；
(2) k=0,FBP算法重建出图像得到预处理图像
(3) 对进行投影，得到投影值Pk；
(4) 根据Pk，利用FBP重建图像，并将所得到的结果乘以β即为Ir；
(5) 修正图像
(6) 令
(7) 如果Δ>ε，跳到步骤(3),否则，结束循环。
为验证该算法的性能，进行以下实验，在实验中，β的取值均为1。
2.1 平行束投影
通过MATLAB模拟一幅大小为256×256的头部幻影图，投影角度范围θ=0∶0.2∶179.8， 使用Radon函数获取图像的投影值，然后分别使用FBP算法和本文提出的算法MPIPM进行重建。最后，使用均方差(mean square error,MSE)、通用图像质量指标(universal quality index,UQI)、图像互信息(mutual information,MI)3种评价准则来评估两种重建算法的性能[15-16]，其中，MSE越小,UQI和MI越接近于1，说明重建图像越接近于参考图像，重建效果越好。
表1是分别使用FBP算法和MPIPM算法重建图像，并使用3种质量评估方法对重建结果进行评估得到的结果，从表中可以看出，MPIPM算法重建的图像质量比FBP算法重建的图像质量高，随着MPIPM算法的迭代次数增加，MSE、MI和UQI相应的变化如图1所示，可以发现MPIPM算法具有较快的收敛速度；图2是分别使用FBP算法和MPIPM算法获得的重建结果，
通过对比图像可以发现，对于平行束投影的完整数据投影，MPIPM算法重建出的图像具有更高的分辨率，有效地提高了重建效果。正则化的约束条件
2.2 稀疏矩阵投影
通过MATLAB模拟一幅大小为128×128的头部幻影图，投影角度范围θ=0∶2∶178，然后分别使用MPIPM、SART、SIRT以及有序子集期望最大算法(ordered subsets expectation maximization,OSEM)对图像进行重建，并通过MSE、UQI、MI以及算法运行时间比较各算法的性能。在该实验中，皆以FBP算法得到的图像作为初始图像进行计算，实验结果见表2和图3。
对比表2的数据可以得出MPIPM算法重建图像耗时少且图像效果较好。图3是不同IR算法重建图像的MSE随迭代次数的变化，从图中可以看到，SIRT的收敛速度比较慢，SART和OSEM算法虽然收敛速度较SIRT快，但是从表2可见其耗时较长，这是因为它们还需要计算投影矩阵，而且投影矩阵都比较大，消耗较多的时间，而MPIPM算法可以避免投影矩阵的计算，减少时间消耗。综合表2和图3可以得出，MPIPM算法在重建时间和重建效果上都优于其他3种方法。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。