简介面向医学影像分析与深度学习初学者提供肺结节检测与分类的完整Python实现以Jupyter Notebook交互形式呈现覆盖从CT影像读入到结节分类的算法流程。压缩包共439个文件大小约195.91MB主要包含397个DICOM原始影像、13个Notebook代码文件、CSV标注与元数据表、模型权重及可视化图表功能模块涵盖数据划分、预处理、模型搭建、训练评估与预测验证。已有2415人学习下载。项目经过两周调试前三个代码文件可顺利运行首个文件注释详尽能帮助阅读者快速理解卷积神经网络在医学影像分类中的应用。整体代码结构清晰从LIDC-IDRI数据集处理到结果可视化均有对应实现适合作为算法复现、实验扩展或毕业设计的参考素材也可迁移至其他医学影像分类任务。1. 肺结节分类算法代码先想清楚你要做的是哪一版分类再动手写代码肺结节分类算法代码这个标题听着像是一个可以直接拿来跑的东西但实际动手时会发现找错任务边界才是大多数项目翻车的起点。临床上最常见的肺结节分类是良恶性二分类输入可以是医生标注好的结节区域也可以是从检测算法输出的候选框里截出来的ROI。前者是纯分类问题后者是检测加分类的级联问题两者在数据组织、标签含义和评估方式上差别很大。这篇文章只做一件事把从DICOM原始数据到得到可评估分类结果的每一段代码怎么组织、参数怎么设、坑在哪讲清楚覆盖特征工程加XGBoost和2D/3D深度学习两条主流路线适合刚接触医学影像的算法工程师和正在复现论文的学生。先别急着选模型先把数据管线和评估口径定下来后面才有得聊。2. 预处理是分类精度的第一道闸门DICOM读取、重采样与窗位窗宽任何肺结节分类算法的代码都是从预处理开始的。CT原始DICOM序列并不是一个可以直接丢给神经网络的张量它的像素值是存储值不是物理意义上的CT值层间间距和面内间距不一致不同扫描协议的窗位窗宽也不同。这三项不处理后面喂给模型的数据就自带偏差分类精度上不去时你甚至会怀疑是模型选型的问题。常见做法是先做四道工序读取DICOM并换算HU、重采样到各向同性、按标注中心裁剪ROI、加窗加归一化。每一步都有参数讲究不能只调库不思考。2.1 DICOM体数据读取先把存储值换算成真实的CT值HU读取DICOM的常用库是pydicom和SimpleITK。我一般用pydicom手动读因为需要同时处理RescaleSlope和RescaleIntercept这两个字段手动读不容易漏。CT图像在PACS里存的通常是整数存储值真实CT值亨氏单位HU必须通过存储值 * RescaleSlope RescaleIntercept换算这一步漏了肺窗和纵隔窗显示都不对更别提给模型训练了。import os import glob import numpy as np import pydicom def load_ct_volume(dicom_dir): 读取一个序列的DICOM切片按空间位置排序成3D体数据返回HU数组和间距。 slices [] for dcm_path in sorted(glob.glob(os.path.join(dicom_dir, *.dcm))): ds pydicom.dcmread(dcm_path) # 关键存储值转成真实CT值(HU)不乘RescaleSlope会导致HU整体偏小 hu ds.pixel_array.astype(np.float32) * float(ds.RescaleSlope) float(ds.RescaleIntercept) pos_z float(ds.ImagePositionPatient[2]) slices.append((pos_z, hu)) slices.sort(keylambda item: item[0]) volume np.stack([item[1] for item in slices], axis0) # (Z, H, W) spacing np.array([ float(ds.PixelSpacing[0]), # 行间距单位mm float(ds.PixelSpacing[1]), # 列间距单位mm float(slices[-1][0] - slices[0][0]) / max(len(slices) - 1, 1), # 层间距 ], dtypenp.float32) return volume, spacing这段代码有三个地方需要说明。第一排序用的ImagePositionPatient[2]是切片在病人坐标系里的Z坐标比直接按文件名排序可靠因为扫描方向可能从头到脚也可能反过来。第二pixel_array是uint16原始存储值先转float32再计算避免溢出。第三返回的spacing顺序和volume轴的顺序对应volume是(Z,H,W)spacing对应(层间距, 行间距, 列间距)后面重采样和裁剪都依赖这个约定。如果个别DICOM没有RescaleSlope字段代码会抛异常处理方法是在循环里用getattr(ds, RescaleSlope, 1)兜底。2.2 重采样与ROI裁剪把不同扫描参数的结节统一到同一个坐标网格不同医院的CT扫描参数差异很大层间距有的0.5mm有的2mm面内分辨率也不一样。如果不做重采样同一物理尺寸的结节在不同样本里占的体素数不同卷积核的感受野覆盖范围也跟着变。常见做法是用scipy.ndimage的zoom函数把体数据线性重采样到各向同性1mm。插值阶数方面HU体数据用三阶样条order3但后续如果涉及mask标注一定要用最近邻否则边界会被磨掉。import numpy as np from scipy.ndimage import zoom def resample_volume(volume, orig_spacing, target_spacing(1.0, 1.0, 1.0)): 把非各向同性CT体数据重采样到目标间距默认1x1x1 mm。 factors np.array(orig_spacing, dtypenp.float32) / np.array(target_spacing, dtypenp.float32) return zoom(volume, factors, order3)重采样后再裁剪就省心很多因为三个方向的间距都是1mm不用再分别算半宽。以标注中心为原点裁固定物理尺寸比如48mm那么体素尺寸就是48x48x48def crop_roi(volume, center_voxel, size_mm48): 按标注中心裁剪固定物理尺寸的ROI适用于重采样到1mm后的数据。 half size_mm // 2 z, y, x [int(round(c)) for c in center_voxel] z0, z1 max(z - half, 0), min(z half, volume.shape[0]) y0, y1 max(y - half, 0), min(y half, volume.shape[1]) x0, x1 max(x - half, 0), min(x half, volume.shape[2]) roi volume[z0:z1, y0:y1, x0:x1] # 边缘不足时用最小像素值填充避免直接裁回小尺寸 if roi.shape ! (size_mm, size_mm, size_mm): pad [(half, half), (half, half), (half, half)] roi np.pad(roi, pad, modeconstant, constant_values-1000)[ z-half:zhalf, y-half:yhalf, x-half:xhalf ] return roi这段话里的重点在最后那个填充分支。靠近肺边缘或扫描边缘的结节裁剪窗口可能超出体数据范围直接切片会让样本之间尺寸不一致DataLoader没法打包。用-1000填充是医学影像里的惯例因为空气的HU值就是-1000。注意center_voxel是在重采样后的坐标空间如果是原始坐标需要先按重采样倍率换算或者干脆先重采样完再用代码里的坐标变换函数统一处理这也是为什么我建议把重采样放在ROI裁剪之前。2.3 窗位窗宽与归一化你看到的图和你喂给网络的数不是一回事CT值是16bit范围的整数直接归一化到[0,1]会损失软组织对比度。临床上医生看肺部结节用的是肺窗窗位约-600HU窗宽约1500HU也就是只保留-1350到150HU这个区间。算法代码里同样应该先加窗再归一化这比全局min-max归一化更符合读片逻辑。窗位窗宽的选择本身就是调参点找磨玻璃结节GGO和找实性结节合适的窗不完全一样。def apply_window(volume, center-600, width1500): 肺窗窗位窗宽调整常见窗口肺窗(-600,1500)纵隔窗(40,400)。 lower center - width / 2.0 upper center width / 2.0 return np.clip(volume, lower, upper) def normalize_to_unit_range(volume): 把加窗后的数据线性归一化到[0,1]避免模型输入数值范围不稳定。 vmin, vmax float(volume.min()), float(volume.max()) return (volume - vmin) / (vmax - vmin 1e-6) # 完整预处理链路的组合调用 def preprocess_pipeline(dicom_dir, center_voxel): vol, spacing load_ct_volume(dicom_dir) vol resample_volume(vol, spacing) # 1mm各向同性 roi crop_roi(vol, center_voxel, size_mm48) # 48mm立方体 roi apply_window(roi, center-600, width1500) # 肺窗 roi normalize_to_unit_range(roi) # [0,1] return roi.astype(np.float32)我在数据集实验里对比过两组设置一组只做min-max归一化另一组先加肺窗再归一化在同一个XGBoost特征管线里AUC大约有3到5个百分点的差距。这里想强调的不是绝对值而是加窗这个操作相当于给模型提供了先验知识告诉它哪些灰度范围是有效的。对于要做2D切片的项目一般从ROI体积里取中间切片或最大面积切片我个人习惯取中间切片配三视图横断位、冠状位、矢状位做三通道输入这个在第四章会展开。3. 特征工程加XGBoost小样本场景下最稳的肺结节分类方案很多入门项目一上来就堆3D卷积网络但肺结节标注样本通常只有几百到一两千个深度学习在这个量级很容易过拟合。特征工程加XGBoost是数据量不足时更稳的一条路也是我建议先做的基线方案。基线跑通后深度学习的增量才有对比意义。这一章的路线是从预处理后的ROI里提取形态和纹理特征组织成表格数据再用XGBoost做二分类最后用AUC和敏感度评估。整段代码可以直接跑通但需要说明的是特征提取的参数直接决定分类上限。3.1 手工特征提取形态、直方图、GLCM和LBP怎么落到代码里肺结节的特征可以分成三类形态特征直径、体积、球度、边缘毛刺、灰度统计特征均值、方差、偏度、峰度、纹理特征GLCM、LBP、HOG。医学影像社区里有个pyradiomics库可以直接提上千个特征但这里我用更透明的手工提取方式方便排查特征含义和复现。注意GLCM和LBP都是基于2D切片的所以先取ROI的中间切片再转成uint8。from skimage.feature import local_binary_pattern, graycomatrix, graycoprops def extract_texture_features(roi_2d): 从2D结节切片提取LBP和GLCM纹理特征返回dict。 feats {} # 转uint8GLCM和LBP都要求输入是整数灰度图 img8 np.clip(roi_2d * 255.0, 0, 255).astype(np.uint8) # LBP局部二值模式uniform模式P8,R1输出9维直方图 lbp local_binary_pattern(img8, P8, R1, methoduniform) lbp_hist, _ np.histogram(lbp, binslist(range(10)), densityTrue) for i, v in enumerate(lbp_hist): feats[flbp_bin_{i}] float(v) # GLCM灰度共生矩阵统计距离1和3角度0/45/90/135度 glcm graycomatrix( img8, distances[1, 3], angles[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4], levels256, symmetricTrue, normedTrue, ) feats[glcm_contrast] float(graycoprops(glcm, contrast).mean()) feats[glcm_homogeneity] float(graycoprops(glcm, homogeneity).mean()) feats[glcm_energy] float(graycoprops(glcm, energy).mean()) feats[glcm_correlation] float(graycoprops(glcm, correlation).mean()) return feats def extract_shape_features(roi_mask, voxel_size_mm1.0): 从二值mask提取形态特征假设mask来自结节分割标注。 feats {} n_voxels roi_mask.sum() feats[volume_mm3] float(n_voxels * voxel_size_mm ** 3) feats[diameter_max_mm] float(np.max(roi_mask.sum(axis2).sum(axis1))) * voxel_size_mm feats[sphericity] float(3.0 * (n_voxels * 4.19) ** (2.0 / 3.0) / (roi_mask.sum(axis2).sum() * 6.0 1e-6)) return feats这里有几个参数值得说明。LBP的P和R分别是邻域点数和半径(8,1)是纹理识别的常用起点想捕获更粗的纹理可以把R调到2但特征维度不变bin数量只受uniform模式影响。GLCM的distances[1,3]表示统计1像素和3像素间距的共生关系前者刻画细纹理后者刻画粗纹理levels256是标准灰度级数。要注意的是这些特征如果直接扔给XGBoost值域差异大体积几百上千LBP密度零点几XGBoost对量纲不敏感但如果是给线性模型或SVM用必须先标准化。3.2 用XGBoost做良恶性分类一版可以直接跑通的最小训练代码特征矩阵准备好后训练部分用XGBoost。这里重点讲三个参数scale_pos_weight处理类别不平衡eval_metricauc让早停基于AUC而非准确率early_stopping_rounds防止过拟合。标签按通常约定恶性为1良性为0。import numpy as np import xgboost as xgb from sklearn.model_selection import StratifiedKFold from sklearn.metrics import roc_auc_score, confusion_matrix def train_xgboost_nodule_classifier(features, labels): 使用5折分层交叉验证训练XGBoost肺结节良恶性分类器。 skf StratifiedKFold(n_splits5, shuffleTrue, random_state42) oof_pred np.zeros(len(labels)) auc_list [] for fold, (tr_idx, va_idx) in enumerate(skf.split(features, labels)): X_tr, X_va features[tr_idx], features[va_idx] y_tr, y_va labels[tr_idx], labels[va_idx] # scale_pos_weight的意义是负样本/正样本数量比让模型关注少数类(恶性) pos_weight float(np.sum(y_tr 0)) / max(np.sum(y_tr 1), 1) model xgb.XGBClassifier( n_estimators200, max_depth4, # 深度4-6足够过深在小样本上容易过拟合 learning_rate0.05, subsample0.8, colsample_bytree0.8, scale_pos_weightpos_weight, eval_metricauc, use_label_encoderFalse, n_jobs8, ) model.fit( X_tr, y_tr, eval_set[(X_va, y_va)], early_stopping_rounds20, verboseFalse, ) oof_pred[va_idx] model.predict_proba(X_va)[:, 1] auc_list.append(roc_auc_score(y_va, oof_pred[va_idx])) print(f5-fold AUC: {np.mean(auc_list):.4f} ± {np.std(auc_list):.4f}) return model, oof_pred这段示例代码的要点在eval_metricauc和early_stopping_rounds20的配合验证集AUC连续20轮不涨就停返回的是验证集上最优的那棵树。scale_pos_weight不是拍脑袋设的它等于负样本数除以正样本数代码里用训练集动态算每个fold略有差异这是刻意为之——避免把测试集信息引入训练。max_depth4对这个量级的特征维度是够的我见过有人把深度调到10训练集AUC到0.99验证集只有0.7典型的过拟合。3.3 类别不平衡下的评估只看准确率会被它骗得很惨肺结节的良恶性分布天然不平衡公开数据集里良性比例通常高于恶性自己从医院拿的数据更夸张。如果直接用准确率评估模型把所有结节都判成良性准确率也可能有70%以上看起来不差实际毫无临床价值。正确的评估口径是AUC、敏感度、特异度和F1。下面的代码在OOF预测结果上同时输出这四个指标以及按约登指数选出的最优阈值。from sklearn.metrics import roc_curve def evaluate_nodule_classifier(y_true, y_pred_proba): 输出AUC、敏感度、特异度以及约登指数对应的最优阈值。 auc roc_auc_score(y_true, y_pred_proba) fpr, tpr, thresholds roc_curve(y_true, y_pred_proba) # 约登指数 敏感度 特异度 - 1取最大值对应的阈值 youden tpr - fpr best_idx np.argmax(youden) best_thr thresholds[best_idx] tn, fp, fn, tp confusion_matrix( y_true, (y_pred_proba best_thr).astype(int) ).ravel() sensitivity tp / max(tp fn, 1) specificity tn / max(tn fp, 1) print(fAUC{auc:.4f}, 最优阈值{best_thr:.3f}, f敏感度{sensitivity:.3f}, 特异度{specificity:.3f}) return auc, best_thr注意这个评估脚本也用于第4章的深度学习模型只不过输入换成网络输出的概率。在类别不平衡场景里AUC比准确率高一个数量级的信息量因为它不依赖阈值。还有一个细节约登指数假设敏感度和特异度同等重要肺结节筛查场景可能更在意敏感度这时候就不该选约登指数而是直接取敏感度大于0.95里的最高特异度对应的阈值。代码里把阈值选择逻辑抽出来就是为了方便按业务诉求改口径。4. 深度学习路线2D与3D卷积分类的代码骨架和参数敏感点当样本量超过一两千且算力充足时深度学习通常能超过手工特征加XGBoost。这条路线又分2D和3D2D方案用结节切面图配合预训练图像分类模型代码简单、迁移学习效果好3D方案保留完整的空间结构但参数量大、容易过拟合。这一章不讨论哪个绝对更好而是把两条路的代码骨架都摆出来参数怎么调、输入张量怎么组直接看代码就能复现。先提个醒医学影像任务里模型的选择空间没有预处理和损失函数的调参空间大别在backbone上反复横跳把2D预训练模型跑透再考虑3D。4.1 2D方案用MobileNetV2在单张轴向切片上做二分类2D方案里最常见的做法是对每个结节取横断位中间切片或者取三视图横断位、冠状位、矢状位拼接成三通道。取三视图时要注意三个方向的切片必须从同一个重采样后的ROI体积中切保持窗位窗宽一致否则通道间的数值分布不同预训练模型的BN层统计量会被打乱。下面用MobileNetV2做backbone原因是它在ImageNet上的预训练权重可复用、参数少、在小数据集上不容易过拟合。我一般不用ResNet50起步参数量差距在医学图像小数据集上很现实。import torch import torch.nn as nn import torchvision.models as models class NoduleClassifier2D(nn.Module): 2D肺结节分类器MobileNetV2骨干 自定义分类头。 def __init__(self, num_classes2, dropout0.3): super().__init__() # 使用ImageNet预训练权重迁移学习中冻结前几层是常见做法 self.backbone models.mobilenet_v2(weightsmodels.MobileNet_V2_Weights.DEFAULT) in_features self.backbone.classifier[1].in_features self.backbone.classifier nn.Sequential( nn.Dropout(pdropout), nn.Linear(in_features, num_classes), ) def forward(self, x): # x: (B, 3, H, W)三视图输入时H,W一般为224或256 return self.backbone(x)def build_2d_input(roi_volume): 从48mm ROI体积中取横断位、冠状位、矢状位三个中间切片堆成3通道。 z_mid roi_volume.shape[0] // 2 y_mid roi_volume.shape[1] // 2 x_mid roi_volume.shape[2] // 2 slice_axial roi_volume[z_mid, :, :] # 横断位 slice_coronal roi_volume[:, y_mid, :] # 冠状位 slice_sagittal roi_volume[:, :, x_mid] # 矢状位 # 归一化已经在预处理阶段做过这里直接堆叠 img np.stack([slice_axial, slice_coronal, slice_sagittal], axis0) return torch.from_numpy(img).float().unsqueeze(0)这两个代码块的配合逻辑是先由预处理管线得到48x48x48的ROI体积build_2d_input提取三个正交方向的中间切片并堆成(1,3,48,48)张量再缩放到224后进模型。值得注意的是MobileNetV2的输入规格是224x224中间的尺寸缩放建议用双线性插值在DataLoader里做而不是离线存224的图。离线缩放的后悔药少因为一旦存了固定尺寸之后想换尺度就得重新生成数据。在线缩放只多几行transform代码却能保留原始分辨率信息。4.2 3D方案把完整结节体积喂给3D卷积网络3D方案放弃预训练模型的便利从零训练轻量3D卷积网络。这里我展示的是一个四层3D卷积加BatchNorm的紧凑结构参数量比标准3D ResNet小得多适合在几百个样本的场景下训练。3D网络输入张量的形状是(B, 1, D, H, W)D是深度方向与第2章里volume的轴顺序一致。import torch import torch.nn as nn class NoduleClassifier3D(nn.Module): 轻量3D卷积分类器适合几百到一两千样本的小规模肺结节数据。 def __init__(self, num_classes2): super().__init__() self.features nn.Sequential( nn.Conv3d(1, 32, kernel_size3, padding1), nn.BatchNorm3d(32), nn.ReLU(inplaceTrue), nn.MaxPool3d(2), nn.Conv3d(32, 64, kernel_size3, padding1), nn.BatchNorm3d(64), nn.ReLU(inplaceTrue), nn.MaxPool3d(2), nn.Conv3d(64, 128, kernel_size3, padding1), nn.BatchNorm3d(128), nn.ReLU(inplaceTrue), nn.MaxPool3d(2), ) self.classifier nn.Sequential( nn.AdaptiveAvgPool3d((1, 1, 1)), nn.Flatten(), nn.Linear(128, num_classes), ) def forward(self, x): return self.classifier(self.features(x))def build_3d_input(roi_volume): 把48mm ROI体积整理成模型需要的(B,1,D,H,W)形状。 vol roi_volume[np.newaxis, np.newaxis, ...] # (1, 1, D, H, W) return torch.from_numpy(vol).float()3D卷积有个参数敏感点要专门说感受野与ROI尺寸。48mm体积经过三次MaxPool3d(2)特征图缩到6x6x6每个体素对应的感受野刚好覆盖48mm范围所以小结节的信息没有被压掉。如果你的ROI只裁了32mm三层下采样后是4x4x4感受野偏细碎分类性能可能反而不如2D方案。另一个常见误区是直接对ROI做零填充到64x64x64来凑batch这是给网络注水不如先把ROI裁准。样本量不够时3D网络很容易过拟合我的建议是配合强数据增广随机旋转、随机缩放、随机裁剪这在4.3节里一并说。4.3 训练流程的通用设置学习率、早停、评价指标怎么配深度模型的训练代码套路相对固定但参数选择决定了结果上限。学习率方面2D迁移学习用1e-4是安全起点3D从头训练建议降到5e-5优化器用AdamW配合weight_decay1e-4学习率调度用余弦退火目标是让后期loss震荡小一点。损失函数在类别不平衡时用带权重的CrossEntropyLoss权重直接用label_rad 1 - (label_num / total_num)这类逆频率公式。下面给出2D训练循环的完整代码3D只需把DataLoader的返回张量换成4D即可。from torch.utils.data import DataLoader, TensorDataset from sklearn.metrics import roc_auc_score def train_2d_model(train_images, train_labels, val_images, val_labels, epochs60): 训练2D肺结节分类网络每个epoch结束在验证集上计算AUC并保存最优权重。 device torch.device(cuda if torch.cuda.is_available() else cpu) train_ds TensorDataset(torch.from_numpy(train_images), torch.from_numpy(train_labels)) val_ds TensorDataset(torch.from_numpy(val_images), torch.from_numpy(val_labels)) train_loader DataLoader(train_ds, batch_size32, shuffleTrue, num_workers4) val_loader DataLoader(val_ds, batch_size64, shuffleFalse) model NoduleClassifier2D(num_classes2).to(device) # 按训练集标签统计逆频率作为类别权重降低多数类对loss的主导 neg int((train_labels 0).sum()) pos int((train_labels 1).sum()) weights torch.tensor([pos / (neg pos), neg / (neg pos)], dtypetorch.float32).to(device) criterion nn.CrossEntropyLoss(weightweights) optimizer torch.optim.AdamW(model.parameters(), lr1e-4, weight_decay1e-4) scheduler torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_maxepochs) best_auc 0.0 for epoch in range(epochs): model.train() for inputs, labels in train_loader: inputs, labels inputs.to(device), labels.to(device) logits model(inputs) loss criterion(logits, labels) optimizer.zero_grad() loss.backward() optimizer.step() scheduler.step() model.eval() val_proba [] with torch.no_grad(): for inputs, _ in val_loader: val_proba.append(torch.softmax(model(inputs.to(device)), dim1)[:, 1].cpu().numpy()) val_auc roc_auc_score(val_labels, np.concatenate(val_proba)) print(fEpoch {epoch1}/{epochs}, val AUC: {val_auc:.4f}) if val_auc best_auc: best_auc val_auc torch.save(model.state_dict(), best_nodule_2d.pt) return best_auc这段示例代码里有几处直接影响结果的地方需要展开。weights的赋值逻辑很多人会写反类别权重应该给少数类更大的权重所以恶性样本pos占比低时它的权重设置为neg/(negpos)来放大惩罚代码里正是这个写法。CosineAnnealingLR配合早停实际使用时我见过最多的问题是学习率还没降到低位验证集就已经过拟合停了所以最好把epochs设得大一些80到100让余弦退火完整走完一个周期。还有一点是best_auc只保存模型状态字典验证时如果要复现最好连同预处理参数一起记录否则不同版本的归一化参数会让模型表现漂移。3D训练只需要把模型类换成NoduleClassifier3Dbatch_size建议减半显存不够时不要强行加batch改用梯度累积。5. 肺结节分类常见踩坑与排查样本泄漏、DICOM窗口、过拟合代码跑通只是第一步结果可靠才是关键。这一章集中写我在这个方向上踩过的几个真实的坑每一条都是现象、原因、解决三步。如果你发现验证集AUC异常高但实际临床数据表现差或者训练loss不断下降但验证指标震荡大概率能在下面几条里找到对应的问题。5.1 训练验证样本泄漏同一病人的多个切片同时出现在两边现象验证集AUC高达0.98但换了新数据或者换到医院的另一个批次AUC掉到0.8以下。原因同一个病人的CT扫描被重采样后切成多个结节ROI或者同一个结节在增广时生成了多个样本它们被随机分到了训练集和验证集。模型其实在认病人特征比如钙化形态、扫描设备伪影而不是在认结节良恶性。解决按病人ID做分组划分保证同一个病人的所有ROI只出现在同一个fold里。# 错误的随机划分同一个病人可能同时出现在两边 # train_idx, val_idx train_test_split(np.arange(n_samples), test_size0.2) # 正确的做法先按病人ID分组再对组做划分 def split_by_patient(patient_ids, labels, test_ratio0.2, seed42): 病人级划分同一病人的所有结点只进训练集或验证集。 unique_patients np.unique(patient_ids) rng np.random.RandomState(seed) test_patients set(rng.choice(unique_patients, sizeint(len(unique_patients) * test_ratio), replaceFalse)) train_idx [i for i, pid in enumerate(patient_ids) if pid not in test_patients] val_idx [i for i, pid in enumerate(patient_ids) if pid in test_patients] return train_idx, val_idx这段代码的判断逻辑就是按patient_ids过滤索引而不是按样本索引随机切。另一个容易忽略的泄漏点来自数据增强如果增广发生在划分之前同一个母样本生成的两个副本可能分到两侧一定要先划分再增广。5.2 DICOM读取少做一步全图发黑或发白十有八九是rescale没处理现象把ROI可视化出来肺窗下背景不是黑色结节区域亮成一团模型训练loss不下降。原因原始CT存储值是uint16很多扫描设备用不同的RescaleSlope和RescaleIntercept漏掉这一步等于把原始存储值直接当HU用。不同设备的存储值范围差异很大有的设备斜率是2有的是1截距也各不相同。解决读取时严格应用第2章的换算公式并在预处理后做一次可视化抽查。def sanity_check_hu_range(volume): 打印体数据的HU范围快速确认rescale是否正确。 print(fHU min{volume.min():.1f}, max{volume.max():.1f}) # 人体软组织的HU一般在-200~100之间肺实质在-900~-500之间 # 如果看到min为负几千、max为几万说明存储值没有正确换算我在代码里加这个函数是因为肉眼检查DICOM显示器上的窗位是一回事检查数据数组是另一回事。实际项目里有个血泪经验某个批次的数据扫描协议改了新数据RescaleSlope是2旧代码写死了只乘1跑了一个星期才发现。所以除了读字段最好在预处理链路的入口打印一次HU范围保持对数据分布变动的敏感度。5.3 小样本过拟合先跑特征工程基线再上深度学习现象深度学习模型训练集AUC到0.99验证集只有0.85且每次随机种子不一样验证指标浮动超过5个百分点。原因样本量只有几百模型容量过大网络把训练集的噪声直接背下来了。解决先用第3章的特征工程加XGBoost跑基线如果基线AUC已经在0.85以上深度学习至少要做到0.90才有替换价值。如果非要上深度学习2D预训练模型配合Migration Learning比3D从头训练更抗过拟合。另一个被低估的干预手段是数据增强。肺结节分类里我常用的增广包括随机旋转10度、随机缩放0.9到1.1倍、随机平移几个像素、轴向翻转。注意不要做颜色抖动或高斯噪声那是自然图像增广的用法CT值本身有物理意义随意改变数值分布会破坏窗位窗宽的信息。增广强度也要控制旋转角度太大可能让结节边缘形态失真10度已经够用。5.4 重采样插值选错mask用线性插值会磨掉边界现象结节分割的mask在重采样后边缘糊成一片ROI裁剪位置偏移分类模型训练正常但推理时结果波动大。原因对二值mask用了order3或order1插值边界体素被平均成中间值mask不再是二值。解决mask重采样强制order0最近邻。而且mask的坐标轴排序必须和体数据一致否则重采样结果会出现空间错位。这也是为什么我建议在预处理阶段把mask和CT放在同一个函数里处理用同一套factors。CT用三阶插值保留灰度细节mask用零阶保留边界硬属性两者分开调用zoom但使用相同的重采样倍率。def resample_ct_and_mask(ct, mask, orig_spacing, target_spacing(1.0, 1.0, 1.0)): CT用三阶样条插值mask用最近邻插值共享同一倍率。 factors np.array(orig_spacing, dtypenp.float32) / np.array(target_spacing, dtypenp.float32) ct_resampled zoom(ct, factors, order3) mask_resampled zoom(mask, factors, order0) # 二值mask绝不能用order1 return ct_resampled, mask_resampled这个坑的隐蔽之处在于mask出错往往不会让程序报错只是分类指标悄悄变差。如果发现同一套CT数据在不同软件里读出来的结节直径不一致先检查重采样插值方式这个问题在公开数据集上也有因为有的预处理代码对mask和CT统一用线性插值导致标注边界系统性膨胀。排查方法很简单重采样后输出mask的取值集合如果只有{0,1}说明没问题如果出现了0.5这种中间值就是插值选错了。6. 验证与进阶技巧先跑通小规模代码再用FROC和可解释性验收最后一章不讲大而全的部署方案分享三个我实际在用的验证和进阶技巧。它们能帮你在模型训练完成之后用更可靠的方式证明它真的可用也能帮你在向同事或导师汇报时用更直观的证据说服对方。6.1 冷启动先用一百个结节样本跑通全流程第一次在高成本数据集上实验时我的习惯是先抽一百个样本把预处理、训练、评估整个代码路径完整跑一遍。这不是浪费时间而是用最小成本暴露代码问题。比如DICOM读取的字段缺失、DataLoader的padding尺寸不一致、显存溢出在全量数据上排查一次的时间成本远高于小样本。一百个样本上即使结果没有统计意义也能确认pipeline是通畅的。我会在这一步固定所有随机种子PyTorch、NumPy、Python内置random保证后续全量实验可以复现。另外把预处理中每个阶段的输出shape打印出来比对预期这比在训练中断后再逐步debug高效得多。6.2 用FROC曲线评估检测加分类的联合性能如果你的输入不只是医生标注好的ROI而是先有一个检测器给出候选框再对候选框做分类那么单纯看分类AUC不够因为候选框本身的定位误差会直接影响分类输入。这时候应该用FROC曲线横轴是平均每例CT的假阳性个数纵轴是检测灵敏度它把检测和分类的级联误差统一在一次评估里。LUNA16公开数据集上的评估标准就是FROC我一般按七个预设假阳性率0.125到8计算平均灵敏度作为模型的中期验收指标。平均假阳性数/例0.1250.250.51248敏感性0.620.680.740.790.840.870.90上面是一组示意数据临床项目中每个假阳性率下的敏感性来自检测加分类整条链路在测试集上的输出不是单独给分类器算的。如果你的分类任务是纯ROI输入用AUC就够了FROC不需要强行引入。但如果你未来要做端到端的结节筛查系统FROC是躲不掉的验收方式。6.3 用Grad-CAM看模型到底在盯结节的哪块区域深度学习模型在医学影像上训练完最常被问的问题是它凭什么做判断。Grad-CAM是最快的可视化手段把最后一层卷积的特征图按梯度加权叠加到原图上能看到模型关注的区域。对肺结节分类来说合理的关注区域应该是结节本体和边缘毛刺部分而不是周围血管或扫描噪声。如果高亮区域飘到背景上说明模型学到的是伪相关第一反应不是调模型而是回查预处理和划分是否有泄漏。下面给一个不依赖三方库的最小Grad-CAM实现。def grad_cam_2d(model, input_tensor, target_class1): 用钩子提取2D卷积网络最后一层特征和梯度生成Grad-CAM热力图。 activations {} gradients {} def save_act(name): def hook(module, inp, out): activations[name] out return hook def save_grad(name): def hook(module, inp, out): gradients[name] out[0] return hook # MobileNetV2的最后一层卷积通常在features.18 layer_name backbone.features.18 target_layer dict(model.named_modules())[layer_name] target_layer.register_forward_hook(save_act(layer_name)) target_layer.register_full_backward_hook(save_grad(layer_name)) model.eval() model.zero_grad() logits model(input_tensor) logits[0, target_class].backward() act activations[layer_name][0] # (C, H, W) grad gradients[layer_name][0] weights grad.mean(dim(1, 2), keepdimTrue) # 全局平均池化得到每通道权重 cam torch.relu((weights * act).sum(dim0)) # 加权求和再ReLU cam cam / (cam.max() 1e-6) # 归一化到[0,1] return cam.cpu().detach().numpy()这个实现的原理是用前向钩子拿激活、用后向钩子拿梯度然后做全局平均池化得到通道权重本质上就是Grad-CAM原文的做法。backbone.features.18这个层名只适用于MobileNetV2换backbone时要换成对应最后一层卷积的路径。注意register_full_backward_hook拿到的gradient和PyTorch老版本里register_backward_hook的行为有差异新版请用前者。可视化时把热力图用cv2.applyColorMap转成伪彩色以0.4的透明度叠加到原切片上观察热点分布是否吻合结节边缘和内部实性部分。我个人的习惯是每个模型迭代都保留一张Grad-CAM截图和OOF预测概率一起存档。时间久了回头看这些截图能帮你快速判断某个实验改动到底是让模型更关注病灶还是让它学会了别的捷径。医学影像AI这条路数据和验证口径的重要性远大于模型结构上的炫技先让代码可复现、指标可解释比追一两分AUC更值得投入。希望这篇代码笔记能帮你顺利跑通自己的肺结节分类实验。本文还有配套的精品资源点击获取
简介面向医学影像分析与深度学习初学者提供肺结节检测与分类的完整Python实现以Jupyter Notebook交互形式呈现覆盖从CT影像读入到结节分类的算法流程。压缩包共439个文件大小约195.91MB主要包含397个DICOM原始影像、13个Notebook代码文件、CSV标注与元数据表、模型权重及可视化图表功能模块涵盖数据划分、预处理、模型搭建、训练评估与预测验证。已有2415人学习下载。项目经过两周调试前三个代码文件可顺利运行首个文件注释详尽能帮助阅读者快速理解卷积神经网络在医学影像分类中的应用。整体代码结构清晰从LIDC-IDRI数据集处理到结果可视化均有对应实现适合作为算法复现、实验扩展或毕业设计的参考素材也可迁移至其他医学影像分类任务。1. 肺结节分类算法代码先想清楚你要做的是哪一版分类再动手写代码肺结节分类算法代码这个标题听着像是一个可以直接拿来跑的东西但实际动手时会发现找错任务边界才是大多数项目翻车的起点。临床上最常见的肺结节分类是良恶性二分类输入可以是医生标注好的结节区域也可以是从检测算法输出的候选框里截出来的ROI。前者是纯分类问题后者是检测加分类的级联问题两者在数据组织、标签含义和评估方式上差别很大。这篇文章只做一件事把从DICOM原始数据到得到可评估分类结果的每一段代码怎么组织、参数怎么设、坑在哪讲清楚覆盖特征工程加XGBoost和2D/3D深度学习两条主流路线适合刚接触医学影像的算法工程师和正在复现论文的学生。先别急着选模型先把数据管线和评估口径定下来后面才有得聊。2. 预处理是分类精度的第一道闸门DICOM读取、重采样与窗位窗宽任何肺结节分类算法的代码都是从预处理开始的。CT原始DICOM序列并不是一个可以直接丢给神经网络的张量它的像素值是存储值不是物理意义上的CT值层间间距和面内间距不一致不同扫描协议的窗位窗宽也不同。这三项不处理后面喂给模型的数据就自带偏差分类精度上不去时你甚至会怀疑是模型选型的问题。常见做法是先做四道工序读取DICOM并换算HU、重采样到各向同性、按标注中心裁剪ROI、加窗加归一化。每一步都有参数讲究不能只调库不思考。2.1 DICOM体数据读取先把存储值换算成真实的CT值HU读取DICOM的常用库是pydicom和SimpleITK。我一般用pydicom手动读因为需要同时处理RescaleSlope和RescaleIntercept这两个字段手动读不容易漏。CT图像在PACS里存的通常是整数存储值真实CT值亨氏单位HU必须通过存储值 * RescaleSlope RescaleIntercept换算这一步漏了肺窗和纵隔窗显示都不对更别提给模型训练了。import os import glob import numpy as np import pydicom def load_ct_volume(dicom_dir): 读取一个序列的DICOM切片按空间位置排序成3D体数据返回HU数组和间距。 slices [] for dcm_path in sorted(glob.glob(os.path.join(dicom_dir, *.dcm))): ds pydicom.dcmread(dcm_path) # 关键存储值转成真实CT值(HU)不乘RescaleSlope会导致HU整体偏小 hu ds.pixel_array.astype(np.float32) * float(ds.RescaleSlope) float(ds.RescaleIntercept) pos_z float(ds.ImagePositionPatient[2]) slices.append((pos_z, hu)) slices.sort(keylambda item: item[0]) volume np.stack([item[1] for item in slices], axis0) # (Z, H, W) spacing np.array([ float(ds.PixelSpacing[0]), # 行间距单位mm float(ds.PixelSpacing[1]), # 列间距单位mm float(slices[-1][0] - slices[0][0]) / max(len(slices) - 1, 1), # 层间距 ], dtypenp.float32) return volume, spacing这段代码有三个地方需要说明。第一排序用的ImagePositionPatient[2]是切片在病人坐标系里的Z坐标比直接按文件名排序可靠因为扫描方向可能从头到脚也可能反过来。第二pixel_array是uint16原始存储值先转float32再计算避免溢出。第三返回的spacing顺序和volume轴的顺序对应volume是(Z,H,W)spacing对应(层间距, 行间距, 列间距)后面重采样和裁剪都依赖这个约定。如果个别DICOM没有RescaleSlope字段代码会抛异常处理方法是在循环里用getattr(ds, RescaleSlope, 1)兜底。2.2 重采样与ROI裁剪把不同扫描参数的结节统一到同一个坐标网格不同医院的CT扫描参数差异很大层间距有的0.5mm有的2mm面内分辨率也不一样。如果不做重采样同一物理尺寸的结节在不同样本里占的体素数不同卷积核的感受野覆盖范围也跟着变。常见做法是用scipy.ndimage的zoom函数把体数据线性重采样到各向同性1mm。插值阶数方面HU体数据用三阶样条order3但后续如果涉及mask标注一定要用最近邻否则边界会被磨掉。import numpy as np from scipy.ndimage import zoom def resample_volume(volume, orig_spacing, target_spacing(1.0, 1.0, 1.0)): 把非各向同性CT体数据重采样到目标间距默认1x1x1 mm。 factors np.array(orig_spacing, dtypenp.float32) / np.array(target_spacing, dtypenp.float32) return zoom(volume, factors, order3)重采样后再裁剪就省心很多因为三个方向的间距都是1mm不用再分别算半宽。以标注中心为原点裁固定物理尺寸比如48mm那么体素尺寸就是48x48x48def crop_roi(volume, center_voxel, size_mm48): 按标注中心裁剪固定物理尺寸的ROI适用于重采样到1mm后的数据。 half size_mm // 2 z, y, x [int(round(c)) for c in center_voxel] z0, z1 max(z - half, 0), min(z half, volume.shape[0]) y0, y1 max(y - half, 0), min(y half, volume.shape[1]) x0, x1 max(x - half, 0), min(x half, volume.shape[2]) roi volume[z0:z1, y0:y1, x0:x1] # 边缘不足时用最小像素值填充避免直接裁回小尺寸 if roi.shape ! (size_mm, size_mm, size_mm): pad [(half, half), (half, half), (half, half)] roi np.pad(roi, pad, modeconstant, constant_values-1000)[ z-half:zhalf, y-half:yhalf, x-half:xhalf ] return roi这段话里的重点在最后那个填充分支。靠近肺边缘或扫描边缘的结节裁剪窗口可能超出体数据范围直接切片会让样本之间尺寸不一致DataLoader没法打包。用-1000填充是医学影像里的惯例因为空气的HU值就是-1000。注意center_voxel是在重采样后的坐标空间如果是原始坐标需要先按重采样倍率换算或者干脆先重采样完再用代码里的坐标变换函数统一处理这也是为什么我建议把重采样放在ROI裁剪之前。2.3 窗位窗宽与归一化你看到的图和你喂给网络的数不是一回事CT值是16bit范围的整数直接归一化到[0,1]会损失软组织对比度。临床上医生看肺部结节用的是肺窗窗位约-600HU窗宽约1500HU也就是只保留-1350到150HU这个区间。算法代码里同样应该先加窗再归一化这比全局min-max归一化更符合读片逻辑。窗位窗宽的选择本身就是调参点找磨玻璃结节GGO和找实性结节合适的窗不完全一样。def apply_window(volume, center-600, width1500): 肺窗窗位窗宽调整常见窗口肺窗(-600,1500)纵隔窗(40,400)。 lower center - width / 2.0 upper center width / 2.0 return np.clip(volume, lower, upper) def normalize_to_unit_range(volume): 把加窗后的数据线性归一化到[0,1]避免模型输入数值范围不稳定。 vmin, vmax float(volume.min()), float(volume.max()) return (volume - vmin) / (vmax - vmin 1e-6) # 完整预处理链路的组合调用 def preprocess_pipeline(dicom_dir, center_voxel): vol, spacing load_ct_volume(dicom_dir) vol resample_volume(vol, spacing) # 1mm各向同性 roi crop_roi(vol, center_voxel, size_mm48) # 48mm立方体 roi apply_window(roi, center-600, width1500) # 肺窗 roi normalize_to_unit_range(roi) # [0,1] return roi.astype(np.float32)我在数据集实验里对比过两组设置一组只做min-max归一化另一组先加肺窗再归一化在同一个XGBoost特征管线里AUC大约有3到5个百分点的差距。这里想强调的不是绝对值而是加窗这个操作相当于给模型提供了先验知识告诉它哪些灰度范围是有效的。对于要做2D切片的项目一般从ROI体积里取中间切片或最大面积切片我个人习惯取中间切片配三视图横断位、冠状位、矢状位做三通道输入这个在第四章会展开。3. 特征工程加XGBoost小样本场景下最稳的肺结节分类方案很多入门项目一上来就堆3D卷积网络但肺结节标注样本通常只有几百到一两千个深度学习在这个量级很容易过拟合。特征工程加XGBoost是数据量不足时更稳的一条路也是我建议先做的基线方案。基线跑通后深度学习的增量才有对比意义。这一章的路线是从预处理后的ROI里提取形态和纹理特征组织成表格数据再用XGBoost做二分类最后用AUC和敏感度评估。整段代码可以直接跑通但需要说明的是特征提取的参数直接决定分类上限。3.1 手工特征提取形态、直方图、GLCM和LBP怎么落到代码里肺结节的特征可以分成三类形态特征直径、体积、球度、边缘毛刺、灰度统计特征均值、方差、偏度、峰度、纹理特征GLCM、LBP、HOG。医学影像社区里有个pyradiomics库可以直接提上千个特征但这里我用更透明的手工提取方式方便排查特征含义和复现。注意GLCM和LBP都是基于2D切片的所以先取ROI的中间切片再转成uint8。from skimage.feature import local_binary_pattern, graycomatrix, graycoprops def extract_texture_features(roi_2d): 从2D结节切片提取LBP和GLCM纹理特征返回dict。 feats {} # 转uint8GLCM和LBP都要求输入是整数灰度图 img8 np.clip(roi_2d * 255.0, 0, 255).astype(np.uint8) # LBP局部二值模式uniform模式P8,R1输出9维直方图 lbp local_binary_pattern(img8, P8, R1, methoduniform) lbp_hist, _ np.histogram(lbp, binslist(range(10)), densityTrue) for i, v in enumerate(lbp_hist): feats[flbp_bin_{i}] float(v) # GLCM灰度共生矩阵统计距离1和3角度0/45/90/135度 glcm graycomatrix( img8, distances[1, 3], angles[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4], levels256, symmetricTrue, normedTrue, ) feats[glcm_contrast] float(graycoprops(glcm, contrast).mean()) feats[glcm_homogeneity] float(graycoprops(glcm, homogeneity).mean()) feats[glcm_energy] float(graycoprops(glcm, energy).mean()) feats[glcm_correlation] float(graycoprops(glcm, correlation).mean()) return feats def extract_shape_features(roi_mask, voxel_size_mm1.0): 从二值mask提取形态特征假设mask来自结节分割标注。 feats {} n_voxels roi_mask.sum() feats[volume_mm3] float(n_voxels * voxel_size_mm ** 3) feats[diameter_max_mm] float(np.max(roi_mask.sum(axis2).sum(axis1))) * voxel_size_mm feats[sphericity] float(3.0 * (n_voxels * 4.19) ** (2.0 / 3.0) / (roi_mask.sum(axis2).sum() * 6.0 1e-6)) return feats这里有几个参数值得说明。LBP的P和R分别是邻域点数和半径(8,1)是纹理识别的常用起点想捕获更粗的纹理可以把R调到2但特征维度不变bin数量只受uniform模式影响。GLCM的distances[1,3]表示统计1像素和3像素间距的共生关系前者刻画细纹理后者刻画粗纹理levels256是标准灰度级数。要注意的是这些特征如果直接扔给XGBoost值域差异大体积几百上千LBP密度零点几XGBoost对量纲不敏感但如果是给线性模型或SVM用必须先标准化。3.2 用XGBoost做良恶性分类一版可以直接跑通的最小训练代码特征矩阵准备好后训练部分用XGBoost。这里重点讲三个参数scale_pos_weight处理类别不平衡eval_metricauc让早停基于AUC而非准确率early_stopping_rounds防止过拟合。标签按通常约定恶性为1良性为0。import numpy as np import xgboost as xgb from sklearn.model_selection import StratifiedKFold from sklearn.metrics import roc_auc_score, confusion_matrix def train_xgboost_nodule_classifier(features, labels): 使用5折分层交叉验证训练XGBoost肺结节良恶性分类器。 skf StratifiedKFold(n_splits5, shuffleTrue, random_state42) oof_pred np.zeros(len(labels)) auc_list [] for fold, (tr_idx, va_idx) in enumerate(skf.split(features, labels)): X_tr, X_va features[tr_idx], features[va_idx] y_tr, y_va labels[tr_idx], labels[va_idx] # scale_pos_weight的意义是负样本/正样本数量比让模型关注少数类(恶性) pos_weight float(np.sum(y_tr 0)) / max(np.sum(y_tr 1), 1) model xgb.XGBClassifier( n_estimators200, max_depth4, # 深度4-6足够过深在小样本上容易过拟合 learning_rate0.05, subsample0.8, colsample_bytree0.8, scale_pos_weightpos_weight, eval_metricauc, use_label_encoderFalse, n_jobs8, ) model.fit( X_tr, y_tr, eval_set[(X_va, y_va)], early_stopping_rounds20, verboseFalse, ) oof_pred[va_idx] model.predict_proba(X_va)[:, 1] auc_list.append(roc_auc_score(y_va, oof_pred[va_idx])) print(f5-fold AUC: {np.mean(auc_list):.4f} ± {np.std(auc_list):.4f}) return model, oof_pred这段示例代码的要点在eval_metricauc和early_stopping_rounds20的配合验证集AUC连续20轮不涨就停返回的是验证集上最优的那棵树。scale_pos_weight不是拍脑袋设的它等于负样本数除以正样本数代码里用训练集动态算每个fold略有差异这是刻意为之——避免把测试集信息引入训练。max_depth4对这个量级的特征维度是够的我见过有人把深度调到10训练集AUC到0.99验证集只有0.7典型的过拟合。3.3 类别不平衡下的评估只看准确率会被它骗得很惨肺结节的良恶性分布天然不平衡公开数据集里良性比例通常高于恶性自己从医院拿的数据更夸张。如果直接用准确率评估模型把所有结节都判成良性准确率也可能有70%以上看起来不差实际毫无临床价值。正确的评估口径是AUC、敏感度、特异度和F1。下面的代码在OOF预测结果上同时输出这四个指标以及按约登指数选出的最优阈值。from sklearn.metrics import roc_curve def evaluate_nodule_classifier(y_true, y_pred_proba): 输出AUC、敏感度、特异度以及约登指数对应的最优阈值。 auc roc_auc_score(y_true, y_pred_proba) fpr, tpr, thresholds roc_curve(y_true, y_pred_proba) # 约登指数 敏感度 特异度 - 1取最大值对应的阈值 youden tpr - fpr best_idx np.argmax(youden) best_thr thresholds[best_idx] tn, fp, fn, tp confusion_matrix( y_true, (y_pred_proba best_thr).astype(int) ).ravel() sensitivity tp / max(tp fn, 1) specificity tn / max(tn fp, 1) print(fAUC{auc:.4f}, 最优阈值{best_thr:.3f}, f敏感度{sensitivity:.3f}, 特异度{specificity:.3f}) return auc, best_thr注意这个评估脚本也用于第4章的深度学习模型只不过输入换成网络输出的概率。在类别不平衡场景里AUC比准确率高一个数量级的信息量因为它不依赖阈值。还有一个细节约登指数假设敏感度和特异度同等重要肺结节筛查场景可能更在意敏感度这时候就不该选约登指数而是直接取敏感度大于0.95里的最高特异度对应的阈值。代码里把阈值选择逻辑抽出来就是为了方便按业务诉求改口径。4. 深度学习路线2D与3D卷积分类的代码骨架和参数敏感点当样本量超过一两千且算力充足时深度学习通常能超过手工特征加XGBoost。这条路线又分2D和3D2D方案用结节切面图配合预训练图像分类模型代码简单、迁移学习效果好3D方案保留完整的空间结构但参数量大、容易过拟合。这一章不讨论哪个绝对更好而是把两条路的代码骨架都摆出来参数怎么调、输入张量怎么组直接看代码就能复现。先提个醒医学影像任务里模型的选择空间没有预处理和损失函数的调参空间大别在backbone上反复横跳把2D预训练模型跑透再考虑3D。4.1 2D方案用MobileNetV2在单张轴向切片上做二分类2D方案里最常见的做法是对每个结节取横断位中间切片或者取三视图横断位、冠状位、矢状位拼接成三通道。取三视图时要注意三个方向的切片必须从同一个重采样后的ROI体积中切保持窗位窗宽一致否则通道间的数值分布不同预训练模型的BN层统计量会被打乱。下面用MobileNetV2做backbone原因是它在ImageNet上的预训练权重可复用、参数少、在小数据集上不容易过拟合。我一般不用ResNet50起步参数量差距在医学图像小数据集上很现实。import torch import torch.nn as nn import torchvision.models as models class NoduleClassifier2D(nn.Module): 2D肺结节分类器MobileNetV2骨干 自定义分类头。 def __init__(self, num_classes2, dropout0.3): super().__init__() # 使用ImageNet预训练权重迁移学习中冻结前几层是常见做法 self.backbone models.mobilenet_v2(weightsmodels.MobileNet_V2_Weights.DEFAULT) in_features self.backbone.classifier[1].in_features self.backbone.classifier nn.Sequential( nn.Dropout(pdropout), nn.Linear(in_features, num_classes), ) def forward(self, x): # x: (B, 3, H, W)三视图输入时H,W一般为224或256 return self.backbone(x)def build_2d_input(roi_volume): 从48mm ROI体积中取横断位、冠状位、矢状位三个中间切片堆成3通道。 z_mid roi_volume.shape[0] // 2 y_mid roi_volume.shape[1] // 2 x_mid roi_volume.shape[2] // 2 slice_axial roi_volume[z_mid, :, :] # 横断位 slice_coronal roi_volume[:, y_mid, :] # 冠状位 slice_sagittal roi_volume[:, :, x_mid] # 矢状位 # 归一化已经在预处理阶段做过这里直接堆叠 img np.stack([slice_axial, slice_coronal, slice_sagittal], axis0) return torch.from_numpy(img).float().unsqueeze(0)这两个代码块的配合逻辑是先由预处理管线得到48x48x48的ROI体积build_2d_input提取三个正交方向的中间切片并堆成(1,3,48,48)张量再缩放到224后进模型。值得注意的是MobileNetV2的输入规格是224x224中间的尺寸缩放建议用双线性插值在DataLoader里做而不是离线存224的图。离线缩放的后悔药少因为一旦存了固定尺寸之后想换尺度就得重新生成数据。在线缩放只多几行transform代码却能保留原始分辨率信息。4.2 3D方案把完整结节体积喂给3D卷积网络3D方案放弃预训练模型的便利从零训练轻量3D卷积网络。这里我展示的是一个四层3D卷积加BatchNorm的紧凑结构参数量比标准3D ResNet小得多适合在几百个样本的场景下训练。3D网络输入张量的形状是(B, 1, D, H, W)D是深度方向与第2章里volume的轴顺序一致。import torch import torch.nn as nn class NoduleClassifier3D(nn.Module): 轻量3D卷积分类器适合几百到一两千样本的小规模肺结节数据。 def __init__(self, num_classes2): super().__init__() self.features nn.Sequential( nn.Conv3d(1, 32, kernel_size3, padding1), nn.BatchNorm3d(32), nn.ReLU(inplaceTrue), nn.MaxPool3d(2), nn.Conv3d(32, 64, kernel_size3, padding1), nn.BatchNorm3d(64), nn.ReLU(inplaceTrue), nn.MaxPool3d(2), nn.Conv3d(64, 128, kernel_size3, padding1), nn.BatchNorm3d(128), nn.ReLU(inplaceTrue), nn.MaxPool3d(2), ) self.classifier nn.Sequential( nn.AdaptiveAvgPool3d((1, 1, 1)), nn.Flatten(), nn.Linear(128, num_classes), ) def forward(self, x): return self.classifier(self.features(x))def build_3d_input(roi_volume): 把48mm ROI体积整理成模型需要的(B,1,D,H,W)形状。 vol roi_volume[np.newaxis, np.newaxis, ...] # (1, 1, D, H, W) return torch.from_numpy(vol).float()3D卷积有个参数敏感点要专门说感受野与ROI尺寸。48mm体积经过三次MaxPool3d(2)特征图缩到6x6x6每个体素对应的感受野刚好覆盖48mm范围所以小结节的信息没有被压掉。如果你的ROI只裁了32mm三层下采样后是4x4x4感受野偏细碎分类性能可能反而不如2D方案。另一个常见误区是直接对ROI做零填充到64x64x64来凑batch这是给网络注水不如先把ROI裁准。样本量不够时3D网络很容易过拟合我的建议是配合强数据增广随机旋转、随机缩放、随机裁剪这在4.3节里一并说。4.3 训练流程的通用设置学习率、早停、评价指标怎么配深度模型的训练代码套路相对固定但参数选择决定了结果上限。学习率方面2D迁移学习用1e-4是安全起点3D从头训练建议降到5e-5优化器用AdamW配合weight_decay1e-4学习率调度用余弦退火目标是让后期loss震荡小一点。损失函数在类别不平衡时用带权重的CrossEntropyLoss权重直接用label_rad 1 - (label_num / total_num)这类逆频率公式。下面给出2D训练循环的完整代码3D只需把DataLoader的返回张量换成4D即可。from torch.utils.data import DataLoader, TensorDataset from sklearn.metrics import roc_auc_score def train_2d_model(train_images, train_labels, val_images, val_labels, epochs60): 训练2D肺结节分类网络每个epoch结束在验证集上计算AUC并保存最优权重。 device torch.device(cuda if torch.cuda.is_available() else cpu) train_ds TensorDataset(torch.from_numpy(train_images), torch.from_numpy(train_labels)) val_ds TensorDataset(torch.from_numpy(val_images), torch.from_numpy(val_labels)) train_loader DataLoader(train_ds, batch_size32, shuffleTrue, num_workers4) val_loader DataLoader(val_ds, batch_size64, shuffleFalse) model NoduleClassifier2D(num_classes2).to(device) # 按训练集标签统计逆频率作为类别权重降低多数类对loss的主导 neg int((train_labels 0).sum()) pos int((train_labels 1).sum()) weights torch.tensor([pos / (neg pos), neg / (neg pos)], dtypetorch.float32).to(device) criterion nn.CrossEntropyLoss(weightweights) optimizer torch.optim.AdamW(model.parameters(), lr1e-4, weight_decay1e-4) scheduler torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_maxepochs) best_auc 0.0 for epoch in range(epochs): model.train() for inputs, labels in train_loader: inputs, labels inputs.to(device), labels.to(device) logits model(inputs) loss criterion(logits, labels) optimizer.zero_grad() loss.backward() optimizer.step() scheduler.step() model.eval() val_proba [] with torch.no_grad(): for inputs, _ in val_loader: val_proba.append(torch.softmax(model(inputs.to(device)), dim1)[:, 1].cpu().numpy()) val_auc roc_auc_score(val_labels, np.concatenate(val_proba)) print(fEpoch {epoch1}/{epochs}, val AUC: {val_auc:.4f}) if val_auc best_auc: best_auc val_auc torch.save(model.state_dict(), best_nodule_2d.pt) return best_auc这段示例代码里有几处直接影响结果的地方需要展开。weights的赋值逻辑很多人会写反类别权重应该给少数类更大的权重所以恶性样本pos占比低时它的权重设置为neg/(negpos)来放大惩罚代码里正是这个写法。CosineAnnealingLR配合早停实际使用时我见过最多的问题是学习率还没降到低位验证集就已经过拟合停了所以最好把epochs设得大一些80到100让余弦退火完整走完一个周期。还有一点是best_auc只保存模型状态字典验证时如果要复现最好连同预处理参数一起记录否则不同版本的归一化参数会让模型表现漂移。3D训练只需要把模型类换成NoduleClassifier3Dbatch_size建议减半显存不够时不要强行加batch改用梯度累积。5. 肺结节分类常见踩坑与排查样本泄漏、DICOM窗口、过拟合代码跑通只是第一步结果可靠才是关键。这一章集中写我在这个方向上踩过的几个真实的坑每一条都是现象、原因、解决三步。如果你发现验证集AUC异常高但实际临床数据表现差或者训练loss不断下降但验证指标震荡大概率能在下面几条里找到对应的问题。5.1 训练验证样本泄漏同一病人的多个切片同时出现在两边现象验证集AUC高达0.98但换了新数据或者换到医院的另一个批次AUC掉到0.8以下。原因同一个病人的CT扫描被重采样后切成多个结节ROI或者同一个结节在增广时生成了多个样本它们被随机分到了训练集和验证集。模型其实在认病人特征比如钙化形态、扫描设备伪影而不是在认结节良恶性。解决按病人ID做分组划分保证同一个病人的所有ROI只出现在同一个fold里。# 错误的随机划分同一个病人可能同时出现在两边 # train_idx, val_idx train_test_split(np.arange(n_samples), test_size0.2) # 正确的做法先按病人ID分组再对组做划分 def split_by_patient(patient_ids, labels, test_ratio0.2, seed42): 病人级划分同一病人的所有结点只进训练集或验证集。 unique_patients np.unique(patient_ids) rng np.random.RandomState(seed) test_patients set(rng.choice(unique_patients, sizeint(len(unique_patients) * test_ratio), replaceFalse)) train_idx [i for i, pid in enumerate(patient_ids) if pid not in test_patients] val_idx [i for i, pid in enumerate(patient_ids) if pid in test_patients] return train_idx, val_idx这段代码的判断逻辑就是按patient_ids过滤索引而不是按样本索引随机切。另一个容易忽略的泄漏点来自数据增强如果增广发生在划分之前同一个母样本生成的两个副本可能分到两侧一定要先划分再增广。5.2 DICOM读取少做一步全图发黑或发白十有八九是rescale没处理现象把ROI可视化出来肺窗下背景不是黑色结节区域亮成一团模型训练loss不下降。原因原始CT存储值是uint16很多扫描设备用不同的RescaleSlope和RescaleIntercept漏掉这一步等于把原始存储值直接当HU用。不同设备的存储值范围差异很大有的设备斜率是2有的是1截距也各不相同。解决读取时严格应用第2章的换算公式并在预处理后做一次可视化抽查。def sanity_check_hu_range(volume): 打印体数据的HU范围快速确认rescale是否正确。 print(fHU min{volume.min():.1f}, max{volume.max():.1f}) # 人体软组织的HU一般在-200~100之间肺实质在-900~-500之间 # 如果看到min为负几千、max为几万说明存储值没有正确换算我在代码里加这个函数是因为肉眼检查DICOM显示器上的窗位是一回事检查数据数组是另一回事。实际项目里有个血泪经验某个批次的数据扫描协议改了新数据RescaleSlope是2旧代码写死了只乘1跑了一个星期才发现。所以除了读字段最好在预处理链路的入口打印一次HU范围保持对数据分布变动的敏感度。5.3 小样本过拟合先跑特征工程基线再上深度学习现象深度学习模型训练集AUC到0.99验证集只有0.85且每次随机种子不一样验证指标浮动超过5个百分点。原因样本量只有几百模型容量过大网络把训练集的噪声直接背下来了。解决先用第3章的特征工程加XGBoost跑基线如果基线AUC已经在0.85以上深度学习至少要做到0.90才有替换价值。如果非要上深度学习2D预训练模型配合Migration Learning比3D从头训练更抗过拟合。另一个被低估的干预手段是数据增强。肺结节分类里我常用的增广包括随机旋转10度、随机缩放0.9到1.1倍、随机平移几个像素、轴向翻转。注意不要做颜色抖动或高斯噪声那是自然图像增广的用法CT值本身有物理意义随意改变数值分布会破坏窗位窗宽的信息。增广强度也要控制旋转角度太大可能让结节边缘形态失真10度已经够用。5.4 重采样插值选错mask用线性插值会磨掉边界现象结节分割的mask在重采样后边缘糊成一片ROI裁剪位置偏移分类模型训练正常但推理时结果波动大。原因对二值mask用了order3或order1插值边界体素被平均成中间值mask不再是二值。解决mask重采样强制order0最近邻。而且mask的坐标轴排序必须和体数据一致否则重采样结果会出现空间错位。这也是为什么我建议在预处理阶段把mask和CT放在同一个函数里处理用同一套factors。CT用三阶插值保留灰度细节mask用零阶保留边界硬属性两者分开调用zoom但使用相同的重采样倍率。def resample_ct_and_mask(ct, mask, orig_spacing, target_spacing(1.0, 1.0, 1.0)): CT用三阶样条插值mask用最近邻插值共享同一倍率。 factors np.array(orig_spacing, dtypenp.float32) / np.array(target_spacing, dtypenp.float32) ct_resampled zoom(ct, factors, order3) mask_resampled zoom(mask, factors, order0) # 二值mask绝不能用order1 return ct_resampled, mask_resampled这个坑的隐蔽之处在于mask出错往往不会让程序报错只是分类指标悄悄变差。如果发现同一套CT数据在不同软件里读出来的结节直径不一致先检查重采样插值方式这个问题在公开数据集上也有因为有的预处理代码对mask和CT统一用线性插值导致标注边界系统性膨胀。排查方法很简单重采样后输出mask的取值集合如果只有{0,1}说明没问题如果出现了0.5这种中间值就是插值选错了。6. 验证与进阶技巧先跑通小规模代码再用FROC和可解释性验收最后一章不讲大而全的部署方案分享三个我实际在用的验证和进阶技巧。它们能帮你在模型训练完成之后用更可靠的方式证明它真的可用也能帮你在向同事或导师汇报时用更直观的证据说服对方。6.1 冷启动先用一百个结节样本跑通全流程第一次在高成本数据集上实验时我的习惯是先抽一百个样本把预处理、训练、评估整个代码路径完整跑一遍。这不是浪费时间而是用最小成本暴露代码问题。比如DICOM读取的字段缺失、DataLoader的padding尺寸不一致、显存溢出在全量数据上排查一次的时间成本远高于小样本。一百个样本上即使结果没有统计意义也能确认pipeline是通畅的。我会在这一步固定所有随机种子PyTorch、NumPy、Python内置random保证后续全量实验可以复现。另外把预处理中每个阶段的输出shape打印出来比对预期这比在训练中断后再逐步debug高效得多。6.2 用FROC曲线评估检测加分类的联合性能如果你的输入不只是医生标注好的ROI而是先有一个检测器给出候选框再对候选框做分类那么单纯看分类AUC不够因为候选框本身的定位误差会直接影响分类输入。这时候应该用FROC曲线横轴是平均每例CT的假阳性个数纵轴是检测灵敏度它把检测和分类的级联误差统一在一次评估里。LUNA16公开数据集上的评估标准就是FROC我一般按七个预设假阳性率0.125到8计算平均灵敏度作为模型的中期验收指标。平均假阳性数/例0.1250.250.51248敏感性0.620.680.740.790.840.870.90上面是一组示意数据临床项目中每个假阳性率下的敏感性来自检测加分类整条链路在测试集上的输出不是单独给分类器算的。如果你的分类任务是纯ROI输入用AUC就够了FROC不需要强行引入。但如果你未来要做端到端的结节筛查系统FROC是躲不掉的验收方式。6.3 用Grad-CAM看模型到底在盯结节的哪块区域深度学习模型在医学影像上训练完最常被问的问题是它凭什么做判断。Grad-CAM是最快的可视化手段把最后一层卷积的特征图按梯度加权叠加到原图上能看到模型关注的区域。对肺结节分类来说合理的关注区域应该是结节本体和边缘毛刺部分而不是周围血管或扫描噪声。如果高亮区域飘到背景上说明模型学到的是伪相关第一反应不是调模型而是回查预处理和划分是否有泄漏。下面给一个不依赖三方库的最小Grad-CAM实现。def grad_cam_2d(model, input_tensor, target_class1): 用钩子提取2D卷积网络最后一层特征和梯度生成Grad-CAM热力图。 activations {} gradients {} def save_act(name): def hook(module, inp, out): activations[name] out return hook def save_grad(name): def hook(module, inp, out): gradients[name] out[0] return hook # MobileNetV2的最后一层卷积通常在features.18 layer_name backbone.features.18 target_layer dict(model.named_modules())[layer_name] target_layer.register_forward_hook(save_act(layer_name)) target_layer.register_full_backward_hook(save_grad(layer_name)) model.eval() model.zero_grad() logits model(input_tensor) logits[0, target_class].backward() act activations[layer_name][0] # (C, H, W) grad gradients[layer_name][0] weights grad.mean(dim(1, 2), keepdimTrue) # 全局平均池化得到每通道权重 cam torch.relu((weights * act).sum(dim0)) # 加权求和再ReLU cam cam / (cam.max() 1e-6) # 归一化到[0,1] return cam.cpu().detach().numpy()这个实现的原理是用前向钩子拿激活、用后向钩子拿梯度然后做全局平均池化得到通道权重本质上就是Grad-CAM原文的做法。backbone.features.18这个层名只适用于MobileNetV2换backbone时要换成对应最后一层卷积的路径。注意register_full_backward_hook拿到的gradient和PyTorch老版本里register_backward_hook的行为有差异新版请用前者。可视化时把热力图用cv2.applyColorMap转成伪彩色以0.4的透明度叠加到原切片上观察热点分布是否吻合结节边缘和内部实性部分。我个人的习惯是每个模型迭代都保留一张Grad-CAM截图和OOF预测概率一起存档。时间久了回头看这些截图能帮你快速判断某个实验改动到底是让模型更关注病灶还是让它学会了别的捷径。医学影像AI这条路数据和验证口径的重要性远大于模型结构上的炫技先让代码可复现、指标可解释比追一两分AUC更值得投入。希望这篇代码笔记能帮你顺利跑通自己的肺结节分类实验。本文还有配套的精品资源点击获取
