Healpix C++  3.50
healpix_map.h
Go to the documentation of this file.
1 /*
2  * This file is part of Healpix_cxx.
3  *
4  * Healpix_cxx is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * Healpix_cxx is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with Healpix_cxx; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17  *
18  * For more information about HEALPix, see http://healpix.sourceforge.net
19  */
20 
21 /*
22  * Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
23  * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
24  * (DLR).
25  */
26 
27 /*! \file healpix_map.h
28  * Copyright (C) 2003-2013 Max-Planck-Society
29  * \author Martin Reinecke
30  */
31 
32 #ifndef HEALPIX_MAP_H
33 #define HEALPIX_MAP_H
34 
35 #include "arr.h"
36 #include "healpix_base.h"
37 
38 //! Healpix value representing "undefined"
39 const double Healpix_undef=-1.6375e30;
40 
41 /*! A HEALPix map of a given datatype */
42 template<typename T> class Healpix_Map: public Healpix_Base
43  {
44  private:
45  arr<T> map;
46 
47  public:
48  /*! Constructs an unallocated map. */
50  /*! Constructs a map with a given \a order and the ordering
51  scheme \a scheme. */
53  : Healpix_Base (order, scheme), map(npix_) {}
54  /*! Constructs a map with a given \a nside and the ordering
55  scheme \a scheme. */
56  Healpix_Map (int nside, Healpix_Ordering_Scheme scheme, const nside_dummy)
57  : Healpix_Base (nside, scheme, SET_NSIDE), map(npix_) {}
58  /*! Constructs a map from the contents of \a data and sets the ordering
59  scheme to \a Scheme. The size of \a data must be a valid HEALPix
60  map size. */
62  : Healpix_Base (npix2nside(data.size()), scheme, SET_NSIDE), map(data) {}
63 
64  /*! Deletes the old map, creates a map from the contents of \a data and
65  sets the ordering scheme to \a scheme. The size of \a data must be a
66  valid HEALPix map size.
67  \note On exit, \a data is zero-sized! */
68  void Set (arr<T> &data, Healpix_Ordering_Scheme scheme)
69  {
70  Healpix_Base::SetNside(npix2nside (data.size()), scheme);
71  map.transfer(data);
72  }
73 
74  /*! Deletes the old map and creates a new map with a given \a order
75  and the ordering scheme \a scheme. */
76  void Set (int order, Healpix_Ordering_Scheme scheme)
77  {
78  Healpix_Base::Set(order, scheme);
79  map.alloc(npix_);
80  }
81  /*! Deletes the old map and creates a new map with a given \a nside
82  and the ordering scheme \a scheme. */
83  void SetNside (int nside, Healpix_Ordering_Scheme scheme)
84  {
85  Healpix_Base::SetNside(nside, scheme);
86  map.alloc(npix_);
87  }
88 
89  /*! Fills the map with \a val. */
90  void fill (const T &val)
91  { map.fill(val); }
92 
93  /*! Imports the map \a orig into the current map, adjusting the
94  ordering scheme. \a orig must have the same resolution as the
95  current map. */
96  void Import_nograde (const Healpix_Map<T> &orig)
97  {
99  "Import_nograde: maps have different nside");
100  if (orig.scheme_ == scheme_)
101  for (int m=0; m<npix_; ++m) map[m] = orig.map[m];
102  else
103  {
104  swapfunc swapper = (scheme_ == NEST) ?
106 #pragma omp parallel
107 {
108  int m;
109 #pragma omp for schedule (dynamic,5000)
110  for (m=0; m<npix_; ++m) map[(this->*swapper)(m)] = orig.map[m];
111 }
112  }
113  }
114 
115  /*! Imports the map \a orig into the current map, adjusting the
116  ordering scheme and the map resolution. \a orig must have lower
117  resolution than the current map, and \a this->Nside() must be an
118  integer multiple of \a orig.Nside(). */
119  void Import_upgrade (const Healpix_Map<T> &orig)
120  {
121  planck_assert(nside_>orig.nside_,"Import_upgrade: this is no upgrade");
122  int fact = nside_/orig.nside_;
123  planck_assert (nside_==orig.nside_*fact,
124  "the larger Nside must be a multiple of the smaller one");
125 
126 #pragma omp parallel
127 {
128  int m;
129 #pragma omp for schedule (dynamic,5000)
130  for (m=0; m<orig.npix_; ++m)
131  {
132  int x,y,f;
133  orig.pix2xyf(m,x,y,f);
134  for (int j=fact*y; j<fact*(y+1); ++j)
135  for (int i=fact*x; i<fact*(x+1); ++i)
136  {
137  int mypix = xyf2pix(i,j,f);
138  map[mypix] = orig.map[m];
139  }
140  }
141 }
142  }
143 
144  /*! Imports the map \a orig into the current map, adjusting the
145  ordering scheme and the map resolution. \a orig must have higher
146  resolution than the current map, and \a orig.Nside() must be an
147  integer multiple of \a this->Nside().
148  \a pessimistic determines whether or not
149  pixels are set to \a Healpix_undef when not all of the corresponding
150  high-resolution pixels are defined.
151 
152  This method is instantiated for \a float and \a double only. */
153  void Import_degrade (const Healpix_Map<T> &orig, bool pessimistic=false);
154 
155  /*! Imports the map \a orig into the current map, adjusting the
156  ordering scheme and the map resolution if necessary.
157  When downgrading, \a pessimistic determines whether or not
158  pixels are set to \a Healpix_undef when not all of the corresponding
159  high-resolution pixels are defined.
160 
161  This method is instantiated for \a float and \a double only. */
162  void Import (const Healpix_Map<T> &orig, bool pessimistic=false)
163  {
164  if (orig.nside_ == nside_) // no up/degrading
165  Import_nograde(orig);
166  else if (orig.nside_ < nside_) // upgrading
167  Import_upgrade(orig);
168  else
169  Import_degrade(orig,pessimistic);
170  }
171 
172  /*! Returns a constant reference to the pixel with the number \a pix. */
173  const T &operator[] (int pix) const { return map[pix]; }
174  /*! Returns a reference to the pixel with the number \a pix. */
175  T &operator[] (int pix) { return map[pix]; }
176 
177  /*! Swaps the map ordering from RING to NEST and vice versa.
178  This is done in-place (i.e. with negligible space overhead). */
179  void swap_scheme()
180  {
181  swapfunc swapper = (scheme_ == NEST) ?
183 
184  arr<int> cycle=swap_cycles();
185 
186  for (tsize m=0; m<cycle.size(); ++m)
187  {
188  int istart = cycle[m];
189 
190  T pixbuf = map[istart];
191  int iold = istart, inew = (this->*swapper)(istart);
192  while (inew != istart)
193  {
194  map[iold] = map[inew];
195  iold = inew;
196  inew = (this->*swapper)(inew);
197  }
198  map[iold] = pixbuf;
199  }
200  scheme_ = (scheme_==RING) ? NEST : RING;
201  }
202 
203  /*! performs the actual interpolation using \a pix and \a wgt. */
205  const fix_arr<double,4> &wgt) const
206  {
207  double wtot=0;
208  T res=T(0);
209  for (tsize i=0; i<4; ++i)
210  {
211  T val=map[pix[i]];
212  if (!approx<double>(val,Healpix_undef))
213  { res+=T(val*wgt[i]); wtot+=wgt[i]; }
214  }
215  return (wtot==0.) ? T(Healpix_undef) : T(res/wtot);
216  }
217  /*! Returns the interpolated map value at \a ptg */
218  T interpolated_value (const pointing &ptg) const
219  {
220  fix_arr<int,4> pix;
221  fix_arr<double,4> wgt;
222  get_interpol (ptg, pix, wgt);
223  return interpolation (pix, wgt);
224  }
225 
226  /*! Returns a constant reference to the map data. */
227  const arr<T> &Map() const { return map; }
228 
229  /*! Returns the minimum and maximum value of the map in
230  \a Min and \a Max.
231 
232  This method is instantiated for \a float and \a double only. */
233  void minmax (T &Min, T &Max) const;
234 
235  /*! Swaps the contents of two Healpix_Map objects. */
236  void swap (Healpix_Map &other)
237  {
238  Healpix_Base::swap(other);
239  map.swap(other.map);
240  }
241 
242  /*! Returns the average of all defined map pixels. */
243  double average() const
244  {
245  kahan_adder<double> adder;
246  int pix=0;
247  for (int m=0; m<npix_; ++m)
248  if (!approx<double>(map[m],Healpix_undef))
249  { ++pix; adder.add(map[m]); }
250  return (pix>0) ? adder.result()/pix : Healpix_undef;
251  }
252 
253  /*! Adds \a val to all defined map pixels. */
254  void Add (T val)
255  {
256  for (int m=0; m<npix_; ++m)
257  if (!approx<double>(map[m],Healpix_undef))
258  { map[m]+=val; }
259  }
260 
261  /*! Multiplies all defined map pixels by \a val. */
262  void Scale (T val)
263  {
264  for (int m=0; m<npix_; ++m)
265  if (!approx<double>(map[m],Healpix_undef))
266  { map[m]*=val; }
267  }
268 
269  /*! Returns the root mean square of the map, not counting undefined
270  pixels. */
271  double rms() const
272  {
273  using namespace std;
274 
275  double result=0;
276  int pix=0;
277  for (int m=0; m<npix_; ++m)
278  if (!approx<double>(map[m],Healpix_undef))
279  { ++pix; result+=map[m]*map[m]; }
280  return (pix>0) ? sqrt(result/pix) : Healpix_undef;
281  }
282  /*! Returns the maximum absolute value in the map, ignoring undefined
283  pixels. */
284  T absmax() const
285  {
286  using namespace std;
287 
288  T result=0;
289  for (int m=0; m<npix_; ++m)
290  if (!approx<double>(map[m],Healpix_undef))
291  { result = max(result,abs(map[m])); }
292  return result;
293  }
294  /*! Returns \a true, if no pixel has the value \a Healpix_undef,
295  else \a false. */
296  bool fullyDefined() const
297  {
298  for (int m=0; m<npix_; ++m)
299  if (approx<double>(map[m],Healpix_undef))
300  return false;
301  return true;
302  }
303  /*! Sets all pixels with the value \a Healpix_undef to 0, and returns
304  the number of modified pixels. */
306  {
307  tsize res=0;
308  for (int m=0; m<npix_; ++m)
309  if (approx<double>(map[m],Healpix_undef))
310  { map[m]=0.; ++res; }
311  return res;
312  }
313  };
314 
315 #endif
bool fullyDefined() const
Definition: healpix_map.h:296
void SetNside(int nside, Healpix_Ordering_Scheme scheme)
Definition: healpix_map.h:83
void swap(Healpix_Map &other)
Definition: healpix_map.h:236
T absmax() const
Definition: healpix_map.h:284
void fill(const T &val)
double average() const
Definition: healpix_map.h:243
I nest2ring(I pix) const
void Set(arr< T > &data, Healpix_Ordering_Scheme scheme)
Definition: healpix_map.h:68
Healpix_Ordering_Scheme
void transfer(arrT &other)
Healpix_Map(const arr< T > &data, Healpix_Ordering_Scheme scheme)
Definition: healpix_map.h:61
std::size_t tsize
T interpolated_value(const pointing &ptg) const
Definition: healpix_map.h:218
void Import(const Healpix_Map< T > &orig, bool pessimistic=false)
Definition: healpix_map.h:162
void Scale(T val)
Definition: healpix_map.h:262
tsize replaceUndefWith0()
Definition: healpix_map.h:305
void SetNside(I nside, Healpix_Ordering_Scheme scheme)
void swap(T_Healpix_Base &other)
I ring2nest(I pix) const
void Set(int order, Healpix_Ordering_Scheme scheme)
Definition: healpix_map.h:76
void Set(int order, Healpix_Ordering_Scheme scheme)
void swap_scheme()
Definition: healpix_map.h:179
static I npix2nside(I npix)
Definition: healpix_base.cc:43
#define planck_assert(testval, msg)
const double Healpix_undef
Healpix value representing "undefined".
Definition: healpix_map.h:39
Healpix_Map(int order, Healpix_Ordering_Scheme scheme)
Definition: healpix_map.h:52
void Import_nograde(const Healpix_Map< T > &orig)
Definition: healpix_map.h:96
T interpolation(const fix_arr< int, 4 > &pix, const fix_arr< double, 4 > &wgt) const
Definition: healpix_map.h:204
void Add(T val)
Definition: healpix_map.h:254
tsize size() const
Healpix_Map(int nside, Healpix_Ordering_Scheme scheme, const nside_dummy)
Definition: healpix_map.h:56
double rms() const
Definition: healpix_map.h:271
const arr< T > & Map() const
Definition: healpix_map.h:227
void Import_upgrade(const Healpix_Map< T > &orig)
Definition: healpix_map.h:119
void fill(const T &val)
Definition: healpix_map.h:90
Healpix_Ordering_Scheme scheme_
Definition: healpix_base.h:56

Generated on Mon Dec 10 2018 10:24:22 for Healpix C++