CN103745440B - Ct系统金属伪影校正方法 - Google Patents

Abstract

本发明提供了一种CT系统金属伪影校正方法，根据初始化参数创建所述CT系统的系统矩阵，基于所述系统矩阵计算CT图像重建的过程中的物体衰减系数f，再优化所述物体衰减系数f，实现对CT系统的金属伪影校正。本发明提供的金属伪影校正算法，能够在SART算法的迭代重建过程中，通过逐网格计算衰减系数f的方法，抑制因投影数据跃变和射束硬化引起的金属伪影，提高重建效果。

Description

CT系统金属伪影校正方法

技术领域

本发明涉及CT系统图像重建方法技术领域，尤其是涉及一种CT系统金属伪影校正方法。

背景技术

CT成像技术发展迅速，其具有空间分辨率高、扫描速度快、射线利用率高和图像具有各向同性等优势，目前广泛应用于人体器官成像、疾病诊断、术后评估等方面。但是，CT成像过程中，若被扫描物体存在高衰减系数的物质(如人体内的金属器件、病人必须携带的活检剂和芯片中的金属等)会导致严重的金属伪影。金属伪影在CT图像中表现为条纹状亮纹、暗带或阴影。金属伪影严重破坏了CT图像的真实性，极大地限制了CT系统的检测精度和临床应用。因此，非常有必要去除或最大程度地校正金属伪影。

目前，金属伪影校正常用三种方法：射束过滤、双能系统和算法校正。基于迭代重建的算法校正是目前最常用的金属伪影校正方法。但是，上述基于迭代重建的金属伪影校正算法需要基于某些先验信息，很多情况下，这些先验信息是未知的。

发明内容

本发明的目的是：提供一种无需先验信息的CT系统金属伪影校正方法。

本发明的技术方案是：一种CT系统金属伪影校正方法，包括下述步骤：

初始化参数，其中，所述参数包括物体的衰减系数、最大投影角度Θ、重建图像大小N_x×N_y×N_z、扫描物体大小为L_x×L_y×L_z、CT探测器数量D及单个探测单元的大小D_s；

创建所述CT系统的系统矩阵θ为投影角度，t是投影角度的序号，i是该投影角度下第i条X射线；

基于所述系统矩阵计算CT图像重建的过程中的物体衰减系数f；

优化所述物体衰减系数f。

在一些较佳的实施例中，其中，初始化参数包括下述步骤：

根据投影角度，确定所述最大投影角度Θ和投影角度数量n；

根据所述探测器单元数量D及单个探测单元的大小Ds确定需要计算的X射线条数；

根据所述重建图像大小N_x×N_y×N_z和扫描物体大小为L_x×L_y×L_z，确定网格的大小和数量；

初始化重建前的物体衰减系数。

在一些较佳的实施例中，其中，创建所述CT系统的矩阵包括下述步骤：

输入所述初始化参数；

确定投影角度θt；

将需要计算衰减系数的射线记作l_ti(1≤t≤n,1≤i≤D)，其中，t是投影角度的序号，n是投影角度数量，j是该投影角度下第j条X射线，D是探测器单元数量；

计算X射线l_ti经过的网格以及所述l_ti与所述网格的交线长度，得到所述CT系统的系统矩阵

所述网格的网格序号保存在内存M_ij中，所述交线长度保存在内存M'_ij中，所述M_ij和M'_ij中的数据一一对应；

若t=n，完成所述CT系统的矩阵若t≠n，则在投影角度θ_t+1下重复执行上述步骤。

在一些较佳的实施例中，其中，基于所述系统矩阵计算CT图像重建的过程中的物体衰减系数f，包括下述步骤：

在当前投影角度θ_t时，读取探测器i探测到的X射线l_ti透射后的能量；

根据通过探测器探测到的能量计算X射线经过的网格的衰减系数；

在所述X射线l_ti经过的路径上，计算l_ti穿过所述网格g_i,j后的能量，并根据所述能量计算g_i,j的衰减系数，计算的公式为其中，E_j+1是射线穿过网格j+1透射后的能量，E₀是射线的初始能量，f_j’是物体网格g_i,j的衰减系数，η是变换系数；

所述X射线l_ti穿过所述网格g_i,j后，计算透射后X射线的能量；

将穿过网格g_i,j的X射线能量作为下一个网格g_i,j+1的入射X射线能量，计算所述X射线穿过g_i,j+1后的能量和g_i,j+1的衰减系数；

重复上述步骤，直至X射线完全穿透物体，计算最终透射X射线的能量和X射线路径上最后一个网格的衰减系数。

在一些较佳的实施例中，其中，根据所述能量计算g_i,j的衰减系数，为采用下述公式计算g_i,j的衰减系数：

$f_i^{(k+1)}=f_j^{\left(k\right)}+\mathrm{\λ}\frac{\underset{i}{\mathrm{\Σ}}\left(\frac{p_i^{\mathrm{\θ}}-\underset{j}{\mathrm{\Σ}}{M}_{i,j}^{\mathrm{\θ}}\left(\frac{f_j^{\left(k\right)}}{{(1-\mathrm{\α}\underset{j^{\′}=0}{\overset{j-1}{\mathrm{\Σ}}}{f}_{j^{\′}}^{\left(k\right)})}^{\mathrm{\β}}}\right)}{\underset{j}{\mathrm{\Σ}}{M}_{i,j}^{\mathrm{\θ}}}\right)w_{\mathrm{ij}}}{\underset{j}{\mathrm{\Σ}}{M}_{i,j}^{\mathrm{\θ}}}$

上式中，表示系统矩阵，p表示探测器采集的投影数据，f表示计算得到的衰减系数，θ表示投影角度，i表示第i条X射线，即投射到探测器i上的射线；j表示X射线i穿过的第j个网格，k表示迭代次数,α∈[0,1]是校正强度，β∈[2.5,3.5]。在一些较佳的实施例中，其中，优化所述物体衰减系数f，包括下述步骤：

负值校正：将CT图像重建的过程中的物体衰减系数f小于0的衰减系数设置为0；

差异图像计算：计算负值校正后的衰减系数和上一次迭代后经过衰减系数优化后的衰减系数之间的差异；

采用梯度下降流优化计算衰减系数的优化值。

在一些较佳的实施例中，其中，采用梯度下降流优化计算衰减系数的优化值，包括下述步骤：

步骤A：定义k₁=1，k₃=1，计算下述公式

${\stackrel{\→}{v}}_{s,t}^{k_3}=\frac{\∂{\left|\right|f_{s,t}^{k_3}\left|\right|}_{\mathrm{TV}}}{\∂f_{s,t}^{k_3}}$

${\hat{v}}_{s,t}^{k_3}=\frac{{\stackrel{\→}{v}}_{s,t}^{k_3}}{\left|{\stackrel{\→}{v}}_{s,t}^{k_3}\right|}$

$f_{s,t}^{k_3}=f_{s,t}^{k_3}-{\mathrm{ad}}^{k_1}{\hat{v}}_{s,t}^{k_3}$

其中，梯度下降流的迭代次数为k₃，最大迭代次数为K₃；

步骤B：定义k₃:=k₃+1，重复上述步骤A，直至k₃=K₃；

步骤C：定义k₁:=k₁+1，转至步骤A，直至k₁=K₁，算法结束，得到衰减系数的优化值。

在一些较佳的实施例中，初始化后的物体的衰减系数为0；或使用滤波反投影算法重建CT图像，将滤波反投影重建后的结果作为衰减系数的初始值。

本发明的优点是：

本发明提供了一种基于联合代数重建（Simultaneous AlgebraicReconstruction Technique,SART）的金属伪影校正算法，能够在SART算法的迭代重建过程中，通过迭代方式和优化搜索方法改善数据跃变部分的重建效果。在SART算法中逐个依次计算某一条X射线经过的每个体素的衰减系数，对于每个体素，可计算得到经过该体素并衰减后的X射线能量，并将此作为下一个体素的入射射线能量，以此逼近物体特别是金属物体的衰减系数，从而模拟真实的射束硬化现象，降低金属伪影。

附图说明

图1为本发明提供的CT系统金属伪影校正方法的步骤流程图；

图2为本发明一实施例提供的初始化参数的步骤流程图；

图3为本发明一实施例提供的创建所述CT系统的系统矩阵的步骤流程图；

图4为本发明一实施例提供的系统矩阵示意图；

图5为本发明一实施例提供的基于所述系统矩阵计算CT图像重建过程中的物体衰减系数f的步骤流程图；

图6为本发明一实施例提供的本发明计算X射线光子经过物体网格衰减后的光子能量的示意图；

图7为本发明一实施例提供的优化所述物体衰减系数f的步骤流程图。

具体实施方式

请参考图1，图1为本发明实施例提供的CT系统金属伪影校正方法100的步骤流程图，包括下述步骤：

步骤S110：初始化参数，其中，所述参数包括物体的衰减系数、最大投影角度Θ、重建图像大小N_x×N_y×N_z、扫描物体大小为L_x×L_y×L_z、CT探测器数量D及单个探测单元的大小Ds。

请参阅图2，为初始化参数的步骤流程图，包括下述步骤：

步骤S111：根据投影角度，确定所述最大投影角度Θ和投影角度数量n；

步骤S112：根据所述探测器单元数量D及单个探测单元的大小Ds确定需要计算的X射线条数；

步骤S113：根据所述重建图像大小N_x×N_y×N_z和扫描物体大小L_x×L_y×L_z，确定网格的大小和数量；

步骤S114：初始化重建前的物体衰减系数。

可以理解，初始化后的物体的衰减系数为0；或使用滤波反投影算法重建CT图像，将滤波反投影重建后的结果作为衰减系数的初始值。

步骤S120：创建所述CT系统的系统矩阵θ为投影角度，i是该投影角度下第i条X射线；

请参阅图3，为创建所述CT系统的系统矩阵的步骤流程图，包括下述步骤：

步骤S121：输入所述初始化参数；

步骤S122：确定投影角度θ_t；

步骤S123：将需要计算衰减系数的射线记作l_ti(1≤t≤n,1≤i≤D)，其中，t是投影角度的序号，i是该投影角度下第i条X射线；

步骤S124：计算X射线l_ti经过的图像网格以及所述l_ti与所述图像网格的交线长度，得到所述CT系统的系统矩阵

步骤S125：所述图像网格的网格序号保存在内存M_ij中，所述交线长度保存在内存M'_ij中，所述M_ij和M'_ij中的数据一一对应；

步骤S126：若t=n，完成所述CT系统的系统矩阵若t≠n，则在投影角度θ_t+1下重复执行上述步骤。

请参阅图4，为系统矩阵示意图，在图4中，被扫描物体被网格化为6×6的网格，射线l_ti经过了网格f₁、f₇、f₈、f₉、f₁₅、f₁₆、f₁₇、f₂₃和f₂₄，计算X射线经过的网格，将网格序号保存在内存M_ij中；计算X射线l_ij经过这些网格的路径长度（即系统矩阵），并保存在内存M’_ij中。M_ij和M’_ij中的数据存在一一对应关系。

步骤S130：基于所述系统矩阵计算CT图像重建的过程中的物体衰减系数f；

请参阅图5，为基于所述系统矩阵计算CT图像重建的过程中的物体衰减系数f的步骤流程图，包括下述步骤：

步骤S131：在当前投影角度θ_t时，读取探测器i探测到的X射线l_ti透射后的能量；

步骤S132：根据通过探测器探测到的能量计算X射线经过的物体网格的衰减系数；

其中，根据所述能量计算g_i,j的衰减系数，为采用下述公式计算g_i,j的衰减系数：

$f_i^{(k+1)}=f_j^{\left(k\right)}+\mathrm{\λ}\frac{\underset{i}{\mathrm{\Σ}}\left(\frac{p_i^{\mathrm{\θ}}-\underset{j}{\mathrm{\Σ}}{M}_{i,j}^{\mathrm{\θ}}\left(\frac{f_j^{\left(k\right)}}{{(1-\mathrm{\α}\underset{j^{\′}=0}{\overset{j-1}{\mathrm{\Σ}}}{f}_{j^{\′}}^{\left(k\right)})}^{\mathrm{\β}}}\right)}{\underset{j}{\mathrm{\Σ}}{M}_{i,j}^{\mathrm{\θ}}}\right)w_{\mathrm{ij}}}{\underset{j}{\mathrm{\Σ}}{M}_{i,j}^{\mathrm{\θ}}}$

上式中，表示系统矩阵，p表示探测器采集的投影数据，f表示计算得到的衰减系数，θ表示投影角度，i表示第i条X射线，即投射到探测器i上的射线；j表示X射线i穿过的第j个网格，k表示迭代次数，,α∈[0,1]是校正强度，β∈[2.5,3.5]。。

步骤S133：在所述X射线l_ti经过的路径上，计算l_ti穿过所述网格g_i,j后的能量，并根据所述能量计算g_i,j的衰减系数，计算的公式为其中，E_j+1是射线穿过网格j+1透射后的能量，E₀是射线的初始能量，f_j’是物体网格g_i,j的衰减系数，η是变换系数；

请参阅图6为本发明计算X射线光子经过网格衰减后的光子能量的示意图。

步骤S134：所述X射线l_ti穿所述网格g_i,j后，计算透射后X射线的能量；

步骤S135：将穿过网格g_i,j的X射线能量作为下一个网格g_i,j+1的入射X射线能量，计算所述X射线穿过g_i,j+1后的能量和g_i,j+1的衰减系数；

步骤S136：重复上述步骤，直至X射线完全穿透物体，计算最终透射X射线的能量和X射线路径上最后一个网格的衰减系数。

步骤S140：优化所述物体衰减系数f。

请参阅图7，为优化所述物体衰减系数f的步骤流程图，包括下述步骤：

步骤S141：负值校正：将CT图像重建过程中的物体衰减系数f小于0的衰减系数设置为0；

步骤S142：差异图像计算：计算负值校正后的衰减系数和上一次迭代后经过衰减系数优化后的衰减系数之间的差异；

步骤S143：采用梯度下降流优化计算衰减系数的优化值。

具体地，采用梯度下降流优化计算衰减系数的优化值，包括下述步骤：

步骤A：定义k₃=1，计算下述公式

${\stackrel{\→}{v}}_{s,t}^{k_3}=\frac{\∂{\left|\right|f_{s,t}^{k_3}\left|\right|}_{\mathrm{TV}}}{\∂f_{s,t}^{k_3}}$

${\hat{v}}_{s,t}^{k_3}=\frac{{\stackrel{\→}{v}}_{s,t}^{k_3}}{\left|{\stackrel{\→}{v}}_{s,t}^{k_3}\right|}$

$f_{s,t}^{k_3}=f_{s,t}^{k_3}-{\mathrm{ad}}^{k_1}{\hat{v}}_{s,t}^{k_3}$

其中，梯度下降流的迭代次数为k₃，最大迭代次数为K₃；

步骤B：定义k₃:=k₃+1，重复上述步骤A，直至k₃=K₃；

步骤C：定义k₁:=k₁+1，转至步骤A，直至k₁=K₁，算法结束，得到衰减系数的优化值。

可以理解，通过上述步骤S110至步骤S140可以实现对CT系统金属伪影校正。

本发明提供了一种基于联合代数重建（Simultaneous AlgebraicReconstruction Technique,SART）的金属伪影校正算法，能够在SART算法的迭代重建过程中，通过迭代方式和优化搜索方法改善数据跃变部分的重建效果。在SART算法中逐个依次计算某一条X射线经过的每个网格的衰减系数，对于每个网格，可计算得到经过该网格衰减后的X射线能量，并将此作为下一个网格的入射射线能量，以此逼近物体特别是金属物体的衰减系数，从而模拟真实的射束硬化现象，降低金属伪影。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (7)

1.一种CT系统金属伪影校正方法，其特征在于，包括下述步骤：

初始化参数，其中，所述参数包括物体的衰减系数f、最大投影角度Θ、重建图像大小N_x×N_y×N_z、扫描物体大小为L_x×L_y×L_z、CT探测器数量D及单个探测单元的大小D_s；

创建所述CT系统的系统矩阵θ为投影角度，i是该投影角度下第i条X射线；

基于所述系统矩阵计算CT图像重建过程中的物体衰减系数f；

优化所述物体衰减系数f；

其中，基于所述系统矩阵计算CT图像重建过程中的物体衰减系数f，包括下述步骤：

在当前投影角度θ_t时，读取探测器i探测到的X射线l_ti透射后的能量，其中，t是投影角度的序号；

根据通过探测器探测到的能量计算X射线经过的网格的衰减系数；

在所述X射线l_ti经过的路径上，计算l_ti穿过所述网格g_i,j后的能量，并根据所述能量计算g_i,j的衰减系数，其中，计算的公式为其中，E_j+1是射线穿过网格j+1透射后的能量，E₀是射线的初始能量，f_j’是物体网格g_i,j的衰减系数，η是变换系数；

所述X射线l_ti穿过所述网格g_i,j后，计算透射后X射线的能量；

将穿过网格g_i,j的X射线能量作为下一个网格g_i,j+1的入射X射线能量，计算所述X射线穿过g_i,j+1后的能量和g_i,j+1的衰减系数；

重复上述步骤，直至X射线完全穿透物体，计算最终透射X射线的能量和X射线路径上最后一个网格的衰减系数。

2.根据权利要求1所述的CT系统金属伪影校正方法，其特征在于，其中，初始化参数包括下述步骤：

根据投影角度，确定所述最大投影角度Θ和投影角度数量n；

根据所述探测器数量D及单个探测单元的大小D_s确定需要计算的X射线条数；

根据所述重建图像大小N_x×N_y×N_z和扫描物体大小L_x×L_y×L_z，确定网格的大小和数量；

初始化重建前的物体衰减系数。

3.根据权利要求1所述的CT系统金属伪影校正方法，其特征在于，其中，创建所述CT系统的系统矩阵包括下述步骤：

输入所述初始化参数；

确定投影角度θ_t；

将需要计算衰减系数的射线记作l_ti，1≤t≤n,1≤i≤D，其中，t是投影角度的序号，n是投影角度数量，i是该投影角度下第i条X射线，D是探测器数量；

计算X射线l_ti经过的网格以及所述l_ti与所述网格的交线长度，得到所述CT系统的系统矩阵

所述网格的网格序号保存在内存M_ij中，所述交线长度保存在内存M'_ij中，所述M_ij和M'_ij中的数据一一对应；

若t＝n，完成所述CT系统的系统矩阵的创建；若t≠n，则在投影角度θ_t+1下重复执行上述步骤。

4.根据权利要求1所述的CT系统金属伪影校正方法，其特征在于，其中，根据所述能量计算g_i,j的衰减系数，为采用下述公式计算g_i,j的衰减系数：

$f_j^{(k+1)}=f_j^{\left(k\right)}+\mathrm{\λ}\frac{\underset{i}{\Σ}\left(\frac{p_i^{\mathrm{\θ}}-\underset{j}{\Σ}{M}_{i,j}^{\mathrm{\θ}}\left(\frac{f_j^{\left(k\right)}}{{(1-\mathrm{\α}\underset{j^{\′}=0}{\overset{j-1}{\Σ}}{f}_{j^{\′}}^{\left(k\right)})}^{\mathrm{\β}}}\right)}{\underset{j}{\Σ}{M}_{i,j}^{\mathrm{\θ}}}\right)w_{ij}}{\underset{j}{\Σ}{M}_{i,j}^{\mathrm{\θ}}}$

上式中，表示系统矩阵，p表示探测器采集的投影数据，f表示计算得到的衰减系数，θ表示投影角度，i表示第i条X射线，即投射到探测器i上的射线；j表示X射线i穿过的第j个网格，k表示迭代次数，α∈(0,1)是校正强度，β∈(2.5,3.5)。

5.根据权利要求1所述的CT系统金属伪影校正方法，其特征在于，其中，优化所述物体衰减系数f，包括下述步骤：

负值校正：将CT图像重建的过程中的物体衰减系数f小于0的衰减系数设置为0；

差异图像计算：计算负值校正后的衰减系数和上一次迭代后经过衰减系数优化后的衰减系数之间的差异；

采用梯度下降流优化计算衰减系数的优化值。

6.根据权利要求5所述的CT系统金属伪影校正方法，其特征在于，其中，采用梯度下降流优化计算衰减系数的优化值，包括下述步骤：

步骤A：定义k₁＝1，k₃＝1，计算下述公式

$f_{s,t}^{k_3}=f_{s,t}^{k_3}-{\mathrm{ad}}^{k_1}{\hat{v}}_{s,t}^{k_3}$

其中，梯度下降流的迭代次数为k₃，最大迭代次数为K₃；

步骤B：定义k₃:＝k₃+1，重复上述步骤A，直至k₃＝K₃；

步骤C：定义k₁:＝k₁+1，转至步骤A，直至k₁＝K₁，算法结束，得到衰减系数的优化值。

7.根据权利要求1所述的CT系统金属伪影校正方法，其特征在于，初始化后的物体的衰减系数为0；或使用滤波反投影算法重建CT图像，将滤波反投影重建后的结果作为衰减系数的初始值。