DICOM (CT and RTDOSE) mapping in matplotlib

ct-rtdose-plot.py

"""
Step 1 - Read Dose and CT .dcm files
Step 2 - Resample dose to CT
Step 3 - Plot
"""
import pdb
import pydicom
import traceback
import numpy as np
import SimpleITK as sitk
from pathlib import Path

if __name__ == "__main__":

    try:
        
        pathCT     = Path('...', '{CTFolder}')
        pathDoseA5 = Path('...', '{RTDOSE}.dcm')
        sliceId = 94 # [77, 94, 115]

        if Path(pathCT).exists() and Path(pathDoseA5).exists():
            
            doseArray, ctArray     = [], []
            doseSpacing, ctSpacing = [], []
            doseImagePositionPatient, ctImagePositionPatientMin, ctImagePositionPatientMax = [], [], []

            # Step 1.1 - Read Dose .dcm
            try:
                dsDose      = pydicom.dcmread(pathDoseA5)
                doseArray   = dsDose.pixel_array * dsDose.DoseGridScaling
                doseSpacing = [float(each) for each in dsDose.PixelSpacing] + [float(dsDose.SliceThickness)]
                assert doseArray.shape[0] == float(dsDose.NumberOfFrames)
                doseImagePositionPatient = [float(each) for each in dsDose.ImagePositionPatient]
            except:
                traceback.print_exc()
                pdb.set_trace()
            
            # Step 1.2 - Read CT .dcms
            try:
                dsCTs = []
                for pathCTDicom in Path(pathCT).iterdir(): dsCTs.append(pydicom.dcmread(pathCTDicom))
                ctSpacing = [float(each) for each in dsCTs[0].PixelSpacing] + [float(dsCTs[0].SliceThickness)]
                dsCTs.sort(key=lambda x: x.InstanceNumber)
                ctImagePositionPatientMin = [float(each) for each in dsCTs[0].ImagePositionPatient]
                ctImagePositionPatientMax = [float(each) for each in dsCTs[-1].ImagePositionPatient]
                for dsCT in dsCTs:
                    ctArray.append(dsCT.pixel_array * dsCT.RescaleSlope + dsCT.RescaleIntercept)
                ctArray = np.array(ctArray[::-1])
            except:
                traceback.print_exc()
                pdb.set_trace()
            
            if len(doseArray) and len(ctArray):
                print ('')
                print (f' - doseArray  : {doseArray.shape} | ctArray: {ctArray.shape} ') # [Slices,H,W]
                print (f' ---- Dose Max: {np.max(doseArray)} | Dose Min: {np.min(doseArray)} ')
                print (f' - doseSpacing: {doseSpacing} | ctSpacing: {ctSpacing} ') # [H,W,Slices]
                print (f' - doseImagePositionPatient: {doseImagePositionPatient} ')
                print (f' - ctImagePositionPatientMin: {ctImagePositionPatientMin}, ctImagePositionPatientMax: {ctImagePositionPatientMax}')
                print ('')
                
                # Step 2 - Resample dose to CT
                doseImage = sitk.GetImageFromArray(doseArray)
                doseImage.SetSpacing(doseSpacing)
                resampler = sitk.ResampleImageFilter()
                resampler.SetOutputSpacing(ctSpacing)
                resampler.SetSize(ctArray.shape[::-1]) # SimpleITK convention: [H,W,Slices], numpy convention: [Slices,H,W]
                resampled_image = resampler.Execute(doseImage)
                doseArrayResampled = sitk.GetArrayFromImage(resampled_image)
                print (f' - doseArrayResampled  : {doseArrayResampled.shape} | ctArray: {ctArray.shape} ')

                dx, dy, dz = ((np.array(doseImagePositionPatient) - np.array(ctImagePositionPatientMax)) / np.array(ctSpacing)).astype(int)
                print (f' - dx, dy, dz: {dx, dy, dz} ')
                from scipy.ndimage.interpolation import shift
                doseArrayResampled = shift(doseArrayResampled, (dz, dy, dx))
                
                # Step 3 - Plot
                import matplotlib.pyplot as plt
                import matplotlib.colors as mcolors
                boundaries = [0      , 13.56       , 27.12        , 37.97        , 43.40      , 48.82    , 51.54     , 54.25   , 58.05        , 66.50        , 70.00      , 74.90]
                rgbColors = np.array([[0,0,0,0], [192,192,192,128], [192,192,255,128], [128,128,255,128], [64,64,192,128], [0,0,128,128], [0,255,0,128], [0,128,0,128], [255,255,255,128], [255,255,128,128], [255,192,0,128], [255,0,0,128]])/255.0
                cmap = mcolors.ListedColormap(rgbColors)
                norm = mcolors.BoundaryNorm(boundaries, cmap.N, clip=True)

                sliceDoseA = doseArrayResampled[sliceId-1]
                sliceCTA   = ctArray[sliceId-1]
                sliceDoseB = doseArrayResampled[sliceId]
                sliceCTB   = ctArray[sliceId]
                sliceDoseC = doseArrayResampled[sliceId+1]
                sliceCTC   = ctArray[sliceId+1]

                X = np.linspace(0, sliceDoseA.shape[1]-1, sliceDoseA.shape[1])
                Y = np.linspace(0, sliceDoseA.shape[0]-1, sliceDoseA.shape[0])
                X, Y = np.meshgrid(X, Y)
                
                f,axarr = plt.subplots(1,3)
                axarr[0].imshow(sliceCTA, cmap='gray', interpolation='bicubic')
                axarr[0].contourf(X,Y,sliceDoseA, levels=boundaries, colors=rgbColors)
                axarr[0].set_title(f'Slice={sliceId}')
                
                axarr[1].imshow(sliceCTB, cmap='gray', interpolation='bicubic')
                axarr[1].contourf(X,Y,sliceDoseB, levels=boundaries, colors=rgbColors)
                axarr[1].set_title(f'Slice={sliceId+1}')

                axarr[2].imshow(sliceCTC, cmap='gray', interpolation='bicubic')
                doseIm = axarr[2].contourf(X,Y,sliceDoseC, levels=boundaries, colors=rgbColors)
                axarr[2].set_title(f'Slice={sliceId+2}')
                
                from mpl_toolkits.axes_grid1 import make_axes_locatable
                divider = make_axes_locatable(axarr[2])
                cax = divider.append_axes("right", size="5%", pad=0.05)
                f.colorbar(doseIm, cax=cax, boundaries=boundaries, norm=norm, ticks=boundaries)

                plt.show(block=False)
                pdb.set_trace()

        else:
            print (f' - ERROR: pathCT: {Path(pathCT).exists()} | pathDoseA5: {Path(pathDoseA5).exists()} ')
            pdb.set_trace()
    
    except:
        traceback.print_exc()
        pdb.set_trace()