尧图建网站 尧图建网站 YAOTU WEB BUILD 免费咨询
ARTICLE DETAIL

资讯详情

深耕网站建设与建站编程的一线实战洞察。

Python实现CT DICOM三维重建:从pydicom读取到VTK体绘制

Python实现CT DICOM三维重建:从pydicom读取到VTK体绘制 简介本资源是一份面向医学影像处理初学者与Python可视化开发者的三维重建实践代码包聚焦CT DICOM数据的读取、体素重建与VTK交互式3D渲染全流程。资源包含2个核心文件1个标准CT切片DICOM样本.dcm用于真实数据输入1个完整可运行的Python脚本.py内建pydicom解析逻辑、vtkImageData体素构建、vtkVolumeMapper体绘制及RenderWindow交互控制覆盖从元数据提取、空间坐标对齐到光照透明度调节等关键环节。压缩包仅1.01MB结构精简无冗余依赖开箱即用。已有4286人学习下载适合希望快速掌握医学图像三维可视化技术栈的开发者——不仅能直接复现临床级CT体绘制效果还可作为二次开发基础模板延伸支持色彩映射、GPU加速或病灶标注模块集成。 做了几年医学影像三维可视化我注意到一个经常误导新人的点一提到三维重建很多人第一反应是去算焦距、相机内参、对极几何。那是摄影测量和SLAM的路子。医学CT的DICOM数据里切片本身已经携带了精确的物理空间信息——像素间距、层厚、病人坐标——所以用Python做CT三维重建核心工作其实是两件事把一组DICOM切片读成三维体数据再用VTK把这个体数据渲染出来。这套pydicom numpy VTK的组合我前后在好几个项目里跑过从本地DICOM文件夹读数据到渲染出可交互的三维体绘制整个过程完整跑通并不复杂但坑确实不少。本文会把环境选型、核心代码、常见元数据陷阱和交互定制全部拆开讲适合刚接触医学影像处理的Python开发者也适合想把手里的DICOM数据快速看成一个三维模型的医生或科研人员。1. 在动手写代码之前DICOM文件到底能不能直接三维重建1.1 单个dcm是二维切片一组dcm才是三维数据DICOM里一张图像就是一张二维灰度图CT扫描天生是断层成像单张dcm里只有一个断层的像素矩阵。想重建三维必须拿到同一个扫描序列Series下的所有切片通常是几十到上千张按顺序叠起来才构成完整的三维体数据。很多人第一次下载DICOM发现文件夹里几百个文件文件名乱七八糟就开始懵。其实DICOM的组织逻辑非常规整患者Patient下可以有多个检查Study一个检查下可以有多个序列Series一个序列下有若干实例Instance即一张张切片我们做三维重建实际操作对象是一个Series。同一Series内部的切片在空间上是连续排列的这一点非常关键。1.2 为什么说CT重建不需要算焦距和内外参CT机和普通相机有本质区别。相机成像靠镜头把三维世界投影到二维传感器上所以从多张二维照片恢复三维结构需要标定焦距、畸变、位姿这是摄影测量和视觉SLAM做的事情。而CT数据在采集阶段就已经完成了从投影数据到体素分布的重建DICOM文件里存的不是投影图而是人体组织对X射线衰减系数的三维分布。这意味着CT数据自带一套物理坐标系。每一层的ImagePositionPatient0020,0032记录了图像左上角在病人坐标系中的位置ImageOrientationPatient0020,0037给出了方向余弦。我们要做的不是反推几何而是把这些元数据读出来告诉VTK每个体素在真实空间里占多大、间隔多少让它按真实物理比例渲染。所以在医学影像这个语境里所谓三维重建准确说是三维可视化——把已经存在的体数据用体绘制或表面重建的方式呈现出来。这个认知很重要否则你会在错误的路上浪费大量时间。1.3 数据从哪里来文件夹、网闸拷贝还是DICOM服务医院环境一般从PACS影像归档与通信系统导出或者直接把整个Series目录拷贝到本地。离线实验阶段我建议直接读文件夹不要一开始就碰DICOM网络协议。但如果你在真实项目里需要从PACS拉数据那就要走DICOM C-STORE / C-GET这类网络服务了。常见的服务端工具链有dcm4che、OrthancPython侧可以用pynetdicom写SCU。等本地渲染流程跑通了再往这个方向扩展也不迟。有一段时期大家还会提DICOM虚拟打印服务器主要是为了把影像推给胶片打印机或者做Web端预览属于比较老的方案现在日常开发很少用了。2. 环境准备与版本选型这一套组合我稳定跑了一年2.1 为什么选pydicom VTK而不是VTK自带的DICOM readerVTK自带了一个vtkDICOMImageReader表面上看一行代码就能读目录reader vtk.vtkDICOMImageReader() reader.SetDirectoryName(dcm_dir) reader.Update()但实际用下来这个类在真实临床数据上有几个比较头疼的问题。对比项vtkDICOMImageReaderpydicom 手动读取多Series识别目录里多个序列时容易串完全可控按SeriesInstanceUID筛选Rescale Intercept/Slope部分版本支持不稳定自己读取显式转换层间距计算经常直接用SliceThickness不校验根据ImagePositionPatient精确计算预处理不灵活无法任意裁剪/归一化先拿numpy数组想怎么处理都行窗宽窗位调节需要额外配置Transfer Function可控性极强说白了vtkDICOMImageReader适合快速看个结果但遇到真实的Hounsfield UnitCT值转换、层间距不准确、多序列目录这些情况它就力不从心了。自己用pydicom读然后把numpy数组喂给VTK整个过程每一步都清清楚楚出了问题也容易定位。2.2 安装步骤与验证直接pip安装三个库就行pip install pydicom vtk numpy我当前的组合是Python 3.9 VTK 9.2 pydicom 2.3 numpy 1.24跑了很久没出过兼容性问题。VTK 9.x的Python绑定已经相当成熟不需要自己编译。安装完先验证一下python -c import vtk; print(vtk.vtkVersion.GetVTKVersion()) python -c import pydicom; print(pydicom.__version__)如果第一行报ImportError但pip又显示已安装通常是OpenGL相关动态库缺失。Linux下可以这样修复sudo apt-get install libgl1-mesa-dev libglu1-mesa-devWindows下如果提示缺少OpenGL DLL先更新显卡驱动一般能解决。2.3 源码编译的坑vtk::mpi链接错误与PCL场景提醒如果你是pip安装可以跳过这一节。但如果你是从源码编译VTK尤其是作为一个依赖库被PCL点云库拉进来一起编那就会遇到一些奇葩的CMake错误。热搜词里那个the link interface of target vtk::mpi contains: mpi::mpi_c就是典型的MPI模块链接问题。这个报错的本质是VTK在配置时检测到MPI库把MPI相关target暴露到了链接接口里但你的编译器/链接器找不到对应的MPI C库路径。解决办法有两个方向如果你根本不需要MPI支持直接关闭它在CMake配置时加-DVTK_GROUP_ENABLE_MPINO如果你确实需要MPI指定正确的MPI路径-DMPI_C_COMPILER/path/to/mpicc -DMPI_CXX_COMPILER/path/to/mpicxx另外从PCL那条路过来的朋友要注意PCL编译时链接的VTK版本必须和运行时一致否则会出现Qt符号版本冲突、OpenGL上下文初始化失败之类的问题而且这类问题通常极其难排查。我的建议是能用预编译二进制就用预编译不要自己折腾VTK源码编译投入产出比太低了。3. 核心实现拆解从pydicom读取到体绘制的完整流程3.1 读取DICOM并提取关键元数据读一个DICOM文件用pydicom非常简单import pydicom ds pydicom.dcmread(path/to/slice.dcm, forceTrue) print(ds.PatientName) print(ds.PixelSpacing) # (行间距, 列间距)单位mm print(ds.SliceThickness) # 层厚单位mm print(ds.InstanceNumber) # 切片序号 print(ds.RescaleIntercept) # 通常为 -1024 print(ds.RescaleSlope) # 通常为 1这里有几个关键Tag是三维重建必须关注的我列了个表Tag分组标签名含义用途(0028,0030)PixelSpacing像素在行/列方向的物理间距设置体素在X/Y方向的尺寸(0018,0050)SliceThickness层厚层间距离的参考值(0020,0032)ImagePositionPatient图像左上角世界坐标计算层间距和图像原点(0020,0037)ImageOrientationPatient图像方向余弦判断切片是否倾斜(0028,1052)RescaleIntercept灰度重采样截距CT值HU转换(0028,1053)RescaleSlope灰度重采样斜率CT值HU转换(0020,0013)InstanceNumber实例序号切片排序CT的pixel_array一般是uint16甚至int16但真实物理含义需要经过公式转换hu ds.pixel_array * ds.RescaleSlope ds.RescaleIntercept转换后得到的才是Hounsfield UnitCT值空气约-1000水约0骨骼通常几百到一千以上。做体绘制时颜色函数和透明度函数直接作用在这个值上。3.2 切片堆叠成三维体数据的正确姿势读多个切片后堆叠顺序是个容易被忽视的细节。很多人按文件名排序但DICOM文件名是生成系统定义的和扫描顺序完全没有关系。正确做法是按InstanceNumber排序import os import pydicom import numpy as np def load_dicom_series(dcm_dir): ds_list [] for fname in os.listdir(dcm_dir): fpath os.path.join(dcm_dir, fname) try: ds pydicom.dcmread(fpath, forceTrue) # 只保留CT类或至少是图像类文件 if hasattr(ds, pixel_array): ds_list.append(ds) except Exception: continue ds_list.sort(keylambda x: x.InstanceNumber if hasattr(x, InstanceNumber) else 0) slices np.stack([ds.pixel_array for ds in ds_list]) # 此时 slices.shape (num_slices, rows, cols)对应 (z, y, x) return slices, ds_list更稳健的做法是按ImagePositionPatient的第三个分量Z坐标排序。有些扫描序列InstanceNumber可能不连续或缺失但每一层的空间Z坐标永远是准确的。def load_dicom_series_by_position(dcm_dir): ds_list [...] # 取每张切片的 ImagePositionPatient 的 Z 值 ds_list.sort(keylambda x: float(x.ImagePositionPatient[2])) ...这里还要注意层间距不一定是SliceThickness。有些扫描设备层间隔和层厚不一致最稳妥的计算方式是取相邻两层的ImagePositionPatient的Z坐标之差z_values [float(ds.ImagePositionPatient[2]) for ds in ds_list] slice_spacing abs(z_values[1] - z_values[0])只有拿到这个值VTK渲染出来的模型才不会被压扁或拉长。3.3 体绘制管线颜色函数、透明度函数、mapper的选择体绘制Volume Rendering本质上是让光线穿过三维体数据在路径上根据每个体素的颜色和不透明度累积颜色。VTK里最常用的是vtkGPUVolumeRayCastMapper它把光线投射算法跑在GPU上交互旋转时能保持流畅。管线组成如下vtkImageData存放三维体数据vtkGPUVolumeRayCastMapper把体数据映射为渲染所需的图元vtkColorTransferFunction把CT值映射为RGB颜色vtkPiecewiseFunction把CT值映射为不透明度vtkVolumeProperty把颜色和不透明函数绑定到Volume对象vtkVolume最终添加到渲染器的三维对象关键点是颜色函数和不透明函数的设置。CT值范围大约是-1024到3000但不同组织分布在不同区间空气约-1000肺组织约-800到-500脂肪约-100到-50软组织/肌肉约20到80骨骼300以上体绘制的效果基本完全取决于你在这几个CT值节点上怎么分配颜色和透明度。后面第4节会展开讲。3.4 完整可运行源码下面是可以直接复制运行的完整脚本我把注释写得比较详细。import os import pydicom import numpy as np import vtk from vtk.util.numpy_support import numpy_to_vtk def load_dicom_series(dcm_dir): 读取一个DICOM目录返回体数据numpy数组、间距和原始DS列表 ds_list [] for fname in os.listdir(dcm_dir): fpath os.path.join(dcm_dir, fname) try: ds pydicom.dcmread(fpath, forceTrue) if hasattr(ds, pixel_array): ds_list.append(ds) except Exception: continue if not ds_list: raise ValueError(没有找到可用的DICOM文件) ds_list.sort(keylambda x: float(x.ImagePositionPatient[2])) slices np.stack([ds.pixel_array for ds in ds_list]) # 间距 px, py [float(v) for v in ds_list[0].PixelSpacing] z0 float(ds_list[0].ImagePositionPatient[2]) z1 float(ds_list[1].ImagePositionPatient[2]) pz abs(z1 - z0) return slices, (px, py, pz), ds_list def numpy_to_vtk_image_data(slices, spacing): 把numpy体数据转成VTK ImageData # slices的shape是 (n_slices, rows, cols)对应 (z, y, x) # VTK希望内存布局是x最快变化所以转置为 (cols, rows, slices) # 注意这一步直接影响渲染方向做错会看到左右颠倒或上下颠倒 vol slices.transpose(2, 1, 0).astype(np.float32) vtk_data numpy_to_vtk(vol.ravel(), deepTrue, array_typevtk.VTK_FLOAT) image vtk.vtkImageData() cols, rows, n_slices vol.shape image.SetDimensions(cols, rows, n_slices) image.SetSpacing(spacing[0], spacing[1], spacing[2]) image.GetPointData().SetScalars(vtk_data) return image def make_volume(image): 创建体绘制Volume对象 # 颜色传递函数 color_func vtk.vtkColorTransferFunction() color_func.AddRGBPoint(-1024, 0.0, 0.0, 0.0) # 空气 color_func.AddRGBPoint(-300, 0.4, 0.2, 0.1) # 低密度组织 color_func.AddRGBPoint(40, 0.8, 0.5, 0.3) # 软组织 color_func.AddRGBPoint(300, 0.9, 0.8, 0.7) # 高密度 color_func.AddRGBPoint(1500, 1.0, 1.0, 1.0) # 骨 # 不透明度传递函数 opacity_func vtk.vtkPiecewiseFunction() opacity_func.AddPoint(-1024, 0.0) opacity_func.AddPoint(-300, 0.0) opacity_func.AddPoint(40, 0.4) opacity_func.AddPoint(300, 0.6) opacity_func.AddPoint(1500, 0.9) prop vtk.vtkVolumeProperty() prop.SetColor(color_func) prop.SetScalarOpacity(opacity_func) prop.ShadeOn() # 开启光照让表面有立体感 prop.SetInterpolationTypeToLinear() mapper vtk.vtkGPUVolumeRayCastMapper() mapper.SetInputData(image) mapper.SetBlendModeToComposite() volume vtk.vtkVolume() volume.SetMapper(mapper) volume.SetProperty(prop) return volume def main(): dcm_dir path/to/dicom_series slices, spacing, ds_list load_dicom_series(dcm_dir) print(f体数据 shape: {slices.shape}, spacing: {spacing}) image numpy_to_vtk_image_data(slices, spacing) volume make_volume(image) renderer vtk.vtkRenderer() renderer.AddVolume(volume) renderer.SetBackground(0.05, 0.05, 0.08) render_window vtk.vtkRenderWindow() render_window.AddRenderer(renderer) render_window.SetSize(900, 700) render_window.SetWindowName(CT DICOM 3D Volume Rendering) interactor vtk.vtkRenderWindowInteractor() interactor.SetRenderWindow(render_window) style vtk.vtkInteractorStyleTrackballCamera() interactor.SetInteractorStyle(style) interactor.Initialize() render_window.Render() interactor.Start() if __name__ __main__: main()这个脚本就是最简版的可运行Demo。跑起来之后鼠标左键拖动旋转视角滚轮缩放中键平移基本交互都有了。3.5 跑出来效果不佳怎么调第一次跑出来效果通常不会太好。常见情况是整体灰蒙蒙、目标组织不清晰或者背景噪点很多。这基本不是代码问题而是Transfer Function设置的问题。我的经验是先定CT值范围再调整不透明度函数最后才谈颜色。比如你想优先看骨骼就把40以下的CT值不透明度压到0300以上的抬到接近1你想看软组织就把-100到100区间的梯度拉大。体绘制调参没有银弹一份数据一套参数多试几次就熟了。4. 重建显示出问题先检查这三个看不见的元数据4.1 体素间距不设置人像会变细长或变扁平这是新手最容易踩的坑。如果不设置SpacingVTK默认每个体素是1x1x1的立方体而实际腹部CT的层间距通常是2.5mm到5mm平面内像素间距却只有0.6mm到0.8mm。这就意味着你看到的模型在Z方向被拉伸了好几倍整个人像一根竹竿。正确设置后体数据在物理空间中的尺寸才是真实的。以512x512x300的体数据为例像素间距0.7mm、层间距2.5mm那么物理尺寸大约是X 512 * 0.7 358.4mm Y 512 * 0.7 358.4mm Z 300 * 2.5 750mm如果不设置SpacingVTK会把它当512x512x300的立方体渲染Z方向被拉得很夸张。判断比例是否正确有个小技巧渲染出来观察人体的头脚方向和左右方向的比例正常成年人躯干左右径和前后径相似但头脚方向远长于左右径。如果你看到一个人变成了横向的板砖那就是层间距和平面间距没有分清楚。4.2 窗宽窗位不处理渲染出来一片白CT值的有效范围很宽从空气的-1024到高密度骨骼的2000以上。如果不做映射直接把CT值除以最大值转成0到1那么大部分软组织的值会落在0.1附近整体看起来就是一团黑骨骼附近又一片白。这就是窗宽窗位Window Level / Window Width的核心作用。体绘制里的Transfer Function本质上就是窗宽窗位的高级版本——你不仅要决定显示范围还要给不同范围的CT值分配颜色和透明度。不同扫描部位常用的窗宽窗位参考窗类型窗位 (WL)窗宽 (WW)适用场景肺窗-6001500观察肺实质、气道腹部窗40400观察肝、脾、肾等软组织骨窗6002000观察骨骼结构脑窗4080观察脑组织在体绘制里我会先根据目标组织确定要显示的CT值区间然后在不透明度函数里把区间外的值压到0区间内的值按梯度提升。这样做出来的模型才干净。4.3 大序列内存撑爆怎么低开销加载医学影像数据量不容小觑。512x512x500的16位CT序列裸数据就是512x512x500x2字节约256MB。如果读进来后转成float64直接飙到1GB以上。很多人的机器不是撑不住渲染而是在转数据类型这一步就卡死了。推荐的处理策略读取时保持原始dtype用int16或uint16需要做HU转换时用float32而不是float64渲染前如果不需要精确HU值直接归一化到uint8使用numpy的memmap做内存映射读取避免一次性载入一个更激进的做法是直接以uint16类型喂给VTK颜色函数和透明度函数都基于原始值定义。这样内存占用最小而且不会丢失CT值的精度。# 内存优化版本 slices np.stack([ds.pixel_array for ds in ds_list]) # 保持原始dtype vol slices.transpose(2, 1, 0) # 不做astype保持原样 vtk_data numpy_to_vtk(vol.ravel(), deepTrue, array_typevtk.VTK_SHORT)5. 交互相当重要鼠标坐标读取、快捷键与观察体验5.1 默认交互器已经够用但需要调整vtkInteractorStyleTrackballCamera是三维观察的默认选择左键旋转、滚轮缩放、中键平移跟大多数三维软件一致。但医学影像场景有点特殊医生习惯了PACS系统里那种横断面/冠状面/矢状面切换的阅片模式直接甩一个三维旋转视角过去很多人会迷路。实际开发中我一般同时提供两个模式三维浏览模式vtkInteractorStyleTrackballCamera给诊断和展示用切片浏览模式vtkInteractorStyleImage配合鼠标滚轮逐层切换切片切换交互器样式在VTK里非常直接interactor.SetInteractorStyle(vtk.vtkInteractorStyleImage())5.2 获取鼠标坐标的两种实现热搜词里vtk获取鼠标坐标出现频率很高这里展开聊聊。VTK里获取鼠标位置有很多种picker最常用的是vtkWorldPointPicker和vtkPropPicker。两者的核心区别在于vtkWorldPointPicker从屏幕像素位置出发沿着视锥方向求与深度缓冲区的交点返回世界坐标。不需要点击到某个Actor上就能获取位置适合读任意位置的空间点vtkPropPicker只在场景中的Actor表面拾取能判断你点到了哪个Actor适合做交互选择class MedicalInteractorStyle(vtk.vtkInteractorStyleTrackballCamera): def __init__(self): super().__init__() self.picker vtk.vtkWorldPointPicker() self.renderer None def OnLeftButtonDown(self): x, y self.GetInteractor().GetEventPosition() self.picker.Pick(x, y, 0, self.renderer) world_pos self.picker.GetPickPosition() print(世界坐标:, world_pos) super().OnLeftButtonDown()拿到了世界坐标如果想换算成体数据索引体素坐标公式很简单index [ int(round((world_pos[0] - origin[0]) / spacing[0])), int(round((world_pos[1] - origin[1]) / spacing[1])), int(round((world_pos[2] - origin[2]) / spacing[2])) ]这里的origin是体数据的原点坐标。如果直接用vtkImageImport而没有显式设置Origin默认是0,0,0。不过我提醒一句vtkWorldPointPicker返回的是渲染面上的三维点不是光线进入体数据的入口。如果要做精确的体素拾取需要自己做ray casting求交或者用vtkCellPicker拾取一个代理几何体来表达位置。这是很多人实现点击体数据某个位置反馈CT值时踩的深坑。5.3 自定制快捷键窗宽窗位、复位视角为了让渲染结果更像一个临床工具我通常会绑定几个快捷键。最实用的是用加减号调窗宽、左右括号调窗位R键复位视角。class MedicalInteractorStyle(vtk.vtkInteractorStyleTrackballCamera): def __init__(self, volume, color_func, opacity_func): super().__init__() self.volume volume self.color_func color_func self.opacity_func opacity_func self.window 400 self.level 40 def OnKeyPress(self): key self.GetInteractor().GetKeySym() if key plus or key equal: self.window 50 self._update_transfer_func() elif key minus: self.window - 50 self._update_transfer_func() elif key bracketleft: self.level - 50 self._update_transfer_func() elif key bracketright: self.level 50 self._update_transfer_func() elif key r or key R: self.GetInteractor().GetRenderWindow().GetRenderers().GetFirstRenderer().ResetCamera() self.GetInteractor().GetRenderWindow().Render() super().OnKeyPress()更新颜色函数和不透明度函数时要小心不要每次都重建整条传递函数曲线。比较高效的做法是只更新关键节点的位置然后触发重新渲染。def _update_transfer_func(self): lower self.level - self.window / 2 upper self.level self.window / 2 # 这里根据新的窗宽窗位重建颜色/不透明度节点 # 核心是让当前显示范围映射到可见区间 self.GetInteractor().GetRenderWindow().Render()6. 二次开发方向表面重建、PACS联动与批处理6.1 从体绘制切换到表面重建体绘制适合观察整体密度分布但如果你要导出一个模型去做3D打印或者定量测量就必须做表面重建。VTK里最常用的是Marching Cubes算法VTK 8.2之后推荐用vtkFlyingEdges3D速度更快、内存更省。以提取骨骼表面为例import vtk iso vtk.vtkFlyingEdges3D() iso.SetInputData(image) iso.SetValue(0, 300) # 提取CT值300以上的等值面 iso.ComputeNormalsOn() mapper vtk.vtkPolyDataMapper() mapper.SetInputConnection(iso.GetOutputPort()) actor vtk.vtkActor() actor.SetMapper(mapper) # 输出为STL文件用于3D打印或导入其他软件 stl_writer vtk.vtkSTLWriter() stl_writer.SetFileName(bone.stl) stl_writer.SetInputConnection(iso.GetOutputPort()) stl_writer.Write()阈值的选择直接影响重建质量。阈值太高骨头会碎成渣太低会把软组织也包进来。我一般会先通过体绘制观察一下数据的大致分布再结合直方图选择一个合理的阈值。不同部位、不同扫描参数下骨骼的CT值可能有差异但通常在200到400之间。6.2 和PACS/本地DICOM服务的对接思路如果你不满足于读文件夹想直接从PACS拉数据那就需要走DICOM网络协议。最实用的方案是用pynetdicom实现一个C-STORE SCU向远程PACS发起查询和获取请求。from pynetdicom import AE, debug_logger from pynetdicom.sop_class import CTImageStorage # debug_logger() ae AE() ae.add_requested_context(CTImageStorage) assoc ae.associate(127.0.0.1, 11112, ae_titleMY_STORE_SCU) if assoc.is_established: # 发送DICOM文件 ds pydicom.dcmread(slice.dcm) status assoc.send_c_store(ds) assoc.release()这个方向适合做数据自动化流转比如新扫描一出来自动把DICOM数据拉到后处理服务器上做三维重建。另外Orthanc这个轻量级DICOM服务器也值得关注它自带REST API可以把DICOM实例转成JSON/PNG甚至直接集成VTK做Web端预览极大降低和医院系统对接的成本。6.3 批量处理与工程化建议当你的DICOM目录里同时存在多个Series时必须先按SeriesInstanceUID分组再对每个Series分别做重建。这一步很多人会漏掉结果就是把两个不同扫描序列混在一起堆叠渲染出来的东西完全错乱。from collections import defaultdict series_dict defaultdict(list) for fname in os.listdir(dcm_dir): ds pydicom.dcmread(os.path.join(dcm_dir, fname), forceTrue) if hasattr(ds, pixel_array): series_dict[ds.SeriesInstanceUID].append(ds) for series_uid, ds_list in series_dict.items(): ds_list.sort(keylambda x: float(x.ImagePositionPatient[2])) slices np.stack([ds.pixel_array for ds in ds_list]) # 对每个Series进行重建和保存数据量大时我习惯先把处理好的体数据保存为.npy或者.vti格式后续调试渲染参数时直接加载不用每次重新解析几百个DICOM文件。VTK自带的vtkXMLImageDataWriter写.vti文件读取用vtkXMLImageDataReader效率很高。最后再说一个实际开发中的体会。很多人觉得医院影像工作站里那种高大上的三维效果靠的是神秘算法其实不是大部分效果差异来自窗宽窗位的把握、Transfer Function的细腻调节和交互细节的打磨。你如果在本地能把体绘制参数调到得心应手后面接PACS服务、做自动分割、训练分割模型整个思路都会顺畅很多。这套Python VTK的流程前期投入不大后续扩展空间却非常广值得好好花点时间跑通它。本文还有配套的精品资源点击获取
返回列表