iDeformationElastic.cc

Go to the documentation of this file.

#include "iDeformationElastic.hh"

/*
  \fn public iDeformationElastic::iDeformationElastic()
  \brief Constructor of iDeformationElastic. Simply set all data members to default values.
*/
iDeformationElastic::iDeformationElastic() : vDeformation()
{
  m_pathToFwdDeformationFiles.clear();
  m_pathToBwdDeformationFiles.clear();

  //file_parser = ParserType::New();
  #ifdef CASTOR_ELASTIX
  ITK_inputImage  = ImageType::New();
  ITK_outputImage = ImageType::New();
  #endif
}



/*
  \fn public iDeformationElastic::iDeformationElastic()
  \brief Destructor of iDeformationElastic. Free memory from all allocated tabs.
*/
iDeformationElastic::~iDeformationElastic() 
{
  if( m_initialized)
  {
    m_pathToFwdDeformationFiles.clear();
    m_pathToBwdDeformationFiles.clear();
  }
}



/*
  \fn public void ShowHelp()
  \brief This function is an implementation of the pure virtual mother function. It is used to
         print out specific help about the deformation model and its options.
*/
void iDeformationElastic::ShowHelp()
{
  cout << "-- Motion correction based on non-rigid deformations calculated using the Elastix software (v.4.8) --" << endl;
  cout << "For the current implementation, a text file containing the filename of the Elastix deformations is required:" << endl; 
  cout << "The 1st line indicates the path to the forward transformation parameters, while the second line indicates the path to backward transformation parameters." << endl;
  cout << "The path should contain '%d' instead of the number of gates in order to adequately load the parameter files. Numbering should start from 1 (not 0)" << endl;
  cout << "(ex 'pathTo/parameterFiles_fileNumber.xxx' becomes 'pathTo/parameterFiles_%d.xxx' )" << endl;
}


/*
  \fn public virtual int ReadAndCheckConfigurationFile()
  \param const string& a_configurationFile
  \brief This function is an implementation of the pure virtual mother function. It is used to read options from a configuration file.
  \return 0 if success, other value otherwise.
*/
int iDeformationElastic::ReadAndCheckConfigurationFile(const string& a_fileOptions)
{
  // Implement here the reading of any options specific to this deformation model (i.e : parameters or path to deformation files), through a configuration file
  // The ReadDataASCIIFile() functions could be helpful to recover data from a file
  
  ifstream in_file(a_fileOptions.c_str(), ios::in);
  
  if(!in_file)
  {
    Cerr("***** iDeformationElastic::ReadAndCheckConfigurationFile -> Error while trying to read configuration file : " << a_fileOptions << endl);
    return 1;
  }
  
  string line;
  char* str_tmp;
  
  // Get forward files basename
  getline(in_file, line);
  str_tmp = new char[line.size()];

  for(int t=0 ; t<m_nbTransformations ; t++)
  {
    //TODO : error handling (incorrect string)
    sprintf(str_tmp, line.c_str(), t+1);
    m_pathToFwdDeformationFiles.push_back(str_tmp);
  }

  // Get backward files basename
  getline(in_file, line);
  
  delete str_tmp;
  str_tmp = new char[line.size()];

  for(int t=0 ; t<m_nbTransformations ; t++)
  {
    //TODO : error handling (incorrect string)
    sprintf(str_tmp, line.c_str(), t+1);
    //sprintf(str_tmp, line.c_str(), 1);
    m_pathToBwdDeformationFiles.push_back(str_tmp);
  }
  
  delete str_tmp;
  in_file.close();
  
  // Check if files exist
  for(int t=0 ; t<m_nbTransformations ; t++)
  {
    ifstream in_file_fwd(m_pathToFwdDeformationFiles.at(t).c_str(), ios::in);
    ifstream in_file_bwd(m_pathToFwdDeformationFiles.at(t).c_str(), ios::in);
    if(!in_file_fwd)
    {
      Cerr("***** iDeformationElastic::ReadAndCheckConfigurationFile -> Error while trying to read forward deformation vectors file number" << t <<  ": " 
           << m_pathToFwdDeformationFiles.at(t) << endl);
      return 1;
    }
    if(!in_file_bwd)
    {
      Cerr("***** iDeformationElastic::ReadAndCheckConfigurationFile -> Error while trying to read backward deformation vectors file number" << t <<  ": " 
           << m_pathToBwdDeformationFiles.at(t) << endl);
      return 1;
    }
    in_file_fwd.close();
    in_file_bwd.close();
  }

  return 0;
}


/*
  \fn public virtual int ReadAndCheckOptionsList()
  \param const string& a_optionsList
  \brief This function is an implementation of the pure virtual mother function. It is used to read options from a list of options.
  \return 0 if success, other value otherwise.
*/
int iDeformationElastic::ReadAndCheckOptionsList(const string& a_listOptions)
{
  // Implement here the reading of any options specific to this deformation model, through a list of options separated by commas
  // The ReadStringOption() function could be helpful to parse the list of parameters in an array
  
  // Normal end
  return 0;
}



/*
  \fn CheckSpecificParameters
  \brief This function is an implementation of the pure virtual mother function. It is used to
         check parameters of the child deformation model before initialization.
  \return 0 if success, other value otherwise.
*/
int iDeformationElastic::CheckSpecificParameters()
{
  // Implement here checks over parameters which should be read using either ReadAndCheckConfigurationFile() and ReadAndCheckOptionsList() functions
  
  
  // All the variables checked here are member variables of the mother deformation class 
  
  // Check image dimensions
  if (mp_ID==NULL)
  {
    Cerr("***** iDeformationElastic::CheckParameters() -> No image dimensions provided !" << endl);
    return 1;
  }
  
  // Check resp gates
  if (m_nbTransformations<0)
  {
    Cerr("***** iDeformationElastic::CheckParameters() -> Wrong number of gates !" << endl);
    return 1;
  }
  
  // Check verbosity
  if (m_verbose<0)
  {
    Cerr("***** iDeformationElastic::CheckParameters() -> Wrong verbosity level provided !" << endl);
    return 1;
  }

  // Check vector files
  if (m_pathToFwdDeformationFiles.empty() || m_pathToBwdDeformationFiles.empty())
  {
    Cerr("***** iDeformationElastic::CheckParameters() -> Vector files not initialized !" << endl);
    return 1;
  }
  
  // Normal end
  m_checked = true;
  return 0;
}


/*
  \fn public int Initialize()
  \brief This function is an implementation of the pure virtual mother function. It is used to
         initialize specific stuff to the child deformation model.
  \return 0 if success, other value otherwise.
*/
int iDeformationElastic::Initialize()
{
  // Implement here the allocation/initialization of whatever member variables specifically used by this deformation model
  
  
  // Check here if Castor has been compiled with CASTOR_ELASTIX option. Throw error otherwise.
  #ifndef CASTOR_ELASTIX
    Cerr("***** iDeformationElastic::Initialize() -> This specific deformation class should be compiled with the CASTOR_ELASTIX flag on in CMake !" << endl);
    Cerr("*****                                      It requires ITK and the elastix library !" << endl);
    return 1;
  #endif


  if (!m_checked)
  {
    Cerr("***** iDeformationElastic::Initialize() -> Parameters should be checked before Initialize() !" << endl);
    return 1;
  }
  
  // Normal end
  m_initialized = true;
  return 0;
}





/*
  \fn public virtual int  iDeformationElastic::ApplyDeformations()
  \param FLTNB* ap_inputImage : input image to deform
  \param FLTNB* ap_outputImage : image in which the output of the deformation should be recovered
  \param FLTNB* a_direction : a direction for the deformation to perform (forward or backward)
  \param int a_defIdx : index of the deformation
  \brief This function is an implementation of the pure virtual mother function. The actual deformation should be implemented here 
  \return 0 if success, other value otherwise.
*/
int iDeformationElastic::ApplyDeformations(FLTNB* ap_inputImage, FLTNB* ap_outputImage, int a_direction, int a_defIdx)
{
  // The deformation model should be implemented here, with the help of any private functions if required
  
  /* The 'a_defIdx' parameter defines the deformation index of the transformation
  * The 'a_direction' parameter is an integer which indicates the direction of the deformation to perform, i.e :
  * - FORWARD_DEFORMATION (from the reference position to the 'a_defIdx' position) 
  * - BACKWARD_DEFORMATION (from the 'a_defIdx' position to the reference position) 
  * The integers FORWARD_DEFORMATION & BACKWARD_DEFORMATION are macros defined in the beginning of vDeformation.hh
  */
  
  // The deformation should be applied to theap_inputImage matrix, and the resulting image should be recovered in ap_outputImage matrix

  /* IMAGE DIMENSIONS :
  * For code efficiency and readability, the spatial index of a voxel is a cumulative 1D index. That is to say, given a voxel [indexX,indexY,indexZ],
  * its cumulative 1D index is computed by 'index = indexZ*nbVoxXY + indexY*nbVoxX + indexX'.
  *
  * The image dimensions can be recovered from the mp_ID class
  * Total number of voxels         : mp_ID->GetNbVoxXYZ()
  * Number of voxels in a slice    : mp_ID->GetNbVoxXY()
  * Number of voxels on the X-axis : mp_ID->GetNbVoxX()
  */
    
  string path_to_deformation_file;
  
  if (a_direction==FORWARD_DEFORMATION) path_to_deformation_file = m_pathToFwdDeformationFiles.at(a_defIdx);
  else if (a_direction==BACKWARD_DEFORMATION) path_to_deformation_file = m_pathToBwdDeformationFiles.at(a_defIdx);
  else
  {
    Cerr("*****iDeformationElastic::ApplyDeformationToImage()->Error, unknown direction type( must be 'forward' or 'backward' # " << endl);
    return 1;
  }
   
  
  #ifdef CASTOR_ELASTIX

  //Convert the input image from the binary to the ITK format

  m_RegionSizeInputImage[0] = mp_ID->GetNbVoxX();
  m_RegionSizeInputImage[1] = mp_ID->GetNbVoxY();
  m_RegionSizeInputImage[2] = mp_ID->GetNbVoxZ();

  m_RegionStartInputImage[0] = 0;
  m_RegionStartInputImage[1] = 0;
  m_RegionStartInputImage[2] = 0;

  m_RegionInputImage.SetSize( m_RegionSizeInputImage );
  m_RegionInputImage.SetIndex( m_RegionStartInputImage );

  ITK_inputImage->SetRegions(m_RegionInputImage);
  ITK_inputImage->Allocate();
  ITK_inputImage->FillBuffer(0);

  m_RegionSpacingInputImage[0] = mp_ID->GetVoxSizeX();
  m_RegionSpacingInputImage[1] = mp_ID->GetVoxSizeY();
  m_RegionSpacingInputImage[2] = mp_ID->GetVoxSizeZ();

  ITK_inputImage->SetSpacing( m_RegionSpacingInputImage );

  
  for(int i = 0; i< mp_ID->GetNbVoxX(); i++)
   for(int j = 0; j< mp_ID->GetNbVoxY(); j++)
    for(int k = 0; k< mp_ID->GetNbVoxZ(); k++)
    {
      m_PixelIndex[0] = i;
      m_PixelIndex[1] = j;
      m_PixelIndex[2] = k;
     
      ITK_inputImage->SetPixel( m_PixelIndex, (float)ap_inputImage[i + j*mp_ID->GetNbVoxX() + k*mp_ID->GetNbVoxXY()] );
    }

  /*WriterType::Pointer writer = WriterType::New();
  writer->SetFileName( "/home/fred/ProjetC/ElasticRegistration_Elastix-v4.8/test_correctionmouvement/test.mhd" );
  writer->SetInput( ITK_inputImage );
  try
      {
      writer->Update();
      }
    catch( itk::ExceptionObject & err )
      {
      std::cerr << "ExceptionObject caught !" << std::endl;
      std::cerr << err << std::endl;
      return EXIT_FAILURE;
      }  */


  // Try parsing transform parameters text file.
  ParserType::Pointer file_parser = ParserType::New();

  file_parser->SetParameterFileName( path_to_deformation_file );
  
  try
      {
  file_parser->ReadParameterFile();
  }
  catch( itk::ExceptionObject & e )
  {
  std::cout << e.what() << std::endl;
  
  // Do some error handling!
  Cout("iDeformationElastic::ApplyDeformationToImage()-> ERROR Reading deformation coefficients # " << a_defIdx << ", direction: "<< a_direction << endl);
  return 1;
  }

  // Retrieve parameter settings as map.
  RegistrationParametersType transform_parameters = file_parser->GetParameterMap();

  // Create the Fixed Image

   /* ImageType::Pointer Fixed_Image = ImageType::New();

    string FixedImageName = "/home/fred/ProjetC/ElasticRegistration_Elastix-v4.8/test_correctionmouvement/test_gate1_recon_it7.mhd";
    cout << "Loading Fixed Image: " << "/home/fred/ProjetC/ElasticRegistration_Elastix-v4.8/test_correctionmouvement/test_gate1_recon_it7.mhd" << " ...\n";

    ReaderType::Pointer reader_fixed = ReaderType::New();

    reader_fixed->SetFileName( FixedImageName );
    reader_fixed->Update();

    Fixed_Image = reader_fixed->GetOutput();*/


  // We are ready to call the warping function  //
  TRANSFORMIX* transformix = new TRANSFORMIX();
  int error = 0;

  string results_dir = "/tmp";
  
  if(m_verbose>1)
      {
  Cout("iDeformationElastic::ApplyDeformationToImage()->Apply coefficients # " << a_defIdx << ", direction: "<< a_direction << endl);
  Cout("iDeformationElastic::Reading_coefficients # " << a_defIdx << ": "<< path_to_deformation_file << endl<<endl);
      }

  try
     {
  error = transformix->TransformImage(
    static_cast<typename itk::DataObject::Pointer>( ITK_inputImage.GetPointer() ),
          //static_cast<typename itk::DataObject::Pointer>( Fixed_Image.GetPointer() ),
     transform_parameters,    // Parameter map read in previous code
    results_dir,
          false,    // boolean write_log_file, Enable/disable writing of elastix.log
    false     // boolean output_to_console, Enable/disable output to console
    );
     }
     catch( itk::ExceptionObject &err )
     {  
    // Do some error handling.
    Cout("iDeformationElastic:: ApplyDeformationToImage()->ERROR Applying coefficients # " << a_defIdx << ", direction: "<< a_direction << endl);
    return 1;
     }

  if( error == 0 )
     {
  if( transformix->GetResultImage().IsNotNull() )
        {
    ITK_outputImage = static_cast<ImageType *>( transformix->GetResultImage().GetPointer() );

    //Convert the Output image from the ITK format to the binary format
    for(int i = 0; i< mp_ID->GetNbVoxX(); i++)
           for(int j = 0; j< mp_ID->GetNbVoxY(); j++)
            for(int k = 0; k< mp_ID->GetNbVoxZ(); k++)
       {
         m_PixelIndex[0] = i;
         m_PixelIndex[1] = j;
         m_PixelIndex[2] = k;

         if (ITK_outputImage->GetPixel( m_PixelIndex ) < 0.0) ap_outputImage[i + j*mp_ID->GetNbVoxX() + k*mp_ID->GetNbVoxXY()] = 0.0;
         else ap_outputImage[i + j*mp_ID->GetNbVoxX() + k*mp_ID->GetNbVoxXY()] = ITK_outputImage->GetPixel( m_PixelIndex );
    //ap_outputImage[i + j*mp_ID->GetNbVoxX() + k*mp_ID->GetNbVoxXY()] = ITK_inputImage->GetPixel( m_PixelIndex );
       }
        }
      else
      {
    // Registration failure. Do some error handling.
    Cout("iDeformationElastic:: ApplyDeformationToImage()->ERROR in saving the warped image after applying coefficients # " << a_defIdx << ", direction: "<< a_direction << endl);
    return 1;
        }
     }

  // Typedef the ITKImageType first...
  /*stringstream ss;
  ss << a_defIdx;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName( "/home/fred/ProjetC/ElasticRegistration_Elastix-v4.8/test_correctionmouvement/test_warping_"+ss.str()+"_"+a_direction+".mhd" );
  writer->SetInput( ITK_outputImage );
  try
      {
      writer->Update();
      }
    catch( itk::ExceptionObject & err )
      {
      std::cerr << "ExceptionObject caught !" << std::endl;
      std::cerr << err << std::endl;
      return EXIT_FAILURE;
      }*/
  

  // Clean up memory.
  delete transformix;

  #endif
  
  
  
  // Any error should return a value >0.
  
  // Normal end
  return 0;
}