LCOV - code coverage report
Current view: top level - src - src_placer.cxx (source / functions) Hit Total Coverage
Test: code_coverage_filter.info Lines: 109 193 56.5 %
Date: 2024-03-28 16:04:17 Functions: 14 22 63.6 %
Branches: 127 354 35.9 %

           Branch data     Line data    Source code
       1                 :            : /** 
       2                 :            :  * \file src_placer.cxx \brief Source Placement Utilities Implementation
       3                 :            :  *
       4                 :            :  * This source files implements the different methods for the src_placer class
       5                 :            :  */
       6                 :            : 
       7                 :            : #include "src_placer.h"
       8                 :            : 
       9                 :            : //! Default Constructor 
      10                 :       1512 : src_placer::src_placer()
      11                 :            : {
      12                 :        756 :     _power_uncertainty = 0;
      13                 :        756 : }
      14                 :            : 
      15                 :            : //! Destructor
      16                 :       1302 : src_placer::~src_placer()
      17                 :            : {
      18                 :        651 : }
      19                 :            : 
      20                 :            : //! This method reads the output vtk file of Fullmonte and stores the fluence at each element
      21                 :        105 : void src_placer::read_vtk_file(string file_path, unsigned int mesh_size)
      22                 :            : {
      23         [ +  - ]:        105 :     fprintf(stderr, "\033[1;%dm\nReading Output VTK file\033[0m\n", 33);
      24                 :            : 
      25                 :            : 
      26 [ +  - ][ +  - ]:        105 :     cout << file_path << endl;
      27                 :        105 :     _list_of_elements_fluence.clear();
      28                 :            : 
      29         [ +  - ]:        210 :     ifstream fin;
      30         [ +  - ]:        105 :     fin.open(file_path);
      31                 :            : 
      32 [ +  - ][ -  + ]:        105 :     if (!fin)
      33                 :            :     {
      34                 :          0 :         fprintf(stderr, "\033[1;%dmsrc_placer::read_vtk_file: Could not read "
      35         [ #  # ]:          0 :                         "VTK file %s\033[0m\n", 31, file_path.c_str());
      36                 :            :         
      37                 :          0 :         exit(-1);
      38                 :            :     }
      39                 :            : 
      40                 :        210 :     string line;
      41                 :        105 :     bool found_fluence = false;
      42         [ +  - ]:        210 :     vector<float> previous_fluence(mesh_size, 0); // adding this vector to account for accumulation error in Fullmonte
      43                 :        105 :     int fluence_num = 0;
      44                 :        105 :     int index = 0; 
      45 [ +  - ][ +  - ]:  411448842 :     while(getline(fin, line))
                 [ +  + ]
      46                 :            :     {
      47 [ +  - ][ +  - ]: 1234346106 :         if (line.substr(0, string("Tissue").size()) == "Tissue" || 
         [ +  - ][ +  + ]
         [ -  + ][ +  - ]
         [ +  - ][ +  - ]
           [ +  +  #  #  
             #  #  #  # ]
      48 [ +  - ][ +  - ]:  822897369 :             line.substr(0, string("Region").size()) == "Region")
         [ +  - ][ +  + ]
         [ +  + ][ +  + ]
           [ #  #  #  #  
                   #  # ]
      49                 :        105 :             found_fluence = false;
      50                 :            :         
      51 [ +  + ][ +  - ]:  411448737 :         if (found_fluence && line.substr(0, string("Fluence").size()) == "Fluence")
         [ +  - ][ +  - ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
           [ +  +  #  #  
             #  #  #  # ]
      52                 :            :         {
      53                 :            :             // Appending the fiber element number
      54 [ +  - ][ +  - ]:       6594 :             _fiber_elements.push_back(atoi(line.substr(string("Fluence").size(), line.find("_") - string("Fluence").size()).c_str()));
         [ +  - ][ +  - ]
                 [ +  - ]
      55                 :            : 
      56                 :            :             //_source_element_fluence_matrix.push_back(_list_of_elements_fluence);
      57         [ +  - ]:       6594 :             append_fluence_vector_to_matrix();
      58                 :       6594 :             _list_of_elements_fluence.clear();
      59                 :            : 
      60                 :       6594 :             index = 0; 
      61                 :            :             
      62                 :       6594 :             continue;
      63                 :            :         }
      64                 :            : 
      65         [ +  + ]:  411442143 :         if (found_fluence)
      66                 :            :         {
      67         [ +  - ]:  630268716 :             stringstream str(line);
      68                 :            :             float fluence;
      69                 :            : 
      70 [ +  - ][ +  - ]: 3151330182 :             while (str >> fluence)
                 [ +  + ]
      71                 :            :             {
      72         [ +  - ]: 2836195824 :                 _list_of_elements_fluence.push_back(fluence - previous_fluence[index]); // TODO: remove this
      73         [ -  + ]: 2836195824 :                 if (fluence-previous_fluence[index] < 0)
      74                 :            :                 {
      75 [ #  # ][ #  # ]:          0 :                     cerr << "Negative " << fluence << " " << previous_fluence[index] << " " << index << " " << fluence_num << endl;
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  # ]
                 [ #  # ]
      76                 :          0 :                     exit(-1);
      77                 :            :                 }
      78         [ -  + ]: 2836195824 :                 else if (index >= static_cast<int>(mesh_size)) 
      79                 :            :                 {
      80 [ #  # ][ #  # ]:          0 :                     cerr << "src_placer::read_vtk_file: Index " << index << " is out of bounds. Exiting! " << endl;
         [ #  # ][ #  # ]
      81                 :          0 :                     exit(-1);
      82                 :            :                 }
      83                 : 2836195824 :                 previous_fluence[index] = fluence;
      84                 : 2836195824 :                 index++;
      85                 :            :             }
      86                 :            : 
      87                 :  315134358 :             fluence_num++;
      88                 :            :         }
      89                 :            : 
      90 [ +  - ][ +  - ]:  411442143 :         if (line.substr(0, string("Fluence").size()) == "Fluence")
         [ +  - ][ +  + ]
      91                 :            :         {
      92                 :            :             // Appending the fiber element number
      93 [ +  - ][ +  - ]:        105 :             _fiber_elements.push_back(atoi(line.substr(string("Fluence").size(), line.find("_") - string("Fluence").size()).c_str()));
         [ +  - ][ +  - ]
                 [ +  - ]
      94                 :            : 
      95                 :        105 :             found_fluence = true;
      96                 :            :         }
      97                 :            :     }
      98         [ +  - ]:        105 :     append_fluence_vector_to_matrix();
      99         [ +  - ]:        105 :     fprintf(stderr, "\033[1;%dmDone Reading Output!\033[0m\n\n", 36);
     100                 :        105 : }
     101                 :            : 
     102                 :            : //! This method returns a vector of boolean indicating ignored rows corresponding to tetras not affected by the optimizer
     103                 :        284 : void src_placer::get_ignored_tetras(const unordered_set<unsigned> & tumor_region_number, 
     104                 :            :                 mesh* data, const vector<double> & d_max_tissues,
     105                 :            :                 unsigned int & actual_mesh_size, vector<bool>& tetras_to_consider)
     106                 :            : {
     107                 :            :     /** 
     108                 :            :      * @details
     109                 :            :      * The way we ignore tetra is by looking at the maximum fluence at that tetra from all 
     110                 :            :      * sources, and if it is less than D_max of that tetra then ignore (since we assume that
     111                 :            :      * the total power limit is one, so the worst case would have all power at that source
     112                 :            :      */
     113                 :            : 
     114                 :        284 :     actual_mesh_size = 0;
     115                 :            : 
     116                 :        284 :         tetras_to_consider.clear(); 
     117                 :        284 :     tetras_to_consider.resize(_source_element_fluence_matrix[0].size()); // false if ignored, true if not
     118                 :            : 
     119         [ +  + ]:  120239068 :     for (unsigned int i = 0; i < _source_element_fluence_matrix[0].size(); i++)
     120                 :            :     {
     121                 :            : 
     122                 :  120238784 :         int mat_id = data->get_list_of_tetrahedra()[i]->get_material_id();
     123   [ +  +  +  + ]:  480954568 :         if (mat_id != 0 && 
                 [ +  + ]
     124 [ +  - ][ +  + ]:  360715784 :             tumor_region_number.find(mat_id) == tumor_region_number.end())
         [ +  + ][ +  + ]
           [ #  #  #  #  
                   #  # ]
     125                 :            :         {
     126                 :  118566024 :             double max_fluence = 0;
     127         [ +  + ]: 3139494720 :             for (unsigned int j = 0; j < _source_element_fluence_matrix.size(); j++)
     128                 :            :             {
     129                 : 3020928696 :                 max_fluence = max(max_fluence, (double)_source_element_fluence_matrix[j][i]);
     130                 :            :             }
     131                 :            : 
     132         [ +  + ]:  118566024 :             if (max_fluence + DOUBLE_ACCURACY >= d_max_tissues[mat_id - 1])
     133                 :            :             {
     134                 :   38577535 :                 actual_mesh_size++;
     135                 :            : 
     136                 :   38577535 :                 tetras_to_consider[i] = true;
     137                 :            :             }
     138                 :            :             else
     139                 :  118566024 :                 tetras_to_consider[i] = false;
     140                 :            :         }
     141 [ +  - ][ +  + ]:    1672760 :         else if (tumor_region_number.find(mat_id) != tumor_region_number.end())
     142                 :            :         {
     143                 :    1672476 :             tetras_to_consider[i] = true;
     144                 :    1672476 :             actual_mesh_size++;
     145                 :            :         }
     146                 :            :         else
     147                 :        284 :             tetras_to_consider[i] = false;
     148                 :            :     }
     149                 :        284 : }
     150                 :            : 
     151                 :            : //!
     152                 :            : //! \deprecated
     153                 :            : //! This method runs the tcl script specified
     154                 :            : //!
     155                 :          0 : void src_placer::run_tcl_script(string tcl_file_path)
     156                 :            : {
     157                 :          0 :     fprintf(stderr, "\033[1;%dmRunning TCL script: %s\033[0m\n", 33, tcl_file_path.c_str());
     158                 :            : 
     159         [ #  # ]:          0 :     int command = system(("tclmonte.sh " + tcl_file_path).c_str());
     160                 :            : 
     161         [ #  # ]:          0 :     if (command == 0)
     162                 :          0 :         fprintf(stderr, "\033[1;%dmDone running script!\033[0m\n\n", 36);
     163                 :          0 : }
     164                 :            : 
     165                 :            : //!
     166                 :            : //! \deprecated
     167                 :            : //! This method reads the output of the IPOPT optimization (as a vector of source powers)
     168                 :            : //! It then generates a script of the composite sources to run FullMonte
     169                 :            : //!
     170                 :          0 : void src_placer::read_optimization_output(string filename, const vector<tetrahedron*>& sources, 
     171                 :            :                                          string mesh_path, string script_name, mesh* /*data*/,
     172                 :            :                                          const unordered_set<unsigned> & tumor_region_number, 
     173                 :            :                                          bool read_data_from_vtk)
     174                 :            : {
     175                 :            : // TODO: Deprecate this 
     176         [ #  # ]:          0 :     fprintf(stderr, "\033[1;%dmReading Optimization Output File: %s\033[0m\n", 33, filename.c_str());
     177                 :            :     
     178         [ #  # ]:          0 :     ifstream fin;
     179         [ #  # ]:          0 :     fin.open(filename);
     180                 :            : 
     181 [ #  # ][ #  # ]:          0 :     if (!fin)
     182                 :            :     {
     183                 :          0 :         fprintf(stderr, "\033[1;%dmsrc_placer::read_optimization_output: Could not read"
     184         [ #  # ]:          0 :                         "Opt. file %s\033[0m\n", 31, filename.c_str());
     185                 :            :     
     186                 :          0 :         exit(-1);
     187                 :            :     }
     188                 :            : 
     189 [ #  # ][ #  # ]:          0 :     vector<float> source_powers(sources.size());
     190                 :          0 :     unsigned int index = 0;
     191         [ #  # ]:          0 :     string line;
     192 [ #  # ][ #  # ]:          0 :     while(getline(fin, line))
                 [ #  # ]
     193                 :            :     {
     194 [ #  # ][ #  # ]:          0 :         if (line.substr(0, string("x[").size()) == "x[")
         [ #  # ][ #  # ]
     195                 :            :         {
     196         [ #  # ]:          0 :             stringstream str(line);
     197                 :          0 :             string word;
     198 [ #  # ][ #  # ]:          0 :             while(str >> word);
                 [ #  # ]
     199                 :            : 
     200         [ #  # ]:          0 :             float source_power = stof(word);
     201                 :            : 
     202                 :          0 :             source_powers[index] = source_power;
     203                 :            : 
     204                 :          0 :             index++;
     205                 :            :         }
     206                 :            :     }
     207                 :            : 
     208                 :            :     assert(sources.size() == index);
     209                 :            : 
     210         [ #  # ]:          0 :     fprintf(stderr, "\033[1;%dmDone Reading Optimization Output!\033[0m\n\n", 36);
     211                 :            : 
     212                 :            :     // Generating script
     213         [ #  # ]:          0 :     vector<tetrahedron*> sources_to_run; 
     214         [ #  # ]:          0 :     vector<float> source_powers_to_run;
     215                 :          0 :     int count_sources_to_run = 0; 
     216         [ #  # ]:          0 :     for (unsigned int i = 0; i < sources.size(); i++)
     217                 :            :     {
     218 [ #  # ][ #  # ]:          0 :         if (tumor_region_number.find((unsigned) sources[i]->get_material_id()) != tumor_region_number.end() && 
           [ #  #  #  # ]
         [ #  # ][ #  # ]
         [ #  # ][ #  #  
          #  #  #  #  #  
                      # ]
     219                 :          0 :                         source_powers[i] > SOURCE_POWER_MIN)
     220                 :            :         {
     221         [ #  # ]:          0 :             sources_to_run.push_back(sources[i]);
     222         [ #  # ]:          0 :             source_powers_to_run.push_back(source_powers[i]);
     223                 :            : 
     224                 :            : #ifdef VERBOSE_PLACER 
     225                 :            :             cout << "Source " << count_sources_to_run << " Coordinates: " << endl;
     226                 :            :             
     227                 :            :             cout << "   P1: " << data->get_list_of_points()[sources[i]->get_p1_id() - 1]->get_xcoord() << " "
     228                 :            :                         << data->get_list_of_points()[sources[i]->get_p1_id() - 1]->get_ycoord() << " "
     229                 :            :                         << data->get_list_of_points()[sources[i]->get_p1_id() - 1]->get_zcoord() << endl;
     230                 :            :             cout << "   P2: " << data->get_list_of_points()[sources[i]->get_p2_id() - 1]->get_xcoord() << " "
     231                 :            :                         << data->get_list_of_points()[sources[i]->get_p2_id() - 1]->get_ycoord() << " "
     232                 :            :                         << data->get_list_of_points()[sources[i]->get_p2_id() - 1]->get_zcoord() << endl;
     233                 :            :             cout << "   P3: " << data->get_list_of_points()[sources[i]->get_p3_id() - 1]->get_xcoord() << " "
     234                 :            :                         << data->get_list_of_points()[sources[i]->get_p3_id() - 1]->get_ycoord() << " "
     235                 :            :                         << data->get_list_of_points()[sources[i]->get_p3_id() - 1]->get_zcoord() << endl;
     236                 :            :             cout << "   P4: " << data->get_list_of_points()[sources[i]->get_p4_id() - 1]->get_xcoord() << " "
     237                 :            :                         << data->get_list_of_points()[sources[i]->get_p4_id() - 1]->get_ycoord() << " "
     238                 :            :                         << data->get_list_of_points()[sources[i]->get_p4_id() - 1]->get_zcoord() << endl;
     239                 :            : 
     240                 :            :             cout << "   Region: " << sources[i]->get_material_id() << endl;
     241                 :            : #endif
     242                 :            : 
     243                 :          0 :             count_sources_to_run++;
     244                 :            :         }
     245                 :            :     }
     246                 :            : 
     247 [ #  # ][ #  # ]:          0 :     cerr << "Number of sources to use is: " << count_sources_to_run << endl;
                 [ #  # ]
     248                 :            : 
     249                 :            :     // TODO;
     250         [ #  # ]:          0 :     return;
     251                 :            : 
     252                 :            :     if (read_data_from_vtk)
     253                 :            :         make_tcl_script_absorption_multiple_sources(script_name, 
     254                 :            :                                             mesh_path, sources_to_run, source_powers_to_run,
     255                 :            :                                             true);
     256                 :            :     else
     257                 :            :         make_tcl_script_absorption_multiple_sources(script_name, 
     258                 :            :                                             mesh_path, sources_to_run, source_powers_to_run,
     259                 :            :                                             false);
     260                 :            :     
     261                 :            :     run_tcl_script("TCL_scripts/"+script_name);
     262                 :            : }
     263                 :            : 
     264                 :            : //!
     265                 :            : //! \deprecated
     266                 :            : //! This function is modified to specify the candidate sources we want
     267                 :            : //!
     268                 :          0 : vector<tetrahedron*> src_placer::get_candidate_sources(mesh* data, const unordered_set<unsigned> & tumor_region_number,
     269                 :            :                                                 unsigned int number_of_tetra_to_skip)
     270                 :            : {
     271                 :            : 
     272                 :            :     /*
     273                 :            :      * Right now, the tumor is a cube from -1 to 1 in all directions
     274                 :            :      * The candidate sources are chosen to be every other 10th element 
     275                 :            :      * in the tumor region with an additional 0.5 units in all directions
     276                 :            :      */
     277                 :            : 
     278         [ #  # ]:          0 :     TRegion* reg = new TRegion;
     279                 :            : 
     280                 :            :     // Defining the tumor region with additional boundary
     281                 :          0 :     reg->x_min = (float) 13; //-2.7;
     282                 :          0 :     reg->y_min = (float) 64; //-2.7;
     283                 :          0 :     reg->z_min = (float) 3;  //-2.7;
     284                 :          0 :     reg->x_max = (float) 20; //-1.3;
     285                 :          0 :     reg->y_max = (float) 73; //-1.3;
     286                 :          0 :     reg->z_max = (float) 9;  //-1.3;
     287                 :            : //    13; //
     288                 :            : //    64; //
     289                 :            : //    3;  // 
     290                 :            : //    20; //
     291                 :            : //    73; //
     292                 :            : //    9;  //
     293                 :            :     tetrahedron * candidate;
     294                 :            : //     SP_Point *c_p1, *c_p2, *c_p3, *c_p4;
     295                 :          0 :     vector<tetrahedron*> sources;
     296                 :            : 
     297                 :          0 :     int current_num = 0;
     298 [ #  # ][ #  # ]:          0 :     for (unsigned int i = 0; i < data->get_mesh_size(); i++)
     299                 :            :     {
     300         [ #  # ]:          0 :         candidate = data->get_list_of_tetrahedra()[i];
     301                 :            : 
     302                 :            : //         c_p1 = data->get_list_of_points()[candidate->get_p1_id() - 1];
     303                 :            : //         c_p2 = data->get_list_of_points()[candidate->get_p2_id() - 1];
     304                 :            : //         c_p3 = data->get_list_of_points()[candidate->get_p3_id() - 1];
     305                 :            : //         c_p4 = data->get_list_of_points()[candidate->get_p4_id() - 1];
     306                 :            : 
     307                 :            : //         if (check_if_point_in_tumor_plus_region(c_p1, reg) || 
     308                 :            : //             check_if_point_in_tumor_plus_region(c_p2, reg) ||
     309                 :            : //             check_if_point_in_tumor_plus_region(c_p3, reg) || 
     310                 :            : //             check_if_point_in_tumor_plus_region(c_p4, reg))
     311 [ #  # ][ #  # ]:          0 :         if (tumor_region_number.find((unsigned)candidate->get_material_id()) != tumor_region_number.end())
                 [ #  # ]
     312                 :            :         {
     313         [ #  # ]:          0 :             if (current_num % number_of_tetra_to_skip == 0)
     314                 :            :             {
     315         [ #  # ]:          0 :                 sources.push_back(candidate);
     316                 :            :             }
     317                 :          0 :             current_num++;
     318                 :            :         }
     319                 :            :     }
     320                 :            : 
     321                 :          0 :     _num_tumor_elements = current_num;
     322                 :            : 
     323                 :          0 :     return sources;
     324                 :            : }
     325                 :            : 
     326                 :            : 
     327                 :            : // Getters
     328                 :            : 
     329                 :            : //! This method returns the fluence at each element in the mesh
     330                 :          0 : vector<float>& src_placer::get_list_of_elements_fluence()
     331                 :            : {
     332                 :          0 :     return _list_of_elements_fluence;
     333                 :            : }
     334                 :            : 
     335                 :            : //! This method returns the matrix of the fluence due to each light source at each element
     336                 :       2137 : vector<vector<float>> & src_placer::get_source_element_fluence_matrix()
     337                 :            : {
     338                 :       2137 :     return _source_element_fluence_matrix;
     339                 :            : }
     340                 :            : 
     341                 :            : //! This method returns the vector fiber elements that specifies the fiber each element belongs to
     342                 :         21 : vector<unsigned int> & src_placer::get_fiber_elements()
     343                 :            : {
     344                 :         21 :     return _fiber_elements;
     345                 :            : }
     346                 :            : 
     347                 :            : //! This method returns the number of tumor elements in the mesh
     348                 :        284 : unsigned int src_placer::get_num_tumor_elements()
     349                 :            : {
     350                 :        284 :     return _num_tumor_elements;
     351                 :            : }
     352                 :            : 
     353                 :            : //! This method returns the total energy specified by the user
     354                 :          0 : double src_placer::get_total_energy()
     355                 :            : {
     356                 :          0 :     return _total_energy;
     357                 :            : }
     358                 :            : 
     359                 :            : //! This method returns the number of photon packets used in FullMonte
     360                 :          0 : double src_placer::get_num_packets()
     361                 :            : {
     362                 :          0 :     return _num_packets;
     363                 :            : }
     364                 :            : 
     365                 :            : // Setters
     366                 :            : //! This method appends the computed vector of fluence at each element due to the last source added
     367                 :       6699 : void src_placer::append_fluence_vector_to_matrix()
     368                 :            : {
     369                 :       6699 :     long unsigned int current_num_rows = _source_element_fluence_matrix.size();
     370                 :            : 
     371                 :       6699 :     _source_element_fluence_matrix.resize(current_num_rows + 1);
     372                 :            : 
     373                 :       6699 :     _source_element_fluence_matrix[current_num_rows] = _list_of_elements_fluence;
     374                 :       6699 : }
     375                 :            : 
     376                 :            : //! This method sets the total energy
     377                 :       1512 : void src_placer::set_total_energy(double total_energy)
     378                 :            : {
     379                 :       1512 :     _total_energy = total_energy;
     380                 :       1512 : }
     381                 :            : 
     382                 :            : //! This method sets the number of photons used in FullMonte
     383                 :       1512 : void src_placer::set_num_packets(double num_packets)
     384                 :            : {
     385                 :       1512 :     _num_packets = num_packets;
     386                 :       1512 : }
     387                 :            : 
     388                 :            : 
     389                 :            : //! This method sets the power uncertainty factor (variance)
     390                 :          0 : void src_placer::set_power_uncertainty(double sigma) 
     391                 :            : {
     392                 :          0 :     _power_uncertainty = sigma;
     393                 :          0 : }
     394                 :            : 
     395                 :            : // Helper private methods
     396                 :            : //!
     397                 :            : //! This method checks if the given point is in the given tumor region in addition to some boundary
     398                 :            : //! This method is used in defining the candidate sources
     399                 :            : //!
     400                 :          0 : bool src_placer::check_if_point_in_tumor_plus_region(SP_Point* p, TRegion* reg)
     401                 :            : {
     402                 :          0 :     float upper_xlim = reg->x_max,
     403                 :          0 :           lower_xlim = reg->x_min,
     404                 :          0 :           upper_ylim = reg->y_max,
     405                 :          0 :           lower_ylim = reg->y_min,
     406                 :          0 :           upper_zlim = reg->z_max,
     407                 :          0 :           lower_zlim = reg->z_min;
     408                 :            : 
     409   [ #  #  #  # ]:          0 :     if (p->get_xcoord() < upper_xlim && p->get_xcoord() > lower_xlim &&
                 [ #  # ]
     410   [ #  #  #  # ]:          0 :         p->get_ycoord() < upper_ylim && p->get_ycoord() > lower_ylim &&
     411 [ #  # ][ #  # ]:          0 :         p->get_zcoord() < upper_zlim && p->get_zcoord() > lower_zlim
     412                 :            :        )
     413                 :          0 :         return true;
     414                 :            : 
     415                 :          0 :     return false;
     416                 :            : }
     417                 :            : 
     418                 :            : /**
     419                 :            :  * This method multiples the matrix of the fluence with the current source power  
     420                 :            :  * To find the current fluence at each element (Note that the matrix is transposed)
     421                 :            :  * @param mesh_size             [in] number of tetrahedra in the mesh
     422                 :            :  * @param current_source_powers [in] current power allocation of the sources 
     423                 :            :  * @param current_fluence_out   [out] stores the result
     424                 :            :  */
     425                 :        387 : void src_placer::compute_current_fluence(unsigned int mesh_size, 
     426                 :            :                                          const vector<double>& current_source_powers,
     427                 :            :                                                vector<double>& current_fluence_out)
     428                 :            : {
     429                 :        387 :     current_fluence_out.clear();
     430                 :        387 :     current_fluence_out.resize(mesh_size);
     431                 :            : 
     432                 :            : 
     433         [ +  + ]:  163846899 :     for (unsigned int i = 0; i < mesh_size; i++)
     434                 :            :     {
     435                 :  163846512 :         current_fluence_out[i] = 0; 
     436         [ +  + ]: 3358218432 :         for (unsigned int j = 0; j < current_source_powers.size(); j++)
     437                 :            :         {
     438         [ -  + ]: 3194371920 :             if (_source_element_fluence_matrix[0].size() == mesh_size + 1)
     439                 :          0 :                 current_fluence_out[i] += _source_element_fluence_matrix[j][i+1]*current_source_powers[j];
     440                 :            :             else
     441                 : 3194371920 :                 current_fluence_out[i] += _source_element_fluence_matrix[j][i]*current_source_powers[j];
     442                 :            :         }
     443                 :            :     }
     444                 :        387 : }
     445                 :            : 
     446                 :            : 
     447                 :            : /** 
     448                 :            :  * This method computes the v_alpha of the different tissue types
     449                 :            :  * @param current_fluences     [in] A vector of the current fluence distribution at each tetrahedron
     450                 :            :  * @param tetra_types          [in] Region ID at each tetrahedron
     451                 :            :  * @param tetra_volumes        [in] volume of each tetrahedron
     452                 :            :  * @param dose_thresholds      [in] PDT dose threshold value at each tetrahedron
     453                 :            :  * @param tumor_region_index   [in] Region ID of the tumor tissue(s)
     454                 :            :  * @param tumor_V_alpha_target [in] The percentage of the dose at which to compute the volume receiving that dose for each tissue
     455                 :            :  */
     456                 :        387 : vector<double> src_placer::compute_v_alpha_tissues(const vector<double>& current_fluences, 
     457                 :            :                                    const vector<unsigned int>& tetra_types,
     458                 :            :                                    const vector<double>& tetra_volumes,
     459                 :            :                                    const vector<double>& dose_thresholds,
     460                 :            :                                    const unordered_set<unsigned>& tumor_region_index, 
     461                 :            :                                    double tumor_V_alpha_target)
     462                 :            : {
     463                 :            : 
     464                 :            : //     double tumor_V_alpha_target = 0.9; // (0, 1) mapping to (0%, 100%)
     465                 :            : 
     466                 :        387 :     unsigned int num_regions = (unsigned int)dose_thresholds.size(); 
     467                 :            : 
     468         [ +  - ]:        387 :     vector<double> v90_tissues(num_regions);
     469                 :            : 
     470         [ +  - ]:        774 :     vector<double> region_volumes(num_regions ,0.0); 
     471         [ +  + ]:  163846899 :     for (unsigned int i = 0; i < tetra_types.size(); i++)
     472                 :            :     {
     473         [ +  + ]:  163846512 :         if (tetra_types[i] != 0)
     474                 :            :         {
     475                 :  163846125 :             region_volumes[tetra_types[i] - 1] += tetra_volumes[i];
     476                 :            :         }
     477                 :            :     }
     478                 :            :     
     479                 :            : 
     480                 :            :     // Checking which tetra exceeded its dose
     481         [ +  - ]:        774 :     vector<double> violated_volume_by_region(num_regions, 0.0);
     482         [ +  + ]:  163846899 :     for (unsigned int i = 0; i < tetra_types.size(); i++)
     483                 :            :     {
     484         [ +  + ]:  163846512 :         if (tetra_types[i] != 0)
     485                 :            :         {
     486 [ +  - ][ +  + ]:  163846125 :             if (tumor_region_index.find(tetra_types[i]) != tumor_region_index.end())
     487                 :            :             {
     488         [ +  + ]:    2279043 :                 if (current_fluences[i] < tumor_V_alpha_target*dose_thresholds[tetra_types[i] - 1])
     489                 :      47479 :                     violated_volume_by_region[tetra_types[i] - 1] += tetra_volumes[i];
     490                 :            :             }
     491                 :            :             else
     492                 :            :             {
     493         [ +  + ]:  161567082 :                 if (current_fluences[i] > tumor_V_alpha_target*dose_thresholds[tetra_types[i] - 1]) // didnt multiply by 0.9 cuz an additional guardband of 10% on the healthy
     494                 :            :                                                                                     // tissues was taken
     495                 :            :                 {
     496                 :   12302194 :                     violated_volume_by_region[tetra_types[i] - 1] += tetra_volumes[i];
     497                 :            :                 }
     498                 :            :             }
     499                 :            :         }
     500                 :            :     }
     501                 :            : 
     502         [ +  + ]:       2322 :     for (unsigned int i = 0; i < num_regions; i++)
     503                 :            :     {
     504                 :            :         // show tumor as a coverage percentage and healthy tissues as damage volume
     505 [ +  - ][ +  + ]:       1935 :         if (tumor_region_index.find(i+1) != tumor_region_index.end())
     506                 :        387 :             v90_tissues[i] = 100*(1 - (violated_volume_by_region[i]/region_volumes[i]));
     507                 :            :         else
     508                 :       1548 :             v90_tissues[i] = violated_volume_by_region[i]*1e-3; // to convert to cm^3
     509                 :            :     }
     510                 :            : 
     511                 :        774 :     return v90_tissues;
     512 [ +  - ][ +  - ]:         84 : }
     513                 :            : 
     514                 :            : 

Generated by: LCOV version 1.12