#******************************************************************************** # Title: autoPlanBrainArc.py # Author: John Eley # Date: 7 November 2014 # # Purpose: Automatically generate beams for arc plan and start optimization for single-field uniform dose # # Input: See required input below in USER INPUT # # Notes: Edited 12/09/14 for use for the proton arc therapy paper # Target doses are set to be robust, but healthy brain is not, that would likely be a conflicting objective # # Isocenter must be specified manually, and some permutation is needed and tricky (hint, record script adding new beam to get them) # # New RayStation Version slightly changed name of plan->beamset->TreatmentTechnique to be "ProtonPencilBeamScanning" instead of "PencilBeamScanning" # # Plan Name Prefix Nomenclature # Plan type: "UD" or "IM" (uniform dose or intensity modulated (pair)) # Robust: "r0" = no, "r1" = yes # Grid size: "g0.5" = 0.5 cm dose grid # Example: "UD_r1_g0.5_" (do end with trailing underscore) # # 060315 - Updated for RayStation 4.6 (means it will no longer be compatible with 4.5, changes to ... # TreatmentTechnique="ProtonPencilBeamScanning" (previously was "PencilBeamScanning") # In the "CreatePbsIonBeam" function, the treatment technique line is deleted now # In "RunOpimization" the reduceOARdose thing is no longer needed #******************************************************************************** # # USER INPUT # # Run specific variables (unique for each run) pcPlanNamePrefix = 'IM_r1_n13_' iFlagIMPTArcPair = 1 # 0=single beam per fraction, unif dose, 1=IMPTpair per fraction, together provide uniform dose but allows IMPT bw pair # Variables common to study dPrescribedDose = 7800 # cGy target dose dThetaNaught = 0 # deg, gantry angle of first (existing) beam, which should be named "f1" dArcAngle = 360 # deg, extent of entire arc (e.g., 180 deg or 360 deg) nPlans = 13 # n, this is also nBeams (for the single beam arc, IMPT pairs have nBeams = 2 * nPlans), use 16 for AAPM abstract 030915, use 4 and 13 for paper (4 gives MGH style, 13 is the pseudoarc) dDeltaThetaPair = 45 # deg, angular spacing bw beams in IMPT pairs dDoseGridSize = 0.5 # cm, voxel width (must be within 0.1 - 0.5 cm, 0.5 cm for testing, 0.2 mm for final paper), cannot run less than 0.4 with robustness and 16 plans and you get crash pcIsRobust = "True" # "True" = use robust opt (2-mm xyz, 3.5% range), "False" = no robust opt (CTV plan only), greatly affects speed targetName = "CTV" iFlagAvoidanceVolume = 0 # if 1, add an avoidance constraint to objective function (dmax = 0, w = 1), this is okay even for UD b/c I have a hard constraint on target coverage avoidanceVolumeName = "" # only used if iFlagAvoidanceVolume == 1 brainName = "External" # This is used only to encourage doses outside the target to be low, it can be brain or external dXIso = 1.15 # Must be determined for each patient and target, hint record script adding a beam to find this dYIso = 11.32 dZIso = 16.33 pcRangeShifter = "rs57" # # END USER INPUT # # Occasionally changed variables nMaxNumberOfIterations = 200 # I think 200 is safe, 70-80 seems typical for convergence pcBeamName = "f1" # name of first beam in plan pcBeamName2 = "f2" # for IMPT arc pair dPlanPrescriptionDose = float(dPrescribedDose) / float(nPlans) # Get libraries and link to RayStation data import time import math import string from connect import * patient = get_current("Patient") #plan = get_current("Plan") #beam_set = get_current("BeamSet") startTime = time.time() # Calculate the interbeam spacing if dArcAngle == 360: dDeltaTheta = float(dArcAngle) / float(nPlans) # This if statment is needed so we don't place 2 beams at 0 deg (since 0 and 360 would each get one) else: dDeltaTheta = float(dArcAngle) / (float(nPlans) - 1) # This makes the first beam at the start and the last beam at the end of the arc # # Add plans for arc # for iPlan in range(1,nPlans + 1): print(" ") print("Plan: " + str(iPlan)) print(" ") pcPlanName = pcPlanNamePrefix + str(iPlan) dThetaPrime = float(dThetaNaught) + (float(iPlan) - 1) * float(dDeltaTheta) # Handle wrapping of gantry angle if dThetaPrime >= 360.0: dTheta = float(dThetaPrime) - 360.0 else: dTheta = float(dThetaPrime) with CompositeAction('Add Treatment plan'): plan = patient.AddNewPlan(PlanName=pcPlanName, PlannedBy="", Comment="", ExaminationName="CT 1", AllowDuplicateNames=False) beamset = plan.AddNewBeamSet(Name=pcPlanName, ExaminationName="CT 1", MachineName="ProBeam", NominalEnergy=0, Modality="Protons", TreatmentTechnique="ProtonPencilBeamScanning", PatientPosition="HeadFirstSupine", NumberOfFractions=1, CreateSetupBeams=False, UseLocalizationPointAsSetupIsocenter=False, Comment="") beamset.AddDosePrescriptionToRoi(RoiName=targetName, DoseVolume=0, PrescriptionType="AverageDose", DoseValue=dPlanPrescriptionDose, RelativePrescriptionLevel=1, AutoScaleDose=False) with CompositeAction('Add beam (' + pcBeamName + ', Beam Set: ' + pcPlanName + ')'): beamset.CreatePbsIonBeam(SnoutId="S1", SpotTuneId="4.0", RangeShifter=pcRangeShifter, MinimumAirGap=None, Isocenter={ 'x': dXIso, 'y': dYIso, 'z': dZIso }, Name=pcBeamName, Description="", GantryAngle=dTheta, CouchAngle=0, CollimatorAngle=0, ApertureBlock=None) with CompositeAction('Set default grid'): plan.SetDefaultDoseGrid(VoxelSize={ 'x': dDoseGridSize, 'y': dDoseGridSize, 'z': dDoseGridSize }) plan.ComputeRoiVoxelVolumes() # If IMPT arc pair is flagged, add another beam orthogonal for IMPT if iFlagIMPTArcPair == 1: dTheta2Prime = float(dTheta) + dDeltaThetaPair if dTheta2Prime >= 360.0: dTheta2 = float(dTheta2Prime) - 360.0 else: dTheta2 = float(dTheta2Prime) beamset.CopyBeam(BeamName=pcBeamName) beamset.Beams[1].Name = pcBeamName2 beamset.Beams[1].GantryAngle = dTheta2 # Add objective functions dMinDose = float(dPlanPrescriptionDose) * 0.95 dMinDose93 = float(dPlanPrescriptionDose) * 0.93 dMaxDose = float(dPlanPrescriptionDose) * 1.05 dMaxDose107 = float(dPlanPrescriptionDose) * 1.07 dMaxDose115 = float(dPlanPrescriptionDose) * 1.15 dMaxBrainDose = float(dPlanPrescriptionDose) * 0.95 with CompositeAction('Add Optimization Function'): plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="UniformDose", RoiName=targetName, IsConstraint=False, RestrictAllBeamsIndividually=False, RestrictToBeam=None, IsRobust=False, RestrictToBeamSet=None) plan.PlanOptimizations[0].Objective.ConstituentFunctions[0].DoseFunctionParameters.DoseLevel = dPlanPrescriptionDose plan.PlanOptimizations[0].Objective.ConstituentFunctions[0].DoseFunctionParameters.Weight = 1 with CompositeAction('Add Optimization Function'): plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="MinDose", RoiName=targetName, IsConstraint=True, RestrictAllBeamsIndividually=False, RestrictToBeam=None, IsRobust=pcIsRobust, RestrictToBeamSet=None) plan.PlanOptimizations[0].Constraints[0].DoseFunctionParameters.DoseLevel = dMinDose93 # This forces the absolute min dose to be >= 93% of Rx with CompositeAction('Add Optimization Function'): plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="MinDvh", RoiName=targetName, IsConstraint=True, RestrictAllBeamsIndividually=False, RestrictToBeam=None, IsRobust=pcIsRobust, RestrictToBeamSet=None) plan.PlanOptimizations[0].Constraints[1].DoseFunctionParameters.DoseLevel = dMinDose # The makes sure 99% of target is covered by >= 95% of Rx plan.PlanOptimizations[0].Constraints[1].DoseFunctionParameters.PercentVolume = 99 # ELEY 082815 Changed this hard constraint from 107% to 115% because chordoma1 could not meet these constraints and the plan did poorly in optimization with CompositeAction('Add Optimization Function'): plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="MaxDose", RoiName=targetName, IsConstraint=True, RestrictAllBeamsIndividually=False, RestrictToBeam=None, IsRobust=pcIsRobust, RestrictToBeamSet=None) plan.PlanOptimizations[0].Constraints[2].DoseFunctionParameters.DoseLevel = dMaxDose107 # ELEY 082815 Changed this one to be an objective rather than a constraint with CompositeAction('Add Optimization Function'): plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="MaxDvh", RoiName=targetName, IsConstraint=False, RestrictAllBeamsIndividually=False, RestrictToBeam=None, IsRobust=False, RestrictToBeamSet=None) #plan.PlanOptimizations[0].Constraints[3].DoseFunctionParameters.DoseLevel = dMaxDose #plan.PlanOptimizations[0].Constraints[3].DoseFunctionParameters.PercentVolume = 1 plan.PlanOptimizations[0].Objective.ConstituentFunctions[1].DoseFunctionParameters.DoseLevel = dMaxDose plan.PlanOptimizations[0].Objective.ConstituentFunctions[1].DoseFunctionParameters.PercentVolume = 1 plan.PlanOptimizations[0].Objective.ConstituentFunctions[1].DoseFunctionParameters.Weight = 1 # This seems to cause a crash for whatever reason # Eley 082815 added this, if no dose over 107 is allowed in the target its not allowed outside either #with CompositeAction('Add Optimization Function'): # plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="MaxDose", # RoiName=brainName, # IsConstraint=True, # RestrictAllBeamsIndividually=False, # RestrictToBeam=None, # IsRobust=pcIsRobust, # RestrictToBeamSet=None) # plan.PlanOptimizations[0].Constraints[3].DoseFunctionParameters.DoseLevel = dMaxDose107 with CompositeAction('Add Optimization Function'): plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="MaxDose", RoiName=brainName, IsConstraint=False, RestrictAllBeamsIndividually=False, RestrictToBeam=None, IsRobust=False, RestrictToBeamSet=None) plan.PlanOptimizations[0].Objective.ConstituentFunctions[2].DoseFunctionParameters.DoseLevel = dMaxBrainDose plan.PlanOptimizations[0].Objective.ConstituentFunctions[2].DoseFunctionParameters.Weight = 100 with CompositeAction('Add Optimization Function'): plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="MaxDose", RoiName=brainName, IsConstraint=False, RestrictAllBeamsIndividually=False, RestrictToBeam=None, IsRobust=False, RestrictToBeamSet=None) plan.PlanOptimizations[0].Objective.ConstituentFunctions[3].DoseFunctionParameters.DoseLevel = 0 plan.PlanOptimizations[0].Objective.ConstituentFunctions[3].DoseFunctionParameters.Weight = 0.01 if iFlagAvoidanceVolume == 1: with CompositeAction('Add Optimization Function'): plan.PlanOptimizations[0].AddOptimizationFunction(FunctionType="MaxDose", RoiName=avoidanceVolumeName, IsConstraint=False, RestrictAllBeamsIndividually=False, RestrictToBeam=None, IsRobust=False, RestrictToBeamSet=None) plan.PlanOptimizations[0].Objective.ConstituentFunctions[4].DoseFunctionParameters.DoseLevel = 0 plan.PlanOptimizations[0].Objective.ConstituentFunctions[4].DoseFunctionParameters.Weight = 1 # Set robust opt parameters (2-mm positional uncertainty w 3.5% range uncertainty) plan.PlanOptimizations[0].OptimizationParameters.SaveRobustnessParameters(PositionUncertaintyAnterior=0.2, PositionUncertaintyPosterior=0.2, PositionUncertaintySuperior=0.2, PositionUncertaintyInferior=0.2, PositionUncertaintyLeft=0.2, PositionUncertaintyRight=0.2, DensityUncertainty=0.035, IndependentBeams=False, ComputeExactScenarioDoses=False) # Start optimization plan.PlanOptimizations[0].OptimizationParameters.Algorithm.MaxNumberOfIterations = nMaxNumberOfIterations plan.PlanOptimizations[0].RunOptimization() # # END for iPlan # print("Execution time: " + str(time.time() - startTime) + " seconds")