MosMedData: 新冠肺炎胸部 CT扫描数据集上基于3D-CNN实现二分类¶

!pip install nibabel
!pip install scipy

import os
import zipfile
import numpy as np
import paddle
import paddle.nn as nn
from paddle.io import Dataset
import random
import matplotlib.pyplot as plt
import sys 
from scipy import ndimage
import nibabel as nib

paddle.__version__

'2.3.2'

# 制作一个目录来存储数据。
os.makedirs("MosMedData")

# 在新创建的目录中解压缩数据。
with zipfile.ZipFile("data/data112472/ct-0.zip", "r") as z_fp:
    z_fp.extractall("./MosMedData/")

with zipfile.ZipFile("data/data112472/ct-23.zip", "r") as z_fp:
    z_fp.extractall("./MosMedData/")

def read_nifti_file(filepath):
    """读取和加载数据"""
    # 读取文件
    scan = nib.load(filepath)
    # 获取原始数据
    scan = scan.get_fdata()
    return scan

def normalize(volume):
    """数据归一化"""
    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):
    """跨 z 轴调整大小"""
    # 设置所需的深度
    desired_depth = 64
    desired_width = 128
    desired_height = 128
    # 获取当前深度
    current_depth = img.shape[-1]
    current_width = img.shape[0]
    current_height = img.shape[1]
    # 计算深度因子
    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
    # 旋转
    img = ndimage.rotate(img, 90, reshape=False)
    # 跨z轴调整大小
    img = ndimage.zoom(img, (width_factor, height_factor, depth_factor), order=1)
    return img

def process_scan(path):
    """读取和调整数据大小"""
    # 读取扫描文件
    volume = read_nifti_file(path)
    # 归一化
    volume = normalize(volume)
    # 调整宽度、高度和深度
    volume = resize_volume(volume)
    return volume

# 文件夹“CT-0”包括具有正常肺组织的CT扫描，
# 没有病毒性肺炎的CT。
normal_scan_paths = [
    os.path.join(os.getcwd(), "MosMedData/CT-0", x)
    for x in os.listdir("MosMedData/CT-0")
]

#文件夹“CT-23”包括具有感染的CT扫描，

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)))

CT scans with normal lung tissue: 100
CT scans with abnormal lung tissue: 100

abnormal_path = "abnormal_scans.npy"
normal_path = "normal_scans.npy"
# 因此处转换比较耗时，如果关闭此次运行，再次运行又需重新转换，所以我们可以缓存一下转换结果。
if os.path.exists(abnormal_path) and os.path.exists(normal_path):
    abnormal_scans = np.load(abnormal_path)
    normal_scans = np.load(normal_path)
else:
    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])
    np.save(abnormal_path,abnormal_scans)
    np.save(normal_path,normal_scans)

# 肺部扫描数据集
class dataset(Dataset):
    """
        步骤一：继承paddle.io.Dataset类
    """
    def __init__(self, abnormal_scans,normal_scans, mode='train'):
        """
        步骤二：实现构造函数，
        Args:
            abnormal_scans: 异常CT
            normal_scans: 正常CT
            mode:模式 train、val
        """
        super().__init__()
        self.imgs = []
        self.labels = []
        self.train = False
        abnormal_labels = np.array([1 for _ in range(len(abnormal_scans))])
        normal_labels = np.array([0 for _ in range(len(normal_scans))])
        if mode == 'train':
            x_train = np.concatenate((abnormal_scans[:70], normal_scans[:70]), axis=0)
            y_train = np.concatenate((abnormal_labels[:70], normal_labels[:70]), axis=0)
            self.imgs = x_train
            self.labels = y_train
            self.train = True
        else:
            x_val = np.concatenate((abnormal_scans[70:], normal_scans[70:]), axis=0)
            y_val = np.concatenate((abnormal_labels[70:], normal_labels[70:]), axis=0)
            self.imgs = x_val
            self.labels = y_val
            
    def __getitem__(self, index):
        """
        步骤三：实现__getitem__方法，定义指定index时如何获取数据，并返回单条数据（训练数据/测试数据，对应的标签）
        """
        img = self.imgs[index]
        lab = self.labels[index]
        if self.train:
            img = self.rotate(img)
        img = np.array(img).astype('float32')
        lab = np.array([lab], dtype="int64")
        img = np.expand_dims(img,3)
        return img, lab

    def __len__(self):
        """
        步骤四：实现__len__方法，返回数据集总数目
        """
        return len(self.labels)

    def rotate(self,volume):
        # 定义增强旋转角度
        angles = [-20, -10, -5, 5, 10, 20]
        # 随机选择一个角度
        angle = random.choice(angles)
        # 旋转CT
        volume = ndimage.rotate(volume, angle, reshape=False)
        volume[volume < 0] = 0
        volume[volume > 1] = 1
        return volume

train_dataset = dataset(abnormal_scans,normal_scans)
validation_dataset= dataset(abnormal_scans,normal_scans,mode='val')

train_loader = paddle.io.DataLoader(train_dataset, batch_size=10, shuffle=True)
validation_loader = paddle.io.DataLoader(validation_dataset, batch_size = 10, shuffle=False)

images, labels  = train_dataset.__getitem__(0)
print(images.shape)
image =np.squeeze(images, 3)
print("Dimension of the CT scan is:", image.shape)
plt.imshow(np.squeeze(image[:, :, 30]), cmap="gray")
plt.show()

(128, 128, 64, 1)
Dimension of the CT scan is: (128, 128, 64)

def plot_slices(num_rows, num_columns, width, height, data):
    """绘制CT的横断面"""
    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 = 6.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()


# 绘制CT横断面.
plot_slices(8, 8, 128, 128, image[:, :, :64])

def Accuracy(pre,label):
    a = np.squeeze(pre.cpu().numpy(), -1)
    a = a<0.5
    b = np.squeeze(label.cpu().numpy(), -1)==0
    c = np.sum(a==b)
    return c / len(a)

his_trainloss = []
his_trainacc = []
his_valloss = []
his_valacc = []

def trainer():
    scheduler = paddle.optimizer.lr.ExponentialDecay(learning_rate=0.0001, gamma=0.96, verbose=True)
    optim = paddle.optimizer.AdamW(learning_rate=scheduler,parameters=cnn.parameters())
    criterion = paddle.nn.BCELoss()
    num_epochs = 25
    max_accuracy = 0
    save_path ="best3dcnn"
    
    for epoch in range(num_epochs):
        
        traintotal,trainacc,trainloss = 0, 0, 0
        cnn.train()
        for i, (img, labels) in enumerate(train_loader):     
            predict = cnn(img)
            loss = criterion(predict, labels.astype(np.float32))
            optim.clear_grad()
            loss.backward()
            optim.minimize(loss)
            acc = Accuracy(predict, labels)
            trainacc += acc
            trainloss += loss.numpy()
            if i%10==0:
                print(f"epoch {epoch} iter {i} loss: {loss.numpy()}, accuracy {acc}")
            traintotal += 1
        valtotal, valacc,valloss = 0, 0, 0
        cnn.eval()
        for i ,(img, labels) in enumerate(validation_loader):
            with paddle.no_grad():
                predict = cnn(img)
                loss = criterion(predict, labels.astype(np.float32))
            valacc +=Accuracy(predict, labels)
            valloss += loss.numpy()
            valtotal+=1
        print(f"epoch {epoch} val Loss {valloss/valtotal} Accuray {valacc/valtotal} ")   
        # 只保存最优模型
        if valacc/valtotal>max_accuracy:
            max_accuracy = valacc/valtotal
            model_state = cnn.state_dict()
            paddle.save(model_state, save_path)
            print("max accuracy {}".format(max_accuracy))
        scheduler.step()
        his_trainloss.append(trainloss/traintotal)
        his_trainacc.append(trainacc/traintotal)
        his_valloss.append(valloss/valtotal)
        his_valacc.append(valacc/valtotal)

trainer()

Epoch 0: ExponentialDecay set learning rate to 0.0001.
epoch 0 iter 0 loss: [0.6301036], accuracy 0.5
epoch 0 iter 10 loss: [0.7426481], accuracy 0.4
epoch 0 val Loss [1.1324016] Accuray 0.5 
max accuracy 0.5
Epoch 1: ExponentialDecay set learning rate to 9.6e-05.
epoch 1 iter 0 loss: [0.7325362], accuracy 0.5
epoch 1 iter 10 loss: [0.5694233], accuracy 0.6
epoch 1 val Loss [1.0659176] Accuray 0.5 
Epoch 2: ExponentialDecay set learning rate to 9.216e-05.
epoch 2 iter 0 loss: [0.5374714], accuracy 0.6
epoch 2 iter 10 loss: [0.40864113], accuracy 0.8
epoch 2 val Loss [0.7484214] Accuray 0.5166666666666667 
max accuracy 0.5166666666666667
Epoch 3: ExponentialDecay set learning rate to 8.847359999999999e-05.
epoch 3 iter 0 loss: [0.4230544], accuracy 0.8
epoch 3 iter 10 loss: [0.4823562], accuracy 0.9
epoch 3 val Loss [0.6880954] Accuray 0.6166666666666667 
max accuracy 0.6166666666666667
Epoch 4: ExponentialDecay set learning rate to 8.493465599999999e-05.
epoch 4 iter 0 loss: [0.6360772], accuracy 0.6
epoch 4 iter 10 loss: [0.6532255], accuracy 0.6
epoch 4 val Loss [0.5975937] Accuray 0.7000000000000001 
max accuracy 0.7000000000000001
Epoch 5: ExponentialDecay set learning rate to 8.153726975999998e-05.
epoch 5 iter 0 loss: [0.676396], accuracy 0.6
epoch 5 iter 10 loss: [0.5197036], accuracy 0.7
epoch 5 val Loss [0.5965887] Accuray 0.7000000000000001 
Epoch 6: ExponentialDecay set learning rate to 7.827577896959998e-05.
epoch 6 iter 0 loss: [0.5222303], accuracy 0.7
epoch 6 iter 10 loss: [0.28494963], accuracy 0.9
epoch 6 val Loss [0.62934595] Accuray 0.6833333333333332 
Epoch 7: ExponentialDecay set learning rate to 7.514474781081598e-05.
epoch 7 iter 0 loss: [0.5931848], accuracy 0.7
epoch 7 iter 10 loss: [0.34113288], accuracy 0.9
epoch 7 val Loss [0.5483835] Accuray 0.7000000000000001 
Epoch 8: ExponentialDecay set learning rate to 7.213895789838334e-05.
epoch 8 iter 0 loss: [0.4956868], accuracy 0.8
epoch 8 iter 10 loss: [0.8255708], accuracy 0.5
epoch 8 val Loss [0.67082316] Accuray 0.5833333333333334 
Epoch 9: ExponentialDecay set learning rate to 6.9253399582448e-05.
epoch 9 iter 0 loss: [0.5112359], accuracy 0.8
epoch 9 iter 10 loss: [0.5681261], accuracy 0.6
epoch 9 val Loss [0.6324358] Accuray 0.65 
Epoch 10: ExponentialDecay set learning rate to 6.648326359915008e-05.
epoch 10 iter 0 loss: [0.4329744], accuracy 0.9
epoch 10 iter 10 loss: [0.7594853], accuracy 0.5
epoch 10 val Loss [0.5622995] Accuray 0.7000000000000001 
Epoch 11: ExponentialDecay set learning rate to 6.382393305518408e-05.
epoch 11 iter 0 loss: [0.3882313], accuracy 0.8
epoch 11 iter 10 loss: [0.48908353], accuracy 0.8
epoch 11 val Loss [0.54760677] Accuray 0.7666666666666666 
max accuracy 0.7666666666666666
Epoch 12: ExponentialDecay set learning rate to 6.12709757329767e-05.
epoch 12 iter 0 loss: [0.348113], accuracy 0.8
epoch 12 iter 10 loss: [0.41176265], accuracy 0.9
epoch 12 val Loss [0.5595346] Accuray 0.7833333333333332 
max accuracy 0.7833333333333332
Epoch 13: ExponentialDecay set learning rate to 5.882013670365765e-05.
epoch 13 iter 0 loss: [0.3314957], accuracy 1.0
epoch 13 iter 10 loss: [0.32749084], accuracy 1.0
epoch 13 val Loss [0.5549964] Accuray 0.7333333333333334 
Epoch 14: ExponentialDecay set learning rate to 5.6467331235511337e-05.
epoch 14 iter 0 loss: [0.36251408], accuracy 0.9
epoch 14 iter 10 loss: [0.23271704], accuracy 1.0
epoch 14 val Loss [0.5716041] Accuray 0.7333333333333334 
Epoch 15: ExponentialDecay set learning rate to 5.4208637986090884e-05.
epoch 15 iter 0 loss: [0.35667133], accuracy 0.9
epoch 15 iter 10 loss: [0.38534072], accuracy 0.9
epoch 15 val Loss [0.55513054] Accuray 0.7333333333333334 
Epoch 16: ExponentialDecay set learning rate to 5.2040292466647244e-05.
epoch 16 iter 0 loss: [0.37971848], accuracy 0.9
epoch 16 iter 10 loss: [0.37076467], accuracy 0.8
epoch 16 val Loss [0.546179] Accuray 0.6833333333333332 
Epoch 17: ExponentialDecay set learning rate to 4.9958680767981346e-05.
epoch 17 iter 0 loss: [0.4073751], accuracy 0.8
epoch 17 iter 10 loss: [0.35747206], accuracy 0.8
epoch 17 val Loss [0.50928646] Accuray 0.7666666666666667 
Epoch 18: ExponentialDecay set learning rate to 4.7960333537262095e-05.
epoch 18 iter 0 loss: [0.26402354], accuracy 1.0
epoch 18 iter 10 loss: [0.3062905], accuracy 0.9
epoch 18 val Loss [0.57567877] Accuray 0.7166666666666667 
Epoch 19: ExponentialDecay set learning rate to 4.6041920195771606e-05.
epoch 19 iter 0 loss: [0.19812493], accuracy 1.0
epoch 19 iter 10 loss: [0.37852132], accuracy 0.9
epoch 19 val Loss [0.55201167] Accuray 0.6666666666666666 
Epoch 20: ExponentialDecay set learning rate to 4.4200243387940746e-05.
epoch 20 iter 0 loss: [0.2532389], accuracy 1.0
epoch 20 iter 10 loss: [0.19800806], accuracy 1.0
epoch 20 val Loss [0.6551569] Accuray 0.7000000000000001 
Epoch 21: ExponentialDecay set learning rate to 4.243223365242311e-05.
epoch 21 iter 0 loss: [0.26303864], accuracy 1.0
epoch 21 iter 10 loss: [0.36150858], accuracy 0.8
epoch 21 val Loss [0.5345089] Accuray 0.6833333333333332 
Epoch 22: ExponentialDecay set learning rate to 4.073494430632618e-05.
epoch 22 iter 0 loss: [0.18942678], accuracy 1.0
epoch 22 iter 10 loss: [0.51950914], accuracy 0.7
epoch 22 val Loss [0.56237996] Accuray 0.6333333333333333 
Epoch 23: ExponentialDecay set learning rate to 3.910554653407313e-05.
epoch 23 iter 0 loss: [0.19643877], accuracy 0.9
epoch 23 iter 10 loss: [0.2246792], accuracy 1.0
epoch 23 val Loss [0.61442226] Accuray 0.7000000000000001 
Epoch 24: ExponentialDecay set learning rate to 3.754132467271021e-05.
epoch 24 iter 0 loss: [0.30469954], accuracy 1.0
epoch 24 iter 10 loss: [0.08427861], accuracy 1.0
epoch 24 val Loss [0.51900846] Accuray 0.7333333333333333 
Epoch 25: ExponentialDecay set learning rate to 3.60396716858018e-05.

fig, ax = plt.subplots(2, 1, figsize=(8, 12))
ax = ax.ravel()


ax[0].plot(his_trainloss)
ax[0].plot(his_valloss)
ax[0].set_title("Model {}".format("Loss"))
ax[0].set_xlabel("epochs")
ax[0].legend(["train", "val"])

ax[1].plot(his_trainacc)
ax[1].plot(his_valacc)
ax[1].set_title("Model {}".format("Accuracy"))
ax[1].set_xlabel("epochs")
# ax[1].set_ylabel(metric)
ax[1].legend(["train", "val"])

<matplotlib.legend.Legend at 0x7efe65dae250>

param_state_dict = paddle.load( "best3dcnn")
cnn.set_dict(param_state_dict)
cnn.eval() #预测模式
def predictor(index):
    imagen,label = validation_dataset.__getitem__(index)
    image = paddle.to_tensor(np.expand_dims(imagen, axis=0), dtype=paddle.float32)
    prediction = cnn(image)
    scores = [1 - prediction[0][0], prediction[0][0]]
    class_names = ["normal", "abnormal"]
    for score, name in zip(scores, class_names):
        print(
            "模型预测： %.2f 概率，这张CT %s"
            % ((100 * score), name)
        )
    print("真实标签为：",class_names[int(label[0])])
    data = np.rot90(imagen)
    data = np.transpose(data)
    data = np.reshape(data, (8, 8, 128, 128))
    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 = 6.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.text(-930,-850 , "The true label is:"+class_names[int(label[0])], size = 15,\
            color = "yellow", style = "italic", weight = "light",\
            bbox = dict(facecolor = "black", alpha = 0.3))  
    plt.text(-930,-790 , "This model is %.2f percent confident that CT scan is %s"% ((100 * scores[0]),class_names[0]), size = 15,\
            color = "w", style = "italic", weight = "light",\
            bbox = dict(facecolor = "black", alpha = 0.3))  
    plt.text(-930,-730 , "This model is %.2f percent confident that CT scan is %s"% ((100 * scores[1]),class_names[1]), size = 15,\
            color = "r", style = "italic", weight = "light",\
            bbox = dict(facecolor = "black", alpha = 0.3)) 
    plt.show()

predictor(2)

模型预测： 43.85 概率，这张CT normal
模型预测： 56.15 概率，这张CT abnormal
真实标签为： abnormal

predictor(11)

模型预测： 32.01 概率，这张CT normal
模型预测： 67.99 概率，这张CT abnormal
真实标签为： abnormal

predictor(55)

模型预测： 79.78 概率，这张CT normal
模型预测： 20.22 概率，这张CT abnormal
真实标签为： normal