![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00016 #include "gVariables.hh" 00017 #include "gOptions.hh" 00018 #include "iDataFilePET.hh" 00019 #include "sOutputManager.hh" 00020 #include "sScannerManager.hh" 00021 #include <iomanip> //std::setprecision 00022 00023 #ifdef CASTOR_ROOT 00024 #include "TROOT.h" 00025 #include "TApplication.h" 00026 #include "TGClient.h" 00027 #include "TCanvas.h" 00028 #include "TSystem.h" 00029 #include "TTree.h" 00030 #include "TBranch.h" 00031 #include "TFile.h" 00032 #endif 00033 00038 void ShowHelp(int a_returnCode) 00039 { 00040 cout << endl; 00041 cout << "NOTE : This program is an example code providing guidance regarding how to generate a PET CASToR datafile from any system datafile " << endl; 00042 cout << " It does NOT perform any conversion by itself and must be adjusted to the conversion of any system dataset !" << endl; 00043 cout << endl; 00044 cout << endl; 00045 cout << "Usage: castor-datafileConversionEx -ih input_datafile : (if datafile is histogram mode)" << endl; 00046 cout << " -il input_datafile : (if datafile is list-mode)" << endl; 00047 cout << " -o output_datafile" << endl; 00048 cout << " -s scanner_alias "<< endl; 00049 cout << " [-nc norm_factors_file" << endl; 00050 cout << " [-sc scat_factors_file" << endl; 00051 cout << " [-rc rdm_factors_file" << endl; 00052 cout << " [-ac atn_factors_file" << endl; 00053 cout << " [-cf calibration factor" << endl; 00054 cout << " [-ist isotope_alias" << endl; 00055 00056 //cout << " [-poi flag for POI]" << endl; 00057 //cout << " [-tof flag for TOF]" << endl; 00058 00059 cout << " [-vb verbosity]" << endl; 00060 cout << endl; 00061 00062 cout << endl; 00063 cout << "[Main settings]:" << endl; 00064 cout << " -ih path_to_histo_datafile : provide an input histogram datafile to convert" << endl; 00065 cout << " -il path_to_listm_datafile : provide an input list-mode datafile to convert" << endl; 00066 cout << " -o path_to_cstr_file.cdh : provide the path to the output file will be created inside this folder (no default)" << endl; 00067 cout << " -s scanner_alias : provide the name of the scanner used for to acquire the original data" << endl; 00068 cout << " : Must correspond to a .geom or .hscan file in the config/scanner repository." << endl; 00069 cout << endl; 00070 cout << "[Optional settings]:" << endl; 00071 cout << " -nc norm_factors_file : provide a file containing normalization correction factors" << endl; 00072 cout << " -sc scat_factors_file : provide a file containing scatter correction factors" << endl; 00073 cout << " -rc rdm_factors_file : provide a file containing random correction factors" << endl; 00074 cout << " -ac atn_factors_file : provide a file containing attenuation correction factors" << endl; 00075 cout << " -cf calibration factor : provide a calibration factor" << endl; 00076 cout << " -ist isotope_alias : provide alias of the isotope used in the input datafile"<< endl; 00077 cout << " Supported PET isotopes and their parameters are listed in config/misc/isotopes_pet"<< endl; 00078 cout << " Other isotopes can be added in the same file"<< endl; 00079 cout << endl; 00080 cout << "[Miscellaneous settings]:" << endl; 00081 cout << " -vb : give the verbosity level, from 0 (no verbose) to above 5 (at the event level) (default: 1)." << endl; 00082 cout << endl; 00083 00084 Exit(a_returnCode); 00085 } 00086 00087 00088 /* 00089 Main program 00090 */ 00091 int main(int argc, char** argv) 00092 { 00093 // No argument, then show help 00094 if (argc==1) ShowHelp(0); 00095 00096 // ============================================================================================================ 00097 // Parameterized variables with their default values 00098 // ============================================================================================================ 00099 00100 // String containing the path to the input data filename provided by the user 00101 string path_to_input_file = ""; 00102 00103 // Datafile mode (histogram or list (default) 00104 int data_mode = MODE_LIST; 00105 00106 // Path to user provided data 00107 string path_to_out_file = ""; 00108 string path_to_data_file = ""; 00109 string path_to_header_file = ""; 00110 string path_to_norm_file = ""; 00111 string path_to_scat_file = ""; 00112 string path_to_rdm_file = ""; 00113 string path_to_acf_file = ""; 00114 FLTNB calib_factor = 1.; 00115 00116 string scanner_alias = ""; 00117 string istp_alias = ""; 00118 00119 // Verbosity 00120 int vb = 1; 00121 00122 // Input datafile is histogram 00123 bool is_histogram_flag = 0; 00124 00125 00126 // ============================================================================================================ 00127 // Read command-line parameters 00128 // ============================================================================================================ 00129 00130 for (int i=1; i<argc; i++) 00131 { 00132 string option = (string)argv[i]; 00133 00134 // Provide an histogram datafile to convert 00135 if (option=="-ih") 00136 { 00137 if (!path_to_input_file.empty()) 00138 { 00139 Cerr("***** castor-datafileConversionEx :: the following file names have already been provided (either -ih or -il option must be used): " << endl); 00140 Cout(path_to_input_file << endl); 00141 Exit(EXIT_FAILURE); 00142 } 00143 else 00144 { 00145 if (argv[i+1] == NULL) 00146 { 00147 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00148 Exit(EXIT_FAILURE); 00149 } 00150 else 00151 path_to_input_file = (string)argv[i+1]; 00152 00153 data_mode = MODE_HISTOGRAM; 00154 i++; 00155 } 00156 } 00157 00158 // Provide a list-mode datafile to convert 00159 if (option=="-il") 00160 { 00161 if (!path_to_input_file.empty()) 00162 { 00163 Cerr("***** castor-datafileConversionEx :: the following file names have already been provided (either -ih or -il option must be used): " << endl); 00164 Cout(path_to_input_file << endl); 00165 Exit(EXIT_FAILURE); 00166 } 00167 else 00168 { 00169 if (argv[i+1] == NULL) 00170 { 00171 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00172 Exit(EXIT_FAILURE); 00173 } 00174 else 00175 path_to_input_file = (string)argv[i+1]; 00176 00177 data_mode = MODE_LIST; 00178 i++; 00179 } 00180 } 00181 00182 // Output CASToR datafile 00183 else if (option=="-o") 00184 { 00185 if (argv[i+1] == NULL) 00186 { 00187 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00188 Exit(EXIT_FAILURE); 00189 } 00190 else 00191 path_to_out_file = (string)argv[i+1]; 00192 i++; 00193 } 00194 00195 // Scanner alias 00196 else if (option=="-s") 00197 { 00198 if (argv[i+1] == NULL) 00199 { 00200 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00201 Exit(EXIT_FAILURE); 00202 } 00203 else 00204 scanner_alias = (string)argv[i+1]; 00205 i++; 00206 } 00207 00208 // Provide a norm correction file 00209 else if (option=="-nc") 00210 { 00211 if (argv[i+1] == NULL) 00212 { 00213 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00214 Exit(EXIT_FAILURE); 00215 } 00216 else 00217 path_to_norm_file = (string)argv[i+1]; 00218 00219 i++; 00220 } 00221 00222 // Provide a scatter correction file 00223 else if (option=="-sc") 00224 { 00225 if (argv[i+1] == NULL) 00226 { 00227 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00228 Exit(EXIT_FAILURE); 00229 } 00230 else 00231 path_to_scat_file = (string)argv[i+1]; 00232 00233 i++; 00234 } 00235 00236 // Provide a rdm correction file 00237 else if (option=="-rc") 00238 { 00239 if (argv[i+1] == NULL) 00240 { 00241 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00242 Exit(EXIT_FAILURE); 00243 } 00244 else 00245 path_to_rdm_file = (string)argv[i+1]; 00246 00247 i++; 00248 } 00249 00250 // Provide an acf correction file 00251 else if (option=="-ac") 00252 { 00253 if (argv[i+1] == NULL) 00254 { 00255 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00256 Exit(EXIT_FAILURE); 00257 } 00258 else 00259 path_to_acf_file = (string)argv[i+1]; 00260 00261 i++; 00262 } 00263 00264 // Provide a calibration factor 00265 else if (option=="-cf") 00266 { 00267 if (argv[i+1] == NULL) 00268 { 00269 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00270 Exit(EXIT_FAILURE); 00271 } 00272 else 00273 { 00274 string val_str = (string)argv[i+1]; 00275 if(ConvertFromString(val_str, &calib_factor) ) 00276 { 00277 Cerr("***** castor-datafileConversionEx :: Exception when trying to read'" << val_str << " for option: " << option << endl); 00278 Exit(EXIT_FAILURE); 00279 } 00280 } 00281 i++; 00282 } 00283 00284 // Provide an isotope alias 00285 else if (option=="-ist") 00286 { 00287 if (argv[i+1] == NULL) 00288 { 00289 Cerr("***** castor-datafileConversionEx :: argument missing for option: " << option << endl); 00290 Exit(EXIT_FAILURE); 00291 } 00292 else 00293 istp_alias = (string)argv[i+1]; 00294 00295 i++; 00296 } 00297 00298 // Verbosity level 00299 else if (option=="-vb") 00300 { 00301 if (i>=argc-1) 00302 { 00303 Cerr("***** castor-datafileConversionEx :: Argument missing for option: " << option << endl); 00304 Exit(EXIT_FAILURE); 00305 } 00306 vb = atoi(argv[i+1]); 00307 i++; 00308 } 00309 00310 else 00311 { 00312 Cerr("***** castor-datafileConversionEx :: Unknown option '" << option << "' !" << endl); 00313 Exit(EXIT_FAILURE); 00314 } 00315 } 00316 00317 00318 // ============================================================================================================ 00319 // Basic initialization checks (minimal initializations mandatory for the next steps) 00320 // ============================================================================================================ 00321 00322 // data files 00323 if (path_to_input_file.empty() ) 00324 { 00325 Cerr("***** castor-datafileConversionEx :: Please provide at least one data filename (-ih for histogram, -il for list-mode)" << endl); 00326 cout << " -input filename : give an input datafile (no default)." << endl; 00327 Exit(EXIT_FAILURE); 00328 } 00329 00330 // output directory 00331 if (path_to_out_file.empty() ) 00332 { 00333 Cerr("***** castor-datafileConversionEx :: Please provide the output file name (-o)" << endl); 00334 Exit(EXIT_FAILURE); 00335 } 00336 00337 // scanner 00338 if (scanner_alias.empty()) 00339 { 00340 Cerr("***** castor-datafileConversionEx :: Please provide a scanner alias (-s) :" << endl); 00341 Cout(" It must correspond to a .geom or .hscan (user-made LUT) file in the config/scanner directory." << endl); 00342 ShowHelp(0); 00343 Exit(EXIT_FAILURE); 00344 } 00345 00346 // Throw error by default. 00347 // These two lines must be erased when implementing any datafile converter! 00348 //Cerr("castor-GATERootToCastor :: This script is only for demonstration purposes " << endl); 00349 //Exit(EXIT_FAILURE); 00350 00351 00352 // ============================================================================================================ 00353 // SOutputManager object initialisation: 00354 // (Used for log file) 00355 // ============================================================================================================ 00356 00357 sOutputManager* p_outputManager = sOutputManager::GetInstance(); 00358 // Set verbose level 00359 p_outputManager->SetVerbose(vb); 00360 // Set MPI rank 00361 p_outputManager->SetMPIRank(0); 00362 00363 // Set path to the config directory 00364 if (p_outputManager->CheckConfigDir("")) 00365 { 00366 Cerr("***** castor-datafileConversionEx :: A problem occured while checking for the config directory path !" << endl); 00367 Exit(EXIT_FAILURE); 00368 } 00369 // Initialize output directory and base name 00370 if (p_outputManager->InitOutputDirectory(path_to_out_file, "")) 00371 { 00372 Cerr("***** castor-datafileConversionEx :: A problem occured while initializing output directory !" << endl); 00373 00374 Exit(EXIT_FAILURE); 00375 } 00376 // Log command line 00377 if (p_outputManager->LogCommandLine(argc,argv)) 00378 { 00379 Cerr("***** castor-datafileConversionEx :: A problem occured while logging command line arguments !" << endl); 00380 Exit(EXIT_FAILURE); 00381 } 00382 00383 // Output progression options 00384 cout << std::fixed; 00385 cout << std::setprecision(2); 00386 00387 // ============================================================================================================ 00388 // ScannerManager object initialization 00389 // (Basic initialization to locate scanner geometry files) 00390 // ============================================================================================================ 00391 00392 sScannerManager* p_scannermanager = sScannerManager::GetInstance(); 00393 p_scannermanager->SetVerbose(vb); 00394 00395 // Check if the scanner exists and get the name from ScannerManager 00396 scanner_alias = (scanner_alias.find(OS_SEP)) ? 00397 scanner_alias.substr(scanner_alias.find_last_of(OS_SEP)+1) : 00398 scanner_alias; 00399 00400 if(p_scannermanager->FindScannerSystem(scanner_alias) ) 00401 { 00402 Cerr("**** castor-datafileConversionEx :: A problem occurred while searching for scanner system !" << endl); 00403 Exit(EXIT_FAILURE); 00404 } 00405 00406 // Build output file names 00407 path_to_data_file = path_to_out_file + ".Cdf"; 00408 path_to_header_file = path_to_out_file + ".Cdh"; 00409 00410 // ============================================================================================================ 00411 // Instanciate/Initialize CASToR Datafile 00412 // (Enable/Disable optionnal field to be written) 00413 // ============================================================================================================ 00414 00415 // Instantiate & Initialize oImageDimensionsAndQuantification object, required for datafile generation (number of threads) 00416 oImageDimensionsAndQuantification* p_ID = new oImageDimensionsAndQuantification(); 00417 00418 // Instanciate & Initialize iDataFilePET and Event objects 00419 iDataFilePET* Out_data_file = NULL; 00420 vEvent* Event = NULL; 00421 00422 // ===> This variable set the maximum number of LORs contained in one event 00423 // ===> in the whole dataset (for compression purposes) 00424 // ===> This depends on the dataset to convert 00425 uint16_t max_nb_lines_per_event = 1; 00426 00427 // Create datafile and set variables 00428 Out_data_file = new iDataFilePET(); 00429 00430 // Set the image oImageDimensionsAndQuantification object (required for initialization) 00431 Out_data_file->SetImageDimensionsAndQuantification(p_ID); 00432 00433 // Set path to header datafile 00434 Out_data_file->SetHeaderDataFileName(path_to_header_file); 00435 00436 // Just required for initialization 00437 Out_data_file->SetPercentageLoad(0); 00438 00439 // Set verbosity 00440 Out_data_file->SetVerbose(vb); 00441 00442 // Set data mode (list/hitogram) 00443 Out_data_file->SetDataMode(data_mode); 00444 00445 // Set modality 00446 Out_data_file->SetDataType(TYPE_PET); 00447 00448 // Set maximum number of LORs contained in one event 00449 Out_data_file->SetMaxNberLinesPerEvent(max_nb_lines_per_event); 00450 00451 // Set calibration factor 00452 Out_data_file->SetCalibrationFactor(calib_factor); 00453 00454 // Set isotope 00455 if(!istp_alias.empty()) 00456 { 00457 if(p_ID->SetPETIsotope(0, istp_alias) ) 00458 { 00459 Cerr("**** castor-datafileConversionEx :: A problem occurred while checking isotope name !" << endl); 00460 Exit(EXIT_FAILURE); 00461 } 00462 Out_data_file->SetIsotope(istp_alias); 00463 } 00464 00465 // Allocate the Event structure (containing all information for a listmode 00466 // or histogram event, depending on data mode) 00467 if(data_mode == MODE_HISTOGRAM) 00468 { 00469 // Instanciate list-mode Event 00470 Event = new iEventHistoPET(); 00471 Event->SetNbLines(max_nb_lines_per_event); 00472 Event->AllocateID(); 00473 } 00474 else 00475 { 00476 // Instanciate histogram Event 00477 Event = new iEventListPET(); 00478 Event->SetNbLines(max_nb_lines_per_event); 00479 Event->AllocateID(); 00480 } 00481 00482 00483 // ============================================================================================================ 00484 // ===> Read sinogram correction file(s) (if any) and set datafile flag for each correction 00485 // ============================================================================================================ 00486 00487 // ===> Compute sinogram dimensions. 00488 // ===> Values must be adjusted to the dataset to convert. 00489 uint32_t nbins = 1; 00490 uint32_t nangles = 1; 00491 uint32_t nsinos = 1; 00492 uint32_t nb_cor_factors = nbins 00493 * nangles 00494 * nsinos; 00495 00496 // --- Normalization correction factors --- 00497 00498 // FLTNBDATA is either a float (default) or double type. 00499 // Defined in /management/gVariables.hh 00500 FLTNBDATA *p_norm = NULL; 00501 bool norm_flag = false; 00502 00503 if( !path_to_norm_file.empty() ) 00504 { 00505 if(vb>=2) Cout("--> Reading the norm correction factors file = " << path_to_norm_file << "..." << endl); 00506 ifstream ifile(path_to_norm_file.c_str(), ios::binary | ios::in); 00507 00508 if (!ifile) 00509 { 00510 Cerr("***** castor-datafileConversionEx :: Error reading normalization factor file: "<< path_to_norm_file <<" !" << endl); 00511 Exit(EXIT_FAILURE); 00512 } 00513 else 00514 { 00515 // Allocate buffer storing the normalisation 00516 p_norm = new FLTNBDATA [nb_cor_factors * sizeof( FLTNBDATA )]; 00517 ifile.read(reinterpret_cast < char* > (p_norm), nb_cor_factors * sizeof( FLTNBDATA)); 00518 ifile.close(); 00519 } 00520 00521 // ===> Set normalization correction flag in datafile objet 00522 // ===> (required to indicate that normalization correction will be provided) 00523 Out_data_file->SetNormCorrectionFlagOn(); 00524 norm_flag = true; 00525 } 00526 00527 // --- Scatter correction factors file --- 00528 FLTNBDATA *p_scat = NULL; 00529 bool scat_flag = false; 00530 00531 if( !path_to_scat_file.empty() ) 00532 { 00533 if(vb>=2) Cout("--> Reading the scatter correction factors file = " << path_to_scat_file << "..." << endl); 00534 ifstream ifile(path_to_scat_file.c_str(), ios::binary | ios::in); 00535 00536 if (!ifile) 00537 { 00538 Cerr("***** castor-datafileConversionEx :: Error reading scatter factor file: "<< path_to_scat_file <<" !" << endl); 00539 Exit(EXIT_FAILURE); 00540 } 00541 else 00542 { 00543 // Allocating buffer storing the scatter 00544 p_scat = new FLTNBDATA [nb_cor_factors * sizeof( FLTNBDATA )]; 00545 ifile.read(reinterpret_cast < char* > (p_scat), nb_cor_factors * sizeof( FLTNBDATA)); 00546 ifile.close(); 00547 } 00548 00549 // ===> Set scatter correction flag in datafile objet 00550 // ===> (required to indicate that scatter correction will be provided) 00551 Out_data_file->SetScatterCorrectionFlagOn(); 00552 scat_flag = true; 00553 } 00554 00555 // --- Random correction factors file --- 00556 FLTNBDATA *p_rdm = NULL; 00557 bool rdm_flag = false; 00558 00559 if( !path_to_rdm_file.empty() ) 00560 { 00561 if(vb>=2) Cout("--> Reading the random correction factors file = " << path_to_rdm_file << "..." << endl); 00562 ifstream ifile(path_to_rdm_file.c_str(), ios::binary | ios::in); 00563 00564 if (!ifile) 00565 { 00566 Cerr("***** castor-datafileConversionEx :: Error reading random correction factor file: "<< path_to_rdm_file <<" !" << endl); 00567 Exit(EXIT_FAILURE); 00568 } 00569 else 00570 { 00571 // Allocating buffer storing the random 00572 p_rdm = new FLTNBDATA [nb_cor_factors * sizeof( FLTNBDATA )]; 00573 ifile.read(reinterpret_cast < char* > (p_rdm), nb_cor_factors * sizeof( FLTNBDATA)); 00574 ifile.close(); 00575 } 00576 00577 // ===> Set random correction flag in datafile objet 00578 // ===> (required to indicate that random correction will be provided) 00579 Out_data_file->SetRandomCorrectionFlagOn(); 00580 rdm_flag = true; 00581 } 00582 00583 // --- Attenuation correction factors file --- 00584 FLTNBDATA *p_acf = NULL; 00585 bool acf_flag = false; 00586 00587 if( !path_to_acf_file.empty() ) 00588 { 00589 if(vb>=2) Cout("--> Reading the attenuation correction factors file = " << path_to_acf_file << "..." << endl); 00590 ifstream ifile(path_to_acf_file.c_str(), ios::binary | ios::in); 00591 00592 if (!ifile) 00593 { 00594 Cerr("***** castor-datafileConversionEx :: Error reading attenuation correction factor file: "<< path_to_acf_file <<" !" << endl); 00595 Exit(EXIT_FAILURE); 00596 } 00597 else 00598 { 00599 // Allocating buffer storing the acf 00600 p_acf = new FLTNBDATA [nb_cor_factors * sizeof( FLTNBDATA )]; 00601 ifile.read(reinterpret_cast < char* > (p_acf), nb_cor_factors* sizeof( FLTNBDATA)); 00602 ifile.close(); 00603 } 00604 00605 // ===> Set attenuation correction flag in datafile objet 00606 // ===> (required to indicate that attenuation correction will be provided) 00607 Out_data_file->SetAtnCorrectionFlagOn(); 00608 acf_flag = true; 00609 } 00610 00611 // Initialize datafile object 00612 Out_data_file->PROJ_InitFile(); 00613 00614 // Compute the event size relatively to the mandatory 00615 // and optionnal (correction factors) elements 00616 Out_data_file->ComputeSizeEvent(); 00617 00618 // Init Datafile 00619 Out_data_file->PrepareDataFile(); 00620 00621 00622 // ============================================================================================================ 00623 // *********************************************** CONVERSION ************************************************* 00624 // ============================================================================================================ 00625 Cout(" --- Start conversion of datafile : " << path_to_input_file << endl 00626 << " CASToR output header datafile: " << path_to_header_file << endl 00627 << " CASToR output binary datafile: " << path_to_data_file << endl << endl); 00628 00629 00630 // ============================================================================================================ 00631 // Init variables 00632 // ============================================================================================================ 00633 00634 // Index for progression printing 00635 uint32_t printing_index = 0; 00636 00637 // ===> Recover the total number of events (number of coincidences or 00638 // ===> histogram bins) from the dataset metadata 00639 uint32_t nEvents = 0; 00640 00641 00642 // ===> Acquisition duration and start time in seconds 00643 // ===> This should be recovered from the dataset metadata) 00644 FLTNB start_time_sec = 0.; 00645 FLTNB stop_time_sec = 1.; 00646 FLTNB duration_sec = stop_time_sec - start_time_sec; 00647 00648 // ===> Scanner variables (blocs, modules, submodules, crystals etc..) 00649 // ===> depending on the system. 00650 uint32_t nRsectorsPerRing = 70; 00651 uint32_t nb_crystals_trs = 9, 00652 nb_crystals_axl = 6; 00653 uint32_t nb_crystals_per_ring = nRsectorsPerRing 00654 * nb_crystals_trs; 00655 00656 // ===> Variables to recover indices of the two crystals for each potential 00657 // ===> LORs contributing to the event 00658 uint32_t *castor_id1 = new uint32_t[max_nb_lines_per_event]; 00659 uint32_t *castor_id2 = new uint32_t[max_nb_lines_per_event]; 00660 00661 // -------------------------------- 00662 // Initialization of root variables 00663 // The following lines are just required for the demonstration dataset used in this sample code, 00664 // Must be erased when adjusting this code to the conversion of any dataset 00665 00666 // ROOT objects/variables 00667 int32_t rSectorID1, rSectorID2; 00668 int32_t moduleID1, moduleID2; 00669 int32_t crystalID1, crystalID2; 00670 00671 #ifdef CASTOR_ROOT 00672 TApplication *Tapp = new TApplication("tapp",0,0); 00673 TTree* Coincidences = new TTree; 00674 TFile* Tfile_root = new TFile(path_to_input_file.c_str(),"READ","ROOT file with histograms"); 00675 Coincidences = (TTree*)Tfile_root->Get("Coincidences"); 00676 nEvents += Coincidences->GetEntries(); 00677 00678 //------------- ROOT branches 00679 double_t time1=0; 00680 Coincidences->SetBranchAddress("time1",&time1); 00681 Coincidences->SetBranchAddress("rsectorID1",&rSectorID1); 00682 Coincidences->SetBranchAddress("moduleID1",&moduleID1); 00683 Coincidences->SetBranchAddress("crystalID1",&crystalID1); 00684 00685 double_t time2=0; 00686 Coincidences->SetBranchAddress("time2",&time2); 00687 Coincidences->SetBranchAddress("rsectorID2",&rSectorID2); 00688 Coincidences->SetBranchAddress("moduleID2",&moduleID2); 00689 Coincidences->SetBranchAddress("crystalID2",&crystalID2); 00690 00691 // Recover acquisition start time from the first event 00692 Coincidences->GetEntry(0); 00693 #endif 00694 // -------------------------------- 00695 00696 // Just for progression output 00697 uint32_t printing_ratio = (nEvents>1000) ? 1000 : nEvents/10; 00698 00699 // Loop on the data elements (histogram bins or LORs) 00700 for (uint32_t i=0; i<nEvents ; i++) 00701 { 00702 // ===> Recover event value (list-mode, or dynamic histogram dataset) 00703 uint32_t time_ms = 0; 00704 stop_time_sec = (FLTNB)time_ms/1000; 00705 00706 // ===> Recover the number of lines in the event (if compression). 00707 int nb_lines_in_event = 1; 00708 00709 // ----------------------------------------------------------------- 00710 // The following lines are just required for the demonstration dataset used in this sample code. 00711 // Must be erased/adapted when adjusting this code to the conversion of any dataset 00712 #ifdef CASTOR_ROOT 00713 Coincidences->GetEntry(i); 00714 time_ms = time1; 00715 #endif 00716 00717 // ===> Compute CASToR indexation depending on the system layout 00718 // ===> This directly depends on the type of scanner file used to 00719 // ===> described the system geometry. 00720 // ===> See documentation for more information 00721 uint32_t crystals_trs_id1 = (uint32_t)(crystalID1/nb_crystals_trs); 00722 00723 uint32_t ringID1 = moduleID1*nb_crystals_axl 00724 + crystals_trs_id1; 00725 00726 uint32_t crystals_trs_id2 = (uint32_t)(crystalID2/nb_crystals_trs); 00727 00728 uint32_t ringID2 = moduleID2*nb_crystals_axl 00729 + crystals_trs_id2; 00730 00731 // Loop on the LORs contained in the event 00732 for(int l=0 ; l<nb_lines_in_event ; l++) 00733 { 00734 // Recover the CASToR ID position of the two crystals according to 00735 // the LUT corresponding to this scanner 00736 castor_id1[l] = ringID1*nb_crystals_per_ring 00737 + rSectorID1*nb_crystals_trs 00738 + crystalID1%nb_crystals_trs ; 00739 00740 castor_id2[l] = ringID2*nb_crystals_per_ring 00741 + rSectorID2*nb_crystals_trs 00742 + crystalID2%nb_crystals_trs ; 00743 } 00744 // ----------------------------------------------------------------- 00745 00746 // --- Histo-mode event --- 00747 if(is_histogram_flag) // histogram mode event 00748 { 00749 // Set the number of lines in this event 00750 Event->SetNbLines(nb_lines_in_event); 00751 00752 // Set CASToR IDs in the event structure for each contributing line 00753 for(int l=0 ; l<nb_lines_in_event ; l++) 00754 { 00755 Event->SetID1(l, castor_id1[l]); // 1st argument : line, 2nd argument : index 00756 Event->SetID2(l, castor_id2[l]); // 1st argument : line, 2nd argument : index 00757 } 00758 00759 // Set information related to each TOF bin 00760 int nb_tof_bins = 1; 00761 for (int b=0; b<nb_tof_bins; b++) 00762 { 00763 // ===> Recover and set event value (histogram dataset) 00764 FLTNB event_value = 1.; 00765 Event->SetEventValue(b, event_value); // 1st argument : TOF bin, 2nd argument : event value 00766 00767 // ===> Set scatter correction factor (if any) 00768 // ===> Require computation of the scatter correction value for 00769 // ===> this LOR from the factors loaded in the "p_scat" buffer 00770 if(scat_flag) 00771 { 00772 FLTNB scat_correction_coeff = 0.; // computed from p_scat 00773 ((iEventHistoPET*)Event)->SetScatterRate(b, scat_correction_coeff); 00774 } 00775 } 00776 00777 // ===> Set timestamp of the event in ms (global timestamp for histogram) 00778 // ===> i.e : 0, or timestamp of a dynamic frame 00779 Event->SetTimeInMs(time_ms); // uint32_t 00780 00781 // ===> Set random correction factor (if any) 00782 // ===> Require computation of the random correction value for 00783 // ===> this LOR from the factors loaded in the "p_rdm" buffer 00784 if(rdm_flag) 00785 { 00786 FLTNB rdm_correction_coeff = 0.; // computed from p_rdm 00787 ((iEventPET*)Event)->SetRandomRate(rdm_correction_coeff); 00788 } 00789 00790 // ===> Set normalization correction factor (if any) 00791 // ===> Require computation of the norm correction value for 00792 // ===> this LOR from the factors loaded in the "p_norm" buffer 00793 if(norm_flag) 00794 { 00795 FLTNB norm_correction_coeff = 0.; // computed from p_norm 00796 ((iEventPET*)Event)->SetNormalizationFactor(norm_correction_coeff); 00797 } 00798 00799 // ===> Set attenuation correction factor (if any) 00800 // ===> Require computation of the atn correction value for 00801 // ===> this LOR from the factors loaded in the "p_acf" buffer 00802 if(acf_flag) 00803 { 00804 FLTNB atn_correction_coeff = 0.; // computed from p_acf 00805 ((iEventPET*)Event)->SetAttenuationCorrectionFactor(atn_correction_coeff); 00806 } 00807 00808 // Write event in the datafile 00809 // ('0' argument is just thread-related. No change required) 00810 Out_data_file->PROJ_WriteEvent(Event, 0); 00811 } 00812 00813 // --- List-mode event --- 00814 else 00815 { 00816 // Set the number of lines in this event 00817 Event->SetNbLines(nb_lines_in_event); 00818 00819 // Set CASToR ID in the event structure 00820 for(int l=0 ; l<nb_lines_in_event ; l++) 00821 { 00822 Event->SetID1(l, castor_id1[l]); // 1st argument : line, 2nd argument :index 00823 Event->SetID2(l, castor_id2[l]); // 1st argument : line, 2nd argument :index 00824 } 00825 00826 // ===> Set scatter correction factor (if any) 00827 // ===> Require computation of the scatter correction value for 00828 // ===> this LOR from the factors loaded in the "p_scat" buffer 00829 if(scat_flag) 00830 { 00831 FLTNB scat_correction_coeff = 0.; // computed from p_scat 00832 // First argument '0' by default for list-mode (it is used for histograms with TOF bins). 00833 ((iEventListPET*)Event)->SetScatterRate(0, scat_correction_coeff); 00834 } 00835 00836 // ===> Set timestamp of the event in ms 00837 Event->SetTimeInMs(time_ms); // uint32_t 00838 00839 // ===> Set random correction factor (if any) 00840 // ===> Require computation of the random correction value for 00841 // ===> this LOR from the factors loaded in the "p_rdm" buffer 00842 if(rdm_flag) 00843 { 00844 FLTNB rdm_correction_coeff = 0.; // computed from p_rdm 00845 ((iEventPET*)Event)->SetRandomRate(rdm_correction_coeff); 00846 } 00847 00848 // ===> Set normalization correction factor (if any) 00849 // ===> Require computation of the norm correction value for 00850 // ===> this LOR from the factors loaded in the "p_norm" buffer 00851 if(norm_flag) 00852 { 00853 FLTNB norm_correction_coeff = 0.; // computed from p_norm 00854 ((iEventPET*)Event)->SetNormalizationFactor(norm_correction_coeff); 00855 } 00856 00857 // ===> Set attenuation correction factor (if any) 00858 // ===> Require computation of the atn correction value for 00859 // ===> this LOR from the factors loaded in the "p_acf" buffer 00860 if(acf_flag) 00861 { 00862 FLTNB atn_correction_coeff = 0.; // computed from p_acf 00863 ((iEventPET*)Event)->SetAttenuationCorrectionFactor(atn_correction_coeff); 00864 } 00865 00866 // Write event in the datafile 00867 // ('0' argument is just thread-related. No change required) 00868 Out_data_file->PROJ_WriteEvent(Event, 0); 00869 } 00870 00871 // Progression output 00872 if (printing_index%printing_ratio==0) 00873 { 00874 FLTNB percent = ( ((FLTNB)(printing_index+1))/((FLTNB)nEvents) ) * ((FLTNB)100); 00875 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 " 00876 << percent << " \% "; 00877 } 00878 printing_index++; 00879 } 00880 00881 Cout("The simulated dataset contained " << nEvents << " coincidences / bins" << endl); 00882 00883 // Set duration and nb of events in datafile 00884 Out_data_file->SetStartTime(start_time_sec); 00885 Out_data_file->SetDuration(duration_sec); 00886 Out_data_file->SetNbEvents(nEvents); 00887 00888 // Write datafile header 00889 Out_data_file->PROJ_WriteHeader(); 00890 00891 // Write raw datafile 00892 Out_data_file->PROJ_WriteData(); 00893 00894 // Delete tmp writing file (if any) 00895 Out_data_file->PROJ_DeleteTmpDatafile(); 00896 00897 cout << endl; 00898 00899 // ============================================================================================================ 00900 // End 00901 // ============================================================================================================ 00902 00903 // Free Root objects 00904 #ifdef CASTOR_ROOT 00905 delete Coincidences; 00906 delete Tfile_root; 00907 delete Tapp; 00908 #endif 00909 00910 delete[] castor_id1; 00911 delete[] castor_id2; 00912 00913 if(p_acf != NULL) delete[] p_acf; 00914 if(p_rdm != NULL) delete[] p_rdm; 00915 if(p_scat != NULL) delete[] p_scat; 00916 if(p_norm != NULL) delete[] p_norm; 00917 00918 // Delete objects 00919 if (Event) 00920 { 00921 if (is_histogram_flag) 00922 delete (iEventHistoPET*)Event; 00923 else 00924 delete (iEventListPET*)Event; 00925 } 00926 00927 delete Out_data_file; 00928 delete p_ID; 00929 00930 return 0; 00931 }