add Healpix_Map and tests

Healpix_cxx/healpix_base.h

@@ -86,7 +86,8 @@ template<typename I> class T_Healpix_Base: public Healpix_Tables
 
     I nest_peano_helper (I pix, int dir) const;
 
-    using swapfunc = I(I) const;
+    typedef I (T_Healpix_Base::*swapfunc)(I pix) const;
+//    using swapfunc = I(I) const;
 
   public:
     enum {order_max=Orderhelper__<I>::omax};

Healpix_cxx/healpix_map.cc

+/*
+ *  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).
+ */
+
+/*
+ *  Copyright (C) 2003-2020 Max-Planck-Society
+ *  Author: Martin Reinecke
+ */
+
+#include "healpix_map.h"
+
+using namespace std;
+
+namespace mr {
+
+namespace healpix {
+
+// template<typename T> void Healpix_Map<T>::Import_degrade
+//   (const Healpix_Map<T> &orig, bool pessimistic)
+//   {
+//   MR_assert(nside_<orig.nside_,"Import_degrade: this is no degrade");
+//   int fact = orig.nside_/nside_;
+//   MR_assert (orig.nside_==nside_*fact,
+//     "the larger Nside must be a multiple of the smaller one");
+//
+//   int minhits = pessimistic ? fact*fact : 1;
+// #pragma omp parallel
+// {
+//   int m;
+// #pragma omp for schedule (static)
+//   for (m=0; m<npix_; ++m)
+//     {
+//     int x,y,f;
+//     pix2xyf(m,x,y,f);
+//     int hits = 0;
+//     kahan_adder<double> adder;
+//     for (int j=fact*y; j<fact*(y+1); ++j)
+//       for (int i=fact*x; i<fact*(x+1); ++i)
+//         {
+//         int opix = orig.xyf2pix(i,j,f);
+//         if (!approx<double>(orig.map[opix],Healpix_undef))
+//           {
+//           ++hits;
+//           adder.add(orig.map[opix]);
+//           }
+//         }
+//     map[m] = T((hits<minhits) ? Healpix_undef : adder.result()/hits);
+//     }
+// }
+//   }
+//
+// template void Healpix_Map<float>::Import_degrade
+//   (const Healpix_Map<float> &orig, bool pessimistic);
+// template void Healpix_Map<double>::Import_degrade
+//   (const Healpix_Map<double> &orig, bool pessimistic);
+
+template<typename T> void Healpix_Map<T>::minmax (T &Min, T &Max) const
+  {
+  Min = T(1e30); Max = T(-1e30);
+  for (int m=0; m<npix_; ++m)
+    {
+    T val = map[m];
+    if (!approx<double>(val,Healpix_undef))
+      {
+      if (val>Max) Max=val;
+      if (val<Min) Min=val;
+      }
+    }
+  }
+
+template void Healpix_Map<float>::minmax (float &Min, float &Max) const;
+template void Healpix_Map<double>::minmax (double &Min, double &Max) const;
+
+namespace {
+
+template<typename Iterator> typename iterator_traits<Iterator>::value_type
+  mymedian(Iterator first, Iterator last)
+  {
+  Iterator mid = first+(last-first-1)/2;
+  nth_element(first,mid,last);
+  if ((last-first)&1) return *mid;
+  return typename iterator_traits<Iterator>::value_type
+    (0.5*((*mid)+(*min_element(mid+1,last))));
+  }
+
+} // unnamed namespace
+
+template<typename T> Healpix_Map<T> Healpix_Map<T>::median (double rad) const
+  {
+  Healpix_Map<T> out(Nside(), Scheme(), SET_NSIDE);
+
+  rangeset<int> pixset;
+  vector<T> list;
+  for (int m=0; m<Npix(); ++m)
+    {
+    query_disc(pix2ang(m),rad,pixset);
+    list.resize(pixset.nval());
+    size_t cnt=0;
+    for (size_t j=0; j<pixset.nranges(); ++j)
+      for (int i=pixset.ivbegin(j); i<pixset.ivend(j); ++i)
+        if (!approx(map[i], T(Healpix_undef)))
+          list[cnt++] = map[i];
+    out[m] = (cnt>0) ? mymedian(list.begin(),list.begin()+cnt)
+                     : Healpix_undef;
+    }
+  return out;
+  }
+
+template Healpix_Map<float> Healpix_Map<float>::median (double rad) const;
+template Healpix_Map<double> Healpix_Map<double>::median (double rad) const;
+
+}}

Healpix_cxx/healpix_map.h

+/*
+ *  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>
+#include "mr_util/error_handling.h"
+#include "mr_util/threading.h"
+#include "Healpix_cxx/healpix_base.h"
+
+namespace mr {
+
+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. */
+//     double average() const
+//       {
+//       kahan_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;
+//       }
+
+    /*! 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

Makefile.am

@@ -38,7 +38,9 @@ libmrutil_la_SOURCES = \
   Healpix_cxx/healpix_tables.h \
   Healpix_cxx/healpix_tables.cc \
   Healpix_cxx/healpix_base.h \
-  Healpix_cxx/healpix_base.cc
+  Healpix_cxx/healpix_base.cc \
+  Healpix_cxx/healpix_map.h \
+  Healpix_cxx/healpix_map.cc
 
 # format is "current:revision:age"
 # any change: increase revision
@@ -77,11 +79,13 @@ nobase_include_HEADERS = \
 
 EXTRA_DIST = test/test_libsharp.sh test/test_space_filling.sh
 
-check_PROGRAMS = sharp2_testsuite space_filling_test
+check_PROGRAMS = sharp2_testsuite space_filling_test hpxtest
 sharp2_testsuite_SOURCES = test/sharp2_testsuite.cc
 sharp2_testsuite_LDADD = libmrutil.la
 space_filling_test_SOURCES = test/space_filling_test.cc
 space_filling_test_LDADD = libmrutil.la
+hpxtest_SOURCES = test/hpxtest.cc
+hpxtest_LDADD = libmrutil.la
 
 TESTS = test/test_libsharp.sh test/test_space_filling.sh