diff --git a/examples/PlanOptimization/run_4DprotonEvaluation.py b/examples/PlanEvaluation/run_4DprotonEvaluation.py similarity index 100% rename from examples/PlanOptimization/run_4DprotonEvaluation.py rename to examples/PlanEvaluation/run_4DprotonEvaluation.py diff --git a/examples/PlanOptimization/run_evaluatePhotonRobustness.py b/examples/PlanEvaluation/run_evaluatePhotonRobustness.py similarity index 100% rename from examples/PlanOptimization/run_evaluatePhotonRobustness.py rename to examples/PlanEvaluation/run_evaluatePhotonRobustness.py diff --git a/examples/PlanOptimization/run_evaluateProtonRobustness.py b/examples/PlanEvaluation/run_evaluateProtonRobustness.py similarity index 100% rename from examples/PlanOptimization/run_evaluateProtonRobustness.py rename to examples/PlanEvaluation/run_evaluateProtonRobustness.py diff --git a/examples/PlanOptimization/run_4DProtonOptimization.py b/examples/PlanOptimization/run_4DProtonOptimization.py index dd8bb18..5f3538f 100644 --- a/examples/PlanOptimization/run_4DProtonOptimization.py +++ b/examples/PlanOptimization/run_4DProtonOptimization.py @@ -36,7 +36,6 @@ from opentps.core.data.plan._protonPlanDesign import ProtonPlanDesign 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.dataLoader import readData from opentps.core.io.scannerReader import readScanner @@ -47,6 +46,7 @@ from opentps.core.processing.planOptimization.planOptimization import IntensityModulationOptimizer from opentps.core.processing.imageProcessing.resampler3D import resampleImage3DOnImage3D from opentps.core.io.dicomIO import writeRTDose, readDicomDose +import opentps.core.processing.planOptimization.objectives.dosimetricObjectives as doseObj logger = logging.getLogger(__name__) @@ -68,8 +68,8 @@ #%% #CT and ROI creation #--------------------------- -# Generic example: 4DCT composed of 3 CTs : 2 phases and the MidP. -# The anatomy consists of a square target moving vertically, with an organ at risk and soft tissue (muscle) in front of it. +# Generic example: 4DCT composed of 3 CTs : 2 phases and the MidP. +# The anatomy consists of a square target moving vertically, with an organ at risk and soft tissue (muscle) in front of it. CT4D = [] ROI4D = [] for i in range(0, 3): @@ -87,7 +87,7 @@ Patient.id = f'12082024' Patient.birthDate = dt.strftime('%Y%m%d') patient.sex = "" - + ctSize = 150 ct = CTImage(seriesInstanceUID=ctSeriesInstanceUID, frameOfReferenceUID=frameOfReferenceUID) ct.name = f'CT_Phase_{i}' @@ -130,7 +130,7 @@ data[25:45, 50:80, 65:85] = True TV.imageArray = data ROI.append(TV) - + # Muscle Muscle = ROIMask() Muscle.patient = patient @@ -225,8 +225,8 @@ # planDesign.robustness.RandomPeriodError = 5.0 # % # planDesign.robustness.Breathing_period = 1 # x100% # default value - planDesign.spotSpacing = 10.0 - planDesign.layerSpacing = 10.0 + planDesign.spotSpacing = 10.0 + planDesign.layerSpacing = 10.0 planDesign.targetMargin = 7 # Enough to encompass target motion planDesign.defineTargetMaskAndPrescription(target = RefTV, targetPrescription = 60.) @@ -234,7 +234,7 @@ plan = planDesign.buildPlan() plan.rtPlanName = f"RobustPlan_4D" - # refIndex : + # refIndex : # ACCUMULATED -> Index of the Image in the 4DCT one wish we will accumulate the dose. ## SYSTEMATIC -> Index of the Image in the 4DCT who will be used as the nominal. So the one closer to the MidP. Or the Midp. @@ -248,10 +248,10 @@ #%% # Set objectives #---------------------- -plan.planDesign.objectives.addFidObjective(RefTV, FidObjective.Metrics.DMAX, limitValue = 63.0, weight = 100.0, robust=True) -plan.planDesign.objectives.addFidObjective(RefTV, FidObjective.Metrics.DMIN, limitValue = 60.0, weight = 100.0, robust=True) -plan.planDesign.objectives.addFidObjective(RefOAR, FidObjective.Metrics.DMAX, limitValue = 40.0, weight = 80.0) -plan.planDesign.objectives.addFidObjective(RefBody, FidObjective.Metrics.DMAX, limitValue = 40.0, weight = 80.0) +plan.planDesign.objectives.addObjective(doseObj.DMax(RefTV,63, weight=100.0, robust=True)) +plan.planDesign.objectives.addObjective(doseObj.DMin(RefTV,60, weight=100.0, robust=True)) +plan.planDesign.objectives.addObjective(doseObj.DMax(RefOAR, 40, weight=80.0)) +plan.planDesign.objectives.addObjective(doseObj.DMax(RefBody, 40, weight=80.0)) #%% # Optimize treatment plan diff --git a/examples/PlanOptimization/run_SimpleOptimizationPhoton.py b/examples/PlanOptimization/run_SimpleOptimizationPhoton.py index 4327ac1..6d387db 100644 --- a/examples/PlanOptimization/run_SimpleOptimizationPhoton.py +++ b/examples/PlanOptimization/run_SimpleOptimizationPhoton.py @@ -41,7 +41,6 @@ from opentps.core.data.images import ROIMask from opentps.core.data import DVH from opentps.core.data import Patient -from opentps.core.data.plan import FidObjective from opentps.core.io.scannerReader import readScanner from opentps.core.io.serializedObjectIO import saveRTPlan, loadRTPlan from opentps.core.processing.doseCalculation.doseCalculationConfig import DoseCalculationConfig @@ -49,6 +48,7 @@ from opentps.core.processing.planOptimization.planOptimization import IntensityModulationOptimizer from opentps.core.processing.doseCalculation.photons.cccDoseCalculator import CCCDoseCalculator from opentps.core.data.plan import PhotonPlanDesign +import opentps.core.processing.planOptimization.objectives.dosimetricObjectives as doseObj logger = logging.getLogger(__name__) @@ -88,6 +88,12 @@ data[100:120, 100:120, 100:120] = True roi.imageArray = data +body = roi.copy() +body.name = 'Body' +body.dilateMask(20) +body.imageArray = np.logical_xor(body.imageArray, roi.imageArray).astype(bool) + + #%% #Output path #----------- @@ -127,21 +133,38 @@ planDesign.xBeamletSpacing_mm = 5 planDesign.yBeamletSpacing_mm = 5 planDesign.targetMargin = 5.0 - planDesign.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) - - plan = planDesign.buildPlan() + planDesign.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) - beamlets = ccc.computeBeamlets(ct, plan) + plan = planDesign.buildPlan() + + beamlets = ccc.computeBeamlets(ct, plan) doseInfluenceMatrix = copy.deepcopy(beamlets) - + plan.planDesign.beamlets = beamlets beamlets.storeOnFS(os.path.join(output_path, "BeamletMatrix_" + plan.seriesInstanceUID + ".blm")) # Save plan with initial spot weights in serialized format (OpenTPS format) saveRTPlan(plan, plan_file) +plan.planDesign.objectives.addObjective(doseObj.DMax(body,5, weight=1.0)) + +plan.planDesign.objectives.addObjective(doseObj.DMax(roi, 21, weight=10.0)) +plan.planDesign.objectives.addObjective(doseObj.DMin(roi, 20, weight=20.0)) + +# Other examples of objectives + +# plan.planDesign.objectives.addObjective(doseObj.DUniform(roi, 20, weight=10)) +# +# plan.planDesign.objectives.addObjective(doseObj.DMaxMean(roi,21,weight=5)) +# plan.planDesign.objectives.addObjective(doseObj.DMinMean(roi,20,weight=5)) +# plan.planDesign.objectives.addObjective(doseObj.DFallOff(roi,oar,10,5,15,weight=5)) +# +# plan.planDesign.objectives.addObjective(doseObj.DVHMax(roi, 21,0.05, weight=5)) +# plan.planDesign.objectives.addObjective(doseObj.DVHMin(roi, 18,0.90, weight=5)) +# +# plan.planDesign.objectives.addObjective(doseObj.EUDMin(roi, 19, 1, weight=50)) +# plan.planDesign.objectives.addObjective(doseObj.EUDMax(roi, 21, 1, weight=50)) +# plan.planDesign.objectives.addObjective(doseObj.EUDUniform(roi, 20, 1, weight=50)) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMAX, 20.0, 1.0) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMIN, 20.5, 1.0) plan.numberOfFractionsPlanned = 30 @@ -166,6 +189,7 @@ #--------------------------------- target_DVH = DVH(roi, doseImage) +body_DVH = DVH(body, doseImage) print('D5 - D95 = {} Gy'.format(target_DVH.D5 - target_DVH.D95)) clinROI = [roi.name, roi.name] clinMetric = ["Dmin", "Dmax"] @@ -178,13 +202,14 @@ # center of mass #--------------- roi = resampleImage3DOnImage3D(roi, ct) +body = resampleImage3DOnImage3D(body,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_mask = roi.imageArray[:, :, Z_coord].transpose(1, 0) +img_body = body.imageArray[:, :, Z_coord].transpose(1, 0) img_dose = resampleImage3DOnImage3D(doseImage, ct) img_dose = img_dose.imageArray[:, :, Z_coord].transpose(1, 0) @@ -195,10 +220,13 @@ 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 +ax[0].contour(img_body,[0.5],colors='green') # body +ax[0].contour(img_mask,[0.5],colors='red') # 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].plot(target_DVH.histogram[0], target_DVH.histogram[1], label=target_DVH.name,color='red') +ax[1].plot(body_DVH.histogram[0], body_DVH.histogram[1], label=body_DVH.name,color='green') ax[1].set_xlabel("Dose (Gy)") ax[1].set_ylabel("Volume (%)") ax[1].grid(True) diff --git a/examples/PlanOptimization/run_SimpleOptimizationProton.py b/examples/PlanOptimization/run_SimpleOptimizationProton.py index 3fceeda..e293da0 100644 --- a/examples/PlanOptimization/run_SimpleOptimizationProton.py +++ b/examples/PlanOptimization/run_SimpleOptimizationProton.py @@ -30,14 +30,13 @@ #%% #import the needed opentps.core packages - +import opentps.core.processing.planOptimization.objectives.dosimetricObjectives as doseObj 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 ProtonPlanDesign 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.io.serializedObjectIO import saveRTPlan, loadRTPlan @@ -81,6 +80,20 @@ data[100:120, 100:120, 100:120] = True roi.imageArray = data +# roi2 = ROIMask() +# roi2.patient = patient +# roi2.name = 'TV2' +# roi2.color = (0, 0, 255) # blue +# data = np.zeros((ctSize, ctSize, ctSize)).astype(bool) +# data[100:120, 70:100, 70:120] = True +# roi2.imageArray = data + + +body = roi.copy() +body.name = 'Body' +body.dilateMask(20) +body.imageArray = np.logical_xor(body.imageArray, roi.imageArray).astype(bool) + #%% #Configure dose engine #--------------------- @@ -91,10 +104,10 @@ mc2.ctCalibration = ctCalibration mc2._independentScoringGrid = True -scoringSpacing = [2, 2, 2] +scoringSpacing = [2,2,2] mc2._scoringVoxelSpacing = scoringSpacing -#%% +#%%# #Design plan #----------- @@ -108,48 +121,55 @@ planInit.beamNames = beamNames planInit.couchAngles = couchAngles planInit.calibration = ctCalibration -planInit.spotSpacing = 6.0 -planInit.layerSpacing = 6.0 -planInit.targetMargin = 0.0 +planInit.spotSpacing = 5.0 +planInit.layerSpacing = 5.0 +planInit.targetMargin = 2.0 planInit.setScoringParameters(scoringSpacing=[2, 2, 2], adapt_gridSize_to_new_spacing=True) # needs to be called after scoringGrid settings but prior to spot placement -planInit.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) - +planInit.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) + plan = planInit.buildPlan() # Spot placement plan.PlanName = "NewPlan" -beamlets = mc2.computeBeamlets(ct, plan, roi=[roi]) +beamlets = mc2.computeBeamlets(ct, plan, roi=[roi,body]) plan.planDesign.beamlets = beamlets # doseImageRef = beamlets.toDoseImage() #%% #objectives #---------- +plan.planDesign.objectives.addObjective(doseObj.DMax(body,5, weight=1.0)) + +plan.planDesign.objectives.addObjective(doseObj.DMax(roi, 21, weight=10.0)) +plan.planDesign.objectives.addObjective(doseObj.DMin(roi, 20, weight=20.0)) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMAX, 20.0, 1.0) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMIN, 20.5, 1.0) # Other examples of objectives -# plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMEAN, 20, 1.0) -# plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DUNIFORM, 20, 1.0) -# plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DVHMIN, 19, 1.0, volume = 95) -# plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DVHMAX, 21, 1.0, volume = 5) -# plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.EUDMIN, 19.5, 1.0, EUDa = 0.2) -# plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.EUDMAX, 20, 1.0, EUDa = 1) -# plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.EUDUNIFORM, 20.5, 1.0, EUDa = 0.5) -# plan.planDesign.objectives.addFidObjective(BODY, FidObjective.Metrics.DFALLOFF, weight=10, fallOffDistance=1, fallOffLowDoseLevel=0, fallOffHighDoseLevel=21) + +# plan.planDesign.objectives.addObjective(doseObj.DUniform(roi, 20, weight=10)) +# +# plan.planDesign.objectives.addObjective(doseObj.DMaxMean(roi,21,weight=5)) +# plan.planDesign.objectives.addObjective(doseObj.DMinMean(roi,20,weight=5)) +# plan.planDesign.objectives.addObjective(doseObj.DFallOff(roi,oar,10,5,15,weight=5)) +# +# plan.planDesign.objectives.addObjective(doseObj.DVHMax(roi, 21,0.05, weight=5)) +# plan.planDesign.objectives.addObjective(doseObj.DVHMin(roi, 18,0.90, weight=5)) +# +# plan.planDesign.objectives.addObjective(doseObj.EUDMin(roi, 19, 1, weight=50)) +# plan.planDesign.objectives.addObjective(doseObj.EUDMax(roi, 21, 1, weight=50)) +# plan.planDesign.objectives.addObjective(doseObj.EUDUniform(roi, 20, 1, weight=50)) #%% #Optimize plan #------------- -solver = IntensityModulationOptimizer(method='Scipy_L-BFGS-B', plan=plan, maxiter=50) +solver = IntensityModulationOptimizer(method='Scipy_L-BFGS-B', plan=plan, maxiter=50, hardwareAcceleration = None) doseImage, ps = solver.optimize() #%% #Final dose computation #---------------------- -mc2.nbPrimaries = 1e7 +mc2.nbPrimaries = 1e6 doseImage = mc2.computeDose(ct, plan) #%% @@ -162,16 +182,19 @@ 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)) +bodyResampled = resampleImage3D(body, origin=ct.origin, spacing=scoringSpacing) +body_DVH = DVH(bodyResampled, doseImage) # center of mass roi = resampleImage3DOnImage3D(roi, ct) +body = resampleImage3DOnImage3D(body,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_mask = roi.imageArray[:, :, Z_coord].transpose(1, 0) +img_body = body.imageArray[:, :, Z_coord].transpose(1, 0) img_dose = resampleImage3DOnImage3D(doseImage, ct) img_dose = img_dose.imageArray[:, :, Z_coord].transpose(1, 0) @@ -183,10 +206,12 @@ # Display dose fig, ax = plt.subplots(1, 2, figsize=(12, 5)) ax[0].imshow(img_ct, cmap='gray') -ax[0].imshow(img_mask, alpha=.2, cmap='binary') # PTV +ax[0].contour(img_body,[0.5],colors='green') # Body +ax[0].contour(img_mask,[0.5],colors='red') # 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].plot(target_DVH.histogram[0], target_DVH.histogram[1], label=target_DVH.name,color='red') +ax[1].plot(body_DVH.histogram[0], body_DVH.histogram[1], label=body_DVH.name,color='green') ax[1].set_xlabel("Dose (Gy)") ax[1].set_ylabel("Volume (%)") plt.grid(True) diff --git a/examples/PlanOptimization/run_beamletFreeOpti.py b/examples/PlanOptimization/run_beamletFreeOpti.py index 5362257..a6a9280 100644 --- a/examples/PlanOptimization/run_beamletFreeOpti.py +++ b/examples/PlanOptimization/run_beamletFreeOpti.py @@ -34,12 +34,12 @@ from opentps.core.data.plan._protonPlanDesign import ProtonPlanDesign 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.protons.mcsquareDoseCalculator import MCsquareDoseCalculator from opentps.core.processing.imageProcessing.resampler3D import resampleImage3DOnImage3D +import opentps.core.processing.planOptimization.objectives.dosimetricObjectives as doseObj #%% #Output path @@ -90,6 +90,11 @@ data[100:120, 100:120, 100:120] = True roi.imageArray = data +body = roi.copy() +body.name = 'body' +body.dilateMask(20) +body.imageArray = np.logical_xor(body.imageArray, roi.imageArray).astype(bool) + image = plt.imshow(ct.imageArray[110,:,:],cmap='Blues') plt.colorbar(image) plt.contour(roi.imageArray[110,:,:],colors="red") @@ -129,31 +134,33 @@ planDesign.beamNames = beamNames planDesign.couchAngles = couchAngles planDesign.calibration = ctCalibration -planDesign.spotSpacing = 5.0 -planDesign.layerSpacing = 5.0 -planDesign.targetMargin = 5.0 +planDesign.spotSpacing = 3.0 +planDesign.layerSpacing = 3.0 +planDesign.targetMargin = 1.0 # needs to be called prior to spot placement -planDesign.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) - +planDesign.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) + plan = planDesign.buildPlan() # Spot placement plan.PlanName = "NewPlan" -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.objectives.addObjective(doseObj.DMax(body,5, weight=1.0)) +plan.planDesign.objectives.addObjective(doseObj.DMax(roi, 20, weight=10.0)) +plan.planDesign.objectives.addObjective(doseObj.DMin(roi, 20, weight=10.0)) #%% #Mcsquare beamlet free planOptimization #-------------------------------------- #Now that we have every needed objects we can compute the optimization through MCsquare. :warning: It may take some time to compute. -doseImage = mc2.optimizeBeamletFree(ct, plan, [roi]) +doseImage = mc2.optimizeBeamletFree(ct, plan, [roi,body]) #%% #Dose volume histogram #--------------------- target_DVH = DVH(roi, doseImage) +body_DVH = DVH(body, 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)) @@ -169,8 +176,8 @@ 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_mask = roi.imageArray[:, :, Z_coord].transpose(1, 0) +img_body = body.imageArray[:, :, Z_coord].transpose(1, 0) img_dose = resampleImage3DOnImage3D(doseImage, ct) img_dose = img_dose.imageArray[:, :, Z_coord].transpose(1, 0) @@ -182,10 +189,12 @@ 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 +ax[0].contour(img_body,[0.5],colors='green') # Body +ax[0].contour(img_mask,[0.5],colors='red') # 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].plot(target_DVH.histogram[0], target_DVH.histogram[1], label=target_DVH.name,color='red') +ax[1].plot(body_DVH.histogram[0], body_DVH.histogram[1], label=body_DVH.name,color='green') ax[1].set_xlabel("Dose (Gy)") ax[1].set_ylabel("Volume (%)") ax[0].set_title("Computed dose") diff --git a/examples/PlanOptimization/run_boundConstraintsOpti.py b/examples/PlanOptimization/run_boundConstraintsOpti.py index 2eeadfd..f9d869e 100644 --- a/examples/PlanOptimization/run_boundConstraintsOpti.py +++ b/examples/PlanOptimization/run_boundConstraintsOpti.py @@ -34,14 +34,13 @@ #%% #import the needed opentps.core packages - +import opentps.core.processing.planOptimization.objectives.dosimetricObjectives as doseObj 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 ProtonPlanDesign 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.io.serializedObjectIO import saveRTPlan, loadRTPlan @@ -97,6 +96,11 @@ data[100:120, 100:120, 100:120] = True roi.imageArray = data +body = roi.copy() +body.name = 'Body' +body.dilateMask(20) +body.imageArray = np.logical_xor(body.imageArray, roi.imageArray).astype(bool) + #%% #Configuration of Mcsquare #------------------------- @@ -111,7 +115,7 @@ #Plan Creation #------------- -# Design plan +# Design plan beamNames = ["Beam1"] gantryAngles = [0.] couchAngles = [0.] @@ -134,12 +138,12 @@ planInit.targetMargin = 5.0 planInit.setScoringParameters(scoringSpacing=[2, 2, 2], adapt_gridSize_to_new_spacing=True) # needs to be called after scoringGrid settings but prior to spot placement - planInit.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) - + planInit.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) + plan = planInit.buildPlan() # Spot placement plan.PlanName = "NewPlan" - beamlets = mc2.computeBeamlets(ct, plan, roi=[roi]) + beamlets = mc2.computeBeamlets(ct, plan, roi=[roi,body]) plan.planDesign.beamlets = beamlets beamlets.storeOnFS(os.path.join(output_path, "BeamletMatrix_" + plan.seriesInstanceUID + ".blm")) @@ -147,8 +151,10 @@ saveRTPlan(plan, plan_file) # Set objectives (attribut is already initialized in planDesign object) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMAX, 20.0, 1.0) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMIN, 20.5, 1.0) +plan.planDesign.objectives.addObjective(doseObj.DMax(body,5, weight=1.0)) + +plan.planDesign.objectives.addObjective(doseObj.DMax(roi, 21, weight=10.0)) +plan.planDesign.objectives.addObjective(doseObj.DMin(roi, 20, weight=20.0)) solver = BoundConstraintsOptimizer(method='Scipy_L-BFGS-B', plan=plan, maxiter=50, bounds=(0.2, 50)) @@ -162,6 +168,7 @@ #--------------------- target_DVH = DVH(roi, doseImage) +body_DVH = DVH(body, 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)) @@ -172,13 +179,14 @@ #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) +body = resampleImage3DOnImage3D(body,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_mask = roi.imageArray[:, :, Z_coord].transpose(1, 0) +img_body = body.imageArray[:, :, Z_coord].transpose(1, 0) img_dose = resampleImage3DOnImage3D(doseImage, ct) img_dose = img_dose.imageArray[:, :, Z_coord].transpose(1, 0) @@ -188,14 +196,15 @@ fig, ax = plt.subplots(1, 2, figsize=(12, 5)) ax[0].imshow(img_ct, cmap='gray') -ax[0].imshow(img_mask, alpha=.2, cmap='binary') +ax[0].contour(img_body,[0.5],colors='green') # Body +ax[0].contour(img_mask,[0.5],colors='red') # 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].plot(target_DVH.histogram[0], target_DVH.histogram[1], label=target_DVH.name,color='red') +ax[1].plot(body_DVH.histogram[0], body_DVH.histogram[1], label=body_DVH.name,color='green') ax[1].set_xlabel("Dose (Gy)") ax[1].set_ylabel("Volume (%)") plt.grid(True) plt.legend() -plt.show() -plt.savefig(f'{output_path}/Dose_BoundContraintOpti.png', format = 'png') +plt.savefig(os.path.join(output_path, 'SimpleOpti1.png'),format = 'png') plt.show() \ No newline at end of file diff --git a/examples/PlanOptimization/run_robustOptimizationProtons.py b/examples/PlanOptimization/run_robustOptimizationProtons.py index 8c83bdf..075323c 100644 --- a/examples/PlanOptimization/run_robustOptimizationProtons.py +++ b/examples/PlanOptimization/run_robustOptimizationProtons.py @@ -30,13 +30,13 @@ #%% # import the needed opentps.core packages +import opentps.core.processing.planOptimization.objectives.dosimetricObjectives as doseObj from opentps.core.data.images import CTImage from opentps.core.data.images import ROIMask from opentps.core.data.plan._protonPlanDesign import ProtonPlanDesign from opentps.core.data.plan import RobustnessProton 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.io.serializedObjectIO import loadRTPlan, saveRTPlan @@ -88,6 +88,11 @@ data[100:120, 100:120, 100:120] = True roi.imageArray = data +body = roi.copy() +body.name = 'Body' +body.dilateMask(20) +body.imageArray = np.logical_xor(body.imageArray, roi.imageArray).astype(bool) + #%% # Design plan #---------------- @@ -136,7 +141,7 @@ # All scenarios (includes diagonals on sphere) # planDesign.robustness.selectionStrategy = planDesign.robustness.Strategies.ALL - # Random scenario sampling + # Random scenario sampling # planDesign.robustness.selectionStrategy = planDesign.robustness.Strategies.RANDOM planDesign.robustness.numScenarios = 5 # specify how many random scenarios to simulate, default = 100 @@ -149,11 +154,11 @@ plan = planDesign.buildPlan() # Spot placement plan.PlanName = "RobustPlan" - nominal, scenarios = mc2.computeRobustScenarioBeamlets(ct, plan, roi=[roi], storePath=output_path) + nominal, scenarios = mc2.computeRobustScenarioBeamlets(ct, plan, roi=[roi,body], storePath=output_path) plan.planDesign.beamlets = nominal plan.planDesign.robustness.scenarios = scenarios plan.planDesign.robustness.numScenarios = len(scenarios) - + #saveRTPlan(plan, plan_file) @@ -164,8 +169,10 @@ #%% # Set objectives (attribut is already initialized in planDesign object) #---------------------- -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMAX, 20.0, 1.0, robust=True) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMIN, 20.5, 1.0, robust=True) +plan.planDesign.objectives.addObjective(doseObj.DMax(body,5, weight=1.0)) + +plan.planDesign.objectives.addObjective(doseObj.DMax(roi, 21, weight=10.0)) +plan.planDesign.objectives.addObjective(doseObj.DMin(roi, 20, weight=20.0)) solver = IntensityModulationOptimizer(method='Scipy_L-BFGS-B', plan=plan, maxiter=50) @@ -177,10 +184,15 @@ plan_file = os.path.join(output_path, "Plan_Proton_WaterPhantom_cropped_optimized.tps") saveRTPlan(plan, plan_file, unloadBeamlets=False) +mc2.nbPrimaries = 1e6 +doseImage = mc2.computeDose(ct, plan) + #%% # Compute DVH #---------------------- target_DVH = DVH(roi, doseImage) +body_DVH = DVH(body, 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)) @@ -189,13 +201,14 @@ # Center of mass #---------------------- roi = resampleImage3DOnImage3D(roi, ct) +body = resampleImage3DOnImage3D(body,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_mask = roi.imageArray[:, :, Z_coord].transpose(1, 0) +img_body = body.imageArray[:, :, Z_coord].transpose(1, 0) img_dose = resampleImage3DOnImage3D(doseImage, ct) img_dose = img_dose.imageArray[:, :, Z_coord].transpose(1, 0) @@ -206,10 +219,12 @@ 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 +ax[0].contour(img_body,[0.5],colors='green') # Body +ax[0].contour(img_mask,[0.5],colors='red') # 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].plot(target_DVH.histogram[0], target_DVH.histogram[1], label=target_DVH.name,color='red') +ax[1].plot(body_DVH.histogram[0], body_DVH.histogram[1], label=body_DVH.name,color='green') ax[1].set_xlabel("Dose (Gy)") ax[1].set_ylabel("Volume (%)") plt.grid(True) diff --git a/examples/PlanOptimization/run_simpleOptimization_createDicomStudy.py b/examples/PlanOptimization/run_simpleOptimization_createDicomStudy.py index 8dea16e..a0de5e8 100644 --- a/examples/PlanOptimization/run_simpleOptimization_createDicomStudy.py +++ b/examples/PlanOptimization/run_simpleOptimization_createDicomStudy.py @@ -34,7 +34,7 @@ #%% #import the needed opentps.core packages - +import opentps.core.processing.planOptimization.objectives.dosimetricObjectives as doseObj from opentps.core.io.dicomIO import writeRTPlan, writeDicomCT, writeRTDose, writeRTStruct from opentps.core.processing.planOptimization.tools import evaluateClinical from opentps.core.data.images import CTImage, DoseImage @@ -44,7 +44,6 @@ from opentps.core.data import DVH from opentps.core.data import Patient from opentps.core.data import RTStruct -from opentps.core.data.plan import FidObjective from opentps.core.io import mcsquareIO from opentps.core.io.scannerReader import readScanner from opentps.core.io.serializedObjectIO import saveRTPlan, loadRTPlan @@ -123,6 +122,11 @@ data[100:120, 100:120, 100:120] = True roi.imageArray = data +body = roi.copy() +body.name = 'Body' +body.dilateMask(20) +body.imageArray = np.logical_xor(body.imageArray, roi.imageArray).astype(bool) + #%% #Configuration of Mcsquare #------------------------- @@ -161,11 +165,11 @@ planDesign.targetMargin = 5.0 planDesign.setScoringParameters(scoringSpacing=[2, 2, 2], adapt_gridSize_to_new_spacing=True) planDesign.defineTargetMaskAndPrescription(target = roi, targetPrescription = 20.) # needs to be called prior spot placement - + plan = planDesign.buildPlan() # Spot placement plan.rtPlanName = "Simple_Patient" - beamlets = mc2.computeBeamlets(ct, plan, roi=[roi]) + beamlets = mc2.computeBeamlets(ct, plan, roi=[roi,body]) plan.planDesign.beamlets = beamlets beamlets.storeOnFS(os.path.join(output_path, "BeamletMatrix_" + plan.seriesInstanceUID + ".blm")) # Save plan with initial spot weights in serialized format (OpenTPS format) @@ -177,8 +181,10 @@ #---------- # Set objectives (attribut is already initialized in planDesign object) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMAX, 20.0, 1.0) -plan.planDesign.objectives.addFidObjective(roi, FidObjective.Metrics.DMIN, 20.5, 1.0) +plan.planDesign.objectives.addObjective(doseObj.DMax(body,5, weight=1.0)) + +plan.planDesign.objectives.addObjective(doseObj.DMax(roi, 21, weight=10.0)) +plan.planDesign.objectives.addObjective(doseObj.DMin(roi, 20, weight=20.0)) #%% #Optimize plan @@ -205,6 +211,7 @@ #--------------------- target_DVH = DVH(roi, doseImage) +body_DVH = DVH(body, doseImage) print('D5 - D95 = {} Gy'.format(target_DVH.D5 - target_DVH.D95)) clinROI = [roi.name, roi.name] clinMetric = ["Dmin", "Dmax"] @@ -217,7 +224,7 @@ doseImage.referenceCT = ct doseImage.patient = patient doseImage.studyInstanceUID = studyInstanceUID -doseImage.frameOfReferenceUID = frameOfReferenceUID +doseImage.frameOfReferenceUID = frameOfReferenceUID doseImage.sopClassUID = '1.2.840.10008.5.1.4.1.1.481.2' doseImage.mediaStorageSOPClassUID = '1.2.840.10008.5.1.4.1.1.481.2' doseImage.sopInstanceUID = pydicom.uid.generate_uid() @@ -242,8 +249,8 @@ 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_mask = roi.imageArray[:, :, Z_coord].transpose(1, 0) +img_body = body.imageArray[:, :, Z_coord].transpose(1, 0) img_dose = resampleImage3DOnImage3D(doseImage, ct) img_dose = img_dose.imageArray[:, :, Z_coord].transpose(1, 0) @@ -255,10 +262,12 @@ 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 +ax[0].contour(img_body,[0.5],colors='green') # Body +ax[0].contour(img_mask,[0.5],colors='red') 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].plot(target_DVH.histogram[0], target_DVH.histogram[1], label=target_DVH.name,color='red') +ax[1].plot(body_DVH.histogram[0], body_DVH.histogram[1], label=body_DVH.name,color='green') ax[1].set_xlabel("Dose (Gy)") ax[1].set_ylabel("Volume (%)") ax[1].grid(True)