vDeformation.cc

Go to the documentation of this file.

/*
  Implementation of class vDeformation

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

#include "vDeformation.hh"
#include "oImageDimensionsAndQuantification.hh"
#include "oImageSpace.hh"


// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn vDeformation
  \brief Constructor of vDeformation. Simply set all data members to default values.
*/
vDeformation::vDeformation() 
{
  mp_ID = NULL;
  m_nbTransformations = -1;
  m_verbose = -1;
  m_checked = false;
  m_initialized = false;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ~vDeformation
  \brief Destructor of vDeformation.
*/
vDeformation::~vDeformation() {}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
int vDeformation::CheckParameters()
{
  if(m_verbose>=2) Cout("vDeformation::CheckParameters ..."<< endl); 
    
  // Check image dimensions
  if (mp_ID==NULL)
  {
    Cerr("***** vDeformation::CheckParameters() -> No image dimensions provided !" << endl);
    return 1;
  }
  
  // Check verbosity
  if (m_verbose<0)
  {
    Cerr("***** vDeformation::CheckParameters() -> Wrong verbosity level provided !" << endl);
    return 1;
  }

  // Check number of basis functions
  if (m_nbTransformations <0)
  {
    Cerr("***** vDeformation::CheckParameters() -> Number of transformations in the deformation has not been initialized !" << endl);
    return 1;
  }
  
  // Check parameters of the child class (if this function is overloaded)
  if (CheckSpecificParameters())
  {
    Cerr("***** vDeformation::CheckParameters() -> An error occurred while checking parameters of the child dynamic class !" << endl);
    return 1;
  }

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




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ApplyDeformationsToForwardImage
  \param ap_Image : required to access the forward image
  \param a_defIdx : index of the deformation
  \param fr : frame index
  \param rimg : respiratory image index
  \param cimg : cardiac image index
  \brief Apply initial deformations on the forwardImage before loop on events
  \todo May not be required anymore for the current implementation (always reference position at start of subset)
  \return 0 if success, positive value otherwise
*/
int vDeformation::ApplyDeformationsToForwardImage(oImageSpace* ap_Image, int a_defIdx, int fr, int rimg, int cimg)
{
  if (m_verbose >= 2) Cout("vDeformation::ApplyDeformationsToForwardImage ... " <<endl);
  
  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** vDeformation::ApplyDeformationsToForwardImage() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif
  
  if (ApplyDeformations(ap_Image->m4p_forwardImage[fr][rimg][cimg], ap_Image->m4p_forwardImage[fr][rimg][cimg], FORWARD_DEFORMATION, a_defIdx) )
  {
    Cerr("***** vDeformation::ApplyDeformationsToForwardImage()-> An error occured while performing initial deformation of the forward image !" << endl);
    return 1;
  }

  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ApplyDeformationsToBackwardImage
  \param ap_Image : required to access the backward image and its deformation backup matrice
  \param a_defIdx : index of the deformation
  \brief Apply backward transformation of the backward image to the reference position
  \details Loop on frames
           Recover any potential data stored in the backup matrice m2p_defTmpBackwardImage
  \return 0 if success, positive value otherwise
*/
int vDeformation::ApplyDeformationsToBackwardImage(oImageSpace* ap_Image, int a_defIdx)
{
  if(m_verbose >= 2) Cout("vDeformation::ApplyDeformationsToBackwardImage ... " <<endl);

  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** vDeformation::ApplyDeformationsToBackwardImage() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif
  
  for(int bimg=0; bimg<ap_Image->GetNbBackwardImages(); bimg++)
    for(int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++)
      for(int rimg=0; rimg<mp_ID->GetNbRespGates(); rimg++)
        for(int cimg=0; cimg<mp_ID->GetNbCardGates(); cimg++)
        {
          // Perform backward deformation
          if(ApplyDeformations(ap_Image->m6p_backwardImage[bimg][0][fr][rimg][cimg], ap_Image->m6p_backwardImage[bimg][0][fr][rimg][cimg], BACKWARD_DEFORMATION, a_defIdx) )
          {
            Cerr("***** vDeformation::ApplyDeformationsToBackwardImage()-> An error occured while performing backward deformation of the backward image !" << endl);
          Cerr("*****                                                      frame index " << fr << " respiratory image index " << rimg<< " cardiac image index " << cimg<<  " !" << endl);
            return 1;
          }
          
          // Recover the content of the temporary backup deformation image to the backward image
          for(int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
            ap_Image->m6p_backwardImage[bimg][0][fr][rimg][cimg][v] += ap_Image->m5p_defTmpBackwardImage[bimg][fr][rimg][cimg][v] ;
        }
  
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ApplyDeformationsToHistoSensitivityImage
  \param ap_Image : required to access the backward image and its deformation backup matrice
  \param a_defIdx : index of the deformation
  \brief Apply backward transformations of the sensitivity image to the reference position (histogram mode)
  \details Loop on frames
          Recover any potential data stored in the backup matrice m4p_defTmpSensitivityImage
  \return 0 if success, positive value otherwise
*/
int vDeformation::ApplyDeformationsToHistoSensitivityImage(oImageSpace* ap_Image, int a_defIdx)
{
  if(m_verbose >= 2) Cout("vDeformation::ApplyDeformationsToHistoSensitivityImage ... " <<endl);
  
  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** vDeformation::ApplyDeformationsToHistoSensitivityImage() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif
  
  for(int fr=0; fr<mp_ID->GetNbTimeFrames(); fr++)
    for(int rimg=0; rimg<mp_ID->GetNbRespGates() ; rimg++)
      for(int cimg=0; cimg<mp_ID->GetNbCardGates() ; cimg++)
      {
        // Perform backward deformation
        if( ApplyDeformations(ap_Image->m5p_sensitivity[0][fr][rimg][cimg], ap_Image->m5p_sensitivity[0][fr][rimg][cimg], BACKWARD_DEFORMATION, a_defIdx) )
        {
          Cerr("***** vDeformation::ApplyDeformationsToHistoSensitivityImage()-> An error occured while performing backward deformation of the backward image !" << endl);
          Cerr("*****                                                            frame index " << fr << " respiratory image index " << rimg<< " cardiac image index " << cimg<<  " !" << endl);
          return 1;
        }
        
        // Recover the content of the temporary backup deformation image to the sensitivity image
        for(int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
          ap_Image->m5p_sensitivity[0][fr][rimg][cimg][v] += ap_Image->m4p_defTmpSensitivityImage[fr][rimg][cimg][v];
      }

  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PerformDeformation
  \param ap_Image : required to access oImageSpace image matrices
  \param a_defIdx : index of the deformation
  \param fr : frame index
  \param rimg : respiratory image index
  \param cimg : cardiac image index
  \brief Apply deformations during reconstruction
  \details 1. Recover all the data of the multithreaded backward image matrices in the first one (thread index 0)
           2. Perform backward deformation of the backward image to the reference position with defIdx-1
           3. Add coefficients of the backward image matrice to the temporary backup image matrice & reset backward image
           4. Apply forward deformation of the forward image (backward deformation to the reference position with defIdx-1, and forward deformation with defIdx)
  \return 0 if success, positive value otherwise
*/
int vDeformation::PerformDeformation(oImageSpace* ap_Image, int a_defIdx, int fr, int rimg, int cimg)
{
  if(m_verbose >= 2) Cout("vDeformation::PerformDeformation ... " <<endl);
  
  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** vDeformation::PerformDeformation() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif
  
  // REDUCE
  for (int bimg=0; bimg<ap_Image->GetNbBackwardImages(); bimg++)
  {
    for (int th=1 ; th<mp_ID->GetNbThreadsForProjection() ; th++) //TODO add loops
    {
      // Synchronisation of the multi-threaded backwardImages inside the first image Reduce)
      for(int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
        ap_Image->m6p_backwardImage[bimg][0][fr][rimg][cimg][v] += ap_Image->m6p_backwardImage[bimg][th][fr][rimg][cimg][v];
    }

    // BACKWARD DEFORMATION of the backward image  //TODO add loops
    if (ApplyDeformations(ap_Image->m6p_backwardImage[bimg][0][fr][rimg][cimg], ap_Image->m6p_backwardImage[bimg][0][fr][rimg][cimg], BACKWARD_DEFORMATION, a_defIdx-1) )
    {
      Cerr("***** vDeformation::ApplyDeformations()-> An error occured while performing backward deformation of the backward image !" << endl);
      Cerr("*****                                     frame index " << fr << " respiratory image index " << rimg<< " cardiac image index " << cimg<<  " !" << endl);
      return 1;
    }

    // Store Backward deformation update coefficients in the temporary backup image  //TODO add loops
    for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
    {
      ap_Image->m5p_defTmpBackwardImage[bimg][fr][rimg][cimg][v] += ap_Image->m6p_backwardImage[bimg][0][fr][rimg][cimg][v];
    }
  }
  
  // Reset the backward images (i.e set all voxels to the specific fr/rimg/cimg to 0)
  for (int bimg=0; bimg<ap_Image->GetNbBackwardImages(); bimg++)
    for(int th=0 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
      for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
        ap_Image->m6p_backwardImage[bimg][th][fr][rimg][cimg][v] = 0.;
        

  // BACKWARD DEFORMATION of the Forward image //TODO add loops
  if (ApplyDeformations(ap_Image->m4p_forwardImage[fr][rimg][cimg], ap_Image->m4p_forwardImage[fr][rimg][cimg], BACKWARD_DEFORMATION, a_defIdx-1) )
  {
    Cerr("***** vDeformation::ApplyDeformations()-> An error occured while performing backward deformation of the forward image !" << endl);
    Cerr("*****                                     frame index " << fr << " respiratory image index " << rimg<< " cardiac image index " << cimg<<  " !" << endl);
    return 1;
  }
  // FORWARD DEFORMATION of the Forward image
  if (ApplyDeformations(ap_Image->m4p_forwardImage[fr][rimg][cimg], ap_Image->m4p_forwardImage[fr][rimg][cimg], FORWARD_DEFORMATION, a_defIdx) )
  {
    Cerr("***** vDeformation::ApplyDeformations()-> An error occured while performing forward deformation of the forxar image !" << endl);
    Cerr("*****                                     frame index " << fr << " respiratory image index " << rimg<< " cardiac image index " << cimg<<  " !" << endl);
    return 1;
  }
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PerformHistoSensitivityDeformation
  \param ap_Image : required to access oImageSpace image matrices
  \param a_defIdx : index of the deformation
  \param fr : frame index
  \param rimg : respiratory image index
  \param cimg : cardiac image index
  \brief Apply deformations on the sensitivity image during reconstruction in histogram mode
  \details 1. Recover all the data of the multithreaded sensitivity image matrice in the first one (thread index 0)
           2. Perform backward deformation of the sensitivity image to the reference position with defIdx-1
           3. Add coefficients of the sensitivity image matrice to the temporary backup image matrice & reset sensitivity image
  \return 0 if success, positive value otherwise
*/
int vDeformation::PerformHistoSensitivityDeformation(oImageSpace* ap_Image, int a_defIdx, int fr, int rimg, int cimg)
{
  if(m_verbose >= 2) Cout("vDeformation::PerformHistoSensitivityDeformation ... " <<endl);
  
  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** vDeformation::PerformHistoSensitivityDeformation() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif
  
  // REDUCE
  for (int th=1 ; th<mp_ID->GetNbThreadsForProjection() ; th++)
  {
    // Synchronisation of the multi-threaded sensitivity images inside the first image )
    for(int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
      ap_Image->m5p_sensitivity[0][fr][rimg][cimg][v] += ap_Image->m5p_sensitivity[th][fr][rimg][cimg][v];
  }
  
  // BACKWARD DEFORMATION of the sensitivity image
  if (ApplyDeformations(ap_Image->m5p_sensitivity[0][fr][rimg][cimg], ap_Image->m5p_sensitivity[0][fr][rimg][cimg], BACKWARD_DEFORMATION, a_defIdx-1) )
  {
    Cerr("***** vDeformation::PerformHistoSensitivityDeformation()-> An error occured while performing backward deformation of the sensitivity image !" << endl);
    Cerr("*****                                                      frame index " << fr << " respiratory image index " << rimg<< " cardiac image index " << cimg<<  " !" << endl);
    return 1;
  }
  
  // Store Backward deformation update coefficients in temporary image
  for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
  {
    ap_Image->m4p_defTmpSensitivityImage[fr][rimg][cimg][v] += ap_Image->m5p_sensitivity[0][fr][rimg][cimg][v];
  }

  // Reset the sensitivity images (i.e set all voxels to the specific fr/rimg/cimg to 0) 
  for (int th=0; th<mp_ID->GetNbThreadsForProjection(); th++)
    for (int v=0; v<mp_ID->GetNbVoxXYZ(); v++)
      ap_Image->m5p_sensitivity[th][fr][rimg][cimg][v] = 0.;

  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PerformSensitivityDeformation
  \param ap_Image : required to access oImageSpace image matrices
  \param a_defDirection : a direction for the deformation to perform (forward or backward)
  \param a_defIdx : index of the deformation
  \param fr : frame index
  \param rg : respiratory gate index
  \param cg : cardiac gate index
  \brief Apply image deformations during sensitivity image generation for list-mode
  \details Depending on the deformation direction (forward or backward):
           Forward : Perform forward deformation of the forward image to the deformation index position
           Backward: Perform backward deformation of the backward image to the reference position
  \return 0 if success, positive value otherwise
*/
int vDeformation::PerformSensitivityDeformation(oImageSpace* ap_Image, int a_defDirection, int a_defIdx, int fr, int rg, int cg)
{
  if(m_verbose >= 2) Cout("vDeformation::PerformSensitivityDeformation ... " <<endl);
  
  #ifdef CASTOR_DEBUG
  if (!m_initialized)
  {
    Cerr("***** vDeformation::PerformSensitivityDeformation() -> Called while not initialized !" << endl);
    Exit(EXIT_DEBUG);
  }
  #endif
  
  if (a_defDirection == FORWARD_DEFORMATION)
  {
    if (ApplyDeformations(ap_Image->m4p_forwardImage[fr][rg][cg], ap_Image->m4p_forwardImage[fr][rg][cg], a_defDirection, a_defIdx) )
    {
      Cerr("***** vDeformation::PerformSensitivityDeformation -> An error occured while performing forward deformation !" << endl);
      Cerr("*****                                                frame index " << fr << " respiratory gate index " << rg<< " cardiac gate index " << cg<<  " !" << endl);
      return 1;
    }
  }
  else if (a_defDirection == BACKWARD_DEFORMATION)
  {
    if (ApplyDeformations(ap_Image->m6p_backwardImage[0][0][fr][rg][cg], ap_Image->m6p_backwardImage[0][0][fr][rg][cg], a_defDirection, a_defIdx) )
    {
      Cerr("***** vDeformation::PerformSensitivityDeformation -> An error occured while performing backward deformation !" << endl);
      Cerr("*****                                                frame index " << fr << " respiratory gate index " << rg<< " cardiac gate index " << cg<<  " !" << endl);
      return 1;
    }
  }
  else 
  {
    Cerr("***** vDeformation::PerformDeformation -> Unknown type of deformation !" << endl);
    return 1;
  }
  
  return 0;
}