CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00008 #include "gVariables.hh" 00009 #include "gOptions.hh" 00010 #include "iDataFilePET.hh" 00011 #include "iDataFileSPECT.hh" 00012 #include "sOutputManager.hh" 00013 #include "sScannerManager.hh" 00014 #include "gDataConversionUtilities.hh" 00015 #include <iomanip> //std::setprecision 00016 00017 #ifdef CASTOR_ROOT 00018 #include "TROOT.h" 00019 #include "TApplication.h" 00020 #include "TGClient.h" 00021 #include "TCanvas.h" 00022 #include "TSystem.h" 00023 #include "TTree.h" 00024 #include "TBranch.h" 00025 #include "TFile.h" 00026 #endif 00027 00028 00029 #include "oImageSpace.hh" 00030 #include "oProjectorManager.hh" 00031 00032 00038 void ShowHelp(int a_returnCode) 00039 { 00040 cout << endl; 00041 cout << "Usage: castor-GATERootToCastor -i path_to_ifile.root -OR- -il path_to_ifile.txt "<< endl; 00042 cout << " -s scanner_alias "<< endl; 00043 cout << " -o path_to_out_file "<< endl; 00044 cout << " -m path_to_macro.mac " << endl; 00045 cout << " [-t only convert true prompts]" << endl; 00046 cout << " [-oh histogram datafile output]" << endl; 00047 cout << " [-atn path_to_atn_image(cm-1)]" << endl; 00048 cout << " [-k recover coincidence kind]" << endl; 00049 cout << " [-ist isotope_alias" << endl; 00050 cout << " [-ot time offsets in s]" << endl; 00051 cout << " [-vb verbosity]" << endl; 00052 cout << endl; 00053 cout << endl; 00054 cout << "[Main settings]:" << endl; 00055 cout << " -i path_to_file.root : give an input root datafile to convert" << endl; 00056 cout << " -il path_to_file.txt : give an input text file containing a list of root files to convert" << endl; 00057 cout << " : they will be concatenated into 1 CASToR datafile" << endl; 00058 cout << " : the path to each datafile to convert must be entered on a newline" << endl; 00059 cout << " -m path_to_macro.mac : gives the input GATE macro file used for the GATE simulation" << endl; 00060 cout << " -o path_to_out_file : give the path to the output file which will be created inside this folder (no default)" << endl; 00061 cout << " -s scanner_alias : provide the name of the scanner used for to acquire the original data" << endl; 00062 cout << " : Must correspond to a .geom or .hscan file in the config/scanner repository." << endl; 00063 cout << " : A geom file can be created from the macro files using the facultative option -geo (see below)" << endl; 00064 cout << endl; 00065 cout << "[Optional settings]:" << endl; 00066 cout << " -t : only the true prompts will be converted" << endl; 00067 cout << " -oh : Indicate if the output datafile should be written in histogram format (default : List-mode)" << endl; 00068 cout << " -atn path_to_atn_image : For histogram output, provide an attenuation image (cm-1) related to the acquisition." << endl; 00069 cout << " Analytic projection will be performed during the data conversion in order to estimate attenuation correction factors (acf) for each histogram event" << endl; 00070 cout << " -k : For List-mode output, write kind of coincidence (true/scatter/rdm/...) in the output datafile (disabled by default)" << endl; 00071 cout << " -ist isotope_alias : provide alias of the isotope used during the simulation"<< endl; 00072 cout << " Aliases for supported PET isotopes and their parameters are listed in config/misc/isotopes_pet"<< endl; 00073 cout << " Other isotopes can be added in the same file"<< endl; 00074 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; 00075 cout << " The option must start with the path to the output image which will be generated." << endl; 00076 cout << " Dimensions and voxel sizes of the image must be provided using commas as separators, as in the following template:" << endl; 00077 cout << " path/to/image:dimX,dimY,dimZ,voxSizeX,voxSizeY,voxSizeZ"<< endl; 00078 cout << " -geo : Generate a CASToR geom file from the provided GATE macro file(s)"<< endl; 00079 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; 00080 //cout << " -spect_b nbinsT,nbinsA : Option specific to SPECThead simulation with root output."<< endl; 00081 //cout << " Give transaxial and axial number of bins for projections, separated by a comma."<< endl; 00082 //cout << " Pixel sizes will be computed from the crystal surface and the transaxial/axial number of bins" << endl; 00083 cout << " -ot list : Provide a serie of time offset in seconds to apply before each input file"<< endl; 00084 cout << " (In the case of converting several datafiles of a dynamic acquisition, timestamps of events may be resetted for each file" << endl; 00085 cout << " This variable allows to manually increment the time between each datafile(s) if required" << endl; 00086 cout << " The number of time offsets must be equal to the number of input files, provided by -i or -if options." << endl; 00087 cout << " 'list' is a list of time offsets, separated by ','" << endl; 00088 cout << endl; 00089 cout << "[Miscellaneous settings]:" << endl; 00090 cout << " -vb : give the verbosity level, from 0 (no verbose) to above 5 (at the event level) (default: 1)." << endl; 00091 cout << endl; 00092 00093 Exit(a_returnCode); 00094 } 00095 00096 00097 /* 00098 Main program 00099 */ 00100 00101 int main(int argc, char** argv) 00102 { 00103 00104 // No argument, then show help 00105 if (argc==1) ShowHelp(0); 00106 00107 // ============================================================================================================ 00108 // Parameterized variables with their default values 00109 // ============================================================================================================ 00110 00111 // String which gathers the path to the input data filename provided by the user. no default 00112 string input_file = ""; 00113 vector<string> path_to_input_file; 00114 00115 // Path to user provided data and output files/images 00116 string path_to_out_file = ""; 00117 string path_to_data_filename = ""; 00118 string path_to_header_filename = ""; 00119 string path_to_mac_file = ""; 00120 string path_to_atn_image = ""; 00121 string path_to_src_image = ""; 00122 00123 // Scanner name provided by the user 00124 string scanner_name = ""; 00125 00126 // GATE system type 00127 int GATE_system_type = GATE_SYS_UNKNOWN; 00128 00129 // Only recover true GEvents 00130 bool true_only_flag = false; 00131 00132 // Verbosity 00133 int vb = 0; 00134 00135 // Histogram output 00136 bool histo_out_flag = false; 00137 00138 // Estimate acf (histogram output) 00139 bool estimate_acf_flag = false; 00140 00141 // Recover kind of coincidence (list-mode output) 00142 bool kind_flag = false; 00143 00144 // Isotope alias 00145 string istp_alias = ""; 00146 00147 // Variables related to the source image 00148 bool src_img_flag = false; 00149 INTNB dim_src_img[3]; 00150 FLTNB vox_src_img[3]; 00151 FLTNB* p_src_img = NULL; 00152 00153 // Generate geom file 00154 bool geom_flag = false; 00155 00156 // Input is interfile (for SPECT simulation) 00157 bool input_is_intf = false; 00158 00159 // SPECT bins 00160 uint16_t spect_bin_axl = 0, 00161 spect_bin_trs = 0; 00162 00163 // Time offsets 00164 string offset_time_str = ""; 00165 uint32_t* offset_time_list = NULL; 00166 00167 // ============================================================================================================ 00168 // Read command-line parameters 00169 // ============================================================================================================ 00170 for (int i=1; i<argc; i++) 00171 { 00172 string option = (string)argv[i]; 00173 00174 if (option=="-h" || option=="--help" || option=="-help") ShowHelp(0); 00175 00176 // Just one file is provided 00177 if (option=="-i") 00178 { 00179 if (path_to_input_file.size() > 0) 00180 { 00181 Cerr("***** castor-GATERootToCastor :: the following file names have already been provided (-i/-il options must be used ONCE): " << endl); 00182 for (size_t i=0 ; i<path_to_input_file.size() ; i++) 00183 Cout(path_to_input_file[i] << endl); 00184 Exit(EXIT_FAILURE); 00185 } 00186 else 00187 { 00188 if (argv[i+1] == NULL) 00189 { 00190 Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl); 00191 Exit(EXIT_FAILURE); 00192 } 00193 else 00194 { 00195 input_file = argv[i+1]; 00196 path_to_input_file.push_back(input_file); 00197 } 00198 i++; 00199 } 00200 } 00201 00202 // Read a text file containing a list of root datafiles to read 00203 else if (option=="-il") 00204 { 00205 if (path_to_input_file.size() > 0) 00206 { 00207 Cerr("***** castor-GATERootToCastor :: the following file names have already been provided (-i/-il options must be used ONCE) " << endl); 00208 for (size_t i=0 ; i<path_to_input_file.size() ; i++) 00209 Cout(path_to_input_file[i] << endl); 00210 Exit(EXIT_FAILURE); 00211 } 00212 else 00213 { 00214 if (argv[i+1] == NULL) 00215 { 00216 Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl); 00217 Exit(EXIT_FAILURE); 00218 } 00219 else 00220 { 00221 ifstream ifile(argv[i+1], ios::in); 00222 string line; 00223 00224 // Read list of rootfgiles 00225 if(ifile) 00226 { 00227 // Get the position of the list_file, then append the name of the content datafiles to this position. 00228 input_file = argv[i+1]; 00229 00230 // Recover the files 00231 while (getline(ifile, line)) 00232 { 00233 string path_to_datafile = GetPathOfFile(input_file) + OS_SEP + line; 00234 path_to_input_file.push_back(path_to_datafile); 00235 } 00236 } 00237 else 00238 { 00239 Cerr("***** castor-GATERootToCastor :: Error, can't read txt file: " << argv[i+1] << endl); 00240 Exit(EXIT_FAILURE); 00241 } 00242 00243 ifile.close(); 00244 } 00245 i++; 00246 } 00247 } 00248 // Mac file 00249 else if (option=="-m") 00250 { 00251 path_to_mac_file = (string)argv[i+1]; 00252 i++; 00253 } 00254 // Scanner alias 00255 else if (option=="-s") 00256 { 00257 if (argv[i+1] == NULL) 00258 { 00259 Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl); 00260 Exit(EXIT_FAILURE); 00261 } 00262 else 00263 scanner_name = argv[i+1]; 00264 i++; 00265 } 00266 // Output CASToR datafile 00267 else if (option=="-o") 00268 { 00269 if (argv[i+1] == NULL) 00270 { 00271 Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl); 00272 Exit(EXIT_FAILURE); 00273 } 00274 else 00275 path_to_out_file = argv[i+1]; 00276 i++; 00277 } 00278 // Recover only trues 00279 else if (option=="-t") 00280 { 00281 #ifdef CASTOR_ROOT 00282 true_only_flag = true; 00283 #else 00284 Cerr("***** castor-GATERootToCastor :: Option: " << option << " is only available for dataset generated with Gate in a Root format." << endl); 00285 Cerr("***** castor-GATERootToCastor :: Root support is currently not enabled (CASTOR_ROOT environnement variable should be set before compilation)" << endl); 00286 Exit(EXIT_FAILURE); 00287 #endif 00288 } 00289 00290 // Time offsets 00291 else if (option=="-ot") 00292 { 00293 if (i>=argc-1) 00294 { 00295 Cerr("***** castor-GATERootToCastor :: Argument missing for option: " << option << endl); 00296 Exit(EXIT_FAILURE); 00297 } 00298 offset_time_str = (string)argv[i+1]; 00299 i++; 00300 } 00301 00302 // Output datafile in histogram mode 00303 else if (option=="-oh") 00304 { 00305 histo_out_flag = true; 00306 } 00307 00308 else if (option=="-atn") 00309 { 00310 if (i>=argc-1) 00311 { 00312 Cerr("***** castor-GATERootToCastor :: Argument missing for option: " << option << endl); 00313 Exit(EXIT_FAILURE); 00314 } 00315 path_to_atn_image = argv[i+1]; 00316 estimate_acf_flag = true; 00317 i++; 00318 } 00319 00320 // Output datafile in histogram mode 00321 else if (option=="-k") 00322 { 00323 kind_flag = true; 00324 } 00325 00326 // Provide an isotope alias 00327 else if (option=="-ist") 00328 { 00329 if (argv[i+1] == NULL) 00330 { 00331 Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl); 00332 Exit(EXIT_FAILURE); 00333 } 00334 else 00335 istp_alias = argv[i+1]; 00336 00337 i++; 00338 } 00339 00340 // Generate image from the sourceID 00341 else if (option=="-isrc") 00342 { 00343 if (argv[i+1] == NULL) 00344 { 00345 Cerr("***** castor-GATERootToCastor :: argument missing for option: " << option << endl); 00346 Exit(EXIT_FAILURE); 00347 } 00348 else 00349 { 00350 // Recover the path to the output image 00351 string input = argv[i+1]; 00352 int pos_colon = input.find_first_of(":"); 00353 path_to_src_image = input.substr(0,pos_colon); 00354 input = input.substr(pos_colon + 1); 00355 00356 // Get string section related to dimensions 00357 string p_dims_str[6]; 00358 if (ReadStringOption(input.c_str(), p_dims_str, 6, ",", option)) 00359 { 00360 Cerr("***** castor-GATERootToCastor :: Invalid argument " << argv[i+1] << " for option " << option << " !" << endl); 00361 Exit(EXIT_FAILURE); 00362 } 00363 00364 // Recover dimensions 00365 for(int i=0 ; i<3 ; i++) 00366 if (ConvertFromString(p_dims_str[i], &dim_src_img[i])) 00367 { 00368 Cerr("***** castor-GATERootToCastor :: Conversion error for elt " << p_dims_str[i] << " for option " << option << " !" << endl); 00369 Exit(EXIT_FAILURE); 00370 } 00371 00372 // Recover dimensions 00373 for(int i=0 ; i<3 ; i++) 00374 if (ConvertFromString(p_dims_str[i+3], &vox_src_img[i])) 00375 { 00376 Cerr("***** castor-GATERootToCastor :: Conversion error for elt " << p_dims_str[i+3] << " for option " << option << " !" << endl); 00377 Exit(EXIT_FAILURE); 00378 } 00379 00380 // Initilize source image 00381 p_src_img = new FLTNB[dim_src_img[0]*dim_src_img[1]*dim_src_img[2]]; 00382 00383 for(int v=0 ; v<dim_src_img[0]*dim_src_img[1]*dim_src_img[2] ; v++) 00384 p_src_img[v] = 0; 00385 00386 src_img_flag = true; 00387 } 00388 i++; 00389 } 00390 00391 // Generate geom file 00392 else if (option=="-geo") 00393 { 00394 geom_flag = true; 00395 } 00396 00397 // SPECT bins 00398 else if (option=="-spect_b") 00399 { 00400 string input = argv[i+1]; 00401 uint16_t s_bins[2]; 00402 if (ReadStringOption(input.c_str(), s_bins, 2, ",", option)) 00403 { 00404 Cerr("***** castor-GATERootToCastor :: Invalid argument " << argv[i+1] << " for option " << option << " !" << endl); 00405 Exit(EXIT_FAILURE); 00406 } 00407 00408 spect_bin_trs = s_bins[0]; 00409 spect_bin_axl = s_bins[1]; 00410 00411 i++; 00412 } 00413 00414 00415 // Verbosity level 00416 else if (option=="-vb") 00417 { 00418 if (i>=argc-1) 00419 { 00420 Cerr("***** castor-GATERootToCastor :: Argument missing for option: " << option << endl); 00421 Exit(EXIT_FAILURE); 00422 } 00423 vb = atoi(argv[i+1]); 00424 i++; 00425 } 00426 00427 else 00428 { 00429 Cerr("***** castor-GATERootToCastor :: Unknown option '" << option << "' !" << endl); 00430 Exit(EXIT_FAILURE); 00431 } 00432 } 00433 00434 00435 // ============================================================================================================ 00436 // Mandatory checks: 00437 // ============================================================================================================ 00438 00439 // Basic initialization checks (minimal initializations mandatory for the next steps) 00440 00441 // data files 00442 if (path_to_input_file.empty() ) 00443 { 00444 Cerr("***** castor-GATERootToCastor :: Please provide at least one data filename (-i or -if)" << endl); 00445 ShowHelp(0); 00446 Exit(EXIT_FAILURE); 00447 } 00448 else 00449 { 00450 // Check if we have interfile 00451 if(path_to_input_file.size() == 1) 00452 { 00453 string check; 00454 int rvalue = 0; 00455 rvalue = IntfKeyGetValueFromFile(path_to_input_file[0], "interfile", &check, 1, KEYWORD_OPTIONAL); 00456 00457 if(rvalue == 1) 00458 { 00459 // error 00460 Cerr("***** castor-GATERootToCastor :: Errot when trying to read file: " << path_to_input_file[0] << "!" << endl); 00461 Exit(EXIT_FAILURE); 00462 } 00463 else if(rvalue == 0) 00464 // key found, we have an interfile as input 00465 input_is_intf = true; 00466 } 00467 00468 if(vb >= 2) 00469 { 00470 Cout(" Selected root data-file(s) to convert: " << endl); 00471 for (size_t i=0 ; i<path_to_input_file.size() ; i++) 00472 Cout(path_to_input_file[i] << endl); 00473 } 00474 } 00475 00476 // Check ROOT is enabled if the input file is not interfile (SPECT) 00477 #ifndef CASTOR_ROOT 00478 if(!input_is_intf) 00479 { 00480 Cerr("***** castor-GATERootToCastor :: CASToR must be compiled with ROOT to read input root files (check installation instructions)" << endl); 00481 Exit(EXIT_FAILURE); 00482 } 00483 #endif 00484 00485 00486 // output directory 00487 if (path_to_out_file.empty() ) 00488 { 00489 Cerr("***** castor-GATERootToCastor :: Please provide the output file name (-o)" << endl); 00490 ShowHelp(0); 00491 Exit(EXIT_FAILURE); 00492 } 00493 else 00494 if(vb >= 2) Cout(" selected output file:" << path_to_out_file << endl); 00495 00496 // macro 00497 if (path_to_mac_file.empty()) 00498 { 00499 Cerr("***** castor-GATERootToCastor :: Please provide the macro file associated to the GATE root datafile (-m) :" << endl); 00500 ShowHelp(0); 00501 Exit(EXIT_FAILURE); 00502 } 00503 else 00504 if(vb >= 2) Cout(" selected macro file: " << path_to_mac_file << endl); 00505 00506 // scanner 00507 if (scanner_name.empty()) 00508 { 00509 Cerr("***** castor-GATERootToCastor :: Please provide a scanner alias (-s) :" << endl); 00510 ShowHelp(0); 00511 Exit(EXIT_FAILURE); 00512 } 00513 else 00514 if(vb >= 2) Cout(" selected scanner alias: " << scanner_name << endl); 00515 00516 00517 // (Required for options using interfile I/O) 00518 GetUserEndianness(); 00519 00520 00521 // ============================================================================================================ 00522 // SOutputManager object initialisation: 00523 // ============================================================================================================ 00524 00525 sOutputManager* p_outputManager = sOutputManager::GetInstance(); 00526 // Set verbose level 00527 p_outputManager->SetVerbose(vb); 00528 // Set MPI rank 00529 p_outputManager->SetMPIRank(0); 00530 00531 // Set path to the config directory 00532 if (p_outputManager->CheckConfigDir("")) 00533 { 00534 Cerr("***** castor-GATERootToCastor :: A problem occured while checking for the config directory path !" << endl); 00535 Exit(EXIT_FAILURE); 00536 } 00537 // Initialize output directory and base name 00538 if (p_outputManager->InitOutputDirectory(path_to_out_file, "")) 00539 { 00540 Cerr("*****castor-GATERootToCastor :: A problem occured while initializing output directory !" << endl); 00541 Exit(EXIT_FAILURE); 00542 } 00543 // Log command line 00544 if (p_outputManager->LogCommandLine(argc,argv)) 00545 { 00546 Cerr("***** castor-GATERootToCastor :: A problem occured while logging command line arguments !" << endl); 00547 Exit(EXIT_FAILURE); 00548 } 00549 00550 // Output progression options 00551 cout << std::fixed; 00552 cout << std::setprecision(2); 00553 00554 00555 // ============================================================================================================ 00556 // Check system type from the macro file 00557 // ============================================================================================================ 00558 GATE_system_type = GetGATESystemType(path_to_mac_file); 00559 00560 if(GATE_system_type<0) 00561 { 00562 // Unknown system 00563 Cerr("***** castor-GATERootToCastor :: GATE system type not found : " << endl); 00564 Cerr(" This script only supports conversion for cylindricalPET and ecat systems" << endl); 00565 Cerr(" The system type is recovered from the lines '/gate/systems/...'" << endl); 00566 Exit(EXIT_FAILURE); 00567 } 00568 00569 00570 // ============================================================================================================ 00571 // Generate the geom file from the mac file(s) is required 00572 // ============================================================================================================ 00573 if(geom_flag) 00574 { 00575 string scanner_repository = sOutputManager::GetInstance()->GetPathToConfigDir() + "scanner" + OS_SEP; 00576 string path_to_geom = scanner_repository + scanner_name + ".geom"; 00577 00578 // Check system type 00579 switch ( GATE_system_type ) 00580 { 00581 case GATE_SYS_CYLINDRICAL: 00582 if( vb>=2 )Cout(endl << " --- CylindricalPET system detected. Proceeding to conversion... --- " << endl << endl); 00583 00584 if(CreateGeomWithCylindrical(path_to_mac_file , path_to_geom) ) 00585 { 00586 Cerr("***** castor-GATERootToCastor :: An error occured while trying to process mac file for cylindrical system: " << path_to_mac_file << endl); 00587 Exit(EXIT_FAILURE); 00588 } 00589 break; 00590 00591 case GATE_SYS_ECAT: 00592 if( vb>=2 )Cout(endl << " --- ECAT system detected. Proceeding to conversion... --- " << endl << endl); 00593 00594 if(CreateGeomWithECAT(path_to_mac_file , path_to_geom) ) 00595 { 00596 Cerr("***** castor-GATEMacToGeom :: An error occured while trying to process mac file for ecat system: " << path_to_mac_file << endl); 00597 Exit(EXIT_FAILURE); 00598 } 00599 break; 00600 /* // Not yet validated 00601 case GATE_SYS_SPECT: 00602 if( vb>=2 )Cout(endl << " --- SPECThead system detected. Proceeding to conversion... --- " << endl << endl); 00603 00604 if(CreateGeomWithSPECT(path_to_mac_file , path_to_geom) ) 00605 { 00606 Cerr("***** castor-GATEMacToGeom :: An error occured while trying to process mac file for SPECT system: " << path_to_mac_file << endl); 00607 Exit(EXIT_FAILURE); 00608 } 00609 break;*/ 00610 00611 default: // Unknown system 00612 Cerr("***** castor-GATERootToCastor :: System type not found : " << endl); 00613 Cerr(" This script only supports conversion for cylindricalPET and ecat systems" << endl); 00614 //Cerr(" This script only supports conversion for cylindricalPET ecat and SPECThead systems" << endl); 00615 Cerr(" The system type is recovered from the lines '/gate/systems/...'" << endl); 00616 Exit(EXIT_FAILURE); 00617 break; 00618 } 00619 00620 if( vb>=2 )Cout(endl << " --- Conversion completed --- " << endl << endl); 00621 } 00622 00623 00624 // ============================================================================================================ 00625 // ScannerManager object initialisation: 00626 // ============================================================================================================ 00627 00628 sScannerManager* p_scannermanager = sScannerManager::GetInstance(); 00629 p_scannermanager->SetVerbose(vb); 00630 00631 // Check if the scanner exists and get the name from ScannerManager 00632 scanner_name = (scanner_name.find(OS_SEP)) ? 00633 scanner_name.substr(scanner_name.find_last_of(OS_SEP)+1) : 00634 scanner_name; 00635 00636 if(p_scannermanager->FindScannerSystem(scanner_name) ) 00637 { 00638 Cerr("**** castor-GATERootToCastor :: A problem occurred while searching for scanner system !" << endl); 00639 Exit(EXIT_FAILURE); 00640 } 00641 00642 // Build output file names 00643 path_to_data_filename = path_to_out_file + ".Cdf"; 00644 path_to_header_filename = path_to_out_file + ".Cdh"; 00645 00646 00647 // ============================================================================================================ 00648 // Instanciate/Initialize CASToR Datafile 00649 // ============================================================================================================ 00650 00651 // Instantiate & Initialize oImageDimensionsAndQuantification object, required for datafile generation (number of threads) 00652 oImageDimensionsAndQuantification* p_ID = new oImageDimensionsAndQuantification(); 00653 00654 // Set isotope if provided 00655 if(!istp_alias.empty()) 00656 { 00657 if(p_ID->SetPETIsotope(0, istp_alias) ) 00658 { 00659 Cerr("**** castor-GATERootToCastor :: A problem occurred while checking isotope name !" << endl); 00660 Exit(EXIT_FAILURE); 00661 } 00662 } 00663 00664 // Instanciate & Initialize iDataFilePET and Event objects 00665 vDataFile* Out_data_file = NULL; 00666 vEvent* Event = NULL; 00667 00668 uint16_t max_nb_lines_per_event = 1; // No compression for root files 00669 00670 if(GATE_system_type == GATE_SYS_SPECT) 00671 { 00672 Out_data_file = new iDataFileSPECT(); 00673 iDataFileSPECT* p_datafile = (dynamic_cast<iDataFileSPECT*>(Out_data_file)); 00674 p_datafile->SetDataType(TYPE_SPECT); 00675 p_datafile->SetIsotope(istp_alias); 00676 histo_out_flag = true; // force histogram output for SPECT 00677 } 00678 else 00679 { 00680 Out_data_file = new iDataFilePET(); 00681 iDataFilePET* p_datafile = (dynamic_cast<iDataFilePET*>(Out_data_file)); 00682 p_datafile->SetDataType(TYPE_PET); 00683 p_datafile->SetIsotope(istp_alias); 00684 00685 // ACF computed 00686 if(estimate_acf_flag) 00687 p_datafile->SetAtnCorrectionFlagOn(); 00688 00689 p_datafile->SetMaxNberLinesPerEvent(max_nb_lines_per_event); 00690 } 00691 00692 Out_data_file->SetImageDimensionsAndQuantification(p_ID); 00693 Out_data_file->SetHeaderDataFileName(path_to_header_filename); 00694 Out_data_file->SetPercentageLoad(0); // 0 (default) 00695 Out_data_file->SetVerbose(0); 00696 00697 00698 // Init Histogram-mode Event 00699 if(histo_out_flag) 00700 { 00701 Out_data_file->SetDataMode(MODE_HISTOGRAM); 00702 00703 // Instanciate histogram Event depending on modality 00704 if(GATE_system_type == GATE_SYS_SPECT) 00705 Event = new iEventHistoSPECT(); 00706 else 00707 Event = new iEventHistoPET(); 00708 00709 Event->SetNbLines(max_nb_lines_per_event); 00710 Event->AllocateID(); 00711 } 00712 00713 // Init List-mode Event 00714 else 00715 { 00716 Out_data_file->SetDataMode(MODE_LIST); 00717 00718 // Instanciate histogram Event depending on modality 00719 if(GATE_system_type == GATE_SYS_SPECT) 00720 { 00721 Event = new iEventListSPECT(); 00722 // record coincidence kind or not 00723 if(kind_flag) 00724 ((iDataFileSPECT*)Out_data_file)->SetEventKindFlagOn(); 00725 } 00726 else 00727 { 00728 Event = new iEventListPET(); 00729 // record coincidence kind or not 00730 if(kind_flag) 00731 ((iDataFilePET*)Out_data_file)->SetEventKindFlagOn(); 00732 } 00733 00734 Event->SetNbLines(max_nb_lines_per_event); 00735 Event->AllocateID(); 00736 } 00737 00738 Out_data_file->PROJ_InitFile(); 00739 Out_data_file->ComputeSizeEvent(); 00740 Out_data_file->PrepareDataFile(); 00741 00742 00743 // ============================================================================================================ 00744 // If acf estimation is enabled for histogram output, initialize all required objects for analytic projection 00745 // (Image space, projector and scanner geometry) 00746 // ============================================================================================================ 00747 oProjectorManager* p_ProjectorManager = new oProjectorManager(); 00748 oImageSpace* p_ImageSpace = new oImageSpace(); 00749 00750 00751 if(estimate_acf_flag) 00752 { 00753 if(!histo_out_flag) 00754 { 00755 Cerr("***** castor-GATERootToCastor :: Estimation of acf from an attenuation image (-atn option) is only available for histogram datafile output !" << endl 00756 << " Add the (-oh) option to the command line to enable this option." << endl); 00757 Exit(1); 00758 } 00759 00760 Intf_fields IF; 00761 IntfKeyInitFields(&IF); 00762 if(IntfReadHeader(path_to_atn_image, &IF, vb) ) 00763 { 00764 Cerr("***** castor-GATERootToCastor :: An error occurred while trying to read the interfile header of attenuation file " << path_to_atn_image << " !" << endl); 00765 Exit(1); 00766 } 00767 00768 // --- oImageDimensionsAndQuantification initialization --- 00769 p_ID->SetNbVoxX(IF.mtx_size[0]); 00770 p_ID->SetNbVoxY(IF.mtx_size[1]); 00771 p_ID->SetNbVoxZ(IF.mtx_size[2]); 00772 p_ID->SetNbThreads("1"); 00773 p_ID->SetNbBeds(1); 00774 p_ID->SetVoxSizeX(IF.vox_size[0]); 00775 p_ID->SetVoxSizeY(IF.vox_size[1]); 00776 p_ID->SetVoxSizeZ(IF.vox_size[2]); 00777 p_ID->SetFOVOutMasking(0., 0); 00778 p_ID->SetFOVSizeX(-1.); 00779 p_ID->SetFOVSizeY(-1.); 00780 p_ID->SetFOVSizeZ(-1.); 00781 p_ID->SetOffsetX(0); 00782 p_ID->SetOffsetY(0); 00783 p_ID->SetOffsetZ(0); 00784 p_ID->SetVerbose(vb); 00785 p_ID->SetNbRespGates(1); 00786 p_ID->SetNbCardGates(1); 00787 p_ID->SetFrames(""); 00788 00789 if (p_ID->CheckParameters()) 00790 { 00791 Cerr("***** castor-GATERootToCastor :: A problem occured while checking image dimensions parameters !" << endl); 00792 Exit(1); 00793 } 00794 if (p_ID->Initialize()) 00795 { 00796 Cerr("***** castor-GATERootToCastor :: A problem occured while initializing image dimensions !" << endl); 00797 Exit(1); 00798 } 00799 00800 // Initialization of DynamicDataManager class, related 4D data splitting management 00801 if (p_ID->InitDynamicData( "", 0, 0, 0, 1, 1 ) ) 00802 { 00803 Cerr("***** castor-GATERootToCastor :: A problem occured while initializing Dynamic data manager's class !" << endl); 00804 Exit(EXIT_FAILURE); 00805 } 00806 00807 00808 // --- Image space initialization --- 00809 p_ImageSpace->SetImageDimensionsAndQuantification(p_ID); 00810 p_ImageSpace->SetVerbose(vb); 00811 00812 // Allocate memory for main image 00813 p_ImageSpace->LMS_InstantiateImage(); 00814 00815 // Read attenuation image 00816 if(p_ImageSpace->PROJ_LoadInitialImage(path_to_atn_image) ) 00817 { 00818 Cerr("***** castor-GATERootToCastor :: Error during image initialization !" << endl); 00819 Exit(EXIT_FAILURE); 00820 } 00821 00822 00823 // --- Build Scanner geometry --- 00824 if(p_scannermanager->BuildScannerObject() ) 00825 { 00826 Cerr("**** castor-GATERootToCastor :: A problem occurred during scanner object construction !" << endl); 00827 Exit(EXIT_FAILURE); 00828 } 00829 00830 if(p_scannermanager->InstantiateScanner() ) 00831 { 00832 Cerr("**** castor-GATERootToCastor :: A problem occurred while creating Scanner object !" << endl); 00833 Exit(EXIT_FAILURE); 00834 } 00835 00836 if(p_scannermanager->BuildLUT() ) 00837 { 00838 Cerr("***** castor-GATERootToCastor :: A problem occurred while generating/reading the LUT !" << endl); 00839 Exit(EXIT_FAILURE); 00840 } 00841 00842 // Check the scanner manager parameters and initialize the scanner 00843 if (p_scannermanager->CheckParameters()) 00844 { 00845 Cerr("***** castor-GATERootToCastor :: A problem occured while checking scanner manager parameters !" << endl); 00846 Exit(1); 00847 } 00848 00849 if (p_scannermanager->Initialize()) 00850 { 00851 Cerr("***** castor-GATERootToCastor :: A problem occured while initializing scanner !" << endl); 00852 Exit(1); 00853 } 00854 00855 00856 // --- Projector Manager initialization--- 00857 p_ProjectorManager->SetScanner(p_scannermanager->GetScannerObject()); 00858 p_ProjectorManager->SetImageDimensionsAndQuantification(p_ID); 00859 p_ProjectorManager->SetDataFile(Out_data_file); 00860 p_ProjectorManager->SetComputationStrategy(FIXED_LIST_COMPUTATION_STRATEGY); 00861 p_ProjectorManager->SetOptionsForward("incrementalSiddon"); 00862 p_ProjectorManager->SetOptionsBackward("incrementalSiddon"); 00863 p_ProjectorManager->SetOptionsCommon(""); 00864 p_ProjectorManager->SetVerbose(vb); 00865 00866 // Check parameters 00867 if (p_ProjectorManager->CheckParameters()) 00868 { 00869 Cerr("***** castor-GATERootToCastor :: A problem occured while checking projector manager's parameters !" << endl); 00870 Exit(EXIT_FAILURE); 00871 } 00872 00873 // Initialize projector manager 00874 if (p_ProjectorManager->Initialize()) 00875 { 00876 Cerr("***** castor-GATERootToCastor :: A problem occured while initializing projector manager !" << endl); 00877 Exit(EXIT_FAILURE); 00878 } 00879 } 00880 00881 00882 // ============================================================================================================ 00883 // Parse Time offsets 00884 // ============================================================================================================ 00885 offset_time_list = new uint32_t[path_to_input_file.size()]; 00886 for(size_t f=0 ; f<path_to_input_file.size() ; f++) 00887 offset_time_list[f] = 0; 00888 00889 // Parse offset_time_str, if it contains any data 00890 if(offset_time_str != "") 00891 { 00892 vector<string> offsets; 00893 size_t comma_pos = 0; 00894 while ((comma_pos=offset_time_str.find_first_of(",")) != string::npos) 00895 { 00896 string offset = offset_time_str.substr(0,comma_pos); 00897 offsets.push_back(offset); 00898 offset_time_str = offset_time_str.substr(comma_pos+1); 00899 } 00900 00901 // Check we have a correct number of offsets 00902 if(offsets.size() != path_to_input_file.size()) 00903 { 00904 Cerr("**** castor-GATERootToCastor :: Unmatching number of offsets with -ot option ! " 00905 << offsets.size() << " have been found, while "<< path_to_input_file.size() <<"input file have been provided !" << endl); 00906 Exit(EXIT_FAILURE); 00907 } 00908 00909 for(size_t o=0 ; o<offsets.size() ; o++) 00910 { 00911 double offset; 00912 if(ConvertFromString(offsets[o] , &offset) ) 00913 { 00914 Cerr("**** castor-GATERootToCastor :: Error while trying to convert offset : "<< offsets[o] <<" in ms ! " << endl); 00915 Exit(EXIT_FAILURE); 00916 } 00917 00918 // convert in ms 00919 offset_time_list[o] = (uint32_t)offset*1000; 00920 } 00921 } 00922 00923 00924 // ============================================================================================================ 00925 // *********************************************** CONVERSION ************************************************* 00926 // ============================================================================================================ 00927 Cout(" --- Start conversion of datafile(s) : " << input_file << endl 00928 << " using mac file: " << path_to_mac_file << endl 00929 << " CASToR output header datafile: " << path_to_header_filename << endl 00930 << " CASToR output binary datafile: " << path_to_data_filename << endl << endl); 00931 00932 // ============================================================================================================ 00933 // Variables initialization 00934 // ============================================================================================================ 00935 00936 // Counter for the number of events (and the number of trues, scatters and randoms for simulated data) 00937 uint64_t nLORs_tot = 0, 00938 nLORs_trues = 0, 00939 nLORs_rdms = 0, 00940 nLORs_scatters =0, 00941 nLORs_unknown =0, 00942 nBINs = 0; 00943 00944 // Scanner variables (PET) 00945 uint32_t nCrystalsTot = 0; 00946 uint32_t nRsectorsPerRing = 1; 00947 uint32_t nModulesTransaxial = 1, nModulesAxial = 1; 00948 uint32_t nSubmodulesTransaxial = 1, nSubmodulesAxial = 1; 00949 uint32_t nCrystalsTransaxial = 1, nCrystalsAxial = 1; 00950 uint32_t nBlocksLine = 1, nBlocksPerRing = 1; 00951 uint8_t nLayers = 1; 00952 uint32_t nb_crystal_per_layer[4] = {0,0,0,0}; 00953 00954 // Scanner variables (SPECT) 00955 uint32_t nProjectionsByHead =1; 00956 uint32_t nProjectionsTot = 1; 00957 uint32_t nHeads =1; 00958 uint32_t nPixelsAxial = 1; 00959 uint32_t nPixelsTransaxial = 1; 00960 uint32_t nbSimulatedPixels = 1; 00961 float_t distToDetector = 0.; 00962 int headRotDirection = GEO_ROT_CW; 00963 float_t head1stAngleDeg = 0; 00964 float_t headAngPitchDeg = -1; 00965 float_t headAngStepDeg = -1; 00966 float_t crystalSizeAxl=-1., 00967 crystalSizeTrs=-1.; 00968 FLTNB* p_proj_spect_image=NULL; 00969 00970 // layer repeaters 00971 vector<uint32_t> nLayersRptTransaxial, nLayersRptAxial; 00972 00973 // Castor variables 00974 uint8_t** p_bins = NULL; 00975 uint32_t bins_elts = 0; 00976 uint32_t start_time_ms = 0; // 0 as default initialization 00977 uint32_t duration_ms = 1000; // 1 second as default initialization 00978 00979 00980 #ifdef CASTOR_ROOT 00981 uint32_t castorID1=0; 00982 uint32_t castorID2=0; 00983 uint8_t kind; 00984 00985 // ROOT data variables 00986 int32_t crystalID1, crystalID2; 00987 00988 // cylindricalPET specific 00989 int32_t rsectorID1, rsectorID2; 00990 int32_t moduleID1, moduleID2; 00991 int32_t submoduleID1, submoduleID2; 00992 int32_t layerID1, layerID2; 00993 00994 // ecat specific 00995 int32_t blockID1, blockID2; 00996 00997 // SPECThead specific 00998 float_t rotAngle; 00999 float_t gPosX, gPosY, gPosZ; 01000 int32_t headID; 01001 int32_t pixelID; 01002 01003 01004 // others 01005 int32_t eventID1= 0, eventID2= 0; 01006 int32_t comptonPhantomID1 = 0, comptonPhantomID2= 0; 01007 int32_t rayleighPhantomID1 = 0,rayleighPhantomID2 = 0; 01008 double_t time1= 0, time2= 0; 01009 int32_t sourceID1, sourceID2; 01010 float_t sourcePosX1, sourcePosY1, sourcePosZ1, 01011 sourcePosX2, sourcePosY2, sourcePosZ2; 01012 01013 01014 // ============================================================================================================ 01015 // ROOT objects declaration 01016 // ============================================================================================================ 01017 01018 TApplication *Tapp = new TApplication("tapp",0,0); 01019 TTree** GEvents = new TTree *[path_to_input_file.size()]; 01020 TFile** Tfile_root = new TFile*[path_to_input_file.size()]; 01021 01022 if(!input_is_intf) 01023 { 01024 // Compute the total number of LORs in the dataset 01025 for (size_t iFic=0 ; iFic<path_to_input_file.size() ; iFic++) 01026 { 01027 Tfile_root[iFic] = new TFile(path_to_input_file[iFic].c_str(),"READ","ROOT file with histograms"); 01028 if(GATE_system_type == GATE_SYS_SPECT) 01029 GEvents[iFic] = (TTree*)Tfile_root[iFic]->Get("Singles"); 01030 else 01031 GEvents[iFic] = (TTree*)Tfile_root[iFic]->Get("Coincidences"); 01032 01033 nLORs_tot += GEvents[iFic]->GetEntries(); 01034 } 01035 } 01036 #endif 01037 01038 01039 // Indexes for progression output 01040 uint64_t printing_index = 0; 01041 uint64_t printing_ratio = (nLORs_tot>10000) ? 10000 : nLORs_tot/10; 01042 01043 // ============================================================================================================ 01044 // Recover system geometric elements values 01045 // ============================================================================================================ 01046 01047 // SPECThead system 01048 if(GATE_system_type == GATE_SYS_SPECT) 01049 { 01050 // Data provided as an interfile projection image 01051 if(input_is_intf) 01052 { 01053 if(ReadIntfSPECT( path_to_input_file[0], 01054 distToDetector, 01055 nHeads, 01056 nPixelsAxial, 01057 nPixelsTransaxial, 01058 crystalSizeAxl, 01059 crystalSizeTrs, 01060 nProjectionsTot, 01061 nProjectionsByHead, 01062 head1stAngleDeg, 01063 headAngPitchDeg, 01064 headAngStepDeg, 01065 headRotDirection, 01066 start_time_ms, 01067 duration_ms, 01068 vb) ) 01069 { 01070 Cerr("**** castor-GATERootToCastor :: Error when reading mac file : "<< path_to_mac_file <<" !" << endl); 01071 Exit(EXIT_FAILURE); 01072 } 01073 01074 nCrystalsTot = nPixelsAxial * nPixelsTransaxial; 01075 01076 // Read Interfile 01077 Intf_fields IF_proj_spect_data; 01078 p_proj_spect_image = new FLTNB[nHeads*nProjectionsByHead*nCrystalsTot]; 01079 01080 if( IntfReadProjectionImage(path_to_input_file[0], p_proj_spect_image, &IF_proj_spect_data, vb, false) ) 01081 { 01082 Cerr("**** castor-GATERootToCastor :: Error when trying to read image : "<< path_to_input_file[0] <<" !" << endl); 01083 Exit(EXIT_FAILURE); 01084 } 01085 } 01086 01087 // Data provided as a root file 01088 else 01089 { 01090 if(ReadMacSPECT(path_to_mac_file, 01091 distToDetector, 01092 nHeads, 01093 nPixelsAxial, 01094 nPixelsTransaxial, 01095 crystalSizeAxl, 01096 crystalSizeTrs, 01097 nProjectionsTot, 01098 nProjectionsByHead, 01099 head1stAngleDeg, 01100 headAngPitchDeg, 01101 headAngStepDeg, 01102 headRotDirection, 01103 start_time_ms, 01104 duration_ms, 01105 vb) ) 01106 { 01107 Cerr("**** castor-GATERootToCastor :: Error when reading mac file : "<< path_to_mac_file <<" !" << endl); 01108 Exit(EXIT_FAILURE); 01109 } 01110 01111 nbSimulatedPixels = nPixelsAxial*nPixelsTransaxial; 01112 01113 if((spect_bin_trs>0 || spect_bin_axl>0) && nbSimulatedPixels > 1 ) 01114 { 01115 Cerr("**** castor-GATERootToCastor :: WARNING : Spect bins have been initialized, but the simulation already provide a specific number of pixels (="<< nPixelsAxial*nPixelsTransaxial <<") !"<< endl << 01116 " Pixel matrix used by default !" << endl); 01117 Exit(EXIT_FAILURE); 01118 } 01119 else // Check bins have been provided, and initialize pixel sizes with provided values 01120 { 01121 if(spect_bin_trs == 0 && spect_bin_axl == 0 && nbSimulatedPixels==1) 01122 { 01123 Cerr("**** castor-GATERootToCastor :: Error : Axial and transaxial bins values expected (use option -spect_b) !"<< endl); 01124 Exit(EXIT_FAILURE); 01125 } 01126 01127 // check crystal sizes have been correctly read 01128 if(crystalSizeAxl<0 || crystalSizeTrs<0) 01129 { 01130 Cerr("**** castor-GATERootToCastor :: Crystal dimensions not correctly read in the mac files !" << endl); 01131 Exit(EXIT_FAILURE); 01132 } 01133 01134 nPixelsTransaxial = spect_bin_trs; 01135 nPixelsAxial = spect_bin_axl; 01136 01137 if(vb>=2) Cout("Transaxial/Axial nb pixels : " << nPixelsTransaxial << " , " << nPixelsAxial << endl << 01138 "Transaxial/Axial pixel sizes : " << crystalSizeTrs/nPixelsTransaxial << " , " << crystalSizeAxl/nPixelsAxial << endl); 01139 } 01140 } 01141 01142 nCrystalsTot = nPixelsAxial * nPixelsTransaxial; 01143 01144 if(vb>=2) Cout("Number of Projections: " << nProjectionsByHead << endl << 01145 "Detected number of crystals in the system : " << nCrystalsTot << endl); 01146 01147 01148 // Histogram bin vector initialization using the total number of projections & pixels 01149 if(histo_out_flag) 01150 { 01151 bins_elts = nHeads*nProjectionsByHead; 01152 01153 p_bins = new uint8_t*[bins_elts]; 01154 for (size_t p=0; p<bins_elts ; p++) 01155 { 01156 p_bins[p] = new uint8_t[nCrystalsTot]; 01157 01158 for (size_t c=0; c<nCrystalsTot ; c++) 01159 { 01160 p_bins[p][c] = input_is_intf ? (uint8_t)p_proj_spect_image[p*nCrystalsTot+c] : 0 ; 01161 01162 // Get number of singles if input is interfile proj image 01163 if(input_is_intf) 01164 { 01165 nLORs_tot += p_proj_spect_image[p*nCrystalsTot+c]; 01166 nLORs_unknown += p_proj_spect_image[p*nCrystalsTot+c]; 01167 } 01168 } 01169 } 01170 } 01171 01172 delete[] p_proj_spect_image; 01173 } 01174 01175 else // PET systems 01176 { 01177 // cylindricalPET system 01178 if(GATE_system_type == GATE_SYS_CYLINDRICAL) 01179 { 01180 if(ReadMacCylindrical(path_to_mac_file, 01181 nLayers, 01182 nb_crystal_per_layer, 01183 nCrystalsTot, 01184 nCrystalsAxial, 01185 nCrystalsTransaxial, 01186 nLayersRptAxial, 01187 nLayersRptTransaxial, 01188 nSubmodulesAxial, 01189 nSubmodulesTransaxial, 01190 nModulesAxial, 01191 nModulesTransaxial, 01192 nRsectorsPerRing, 01193 start_time_ms, 01194 duration_ms, 01195 vb) ) 01196 { 01197 Cerr("**** castor-GATERootToCastor :: Error when reading mac file : "<< path_to_mac_file <<" !" << endl); 01198 Exit(EXIT_FAILURE); 01199 } 01200 } 01201 01202 else // ECAT 01203 { 01204 // Reading the macro file 01205 if(ReadMacECAT(path_to_mac_file, 01206 nCrystalsTot, 01207 nCrystalsAxial, 01208 nCrystalsTransaxial, 01209 nBlocksLine, 01210 nBlocksPerRing, 01211 start_time_ms, 01212 duration_ms, 01213 vb) ) 01214 { 01215 Cerr("**** castor-GATERootToCastor :: Error when reading mac file : "<< path_to_mac_file <<" !" << endl); 01216 Exit(EXIT_FAILURE); 01217 } 01218 01219 } 01220 01221 if(vb>=2) Cout("Detected number of crystals in the system : " << nCrystalsTot << endl); 01222 01223 // Histogram bin vector initialization using the total number of crystals 01224 if(histo_out_flag) 01225 { 01226 Cout(" Allocating memory for bins... " << endl << 01227 " Warning : this step could required huge amount of memory if the system contains a high number of crystals !" << endl); 01228 01229 bins_elts = nCrystalsTot; 01230 01231 p_bins = new uint8_t*[bins_elts]; 01232 for (size_t c=0; c<bins_elts ; c++) 01233 { 01234 p_bins[c] = new uint8_t[nCrystalsTot-c]; 01235 01236 for (size_t c2=0; c2<nCrystalsTot-c ; c2++) 01237 p_bins[c][c2] = 0; 01238 } 01239 01240 Cout(" Memory allocation for bins completed " << endl << endl); 01241 } 01242 01243 } // end of PET section 01244 01245 01246 01247 01248 // ============================================================================================================ 01249 // Loop on the input datafile(s) and process event by event 01250 // ============================================================================================================ 01251 if(!input_is_intf) // SPECT projection interfile : data already processed 01252 for (size_t iFic=0 ; iFic<path_to_input_file.size() ; iFic++) 01253 { 01254 if(vb>=2) Cout(endl << "Processing datafile " << iFic << " : " << path_to_input_file[iFic] << "..." << endl); 01255 01256 // Set variables of the root tree 01257 // If we got a cylindricalPET system 01258 #ifdef CASTOR_ROOT 01259 if(GATE_system_type == GATE_SYS_CYLINDRICAL) 01260 { 01261 GEvents[iFic]->SetBranchAddress("time1",&time1); 01262 GEvents[iFic]->SetBranchAddress("rsectorID1",&rsectorID1); 01263 GEvents[iFic]->SetBranchAddress("moduleID1",&moduleID1); 01264 GEvents[iFic]->SetBranchAddress("submoduleID1",&submoduleID1); 01265 GEvents[iFic]->SetBranchAddress("crystalID1",&crystalID1); 01266 GEvents[iFic]->SetBranchAddress("layerID1", &layerID1); 01267 GEvents[iFic]->SetBranchAddress("comptonPhantom1",&comptonPhantomID1); 01268 GEvents[iFic]->SetBranchAddress("RayleighPhantom1",&rayleighPhantomID1); 01269 GEvents[iFic]->SetBranchAddress("eventID1",&eventID1); 01270 GEvents[iFic]->SetBranchAddress("sourceID1",&sourceID1); 01271 01272 GEvents[iFic]->SetBranchAddress("time2",&time2); 01273 GEvents[iFic]->SetBranchAddress("rsectorID2",&rsectorID2); 01274 GEvents[iFic]->SetBranchAddress("moduleID2",&moduleID2); 01275 GEvents[iFic]->SetBranchAddress("submoduleID2",&submoduleID2); 01276 GEvents[iFic]->SetBranchAddress("crystalID2",&crystalID2); 01277 GEvents[iFic]->SetBranchAddress("layerID2", &layerID2); 01278 GEvents[iFic]->SetBranchAddress("comptonPhantom2",&comptonPhantomID2); 01279 GEvents[iFic]->SetBranchAddress("RayleighPhantom2",&rayleighPhantomID2); 01280 GEvents[iFic]->SetBranchAddress("eventID2",&eventID2); 01281 GEvents[iFic]->SetBranchAddress("sourceID2",&sourceID2); 01282 01283 GEvents[iFic]->SetBranchAddress("sourcePosX1",&sourcePosX1); 01284 GEvents[iFic]->SetBranchAddress("sourcePosY1",&sourcePosY1); 01285 GEvents[iFic]->SetBranchAddress("sourcePosZ1",&sourcePosZ1); 01286 } 01287 else if(GATE_system_type == GATE_SYS_SPECT) 01288 { 01289 GEvents[iFic]->SetBranchAddress("time",&time1); 01290 GEvents[iFic]->SetBranchAddress("headID", &headID); 01291 GEvents[iFic]->SetBranchAddress("crystalID",&crystalID1); 01292 GEvents[iFic]->SetBranchAddress("pixelID", &pixelID); 01293 GEvents[iFic]->SetBranchAddress("rotationAngle", &rotAngle); 01294 GEvents[iFic]->SetBranchAddress("globalPosX",&gPosX); 01295 GEvents[iFic]->SetBranchAddress("globalPosY",&gPosY); 01296 GEvents[iFic]->SetBranchAddress("globalPosZ",&gPosZ); 01297 GEvents[iFic]->SetBranchAddress("comptonPhantom",&comptonPhantomID1); 01298 GEvents[iFic]->SetBranchAddress("RayleighPhantom",&rayleighPhantomID1); 01299 GEvents[iFic]->SetBranchAddress("sourcePosX",&sourcePosX1); 01300 GEvents[iFic]->SetBranchAddress("sourcePosY",&sourcePosY1); 01301 GEvents[iFic]->SetBranchAddress("sourcePosZ",&sourcePosZ1); 01302 } 01303 else 01304 { 01305 GEvents[iFic]->SetBranchAddress("time1",&time1); 01306 GEvents[iFic]->SetBranchAddress("blockID1",&blockID1); 01307 GEvents[iFic]->SetBranchAddress("crystalID1",&crystalID1); 01308 GEvents[iFic]->SetBranchAddress("comptonPhantom1",&comptonPhantomID1); 01309 GEvents[iFic]->SetBranchAddress("RayleighPhantom1",&rayleighPhantomID1); 01310 GEvents[iFic]->SetBranchAddress("eventID1",&eventID1); 01311 GEvents[iFic]->SetBranchAddress("sourceID1",&sourceID1); 01312 01313 GEvents[iFic]->SetBranchAddress("time2",&time2); 01314 GEvents[iFic]->SetBranchAddress("blockID2",&blockID2); 01315 GEvents[iFic]->SetBranchAddress("crystalID2",&crystalID2); 01316 GEvents[iFic]->SetBranchAddress("comptonPhantom2",&comptonPhantomID2); 01317 GEvents[iFic]->SetBranchAddress("RayleighPhantom2",&rayleighPhantomID2); 01318 GEvents[iFic]->SetBranchAddress("eventID2",&eventID2); 01319 GEvents[iFic]->SetBranchAddress("sourceID2",&sourceID2); 01320 GEvents[iFic]->SetBranchAddress("sourcePosX2",&sourcePosX2); 01321 GEvents[iFic]->SetBranchAddress("sourcePosY2",&sourcePosY2); 01322 GEvents[iFic]->SetBranchAddress("sourcePosZ2",&sourcePosZ2); 01323 } 01324 01325 01326 01327 01328 //---------------------------------------------- /ROOT data variables 01329 01330 // Loop on the GEvents in the current datafile 01331 for (int i=0; i<GEvents[iFic]->GetEntries() ; i++) 01332 { 01333 GEvents[iFic]->GetEntry(i); 01334 01335 printing_index++; 01336 01337 // ID Conversions 01338 if (GATE_system_type == GATE_SYS_CYLINDRICAL) 01339 { 01340 castorID1 = ConvertIDcylindrical(nRsectorsPerRing, 01341 nModulesTransaxial, nModulesAxial, 01342 nSubmodulesTransaxial, nSubmodulesAxial, 01343 nCrystalsTransaxial, nCrystalsAxial, 01344 nLayers, nb_crystal_per_layer, 01345 nLayersRptTransaxial, nLayersRptAxial, 01346 layerID1, crystalID1, submoduleID1, moduleID1, rsectorID1); 01347 01348 castorID2 = ConvertIDcylindrical(nRsectorsPerRing, 01349 nModulesTransaxial, nModulesAxial, 01350 nSubmodulesTransaxial, nSubmodulesAxial, 01351 nCrystalsTransaxial, nCrystalsAxial, 01352 nLayers, nb_crystal_per_layer, 01353 nLayersRptTransaxial, nLayersRptAxial, 01354 layerID2, crystalID2, submoduleID2, moduleID2, rsectorID2); 01355 01356 if(vb >= 4) 01357 { 01358 Cout("File#" << iFic << ", event#" << i << endl;); 01359 Cout("Crystal 1 : RsectorID: " << rsectorID1 << ", moduleID: " << moduleID1 << ", submoduleID: " << submoduleID1 << ", crystalID: " << crystalID1 01360 << ", layerID: " << layerID1 << " --> castorID1:" << castorID1 << endl;); 01361 Cout("Crystal 2 :RsectorID: " << rsectorID2 << ", moduleID2: " << moduleID2 << ", submoduleID: " << submoduleID2 << ", crystalID: " << crystalID2 01362 << ", layerID: " << layerID2 << " --> castorID2:" << castorID2 << endl;); 01363 } 01364 } 01365 01366 else if (GATE_system_type == GATE_SYS_SPECT) 01367 { 01368 // Get projection ID 01369 castorID1 = ConvertIDSPECTRoot1(headID, 01370 rotAngle, 01371 headAngStepDeg, 01372 nProjectionsByHead); 01373 01374 // Get crystal ID 01375 castorID2 = ConvertIDSPECTRoot2(nbSimulatedPixels, 01376 nPixelsAxial, 01377 nPixelsTransaxial, 01378 headID, 01379 crystalID1, 01380 pixelID, 01381 rotAngle, 01382 headAngPitchDeg, 01383 crystalSizeAxl, 01384 crystalSizeTrs, 01385 gPosX, 01386 gPosY, 01387 gPosZ); 01388 01389 if(vb >= 4) 01390 { 01391 Cout("File#" << iFic << ", event#" << i << endl;); 01392 Cout("Projection ID : headID: " << headID << ", rotation angle (deg): " << rotAngle + headID*headAngPitchDeg << ", computed angular step between projections: " << headAngStepDeg << " --> castorID1:" << castorID1 << endl;); 01393 Cout("Crystal ID : crystalID: " << crystalID1 << ", pixelID: " << pixelID << ", rotAngle: " << rotAngle << 01394 ", globalPosX " << gPosX << ", globalPosY " << gPosY << ", globalPosZ " << gPosZ << " --> castorID2:" << castorID2 << endl;); 01395 } 01396 01397 } 01398 01399 else if (GATE_system_type == GATE_SYS_ECAT) 01400 { 01401 castorID1 = ConvertIDecat(nBlocksPerRing, nBlocksLine, nCrystalsTransaxial, nCrystalsAxial, crystalID1, blockID1); 01402 castorID2 = ConvertIDecat(nBlocksPerRing, nBlocksLine, nCrystalsTransaxial, nCrystalsAxial, crystalID2, blockID2); 01403 } 01404 01405 // Find out the kind of coincidence (true, scatter, random) 01406 kind = ComputeKindGATEEvent(eventID1, eventID2, comptonPhantomID1, comptonPhantomID2, rayleighPhantomID1, rayleighPhantomID2); 01407 01408 // Count nb LORs according to kind 01409 switch (kind) 01410 { 01411 case 1: 01412 nLORs_trues++; 01413 break; 01414 01415 case 2: 01416 nLORs_scatters++; 01417 break; 01418 01419 case 3: 01420 nLORs_scatters++; 01421 break; 01422 01423 case 4: 01424 nLORs_rdms++; 01425 break; 01426 01427 default: 01428 nLORs_unknown++; 01429 } 01430 01431 // Skip next step if event is not true if we only recover true 01432 if (true_only_flag==true && kind != KIND_TRUE) 01433 continue; 01434 01435 01436 // --- Write Event --- 01437 01438 // SPECT event 01439 if(GATE_system_type == GATE_SYS_SPECT) 01440 { 01441 // HISTOGRAM mode 01442 if(histo_out_flag) 01443 p_bins[castorID1][castorID2]++; 01444 else 01445 // LIST-mode 01446 { 01447 // Write event in the datafile 01448 uint32_t time_in_ms = time1*1000; 01449 Event->SetNbLines(1); 01450 Event->SetID1(0, castorID1); // 1st argument : line, 2nd argument :index 01451 Event->SetID2(0, castorID2); // 1st argument : line, 2nd argument :index 01452 Event->SetTimeInMs(time_in_ms); 01453 01454 ((iEventListSPECT*)Event)->SetKind(kind); 01455 01456 Out_data_file->PROJ_WriteEvent(Event, 0); 01457 } 01458 } 01459 else // PET event 01460 { 01461 // HISTOGRAM mode 01462 if(histo_out_flag) 01463 { 01464 if (castorID1 < castorID2) 01465 p_bins[castorID1][castorID2-castorID1-1]++; 01466 else 01467 p_bins[castorID2][castorID1-castorID2-1]++; 01468 } 01469 // LIST-mode 01470 else 01471 { 01472 // Write event in the datafile 01473 uint32_t time_in_ms = time1*1000; 01474 int nb_lines_in_event = 1; // 1 by default for GATE root files 01475 01476 Event->SetNbLines(nb_lines_in_event); 01477 Event->SetID1(0, castorID1); // 1st argument : line, 2nd argument :index 01478 Event->SetID2(0, castorID2); // 1st argument : line, 2nd argument :index 01479 Event->SetTimeInMs(time_in_ms); 01480 01481 ((iEventListPET*)Event)->SetKind(kind); 01482 01483 Out_data_file->PROJ_WriteEvent(Event, 0); 01484 } 01485 } 01486 01487 01488 // Source image (if no rdm event) 01489 if(src_img_flag && 01490 eventID1 == eventID2) 01491 { 01492 01493 INTNB x,y,z; 01494 x = (int)(( sourcePosX1 + dim_src_img[0]/2*vox_src_img[0]) / vox_src_img[0]); 01495 // CASToR Y axis inverted in comparison with GATE (left-handed coordinate) 01496 y = (int)(( dim_src_img[1]/2*vox_src_img[1] -sourcePosY1 ) / vox_src_img[1]); 01497 z = (int)(( sourcePosZ1 + dim_src_img[2]/2*vox_src_img[2]) / vox_src_img[2]); 01498 01499 if(x >= 0 && x < dim_src_img[0] && 01500 y >= 0 && y < dim_src_img[1] && 01501 z >= 0 && z < dim_src_img[2] ) 01502 p_src_img[z*dim_src_img[1]*dim_src_img[0] + y*dim_src_img[0] + x]++; 01503 } 01504 01505 01506 // Progression 01507 if (printing_index%(nLORs_tot/printing_ratio) == 0) 01508 { 01509 FLTNB percent = ( ((FLTNB)(printing_index+1))/((FLTNB)nLORs_tot) ) * ((FLTNB)100); 01510 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 " 01511 << percent << " \% "; 01512 } 01513 01514 01515 } // end of loop on events in input datafile iFic 01516 #endif 01517 01518 if(vb>=2) Cout(endl << "Datafile " << iFic << " : " << path_to_input_file[iFic] << " process OK" << endl); 01519 01520 } // end of loop on input datafiles 01521 01522 // Free Root objects 01523 #ifdef CASTOR_ROOT 01524 delete[] Tfile_root; 01525 delete[] GEvents; 01526 delete Tapp; 01527 #endif 01528 01529 01530 01531 // ============================================================================================================ 01532 // Write the HISTOGRAM datafile and header 01533 // ============================================================================================================ 01534 if(GATE_system_type == GATE_SYS_SPECT) 01535 { 01536 if (histo_out_flag) 01537 { 01538 printing_index=0; 01539 01540 // Writing the datafile 01541 Cout(endl << "Generate the histogram datafile" << endl); 01542 01543 uint32_t nb_bins = bins_elts*nCrystalsTot; 01544 printing_ratio = (nb_bins>1000) ? 1000 : nb_bins/10; 01545 01546 // Loop on the crystal IDs 01547 for (size_t id1=0 ; id1<bins_elts ; id1++) 01548 for (size_t id2=0; id2<nCrystalsTot ; id2++) 01549 { 01550 uint32_t nb_events_in_bin = p_bins[id1][id2]; 01551 01552 int nb_lines = 1; // 1 by default for GATE root file; 01553 Event->SetNbLines(nb_lines); 01554 Event->SetID1(0, id1); // 1st argument : line, 2nd argument :index 01555 Event->SetID2(0, id2); // 1st argument : line, 2nd argument :index 01556 Event->SetEventValue(0, nb_events_in_bin); 01557 01558 // Write event in Datafile 01559 Out_data_file->PROJ_WriteEvent(Event, 0); 01560 nBINs++; 01561 01562 // Progression 01563 01564 if (printing_index%((nb_bins)/printing_ratio) == 0) 01565 { 01566 FLTNB percent = ( ((FLTNB)(printing_index+1))/((FLTNB)nb_bins) ) * ((FLTNB)100); 01567 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 " 01568 << percent << " \% "; 01569 } 01570 01571 printing_index++; 01572 } 01573 01574 Out_data_file->SetStartTime((FLTNB)(start_time_ms)/1000); // start time in s 01575 Out_data_file->SetDuration((FLTNB)(duration_ms)/1000); // duration in s 01576 Out_data_file->SetNbEvents( nBINs ); 01577 01578 Cout(endl << "The output histogram contains " << nBINs << " events." << endl); 01579 } 01580 // Write the List-mode datafile header 01581 else 01582 { 01583 01584 Out_data_file->SetStartTime((FLTNB)(start_time_ms)/1000); // start time in s 01585 Out_data_file->SetDuration((FLTNB)(duration_ms)/1000); // duration in s 01586 //uint32_t nLORs = (true_only_flag) ? nLORs_trues : nLORs_tot; 01587 Out_data_file->SetNbEvents( (true_only_flag) ? 01588 nLORs_trues : 01589 nLORs_tot ) ; 01590 } 01591 01592 FLTNB* p_angles = new FLTNB[nProjectionsTot]; 01593 FLTNB* p_distToDetector = new FLTNB[nProjectionsTot]; 01594 01595 for(size_t h=0 ; h<nHeads ; h++) 01596 for(size_t p=0 ; p<nProjectionsByHead ; p++) 01597 { 01598 int idx_proj = h*nProjectionsByHead+p; 01599 p_angles[idx_proj] = head1stAngleDeg // initial angular position of the head 01600 + p*headAngStepDeg // projection angular position 01601 + h*headAngPitchDeg; // angular shift for each head 01602 01603 // same distance for each head with GATE 01604 p_distToDetector[idx_proj] = distToDetector; 01605 } 01606 01607 ((iDataFileSPECT*)Out_data_file)->SetNbBins(nPixelsTransaxial, nPixelsAxial); 01608 ((iDataFileSPECT*)Out_data_file)->SetNbProjections(nProjectionsTot); 01609 01610 if( ((iDataFileSPECT*)Out_data_file)->InitAngles(p_angles) ) 01611 { 01612 Cerr("**** castor-GATERootToCastor :: Error when trying to set projection angles values !" << endl); 01613 Exit(EXIT_FAILURE); 01614 } 01615 01616 if( ((iDataFileSPECT*)Out_data_file)->InitCorToDetectorDistance(p_distToDetector) ) 01617 { 01618 Cerr("**** castor-GATERootToCastor :: Error when trying to set distance between center of rotation and detectors !" << endl); 01619 Exit(EXIT_FAILURE); 01620 } 01621 01622 ((iDataFileSPECT*)Out_data_file)->SetHeadRotDirection(headRotDirection); 01623 01624 delete[] p_angles; 01625 delete[] p_distToDetector; 01626 01627 Out_data_file->PROJ_WriteHeader(); 01628 } // end of SPECT section 01629 01630 else 01631 { 01632 if (histo_out_flag) 01633 { 01634 printing_index=0; 01635 01636 // Writing the datafile 01637 Cout(endl << "Generate the histogram datafile" << endl); 01638 uint32_t nb_bins = (nCrystalsTot*nCrystalsTot - nCrystalsTot)/2; 01639 printing_ratio = (nb_bins>1000) ? 1000 : nb_bins/10; 01640 01641 // Loop on the crystal IDs 01642 for (size_t id1=0 ; id1 <bins_elts ; id1++) 01643 for (size_t id2 = id1+1; id2 < nCrystalsTot;id2++) 01644 { 01645 uint32_t nb_events_in_bin = p_bins[id1][id2-id1-1]; 01646 01647 int nb_lines = 1; // 1 by default for GATE root file; 01648 Event->SetNbLines(nb_lines); 01649 Event->SetID1(0, id1); // 1st argument : line, 2nd argument :index 01650 Event->SetID2(0, id2); // 1st argument : line, 2nd argument :index 01651 Event->SetEventValue(0, nb_events_in_bin); 01652 01653 // Estimate acf for the bin 01654 if(estimate_acf_flag) 01655 { 01656 // Compute the system matrix elements for the two crystals 01657 oProjectionLine* line = p_ProjectorManager->ComputeProjectionLine(Event, 0); 01658 01659 // Compute forward projection of attenuation image 01660 FLTNB fp = 0.; 01661 if (line->NotEmptyLine()) 01662 fp = line->ForwardProject(p_ImageSpace->m4p_image[0][0][0]) * 0.1; // 0.1 -> convert in mm-1 01663 01664 // Write atn correction factor in Event 01665 ((iEventPET*)Event)->SetAttenuationCorrectionFactor(1/std::exp(-fp)); 01666 } 01667 01668 // Write event in Datafile 01669 Out_data_file->PROJ_WriteEvent(Event, 0); 01670 nBINs++; 01671 01672 // Progression 01673 if (printing_index%((nb_bins)/printing_ratio) == 0) 01674 { 01675 FLTNB percent = ( ((FLTNB)(printing_index+1))/((FLTNB)nb_bins) ) * ((FLTNB)100); 01676 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 " 01677 << percent << " \% "; 01678 } 01679 printing_index++; 01680 } 01681 01682 Out_data_file->SetStartTime((FLTNB)(start_time_ms)/1000); // start time in s 01683 Out_data_file->SetDuration((FLTNB)(duration_ms)/1000); // duration in s 01684 Out_data_file->SetNbEvents( nBINs ); 01685 Out_data_file->PROJ_WriteHeader(); 01686 01687 Cout(endl << "The output histogram contains " << nBINs << " events." << endl); 01688 cout << "trest " << endl; 01689 } 01690 // Write the List-mode datafile header 01691 else 01692 { 01693 Out_data_file->SetStartTime((FLTNB)(start_time_ms)/1000); // start time in s 01694 Out_data_file->SetDuration((FLTNB)(duration_ms)/1000); // duration in s 01695 //uint32_t nLORs = (true_only_flag) ? nLORs_trues : nLORs_tot; 01696 Out_data_file->SetNbEvents( (true_only_flag) ? 01697 nLORs_trues : 01698 nLORs_tot ); 01699 Out_data_file->PROJ_WriteHeader(); 01700 } 01701 } // end of PET section 01702 01703 01704 // Write the source image 01705 if(src_img_flag) 01706 { 01707 if (vb>=2) Cout("Writing source image ..." << endl); 01708 01709 // Initialize Intf_fields object with the source image dimensions 01710 Intf_fields IF; 01711 IntfKeyInitFields(&IF); 01712 01713 IF.mtx_size[0] = dim_src_img[0]; 01714 IF.mtx_size[1] = dim_src_img[1]; 01715 IF.mtx_size[2] = dim_src_img[2]; 01716 IF.vox_size[0] = vox_src_img[0]; 01717 IF.vox_size[1] = vox_src_img[1]; 01718 IF.vox_size[2] = vox_src_img[2]; 01719 IF.nb_total_imgs = dim_src_img[2]; 01720 01721 if (IntfWriteImgFile(path_to_src_image.append("_src"), p_src_img, IF, vb) ) 01722 { 01723 Cerr("***** castor-GATERootToCastor :: Error writing Interfile of output image !" << endl); 01724 Exit(EXIT_FAILURE); 01725 } 01726 } 01727 01728 Cout(endl << "The simulated dataset contained " << nLORs_tot << " coincidences/singles, with: " << endl 01729 << " " << nLORs_trues << " trues (" << 100*nLORs_trues/nLORs_tot << " %), " << endl 01730 << " " << nLORs_scatters << " scatters (" << 100*nLORs_scatters/nLORs_tot << " %)," << endl 01731 << " " << nLORs_rdms << " randoms (" << 100*nLORs_rdms/nLORs_tot << " %)," << endl 01732 << " " << nLORs_unknown << " unknown (" << 100*nLORs_unknown/nLORs_tot << " %)." << endl); 01733 01734 if (vb>=2) Cout("Writing raw datafile ..." << endl); // todo just change name if only one datafile 01735 01736 01737 Out_data_file->PROJ_WriteData(); 01738 Out_data_file->PROJ_DeleteTmpDatafile(); 01739 01740 Cout(endl << " --- Conversion completed --- " << endl); 01741 01742 // ============================================================================================================ 01743 // End 01744 // ============================================================================================================ 01745 01746 if (vb>=2) Cout(" Deallocating objects ..." << endl); 01747 01748 // Delete objects 01749 if(histo_out_flag) 01750 { 01751 for(size_t b=0; b<bins_elts ; b++) 01752 delete[] p_bins[b]; 01753 01754 delete[]p_bins; 01755 01756 if(Event) 01757 { 01758 if(GATE_system_type == GATE_SYS_SPECT) 01759 delete (iEventHistoSPECT*)Event; 01760 else 01761 delete (iEventHistoPET*)Event; 01762 } 01763 } 01764 else 01765 if(Event) 01766 { 01767 if(GATE_system_type == GATE_SYS_SPECT) 01768 delete (iEventListSPECT*)Event; 01769 else 01770 delete (iEventListPET*)Event; 01771 } 01772 01773 delete[]offset_time_list; 01774 delete Out_data_file; 01775 01776 if(src_img_flag && p_src_img) 01777 delete[] p_src_img; 01778 01779 // Free objects created for analytic projection (acf estimation) 01780 if(estimate_acf_flag) 01781 p_ImageSpace->LMS_DeallocateImage(); 01782 if(p_ImageSpace) delete p_ImageSpace; 01783 if(p_ProjectorManager) delete p_ProjectorManager; 01784 if(p_ID) delete p_ID; 01785 01786 Cout(" --- END --- " << endl << endl); 01787 01788 return 0; 01789 }