CN108492269B - 基于梯度正则卷积神经网络的低剂量ct图像去噪方法 - Google Patents

Abstract

本发明公开了一种低剂量CT图像去噪方法，主要解决现有技术中仅考虑恢复图像灰度信息而导致的图像边缘模糊、细节丢失的问题。其实现方案：1)从人体同一部位获取多张低剂量和全剂量CT图像；2)对获得的CT图像数据集进行扩充和归一化，再在每对CT图像上进行取块，得到CT图像块数据集；3)建立一个12层的全卷积去噪神经网络，使用CT图像块数据集作为网络训练数据，利用带动量项的小批量梯度下降算法优化该网络，得到训练好的网络；4)输入一张完整的低剂量CT图像到该网络中，即可输出对应的去噪CT图像。本发明能在图像去噪的同时，较好地保持图像的边缘和细节，可用于对低剂量CT图像的增强，便于医生对CT图像病例的识别。

Description

基于梯度正则卷积神经网络的低剂量CT图像去噪方法

技术领域

本发明属于图像处理技术领域，特别涉及一种CT图像的去噪方法，可用于CT设备成像时辐射量为传统CT成像的1/6到1/4的低剂量CT图像进行增强，提高图像质量和成像器官、组织的可视效果。

背景技术

传统的CT成像由于使用了较高的X射线辐射，根据成像器官或组织的不同，一般为1毫希弗到12毫希弗，对病人具有一定的辐射伤害，而且辐射可能导致的癌变和基因突变也引起了人们的担忧，因此降低成像时X射线剂量，一般为传统CT成像使用剂量的1/6到1/4，进行低剂量CT成像成为了一种新的选择。但是，这种做法会导致成像的结果质量很差，图像中会有大量的噪声和伪影，严重影响了医生对病变区域和重要器官、组织的诊断和分析。所以，急需一种有效的方法以提高低剂量CT图像的质量。

根据处理数据类型的不同，低剂量CT图像去噪主要有两种方式。其一，是对CT成像过程中成像光子的正弦信号进行预处理，滤除噪声，然后再使用反向投影算法进行成像。这种方法利用了CT成像过程中噪声的生成机制，即成像光子在统计上服从复合泊松分布，代表算法包括双边滤波、惩罚加权最小二乘法、静态小波变换和最大似然估计等等。这类算法简单高效，但对原始信号的依赖较大，空间分辨率不高。其二，是对反向投影算法重建后的图像进行处理，直接在图像域滤除噪声，对原始信号依赖较低，但无法有效地确定噪声的分布。常见方法包括小波分解、非局部均值BM3D、稀疏表示与字典学习KSVD等等。

近年来，深度学习在图像处理和分析领域掀起了一阵热浪，从低层次的图像去噪、超分辨任务到高层次的图像检测与分类任务都有涉及。这种类脑方式的信息处理机制模仿了人类视觉神经系统，对图像处理和分析十分有效，已经有一些研究工作将深度学习方法应用到了图像去噪方面，取得了不错的结果。在低剂量CT图像去噪的应用上，Kang等人首先对低剂量CT图像进行小波分解，利用得到的小波系数作为输入，然后通过卷积神经网络去除噪声，重建CT图像。Chen等人提出了一个对称的卷积和转置卷积神经网络，先对输入图像进行特征提取和编码，再经过网络对编码信号进行解码重建图像，在编解码的过程中滤除噪声。Wolterink等人率先使用了生成对抗网络进行去噪，通过生成网络来模拟低剂量CT图像的噪声生成机制，结果发现使用卷积神经网络对低剂量CT图像去噪是可行的，使用生成对抗策略可以生成更接近于全剂量的CT图像。但这类算法在重建图像的过程中，仅仅考虑了对图像灰度信息的恢复，没有过多地考虑图像的其他信息，很容易导致恢复的图像边缘模糊，细节损失严重。

发明内容

本发明的目的在于针对上述深度学习在图像重建过程中的不足，提出了一种基于梯度正则卷积神经网络的低剂量CT图像去噪方法，以保留更完整的图像边缘和细节，提高恢复图像的质量。

为实现上述目的，本发明在图像恢复的过程中，同时保持灰度信息的和梯度信息，实现恢复图像细节的保留，减小模糊边缘和损失细节的现象，其实现方案包括如下：

(1)使用CT设备同时对人体同一部位进行全剂量和低剂量成像，得到多张成对的低剂量CT图像X和全剂量CT图像Y，记为{X,Y}，其中，全剂量CT成像时X射线管电压为120千伏，管电流为200毫安，辐射剂量约为3毫希弗；低剂量CT图像时X射线管电压为120千伏，管电流为50毫安，辐射剂量为0.75毫希弗；采集的CT图像数据集记为D；

(2)对采集的CT图像数据集D进行扩充，得到扩充CT图像数据集D_A；

(3)对D_A中每张低剂量CT和全剂量CT图像进行归一化操作，将灰度值映射到区间[0,1]；

(4)对数据集D_A中每对CT图像进行取块，建立图像块数据集D_P；

(5)根据得到的图像块数据集D_P，搭建一个全卷积去噪神经网络N，网络层数通过交叉验证设为12，网络的输入为低剂量CT图像，输出为去噪后的CT图像，使用MSRA方法初始化网络的权值W，网络的偏置b均初始化为数值0；

(6)使用图像块数据集D_P，训练全卷积去噪网络N：

(6a)打乱D_P中图像块对的顺序，依次从D_P中选择256对CT图像块数据，记为一个训练批次{X_PB,Y_PB}，其中低剂量图像块X_PB作为网络的输入，全剂量图像块Y_PB用来计算网络的损失；

(6b)将X_PB输入到网络N中，经过前向传播，得到网络的输出A_PB，计算A_PB和Y_PB之间的均方误差，得到网络N的灰度损失L_V；

(6c)使用梯度算子g分别卷积A_PB和Y_PB，得到对应的梯度图G_A和G_Y，计算G_A和G_Y之间的均方误差，得到网络N的梯度损失L_G；

(6d)对网络N的灰度损失L_V和梯度损失L_G进行线性加权，得到网络的整体损失L＝L_V+λ×L_G，其中λ为加权系数；

(6e)使用带动量项的批量梯度下降算法优化整体损失L，初始学习速率为0.0001，动量项参数为0.9，权值衰减因子为0.0001，对网络的权值W和偏置b进行更新；

(6f)重复以上步骤(6a)到(6e)，直到达到网络训练的最大迭代次数T＝100000，得到训练好的去噪网络N_T；

(7)将一张完整的低剂量CT图像X_T输入到训练好的去噪网络N_T中，经过该去噪网络N_T的前向传播，得到去噪后的CT图像A_T。

本发明与现有技术相比具有以下优点：

1、本发明在图像恢复的过程中既考虑了图像灰度信息的重建，也考虑了图像梯度的保持，使得网络既有好的去噪效果，也能保留图像的大部分细节。

2、本发明是直接在图像域进行去噪，独立于具体的噪声分布和CT仪器参数，应用范围更广泛。

附图说明

图1是本发明的实现总流程图；

图2是本发明中所使用的低剂量CT图像；

图3是本发明中所使用的全剂量CT图像；

图4是本发明中所使用的低剂量CT图像扩充样例图；

图5是本发明中所使用的全剂量CT图像扩充样例图；

图6是本发明中的图像取块示例图；

图7是本发明中所构建的全卷积神经网络结构图；

图8是用本发明对低剂量CT图像去噪后的结果图。

具体实施方式

以下结合附图，对本发明的具体实施方案和效果作进一步的解释和说明：

参照图1，本发明基于梯度正则化的低剂量CT图像去噪方法，其实现步骤如下：

步骤1：数据准备。

1a)使用CT设备同时对人体同一部位进行全剂量和低剂量成像，其中，全剂量CT成像时X射线管电压为120千伏，管电流为200毫安，辐射剂量约为3毫希弗；低剂量CT图像时X射线管电压为120千伏，管电流为50毫安，辐射剂量为0.75毫希弗；

1b)将得到成对的低剂量CT图像X和全剂量CT图像Y，记为{X,Y}，其中低剂量CT图像X，如图2所示，全剂量CT图像Y如图3所示，两种图像的大小均为512*512，将采集的CT图像数据集记为D；

1c)对采集的CT图像数据集D进行扩充，得到扩充CT图像数据集D_A：

1c1)对图像对{X,Y}中的低剂量CT图像X进行扩充：

参照图4，本步骤的具体实现如下：首先，对低剂量CT图像，如图4(a)所示，沿其垂直中轴线左右翻转，得到第2张低剂量CT图像，如图4(e)所示；然后，对这2张低剂量CT图像分别沿逆时针方向旋转90度，180度和270度，得到6张新的不同角度的低剂量CT图像，如图4(b)～4(d)和4(f)～4(h)所示，即1张低剂量CT图像经扩充后共计得到8张低剂量CT图像；

1c2)对图像对{X,Y}中的全剂量CT图像Y进行扩充：

参照图5，本步骤的具体实现为：首先，对全剂量CT图像如图5(a)所示，沿其垂直中轴线左右翻转，得到第2张全剂量CT图像，如图5(e)所示；然后，对这2张全剂量CT图像分别沿逆时针方向旋转90度，180度和270度，得到6张新的不同角度的全剂量CT图像，如图5(b)～5(d)和5(f)～5(h)所示，即1张全剂量CT图像经扩充后共计得到8张全剂量CT图像；

1c3)对CT图像数据集D中每一对CT图像对{X,Y}，都进行上述步骤1c1)和1c2)的扩充操作，得到扩充CT图像数据集D_A；

1d)对D_A中每张低剂量CT和全剂量CT图像进行归一化操作，将灰度值映射到区间[0,1]，归一化公式如下：

其中，X(r,c)表示CT图像X的在第r行第c列的像素灰度值，max(X)表示图像中所有像素的最大值，min(X)表示图像中所有像素的最小值，

表示归一化后的CT图像在第r行第c列的数值；

1e)对数据集D_A中归一化后的每对CT图像进行取块，建立图像块数据集D_P；

参照图6，本步骤的具体实现如下：

1e1)将CT图像的中心区域定义为距CT图像四周像素距离至少为48的图像区域；

1e2)对图像对{X,Y}中的低剂量CT图像X进行取块：

在低剂量CT图像上的中心区域随机选择位置，记录此位置，并以此位置从左上角向右下方截取尺寸为48*48的图像块，记为X_P，在一张低剂量CT图像上共计取32个图像块；

1e3)对图像对{X,Y}中的全剂量CT图像Y进行取块：

在全剂量CT图像上，根据在低剂量CT图像上记录的取块位置，并以此位置从左上角向右下方截取尺寸为48*48的图像块Y_P，在一张全剂量CT图像上共计取32个图像块；

1e4)将截取的一一对应的图像块X_P和Y_P记为CT图像块对{X_P,Y_P}，对数据集D_A中每一对CT图像对{X,Y}都进行上述步骤1e2)和1e3)的取块操作，得到CT图像块数据集D_P。

步骤2：根据得到的图像块数据集D_P，搭建一个全卷积去噪神经网络N。

参照图7，本步骤的具体实现如下：

2a)将CT图像作为网络N的输入，网络N的输出结果为去噪后的CT图像，输入输出图像大小保持不变；

2b)通过交叉验证，将网络N的层数设为12，包括输入层、中间层和输出层；

2c)设第1层为输入层，其包含一个1个卷积层Conv和1个修正线性单元激活层ReLU，该卷积层含有64个卷积核，每个卷积核大小为3*3，如图7中黄色模块所示；

2d)设第2层至第11层为中间层，该中间层的每一层都由同一种模块构造，该模块按连接顺序依次含有1个卷积层、1个ReLU层和1个批规范化层BN，该卷积层含有64个卷积核，每个卷积核大小为3*3，如图7中绿色模块所示；

2e)设第12层为输出层，其包含1个卷积层，该卷积层只包含1个卷积核，卷积核大小为3*3，如图7中蓝色模块所示。

上述网络N中的卷积层，其数学形式如下：

其中，

表示网络第l层的第i个特征图，当l＝0时，F⁰表示网络输入的低剂量CT图像块，

表示网络第l层的第i个卷积核的权值，

表示网络第l层的第i个卷积核的偏置，n_l表示网络第l层的卷积核数量，Σ表示求和运算，

表示图像卷积运算，该卷积运算采用“same”方式，以保持卷积前后图像大小不变。

上述网络N中的修正线性单元激活层ReLU，其数学形式如下：

其中，x表示输入数据；

上述网络N中的批规范化层BN，其数学形式如下：

其中，x表示输入数据，

表示其均值，σ表示其标准差，γ表示BN层的斜率参数，β表示BN层的偏置参数。

步骤3，对网络N的权值W和偏置b进行初始化。

3a)使用MSRA方法对网络N的权值W进行初始化，其公式如下：

其中，W表示网络的权值，N(·,·)表示高斯分布，即网络的权值W服从均值为0，标准差为

的高斯分布；

3b)将网络N的偏置b均初始化为数值0。

步骤4：使用图像块数据集D_P，训练全卷积去噪网络N。

4a)打乱D_P中图像块对的顺序，依次从D_P中选择m对CT图像块数据，记为一个训练批次{X_PB,Y_PB}，其中低剂量图像块

作为网络的输入，全剂量图像块

用来计算网络的损失，i＝1,2,...,m，m表示批量大小，本实例中取值为m＝256；

4b)将X_PB输入到网络N中，经过前向传播，得到网络的输出A_PB，计算A_PB和Y_PB之间的均方误差，得到网络N的灰度损失L_V：

其中，row和col分别表示输入的图像块的高度和宽度，取值均为48，

表示第i个全剂量CT图像块的第r行第c列的像素灰度值，

表示第i个网络输出去噪图像块的第r行第c列的像素灰度值，||·||²表示平方运算，∑表示求和运算；

4c)选择梯度卷积算子g，本实例中选取如下：

4d)使用梯度算子g分别卷积A_PB和Y_PB，得到对应的梯度图G_A和G_Y：

其中，

表示图像卷积运算，该卷积运算采用“same”方式，以保持卷积前后图像大小不变。

4e)计算G_A和G_Y之间的均方误差，得到网络N的梯度损失L_G：

其中，

表示第i个全剂量CT图像块对应的梯度图在第r行第c列的数值，

表示第i个网络输出去噪图像块对应的梯度图在第r行第c列的数值；

4f)对网络N的灰度损失L_V和梯度损失L_G进行线性加权，得到网络的总体损失L：

L＝L_V+λ×L_G，

其中，λ为加权系数，通过交叉验证取值为2；

4g)使用带动量项的批量梯度下降算法优化整体损失L，初始学习速率为0.0001，动量项参数为0.9，权值衰减因子为0.0001，对网络的权值W进行更新；

4h)重复以上步骤(4a)到(4g)，直到达到网络训练的最大迭代次数T＝100000，得到训练好的去噪网络N_T。

步骤5：使用训练好的去噪网络N_T对低剂量CT图像进行去噪。

将一张完整的低剂量CT图像X_T，如图8(a)所示，输入到训练好的去噪网络N_T中，经过该去噪网络N_T的前向传播，得到去噪后的CT图像A_T，如图8(b)所示。

去噪网络虽然是利用图像块去训练的，但其是一个全卷积网络，不受输入图像尺寸的影响，因此可以直接输入完整的CT图像进行去噪。使用图像块进行训练是因为完整的CT图像数量有限，可能导致训练的网络过拟合，且使用图像块训练可以加快网络的训练速度。

从图8可以看出，在图像恢复的过程中同时考虑图像灰度信息的重建和图像梯度的保持，得到的去噪的后图像细节更丰富、清晰，可视性更好。

以上描述仅是本发明的一个具体实例，不构成对本发明的任何限制，显然对于本领域的专业人员来说，在了解了本发明内容和原理后，都可能在不背离本发明原理、结构的情况下，进行形式和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。

Claims (10)

1.一种基于梯度正则卷积神经网络的低剂量CT图像去噪方法，包括：

(1)使用CT设备同时对人体同一部位进行全剂量和低剂量成像，得到多张成对的低剂量CT图像X和全剂量CT图像Y，记为{X,Y}，其中，全剂量CT成像时X射线管电压为120千伏，管电流为200毫安，辐射剂量约为3毫希弗；低剂量CT成像时X射线管电压为120千伏，管电流为50毫安，辐射剂量为0.75毫希弗；采集的CT图像数据集记为D；

(2)对采集的CT图像数据集D进行扩充，得到扩充CT图像数据集D_A；

(3)对D_A中每张低剂量CT和全剂量CT图像进行归一化操作，将灰度值映射到区间[0,1]；

(4)对数据集D_A中每对CT图像进行取块，建立图像块数据集D_P；

(5)根据得到的图像块数据集D_P，搭建一个全卷积去噪神经网络N，网络层数通过交叉验证设为12，网络的输入为低剂量CT图像，输出为去噪后的CT图像，使用MSRA方法初始化网络的权值W，网络的偏置b均初始化为数值0；

(6)使用图像块数据集D_P，训练全卷积去噪神经网络N：

(6a)打乱D_P中图像块对的顺序，依次从D_P中选择256对CT图像块数据，记为一个训练批次{X_PB,Y_PB}，其中低剂量图像块X_PB作为网络的输入，全剂量图像块Y_PB用来计算网络的损失；

(6b)将X_PB输入到网络N中，经过前向传播，得到网络的输出A_PB，计算A_PB和Y_PB之间的均方误差，得到网络N的灰度损失L_V；

(6c)使用梯度算子g分别卷积A_PB和Y_PB，得到对应的梯度图G_A和G_Y，计算G_A和G_Y之间的均方误差，得到网络N的梯度损失L_G；

(6d)对网络N的灰度损失L_V和梯度损失L_G进行线性加权，得到网络的整体损失L＝L_V+λ×L_G，其中λ为加权系数；

(6e)使用带动量项的批量梯度下降算法优化整体损失L，初始学习速率为0.0001，动量项参数为0.9，权值衰减因子为0.0001，对网络的权值W和偏置b进行更新；

(6f)重复以上步骤(6a)到(6e)，直到达到网络训练的最大迭代次数T＝100000，得到训练好的去噪网络N_T；

(7)将一张完整的低剂量CT图像X_T输入到训练好的去噪网络N_T中，经过该去噪网络N_T的前向传播，得到去噪后的CT图像A_T。

2.根据权利要求1所述的方法，其中步骤(2)对采集的CT图像数据集D进行扩充，得到扩充CT图像数据集D_A，按如下步骤进行：

(2a)对图像对中低剂量CT图像X进行扩充：

首先，对低剂量CT图像沿其垂直中轴线左右翻转，得到第2张低剂量CT图像，然后对这2张低剂量CT图像分别逆时针旋转90度，180度和270度，得到6张新的低剂量CT图像，即对1张低剂量CT图像经扩充后共计有8张低剂量CT图像；

(2b)对图像对中全剂量CT图像Y进行扩充：

首先，对全剂量CT图像沿其垂直中轴线左右翻转，得到第2张全剂量CT图像，然后对这2张全剂量CT图像分别逆时针旋转90度，180度和270度，得到6张新的全剂量CT图像，即对1张全剂量CT图像经扩充后共计有8张全剂量CT图像；

(2c)重复操作(2a)和(2b)，依次处理CT图像数据集D中每一对CT图像，得到扩充CT图像数据集D_A。

3.根据权利要求1所述的方法，其中步骤(3)对每张低剂量CT图像和全剂量CT图像归一化到[0,1]，按如下公式进行：

其中X为输入的CT图像，min(X)表示CT图像灰度值的最小值，max(X)表示CT图像灰度值的最大值，

表示归一化后的图像。

4.根据权利要求1所述的方法，其中步骤(4)对数据集D_A中每对CT图像进行取块，建立图像块数据集D_P，按如下步骤进行：

(4a)对图像对中低剂量CT图像X进行取块：

在低剂量CT图像上的中心区域随机选择位置，记录此位置，并以此位置从左上角向右下方截取尺寸为48*48的图像块X_P，共计取32个；

(4b)对图像对中全剂量CT图像Y进行取块：

在全剂量CT图像上，根据在低剂量CT图像上记录的取块位置，并以此位置从左上角向右下方截取尺寸为48*48的图像块Y_P，共计取32个；

(4c)截取的图像块X_P和Y_P，记为CT图像块对{X_P,Y_P}，重复操作(4a)和(4b)，依次处理D_A中每一对CT图像，得到CT图像块数据集D_P。

5.根据权利要求1所述的方法，其中步骤(5)中根据图像块数据集D_P，搭建的12层全卷积去噪网络N，包括输入层、中间层和输出层，其中：

第1层是输入层，其包含一个1个卷积层和1个修正线性单元激活层，该卷积层含有64个卷积核，每个卷积核大小为3*3；

第2层至第11层是中间层，每一层都由同一种模块构造，该模块按连接顺序依次含有1个卷积层、1个修正线性单元激活层和1个批规范化层，该卷积层含有64个卷积核，每个卷积核大小为3*3；

第12层是输出层，其包含1个卷积层，该卷积层只包含1个卷积核，卷积核大小为3*3。

6.根据权利要求1所述的方法，其中步骤(5)对12层全卷积去噪网络N使用MSRA方法初始化权值，按如下公式进行：

其中，W表示网络N的权值，N(·,·)表示高斯分布，即网络N的权值W服从均值为0，标准差为

的高斯分布。

7.根据权利要求1所述的方法，其中步骤(6b)中计算网络N的灰度损失L_V，按如下公式进行：

其中，m为训练批次大小，取值为256，||·||²表示平方运算，Σ表示求和运算，Y_PB表示全剂量CT图像块，A_PB表示网络N输出的去噪图像块。

8.根据权利要求1所述的方法，其中步骤(6c)中的梯度算子g，全剂量CT图像块的梯度图G_Y和网络输出图像块的梯度图G_A，按如下公式确定：

其中，Y_PB表示全剂量CT图像块，A_PB表示网络输出的去噪图像块，

表示图像卷积运算，该卷积运算采用“same”方式，以保持卷积前后图像大小不变。

9.根据权利要求1所述的方法，其中步骤(6c)中计算网络N的梯度损失L_G，按如下公式进行：

其中，m为训练批次大小，取值为256，||·||²表示平方运算，Σ表示求和运算，G_Y表示全剂量CT图像块对应的梯度图，G_A表示网络N输出的图像块对应的梯度图。

10.根据权利要求1所述的方法，其中步骤(6d)中计算网络的整体损失L，按如下公式进行：

L＝L_V+λ×L_G，

其中，L_V为灰度损失，L_G为梯度损失，λ为加权系数，取值为2。

Images