healpix_map.h 10.5 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
/*
 *  This file is part of Healpix_cxx.
 *
 *  Healpix_cxx is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  Healpix_cxx is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with Healpix_cxx; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 *  For more information about HEALPix, see http://healpix.sourceforge.net
 */

/*
 *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
 *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
 *  (DLR).
 */

/*! \file healpix_map.h
 *  Copyright (C) 2003-2020 Max-Planck-Society
 *  \author Martin Reinecke
 */

#ifndef HEALPIX_MAP_H
#define HEALPIX_MAP_H

#include <vector>
#include <algorithm>
Martin Reinecke's avatar
Martin Reinecke committed
37
38
#include "ducc0/error_handling.h"
#include "ducc0/threading.h"
Martin Reinecke's avatar
Martin Reinecke committed
39
40
#include "Healpix_cxx/healpix_base.h"

Martin Reinecke's avatar
Martin Reinecke committed
41
namespace ducc0 {
Martin Reinecke's avatar
Martin Reinecke committed
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242

namespace healpix {

//! Healpix value representing "undefined"
const double Healpix_undef=-1.6375e30;

/*! A HEALPix map of a given datatype */
template<typename T> class Healpix_Map: public Healpix_Base
  {
  private:
    std::vector<T> map;

  public:
    /*! Constructs an unallocated map. */
    Healpix_Map () {}
    /*! Constructs a map with a given \a order and the ordering
        scheme \a scheme. */
    Healpix_Map (int order, Ordering_Scheme scheme)
      : Healpix_Base (order, scheme), map(npix_) {}
    /*! Constructs a map with a given \a nside and the ordering
        scheme \a scheme. */
    Healpix_Map (int nside, Ordering_Scheme scheme, const nside_dummy)
      : Healpix_Base (nside, scheme, SET_NSIDE), map(npix_) {}
    /*! Constructs a map from the contents of \a data and sets the ordering
        scheme to \a Scheme. The size of \a data must be a valid HEALPix
        map size. */
    Healpix_Map (const std::vector<T> &data, Ordering_Scheme scheme)
      : Healpix_Base (npix2nside(data.size()), scheme, SET_NSIDE), map(data) {}

    /*! Deletes the old map, creates a map from the contents of \a data and
        sets the ordering scheme to \a scheme. The size of \a data must be a
        valid HEALPix map size.
        \note On exit, \a data is zero-sized! */
    void Set (std::vector<T> &data, Ordering_Scheme scheme)
      {
      Healpix_Base::SetNside(npix2nside (data.size()), scheme);
      map.transfer(data);
      }

    /*! Deletes the old map and creates a new map  with a given \a order
        and the ordering scheme \a scheme. */
    void Set (int order, Ordering_Scheme scheme)
      {
      Healpix_Base::Set(order, scheme);
      map.alloc(npix_);
      }
    /*! Deletes the old map and creates a new map  with a given \a nside
        and the ordering scheme \a scheme. */
    void SetNside (int nside, Ordering_Scheme scheme)
      {
      Healpix_Base::SetNside(nside, scheme);
      map.alloc(npix_);
      }

    /*! Fills the map with \a val. */
    void fill (const T &val)
      { std::fill(map.begin(),map.end(),val); }

    /*! Imports the map \a orig into the current map, adjusting the
        ordering scheme. \a orig must have the same resolution as the
        current map. */
    void Import_nograde (const Healpix_Map<T> &orig)
      {
      MR_assert (nside_==orig.nside_,
        "Import_nograde: maps have different nside");
      if (orig.scheme_ == scheme_)
        for (int m=0; m<npix_; ++m) map[m] = orig.map[m];
      else
        {
        swapfunc swapper = (scheme_ == NEST) ?
          &Healpix_Base::ring2nest : &Healpix_Base::nest2ring;
#pragma omp parallel for schedule (dynamic,5000)
        for (int m=0; m<npix_; ++m) map[(this->*swapper)(m)] = orig.map[m];
        }
      }

    /*! Imports the map \a orig into the current map, adjusting the
        ordering scheme and the map resolution. \a orig must have lower
        resolution than the current map, and \a this->Nside() must be an
        integer multiple of \a orig.Nside(). */
    void Import_upgrade (const Healpix_Map<T> &orig)
      {
      MR_assert(nside_>orig.nside_,"Import_upgrade: this is no upgrade");
      int fact = nside_/orig.nside_;
      MR_assert (nside_==orig.nside_*fact,
        "the larger Nside must be a multiple of the smaller one");

#pragma omp parallel for schedule (dynamic,5000)
      for (int m=0; m<orig.npix_; ++m)
        {
        int x,y,f;
        orig.pix2xyf(m,x,y,f);
        for (int j=fact*y; j<fact*(y+1); ++j)
          for (int i=fact*x; i<fact*(x+1); ++i)
            {
            int mypix = xyf2pix(i,j,f);
            map[mypix] = orig.map[m];
            }
        }
      }

    /*! Imports the map \a orig into the current map, adjusting the
        ordering scheme and the map resolution. \a orig must have higher
        resolution than the current map, and \a orig.Nside() must be an
        integer multiple of \a this->Nside().
        \a pessimistic determines whether or not
        pixels are set to \a Healpix_undef when not all of the corresponding
        high-resolution pixels are defined.

        This method is instantiated for \a float and \a double only. */
    void Import_degrade (const Healpix_Map<T> &orig, bool pessimistic=false);

    /*! Imports the map \a orig into the current map, adjusting the
        ordering scheme and the map resolution if necessary.
        When downgrading, \a pessimistic determines whether or not
        pixels are set to \a Healpix_undef when not all of the corresponding
        high-resolution pixels are defined.

        This method is instantiated for \a float and \a double only. */
    void Import (const Healpix_Map<T> &orig, bool pessimistic=false)
      {
      if (orig.nside_ == nside_) // no up/degrading
        Import_nograde(orig);
      else if (orig.nside_ < nside_) // upgrading
        Import_upgrade(orig);
      else
        Import_degrade(orig,pessimistic);
      }

    /*! Returns a constant reference to the pixel with the number \a pix. */
    const T &operator[] (int pix) const { return map[pix]; }
    /*! Returns a reference to the pixel with the number \a pix. */
    T &operator[] (int pix) { return map[pix]; }

    /*! Swaps the map ordering from RING to NEST and vice versa.
        This is done in-place (i.e. with negligible space overhead). */
    void swap_scheme()
      {
      swapfunc swapper = (scheme_ == NEST) ?
        &Healpix_Base::ring2nest : &Healpix_Base::nest2ring;

      std::vector<int> cycle=swap_cycles();

#pragma omp parallel for schedule(dynamic,1)
      for (size_t m=0; m<cycle.size(); ++m)
        {
        int istart = cycle[m];

        T pixbuf = map[istart];
        int iold = istart, inew = (this->*swapper)(istart);
        while (inew != istart)
          {
          map[iold] = map[inew];
          iold = inew;
          inew = (this->*swapper)(inew);
          }
        map[iold] = pixbuf;
        }
      scheme_ = (scheme_==RING) ? NEST : RING;
      }

    /*! performs the actual interpolation using \a pix and \a wgt. */
    T interpolation (const std::array<int,4> &pix,
      const std::array<double,4> &wgt) const
      {
      double wtot=0;
      T res=T(0);
      for (size_t i=0; i<4; ++i)
        {
        T val=map[pix[i]];
        if (!approx<double>(val,Healpix_undef))
          { res+=T(val*wgt[i]); wtot+=wgt[i]; }
        }
      return (wtot==0.) ? T(Healpix_undef) : T(res/wtot);
      }
    /*! Returns the interpolated map value at \a ptg */
    T interpolated_value (const pointing &ptg) const
      {
      std::array<int,4> pix;
      std::array<double,4> wgt;
      get_interpol (ptg, pix, wgt);
      return interpolation (pix, wgt);
      }

    /*! Returns a constant reference to the map data. */
    const std::vector<T> &Map() const { return map; }

    /*! Returns the minimum and maximum value of the map in
        \a Min and \a Max.

        This method is instantiated for \a float and \a double only. */
    void minmax (T &Min, T &Max) const;

    /*! Swaps the contents of two Healpix_Map objects. */
    void swap (Healpix_Map &other)
      {
      Healpix_Base::swap(other);
      map.swap(other.map);
      }

    /*! Returns the average of all defined map pixels. */
Martin Reinecke's avatar
Martin Reinecke committed
243
244
245
246
247
248
249
250
251
    double average() const
      {
      tree_adder<double> adder;
      int pix=0;
      for (int m=0; m<npix_; ++m)
        if (!approx<double>(map[m],Healpix_undef))
          { ++pix; adder.add(map[m]); }
      return (pix>0) ? adder.result()/pix : Healpix_undef;
      }
Martin Reinecke's avatar
Martin Reinecke committed
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321

    /*! Adds \a val to all defined map pixels. */
    void Add (T val)
      {
      for (int m=0; m<npix_; ++m)
        if (!approx<double>(map[m],Healpix_undef))
          { map[m]+=val; }
      }

    /*! Multiplies all defined map pixels by \a val. */
    void Scale (T val)
      {
      for (int m=0; m<npix_; ++m)
        if (!approx<double>(map[m],Healpix_undef))
          { map[m]*=val; }
      }

    /*! Returns the root mean square of the map, not counting undefined
        pixels. */
    double rms() const
      {
      using namespace std;

      double result=0;
      int pix=0;
      for (int m=0; m<npix_; ++m)
        if (!approx<double>(map[m],Healpix_undef))
          { ++pix; result+=map[m]*map[m]; }
      return (pix>0) ? sqrt(result/pix) : Healpix_undef;
      }
    /*! Returns the maximum absolute value in the map, ignoring undefined
        pixels. */
    T absmax() const
      {
      using namespace std;

      T result=0;
      for (int m=0; m<npix_; ++m)
        if (!approx<double>(map[m],Healpix_undef))
          { result = max(result,abs(map[m])); }
      return result;
      }
    /*! Returns \a true, if no pixel has the value \a Healpix_undef,
        else \a false. */
    bool fullyDefined() const
      {
      for (int m=0; m<npix_; ++m)
        if (approx<double>(map[m],Healpix_undef))
          return false;
      return true;
      }
    /*! Sets all pixels with the value \a Healpix_undef to 0, and returns
        the number of modified pixels. */
    size_t replaceUndefWith0()
      {
      size_t res=0;
      for (int m=0; m<npix_; ++m)
        if (approx<double>(map[m],Healpix_undef))
          { map[m]=0.; ++res; }
      return res;
      }

    /*! Returns a map that contains at each pixel the median value of all
        pixels within the radius \a rad (in radians) around the pixel center. */
    Healpix_Map median(double rad) const;
  };

}}

#endif