基于 Luna16 数据集的肺部结节处理与模型训练实战

Jack电子实验室

已于 2024-12-15 21:01:33 修改

阅读量3k

点赞数 43

文章标签：网络深度学习 python pytorch 人工智能 conda 医学图像分割

于 2024-12-15 21:01:04 首次发布

本文链接：https://blog.csdn.net/qq_63860529/article/details/144492586

版权

基于 Luna16 数据集的肺部结节处理与模型训练实战

在医学影像分析领域，对肺部结节的准确检测和分割具有重要的临床意义，能够为早期肺部疾病的诊断提供有力支持。Luna 数据集作为肺部影像研究中的关键资源，为我们开展相关工作提供了丰富的数据基础。本文将详细介绍我对 Luna 数据集的处理流程，以及基于该数据集进行肺部结节分割所使用的模型训练过程，希望能为同行们提供一些参考和启发。

一、Luna 数据集简介

Luna 数据集主要来源于肺部的 CT 扫描影像，包含了大量的肺部病例数据，这些数据对于肺部结节的研究至关重要。数据集中的每一个病例都包含了详细的肺部 CT 图像信息，并且对于存在的肺部结节，还提供了精确的标注信息，如结节的位置、大小、形状等，这些标注信息为我们进行模型训练提供了准确的目标参考，使得我们能够开发出针对肺部结节检测和分割的有效模型，从而辅助医生更准确地诊断肺部疾病。

二、数据集处理过程

（一）数据预处理

数据读取与格式转换
首先，使用 SimpleITK 库来读取 Luna 数据集中的 .mhd 格式的图像文件，并将其转换为 numpy 数组格式，以便后续进行处理。例如：

import SimpleITK as sitk
import numpy as np

# 读取图像文件
itk_img = sitk.ReadImage("example.mhd")
img_array = sitk.GetArrayFromImage(itk_img)

同时，对于一些可能需要的图像格式转换操作，如将原始的医学影像格式转换为常见的 .png 或 .jpg 格式，以方便在某些特定场景下的查看和处理，也可以在这一步骤中进行相应的代码实现（若实际情况有此需求）。
2. 像素值归一化
为了使图像数据在模型训练中具有更好的稳定性和一致性，对图像的像素值进行归一化处理是必不可少的。通过定义 normalizePlanes 函数，将 CT 图像的 HU 值（Hounsfield Units）映射到 [0, 255] 的范围：

def normalizePlanes(npzarray):
    maxHU = 400
    minHU = -1000
    npzarray = (npzarray - minHU) / (maxHU - minHU)
    npzarray[npzarray > 1] = 1
    npzarray[npzarray < 0] = 0
    npzarray *= 255
    return (npzarray.astype(int))

这样的归一化操作能够使不同病例的图像数据在数值上具有可比性，有利于模型的学习和训练。

（二）制作肺部结节掩码

根据数据集中提供的肺部结节标注信息，通过 make_mask 函数来生成对应的肺部结节掩码。掩码是一个二维数组，用于精确标记出结节在图像中的位置，其实现过程如下：

def make_mask(center, diam, z, width, height, spacing, origin):
    mask = np.zeros([height, width])  # 匹配图像

    # 定义结节所在的体素范围
    v_center = (center - origin) / spacing
    v_diam = int(diam / spacing[0] + 5)
    v_xmin = np.max([0, int(v_center[0] - v_diam) - 5])
    v_xmax = np.min([width - 1, int(v_center[0] + v_diam) + 5])
    v_ymin = np.max([0, int(v_center[1] - v_diam) - 5])
    v_ymax = np.min([height - 1, int(v_center[1] + v_diam) + 5])

    v_xrange = range(v_xmin, v_xmax + 1)
    v_yrange = range(v_ymin, v_ymax + 1)

    x_data = [x * spacing[0] + origin[0] for x in range(width)]
    y_data = [x * spacing[1] + origin[1] for x in range(height)]

    # 结节周围全都填充 1
    for v_x in v_xrange:
        for v_y in v_yrange:
            p_x = spacing[0] * v_x + origin[0]
            p_y = spacing[1] * v_y + origin[1]
            if np.linalg.norm(center - np.array([p_x, p_y, z])) <= diam:
                mask[int((p_y - origin[1]) / spacing[1]), int((p_x - origin[0]) / spacing[0])] = 1.0

    return (mask)

在这个函数中，首先根据结节的中心坐标、直径以及图像的空间信息（如原点、间距等）计算出结节在体素坐标系下的大致范围，然后通过遍历这个范围内的体素，根据与结节中心的距离判断是否在结节内，从而在掩码中相应位置进行标记。

（三）数据划分与保存

经过上述预处理和掩码制作步骤后，将数据集按照一定的比例划分为训练集、测试集，并将处理后的图像数据和掩码数据保存为 .npy 格式，以便后续快速加载和使用。例如：

# 假设已经完成了数据预处理和掩码制作，得到了 imgs 和 masks 数组
np.save("trainImages.npy", imgs_train)
np.save("trainMasks.npy", masks_train)
np.save("testImages.npy", imgs_test)
np.save("testMasks.npy", masks_test)

这样，我们就完成了对 Luna 数据集的初步处理，得到了适合模型训练和测试的数据集格式。

三、模型构建与训练

（一）UNet 模型架构

选择 UNet 作为我们的分割模型，它是一种经典的卷积神经网络架构，在图像分割任务中表现出色。以下是 UNet 模型的具体实现代码：

from keras.models import *
from keras.layers import *
from keras.optimizers import *

def unet(pretrained_weights=None, input_size=(512, 512, 1)):
    inputs = Input(input_size)
    conv1 = Conv2D(64, 3, activation='relu', padding='same', kernel_initializer='he_normal')(inputs)
    conv1 = Conv2D(64, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv1)
    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)
    conv2 = Conv2D(128, 3, activation='relu', padding='same', kernel_initializer='he_normal')(pool1)
    conv2 = Conv2D(128, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv2)
    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)
    conv3 = Conv2D(256, 3, activation='relu', padding='same', kernel_initializer='he_normal')(pool2)
    conv3 = Conv2D(256, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv3)
    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)
    conv4 = Conv2D(512, 3, activation='relu', padding='same', kernel_initializer='he_normal')(pool3)
    conv4 = Conv2D(512, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv4)
    drop4 = Dropout(0.5)(conv4)
    pool4 = MaxPooling2D(pool_size=(2, 2))(drop4)

    conv5 = Conv2D(1024, 3, activation='relu', padding='same', kernel_initializer='he_normal')(pool4)
    conv5 = Conv2D(1024, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv5)
    drop5 = Dropout(0.5)(conv5)

    up6 = Conv2D(512, 2, activation='relu', padding='same', kernel_initializer='he_normal')(UpSampling2D(size=(2, 2))(drop5))
    merge6 = concatenate([drop4, up6], axis=3)
    conv6 = Conv2D(512, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge6)
    conv6 = Conv2D(512, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv6)

    up7 = Conv2D(256, 2, activation='relu', padding='same', kernel_initializer='he_normal')(UpSampling2D(size=(2, 2))(conv6))
    merge7 = concatenate([conv3, up7], axis=3)
    conv7 = Conv2D(256, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge7)
    conv7 = Conv2D(256, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv7)

    up8 = Conv2D(128, 2, activation='relu', padding='same', kernel_initializer='he_normal')(UpSampling2D(size=(2, 2))(conv7))
    merge8 = concatenate([conv2, up8], axis=3)
    conv8 = Conv2D(128, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge8)
    conv8 = Conv2D(128, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge8)

    up9 = Conv2D(64, 2, activation='relu', padding='same', kernel_initializer='he_normal')(UpSampling2D(size=(2, 2))(conv8))
    merge9 = concatenate([conv1, up9], axis=3)
    conv9 = Conv2D(64, 3, activation='relu', padding='same', kernel_initializer='he_normal')(merge9)
    conv9 = Conv2D(64, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv9)
    conv9 = Conv2D(2, 3, activation='relu', padding='same', kernel_initializer='he_normal')(conv9)
    conv10 = Conv2D(1, 1, activation='sigmoid')(conv9)

    model = Model(inputs=inputs, outputs=conv10)

    model.compile(optimizer=Adam(lr=1e-4), loss='binary_crossentropy', metrics=['accuracy'])
    
    if(pretrained_weights):
        model.load_weights(pretrained_weights)
        
    return model

该模型通过不断的下采样（卷积和池化操作）提取图像特征，然后通过上采样并结合相应层次的下采样特征进行特征融合，最终输出分割结果。这种结构能够有效地捕捉图像中的细节信息，对于肺部结节这种小目标的分割具有较好的效果。

（二）模型训练与评估指标

训练过程
在模型训练过程中，首先加载之前处理好并保存的训练集和测试集数据，包括图像数据和对应的掩码数据，并对训练数据进行标准化处理，以提高模型的训练效果。然后创建 UNet 模型实例，并设置模型检查点（ModelCheckpoint），用于保存训练过程中损失最小的模型权重到 unet.hdf5 文件。使用训练数据对模型进行训练，训练的批次大小为 16，共训练 100 个轮次，并且在训练过程中打乱数据顺序，以防止模型过拟合。

def train_and_predict(use_existing):
    print('-'*30)
    print('Loading and preprocessing train data...')
    print('-'*30)
    imgs_train = np.load(working_path+"trainImages.npy").astype(np.float32)
    imgs_train = np.transpose(imgs_train, (0, 2, 3, 1))
    imgs_mask_train = np.load(working_path+"trainMasks.npy").astype(np.float32)
    imgs_mask_train = np.transpose(imgs_mask_train, (0, 2, 3, 1))

    imgs_test = np.load(working_path+"testImages.npy").astype(np.float32)
    imgs_test = np.transpose(imgs_test, (0, 2, 3, 1))
    imgs_mask_test_true = np.load(working_path+"testMasks.npy").astype(np.float32)
    imgs_mask_test_true = np.transpose(imgs_mask_test_true, (0, 2, 3, 1))

    mean = np.mean(imgs_train)  # mean for data centering
    std = np.std(imgs_train)  # std for data normalization

    imgs_train -= mean  # images should already be standardized, but just in case
    imgs_train /= std

    print('-'*30)
    print('Creating and compiling model...')
    print('-'*30)
    model = unet()
    model_checkpoint = ModelCheckpoint('unet.hdf5', monitor='loss', save_best_only=True)
    if use_existing:
        model.load_weights('./unet.hdf5')
    print('-'*30)
    print('Fitting model...')
    print('-'*30)
    model.fit(imgs_train, imgs_mask_train, batch_size=16, epochs=100, verbose=1, shuffle=True,
              callbacks=[model_checkpoint])

评估指标
为了准确评估模型的分割性能，定义了 Dice 系数作为主要的评估指标。Dice 系数常用于衡量图像分割结果与真实标签之间的相似度，取值范围在 0 到 1 之间，越接近 1 表示分割效果越好。以下是基于 Keras 后端和 Numpy 分别实现的 Dice 系数计算函数：

from keras import backend as K

def dice_coef(y_true, y_pred):
    smooth = 1.
    y_true_f = K.flatten(y_true)
    y_pred_f = K.flatten(y_pred)
    intersection = K.sum(y_true_f * y_pred_f)
    return (2. * intersection + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)

def dice_coef_np(y_true, y_pred):
    y_true_f = y_true.flatten()
    y_pred_f = y_pred.flatten()
    intersection = np.sum(y_true_f * y_pred_f)
    return (2. * intersection + smooth) / (np.sum(y_true_f) + np.sum(y_pred_f) + smooth)

同时，为了在模型训练中直接使用 Dice 系数相关的损失来优化模型，还定义了 Dice 系数损失函数：

def dice_coef_loss(y_true, y_pred):
    return -dice_coef(y_true, y_pred)

（三）模型预测与结果可视化

预测过程
在模型训练完成后，加载保存的最佳模型权重，对测试集数据进行预测，得到预测的掩码结果并保存为 masksTestPredicted.npy 文件。

    print('-'*30)
    print('Loading saved weights...')
    print('-'*30)
    model.load_weights('./unet.hdf5')

    print('-'*30)
    print('Predicting masks on test data...')
    print('-'*30)
    num_test = len(imgs_test)
    imgs_mask_test = np.ndarray([num_test,512,512,1],dtype=np.float32)
    for i in range(num_test):
        imgs_mask_test[i] = model.predict([imgs_test[i:i+1]], verbose=0)[0]
    np.save('masksTestPredicted.npy', imgs_mask_test)

结果可视化
通过定义 visualize_results 函数，对测试集上的原始图像、真实掩码以及预测掩码进行可视化展示，以便直观地观察模型的分割效果。通过循环展示前 num_images（默认为 5）个测试样本的结果，能够清晰地看到模型在不同病例上的表现情况。

import matplotlib.pyplot as plt

def visualize_results(imgs_test, imgs_mask_test_true, imgs_mask_test, num_images=5):
    plt.figure(figsize=(15, 5 * num_images))
    
    for i in range(num_images):
        plt.subplot(num_images, 3, i * 3 + 1)
        plt.imshow(imgs_test[i].squeeze(), cmap='gray')
        plt.title('Original Image')
        plt.axis('off')

        plt.subplot(num_images, 3, i * 3 + 2)
        plt.imshow(imgs_mask_test_true[i,..., 0], cmap='gray')  
        plt.title('True Mask')
        plt.axis('off')

        plt.sub