LUNA16三维CT肺结节检测的数据预处理与PyTorch加载链路

简介：本资源是一套基于Python与PyTorch实现的3D CT肺结节检测完整项目，面向人工智能、医学影像、生物信息等方向的高校学生、科研人员及初学者，聚焦医学图像分析中的关键任务——肺部小结节自动识别与定位。项目以国际公开LUNA16数据集为基准，涵盖数据预处理（dicom2raw、full_prep）、候选区域检测（net_detector、test_detect）、分类判别（net_classifier、test_classifier）及结果可视化全流程，配套设计文档、训练/提交配置脚本及Jupyter演示笔记，结构清晰、模块解耦度高。压缩包共53个文件，含38个核心Python源码（含网络定义、数据加载、训练逻辑）、4个CSV标注文件、2个Numpy中间数据、2个PNG效果示意图及Shell运行脚本等，整体9.62MB，轻量易部署。目前已有96人学习下载，提供可直接运行的完整代码、详细注释与分步说明，支持毕设开题、课程设计实践及深度学习进阶训练，亦可作为医学AI入门项目的优质参考范例。

1. 这不是“跑通一个模型”那么简单：3D CT肺结节检测真正卡在数据预处理、三维卷积建模与LUNA16标注一致性上

你下载了医学影像分析-3DCT肺结节检测（Python+PyTorch+LUNA16数据集）.zip，解压后看到一堆.mhd/.raw文件和annotations.csv，trainset.py里写着torch.nn.Conv3d，但train.py一运行就报RuntimeError: expected scalar type Float but found Byte——这不是环境没装好，而是 LUNA16 的原始 DICOM 转 MHD 流程中，窗宽窗位（WW/WL）未归一化、体素间距未重采样、结节标注半径未按实际 mm 映射到体素坐标，导致输入张量分布崩坏、GT mask 尺寸错位。本篇不讲“如何安装 PyTorch”，而是聚焦LUNA16 数据集在 PyTorch 中的三维可复现加载链路：从原始.mhd文件解析、HU 值截断与线性归一化、各向异性体素重采样、结节中心点与直径的体素空间对齐，到最终构建torch.utils.data.Dataset返回(C, D, H, W)张量与(D, H, W)二值 mask。适合已配好python 3.10+和pytorch 2.0+环境、能跑通 MNIST 但第一次处理医学三维影像的工程师；也适合熟悉 PyTorch 训练循环、却在 LUNA16 上反复卡在DataLoader返回 shape 不一致的老手。

2. 解析 LUNA16 的 MHD+RAW 结构：用 SimpleITK 读取并校验体素物理尺寸与 HU 值范围

LUNA16 数据集以.mhd（元数据头文件）+.raw（二进制体素数据）形式组织，其关键元信息不在 CSV 标注里，而在.mhd文件中。直接np.fromfile()会丢失方向、间距、原点等物理坐标信息，导致后续重采样和标注映射完全失效。必须用SimpleITK（非nibabel或pydicom）解析，因其完整支持 LUNA16 官方转换脚本生成的 MHD 格式，并能自动处理TransformMatrix和ElementSpacing。

2.1 读取单例 CT 并提取关键元数据

import SimpleITK as sitk import numpy as np def load_ct_scan(mhd_path): # 用 SimpleITK 读取，保留全部元数据 itk_img = sitk.ReadImage(mhd_path) # 获取原始体素数组（注意：sitk.GetArrayFromImage 返回 (D, H, W)，Z轴在前） img_array = sitk.GetArrayFromImage(itk_img).astype(np.float32) # 提取物理参数（单位：mm） spacing = np.array(itk_img.GetSpacing()) # [z_spacing, y_spacing, x_spacing] origin = np.array(itk_img.GetOrigin()) # [z_origin, y_origin, x_origin] direction = np.array(itk_img.GetDirection()) # 9x1 向量，描述坐标系旋转 print(f"原始形状: {img_array.shape}, 体素间距: {spacing}, 原点: {origin}") print(f"HU 值范围: [{img_array.min():.1f}, {img_array.max():.1f}]") return img_array, spacing, origin, direction # 示例调用 img, spc, org, dir_ = load_ct_scan("subset0/1.3.6.1.4.1.14519.5.2.1.6279.6001.105756658031515062000744821260.mhd")

提示：LUNA16 的.mhd中ElementSpacing通常为[1.0, 0.585938, 0.585938]（即 Z 轴 1mm，XY 轴约 0.586mm），但部分病例存在0.78125或0.625，必须逐例读取，不可硬编码。img_array.min()常为-1024（空气），max()可达3071（金属伪影），但有效肺实质 HU 区间为[-1000, 400]，需截断。

2.2 HU 值标准化：从原始整型到网络可用浮点张量

CT 图像本质是 Hounsfield Unit（HU），其物理意义明确：水为 0，空气为 -1000，骨约为 +1000。深度学习模型无法直接处理[-1000, 3000]的宽动态范围，且 LUNA16 中存在超出范围的噪声值（如-3000）。标准做法是截断 + 线性缩放：

def normalize_hu(image_array, min_hu=-1000, max_hu=400): """将 HU 值截断并归一化到 [0, 1]""" image_array = np.clip(image_array, min_hu, max_hu) # 截断 image_array = (image_array - min_hu) / (max_hu - min_hu) # 线性归一化 return image_array.astype(np.float32) # 应用 normalized_img = normalize_hu(img) # shape: (D, H, W), dtype: float32 print(f"归一化后范围: [{normalized_img.min():.3f}, {normalized_img.max():.3f}]") # 恒为 [0.0, 1.0]

2.2.1 为什么不用 Z-score？

LUNA16 各病例 HU 分布差异极大（如含金属植入物的病例均值偏高），全局mean/std会导致不同样本归一化失真。min/max截断基于解剖学先验（肺组织 HU ∈ [-1000, 400]），是领域共识，PyTorch 官方医学教程（如 MONAI）亦采用此法。

2.2.2 关键参数表：LUNA16 典型 HU 截断阈值

3. 构建三维训练样本：重采样、结节标注映射与 Patch 提取

LUNA16 的标注文件annotations.csv提供结节中心(x, y, z)（世界坐标，单位 mm）和直径d（单位 mm）。但img_array是体素坐标，二者需通过origin和spacing转换。若跳过此步直接用 CSV 中的(x,y,z)作索引，会导致结节中心偏移 3~5 个体素——对 3D CNN 而言，这足以让正样本丢失。

3.1 将世界坐标（mm）映射到体素坐标（index）

def world_to_voxel(world_coord, origin, spacing, direction=None): """ 将世界坐标 (x,y,z) mm 转为体素坐标 (i,j,k) index 注意：LUNA16 .mhd 中 world_coord = (x,y,z)，但 img_array 索引为 (k,j,i) 即 (z,y,x) """ # world_coord 是 [x, y, z]，origin 是 [x0, y0, z0]，spacing 是 [sx, sy, sz] # 公式：voxel_index = (world_coord - origin) / spacing voxel_coord = (np.array(world_coord) - np.array(origin)) / np.array(spacing) # 四舍五入取整，确保为整数索引 return np.round(voxel_coord).astype(int)[[2, 1, 0]] # 调换顺序：(z,y,x) -> (k,j,i) # 示例：从 annotations.csv 读取一行 # 1.3.6.1.4.1.14519.5.2.1.6279.6001.105756658031515062000744821260,123.45,67.89,101.23,6.5 seriesuid, x, y, z, diameter = "1.3.6.1.4.1.14519.5.2.1.6279.6001.105756658031515062000744821260", 123.45, 67.89, 101.23, 6.5 world_center = [float(x), float(y), float(z)] # [x,y,z] in mm voxel_center = world_to_voxel(world_center, org, spc) # 返回 [k,j,i] 即 [z_idx, y_idx, x_idx] print(f"世界中心 {world_center} -> 体素中心 {voxel_center}")

注意：world_to_voxel输出[k,j,i]，而img_array[k,j,i]才是该点灰度值。若误用[i,j,k]，则坐标全错。

3.2 各向异性体素重采样：统一为 1mm³ 立方体素

LUNA16 原始体素间距不一致（如[1.0, 0.585938, 0.585938]），直接输入 3D CNN 会导致 Z 轴感受野远大于 XY 轴。必须重采样为各向同性（如1.0mm）。使用SimpleITK的ResampleImageFilter：

def resample_image(itk_img, new_spacing=[1.0, 1.0, 1.0], interpolator=sitk.sitkLinear): original_spacing = itk_img.GetSpacing() original_size = itk_img.GetSize() # 计算新尺寸：new_size = old_size * old_spacing / new_spacing new_size = [ int(round(original_size[0] * original_spacing[0] / new_spacing[0])), int(round(original_size[1] * original_spacing[1] / new_spacing[1])), int(round(original_size[2] * original_spacing[2] / new_spacing[2])) ] resample = sitk.ResampleImageFilter() resample.SetOutputSpacing(new_spacing) resample.SetSize(new_size) resample.SetOutputDirection(itk_img.GetDirection()) resample.SetOutputOrigin(itk_img.GetOrigin()) resample.SetInterpolator(interpolator) resample.SetDefaultPixelValue(itk_img.GetPixelIDValue()) return resample.Execute(itk_img) # 应用 resampled_itk = resample_image(itk_img, new_spacing=[1.0, 1.0, 1.0]) resampled_array = sitk.GetArrayFromImage(resampled_itk).astype(np.float32) resampled_array = normalize_hu(resampled_array) # 重采样后需重新归一化 print(f"重采样后形状: {resampled_array.shape}") # 如 (512, 512, 512) 或类似立方体

3.3 构建训练 Patch：以结节为中心裁剪 64×64×64 体素块

为平衡显存与感受野，主流做法是提取固定尺寸 Patch（如64×64×64），而非全图输入。需确保：

• 正样本：Patch 中心为结节中心，mask 中对应位置为 1；
• 负样本：从无结节区域随机采样，mask 全 0。

def get_patch_and_mask(image_array, voxel_center, patch_size=64, radius_factor=1.0): """ 提取以 voxel_center 为中心的 patch，及对应二值 mask radius_factor: 结节半径放大倍数，用于扩大正样本区域（如 1.5 倍） """ d, h, w = image_array.shape z, y, x = voxel_center # 计算 patch 边界（确保不越界） z_min = max(0, z - patch_size//2) z_max = min(d, z + patch_size//2) y_min = max(0, y - patch_size//2) y_max = min(h, y + patch_size//2) x_min = max(0, x - patch_size//2) x_max = min(w, x + patch_size//2) # 提取图像 patch patch = image_array[z_min:z_max, y_min:y_max, x_min:x_max] # 补零至目标尺寸 pad_z = patch_size - (z_max - z_min) pad_y = patch_size - (y_max - y_min) pad_x = patch_size - (x_max - x_min) patch = np.pad(patch, ((0, pad_z), (0, pad_y), (0, pad_x)), mode='constant', constant_values=0) # 构建 mask：以 voxel_center 为球心，半径 = diameter/2 * radius_factor mask = np.zeros((patch_size, patch_size, patch_size), dtype=np.float32) # 计算 patch 内结节中心偏移 center_in_patch = (patch_size//2, patch_size//2, patch_size//2) # 生成球形 mask（简化为立方体也可，但球形更符合解剖） z_grid, y_grid, x_grid = np.ogrid[:patch_size, :patch_size, :patch_size] dist_sq = (z_grid - center_in_patch[0])**2 + \ (y_grid - center_in_patch[1])**2 + \ (x_grid - center_in_patch[2])**2 radius_voxel = int(round(diameter / 2.0 / 1.0 * radius_factor)) # 直径转体素半径，假设重采样后 spacing=1.0 mask[dist_sq <= radius_voxel**2] = 1.0 return patch, mask # 示例 patch, mask = get_patch_and_mask(resampled_array, voxel_center, patch_size=64, radius_factor=1.5) print(f"Patch shape: {patch.shape}, Mask sum: {mask.sum():.0f}") # 正样本 mask.sum() > 0

4. PyTorch Dataset 实现：支持多结节、负样本平衡与内存优化

LUNA16 每例 CT 可含多个结节（annotations.csv中同一seriesuid多行），且负样本数量远超正样本。Dataset 必须：

• 支持单例多结节迭代；
• 按 1:1 或 1:3 比例混合正负样本；
• 避免一次性加载所有.raw到内存（LUNA16 单例可达 500MB）。

4.1 构建结节级索引表（非病例级）

import pandas as pd from pathlib import Path def build_nodule_index(annotations_csv, subset_dirs): """ 构建 (seriesuid, x, y, z, diameter) 索引列表，每行一个结节 """ ann_df = pd.read_csv(annotations_csv) nodule_list = [] for _, row in ann_df.iterrows(): seriesuid = row['seriesuid'] # 查找对应 .mhd 文件路径 mhd_path = None for subset_dir in subset_dirs: p = Path(subset_dir) / f"{seriesuid}.mhd" if p.exists(): mhd_path = str(p) break if mhd_path: nodule_list.append({ 'mhd_path': mhd_path, 'x': row['coordX'], 'y': row['coordY'], 'z': row['coordZ'], 'diameter': row['diameter_mm'] }) return nodule_list # 使用 subset_dirs = ["subset0", "subset1"] # 实际需列出全部 subset0-9 nodule_index = build_nodule_index("annotations.csv", subset_dirs) print(f"共 {len(nodule_index)} 个结节标注")

4.2 自定义 Dataset：惰性加载 + 负样本在线采样

import torch from torch.utils.data import Dataset class LUNADataset(Dataset): def __init__(self, nodule_index, patch_size=64, neg_ratio=3, transform=None): self.nodule_index = nodule_index self.patch_size = patch_size self.neg_ratio = neg_ratio # 每个正样本配几个负样本 self.transform = transform # 正样本索引（每个结节一行） self.positive_samples = [(i, 'pos') for i in range(len(nodule_index))] # 负样本索引：随机选病例，再随机选无结节位置 self.negative_samples = [(i, 'neg') for i in range(len(nodule_index) * neg_ratio)] self.all_samples = self.positive_samples + self.negative_samples def __len__(self): return len(self.all_samples) def __getitem__(self, idx): sample_type = self.all_samples[idx][1] if sample_type == 'pos': # 加载正样本 nodule_info = self.nodule_index[self.all_samples[idx][0]] itk_img = sitk.ReadImage(nodule_info['mhd_path']) img_array = sitk.GetArrayFromImage(itk_img).astype(np.float32) img_array = normalize_hu(img_array) # 重采样（此处为演示，实际应缓存或预处理） resampled_itk = resample_image(itk_img, new_spacing=[1.0, 1.0, 1.0]) resampled_array = sitk.GetArrayFromImage(resampled_itk).astype(np.float32) resampled_array = normalize_hu(resampled_array) # 映射中心点 world_center = [nodule_info['x'], nodule_info['y'], nodule_info['z']] spacing = np.array(resampled_itk.GetSpacing()) origin = np.array(resampled_itk.GetOrigin()) voxel_center = world_to_voxel(world_center, origin, spacing) # 提取 patch 和 mask patch, mask = get_patch_and_mask( resampled_array, voxel_center, patch_size=self.patch_size, radius_factor=1.5 ) else: # 'neg'：随机采样负样本 # 随机选一个病例 rand_idx = np.random.randint(0, len(self.nodule_index)) nodule_info = self.nodule_index[rand_idx] itk_img = sitk.ReadImage(nodule_info['mhd_path']) resampled_itk = resample_image(itk_img, new_spacing=[1.0, 1.0, 1.0]) resampled_array = sitk.GetArrayFromImage(resampled_itk).astype(np.float32) resampled_array = normalize_hu(resampled_array) # 随机选一个无结节位置（避开所有已知结节） d, h, w = resampled_array.shape # 简化：随机选中心，确保距所有结节 > 3*radius z, y, x = np.random.randint(32, d-32), np.random.randint(32, h-32), np.random.randint(32, w-32) patch = resampled_array[ z-32:z+32, y-32:y+32, x-32:x+32 ] mask = np.zeros((64,64,64), dtype=np.float32) # 转为 torch tensor，添加通道维度 image_tensor = torch.from_numpy(patch).unsqueeze(0) # (1, 64, 64, 64) mask_tensor = torch.from_numpy(mask).unsqueeze(0) # (1, 64, 64, 64) if self.transform: image_tensor = self.transform(image_tensor) return image_tensor, mask_tensor # 实例化 dataset = LUNADataset(nodule_index, patch_size=64, neg_ratio=3) loader = torch.utils.data.DataLoader(dataset, batch_size=4, shuffle=True, num_workers=4)

提示：实际训练中，resample_image和world_to_voxel计算开销大，应预处理并缓存重采样后的.nii.gz文件，或使用MONAI CacheDataset。此处为展示逻辑完整性。

5. 验证数据加载正确性：可视化 Patch、Mask 与体素坐标对齐

在train.py开头加入验证逻辑，避免训 10 小时才发现 mask 全黑或中心偏移。核心是用 matplotlib 画出中间层 Z 切片，叠加 mask 轮廓：

import matplotlib.pyplot as plt def visualize_patch_and_mask(image_tensor, mask_tensor, save_path=None): """ 可视化 (1, D, H, W) 图像与 mask 的中间切片 """ img = image_tensor[0].numpy() # (D, H, W) mask = mask_tensor[0].numpy() # (D, H, W) mid_z = img.shape[0] // 2 fig, axes = plt.subplots(1, 2, figsize=(10, 4)) # 图像切片 axes[0].imshow(img[mid_z], cmap='gray') axes[0].set_title(f'Image Slice Z={mid_z}') axes[0].axis('off') # Mask 叠加（红色轮廓） axes[1].imshow(img[mid_z], cmap='gray') # 仅显示 mask 边界（用 skimage.measure.find_contours） from skimage import measure contours = measure.find_contours(mask[mid_z], 0.5) for contour in contours: axes[1].plot(contour[:, 1], contour[:, 0], linewidth=1.5, color='red') axes[1].set_title('Mask Contour Overlay') axes[1].axis('off') if save_path: plt.savefig(save_path, bbox_inches='tight', dpi=150) plt.close() else: plt.show() # 在 DataLoader 循环中验证 for i, (x, y) in enumerate(loader): print(f"Batch {i}: image shape {x.shape}, mask shape {y.shape}") visualize_patch_and_mask(x[0], y[0], save_path=f"debug_batch_{i}.png") if i >= 2: # 只看前3个 batch break

5.1 正确输出特征：

• 左图：CT 切片可见肺纹理、血管、可能的结节高亮区；
• 右图：红色轮廓紧密包裹左图中的高密度区域，且轮廓中心与结节视觉中心一致；
• 若轮廓漂移、破碎或完全不重合，则world_to_voxel或resample_image参数有误。

5.2 常见失败模式与修复指令：

执行python debug_loader.py后，若生成的debug_batch_0.png中红色轮廓精准套住结节，则数据加载链路已打通。此时可安全接入nn.Conv3d模型，进入训练阶段。

本文还有配套的精品资源，点击获取