使用pydicom从轴向视图中提取矢状和冠状切口 [英] Extract sagittal and coronal cuts from axial view using pydicom

问题描述

我正在尝试读取一系列.dcm文件，这些文件默认情况下显示轴向视图.下面是代码:

I am trying to read a series of .dcm files which are by default show axial view. Below is the code:

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. 所以我的问题是，在轴向切开的情况下，如何获得矢状切面和冠状切面并绘制它们?

As you can see from the code I could not figure out the relevant idx and idy for particular dcm file. So my question is how to get sagittal and coronal cuts and plot them, given the axial cuts?

先谢谢了.

正如@ColonelFazackerley回答得很完美.我在下面添加一行，只是为了说明如何使用它.

As @ColonelFazackerley answered perfectly. I am adding below line just to show how I used it.

# 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

tested against series 2 from this example 3D CT data http://www.pcir.org/researchers/54879843_20060101.html

编辑说明:接受作为pydicom项目的示例 https://github.com/pydicom/pydicom/blob/master/examples/image_processing/reslice.py

edit note: accepted as an example into the pydicom project https://github.com/pydicom/pydicom/blob/master/examples/image_processing/reslice.py

这篇关于使用pydicom从轴向视图中提取矢状和冠状切口的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！