Search code examples
pythondicompydicom

Extract sagittal and coronal cuts from axial view using pydicom


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)

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?

Thanks in advance.

Edit: 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)  

Solution

  • """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()
    

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

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