1. 背景

1.影像组学，常见的套路对影像图像的病灶进行勾画，制作mask文件。通过特征软件或者pyradiomics库对原始数据文件结合mask的范围，提取几百到几千不等的特征。然后通过统计学对特征进行筛选，然后通过机器学习对筛选后的进行建模，进行分类或者生存预测等任务。例如这篇文献

《 Development and Validation of a Radiomics Nomogram
for Differentiating Mycoplasma Pneumonia and
Bacterial Pneumonia》

2.也有作者使用深度学习的方法，使用ImageNet数据库上进行了预训练的2D分类模型，删除了网络顶部的最后一个完全连接的层，并使用全局最大池化来取特征图每一层的最大值，以将特征图转换为原始值。然后将这些提取的特征输入到机器学习进行建模。例如这篇文献

《Computed tomography-based deep-learning prediction of neoadjuvant
chemoradiotherapy treatment response in esophageal squamous cell
carcinoma

3.也有作者使用深度学习3D分类模型，使用大量医学数据库上进行了预训练的3D分类模型，对自己的数据进行微调，然后提取全连接的特征，再去建模。例如这篇。

《Deep learning and radiomics analysis for prediction of placenta
invasion based on T2WI》

4.影像数据一般比较少，使用深度学习的时候可以使用迁移学习，但是像ImageNet这种规模的数据集， 和医学数据差异比较大。有作者就提出使用大量医学数据训练过的模型，对新的数据进行微调，能否提高微调结果的性能。

《Med3D: Transfer Learning for 3D Medical Image Analysis 》

这篇文章作者使用8种分割医疗数据进行训练。得到训练的好的编码器，作为其它分类、检测、分割任务进行微调得到不错的效果。

5.本项目也使用MedicalNet训练好的权重，因为权重文件是pytorch格式，通过PaddlePaddle的X2Paddle ，在本地转换好PaddlePaddle格式的权重文件。对MedicalNet模型尾部的分割头去掉，改成全连接层。然后加载权重，固定权重。相当于使用预训练模型的权重提取3D图像的特征，然后通过全连接层对特征进行分类。模型结构的修改参考“brain_mri_medicalnet_classification” 的项目。在这里也感谢realty66 的知识分享。

6.整个项目的示意图如下所示：

2. 数据处理

数据来源：

BraTS20，HGG有100例，LGG有76例。每一例中都有五个文件，分别是T1、T2、T1CE、Flair、Seg，四个模态文件和一个分割文件。

数据处理

每个原始数据都带有一个对应的mask数据。在输入3D分类模型前，为了让模型更好对病灶区域进行特征提取。根据mask文件对病灶区域进行裁剪，首先根据mask的范围对3D数据三个维度上进行裁剪。然后对裁剪后的数据非前景的区域用0进行填充（这是对磁共振数据的处理，假如是CT数据应该设置窗宽窗位，然后缩放到0~255，再对背景填充0），这样得到的数据只有前景的信息。

#安装 SimpleITK，用来读取医疗数据
!pip install SimpleITK

!unzip -o /home/aistudio/data/data164244/MedicalNet.zip -d /home/aistudio/
!unzip -o /home/aistudio/data/data158320/brats176.zip -d /home/aistudio/work

#生成txt文档，每一行:xxx.nii 空格 类别
import os
import random
random.seed(2022)
classifier = {"LGG":0,"HGG":1}

path_origin = '/home/aistudio/work/brats176/'
train_txt = open(os.path.join(path_origin,'train_list.txt'),'w')
val_txt = open(os.path.join(path_origin,'val_list.txt'),'w')

for k in classifier.keys():
    classifier_path = os.path.join(path_origin,k)
    fileList = []
    for root, dirs, files in os.walk(classifier_path):
            for fileObj in files:
                fileList.append(os.path.join(root, fileObj))
    #只使用T2模态
    files = list(filter(lambda x: x.endswith('t2.nii'), fileList))
    random.shuffle(files)
    rate = int(len(files) * 0.8)#训练集和测试集8：2
    for i,f_path in enumerate(files):
        if i < rate:
            train_txt.write(f_path + ' ' + str(classifier[k])+ '\n')
        else:
            val_txt.write(f_path + ' ' + str(classifier[k])+ '\n')
train_txt.close()
val_txt.close()
print('完成')

完成

3. Dataset

DataSet主要构建中，已经包括了对原始数据的裁剪，并且Transform中使用了对前景进行归一化，和重采样到32x64x64

#创建DataSet 和Transform 用来读取医疗数据，并转换成numpy，归一化，重采样合适的大小
from paddle.io import Dataset 
import numpy as np
import os
import SimpleITK as sitk
from transforms import transforms as T

SIZE = (32,64,64)
train_transforms =T.Compose( [
    T.MRINormalize(),
    T.Resize3D(target_size=SIZE),#重采样
])

val_transforms = T.Compose([
    T.MRINormalize(),
    T.Resize3D(target_size=SIZE),
])
class MyDataset(Dataset):
    def __init__(self, data_dir, mode='train',transform=None):
        super().__init__()
        self.data_dir = data_dir
        self.transform = transform
        self.data_paths = []
        self.labels = []

        if mode == 'train':
            with open(os.path.join(self.data_dir, "train_list.txt"), "r", encoding="utf-8") as f:
                info = f.readlines()
            for img_info in info:
                img_path, label = img_info.split()
                self.data_paths.append(img_path)
                self.labels.append(int(label))

        elif mode == 'val':
            with open(os.path.join(self.data_dir, "val_list.txt"), "r", encoding="utf-8") as f:
                info = f.readlines()
            for img_info in info:
                img_path, label = img_info.split()
                self.data_paths.append(img_path)
                self.labels.append(int(label))

    def __getitem__(self, index):
        img_path = self.data_paths[index]
        mask_path = img_path.replace('t2.nii','seg.nii')
        data = self.read_Nifit(img_path)
        mask = self.read_Nifit(mask_path)
        data = self.maskcroppingbox(data,mask)
        if data.ndim != 3:
            raise("  ")
        if self.transform is not None:
            data = self.transform(data)
        data = np.expand_dims(data, axis=0)
        label = self.labels[index]
        label = np.array([label], dtype="int32")
        #data = np.ascontiguousarray(data)
        return data, label

    def read_Nifit(self,path):
        sitkImg = sitk.ReadImage(path)
        npImg = sitk.GetArrayFromImage(sitkImg)
        npImg = npImg.astype('float32')
        return npImg

    def maskcroppingbox(self,img,mask, use2D=False):
        #根据mask范围进行裁剪
        mask_2 = np.argwhere(mask)
        (zstart, ystart, xstart), (zstop, ystop, xstop) = mask_2.min(axis=0), mask_2.max(axis=0) + 1
        
        roi_image = img[zstart-1:zstop+1,ystart:ystop,xstart:xstop]
        roi_mask = mask[zstart-1:zstop+1,ystart:ystop,xstart:xstop]
        roi_image[roi_mask<1]=0
        return roi_image

    def __len__(self):
        return len(self.data_paths)

train_dataset = MyDataset('./work/brats176/',mode='train',transform=train_transforms)
val_dataset = MyDataset('./work/brats176/',mode='val',transform=train_transforms)

#展示处理后的数据中某一层
import matplotlib.pyplot as plt

img,label = train_dataset[1]
d,h,w = img[0].shape
print(img[0].shape)
plt.imshow(img[0][int(d/2),:,:],'gray')
plt.show()

(32, 64, 64)

4.建模

使用PaddlePaddle的高层API，创建模型和加载预训练权重，模型结构中已经对卷积权重进行固定，不再反向传播。只用最后的全连接层进行反向重播。

#创建Resnet3D分类模型并加载预训练模型，
from model.model import generate_model
import os
import paddle
layers = 50
load_layer_state_dict = paddle.load( "/home/aistudio/MedicalNet/resnet_"+str(layers)+".pdiparams")
resnet3d = generate_model(layers,checkpoint=load_layer_state_dict)
model = paddle.Model(resnet3d)

model.summary((4,1,32,64, 64))

#使用较小的学习率进行微调
opt = paddle.optimizer.Adam(learning_rate=0.0001, parameters=model.parameters())
model.prepare(opt,
              paddle.nn.CrossEntropyLoss(),
              paddle.metric.Accuracy())
vdl_callback = paddle.callbacks.VisualDL(log_dir='log')
#callbacks的使用可以学习“nanting03”大佬的“一文搞懂Paddle2.0中的Callbacks”
#一文搞懂Paddle2.0中的Callbacks：https://aistudio.baidu.com/aistudio/projectdetail/1318174
class SaveBestModel(paddle.callbacks.Callback):
    def __init__(self, target=0.5, path='./best_model', verbose=0):
        self.target = target
        self.epoch = None
        self.path = path

    def on_epoch_end(self, epoch, logs=None):
        self.epoch = epoch

    def on_eval_end(self, logs=None):
        if logs.get('acc') > self.target:
            self.target = logs.get('acc')
            self.model.save(self.path)
            print('best acc is {} at epoch {}'.format(self.target, self.epoch))

save_dir = 'save_models'+str(layers)
callback_savebestmodel = SaveBestModel(target=0.5, path=save_dir+'/best_model')
model.fit(
    train_data=train_dataset, 
    eval_data=val_dataset, 
    batch_size=4, 
    epochs=20, 
    eval_freq=1, 
    log_freq=20, 
    save_dir=save_dir, 
    save_freq=3, 
    verbose=1, 
    drop_last=False, 
    shuffle=True,
    num_workers=1, 
    callbacks=[vdl_callback,callback_savebestmodel]
    )

5. 预测数据并绘制图像

对验证集进行预测并绘制ROC曲线、DCA曲线、混淆矩阵图和校准曲线。并且对比了使用Resnet3D_18和Resnet3D_34和Resnet3D_50预训练权重的效果。

import pandas as pd
from model.model import generate_model
import paddle
Layers = [18, 34, 50]
for layer in Layers:
    resnet3d = generate_model(layer)
    model = paddle.Model(resnet3d)
    #加载训练好的模型
    model.load('/home/aistudio/save_models'+str(layer)+'/best_model.pdparams')
    classes = [0,1]#类别
    ids = list() #存储行名
    labels = list() # 存储真实标签
    pred_labels = list() #存储预测标签

    pred_scores = list() #存储预测概率
    for id,data in enumerate(val_dataset) :
        img,label = data
        img= np.expand_dims(img, axis=0)
        result = model.predict_batch(img)
        #经过Softmax 转换成【0,1】范围的概率，每个类别的概率相加等于1
        pred_score = paddle.squeeze(paddle.nn.Softmax()(paddle.to_tensor(result))).numpy()
        index = np.argmax(result)
        predict_label = classes[index]

        ids.append(id)
        labels.append(label[0])
        pred_labels.append(predict_label)
        pred_scores.append(list(pred_score))

    data ={
        'ids':ids,
        'labels':labels,
        'pred_labels':pred_labels,
        'pred_scores':pred_scores
    }
    #存储成DataFrame方便后续操作
    if layer == 18:
        df18 = pd.DataFrame(data)
    elif layer == 34:
        df34 = pd.DataFrame(data)
    elif layer == 50:
        df50 = pd.DataFrame(data)
#展示某个DataFrame数据
df34.head()

W0821 20:02:50.111349   171 gpu_resources.cc:61] Please NOTE: device: 0, GPU Compute Capability: 7.0, Driver API Version: 11.2, Runtime API Version: 10.1
W0821 20:02:50.115998   171 gpu_resources.cc:91] device: 0, cuDNN Version: 7.6.

	ids	labels	pred_labels	pred_scores
0	0	0	0	[0.7040171, 0.29598293]
1	1	0	0	[0.7981601, 0.20183994]
2	2	0	0	[0.5934987, 0.4065013]
3	3	0	0	[0.6932341, 0.30676588]
4	4	0	0	[0.8255727, 0.17442727]

5.1绘制ROC曲线(区分度)

评价预测模型的区分能力的指标可以使用ROC曲线在面积（AUC），AUC越大，说明模型的判别区分能力越好。

可以看到Resnet3D-34的AUC为0.88，三个模型中最高，Resnet3D-50AUC最低，但是三个模型都在0.8以上。

import sklearn
from sklearn.metrics import roc_curve,auc,accuracy_score,confusion_matrix
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings("ignore")
def fun(x,index):
    #返回列表中第几个元素
    return x[index]
#绘制roc曲线
plt.figure(figsize=(8,8))
colors=['red','yellow','blue']
for index, df in enumerate([df18,df34,df50]) :
    y_true = df['labels']
    y_pred = df['pred_labels']
    y_score = df['pred_scores'].apply(lambda x:fun(x,1))
    #计算roc和 auc
    fpr, tpr, thresholds=roc_curve(y_true,y_score,pos_label=1)
    ACC = accuracy_score(y_true,y_pred)
    AUC = auc(fpr, tpr) 
    plt.plot(fpr,tpr,colors[index],label='Resnet3D-'+str(Layers[index])+' ACC=%0.2f AUC = %0.2f'% (ACC ,AUC))
plt.legend(loc='lower right',fontsize = 12)
plt.plot([0,1],[0,1],color='black',linestyle='dashed')
plt.ylabel('True Positive Rate',fontsize = 14)
plt.xlabel('Flase Positive Rate',fontsize = 14)
plt.show()

5.2 绘制混淆矩阵

从混淆矩阵中，可以直观看到那些类别是模型分类对，那些分类错误的。
从结果中看到Resnet3D-34对LGG类别的预测效果最好，Resnet3D-50对HGG类别预测效果最好.
综合来说Resnet3D-34对LGG和HGG类别预测效果比其它三个模型都要好。

plt.figure(figsize=(15,5))
for index, df in enumerate([df18,df34,df50]) :
    y_true = df['labels']
    y_pred = df['pred_labels']
    plt.subplot(1,3,index+1)
    sns.heatmap(confusion_matrix(y_true, y_pred),annot=True, cmap='Blues')
    plt.xlabel('Predicted Label')
    plt.ylabel('True Label')
    plt.title("Resnet3D-"+str(Layers[index]))
plt.show()

5.3 绘制DCA曲线

判断模型的区分能力可以使用ROC曲线，但是ROC曲线体现模型的准确率、灵敏度等指标，并没有关注不同模型中不同切点所带来的获益与风险的关系。而决策曲线分析（decision curve analysis, DCA）是衡量临床实用性的一种广泛使用的方法。
DCA曲线中，三个模型都比较接近，但是可以看出黄色曲线(Resnet3D-34)最靠近右上角。

#之前都是用R绘制DCA曲线，现在在网上找到用python绘制DCA方法，感谢作者的分享
#用python绘制DCA曲线参考资料http://t.csdn.cn/dQv8s

def calculate_net_benefit_model(thresh_group, y_pred_score, y_label):
    net_benefit_model = np.array([])
    for thresh in thresh_group:
        y_pred_label = y_pred_score > thresh
        tn, fp, fn, tp = confusion_matrix(y_label, y_pred_label).ravel()
        n = len(y_label)
        net_benefit = (tp / n) - (fp / n) * (thresh / (1 - thresh))
        net_benefit_model = np.append(net_benefit_model, net_benefit)
    return net_benefit_model

def calculate_net_benefit_all(thresh_group, y_label):
    net_benefit_all = np.array([])
    tn, fp, fn, tp = confusion_matrix(y_label, y_label).ravel()
    total = tp + tn
    for thresh in thresh_group:
        net_benefit = (tp / total) - (tn / total) * (thresh / (1 - thresh))
        net_benefit_all = np.append(net_benefit_all, net_benefit)
    return net_benefit_all

def plot_DCA(ax, thresh_group, net_benefit_model, net_benefit_all,color='crimson',index=18):

    ax.plot(thresh_group, net_benefit_model, color = color, label = 'Resnet3D-'+str(index))

    #Figure Configuration， 美化一下细节
    ax.set_xlim(0,1)
    ax.set_ylim(net_benefit_model.min() - 0.15, net_benefit_model.max() + 0.15)#adjustify the y axis limitation
    ax.set_xlabel(
        xlabel = 'High Risk Threshold', 
        fontdict= {'family': 'Times New Roman', 'fontsize': 13}
        )
    ax.set_ylabel(
        ylabel = 'Net Benefit', 
        fontdict= {'family': 'Times New Roman', 'fontsize': 13}
        )
    ax.grid('off')
    ax.spines['right'].set_color((0.8, 0.8, 0.8))
    ax.spines['top'].set_color((0.8, 0.8, 0.8))
    ax.legend(loc = 'upper right')
    return ax

fig, ax = plt.subplots()
fig.set_size_inches(10, 8)
for index, df in enumerate([df18,df34,df50]) :
    y_label = df['labels']
    y_pred_score = df['pred_scores'].apply(lambda x:fun(x,1))

    thresh_group = np.arange(0,1,0.01)
    net_benefit_model = calculate_net_benefit_model(thresh_group, y_pred_score, y_label)
    net_benefit_all = calculate_net_benefit_all(thresh_group, y_label)
    ax = plot_DCA(ax, thresh_group, net_benefit_model, net_benefit_all,color=colors[index],index= Layers[index])
ax.plot(thresh_group, net_benefit_all, color = 'black',label = 'ALL')
ax.plot((0, 1), (0, 0), color = 'black', linestyle = ':', label = 'None')
plt.show()

5.4绘制校准曲线（校准度）

一个模型区分度好（AUC高），但是未必能提现预测概率和观测概率的一致性问题。校准曲线就能体现，横坐标是预测的概率，纵坐标是观测的概率，如果校准曲线靠近中间的虚线，代表模型的校准度好。
可以到Resnet3D_34模型的校准曲线最接近黑色虚线。

from sklearn.calibration import calibration_curve
for index, df in enumerate([df18,df34,df50]) :
    y_true = df['labels']
    y_score = df['pred_scores'].apply(lambda x:fun(x,1))
    prob_true, prob_pred = calibration_curve(list(y_true), list(y_score), n_bins=5)
    plt.plot(prob_pred, prob_true, linewidth=1, color=colors[index], marker="o",label="Mean value")
plt.plot([0,1],[0,1],color='black',linestyle='dashed')
plt.legend(["Resnet3D-18","Resnet3D-34","Resnet3D-50"],loc="upper left")
plt.xlabel('Predicted Probability of hypertension')
plt.ylabel('Actual probability')

Text(0,0.5,'Actual probability')

6.总结

现在搞影像组学的方法多种多样，有纯机器学习的，也有深度学习加机器学习的。这个项目使用Resnet3D网络，使用经过医学图像训练好的权重，然后固定卷积权重，然后模型结尾增加全连接层。这样相当于用3D卷积提取特征，然后全连接层进行训练、分类。项目比较了三个模型，最终发现Resnet3D-34的模型对胶质瘤的分类效果比较好

此文章为搬运
原项目链接