healpix_map.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_map.h
00028  *  Copyright (C) 2003, 2004, 2005, 2006 Max-Planck-Society
00029  *  \author Martin Reinecke
00030  */
00031 
00032 #ifndef HEALPIX_MAP_H
00033 #define HEALPIX_MAP_H
00034 
00035 #include "arr.h"
00036 #include "healpix_base.h"
00037 
00038 /*! A HEALPix map of a given datatype */
00039 template<typename T> class Healpix_Map: public Healpix_Base
00040   {
00041   private:
00042     arr<T> map;
00043 
00044   public:
00045     /*! Constructs an unallocated map. */
00046     Healpix_Map () {}
00047     /*! Constructs a map with a given \a order and the ordering
00048         scheme \a scheme. */
00049     Healpix_Map (int order, Healpix_Ordering_Scheme scheme)
00050       : Healpix_Base (order, scheme), map(npix_) {}
00051     /*! Constructs a map with a given \a nside and the ordering
00052         scheme \a scheme. */
00053     Healpix_Map (int nside, Healpix_Ordering_Scheme scheme, const nside_dummy)
00054       : Healpix_Base (nside, scheme, SET_NSIDE), map(npix_) {}
00055     /*! Constructs a map from the contents of \a data and sets the ordering
00056         scheme to \a Scheme. The size of \a data must be a valid HEALPix
00057         map size. */
00058     Healpix_Map (const arr<T> &data, Healpix_Ordering_Scheme scheme)
00059       : Healpix_Base (npix2nside(data.size()), scheme, SET_NSIDE), map(data) {}
00060 
00061     /*! Deletes the old map, creates a map from the contents of \a data and
00062         sets the ordering scheme to \a scheme. The size of \a data must be a
00063         valid HEALPix map size. */
00064     void Set (arr<T> &data, Healpix_Ordering_Scheme scheme)
00065       {
00066       Healpix_Base::SetNside(npix2nside (data.size()), scheme);
00067       map.transfer(data);
00068       }
00069 
00070     /*! Deletes the old map and creates a new map  with a given \a order
00071         and the ordering scheme \a scheme. */
00072     void Set (int order, Healpix_Ordering_Scheme scheme)
00073       {
00074       Healpix_Base::Set(order, scheme);
00075       map.alloc(npix_);
00076       }
00077     /*! Deletes the old map and creates a new map  with a given \a nside
00078         and the ordering scheme \a scheme. */
00079     void SetNside (int nside, Healpix_Ordering_Scheme scheme)
00080       {
00081       Healpix_Base::SetNside(nside, scheme);
00082       map.alloc(npix_);
00083       }
00084 
00085     /*! Fills the map with \a val. */
00086     void fill (const T &val)
00087       { map.fill(val); }
00088 
00089     /*! Imports the map \a orig into the current map, adjusting the
00090         ordering scheme. \a orig must have the same resolution as the
00091         current map. */
00092     void Import_nograde (const Healpix_Map<T> &orig)
00093       {
00094       planck_assert (nside_==orig.nside_,
00095         "Import_nograde: maps have different nside");
00096       if (orig.scheme_ == scheme_)
00097         for (int m=0; m<npix_; ++m) map[m] = orig.map[m];
00098       else
00099         {
00100         swapfunc swapper = (scheme_ == NEST) ?
00101           &Healpix_Base::ring2nest : &Healpix_Base::nest2ring;
00102 #pragma omp parallel
00103 {
00104         int m;
00105 #pragma omp for schedule (dynamic,5000)
00106         for (m=0; m<npix_; ++m) map[(this->*swapper)(m)] = orig.map[m];
00107 }
00108         }
00109       }
00110 
00111     /*! Imports the map \a orig into the current map, adjusting the
00112         ordering scheme and the map resolution. \a orig must have higher
00113         resolution than the current map. */
00114     void Import_upgrade (const Healpix_Map<T> &orig)
00115       {
00116       planck_assert(nside_>orig.nside_,"Import_upgrade: this is no upgrade");
00117       int fact = nside_/orig.nside_;
00118       planck_assert (nside_==orig.nside_*fact,
00119         "the larger Nside must be a multiple of the smaller one");
00120       pix2xyf to_xyf = (orig.scheme_==RING) ?
00121         &Healpix_Map::ring2xyf : &Healpix_Map::nest2xyf;
00122       xyf2pix from_xyf = (scheme_==RING) ?
00123         &Healpix_Map::xyf2ring : &Healpix_Map::xyf2nest;
00124 
00125 #pragma omp parallel
00126 {
00127       int m;
00128 #pragma omp for schedule (dynamic,5000)
00129       for (m=0; m<orig.npix_; ++m)
00130         {
00131         int x,y,f;
00132         (orig.*to_xyf)(m,x,y,f);
00133         for (int j=fact*y; j<fact*(y+1); ++j)
00134           for (int i=fact*x; i<fact*(x+1); ++i)
00135             {
00136             int mypix = (this->*from_xyf)(i,j,f);
00137             map[mypix] = orig.map[m];
00138             }
00139         }
00140 }
00141       }
00142 
00143     /*! Imports the map \a orig into the current map, adjusting the
00144         ordering scheme and the map resolution. \a orig must have higher
00145         resolution than the current map.
00146         \a pessimistic determines whether or not
00147         pixels are set to \a Healpix_undef when not all of the corresponding
00148         high-resolution pixels are defined.
00149 
00150         This method is instantiated for \a float and \a double only. */
00151     void Import_degrade (const Healpix_Map<T> &orig, bool pessimistic=false);
00152 
00153     /*! Imports the map \a orig into the current map, adjusting the
00154         ordering scheme and the map resolution if necessary.
00155         When downgrading, \a pessimistic determines whether or not
00156         pixels are set to \a Healpix_undef when not all of the corresponding
00157         high-resolution pixels are defined.
00158 
00159         This method is instantiated for \a float and \a double only. */
00160     void Import (const Healpix_Map<T> &orig, bool pessimistic=false)
00161       {
00162       if (orig.nside_ == nside_) // no up/degrading
00163         Import_nograde(orig);
00164       else if (orig.nside_ < nside_) // upgrading
00165         Import_upgrade(orig);
00166       else
00167         Import_degrade(orig,pessimistic);
00168       }
00169 
00170     /*! Returns a constant reference to the pixel with the number \a pix. */
00171     const T &operator[] (int pix) const { return map[pix]; }
00172     /*! Returns a reference to the pixel with the number \a pix. */
00173     T &operator[] (int pix) { return map[pix]; }
00174 
00175     /*! Swaps the map ordering from RING to NEST and vice versa.
00176         This is done in-place (i.e. with negligible space overhead). */
00177     void swap_scheme()
00178       {
00179       static const int clen[] = { 0,7,5,4,12,10,13,18,14,19,18,17,27,21 };
00180       static const int cycle[][30] = {
00181         { },
00182         { 0,1,8,12,16,21,40 },
00183         { 0,1,2,40,114 },
00184         { 0,4,160,263 },
00185         { 0,4,30,49,51,87,526,1027,1105,1387,1807,2637 },
00186         { 0,8,10,18,39,74,146,307,452,4737 },
00187         { 0,1,2,7,9,17,80,410,1526,1921,32859,33566,38931 },
00188         { 0,5,6,10,12,24,27,95,372,494,924,1409,3492,4248,9137,66043,103369,
00189           156899 },
00190         { 0,1,2,3,4,45,125,351,697,24337,102940,266194,341855,419857 },
00191         { 0,1,2,3,9,16,1705,2082,2126,8177,12753,15410,52642,80493,83235,
00192           88387,99444,1675361,2495125 },
00193         { 0,2,6,8,9,11,20,50,93,152,183,2137,13671,44794,486954,741908,
00194           4803258,5692573 },
00195         { 0,1,5,6,44,53,470,2847,3433,4906,13654,14710,400447,1797382,
00196           2744492,18775974,23541521 },
00197         { 0,4,9,10,16,33,83,117,318,451,5759,10015,128975,171834,211256,
00198           347608,1278690,2154097,2590798,3427694,5581717,21012301,27023976,
00199           72522811,95032729,139166747,171822389 },
00200         { 0,5,10,267,344,363,2968,3159,9083,18437,76602,147614,1246902,
00201           1593138,2035574,6529391,9511830,11340287,29565945,281666026,
00202           677946848 } };
00203 
00204       swapfunc swapper = (scheme_ == NEST) ?
00205         &Healpix_Base::ring2nest : &Healpix_Base::nest2ring;
00206 
00207       planck_assert (order_>=0, "swap_scheme(): need hierarchical map");
00208 
00209       for (int m=0; m<clen[order_]; ++m)
00210         {
00211         int istart = cycle[order_][m];
00212 
00213         T pixbuf = map[istart];
00214         int iold = istart, inew = (this->*swapper)(istart);
00215         while (inew != istart)
00216           {
00217           map[iold] = map[inew];
00218           iold = inew;
00219           inew = (this->*swapper)(inew);
00220           }
00221         map[iold] = pixbuf;
00222         }
00223       scheme_ = (scheme_==RING) ? NEST : RING;
00224       }
00225 
00226     /*! performs the actual interpolation using \a pix and \a wgt. */
00227     T interpolation (const fix_arr<int,4> &pix,
00228       const fix_arr<double,4> &wgt) const
00229       {
00230       return map[pix[0]]*wgt[0] + map[pix[1]]*wgt[1]
00231            + map[pix[2]]*wgt[2] + map[pix[3]]*wgt[3];
00232       }
00233     /*! Returns the interpolated map value at \a ptg */
00234     T interpolated_value (const pointing &ptg) const
00235       {
00236       fix_arr<int,4> pix;
00237       fix_arr<double,4> wgt;
00238       get_interpol (ptg, pix, wgt);
00239       return interpolation (pix, wgt);
00240       }
00241 
00242     /*! Returns a constant reference to the map data. */
00243     const arr<T> &Map() const { return map; }
00244 
00245     /*! Returns the minimum and maximum value of the map in
00246         \a Min and \a Max.
00247 
00248         This method is instantiated for \a float and \a double only. */
00249     void minmax (T &Min, T &Max) const;
00250 
00251     /*! Swaps the contents of two Healpix_Map objects. */
00252     void swap (Healpix_Map &other)
00253       {
00254       Healpix_Base::swap(other);
00255       map.swap(other.map);
00256       }
00257 
00258     /*! Returns the average of all defined map pixels. */
00259     double average() const
00260       {
00261       double avg=0;
00262       int pix=0;
00263       for (int m=0; m<npix_; ++m)
00264         if (!approx<double>(map[m],Healpix_undef))
00265           { ++pix; avg+=map[m]; }
00266       return avg/pix;
00267       }
00268 
00269     /*! Adds \a val to all defined map pixels. */
00270     void add (T val)
00271       {
00272       for (int m=0; m<npix_; ++m)
00273         if (!approx<double>(map[m],Healpix_undef))
00274           { map[m]+=val; }
00275       }
00276 
00277     /*! Returns the root mean square of the map, not counting undefined
00278         pixels. */
00279     double rms() const
00280       {
00281       using namespace std;
00282 
00283       double result=0;
00284       int pix=0;
00285       for (int m=0; m<npix_; ++m)
00286         if (!approx<double>(map[m],Healpix_undef))
00287           { ++pix; result+=map[m]*map[m]; }
00288       return sqrt(result/pix);
00289       }
00290     /*! Returns the maximum absolute value in the map, ignoring undefined
00291         pixels. */
00292     T absmax() const
00293       {
00294       using namespace std;
00295 
00296       T result=0;
00297       for (int m=0; m<npix_; ++m)
00298         if (!approx<double>(map[m],Healpix_undef))
00299           { result = max(result,abs(map[m])); }
00300       return result;
00301       }
00302   };
00303 
00304 #endif