Python中大型3D DICOM CT图像带权重动态滤波下采样的高效方法

我来帮你梳理一下针对CT DICOM序列进行Z轴整数比例下采样的实现思路，结合你提到的π相关线性权重映射需求，下面是具体的步骤和代码示例：

1. 先加载DICOM序列为3D数组 #

首先需要把一系列2D DICOM切片转换成标准的3D numpy数组，这里用pydicom和numpy实现最方便：

import pydicom
import numpy as np
from pathlib import Path

def load_dicom_series(dicom_dir):
    # 读取文件夹下所有DICOM文件并按Z轴位置排序
    dicom_paths = sorted(list(Path(dicom_dir).glob("*.dcm")))
    slices = [pydicom.dcmread(path) for path in dicom_paths]
    slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))
    
    # 堆叠成(H, W, Z)形状的3D数组（和你描述的origimg[x,y,z]维度匹配）
    img_3d = np.stack([s.pixel_array for s in slices], axis=-1)
    return img_3d, slices[0]  # 返回图像数组和元数据模板，方便后续保存

2. 核心：带线性权重的Z轴下采样逻辑 #

针对整数比例下采样（比如160→80，步长为2；160→40，步长为4），我们封装一个通用函数，支持自定义权重规则：

def z_downsample_with_weight(img_3d, target_z, weight_func):
    """
    对3D DICOM图像进行Z轴整数比例下采样
    参数:
        img_3d: 输入3D数组，形状为 (H, W, orig_z)
        target_z: 目标Z轴长度，需满足 orig_z % target_z == 0（整数比例）
        weight_func: 权重函数，输入原Z层索引窗口，输出对应权重数组
    返回:
        downsampled_img: 下采样后的3D数组，形状为 (H, W, target_z)
    """
    orig_h, orig_w, orig_z = img_3d.shape
    step = orig_z // target_z
    downsampled_img = np.zeros((orig_h, orig_w, target_z), dtype=img_3d.dtype)
    
    for new_z in range(target_z):
        # 定位原序列中对应的Z层窗口
        start_idx = new_z * step
        orig_z_window = np.arange(start_idx, start_idx + step)
        # 获取自定义权重并归一化（避免亮度偏差）
        weights = weight_func(orig_z_window)
        weights = weights / np.sum(weights)
        # 加权计算当前目标Z层的像素值
        downsampled_img[:, :, new_z] = np.sum(img_3d[:, :, orig_z_window] * weights, axis=-1)
    
    return downsampled_img

3. 实现π相关的线性权重函数 #

你提到权重和π相关，这里举一个基于三角函数的线性权重示例（可以根据你的具体映射规则调整）：

def pi_based_linear_weight(z_window):
    # 计算窗口内的相对位置（0到1）
    rel_pos = np.arange(len(z_window)) / (len(z_window)-1) if len(z_window) > 1 else np.array([1.0])
    # 基于π的正弦加权，越靠近窗口后半段权重越高（可根据需求修改公式）
    weights = np.sin(rel_pos * np.pi)
    return weights

如果是简单的相邻层线性平均（比如160→80时，每层取原两层各50%权重），也可以写成：

def simple_linear_weight(z_window):
    return np.full(len(z_window), 1/len(z_window))

4. 完整流程示例 #

# 加载原始DICOM序列
dicom_dir = "./your_ct_dicom_folder"
img_3d, ref_dicom = load_dicom_series(dicom_dir)
print(f"原始图像形状: {img_3d.shape}")  # 输出类似 (512, 512, 160)

# 下采样到Z=80
target_z = 80
downsampled_img = z_downsample_with_weight(img_3d, target_z, pi_based_linear_weight)
print(f"下采样后图像形状: {downsampled_img.shape}")  # 输出 (512, 512, 80)

# 可选：保存下采样后的DICOM序列
def save_dicom_series(img_3d, ref_dicom, output_dir):
    Path(output_dir).mkdir(exist_ok=True)
    slice_thickness = float(ref_dicom.SliceThickness)
    step = orig_z // target_z
    for z in range(img_3d.shape[-1]):
        ds = ref_dicom.copy()
        ds.PixelData = img_3d[:, :, z].tobytes()
        ds.Rows, ds.Columns = img_3d.shape[:2]
        # 更新Z轴位置信息，保持物理空间连续性
        ds.ImagePositionPatient = list(ds.ImagePositionPatient[:2]) + [float(z) * slice_thickness * step]
        ds.save_as(Path(output_dir) / f"slice_{z:03d}.dcm")

save_dicom_series(downsampled_img, ref_dicom, "./downsampled_ct")

关键注意点 #

• 确保target_z是原Z轴长度的约数，否则需要额外处理非整数比例的插值逻辑
• 权重函数可以完全根据你的π相关映射需求修改，只要输出和输入窗口长度匹配的权重数组即可
• DICOM的像素值通常对应HU值，加权时建议保留原数据类型（比如int16），避免精度损失

内容的提问来源于stack exchange，提问作者Sean