CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
castor-GATERootToCastor.cc
Go to the documentation of this file.
1 /*
2 This file is part of CASToR.
3 
4  CASToR is free software: you can redistribute it and/or modify it under the
5  terms of the GNU General Public License as published by the Free Software
6  Foundation, either version 3 of the License, or (at your option) any later
7  version.
8 
9  CASToR is distributed in the hope that it will be useful, but WITHOUT ANY
10  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  details.
13 
14  You should have received a copy of the GNU General Public License along with
15  CASToR (in file GNU_GPL.TXT). If not, see <http://www.gnu.org/licenses/>.
16 
17 Copyright 2017-2018 all CASToR contributors listed below:
18 
19  --> current contributors: Thibaut MERLIN, Simon STUTE, Didier BENOIT, Claude COMTAT, Marina FILIPOVIC, Mael MILLARDET
20  --> past contributors: Valentin VIELZEUF
21 
22 This is CASToR version 2.0.
23 */
24 
31 #include "gVariables.hh"
32 #include "gOptions.hh"
33 #include "iDataFilePET.hh"
34 #include "iDataFileSPECT.hh"
35 #include "sOutputManager.hh"
36 #include "sScannerManager.hh"
38 #include <iomanip> //std::setprecision
39 
40 #ifdef CASTOR_ROOT
41  #include "TROOT.h"
42  #include "TApplication.h"
43  #include "TGClient.h"
44  #include "TCanvas.h"
45  #include "TSystem.h"
46  #include "TTree.h"
47  #include "TBranch.h"
48  #include "TFile.h"
49 #endif
50 
51 
52 #include "oImageSpace.hh"
53 #include "oProjectorManager.hh"
54 
55 
61 void ShowHelp(int a_returnCode)
62 {
63  cout << endl;
64  cout << "Usage: castor-GATERootToCastor -i path_to_ifile.root -OR- -il path_to_ifile.txt "<< endl;
65  cout << " -s scanner_alias "<< endl;
66  cout << " -o path_to_out_file "<< endl;
67  cout << " -m path_to_macro.mac " << endl;
68  cout << " [-t only convert true prompts]" << endl;
69  cout << " [-oh histogram datafile output]" << endl;
70  cout << " [-atn path_to_atn_image(cm-1)]" << endl;
71  cout << " [-k recover coincidence kind]" << endl;
72  cout << " [-ist isotope_alias" << endl;
73  cout << " [-ot time offsets in s]" << endl;
74  cout << " [-vb verbosity]" << endl;
75  cout << endl;
76  cout << endl;
77  cout << "[Main settings]:" << endl;
78  cout << " -i path_to_file.root : give an input root datafile to convert" << endl;
79  cout << " -il path_to_file.txt : give an input text file containing a list of root files to convert" << endl;
80  cout << " : they will be concatenated into 1 CASToR datafile" << endl;
81  cout << " : the path to each datafile to convert must be entered on a newline" << endl;
82  cout << " -m path_to_macro.mac : gives the input GATE macro file used for the GATE simulation" << endl;
83  cout << " -o path_to_out_file : give the path to the output file which will be created inside this folder (no default)" << endl;
84  cout << " -s scanner_alias : provide the name of the scanner used for to acquire the original data" << endl;
85  cout << " : Must correspond to a .geom or .hscan file in the config/scanner repository." << endl;
86  cout << " : A geom file can be created from the macro files using the facultative option -geo (see below)" << endl;
87  cout << endl;
88  cout << "[Optional settings]:" << endl;
89  cout << " -t : only the true prompts will be converted" << endl;
90  cout << " -oh : Indicate if the output datafile should be written in histogram format (default : List-mode)" << endl;
91  cout << " -atn path_image:proj : For PET histogram output, provide an attenuation image (cm-1) related to the acquisition." << endl;
92  cout << " Analytic projection will be performed during the data conversion in order to estimate PET attenuation correction factors (acf) for each histogram event" << endl;
93  cout << " path_image : path to the attenuation image" << endl;
94  cout << " proj : (optional) projector to use for analytic projection. Defaut projector = Incremental siddon" << endl;
95  cout << " -k : For List-mode output, write kind of coincidence (true/scatter/rdm/...) in the output datafile (disabled by default)" << endl;
96  cout << " -ist isotope_alias : provide alias of the isotope used during the simulation"<< endl;
97  cout << " Aliases for supported PET isotopes and their parameters are listed in config/misc/isotopes_pet"<< endl;
98  cout << " Other isotopes can be added in the same file"<< endl;
99  cout << " -isrc path_to_img:dims : Provide name and dimensions (separated by a colon) of an image generated with the source (annihilation) XYZ position"<< endl;
100  cout << " The option must start with the path to the output image which will be generated." << endl;
101  cout << " Dimensions and voxel sizes of the image must be provided using commas as separators, as in the following template:" << endl;
102  cout << " path/to/image:dimX,dimY,dimZ,voxSizeX,voxSizeY,voxSizeZ"<< endl;
103  cout << " -geo : Generate a CASToR geom file from the provided GATE macro file(s)"<< endl;
104  cout << " A geom file with the 'scanner_alias' (from the -s option) basename will be generated in the scanner repository (default location : /config/scanner)" << endl;
105  cout << " -sp_bins nbinsT,nbinsA : Option specific to simulation using SPECThead systems, with root output."<< endl;
106  cout << " Give transaxial and axial number of bins for projections, separated by a comma."<< endl;
107  cout << " Pixel sizes will be computed from the crystal surface and the transaxial/axial number of bins" << endl;
108  cout << " -ot list : Provide a serie of time offset in seconds to apply before each input file"<< endl;
109  cout << " (In the case of converting several datafiles of a dynamic acquisition, timestamps of events may be resetted for each file" << endl;
110  cout << " This variable allows to manually increment the time between each datafile(s) if required" << endl;
111  cout << " The number of time offsets must be equal to the number of input files, provided by -i or -if options." << endl;
112  cout << " 'list' is a list of time offsets, separated by ','" << endl;
113  cout << endl;
114  cout << "[Miscellaneous settings]:" << endl;
115  cout << " -vb : give the verbosity level, from 0 (no verbose) to above 5 (at the event level) (default: 1)." << endl;
116  cout << endl;
117  #ifdef CASTOR_VERSION
118  cout << " This program is part of the CASToR release version " << CASTOR_VERSION << "." << endl;
119  cout << endl;
120  #endif
121  Exit(a_returnCode);
122 }
123 
124 
125 /*
126  Main program
127 */
128 
129 int main(int argc, char** argv)
130 {
131 
132  // No argument, then show help
133  if (argc==1) ShowHelp(0);
134 
135  // ============================================================================================================
136  // Parameterized variables with their default values
137  // ============================================================================================================
138 
139  // String which gathers the path to the input data filename provided by the user. no default
140  string input_file = "";
141  vector<string> path_to_input_file;
142 
143  // Path to user provided data and output files/images
144  string path_to_out_file = "";
145  string path_to_data_filename = "";
146  string path_to_header_filename = "";
147  string path_to_mac_file = "";
148  string path_to_atn_image = "";
149  string path_to_src_image = "";
150 
151  // Scanner name provided by the user
152  string scanner_name = "";
153 
154  // GATE system type
155  int GATE_system_type = GATE_SYS_UNKNOWN;
156 
157  // Only recover true GEvents
158  bool true_only_flag = false;
159 
160  // Verbosity
161  int vb = 0;
162 
163  // Histogram output
164  bool histo_out_flag = false;
165 
166  // Estimate acf (histogram output)
167  bool estimate_acf_flag = false;
168  // Projector for analytic projection
169  string options_projector = "incrementalSiddon";
170 
171  // Recover kind of coincidence (list-mode output)
172  bool kind_flag = false;
173 
174  // Isotope alias
175  string istp_alias = "";
176 
177  // Variables related to the source image
178  bool src_img_flag = false;
179  INTNB dim_src_img[3];
180  FLTNB vox_src_img[3];
181  FLTNB* p_src_img = NULL;
182 
183  // Generate geom file
184  bool geom_flag = false;
185 
186  // Input is interfile (for SPECT simulation) -> not supported (def = false)
187  bool input_is_intf = false;
188 
189  // SPECT bins
190  uint16_t spect_bin_axl = 0,
191  spect_bin_trs = 0;
192 
193  // Time offsets
194  string offset_time_str = "";
195  uint32_t* offset_time_list = NULL;
196 
197  // ============================================================================================================
198  // Read command-line parameters
199  // ============================================================================================================
200  for (int i=1; i<argc; i++)
201  {
202  string option = (string)argv[i];
203 
204  if (option=="-h" || option=="--help" || option=="-help") ShowHelp(0);
205 
206  // Just one file is provided
207  if (option=="-i")
208  {
209  if (path_to_input_file.size() > 0)
210  {
211  Cerr("***** castor-GATERootToCastor :: the following file names have already been provided (-i/-il options must be used ONCE): " << endl);
212  for (size_t i=0 ; i<path_to_input_file.size() ; i++)
213  Cout(path_to_input_file[i] << endl);
214  Exit(EXIT_FAILURE);
215  }
216  else
217  {
218  if (argv[i+1] == NULL)
219  {
220  Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl);
221  Exit(EXIT_FAILURE);
222  }
223  else
224  {
225  input_file = argv[i+1];
226  path_to_input_file.push_back(input_file);
227  }
228  i++;
229  }
230  }
231 
232  // Read a text file containing a list of root datafiles to read
233  else if (option=="-il")
234  {
235  if (path_to_input_file.size() > 0)
236  {
237  Cerr("***** castor-GATERootToCastor :: the following file names have already been provided (-i/-il options must be used ONCE) " << endl);
238  for (size_t i=0 ; i<path_to_input_file.size() ; i++)
239  Cout(path_to_input_file[i] << endl);
240  Exit(EXIT_FAILURE);
241  }
242  else
243  {
244  if (argv[i+1] == NULL)
245  {
246  Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl);
247  Exit(EXIT_FAILURE);
248  }
249  else
250  {
251  ifstream ifile(argv[i+1], ios::in);
252  string line;
253 
254  // Read list of rootfgiles
255  if(ifile)
256  {
257  // Get the position of the list_file, then append the name of the content datafiles to this position.
258  input_file = argv[i+1];
259 
260  // Recover the files
261  while (getline(ifile, line))
262  {
263  string path_to_datafile = GetPathOfFile(input_file) + OS_SEP + line;
264  path_to_input_file.push_back(path_to_datafile);
265  }
266  }
267  else
268  {
269  Cerr("***** castor-GATERootToCastor :: Error, can't read txt file: " << argv[i+1] << endl);
270  Exit(EXIT_FAILURE);
271  }
272 
273  ifile.close();
274  }
275  i++;
276  }
277  }
278  // Mac file
279  else if (option=="-m")
280  {
281  path_to_mac_file = (string)argv[i+1];
282  i++;
283  }
284  // Scanner alias
285  else if (option=="-s")
286  {
287  if (argv[i+1] == NULL)
288  {
289  Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl);
290  Exit(EXIT_FAILURE);
291  }
292  else
293  scanner_name = argv[i+1];
294  i++;
295  }
296  // Output CASToR datafile
297  else if (option=="-o")
298  {
299  if (argv[i+1] == NULL)
300  {
301  Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl);
302  Exit(EXIT_FAILURE);
303  }
304  else
305  path_to_out_file = argv[i+1];
306  i++;
307  }
308  // Recover only trues
309  else if (option=="-t")
310  {
311  #ifdef CASTOR_ROOT
312  true_only_flag = true;
313  #else
314  Cerr("***** castor-GATERootToCastor :: Option: " << option << " is only available for dataset generated with Gate in a Root format." << endl);
315  Cerr("***** castor-GATERootToCastor :: Root support is currently not enabled (CASTOR_ROOT environnement variable should be set before compilation)" << endl);
316  Exit(EXIT_FAILURE);
317  #endif
318  }
319 
320  // Time offsets
321  else if (option=="-ot")
322  {
323  if (i>=argc-1)
324  {
325  Cerr("***** castor-GATERootToCastor :: Argument missing for option: " << option << endl);
326  Exit(EXIT_FAILURE);
327  }
328  offset_time_str = (string)argv[i+1];
329  i++;
330  }
331 
332  // Output datafile in histogram mode
333  else if (option=="-oh")
334  {
335  histo_out_flag = true;
336  }
337 
338  else if (option=="-atn")
339  {
340  if (i>=argc-1)
341  {
342  Cerr("***** castor-GATERootToCastor :: Argument missing for option: " << option << endl);
343  Exit(EXIT_FAILURE);
344  }
345 
346  path_to_atn_image = argv[i+1];
347 
348  // check first if a projector has been provided (':' included)
349  size_t pos = path_to_atn_image.find_first_of(":");
350 
351  if (pos != string::npos)
352  {
353  options_projector = path_to_atn_image.substr(pos+1);
354  path_to_atn_image = path_to_atn_image.substr(0,pos);
355  }
356 
357  estimate_acf_flag = true;
358  i++;
359  }
360 
361  // Output datafile in histogram mode
362  else if (option=="-k")
363  {
364  kind_flag = true;
365  }
366 
367  // Provide an isotope alias
368  else if (option=="-ist")
369  {
370  if (argv[i+1] == NULL)
371  {
372  Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl);
373  Exit(EXIT_FAILURE);
374  }
375  else
376  istp_alias = argv[i+1];
377 
378  i++;
379  }
380 
381  // Generate image from the sourceID
382  else if (option=="-isrc")
383  {
384  if (argv[i+1] == NULL)
385  {
386  Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl);
387  Exit(EXIT_FAILURE);
388  }
389  else
390  {
391  // Recover the path to the output image
392  string input = argv[i+1];
393  int pos_colon = input.find_first_of(":");
394  path_to_src_image = input.substr(0,pos_colon);
395  input = input.substr(pos_colon + 1);
396 
397  // Get string section related to dimensions
398  string p_dims_str[6];
399  if (ReadStringOption(input.c_str(), p_dims_str, 6, ",", option))
400  {
401  Cerr("***** castor-GATERootToCastor :: Invalid argument " << argv[i+1] << " for option " << option << " !" << endl);
402  Exit(EXIT_FAILURE);
403  }
404 
405  // Recover dimensions
406  for(int i=0 ; i<3 ; i++)
407  if (ConvertFromString(p_dims_str[i], &dim_src_img[i]))
408  {
409  Cerr("***** castor-GATERootToCastor :: Conversion error for elt " << p_dims_str[i] << " for option " << option << " !" << endl);
410  Exit(EXIT_FAILURE);
411  }
412 
413  // Recover dimensions
414  for(int i=0 ; i<3 ; i++)
415  if (ConvertFromString(p_dims_str[i+3], &vox_src_img[i]))
416  {
417  Cerr("***** castor-GATERootToCastor :: Conversion error for elt " << p_dims_str[i+3] << " for option " << option << " !" << endl);
418  Exit(EXIT_FAILURE);
419  }
420 
421  // Initilize source image
422  p_src_img = new FLTNB[dim_src_img[0]*dim_src_img[1]*dim_src_img[2]];
423 
424  for(int v=0 ; v<dim_src_img[0]*dim_src_img[1]*dim_src_img[2] ; v++)
425  p_src_img[v] = 0;
426 
427  src_img_flag = true;
428  }
429  i++;
430  }
431 
432  // Generate geom file
433  else if (option=="-geo")
434  {
435  geom_flag = true;
436  }
437 
438  // SPECT bins
439  else if (option=="-sp_bins")
440  {
441  string input = argv[i+1];
442  uint16_t s_bins[2];
443  if (ReadStringOption(input.c_str(), s_bins, 2, ",", option))
444  {
445  Cerr("***** castor-GATERootToCastor :: Invalid argument " << argv[i+1] << " for option " << option << " !" << endl);
446  Exit(EXIT_FAILURE);
447  }
448 
449  spect_bin_trs = s_bins[0];
450  spect_bin_axl = s_bins[1];
451 
452  i++;
453  }
454 
455 
456  // Verbosity level
457  else if (option=="-vb")
458  {
459  if (i>=argc-1)
460  {
461  Cerr("***** castor-GATERootToCastor :: Argument missing for option: " << option << endl);
462  Exit(EXIT_FAILURE);
463  }
464  vb = atoi(argv[i+1]);
465  i++;
466  }
467 
468  else
469  {
470  Cerr("***** castor-GATERootToCastor :: Unknown option '" << option << "' !" << endl);
471  Exit(EXIT_FAILURE);
472  }
473  }
474 
475 
476  // ============================================================================================================
477  // Mandatory checks:
478  // ============================================================================================================
479 
480  // Basic initialization checks (minimal initializations mandatory for the next steps)
481 
482  // data files
483  if (path_to_input_file.empty() )
484  {
485  Cerr("***** castor-GATERootToCastor :: Please provide at least one data filename (-i or -if)" << endl);
486  ShowHelp(0);
487  Exit(EXIT_FAILURE);
488  }
489  else
490  {
491  /*
492  // Check if we have interfile
493  if(path_to_input_file.size() == 1)
494  {
495  string check;
496  int rvalue = 0;
497  rvalue = IntfKeyGetValueFromFile(path_to_input_file[0], "interfile", &check, 1, KEYWORD_OPTIONAL);
498 
499  if(rvalue == 1)
500  {
501  // error
502  Cerr("***** castor-GATERootToCastor :: Error when trying to read file: " << path_to_input_file[0] << "!" << endl);
503  Exit(EXIT_FAILURE);
504  }
505  else if(rvalue == 0)
506  // key found, we have an interfile as input
507  input_is_intf = true;
508  }
509  */
510  if(vb >= 2)
511  {
512  Cout(" Selected root data-file(s) to convert: " << endl);
513  for (size_t i=0 ; i<path_to_input_file.size() ; i++)
514  Cout(path_to_input_file[i] << endl);
515  }
516  }
517 
518 
519 
520  if(!input_is_intf)
521  {
522  // Check ROOT is enabled if the input file is not interfile (SPECT)
523  #ifndef CASTOR_ROOT
524  Cerr("***** castor-GATERootToCastor :: CASToR must be compiled with ROOT to read input root files (check installation instructions)" << endl);
525  Exit(EXIT_FAILURE);
526  #endif
527  }
528  else if(src_img_flag) // intf input + image of the source -> error
529  {
530  Cerr("***** castor-GATERootToCastor :: Can't use -isrc with interfile dataset ");
531  Cerr(" (image of the source can only be generated from root file) !" << endl);
532  Exit(EXIT_FAILURE);
533  }
534 
535  // Check file(s) existence
536  for (size_t i=0 ; i<path_to_input_file.size() ; i++)
537  {
538  ifstream input_file(path_to_input_file[i].c_str(), ios::in);
539  if (!input_file)
540  {
541  Cerr("***** gOptions::ReadDataASCIIFile() -> Couldn't find or read data-file '"<< path_to_input_file[i] << "' !" << endl);
542  Exit(EXIT_FAILURE);
543  }
544  }
545 
546  // output directory
547  if (path_to_out_file.empty() )
548  {
549  Cerr("***** castor-GATERootToCastor :: Please provide the output file name (-o)" << endl);
550  ShowHelp(0);
551  Exit(EXIT_FAILURE);
552  }
553  else
554  if(vb >= 2) Cout(" selected output file:" << path_to_out_file << endl);
555 
556  // macro
557  if (path_to_mac_file.empty())
558  {
559  Cerr("***** castor-GATERootToCastor :: Please provide the macro file associated to the GATE root datafile (-m) :" << endl);
560  ShowHelp(0);
561  Exit(EXIT_FAILURE);
562  }
563  else
564  if(vb >= 2) Cout(" selected macro file: " << path_to_mac_file << endl);
565 
566  // scanner
567  if (scanner_name.empty())
568  {
569  Cerr("***** castor-GATERootToCastor :: Please provide a scanner alias (-s) :" << endl);
570  ShowHelp(0);
571  Exit(EXIT_FAILURE);
572  }
573  else
574  if(vb >= 2) Cout(" selected scanner alias: " << scanner_name << endl);
575 
576 
577  // (Required for options using interfile I/O)
579 
580 
581  // ============================================================================================================
582  // SOutputManager object initialisation:
583  // ============================================================================================================
584 
585  sOutputManager* p_outputManager = sOutputManager::GetInstance();
586  // Set verbose level
587  p_outputManager->SetVerbose(vb);
588  // Set MPI rank
589  p_outputManager->SetMPIRank(0);
590 
591  // Set path to the config directory
592  if (p_outputManager->CheckConfigDir(""))
593  {
594  Cerr("***** castor-GATERootToCastor :: A problem occured while checking for the config directory path !" << endl);
595  Exit(EXIT_FAILURE);
596  }
597  // Initialize output directory and base name
598  if (p_outputManager->InitOutputDirectory(path_to_out_file, ""))
599  {
600  Cerr("*****castor-GATERootToCastor :: A problem occured while initializing output directory !" << endl);
601  Exit(EXIT_FAILURE);
602  }
603  // Log command line
604  if (p_outputManager->LogCommandLine(argc,argv))
605  {
606  Cerr("***** castor-GATERootToCastor :: A problem occured while logging command line arguments !" << endl);
607  Exit(EXIT_FAILURE);
608  }
609 
610  // Output progression options
611  cout << std::fixed;
612  cout << std::setprecision(2);
613 
614 
615  // ============================================================================================================
616  // Check system type from the macro file
617  // ============================================================================================================
618  GATE_system_type = GetGATESystemType(path_to_mac_file);
619 
620  if(GATE_system_type<0)
621  {
622  // Unknown system
623  Cerr("***** castor-GATERootToCastor :: GATE system type not found : " << endl);
624  Cerr(" This script only supports conversion for cylindricalPET, ecat and SPECThead systems" << endl);
625  Cerr(" The system type is recovered from the lines '/gate/systems/...'" << endl);
626  Exit(EXIT_FAILURE);
627  }
628 
629 
630  // ============================================================================================================
631  // Generate the geom file from the mac file(s) is required
632  // ============================================================================================================
633  if(geom_flag)
634  {
635  string scanner_repository = sOutputManager::GetInstance()->GetPathToConfigDir() + "scanner" + OS_SEP;
636  string path_to_geom = scanner_repository + scanner_name + ".geom";
637 
638  // Check system type
639  switch ( GATE_system_type )
640  {
642  if( vb>=2 )Cout(endl << " --- CylindricalPET system detected. Proceeding to conversion... --- " << endl << endl);
643 
644  if(CreateGeomWithCylindrical(path_to_mac_file , path_to_geom) )
645  {
646  Cerr("***** castor-GATERootToCastor :: An error occured while trying to process mac file for cylindrical system: " << path_to_mac_file << endl);
647  Exit(EXIT_FAILURE);
648  }
649  break;
650 
651  case GATE_SYS_ECAT:
652  if( vb>=2 )Cout(endl << " --- ECAT system detected. Proceeding to conversion... --- " << endl << endl);
653 
654  if(CreateGeomWithECAT(path_to_mac_file , path_to_geom) )
655  {
656  Cerr("***** castor-GATEMacToGeom :: An error occured while trying to process mac file for ecat system: " << path_to_mac_file << endl);
657  Exit(EXIT_FAILURE);
658  }
659  break;
660  // TODO
661  case GATE_SYS_SPECT:
662  if( vb>=2 )Cout(endl << " --- SPECThead system detected. Proceeding to conversion... --- " << endl << endl);
663 
664  if(CreateGeomWithSPECT(path_to_mac_file , path_to_geom) )
665  {
666  Cerr("***** castor-GATEMacToGeom :: An error occured while trying to process mac file for SPECT system: " << path_to_mac_file << endl);
667  Exit(EXIT_FAILURE);
668  }
669  break;
670 
671  default: // Unknown system
672  Cerr("***** castor-GATERootToCastor :: System type not found : " << endl);
673  Cerr(" This script only supports conversion for cylindricalPET ecat and SPECThead systems" << endl);
674  Cerr(" The system type is recovered from the lines '/gate/systems/...'" << endl);
675  Exit(EXIT_FAILURE);
676  break;
677  }
678 
679  if( vb>=2 )Cout(endl << " --- Conversion completed --- " << endl << endl);
680  }
681 
682 
683  // ============================================================================================================
684  // ScannerManager object initialisation:
685  // ============================================================================================================
686 
687  sScannerManager* p_scannermanager = sScannerManager::GetInstance();
688  p_scannermanager->SetVerbose(vb);
689 
690  // Check if the scanner exists and get the name from ScannerManager
691  scanner_name = (scanner_name.find(OS_SEP)) ?
692  scanner_name.substr(scanner_name.find_last_of(OS_SEP)+1) :
693  scanner_name;
694 
695  if(p_scannermanager->FindScannerSystem(scanner_name) )
696  {
697  Cerr("**** castor-GATERootToCastor :: A problem occurred while searching for scanner system !" << endl);
698  Exit(EXIT_FAILURE);
699  }
700 
701  // Build output file names
702  path_to_data_filename = path_to_out_file + ".Cdf";
703  path_to_header_filename = path_to_out_file + ".Cdh";
704 
705 
706  // ============================================================================================================
707  // Instanciate/Initialize CASToR DataFile
708  // ============================================================================================================
709 
710  // Instantiate & Initialize oImageDimensionsAndQuantification object, required for datafile generation (number of threads)
712 
713  // Set isotope if provided
714  if(!istp_alias.empty())
715  {
716  if(p_ID->SetPETIsotope(0, istp_alias) )
717  {
718  Cerr("**** castor-GATERootToCastor :: A problem occurred while checking isotope name !" << endl);
719  Exit(EXIT_FAILURE);
720  }
721  }
722 
723  // Instanciate & Initialize iDataFilePET and Event objects
724  vDataFile* Out_data_file = NULL;
725  vEvent* Event = NULL;
726 
727  uint16_t max_nb_lines_per_event = 1; // No compression for root files
728 
729  if(GATE_system_type == GATE_SYS_SPECT)
730  {
731  Out_data_file = new iDataFileSPECT();
732  iDataFileSPECT* p_datafile = (dynamic_cast<iDataFileSPECT*>(Out_data_file));
733  p_datafile->SetDataType(TYPE_SPECT);
734  p_datafile->SetIsotope(istp_alias);
735  histo_out_flag = true; // force histogram output for SPECT
736  }
737  else
738  {
739  Out_data_file = new iDataFilePET();
740  iDataFilePET* p_datafile = (dynamic_cast<iDataFilePET*>(Out_data_file));
741  p_datafile->SetDataType(TYPE_PET);
742  p_datafile->SetIsotope(istp_alias);
743 
744  // ACF computed
745  if(estimate_acf_flag)
746  p_datafile->SetAtnCorrectionFlagOn();
747 
748  p_datafile->SetMaxNumberOfLinesPerEvent(max_nb_lines_per_event);
749  }
750 
751  Out_data_file->SetImageDimensionsAndQuantification(p_ID);
752  Out_data_file->SetHeaderDataFileName(path_to_header_filename);
753  Out_data_file->SetVerbose(0);
754 
755 
756  // Init Histogram-mode Event
757  if(histo_out_flag)
758  {
759  Out_data_file->SetDataMode(MODE_HISTOGRAM);
760 
761  // Instanciate histogram Event depending on modality
762  if(GATE_system_type == GATE_SYS_SPECT)
763  Event = new iEventHistoSPECT();
764  else
765  Event = new iEventHistoPET();
766 
767  Event->SetNbLines(max_nb_lines_per_event);
768  Event->AllocateID();
769  }
770 
771  // Init List-mode Event
772  else
773  {
774  Out_data_file->SetDataMode(MODE_LIST);
775 
776  // Instanciate histogram Event depending on modality
777  if(GATE_system_type == GATE_SYS_SPECT)
778  {
779  Event = new iEventListSPECT();
780  // record coincidence kind or not
781  if(kind_flag)
782  ((iDataFileSPECT*)Out_data_file)->SetEventKindFlagOn();
783  }
784  else
785  {
786  Event = new iEventListPET();
787  // record coincidence kind or not
788  if(kind_flag)
789  ((iDataFilePET*)Out_data_file)->SetEventKindFlagOn();
790  }
791 
792  Event->SetNbLines(max_nb_lines_per_event);
793  Event->AllocateID();
794  }
795 
796  Out_data_file->PROJ_InitFile();
797  Out_data_file->ComputeSizeEvent();
798  Out_data_file->PrepareDataFile();
799 
800 
801  // ============================================================================================================
802  // If acf estimation is enabled for histogram output, initialize all required objects for analytic projection
803  // (Image space, projector and scanner geometry)
804  // ============================================================================================================
805  oProjectorManager* p_ProjectorManager = new oProjectorManager();
806  oImageSpace* p_ImageSpace = new oImageSpace();
807 
808 
809  if(estimate_acf_flag)
810  {
811  // Check if histogram output has been enabled, throw error otherwise
812  if(!histo_out_flag)
813  {
814  Cerr("***** castor-GATERootToCastor :: Estimation of acf from an attenuation image (-atn option) is only available for histogram datafile output !" << endl
815  << " Add the (-oh) option to the command line to enable this option." << endl);
816  Exit(1);
817  }
818 
819  // Check if system is SPECT, throw error in this case
820  if(GATE_system_type == GATE_SYS_SPECT)
821  {
822  Cerr("***** castor-GATERootToCastor :: Estimation of acf from an attenuation image (-atn option) only available for PET systems ! (detected system is SPECT)" << endl);
823  Exit(1);
824  }
825 
826  Intf_fields IF;
827  IntfKeyInitFields(&IF);
828  if(IntfReadHeader(path_to_atn_image, &IF, vb) )
829  {
830  Cerr("***** castor-GATERootToCastor :: An error occurred while trying to read the interfile header of attenuation file " << path_to_atn_image << " !" << endl);
831  Exit(1);
832  }
833 
834  // --- oImageDimensionsAndQuantification initialization ---
835  p_ID->SetNbVoxX(IF.mtx_size[0]);
836  p_ID->SetNbVoxY(IF.mtx_size[1]);
837  p_ID->SetNbVoxZ(IF.mtx_size[2]);
838  p_ID->SetNbThreads("1");
839  p_ID->SetNbBeds(1);
840  p_ID->SetVoxSizeX(IF.vox_size[0]);
841  p_ID->SetVoxSizeY(IF.vox_size[1]);
842  p_ID->SetVoxSizeZ(IF.vox_size[2]);
843  p_ID->SetFOVOutMasking(0., 0);
844  p_ID->SetFOVSizeX(-1.);
845  p_ID->SetFOVSizeY(-1.);
846  p_ID->SetFOVSizeZ(-1.);
847  p_ID->SetOffsetX(0);
848  p_ID->SetOffsetY(0);
849  p_ID->SetOffsetZ(0);
850  p_ID->SetVerbose(vb);
851  p_ID->SetNbRespGates(1);
852  p_ID->SetNbCardGates(1);
853  p_ID->SetFrames("");
854 
855  if (p_ID->CheckParameters())
856  {
857  Cerr("***** castor-GATERootToCastor :: A problem occured while checking image dimensions parameters !" << endl);
858  Exit(1);
859  }
860  if (p_ID->Initialize())
861  {
862  Cerr("***** castor-GATERootToCastor :: A problem occured while initializing image dimensions !" << endl);
863  Exit(1);
864  }
865 
866  // Initialization of DynamicDataManager class, related 4D data splitting management
867  if (p_ID->InitDynamicData( "", 0, 0, 0, 0, 1, 1 ) )
868  {
869  Cerr("***** castor-GATERootToCastor :: A problem occured while initializing Dynamic data manager's class !" << endl);
870  Exit(EXIT_FAILURE);
871  }
872 
873 
874  // --- Image space initialization ---
875  p_ImageSpace->SetImageDimensionsAndQuantification(p_ID);
876  p_ImageSpace->SetVerbose(vb);
877 
878  // Allocate memory for main image
879  p_ImageSpace->LMS_InstantiateImage();
880 
881  // Read attenuation image
882  if(p_ImageSpace->PROJ_LoadInitialImage(path_to_atn_image) )
883  {
884  Cerr("***** castor-GATERootToCastor :: Error during image initialization !" << endl);
885  Exit(EXIT_FAILURE);
886  }
887 
888 
889  // --- Build Scanner geometry ---
890  if(p_scannermanager->BuildScannerObject() )
891  {
892  Cerr("**** castor-GATERootToCastor :: A problem occurred during scanner object construction !" << endl);
893  Exit(EXIT_FAILURE);
894  }
895 
896  if(p_scannermanager->InstantiateScanner() )
897  {
898  Cerr("**** castor-GATERootToCastor :: A problem occurred while creating Scanner object !" << endl);
899  Exit(EXIT_FAILURE);
900  }
901 
902  if(p_scannermanager->BuildLUT() )
903  {
904  Cerr("***** castor-GATERootToCastor :: A problem occurred while generating/reading the LUT !" << endl);
905  Exit(EXIT_FAILURE);
906  }
907 
908  // Check the scanner manager parameters and initialize the scanner
909  if (p_scannermanager->CheckParameters())
910  {
911  Cerr("***** castor-GATERootToCastor :: A problem occured while checking scanner manager parameters !" << endl);
912  Exit(1);
913  }
914 
915  if (p_scannermanager->Initialize())
916  {
917  Cerr("***** castor-GATERootToCastor :: A problem occured while initializing scanner !" << endl);
918  Exit(1);
919  }
920 
921 
922  // --- Projector Manager initialization---
923  p_ProjectorManager->SetScanner(p_scannermanager->GetScannerObject());
924  p_ProjectorManager->SetImageDimensionsAndQuantification(p_ID);
925  p_ProjectorManager->SetDataFile(Out_data_file);
927  p_ProjectorManager->SetOptionsForward(options_projector);
928  p_ProjectorManager->SetOptionsBackward(options_projector);
929  p_ProjectorManager->SetOptionsCommon("");
930  p_ProjectorManager->SetVerbose(vb);
931 
932  // Check parameters
933  if (p_ProjectorManager->CheckParameters())
934  {
935  Cerr("***** castor-GATERootToCastor :: A problem occured while checking projector manager's parameters !" << endl);
936  Exit(EXIT_FAILURE);
937  }
938 
939  // Initialize projector manager
940  if (p_ProjectorManager->Initialize())
941  {
942  Cerr("***** castor-GATERootToCastor :: A problem occured while initializing projector manager !" << endl);
943  Exit(EXIT_FAILURE);
944  }
945  }
946 
947 
948  // ============================================================================================================
949  // Parse Time offsets
950  // ============================================================================================================
951  offset_time_list = new uint32_t[path_to_input_file.size()];
952  for(size_t f=0 ; f<path_to_input_file.size() ; f++)
953  offset_time_list[f] = 0;
954 
955  // Parse offset_time_str, if it contains any data
956  if(offset_time_str != "")
957  {
958  vector<string> offsets;
959  size_t comma_pos = 0;
960  while ((comma_pos=offset_time_str.find_first_of(",")) != string::npos)
961  {
962  string offset = offset_time_str.substr(0,comma_pos);
963  offsets.push_back(offset);
964  offset_time_str = offset_time_str.substr(comma_pos+1);
965  }
966 
967  // Check we have a correct number of offsets
968  if(offsets.size() != path_to_input_file.size())
969  {
970  Cerr("**** castor-GATERootToCastor :: Unmatching number of offsets with -ot option ! "
971  << offsets.size() << " have been found, while "<< path_to_input_file.size() <<"input file have been provided !" << endl);
972  Exit(EXIT_FAILURE);
973  }
974 
975  for(size_t o=0 ; o<offsets.size() ; o++)
976  {
977  double offset;
978  if(ConvertFromString(offsets[o] , &offset) )
979  {
980  Cerr("**** castor-GATERootToCastor :: Error while trying to convert offset : "<< offsets[o] <<" in ms ! " << endl);
981  Exit(EXIT_FAILURE);
982  }
983 
984  // convert in ms
985  offset_time_list[o] = (uint32_t)offset*1000;
986  }
987  }
988 
989 
990  // ============================================================================================================
991  // *********************************************** CONVERSION *************************************************
992  // ============================================================================================================
993  Cout(" --- Start conversion of datafile(s) : " << input_file << endl
994  << " using mac file: " << path_to_mac_file << endl
995  << " CASToR output header datafile: " << path_to_header_filename << endl
996  << " CASToR output binary datafile: " << path_to_data_filename << endl << endl);
997 
998  // ============================================================================================================
999  // Variables initialization
1000  // ============================================================================================================
1001 
1002  // Counter for the number of events (and the number of trues, scatters and randoms for simulated data)
1003  uint64_t nLORs_tot = 0,
1004  nLORs_trues = 0,
1005  nLORs_rdms = 0,
1006  nLORs_scatters =0,
1007  nLORs_unknown =0,
1008  nBINs = 0;
1009 
1010  // Scanner variables (PET)
1011  uint32_t nCrystalsTot = 0;
1012  uint32_t nRsectorsAngPos = 1, nRsectorsAxial = 1;
1013  // index order of rsectors (transaxial/axial) for cylindricalPET depending on the GATE macro
1014  int rsector_id_order = 0;
1015  uint32_t nModulesTransaxial = 1, nModulesAxial = 1;
1016  uint32_t nSubmodulesTransaxial = 1, nSubmodulesAxial = 1;
1017  uint32_t nCrystalsTransaxial = 1, nCrystalsAxial = 1;
1018  uint32_t nBlocksLine = 1, nBlocksPerRing = 1;
1019  uint8_t nLayers = 1;
1020  uint32_t nb_crystal_per_layer[4] = {0,0,0,0};
1021 
1022  // Scanner variables (SPECT)
1023  uint32_t nProjectionsByHead =1;
1024  uint32_t nProjectionsTot = 1;
1025  uint32_t nHeads =1;
1026  uint32_t nPixelsAxial = 1;
1027  uint32_t nPixelsTransaxial = 1;
1028  uint32_t nbSimulatedPixels = 1;
1029  float_t distToDetector = 0.;
1030  int headRotDirection = GEO_ROT_CW;
1031  float_t head1stAngleDeg = 0;
1032  float_t headAngPitchDeg = -1;
1033  float_t headAngStepDeg = -1;
1034  float_t crystalSizeAxl=-1.,
1035  crystalSizeTrs=-1.;
1036  FLTNB* p_proj_spect_image=NULL;
1037 
1038  // layer repeaters
1039  vector<uint32_t> nLayersRptTransaxial, nLayersRptAxial;
1040 
1041  // Castor variables
1042  uint8_t** p_bins = NULL;
1043  uint32_t bins_elts = 0;
1044  uint32_t start_time_ms = 0; // 0 as default initialization
1045  uint32_t duration_ms = 1000; // 1 second as default initialization
1046 
1047 
1048  #ifdef CASTOR_ROOT
1049  uint32_t castorID1=0;
1050  uint32_t castorID2=0;
1051  uint8_t kind;
1052 
1053  // ROOT data variables
1054  int32_t crystalID1=0, crystalID2=0;
1055 
1056  // cylindricalPET specific
1057  int32_t rsectorID1=0 , rsectorID2=0;
1058  int32_t moduleID1=0 , moduleID2=0;
1059  int32_t submoduleID1=0, submoduleID2=0;
1060  int32_t layerID1=0 , layerID2=0;
1061 
1062  // ecat specific
1063  int32_t blockID1=0, blockID2=0;
1064 
1065  // SPECThead specific
1066  float_t rotAngle;
1067  float_t gPosX, gPosY, gPosZ;
1068  int32_t headID=0;
1069  int32_t pixelID=0;
1070 
1071 
1072  // others
1073  int32_t eventID1= 0, eventID2= 0;
1074  int32_t comptonPhantomID1 = 0, comptonPhantomID2= 0;
1075  int32_t rayleighPhantomID1 = 0,rayleighPhantomID2 = 0;
1076  double_t time1= 0, time2= 0;
1077  int32_t sourceID1=0, sourceID2=0;
1078  float_t sourcePosX1=0., sourcePosY1=0., sourcePosZ1=0.;
1079 
1080 
1081  // ============================================================================================================
1082  // ROOT objects declaration
1083  // ============================================================================================================
1084 
1085  TApplication *Tapp = new TApplication("tapp",0,0);
1086  TTree** GEvents = new TTree *[path_to_input_file.size()];
1087  TFile** Tfile_root = new TFile*[path_to_input_file.size()];
1088 
1089  if(!input_is_intf)
1090  {
1091  // Compute the total number of LORs in the dataset
1092  for (size_t iFic=0 ; iFic<path_to_input_file.size() ; iFic++)
1093  {
1094  Tfile_root[iFic] = new TFile(path_to_input_file[iFic].c_str(),"READ","ROOT file with histograms");
1095  if(GATE_system_type == GATE_SYS_SPECT)
1096  GEvents[iFic] = (TTree*)Tfile_root[iFic]->Get("Singles");
1097  else
1098  GEvents[iFic] = (TTree*)Tfile_root[iFic]->Get("Coincidences");
1099 
1100  nLORs_tot += GEvents[iFic]->GetEntries();
1101  }
1102  }
1103  #endif
1104 
1105 
1106  // Indexes for progression output
1107  uint64_t printing_index = 0;
1108  uint64_t printing_ratio = (nLORs_tot>10000) ? 10000 : nLORs_tot/10;
1109 
1110  // ============================================================================================================
1111  // Recover system geometric elements values
1112  // ============================================================================================================
1113 
1114  // SPECThead system
1115  if(GATE_system_type == GATE_SYS_SPECT)
1116  {
1117  duration_ms =0.;
1118 
1119  // Data provided as an interfile projection image
1120  if(input_is_intf)
1121  {
1122  if(ReadIntfSPECT( path_to_input_file[0],
1123  distToDetector,
1124  nHeads,
1125  nPixelsAxial,
1126  nPixelsTransaxial,
1127  crystalSizeAxl,
1128  crystalSizeTrs,
1129  nProjectionsTot,
1130  nProjectionsByHead,
1131  head1stAngleDeg,
1132  headAngPitchDeg,
1133  headAngStepDeg,
1134  headRotDirection,
1135  start_time_ms,
1136  duration_ms,
1137  vb) )
1138  {
1139  Cerr("**** castor-GATERootToCastor :: Error when reading mac file : "<< path_to_mac_file <<" !" << endl);
1140  Exit(EXIT_FAILURE);
1141  }
1142 
1143  nCrystalsTot = nPixelsAxial * nPixelsTransaxial;
1144 
1145  // Read Interfile
1146  Intf_fields IF_proj_spect_data;
1147  p_proj_spect_image = new FLTNB[nHeads*nProjectionsByHead*nCrystalsTot];
1148 
1149  if( IntfReadProjectionImage(path_to_input_file[0], p_proj_spect_image, &IF_proj_spect_data, vb, false) )
1150  {
1151  Cerr("**** castor-GATERootToCastor :: Error when trying to read image : "<< path_to_input_file[0] <<" !" << endl);
1152  Exit(EXIT_FAILURE);
1153  }
1154  }
1155 
1156  // Data provided as a root file
1157  else
1158  {
1159  if(ReadMacSPECT(path_to_mac_file,
1160  distToDetector,
1161  nHeads,
1162  nPixelsAxial,
1163  nPixelsTransaxial,
1164  crystalSizeAxl,
1165  crystalSizeTrs,
1166  nProjectionsTot,
1167  nProjectionsByHead,
1168  head1stAngleDeg,
1169  headAngPitchDeg,
1170  headAngStepDeg,
1171  headRotDirection,
1172  start_time_ms,
1173  duration_ms,
1174  vb) )
1175  {
1176  Cerr("**** castor-GATERootToCastor :: Error when reading mac file : "<< path_to_mac_file <<" !" << endl);
1177  Exit(EXIT_FAILURE);
1178  }
1179 
1180  nbSimulatedPixels = nPixelsAxial*nPixelsTransaxial;
1181 
1182  if((spect_bin_trs>0 || spect_bin_axl>0) && nbSimulatedPixels > 1 )
1183  {
1184  Cerr("**** castor-GATERootToCastor :: WARNING : Spect bins have been initialized, but the simulation already provide a specific number of pixels (="<< nPixelsAxial*nPixelsTransaxial <<") !"<< endl <<
1185  " Pixel matrix used by default !" << endl);
1186  Exit(EXIT_FAILURE);
1187  }
1188  else // Check bins have been provided, and initialize pixel sizes with provided values
1189  {
1190  if(spect_bin_trs == 0 && spect_bin_axl == 0 && nbSimulatedPixels==1)
1191  {
1192  Cerr("**** castor-GATERootToCastor :: Error : Axial and transaxial bins values expected (use option -spect_b) !"<< endl);
1193  Exit(EXIT_FAILURE);
1194  }
1195 
1196  // check crystal sizes have been correctly read
1197  if(crystalSizeAxl<0 || crystalSizeTrs<0)
1198  {
1199  Cerr("**** castor-GATERootToCastor :: Crystal dimensions not correctly read in the mac files !" << endl);
1200  Exit(EXIT_FAILURE);
1201  }
1202 
1203  nPixelsTransaxial = spect_bin_trs;
1204  nPixelsAxial = spect_bin_axl;
1205 
1206  if(vb>=2) Cout("Transaxial/Axial nb pixels : " << nPixelsTransaxial << " , " << nPixelsAxial << endl <<
1207  "Transaxial/Axial pixel sizes : " << crystalSizeTrs/nPixelsTransaxial << " , " << crystalSizeAxl/nPixelsAxial << endl);
1208  }
1209  }
1210 
1211  nCrystalsTot = nPixelsAxial * nPixelsTransaxial;
1212 
1213  if(vb>=2) Cout("Number of Projections: " << nProjectionsByHead << endl <<
1214  "Detected number of crystals in the system : " << nCrystalsTot << endl);
1215 
1216 
1217  // Histogram bin vector initialization using the total number of projections & pixels
1218  if(histo_out_flag)
1219  {
1220  bins_elts = nHeads*nProjectionsByHead;
1221 
1222  p_bins = new uint8_t*[bins_elts];
1223  for (size_t p=0; p<bins_elts ; p++)
1224  {
1225  p_bins[p] = new uint8_t[nCrystalsTot];
1226 
1227  for (size_t c=0; c<nCrystalsTot ; c++)
1228  {
1229  p_bins[p][c] = input_is_intf ? (uint8_t)p_proj_spect_image[p*nCrystalsTot+c] : 0 ;
1230 
1231  // Get number of singles if input is interfile proj image
1232  if(input_is_intf)
1233  {
1234  nLORs_tot += p_proj_spect_image[p*nCrystalsTot+c];
1235  nLORs_unknown += p_proj_spect_image[p*nCrystalsTot+c];
1236  }
1237  }
1238  }
1239  }
1240 
1241  delete[] p_proj_spect_image;
1242  }
1243 
1244  else // PET systems
1245  {
1246  // cylindricalPET system
1247  if(GATE_system_type == GATE_SYS_CYLINDRICAL)
1248  {
1249  if(ReadMacCylindrical(path_to_mac_file,
1250  nLayers,
1251  nb_crystal_per_layer,
1252  nCrystalsTot,
1253  nCrystalsAxial,
1254  nCrystalsTransaxial,
1255  nLayersRptAxial,
1256  nLayersRptTransaxial,
1257  nSubmodulesAxial,
1258  nSubmodulesTransaxial,
1259  nModulesAxial,
1260  nModulesTransaxial,
1261  nRsectorsAxial,
1262  nRsectorsAngPos,
1263  rsector_id_order,
1264  start_time_ms,
1265  duration_ms,
1266  vb) )
1267  {
1268  Cerr("**** castor-GATERootToCastor :: Error when reading mac file : "<< path_to_mac_file <<" !" << endl);
1269  Exit(EXIT_FAILURE);
1270  }
1271  }
1272 
1273  else // ECAT
1274  {
1275  // Reading the macro file
1276  if(ReadMacECAT(path_to_mac_file,
1277  nCrystalsTot,
1278  nCrystalsAxial,
1279  nCrystalsTransaxial,
1280  nBlocksLine,
1281  nBlocksPerRing,
1282  start_time_ms,
1283  duration_ms,
1284  vb) )
1285  {
1286  Cerr("**** castor-GATERootToCastor :: Error when reading mac file : "<< path_to_mac_file <<" !" << endl);
1287  Exit(EXIT_FAILURE);
1288  }
1289 
1290  }
1291 
1292  if(vb>=2) Cout("Detected number of crystals in the system : " << nCrystalsTot << endl);
1293 
1294  // Histogram bin vector initialization using the total number of crystals
1295  if(histo_out_flag)
1296  {
1297  Cout(" Allocating memory for bins... " << endl <<
1298  " Warning : this step can require huge amount of memory if the system contains a high number of crystals !" << endl);
1299 
1300  bins_elts = nCrystalsTot;
1301 
1302  p_bins = new uint8_t*[bins_elts];
1303  for (size_t c=0; c<bins_elts ; c++)
1304  {
1305  p_bins[c] = new uint8_t[nCrystalsTot-c];
1306 
1307  for (size_t c2=0; c2<nCrystalsTot-c ; c2++)
1308  p_bins[c][c2] = 0;
1309  }
1310 
1311  Cout(" Memory allocation for bins completed " << endl << endl);
1312  }
1313 
1314  } // end of PET section
1315 
1316 
1317 
1318 
1319  // ============================================================================================================
1320  // Loop on the input datafile(s) and process event by event
1321  // ============================================================================================================
1322  if(!input_is_intf) // SPECT projection interfile : data already processed
1323  for (size_t iFic=0 ; iFic<path_to_input_file.size() ; iFic++)
1324  {
1325  if(vb>=2) Cout(endl << "Processing datafile " << iFic << " : " << path_to_input_file[iFic] << "..." << endl);
1326 
1327  // Set variables of the root tree
1328  // If we got a cylindricalPET system
1329  #ifdef CASTOR_ROOT
1330  if(GATE_system_type == GATE_SYS_CYLINDRICAL)
1331  {
1332  GEvents[iFic]->SetBranchAddress("time1",&time1);
1333  GEvents[iFic]->SetBranchAddress("rsectorID1",&rsectorID1);
1334  GEvents[iFic]->SetBranchAddress("moduleID1",&moduleID1);
1335  GEvents[iFic]->SetBranchAddress("submoduleID1",&submoduleID1);
1336  GEvents[iFic]->SetBranchAddress("crystalID1",&crystalID1);
1337  GEvents[iFic]->SetBranchAddress("layerID1", &layerID1);
1338  GEvents[iFic]->SetBranchAddress("comptonPhantom1",&comptonPhantomID1);
1339  GEvents[iFic]->SetBranchAddress("RayleighPhantom1",&rayleighPhantomID1);
1340  GEvents[iFic]->SetBranchAddress("eventID1",&eventID1);
1341  GEvents[iFic]->SetBranchAddress("sourceID1",&sourceID1);
1342 
1343  GEvents[iFic]->SetBranchAddress("time2",&time2);
1344  GEvents[iFic]->SetBranchAddress("rsectorID2",&rsectorID2);
1345  GEvents[iFic]->SetBranchAddress("moduleID2",&moduleID2);
1346  GEvents[iFic]->SetBranchAddress("submoduleID2",&submoduleID2);
1347  GEvents[iFic]->SetBranchAddress("crystalID2",&crystalID2);
1348  GEvents[iFic]->SetBranchAddress("layerID2", &layerID2);
1349  GEvents[iFic]->SetBranchAddress("comptonPhantom2",&comptonPhantomID2);
1350  GEvents[iFic]->SetBranchAddress("RayleighPhantom2",&rayleighPhantomID2);
1351  GEvents[iFic]->SetBranchAddress("eventID2",&eventID2);
1352  GEvents[iFic]->SetBranchAddress("sourceID2",&sourceID2);
1353 
1354  GEvents[iFic]->SetBranchAddress("sourcePosX1",&sourcePosX1);
1355  GEvents[iFic]->SetBranchAddress("sourcePosY1",&sourcePosY1);
1356  GEvents[iFic]->SetBranchAddress("sourcePosZ1",&sourcePosZ1);
1357  }
1358  else if(GATE_system_type == GATE_SYS_SPECT)
1359  {
1360  GEvents[iFic]->SetBranchAddress("time",&time1);
1361  GEvents[iFic]->SetBranchAddress("headID", &headID);
1362  GEvents[iFic]->SetBranchAddress("crystalID",&crystalID1);
1363  GEvents[iFic]->SetBranchAddress("pixelID", &pixelID);
1364  GEvents[iFic]->SetBranchAddress("rotationAngle", &rotAngle);
1365  GEvents[iFic]->SetBranchAddress("globalPosX",&gPosX);
1366  GEvents[iFic]->SetBranchAddress("globalPosY",&gPosY);
1367  GEvents[iFic]->SetBranchAddress("globalPosZ",&gPosZ);
1368  GEvents[iFic]->SetBranchAddress("comptonPhantom",&comptonPhantomID1);
1369  GEvents[iFic]->SetBranchAddress("RayleighPhantom",&rayleighPhantomID1);
1370  GEvents[iFic]->SetBranchAddress("sourcePosX",&sourcePosX1);
1371  GEvents[iFic]->SetBranchAddress("sourcePosY",&sourcePosY1);
1372  GEvents[iFic]->SetBranchAddress("sourcePosZ",&sourcePosZ1);
1373  }
1374  else
1375  {
1376  GEvents[iFic]->SetBranchAddress("time1",&time1);
1377  GEvents[iFic]->SetBranchAddress("blockID1",&blockID1);
1378  GEvents[iFic]->SetBranchAddress("crystalID1",&crystalID1);
1379  GEvents[iFic]->SetBranchAddress("comptonPhantom1",&comptonPhantomID1);
1380  GEvents[iFic]->SetBranchAddress("RayleighPhantom1",&rayleighPhantomID1);
1381  GEvents[iFic]->SetBranchAddress("eventID1",&eventID1);
1382  GEvents[iFic]->SetBranchAddress("sourceID1",&sourceID1);
1383 
1384  GEvents[iFic]->SetBranchAddress("time2",&time2);
1385  GEvents[iFic]->SetBranchAddress("blockID2",&blockID2);
1386  GEvents[iFic]->SetBranchAddress("crystalID2",&crystalID2);
1387  GEvents[iFic]->SetBranchAddress("comptonPhantom2",&comptonPhantomID2);
1388  GEvents[iFic]->SetBranchAddress("RayleighPhantom2",&rayleighPhantomID2);
1389  GEvents[iFic]->SetBranchAddress("eventID2",&eventID2);
1390  GEvents[iFic]->SetBranchAddress("sourceID2",&sourceID2);
1391  GEvents[iFic]->SetBranchAddress("sourcePosX1",&sourcePosX1);
1392  GEvents[iFic]->SetBranchAddress("sourcePosY1",&sourcePosY1);
1393  GEvents[iFic]->SetBranchAddress("sourcePosZ1",&sourcePosZ1);
1394  }
1395 
1396 
1397 
1398 
1399  //---------------------------------------------- /ROOT data variables
1400 
1401  // Loop on the GEvents in the current datafile
1402  for (int i=0; i<GEvents[iFic]->GetEntries() ; i++)
1403  {
1404  GEvents[iFic]->GetEntry(i);
1405  printing_index++;
1406 
1407  if(vb >= 4)
1408  Cout("File#" << iFic << ", event#" << i << endl;);
1409 
1410  // ID Conversions
1411  if (GATE_system_type == GATE_SYS_CYLINDRICAL)
1412  {
1413  if(vb >= 4)
1414  {
1415  Cout("Crystal 1 : RsectorID: " << rsectorID1 << ", moduleID: " << moduleID1 << ", submoduleID: " << submoduleID1 << ", crystalID: " << crystalID1
1416  << ", layerID: " << layerID1 << endl;);
1417  Cout("Crystal 2 :RsectorID: " << rsectorID2 << ", moduleID2: " << moduleID2 << ", submoduleID: " << submoduleID2 << ", crystalID: " << crystalID2
1418  << ", layerID: " << layerID2 << endl;);
1419  }
1420 
1421  castorID1 = ConvertIDcylindrical(nRsectorsAngPos, nRsectorsAxial, rsector_id_order,
1422  nModulesTransaxial, nModulesAxial,
1423  nSubmodulesTransaxial, nSubmodulesAxial,
1424  nCrystalsTransaxial, nCrystalsAxial,
1425  nLayers, nb_crystal_per_layer,
1426  nLayersRptTransaxial, nLayersRptAxial,
1427  layerID1, crystalID1, submoduleID1, moduleID1, rsectorID1);
1428 
1429  castorID2 = ConvertIDcylindrical(nRsectorsAngPos, nRsectorsAxial, rsector_id_order,
1430  nModulesTransaxial, nModulesAxial,
1431  nSubmodulesTransaxial, nSubmodulesAxial,
1432  nCrystalsTransaxial, nCrystalsAxial,
1433  nLayers, nb_crystal_per_layer,
1434  nLayersRptTransaxial, nLayersRptAxial,
1435  layerID2, crystalID2, submoduleID2, moduleID2, rsectorID2);
1436 
1437  if(vb >= 4)
1438  {
1439  Cout("--> castorID1: " << castorID1 << endl;);
1440  Cout("--> castorID2: " << castorID2 << endl;);
1441  }
1442  }
1443 
1444  else if (GATE_system_type == GATE_SYS_SPECT)
1445  {
1446  if(vb >= 4)
1447  {
1448  Cout("Projection ID : headID: " << headID << ", rotation angle (deg): " << rotAngle << endl;);
1449  Cout("Crystal ID : crystalID: " << crystalID1 << ", pixelID: " << pixelID << ", globalPosX " << gPosX << ", globalPosY " << gPosY << ", globalPosZ " << gPosZ << endl;);
1450  }
1451 
1452  // Get projection ID
1453  castorID1 = ConvertIDSPECTRoot1(headID,
1454  rotAngle,
1455  headAngStepDeg,
1456  nProjectionsByHead);
1457 
1458  // Get crystal ID
1459  castorID2 = ConvertIDSPECTRoot2(nbSimulatedPixels,
1460  nPixelsTransaxial,
1461  nPixelsAxial,
1462  headID,
1463  crystalID1,
1464  pixelID,
1465  rotAngle,
1466  headAngPitchDeg,
1467  crystalSizeAxl,
1468  crystalSizeTrs,
1469  gPosX,
1470  gPosY,
1471  gPosZ);
1472 
1473  if(vb >= 4)
1474  {
1475  Cout("--> castorID1: " << castorID1 << endl;);
1476  Cout("--> castorID2: " << castorID2 << endl;);
1477  }
1478 
1479  }
1480 
1481  else if (GATE_system_type == GATE_SYS_ECAT)
1482  {
1483  if(vb >= 4)
1484  {
1485  Cout("Crystal 1 : BlockID: " << blockID1 << ", crystalID: " << crystalID1 << endl;);
1486  Cout("Crystal 2 : BlockID: " << blockID2 << ", crystalID: " << crystalID2 << endl;);
1487  }
1488 
1489  castorID1 = ConvertIDecat(nBlocksPerRing, nBlocksLine, nCrystalsTransaxial, nCrystalsAxial, crystalID1, blockID1);
1490  castorID2 = ConvertIDecat(nBlocksPerRing, nBlocksLine, nCrystalsTransaxial, nCrystalsAxial, crystalID2, blockID2);
1491 
1492  if(vb >= 4)
1493  {
1494  Cout("--> castorID1: " << castorID1 << endl;);
1495  Cout("--> castorID2: " << castorID2 << endl;);
1496  }
1497 
1498  }
1499 
1500  // Find out the kind of coincidence (true, scatter, random)
1501  kind = ComputeKindGATEEvent(eventID1, eventID2, comptonPhantomID1, comptonPhantomID2, rayleighPhantomID1, rayleighPhantomID2);
1502 
1503  // Count nb LORs according to kind
1504  switch (kind)
1505  {
1506  case 1:
1507  nLORs_trues++;
1508  break;
1509 
1510  case 2:
1511  nLORs_scatters++;
1512  break;
1513 
1514  case 3:
1515  nLORs_scatters++;
1516  break;
1517 
1518  case 4:
1519  nLORs_rdms++;
1520  break;
1521 
1522  default:
1523  nLORs_unknown++;
1524  }
1525 
1526  // Skip next step if event is not true if we only recover true
1527  if (true_only_flag==true && kind != KIND_TRUE)
1528  continue;
1529 
1530 
1531  // --- Write Event ---
1532 
1533  // SPECT event
1534  if(GATE_system_type == GATE_SYS_SPECT)
1535  {
1536  // HISTOGRAM mode
1537  if(histo_out_flag)
1538  p_bins[castorID1][castorID2]++;
1539  else
1540  // LIST-mode
1541  {
1542  // Write event in the datafile
1543  uint32_t time_in_ms = time1*1000;
1544  Event->SetNbLines(1);
1545  Event->SetID1(0, castorID1); // 1st argument : line, 2nd argument :index
1546  Event->SetID2(0, castorID2); // 1st argument : line, 2nd argument :index
1547  Event->SetTimeInMs(time_in_ms);
1548 
1549  ((iEventListSPECT*)Event)->SetKind(kind);
1550 
1551  Out_data_file->WriteEvent(Event, 0);
1552  }
1553  }
1554  else // PET event
1555  {
1556  // HISTOGRAM mode
1557  if(histo_out_flag)
1558  {
1559  if (castorID1 < castorID2)
1560  p_bins[castorID1][castorID2-castorID1-1]++;
1561  else
1562  p_bins[castorID2][castorID1-castorID2-1]++;
1563  }
1564  // LIST-mode
1565  else
1566  {
1567  // Write event in the datafile
1568  uint32_t time_in_ms = time1*1000;
1569  int nb_lines_in_event = 1; // 1 by default for GATE root files
1570 
1571  Event->SetNbLines(nb_lines_in_event);
1572  Event->SetID1(0, castorID1); // 1st argument : line, 2nd argument :index
1573  Event->SetID2(0, castorID2); // 1st argument : line, 2nd argument :index
1574  Event->SetTimeInMs(time_in_ms);
1575 
1576  ((iEventListPET*)Event)->SetKind(kind);
1577 
1578  Out_data_file->WriteEvent(Event, 0);
1579  }
1580  }
1581 
1582 
1583  // Source image (if no rdm event)
1584  if(src_img_flag &&
1585  eventID1 == eventID2)
1586  {
1587  INTNB x,y,z;
1588  x = (int)(( sourcePosX1 + dim_src_img[0]/2*vox_src_img[0]) / vox_src_img[0]);
1589  // CASToR Y axis inverted in comparison with GATE (left-handed coordinate)
1590  y = (int)(( sourcePosY1 + dim_src_img[1]/2*vox_src_img[1]) / vox_src_img[1]);
1591  z = (int)(( sourcePosZ1 + dim_src_img[2]/2*vox_src_img[2]) / vox_src_img[2]);
1592 
1593  if(x >= 0 && x < dim_src_img[0] &&
1594  y >= 0 && y < dim_src_img[1] &&
1595  z >= 0 && z < dim_src_img[2] )
1596  p_src_img[z*dim_src_img[1]*dim_src_img[0] + y*dim_src_img[0] + x]++;
1597  }
1598 
1599 
1600  // Progression
1601  if (printing_index%(nLORs_tot/printing_ratio) == 0)
1602  {
1603  FLTNB percent = ( ((FLTNB)(printing_index+1))/((FLTNB)nLORs_tot) ) * ((FLTNB)100);
1604  cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b "
1605  << percent << " % ";
1606  }
1607 
1608 
1609  } // end of loop on events in input datafile iFic
1610  #endif
1611 
1612  if(vb>=2) Cout(endl << "DataFile " << iFic << " : " << path_to_input_file[iFic] << " process OK" << endl);
1613 
1614  } // end of loop on input datafiles
1615 
1616  // Free Root objects
1617  #ifdef CASTOR_ROOT
1618  delete[] Tfile_root;
1619  delete[] GEvents;
1620  delete Tapp;
1621  #endif
1622 
1623 
1624 
1625  // ============================================================================================================
1626  // Write the HISTOGRAM datafile and header
1627  // ============================================================================================================
1628  if(GATE_system_type == GATE_SYS_SPECT)
1629  {
1630  if (histo_out_flag)
1631  {
1632  printing_index=0;
1633 
1634  // Writing the datafile
1635  Cout(endl << "Generate the histogram datafile" << endl);
1636 
1637  uint32_t nb_bins = bins_elts*nCrystalsTot;
1638  printing_ratio = (nb_bins>1000) ? 1000 : nb_bins/10;
1639 
1640  // Loop on the crystal IDs
1641  for (size_t id1=0 ; id1<bins_elts ; id1++)
1642  for (size_t id2=0; id2<nCrystalsTot ; id2++)
1643  {
1644  uint32_t nb_events_in_bin = p_bins[id1][id2];
1645 
1646  int nb_lines = 1; // 1 by default for GATE root file;
1647  Event->SetNbLines(nb_lines);
1648  Event->SetID1(0, id1); // 1st argument : line, 2nd argument :index
1649  Event->SetID2(0, id2); // 1st argument : line, 2nd argument :index
1650  Event->SetEventValue(0, nb_events_in_bin);
1651 
1652  // Write event in DataFile
1653  Out_data_file->WriteEvent(Event, 0);
1654  nBINs++;
1655 
1656  // Progression
1657 
1658  if (printing_index%((nb_bins)/printing_ratio) == 0)
1659  {
1660  FLTNB percent = ( ((FLTNB)(printing_index+1))/((FLTNB)nb_bins) ) * ((FLTNB)100);
1661  cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b "
1662  << percent << " % ";
1663  }
1664 
1665  printing_index++;
1666  }
1667 
1668  Out_data_file->SetStartTime((FLTNB)(start_time_ms)/1000); // start time in s
1669  Out_data_file->SetDuration((FLTNB)(duration_ms)/1000); // duration in s
1670  Out_data_file->SetNbEvents( nBINs );
1671 
1672  Cout(endl << "The output histogram contains " << nBINs << " events." << endl);
1673  }
1674  // Write the List-mode datafile header
1675  else
1676  {
1677 
1678  Out_data_file->SetStartTime((FLTNB)(start_time_ms)/1000); // start time in s
1679  Out_data_file->SetDuration((FLTNB)(duration_ms)/1000); // duration in s
1680  //uint32_t nLORs = (true_only_flag) ? nLORs_trues : nLORs_tot;
1681  Out_data_file->SetNbEvents( (true_only_flag) ?
1682  nLORs_trues :
1683  nLORs_tot ) ;
1684  }
1685 
1686  FLTNB* p_angles = new FLTNB[nProjectionsTot];
1687  FLTNB* p_distToDetector = new FLTNB[nProjectionsTot];
1688 
1689  for(size_t h=0 ; h<nHeads ; h++)
1690  for(size_t p=0 ; p<nProjectionsByHead ; p++)
1691  {
1692  int idx_proj = h*nProjectionsByHead+p;
1693  p_angles[idx_proj] = head1stAngleDeg // initial angular position of the head
1694  + p*headAngStepDeg // projection angular position
1695  + h*headAngPitchDeg; // angular shift for each head
1696 
1697  // same distance for each head with GATE
1698  p_distToDetector[idx_proj] = distToDetector;
1699  }
1700 
1701  ((iDataFileSPECT*)Out_data_file)->SetNbBins(nPixelsTransaxial, nPixelsAxial);
1702  ((iDataFileSPECT*)Out_data_file)->SetNbProjections(nProjectionsTot);
1703 
1704  if( ((iDataFileSPECT*)Out_data_file)->InitAngles(p_angles) )
1705  {
1706  Cerr("**** castor-GATERootToCastor :: Error when trying to set projection angles values !" << endl);
1707  Exit(EXIT_FAILURE);
1708  }
1709 
1710  if( ((iDataFileSPECT*)Out_data_file)->InitCorToDetectorDistance(p_distToDetector) )
1711  {
1712  Cerr("**** castor-GATERootToCastor :: Error when trying to set distance between center of rotation and detectors !" << endl);
1713  Exit(EXIT_FAILURE);
1714  }
1715 
1716  ((iDataFileSPECT*)Out_data_file)->SetHeadRotDirection(headRotDirection);
1717 
1718  delete[] p_angles;
1719  delete[] p_distToDetector;
1720 
1721  Out_data_file->WriteHeader();
1722  } // end of SPECT section
1723 
1724  else
1725  {
1726  if (histo_out_flag)
1727  {
1728  printing_index=0;
1729 
1730  // Writing the datafile
1731  Cout(endl << "Generate the histogram datafile" << endl);
1732  uint32_t nb_bins = (nCrystalsTot*nCrystalsTot - nCrystalsTot)/2;
1733  printing_ratio = (nb_bins>1000) ? 1000 : nb_bins/10;
1734 
1735  // Loop on the crystal IDs
1736  for (size_t id1=0 ; id1 <bins_elts ; id1++)
1737  for (size_t id2 = id1+1; id2 < nCrystalsTot;id2++)
1738  {
1739  uint32_t nb_events_in_bin = p_bins[id1][id2-id1-1];
1740 
1741  int nb_lines = 1; // 1 by default for GATE root file;
1742  Event->SetNbLines(nb_lines);
1743  Event->SetID1(0, id1); // 1st argument : line, 2nd argument :index
1744  Event->SetID2(0, id2); // 1st argument : line, 2nd argument :index
1745  Event->SetEventValue(0, nb_events_in_bin);
1746 
1747  // Estimate acf for the bin
1748  if(estimate_acf_flag)
1749  {
1750  // Compute the system matrix elements for the two crystals
1751  oProjectionLine* line = p_ProjectorManager->ComputeProjectionLine(Event, 0);
1752 
1753  // Compute forward projection of attenuation image
1754  FLTNB fp = 0.;
1755  if (line->NotEmptyLine())
1756  fp = line->ForwardProject(p_ImageSpace->m4p_image[0][0][0]) * 0.1; // 0.1 -> convert in mm-1
1757 
1758  // Write atn correction factor in Event
1759  ((iEventPET*)Event)->SetAttenuationCorrectionFactor(1/std::exp(-fp));
1760  }
1761 
1762  // Write event in DataFile
1763  Out_data_file->WriteEvent(Event, 0);
1764  nBINs++;
1765 
1766  // Progression
1767  if (printing_index%((nb_bins)/printing_ratio) == 0)
1768  {
1769  FLTNB percent = ( ((FLTNB)(printing_index+1))/((FLTNB)nb_bins) ) * ((FLTNB)100);
1770  cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b "
1771  << percent << " % ";
1772  }
1773  printing_index++;
1774  }
1775 
1776  Out_data_file->SetStartTime((FLTNB)(start_time_ms)/1000); // start time in s
1777  Out_data_file->SetDuration((FLTNB)(duration_ms)/1000); // duration in s
1778  Out_data_file->SetNbEvents( nBINs );
1779  Out_data_file->WriteHeader();
1780 
1781  Cout(endl << "The output histogram contains " << nBINs << " events." << endl);
1782  }
1783  // Write the List-mode datafile header
1784  else
1785  {
1786  Out_data_file->SetStartTime((FLTNB)(start_time_ms)/1000); // start time in s
1787  Out_data_file->SetDuration((FLTNB)(duration_ms)/1000); // duration in s
1788  //uint32_t nLORs = (true_only_flag) ? nLORs_trues : nLORs_tot;
1789  Out_data_file->SetNbEvents( (true_only_flag) ?
1790  nLORs_trues :
1791  nLORs_tot );
1792  Out_data_file->WriteHeader();
1793  }
1794  } // end of PET section
1795 
1796 
1797  // Write the source image
1798  if(src_img_flag)
1799  {
1800  if (vb>=2) Cout("Writing source image ..." << endl);
1801 
1802  // Initialize Intf_fields object with the source image dimensions
1803  Intf_fields IF;
1804  IntfKeyInitFields(&IF);
1805 
1806  IF.mtx_size[0] = dim_src_img[0];
1807  IF.mtx_size[1] = dim_src_img[1];
1808  IF.mtx_size[2] = dim_src_img[2];
1809  IF.vox_size[0] = vox_src_img[0];
1810  IF.vox_size[1] = vox_src_img[1];
1811  IF.vox_size[2] = vox_src_img[2];
1812  IF.nb_total_imgs = dim_src_img[2];
1813 
1814  if (IntfWriteImgFile(path_to_src_image.append("_src"), p_src_img, IF, vb) )
1815  {
1816  Cerr("***** castor-GATERootToCastor :: Error writing Interfile of output image !" << endl);
1817  Exit(EXIT_FAILURE);
1818  }
1819  }
1820 
1821  Cout(endl << "The simulated dataset contained " << nLORs_tot << " coincidences/singles, with: " << endl
1822  << " " << nLORs_trues << " trues (" << 100*nLORs_trues/nLORs_tot << " %), " << endl
1823  << " " << nLORs_scatters << " scatters (" << 100*nLORs_scatters/nLORs_tot << " %)," << endl
1824  << " " << nLORs_rdms << " randoms (" << 100*nLORs_rdms/nLORs_tot << " %)," << endl
1825  << " " << nLORs_unknown << " unknown (" << 100*nLORs_unknown/nLORs_tot << " %)." << endl);
1826 
1827  if (vb>=2) Cout("Writing raw datafile ..." << endl); // todo just change name if only one datafile
1828 
1829 
1830  Out_data_file->PROJ_WriteData();
1831  Out_data_file->PROJ_DeleteTmpDataFile();
1832 
1833  Cout(endl << " --- Conversion completed --- " << endl);
1834 
1835  // ============================================================================================================
1836  // End
1837  // ============================================================================================================
1838 
1839  if (vb>=2) Cout(" Deallocating objects ..." << endl);
1840 
1841  // Delete objects
1842  if(histo_out_flag)
1843  {
1844  for(size_t b=0; b<bins_elts ; b++)
1845  delete[] p_bins[b];
1846 
1847  delete[]p_bins;
1848 
1849  if(Event)
1850  {
1851  if(GATE_system_type == GATE_SYS_SPECT)
1852  delete (iEventHistoSPECT*)Event;
1853  else
1854  delete (iEventHistoPET*)Event;
1855  }
1856  }
1857  else
1858  if(Event)
1859  {
1860  if(GATE_system_type == GATE_SYS_SPECT)
1861  delete (iEventListSPECT*)Event;
1862  else
1863  delete (iEventListPET*)Event;
1864  }
1865 
1866  delete[]offset_time_list;
1867  delete Out_data_file;
1868 
1869  if(src_img_flag && p_src_img)
1870  delete[] p_src_img;
1871 
1872  // Free objects created for analytic projection (acf estimation)
1873  if(estimate_acf_flag)
1874  p_ImageSpace->LMS_DeallocateImage();
1875  if(p_ImageSpace) delete p_ImageSpace;
1876  if(p_ProjectorManager) delete p_ProjectorManager;
1877  if(p_ID) delete p_ID;
1878 
1879  Cout(" --- END --- " << endl << endl);
1880 
1881  return 0;
1882 }
void SetDataMode(int a_dataMode)
set the data mode
Definition: vDataFile.hh:406
uint32_t ConvertIDSPECTRoot1(int32_t a_headID, float_t a_rotAngle, float_t a_angStep, uint32_t a_nProjectionsByHead)
Compute a CASToR projection index of a GATE SPECThead system.
This class is designed to be a mother virtual class for DataFile.
Definition: vDataFile.hh:103
void LMS_InstantiateImage()
Allocate memory for the main image matrices (for list-mode sensitivity generation) ...
This header file is mainly used to declare some macro definitions and all includes needed from the st...
int SetPETIsotope(int a_bed, const string &a_isotope)
Set the PET isotope for the provided bed.
static sScannerManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
int InitDynamicData(string a_pathTo4DDataSplittingFile, int a_respMotionCorrectionFlag, int a_cardMotionCorrectionFlag, int a_doubleMotionCorrectionFlag, int a_invMotionCorrectionFlag, int a_nbRespGates, int a_nbCardGates)
Call the eponym function from the oDynamicDataManager object in order to initialize its data...
#define MODE_HISTOGRAM
Definition: vDataFile.hh:59
#define TYPE_PET
Definition: vDataFile.hh:74
void SetNbLines(uint16_t a_value)
Set the number of lines of the Event.
Definition: vEvent.hh:154
void SetVerbose(int a_verbose)
set verbosity
void SetAtnCorrectionFlagOn()
set to true the flag indicating the presence of attenuation correction factors in the datafile ...
void SetFOVSizeZ(FLTNB a_fovSizeZ)
Set the FOV's size along the Z axis, in mm.
#define MODE_LIST
Definition: vDataFile.hh:57
#define GATE_SYS_UNKNOWN
FLTNB vox_size[3]
#define FLTNB
Definition: gVariables.hh:81
This file gathers various function dedicated to data conversion in order to convert various type of G...
uint32_t nb_total_imgs
void SetComputationStrategy(int a_computationStrategy)
Set the computation strategy for the system matrix elements storage.
int Initialize()
A function used to initialize the manager and the projectors or system matrices it manages...
virtual int WriteHeader()=0
This function is implemented in child classes. Generate a header file according to the data output ...
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
Set the image dimensions in use.
int FindScannerSystem(string a_scannerName)
Look for a file matching with the scanner name in parameter inside the scanner repository.
void SetFOVSizeY(FLTNB a_fovSizeY)
Set the FOV's size along the Y axis, in mm.
void SetNbVoxZ(INTNB a_nbVoxZ)
Set the number of voxels along the Z axis.
int ReadMacSPECT(string a_pathMac, float_t &distToDetector, uint32_t &nHeads, uint32_t &nPixAxl, uint32_t &nPixTrs, float_t &crystalSizeAxl, float_t &crystalSizeTrs, uint32_t &nProjectionsTot, uint32_t &nProjectionsByHead, float_t &head1stAngle, float_t &headAngPitch, float_t &headAngStepDeg, int &headRotDirection, uint32_t &start_time_ms, uint32_t &duration_ms, int vb)
Recover informations about the scanner element of an ECAT system, and acquisition duration...
int BuildScannerObject()
Instantiate the specific scanner object related to the modality, and set verbosity of scanner object...
#define GATE_SYS_ECAT
void SetNbRespGates(int a_nbRespGates)
Set the number of respiratory gates.
#define KIND_TRUE
void SetOffsetY(FLTNB a_offsetY)
Set the image offset along the Y axis, in mm.
virtual int WriteEvent(vEvent *ap_Event, int a_th=0)=0
This function is implemented in child classes Write event according to the chosen type of data...
void ShowHelp(int a_returnCode)
void SetFOVOutMasking(FLTNB a_fovOutPercent, INTNB a_nbSliceOutMask)
Set the output FOV masking settings: transaxial unmasked FOV percent and number of extrem slices to r...
Inherit from iEventSPECT. Class for SPECT histogram mode events.
void SetVerbose(int a_verboseLevel)
set verbosity
Definition: vDataFile.hh:439
void SetOptionsForward(const string &a_optionsForward)
Set the forward projection options contained in the provided string.
void SetDuration(FLTNB a_value)
Definition: vDataFile.hh:516
virtual int ComputeSizeEvent()=0
This function is implemented in child classes Computation of the size of each event according to th...
#define GATE_SYS_SPECT
virtual void SetEventValue(int a_bin, FLTNBDATA a_value)=0
This function is implemented by child classes.
Definition: vEvent.hh:140
static sOutputManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
void SetOffsetZ(FLTNB a_offsetZ)
Set the image offset along the Z axis, in mm.
#define TYPE_SPECT
Definition: vDataFile.hh:76
Declaration of class iDataFilePET.
void SetOptionsCommon(const string &a_optionsCommon)
Set the common projection options contained in the provided string.
int CheckParameters()
Check if all parameters have been correctly initialized, and call the CheckParameters function of the...
void SetNbVoxX(INTNB a_nbVoxX)
Set the number of voxels along the X axis.
int SetNbThreads(const string &a_nbThreads)
Set the number of threads.
uint32_t ConvertIDcylindrical(uint32_t nRsectorsAngPos, uint32_t nRsectorsAxial, int a_rsectorIdOrder, uint32_t nModulesTransaxial, uint32_t nModulesAxial, uint32_t nSubmodulesTransaxial, uint32_t nSubmodulesAxial, uint32_t nCrystalsTransaxial, uint32_t nCrystalsAxial, uint8_t nLayers, uint32_t *nCrystalPerLayer, vector< uint32_t > nLayersRptTransaxial, vector< uint32_t > nLayersRptAxial, int32_t layerID, int32_t crystalID, int32_t submoduleID, int32_t moduleID, int32_t rsectorID)
Compute a CASToR crystal index of a GATE cylindricalPET system from its indexes (rsector/module/submo...
Declaration of class iDataFileSPECT.
void Exit(int code)
void SetOptionsBackward(const string &a_optionsBackward)
Set the backward projection options contained in the provided string.
int IntfReadHeader(const string &a_pathToHeaderFile, Intf_fields *ap_IntfFields, int vb)
Read an Interfile header.
int InstantiateScanner()
Instantiate scanner using the related function in the scanner classes.
int LogCommandLine(int argc, char **argv)
Write log file header with the provided command line options and different informations.
int PROJ_LoadInitialImage(const string &a_pathToImage)
Load the initial image for the analytical projection.
void SetIsotope(string a_value)
initialize the isotope string value
int PROJ_WriteData()
Write/Merge chunk of data in a general data file.
Definition: vDataFile.cc:715
int ConvertFromString(const string &a_str, string *a_result)
Copy the 'a_str' string in the position pointed by 'a_result'.
Definition: gOptions.cc:771
void SetIsotope(string a_value)
initialize the isotope string value
#define FIXED_LIST_COMPUTATION_STRATEGY
void SetNbVoxY(INTNB a_nbVoxY)
Set the number of voxels along the Y axis.
int CheckConfigDir(const string &a_path)
Set the path to the CASTOR config directory from the given path if not empty or through the existence...
Inherit from iEventPET. Class for PET list-mode events.
int CheckParameters()
A function used to check the parameters settings.
void SetNbCardGates(int a_nbCardGates)
Set the number of cardiac gates.
#define Cerr(MESSAGE)
int ReadMacCylindrical(string a_pathMac, uint8_t &nLayers, uint32_t *nb_crystal_per_layer, uint32_t &nCrystalsTot, uint32_t &nCrystalsAxial, uint32_t &nCrystalsTransaxial, vector< uint32_t > &nLayersRptAxial, vector< uint32_t > &nLayersRptTransaxial, uint32_t &nSubmodulesAxial, uint32_t &nSubmodulesTransaxial, uint32_t &nModulesAxial, uint32_t &nModulesTransaxial, uint32_t &nRsectorsAxial, uint32_t &nRsectorsAngPos, int &rsector_id_order, uint32_t &start_time_ms, uint32_t &duration_ms, int vb)
Recover informations about the scanner element of a cylindricalPET system and acquisition duration...
Inherit from iEventPET. Class for PET histogram mode events.
FLTNB **** m4p_image
Definition: oImageSpace.hh:81
int BuildLUT()
Call the eponym function of the scanner class.
int ComputeKindGATEEvent(uint32_t eventID1, uint32_t eventID2, int comptonPhantom1, int comptonPhantom2, int rayleighPhantom1, int rayleighPhantom2)
Determine kind of a given coincidence event, from its attributes.
void SetStartTime(FLTNB a_value)
Definition: vDataFile.hh:508
int Initialize()
A function used to initialize all that is needed.
Singleton class that manages output writing on disk (images, sinograms, etc). It also manages loggi...
const string & GetPathToConfigDir()
Return the path to the CASTOR config directory.
void SetFrames(const string &a_frameList)
Set the frame list (a string that will be parsed by the InitializeFramingAndQuantification function) ...
void SetFOVSizeX(FLTNB a_fovSizeX)
Set the FOV's size along the X axis, in mm.
void SetDataType(int a_dataType)
set the data type
Definition: vDataFile.hh:413
Singleton class that Instantiate and initialize the scanner object.
int IntfReadProjectionImage(const string &a_pathToHeaderFile, FLTNB *ap_ImgMatrix, Intf_fields *ap_IF, int vb, bool a_lerpFlag)
Main function which reads a projection Interfile 3D projection image and store its content in the pro...
void LMS_DeallocateImage()
Free memory for the main image matrices (for list-mode sensitivity generation)
int ReadMacECAT(string a_pathMac, uint32_t &nCrystalsTot, uint32_t &nCrystalsAxial, uint32_t &nCrystalsTransaxial, uint32_t &nBlocksLine, uint32_t &nBlocksPerRing, uint32_t &start_time_ms, uint32_t &duration_ms, int vb)
Recover informations about the scanner element of an ECAT system and acquisition duration, from a GATE macro file.
Declaration of class sScannerManager.
void SetVoxSizeY(FLTNB a_voxSizeY)
Set the voxel's size along the Y axis, in mm.
int CreateGeomWithECAT(string a_pathMac, string a_pathGeom)
Read a GATE macro file containing the description of an ecat system, and convert it to a geom file...
virtual int PrepareDataFile()=0
This function is implemented in child classes Store different kind of information inside arrays (da...
oProjectionLine * ComputeProjectionLine(vEvent *ap_Event, int a_th)
This function is used to compute system matrix elements from the associated projector or pre-computed...
Inherit from vDataFile. Class that manages the reading of a SPECT input file (header + data)...
void SetOffsetX(FLTNB a_offsetX)
Set the image offset along the X axis, in mm.
void SetHeaderDataFileName(const string &a_headerFileName)
set the data header file name
Definition: vDataFile.hh:476
vScanner * GetScannerObject()
void SetVoxSizeX(FLTNB a_voxSizeX)
Set the voxel's size along the X axis, in mm.
void IntfKeyInitFields(Intf_fields *ap_IF)
Init the file of an Interfile fields structure passed in parameter to their default values...
int CheckParameters()
A function used to check the parameters settings.
uint32_t mtx_size[7]
void SetMPIRank(int a_mpiRank)
Initialize the machine index for MPI.
void SetVerbose(int a_verboseLevel)
Set the verbose level.
int AllocateID()
Instantiate the mp_ID1 and mp_ID2 indices arrays.
Definition: vEvent.cc:69
bool NotEmptyLine()
This function is used to know if the line contains any voxel contribution.
void SetDataFile(vDataFile *ap_DataFile)
Set a data file in use to later recover some information from it.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
set the pointer to the oImageDimensionsAndQuantification object
Definition: oImageSpace.hh:624
#define GATE_SYS_CYLINDRICAL
Inherit from vEvent. Main PET class for the Event objects.
Definition: iEventPET.hh:42
void SetNbEvents(int64_t a_value)
initialize the number of events with a int64_t value
Definition: vDataFile.hh:500
int Initialize()
Initialization : .
void SetVerbose(int a_verbose)
Set the member m_verboseLevel to the provided value.
void SetID2(int a_line, uint32_t a_value)
Set the indice associated with the line index for the 2nd ID of the Event.
Definition: vEvent.hh:170
#define INTNB
Definition: gVariables.hh:92
#define OS_SEP
This class is designed to manage and store system matrix elements associated to a vEvent...
string GetPathOfFile(const string &a_pathToFile)
Simply return the path to the directory of a file path string passed in parameter.
Definition: gOptions.cc:1166
void SetVoxSizeZ(FLTNB a_voxSizeZ)
Set the voxel's size along the Z axis, in mm.
Declaration of class oImageSpace.
void SetVerbose(int a_verboseLevel)
set verbosity
Definition: oImageSpace.hh:617
This class is designed to manage the projection part of the reconstruction.
Interfile fields. This structure contains all the Interfile keys currently managed by CASToR Decl...
void SetTimeInMs(uint32_t a_value)
Set the timestamp of the Event.
Definition: vEvent.hh:139
Declaration of class sOutputManager.
This class holds all the matrices in the image domain that can be used in the algorithm: image...
Definition: oImageSpace.hh:61
Mother class for the Event objects.
Definition: vEvent.hh:43
int InitOutputDirectory(const string &a_pathFout, const string &a_pathDout)
Create the output directory if any, extract the base name and create the log file.
FLTNB ForwardProject(FLTNB *ap_image=NULL)
Simply forward projects the provided image if not null, or else 1, for the current TOF bin...
int main(int argc, char **argv)
This class is designed to manage all dimensions and quantification related stuff. ...
void SetNbBeds(int a_nbBeds)
Set the number of bed positions and allocate the bed positions if not already done.
This file is used for all kind of different functions designed for options parsing and ASCII file rea...
void SetVerbose(int a_verboseLevel)
set verbosity
#define Cout(MESSAGE)
uint32_t ConvertIDecat(int32_t nBlocksPerRing, int32_t nBlocksLine, int32_t nCrystalsTransaxial, int32_t nCrystalsAxial, int32_t crystalID, int32_t blockID)
Compute a CASToR crystal index of a GATE ecat system from its indexes (block/crystal) and the system ...
int PROJ_DeleteTmpDataFile()
Delete temporary datafile used for multithreaded output writing if needed.
Definition: vDataFile.cc:786
#define CASTOR_VERSION
Definition: gVariables.hh:70
void SetScanner(vScanner *ap_Scanner)
Set the scanner in use.
int ReadIntfSPECT(string a_pathIntf, float_t &ap_distToDetector, uint32_t &ap_nHeads, uint32_t &ap_nPixAxl, uint32_t &ap_nPixTrs, float_t &ap_crystalSizeAxl, float_t &ap_crystalSizeTrs, uint32_t &ap_nProjectionsTot, uint32_t &ap_nProjectionsByHead, float_t &ap_head1stAngle, float_t &ap_headAngPitch, float_t &headAngStepDeg, int &ap_headRotDirection, uint32_t &ap_start_time_ms, uint32_t &ap_duration_ms, int vb)
Recover informations about the scanner element of an ECAT system, and acquisition duration...
Declaration of class oProjectorManager.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
set the pointer to the oImageDimensionsAndQuantification object
Definition: vDataFile.hh:453
Inherit from vDataFile. Class that manages the reading of a PET input file (header + data)...
Definition: iDataFilePET.hh:44
int ReadStringOption(const string &a_input, T *ap_return, int a_nbElts, const string &sep, const string &a_option)
Parse the 'a_input' string corresponding to the 'a_option' into 'a_nbElts' elements, using the 'sep' separator. The results are returned in the templated 'ap_return' dynamic templated array. Call "ConvertFromString()" to perform the correct conversion depending on the type of the data to convert.
Definition: gOptions.cc:62
int GetGATESystemType(const string &a_pathMac)
Read a GATE macro file and identify the system type from the 'gate/systems/' command lines...
void SetMaxNumberOfLinesPerEvent(uint16_t a_value)
set the max number of line per event in the datafile
int CreateGeomWithCylindrical(string a_pathMac, string a_pathGeom)
Read a GATE macro file containing the description of a cylindricalPET system, and convert it to a geo...
virtual int PROJ_InitFile()=0
This function is implemented in child classes Initialize the fstream objets for output writing as w...
int IntfWriteImgFile(const string &a_pathToImg, FLTNB *ap_ImgMatrix, const Intf_fields &ap_IntfF, int vb)
Main function dedicated to Interfile 3D image writing. Recover image information from a provided In...
#define GEO_ROT_CW
Definition: vScanner.hh:47
int CreateGeomWithSPECT(string a_pathMac, string a_pathGeom)
Read a GATE macro file containing the description of a SPECThead system, and convert it to a geom fil...
void GetUserEndianness()
Check user/host computer endianness and write it to the global variable User_Endianness.
uint32_t ConvertIDSPECTRoot2(uint32_t a_nbSimulatedPixels, uint32_t a_nPixTrs, uint32_t a_nPixAxl, int32_t a_headID, int32_t a_crystalID, int32_t a_pixelID, float_t a_rotAngle, float_t a_headAngPitch, float_t a_crystalSizeAxl, float_t a_crystalSizeTrs, float_t a_gPosX, float_t a_gPosY, float_t a_gPosZ)
Compute a CASToR crystal index of a GATE SPECThead system.
Inherit from iEventSPECT. Class for SPECT list-mode events.
void SetID1(int a_line, uint32_t a_value)
Set the indice associated with the line index for the 1st ID of the Event.
Definition: vEvent.hh:162