healpix_base2.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_base2.h
00028  *  Copyright (C) 2003, 2004, 2005, 2006 Max-Planck-Society
00029  *  \author Martin Reinecke
00030  */
00031 
00032 #ifndef HEALPIX_BASE2_H
00033 #define HEALPIX_BASE2_H
00034 
00035 #include "healpix_base.h"
00036 #include "datatypes.h"
00037 
00038 /*! Functionality related to the HEALPix pixelisation. Supports resolutions up to
00039     N_side = 2^29. */
00040 class Healpix_Base2
00041   {
00042   protected:
00043     enum { order_max=29 };
00044 
00045     class Tablefiller
00046       {
00047       public:
00048         Tablefiller();
00049       };
00050     static Tablefiller Filler;
00051     friend class Tablefiller;
00052 
00053     static short ctab[0x100], utab[0x100];
00054 
00055     static const int jrll[];
00056     static const int jpll[];
00057 
00058     /*! The order of the map; -1 for nonhierarchical map. */
00059     int order_;
00060     /*! The N_side parameter of the map; 0 if not allocated. */
00061     int64 nside_;
00062     int64 npface_, ncap_, npix_;
00063     double fact1_, fact2_;
00064     /*! The map's ordering scheme. */
00065     Healpix_Ordering_Scheme scheme_;
00066 
00067     inline int64 ring_above (double z) const;
00068     void in_ring (int64 iz, double phi0, double dphi,
00069       std::vector<int64> &listir) const;
00070 
00071     int64 xyf2nest(int ix, int iy, int face_num) const;
00072     void nest2xyf(int64 pix, int &ix, int &iy, int &face_num) const;
00073     int64 xyf2ring(int ix, int iy, int face_num) const;
00074     void ring2xyf(int64 pix, int &ix, int &iy, int &face_num) const;
00075 
00076     typedef int64 (Healpix_Base2::*swapfunc)(int64 pix) const;
00077     typedef void (Healpix_Base2::*pix2xyf)
00078                  (int64 pix, int &x, int &y, int &f) const;
00079     typedef int64 (Healpix_Base2::*xyf2pix) (int x, int y, int f) const;
00080 
00081   public:
00082     /*! Calculates the map order from its \a N_side parameter.
00083         Returns -1 if \a nside is not a power of 2.
00084         \param nside the \a N_side parameter */
00085     static int nside2order (int64 nside)
00086       {
00087       planck_assert (nside>0, "invalid value for Nside");
00088       if ((nside)&(nside-1)) return -1;
00089       return ilog2(nside);
00090       }
00091     /*! Calculates the \a N_side parameter from the number of pixels.
00092         \param npix the number of pixels */
00093     static int64 npix2nside (int64 npix);
00094     /*! Constructs an unallocated object. */
00095     Healpix_Base2 ()
00096       : order_(-1), nside_(0), npface_(0), ncap_(0), npix_(0),
00097         fact1_(0), fact2_(0), scheme_(RING) {}
00098     /*! Constructs an object with a given \a order and the ordering
00099         scheme \a scheme. */
00100     Healpix_Base2 (int order, Healpix_Ordering_Scheme scheme)
00101       { Set (order, scheme); }
00102     /*! Constructs an object with a given \a nside and the ordering
00103         scheme \a scheme. The \a nside_dummy parameter must be set to
00104         SET_NSIDE. */
00105     Healpix_Base2 (int64 nside, Healpix_Ordering_Scheme scheme,
00106         const nside_dummy)
00107       { SetNside (nside, scheme); }
00108 
00109     /* Adjusts the object to \a order and \a scheme. */
00110     void Set (int order, Healpix_Ordering_Scheme scheme)
00111       {
00112       planck_assert ((order>=0)&&(order<=order_max), "bad order");
00113       order_  = order;
00114       nside_  = int64(1)<<order;
00115       npface_ = nside_<<order_;
00116       ncap_   = (npface_-nside_)<<1;
00117       npix_   = 12*npface_;
00118       fact2_  = 4./npix_;
00119       fact1_  = (nside_<<1)*fact2_;
00120       scheme_ = scheme;
00121       }
00122     /* Adjusts the object to \a nside and \a scheme. */
00123     void SetNside (int64 nside, Healpix_Ordering_Scheme scheme)
00124       {
00125       order_  = nside2order(nside);
00126       planck_assert ((scheme!=NEST) || (order_>=0),
00127         "SetNside: nside must be power of 2 for nested maps");
00128       nside_  = nside;
00129       npface_ = nside_*nside_;
00130       ncap_   = (npface_-nside_)<<1;
00131       npix_   = 12*npface_;
00132       fact2_  = 4./npix_;
00133       fact1_  = (nside_<<1)*fact2_;
00134       scheme_ = scheme;
00135       }
00136 
00137     /*! Returns the z-coordinate of the ring \a ring. This also works
00138         for the (not really existing) rings 0 and 4*nside. */
00139     double ring2z (int64 ring) const;
00140     /*! Returns the number of the ring in which \a pix lies. */
00141     int64 pix2ring (int64 pix) const;
00142 
00143     /*! Translates a pixel number from NEST to RING. */
00144     int64 nest2ring (int64 pix) const;
00145     /*! Translates a pixel number from RING to NEST. */
00146     int64 ring2nest (int64 pix) const;
00147     /*! Translates a pixel number from NEST to its Peano index. */
00148     int64 nest2peano (int64 pix) const;
00149     /*! Translates a pixel number from its Peano index to NEST. */
00150     int64 peano2nest (int64 pix) const;
00151 
00152     int64 ang2pix_z_phi (double z, double phi) const;
00153 
00154     /*! Returns the number of the pixel which contains the angular coordinates
00155         \a ang. */
00156     int64 ang2pix (const pointing &ang) const
00157       { return ang2pix_z_phi (cos(ang.theta), ang.phi); }
00158     /*! Returns the number of the pixel which contains the vector \a vec
00159         (\a vec is normalized if necessary). */
00160     int64 vec2pix (const vec3 &vec) const
00161       { return ang2pix_z_phi (vec.z/vec.Length(), safe_atan2(vec.y,vec.x)); }
00162 
00163     void pix2ang_z_phi (int64 pix, double &z, double &phi) const;
00164 
00165     /*! Returns the angular coordinates of the center of the pixel with
00166         number \a pix. */
00167     pointing pix2ang (int64 pix) const
00168       {
00169       double z, phi;
00170       pix2ang_z_phi (pix,z,phi);
00171       return pointing(acos(z),phi);
00172       }
00173     /*! Returns the vector to the center of the pixel with number \a pix. */
00174     vec3 pix2vec (int64 pix) const
00175       {
00176       double z, phi;
00177       pix2ang_z_phi (pix,z,phi);
00178       vec3 res;
00179       res.set_z_phi (z, phi);
00180       return res;
00181       }
00182 
00183     /*! Returns the numbers of all pixels whose centers lie within \a radius
00184         of \a dir in \a listpix.
00185         \param dir the angular coordinates of the disc center
00186         \param radius the radius (in radians) of the disc
00187         \param listpix a vector containing the numbers of all pixels within
00188                the disc
00189         \note This method is more efficient in the RING scheme. */
00190     void query_disc (const pointing &dir, double radius,
00191       std::vector<int64> &listpix) const;
00192     /*! Returns the numbers of all pixels that lie at least partially within
00193         \a radius of \a dir in \a listpix. It may also return a few pixels
00194         which do not lie in the disk at all.
00195         \param dir the angular coordinates of the disc center
00196         \param radius the radius (in radians) of the disc
00197         \param listpix a vector containing the numbers of all pixels within
00198                the disc
00199         \note This method works in both RING and NEST schemes, but is
00200           considerably faster in the RING scheme. */
00201 //FIXME: factor 1.362 not OK!
00202     void query_disc_inclusive (const pointing &dir, double radius,
00203       std::vector<int64> &listpix) const
00204         { query_disc (dir,radius+1.362*pi/(4*nside_),listpix); }
00205 
00206     /*! Returns useful information about a given ring of the map.
00207         \param ring the ring number (the number of the first ring is 1)
00208         \param startpix the number of the first pixel in the ring
00209         \param ringpix the number of pixels in the ring
00210         \param costheta the cosine of the colatitude (in radians) of the ring
00211         \param sintheta the sine of the colatitude (in radians) of the ring
00212         \param shifted if \a true, the center of the first pixel is not at
00213                \a phi=0 */
00214     void get_ring_info (int64 ring, int64 &startpix, int64 &ringpix,
00215       double &costheta, double &sintheta, bool &shifted) const;
00216     /*! Returns useful information about a given ring of the map.
00217         \param ring the ring number (the number of the first ring is 1)
00218         \param startpix the number of the first pixel in the ring
00219         \param ringpix the number of pixels in the ring
00220         \param theta the colatitude (in radians) of the ring
00221         \param shifted if \a true, the center of the first pixel is not at
00222                \a phi=0 */
00223     void get_ring_info2 (int64 ring, int64 &startpix, int64 &ringpix,
00224       double &theta, bool &shifted) const;
00225 
00226     /*! Returns the neighboring pixels of \a pix in \a result.
00227         On exit, \a result contains (in this order)
00228         the pixel numbers of the SW, W, NW, N, NE, E, SE and S neighbor
00229         of \a pix. If a neighbor does not exist (this can only be the case
00230         for the W, N, E and S neighbors), its entry is set to -1.
00231 
00232         \note This method works in both RING and NEST schemes, but is
00233           considerably faster in the NEST scheme. */
00234     void neighbors (int64 pix, fix_arr<int64,8> &result) const;
00235     /*! Returns interpolation information for the direction \a ptg.
00236         The surrounding pixels are returned in \a pix, their corresponding
00237         weights in \a wgt.
00238         \note This method works in both RING and NEST schemes, but is
00239           considerably faster in the RING scheme. */
00240     void get_interpol (const pointing &ptg, fix_arr<int64,4> &pix,
00241                        fix_arr<double,4> &wgt) const;
00242 
00243     /*! Returns the order parameter of the object. */
00244     int Order() const { return order_; }
00245     /*! Returns the \a N_side parameter of the object. */
00246     int64 Nside() const { return nside_; }
00247     /*! Returns the number of pixels of the object. */
00248     int64 Npix() const { return npix_; }
00249     /*! Returns the ordering scheme of the object. */
00250     Healpix_Ordering_Scheme Scheme() const { return scheme_; }
00251 
00252     /*! Returns \a true, if both objects have the same nside and scheme,
00253         else  \a false. */
00254     bool conformable (const Healpix_Base2 &other) const
00255       { return ((nside_==other.nside_) && (scheme_==other.scheme_)); }
00256 
00257     /*! Swaps the contents of two Healpix_Base objects. */
00258     void swap (Healpix_Base2 &other)
00259       {
00260       std::swap(order_,other.order_);
00261       std::swap(nside_,other.nside_);
00262       std::swap(npface_,other.npface_);
00263       std::swap(ncap_,other.ncap_);
00264       std::swap(npix_,other.npix_);
00265       std::swap(fact1_,other.fact1_);
00266       std::swap(fact2_,other.fact2_);
00267       std::swap(scheme_,other.scheme_);
00268       }
00269 
00270     /*! Returns the maximum angular distance (in radian) between any pixel
00271         center and its corners. */
00272     double max_pixrad() const;
00273   };
00274 
00275 #endif