跳至主要内容

使用深度学习（卷积自编码器）重建脑部 MRI 图像

在本教程中，您将学习并理解如何读取 NIfTI 格式的脑部磁共振成像（MRI）图像，并使用卷积自编码器对其进行重建。
· 13分钟 阅读

使用 AI 探索

您将使用 3T 脑部 MRI 数据集来训练您的网络。为观察模型的有效性，您将把模型测试在：

  • 未见过的 3T MRI 图像，
  • 带噪声的 3T MRI 图像，以及
  • 使用定性指标：峰值信噪比（PSNR）来评估重建图像的性能。

本教程不会涉及医学影像的复杂细节，而是专注于深度学习部分！注意：本教程主要涵盖卷积自编码器的实践实现。因此，如果您尚不了解卷积神经网络（CNN）和自编码器，您可能需要先查看 CNN 和 自编码器 教程。 本教程的最大亮点在于，您将把 3D 体数据以 2D 图像的形式加载并输入模型。 总之，今天的教程将涵盖以下主题：

  • 首先，您将简要了解磁共振成像（MRI），
  • 然后您将理解脑部 MRI 数据集：它包含何种图像、如何导入模块、如何读取图像、创建图像数组、对脑部 MRI 图像进行预处理以便输入模型，最后探索这些图像。
  • 在卷积自编码器的实现中：您将把预处理后的数据拟合到模型中，可视化训练与验证损失曲线，保存已训练的模型，并在测试集上进行预测。
  • 接着，您将通过向测试图像中添加噪声来检验预训练模型的鲁棒性，并定量观察模型表现如何。
  • 最后，您将使用定量指标峰值信噪比 (PSNR) 来评估预测结果并衡量模型性能。

MR 图像简要介绍

医学成像中使用了多种系统，从磁场强度为 0.3 特斯拉（T）的开放式 MRI，到磁场强度可达 1.0 T 的四肢 MRI 系统，再到临床使用的全身扫描仪，其磁场强度可达 3.0 T。特斯拉是衡量 MR 图像磁场定量强度的单位。高场 MR 扫描仪（7T、11.5T）即便在体素（一个三维小块或网格）尺寸更小的情况下也能提供更高的 SNR（信噪比），因此更有利于获得更准确的诊断。

更小的体素尺寸带来更高的分辨率，进而有助于临床诊断。然而，所用 MR 扫描仪的磁场强度会对体素尺寸设置下限，以在保持 MR 图像细节的同时维持良好的信噪比（SNR）。
尽管 7T 和 11.5T 图像质量更佳，但由于成本限制，它们很少在生产环境中部署。

根据最新论文报道，3T 扫描仪约有 20,000 台，而 7T 扫描仪仅约 40 台。

理解脑部 MRI 3T 数据集

脑部 MRI 数据集由 3D 体数据组成，每个体数据包含总计 207 张/层在大脑不同切片位置获取的 MRI 图像。每一层的尺寸为 173 x 173。图像为单通道灰度图。共有 30 个受试者，每个受试者包含一名患者的 MRI 扫描。图像格式不是 jpeg、png 等，而是 NIfTI 格式。您将在后文看到如何读取 NIfTI 格式的图像。

该数据集包含 T1 模态 MR 图像，T1 序列传统上被认为适合解剖结构评估。今天您将使用的数据集为 3T 脑部 MRI。

该数据集是公开的，可从此来源下载。

提示：如果您想学习如何在 MNIST 数据集上实现用于分类任务的多层感知器（MLP），请查看本教程。

注意：在开始之前请注意，模型将于配备 Nvidia 1080 Ti GPU、Xeon e5 GeForce 处理器和 32GB 内存的系统上训练。如果您使用 Jupyter Notebook，您需要增加三行代码，通过名为 os 的模块指定 CUDA 设备顺序和可见设备。

在下面的代码中，您基本上是通过 os.environ 在笔记本中设置环境变量。建议在初始化 Keras 之前执行以下操作，以限制 Keras 后端 TensorFlow 使用第一块 GPU。如果您训练所用机器的 GPU 编号为 0，请确保使用 0 而非 1。您可以在终端中运行简单命令检查，例如：nvidia-smi

import os
os.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 os
import cv2
from keras.layers import Input,Dense,Flatten,Dropout,merge,Reshape,Conv2D,MaxPooling2D,UpSampling2D,Conv2DTranspose
from keras.layers.normalization import BatchNormalization
from keras.models import Model,Sequential
from keras.callbacks import ModelCheckpoint
from keras.optimizers import Adadelta, RMSprop,SGD,Adam
from keras import regularizers
from keras import backend as K

Using TensorFlow backend.

import numpy as np
import scipy.misc
import numpy.random as rng
from PIL import Image, ImageDraw, ImageFont
from sklearn.utils import shuffle
import nibabel as nib #reading MR images
from sklearn.cross_validation import train_test_split
import math
import glob
from 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)

数据预处理

由于 images 是一个列表，您将使用 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_split
train_X,valid_X,train_ground,valid_ground = train_test_split(images,
                                                             images,
                                                             test_size=0.2,
                                                             random_state=13)

注意，对于本任务，您不需要训练与测试标签。这就是您要把训练图像传入两次的原因。训练图像同时作为输入与类似分类任务中“标签”的真实值。

数据探索

现在让我们分析一下数据集中图像的样子，并再次查看图像的尺寸，因为您已经通过添加三行三列进行了填充。可使用 NumPy 数组的 .shape 属性：

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

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

从上述输出可见，数据的形状为 6210 x 176 x 176，因为共有 6210 个样本，每个样本是 176 x 176 x 1 的矩阵。

现在，让我们看看数据集中部分训练与验证图像：

plt.figure(figsize=[5,5])

# Display the first image in training data
plt.subplot(121)
curr_img = np.reshape(train_X[0], (176,176))
plt.imshow(curr_img, cmap='gray')

# Display the first image in testing data
plt.subplot(122)
curr_img = np.reshape(valid_X[0], (176,176))
plt.imshow(curr_img, cmap='gray')

<matplotlib.image.AxesImage at 0x7ff67d5d59b0>
images

上述两幅图分别来自训练集与验证集。您可以看到二者不同。看看卷积自编码器是否能够学习到特征并正确重建这些图像将很有意思。

现在一切就绪，可以定义网络并将数据送入网络了。事不宜迟，让我们进入下一步！

卷积自编码器！

图像尺寸为 176 x 176 x 1，或 30976 维向量。您将图像矩阵转换为数组，将其缩放到 0 到 1 之间，重塑为 176 x 176 x 1 的形状，并将其作为输入馈入网络。

此外，您将使用 128 的批大小；使用更大的批大小如 256 或 512 也可以，取决于您训练模型的硬件。这对学习参数的确定及预测准确性影响很大。您将把网络训练 50 个周期。

batch_size = 128
epochs = 300
inChannel = 1
x, y = 176, 176
input_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 32
    conv1 = 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 32
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1) #14 x 14 x 64
    conv2 = 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 64
    conv3 = 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)


    #decoder
    conv4 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv3) #7 x 7 x 128
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    up1 = UpSampling2D((2,2))(conv4) # 14 x 14 x 128
    conv5 = Conv2D(32, (3, 3), activation='relu', padding='same')(up1) # 14 x 14 x 64
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    up2 = UpSampling2D((2,2))(conv5) # 28 x 28 x 64
    decoded = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(up2) # 28 x 28 x 1
    return decoded

创建模型后，您需要编译它，优化器使用 RMSProp。

同时还需通过 loss 参数指定损失类型。在本例中为均方误差，因为每个批次后都会逐像素计算预测输出与真实值之间的均方误差：

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

使用 summary 函数可视化您在上一步创建的各层，这将显示每层的参数（权重与偏置）数量以及模型的总参数量。

autoencoder.summary()

_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
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,713
Trainable params: 426,433
Non-trainable params: 1,280
_________________________________________________________________

终于可以使用 Keras 的 fit() 函数来训练模型了！模型训练 50 个周期。fit() 函数将返回一个 history 对象；通过将其结果存储在 fashion_train 中，您可以稍后绘制训练与验证损失曲线图，从而直观分析模型表现。

训练模型

autoencoder_train = autoencoder.fit(train_X, train_ground, batch_size=batch_size,epochs=epochs,verbose=1,validation_data=(valid_X, valid_ground))

Train on 1224 samples, validate on 306 samples
Epoch 1/300
1224/1224 [==============================] - 7s - loss: 0.1201 - val_loss: 0.0838
Epoch 2/300
1224/1224 [==============================] - 7s - loss: 0.0492 - val_loss: 0.0534
...
Epoch 299/300
1224/1224 [==============================] - 7s - loss: 1.3101e-04 - val_loss: 6.1086e-04
Epoch 300/300
1224/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()
training and validation loss

可以看到，验证损失与训练损失走势一致。这表明您的模型没有过拟合：验证损失在下降而非上升，且整个训练阶段训练与验证损失之间几乎没有明显差距。

因此，可以认为模型的泛化能力较好。

最后，是时候使用 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
Test Images
Reconstruction of Test Images
Reconstruction of Test Images

从上图可以观察到，您的模型在重建测试图像方面表现出色。至少从定性角度看，测试图像与重建图像几乎完全相同。

也许在重建原始 3T 图像中的某些局部细节方面还能做得更好。

在带噪 3T 图像上进行预测

首先我们在验证图像中添加一些噪声，均值为 0、标准差为 0.03。

[a,b,c,d]= np.shape(valid_X)
mean = 0
sigma = 0.03
gauss = 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
noisy test images
Reconstruction of 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 的 NIfTI 格式图像，对其进行分析、预处理，并使用 3T 脑部 MRI 数据集将其输入模型。它向您展示了自编码器一个很不错的实际应用。如果您能轻松跟上，或稍加努力也完成了，做得很棒！

您可以尝试调整架构，从定量与定性两个方面改进预测效果。也许可以增加更多层、或者延长训练时间？尝试这些组合，看看是否有帮助。

仍有许多内容值得探索，何不尝试 DataCamp 的Python 深度学习入门课程？如果您还没有学过，建议现在开始。您将从基础学起，逐步迈向深度学习领域的进阶，这无疑会成为您学习如何在 Python 中使用卷积神经网络、进行人脸与对象检测等任务时不可或缺的资源。

主题

Python深度学习