Skip to content

GitLab

Commit 714edde6 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

add more Healpix functionality

parent 39552dd8
Expand all Hide whitespace changes
Inline Side-by-side
Healpix_cxx/moc.h 0 → 100644
View file @ 714edde6
/*
* 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 moc.h
* Copyright (C) 2014-2020 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef HEALPIX_MOC_H
#define HEALPIX_MOC_H
#include "Healpix_cxx/healpix_base.h"
#include "mr_util/compress_utils.h"
namespace mr {
namespace healpix {
template<typename I> class Moc
{
public:
enum{maxorder=T_Healpix_Base<I>::order_max};
private:
rangeset<I> rs;
static Moc fromNewRangeSet(const rangeset<I> &rngset)
{
Moc res;
res.rs = rngset;
return res;
}
public:
const rangeset<I> &Rs() { return rs; }
size_t maxOrder() const
{
I combo=0;
for (size_t i=0; i<rs.nranges(); ++i)
combo|=rs.ivbegin(i)|rs.ivend(i);
return maxorder-(trailingZeros(combo)>>1);
}
Moc degradedToOrder (int order, bool keepPartialCells) const
{
int shift=2*(maxorder-order);
I ofs=(I(1)<<shift)-1;
I mask = ~ofs;
I adda = keepPartialCells ? I(0) : ofs,
addb = keepPartialCells ? ofs : I(0);
rangeset<I> rs2;
for (size_t i=0; i<rs.nranges(); ++i)
{
I a=(rs.ivbegin(i)+adda)&mask;
I b=(rs.ivend (i)+addb)&mask;
if (b>a) rs2.append(a,b);
}
return fromNewRangeSet(rs2);
}
void addPixelRange (int order, I p1, I p2)
{
int shift=2*(maxorder-order);
rs.add(p1<<shift,p2<<shift);
}
void appendPixelRange (int order, I p1, I p2)
{
int shift=2*(maxorder-order);
rs.append(p1<<shift,p2<<shift);
}
void appendPixel (int order, I p)
{ appendPixelRange(order,p,p+1); }
/*! Returns a new Moc that contains the union of this Moc and \a other. */
Moc op_or (const Moc &other) const
{ return fromNewRangeSet(rs.op_or(other.rs)); }
/*! Returns a new Moc that contains the intersection of this Moc and
\a other. */
Moc op_and (const Moc &other) const
{ return fromNewRangeSet(rs.op_and(other.rs)); }
Moc op_xor (const Moc &other) const
{ return fromNewRangeSet(rs.op_xor(other.rs)); }
/*! Returns a new Moc that contains all parts of this Moc that are not
contained in \a other. */
Moc op_andnot (const Moc &other) const
{ return fromNewRangeSet(rs.op_andnot(other.rs)); }
/*! Returns the complement of this Moc. */
Moc complement() const
{
rangeset<I> full; full.append(I(0),I(12)*(I(1)<<(2*maxorder)));
return fromNewRangeSet(full.op_andnot(rs));
}
/** Returns \c true if \a other is a subset of this Moc, else \c false. */
bool contains(const Moc &other) const
{ return rs.contains(other.rs); }
/** Returns \c true if the intersection of this Moc and \a other is not
empty. */
bool overlaps(const Moc &other) const
{ return rs.overlaps(other.rs); }
/** Returns a vector containing all HEALPix pixels (in ascending NUNIQ
order) covered by this Moc. The result is well-formed in the sense that
every pixel is given at its lowest possible HEALPix order. */
std::vector<I> toUniq() const
{
std::vector<I> res;
std::vector<std::vector<I> > buf(maxorder+1);
for (size_t i=0; i<rs.nranges(); ++i)
{
I start=rs.ivbegin(i), end=rs.ivend(i);
while(start!=end)
{
unsigned logstep=std::min<int>(maxorder,trailingZeros(start)>>1);
logstep=std::min(logstep,ilog2(end-start)>>1);
buf[maxorder-logstep].push_back(start);
start+=I(1)<<(2*logstep);
}
}
for (int o=0; o<=maxorder; ++o)
{
I ofs=I(1)<<(2*o+2);
int shift=2*(maxorder-o);
for (size_t j=0; j<buf[o].size(); ++j)
res.push_back((buf[o][j]>>shift)+ofs);
}
return res;
}
static Moc fromUniq (const std::vector<I> &vu)
{
rangeset<I> r, rtmp;
int lastorder=0;
int shift=2*maxorder;
for (size_t i=0; i<vu.size(); ++i)
{
int order = ilog2(vu[i]>>2)>>1;
if (order!=lastorder)
{
r=r.op_or(rtmp);
rtmp.clear();
lastorder=order;
shift=2*(maxorder-order);
}
I pix = vu[i]-(I(1)<<(2*order+2));
rtmp.append (pix<<shift,(pix+1)<<shift);
}
r=r.op_or(rtmp);
return fromNewRangeSet(r);
}
static void uniq_nest2peano(std::vector<I> &vu)
{
using namespace std;
if (vu.empty()) return;
size_t start=0;
int order=-1;
T_Healpix_Base<I> base;
I offset=0;
for (size_t j=0; j<vu.size(); ++j)
{
int neworder=ilog2(vu[j]>>2)>>1;
if (neworder>order)
{
sort(vu.begin()+start,vu.begin()+j);
order=neworder;
start=j;
base.Set(order,NEST);
offset=I(1)<<(2*order+2);
}
vu[j]=base.nest2peano(vu[j]-offset)+offset;
}
sort(vu.begin()+start,vu.end());
}
static void uniq_peano2nest(std::vector<I> &vu)
{
using namespace std;
if (vu.empty()) return;
size_t start=0;
int order=-1;
T_Healpix_Base<I> base;
I offset=0;
for (size_t j=0; j<vu.size(); ++j)
{
int neworder=ilog2(vu[j]>>2)>>1;
if (neworder>order)
{
sort(vu.begin()+start,vu.begin()+j);
order=neworder;
start=j;
base.Set(order,NEST);
offset=I(1)<<(2*order+2);
}
vu[j]=base.peano2nest(vu[j]-offset)+offset;
}
sort(vu.begin()+start,vu.end());
}
std::vector<uint8_t> toCompressed() const
{
obitstream obs;
interpol_encode(rs.data().begin(),rs.data().end(),obs);
return obs.state();
}
static Moc fromCompressed(const std::vector<uint8_t> &data)
{
ibitstream ibs(data);
std::vector<I> v;
interpol_decode(v,ibs);
Moc out;
out.rs.setData(v);
return out;
}
bool operator==(const Moc &other) const
{
if (this == &other)
return true;
return rs==other.rs;
}
size_t nranges() const
{ return rs.nranges(); }
I nval() const
{ return rs.nval(); }
};
}}
#endif
Makefile.am
View file @ 714edde6
......@@ -33,6 +33,8 @@ libmrutil_la_SOURCES = \
mr_util/unity_roots.h \
mr_util/useful_macros.h \
mr_util/rangeset.h \
mr_util/crangeset.h \
mr_util/compress_utils.h \
mr_util/geom_utils.h \
mr_util/geom_utils.cc \
Healpix_cxx/healpix_tables.h \
......@@ -40,7 +42,8 @@ libmrutil_la_SOURCES = \
Healpix_cxx/healpix_base.h \
Healpix_cxx/healpix_base.cc \
Healpix_cxx/healpix_map.h \
Healpix_cxx/healpix_map.cc
Healpix_cxx/healpix_map.cc \
Healpix_cxx/moc.h
# format is "current:revision:age"
# any change: increase revision
......
mr_util/compress_utils.h 0 → 100644
View file @ 714edde6
/*
* This file is part of libcxxsupport.
*
* libcxxsupport 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.
*
* libcxxsupport 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 libcxxsupport; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file compress_utils.h
* Support for compression of integer arrays.
*
* Copyright (C) 2013-2020 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef MRUTIL_COMPRESS_UTILS_H
#define MRUTIL_COMPRESS_UTILS_H
#include <limits>
#include <vector>
#include <iterator>
#include "mr_util/error_handling.h"
#include "mr_util/math_utils.h"
namespace mr {
class obitstream
{
private:
std::vector<uint8_t> data;
size_t bitpos;
template<typename T> void put_internal (const T &val, uint8_t bits)
{
int bitsleft = 8-(bitpos&7);
if (bitsleft==8) data.push_back(0);
if (bits<=bitsleft) // wbits==bits
{
data.back() |= ((val&((T(1)<<bits)-1))<<(bitsleft-bits));
bitpos+=bits;
}
else // wbits==bitsleft
{
data.back() |= ((val>>(bits-bitsleft))&((1<<bitsleft)-1));
bitpos+=bitsleft;
put_internal(val,bits-bitsleft);
}
}
public:
obitstream() : bitpos(0) {}
template<typename T> void put (const T &val, uint8_t bits)
{
if (bits==0) return;
MR_assert(bits<=(sizeof(T)<<3),"too many bits for this data type");
// MR_assert(val-(val&((T(1)<<bits)-1))==0,"value has too many bits");
put_internal(val,bits);
}
const std::vector<uint8_t> &state() const
{ return data; }
};
class ibitstream
{
private:
const std::vector<uint8_t> data;
size_t bitpos;
template<typename T> void get_internal (T &val, int bits)
{
int bitsleft = 8-(bitpos&7);
if (bits<=bitsleft)
{
val |= ((data[bitpos>>3]>>(bitsleft-bits))&((1<<bits)-1));
bitpos+=bits;
}
else
{
val |= T((data[bitpos>>3])&((1<<bitsleft)-1))<<(bits-bitsleft);
bitpos+=bitsleft;
get_internal(val,bits-bitsleft);
}
}
public:
ibitstream(const std::vector<uint8_t> &indata) : data(indata), bitpos(0) {}
template<typename T> T get (uint8_t bits)
{
if (bits==0) return T(0);
MR_assert(bits<=(sizeof(T)<<3),"too many bits for this data type");
MR_assert((bitpos+bits)<=8*data.size(),"reading past end of stream");
T res=T(0);
get_internal(res,bits);
return res;
}
void rewind (uint8_t bits) {bitpos-=bits;}
};
template<typename T,bool> struct notNegativeHelper__
{ notNegativeHelper__(const T &) {} };
template<typename T> struct notNegativeHelper__<T,true>
{
notNegativeHelper__(const T &val)
{ MR_assert(val>=T(0),"numbers must be nonnegative");}
};
template<typename T> void assertNotNegative(const T &val)
{ notNegativeHelper__<T,std::numeric_limits<T>::is_signed> dummy(val); }
template<typename Iter> void interpol_encode2 (Iter l, Iter r, obitstream &obs,
int shift)
{
if (r-l<=1) return;
typedef std::iterator_traits<Iter> traits;
typedef typename traits::value_type T;
Iter m=l+(r-l)/2;
T nval = ((*r-*l)>>shift) - (r-l) + 1;
if (nval<=1) return;
uint8_t nb = 1+ilog2(nval-1);
T val = ((*m)>>shift)-(((*l)>>shift)+(m-l));
T nshort=(T(1)<<nb)-nval;
#if 0
// optional rotation
T nrot=nval-(nshort>>1);
if (val>=nrot)
val-=nrot;
else
val+=(nval-nrot);
#endif
if (val<nshort)
obs.put(val,nb-1);
else
obs.put(val+nshort,nb);
interpol_encode2(l,m,obs,shift);
interpol_encode2(m,r,obs,shift);
}
template<typename Iter> void interpol_encode (Iter l, Iter r, obitstream &obs)
{
typedef std::iterator_traits<Iter> traits;
typedef typename traits::value_type T;
if (l==r) // empty range
{ obs.put(0,8); return; }
assertNotNegative(*l);
T combo=*l;
for (Iter i=l+1; i!=r; ++i)
{
MR_assert(*i>*(i-1),"numbers not strictly increasing");
combo|=*i;
}
int shift = trailingZeros(combo);
T maxnum=(*(r-1))>>shift;
if (T(r-l)>maxnum) maxnum=T(r-l);
uint8_t maxbits=1+ilog2(maxnum);
obs.put(maxbits,8);
obs.put(shift,8);
obs.put(r-l,maxbits);
obs.put((*l)>>shift,maxbits);
if (r-l==1) return;
obs.put((*(r-1))>>shift,maxbits);
interpol_encode2(l,r-1,obs,shift);
}
template<typename Iter> void interpol_decode2 (Iter l, Iter r, ibitstream &ibs,
int shift)
{
if (r-l<=1) return;
typedef std::iterator_traits<Iter> traits;
typedef typename traits::value_type T;
Iter m=l+(r-l)/2;
T nval = ((*r-*l)>>shift) - (r-l) + 1;
T val=0;
if (nval>1)
{
uint8_t nb = 1+ilog2(nval-1);
T nshort=(T(1)<<nb)-nval;
val=ibs.get<T>(nb-1);
if (val>=nshort)
val=(val<<1)+ ibs.get<T>(1) - nshort;
#if 0
// optional rotation
T nrot=nval-(nshort>>1);
if (val<(nval-nrot))
val+=nrot;
else
val-=(nval-nrot);
#endif
}
*m=*l+(((m-l)+val)<<shift);
interpol_decode2(l,m,ibs,shift);
interpol_decode2(m,r,ibs,shift);
}
template<typename T> void interpol_decode (std::vector<T> &v, ibitstream &ibs)
{
uint8_t maxbits=ibs.get<uint8_t>(8);
if (maxbits==0) { v.clear(); return; }
int shift = ibs.get<int>(8);
v.resize(ibs.get<size_t>(maxbits));
v[0]=ibs.get<T>(maxbits)<<shift;
if (v.size()==1) return;
v[v.size()-1]=ibs.get<T>(maxbits)<<shift;
interpol_decode2(v.begin(),v.end()-1,ibs,shift);
}
}
#endif
mr_util/crangeset.h 0 → 100644
View file @ 714edde6
/*
* This file is part of libcxxsupport.
*
* libcxxsupport 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.
*
* libcxxsupport 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 libcxxsupport; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file crangeset.h
* Class for storing sets of ranges of integer numbers
*
* Copyright (C) 2015-2020 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef MRUTIL_CRANGESET_H
#define MRUTIL_CRANGESET_H
#include <algorithm>
#include <vector>
#include <utility>
#include <iostream>
#include "mr_util/error_handling.h"
#include "mr_util/math_utils.h"
namespace mr {
/*
compact rangeset (CRS)
The underlying type T must be signed.
BUT: All "on" ranges must lie entirely in N_0.
A nonnegative number indicates the start of an "on" range at this number.
A negative number indicates the start of an "off" range at the absolute value
of this number.
All numbers r are sorted in a way that |r_i| < |r_{i+1}|.
If two consecutive numbers are both nonnegative or negative, it means that the
interval between them contains in fact _two_ ranges: one of length 1 and the
other filling the remainder.
Consequences:
- The first number in CRS must be nonnegative
- numbers r_i and r_{i+1} in the CRS have the property
|r_{i+1}| > |r_i| + 1 (because of the special treatment of short intervals)
Example:
The range set
[5;7[; [10;15[; [16;22[; [23;24[; [25;30[; [35;36[; [40;45[
would be encoded as
5 -7 10 -15 -22 -24 -30 35 40 -45
*/
/*! Class for storing sets of ranges of integer numbers
T must be a signed integer type, but all numbers entered into the range set
must be nonnegative! */
template<typename T> class crangeset
{
private:
typedef std::vector<T> rtype;
rtype r;
struct abscomp
{
bool operator()(const T &a, const T &b)
{
using namespace std;
return abs(a)<abs(b);
}
};
class RsIter
{
private:
size_t idx;
bool extra;
const crangeset &ref;
T val;
bool singleVal() const
{
if (idx==ref.r.size()-1)
return (ref.r[idx]>=0);
else
return ((ref.r[idx]^ref.r[idx+1])>=0);
}
public:
RsIter(const crangeset &ref_) : idx(0), extra(false), ref(ref_), val(0)
{
using namespace std;
if (!atEnd())
val=abs(ref.r[0]);