Commit 09902695 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

add pyHealpix; still need to fix warnings

parent d7af95c1
/*
* 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_cxx/healpix_base.h"
#include "mr_util/geom_utils.h"
#include "mr_util/constants.h"
#include "mr_util/mav.h"
#include "mr_util/space_filling.h"
using namespace std;
namespace mr {
namespace healpix {
namespace {
template<typename T> inline T spread(int v);
template<> inline int spread(int v)
{ return spread_bits_2D_32(v); }
template<> inline int64_t spread(int v)
{ return spread_bits_2D_64(v); }
template<typename T> inline T interleave(int x, int y);
template<> inline int interleave(int x, int y)
{ return coord2morton2D_32({uint32_t(x), uint32_t(y)}); }
template<> inline int64_t interleave(int x, int y)
{ return coord2morton2D_64({uint32_t(x), uint32_t(y)}); }
template<typename T> inline void deinterleave(T v, int &x, int &y);
template<> inline void deinterleave(int v, int &x, int &y)
{
auto res = morton2coord2D_32(v);
x = res[0];
y = res[1];
}
template<> inline void deinterleave(int64_t v, int &x, int &y)
{
auto res = morton2coord2D_64(v);
x = res[0];
y = res[1];
}
}
template<typename I> int T_Healpix_Base<I>::nside2order (I nside)
{
MR_assert (nside>I(0), "invalid value for Nside");
return ((nside)&(nside-1)) ? -1 : ilog2(nside);
}
template<typename I> I T_Healpix_Base<I>::npix2nside (I npix)
{
I res=isqrt(npix/I(12));
MR_assert (npix==res*res*I(12), "invalid value for npix");
return res;
}
template<typename I> I T_Healpix_Base<I>::ring_above (double z) const
{
double az=abs(z);
if (az<=twothird) // equatorial region
return I(nside_*(2-1.5*z));
I iring = I(nside_*sqrt(3*(1-az)));
return (z>0) ? iring : 4*nside_-iring-1;
}
namespace {
/* Short note on the "zone":
zone = 0: pixel lies completely outside the queried shape
1: pixel may overlap with the shape, pixel center is outside
2: pixel center is inside the shape, but maybe not the complete pixel
3: pixel lies completely inside the shape */
template<typename I, typename I2> inline void check_pixel (int o, int order_,
int omax, int zone, rangeset<I2> &pixset, I pix, vector<pair<I,int> > &stk,
bool inclusive, int &stacktop)
{
if (zone==0) return;
if (o<order_)
{
if (zone>=3)
{
int sdist=2*(order_-o); // the "bit-shift distance" between map orders
pixset.append(pix<<sdist,(pix+1)<<sdist); // output all subpixels
}
else // (1<=zone<=2)
for (int i=0; i<4; ++i)
stk.push_back(make_pair(4*pix+3-i,o+1)); // add children
}
else if (o>order_) // this implies that inclusive==true
{
if (zone>=2) // pixel center in shape
{
pixset.append(pix>>(2*(o-order_))); // output the parent pixel at order_
stk.resize(stacktop); // unwind the stack
}
else // (zone==1): pixel center in safety range
{
if (o<omax) // check sublevels
for (int i=0; i<4; ++i) // add children in reverse order
stk.push_back(make_pair(4*pix+3-i,o+1));
else // at resolution limit
{
pixset.append(pix>>(2*(o-order_))); // output the parent pixel at order_
stk.resize(stacktop); // unwind the stack
}
}
}
else // o==order_
{
if (zone>=2)
pixset.append(pix);
else if (inclusive) // and (zone>=1)
{
if (order_<omax) // check sublevels
{
stacktop=stk.size(); // remember current stack position
for (int i=0; i<4; ++i) // add children in reverse order
stk.push_back(make_pair(4*pix+3-i,o+1));
}
else // at resolution limit
pixset.append(pix); // output the pixel
}
}
}
template<typename I> bool check_pixel_ring (const T_Healpix_Base<I> &b1,
const T_Healpix_Base<I> &b2, I pix, I nr, I ipix1, int fct,
double cz, double cphi, double cosrp2, I cpix)
{
if (pix>=nr) pix-=nr;
if (pix<0) pix+=nr;
pix+=ipix1;
if (pix==cpix) return false; // disk center in pixel => overlap
int px,py,pf;
b1.pix2xyf(pix,px,py,pf);
for (int i=0; i<fct-1; ++i) // go along the 4 edges
{
I ox=fct*px, oy=fct*py;
double pz,pphi;
b2.pix2zphi(b2.xyf2pix(ox+i,oy,pf),pz,pphi);
if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2) // overlap
return false;
b2.pix2zphi(b2.xyf2pix(ox+fct-1,oy+i,pf),pz,pphi);
if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2) // overlap
return false;
b2.pix2zphi(b2.xyf2pix(ox+fct-1-i,oy+fct-1,pf),pz,pphi);
if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2) // overlap
return false;
b2.pix2zphi(b2.xyf2pix(ox,oy+fct-1-i,pf),pz,pphi);
if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2) // overlap
return false;
}
return true;
}
} // unnamed namespace
template<typename I> template<typename I2>
void T_Healpix_Base<I>::query_disc_internal
(pointing ptg, double radius, int fact, rangeset<I2> &pixset) const
{
bool inclusive = (fact!=0);
pixset.clear();
ptg.normalize();
if (scheme_==RING)
{
I fct=1;
if (inclusive)
{
MR_assert (((I(1)<<order_max)/nside_)>=fact,
"invalid oversampling factor");
fct = fact;
}
T_Healpix_Base b2;
double rsmall, rbig;
if (fct>1)
{
b2.SetNside(fct*nside_,RING);
rsmall = radius+b2.max_pixrad();
rbig = radius+max_pixrad();
}
else
rsmall = rbig = inclusive ? radius+max_pixrad() : radius;
if (rsmall>=pi)
{ pixset.append(0,npix_); return; }
rbig = min(pi,rbig);
double cosrsmall = cos(rsmall);
double cosrbig = cos(rbig);
double z0 = cos(ptg.theta);
double xa = 1./sqrt((1-z0)*(1+z0));
I cpix=zphi2pix(z0,ptg.phi);
double rlat1 = ptg.theta - rsmall;
double zmax = cos(rlat1);
I irmin = ring_above (zmax)+1;
if ((rlat1<=0) && (irmin>1)) // north pole in the disk
{
I sp,rp; bool dummy;
get_ring_info_small(irmin-1,sp,rp,dummy);
pixset.append(0,sp+rp);
}
if ((fct>1) && (rlat1>0)) irmin=max(I(1),irmin-1);
double rlat2 = ptg.theta + rsmall;
double zmin = cos(rlat2);
I irmax = ring_above (zmin);
if ((fct>1) && (rlat2<pi)) irmax=min(4*nside_-1,irmax+1);
for (I iz=irmin; iz<=irmax; ++iz)
{
double z=ring2z(iz);
double x = (cosrbig-z*z0)*xa;
double ysq = 1-z*z-x*x;
double dphi=-1;
if (ysq<=0) // no intersection, ring completely inside or outside
dphi = (fct==1) ? 0: pi-1e-15;
else
dphi = atan2(sqrt(ysq),x);
if (dphi>0)
{
I nr, ipix1;
bool shifted;
get_ring_info_small(iz,ipix1,nr,shifted);
double shift = shifted ? 0.5 : 0.;
I ipix2 = ipix1 + nr - 1; // highest pixel number in the ring
I ip_lo = ifloor<I>(nr*inv_twopi*(ptg.phi-dphi) - shift)+1;
I ip_hi = ifloor<I>(nr*inv_twopi*(ptg.phi+dphi) - shift);
if (fct>1)
{
while ((ip_lo<=ip_hi) && check_pixel_ring
(*this,b2,ip_lo,nr,ipix1,fct,z0,ptg.phi,cosrsmall,cpix))
++ip_lo;
while ((ip_hi>ip_lo) && check_pixel_ring
(*this,b2,ip_hi,nr,ipix1,fct,z0,ptg.phi,cosrsmall,cpix))
--ip_hi;
}
if (ip_lo<=ip_hi)
{
if (ip_hi>=nr)
{ ip_lo-=nr; ip_hi-=nr; }
if (ip_lo<0)
{
pixset.append(ipix1,ipix1+ip_hi+1);
pixset.append(ipix1+ip_lo+nr,ipix2+1);
}
else
pixset.append(ipix1+ip_lo,ipix1+ip_hi+1);
}
}
}
if ((rlat2>=pi) && (irmax+1<4*nside_)) // south pole in the disk
{
I sp,rp; bool dummy;
get_ring_info_small(irmax+1,sp,rp,dummy);
pixset.append(sp,npix_);
}
}
else // scheme_==NEST
{
if (radius>=pi) // disk covers the whole sphere
{ pixset.append(0,npix_); return; }
int oplus = 0;
if (inclusive)
{
MR_assert ((I(1)<<(order_max-order_))>=fact,
"invalid oversampling factor");
MR_assert ((fact&(fact-1))==0,
"oversampling factor must be a power of 2");
oplus=ilog2(fact);
}
int omax=order_+oplus; // the order up to which we test
vec3 vptg(ptg);
vector<T_Healpix_Base<I> > base(omax+1);
vector<double> crpdr(omax+1), crmdr(omax+1);
double cosrad=cos(radius);
for (int o=0; o<=omax; ++o) // prepare data at the required orders
{
base[o].Set(o,NEST);
double dr=base[o].max_pixrad(); // safety distance
crpdr[o] = (radius+dr>pi) ? -1. : cos(radius+dr);
crmdr[o] = (radius-dr<0.) ? 1. : cos(radius-dr);
}
vector<pair<I,int> > stk; // stack for pixel numbers and their orders
stk.reserve(12+3*omax); // reserve maximum size to avoid reallocation
for (int i=0; i<12; ++i) // insert the 12 base pixels in reverse order
stk.push_back(make_pair(I(11-i),0));
int stacktop=0; // a place to save a stack position
while (!stk.empty()) // as long as there are pixels on the stack
{
// pop current pixel number and order from the stack
I pix=stk.back().first;
int o=stk.back().second;
stk.pop_back();
double z,phi;
base[o].pix2zphi(pix,z,phi);
// cosine of angular distance between pixel center and disk center
double cangdist=cosdist_zphi(vptg.z,ptg.phi,z,phi);
if (cangdist>crpdr[o])
{
int zone = (cangdist<cosrad) ? 1 : ((cangdist<=crmdr[o]) ? 2 : 3);
check_pixel (o, order_, omax, zone, pixset, pix, stk, inclusive,
stacktop);
}
}
}
}
template<typename I> void T_Healpix_Base<I>::query_disc
(pointing ptg, double radius, rangeset<I> &pixset) const
{
query_disc_internal (ptg, radius, 0, pixset);
}
template<typename I> void T_Healpix_Base<I>::query_disc_inclusive
(pointing ptg, double radius, rangeset<I> &pixset, int fact) const
{
MR_assert(fact>0,"fact must be a positive integer");
if ((sizeof(I)<8) && (((I(1)<<order_max)/nside_)<fact))
{
T_Healpix_Base<int64_t> base2(nside_,scheme_,SET_NSIDE);
base2.query_disc_internal(ptg,radius,fact,pixset);
return;
}
query_disc_internal (ptg, radius, fact, pixset);
}
template<typename I> template<typename I2>
void T_Healpix_Base<I>::query_multidisc (const vector<vec3> &norm,
const vector<double> &rad, int fact, rangeset<I2> &pixset) const
{
bool inclusive = (fact!=0);
size_t nv=norm.size();
MR_assert(nv==rad.size(),"inconsistent input arrays");
pixset.clear();
if (scheme_==RING)
{
I fct=1;
if (inclusive)
{
MR_assert (((I(1)<<order_max)/nside_)>=fact,
"invalid oversampling factor");
fct = fact;
}
T_Healpix_Base b2;
double rpsmall, rpbig;
if (fct>1)
{
b2.SetNside(fct*nside_,RING);
rpsmall = b2.max_pixrad();
rpbig = max_pixrad();
}
else
rpsmall = rpbig = inclusive ? max_pixrad() : 0;
I irmin=1, irmax=4*nside_-1;
vector<double> z0,xa,cosrsmall,cosrbig;
vector<pointing> ptg;
vector<I> cpix;
for (size_t i=0; i<nv; ++i)
{
double rsmall=rad[i]+rpsmall;
if (rsmall<pi)
{
double rbig=min(pi,rad[i]+rpbig);
pointing pnt=pointing(norm[i]);
cosrsmall.push_back(cos(rsmall));
cosrbig.push_back(cos(rbig));
double cth=cos(pnt.theta);
z0.push_back(cth);
if (fct>1) cpix.push_back(zphi2pix(cth,pnt.phi));
xa.push_back(1./sqrt((1-cth)*(1+cth)));
ptg.push_back(pnt);
double rlat1 = pnt.theta - rsmall;
double zmax = cos(rlat1);
I irmin_t = (rlat1<=0) ? 1 : ring_above (zmax)+1;
if ((fct>1) && (rlat1>0)) irmin_t=max(I(1),irmin_t-1);
double rlat2 = pnt.theta + rsmall;
double zmin = cos(rlat2);
I irmax_t = (rlat2>=pi) ? 4*nside_-1 : ring_above (zmin);
if ((fct>1) && (rlat2<pi)) irmax_t=min(4*nside_-1,irmax_t+1);
if (irmax_t < irmax) irmax=irmax_t;
if (irmin_t > irmin) irmin=irmin_t;
}
}
for (I iz=irmin; iz<=irmax; ++iz)
{
double z=ring2z(iz);
I ipix1,nr;
bool shifted;
get_ring_info_small(iz,ipix1,nr,shifted);
double shift = shifted ? 0.5 : 0.;
rangeset<I2> tr;
tr.append(ipix1,ipix1+nr);
for (size_t j=0; j<z0.size(); ++j)
{
double x = (cosrbig[j]-z*z0[j])*xa[j];
double ysq = 1.-z*z-x*x;
double dphi = (ysq<=0) ? pi-1e-15 : atan2(sqrt(ysq),x);
I ip_lo = ifloor<I>(nr*inv_twopi*(ptg[j].phi-dphi) - shift)+1;
I ip_hi = ifloor<I>(nr*inv_twopi*(ptg[j].phi+dphi) - shift);
if (fct>1)
{
while ((ip_lo<=ip_hi) && check_pixel_ring
(*this,b2,ip_lo,nr,ipix1,fct,z0[j],ptg[j].phi,cosrsmall[j],cpix[j]))
++ip_lo;
while ((ip_hi>ip_lo) && check_pixel_ring
(*this,b2,ip_hi,nr,ipix1,fct,z0[j],ptg[j].phi,cosrsmall[j],cpix[j]))
--ip_hi;
}
if (ip_hi>=nr)
{ ip_lo-=nr; ip_hi-=nr;}
if (ip_lo<0)
tr.remove(ipix1+ip_hi+1,ipix1+ip_lo+nr);
else
tr.intersect(ipix1+ip_lo,ipix1+ip_hi+1);
}
pixset.append(tr);
}
}
else // scheme_ == NEST
{
int oplus = 0;
if (inclusive)
{
MR_assert ((I(1)<<(order_max-order_))>=fact,
"invalid oversampling factor");
MR_assert ((fact&(fact-1))==0,
"oversampling factor must be a power of 2");
oplus=ilog2(fact);
}
int omax=order_+oplus; // the order up to which we test
// TODO: ignore all disks with radius>=pi
vector<T_Healpix_Base<I> > base(omax+1);
mav<double,3> crlimit({size_t(omax)+1,nv,3});
for (int o=0; o<=omax; ++o) // prepare data at the required orders
{
base[o].Set(o,NEST);
double dr=base[o].max_pixrad(); // safety distance
for (size_t i=0; i<nv; ++i)
{
crlimit(o,i,0) = (rad[i]+dr>pi) ? -1. : cos(rad[i]+dr);
crlimit(o,i,1) = (o==0) ? cos(rad[i]) : crlimit(0,i,1);
crlimit(o,i,2) = (rad[i]-dr<0.) ? 1. : cos(rad[i]-dr);
}
}
vector<pair<I,int> > stk; // stack for pixel numbers and their orders
stk.reserve(12+3*omax); // reserve maximum size to avoid reallocation
for (int i=0; i<12; ++i) // insert the 12 base pixels in reverse order
stk.push_back(make_pair(I(11-i),0));
int stacktop=0; // a place to save a stack position
while (!stk.empty()) // as long as there are pixels on the stack
{
// pop current pixel number and order from the stack
I pix=stk.back().first;
int o=stk.back().second;
stk.pop_back();
vec3 pv(base[o].pix2vec(pix));
size_t zone=3;
for (size_t i=0; i<nv; ++i)
{
double crad=dotprod(pv,norm[i]);
for (size_t iz=0; iz<zone; ++iz)
if (crad<crlimit(o,i,iz))
if ((zone=iz)==0) goto bailout;
}
check_pixel (o, order_, omax, zone, pixset, pix, stk, inclusive,
stacktop);
bailout:;
}
}
}
template<typename I> void T_Healpix_Base<I>::query_multidisc_general
(const vector<vec3> &norm, const vector<double> &rad, bool inclusive,
const vector<int> &cmds, rangeset<I> &pixset) const
{
size_t nv=norm.size();
MR_assert(nv==rad.size(),"inconsistent input arrays");
pixset.clear();
if (scheme_==RING)
{
MR_fail ("not yet implemented");
}
else // scheme_ == NEST
{
int oplus=inclusive ? 2 : 0;
int omax=min<int>(order_max,order_+oplus); // the order up to which we test
// TODO: ignore all disks with radius>=pi
vector<T_Healpix_Base<I> > base(omax+1);
mav<double,3> crlimit({size_t(omax+1),nv,3});
for (int o=0; o<=omax; ++o) // prepare data at the required orders
{
base[o].Set(o,NEST);
double dr=base[o].max_pixrad(); // safety distance
for (size_t i=0; i<nv; ++i)
{
crlimit(o,i,0) = (rad[i]+dr>pi) ? -1. : cos(rad[i]+dr);
crlimit(o,i,1) = (o==0) ? cos(rad[i]) : crlimit(0,i,1);
crlimit(o,i,2) = (rad[i]-dr<0.) ? 1. : cos(rad[i]-dr);
}
}
vector<pair<I,int> > stk; // stack for pixel numbers and their orders
stk.reserve(12+3*omax); // reserve maximum size to avoid reallocation
for (int i=0; i<12; ++i) // insert the 12 base pixels in reverse order
stk.push_back(make_pair(I(11-i),0));
int stacktop=0; // a place to save a stack position
vector<size_t> zone(nv);
vector<size_t> zstk; zstk.reserve(cmds.size());
while (!stk.empty()) // as long as there are pixels on the stack
{
// pop current pixel number and order from the stack
I pix=stk.back().first;
int o=stk.back().second;
stk.pop_back();
vec3 pv(base[o].pix2vec(pix));
for (size_t i=0; i<nv; ++i)