All the examples might not be up to date on the website, please go to our GitLab repository and download latest version of the notebook.
Simple dose computation
In this example we are going to create a generic CT and use the MCsquare dose calculator to compute the dose image
#imports
import numpy as np
from matplotlib import pyplot as plt
import math
from opentps.core.data.images import CTImage
from opentps.core.data.images import ROIMask
from opentps.core.data.plan import PlanDesign
from opentps.core.data import DVH
from opentps.core.data import Patient
from opentps.core.io import mcsquareIO
from opentps.core.io.scannerReader import readScanner
from opentps.core.processing.doseCalculation.doseCalculationConfig import DoseCalculationConfig
from opentps.core.processing.doseCalculation.mcsquareDoseCalculator import MCsquareDoseCalculator
from opentps.core.processing.imageProcessing.resampler3D import resampleImage3DOnImage3D,resampleImage3D
02/08/2023 04:47:30 PM - root - INFO - Loading logging configuration: C:\Users\romai\AppData\Roaming\Python\Python39\site-packages\opentps\core\config\logger\logging_config.json
02/08/2023 04:47:30 PM - opentps.core._loggingConfig - INFO - Log level set: INFO
02/08/2023 04:47:31 PM - opentps.core.processing.imageProcessing.cupyImageProcessing - WARNING - Cannot import Cupy module
02/08/2023 04:47:32 PM - opentps.core.processing.registration.registrationMorphons - WARNING - cupy not found.
02/08/2023 04:47:32 PM - opentps.core.processing.C_libraries.libInterp3_wrapper - WARNING - cupy not found.
Generic CT creation
We will first create a generic CT of a box filled with water and air
ctCalibration = readScanner(DoseCalculationConfig().scannerFolder)
bdl = mcsquareIO.readBDL(DoseCalculationConfig().bdlFile)
patient = Patient()
patient.name = 'Patient'
ctSize = 150
ct = CTImage()
ct.name = 'CT'
ct.patient = patient
huAir = -1024.
huWater = ctCalibration.convertRSP2HU(1.)
data = huAir * np.ones((ctSize, ctSize, ctSize))
data[:, 50:, :] = huWater
ct.imageArray = data
Region of interest
We will now create a region of interest wich is a small 3D box of size 20*20*20
roi = ROIMask()
roi.patient = patient
roi.name = 'TV'
roi.color = (255, 0, 0) # red
data = np.zeros((ctSize, ctSize, ctSize)).astype(bool)
data[65:85, 65:85, 65:85] = True
roi.imageArray = data
Configuration of MCsquare
To configure the MCsquare calculator we need to calibrate it with the CT calibration obtained above
# Configure MCsquare
mc2 = MCsquareDoseCalculator()
mc2.beamModel = bdl
mc2.ctCalibration = ctCalibration
mc2.nbPrimaries = 1e7
Plan creation
We will now create a plan and create one beam
# Design plan
beamNames = ["Beam1","Beam2","Beam3"]
gantryAngles = [0.,90.,270.]
couchAngles = [0.,0.,0.]
# Generate new plan
planDesign = PlanDesign()
planDesign.ct = ct
planDesign.targetMask = roi
planDesign.gantryAngles = gantryAngles
planDesign.beamNames = beamNames
planDesign.couchAngles = couchAngles
planDesign.calibration = ctCalibration
planDesign.spotSpacing = 5.0
planDesign.layerSpacing = 5.0
planDesign.targetMargin = 5.0
plan = planDesign.buildPlan() # Spot placement
plan.PlanName = "NewPlan"
02/08/2023 04:47:49 PM - opentps.core.data.plan._planDesign - INFO - Building plan ...
02/08/2023 04:47:49 PM - opentps.core.processing.planOptimization.planInitializer - INFO - Target is dilated using a margin of 5.0 mm. This process might take some time.
02/08/2023 04:47:49 PM - opentps.core.processing.imageProcessing.roiMasksProcessing - INFO - Using SITK to dilate mask.
02/08/2023 04:47:58 PM - opentps.core.data.plan._planDesign - INFO - New plan created in 9.920757055282593 sec
02/08/2023 04:47:58 PM - opentps.core.data.plan._planDesign - INFO - Number of spots: 913
Center of mass
Here we look at the part of the 3D CT image where “stuff is happening” by getting the CoM. We use the function resampleImage3DOnImage3D to the same array size for both images.
roi = resampleImage3DOnImage3D(roi, ct)
COM_coord = roi.centerOfMass
COM_index = roi.getVoxelIndexFromPosition(COM_coord)
Z_coord = COM_index[2]
img_ct = ct.imageArray[:, :, Z_coord].transpose(1, 0)
contourTargetMask = roi.getBinaryContourMask()
img_mask = contourTargetMask.imageArray[:, :, Z_coord].transpose(1, 0)
image = plt.imshow(img_ct,cmap='Blues')
plt.colorbar(image)
plt.contour(img_mask,colors="red")
plt.title("Created CT with ROI")
plt.text(5,40,"Air",color= 'black')
plt.text(5,100,"Water",color = 'white')
plt.text(71,77,"TV",color ='red')
plt.savefig('img/SimpleDose1.png',format = 'png')
plt.close()
02/08/2023 04:47:59 PM - opentps.core.processing.imageProcessing.roiMasksProcessing - INFO - Using SITK to dilate mask.
Dose computation
We now use the MCsquared dose calculator to compute the dose of the created plan
doseImage = mc2.computeDose(ct, plan)
02/08/2023 04:47:59 PM - opentps.core.processing.doseCalculation.mcsquareDoseCalculator - INFO - Prepare MCsquare Dose calculation
02/08/2023 04:47:59 PM - opentps.core.io.mhdIO - INFO - Write MHD file: C:\Users\romai\openTPS_workspace\Simulations\MCsquare_simulation\CT.mhd
02/08/2023 04:47:59 PM - opentps.core.io.mcsquareIO - INFO - Write plan: C:\Users\romai\openTPS_workspace\Simulations\MCsquare_simulation\PlanPencil.txt
02/08/2023 04:48:01 PM - opentps.core.processing.doseCalculation.mcsquareDoseCalculator - INFO - Start MCsquare simulation
img_dose = resampleImage3DOnImage3D(doseImage, ct)
img_dose = img_dose.imageArray[:, :, Z_coord].transpose(1, 0)
scoringSpacing = [2, 2, 2]
scoringGridSize = [int(math.floor(i / j * k)) for i, j, k in zip([150,150,150], scoringSpacing, [1,1,1])]
roiResampled = resampleImage3D(roi, origin=ct.origin, gridSize=scoringGridSize, spacing=scoringSpacing)
target_DVH = DVH(roiResampled, doseImage)
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
#ax[0].axes.get_xaxis().set_visible(False)
#ax[0].axes.get_yaxis().set_visible(False)
ax[0].imshow(img_ct, cmap='gray')
ax[0].imshow(img_mask, alpha=.2, cmap='binary') # PTV
dose = ax[0].imshow(img_dose, cmap='jet', alpha=.2)
plt.colorbar(dose, ax=ax[0])
ax[1].plot(target_DVH.histogram[0], target_DVH.histogram[1], label=target_DVH.name)
ax[1].set_xlabel("Dose (Gy)")
ax[1].set_ylabel("Volume (%)")
plt.grid(True)
plt.legend()
plt.savefig('img/SimpleDose2.png',format = 'png')
plt.close()
print('D95 = ' + str(target_DVH.D95) + ' Gy')
print('D5 = ' + str(target_DVH.D5) + ' Gy')
print('D5 - D95 = {} Gy'.format(target_DVH.D5 - target_DVH.D95))
D95 = 10.667724609375 Gy
D5 = 15.92838062959559 Gy
D5 - D95 = 5.2606560202205905 Gy
Download this notebook via our GitLab repository