Changeset 2978 in Sophya for trunk


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

Location:
trunk/SophyaLib/SkyMap
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/SophyaLib/SkyMap/fiospherehealpix.cc

    r2869 r2978  
    8181uint_8 itab[3];
    8282is.Get(itab, 3);
     83 
     84//  Juin 2006 : NumVer 1->2 avec mise en place de HEALPix NESTED
     85//  Si itab[2]==1 --> NESTED
     86bool fgring = true;
     87if (itab[2] == 1) fgring = false;
     88
    8389// Let's Read the SphereCoordSys object  -- ATTENTIOn - $CHECK$
    8490  FIO_SphereCoordSys fio_scs( dobj_->GetCoordSys());
     
    100106  double Omega;
    101107  is.GetR8(Omega);
    102   dobj_->setParameters(nSide,nPix, Omega);
     108  dobj_->setParameters(nSide,nPix, Omega, fgring);
    103109
    104110// On lit les DataBlocks;
     
    124130//  On ecrit 3 uint_8
    125131//  0 : Numero de version,  1 : Size index,  2  reserve a l
     132//  Juin 2006 : NumVer 1->2 avec mise en place de HEALPix NESTED
     133//  totalement back-compatible , on utilise itab[2] qu'on met a 1
     134//  pour NESTED
    126135uint_8 itab[3];
    127 itab[0] = 1;
     136itab[0] = 2;
    128137itab[1] = dobj_->SizeIndex();
    129 itab[2] = 0;
     138itab[2] = (dobj_->IfRING()) ? 0 : 1;
    130139os.Put(itab, 3);
    131140// Let's write the SphereCoordSys object
     
    164173
    165174#ifdef __CXX_PRAGMA_TEMPLATES__
     175#pragma define_template FIO_SphereHEALPix<uint_2>
    166176#pragma define_template FIO_SphereHEALPix<int_4>
    167177#pragma define_template FIO_SphereHEALPix<r_8>
     
    172182#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
    173183namespace SOPHYA {
     184template class FIO_SphereHEALPix<uint_2>;
    174185template class FIO_SphereHEALPix<int_4>;
    175186template class FIO_SphereHEALPix<r_8>;
  • trunk/SophyaLib/SkyMap/pixelmap.h

    r2973 r2978  
    161161{
    162162os << "PixelMap<T>::Show() class: " << typeid(*this).name()
    163    << " T= " << DataTypeInfo<T>::getTypeName() << endl;
     163   << " T= " << DataTypeInfo<T>::getTypeName()
     164   << " TypeOfMap= " << TypeOfMap() << endl;
    164165double solang = 0.;
    165166if (NbPixels() > 0) solang = PixSolAngle(NbPixels()/2);
  • trunk/SophyaLib/SkyMap/skymapinit.cc

    r2615 r2978  
    5454
    5555  // ---------- Les SphereHEALPix ---------
     56
     57  PPRegister(FIO_SphereHEALPix<uint_2>); 
     58  DObjRegister(FIO_SphereHEALPix<uint_2>, SphereHEALPix<uint_2>);
    5659
    5760  PPRegister(FIO_SphereHEALPix<int_4>); 
  • trunk/SophyaLib/SkyMap/sphereecp.h

    r2973 r2978  
    2525//! return the size index value corresponding to resolution res (in radian)
    2626static inline int_4 ResolToSizeIndex(double res)
    27   { return (M_PI/res); }
     27  { return (int_4)((M_PI/res)+0.5); }
    2828
    2929  // Constructeur par defaut
  • 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    {
  • trunk/SophyaLib/SkyMap/spherehealpix.h

    r2973 r2978  
    5252         { return  HEALPix::ResolToSizeIndex(res); }
    5353
    54 SphereHEALPix();
    55 SphereHEALPix(int_4 m);
     54SphereHEALPix(bool fgring=true);
     55SphereHEALPix(int_4 m, bool fgring=true);
    5656SphereHEALPix(const SphereHEALPix<T>& s, bool share);
    5757SphereHEALPix(const SphereHEALPix<T>& s);
     
    9393virtual T SetPixels(T v);
    9494
    95 /*! Pixel Solid angle  (steradians)
     95/*! \brief Pixel Solid angle  (steradians)
    9696
    9797    All the pixels have the same solid angle. The dummy argument is
     
    108108
    109109virtual void Resize(int_4 m);
    110 
    111 /*!
    112 
    113 \return type of storage of the map : RING or NESTED
    114 
    115 at the moment, always RING
    116 */
    117 inline virtual string TypeOfMap() const {return string("RING");};
    118 
    119 
     110virtual string TypeOfMap() const;
     111
     112inline bool IfRING() const { return fgring_; }
     113inline bool IfNESTED() const { return ( (fgring_) ? false : true ); }
     114
     115/*
    120116virtual T& PixValNest(int_4 k);
    121117virtual T const& PixValNest(int_4 k) const;
    122 
    123118virtual int_4 PixIndexSphNest(double theta,double phi) const;
    124 
    125119virtual void PixThetaPhiNest(int_4 k,double& theta,double& phi) const;
     120*/
    126121
    127122void Pixelize(int_4);
     
    138133inline void Print(ostream& os) const { print(os); }
    139134
    140 
    141 
    142  
    143 
    144 
    145 
    146 // Operations diverses  = , +=, ...
    147 
     135//--------------------- Operations diverses  = , +=, ...
    148136
    149137SphereHEALPix<T>& Set(const SphereHEALPix<T>& a);
     
    153141// A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)
    154142
    155   //! Fill SphereHEALPix with all elements equal to \b x
     143//! Fill SphereHEALPix with all elements equal to \b x
    156144virtual SphereHEALPix<T>& SetT(T a);
    157145inline  SphereHEALPix<T>& operator = (T a) {return SetT(a);}
     
    202190void SetThetaSlices();
    203191
    204 inline void setParameters(int_4 nside, int_4 nbpixels, double solangle)
    205   {
    206     nSide_= nside;
    207     nPix_= nbpixels;
    208     omeg_= solangle;
    209   }
     192inline void setParameters(int_4 nside, int_4 nbpixels, double solangle,
     193                          bool fgring)
     194{
     195  nSide_= nside;
     196  nPix_= nbpixels;
     197  omeg_= solangle;
     198  fgring_ = fgring;
     199}
    210200
    211201// ------------- variables internes -----------------------
     
    214204int_4 nPix_;
    215205double omeg_;
     206bool fgring_;  // true -> RING pixelisation , false -> NESTED
    216207
    217208NDataBlock<T> pixels_;
    218 NDataBlock<int_4> sliceBeginIndex_;         // Rationalisation Mac. D.Y.
     209NDataBlock<int_4> sliceBeginIndex_;    // Indices in RING scheme
    219210NDataBlock<int_4> sliceLenght_;
    220211
  • trunk/SophyaLib/SkyMap/spherethetaphi.h

    r2973 r2978  
    2929//! return the size index value corresponding to resolution res (in radian)
    3030static inline int_4 ResolToSizeIndex(double res)
    31   { return (M_PI/2./res); }
     31  { return (int_4)((M_PI/2./res)+0.5); }
    3232
    3333SphereThetaPhi();
Note: See TracChangeset for help on using the changeset viewer.