-
Couldn't load subscription status.
- Fork 61
Inverse Planning Support and pyRadPlan Integration in SlicerRT's External Beam Planning Module #275
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: master
Are you sure you want to change the base?
Conversation
|
Congratulations for this awesome work!!! I have one question:
Shouldn't there be a separate table and "C++ struct" for constraints and another for objectives? |
|
Yes, thank you for pointing this out. I forgot to mention this in the "Future Works". Currently pyRadPlan is not handling constraints, so we set the focus on the objectives for now. But yes, there should be a similar structure/table for constraints in the future. |
|
Thanks for clarifying :) One observation, we could maybe consider to have a dictionary of predefined standard optimization functions, stored in a dict on the C++ side, and shown as a dropdown menu in the GUI, following a bit https://youtu.be/FIWYQRacRDI?feature=shared&t=326 Other suggestions for "future work", would be to have a "Robust" checkbox (even if not support at this stage, it would be good to leave room for it in the GUI); and to include a optimization-progress feedback in the GUI / progress bar or so. Another comment: Weight instead of Penalty might be more physician-friendly as label. |
|
Just letting everyone know that I'm back from vacation, and I'll look at this ASAP! This is very exciting, thank you so much for the hard work! |
|
This 20-proton-patients dataset might be helpful for stress-testing this PR: https://consigna.uv.es/?recoge:ca:3245_9718_9829_7607 |
|
Hi everyone, sorry for the late reply. I also have more time now and will be available to actively work on the pull request in the upcoming weeks. To comment on the optimization/objective functions: I'll add the robust checkbox and try including an optimization progress bar. And yes, I would greatly appreciate feedback on any naming conventions. I mostly oriented myself on pyRadPlan's conventions here, which might not always be transferable. |
|
Thanks a lot! Regarding names, etc. I leave here a link https://consigna.uv.es/?recoge:ca:7746_1317_7296_5505 with screenshots of GUIs to get inspiration from. Maybe it's good if we can make naming match the most commonly used conventions. Regarding checkbox: would it also make sense to create a "Mean" checkbox? Or do we prefer having separate "MaxDose" and "MeanMaxDose" objectives?
Or maybe instead, one has to choose first the optimizer, and then the content of the dropdown menu is dynamically adjusted? |
|
@LiBu981 I did a quick (2-hour) review round of this PR and it looks quite good! Thank you very much for the work you've done on this! I only have one concrete comment currently: I'd rename the Also a question: in your branch does forward planning still work as before? I'll continue the review in the following days. After that maybe we could have a meeting, even with @ferdymercury or @wahln to discuss some details, facilitating the review process? |
|
Maybe a (not so) short comment from my side on the design choice of letting the optimizer define the objectives that can be used itself. We first thought about providing a general objective & constraint interface (i.e. similar to the dose engines or the optimizer, objectives and, in the future, constraints can be defined in Python or C++). However, we were concerned about the computational overhead of this, as then data would need to be transfered during each iteration of the optimization if you mix a python optimizer with a C++ objective. Sometimes RT optimizers also have dedicated implementations (they need 2nd order derivatives, or use some tricks in evaluating linear constraints, and so on) that would be difficult to represent in a general interface. Also, what if you have a very specialized optimization routine doing something really different, like working with dose rate, biological modeling, or LET on top - this would also already create a headache. At some point you would end up with defining so many flags / weird class hierarchy situtations to make sure an optimizer understands capabilities of the objectives and so on... So in the end, we thought it would be better if the interfaced "optimizer" provides a list of what objectives are available and what configuration parameters they have, and SlicerRT just handles it on this level for display and returning the information to the optimizer. This also makes it very simple to interface a custom planning routine that works with some crazy objective prototypes. Every idea is welcome here on how to have the most user- and developer-friendly situation here. |
|
@LiBu981 I remember the weird bug of SlicerRT crashing when the plan was created too quickly after starting the EBP Module. Is this still a thing? |
Yes, it still occasionally crashes on creating a new plan. |
Thanks Niklas for the feedback! I understand your point and motivation. The drawback (from a user-friendly perspective) of delegating everything to the optimizer would be that, if I now want to change from EngineA or EngineB to test which one gives faster or better results, one has to redefine all objectives and constraints. This makes sense if EngineA does LET optimization and EngineB does dose rate optimization. However, for more 'common' examples, maybe the user just defined MinDose and MaxDose objectives, since they are not even aware of those other crazy objectives. Now let's imagine that both EngineA and EngineB implement exactly the same number and type of objectives, but since they were developed by different persons, one calls these objectives Maybe this is also a question of whether we want this to be optimization-research-oriented or newbie-user-oriented. Tradeoff and compromises. Maybe there is a clever way of doing this, such as follows:
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I did a first round of review. Please fix the obvious ones and let's discuss anytihng that is not obvious. Thank you very much!
Beams/MRML/vtkMRMLRTBeamNode.cxx
Outdated
| this->CollimatorAngle = 0.0; | ||
| this->CouchAngle = 0.0; | ||
|
|
||
| this->DoseGridDim[0] = -1; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Where are these (dimensions and spacing) set? I cannot find it. Also, dimensions and spacing are not enough to define a volume geometry, we also need orientation. I understand that for the Dij matrix it may not be needed (manybe neither spacing, not sure), but this seems to be hanging in the air here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Currently it is saved here , if passed as input when a dose influence matrix is saved in the beamNode.
This is used for example in the pyRadPlanEngine.py whereby the dij.dose_grid is defined in the prepareRTDataset.py.
I added this to be able to resample the dose in the MockOptimizer. Maybe this can be handled differently though?
Beams/MRML/vtkMRMLRTBeamNode.h
Outdated
| DoseInfluenceMatrixType DoseInfluenceMatrix; | ||
|
|
||
| /// Dose grid dimensions (on which the dose influence matrix is defined) | ||
| double DoseGridDim[3]; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This should not be double. The get/set functions use int. I'd use unsigned int probably.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also now we can initialize values in the header like:
unsigned int DoseGridDim {5, 5, 5};
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I changed it to int now, as I initialize here with {-1, -1, -1} to check weather a DoseGridDim (same with DoseGridSpacing) was set. So unsigned int wouldn't work with this logic.
(e.g. in pyRadPlanUtils when called from the pyRadPlan optimizer )
| vtkMRMLRTObjectiveNode(const vtkMRMLRTObjectiveNode&); | ||
| void operator=(const vtkMRMLRTObjectiveNode&); | ||
|
|
||
| char* Name; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just use std::string, much easier.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I used char* so that vtkGetStringMacro(Name) works? Is there a way to do this with std::string?
(I oriented myself on the beam node.)
|
|
||
| char* Name; | ||
| /// Name of the segmentation assigned to this objective | ||
| std::string Segmentation; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Node reference, very important!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See comment :)
| qSlicerScriptedDoseEngine.h | ||
| qSlicerAbstractPlanOptimizer.h | ||
| qSlicerPlanOptimizerPluginHandler.h | ||
| qSlicerPlanOptimizerLogic.h |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It is strange to have classes called Logic in the Widgets project. I'm sure there is a good reason, let's discuss this later.
|
|
||
| // Connect plan modified events to population of the table | ||
| qvtkReconnect(d->PlanNode, planNode, vtkCommand::ModifiedEvent, this, SLOT(updateObjectivesTable())); | ||
|
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm struggling a bit to add a connection between the SegmentationNode in the Plan and in the Objectives Table, so that on modifications of the selected segmentation (e.g. name of a segment) the selectable segments in the objectives table are modified as well. (I tried looking into the qMRMLBeamsTableView.cxx for inspiration.)
Does someone have ideas?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you don't know the type of modification (for example segment has been added, deleted or renamed) then just clean a table and refill it with new information for each segment.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If position of the segment ID with segmentation node is unchanged, then just update corresponding table items according to the position of the segment within table.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In general I agree with @MichaelColonel in that in case of such custom table widgets, the simplest is to re-populate the table on relevant changes. The table will never contain more than a few dozen rows, so speed is not an issue. In general there are two ways to implement tables that sync with the MRML scene:
- Qt model-view pattern: see
qMRMLSubjectHierarchyModelandqMRMLSubjectHierarchyTreeViewfor example. It seems complicated, but it is the correct way - However, sometimes a simple table widget is enough, especially if there is not a lot of data to handle (unlike the whole MRML scene can be). See functions starting with
populatefor example. I think even the SegmentsTableView worked like that in the beginning.
Your solution seems to be halfway between the two, in that you handle addition, deletion, etc. So either simplify it to always re-populate (but it's work done twice, and more time), or do what @MichaelColonel suggests second, i.e. observe the segment modified event of the segmentation node, then look for its row in the table, and update the item text. I hope this helped a bit. If not, happy to meet Friday morning to discuss (this and other questions too).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay, thank you so much for all the feedback. I finally managed to get to this and implemented reactions to changes in the segmentationNode and the planNode and it seems to work properly, though I would appreciate a review (see commit 4215189). I hope this corresponds to your comments :) Thanks in advance!
| @@ -0,0 +1,680 @@ | |||
| /*============================================================================== | |||
| Copyright (c) Laboratory for Percutaneous Surgery (PerkLab) | |||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Feel free to change the copyright in new files to your institution
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See commit c803da7
| qCritical() << Q_FUNC_INFO << ": Segment ID is empty!"; | ||
| return; | ||
| } | ||
| vtkMRMLSegmentationNode* segmentationNode = d->PlanNode->GetSegmentationNode(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fix formatting here
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sorry, I thought I had fixed this already, but some places it were still tabs instead of spaces. Should be fixed now.
(See commit 5b22bd3)
|
Hi @LiBu981 I tried to build this branch but I got build errors. When I tried to rebase it to the latest SlicerRT I immediately got conflicts, and since you have dozens of commits and I'm not that familiar with the content and the reason of the choices you made, I'd like to kindly ask you to do the rebase on the latest SlicerRT so that the branch is up-to-date. Thank you! |
|
Is rebased now :) |
|
This should fix the crash: commit 2dc6206 I haven't had any crashes now. |
|
Thank you so much @LiBu981 ! I'll build it today and start reviewing (and understanding) it in more detail |
Hmmm the problem might actually be that your repository is a fork of https://github.com/e0404/SlicerRT rather than a fork of https://github.com/SlicerRt/SlicerRT I think the easiest would be to create a local backup, erase your repo from GH, and fork again cleanly from the SlicerRt, and then push back from local copy your branches to upstream. But I think there are other ways such as cli/cli#10956 |
|
Thanks a lot @LiBu981 I was able to push this commit now! Strange, because when we do a PR in Slicer, the people with push rights (JC, Andras, ...) can push directly to the topic branches (which are from forks). Anyway, good for now :) Thanks again |
yep. Maybe due to of the fork-of-fork structure (see above). |
|
Ah OK I didn't catch that. Thanks @ferdymercury |
| /// Structure to represent objective functions (only name and parameters) | ||
| struct ObjectiveStruct | ||
| { | ||
| std::string name; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As a Qt class, it would make sense to use QString and QMap. Is there any reason not to use those?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I thought I had problems before with some conversions here, but I've changed it to Qt types now without any issues. So not entirely sure why I didn't use it before.... (see commit 7acb67b)
| // Optimization calculation related functions | ||
| public: | ||
| /// Perform Optimization calculation for a single beam | ||
| /// \param Beam node for which the Optimization is calculated |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This line seems to be wrong. It says beam node but the argument is plan node.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for spotting! It was left over from the Dose Engine. I hope I've caught any wrong descriptions now (see commit d7e8297).
|
|
||
| // Optimization calculation related functions | ||
| public: | ||
| /// Perform Optimization calculation for a single beam |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"for a single beam" => seems otherwise (i.e. whole plan)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As above :) (see commit d7e8297)
| from Python.PyRadPlanUtils import prepareCt, prepareCst, preparePln | ||
|
|
||
| class pyRadPlanEngine(AbstractScriptedDoseEngine): | ||
| """ Mock python dose engine to test python interface for External Beam Planning |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Fix description please
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See commit 1b4bfcf
| def calculateDoseUsingEngine(self, beamNode, resultDoseVolumeNode): | ||
|
|
||
| #################################### Import pyRadPlan libraries ###################################### | ||
| import logging |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As far as I know the good practice is to import everything on the top of the file
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I think I can put everything except the pyRadPlan internal imports outside/ at the top of the file. Otherwise, we will get an error on opening Slicer before pyRadPlan is installed the first time (with the popup install request).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
See commit 1b4bfcf
|
Finishing for today. Sorry for the many comments here. I also added some in the code, they are minor fixes. One thing that I'd like to ask in general is to review the usage of smart pointers. In some cases it is not needed, for example in
Thanks a lot again! I'll finish reviewing the code probably tomorrow and will finally try it out in action. I think we are close to integration! |
|
I've reviewed the use of the vtkSmartPointer and hope I've caught any wrong uses (see commit b2a87). Though, if I saw correctly, we normally use a vtkSmartPointer for the vtkImageData before passing it to SetAndObserveImageData() (e.g. in |
| // Give default name for result node (engine can give it a more meaningful name) | ||
| std::string resultOptimizationNodeName = std::string(planNode->GetName()) + "_Optimization"; | ||
| resultOptimizationVolumeNode->SetName(resultOptimizationNodeName.c_str()); | ||
| resultOptimizationVolumeNode = vtkMRMLScalarVolumeNode::New(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This will result in a memory leak. Please use vtkNew, otherwise this reference will be dangling when the scene is cleared and the object will not be deleted from memory.
|
|
||
| // Create total dose image data from reference volume | ||
| vtkImageData* totalDoseImageData = vtkImageData::New(); | ||
| vtkSmartPointer<vtkImageData> totalDoseImageData = vtkSmartPointer<vtkImageData>::New(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Better to use vtkNew here as well.
|
I just pushed a really minor commit. I hope I'll have more time for this in the next week - I haven't been lucky lately... Thanks for your patience. Let me know if there is anything in particular I should look at. Thank you! |
|
Thank you a lot @cpinter! I think a review on the objective table's reactions to changes in the plan and segmentation node see here would be good, as I'm not 100% sure if I implemented everything correctly. There are a few minor issues I wanted to fix the coming days as well as integrate your comments, though from my side this was the biggest ToDo still open. |
…dangling reference when scene is deleted)

Summary
This pull request extends SlicerRT's External Beam Planning (EBP) module to support a more modular architecture for treatment planning by separating dose calculation and optimization for inverse planning, while focusing on integrating the treatment planning toolkit pyRadPlan.
Key Changes
Separation of Dose Calculation and Optimization
Code under
ExternalBeamPlanning\Widgets\Pythonhas been restructured intoDoseEnginesPlanOptimizersWhen "Inverse Treatment Planning?" is selected, dose calculation and plan optimization are handled as two distinct steps.
This separation allows different engines to be used for each task independently.
Inverse Planning
Each dose engine reports whether it supports inverse planning (
isInverseflag).If supported and "Inverse Treatment Planning?" is selected:
beamNode.The plan optimizer:
planNodeand selectedobjectives.beamNodesunder theplanNode(as sparse matrix in C++ or as field data in python).resultOptimizationVolumeNode.optimizePlanUsingOptimizer(planNode, objectives, resultOptimizationVolumeNode)Objectives Table
A new UI table (based on qMRMLBeamsTableView) allowing the user to define objectives for inverse planning.
For each objective (each row) the user selects:
An
objectiveNodeis created for each row, holding...... and is saved in the optimizer
For selection available objectives are defined by the optimizer and stored as
std::vector<ObjectiveStruct>, where:qSlicerSquaredDeviationObjectiveis provided and registered as an available objective for the Mock Optimizer, though it is not yet used during the optimization process.Integration of pyRadPlan
On opening the EBP module, the system checks for a pyRadPlan installation and prompts the user to install if needed.
The following Python-based components are integrated:
pyRadPlanEngine.py--> dose enginepyRadPlanPlanOptimizer.py--> plan optimizerPyRadPlanUtils.py--> constructsct,cstandplndata structures from thebeamNode,referenceVolumeandplanNode(pyRadPlan applies stable validation and serialization with pydantic)
Dose calculation (inverse planning):
stfanddijdata structures are computed fromct,cstandpln.dij) is saved in thebeamNode, along with the corresponding dose grid dimension and spacing.Optimization:
dij.cst.dij.ScalarVolumeNodes.Known Limitations / Future Works
planNodeon startup - this requires further investigation and handling.