CBDose C++ Source Code

The C++ code for the CBDose implementation consists of a header file (CBDose.hpp) and a source file (CBDose.cpp). The CBDose class calculates dose from external beam radiation therapy using a corrections-based method.

Header File (CBDose.hpp)

CBDose.hpp - Click here to download
/******************************************************************************/
/*                                                                            */
/* Copyright 2016-2017 Steven Dolly                                           */
/*                                                                            */
/* Licensed under the Apache License, Version 2.0 (the "License");            */
/* you may not use this file except in compliance with the License.           */
/* You may obtain a copy of the License at:                                   */
/*                                                                            */
/*     http://www.apache.org/licenses/LICENSE-2.0                             */
/*                                                                            */
/* Unless required by applicable law or agreed to in writing, software        */
/* distributed under the License is distributed on an "AS IS" BASIS,          */
/* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   */
/* See the License for the specific language governing permissions and        */
/* limitations under the License.                                             */
/*                                                                            */
/******************************************************************************/

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// CBDose.hpp                                                                 //
// Corrections-Based Dose Calculation Class                                   //
// Created September 13, 2017 (Steven Dolly)                                  //
//                                                                            //
// This header file defines a class for dose calculation for a radiotherapy   //
// linac, using the corrections-based methodology. The CBDose class reads in  //
// beam data and calculates the dose to a point for a given beam.             //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

// Header guards
#ifndef CBDOSE_HPP
#define CBDOSE_HPP

// Standard C++ header files
#include <string>
#include <vector>

namespace solutio
{
  // Class to represent a linac beam
  class LinacBeam {
    public:
      // Set and get functions
      void SetFieldSize(float x, float y);
      void SetFieldSize(float x1, float x2, float y1, float y2);
      void SetSSD(float ssd){ SSD = ssd; }
      float GetX1(){ return X1; }
      float GetX2(){ return X2; }
      float GetY1(){ return Y1; }
      float GetY2(){ return Y2; }
      float GetX(){ return (X1-X2); }
      float GetY(){ return (Y1-Y2); }
      float GetSSD(){ return SSD; }
    private:
      // Beam setup data
      float X1; // X1 jaw position
      float X2; // X2 jaw position
      float Y1; // Y1 jaw position
      float Y2; // Y2 jaw position
      float SSD;
  };
  
  // Class to represent calculation point
  class CalcPoint {
    public:
      void SetPoint(float d, float doa);
      float GetDepth(){ return depth; }
      float GetOAD(){ return off_axis_distance; }
    private:
      float depth;
      float off_axis_distance;
  };
  
  // Utility calculation equations
  float SquareField(float a, float b);
  float SquareField(float r);
  float MayneordF(float f_1, float f_2, float d_0, float d);
  float AnalyticPenumbraModel(float oad, float field_size);
  
  // CBDose algorithm class:
  // 1) Loads/stores linac beam data
  // 2) Calculates dose for a given LinacBeam at a given CalcPoint, using the
  //    corrections-based formalism
  class CBDose {
    public:
      CBDose();
      // Load beam data from text file
      void LoadData(std::string file_name);
      // Get data from memory
      float Getk(){ return k; };
      float Getd_0(){ return d_0; };
      float GetSSD_0(){ return SSD_0; };
      float GetSAD(){ return SAD; };
      float GetS_c(float r);
      float GetS_p(float r);
      float GetPDD(float d, float r, float f);
      float GetTPR(float d, float r);
      float GetOAR(float d, float oad);
      // Calculation functions
      float PDDToTPR(float d, float r_d);
      float CalcDose(float mu, LinacBeam &beam, CalcPoint &point, 
          std::string type = "SAD");
      float CalcMU(float dose, LinacBeam &beam, CalcPoint &point, 
          std::string type = "SAD");
    private:
      float k; // Calibration constant in cGy/MU
      float d_0; // Depth of calibration in cm
      float SSD_0; // Calibration SSD, in cm
      float SAD; // Source-to-isocenter distance, in cm (usually 100 cm)
      float SSD_PDD; // SSD for PDD measurements, in cm
      
      std::vector<float> r_scatter;
      std::vector<float> S_c_data;
      std::vector<float> S_p_data;
      
      std::vector<float> r_pdd;
      std::vector<float> d_pdd;
      std::vector< std::vector<float> > pdd_data;
      
      std::vector<float> r_tpr;
      std::vector<float> d_tpr;
      std::vector< std::vector<float> > tpr_data;
      
      std::vector<float> oad_oar;
      std::vector<float> d_oar;
      std::vector< std::vector<float> > oar_data;
  };
}

#endif

Main File (CBDose.cpp)

CBDose.cpp - Click here to download
/******************************************************************************/
/*                                                                            */
/* Copyright 2016-2017 Steven Dolly                                           */
/*                                                                            */
/* Licensed under the Apache License, Version 2.0 (the "License");            */
/* you may not use this file except in compliance with the License.           */
/* You may obtain a copy of the License at:                                   */
/*                                                                            */
/*     http://www.apache.org/licenses/LICENSE-2.0                             */
/*                                                                            */
/* Unless required by applicable law or agreed to in writing, software        */
/* distributed under the License is distributed on an "AS IS" BASIS,          */
/* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   */
/* See the License for the specific language governing permissions and        */
/* limitations under the License.                                             */
/*                                                                            */
/******************************************************************************/

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// CBDose.cpp                                                                 //
// Corrections-Based Dose Calculation Class                                   //
// Created September 13, 2017 (Steven Dolly)                                  //
//                                                                            //
// This main file contains a class for dose calculation for a radiotherapy    //
// linac, using the corrections-based methodology. The CBDose class reads in  //
// beam data and calculates the dose to a point for a given beam.             //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

// Class header
#include "CBDose.hpp"

// Standard C++ headers
#include <iostream>
#include <fstream>
#include <sstream>

// Standard C headers
#include <cmath>

// Solutio C++ headers
#include "Utilities/DataInterpolation.hpp"

namespace solutio
{
  /////////////////////////////////////
  // Class to manage beam setup data //
  /////////////////////////////////////
  
  void LinacBeam::SetFieldSize(float x, float y)
  {
    X1 = x;
    X2 = -x;
    Y1 = y;
    Y2 = -y;
  }
  void LinacBeam::SetFieldSize(float x1, float x2, float y1, float y2)
  {
    X1 = x1;
    X2 = x2;
    Y1 = y1;
    Y2 = y2;
  }
  
  ////////////////////////////////////////////
  // Class to manage calculation point data //
  ////////////////////////////////////////////
  
  void CalcPoint::SetPoint(float d, float doa)
  {
    depth = d;
    off_axis_distance = doa;
  }
  
  ///////////////////////
  // Utility functions //
  ///////////////////////
  
  // Calculate equivalent square field sizes
  float SquareField(float a, float b){ return ( (2*a*b) / (a+b) ); }
  float SquareField(float r)
  {
    
  }
  
  // Mayneord F factor for PDD conversion to different SSD
  float MayneordF(float f_1, float f_2, float d_0, float d)
  {
    return (pow(((f_2+d_0)/(f_2+d)), 2.0 ) * pow(((f_1+d)/(f_1+d_0)), 2.0 ));
  }
  
  float AnalyticPenumbraModel(float oad, float field_size)
  {
    float A = 0.173;
    float B1 = 0.456;
    float B2 = 2.892;
    float T = 0.01;
    return (T + (1.0-T)*( A*((erf(B1*(field_size-oad))+1.0)/2.0) + 
        (1-A)*((erf(B2*(field_size-oad))+1.0)/2.0) ) );
  }
  
  ///////////////////////////////
  // Class to manage beam data //
  ///////////////////////////////
  
  CBDose::CBDose()
  {
    SAD = 100.0;
  }
  
  void CBDose::LoadData(std::string file_name)
  {
    // Initialization and open file
    std::ifstream fin;
    std::string input, str;
    size_t p1, p2;
    bool reading = true;
    float temp;
    
    fin.open(file_name.c_str());
    
    // Get first line and display
    std::getline(fin, input);
    std::cout << input << '\n';
    for(int n = 0; n < 3; n++){ std::getline(fin, input); }
    
    // Get calibration constant, SSD, and depth
    std::getline(fin, input);
    p1 = input.find(':'); p1++; p2 = input.find('c');
    str = input.substr(p1,p2-p1);
    std::stringstream(str) >> k;
    
    std::getline(fin, input);
    p1 = input.find(':'); p1++; p2 = input.find('c');
    str = input.substr(p1,p2-p1);
    std::stringstream(str) >> d_0;
    
    std::getline(fin, input);
    p1 = input.find(':'); p1++; p2 = input.find('c');
    str = input.substr(p1,p2-p1);
    std::stringstream(str) >> SSD_0;
    
    std::getline(fin, input);
    p1 = input.find(':'); p1++; p2 = input.find('c');
    str = input.substr(p1,p2-p1);
    std::stringstream(str) >> SSD_PDD;
    
    std::cout << "Calibration constant (k) = " << k << " cGy/MU @ " << d_0 << " cm\n";
    std::cout << "PDD measured using " << SSD_PDD << " cm SSD\n";
    for(int n = 0; n < 4; n++){ std::getline(fin, input); }
    
    // Get scatter factor tables for S_c and S_p
    int ind = 0;
    while(reading)
    {
      std::getline(fin, input);
      if(input == "end scatter factors")
      {
        reading = false;
      }
      else
      {
        p1 = input.find(' ');
        std::stringstream(input.substr(0,p1)) >> temp;
        r_scatter.push_back(temp);
        
        p2 = input.find(' ',p1+1);
        std::stringstream(input.substr(p1,p2-p1)) >> temp;
        S_c_data.push_back(temp);
        
        std::stringstream(input.substr(p2)) >> temp;
        S_p_data.push_back(temp);
      }
    }
    for(int n = 0; n < 3; n++){ std::getline(fin, input); }
    
    // Get PDD table
    std::getline(fin, input);
    p1 = 0; p2 = input.find(' ');
    while(p1 != std::string::npos){
      std::stringstream(input.substr(p1,p2-p1)) >> temp;
      r_pdd.push_back(temp);
      p1 = p2; p2 = input.find(' ',p1+1);
    }
    
    reading = true;
    while(reading)
    {
      std::getline(fin, input);  
      if(input == "end pdd table")
      {
        reading = false;
        continue;
      }
      
      std::vector<float> buffer;
      p1 = 0; p2 = input.find(' ');
      std::stringstream(input.substr(p1,p2-p1)) >> temp;
      d_pdd.push_back(temp);
      p1 = p2; p2 = input.find(' ',p1+1);
      while(p1 != std::string::npos)
      {
        std::stringstream(input.substr(p1,p2-p1)) >> temp;
        buffer.push_back(temp);
        p1 = p2; p2 = input.find(' ',p1+1);
      }
      
      pdd_data.push_back(buffer);
    }
    for(int n = 0; n < 2; n++){ std::getline(fin, input); }
    
    // Get TPR table
    
    // Read if TPR data is available, else calculate from PDD
    if(input == "TPR Table")
    {
      for(int n = 0; n < 2; n++){ std::getline(fin, input); }
      p1 = 0; p2 = input.find(' ');
      while(p1 != std::string::npos)
      {
        std::stringstream(input.substr(p1,p2-p1)) >> temp;
        r_tpr.push_back(temp);
        p1 = p2; p2 = input.find(' ',p1+1);
      }
      
      reading = true;
      while(reading)
      {
        std::getline(fin, input);
        if(input == "end tpr table")
        {
          reading = false;
          continue;
        }
        
        std::vector<float> buffer;
        p1 = 0; p2 = input.find(' ');
        std::stringstream(input.substr(p1,p2-p1)) >> temp;
        d_tpr.push_back(temp);
        p1 = p2; p2 = input.find(' ',p1+1);
        while(p1 != std::string::npos)
        {
          std::stringstream(input.substr(p1,p2-p1)) >> temp;
          buffer.push_back(temp);
          p1 = p2; p2 = input.find(' ',p1+1);
        }
        
        tpr_data.push_back(buffer);
      }
      for(int n = 0; n < 3; n++){ std::getline(fin, input); }
    }
    else if(input == "no tpr table")
    {
      std::cout << "No TPR data found, calculating from PDD data...\n";
      for(int n = 1; n < r_pdd.size(); n++) r_tpr.push_back(r_pdd[n]);
      for(int n = 0; n < d_pdd.size(); n++) d_tpr.push_back(d_pdd[n]);
      for(int n_d = 0; n_d < d_tpr.size(); n_d++)
      {
        std::vector<float> buffer;
        for(int n_r = 0; n_r < r_tpr.size(); n_r++)
        {
          buffer.push_back(PDDToTPR(d_tpr[n_d], r_tpr[n_r]));
        }
        tpr_data.push_back(buffer);
      }
      for(int n = 0; n < 4; n++){ std::getline(fin, input); }
    }
    else
    {
      std::cout << "Error in reading/calculating TPR data!\n";
      for(int n = 0; n < 2; n++){ std::getline(fin, input); }
    }
    
    // Get OAR table
    std::getline(fin, input);
    std::cout << input << '\n';
    p1 = 0; p2 = input.find(' ');
    while(p1 != std::string::npos)
    {
      std::stringstream(input.substr(p1,p2-p1)) >> temp;
      oad_oar.push_back(temp);
      p1 = p2; p2 = input.find(' ',p1+1);
    }
    
    reading = true;
    while(reading)
    {
      std::getline(fin, input);
      if(input == "end oar table")
      {
        reading = false;
        continue;
      }
      
      std::vector<float> buffer;
      p1 = 0; p2 = input.find(' ');
      std::stringstream(input.substr(p1,p2-p1)) >> temp;
      d_oar.push_back(temp);
      p1 = p2; p2 = input.find(' ',p1+1);
      while(p1 != std::string::npos)
      {
        std::stringstream(input.substr(p1,p2-p1)) >> temp;
        buffer.push_back(temp);
        p1 = p2; p2 = input.find(' ',p1+1);
      }
      
      oar_data.push_back(buffer);
    }
    
    // Close file
    fin.close();
  }
  
  // Get data from tables using linear interpolation
  float CBDose::GetS_c(float r)
  {
    return LinearInterpolation(r_scatter, S_c_data, r);
  }
  float CBDose::GetS_p(float r)
  {
    return LinearInterpolation(r_scatter, S_p_data, r);
  }
  float CBDose::GetPDD(float d, float r, float f)
  {
    float pdd_1 = LinearInterpolation(d_pdd, r_pdd, pdd_data, d, r);
    float pdd_2;
    if(f == SSD_PDD) pdd_2 = pdd_1;
    else
    {
      float r_1 = r*((SSD_PDD + d)/SSD_PDD);
      float r_2 = r*((f+d)/f);
      float r_10 = r*((SSD_PDD + d_0)/SSD_PDD);
      float r_20 = r*((f+d_0)/f);
      float tpr_ratio = GetTPR(d, r_2) / GetTPR(d, r_1);
      float Sp_ratio = (GetS_p(r_10)/GetS_p(r_1)) * (GetS_p(r_2)/GetS_p(r_20));
      pdd_2 = pdd_1 * MayneordF(SSD_PDD, f, d_0, d) * tpr_ratio * Sp_ratio;
    }
    return pdd_2;
  }
  float CBDose::GetTPR(float d, float r)
  {
    return LinearInterpolation(d_tpr, r_tpr, tpr_data, d, r);
  }
  
  float CBDose::GetOAR(float d, float oad)
  {
    return LinearInterpolation(d_oar, oad_oar, oar_data, d, oad);
  }
  // Convert PDD(d, r, f) to TPR(d, r_d)
  float CBDose::PDDToTPR(float d, float r_d)
  {
    float r = r_d*(SSD_PDD/(SSD_PDD+d));
    float r_d0 = r*((SSD_PDD+d_0)/SSD_PDD);
    return ( (GetPDD(d,r,SSD_PDD)/100.0) * pow(((SSD_PDD+d)/(SSD_PDD+d_0)),2.0) * (GetS_p(r_d0)/GetS_p(r_d)) );
  }
  // Calculate dose or monitor units, depending on variable "type"
  float CBDose::CalcDose(float mu, LinacBeam &beam, CalcPoint &point, 
      std::string type)
  {
    // Calculate source to point distance
    float SPD = beam.GetSSD() + point.GetDepth();
    // Calculate field sizes
    float r_c = SquareField(beam.GetX(), beam.GetY());
    float r = r_c*(beam.GetSSD() / GetSAD());
    float r_0 = r_c*((beam.GetSSD() + Getd_0()) / GetSAD());
    float r_d = r_c*(SPD / GetSAD());
    // Get scatter factors
    float S_c = GetS_c(r_c);
    float S_p;
    if(type == "SAD") S_p = GetS_p(r_d);
    else GetS_p(r_0);
    // Get PDD/TPR
    float depth_dose;
    if(type == "SAD") depth_dose = GetTPR(point.GetDepth(), r_d);
    else depth_dose = GetPDD(point.GetDepth(), r, beam.GetSSD()) / 100.0;
    // Calculate inverse square factor
    float isf;
    if(type == "SAD") isf = pow(((GetSSD_0()+Getd_0())/SPD),2.0);
    else isf = pow(((GetSSD_0()+Getd_0()) / (beam.GetSSD()+Getd_0())), 2.0);
    // Get off-axis factor
    float OAR = GetOAR(point.GetDepth(), point.GetOAD());
    // Calculate dose
    return ( mu * (Getk()*S_c*S_p*depth_dose*isf*OAR) );
  }
  
  float CBDose::CalcMU(float dose, LinacBeam &beam, CalcPoint &point, 
      std::string type)
  {
    return ( dose / CalcDose(1.0, beam, point, type) );
  }
}

Please do not hesistate to submit any questions or comments you may have about the site by using the Contact page.