topasmc / dicom-interface

DICOM RT Interface for TOPAS
MIT License
23 stars 9 forks source link

How is it implemented? #9

Closed nicverbeek closed 8 months ago

nicverbeek commented 4 years ago

Hello,

first of all, many thanks for your outstanding work on the dicom interface. I would love to use it to integrate our clinical beamline at the West German Proton Therapy Center Essen. I do have a few questions, though:

  1. How do you interpolate between energies? I would run MC simulations in 5 MeV steps to adjust the TOPAS energies to our actual facility ranges measured in water. To do this I would need to know if this step size is sufficient in terms of accuracy, or if it can even be chosen rougher. So how does the interpolation between the energies and the energy spread work?

  2. the same for the Fermi-Eyges source parameters (PHSP source) how does the interpolation between the Fermi-Eyges source parameters work?

3.how is the source specified? Are the beam parameters specified in the isocenter? If the beam parameters are specified in the isocenter, how are they transported back? Is the source transported back in vacuum?

  1. are different Range Shifter positions supported?

  2. how is the scanning realized? Is it simulated by magnetic fields in TOPAS?

  3. what is taken from the planning system? I define the source with the energies and how is the conversion between the energies I define and the energies of the planning system done? To be more precise: In the treatment room we have a nominal energy of 160 MeV with a corresponding range. To achieve this range in RayStation, 160.25 MeV is applied in the RS MC software. To achieve this range in TOPAS, I have to apply 159.93 MeV in TOPAS.

How is this handled in your software?

  1. why is there a Stopping Power file, do you not use the materials stored in TOPAS? How does this work in general with the material definition, Can I define my Range Shifter material in the TOPAS file where I access RTI?

  2. i had read in the Issues section in github that there is a source to be programmed at the moment, is that still correct? Because the files which are defined for an own source 'rbe_1p1.hpp' and edited 'treatment_session.hpp' are already there? Is not only the machine description file 'rti/treatment_machines/pbs/rbe_1p1.table' missing? I would definitely want to use a bivariate source, which is why it is not clear to me what I have to create from scratch. If this is now done via the .table file, where do I put it and how do I call it in TOPAS? Because I cannot find the location, where the .table file is called in rbe_1p1.hpp.

I'm sorry that these are so many questions, but I would really love to use your DICOM interface to build the treatment engine for our facility. I really appreciate your work, and would be happy if you could support me.

Best,

Nico

ferdymercury commented 4 years ago

Hi, here answers to some of your questions:

  1. 5 MeV steps sounds reasonable. The user can decide what interpolator to use. This is done when you define your treatment machine, in the characterize_beamlet() function. There, you can call the linear interpolator see https://github.com/topasmc/dicom-interface/blob/f80227e93ab8a51cb1f18984be88c58eeac5db6b/rti/base/rti_utils.hpp#L45, or write a custom one in the code. One example of use here: https://github.com/topasmc/dicom-interface/blob/f80227e93ab8a51cb1f18984be88c58eeac5db6b/rti/base/rti_stopping_power.hpp#L27 If you want, you can also use a 2D Table for 2D linear interpolation, see https://github.com/topasmc/dicom-interface/blob/f80227e93ab8a51cb1f18984be88c58eeac5db6b/rti/base/rti_utils.hpp#L110

One side comment, you might also need to use a correction function for the MU to number of proton conversion. If you use RayStation, this conversion is normally easily off by 10% compared to TOPAS.

  1. Yes, beam lateral spread are at isocenter. The beam starts at a predefined distance from the iso, magnets are not included. You specify this value via a TOPAS parameter. The default value is 39 cm, but can be overriden. There is a potential future option to read this value automatically from the SAD of the DICOM RTplan file, see https://github.com/topasmc/dicom-interface/issues/5

  2. Range shifter position is read from the DICOM RTplan

  3. Beam is pointed at X-Y scanning position on isocenter. No magnetic scanning is simulated AFAIK.

  4. We usually define the conversion of nominal energy to TOPAS energy with a table. The nominal energy is read from the RTplan, and then you apply a desired table conversion. Alternatively you can define it with a analytical equation. You specify then in the characterize_beamlet() function what you want to use. For example:

    int nEnergies = 3;
    int nRows = 2;
    float eneConversion[nRows][nEnergies] =
        {
            // Nominal E in MeV
            { 100, 105, 110},
            // Corrected TOPAS Energy in MeV
            {100.57, 105.66, 110.62}
        };
    const float pEne = rti::TableInterpolation(&(eneConversion[0][0]), &(eneConversion[1][0]), s.e, nEnergies);
    auto energy = new rti::norm_1d<T>({pEne},{energySpread});

    If you have more calibrated things like energy spread, then add one more row to this table.

The other questions I leave for the author of this software. Basically, the rbe1p1 is just one example. You have to copy this file and create your own version of it, and implement your tables. There is an idea to do this just based on 'ini files', but is 'future work' https://github.com/topasmc/dicom-interface/issues/10.

ferdymercury commented 4 years ago

See also https://github.com/topasmc/dicom-interface/wiki/How-to-create-your-treatment-machine

jungwookshin commented 4 years ago

@ferdymercury Thank you very much!

ferdymercury commented 4 years ago

you are welcome. @jungwookshin can you answer to questions 2, 7 and 8?

jungwookshin commented 4 years ago

Reply to 2: currently I didn't support direct twiss parameters (alpha, beta, gamma and emittance). Instead, what you can set up for a phase-space source are x, x', y, y' and their correlations, x-x' and y-y'. This is equivalent of "BiGaussian" in topas, https://topas.readthedocs.io/en/latest/parameters/source/emittance.html

Reply to 7: Having stopping power is not critical for the RTI project. I used to correct beam energy at some distance from the isocenter. For example, if you calibrated your beam's energy at iso-center and then want to know the beam energy at 50 cm or something, you may need to adjust beam energy. Using relative stopping power of water and air, the energy differences can be calculated. I used this but it's totally up to you... Or you can use simple WET_air value, e.g., 0.1 mm wet for 10 cm air.

Reply to 8: The file was missing. Thank you for picking up. As I have managed two repositories for this project, one for private and the other for public (here) as I need to keep my institution's data private. I will merge the source code.

@nicverbeek hope answers from @ferdymercury and myself help you to understand. Let us know. Then, I would like to close this thread or we can create FAQ page in wiki based on your questions.

ferdymercury commented 1 year ago

@nicverbeek are your questions solved?