Geometry - topasmc/dicom-interface GitHub Wiki

Vector and matrix

RTI's own vector and matrix classes were implemented.

vec2<T> v2;
vec3<T> v3;
vec4<T> v4;
mat3x3<T>  m3;
mat4x4<T>  m4;
v2.x; v2.y
v3.x; v3.y; v3.z
v4.x; v4.y; v4.z; v4.s //s is scale

Coordinate transform

Coordinate transform is necessary when the phase-space variables (\x,y,z,x',y',z') from a beamlet coordinate system are transformed to global (patient or IEC) coordinate system. For this purpose, the rti::coordiante_transform class is introduced and consists of one translational shift and four rotations, i.e., collimator, gantry, couch, and iec2dicom. See below.

template <typename T>
class coordinate_transform{
public:
    /// Final rotation matrix and translation vector
    rti::mat3x3<T>  rotation    ;
    rti::vec3<T>    translation ;

    /// rotation matrices 
    rti::mat3x3<T>  collimator      ; ///< Rotation due to collimator angle
    rti::mat3x3<T>  gantry          ; ///< Rotation due to gantry angle
    rti::mat3x3<T>  patient_support ; ///< Rotation due to couch
    rti::mat3x3<T>  iec2dicom       ; ///-90 deg (iec2dicom), 90 deg (dicom2iec)

The final rotation, to which phase-space variables is multiplied, is then

rotation = iec2dicom * patient_support *  gantry * collimator    ;

For how to pass this coordinate system to beamlet, see Beamsource.

Following example shows how to create a new coordinate system from TOPAS parameter file. RTI takes radian as angle unit.

std::array<float, 4> rotations;
rotations[0] = static_cast<float>(fPm->GetDoubleParameter(GetFullParmName("RotCollimator"), "Angle") / rad);
rotations[1] = static_cast<float>(fPm->GetDoubleParameter(GetFullParmName("RotGantry") "Angle") / rad);
rotations[2] = static_cast<float>(fPm->GetDoubleParameter(GetFullParmName("RotPatientSupport"), "Angle") / rad);
rotations[3] = static_cast<float>(fPm->GetDoubleParameter(GetFullParmName("RotIEC2DICOM"), "Angle") / rad);

rti::vec3<float> p_xyz; //iso-center shift
rt_coordinate_topas_ = rti::coordinate_transform<float>(rotations, p_xyz);

CT, rtdose, and DVF

RTI also has box type geometries. A template class, rect3d, is implemented and then inherited by other similar box objects, such as RTDOSE, CT, and DVF (deformation vector field).

    template<typename T, typename R> class rect3d

    template<typename R> class ct : public rect3d<int16_t, R>

    template<typename R> class rtdose : public rect3d<float,R>

    template<typename S, typename R> class dvf : public rect3d<rti::vec3<S>, R>

Using these objects, we can read DICOM data and output into binary format.

note: reading mha file is supported for DVF data only (2019 July).

    //Read CT dir
    rti::ct<float> myct("ct_dir");
    myct.load_data(); //Load CT contents

    //Read RTDOSE file
    rti::rtdose<float> mydose("RTDOSE_FILE.dcm");
    mydose.load_data(); //Load RTDOSE contents

    //Read DVF (deformation vector field) data
    rti::dvf<float, float> mydvf("DVF_file.dcm");
    mydvf.load_data();

    //Output to binary format for checking purpose
    myct.write_data("ct.raw"); 
    mydose.write_data("dose.raw"); 
    mydvf.write_data("dvf.raw");

We can create an empty rectangle and copy the point's structure from others.

//Create empty classes
rti::rect3d<int16_t,float>  dose0;
rti::rect3d<float,float>  dose1;
rti::dvf<float,float>  dvf5to0;
rti::clone_structure(myct,  dose0);   //make dose0 grid structure same with CT
rti::clone_structure(mydose,  dose1); //make dose1 grid structure same with RTDOSE
rti::clone_structure(mydose, dvf5to0);  //make DVF grid structure same with RTDOSE

Then, we can fill values into created ones based on trilinear interpolation.

  • Interpolate CT and fill -1000 for points outside of CT grid
    int16_t fill_air = -1000; 
    rti::interpolate(myct, dose0, fill_air);
  • Interpolate RTDOSE and fill 0 for points outside of dose grid
    float fill_value = 0.0;
    rti::interpolate(mydose, dose1, fill_value);
  • Interpolate DVF and fill (0,0,0) for points outside of DVF grid
    rti::vec3<float> fill_vector(0.0,0.0,0.0);
    rti::interpolate(mydvf, dvf5to0,  fill_vector);
  • Deform dose from moving image (moving) to reference image (dose0) using DVF image
    rti::warp_linear(mydose, dose0, mydvf);

Beam limiting devices

beam limiting devices in DICOM-RT are rangeshifter, block, snout, wedges, bolus, applicator, Lateral spreading device, range modulator, compensator, and MLC for RT-ion, from http://dicom.nema.org/dicom/2013/output/chtml/part03/sect_C.8.html#sect_C.8.8.25. Currently we support only rangeshifter and block (aperture).

Rangeshifter

RTIs rangeshifter class is defined as

class rangeshifter : public geometry{
public:
    const bool        is_rectangle; ///< type of volume
    const rti::vec3<float> volume;  ///< volume dimension

    /// Constructor
    rangeshifter(
        rti::vec3<float>& v,    ///< x,y,z or r, theta, thickness
        rti::vec3<float>& p,    ///< position
        rti::mat3x3<float>& r,  ///< rotation matrix
        bool is_rect=true)

A is_rect defines the shape of a rangeshifter. If is_rect is false, the volume parameter, v.x, v.y and v.z should be used for r, theta, and thickness. Otherwise they stand length of x,y, and z.

Usually, the rangeshifter's thickness information is expressed in RangeShifterID (300a,0318). For example, RS id '33mm RS' may represent 33 mm thickness of range shifter. Besides, the thickness is expressed in Water equivalent thickenss depending on the institute. Like to aperture, there is no universial rule to parse rangeshifter id to the thickness. The position of rangeshiter can be calculated from snout position or is specified in a DICOM tag (Isocenter to range shifter distance,300a,0364). So again,treatment_machine_ion::characterize_rangeshifter method is provided to users to implement their own rule by class overriding.

Aperture (Block)

The rti's aperture class is defined as follows

class aperture : public geometry{
public:
    /// whehter aperture shape is box or cylinder type
    const bool is_rectangle ;

    /// aperture dimension, (Lx, Ly, Lz) for box or (R, H, dummy) for cylinder
    const rti::vec3<float> volume;
    
    /// X-Y opening points. Divergence is not considered.
    const std::vector<std::vector<std::array<float,2>>> block_data;
public:
    /// \param xypts a set of (x,y) points. 
    /// \param v a volume, e.g., (Lx,Ly,Lz) or (R, thickness, ignored). 
    /// \param p a center position of the aperture. 
    /// \param is_rect tells whether aperture is box shape or cylinder shape.
    aperture(
        std::vector<std::vector<std::array<float,2>>> xypts,
        rti::vec3<float>& v,
        rti::vec3<float>& p,
        rti::mat3x3<float>& r,
        bool is_rect = true)

All parameters are same with rangeshifter except xypts which is a x-y points vector for aperture holes. It is difficult to generalize all the variants of apertures too. For example, aperture thickness is usually specified by a tag (300a,0100) but some hospitals determine the thickness depending on the beam energy. The higher beam energy, the thicker the aperture. Besides, there are no tags are available in DICOM to tell the shape of aperture whether circular or rectangular. So, treatment_machine_ion::characterize_aperture method is provided to implement their own rule by class overriding.

etc

The issues with beam limiting devices are (correct me if I'm wrong):

  • Snout: the geometry of a snout will not be created for MC simulation at least IMPT simulation. Snout ID (300a,030f) will be used to determine the cross sectional size of apertures and rangeshifter.
  • Compensator: nice to have but not a high priority task.
  • Wedges: No plan to implement yet.
  • Bolus: No plan to implement yet.
  • Range modulator: No plan to implement yet.
  • MLC: No plan to implement yet.
⚠️ **GitHub.com Fallback** ⚠️