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.

Beamlet free optimization

This example will present the basis of beamlet optimization with the OpenTPS core.

#imports 

import numpy as np

from matplotlib import pyplot as plt

from opentps.core.data.images import CTImage
from opentps.core.data.images import ROIMask
from opentps.core.data.plan import ObjectivesList
from opentps.core.data.plan import PlanDesign
from opentps.core.data import DVH
from opentps.core.data import Patient
from opentps.core.data.plan import FidObjective
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

02/08/2023 02:04:53 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 02:04:53 PM - opentps.core._loggingConfig - INFO - Log level set: INFO
02/08/2023 02:04:54 PM - opentps.core.processing.imageProcessing.cupyImageProcessing - WARNING - Cannot import Cupy module
02/08/2023 02:04:55 PM - opentps.core.processing.registration.registrationMorphons - WARNING - cupy not found.
02/08/2023 02:04:55 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

#calibratioin of the CT
ctCalibration = readScanner(DoseCalculationConfig().scannerFolder)
bdl = mcsquareIO.readBDL(DoseCalculationConfig().bdlFile)

#creation of the patient object
patient = Patient()
patient.name = 'Patient'

#size of the 3D box
ctSize = 150

#creation of the CTImage object
ct = CTImage()
ct.name = 'CT'
ct.patient = patient

huAir = -1024. #Hounsfield unit of water
huWater = ctCalibration.convertRSP2HU(1.) #convert a stopping power of 1. to HU units
data = huAir * np.ones((ctSize, ctSize, ctSize))
data[:, 50:, :] = huWater
ct.imageArray = data #the CT generic image is created

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[100:120, 100:120, 100:120] = True
roi.imageArray = data

image = plt.imshow(ct.imageArray[110,:,:],cmap='Blues')
plt.colorbar(image)
plt.contour(roi.imageArray[110,:,:],colors="red")
plt.title("Created CT with ROI")
plt.text(5,40,"Air",color= 'black')
plt.text(5,100,"Water",color = 'white')
plt.text(106,111,"TV",color ='red')
plt.savefig('img/beamFree1.png',format = 'png')
plt.close()

png

Configuration of MCsquare

To configure the MCsquare calculator we need to calibrate it with the CT calibration obtained above

mc2 = MCsquareDoseCalculator()
mc2.beamModel = bdl
mc2.ctCalibration = ctCalibration
mc2.nbPrimaries = 1e7

Plan creation and design

We will now create a plan and set objectives for the optimization and set a goal of 20 Gy to the target

# Design plan
beamNames = ["Beam1"]
gantryAngles = [0.]
couchAngles = [0.]

# Generate new plan
planDesign = PlanDesign() #create a new plan
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"


plan.planDesign.objectives = ObjectivesList() #create a new objective set
plan.planDesign.objectives.setTarget(roi.name, 20.0) #setting a target of 20 Gy for the target
plan.planDesign.objectives.fidObjList = []
plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMAX, 20.0, 1.0)
plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMIN, 20.0, 1.0)
plan.planDesign.defineTargetMaskAndPrescription()

02/08/2023 02:05:10 PM - opentps.core.data.plan._planDesign - INFO - Building plan ...
02/08/2023 02:05:11 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 02:05:11 PM - opentps.core.processing.imageProcessing.roiMasksProcessing - INFO - Using SITK to dilate mask.
02/08/2023 02:05:16 PM - opentps.core.data.plan._planDesign - INFO - New plan created in 5.482358455657959 sec
02/08/2023 02:05:16 PM - opentps.core.data.plan._planDesign - INFO - Number of spots: 317

MCsquare beamlet free planOptimization

Now that we have every needed object defined we can compute the optimization through MCsquare. /!\ It may take some time to compute

doseImage = mc2.optimizeBeamletFree(ct, plan, [roi])

02/08/2023 02:05:16 PM - opentps.core.io.mhdIO - INFO - Write MHD file: C:\Users\romai\openTPS_workspace\Simulations\MCsquare_simulation\CT.mhd
02/08/2023 02:05:16 PM - opentps.core.io.mcsquareIO - INFO - Write plan: C:\Users\romai\openTPS_workspace\Simulations\MCsquare_simulation\PlanPencil.txt
02/08/2023 02:05:17 PM - opentps.core.io.mcsquareIO - INFO - Write plan objectives: C:\Users\romai\openTPS_workspace\Simulations\MCsquare_simulation\PlanObjectives.txt
02/08/2023 02:05:17 PM - opentps.core.io.mhdIO - INFO - Write MHD file: C:\Users\romai\openTPS_workspace\Simulations\MCsquare_simulation\structs\TV.mhd
02/08/2023 02:05:17 PM - opentps.core.processing.doseCalculation.mcsquareDoseCalculator - INFO - Start MCsquare simulation

Dose volume histogram

target_DVH = DVH(roi, doseImage)
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 = 17.891438802083336 Gy
D5 = 21.97541267641129 Gy
D5 - D95 =  4.083973874327953 Gy

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)
img_dose = resampleImage3DOnImage3D(doseImage, ct)
img_dose = img_dose.imageArray[:, :, Z_coord].transpose(1, 0)

02/08/2023 02:09:37 PM - opentps.core.processing.imageProcessing.roiMasksProcessing - INFO - Using SITK to dilate mask.

Plot of the dose

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 (%)")
ax[0].set_title("Computed dose")
ax[1].set_title("DVH")
plt.grid(True)
plt.legend()
plt.savefig('img/beamFree2.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 = 17.891438802083336 Gy
D5 = 21.97541267641129 Gy
D5 - D95 =  4.083973874327953 Gy

png

Download this notebook via our GitLab repository