Python肝脏CT图像分割与三维重建全流程实战详解

简介：一份基于Python的肝脏CT图像分割与三维重建项目资源，面向医学图像处理、计算机视觉及相关专业的学生和开发者，可应用于毕业设计、课程设计、期末大作业或竞赛立项演示。资源包含完整源码与预训练模型，覆盖从肝脏CT预处理、分割到三维表面重建与可视化的关键流程。压缩包共124个文件，以42个Python脚本和42张效果图像为主，另有模型权重、VTK三维网格文件、XML配置、文本说明及可视化动图等，整体约105.2MB，便于对照源码查看各阶段输出；目前已有180人浏览学习。项目自带说明文档和示例图片，能帮助快速理解分割与重建原理；代码结构清晰、模块划分明确，适合作为图像处理方向的入门进阶范例，也支持在此基础上修改扩展，完成二次开发。

1. 为什么“Python+肝脏CT分割+三维重建”值得你花一个周末跑通

做肝脏CT分割的时候，大家第一反应往往是“有没有现成的模型可以直接用”。实际上下载到手的这套基于Python的肝脏CT图像分割及三维重建源码+模型.zip，最常见的工作流是：用DICOM读取肝脏CT，预处理后交给一个训练好的分割网络，得到肝脏在每一层CT上的掩码，再通过三维重建生成光滑的肝脏3D模型文件。对医工交叉入门者、论文需要对照实验的在校同学，以及想快速验证算法可行性的工程师来说，这是最省事的起点；但前提是你得把数据流、模型输入规范和输出约定都理顺，否则很容易在DICOM转NIfTI这一步就翻车。这篇笔记就是顺着这条管线，把每步的代码、参数和踩过的坑讲透。

2. 数据入口：把CT DICOM序列转成体素矩阵，这步决定后面所有结果

2.1 用pydicom把CT目录读成连续3D数组

DICOM序列文件夹看起来是一堆 .dcm 文件，但能不能正确堆成三维体数据，取决于两件事：切片顺序对不对，rescale参数有没有乘上。很多人直接遍历文件夹读 pixel_array，结果z轴顺序混乱，生成的mask前后颠倒，三维重建出来肝脏是镜像的。

我一般这样写读取逻辑：

import os
import numpy as np
import pydicom
from pydicom.filereader import dcmread

def load_dicom_series(dcm_dir):
    dcm_files = [f for f in os.listdir(dcm_dir) if f.lower().endswith('.dcm')]
    slices = []
    for f in dcm_files:
        ds = dcmread(os.path.join(dcm_dir, f))
        slices.append(ds)

    # 按 ImagePositionPatient 的第3个分量排序，而不是按文件名
    slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))

    # 先验证层厚：从相邻两层的位置差推算
    if len(slices) < 2:
        raise ValueError("切片数量小于2，无法构成三维体数据")
    spacing = (
        float(slices[0].PixelSpacing[0]),
        float(slices[0].PixelSpacing[1]),
        abs(float(slices[1].ImagePositionPatient[2]) - float(slices[0].ImagePositionPatient[2]))
    )

    # 把像素值转成HU（亨氏单位）
    volume = np.zeros((len(slices), slices[0].Rows, slices[0].Columns), dtype=np.int16)
    for i, ds in enumerate(slices):
        slope = float(ds.RescaleSlope) if ds.RescaleSlope is not None else 1.0
        intercept = float(ds.RescaleIntercept) if ds.RescaleIntercept is not None else 0.0
        volume[i] = ds.pixel_array * slope + intercept

    return volume, spacing

这段代码看起来简单，但有两个关键点：一是排序必须用 ImagePositionPatient 的z坐标，很多DICOM文件名是乱序编号，按文件名排序会在z轴产生断层；二是 pixel_array 读出来是原始存储值，必须乘 RescaleSlope 再加 RescaleIntercept 才是真正的HU值。肝脏CT的典型HU范围是 -50 到 250，如果你直接拿原始值当输入，后续归一化就会偏。

2.2 窗宽窗位与重采样：两个决定模型胜负的预热参数

CT图像是12位或16位灰度，窗宽窗位决定了肉眼能看到什么，但对于深度学习模型来说，更重要的是归一化策略。很多开源肝脏分割模型在训练时用 clip = [-100, 200] 或 [-75, 300] 这类窗，然后把像素值线性缩放到 [0, 1] 或 z-score 标准化。

我常用的做法是：

def preprocess_ct(volume_np, window_min=-100, window_max=200):
    clipped = np.clip(volume_np, window_min, window_max)
    # 线性缩放到 0-1
    normalized = (clipped - window_min) / (window_max - window_min)
    return normalized.astype(np.float32)

注意，这个窗不是随便拍的。肝脏实质的CT值通常在40-80HU，肝脏肿瘤在低密度增强形势下可能落在20-50HU，窗底用-100是为了保留肝内血管和胆管的低密度信息，窗顶用200是为了不把钙化点和增强后的肝实质压成纯白。推理时必须和训练时的窗保持一致，否则模型输入分布偏移，分割结果会突然退化。

重采样是第二个隐藏坑。原始DICOM的spacing可能是 (0.7, 0.7, 5)mm，也可能是各向异性 (0.5, 0.5, 1.5)mm。肝脏形态不规则，三维重建的过程对spacing很敏感。常见做法是把各向异性数据线性插值成各向同性体素，例如统一到 1.5mm³：

import SimpleITK as sitk

def resample_to_isotropic(volume_np, spacing, target_spacing=1.5):
    # 转成 SimpleITK 图像再重采样
    sitk_img = sitk.GetImageFromArray(volume_np)
    sitk_img.SetSpacing((spacing[0], spacing[1], spacing[2]))
    original_size = sitk_img.GetSize()
    new_size = [
        int(round(original_size[0] * spacing[0] / target_spacing)),
        int(round(original_size[1] * spacing[1] / target_spacing)),
        int(round(original_size[2] * spacing[2] / target_spacing))
    ]
    resampler = sitk.ResampleImageFilter()
    resampler.SetSize(new_size)
    resampler.SetOutputSpacing((target_spacing, target_spacing, target_spacing))
    resampler.SetInterpolator(sitk.sitkLinear)
    resampled = resampler.Execute(sitk_img)
    return sitk.GetArrayFromImage(resampled), (target_spacing, target_spacing, target_spacing)

失败时看什么：如果重采样后体积明显变化，先检查原始spacing是否被正确读取——尤其是层厚和层间距混在一起的情况，有些DICOM里 SpacingBetweenSlices 和 PixelSpacing 的第三个分量并不等价。

2.3 统一输出NIfTI：给模型和三维重建交一份标准数据

预处理最好统一到一个中间格式，我一般倾向输出 .nii.gz，原因有三个：保留spacing元数据，方便3D Slicer直接查看，避免反复读写DICOM时把方向信息搞丢。

def array_to_nifti(volume_np, spacing, output_path):
    sitk_img = sitk.GetImageFromArray(volume_np)
    sitk_img.SetSpacing((spacing[0], spacing[1], spacing[2]))
    sitk.WriteImage(sitk_img, output_path)

这里有一个容易忽略的点：如果有横断面扫描方向倾斜的大角度定位序列，最好直接用 SimpleITK 读取原始DICOM序列（ sitk.ImageSeriesReader ），它能自动处理方向cosine，再转成NIfTI。对于大多数常规腹部平扫，手工构造的坐标已经够用，但遇到定位像、增强多期时，用 ImageSeriesReader 更稳。统一成NIfTI后，后面的分割模型和三维重建都只跟这个文件打交道，省得每个模块都重复解析DICOM。

3. 分割选型与推理：为什么肝脏任务默认用3D U-Net，下载的权重怎么跑

3.1 模型骨架选型：从U-Net到nnU-Net的取舍

肝脏分割这个任务，目前主流方案基本都绕不开U-Net家族。2D U-Net可以逐层分割，但切片之间缺乏上下文，在肝脏与胆囊、胃壁粘连的层面经常产生不连续的结果。3D U-Net直接吃整个体数据块，能利用z轴相邻信息，但显存开销大；像LiTS挑战赛上的Top方案，很多是基于3D U-Net加级联结构。源码包里最常见的模型组织方式也就是这两种：一种是3D U-Net，另一种是nnU-Net的推理脚本。

我的建议是先从3D U-Net起步。nnU-Net虽然效果好，但它的预处理管线很重，包括自动spacing归一化、前景裁剪、数据增强策略，整套调完可能要折腾好几天。而3D U-Net的训练和推理相对轻量，对于只想复现分割结果的人来说，理解成本低得多。

别人权重的加载方式要看它的state_dict结构，第一件事是把模型打印出来：

model = UNet3D(in_channels=1, out_channels=2, base_channels=16)
state = torch.load('liver_seg_model.pt', map_location='cpu')
model.load_state_dict(state)

很多翻车现场都是因为训练时用了 nn.DataParallel ，state_dict的key前面带着 module. 前缀，直接load会报错。这时用 torch.load 后过滤一下key即可。

3.2 推理脚本：加载权重、滑动窗口、概率图后处理

肝脏CT体数据往往有 512×512×300 以上尺寸，一块整的3D输入放不进去，滑动窗口推理是最常用的方案。下面给一个完整最小脚本：

import torch
import numpy as np
import SimpleITK as sitk

def predict_volume(model, volume_np, device='cuda', patch_size=96, step=48):
    # volume_np: [D, H, W] 预处理后的0-1浮点数据
    model.eval()
    model.to(device)
    D, H, W = volume_np.shape[:3]
    output_prob = np.zeros((2, D, H, W), dtype=np.float32)
    count_map = np.zeros((D, H, W), dtype=np.float32)

    for z in range(0, D - patch_size + 1, step):
        for y in range(0, H - patch_size + 1, step):
            for x in range(0, W - patch_size + 1, step):
                patch = volume_np[z:z+patch_size, y:y+patch_size, x:x+patch_size]
                if patch.shape != (patch_size, patch_size, patch_size):
                    continue
                inp = torch.from_numpy(patch).unsqueeze(0).unsqueeze(0).float().to(device)
                with torch.no_grad():
                    logits = model(inp)  # [1, 2, patch_size, patch_size, patch_size]
                    prob = torch.softmax(logits, dim=1)[0].cpu().numpy()
                output_prob[:, z:z+patch_size, y:y+patch_size, x:x+patch_size] += prob
                count_map[z:z+patch_size, y:y+patch_size, x:x+patch_size] += 1.0

    # 重叠区域取平均概率，避免patch边缘出现接缝
    valid = count_map > 0
    output_prob[:, valid] /= count_map[valid][np.newaxis, ...]
    return output_prob

patch_size=96, step=48 意味着相邻patch有50%重叠，重叠区域做加权平均，这样消除边缘伪影。如果你显存只有8G，patch降到64，step降到32；如果显存充足（24G以上），patch可以到128，step也可以放大到64，减少一半推理时间。

注意这里 output_prob 是双通道概率图，通道1对应肝前景。不要直接用 argmax 出结果，后面还有全局阈值和连通域处理。

3.3 阈值与连通域：把模型输出变成医生能看懂的掩码

网络输出的是概率，转成二值mask时阈值一般取0.5，但肝脏分割有一个常见现象：肝脏占整个腹部体积的比例只有10%-15%，如果模型在部分层上不自信，概率值在0.3-0.6之间震荡，这时固定阈值很吃亏。我建议先看一遍概率分布，再决定阈值：

import matplotlib.pyplot as plt

prob_map = output_prob[1]
plt.hist(prob_map.flatten()[::1000], bins=50)
plt.title('liver foreground probability')
plt.show()

一个常用的经验值：如果模型训练充分，概率图会出现明显的双峰分布，取波谷位置作为阈值通常比0.5更稳。如果训练不足，概率图会挤在0.1以下，这时候返回去检查预处理是否匹配训练时的窗宽窗位。

拿到二值mask后立刻做最大连通域和孔洞填充，是医疗分割任务少不了的步骤：

from scipy.ndimage import label, binary_fill_holes

def postprocess_mask(prob_map, threshold=0.5):
    binary = (prob_map > threshold).astype(np.uint8)
    # 只保留最大连通域，过滤脾脏、胃壁误检
    structure = np.ones((3, 3, 3), dtype=np.uint8)
    labeled, num = label(binary, structure=structure)
    if num > 0:
        sizes = np.bincount(labeled.ravel())
        largest_label = np.argmax(sizes[1:]) + 1
        binary = (labeled == largest_label).astype(np.uint8)
    # 孔洞填充：肝内血管横断面在mask里可能形成空洞
    binary = binary_fill_holes(binary).astype(np.uint8)
    return binary

先提最大连通域是为了排除脾脏、胆囊壁等位置的假阳性，孔洞填充是为了让肝内血管不会在三维重建时形成网眼状空洞。但要注意，如果肿瘤本身在CT上呈低密度并且包膜完整，孔洞填充会把肿瘤内部空腔也填平，风险是肿瘤体积被低估。一个折中做法是只做闭运算（ binary_closing ），不直接fill holes。

4. 三维重建：从掩码到光滑STL模型的操作步骤

4.1 后处理：取最大连通域、填孔、闭运算

三维重建输入是上一步的肝脏mask，但直接把mask丢进重建算法，出来的网格会是毛糙的，表面布满钉刺和各种凹陷。原因很简单：逐体素的二值mask在斜面处有锯齿，肝脏与周围组织分割的边缘本身也很不规则。所以在重建前先把mask打磨一下：

from scipy.ndimage import binary_closing, binary_opening

def refine_mask(mask_np):
    # 先闭运算填补小的凹陷，再开运算去除表面毛刺
    struct = np.ones((3, 3, 3), dtype=np.uint8)
    refined = binary_closing(mask_np, structure=struct, iterations=2)
    refined = binary_opening(refined, structure=struct, iterations=1)
    return refined.astype(np.uint8)

这一小步对最终视觉效果影响极大。闭合运算可以想象成把“外壳”上的小凹痕填平，开运算像把表面突起磨掉一层。迭代次数不是越多越好，肝脏和胆囊分界处的细长结构会被过度开运算磨没，我一般 closing=2, opening=1 起步，看结果再调。

4.2 用Marching Cubes生成曲面并导出STL/OBJ

Marching Cubes是三维重建里最成熟的标准算法，scikit-image的实现就够用：

from skimage import measure
import trimesh

def mask_to_mesh(mask_np, spacing, level=0.5):
    verts, faces, normals, values = measure.marching_cubes(
        mask_np, level=level, spacing=spacing, gradient_direction='ascent'
    )
    mesh = trimesh.Trimesh(vertices=verts, faces=faces, vertex_normals=normals)
    mesh.fix_normals()
    return mesh

# 使用示例
mesh = mask_to_mesh(mask, spacing=(1.5, 1.5, 1.5))
mesh.export('liver.stl')

实现上的关键是 level=0.5 ：Marching Cubes等值面在插值时，0.5表示体素值为0和1的边界中点；如果你把mask的体素值归一化到 0.0~1.0，这个参数就是标准二值边界。spacing参数必须与预处理时的体素间距一致，否则模型会被拉长或压扁。

trimesh帮我们省了很多事：它能自动修复法线方向，导出STL、PLY、OBJ都支持。如果你想在ParaView或3D Slicer里做更精细的后续处理，导成 .ply 格式通常比STL更好，因为PLY可以存储顶点颜色和法线信息。

4.3 可选的表面平滑：防锯齿与小凹痕处理

Marching Cubes出来的网格通常有几十万到百万面片，体积大、渲染卡，而且表面略带台阶感。一个简单的拉普拉斯平滑能明显改善视觉效果：

mesh = mesh.simplify_quadric_decimation(target_count=200000)
mesh = mesh.laplacian_smooth(iterations=10, lamb=0.5)
mesh.export('liver_smooth.stl')

注意，这里有一个需要权衡的地方。平滑会改变体积，拉普拉斯算法会让整体表面向内收缩一点，迭代次数越多收缩越明显。如果論文里要报体积定量结果，用未经平滑或仅极小平滑的网格算；如果只是给临床医生看三维形态，用平滑后的模型更直观。我通常的做法是同时导两个文件， liver_raw.obj 用于量化和距离评估， liver_smooth.obj 用于视觉汇报。

如果对表面质量要求更高，可以使用vtk的 vtkWindowedSincPolyDataFilter 做体积保留型平滑，它的频率域滤波能减少体积收缩现象，但代价是计算量更大。对于肝脏这种形状相对规整的器官，trimesh的laplacian平滑已经足够。

5. 避坑笔记：肝脏CT分割与三维重建里5个常见踩坑现场

5.1 现象：DICOM读出来的HU值全是负数

第一次跑通代码的人最容易蒙圈：打印出来的pixel_array值在-1000左右徘徊，看起来不像组织。原因其实不复杂：CT的原始存储值不是HU，必须经过rescale变换。具体讲，有些设备的存储像素值是加了偏移的整数，比如把空气标成0而不是-1000。只要漏掉 RescaleSlope 和 RescaleIntercept ，后续所有归一化全部报废。解决方法是按前面代码那样，在读取时就完成 pixel_array * slope + intercept ，并且在预处理阶段用真实的HU值做窗截断。顺手检查一下 ds.RescaleIntercept 是否等于-1024，很多设备都是这个常数。

5.2 现象：预测全是背景，或肝脏区域漏了一半

推理时模型输出概率图整体偏低，阈值0.5切出来几乎全黑。我遇到过两次，一次是训练时用 [-200, 150] 窗，我推理时用了 [-100, 200] ，分布错位导致的；另一次是模型对某个增强期相的CT不适应，训练数据全是门静脉期，推理却给了动脉期。解决路线：先打印输入统计值，看min/max是否落在训练时的数值范围，再对比窗宽窗位设置，最后再考虑降阈值。这里特别提醒，如果是多期增强CT，不同期相的肝脏CT值范围差异很大，最好训练时就按期相拆分数据，或者把期相信息做成模型的额外输入。

5.3 现象：同一个病例在两个工作站上体积差20%

体积对不上，最常见的是spacing不一致。同一台CT扫描仪，重建间隔可能是1mm、2mm、5mm，不同中心的重建参数完全不同。如果用原始spacing直接mesh，两个来源的mask体积自然差一大截。解决方法是统一重采样到各向同性spacing再做重建。另一个隐蔽原因是：线性插值重采样会把边缘体素的值混到软阴影，导致分割结果边缘偏移半个体素，累积起来就是几个百分点的体积差异。这种情况推荐改用 sitk.sitkNearestNeighbor 对mask做重采样，保留硬边界。

5.4 现象：生成STL后模型被拉长或翻转

这是三维重建阶段最常见的翻车现场。Marching Cubes的spacing三个分量没有对应到X/Y/Z轴，或者DICOM读取时切片顺序颠倒，都会造成三维模型变形。另外还有一个极容易忽略的点：SimpleITK的 GetArrayFromImage 给出的数组顺序是 [z, y, x]，而多数三维重建库期望的坐标顺序是 [x, y, z]。如果直接把数组当作 [x, y, z] 输入，模型就会绕对角线转90度。我个人习惯是统一把所有操作放在SimpleITK的坐标系统里，以NIfTI文件的spacing和direction为基准，只在导出的最后一步做轴转换。

5.5 现象：GPU显存不够，推理直接OOM

512×512×300的体数据直接推理，即使8G显存也撑不住。滑动窗口是标准解法，但窗口大小和步长又要重新调参。我踩过一次教训：把patch_size设成96，显存刚好够，但推理速度和结果都正常；后来换了更大尺寸的体数据，忘了修改step参数，导致部分区域重叠不足，在z轴方向留下明显接缝。更稳的做法是记录valid count做加权平均，同时给三个方向设置循环边界保护，窗口别越过边界。如果你是CPU推理，patch大小降到64会更现实，速度虽然慢，但至少不崩溃。

6. 评估与验收：Dice、HD95和把STL放回CT原图核对

分割与重建的最终目标是拿到一个能用的模型文件，怎么算“能用”？两个硬指标加一个软检查。

硬指标一：Dice系数。把预测mask和标注mask做逐体素对比：

def dice_score(gt_np, pred_np):
    intersection = np.sum((gt_np > 0) & (pred_np > 0))
    return 2.0 * intersection / (np.sum(gt_np > 0) + np.sum(pred_np > 0))

肝脏分割公开数据集上的SOTA大致能到0.95左右；自己训练的模型，测试集上Dice稳定在0.90以上就算及格。硬指标二：HD95（95%分位Hausdorff距离），它衡量两个表面之间的最大空隙，用SimpleITK实现：

def surface_distance(gt_np, pred_np, input_spacing=(1.5, 1.5, 1.5)):
    gt_sitk = sitk.GetImageFromArray((gt_np > 0).astype(np.uint8))
    pred_sitk = sitk.GetImageFromArray((pred_np > 0).astype(np.uint8))
    gt_sitk.SetSpacing(input_spacing)
    pred_sitk.SetSpacing(input_spacing)
    # signed distance map需要float类型
    gt_dist = sitk.SignedMaurerDistanceMap(gt_sitk, useImageSpacing=True, insideIsPositive=False)
    pred_sitk_f = sitk.Cast(pred_sitk, sitk.sitkFloat32)
    pred_dist = sitk.SignedMaurerDistanceMap(pred_sitk_f, useImageSpacing=True, insideIsPositive=False)
    return sitk.HausdorffDistanceImageFilter().Execute(gt_dist > 0, pred_dist > 0).GetHausdorffDistance()

我一般要求HD95小于5mm，超过这个值说明边缘有系统性偏差，大概率是预处理阶段的问题。

软检查：导出STL后，在3D Slicer里同时载入原始CT和STL模型，检查模型是否贴合肝实质边界。这一步能发现Dice报告不出来的问题——比如模型整体偏上一两个体素、胆囊被吞进来一半、肝左叶被削掉一截。如果模型空间位置对不上，不要急着调表面平滑参数，先回头看重采样时有没有对齐origin。

我现在的习惯是把整条管线做成三个固定接口： preprocess(dicom_dir) 输出NIfTI、 predict(nii_path) 输出mask、 reconstruct(mask_path) 输出STL。新到一批数据，先用同一个preprocess跑一遍，再在Slicer里目测，程序不报错只代表数据流没问题，形态对得上才是真正验收通过。三维重建是个“结果一眼能看到好坏”的活，把前两个阶段的坑填平了，最后生成的模型大概率不会辜负你。希望帮到你。

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