CN104899907B - 基于伽马先验的稀疏角度ct图像重建方法 - Google Patents

Abstract

本发明公开了一种基于伽马先验的稀疏角度CT图像重建方法，能够直接应用于医用CT投影数据，解决了在稀疏角度CT成像中，易受到条状伪影的影响，CT图像质量发生显著下降的问题。本发明方法中，首先按照等角度对物体进行稀疏角度采样，获取稀疏投影数据；然后对所采集的投影数据利用伽马先验模型进行迭代重建；最后对重建后的CT图像按照合适的窗宽和窗位进行显示，临床数据的实验结果证明了本发明方法能有效的抑制稀疏角度下CT图像中的条状伪影，使稀疏角度下CT的图像质量接近正常角度CT图像的临床质量水平，本发明可以在目前的CT设备中获得很好的应用和推广。

Description

基于伽马先验的稀疏角度CT图像重建方法

技术领域

本发明涉及一种基于伽马先验的稀疏角度CT图像重建方法，属于计算机断层显像技术领域内的稀疏角度CT成像问题。

背景技术

作为目前一种常规有效的临床医学诊断工具，X射线计算机断层成像(X-rayComputerized Tomography，CT)能够在毫米尺度上清晰的获得人体不同组织对于X射线的衰减信息，从而为临床医生的诊断和预防提供丰富的三维人体器官组织信息。今天，作为一种成熟的和在临床上普遍认可的检查方法，CT已经成为放射诊断领域内不可缺少的主要工具之一。然而，随着CT断层扫描在临床诊断尤其是常规检查中的普及，CT扫描中的辐射剂量问题已经引起了人们越来越多的关注，大量的临床研究表明超过正常范围的CT辐射剂量易诱发人体新陈代谢异常乃至癌症等疾病。

当前降低CT剂量的算法主要分为降低发射管电压(或电流)和稀疏角度扫描两大类，前者主要是基于降低每次发射的光子数的，这样可以减少每次的剂量，但是不可避免的带来噪声。例如使用当前在临床CT重建中普遍使用的滤波反投影(FilteredBackprojection，FBP)算法进行处理，往往带来非常严重的噪声。后者基于稀疏角度的重建，一般并不降低正常所需的电压和电流，但是在每次扫描时，尽量减少扫描角度，从而达到降低总的剂量。然而，在目前的CT设备中，通过减少扫描角度容易产生条状伪影，降低CT重建的质量，从而影响临床医生对异常组织的确诊率。

随着压缩感知(compressed sensing(CS))的兴起，基于全变分(total variation(TV))的重建算法取得的成功的进展。此方法能够利用CT图像域的局部梯度稀疏信息先验来获得更好的重建效果。其原理在E.Candes，J.Romberg等的“Robust uncertaintyprinciples:Exact signal reconstruction from highly incomplete frequencyinformation”和D.Donoho的“Compressed sensing”等的压缩感知的参考文献中进行了详细的分析。基于TV准则的稀疏角度重建在角度降低的一定程度后，仍然会出现明显的条状伪影。

发明内容

发明目的：针对现有技术中存在的问题，本发明提出一种基于伽马先验的稀疏角度CT重建方法，通过对传统的迭代方法中TV先验进行改进，抑制稀疏角度下CT图像内的条状伪影，提高稀疏角度CT图像质量，使其达到符合临床医生诊断的质量水平。

技术方案：为实现上述发明目的，本发明采用如下技术方案：

一种基于伽马先验的稀疏角度CT重建方法，包括如下步骤：

(1)对目标物进行等角度稀疏角度扫描获取投影数据；

(2)基于伽马先验模型对所获取的投影数据进行迭代重建，得到重建后的图像；所述伽马先验模型为其中f表示重建后的图像，G表示成像系统矩阵，y表示投影数据，λ表示保真项与约束项Ψ_Γ(f)的平衡参数，Ψ_Γ(f)表示以伽马稀疏项作为约束的函数；

(3)显示所述重建后的图像。

优选地，所述Ψ_Γ(f)定义为：

其中，α表示伽马函数的形状参数，β表示伽马函数的尺度参数，m,n分别表示重建后图像的两个维数，Γ(α)表示伽马函数在α处的值，表示在像素f_i,j处水平方向梯度的模，表示在像素f_i,j处竖直方向梯度的模。

优选地，所述步骤(2)的迭代重建利用共轭梯度法进行迭代更新，具体步骤包括：

(2.1)设置迭代参数，所述参数包括平衡参数、伽马形状参数、伽马尺度参数、步长选择参数和最大迭代次数；

(2.2)初始化迭代变量，所变量包括当前迭代次数、迭代初始图像、初始梯度方向和初始搜索方向；

(2.3)按照共轭梯度法进行迭代更新，直到梯度方向满足设定的最小幅度阈值条件或超过最大迭代次数。

优选地，所述步骤(2.3)中采用Armijo回溯搜索方法进行迭代步长的更新。

有益效果：本发明通过对压缩感知理论的分析，得出其原因在于先验项不够“稀疏”，也就是不够接近压缩感知理论中的提到了L0范数。基于此，借助概率统计的知识，提出了伽马先验。伽马先验能够较TV先验更加逼近L0范数。实验结果验证了在投影角度降低到一定角度后，伽马先验能得到较TV先验更佳图像，即能够有效的抑制CT图像中在稀疏条件下易出现的条状伪影，使稀疏角度下CT的图像质量接近正常角度CT图像的临床质量水平，本发明可以在目前的CT设备中获得很好的应用和推广。

附图说明

图1为本发明实施例的方法流程图。

图2为本发明对比实验中临床骨盆投影数据图。

图3为本发明对比实验中全角度重建的结果图。

图4为本发明对比实验中TV方法重建的结果图。

图5为本发明对比实验中伽马方法重建的结果图。

图6为本发明对比实验中重建图像的局部放大对比图(a、b、c分别对应全角度、TV和伽马三种重建方法)。

图7为本发明对比实验中重建图像的局部放大对比图(a、b、c分别对应全角度、TV和伽马三种重建方法)。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

定义f表示所要重建的图像，TV方法的模型如下：

本发明的基于伽马先验的稀疏角度CT重建方法，对TV方法进行了改进，建立基于伽马先验迭代模型，该模型以把稀疏角度投影数据重建为尽可能接近与之相对应的全角度投影数据重建的图像为目标，以伽马稀疏项作为约束，模型可表示为：

其中，G表示成像系统矩阵，根据成像系统自动生成的，B.D.Man和S.Basu，的“Distance-driven projection and backprojection in three dimensions”文献中有关于系统矩阵的详细描述，f表示所要重建的图像，Φ(f)表示所建立的最小化能量函数模型，λ表示保真项与约束项Ψ_Γ(f)的平衡参数。Γ(α)表示伽马函数在α处的值。表示在像素f_i，j处水平方向梯度的模。表示在像素f_i，j处竖直方向梯度的模。α表示伽马函数的形状参数，β表示伽马函数的尺度参数，m，n分别表示重建图像的两个维数。

如图1所示，本发明实施例公开的一种基于伽马先验的稀疏角度CT图像重建方法，主要包括：首先，建立系统矩阵G，对目标物进行等角度稀疏角度扫描获取投影数据y；然后，调整伽马先验向里面的形状参数α和尺度参数β，基于伽马先验模型对所获取的投影数据进行迭代重建，迭代过程中的步长选取方法使用的是Armijo回溯方法；最后，对所得到的重建后的图像选择合适的窗宽和窗位进行显示。其中具体迭代重建的步骤为：

(1)设置迭代参数的值，包括平衡参数λ，伽马形状参数α，伽马尺度参数β，以及步长选择参数η，ρ，τ，最大迭代次数N₀。

(2)初始化迭代变量，包括当前迭代次数k＝1，迭代初始图像f₀，初始梯度方向g₀，和初始搜索方向d₀。

(3)按照共轭梯度法进行搜索步长迭代，直到迭代停止条件。

上述步骤的算法流程的伪代码如下：

上述算法中参数具体指代内容，列在下表中。

表1

效果评估：我们对获得的一个临床骨盆投影数据(图2)进行了等间隔采样。此次扫描的具体参数是电压120kVp，电流100mA，306个投影角度，且探测器的个数地512。图3为全角度重建的结果。图4和图5分别TV方法和本发明的伽马方法等间隔90个采样角度的重建结果。

通过观察图3至图5(分别对应的是全角度重建图像，TV先验稀疏重建的CT图像，伽马先验稀疏重建的图像)以及局部放大图(图6、图7)，可以看到伽马先验方法较经典的TV先验可以有效抑制条状伪影。

为了量化的验证本发明算法的有效性，我们计算了两个选定区域(图3中实线框及虚线框所示区域)相对于全角度重建图像(图3)的PSNR，这里PSNR的定义为：

式中，I表示重建图像，P表示参考图像。m和n表示图像大小。P_max表示像素灰度的最大值。全图以及局部图的PSNR值列在下表2中。结果表明所提出的方法，较经典的TV约束无论是数量还是视觉质量上进行对比，都具有较好的可比性。

表2

方法	全图像	实线框区域	虚线框区域
				TV	39.84	38.47	35.63
伽马	40.80	39.11	37.31

Claims (3)

1.一种基于伽马先验的稀疏角度CT图像重建方法，其特征在于，包括如下步骤：

(1)对目标物进行等角度稀疏角度扫描获取投影数据；

(2)基于伽马先验模型对所获取的投影数据进行迭代重建，得到重建后的图像；所述伽马先验模型为其中f表示重建后的图像，G表示成像系统矩阵，y表示投影数据，λ表示保真项与约束项Ψ_Γ(f)的平衡参数，Ψ_Γ(f)表示以伽马稀疏项作为约束的函数；Ψ_Γ(f)定义为：

${\mathrm{\Ψ}}_{\mathrm{\Γ}}\left(f\right)={\mathrm{\Σ}}_{l=1}^m{\mathrm{\Σ}}_{j=1}^n{\∫}_0^{|\▿f_{i,j}^v|}\frac{t^{\mathrm{\α}-1}{\mathrm{\β}}^{\mathrm{\α}}{e}^{-\mathrm{\β}t}}{\mathrm{\Γ}\left(\mathrm{\α}\right)}dt+{\mathrm{\Σ}}_{i=1}^m{\mathrm{\Σ}}_{j=1}^n{\∫}_0^{|\▿f_{i,j}^h|}\frac{t^{\mathrm{\α}-1}{\mathrm{\β}}^{\mathrm{\α}}{e}^{-\mathrm{\β}t}}{\mathrm{\Γ}\left(\mathrm{\α}\right)}dt$

其中，α表示伽马函数的形状参数，β表示伽马函数的尺度参数，m,n分别表示重建后图像的两个维数，Γ(α)表示伽马函数在α处的值，表示在像素f_i,j处水平方向梯度的模，表示在像素f_i,j处竖直方向梯度的模；

(3)显示所述重建后的图像。

2.根据权利要求1所述的基于伽马先验的稀疏角度CT图像重建方法，其特征在于，所述步骤(2)的迭代重建利用共轭梯度法进行迭代更新，具体步骤包括：

(2.1)设置迭代参数，所述参数包括平衡参数、伽马形状参数、伽马尺度参数、步长选择参数和最大迭代次数；

(2.2)初始化迭代变量，所变量包括当前迭代次数、迭代初始图像、初始梯度方向和初始搜索方向；

(2.3)按照共轭梯度法进行迭代更新，直到梯度方向满足设定的最小幅度阈值条件或超过最大迭代次数。

3.根据权利要求2所述的基于伽马先验的稀疏角度CT图像重建方法，其特征在于，所述步骤(2.3)中采用Armijo回溯搜索方法进行迭代步长的更新。