CN107274459B - 一种用于加快锥形束ct图像迭代重建的预条件方法 - Google Patents

Abstract

本发明涉及一种用于加快锥形束CT图像迭代重建的预条件方法，通过构建一个3D块三对角预条件算子

，考虑到图像层间的正投影耦合关系，较好的反应正投影算子的属性，使得迭代快速收敛，解预条件方程的浮点计算量仅仅是每个迭代步骤的图像尺寸的3倍。本发明的3D块三对角预条件算子

能更好地反应系统矩阵的结构，显著提高收敛速度。

Description

一种用于加快锥形束CT图像迭代重建的预条件方法

技术领域

本发明属于一般的医疗图像生成的技术领域，特别涉及一种可以有效减少迭代步数、加速锥形束CT图像迭代重建的用于加快锥形束CT图像迭代重建的预条件方法。

背景技术

随着大众对非创伤诊断手段准确性期望的提高，当代医用CT对图像的精确度有越来越高的要求。事实上CT图像的精确性受到许多因素的制约，比如射线硬化，焦外辐射，散射，探测器有限尺寸等因素。所有这些因素都要做适当的校正和理想化近似。然后对校正后的投影数据做图像重建。见“CT原理图”。现在常用医疗CT图像重建方法是滤波反投影(FBP)。这一方法的优点是简单快速。缺点是它对数据有严格的假设题条件。比如数据应该是确定性的而不是具有统计特征的，数据是单色的，即做过精确度射线硬化校正等。这些条件很难精确满足，因此FBP重建图像在有些情况下会有一些伪影。而迭代重建方法的输入数据可以是没有做过射线硬化校正的，也可以是校正后仍然具有统计特性的数据。迭代重建可以建立相应的模型来重建图像。

在X射线计算机断层摄影(CT)成像中，迭代图像重建是降低辐射剂量并保持良好图像质量以满足放射科医师用于诊断目的的有效途径。迭代重建可以在低剂量条件下，产生较高分辨率和较少图像伪影的图像。

然而，在实际的迭代重建计算过程中，迭代重建的计算速度很慢，经常大于半小时，而医院的临床医生希望得到实时图像。不仅是为了提高设备使用效率也是急诊的需要。速度慢阻止了迭代重建替代常规的滤波反投影(FBP)分析方法，成为广泛使用的方法，而众所周知，统计迭代重建的图像质量要远优于滤波反投影图像重建的图像质量。迭代重建常用方法是共轭梯度法(CG)，见图“CG迭代算法”。因为要反复迭代许多步才收敛，而且每一步的迭代都需要做大量计算来进行正投影和反投影。相比而言FBP的主要计算就一次反投影。因此如何加快迭代重建是这一领域的重点和难点。

当然使用高级并行计算硬件是可以加速的，但也是有限度的。另一个途径是设计算法加速。一方面可以用算法加速正投影和反投影的计算，另一方面是用新算法减少迭代步数。算法加速不仅节省设备成本而且避免高级硬件所需要的电力能耗。本发明用于减少迭代步数。

有几个尝试通过预条件加速收敛的数值算法已经被提出。图像重建问题早期设计的预条件算子是平移不变的，它适用于某些特殊情况，但对移变投影无效，特别是当权重涉及统计优化问题时没有显著加速收敛的效果。对于2D重建的问题已经有了有效的移变预条件算法。然而，锥束3D CT重建的预条件算子设计具有更多的挑战。本发明构建了一个3D块三对角预条件算子。所提出的预条件算子是具有移变特性，计算量小而且加速效果显著。

发明内容

现有技术中，迭代重建的收敛速度很慢，而导致的阻止了其替代常规的分析方法成为广泛使用的方法，而众所周知，统计迭代重建的图像质量要远优于滤波反投影图像重建的图像质量。本发明解决的技术问题是，提供了一种优化的用于加快锥形束CT图像迭代重建收敛的预条件方法。

本发明所采用的技术方案是，一种用于加快锥形束CT图像迭代重建的预条件方法，所述方法包括以下步骤：

步骤1.1：采集CT原始投影数据y，校正，通过迭代重建求解系统方程(P^TP+μC^TC)x＝P^Ty，重建3D图像x.其中，P是CT图像的正投影算子，C是微分算子，μ是预先设定的正则化参数；

步骤1.2：令矩阵K＝(P^TP+μC^TC)；

步骤1.3：分解K的子块

其中，F_2D是二维傅里叶变换算子，Λ_ij为对角矩阵；构造K的块三对角预条件矩阵M，

其中，E_k＝Λ_k,k+1＝Λ_k+1,k,；

步骤1.4：令

分解S＝LDL^T，其中，

步骤1.5：在每一次迭代步中，利用傅立叶变换和LDL^T分解来求解预条件方程Mz＝r，其中，r为一个已知的3D图像，求解z；

步骤1.6：输出z。

优选地，所述CT的X射线为三维锥形束，相邻图像层正投影耦合。

优选地，所述步骤1.3中，E_k-1和Λ_k的计算包括以下步骤：

步骤3.1：生成冲击函数△_k，满足第k层除了中心为1，其余元素为0；

步骤3.2：得到Λ_k＝diag{|F_2D·χ_k·K·χ_k·△_k|}；其中，χ_k限制3D图像在第k层的2D图像；

步骤3.2：得到E_k-1＝diag{|F_2D·χ_k-1·K·χ_k·△_k|}；其中，χ_k限制3D图像在第k层的2D图像，χ_k-1限制3D图像在第k-1层的2D图像。

优选地，所述步骤3.1中，第k层包括J*J个元素；当J为偶数时，则第

行的第

个元素为1，其余元素为0；当J为奇数，则第

行的第

个元素为1，其余元素为0。

优选地，所述D₁＝Λ₁。

优选地，当k＝2,…N时，所述L_k-1＝E_k-1/D_k-1，D_k＝Λ_k-E_k-1L_k-1。

优选地，所述步骤1.5中，在每一次迭代步包括以下步骤：

步骤5.1：设置预条件方程Mz＝r；

步骤5.2：建立三个方程

Dx₂＝x₁、Lx₃＝x₂，通过分别求解带宽上块三对角系统L^T、块对角系统D和带宽下块三对角系统L来完成预条件方程Mz＝r的求解；

步骤5.3：解得x₃后，以

进行傅里叶变换，得到z。

本发明提供了一种优化的用于加快锥形束CT图像迭代重建的预条件方法，通过构建一个3D块三对角预条件算子M，考虑到图像层间的正投影耦合关系，较好的反应正投影算子的属性，使得迭代快速收敛，解预条件方程的浮点计算量仅仅是每个迭代步骤的图像尺寸的3倍。本发明的3D块三对角预条件算子M能更好地反应系统矩阵的结构，显著提高收敛速度。

附图说明

图1为从头部模体的轴向锥形束扫描数据测试，图像层间隔0.625mm、每个图像层具有512*512元素、重建FOV设为250mm时，现有技术的共轭梯度法1和本发明方法2的9次迭代的标准误差；其中，线条1为共轭梯度法的数据，线条2为本发明方法的数据，纵坐标数据为标准误差，横坐标数据为迭代次数。

具体实施方式

下面结合实施例对本发明做进一步的详细描述，但本发明的保护范围并不限于此。

本发明涉及一种用于加快锥形束CT图像迭代重建的预条件方法，所述方法包括以下步骤：

步骤1.1：采集CT原始投影数据y，校正，通过迭代重建求解系统方程(P^TP+μC^TC)x＝P^Ty，重建3D图像x.其中，P是CT图像的正投影算子，C是微分算子，μ是预先设定的正则化参数。

本发明中，若没有正则化项，则系统具有非常小的特征值，正则项使系统不仅保持对称正定，且在数值上具有良好的条件数特性。

本发明中，校正主要包括暗电流校正、空气校正、串扰校正、射线硬化校正和探测器均匀性校正。校正项为本领域技术人员容易理解的内容，可以依本领域技术人员的需求自行处理。

步骤1.2：令矩阵K＝(P^TP+μC^TC)。

步骤1.3：分解K的子块

其中，F_2D是二维傅里叶变换算子，Λ_ij为对角矩阵；构造K的块三对角预条件矩阵M，

其中，E_k＝Λ_k,k+1＝Λ_k+1,k,。

本发明中，矩阵K可以被拆分为

K的主要部分是对角带宽部分，且K_kk≈K_k+1k+1，但，当k和j不接近时，K_kk与K_jj不接近。分解K的子块

其中，F_2D是二维傅里叶变换算子，Λ_ij为对角矩阵；计算得到

其中，E_k＝Λ_k,k+1＝Λ_k+1,k。

本发明中，扩展块对角预条件算子，包括非对角线块，以反映CT图像的正投影算子的空间特性。

本发明中，步骤1.3利用预条件算子M的Toeplitz近似结构做块循环矩阵逼近，从而可以用二维Fourier变换加速计算。

所述步骤1.3中，E_k-1和Λ_k的计算包括以下步骤：

步骤3.1：生成冲击函数△_k，满足第k层除了中心为1，其余元素为0；

步骤3.2：得到Λ_k＝diag{|F_2D·χ_k·K·χ_k·△_k|}；其中，χ_k限制3D图像在第k层的2D图像；

步骤3.2：得到E_k-1＝diag{|F_2D·χ_k-1·K·χ_k·△_k|}；其中，χ_k限制3D图像在第k层的2D图像，χ_k-1限制3D图像在第k-1层的2D图像。

所述步骤3.1中，第k层包括J*J个元素；当J为偶数时，则第

行的第

个元素为1，其余元素为0；当J为奇数，则第

行的第

个元素为1，其余元素为0。

本发明中，预条件算子M的构造过程中，应用了稀疏点源加插值的办法来降低计算量，不需要独立计算每个块，而是可以沿着Z轴，在具有适当间隔的图像层处，分配冲击函数△_k，只进行K运算，不进行χ运算，创建代表块后，使用插值形成其他块。

本发明中，冲击函数△_k的中心为1，只要靠近中心的点为1即可。

本发明中，充分利用系统的对称正定和TOEPLITZ结构特点开发具有实质加速效果而且计算简单的预条件算子M，M的子块是随空间位置而变化的，使用块三对角结构使得预条件算子M精巧地体现了系统矩阵K的结构特点，所以加速效果好而且求解快速简单。

步骤1.4：令

分解S＝LDL^T，其中，

所述D₁＝Λ₁。

当k＝2,…N时，所述L_k-1＝E_k-1/D_k-1，D_k＝Λ_k-E_k-1L_k-1。

本发明中，S＝LDL^T为块三角分解。

本发明中，E_k-1/D_k-1是对角元素的逐个相除。

步骤1.5：在每一次迭代步中，利用傅立叶变换和S＝LDL^T分解来求解预条件方程Mz＝r，其中，r为一个已知的3D图像，求解z。

本发明中，在迭代的每个预条件步骤中，都需要求解预条件方程Mz＝r，根据其结构特点，M可以分解成几个简单矩阵相乘，特别是对它的核心部分S的分解得到S＝LDL^T，随后在每个迭代步只需要求解几个简单的方程，最终以一次傅里叶变换完成z的求解。

步骤1.6：输出z。

所述CT的X射线为三维锥形束，相邻图像层正投影耦合。

本发明中，S＝LDL^T的分解在迭代前完成，并在所有的迭代步骤中反复使用。当要重建的3D图像有N_t＝J*J*N个像素时，求解预条件方程的成本为3倍像素的存储量和浮点计算。具体地讲，仅为3N_t+2*N*log₂(J*J)浮点计算，不需要额外的存储。

本发明不仅可以用于加速锥束CT迭代重建，也可用于加速其它类型CT迭代重建，以及PET，SPECT等医疗图像迭代重建。

本发明中，采用MinFound CT64从头部模体的轴向锥形束扫描数据测试本方法，图像层间隔0.625mm，每个图像层具有512*512元素，重建FOV设为250mm。比较通用的共轭梯度法1和本发明方法2的收敛速度，以200次迭代的图像作为参考图像，如图1，说明了9次迭代的标准误差，显示预条件算子M的有效性。

本发明解决了现有技术中，迭代重建的计算速度很慢的问题所周知，统计迭代重建的图像质量要远优于滤波反投影图像重建的图像质量的问题。本发明通过构建一个3D块三对角预条件算子M，考虑到图像层间的正投影耦合关系，较好的反应正投影算子的属性，使得迭代快速收敛，每个迭代步骤的预条件计算成本仅仅是图像尺寸的3倍。本发明的3D块三对角预条件算子M能更好地反应系统矩阵的结构，显著提高收敛速度。

Claims (7)

1.一种用于加快锥形束CT图像迭代重建的预条件方法，其特征在于：所述方法包括以下步骤：

步骤1.1：采集CT原始投影数据y，校正，通过迭代重建求解系统方程(P^TP+μC^TC）x＝P^Ty，重建3D图像x.其中，P是CT图像的正投影算子，C是微分算子，μ是预先设定的正则化参数；

步骤1.2：令矩阵K＝(P^TP+μC^TC)；

步骤1.3：分解K的子块

其中，F_2D是二维傅里叶变换算子，Λ_ij为对角矩阵；构造K的块三对角预条件矩阵M，

其中，E_k＝Λ_k，k+1＝Λ_k+1，k；

步骤1.4：令

分解S＝LDL^T，其中，

步骤1.5：在每一次迭代步中，利用傅立叶变换和LDL^T分解来求解预条件方程Mz＝r，其中，r为一个已知的3D图像，求解z；

步骤1.6：输出z。

2.根据权利要求1所述的一种用于加快锥形束CT图像迭代重建的预条件方法，其特征在于：所述CT的X射线为三维锥形束，相邻图像层正投影耦合。

3.根据权利要求2所述的一种用于加快锥形束CT图像迭代重建的预条件方法，其特征在于：所述步骤1.3中，E_k-1和Λ_k的计算包括以下步骤：

步骤3.1：生成冲击函数Δ_k，满足第k层除了中心为1，其余元素为0；

步骤3.2：得到Λ_k＝diag{|F_2D·x_k·K·x_k·Δ_k|}；其中，x_k限制3D图像在第k层的2D图像；

步骤3.2：得到E_k-1＝diag{|F_2D·x_k-1·K·x_k·Δ_k|}；其中，x_k限制3D图像在第k层的2D图像，x_k-1限制3D图像在第k-1层的2D图像。

4.根据权利要求3所述的一种用于加快锥形束CT图像迭代重建的预条件方法，其特征在于：所述步骤3.1中，第k层包括J*J个元素；当J为偶数时，则第

行的第

个元素为1，其余元素为0；当J为奇数，则第

行的第

个元素为1，其余元素为0。

5.根据权利要求1所述的一种用于加快锥形束CT图像迭代重建的预条件方法，其特征在于：所述D₁＝Λ₁。

6.根据权利要求1所述的一种用于加快锥形束CT图像迭代重建的预条件方法，其特征在于：当k＝2，…N时，所述L_k-1＝E_k-1/D_k-1，D_k＝Λ_k-E_k-1L_k-1。

7.根据权利要求1所述的一种用于加快锥形束CT图像迭代重建的预条件方法，其特征在于：所述步骤1.5中，在每一次迭代步包括以下步骤：

步骤5.1：设置预条件方程Mz＝r；

步骤5.2：建立三个方程

Dx₂＝x₁、Lx₃＝x₂，通过分别求解带宽上块三对角系统L^T、块对角系统D和带宽下块三对角系统L来完成预条件方程Mz＝r的求解；

步骤5.3：解得x₃后，以

进行傅里叶变换，得到z。

Images