iScannerSPECTConv.cc

Go to the documentation of this file.

/*
  Implementation of class iScannerSPECTConv

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

#include "iScannerSPECTConv.hh"
#include "sOutputManager.hh"
#include "sScannerManager.hh"





// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \brief iScannerSPECTConv constructor. 
         Initialize the member variables to their default values.
*/
iScannerSPECTConv::iScannerSPECTConv() : vScanner() 
{
  // Set member variables to default values
  m_scannerType = 2;
  m_nbCrystals = -1;
  m_nbHeads = -1;
  m_nbPixelsTrans = 0;
  m_pixelsSizeTrans = -1.;
  m_gapSizeTrans = -1.;

  m_nbPixelsAxial = 0;
  m_pixelsSizeAxial = -1.;
  m_gapSizeAxial = -1.;

  m_vPixelsSizeTrans = -1.;
  m_vPixelsSizeAxial = -1.;

  m_vNbPixelsTrans = 0;
  m_vNbPixelsAxial = 0;

  mp_nbOfBins[0] = 0;
  mp_nbOfBins[1] = 0;

  m_crystalDepth = -1.;

  mp_focalModelTrans = NULL;
  mp_nbCoefModelTrans = NULL;
  m2p_transFocalParameters = NULL;    

  mp_focalModelAxial = NULL;
  mp_nbCoefModelAxial = NULL;
  m2p_axialFocalParameters = NULL;  

  m_nbOfProjections = 0;
  mp_projectionAngles = NULL;
  mp_radius = NULL;
  mp_CORtoDetectorDistance = NULL; 
  m_rotDirection = GEO_ROT_CW;

  mp_crystalCentralPositionX = NULL; 
  mp_crystalCentralPositionY = NULL; 
  mp_crystalCentralPositionZ = NULL; 

  mp_crystalOrientationX = NULL; 
  mp_crystalOrientationY = NULL; 
  mp_crystalOrientationZ = NULL; 
 
  mp_crystalFocalPositionX = NULL; 
  mp_crystalFocalPositionY = NULL; 
  mp_crystalFocalPositionZ = NULL; 

  mp_positionMatrix_ref = NULL;
  mp_positionMatrix_out = NULL;
  mp_rotationMatrix = NULL;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \brief iScannerSPECTConv destructor. 
*/
iScannerSPECTConv::~iScannerSPECTConv() 
{
  if( mp_focalModelTrans ) delete[] mp_focalModelTrans;
  if( mp_nbCoefModelTrans ) delete[] mp_nbCoefModelTrans;

  if( m2p_transFocalParameters )
  {
    for( int i = 0; i < m_nbHeads; ++i )
      delete m2p_transFocalParameters[ i ];
    
    delete[] m2p_transFocalParameters;
  }
    

  if( mp_focalModelAxial ) delete[] mp_focalModelAxial;
  if( mp_nbCoefModelAxial ) delete[] mp_nbCoefModelAxial;

  if( m2p_axialFocalParameters )
  {
    for( int i = 0; i < m_nbHeads; ++i )
      delete m2p_axialFocalParameters[ i ];
    
    delete[] m2p_axialFocalParameters;
  }
  
  if( mp_projectionAngles ) delete[] mp_projectionAngles;
  if( mp_radius ) delete[] mp_radius;
  if( mp_CORtoDetectorDistance ) delete[] mp_CORtoDetectorDistance; 

  if (mp_positionMatrix_ref) delete mp_positionMatrix_ref;
  if (mp_positionMatrix_out) delete mp_positionMatrix_out;
  if (mp_rotationMatrix) delete mp_rotationMatrix;


  if(mp_crystalCentralPositionX) delete mp_crystalCentralPositionX;
  if(mp_crystalCentralPositionY) delete mp_crystalCentralPositionY;
  if(mp_crystalCentralPositionZ) delete mp_crystalCentralPositionZ;

  if(mp_crystalOrientationX) delete mp_crystalOrientationX;
  if(mp_crystalOrientationY) delete mp_crystalOrientationY;
  if(mp_crystalOrientationZ) delete mp_crystalOrientationZ;
 
  if(mp_crystalFocalPositionX) delete mp_crystalFocalPositionX;
  if(mp_crystalFocalPositionY) delete mp_crystalFocalPositionY;
  if(mp_crystalFocalPositionZ) delete mp_crystalFocalPositionZ;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iScannerSPECTConv::Instantiate(bool a_scannerFileIsLUT)
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::Instantiate ..."<< endl); 
    
  sScannerManager* p_scannerManager; 
  p_scannerManager = sScannerManager::GetInstance();  
  
  // Get the number of layers in the scanner
  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "number of detector heads", &m_nbHeads, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read the number of SPECT heads !" << endl);
    return 1;
  }
  
  mp_radius = new FLTNB[m_nbHeads];
  
  // Get default radius
  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "scanner radius", mp_radius, m_nbHeads, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read the number of SPECT heads !" << endl);
    return 1;
  }

  for( int i = 0; i < m_nbHeads; i++ )
    if(mp_radius[i] == 0)
    {
      Cerr("***** iScannerSPECTConv::Instantiate() -> scanner radius == 0 ? really ? :) ... " << endl);
      return 1;
    } 

  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "trans number of pixels", &m_nbPixelsTrans, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read the transaxial number of pixels !" << endl);
    return 1;
  }
 
  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "trans pixel size", &m_pixelsSizeTrans, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read the transaxial pixel size !" << endl);
    return 1;
  }
  
  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "axial number of pixels", &m_nbPixelsAxial, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read the axial number of pixels !" << endl);
    return 1;
  } 
  
  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "axial pixel size", &m_pixelsSizeAxial, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read the axial pixel size !" << endl);
    return 1;
  } 

  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "detector depth", &m_crystalDepth, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read detector depth !" << endl);
    return 1;
  } 
  

  string *head_name = new string[ m_nbHeads + 1 ];
  for( int i = 0; i < m_nbHeads; i++ )
  {
    ostringstream oss( ostringstream::out );
    oss << "head" << i+1;
    head_name[ i ] = oss.str();
  }
  head_name[ m_nbHeads ] = "eof";


  if(m_verbose>=2) Cout("iScannerSPECTConv::Allocating memory ..."<< endl); 
  
  // Allocating memories
  mp_focalModelTrans = new string[ m_nbHeads ];
  mp_focalModelAxial = new string[ m_nbHeads ];
  mp_nbCoefModelTrans = new uint8_t[ m_nbHeads ];
  mp_nbCoefModelAxial = new uint8_t[ m_nbHeads ];
  m2p_transFocalParameters = new FLTNB*[ m_nbHeads ];
  m2p_axialFocalParameters = new FLTNB*[ m_nbHeads ];
  
  
  
  // Read informations about SPECT head
  for( int i = 0; i < m_nbHeads; i++ )
  {
    if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "trans focal model", &mp_focalModelTrans[ i ], 1, KEYWORD_MANDATORY, head_name[ i ], head_name[ i + 1 ]))
    {
      Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read transaxial focal model of head " << i << " !" << endl);
      return 1;
    }

    if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "trans number of coef model", &mp_nbCoefModelTrans[ i ], 1, KEYWORD_MANDATORY, head_name[ i ], head_name[ i + 1 ]))
    {
      Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read transaxial number of coef model of head " << i << " !" << endl);
      return 1;
    }

    if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "axial focal model", &mp_focalModelAxial[ i ], 1, KEYWORD_MANDATORY, head_name[ i ], head_name[ i + 1 ]))
    {
      Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read axial focal model of head " << i << " !" << endl);
      return 1;
    }
    
    if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "axial number of coef model", &mp_nbCoefModelAxial[ i ], 1, KEYWORD_MANDATORY, head_name[ i ], head_name[ i + 1 ]))
    {
      Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read transaxial focal model of head " << i << " !" << endl);
      return 1;
    }
    
    m2p_transFocalParameters[i] = new FLTNB[mp_nbCoefModelTrans[ i ]];
    m2p_axialFocalParameters[i] = new FLTNB[mp_nbCoefModelAxial[ i ]]; 
    
    if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "trans parameters", m2p_transFocalParameters[i], mp_nbCoefModelTrans[ i ], KEYWORD_MANDATORY, head_name[ i ], head_name[ i + 1 ]))
    {
      Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read transaxial focal model of head " << i << " !" << endl);
      return 1;
    }
    
    if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "axial parameters", m2p_axialFocalParameters[i], mp_nbCoefModelAxial[ i ], KEYWORD_MANDATORY, head_name[ i ], head_name[ i + 1 ]))
    {
      Cerr("***** iScannerSPECTConv::Instantiate() -> An error occurred while trying to read transaxial focal model of head " << i << " !" << endl);
      return 1;
    }
  }

  // Instanciation of objects for matrix calculation
  mp_positionMatrix_ref = new oMatrix(3,1);
  mp_positionMatrix_out = new oMatrix(3,1);
  mp_rotationMatrix = new oMatrix(3,3);  
  
  delete[] head_name;
   
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn BuildLUT
  \param a_scannerFileIsLUT : boolean indicating if the file describing 
                              the SPECT camera is a generic file (0) or custom Look-up-table (1)
  \brief Call the functions to generate the LUT or read the user-made LUT depending on the user choice
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::BuildLUT(bool a_scannerFileIsLUT)
{  
  if(m_verbose>=2) Cout("iScannerSPECTConv::BuildLUT ..."<< endl);
  
  // Initialize nb_crystals
  if(mp_nbOfBins[0] == 0 && mp_nbOfBins[1] == 0) // number of bins not initialized in the datafile, then we read the number of crystals from the geom file
  {
    if(m_nbPixelsTrans==0 || m_nbPixelsAxial==0)
    {
     Cerr("***** iScannerSPECTConv::BuildLUT() -> Error, transaxial or axial number of pixels in the geometric file should be >0 !" << endl);
     return 1;
    }
    
    m_nbCrystals = m_nbPixelsTrans*m_nbPixelsAxial;
  }
  else // initialization with number of bins
  {
    m_nbCrystals = mp_nbOfBins[0]*mp_nbOfBins[1];
  }
  
  mp_crystalCentralPositionX = new FLTNB[m_nbOfProjections*m_nbCrystals];
  mp_crystalCentralPositionY = new FLTNB[m_nbOfProjections*m_nbCrystals];
  mp_crystalCentralPositionZ = new FLTNB[m_nbOfProjections*m_nbCrystals];

  mp_crystalOrientationX = new FLTNB[m_nbOfProjections*m_nbCrystals];
  mp_crystalOrientationY = new FLTNB[m_nbOfProjections*m_nbCrystals];
  mp_crystalOrientationZ = new FLTNB[m_nbOfProjections*m_nbCrystals];
  
  mp_crystalFocalPositionX = new FLTNB[m_nbOfProjections*m_nbCrystals];
  mp_crystalFocalPositionY = new FLTNB[m_nbOfProjections*m_nbCrystals];
  mp_crystalFocalPositionZ = new FLTNB[m_nbOfProjections*m_nbCrystals]; 
    
  // Either generate the LUT from a generic file, or load the user precomputed LUT
  if (!a_scannerFileIsLUT)
  {
    if (ComputeLUT())
    {
     Cerr("***** iScannerSPECTConv::BuildLUT() -> A problem occured while generating scanner LUT !" << endl);
     return 1;
    }
  }
  else 
  {
    if (LoadLUT())
    {
      Cerr("***** iScannerSPECTConv::BuildLUT() -> A problem occured while loading scanner LUT !" << endl);
      return 1;
    }
  }
  
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn CheckParameters
  \brief Check that all parameters have been correctly initialized.
  \return 0 if success, positive value otherwise
  \todo Keep the check on crystal depth ?
*/
int iScannerSPECTConv::CheckParameters()
{  
  if(m_verbose>=2) Cout("iScannerSPECTConv::CheckParameters ..."<< endl);
  
  // Check if all parameters have been correctly initialized. Return error otherwise

  if(m_nbCrystals<0)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, number of crystals has not been initialized !" <<endl);
    return 1;
  }

  if(m_nbHeads<0)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, number of heads in the SPECT system has not been initialized !" <<endl);
    return 1;
  }

  if(m_nbOfProjections<=0)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, number of projection angles has not been initialized !" <<endl);
    return 1;
  }
  
  if(m_nbPixelsTrans<=0 || m_nbPixelsAxial<=0 )
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, number of transaxial/axial pixels have not correctly been initialized ! (should be >0)" <<endl);
    return 1;
  }
    
  if(m_pixelsSizeTrans<=0 || m_pixelsSizeAxial<=0 )
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, transaxial/axial pixel sizes have not correctly been initialized ! (should be >0)" <<endl);
    return 1;
  }


  if(m_vNbPixelsTrans<=0 || m_vNbPixelsAxial<=0 )
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, transaxial/axial number of virtual pixels have not correctly been initialized ! (should be >0)" <<endl);
    return 1;
  }


  if(m_vPixelsSizeTrans<=0 || m_vPixelsSizeAxial<=0 )
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, transaxial/axial virtual pixel sizes have not correctly been initialized ! (should be >0)" <<endl);
    return 1;
  }


  if(m_gapSizeTrans<0 || m_gapSizeAxial<0 )
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, transaxial/axial pixel gap sizes have not correctly been initialized ! (should be >0)" <<endl);
    return 1;
  }

  // todo : maybe delete this as we don't used it
  if(m_crystalDepth<=0)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, crystal depth has not correctly been initialized ! (should be >0)" <<endl);
    return 1;
  }

  if(mp_focalModelTrans == NULL || mp_focalModelAxial == NULL)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, transaxial/axial focal models have not correctly been initialized !" <<endl);
    return 1;
  }

  if(mp_nbCoefModelTrans == NULL || mp_nbCoefModelAxial == NULL)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, number of coefficients for the transaxial/axial focal models have not correctly been initialized !" <<endl);
    return 1;
  }
  
  if(m2p_transFocalParameters == NULL || m2p_axialFocalParameters == NULL)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, parameters for the transaxial/axial focal models have not correctly been initialized !" <<endl);
    return 1;
  }
  
  if(mp_projectionAngles == NULL)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, projection angles have not correctly been initialized !" <<endl);
    return 1;
  }

  if(mp_radius == NULL)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, default scanner radius for each detector heads (scanner file) not correctly initialized !" <<endl);
    return 1;
  }
  
  if(mp_CORtoDetectorDistance == NULL)
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> Error, distances between center of rotation and each detector heads have not correctly been initialized !" <<endl);
    return 1;
  }
  
  if(mp_crystalCentralPositionX == NULL ||
     mp_crystalCentralPositionY == NULL ||
     mp_crystalCentralPositionZ == NULL  )
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> LUT elements (crystal central positions) have not correctly been initialized !" <<endl);
    return 1;
  }
    
  if(mp_crystalOrientationX == NULL ||
     mp_crystalOrientationY == NULL ||
     mp_crystalOrientationZ == NULL  )
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> LUT elements (crystal orientations) have not correctly been initialized !" <<endl);
    return 1;
  }

  if(mp_crystalFocalPositionX == NULL ||
     mp_crystalFocalPositionY == NULL ||
     mp_crystalFocalPositionZ == NULL  )
  {
    Cerr("*****iScannerSPECTConv::CheckParameters()-> LUT elements (crystal focal positions) have not correctly been initialized !" <<endl);
    return 1;
  }

  // No need to check mp_nbOfBins, Default values = 0

  m_allParametersChecked = true;
  
  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn Initialize
  \brief Check general initialization and set several parameters to their default value
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::Initialize()
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::Initialize ..."<< endl); 
  
  if(m_allParametersChecked == false)
  {
    Cerr("***** iScannerSPECTConv::Initialize() -> Parameters have not been checked !" << endl);
    return 1;
  }
  
  // Any initialization
  
  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn LoadLUT
  \brief Load a precomputed scanner LUT
  \details Read mandatory data from the header of the LUT. Then load the LUT elements for each crystal
  \todo Not yet implemented
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::LoadLUT()
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::LoadLUT ..."<< endl); 
  
  Cerr("iScannerSPECTConv::LoadLUT() -> Not yet implemented !" << endl);
  return 1;
  // Read LUT according to SPECT camera hscan/LUT files. Check the Castor presentation file to see mandatory/optional fields to read
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ComputeLUT
  \brief Computes the LUT of the scanner from a generic (.geom) file.
  \details Read mandatory data from the geom file. Then compute the LUT elements for each crystal from the geometry described in the file
           Compute the look-up-tables of the system containing the locations of the scanner elements center in space and their orientations
  \todo  center of rotation & head first angles : get this from acquisition header file
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::ComputeLUT()
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::ComputeLUT() -> Start LUT generation" << endl);
  
  // Compute LUT according to SPECT camera geom files. Check the Castor presentation file for more details about the related  mandatory/optional fields to be read

  // ============================================================================================================
  // Parameters declarations
  // ============================================================================================================
  
  // todo center of rotation & head first angles : get this from acquisition header file
  FLTNB *CORx, *CORy, *CORz;
  FLTNB *head_angleX, *head_angleY, *head_angleZ;
  
  CORx = new FLTNB[m_nbHeads];
  CORy = new FLTNB[m_nbHeads]; 
  CORz = new FLTNB[m_nbHeads];
  head_angleX = new FLTNB[m_nbHeads];
  head_angleY = new FLTNB[m_nbHeads];
  head_angleZ = new FLTNB[m_nbHeads];
    
  //  for(int hId=0 ; hId<m_nbHeads ; hId++)
    // todo Initialisation COR and head_angle

  // ============================================================================================================
  // Generate the LUT
  // ============================================================================================================
  
  // oMatrix crystal_center_ref will be used to gather cartesian positions of the crystals center (on the surface of the crystals) 
  // in the reference head (directly above isocenter)
  oMatrix** crystal_center_ref = new oMatrix *[m_nbCrystals];
  // oMatrix crystal_center_out will be used to recover positions of each crystal after each rotation (each projection angles)
  oMatrix* crystal_center_out = new oMatrix(3,1);
  
  // same method for matrices recovering the positions of focal point
  oMatrix* focal_projection_position_mtx = new oMatrix(3,3);
  oMatrix* focal_projection_position_mtx_output = new oMatrix(3,3);
  
  for (int c = 0; c<m_nbCrystals; c++)
    crystal_center_ref[c]  = new oMatrix(3,1);


  for(int i=0 ; i<3 ; i++)
  {
    crystal_center_out->SetMatriceElt(i,0,0);
    
    for(int j=0 ; j<3 ; j++)
    {
      focal_projection_position_mtx->SetMatriceElt(i,j,0);
      focal_projection_position_mtx_output->SetMatriceElt(i,j,0);
    }
  }

  // Generation of the rotation matrix allowing to compute the position of all the rsectors. 
  oMatrix** rotation_mtx = new oMatrix*[m_nbOfProjections];

  for(int i=0; i<m_nbOfProjections; i++)
  {
    rotation_mtx[i] = new oMatrix(3,3);
    
    int dir = (m_rotDirection == GEO_ROT_CCW) ? -1 : 1;
    
    rotation_mtx[i]->SetMatriceElt(0,0, cos(mp_projectionAngles[i] * M_PI/180.) );
    rotation_mtx[i]->SetMatriceElt(1,0, -dir*sin(mp_projectionAngles[i] * M_PI/180.) );
    rotation_mtx[i]->SetMatriceElt(2,0,0);
    rotation_mtx[i]->SetMatriceElt(0,1, dir*sin(mp_projectionAngles[i] * M_PI/180.) );
    rotation_mtx[i]->SetMatriceElt(1,1, cos(mp_projectionAngles[i] * M_PI/180.) );
    rotation_mtx[i]->SetMatriceElt(2,1,0);
    rotation_mtx[i]->SetMatriceElt(0,2,0);
    rotation_mtx[i]->SetMatriceElt(1,2,0);
    rotation_mtx[i]->SetMatriceElt(2,2,1);
  }

  // Initialization of pixel size
  FLTNB size_pixel_trans = 0., size_pixel_axial = 0.;
  int   nb_pixel_trans=0, nb_pixel_axial=0;

  // nb of bins have been recovered from the datafile (>0)
  if(mp_nbOfBins[0]>0 && mp_nbOfBins[1]>0)
  {
    size_pixel_trans = m_pixelsSizeTrans/mp_nbOfBins[0];
    size_pixel_axial = m_pixelsSizeAxial/mp_nbOfBins[1];
    
    // Initialize pixel value
    nb_pixel_trans = mp_nbOfBins[0];
    nb_pixel_axial = mp_nbOfBins[1];
    
    // Set the member variable recovering the pixel value
    m_vNbPixelsTrans = mp_nbOfBins[0];
    m_vNbPixelsAxial = mp_nbOfBins[1];
  }
  else // Otherwise, we recover these informations from the geom file
  {
    size_pixel_trans = m_pixelsSizeTrans;
    size_pixel_axial = m_pixelsSizeAxial;
    nb_pixel_trans = m_nbPixelsTrans;
    nb_pixel_axial = m_nbPixelsAxial;

    // Set the member variable recovering the pixel value
    m_vNbPixelsTrans = m_nbPixelsTrans;
    m_vNbPixelsAxial = m_nbPixelsAxial;
  }

  // Recover the actual trans/axial size of the pixels in member variable, in case we need them during reconstruction
  m_vPixelsSizeTrans = size_pixel_trans;
  m_vPixelsSizeAxial = size_pixel_axial;
  

  // Recover the trans/axial gap size from the geom file
  sScannerManager* p_scannerManager; 
  p_scannerManager = sScannerManager::GetInstance();

  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "trans gap size", &m_gapSizeTrans, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::ComputeLUT() -> An error occurred while trying to read the transaxial gap size !" << endl);
    return 1;
  } 
  
  if(ReadDataASCIIFile(p_scannerManager->GetPathToScannerFile(), "axial gap size", &m_gapSizeAxial, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::ComputeLUT() -> An error occurred while trying to read the axial gap size !" << endl);
    return 1;
  } 
  
  
  // Generate cartesian coordinates of the crystal centers for the SPECT at the reference position (directly above isocenter)
  // Starting position of the SPECT camera for all axis
  FLTNB x_start = -(nb_pixel_trans*size_pixel_trans + (nb_pixel_trans-1)*m_gapSizeTrans) / 2 + (size_pixel_trans/2);
  FLTNB z_start = -(nb_pixel_axial*size_pixel_axial + (nb_pixel_axial-1)*m_gapSizeAxial) / 2 + (size_pixel_axial/2);

  for( int jj = 0; jj < nb_pixel_axial; ++jj )
  {
    FLTNB Zcrist = z_start + jj * (size_pixel_axial + m_gapSizeAxial);
    for( int ii = 0; ii < nb_pixel_trans; ++ ii )
    {
      FLTNB Xcrist = x_start + ii * (size_pixel_trans + m_gapSizeTrans);

      crystal_center_ref[ii + nb_pixel_trans * jj]->SetMatriceElt(0,0,Xcrist);
      crystal_center_ref[ii + nb_pixel_trans * jj]->SetMatriceElt(1,0,     0); // This value will be set for each projection angles
      crystal_center_ref[ii + nb_pixel_trans * jj]->SetMatriceElt(2,0,Zcrist);
    }
  }

  // ============================================================================================================
  // Loop over all the projection angles
  // Positions of the scanner elements are progressively stored in the LUT file.
  // ============================================================================================================
  if (m_verbose>=2) Cout("  --> Generate positions for each projection angle"<< endl); 
  
  for (int a=0 ; a<m_nbOfProjections ; a++)
  {
    // Set y-position of the crystal reference matrix according to the distance between SPECT camera and center of rotation of the specific head.(mp_CORtoDetectorDistance)
    int hId = a*m_nbHeads/m_nbOfProjections;

    // Get surface crystal position (and not center crystal positions for SPECT), and set the reference crystal matrix
    FLTNB Ycrist = mp_CORtoDetectorDistance[hId];

    for (int c=0 ; c < m_nbCrystals ; c++) 
      crystal_center_ref[c]->SetMatriceElt(1,0,Ycrist);

    for (int c=0 ; c<m_nbCrystals ; c++)
    {
      // crystal indexation
      int cryID = a*m_nbCrystals + c;

      // SET CARTESIAN COORDINATES OF THE CRYSTAL SURFACE CENTER POSITIONS      
      // Rotation
      rotation_mtx[a]->Multiplication(crystal_center_ref[c], crystal_center_out);

      // Get crystal surface positions
      mp_crystalCentralPositionX[cryID] = crystal_center_out->GetMatriceElt(0,0);
      mp_crystalCentralPositionY[cryID] = crystal_center_out->GetMatriceElt(1,0);
      mp_crystalCentralPositionZ[cryID] = crystal_center_out->GetMatriceElt(2,0);

      // COMPUTE FOCAL POSITIONS AND SET CARTESIAN COORDINATES OF THE CRYSTAL FOCAL POSITIONS
      if(ComputeFocalPositions(crystal_center_ref[c]->GetMatriceElt(0,0), 
                               crystal_center_ref[c]->GetMatriceElt(1,0), 
                               crystal_center_ref[c]->GetMatriceElt(2,0), 
                                                                    hId, 
                                                                  cryID) )
      {
        Cerr("***** iScannerSPECTConv::ComputeLUT() -> An error occurred while computing the focal positions ! " << endl);
        return 1;
      }
      
      // Set elements of the focal matrix
      focal_projection_position_mtx->SetMatriceElt(0,0,mp_crystalFocalPositionX[cryID] );
      focal_projection_position_mtx->SetMatriceElt(1,0,mp_crystalFocalPositionY[cryID] );
      focal_projection_position_mtx->SetMatriceElt(2,0,mp_crystalFocalPositionZ[cryID] );
      focal_projection_position_mtx->SetMatriceElt(0,1,0);
      focal_projection_position_mtx->SetMatriceElt(1,1,0);
      focal_projection_position_mtx->SetMatriceElt(2,1,0);
      focal_projection_position_mtx->SetMatriceElt(0,2,0);
      focal_projection_position_mtx->SetMatriceElt(1,2,0);
      focal_projection_position_mtx->SetMatriceElt(2,2,1);

      // Rotate and compute projections of focal position for this crystal
      rotation_mtx[a]->Multiplication(focal_projection_position_mtx, focal_projection_position_mtx_output);
      mp_crystalFocalPositionX[cryID] = focal_projection_position_mtx_output->GetMatriceElt(0,0);
      mp_crystalFocalPositionY[cryID] = focal_projection_position_mtx_output->GetMatriceElt(1,0); 
      mp_crystalFocalPositionZ[cryID] = focal_projection_position_mtx_output->GetMatriceElt(2,0); 

      // GET CRYSTAL ORIENTATIONS
      mp_crystalOrientationX[cryID] = rotation_mtx[a]->GetMatriceElt(0,0);
      mp_crystalOrientationY[cryID] = rotation_mtx[a]->GetMatriceElt(0,1);
      mp_crystalOrientationZ[cryID] = 0;
    }
  }

  for(int i=0; i<m_nbOfProjections; i++)
    delete rotation_mtx[i];

  delete[] rotation_mtx;


  for (int c=0; c<m_nbCrystals; c++)
    delete crystal_center_ref[c];

  delete[] crystal_center_ref;
  delete crystal_center_out;
  delete focal_projection_position_mtx;
  delete focal_projection_position_mtx_output;

  delete[] head_angleX;
  delete[] head_angleY;
  delete[] head_angleZ;
  delete[] CORx;
  delete[] CORy;
  delete[] CORz;

  if (m_verbose>=2) Cout("  --> LUT generation completed" << endl);

  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ComputeFocalPositions
  \param a_posX : cartesian position of the crystal on the x-axis 
  \param a_posY : cartesian position of the crystal on the y-axis 
  \param a_posZ : cartesian position of the crystal on the z-axis 
  \param a_headID : head index of the crystal 
                   (required to select the related focal model parameters
                    and distance to center of rotation)
  \param a_cryID : crystal index in the LUT
  \brief Compute focal positions for a specific crystal ID
  \details Compute the focal positions using the implemented "constant", 
           "polynomial", "slanthole" focal models related to the gamma cameras
           The "custom" focal model is dedicated to user-made focal model
           and should be implemented by the user
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::ComputeFocalPositions(FLTNB a_posX, FLTNB a_posY, FLTNB a_posZ, int a_headID, int a_cryID)
{
  #ifdef CASTOR_DEBUG
  if(m_verbose>=4) Cout("iScannerSPECTConv::ComputeFocalPositions ..." << endl);
  #endif
  
  mp_crystalFocalPositionY[a_cryID] = -a_posY;

  // ===================================================================
  // AXIAL FOCAL POSITIONS
  // ===================================================================
  // constant
  if(mp_focalModelAxial[a_headID] == "constant")
  {
    mp_crystalFocalPositionZ[a_cryID] = a_posZ;
  }
  
  // polynomial
  else if(mp_focalModelAxial[a_headID] == "polynomial" ) 
  {
    if (mp_nbCoefModelAxial[a_headID] >= 1 && mp_nbCoefModelAxial[a_headID] <= 3)
    {
      FLTNB focal_posZ = 0;
      
      for(int coef=0 ; coef<mp_nbCoefModelAxial[a_headID] ; coef++)
        focal_posZ += m2p_axialFocalParameters[a_headID][coef]*abs(pow(a_posZ,coef));
        
      // Compute the projection of the focal position on the Y plane (thales) in order to be sure that the focal positions are not in the field of view (LOR would not go through the whole FOV) 
      mp_crystalFocalPositionZ[a_cryID] = -(2*mp_CORtoDetectorDistance[a_headID]-focal_posZ)*a_posZ/focal_posZ;
      
    }
    else if (mp_nbCoefModelAxial[a_headID] > 3) // Error (nb parameters >3)
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Error, the number of coeffs for the axial model should be <4 (max : polynom order 2 = 3 coeffs) ! " << endl);
      return 1;
    }
    else // Error (nb parameters == 0)
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Error, the number of coeffs for the axial model should be >0 ('axial number of coef model' in geom file) ! " << endl);
      return 1;
    }
  }

  // slanthole
  else if(mp_focalModelAxial[a_headID] == "slanthole" )
  {
    if (mp_nbCoefModelAxial[a_headID] == 1) 
    {
      FLTNB cst = a_posZ - m2p_axialFocalParameters[a_headID][0]*mp_CORtoDetectorDistance[a_headID];
      
      // Compute the projection of the focal position on the Y plane
      mp_crystalFocalPositionZ[a_cryID] = -m2p_axialFocalParameters[a_headID][0]*mp_CORtoDetectorDistance[a_headID] + cst;
    } 
    else // Error (nb parameters != 1)
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Error, the number of coeffs for a slanthole model should be equal to 1 (slope) : 'axial number of coef model' in geom file ! " << endl);
      return 1;
    }
  }
  
  // custom
  else if(mp_focalModelAxial[a_headID] == "custom" )
  {
    if (mp_nbCoefModelAxial[a_headID] == 1) // Uni-Parameter
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Custom model should be implemented by the user (in iScannerSPECTConv::ComputeFocalPositions()) ! " << endl);
      return 1;
    } 
    else if (mp_nbCoefModelAxial[a_headID] > 1) // Multi-Parameters
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Custom model should be implemented by the user (in iScannerSPECTConv::ComputeFocalPositions()) ! " << endl);
      return 1;
    }
    else // Error (nb parameters == 0)
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Custom model should be implemented by the user (in iScannerSPECTConv::ComputeFocalPositions()) ! " << endl);
      return 1;
    }
  }
  else  // Error, unknown model
  {
    Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Error, current model " << mp_focalModelAxial[a_headID] << " is unknown !" << endl);
    Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Should be either 'constant' (parallel), 'polynomial', 'hyperbolic, or 'custom'" << endl);
    return 1;
  }
  
  
  
  // ===================================================================
  // TRANSAXIAL FOCAL POSITIONS
  // ===================================================================
  // constant
  if(mp_focalModelTrans[a_headID] == "constant")
  {
    mp_crystalFocalPositionX[a_cryID] = a_posX;
  }
  
  // polynomial
  else if(mp_focalModelTrans[a_headID] == "polynomial" ) 
  {
    if (mp_nbCoefModelTrans[a_headID] >= 1 && mp_nbCoefModelTrans[a_headID] <= 3)
    {
      FLTNB focal_posX = 0;
      
      for(int coef=0 ; coef<mp_nbCoefModelTrans[a_headID] ; coef++)
        focal_posX += m2p_transFocalParameters[a_headID][coef]*abs(pow(a_posX,coef));
        
      // Compute the projection of the focal position on the Y plane (thales) in order to be sure that the focal positions are not in the field of view (LOR would not go through the whole FOV) 
      mp_crystalFocalPositionX[a_cryID] = -(2*mp_CORtoDetectorDistance[a_headID]-focal_posX)*a_posX/focal_posX;
      
    } 
    else if (mp_nbCoefModelTrans[a_headID] > 3)  // Error (nb parameters >3)
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Error, the number of coeffs for the trans model should be <4 (max : polynom order 2 = 3 coeffs) ! " << endl);
      return 1;
    }
    else // Error (nb parameters == 0)
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Error, the number of coeffs for the axial model should be >0 ('trans number of coef model' in geom file) ! " << endl);
      return 1;
    }
  }
  
  // slanthole
  else if(mp_focalModelTrans[a_headID] == "slanthole" ) 
  {
    if (mp_nbCoefModelTrans[a_headID] == 1) 
    {
      FLTNB cst = a_posX - m2p_transFocalParameters[a_headID][0]*mp_CORtoDetectorDistance[a_headID];
      
      // Compute the projection of the focal position on the Y plane
      mp_crystalFocalPositionX[a_cryID] = -m2p_transFocalParameters[a_headID][0]*mp_CORtoDetectorDistance[a_headID] + cst;
      
    } 
    else // Error (nb parameters != 1
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Error, the number of coeffs for a slanthole model should be equal to 1 (slope) : 'trans number of coef model' in geom file ! " << endl);
      return 1;
    }
  }
  
  // custom
  else if(mp_focalModelTrans[a_headID] == "custom" )
  {
    if (mp_nbCoefModelTrans[a_headID] == 1) // Uni-Parameter
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Custom model should be implemented by the user (in iScannerSPECTConv::ComputeFocalPositions()) ! " << endl);
      return 1;
    } 
    else if (mp_nbCoefModelTrans[a_headID] > 1) // Multi-Parameters
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Custom model should be implemented by the user (in iScannerSPECTConv::ComputeFocalPositions()) ! " << endl);
      return 1;
    }
    else // Error (nb parameters == 0)
    {
      Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Custom model should be implemented by the user (in iScannerSPECTConv::ComputeFocalPositions()) ! " << endl);
      return 1;
    }
  }
  else  // Error, unknown model
  {
    Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Error, current model " << mp_focalModelTrans[a_headID] << " is unknown !" << endl);
    Cerr("***** iScannerSPECTConv::ComputeFocalPositions() -> Should be either 'constant' (parallel), 'polynomial', 'hyperbolic, or 'custom'" << endl);
    return 1;
  }
  
  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn GetPositionsAndOrientations
  \param a_index1 : 1st index of the event (projection angle)
  \param a_index2 : 2nd index of the event (crystal index in the gamma camera)
  \param ap_Position1[3] : x,y,z cartesian position of the focal point
  \param ap_Position2[3] : x,y,z cartesian position of the crystal
  \param ap_Orientation1[3] : return -1 by default (no orientation components required for the focal point)
  \param ap_Orientation2[3] : x,y,z components of the orientation vector related to the crystal
  \param ap_POI1 : ignored in SPECT (no POI component required for the focal point)
  \param ap_POI2 : x,y,z components of the Point Of Interation related to the crystal.
                   Currently ignored (POI management not implemented in SPECT)
  \brief Get the central positions and orientations of the scanner elements from their indices.
  \todo : Point Of Interactions management is NOT implemented and ignored
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::GetPositionsAndOrientations( int a_index1, int a_index2,
                                               FLTNB ap_Position1[3], FLTNB ap_Position2[3],
                                               FLTNB ap_Orientation1[3], FLTNB ap_Orientation2[3],
                                               FLTNB* ap_POI1, FLTNB* ap_POI2 )
{
  #ifdef CASTOR_VERBOSE
  if(m_verbose>=4) Cout("iScannerSPECTConv::GetPositionsAndOrientations ..." << endl);
  #endif

  // SPECT indexes : 1st for the projection, 2nd for the crystal

  // First check projection angle existency
  if (a_index1<0 || a_index1>=m_nbOfProjections)
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> Projection index (" << a_index1 << ") out of range [0:" << m_nbOfProjections-1 << "] !" << endl);
    return 1;
  }
  
  // Second check crystals existency
  if (a_index2<0 || a_index2>=m_nbCrystals)
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> Crystal index (" << a_index2 << ") out of range [0:" << m_nbCrystals-1 << "] !" << endl);
    return 1;
  }

  // Get crystal index related to the projection
  int index = a_index1*m_nbCrystals + a_index2;
  
  // Get position of the focal point related to the crystal
  ap_Position1[0] = mp_crystalFocalPositionX[index];
  ap_Position1[1] = mp_crystalFocalPositionY[index];
  ap_Position1[2] = mp_crystalFocalPositionZ[index];

  // todo : Due to the current implementation of SPECT projection, POI 
  //        and DOI are not handled and ignored.
  //        An error is returned if POI are provided

  // Case when POI is not provided
  if (ap_POI2==NULL)
  {
    //FLTNB depth = mp_meanDepthOfInteraction[GetLayer(a_index2)] - mp_sizeCrystalDepth[GetLayer(a_index2)]/2;
    //ap_Position2[0] = mp_crystalCentralPositionX[a_index2] + depth*mp_crystalOrientationX[a_index2];
    //ap_Position2[1] = mp_crystalCentralPositionY[a_index2] + depth*mp_crystalOrientationY[a_index2];
    //ap_Position2[2] = mp_crystalCentralPositionZ[a_index2] + depth*mp_crystalOrientationZ[a_index2];
    ap_Position2[0] = mp_crystalCentralPositionX[index];
    ap_Position2[1] = mp_crystalCentralPositionY[index];
    ap_Position2[2] = mp_crystalCentralPositionZ[index];
  }
  // Case when POI[2] is negative (meaning we only have POI[0] or POI[1] specified and to be taken into account)
  else if (ap_POI2[2]<0.)
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> POI management not implemented yet for SPECT !" << endl);
    return 1;
  }
  // Case when only the DOI is provided
  else if (ap_POI2[0]==0. && ap_POI2[1]==0.)
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> POI management not implemented yet for SPECT !" << endl);
    return 1;
  }
  // Case when the full POI is taken into account
  else
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> POI management not implemented yet for SPECT !" << endl);
    return 1;
  }

  // Get orientations
  ap_Orientation1[0] = -mp_crystalOrientationX[a_index2];
  ap_Orientation1[1] = -mp_crystalOrientationY[a_index2];
  ap_Orientation1[2] = -mp_crystalOrientationZ[a_index2];
  ap_Orientation2[0] = mp_crystalOrientationX[a_index2];
  ap_Orientation2[1] = mp_crystalOrientationY[a_index2];
  ap_Orientation2[2] = mp_crystalOrientationZ[a_index2];

  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn GetRdmPositionsAndOrientations()
  \param a_index1 : 1st index of the event (projection angle)
  \param a_index2 : 2nd index of the event (crystal index in the gamma camera)
  \param ap_Position1[3] : x,y,z cartesian position of the focal point
  \param ap_Position2[3] : x,y,z cartesian position of the crystal
  \param ap_Orientation1[3] : return -1 by default (no orientation components required for the focal point)
  \param ap_Orientation2[3] : x,y,z components of the orientation vector related to the crystal
  \brief Get the focal point and random positions on the crystal surface and its orientations from the event indices.
  \details - Computed the LUT index described by the projection angle and crystal index passed in parameters. 
           - Compute random positions on the surface of the crystal
           - Write the corresponding random cartesian coordinates in the positions parameters.
  \todo This implementation has to be checked and adapted,
        as the current implementation of SPECT projection assumes crystal position located on the center of the crystal surface
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::GetRdmPositionsAndOrientations(int a_index1, int a_index2,
                                                FLTNB ap_Position1[3], FLTNB ap_Position2[3],
                                                FLTNB ap_Orientation1[3], FLTNB ap_Orientation2[3] )
{
  #ifdef CASTOR_VERBOSE
  if(m_verbose>=4) Cout("iScannerSPECTConv::GetRdmPositionsAndOrientations ..." << endl);
  #endif

  // SPECT indexes : 1st for the projection, 2nd for the crystal

  // First check projection angle existency
  if (a_index1<0 || a_index1>=m_nbOfProjections)
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> Projection index (" << a_index1 << ") out of range [0:" << m_nbOfProjections-1 << "] !" << endl);
    return 1;
  }
  
  // Second check crystals existency
  if (a_index2<0 || a_index2>=m_nbCrystals)
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> Crystal index (" << a_index2 << ") out of range [0:" << m_nbCrystals-1 << "] !" << endl);
    return 1;
  }
  
  // Get crystal index related to the projection
  int index = a_index1*m_nbCrystals + a_index2;

  // Get position of the focal point related to the crystal
  ap_Position1[0] = mp_crystalFocalPositionX[index];
  ap_Position1[1] = mp_crystalFocalPositionY[index];
  ap_Position1[2] = mp_crystalFocalPositionZ[index];
  
  // Get instance of random number generator
  sRNG* p_RNG = sRNG::GetInstance(); 

  // Get random numbers for the first crystal
  FLTNB axial = (p_RNG->GenerateRdmNber()-0.5) * m_pixelsSizeAxial;
  FLTNB trans = (p_RNG->GenerateRdmNber()-0.5) * m_pixelsSizeTrans;
  // Do not consider random depth (position on the surface of the crystal)
  //FLTNB depth = (p_RNG->GenerateRdmNber()-0.5) * m_crystalDepth;
  
  ap_Position2[0] = mp_crystalCentralPositionX[index] + trans*mp_crystalOrientationY[index] + axial*mp_crystalOrientationX[index]*mp_crystalOrientationZ[index];
  ap_Position2[1] = mp_crystalCentralPositionY[index] + trans*mp_crystalOrientationX[index] + axial*mp_crystalOrientationY[index]*mp_crystalOrientationZ[index];
  ap_Position2[2] = mp_crystalCentralPositionZ[index] + axial*sqrt(1-mp_crystalOrientationZ[index]*mp_crystalOrientationZ[index]);

  // Get orientations
  ap_Orientation1[0] = -1.;
  ap_Orientation1[1] = -1.;
  ap_Orientation1[2] = -1.;
  ap_Orientation2[0] = mp_crystalOrientationX[a_index2];
  ap_Orientation2[1] = mp_crystalOrientationY[a_index2];
  ap_Orientation2[2] = mp_crystalOrientationZ[a_index2];
  
  return 0;
}


  
// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn GetPositionWithRandomDepth
  \param a_index1 :
  \param a_index2 : index of the crystal
  \param ap_Position1[3] :
  \param ap_Position2[3] : x,y,z cartesian position of the point related to the crystal
  \brief Get the positions and orientations of scanner elements from their indices, with a random depth.
  \todo Not yet implemented
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::GetPositionWithRandomDepth(int a_index1, int a_index2, FLTNB ap_Position1[3], FLTNB ap_Position2[3])
{
  #ifdef CASTOR_VERBOSE
  if(m_verbose>=4) Cout("iScannerSPECTConv::GetPositionWithRandomDepth ..." << endl);
  #endif
  
  // This function was first implemented for PET testing purposed. Not implemented yet.
  Cerr("***** iScannerSPECTConv::GetPositionWithRandomDepth() -> This function was implemented for PET testing purposes. Not implemented for SPECT !" << endl);
  return 1;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn GetTwoCorners
  \param a_index1 : index of the projection angle
  \param a_index2 : index of the  crystal 
  \param ap_CornerInf1[3]
  \param ap_CornerSup1[3]
  \param ap_CornerInf2[3]
  \param ap_CornerSup2[3]
  \brief Get the cartesian coordinaters of the two opposite corners of a scanner element.
  \todo Not yet implemented
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::GetTwoCorners(int a_index1, int a_index2,
                               FLTNB ap_CornerInf1[3], FLTNB ap_CornerSup1[3],
                               FLTNB ap_CornerInf2[3], FLTNB ap_CornerSup2[3])
{
  #ifdef CASTOR_VERBOSE
  if(m_verbose>=4) Cout("iScannerSPECTConv::GetTwoCorners ..." << endl);
  #endif
  
  // First check projection angle existency
  if (a_index1<0 || a_index1>=m_nbOfProjections)
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> Projection index (" << a_index1 << ") out of range [0:" << m_nbOfProjections-1 << "] !" << endl);
    return 1;
  }
  
  // Second check crystals existency
  if (a_index2<0 || a_index2>=m_nbCrystals)
  {
    Cerr("***** iScannerSPECTConv::GetPositionsAndOrientations() -> Crystal index (" << a_index2 << ") out of range [0:" << m_nbCrystals-1 << "] !" << endl);
    return 1;
  }
  
  Cerr("***** iScannerSPECTConv::GetTwoCorners() -> Not implemented yet !" << endl);
  return 1;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn GetGeometricInfoFromDatafile
  \param a_pathToDF : string containing the path to datafile header
  \brief Recover geometric informations specific to the scanner class from the datafile header
  \details -Recover nb of bins and projections
           -Recover the projection angles from the header 
           (directly read from the datafile, or extrapolated from a first and last angle)
           -Recover the distance between the gamma camera detector surfaces and the center of rotation
           (directly read from the datafile, or extracted from the gamma camera configuratino file)
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::GetGeometricInfoFromDatafile(string a_pathToDF)
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::GetGeometricInfoFromDatafile ..." << endl);
  
  if (ReadDataASCIIFile(a_pathToDF , "Number of bins", mp_nbOfBins, 2, KEYWORD_OPTIONAL) == 1)
  {
    Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> Error while reading number of bins in the header data file " << endl);
    return 1;
  }

  if (ReadDataASCIIFile(a_pathToDF , "Number of projections", &m_nbOfProjections, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> Error while reading number of projections in the header data file " << endl);
    return 1;
  }

  string rotation_direction = "";
  if (ReadDataASCIIFile(a_pathToDF , "Head rotation direction", &rotation_direction, 1, KEYWORD_OPTIONAL) == 1)
  {
    Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> Error while reading head rotation orientation in the header data file " << endl);
    return 1;
  }
 
  if (PROJ_SetSPECTRotDirection(rotation_direction) )
  {
    Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> Error occured while trying to initialize head rotation orientation " << endl);
    return 1;
  }
  
  // Allocate a one-dimensional vector to retrieve angles from the datafile, then copy it in the member variables
  FLTNB* angles = new FLTNB[m_nbOfProjections];
  FLTNB first_and_last_angles[2] = {-1.,-1.};

  // Two possible initializations for projection angles :
  // - All angles are provided with the keyword "Angles"
  // - The first and last angles are provided, and the intermediate angles are extrapolated from the number of projections
  
  // 'First/Last angles' tags : checking issue during data reading/conversion (==1)
  if (ReadDataASCIIFile(a_pathToDF, "First and last projection angles", first_and_last_angles, 2, KEYWORD_OPTIONAL) == 1)
  {
    Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> Error while reading Angle mandatory field in the header data file '" << endl);
    return 1;
  }
  
  // Check for 'Projection angles' tag
  int rvalue = ReadDataASCIIFile(a_pathToDF, "Projection angles", angles, m_nbOfProjections, KEYWORD_OPTIONAL);
  
  // Error while reading "Angles" tag (==1)
  if(rvalue==1)
  {
    Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> Error while reading Angles field in the header data file !'" << endl);
    return 1;
  }

  // Check if information on projection angles has been provided
  if ( rvalue>=2 &&  // "Angles" tag not found
      (first_and_last_angles[0] <0 || first_and_last_angles[1] <0) ) // Tags first/last angles not found)
  {
    Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> No information on projection angles provided in the datafile !'" << endl);
    Cerr("                                                           This information should be provided using either the 'Angles' tag, or both 'First angles', 'Last angles' tags !'" << endl);
    return 1;
  }
  else if (rvalue>=2) // "Angles" tag not found, but first angle and last angle provided
  {
    for(int a=0 ; a<m_nbOfProjections ; a++)
      angles[a] = first_and_last_angles[0] + a*(first_and_last_angles[1]-first_and_last_angles[0]) / (m_nbOfProjections-1);
  }
  // else : Angles have been recovered using ReadDataASCIIFile() above


  // Instanciate here the projection angles variable
  mp_projectionAngles = new FLTNB[m_nbOfProjections];
    
  for(int a=0 ; a<m_nbOfProjections ; a++)
    mp_projectionAngles[a] = angles[a];
    
    
    
  // Recover distance between the detectors and scanner center of rotation
  // Allocate with the number of projections by default
  mp_CORtoDetectorDistance = new FLTNB[m_nbOfProjections]; 
  
  // Initialize by default with the scanner radius
  for(int hId=0 ; hId<m_nbHeads ; hId++)
    for(int a=0 ; a<m_nbOfProjections/m_nbHeads ; a++)
      mp_CORtoDetectorDistance[a+hId*m_nbOfProjections/m_nbHeads] = mp_radius[hId];
  
  // Read datafile value if any. First case : we have a radius specific to each projections
  int read_flag = 0;
  read_flag = ReadDataASCIIFile(a_pathToDF, "Distance camera surface to COR", mp_CORtoDetectorDistance, m_nbOfProjections, KEYWORD_OPTIONAL);
  
  if(read_flag==1)
  {
    // Error during reading
    Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> Error while reading the distance between the camera detectors to the center of rotation in the header data file " << endl);
    return 1;
  } 
  // Check the second case : we have a global radius for each projection
  else if (read_flag==2)
  {
    read_flag = ReadDataASCIIFile(a_pathToDF, "Global distance camera surface to COR", mp_CORtoDetectorDistance, 1, KEYWORD_OPTIONAL);

    if(read_flag==0) 
    {
      // Field was found : Initialize the distance for each projection angle with the global value
      for(int a=1 ; a<m_nbOfProjections ; a++)
        mp_CORtoDetectorDistance[a] = mp_CORtoDetectorDistance[0];
    }
    else if(read_flag==1)
    {
      // Error during reading
      Cerr("***** iScannerSPECTConv::GetGeometricInfoFromDatafile() -> Error while reading the global distance between the camera detectors to the center of rotation in the header data file " << endl);
      return 1;
    }
    else if (read_flag==2)
    {
      // Initialization with default values from the scanner file
      for(int hId=0 ; hId<m_nbHeads ; hId++)
        for(int a=0 ; a<m_nbOfProjections/m_nbHeads ; a++)
          mp_CORtoDetectorDistance[a+hId*m_nbOfProjections/m_nbHeads] = mp_radius[hId];
    }
  }

  delete[] angles;
  
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn GetSPECTSpecificParameters
  \param ap_nbOfProjections : total number of views
  \param ap_nbHeads : number of heads in the system
  \param ap_nbOfBins : 2 elements array containing transaxial/axial number of pixels
  \param ap_pixSizeXY : 2 elements array containing transaxial/axial pixel sizes
  \param ap_angles : Array containing angles of each projection view
  \param ap_CORtoDetectorDistance : Radius (distance between FOV center and detector)
  \param ap_headRotDirection : head rotation direction
  \brief Set pointers passed in argument with the related SPECT specific variables
         This function is used to recover these values in the datafile object
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::GetSPECTSpecificParameters(uint16_t* ap_nbOfProjections, 
                                                  uint16_t* ap_nbHeads, 
                                                  uint16_t* ap_nbOfBins,
                                                   FLTNB*  ap_pixSizeXY, 
                                                   FLTNB*& ap_angles, 
                                                   FLTNB*& ap_CORtoDetectorDistance,
                                                   int*    ap_headRotDirection)
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::GetSPECTSpecificParameters ..." << endl);
    
  // Verify that all parameters have been correctly checked
  if(m_allParametersChecked == false)
  {
    Cerr("***** iScannerSPECTConv::GetSPECTSpecificParameters() -> Parameters have not been checked !" << endl);
    return 1;
  }

  *ap_nbOfProjections = m_nbOfProjections;
  *ap_nbHeads = m_nbHeads;
  ap_nbOfBins[0] = mp_nbOfBins[0];
  ap_nbOfBins[1] = mp_nbOfBins[1];
  ap_pixSizeXY[0] = m_vPixelsSizeTrans;
  ap_pixSizeXY[1] = m_vPixelsSizeAxial;
  ap_angles = mp_projectionAngles;
  ap_CORtoDetectorDistance = mp_CORtoDetectorDistance;
  *ap_headRotDirection = m_rotDirection;
  
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_SetSPECTAngles
  \param ap_projectionAngles : an array containing the projection angles
  \brief Set the projection angles with the array provided in parameter
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::PROJ_SetSPECTAngles(FLTNB* ap_projectionAngles)
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::PROJ_SetSPECTAngles ..." << endl);

  // Check initialization of the number of projections
  if(m_nbOfProjections <= 0)
  {
    Cerr("***** iScannerSPECTConv::PROJ_SetSPECTAngles -> Error number of projection should be >0 ! '" << endl);
    return 1;
  }
  else
  {
    mp_projectionAngles = new FLTNB[m_nbOfProjections];
    
    for(int a=0 ; a<m_nbOfProjections ; a++)
      mp_projectionAngles[a] = ap_projectionAngles[a];
  }  
  
  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_SetSPECTCORtoDetectorDistance
  \param a_distance
  \brief Set distance between the center of rotation and SPECT detectors if arg value>0, 
         Set with the geometric information in the scanner configuration file otherwise
  \return 0 if success, positive value otherwise
*/
int iScannerSPECTConv::PROJ_SetSPECTCORtoDetectorDistance(FLTNB a_distance)
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::PROJ_SetSPECTCORtoDetectorDistance ..." << endl);
  
  // Check initialization of the number of projections
  if(m_nbOfProjections <= 0)
  {
    Cerr("***** iScannerSPECTConv::PROJ_SetSPECTCORtoDetectorDistance -> Error number of projection should be >0 ! '" << endl);
    return 1;
  }
  else
  {
    mp_CORtoDetectorDistance = new FLTNB[m_nbOfProjections];
    
    if(a_distance>0)
    {
      // Set all radius to the provided distance
      for(int a=0 ; a<m_nbOfProjections ; a++)
        mp_CORtoDetectorDistance[a] = a_distance;
    }
    else
    {
      // Set all distance according to scanner geometric informations (default)
      for(int hId=0 ; hId<m_nbHeads ; hId++)
        for(int a=0 ; a<m_nbOfProjections/m_nbHeads ; a++)
          mp_CORtoDetectorDistance[a+hId*m_nbOfProjections/m_nbHeads] = mp_radius[hId];
    }
  }
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_GetSPECTNbBins
  \brief Get the number of SPECT heads in the pointer provided in parameter
  \return 0 by default (no error)
*/
int iScannerSPECTConv::PROJ_GetSPECTNbBins(uint16_t* ap_nbOfBins) 
{
  ap_nbOfBins = mp_nbOfBins;
  return 0;
}



// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_SetSPECTNbBins
  \param ap_nbOfBins
  \brief Set number of bins
  \return 0 by default (no error)
*/
int iScannerSPECTConv::PROJ_SetSPECTNbBins(uint16_t* ap_nbOfBins)
{
  for(int i=0 ; i<2 ; i++)
    mp_nbOfBins[i]=ap_nbOfBins[i]; 
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_SetSPECTNbProjections
  \param a_nbOfProjections
  \brief Set number of projections
  \return 0 by default (no error)
*/
int iScannerSPECTConv::PROJ_SetSPECTNbProjections(uint32_t a_nbOfProjections)
{
  m_nbOfProjections = a_nbOfProjections; 
  return 0;
}




// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn PROJ_SetSPECTRotOrientation
  \param a_rotDirection
  \brief Set head rotation orientation
  \details This function is surcharged by the SPECT scanner daughter classes
           Returns an error by default.
  \return 1 (error) if not surcharged by a daughter class
*/
int iScannerSPECTConv::PROJ_SetSPECTRotDirection( string a_rotDirection )
{
  if(a_rotDirection == "CCW" ||
          a_rotDirection == "Ccw" ||
          a_rotDirection == "ccw" )
    m_rotDirection = GEO_ROT_CCW ;
    
  else if(a_rotDirection == ""   || // Default
          a_rotDirection == "CW" ||
          a_rotDirection == "Cw" ||
          a_rotDirection == "cw" )
    m_rotDirection = GEO_ROT_CW ;
    
  else
  {
    Cerr("***** iScannerSPECTConv::PROJ_SetSPECTRotDirection -> Error while initializing head rotation direction !" << endl);
    Cerr("     "<< a_rotDirection <<"' is unknown. Direction must be  'CW' (clockwise) or 'CCW' (counter-clockwise).");
    return 1;
  }
  
  return 0;
}


    
    
// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================
/*
  \fn ShowHelp
  \brief Display help
  \todo Provide informations about SPECT system initialization ?
*/
void iScannerSPECTConv::ShowHelp()
{
  if(m_verbose>=2) Cout("iScannerSPECTConv::ShowHelp ..." << endl);
  
  cout << "This scanner class is dedicated to the description of parallel, convergent and multi-convergent SPECT systems." << endl;
}