Direct implementation - topasmc/dicom-interface GitHub Wiki

For direct class implementation, one file (rbe_1p1.hpp) needs to be created to describe the treatment machine and one existing file (treatment_session.hpp) needs to be edited.

Create a new treatment machine file, rbe_1p1.hpp and put somewhere (see below) in your system.

$ cd /<your_sw_path>/rti.git/rti/treatment_machine/rbe/ 
$ vi rbe_1p1.hpp

note: complete file is in /<your_sw_path>/rti.git/rti/treatment_machine/rbe/rbe_1p1_complete.hpp

The new treatment machine rbe_1p1 inherits rti_treatment_machine_ion and implements site-specific information handling in overriding methods, such as, characterize_history, characterize_beamlet, characterize_rangeshifter, and characterize_aperture. For the more detail explanation with the code, see below.

Define a class as following:

#include <rti/base/rti_treatment_machine_ion.hpp>

template <typename T>
class rbe_1p1 : public treatment_machine_ion<T> 

And you need to implement the following methods:

• rbe_1p1 : in this default constructor, you can define SAD.

  rbe_1p1()
  {
      treatment_machine_ion<T>::SAD_[0] = 1800.0; 
      treatment_machine_ion<T>::SAD_[1] = 2000.0;
  }

• characterize_history is a method where you put the conversion rule from the spot's meterset (MU or NP) to number of histories to be simulated. scale is a factor for particles per history. This factor is important because the real number of protons to be simulated is too high to complete in a reasonable amount of time. This example assumed that the meterset is the number of particles (NP).

    virtual size_t
    characterize_history(
        const rti::beam_module_ion::spot& s, 
        float scale)
    {
        return s.meterset / scale;
    }

• characterize_beamlet is a method to convert spot properties to a PDF distribution (in RTI, we call it beamlet) for energy and phase-space variables. For detailed information on these distribution functions, see Beam source->Distributions

  virtual
  rti::beamlet<T>
  characterize_beamlet(const rti::beam_module_ion::spot& s)
  {
  
      // 1. energy distribution 
      // constant energy
      auto energy = new rti::const_1d<T>({s.e},{0});
  
      // 2. X-Y position at z
      // DIR vector at beam(x,y,z)
      rti::vec3<T> iso(s.x, s.y, 0);
      rti::vec3<T> beam = this->beam_starting_position(iso,rti::treatment_machine_ion<T>::source_to_isocenter_mm_);
      rti::vec3<T> dir = iso - beam;
      dir.normalize(); 
      
      // 3. Complete fluence distribution
      // this samples x, x', y, y', z, z'
      std::array<T,6> spot_mean = {beam.x, beam.y, beam.z, dir.x, dir.y, dir.z};
      std::array<T,6> spot_sigm = {s.fwhm_x/T(2.355), s.fwhm_y/T(2.355), 0.0, 0.0, 0.0, 0.0};
      std::array<T,2> corr      = {0.0, 0.0};
      auto fluence= new rti::phsp_6d<T>(spot_mean, spot_sigm, corr);
  
      return rti::beamlet<T>(energy, fluence);
  }

• characterize_rangeshifter is a method to configure a rangeshifter from DICOM information. This was very difficult to generalize the conversion from DICOM information to a rangeshifter's dimension, shape, and a position because some hospitals imply the thickness in range shifter's ID, which is string, e.g., "RS 10mm" or "RS10" instead of using the thickness in WET DICOM tag. In addition, the rangeshifter's positions can be defined by offset from snout position instead of tray distance DICOM tag. So this method offers users to configure a rangeshifter based on their DICOM specification.

  rti::rangeshifter*
  characterize_rangeshifter(
      const rti::dataset* ds,
      rti::modality_type m)
  {
      //0. Get a Rangeshifter sequence for given modality.
      auto seq_tags = &rti::seqtags_per_modality.at(m);
  
      //1. Get a rangeshifter dataset from a rangeshifter sequence.
      auto  rs_ds = (*ds)( seq_tags->at("rs")) ; 
      assert(rs_ds.size() >=1);
  
      //2. Snout position from control point 0 of ControlPointSequence 0,
      std::vector<float> ftmp;
      auto layer0    = (*ds)(seq_tags->at("ctrl"))[0]; //layer0 for snout position
      layer0->get_values( "SnoutPosition", ftmp);
      
      rti::vec3<float>    lxyz(300.0, 300.0, 0.0) ;  //x-y sizes
      rti::vec3<float>    pxyz(0.0, 0.0, ftmp[0]) ;  //position
      rti::mat3x3<float>  rxyz(0.0, 0.0, 0.0) ;      //rotation
  
      //3. check RangeshifterID to detiner shifter thickness.
      for(auto r : rs_ds){
          std::vector<std::string> rs_id(0);
          r->get_values("RangeShifterID" , rs_id);
          
          if ( !rs_id[0].compare("RS10")){
              lxyz.z += 10.0;
          }else if ( !rs_id[0].compare("RS20")){
              lxyz.z += 20.0;
          }
          assert(lxyz.z>0);
      }
  
      pxyz.z -= lxyz.z ;
      std::cout<<"Range shifter thickness: " << lxyz.z 
      << " (mm) and position: " << pxyz.z <<" (mm)" << std::endl;
  
      return new rti::rangeshifter(lxyz, pxyz, rxyz);
      
  }

• characterize_aperture is a method to configure an aperture from DICOM information. In this method, users can configure site-specific aperture. Multiple-holes aperture is supported.

  rti::aperture*
  characterize_aperture(
      const rti::dataset* ds,
      rti::modality_type m
  ){
      //1. aperture opening points. xypts is a vector of holes.
      auto xypts = this->characterize_aperture_opening(ds,m);
  
      //2. get a block sequence from a given modality
      auto seq_tags = &rti::seqtags_per_modality.at(m);
      auto  blk_ds = (*ds)( seq_tags->at("blk")) ; 
      assert(blk_ds.size() >=1);
  
      std::vector<float> ftmp;
  
      blk_ds[0]->get_values( "BlockThickness", ftmp); //thickness was given
      rti::vec3<float> lxyz(400.0, 400.0, ftmp[0]);   //x-y sizes are not given from DICOM. so nominal values are used.
  
      blk_ds[0]->get_values("IsocenterToBlockTrayDistance", ftmp); //distance from the isocenter
      rti::vec3<float> pxyz(0.0, 0.0, ftmp[0]);
  
      rti::mat3x3<float> rxyz(0.0, 0.0, 0.0);
  
      return new rti::aperture(xypts, lxyz, pxyz, rxyz);
  }

This is a step to register a new treatment machine rbe_1p1 in treatment_session, especially in create_machine. So when a RTIP or RTIBTR dicom has its machine name as "RBE" and treatment machine name as "1.1", treatement_session instantiates rbe_1p1 and calls rbe_1p1.create_beamline and rbe_1p1.create_source upon requests from MC engine.

See below

#include <rti/treatment_machine/rbe/rbe_1p1.hpp>

bool
create_machine(
    std::string machine_name, 
    std::string mc_code)
{

    if( !machine_name.compare("rbe:1.1")){ //lower cases only.
        tx_machine_ = new rti::rbe::rbe_1p1<T>;
        return true;
    }
    return false;   
}