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)