基于Python的肝脏CT分割与三维重建：从DICOM到3D打印模型全流程

简介本资源为基于Python的肝脏CT图像分割及三维重建完整项目包面向计算机、人工智能、生物医学工程等专业的在校学生与从业人员可用于毕业设计、课程设计、期末大作业及竞赛初期项目立项演示。项目围绕医学影像处理展开涵盖CT图像预处理、肝脏区域分割、三维体绘制与可视化重建等核心环节具有较强的代表性与学习借鉴价值。压缩包共124个文件约105.2MB以42个py源码文件为主体辅以42个png结果图、15个pyc编译文件、10个txt说明文档及5个tar模型包另含xml配置、gif演示与vtk可视化文件结构完整、便于按模块查阅。目前已有180人浏览学习。源码经本地运行与功能测试读者可据此掌握分割算法实现思路、模型调用方式与三维重建流程并可在原有基础上进行二次开发与功能扩展。1. 肝脏CT分割与三维重建从DICOM到可打印模型的完整链路拿到一份腹部增强CT放射科医生要在几百层切片里逐层勾出肝脏边界再凭空间想象力判断肿瘤位置和残肝体积——这个过程在临床上是常态但效率极低。基于Python的肝脏CT图像分割及三维重建要解决的就是把这条链路自动化输入DICOM序列输出肝脏掩膜和可交互的三维网格模型。它适合影像科工程师、医学图像处理方向的学生以及需要做术前规划或科研量化的从业者。整套流程的核心环节只有三个窗宽窗位预处理、分割网络推理、面绘制重建。源码和模型的价值在于把这三个环节的参数和接口固定下来让你不用从零调通。下面按实际落地顺序拆开讲每一步都给可复现的命令和参数。2. 数据预处理DICOM读取、窗宽窗位与HU值归一化2.1 为什么不能直接把DICOM像素丢给模型CT图像的原始像素值叫HUHounsfield Unit范围通常在-1024到3071之间。肝脏实质的HU值大约在40到70肿瘤可能低到20或高到100以上而骨骼能到1000以上。如果直接把原始HU送进网络骨骼和对比剂的极端值会主导梯度肝脏区域的微弱差异被淹没。常见做法是先把HU截断到[-200, 300]这个腹部窗再线性映射到[0, 1]。这个区间能覆盖肝脏、脾脏、肾脏和大部分病灶同时排除骨骼和空气的干扰。另一个坑是不同扫描仪的斜率RescaleSlope和截距RescaleIntercept不同。DICOM里的像素值需要先做HU pixel * slope intercept才能得到真实HU。很多开源代码直接读pixel_array在部分设备上会得到完全错误的结果。import pydicom import numpy as np def load_dicom_series(dicom_dir): 读取DICOM序列并按InstanceNumber排序 slices [] for fname in os.listdir(dicom_dir): ds pydicom.dcmread(os.path.join(dicom_dir, fname)) # 关键用ImagePositionPatient的z坐标排序比InstanceNumber可靠 slices.append(ds) slices.sort(keylambda s: float(s.ImagePositionPatient[2])) return slices def hu_to_normalized(ds, hu_min-200, hu_max300): HU值截断并归一化到[0,1] pixel ds.pixel_array.astype(np.float32) hu pixel * float(ds.RescaleSlope) float(ds.RescaleIntercept) hu np.clip(hu, hu_min, hu_max) return (hu - hu_min) / (hu_max - hu_min)load_dicom_series里用ImagePositionPatient[2]排序而不是InstanceNumber是因为部分设备在多序列拼接时InstanceNumber会重复或跳号。hu_to_normalized的两个参数hu_min和hu_max是腹部窗的典型值如果你做的是肝脏专门分割可以收紧到[-100, 200]让肝脏对比度更高。归一化后的数组直接堆叠成三维体数据形状为(D, H, W)D是层数。2.2 层厚不一致与各向异性重采样CT扫描的层厚常见有1mm、1.25mm、2.5mm、5mm。层厚5mm的序列在z轴方向只有几十层直接做三维重建会得到阶梯状表面。标准做法是把体数据重采样到各向同性间距比如1mm×1mm×1mm。用scipy.ndimage.zoom或SimpleITK都可以但要注意插值顺序图像用三阶样条掩膜用最近邻否则标签会被插值成小数。import SimpleITK as sitk def resample_isotropic(image, mask, target_spacing(1.0, 1.0, 1.0)): 将图像和掩膜重采样到各向同性间距 original_spacing image.GetSpacing() original_size image.GetSize() # 计算重采样后的尺寸 new_size [ int(round(original_size[i] * original_spacing[i] / target_spacing[i])) for i in range(3) ] resampler sitk.ResampleImageFilter() resampler.SetSize(new_size) resampler.SetOutputSpacing(target_spacing) resampler.SetOutputDirection(image.GetDirection()) resampler.SetOutputOrigin(image.GetOrigin()) # 图像用线性插值掩膜用最近邻 resampler.SetInterpolator(sitk.sitkLinear) image_resampled resampler.Execute(image) resampler.SetInterpolator(sitk.sitkNearestNeighbor) mask_resampled resampler.Execute(mask) return image_resampled, mask_resampledtarget_spacing设为(1.0, 1.0, 1.0)是通用选择如果你的GPU显存有限可以放宽到(1.5, 1.5, 1.5)。SetInterpolator对图像和掩膜分别设置是必须的掩膜用线性插值会产生0.3、0.7这样的标签值后续计算Dice时直接报错。重采样后的体数据再切成256×256或512×512的patch送入网络。3. 分割模型选型与推理U-Net、nnU-Net还是MONAI3.1 肝脏分割的模型边界在哪里肝脏分割在医学图像领域已经比较成熟公开数据集LiTS和Sliver07上的Dice能到0.95以上。但实际落地时模型面临的挑战不是肝脏本身而是边界模糊区域肝脏与胃壁、脾脏、心脏相邻处CT值接近梯度信息弱。另外肿瘤浸润区域肝脏边界会变形模型容易把肿瘤漏在肝脏外。选型上U-Net是基线nnU-Net是当前公认的强基线MONAI是工程化封装。如果你要快速跑通用MONAI的UNet加预训练权重最省事如果要刷指标nnU-Net的自适应配置自动选择patch size、归一化方式、损失函数很难被手工调参超过。源码包里如果带的是自定义U-Net重点看它的损失函数Dice BCE组合是标配但肝脏分割里Dice权重通常要高于BCE因为前景占比小。import torch import monai from monai.networks.nets import UNet # MONAI UNet配置4层下采样通道数32起步 model UNet( spatial_dims3, in_channels1, out_channels2, # 背景肝脏 channels(32, 64, 128, 256, 512), strides(2, 2, 2, 2), num_res_units2, normbatch, dropout0.1, ) model model.cuda() model.eval() # 推理滑动窗口高斯权重融合 from monai.inferers import SlidingWindowInferer inferer SlidingWindowInferer( roi_size(128, 128, 128), sw_batch_size4, overlap0.5, modegaussian, ) with torch.no_grad(): logits inferer(input_tensor, model) pred torch.argmax(logits, dim1)channels从32到512是显存和精度的折中如果显存只有8GB把第一层降到16去掉最后一层。num_res_units2表示每个下采样块里有两个残差单元能缓解梯度消失。SlidingWindowInferer的overlap0.5是经验值重叠太少会在拼接处出现接缝太多则推理时间翻倍。modegaussian让窗口边缘权重低、中心权重高融合后边界更平滑。3.2 后处理连通域与形态学修补网络输出的掩膜经常有孤立小区域和内部空洞。肝脏是最大的连通域所以保留最大连通域能去掉大部分假阳性。内部空洞用二值填充补上边缘的毛刺用形态学开运算平滑。from scipy import ndimage def postprocess_liver_mask(mask): 保留最大连通域并填充内部空洞 # 保留最大连通域 labeled, num ndimage.label(mask) if num 1: sizes ndimage.sum(mask, labeled, range(1, num 1)) largest np.argmax(sizes) 1 mask (labeled largest) # 填充内部空洞 mask ndimage.binary_fill_holes(mask) # 开运算平滑边缘 mask ndimage.binary_opening(mask, structurenp.ones((3, 3, 3))) return mask.astype(np.uint8)ndimage.label默认是6连通对肝脏这种块状结构够用。binary_fill_holes会把肝脏内部所有封闭空洞填上包括血管造成的低密度区这是符合预期的——三维重建时内部空洞会导致网格破碎。binary_opening的结构元素用3×3×3再大就会侵蚀肝脏边缘。后处理之后用skimage.measure.marching_cubes提取等值面得到顶点和面片。4. 三维重建marching cubes参数与网格导出4.1 从掩膜到网格的四个关键参数Marching cubes是面绘制的标准算法输入是三维二值体数据输出是三角网格。skimage的实现有三个参数直接影响结果level、spacing和step_size。level设为0.5因为掩膜是0/1二值0.5是等值面位置。spacing要设成体数据的实际物理间距否则重建出来的模型在z轴方向会被拉伸或压缩。step_size控制采样步长设为1表示逐体素计算精度最高但速度慢设为2会跳过部分体素速度快但可能丢失小结构。from skimage import measure import trimesh def mask_to_mesh(mask, spacing(1.0, 1.0, 1.0)): marching cubes提取网格并导出STL verts, faces, normals, values measure.marching_cubes( mask, level0.5, spacingspacing, step_size1, allow_degenerateFalse, ) mesh trimesh.Trimesh(verticesverts, facesfaces, normalsnormals) # 网格简化保留90%的顶点减少面片数 mesh mesh.simplify_quadric_decimation(0.9) # 平滑拉普拉斯平滑迭代10次 mesh mesh.smooth_laplacian(lamb0.5, iterations10) return mesh mesh mask_to_mesh(liver_mask, spacing(1.0, 1.0, 1.0)) mesh.export(liver_model.stl)allow_degenerateFalse会去掉零面积三角形避免后续3D打印切片报错。simplify_quadric_decimation(0.9)把面片数降到原来的90%对视觉影响很小但文件体积明显减小。smooth_laplacian的lamb0.5是平滑强度迭代10次能去掉阶梯感但迭代太多会让肝脏的锐利边缘变圆。导出STL后可以用MeshLab或Blender打开检查重点看有没有自交面和法线翻转。4.2 体积计算与术前规划指标三维重建不只是为了看还要算。肝脏体积、肿瘤体积、残肝体积比Future Liver Remnant, FLR是术前规划的核心指标。体积计算用体素计数乘以单个体素的物理体积比网格体积更准因为网格简化会引入误差。def compute_volumes(liver_mask, tumor_mask, spacing): 计算肝脏、肿瘤和残肝体积毫升 voxel_volume spacing[0] * spacing[1] * spacing[2] / 1000.0 # mm^3转ml liver_vol np.sum(liver_mask) * voxel_volume tumor_vol np.sum(tumor_mask) * voxel_volume # 残肝 肝脏 - 肿瘤假设肿瘤在肝脏内 remnant_vol liver_vol - tumor_vol flr_ratio remnant_vol / liver_vol return { liver_ml: round(liver_vol, 1), tumor_ml: round(tumor_vol, 1), remnant_ml: round(remnant_vol, 1), flr_ratio: round(flr_ratio, 3), }spacing的单位是毫米除以1000转成毫升。FLR比值低于0.25时术后肝衰竭风险显著升高这是临床上的硬指标。如果你的分割掩膜里肿瘤和肝脏有重叠先做tumor_mask tumor_mask liver_mask再计算。体积计算的结果和商业软件如Myrian、IntelliSpace对比误差通常在5%以内主要来源是层厚和分割边界。5. 避坑与排查DICOM方向、显存溢出与网格破碎5.1 现象重建模型上下颠倒或左右镜像原因DICOM的ImageOrientationPatient定义了图像坐标系但不同设备厂商的坐标系约定不同。直接按数组索引重建在部分数据上会得到镜像结果。解决用ImagePositionPatient和ImageOrientationPatient构建仿射矩阵通过SimpleITK的GetDirection()获取方向余弦在重采样时保留方向信息。导出网格前用trimesh的apply_transform做一次坐标变换确保模型在RAS坐标系下。5.2 现象推理时CUDA out of memory原因3D U-Net的显存占用和patch size的三次方成正比。128×128×128的patch在FP32下大约占4GB加上梯度缓存和中间特征图8GB显存很容易爆。解决把patch降到96×96×96或者用混合精度torch.cuda.amp。如果还不行改用滑动窗口推理每次只送一个窗口窗口之间重叠0.5。注意sw_batch_size不要设太大它控制并行窗口数设成2或4就够。5.3 现象STL文件导入3D打印机后切片失败原因marching cubes输出的网格可能有自交面、非流形边或法线不一致。3D打印切片软件对网格水密性要求严格。解决用trimesh的mesh.fill_holes()补洞mesh.fix_normals()统一法线mesh.remove_degenerate_faces()去退化面。如果还有问题用mesh.split()拆成多个连通分量只保留最大的那个。导出前用mesh.is_watertight检查返回True才能保证切片成功。5.4 现象分割结果在肝脏顶部或底部缺失原因CT序列的两端层数少肝脏在z轴方向被截断网络在边界处没有足够上下文。解决在预处理时对z轴做镜像padding两端各补10层推理完再裁掉。或者用nnU-Net的自适应patch它会根据图像尺寸自动调整。如果源码包里没有这个逻辑手动在load_dicom_series之后加一步np.padmode用reflect。5.5 现象Dice系数在验证集上很高但视觉效果差原因Dice对大面积区域敏感肝脏内部的小肿瘤或血管分割错误对Dice影响很小。解决加看HD9595% Hausdorff Distance和ASSDAverage Symmetric Surface Distance这两个指标对边界误差更敏感。另外把肿瘤区域单独算Dice不要混在肝脏整体里。如果源码包只给了Dice自己用medpy.metric.binary.hd95补上。6. 进阶技巧用PyTorch 2.0编译加速与多模态融合6.1 推理速度翻倍torch.compile的实际收益PyTorch 2.0的torch.compile对3D U-Net的推理加速在1.3到1.8倍之间取决于patch size和GPU型号。用法很简单在模型加载后加一行model torch.compile(model, modereduce-overhead)modereduce-overhead适合推理场景它会用CUDA Graph减少kernel启动开销。第一次推理会触发编译耗时几十秒之后每次推理都快。注意torch.compile对动态shape支持有限如果你的patch size每次都变编译会反复触发反而更慢。所以推理时固定patch size用滑动窗口处理不同尺寸的输入。6.2 多模态融合CTMRI的通道拼接如果手头有同一患者的MRI数据可以把T2加权像配准到CT空间作为第二通道送入网络。肝脏在MRI上的边界比CT更清晰尤其是肝硬化患者。配准用SimpleITK的Elastix或ANTs刚性配准加B样条形变。融合时把MRI归一化到和CT相同的[0,1]区间通道维度从1变成2。注意MRI的强度不均匀先做N4偏置场校正否则网络会学到伪影。# 多模态输入CT MRI双通道 ct_tensor torch.from_numpy(ct_normalized).unsqueeze(0).unsqueeze(0) # (1,1,D,H,W) mr_tensor torch.from_numpy(mr_normalized).unsqueeze(0).unsqueeze(0) input_tensor torch.cat([ct_tensor, mr_tensor], dim1) # (1,2,D,H,W)对应的模型in_channels要改成2。如果只有CT不要用零填充假装双通道网络会学到无用的零通道。多模态训练时MRI缺失的样本要单独处理不能直接补零。6.3 一个我踩过的坑模型权重加载的键名不匹配源码包里的.pth文件如果用了DataParallel或DistributedDataParallel训练保存的state_dict键名会带module.前缀。直接model.load_state_dict()会报missing keys。解决state_dict torch.load(model.pth, map_locationcpu) # 去掉module.前缀 new_state_dict {k.replace(module., ): v for k, v in state_dict.items()} model.load_state_dict(new_state_dict, strictTrue)strictTrue会严格检查键名如果有缺失或多余会直接报错比strictFalse安全。如果键名不匹配是因为模型结构改了那就只能逐层对比用model.named_parameters()打印出来和state_dict的键做差集。6.4 验证重建精度的土办法没有体模的情况下用同一患者的复查CT做交叉验证。第一次扫描重建模型第二次扫描再重建两个模型配准后算表面距离。如果平均距离小于2mm说明分割和重建的重复性够用。另一个办法是打印出来用卡尺量但精度受打印机限制。我一般会在导出STL前把网格顶点坐标存成CSV和商业软件的结果做点对点对比这样能定位到具体是哪个区域偏差大。这套流程从DICOM读取到STL导出核心代码不超过300行但参数和边界条件很多。我自己的习惯是每换一个数据集先跑一遍预处理可视化确认HU窗口和重采样没问题再动模型。分割结果出来后不要只看Dice把掩膜叠加到原始CT上逐层翻一遍尤其是肝脏与周围器官的交界处。三维重建的模型导出后用MeshLab检查法线和自交面别等到切片失败才回头找原因。希望帮到你。本文还有配套的精品资源点击获取