CN112656438B - 一种基于曲面全变差的低剂量ct投影域去噪及重建方法 - Google Patents

Abstract

本发明涉及一种基于曲面全变差的低剂量CT投影域去噪及重建方法，步骤一：根据投影域数据噪声的产生机制，建立优化泛函数据项Ψ；建立包括基于l_1/2范数的曲面全变差正则项的优化泛函正则项Υ，将优化泛函数据项和优化泛函正则项γ之和作为去噪优化泛函；步骤二：对步骤一中的去噪优化泛函进行求解，得到降噪后的CT投影域数据；步骤三：基于连续积分方程模型，根据步骤二中降噪后的CT投影域数据进行图像重建。本方法在进行投影域去噪的时候，建立基于l_1/2范数的曲面全变差正则项来提高去噪的精度。在完成投影域数据去噪后，采用连续积分方程模型来进行图像重建，加快了图像重建的速度和提高图像重建的精度。

Description

一种基于曲面全变差的低剂量CT投影域去噪及重建方法

技术领域

本发明涉及机CT图像处理方法领域，更具体地，涉及一种基于曲面全变差的低剂量CT投影域去噪及重建方法。

背景技术

计算机断层扫描(ComputedTomography,CT)能够以非侵入的方式得到人体内部的解剖图像，清晰地显示各部位的结构组织及病灶，帮助医生进行临床诊断，是一种广泛应用于临床检查的成像技术。然而，CT扫描过程中需要利用X射线穿透人体来采集各个角度的投影数据。许多研究表明，CT扫描过程中产生的X射线辐射会诱发人体的新陈代谢异常甚至有致癌风险，其对人体健康的影响不容忽视。为了应对标准CT扫描所产生的医源性辐射伤害问题，Naidich等人在20世纪90年代提出了低剂量CT(LowDoseCT,LDCT)。LDCT通过降低CT扫描设置中的毫安秒(mAs)设置，即降低通过X射线管的电流或缩短曝光时间来减少X射线的辐射剂量。然而，降低mAs设置会使到达探测器的光子所服从的分布由高斯分布退化为泊松分布。与标准CT相比，LDCT的原始数据生成过程中会受到相对更高水平噪声的污染。

目前在图像降噪领域被广泛采用的全变差(Total Variation，TV)正则化假设真实图像是分片常数的，如公开号为“CN104574416A”，公开日为2015年4月29日的专利公布文件中公开了一种低剂量能谱CT图像去噪方法。但CT投影域数据通常并不满足这一假设，因此直接将全变差正则化方法应用到投影域降噪问题中容易产生阶梯伪影，影响了去噪的精度。同时，在图像降噪后需要对图像进行重建，目前图像重建方法大部分是基于离散模型的算法，其所求解的图像函数解属于分片常数函数的范畴，所需要的重建时间也比较长。而且由于离散化的过程会引入误差，使用这种方法导出的优化模型天然地存在模型误差。

发明内容

本发明为克服上述现有技术中低剂量CT去噪存在阶梯伪影导致去噪精度低和重建时间长的问题，提供一种基于曲面全变差的低剂量CT投影域去噪及重建方法，实现高精度去噪和加快图像重建速度。

为解决上述技术问题，本发明采用的技术方案是：一种基于曲面全变差的低剂量CT投影域去噪及重建方法，包括如下步骤：

步骤一：根据投影域数据噪声的产生机制，建立优化泛函数据项Ψ；建立包括基于l_1/2范数的曲面全变差正则项的优化泛函正则项Υ，将优化泛函数据项和优化泛函正则项Υ之和作为优化泛函；

步骤二：对步骤一中的优化泛函进行求解，得到降噪后的CT投影域数据；

步骤三：基于连续积分方程模型，根据步骤二中降噪后的CT投影域数据进行图像重建。

在上述的技术方案中，通过建立包括正则项和数据项的优化泛函并对其进行优化求解，将输入的投影域数据中的噪声去除，实现CT投影域的去噪。其中，曲面全变差正则项引入曲面的主曲率方向，用主曲率方向代替梯度方向，并结合投影域本身的性质，主曲率方向同时包含了一阶梯度和二阶梯度信息，判断图像某点的变化强度大的方向和平滑方向比一阶梯度方向更为精确，避免图像中出现伪影，令图像去噪更加准确。在图像重建的时候采用连续积分方程模型来设计CT图像重建算法，连续积分方程基于实际物理模型，具有更高的精度和重建速度。

优选的，在所述步骤一中，优化泛函正则项Υ还包括基于l₁范数的曲波变换系数正则项。根据CT投影域数据曲波系数的稀疏性，基于l₁范数的曲波变换系数正则项能够改善去噪效果。

优选的，曲面全变差正则项建立方法为：将CT投影域数据视为

空间中的曲面，约束曲面在主曲率方向上的跳变，将主曲率方向上梯度的

范数作为曲面全变差正则项，具体表示为：

式中，Y为带噪声的输入投影域数据；B为交错差分算子；S为主曲率方向算子；

为横向差分算子；

为纵向差分算子；L是lipschitz连续系数；M为Y的行数。

具体的，定义主曲率方向算子

表示为：

式中，v₁和v₂分别为Y的最大主曲率和最小主曲率方向；D_Yρ为主曲率绝对差。

优选的，曲波变换系数正则项建立方法为对CT投影域数据Y进行离散曲波变换，得到系数CY，具体表示为：Υ₂(Y)＝||CY||₁。

优选的，步骤一中的优化泛函数据项Ψ具体为：

其中，Y为带噪声的输入投影域数据；Q为到达探测板的光子数；P为探测板所接收到的光子数；I₀为X射线入射强度。

优选的，步骤一中的优化泛函具体为：

式中，I₀为X射线入射的强度；λ₁和λ₂为两种正则化项的惩罚系数，取为正实数；C为曲波变换算子；S为主曲率方向算子；B为交错差分算子。

优选的，所述步骤三中的具体流程为：

S3.1：建立图像重建的优化泛函如下：

式中，f是图像重建函数；g为平行束投影域函数；

为Radon变换；λ是正则化项系数；Φ是恒大于零且L-lipschitz连续的函数；B表示线性变换；

S3.2：求解图像重建的优化泛函；

S3.3：对图像重建的优化泛函进行迭代求解，求得图像重建函数f。

在连续积分方程模型的基础上，建立图像重建的优化泛函中引入了Radon变换，通过Radon变换的卷积性质对后续的迭代算法进行加速，使得图像重建的速度更加快，重建时间更短。

优选的，在所述步骤S3.2中，采用邻近梯度法求解图像重建的优化泛函，具体的流程为：

首先计算

的梯度。设

则

其中

是Radon变换的共轭算子。利用

的卷积性质与傅里叶变换完成加速计算。

计算

是f在L²空间上某组基函数，z_i是基函数的系数，

是探测板的探测范围，

由此可以解得

其中利用卷积性质可以加速

的计算，设矩阵

其元素

由于W是循环矩阵，因此

w为矩阵W的第一行，z＝[z₁,z₂,…,z_L]^T。

计算

其中

表示g(θ,u)关于u分量的傅里叶变换。

优选的，在所述步骤S3.3中，采用迭代公式对图像重建的优化泛函进行迭代求解，迭代公式具体为：

其中，

表示傅里叶变换；

表示反傅里叶变换；prox_Φ是Φ的临近梯度算子；L是Φ的lipschitz连续系数；⊙表示向量按元素相乘；b^(k+1)表示第k+1步迭代的中间结果，z^{(k +1)}表示第k+1步的图像函数。

优选的，在所述步骤二中，使用带有随机步长的ADMM优化算法对步骤一中的优化泛函进行求解。采用随机步长可以加快优化算法的收敛速度。

与现有技术相比，本发明的有益效果是：本方法在进行投影与去噪的时候，建立基于l_1/2范数的曲面全变差正则项，引入曲面的主曲率方向，用主曲率方向代替梯度方向，并结合投影域本身的性质来避免图像中微影的出现，从而提高去噪的精度。在完成投影域数据去噪后，采用连续积分方程模型来进行图像重建，加快了图像重建的速度和提高图像重建的精度。

附图说明

图1是本发明的一种基于曲面全变差的低剂量CT投影域去噪及重建方法的整体流程示意图；

图2是有噪声的CT图像；

图3是高剂量的CT图像；

图4是经过本方法处理图2后的CT图像。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本专利的限制。

下面通过具体实施例，并结合附图，对本发明的技术方案作进一步的具体描述：

实施例1

如图1所示为一种基于曲面全变差的低剂量CT投影域去噪及重建方法的实施例，包括如下步骤：

步骤一：根据投影域数据噪声的产生机制，建立优化泛函数据项Ψ；建立包括基于l_1/2范数的曲面全变差正则项和于l₁范数的曲波变换系数正则项的优化泛函正则项Υ，将优化泛函数据项和优化泛函正则项Υ之和作为优化泛函；

其中，曲面全变差正则项建立方法为：将CT投影域数据视为

空间中的曲面，约束曲面在主曲率方向上的跳变，将主曲率方向上梯度的

范数作为曲面全变差正则项，具体表示为：

式中，Y为带噪声的输入投影域数据；B为交错差分算子；S为主曲率方向算子；

为横向差分算子；

为纵向差分算子；L是lipschitz连续系数；M为Y的行数。

曲波变换系数正则项建立方法为对CT投影域数据Y进行离散曲波变换，得到系数CY，具体表示为：Υ₂(Y)＝||CY||₁。

优化泛函为优化泛函数据项Ψ、曲面全变差正则项Υ₁(Y)和曲波变换系数正则项Υ₂(Y)三者之和，具体为：

其中，Y为带噪声的输入投影域数据；Q为到达探测板的光子数；P为探测板所接收到的光子数；I₀为X射线入射强度。

步骤二：使用带有随机步长的ADMM优化算法对步骤一中的优化泛函进行求解，得到降噪后的CT投影域数据；

步骤三：基于连续积分方程模型，根据步骤二中降噪后的CT投影域数据进行图像重建，具体流程为：

S3.1：建立图像重建的优化泛函如下：

式中，f是图像重建函数；g为平行束投影域函数；

为Radon变换；λ是正则化项系数；Φ是恒大于零且L-lipschitz连续的函数；A表示线性变换；

S3.2：采用邻近梯度法求解图像重建的优化泛函；

S3.3：对图像重建的优化泛函进行迭代求解，求得图像重建函数f，具体的迭代公式为：

其中，

表示傅里叶变换；

表示反傅里叶变换；prox_Φ是Φ的临近梯度算子；L是Φ的lipschitz连续系数；⊙表示向量按元素相乘；b^(k+1)表示第k+1步迭代的中间结果；z^{(k +1)}表示第k+1步的图像函数。

具体的，在所述步骤S3.2中，采用邻近梯度法求解图像重建的优化泛函，具体的流程为：

首先计算

的梯度。设

则

其中

是Radon变换的共轭算子。利用

的卷积性质与傅里叶变换完成加速计算。

计算

是f在L²空间上某组基函数，z_i是基函数的系数，

是探测板的探测范围，

由此可以解得

其中利用卷积性质可以加速

的计算，设矩阵

其元素

由于W是循环矩阵，因此

w为矩阵W的第一行，z＝[z₁,z₂,…,z_L]^T。

计算

其中

表示g(θ,u)关于u分量的傅里叶变换。

本实施例的工作原理：通过建立包括正则项和数据项的优化泛函并对其进行优化求解，将输入的投影域数据中的噪声去除，实现CT投影域的去噪。其中，曲面全变差正则项引入曲面的主曲率方向，用主曲率方向代替梯度方向，并结合投影域本身的性质，主曲率方向同时包含了一阶梯度和二阶梯度信息，判断图像某点的变化强度大的方向和平滑方向比一阶梯度方向更为精确，避免图像中出现伪影，令图像去噪更加准确。在图像重建的时候采用连续积分方程模型来设计CT图像重建算法，连续积分方程基于实际物理模型，具有更高的精度和重建速度。

本实施例的有益效果:本方法在进行投影与去噪的时候，建立基于l_1/2范数的曲面全变差正则项，引入曲面的主曲率方向，用主曲率方向代替梯度方向，并结合投影域本身的性质来避免图像中微影的出现，从而提高去噪的精度。在完成投影域数据去噪后，采用连续积分方程模型来进行图像重建，加快了图像重建的速度和提高图像重建的精度。如图4所示，经过本方法对图2低剂量有噪声的CT图像进行去噪重建，得到接近于图3的CT图像。图3的CT图像为高剂量的CT图像，是低剂量(有噪声)图像降噪重建的目标。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (8)

1.一种基于曲面全变差的低剂量CT投影域去噪及重建方法，其特征在于，包括如下步骤：

步骤一：根据投影域数据噪声的产生机制，建立优化泛函数据项Ψ；建立包括基于l_1/2范数的曲面全变差正则项的优化泛函正则项Υ，将优化泛函数据项和优化泛函正则项Υ之和作为去噪优化泛函；

建立曲面全变差正则项的方法为将CT投影域数据视为

空间中的曲面，约束曲面在主曲率方向上的跳变，将主曲率方向上梯度的

范数作为曲面全变差正则项；具体表示为：

式中，Y为带噪声的输入投影域数据；B为交错差分算子；S为主曲率方向算子；

为横向差分算子；

为纵向差分算子；M为Y的行数；

步骤二：对步骤一中的去噪优化泛函进行求解，得到降噪后的CT投影域数据；

步骤三：基于连续积分方程模型，根据步骤二中降噪后的CT投影域数据进行图像重建。

2.根据权利要求1所述的一种基于曲面全变差的低剂量CT投影域去噪及重建方法，其特征在于，在所述步骤一中，优化泛函正则项Υ还包括基于l₁范数的曲波变换系数正则项。

3.根据权利要求2所述的一种基于曲面全变差的低剂量CT投影域去噪及重建方法，其特征在于，建立曲波变换系数正则项的方法为对CT投影域数据Y进行离散曲波变换，得到系数CY，将曲波变换系数正则项表示为：Υ₂(Y)＝||CY||₁，其中C为曲波变换算子。

4.根据权利要求3所述的一种基于曲面全变差的低剂量CT投影域去噪及重建方法，其特征在于，步骤一中的优化泛函数据项Ψ具体为：

其中，Y为带噪声的输入投影域数据；Q为到达探测板的光子数；P为探测板所接收到的光子数；I₀为X射线入射强度。

5.根据权利要求4所述的一种基于曲面全变差的低剂量CT投影域去噪及重建方法，其特征在于，步骤一中的去噪优化泛函具体为：

式中，I₀为X射线入射的强度；λ₁和λ₂为两种正则化项的惩罚系数，取为正实数；C为曲波变换算子；S为主曲率方向算子；B为交错差分算子。

6.根据权利要求2所述的一种基于曲面全变差的低剂量CT投影域去噪及重建方法，其特征在于，所述步骤三中的具体流程为：

S3.1：建立图像重建的优化泛函如下：

式中，f是图像重建函数；g为平行束投影域函数；

为Radon变换；λ是正则化项系数；Φ是恒大于零且L-lipschitz连续的函数；A表示线性变换；

S3.2：求解图像重建的优化泛函；

S3.3：对图像重建的优化泛函进行迭代求解，求得图像重建函数f。

7.根据权利要求6所述的一种基于曲面全变差的低剂量CT投影域去噪及重建方法，其特征在于，在所述步骤S3.2中，采用邻近梯度法求解图像重建的优化泛函。

8.根据权利要求1-7任一所述的一种基于曲面全变差的低剂量CT投影域去噪及重建方法，其特征在于，在所述步骤二中，使用带有随机步长的ADMM优化算法对步骤一中的优化泛函进行求解。

Images