Ignore:
Timestamp:
Jun 21, 2006, 10:38:26 AM (19 years ago)
Author:
ansari
Message:

1/ Mise en place de la prise en charge du schema de pixelisation NESTED dans
SphereHEALPix<T>
2/ petites corrections sur PixelMap<T> et cartes speheriques

Reza 21/06/2006

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/SkyMap/spherehealpix.cc

    r2973 r2978  
    6262/* --Methode-- */
    6363
    64 template<class T>
    65 SphereHEALPix<T>::SphereHEALPix() : pixels_(), sliceBeginIndex_(),
    66                                                 sliceLenght_()
    67 
     64//! Default constructor - optional pixelisation scheme parameter
     65template<class T>
     66SphereHEALPix<T>::SphereHEALPix(bool fgring) : fgring_(fgring), pixels_(),
     67                                    sliceBeginIndex_(), sliceLenght_()
    6868
    6969{
    … …  
    7171}
    7272
    73 /*! \fn   SOPHYA::SphereHEALPix::SphereHEALPix(int_4 m)
    74  
     73/*! \brief Constructor with specification of nside and optional pixelisation scheme
     74  
    7575  \param <m> : "nside" of the Healpix algorithm
     76  \param <fgring> : if true -> RING pixelisation (default), if not NESTED
    7677
    7778  The total number of pixels will be Npix =  12*nside**2
    … …  
    8182
    8283template<class T>
    83 SphereHEALPix<T>::SphereHEALPix(int_4 m)
    84 
    85 {
     84SphereHEALPix<T>::SphereHEALPix(int_4 m, bool fgring)
     85
     86{
     87  fgring_ = fgring;
    8688  InitNul();
    8789  Pixelize(m);
    … …  
    98100  nPix_ = s.nPix_;
    99101  omeg_ = s.omeg_;
     102  fgring_ = s.fgring_;
    100103  if(s.mInfo_) this->mInfo_= new DVList(*s.mInfo_);
    101104}
    … …  
    111114  nPix_ = s.nPix_;
    112115  omeg_ = s.omeg_;
     116  fgring_ = s.fgring_;
    113117  if(s.mInfo_) this->mInfo_= new DVList(*s.mInfo_);
    114118  //  CloneOrShare(s);
    … …  
    123127  nPix_ = a.nPix_;
    124128  omeg_ = a.omeg_;
     129  fgring_ = a.fgring_;
    125130  pixels_.CloneOrShare(a.pixels_);
    126131  sliceBeginIndex_.CloneOrShare(a.sliceBeginIndex_);
    … …  
    137142  nPix_ = a.nPix_;
    138143  omeg_ = a.omeg_;
     144  fgring_ = a.fgring_;
    139145  pixels_.Share(a.pixels_);
    140146  sliceBeginIndex_.Share(a.sliceBeginIndex_);
    … …  
    148154SphereHEALPix<T>& SphereHEALPix<T>::Set(const SphereHEALPix<T>& a)
    149155{
    150   if (this != &a)
    151     {
    152      
     156  if (this != &a)   {
    153157      if (a.NbPixels() < 1)
    154158        throw RangeCheckError("SphereHEALPix<T>::Set(a ) - SphereHEALPix a not allocated ! ");
    … …  
    156160      else CopyElt(a);
    157161     
    158      
    159162      if (this->mInfo_) delete this->mInfo_;
    160163      this->mInfo_ = NULL;
    161164      if (a.mInfo_) this->mInfo_ = new DVList(*(a.mInfo_));
    162     }
     165  }
    163166  return(*this);
    164167}
    … …  
    168171{
    169172  if (NbPixels() < 1)
    170     throw RangeCheckError("SphereHEALPix<T>::CopyElt(const SphereHEALPix<T>& )  - Not Allocated SphereHEALPix ! ");
     173    throw RangeCheckError("SphereHEALPix<T>::CopyElt(a)  - Not Allocated SphereHEALPix ! ");
    171174  if (NbPixels() != a.NbPixels())
    172     throw(SzMismatchError("SphereHEALPix<T>::CopyElt(const SphereHEALPix<T>&) SizeMismatch")) ;
     175    throw(SzMismatchError("SphereHEALPix<T>::CopyElt(a) SizeMismatch")) ;
    173176  nSide_= a.nSide_;
    174177  nPix_ = a.nPix_;
    175178  omeg_ = a.omeg_;
    176179  int k;
    177   for (k=0; k< nPix_; k++) pixels_(k) = a.pixels_(k);
     180  if (fgring_ != a.fgring_)
     181    for (k=0; k< nPix_; k++) pixels_(k) = a.pixels_(k);
     182  else {
     183    if (fgring_)  for (k=0; k< nPix_; k++)
     184      pixels_(k) = a.pixels_(ring2nest(nSide_, k));
     185    else  for (k=0; k< nPix_; k++)
     186      pixels_(k) = a.pixels_(nest2ring(nSide_, k));                         
     187  }
    178188  for (k=0; k< a.sliceBeginIndex_.Size(); k++) sliceBeginIndex_(k) = a.sliceBeginIndex_(k);
    179189  for (k=0; k< a.sliceLenght_.Size(); k++) sliceLenght_(k) = a.sliceLenght_(k);
    … …  
    181191}
    182192
    183 
    184 
    185 //++
    186 // Titre        Destructor
    187 //--
    188 //++
    189193template<class T>
    190194SphereHEALPix<T>::~SphereHEALPix()
    191 
    192 //--
    193 {
    194 }
    195 
    196 //++
    197 // Titre        Public Methods
    198 //--
     195{
     196}
    199197
    200198/*!  \fn   void SOPHYA::SphereHEALPix::Resize(int_4 m)
    … …  
    207205template<class T>
    208206void SphereHEALPix<T>::Resize(int_4 m)
    209 
    210 
    211207{
    212208  if ((m <= 0 && nSide_ > 0))  {
    … …  
    219215}
    220216
     217
     218//! return type of the map pixelisation : RING or NESTED
     219template<class T>
     220string SphereHEALPix<T>::TypeOfMap() const
     221{
     222  if (fgring_)  return string("RING");
     223  else return string("NESTED");
     224}
     225
    221226template<class T>
    222227void  SphereHEALPix<T>::Pixelize( int_4 m)
    223 
    224228//    prépare la pixelisation Gorski (m a la même signification
    225229//    que pour le constructeur)
    226 //
    227 //   
    228230{
    229231  if (m<=0 || m> 8192) {
    … …  
    291293{
    292294  uint_4 nbslices = uint_4(4*nSide_-1);
    293   if (nSide_<=0)
    294     {
    295       nbslices = 0;
    296       throw PException(" sphere not pixelized, NbSlice=0 ");
    297     }
     295  if (nSide_<=0)  {
     296    nbslices = 0;
     297    throw PException(" sphere not pixelized, NbSlice=0 ");
     298  }
    298299  return nbslices;
    299300}
    … …  
    338339template<class T>
    339340void SphereHEALPix<T>::GetThetaSlice(int_4 index,r_8& theta,TVector<r_8>& phi,TVector<T>& value) const
    340 
    341 {
    342 
    343   if (index<0 || index >= NbThetaSlices())
    344     {
     341{
     342  if (index<0 || index >= NbThetaSlices())  {
    345343      cout << " SphereHEALPix::GetThetaSlice : Pixel index out of range" <<endl;
    346344      throw RangeCheckError(" SphereHEALPix::GetThetaSlice()  index out of range");
    347     }
    348 
     345  }
    349346
    350347  int_4 iring= sliceBeginIndex_(index);
    351348  int_4 lring  =  sliceLenght_(index);
    352349
    353    phi.ReSize(lring);
    354    value.ReSize(lring);
     350  phi.ReSize(lring);
     351  value.ReSize(lring);
    355352
    356353  double TH= 0.;
    357354  double FI= 0.;
    358   for(int_4 kk = 0; kk < lring;kk++)
    359     {
     355  if (fgring_) {  // RING pixelisation scheme
     356    for(int_4 kk = 0; kk < lring;kk++)  {
    360357      PixThetaPhi(kk+iring,TH,FI);
    361358      phi(kk)= FI;
    362       value(kk)= PixVal(kk+iring);
    363     }
     359      value(kk)= pixels_(kk+iring);
     360    }
     361  }
     362  else {  // NESTED pixelisation scheme
     363    for(int_4 kk = 0; kk < lring;kk++)  {
     364      int kkn = ring2nest(nSide_, kk+iring);
     365      PixThetaPhi(kkn,TH,FI);
     366      phi(kk)= FI;
     367      value(kk)= pixels_(kkn);
     368    }
     369  }
    364370  theta= TH;
    365371}
    … …  
    394400  value.ReSize(lring);
    395401
    396   for(int_4 kk = 0; kk < lring;kk++)
    397     {
     402  if (fgring_) {  // RING pixelisation scheme
     403    for(int_4 kk = 0; kk < lring;kk++)  {
    398404      pixelIndices(kk)= kk+iring;
    399       value(kk)= PixVal(kk+iring);
    400     }
    401   PixThetaPhi(iring, theta, phi0);
     405      value(kk)= pixels_(kk+iring);
     406    }
     407    PixThetaPhi(iring, theta, phi0);
     408  }
     409  else {  // NESTED pixelisation scheme
     410    for(int_4 kk = 0; kk < lring;kk++)  {
     411      int_4 kkn = ring2nest(nSide_, kk+iring);
     412      pixelIndices(kk)= kkn;
     413      value(kk)= pixels_(kkn);
     414    }
     415    PixThetaPhi(ring2nest(nSide_,iring), theta, phi0);   
     416  }
    402417}
    403418
    … …  
    427442
    428443
    429 /*! \fn T& SOPHYA::SphereHEALPix::PixVal(int_4 k)
    430 
    431    \return value of  pixel with "RING" index k
    432 */
     444
     445//!   \return value of the \b k th pixel
    433446template<class T>
    434447T& SphereHEALPix<T>::PixVal(int_4 k)
    … …  
    436449{
    437450  if((k < 0) || (k >= nPix_))
    438     {
    439       throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range");
    440     }
     451      throw RangeCheckError("SphereHEALPix::PixVal() Pixel index out of range");
    441452  return pixels_(k);
    442453}
    443454
    444 /*! \fn T const& SOPHYA::SphereHEALPix::PixVal(int_4 k) const
    445 
    446    \return value of  pixel with "RING" index k
    447 */
     455//!   \return value of the \b k th pixel
    448456template<class T>
    449457T const& SphereHEALPix<T>::PixVal(int_4 k) const
    … …  
    451459{
    452460  if((k < 0) || (k >= nPix_))
    453     {
    454       throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range");
    455     }
     461    throw RangeCheckError("SphereHEALPix::PIxVal Pixel index out of range");
    456462  return *(pixels_.Data()+k);
    457463}
    … …  
    460466
    461467   \return value of  pixel with "NESTED" index k
    462 */
     468
    463469template<class T>
    464470T& SphereHEALPix<T>::PixValNest(int_4 k)
    … …  
    472478  return pixels_(nest2ring(nSide_,k));
    473479}
    474 
     480*/
    475481/*!  \fn T const& SOPHYA::SphereHEALPix::PixValNest(int_4 k) const
    476482
    477483  \return value of  pixel with "NESTED" index k
    478 */
    479484
    480485template<class T>
    … …  
    489494  return *(pixels_.Data()+pix);
    490495}
     496*/
    491497
    492498/*! \fn bool SOPHYA::SphereHEALPix::ContainsSph(double theta, double phi) const
    … …  
    508514
    509515{
    510   return ang2pix_ring(nSide_,theta,phi);
     516  if (fgring_) return ang2pix_ring(nSide_,theta,phi);
     517  else return ang2pix_nest(nSide_,theta,phi);
    511518}
    512519
    … …  
    514521
    515522 \return "NESTED" index of the pixel corresponding to direction (theta, phi).
    516  */
     523
    517524template<class T>
    518525int_4 SphereHEALPix<T>::PixIndexSphNest(double theta,double  phi) const
    … …  
    521528  return ang2pix_nest(nSide_,theta,phi);
    522529}
    523 
    524 
    525 /* --Methode-- */
    526 /*! \fn void SOPHYA::SphereHEALPix::PixThetaPhi(int_4 k,double& theta,double& phi) const
    527  \return (theta,phi) coordinates of middle of  pixel with "RING" index k
    528  */
     530*/
     531
     532//!  \return (theta,phi) coordinates of middle of  pixel with "RING" index k
    529533template<class T>
    530534void SphereHEALPix<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const
    531535{
    532   pix2ang_ring(nSide_,k,theta,phi);
    533 }
    534 
    535 /*! \fn T SOPHYA::SphereHEALPix::SetPixels(T v)
    536 
    537 Set all pixels to value v
    538 */
     536  if (fgring_) pix2ang_ring(nSide_,k,theta,phi);
     537  else pix2ang_nest(nSide_,k,theta,phi);
     538}
     539
     540//! Set all pixels to value v
    539541template <class T>
    540542T SphereHEALPix<T>::SetPixels(T v)
    … …  
    549551
    550552  \return (theta,phi) coordinates of middle of  pixel with "NESTED" index k
    551  */
    552553template<class T>
    553554void SphereHEALPix<T>::PixThetaPhiNest(int_4 k,double& theta,double& phi) const
    … …  
    556557  pix2ang_nest(nSide_,k,theta,phi);
    557558}
    558 
    559 
    560 /*! \fn int_4 SOPHYA::SphereHEALPix::NestToRing(int_4 k) const
    561 
    562    translation from NESTED index  into RING index
    563 */
     559*/
     560
     561//!  Conversion from NESTED index  into RING index
    564562template<class T>
    565563int_4 SphereHEALPix<T>::NestToRing(int_4 k) const
    566 
    567564{
    568565  return  nest2ring(nSide_,k);
    569566}
    570567
    571 
    572 /*!  \fn  int_4 SOPHYA::SphereHEALPix::RingToNest(int_4 k) const
    573 
    574  translation from  RING index  into NESTED index
    575 */
     568//!  Conversion from  RING index  into NESTED index
    576569template<class T>
    577570int_4 SphereHEALPix<T>::RingToNest(int_4 k) const
    … …  
    592585}
    593586
    594 /*! Add a constant value \b x to a SphereHEALPix */
     587//! Add a constant value \b x to a SphereHEALPix
    595588template <class T>
    596589SphereHEALPix<T>& SphereHEALPix<T>::Add(T a)
    … …  
    641634{
    642635  if (NbPixels() != a.NbPixels() )
    643     {
    644     throw(SzMismatchError("SphereHEALPix<T>::AddElt(const SphereHEALPix<T>&) SizeMismatch")) ;
    645     }
     636    throw(SzMismatchError("SphereHEALPix<T>::AddElt(a) SizeMismatch")) ;
     637  if (fgring_ != a.fgring_)
     638    throw(ParmError("SphereHEALPix<T>::AddElt(a) different pixelisation RING<>NESTED")) ;
     639
    646640  pixels_ += a.pixels_;
    647641  return (*this);
    … …  
    653647{
    654648  if (NbPixels() != a.NbPixels() )
    655     {
    656     throw(SzMismatchError("SphereHEALPix<T>::SubElt(const SphereHEALPix<T>&) SizeMismatch")) ;
    657     }
     649    throw(SzMismatchError("SphereHEALPix<T>::SubElt(a) SizeMismatch")) ;
     650  if (fgring_ != a.fgring_)
     651    throw(ParmError("SphereHEALPix<T>::SubElt(a) different pixelisation RING<>NESTED")) ;
     652
    658653  pixels_ -= a.pixels_;
    659654  return (*this);
    … …  
    665660{
    666661  if (NbPixels() != a.NbPixels() )
    667     {
    668     throw(SzMismatchError("SphereHEALPix<T>::MulElt(const SphereHEALPix<T>&) SizeMismatch")) ;
    669     }
     662    throw(SzMismatchError("SphereHEALPix<T>::MulElt(a) SizeMismatch")) ;
     663  if (fgring_ != a.fgring_)
     664    throw(ParmError("SphereHEALPix<T>::MulElt(a) different pixelisation RING<>NESTED")) ;
     665
    670666  pixels_ *= a.pixels_;
    671667  return (*this);
    … …  
    677673{
    678674  if (NbPixels() != a.NbPixels() )
    679     {
    680     throw(SzMismatchError("SphereHEALPix<T>::DivElt(const SphereHEALPix<T>&) SizeMismatch")) ;
    681     }
     675    throw(SzMismatchError("SphereHEALPix<T>::DivElt(a) SizeMismatch")) ;
     676  if (fgring_ != a.fgring_)
     677    throw(ParmError("SphereHEALPix<T>::DivElt(a) different pixelisation RING<>NESTED")) ;
    682678  pixels_ /= a.pixels_;
    683679  return (*this);
    … …  
    691687{
    692688  Show(os);
     689  os << "SphereHEALPix<T>(" << TypeOfMap() << ") NSide= "
     690     << nSide_ << " nPix_   = " << nPix_  << " omeg_ = " << omeg_   << endl;
     691
    693692  if(this->mInfo_) os << "  DVList Info= " << *(this->mInfo_) << endl;
    694   //
    695   os << " nSide_ = " << nSide_  << endl;
    696   os << " nPix_   = " << nPix_   << endl;
    697   os << " omeg_ = " << omeg_   << endl;
    698 
    699   os << " content of pixels : ";
     693  os << "... Pixel values : ";
    700694  for(int i=0; i < nPix_; i++)
    701695    {
Note: See TracChangeset for help on using the changeset viewer.