使用Python、Numpy和Matplotlib绘制3D医学图像失败求助

我来帮你排查下这个3D医学图像渲染失败的问题，从你的代码和需求来看，主要可能有几个常见的坑，咱们一步步来解决：

1. 先确认DICOM数据读取的正确性 #

你用imageio.volread读取DICOM序列虽然可行，但很容易踩文件排序错误的坑：os.listdir返回的文件名顺序不一定是DICOM切片的实际扫描顺序（比如按切片z轴位置排序），这会导致读出来的3D体素数据混乱，最终渲染不出正确结构。

更可靠的方式是用pydicom来读取并正确排序切片，示例代码如下：

import pydicom
import numpy as np
import matplotlib.pyplot as plt
from pydicom.fileset import FileSet

def make_ax(grid=False):
    fig = plt.figure(figsize=(10, 8))
    # 替换弃用的gca写法，用add_subplot创建3D轴
    ax = fig.add_subplot(projection='3d')
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    ax.grid(grid)
    return ax

# 读取并排序DICOM序列
dicom_path = 'Path to .dcm images'
fs = FileSet(dicom_path)
slices = list(fs)
# 按切片的z轴位置排序，确保体素顺序正确
slices.sort(key=lambda x: x.ImagePositionPatient[2])
# 转换为3D numpy数组
vol = np.stack([slice_.pixel_array for slice_ in slices])
print(f"体素数据形状: {vol.shape}")

2. 优化voxels函数的使用 #

ax.voxels默认只渲染非零体素，而医学图像大部分区域都是背景（像素值为0），直接渲染可能看不到任何结构。另外，大体积数据会导致渲染性能极差，甚至崩溃。可以做这两个优化：

• 设置阈值过滤，只渲染有意义的组织区域
• 对数据降采样，减少体素数量提升渲染速度

示例代码：

ax = make_ax(True)
# 设置阈值，过滤掉背景（根据你的图像调整数值）
threshold = 50
# 降采样，每隔2个取一个体素（可根据数据大小调整步长）
downsampled_vol = vol[::2, ::2, ::2]
# 渲染大于阈值的体素，添加面颜色让结构更明显
ax.voxels(downsampled_vol > threshold, 
          edgecolors='gray', 
          facecolors='lightblue',
          alpha=0.5)  # 透明度参数，避免遮挡
plt.show()

3. 检查维度兼容性 #

确认vol的维度是(z, y, x)，这是ax.voxels要求的输入形状（对应3D空间的z轴切片数、y轴高度、x轴宽度）。如果你的数据维度是(x, y, z)，需要转置调整：

vol = vol.transpose(2, 0, 1)

整合后的完整代码 #

把上面的优化点整合起来，最终代码如下：

import pydicom
import numpy as np
import matplotlib.pyplot as plt
from pydicom.fileset import FileSet

def make_ax(grid=False):
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(projection='3d')
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    ax.grid(grid)
    return ax

# 读取DICOM序列
dicom_path = 'Path to .dcm images'
fs = FileSet(dicom_path)
slices = list(fs)
slices.sort(key=lambda x: x.ImagePositionPatient[2])
vol = np.stack([slice_.pixel_array for slice_ in slices])

# 预处理数据
threshold = 50
downsampled_vol = vol[::2, ::2, ::2]

# 渲染3D体素
ax = make_ax(True)
ax.voxels(downsampled_vol > threshold, 
          edgecolors='gray', 
          facecolors='lightblue',
          alpha=0.5)
plt.show()

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