See this page for information about how to run CLI module from a Python script: https://www.slicer.org/wiki/Documentation/Nightly/Developers/Python_scripting#Running_a_CLI_from_Python
1 Like
It work now, thank you
1 Like
Posting here a script I made using pydicom and SimpleITK by inputting CT and RTDOSE .dcms. Sample Output
"""
Step 1 - Read Dose and CT .dcm files
Step 2 - Resample dose to CT
Step 2.1 - Shift dose to match 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} ')
# Step 2.1 - Shift dose to match CT
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()