source: Sophya/trunk/SophyaLib/Samba/spherethetaphi.cc@ 372

//
// 
#include "spherethetaphi.h"
#include "nbmath.h"
#include <iostream.h>
//***************************************************************
//++ 
// Class        SphereThetaPhi
// 
//
// include      spherethetaphi.h nbmath.h 
//
//    Découpage de la sphère selon theta et phi, chaque 
//    hémisphère étant découpé en (m-1) parallèles (l'équateur compte pour du 
//    beurre), chacune des m bandes de theta ainsi définies étant découpée par 
//    des méridiens équirepartis, ce découpage étant fait de sorte que tous 
//    les pixels aient la même surface et soient le plus carré possible.
//    On commence par découper l'hémisphère de z positif en partant du pôle et 
//    en allant vers l'équateur. Le premier pixel est la calotte polaire, 
//    il est circulaire et centré sur theta=0.
//--
//++
//
// Links        Parents
//
//    SphericalMap 
//-- 
//++
//
// Links        Descendants
//
//     
//-- 

/* --Methode-- */
//++
// Titre        Constructeurs
//--
//++


SphereThetaPhi::SphereThetaPhi()

//--
{
  InitNul();
}

/* --Methode-- */
//++
SphereThetaPhi::SphereThetaPhi(char* flnm)

//    Constructeur : charge une image à partir d'un fichier 
//--
{
  InitNul();
  PInPersist s(flnm);
  Read(s);
}

/* --Methode-- */

//++
SphereThetaPhi::SphereThetaPhi(int_4 m, int_4 pet)

//    Constructeur : m est le nombre de bandes en theta sur un hémisphère 
//    (la calotte constituant la premiere bande).
//    pet est le nombre de pixels (pétales) de la bande en contact avec la 
//    calotte polaire. Pour l'instant pet est inopérant! 
//--
{
  InitNul();
  Pixelize(m,pet);
}


//++
// Titre        Destructeur
//--
//++
SphereThetaPhi::~SphereThetaPhi()

//--
{
  Clear();
}

//++
// Titre        Méthodes
//--
void SphereThetaPhi::InitNul()
//
//    initialise à zéro les variables de classe pointeurs
{
  mTheta = NULL;
  mNphi = NULL;
  mTNphi = NULL;
  // mPix = NULL;
  _pixel=NULL;
  mNTheta = mNPet = 0;
  mNPix = 0;
}

/* --Methode-- */
void SphereThetaPhi::Clear()

{
  if (mTheta)  delete[] mTheta;
  if (mNphi)  delete[] mNphi;
  if (mTNphi) delete[] mTNphi;
  //if (mPix)   delete[] mPix;
  if (_pixel) _pixel->Detach();
  mTheta = NULL;
  mNphi = NULL;
  mTNphi = NULL;
  //mPix = NULL;
  _pixel=NULL;
  mNTheta = mNPet = 0;
  mNPix = 0;
}

/* --Methode-- */
//++
void SphereThetaPhi::WriteSelf(POutPersist& s) const

//    créer un fichier image
//--
{
  char strg[256];
  if (mInfo_) {sprintf(strg, "SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo", (int_4)mNTheta, (int_4)mNPix);
  }
  else { sprintf(strg, "SphereThetaPhi: NSlices=%6d  NPix=%9d ", 
               (int_4)mNTheta, (int_4)mNPix);
  }
  s.PutLine(strg);
  if (mInfo_)  s << (*mInfo_);
  s.PutI4(mNTheta);
  s.PutI4(mNPet);
  s.PutI4(mNPix);
  s.PutR8(mOmeg);
  s.PutR4s(mTheta, mNTheta+1);
  s.PutI4s(mNphi, mNTheta);
  s.PutI4s(mTNphi, mNTheta+1);
  s.PutR8s(_pixel->Data(), mNPix);
  //s.PutR8s(mPix, mNPix);
  return;
}

/* --Methode-- */
//++
void SphereThetaPhi::ReadSelf(PInPersist& s)

//    relit un fichier d'image
//--
{
  Clear();
  char strg[256];
  s.GetLine(strg, 255);
// Pour savoir s'il y avait un DVList Info associe
  bool hadinfo = false;
  if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo = true;
  if (hadinfo) {    // Lecture eventuelle du DVList Info
    if (mInfo_ == NULL)  mInfo_ = new DVList;
    s >> (*mInfo_);
  }
  s.GetI4(mNTheta);
  s.GetI4(mNPet);
  s.GetI4(mNPix);
  s.GetR8(mOmeg);
  mTheta = new r_4[mNTheta+1];
  mNphi = new int_4[mNTheta];
  mTNphi = new int_4[mNTheta+1];
  //mPix = new r_8[mNPix+1];
  _pixel=new PDataArray<r_8>(mNPix,true);
  s.GetR4s(mTheta, mNTheta+1);
  s.GetI4s(mNphi, mNTheta);
  s.GetI4s(mTNphi, mNTheta+1);
  s.GetR8s(_pixel->Data(), mNPix);
  //s.GetR8s(mPix, mNPix);
  return;
}


/* --Methode-- */
//++
int_4 SphereThetaPhi::NbPixels() const

//    Retourne le nombre de pixels du découpage
//--
{
  return(mNPix);
}


static r_8 dummy_pixel = 0;

/* --Methode-- */
//++
r_8& SphereThetaPhi::PixVal(int_4 k)

//    Retourne la valeur du contenu du pixel d'indice k
//--
{
  if ( (k<0) || (k >= mNPix) ) {
    //  THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
    cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
  THROW(rangeCheckErr);

    //throw "SphereThetaPhi::PIxVal Pixel index out of range";
  }
  return (*(_pixel->Data()+k));
  //return(mPix[k]);
}
//++
r_8 const& SphereThetaPhi::PixVal(int_4 k) const

//    Retourne la valeur du contenu du pixel d'indice k
//--
{
  if ( (k<0) || (k >= mNPix) ) {
    cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
  //    THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));

    throw "SphereThetaPhi::PIxVal Pixel index out of range"; 
  } 
  return *(_pixel->Data()+k);
  //return(mPix[k]);
}


/* --Methode-- */
//++
int_4 SphereThetaPhi::PixIndexSph(r_4 theta, r_4 phi) const

//    Retourne l'indice du pixel vers lequel pointe une direction définie par 
//    ses coordonnées sphériques
//--
{
  r_4 dphi;
  int_4 i,j,k;
  bool fgzn = false;

  if( (theta > (r_4)Pi) || (theta < 0. ) ) return(-1);
  if( (phi < 0.) || (phi > DeuxPi) )  return(-1);
  if( theta > (r_4)Pi*0.5)  { fgzn = true;  theta = Pi-theta; }
  
  // La bande d'indice kt est limitée par les valeurs de theta contenues dans
  // mTheta[kt] et mTheta[kt+1]
  for( i=1; i< mNTheta; i++ ) 
    if( theta < mTheta[i] ) break;
  
  dphi= (r_4)DeuxPi/(r_4)mNphi[i-1];
  
  if (fgzn)  k= mNPix-mTNphi[i]+(int_4)(phi/dphi);
  else k= mTNphi[i-1]+(int_4)(phi/dphi); 
  
  return(k);
}


/* --Methode-- */
//++
void SphereThetaPhi::PixThetaPhi(int_4 k, r_4& theta, r_4& phi) const

//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
//--
{
  int_4 i; 
  bool fgzn = false;
  
  if( (k < 0) || (k >= mNPix) ) {theta = -99999.; phi = -99999.;  return; }
  if( k >= mNPix/2)  {fgzn = true;  k = mNPix-1-k; }

  // recupère l'indice i de la tranche qui contient le pixel k
  for( i=0; i< mNTheta; i++ ) 
    if( k < mTNphi[i+1] ) break;

  // angle theta
  theta= 0.5*(mTheta[i]+mTheta[i+1]);
  if (fgzn) theta= Pi-theta;
  
  // angle phi
  k -= mTNphi[i];
  phi= (r_4)DeuxPi/(r_4)mNphi[i]*(r_4)(k+.5);
  if (fgzn) phi= DeuxPi-phi;

  return;
}


//++
r_8         SphereThetaPhi::PixSolAngle(int_4 dummy) const
//    Pixel Solid angle  (steradians)
//    All the pixels have the same solid angle. The dummy argument is
//    for compatibility with eventual pixelizations which would not 
//    fulfil this requirement.
//--
{
  return mOmeg;
}


/* --Methode-- */
//++
void SphereThetaPhi::Limits(int_4 k, r_4& tetMin, r_4& tetMax, 
                           r_4& phiMin, r_4& phiMax) 

//    Retourne les valeurs de theta et phi limitant le pixel d'indice k
//--
  {
  int_4 j;
  r_4 dphi;
  bool fgzn= false;
  
  if( (k< 0) || (k >= mNPix) ) {
    tetMin= -99999.; 
    phiMin= -99999.;  
    tetMax= -99999.; 
    phiMax= -99999.;  
    return;
  }
  
  // si on se trouve dans l'hémisphère Sud
  if( k >= mNPix/2 ) { 
    fgzn= true; 
    k= mNPix-1-k; 
  }
  
  // on recupere l'indice i de la tranche qui contient le pixel k
  int i;
  for( i=0; i< mNTheta; i++ ) 
    if( k< mTNphi[i+1] ) break;
  
  // valeurs limites de theta dans l'hemisphere Nord
  tetMin= mTheta[i];
  tetMax= mTheta[i+1];
  // valeurs limites de theta dans l'hemisphere Sud
  if (fgzn) {
    tetMin= Pi-mTheta[i+1];
    tetMax= Pi-mTheta[i];
  }
  
  // pixel correspondant dans l'hemisphere Nord
  if (fgzn) k= mTNphi[i+1]-k+mTNphi[i]-1;
  
  // indice j de discretisation ( phi= j*dphi )
  j= k-mTNphi[i];
  dphi= (r_4)DeuxPi/(r_4)mNphi[i];
  
  // valeurs limites de phi 
  phiMin= dphi*(r_4)(j);
  phiMax= dphi*(r_4)(j+1);
  return;
}

/* --Methode-- */
//++
int_4 SphereThetaPhi::NbThetaSlices() const

//    Retourne le nombre de tranches en theta sur la sphere
//--
{
  int_4 nbslices;
  nbslices= 2*mNTheta;
  return(nbslices);
}

/* --Methode-- */
//++
int_4 SphereThetaPhi::NPhi(int_4 kt) const

//    Retourne le nombre de pixels en phi de la tranche kt
//-- 
{
  
  int_4 nbpix;
  
  // verification
  if( (kt< 0) || (kt>= 2*mNTheta) ) return(-1);
  
  // si on se trouve dans l'hemisphere Sud
  if( kt >= mNTheta ) {
    kt= 2*mNTheta-1-kt; 
  }
  
  // nombre de pixels
  nbpix= mNphi[kt];
  return(nbpix);
}


/* --Methode-- */
//++
void SphereThetaPhi::Theta(int_4 kt,r_4& tetMin,r_4& tetMax) 

//    Retourne les valeurs de theta limitant la tranche kt
//--
{
  
  bool fgzn= false;
  
  // verification
  if( (kt< 0) || (kt>= 2*mNTheta) ) {
    tetMin= -99999.; 
    tetMax= -99999.; 
    return;
  }

  // si on se trouve dans l'hemisphere Sud
  if( kt >= mNTheta ) {
    fgzn= true;
    kt= 2*mNTheta-1-kt; 
  }
  
  // valeurs limites de theta dans l'hemisphere Nord
  tetMin= mTheta[kt];
  tetMax= mTheta[kt+1];
  // valeurs limites de theta dans l'hemisphere Sud
  if (fgzn) {
    tetMin= Pi-mTheta[kt+1];
    tetMax= Pi-mTheta[kt];
  }  
}

/* --Methode-- */
//++
void SphereThetaPhi::Phi(int_4 kt,int_4 jp,r_4& phiMin,r_4& phiMax) 

//    Retourne les valeurs de phi limitant le pixel jp de la tranche kt
//--
{
  
  r_4 dphi;
  
  // verification
  if( (kt< 0) || (kt>= 2*mNTheta) ) {
    phiMin= -99999.; 
    phiMax= -99999.; 
    return;
  }
  
  // si on se trouve dans l'hemisphere Sud
  if( kt >= mNTheta ) kt= 2*mNTheta-1-kt;
  
  // verifie si la tranche kt contient au moins jp pixels
  if( (jp< 0) || (jp >= mNphi[kt]) ) {
    phiMin= -88888.; 
    phiMax= -88888.; 
    return;
  }
  
  dphi= (r_4)DeuxPi/(r_4)mNphi[kt];
  phiMin= dphi*(r_4)(jp);
  phiMax= dphi*(r_4)(jp+1);
  return;
}

/* --Methode-- */
//++
int_4 SphereThetaPhi::Index(int_4 kt,int_4 jp) const

//    Retourne l'indice du pixel d'indice jp dans la tranche kt 
//--
{
  
  int_4 k;
  bool fgzn= false;
  
  // si on se trouve dans l'hemisphere Sud
  if( kt >= mNTheta ) {
    fgzn= true; 
    kt= 2*mNTheta-1-kt;
  }
  
  // si la tranche kt contient au moins i pixels
  if( (jp>=0) && (jp<mNphi[kt]) ) {
    // dans l'hemisphere Sud
    if (fgzn) k= mNPix-mTNphi[kt+1]+jp;
    // dans l'hemisphere Nord
    else k= mTNphi[kt]+jp;
  }
  else{
    k= 9999;
    printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
  }
  return(k);
}

/* --Methode-- */
//++
void SphereThetaPhi::ThetaPhiIndex(int_4 k,int_4& kt,int_4& jp) 

//    Retourne les indices kt et jp du pixel d'indice k
//--
{
  
  bool fgzn= false;
  
  // si on se trouve dans l'hemisphere Sud
  if( k >= mNPix/2 ) { 
    fgzn= true; 
    k= mNPix-1-k; 
  }
  
  // on recupere l'indice kt de la tranche qui contient le pixel k
  int i;
  for( i=0; i< mNTheta; i++ ) 
    if( k< mTNphi[i+1] ) break;
  
  // indice  kt de tranche
  if (fgzn) kt= 2*mNTheta-1-i;
  else kt= i;
  
  // indice jp de pixel
  if (fgzn) jp= mTNphi[i+1]-k-1;
  else jp= k-mTNphi[i];
  
}
//++
void  SphereThetaPhi::Pixelize( int_4 m, int_4 pet) 

//    effectue le découpage en pixels (m et pet ont la même signification 
//    que pour le constructeur)
//
//    Chaque bande de theta sera découpée en partant de phi=0 ...
//    L'autre hémisphère est parcourue dans le même sens en phi et de
//    l'équateur vers le pôle (le pixel qui suit le dernier de la bande la plus
//    proche de l'équateur a z>0 est celui de plus petit phi de la bande la 
//    plus proche de l'equateur a z<0).
//--
{
  int_4 ntotpix,i,j;
  //r_8 x, omeg, omeg2pi, teti, tetm, tet, nphi, costeti;
  r_8 omeg2pi;
  // Decodage et controle des arguments d'appel :
  // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
  if (m < 2) m = 2;
  if (m > 16384) m = 16384;
  
  // Au moins 4 et au plus 256 pixels dans la premiere bande decoupee en phi
  if (pet < 4) pet = 4;
  if (pet > 256) pet = 256;
  
  // On memorise les arguments d'appel
  mNTheta = m;  mNPet = pet;
  
  // On commence par decouper l'hemisphere z>0.
  // Creation des vecteurs contenant :
  // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
  mTheta = new r_4[m+1];
  
  // Le nombre de pixels en phi de chacune des bandes en theta
  mNphi = new int_4[m];
  
  // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
  // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
  mTNphi = new int_4[m+1];
  
  // Calcul du nombre total de pixels dans chaque bande pour optimiser
  // le rapport largeur/hauteur des pixels 
  
  //calotte polaire
  mTNphi[0]=0;
  mNphi[0] = 1;
  
  //bandes jusqu'a l'equateur
  for(j=1; j < m; j++) {
    mTNphi[j] = mTNphi[j-1]+mNphi[j-1];
    mNphi[j]  = (int_4)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
  }
  
  // Nombre total de pixels sur l'hemisphere
  ntotpix = mTNphi[m-1]+mNphi[m-1];
  mTNphi[m] = ntotpix;
  // et sur la sphere entiere
  mNPix=2*ntotpix;
  
  // Creation et initialisation du vecteur des contenus des pixels 
  //mPix = new r_8[mNPix];
  //for(i=0; i<mNPix; i++)  mPix[i] = 0.;
  _pixel=new PDataArray<r_8>(mNPix,true);
  for(i=0; i<mNPix; i++)  *(_pixel->Data()+i)=0.;
  // Determination des limites des bandes en theta :
  // omeg est l'angle solide couvert par chaque pixel,
  // une bande donnee kt couvre un angle solide mNphi[kt]*omeg
  // egal a 2* Pi*(cos mTheta[kt]-cos mTheta[kt+1]). De meme, l'angle solide 
  //de la
  // calotte allant du pole a la limite haute de la bande kt vaut
  // 2* Pi*(1.-cos mTheta[kt+1])= mTNphi[kt]*omeg...
  
  omeg2pi = 1./(r_8)ntotpix;
  mOmeg = omeg2pi*DeuxPi;
  
  for(j=0; j <= m; j++) {
    mTheta[j] = acos(1.-(double)mTNphi[j]*omeg2pi);
  }
}


//++
void  SphereThetaPhi::GetThetaSlice(int_4 index, r_4& theta, Vector& phi, Vector& value) const 

//    Retourne, pour la tranche en theta d'indice 'index' le theta 
//    correspondant, un vecteur (Peida) contenant les phi des pixels de 
//    la tranche, un vecteur (Peida) contenant les valeurs de pixel 
//    correspondantes
//--

{
  cout << "entree GetThetaSlice, couche no " << index << endl;
  if (index<0 || index > NbThetaSlices()) {
    // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 
    cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
  THROW(rangeCheckErr);
  }
  int_4 iring=Index(index,0);
  int_4 bid=this->NPhi(index);
  int lring=bid;
  cout << " iring= " << iring << " lring= " << lring << endl;
  phi.Realloc(lring);
  value.Realloc(lring);
  float T=0.;
  float F=0.;
  for (int kk=0; kk<lring;kk++) {
    PixThetaPhi(kk+iring,T,F);
    phi(kk)=F;
    value(kk)=PixVal(kk+iring);
  }
  theta=T;

}


/* --Methode-- */