castor-mergeBedPositions.cc

Go to the documentation of this file.

#include "gVariables.hh"
#include "gOptions.hh"
#include "oIterativeAlgorithm.hh"
#include "oSensitivityGenerator.hh"
#include "oImageDimensionsAndQuantification.hh"
#include "iDataFilePET.hh"
#include "iDataFileSPECT.hh"
#include "sOutputManager.hh"
#include "sScannerManager.hh"
#include "sRandomNumberGenerator.hh"
#include "iScannerPET.hh"
#include "sAddonManager.hh"

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


void ShowHelp()
{
  // Show help
  cout << endl;
  cout << "Usage:  castor-mergeBedPositions  -img image1.hdr  -img image2.hdr  -(f/d)out output  [settings]" << endl;
  cout << endl;
  cout << "This program can be used to merge several images corresponding to multiple bed positions of the same acquisition. The slices of any image" << endl;
  cout << "are weighted during merging. By default, these weights are uniform (1 for all slices of all images). The -sens option can be used to provide" << endl;
  cout << "a sensitivity image for each image to be merged, where the voxels' value of the sensitivity images are used as weights." << endl;
  cout << "The bed displacement between two successive bed positions is recovered from the image headers as 'multiple bed displacement'." << endl;
  cout << "If not found, the default one from the scanner is used." << endl;
  cout << "Important note: this program assumes the provided images exactly match the axial FOV of the scanner; if bigger or smaller results will be wrong." << endl;
  cout << endl;
  cout << "[Mandatory parameters]:" << endl;
  cout << "  -img imageX.hdr      : Give an image corresponding to a bed position (must give at least two images; the provided order defines the order" << endl;
  cout << "                         in which bed positions are merged)" << 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 << "  -sens sensX.hdr      : Give the sensitivity image to weight the imageX.hdr while merging the bed positions (must provide as many sensitivity" << endl;
  cout << "                         images as images to be merged; the association of sensitivity images to the images is done with respect to the order" << endl;
  cout << "                         both images are provided)" << endl;
  cout << "  -oweight             : Also save the global weights" << endl;
  cout << "  -inv                 : Virtually inverse the order of the provided images (and sensitivity images)" << endl;
  cout << "  -vb                  : Give the verbosity level, from 0 (no verbose) to 2 (default: 1)" << endl;
  cout << "  -conf                : Give the path to the CASToR config directory (default: located through the CASTOR_CONFIG environment variable)." << 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
  // --------------------------------------------------------------------------------

  // Number of bed position 
  uint32_t nb_beds = 0;
  // Vector containing string pointing to the images to be merged
  vector<string> path_to_image_filename;
  // Vector containing string pointing to the sensitivity images used as weights
  vector<string> path_to_sens_filename;
  // Invert bed datafile names order
  bool invert_order_flag = false;

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

  // Output directory name.
  string path_dout = "";
  // Or root name
  string path_fout = "";
  // Save global weight
  bool save_weights = false;

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

  // Verbose level
  int verbose = 1;
  // Path to config directory
  string path_to_config_dir = "";

  // ============================================================================================================
  // 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-mergeBedPositions() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      if (ConvertFromString(argv[i+1], &verbose))
      {
        Cerr("***** castor-mergeBedPositions() -> Exception when trying to read provided verbosity level '" << verbose << " for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      i++;
    }
    // Path to config directory
    else if (option=="-conf")
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-mergeBedPositions() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      path_to_config_dir = (string)argv[i+1];
      i++;
    }

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

    // Images
    else if (option=="-img") // This is a mandatory option
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-mergeBedPositions() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      string file_name = (string)argv[i+1];
      path_to_image_filename.push_back(file_name);
      nb_beds++;
      i++;
    }
    // Sensitivity images
    else if (option=="-sens")
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-mergeBedPositions() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      string file_name = (string)argv[i+1];
      path_to_sens_filename.push_back(file_name);
      i++;
    }
    // Invert bed positions order
    else if (option=="-inv")
    {
      invert_order_flag = true;
    }

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

    // Name of the output directory
    else if (option=="-dout") // This is a mandatory option
    {
      if (i>=argc-1)
      {
        Cerr("***** castor-mergeBedPositions() -> 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-mergeBedPositions() -> Argument missing for option: " << option << endl);
        Exit(EXIT_FAILURE);
      }
      path_fout = argv[i+1];
      i++;
    }
    // Save global weights
    else if (option=="-oweight")
    {
      save_weights = true;
    }

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

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

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

  // Data files
  if (nb_beds < 2)
  {
    Cerr("***** castor-mergeBedPositions() -> Please provide at least two image filename !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Output files
  if (path_fout.empty() && path_dout.empty())
  {
    Cerr("***** castor-mergeBedPositions() -> 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-mergeBedPositions() -> Please provide either output option -fout or -dout but not both !" << endl);
    Exit(EXIT_FAILURE);
  }
  // If sensitivity is provided, then check that there are as many files as input images
  if (path_to_sens_filename.size()>0 && path_to_sens_filename.size()!=path_to_image_filename.size())
  {
    Cerr("***** castor-mergeBedPositions() -> If using sensitivity images as weights, please provide as many sensitivity images as input images !" << 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);
  // Set path to the config directory
  if (p_OutputManager->CheckConfigDir(path_to_config_dir))
  {
    Cerr("***** castor-mergeBedPositions() -> A problem occured while checking for the config directory path !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Initialize output directory and base name
  if (p_OutputManager->InitOutputDirectory(path_fout, path_dout))
  {
    Cerr("***** castor-mergeBedPositions() -> A problem occured while initializing output directory !" << endl);
    Exit(EXIT_FAILURE);
  }
  // Log command line
  if (p_OutputManager->LogCommandLine(argc,argv))
  {
    Cerr("***** castor-mergeBedPositions() -> A problem occured while logging command line arguments !" << endl);
    Exit(EXIT_FAILURE);
  }

  // ============================================================================================================
  // Create sScannerManager
  // ============================================================================================================
  sScannerManager* p_ScannerManager = sScannerManager::GetInstance();  
  p_ScannerManager->SetVerbose(verbose);

  // ============================================================================================================
  // Read all input images
  // ============================================================================================================

  // Verbose
  if (verbose>=1) Cout("castor-mergeBedPositions() -> Load input images" << endl);

  // Get user endianness (interfile I/O)
  GetUserEndianness();

  // Create interfile image fields and image pointers
  Intf_fields** p_image_fields = (Intf_fields**)malloc(nb_beds*sizeof(Intf_fields*));
  FLTNB** p_images = (FLTNB**)malloc(nb_beds*sizeof(FLTNB*));

  // Read interfile images
  for (uint32_t bed=0; bed<nb_beds; bed++)
  {
    // Create the image fields
    p_image_fields[bed] = new Intf_fields;
    // Read the image
    if (invert_order_flag) p_images[bed] = IntfLoadImageFromScratch(path_to_image_filename[nb_beds-1-bed],p_image_fields[bed],verbose);
    else p_images[bed] = IntfLoadImageFromScratch(path_to_image_filename[bed],p_image_fields[bed],verbose);
    // Check reading
    if (p_images[bed]==NULL)
    {
      Cerr("***** castor-mergeBedPositions() -> A problem occured while reading image for bed " << bed+1 << " !" << endl);
      return 1;
    }
  }

  // Check dimensions consistency between all image fields, as well as originating system and bed displacement
  for (uint32_t bed=1; bed<nb_beds; bed++)
  {
    if (IntfCheckDimensionsConsistency(*p_image_fields[0],*p_image_fields[bed]))
    {
      Cerr("***** castor-mergeBedPositions() -> Inconsistency between image dimensions !" << endl);
      Exit(EXIT_FAILURE);
    }
  }

  // Get dimensions locally
  uint32_t dim_x = p_image_fields[0]->mtx_size[0];
  uint32_t dim_y = p_image_fields[0]->mtx_size[1];
  uint32_t dim_z = p_image_fields[0]->mtx_size[2];
  uint32_t dim_xyz = dim_x * dim_y * dim_z;

  // ============================================================================================================
  // Recover the bed displacement if not present in the image headers
  // ============================================================================================================

  // The bed displacement
  FLTNB multiple_bed_displacement = p_image_fields[0]->multiple_bed_displacement;

  // If not found in the interfile header, then try into the scanner configuration file
  if (multiple_bed_displacement<=0.)
  {
    // Verbose
    if (verbose>=1) Cout("castor-mergeBedPositions() -> Bed displacement not found in image headers, search for it in the scanner configuration file" << endl);
    // Find scanner configuration file
    if (p_ScannerManager->FindScannerSystem(p_image_fields[0]->originating_system) )
    {
      Cerr("***** castor-mergeBedPositions() -> A problem occurred while searching for scanner system !" << endl);
      Exit(EXIT_FAILURE);
    }
    // Get the bed displacement from the scanner configuration file
    if (ReadDataASCIIFile(p_ScannerManager->GetPathToScannerFile(), "multiple bed displacement", &multiple_bed_displacement, 1, KEYWORD_MANDATORY))
    {
      Cerr("***** castor-mergeBedPositions() -> Multiple bed displacement field not found in the scanner configuration file '" <<
           p_ScannerManager->GetPathToScannerFile() << " !" << endl);
      Exit(EXIT_FAILURE);
    }
  }

  // Verbose
  if (verbose>=1) Cout("castor-mergeBedPositions() -> Use a bed displacement of " << multiple_bed_displacement << " mm" << endl);

  // ============================================================================================================
  // Check that the bed displacement is compatible with the slice thickness
  // ============================================================================================================

  // Define the tolerance associated to the rounding
  FLTNB tolerance = 1.e-4;
  // Compute floating point number of slices
  FLTNB nb_slices_fltnb = multiple_bed_displacement / p_image_fields[0]->vox_size[2];
  // Compute integer number of slices
  FLTNB nb_slices_intnb = ((FLTNB)((uint32_t)nb_slices_fltnb));
  // The number of slices defining the displacement
  uint32_t displacement_slices = 0;
  // Case where the floating number minus the integer number is less than 0.5; this means that
  // the closest integer to the floating number of slices is below it.
  if (nb_slices_fltnb-nb_slices_intnb<0.5)
  {
    // If the difference is below the tolerance, we are ok taking the 'integerized' number of slices as the displacement
    if (nb_slices_fltnb-nb_slices_intnb<tolerance) displacement_slices = ((uint32_t)nb_slices_fltnb);
    // Otherwise, the displacement is not compatible with the slice thickness
    else
    {
      Cerr("***** castor-mergeBedPositions() -> Bed displacement is not compatible with the slice thickness of " << p_image_fields[0]->vox_size[2] << " mm !" << endl);
      Exit(EXIT_FAILURE);
    }
  }
  // Case where the floating number minus the integer number is more than 0.5; this means that
  // the closest integer to the floating number of slices is above it.
  else
  {
    // If the minus the difference + 1 is below the tolerance, we are of taking the 'integerized' number of slices + 1 as the displacement
    if (nb_slices_intnb+1.-nb_slices_fltnb<tolerance) displacement_slices = ((uint32_t)nb_slices_fltnb) + 1;
    // Otherwise, the displacement is not compatible with the slice thickness
    else
    {
      Cerr("***** castor-mergeBedPositions() -> Bed displacement is not compatible with the slice thickness of " << p_image_fields[0]->vox_size[2] << " mm !" << endl);
      Exit(EXIT_FAILURE);
    }
  }

  // Still check that the displacement slices makes sense
  if (displacement_slices==0 || displacement_slices>dim_z)
  {
    Cerr("***** castor-mergeBedPositions() -> Computed displacement of " << displacement_slices << " slices makes no sense !" << endl);
    Exit(EXIT_FAILURE);
  }

  // Verbose
  if (verbose>=1) Cout("castor-mergeBedPositions() -> Bed displacement corresponds to " << displacement_slices << " slices" << endl);

  // Compute the number of slices (this is one full image + the number of beds - 1 times the displacement)
  uint32_t whole_dim_z = dim_z + ((uint32_t)(nb_beds-1)) * displacement_slices;
  uint32_t whole_dim_xyz = dim_x * dim_y * whole_dim_z;

  // Verbose
  if (verbose>=1) Cout("castor-mergeBedPositiosn() -> Output merged image has " << whole_dim_z << " slices" << endl);

  // ============================================================================================================
  // If sensitivity, then load it as weights, otherwise set everything to 1
  // ============================================================================================================

  // Create the weights matrix
  FLTNB** p_weights = (FLTNB**)malloc(nb_beds*sizeof(FLTNB*));
  for (uint32_t bed=0; bed<nb_beds; bed++) p_weights[bed] = NULL;

  // If we have sensitivity
  if (path_to_sens_filename.size()>0)
  {
    // Verbose
    if (verbose>=1) Cout("castor-mergeBedPositions() -> Load sensitivity images" << endl);
    // Create interfile image fields and image pointers for sensitivity
    Intf_fields** p_sens_fields = (Intf_fields**)malloc(nb_beds*sizeof(Intf_fields*));
    // Read interfile images
    for (uint32_t bed=0; bed<nb_beds; bed++)
    {
      // Create the image fields
      p_sens_fields[bed] = new Intf_fields;
      // Read the image
      if (invert_order_flag) p_weights[bed] = IntfLoadImageFromScratch(path_to_sens_filename[nb_beds-1-bed],p_sens_fields[bed],verbose);
      else p_weights[bed] = IntfLoadImageFromScratch(path_to_sens_filename[bed],p_sens_fields[bed],verbose);
      // Check reading
      if (p_weights[bed]==NULL)
      {
        Cerr("***** castor-mergeBedPositions() -> A problem occured while reading sensitivity for bed " << bed+1 << " !" << endl);
        return 1;
      }
    }
    // Check dimensions consistency between all image fields, as well as originating system and bed displacement
    for (uint32_t bed=1; bed<nb_beds; bed++)
    {
      if (IntfCheckDimensionsConsistency(*p_sens_fields[0],*p_sens_fields[bed]))
      {
        Cerr("***** castor-mergeBedPositions() -> Inconsistency between sensitivity dimensions !" << endl);
        Exit(EXIT_FAILURE);
      }
    }
    // Check dimensions consistency, etc, between input images and sensitivity images
    for (uint32_t bed=0; bed<nb_beds; bed++)
    {
      if (IntfCheckDimensionsConsistency(*p_image_fields[bed],*p_sens_fields[bed]))
      {
        Cerr("***** castor-mergeBedPositions() -> Inconsistency between input image and sensitivity dimensions for bed " << bed+1 << " !" << endl);
        Exit(EXIT_FAILURE);
      }
    }
  }
  // If we do not have sensitivity
  else
  {
    // Verbose
    if (verbose>=1) Cout("castor-mergeBedPositions() -> Initialize weights to 1." << endl);
    // We set everything to 1.
    for (uint32_t bed=0; bed<nb_beds; bed++)
    {
      p_weights[bed] = (FLTNB*)malloc(dim_xyz*sizeof(FLTNB));
      for (uint32_t v=0; v<dim_xyz; v++) p_weights[bed][v] = 1.;
    }
  }

  // ============================================================================================================
  // Compute the whole body image
  // ============================================================================================================

  // Verbose
  if (verbose>=1) Cout("castor-mergeBedPositiosn() -> Compute whole image" << endl);
  // Allocate whole image
  FLTNB* p_whole_image = (FLTNB*)malloc(whole_dim_xyz*sizeof(FLTNB));
  for (uint32_t v=0; v<whole_dim_xyz; v++) p_whole_image[v] = 0.;
  // Allocate whole weights
  FLTNB* p_whole_weight = (FLTNB*)malloc(whole_dim_xyz*sizeof(FLTNB));
  for (uint32_t v=0; v<whole_dim_xyz; v++) p_whole_weight[v] = 0.;

  // Loop on all beds
  for (uint32_t bed=0; bed<nb_beds; bed++)
  {
    // Loop on all slices
    for (uint32_t z=0; z<dim_z; z++)
    {
      // Compute whole z index
      uint32_t whole_z = z + bed * displacement_slices;
      // For efficiency
      uint32_t base_z = z * dim_x * dim_y;
      uint32_t whole_base_z = whole_z * dim_x * dim_y;
      // Loop on all voxels inside the slice
      for (uint32_t xy=0; xy<dim_x*dim_y; xy++)
      {
        // Add the contribution of this bed to the whole image
        p_whole_image[whole_base_z+xy] += p_images[bed][base_z+xy] * p_weights[bed][base_z+xy];
        // Add the weights to the global weights
        p_whole_weight[whole_base_z+xy] += p_weights[bed][base_z+xy];
      }
    }
  }
  // Normalize the whole image with respect to the weights
  for (uint32_t v=0; v<whole_dim_xyz; v++)
  {
    if (p_whole_weight[v]>0.) p_whole_image[v] /= p_whole_weight[v];
    else p_whole_image[v] = 0.;
  }

  // ============================================================================================================
  // Save output image
  // ============================================================================================================

  // Verbose
  if (verbose>=1) Cout("castor-mergeBedPositions() -> Write output image(s)" << endl);

  // Build interfile fields (we use first one of the input images and modify the relevant fields)
  // Set header metadata using Image Dimensions object
  p_image_fields[0]->mtx_size[2] = whole_dim_z;
  p_image_fields[0]->nb_total_imgs = whole_dim_z;
  p_image_fields[0]->nb_bytes_pixel = sizeof(FLTNB);
  /* TODO: UPDATE TIME START AND DURATION BUT WILL HAVE TO WRITE THEM IN THE HEADERS OF IMAGES WRITTEN DURING ITERATIONS
  ap_IF->study_duration = ap_ID->GetFinalTimeStopInSec(0) -
                          ap_ID->GetFrameTimeStartInSec(0,0);
  for(int fr=0 ; fr<ap_ID->GetNbTimeFrames() ; fr++)
  {
    ap_IF->image_duration.push_back(ap_ID->GetFrameDurationInSec(0, fr));
    ap_IF->frame_group_pause.push_back((fr == 0) ? 0 
                                                 : ap_ID->GetFrameTimeStartInSec(0,fr) - ap_ID->GetFrameTimeStopInSec(0,fr-1));
  }
  */
  // File name
  string whole_image_filename = p_OutputManager->GetPathName() + p_OutputManager->GetBaseName() + "_whole";
  // Write the image
  if (IntfWriteImageFromIntfFields(whole_image_filename, p_whole_image, *p_image_fields[0], verbose))
  {
    Cerr("***** castor-mergeBedPositions() -> A problem occured while writing output whole image !" << endl);
    Exit(EXIT_FAILURE);
  }

  // Save global weights if asked for
  if (save_weights)
  {
    // File name
    string whole_weight_filename = p_OutputManager->GetPathName() + p_OutputManager->GetBaseName() + "_weights";
    // Write the image
    if (IntfWriteImageFromIntfFields(whole_weight_filename, p_whole_weight, *p_image_fields[0], verbose))
    {
      Cerr("***** castor-mergeBedPositions() -> A problem occured while writing global weights image !" << endl);
      Exit(EXIT_FAILURE);
    }
  }

  // ============================================================================================================
  // Free memory properly
  // ============================================================================================================

  // Input images
  if (p_images)
  {
    for (uint32_t bed=0; bed<nb_beds; bed++) if (p_images[bed]) free(p_images[bed]);
    free(p_images);
  }
  // Weights
  if (p_weights)
  {
    for (uint32_t bed=0; bed<nb_beds; bed++) if (p_weights[bed]) free(p_weights[bed]);
    free(p_weights);
  }
  // Whole image
  if (p_whole_image) free(p_whole_image);
  // Whole image
  if (p_whole_weight) free(p_whole_weight);

  // Ending
  if (verbose>=1) Cout(endl);
  #ifdef CASTOR_MPI
  MPI_Finalize();
  #endif
  return EXIT_SUCCESS;
}