Hello everyone, I am “K student ah”!

It seems that I have not updated for some time. There are too many things in this period, and I am a little lazy, but I still insist on it. When the school starts, the update frequency may be stable.

The end of the conversation, into the topic, some time ago to help others to do a 3D classification of the live, thinking how also have to sort out a blog for everyone, so there is this article. The project described in this article has the following improvements over earlier projects:

  • 1. Set the dynamic learning rate
  • 2. Added early stop policy
  • 3. Save time of model is more “intelligent”
  • 4. There are also obvious optimizations in data loading

🚀 My environment:

  • Language: Python3.6.5
  • Compiler: Jupyter Notebook
  • Deep learning environment: TensorFlow2.4.1
  • Graphics card: NVIDIA GeForce RTX 3080
  • Data link: pan.baidu.com/s/1_8T3PthA… Extraction code: UEKA

🚀 from the column:100 Examples of Deep Learning

If you are an expert in deep learning, you can take a look at the column I wrote for you: Introduction to Deep Learning.

CT scan 3D image classification

  • Creation time: The afternoon of August 23, 2021
  • Task description: Train a 3D convolutional neural network to predict the presence of pneumonia.

I. Preliminary work

1. Introduction

This case will demonstrate whether viral pneumonia can be predicted in COMPUTED tomography (CT) by constructing 3D convolutional neural network (CNN). 2D CNN is usually used to process RGB images (3 channels). The CNN of 3D is only the 3D equivalent, and we can simply understand the 3D image as the superposition of 2D image. 3D CNN can be understood as a powerful model for learning stereoscopic data.

import os,zipfile
import numpy as np
from tensorflow import keras
from tensorflow.keras import layers

import tensorflow as tf
gpus = tf.config.list_physical_devices("GPU")

if gpus:
    tf.config.experimental.set_memory_growth(gpus[0].True)  Set GPU memory usage as required
    tf.config.set_visible_devices([gpus[0]],"GPU")
    
Print the graphics card information and confirm that the GPU is available
print(gpus)
Copy the code
[PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')]
Copy the code

2. Load and preprocess data

The data file is Nifti with the.nii extension. I use the nibabel package to read the files, which you can install by PIP install nibabel.

Data preprocessing steps:

  1. First turn the volume 90 degrees to make sure the orientation is fixed
  2. Scale the HU value between 0 and 1.
  3. Adjust width, height, and depth.

I have defined several helper functions to handle the data, which will be used when building training and validation datasets.

import nibabel as nib
from scipy import ndimage

def read_nifti_file(filepath) :
    # read file
    scan = nib.load(filepath)
    # Fetch data
    scan = scan.get_fdata()
    return scan

def normalize(volume) :
    """ "normalization """
    min = -1000
    max = 400
    volume[volume < min] = min
    volume[volume > max] = max
    volume = (volume - min)/(max - min)
    volume = volume.astype("float32")
    return volume

def resize_volume(img) :
    """ Modify image size """
    # Set the desired depth
    desired_depth = 64
    desired_width = 128
    desired_height = 128
    # Get current depth
    current_depth = img.shape[-1]
    current_width = img.shape[0]
    current_height = img.shape[1]
    # Compute depth factor
    depth = current_depth / desired_depth
    width = current_width / desired_width
    height = current_height / desired_height
    depth_factor = 1 / depth
    width_factor = 1 / width
    height_factor = 1 / height
    # rotation
    img = ndimage.rotate(img, 90, reshape=False)
    # Data adjustment
    img = ndimage.zoom(img, (width_factor, height_factor, depth_factor), order=1)
    return img

def process_scan(path) :
    # read file
    volume = read_nifti_file(path)
    # normalization
    volume = normalize(volume)
    Adjust width, height and depth
    volume = resize_volume(volume)
    return volume
Copy the code

Path to read the CT scan file

# "CT-0" folder is a CT scan of normal lung tissue
normal_scan_paths = [
    os.path.join(os.getcwd(), "MosMedData/CT-0", x)
    for x in os.listdir("MosMedData/CT-0")]# "CT-23" folder contains CT scans of people with pneumonia
abnormal_scan_paths = [
    os.path.join(os.getcwd(), "MosMedData/CT-23", x)
    for x in os.listdir("MosMedData/CT-23")]print("CT scans with normal lung tissue: " + str(len(normal_scan_paths)))
print("CT scans with abnormal lung tissue: " + str(len(abnormal_scan_paths)))
Copy the code
CT scans with normal lung tissue: 100
CT scans with abnormal lung tissue: 100
Copy the code
Read the data and preprocess it
abnormal_scans = np.array([process_scan(path) for path in abnormal_scan_paths])
normal_scans = np.array([process_scan(path) for path in normal_scan_paths])

# Tag digitization
abnormal_labels = np.array([1 for _ in range(len(abnormal_scans))])
normal_labels = np.array([0 for _ in range(len(normal_scans))])
Copy the code

Second, build training and validation sets

Read scan and assign labels from the class directory. The scan is downsampled to have a 128x128x64 shape. Readjust the original HU value to a range of 0 to 1. Finally, the data set is split into training and validation subsets.

Divide the training set and verification set in a 7:3 ratio
x_train = np.concatenate((abnormal_scans[:70], normal_scans[:70]), axis=0)
y_train = np.concatenate((abnormal_labels[:70], normal_labels[:70]), axis=0)
x_val = np.concatenate((abnormal_scans[70:], normal_scans[70:]), axis=0)
y_val = np.concatenate((abnormal_labels[70:], normal_labels[70:]), axis=0)
print(
    "Number of samples in train and validation are %d and %d."
    % (x_train.shape[0], x_val.shape[0]))Copy the code
Number of samples in train and validation are 140 and 60.
Copy the code

Third, data enhancement

CT scans also enhance the data by rotating at random angles during training. Since the data is stored in the shape of Rank-3 (sample, height, width, depth), we add a dimension of size 1 at axis 4 to enable 3D convolution of the data. Therefore, the new shape (sample, height, width, depth, 1). There are different types of preprocessing and enhancement techniques out there, and this example shows some simple beginnings.

import random
from scipy import ndimage

@tf.function
def rotate(volume) :
    """ Different degrees of rotation. """
    def scipy_rotate(volume) :
        Define some rotation angles
        angles = [-20, -10, -5.5.10.20]
        # Select an Angle at random
        angle = random.choice(angles)

        volume = ndimage.rotate(volume, angle, reshape=False)
        volume[volume < 0] = 0
        volume[volume > 1] = 1
        return volume

    augmented_volume = tf.numpy_function(scipy_rotate, [volume], tf.float32)
    return augmented_volume

def train_preprocessing(volume, label) :
    volume = rotate(volume)
    volume = tf.expand_dims(volume, axis=3)
    return volume, label

def validation_preprocessing(volume, label) :
    volume = tf.expand_dims(volume, axis=3)
    return volume, label
Copy the code

While defining the training and validation data loaders, the training data will be randomly rotated at different angles. Both training and validation data have been readjusted to have values between 0 and 1.

Define the data loader
train_loader = tf.data.Dataset.from_tensor_slices((x_train, y_train))
validation_loader = tf.data.Dataset.from_tensor_slices((x_val, y_val))

batch_size = 2

train_dataset = (
    train_loader.shuffle(len(x_train))
    .map(train_preprocessing)
    .batch(batch_size)
    .prefetch(2)
)

validation_dataset = (
    validation_loader.shuffle(len(x_val))
    .map(validation_preprocessing)
    .batch(batch_size)
    .prefetch(2))Copy the code

4. Data visualization

import matplotlib.pyplot as plt

data = train_dataset.take(1)
images, labels = list(data)[0]
images = images.numpy()
image = images[0]
print("Dimension of the CT scan is:", image.shape)
plt.imshow(np.squeeze(image[:, :, 30]), cmap="gray")
Copy the code
Dimension of the CT scan is: (128, 128, 64, 1)
Copy the code

def plot_slices(num_rows, num_columns, width, height, data) :
    """Plot a montage of 20 CT slices"""
    data = np.rot90(np.array(data))
    data = np.transpose(data)
    data = np.reshape(data, (num_rows, num_columns, width, height))
    rows_data, columns_data = data.shape[0], data.shape[1]
    heights = [slc[0].shape[0] for slc in data]
    widths = [slc.shape[1] for slc in data[0]]
    fig_width = 12.0
    fig_height = fig_width * sum(heights) / sum(widths)
    f, axarr = plt.subplots(
        rows_data,
        columns_data,
        figsize=(fig_width, fig_height),
        gridspec_kw={"height_ratios": heights},
    )
    for i in range(rows_data):
        for j in range(columns_data):
            axarr[i, j].imshow(data[i][j], cmap="gray")
            axarr[i, j].axis("off")
    plt.subplots_adjust(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1)
    plt.show()

# Visualize montage of slices.
# 4 rows and 10 columns for 100 slices of the CT scan.
plot_slices(4.10.128.128, image[:, :, :40])
Copy the code

5. Build 3D convolutional neural network model

To make the model easier to understand, I built it into blocks.

def get_model(width=128, height=128, depth=64) :
    Construction of 3D Convolutional Neural Network Model

    inputs = keras.Input((width, height, depth, 1))

    x = layers.Conv3D(filters=64, kernel_size=3, activation="relu")(inputs)
    x = layers.MaxPool3D(pool_size=2)(x)
    x = layers.BatchNormalization()(x)

    x = layers.Conv3D(filters=64, kernel_size=3, activation="relu")(x)
    x = layers.MaxPool3D(pool_size=2)(x)
    x = layers.BatchNormalization()(x)

    x = layers.Conv3D(filters=128, kernel_size=3, activation="relu")(x)
    x = layers.MaxPool3D(pool_size=2)(x)
    x = layers.BatchNormalization()(x)

    x = layers.Conv3D(filters=256, kernel_size=3, activation="relu")(x)
    x = layers.MaxPool3D(pool_size=2)(x)
    x = layers.BatchNormalization()(x)

    x = layers.GlobalAveragePooling3D()(x)
    x = layers.Dense(units=512, activation="relu")(x)
    x = layers.Dropout(0.3)(x)

    outputs = layers.Dense(units=1, activation="sigmoid")(x)

    # Define model
    model = keras.Model(inputs, outputs, name="3dcnn")
    return model

# Build a model
model = get_model(width=128, height=128, depth=64)
model.summary()
Copy the code
Model: "3dcnn" _________________________________________________________________ Layer (type) Output Shape Param # ================================================================= input_1 (InputLayer) [(None, 128, 128, 64, 1)] 0 _________________________________________________________________ conv3d (Conv3D) (None, 126, 126, 62, 64) 1792 _________________________________________________________________ max_pooling3d (MaxPooling3D) (None, 63, 63, 31, 64) 0 _________________________________________________________________ batch_normalization (BatchNo (None, 63, 63, 31, 64) 256 _________________________________________________________________ conv3d_1 (Conv3D) (None, 61, 61, 29, 64) 110656 _________________________________________________________________ max_pooling3d_1 (MaxPooling3 (None, 30, 30, 14, 64) 0 _________________________________________________________________ batch_normalization_1 (Batch (None, 30, 30, 14, 64) 256 _________________________________________________________________ conv3d_2 (Conv3D) (None, 28, 28, 12, 128) 221312 _________________________________________________________________ max_pooling3d_2 (MaxPooling3 (None, 14, 14, 6, 128) 0 _________________________________________________________________ batch_normalization_2 (Batch (None, 14, 14, 6, 128) 512 _________________________________________________________________ conv3d_3 (Conv3D) (None, 12, 12, 4, 256) 884992 _________________________________________________________________ max_pooling3d_3 (MaxPooling3 (None, 6, 6, 2, 256) 0 _________________________________________________________________ batch_normalization_3 (Batch (None, 6, 6, 2, 256) 1024 _________________________________________________________________ global_average_pooling3d (Gl (None, 256) 0 _________________________________________________________________ dense (Dense) (None, 512) 131584 _________________________________________________________________ dropout (Dropout) (None, 512) 0 _________________________________________________________________ dense_1 (Dense) (None, 1) 513 = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = Total params: 1352897 Trainable params: 1351873 Non - trainable params: 1024 _________________________________________________________________Copy the code

6. Training model

# Set dynamic learning rate
initial_learning_rate = 1e-4
lr_schedule = keras.optimizers.schedules.ExponentialDecay(
    initial_learning_rate, decay_steps=30, decay_rate=0.96, staircase=True
)
# compiler
model.compile(
    loss="binary_crossentropy",
    optimizer=keras.optimizers.Adam(learning_rate=lr_schedule),
    metrics=["acc"],)# Save the model
checkpoint_cb = keras.callbacks.ModelCheckpoint(
    "3d_image_classification.h5", save_best_only=True
)
# Define an early stop policy
early_stopping_cb = keras.callbacks.EarlyStopping(monitor="val_acc", patience=15)

epochs = 100
model.fit(
    train_dataset,
    validation_data=validation_dataset,
    epochs=epochs,
    shuffle=True,
    verbose=2,
    callbacks=[checkpoint_cb, early_stopping_cb],
)
Copy the code
Epoch 1/100 70/70-18S-Loss: 0.7115 - ACC: 0.5143-val_Loss: 0.8686 - Val_ACC: Epoch 2/100 70/70-15S-loss: 0.6595 - ACC: 0.6643 - val_loss: 0.8958 - val_ACC: 0.5000 Epoch 3/100 70/70-15S-loss: 0.6623-ACC: 0.6214-val_loss: 0.8238-val_ACC: 0.5000 Epoch 4/100...... 70/70-15S-loss: 0.5098-ACC: 0.7571-val_loss: 0.5803-VAL_ACC: 0.6667 Epoch 43/100 70/70-15S-loss: 0.5083-ACC: 0.7857-VAL_loss: 0.6009-Val_ACC: 0.6333 Epoch 44/100 70/70-15S-loss: 0.5276-ACC: 0.7429-VAL_loss: 0.5300-Val_ACC: 0.7000 Epoch 45/100 70/70-15S-loss: 0.5158-ACC: 0.7500-Val_loss: 0.5504-VAL_ACC: 0.6833 Epoch 46/100 70/70-15S-loss: 0.4648-ACC: 0.8357-Val_Loss: 0.5724-Val_ACC: 0.6833 < tensorflow. Python. Keras. Callbacks. History at 0 x1becc8326d8 >Copy the code

Note: Since the sample size is very small (only 200) and I did not specify random seeds. Therefore, you can expect that the results may differ significantly.

7. Visual model performance

fig, ax = plt.subplots(1.2, figsize=(20.3))
ax = ax.ravel()

for i, metric in enumerate(["acc"."loss"]):
    ax[i].plot(model.history.history[metric])
    ax[i].plot(model.history.history["val_" + metric])
    ax[i].set_title("Model {}".format(metric))
    ax[i].set_xlabel("epochs")
    ax[i].set_ylabel(metric)
    ax[i].legend(["train"."val"])
Copy the code

8. Prediction of single CT scan

# Load model
model.load_weights("3d_image_classification.h5")
prediction = model.predict(np.expand_dims(x_val[0], axis=0))0]
scores = [1 - prediction[0], prediction[0]]

class_names = ["normal"."abnormal"]
for score, name in zip(scores, class_names):
    print(
        "This model is %.2f percent confident that CT scan is %s"
        % ((100 * score), name)
    )
Copy the code
This model is 27.88 percent confident that CT scan is normal
This model is 72.12 percent confident that CT scan is abnormal
Copy the code

Note: This article refers to the official code, and some changes have been made.


🚀 from the column:100 Examples of Deep Learning

💖 praise after see, collect again, form a good habit! 💖