#******************************************************************************** # Title: interplaySpots.py # Author: John Eley # Date: 24 February 2015 # # Purpose: Simulate interplay for single fraction # # Notes: # # I create nMotionState copies of the original plan, calculate the state dose for each, warp back to the ref phase, then I sum the # total 4D dose # # Before running script plan is already created and optimized for a single # reference plan and deformation fields are # already calculated for every motion state wrt the reference motion state. # # I always index registration 0, this must be upgraded if the user wants to have more than one registration for a patient # # Currently I assume a fixed dose rate of 16.5 MU / min = 0.275 MU/s = 0.000275 MU / ms , # which I determined by optimizing a plan using the ProBeam model in RS for "zzz water phantom" 2Gy to 10x10x10cm3 target required 16.5 MU, and I know this should take 60s on our system # # Note to self, Ok, here is what I found (1/28/14), if you sum up all the spot weights in the entire beam # over all IESs (I only test single beam so far), you get 1. So # this tells me that also the relative weight of the IES slice (Segments.RelativeWeight) is equal to the sum # of the spot weights in that IES, so this also means that # (spotMU = spotWeight * beamMU), nice and simple # # The sum of Beams[0].Segments[iIES].RelativeWeights over all iIES is 1 for the input reference plan (I guess for each beam but didn't test multi beams) # # I made dTInterbeamPause constant as 4 minutes (240000 ms), which is really just a guess at the time between beams for a patient # # Currently, only supports iFraction = 0 (1st fraction) # # I am not 10)% sure if my registrations transform the dose the "right way" maybe vectors need to be reversed, should check state doses to confirm # # It is assumed that all CT states listed in ppcMotionStates are (1) in order temporally, (2) equal phase width, i.e., equ-phase-binned, and (3) cover 2*pi, non-repetititive (i.e., first does not equal last) # # Warning: I do not guarantee state dose or 4D dose if you recalculate using Raystation buttons, as I overwrite dose values here. Recalculating # any of the state or 4D plans manually in RayStation may invalidate 4D dose or state dose. # # Well after looking into it again it seems (4/6/15) that the reason I gave up on eval doses was because # I cannot sum them automatically. Therefore, maybe I need to stick with the multiplan idea # but maybe I can give the user a flag to auto delete state plans if so desired. # # Update: now seems to work when planned on average CT too # # Limitations are still only 1st fraction supported (no summing of fractions yet) # and head first supine only (search for "supine" to see dependencies) # # To be compatible with v4.5 or v4.6, TreatmentTechnique should be either "PencilBeamScanning" or "ProtonPencilBeamScanning", respectively #******************************************************************************** # USER INPUT outputPlanNamePrefix = '' # prefix/nickname for 4D dose trials (e.g., "t1" for test 1) refPlanName = 'eleyNB' # This should be the reference 3D plan (3D static as Christoph would call it) on the reference state CT, can either be a plan on a single phase of 4DCT or even a plan on average CT or freebreathing CT ppcMotionStates = ['1', '2', '3', '4', '5', '6', '7', '8', '9', '10'] # List the numbers of all state CTs (phase CTs) in your preferred order (and covering 2 pi cycle) including the reference phase, e.g., include "4" for "CT 4", this is not the same as the index below iMotionState, iMotionState goes from 0 to nMotionStates - 1, ppcMotionStates refer to the CT names and can be in any order or even skip numbers, they are just names pcRefMotionState = '2' # number matching the name of the reference phase CT, corresponds to ppcMotionStates, 4D dose will be displayed here, if planning is done on a 4DCT phase, it is assumed to be this one dStartingPhase = 0 # deg (out of 360), phase of respiration when irradiation starts dRespiratoryPeriod = 4000 # ms, respiratory period iFlagRefPlanOutside4dct = 1 # 0 = assumes that the reference plan was made on the reference state (of the 4DCT) above, 1 = assumes the plan was made outside the 4DCT (e.g., on an average CT, MIP, etc.) iFlagDoseRateThrottling = 0 # 0 = assume constant dose rate for all IES (16.5 MU / min), 1 = throttle dose rate individually for each IES so that minimum spot MU takes 3 ms iFlagDebugWStateTree = 0 iDefRegistration = 0 # if there is only 1 deformable registration, this should be set to 0, otherwise, one must know which to retrieve # END USER INPUT # Load libraries and variables from RayStation from connect import * import copy import math import string # string library used for string manipulation import json # useful library to manipulate strings to lists and back and forth patient = get_current("Patient") plan = patient.TreatmentPlans[refPlanName] beamSet = plan.BeamSets[refPlanName] #plan = get_current("Plan") #beamSet = get_current("BeamSet") #refPlanName = plan.Name # This works, but requires the reference plan be selected as the active plan, I prefer to just manually specify this above in USER INPUT # Debugging if iFlagDebugWStateTree == 1: import statetree statetree.RunStateTree() # # ## ### Determine the time of irradiation for each spot in each IES in each beam ## # # print(" ") print("- ") print("-- ") print("--- Calculating Irradiation Times for All Spots") print("-- ") print("- ") print(" ") # Determine number of beams (Always when finding n I use complex string counting arguments bc RS doesn't just store a simple variable I can access like nBeams ppcTempBeamNames = beamSet.Beams._ # This makes a list of the beam names nBeams = len(ppcTempBeamNames.split()) # This counts the number of separate beam names print("Testing nBeams = " + str(nBeams)) # Initialize a nested list to hold the number of spots in each beam, each IES piNIES = [0] * nBeams # n IES in each beam, access via: piNIES[iBeam] ppiNSpots = [0] * nBeams # n spots in this IES in this beam, access via: ppiNSpots[iBeam][iIES] # Initialize a nested list to hold the dose rate for each IES in each beam ppdDoseRate = [0] * nBeams # MU / ms, access via: ppdDoseRate[iBeam][iIES] # Initialize a nested list to hold the start times for each spot in each IES in each beam pppdTStartIrr = [0] * nBeams # ms, start times for each spot relative to beam-on (at t = 0), access via pppdTStartIrr[iBeam][iIES][iSpot] dTSpotStartIrr = 0 # ms, variable to keep track of current time that spot irradiation starts for iBeam in range(0,nBeams): # Determine the number of iso-energy slices (IES) in this beam, btw RayStation calls these "segments" instead of IES or energy layers pcTempIESIndices = beamSet.Beams[iBeam].Segments._ # This grabs a string containing the segment indices, but requires some manipulation to actually count the segments (cleaning up the string) pcIESIndices = pcTempIESIndices.translate(None, '_') # remove underscores from string piIESIndices = json.loads('[' + pcIESIndices + ']') # this line relies on the function json.loads that converts strings into lists/arrays of integers piNIES[iBeam] = len(piIESIndices) # Determine the number of spots in each IES for this beam ppiNSpots[iBeam] = [0] * piNIES[iBeam] # update / reinitialize for the correct the number of IES in each beam for iIES in range(0,piNIES[iBeam]): ppiNSpots[iBeam][iIES] = len(beamSet.Beams[iBeam].Segments[iIES].Spots.Weights) # # Specify beam temporal characteristics of Varian cyclotron (MU per second as a function of energy) and inter IES pause time, inter spot time # dTInterIESPause = 1000 # ms, pause bw energy layers of 1000 ms dTInterSpotPause = 0.05 # ms, pause bw spots dTInterBeamPause = 240000 # ms, pause bw beams (currently this is arbitrarily set to 4 minutes) # Determine "dose rate" (really MU rate) for each energy layer (IES) in this beam ppdDoseRate[iBeam] = [0] * piNIES[iBeam] # MU / ms, update / reinitialize for the correct the number of IES in each beam for iIES in range(0,piNIES[iBeam]): #dE = beamSet.Beams[iBeam].Segments[iIES].NominalEnergy if iFlagDoseRateThrottling == 1: # Looks for the lowest weighted spot in the IES and sets the dose rate such that that spot takes 3 ms (dose rate throttling), which is what our system is actually going to do dMinMU = min(beamSet.Beams[iBeam].Segments[iIES].Spots.Weights) if dMinMU == 0: print("iIES: " + str(iIES)) print " Minimum MU in this IES is zero, there are going to be problems with this. This method only provides reasonable results if an appropriate min MU is set in optimization. - JE" ppdDoseRate[iBeam][iIES] = dMinMU / 3 # MU / ms, dose rate is throttled such that minMU for this IES is delivered in 3 ms else: # Assume constant dose rate for all IES ppdDoseRate[iBeam][iIES] = 0.000275 # MU / ms (this rough approximation gives 2 Gy to 1 L volume in 1 minute, see notes above) # END for iIES # END specify beam temporal characteristics # # Build a list giving the irradiation start time for each spot in each IES # I assume irradiation order is the same as the order of IES and spots given by RayStation # which seems to be highest to lowest energy in that order # pppdTStartIrr[iBeam] = [0] * piNIES[iBeam] # ms, reinitialize for this beam for this nIES for iIES in range(0,piNIES[iBeam]): pppdTStartIrr[iBeam][iIES] = [0] * ppiNSpots[iBeam][iIES] # Calculate the irradiation time for each spot in this IES in this beam for iSpot in range(0,ppiNSpots[iBeam][iIES]): pppdTStartIrr[iBeam][iIES][iSpot] = dTSpotStartIrr; # ms dSpotMU = beamSet.Beams[iBeam].Segments[iIES].Spots.Weights[iSpot] * \ beamSet.Beams[iBeam].BeamMU # MU for this spot (remember Weights are relative to the whole beam and sum to 1, hence I multiple by BeamMU) dTSpotIrr = dSpotMU / ppdDoseRate[iBeam][iIES] # ms, time for this spot to complete irradiation dTSpotStartIrr = dTSpotStartIrr + dTSpotIrr + dTInterSpotPause # ms, update to main time variable (add elapsed time) # END for iSpot dTSpotStartIrr = dTSpotStartIrr + dTInterIESPause # add pause bw energy layers # END for iIES dTSpotStartIrr = dTSpotStartIrr + dTInterBeamPause # add pause bw beam delivery # END for iBeam # NOTE: I would like to add a function here to plot the spot MU and irradiation time x is time, y is MU) # # ## ### Determine the phase (motion state) of irradiation for each spot in each IES in each beam (see notes from 1/27/15) ## # # print(" ") print("- ") print("-- ") print("--- Sorting Beam Spots into Motion States (Respiratory Phases)") print("-- ") print("- ") print(" ") nMotionStates = len(ppcMotionStates) dDeltaT = dRespiratoryPeriod / nMotionStates # ms, time between adjacent motion states pppiMotionStateForSpot = [0] * nBeams # initialize for iBeam in range(0,nBeams): pppiMotionStateForSpot[iBeam] = [0] * piNIES[iBeam] # update / reinitialize for iIES in range(0,piNIES[iBeam]): pppiMotionStateForSpot[iBeam][iIES] = [0] * ppiNSpots[iBeam][iIES] # update / reinitialize for iSpot in range(0,ppiNSpots[iBeam][iIES]): dTShift = pppdTStartIrr[iBeam][iIES][iSpot] + \ (dStartingPhase * dRespiratoryPeriod / 360) # time that spot irradiation starts, shifted to allow simulation of variable patient respiratory phase at "beam on" nPeriodsElapsed = math.floor(dTShift / dRespiratoryPeriod) # n, number of respiratory periods completed since start of simulation (simulation always starts with patient at phase = 0, beam on can be shifted to later phases using the startingPhase variable dTPrime = dTShift - (nPeriodsElapsed * dRespiratoryPeriod) # subtract repetitive cycles iMotionStateForSpot = int(math.floor(dTPrime / dDeltaT)) # this is the motion state in which this spot will start (and I assume will also end) irradiation (I do not handle spots that cross motio states) pppiMotionStateForSpot[iBeam][iIES][iSpot] = iMotionStateForSpot # END for iSpot # END for iIES # End for iBeam # # ## ### Create state plans and calculate dose for each motion state, warp doses to reference motion state ## # # # Maybe this next line will become obsolete? dSumAllIESWeights = 1 # (I think this is always 1) This is the sum of Segments[iIES].RelativeWeight and is also the sum of all spot weights in a beam # Just keep in mind that iMotionState != pcMotionState, pcMotionState is the number given the CT phase, it is really a name not an index, and does not have to go from 0 to n-1 piDoseEvaluationIndex = [0] * nMotionStates # initialize iRefMotionState = -1 # initialize to -1 so it does not match any accidental CT piWarpDoseEvaluationIndex = [-1] * nMotionStates # initialize to -1 so it does not match any indices accidentally ppcStatePlanNames = [0] * nMotionStates # Initialize a list to hold state plan names piFlagMotionStateHasDose = [0] * nMotionStates # logic flags, if motion state has no dose, don't try to sum it later or it may cause probelms nMotionStatesWithDose = 0 for iMotionState in range(0,nMotionStates): print(" ") print("- ") print("-- ") print("--- Calculating Motion State: " + str(iMotionState)) print("-- ") print("- ") print(" ") pcMotionState = ppcMotionStates[iMotionState] if pcMotionState == pcRefMotionState: iRefMotionState = iMotionState #statePlanName = str(outputPlanNamePrefix) + '_ref_' + refPlanName + 't' + str(dRespiratoryPeriod) + 'p' + str(dStartingPhase) + '_ms' + str(pcMotionState) statePlanName = str(outputPlanNamePrefix) + refPlanName + '_p' + str(dStartingPhase) + '_ms' + str(pcMotionState) + '_x' statePlanNameFinal = str(outputPlanNamePrefix) + refPlanName + '_p' + str(dStartingPhase) + '_ms' + str(pcMotionState) #print("Testing statePlanName: " + statePlanName) ppcStatePlanNames[iMotionState] = statePlanNameFinal # name used later for summing all state plans # # Create a new plan for each motion state # patient.CopyPlan(PlanName=refPlanName, NewPlanName=statePlanName) # # Calculate the state dose on the refPlanName CT, # this is just a placeholder and for now is an identical copy of the reference plan, I will overwrite this later with the warped dose # but I want to put it in the "real" plan dose location not an evaluation dose # patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName].ComputeDose(ComputeBeamDoses=True, DoseAlgorithm="SpotWeightPencilBeam", ForceRecompute=True) # ELEY 040615, this line above might not actually be needed ???, just leave it in for now until better understood # # Loop over each beam, each IES, each spot - make a local copy of the spot weights, set spots to zero if they don't occur in this motion state, rescale beamMU and rel spot weight as needed, then tie back into main spot lists, # pdSumBeamSpotWeights = [0] * nBeams # initialize array to contain the sum of spot weights for each beam for this motion state (used later to scale BeamMU) for iBeam in range(0,nBeams): # Test if beam has any spots in this motion state, if not, just set the entire beam MU to essentially zero iBeamHasSpots = 0 for iIES in range(0,piNIES[iBeam]): for iSpot in range(0,ppiNSpots[iBeam][iIES]): if pppiMotionStateForSpot[iBeam][iIES][iSpot] == iMotionState: iBeamHasSpots = 1 piFlagMotionStateHasDose[iMotionState] = 1 # importantly, this is initialized outside of beams and even MS loops, I want to count if any beams give dose to this MS # END if # END for iSpot # End for iIES if iBeamHasSpots == 0: # Set beam MU to nearly zero (1/10000 original MU), RS doesn't allow actually setting to zero (and call it a day) patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName] \ .Beams[iBeam].BeamMU = patient \ .TreatmentPlans[statePlanName] \ .BeamSets[statePlanName] \ .Beams[iBeam].BeamMU / 10000 # You can change the beam MU but it cannot be zero or RayStation gives error so I just set it to 1/10000 of the initial weight print(" ") print(": State plan " + str(pcMotionState) + " has no dose contribution from beam " + str(iBeam) + ".") print(" ") else: # Loop over all spots, sum spot weights for this beam for this motion state for all IES for iIES in range(0,piNIES[iBeam]): # Make a local copy of the spot weights in this IES pdSpotWeights = 0 # Initialize / erase entire list (from previous iterations) pdSpotWeights = patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName] \ .Beams[iBeam].Segments[iIES] \ .Spots.Weights # Find sum of spot weights for this beam for all IES that do occur in this motion state for iSpot in range(0,ppiNSpots[iBeam][iIES]): if pppiMotionStateForSpot[iBeam][iIES][iSpot] == iMotionState: pdSumBeamSpotWeights[iBeam] = pdSumBeamSpotWeights[iBeam] + pdSpotWeights[iSpot] # Keep a running tab of the sum spot weights that remain in this motion state for each beam # END if # END for iSpot # END for iIES # Rescale beam spot weights and beam MU for iIES in range(0,piNIES[iBeam]): # Again, Make a local copy of the spot weights in this IES (For whatever reason, before I had problems to overwrite spot weights directly, but if I make a local copy I can change them, and then I can pass the whole local copy list back into the IES with, surprisingly no problems) pdSpotWeights = 0 # Initialize / erase entire list (from previous iterations) pdSpotWeights = patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName] \ .Beams[iBeam].Segments[iIES] \ .Spots.Weights for iSpot in range(0,ppiNSpots[iBeam][iIES]): if pppiMotionStateForSpot[iBeam][iIES][iSpot] == iMotionState: pdSpotWeights[iSpot] = pdSpotWeights[iSpot] / pdSumBeamSpotWeights[iBeam] # Scale up so that the sum of spot weights is always 1 for the beam else: pdSpotWeights[iSpot] = 0 # Zero weight to spots that do not fall in this motion state # END if / else # END for iSpot # Here pass back in the modified local copy of the spots to the IES patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName] \ .Beams[iBeam].Segments[iIES] \ .Spots.Weights = pdSpotWeights # END for iIES # And here finally reset beam MU (scale down since I scaled the spot weights up, see above) patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName] \ .Beams[iBeam].BeamMU = patient \ .TreatmentPlans[statePlanName] \ .BeamSets[statePlanName] \ .Beams[iBeam].BeamMU * pdSumBeamSpotWeights[iBeam] # END if iBeamHasSpots (in this motion state) # End for iBeam # Keep track of the number of states with dose, this is needed for the final 4D sum if piFlagMotionStateHasDose[iMotionState] == 1: nMotionStatesWithDose = nMotionStatesWithDose + 1 # # If this is not the reference state OR if this is the reference state, but the ref plan was created on another CT image (other than the ref state of the 4DCT) # then compute state dose (ComputeDoseOnAdditionalSets), furthermore if it is not the ref CT state, then warp it onto the ref CT state # # This else is entered only if (1) this iMotionState is not the reference motion state, or (2) the reference plan was not planned on the 4DCT ref state if pcMotionState != pcRefMotionState or iFlagRefPlanOutside4dct == 1: print("Testing ComputeDoseOnAdditionalSets for pcMotionState: " + str(pcMotionState)) print("Testing: If plan was made on average CT, you should see this execute for the ref CT state.") # # Calculate the state dose on the correct state CT # patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName] \ .ComputeDoseOnAdditionalSets(OnlyOneDosePerImageSet=False, AllowGridExpansion=False, ExaminationNames=["CT " + str(pcMotionState)], FractionNumbers=[0], ComputeBeamDoses=True) # Take care of some routing and indexing issues, I need to find an exam index where the current state CT name is stored ppcTempExaminationNames = patient.TreatmentDelivery.FractionEvaluations[0].DoseOnExaminations._ # This makes a list of the exam names nExams = len(ppcTempExaminationNames.split()) # This counts the number of separate exam names iExamToRetrieve = -1 for iExam in range(0,nExams): if (patient.TreatmentDelivery.FractionEvaluations[0].DoseOnExaminations[iExam].OnExamination.Name) == ('CT ' + str(pcMotionState)): iExamToRetrieve = iExam print("Motion State: " + str(iMotionState)) print("iExamToRetrieve: " + str(iExamToRetrieve)) # END if if (patient.TreatmentDelivery.FractionEvaluations[0].DoseOnExaminations[iExam].OnExamination.Name) == ('CT ' + str(pcRefMotionState)): iRefExam = iExam print("Motion State: " + str(iMotionState)) print("iRefExam: " + str(iRefExam)) # END if # END for # Also, I need to determine which dose evaluation index I will use, since I have just previously created a new dose evaluation (that is the state dose on the state CT), I think I am safe to always grab the highest index of a dose evaluation (e.g., if there were previously 4 dose evaluations, now there will be 5 and I want the 5th one) pcTempDoseEvaluationIndices = patient.TreatmentDelivery.FractionEvaluations[0].DoseOnExaminations[iExamToRetrieve].DoseEvaluations._ # This grabs a string containing the dose exam indices, but requires some manipulation to actually count the exams pcDoseEvaluationIndices = pcTempDoseEvaluationIndices.translate(None, '_') # remove underscores from string piDoseEvaluationIndices = json.loads('[' + pcDoseEvaluationIndices + ']') # this line relies on the function json.loads that converts strings into lists/arrays of integers iDoseEvaluation = len(piDoseEvaluationIndices) - 1 # grab the highest index print("Test1") print("pcMotionState: " + str(pcMotionState)) print("iDoseEvaluation: " + str(iDoseEvaluation)) #if iMotionState == iRefMotionState: # iRefDoseEvaluation = iDoseEvaluation # Store this to retrieve the state dose in the reference state CT # If this is not the reference 4DCT phase, I need to deform to that if pcMotionState != pcRefMotionState: # Also, I need to find the name (index) of the deformation vectors for this motion state pcTempDeformations = patient.Registrations[iDefRegistration].StructureRegistrations._ nDeformations = len(pcTempDeformations.split()) # split() divides the words in a string into a list that I can count using len iFoundDeformation = 0 for iDeformation in range (0,nDeformations): # Testing print(" ") print("TESTING - Retrieval of Deformation Matrix") print("iDeformation: " + str(iDeformation)) print("Register to: " + patient.Registrations[iDefRegistration].StructureRegistrations[iDeformation].ToExamination.Name) print("Register from: " + patient.Registrations[iDefRegistration].StructureRegistrations[iDeformation].FromExamination.Name) if patient.Registrations[iDefRegistration].StructureRegistrations[iDeformation].ToExamination.Name == 'CT ' + str(pcMotionState) and \ patient.Registrations[iDefRegistration].StructureRegistrations[iDeformation].FromExamination.Name == 'CT ' + str(pcRefMotionState): iStructRegistrationToRetrieve = iDeformation # Testing print("I found it at iDeformation: " + str(iDeformation)) print("BTW This is motion state: " + str(iMotionState)) iFoundDeformation = 1 # END if # END for if iFoundDeformation == 0: print(" ") print(" Registration not found.") print(" ") # Warp the state doses onto the reference state CT (I am not 100% sure I am warping the right way, need to validate) with CompositeAction('Deform dose(s)'): patient.MapDose(FractionNumber=0, SetTotalDoseEstimateReference=True, DoseDistribution=patient.TreatmentDelivery.FractionEvaluations[0].DoseOnExaminations[iExamToRetrieve].DoseEvaluations[iDoseEvaluation], StructureRegistration=patient.Registrations[iDefRegistration].StructureRegistrations[iStructRegistrationToRetrieve]) # Keep track of the dose evaluation index of the warped state doses for later retrieval (warped on the reference CT) pcTempWarpDoseEvaluationIndices = patient.TreatmentDelivery.FractionEvaluations[0].DoseOnExaminations[iRefExam].DoseEvaluations._ # This grabs a string containing the dose exam indices, but requires some manipulation to actually count the exams pcWarpDoseEvaluationIndices = pcTempWarpDoseEvaluationIndices.translate(None, '_') # remove underscores from string piWarpDoseEvaluationIndices = json.loads('[' + pcWarpDoseEvaluationIndices + ']') # this line relies on the function json.loads that converts strings into lists/arrays of integers iWarpDoseEvaluation = len(piWarpDoseEvaluationIndices) - 1 # grab the highest index piWarpDoseEvaluationIndex[iMotionState] = iWarpDoseEvaluation # store this for later use # END if this is not the ref 4DCT state # # Write the new dose back into the primary dose grid for this plan # # Determine index of dose evaluation to retrieve if pcMotionState != pcRefMotionState: iStateDoseEvaluationIndex = piWarpDoseEvaluationIndex[iMotionState] else: iStateDoseEvaluationIndex = iDoseEvaluation print("Test2") print("pcMotionState: " + str(pcMotionState)) print("iStateDoseEvaluationIndex: " + str(iStateDoseEvaluationIndex)) # Grab a local copy ppdDoseTemp = 0 # Re-initialize each time ppdDoseTemp = patient.TreatmentDelivery.FractionEvaluations[0] \ .DoseOnExaminations[iRefExam] \ .DoseEvaluations[iStateDoseEvaluationIndex] \ .DoseValues.DoseData if ppdDoseTemp == 0: print(" Cannot find warped dose.") # Write it to primary dose grid patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName].FractionDose.SetDoseValues(Array=ppdDoseTemp, CalculationInfo="State Dose on Ref CT") refCtName = 'CT ' + str(pcRefMotionState) plan = patient.AddNewPlan(PlanName=statePlanNameFinal, PlannedBy="", Comment="", ExaminationName=refCtName, AllowDuplicateNames=False) # Prepare an empty dose grid w same dimensions as the input plan with CompositeAction('Set default grid'): plan.SetDefaultDoseGrid(VoxelSize={ 'x': 0.5, 'y': 0.5, 'z': 0.5 }) dXCorner = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.Corner.x dYCorner = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.Corner.y dZCorner = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.Corner.z dVoxelSizeX = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.VoxelSize.x dVoxelSizeY = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.VoxelSize.y dVoxelSizeZ = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.VoxelSize.z iNVoxelsX = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.NrVoxels.x iNVoxelsY = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.NrVoxels.y iNVoxelsZ = patient.TreatmentPlans[statePlanName].BeamSets[statePlanName].FractionDose.InDoseGrid.NrVoxels.z patient.TreatmentPlans[statePlanNameFinal].UpdateDoseGrid(Corner={ 'x': dXCorner, 'y': dYCorner, 'z': dZCorner }, VoxelSize={ 'x': dVoxelSizeX, 'y': dVoxelSizeY, 'z': dVoxelSizeZ }, NumberOfVoxels={ 'x': iNVoxelsX, 'y': iNVoxelsY, 'z': iNVoxelsZ }) beamset = plan.AddNewBeamSet(Name=statePlanNameFinal, ExaminationName=refCtName, MachineName="ProBeam", NominalEnergy=0, Modality="Protons", TreatmentTechnique="ProtonPencilBeamScanning", PatientPosition="HeadFirstSupine", NumberOfFractions=1, CreateSetupBeams=False, UseLocalizationPointAsSetupIsocenter=False, Comment="") # Write it to primary dose grid patient.TreatmentPlans[statePlanNameFinal] \ .BeamSets[statePlanNameFinal].FractionDose.SetDoseValues(Array=ppdDoseTemp, CalculationInfo="State Dose on Ref CT") # Can't do it directly like this, this is just kept for future notice as a bad way to try # patient.TreatmentPlans[statePlanNameFinal] \ # .BeamSets[statePlanNameFinal].FractionDose = patient.TreatmentPlans[statePlanName] \ # .BeamSets[statePlanName].FractionDose # This else is entered only if (1) this iMotionState is the reference motion state, and (2) the reference plan was planned on this reference motion state else: # Here I only need to do a recalc for the ref state, since I updated the spot weights and beam MU patient.TreatmentPlans[statePlanName] \ .BeamSets[statePlanName].ComputeDose(ComputeBeamDoses=True, DoseAlgorithm="SpotWeightPencilBeam", ForceRecompute=True) # Also, need to rename it to the final name (This is new and not fully tested 040615, needed to have compatibility with new framework) patient.TreatmentPlans[statePlanName].Name = statePlanNameFinal # END if # END for iMotionState # # ## ### Calculate the total 4D dose ## # # # Calculate the total 4D dose (sum the warped state doses), it looks like the summation may be an unscriptable task, might have to take the last step manually, or maybe I can create a temp evaluation dose and overwrite voxel by voxel with the sum of the last number of deformed state doses # If I can grab doses voxel by voxel, I could even just do this by creating a simple copy of the reference plan, name the new plan "plan_4D", calculate first the regular 3D dose, and then overwrite everything # That would give me a nice "4D Plan", as it is nicer to see the plan than the evaluation dose (which e.g., may not have dicom export or whatever feature) print(" ") print("- ") print("-- ") print("--- Calculating Total 4D Dose") print("-- ") print("- ") print(" ") # ELEY 040615 - Here I need to recode so that the new plans use the ref CT state, not the ave CT final4dSumPlanName = str(outputPlanNamePrefix) + refPlanName + '_4D' plan = patient.AddNewPlan(PlanName=final4dSumPlanName, PlannedBy="", Comment="", ExaminationName=refCtName, AllowDuplicateNames=False) # Prepare an empty dose grid w same dimensions as the input plan with CompositeAction('Set default grid'): plan.SetDefaultDoseGrid(VoxelSize={ 'x': 0.5, 'y': 0.5, 'z': 0.5 }) dXCorner = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.Corner.x dYCorner = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.Corner.y dZCorner = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.Corner.z dVoxelSizeX = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.VoxelSize.x dVoxelSizeY = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.VoxelSize.y dVoxelSizeZ = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.VoxelSize.z iNVoxelsX = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.NrVoxels.x iNVoxelsY = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.NrVoxels.y iNVoxelsZ = patient.TreatmentPlans[refPlanName].BeamSets[refPlanName].FractionDose.InDoseGrid.NrVoxels.z patient.TreatmentPlans[final4dSumPlanName].UpdateDoseGrid(Corner={ 'x': dXCorner, 'y': dYCorner, 'z': dZCorner }, VoxelSize={ 'x': dVoxelSizeX, 'y': dVoxelSizeY, 'z': dVoxelSizeZ }, NumberOfVoxels={ 'x': iNVoxelsX, 'y': iNVoxelsY, 'z': iNVoxelsZ }) beamset = plan.AddNewBeamSet(Name=final4dSumPlanName, ExaminationName=refCtName, MachineName="ProBeam", NominalEnergy=0, Modality="Protons", TreatmentTechnique="ProtonPencilBeamScanning", PatientPosition="HeadFirstSupine", NumberOfFractions=1, CreateSetupBeams=False, UseLocalizationPointAsSetupIsocenter=False, Comment="") # Modify list to only include motion states that have any dose (trying to fix the "[state plan] has no computed fraction dose" problem) ppcStatePlanNamesWithDose = [0] * nMotionStatesWithDose iCount = 0 for iMotionState in range(0,nMotionStates): if piFlagMotionStateHasDose[iMotionState] == 1: ppcStatePlanNamesWithDose[iCount] = ppcStatePlanNames[iMotionState] iCount = iCount + 1 print ppcStatePlanNamesWithDose # Sum state plan doses to calc total 4D dose and overwrite dWeight = 1 # 1 gives single fraction 4D dose (this makes most sense), but could also scale to nFractions if you want to get total RxDose pdWeights = [dWeight] * nMotionStatesWithDose print pdWeights patient.TreatmentPlans[final4dSumPlanName] \ .BeamSets[final4dSumPlanName].FractionDose.SetWeightedDose(PlanNames=ppcStatePlanNamesWithDose, Weights=pdWeights)