第一课第三周大作业-- MRI 脑肿瘤自动分割教程
Posted Tina姐
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了第一课第三周大作业-- MRI 脑肿瘤自动分割教程相关的知识，希望对你有一定的参考价值。
文章目录
作业文件
吴恩达-医学图像AI 专项课程-作业/第一课/第一课大作业/week3seg/C1M3_Assignmentlast.ipynb
欢迎来到“医学人工智能”课程 1 的最后一部分！
您将学习如何使用 MRI数据构建神经网络以自动分割大脑中的肿瘤区域。
你将学习：
- MR 图像中有什么
- MRI 数据集的标准数据准备技术
- 分割的度量和损失函数
- 可视化和评估分割模型
作业目录
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 表示完美匹配 。
DSC
(
A
,
B
)
=
2
×
∣
A
∩
B
∣
∣
A
∣
+
∣
B
∣
.
\\text{DSC}(A, B) = \\frac{2 \\times |A \\cap B|}{|A| + |B|}.
DSC(A,B)=∣A∣+∣B∣2×∣A∩B∣.
这里我们可以解读 𝐴 和 𝐵 作为体素集，𝐴 作为预测的肿瘤区域和 𝐵 是 ground truth。
我们的模型会将每个体素映射到 0 或 1
- 0 表示它是背景体素
- 1 表示它是分割区域的一部分。
DSC ( f , x , y ) = 2 × ∑ i , j f ( x ) i j × y i j + ϵ ∑ i , j f ( x ) i j + ∑ i , j y i j + ϵ \\text{DSC}(f, x, y) = \\frac{2 \\times \\sum_{i, j} f(x)_{ij} \\times y_{ij} + \\epsilon}{\\sum_{i,j} f(x)_{ij} + \\sum_{i, j} y_{ij} + \\epsilon} DSC(f,x,y)=∑i,jf(x)ij+∑i,jyij+ϵ2×∑i,jf(x)ij×yij+ϵ
- x x x : the input image
- f ( x ) f(x) f(x) : the model output (prediction)
- y y y : the label (actual ground truth)
-
ϵ
\\epsilon
ϵ : 是为了避免被零除而添加的小数
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
L D i c e ( p , q ) = 1 − 2 × ∑ i , j p i j q i j + ϵ ( ∑ i , j p i j 2 ) + ( ∑ i , j q i j 2 ) + ϵ \\mathcal{L}_{Dice}(p, q) = 1 - \\frac{2\\times\\sum_{i, j} p_{ij}q_{ij} + \\epsilon}{\\left(\\sum_{i, j} p_{ij}^2 \\right) + \\left(\\sum_{i, j} q_{ij}^2 \\right) + \\epsilon} LDice(p,q)=1−(∑i,jpij2)+(∑i,jqij2)+ϵ2×∑i,jpijqij+ϵ
- p p p is our predictions
- q q q is the ground truth
- q i q_i qi 是 0 or 1.
soft dice loss 就是 1-DSC
如果是多类别，就把每个类别的 DSC 累加起来：
L
D
i
c
e
(
p
,
q
)
=
1
−
1
N
∑
c
=
1
C
2
×
∑
i
,
j
p
c
i
j
q
c
i
j
+
ϵ
(
∑
i
,
j
p
c
i
j
2
)
+
(
∑
i
,
j
q
c
i
j
2
)
+
ϵ
\\mathcal{L}_{Dice}(p, q) = 1 - \\frac{1}{N} \\sum_{c=1}^{C} \\frac{2\\times\\sum_{i, j} p_{cij}q_{cij} + \\epsilon}{\\left(\\sum_{i, j} p_{cij}^2 \\right) + \\left(\\sum_{i, j} q_{cij}^2 \\right) + \\epsilon}
LDice(p,q)=1−N1c=1∑C(∑i,j第一课第三周3-4节--2D和3D分割的优缺点以及如何构建UNet模型