CN102314698A - 基于阿尔法散度约束的全变分最小化剂量ct重建方法 - Google Patents

Abstract

本发明公开了一种基于阿尔法散度约束的全变分最小化剂量CT重建方法，包括以下步骤：(1)利用CT成像设备采集原始投影数据，对采集到的投影数据进行CT图像预重建；(2)将步骤(1)的结果作为初值，并构建CT重建模型；(3)采用迭代过程求解步骤(2)中CT重建模型，对每步迭代后的结果进行全变分图像恢复，同时判断每步迭代后的结果是否满足终止迭代条件：若否，将全变分恢复得到的图像数据作为下一步迭代初值并继续迭代过程；若是，则将当前全变分恢复得到的图像数据作为最终重建结果。该方法可实现根据低剂量的CT成像设备投影数据即可重建并得到高质量的CT重建图像。

Description

基于阿尔法散度约束的全变分最小化剂量CT重建方法

技术领域

本发明涉及一种医用CT图像处理的方法，具体来说涉及一种基于阿尔法(α)散度约束的全变分最小化低剂量CT重建方法。 

背景技术

临床CT成像中，X线辐射会对人体的潜在风险，如增加癌症发生率及余生肿瘤致死率，因此一直受到广泛关注。但根据目前的图像重建方法，为了得到更好的重建图像，多排CT、双源CT已在临床中广泛应用，新的CT设备的使用造成更大的X射线使用剂量，也对人体造成更大危害，因此，如何在保证图像质量的前提下最大限度地降低X射线使用剂量已成为医学CT成像领域的迫切需要。 

为达到上述目的，当前方法中除了先进的硬件方法外，通过降低管电流实现低剂量成像是最常用的方法，但该方法重建的图像质量存在严重退化现象，难以满足临床诊断需求。大量研究表明，优质的低剂量CT图像重建方法是大幅降低X射线辐射剂量的有效途径。其中，基于图像域和投影域的滤波方法在一定条件下能取得较好的图像质量，但其只能从图像和投影数据自身获取信息，其处理的精度依赖于图像数据的噪声和伪影特性。当辐射剂量降低或投影数据采集不足时，图像和投影数据的噪声和伪影特性极为复杂，重建出的图像质量就会很差。 

发明内容

本发明的目的在于提供一种基于阿尔法(α)散度约束的全变分最小化低剂量CT重建方法，该方法可实现根据低剂量的CT成像设备投影数据即可重建并得到高质量的CT重建图像。 

本发明的目的可通过以下的技术措施来实现： 

一种基于阿尔法(α)散度约束的全变分最小化低剂量CT重建方法，包括以下步骤： 

(1)利用CT成像设备采集原始投影数据，对采集到的投影数据进行CT图像预重建； 

(2)将步骤(1)中的CT图像预重建的结果作为初值，并构建CT重建模型； 

(3)采用迭代过程求解步骤(2)中CT重建模型，对每步迭代后的结果进行全变分图像恢复，同时判断每步迭代后的结果是否满足终止迭代条件：若否，将全变分恢复得到的图像数据作为下一步迭代初值并继续迭代过程；若是，则将当前全变分恢复得到的图像数据作为最终重建结果。 

所述步骤(1)中的CT图像预重建过程为：对CT成像设备采集的原始投影数据采用滤波反投影方法进行图像预重建处理或采用快速的迭代算法进行图像预处理重建。 

所述的快速的迭代算法为最大似然期望最大算法或最大后验估计算法。 

所述步骤(2)中的CT重建模型为： 

其中：f为采集到的原始投影数据；Kμ为迭代过程中生成的投影数据，K＝{K_i，j}为系统矩阵，i，j分别为待重建CT图像的像素个数和CT探测器探测元的个数；μ表示待重建的衰减图像；D_α(f‖Kμ)为数据保真项，α∈[0.6，2]；|μ|_TV为全变分正则化项。 

所述步骤(3)中的迭代过程采用期望最大化(EM)迭代方法进行计算。 

所述系统矩阵是从CT成像设备的系统得到的基于面积加权的系统矩阵，或者，基于体素加权的系统矩阵。 

所述步骤(3)中全变分图像恢复的计算过程为：  ${\mathrm{\μ}}_{k+1}^{\mathrm{\α}}={\mathrm{\μ}}_{k+.5}^{\mathrm{\α}}-\mathrm{\α}\stackrel{\‾}{\mathrm{\λ}}{\mathrm{\μ}}_k^{\mathrm{\α}}\mathrm{div}\left(w\right),$ 其中 $\stackrel{\‾}{\mathrm{\λ}}=\mathrm{\λ}/{\mathrm{\Σ}}_i{K}_{i,j};$ μ_k+1为全变分恢复后的图像数据，div为散度算子，w为二元变量。 

所述步骤(4)中迭代终止条件为：相邻两次全变分图像恢复后的重建模型的差小于阈值。所述阈值为0.001。 

本发明对比现有技术，有如下优点： 

1、由于α散度应用于统计推测时对噪声和离群值具有较好的鲁棒性，本方法利 用α散度来测定低剂量投影数据和估计的投影数据之间的偏差，而全变分在图像降噪的同时保持图像的边缘和纹理结构信息。本方法结合以上两点，从而能够实现低剂量CT图像的优质重建； 

2、试验结果显示，本方法能较好的保持图像边缘，并且抑制匀质区域的噪声，同时具有较高的信噪比。 

附图说明

图1是本发明的基于阿尔法散度约束的全变分最小化低剂量CT重建方法的处理流程图； 

图2(a)是Shepp-Logan体模图像； 

图2(b)是对Shepp-Logan体模低剂量投影数据采用Ramp滤波后的FBP重建图像； 

图2(c)是对Shepp-Logan体模低剂量投影数据采用Hann窗滤波后的FBP重建图像； 

图2(d)是对Shepp-Logan体模低剂量投影数据采用本发明方法的重建图像，其中α取值0.6，λ取值0.2； 

图2(e)是对Shepp-Logan体模低剂量投影数据采用本发明方法的重建图像，其中α取值1.0，λ取值0.2； 

图2(f)是对Shepp-Logan体模低剂量投影数据采用本发明方法的重建图像，其中α取值2.0，λ取值0.2的重建图像； 

图3(a)至图3(d)分别是图2(c)至图2(f)重建图像的水平轮廓图； 

图4为性能体模图像； 

图5(a)为图2(d)对应重建过程中目标函数随迭代次数的变化曲线； 

图5(b)为图2(e)对应重建过程中目标函数随迭代次数的变化曲线； 

图5(c)为图2(f)对应重建过程中目标函数随迭代次数的变化曲线。 

具体实施方式

图1示出了本发明的基于阿尔法(α)散度约束的全变分最小化低剂量CT重建方法(简称：αD-TV)的处理流程图，具体过程如下： 

1、利用CT设备采集低剂量投影数据，由于临床CT扫描中，不同的成像部位，剂量使用不同，本方法在使用过程中，采用一般剂量的1/10左右即可达到良好效果， 即使用最小剂量即可。对采集到的投影数据进行CT图像预重建，CT图像预重建过程为：对CT成像设备采集原始投影数据采用滤波反投影方法进行CT图像预重建或采用快速的迭代算法进行CT图像预重建。快速的迭代算法为最大似然期望最大算法或最大后验估计算法，也可采用其他方法进行CT图像预重建。 

接下来，构建基于α散度约束的全变分最小化CT重建模型，并将图像重建预处理后的结果作为CT重建模型优化求解的初值进行求解，本方法优选如下CT重建模型：  $\underset{\mathrm{\μ}\≥0}{\min }{D}_{\mathrm{\α}}\left(f\right|\left|\mathrm{K\μ}\right)+\mathrm{\λ}{\left|\mathrm{\μ}\right|}_{\mathrm{TV}}$ 其中： 

D_α(f‖Kμ)为数据保真项，其形式为：  $D_{\mathrm{\α}}(f,\mathrm{K\μ})=\frac{1}{\mathrm{\α}(1-\mathrm{\α})}{\mathrm{\Σ}}_{i=1}^I[\mathrm{\α}{f}_i+(1-\mathrm{\α}){\left(\mathrm{K\μ}\right)}_i-f_i^{\mathrm{\α}}{\left(\mathrm{K\μ}\right)}_i^{1-\mathrm{\α}}],$ 其中，α∈[0.6，2]；K＝{K_i，j}为系统矩阵，i，j分别为待重建CT图像的像素个数和CT探测器探测元的个数；μ表示待重建的衰减图像；|μ|_TV为全变分正则化项，其表达式为：  ${\left|u\right|}_{\mathrm{TV}}=\underset{w\∈C_0^{\∞}\left(R^N\right),{\left|\right|w\left|\right|}_{\∞}\≤1}{\sup }{\∫}_{R^N}\mathrm{\μ}\mathrm{div}\left(w\right),$ div为散度算子，R，R^N分别表示1维实数空间和N维实数空间；w为二元变量，初始化w⁰为0，对于任意的n≥0，其计算方法为： $w^{n+1}=\frac{w^n+t\▿(\mathrm{\α}\sqrt{2\mathrm{\α}-1}\stackrel{\‾}{\mathrm{\λ}}{\mathrm{\μ}}_k^{\mathrm{\α}}\mathrm{div}\left(w\right)-\frac{\mathrm{\α}}{\sqrt{2\mathrm{\α}-1}}{\mathrm{\μ}}_{k+.5}^{\mathrm{\α}})}{1+t|\▿(\mathrm{\α}\sqrt{2\mathrm{\α}-1}\stackrel{\‾}{\mathrm{\λ}}{\mathrm{\μ}}_k^{\mathrm{\α}}\mathrm{div}\left(w\right)-\frac{\mathrm{\α}}{\sqrt{2\mathrm{\α}-1}}{\mathrm{\μ}}_{k+.5}^{\mathrm{\α}})|},$ 其中，t为利用非单调的梯度投影算法自适应估计的时间步长； 

λ＞0为正则化参数，用于刻画TV正则化的强度，在具体实施过程中根据实际需求设定。上述CT重建模型为最优实施方式，根据需要也可以采用其他重建模型，并结合下面的求解处理过程来得到重建图像。 

3、采用改进的期望最大化(Expectation Maximization，EM)迭代算法求解上述CT重建模型，且对每步迭代后的结果进行全变分图像恢复并得到下一步迭代初 值，同时判断每步迭代后的结果是否满足终止迭代条件：若否，继续迭代过程；若是，则将当前全变分恢复得到的图像数据作为CT重建模型的最终重建结果。 

改进的EM迭代算法形式为： ${\mathrm{\μ}}_{k+.5}^{\mathrm{\α}}={\mathrm{\μ}}_k^{\mathrm{\α}}{\mathrm{\Σ}}_i{K}_{i,j}{(f/{\left(\mathrm{K\μ}\right)}_i)}^{\mathrm{\α}}/{\mathrm{\Σ}}_i{K}_{i,j},$ 其中μ_k是第k步迭代初值，所述步骤(1)得到的结果即为μ₀或μ_init，μ_k+.5是第k步迭代后的结果，f为CT成像设备采集的原始投影数据，K＝{K_i，j}为系统矩阵，i，j分别待重建CT图像的像素个数和CT探测器探测元的个数。全变分图像恢复的形式为： ${\mathrm{\μ}}_{k+1}^{\mathrm{\α}}={\mathrm{\μ}}_{k+.5}^{\mathrm{\α}}-\mathrm{\α}\stackrel{\‾}{\mathrm{\λ}}{\mathrm{\μ}}_k^{\mathrm{\α}}\mathrm{div}\left(w\right),$ 其中 $\stackrel{\‾}{\mathrm{\λ}}=\mathrm{\λ}/{\mathrm{\Σ}}_i{K}_{i,j},$ μ_k+5为待恢复处理的图像；μ_k+1为恢复后的图像；w为二元变量，初始化w⁰为0，对于任意的n≥0，其计算方法为： 

$w^{n+1}=\frac{w^n+t\▿(\mathrm{\α}\sqrt{2\mathrm{\α}-1}\stackrel{\‾}{\mathrm{\λ}}{\mathrm{\μ}}_k^{\mathrm{\α}}\mathrm{div}\left(w\right)-\frac{\mathrm{\α}}{\sqrt{2\mathrm{\α}-1}}{\mathrm{\μ}}_{k+.5}^{\mathrm{\α}})}{1+t|\▿(\mathrm{\α}\sqrt{2\mathrm{\α}-1}\stackrel{\‾}{\mathrm{\λ}}{\mathrm{\μ}}_k^{\mathrm{\α}}\mathrm{div}\left(w\right)-\frac{\mathrm{\α}}{\sqrt{2\mathrm{\α}-1}}{\mathrm{\μ}}_{k+.5}^{\mathrm{\α}})|}$

并利用非单调的梯度投影算法自适应估计时间步长t来加速算法的全局收敛。 

4、判断全变分恢复后的图像是否满足迭代终止条件，若否，使用全变分恢复后的图像数据做为下一步迭代初值，并继续迭代过程，若是，则将全变分恢复后的图像数据作为本重建方法的最终重建结果。迭代终止条件为：相邻两次全变分图像恢复后的重建模型的差小于阈值0.001，阈值的取值根据实际需求设定。 

下面对具体数据采用本方法进行CT图像的重建： 

首先采用如图2(a)所示的修定的shepp-logan体模图像作为本发明的模拟实验对象。将shepp-logan体模图像作为待重建的目标CT图像，由它仿真生成原始数据。体模大小设为512×512，射源到旋转中心和探测器的距离分别为541mm和949mm，旋转角在[0，2π]间采样值为984，探测器单元为1.00mm。通过系统矩阵K得到投影数据(sinogram)。重建过程中使用滤波反投影算法，滤波过程采用汉宁(Hanning)窗进行，其截止频率设为奈奎斯特频率的80％。 

图2(b)至图2(f)描述了本发明方法和其他方法的效果对比。图2(b)为低剂量数据采用Ramp滤波后的FBP重建图像，图2(c)低剂量数据采用hann窗滤波后的FBP重建图像，图2(d)至图2(f)为本法分别采用参数为：(α＝0.6，λ＝0.2)、(α＝1.0，λ＝0.2)、(α＝2.0，λ＝0.2)得到的重建图像。可以对比看出本发明方法保持了良好的边缘并且在计量衰减大的区域有更佳的抑制噪声的效果。 

从图3(a)至3(d)描述了对低剂量投影数据重建图像上述图2(c)至(f)的水平剖面图，可以看出，使用本方法的重建方法所重建的图像的侧面轮廓图无论在背景区域或者还是在边缘区域都更接近于真实图像的侧面轮廓图，使用本方法的重建方法能够更好的克服重建中的病态问题，重建出更接近于真实体模图像的重建图像。 

表1列出了图2(c)至(f)的图像重建时间，可以看出，本方法在α＝0.6时相对于其他两种情况具有较短的重建时间。 

表1 

表2显示了性能体模中两个感兴趣区域(如图4所示)采用Ramp滤波后的FBP重建图像和本方法重建的对比噪声比(contrast-to-noise ratio，CNR)的结果。可以看出，本方法的CNR值较高，说明了病灶与背景区域的对比明显，利于临床诊断。 

表2 

图5显示了本发明方法的单调性特征，表明本方法能够快速的收敛到稳定的值，这对较大CT图像的重建是非常有利的。 

本发明的实施方式不限于此，在本发明上述基本技术思想前提下，按照本领域的普通技术知识和惯用手段对本发明内容所做出其它多种形式的修改、替换或变更，均落在本发明权利保护范围之内。 

Claims (9)

1.一种基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于包括以下步骤：

(1)利用CT成像设备采集原始投影数据，对采集到的投影数据进行CT图像预重建；

(2)将步骤(1)中的CT图像预重建的结果作为初值，并构建CT重建模型；

(3)采用迭代过程求解步骤(2)中CT重建模型，对每步迭代后的结果进行全变分图像恢复，同时判断每步迭代后的结果是否满足终止迭代条件：若否，将全变分恢复得到的图像数据作为下一步迭代初值并继续迭代过程；若是，则将当前全变分恢复得到的图像数据作为最终重建结果。

2.根据权利要求1所述的基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于：所述步骤(1)中的CT图像预重建过程为：对CT成像设备采集的原始投影数据采用滤波反投影方法进行图像预重建处理或采用快速的迭代算法进行图像预处理重建。

3.根据权利要求2所述的基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于：所述的快速的迭代算法为最大似然期望最大算法或最大后验估计算法。

4.根据权利要求1至3任意一项所述的基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于：所述步骤(2)中的CT重建模型为：

其中：f为采集到的原始投影数据；Kμ为迭代过程中生成的投影数据，K＝{K_i，j}为系统矩阵，i，j分别为待重建CT图像的像素个数和CT探测器探测元的个数；μ表示待重建的衰减图像；D_α(f‖Kμ)为数据保真项，α∈[0.6，2]；|μ|_TV为全变分正则化项。

5.根据权利要求4所述的基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于：所述步骤(3)中的迭代过程采用期望最大化迭代方法进行计算。

6.根据权利要求5所述的基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于：所述系统矩阵是从CT成像设备的系统得到的基于面积加权的系统矩阵，或者，基于体素加权的系统矩阵。

7.根据权利要求6所述的基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于：所述步骤(3)中全变分图像恢复的计算过程为： ${\mathrm{\μ}}_{k+1}^{\mathrm{\α}}={\mathrm{\μ}}_{k+.5}^{\mathrm{\α}}-\mathrm{\α}\stackrel{\‾}{\mathrm{\λ}}{\mathrm{\μ}}_k^{\mathrm{\α}}\mathrm{div}\left(w\right),$ 其中 $\stackrel{\‾}{\mathrm{\λ}}=\mathrm{\λ}/{\mathrm{\Σ}}_i{K}_{i,j};$ μ_k+1为全变分恢复后的图像数据，div为散度算子，w为二元变量。

8.根据权利要求7所述的基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于：所述步骤(4)中迭代终止条件为：相邻两次全变分图像恢复后的重建模型的差小于阈值。

9.根据权利要求8所述的基于阿尔法散度约束的全变分最小化低剂量CT重建方法，其特征在于：所述阈值为0.001。

Images