33 #include "geom_utils.h" 34 #include "lsconstants.h" 40 planck_assert (nside>I(0),
"invalid value for Nside");
41 return ((nside)&(nside-1)) ? -1 : ilog2(nside);
45 I res=isqrt(npix/I(12));
46 planck_assert (npix==res*res*I(12),
"invalid value for npix");
54 return I(nside_*(2-1.5*z));
55 I iring = I(nside_*sqrt(3*(1-az)));
56 return (z>0) ? iring : 4*nside_-iring-1;
67 template<
typename I,
typename I2>
inline void check_pixel (
int o,
int order_,
68 int omax,
int zone, rangeset<I2> &pixset, I pix, vector<pair<I,int> > &stk,
69 bool inclusive,
int &stacktop)
77 int sdist=2*(order_-o);
78 pixset.append(pix<<sdist,(pix+1)<<sdist);
81 for (
int i=0; i<4; ++i)
82 stk.push_back(make_pair(4*pix+3-i,o+1));
88 pixset.append(pix>>(2*(o-order_)));
94 for (
int i=0; i<4; ++i)
95 stk.push_back(make_pair(4*pix+3-i,o+1));
98 pixset.append(pix>>(2*(o-order_)));
112 for (
int i=0; i<4; ++i)
113 stk.push_back(make_pair(4*pix+3-i,o+1));
123 double cz,
double cphi,
double cosrp2, I cpix)
125 if (pix>=nr) pix-=nr;
128 if (pix==cpix)
return false;
130 b1.pix2xyf(pix,px,py,pf);
131 for (
int i=0; i<fct-1; ++i)
133 I ox=fct*px, oy=fct*py;
135 b2.
pix2zphi(b2.xyf2pix(ox+i,oy,pf),pz,pphi);
136 if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2)
138 b2.
pix2zphi(b2.xyf2pix(ox+fct-1,oy+i,pf),pz,pphi);
139 if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2)
141 b2.
pix2zphi(b2.xyf2pix(ox+fct-1-i,oy+fct-1,pf),pz,pphi);
142 if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2)
144 b2.
pix2zphi(b2.xyf2pix(ox,oy+fct-1-i,pf),pz,pphi);
145 if (cosdist_zphi(pz,pphi,cz,cphi)>cosrp2)
153 template<
typename I>
template<
typename I2>
155 (pointing ptg,
double radius,
int fact, rangeset<I2> &pixset)
const 157 bool inclusive = (fact!=0);
166 planck_assert (((I(1)<<order_max)/nside_)>=fact,
167 "invalid oversampling factor");
176 rbig = radius+max_pixrad();
179 rsmall = rbig = inclusive ? radius+max_pixrad() : radius;
182 { pixset.append(0,npix_);
return; }
186 double cosrsmall = cos(rsmall);
187 double cosrbig = cos(rbig);
189 double z0 = cos(ptg.theta);
190 double xa = 1./sqrt((1-z0)*(1+z0));
192 I cpix=zphi2pix(z0,ptg.phi);
194 double rlat1 = ptg.theta - rsmall;
195 double zmax = cos(rlat1);
196 I irmin = ring_above (zmax)+1;
198 if ((rlat1<=0) && (irmin>1))
201 get_ring_info_small(irmin-1,sp,rp,dummy);
202 pixset.append(0,sp+rp);
205 if ((fct>1) && (rlat1>0)) irmin=max(I(1),irmin-1);
207 double rlat2 = ptg.theta + rsmall;
208 double zmin = cos(rlat2);
209 I irmax = ring_above (zmin);
211 if ((fct>1) && (rlat2<pi)) irmax=min(4*nside_-1,irmax+1);
213 for (I iz=irmin; iz<=irmax; ++iz)
216 double x = (cosrbig-z*z0)*xa;
217 double ysq = 1-z*z-x*x;
220 dphi = (fct==1) ? 0: pi-1e-15;
222 dphi = atan2(sqrt(ysq),x);
227 get_ring_info_small(iz,ipix1,nr,shifted);
228 double shift = shifted ? 0.5 : 0.;
230 I ipix2 = ipix1 + nr - 1;
232 I ip_lo = ifloor<I>(nr*inv_twopi*(ptg.phi-dphi) - shift)+1;
233 I ip_hi = ifloor<I>(nr*inv_twopi*(ptg.phi+dphi) - shift);
237 while ((ip_lo<=ip_hi) && check_pixel_ring
238 (*
this,b2,ip_lo,nr,ipix1,fct,z0,ptg.phi,cosrsmall,cpix))
240 while ((ip_hi>ip_lo) && check_pixel_ring
241 (*
this,b2,ip_hi,nr,ipix1,fct,z0,ptg.phi,cosrsmall,cpix))
248 { ip_lo-=nr; ip_hi-=nr; }
251 pixset.append(ipix1,ipix1+ip_hi+1);
252 pixset.append(ipix1+ip_lo+nr,ipix2+1);
255 pixset.append(ipix1+ip_lo,ipix1+ip_hi+1);
259 if ((rlat2>=pi) && (irmax+1<4*nside_))
262 get_ring_info_small(irmax+1,sp,rp,dummy);
263 pixset.append(sp,npix_);
269 { pixset.append(0,npix_);
return; }
274 planck_assert ((I(1)<<(order_max-order_))>=fact,
275 "invalid oversampling factor");
276 planck_assert ((fact&(fact-1))==0,
277 "oversampling factor must be a power of 2");
280 int omax=order_+oplus;
283 arr<T_Healpix_Base<I> > base(omax+1);
284 arr<double> crpdr(omax+1), crmdr(omax+1);
285 double cosrad=cos(radius);
286 for (
int o=0; o<=omax; ++o)
289 double dr=base[o].max_pixrad();
290 crpdr[o] = (radius+dr>pi) ? -1. : cos(radius+dr);
291 crmdr[o] = (radius-dr<0.) ? 1. : cos(radius-dr);
293 vector<pair<I,int> > stk;
294 stk.reserve(12+3*omax);
295 for (
int i=0; i<12; ++i)
296 stk.push_back(make_pair(I(11-i),0));
303 I pix=stk.back().first;
304 int o=stk.back().second;
308 base[o].pix2zphi(pix,z,phi);
310 double cangdist=cosdist_zphi(vptg.z,ptg.phi,z,phi);
312 if (cangdist>crpdr[o])
314 int zone = (cangdist<cosrad) ? 1 : ((cangdist<=crmdr[o]) ? 2 : 3);
316 check_pixel (o, order_, omax, zone, pixset, pix, stk, inclusive,
324 (pointing ptg,
double radius, rangeset<I> &pixset)
const 326 query_disc_internal (ptg, radius, 0, pixset);
330 (pointing ptg,
double radius, rangeset<I> &pixset,
int fact)
const 332 planck_assert(fact>0,
"fact must be a positive integer");
333 if ((
sizeof(I)<8) && (((I(1)<<order_max)/nside_)<fact))
336 base2.query_disc_internal(ptg,radius,fact,pixset);
339 query_disc_internal (ptg, radius, fact, pixset);
342 template<
typename I>
template<
typename I2>
344 const arr<double> &rad,
int fact, rangeset<I2> &pixset)
const 346 bool inclusive = (fact!=0);
347 tsize nv=norm.size();
348 planck_assert(nv==rad.size(),
"inconsistent input arrays");
356 planck_assert (((I(1)<<order_max)/nside_)>=fact,
357 "invalid oversampling factor");
361 double rpsmall, rpbig;
366 rpbig = max_pixrad();
369 rpsmall = rpbig = inclusive ? max_pixrad() : 0;
371 I irmin=1, irmax=4*nside_-1;
372 vector<double> z0,xa,cosrsmall,cosrbig;
373 vector<pointing> ptg;
375 for (tsize i=0; i<nv; ++i)
377 double rsmall=rad[i]+rpsmall;
380 double rbig=min(pi,rad[i]+rpbig);
381 pointing pnt=pointing(norm[i]);
382 cosrsmall.push_back(cos(rsmall));
383 cosrbig.push_back(cos(rbig));
384 double cth=cos(pnt.theta);
386 if (fct>1) cpix.push_back(zphi2pix(cth,pnt.phi));
387 xa.push_back(1./sqrt((1-cth)*(1+cth)));
390 double rlat1 = pnt.theta - rsmall;
391 double zmax = cos(rlat1);
392 I irmin_t = (rlat1<=0) ? 1 : ring_above (zmax)+1;
394 if ((fct>1) && (rlat1>0)) irmin_t=max(I(1),irmin_t-1);
396 double rlat2 = pnt.theta + rsmall;
397 double zmin = cos(rlat2);
398 I irmax_t = (rlat2>=pi) ? 4*nside_-1 : ring_above (zmin);
400 if ((fct>1) && (rlat2<pi)) irmax_t=min(4*nside_-1,irmax_t+1);
402 if (irmax_t < irmax) irmax=irmax_t;
403 if (irmin_t > irmin) irmin=irmin_t;
407 for (I iz=irmin; iz<=irmax; ++iz)
412 get_ring_info_small(iz,ipix1,nr,shifted);
413 double shift = shifted ? 0.5 : 0.;
415 tr.append(ipix1,ipix1+nr);
416 for (tsize j=0; j<z0.size(); ++j)
418 double x = (cosrbig[j]-z*z0[j])*xa[j];
419 double ysq = 1.-z*z-x*x;
420 double dphi = (ysq<=0) ? pi-1e-15 : atan2(sqrt(ysq),x);
421 I ip_lo = ifloor<I>(nr*inv_twopi*(ptg[j].phi-dphi) - shift)+1;
422 I ip_hi = ifloor<I>(nr*inv_twopi*(ptg[j].phi+dphi) - shift);
425 while ((ip_lo<=ip_hi) && check_pixel_ring
426 (*
this,b2,ip_lo,nr,ipix1,fct,z0[j],ptg[j].phi,cosrsmall[j],cpix[j]))
428 while ((ip_hi>ip_lo) && check_pixel_ring
429 (*
this,b2,ip_hi,nr,ipix1,fct,z0[j],ptg[j].phi,cosrsmall[j],cpix[j]))
433 { ip_lo-=nr; ip_hi-=nr;}
435 tr.remove(ipix1+ip_hi+1,ipix1+ip_lo+nr);
437 tr.intersect(ipix1+ip_lo,ipix1+ip_hi+1);
447 planck_assert ((I(1)<<(order_max-order_))>=fact,
448 "invalid oversampling factor");
449 planck_assert ((fact&(fact-1))==0,
450 "oversampling factor must be a power of 2");
453 int omax=order_+oplus;
457 arr<T_Healpix_Base<I> > base(omax+1);
458 arr3<double> crlimit(omax+1,nv,3);
459 for (
int o=0; o<=omax; ++o)
462 double dr=base[o].max_pixrad();
463 for (tsize i=0; i<nv; ++i)
465 crlimit(o,i,0) = (rad[i]+dr>pi) ? -1. : cos(rad[i]+dr);
466 crlimit(o,i,1) = (o==0) ? cos(rad[i]) : crlimit(0,i,1);
467 crlimit(o,i,2) = (rad[i]-dr<0.) ? 1. : cos(rad[i]-dr);
471 vector<pair<I,int> > stk;
472 stk.reserve(12+3*omax);
473 for (
int i=0; i<12; ++i)
474 stk.push_back(make_pair(I(11-i),0));
481 I pix=stk.back().first;
482 int o=stk.back().second;
485 vec3 pv(base[o].pix2vec(pix));
488 for (tsize i=0; i<nv; ++i)
490 double crad=dotprod(pv,norm[i]);
491 for (tsize iz=0; iz<zone; ++iz)
492 if (crad<crlimit(o,i,iz))
493 if ((zone=iz)==0)
goto bailout;
496 check_pixel (o, order_, omax, zone, pixset, pix, stk, inclusive,
504 (
const arr<vec3> &norm,
const arr<double> &rad,
bool inclusive,
505 const vector<int> &cmds, rangeset<I> &pixset)
const 507 tsize nv=norm.size();
508 planck_assert(nv==rad.size(),
"inconsistent input arrays");
513 planck_fail (
"not yet implemented");
517 int oplus=inclusive ? 2 : 0;
518 int omax=min<int>(order_max,order_+oplus);
522 arr<T_Healpix_Base<I> > base(omax+1);
523 arr3<double> crlimit(omax+1,nv,3);
524 for (
int o=0; o<=omax; ++o)
527 double dr=base[o].max_pixrad();
528 for (tsize i=0; i<nv; ++i)
530 crlimit(o,i,0) = (rad[i]+dr>pi) ? -1. : cos(rad[i]+dr);
531 crlimit(o,i,1) = (o==0) ? cos(rad[i]) : crlimit(0,i,1);
532 crlimit(o,i,2) = (rad[i]-dr<0.) ? 1. : cos(rad[i]-dr);
536 vector<pair<I,int> > stk;
537 stk.reserve(12+3*omax);
538 for (
int i=0; i<12; ++i)
539 stk.push_back(make_pair(I(11-i),0));
544 vector<tsize> zstk; zstk.reserve(cmds.size());
549 I pix=stk.back().first;
550 int o=stk.back().second;
553 vec3 pv(base[o].pix2vec(pix));
555 for (tsize i=0; i<nv; ++i)
558 double crad=dotprod(pv,norm[i]);
559 for (tsize iz=0; iz<zone[i]; ++iz)
560 if (crad<crlimit(o,i,iz))
564 for (tsize i=0; i<cmds.size(); ++i)
570 tmp=zstk.back(); zstk.pop_back();
571 zstk.back() = max(zstk.back(),tmp);
574 tmp=zstk.back(); zstk.pop_back();
575 zstk.back() = min(zstk.back(),tmp);
578 zstk.push_back(zone[cmds[i]]);
581 planck_assert(zstk.size()==1,
"inconsistent commands");
582 tsize zn=zstk[0]; zstk.pop_back();
584 check_pixel (o, order_, omax, zn, pixset, pix, stk, inclusive,
593 {
return utab[v&0xff] | (utab[(v>>8)&0xff]<<16); }
596 return int64(utab[ v &0xff]) | (int64(utab[(v>> 8)&0xff])<<16)
597 | (int64(utab[(v>>16)&0xff])<<32) | (int64(utab[(v>>24)&0xff])<<48);
602 int raw = (v&0x5555) | ((v&0x55550000)>>15);
603 return ctab[raw&0xff] | (ctab[raw>>8]<<4);
607 int64 raw = v&0x5555555555555555ull;
609 return ctab[ raw &0xff] | (ctab[(raw>> 8)&0xff]<< 4)
610 | (ctab[(raw>>32)&0xff]<<16) | (ctab[(raw>>40)&0xff]<<20);
615 #include <x86intrin.h> 618 {
return _pdep_u32(v,0x55555555u); }
620 {
return _pdep_u64(v,0x5555555555555555ull); }
622 {
return _pext_u32(v,0x55555555u); }
624 {
return _pext_u64(v,0x5555555555555555ull); }
629 int &iy,
int &face_num)
const 631 face_num = pix>>(2*order_);
633 ix = compress_bits(pix);
634 iy = compress_bits(pix>>1);
639 {
return (I(face_num)<<(2*order_)) + spread_bits(ix) + (spread_bits(iy)<<1); }
644 template<
typename I>
inline I special_div(I a, I b)
648 return (t<<1)+(a>=b);
656 I iring, iphi, kshift, nr;
661 iring = (1+isqrt(1+2*pix))>>1;
662 iphi = (pix+1) - 2*iring*(iring-1);
665 face_num=special_div(iphi-1,nr);
667 else if (pix<(npix_-ncap_))
670 I tmp = (order_>=0) ? ip>>(order_+2) : ip/(4*nside_);
672 iphi = ip-tmp*4*nside_ + 1;
673 kshift = (iring+nside_)&1;
677 I ifm = iphi - (ire>>1) + nside_ -1,
678 ifp = iphi - (irm>>1) + nside_ -1;
680 { ifm >>= order_; ifp >>= order_; }
682 { ifm /= nside_; ifp /= nside_; }
683 face_num = (ifp==ifm) ? (ifp|4) : ((ifp<ifm) ? ifp : (ifm+8));
688 iring = (1+isqrt(2*ip-1))>>1;
689 iphi = 4*iring + 1 - (ip - 2*iring*(iring-1));
693 face_num=special_div(iphi-1,nr)+8;
696 I irt = iring - ((2+(face_num>>2))*nside_) + 1;
697 I ipt = 2*iphi- jpll[face_num]*nr - kshift -1;
699 if (ipt>=nl2) ipt-=8*nside_;
709 I jr = (jrll[face_num]*nside_) - ix - iy - 1;
711 I nr, kshift, n_before;
714 get_ring_info_small(jr,n_before,nr,shifted);
717 I jp = (jpll[face_num]*nr + ix - iy + 1 + kshift) / 2;
718 planck_assert(jp<=4*nr,
"must not happen");
721 return n_before + jp - 1;
725 : order_(-1), nside_(0), npface_(0), ncap_(0), npix_(0),
726 fact1_(0), fact2_(0), scheme_(
RING) {}
731 planck_assert ((order>=0)&&(order<=order_max),
"bad order");
733 nside_ = I(1)<<order;
734 npface_ = nside_<<order_;
735 ncap_ = (npface_-nside_)<<1;
738 fact1_ = (nside_<<1)*fact2_;
744 order_ = nside2order(nside);
745 planck_assert ((scheme!=
NEST) || (order_>=0),
746 "SetNside: nside must be power of 2 for nested maps");
748 npface_ = nside_*nside_;
749 ncap_ = (npface_-nside_)<<1;
752 fact1_ = (nside_<<1)*fact2_;
759 return 1 - ring*ring*fact2_;
761 return (2*nside_-ring)*fact1_;
762 ring=4*nside_ - ring;
763 return ring*ring*fact2_ - 1;
771 return (1+I(isqrt(1+2*pix)))>>1;
772 else if (pix<(npix_-ncap_))
773 return (pix-ncap_)/(4*nside_) + nside_;
775 return 4*nside_-((1+I(isqrt(2*(npix_-pix)-1)))>>1);
779 int face_num, ix, iy;
780 nest2xyf(pix,ix,iy,face_num);
781 return (I(jrll[face_num])<<order_) - ix - iy - 1;
787 planck_assert(order_>=0,
"hierarchical map required");
788 int ix, iy, face_num;
789 nest2xyf (pix, ix, iy, face_num);
790 return xyf2ring (ix, iy, face_num);
795 planck_assert(order_>=0,
"hierarchical map required");
796 int ix, iy, face_num;
797 ring2xyf (pix, ix, iy, face_num);
798 return xyf2nest (ix, iy, face_num);
802 (I pix,
int dir)
const 804 int face = (pix>>(2*order_));
806 int state = ((peano_face2path[dir][face]<<4))|(dir<<7);
807 int shift=2*order_-4;
808 for (; shift>=0; shift-=4)
810 state=peano_arr2[(state&0xF0) | ((pix>>shift)&0xF)];
811 result = (result<<4) | (state&0xF);
815 state=peano_arr[((state>>2)&0xFC) | (pix&0x3)];
816 result = (result<<2) | (state&0x3);
819 return result + (I(peano_face2face[dir][face])<<(2*order_));
823 {
return nest_peano_helper(pix,0); }
826 {
return nest_peano_helper(pix,1); }
829 double sth,
bool have_sth)
const 832 double tt = fmodulo(phi*inv_halfpi,4.0);
839 double temp1 = nside_*(0.5+tt);
840 double temp2 = nside_*z*0.75;
841 I jp = I(temp1-temp2);
842 I jm = I(temp1+temp2);
845 I ir = nside_ + 1 + jp - jm;
848 I t1 = jp+jm-nside_+kshift+1+nl4+nl4;
850 (t1>>1)&(nl4-1) : ((t1>>1)%nl4);
852 return ncap_ + (ir-1)*nl4 + ip;
856 double tp = tt-I(tt);
857 double tmp = ((za<0.99)||(!have_sth)) ?
858 nside_*sqrt(3*(1-za)) :
859 nside_*sth/sqrt((1.+za)/3.);
862 I jm = I((1.0-tp)*tmp);
866 planck_assert((ip>=0)&&(ip<4*ir),
"must not happen");
869 return (z>0) ? 2*ir*(ir-1) + ip : npix_ - 2*ir*(ir+1) + ip;
876 double temp1 = nside_*(0.5+tt);
877 double temp2 = nside_*(z*0.75);
878 I jp = I(temp1-temp2);
879 I jm = I(temp1+temp2);
880 I ifp = jp >> order_;
881 I ifm = jm >> order_;
882 int face_num = (ifp==ifm) ? (ifp|4) : ((ifp<ifm) ? ifp : (ifm+8));
884 int ix = jm & (nside_-1),
885 iy = nside_ - (jp & (nside_-1)) - 1;
886 return xyf2nest(ix,iy,face_num);
890 int ntt = min(3,
int(tt));
892 double tmp = ((za<0.99)||(!have_sth)) ?
893 nside_*sqrt(3*(1-za)) :
894 nside_*sth/sqrt((1.+za)/3.);
897 I jm = I((1.0-tp)*tmp);
901 xyf2nest(nside_-jm -1,nside_-jp-1,ntt) : xyf2nest(jp,jm,ntt+8);
907 double &phi,
double &sth,
bool &have_sth)
const 914 I iring = (1+I(isqrt(1+2*pix)))>>1;
915 I iphi = (pix+1) - 2*iring*(iring-1);
917 double tmp=(iring*iring)*fact2_;
919 if (z>0.99) { sth=sqrt(tmp*(2.-tmp)); have_sth=
true; }
920 phi = (iphi-0.5) * halfpi/iring;
922 else if (pix<(npix_-ncap_))
926 I tmp = (order_>=0) ? ip>>(order_+2) : ip/nl4;
927 I iring = tmp + nside_,
930 double fodd = ((iring+nside_)&1) ? 1 : 0.5;
932 z = (2*nside_-iring)*fact1_;
933 phi = (iphi-fodd) * pi*0.75*fact1_;
938 I iring = (1+I(isqrt(2*ip-1)))>>1;
939 I iphi = 4*iring + 1 - (ip - 2*iring*(iring-1));
941 double tmp=(iring*iring)*fact2_;
943 if (z<-0.99) { sth=sqrt(tmp*(2.-tmp)); have_sth=
true; }
944 phi = (iphi-0.5) * halfpi/iring;
949 int face_num, ix, iy;
950 nest2xyf(pix,ix,iy,face_num);
952 I jr = (I(jrll[face_num])<<order_) - ix - iy - 1;
958 double tmp=(nr*nr)*fact2_;
960 if (z>0.99) { sth=sqrt(tmp*(2.-tmp)); have_sth=
true; }
962 else if (jr > 3*nside_)
965 double tmp=(nr*nr)*fact2_;
967 if (z<-0.99) { sth=sqrt(tmp*(2.-tmp)); have_sth=
true; }
972 z = (2*nside_-jr)*fact1_;
975 I tmp=I(jpll[face_num])*nr+ix-iy;
976 planck_assert(tmp<8*nr,
"must not happen");
977 if (tmp<0) tmp+=8*nr;
978 phi = (nr==nside_) ? 0.75*halfpi*tmp*fact1_ :
983 template<
typename I>
template<
typename I2>
985 (
const vector<pointing> &vertex,
int fact, rangeset<I2> &pixset)
const 987 bool inclusive = (fact!=0);
988 tsize nv=vertex.size();
989 tsize ncirc = inclusive ? nv+1 : nv;
990 planck_assert(nv>=3,
"not enough vertices in polygon");
992 for (tsize i=0; i<nv; ++i)
993 vv[i]=vertex[i].to_vec3();
994 arr<vec3> normal(ncirc);
996 for (tsize i=0; i<nv; ++i)
998 normal[i]=crossprod(vv[i],vv[(i+1)%nv]).Norm();
999 double hnd=dotprod(normal[i],vv[(i+2)%nv]);
1000 planck_assert(abs(hnd)>1e-10,
"degenerate corner");
1002 flip = (hnd<0.) ? -1 : 1;
1004 planck_assert(flip*hnd>0,
"polygon is not convex");
1007 arr<double> rad(ncirc,halfpi);
1011 find_enclosing_circle (vv, normal[nv], cosrad);
1012 rad[nv]=acos(cosrad);
1014 query_multidisc(normal,rad,fact,pixset);
1018 (
const vector<pointing> &vertex, rangeset<I> &pixset)
const 1020 query_polygon_internal(vertex, 0, pixset);
1024 (
const vector<pointing> &vertex, rangeset<I> &pixset,
int fact)
const 1026 planck_assert(fact>0,
"fact must be a positive integer");
1027 if ((
sizeof(I)<8) && (((I(1)<<order_max)/nside_)<fact))
1030 base2.query_polygon_internal(vertex,fact,pixset);
1033 query_polygon_internal(vertex, fact, pixset);
1037 (
double theta1,
double theta2,
bool inclusive, rangeset<I> &pixset)
const 1041 I ring1 = max(I(1),1+ring_above(cos(theta1))),
1042 ring2 = min(4*nside_-1,ring_above(cos(theta2)));
1045 ring1 = max(I(1),ring1-1);
1046 ring2 = min(4*nside_-1,ring2+1);
1051 get_ring_info_small(ring1,sp1,rp1,dummy);
1052 get_ring_info_small(ring2,sp2,rp2,dummy);
1055 if (pix1<=pix2) pixset.append(pix1,pix2);
1058 planck_fail(
"query_strip not yet implemented for NESTED");
1062 double theta2,
bool inclusive, rangeset<I> &pixset)
const 1067 query_strip_internal(theta1,theta2,inclusive,pixset);
1070 query_strip_internal(0.,theta2,inclusive,pixset);
1072 query_strip_internal(theta1,pi,inclusive,ps2);
1078 (I ring, I &startpix, I &ringpix,
bool &shifted)
const 1084 startpix = 2*ring*(ring-1);
1086 else if (ring < 3*nside_)
1088 shifted = ((ring-nside_) & 1) == 0;
1090 startpix = ncap_ + (ring-nside_)*ringpix;
1095 I nr= 4*nside_-ring;
1097 startpix = npix_-2*nr*(nr+1);
1102 I &ringpix,
double &costheta,
double &sintheta,
bool &shifted)
const 1104 I northring = (ring>2*nside_) ? 4*nside_-ring : ring;
1105 if (northring < nside_)
1107 double tmp = northring*northring*fact2_;
1109 sintheta = sqrt(tmp*(2-tmp));
1110 ringpix = 4*northring;
1112 startpix = 2*northring*(northring-1);
1116 costheta = (2*nside_-northring)*fact1_;
1117 sintheta = sqrt((1+costheta)*(1-costheta));
1119 shifted = ((northring-nside_) & 1) == 0;
1120 startpix = ncap_ + (northring-nside_)*ringpix;
1122 if (northring != ring)
1124 costheta = -costheta;
1125 startpix = npix_ - startpix - ringpix;
1129 I &startpix, I &ringpix,
double &theta,
bool &shifted)
const 1131 I northring = (ring>2*nside_) ? 4*nside_-ring : ring;
1132 if (northring < nside_)
1134 double tmp = northring*northring*fact2_;
1135 double costheta = 1 - tmp;
1136 double sintheta = sqrt(tmp*(2-tmp));
1137 theta = atan2(sintheta,costheta);
1138 ringpix = 4*northring;
1140 startpix = 2*northring*(northring-1);
1144 theta = acos((2*nside_-northring)*fact1_);
1146 shifted = ((northring-nside_) & 1) == 0;
1147 startpix = ncap_ + (northring-nside_)*ringpix;
1149 if (northring != ring)
1152 startpix = npix_ - startpix - ringpix;
1157 fix_arr<I,8> &result)
const 1159 int ix, iy, face_num;
1161 ring2xyf(pix,ix,iy,face_num) : nest2xyf(pix,ix,iy,face_num);
1163 const I nsm1 = nside_-1;
1164 if ((ix>0)&&(ix<nsm1)&&(iy>0)&&(iy<nsm1))
1167 for (
int m=0; m<8; ++m)
1168 result[m] = xyf2ring(ix+nb_xoffset[m],iy+nb_yoffset[m],face_num);
1171 I fpix = I(face_num)<<(2*order_),
1172 px0=spread_bits(ix ), py0=spread_bits(iy )<<1,
1173 pxp=spread_bits(ix+1), pyp=spread_bits(iy+1)<<1,
1174 pxm=spread_bits(ix-1), pym=spread_bits(iy-1)<<1;
1176 result[0] = fpix+pxm+py0; result[1] = fpix+pxm+pyp;
1177 result[2] = fpix+px0+pyp; result[3] = fpix+pxp+pyp;
1178 result[4] = fpix+pxp+py0; result[5] = fpix+pxp+pym;
1179 result[6] = fpix+px0+pym; result[7] = fpix+pxm+pym;
1184 for (
int i=0; i<8; ++i)
1186 int x=ix+nb_xoffset[i], y=iy+nb_yoffset[i];
1189 { x+=nside_; nbnum-=1; }
1191 { x-=nside_; nbnum+=1; }
1193 { y+=nside_; nbnum-=3; }
1195 { y-=nside_; nbnum+=3; }
1197 int f = nb_facearray[nbnum][face_num];
1200 int bits = nb_swaparray[nbnum][face_num>>2];
1201 if (bits&1) x=nside_-x-1;
1202 if (bits&2) y=nside_-y-1;
1203 if (bits&4) std::swap(x,y);
1204 result[i] = (scheme_==
RING) ? xyf2ring(x,y,f) : xyf2nest(x,y,f);
1213 fix_arr<I,4> &pix, fix_arr<double,4> &wgt)
const 1215 planck_assert((ptg.theta>=0)&&(ptg.theta<=pi),
"invalid theta value");
1216 double z = cos (ptg.theta);
1217 I ir1 = ring_above(z);
1219 double theta1, theta2, w1, tmp, dphi;
1225 get_ring_info2 (ir1, sp, nr, theta1, shift);
1227 tmp = (ptg.phi/dphi - .5*shift);
1228 i1 = (tmp<0) ? I(tmp)-1 : I(tmp);
1229 w1 = (ptg.phi-(i1+.5*shift)*dphi)/dphi;
1232 if (i2>=nr) i2 -=nr;
1233 pix[0] = sp+i1; pix[1] = sp+i2;
1234 wgt[0] = 1-w1; wgt[1] = w1;
1238 get_ring_info2 (ir2, sp, nr, theta2, shift);
1240 tmp = (ptg.phi/dphi - .5*shift);
1241 i1 = (tmp<0) ? I(tmp)-1 : I(tmp);
1242 w1 = (ptg.phi-(i1+.5*shift)*dphi)/dphi;
1245 if (i2>=nr) i2 -=nr;
1246 pix[2] = sp+i1; pix[3] = sp+i2;
1247 wgt[2] = 1-w1; wgt[3] = w1;
1252 double wtheta = ptg.theta/theta2;
1253 wgt[2] *= wtheta; wgt[3] *= wtheta;
1254 double fac = (1-wtheta)*0.25;
1255 wgt[0] = fac; wgt[1] = fac; wgt[2] += fac; wgt[3] +=fac;
1256 pix[0] = (pix[2]+2)&3;
1257 pix[1] = (pix[3]+2)&3;
1259 else if (ir2==4*nside_)
1261 double wtheta = (ptg.theta-theta1)/(pi-theta1);
1262 wgt[0] *= (1-wtheta); wgt[1] *= (1-wtheta);
1263 double fac = wtheta*0.25;
1264 wgt[0] += fac; wgt[1] += fac; wgt[2] = fac; wgt[3] =fac;
1265 pix[2] = ((pix[0]+2)&3)+npix_-4;
1266 pix[3] = ((pix[1]+2)&3)+npix_-4;
1270 double wtheta = (ptg.theta-theta1)/(theta2-theta1);
1271 wgt[0] *= (1-wtheta); wgt[1] *= (1-wtheta);
1272 wgt[2] *= wtheta; wgt[3] *= wtheta;
1276 for (tsize m=0; m<pix.size(); ++m)
1277 pix[m] = ring2nest(pix[m]);
1282 std::swap(order_,other.
order_);
1283 std::swap(nside_,other.
nside_);
1284 std::swap(npface_,other.npface_);
1285 std::swap(ncap_,other.ncap_);
1286 std::swap(npix_,other.npix_);
1287 std::swap(fact1_,other.fact1_);
1288 std::swap(fact2_,other.fact2_);
1289 std::swap(scheme_,other.
scheme_);
1295 va.set_z_phi (2./3., pi/(4*nside_));
1296 double t1 = 1.-1./nside_;
1298 vb.set_z_phi (1-t1/3, 0);
1299 return v_angle(va,vb);
1304 if (ring>=2*nside_) ring=4*nside_-ring;
1305 double z=ring2z(ring), z_up=ring2z(ring-1);
1307 uppos.set_z_phi(z_up,0);
1310 mypos.set_z_phi(z,pi/(4*ring));
1311 double v1=v_angle(mypos,uppos);
1312 if (ring!=1)
return v1;
1313 uppos.set_z_phi(ring2z(ring+1),pi/(4*(min(nside_,ring+1))));
1314 return max(v1,v_angle(mypos,uppos));
1316 mypos.set_z_phi(z,0);
1317 double vdist=v_angle(mypos,uppos);
1318 double hdist=sqrt(1.-z*z)*pi/(4*nside_);
1319 return max(hdist,vdist);
1323 int face,
double &z,
double &phi,
double &sth,
bool &have_sth)
const 1326 double jr = jrll[face] - x - y;
1331 double tmp = nr*nr/3.;
1335 sth = std::sqrt(tmp*(2.0-tmp));
1342 double tmp = nr*nr/3.;
1346 sth = std::sqrt(tmp*(2.-tmp));
1356 double tmp=jpll[face]*nr+x-y;
1359 phi = (nr<1e-15) ? 0 : (0.5*halfpi*tmp)/nr;
1364 vec3 locToVec3 (
double z,
double phi,
double sth,
bool have_sth)
1367 return vec3(sth*cos(phi),sth*sin(phi),z);
1371 res.set_z_phi (z, phi);
1379 vector<vec3> &out)
const 1383 pix2xyf(pix, ix, iy, face);
1384 double dc = 0.5 / nside_;
1385 double xc = (ix + 0.5)/nside_, yc = (iy + 0.5)/nside_;
1386 double d = 1.0/(step*nside_);
1387 for (tsize i=0; i<step; ++i)
1391 xyf2loc(xc+dc-i*d, yc+dc, face, z, phi, sth, have_sth);
1392 out[i] = locToVec3(z, phi, sth, have_sth);
1393 xyf2loc(xc-dc, yc+dc-i*d, face, z, phi, sth, have_sth);
1394 out[i+step] = locToVec3(z, phi, sth, have_sth);
1395 xyf2loc(xc-dc+i*d, yc-dc, face, z, phi, sth, have_sth);
1396 out[i+2*step] = locToVec3(z, phi, sth, have_sth);
1397 xyf2loc(xc+dc, yc-dc+i*d, face, z, phi, sth, have_sth);
1398 out[i+3*step] = locToVec3(z, phi, sth, have_sth);
1404 planck_assert(order_>=0,
"need hierarchical map");
1405 planck_assert(order_<=13,
"map too large");
1406 arr<int> result(swap_clen[order_]);
1408 for (
int m=0; m<order_;++m) ofs+=swap_clen[m];
1409 for (tsize m=0; m<result.size();++m) result[m]=swap_cycle[m+ofs];
I nest2peano(I pix) const
void query_disc(pointing ptg, double radius, rangeset< I > &pixset) const
void get_ring_info2(I ring, I &startpix, I &ringpix, double &theta, bool &shifted) const
void query_polygon_inclusive(const std::vector< pointing > &vertex, rangeset< I > &pixset, int fact=1) const
void query_polygon(const std::vector< pointing > &vertex, rangeset< I > &pixset) const
static int nside2order(I nside)
void get_interpol(const pointing &ptg, fix_arr< I, 4 > &pix, fix_arr< double, 4 > &wgt) const
void SetNside(I nside, Healpix_Ordering_Scheme scheme)
void swap(T_Healpix_Base &other)
void Set(int order, Healpix_Ordering_Scheme scheme)
void query_strip(double theta1, double theta2, bool inclusive, rangeset< I > &pixset) const
void pix2zphi(I pix, double &z, double &phi) const
void boundaries(I pix, tsize step, std::vector< vec3 > &out) const
static I npix2nside(I npix)
void get_ring_info(I ring, I &startpix, I &ringpix, double &costheta, double &sintheta, bool &shifted) const
double max_pixrad() const
I ring_above(double z) const
void query_disc_inclusive(pointing ptg, double radius, rangeset< I > &pixset, int fact=1) const
void get_ring_info_small(I ring, I &startpix, I &ringpix, bool &shifted) const
void neighbors(I pix, fix_arr< I, 8 > &result) const
I peano2nest(I pix) const
double ring2z(I ring) const
Healpix_Ordering_Scheme scheme_