Reformat axial CT sequence as sagittal and coronal

[alancam]

Hello,

I have a thin slice axial CT that I would like to reformat in the orthogonal sagittal and coronal planes.

I've been trying to combine snippets from multiple places using pydicom but I just can't get success. I'm very much a beginner!

The axial images a a file per slice, which is what I would also want my reformats/reconstructions to be.

I've gotten to the below so far but receive an error when I try to save the slice files using dcmwrite ('DataElement' object has no attribute 'is_implicit_VR'). If I try to change this line to
pydicom.filewriter.dcmwrite(os.path.join(dir_path, filename), coronal_ds[:,slice,:], write_like_original=False)
(as I'm assuming I am writing the wrong thing) I also get an error (KeyError "'(0008, 0008) Image Type CS: ['Derived', 'Secondary']'" TypeError: '>' not supported between instances of 'DataElement' and 'int').

I think it's patently obvious that my understanding is very minimal!

Any guidance would be much appreciated.

import pydicom
import numpy as np
import sys
import os

# Load the DICOM files in the directory
dir_path = 'G:/RAIQC Nodules/Training/TG60 - Copy/20190208/CT/'
files = [pydicom.dcmread(dir_path + f) for f in os.listdir(dir_path) if f.endswith('.dcm')]
print("file count: {}".format(len(files)))

print("file count: {}".format(len(files)))

# skip files with no SliceLocation (eg scout views)
axial_data = []
skipcount = 0
for f in files:
    if hasattr(f, 'SliceLocation'):
        axial_data.append(f)
    else:
        skipcount = skipcount + 1

print("skipped, no SliceLocation: {}".format(skipcount))

# ensure they are in the correct order
axial_data = sorted(axial_data, key=lambda s: s.SliceLocation)

# pixel aspects, assuming all slices are the same
ps = axial_data[0].PixelSpacing
ss = axial_data[0].SliceThickness
ax_aspect = ps[1]/ps[0]
sag_aspect = ps[1]/ss
cor_aspect = ss/ps[0]

# create 3D array
img_shape = list(axial_data[0].pixel_array.shape)
img_shape.append(len(axial_data))
img3d = np.zeros(img_shape)

# fill 3D array with the images from the files
for i, s in enumerate(axial_data):
    img2d = s.pixel_array
    img3d[:, :, i] = img2d


#sag
sagittal_ds = files[0].copy()
sagittal_ds.PixelData = img3d[:, img_shape[1]//2, :]
#cor
coronal_ds = files[0].copy()
coronal_ds.PixelData = img3d[img_shape[0]//2, :, :].T

# Get an array ready for images
images = []
i = 0

# Save coronal data

coronal_ds.Rows, coronal_ds.Columns = coronal_ds.PixelData.shape
coronal_ds.ImagePositionPatient = [0, 0, 0]
coronal_ds.SliceThickness = axial_data[0].SliceThickness
coronal_ds.PhotometricInterpretation = 'RGB'
coronal_ds.SOPInstanceUID = pydicom.uid.generate_uid()
coronal_ds.ImageType = ['Derived', 'Secondary']

for i, slice in enumerate(coronal_ds):
    # Set the filename for the new DICOM file
    filename = f'coronal_{i+1}.dcm'

    # Save the DICOM file to disk
    pydicom.filewriter.dcmwrite(os.path.join(dir_path, filename), coronal_ds[:,slice,:], write_like_original=False)




# Reset array ready for images
images = []
i = 0

# Save sagittal data

sagittal_ds.Rows, sagittal_ds.Columns = sagittal_ds.PixelData.shape
sagittal_ds.ImagePositionPatient = [0, 0, 0]
sagittal_ds.SliceThickness = axial_data[0].SliceThickness
sagittal_ds.PhotometricInterpretation = 'RGB'
sagittal_ds.SOPInstanceUID = pydicom.uid.generate_uid()
sagittal_ds.ImageType = ['Derived', 'Secondary']

for i, slice in enumerate(sagittal_ds):
    # Set the filename for the new DICOM file
    filename = f'sagittal_{i+1}.dcm'

    # Save the DICOM file to disk
    pydicom.filewriter.dcmwrite(os.path.join(dir_path, filename), sagittal_ds[:,slice,:], write_like_original=False)