YOLOv8检测LUNA16肺结节实战（一）：数据预处理

说明：

肺结节检测已经相对成熟，这里将尝试用最新的YOLOv8检测LUNA16数据集中的肺结节。这里的代码参考了CSDN的多篇博客，然后通过调试，把整个预处理过程的代码整合到一个文件（preprocess.py）。代码已经成功运行，如有错误，欢迎指正。

注：1）该代码生成的是2维的图片（非3D）；2）未进行flip（不知道为什么需要该操作，欢迎告知）。

LUNA16数据集：https://luna16.grand-challenge.org/Download/

主要参考：
图像预处理+生成训练图片（.jpg）

提供了详细的代码（参考了.xml标签文件生成函数）

提供了LUNA16数据集的基础知识

1、安装YOLOv8：

# 创建虚拟环境：
conda create -n yolov8 python=3.8
conda activate yolov8

# 安装PyTorch:
pip install torch==1.9.0+cu111 torchvision==0.10.0+cu111  -f https://download.pytorch.org/whl/torch_stable.html

# 安装YOLOv8:
pip install ultralytics -i https://pypi.tuna.tsinghua.edu.cn/simple

# 安装其它依赖项（未列举完全）：
pip install SimpleITK scikit-learn scikit-image -i https://pypi.tuna.tsinghua.edu.cn/simple/

2、设置基础路径：

luna_base_path = "/media/bsuo/Elements/CT_image/LUNA16/subset_data/"
annotation_path = "/media/bsuo/Elements/CT_image/LUNA16/CSVFILES/annotations.csv"
output_npy_base_path = "/media/bsuo/Seagate/CT_image/LUNA_YOLOv8/npy/"
output_data_base_path = "/media/bsuo/Seagate/CT_image/LUNA_YOLOv8/data/"

3、CT扫描像素值归一化处理函数：

def normalizePlanes(npzarray):
    maxHU = 400
    minHU = -1000
    npzarray = (npzarray - minHU) / (maxHU - minHU)  # 像素直归一化处理
    npzarray[npzarray > 1] = 1
    npzarray[npzarray < 0] = 0
    npzarray *= 255  # 所有归一的像素乘以255，即1-255范围

    return npzarray.astype(int)

4、生成肺结节在CT扫描中掩码的函数：

def make_mask(center, diam, z, width, height, spacing, origin):
    mask = np.zeros([height, width])  # 匹配图像，生成一个确定高和宽的，元素全为0的掩码

    # 定义结节所在的体素范围
    v_center = (center - origin) / spacing  # 将肺结节中心坐标从世界坐标系转换为体素坐标系
    v_diam = int(diam / spacing[0] + 5)  # 计算肺结节直径在体素坐标系中的大小
    v_xmin = np.max([0, int(v_center[0] - v_diam) - 5])  # 计算肺结节所在的 x 方向上的最小体素坐标
    v_xmax = np.min([width - 1, int(v_center[0] + v_diam) + 5])  # 计算肺结节所在的 x 方向上的最大体素坐标
    v_ymin = np.max([0, int(v_center[1] - v_diam) - 5])  # 计算肺结节所在的 y 方向上的最小体素坐标
    v_ymax = np.min([height - 1, int(v_center[1] + v_diam) + 5])  # 计算肺结节所在的 y 方向上的最大体素坐标

    v_xrange = range(v_xmin, v_xmax + 1)  # 表示 x 方向上的体素坐标范围
    v_yrange = range(v_ymin, v_ymax + 1)  # 表示 y 方向上的体素坐标范围

    x_data = [x * spacing[0] + origin[0] for x in range(width)]  # 通过列表推导式，创建一个列表 x_data，其中包含了图像每个像素点在 x 方向上的世界坐标
    y_data = [x * spacing[1] + origin[1] for x in range(height)]  # 通过列表推导式，创建一个列表 y_data，其中包含了图像每个像素点在 y 方向上的世界坐标

    # 结节周围全都填充1，用于在图像中标记肺结节位置
    for v_x in v_xrange:
        for v_y in v_yrange:
            p_x = spacing[0] * v_x + origin[0]  # 计算了当前体素在世界坐标系中的 x 坐标
            p_y = spacing[1] * v_y + origin[1]  # 计算了当前体素在世界坐标系中的 y 坐标
            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

 5、生成标签文件的函数：

# 美化.xml标签文件格式，提高可读性
def beatau(e, level=0):
    if len(e) > 0:
        e.text = '\n' + '\t' * (level + 1)
        for child in e:
            beatau(child, level + 1)
        child.tail = child.tail[:-1]
    e.tail = '\n' + '\t' * level

# 生成.xml标签文件的主要函数
def csvtoxml(name, x, y, w, h):
    root = Element('annotation')  # 根节点
    erow1 = Element('folder')  # 节点1
    erow1.text = "VOC"

    erow2 = Element('filename')  # 节点2
    erow2.text = str(name)

    erow3 = Element('size')  # 节点3
    erow31 = Element('width')
    erow31.text = "512"
    erow32 = Element('height')
    erow32.text = "512"
    erow33 = Element('depth')
    erow33.text = "3"
    erow3.append(erow31)
    erow3.append(erow32)
    erow3.append(erow33)

    erow4 = Element('object')
    erow41 = Element('name')
    erow41.text = 'nodule'
    erow42 = Element('bndbox')
    erow4.append(erow41)
    erow4.append(erow42)
    erow421 = Element('xmin')
    erow421.text = str(x - np.round(w / 2).astype(int))
    erow422 = Element('ymin')
    erow422.text = str(y - np.round(h / 2).astype(int))
    erow423 = Element('xmax')
    erow423.text = str(x + np.round(w / 2).astype(int))
    erow424 = Element('ymax')
    erow424.text = str(y + np.round(h / 2).astype(int))
    erow42.append(erow421)
    erow42.append(erow422)
    erow42.append(erow423)
    erow42.append(erow424)

    root.append(erow1)
    root.append(erow2)
    root.append(erow3)
    root.append(erow4)
    beatau(root)

    return ElementTree(root)

6、主要处理函数：提取图像（.npy）、生成掩码（.npy）、生成训练用标签（.xml）

# 历遍所有子集，读取.mhd文件
def process_subset(subset_path):
    # 获取文件列表
    file_list = glob(os.path.join(subset_path, "*.mhd"))

    # 获取文件名与案例相关的路径
    def get_filename(case):
        nonlocal file_list
        for f in file_list:
            if case in f:
                return f

    # 读取包含结节注释的CSV文件
    df_node = pd.read_csv(annotation_path)
    df_node["file"] = df_node["seriesuid"].apply(get_filename)  # 对于 "seriesuid" 中的每个值，找到对应的文件路径
    df_node = df_node.dropna()  # 从数据框中删除包含 NaN 的行，确保数据框中不包含缺失文件路径的记录

    # 遍历所有图像文件
    for fcount, img_file in enumerate(tqdm(file_list)):
        print("获取图像文件 %s 的掩膜" % img_file.replace(subset_path, ""))
        mini_df = df_node[df_node["file"] == img_file]  # mini_df 包含了与当前图像文件相关联的所有结节注释
        if len(mini_df) > 0:  # 跳过没有结节的文件
            # 读取图像数据
            itk_img = sitk.ReadImage(img_file)
            img_array = sitk.GetArrayFromImage(itk_img)  # 将SimpleITK图像对象转换为NumPy三维数组(z, y, x)
            num_z, height, width = img_array.shape  # 图像的深度（z方向上的切片数量）、高度和宽度
            origin = np.array(itk_img.GetOrigin())  # 世界坐标系下的 x, y, z (mm)
            spacing = np.array(itk_img.GetSpacing())  # 世界坐标中的体素间隔 (mm)

            # 遍历所有结节的注释
            for node_idx, cur_row in mini_df.iterrows():
                node_x = cur_row["coordX"]
                node_y = cur_row["coordY"]
                node_z = cur_row["coordZ"]
                diam = cur_row["diameter_mm"]

                # 保留三个切片
                imgs = np.ndarray([3, height, width], dtype=np.float32)  # 创建一个形状为 (3, height, width) 的三维数组 imgs，用于存储图像数据
                masks = np.ndarray([3, height, width], dtype=np.uint8)  # 创建一个形状为 (3, height, width) 的三维数组 masks，用于存储掩膜数据
                center = np.array([node_x, node_y, node_z])  # 表示结节在图像中的中心位置
                v_center = np.rint((center - origin) / spacing)  # 体素坐标系的结节中心 (x, y, z)

                for i, i_z in enumerate(np.arange(int(v_center[2]) - 1, int(v_center[2]) + 2).clip(0, num_z - 1)):  # 防止超出 z
                    # 用于生成.xml标签labels的值
                    name = "%04d_%04d.jpg" % (fcount, node_idx)
                    x = v_center[0]
                    y = v_center[1]
                    w = int(diam / spacing[0] + 5)  # 计算肺结节直径在体素坐标系中的大小
                    h = int(diam / spacing[0] + 5)  # 计算肺结节直径在体素坐标系中的大小
                    label = csvtoxml(name, x, y, w, h)
                    # 生成掩码（.npy）
                    mask = make_mask(center, diam, i_z * spacing[2] + origin[2], width, height, spacing, origin)  # 调用 make_mask 函数生成结节的掩膜
                    masks[i] = mask  # 将生成的掩膜保存到 masks 数组的相应位置
                    # 生成图片（.npy）
                    imgs[i] = normalizePlanes(img_array[i_z])  # 将规范化后的数据保存到 imgs 数组的相应位置

                # 拼接输出路径
                output_npy_path = os.path.join(output_npy_base_path, os.path.basename(subset_path))
                if not os.path.exists(output_npy_path):
                    os.makedirs(output_npy_path)
                output_data_path = os.path.join(output_data_base_path, os.path.basename(subset_path), "labels")
                if not os.path.exists(output_data_path):
                    os.makedirs(output_data_path)

                np.save(os.path.join(output_npy_path, "images_%04d_%04d.npy" % (fcount, node_idx)), imgs)  # 将处理后的图像数据 imgs 保存为一个 .npy 文件
                np.save(os.path.join(output_npy_path, "masks_%04d_%04d.npy" % (fcount, node_idx)), masks)  # 将处理后的结节掩膜 masks 保存为一个 .npy 文件
                label.write(os.path.join(output_data_path, "%04d_%04d_1.xml" % (fcount, node_idx)))
                label.write(os.path.join(output_data_path, "%04d_%04d_2.xml" % (fcount, node_idx)))
                label.write(os.path.join(output_data_path, "%04d_%04d_3.xml" % (fcount, node_idx)))

7、获取肺实质函数：

# 获取肺实质
def get_lungmask(img_list):
    for img_file in img_list:
        # 加载图像数据
        imgs_to_process = np.load(img_file).astype(np.float64)
        print("处理图像", img_file)

        # 处理每个切片
        for i in range(len(imgs_to_process)):
            # print(imgs_to_process.shape)
            img = imgs_to_process[i]

            # 标准化
            mean = np.mean(img)
            std = np.std(img)
            img = (img - mean) / std

            # 寻找肺部附近的平均像素，以重新调整过度曝光的图像
            middle = img[100:400, 100:400]

            # 使用 Kmeans 算法将前景（放射性不透明组织）和背景（放射性透明组织，即肺部）分离。
            # 仅在图像中心进行此操作，以尽可能避免图像的非组织部分。
            kmeans = KMeans(n_clusters=2, n_init=10).fit(np.reshape(middle, [np.prod(middle.shape), 1]))
            centers = sorted(kmeans.cluster_centers_.flatten())
            threshold = np.mean(centers)
            thresh_img = np.where(img < threshold, 1.0, 0.0)  # 阈值化图像，二值化处理

            # 腐蚀和膨胀
            eroded = morphology.erosion(thresh_img, np.ones([4, 4]))
            dilation = morphology.dilation(eroded, np.ones([10, 10]))

            # 对二值图像进行标记，标记连通区域
            labels = measure.label(dilation)
            label_vals = np.unique(labels)
            regions = measure.regionprops(labels)

            good_labels = []
            for prop in regions:
                B = prop.bbox  # 边界框
                if B[2] - B[0] < 475 and B[3] - B[1] < 475 and B[0] > 40 and B[2] < 472:
                    good_labels.append(prop.label)

            # 创建肺部掩码lungmask
            lungmask = np.zeros([512, 512], dtype=np.int8)
            for N in good_labels:
                lungmask = lungmask + np.where(labels == N, 1, 0)
                # print(lungmask.shape)
            lungmask = morphology.dilation(lungmask, np.ones([10, 10]))
            imgs_to_process[i] = lungmask

        # 将处理后的 lungmask 保存
        # print(imgs_to_process.shape)
        np.save(img_file.replace("images", "lungmasks"), imgs_to_process)

8、生成训练用图片函数：分割肺实质后保存成.jpg图片

def npytojpg(img_list):
    for i in range(len(img_list)):
        # 加载原始图像和肺部 mask
        origin_img = np.load(img_list[i])
        lung_mask = np.load(img_list[i].replace("images", "lungmasks"))
        # 使用os.path.basename()获取文件名
        file_name_npy = os.path.basename(img_list[i])
        # 使用os.path.splitext()分割文件名和后缀
        file_name, file_extension = os.path.splitext(file_name_npy)

        for j in range(len(origin_img)):
            # 获取当前切片的原始图像和肺部 mask
            img = origin_img[j]
            lung_img = lung_mask[j]

            # 对原始图像进行肺部分割
            segment = img * lung_img

            # 将分割结果转为 PIL Image，并保存为 JPEG 文件
            segment = Image.fromarray(segment)
            segment = segment.convert('RGB')
            new_file_name = file_name.replace("images_", "")
            new_file_path = os.path.join(save_path, f"{new_file_name}_{j + 1}.jpg")
            segment.save(new_file_path)

9、依次调用以上函数：

# 遍历所有子集，生成images*.npy、masks*.npy和*.xml文件
for subset_number in range(10):
    subset_path = os.path.join(luna_base_path, f"subset{subset_number}")
    process_subset(subset_path)

# 遍历所有子集，生成lungmasks*.npy文件
for i in range(10):
    img_list = glob(os.path.join(output_npy_base_path, f"subset{i}/", "images*.npy"))
    get_lungmask(img_list)

# 历遍所有子集，生成*.jpg文件
for subset_num in range(10):
    # 获取文件列表
    img_list = glob(os.path.join(output_npy_base_path, "subset{}".format(subset_num), "images*.npy"))
    print(img_list)
    # 指定保存路径
    save_path = os.path.join(output_data_base_path, "subset{}/".format(subset_num), "images")
    # 如果保存路径不存在，则创建
    if not os.path.exists(save_path):
        os.makedirs(save_path)
    npytojpg(img_list)