CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
castor-recon.cc
Go to the documentation of this file.
1 /*
2 This file is part of CASToR.
3 
4  CASToR is free software: you can redistribute it and/or modify it under the
5  terms of the GNU General Public License as published by the Free Software
6  Foundation, either version 3 of the License, or (at your option) any later
7  version.
8 
9  CASToR is distributed in the hope that it will be useful, but WITHOUT ANY
10  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  details.
13 
14  You should have received a copy of the GNU General Public License along with
15  CASToR (in file GNU_GPL.TXT). If not, see <http://www.gnu.org/licenses/>.
16 
17 Copyright 2017-2018 all CASToR contributors listed below:
18 
19  --> current contributors: Thibaut MERLIN, Simon STUTE, Didier BENOIT, Claude COMTAT, Marina FILIPOVIC, Mael MILLARDET
20  --> past contributors: Valentin VIELZEUF
21 
22 This is CASToR version 2.0.
23 */
24 
36 #include "gVariables.hh"
37 #include "gOptions.hh"
38 #include "oIterativeAlgorithm.hh"
39 #include "oSensitivityGenerator.hh"
41 #include "iDataFilePET.hh"
42 #include "iDataFileSPECT.hh"
43 #include "iDataFileCT.hh"
44 #include "sOutputManager.hh"
45 #include "sScannerManager.hh"
47 #include "iScannerPET.hh"
48 #include "sAddonManager.hh"
49 #include "sChronoManager.hh"
50 
51 // =============================================================================================================================================
52 // =============================================================================================================================================
53 // =============================================================================================================================================
54 // H E L P F U N C T I O N S
55 // =============================================================================================================================================
56 // =============================================================================================================================================
57 // =============================================================================================================================================
58 
59 
64 void ShowHelp()
65 {
66  // Return when using MPI and mpi_rank is not 0
67  #ifdef CASTOR_MPI
68  int mpi_rank = 0;
69  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
70  if (mpi_rank!=0) return;
71  #endif
72  // Show help
73  cout << endl;
74  cout << "Usage: castor-recon.exe -df file.cdh -(f/d)out output -it iter [settings]" << endl;
75  cout << endl;
76  cout << "[Main options]:" << endl;
77  cout << " -df file.cdh : Give an input CASTOR datafile header (no default)." << endl;
78 //MULTIFILE cout << " -df file.cdh : Give an input CASTOR datafile header or several ones separated by commas (no default)." << endl;
79 //MULTIBED cout << " Can use this option multiple times to specify multiple bed positions." << endl;
80  cout << " -fout name : Give the root name for all output files (no default, alternative to -dout)" << endl;
81  cout << " -dout name : Give the name of the output directory where all output files will be written (no default, alternative to -fout)" << endl;
82  cout << " -it list : Give the sequence of iterations:subsets separated by commas (no default)." << endl;
83  cout << " -dim x,y,z : Give the number of voxels in each dimension (default: those of the scanner)." << endl;
84  cout << " -fov x,y,z : Give the size of the field-of-view in each dimension, in mm (default: those of the scanner, or calculated from" << endl;
85  cout << " the voxel sizes provided using '-vox')." << endl;
86  cout << " -vox x,y,z : Give the size of the voxels in each dimension, in mm (default: those of the scanner, or calculated from the fov" << endl;
87  cout << " if a fov value is provided using '-fov')." << endl;
88  cout << " -vb : Give the general verbosity level, from 0 (no verbose) to 5 and above (at the event level) (default: 1)." << endl;
89  cout << endl;
90  cout << "[Specific options]:" << endl;
91  cout << " -help-dim : Print out specific help about dimensions settings." << endl; // managed by main
92  cout << " -help-in : Print out specific help about input settings." << endl; // managed by main
93  cout << " -help-out : Print out specific help about output settings." << endl; // managed by main
94  cout << " -help-algo : Print out specific help about reconstruction optimizer algorithm settings." << endl; // managed by main
95  cout << " -help-proj : Print out specific help about projection operators." << endl; // managed by main
96  cout << " -help-dynamic : Print out specific help about dynamic methodologies settings." << endl; // managed by main
97  cout << " -help-imgp : Print out specific help about image processing modules." << endl; // managed by main
98  cout << " -help-comp : Print out specific help about computing settings." << endl; // managed by main
99  cout << " -help-corr : Print out specific help about all corrections that can be disabled." << endl;
100  cout << " -help-misc : Print out specific help about miscellaneous and verbose settings." << endl; // managed by main
101  cout << endl;
102  cout << "[Implemented Modules]:" << endl;
103  cout << " -help-scan : Show the list of all scanners from the configuration directory." << endl; // managed by sScannerManager
104  cout << " -help-projm : Show the list and description of all implemented projectors." << endl; // managed by oProjectorManager
105  cout << " -help-opti : Show the list and description of all implemented optimizer algorithms." << endl; // managed by oOptimizerManager
106  //cout << " -help-pnlt : Show the list and description of all implemented penalties for optimizers." << endl; // managed by oOptimizerManager
107  //cout << " -help-motion : Show the list and description of all implemented image-based deformation models." << endl; // managed by oImageDeformationManager
108  //cout << " -help-dynm : Show the list and description of all implemented dynamic models." << endl; // managed by oDynamicModelManager
109  cout << " -help-conv : Show the list and description of all implemented image convolvers." << endl; // managed by oImageConvolverManager
110  cout << " -help-proc : Show the list and description of all implemented image processing modules." << endl; // managed by oImageProcessingManager
111  cout << endl;
112  cout << endl;
113  cout << " --help,-h,-help : Print out this help page." << endl; // managed by main
114  cout << endl;
115  #ifdef CASTOR_MPI
116  cout << " Compiled with MPI" << endl;
117  #endif
118  #ifdef CASTOR_OMP
119  cout << " Compiled with OpenMP" << endl;
120  #endif
121  #ifdef CASTOR_GPU
122  cout << " Compiled for GPU" << endl;
123  #endif
124  #if defined(CASTOR_OMP) || defined(CASTOR_MPI) || defined(CASTOR_GPU)
125  cout << endl;
126  #endif
127  #ifdef BUILD_DATE
128  cout << " Build date: " << BUILD_DATE << endl;
129  cout << endl;
130  #endif
131  #ifdef CASTOR_VERSION
132  cout << " This program is part of the CASToR release version " << CASTOR_VERSION << "." << endl;
133  cout << endl;
134  #endif
135 }
136 
137 // =====================================================================
138 // ---------------------------------------------------------------------
139 // ---------------------------------------------------------------------
140 // =====================================================================
146 {
147  // Return when using MPI and mpi_rank is not 0
148  #ifdef CASTOR_MPI
149  int mpi_rank = 0;
150  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
151  if (mpi_rank!=0) return;
152  #endif
153  // Show help
154  cout << endl;
155  cout << "[Input settings]:" << endl;
156  cout << endl;
157 //MULTIBED cout << " -df-inv : Invert the order of provided datafiles corresponding to the different bed positions." << endl;
158 //MULTIBED cout << endl;
159 //MULTIBED cout << " -df file.Cdf : Give an input CASTOR datafile (no default). Can use this option multiple times to specify multiple bed positions." << endl;
160 //MULTIFILE cout << " -df file.Cdf : Give an input CASTOR datafile header, or several ones separated by commas (no default)." << endl;
161 //MULTIFILE cout << " Several datafiles can be provided only in histogram mode, when these multiple datafiles represent different time" << endl;
162 //MULTIFILE cout << " frames (identified by the time start and duration inside the header) or different respiratory or cardiac gates" << endl;
163 //MULTIFILE cout << " (identified by the respiratory or cardiac gate index provided in the header)." << endl;
164  cout << " -df file.Cdf : Give an input CASTOR datafile header (no default)." << endl;
165  cout << endl;
166  cout << " -img file.hdr : Give an input image as the initialization of the algorithm (default: uniform value)." << endl;
167  cout << endl;
168  cout << " -sens file.hdr : Provide the sensitivity image (default: sensitivity image is computed before reconstruction)." << endl;
169  cout << " The image file should integrate all sensitivity images if more than one are required. If dual-gating is enabled and if it" << endl;
170  cout << " requires sensitivity images for each gate, the image should integrate nb_resp_gates * nb_card_gates sensitivity images" << endl;
171  cout << " (all cardiac-gated based sensitivity images for each one of the respiratory gates)." << endl;
172  cout << endl;
173  cout << " -multimodal file.hdr : Provide additional images, from other modalities (anatomical, functional), or processed, for use in constrained recon " << endl;
174  cout << " To provide several additional images for a single reconstructed image, the option should be used several times, new option for new image file " << endl;
175  cout << endl;
176  cout << " -mask file.hdr : Provide a mask image. The mask image is used in the projection, which requires the application of this same mask to the sensitivity, to avoid border issues. " << endl;
177  cout << endl;
178  cout << " -norm file.cdh : For list-mode data, provide a normalization data file for sensitivity computation (default: use the scanner LUT and" << endl;
179  cout << " assume all LORs with a weight of 1.). This restricts the computation of the sensitivity to the LORs provided in the" << endl;
180  cout << " normalization file and associated normalization factors and/or attenuation factors." << endl;
181  cout << " For dynamic reconstructions with multiple frames, as many normalization files as frames can be supplied, their names" << endl;
182  cout << " separated by commas. This is useful when dead-times correction is included in the normalization factors." << endl;
183 //MULTIBED cout << " Can also use this option multiple times when multiple bed positions are reconstructed at once." << endl;
184  cout << endl;
185  cout << " -atn file.hdr : Give an input attenuation image (unit has to be cm-1) for sensitivity image generation or SPECT attenuation correction." << endl;
186  cout << endl;
187  cout << " -help-in : Print out this help." << endl;
188  cout << endl;
189 }
190 
191 // =====================================================================
192 // ---------------------------------------------------------------------
193 // ---------------------------------------------------------------------
194 // =====================================================================
200 {
201  // Return when using MPI and mpi_rank is not 0
202  #ifdef CASTOR_MPI
203  int mpi_rank = 0;
204  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
205  if (mpi_rank!=0) return;
206  #endif
207  // Show help
208  cout << endl;
209  cout << "[Output settings]:" << endl;
210  cout << endl;
211  cout << " -fout name : Give the root name for all output files. All output files will be written as 'name_suffix.ext'." << endl;
212  cout << " So the provided name should not end with '.' or '/' character. (no default, alternative to -dout)" << endl;
213  cout << " -dout name : Give the name of the output directory where all output files will be written. All files will also" << endl;
214  cout << " be prefixed by the name of the directory. The provided name should not end with '.' or '/' character." << endl;
215  cout << " (no default, alternative to -fout)" << endl;
216  cout << endl;
217  cout << " -oit list : Give the sequence of output iterations as a list of 'a:b' pairs separated by commas. This will output one" << endl;
218  cout << " iteration over 'a' until 'b' is reached, then it goes to the next pair of setting." << endl;
219  cout << " Set '-1' to save only the last iteration. (default: save all iterations)" << endl;
220  cout << endl;
221  cout << " -onbp prec : By default, numbers are displayed using scientific format. This option allows to customize the format and precision" << endl;
222  cout << " : The format is format,precision. f is the format (s=scientific, f=fixed), p is the precision" << endl;
223  cout << " eg: -onbp f,5 --> fixed with 5 digits precision, -onbp --> scientifici with max precision." << endl;
224  cout << endl;
225  cout << " -omd : (M)erge (D)ynamic images. Indicate if a dynamic serie of 3D images should be written on disk in one file" << endl;
226  cout << " instead of a serie of 3D images associated with an interfile metaheader" << endl;
227  cout << endl;
228 /* TO BE IMPLEMENTED
229  cout << " -otb : Flag to say that we want to save the time basis images too (default: only the frames are saved)." << endl;
230  cout << endl;
231 */
232  cout << " -osens : Flag to say that we want to save the sensitivity image of each subset/iteration, when in histogram mode." << endl;
233  cout << endl;
234  cout << " -osub : Flag to say that we want to save the image after each subset." << endl;
235  cout << endl;
236  cout << " -olut : If a scanner LUT (geometry information) is computed from a .geom file, it will be save on disk in the scanner repository" << endl;
237  cout << endl;
238  cout << " -sens-histo : If input file is a histogram, compute the global sensitivity from it (and still proceed to normal reconstruction after)" << endl;
239  cout << endl;
240  cout << " -sens-only : For list-mode data, exit directly after computing and saving the sensitivity." << endl;
241  cout << endl;
242  cout << " -help-out : Print out this help." << endl;
243  cout << endl;
244 }
245 
246 // =====================================================================
247 // ---------------------------------------------------------------------
248 // ---------------------------------------------------------------------
249 // =====================================================================
255 {
256  // Return when using MPI and mpi_rank is not 0
257  #ifdef CASTOR_MPI
258  int mpi_rank = 0;
259  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
260  if (mpi_rank!=0) return;
261  #endif
262  // Show help
263  cout << endl;
264  cout << "[Dimensions options]:" << endl;
265  cout << endl;
266  cout << " -dim x,y,z : Give the number of voxels in each dimension (default: those of the scanner)." << endl;
267  cout << endl;
268  cout << " -fov x,y,z : Give the size of the field-of-view in each dimension, in mm (default: those of the scanner)." << endl;
269  cout << endl;
270  cout << " -vox x,y,z : Give the size of the voxels in each dimension, in mm (default: those of the scanner, or calculated from the fov if a fov value is provided using '-fov')." << endl;
271  cout << endl;
272  cout << " -off x,y,z : Give the offset of the field-of-view in each dimension, in mm (default: 0.,0.,0.)." << endl;
273  cout << " (note this has no effect when using a pre-computed system matrix)" << endl;
274  cout << endl;
275  cout << " -fov-out percent : Give the percentage of the eliptical transaxial FOV to be kept while saving the image (default: no making)" << endl;
276  cout << endl;
277  cout << " -slice-out value : Give the number of axial slices to be masked at each border of the axial field-of-view (default: 0)" << endl;
278  cout << endl;
279  cout << " -flip-out value : Flip the image before saving it (not done in the computation); specify the axis (e.g. 'X', 'XY', 'YZ') (default: no flip)" << endl;
280  cout << endl;
281  cout << " -help-dim : Print out this help." << endl;
282  cout << endl;
283 }
284 
285 // =====================================================================
286 // ---------------------------------------------------------------------
287 // ---------------------------------------------------------------------
288 // =====================================================================
294 {
295  // Return when using MPI and mpi_rank is not 0
296  #ifdef CASTOR_MPI
297  int mpi_rank = 0;
298  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
299  if (mpi_rank!=0) return;
300  #endif
301  // Show help
302  cout << "[Algorithm settings]:" << endl;
303  cout << endl;
304  cout << " -opti param : Give the algorithm to be used, along with a configuration file (algo:file.conf) or the list of parameters" << endl;
305  cout << " associated to the algorithm (algo,param1,param2,...). If the algorithm only is specified, the default" << endl;
306  cout << " configuration file is used if any. By default, the MLEM algorithm is used. For specific help, use option -help-opti." << endl;
307  cout << endl;
308  cout << " -opti-fom : Flag to say that we want to compute and print figures-of-merit (likelihood and RMSE) in the data-space." << endl;
309  cout << endl;
310  cout << " -opti-stat : Flag to say that we want to compute and print basic statistics about the image update." << endl;
311  cout << endl;
312 /* NOT YET IMPLEMENTED
313  cout << " -pnlt param : Give the penalty to be used with the algorithm (if the latter allows to do so), along with a configuration file (penalty:file.conf)" << endl;
314  cout << " or the list of parameters associated to the penalty (penalty,param1,param2,...). If the penalty only is specified, the default" << endl;
315  cout << " configuration file is used if any. By default, no penalty is used. For specific help, use option -help-penalty." << endl;
316  cout << endl;
317 */
318  cout << " -help-opti : Print out specific help about algorithm settings." << endl; // managed by oOptimizerManager
319  cout << endl;
320 /* NOT YET IMPLEMENTED
321  cout << " -help-pnlt : Print out specific help about penalty settings." << endl; // managed by oOptimizerManager
322  cout << endl;
323 */
324 }
325 
326 // =====================================================================
327 // ---------------------------------------------------------------------
328 // ---------------------------------------------------------------------
329 // =====================================================================
335 {
336  // Return when using MPI and mpi_rank is not 0
337  #ifdef CASTOR_MPI
338  int mpi_rank = 0;
339  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
340  if (mpi_rank!=0) return;
341  #endif
342  // Show help
343  cout << "[Projector settings]:" << endl;
344  cout << endl;
345  cout << " -proj param : Give the projector to be used for both forward and backward projections, along with a configuration file (proj:file.conf)" << endl;
346  cout << " or the list of parameters associated to the projector (proj,param1,param2,...). If the projector only is specified, the" << endl;
347  cout << " default configuration file is used. By default, the Siddon projector is used. For specific help, use option -help-proj." << endl;
348  cout << endl;
349  cout << " -projF param : Give the projector to be used for forward projections. See option -proj for details." << endl;
350  cout << endl;
351  cout << " -projB param : Give the projector to be used for backward projections. See option -proj for details." << endl;
352  cout << endl;
353  cout << " -proj-common : Give parameterization of different common projector-related options." << endl;
354  cout << endl;
355 /* TO BE IMPLEMENTED
356  cout << " -ignore-POI : Flag to say that we want to ignore any potential POI information (default: use it if present in the datafile)." << endl;
357  cout << endl;
358  * */
359  cout << " -ignore-TOF : Flag to say that we want to ignore any potential TOF information (default: use it if present in the datafile)." << endl;
360  cout << endl;
361  cout << " -help-projm : Print out specific help about projector settings." << endl; // managed by oProjectorManager
362  cout << endl;
363 }
364 
365 // =====================================================================
366 // ---------------------------------------------------------------------
367 // ---------------------------------------------------------------------
368 // =====================================================================
374 {
375  // Return when using MPI and mpi_rank is not 0
376  #ifdef CASTOR_MPI
377  int mpi_rank = 0;
378  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
379  if (mpi_rank!=0) return;
380  #endif
381  // Show help
382  cout << "[Image processing settings]:" << endl;
383  cout << endl;
384  cout << " -conv param;when : Give an image convolver model to be used within the algorithm, along with a configuration file (conv:file.conf::when) or the" << endl;
385  cout << " list of parameters associated to the convolver (conv,param1,param2,...::when). If the convolver only is specified, its default" << endl;
386  cout << " configuration file is used. By default, no convolver is applied. Multiple convolvers can be combined simply by repeating this" << endl;
387  cout << " this option. The mandatory 'when' parameter specifies when the convolver is applied (psf, sieve, forward, backward, post, intra)." << endl;
388  cout << " For more specific help, use option -help-conv." << endl;
389  cout << endl;
390  cout << " -help-conv : Print out specific help about the image convolver settings." << endl; // managed by oImageConvolverManager
391  cout << endl;
392  cout << " -proc param;when : Give an image processing module to be used within the algorithm, along with a configuration file (proc:file.conf;when) or the" << endl;
393  cout << " list of parameters associated to the module (proc,param1,param2,...;when). If the module only is specified, its default" << endl;
394  cout << " configuration file is used. By default, no image processing module is applied. Multiple modules can be combined simply by" << endl;
395  cout << " repeating this this option. The mandatory 'when' parameter specifies when the module is applied (forward, backward, post, intra)." << endl;
396  cout << " For more specific help, use option -help-proc." << endl;
397  cout << endl;
398  cout << " -help-proc : Print out specific help about the image processing module settings." << endl; // managed by oImageProcessingManager
399  cout << endl;
400 }
401 
402 // =====================================================================
403 // ---------------------------------------------------------------------
404 // ---------------------------------------------------------------------
405 // =====================================================================
411 {
412  // Return when using MPI and mpi_rank is not 0
413  #ifdef CASTOR_MPI
414  int mpi_rank = 0;
415  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
416  if (mpi_rank!=0) return;
417  #endif
418  // Show help
419  cout << endl;
420  cout << "[Dynamic settings]:" << endl;
421  cout << endl;
422  cout << " -frm list : Give the framing for the reconstruction where 'list' is a list of all frames durations in seconds separated" << endl;
423  cout << " by commas. To specify a dead frame, add 'df' concatenated to the frame duration. Milliseconds maximum precision." << endl;
424  cout << " (default: 1 frame of the whole input file duration)" << endl;
425  cout << endl;
426  cout << " -g path_to_file : Provide text file for dynamic data splitting associated to respiratory/cardiac gating or involuntary patient motion correction" << endl;
427  cout << endl;
428 /* TO BE VALIDATED
429  cout << " -time-basis nb:file : Give the number of time basis functions related to the frames given by option -frm." << endl;
430  cout << " The file should contain one line for each time-basis function, and as many numbers as frames." << endl;
431  cout << " -resp-basis nb:file : Give the number of intrinsic respiratory basis functions related to the respiratory gates" << endl;
432  cout << " The file should contain one line for each respiratory basis function, and as many numbers as respiratory gates." << endl;
433  cout << " -card-basis nb:file : Give the number of intrinsic cardiac basis functions related to the cardiac gates" << endl;
434  cout << " The file should contain one line for each cardiac basis function, and as many numbers as cardiac gates." << endl;
435  cout << endl;
436  cout << " -dynamic-model param : Dynamic model applied to either the frames of a dynamic acquisition, respiratory-gated frames, cardiac-gated frames, or simultaneously between these datasets." << endl;
437  cout << " Select the dynamic model to be used, along with a configuration file (model:file) or the list of parameters associated to the model (model_name,param1,param2,...)." << endl;
438  cout << endl;
439  cout << " -rm param : Provide an image-based deformation model to be used for respiratory motion correction, along with a configuration file (deformation:file)" << endl;
440  cout << " or the list of parameters associated to the projector (deformation,param1,param2,...)." << endl;
441  cout << endl;
442  cout << " -cm param : Provide an image-based deformation model to be used for cardiac motion correction, along with a configuration file (deformation:file)" << endl;
443  cout << " or the list of parameters associated to the projector (deformation,param1,param2,...)." << endl;
444  cout << endl;
445  cout << " -rcm param : Provide an image-based deformation model to be used for both respiratory and cardiac motion corrections, along with a configuration file (deformation:file)" << endl;
446  cout << " or the list of parameters associated to the projector (deformation,param1,param2,...)." << endl;
447  cout << endl;
448  cout << " -im param : Provide an image-based deformation model to be used for involuntary motion correction, along with a configuration file (deformation:file)" << endl;
449  cout << " or the list of parameters associated to the projector (deformation,param1,param2,...)." << endl;
450  cout << endl;
451  cout << " -qdyn file : Provide a text file containing quantitative factors specific to dynamic frames, respiratory or cardiac gates." << endl;
452  cout << " The file should provide factors with the keywords 'FRAME_QUANTIFICATION_FACTORS' and 'GATE_QUANTIFICATION_FACTORS' " << endl;
453  cout << " The number of factors must be consistent with the number of frames/gates " << endl;
454  cout << " If the data contains several frames and gates, the gate quantification factors should be entered on a separate line for each frame " << endl;
455  cout << endl;
456  cout << " -help-dynm : Print out specific help about dynamic model." << endl; // managed by oDynamicModelManager
457  cout << endl;
458  cout << " -help-motion : Print out specific help about deformation in image space for motion correction." << endl; // managed by oImageDeformationManager
459  cout << endl;
460  */
461  cout << " -help-dynamic : Print out this help." << endl;
462  cout << endl;
463 }
464 
465 // =====================================================================
466 // ---------------------------------------------------------------------
467 // ---------------------------------------------------------------------
468 // =====================================================================
474 {
475  // Return when using MPI and mpi_rank is not 0
476  #ifdef CASTOR_MPI
477  int mpi_rank = 0;
478  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
479  if (mpi_rank!=0) return;
480  #endif
481  // Show help
482  cout << endl;
483  cout << "[Computation settings]:" << endl;
484  cout << endl;
485  #ifdef CASTOR_OMP
486  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;
487  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;
488  cout << endl;
489  #endif
490  #ifdef CASTOR_GPU
491  cout << " -gpu : Flag to say that we want to use the GPU device (default: use the CPU only)." << endl;
492  cout << endl;
493  #endif
494  cout << " -load size : Give a number between 0 and 100 to specify the number of events which will be loaded in the RAM." << endl;
495  cout << " It will be interpreted as a percent ratio (default: 100 means all data buffered. 0 means no data buffered.)" << endl;
496  cout << endl;
497  cout << " -proj-comp : Give the strategy for projection line computation. Here are the three different strategies that can be used:" << endl;
498  cout << " 1 : Image-based system matrix elements storage: The voxels weights are added in a matrix representing the whole image, so" << endl;
499  cout << " the addition of a new line to the previous ones is straightforward only by adding the weights to the corresponding voxels." << endl;
500  cout << " As it is insanely long, it can possibly be used for example with extremely complex projectors that makes use of huge number" << endl;
501  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;
502  cout << " This strategy is not compatible with SPECT reconstruction including attenuation correction." << endl;
503  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;
504  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;
505  cout << " has a fixed size which is provided by the EstimateMaxNumberOfVoxelsPerLine() function from the vProjector class. There are no" << endl;
506  cout << " ckecks at all for possible buffer overflows. This is the fastest strategy and default one." << endl;
507  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;
508  cout << " upgraded if the current number of contributing voxels exceed the list's size. The first allocated size corresponds to the diagonal" << endl;
509  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;
510  cout << " is a bit slower than the fixed-list strategy during the first iteration, but is optimal with respect to RAM usage." << endl;
511  cout << endl;
512  cout << " -rng-seed : Give a seed for the random number generator (should be >=0)" << endl;
513  cout << endl;
514  cout << " -rng-extra : Give the number of additional random number generators, if needed (in addition to one rng per thread) " << endl;
515  cout << endl;
516  cout << " -help-comp : Print out this help." << endl;
517  cout << endl;
518  #ifdef CASTOR_MPI
519  cout << " Compiled with MPI" << endl;
520  #endif
521  #ifdef CASTOR_OMP
522  cout << " Compiled with OpenMP" << endl;
523  #endif
524  #ifdef CASTOR_GPU
525  cout << " Compiled for GPU" << endl;
526  #endif
527  #if defined(CASTOR_OMP) || defined(CASTOR_MPI) || defined(CASTOR_GPU)
528  cout << endl;
529  #endif
530 }
531 
532 // =====================================================================
533 // ---------------------------------------------------------------------
534 // ---------------------------------------------------------------------
535 // =====================================================================
541 {
542  // Return when using MPI and mpi_rank is not 0
543  #ifdef CASTOR_MPI
544  int mpi_rank = 0;
545  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
546  if (mpi_rank!=0) return;
547  #endif
548  // Show help
549  cout << endl;
550  cout << "[Miscellaneous settings]:" << endl;
551  cout << endl;
552  cout << " -vb : Give the general verbosity level, from 0 (no verbose) to 5 and above (at the event level) (default: 1)." << endl;
553  cout << " -vb-algo : Give the verbose level specific to the algorithm (default: same as general verbose level)." << endl;
554  cout << " -vb-opti : Give the verbose level specific to the optimizer (default: same as general verbose level)." << endl;
555  cout << " -vb-proj : Give the verbose level specific to the projector (default: same as general verbose level)." << endl;
556  cout << " -vb-conv : Give the verbose level specific to the image convolver (default: same as general verbose level)." << endl;
557  cout << " -vb-proc : Give the verbose level specific to the image processing (default: same as general verbose level)." << endl;
558  cout << " -vb-scan : Give the verbose level specific to the scanner (default: same as general verbose level)." << endl;
559  cout << " -vb-data : Give the verbose level specific to the data and image management (default: same as general verbose level)." << endl;
560  cout << " -vb-defo : Give the verbose level specific to the deformation (default: same as general verbose level)." << endl;
561  cout << " -vb-dyna : Give the verbose level specific to the dynamic model (default: same as general verbose level)." << endl;
562  cout << " -vb-sens : Give the verbose level specific to the sensitivity computation (default: same as general verbose level)." << endl;
563  cout << endl;
564  cout << " -conf : Give the path to the CASToR config directory (default: located through the CASTOR_CONFIG environment variable)." << endl;
565  cout << endl;
566  cout << " -help-misc : Print out this help." << endl;
567  cout << endl;
568 }
569 
570 // =====================================================================
571 // ---------------------------------------------------------------------
572 // ---------------------------------------------------------------------
573 // =====================================================================
579 {
580  // Return when using MPI and mpi_rank is not 0
581  #ifdef CASTOR_MPI
582  int mpi_rank = 0;
583  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
584  if (mpi_rank!=0) return;
585  #endif
586  // Show help
587  cout << endl;
588  cout << "[Correction settings]:" << endl;
589  cout << endl;
590  cout << " -ignore-corr list : Give the list of corrections that should be ignored, separated by commas (default: all corrections applied if present)." << endl;
591  cout << " Here is a list of all corrections that can be ignored:" << endl;
592  cout << " - attn: to ignore the attenuation correction (emission only)" << endl;
593  cout << " - norm: to ignore the normalization correction (emission only)" << endl;
594  cout << " - rand: to ignore the random correction (PET only)" << endl;
595  cout << " - scat: to ignore the scatter correction" << endl;
596  cout << " - deca: to ignore the decay correction (emission only)" << endl;
597  cout << " - brat: to ignore the branching ratio correction (emission only)" << endl;
598  cout << " - fdur: to ignore the frame duration correction (emission only)" << endl;
599  cout << " - cali: to ignore the calibration correction (emission only)" << endl;
600  cout << endl;
601 }
602 
603 // =============================================================================================================================================
604 // =============================================================================================================================================
605 // =============================================================================================================================================
606 // M A I N P R O G R A M
607 // =============================================================================================================================================
608 // =============================================================================================================================================
609 // =============================================================================================================================================
610 
611 int main(int argc, char** argv)
612 {
613  // ============================================================================================================
614  // MPI stuff
615  // ============================================================================================================
616  int mpi_rank = 0;
617  int mpi_size = 1;
618  #ifdef CASTOR_MPI
619  MPI_Init(&argc, &argv);
620  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
621  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
622  #endif
623 
624  // No argument, then show help
625  if (argc==1)
626  {
627  ShowHelp();
628  Exit(EXIT_SUCCESS);
629  }
630 
631  // ============================================================================================================
632  // Parameterized variables with their default values
633  // ============================================================================================================
634 
635  // --------------------------------------------------------------------------------
636  // Dimensions settings
637  // --------------------------------------------------------------------------------
638 
639  // Initialization of the voxel and field-of-view values in each spatial dimensions. default values depend of the scanner.
640  INTNB nb_voxX=-1, nb_voxY=-1, nb_voxZ=-1;
641  FLTNB fov_sizeX=-1., fov_sizeY=-1., fov_sizeZ=-1.;
642  FLTNB vox_sizeX=-1., vox_sizeY=-1., vox_sizeZ=-1.;
643  // Initialization offset values of the field-of-view in each spatial dimensions
644  FLTNB offsetX = 0., offsetY = 0., offsetZ = 0.;
645  FLTNB fov_out = 0.;
646  INTNB slice_out = 0;
647  // Output flip
648  string flip_out = "";
649 
650  // --------------------------------------------------------------------------------
651  // Dynamic settings
652  // --------------------------------------------------------------------------------
653 
654  // Frame list descriptor
655  string frame_list = "";
656  // Number of respiratory gates in the data. default : 1
657  int nb_resp_gates = 1;
658  // Number of cardiac gates in the data. default : 1
659  int nb_card_gates = 1;
660  // Path to file containing the respiratory/cardiac gate splitting of the data
661  string path_to_4D_data_splitting_file = "";
662  // String gathering the name of the used dynamic model applied to the frames/resp gates/card gates of the dynamic acquisition (default : none), corresponding file, and corresponding string gathering parameters
663  string dynamic_model_options = "";
664  // String gathering the name of the used deformation model for respiratory motion (default : none), corresponding file, and corresponding string gathering parameters
665  string resp_motion_options = "";
666  // String gathering the name of the used deformation model for cardiac motion (default : none), corresponding file, and corresponding string gathering parameters
667  string card_motion_options = "";
668  // String gathering the name of the used deformation model for both respiratory and cardiac motions (default : none), corresponding file, and corresponding string gathering parameters
669  string double_motion_options = "";
670  // String gathering the name of the used deformation model for involuntary motion (default : none), corresponding file, and corresponding string gathering parameters
671  string ipat_motion_options = "";
672  // Number of time basis functions and associated file
673  int nb_time_basis = 1;
674  string path_to_time_basis_coef = "";
675  // Number of respiratory basis functions and associated file
676  int nb_resp_basis = 1;
677  string path_to_resp_basis_coef = "";
678  // Number of cardiac basis functions and associated file
679  int nb_card_basis = 1;
680  string path_to_card_basis_coef = "";
681  // Number of cardiac basis functions and associated file
682  string path_to_dynamic_quantification_file= "";
683 
684  // --------------------------------------------------------------------------------
685  // Input settings
686  // --------------------------------------------------------------------------------
687 
688  // Number of bed positions
689  int nb_beds = 0;
690  // Vector (for bed positions) containing string which gather the path to the data filenames provided by the user. no default
691  vector<string> path_to_data_filename;
692  // Invert bed datafile names order (bed positions)
693  bool invert_datafile_bed_order_flag = false;
694  // Path to an image used as initialization. no default (uniform image is used in this case)
695  string path_to_initial_img = "";
696  // Path to an image used as initialization for the sensitivity image. no default (sensitivity image is computed in this case)
697  string path_to_sensitivity_img;
698  // Path to a normalization data file for sensitivity computation with list-mode data
699  vector<string> path_to_normalization_filename;
700  // Path to an image used for the attenuation. default : uniform image
701  string path_to_attenuation_img;
702  // Paths to additional images
703  vector<string> path_to_multimodal_img;
704  // Path to an image used as mask
705  string path_to_mask_img = "";
706  // Number of resp and card atn images for sensitivity image generation
707  int nb_atn_resp_imgs = 1;
708  int nb_atn_card_imgs = 1;
709  // Ignored corrections (by default, if present, there are applied)
710  string ignored_corrections = "";
711 
712  // --------------------------------------------------------------------------------
713  // Output settings
714  // --------------------------------------------------------------------------------
715 
716  // Output directory name.
717  string path_dout = "";
718  // Or root name
719  string path_fout = "";
720  // Iterations to be saved
721  string output_iterations = "";
722  // Merge output dynamic images on disk into one file or not
723  bool merge_dynamic_imgs_flag = false;
724  // Write scanner LUT generated by a geom file on disk
725  bool save_LUT_flag = false;
726  // Save the sensitivity image in histogram mode for each subset/iteration
727  bool save_sens_histo = false;
728  // Save the image after each subset
729  bool save_subset_image = false;
730  // Save the time basis functions (flag). default : no saving
731  //bool save_time_basis_functions_flag = false; //TODO
732  // Exit directly after computing the sensitivity
733  bool exit_after_sensitivity = false;
734  // Compute sensitivity with a histogram file as input
735  bool sensitivity_from_histogram = false;
736  // Precision for output number display
737  string onb_prec = "s,0";
738 
739  // --------------------------------------------------------------------------------
740  // Projector settings
741  // --------------------------------------------------------------------------------
742 
743  // String gathering the name of the projector for forward/backward projection with specific options (default: Native Siddon for both)
744  string options_projector = "incrementalSiddon";
745  string options_projectorF = "incrementalSiddon";
746  string options_projectorB = "incrementalSiddon";
747  string options_projector_common = "";
748  // Default projector computation strategy
749  int projector_computation_strategy = FIXED_LIST_COMPUTATION_STRATEGY;
750  // POI & TOF flags for the use of such information in the reconstruction
751  bool ignore_POI = false;
752  bool ignore_TOF = false;
753 
754  // --------------------------------------------------------------------------------
755  // Algorithm settings
756  // --------------------------------------------------------------------------------
757 
758  // Numbers of iterations and associated numbers of subsets (default : 1 for both)
759  string nb_iterations_subsets = "";
760  // String gathering the name of the optimizer with specific options (default: MLEM)
761  string options_optimizer = "MLEM";
762  // Boolean to say if we want to compute FOM in the data-space
763  bool optimizer_fom = false;
764  // Boolean to say if we want to compute image update basic statistics
765  bool optimizer_stat = false;
766  // String gathering the name of the penalty with specific options (default: no penalty)
767  string options_penalty = "";
768 
769  // --------------------------------------------------------------------------------
770  // Image convolvers and processing modules
771  // --------------------------------------------------------------------------------
772 
773  // String vector gathering the name of the image convolvers with specific options (default: no image convolver)
774  vector<string> options_image_convolver;
775  // String vector gathering the name of the image processing modules with specific options (default: no image processing module)
776  vector<string> options_image_processing;
777 
778  // --------------------------------------------------------------------------------
779  // Computation settings
780  // --------------------------------------------------------------------------------
781 
782  // Using GPU (flag) ->NOTE : default : only CPU
783  bool gpu_flag = false;
784  // Number of threads
785  string nb_threads = "1";
786 
787  // --------------------------------------------------------------------------------
788  // Miscellaneous settings
789  // --------------------------------------------------------------------------------
790 
791  // General verbose level
792  int verbose_general = 1;
793  // Specific verbose levels
794  int verbose_algo = -1;
795  int verbose_opti = -1;
796  int verbose_proj = -1;
797  int verbose_conv = -1;
798  int verbose_proc = -1;
799  int verbose_scan = -1;
800  int verbose_data = -1;
801  int verbose_defo = -1;
802  int verbose_dyna = -1;
803  int verbose_sens = -1;
804  // Path to config directory
805  string path_to_config_dir = "";
806  // Initial seed for random number generator
807  int64_t random_generator_seed = -1;
808  int nb_extra_random_generators = 0;
809 
810  // ============================================================================================================
811  // Read command-line parameters
812  // ============================================================================================================
813 
814  // TODO : replace remaining atoi, atof, etc, by ConvertFromString()
815 
816  // Must manually increment the option index when an argument is needed after an option
817  for (int i=1; i<argc; i++)
818  {
819  // Get the option as a string
820  string option = (string)argv[i];
821 
822  // --------------------------------------------------------------------------------
823  // Miscellaneous settings
824  // --------------------------------------------------------------------------------
825 
826  // Show help
827  if (option=="-h" || option=="--help" || option=="-help")
828  {
829  ShowHelp();
830  Exit(EXIT_SUCCESS);
831  }
832  // Specific help for integrated scanners (managed by scanner manager)
833  else if (option=="-help-scan")
834  {
835  if(sScannerManager::GetInstance()->ShowScannersDescription())
836  Cerr("***** castor-recon() -> An error occured when trying to output the available scanners from the scanner repository !'" << endl;);
837  Exit(EXIT_SUCCESS);
838  }
839  // Specific help for optimizer settings (managed by optimizer children)
840  else if (option=="-help-opti")
841  {
843  Exit(EXIT_SUCCESS);
844  }
845  // Specific help for penalty settings (managed by penalty children)
846  else if (option=="-help-pnlt")
847  {
849  Exit(EXIT_SUCCESS);
850  }
851  // Specific help for projector settings (managed by projector children)
852  else if (option=="-help-projm")
853  {
854  // Call specific showHelp function from vProjector children
856  // Call the static showHelp function from vProjector
858  Exit(EXIT_SUCCESS);
859  }
860  // Specific help for image convolver settings (managed by image convolver children)
861  else if (option=="-help-conv")
862  {
863  // Call specific showHelp function from vImageConvolver children
865  // Call the static showHelp function from oImageConvolverManager
867  Exit(EXIT_SUCCESS);
868  }
869  // Specific help for image processing settings (managed by image processing children)
870  else if (option=="-help-proc")
871  {
872  // Call specific showHelp function from vImageProcessingModule children
874  // Call the static showHelp function from oImageProcessingManager
876  Exit(EXIT_SUCCESS);
877  }
878  // Specific help for dynamic model (managed by dynamic model children)
879  else if (option=="-help-dynm")
880  {
882  Exit(EXIT_SUCCESS);
883  }
884  // Specific help for image deformation settings (managed by projector children)
885  else if (option=="-help-motion")
886  {
888  Exit(EXIT_SUCCESS);
889  }
890  // Specific help for input settings (managed by main)
891  else if (option=="-help-in")
892  {
893  ShowHelpInput();
894  Exit(EXIT_SUCCESS);
895  }
896  // Specific help for output settings (managed by main)
897  else if (option=="-help-out")
898  {
899  ShowHelpOutput();
900  Exit(EXIT_SUCCESS);
901  }
902  // Specific help for dimensions settings (managed by main)
903  else if (option=="-help-dim")
904  {
906  Exit(EXIT_SUCCESS);
907  }
908  // Specific help for optimizer settings (managed by main)
909  else if (option=="-help-algo")
910  {
911  ShowHelpAlgo();
912  Exit(EXIT_SUCCESS);
913  }
914  // Specific help for projector settings (managed by main)
915  else if (option=="-help-proj")
916  {
917  ShowHelpProj();
918  Exit(EXIT_SUCCESS);
919  }
920  // Specific help for image processing settings (managed by main)
921  else if (option=="-help-imgp")
922  {
923  ShowHelpImgp();
924  Exit(EXIT_SUCCESS);
925  }
926  // Specific help for computation settings (managed by main)
927  else if (option=="-help-comp")
928  {
930  Exit(EXIT_SUCCESS);
931  }
932  // Specific help for miscellaneous settings (managed by main)
933  else if (option=="-help-misc")
934  {
936  Exit(EXIT_SUCCESS);
937  }
938  // Specific help for dynamic settings (managed by main)
939  else if (option=="-help-dynamic")
940  {
941  ShowHelpDynamic();
942  Exit(EXIT_SUCCESS);
943  }
944  // Specific help for correction settings (managed by main)
945  else if (option=="-help-corr")
946  {
948  Exit(EXIT_SUCCESS);
949  }
950  // General verbosity level
951  else if (option=="-vb")
952  {
953  if (i>=argc-1)
954  {
955  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
956  Exit(EXIT_FAILURE);
957  }
958  if (ConvertFromString(argv[i+1], &verbose_general))
959  {
960  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_general << " for option: " << option << endl);
961  Exit(EXIT_FAILURE);
962  }
963  i++;
964  }
965  // Algorithm verbosity level
966  else if (option=="-vb-algo")
967  {
968  if (i>=argc-1)
969  {
970  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
971  Exit(EXIT_FAILURE);
972  }
973  if (ConvertFromString(argv[i+1], &verbose_algo))
974  {
975  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_algo << " for option: " << option << endl);
976  Exit(EXIT_FAILURE);
977  }
978  i++;
979  }
980  // Optimizer verbosity level
981  else if (option=="-vb-opti")
982  {
983  if (i>=argc-1)
984  {
985  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
986  Exit(EXIT_FAILURE);
987  }
988  if (ConvertFromString(argv[i+1], &verbose_opti))
989  {
990  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_opti << " for option: " << option << endl);
991  Exit(EXIT_FAILURE);
992  }
993  i++;
994  }
995  // Projector verbosity level
996  else if (option=="-vb-proj")
997  {
998  if (i>=argc-1)
999  {
1000  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1001  Exit(EXIT_FAILURE);
1002  }
1003  if (ConvertFromString(argv[i+1], &verbose_proj))
1004  {
1005  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_proj << " for option: " << option << endl);
1006  Exit(EXIT_FAILURE);
1007  }
1008  i++;
1009  }
1010  // Image convolver verbosity level
1011  else if (option=="-vb-conv")
1012  {
1013  if (i>=argc-1)
1014  {
1015  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1016  Exit(EXIT_FAILURE);
1017  }
1018  if (ConvertFromString(argv[i+1], &verbose_conv))
1019  {
1020  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_conv << " for option: " << option << endl);
1021  Exit(EXIT_FAILURE);
1022  }
1023  i++;
1024  }
1025  // Image processing verbosity level
1026  else if (option=="-vb-proc")
1027  {
1028  if (i>=argc-1)
1029  {
1030  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1031  Exit(EXIT_FAILURE);
1032  }
1033  if (ConvertFromString(argv[i+1], &verbose_proc))
1034  {
1035  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_proc << " for option: " << option << endl);
1036  Exit(EXIT_FAILURE);
1037  }
1038  i++;
1039  }
1040  // Scanner verbosity level
1041  else if (option=="-vb-scan")
1042  {
1043  if (i>=argc-1)
1044  {
1045  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1046  Exit(EXIT_FAILURE);
1047  }
1048  if (ConvertFromString(argv[i+1], &verbose_scan))
1049  {
1050  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_scan << " for option: " << option << endl);
1051  Exit(EXIT_FAILURE);
1052  }
1053  i++;
1054  }
1055  // Data and image verbosity level
1056  else if (option=="-vb-data")
1057  {
1058  if (i>=argc-1)
1059  {
1060  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1061  Exit(EXIT_FAILURE);
1062  }
1063  if (ConvertFromString(argv[i+1], &verbose_data))
1064  {
1065  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_data << " for option: " << option << endl);
1066  Exit(EXIT_FAILURE);
1067  }
1068  i++;
1069  }
1070  // Deformation verbosity level
1071  else if (option=="-vb-defo")
1072  {
1073  if (i>=argc-1)
1074  {
1075  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1076  Exit(EXIT_FAILURE);
1077  }
1078  if (ConvertFromString(argv[i+1], &verbose_defo))
1079  {
1080  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_defo << " for option: " << option << endl);
1081  Exit(EXIT_FAILURE);
1082  }
1083  i++;
1084  }
1085  // Dynamic verbosity level
1086  else if (option=="-vb-dyna")
1087  {
1088  if (i>=argc-1)
1089  {
1090  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1091  Exit(EXIT_FAILURE);
1092  }
1093  if (ConvertFromString(argv[i+1], &verbose_dyna))
1094  {
1095  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_dyna << " for option: " << option << endl);
1096  Exit(EXIT_FAILURE);
1097  }
1098  i++;
1099  }
1100  // Sensitivity verbosity level
1101  else if (option=="-vb-sens")
1102  {
1103  if (i>=argc-1)
1104  {
1105  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1106  Exit(EXIT_FAILURE);
1107  }
1108  if (ConvertFromString(argv[i+1], &verbose_sens))
1109  {
1110  Cerr("***** castor-recon() -> Exception when trying to read provided verbosity level '" << verbose_sens << " for option: " << option << endl);
1111  Exit(EXIT_FAILURE);
1112  }
1113  i++;
1114  }
1115  // RNG seed
1116  else if (option=="-rng-seed")
1117  {
1118  if (i>=argc-1)
1119  {
1120  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1121  Exit(EXIT_FAILURE);
1122  }
1123  if(ConvertFromString(argv[i+1], &random_generator_seed))
1124  {
1125  Cerr("***** castor-recon() -> Exception when trying to read provided number '" << random_generator_seed << " for option: " << option << endl);
1126  Exit(EXIT_FAILURE);
1127  }
1128  i++;
1129  }
1130  // number of extra RNG
1131  else if (option=="-rng-extra")
1132  {
1133  if (i>=argc-1)
1134  {
1135  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1136  Exit(EXIT_FAILURE);
1137  }
1138  if(ConvertFromString(argv[i+1], &nb_extra_random_generators))
1139  {
1140  Cerr("***** castor-recon() -> Exception when trying to read provided number '" << nb_extra_random_generators << " for option: " << option << endl);
1141  Exit(EXIT_FAILURE);
1142  }
1143  i++;
1144  }
1145  // Path to config directory
1146  else if (option=="-conf")
1147  {
1148  if (i>=argc-1)
1149  {
1150  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1151  Exit(EXIT_FAILURE);
1152  }
1153  path_to_config_dir = (string)argv[i+1];
1154  i++;
1155  }
1156 
1157  // --------------------------------------------------------------------------------
1158  // Dimensions settings
1159  // --------------------------------------------------------------------------------
1160 
1161  // Dimensions: number of voxels
1162  else if (option=="-dim")
1163  {
1164  if (i>=argc-1)
1165  {
1166  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1167  Exit(EXIT_FAILURE);
1168  }
1169  INTNB input[3];
1170  if (ReadStringOption(argv[i+1], input, 3, ",", option))
1171  {
1172  Cerr("***** castor-recon() -> Invalid argument " << argv[i+1] << " for option " << option << " !" << endl);
1173  Exit(EXIT_FAILURE);
1174  }
1175  nb_voxX = input[0];
1176  nb_voxY = input[1];
1177  nb_voxZ = input[2];
1178  i++;
1179  }
1180  // Dimensions: size of the field-of-view in mm
1181  else if (option=="-fov")
1182  {
1183  if (i>=argc-1)
1184  {
1185  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1186  Exit(EXIT_FAILURE);
1187  }
1188  FLTNB input[3];
1189  if (ReadStringOption(argv[i+1], input, 3, ",", option))
1190  {
1191  Cerr("***** castor-recon() -> Invalid argument " << argv[i+1] << " for option " << option << " !" << endl);
1192  Exit(EXIT_FAILURE);
1193  }
1194  fov_sizeX = input[0];
1195  fov_sizeY = input[1];
1196  fov_sizeZ = input[2];
1197  i++;
1198  }
1199  // Dimensions: size of the voxels in mm
1200  else if (option=="-vox")
1201  {
1202  if (i>=argc-1)
1203  {
1204  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1205  Exit(EXIT_FAILURE);
1206  }
1207  FLTNB input[3];
1208  if (ReadStringOption(argv[i+1], input, 3, ",", option))
1209  {
1210  Cerr("***** castor-recon() -> Invalid argument " << argv[i+1] << " for option " << option << " !" << endl);
1211  Exit(EXIT_FAILURE);
1212  }
1213  vox_sizeX = input[0];
1214  vox_sizeY = input[1];
1215  vox_sizeZ = input[2];
1216  i++;
1217  }
1218  // Dimensions: offset of the image in mm
1219  else if (option=="-off")
1220  {
1221  if (i>=argc-1)
1222  {
1223  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1224  Exit(EXIT_FAILURE);
1225  }
1226  FLTNB input[3];
1227  if (ReadStringOption(argv[i+1], input, 3, ",", option))
1228  {
1229  Cerr("***** castor-recon() -> Invalid argument " << argv[i+1] << " for option " << option << " !" << endl);
1230  Exit(EXIT_FAILURE);
1231  }
1232  offsetX = input[0];
1233  offsetY = input[1];
1234  offsetZ = input[2];
1235  i++;
1236  }
1237  // Size of the output transaxial FOV in percentage
1238  else if (option=="-fov-out")
1239  {
1240  if (i>=argc-1)
1241  {
1242  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1243  Exit(EXIT_FAILURE);
1244  }
1245  fov_out = atof(argv[i+1]);
1246  i++;
1247  }
1248  // Number of slices to be masked at output
1249  else if (option=="-slice-out")
1250  {
1251  if (i>=argc-1)
1252  {
1253  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1254  Exit(EXIT_FAILURE);
1255  }
1256  slice_out = ((INTNB)(atoi(argv[i+1])));
1257  i++;
1258  }
1259  // Output flip
1260  else if (option=="-flip-out")
1261  {
1262  if (i>=argc-1)
1263  {
1264  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1265  Exit(EXIT_FAILURE);
1266  }
1267  flip_out = (string)(argv[i+1]);
1268  i++;
1269  }
1270 
1271  // --------------------------------------------------------------------------------
1272  // Dynamic settings
1273  // --------------------------------------------------------------------------------
1274 
1275  // Dimensions: number of frames
1276  else if (option=="-frm") // TODO: more checks to be done in the option reading
1277  {
1278  if (i>=argc-1)
1279  {
1280  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1281  Exit(EXIT_FAILURE);
1282  }
1283  frame_list = (string)argv[i+1];
1284  i++;
1285  }
1286  // Dimensions: number of respiratory gates
1287  else if (option=="-g")
1288  {
1289  if (i>=argc-1)
1290  {
1291  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1292  Exit(EXIT_FAILURE);
1293  }
1294  // Path to the dynamic file
1295  path_to_4D_data_splitting_file = ((string)argv[i+1]);
1296  // Recover number of respiratory gates from file
1297  if (ReadDataASCIIFile(path_to_4D_data_splitting_file, "nb_respiratory_gates", &nb_resp_gates, 1, KEYWORD_OPTIONAL) == KEYWORD_OPTIONAL_ERROR)
1298  {
1299  Cerr("***** castor-recon() -> Error when trying to read the number of respiratory gates in the file " << path_to_4D_data_splitting_file << " for option " << option << endl);
1300  Exit(EXIT_FAILURE);
1301  }
1302  // Recover number of cardiac gates from file
1303  if (ReadDataASCIIFile(path_to_4D_data_splitting_file, "nb_cardiac_gates", &nb_card_gates, 1, KEYWORD_OPTIONAL) == KEYWORD_OPTIONAL_ERROR)
1304  {
1305  Cerr("***** castor-recon() -> Error when trying to read the number of cardiac gates in the file " << path_to_4D_data_splitting_file << " for option " << option << endl);
1306  Exit(EXIT_FAILURE);
1307  }
1308  // Check for wrong initialization
1309  if (nb_resp_gates<1 || nb_card_gates <1)
1310  {
1311  Cerr("***** castor-recon() -> Incorrect initialization of the number of gates for the option: " << option << ". This number should be >= 1" << endl);
1312  Exit(EXIT_FAILURE);
1313  }
1314  i++;
1315  }
1316  // Time basis functions
1317  else if (option=="-time-basis")
1318  {
1319  if (i>=argc-1)
1320  {
1321  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1322  Exit(EXIT_FAILURE);
1323  }
1324  string input = argv[i+1];
1325  size_t column_pos = input.find_first_of(":");
1326  if (column_pos==string::npos)
1327  {
1328  Cerr("***** castor-recon() -> Incorrect argument after option " << option << ", ':' sign is missing !" << endl);
1329  Exit(EXIT_FAILURE);
1330  }
1331  string str = input.substr(0, column_pos);
1332  nb_time_basis = atoi(str.c_str());
1333  path_to_time_basis_coef = input.substr(column_pos+1);
1334  i++;
1335  }
1336  // Respiratory basis functions
1337  else if (option=="-resp-basis")
1338  {
1339  if (i>=argc-1)
1340  {
1341  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1342  Exit(EXIT_FAILURE);
1343  }
1344  string input = argv[i+1];
1345  size_t column_pos = input.find_first_of(":");
1346  if (column_pos==string::npos)
1347  {
1348  Cerr("***** castor-recon() -> Incorrect argument after option " << option << ", ':' sign is missing !" << endl);
1349  Exit(EXIT_FAILURE);
1350  }
1351  string str = input.substr(0, column_pos);
1352  nb_resp_basis = atoi(str.c_str());
1353  path_to_resp_basis_coef = input.substr(column_pos+1);
1354  i++;
1355  }
1356  // Cardiac basis functions
1357  else if (option=="-card-basis")
1358  {
1359  if (i>=argc-1)
1360  {
1361  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1362  Exit(EXIT_FAILURE);
1363  }
1364  string input = argv[i+1];
1365  size_t column_pos = input.find_first_of(":");
1366  if (column_pos==string::npos)
1367  {
1368  Cerr("***** castor-recon() -> Incorrect argument after option " << option << ", ':' sign is missing !" << endl);
1369  Exit(EXIT_FAILURE);
1370  }
1371  string str = input.substr(0, column_pos);
1372  nb_card_basis = atoi(str.c_str());
1373  path_to_card_basis_coef = input.substr(column_pos+1);
1374  i++;
1375  }
1376  // Dynamic model applied to the frames/respiratory gates/cardiac gates of the dynamic acquisition
1377  else if (option=="-dynamic-model")
1378  {
1379  if (i>=argc-1)
1380  {
1381  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1382  Exit(EXIT_FAILURE);
1383  }
1384  dynamic_model_options = (string)argv[i+1];
1385  i++;
1386  }
1387  // Data for respiratory motion correction based on deformation
1388  else if (option=="-rm")
1389  {
1390  if (i>=argc-1)
1391  {
1392  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1393  Exit(EXIT_FAILURE);
1394  }
1395  resp_motion_options = (string)argv[i+1];
1396  i++;
1397  }
1398  // Data for cardiac motion correction based on deformation
1399  else if (option=="-cm")
1400  {
1401  if (i>=argc-1)
1402  {
1403  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1404  Exit(EXIT_FAILURE);
1405  }
1406  card_motion_options = (string)argv[i+1];
1407  i++;
1408  }
1409  // Data for respiratory and cardiac motion corrections based on deformation
1410  else if (option=="-rcm")
1411  {
1412  if (i>=argc-1)
1413  {
1414  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1415  Exit(EXIT_FAILURE);
1416  }
1417  double_motion_options = (string)argv[i+1];
1418  i++;
1419  }
1420  // Data for involuntary motion correction based on deformation
1421  else if (option=="-im")
1422  {
1423  if (i>=argc-1)
1424  {
1425  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1426  Exit(EXIT_FAILURE);
1427  }
1428  ipat_motion_options = (string)argv[i+1];
1429  i++;
1430  }
1431  // Data for involuntary motion correction based on deformation
1432  else if (option=="-qdyn")
1433  {
1434  if (i>=argc-1)
1435  {
1436  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1437  Exit(EXIT_FAILURE);
1438  }
1439  path_to_dynamic_quantification_file = (string)argv[i+1];
1440  i++;
1441  }
1442 
1443  // --------------------------------------------------------------------------------
1444  // Input settings
1445  // --------------------------------------------------------------------------------
1446 
1447  // DataFiles
1448  else if (option=="-df") // This is a mandatory option
1449  {
1450  if (i>=argc-1)
1451  {
1452  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1453  Exit(EXIT_FAILURE);
1454  }
1455  string file_name = (string)argv[i+1];
1456  path_to_data_filename.push_back(file_name);
1457  nb_beds++;
1458  i++;
1459  }
1460  // Invert datafile order
1461  else if (option=="-df-inv")
1462  {
1463  invert_datafile_bed_order_flag = true;
1464  }
1465  // Image for the initialisation of the algorithm : What should we do in case of multiple beds ? or multiple frames ? TODO
1466  else if (option=="-img")
1467  {
1468  if (i>=argc-1)
1469  {
1470  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1471  Exit(EXIT_FAILURE);
1472  }
1473  path_to_initial_img = argv[i+1];
1474  i++;
1475  }
1476  // Sensitivity image
1477  else if (option=="-sens")
1478  {
1479  if (i>=argc-1)
1480  {
1481  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1482  Exit(EXIT_FAILURE);
1483  }
1484 
1485  path_to_sensitivity_img = argv[i+1];
1486  i++;
1487  }
1488  // Normalization data files
1489  else if (option=="-norm")
1490  {
1491  if (i>=argc-1)
1492  {
1493  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1494  Exit(EXIT_FAILURE);
1495  }
1496  string file_name = (string)argv[i+1];
1497  path_to_normalization_filename.push_back(file_name);
1498  i++;
1499  }
1500 
1501  // Image for the attenuation
1502  else if (option=="-atn")
1503  {
1504  if (i>=argc-1)
1505  {
1506  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1507  Exit(EXIT_FAILURE);
1508  }
1509 
1510  path_to_attenuation_img = (string)argv[i+1];
1511 
1512  if (path_to_attenuation_img != "") // parse
1513  {
1514  // Retrieve nb of gated images from header if required
1515  Intf_fields IF;
1516  IntfKeyInitFields(&IF);
1517  if(IntfReadHeader(path_to_attenuation_img, &IF, 0) )
1518  {
1519  Cerr("***** castor-recon() -> An error occurred while trying to read the interfile header of attenuation file " << path_to_attenuation_img << " !" << endl);
1520  Exit(EXIT_FAILURE);
1521  }
1522  nb_atn_resp_imgs = IF.nb_resp_gates;
1523  nb_atn_card_imgs = IF.nb_card_gates;
1524  }
1525  i++;
1526  }
1527  // Multimodal images
1528  else if (option=="-multimodal")
1529  {
1530  if (i>=argc-1)
1531  {
1532  Cerr("***** Argument missing for option: " << option << endl);
1533  Exit(EXIT_FAILURE);
1534  }
1535  string file_name = (string)argv[i+1];
1536  path_to_multimodal_img.push_back(file_name);
1537  i++;
1538  }
1539  // Mask image
1540  else if (option=="-mask")
1541  {
1542  if (i>=argc-1)
1543  {
1544  Cerr("***** Argument missing for option: " << option << endl);
1545  Exit(EXIT_FAILURE);
1546  }
1547 
1548  path_to_mask_img = argv[i+1];
1549  i++;
1550  }
1551 
1552  // --------------------------------------------------------------------------------
1553  // Output settings
1554  // --------------------------------------------------------------------------------
1555 
1556  // Name of the output directory
1557  else if (option=="-dout") // This is a mandatory option
1558  {
1559  if (i>=argc-1)
1560  {
1561  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1562  Exit(EXIT_FAILURE);
1563  }
1564  path_dout = argv[i+1];
1565  i++;
1566  }
1567  // Base name of the output files
1568  else if (option=="-fout") // This is a mandatory option
1569  {
1570  if (i>=argc-1)
1571  {
1572  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1573  Exit(EXIT_FAILURE);
1574  }
1575  path_fout = argv[i+1];
1576  i++;
1577  }
1578  // List of iterations to be outputed
1579  else if (option=="-oit")
1580  {
1581  if (i>=argc-1)
1582  {
1583  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1584  Exit(EXIT_FAILURE);
1585  }
1586  output_iterations = (string)argv[i+1];
1587  i++;
1588  }
1589  // Output number precision
1590  else if (option=="-onbp")
1591  {
1592  if (i>=argc-1)
1593  {
1594  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1595  Exit(EXIT_FAILURE);
1596  }
1597  onb_prec = argv[i+1];
1598  i++;
1599  }
1600  // Flag to say that we want to save time basis functions too
1601  else if (option=="-omd")
1602  {
1603  merge_dynamic_imgs_flag = true;
1604  }
1605  // Flag to say that we want to save time basis functions too
1606  else if (option=="-olut")
1607  {
1608  save_LUT_flag = true;
1609  }
1610  // Flag to say that we want to save the sensitivity image for each subset/iteration in histogram mode
1611  else if (option=="-osens")
1612  {
1613  save_sens_histo = true;
1614  }
1615  // Flag to say that we want to save the image after each subset
1616  else if (option=="-osub")
1617  {
1618  save_subset_image = true;
1619  }
1620  // Flag to say that we want to save time basis functions too
1621  else if (option=="-otb")
1622  {
1623  Cerr("!!!!! castor-recon() -> Warning: option -otb not yet implemented !" << endl);
1624  }
1625  // Flag to say that we exit directly after computing and saving the sensitivity
1626  else if (option=="-sens-only")
1627  {
1628  exit_after_sensitivity = true;
1629  }
1630  // Flag to say that we still want to compute global sensitivity from the input histogram file
1631  else if (option=="-sens-histo")
1632  {
1633  sensitivity_from_histogram = true;
1634  }
1635 
1636  // --------------------------------------------------------------------------------
1637  // Algorithm settings
1638  // --------------------------------------------------------------------------------
1639 
1640  // List of iterations/subsets
1641  else if (option=="-it") // This is a mandatory option
1642  {
1643  if (i>=argc-1)
1644  {
1645  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1646  Exit(EXIT_FAILURE);
1647  }
1648  nb_iterations_subsets = (string)argv[i+1];
1649  i++;
1650  }
1651  // Optimizer settings
1652  else if (option=="-opti")
1653  {
1654  if (i>=argc-1)
1655  {
1656  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1657  Exit(EXIT_FAILURE);
1658  }
1659  options_optimizer = (string)argv[i+1];
1660  i++;
1661  }
1662  // Optimizer FOM
1663  else if (option=="-opti-fom")
1664  {
1665  optimizer_fom = true;
1666  }
1667  // Optimizer image update stat
1668  else if (option=="-opti-stat")
1669  {
1670  optimizer_stat = true;
1671  }
1672  // Penalty settings
1673  else if (option=="-penalty")
1674  {
1675  if (i>=argc-1)
1676  {
1677  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1678  Exit(EXIT_FAILURE);
1679  }
1680  options_penalty = (string)argv[i+1];
1681  i++;
1682  }
1683  // Image convolver settings
1684  else if (option=="-conv")
1685  {
1686  if (i>=argc-1)
1687  {
1688  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1689  Exit(EXIT_FAILURE);
1690  }
1691  string convolver = (string)argv[i+1];
1692  options_image_convolver.push_back(convolver);
1693  i++;
1694  }
1695  // Image processing settings
1696  else if (option=="-proc")
1697  {
1698  if (i>=argc-1)
1699  {
1700  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1701  Exit(EXIT_FAILURE);
1702  }
1703  string module = (string)argv[i+1];
1704  options_image_processing.push_back(module);
1705  i++;
1706  }
1707 
1708  // --------------------------------------------------------------------------------
1709  // Projection settings
1710  // --------------------------------------------------------------------------------
1711 
1712  // Projector settings
1713  else if (option=="-proj")
1714  {
1715  if (i>=argc-1)
1716  {
1717  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1718  Exit(EXIT_FAILURE);
1719  }
1720  options_projectorF = (string)argv[i+1];
1721  options_projectorB = (string)argv[i+1];
1722  i++;
1723  }
1724  // Forward projector settings
1725  else if (option=="-projF")
1726  {
1727  if (i>=argc-1)
1728  {
1729  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1730  Exit(EXIT_FAILURE);
1731  }
1732  options_projectorF = (string)argv[i+1];
1733  i++;
1734  }
1735  // Backward projector settings
1736  else if (option=="-projB")
1737  {
1738  if (i>=argc-1)
1739  {
1740  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1741  Exit(EXIT_FAILURE);
1742  }
1743  options_projectorB = (string)argv[i+1];
1744  i++;
1745  }
1746  // Common projector settings
1747  else if (option=="-proj-common")
1748  {
1749  if (i>=argc-1)
1750  {
1751  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1752  Exit(EXIT_FAILURE);
1753  }
1754  options_projector_common = (string)argv[i+1];
1755  i++;
1756  }
1757  // Ignore TOF flag
1758  else if (option=="-ignore-TOF")
1759  {
1760  ignore_TOF = true;
1761  }
1762  // Ignore POI flag
1763  else if (option=="-ignore-POI")
1764  {
1765  ignore_POI = true;
1766  }
1767  // Projection line computation strategy
1768  else if (option=="-proj-comp")
1769  {
1770  if (i>=argc-1)
1771  {
1772  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1773  Exit(EXIT_FAILURE);
1774  }
1775  projector_computation_strategy = atoi(argv[i+1]);
1776  i++;
1777  }
1778 
1779  // --------------------------------------------------------------------------------
1780  // Quantification settings
1781  // --------------------------------------------------------------------------------
1782 
1783  // Corrections settings
1784  else if (option=="-ignore-corr")
1785  {
1786  if (i>=argc-1)
1787  {
1788  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1789  Exit(EXIT_FAILURE);
1790  }
1791  ignored_corrections = (string)argv[i+1];
1792  i++;
1793  }
1794 
1795  // --------------------------------------------------------------------------------
1796  // Computation settings
1797  // --------------------------------------------------------------------------------
1798 
1799  // Flag to say that we want to use the GPU
1800  #ifdef CASTOR_GPU
1801  else if (option=="-gpu")
1802  {
1803  gpu_flag = 1;
1804  }
1805  #endif
1806  // Number of threads
1807  #ifdef CASTOR_OMP
1808  else if (option=="-th")
1809  {
1810  if (i>=argc-1)
1811  {
1812  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1813  Exit(EXIT_FAILURE);
1814  }
1815  nb_threads = (string)argv[i+1];
1816  i++;
1817  }
1818  #else
1819  else if (option=="-th")
1820  {
1821  if (i>=argc-1)
1822  {
1823  Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
1824  Exit(EXIT_FAILURE);
1825  }
1826  Cerr("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
1827  Cerr("!!!!! castor-recon() -> Option -th is available only if the code is compiled using the CASTOR_OMP environment variable set to 1 !" << endl);
1828  Cerr("!!!!! The execution will continue BUT WITH ONLY ONE THREAD !" << endl);
1829  Cerr("!!!!! We strongly advice to compile CASToR with OpenMP !" << endl);
1830  Cerr("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
1831  i++;
1832  }
1833  #endif
1834 
1835  // --------------------------------------------------------------------------------
1836  // Unknown option!
1837  // --------------------------------------------------------------------------------
1838 
1839  else
1840  {
1841  Cerr("***** castor-recon() -> Unknown option '" << option << "' !" << endl);
1842  Exit(EXIT_FAILURE);
1843  }
1844  }
1845 
1846  // Synchronize MPI processes
1847  #ifdef CASTOR_MPI
1848  MPI_Barrier(MPI_COMM_WORLD);
1849  #endif
1850 
1851  // Affect specific verbose leves if not set
1852  if (verbose_algo==-1) verbose_algo = verbose_general;
1853  if (verbose_opti==-1) verbose_opti = verbose_general;
1854  if (verbose_proj==-1) verbose_proj = verbose_general;
1855  if (verbose_conv==-1) verbose_conv = verbose_general;
1856  if (verbose_proc==-1) verbose_proc = verbose_general;
1857  if (verbose_scan==-1) verbose_scan = verbose_general;
1858  if (verbose_data==-1) verbose_data = verbose_general;
1859  if (verbose_defo==-1) verbose_defo = verbose_general;
1860  if (verbose_dyna==-1) verbose_dyna = verbose_general;
1861  if (verbose_sens==-1) verbose_sens = verbose_general;
1862 
1863  // ============================================================================================================
1864  // Some checks
1865  // ============================================================================================================
1866 
1867  // Data files
1868  if (nb_beds < 1)
1869  {
1870  Cerr("***** castor-recon() -> Please provide at least one data filename !" << endl);
1871  Exit(EXIT_FAILURE);
1872  }
1873  // Output files
1874  if (path_fout.empty() && path_dout.empty())
1875  {
1876  Cerr("***** castor-recon() -> Please provide an output option for output files (-fout or -dout) !" << endl);
1877  Exit(EXIT_FAILURE);
1878  }
1879  // Check that only one option has been provided
1880  if (!path_fout.empty() && !path_dout.empty())
1881  {
1882  Cerr("***** castor-recon() -> Please provide either output option -fout or -dout but not both !" << endl);
1883  Exit(EXIT_FAILURE);
1884  }
1885 
1886  // Check if gated reconstruction is enabled but no file describing the gating of the data has been provided
1887  // TODO: maybe do it in the dynamic data manager if possible to clean this main as much as possible
1888  if ( (nb_resp_gates>1 || nb_card_gates>1) && path_to_4D_data_splitting_file.empty() )
1889  {
1890  Cerr("***** castor-recon() -> gating is enabled, but no file describing the splitting of the data has been provided (-gating option) !" << endl);
1891  Exit(EXIT_FAILURE);
1892  }
1893 
1894  // ============================================================================================================
1895  // Singletons initialization: create here all needed singletons
1896  // ============================================================================================================
1897 
1898  if (verbose_general>=5) Cout("----- Singletons initializations ... -----" << endl);
1899 
1900  // Get user endianness (interfile I/O)
1902 
1903  // ----------------------------------------------------------------------------------------
1904  // Create sOutputManager
1905  // ----------------------------------------------------------------------------------------
1906  sOutputManager* p_outputManager = sOutputManager::GetInstance();
1907  // Set verbose level
1908  p_outputManager->SetVerbose(verbose_general);
1909  // Set MPI rank
1910  p_outputManager->SetMPIRank(mpi_rank);
1911  // Set output dynamic image policy
1912  p_outputManager->SetMergeDynImagesFlag(merge_dynamic_imgs_flag);
1913  // Set datafile name(s) in order to be able to recover them from interfile classes
1914  p_outputManager->SetDataFileName(path_to_data_filename);
1915  // Set output number precision
1916  p_outputManager->SetOutNbPrec(onb_prec);
1917 
1918 
1919  // Set path to the config directory
1920  if (p_outputManager->CheckConfigDir(path_to_config_dir))
1921  {
1922  Cerr("***** castor-recon() -> A problem occured while checking for the config directory path !" << endl);
1923  Exit(EXIT_FAILURE);
1924  }
1925  // Initialize output directory and base name
1926  if (p_outputManager->InitOutputDirectory(path_fout, path_dout))
1927  {
1928  Cerr("***** castor-recon() -> A problem occured while initializing output directory !" << endl);
1929  Exit(EXIT_FAILURE);
1930  }
1931  // General call verbose
1932  #ifdef CASTOR_VERSION
1933  Cout("castor-recon() -> Launch reconstruction from CASToR version " << CASTOR_VERSION << "." << endl);
1934  #endif
1935  // Log command line
1936  if (p_outputManager->LogCommandLine(argc,argv))
1937  {
1938  Cerr("***** castor-recon() -> A problem occured while logging command line arguments !" << endl);
1939  Exit(EXIT_FAILURE);
1940  }
1941 
1942  // ----------------------------------------------------------------------------------------
1943  // Create oImageDimensionsAndQuantification
1944  // ----------------------------------------------------------------------------------------
1945 
1946  oImageDimensionsAndQuantification* p_ImageDimensionsAndQuantification = new oImageDimensionsAndQuantification();
1947 
1948  // ----------------------------------------------------------------------------------------
1949  // Create sScannerManager
1950  // ----------------------------------------------------------------------------------------
1951  sScannerManager* p_ScannerManager = sScannerManager::GetInstance();
1952  p_ScannerManager->SetVerbose(verbose_scan);
1953  p_ScannerManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
1954  p_ScannerManager->SetSaveLUTFlag(save_LUT_flag);
1955 
1956  if (verbose_general>=5) Cout("----- Geometry Initialization ... -----" << endl);
1957 
1958  // TODO: put all that stuff into only one function taking the path_to_data_filename[0] as the only parameter
1959 
1960  // Get system name from the dataFile
1961  string scanner_name = "";
1962  if (ReadDataASCIIFile(path_to_data_filename[0], "Scanner name", &scanner_name, 1, KEYWORD_MANDATORY))
1963  {
1964  Cerr("***** castor-recon() -> A problem occured while trying to find the system name in the datafile header !" << endl);
1965  Exit(EXIT_FAILURE);
1966  }
1967  if (p_ScannerManager->FindScannerSystem(scanner_name) )
1968  {
1969  Cerr("***** castor-recon() -> A problem occurred while searching for scanner system !" << endl);
1970  Exit(EXIT_FAILURE);
1971  }
1972  if (p_ScannerManager->BuildScannerObject() )
1973  {
1974  Cerr("***** castor-recon() -> A problem occurred during scanner object construction ! !" << endl);
1975  Exit(EXIT_FAILURE);
1976  }
1977  if (p_ScannerManager->InstantiateScanner() )
1978  {
1979  Cerr("***** castor-recon() -> A problem occurred while creating Scanner object !" << endl);
1980  Exit(EXIT_FAILURE);
1981  }
1982  if (p_ScannerManager->GetGeometricInfoFromDataFile(path_to_data_filename[0]))
1983  {
1984  Cerr("***** castor-recon() -> A problem occurred while retrieving scanner fields from the datafile header !" << endl);
1985  Exit(EXIT_FAILURE);
1986  }
1987  if (p_ScannerManager->BuildLUT() )
1988  {
1989  Cerr("***** castor-recon() -> A problem occurred while generating/reading the LUT !" << endl);
1990  Exit(EXIT_FAILURE);
1991  }
1992  // Check the scanner manager parameters and initialize the scanner
1993  if (p_ScannerManager->CheckParameters())
1994  {
1995  Cerr("***** castor-recon() -> A problem occured while checking scanner manager parameters !" << endl);
1996  Exit(EXIT_FAILURE);
1997  }
1998  if (p_ScannerManager->Initialize())
1999  {
2000  Cerr("***** castor-recon() -> A problem occured while initializing scanner !" << endl);
2001  Exit(EXIT_FAILURE);
2002  }
2003 
2004  // If no number of voxels provided, then get the default ones from the scanner
2005  if (nb_voxX<=0 || nb_voxY<=0 || nb_voxZ<=0)
2006  {
2007  if (ReadDataASCIIFile(p_ScannerManager->GetPathToScannerFile(), "voxels number transaxial", &nb_voxX, 1, KEYWORD_MANDATORY))
2008  {
2009  Cerr("***** castor-recon() -> A problem occured while reading for default number of transaxial voxels !" << endl);
2010  Exit(EXIT_FAILURE);
2011  }
2012  if (ReadDataASCIIFile(p_ScannerManager->GetPathToScannerFile(), "voxels number transaxial", &nb_voxY, 1, KEYWORD_MANDATORY))
2013  {
2014  Cerr("***** castor-recon() -> A problem occured while reading for default number of transaxial voxels !" << endl);
2015  Exit(EXIT_FAILURE);
2016  }
2017  if (ReadDataASCIIFile(p_ScannerManager->GetPathToScannerFile(), "voxels number axial", &nb_voxZ, 1, KEYWORD_MANDATORY))
2018  {
2019  Cerr("***** castor-recon() -> A problem occured while reading for default number of axial voxels !" << endl);
2020  Exit(EXIT_FAILURE);
2021  }
2022  }
2023  // If no FOV nor VOX size provided, then get the default one
2024  if ( (fov_sizeX<=0 || fov_sizeY<=0 || fov_sizeZ<=0) && (vox_sizeX<=0 || vox_sizeY<=0 || vox_sizeZ<=0) )
2025  {
2026  if (ReadDataASCIIFile(p_ScannerManager->GetPathToScannerFile(), "field of view transaxial", &fov_sizeX, 1, KEYWORD_MANDATORY))
2027  {
2028  Cerr("***** castor-recon() -> A problem occured while reading for default transaxial FOV size !" << endl);
2029  Exit(EXIT_FAILURE);
2030  }
2031  if (ReadDataASCIIFile(p_ScannerManager->GetPathToScannerFile(), "field of view transaxial", &fov_sizeY, 1, KEYWORD_MANDATORY))
2032  {
2033  Cerr("***** castor-recon() -> A problem occured while reading for default transaxial FOV size !" << endl);
2034  Exit(EXIT_FAILURE);
2035  }
2036  if (ReadDataASCIIFile(p_ScannerManager->GetPathToScannerFile(), "field of view axial", &fov_sizeZ, 1, KEYWORD_MANDATORY))
2037  {
2038  Cerr("***** castor-recon() -> A problem occured while reading for default axial FOV size !" << endl);
2039  Exit(EXIT_FAILURE);
2040  }
2041  }
2042 
2043  if (verbose_general>=5) Cout("----- Geometry Initialization OK -----" << endl);
2044 
2045  // ----------------------------------------------------------------------------------------
2046  // Create oImageDimensionsAndQuantification
2047  // ----------------------------------------------------------------------------------------
2048 
2049  if (verbose_general>=5) Cout("----- Image dimensions initialization ... -----" << endl);
2050 
2051  if (p_ImageDimensionsAndQuantification->SetNbThreads(nb_threads))
2052  {
2053  Cerr("***** castor-recon() -> A problem occured while setting the number of threads !" << endl);
2054  Exit(EXIT_FAILURE);
2055  }
2056  p_ImageDimensionsAndQuantification->SetNbBeds(nb_beds);
2057  p_ImageDimensionsAndQuantification->SetNbVoxX(nb_voxX);
2058  p_ImageDimensionsAndQuantification->SetNbVoxY(nb_voxY);
2059  p_ImageDimensionsAndQuantification->SetNbVoxZ(nb_voxZ);
2060  p_ImageDimensionsAndQuantification->SetVoxSizeX(vox_sizeX);
2061  p_ImageDimensionsAndQuantification->SetVoxSizeY(vox_sizeY);
2062  p_ImageDimensionsAndQuantification->SetVoxSizeZ(vox_sizeZ);
2063  p_ImageDimensionsAndQuantification->SetFOVSizeX(fov_sizeX);
2064  p_ImageDimensionsAndQuantification->SetFOVSizeY(fov_sizeY);
2065  p_ImageDimensionsAndQuantification->SetFOVSizeZ(fov_sizeZ);
2066  p_ImageDimensionsAndQuantification->SetFOVOutMasking(fov_out,slice_out);
2067  p_ImageDimensionsAndQuantification->SetOffsetX(offsetX);
2068  p_ImageDimensionsAndQuantification->SetOffsetY(offsetY);
2069  p_ImageDimensionsAndQuantification->SetOffsetZ(offsetZ);
2070  p_ImageDimensionsAndQuantification->SetMPIRankAndSize(mpi_rank, mpi_size);
2071  p_ImageDimensionsAndQuantification->SetVerbose(verbose_data);
2072  p_ImageDimensionsAndQuantification->SetIgnoredCorrections(ignored_corrections);
2073  p_ImageDimensionsAndQuantification->SetFrames(frame_list);
2074  p_ImageDimensionsAndQuantification->SetNbTimeBasisFunctions(nb_time_basis);
2075  p_ImageDimensionsAndQuantification->SetTimeBasisFunctionsFile(path_to_time_basis_coef);
2076  p_ImageDimensionsAndQuantification->SetNbMultiModalImages(path_to_multimodal_img.size());
2077  if (p_ImageDimensionsAndQuantification->SetFlipOut(flip_out))
2078  {
2079  Cerr("***** castor-recon() -> A problem occured while setting the output flip option !" << endl);
2080  Exit(EXIT_FAILURE);
2081  }
2082  if (resp_motion_options=="" && double_motion_options=="")
2083  {
2084  p_ImageDimensionsAndQuantification->SetNbRespGates(nb_resp_gates);
2085  p_ImageDimensionsAndQuantification->SetNbRespBasisFunctions(nb_resp_basis);
2086  p_ImageDimensionsAndQuantification->SetRespBasisFunctionsFile(path_to_resp_basis_coef);
2087  }
2088  else
2089  {
2090  if (path_to_resp_basis_coef!="")
2091  {
2092  Cerr("***** castor-recon() -> Cannot use both respiratory motion correction and respiratory basis functions, it has no sense !" << endl);
2093  Exit(EXIT_FAILURE);
2094  }
2095  // Set only one gate here because we correct for motion, so we reconstruct only one image
2096  p_ImageDimensionsAndQuantification->SetNbRespGates(1);
2097  }
2098  if (card_motion_options=="" && double_motion_options=="")
2099  {
2100  p_ImageDimensionsAndQuantification->SetNbCardGates(nb_card_gates);
2101  p_ImageDimensionsAndQuantification->SetNbCardBasisFunctions(nb_card_basis);
2102  p_ImageDimensionsAndQuantification->SetCardBasisFunctionsFile(path_to_card_basis_coef);
2103  }
2104  else
2105  {
2106  if (path_to_card_basis_coef!="")
2107  {
2108  Cerr("***** castor-recon() -> Cannot use both cardiac motion correction and cardiac basis functions, it has no sense !" << endl);
2109  Exit(EXIT_FAILURE);
2110  }
2111  // Set only one gate here because we correct for motion, so we reconstruct only one image
2112  p_ImageDimensionsAndQuantification->SetNbCardGates(1);
2113  }
2114  if (p_ImageDimensionsAndQuantification->CheckParameters())
2115  {
2116  Cerr("***** castor-recon() -> A problem occured while checking image dimensions parameters !" << endl);
2117  Exit(EXIT_FAILURE);
2118  }
2119  if (p_ImageDimensionsAndQuantification->Initialize())
2120  {
2121  Cerr("***** castor-recon() -> A problem occured while initializing image dimensions !" << endl);
2122  Exit(EXIT_FAILURE);
2123  }
2124  // Initialization of DynamicDataManager class, related 4D data splitting management
2125  if (p_ImageDimensionsAndQuantification->InitDynamicData(path_to_4D_data_splitting_file,
2126  !resp_motion_options.empty(),
2127  !card_motion_options.empty(),
2128  !double_motion_options.empty(),
2129  !ipat_motion_options.empty(),
2130  nb_resp_gates,
2131  nb_card_gates) )
2132  {
2133  Cerr("***** castor-recon() -> A problem occured while initializing Dynamic data manager's class !" << endl);
2134  Exit(EXIT_FAILURE);
2135  }
2136  // Get the number of events in the data from the header (in order to check consistency between
2137  // the number of events in the datafile and in the dynamic files, if any gating is enabled)
2138  int64_t nb_events = 0;
2139  ReadDataASCIIFile(path_to_data_filename[0], "Number of events", &nb_events, 1, KEYWORD_MANDATORY);
2140  // Check dynamic parameters
2141  if (p_ImageDimensionsAndQuantification->CheckDynamicParameters(nb_events) )
2142  {
2143  Cerr("***** castor-recon() -> A problem occured while checking Dynamic data manager's parameters !" << endl);
2144  Exit(EXIT_FAILURE);
2145  }
2146  // Initialize dynamic specific quantitative factors
2147  if (p_ImageDimensionsAndQuantification->SetDynamicSpecificQuantificationFactors(path_to_dynamic_quantification_file) )
2148  {
2149  Cerr("***** castor-recon() -> A problem occured while initializing specific dynamic quantification factors!" << endl);
2150  Exit(EXIT_FAILURE);
2151  }
2152 
2153  if (verbose_general>=5) Cout("----- Image dimensions initialization OK -----" << endl);
2154 
2155  // ----------------------------------------------------------------------------------------
2156  // Create the image space
2157  // ----------------------------------------------------------------------------------------
2158 
2159  oImageSpace* p_ImageSpace = new oImageSpace();
2160  p_ImageSpace->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2161  p_ImageSpace->SetVerbose(verbose_data);
2162 
2163  // ----------------------------------------------------------------------------------------
2164  // Create sChronoManager
2165  // ----------------------------------------------------------------------------------------
2166 
2167  if (verbose_general>=5) Cout("----- Profiling manager initialization ... -----" << endl);
2168 
2169  sChronoManager* p_ChronoManager = sChronoManager::GetInstance();
2170  p_ChronoManager->SetNbThreads( p_ImageDimensionsAndQuantification->GetNbThreadsForProjection(),
2171  p_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation() );
2172  p_ChronoManager->SetVerbose(0);
2173  if (p_ChronoManager->CheckParameters())
2174  {
2175  Cerr("***** castor-recon() -> A problem occured while checking parameters for the chrono manager !" << endl);
2176  Exit(EXIT_FAILURE);
2177  }
2178  if (p_ChronoManager->Initialize())
2179  {
2180  Cerr("***** castor-recon() -> A problem occured while initializing the chrono manager !" << endl);
2181  Exit(EXIT_FAILURE);
2182  }
2183 
2184  if (verbose_general>=5) Cout("----- Profiling manager initialization OK -----" << endl);
2185 
2186  // ----------------------------------------------------------------------------------------
2187  // Random Number Generator initialization: (we first require to know the number of threads to use from p_ImageDimensionsAndQuantification)
2188  // ----------------------------------------------------------------------------------------
2189 
2190  if (verbose_general>=5) Cout("----- Random number generator initialization ... -----" << endl);
2191 
2192  sRandomNumberGenerator* p_RandomNumberGenerator = sRandomNumberGenerator::GetInstance();
2193  p_RandomNumberGenerator->SetVerbose(verbose_general);
2194  // Use a user-provided seed to initialize the RNG if one has been provided. Use random number otherwise.
2195  if (random_generator_seed>=0) p_RandomNumberGenerator->Initialize(random_generator_seed, p_ImageDimensionsAndQuantification->GetNbThreadsMax(), nb_extra_random_generators);
2196  else p_RandomNumberGenerator->Initialize(p_ImageDimensionsAndQuantification->GetNbThreadsMax(), nb_extra_random_generators);
2197 
2198  if (verbose_general >=5) Cout("----- Random number generator initialization OK -----" << endl);
2199 
2200  // ----------------------------------------------------------------------------------------
2201  // Create vDataFile
2202  // ----------------------------------------------------------------------------------------
2203 
2204  if (verbose_general>=5) Cout("----- DataFile initialization ... -----" << endl);
2205 
2206  vDataFile** p_DataFile = new vDataFile*[nb_beds];
2207 
2208  if (p_ScannerManager->GetScannerType() == SCANNER_PET)
2209  {
2210  // Create specific data file
2211  for (int i=0 ; i<nb_beds ; i++)
2212  {
2213  p_DataFile[i] = new iDataFilePET();
2214  (dynamic_cast<iDataFilePET*>(p_DataFile[i]))->SetIgnoreTOFFlag(ignore_TOF);
2215  }
2216  }
2217  else if (p_ScannerManager->GetScannerType() == SCANNER_SPECT_CONVERGENT)
2218  {
2219  // Create specific data file
2220  for (int i=0 ; i<nb_beds ; i++)
2221  {
2222  p_DataFile[i] = new iDataFileSPECT();
2223  }
2224  }
2225  else if (p_ScannerManager->GetScannerType() == SCANNER_CT)
2226  {
2227  // Create specific data file
2228  for (int i=0 ; i<nb_beds ; i++)
2229  {
2230  p_DataFile[i] = new iDataFileCT();
2231  }
2232  }
2233  // Unknown scanner
2234  else
2235  {
2236  Cerr("***** castor-recon() -> Unknown scanner type (" << p_ScannerManager->GetScannerType() << ") for datafile construction ! Abort." << endl);
2237  Exit(EXIT_FAILURE);
2238  }
2239 
2240  // Load raw data in memory and do other stuff if needed.
2241  for (int bed=0 ; bed<nb_beds ; bed++)
2242  {
2243  // If it is asked to revert the order of the bed positions, then we proceed here
2244  if (invert_datafile_bed_order_flag) p_DataFile[bed]->SetHeaderDataFileName(path_to_data_filename.at(nb_beds-1-bed));
2245  else p_DataFile[bed]->SetHeaderDataFileName(path_to_data_filename.at(bed));
2246  p_DataFile[bed]->SetBedIndex(bed);
2247  p_DataFile[bed]->SetVerbose(verbose_data);
2248  p_DataFile[bed]->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2249  p_DataFile[bed]->SetIgnorePOIFlag(ignore_POI);
2250  if (p_DataFile[bed]->ReadInfoInHeader())
2251  {
2252  Cerr("***** castor-recon() -> A problem occurred during datafile header reading ! Abort." << endl);
2253  Exit(EXIT_FAILURE);
2254  }
2255  if (p_DataFile[bed]->CheckParameters())
2256  {
2257  Cerr("***** castor-recon() -> A problem occurred while checking datafile parameters ! Abort." << endl);
2258  Exit(EXIT_FAILURE);
2259  }
2260  if (p_DataFile[bed]->ComputeSizeEvent())
2261  {
2262  Cerr("***** castor-recon() -> A problem occurred in datafile initialization ! Abort." << endl);
2263  Exit(EXIT_FAILURE);
2264  }
2265  if (p_DataFile[bed]->InitializeMappedFile())
2266  {
2267  Cerr("***** castor-recon() -> A problem occurred in datafile initialization ! Abort." << endl);
2268  Exit(EXIT_FAILURE);
2269  }
2270  if (p_DataFile[bed]->PrepareDataFile())
2271  {
2272  Cerr("***** castor-recon() -> A problem occured in datafile preparation ! Abort." << endl);
2273  Exit(EXIT_FAILURE);
2274  }
2275  }
2276  // Check consistency between all datafiles; to do that we compare the first with all the others
2277  for (int bed=1; bed<nb_beds; bed++)
2278  {
2279  if (p_DataFile[0]->CheckConsistencyWithAnotherBedDataFile(p_DataFile[bed]))
2280  {
2281  int bed_index_problem = bed + 1;
2282  if (invert_datafile_bed_order_flag) bed_index_problem = nb_beds - bed;
2283  Cerr("***** castor-recon() -> A problem occured while checking consistency between first bed and bed " << bed_index_problem << " !" << endl);
2284  Exit(EXIT_FAILURE);
2285  }
2286  }
2287  // And finally, deal with the multiple bed positions
2288  if (p_ImageDimensionsAndQuantification->DealWithBedPositions(p_DataFile))
2289  {
2290  Cerr("***** castor-recon() -> A problem occured while dealing with the bed positions !" << endl);
2291  Exit(EXIT_FAILURE);
2292  }
2293 
2294  if (verbose_general>=5) Cout("----- DataFile initialization OK -----" << endl);
2295 
2296  // ----------------------------------------------------------------------------------------
2297  // Create Projector Manager
2298  // ----------------------------------------------------------------------------------------
2299 
2300  // Verbose
2301  if (verbose_general>=5) Cout("----- Projector initialization ... -----" << endl);
2302  // Create object
2303  oProjectorManager* p_ProjectorManager = new oProjectorManager();
2304  // Set all parameters
2305  p_ProjectorManager->SetScanner(p_ScannerManager->GetScannerObject());
2306  p_ProjectorManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2307  p_ProjectorManager->SetDataFile(p_DataFile[0]);
2308  p_ProjectorManager->SetComputationStrategy(projector_computation_strategy);
2309  p_ProjectorManager->SetOptionsForward(options_projectorF);
2310  p_ProjectorManager->SetOptionsBackward(options_projectorB);
2311  p_ProjectorManager->SetOptionsCommon(options_projector_common);
2312 
2313  if (!path_to_mask_img.empty())
2314  {
2315  if (p_ImageSpace->InitMaskImage(path_to_mask_img))
2316  {
2317  Cerr("The option for using a mask image for projection is set, but the mask image is not provided or could not be loaded !" << endl);
2318  Exit(EXIT_FAILURE);
2319  }
2320  p_ProjectorManager->ProcessAndSetMask(p_ImageSpace->mp_maskImage);
2321  p_ImageSpace->DeallocateMaskImage();
2322  }
2323 
2324  p_ProjectorManager->SetVerbose(verbose_proj);
2325  // Check parameters
2326  if (p_ProjectorManager->CheckParameters())
2327  {
2328  Cerr("***** castor-recon() -> A problem occured while checking projector manager's parameters !" << endl);
2329  Exit(EXIT_FAILURE);
2330  }
2331  // Initialize projector manager
2332  if (p_ProjectorManager->Initialize())
2333  {
2334  Cerr("***** castor-recon() -> A problem occured while initializing projector manager !" << endl);
2335  Exit(EXIT_FAILURE);
2336  }
2337  // Check specific requirements for SPECT with attenuation correction
2338  if (p_ProjectorManager->CheckSPECTAttenuationCompatibility(path_to_attenuation_img))
2339  {
2340  Cerr("***** castor-recon() -> A problem occured while checking projector's compatibility with SPECT and attenuation correction !" << endl);
2341  Exit(EXIT_FAILURE);
2342  }
2343  // Verbose
2344  if (verbose_general>=5) Cout("----- Projector initialization OK -----" << endl);
2345 
2346  // ----------------------------------------------------------------------------------------
2347  // Create Optimizer Manager
2348  // ----------------------------------------------------------------------------------------
2349 
2350  // Verbose
2351  if (verbose_general>=5) Cout("----- Optimizer initialization ... -----" << endl);
2352  // Create object
2353  oOptimizerManager* p_OptimizerManager = new oOptimizerManager();
2354  // Set all parameters
2355  p_OptimizerManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2356  p_OptimizerManager->SetImageSpace(p_ImageSpace);
2357  p_OptimizerManager->SetDataMode(p_DataFile[0]->GetDataMode());
2358  p_OptimizerManager->SetDataType(p_DataFile[0]->GetDataType());
2359  p_OptimizerManager->SetDataSpec(p_DataFile[0]->GetDataSpec());
2360  p_OptimizerManager->SetOptionsOptimizer(options_optimizer);
2361  p_OptimizerManager->SetOptimizerFOMFlag(optimizer_fom);
2362  p_OptimizerManager->SetOptimizerImageStatFlag(optimizer_stat);
2363  p_OptimizerManager->SetOptionsPenalty(options_penalty);
2364  p_OptimizerManager->SetNbTOFBins(p_ProjectorManager->GetNbTOFBins());
2365  p_OptimizerManager->SetVerbose(verbose_opti);
2366  // Check parameters
2367  if (p_OptimizerManager->CheckParameters())
2368  {
2369  Cerr("***** castor-recon() -> A problem occured while checking optimizer manager's parameters !" << endl);
2370  Exit(EXIT_FAILURE);
2371  }
2372  // Initialize optimizer manager
2373  if (p_OptimizerManager->Initialize())
2374  {
2375  Cerr("***** castor-recon() -> A problem occured while initializing optimizer manager !" << endl);
2376  Exit(EXIT_FAILURE);
2377  }
2378  // Verbose
2379  if (verbose_general>=5) Cout("----- Optimizer initialization OK -----" << endl);
2380 
2381  // ----------------------------------------------------------------------------------------
2382  // Create Image Convolver Manager
2383  // ----------------------------------------------------------------------------------------
2384 
2385  // Verbose
2386  if (verbose_general>=5) Cout("----- Image Convolver initialization (if any) ... -----" << endl);
2387  // Create object
2388  oImageConvolverManager* p_ImageConvolverManager = new oImageConvolverManager();
2389  // Set all parameters
2390  p_ImageConvolverManager->SetVerbose(verbose_conv);
2391  p_ImageConvolverManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2392  p_ImageConvolverManager->SetOptions(options_image_convolver);
2393  // Check parameters
2394  if (p_ImageConvolverManager->CheckParameters())
2395  {
2396  Cerr("***** castor-recon() -> A problem occured while checking image convolver manager's parameters !" << endl);
2397  Exit(EXIT_FAILURE);
2398  }
2399  // Initialize image convolver manager
2400  if (p_ImageConvolverManager->Initialize())
2401  {
2402  Cerr("***** castor-recon() -> A problem occured while initializing image convolver manager !" << endl);
2403  Exit(EXIT_FAILURE);
2404  }
2405  // Verbose
2406  if (verbose_general>=5) Cout("----- Image Convolver initialization OK -----" << endl);
2407 
2408  // ----------------------------------------------------------------------------------------
2409  // Create Image Processing Manager
2410  // ----------------------------------------------------------------------------------------
2411 
2412  // Verbose
2413  if (verbose_general>=5) Cout("----- Image Processing initialization (if any) ... -----" << endl);
2414  // Create object
2415  oImageProcessingManager* p_ImageProcessingManager = new oImageProcessingManager();
2416  // Set all parameters
2417  p_ImageProcessingManager->SetVerbose(verbose_proc);
2418  p_ImageProcessingManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2419  p_ImageProcessingManager->SetOptions(options_image_processing);
2420  // Check parameters
2421  if (p_ImageProcessingManager->CheckParameters())
2422  {
2423  Cerr("***** castor-recon() -> A problem occured while checking image processing manager's parameters !" << endl);
2424  Exit(EXIT_FAILURE);
2425  }
2426  // Initialize image processing manager
2427  if (p_ImageProcessingManager->Initialize())
2428  {
2429  Cerr("***** castor-recon() -> A problem occured while initializing image processing manager !" << endl);
2430  Exit(EXIT_FAILURE);
2431  }
2432  // Verbose
2433  if (verbose_general>=5) Cout("----- Image Processing initialization OK -----" << endl);
2434 
2435  // ----------------------------------------------------------------------------------------
2436  // Create Dynamic Model Manager
2437  // ----------------------------------------------------------------------------------------
2438 
2439  // Verbose
2440  if (verbose_general>=5) Cout("----- Dynamic model initialization (if any) ... -----" << endl);
2441  // Create object
2442  oDynamicModelManager* p_DynamicModelManager = new oDynamicModelManager();
2443  // Set all parameters
2444  p_DynamicModelManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2445  p_DynamicModelManager->SetOptions(dynamic_model_options);
2446  p_DynamicModelManager->SetVerbose(verbose_dyna);
2447  // Check parameters
2448  if (p_DynamicModelManager->CheckParameters())
2449  {
2450  Cerr("***** castor-recon() -> A problem occured while checking dynamic model manager's parameters !" << endl);
2451  Exit(EXIT_FAILURE);
2452  }
2453  // Initialize optimizer manager
2454  if (p_DynamicModelManager->Initialize())
2455  {
2456  Cerr("***** castor-recon() -> A problem occured while initializing dynamic model manager !" << endl);
2457  Exit(EXIT_FAILURE);
2458  }
2459  // Verbose
2460  if (verbose_general>=5) Cout("----- Dynamic model initialization OK -----" << endl);
2461 
2462  // ----------------------------------------------------------------------------------------
2463  // Create Deformation Manager
2464  // ----------------------------------------------------------------------------------------
2465 
2466  // Verbose
2467  if (verbose_general>=5) Cout("----- Image deformation initialization (if any) ... -----" << endl);
2468  // Create object
2469  oDeformationManager* p_DeformationManager = new oDeformationManager();
2470  // Set all parameters
2471  p_DeformationManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2472  p_DeformationManager->SetDataMode(p_DataFile[0]->GetDataMode()); // required to know if sensitivity image deformation should be enabled
2473 
2474  if(resp_motion_options != "")
2475  {
2476  p_DeformationManager->SetOptions(resp_motion_options);
2477  p_DeformationManager->SetNbTransformations(nb_resp_gates-1);
2478  p_DeformationManager->SetMotionType(DEF_RESP_MOT);
2479  }
2480  else if (card_motion_options != "")
2481  {
2482  p_DeformationManager->SetOptions(card_motion_options);
2483  p_DeformationManager->SetNbTransformations(nb_card_gates-1);
2484  p_DeformationManager->SetMotionType(DEF_CARD_MOT);
2485  }
2486  else if (double_motion_options != "")
2487  {
2488  p_DeformationManager->SetOptions(double_motion_options);
2489  p_DeformationManager->SetNbTransformations(nb_resp_gates*nb_card_gates-1);
2490  p_DeformationManager->SetMotionType(DEF_DUAL_MOT);
2491  }
2492  else if (ipat_motion_options != "")
2493  {
2494  p_DeformationManager->SetOptions(ipat_motion_options);
2495  p_DeformationManager->SetNbTransformations(p_ImageDimensionsAndQuantification->GetNbIPatMotionSubsets() - 1);
2496  p_DeformationManager->SetMotionType(DEF_IPAT_MOT);
2497  }
2498  else
2499  p_DeformationManager->SetNbTransformations(0);
2500 
2501  p_DeformationManager->SetVerbose(verbose_defo);
2502  // Check parameters
2503  if (p_DeformationManager->CheckParameters())
2504  {
2505  Cerr("***** castor-recon() -> A problem occured while checking image deformation manager's parameters !" << endl);
2506  Exit(EXIT_FAILURE);
2507  }
2508  // Initialize optimizer manager
2509  if (p_DeformationManager->Initialize())
2510  {
2511  Cerr("***** castor-recon() -> A problem occured while initializing image deformation manager !" << endl);
2512  Exit(EXIT_FAILURE);
2513  }
2514  // Verbose
2515  if (verbose_general >=5) Cout("----- Image deformation initialization OK -----" << endl);
2516 
2517  // ============================================================================================================
2518  // Check for supported/implemented combinations of options, after creating managers/input data and before launching the algorithm
2519  // ============================================================================================================
2520 
2521  // ============================================================================================================
2522  // Algorithm initialization: create here sensitivity computation and launch algorithm
2523  // ============================================================================================================
2524 
2525  // --------------------------------------------------------------------------------------------
2526  // Create sensitivity if one of the following conditions is met
2527  // 1. Input file is a list-mode file and a sensitivity image was not provided
2528  // 2. Input file is a histogram and the -sens-histo option is provided (i.e. we explicitely
2529  // ask to compute the global sensitivity from the histogram datafile)
2530  // --------------------------------------------------------------------------------------------
2531  if ( (path_to_sensitivity_img.empty() && p_DataFile[0]->GetDataMode()==MODE_LIST) ||
2532  (p_DataFile[0]->GetDataMode()==MODE_HISTOGRAM && sensitivity_from_histogram)
2533  )
2534  {
2535  // Verbose
2536  if (verbose_general>=5) Cout("----- Image Sensitivity generation initialization ... -----" << endl);
2537  // Create object
2538  oSensitivityGenerator* p_Sensitivity = new oSensitivityGenerator();
2539  // Set parameters
2540  p_Sensitivity->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2541  p_Sensitivity->SetImageSpace(p_ImageSpace);
2542  p_Sensitivity->SetScanner(p_ScannerManager->GetScannerObject());
2543  p_Sensitivity->SetProjectorManager(p_ProjectorManager);
2544  p_Sensitivity->SetImageConvolverManager(p_ImageConvolverManager);
2545  p_Sensitivity->SetDeformationManager(p_DeformationManager);
2546  p_Sensitivity->SetGPUflag(gpu_flag);
2547  p_Sensitivity->SetPathToAttenuationImage(path_to_attenuation_img);
2548  p_Sensitivity->SetPathToMaskImage(path_to_mask_img);
2549  p_Sensitivity->SetNumberOfAtnGateImages(nb_atn_resp_imgs, nb_atn_card_imgs);
2550  p_Sensitivity->SetPathToNormalizationFileName(path_to_normalization_filename,invert_datafile_bed_order_flag);
2551  p_Sensitivity->SetDataFile(p_DataFile);
2552  p_Sensitivity->SetComputeFromHistogramFlag(sensitivity_from_histogram);
2553  p_Sensitivity->SetVerbose(verbose_sens);
2554  // Check parameters
2555  if (p_Sensitivity->CheckParameters())
2556  {
2557  Cerr("***** castor-recon() -> A problem occured while checking parameters of the sensitivity generator !" << endl);
2558  Exit(EXIT_FAILURE);
2559  }
2560  // Initialize the sensitivity generator
2561  if (p_Sensitivity->Initialize())
2562  {
2563  Cerr("***** castor-recon() -> A problem occured while initializing the sensitivity generator !" << endl);
2564  Exit(EXIT_FAILURE);
2565  }
2566  // Set the sensitivity mode ON for the projector manager
2567  p_ProjectorManager->SetSensitivityModeOn();
2568  // Launch the computation
2569  if (p_Sensitivity->Launch())
2570  {
2571  Cerr("***** castor-recon() -> A problem occured while computing the sensitivity !" << endl);
2572  Exit(EXIT_FAILURE);
2573  }
2574  // Set the sensitivity mode OFF for the projector manager
2575  p_ProjectorManager->SetSensitivityModeOff();
2576  // Get the path to the sensitivity image (will be given to the algorithm if input file is a list-mode)
2577  if (p_DataFile[0]->GetDataMode()==MODE_LIST) path_to_sensitivity_img = p_Sensitivity->GetPathToSensitivityImage();
2578  else path_to_sensitivity_img = "";
2579  // Delete the generator
2580  delete p_Sensitivity;
2581  // Exit now if asked for
2582  if (exit_after_sensitivity)
2583  {
2584  if (verbose_general>=VERBOSE_LIGHT) Cout("castor-recon() -> Asked to exit after sensitivity computation." << endl);
2585  // Delete objects in the inverse order in which they were created
2586  delete p_ImageSpace;
2587  delete p_DeformationManager;
2588  delete p_DynamicModelManager;
2589  delete p_ImageProcessingManager;
2590  delete p_ImageConvolverManager;
2591  delete p_OptimizerManager;
2592  delete p_ProjectorManager;
2593  for (int i=0 ; i<nb_beds ; i++) delete p_DataFile[i];
2594  delete[] p_DataFile;
2595  delete p_ImageDimensionsAndQuantification;
2596  // And exit
2597  return 0;
2598  }
2599  }
2600 
2601  // ----------------------------------------------------------------------------------------
2602  // Create algorithm
2603  // ----------------------------------------------------------------------------------------
2604 
2605  // If the number of events in a datafile is below the number of threads for projections, then we must reduce this number of threads, otherwise
2606  // the datafile reading cannot work (and it has no sense anyway). This is done by the following function
2607  p_ImageDimensionsAndQuantification->CheckNumberOfProjectionThreadsConsistencyWithDataFileSize(p_DataFile);
2608  // Verbose
2609  if (verbose_general>=5) Cout("----- Iterative reconstruction algorithm initialization ... -----" << endl);
2610  // Create object
2611  oIterativeAlgorithm* p_Algorithm = new oIterativeAlgorithm();
2612  // Set parameters
2613  p_Algorithm->SetImageConvolverManager(p_ImageConvolverManager);
2614  p_Algorithm->SetImageProcessingManager(p_ImageProcessingManager);
2615  p_Algorithm->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
2616  p_Algorithm->SetImageSpace(p_ImageSpace);
2617  p_Algorithm->SetProjectorManager(p_ProjectorManager);
2618  p_Algorithm->SetDynamicModelManager(p_DynamicModelManager);
2619  p_Algorithm->SetDeformationManager(p_DeformationManager);
2620  p_Algorithm->SetOptimizerManager(p_OptimizerManager);
2621  p_Algorithm->SetDataFile(p_DataFile);
2622  p_Algorithm->SetGPUflag(gpu_flag);
2623  p_Algorithm->SetVerbose(verbose_algo);
2624  p_Algorithm->SetNbBeds(nb_beds);
2625  p_Algorithm->SetPathInitImage(path_to_initial_img);
2626  p_Algorithm->SetPathToAttenuationImage(path_to_attenuation_img);
2627  p_Algorithm->SetPathToSensitivityImage(path_to_sensitivity_img);
2628  p_Algorithm->SetPathToMultiModalImage(path_to_multimodal_img);
2629  p_Algorithm->SetPathToMaskImage(path_to_mask_img);
2630  p_Algorithm->SetSaveSensitivityHistoFlag(save_sens_histo);
2631  p_Algorithm->SetSaveSubsetImageFlag(save_subset_image);
2632  if (p_Algorithm->SetNbIterationsAndSubsets(nb_iterations_subsets))
2633  {
2634  Cerr("***** castor-recon() -> Error while setting the numbers of iterations and subsets !" << endl);
2635  Exit(EXIT_FAILURE);
2636  }
2637  if (p_Algorithm->SetOutputIterations(output_iterations))
2638  {
2639  Cerr("***** castor-recon() -> Error while setting the selected output iterations !" << endl);
2640  Exit(EXIT_FAILURE);
2641  }
2642  if (p_Algorithm->Iterate())
2643  {
2644  Cerr("***** castor-recon() -> Error while performing the reconstruction" << endl);
2645  Exit(EXIT_FAILURE);
2646  }
2647 
2648  // Display profiling results
2649  p_ChronoManager->Display();
2650 
2651  // ============================================================================================================
2652  // End
2653  // ============================================================================================================
2654 
2655  // Delete objects in the inverse order in which they were created
2656  if (verbose_general>=5) Cout("----- Deleting CASToR objects ... -----" << endl);
2657  delete p_Algorithm;
2658  delete p_ImageSpace;
2659  delete p_DeformationManager;
2660  delete p_DynamicModelManager;
2661  delete p_ImageProcessingManager;
2662  delete p_ImageConvolverManager;
2663  delete p_OptimizerManager;
2664  delete p_ProjectorManager;
2665  for (int i=0 ; i<nb_beds ; i++) delete p_DataFile[i];
2666  delete[] p_DataFile;
2667  delete p_ImageDimensionsAndQuantification;
2668  if (verbose_general>=5) Cout("----- CASToR objects successfully deleted -----" << endl);
2669 
2670  // Ending
2671  if (verbose_general>=1) Cout(endl);
2672  #ifdef CASTOR_MPI
2673  MPI_Finalize();
2674  #endif
2675  return EXIT_SUCCESS;
2676 }
void SetNbThreads(int a_nbThreadsForProjection, int a_nbThreadsForImageComputation)
Set the number of threads for both projection and image computations.
void ShowHelpProj()
Display command line options related to the Projector module for castor-recon.
void SetIgnorePOIFlag(bool a_ignorePOIFlag)
Set a boolean that that if we ignore POI information or not.
Definition: vDataFile.hh:468
void SetRespBasisFunctionsFile(const string &a_respBasisFunctionsFile)
Set the file name containing the respiratory basis functions coefficients.
This class is designed to be a mother virtual class for DataFile.
Definition: vDataFile.hh:103
This header file is mainly used to declare some macro definitions and all includes needed from the st...
static sScannerManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
Set the member mp_ImageDimensionsAndQuantification to the provided value.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
Set the Image Dimensions and Quantification Object.
static sRandomNumberGenerator * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
Declaration of class oImageDimensionsAndQuantification.
void SetPathToAttenuationImage(string a_pathToAttenuationImage)
This function is used to set the path to the attenuation image.
int InitDynamicData(string a_pathTo4DDataSplittingFile, int a_respMotionCorrectionFlag, int a_cardMotionCorrectionFlag, int a_doubleMotionCorrectionFlag, int a_invMotionCorrectionFlag, int a_nbRespGates, int a_nbCardGates)
Call the eponym function from the oDynamicDataManager object in order to initialize its data...
#define MODE_HISTOGRAM
Definition: vDataFile.hh:59
void SetNumberOfAtnGateImages(int a_nbAtnRespGateImages, int a_nbAtnCardGateImages)
int Iterate()
Just call either the IterateCPU or the IterateGPU function as asked for.
int SetNbIterationsAndSubsets(const string &a_nbIterationsSubsets)
Set the number of iterations and subsets.
void SetVerbose(int a_verbose)
set verbosity
void SetGPUflag(bool a_flagGPU)
Set the GPU flag.
void SetFOVSizeZ(FLTNB a_fovSizeZ)
Set the FOV's size along the Z axis, in mm.
#define MODE_LIST
Definition: vDataFile.hh:57
static void ShowCommonHelp()
This function is used to print out some help about the use of options common to all projectors...
Definition: vProjector.cc:112
void SetPathToMultiModalImage(vector< string > a_pathToMultiModalImage)
Set path to multimodal images.
void SetIgnoredCorrections(const string &a_ignoredCorrectionsList)
Set the string specifying the corrections that will be ignored.
#define FLTNB
Definition: gVariables.hh:81
int InitMaskImage(const string &a_pathToImage)
Memory allocation and initialization for the mask image.
Definition: oImageSpace.cc:675
void SetOptimizerImageStatFlag(bool a_optimizerImageStatFlag)
Set the optimizer image stat flag that specifies if some basic statistics about image update is compu...
void ShowHelpDynamicModel()
Show help about all implemented dynamic models.
void ShowHelpDeformation()
Show help about all implemented deformations.
void ShowHelp()
Display main command line options for castor-recon.
Definition: castor-recon.cc:64
void SetComputationStrategy(int a_computationStrategy)
Set the computation strategy for the system matrix elements storage.
void SetProjectorManager(oProjectorManager *ap_ProjectorManager)
Set the Projector Manager Object.
int Initialize()
A function used to initialize the manager and the projectors or system matrices it manages...
static void ShowCommonHelp()
This function does not take any parameter and is used to display some help about the syntax of the op...
void SetSaveSensitivityHistoFlag(bool a_saveSensitivityHistoFlag)
Set the flag that specifies if the sensitivity image in histogram mode has to be saved for each subse...
int Initialize()
A public function used to initialize the sensitivity generator.
void SetVerbose(int a_verboseLevel)
Set the member m_verboseLevel to the provided value.
void ShowHelpImgp()
Display command line options related to the Image Processing module for castor-recon.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
Set the image dimensions in use.
int FindScannerSystem(string a_scannerName)
Look for a file matching with the scanner name in parameter inside the scanner repository.
void SetDeformationManager(oDeformationManager *ap_DeformationManager)
Set the Deformation Manager Object.
void SetFOVSizeY(FLTNB a_fovSizeY)
Set the FOV's size along the Y axis, in mm.
void SetNbVoxZ(INTNB a_nbVoxZ)
Set the number of voxels along the Z axis.
void SetNbTimeBasisFunctions(int a_nbTimeBasisFunctions)
Set the number of time basis functions.
void SetOptions(const string &a_options)
Set the respiratory motion options contained in the provided string.
int CheckParameters()
Check validity of all parameters.
void SetVerbose(int a_verboseLevel)
Set the verbose level.
int BuildScannerObject()
Instantiate the specific scanner object related to the modality, and set verbosity of scanner object...
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
Set the image dimensions in use.
string GetPathToSensitivityImage()
This function return the path to the sensitivity image.
void SetNbRespGates(int a_nbRespGates)
Set the number of respiratory gates.
void ShowHelpComputation()
Display command line options related to the computation settings for castor-recon.
int CheckParameters()
A function used to check the parameters settings.
void SetPathToSensitivityImage(string a_pathToSensitivityImage)
Set path to the sensitivity image.
This class is designed to manage the optimization part of an iterative reconstruction.
void SetOffsetY(FLTNB a_offsetY)
Set the image offset along the Y axis, in mm.
void CheckNumberOfProjectionThreadsConsistencyWithDataFileSize(vDataFile **a2p_DataFile)
void SetFOVOutMasking(FLTNB a_fovOutPercent, INTNB a_nbSliceOutMask)
Set the output FOV masking settings: transaxial unmasked FOV percent and number of extrem slices to r...
void ShowHelpImageProcessingModule()
Show help about all implemented image processing modules.
void SetVerbose(int a_verboseLevel)
set verbosity
Definition: vDataFile.hh:439
void ShowHelpMiscellaneous()
Display command line options related to miscellaneous settings for castor-recon.
int DealWithBedPositions(vDataFile **a2p_DataFile)
Deal with provided or default bed relative positions.
void SetDeformationManager(oDeformationManager *ap_DeformationManager)
This function is used to set the pointer to the oDeformationManager in use.
void SetMotionType(int a_motionType)
Set the nature of motion correction (Deformation type macro)
void SetOptionsForward(const string &a_optionsForward)
Set the forward projection options contained in the provided string.
void SetDataMode(int a_dataMode)
Set the mode of reconstruction.
void SetSensitivityModeOff()
Say that the projector will no longer be used to compute the global sensitivity.
int SetOutputIterations(const string &a_outputIterations)
Set the selected output iterations.
void SetVerbose(int a_verboseLevel)
Set Verbosity.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
Set the member mp_ImageDimensionsAndQuantification to the provided value.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
Set the image dimensions in use.
This is the main class for iterative reconstructions, that manages the iteration loops. This class manages an iterative reconstruction of any kind, using a vDataFile, and through the use of an oProjector, an oOptimizer, a oConvolver, a oImageSpace.
int Launch()
A public function used to launch the sensitivity generator (compute the sensitivity image) ...
int Initialize()
A function used to initialize the manager and all image processing modules it manages.
void SetDataSpec(int a_dataSpec)
Set the specificity of the data: EMISSION or TRANSMISSION.
string GetPathToScannerFile()
void SetVerbose(int a_verboseLevel)
Set verbosity level.
void ShowHelpPenalty()
Show help about all implemented penalties.
int Initialize()
Set the dynamic model flag and instanciate/initialize model objects through the ParseOptionsAndInitia...
static sOutputManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
void SetOffsetZ(FLTNB a_offsetZ)
Set the image offset along the Z axis, in mm.
#define SCANNER_PET
Declaration of class iDataFileCT.
void SetVerbose(int a_verboseLevel)
void SetImageConvolverManager(oImageConvolverManager *ap_ImageConvolverManager)
This function is used to set the pointer to the oImageConvolverManager in use.
void SetSaveLUTFlag(bool a_flag)
Set to on the flag indicating a LUT generated by a geom file should be written on disk or not...
Declaration of class iDataFilePET.
void SetOptions(vector< string > a_options)
Set the member m_options to the provided value.
Declaration of class oIterativeAlgorithm.
void SetOptionsCommon(const string &a_optionsCommon)
Set the common projection options contained in the provided string.
void ShowHelpDynamic()
Display command line options related to the dynamic features for castor-recon.
int Initialize()
A function used to initialize the manager and the optimizer it manages.
int CheckParameters()
Check if all parameters have been correctly initialized, and call the CheckParameters function of the...
void SetNbVoxX(INTNB a_nbVoxX)
Set the number of voxels along the X axis.
int SetNbThreads(const string &a_nbThreads)
Set the number of threads.
uint16_t nb_resp_gates
Declaration of class iDataFileSPECT.
void Exit(int code)
void SetOptionsBackward(const string &a_optionsBackward)
Set the backward projection options contained in the provided string.
void DeallocateMaskImage()
Free memory for the mask image.
Definition: oImageSpace.cc:716
uint16_t nb_card_gates
Declaration of class iScannerPET.
int IntfReadHeader(const string &a_pathToHeaderFile, Intf_fields *ap_IntfFields, int vb)
Read an Interfile header.
void SetImageSpace(oImageSpace *ap_ImageSpace)
Set the image space in use.
int SetOutNbPrec(string a_format)
Set the output format and precision used for numeric display.
int InstantiateScanner()
Instantiate scanner using the related function in the scanner classes.
int LogCommandLine(int argc, char **argv)
Write log file header with the provided command line options and different informations.
void SetPathToNormalizationFileName(vector< string > ap_pathToNormalizationFileName, bool a_inverseDataFileOrderFlag)
This function is used to set the path to the normalization file names, and the flag saying if their o...
static sAddonManager * GetInstance()
void SetSaveSubsetImageFlag(bool a_saveImageAfterSubsets)
Set the flag that specifies if the image has to be saved for each subset.
void SetComputeFromHistogramFlag(bool a_computeFromHistogramFlag)
This function is used to set the m_computeFromHistogramFlag.
int ConvertFromString(const string &a_str, string *a_result)
Copy the 'a_str' string in the position pointed by 'a_result'.
Definition: gOptions.cc:771
#define DEF_DUAL_MOT
#define FIXED_LIST_COMPUTATION_STRATEGY
void SetNbVoxY(INTNB a_nbVoxY)
Set the number of voxels along the Y axis.
void SetNbRespBasisFunctions(int a_nbRespBasisFunctions)
Set the number of respiratory basis functions.
int CheckConfigDir(const string &a_path)
Set the path to the CASTOR config directory from the given path if not empty or through the existence...
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
Set the image dimensions in use.
void SetVerbose(int a_verboseLevel)
Set the verbose level.
int CheckParameters()
A function used to check the parameters settings.
void SetNbCardGates(int a_nbCardGates)
Set the number of cardiac gates.
void SetDataMode(int a_dataMode)
Set the mode of the data (histogram, list-mode)
#define Cerr(MESSAGE)
#define SCANNER_SPECT_CONVERGENT
int BuildLUT()
Call the eponym function of the scanner class.
void SetDataType(int a_dataType)
Set the type of the data (pet, spect, etc)
int GetNbIPatMotionSubsets()
call the eponym function from the oDynamicDataManager object
int Initialize()
A function used to initialize all that is needed.
Singleton class that manages output writing on disk (images, sinograms, etc). It also manages loggi...
void SetMPIRankAndSize(int a_mpiRank, int a_mpiSize)
Set the MPI rank of the MPI instance, and the MPI size (the number of instances)
#define DEF_RESP_MOT
void SetOptionsOptimizer(const string &a_optionsOptimizer)
Set the optimizer projection options contained in the provided string.
This class is designed to manage the use of dynamic model in the reconstruction.
void SetFrames(const string &a_frameList)
Set the frame list (a string that will be parsed by the InitializeFramingAndQuantification function) ...
void SetFOVSizeX(FLTNB a_fovSizeX)
Set the FOV's size along the X axis, in mm.
void SetNbBeds(int a_nbBeds)
Set number of beds (bed positions)
void ShowHelpOutput()
Display command line options related to output settings for castor-recon.
#define VERBOSE_LIGHT
int ReadDataASCIIFile(const string &a_file, const string &a_keyword, T *ap_return, int a_nbElts, bool a_mandatoryFlag)
Look for "a_nbElts" elts in the "a_file" file matching the "a_keyword" string passed as parameter a...
Definition: gOptions.cc:123
int GetNbThreadsMax()
Get the maximum between the number of threads used for projections and image operations.
Singleton class that Instantiate and initialize the scanner object.
int Initialize(int a_nbThreads, int a_nbExtra)
Instantiate pseudo random number generators, one per thread by default, and additional extra ones if ...
void SetImageProcessingManager(oImageProcessingManager *ap_ImageProcessingManager)
Set the Image Processing Manager Object.
void SetOptions(vector< string > a_options)
Set the member m_options to the provided value.
void SetOptionsPenalty(const string &a_optionsPenalty)
Set the penalty projection options contained in the provided string.
int CheckParameters()
A function used to check the parameters settings.
Declaration of class sScannerManager.
void SetPathInitImage(string a_pathToInitialImage)
Set path to an initial image.
void SetVoxSizeY(FLTNB a_voxSizeY)
Set the voxel's size along the Y axis, in mm.
void SetPathToMaskImage(string a_pathToMaskImage)
Set path to a mask image.
Inherit from vDataFile. Class that manages the reading of a SPECT input file (header + data)...
This class is designed to manage the different image convolvers and to apply them.
void SetOffsetX(FLTNB a_offsetX)
Set the image offset along the X axis, in mm.
void SetHeaderDataFileName(const string &a_headerFileName)
set the data header file name
Definition: vDataFile.hh:476
FLTNB * mp_maskImage
Definition: oImageSpace.hh:121
vScanner * GetScannerObject()
void SetVoxSizeX(FLTNB a_voxSizeX)
Set the voxel's size along the X axis, in mm.
void IntfKeyInitFields(Intf_fields *ap_IF)
Init the file of an Interfile fields structure passed in parameter to their default values...
int CheckParameters()
A function used to check the parameters settings.
void SetSensitivityModeOn()
Say that the projector will be used to compute the global sensitivity.
void SetMergeDynImagesFlag(bool a_flag)
Set to on the flag indicating that a dynamic serie of 3D images should be written on disk in one file...
void SetPathToAttenuationImage(string a_pathToAttenuationImage)
This function is used to set the path to the attenuation image.
void SetMPIRank(int a_mpiRank)
Initialize the machine index for MPI.
void SetVerbose(int a_verboseLevel)
Set the verbose level.
void SetTimeBasisFunctionsFile(const string &a_timeBasisFunctionsFile)
Set the file name containing the time basis functions coefficients.
void SetDataFile(vDataFile **a2p_DataFile)
This function is used to set the pointer to the vDataFile array in use.
void SetDataFile(vDataFile *ap_DataFile)
Set a data file in use to later recover some information from it.
#define DEF_CARD_MOT
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
set the pointer to the oImageDimensionsAndQuantification object
Definition: oImageSpace.hh:624
int main(int argc, char **argv)
void ShowHelpProjector()
Show help about all implemented projectors.
#define KEYWORD_MANDATORY
Definition: gOptions.hh:48
int CheckParameters()
A function used to check the parameters settings.
#define KEYWORD_OPTIONAL
Definition: gOptions.hh:50
static void ShowCommonHelp()
This function does not take any parameter and is used to display some help about the syntax of the op...
void SetImageSpace(oImageSpace *ap_ImageSpace)
This function is used to set the pointer to the oImageSpace in use.
void SetPathToMaskImage(string a_pathToMaskImage)
Set path to a mask image.
int CheckParameters()
This function is used to check parameters after the latter have been all set using Set functions...
Singleton class that generate a thread-safe random generator number for openMP As singleton...
int ProcessAndSetMask(FLTNB *ap_maskImage)
Process and set the provided mask image for projector masking.
int CheckParameters()
This function is used to check parameters after the latter have been all set using Set functions...
int Initialize()
Initialization : .
Declaration of class sRandomNumberGenerator.
void SetNbTOFBins(int a_nbTOFBins)
Set the number of TOF bins in use.
void SetNbMultiModalImages(int a_nbMultiModalImages)
Set the number of additional multimodal images.
void SetScanner(vScanner *ap_Scanner)
This function is used to set the pointer to the vScanner in use.
void SetVerbose(int a_verbose)
Set the member m_verboseLevel to the provided value.
#define INTNB
Definition: gVariables.hh:92
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ID)
Set the pointer to the image dimensions and quantification object.
void SetVerbose(int a_verboseLevel)
Set the verbose level.
int Initialize()
Initialize all thread-safe buffers for profiling.
This class is designed to manage the image-based deformation part of the reconstruction.
void SetVoxSizeZ(FLTNB a_voxSizeZ)
Set the voxel's size along the Z axis, in mm.
int CheckDynamicParameters(int64_t a_nbEvents)
Call the eponym function from the oDynamicDataManager object in order to check its parameters...
void SetVerbose(int a_verboseLevel)
set verbosity
Definition: oImageSpace.hh:617
Declaration of class oSensitivityGenerator.
This class is designed to manage the different image processing modules and to apply them...
This class is designed to manage the projection part of the reconstruction.
int GetNbTOFBins()
Get the number of TOF bins associated to the projector.
Interfile fields. This structure contains all the Interfile keys currently managed by CASToR Decl...
void SetBedIndex(int a_bedIndex)
set the bed index corresponding to this data file
Definition: vDataFile.hh:420
#define DEF_IPAT_MOT
int Initialize()
Set the flags for the different motion types and instanciate/initialize deformation objects through t...
Declaration of class sOutputManager.
This class holds all the matrices in the image domain that can be used in the algorithm: image...
Definition: oImageSpace.hh:61
void SetImageSpace(oImageSpace *ap_ImageSpace)
Set the Image Space Object.
void SetImageConvolverManager(oImageConvolverManager *ap_ImageConvolverManager)
Set the Image Convolver Manager Object.
int GetNbThreadsForImageComputation()
Get the number of threads used for image operations.
int GetDataMode()
Definition: vDataFile.hh:300
int InitOutputDirectory(const string &a_pathFout, const string &a_pathDout)
Create the output directory if any, extract the base name and create the log file.
int CheckParameters()
A public function used to check the parameters settings.
Inherit from vDataFile. Class that manages the reading of a CT input file (header + data)...
Definition: iDataFileCT.hh:48
This class is designed to manage all dimensions and quantification related stuff. ...
void ShowHelpInput()
Display command line options related to input settings for castor-recon.
int GetNbThreadsForProjection()
Get the number of threads used for projections.
void SetNbBeds(int a_nbBeds)
Set the number of bed positions and allocate the bed positions if not already done.
void SetNbCardBasisFunctions(int a_nbCardBasisFunctions)
Set the number of cardiac basis functions.
This file is used for all kind of different functions designed for options parsing and ASCII file rea...
Declaration of class sChronoManager.
void SetProjectorManager(oProjectorManager *ap_ProjectorManager)
This function is used to set the pointer to the oProjectorManager in use.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
This function is used to set the pointer to the oImageDimensionsAndQuantification in use...
void ShowHelpDimensions()
Display command line options related to image dimensions for castor-recon.
void SetGPUflag(bool a_flagGPU)
This function is used to set the GPU flag; do we use GPU or not.
void SetVerbose(int a_verboseLevel)
Set the member m_verboseLevel to the provided value.
void SetVerbose(int a_verboseLevel)
set verbosity
#define Cout(MESSAGE)
int CheckSPECTAttenuationCompatibility(const string &a_pathToAttenuationImage)
A function used to check specific compatibility with SPECT and attenuation correction.
This class is designed to manage some profiling of the code.
#define CASTOR_VERSION
Definition: gVariables.hh:70
void SetDataFileName(vector< string > ap_dataFileName)
void SetCardBasisFunctionsFile(const string &a_cardBasisFunctionsFile)
Set the file name containing the cardiac basis functions coefficients.
void SetScanner(vScanner *ap_Scanner)
Set the scanner in use.
#define SCANNER_CT
void ShowHelpAlgo()
Display command line options related to the Optimization algorithm module for castor-recon.
void ShowHelpOptimizer()
Show help about all implemented optimizers.
void SetOptimizerManager(oOptimizerManager *ap_OptimizerManager)
Set the Optimizer Manager Object.
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
set the pointer to the oImageDimensionsAndQuantification object
Definition: vDataFile.hh:453
Inherit from vDataFile. Class that manages the reading of a PET input file (header + data)...
Definition: iDataFilePET.hh:44
int GetGeometricInfoFromDataFile(string a_pathToDataFilename)
Call the specialized function of the scanner object in order to get geometric informations from the d...
void SetNbTransformations(int a_nbTransformations)
Set the total number of transformations/deformations.
void SetVerbose(int a_verbose)
Set the verbose level.
int ReadStringOption(const string &a_input, T *ap_return, int a_nbElts, const string &sep, const string &a_option)
Parse the 'a_input' string corresponding to the 'a_option' into 'a_nbElts' elements, using the 'sep' separator. The results are returned in the templated 'ap_return' dynamic templated array. Call "ConvertFromString()" to perform the correct conversion depending on the type of the data to convert.
Definition: gOptions.cc:62
void SetDataFile(vDataFile **a2p_DataFile)
Set the list of DataFile.
void SetOptions(const string &a_options)
Set the motion options contained in the provided string, and the related number of gates...
static sChronoManager * GetInstance()
Instantiate the singleton if not already done, then return the pointer to its instance.
void SetDynamicModelManager(oDynamicModelManager *ap_DynamicModelManager)
Set the Dynamic Model Manager Object.
This class is designed to manage the computation of the sensitivity image.
Declaration of class sAddonManager.
void GetUserEndianness()
Check user/host computer endianness and write it to the global variable User_Endianness.
int SetFlipOut(const string &a_flipOut)
Set the output flip options, the parameter being a string potentially containing the letters x...
int Initialize()
A function used to initialize the manager and all image convolvers it manages.
#define KEYWORD_OPTIONAL_ERROR
Definition: gOptions.hh:56
int SetDynamicSpecificQuantificationFactors(const string &a_quantificationFile)
Apply specific quantification factors manually provided as an option.
void ShowHelpCorrection()
Display command line options related to correction settings for castor-recon.
void ShowHelpImageConvolver()
Show help about all implemented image convolvers.
void Display()
Display the results of the duration buffers.
void SetOptimizerFOMFlag(bool a_optimizerFOMFlag)
Set the optimizer FOM flag that specifies if some figures-of-merit (FOM) will be computed in the data...