[人工智能]第一课第三周大作业-- MRI 脑肿瘤自动分割教程

吴恩达-医学图像AI 专项课程-作业/第一课/第一课大作业/week3seg/C1M3_Assignmentlast.ipynb

欢迎来到“医学人工智能”课程 1 的最后一部分！

您将学习如何使用 MRI数据构建神经网络以自动分割大脑中的肿瘤区域。

你将学习：

• MR 图像中有什么
• MRI 数据集的标准数据准备技术
• 分割的度量和损失函数
• 可视化和评估分割模型

作业目录

• 1. Dataset
  • 1.1 What is an MRI?
  • 1.2 MRI Data Processing
  • 1.3 Exploring the Dataset
  • 1.4 Data Preprocessing
    • 1.4.1 Sub-volume Sampling
    • 1.4.2 Standardization
• 2. Model: 3D U-Net
• 3. Metrics
  • 3.1 Dice Coefficient
  • 3.2 Soft Dice Loss
• 4. Training
• 5. Evaluation
  • 5.1 Overall Performance
  • 5.2 Patch-level Predictions
  • 5.3 Running on Entire Scans

2.1 MRI data

在这个作业中，我们将构建一个多类分割模型。我们将在每个图像中识别 3 种不同的异常：水肿、非增强肿瘤和增强肿瘤。

可以使用ITK-SNAP打开来看一下图像和标签

使用的MRI数据来自

Medical Segmentation Decathlon

这是一个极其丰富的数据集，提供 3D 表示中每个点（体素）相关联的标签。

2.2 探索数据

我们的数据集以NifTI-1 格式存储，我们将使用 NiBabel库与文件进行交互。每个训练样本由两个单独的文件组成：

第一个文件是一个图像文件，包含形状为 (240, 240, 155, 4) 的 MR 图像的 4D 阵列。

• 前 3 个维度是 3D 体积中每个点的 X、Y 和 Z 值，通常称为体素。
• 第 4 维是 4 个不同序列的值
  • 0：FLAIR序列：“流体衰减反转恢复”（FLAIR）
  • 1：T1w：“T1 加权”
  • 2：t1gd：“T1 加权钆对比度增强”（T1-Gd）
  • 3：T2w：“T2 加权”

每个训练示例中的第二个文件是一个标签文件，其中包含一个形状为 (240, 240, 155) 的 3D 数组。

• 此数组中的整数值表示相应图像文件中每个体素的“标签”：
  • 0：背景
  • 1：水肿
  • 2：非强化肿瘤
  • 3：增强肿瘤

我们可以访问总共 484 张训练图像，我们将这些图像分成训练 (80%) 和验证 (20%) 数据集。

让我们从查看一个案例并可视化数据开始！您可以通过此笔记本访问 10 个不同的案例，我们强烈建议您自行进一步探索数据。

分别按横断位、冠状位、矢状位显示图像，并且把标签显示在图像上

image, label = load_case(DATA_DIR + "imagesTr/BRATS_003.nii.gz", DATA_DIR + "labelsTr/BRATS_003.nii.gz")
image = util.get_labeled_image(image, label)

util.plot_image_grid(image)

我们还编写了一个实用函数，该函数生成一个GIF，显示在每个轴上迭代的效果

image, label = load_case(DATA_DIR + "imagesTr/BRATS_003.nii.gz", DATA_DIR + "labelsTr/BRATS_003.nii.gz")
util.visualize_data_gif(util.get_labeled_image(image, label))

2.4 使用 patch 进行数据预处理

虽然我们的数据集是在配准后以 NIfTI 格式提供给我们的，但在将数据提供给我们的模型之前，我们仍然需要做一些小的预处理。

生成子体积（sub-volumes）

我们将首先生成数据的“patch 补丁”，您可以将其视为整个 MR 图像的子体积。

我们生成补丁的原因是因为原始图像太大，容易超出内存。

因此，我们将使用这种通用技术来生成空间一致的数据子体积，这些子体积可以输入到我们的网络中。

具体来说，我们将从我们的图像中生成形状为 [160, 160, 16] 的随机采样子体积。

此外，由于随机采样子体积，可能会导致大部分子体积只是脑组织或黑色背景，没有任何肿瘤，我们希望生成的子体积中至少包含一定数量肿瘤体素。

因此，我们将只选择具有最多 95% 非肿瘤区域（或至少 5% 肿瘤）的补丁（这是生成的子体积中，肿瘤和非肿瘤体素的占比）

我们通过根据背景标签中存在的值过滤来做到这一点。

def get_sub_volume(image, label, 
                   orig_x = 240, orig_y = 240, orig_z = 155, 
                   output_x = 160, output_y = 160, output_z = 16,
                   num_classes = 4, max_tries = 1000, 
                   background_threshold=0.95):
    """
    Extract random sub-volume from original images.

    Args:
        image (np.array): original image, 
            of shape (orig_x, orig_y, orig_z, num_channels)
        label (np.array): original label. 
            labels coded using discrete values rather than
            a separate dimension, 
            so this is of shape (orig_x, orig_y, orig_z)
        orig_x (int): x_dim of input image
        orig_y (int): y_dim of input image
        orig_z (int): z_dim of input image
        output_x (int): desired x_dim of output
        output_y (int): desired y_dim of output
        output_z (int): desired z_dim of output
        num_classes (int): number of class labels
        max_tries (int): maximum trials to do when sampling
        background_threshold (float): limit on the fraction 
            of the sample which can be the background

    returns:
        X (np.array): sample of original image of dimension 
            (num_channels, output_x, output_y, output_z)
        y (np.array): labels which correspond to X, of dimension 
            (num_classes, output_x, output_y, output_z)
    """
    # Initialize features and labels with `None`
    X = None
    y = None

    ### START CODE HERE (REPLACE INSTANCES OF 'None' with your code) ###
    
    tries = 0
    
    while tries < max_tries:
        # randomly sample sub-volume by sampling the corner voxel
        # hint: make sure to leave enough room for the output dimensions!
        start_x = np.random.randint(0,orig_x - output_x + 1)
        start_y = np.random.randint(0,orig_y - output_y + 1)
        start_z = np.random.randint(0,orig_z - output_z + 1)

        # extract relevant area of label
        y = label[start_x: start_x + output_x,
                  start_y: start_y + output_y,
                  start_z: start_z + output_z]
        
        # One-hot encode the categories.
        # This adds a 4th dimension, 'num_classes'
        # (output_x, output_y, output_z, num_classes)
        y = keras.utils.to_categorical(y,num_classes= num_classes)


        # compute the background ratio
        bgrd_ratio = np.sum(y[:, :, :, 0])/(output_x * output_y * output_z)

        # increment tries counter
        tries += 1

        # if background ratio is below the desired threshold,
        # use that sub-volume.
        # otherwise continue the loop and try another random sub-volume
        if bgrd_ratio < background_threshold:

            # make copy of the sub-volume
            X = np.copy(image[start_x: start_x + output_x,
                              start_y: start_y + output_y,
                              start_z: start_z + output_z, :])
            
            # change dimension of X
            # from (x_dim, y_dim, z_dim, num_channels)
            # to (num_channels, x_dim, y_dim, z_dim)
            X = np.moveaxis(X,3,0)

            # change dimension of y
            # from (x_dim, y_dim, z_dim, num_classes)
            # to (num_classes, x_dim, y_dim, z_dim)
            y = np.moveaxis(y,3,0)


            ### END CODE HERE ###
            
            # take a subset of y that excludes the background class
            # in the 'num_classes' dimension
            y = y[1:, :, :, :]
    
            return X, y

    # if we've tried max_tries number of samples
    # Give up in order to avoid looping forever.
    print(f"Tried {tries} times to find a sub-volume. Giving up...")

标准化

最后，鉴于 MR 图像中的值涵盖的范围非常广泛，我们将这些值标准化为均值为 0，标准差为 1。

这是深度图像处理中的常用技术，因为标准化使网络更容易学习。

# UNQ_C2 (UNIQUE CELL IDENTIFIER, DO NOT EDIT)
def standardize(image):
    """
    Standardize mean and standard deviation 
        of each channel and z_dimension.

    Args:
        image (np.array): input image, 
            shape (num_channels, dim_x, dim_y, dim_z)

    Returns:
        standardized_image (np.array): standardized version of input image
    """
    
    ### START CODE HERE (REPLACE INSTANCES OF 'None' with your code) ###
    
    # initialize to array of zeros, with same shape as the image
    standardized_image = np.zeros(image.shape)

    # iterate over channels
    for c in range(image.shape[0]):
        # iterate over the `z` dimension
        for z in range(image.shape[3]):
            # get a slice of the image 
            # at channel c and z-th dimension `z`
            image_slice = image[c,:,:,z]

            # subtract the mean from image_slice
            mean=np.mean(image_slice)
            centered = image_slice-mean
            
            # divide by the standard deviation (only if it is different from zero)
            std=np.std(centered)
            if std != 0:
                centered_scaled = centered/std

            # update  the slice of standardized image
            # with the scaled centered and scaled image
            standardized_image[c, :, :, z] = centered_scaled

    ### END CODE HERE ###

    return standardized_image

3D U-Net

模型使用的是3D unet

4.1 骰子相似系数

除了架构之外，任何深度学习方法中最重要的元素之一就是损失函数的选择。

您可能熟悉的一个自然选择是交叉熵损失函数。

然而，由于严重的类别不平衡（通常没有很多正区域），这种损失函数对于分割任务并不理想。

分割任务的一个更常见的损失是 Dice 相似系数，它是衡量两个轮廓重叠程度的指标。

Dice = 0 表示完全不匹， Dice =1 表示完美匹配 。

$$\text{DSC}(A,B)=\frac{2\times |A\cap B|}{|A|+|B|}.$$
这里我们可以解读 𝐴 和 𝐵 作为体素集，𝐴 作为预测的肿瘤区域和 𝐵 是 ground truth。

我们的模型会将每个体素映射到 0 或 1

• 0 表示它是背景体素
• 1 表示它是分割区域的一部分。

$$\text{DSC}(f,x,y)=\frac{2\times {\sum }_{i,j}f(x{)}_{ij}\times y_{ij}+?}{{\sum }_{i,j}f(x{)}_{ij}+{\sum }_{i,j}{y}_{ij}+?}$$

• $x$ : the input image
• $f(x)$ : the model output (prediction)
• $y$ : the label (actual ground truth)
• $?$ : 是为了避免被零除而添加的小数
  

def single_class_dice_coefficient(y_true, y_pred, axis=(0, 1, 2), 
                                  epsilon=0.00001):
    ### START CODE HERE (REPLACE INSTANCES OF 'None' with your code) ###
    
    dice_numerator = 2 * K.sum(y_true * y_pred) + epsilon
    dice_denominator = K.sum(y_true) + K.sum(y_pred) + epsilon
    dice_coefficient = dice_numerator / dice_denominator
    return dice_coefficient

4.2 soft dice loss

$${\mathcal{L}}_{Dice}(p,q)=1?\frac{2\times {\sum }_{i,j}{p}_{ij}{q}_{ij}+?}{\left({\sum }_{i,j}{p}_{ij}^2\right)+\left({\sum }_{i,j}{q}_{ij}^2\right)+?}$$

• $p$ is our predictions
• $q$ is the ground truth
• $q_i$ 是 0 or 1.

soft dice loss 就是 1-DSC

如果是多类别，就把每个类别的 DSC 累加起来：
$${\mathcal{L}}_{Dice}(p,q)=1?\frac{1}{N}\sum _{c=1}^C\frac{2\times {\sum }_{i,j}{p}_{cij}{q}_{cij}+?}{\left({\sum }_{i,j}{p}_{cij}^2\right)+\left({\sum }_{i,j}{q}_{cij}^2\right)+?}$$

def soft_dice_loss(y_true, y_pred, axis=(1, 2, 3), 
                   epsilon=0.00001):
    dice_numerator = 2. * K.sum(y_true * y_pred, axis = axis) + epsilon
    dice_denominator = K.sum(y_true**2, axis = axis) + K.sum(y_pred**2, axis = axis) + epsilon
    dice_loss = 1 - K.mean(dice_numerator / dice_denominator)
    return dice_loss

5 创建模型

model = util.unet_model_3d(loss_function=soft_dice_loss, metrics=[dice_coefficient])

模型在 utils.py 文件中提供

5.1 训练

如果你想训练模型，将5.1 的代码运行一遍

6 评估

训练好的模型，没有给到大家，你可以自己训练一个模型，再评估。或者直接阅读以下评估方法，不用实践。

对于分割模型，我们同样可以计算 敏感度，特异性

还可以可视化预测结果

文章持续更新，可以关注微信公众号【医学图像人工智能实战营】获取最新动态，一个关注于医学图像处理领域前沿科技的公众号。坚持已实践为主，手把手带你做项目，打比赛，写论文。凡原创文章皆提供理论讲解，实验代码，实验数据。只有实践才能成长的更快，关注我们，一起学习进步~

我是Tina, 我们下篇博客见~

白天工作晚上写文，呕心沥血

觉得写的不错的话最后，求点赞，评论，收藏。或者一键三连