来自伦敦帝国理工学院的 Martin Rajchl、S. Ira Ktena 和 Nick Pawlowski 的客座文章
DLTK，即医学影像深度学习工具包，扩展了 TensorFlow 以支持生物医学图像的深度学习。它提供专业操作和函数、模型实现、教程（如本博文中所用）以及典型应用的代码示例.

这篇文章是对深度学习在生物医学图像领域中应用的简要介绍，我们将展示当前工程问题中的一些问题和解决方案，并向您展示如何使用原型解决您的问题。

网站: https://dltk.github.io; 源代码: https://github.com/DLTK/DLTK;

概述

什么是生物医学图像分析，为什么需要它？ 生物医学图像是在不同尺度（例如显微镜、宏观等）上对人体的测量。它们来自各种成像模式（例如 CT 扫描仪、超声仪等），并测量人体的物理特性（例如射线密度、对 X 射线的透明度）。这些图像由领域专家（例如放射科医师）解释用于临床任务（例如诊断），并对医生的决策产生重大影响。

医学图像示例（从左上到右下）：多序列脑部 MRI：T1 加权、T1 反转恢复和 T2 FLAIR 通道；拼接全身 MRI；平面心脏超声；胸部 X 光；心脏电影 MRI。

生物医学图像通常是体积图像（3D），有时还具有额外的时域维度（4D）和/或多个通道（4-5D）（例如多序列 MR 图像）。生物医学图像的差异与自然图像（例如照片）有很大不同，因为临床协议旨在分层说明图像的获取方式（例如患者仰卧，头部没有倾斜等）。在分析中，我们的目标是检测细微差异（即指示异常发现的一些小区域）。

为什么使用计算机视觉和机器学习？ 长期以来，计算机视觉方法一直被用于自动分析生物医学图像。深度学习的最新进展已取代了许多其他机器学习方法，因为它避免了手工设计特征的创建，从而消除了该过程中一个关键的错误来源。此外，GPU 加速全连接网络的快速推理速度使我们能够将分析扩展到前所未有的数据量（例如 10⁶ 个受试者图像）。

我们可以直接将深度学习库用于生物医学成像吗？为什么要创建 DLTK？
创建 DLTK 的主要原因是将该领域的专业工具包含在开箱即用的环境中。虽然许多深度学习库向开发人员公开了低级操作（例如张量乘法等），但许多高级专业操作在体积图像上使用时是缺失的（例如可微 3D 上采样层等），并且由于图像的额外空间维度，我们可能会遇到内存问题（例如存储 1k CT 图像数据库的单个副本，图像维度为 512x512x256 个体素，在 float32 中约为 268 GB）。由于采集的性质不同，一些图像需要特殊的预处理（例如强度归一化、偏差场校正、去噪、空间归一化/配准等）。

文件格式、标题和读取图像

虽然许多成像模式供应商以DICOM标准格式生成图像，将体积保存为一系列 2D 切片，但许多分析库依赖于更适合计算和与医学图像交互的格式。我们使用NifTI （或 .nii 格式），最初是为脑成像开发的，但广泛用于DLTK 和本教程中的大多数其他体积图像。这和其他格式保存了重建图像容器并将其在物理空间中定向所需的信息。

为此，它需要专业的标题信息，我们将介绍一些需要考虑的用于深度学习的属性。

• 维度和大小存储有关如何重建图像的信息（例如将体积重建为三个维度，并使用大小向量）。
  
• 数据类型
  
• 体素间距（也称为体素的物理尺寸，通常以毫米为单位）
  
• 物理坐标系原点
  
• 方向

为什么这些属性很重要？ 网络将在体素空间中进行训练，这意味着我们将创建形状和维度为 [batch_size, dx, dy, dz, channels/features] 的张量，并将它们馈送到网络。网络将在该体素空间中进行训练，并假设所有图像（包括看不见的测试图像）都在该空间中进行归一化，或者可能存在泛化问题。在该体素空间中，特征提取器（例如卷积层）将假设体素维度是各向同性的（即在每个维度中都相同），并且所有图像都以相同的方式定向。

但是，由于大多数图像都描绘了物理空间，因此我们需要将该物理空间转换为常见的体素空间。

如果所有图像都以相同的方式定向（有时我们需要配准以在空间上对图像进行归一化：查看 MIRTK），我们可以通过以下方式计算从物理空间到体素空间的缩放变换：

phys_coords = origin + voxel_spacing * voxel_coord

其中所有这些信息都是存储在 .nii 标题中的向量。

读取 .nii 图像： 有几个库可以读取 .nii 文件并访问标题信息并对其进行解析以获得重建的图像容器作为numpy数组。我们选择了SimpleITK，它是ITK库的 Python 包装器，它允许我们导入额外的图像过滤器以进行预处理和其他任务。

import SimpleITK as sitk
import numpy as np

# A path to a T1-weighted brain .nii image:
t1_fn = './brain_t1_0001.nii'

# Read the .nii image containing the volume with SimpleITK:
sitk_t1 = sitk.ReadImage(t1_fn)

# and access the numpy array:
t1 = sitk.GetArrayFromImage(sitk_t1)

数据 I/O 注意事项

根据训练数据库的大小，有多种选择可以将 .nii 图像数据馈送到网络图。每种方法在速度方面都有特定的权衡，并且在训练过程中可能成为瓶颈。我们将介绍并解释三种选择。

内存中和馈送字典： 我们可以为网络图创建一个 tf.placeholder，并在训练期间通过 feed_dict 馈送它。我们从磁盘读取所有 .nii 文件，在 Python 中处理它们（参见load_data()），并将所有训练示例存储在内存中，并在其中进行馈送。

# Load all data into memory
data = load_data(all_filenames, tf.estimator.ModeKeys.TRAIN, reader_params)

# Create placeholder variables and define their shapes (here, 
# we input a volume image of size [128, 224, 244] and a single
# channel (i.e. greyscale):
x = tf.placeholder(reader_example_dtypes['features']['x'], 
                   [None, 128, 224, 224, 1])
y = tf.placeholder(reader_example_dtypes['labels']['y'], 
                   [None, 1])

# Create a tf.data.Dataset
dataset = tf.data.Dataset.from_tensor_slices((x, y))
dataset = dataset.repeat(None)
dataset = dataset.batch(batch_size)
dataset = dataset.prefetch(1)

# Create an iterator
iterator = dataset.make_initializable_iterator()
nx = iterator.get_next()

with tf.train.MonitoredTrainingSession() as sess_dict:
    
    sess_dict.run(iterator.initializer, 
               feed_dict={x: data['features'], y: data['labels']})
    
    for i in range(iterations):
        # Get next features/labels pair
        dict_batch_feat, dict_batch_lbl = sess_dict.run(nx)

TL;DR：这种直接方法通常是最快且最容易实现的方法，因为它避免了持续从磁盘读取数据，但是需要将整个训练示例数据库（以及验证示例）保存在内存中，对于更大的数据库或更大的图像文件来说这是不可行的。

使用TFRecords数据库： 对于大多数关于图像体积的深度学习问题，训练示例数据库太大，无法容纳在内存中。TFRecords格式允许序列化训练示例并将它们存储在磁盘上，并具有快速写入访问权限（即并行数据读取）。

def _int64_feature(value):
    return tf.train.Feature(int64_list=tf.train.Int64List(value=[value]))

def _float_feature(value):
    return tf.train.Feature(float_list=tf.train.FloatList(value=value))

# path to save the TFRecords file
train_filename = 'train.tfrecords'

# open the file
writer = tf.python_io.TFRecordWriter(train_filename)

# iterate through all .nii files:
for meta_data in all_filenames:

    # Load the image and label
    img, label = load_img(meta_data, reader_params)
    
    # Create a feature
    feature = {'train/label': _int64_feature(label),
               'train/image': _float_feature(img.ravel())}
               
    # Create an example protocol buffer
    example = tf.train.Example(features=tf.train.Features(feature=feature))
    
    # Serialize to string and write on the file
    writer.write(example.SerializeToString())
    
writer.close()

该格式可以与TensorFlow直接交互，并可以直接集成到 tf.graph 中的训练循环中。

def decode(serialized_example):
    # Decode examples stored in TFRecord
    # NOTE: make sure to specify the correct dimensions for the images
    features = tf.parse_single_example(
        serialized_example,
        features={'train/image': tf.FixedLenFeature([128, 224, 224, 1], tf.float32),
                  'train/label': tf.FixedLenFeature([], tf.int64)})

    # NOTE: No need to cast these features, as they are already `tf.float32` values.
    return features['train/image'], features['train/label']

dataset = tf.data.TFRecordDataset(train_filename).map(decode)
dataset = dataset.repeat(None)
dataset = dataset.batch(batch_size)
dataset = dataset.prefetch(1)

iterator = dataset.make_initializable_iterator()
features, labels = iterator.get_next()
nx = iterator.get_next()

with tf.train.MonitoredTrainingSession() as sess_rec:
    sess_rec.run(iterator.initializer)

    for i in range(iterations):
        try:
            # Get next features-labels pair
            rec_batch_feat, rec_batch_lbl = sess_rec.run([features, labels])

        except tf.errors.OutOfRangeError:
            pass

TL;DR：TFRecords是从磁盘访问文件的快速方法，但是需要存储整个训练数据库的另一个副本。如果我们的目标是使用数 TB 大小的数据库，这将是不可接受的。

使用本机 Python 生成器： 最后，我们可以使用 Python 生成器，创建一个read_fn()以直接加载图像数据……

def read_fn(file_references, mode, params=None):
    
    # We define a `read_fn` and iterate through the `file_references`, which
    # can contain information about the data to be read (e.g. a file path):
    for meta_data in file_references:
        
        # Here, we parse the `subject_id` to construct a file path to read
        # an image from.
        subject_id = meta_data[0]
        data_path = '../../data/IXI_HH/1mm'
        t1_fn = os.path.join(data_path, '{}/T1_1mm.nii.gz'.format(subject_id))
        
        # Read the .nii image containing a brain volume with SimpleITK and get 
        # the numpy array:
        sitk_t1 = sitk.ReadImage(t1_fn)
        t1 = sitk.GetArrayFromImage(sitk_t1)

        # Normalise the image to zero mean/unit std dev:
        t1 = whitening(t1)
        
        # Create a 4D Tensor with a dummy dimension for channels
        t1 = t1[..., np.newaxis]
        
        # If in PREDICT mode, yield the image (because there will be no label
        # present). Additionally, yield the sitk.Image pointer (including all
        # the header information) and some metadata (e.g. the subject id),
        # to facilitate post-processing (e.g. reslicing) and saving.
        # This can be useful when you want to use the same read function as 
        # python generator for deployment.
        if mode == tf.estimator.ModeKeys.PREDICT:
            yield {'features': {'x': t1}}

        # Labels: Here, we parse the class *sex* from the file_references 
        # \in [1,2] and shift them to \in [0,1] for training:
        sex = np.int32(meta_data[1]) - 1
        y = sex
        
        # If training should be done on image patches for improved mixing, 
        # memory limitations or class balancing, call a patch extractor
        if params['extract_examples']:
            images = extract_random_example_array(
                t1,
                example_size=params['example_size'],
                n_examples=params['n_examples'])
            
            # Loop the extracted image patches and yield
            for e in range(params['n_examples']):
                yield {'features': {'x': images[e].astype(np.float32)},
                       'labels': {'y': y.astype(np.int32)}}
                     
        # If desired (i.e. for evaluation, etc.), return the full images
        else:
            yield {'features': {'x': images},
                   'labels': {'y': y.astype(np.int32)}}

    return

以及tf.data.Dataset.from_generator() 来排队示例。

# Generator function
def f():
    fn = read_fn(file_references=all_filenames,
                 mode=tf.estimator.ModeKeys.TRAIN, 
                 params=reader_params)
    
    ex = next(fn)
    # Yield the next image
    yield ex
    
# Timed example with generator io
dataset = tf.data.Dataset.from_generator(
    f, reader_example_dtypes, reader_example_shapes)
dataset = dataset.repeat(None)
dataset = dataset.batch(batch_size)
dataset = dataset.prefetch(1)

iterator = dataset.make_initializable_iterator()
next_dict = iterator.get_next()

with tf.train.MonitoredTrainingSession() as sess_gen:
    # Initialize generator
    sess_gen.run(iterator.initializer)

    with Timer('Generator'):
        for i in range(iterations):
            # Fetch the next batch of images
            gen_batch_feat, gen_batch_lbl = sess_gen.run([next_dict['features'], next_dict['labels']])

TL;DR：它避免创建图像数据库的额外副本，但由于生成器无法并行读取和映射函数，因此比TFRecords慢得多。 速度基准测试和选择方法： 我们运行了这三种将 .nii 文件读入TensorFlow 的方法，并比较了加载和馈送固定大小的示例数据库所需的时间。所有代码和结果都可以在此处找到。

最明显最快的方法是通过占位符从内存中馈送，耗时 5.6 秒，其次是TFRecords，耗时 31.1 秒，而使用 Python 生成器从磁盘读取未经优化的速度为 123.5 秒。但是，只要训练过程中的正向/反向传递是计算瓶颈，数据 I/O 的速度就可以忽略不计。

数据归一化

与自然图像一样，我们可以对生物医学图像数据进行归一化，但方法可能略有不同。归一化的目的是消除数据中的一些变化（例如不同的受试者姿势或图像对比度的差异等），这些变化是已知的，因此简化了我们感兴趣的细微差异的检测（例如病理的存在）。在这里，我们将介绍最常见的归一化形式。

体素强度的归一化： 这种形式高度依赖于获取数据的成像模式。典型的零均值，单位方差归一化是定性图像（例如加权脑部 MR 图像，其中对比度高度依赖于采集参数，通常由专家设置）的标准。如果我们使用这种统计方法，我们使用来自完整单一体积的统计信息，而不是整个数据库。与之相反，定量成像测量物理量（例如 CT 成像中的射线密度，其中强度在不同扫描仪之间具有可比性），并受益于裁剪和/或重新缩放，例如 简单范围归一化（例如到 [-1,1]）。

强度归一化方法示例

空间归一化： 对图像方向进行归一化可以避免模型不得不学习所有可能的定向，这在很大程度上减少了所需的训练图像数量（见标题属性的重要性，可以了解图像的定向）。我们还会考虑体素间距，即使是在从同一台扫描仪获取的图像之间，体素间距也可能有所不同。可以通过重新采样到各向同性分辨率来实现这一点。

def resample_img(itk_image, out_spacing=[2.0, 2.0, 2.0], is_label=False):
    
    # Resample images to 2mm spacing with SimpleITK
    original_spacing = itk_image.GetSpacing()
    original_size = itk_image.GetSize()

    out_size = [
        int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),
        int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),
        int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))]

    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size)
    resample.SetOutputDirection(itk_image.GetDirection())
    resample.SetOutputOrigin(itk_image.GetOrigin())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(itk_image.GetPixelIDValue())

    if is_label:
        resample.SetInterpolator(sitk.sitkNearestNeighbor)
    else:
        resample.SetInterpolator(sitk.sitkBSpline)

    return resample.Execute(itk_image)

# Assume to have some sitk image (itk_image) and label (itk_label)
resampled_sitk_img = resample_img(itk_image, out_spacing=[2.0, 2.0, 2.0], is_label=False)
resampled_sitk_lbl = resample_img(itk_label, out_spacing=[2.0, 2.0, 2.0], is_label=True)

如果需要进一步标准化，我们可以使用医学图像配准软件包（例如 *MIRTK* 等），并将图像配准到相同的空间，以便图像之间的体素位置相互对应。分析结构性脑部 MR 图像（例如 T1 加权 MR 图像）的典型步骤是将训练数据库中的所有图像配准到参考标准，例如平均图谱（例如 *MNI 305* 图谱）。根据配准方法的自由度，这也可以对大小（仿射配准）或形状（可变形配准）进行标准化。这两种变体很少使用，因为它们会删除图像中的一些信息（即形状信息或大小信息），而这些信息可能对分析很重要（例如，心脏较大可能预示着心脏病）。

数据增强

通常，可用数据的数量有限，并且某些变异没有被涵盖。以下是一些示例：

• 软组织器官，存在多种正常的形状
  
• 病理，例如癌症病灶，其形状和位置差异很大
  
• 自由手超声图像，可能存在大量可能的视角
  

为了对未见测试用例进行适当的泛化，我们通过模拟我们希望对其保持稳健的数据变异来增强训练图像。与标准化方法类似，我们区分强度增强和空间增强

强度增强的示例

• 向训练图像添加噪声可以泛化到噪声图像
  
• 添加随机偏移或对比度以处理图像之间的差异
  

空间增强的示例

• 在预期对称的方向上翻转图像张量（例如，对脑部扫描进行左右翻转）
  
• 随机变形（例如，用于模拟器官形状的差异）
  
• 沿轴旋转（例如，用于模拟不同的超声视图角度）
  
• 随机裁剪和对补丁进行训练
  

强度和空间增强技术的示例

有关增强和数据 I/O 的重要说明：根据需要或有帮助的增强，某些操作仅在 python 中可用（例如，随机变形），这意味着如果使用使用原始 *TensorFlow*（即 *TFRecords* 或 *tf.placeholder*）的读取方法，则需要预先计算并存储到磁盘，从而大大增加训练数据库的大小。

类别平衡

领域专家解释（例如，手动分割或疾病类别）是监督学习医学图像的必要条件。通常，图像级（例如，疾病类别）或体素级（即分割）标签的比例并不相同，这意味着网络在训练期间不会看到来自每个类别的相同数量的示例。如果类别比率大致相似（例如，对于二元分类情况，30/70），这不会对准确性产生很大影响。但是，由于大多数损失是整个批次的平均成本，因此网络将首先学习正确预测最常出现的类别（例如，背景或正常情况，通常有更多可用的示例）。

训练期间的类别不平衡将对罕见现象（例如，图像分割中的小病灶）产生更大的影响，并严重影响测试准确性。

为了避免这种下降，有两种典型的方法可以解决数据集中类别不平衡的问题

• 通过采样进行类别平衡：在这里，我们的目标是在采样期间纠正看到的示例的频率。这可以通过以下方式实现：a）从每个类别中采样相同数量的样本，b）对过度表示的类别进行欠采样，或 c）对频率较低的类别进行过采样。在 *DLTK* 中，我们有对 a）的实现，可以找到 这里。我们对图像体积中的随机位置进行采样，如果提取的示例包含我们正在查找的类别，则考虑该示例。
  
• 通过损失函数进行类别平衡：与典型的体素级平均损失（例如，分类交叉熵、L2 等）相比，我们可以 a）使用本质上是平衡的损失函数（例如，平滑骰子损失，它是所有类别的平均骰子系数）或 b）通过类别频率重新加权每个预测的损失（例如，中位数频率重新加权交叉熵）。

示例应用程序亮点

有了这篇博文中提供的所有基本知识，我们现在可以研究使用 *TensorFlow* 构建医学图像深度学习的完整应用程序。我们已经使用深度神经网络实现了几个典型的应用程序，并将介绍其中的一些应用程序，让您了解现在可以尝试解决哪些问题。

注意：这些示例应用程序学习了一些有意义的东西，但它们是为了演示目的而构建的，而不是为了高性能实现。

示例数据集

我们提供以下所有示例的 下载和预处理脚本。对于大多数情况（包括上面的演示），我们使用了 *IXI 脑部数据库*。对于图像分割，我们下载了 *MRBrainS13 挑战数据库*，您需要注册才能下载。

多通道脑部 MR 图像的图像分割

多序列图像输入、目标标签和预测的 Tensorboard 可视化

此图像分割应用程序学习从小型（N=5）*MRBrainS* 挑战数据集的多序列 MR 图像（T1 加权、T1 反转恢复和 T2 FLAIR）中预测脑组织和白质病变。它使用具有残差单元作为特征提取器的 3D U-Net 类网络，并在 *TensorBoard* 中跟踪每个标签的骰子系数准确性。

代码和说明可以找到 这里。

从 T1 加权脑部 MR 图像中进行年龄回归和性别分类

回归和分类的示例输入 T1 加权脑部 MR 图像

两个类似的应用程序采用可扩展的 3D ResNet 架构，从 *IXI* 数据库的 T1 加权脑部 MR 图像中学习预测受试者的年龄（回归）或受试者的性别（分类）。这些应用程序之间的主要区别在于损失函数：当我们训练回归网络以使用 L2 损失预测年龄作为连续变量（预测年龄和实际年龄之间的均方差）时，我们使用分类交叉熵损失来预测性别的类别。

这些应用程序的代码和说明可以在这里找到：分类，回归。

对 3T 多通道脑部 MR 图像进行表征学习

使用深度卷积自动编码器网络的测试图像和重建

在这里，我们演示了深度卷积自动编码器架构的使用，这是一种强大的表征学习工具：网络以多序列 MR 图像作为输入，并试图重建它们。通过这样做，它将整个训练数据库的信息压缩到它的潜在变量中。经过训练的权重也可以用于迁移学习或信息压缩。请注意，重建的图像非常平滑：这可能是由于这个应用程序使用了 L2 损失函数，或者网络太小而无法正确编码详细的信息。

代码和说明可以找到 这里。

T1w 脑部 MR 图像的简单图像超分辨率

图像超分辨率：原始目标图像；降采样输入图像；线性上采样图像；预测的超分辨率；

单图像超分辨率旨在学习如何从低分辨率输入中上采样和重建高分辨率图像。这个简单的实现创建了一个图像的低分辨率版本，超分辨率网络学习将图像上采样到其原始分辨率（这里上采样因子为 [4,4,4]）。此外，我们计算了一个线性上采样版本，以显示与重建图像的差异。

代码和说明可以找到 这里。

最后…

我们希望本教程能帮助您轻松了解生物医学图像的深度学习。如果您觉得它有用，我们感谢您分享它并在 *github* 上关注 *DLTK*。如果您需要类似问题的帮助，请访问我们的 *gitter.io 聊天* 并向我们提问。也许有一天，我们可以将您的应用程序托管在 *DLTK* 的 模型库 中。感谢您的阅读！

https://twitter.com/dltk_; https://dltk.github.io; https://gitter.im/DLTK/DLTK;

资源

教程代码，示例应用程序，DLTK 源代码