使用深度学习(卷积自动编码器)重建大脑MRI图像 _卷积

本文概述

了解脑MRI 3T数据集
训练模型
保存模型
在嘈杂的3T图像上预测
定量指标：峰值信噪比(PSNR)

你将使用3T脑MRI数据集来训练你的网络。为了观察模型的有效性, 你将在以下位置测试模型：

看不见的3T MRI图像,
嘈杂的3T MRI图像和
使用定性度量：峰值信噪比(PSNR)评估重建图像的性能。
本教程将不会解决医学成像的复杂性, 而是将重点放在深度学习方面！注意：本教程将主要介绍卷积自动编码器的实际实现。因此, 如果你尚不了解卷积神经网络(CNN)和自动编码器, 则可能需要查看CNN和自动编码器教程。关于本教程的最佳部分是, 你将以2D图像的形式加载3D体积并将其馈送到模型中。简而言之, 你将在今天的教程中解决以下主题：
一开始, 你将了解磁共振成像(MRI),
然后, 你将了解脑MRI数据集：它具有什么样的图像, 导入模块, 如何读取图像, 创建图像阵列, 对脑MRI图像进行预处理以能够将其输入模型中, 最后探索大脑MRI图像。
在卷积自动编码器的实现中：将预处理的数据拟合到模型中, 可视化训练和验证损失图, 通过训练后的模型, 最后在测试集上进行预测。
接下来, 你将通过在测试图像中添加噪声来测试预训练模型的鲁棒性, 并查看模型的定量性能。
最后, 你将使用定量度量峰值信噪比(PSNR)来测试你的预测, 并测量模型的性能。
在医学成像中使用了多种系统, 从磁场强度为0.3特斯拉(T)的开放式MRI装置到场强高达1.0 T的四肢MRI系统以及场强高达3.0 T的全身扫描仪(临床用)。 Tesla是测量MR图像磁场定量强度的单位。高场MR扫描仪(7T, 11.5T)即使具有较小的体素(3维贴片或网格)尺寸, 也可以产生更高的SNR(信噪比), 因此是更准确诊断的首选。
较小的体素尺寸可导致更好的分辨率, 进而有助于临床诊断。但是, MR扫描仪中使用的磁场强度对体素大小设置了下限, 以保持良好的信噪比(SNR), 以便保留MR图像细节。
尽管7T和11.5T具有出色的图像质量, 但由于成本限制, 它们很少在生产中部署。
根据最近的论文, 报告的3T扫描仪数量约为20, 000, 而只有40台7T扫描仪。
了解脑MRI 3T数据集大脑MRI数据集由3D体积组成, 每个体积在大脑的不同切片上拍摄的MRI图像的总数为207个切片/图像。每个切片的尺寸为173 x173。图像是单通道灰度图像。共有30个对象, 每个对象都包含患者的MRI扫描。图像格式不是jpeg, png等, 而是nifti格式。你将在下一节中看到如何读取nifti格式的图像。
该数据集由T1模态MR图像组成, 传统上认为T1序列可用于评估解剖结构。你今天将要使用的数据集包括3T脑MRI。
该数据集是公开的, 可以从此源下载。
提示：如果你想学习如何使用MNIST数据集为分类任务实现多层感知器(MLP), 请查看本教程。
注意：开始之前, 请注意, 该模型将在配备32GB RAM的Nvidia 1080 Ti GPU Xeon e5 GeForce处理器的系统上进行训练。如果你正在使用Jupyter Notebook, 则将需要再添加三行代码, 在其中使用称为os的模块指定CUDA设备顺序和CUDA可见设备。
在下面的代码中, 你基本上使用os.environ在笔记本中设置了环境变量。在初始化Keras以限制Keras后端TensorFlow使用第一个GPU之前, 最好执行以下操作。如果你要在其上训练的机器的GPU为0, 请确保使用0而不是1。你可以通过在终端上运行一个简单的命令来进行检查：例如nvidia-smi

import osos.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"os.environ["CUDA_VISIBLE_DEVICES"]="1" #model will be trained on GPU 1

首先, 导入所有必需的模块, 例如cv2, numpy, matplotlib和最重要的keras, 因为你将在今天的教程中使用该框架！
为了读取nifti格式的图像, 你还必须导入一个名为nibabel的模块。

import osimport cv2from keras.layers import Input, Dense, Flatten, Dropout, merge, Reshape, Conv2D, MaxPooling2D, UpSampling2D, Conv2DTransposefrom keras.layers.normalization import BatchNormalizationfrom keras.models import Model, Sequentialfrom keras.callbacks import ModelCheckpointfrom keras.optimizers import Adadelta, RMSprop, SGD, Adamfrom keras import regularizersfrom keras import backend as K

Using TensorFlow backend.

import numpy as npimport scipy.miscimport numpy.random as rngfrom PIL import Image, ImageDraw, ImageFontfrom sklearn.utils import shuffleimport nibabel as nib #reading MR imagesfrom sklearn.cross_validation import train_test_splitimport mathimport globfrom matplotlib import pyplot as plt%matplotlib inline

/usr/local/lib/python3.5/dist-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20."This module will be removed in 0.20.", DeprecationWarning)

你将使用glob模块, 该模块将返回一个列表, 其中包含你指定的文件夹中的所有卷！

ff = glob.glob('ground3T/*')

让我们打印列表的第一个元素, 并检查列表的长度：在我们的例子中应该是30。

ff[0]

'ground3T/181232.nii.gz'

len(ff)

30

现在, 你都已准备好使用nibabel加载3D卷。请注意, 当你加载Nifti格式卷时, Nibabel不会加载图像阵列。它一直等到你要求数组数据。要求数组数据的通常方法是调用get_data()方法。
由于你要使用2D切片而不是3D, 因此将初始化一个列表, 其中：每次阅读卷时, 你都将遍历3D卷的全部207个切片, 并将每个切片一个接一个地追加到列表中。

images = []

我们还要打印3D体积之一的形状, 它的形状应为173 x 207 x 173(x, y, z；坐标)。
注意：你将仅使用大脑的中间51个切片, 而不是全部207个切片。因此, 我们还看看如何仅使用中心切片并加载它们。

for f in range(len(ff)):a = nib.load(ff[f])a = a.get_data()a = a[:, 78:129, :]for i in range(a.shape[1]):images.append((a[:, i, :]))print (a.shape)

(173, 51, 173)

让我们分析207个切片中一个切片的形状。

a[:, 0, :].shape

(173, 173)

由于图像是列表, 因此你将使用numpy模块将列表转换为numpy数组。

images = np.asarray(images)

现在该检查numpy数组的形状了, 该数组的第一个尺寸应分别为207 x 30 = 6210, 其余两个尺寸应为173 x 173。

images.shape

(1530, 173, 173)

数据集的图像确实是尺寸为173 x 173的灰度图像, 因此在将数据输入模型之前, 对其进行预处理非常重要。首先, 你将每个173 x 173图像转换为173 x 173 x 1的矩阵, 你可以将其输入网络：

images = images.reshape(-1, 173, 173, 1)

images.shape

(1530, 173, 173, 1)

接下来, 使用最大-最小归一化技术重新缩放数据：

m = np.max(images)mi = np.min(images)

m, mi

(3599.0959, -341.83853)

images = (images - mi) / (m - mi)

重新缩放后, 让我们验证数据的最小值和最大值, 分别为0.0和1.0！

np.min(images), np.max(images)

(0.0, 1.0)

这是重要的一步, 这里你将在边界处用零填充图像, 以使图像的尺寸均匀, 并且在将图像传递通过模型时, 更容易对图像进行二分下采样。让我们在三行三列中添加零, 使尺寸为176 x 176

temp = np.zeros([1530, 176, 176, 1])

temp[:, 3:, 3:, :] = images

images = temp

完成所有这些之后, 对数据进行分区很重要。为了使模型更好地泛化, 你将数据分为两部分：训练和验证集。你将在80％的数据上训练模型, 并在剩余训练数据的20％上验证模型。
这也将帮助你减少过度拟合的机会, 因为你将根据训练阶段无法看到的数据来验证模型。
你可以使用scikit-learn的train_test_split模块正确地划分数据：

from sklearn.model_selection import train_test_splittrain_X, valid_X, train_ground, valid_ground = train_test_split(images, images, test_size=0.2, random_state=13)

请注意, 对于此任务, 你不需要培训和测试标签。这就是为什么你将两次通过训练图像。与分类任务中的标签类似, 你的训练图像将既是输入也是基础事实。
现在, 让我们分析数据集中的图像外观, 并再次查看图像的尺寸, 因为借助NumPy array属性.shape在图像上添加了三行和三列：

# Shapes of training setprint("Dataset (images) shape: {shape}".format(shape=images.shape))

Dataset (images) shape: (1530, 176, 176, 1)

从上面的输出中, 你可以看到数据的形状为6210 x 176 x 176, 这是因为存在176 x 176 x 1维矩阵的6210个样本。
现在, 让我们看一下数据集中的几个训练和验证图像：

plt.figure(figsize=[5, 5])# Display the first image in training dataplt.subplot(121)curr_img = np.reshape(train_X[0], (176, 176))plt.imshow(curr_img, cmap='gray')# Display the first image in testing dataplt.subplot(122)curr_img = np.reshape(valid_X[0], (176, 176))plt.imshow(curr_img, cmap='gray')

< matplotlib.image.AxesImage at 0x7ff67d5d59b0>

文章图片
以上两个图的输出来自训练和验证集。你会看到训练图像和验证图像都不同, 看看卷积自动编码器是否能够学习特征并正确重建这些图像将很有趣。
现在你已经准备好定义网络并将数据馈入网络。因此, 事不宜迟, 让我们跳到下一步！
图像的尺寸为176 x 176 x 1或30976维矢量。你将图像矩阵转换为数组, 在0到1之间调整比例, 调整形状以使其尺寸为176 x 176 x 1, 并将其作为输入提供给网络。
同样, 你将使用128的批次大小, 最好使用256或512的较高批次大小, 这完全取决于你训练模型的系统。它在确定学习参数方面做出了巨大贡献, 并影响了预测准确性。你将训练你的网络50个纪元。

batch_size = 128epochs = 300inChannel = 1x, y = 176, 176input_img = Input(shape = (x, y, inChannel))

你可能已经知道, 自动编码器分为两个部分：有一个编码器和一个解码器。
编码器：它具有3个卷积块, 每个块都有一个卷积层, 后跟一个批处理归一化层。最大卷积层用于第一和第二卷积块之后。

第一个卷积块将包含32个大小为3 x 3的滤波器, 然后是下采样(最大合并)层,
第二个块将包含64个大小为3 x 3的滤波器, 然后是另一个下采样层,
编码器的最后一块将具有128个大小为3 x 3的滤波器。

解码器：它具有2个卷积块, 每个块都有一个卷积层, 后跟一个批处理归一化层。在第一和第二卷积块之后使用上采样层。

第一块将具有128个大小为3 x 3的滤波器, 其后是一个上采样层,
第二块将具有64个尺寸为3 x 3的滤波器, 然后是另一个上采样层,
编码器的最后一层将具有1个大小为3 x 3的滤波器, 它将重建具有单个通道的输入。

每次使用时, 最大池化层将对输入进行两次下采样, 而每次使用时, 向上采样层将对输入进行两次下采样。
注意：滤波器的数量, 滤波器的大小, 层数, 训练模型的时期数都是超参数, 应根据自己的直觉来决定, 你可以通过调整这些超参数和衡量模型的性能。这就是你将如何慢慢学习深度学习的艺术！

def autoencoder(input_img):#encoder#input = 28 x 28 x 1 (wide and thin)conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(input_img) #28 x 28 x 32conv1 = BatchNormalization()(conv1)conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv1)conv1 = BatchNormalization()(conv1)pool1 = MaxPooling2D(pool_size=(2, 2))(conv1) #14 x 14 x 32conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1) #14 x 14 x 64conv2 = BatchNormalization()(conv2)conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)conv2 = BatchNormalization()(conv2)pool2 = MaxPooling2D(pool_size=(2, 2))(conv2) #7 x 7 x 64conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2) #7 x 7 x 128 (small and thick)conv3 = BatchNormalization()(conv3)conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)conv3 = BatchNormalization()(conv3)#decoderconv4 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv3) #7 x 7 x 128conv4 = BatchNormalization()(conv4)conv4 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv4)conv4 = BatchNormalization()(conv4)up1 = UpSampling2D((2, 2))(conv4) # 14 x 14 x 128conv5 = Conv2D(32, (3, 3), activation='relu', padding='same')(up1) # 14 x 14 x 64conv5 = BatchNormalization()(conv5)conv5 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv5)conv5 = BatchNormalization()(conv5)up2 = UpSampling2D((2, 2))(conv5) # 28 x 28 x 64decoded = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(up2) # 28 x 28 x 1return decoded

创建模型后, 必须使用优化器将其编译为RMSProp。
注意, 你还必须通过参数loss指定损失类型。在这种情况下, 这就是均方误差, 因为将使用逐像素均方误差来计算每批预测输出和地面真实值之间的每批损失：

autoencoder = Model(input_img, autoencoder(input_img))autoencoder.compile(loss='mean_squared_error', optimizer = RMSprop())

让我们使用摘要功能来可视化在上一步中创建的图层, 这将显示每个图层中的参数数量(权重和偏差)以及模型中的总参数。

autoencoder.summary()

_________________________________________________________________Layer (type)Output ShapeParam #=================================================================input_2 (InputLayer)(None, 176, 176, 1)0_________________________________________________________________conv2d_18 (Conv2D)(None, 176, 176, 32)320_________________________________________________________________batch_normalization_11 (Batc (None, 176, 176, 32)128_________________________________________________________________conv2d_19 (Conv2D)(None, 176, 176, 32)9248_________________________________________________________________batch_normalization_12 (Batc (None, 176, 176, 32)128_________________________________________________________________max_pooling2d_5 (MaxPooling2 (None, 88, 88, 32)0_________________________________________________________________conv2d_20 (Conv2D)(None, 88, 88, 64)18496_________________________________________________________________batch_normalization_13 (Batc (None, 88, 88, 64)256_________________________________________________________________conv2d_21 (Conv2D)(None, 88, 88, 64)36928_________________________________________________________________batch_normalization_14 (Batc (None, 88, 88, 64)256_________________________________________________________________max_pooling2d_6 (MaxPooling2 (None, 44, 44, 64)0_________________________________________________________________conv2d_22 (Conv2D)(None, 44, 44, 128)73856_________________________________________________________________batch_normalization_15 (Batc (None, 44, 44, 128)512_________________________________________________________________conv2d_23 (Conv2D)(None, 44, 44, 128)147584_________________________________________________________________batch_normalization_16 (Batc (None, 44, 44, 128)512_________________________________________________________________conv2d_24 (Conv2D)(None, 44, 44, 64)73792_________________________________________________________________batch_normalization_17 (Batc (None, 44, 44, 64)256_________________________________________________________________conv2d_25 (Conv2D)(None, 44, 44, 64)36928_________________________________________________________________batch_normalization_18 (Batc (None, 44, 44, 64)256_________________________________________________________________up_sampling2d_5 (UpSampling2 (None, 88, 88, 64)0_________________________________________________________________conv2d_26 (Conv2D)(None, 88, 88, 32)18464_________________________________________________________________batch_normalization_19 (Batc (None, 88, 88, 32)128_________________________________________________________________conv2d_27 (Conv2D)(None, 88, 88, 32)9248_________________________________________________________________batch_normalization_20 (Batc (None, 88, 88, 32)128_________________________________________________________________up_sampling2d_6 (UpSampling2 (None, 176, 176, 32)0_________________________________________________________________conv2d_28 (Conv2D)(None, 176, 176, 1)289=================================================================Total params: 427, 713Trainable params: 426, 433Non-trainable params: 1, 280_________________________________________________________________

终于是时候用Keras的fit()函数训练模型了！该模型训练了50个纪元。 fit()函数将返回一个历史对象；通过在fashion_train中讲述该函数的结果, 你以后可以使用它在训练和验证之间绘制损失函数图, 这将帮助你直观地分析模型的性能。
训练模型

autoencoder_train = autoencoder.fit(train_X, train_ground, batch_size=batch_size, epochs=epochs, verbose=1, validation_data=http://www.srcmini.com/(valid_X, valid_ground))

Train on 1224 samples, validate on 306 samplesEpoch 1/3001224/1224 [==============================] - 7s - loss: 0.1201 - val_loss: 0.0838Epoch 2/3001224/1224 [==============================] - 7s - loss: 0.0492 - val_loss: 0.0534...Epoch 299/3001224/1224 [==============================] - 7s - loss: 1.3101e-04 - val_loss: 6.1086e-04Epoch 300/3001224/1224 [==============================] - 7s - loss: 1.0711e-04 - val_loss: 3.9641e-04

最后！你在指纹数据集上训练了200个时期的模型, 现在, 让我们在训练和验证数据之间绘制损失图, 以可视化模型的性能。

loss = autoencoder_train.history['loss']val_loss = autoencoder_train.history['val_loss']epochs = range(300)plt.figure()plt.plot(epochs, loss, 'bo', label='Training loss')plt.plot(epochs, val_loss, 'b', label='Validation loss')plt.title('Training and validation loss')plt.legend()plt.show()

文章图片
最后, 你可以看到验证损失和训练损失都是同步的。它表明你的模型不是过拟合的：验证损失正在减少并且没有增加, 并且在整个训练阶段, 培训和验证损失之间几乎没有任何差距。
因此, 可以说模型的泛化能力很好。
最后, 是时候使用Keras的predict()函数重建测试图像了, 看看你的模型在测试数据上的重建能力如何。
保存模型现在让我们保存经过训练的模型。这是你使用深度学习的重要一步, 因此, 权重是你所要解决的问题的解决方案的核心！
你可以随时将保存的权重加载到同一模型中, 并从停止训练的地方对其进行训练。例如：如果再次训练以上模型, 则权重, 偏差, 损失函数等参数将不会从头开始, 也不再是新鲜的训练。
只需一行代码, 你就可以保存权重并将其加载回模型中。

autoencoder = autoencoder.save_weights('autoencoder_mri.h5')

autoencoder = Model(input_img, autoencoder(input_img))

autoencoder.load_weights('autoencoder_mri.h5')

从这里开始, 你没有测试数据。让我们使用验证数据对你刚刚训练的模型进行预测。
你将在剩余的306个验证图像上预测经过训练的模型, 并绘制少量重建图像以可视化模型能够多么有效地重建验证图像。

pred = autoencoder.predict(valid_X)

plt.figure(figsize=(20, 4))print("Test Images")for i in range(5):plt.subplot(1, 5, i+1)plt.imshow(valid_ground[i, ..., 0], cmap='gray')plt.show()plt.figure(figsize=(20, 4))print("Reconstruction of Test Images")for i in range(5):plt.subplot(1, 5, i+1)plt.imshow(pred[i, ..., 0], cmap='gray')plt.show()

Test Images

文章图片

Reconstruction of Test Images

文章图片
从上图可以看出, 你的模型在重建使用该模型预测的测试图像方面做得非常出色。至少在质量上, 测试和重建的图像看起来几乎完全相似。
也许它可以在重构原始3T图像中存在的一些局部细节方面做得更好。
在嘈杂的3T图像上预测首先, 让我们在验证图像中添加一些噪声, 平均值为零, 标准偏差为0.03。

[a, b, c, d]= np.shape(valid_X)mean = 0sigma = 0.03gauss = np.random.normal(mean, sigma, (a, b, c, d))noisy_images = valid_X + gauss

是时候对嘈杂的验证图像进行预测了。让我们看看你的模型尽管在噪点图像上没有经过训练, 但在噪点图像上的表现如何。

pred_noisy = autoencoder.predict(noisy_images)

plt.figure(figsize=(20, 4))print("Noisy Test Images")for i in range(5):plt.subplot(1, 5, i+1)plt.imshow(noisy_images[i, ..., 0], cmap='gray')plt.show()plt.figure(figsize=(20, 4))print("Reconstruction of Noisy Test Images")for i in range(5):plt.subplot(1, 5, i+1)plt.imshow(pred_noisy[i, ..., 0], cmap='gray')plt.show()

Noisy Test Images

文章图片

Reconstruction of Noisy Test Images

文章图片
该模型似乎做得很好。是不是毕竟, 重建后的图像看起来要比嘈杂的图像好。而且, 该模型在训练过程中从未真正看到过嘈杂的图像。
定量指标：峰值信噪比(PSNR) PSNR块计算两个图像之间的峰值信噪比, 以分贝(dB)为单位。该比率通常用作原始图像和重建图像之间的质量度量。 PSNR越高, 重建图像的质量越好。
因此, 首先让我们计算验证和重构的验证图像的性能。

valid_pred = autoencoder.predict(valid_X)mse =np.mean((valid_X - valid_pred) ** 2)psnr = 20 * math.log10( 1.0 / math.sqrt(mse))

print('PSNR of reconstructed validation images: {psnr}dB'.format(psnr=np.round(psnr, 2)))

PSNR of reconstructed validation images: 34.02dB

接下来, 让我们计算原始验证和预测的噪声图像之间的PSNR

noisy_pred = autoencoder.predict(noisy_images)mse =np.mean((valid_X - noisy_pred) ** 2)psnr_noisy = 20 * math.log10( 1.0 / math.sqrt(mse))

print('PSNR of reconstructed validation images: {psnr}dB'.format(psnr=np.round(psnr_noisy, 2)))

PSNR of reconstructed validation images: 32.48dB

好吧, 从数量上看, 没有噪声重建和有噪声重建的PSNR之间仅相差1.54 dB。此外, 你尚未训练模型来处理噪声。这不是很神奇吗？
本教程是了解如何使用脑部MRI 3T数据集读取MRI nifti格式图像, 进行分析, 预处理并将它们供入模型的良好起点。它向你展示了自动编码器的实用应用之一。如果你能够轻松地进行工作, 甚至只需付出更多的努力, 那就做好了！
你可以尝试使用该体系结构, 并尝试从数量和质量上改进预测。可能是增加了更多的层, 还是训练了更多的时间？尝试这些组合, 看看是否有帮助。
【使用深度学习(卷积自动编码器)重建大脑MRI图像】还有很多内容要讲, 为什么不参加srcmini的Python深度学习课程呢？如果你还没有这样做的话。你将从基础知识中学习, 并慢慢进入精通深度学习领域的领域, 当你学习如何在Python中使用卷积神经网络, 如何检测面部, 物体等时, 毫无疑问它将是必不可少的资源。