healpix_base.h

00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*! \file healpix_base.h
00028  *  Copyright (C) 2003, 2004, 2005, 2006 Max-Planck-Society
00029  *  \author Martin Reinecke
00030  */
00031 
00032 #ifndef HEALPIX_BASE_H
00033 #define HEALPIX_BASE_H
00034 
00035 #include <vector>
00036 #include "cxxutils.h"
00037 #include "lsconstants.h"
00038 #include "pointing.h"
00039 template<typename T, unsigned int sz> class fix_arr;
00040 
00041 /*! The two possible ordering schemes of a HEALPix map. */
00042 typedef enum { RING, /*!< RING scheme */
00043                NEST  /*!< NESTED scheme */
00044              }  Healpix_Ordering_Scheme;
00045 
00046 class nside_dummy {};
00047 extern const nside_dummy SET_NSIDE;
00048 
00049 /*! Functionality related to the HEALPix pixelisation. */
00050 class Healpix_Base
00051   {
00052   protected:
00053     enum { order_max=13 };
00054 
00055     class Tablefiller
00056       {
00057       public:
00058         Tablefiller();
00059       };
00060     static Tablefiller Filler;
00061     friend class Tablefiller;
00062 
00063     static short ctab[0x100], utab[0x100];
00064 
00065     static const int jrll[];
00066     static const int jpll[];
00067 
00068     /*! The order of the map; -1 for nonhierarchical map. */
00069     int order_;
00070     /*! The N_side parameter of the map; 0 if not allocated. */
00071     int nside_;
00072     int npface_, ncap_, npix_;
00073     double fact1_, fact2_;
00074     /*! The map's ordering scheme. */
00075     Healpix_Ordering_Scheme scheme_;
00076 
00077     inline int ring_above (double z) const;
00078     void in_ring (int iz, double phi0, double dphi,
00079       std::vector<int> &listir) const;
00080 
00081     int xyf2nest(int ix, int iy, int face_num) const;
00082     void nest2xyf(int pix, int &ix, int &iy, int &face_num) const;
00083     int xyf2ring(int ix, int iy, int face_num) const;
00084     void ring2xyf(int pix, int &ix, int &iy, int &face_num) const;
00085 
00086     typedef int (Healpix_Base::*swapfunc)(int pix) const;
00087     typedef void (Healpix_Base::*pix2xyf)
00088                  (int pix, int &x, int &y, int &f) const;
00089     typedef int (Healpix_Base::*xyf2pix) (int x, int y, int f) const;
00090 
00091   public:
00092     /*! Calculates the map order from its \a N_side parameter.
00093         Returns -1 if \a nside is not a power of 2.
00094         \param nside the \a N_side parameter */
00095     static int nside2order (int nside)
00096       {
00097       planck_assert (nside>0, "invalid value for Nside");
00098       if ((nside)&(nside-1)) return -1;
00099       return ilog2(nside);
00100       }
00101     /*! Calculates the \a N_side parameter from the number of pixels.
00102         \param npix the number of pixels */
00103     static int npix2nside (int npix);
00104     /*! Constructs an unallocated object. */
00105     Healpix_Base ()
00106       : order_(-1), nside_(0), npface_(0), ncap_(0), npix_(0),
00107         fact1_(0), fact2_(0), scheme_(RING) {}
00108     /*! Constructs an object with a given \a order and the ordering
00109         scheme \a scheme. */
00110     Healpix_Base (int order, Healpix_Ordering_Scheme scheme)
00111       { Set (order, scheme); }
00112     /*! Constructs an object with a given \a nside and the ordering
00113         scheme \a scheme. The \a nside_dummy parameter must be set to
00114         SET_NSIDE. */
00115     Healpix_Base (int nside, Healpix_Ordering_Scheme scheme, const nside_dummy)
00116       { SetNside (nside, scheme); }
00117 
00118     /* Adjusts the object to \a order and \a scheme. */
00119     void Set (int order, Healpix_Ordering_Scheme scheme)
00120       {
00121       planck_assert ((order>=0)&&(order<=order_max), "bad order");
00122       order_  = order;
00123       nside_  = 1<<order;
00124       npface_ = nside_<<order_;
00125       ncap_   = (npface_-nside_)<<1;
00126       npix_   = 12*npface_;
00127       fact2_  = 4./npix_;
00128       fact1_  = (nside_<<1)*fact2_;
00129       scheme_ = scheme;
00130       }
00131     /* Adjusts the object to \a nside and \a scheme. */
00132     void SetNside (int nside, Healpix_Ordering_Scheme scheme)
00133       {
00134       order_  = nside2order(nside);
00135       planck_assert ((scheme!=NEST) || (order_>=0),
00136         "SetNside: nside must be power of 2 for nested maps");
00137       nside_  = nside;
00138       npface_ = nside_*nside_;
00139       ncap_   = (npface_-nside_)<<1;
00140       npix_   = 12*npface_;
00141       fact2_  = 4./npix_;
00142       fact1_  = (nside_<<1)*fact2_;
00143       scheme_ = scheme;
00144       }
00145 
00146     /*! Returns the z-coordinate of the ring \a ring. This also works
00147         for the (not really existing) rings 0 and 4*nside. */
00148     double ring2z (int ring) const;
00149     /*! Returns the number of the ring in which \a pix lies. */
00150     int pix2ring (int pix) const;
00151 
00152     /*! Translates a pixel number from NEST to RING. */
00153     int nest2ring (int pix) const;
00154     /*! Translates a pixel number from RING to NEST. */
00155     int ring2nest (int pix) const;
00156     /*! Translates a pixel number from NEST to its Peano index. */
00157     int nest2peano (int pix) const;
00158     /*! Translates a pixel number from its Peano index to NEST. */
00159     int peano2nest (int pix) const;
00160 
00161     int ang2pix_z_phi (double z, double phi) const;
00162 
00163     /*! Returns the number of the pixel which contains the angular coordinates
00164         \a ang. */
00165     int ang2pix (const pointing &ang) const
00166       { return ang2pix_z_phi (cos(ang.theta), ang.phi); }
00167     /*! Returns the number of the pixel which contains the vector \a vec
00168         (\a vec is normalized if necessary). */
00169     int vec2pix (const vec3 &vec) const
00170       { return ang2pix_z_phi (vec.z/vec.Length(), safe_atan2(vec.y,vec.x)); }
00171 
00172     void pix2ang_z_phi (int pix, double &z, double &phi) const;
00173 
00174     /*! Returns the angular coordinates of the center of the pixel with
00175         number \a pix. */
00176     pointing pix2ang (int pix) const
00177       {
00178       double z, phi;
00179       pix2ang_z_phi (pix,z,phi);
00180       return pointing(acos(z),phi);
00181       }
00182     /*! Returns the vector to the center of the pixel with number \a pix. */
00183     vec3 pix2vec (int pix) const
00184       {
00185       double z, phi;
00186       pix2ang_z_phi (pix,z,phi);
00187       vec3 res;
00188       res.set_z_phi (z, phi);
00189       return res;
00190       }
00191 
00192     /*! Returns the numbers of all pixels whose centers lie within \a radius
00193         of \a dir in \a listpix.
00194         \param dir the angular coordinates of the disc center
00195         \param radius the radius (in radians) of the disc
00196         \param listpix a vector containing the numbers of all pixels within
00197                the disc
00198         \note This method is more efficient in the RING scheme. */
00199     void query_disc (const pointing &dir, double radius,
00200       std::vector<int> &listpix) const;
00201     /*! Returns the numbers of all pixels that lie at least partially within
00202         \a radius of \a dir in \a listpix. It may also return a few pixels
00203         which do not lie in the disk at all.
00204         \param dir the angular coordinates of the disc center
00205         \param radius the radius (in radians) of the disc
00206         \param listpix a vector containing the numbers of all pixels within
00207                the disc
00208         \note This method works in both RING and NEST schemes, but is
00209           considerably faster in the RING scheme. */
00210     void query_disc_inclusive (const pointing &dir, double radius,
00211       std::vector<int> &listpix) const
00212         { query_disc (dir,radius+1.362*pi/(4*nside_),listpix); }
00213 
00214     /*! Returns useful information about a given ring of the map.
00215         \param ring the ring number (the number of the first ring is 1)
00216         \param startpix the number of the first pixel in the ring
00217         \param ringpix the number of pixels in the ring
00218         \param costheta the cosine of the colatitude (in radians) of the ring
00219         \param sintheta the sine of the colatitude (in radians) of the ring
00220         \param shifted if \a true, the center of the first pixel is not at
00221                \a phi=0 */
00222     void get_ring_info (int ring, int &startpix, int &ringpix,
00223       double &costheta, double &sintheta, bool &shifted) const;
00224     /*! Returns useful information about a given ring of the map.
00225         \param ring the ring number (the number of the first ring is 1)
00226         \param startpix the number of the first pixel in the ring
00227         \param ringpix the number of pixels in the ring
00228         \param theta the colatitude (in radians) of the ring
00229         \param shifted if \a true, the center of the first pixel is not at
00230                \a phi=0 */
00231     void get_ring_info2 (int ring, int &startpix, int &ringpix,
00232       double &theta, bool &shifted) const;
00233 
00234     /*! Returns the neighboring pixels of \a pix in \a result.
00235         On exit, \a result contains (in this order)
00236         the pixel numbers of the SW, W, NW, N, NE, E, SE and S neighbor
00237         of \a pix. If a neighbor does not exist (this can only be the case
00238         for the W, N, E and S neighbors), its entry is set to -1.
00239 
00240         \note This method works in both RING and NEST schemes, but is
00241           considerably faster in the NEST scheme. */
00242     void neighbors (int pix, fix_arr<int,8> &result) const;
00243     /*! Returns interpolation information for the direction \a ptg.
00244         The surrounding pixels are returned in \a pix, their corresponding
00245         weights in \a wgt.
00246         \note This method works in both RING and NEST schemes, but is
00247           considerably faster in the RING scheme. */
00248     void get_interpol (const pointing &ptg, fix_arr<int,4> &pix,
00249                        fix_arr<double,4> &wgt) const;
00250 
00251     /*! Returns the order parameter of the object. */
00252     int Order() const { return order_; }
00253     /*! Returns the \a N_side parameter of the object. */
00254     int Nside() const { return nside_; }
00255     /*! Returns the number of pixels of the object. */
00256     int Npix() const { return npix_; }
00257     /*! Returns the ordering scheme of the object. */
00258     Healpix_Ordering_Scheme Scheme() const { return scheme_; }
00259 
00260     /*! Returns \a true, if both objects have the same nside and scheme,
00261         else  \a false. */
00262     bool conformable (const Healpix_Base &other) const
00263       { return ((nside_==other.nside_) && (scheme_==other.scheme_)); }
00264 
00265     /*! Swaps the contents of two Healpix_Base objects. */
00266     void swap (Healpix_Base &other)
00267       {
00268       std::swap(order_,other.order_);
00269       std::swap(nside_,other.nside_);
00270       std::swap(npface_,other.npface_);
00271       std::swap(ncap_,other.ncap_);
00272       std::swap(npix_,other.npix_);
00273       std::swap(fact1_,other.fact1_);
00274       std::swap(fact2_,other.fact2_);
00275       std::swap(scheme_,other.scheme_);
00276       }
00277 
00278     /*! Returns the maximum angular distance (in radian) between any pixel
00279         center and its corners. */
00280     double max_pixrad() const;
00281   };
00282 
00283 #endif