python - 使用 pydicom 从轴向视图中提取矢状和冠状切面
问题描述
我正在尝试阅读一系列默认显示轴向视图的 .dcm 文件。下面是代码:
import os
import numpy as np
import pydicom as dicom
from matplotlib import pyplot as plt
root_dir = 'mydcomDir'
def sortDcm():
print('Given Path to the .dcm directory is: {}'.format(root_dir))
slices = [dicom.read_file(root_dir + '/' + s) for s in os.listdir(root_dir)]
slices.sort(key = lambda x: float(x.ImagePositionPatient[2]))
pos1 = slices[int(len(slices)/2)].ImagePositionPatient[2]
pos2 = slices[(int(len(slices)/2)) + 1].ImagePositionPatient[2]
diff = pos2 - pos1
# if diff > 0:
# slices = np.flipud(slices)
try:
slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
except:
slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
for s in slices:
s.SliceThickness = slice_thickness
# print("from sorted dicom",len(slices))
return slices
dcms = sortDcm()
ref_dicom = dcms[0]
d_array = np.zeros((ref_dicom.Columns,ref_dicom.Rows, len(dcms)), dtype=ref_dicom.pixel_array.dtype)
for dcm in dcms:
d_array[:, :, dcms.index(dcm)] = dcm.pixel_array
# fig = plt.figure(figsize=(12,12))
# plt.subplot(1, 3, 1)
# plt.title("Coronal")
# plt.imshow(np.flipud(d_array[idx , :, :].T))
# plt.subplot(1, 3, 2)
# plt.title("Sagital")
# plt.imshow(np.flipud(d_array[:, idy, :].T))
# plt.subplot(1, 3, 3)
plt.title("axial")
plt.imshow(d_array[:, :, dcms.index(dcm)])
plt.pause(0.001)
正如您从代码中看到的那样,我无法找出特定 dcm 文件的相关 idx 和 idy。所以我的问题是,考虑到轴向切割,如何获得矢状和冠状切割并绘制它们?
提前致谢。
编辑:正如@ColonelFazackerley 完美回答的那样。我在下面添加只是为了展示我是如何使用它的。
# fill 3D array with the images from the files
for i, s in enumerate(slices):
img2d = s.pixel_array
img3d[:,:,i] = img2d
#then to view sagittal and coronal slices for each of the axial slice
for i, s in enumerate(slices):
img2d = s.pixel_array
img3d[:,:,i] = img2d
corId = corId-1
sagId = sagId-1
# plot 3 orthogonal slices
a1 = plt.subplot(1,3,1)
plt.title('Axial')
plt.imshow(img3d[:,:,i],'gray')
a1.set_aspect(ax_aspect)
a2 = plt.subplot(1,3,2)
plt.title('Sagittal')
plt.imshow(np.flipud(img3d[:,sagId,:].T),'gray')
a2.set_aspect(sag_aspect)
a3 = plt.subplot(1,3,3)
plt.imshow(np.flipud(img3d[corId,:,:].T),'gray')
a3.set_aspect(cor_aspect)
plt.title('Coronal')
plt.show()
plt.pause(0.001)
解决方案
"""usage: reslice.py <glob>
where <glob> refers to a set of DICOM image files.
Example: python reslice.py "*.dcm". The quotes are needed to protect the glob
from your system and leave it for the script."""
import pydicom
import numpy as np
import matplotlib.pyplot as plt
import sys
import glob
# load the DICOM files
files=[]
print('glob: {}'.format(sys.argv[1]))
for fname in glob.glob(sys.argv[1], recursive=False):
print("loading: {}".format(fname))
files.append(pydicom.read_file(fname))
print("file count: {}".format(len(files)))
# skip files with no SliceLocation (eg scout views)
slices=[]
skipcount=0
for f in files:
if hasattr(f, 'SliceLocation'):
slices.append(f)
else:
skipcount = skipcount + 1
print("skipped, no SliceLocation: {}".format(skipcount))
# ensure they are in the correct order
slices = sorted(slices, key=lambda s: s.SliceLocation)
# pixel aspects, assuming all slices are the same
ps = slices[0].PixelSpacing
ss = slices[0].SliceThickness
ax_aspect = ps[1]/ps[0]
sag_aspect = ps[1]/ss
cor_aspect = ss/ps[0]
# create 3D array
img_shape = list(slices[0].pixel_array.shape)
img_shape.append(len(slices))
img3d=np.zeros(img_shape)
# fill 3D array with the images from the files
for i, s in enumerate(slices):
img2d = s.pixel_array
img3d[:,:,i] = img2d
# plot 3 orthogonal slices
a1 = plt.subplot(2,2,1)
plt.imshow(img3d[:,:,img_shape[2]//2])
a1.set_aspect(ax_aspect)
a2 = plt.subplot(2,2,2)
plt.imshow(img3d[:,img_shape[1]//2,:])
a2.set_aspect(sag_aspect)
a3 = plt.subplot(2,2,3)
plt.imshow(img3d[img_shape[0]//2,:,:].T)
a3.set_aspect(cor_aspect)
plt.show()
针对此示例 3D CT 数据中的系列 2 进行了测试 http://www.pcir.org/researchers/54879843_20060101.html
编辑说明:接受作为 pydicom 项目的示例 https://github.com/pydicom/pydicom/blob/master/examples/image_processing/reslice.py
推荐阅读
- autoit - 如何点击类下的可见文本?
- java - “WARN osscbcrypt.BCryptPasswordEncoder - Spring Security 中的空编码密码”是什么意思?
- jquery - 如何修复 Chrome 中的页面滚动延迟
- f# - 绑定中的意外符号“<-”
- css - 在正文中在背景图像上方显示带有渐变的日落效果
- javascript - 如何通过 url 共享购物车详细信息
- bash - 我怎样才能得到最后一个在遥控器上提交的人的名字?
- python - python: spyder: 运行 ipdb 时无法在控制台中编辑先前的命令
- python - ValueError:没有返回 HttpResponse 对象。它返回 None 而不是
- java - Java 正则表达式:双引号和字边界问题