CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
vDataFile.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  mp_ID = NULL;
42  m_verbose = -1;
43 
44  // Variables related to the acquisition
45  m2p_dataFile = NULL;
46  m_headerFileName = "";
47  m_dataFileName = "";
48  m_scannerName = "";
52  m_nbEvents = -1;
53  m_bedIndex = -1;
55  m_bedPositionFlag = false;
56  m_startTimeInSec = 0.;
57  m_durationInSec = -1.;
59  m_sizeEvent = -1;
60 
61  // Default POI (meaning we do not have any)
62  mp_POIResolution[0] = -1.;
63  mp_POIResolution[1] = -1.;
64  mp_POIResolution[2] = -1.;
65  mp_POIDirectionFlag[0] = false;
66  mp_POIDirectionFlag[1] = false;
67  mp_POIDirectionFlag[2] = false;
68  m_POIInfoFlag = false;
69  m_ignorePOIFlag = false;
70 
71  // Variable related to Buffer/Container arrays
72  m2p_BufferEvent = NULL;
73  m_mpi1stEvent = -1;
74  m_mpiLastEvent = -1;
75  m_mpiNbEvents = -1;
76  mp_MappedFile = NULL;
77  mp_mappedMemory = NULL;
78 }
79 
80 // =====================================================================
81 // ---------------------------------------------------------------------
82 // ---------------------------------------------------------------------
83 // =====================================================================
84 
86 {
87  // Free the (vEvent) buffer event
88  if (m2p_BufferEvent)
89  {
90  for(int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
91  if (m2p_BufferEvent[th]) delete m2p_BufferEvent[th];
92  delete[] m2p_BufferEvent;
93  }
94  // Close datafiles and delete them
95  if (m2p_dataFile)
96  {
97  for(int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
98  if (m2p_dataFile[th]) m2p_dataFile[th]->close();
99  delete m2p_dataFile;
100  }
101  // Delete the mapped memory object
102  if (mp_MappedFile!=NULL) delete mp_MappedFile;
103 }
104 
105 // =====================================================================
106 // ---------------------------------------------------------------------
107 // ---------------------------------------------------------------------
108 // =====================================================================
109 
110 int vDataFile::ReadInfoInHeader(bool a_affectQuantificationFlag)
111 {
112  // Verbose
114  if (m_verbose>=VERBOSE_DETAIL) Cout("vDataFile::ReadInfoInHeader() -> Read datafile header from '" << m_headerFileName << " ...'" << endl);
115 
116  // Check that the image dimensions and quantification object has been set if the a_affectQuantificationFlag is true
117  if (a_affectQuantificationFlag && mp_ID==NULL)
118  {
119  Cerr("***** vDataFile::ReadInfoInHeader() -> Image dimensions and quantification object has not been set while it is asked to affect" << endl);
120  Cerr(" its parameters from information read into the header !" << endl);
121  return 1;
122  }
123 
124  // Read mandatory general fields in the header, check if errors (either mandatory tag not found or issue during data reading/conversion (==1) )
125  if (ReadDataASCIIFile(m_headerFileName, "Number of events", &m_nbEvents, 1, KEYWORD_MANDATORY) )
126  {
127  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read number of events in the header data file !" << endl);
128  return 1;
129  }
131  {
132  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read scanner name in the header data file !" << endl);
133  return 1;
134  }
135 
136  // Get data file name
138  {
139  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read data filename in the header data file !" << endl);
140  return 1;
141  }
142  // If it is an absolute path (start with a OS_SEP) then let it as is, otherwise past the path of the header file name (because we suppose the two
143  // files are side by side).
144  // This is currently only unix compatible...
145  if (m_dataFileName.substr(0,1)!=OS_SEP && m_headerFileName.find(OS_SEP)!=string::npos)
146  {
147  // Extract the path from the header file name
148  size_t last_slash = m_headerFileName.find_last_of(OS_SEP);
149  // Paste the path to the data file name
150  m_dataFileName = m_headerFileName.substr(0,last_slash+1) + m_dataFileName;
151  }
152 
153  // For data mode, multiple declarations are allowed, so we get the mode as a string and do some checks
154  string data_mode = "";
155  if (ReadDataASCIIFile(m_headerFileName, "Data mode", &data_mode, 1, KEYWORD_MANDATORY) )
156  {
157  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read data mode in the header data file !" << endl);
158  return 1;
159  }
160  if ( data_mode=="LIST" || data_mode=="LISTMODE" || data_mode=="LIST-MODE" ||
161  data_mode=="list" || data_mode=="listmode" || data_mode=="list-mode" || data_mode=="0" ) m_dataMode = MODE_LIST;
162  else if ( data_mode=="HISTOGRAM" || data_mode=="histogram" || data_mode=="Histogram" ||
163  data_mode=="HISTO" || data_mode=="histo" || data_mode=="Histo" || data_mode=="1" ) m_dataMode = MODE_HISTOGRAM;
164  else if ( data_mode=="NORMALIZATION" || data_mode=="normalization" || data_mode=="Normalization" ||
165  data_mode=="NORM" || data_mode=="norm" || data_mode=="Norm" || data_mode=="2" ) m_dataMode = MODE_NORMALIZATION;
166  else
167  {
168  Cerr("***** vDataFile::ReadInfoInHeader() -> Unknown data mode '" << data_mode << "' found in header file !" << endl);
169  return 1;
170  }
171 
172  // For data type, multiple declarations are allowed, so we get the mode as a string and do some checks
173  string data_type = "";
174  if (ReadDataASCIIFile(m_headerFileName, "Data type", &data_type, 1, KEYWORD_MANDATORY) )
175  {
176  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read data type in the header data file !" << endl);
177  return 1;
178  }
179  if ( data_type=="PET" || data_type=="pet" || data_type=="0" )
180  {
183  }
184  else if ( data_type=="SPECT" || data_type=="spect" || data_type=="1" )
185  {
188  }
189  else if ( data_type=="CT" || data_type=="ct" || data_type=="2" )
190  {
193  }
194  else
195  {
196  Cerr("***** vDataFile::ReadInfoInHeader() -> Unknown data type '" << data_type << "' found in header file !" << endl);
197  return 1;
198  }
199 
200  // Get start time and duration (optional for the normalization mode
202  {
204  {
205  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read acquisition start time in the header data file !" << endl);
206  return 1;
207  }
209  {
210  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read acquisition stop time in the header data file !" << endl);
211  return 1;
212  }
213  }
214  else
215  {
216  if (ReadDataASCIIFile(m_headerFileName, "Start time (s)", &m_startTimeInSec, 1, KEYWORD_OPTIONAL) == 1 )
217  {
218  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while reading acquisition start time in the header data file !" << endl);
219  return 1;
220  }
221  if (ReadDataASCIIFile(m_headerFileName, "Duration (s)", &m_durationInSec, 1, KEYWORD_OPTIONAL) == 1 )
222  {
223  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while reading acquisition stop time in the header data file !" << endl);
224  return 1;
225  }
226  }
227 
228  // Set the acquisition timing to the quantification factors
229  if (a_affectQuantificationFlag && mp_ID->SetAcquisitionTime(m_bedIndex, m_startTimeInSec, m_durationInSec))
230  {
231  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while passing acquisition start and stop time to the ImageDimensionsAndQuantification !" << endl);
232  return 1;
233  }
234 
235  // Read optional bed relative position
236  int read_bed_position_status = ReadDataASCIIFile(m_headerFileName, "Horizontal bed relative position (mm)", &m_relativeBedPosition, 1, KEYWORD_OPTIONAL);
237  if (read_bed_position_status==KEYWORD_OPTIONAL_SUCCESS) m_bedPositionFlag = true;
238  else if (read_bed_position_status==KEYWORD_OPTIONAL_SUCCESS)
239  {
240  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while reading the 'Bed relative position (mm)' field in the data file header !" << endl);
241  return 1;
242  }
243 
244  // Read optional fields in the header, check if errors (issue during data reading/conversion (==1) )
245  if (ReadDataASCIIFile(m_headerFileName, "Calibration factor", &m_calibrationFactor, 1, KEYWORD_OPTIONAL) == 1 ||
247  ReadDataASCIIFile(m_headerFileName, "POI correction flag", &m_POIInfoFlag, 3, KEYWORD_OPTIONAL) == 1 )
248  {
249  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while reading optional field in the header data file !" << endl);
250  return 1;
251  }
252 
253  // Read fields specific to the modality (call to the related function implemented in child classes)
254  if (ReadSpecificInfoInHeader(a_affectQuantificationFlag) )
255  {
256  Cerr("***** vDataFile::ReadInfoInHeader() -> Error while trying to read modality-specific informations from the header data file !" << endl);
257  return 1;
258  }
259 
260  // Give the calibration factor to the oImageDimensionsAndQuantification that manages the quantification factors
261  if (a_affectQuantificationFlag && mp_ID->SetCalibrationFactor(m_bedIndex, m_calibrationFactor))
262  {
263  Cerr("***** vDataFile::ReadSpecificInfoInHeader() -> A problem occured while setting the calibration factor to oImageDimensionsAndQuantification !" << endl);
264  return 1;
265  }
266 
267  return 0;
268 }
269 
270 // =====================================================================
271 // ---------------------------------------------------------------------
272 // ---------------------------------------------------------------------
273 // =====================================================================
274 
276 {
277  // Verbose
278  m_verbose = ap_DataFile->GetVerbose();
279  // Variables related to the acquisition
280  m2p_dataFile = NULL;
281  m_scannerName = ap_DataFile->GetScannerName();
282  m_dataMode = ap_DataFile->GetDataMode();
283  m_dataType = ap_DataFile->GetDataType();
284  m_dataSpec = ap_DataFile->GetDataSpec();
285  m_bedIndex = ap_DataFile->GetBedIndex();
287  m_bedPositionFlag = ap_DataFile->GetBedPositionFlag();
288  m_startTimeInSec = ap_DataFile->GetStartTime();
289  m_durationInSec = ap_DataFile->GetDuration();
290  m_calibrationFactor = ap_DataFile->GetCalibrationFactor();
291  m_sizeEvent = ap_DataFile->GetEventSize();
292  // Default POI (meaning we do not have any)
293  m_POIInfoFlag = ap_DataFile->GetPOIInfoFlag();
294  mp_POIResolution[0] = ap_DataFile->GetPOIResolution()[0];
295  mp_POIResolution[1] = ap_DataFile->GetPOIResolution()[1];
296  mp_POIResolution[2] = ap_DataFile->GetPOIResolution()[2];
297  mp_POIDirectionFlag[0] = ap_DataFile->GetPOIDirectionFlag()[0];
298  mp_POIDirectionFlag[1] = ap_DataFile->GetPOIDirectionFlag()[1];
299  mp_POIDirectionFlag[2] = ap_DataFile->GetPOIDirectionFlag()[2];
300  // Call the specific function
301  if (SetSpecificParametersFrom(ap_DataFile))
302  {
303  Cerr("***** vDataFile::SetParametersFrom() -> An error occured while setting specific parameters fro mthe provided datafile !" << endl);
304  return 1;
305  }
306  // End
307  return 0;
308 }
309 
310 // =====================================================================
311 // ---------------------------------------------------------------------
312 // ---------------------------------------------------------------------
313 // =====================================================================
314 
316 {
317  // Verbose
319  if (m_verbose>=VERBOSE_DETAIL) Cout("vDataFile::CheckParameters() -> Check mandatory parameters" << endl);
320  // Check mandatory parameters
321  if (mp_ID == NULL)
322  {
323  Cerr("***** vDataFile::CheckParameters() -> ImageDimensionsAndQuantification object not initialized !" << endl);
324  return 1;
325  }
326  if (m_headerFileName == "")
327  {
328  Cerr("***** vDataFile::CheckParameters() -> String containing path to header file not initialized !" << endl);
329  return 1;
330  }
331  if (m_dataFileName == "")
332  {
333  Cerr("***** vDataFile::CheckParameters() -> String containing path to raw data file not initialized !" << endl);
334  return 1;
335  }
336  if (m_nbEvents<0)
337  {
338  Cerr("***** vDataFile::CheckParameters() -> Number of events incorrectly initialized !" << endl);
339  return 1;
340  }
342  {
343  Cerr("***** vDataFile::CheckParameters() -> Data mode incorrectly initialized !" << endl);
344  return 1;
345  }
347  {
348  Cerr("***** vDataFile::CheckParameters() -> Acquisition duration (s) incorrectly initialized !" << endl);
349  return 1;
350  }
352  {
353  Cerr("***** vDataFile::CheckParameters() -> Data type incorrectly initialized !" << endl);
354  return 1;
355  }
357  {
358  Cerr("***** vDataFile::CheckParameters() -> Data physical property incorrectly initialized !" << endl);
359  return 1;
360  }
361  if (m_bedIndex<0)
362  {
363  Cerr("***** vDataFile::CheckParameters() -> Bed position index incorrectly initialized !" << endl);
364  return 1;
365  }
366  if (m_verbose<0)
367  {
368  Cerr("***** vDataFile::CheckParameters() -> Verbosity incorrectly initialized !" << endl);
369  return 1;
370  }
371  // Call to the related function implemented in child classes
373  {
374  Cerr("***** vDataFile::CheckParameters() -> Error while checking specific parameters !" << endl);
375  return 1;
376  }
377  // End
378  return 0;
379 }
380 
381 // =====================================================================
382 // ---------------------------------------------------------------------
383 // ---------------------------------------------------------------------
384 // =====================================================================
385 
387 {
388  // Verbose
390  if (m_verbose>=VERBOSE_DETAIL) Cout("vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Check consistency between this and the provided datafiles" << endl);
391  // Check data type
392  if (m_dataType!=ap_DataFile->GetDataType())
393  {
394  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Data types are inconsistent !" << endl);
395  return 1;
396  }
397  // Check data mode
398  if (m_dataMode!=ap_DataFile->GetDataMode())
399  {
400  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Data modes are inconsistent !" << endl);
401  return 1;
402  }
403  // For histogram and normalization modes, check the number of events (i.e. the number of histogram bins)
405  {
406  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Total number of events are inconsistent !" << endl);
407  return 1;
408  }
409  // Check the scanner name
410  if (m_scannerName!=ap_DataFile->GetScannerName())
411  {
412  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Scanner names are inconsistent !" << endl);
413  return 1;
414  }
415  // Check the calibration factor
416  if (m_calibrationFactor!=ap_DataFile->GetCalibrationFactor())
417  {
418  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Calibration factors are inconsistent !" << endl);
419  return 1;
420  }
421  // Check event size
422  if (m_sizeEvent!=ap_DataFile->GetEventSize())
423  {
424  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Events sizes are inconsistent !" << endl);
425  return 1;
426  }
427  // Check POI info flag
428  if (m_POIInfoFlag!=ap_DataFile->GetPOIInfoFlag())
429  {
430  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> POI flags are inconsistent !" << endl);
431  return 1;
432  }
433  // Check POI resolutions
434  if ( (mp_POIResolution[0]!=ap_DataFile->GetPOIResolution()[0]) ||
435  (mp_POIResolution[1]!=ap_DataFile->GetPOIResolution()[1]) ||
436  (mp_POIResolution[2]!=ap_DataFile->GetPOIResolution()[2]) )
437  {
438  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> POI resolutions are inconsistent !" << endl);
439  return 1;
440  }
441  // Check the bed position flag
442  if (m_bedPositionFlag!=ap_DataFile->GetBedPositionFlag())
443  {
444  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Bed relative positions must be provided for all beds or not at all !" << endl);
445  return 1;
446  }
447  // Call specific function
449  {
450  Cerr("***** vDataFile::CheckConsistencyWithAnotherBedDataFile() -> Inconsistency detected for specific characteristics !" << endl);
451  return 1;
452  }
453  // End
454  return 0;
455 }
456 
457 // =====================================================================
458 // ---------------------------------------------------------------------
459 // ---------------------------------------------------------------------
460 // =====================================================================
461 
463 {
464  // Verbose
466  if (m_verbose>=VERBOSE_NORMAL) Cout("vDataFile::InitializeMappedFile() -> Map datafile to memory" << endl);
467 
468  // Check file size consistency here (this is a pure virtual function implemented by children)
470  {
471  Cerr("***** vDataFile::InitializeMappedFile() -> A problem occured while checking file size consistency !" << endl);
472  return 1;
473  }
474 
475  // Create mapped file
477  // Open the file
479  {
480  Cerr("***** vDataFile::InitializeMappedFile() -> Failed to open data file '" << m_dataFileName << "' !" << endl);
481  Cerr(" Provided in data file header '" << m_headerFileName << "'." << endl);
482  return 1;
483  }
484 
485  // Get raw pointer to mapped memory
487 
488  // ------------------ Compute here the piece of file that each MPI instance manages --------------------
489 
490  // 1. Compute the number of events that each instance has to manage
491  int64_t instance_size = m_nbEvents / mp_ID->GetMPISize();
492  // All instances manage an equal part of the file but the last one which also manages the few last events
493  if (mp_ID->GetMPIRank()!=mp_ID->GetMPISize()-1) m_mpiNbEvents = instance_size;
494  else m_mpiNbEvents = instance_size + (m_nbEvents - instance_size*mp_ID->GetMPISize());
495  // 2. Compute the first event managed by the instance
496  m_mpi1stEvent = mp_ID->GetMPIRank() * instance_size;
497  // 3. Compute the last event managed by the instance
498  m_mpiLastEvent = m_mpi1stEvent + m_mpiNbEvents - 1;
499 
500  // End
501  return 0;
502 }
503 
504 
505 // =====================================================================
506 // ---------------------------------------------------------------------
507 // ---------------------------------------------------------------------
508 // =====================================================================
509 
511 {
513  if (m_verbose==0) return;
514  // Describe the datafile
515  Cout("vDataFile::Describe() -> Here is some generic content of the datafile" << endl);
516  Cout(" --> Header file name: " << m_headerFileName << endl);
517  Cout(" --> Data file name: " << m_dataFileName << endl);
518  Cout(" --> Data type: " << GetDataTypeToString() << endl);
519  Cout(" --> Data mode: " << GetDataModeToString() << endl);
520  Cout(" --> Data spec: " << GetDataSpecToString() << endl);
521  Cout(" --> Scanner name: " << m_scannerName << endl);
522  Cout(" --> Number of specific events: " << m_nbEvents << endl);
523  Cout(" --> Size of an event: " << m_sizeEvent << " bytes" << endl);
524  Cout(" --> Time start: " << m_startTimeInSec << " sec" << endl);
525  Cout(" --> Duration: " << m_durationInSec << " sec" << endl);
526  Cout(" --> Calibration factor: " << m_calibrationFactor << endl);
527  if (m_bedPositionFlag) Cout(" --> Relative axial bed position: " << m_relativeBedPosition << " mm" << endl);
528  // Call the specific function
530 }
531 
532 // =====================================================================
533 // ---------------------------------------------------------------------
534 // ---------------------------------------------------------------------
535 // =====================================================================
536 
537 int vDataFile::OpenFileForWriting(string a_suffix)
538 {
539  // Check that file is not already open
540  if (m2p_dataFile && m2p_dataFile[0]->is_open())
541  {
542  Cerr("***** vDataFile::OpenFileForWriting() -> Output data file is already open !" << endl);
543  return 1;
544  }
545 
546  // Allocate the m2p_dataFile for each thread
547  m2p_dataFile = new fstream*[1];
548 
549  // Set the name of output files
551  string path_name = p_OutMgr->GetPathName();
552  string base_name = p_OutMgr->GetBaseName();
553  m_headerFileName = path_name + base_name + a_suffix + ".cdh";
554  m_dataFileName = path_name + base_name + a_suffix + ".cdf";
555 
556  // Verbose
557  if (m_verbose>=2) Cout("vDataFile::OpenFileForWriting() -> Output data file header is '" << m_headerFileName << "'" << endl);
558 
559  // Open file
560  m2p_dataFile[0] = new fstream( m_dataFileName.c_str(), ios::binary | ios::out );
561 
562  // Check
563  if (!m2p_dataFile[0]->is_open())
564  {
565  Cerr("***** vDataFile::OpenFileForWriting() -> Failed to create output file '" << m_dataFileName << "' !" << endl);
566  return 1;
567  }
568 
569  // End
570  return 0;
571 }
572 
573 // =====================================================================
574 // ---------------------------------------------------------------------
575 // ---------------------------------------------------------------------
576 // =====================================================================
577 
579 {
580  // Close binary data file
581  if (m2p_dataFile)
582  {
583  // Verbose
584  if (m_verbose>=2) Cout("vDataFile::CloseFile() -> Closing output binary data file" << endl);
585  // Close file
586  if (m2p_dataFile[0]) m2p_dataFile[0]->close();
587  delete m2p_dataFile;
588  }
589  // End
590  return 0;
591 }
592 
593 // =====================================================================
594 // ---------------------------------------------------------------------
595 // ---------------------------------------------------------------------
596 // =====================================================================
597 
598 vEvent* vDataFile::GetEvent(int64_t a_eventIndex, int a_th)
599 {
601  // The event pointer that will be returned
602  vEvent* p_event;
603  // Compute the adress into the buffer
604  char* event_address_in_array = mp_mappedMemory + a_eventIndex*m_sizeEvent;
605  // Get the event
606  p_event = GetEventSpecific(event_address_in_array, a_th);
607  // Return the event
608  return p_event;
609 }
610 
611 
612 
613 
614 // =====================================================================
615 // ---------------------------------------------------------------------
616 // ---------------------------------------------------------------------
617 // =====================================================================
618 
620  // dev tool, but it can be used currently by the datafile
621  // Ignored for now in windows compilation
622 
623 int vDataFile::Shuffle( int64_t nb_events_to_load )
624 {
625  #ifndef _WIN32
626  if (m_verbose >=2) Cout("vDataFile::Shuffle() : Everyday I shuffle in... " << endl);
627  // Buffer storing ordered indices from ( 0 ) to ( nb_events_to_load - 1 )
628  int64_t *rndmIdx = new int64_t[ nb_events_to_load ];
629  std::iota( rndmIdx, rndmIdx + nb_events_to_load, 0 );
630 
631  // Initializing random
632  std::random_device rd;
633  std::mt19937 rndm( rd() );
634  rndm.seed( 1100001001 );
635 
636  // Shuffling the buffer of indices
637  std::shuffle( rndmIdx, rndmIdx + nb_events_to_load, rndm );
638 
639  // Creating a tmp buffer
640  char *mp_arrayEvents_tmp = new char[ nb_events_to_load*m_sizeEvent ];
641 
642  // Copy sorted buffer to shuffled buffer
643  for( int64_t i = 0; i < nb_events_to_load; ++i )
644  {
645  for( int64_t j = 0; j < m_sizeEvent; ++j )
646  {
647  //mp_arrayEvents_tmp[ i * m_sizeEvent + j ] = mp_arrayEvents[ rndmIdx[ i ] * m_sizeEvent + j ];
648  mp_arrayEvents_tmp[ i * m_sizeEvent + j ] = mp_mappedMemory[ rndmIdx[ i ] * m_sizeEvent + j ];
649  }
650  }
651 
652  // Freeing memory
653  delete[] rndmIdx;
654  //delete[] mp_mappedMemory;
655 
656  //mp_arrayEvents = mp_arrayEvents_tmp;
657  mp_mappedMemory = mp_arrayEvents_tmp;
658 
659  if (m_verbose >=2) Cout("vDataFile::Shuffle OK" << endl);
660  #endif
661  return 0;
662 }
663 
664 
665 
666 // =====================================================================
667 // ---------------------------------------------------------------------
668 // ---------------------------------------------------------------------
669 // =====================================================================
670 
671 void vDataFile::GetEventIndexStartAndStop( int64_t* ap_indexStart, int64_t* ap_indexStop, int a_subsetNum, int a_nbSubsets )
672 {
674  // Basically here, the index start for a single MPI instance will be the current subset number.
675  // If multiple instances are used, the whole datafile is split into equal-size concatenated pieces.
676  // So for each instance, we just have to found the first index falling in its range (assuming that
677  // the event index step is always equal to the number of subsets).
678 
679  // For the first machine, index start is the current subset number
680  if (mp_ID->GetMPIRank()==0) *ap_indexStart = a_subsetNum;
681  // For the other machines, we must search for the first index falling in their range belonging to this
682  // subset (a subset being starting at a_subsetNum with a step equal to the number of subsets, a_nbSubsets)
683  else
684  {
685  // Compute the modulo of the first index of this machine minus the subset number with respect to the number of subsets
686  int64_t modulo = (m_mpi1stEvent-a_subsetNum) % a_nbSubsets;
687  // If this modulo is null, then the index start is the first index
688  if (modulo==0) *ap_indexStart = m_mpi1stEvent;
689  // Otherwise, the index start is equal to the first index plus the number of subsets minus the modulo
690  else *ap_indexStart = m_mpi1stEvent + (a_nbSubsets - modulo);
691  }
692 
693  // For index stop, we simply get the last event of the MPI instance (+1 because the for loop is exclusive)
694  *ap_indexStop = m_mpiLastEvent + 1;
695 }
696 
697 // =====================================================================
698 // ---------------------------------------------------------------------
699 // ---------------------------------------------------------------------
700 // =====================================================================
701 
703 {
705  Cerr("*****vDataFile::GetMaxRingDiff() -> This function is not implemented for the used system" << endl);
706  Cerr(" (this error may be prompted if the present function is erroneously called for a SPECT system)" << endl);
707  return -1;
708 }
709 
710 // =====================================================================
711 // ---------------------------------------------------------------------
712 // ---------------------------------------------------------------------
713 // =====================================================================
714 
716 {
717  // Verbose
719 
720  // Close all multithreaded datafiles before merging
721  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
722  if (m2p_dataFile[th]) m2p_dataFile[th]->close();
723 
724  // If the output file doesn't exist yet, simply rename the first temporary file name using the ouput file name
725  if (!ifstream(m_dataFileName))
726  {
727  // If only one projection is required (i.e one threads && no projection of frame or gates)...
728  if(mp_ID->GetNbThreadsForProjection()==1 // No multithreading so no multiple tmp datafiles
729  && mp_ID->GetNbTimeFrames()
731  * mp_ID->GetNb2ndMotImgsForLMS() == 1)
732  {
733  // rename the first temporary file name to the global name
734  string tmp_file_name = m_dataFileName + "_0";
735 
736  // ... then just rename the first temporary file name using the ouput file name
737  rename(tmp_file_name.c_str(),m_dataFileName.c_str());
738  // no need to concatenate, so we leave here.
739  return 0 ;
740  }
741  }
742 
743  // Create the final output file which will concatenate the events inside the thread-specific data files
744  ofstream merged_file(m_dataFileName.c_str(), ios::out | ios::binary | ios::app);
745 
746  // Concatenation : generate input file ("data_file") to read the buffer of the thread-specific data files and store the information in the final output file
747  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
748  {
749  // Build thread file name
750  stringstream ss; ss << th;
751  string file_name = m_dataFileName;
752  file_name.append("_").append(ss.str());
753  // Open it
754  // Note SS: Maybe we can use the m2p_dataFile[th] here by just rewarding to the beginning of the file ?
755  // Note TM: There were some issues involving the use of rdbuf and concatenation of the file in this case (ifstream were needed instead of fstream for some reasons)
756  // But we should have another look on how the projection data writing works with the implementation of the new analytical simulator.
757  ifstream data_file(file_name.c_str(), ios::binary | ios::in);
758 
759  if (!data_file)
760  {
761  Cerr(endl);
762  Cerr("***** vDataFile::PROJ_WriteData() -> Input temporary thread file '" << file_name << "' is missing or corrupted !" << endl);
763  return 1;
764  }
765 
766  // Concatenate it to the merged file
767  merged_file << data_file.rdbuf();
768  // Close file
769  data_file.close();
770 
771  // Re-open datafiles (needed if projecting frame/gates, as the contents of the temporary datafile are copied to the main datafile after each frame/gate)
772  m2p_dataFile[th]->open( file_name.c_str(), ios::binary | ios::out | ios::trunc);
773  }
774 
775  // Close merged file
776  merged_file.close();
777 
778  return 0;
779 }
780 
781 // =====================================================================
782 // ---------------------------------------------------------------------
783 // ---------------------------------------------------------------------
784 // =====================================================================
785 
787 {
788  // Verbose
790  if (m_verbose>=VERBOSE_DETAIL) Cout("vDataFile::PROJ_DeleteTmpDataFile() ..." << endl);
791  // Generate input file ("data_file") to read the buffer of the thread-specific data files and store the information in the final output file
792  for (int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
793  {
794  // Build thread file name
795  stringstream ss; ss << th;
796 
797  if (m2p_dataFile[th]) m2p_dataFile[th]->close();
798 
799  string file_name = m_dataFileName;
800  file_name.append("_").append(ss.str());
801 
802  // Remove temporary file for data output writing (Projection script only)
803  ifstream fcheck(file_name.c_str());
804  if(fcheck.good())
805  {
806  fcheck.close();
807  #ifdef _WIN32
808  string dos_instruction = "del " + file_name;
809  system(dos_instruction.c_str());
810  #else
811  remove(file_name.c_str());
812  #endif
813  }
814  }
815  // End
816  return 0;
817 }
818 
819 // =====================================================================
820 // ---------------------------------------------------------------------
821 // ---------------------------------------------------------------------
822 // =====================================================================
823 
824 vEvent* vDataFile::PROJ_GenerateEvent(int a_idxElt1, int a_idxElt2, int a_th)
825 {
827  // Only 1 line required for projection/sensitivity generation
828  m2p_BufferEvent[a_th]->SetNbLines(1);
829  m2p_BufferEvent[a_th]->SetID1(0, a_idxElt1);
830  m2p_BufferEvent[a_th]->SetID2(0, a_idxElt2);
831  // Return the event
832  return m2p_BufferEvent[a_th];
833 }
834 
835 // =====================================================================
836 // ---------------------------------------------------------------------
837 // ---------------------------------------------------------------------
838 // =====================================================================
839 
841 {
842  // Simply switch between the different data types
843  string data_type = "";
844  if (m_dataType==TYPE_CT) data_type = "CT";
845  else if (m_dataType==TYPE_PET) data_type = "PET";
846  else if (m_dataType==TYPE_SPECT) data_type = "SPECT";
847  else data_type = "unknown";
848  return data_type;
849 }
850 
851 // =====================================================================
852 // ---------------------------------------------------------------------
853 // ---------------------------------------------------------------------
854 // =====================================================================
855 
857 {
858  // Simply switch between the different data modes
859  string data_mode = "";
860  if (m_dataMode==MODE_HISTOGRAM) data_mode = "histogram";
861  else if (m_dataMode==MODE_LIST) data_mode = "list-mode";
862  else if (m_dataMode==MODE_NORMALIZATION) data_mode = "normalization";
863  else data_mode = "unknown";
864  return data_mode;
865 }
866 
867 // =====================================================================
868 // ---------------------------------------------------------------------
869 // ---------------------------------------------------------------------
870 // =====================================================================
871 
873 {
874  // Simply switch between the different data specs
875  string data_spec = "";
876  if (m_dataSpec==SPEC_EMISSION) data_spec = "emission";
877  else if (m_dataSpec==SPEC_TRANSMISSION) data_spec = "transmission";
878  else data_spec = "unknown";
879  return data_spec;
880 }
881 
882 // =====================================================================
883 // ---------------------------------------------------------------------
884 // ---------------------------------------------------------------------
885 // =====================================================================
int64_t GetEventSize()
Definition: vDataFile.hh:339
#define KEYWORD_OPTIONAL_SUCCESS
Definition: gOptions.hh:54
This class is designed to be a mother virtual class for DataFile.
Definition: vDataFile.hh:103
void GetEventIndexStartAndStop(int64_t *ap_indexStart, int64_t *ap_indexStop, int a_subsetNum=0, int a_NbSubsets=1)
Compute the index start and stop of the events loop with respect to the current subset and MPI size a...
Definition: vDataFile.cc:671
#define VERBOSE_DEBUG_EVENT
#define MODE_HISTOGRAM
Definition: vDataFile.hh:59
#define TYPE_PET
Definition: vDataFile.hh:74
void SetNbLines(uint16_t a_value)
Set the number of lines of the Event.
Definition: vEvent.hh:154
int64_t m_sizeEvent
Definition: vDataFile.hh:598
bool mp_POIDirectionFlag[3]
Definition: vDataFile.hh:594
string m_dataFileName
Definition: vDataFile.hh:578
#define TYPE_UNKNOWN
Definition: vDataFile.hh:72
#define MODE_LIST
Definition: vDataFile.hh:57
virtual vEvent * GetEventSpecific(char *ap_buffer, int a_th)=0
This function is implemented in child classes Read an event from the position pointed by 'ap_buffer...
#define MODE_NORMALIZATION
Definition: vDataFile.hh:61
int CheckParameters()
Check the initialization of member variables Call the CheckSpecificParameters() function implemente...
Definition: vDataFile.cc:315
FLTNB GetCalibrationFactor()
Definition: vDataFile.hh:368
bool GetPOIInfoFlag()
Definition: vDataFile.hh:386
int GetNb1stMotImgsForLMS(int a_fr)
call the eponym function from the oDynamicDataManager object
virtual int CheckFileSizeConsistency()=0
int SetParametersFrom(vDataFile *ap_DataFile)
Initialize all parameters from the provided datafile.
Definition: vDataFile.cc:275
oMemoryMapped * mp_MappedFile
Definition: vDataFile.hh:604
int ReadInfoInHeader(bool a_affectQuantificationFlag=true)
Read and check general information from the header datafile Call the ReadSpecificInformationInHeade...
Definition: vDataFile.cc:110
#define VERBOSE_DETAIL
int SetCalibrationFactor(int a_bed, FLTNB a_calibrationFactor)
Set the calibration factor for the provided bed.
FLTNB m_durationInSec
Definition: vDataFile.hh:584
virtual int CheckSpecificConsistencyWithAnotherDataFile(vDataFile *ap_DataFile)=0
Check consistency between 'this' and the provided datafile, for specific characteristics.
int64_t GetSize()
Definition: vDataFile.hh:333
string m_headerFileName
Definition: vDataFile.hh:577
int CheckConsistencyWithAnotherBedDataFile(vDataFile *ap_DataFile)
Check consistency between 'this' and the provided datafile as two bed positions.
Definition: vDataFile.cc:386
string GetDataSpecToString()
Definition: vDataFile.cc:872
int m_dataMode
Definition: vDataFile.hh:580
static sOutputManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
#define TYPE_SPECT
Definition: vDataFile.hh:76
#define TYPE_CT
Definition: vDataFile.hh:78
virtual ~vDataFile()
vDataFile destructor.
Definition: vDataFile.cc:85
virtual int GetMaxRingDiff()
Return an error by default. This function is surcharged by the PET (and CT) scanner daughter class...
Definition: vDataFile.cc:702
int m_dataSpec
Definition: vDataFile.hh:582
FLTNB m_startTimeInSec
Definition: vDataFile.hh:583
int Open(const std::string &filename, size_t mappedBytes=WholeFile, CacheHint hint=Normal)
open file, mappedBytes = 0 maps the whole file
int64_t m_mpi1stEvent
Definition: vDataFile.hh:601
int InitializeMappedFile()
Check the datafile existency, map it to memory and get the raw char* pointer. .
Definition: vDataFile.cc:462
int PROJ_WriteData()
Write/Merge chunk of data in a general data file.
Definition: vDataFile.cc:715
vEvent * GetEvent(int64_t a_eventIndex, int a_th=0)
Definition: vDataFile.cc:598
int CloseFile()
Close as many binary file stream for writing.
Definition: vDataFile.cc:578
int64_t m_mpiLastEvent
Definition: vDataFile.hh:602
#define Cerr(MESSAGE)
const string & GetPathName()
#define VERBOSE_DEBUG_LIGHT
bool m_ignorePOIFlag
Definition: vDataFile.hh:593
virtual int CheckSpecificParameters()=0
This function is implemented in child classes Check specific parameters of child classes...
bool * GetPOIDirectionFlag()
Definition: vDataFile.hh:380
Singleton class that manages output writing on disk (images, sinograms, etc). It also manages loggi...
char * mp_mappedMemory
Definition: vDataFile.hh:605
vDataFile()
vDataFile constructor.
Definition: vDataFile.cc:39
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 GetMPISize()
Get the MPI size (the number of MPI instances)
int m_dataType
Definition: vDataFile.hh:581
everything ... be careful when file is larger than memory
bool GetBedPositionFlag()
Definition: vDataFile.hh:426
int64_t m_nbEvents
Definition: vDataFile.hh:579
oImageDimensionsAndQuantification * mp_ID
Definition: vDataFile.hh:573
virtual int ReadSpecificInfoInHeader(bool a_affectQuantificationFlag=true)=0
This function is implemented in child classes Read and check modality-specific information from the...
fstream ** m2p_dataFile
Definition: vDataFile.hh:599
FLTNB m_relativeBedPosition
Definition: vDataFile.hh:587
#define VERBOSE_NORMAL
FLTNB * GetPOIResolution()
Definition: vDataFile.hh:374
int GetDataType()
Definition: vDataFile.hh:311
int GetBedIndex()
Definition: vDataFile.hh:294
int OpenFileForWriting(string a_suffix="")
Open a binary file stream for writing, with eventually the suffix appended to the file name...
Definition: vDataFile.cc:537
string GetScannerName()
Definition: vDataFile.hh:522
vEvent * PROJ_GenerateEvent(int idx_elt1, int idx_elt2, int a_th)
Generate a standard event and set up its ID Used by the projection, list-mode sensitivity generatio...
Definition: vDataFile.cc:824
#define SPEC_UNKNOWN
Definition: vDataFile.hh:89
int m_bedIndex
Definition: vDataFile.hh:586
virtual void DescribeSpecific()=0
A pure virtual function used to describe the specific parts of the datafile.
#define KEYWORD_MANDATORY
Definition: gOptions.hh:48
#define KEYWORD_OPTIONAL
Definition: gOptions.hh:50
#define SPEC_TRANSMISSION
Definition: vDataFile.hh:93
const unsigned char * GetData() const
raw access
FLTNB mp_POIResolution[3]
Definition: vDataFile.hh:595
void SetID2(int a_line, uint32_t a_value)
Set the indice associated with the line index for the 2nd ID of the Event.
Definition: vEvent.hh:170
#define OS_SEP
bool m_bedPositionFlag
Definition: vDataFile.hh:588
const string & GetBaseName()
int GetDataSpec()
Definition: vDataFile.hh:322
Mother class for the Event objects.
Definition: vEvent.hh:43
int GetDataMode()
Definition: vDataFile.hh:300
int GetNbTimeFrames()
Get the number of time frames.
vEvent ** m2p_BufferEvent
Definition: vDataFile.hh:600
FLTNB GetRelativeBedPosition()
Definition: vDataFile.hh:432
int GetNb2ndMotImgsForLMS()
call the eponym function from the oDynamicDataManager object
bool m_POIInfoFlag
Definition: vDataFile.hh:592
#define SPEC_EMISSION
Definition: vDataFile.hh:91
void Describe()
A function used to describe the generic parts of the datafile.
Definition: vDataFile.cc:510
int SetAcquisitionTime(int a_bed, FLTNB a_timeStartInSec, FLTNB a_durationInSec)
Set the acquisition time if not already set by the SetTimeFrames()
int GetNbThreadsForProjection()
Get the number of threads used for projections.
string GetDataTypeToString()
Definition: vDataFile.cc:840
int64_t m_mpiNbEvents
Definition: vDataFile.hh:603
#define MODE_UNKNOWN
Definition: vDataFile.hh:55
#define DEBUG_VERBOSE(IGNORED1, IGNORED2)
virtual int SetSpecificParametersFrom(vDataFile *ap_DataFile)=0
This function is implemented in child classes Initialize all specific parameters from the provided ...
#define Cout(MESSAGE)
string m_scannerName
Definition: vDataFile.hh:589
Portable read-only memory mapping (Windows and Linux)
FLTNB m_calibrationFactor
Definition: vDataFile.hh:585
int PROJ_DeleteTmpDataFile()
Delete temporary datafile used for multithreaded output writing if needed.
Definition: vDataFile.cc:786
int GetVerbose()
Get the verbose level.
Definition: vDataFile.hh:446
virtual int Shuffle(int64_t)
!!!\ This function has been modified to be used specifically with a
Definition: vDataFile.cc:623
int64_t GetDuration()
Definition: vDataFile.hh:362
int m_verbose
Definition: vDataFile.hh:574
int GetMPIRank()
Get the MPI instance number (the rank)
string GetDataModeToString()
Definition: vDataFile.cc:856
good overall performance
int64_t GetStartTime()
Definition: vDataFile.hh:356
#define VERBOSE_DEBUG_NORMAL
void SetID1(int a_line, uint32_t a_value)
Set the indice associated with the line index for the 1st ID of the Event.
Definition: vEvent.hh:162