Makefile.am

ACLOCAL_AMFLAGS = -I m4

lib_LTLIBRARIES = libhealpix_cxx.la

src_c_utils = \
  c_utils/c_utils.c \
  c_utils/c_utils.h \
  c_utils/walltime_c.c \
  c_utils/walltime_c.h \
  c_utils/trig_utils.c \
  c_utils/trig_utils.h
src_libpocketfft = \
  pocketfft/pocketfft.c \
  pocketfft/pocketfft.h
src_libsharp = \
  libsharp/sharp.c \
  libsharp/sharp.h \
  libsharp/sharp_cxx.h \
  libsharp/sharp_almhelpers.c \
  libsharp/sharp_almhelpers.h \
  libsharp/sharp_complex_hacks.h \
  libsharp/sharp_core.c \
  libsharp/sharp_core_avx.c \
  libsharp/sharp_core.h \
  libsharp/sharp_geomhelpers.c \
  libsharp/sharp_geomhelpers.h \
  libsharp/sharp_internal.h \
  libsharp/sharp_lowlevel.h \
  libsharp/sharp_vecsupport.h \
  libsharp/sharp_vecutil.h \
  libsharp/sharp_ylmgen_c.c \
  libsharp/sharp_ylmgen_c.h \
  libsharp/sharp_legendre_roots.c \
  libsharp/sharp_legendre_roots.h
src_cxxsupport = \
  cxxsupport/announce.cc \
  cxxsupport/geom_utils.cc \
  cxxsupport/string_utils.cc \
  cxxsupport/ls_image.cc \
  cxxsupport/paramfile.cc \
  cxxsupport/pointing.cc \
  cxxsupport/rotmatrix.cc \
  cxxsupport/trafos.cc \
  cxxsupport/walltimer.cc \
  cxxsupport/wigner.cc \
  cxxsupport/error_handling.cc \
  cxxsupport/alloc_utils.h \
  cxxsupport/announce.h \
  cxxsupport/arr.h \
  cxxsupport/bstream.h \
  cxxsupport/datatypes.h \
  cxxsupport/error_handling.h \
  cxxsupport/fftpack_support.h \
  cxxsupport/geom_utils.h \
  cxxsupport/levels_facilities.h \
  cxxsupport/ls_image.h \
  cxxsupport/lsconstants.h \
  cxxsupport/math_utils.h \
  cxxsupport/openmp_support.h \
  cxxsupport/paramfile.h \
  cxxsupport/planck_rng.h \
  cxxsupport/pointing.h \
  cxxsupport/rangeset.h \
  cxxsupport/crangeset.h \
  cxxsupport/rotmatrix.h \
  cxxsupport/safe_cast.h \
  cxxsupport/safe_ptr.h \
  cxxsupport/share_utils.h \
  cxxsupport/sse_utils_cxx.h \
  cxxsupport/string_utils.h \
  cxxsupport/trafos.h \
  cxxsupport/vec3.h \
  cxxsupport/walltimer.h \
  cxxsupport/wigner.h \
  cxxsupport/xcomplex.h \
  cxxsupport/compress_utils.h \
  cxxsupport/colour.h \
  cxxsupport/linear_map.h \
  cxxsupport/sort_utils.h
src_cxxsupport_fits = \
  cxxsupport/fitshandle.cc \
  cxxsupport/fitshandle.h
src_healpix_cxx= \
  Healpix_cxx/alm.cc \
  Healpix_cxx/alm.h \
  Healpix_cxx/alm_healpix_tools.cc \
  Healpix_cxx/alm_healpix_tools.h \
  Healpix_cxx/alm_powspec_tools.cc \
  Healpix_cxx/alm_powspec_tools.h \
  Healpix_cxx/healpix_tables.cc \
  Healpix_cxx/healpix_tables.h \
  Healpix_cxx/healpix_base.cc \
  Healpix_cxx/healpix_base.h \
  Healpix_cxx/healpix_map.cc \
  Healpix_cxx/healpix_map.h \
  Healpix_cxx/powspec.cc \
  Healpix_cxx/powspec.h \
  Healpix_cxx/mask_tools.cc \
  Healpix_cxx/mask_tools.h \
  Healpix_cxx/moc_query.h \
  Healpix_cxx/moc_query.cc \
  Healpix_cxx/weight_utils.h \
  Healpix_cxx/weight_utils.cc
src_healpix_cxx_fits= \
  Healpix_cxx/healpix_data_io.cc \
  Healpix_cxx/healpix_data_io.h \
  Healpix_cxx/healpix_map_fitsio.cc \
  Healpix_cxx/healpix_map_fitsio.h \
  Healpix_cxx/alm_fitsio.cc \
  Healpix_cxx/alm_fitsio.h \
  Healpix_cxx/powspec_fitsio.cc \
  Healpix_cxx/powspec_fitsio.h \
  Healpix_cxx/moc_fitsio.h \
  Healpix_cxx/moc_fitsio.cc

EXTRA_DIST = \
  libsharp/sharp_core_inc0.c \
  libsharp/sharp_core_inc.c \
  libsharp/sharp_core_inc2.c \
  libsharp/sharp_core_inchelper.c \
  cxxsupport/font_data.inc \
  pocketfft/LICENSE.md

include_HEADERS = \
  libsharp/sharp.h \
  libsharp/sharp_lowlevel.h \
  libsharp/sharp_geomhelpers.c \
  libsharp/sharp_geomhelpers.h \
  cxxsupport/arr.h \
  cxxsupport/xcomplex.h \
  cxxsupport/datatypes.h \
  cxxsupport/pointing.h \
  cxxsupport/geom_utils.h \
  cxxsupport/rangeset.h \
  cxxsupport/fitshandle.h \
  cxxsupport/alloc_utils.h \
  cxxsupport/math_utils.h \
  cxxsupport/error_handling.h \
  cxxsupport/lsconstants.h \
  cxxsupport/safe_cast.h \
  cxxsupport/vec3.h \
  cxxsupport/wigner.h \
  cxxsupport/rotmatrix.h \
  cxxsupport/sse_utils_cxx.h \
  cxxsupport/string_utils.h \
  Healpix_cxx/alm.h \
  Healpix_cxx/alm_fitsio.h \
  Healpix_cxx/alm_healpix_tools.h \
  Healpix_cxx/alm_powspec_tools.h \
  Healpix_cxx/healpix_base.h \
  Healpix_cxx/healpix_data_io.h \
  Healpix_cxx/healpix_map.h \
  Healpix_cxx/healpix_map_fitsio.h \
  Healpix_cxx/healpix_tables.h \
  Healpix_cxx/powspec.h \
  Healpix_cxx/powspec_fitsio.h \
  Healpix_cxx/mask_tools.h \
  Healpix_cxx/moc.h \
  Healpix_cxx/moc_fitsio.h \
  Healpix_cxx/moc_query.h \
  Healpix_cxx/weight_utils.h

libhealpix_cxx_la_SOURCES = $(src_c_utils) $(src_libpocketfft) $(src_libsharp) $(src_cxxsupport) $(src_cxxsupport_fits) $(src_healpix_cxx) $(src_healpix_cxx_fits) $(src_healpix_cxx_mod)
libhealpix_cxx_la_LIBADD = $(CFITSIO_LIBS)
libhealpix_cxx_la_LDFLAGS = -version-info 2:1:0

bin_PROGRAMS = alm2map_cxx anafast_cxx calc_powspec hotspots_cxx map2tga \
  median_filter_cxx mult_alm rotalm_cxx smoothing_cxx syn_alm_cxx udgrade_cxx \
  udgrade_harmonic_cxx compute_weights needlet_tool
check_PROGRAMS = hpxtest
alm2map_cxx_SOURCES = Healpix_cxx/alm2map_cxx.cc Healpix_cxx/alm2map_cxx_module.cc
alm2map_cxx_LDADD = libhealpix_cxx.la
anafast_cxx_SOURCES = Healpix_cxx/anafast_cxx.cc Healpix_cxx/anafast_cxx_module.cc
anafast_cxx_LDADD = libhealpix_cxx.la
calc_powspec_SOURCES = Healpix_cxx/calc_powspec.cc Healpix_cxx/calc_powspec_module.cc
calc_powspec_LDADD = libhealpix_cxx.la
hotspots_cxx_SOURCES = Healpix_cxx/hotspots_cxx.cc Healpix_cxx/hotspots_cxx_module.cc
hotspots_cxx_LDADD = libhealpix_cxx.la
map2tga_SOURCES = Healpix_cxx/map2tga.cc Healpix_cxx/map2tga_module.cc
map2tga_LDADD = libhealpix_cxx.la
median_filter_cxx_SOURCES = Healpix_cxx/median_filter_cxx.cc Healpix_cxx/median_filter_cxx_module.cc
median_filter_cxx_LDADD = libhealpix_cxx.la
mult_alm_SOURCES = Healpix_cxx/mult_alm.cc Healpix_cxx/mult_alm_module.cc
mult_alm_LDADD = libhealpix_cxx.la
rotalm_cxx_SOURCES = Healpix_cxx/rotalm_cxx.cc
rotalm_cxx_LDADD = libhealpix_cxx.la
smoothing_cxx_SOURCES = Healpix_cxx/smoothing_cxx.cc Healpix_cxx/smoothing_cxx_module.cc
smoothing_cxx_LDADD = libhealpix_cxx.la
syn_alm_cxx_SOURCES = Healpix_cxx/syn_alm_cxx.cc Healpix_cxx/syn_alm_cxx_module.cc
syn_alm_cxx_LDADD = libhealpix_cxx.la
udgrade_cxx_SOURCES = Healpix_cxx/udgrade_cxx.cc Healpix_cxx/udgrade_cxx_module.cc
udgrade_cxx_LDADD = libhealpix_cxx.la
udgrade_harmonic_cxx_SOURCES = Healpix_cxx/udgrade_harmonic_cxx.cc Healpix_cxx/udgrade_harmonic_cxx_module.cc
udgrade_harmonic_cxx_LDADD = libhealpix_cxx.la
compute_weights_SOURCES = Healpix_cxx/compute_weights.cc Healpix_cxx/compute_weights_module.cc
compute_weights_LDADD = libhealpix_cxx.la
needlet_tool_SOURCES = Healpix_cxx/needlet_tool.cc Healpix_cxx/needlet_tool_module.cc
needlet_tool_LDADD = libhealpix_cxx.la
hpxtest_SOURCES=Healpix_cxx/hpxtest.cc
hpxtest_LDADD = libhealpix_cxx.la

TESTS=hpxtest

AM_CPPFLAGS = $(CFITSIO_CFLAGS) -I$(top_srcdir)/c_utils -I$(top_srcdir) -I$(top_srcdir)/libsharp -I$(top_srcdir)/cxxsupport

pkgconfigdir = $(libdir)/pkgconfig
nodist_pkgconfig_DATA = @PACKAGE_NAME@.pc

DISTCLEANFILES=@PACKAGE_NAME@.pc @PACKAGE_NAME@.pc.in  @PACKAGE_NAME@-uninstalled.pc  @PACKAGE_NAME@-uninstalled.sh

mask_tools.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).
 */

/*! \file mask_tools.cxx
 *  Copyright (C) 2003-2019 Max-Planck-Society
 *  \author Xavier Garrido
 */

#include "mask_tools.h"
#include "lsconstants.h"

#include <numeric>
#include <algorithm>

void maskborder(const Healpix_Map<double> & mask_in,
                Healpix_Map<double> & mask_out,
                size_t & nbordpix)
{
  double minmask = 0;
  double maxmask = 0;
  mask_in.minmax(minmask, maxmask);

  nbordpix = 0;
  planck_assert(minmask == PIXBAD, "No mask region found in map!");
  planck_assert(maxmask == PIXGOOD, "No physical region found in map!");

  for (int p = 0; p < mask_out.Npix(); p++) {
    if (mask_out[p] != PIXBAD) continue;
    fix_arr<int, 8> nb;
    mask_out.neighbors(p, nb);
    for (size_t i = 0; i < 8; i++) {
      if (nb[i] != -1 && mask_out[nb[i]] == PIXGOOD) {
        nbordpix++;
        mask_out[p] = DEFBTAG;
        break;
      }
    } // end over neighbors
  } // end over pixels of the map
}

void dist2holes(const Healpix_Map<double> & mask,
                Healpix_Map<double> & distances)
{
  // Work in NEST pixelisation (faster)
  const Healpix_Ordering_Scheme orig_scheme = distances.Scheme();
  if (orig_scheme == RING) {
    distances.swap_scheme();
  }
  Healpix_Map<double> mask1(mask);
  if (mask1.Scheme() == RING) {
    mask1.swap_scheme();
  }

  size_t nbordpix = 0;
  maskborder(mask, mask1, nbordpix);

  // Treat particular case of no border
  if (nbordpix == 0) {
    distances.fill(0.0);
    return;
  }

  // Find 3D location of those border pixels
  std::vector<vec3> bordpos;
  std::vector<int> bordpix;
  for (size_t p = 0; p < size_t(mask1.Npix()); p++) {
    if (mask1[p] == PIXGOOD) continue;
    distances[p] = 0.0;
    if (mask1[p] == DEFBTAG) {
      bordpos.push_back(mask1.pix2vec(p));
      bordpix.push_back(p);
    }
  }

  // Define low resolution grid
  const size_t nslow = 64;
  const size_t nside = mask1.Nside();
  const size_t nslow_in = std::min(nslow, nside);
  Healpix_Map<double> low_res_map(nslow_in, NEST, SET_NSIDE);
  const size_t nplow = low_res_map.Npix();
  const size_t iratio = std::pow(size_t(nside/nslow_in), 2); // ratio of number of pixels

  // Find low resolution pixels containing border pixels (=non_empty)
  std::vector<int> nb_in_lp(nplow, 0);
  for (size_t p = 0; p < nbordpix; p++) {
    const size_t q = bordpix[p] / iratio;
    nb_in_lp[q]++;// number of border pixels in each low-res pixel
  }
  const size_t sum_in_lp = std::accumulate(nb_in_lp.begin(), nb_in_lp.end(), 0);
  if (sum_in_lp != nbordpix) {
    std::ostringstream msg;
    msg << "Error with low res. pixels (" << sum_in_lp << "!=" << nbordpix << ")";
    planck_assert(true, msg.str());
  }
  const size_t max_in_lp = *std::max_element(nb_in_lp.begin(), nb_in_lp.end());
  if (max_in_lp > iratio) {
    std::ostringstream msg;
    msg << "Error with low res. pixels (" << max_in_lp << ">" << iratio << ")";
    planck_assert(true, msg.str());
  }

  std::vector<std::pair<size_t, size_t> > non_empty_low;
  for (size_t q = 0; q < nplow; q++) {
    if (nb_in_lp[q] > 0) {
      non_empty_low.push_back(std::make_pair(q, nb_in_lp[q]));
    }
  }
  const size_t nploweff = std::count_if(nb_in_lp.begin(), nb_in_lp.end(),
                                        [](size_t i){return i > 0;});
  std::vector<size_t> start_list(nploweff+1, 0);
  for (size_t q = 1; q < nploweff+1; q++) {
    start_list[q] = start_list[q-1] + non_empty_low[q-1].second;
  }
  if (start_list[nploweff] != nbordpix) {
    std::ostringstream msg;
    msg << "Error with low res. pixels (" << start_list[nploweff] << "!=" << nbordpix << ")";
    planck_assert(true, msg.str());
  }
  nb_in_lp.clear();

  // Build location of center for non-empty low resolution "border" pixels
  std::vector<vec3> centlow(nploweff);
  for (size_t qq = 0; qq < nploweff; qq++) {
    const size_t q = non_empty_low[qq].first;
    centlow[qq] = low_res_map.pix2vec(q);
  }

  // Loop on low-resolution pixels
  const double fudge = std::min(4.276/double(nslow_in), pi); // 2.0 x diagonal of longest pixel, at large Nside
  const double sfudge = std::sin(fudge);
  const double cfudge = std::cos(fudge);

  std::vector<double> drange(nploweff);
  std::vector<size_t> p2test(nploweff);
  std::vector<vec3> bordpospack(nbordpix);
  std::vector<vec3> centpix(iratio);
  std::vector<double> dd(iratio);
  size_t n_active = 0;
  size_t n_invalid = 0;
  for (size_t q1 = 0; q1 < nplow; q1++) {
    // Keep only those containing valid small pixel
    const arr<double> & map = mask1.Map();
    const auto start = map.begin()+q1*iratio;
    const auto stop = map.begin()+(q1+1)*iratio;
    if (std::find(start, stop, PIXGOOD) == stop) continue;
    // In C++ 17
    // if (! std::any_of(map.begin()+q1*iratio, map.begin()+(q1+1)*iratio,
    //                   [](int i) {return i == PIXGOOD;})) continue;

    n_active++;
    // Center location of current low-res pixel
    const vec3 & vv = low_res_map.pix2vec(q1);
#pragma omp parallel for schedule (static, 4)
    for (size_t qq = 0; qq < nploweff; qq++) {
      drange[qq] = dotprod(centlow[qq], vv);
    }

    // Min distance (=max scalar product)
    const double distmin = std::min(1.0, *std::max_element(drange.begin(), drange.end()));

    // Find out list of low-res "border" pixels close enough to be checked-out at higher-res
    // radius = smallest distance found above + fudge
    double threshold = -1.0;
    if ((distmin + cfudge) > 0.0) {
      // New radius is less than Pi
      threshold = distmin*cfudge - sfudge*sqrt((1.0-distmin)*(1.0+distmin));
    }
    size_t n2test = 0;
    for (size_t qq = 0; qq < nploweff; qq++) {
      if (drange[qq] >= threshold) {
        // For scalar product
        p2test[n2test] = qq;
        n2test++;
      }
    }

    // Store location of border pixels contained in relevant low-res pixels
    size_t pstart = 0;
    for (size_t q = 0; q < n2test; q++) {
      const size_t qq = p2test[q];
      const size_t i1 = start_list[qq];
      const size_t i2 = start_list[qq+1];
      const size_t n_in_list = i2 -i1;
      std::copy(bordpos.begin()+i1, bordpos.begin()+i2, bordpospack.begin()+pstart);
      pstart += n_in_list;
    }
    const size_t nb = pstart; // Number of border pixels close enough to be considered
    // mean_nbord = mean_nbord + nb

    // Center location for small pixels
    size_t pp = 0;
    for (size_t p = q1*iratio; p < (q1+1)*iratio; p++) {
      if (mask1[p] != PIXGOOD) continue;
      centpix[pp] = mask1.pix2vec(p);
      pp++;
    }
    const size_t np = pp; // Number of valid small pixels in current low-res pixels

    // Scalar product: small pixels * border pixels
#pragma omp parallel for schedule (guided, 4)
    for (size_t i = 0; i < np; i++) {
      double dmin = 1e6;
      for (size_t j = 0; j < nb; j++) {
        const double d = (bordpospack[j] - centpix[i]).SquaredLength();
        if (d < dmin) dmin = d;
      }
      dd[i] = dmin;
    }

    // norm2 -> angular distance (radians)
    pp = 0;
    for (size_t p = q1*iratio; p < (q1+1)*iratio; p++) {
      if (mask1[p] != PIXGOOD) continue;
      if (dd[pp] < 0.0 || dd[pp] > 4.0) n_invalid++;
      distances[p] = 2.0 * std::asin(std::sqrt(dd[pp]) * 0.5);
      pp++;
    }
  }

  // Get back to initial scheme
  if (orig_scheme == RING) {
    distances.swap_scheme();
  }
}

mask_tools.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 mask_tools.h
 *  Copyright (C) 2003-2019 Max-Planck-Society
 *  \author Xavier Garrido
 */

#ifndef MASK_TOOLS_H
#define MASK_TOOLS_H

#include "healpix_map.h"

// Enum for pixel status within a mask map
enum PixelStatus
  {
   PIXBAD = 0,  // Invalid pixel
   PIXGOOD = 1, // Valid pixel
   DEFBTAG = 2  // Default border value
  };

// Mask pixels in the border (adapted from Fortran code by E. Hivon)
// mask_in IN, Healpix map, in NESTED scheme
//    on input: mask of good pixels (mask=1) and bad pixels (mask=0)
//
// mask_out OUT, Healpix map, in NESTED scheme
//    on output: inner border of bad regions takesmask=DEFBTAG=2
//
// nbordpix OUT
//    number of pixels in border
void maskborder(const Healpix_Map<double> & mask_in, Healpix_Map<double> & mask_out, size_t & nbordpix);

// For each valid pixel, returns distance (in radians) from pixel center to
// center of closest invalid pixel.
// The input mask must be NESTED ordered
//
// Algorithm:
//  To begin with, the inner border (=invalid pixels in touch with valid ones)
//  of each invalid region is identified
// Then a small Nside is used to select the closest border pixels, and finally
// the distance of each valid pixel to the selected border pixels is computed
// to find its minimum.
//
// Note: final distances are computed with  2 ASIN (|v1-v2|/2) which is more
// accurate than the usual ACOS(v1.v2) at small separations.
//
// Adapted from Fortran code by E. Hivon
void dist2holes(const Healpix_Map<double> & mask, Healpix_Map<double> & distances);

#endif