MRI-only proton therapy: How to implement calculating Energy of bixels for stf-struct on MRI instead of CT?

[lauratsu]

Dear matRad experts,

I am currently researching on MRI-only proton treatment therapy.
I would like to be able to calculate the energies of the bixels, that are calculated when running "matRad_generateStf(ct,cst,pln)", on MRI instead of CT. So far I don't have an exact plan on how to do it, but I would have to understand how it works for CT first anyway. And for that I need some help and clarification on how the code works, because I am not an expert when it comes to coding.

I am looking for the codeline, where the required energy for a CT beampath is getting calculated. So in general I am looking for the function E = f(CT_beampath) so that i could theoretically exchange it with E = f(MRI_beampath).

I looked into the code and tried to understand how the energy is getting calculated. So far I understand it like this. Please correct me if I am incorrect.

1. In "matRad_generateStf(ct,cst,pln)" -> "generator = matRad_StfGeneratorBase.getGeneratorFromPln(pln);"

Does this select an engine, which is responsible for where a bixel is positioned and what energy it has?
Here is the code I am using.

When I run this code it will select "STFGeneratorParticleIMPT" as engine. But theoretically I would need to create a new engine for my MRI-only workflow right? Because the energy would be calculated differently.
How would I specify what engine I would like to use when running my code. I thought it would get the information from the pln-struct, but when I looked at https://github.com/e0404/matRad/wiki/The-pln-struct I could't find a paramter that specifies the engine.

2. Then next in "matRad_generateStf(ct,cst,pln)" -> "stf = generator.generate(ct,cst);" -> "this.initialize();"

Does it convert the CT values to WET values with HLUT. I would probably have to work around the HLUT for my MRI-only workflow.

3. Then there is "this.createPatientGeometry();"

Does this find the position/voxel of the target volume and afterwards calculates the water equivalent distance to reach that position?

4. Then there is "stf = this.generateSourceGeometry();" -> "beam = this.initRays(beam);"

Does it give us the positions and directions of the rays?

5. Next "beam = this.setBeamletEnergies(beam);" seems to be the interesting part

There is the line "newEnergies = this.availableEnergies(this.availablePeakPos>=targetEntry(k)&this.availablePeakPos<=targetExit(k));" which actually gets the energies. These availableEnergies come from the machine. Does it compare the available Energies from the machine with the energies that were calculated in "this.createPatientGeometry();" and selects the energies that are needed?

It seems that everything is based on the conversion of CT values to WET values. So I would first need to find a solution on how to convert MRI values to WET or rather completely work around HLUT. Does someone has a suggestion on how to achieve that?
I would be happy about any answer to my questions. Any kind of help is appreciated.
Thank you :)

[wahln]

Hi @lauratsu
This is an interesting project. You got quite far "down the rabbit hole", here's a few initial comments. I hope to have some time to look into it in more detail soon, but currently we are all a bit swamped here...

First an explanation about the philosophy behind the stf generation and "StfGenerators" (stf stands for "steering file / information", and is thus describing the irradiation geometry).

• matRad_generateStf(...) is a helper function to select the appropriate generator to generate the beam information. For particles, this prioritizes, as you say, the IMPT Generator for fully-intesity modulated particle therapy. I think currently the only alternative is the "single beamlet generator". There are some helper functions to see which generators are available for your current machine selection, but I assume you want to do IMPT on MRI-only.
• The generators are build in an object-oriented fashion to group reusable code (e.g., photons and ions both need a beamlet grid, and much of the geometry calculations are similar). So navigating through the functions can be a bit complex, because one function may call the implementation from the superclass first until it executes its specific code. This is also the reason why the generation is modular (initialize, createPatientGeometry, etc.).
• The best way to check what a function does is though to "walk up" the hierarchy from StfGeneratorParticleIMPT -> StfGeneratorParticleRayBixelAbstract -> StfGeneratorExternalRayBixelAbstract -> StfGeneratorBase

Now to the functions:

• initialize(...) does mostly do some bookkeeping, like making sure the isocenter is given/obtained as well as some other variables. In StfGeneratorParticleRayBixelAbstract, it adds the computation of the relative stopping power from the look-up table, as this is always needed for particle therapy and energy selection
• createPatientGeometry(...) creates coordinate representations and filters target volumes etc.
• createSourceGeometry(...) creates the beamlet grid and aligns spots etc. "StfGeneratorExternalRayBixelAbstract" provides helper functions for this (i.e. also the initRays you discovered, which sets up the lateral spot grid and is also used for photon beamlets, for example). It also declares a setBeamletEnergies which needs to be finalized in derived classes (e.g. the IMPT generator).

The bemlet energies are selected based on a raytracing step implemented in the setBeamletEnergies implementation in StfGeneratorParticleIMPT. For this, the rSPs are needed which come from calling matRad_calcWaterEqD(...) in the initialize function of StfGeneratorParticleRayBixelAbstract (see above).
Basically, each direction is traced and integrated. Concurrently we check where these traced rays enter and exit targets. then we filter the energies in-between entries and exits, giving us the desired machine energies.

I think as a start you should tinker with this step of computing the rSPs in the latter function. The function tries to find a suitable HU->rSP look-up table and computes a ct.cube from the given ct.cubeHU, otherwise works with the default. Some suggestions I have for this:

• In the dose calculation, we allow to bypass this step and provide a ct.cube already. It might make sense to do extend the generator to allow this as well, which means you could create a synthetic ct.cube from your MRI as you wish. For example, you could start experimenting with simple density overrides. Later this could be better integrated.
• You could also create a look-up table to convert intesities directly into rSP. Don't know what the state of research is on that, but I'm skeptical.
• You could create a full "synthetic CT" and not change anything in the Generators and keep the default look-up table used by calcWaterEqD

So, in the end, you might get away with minimal modifications to the Stf-Generators. If you want to include an active magnetic field during irradiation, this becomes a different beast...

Let us know if you have further questions.