使用多解码器框架实现图像超分辨率

在“训练脚本”部分，您将使用一种新颖的深度学习架构来处理图像超分辨率问题。具体任务是将低场强 3 特斯拉脑部 MR 图像非线性映射到高场强 7 特斯拉脑部 MR 图像。在本“测试脚本”部分，您将使用第 1 部分训练得到的权重，对未见过的训练数据进行预测。您还将学习如何最终将 2D 图像保存为一个合并的体数据。

教程分为两部分：第一部分带您完成训练流程，第二部分讲解测试流程。

第 1 部分：训练脚本
第 2 部分：测试脚本

注意：如果您有兴趣阅读论文，可在此处找到原文。

训练脚本

概括来说，您将在本教程第 1 部分中完成以下内容：

首先导入训练深度学习模型所需的模块，
接着简要介绍 3T 与 7T MRI 数据集，
然后定义初始化参数并加载 3T 与 7T 数据集；在加载过程中您还将即时调整图像大小，
下一步对加载的数据进行预处理：将训练与测试列表转换为 numpy 矩阵，将矩阵类型转换为 float32，使用最大-最小策略进行缩放，重塑数组，最后按 80%/20% 划分训练集与验证集，
随后创建1 编码器 - 3 解码器的架构：包含合并连接与多解码器，
然后定义损失函数，建立三种不同的模型并完成编译，
最后，训练您的合并多解码器模型，在验证数据上测试，并计算定量结果。

Python 模块依赖

在开始学习本教程之前，请确保您的模块版本与下列完全一致：

Keras==2.0.4
tensorflow==1.8.0
scipy==0.19.0
numpy==1.14.5
Pillow==4.1.1
nibabel==2.1.0
scikit_learn==0.18.1

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

在下面的代码中，您通过 os.environ 在笔记本中设置环境变量。最好在初始化 Keras 之前这样做，以限制其后端 TensorFlow 仅使用第一块 GPU。如果您的训练机器的 GPU 在 0 号位置，请确保将 1 改为 0。您可以在终端运行一个简单命令进行检查，例如 nvidia-smi

import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"]="0" #model will be trained on GPU
技术支持

导入模块

首先导入所有必需模块，如 tensorflow、numpy，尤其是 keras 以及将要用到的函数或层，如 Input、Conv2D、MaxPooling2D 等，因为您将使用这些框架来训练模型！
为读取 nifti 格式图像，您还需导入名为 nibabel 的模块。

import os
import numpy as np
import math
import tensorflow as tf
import nibabel as nib
import numpy as np
from keras.layers import Input,Dense,merge,Reshape,Conv2D,MaxPooling2D,UpSampling2D
from keras.layers.normalization import BatchNormalization
from keras.models import Model,Sequential
from keras.callbacks import ModelCheckpoint
from keras.optimizers import RMSprop
from keras import backend as K
import scipy.misc
from sklearn.utils import shuffle
from sklearn.cross_validation import train_test_split
import matplotlib.pyplot as plt
from keras.models import model_from_json
技术支持

Using TensorFlow backend.
/usr/local/lib/python3.5/dist-packages/sklearn/cross_validation.py:41: 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 is different from that of this module. This module will be removed in 0.20.
  "This module will be removed in 0.20.", DeprecationWarning)
技术支持

了解脑部 MRI 3T 与 7T 数据集

脑部 MRI 3T 与 7T 数据集由 3D 体数据组成，每个体数据共包含 207 张/帧脑部 MRI 切片，取自大脑的不同层面。每张切片的尺寸为 173 x 173。图像为单通道灰度图像。数据集中共有 39 名受试者，每位受试者包含一份患者的 MRI 扫描。图像格式不是 jpeg、png 等，而是 nifti 格式。稍后您将看到如何读取 nifti 格式图像。

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

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

使用 28 位受试者的数据进行训练，其余 11 位用于测试。

定义初始化参数

先定义数据的维度。您将在读数阶段将图像尺寸从 173x173 调整到 176x176。同时定义数据目录、训练模型时的批大小、通道数、一个 Input() 层、训练与测试矩阵（列表形式），以及为缩放而加载的包含 MRI 数据集最小值与最大值的文本文件。

注意：进行缩放时，您也可以直接取数据集的最大值与最小值，因为您可能没有 maxANDmin.txt 文件。

x,y = 173,173
full_z = 207
resizeTo=176
batch_size = 32
inChannel = outChannel = 1
input_shape=(x,y,inChannel)
input_img = Input(shape = (resizeTo, resizeTo, inChannel))
inp = "ground3T"
out = "ground7T"     
train_matrix = []
test_matrix = []
min_max = np.loadtxt('maxANDmin.txt')技术支持

加载数据

接下来使用 nibabel 库加载 MRI 数据，并在 x 与 y 维度通过零填充将图像从 173 x 173 调整为 176 x 176。

注意，当您加载一个 Nifti 格式体数据时，Nibabel 并不会立即加载图像数组，而是在您请求数组数据时才加载。标准做法是调用 get_data() 方法。

由于您需要 2D 切片而非 3D 体数据，您将使用之前初始化的 train 和 test 列表；每次读取一个体数据时，遍历该 3D 体数据的全部 207 张切片，并将每张切片逐一追加到列表中。

folder = os.listdir(inp)
技术支持

我们先加载 3T 数据。

for f in folder:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load(inp + f)
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        train_matrix.append(a[:,j,:])
技术支持

然后加载 7T 数据。由于 3T 与 7T MR 体数据数量相等，加载 7T 数据时也复用同一个 folder 变量。

for f in folder:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    b = nib.load(out + f)
    b = b.get_data()
    temp[3:,:,3:] = b
    b = temp
    for j in range(full_z):
        test_matrix.append(b[:,j,:])
技术支持

数据预处理

由于训练与测试矩阵是列表，您将使用 numpy 将其转换为 numpy 数组。

进一步地，将 numpy 数组的 type 转为 float32，并对 输入与真实值同时进行缩放。

train_matrix = np.asarray(train_matrix)
train_matrix = train_matrix.astype('float32')
m = min_max[0]
mi = min_max[1]
train_matrix = (train_matrix - mi) / (m - mi)

test_matrix = np.asarray(test_matrix)
test_matrix = test_matrix.astype('float32')
test_matrix = (test_matrix - mi) / (m - mi)
技术支持

快速打印 train_matrix（3T）与 test_matrix（7T）的形状。它们应各有 28 x 207 = 5796 张图像，每张尺寸为 176 x 176。

train_matrix.shape
技术支持

(5796, 176, 176)
技术支持

test_matrix.shape
技术支持

(5796, 176, 176)
技术支持

接着，创建两个新变量 augmented_images（3T/输入）与 Haugmented_images（7T/真实值），其形状与训练与测试矩阵一致。这将是一个 4D 矩阵：第一维为图像总数，第二与第三维为每张图像的尺寸，最后一维为通道数（此处为 1）。

augmented_images=np.zeros(shape=[(train_matrix.shape[0]),(train_matrix.shape[1]),(train_matrix.shape[2]),(1)])
Haugmented_images=np.zeros(shape=[(train_matrix.shape[0]),(train_matrix.shape[1]),(train_matrix.shape[2]),(1)])
技术支持

然后逐一遍历所有图像，每次将训练与测试矩阵重塑为 176 x 176，并分别追加到 augmented_images（3T/输入）与 Haugmented_images（7T/真实值）。

for i in range(train_matrix.shape[0]):
    augmented_images[i,:,:,0] = train_matrix[i,:,:].reshape(resizeTo,resizeTo)
    Haugmented_images[i,:,:,0] = test_matrix[i,:,:].reshape(resizeTo,resizeTo)
技术支持

完成以上步骤后，划分数据非常重要。为了让模型具备良好的泛化能力，您将数据分为训练集与验证集两部分。使用 80% 的数据进行训练，剩余 20% 用于验证。

这也有助于降低过拟合的可能性，因为您将在模型未见过的数据上进行验证。

您可以使用开头导入的 scikit-learn 中的 train_test_split 模块来正确划分数据：

data,Label = shuffle(augmented_images,Haugmented_images, random_state=2)
X_train, X_test, y_train, y_test = train_test_split(data, Label, test_size=0.2, random_state=2)
技术支持

X_test = np.array(X_test)
y_test = np.array(y_test)
X_test = X_test.astype('float32')
y_test = y_test.astype('float32')
技术支持

X_train = np.array(X_train)
y_train = np.array(y_train)
X_train = X_train.astype('float32')
y_train = y_train.astype('float32')
技术支持

模型：1 个编码器与 3 个解码器！

图：选择性自编码器的反向传播

图片来源于这篇论文。

合并连接

接下来，您在编码器部分按照如下单元定义所示，使用一系列卷积层块后接最大池化层来构建所提架构。为在输出端重建图像的原始尺寸，在每个解码器块中引入上采样层。在上采样过程中，由于解码器的下采样输入缺少细节，可能引入伪影。因此，您将解码器的输入与来自编码器的上采样版本进行拼接，以在解码器上采样时提供细节信息，从而获得更好的重建效果。事实上，加入合并连接会显著提升 PSNR（约 5 dB 的量级）。该架构设置受此论文启发。

多解码器

该方法采用单编码器、多个解码器的架构，并以单通道输入。编码器与三个解码器的每个块中都包含三层卷积，随后接批归一化层以保持数值稳定性。

编码器中的第一个卷积块含 32 个滤波器，之后每个卷积块的滤波器数量翻倍。在所有解码器中，第一个卷积块含 256 个滤波器，之后每个块滤波器数量减半。所有卷积块均使用 3x3 的滤波器大小。

除最后一层外，其余各层均使用 ReLU 作为激活函数。由于数据被归一化到 0 到 1 之间，最终层使用 Sigmoid 激活。

众所周知，不同尺度的局部图像细节对图像重建至关重要。所提架构通过在编码器与解码器中分别使用层级化的下采样（最大池化）与上采样（均为因子 2），在不同尺度上处理图像。经过三次下采样得到的编码表示将高维输入压缩为潜在空间表示。

编码器

def encoder(input_img):
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(input_img)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(pool3)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv4)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='sigmoid', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    return conv5,conv4,conv3,conv2,conv1
技术支持

解码器 1

def decoder_1(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_1 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_1
技术支持

解码器 2

def decoder_2(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_2 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_2
技术支持

解码器 3

def decoder_3(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded_3 = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded_3
技术支持

在以上 4 个代码单元中，您共定义了四个函数：一个编码器与三个解码器。由于它们是函数，您本可以只定义一次 decoder() 并调用三次，但为了更好地理解且避免三个解码器随机性带来的问题，这里分别定义三次。

损失函数

接下来采用均方误差，但会将真实值 y_t 与预测值 y_p 中等于 0 的像素排除在外。

def root_mean_sq_GxGy(y_t, y_p):
    a1=1
    where = tf.not_equal(y_t, 0)
    a_t=tf.boolean_mask(y_t,where,name='boolean_mask')
    a_p=tf.boolean_mask(y_p,where,name='boolean_mask')
    return a1*(K.sqrt(K.mean((K.square(a_t-a_p)))))
技术支持

模型定义与编译

首先调用 encoder 函数并传入输入。由于架构中使用了合并连接，encoder 将返回五个 卷积层的输出，随后分别与三个解码器的输出进行合并。

conv5,conv4,conv3,conv2,conv1 = encoder(input_img)
技术支持

autoencoder_1 = Model(input_img, decoder_1(conv5,conv4,conv3,conv2,conv1))
autoencoder_1.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())

autoencoder_2 = Model(input_img, decoder_2(conv5,conv4,conv3,conv2,conv1))
autoencoder_2.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())

autoencoder_3 = Model(input_img, decoder_3(conv5,conv4,conv3,conv2,conv1))
autoencoder_3.compile(loss=root_mean_sq_GxGy, optimizer = RMSprop())
技术支持

/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:2: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.

/usr/local/lib/python3.5/dist-packages/keras/legacy/layers.py:460: UserWarning: The `Merge` layer is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
  name=name)
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:10: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
  # Remove the CWD from sys.path while we load stuff.
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:18: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.
/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:26: UserWarning: The `merge` function is deprecated and will be removed after 08/2017. Use instead layers from `keras.layers.merge`, e.g. `add`, `concatenate`, etc.


WARNING:tensorflow:From /usr/local/lib/python3.5/dist-packages/keras/backend/tensorflow_backend.py:1257: calling reduce_mean (from tensorflow.python.ops.math_ops) with keep_dims is deprecated and will be removed in a future version.
Instructions for updating:
keep_dims is deprecated, use keepdims instead
技术支持

训练模型

仅当验证数据上的峰值信噪比（PSNR）提升时才保存权重。为此您将定义一个 psnr_gray_channel 列表，并先追加默认值 1。之所以追加 1，是因为需要用一个随机数初始化列表，并且在训练初期 PSNR 也会远高于该值，它仅作为一个占位用的虚拟数值。

将初始学习率设为 1e-3，并采用学习率衰减策略：每 20 个 epoch 将当前学习率降低 10%！

psnr_gray_channel = []
psnr_gray_channel.append(1)
learning_rate = 0.001
j=0
技术支持

注意：下一个完整代码需在同一单元中一次性运行，但为帮助理解整个训练流程，这里会拆分为多个小单元展示。

模型训练 500 个 epoch。初始学习率为上一单元定义的 1e-3。您还将每个 epoch 后保存 7T 与预测结果之间的 PSNR、MSE。使用 K.set_value 函数每 20 个 epoch 将三个模型的学习率统一降低 10%。

for jj in range(500):
    myfile_valid_psnr_7T = open('../1_encoder_3_decoders_complete_slices_single_channel/validation_psnr7T_1encoder_3decoders.txt', 'a')
    myfile_valid_mse_7T = open('../1_encoder_3_decoders_complete_slices_single_channel/validation_mse7T_1encoder_3decoders.txt', 'a')

    K.set_value(autoencoder_1.optimizer.lr, learning_rate)
    K.set_value(autoencoder_2.optimizer.lr, learning_rate)
    K.set_value(autoencoder_3.optimizer.lr, learning_rate)
技术支持

接下来打乱输入的 3T 与真实 7T 图像，以避免网络过拟合，因为不打乱会使模型在每个 epoch 都以相同顺序看到样本。然后基于先前定义的 batch_size 计算将生成的批次数。最后开始遍历 num_batches。

train_X,train_Y = shuffle(X_train,y_train)
print ("Epoch is: %d\n" % j)
print ("Number of batches: %d\n" % int(len(train_X)/batch_size))
num_batches = int(len(train_X)/batch_size)
for batch in range(num_batches):
技术支持

除了存储 PSNR 值，您还将分别存储三个自编码器及其三个解码器的损失。

myfile_ae1_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae1_train_loss_1encoder_3decoders.txt', 'a')
myfile_ae2_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae2_train_loss_1encoder_3decoders.txt', 'a')
myfile_ae3_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/ae3_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec1_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec1_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec2_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec2_train_loss_1encoder_3decoders.txt', 'a')
myfile_dec3_loss = open('../1_encoder_3_decoders_complete_slices_single_channel/dec3_train_loss_1encoder_3decoders.txt', 'a')
技术支持

由于每个批次都希望模型看到接下来的 32（batch_size）个样本，下一单元将完成该操作！

batch_train_X = train_X[batch*batch_size:min((batch+1)*batch_size,len(train_X)),:]
batch_train_Y = train_Y[batch*batch_size:min((batch+1)*batch_size,len(train_Y)),:]
技术支持

为使最小自编码器策略生效，在每个 epoch 后，您使用 Keras 的 test_on_batch 函数在训练数据上测试，得到三个损失，并打印出来。

loss_1 = autoencoder_1.test_on_batch(batch_train_X,batch_train_Y)
loss_2 = autoencoder_2.test_on_batch(batch_train_X,batch_train_Y)
loss_3 = autoencoder_3.test_on_batch(batch_train_X,batch_train_Y)
print ('epoch_num: %d batch_num: %d Test_loss_1: %f\n' % (j,batch,loss_1))
print ('epoch_num: %d batch_num: %d Test_loss_2: %f\n' % (j,batch,loss_2))
print ('epoch_num: %d batch_num: %d Test_loss_3: %f\n' % (j,batch,loss_3))
技术支持

网络中可能出现六种情况：
loss_1：自编码器 1 的损失
loss_2：自编码器 2 的损失
loss_3：自编码器 3 的损失

loss_1 大于 loss_2 与 loss_3。如果为真，则仅使用自编码器 1 训练。将自编码器 2 与自编码器 3 的编码器部分设为 False，只训练其解码器，最后将所有损失写入文本文件。

loss_2 大于 loss_1 与 loss_3。如果为真，则仅使用自编码器 2 训练。将自编码器 1 与自编码器 3 的编码器部分设为 False，只训练其解码器，最后将所有损失写入文本文件。

loss_3 大于 loss_1 与 loss_2。如果为真，则仅使用自编码器 3 训练。将自编码器 1 与自编码器 2 的编码器部分设为 False，只训练其解码器，最后将所有损失写入文本文件。

loss_1 等于 loss_2。如果为真，则使用自编码器 1 或 2 进行训练。将自编码器 3 的编码器部分设为 False，同时将未被选中的 1 与 2 中的那个也设为 False，只训练其解码器，最后将所有损失写入文本文件。

loss_2 等于 loss_3。如果为真，则使用自编码器 2 或 3 进行训练。将自编码器 1 的编码器部分设为 False，同时将未被选中的 2 与 3 中的那个也设为 False，只训练其解码器，最后将所有损失写入文本文件。

loss_3 等于 loss_1。如果为真，则使用自编码器 3 或 1 进行训练。将自编码器 2 的编码器部分设为 False，同时将未被选中的 3 与 1 中的那个也设为 False，只训练其解码器，最后将所有损失写入文本文件。

图：模型架构

图片来源于这篇论文。

if loss_1 < loss_2 and loss_1 < loss_3:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1: %f\n' % (j,batch,train_1))
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_loss_3: %f\n' % (j,batch,train_3))
elif loss_2 < loss_1 and loss_2 < loss_3:
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_2: %f\n' % (j,batch,train_2))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_1.summary()
    #autoencoder_3.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" % (train_1))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3: %f\n' % (j,batch,train_3))
elif loss_3 < loss_1 and loss_3 < loss_2:
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_3: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_3: %f\n' % (j,batch,train_3))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    #autoencoder_1.summary()
    #autoencoder_2.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" %(train_1))
    myfile_dec2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2: %f\n' % (j,batch,train_2))
elif loss_1 == loss_2:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
elif loss_2 == loss_3:
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae2_loss.write("%f \n" % (train_2))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    for layer in autoencoder_1.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False

    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec1_loss.write("%f \n" % (train_1))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
elif loss_3 == loss_1:
    train_1 = autoencoder_1.train_on_batch(batch_train_X,batch_train_Y)
    myfile_ae1_loss.write("%f \n" % (train_1))
    print ('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    #myfile.write('epoch_num: %d batch_num: %d AE_Train_loss_1_equal_state: %f\n' % (j,batch,train_1))
    for layer in autoencoder_2.layers[:34]:
        layer.trainable = False
    for layer in autoencoder_3.layers[:34]:
        layer.trainable = False
    #autoencoder_2.summary()
    #autoencoder_3.summary()
    train_2 = autoencoder_2.train_on_batch(batch_train_X,batch_train_Y)
    train_3 = autoencoder_3.train_on_batch(batch_train_X,batch_train_Y)
    myfile_dec2_loss.write("%f \n" % (train_2))
    myfile_dec3_loss.write("%f \n" % (train_3))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    print ('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_2_equal_state: %f\n' % (j,batch,train_2))
    #myfile.write('epoch_num: %d batch_num: %d Decoder_Train_loss_3_equal_state: %f\n' % (j,batch,train_3))


    myfile_ae1_loss.close()
    myfile_ae2_loss.close()
    myfile_ae3_loss.close()
    myfile_dec1_loss.close()
    myfile_dec2_loss.close()
    myfile_dec3_loss.close()
技术支持

这一步很关键：由于您在上述条件中将部分层设为了 False，这里需要将这些层重新设为 True，以便在使用 test_on_batch 函数时三个自编码器的所有层都可用，且不会在整个训练过程中一直保持为 False。

for layer in autoencoder_1.layers[:34]:
            layer.trainable = True
for layer in autoencoder_2.layers[:34]:
            layer.trainable = True
for layer in autoencoder_3.layers[:34]:
            layer.trainable = True
技术支持

您将通过三种方式保存权重，下面是其中两种：其一每 100 个 epoch 保存一次；其二每个 epoch 都保存（会被覆盖）。

if jj % 100 ==0:
            autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1_" + str(jj)+".h5")
            autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2_" + str(jj)+".h5")
            autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3_" + str(jj)+".h5")


autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1.h5")
autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2.h5")
autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3.h5")
技术支持

在验证数据上测试

您将首先打乱验证数据，然后复用上文定义的六种可能情况。满足哪种情况，就使用相应的自编码器来测试模型。

随后计算 MSE 与 PSNR（预测的 decoded_imgs 与真实值之间）。最后将它们保存到文本文件。

X_test,y_test = shuffle(X_test,y_test)
if loss_1 < loss_2 and loss_1 < loss_3:
    decoded_imgs = autoencoder_1.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)


    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))

    #print (check)
elif loss_2 < loss_1 and loss_2 < loss_3:
    decoded_imgs = autoencoder_2.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))

    #print (check)
elif loss_3 < loss_2 and loss_3 < loss_1:
    decoded_imgs = autoencoder_3.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))
    #print (check)

elif loss_1 == loss_2:
    decoded_imgs = autoencoder_1.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))


elif loss_2 == loss_3:
    decoded_imgs = autoencoder_2.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))


elif loss_3 == loss_1:
    decoded_imgs = autoencoder_3.predict(X_test)
    mse_7T=  np.mean((y_test[:,:,:,0] - decoded_imgs[:,:,:,0]) ** 2)
    check_7T = math.sqrt(mse_7T)
    psnr_7T = 20 * math.log10( 1.0 / check_7T)

    myfile_valid_psnr_7T.write("%f \n" % (psnr_7T))
    myfile_valid_mse_7T.write("%f \n" % (mse_7T))
技术支持

此处仅当预测的 7T 与真实值（7T）之间的 PSNR 相较于此前存储在 psnr_gray_channel 列表中的最大 PSNR 更高时才保存权重。

if max(psnr_gray_channel) < psnr_7T:
            autoencoder_1.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE1_" + str(jj)+".h5")
            autoencoder_2.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE2_" + str(jj)+".h5")
            autoencoder_3.save_weights("../Model/CROSSVAL1/CROSSVAL1_AE3_" + str(jj)+".h5")

    psnr_gray_channel.append(psnr_7T)
技术支持

保存输入、真实值与解码结果：定量结果

您将定义一个大小为 176 x 528 的 numpy 矩阵 temp，因为将在一行中保存 3 张图像，每张尺寸 176 x 176。您将从所有验证集 3T、7T 与预测结果中任取一张图像进行保存，并将矩阵乘以 255，因为图像之前被缩放到了 0 到 1 之间。

最后使用 scipy 将图像保存到文件。

temp = np.zeros([resizeTo,resizeTo*3])
temp[:resizeTo,:resizeTo] = X_test[0,:,:,0]
temp[:resizeTo,resizeTo:resizeTo*2] = y_test[0,:,:,0]
temp[:resizeTo,2*resizeTo:] = decoded_imgs[0,:,:,0]
temp = temp*255
scipy.misc.imsave('../Results/1_encoder_3_decoders_complete_slices_single_channel/' + str(j) + '.jpg', temp)
j +=1
技术支持

最后关闭 PSNR 与 MSE 文件。

myfile_valid_psnr_7T.close()
myfile_valid_mse_7T.close()
技术支持

最后，每 20 个 epoch 将学习率在当前值基础上降低 10%。

if jj % 20 ==0:
        learning_rate = learning_rate - learning_rate * 0.10
技术支持

测试脚本

在本教程第 2 部分，您将完成以下内容：

首先导入训练深度学习模型所需的模块，
接着简要介绍 3T 与 7T MRI 测试数据集，
然后定义初始化参数并加载 3T 与 7T 测试数据集；在加载过程中即时调整图像大小，
接着对加载的数据进行预处理：将训练与测试列表转换为 numpy 矩阵，将矩阵类型转换为 float32，使用最大-最小策略进行缩放，重塑数组，最后按 80%/20% 划分训练集与验证集，
随后创建1 编码器 - 3 解码器的架构：包含合并连接与多解码器，
然后定义：损失函数、三种不同的模型，并加载训练好的权重，
最后，使用合并多解码器模型在未见数据上进行预测，并保存定量与定性结果。您还将学习如何使用 nibabel 将 2D 图像保存为合并体数据。

Python 模块依赖

在开始学习本教程之前，请确保您的模块版本与下列完全一致：

Keras==2.0.4
tensorflow==1.8.0
scipy==0.19.0
numpy==1.14.5
Pillow==4.1.1
nibabel==2.1.0
scikit_learn==0.18.1

在下面的代码中，您基本上通过 os.environ 在笔记本中设置环境变量。最好在初始化 Keras 之前这样做，以限制其后端 TensorFlow 仅使用第一块 GPU。如果您的训练机器的 GPU 在 0 号位置，请确保将 1 改为 0。您可以在终端运行一个简单命令进行检查，例如 nvidia-smi

import os
os.environ["CUDA_DEVICE_ORDER"]="PCI_BUS_ID"
os.environ["CUDA_VISIBLE_DEVICES"]="0" #model will be trained on GPU
技术支持

导入模块

import os
from keras.layers import Input,Dense,Flatten,Dropout,merge,Reshape,Conv2D,MaxPooling2D,UpSampling2D
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
import numpy as np
import scipy.misc
import numpy.random as rng
from sklearn.utils import shuffle
import nibabel as nib
from sklearn.cross_validation import train_test_split
import math
技术支持

了解脑部 MRI 3T 与 7T 数据集

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

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

定义初始化参数

先定义数据维度。您将在读数阶段将图像尺寸从 173x173 调整到 176x176。这里还将定义数据目录、训练时的批大小、通道数、一个 Input() 层、训练与测试矩阵（列表形式），最后为缩放加载包含 MRI 数据集最小值与最大值的文本文件。

x,y = 173,173
full_z = 207
resizeTo=176
inChannel = outChannel = 1
input_shape=(x,y,inChannel)
input_img = Input(shape = (resizeTo, resizeTo, inChannel))
技术支持

train_matrix = []
test_matrix = []
ff = os.listdir("../test_crossval1")
save = "../Result_nii_crossval1"
folder_ground = os.listdir("../test_g_crossval1")
ToPredict_images=[]
predict_matrix=[]
技术支持

ground_images=[]
ground_matrix=[]
min_max=np.loadtxt('../maxANDmin.txt')
技术支持

加载测试体数据

接下来使用 nibabel 库加载 MRI 数据，并在 x 与 y 维度通过零填充将图像从 173 x 173 调整为 176 x 176。

注意，当您加载一个 Nifti 格式体数据时，Nibabel 并不会立即加载图像数组，而是在您请求数组数据时才加载。标准做法是调用 get_data() 方法。

for f in ff:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load("../test_crossval1" + f)
    affine = a.affine
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        predict_matrix.append(a[:,j,:])
技术支持

for f in ff:
    temp = np.zeros([resizeTo,full_z,resizeTo])
    a = nib.load("../test_g_crossval1" + f)
    affine = a.affine
    a = a.get_data()
    temp[3:,:,3:] = a
    a = temp
    for j in range(full_z):
        ground_matrix.append(a[:,j,:])
技术支持

数据预处理

由于 3T 与 7T 测试矩阵为列表，您将使用 numpy 将其转换为 numpy 数组。

进一步地，将 numpy 数组的 type 转为 float32，并对 输入与真实值同时进行缩放。

ToPredict_images = np.asarray(predict_matrix)
ToPredict_images = ToPredict_images.astype('float32')
mx = min_max[0]
mn = min_max[1]
ToPredict_images[:,:,:,0] = (ToPredict_images[:,:,:,0] - mn ) / (mx - mn)
技术支持

ground_images = np.asarray(ground_matrix)
ground_images = ground_images.astype('float32')
ground_images[:,:,:,0] = (ground_images[:,:,:,0] - mn ) / (mx - mn)
技术支持

接着，创建两个新变量 ToPredict_images（3T 测试/输入）与 ground_images（7T 测试/真实值），其形状与训练与测试矩阵一致。这将是一个 4D 矩阵：第一维为图像总数，第二与第三维为每张图像的尺寸，最后一维为通道数（此处为 1）。

ToPredict_images=np.zeros(shape=[(ToPredict_images.shape[0]),(ToPredict_images.shape[1]),(ToPredict_images.shape[2]),(1)])
ground_images=np.zeros(shape=[(ground_images.shape[0]),(ground_images.shape[1]),(ground_images.shape[2]),(1)])
技术支持

然后逐一遍历所有图像，每次将训练与测试矩阵重塑为 176 x 176，并分别追加到 ToPredict_images（3T 测试/输入）与 ground_images（7T 测试/真实值）。

for i in range(ToPredict_images.shape[0]):
    ToPredict_images[i,:,:,0] = ToPredict_images[i,:,:].reshape(resizeTo,resizeTo)
技术支持

for i in range(ground_images.shape[0]):
    ground_images[i,:,:,0] = ground_images[i,:,:].reshape(resizeTo,resizeTo)
技术支持

模型！

图：模型架构

图片来源于这篇论文。

编码器

def encoder(input_img):
    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(input_img)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    conv1 = Conv2D(32, (3,3), activation='relu', padding='same')(conv1)
    conv1 = BatchNormalization()(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)
    conv2 = BatchNormalization()(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)
    conv3 = BatchNormalization()(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(pool3)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)
    conv4 = BatchNormalization()(conv4)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv4)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    conv5 = Conv2D(512, (3, 3), activation='sigmoid', padding='same')(conv5)
    conv5 = BatchNormalization()(conv5)
    return conv5,conv4,conv3,conv2,conv1
技术支持

解码器

def decoder(conv5,conv4,conv3,conv2,conv1):
    up6 = merge([conv5, conv4], mode='concat', concat_axis=3)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)
    conv6 = BatchNormalization()(conv6)
    up7 = UpSampling2D((2,2))(conv6)
    up7 = merge([up7, conv3], mode='concat', concat_axis=3)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)
    conv7 = BatchNormalization()(conv7)
    up8 = UpSampling2D((2,2))(conv7)
    up8 = merge([up8, conv2], mode='concat', concat_axis=3)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)
    conv8 = BatchNormalization()(conv8)
    up9 = UpSampling2D((2,2))(conv8)
    up9 = merge([up9, conv1], mode='concat', concat_axis=3)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)
    conv9 = BatchNormalization()(conv9)
    decoded = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(conv9)
    return decoded
技术支持

损失函数

接下来采用均方误差，但会将真实值 y_t 与预测值 y_p 中等于 0 的像素排除在外。

def root_mean_sq_GxGy(y_t, y_p):
    a1=1
    zero = tf.constant(0, dtype=tf.float32)
    where = tf.not_equal(y_t, zero)
    a_t=tf.boolean_mask(y_t,where,name='boolean_mask')
    a_p=tf.boolean_mask(y_p,where,name='boolean_mask')
    return a1*(K.sqrt(K.mean((K.square(a_t-a_p)))))
技术支持

模型定义并在三个自编码器中加载权重！

首先调用 encoder 函数并传入输入。由于架构中使用了合并连接，encoder 将返回五个 卷积层的输出，随后分别与三个解码器的输出进行合并。

conv5,conv4,conv3,conv2,conv1 = encoder(input_img)
技术支持

现在创建三个不同的模型，并将三个已训练的模型权重加载进去。

autoencoder_1 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_1.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE1.h5")
技术支持

autoencoder_2 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_2.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE2.h5")
技术支持

autoencoder_3 = Model(input_img, decoder(conv5,conv4,conv3,conv2,conv1))

autoencoder_3.load_weights("../Model/CROSSVAL1/CROSSVAL1_AE3.h5")
技术支持

模型在测试体数据上的预测：定量与定性结果！

快速初始化两个大小为 11 x 3 x 1 的 numpy 数组：第一维代表用于测试模型的体数据数量；第二维代表三种对比关系的 MSE 与 PSNR：预测输出 7T MR 图像与真实 7T MR 图像、预测输出与输入 3T MR 图像、输入 3T MR 图像与真实 7T MR 图像；第三维代表输入模型的通道数！

mse= np.zeros([11,3,1])
psnr= np.zeros([11,3,1])
i=0 #for iterating over the slices of all the 11 volumes
技术支持

在下一部分中，您将逐个遍历全部 11 个体数据。每次迭代中，分别使用三个自编码器对该体数据进行预测，最后对预测结果取平均。

在每个体数据中，您将遍历其通道数，并计算上述三种情形下的 PSNR 与 MSE。

随后使用 nibabel 库将预测输出、输入（3T）与真实值（7T）保存为 .nii 格式文件：每个体数据包含 207 张切片。

最后使用 numpy 将 PSNR 矩阵保存到文本文件。

正如论文所述，对预测输出取平均有助于降低噪声影响，同时保留重建图像中的局部特征，因此相较于单个解码器的输出，PSNR 会有所提升。

for j in range(11):
    decoded_imgs_1 = autoencoder_1.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs_2 = autoencoder_2.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs_3 = autoencoder_3.predict(ToPredict_images[i:i+207,:,:,:])
    decoded_imgs = np.mean( np.array([ decoded_imgs_1, decoded_imgs_2,decoded_imgs_3 ]), axis=0 )
    for channel in range(1):
        mse[j,0,channel]=  np.mean((ground_images[i:i+207,:,:,channel] - decoded_imgs[:,:,:,channel]) ** 2)
        psnr[j,0,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,0,channel]))
        mse[j,1,channel]=  np.mean((ground_images[i:i+207,:,:,channel] - ToPredict_images[i:i+207,:,:,channel])** 2)
        psnr[j,1,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,1,channel]))
        mse[j,2,channel] =  np.mean((ToPredict_images[i:i+207,:,:,channel] - decoded_imgs[:,:,:,channel]) ** 2)
        checklt = math.sqrt(mse[j,2,channel])
        psnr[j,2,channel] = 20 * math.log10( 1.0 / math.sqrt(mse[j,2,channel]))
    obj = nib.Nifti1Image(decoded_imgs, affine)
    string =str(j)+'_crossval1.nii'
    nib.save(obj, save + string)
    obj = nib.Nifti1Image(ground_images[i:i+207,:,:,:], affine)
    string =str(j)+'_ground_images_crossval1.nii'
    nib.save(obj, save + string)
    obj = nib.Nifti1Image(ToPredict_images[i:i+207,:,:,:], affine)
    string =str(j)+'_ToPredict_images_crossval1.nii'
    nib.save(obj, save + string)
    i=i+207


np.savetxt('psnr_all_slices.txt',psnr[:,:,0])
技术支持

如果您希望进一步学习 Python，请查看 DataCamp 的Matplotlib 数据可视化入门课程。

也可参阅Keras 教程：用 Python 进行深度学习。

主题

Python 机器学习深度学习

训练脚本

Python 模块依赖

导入模块

了解脑部 MRI 3T 与 7T 数据集

定义初始化参数

加载数据

数据预处理

模型：1 个编码器与 3 个解码器！

合并连接

多解码器

编码器

解码器 1

解码器 2

解码器 3

损失函数

模型定义与编译

训练模型

在验证数据上测试

保存输入、真实值与解码结果：定量结果

测试脚本

Python 模块依赖

导入模块

了解脑部 MRI 3T 与 7T 数据集

定义初始化参数

加载测试体数据

数据预处理

模型！

编码器

解码器

损失函数

模型定义并在三个自编码器中加载权重！

模型在测试体数据上的预测：定量与定性结果！

主题

深入了解 Python 与深度学习

Keras 深度学习入门

使用 Keras 的高级深度学习

Python 深度学习入门

深入了解 Python 与深度学习

Keras 深度学习入门

使用 Keras 的高级深度学习

Python 深度学习入门