CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
oImageDimensionsAndQuantification.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 
32 #include "vDataFile.hh"
33 
34 // =====================================================================
35 // ---------------------------------------------------------------------
36 // ---------------------------------------------------------------------
37 // =====================================================================
38 
40 {
41  // Set all members to default values
42  m_nbVoxX = -1;
43  m_nbVoxY = -1;
44  m_nbVoxZ = -1;
45  m_nbVoxXY = -1;
46  m_nbVoxXYZ = -1;
47  m_voxSizeX = -1.;
48  m_voxSizeY = -1.;
49  m_voxSizeZ = -1.;
50  m_fovSizeX = -1.;
51  m_fovSizeY = -1.;
52  m_fovSizeZ = -1.;
53  m_fovOutPercent = 0.;
54  m_flipOutX = false;
55  m_flipOutY = false;
56  m_flipOutZ = false;
57  m_offsetX = 0.;
58  m_offsetY = 0.;
59  m_offsetZ = 0.;
60  m_frameList = "";
61  m_nbTimeFrames = 1;
66  m2p_frameTimeStopInMs = NULL;
68  m_timeStaticFlag = false;
79  m_nbRespGates = 1;
83  m_respStaticFlag = false;
84  m_nbCardGates = 1;
88  m_cardStaticFlag = false;
91  m_mpiRank = 0;
92  m_mpiSize = 1;
93  m_nbBeds = -1;
94  mp_bedPositions = NULL;
95  m_providedBedPosition = false;
96  m_verbose = 0;
97  m_checked = false;
98  m_initialized = false;
99  // Allocate Dynamic data manager object
102 }
103 
104 // =====================================================================
105 // ---------------------------------------------------------------------
106 // ---------------------------------------------------------------------
107 // =====================================================================
108 
110 {
111  // Free frame duration
113  {
114  for (int bed=0; bed<m_nbBeds; bed++) if (m2p_frameDurationsInMs[bed]) delete m2p_frameDurationsInMs[bed];
115  delete m2p_frameDurationsInMs;
116  }
117  // Free frame time start
119  {
120  for (int bed=0; bed<m_nbBeds; bed++) if (m2p_frameTimeStartInMs[bed]) delete m2p_frameTimeStartInMs[bed];
121  delete m2p_frameTimeStartInMs;
122  }
123  // Free frame time stop
125  {
126  for (int bed=0; bed<m_nbBeds; bed++) if (m2p_frameTimeStopInMs[bed]) delete m2p_frameTimeStopInMs[bed];
127  delete m2p_frameTimeStopInMs;
128  }
129  // Free quantification factors
131  {
132  for (int bed=0; bed<m_nbBeds; bed++)
133  {
134  if (m3p_quantificationFactors[bed])
135  {
136  for (int fr=0; fr<m_nbTimeFrames; fr++)
137  if (m3p_quantificationFactors[bed][fr]) delete m3p_quantificationFactors[bed][fr];
138  }
139  delete m3p_quantificationFactors[bed];
140  }
142  }
143  // Free the bed positions
144  if (mp_bedPositions) free(mp_bedPositions);
145  // Delete the dynamic data manager
146  delete mp_DynamicDataManager;
147 }
148 
149 // =====================================================================
150 // ---------------------------------------------------------------------
151 // ---------------------------------------------------------------------
152 // =====================================================================
153 
155 {
156  // If empty, then put all flags to false and return
157  if (a_flipOut=="")
158  {
159  m_flipOutX = false;
160  m_flipOutY = false;
161  m_flipOutZ = false;
162  return 0;
163  }
164  // Test all possible settings !
165  if (a_flipOut=="x" || a_flipOut=="X")
166  {
167  m_flipOutX = true;
168  m_flipOutY = false;
169  m_flipOutZ = false;
170  }
171  else if (a_flipOut=="y" || a_flipOut=="Y")
172  {
173  m_flipOutX = false;
174  m_flipOutY = true;
175  m_flipOutZ = false;
176  }
177  else if (a_flipOut=="z" || a_flipOut=="Z")
178  {
179  m_flipOutX = false;
180  m_flipOutY = false;
181  m_flipOutZ = true;
182  }
183  else if (a_flipOut=="xy" || a_flipOut=="yx" || a_flipOut=="XY" || a_flipOut=="YX")
184  {
185  m_flipOutX = true;
186  m_flipOutY = true;
187  m_flipOutZ = false;
188  }
189  else if (a_flipOut=="zy" || a_flipOut=="yz" || a_flipOut=="ZY" || a_flipOut=="YZ")
190  {
191  m_flipOutX = false;
192  m_flipOutY = true;
193  m_flipOutZ = true;
194  }
195  else if (a_flipOut=="xz" || a_flipOut=="zx" || a_flipOut=="XZ" || a_flipOut=="ZX")
196  {
197  m_flipOutX = true;
198  m_flipOutY = false;
199  m_flipOutZ = true;
200  }
201  else if ( a_flipOut=="xyz" || a_flipOut=="xzy" || a_flipOut=="yxz" || a_flipOut=="yzx" || a_flipOut=="zxy" || a_flipOut=="zyx" ||
202  a_flipOut=="XYZ" || a_flipOut=="XZY" || a_flipOut=="YXZ" || a_flipOut=="YZX" || a_flipOut=="ZXY" || a_flipOut=="ZYX" )
203  {
204  m_flipOutX = true;
205  m_flipOutY = true;
206  m_flipOutZ = true;
207  }
208  // Otherwise, throw an error
209  else
210  {
211  Cerr("***** oImageDimensionsAndQuantification::SetFlipOut() -> Output flip settings is incorrect !" << endl);
212  return 1;
213  }
214  // Normal end
215  return 0;
216 }
217 
218 // =====================================================================
219 // ---------------------------------------------------------------------
220 // ---------------------------------------------------------------------
221 // =====================================================================
222 
224 {
225  // The number of threads can be given as two separated numbers to distinguish between the number of threads for projection computation
226  // and the number of threads for image computation. Be careful that the thread-safe images are allocated with respect to the number
227  // of projection threads; the number of threads used for image computation are strictly used for computation over images (i.e. when
228  // using a loop over voxels.
229 
230  // Look for a comma and check that there is only one comma without missing parameters
231  size_t first_comma = a_nbThreads.find_first_of(",");
232  size_t last_comma = a_nbThreads.find_last_of(",");
233  if (first_comma!=last_comma || first_comma==0 || first_comma==a_nbThreads.length()-1)
234  {
235  Cerr("***** oImageDimensionsAndQuantification::SetNbThreads() -> Wrong syntax in the thread parameters ! See help." << endl);
236  return 1;
237  }
238 
239  // Case for a single parameter
240  if (first_comma==string::npos)
241  {
242  m_nbThreadsForProjection = atoi( a_nbThreads.c_str() );
244  }
245  // Case for two parameters
246  else
247  {
248  m_nbThreadsForProjection = atoi( (a_nbThreads.substr(0,first_comma)).c_str() );
249  m_nbThreadsForImageComputation = atoi( (a_nbThreads.substr(first_comma+1)).c_str() );
250  }
251 
252  // Checks for negative threads numbers
254  {
255  Cerr("***** oImageDimensionsAndQuantification::SetNbThreads() -> Negative number of threads provided for projection computation !" << endl);
256  return 1;
257  }
259  {
260  Cerr("***** oImageDimensionsAndQuantification::SetNbThreads() -> Negative number of threads provided for image computation !" << endl);
261  return 1;
262  }
263 
264  #ifdef CASTOR_OMP
265  // If number of threads is 0, we automatically determine the maximum number of threads using OMP functions
266  if (m_nbThreadsForProjection==0) m_nbThreadsForProjection = omp_get_max_threads();
267  if (m_nbThreadsForImageComputation==0) m_nbThreadsForImageComputation = omp_get_max_threads();
268  // At this step, set by default the number of threads for image operations
269  omp_set_num_threads(m_nbThreadsForImageComputation);
270  #endif
271 
272  // End
273  return 0;
274 }
275 
276 // =====================================================================
277 // ---------------------------------------------------------------------
278 // ---------------------------------------------------------------------
279 // =====================================================================
280 
282 {
283  // Number of threads
286  // MPI stuff
287  m_mpiSize = 1;
288  m_mpiRank = 0;
289  // Unlock
290  m_checked = true;
291 }
292 
293 
294 // =====================================================================
295 // ---------------------------------------------------------------------
296 // ---------------------------------------------------------------------
297 // =====================================================================
298 
300 {
301  // Number of threads
303  {
304  Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Should provide a strictly positive number of threads !" << endl);
305  return 1;
306  }
307  // TODO: authorize -fov and -vox without -dim, checking that the number of voxels obtained is an integer
308  // Number of voxels
309  if (m_nbVoxX<=0 || m_nbVoxY<=0 || m_nbVoxZ<=0)
310  {
311  Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Should provide strictly positive number of voxels !" << endl);
312  return 1;
313  }
314  // When FOV and voxel sizes are both not provided
315  if ((m_voxSizeX<=0. || m_voxSizeY<=0. || m_voxSizeZ<=0.) && (m_fovSizeX<=0. || m_fovSizeY<=0. || m_fovSizeZ<=0.))
316  {
317  Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Should provide strictly positive voxel or FOV dimensions !" << endl);
318  return 1;
319  }
320  // When both FOV and voxel sizes are provided
321  if (m_voxSizeX>0. && m_voxSizeY>0. && m_voxSizeZ>0. && m_fovSizeX>0. && m_fovSizeY>0. && m_fovSizeZ>0.)
322  {
323  Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Both FOV and voxels dimensions provided, should not provide both !" << endl);
324  return 1;
325  }
326  // Check output transaxial FOV which is in percent (0 means we do not apply any FOV masking)
327  if (m_fovOutPercent<0. || m_fovOutPercent>100.)
328  {
329  Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Output transaxial FOV percentage must be in ]0:100] !" << endl);
330  return 1;
331  }
332  // Check that the number of axial slices to be masked is between 0 and dimZ/2
333  if (m_nbSliceOutMask<0 || m_nbSliceOutMask>m_nbVoxZ/2)
334  {
335  Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Number of output axial slices to be masked is incorrectly set !" << endl);
336  return 1;
337  }
338  // Check whether the DynamicDataManager object has been instanciated or not
340  {
341  Cerr("***** oImageDimensionsAndQuantification::CheckParameters() -> Error : DynamicDataManager object not initialized !" << endl);
342  return 1;
343  }
344  // Unlock
345  m_checked = true;
346  // End
347  return 0;
348 }
349 
350 // =====================================================================
351 // ---------------------------------------------------------------------
352 // ---------------------------------------------------------------------
353 // =====================================================================
354 
356 {
357  // Mandatory check
358  if (!m_checked)
359  {
360  Cerr("***** oImageDimensionsAndQuantification::Initialize() -> Cannot initialize before a call to CheckParameters() !" << endl);
361  return 1;
362  }
363 
364  // Verbose
365  if (m_verbose>=2) Cout("oImageDimensionsAndQuantification::Initialize() -> Initialize image dimensions, basis functions and quantification" << endl);
366 
367  // Precompute number of voxels in a slice (XY) and total number of voxels (XYZ)
370 
371  // If FOV dimensions are provided, then compute the voxel dimensions
372  if (m_fovSizeX>0. && m_fovSizeY>0. && m_fovSizeZ>0.)
373  {
375  m_voxSizeY = m_fovSizeY / ((FLTNB)m_nbVoxY);
376  m_voxSizeZ = m_fovSizeZ / ((FLTNB)m_nbVoxZ);
377  }
378  // If voxel dimensions are provided, then compute the voxel dimensions
379  else if (m_voxSizeX>0. && m_voxSizeY>0. && m_voxSizeZ>0.)
380  {
382  m_fovSizeY = m_voxSizeY * ((FLTNB)m_nbVoxY);
383  m_fovSizeZ = m_voxSizeZ * ((FLTNB)m_nbVoxZ);
384  }
385 
386  // Deal with frame list and quantification factors
388  {
389  Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing framing and quantification tabs !" << endl);
390  return 1;
391  }
392 
393  // Deal with time basis functions
395  {
396  Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing time basis functions !" << endl);
397  return 1;
398  }
399 
400  // Deal with respiratory basis functions
402  {
403  Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing respiratory basis functions !" << endl);
404  return 1;
405  }
406 
407  // Deal with cardiac basis functions
409  {
410  Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing cardiac basis functions !" << endl);
411  return 1;
412  }
413 
414  // Deal with ignored corrections
416  {
417  Cerr("***** oImageDimensionsAndQuantification::Initialize() -> A problem occured while initializing ignored corrections !" << endl);
418  return 1;
419  }
420 
421  // Verbose
422  if (m_verbose>=2)
423  {
424  Cout(" --> Image dimensions: [" << m_nbVoxX << ";" << m_nbVoxY << ";" << m_nbVoxZ << "] voxels of [" << m_voxSizeX << ";" << m_voxSizeY << ";" << m_voxSizeZ << "] mm3" << endl);
425  Cout(" --> FOV size: [" << m_fovSizeX << ";" << m_fovSizeY << ";" << m_fovSizeZ << "] mm3" << endl);
426  if (m_nbTimeFrames>1)
427  {
428  if (m_timeStaticFlag) Cout(" --> Time frames: " << m_nbTimeFrames << endl);
429  else Cout(" --> Time frames: " << m_nbTimeFrames << " | Time basis functions: " << m_nbTimeBasisFunctions << endl);
430  }
431  if (m_respBasisFunctionsFile!="")
432  {
433  if (m_respStaticFlag) Cout(" --> Respiratory gates: " << m_nbRespGates << endl);
434  else Cout(" --> Respiratory gates: " << m_nbRespGates << " | Respiratory basis functions: " << m_nbRespBasisFunctions << endl);
435  }
436  if (m_cardBasisFunctionsFile!="")
437  {
438  if (m_cardStaticFlag) Cout(" --> Cardiac gates: " << m_nbCardGates << endl);
439  else Cout(" --> Cardiac gates: " << m_nbCardGates << " | Cardiac basis functions: " << m_nbCardBasisFunctions << endl);
440  }
442  {
443  if (m_nbThreadsForImageComputation==m_nbThreadsForProjection) Cout(" --> Number of parallel threads: " << m_nbThreadsForProjection << endl);
444  else Cout(" --> Number of parallel threads for projection / image computation: [" << m_nbThreadsForProjection << "/" << m_nbThreadsForImageComputation << "]" << endl);
445  }
446  }
447 
448  // Initialized
449  m_initialized = true;
450 
451  // End
452  return 0;
453 }
454 
455 // =====================================================================
456 // ---------------------------------------------------------------------
457 // ---------------------------------------------------------------------
458 // =====================================================================
459 
461 {
462  // Loop on all datafiles
463  for (int bed=0; bed<m_nbBeds; bed++)
464  {
465  // If the number of events in this datafile is below the number of projection threads, then reduce the number of projection threads
466  // to this number
467  if (a2p_DataFile[bed]->GetSize()<m_nbThreadsForProjection)
468  {
469  m_nbThreadsForProjection = a2p_DataFile[bed]->GetSize();
470  Cerr("!!!!! oImageDimensionsAndQuantification::CheckNumberOfProjectionThreadsConsistencyWithDataFileSize() !!!!!" << endl);
471  Cerr(" --> The number of projection threads was reduced to the provided datafile's number of events: " << m_nbThreadsForProjection << endl);
472  }
473  }
474 }
475 
476 // =====================================================================
477 // ---------------------------------------------------------------------
478 // ---------------------------------------------------------------------
479 // =====================================================================
480 
482 {
483  // Note that the consistency between bed positions should have already been tested when this function is called.
484 
485  // First particular case for a single bed
486  if (m_nbBeds==1)
487  {
488  // Here we just copy the relative bed position if provided by the datafile header. It will be ignored in the projection when
489  // there is only one bed position, so no worries.
490  mp_bedPositions[0] = a2p_DataFile[0]->GetRelativeBedPosition();
491  // Flag to say that the bed positions were provided from the datafile
492  m_providedBedPosition = a2p_DataFile[0]->GetBedPositionFlag();
493  // We directly exit here, because below in this function, if the relative bed position has not been provided in the datafile header,
494  // then the scanner default bed displacement will be used, however, this value is not mandatory as some scanners may not be designed
495  // for wholebody.
496  return 0;
497  }
498 
499  // Case where the bed relative positions were provided for all the datafiles
500  if (a2p_DataFile[0]->GetBedPositionFlag())
501  {
502  // In this case we have to compute the center of mass of the provided bed positions, and move it to the 0-center of the Z axis
503  FLTNB center = 0.;
504  // Compute the center of mass
505  for (int bed=0; bed<m_nbBeds; bed++) center += a2p_DataFile[bed]->GetRelativeBedPosition();
506  center /= ((FLTNB)m_nbBeds);
507  // Compute the new bed positions from their relative values to recenter the mass at 0
508  for (int bed=0; bed<m_nbBeds; bed++) mp_bedPositions[bed] = a2p_DataFile[bed]->GetRelativeBedPosition() - center;
509  // Flag to say that the bed positions were provided from the datafile
510  m_providedBedPosition = true;
511  }
512 
513  // Case where the scanner default bed displacement is used
514  else
515  {
516  // Get the scanner object in use
518  // Check that the default bed displacement is more than 0. in the scanner
519  if (ap_Scanner->GetDefaultBedDisplacementInMm()<=0.)
520  {
521  Cerr("***** oImageDimensionsAndQuantification::DealWithBedPositions() -> Bed displacement between two successive bed positions must be strictly positive !" << endl);
522  return 1;
523  }
524  // Loop on the bed positions
525  for (int bed=0; bed<m_nbBeds; bed++)
526  {
527  // Compute bed offset (remember here that the 0. on the Z axis is at the center of the whole image)
528  FLTNB bed_offset = 0.;
529  // For odd numbers of bed positions
530  if (m_nbBeds%2==1) bed_offset = ((FLTNB)( bed-m_nbBeds/2 )) * ap_Scanner->GetDefaultBedDisplacementInMm();
531  // For even numbers of bed positions
532  else bed_offset = (((FLTNB)( bed-m_nbBeds/2 )) + 0.5) * ap_Scanner->GetDefaultBedDisplacementInMm();
533  // Record the value
534  mp_bedPositions[bed] = bed_offset;
535  }
536  // Flag to say that the bed positions were not provided from the datafile
537  m_providedBedPosition = false;
538  }
539 
540  // Verbose
542  {
543  Cout("oImageDimensionsAndQuantification::DealWithBedPositions() -> Use following relative bed positions:" << endl);
544  for (int bed=0; bed<m_nbBeds; bed++) Cout(" --> Bed " << bed << " | Relative axial position: " << mp_bedPositions[bed] << " mm" << endl);
545  }
546 
547  // End
548  return 0;
549 }
550 
551 // =====================================================================
552 // ---------------------------------------------------------------------
553 // ---------------------------------------------------------------------
554 // =====================================================================
555 
557 {
558 // TODO: for whole body 4D dynamic acquisitions, will need to read a set of frames for each bed where the number of frames and frame durations
559 // should be the same for all beds, except which frame is dead and time start and stop also... Have to think about and manage this.
560  // This function requires that the number of beds be set
561  if (m_nbBeds<=0)
562  {
563  Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Number of beds must be set before setting the framing of the acquisition !" << endl);
564  return 1;
565  }
566 
567  // --------------------------------------------------
568  // Allocate time start, stop and duration, to 1 frame
569  m2p_frameDurationsInMs = (uint32_t**)malloc(m_nbBeds*sizeof(uint32_t*));
570  m2p_frameTimeStartInMs = (uint32_t**)malloc(m_nbBeds*sizeof(uint32_t*));
571  m2p_frameTimeStopInMs = (uint32_t**)malloc(m_nbBeds*sizeof(uint32_t*));
572  m3p_quantificationFactors = (FLTNB***)malloc(m_nbBeds*sizeof(FLTNB**));
573  for (int bed=0; bed<m_nbBeds; bed++)
574  {
575  m2p_frameDurationsInMs[bed] = (uint32_t*)malloc(1*sizeof(uint32_t));
576  m2p_frameTimeStartInMs[bed] = (uint32_t*)malloc(1*sizeof(uint32_t));
577  m2p_frameTimeStopInMs[bed] = (uint32_t*)malloc(1*sizeof(uint32_t));
578  }
579 
580  // ---------------------------------------
581  // Particular case where the list is empty
582  if (m_frameList=="")
583  {
584  // Set number of frames to 1
585  m_nbTimeFrames = 1;
586  // Init them to 0 (the vDataFile will set them to the appropriate value later through the oImageDimensionsAndQuantification::SetAcquisitionTime() function)
587  for (int bed=0; bed<m_nbBeds; bed++)
588  {
589  m2p_frameDurationsInMs[bed][0] = 0;
590  m2p_frameTimeStartInMs[bed][0] = 0;
591  m2p_frameTimeStopInMs[bed][0] = 0;
592  // Allocate quantification factors and set them to 1.
593  m3p_quantificationFactors[bed] = (FLTNB**)malloc(1*sizeof(FLTNB*));
595  for (int g=0; g<m_nbRespGates*m_nbCardGates; g++) m3p_quantificationFactors[bed][0][g] = 1.;
596  }
597  // Exit the function
598  return 0;
599  }
600 
601  // --------------------------------------------------------------------
602  // Otherwise, get the parameter as a non-constant string and process it
603  string frame_list = m_frameList;
604 
605  // Current time start
606  uint32_t time_start = 0;
607  // Init number of frames
608  m_nbTimeFrames = 0;
609  // Loop on all commas found in the list
610  size_t comma_pos = 0;
611  while ((comma_pos=frame_list.find_first_of(","))!=string::npos)
612  {
613  // Extract the parameter before the comma
614  string param = frame_list.substr(0,comma_pos);
615  // Search for a dead-frame sign "df"
616  size_t semi_column_pos = param.find("df");
617  // If does not find then this is an actual frame
618  if (semi_column_pos==string::npos)
619  {
620  // Increment number of frames
621  m_nbTimeFrames++;
622  // Realloc the time start, stop and duration
623  for (int bed=0; bed<m_nbBeds; bed++)
624  {
625  m2p_frameDurationsInMs[bed] = (uint32_t*)realloc( m2p_frameDurationsInMs[bed], m_nbTimeFrames*sizeof(uint32_t) );
626  m2p_frameTimeStartInMs[bed] = (uint32_t*)realloc( m2p_frameTimeStartInMs[bed], m_nbTimeFrames*sizeof(uint32_t) );
627  m2p_frameTimeStopInMs[bed] = (uint32_t*)realloc( m2p_frameTimeStopInMs[bed] , m_nbTimeFrames*sizeof(uint32_t) );
628  }
629  // Extract frame duration
630  FLTNB duration = atof(param.c_str());
631 
632  if (duration<=0.)
633  {
634  Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Found a duration which is negative or null !" << endl);
635  return 1;
636  }
637  // Affect time start, stop and duration
638  for (int bed=0; bed<m_nbBeds; bed++)
639  {
640  m2p_frameDurationsInMs[bed][m_nbTimeFrames-1] = ((uint32_t)(duration*1000.));
641  m2p_frameTimeStartInMs[bed][m_nbTimeFrames-1] = time_start;
642  }
643  time_start += m2p_frameDurationsInMs[0][m_nbTimeFrames-1];
644  for (int bed=0; bed<m_nbBeds; bed++)
645  m2p_frameTimeStopInMs[bed][m_nbTimeFrames-1] = time_start;
646  }
647  // Otherwise, this is a dead frame
648  else
649  {
650  // Check that the keyword "df" if at the end of the parameter
651  if (semi_column_pos!=param.size()-2)
652  {
653  Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Invalid dead frame format '" << param << "' !" << endl);
654  return 1;
655  }
656  // Extract the dead frame duration
657  FLTNB duration = atof( (param.substr(0,semi_column_pos)).c_str() );
658 
659  if (duration<=0.)
660  {
661  Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Found a duration which is negative or null !" << endl);
662  return 1;
663  }
664  // Increment time start
665  time_start += ((uint32_t)(duration*1000.));
666  }
667  // Remove first parameter
668  frame_list = frame_list.substr(comma_pos+1);
669  }
670  // Last parameter to extract; search for a dead-frame sign "df"
671  size_t semi_column_pos = frame_list.find("df");
672  // If does not find then this is an actual frame
673  if (semi_column_pos==string::npos)
674  {
675  // Increment number of frames
676  m_nbTimeFrames++;
677  // Realloc the time start, stop and duration
678  for (int bed=0; bed<m_nbBeds; bed++)
679  {
680  m2p_frameDurationsInMs[bed] = (uint32_t*)realloc(m2p_frameDurationsInMs[bed],m_nbTimeFrames*sizeof(uint32_t));
681  m2p_frameTimeStartInMs[bed] = (uint32_t*)realloc(m2p_frameTimeStartInMs[bed],m_nbTimeFrames*sizeof(uint32_t));
682  m2p_frameTimeStopInMs[bed] = (uint32_t*)realloc(m2p_frameTimeStopInMs[bed] ,m_nbTimeFrames*sizeof(uint32_t));
683  }
684  // Extract frame duration
685  FLTNB duration = atof(frame_list.c_str());
686  if (duration<=0.)
687  {
688  Cerr("***** oImageDimensionsAndQuantification::InitializeFramingAndQuantification() -> Found a duration which is negative or null !" << endl);
689  return 1;
690  }
691  // Affect time start, stop and duration
692  for (int bed=0; bed<m_nbBeds; bed++)
693  {
694  m2p_frameDurationsInMs[bed][m_nbTimeFrames-1] = ((uint32_t)(duration*1000.));
695  m2p_frameTimeStartInMs[bed][m_nbTimeFrames-1] = time_start;
696  }
697  time_start += m2p_frameDurationsInMs[0][m_nbTimeFrames-1];
698  for (int bed=0; bed<m_nbBeds; bed++)
699  m2p_frameTimeStopInMs[bed][m_nbTimeFrames-1] = time_start;
700  }
701  // Otherwise, this is a dead frame, so we just ignore it as it is in last position
702  else
703  {
704  ;
705  }
706 
707  // ----------------------------------------------
708  // Allocate and affect the quantification factors
709  for (int bed=0; bed<m_nbBeds; bed++)
710  {
711  m3p_quantificationFactors[bed] = (FLTNB**)malloc(m_nbTimeFrames*sizeof(FLTNB*));
712  for (int fr=0; fr<m_nbTimeFrames; fr++)
713  {
714  m3p_quantificationFactors[bed][fr] = (FLTNB*)malloc(m_nbRespGates*m_nbCardGates*sizeof(FLTNB));
715  for (int g=0; g<m_nbRespGates*m_nbCardGates; g++) m3p_quantificationFactors[bed][fr][g] = 1.;
716  }
717  }
718 
719  // ---
720  // End
721  return 0;
722 }
723 
724 // =====================================================================
725 // ---------------------------------------------------------------------
726 // ---------------------------------------------------------------------
727 // =====================================================================
728 
730 {
731  // Case 1: no time basis functions file provided
732  if (m_timeBasisFunctionsFile=="")
733  {
734  // Set the number of time basis functions to the number of frames
736  // Allocate the conversion matrix from basis functions to frames
738  for (int tbf=0; tbf<m_nbTimeBasisFunctions; tbf++)
739  {
740  m2p_timeBasisFunctions[tbf] = (FLTNB*)malloc(m_nbTimeFrames*sizeof(FLTNB));
741  for (int fr=0; fr<m_nbTimeFrames; fr++)
742  {
743  // Set diagonal to 1, the rest to 0
744  if (tbf==fr) m2p_timeBasisFunctions[tbf][fr] = 1.;
745  else m2p_timeBasisFunctions[tbf][fr] = 0.;
746  }
747  }
748  // Set the static flag to true
749  m_timeStaticFlag = true;
750  }
751  // Case 2: there is a file
752  else
753  {
754  // Allocate the conversion matrix from basis functions to frames
756  for (int tbf=0; tbf<m_nbTimeBasisFunctions; tbf++) m2p_timeBasisFunctions[tbf] = (FLTNB*)malloc(m_nbTimeFrames*sizeof(FLTNB));
757  // Open file
758  ifstream fin(m_timeBasisFunctionsFile.c_str());
759  if (!fin)
760  {
761  Cerr("***** oImageDimensionsAndQuantification::InitializeTimeBasisFunctions() -> Input time basis functions file '" << m_timeBasisFunctionsFile << "' is missing or corrupted !" << endl);
762  return 1;
763  }
764  // Read values
765  for (int tbf=0; tbf<m_nbTimeBasisFunctions; tbf++) for (int fr=0; fr<m_nbTimeFrames; fr++) fin >> m2p_timeBasisFunctions[tbf][fr];
766  // Close file
767  fin.close();
768  }
769  // End
770  return 0;
771 }
772 
773 // =====================================================================
774 // ---------------------------------------------------------------------
775 // ---------------------------------------------------------------------
776 // =====================================================================
777 
779 {
780  // Case 1: no respiratory basis functions file provided
781  if (m_respBasisFunctionsFile=="")
782  {
783  // Set the number of respiratory basis functions to the number of respiratory gates
785 
786  // Allocate the conversion matrix from basis functions to respiratory gates
788  for (int rbf=0; rbf<m_nbRespBasisFunctions; rbf++)
789  {
790  m2p_respBasisFunctions[rbf] = (FLTNB*)malloc(m_nbRespGates*sizeof(FLTNB));
791  for (int rg=0; rg<m_nbRespGates; rg++)
792  {
793  // Set diagonal to 1, the rest to 0
794  if (rbf==rg) m2p_respBasisFunctions[rbf][rg] = 1.;
795  else m2p_respBasisFunctions[rbf][rg] = 0.;
796  }
797  }
798  // Set the static flag to true
799  m_respStaticFlag = true;
800  }
801  // Case 2: there is a file
802  else
803  {
804  // Allocate the conversion matrix from basis functions to respiratory gates
806  for (int rbf=0; rbf<m_nbRespBasisFunctions; rbf++) m2p_respBasisFunctions[rbf] = (FLTNB*)malloc(m_nbRespGates*sizeof(FLTNB));
807  // Open file
808  ifstream fin(m_respBasisFunctionsFile.c_str());
809  if (!fin)
810  {
811  Cerr("***** oImageDimensionsAndQuantification::InitializeRespBasisFunctions() -> Input respiratory basis functions file '" << m_respBasisFunctionsFile << "' is missing or corrupted !" << endl);
812  return 1;
813  }
814  // Read values
815  for (int rbf=0; rbf<m_nbRespBasisFunctions; rbf++) for (int rg=0; rg<m_nbRespGates; rg++) fin >> m2p_respBasisFunctions[rbf][rg];
816  // Close file
817  fin.close();
818  }
819  // End
820  return 0;
821 }
822 
823 // =====================================================================
824 // ---------------------------------------------------------------------
825 // ---------------------------------------------------------------------
826 // =====================================================================
827 
829 {
830  // Case 1: no cardiac basis functions file provided
831  if (m_cardBasisFunctionsFile=="")
832  {
833  // Set the number of cardiac basis functions to the number of cardiac gates
835  // Allocate the conversion matrix from cardiac basis functions to cardiac gates
837  for (int cbf=0; cbf<m_nbCardBasisFunctions; cbf++)
838  {
839  m2p_cardBasisFunctions[cbf] = (FLTNB*)malloc(m_nbCardGates*sizeof(FLTNB));
840  for (int cg=0; cg<m_nbCardGates; cg++)
841  {
842  // Set diagonal to 1, the rest to 0
843  if (cbf==cg) m2p_cardBasisFunctions[cbf][cg] = 1.;
844  else m2p_cardBasisFunctions[cbf][cg] = 0.;
845  }
846  }
847  // Set the static flag to true
848  m_cardStaticFlag = true;
849  }
850  // Case 2: there is a file
851  else
852  {
853  // Allocate the conversion matrix from basis functions to cardiac gates
855  for (int cbf=0; cbf<m_nbCardBasisFunctions; cbf++) m2p_cardBasisFunctions[cbf] = (FLTNB*)malloc(m_nbCardGates*sizeof(FLTNB));
856  // Open file
857  ifstream fin(m_cardBasisFunctionsFile.c_str());
858  if (!fin)
859  {
860  Cerr("***** oImageDimensionsAndQuantification::InitializeCardBasisFunctions() -> Input cardiac basis functions file '" << m_cardBasisFunctionsFile << "' is missing or corrupted !" << endl);
861  return 1;
862  }
863  // Read values
864  for (int cbf=0; cbf<m_nbCardBasisFunctions; cbf++) for (int cg=0; cg<m_nbCardGates; cg++) fin >> m2p_cardBasisFunctions[cbf][cg];
865  // Close file
866  fin.close();
867  }
868  // End
869  return 0;
870 }
871 
872 // =====================================================================
873 // ---------------------------------------------------------------------
874 // ---------------------------------------------------------------------
875 // =====================================================================
876 
878 {
879  // If the m_ignoredCorrectionsList is empty, force all booleans to false and leave this function
880  if (m_ignoredCorrectionsList=="")
881  {
890  return 0;
891  }
892 
893  // Then, count the number of commas in the m_ignoredCorrectionsList to know the number of keywords
894  size_t nb_keywords = count(m_ignoredCorrectionsList.begin(), m_ignoredCorrectionsList.end(), ',') + 1;
895 
896  // Read the keywords
897  string *p_keywords = new string[nb_keywords];
898 
899  if (ReadStringOption(m_ignoredCorrectionsList, p_keywords, nb_keywords, ",", "Ignored corrections list"))
900  {
901  Cerr("***** oImageDimensionsAndQuantification::InitializeIgnoredCorrections() -> An error occured while reading the list of ignored corrections !" << endl);
902  return 1;
903  }
904 
905  // Process them
906  for (size_t k=0; k<nb_keywords; k++)
907  {
908  // Test the different keywords
909  if (p_keywords[k]=="attn") m_ignoreAttnCorrectionFlag = true;
910  else if (p_keywords[k]=="norm") m_ignoreNormCorrectionFlag = true;
911  else if (p_keywords[k]=="scat") m_ignoreScatCorrectionFlag = true;
912  else if (p_keywords[k]=="rand") m_ignoreRandCorrectionFlag = true;
913  else if (p_keywords[k]=="deca") m_ignoreDecaCorrectionFlag = true;
914  else if (p_keywords[k]=="brat") m_ignoreBratCorrectionFlag = true;
915  else if (p_keywords[k]=="fdur") m_ignoreFdurCorrectionFlag = true;
916  else if (p_keywords[k]=="cali") m_ignoreCaliCorrectionFlag = true;
917  else
918  {
919  Cerr("***** oImageDimensionsAndQuantification::InitializeIgnoredCorrections() -> Unknown keyword '" << p_keywords[k] << "' in the provided ignored corrections list !" << endl);
920  return 1;
921  }
922  }
923 
924  delete[] p_keywords;
925  // Normal end
926  return 0;
927 }
928 
929 // =====================================================================
930 // ---------------------------------------------------------------------
931 // ---------------------------------------------------------------------
932 // =====================================================================
933 
934 int oImageDimensionsAndQuantification::SetAcquisitionTime(int a_bed, FLTNB a_timeStartInSec, FLTNB a_durationInSec)
935 {
936  // Check if initialized
937  if (!m_initialized)
938  {
939  Cerr("***** oImageDimensionsAndQuantification::SetAcquisitionTime() -> Object not initialized !" << endl);
940  return 1;
941  }
942  // If the number of frames is more than one, then it was already specified, so we just update the quantification factors and exit
943  if (m_nbTimeFrames>1)
944  {
946  {
947  for (int fr=0; fr<m_nbTimeFrames; fr++) for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++)
948  m3p_quantificationFactors[a_bed][fr][g] /= (((FLTNB)(m2p_frameDurationsInMs[a_bed][fr]))/1000.);
949  }
950  }
951  // If timeStart, timeStop and duration are 0, then it means that it was not specified yet
952  else if (m2p_frameDurationsInMs[a_bed][0]==0 && m2p_frameTimeStartInMs[a_bed][0]==0 && m2p_frameTimeStopInMs[a_bed][0]==0)
953  {
954  m2p_frameDurationsInMs[a_bed][0] = ((uint32_t)(a_durationInSec*1000.));
955  m2p_frameTimeStartInMs[a_bed][0] = ((uint32_t)(a_timeStartInSec*1000.));
956  m2p_frameTimeStopInMs[a_bed][0] = m2p_frameTimeStartInMs[a_bed][0] + m2p_frameDurationsInMs[a_bed][0];
957  // Apply frame's duration onto quantification factors
959  {
960  for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++)
961  m3p_quantificationFactors[a_bed][0][g] /= a_durationInSec;
962  }
963  }
964  // Verbose
965  if (m_verbose>=2 && (a_bed==m_nbBeds-1))
966  {
967  // Case 1: a static single bed
968  if (m_nbTimeFrames==1 && m_nbBeds==1)
969  {
970  Cout("oImageDimensionsAndQuantification::SetAcquisitionTime() -> Static single bed acquisition with duration [ " << GetFrameTimeStartInSec(0,0) << " : "
971  << GetFrameTimeStopInSec(0,0) << " ] seconds" << endl);
972  }
973  // Case 2: a static multi beds
974  else if (m_nbTimeFrames==1 && m_nbBeds>1)
975  {
976  Cout("oImageDimensionsAndQuantification::SetAcquisitionTime() -> Static " << m_nbBeds << " beds acquisition with following bed durations:" << endl);
977  for (int bed=0; bed<m_nbBeds; bed++)
978  Cout(" --> Bed " << bed+1 << " with duration [ " << GetFrameTimeStartInSec(bed,0) << " : " << GetFrameTimeStopInSec(bed,0) << " ] seconds" << endl);
979  }
980  // Case 3: a dynamic single bed
981  else if (m_nbTimeFrames>1 && m_nbBeds==1)
982  {
983  Cout("oImageDimensionsAndQuantification::SetAcquisitionTime() -> Dynamic single bed acquisition with following " << m_nbTimeFrames << " frame durations:" << endl);
984  for (int fr=0; fr<m_nbTimeFrames; fr++)
985  Cout(" --> Frame " << fr+1 << " with duration [ " << GetFrameTimeStartInSec(0,fr) << " : " << GetFrameTimeStopInSec(0,fr) << " ] seconds" << endl);
986  }
987  // Case 4: a dynamic multi beds
988  else
989  {
990  Cout("oImageDimensionsAndQuantification::SetAcquisitionTime() -> Dynamic " << m_nbBeds << " beds acquistion with following " << m_nbTimeFrames << " frame durations:" << endl);
991  for (int bed=0; bed<m_nbBeds; bed++)
992  {
993  Cout(" --> Bed " << bed+1 << " as following framing:" << endl);
994  for (int fr=0; fr<m_nbTimeFrames; fr++)
995  Cout(" | Frame " << fr+1 << " with duration [ " << GetFrameTimeStartInSec(bed,fr) << " : " << GetFrameTimeStopInSec(bed,fr) << " ] seconds" << endl);
996  }
997  }
998  // Correct for frame duration or not
999  if (m_ignoreFdurCorrectionFlag) Cout(" --> Ignore frame duration correction" << endl);
1000  else Cout(" --> Correct for frame duration" << endl);
1001  }
1002  // End
1003  return 0;
1004 }
1005 
1006 // =====================================================================
1007 // ---------------------------------------------------------------------
1008 // ---------------------------------------------------------------------
1009 // =====================================================================
1010 
1012 {
1013  // Check if initialized
1014  if (!m_initialized)
1015  {
1016  Cerr("***** oImageDimensionsAndQuantification::SetCalibrationFactor() -> Object not initialized !" << endl);
1017  return 1;
1018  }
1019  // Ignore calibration factor if asked for
1021  {
1022  // Verbose
1023  if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout("oImageDimensionsAndQuantification::SetCalibrationFactor() -> Ignore calibration factor correction" << endl);
1024  // Exit the function
1025  return 0;
1026  }
1027  // Check if calibration factor is strictly positive
1028  if (a_calibrationFactor<=0.)
1029  {
1030  Cerr("***** oImageDimensionsAndQuantification::SetCalibrationFactor() -> Provided calibration factor (" << a_calibrationFactor << ") is negative or null !" << endl);
1031  return 1;
1032  }
1033  // Affect quantification factor for all frames for this bed (even though the calibration factor should be the same
1034  // for all beds, we do not check it, it should be self consistent in the input files)
1036  {
1037  for (int fr=0; fr<m_nbTimeFrames; fr++) for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++)
1038  m3p_quantificationFactors[a_bed][fr][g] *= a_calibrationFactor;
1039  }
1040  // Verbose
1041  if (m_verbose>=2 && a_bed==m_nbBeds-1)
1042  Cout("oImageDimensionsAndQuantification::SetCalibrationFactor() -> Correct for following calibration factor: " << a_calibrationFactor << endl);
1043  // End
1044  return 0;
1045 }
1046 
1047 // =====================================================================
1048 // ---------------------------------------------------------------------
1049 // ---------------------------------------------------------------------
1050 // =====================================================================
1051 
1053 {
1054  // Check if initialized
1055  if (!m_initialized)
1056  {
1057  Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> Object not initialized !" << endl);
1058  return 1;
1059  }
1060 
1061  // Provide the DynamicDataManager with a 2D array containing quantitative factors to update during initialization
1063 
1064  // Exit the function if no file provided
1065  if (a_quantificationFile=="") return 0;
1066 
1067  // Currently, it is not possible to reconstruction a multi-bed gated acquisition, so we throw an error in this function
1068  if (m_nbBeds>1)
1069  {
1070  Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> Multi-bed gated acquisitions cannot be reconstructed !" << endl);
1071  return 1;
1072  }
1073 
1074  // Verbose
1075  if (m_verbose>=2) Cout("oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors()-> Processing quantification file '" << a_quantificationFile << "'" << endl);
1076 
1077  // GET QUANTIFICATION CORRECTION FACTORS WITH RESPECT TO FRAMES or GATES
1078 
1079  // READ USER-DEFINED QUANTIFICATION FACTORS
1080 
1081  // TODO : For now we assume the configuration file holds :
1082  // TODO : - frame quantification factors : specific to each bed & frame
1083  // TODO : - gate quantification factors : specific to each frame and gates
1084  // TODO : Might require a new ReadDataASCIIFile function to recover this
1085  //
1086  // TODO : Additionally, we should be careful in documentation and explained which gate quantitfication factors are natively taken into account or not
1087 
1088 
1089  // Quantification factors for each bed/frame/gate, intialized to 1.
1090  FLTNB*** pp_dynamic_quantification_factors = new FLTNB**[m_nbBeds];
1091  for (int bed=0 ; bed<m_nbBeds ; bed++)
1092  {
1093  pp_dynamic_quantification_factors[bed] = new FLTNB*[m_nbTimeFrames];
1094  for (int fr=0; fr<m_nbTimeFrames; fr++)
1095  {
1096  pp_dynamic_quantification_factors[bed][fr] = new FLTNB[m_nbRespGates*m_nbCardGates];
1097  for (int g=0; g<m_nbRespGates*m_nbCardGates; g++) pp_dynamic_quantification_factors[bed][fr][g] = 1.;
1098  }
1099  }
1100 
1101  // Build bed-related key words
1102  string *bed_name = new string[m_nbBeds + 1];
1103 
1104  for(int bed=0 ; bed<m_nbBeds ; bed++)
1105  {
1106  ostringstream oss( ostringstream::out );
1107  oss << "bed" << bed+1;
1108  bed_name[bed] = oss.str();
1109  }
1110 
1111  bed_name[m_nbBeds] = "eof";
1112 
1113  for (int bed=0 ; bed<m_nbBeds ; bed++)
1114  {
1115  int return_value = ReadDataASCIIFile(a_quantificationFile, "QUANTIFICATION_FACTORS", pp_dynamic_quantification_factors[bed], m_nbRespGates*m_nbCardGates, m_nbTimeFrames, KEYWORD_MANDATORY, bed_name[bed], bed_name[bed+1]);
1116  if (return_value<0) // string not found error
1117  {
1118  Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> Didn't found quantitative factors in file " << a_quantificationFile << " !" << endl);
1119  return 1;
1120  }
1121  else if(return_value == 1) // reading error
1122  {
1123  Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> An error occured while trying to recover specific quantitative factors for frame !" << endl);
1124  return 1;
1125  }
1126  else if(return_value == 0) // correct reading
1127  {
1128  for(int bed=0 ; bed<m_nbBeds ; bed++)
1129  for (int fr=0; fr<m_nbTimeFrames; fr++)
1130  for (int g=0; g<m_nbRespGates*m_nbCardGates; g++)
1131  if (pp_dynamic_quantification_factors[bed][fr][g] <= 0)
1132  {
1133  Cerr("***** oImageDimensionsAndQuantification::SetDynamicSpecificQuantificationFactors() -> Provided quantification factor (" << pp_dynamic_quantification_factors[bed][fr][g] << ") is negative or null !" << endl);
1134  return 1;
1135  }
1136  else
1137  {
1138  m3p_quantificationFactors[bed][fr][g] *= pp_dynamic_quantification_factors[bed][fr][g];
1139  }
1140  }
1141  }
1142 
1143  // Delete temporary tabs
1144  for (int bed=0; bed<m_nbBeds; bed++)
1145  {
1146  for (int fr=0; fr<m_nbTimeFrames; fr++) delete pp_dynamic_quantification_factors[bed][fr];
1147  delete[] pp_dynamic_quantification_factors[bed];
1148  }
1149  delete[] pp_dynamic_quantification_factors;
1150  delete[] bed_name;
1151 
1152  // End
1153  return 0;
1154 }
1155 
1156 // =====================================================================
1157 // ---------------------------------------------------------------------
1158 // ---------------------------------------------------------------------
1159 // =====================================================================
1160 
1161 FLTNB oImageDimensionsAndQuantification::GetQuantificationFactor(int a_bed, int a_frame, int a_respGate, int a_cardGate)
1162 {
1163  // Make sure that if respiratory or cardiac motion is enabled, we will recover the factor corresponding to the first image
1164  // TODO TM : perhaps more logical to initialize the oImageDimensions with the actual number of gates even when motion is enabled
1165  // TODO TM : (i.e : m_nbRespGates or m_nbCardGates will be equal to 1, but the sensitivity image (list-mode) will require quantification factor for each gate
1166  // TODO TM : and add GetNbRespImages() functions anywhere in the code we need to know the number of resp/card images to be reconstructed
1167  // SS: If the motion correction is enabled, then all gates are pulled together in the optimization process, so there is no need for any quantification factor specific to each gate...
1168  if (m_nbRespGates == 1) a_respGate = 0;
1169  if (m_nbCardGates == 1) a_cardGate = 0;
1170  return m3p_quantificationFactors[a_bed][a_frame][a_respGate*m_nbCardGates+a_cardGate];
1171 }
1172 
1173 // =====================================================================
1174 // ---------------------------------------------------------------------
1175 // ---------------------------------------------------------------------
1176 // =====================================================================
1177 
1178 int oImageDimensionsAndQuantification::SetSPECTIsotope(int a_bed, const string& a_isotope)
1179 {
1180  // Check if initialized
1181  if (!m_initialized)
1182  {
1183  Cerr("***** oImageDimensionsAndQuantification::SetSPECTIsotope() -> Object not initialized !" << endl);
1184  return 1;
1185  }
1186 
1187  // Not yet implemented
1188 
1189  // Normal end
1190  return 0;
1191 }
1192 
1193 // =====================================================================
1194 // ---------------------------------------------------------------------
1195 // ---------------------------------------------------------------------
1196 // =====================================================================
1197 
1198 int oImageDimensionsAndQuantification::SetPETIsotope(int a_bed, const string& a_isotope)
1199 {
1200  // Check if initialized
1201  if (!m_initialized)
1202  {
1203  Cerr("***** oImageDimensionsAndQuantification::SetPETIsotope() -> Object not initialized !" << endl);
1204  return 1;
1205  }
1206 
1207  // ------------------------------------------------------------------
1208  // Preliminary step is checking if we ignore isotope corrections
1209  // ------------------------------------------------------------------
1210 
1212  {
1213  // Verbose
1214  if (m_verbose>=2 && a_bed==m_nbBeds-1)
1215  Cout("oImageDimensionsAndQuantification::SetPETIsotope() -> Ignore isotope dependent corrections" << endl);
1216  // Exit the function
1217  return 0;
1218  }
1219 
1220  // ------------------------------------------------------------------
1221  // Preliminary step is checking for unknown keyword
1222  // ------------------------------------------------------------------
1223 
1224  // Check if the isotope is named as unknown
1225  if (a_isotope=="UNKNOWN" || a_isotope=="Unknown" || a_isotope=="unknown")
1226  {
1227  // In this case we simply assume perfect branching ratio and infinite half-life
1228  if (m_verbose>=2 && a_bed==m_nbBeds-1)
1229  Cout("oImageDimensionsAndQuantification::SetPETIsotope() -> Un-specified isotope; no decay nor branching ratio correction" << endl);
1230  // Exit the function
1231  return 0;
1232  }
1233 
1234  // ------------------------------------------------------------------
1235  // First step is open the isotopes file and find the provided isotope
1236  // ------------------------------------------------------------------
1237 
1238  // Get oOutputManager instance and config directory
1239  sOutputManager* p_output_manager = sOutputManager::GetInstance();
1240  string config_dir = p_output_manager->GetPathToConfigDir();
1241 
1242  // Open the isotope file based on config directory
1243  string file_name = config_dir + "/misc/isotopes_pet.txt";
1244  ifstream fin(file_name.c_str());
1245 
1246  if (!fin)
1247  {
1248  Cerr("***** oImageDimensionsAndQuantification::SetPETIsotope() -> Failed to open PET isotopes data file '" << file_name << "' !" << endl);
1249  return 1;
1250  }
1251 
1252  // Loop on lines to find the isotope
1253  int line_max_size = 10240;
1254  char *line = new char[line_max_size];
1255  FLTNB half_life = -1.;
1256  FLTNB branching_ratio = -1.;
1257  bool found_it = false;
1258  fin.getline(line,line_max_size);
1259  while (!fin.eof())
1260  {
1261  // For the word position
1262  size_t found_position;
1263  // Cast the line into a string for convenience
1264  string test = (string)line;
1265  // Jump to next line if we find a # character as the first character
1266  if ((found_position=test.find("#"))==0)
1267  {
1268  // Read a new line before continuing
1269  fin.getline(line,line_max_size);
1270  continue;
1271  }
1272  // Check if we see the isotope name
1273  found_position = test.find(a_isotope);
1274  if (found_position!=string::npos)
1275  {
1276  // Each line is organised as follows: (some spaces/tabs) isotope_name (some spaces/tabs) half_life (some spaces/tabs) branching_ratio
1277  // We first remove the isotope name from the line
1278  test = test.substr(found_position+1+a_isotope.length());
1279  // Then we convert it into a string stream to ease the reading of both numbers
1280  istringstream fstr(test);
1281  fstr >> half_life >> branching_ratio;
1282  // We found it
1283  found_it = true;
1284  break;
1285  }
1286  // Read a new line
1287  fin.getline(line,line_max_size);
1288  }
1289 
1290  delete[] line;
1291  // Close file
1292  fin.close();
1293 
1294  // Check if we found it or not
1295  if (found_it)
1296  {
1297  // Check rationality of values
1298  if (branching_ratio<=0. || branching_ratio>1.)
1299  {
1300  Cerr("***** oImageDimensionsAndQuantification::SetPETIsotope() -> Branching ratio (" << branching_ratio << ") is not in the ]0:1] range !" << endl);
1301  return 1;
1302  }
1303  // Verbose
1304  if (m_verbose>=2 && a_bed==m_nbBeds-1)
1305  {
1306  // If negative half-life, then we consider it is infinite (no half-life)
1307  if (half_life<=0.)
1308  Cout("oImageDimensionsAndQuantification::SetPETIsotope() -> Isotope " << a_isotope << " has infinite half life and " << branching_ratio << " branching ratio" << endl);
1309  else
1310  Cout("oImageDimensionsAndQuantification::SetPETIsotope() -> Isotope " << a_isotope << " has " << half_life << " seconds half life and " << branching_ratio << " branching ratio" << endl);
1311  }
1312  }
1313  else
1314  {
1315  // Throw error message
1316  Cerr("***** oImageDimensionsAndQuantification::SetPETIsotope() -> Did not find " << a_isotope << " isotope in the PET isotope data file, please add it !" << endl);
1317  return 1;
1318  }
1319 
1320  // ------------------------------------------------------------------------
1321  // Second step is applying branching ratio and decay wrt to frame durations
1322  // ------------------------------------------------------------------------
1323 
1324  // Branching ratio
1326  {
1327  // Apply correction
1328  for (int fr=0; fr<m_nbTimeFrames; fr++)
1329  {
1330  for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++)
1331  m3p_quantificationFactors[a_bed][fr][g] /= branching_ratio;
1332  }
1333  // Verbose
1334  if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout(" --> Correct for branching ratio" << endl);
1335  }
1336  // Verbose
1337  else if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout(" --> Ignore branching ratio correction" << endl);
1338 
1339  // Apply decay factors if half life is not negative (meaning no half-life)
1340  if (half_life>0.)
1341  {
1342  // We correct for it
1344  {
1345  // Apply correction
1346  for (int fr=0; fr<m_nbTimeFrames; fr++)
1347  {
1348  long double lambda = log(2.0)/half_life;
1349  for(int g=0 ; g<m_nbRespGates*m_nbCardGates ; g++)
1350  {
1351  long double dstart = lambda*GetFrameTimeStartInSec(a_bed,fr);
1352  long double dacq = lambda*GetFrameDurationInSec(a_bed,fr);
1353  // Time start decay correction
1354  m3p_quantificationFactors[a_bed][fr][g] *= exp(dstart);
1355  // Frame duration decay correction
1356  m3p_quantificationFactors[a_bed][fr][g] *= dacq/(1.0-exp(-dacq));
1357 
1358  /* Time start decay correction
1359  m3p_quantificationFactors[a_bed][fr][g] *= exp(lambda*GetFrameTimeStartInSec(a_bed,fr));
1360  // Frame duration decay correction
1361  m3p_quantificationFactors[a_bed][fr][g] *= lambda*GetFrameDurationInSec(a_bed,fr)/(1.0-exp(-lambda*(GetFrameDurationInSec(a_bed,fr))));*/
1362  }
1363  }
1364  // Verbose
1365  if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout(" --> Correct for half-life" << endl);
1366  }
1367  // Verbose
1368  else if (m_verbose>=2 && a_bed==m_nbBeds-1) Cout(" --> Ignore half-life correction" << endl);
1369  }
1370 
1371  // End
1372  return 0;
1373 }
1374 
1375 // =====================================================================
1376 // ---------------------------------------------------------------------
1377 // ---------------------------------------------------------------------
1378 // =====================================================================
1379 
1380 int oImageDimensionsAndQuantification::InitDynamicData( string a_pathTo4DDataSplittingFile,
1381  int a_respMotionCorrectionFlag,
1382  int a_cardMotionCorrectionFlag,
1383  int a_doubleMotionCorrectionFlag,
1384  int a_invMotionCorrectionFlag,
1385  int a_nbRespGates, int a_nbCardGates )
1386 {
1387  if (m_verbose>=5) Cout("oImageDimensionsAndQuantification::InitDynamicData()" << endl);
1390  if (mp_DynamicDataManager->InitDynamicData( a_nbRespGates, a_nbCardGates, a_pathTo4DDataSplittingFile,
1391  a_respMotionCorrectionFlag, a_cardMotionCorrectionFlag, a_doubleMotionCorrectionFlag, a_invMotionCorrectionFlag ))
1392  {
1393  Cerr("***** oImageDimensionsAndQuantification::InitDynamicData() -> A problem occured while initializing the dynamic data from dynamic data manager !" << endl);
1394  return 1;
1395  }
1396  return 0;
1397 }
1398 
1399 // =====================================================================
1400 // ---------------------------------------------------------------------
1401 // ---------------------------------------------------------------------
1402 // =====================================================================
1403 
1405 {
1406  if (m_verbose>=5) Cout("oImageDimensionsAndQuantification::CheckDynamicParameters()" << endl);
1407  return mp_DynamicDataManager->CheckParameters(a_nbEvents);
1408 }
1409 
1410 // =====================================================================
1411 // ---------------------------------------------------------------------
1412 // ---------------------------------------------------------------------
1413 // =====================================================================
int InitializeFramingAndQuantification()
A function used to initialize the framing and quantification tables.
int CheckParameters(int64_t a_nbEvents)
Check all mandatory parameters.
This class is designed to be a mother virtual class for DataFile.
Definition: vDataFile.hh:103
int SetSPECTIsotope(int a_bed, const string &a_isotope)
Set the SPECT isotope for the provided bed.
int SetPETIsotope(int a_bed, const string &a_isotope)
Set the PET isotope for the provided bed.
static sScannerManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
Declaration of class oImageDimensionsAndQuantification.
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...
void SetVerbose(int a_verboseLevel)
set verbosity
oImageDimensionsAndQuantification()
The constructor of oImageDimensionsAndQuantification.
FLTNB GetFrameTimeStartInSec(int a_bed, int a_frame)
Get the frame time start for the given bed, in seconds as a FLTNB.
#define FLTNB
Definition: gVariables.hh:81
This class gathers the information about the dynamic splitting of the data.
int SetCalibrationFactor(int a_bed, FLTNB a_calibrationFactor)
Set the calibration factor for the provided bed.
void CheckNumberOfProjectionThreadsConsistencyWithDataFileSize(vDataFile **a2p_DataFile)
int DealWithBedPositions(vDataFile **a2p_DataFile)
Deal with provided or default bed relative positions.
int64_t GetSize()
Definition: vDataFile.hh:333
static sOutputManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
int SetNbThreads(const string &a_nbThreads)
Set the number of threads.
~oImageDimensionsAndQuantification()
The destructor of oImageDimensionsAndQuantification.
int CheckParameters()
A function used to check the parameters settings.
#define Cerr(MESSAGE)
int Initialize()
A function used to initialize all that is needed.
Singleton class that manages output writing on disk (images, sinograms, etc). It also manages loggi...
const string & GetPathToConfigDir()
Return the path to the CASTOR config directory.
Declaration of class vDataFile.
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 InitializeTimeBasisFunctions()
A function used to initialize the time basis functions.
bool GetBedPositionFlag()
Definition: vDataFile.hh:426
void SetImageDimensionsAndQuantification(oImageDimensionsAndQuantification *ap_ImageDimensionsAndQuantification)
set the pointer to the oImageDimensionsAndQuantification object
FLTNB GetQuantificationFactor(int a_bed, int a_frame, int a_respGate, int a_cardGate)
Get the quantification factor corresponding to the provided bed, frame, respiratory and cardiac gates...
#define VERBOSE_NORMAL
vScanner * GetScannerObject()
int SetDynamicSpecificQuantificationFactors(FLTNB **a2p_quantificationFactors)
Compute gate-specific quantificative factors using the number of events within each gate...
#define KEYWORD_MANDATORY
Definition: gOptions.hh:48
int CheckDynamicParameters(int64_t a_nbEvents)
Call the eponym function from the oDynamicDataManager object in order to check its parameters...
int InitializeRespBasisFunctions()
A function used to initialize the respiratory basis functions.
FLTNB GetDefaultBedDisplacementInMm()
Definition: vScanner.hh:289
FLTNB GetFrameTimeStopInSec(int a_bed, int a_frame)
Get the frame time stop for the given bed, in seconds as a FLTNB.
FLTNB GetRelativeBedPosition()
Definition: vDataFile.hh:432
int SetAcquisitionTime(int a_bed, FLTNB a_timeStartInSec, FLTNB a_durationInSec)
Set the acquisition time if not already set by the SetTimeFrames()
int InitDynamicData(int a_nbRespGates, int a_nbCardGates, const string &a_pathTo4DDataSplittingFile, int a_rmMCorrFlag, int a_cMmCorrFlag, int a_dmCorrFlag, int a_pMotionCorrFlag)
Main function for instanciation and initialization of the member variables and arrays. Call the specific initialization function depending of the type of dataset.
#define Cout(MESSAGE)
int InitializeCardBasisFunctions()
A function used to initialize the cardiac basis functions.
FLTNB GetFrameDurationInSec(int a_bed, int a_frame)
Get the frame duration for the given bed, in seconds as a FLTNB.
void SetDefault()
A function used to set number of threads and MPI instances to 1 and bypass the CheckParameters() func...
int ReadStringOption(const string &a_input, T *ap_return, int a_nbElts, const string &sep, const string &a_option)
Parse the 'a_input' string corresponding to the 'a_option' into 'a_nbElts' elements, using the 'sep' separator. The results are returned in the templated 'ap_return' dynamic templated array. Call "ConvertFromString()" to perform the correct conversion depending on the type of the data to convert.
Definition: gOptions.cc:62
int InitializeIgnoredCorrections()
A function used to initialize the ignored corrections.
Generic class for scanner objects.
Definition: vScanner.hh:62
int SetFlipOut(const string &a_flipOut)
Set the output flip options, the parameter being a string potentially containing the letters x...
int SetDynamicSpecificQuantificationFactors(const string &a_quantificationFile)
Apply specific quantification factors manually provided as an option.