![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00009 #include "gVariables.hh" 00010 #include "gOptions.hh" 00011 #include "oAnalyticProjection.hh" 00012 #include "oImageConvolverManager.hh" 00013 #include "oImageDimensionsAndQuantification.hh" 00014 #include "iDataFilePET.hh" 00015 #include "iDataFileSPECT.hh" 00016 #include "iDataFileTransmission.hh" 00017 #include "sOutputManager.hh" 00018 #include "sScannerManager.hh" 00019 #include "sRNG.hh" 00020 #include "iScannerPET.hh" 00021 00022 00023 // ===================================================================== 00024 // --------------------------------------------------------------------- 00025 // --------------------------------------------------------------------- 00026 // ===================================================================== 00031 void ShowHelp(int a_returnCode, int a_mpiRank) 00032 { 00033 // Return when using MPI and mpi_rank is not 0 00034 #ifdef CASTOR_MPI 00035 int mpi_rank = 0; 00036 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00037 if (mpi_rank!=0) return; 00038 #endif 00039 00040 // Show help 00041 cout << endl; 00042 cout << "Usage: castor-proj.exe -img path_to_file.hdr -sc name -(f/d)out output [settings]" << endl; 00043 cout << endl; 00044 cout << "[Main options]:" << endl; 00045 cout << " -img path_to_img.hdr : Give the interfile header of the image to project (no default)." << endl; 00046 cout << " -sc scanner_alias : Give the scanner model. It must correspond to the name of one of the file in the scanner repository (w/o file extension)" << endl; 00047 cout << " -fout name : Give the root name for all output files (no default, alternative to -dout)" << endl; 00048 cout << " -dout name : Give the name of the output directory where all output files will be written (no default, alternative to -fout)" << endl; 00049 cout << " -fileType type : Give the type of events that will be generated : ('TYPE_PET' ou '0' = PET, 'TYPE_SPECT' ou '1' = SPECT) (default : TYPE_PET)" << endl; 00050 cout << " -vb : Give the verbosity level, from 0 (no verbose) to above 5 (at the event level) (default: 1)." << endl; 00051 cout << endl; 00052 cout << "[Specific options]:" << endl; 00053 cout << " -help-in : Print out specific help about input settings." << endl; // managed by main 00054 cout << " -help-out : Print out specific help about output settings." << endl; // managed by main 00055 cout << " -help-proj : Print out specific help about projector settings." << endl; // managed by main 00056 cout << " -help-comp : Print out specific help about computing settings." << endl; // managed by main 00057 cout << " -help-misc : Print out specific help about miscellaneous and verbose settings." << endl; // managed by main 00058 cout << " -help-dynamic : Print out specific help about dynamic options." << endl; // managed by main 00059 cout << " -help-pet : Print out specific help about PET settings for analytic projection." << endl; // managed by main 00060 cout << " -help-spect : Print out specific help about SPECT settings for analytic projection." << endl; // managed by main 00061 cout << endl; 00062 cout << "[Implemented Modules]:" << endl; 00063 cout << " -help-scan : Show the list of all scanners from the configuration directory." << endl; // managed by sScannerManager 00064 cout << " -help-projm : Show the list and description of all implemented projectors." << endl; // managed by oProjectorManager 00065 cout << endl; 00066 cout << " --help,-h,-help : Print out this help page." << endl; // managed by main 00067 cout << endl; 00068 #ifdef CASTOR_MPI 00069 cout << " Compiled with MPI" << endl; 00070 #endif 00071 #ifdef CASTOR_OMP 00072 cout << " Compiled with OpenMP" << endl; 00073 #endif 00074 #ifdef CASTOR_GPU 00075 cout << " Compiled for GPU" << endl; 00076 #endif 00077 #if defined(CASTOR_OMP) || defined(CASTOR_MPI) || defined(CASTOR_GPU) 00078 cout << endl; 00079 #endif 00080 #ifdef BUILD_DATE 00081 cout << " Build date: " << BUILD_DATE << endl; 00082 cout << endl; 00083 #endif 00084 Exit(a_returnCode); 00085 } 00086 00087 00088 00089 00090 // ===================================================================== 00091 // --------------------------------------------------------------------- 00092 // --------------------------------------------------------------------- 00093 // ===================================================================== 00098 void ShowHelpInput() 00099 { 00100 // Return when using MPI and mpi_rank is not 0 00101 #ifdef CASTOR_MPI 00102 int mpi_rank = 0; 00103 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00104 if (mpi_rank!=0) return; 00105 #endif 00106 // Show help 00107 cout << endl; 00108 cout << "[Input settings]:" << endl; 00109 cout << endl; 00110 cout << " -off x,y,z : Give the offset of the field-of-view in each dimension, in mm (default: 0.,0.,0.)." << endl; 00111 cout << endl; 00112 cout << " -atn path_to_img.hdr : Give the interfile header of an input attenuation image (unit has to be cm-1) (default: uniform value)." << endl; 00113 cout << endl; 00114 cout << " -help-in : Print out this help." << endl; 00115 cout << endl; 00116 } 00117 00118 00119 00120 // ===================================================================== 00121 // --------------------------------------------------------------------- 00122 // --------------------------------------------------------------------- 00123 // ===================================================================== 00128 void ShowHelpOutput() 00129 { 00130 // Return when using MPI and mpi_rank is not 0 00131 #ifdef CASTOR_MPI 00132 int mpi_rank = 0; 00133 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00134 if (mpi_rank!=0) return; 00135 #endif 00136 // Show help 00137 cout << endl; 00138 cout << "[Output settings]:" << endl; 00139 cout << endl; 00140 cout << " -fout name : Give the root name for all output files. All output files will be written as 'name_suffix.ext'." << endl; 00141 cout << " So the provided name should not end with '.' or '/' character. (no default, alternative to -dout)" << endl; 00142 cout << " -dout name : Give the name of the output directory where all output files will be written. All files will also" << endl; 00143 cout << " be prefixed by the name of the directory. The provided name should not end with '.' or '/' character." << endl; 00144 cout << " (no default, alternative to -fout)" << endl; 00145 cout << endl; 00146 cout << " -olut : If a scanner LUT (geometry information) is computed from a .geom file, it will be written on disk in the scanner repository" << endl; 00147 cout << endl; 00148 cout << " -help-out : Print out this help." << endl; 00149 cout << endl; 00150 } 00151 00152 00153 00154 // ===================================================================== 00155 // --------------------------------------------------------------------- 00156 // --------------------------------------------------------------------- 00157 // ===================================================================== 00162 void ShowHelpProj() 00163 { 00164 // Return when using MPI and mpi_rank is not 0 00165 #ifdef CASTOR_MPI 00166 int mpi_rank = 0; 00167 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00168 if (mpi_rank!=0) return; 00169 #endif 00170 // Show help 00171 cout << "[Projector settings]:" << endl; 00172 cout << endl; 00173 cout << " -proj param : Give the projector to be used for both forward and backward projections, along with a configuration file (proj:file.conf)" << endl; 00174 cout << " or the list of parameters associated to the projector (proj,param1,param2,...). If the projector only is specified, the" << endl; 00175 cout << " default configuration file is used. By default, the Siddon projector is used." << endl; 00176 cout << endl; 00177 cout << " -proj-comp : Give the strategy for projection line computation. Here are the three different strategies that can be used:" << endl; 00178 cout << " 1 : Image-based system matrix elements storage: The voxels weights are added in a matrix representing the whole image, so" << endl; 00179 cout << " the addition of a new line to the previous ones is straightforward only by adding the weights to the corresponding voxels." << endl; 00180 cout << " As it is insanely long, it can possibly be used for example with extremely complex projectors that makes use of huge number" << endl; 00181 cout << " of ray tracings for a single event, where the list of contributions can become longer than the number of voxels in the image." << endl; 00182 cout << " 2 : Fixed-size list storage of system matrix elements: The voxels are added one by one in two separated lists, one containing voxel" << endl; 00183 cout << " indices and the other voxel weights. When a voxel is added to the oProjectionLine, it is simply pilled-up to the list. The list" << endl; 00184 cout << " has a fixed size which is provided by the EstimateMaxNumberOfVoxelsPerLine() function from the vProjector class. There are no" << endl; 00185 cout << " ckecks at all for possible buffer overflows. This is the fastest strategy and default one." << endl; 00186 cout << " 3 : Adaptative-size list storage of system matrix elements: This is the same as the fixed-size strategy except that the size can be" << endl; 00187 cout << " upgraded if the current number of contributing voxels exceed the list's size. The first allocated size corresponds to the diagonal" << endl; 00188 cout << " of the image. During the first iteration, this size will be upgraded until it will reach a size suitable for all events. Thus it" << endl; 00189 cout << " is a bit slower than the fixed-list strategy during the first iteration, but is optimal with respect to RAM usage." << endl; 00190 cout << endl; 00191 cout << " -help-projm : Print out specific help about projector settings." << endl; // managed by oProjectorManager 00192 cout << endl; 00193 cout << " -help-proj : Print out this help." << endl; // managed by oProjectorManager 00194 cout << endl; 00195 } 00196 00197 00198 00199 // ===================================================================== 00200 // --------------------------------------------------------------------- 00201 // --------------------------------------------------------------------- 00202 // ===================================================================== 00207 void ShowHelpComputation() 00208 { 00209 // Return when using MPI and mpi_rank is not 0 00210 #ifdef CASTOR_MPI 00211 int mpi_rank = 0; 00212 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00213 if (mpi_rank!=0) return; 00214 #endif 00215 // Show help 00216 cout << endl; 00217 cout << "[Computation settings]:" << endl; 00218 cout << endl; 00219 #ifdef CASTOR_OMP 00220 cout << " -th param : Set the number of threads for parallel computation (default: 1). If 0 is given, then the maximum number of available threads is automatically determined." << endl; 00221 cout << " Can also give two parameters separated by a comma (e.g. 16,4), to distinguish between the number of threads for projection and image operations respectively." << endl; 00222 cout << endl; 00223 #endif 00224 #ifdef CASTOR_GPU 00225 cout << " -gpu : Flag to say that we want to use the GPU device (default: use the CPU only)." << endl; 00226 cout << endl; 00227 #endif 00228 cout << " -conv param : Give an image convolver model to be used within the algorithm, along with a configuration file (conv:file.conf) or the" << endl; 00229 cout << " list of parameters associated to the convolver (conv,param1,param2,...). If the convolver only is specified, its default" << endl; 00230 cout << " configuration file is used. By default, no convolver is applied" << endl; 00231 cout << endl; 00232 cout << " -help-conv : Show the list and description of all implemented image convolvers." << endl; // managed by oImageConvolverManager 00233 cout << endl; 00234 cout << " -noise poisson : Give the noise model to be used. Only Poisson noise is implemented for the moment (default : no noise model)" << endl; 00235 cout << endl; 00236 cout << " -rng : Give a seed for the random number generator (should be >=0)" << endl; 00237 cout << endl; 00238 cout << " -help-comp : Print out this help." << endl; 00239 cout << endl; 00240 #ifdef CASTOR_MPI 00241 cout << " Compiled with MPI" << endl; 00242 #endif 00243 #ifdef CASTOR_OMP 00244 cout << " Compiled with OpenMP" << endl; 00245 #endif 00246 #ifdef CASTOR_GPU 00247 cout << " Compiled for GPU" << endl; 00248 #endif 00249 #if defined(CASTOR_OMP) || defined(CASTOR_MPI) || defined(CASTOR_GPU) 00250 cout << endl; 00251 #endif 00252 } 00253 00254 00255 00256 // ===================================================================== 00257 // --------------------------------------------------------------------- 00258 // --------------------------------------------------------------------- 00259 // ===================================================================== 00264 void ShowHelpDynamic() 00265 { 00266 // Return when using MPI and mpi_rank is not 0 00267 #ifdef CASTOR_MPI 00268 int mpi_rank = 0; 00269 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00270 if (mpi_rank!=0) return; 00271 #endif 00272 // Show help 00273 cout << endl; 00274 cout << "[Dynamic settings]:" << endl; 00275 cout << endl; 00276 cout << " -frm list : Give the framing for the reconstruction where 'list' is a list of all frames durations in seconds separated" << endl; 00277 cout << " by commas. To specify a dead frame, add ';0' concatenated to the frame duration. Milliseconds maximum precision." << endl; 00278 cout << " (default: read from the 'image duration (sec)' key in the interfile image header" << endl; 00279 cout << endl; 00280 cout << " (about gated images) : An image contening several respiratory/cardiac 'gates' can be provided to generate one datafile containing histograms of these gated images." << endl; 00281 cout << " The number of gates being directly read from the 'number of respiratory gates' and 'number of cardiac gates' keys in the interfile header" << endl; 00282 cout << endl; 00283 cout << " -help-dynamic : Print out this help." << endl; 00284 cout << endl; 00285 } 00286 00287 00288 00289 // ===================================================================== 00290 // --------------------------------------------------------------------- 00291 // --------------------------------------------------------------------- 00292 // ===================================================================== 00297 void ShowHelpSPECT() 00298 { 00299 // Return when using MPI and mpi_rank is not 0 00300 #ifdef CASTOR_MPI 00301 int mpi_rank = 0; 00302 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00303 if (mpi_rank!=0) return; 00304 #endif 00305 // Show help 00306 cout << "[SPECT settings]:" << endl; 00307 cout << endl; 00308 cout << " -SPECT_bins par : Give two integer parameters corresponding to the transaxial number of bins, separated by a comma (nbBinsX, nbBinsY). " << endl; 00309 cout << " Default = transaxial number of pixels in the SPECT scanner file " << endl; 00310 cout << endl; 00311 cout << " -SPECT_ang par : Give the total number of projections, followed by the first angle and the angle step in float" << endl; 00312 cout << " Format : nb_projs:first_angle,step_angle" << endl; 00313 cout << endl; 00314 cout << " -SPECT_c_ang par : SPECT custom angles : give the total number of projections, followed by the projections angles in float" << endl; 00315 cout << " Format : nb_projs:angle1,angle2,... " << endl; 00316 cout << endl; 00317 cout << " -SPECT_radius par : Give a distance between center of rotation to detectors which will be used for all heads" << endl; 00318 cout << endl; 00319 cout << " -SPECT_rot par : Give the heads rotation direction. Accept two parameters: CW (clockwise, default), CCW (counter-clockwise)" << endl; 00320 cout << endl; 00321 cout << " -help-spect : Print out this specific help." << endl; // managed by oProjectorManager 00322 cout << endl; 00323 } 00324 00325 00326 00327 // ===================================================================== 00328 // --------------------------------------------------------------------- 00329 // --------------------------------------------------------------------- 00330 // ===================================================================== 00335 void ShowHelpPET() 00336 { 00337 // Return when using MPI and mpi_rank is not 0 00338 #ifdef CASTOR_MPI 00339 int mpi_rank = 0; 00340 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00341 if (mpi_rank!=0) return; 00342 #endif 00343 // Show help 00344 cout << "[PET settings]:" << endl; 00345 cout << endl; 00346 cout << " -PET_maxRingDiff param : Set the maximum number of differences between the ring ID of the two crystals (1 = 2D mode) (default = total number of rings in the scanner)" << endl; 00347 cout << endl; 00348 } 00349 00350 00351 00352 // ===================================================================== 00353 // --------------------------------------------------------------------- 00354 // --------------------------------------------------------------------- 00355 // ===================================================================== 00360 void ShowHelpMiscellaneous() 00361 { 00362 // Return when using MPI and mpi_rank is not 0 00363 #ifdef CASTOR_MPI 00364 int mpi_rank = 0; 00365 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00366 if (mpi_rank!=0) return; 00367 #endif 00368 // Show help 00369 cout << endl; 00370 cout << "[Miscellaneous settings]:" << endl; 00371 cout << endl; 00372 cout << " -vb : Give the general verbosity level, from 0 (no verbose) to 5 and above (at the event level) (default: 1)." << endl; 00373 cout << " -vb-algo : Give the verbose level specific to the analytic projection algorithm (default: same as general verbose level)." << endl; 00374 cout << " -vb-proj : Give the verbose level specific to the projector (default: same as general verbose level)." << endl; 00375 cout << " -vb-conv : Give the verbose level specific to the image convolver (default: same as general verbose level)." << endl; 00376 cout << " -vb-scan : Give the verbose level specific to the scanner (default: same as general verbose level)." << endl; 00377 cout << " -vb-data : Give the verbose level specific to the data and image management (default: same as general verbose level)." << endl; 00378 cout << endl; 00379 cout << " -conf : Give the path to the CASToR config directory (default: located through the CASTOR_CONFIG environment variable)." << endl; 00380 cout << endl; 00381 cout << " -help-misc : Print out this help." << endl; 00382 cout << endl; 00383 } 00384 00385 00386 00387 00388 00389 // ===================================================================== 00390 // --------------------------------------------------------------------- 00391 // --------------------------------------------------------------------- 00392 // ===================================================================== 00393 /* 00394 Main program 00395 */ 00396 int main(int argc, char** argv) 00397 { 00398 // ============================================================================================================ 00399 // MPI management 00400 // ============================================================================================================ 00401 int mpi_rank = 0; 00402 //int mpi_size = 1; 00403 #ifdef CASTOR_MPI 00404 int mpi_size = 1; 00405 MPI_Init(&argc, &argv ); 00406 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); 00407 MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); 00408 #endif 00409 00410 // No argument, then show help 00411 if (argc==1) ShowHelp(0,mpi_rank); 00412 if (mpi_rank==0) cout << endl; // TODO : Most stuff is handled with processor 1 at first 00413 00414 // ============================================================================================================ 00415 // Parameterized variables with their default values 00416 // ============================================================================================================ 00417 00418 // Initialization of the voxel and field-of-view values in each spatial dimensions. 00419 // Read from the interfile header of the image 00420 int nb_VoxX=-1, nb_VoxY=-1, nb_VoxZ=-1; 00421 FLTNB fov_SizeX=-1, fov_SizeY=-1, fov_SizeZ=-1; 00422 FLTNB vox_SizeX=-1, vox_SizeY=-1, vox_SizeZ=-1; 00423 // Initialization offset values of the field-of-view in each spatial dimensions 00424 FLTNB offsetX = 0, offsetY = 0, offsetZ = 0; 00425 00426 // Path to an image used as initialization. no default 00427 string path_to_initial_img = ""; 00428 // Path to an image used for the attenuation. default : uniform image 00429 string path_to_attenuation_img = ""; 00430 // Scanner name provided by the user 00431 string scanner_name = ""; 00432 00433 // DataFile mode 00434 int datafile_mode = MODE_HISTOGRAM; 00435 int datafile_type = TYPE_PET; 00436 00437 // Frame list descriptor 00438 string frame_list = ""; 00439 // Number of resp gates 00440 int nb_resp_gates = 1; 00441 // Number of card gates 00442 int nb_card_gates = 1; 00443 00444 // Number of beds position in the data. default : 1 00445 int nb_beds = 1; 00446 00447 // Maximum ring difference between the crystals in a LOR 00448 int max_ring_diff = -1; 00449 00450 // Output directory name. 00451 string path_dout = ""; 00452 // Or root name 00453 string path_fout = ""; 00454 00455 // Write scanner LUT generated by a geom file on disk 00456 bool save_LUT_flag = false; 00457 00458 // String gathering the name of the used noise model (default : none) 00459 string noise_model = ""; 00460 00461 // String gathering the name of the used projector for forward/backward projection (default : Siddon) 00462 string options_projector = "incrementalSiddon"; 00463 00464 // Default projector computation strategy 00465 int projector_computation_strategy = FIXED_LIST_COMPUTATION_STRATEGY; 00466 00467 // String vector gathering the name of the image convolvers with specific options (default: no image convolver) 00468 vector<string> options_image_convolver; 00469 00470 // Using GPU (flag) ->NOTE : default : only CPU 00471 bool gpu_flag = false; 00472 // Size of the buffer used to read the data file (0 by default for projection as it is not used, but should be initialized to use the Datafile objects) 00473 int data_file_percentage_load = 0; 00474 00475 // Verbose level 00476 int verbose_general = 1, 00477 verbose_algo = -1, 00478 verbose_proj = -1, 00479 verbose_conv = -1, 00480 verbose_scan = -1, 00481 verbose_data = -1; 00482 00483 // Path to config directory 00484 string path_to_config_dir = ""; 00485 00486 // Initial seed for random number generator 00487 int64_t seed_RNG = -1; 00488 // Number of threads 00489 string nb_threads = "1"; 00490 00491 // SPECT specific parameters 00492 00493 // SPECT number of transaxial bins 00494 uint16_t SPECT_nb_bins[2] = {0, 0}; 00495 // SPECT Number of projection angles 00496 uint32_t SPECT_nb_projection_angles = 0; 00497 // SPECT projections angles 00498 FLTNB* SPECT_projection_angles = NULL; 00499 // SPECT first angle and step angle for automatic initialization 00500 FLTNB SPECT_first_angle= -1., SPECT_step_angle=-1.; 00501 // SPECT distance from center of rotation to detectors 00502 FLTNB SPECT_radius = -1.; 00503 // SPECT head rotation orientation 00504 string SPECT_rot = "CW"; 00505 00506 // ============================================================================================================ 00507 // Read command-line parameters 00508 // ============================================================================================================ 00509 00510 for (int i=1; i<argc; i++) 00511 { 00512 string option = (string)argv[i]; 00513 00514 if (option=="-h" || option=="--help" || option=="-help") ShowHelp(0,mpi_rank); 00515 00516 // Specific help for integrated scanners (managed by scanner manager) 00517 else if (option=="-help-scan") 00518 { 00519 if(sScannerManager::GetInstance()->ShowScannersDescription()) 00520 Cerr("***** castor-proj() -> An error occured when trying to output the available scanners from the scanner repository !'" << endl;); 00521 Exit(EXIT_SUCCESS); 00522 } 00523 00524 00525 // Specific help for projector settings (managed by projector children) 00526 else if (option=="-help-projm") 00527 { 00528 // Call specific showHelp function from vProjector children 00529 sAddonManager::GetInstance()->ShowHelpProjector(); 00530 // Call the static showHelp function from vProjector 00531 vProjector::ShowCommonHelp(); 00532 Exit(EXIT_SUCCESS); 00533 } 00534 00535 // Specific help for image convolver settings (managed by image convolver children) 00536 else if (option=="-help-conv") 00537 { 00538 // Call specific showHelp function from vImageConvolver children 00539 sAddonManager::GetInstance()->ShowHelpImageConvolver(); 00540 // Call the static showHelp function from oImageConvolverManager 00541 // Call the static showHelp function from oImageConvolverManager 00542 oImageConvolverManager::ShowCommonHelp(); 00543 Exit(EXIT_SUCCESS); 00544 } 00545 00546 // Specific help for input settings (managed by main) 00547 else if (option=="-help-in") 00548 { 00549 ShowHelpInput(); 00550 Exit(EXIT_SUCCESS); 00551 } 00552 // Specific help for output settings (managed by main) 00553 else if (option=="-help-out") 00554 { 00555 ShowHelpOutput(); 00556 Exit(EXIT_SUCCESS); 00557 } 00558 // Specific help for projector settings (managed by main) 00559 else if (option=="-help-proj") 00560 { 00561 ShowHelpProj(); 00562 Exit(EXIT_SUCCESS); 00563 } 00564 // Specific help for computation settings (managed by main) 00565 else if (option=="-help-comp") 00566 { 00567 ShowHelpComputation(); 00568 Exit(EXIT_SUCCESS); 00569 } 00570 // Specific help for dynamic settings (managed by main) 00571 else if (option=="-help-dynamic") 00572 { 00573 ShowHelpDynamic(); 00574 Exit(EXIT_SUCCESS); 00575 } 00576 // Specific help for PET projection settings (managed by main) 00577 else if (option=="-help-pet") 00578 { 00579 ShowHelpPET(); 00580 Exit(EXIT_SUCCESS); 00581 } 00582 // Specific help for SPECT projection settings (managed by main) 00583 else if (option=="-help-spect") 00584 { 00585 ShowHelpSPECT(); 00586 Exit(EXIT_SUCCESS); 00587 } 00588 // Specific help for miscellaneous settings (managed by main) 00589 else if (option=="-help-misc") 00590 { 00591 ShowHelpMiscellaneous(); 00592 Exit(EXIT_SUCCESS); 00593 } 00594 00595 else if (option=="-off") 00596 { 00597 if (i>=argc-1) 00598 { 00599 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00600 Exit(EXIT_FAILURE); 00601 } 00602 else 00603 { 00604 FLTNB input[3]; 00605 if (ReadStringOption(argv[i+1], input, 3, ",", option)) 00606 { 00607 Cerr("***** castor-proj() -> Invalid argument " << argv[i+1] << " for option " << option << " !" << endl); 00608 Exit(EXIT_FAILURE); 00609 } 00610 offsetX = input[0]; 00611 offsetY = input[1]; 00612 offsetZ = input[2]; 00613 } 00614 i++; 00615 } 00616 00617 else if (option=="-frm") 00618 { 00619 if (i>=argc-1) 00620 { 00621 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00622 Exit(EXIT_FAILURE); 00623 } 00624 frame_list = (string)argv[i+1]; 00625 i++; 00626 } 00627 00628 // Scanner name 00629 else if (option=="-sc") 00630 { 00631 if (i>=argc-1) 00632 { 00633 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00634 Exit(EXIT_FAILURE); 00635 } 00636 scanner_name = argv[i+1]; 00637 i++; 00638 } 00639 00640 // Image for the initialisation of the algorithm : What should we do in case of multiple frames ? 00641 else if (option=="-img") 00642 { 00643 if (i>=argc-1) 00644 { 00645 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00646 Exit(EXIT_FAILURE); 00647 } 00648 path_to_initial_img = argv[i+1]; 00649 i++; 00650 } 00651 00652 00653 else if (option=="-fileType") 00654 { 00655 if (i>=argc-1) 00656 { 00657 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00658 Exit(EXIT_FAILURE); 00659 } 00660 00661 string input_str = argv[i+1]; 00662 00663 if(!input_str.compare("TYPE_PET") ) 00664 datafile_type = TYPE_PET; 00665 else if(!input_str.compare("TYPE_SPECT") ) 00666 datafile_type = TYPE_SPECT; 00667 else if(!input_str.compare("TYPE_TRANSMISSION") ) 00668 datafile_type = TYPE_TRANSMISSION; 00669 else // check for number 00670 { 00671 if(ConvertFromString(input_str , &datafile_type)) 00672 { 00673 Cerr("***** castor-proj() -> Exception when trying to read provided entry '" << input_str << " for option: " << option << endl); 00674 Exit(EXIT_FAILURE); 00675 } 00676 } 00677 00678 i++; 00679 } 00680 00681 // Image for the attenuation 00682 else if (option=="-atn") 00683 { 00684 if (i>=argc-1) 00685 { 00686 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00687 Exit(EXIT_FAILURE); 00688 } 00689 path_to_attenuation_img = argv[i+1]; 00690 i++; 00691 } 00692 00693 // Number of beds 00694 else if (option=="-bed") 00695 { 00696 if (i>=argc-1) 00697 { 00698 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00699 Exit(EXIT_FAILURE); 00700 } 00701 nb_beds = atoi(argv[i+1]); 00702 i++; 00703 } 00704 00705 // maximum number of differences between the ring ID of the two crystals 00706 else if (option=="-PET_maxRingDiff") 00707 { 00708 if (i>=argc-1) 00709 { 00710 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00711 Exit(EXIT_FAILURE); 00712 } 00713 max_ring_diff = atoi(argv[i+1]); 00714 i++; 00715 } 00716 00717 // Name of the output directory 00718 else if (option=="-dout") // This is a mandatory option 00719 { 00720 if (i>=argc-1) 00721 { 00722 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00723 Exit(EXIT_FAILURE); 00724 } 00725 path_dout = argv[i+1]; 00726 i++; 00727 } 00728 // Base name of the output files 00729 else if (option=="-fout") // This is a mandatory option 00730 { 00731 if (i>=argc-1) 00732 { 00733 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00734 Exit(EXIT_FAILURE); 00735 } 00736 path_fout = argv[i+1]; 00737 i++; 00738 } 00739 // Flag to say that we want to save time basis functions too 00740 else if (option=="-olut") 00741 { 00742 save_LUT_flag = true; 00743 } 00744 else if (option=="-noise") 00745 { 00746 noise_model = argv[i+1]; 00747 i++; 00748 } 00749 00750 else if (option=="-proj") 00751 { 00752 if (i>=argc-1) 00753 { 00754 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00755 Exit(EXIT_FAILURE); 00756 } 00757 options_projector = (string)argv[i+1]; 00758 i++; 00759 } 00760 00761 // Projection line computation strategy 00762 else if (option=="-proj-comp") 00763 { 00764 if (i>=argc-1) 00765 { 00766 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00767 Exit(EXIT_FAILURE); 00768 } 00769 projector_computation_strategy = atoi(argv[i+1]); 00770 i++; 00771 } 00772 00773 // Image convolver settings 00774 else if (option=="-conv") 00775 { 00776 if (i>=argc-1) 00777 { 00778 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00779 Exit(EXIT_FAILURE); 00780 } 00781 string convolver = (string)argv[i+1]; 00782 convolver.append("::forward"); // convolution only allowed before projection 00783 options_image_convolver.push_back(convolver); 00784 i++; 00785 } 00786 00787 00788 // --- Verbosities --- 00789 00790 // General verbosity level 00791 else if (option=="-vb") 00792 { 00793 if (i>=argc-1) 00794 { 00795 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00796 Exit(EXIT_FAILURE); 00797 } 00798 if (ConvertFromString(argv[i+1], &verbose_general)) 00799 { 00800 Cerr("***** castor-proj() -> Exception when trying to read provided verbosity level '" << verbose_general << " for option: " << option << endl); 00801 Exit(EXIT_FAILURE); 00802 } 00803 i++; 00804 } 00805 00806 // Algorithm verbosity level 00807 else if (option=="-vb-algo") 00808 { 00809 if (i>=argc-1) 00810 { 00811 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00812 Exit(EXIT_FAILURE); 00813 } 00814 if (ConvertFromString(argv[i+1], &verbose_algo)) 00815 { 00816 Cerr("***** castor-proj() -> Exception when trying to read provided verbosity level '" << verbose_algo << " for option: " << option << endl); 00817 Exit(EXIT_FAILURE); 00818 } 00819 i++; 00820 } 00821 00822 // Projector verbosity level 00823 else if (option=="-vb-proj") 00824 { 00825 if (i>=argc-1) 00826 { 00827 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00828 Exit(EXIT_FAILURE); 00829 } 00830 if (ConvertFromString(argv[i+1], &verbose_proj)) 00831 { 00832 Cerr("***** castor-proj() -> Exception when trying to read provided verbosity level '" << verbose_proj << " for option: " << option << endl); 00833 Exit(EXIT_FAILURE); 00834 } 00835 i++; 00836 } 00837 00838 // Image convolver verbosity level 00839 else if (option=="-vb-conv") 00840 { 00841 if (i>=argc-1) 00842 { 00843 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00844 Exit(EXIT_FAILURE); 00845 } 00846 if (ConvertFromString(argv[i+1], &verbose_conv)) 00847 { 00848 Cerr("***** castor-proj() -> Exception when trying to read provided verbosity level '" << verbose_conv << " for option: " << option << endl); 00849 Exit(EXIT_FAILURE); 00850 } 00851 i++; 00852 } 00853 00854 // Scanner verbosity level 00855 else if (option=="-vb-scan") 00856 { 00857 if (i>=argc-1) 00858 { 00859 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00860 Exit(EXIT_FAILURE); 00861 } 00862 if (ConvertFromString(argv[i+1], &verbose_scan)) 00863 { 00864 Cerr("***** castor-proj() -> Exception when trying to read provided verbosity level '" << verbose_scan << " for option: " << option << endl); 00865 Exit(EXIT_FAILURE); 00866 } 00867 i++; 00868 } 00869 00870 // Data and image verbosity level 00871 else if (option=="-vb-data") 00872 { 00873 if (i>=argc-1) 00874 { 00875 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00876 Exit(EXIT_FAILURE); 00877 } 00878 if (ConvertFromString(argv[i+1], &verbose_data)) 00879 { 00880 Cerr("***** castor-proj() -> Exception when trying to read provided verbosity level '" << verbose_data << " for option: " << option << endl); 00881 Exit(EXIT_FAILURE); 00882 } 00883 i++; 00884 } 00885 00886 00887 // RNG seed 00888 else if (option=="-rng") 00889 { 00890 if (i>=argc-1) 00891 { 00892 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00893 Exit(EXIT_FAILURE); 00894 } 00895 if(ConvertFromString(argv[i+1], &seed_RNG)) 00896 { 00897 Cerr("***** castor-proj() -> Exception when trying to read provided number '" << seed_RNG << " for option: " << option << endl); 00898 Exit(EXIT_FAILURE); 00899 } 00900 i++; 00901 } 00902 00903 // Path to config directory 00904 else if (option=="-conf") 00905 { 00906 if (i>=argc-1) 00907 { 00908 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00909 Exit(EXIT_FAILURE); 00910 } 00911 path_to_config_dir = (string)argv[i+1]; 00912 i++; 00913 } 00914 00915 #ifdef CASTOR_GPU 00916 else if (option=="-gpu") 00917 { 00918 gpu_flag = 1; 00919 } 00920 #endif 00921 #ifdef CASTOR_OMP 00922 else if (option=="-th") 00923 { 00924 if (i>=argc-1) 00925 { 00926 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00927 Exit(EXIT_FAILURE); 00928 } 00929 nb_threads = (string)argv[i+1]; 00930 i++; 00931 } 00932 #endif 00933 00934 else if (option=="-load") 00935 { 00936 if (i>=argc-1) 00937 { 00938 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00939 Exit(EXIT_FAILURE); 00940 } 00941 data_file_percentage_load = atoi(argv[i+1]); 00942 if (data_file_percentage_load>100 || data_file_percentage_load<0) 00943 { 00944 Cerr("***** castor-proj() -> Incorrect initialization of the size data buffer" << endl); 00945 Cerr(" Number provided: " << argv[i+1] << " is not in the expected [0;100] interval" << endl); 00946 Exit(EXIT_FAILURE); 00947 } 00948 i++; 00949 } 00950 else if (option=="-SPECT_bins") 00951 { 00952 if (i>=argc-1) 00953 { 00954 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00955 Exit(EXIT_FAILURE); 00956 } 00957 00958 if(ReadStringOption(argv[i+1], SPECT_nb_bins, 2, ",", option)) 00959 { 00960 Cerr("***** castor-proj() -> Invalid argument " << argv[i+1] << " for option " << option << " !" << endl); 00961 Cerr("***** castor-proj() -> Expected two uint16 parameters separated by a comma" << endl); 00962 Exit(EXIT_FAILURE); 00963 } 00964 00965 i++; 00966 } 00967 00968 // Custom initialization of SPECT projection angles using a number of projections, and an equivalent number of pre-set projection angles 00969 else if (option=="-SPECT_c_ang") 00970 { 00971 if (i>=argc-1) 00972 { 00973 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 00974 Exit(EXIT_FAILURE); 00975 } 00976 00977 // parsing 00978 string input_str = argv[i+1]; 00979 size_t colon = input_str.find_first_of(":"); 00980 00981 if (colon==string::npos) 00982 { 00983 Cerr("***** castor-proj() -> Parsing error : colon missing in option: " << option << endl); 00984 Exit(EXIT_FAILURE); 00985 } 00986 else 00987 { 00988 if(ConvertFromString(input_str.substr(0,colon), &SPECT_nb_projection_angles)) 00989 { 00990 Cerr("***** castor-proj() -> An error occured while trying to read the number of SPECT projection values " << input_str.substr(0,colon) << " for option " << option << " !" << endl); 00991 Exit(EXIT_FAILURE); 00992 } 00993 00994 SPECT_projection_angles = new FLTNB[SPECT_nb_projection_angles]; 00995 00996 if(ReadStringOption(input_str.substr(colon+1), SPECT_projection_angles, SPECT_nb_projection_angles, ",", option)) 00997 { 00998 Cerr("***** castor-proj() -> Invalid argument " << argv[i+1] << " for option " << option << " !" << endl); 00999 Exit(EXIT_FAILURE); 01000 } 01001 } 01002 i++; 01003 } 01004 01005 // Generic initialization of SPECT projection angles using a first angle, an angle step, and a number of projection 01006 else if (option=="-SPECT_ang") 01007 { 01008 if (i>=argc-1) 01009 { 01010 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 01011 Exit(EXIT_FAILURE); 01012 } 01013 01014 // parsing 01015 string input_str = argv[i+1]; 01016 size_t colon = input_str.find_first_of(":"); 01017 01018 if (colon==string::npos) 01019 { 01020 Cerr("***** castor-proj() -> Parsing error : colon missing in option: " << option << endl); 01021 Exit(EXIT_FAILURE); 01022 } 01023 else 01024 { 01025 if(ConvertFromString(input_str.substr(0,colon), &SPECT_nb_projection_angles)) 01026 { 01027 Cerr("***** castor-proj() -> An error occured while trying to read the number of SPECT projection values " << input_str.substr(0,colon) << " for option " << option << " !" << endl); 01028 Exit(EXIT_FAILURE); 01029 } 01030 01031 FLTNB input[2]; 01032 if(ReadStringOption(input_str.substr(colon+1), input, 2, ",", option)) 01033 { 01034 Cerr("***** castor-proj() -> Invalid argument " << argv[i+1] << " for option " << option << " !" << endl); 01035 Exit(EXIT_FAILURE); 01036 } 01037 else 01038 { 01039 SPECT_first_angle = input[0]; 01040 SPECT_step_angle = input[1]; 01041 } 01042 } 01043 i++; 01044 } 01045 01046 else if (option=="-SPECT_radius") 01047 { 01048 if (i>=argc-1) 01049 { 01050 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 01051 Exit(EXIT_FAILURE); 01052 } 01053 01054 else 01055 { 01056 if(ConvertFromString(argv[i+1], &SPECT_radius) ) 01057 { 01058 Cerr("***** castor-proj() -> An error occured while trying to read the SPECT global radius " << argv[i+1] << " for option " << option << " !" << endl); 01059 Exit(EXIT_FAILURE); 01060 } 01061 } 01062 i++; 01063 } 01064 01065 else if (option=="-SPECT_rot") 01066 { 01067 if (i>=argc-1) 01068 { 01069 Cerr("***** castor-proj() -> Argument missing for option: " << option << endl); 01070 Exit(EXIT_FAILURE); 01071 } 01072 01073 else 01074 { 01075 SPECT_rot = argv[i+1]; 01076 01077 if(SPECT_rot != "CW" && SPECT_rot != "CCW") 01078 { 01079 Cerr("***** castor-proj() -> '" << SPECT_rot <<"' is unknown argument for option " << option << "."); 01080 Cerr(" SPECT rotation direction must be either 'CW' (clockwise) or 'CCW' (counter-clockwise) !" << endl); 01081 Exit(EXIT_FAILURE); 01082 } 01083 } 01084 i++; 01085 } 01086 01087 else 01088 { 01089 if (mpi_rank==0) Cerr("***** castor-proj() -> Unknown option '" << option << "' !" << endl); 01090 Exit(EXIT_FAILURE); 01091 } 01092 01093 } 01094 01095 #ifdef CASTOR_MPI 01096 MPI_Barrier(MPI_COMM_WORLD); 01097 #endif 01098 01099 01100 // Affect specific verbose levels if not set 01101 if (verbose_algo==-1) verbose_algo = verbose_general; 01102 if (verbose_proj==-1) verbose_proj = verbose_general; 01103 if (verbose_conv==-1) verbose_conv = verbose_general; 01104 if (verbose_scan==-1) verbose_scan = verbose_general; 01105 if (verbose_data==-1) verbose_data = verbose_general; 01106 01107 01108 01109 // ============================================================================================================ 01110 // Mandatory checks: the idea is that the most likely problems are checked here, so that we do not apply a 01111 // too much secured approach in the rest of the classes in order to not overload the code 01112 // with checks everywhere, which could also possibly slows down its execution. 01113 // ============================================================================================================ 01114 01115 // Basic initialization checks (minimal initializations mandatory for the next steps) 01116 01117 // scanner 01118 if (scanner_name.empty()) 01119 { 01120 Cerr("***** castor-proj() -> Please provide a scanner name (format: Modality_Constuctor_Model) :" << endl); 01121 Cerr(" -sc scanner_alias: Give the scanner model." << endl); 01122 Cerr(" It must correspond to the name of one of the file in the scanner repository (without file extension)" << endl); 01123 Exit(EXIT_FAILURE); 01124 } 01125 else 01126 { 01127 if (verbose_general>=2) Cout(" selected scanner name: " << scanner_name << endl); 01128 } 01129 01130 // Output files 01131 if (path_fout.empty() && path_dout.empty()) 01132 { 01133 Cerr("***** castor-proj() -> Please provide an output option for output files (-fout or -dout) !" << endl); 01134 Exit(EXIT_FAILURE); 01135 } 01136 // Check that only one option has been provided 01137 if (!path_fout.empty() && !path_dout.empty()) 01138 { 01139 Cerr("***** castor-proj() -> Please provide either output option -fout or -dout but not both !" << endl); 01140 Exit(EXIT_FAILURE); 01141 } 01142 01143 // Image source 01144 if (path_to_initial_img.empty()) 01145 { 01146 Cerr("***** castor-proj() -> -img path_to_img.hdr : Give an interfile header of the image to project" << endl); 01147 Exit(EXIT_FAILURE); 01148 } 01149 else 01150 { 01151 if (verbose_general>=2) Cout(" selected input image path: " << path_to_initial_img << endl); 01152 } 01153 01154 Cout(endl); 01155 01156 01157 01158 01159 // ============================================================================================================ 01160 // Singletons initialization: create here all the needed singletons. 01161 // ============================================================================================================ 01162 01163 // Get user endianness (for interfile) 01164 GetUserEndianness(); 01165 01166 // OutputManager 01167 sOutputManager* p_outputManager; 01168 p_outputManager = sOutputManager::GetInstance(); 01169 p_outputManager->SetVerbose(verbose_general); 01170 01171 // Set path to the config directory 01172 if (p_outputManager->CheckConfigDir(path_to_config_dir)) 01173 { 01174 Cerr("***** castor-proj() -> A problem occured while checking for the config directory path !" << endl); 01175 Exit(EXIT_FAILURE); 01176 } 01177 01178 // Initialize output directory and base name 01179 if (p_outputManager->InitOutputDirectory(path_fout, path_dout)) 01180 { 01181 Cerr("***** castor-proj() -> A problem occured while initializing output directory !" << endl); 01182 Exit(EXIT_FAILURE); 01183 } 01184 // Log command line 01185 if (p_outputManager->LogCommandLine(argc,argv)) 01186 { 01187 Cerr("***** castor-proj() -> A problem occured while logging command line arguments !" << endl); 01188 Exit(EXIT_FAILURE); 01189 } 01190 01191 // ScannerManager 01192 sScannerManager* p_scannermanager = sScannerManager::GetInstance(); 01193 p_scannermanager->SetVerbose(verbose_scan); 01194 p_scannermanager->SetSaveLUTFlag(save_LUT_flag); 01195 01196 01197 01198 // ============================================================================================================ 01199 // Main part of the program: create all elements needed for the algorithm, create the algorithm, initialize it 01200 // and launch it. 01201 // ============================================================================================================ 01202 01203 01204 // ============================================================================================================ 01205 // Generate system geometry 01206 // ============================================================================================================ 01207 01208 scanner_name = (scanner_name.find(OS_SEP)) ? 01209 scanner_name.substr(scanner_name.find_last_of(OS_SEP)+1) : 01210 scanner_name; 01211 01212 if(p_scannermanager->FindScannerSystem(scanner_name) ) 01213 { 01214 Cerr("***** castor-proj() -> A problem occurred while searching for scanner system !" << endl); 01215 Exit(EXIT_FAILURE); 01216 } 01217 01218 if(p_scannermanager->BuildScannerObject() ) //TODO verifier si correct d'envoyer uniquement le chemin du premier DataFile (pb si dataFile ne proviennent pas du même scan ou de la même acquisition) 01219 { 01220 Cerr("***** castor-proj() -> A problem occurred during scanner object construction !" << endl); 01221 Exit(EXIT_FAILURE); 01222 } 01223 01224 if(p_scannermanager->InstantiateScanner() ) 01225 { 01226 Cerr("***** castor-proj() -> A problem occurred while creating Scanner object !" << endl); 01227 Exit(EXIT_FAILURE); 01228 } 01229 01230 01231 //Set the specific mandatory geometry information of the modality using the SetXXXSpecificParameters functions of the scannerManager 01232 if (p_scannermanager->GetScannerType() == SCANNER_PET) 01233 { 01234 if(p_scannermanager->PROJ_SetPETSpecificParameters(max_ring_diff) ) 01235 { 01236 Cerr("***** castor-proj() -> A problem occurred while setting PET geometric parameters !" << endl); 01237 Exit(EXIT_FAILURE); 01238 } 01239 } 01240 else if (p_scannermanager->GetScannerType() == SCANNER_SPECT_CONVERGENT || p_scannermanager->GetScannerType() == SCANNER_SPECT_PINHOLE ) 01241 { 01242 if(p_scannermanager->PROJ_SetSPECTSpecificParameters(SPECT_nb_bins, 01243 SPECT_nb_projection_angles, 01244 SPECT_first_angle, 01245 SPECT_step_angle, 01246 SPECT_projection_angles, 01247 SPECT_radius, 01248 SPECT_rot) ) 01249 { 01250 Cerr("***** castor-proj() -> A problem occurred while setting SPECT geometric parameters !" << endl); 01251 Exit(EXIT_FAILURE); 01252 } 01253 } 01254 01255 01256 if(p_scannermanager->BuildLUT() ) 01257 { 01258 Cerr("***** castor-proj() -> A problem occurred while generating/reading the LUT !" << endl); 01259 Exit(EXIT_FAILURE); 01260 } 01261 01262 01263 // Check the scanner manager parameters and initialize the scanner 01264 if (p_scannermanager->CheckParameters()) 01265 { 01266 Cerr("***** castor-proj() -> A problem occured while checking scanner manager parameters !" << endl); 01267 Exit(EXIT_FAILURE); 01268 } 01269 if (p_scannermanager->Initialize()) 01270 { 01271 Cerr("***** castor-proj() -> A problem occured while initializing scanner !" << endl); 01272 Exit(EXIT_FAILURE); 01273 } 01274 01275 01276 01277 // ============================================================================================================ 01278 // Recover image dimensions informations from the interfile header of the image to project 01279 // ============================================================================================================ 01280 01281 Intf_fields IF; 01282 IntfKeyInitFields(&IF); 01283 if(IntfReadHeader(path_to_initial_img, &IF, verbose_data) ) 01284 { 01285 Cerr("***** castor-proj() -> An error occurred while trying to read the interfile header of file reading file " << path_to_initial_img << " !" << endl); 01286 Exit(EXIT_FAILURE); 01287 } 01288 01289 if(verbose_general >= 3) IntfKeyPrintFields(IF); 01290 01291 nb_VoxX = nb_VoxX<0 ? IF.mtx_size[0] : nb_VoxX ; 01292 nb_VoxY = nb_VoxY<0 ? IF.mtx_size[1] : nb_VoxY ; 01293 nb_VoxZ = nb_VoxZ<0 ? IF.mtx_size[2] : nb_VoxZ ; 01294 01295 // TODO : nb_beds ? 01296 // If voxel sizes not provided/found, set them to the default value (1mm) 01297 vox_SizeX = IF.vox_size[0]<0 ? 1. : IF.vox_size[0]; 01298 vox_SizeY = IF.vox_size[1]<0 ? 1. : IF.vox_size[1]; 01299 vox_SizeZ = IF.vox_size[2]<0 ? 1. : IF.vox_size[2]; 01300 01301 01302 // Get the frame duration from the image in the case it is not provided by the user 01303 if (frame_list == "") 01304 { 01305 vector <string> image_filenames; 01306 if ( IntfIsMHD(path_to_initial_img, image_filenames) ) 01307 { 01308 Cerr("***** castor-proj() -> Error : while trying to read the interfile header '"<< path_to_initial_img << "' !" << endl); 01309 Exit(EXIT_FAILURE); 01310 } 01311 01312 if(IF.image_duration.size()>0) // Image durations have been provided 01313 { 01314 if(image_filenames.size() > 1) 01315 { 01316 // If image files are splitted over a serie of 3D image files, we will have image durations for each image, including gates 01317 // As we need to recover the duration for each gate once, we skip gates if the data contains any 01318 for(int fr=0 ; fr<IF.nb_time_frames ; fr++) 01319 { 01320 // Check nb of frames and parse framing 01321 if(IF.image_duration.size() != (uint32_t)IF.nb_time_frames * 01322 IF.nb_resp_gates * 01323 IF.nb_card_gates ) 01324 { 01325 Cerr("***** castor-proj() -> Interfile reading error of the input image :"<<endl); 01326 Cerr(" The number of provided image duration: '"<< IF.image_duration.size() 01327 << "' does not match the actual number of time frames * respiratory gates * cardiac gates ) '"<< IF.nb_time_frames * 01328 IF.nb_resp_gates * 01329 IF.nb_card_gates <<"' !" << endl); 01330 Exit(EXIT_FAILURE); 01331 } 01332 01333 // Recover image duration and dead frames (if any) in the frame list 01334 stringstream ss_id, ss_df; 01335 // Go to the next frame (skip gate if any) 01336 ss_id << IF.image_duration.at(fr * IF.nb_resp_gates* IF.nb_card_gates); 01337 string frame_duration_str = ss_id.str(); 01338 01339 // Check if there is any dead frame (pause between frame groups) 01340 if(IF.frame_group_pause.size()>0) 01341 { 01342 FLTNB dead_frame_sec = IF.frame_group_pause.at(fr * IF.nb_resp_gates* IF.nb_card_gates); 01343 if(dead_frame_sec > 0.) 01344 { 01345 ss_df << dead_frame_sec; 01346 frame_duration_str.append(",").append(ss_df.str()).append(";0"); 01347 } 01348 } 01349 01350 frame_list.append(frame_duration_str); 01351 01352 // add comma if not last frame 01353 if(fr+1 < IF.nb_time_frames) 01354 frame_list.append(","); 01355 } 01356 } 01357 else 01358 for(int fr=0 ; fr<IF.nb_time_frames ; fr++) 01359 { 01360 if(IF.image_duration.size() != IF.nb_time_frames) 01361 { 01362 Cerr("***** castor-proj() -> Interfile reading error of the input image :"<<endl); 01363 Cerr(" The number of provided image duration: '"<< IF.image_duration.size() 01364 << "' does not match the actual number of time frames * respiratory gates * cardiac gates ) '"<< IF.nb_time_frames <<"' !" << endl); 01365 Exit(EXIT_FAILURE); 01366 } 01367 01368 stringstream ss; 01369 ss << IF.image_duration.at(fr); 01370 frame_list.append(ss.str()); 01371 // add comma if not last frame 01372 if(fr+1 < IF.nb_time_frames) 01373 frame_list.append(","); 01374 } 01375 } 01376 } 01377 01378 // Check nb gating 01379 (IF.nb_resp_gates >1) ? nb_resp_gates = IF.nb_resp_gates : 1 ; 01380 (IF.nb_card_gates >1) ? nb_card_gates = IF.nb_card_gates : 1 ; 01381 01382 //TODO : offsets ? 01383 01384 // Default initialization of max ring difference with the max number of rings 01385 // TODO : How to deal with potential geometry containing several layers of crystals with different number of rings ? 01386 // Only for PET 01387 if(p_scannermanager->GetScannerType()==SCANNER_PET && max_ring_diff<0) max_ring_diff = p_scannermanager->GetScannerLayerNbRings(0); 01388 01389 if(nb_VoxX<0 || nb_VoxY<0 || nb_VoxZ<0) 01390 { 01391 ReadDataASCIIFile(p_scannermanager->GetPathToScannerFile(), "voxels number transaxial", &nb_VoxX, 1, KEYWORD_MANDATORY); 01392 ReadDataASCIIFile(p_scannermanager->GetPathToScannerFile(), "voxels number transaxial", &nb_VoxY, 1, KEYWORD_MANDATORY); 01393 ReadDataASCIIFile(p_scannermanager->GetPathToScannerFile(), "voxels number axial", &nb_VoxZ, 1, KEYWORD_MANDATORY); 01394 } 01395 01396 if ( (fov_SizeX<0 || fov_SizeY<0 || fov_SizeZ<0) && (vox_SizeX<0 || vox_SizeY<0 || vox_SizeZ<0) ) 01397 { 01398 ReadDataASCIIFile(p_scannermanager->GetPathToScannerFile(), "field of view transaxial", &fov_SizeX, 1, KEYWORD_MANDATORY); 01399 ReadDataASCIIFile(p_scannermanager->GetPathToScannerFile(), "field of view transaxial", &fov_SizeY, 1, KEYWORD_MANDATORY); 01400 ReadDataASCIIFile(p_scannermanager->GetPathToScannerFile(), "field of view axial", &fov_SizeZ, 1, KEYWORD_MANDATORY); 01401 } 01402 /* 01403 else if ( (fov_SizeX<0 || fov_SizeY<0 || fov_SizeZ<0) && (vox_SizeX>=0 || vox_SizeY>=0 || vox_SizeZ>=0)) 01404 { 01405 fov_SizeX = vox_SizeX*(FLTNB)nb_VoxX; 01406 fov_SizeY = vox_SizeY*(FLTNB)nb_VoxY; 01407 fov_SizeZ = vox_SizeZ*(FLTNB)nb_VoxZ; 01408 cout << vox_SizeX << " " << fov_SizeX << endl; 01409 } 01410 else 01411 { 01412 vox_SizeX = fov_SizeX/(FLTNB)nb_VoxX; 01413 vox_SizeY = fov_SizeY/(FLTNB)nb_VoxY; 01414 vox_SizeZ = fov_SizeZ/(FLTNB)nb_VoxZ; 01415 } 01416 */ 01417 if (verbose_general>=2) 01418 { 01419 Cout(" Image dimensions for analytical projections: " << endl); 01420 Cout(" voxels number (X,Y,Z): " << nb_VoxX << "," << nb_VoxY << "," << nb_VoxZ << endl); 01421 Cout(" voxels size (X,Y,Z): " << vox_SizeX << "," << vox_SizeY << "," << vox_SizeZ << endl); 01422 Cout(" FOV (X,Y,Z): " << fov_SizeX << "," << fov_SizeY << "," << fov_SizeZ << endl); 01423 } 01424 01425 // ============================================================================================================ 01426 // oImageDimensionsAndQuantification settings 01427 // ============================================================================================================ 01428 01429 // Create oImageDimensionsAndQuantification 01430 oImageDimensionsAndQuantification* p_ImageDimensionsAndQuantification = new oImageDimensionsAndQuantification(); 01431 p_ImageDimensionsAndQuantification->SetNbVoxX(nb_VoxX); 01432 p_ImageDimensionsAndQuantification->SetNbVoxY(nb_VoxY); 01433 p_ImageDimensionsAndQuantification->SetNbVoxZ(nb_VoxZ); 01434 p_ImageDimensionsAndQuantification->SetNbThreads(nb_threads); 01435 p_ImageDimensionsAndQuantification->SetNbBeds(nb_beds); 01436 p_ImageDimensionsAndQuantification->SetVoxSizeX(vox_SizeX); 01437 p_ImageDimensionsAndQuantification->SetVoxSizeY(vox_SizeY); 01438 p_ImageDimensionsAndQuantification->SetVoxSizeZ(vox_SizeZ); 01439 p_ImageDimensionsAndQuantification->SetFOVSizeX(fov_SizeX); 01440 p_ImageDimensionsAndQuantification->SetFOVSizeY(fov_SizeY); 01441 p_ImageDimensionsAndQuantification->SetFOVSizeZ(fov_SizeZ); 01442 p_ImageDimensionsAndQuantification->SetFOVOutMasking(0., 0); 01443 p_ImageDimensionsAndQuantification->SetOffsetX(offsetX); 01444 p_ImageDimensionsAndQuantification->SetOffsetY(offsetY); 01445 p_ImageDimensionsAndQuantification->SetOffsetZ(offsetZ); 01446 p_ImageDimensionsAndQuantification->SetVerbose(verbose_data); 01447 p_ImageDimensionsAndQuantification->SetNbRespGates(nb_resp_gates); 01448 p_ImageDimensionsAndQuantification->SetNbCardGates(nb_card_gates); 01449 p_ImageDimensionsAndQuantification->SetFrames(frame_list); 01450 if (p_ImageDimensionsAndQuantification->CheckParameters()) 01451 { 01452 Cerr("***** castor-proj() -> A problem occured while checking image dimensions parameters !" << endl); 01453 Exit(EXIT_FAILURE); 01454 } 01455 if (p_ImageDimensionsAndQuantification->Initialize()) 01456 { 01457 Cerr("***** castor-proj() -> A problem occured while initializing image dimensions !" << endl); 01458 Exit(EXIT_FAILURE); 01459 } 01460 01461 // Initialization of DynamicDataManager class, related 4D data splitting management 01462 if (p_ImageDimensionsAndQuantification->InitDynamicData( "", 0, 0, 0, nb_resp_gates, nb_card_gates ) ) 01463 { 01464 Cerr("***** castor-proj() -> A problem occured while initializing Dynamic data manager's class !" << endl); 01465 Exit(EXIT_FAILURE); 01466 } 01467 01468 01469 // ---------------------------------------------------------------------------------------- 01470 // Random Number Generator initialization: (we first required to know the number of threads to use from p_ImageDimensionsAndQuantification) 01471 // ---------------------------------------------------------------------------------------- 01472 01473 sRNG* p_RNG = sRNG::GetInstance(); 01474 p_RNG->SetVerbose(verbose_general); 01475 // Use a user-provided seed to initialize the RNG if one has been provided. Use random number otherwise. 01476 (seed_RNG>=0 )? 01477 p_RNG->Initialize(seed_RNG, p_ImageDimensionsAndQuantification->GetNbThreadsForProjection()): 01478 p_RNG->Initialize(p_ImageDimensionsAndQuantification->GetNbThreadsForProjection()); 01479 01480 01481 // ============================================================================================================ 01482 // Datafile initialization 01483 // ============================================================================================================ 01484 01485 vDataFile** p_DataFile = new vDataFile*[nb_beds]; 01486 01487 if (p_scannermanager->GetScannerType() < 0) 01488 { 01489 Cerr("***** castor-proj() -> Unknown scanner type !" << endl); 01490 Exit(EXIT_FAILURE); 01491 } 01492 else if (p_scannermanager->GetScannerType() == SCANNER_PET) 01493 { 01494 if (datafile_type!=TYPE_PET) 01495 { 01496 Cerr("***** castor-proj() -> Only PET type events are allowed while using a PET scanner !" << endl); 01497 Cerr("***** castor-proj() -> (use '-fileType' to define modality (0=PET, 1=SPECT, 2=CT) )" << endl); 01498 Exit(EXIT_FAILURE); 01499 } 01500 for (int i=0 ; i<nb_beds ; i++) 01501 { 01502 p_DataFile[i] = new iDataFilePET(); 01503 if(path_to_attenuation_img != "") // Enable correction factor for attenuation 01504 (dynamic_cast<iDataFilePET*>(p_DataFile[i]))->SetAtnCorrectionFlagOn(); 01505 } 01506 } 01507 else if (p_scannermanager->GetScannerType() == SCANNER_SPECT_CONVERGENT) 01508 { 01509 if (datafile_type!=TYPE_SPECT) 01510 { 01511 Cerr("***** castor-proj() -> Only SPECT type events are allowed while using a SPECT scanner !" << endl); 01512 Cerr("***** castor-proj() -> (use '-fileType' to define modality (0=PET, 1=SPECT, 2=CT) )" << endl); 01513 Exit(EXIT_FAILURE); 01514 } 01515 for (int i=0 ; i<nb_beds ; i++) 01516 { 01517 p_DataFile[i] = new iDataFileSPECT(); 01518 } 01519 } 01520 else 01521 { 01522 Cerr("***** castor-proj() -> Unsupported scanner type !" << endl); 01523 Exit(EXIT_FAILURE); 01524 } 01525 01526 // Load raw data in memory and do other stuff if needed. 01527 // TODO TM : Even if the loop is there, no support for multi-bed projection as for now (case we have to project an image as in a multi-bad acquisition) 01528 // It would first require some informations about the FOV of each bed (get this from interfile header of image or for command line options ?) 01529 // In the loop, we need to provide the bed in argument of PROJ_InitFile() in order to correctly initialize the beds 01530 for (int bed=0 ; bed<nb_beds ; bed++) 01531 { 01532 p_DataFile[bed]->SetBedIndex(bed); 01533 p_DataFile[bed]->SetPercentageLoad(data_file_percentage_load); 01534 p_DataFile[bed]->SetVerbose(verbose_data); 01535 p_DataFile[bed]->SetDataMode(datafile_mode); 01536 p_DataFile[bed]->SetDataType(datafile_type); 01537 p_DataFile[bed]->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification); 01538 p_DataFile[bed]->PROJ_GetScannerSpecificParameters(); 01539 p_DataFile[bed]->PROJ_InitFile(); 01540 p_DataFile[bed]->ComputeSizeEvent(); 01541 01542 // The Check parameters function is dedicated to reconstruction, as it checks things such as consistenties between size of the datafile and user-provided information 01543 if(p_DataFile[bed]->CheckParameters()) 01544 { 01545 Cerr("***** castor-proj() -> A problem occurred while checking datafile parameters ! Abort." << endl); 01546 Exit(EXIT_FAILURE); 01547 } 01548 p_DataFile[bed]->PrepareDataFile(); 01549 01550 /* 01551 if(p_DataFile[bed]->SetScannerSpecificParameters()) 01552 { 01553 Cerr("***** A problem occured while setting scanner parameters from the datafile !" << endl); 01554 Exit(EXIT_FAILURE); 01555 } 01556 */ 01557 01558 //p_DataFile[bed]->InitDynamicData(p_ImageDimensionsAndQuantification->GetNbTimeFrames(), nb_resp_gates, nb_card_gates, "", 0, 0, bed); 01559 01560 /* 01561 if (p_DataFile[bed]->InitDynamicData(nb_resp_gates, nb_card_gates, "", 0, 0, 0, p_DataFile[bed]) ) 01562 { 01563 Cerr("***** A problem occured while initializing Dynamic data !" << endl); 01564 Exit(EXIT_FAILURE); 01565 } 01566 */ 01567 } 01568 01569 01570 // ============================================================================================================ 01571 // Projector manager initialization 01572 // ============================================================================================================ 01573 01574 // Create object 01575 oProjectorManager* p_ProjectorManager = new oProjectorManager(); 01576 01577 // Set all parameters 01578 p_ProjectorManager->SetScanner(p_scannermanager->GetScannerObject()); 01579 p_ProjectorManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification); 01580 p_ProjectorManager->SetDataFile(p_DataFile[0]); 01581 p_ProjectorManager->SetComputationStrategy(projector_computation_strategy); 01582 p_ProjectorManager->SetOptionsForward(options_projector); 01583 p_ProjectorManager->SetOptionsBackward(options_projector); 01584 p_ProjectorManager->SetOptionsCommon(""); 01585 p_ProjectorManager->SetVerbose(verbose_proj); 01586 01587 // Check parameters 01588 if (p_ProjectorManager->CheckParameters()) 01589 { 01590 Cerr("***** castor-proj() -> A problem occured while checking projector manager's parameters !" << endl); 01591 Exit(EXIT_FAILURE); 01592 } 01593 01594 // Initialize projector manager 01595 if (p_ProjectorManager->Initialize()) 01596 { 01597 Cerr("***** castor-proj() -> A problem occured while initializing projector manager !" << endl); 01598 Exit(EXIT_FAILURE); 01599 } 01600 01601 01602 01603 01604 // ============================================================================================================ 01605 // Image Convolver manager initialization 01606 // ============================================================================================================ 01607 01608 // Create object 01609 oImageConvolverManager* p_ImageConvolverManager = new oImageConvolverManager(); 01610 // Set all parameters 01611 p_ImageConvolverManager->SetVerbose(verbose_conv); 01612 p_ImageConvolverManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification); 01613 p_ImageConvolverManager->SetOptions(options_image_convolver); 01614 01615 // Check parameters 01616 if (p_ImageConvolverManager->CheckParameters()) 01617 { 01618 Cerr("***** castor-proj() -> A problem occured while checking optimizer manager's parameters !" << endl); 01619 Exit(EXIT_FAILURE); 01620 } 01621 // Initialize optimizer manager 01622 if (p_ImageConvolverManager->Initialize()) 01623 { 01624 Cerr("***** castor-proj() -> A problem occured while initializing optimizer manager !" << endl); 01625 Exit(EXIT_FAILURE); 01626 } 01627 01628 01629 01630 // ============================================================================================================ 01631 // Image Space initialization 01632 // ============================================================================================================ 01633 01634 oImageSpace* p_ImageSpace = new oImageSpace(); 01635 p_ImageSpace->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification); 01636 p_ImageSpace->SetVerbose(verbose_data); 01637 01638 01639 01640 // ============================================================================================================ 01641 // Init alogorithm and proceed to analytical projection 01642 // ============================================================================================================ 01643 01644 oAnalyticProjection* p_Projection = new oAnalyticProjection(); 01645 p_Projection->InitOptimizer(p_ImageDimensionsAndQuantification); 01646 if (p_Projection->InitNoiseModel(noise_model) ) 01647 { 01648 Cerr("***** castor-proj() -> A problem occured while initializing noise model !" << endl); 01649 Exit(EXIT_FAILURE); 01650 } 01651 01652 p_Projection->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification); 01653 p_Projection->SetProjectorManager(p_ProjectorManager); 01654 p_Projection->SetImageConvolverManager(p_ImageConvolverManager); 01655 p_Projection->SetImageSpace(p_ImageSpace); 01656 p_Projection->SetDataFile(p_DataFile); 01657 p_Projection->SetGPUflag(gpu_flag); 01658 p_Projection->SetVerbose(verbose_algo); 01659 p_Projection->SetNbBeds(nb_beds); 01660 p_Projection->SetPathInitImage(path_to_initial_img); 01661 p_Projection->SetPathAtnImage(path_to_attenuation_img); 01662 p_Projection->SetScanner(p_scannermanager->GetScannerObject()); 01663 01664 // Launch analytical projection 01665 if (p_Projection->Launch()) 01666 { 01667 Cerr("***** castor-proj() -> Error occured during Projections" << endl;); 01668 } 01669 01670 01671 // ============================================================================================================ 01672 // End 01673 // ============================================================================================================ 01674 01675 if(SPECT_projection_angles) delete SPECT_projection_angles; 01676 01677 if (verbose_general>=5) Cout("----- Deleting CASToR objects ... -----" << endl); 01678 delete p_Projection; 01679 delete p_ImageSpace; 01680 delete p_ImageConvolverManager; 01681 delete p_ProjectorManager; 01682 01683 01684 // Delete objects 01685 for (int i=0 ; i<nb_beds ; i++) delete p_DataFile[i]; 01686 delete[] p_DataFile; 01687 delete p_ImageDimensionsAndQuantification; 01688 01689 01690 if (verbose_general>=5) Cout("----- CASToR objects successfully deleted -----" << endl); 01691 01692 // Ending 01693 if (verbose_general>=1) Cout(endl); 01694 #ifdef CASTOR_MPI 01695 MPI_Finalize(); 01696 #endif 01697 return EXIT_SUCCESS; 01698 }