Example: TG-119 standard treatment plan optimization
Intro
Welcome to the pyanno4rt example notebook!
In this notebook, we will showcase the core functionality of our package using data from the TG-119 standard case (available from our Github repository as .mat-files). The first part will present a beginner-friendly version of the code-based interface, followed by the UI-based interface in the second part.
Import of the relevant classes
First, we import the base classes. Our package is designed for clarity and ease of use, wherefore it has only one class for initializing a treatment plan and one class for initializing the graphical user interface. Hence, the import statements are:
from pyanno4rt.base import TreatmentPlan
from pyanno4rt.gui import GraphicalUserInterface
Code-based interface
If you prefer to work with the command line interface (CLI) or an interactive development environment (IDE), you can initialize the TreatmentPlan class by hand. The parameter space of this class is divided into three parameter groups:
- configuration parameters: design parameters w.r.t general or external data settings for the treatment plan
- optimization parameters: design parameters w.r.t the components (objectives and constraints), method and solver for treatment plan optimization
- evaluation parameters: design parameters w.r.t the evaluation methods used for treatment plan assessment
Well then, let’s create an instance of the TreatmentPlan class!
Treatment plan initialization
For the sake of readability, we will define the parameter groups one by one (of course, you could also directly specify them in the base class arguments). Our package utilizes Python dictionaries for this purpose, which allow an efficient mapping between parameter names and values per group and promote a transparent setup and passing.
Setting up the configuration dictionary
We decide to label our plan ‘TG-119-example’ and set the minimum logging level to ‘info’, which means that any debugging messages will be suppressed. For the modality and the number of fractions, we stick to the default values ‘photon’ and 30. Since we have some MATLAB files available for the TG-119 case, we provide the corresponding paths to the imaging and dose-influence matrix files (you may adapt them). Post-processing interpolation of the imaging data is not required, so we leave the parameter at None. Finally, we know that the dose-influence matrix has been calculated with a resolution of 6 mm in each dimension, so we set the dose resolution parameter accordingly.
configuration = {
'label': 'TG-119-example', # Unique identifier for the treatment plan
'min_log_level': 'info', # Minimum logging level
'modality': 'photon', # Treatment modality
'number_of_fractions': 30, # Number of fractions
'imaging_path': './TG_119_data.mat', # Path to the CT and segmentation data
'target_imaging_resolution': None, # Imaging resolution for post-processing interpolation of the CT and segmentation data
'dose_matrix_path': './TG_119_photonDij.mat', # Path to the dose-influence matrix
'dose_resolution': [6, 6, 6] # Size of the dose grid in [mm] per dimension
}
Great, we have completely defined the first parameter group 👍
Setting up the optimization dictionary
Next, we need to describe how the TG-119 treatment plan should be optimized. In general, the final plan should apply a reasonably high dose to the target volumes while limiting the dose exposure to relevant organs at risk to prevent post-treatment complications.
To achieve this, we define objective functions for the core (‘Core’), for the outer target (‘OuterTarget’), and for the whole body (‘BODY’), where ‘Squared Overdosing’ refers to a function that penalizes dose values above a maximum, and ‘Squared Deviation’ refers to a function that penalizes upward and downward deviations from a target. The definition of these functions is again based on a dictionary, the components dictionary: it takes the segment names from the imaging data as keys and sub-dictionaries as values, in which the component type (‘objective’ or ‘constraint’) and the component instance (or a list of instances) with the component’s class name and parameters are set.
Once the components have been defined, we find ourselves in a trade-off situation, where a higher degree of fulfillment for one objective is usually accompanied with a lower degree of fulfillment for another. We can handle this by choosing the ‘weighted-sum’ method, which bypasses the multi-objective problem by multiplying each objective value with a weight parameter and then summing them up, effectively merging them into a scalar “total” objective function. This works well with the default solution algorithm, the ‘L-BFGS-B’ algorithm from the ‘scipy’ solver, so we pick that one. For the initialization of the fluence vector (holding the decision variables), we opt for ‘target-coverage’ to start off with a satisfactory dose level for the outer target (alternatively we could have passed ‘warm-start’ and replaced None for the initial fluence vector with an array). We place a lower bound of 0 and no upper bound (None) on the fluence, matching its physical properties. As the final step, we limit the number of iterations to 500 and the tolerance (precision goal) for the objective function value to 0.001.
optimization = {
'components': { # Optimization components for each segment of interest
'Core': {
'type': 'objective',
'instance': {
'class': 'Squared Overdosing',
'parameters': {
'maximum_dose': 25,
'weight': 100
}
}
},
'OuterTarget': {
'type': 'objective',
'instance': {
'class': 'Squared Deviation',
'parameters': {
'target_dose': 60,
'weight': 1000
}
}
},
'BODY': {
'type': 'objective',
'instance': {
'class': 'Squared Overdosing',
'parameters': {
'maximum_dose': 30,
'weight': 800
}
}
}
},
'method': 'weighted-sum', # Single- or multi-criteria optimization method
'solver': 'scipy', # Python package to be used for solving the optimization problem
'algorithm': 'L-BFGS-B', # Solution algorithm from the chosen solver
'initial_strategy': 'target-coverage', # Initialization strategy for the fluence vector
'initial_fluence_vector': None, # User-defined initial fluence vector (only for 'warm-start')
'lower_variable_bounds': 0, # Lower bounds on the decision variables
'upper_variable_bounds': None, # Upper bounds on the decision variables
'max_iter': 500, # Maximum number of iterations for the solvers to converge
'tolerance': 0.001 # Precision goal for the objective function value
}
Yeah, this was a tough piece of work! If you have managed to complete the optimization dictionary, feel free to reward yourself with a cup of tea or coffee, maybe a small snack, and a relaxing short break before moving on ☕
Setting up the evaluation dictionary
It is not actually necessary to set up the evaluation dictionary if you are happy with the default values. However, we will initialize it for reasons of completeness. First, we select the DVH type ‘cumulative’ and request its evaluation at 1000 (evenly-spaced) points. With the parameters ‘reference_volume’ and ‘reference_dose’, we let the package calculate dose and volume quantiles at certain levels. By inserting an empty list for ‘reference_dose’, the levels are automatically determined. The last two parameters, ‘display_segments’ and ‘display_metrics’, can be used to filter the names of the segments and metrics to be displayed later in the treatment plan visualization. We also specify empty lists here to not exclude any segment or metric.
evaluation = {
'dvh_type': 'cumulative', # Type of DVH to be calculated
'number_of_points': 1000, # Number of (evenly-spaced) points for which to evaluate the DVH
'reference_volume': [2, 5, 50, 95, 98], # Reference volumes for which to calculate the inverse DVH values
'reference_dose': [], # Reference dose values for which to calculate the DVH values
'display_segments': [], # Names of the segmented structures to be displayed
'display_metrics': [] # Names of the plan evaluation metrics to be displayed
}
Congratulations, you have successfully set up all parameter dictionaries 🎉
Initializing the base class
Now let’s finally put everything together into a complete TreatmentPlan instance.
tp = TreatmentPlan(configuration, optimization, evaluation)
Treatment plan workflow
In this section, we describe the standard workflow in which the generated treatment plan instance comes into play. Our package equips the instance with one method for each work step, which can be called parameter-free.
Configuring the plan
First, a successfully initialized treatment plan needs to be configured. By calling the configure method, the information from the configuration dictionary is transferred to internal instances of the configuration classes, which perform functional (logging, data management) and I/O tasks (processing of imaging data, preparation of data dictionaries). Note that a plan must be configured before it can be optimized.
tp.configure()
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Initializing logger ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - You are running Python version 3.10.14 ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Initializing datahub ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Initializing patient loader ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Initializing plan generator ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Initializing dose information generator ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Importing CT and segmentation data from MATLAB file ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Generating plan configuration for photon treatment ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Generating dose information for photon treatment ...
2024-05-14 20:59:29 - pyanno4rt - TG-119-example - INFO - Adding dose-influence matrix from MATLAB file ...
Optimizing the plan
Afterwards, the treatment plan is ready for optimization. We call the optimize method, which generates the internal optimization classes, passes the optimization parameters from the dictionary, and at the end triggers the solver run. If machine learning model-based components are used, the model fitting would also take place here. In our example, no such components exist, which means that the optimization process starts immediately. Note that a plan must be optimized before it can be evaluated.
tp.optimize()
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Initializing fluence optimizer ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Removing segment overlaps ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Resizing segments from CT to dose grid ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Setting the optimization components ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Setting objective 'Squared Overdosing' for ['Core'] ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Setting objective 'Squared Deviation' for ['OuterTarget'] ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Setting objective 'Squared Overdosing' for ['BODY'] ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Adjusting dose parameters for fractionation ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Initializing dose projection ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Initializing weighted-sum optimization method ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Initializing fluence initializer ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Initializing fluence vector with respect to target coverage ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Initializing SciPy solver with L-BFGS-B algorithm ...
2024-05-14 20:59:31 - pyanno4rt - TG-119-example - INFO - Solving optimization problem ...
2024-05-14 20:59:32 - pyanno4rt - TG-119-example - INFO - At iterate 0: f=144.1014
2024-05-14 20:59:33 - pyanno4rt - TG-119-example - INFO - At iterate 1: f=119.7357
2024-05-14 20:59:33 - pyanno4rt - TG-119-example - INFO - At iterate 2: f=85.2347
2024-05-14 20:59:33 - pyanno4rt - TG-119-example - INFO - At iterate 3: f=30.4996
2024-05-14 20:59:33 - pyanno4rt - TG-119-example - INFO - At iterate 4: f=14.5176
2024-05-14 20:59:33 - pyanno4rt - TG-119-example - INFO - At iterate 5: f=12.4683
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 6: f=10.9733
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 7: f=10.3196
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 8: f=9.684
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 9: f=9.4383
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 10: f=8.9254
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 11: f=8.3518
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 12: f=8.0128
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 13: f=7.5543
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 14: f=7.3486
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 15: f=6.9583
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 16: f=6.5906
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 17: f=6.3163
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 18: f=6.1272
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 19: f=6.0453
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 20: f=5.9811
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 21: f=5.8058
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 22: f=5.7397
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 23: f=5.6837
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 24: f=5.6219
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 25: f=5.5861
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 26: f=5.5432
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 27: f=5.5025
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 28: f=5.4791
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 29: f=5.4677
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 30: f=5.4307
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 31: f=5.4234
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 32: f=5.4037
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 33: f=5.3914
2024-05-14 20:59:34 - pyanno4rt - TG-119-example - INFO - At iterate 34: f=5.3791
2024-05-14 20:59:35 - pyanno4rt - TG-119-example - INFO - At iterate 35: f=5.3647
2024-05-14 20:59:35 - pyanno4rt - TG-119-example - INFO - At iterate 36: f=5.3593
2024-05-14 20:59:35 - pyanno4rt - TG-119-example - INFO - At iterate 37: f=5.3487
2024-05-14 20:59:35 - pyanno4rt - TG-119-example - INFO - At iterate 38: f=5.3436
2024-05-14 20:59:35 - pyanno4rt - TG-119-example - INFO - Computing 3D dose cube from optimized fluence vector ...
2024-05-14 20:59:35 - pyanno4rt - TG-119-example - INFO - Fluence optimizer took 3.62 seconds (3.26 seconds for problem solving) ...
Evaluating the plan
The penultimate step usually is the evaluation of the treatment plan, and following the previous logic, we have added an evaluate method for this purpose. Internally, this creates objects from the DVH and dosimetrics class, which take the parameters of the evaluation dictionary and trigger the respective evaluation processes.
tp.evaluate()
2024-05-14 20:59:37 - pyanno4rt - TG-119-example - INFO - Initializing DVH evaluator ...
2024-05-14 20:59:37 - pyanno4rt - TG-119-example - INFO - Initializing dosimetrics evaluator ...
2024-05-14 20:59:37 - pyanno4rt - TG-119-example - INFO - Evaluating cumulative DVH with 1000 points for all segments ...
2024-05-14 20:59:37 - pyanno4rt - TG-119-example - INFO - Evaluating dosimetrics for all segments ...
Visualizing the plan
We are now at the end of the standard workflow, and of course we would like to conclude by analyzing the results of the treatment plan optimization and evaluation both qualitatively and quantitatively. Our package features a visual analysis tool that provides three sets of visualizations: optimization problem analysis, data-driven model review, and treatment plan evaluation. By clicking on the activated buttons, you can open the plot windows. The visual analysis tool can easily be launched with the visualize method.
tp.visualize()
2024-05-14 20:59:37 - pyanno4rt - TG-119-example - INFO - Initializing visualizer ...
QPixmap::scaled: Pixmap is a null pixmap
2024-05-14 20:59:38 - pyanno4rt - TG-119-example - INFO - Launching visualizer ...
2024-05-14 20:59:39 - pyanno4rt - TG-119-example - INFO - Opening CT/dose slice plot ...
QPixmap::scaled: Pixmap is a null pixmap
2024-05-14 21:00:11 - pyanno4rt - TG-119-example - INFO - Closing visualizer ...
Ideally, you should now see the window below.
(By the way, the top image in this notebook has been extracted from the CT/Dose slice plot 😉)
Shortcut: composing the plan
Many times you will just run all four of the above methods in sequence. To make this a little more convenient, the treatment plan can also be “composed” in a single step, using the appropriately named compose method (and yeah, we love music ❤️).
tp.compose()
Updating parameter values
One last class functionality is the updating of parameter values with the update method. This comes in handy because each of the configure, optimize and evaluate methods is based on the corresponding parameter dictionary, so that, for example, the evaluate method can be called again after updating an evaluation parameter without repeating all the initialization and prior workflow steps.
The update method takes a dictionary with key-value pairs as input, where the former are from the parameter dictionaries, and the latter are the new parameter values. We do not want to change the plan at this point, so we will just overwrite the modality and the DVH type with the previous values for illustration purposes.
tp.update({
'modality': 'photon',
'dvh_type': 'cumulative'
})
Saving and loading treatment plans
Treatment plans generated within our package can be saved as a snapshot folder and loaded from there as a copycat. You can import the corresponding functions from the tools subpackage.
from pyanno4rt.tools import copycat, snapshot
A snapshot automatically includes a JSON file with the parameter dictionaries, a compiled log file, and, if machine learning model-based components are used, subfolders with model configuration files. Optionally, you can specify whether to add the imaging data, the dose-influence matrix, and the model training data (this allows sharing an instance of TreatmentPlan with all input data). The name of the snapshot folder is specified by the treatment plan label from the configuration dictionary.
Assuming the snapshot is to be saved in the current path, the line below would create the minimum-sized version of a snapshot folder.
snapshot(instance=tp, path='./', include_patient_data=False, include_dose_matrix=False, include_model_data=False)
Conversely, a snapshot that has been saved can be loaded back into a Python variable by calling the copycat function with the base class and the folder path.
tp_copy = copycat(base_class=TreatmentPlan, path='./TG-119-example/')
UI-based interface
Our package can also be accessed from a graphical user interface (GUI) if you prefer this option. There are many good reasons for using the GUI:
- Support with the initialization: the GUI window provides the parameter names in a structured form and at the same time already defines parts of the parameter space.
- Handling of large parameter spaces: the parameter space of a treatment plan can very quickly become high-dimensional, especially with many components (and even more with machine learning model-based components), making the code-based generation of the dictionaries more complex than the UI-based generation.
- Faster generation of treatment plans and cross-instance comparison: due to the first two key points, among others, the GUI allows faster treatment plan initialization and workflow execution, and multiple plans can be saved and switched quickly for comparison.
And, of course, a GUI may also simply look good 😎
GUI initialization
So, how can the GUI be called? Instead of initializing the TreatmentPlan class, we create an object of the GraphicalUserInterface class.
gui = GraphicalUserInterface()
GUI opening
Then, you can open the GUI window directly using the launch method.
gui.launch()
Below you can see the main window of the GUI that should now appear.
Without going into detail, the most important widgets shall be described here:
- Menu bar (upper row): load/save a treatment plan, drop an existing treatment plan instance, instance selector, settings/info windows, exit the GUI
- Composer (left column): tabs for setting the configuration, optimization and evaluation parameters
- Workflow (middle column): action buttons for the workflow steps, toolbox for accessing generated data (e.g. log files)
- Viewer (right column): Axial CT/dose preview, interactive DVH
Alternatively, you can launch the GUI directly with an instance.
gui.launch(tp)
Fetching treatment plans from the GUI
The GUI has an internal dictionary in which objects of the TreatmentPlan class generated from the interface are stored. These can also be retrieved after closing the GUI using the fetch method.
tps_gui = gui.fetch()
Outro
We very much hope that this little example illustrates the basic usage of our package for treatment plan optimization. If you have any questions or suggestions, or in the (hopefully unlikely) event that something does not work, please take a look at the “Help and support” section and drop us a line. We would also be happy if you leave a positive comment or recommend our work to others.
Thank you for using pyanno4rt 😊