castor-norm.cc

Go to the documentation of this file.

#include "gVariables.hh"
#include "gOptions.hh"
#include "iEventHistoPET.hh"
#include "iEventNorm.hh"
#include "iProjectorClassicSiddon.hh"
#include "oImageDimensionsAndQuantification.hh"
#include "oInterfileIO.hh"
#include "oProjectionLine.hh"
#include "oProjectorManager.hh"
#include "vDataFile.hh"
#include "iDataFilePET.hh"
#include "iDataFileSPECT.hh"
#include "iDataFileCT.hh"
#include "sOutputManager.hh"
#include "sRandomNumberGenerator.hh"
#include <cstdint>

// =============================================================================================================================================
// =============================================================================================================================================
// =============================================================================================================================================
//                                                        H E L P     F U N C T I O N S
// =============================================================================================================================================
// =============================================================================================================================================
// =============================================================================================================================================


void ShowHelp()
{
  // Show help
  cout << endl;
  cout << "Usage:  castor-norm -df datafile.cdh -img path_to_img.hdr -sc scanner_alias -(f/d)out output" << endl;
  cout << endl;
  cout << "This program can be used to compute direct normalization factors from a normalization phantom." << endl;
  cout << endl;
  cout << "Direct normalization factors are used to counteract on the variation of sensitivity across all LORs of a PET system.";
  cout << " They are computed from a normalization scan (-df), typically of a symmetric and uniform phantom (-img), which illuminates as evenly as possible the different possible LORs.";
  cout << " The normalization factors are given by" << endl;
  cout << "    NF_uivj=A/C_uivj" << endl;
  cout << "where A is the activity of the uniform phantom, and C_uivj is the event count measured for LOR uivj.";
  cout << " The program outputs a normalization datafile that contains the computed NFs." << endl;
  cout << endl;
  cout << "[Mandatory parameters]:" << endl;
  cout << "  -df datafile.cdh     : Give a CASToR histogram datafile of a normalization scan." << endl;
  cout << "  -img path_to_img.hdr : Give the interfile header of the normalization phantom (no default)" << endl;
  cout << "  -sc  scanner_alias   : Give the scanner model. It must correspond to the name of one of the file in the scanner repository (w/o file extension)" << endl;
  cout << "  -fout name           : Give the root name for all output files (no default, alternative to -dout)" << endl;
  cout << "  -dout name           : Give the name of the output directory where all output files will be written (no default, alternative to -fout)" << endl;
  cout << endl;
  cout << "[Options]:" << endl;
  cout << "  -proj  param         : Give the projector to be used for forward projection, along with a configuration file (proj:file.conf)" << endl;
  cout << "  -proj-common         : Give common projector-related options." << endl;
  cout << "                         or the list of parameters associated to the projector (proj,param1,param2,...). If the projector only is specified, the" << endl;
  cout << "                         default configuration file is used. By default, the Siddon projector is used." << endl;
  cout << "  -proj-comp           : Give the strategy for projection line computation. Here are the three different strategies that can be used:" << endl;
  cout << "                     1 : Image-based system matrix elements storage: The voxels weights are added in a matrix representing the whole image, so" << endl;
  cout << "                         the addition of a new line to the previous ones is straightforward only by adding the weights to the corresponding voxels." << endl;
  cout << "                         As it is insanely long, it can possibly be used for example with extremely complex projectors that makes use of huge number" << endl;
  cout << "                         of ray tracings for a single event, where the list of contributions can become longer than the number of voxels in the image." << endl;
  cout << "                         This strategy is not compatible with SPECT reconstruction including attenuation correction." << endl;
  cout << "                     2 : Fixed-size list storage of system matrix elements: The voxels are added one by one in two separated lists, one containing voxel" << endl;
  cout << "                         indices and the other voxel weights. When a voxel is added to the oProjectionLine, it is simply pilled-up to the list. The list" << endl;
  cout << "                         has a fixed size which is provided by the EstimateMaxNumberOfVoxelsPerLine() function from the vProjector class. There are no" << endl;
  cout << "                         checks at all for possible buffer overflows. This is the fastest strategy and default one." << endl;
  cout << "                     3 : Adaptative-size list storage of system matrix elements: This is the same as the fixed-size strategy except that the size can be" << endl;
  cout << "                         upgraded if the current number of contributing voxels exceed the list's size. The first allocated size corresponds to the diagonal" << endl;
  cout << "                         of the image. During the first iteration, this size will be upgraded until it will reach a size suitable for all events. Thus it" << endl;
  cout << "                         is a bit slower than the fixed-list strategy during the first iteration, but is optimal with respect to RAM usage." << endl;
  cout << "  -vb value            : Give the verbosity level, from 0 (no verbose) to 2 (default: 1)" << endl;
  cout << "  --help,-h,-help      : Print out this help page." << endl; // managed by main
  cout << endl;
  #ifdef BUILD_DATE
  cout << "  Build date: " << BUILD_DATE << endl;
  cout << endl;
  #endif
  #ifdef CASTOR_VERSION
  cout << "  This program is part of the CASToR release version " << CASTOR_VERSION << "." << endl;
  cout << endl;
  #endif
}

// =============================================================================================================================================
// =============================================================================================================================================
// =============================================================================================================================================
//                                                        M A I N     P R O G R A M
// =============================================================================================================================================
// =============================================================================================================================================
// =============================================================================================================================================

int main(int argc, char** argv)
{
  // ============================================================================================================
  // MPI stuff (we make all instances except the first one returning 0 directly)
  // ============================================================================================================
  int mpi_rank = 0;
  #ifdef CASTOR_MPI
  int mpi_size = 1;
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
  if (mpi_rank!=0) return 0;
  #endif

  // No argument, then show help
  if (argc==1)
  {
    ShowHelp();
    Exit(EXIT_SUCCESS);
  }

  // ============================================================================================================
  // Parameterized variables with their default values
  // ============================================================================================================

  // --------------------------------------------------------------------------------
  // Input settings
  // --------------------------------------------------------------------------------

  // Path to the normalization datafile
  string path_to_normalization_datafile = "";

  // Path to the normalization phantom. no default
  string path_to_normalization_phantom = "";

  // Scanner name provided by the user
  string scanner_name = "";

  // --------------------------------------------------------------------------------
  // Output settings
  // --------------------------------------------------------------------------------

  // Output directory name.
  string path_dout = "";
  // Or root name
  string path_fout = "";

  // --------------------------------------------------------------------------------
  // Projector settings
  // --------------------------------------------------------------------------------

  string options_projectorF = "incrementalSiddon";
  string options_projector_common = "3.,1,1";

  // Default projector computation strategy
  int projector_computation_strategy = FIXED_LIST_COMPUTATION_STRATEGY;

  // --------------------------------------------------------------------------------
  // Miscellaneous settings
  // --------------------------------------------------------------------------------

  // Verbose level
  int verbose = 1;

  // ============================================================================================================
  // Read command-line parameters
  // ============================================================================================================

  // Must manually increment the option index when an argument is needed after an option
  for (int i=1; i<argc; i++)
  {
    // Get the option as a string
    string option = (string)argv[i];

    // --------------------------------------------------------------------------------
    // Miscellaneous settings
    // --------------------------------------------------------------------------------

    // Show help
    if (option=="-h" || option=="--help" || option=="-help")
    {
      ShowHelp();
      Exit(EXIT_SUCCESS);
    }
    // General verbosity level
    else if (option=="-vb")
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-norm() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      if (ConvertFromString(argv[i+1], &verbose))
      {
        Cerr("***** castor-norm() -> Exception when trying to read provided verbosity level '" << verbose << " for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      i++;
    }

    // --------------------------------------------------------------------------------
    // Input settings
    // --------------------------------------------------------------------------------

    // Images
    else if (option=="-df") // This is a mandatory option
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-norm() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      path_to_normalization_datafile = (string)argv[i+1];
      i++;
    }
    // Normalization scan
    else if (option=="-img")
    {
    if (i>=argc-1)
    {
     Cerr("***** castor-norm() -> Argument missing for option: " << option << endl);
     Exit(EXIT_FAILURE);
    }
    path_to_normalization_phantom = argv[i+1];
    i++;
    }
    // Scanner name
    else if (option=="-sc")
    {
    if (i>=argc-1)
    {
     Cerr("***** castor-norm() -> Argument missing for option: " << option << endl);
     Exit(EXIT_FAILURE);
    }
    scanner_name = argv[i+1];
    i++;
    }

    // --------------------------------------------------------------------------------
    // Output settings
    // --------------------------------------------------------------------------------

    // Name of the output directory
    else if (option=="-dout") // This is a mandatory option
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-norm() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      path_dout = argv[i+1];
      i++;
    }
    // Base name of the output files
    else if (option=="-fout") // This is a mandatory option
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-norm() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      path_fout = argv[i+1];
      i++;
    }

    // --------------------------------------------------------------------------------
    // Projector settings
    // --------------------------------------------------------------------------------

    // Projector settings
    else if (option=="-proj")
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-norm() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      options_projectorF = (string)argv[i+1];
      i++;
    }
    // Common projector settings
    else if (option=="-proj-common")
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-norm() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      options_projector_common = (string)argv[i+1];
      i++;
    }
    // Projection line computation strategy
    else if (option=="-proj-comp")
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-recon() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      projector_computation_strategy = atoi(argv[i+1]);
      i++;
    }

    // --------------------------------------------------------------------------------
    // Unknown option!
    // --------------------------------------------------------------------------------

    else
    {
      Cerr("***** castor-norm() -> Unknown option '" << option << "' !" << endl);
      Exit(EXIT_FAILURE);
    }
  }

  // ============================================================================================================
  // Some checks
  // ============================================================================================================

  // Output files
  if (path_fout.empty() && path_dout.empty())
  {
    Cerr("***** castor-norm() -> Please provide an output option for output files (-fout or -dout) !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Check that only one option has been provided
  if (!path_fout.empty() && !path_dout.empty())
  {
    Cerr("***** castor-norm() -> Please provide either output option -fout or -dout but not both !" << endl);
    Exit(EXIT_FAILURE);
  }

  // ============================================================================================================
  // Create sOutputManager
  // ============================================================================================================
  sOutputManager *p_OutputManager = sOutputManager::GetInstance();
  // Set verbose level
  p_OutputManager->SetVerbose(verbose);
  // Set MPI rank
  p_OutputManager->SetMPIRank(mpi_rank);
  // Initialize output directory and base name
  if (p_OutputManager->InitOutputDirectory(path_fout, path_dout))
  {
    Cerr("***** castor-norm() -> A problem occurred while initializing output directory !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Log command line
  if (p_OutputManager->LogCommandLine(argc,argv))
  {
    Cerr("***** castor-norm() -> A problem occurred while logging command line arguments !" << endl);
    Exit(EXIT_FAILURE);
  }

  // ============================================================================================================
  // Create sScannerManager
  // ============================================================================================================
  sScannerManager *p_ScannerManager = sScannerManager::GetInstance();
  p_ScannerManager->SetVerbose(verbose);
  if (p_ScannerManager->FindScannerSystem(scanner_name) )
  {
    Cerr("***** castor-norm() -> A problem occurred while searching for scanner system !" << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_ScannerManager->BuildScannerObject() )
  {
    Cerr("***** castor-norm() -> A problem occurred during scanner object construction ! !" << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_ScannerManager->InstantiateScanner() )
  {
    Cerr("***** castor-norm() -> A problem occurred while creating Scanner object !" << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_ScannerManager->BuildLUT() )
  {
    Cerr("***** castor-norm() -> A problem occurred while generating/reading the LUT !" << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_ScannerManager->CheckParameters())
  {
    Cerr("***** castor-norm() -> A problem occurred while checking scanner manager parameters !" << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_ScannerManager->Initialize())
  {
    Cerr("***** castor-norm() -> A problem occurred while initializing scanner !" << endl);
    Exit(EXIT_FAILURE);
  }

  // Get the number of normalization factors
  int64_t nb_elems = p_ScannerManager->GetSystemNbElts();

  // ============================================================================================================
  // Recover image dimensions informations from the interfile header of the image to project
  // ============================================================================================================

  // (Required for options using interfile I/O)
  GetUserEndianness();

  Intf_fields IF;
  IntfKeyInitFields(&IF);

  if(IntfReadHeader(path_to_normalization_phantom, &IF, verbose) )
  {
    Cerr("***** castor-norm() -> An error occurred while trying to read the "
         "interfile header of file reading file "
         << path_to_normalization_phantom << " !" << endl);
    Exit(EXIT_FAILURE);
  }

  if(verbose >= 3) IntfKeyPrintFields(IF);

  int nb_VoxX = IF.mtx_size[0];
  int nb_VoxY = IF.mtx_size[1];
  int nb_VoxZ = IF.mtx_size[2];

  FLTNB vox_SizeX = IF.vox_size[0];
  FLTNB vox_SizeY = IF.vox_size[1];
  FLTNB vox_SizeZ = IF.vox_size[2];

  FLTNB offsetX = IF.vox_offset[0];
  FLTNB offsetY = IF.vox_offset[1];
  FLTNB offsetZ = IF.vox_offset[2];

  int nb_time_basis = IF.nb_time_frames > 1 ? IF.nb_time_frames : 1;
  int nb_resp_basis = IF.nb_resp_gates > 1 ? IF.nb_resp_gates : 1;
  int nb_card_basis = IF.nb_card_gates > 1 ? IF.nb_card_gates : 1;

  if (verbose>=2)
  {
    Cout(" Image dimensions for projections: " << endl);
    Cout(" voxels number (X,Y,Z): " << nb_VoxX << "," << nb_VoxY << "," << nb_VoxZ << endl);
    Cout(" voxels size (X,Y,Z): " << vox_SizeX << "," << vox_SizeY << "," << vox_SizeZ << endl);
    Cout(" offset (X,Y,Z): " << offsetX << "," << offsetY << "," << offsetZ << endl);
  }

    // ============================================================================================================
  // oImageDimensionsAndQuantification settings
  // ============================================================================================================

  // Create oImageDimensionsAndQuantification
  oImageDimensionsAndQuantification* p_ImageDimensionsAndQuantification = new oImageDimensionsAndQuantification();
  p_ImageDimensionsAndQuantification->SetDefault();

  p_ImageDimensionsAndQuantification->SetNbVoxX(nb_VoxX);
  p_ImageDimensionsAndQuantification->SetNbVoxY(nb_VoxY);
  p_ImageDimensionsAndQuantification->SetNbVoxZ(nb_VoxZ);
  p_ImageDimensionsAndQuantification->SetVoxSizeX(vox_SizeX);
  p_ImageDimensionsAndQuantification->SetVoxSizeY(vox_SizeY);
  p_ImageDimensionsAndQuantification->SetVoxSizeZ(vox_SizeZ);
  p_ImageDimensionsAndQuantification->SetOffsetX(offsetX);
  p_ImageDimensionsAndQuantification->SetOffsetY(offsetY);
  p_ImageDimensionsAndQuantification->SetOffsetZ(offsetZ);

  p_ImageDimensionsAndQuantification->SetNbTimeBasisFunctions(nb_time_basis);
  p_ImageDimensionsAndQuantification->SetNbRespBasisFunctions(nb_resp_basis);
  p_ImageDimensionsAndQuantification->SetNbCardBasisFunctions(nb_card_basis);

  p_ImageDimensionsAndQuantification->SetNbBeds(1);

  // Threading
  int nb_threads = 1;
  #ifdef CASTOR_OMP
  nb_threads = omp_get_num_procs();
  if (verbose>=2) Cout("  --> Use " << nb_threads << " threads" << endl);
  #endif
  if(p_ImageDimensionsAndQuantification->SetNbThreads(to_string(nb_threads)))
  {
    Cerr("***** castor-norm() -> A problem occurred while setting the number of threads !" << endl);
    Exit(EXIT_FAILURE);
  }

  if (p_ImageDimensionsAndQuantification->CheckParameters())
  {
    Cerr("***** castor-norm() -> A problem occurred while checking image dimensions parameters !" << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_ImageDimensionsAndQuantification->Initialize())
  {
    Cerr("***** castor-norm() -> A problem occurred while initializing image dimensions !" << endl);
    Exit(EXIT_FAILURE);
  }

  // ============================================================================================================
  // Set image space and initialise image
  // ============================================================================================================

  // --- Image space initialization ---
  oImageSpace* p_ImageSpace = new oImageSpace();

  p_ImageSpace->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
  p_ImageSpace->SetVerbose(verbose);

  // Allocate memory for main image
  p_ImageSpace->InstantiateImage();
  p_ImageSpace->InitImage(path_to_normalization_phantom, 0);

  // ============================================================================================================
  // Create the input datafile
  // ============================================================================================================

  // Compute the expected number of normalization factors
  int64_t nb_nf = (int) (nb_elems * (nb_elems - 1) / 2);

  // Create the datafile based on the data type
  iDataFilePET* p_DataFile = new iDataFilePET(); // do we want to normalize other modalities?
  p_DataFile->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
  p_DataFile->SetHeaderDataFileName(path_to_normalization_datafile);
  p_DataFile->SetIgnoreTOFFlag(true); // force to sum all TOF bins
  p_DataFile->SetVerbose(verbose);
  p_DataFile->SetBedIndex(0);
  bool do_not_affect_quantification = false;
  if (p_DataFile->ReadInfoInHeader(do_not_affect_quantification))
  {
    Cerr("***** castor-norm() -> A problem occurred during datafile header reading ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_DataFile->GetDataType() != TYPE_PET)
  {
    Cerr("***** castor-norm() -> Only PET data is currently supported ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_DataFile->GetDataMode() != MODE_HISTOGRAM)
  {
    Cerr("***** castor-norm() -> Only histogram data mode is currently supported. Aborting." << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_DataFile->CheckParameters())
  {
    Cerr("***** castor-norm() -> A problem occurred while checking datafile parameters ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_DataFile->ComputeSizeEvent())
  {
    Cerr("***** castor-norm() -> A problem occurred in datafile initialization ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_DataFile->GetSize() != nb_nf) {
    Cerr("***** castor-norm() -> The input datafile does not contain expected number of events (got" << p_DataFile->GetSize() << ", expected " << nb_nf << ") ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_DataFile->InitializeMappedFile())
  {
    Cerr("***** castor-norm() -> A problem occurred in datafile initialization ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_DataFile->PrepareDataFile())
  {
    Cerr("***** castor-norm() -> A problem occurred in datafile preparation ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }

  // ============================================================================================================
  // Create the output datafile
  // ============================================================================================================

  // Create output datafile object
  iDataFilePET* p_OutputDataFile = new iDataFilePET();
  p_OutputDataFile->SetVerbose(verbose);
  p_OutputDataFile->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
  p_OutputDataFile->SetNormCorrectionFlagOn();
  p_OutputDataFile->SetDataMode(MODE_NORMALIZATION);
  p_OutputDataFile->SetBedIndex(0);
  if (p_OutputDataFile->PROJ_InitFile())
  {
    Cerr("***** castor-norm() -> An error occurred while opening file for writing !" << endl);
    Exit(EXIT_FAILURE);
  }
  p_OutputDataFile->SetNbEvents(nb_nf);
  if (p_OutputDataFile->ComputeSizeEvent())
  {
    Cerr("***** castor-norm() -> A problem occurred in datafile initialization ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }
  if (p_OutputDataFile->PrepareDataFile())
  {
    Cerr("***** castor-norm() -> A problem occurred in datafile preparation ! Abort." << endl);
    Exit(EXIT_FAILURE);
  }

  // ============================================================================================================
  // Create Projector Manager
  // ============================================================================================================

  // Verbose
  if (verbose>=5) Cout("----- Projector initialization ... -----" << endl);
  // Create object
  oProjectorManager* p_ProjectorManager = new oProjectorManager();
  // Set all parameters
  p_ProjectorManager->SetScanner(p_ScannerManager->GetScannerObject());
  p_ProjectorManager->SetImageDimensionsAndQuantification(p_ImageDimensionsAndQuantification);
  p_ProjectorManager->SetDataFile(p_DataFile);
  p_ProjectorManager->SetComputationStrategy(projector_computation_strategy);
  p_ProjectorManager->SetOptionsForward(options_projectorF);
  p_ProjectorManager->SetOptionsBackward("incrementalSiddon"); // will not be used
  p_ProjectorManager->SetOptionsCommon(options_projector_common);
  p_ProjectorManager->SetVerbose(verbose);
  // Check parameters
  if (p_ProjectorManager->CheckParameters())
  {
    Cerr("***** castor-norm() -> A problem occurred while checking projector manager's parameters !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Initialize projector manager
  if (p_ProjectorManager->Initialize())
  {
    Cerr("***** castor-norm() -> A problem occurred while initializing projector manager !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Verbose
  if (verbose>=5) Cout("----- Projector initialization OK -----" << endl);

  // ============================================================================================================
  // Compute normalization factors
  // ============================================================================================================

  // Verbose
  if (verbose>=1) Cout("castor-norm() -> Compute normalization factors for " << p_DataFile->GetSize() << " events" << endl);

  // Variables used in loop
  int64_t printing_index = 0;
  bool problem = false; // used to report problems in the parallel loop

  // Number of normalization factors that were correctly computed
  int64_t nf_correctly_computed = 0;

  // Create normalization event objects
  iEventNorm** p2_EventNorm = new iEventNorm*[nb_threads];
  for (INTNB th = 0; th < nb_threads; th++)
  {
    iEventNorm* p_EventNorm = new iEventNorm();
    p_EventNorm->SetNbLines(1); // 1 = lines in event
    p_EventNorm->AllocateID();
    p2_EventNorm[th] = p_EventNorm;
  }

  // Compute start and stop indices
  int64_t index = 0;
  int64_t index_start = 0;
  int64_t index_stop  = 0;
  p_DataFile->GetEventIndexStartAndStop(&index_start, &index_stop);

  #pragma omp parallel for private(index) schedule(static), num_threads(nb_threads)
  for (index = index_start; index < index_stop; index++)
  {
    // The thread index
    int th = 0;
    #ifdef CASTOR_OMP
    th = omp_get_thread_num();
    #endif

    // Verbose
    if (th==0 && verbose>=2)
    {
      if (printing_index%((int)(index_stop/100/nb_threads))==0)
      {
        int percent = ((int)( ((FLTNB)(index)) * (FLTNB)nb_threads * 100. / index_stop ));
        cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b      "
             << percent << " %                    " << flush;
      }
      printing_index++;
    }
    // Read the event
    iEventHistoPET* p_EventHistoPET = (iEventHistoPET*) p_DataFile->GetEvent(index,th);
    // id1 always less than id2
    uint32_t id1 = p_EventHistoPET->GetID1(0);
    uint32_t id2 = p_EventHistoPET->GetID2(0);

    // Compute the normalization factor
    FLTNBDATA nf_index = 0.; // Normalization factor default value
    FLTNBDATA event_value = p_EventHistoPET->GetEventValue(0); // 0 = no TOF
    if (event_value > 0.){
      // Compute the projection line
      oProjectionLine* p_ProjectionLine = p_ProjectorManager->ComputeProjectionLine(p_EventHistoPET, th);
      if (p_ProjectionLine == NULL){
        Cerr("***** castor-norm() -> A problem occurred while computing the projection line !" << endl);
        // Specify that there was a problem
        problem = true;
        // We must continue here because we are inside an OpenMP loop
        continue;
      }
      FLTNB activity = 0.;
      if (p_ProjectionLine->NotEmptyLine())
      {
        activity = p_ProjectionLine->ForwardProject(p_ImageSpace->m4p_image[0][0][0]);
        if (activity > 0.)
        {
          nf_index = activity / event_value;
          nf_correctly_computed++;
          // Cout(id1 << "," << id2 << ": activity=" << activity << ", event_value=" << event_value << ", nf_index=" << nf_index << endl); // Debug
        } else if (activity < 0.) {
          Cerr("***** castor-norm() -> A negative projection has been computed !" << endl);
          Cerr("*****" << id1 << "," << id2 << ": activity=" << activity << ", event_value=" << event_value << endl);
          // Specify that there was a problem
          problem = true;
          // We must continue here because we are inside an OpenMP loop
          continue;
        }
      }
      // If this part of code is reached, then !p_ProjectionLine->NonEmptyLine() or activity == 0.
      // In any case, we set the normalization to 0.
      if (activity == 0. && verbose >= 3)
      {
          Cerr("!!!!! castor-norm() -> The projection of the line of response having ID1=" << id1 << " and ID2=" << id2 << " is null, but the corresponding event activity is non-null (" << event_value << ").");
          Cerr(" The normalization correction factor (" << nf_index << ") will be erroneous for this line of response." << endl);
      }
    }

    // Write the normalization factor
    iEventNorm* p_EventNorm = p2_EventNorm[th];
    p_EventNorm->SetID1(0, id1);
    p_EventNorm->SetID2(0, id2);
    p_EventNorm->SetNormalizationFactor(nf_index);
    p_OutputDataFile->WriteEvent(p_EventNorm, th);
  }

  // If a problem was encountered, then report it here
  if (problem)
  {
    Cerr("***** castor-norm() -> A problem occurred inside the parallel loop over events !" << endl);
    Exit(EXIT_FAILURE);
  }

  // Verbose
  if (verbose>=2)
  {
    cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
         << "  --> 100 %                        " << endl;
    Cout("  --> Compute normalization factors" << endl);
  }

  // Summary of normalisation
  if (verbose >= 1){
    float percent = 100. * nf_correctly_computed / nb_nf;
    Cout("castor-norm() -> " << nf_correctly_computed << " out of " << nb_nf << " (" << percent << "%) normalization factors correctly computed." << endl);
    Cout("Please be aware that a low percentage may result in poor reconstruction quality!" << endl);
  }

  // ============================================================================================================
  // Write final output
  // ============================================================================================================

  // Write final file
  if(p_OutputDataFile->PROJ_WriteData())
  {
    Cerr("***** castor-norm() -> An error occurred while writing the final file !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Delete temporary files
  if(p_OutputDataFile->PROJ_DeleteTmpDataFile())
  {
    Cerr("***** castor-norm() -> An error occurred while deleting temporary files !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Write header
  if (p_OutputDataFile->WriteHeader())
  {
    Cerr("***** castor-norm() -> An error occurred while writing output header file !" << endl);
    Exit(EXIT_FAILURE);
  }

  // ============================================================================================================
  // Exit
  // ============================================================================================================

  // Delete objects in the inverse order in which they were created
  if (verbose>=5) Cout("----- Deleting CASToR objects ... -----" << endl);

  for (INTNB th = 0; th < nb_threads; th++)
  {
    delete p2_EventNorm[th];
  }
  delete[] p2_EventNorm;

  delete p_ProjectorManager;
  delete p_OutputDataFile;
  delete p_DataFile;

  p_ImageSpace->DeallocateImage();
  delete p_ImageSpace;

  delete p_ImageDimensionsAndQuantification;
  if (verbose>=5) Cout("----- CASToR objects successfully deleted -----" << endl);

  // Ending
  #ifdef CASTOR_MPI
  MPI_Finalize();
  #endif
  return EXIT_SUCCESS;
}