iDeformationRigid.cc

Go to the documentation of this file.

/*
  Implementation of class iDeformationRigid

  - separators: X
  - doxygen: X
  - default initialization: X
  - CASTOR_DEBUG: X
  - CASTOR_VERBOSE: X
*/

#include "iDeformationRigid.hh"


// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn iDeformationRigid
  \brief Constructor of iDeformationRigid. Simply set all data members to default values.
*/
iDeformationRigid::iDeformationRigid() : vDeformation()
{
  m_pathToFwdDeformationFiles.clear();
  m_pathToBwdDeformationFiles.clear();
  m_OriginalArray = NULL;
  m_OutputArray = NULL;
  mp_vectorX = NULL;
  mp_vectorY = NULL;
  mp_vectorZ = NULL;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ~iDeformationRigid
  \brief Destructor of iDeformationRigid. Free memory from all allocated tabs.
*/
iDeformationRigid::~iDeformationRigid() 
{
  if(m_initialized)
  {
    m_pathToFwdDeformationFiles.clear();
    m_pathToBwdDeformationFiles.clear();
    if(m_OriginalArray) delete[] m_OriginalArray;
    if(m_OutputArray) delete[] m_OutputArray;  
    if(mp_vectorX) delete[] mp_vectorX; 
    if(mp_vectorY) delete[] mp_vectorY; 
    if(mp_vectorZ) delete[] mp_vectorZ; 
  }
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ShowHelp
  \brief This function is used to print out specific help about the deformation and its options.
*/
void iDeformationRigid::ShowHelp()
{
  // todo : really we should rather use param files in binary format here
  cout << "-- Rigid transformation algorithm based on trilinear interpolation --" << endl;
  cout << "It requires a set of ASCII transformation files containing the forward and backward deformation vectors" << endl;
  cout << endl;
  cout << "The path and basename of the transformation files must be provided through a main ASCII header file" <<endl;
  cout << "(1) First provide the path to the file using the -rm (respiratory) -cm (cardiac) or -im (patient) motion correction option :" << endl;
  cout << "-rm/cm/im deformationRigid:path_to/header_init_file.txt" << endl;
  cout << "The header must contain the keywords 'Forward vector files' and 'Backward vector files'," << endl;
  cout << "which provide the basename of the vector files related to the forward and backward deformation respectively" << endl;
  cout << "Ex :" << endl;
  cout << "--------------------------------------------------------" << endl;
  cout << "header_init_file.txt :" << endl;
  cout << "--------------------------------------------------------" << endl;
  cout << "Forward vector files  : pathTo/parameterFiles_Fwd_%d.xxx" << endl;
  cout << "Backward vector files : pathTo/parameterFiles_Bwd_%d.xxx" << endl;
  cout << "--------------------------------------------------------" << endl;
  cout << "In the header file, the gate index should be replaced by '%d' as presented in the exemple above." << endl;
  cout << "The Forward/Backward vector ASCII files must be ordered from 1 (not 0) to the last transformation index." << endl;
  cout << "For each voxel, the x,y,z components of the vector must be provided on a new line, as in the following template :" << endl;
  cout << "Ex :" << endl;
  cout << "--------------------------------------------------------" << endl;
  cout << "parameterFiles_Fwd_1.txt :" << endl;
  cout << "--------------------------------------------------------" << endl;
  cout << "vox0_x" << endl;
  cout << "vox0_y" << endl;
  cout << "vox0_z" << endl;
  cout << "vox1_x" << endl;
  cout << "vox1_y" << endl;
  cout << "vox1_z" << endl;
  cout << "ect ..." << endl;
  cout << "--------------------------------------------------------" << endl;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ReadAndCheckConfigurationFile
  \param a_configurationFile
  \brief This function is used to read options from a configuration file.
  \return 0 if success, other value otherwise.
*/
int iDeformationRigid::ReadAndCheckConfigurationFile(const string& a_fileOptions)
{
  if (m_verbose >=2) Cout("iDeformationRigid::ReadAndCheckConfigurationFile ..."<< endl); 
  
  // First, check the nb of transformations has correctly been initializeds
  if (m_nbTransformations<0)
  {
    Cerr("***** iDeformationRigid::ReadAndCheckConfigurationFile() -> Number of transformations in the deformation has not been initialized  !" << endl);
    return 1;
  }
  
  // Need to know the gate numbers before
  ifstream in_file(a_fileOptions.c_str(), ios::in);
  
  if(!in_file)
  {
    Cerr("***** iDeformationRigid::ReadAndCheckConfigurationFile -> Error while trying to read configuration file : " << a_fileOptions << endl);
    return 1;
  }
  
  string line;
  char* str_tmp;
  
  // --- Recover forward files basename (1st string in the file) --- //
  getline(in_file, line);
  str_tmp = new char[line.size()];

  for(int t=0 ; t<m_nbTransformations ; t++)
  {
    sprintf(str_tmp, line.c_str(), t+1);
    m_pathToFwdDeformationFiles.push_back(str_tmp);
  }


  // --- Recover backward files basename (2nd string in the file) --- //
  getline(in_file, line);
  
  delete str_tmp;
  str_tmp = new char[line.size()];

  for(int t=0 ; t<m_nbTransformations ; t++)
  {
    sprintf(str_tmp, line.c_str(), t+1);
    m_pathToBwdDeformationFiles.push_back(str_tmp);
  }
  
  
  // Delete local objects
  delete str_tmp;
  in_file.close();
  
  
  // --- Check file existence --- //
  
   // \todo Perhaps more checks on the consistency of the file. May do this in initialization though
  
  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_pathToBwdDeformationFiles.at(t).c_str(), ios::in);
    
    if(!in_file_fwd)
    {
      Cerr("***** iDeformationRigid::ReadAndCheckConfigurationFile -> Error while trying to read forward deformation vectors file number" << t <<  ": " 
           << m_pathToFwdDeformationFiles.at(t) << endl);
      return 1;
    }

    if(!in_file_bwd)
    {
      Cerr("***** iDeformationRigid::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 ReadAndCheckOptionsList
  \param a_optionsList
  \brief This function is used to read options from a list of options.
         Throw error by defaut for this method, as a file has to be used for initialization
  \return 0 if success, other value otherwise.
*/
int iDeformationRigid::ReadAndCheckOptionsList(const string& a_listOptions)
{
  if(m_verbose >=2) Cout("iDeformationRigid::ReadAndCheckOptionsList ..."<< endl); 
  
  Cerr("***** iDeformationRigid::ReadAndCheckOptionsList -> This deformation algorithm needs to use a file for its initialization. Use help for more info !" << endl);
  return 1;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn CheckSpecificParameters
  \brief This function is used to check parameters after the latter
         have been all set using Set functions.
  \return 0 if success, positive value otherwise.
*/
int iDeformationRigid::CheckSpecificParameters()
{
  if(m_verbose >=2) Cout("iDeformationRigid::CheckSpecificParameters ..."<< endl); 
  
  // Check vector files
  if (m_pathToFwdDeformationFiles.empty() || 
      m_pathToBwdDeformationFiles.empty() )
  {
    Cerr("***** iDeformationRigid::CheckSpecificParameters() -> Vector files not initialized !" << endl);
    return 1;
  }
  
  // todo consistency of param files here
  
  // Normal end
  m_checked = true;
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn Initialize
  \brief This function is used to initialize specific stuff in the deformation model.
  \details Allocate memory for image matrices, and initialize images containing transformation parameters
  \return 0 if success, other value otherwise.
*/
int iDeformationRigid::Initialize()
{
  if(m_verbose >=2) Cout("iDeformationRigid::Initialize ..."<< endl); 
  
  // Forbid initialization without check
  if (!m_checked)
  {
    Cerr("***** iDeformationRigid::Initialize() -> Parameters should be checked before Initialize() !" << endl);
    return 1;
  }
  
  // Memory allocation
  m_OriginalArray = new double [mp_ID->GetNbVoxXYZ()];
  m_OutputArray   = new double [mp_ID->GetNbVoxXYZ()];
  mp_vectorX = new FLTNB [mp_ID->GetNbVoxXYZ()];
  mp_vectorY = new FLTNB [mp_ID->GetNbVoxXYZ()];
  mp_vectorZ = new FLTNB [mp_ID->GetNbVoxXYZ()];
  
  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
  {
    mp_vectorX[v] = 0;
    mp_vectorY[v] = 0;
    mp_vectorZ[v] = 0;
  }

  // Normal end
  m_initialized = true;
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ApplyDeformations
  \param ap_inputImage : input image to deform
  \param ap_outputImage : image in which the output of the deformation should be recovered
  \param a_direction : a direction for the deformation to perform (forward or backward)
  \param a_defIdx : index of the deformation
  \brief This function prepares the deformation to perform 
  \details 1. Selects the right deformation parameters file according to the direction and deformation index argument
           2. Copy the image to deform in the buffer image of this class
           3. Call the deformation function (WarpImage)
           4. Copy back the output of the deformation image to the output image matrice passed in argument
  \return 0 if success, other value otherwise.
*/
int iDeformationRigid::ApplyDeformations(FLTNB* ap_inputImage, FLTNB* ap_outputImage, int a_direction, int a_defIdx)
{
  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** iDeformationRigid::ApplyDeformations() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif

  #ifdef CASTOR_VERBOSE
  if(m_verbose >=4) Cout("iDeformationRigid::ApplyDeformations ... "<< endl); 
  #endif
  
  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("*****iDeformationRigid::ApplyDeformationToForwardImage()->Error, unknown direction type( must be 'forward' or 'backward' # " << endl);
    return 1;
  }

  for(int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
  {
    m_OriginalArray[v] = (double) ap_inputImage[v];
    m_OutputArray[v] = 0.0;
  }
  
  // Apply all required deformations
  WarpImage(path_to_deformation_file, m_OriginalArray, m_OutputArray); 

  for(int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
  {
    if (m_OutputArray[v] < 0) m_OutputArray[v] = 0;
    ap_outputImage[v] = (FLTNB) m_OutputArray[v];
  }
  
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn WarpImage
  \param a_vectorFileIdx : deformation vector file
  \param ap_inputImage : image on which the deformation should be performed
  \param ap_outputImage : matrice which recovers the image after deformation
  \brief This function performs a rigid transformation using x,y,z displacement vectors and trilinear interpolation
  \details 1. Load the X,Y,Z transformation vectors from the file passed in argument
           2. Loop on the voxels and compute the new voxels value using trilinear interpolation
  \todo : perhaps add options to load all the vector parameters in RAM during initialization, instead on reading files during reconstruction
  \todo : perhaps use padded image in order to avoid if statements
  \return 0 if success, other value otherwise.
*/
int iDeformationRigid::WarpImage(string a_vectorFileIdx, double *ap_inputImage, double *ap_outputImage)
{
  #ifdef CASTOR_VERBOSE
  if(m_verbose >=4)
  {
    Cout("iDeformationRigid::WarpImage using the following coefficients file : "<< endl); 
    Cout(a_vectorFileIdx << endl); 
  }
  #endif
  
  // Get the vector coefficients
  // Todo : perhaps add options to load all the vector parameters in images matrices during initialization
  //        require more memory (x2*nb_transformations), but avoid data reading during reconstruction
  if(LoadVectorFiles(a_vectorFileIdx) )
  {
    Cerr("*****iDeformationRigid::WarpImage()-> An error occured while trying to read a deformation parameters file !" << endl);
    return 1;
  }
  
  for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
  {
    // For each axis component, get floor value, sign and fractional part of the vector
    int DepXEnt = (int) (mp_vectorX[v]) ;
    int DepXSign = (int) ( (mp_vectorX[v]==0) ? 0: mp_vectorX[v]/fabs(mp_vectorX[v]));
    float DX = fabs(mp_vectorX[v] - DepXEnt) ;
    
    int DepYEnt = (int) (mp_vectorY[v]) ;
    int DepYSign = (int) ((mp_vectorY[v]==0) ? 0: mp_vectorY[v]/fabs(mp_vectorY[v]));
    float DY = fabs(mp_vectorY[v] - DepYEnt) ;
    
    int DepZEnt = (int) (mp_vectorZ[v]) ;
    int DepZSign = (int) ((mp_vectorZ[v]==0) ? 0: mp_vectorZ[v]/fabs(mp_vectorZ[v]));
    float DZ = fabs(mp_vectorZ[v] - DepZEnt) ;
            
    // Get position of central voxel
    int q = v
            + DepXEnt 
            + DepYEnt * mp_ID->GetNbVoxX()
            + DepZEnt * mp_ID->GetNbVoxXY();
    
    // Trilinear interpolation
    // Todo : Use padded image in order to avoid 'if' statements ?
    if ((q <  mp_ID->GetNbVoxXYZ()) 
     && (q >= 0)) 
      ap_outputImage[v] = ap_inputImage[q] * (1-DX)*(1-DY)*(1-DZ);
             
    if ((q+DepXSign <  mp_ID->GetNbVoxXYZ()) 
     && (q+DepXSign >= 0))
        ap_outputImage[v] += ap_inputImage[q+DepXSign] * DX*(1-DY)*(1-DZ);
                
    if ((q+DepYSign*mp_ID->GetNbVoxX() <  mp_ID->GetNbVoxXYZ() )
    &&  (q+DepYSign*mp_ID->GetNbVoxX() >= 0))
        ap_outputImage[v] += ap_inputImage[q +DepYSign*mp_ID->GetNbVoxX()] * (1-DX)*DY*(1-DZ);
                                 
    if ((q+DepXSign+DepYSign*mp_ID->GetNbVoxX() <  mp_ID->GetNbVoxXYZ()) 
    && ( q+DepXSign+DepYSign*mp_ID->GetNbVoxX() >= 0))
        ap_outputImage[v] += ap_inputImage[q+DepXSign+DepYSign*mp_ID->GetNbVoxX()] * DX*DY*(1-DZ);
                
    if ((q+DepZSign*mp_ID->GetNbVoxXY()  <  mp_ID->GetNbVoxXYZ()) 
     && (q+DepZSign*mp_ID->GetNbVoxXY()) >= 0)
        ap_outputImage[v] += ap_inputImage[q+DepZSign*mp_ID->GetNbVoxXY()] * (1-DX)*(1-DY)*DZ;

    if ((q+DepXSign+DepZSign*mp_ID->GetNbVoxXY()) <  mp_ID->GetNbVoxXYZ()  
    &&  (q+DepXSign+DepZSign*mp_ID->GetNbVoxXY()) >=0 )
        ap_outputImage[v] += ap_inputImage[q+DepXSign+DepZSign*mp_ID->GetNbVoxXY()] * DX*(1-DY)*DZ;

    if ((q+DepYSign*mp_ID->GetNbVoxX()+DepZSign*mp_ID->GetNbVoxXY() <  mp_ID->GetNbVoxXYZ()) 
     && (q+DepYSign*mp_ID->GetNbVoxX()+DepZSign*mp_ID->GetNbVoxXY() >= 0))
        ap_outputImage[v] += ap_inputImage[q+DepYSign*mp_ID->GetNbVoxX()+DepZSign*mp_ID->GetNbVoxXY()] * (1-DX)*DY*DZ;

    if ((q+DepXSign+DepYSign*mp_ID->GetNbVoxX()+DepZSign*mp_ID->GetNbVoxXY() <  mp_ID->GetNbVoxXYZ()) 
    &&  (q+DepXSign+DepYSign*mp_ID->GetNbVoxX()+DepZSign*mp_ID->GetNbVoxXY() >= 0))
        ap_outputImage[v] += ap_inputImage[q+DepXSign+DepYSign*mp_ID->GetNbVoxX()+DepZSign*mp_ID->GetNbVoxXY()] *DX*DY*DZ;
  }  
  
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn LoadVectorFiles
  \param a_vectorFileIdx : deformation vector file
  \brief Load the transformation parameters from an ASCII file provided in argument
  \return 0 if success, other value otherwise.
*/
int iDeformationRigid::LoadVectorFiles(string a_vectorFileIdx)
{
  #ifdef CASTOR_VERBOSE
  if(m_verbose >=4) Cout("iDeformationRigid::LoadVectorFiles ... "<< endl); 
  #endif
  
  ifstream vector_file(a_vectorFileIdx.c_str(), ios::in);
  
  if(!vector_file) 
  {
    Cerr("*****iDeformationRigid::LoadVectorFiles()-> An error occured while trying to read the file " << a_vectorFileIdx << " !" << endl);
    return 1;
  }
  
  for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
  {
    vector_file >> mp_vectorX[v];
    vector_file >> mp_vectorY[v];
    vector_file >> mp_vectorZ[v];
  }

  return 0;
}