1 /* 2 * This file is part of Healpix Java. 3 * 4 * This code is free software; you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation; either version 2 of the License, or 7 * (at your option) any later version. 8 * 9 * This code is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this code; if not, write to the Free Software 16 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 17 * 18 * For more information about HEALPix, see http://healpix.sourceforge.net 19 */ 20 21 package healpix.essentials; 22 23 import java.util.Arrays; 24 25 /** Basic opersations related to the HEALPix pixelisation. 26 This class is conceptually very similar the the Healpix_Base class 27 of Healpix_cxx. It supports the NESTED scheme for nside parameters which are 28 powers of 2, and the RING scheme for arbitrary positive nside parameters. 29 The maximum supported nside parameter is 2^29. 30 31 @copyright 2011, 2012 Max-Planck-Society 32 @author Martin Reinecke */ 33 public class HealpixBase extends HealpixTables 34 { 35 protected final class Xyf 36 { 37 public int ix, iy, face; Xyf()38 public Xyf () {} Xyf(int x, int y, int f)39 public Xyf (int x, int y, int f) 40 { ix=x; iy=y; face=f; } 41 } 42 43 private final class RingInfoSmall 44 { 45 long startpix, ringpix; 46 boolean shifted; 47 } 48 49 /** Maximum order parameter */ 50 public static final int order_max=29; 51 /** Maximum Nside parameter; equivalent to 2^{@code order_max}. */ 52 public static final long ns_max=1L<<order_max; 53 54 /** The order of the map; -1 when {@code nside} is not a power of 2. */ 55 protected int order; 56 57 /** The Nside parameter. */ 58 protected long nside; 59 60 protected long nl2, nl3, nl4, npface, npix, ncap; 61 protected double fact1, fact2; 62 63 /** The ordering scheme. */ 64 protected Scheme scheme; 65 spread_bits(int v)66 private static long spread_bits (int v) 67 { 68 return (long)(utab[ v &0xff]) | ((long)(utab[(v>>> 8)&0xff])<<16) 69 | ((long)(utab[(v>>>16)&0xff])<<32) | ((long)(utab[(v>>>24)&0xff])<<48); 70 } compress_bits(long v)71 private static int compress_bits (long v) 72 { 73 long raw = v&0x5555555555555555L; 74 raw |= raw>>>15; 75 int raw1 = (int)(raw&0xffffL), raw2 = (int)((raw>>>32)&0xffffL); 76 return ctab[raw1&0xff] | (ctab[raw1>>>8]<< 4) 77 | (ctab[raw2&0xff]<<16) | (ctab[raw2>>>8]<<20); 78 } 79 nest2xyf(long ipix)80 private Xyf nest2xyf(long ipix) 81 { 82 long pix=ipix&(npface-1); 83 return new Xyf (compress_bits(pix), compress_bits(pix>>>1), 84 (int)(ipix>>>(2*order))); 85 } 86 xyf2nest(int ix, int iy, int face_num)87 private long xyf2nest(int ix, int iy, int face_num) 88 { 89 return ((long)(face_num)<<(2*order)) + 90 spread_bits(ix) + (spread_bits(iy)<<1); 91 } 92 xyf2ring(int ix, int iy, int face_num)93 private long xyf2ring(int ix, int iy, int face_num) 94 { 95 long jr = ((long)jrll[face_num]*nside) - (long)ix - (long)iy - 1L; 96 97 RingInfoSmall ris = get_ring_info_small(jr); 98 long nr=ris.ringpix>>>2, kshift=ris.shifted ? 0:1; 99 100 long jp = (jpll[face_num]*nr + (long)ix - (long)iy + 1L + kshift) / 2L; 101 if (jp>nl4) 102 jp-=nl4; 103 else 104 if (jp<1) jp+=nl4; 105 106 return ris.startpix + jp - 1L; 107 } 108 ring2xyf(long pix)109 protected Xyf ring2xyf(long pix) 110 { 111 Xyf ret = new Xyf(); 112 long iring, iphi, kshift, nr; 113 114 if (pix<ncap) // North Polar cap 115 { 116 iring = (1+HealpixUtils.isqrt(1L+2L*pix))>>>1; //counted from North pole 117 iphi = (pix+1) - 2*iring*(iring-1); 118 kshift = 0; 119 nr = iring; 120 ret.face=(int)((iphi-1)/nr); 121 } 122 else if (pix<(npix-ncap)) // Equatorial region 123 { 124 long ip = pix - ncap; 125 long tmp = (order>=0) ? ip>>>(order+2) : ip/nl4; 126 iring = tmp+nside; 127 iphi = ip-tmp*nl4 + 1; 128 kshift = (iring+nside)&1; 129 nr = nside; 130 long ire = iring-nside+1, 131 irm = nl2+2-ire; 132 long ifm = iphi - (ire>>>1) + nside -1, 133 ifp = iphi - (irm>>>1) + nside -1; 134 if (order>=0) 135 { ifm >>>= order; ifp >>>= order; } 136 else 137 { ifm /= nside; ifp /= nside; } 138 ret.face = (int)((ifp==ifm) ? (ifp|4) : ((ifp<ifm) ? ifp : (ifm+8))); 139 } 140 else // South Polar cap 141 { 142 long ip = npix - pix; 143 iring = (1+HealpixUtils.isqrt(2L*ip-1L))>>>1; //counted from South pole 144 iphi = 4L*iring + 1 - (ip - 2L*iring*(iring-1L)); 145 kshift = 0; 146 nr = iring; 147 iring = 2L*nl2-iring; 148 ret.face = 8 + (int)((iphi-1)/nr); 149 } 150 151 long irt = iring - ((long)jrll[ret.face]*nside) + 1L; 152 long ipt = 2L*iphi- (long)jpll[ret.face]*nr - kshift -1L; 153 if (ipt>=nl2) ipt-=8L*nside; 154 155 ret.ix = (int)(( ipt-irt)>>>1); 156 ret.iy = (int)((-ipt-irt)>>>1); 157 158 return ret; 159 } 160 pix2xyf(long pix)161 protected Xyf pix2xyf (long pix) 162 { return (scheme==Scheme.RING) ? ring2xyf(pix) : nest2xyf(pix); } 163 xyf2pix(int ix, int iy, int face_num)164 protected long xyf2pix (int ix, int iy, int face_num) 165 { 166 return (scheme==Scheme.RING) ? 167 xyf2ring(ix,iy,face_num) : xyf2nest(ix,iy,face_num); 168 } xyf2pix(Xyf xyf)169 protected long xyf2pix (Xyf xyf) 170 { 171 return (scheme==Scheme.RING) ? 172 xyf2ring(xyf.ix,xyf.iy,xyf.face) : xyf2nest(xyf.ix,xyf.iy,xyf.face); 173 } 174 175 /** Calculates the map order from its Nside parameter. 176 @param nside the Nside parameter 177 @return the map order corresponding to {@code nside}; -1 if 178 {@code nside} is not a power of 2. */ nside2order(long nside)179 public static int nside2order(long nside) throws Exception 180 { 181 HealpixUtils.check (nside>0,"nside must be positive"); 182 return ((nside&(nside-1))!=0) ? -1 : HealpixUtils.ilog2(nside); 183 } 184 185 /** Calculates the Nside parameter from the number of pixels. 186 @param npix the number of pixels 187 @return the computed Nside parameter */ npix2Nside(long npix)188 public static long npix2Nside(long npix) throws Exception 189 { 190 HealpixUtils.check (npix>0,"npix must be positive"); 191 long nside = HealpixUtils.isqrt(npix/12); 192 HealpixUtils.check(12*nside*nside==npix,"npix is not 12*nside*nside"); 193 HealpixUtils.check(nside<=ns_max,"nside is too large"); 194 return nside; 195 } 196 /** Calculates the number of pixels from the Nside parameter. 197 @param nside the Nside parameter 198 @return the computed number of pixels */ nside2Npix(long nside)199 public static long nside2Npix(long nside) throws Exception 200 { 201 HealpixUtils.check (nside>0,"nside must be positive"); 202 HealpixUtils.check(nside<=ns_max,"nside is too large"); 203 return 12*nside*nside; 204 } 205 206 /** Calculates the number of pixels from the order parameter. 207 @param order the order parameter 208 @return the computed number of pixels */ order2Npix(int order)209 public static long order2Npix(int order) throws Exception 210 { 211 HealpixUtils.check (order>=0,"order must be nonnegative"); 212 HealpixUtils.check(order<=order_max,"order is too large"); 213 return 12*(1L<<(2*order)); 214 } 215 216 /** Adjusts the object to nside_in. 217 @param nside_in the new Nside parameter */ setNside(long nside_in)218 public void setNside (long nside_in) throws Exception 219 { 220 if (nside==nside_in) return; 221 nside=nside_in; 222 HealpixUtils.check(nside<=ns_max && nside>0, 223 "Nside must be between 1 and " + ns_max); 224 225 order = nside2order(nside); 226 if (scheme==Scheme.NESTED && order<0) 227 throw new Exception("Nside must be a power of 2 for NESTED scheme"); 228 nl2 = 2*nside; 229 nl3 = 3*nside; 230 nl4 = 4*nside; 231 npface = nside*nside; 232 ncap = 2*nside*(nside-1); // pixels in each polar cap 233 npix = 12*npface; 234 fact2 = 4.0/npix; 235 fact1 = (nside<<1)*fact2; 236 } 237 238 /** Adjusts the object to nside_in and scheme_in. 239 @param nside_in the new Nside parameter 240 @param scheme_in the new ordering scheme */ setNsideAndScheme(long nside_in, Scheme scheme_in)241 public void setNsideAndScheme (long nside_in, Scheme scheme_in) 242 throws Exception 243 { 244 if ((scheme==scheme_in) && (nside==nside_in)) return; 245 HealpixUtils.check (((nside_in&(nside_in-1))==0)||(scheme_in==Scheme.RING), 246 "Nside must be a power of 2 for NESTED scheme"); 247 scheme=scheme_in; 248 setNside(nside_in); 249 } 250 251 /** Initializes the object to Nside=1 and scheme=NESTED. */ HealpixBase()252 public HealpixBase() 253 { 254 try 255 { nside=0; setNsideAndScheme(1,Scheme.NESTED); } 256 catch(Exception Ex) {/*cannot happen*/} 257 } 258 259 /** Initializes the object to a user-supplied Nside and ordering scheme. 260 @param nside_in the Nside parameter 261 @param scheme_in the ordering scheme */ HealpixBase(long nside_in, Scheme scheme_in)262 public HealpixBase(long nside_in, Scheme scheme_in) throws Exception 263 { nside=nside_in-1; setNsideAndScheme(nside_in, scheme_in); } 264 265 /** Returns the current ordering scheme. 266 @return the current ordering scheme */ getScheme()267 public Scheme getScheme() 268 { return scheme; } 269 270 /** Returns the current Nside parameter. 271 @return the current Nside parameter */ getNside()272 public int getNside() 273 { return (int)nside; } 274 275 /** Returns the total number of pixels in the pixelisation. 276 @return the current total number of pixels */ getNpix()277 public long getNpix() 278 { return npix; } 279 280 /** Adjusts the object to scheme. 281 @param scheme_in the new ordering scheme */ setScheme(Scheme scheme_in)282 public void setScheme(Scheme scheme_in) throws Exception 283 { 284 if (scheme_in==Scheme.NESTED && order<0) 285 throw new Exception("Nside must be a power of 2 for NESTED scheme"); 286 scheme = scheme_in; 287 } 288 289 /** Returns the current order parameter. 290 @return the current map order parameter. */ getOrder()291 public int getOrder() 292 { return order; } 293 294 /** Returns the pixel which contains the supplied Pointing. 295 @param ptg the requested location on the sphere. 296 @return the pixel number containing the location. */ ang2pix(Pointing ptg)297 public long ang2pix(Pointing ptg) throws Exception 298 { return loc2pix(new Hploc(ptg)); } 299 300 /** Returns the Pointing corresponding to the center of the supplied pixel. 301 @param pix the requested pixel number. 302 @return the pixel's center coordinates. */ pix2ang(long pix)303 public Pointing pix2ang(long pix) throws Exception 304 { return pix2loc(pix).toPointing(); } 305 306 /** Returns the pixel which contains the supplied Vec3. 307 @param vec the requested location on the sphere (need not be normalized). 308 @return the pixel number containing the location. */ vec2pix(Vec3 vec)309 public long vec2pix(Vec3 vec) throws Exception 310 { return loc2pix(new Hploc(vec)); } 311 312 /** Returns the normalized 3-vector corresponding to the center of the 313 supplied pixel. 314 @param pix the requested pixel number. 315 @return the pixel's center coordinates. */ pix2vec(long pix)316 public Vec3 pix2vec(long pix) throws Exception 317 { return pix2loc(pix).toVec3(); } 318 319 /** Returns nested pixel number for the supplied ring pixel number. 320 @param ipring the requested pixel number in RING scheme. 321 @return the corresponding pixel number in NESTED scheme. */ ring2nest(long ipring)322 public long ring2nest(long ipring) throws Exception 323 { 324 Xyf xyf = ring2xyf(ipring); 325 return xyf2nest (xyf.ix,xyf.iy, xyf.face); 326 } 327 328 /** Returns ring pixel number for the supplied nested pixel number. 329 @param ipnest the requested pixel number in NESTED scheme. 330 @return the corresponding pixel number in RING scheme. */ nest2ring(long ipnest)331 public long nest2ring(long ipnest) throws Exception 332 { 333 Xyf xyf = nest2xyf(ipnest); 334 return xyf2ring (xyf.ix,xyf.iy, xyf.face); 335 } 336 loc2pix(Hploc loc)337 protected long loc2pix (Hploc loc) 338 { 339 double z=loc.z, phi=loc.phi; 340 341 double za = Math.abs(z); 342 double tt = HealpixUtils.fmodulo((phi*Constants.inv_halfpi),4.0);// in [0,4) 343 344 if (scheme==Scheme.RING) 345 { 346 if (za<=Constants.twothird) // Equatorial region 347 { 348 double temp1 = nside*(0.5+tt); 349 double temp2 = nside*z*0.75; 350 long jp = (long)(temp1-temp2); // index of ascending edge line 351 long jm = (long)(temp1+temp2); // index of descending edge line 352 353 // ring number counted from z=2/3 354 long ir = nside + 1 + jp - jm; // in {1,2n+1} 355 long kshift = 1-(ir&1); // kshift=1 if ir even, 0 otherwise 356 357 long t1 = jp+jm-nside+kshift+1+nl4+nl4; 358 long ip = (order>0) ? 359 (t1>>>1)&(nl4-1) : ((t1>>>1)%nl4); // in {0,4n-1} 360 361 return ncap + (ir-1)*nl4 + ip; 362 } 363 else // North & South polar caps 364 { 365 double tp = tt-(long)(tt); 366 double tmp = ((za<0.99)||(!loc.have_sth)) ? 367 nside*Math.sqrt(3*(1-za)) : 368 nside*loc.sth/Math.sqrt((1.+za)/3.); 369 370 long jp = (long)(tp*tmp); // increasing edge line index 371 long jm = (long)((1.0-tp)*tmp); // decreasing edge line index 372 373 long ir = jp+jm+1; // ring number counted from the closest pole 374 long ip = (long)(tt*ir); // in {0,4*ir-1} 375 assert((ip>=0)&&(ip<4*ir)); 376 377 return (z>0) ? 2*ir*(ir-1) + ip : npix - 2*ir*(ir+1) + ip; 378 } 379 } 380 else // scheme_ == NEST 381 { 382 if (za<=Constants.twothird) // Equatorial region 383 { 384 double temp1 = nside*(0.5+tt); 385 double temp2 = nside*(z*0.75); 386 long jp = (long)(temp1-temp2); // index of ascending edge line 387 long jm = (long)(temp1+temp2); // index of descending edge line 388 long ifp = jp >>> order; // in {0,4} 389 long ifm = jm >>> order; 390 long face_num = (ifp==ifm) ? (ifp|4) : ((ifp<ifm) ? ifp : (ifm+8)); 391 392 long ix = jm & (nside-1), 393 iy = nside - (jp & (nside-1)) - 1; 394 return xyf2nest((int)ix,(int)iy,(int)face_num); 395 } 396 else // polar region, za > 2/3 397 { 398 int ntt = Math.min(3,(int)tt); 399 double tp = tt-ntt; 400 double tmp = ((za<0.99)||(!loc.have_sth)) ? 401 nside*Math.sqrt(3*(1-za)) : 402 nside*loc.sth/Math.sqrt((1.+za)/3.); 403 404 long jp = (long)(tp*tmp); // increasing edge line index 405 long jm = (long)((1.0-tp)*tmp); // decreasing edge line index 406 if (jp>=nside) jp = nside-1; // for points too close to the boundary 407 if (jm>=nside) jm = nside-1; 408 return (z>=0) ? 409 xyf2nest((int)(nside-jm -1),(int)(nside-jp-1),ntt) : 410 xyf2nest((int)jp,(int)jm,ntt+8); 411 } 412 } 413 } 414 415 /** Returns the pixel which contains the supplied Zphi. 416 @param zphi the requested location on the sphere. 417 @return the pixel number containing the location. */ zphi2pix(Zphi zphi)418 public long zphi2pix (Zphi zphi) 419 { return loc2pix (new Hploc(zphi)); } 420 pix2loc(long pix)421 protected Hploc pix2loc (long pix) 422 { 423 Hploc loc = new Hploc(); 424 if (scheme==Scheme.RING) 425 { 426 if (pix<ncap) // North Polar cap 427 { 428 long iring = (1+(HealpixUtils.isqrt(1+2*pix)))>>>1; //from North pole 429 long iphi = (pix+1) - 2*iring*(iring-1); 430 431 double tmp = (iring*iring)*fact2; 432 loc.z = 1.0 - tmp; 433 if (loc.z>0.99) { loc.sth=Math.sqrt(tmp*(2.-tmp)); loc.have_sth=true; } 434 loc.phi = (iphi-0.5) * Constants.halfpi/iring; 435 } 436 else if (pix<(npix-ncap)) // Equatorial region 437 { 438 long ip = pix - ncap; 439 long tmp = (order>=0) ? ip>>>(order+2) : ip/nl4; 440 long iring = tmp + nside, 441 iphi = ip-nl4*tmp+1;; 442 // 1 if iring+nside is odd, 1/2 otherwise 443 double fodd = ((iring+nside)&1)!=0 ? 1 : 0.5; 444 445 loc.z = (nl2-iring)*fact1; 446 loc.phi = (iphi-fodd) * Math.PI*0.75*fact1; 447 } 448 else // South Polar cap 449 { 450 long ip = npix - pix; 451 long iring = (1+HealpixUtils.isqrt(2*ip-1))>>>1; //from South pole 452 long iphi = 4*iring + 1 - (ip - 2*iring*(iring-1)); 453 454 double tmp = (iring*iring)*fact2; 455 loc.z = tmp-1.0; 456 if (loc.z<-0.99) { loc.sth=Math.sqrt(tmp*(2.-tmp)); loc.have_sth=true; } 457 loc.phi = (iphi-0.5) * Constants.halfpi/iring; 458 } 459 } 460 else 461 { 462 Xyf xyf= nest2xyf(pix); 463 464 long jr = ((long)(jrll[xyf.face])<<order) -xyf.ix - xyf.iy - 1; 465 466 long nr; 467 if (jr<nside) 468 { 469 nr = jr; 470 double tmp = (nr*nr)*fact2; 471 loc.z = 1 - tmp; 472 if (loc.z>0.99) { loc.sth=Math.sqrt(tmp*(2.-tmp)); loc.have_sth=true; } 473 } 474 else if (jr>nl3) 475 { 476 nr = nl4-jr; 477 double tmp = (nr*nr)*fact2; 478 loc.z = tmp - 1; 479 if (loc.z<-0.99) { loc.sth=Math.sqrt(tmp*(2.-tmp)); loc.have_sth=true; } 480 } 481 else 482 { 483 nr = nside; 484 loc.z = (nl2-jr)*fact1; 485 } 486 487 long tmp=(long)(jpll[xyf.face])*nr+xyf.ix-xyf.iy; 488 assert(tmp<8*nr); // must not happen 489 if (tmp<0) tmp+=8*nr; 490 loc.phi = (nr==nside) ? 0.75*Constants.halfpi*tmp*fact1 : 491 (0.5*Constants.halfpi*tmp)/nr; 492 } 493 return loc; 494 } 495 496 /** Returns the Zphi corresponding to the center of the supplied pixel. 497 @param pix the requested pixel number. 498 @return the pixel's center coordinates. */ 499 public Zphi pix2zphi (long pix) 500 { return pix2loc(pix).toZphi(); } 501 502 /** Returns the neighboring pixels of ipix. 503 This method works in both RING and NEST schemes, but is 504 considerably faster in the NEST scheme. 505 @param ipix the requested pixel number. 506 @return array with indices of the neighboring pixels. 507 The returned array contains (in this order) 508 the pixel numbers of the SW, W, NW, N, NE, E, SE and S neighbor 509 of ipix. If a neighbor does not exist (this can only happen 510 for the W, N, E and S neighbors), its entry is set to -1. */ 511 public long[] neighbours(long ipix) throws Exception 512 { 513 long[] result = new long[8]; 514 Xyf xyf = (scheme==Scheme.NESTED) ? nest2xyf(ipix) : ring2xyf(ipix); 515 int ix = xyf.ix, iy=xyf.iy, face_num=xyf.face; 516 517 long nsm1 = nside-1; 518 if ((ix>0)&&(ix<nsm1)&&(iy>0)&&(iy<nsm1)) 519 { 520 if (scheme==Scheme.RING) 521 for (int m=0; m<8; ++m) 522 result[m]=(xyf2ring(ix+xoffset[m],iy+yoffset[m],face_num)); 523 else 524 { 525 long fpix = (long)(face_num)<<(2*order), 526 px0=spread_bits(ix ), py0=spread_bits(iy )<<1, 527 pxp=spread_bits(ix+1), pyp=spread_bits(iy+1)<<1, 528 pxm=spread_bits(ix-1), pym=spread_bits(iy-1)<<1; 529 530 result[0]=fpix+pxm+py0; result[1]=fpix+pxm+pyp; 531 result[2]=fpix+px0+pyp; result[3]=fpix+pxp+pyp; 532 result[4]=fpix+pxp+py0; result[5]=fpix+pxp+pym; 533 result[6]=fpix+px0+pym; result[7]=fpix+pxm+pym; 534 } 535 } 536 else 537 { 538 for (int i=0; i<8; ++i) 539 { 540 int x=ix+xoffset[i]; 541 int y=iy+yoffset[i]; 542 int nbnum=4; 543 if (x<0) 544 { x+=nside; nbnum-=1; } 545 else if (x>=nside) 546 { x-=nside; nbnum+=1; } 547 if (y<0) 548 { y+=nside; nbnum-=3; } 549 else if (y>=nside) 550 { y-=nside; nbnum+=3; } 551 552 int f = facearray[nbnum][face_num]; 553 554 if (f>=0) 555 { 556 int bits = swaparray[nbnum][face_num>>>2]; 557 if ((bits&1)>0) x=(int)(nside-x-1); 558 if ((bits&2)>0) y=(int)(nside-y-1); 559 if ((bits&4)>0) { int tint=x; x=y; y=tint; } 560 result[i] = (scheme==Scheme.NESTED) ? 561 xyf2nest(x,y,f) : xyf2ring(x,y,f); 562 } 563 else 564 result[i]=-1; 565 } 566 } 567 return result; 568 } 569 570 /** Returns the maximum angular distance between a pixel center and its 571 corners. 572 @return maximum angular distance between a pixel center and its 573 corners. */ 574 public double maxPixrad() 575 { 576 Vec3 va = new Vec3(new Zphi(2./3., Math.PI/nl4)); 577 double t1 = 1.-1./nside; 578 t1*=t1; 579 Vec3 vb = new Vec3(new Zphi(1-t1/3, 0)); 580 return va.angle(vb); 581 } 582 583 private RingInfoSmall get_ring_info_small (long ring) 584 { 585 RingInfoSmall ret = new RingInfoSmall(); 586 if (ring<nside) 587 { 588 ret.shifted = true; 589 ret.ringpix = 4*ring; 590 ret.startpix = 2*ring*(ring-1); 591 } 592 else if (ring<nl3) 593 { 594 ret.shifted = ((ring-nside) & 1) == 0; 595 ret.ringpix = nl4; 596 ret.startpix = ncap + (ring-nside)*nl4; 597 } 598 else 599 { 600 ret.shifted = true; 601 long nr=nl4-ring; 602 ret.ringpix = 4*nr; 603 ret.startpix = npix-2*nr*(nr+1); 604 } 605 return ret; 606 } 607 608 private long ringAbove (double z) 609 { 610 double az=Math.abs(z); 611 if (az> Constants.twothird) // polar caps 612 { 613 long iring = (long)(nside*Math.sqrt(3*(1-az))); 614 return (z>0) ? iring : nl4-iring-1; 615 } 616 else // equatorial region 617 return (long)(nside*(2.0-1.5*z)); 618 } 619 620 /** Returns the z-coordinate (equivalent to the cosine of the colatitude) 621 for the requested ring. 622 This method also accepts the not-really-existing ring indices 0 and 623 4*Nside, which correspond to North and South Poles, respectively. 624 @param ring ring index: the northernmost ring in the map has index 1; 625 ring indices are increasing towards the South pole. 626 @return z-coordinate of the ring. */ 627 public double ring2z (long ring) 628 { 629 if (ring<nside) 630 return 1 - ring*ring*fact2; 631 if (ring <=nl3) 632 return (nl2-ring)*fact1; 633 ring=nl4 - ring; 634 return ring*ring*fact2 - 1; 635 } 636 /** Returns the colatitude for the requested ring. 637 This method also accepts the not-really-existing ring indices 0 and 638 4*Nside, which correspond to North and South Poles, respectively. 639 @param ring ring index: the northernmost ring in the map has index 1; 640 ring indices are increasing towards the South pole. 641 @return z-coordinate of the ring. */ 642 public double ring2theta (long ring) 643 { 644 if (ring<nside) 645 { 646 double tmp=ring*ring*fact2; 647 return FastMath.atan2 (Math.sqrt(tmp*(2.-tmp)),1-tmp); 648 } 649 if (ring <=nl3) 650 { 651 double cth=(nl2-ring)*fact1; 652 return FastMath.atan2 (Math.sqrt(1.-cth*cth),cth); 653 } 654 ring=nl4 - ring; 655 double tmp=ring*ring*fact2; 656 return FastMath.atan2 (Math.sqrt(tmp*(2.-tmp)),tmp-1); 657 } 658 659 private final class pstack 660 { 661 private long[] p; 662 private int[] o; 663 private int s, m; 664 665 public pstack (int sz) 666 { p=new long[sz]; o=new int[sz]; s=m=0; } 667 public void push (long p_, int o_) 668 { p[s]=p_; o[s]=o_; ++s; } 669 public void pop () 670 { --s; } 671 public void popToMark () 672 { s=m; } 673 public int size () 674 { return s; } 675 public void mark () 676 { m=s; } 677 public int otop () 678 { return o[s-1]; } 679 public long ptop () 680 { return p[s-1]; } 681 } 682 683 private RangeSet queryStripInternal(double theta1, double theta2, 684 boolean inclusive) throws Exception 685 { 686 RangeSet pixset = new RangeSet(1); 687 if (scheme==Scheme.RING) 688 { 689 long ring1 = Math.max(1,1+ringAbove(FastMath.cos(theta1))), 690 ring2 = Math.min(4*nside-1,ringAbove(FastMath.cos(theta2))); 691 if (inclusive) 692 { 693 ring1 = Math.max(1,ring1-1); 694 ring2 = Math.min(4*nside-1,ring2+1); 695 } 696 697 RingInfoSmall ris1 = get_ring_info_small(ring1), 698 ris2 = get_ring_info_small(ring2); 699 long pix1 = ris1.startpix, 700 pix2 = ris2.startpix+ris2.ringpix; 701 if (pix1<pix2) pixset.append(pix1,pix2); 702 } 703 else 704 HealpixUtils.check(false,"queryStrip not yet implemented for NESTED"); 705 706 return pixset; 707 } 708 709 /** Returns a range set of pixels whose centers lie within a given latitude 710 range (if {@code inclusive==false}), or which overlap with this range 711 (if {@code inclusive==true}).<p> 712 The latitude range is defined as follows: 713 <ul> 714 <li>if {@code theta1<theta2}, it is the range between {@code theta1} 715 and {@code theta2}</li> 716 <li>otherwise it is the range between 0 and {@code theta2}, and between 717 {@code theta1} and pi.</li> 718 </ul> 719 This method currently only works in the RING scheme. */ 720 public RangeSet queryStrip(double theta1, double theta2, 721 boolean inclusive) throws Exception 722 { 723 if (theta1<theta2) 724 return queryStripInternal(theta1,theta2,inclusive); 725 RangeSet res = queryStripInternal(0.,theta2,inclusive); 726 return res.union(queryStripInternal(theta1,Math.PI,inclusive)); 727 } 728 729 private boolean checkPixelRing (HealpixBase b2, long pix, long nr, 730 long ipix1, int fct, Zphi czphi, double cosrp2, long cpix) 731 { 732 if (pix>=nr) pix-=nr; 733 if (pix<0) pix+=nr; 734 pix+=ipix1; 735 if (pix==cpix) return false; // disk center in pixel => overlap 736 Xyf xyf=pix2xyf(pix); 737 for (int i=0; i<fct-1; ++i) // go along the 4 edges 738 { 739 int ox=fct*xyf.ix, oy=fct*xyf.iy; 740 Xyf xyf2 = new Xyf (ox,oy,xyf.face); 741 double pz,pphi; 742 xyf2.ix=ox+i; xyf2.iy=oy; 743 if (HealpixUtils.cosdist_zphi(czphi,b2.pix2zphi(b2.xyf2pix(xyf2)))>cosrp2) 744 return false; // overlap 745 xyf2.ix=ox+fct-1; xyf2.iy=oy+i; 746 if (HealpixUtils.cosdist_zphi(czphi,b2.pix2zphi(b2.xyf2pix(xyf2)))>cosrp2) 747 return false; // overlap 748 xyf2.ix=ox+fct-1-i; xyf2.iy=oy+fct-1; 749 if (HealpixUtils.cosdist_zphi(czphi,b2.pix2zphi(b2.xyf2pix(xyf2)))>cosrp2) 750 return false; // overlap 751 xyf2.ix=ox; xyf2.iy=oy+fct-1-i; 752 if (HealpixUtils.cosdist_zphi(czphi,b2.pix2zphi(b2.xyf2pix(xyf2)))>cosrp2) 753 return false; // overlap 754 } 755 return true; 756 } 757 758 private RangeSet queryDiscInternal(Pointing ptg, double radius, int fact) 759 throws Exception 760 { 761 boolean inclusive=(fact!=0); 762 RangeSet pixset = new RangeSet(); 763 if (scheme==Scheme.RING) 764 { 765 int fct=1; 766 if (inclusive) 767 { 768 HealpixUtils.check ((1L<<order_max)/nside>=fact, 769 "invalid oversampling factor"); 770 fct = fact; 771 } 772 HealpixBase b2=null; 773 double rsmall,rbig; 774 if (fct>1) 775 { 776 b2=new HealpixBase(fct*nside,Scheme.RING); 777 rsmall = radius+b2.maxPixrad(); 778 rbig = radius+maxPixrad(); 779 } 780 else 781 rsmall = rbig = inclusive ? radius+maxPixrad() : radius; 782 783 if (rsmall>=Math.PI) 784 { pixset.append(0,npix); return pixset; } 785 786 rbig = Math.min (Math.PI,rbig); 787 788 double cosrsmall = FastMath.cos(rsmall); 789 double cosrbig = FastMath.cos(rbig); 790 791 double z0 = FastMath.cos(ptg.theta); 792 double xa = 1./Math.sqrt((1-z0)*(1+z0)); 793 794 Zphi czphi = new Zphi(z0,ptg.phi); 795 long cpix = zphi2pix(czphi); 796 797 double rlat1 = ptg.theta - rsmall; 798 double zmax = FastMath.cos(rlat1); 799 long irmin = ringAbove (zmax)+1; 800 801 if ((rlat1<=0) && (irmin>1)) // north pole in the disk 802 { 803 RingInfoSmall info =get_ring_info_small(irmin-1); 804 pixset.append(0,info.startpix+info.ringpix); 805 } 806 807 if ((fct>1) && (rlat1>0)) irmin=Math.max(1L,irmin-1); 808 809 double rlat2 = ptg.theta + rsmall; 810 double zmin = FastMath.cos(rlat2); 811 long irmax = ringAbove (zmin); 812 813 if ((fct>1) && (rlat2<Math.PI)) irmax=Math.min(nl4-1,irmax+1); 814 815 for (long iz=irmin; iz<=irmax; ++iz) // rings partially in the disk 816 { 817 double z=ring2z(iz); 818 819 double x = (cosrbig-z*z0)*xa; 820 double ysq = 1-z*z-x*x; 821 double dphi=-1; 822 if (ysq<=0) // no intersection, ring completely inside or outside 823 dphi = (fct==1) ? 0: Math.PI-1e-15; 824 else 825 dphi = FastMath.atan2(Math.sqrt(ysq),x); 826 if (dphi>0) 827 { 828 RingInfoSmall info =get_ring_info_small(iz); 829 long ipix1 = info.startpix, nr=info.ringpix, ipix2=ipix1+nr-1; 830 double shift = info.shifted ? 0.5 : 0.; 831 832 long ip_lo = (long)Math.floor 833 (nr*Constants.inv_twopi*(ptg.phi-dphi) - shift)+1; 834 long ip_hi = (long)Math.floor 835 (nr*Constants.inv_twopi*(ptg.phi+dphi) - shift); 836 837 if (fct>1) 838 { 839 while ((ip_lo<=ip_hi) && checkPixelRing 840 (b2,ip_lo,nr,ipix1,fct,czphi,cosrsmall,cpix)) 841 ++ip_lo; 842 while ((ip_hi>ip_lo) && checkPixelRing 843 (b2,ip_hi,nr,ipix1,fct,czphi,cosrsmall,cpix)) 844 --ip_hi; 845 } 846 847 if (ip_lo<=ip_hi) 848 { 849 if (ip_hi>=nr) 850 { ip_lo-=nr; ip_hi-=nr; } 851 if (ip_lo<0) 852 { 853 pixset.append(ipix1,ipix1+ip_hi+1); 854 pixset.append(ipix1+ip_lo+nr,ipix2+1); 855 } 856 else 857 pixset.append(ipix1+ip_lo,ipix1+ip_hi+1); 858 } 859 } 860 } 861 862 if ((rlat2>=Math.PI) && (irmax+1<nl4)) // south pole in the disk 863 { 864 RingInfoSmall info =get_ring_info_small(irmax+1); 865 pixset.append(info.startpix,npix); 866 } 867 } 868 else // scheme_==NEST 869 { 870 if (radius>=Math.PI) // disk covers the whole sphere 871 { pixset.append(0,npix); return pixset; } 872 873 int oplus=0; 874 if (inclusive) 875 { 876 HealpixUtils.check ((1L<<order_max)>=fact,"invalid oversampling factor"); 877 HealpixUtils.check ((fact&(fact-1))==0, 878 "oversampling factor must be a power of 2"); 879 oplus=HealpixUtils.ilog2(fact); 880 } 881 int omax=Math.min(order_max,order+oplus); // the order up to which we test 882 883 Vec3 vptg = new Vec3(ptg); 884 885 double[] crpdr = new double[omax+1]; 886 double[] crmdr = new double[omax+1]; 887 888 double cosrad=FastMath.cos(radius), 889 sinrad=FastMath.sin(radius); 890 for (int o=0; o<=omax; o++) // prepare data at the required orders 891 { 892 double dr = HealpixProc.mpr[o]; // safety distance 893 double cdr = HealpixProc.cmpr[o]; 894 double sdr = HealpixProc.smpr[o]; 895 crpdr[o] = (radius+dr>Math.PI) ? -1. : cosrad*cdr-sinrad*sdr; 896 crmdr[o] = (radius-dr<0.) ? 1. : cosrad*cdr+sinrad*sdr; 897 } 898 899 pstack stk=new pstack(12+3*omax); 900 901 /* Still experimental, therefore disabled 902 Fxyf fxyf=new Fxyf(vptg); 903 for (int o=order; o>=0; --o) 904 { 905 long nsd=HealpixProc.bn[o].nside; 906 double fx=nsd*fxyf.fx-(int)(nsd*fxyf.fx), 907 fy=nsd*fxyf.fy-(int)(nsd*fxyf.fy); 908 double fmin = Math.min(Math.min(fx,fy),Math.min(1-fx,1-fy)); 909 if (fmin*0.7>nsd*radius) 910 { 911 System.out.println("contained completely at order "+o); 912 stk.push(HealpixProc.bn[o].vec2pix(vptg),o); 913 break; 914 } 915 } 916 if (stk.size()==0) 917 */ 918 919 for (int i=0; i<12; i++) // insert the 12 base pixels in reverse order 920 stk.push(11-i,0); 921 922 while (stk.size()>0) {// as long as there are pixels on the stack 923 // pop current pixel number and order from the stack 924 long pix=stk.ptop(); 925 int o=stk.otop(); 926 stk.pop(); 927 928 Zphi pos=HealpixProc.bn[o].pix2zphi(pix); 929 // cosine of angular distance between pixel center and disk center 930 double cangdist=HealpixUtils.cosdist_zphi(vptg.z,ptg.phi,pos.z,pos.phi); 931 932 if (cangdist>crpdr[o]) 933 { 934 int zone = (cangdist<cosrad) ? 1 : ((cangdist<=crmdr[o]) ? 2 : 3); 935 936 check_pixel (o, omax, zone, pixset, pix, stk, inclusive); 937 } 938 } 939 } 940 return pixset; 941 } 942 943 /** Returns a range set of pixels whose centers lie within a given disk. <p> 944 This method is more efficient in the RING scheme. 945 @param ptg the angular coordinates of the disk center 946 @param radius the radius (in radians) of the disk 947 @return the requested set of pixel number ranges */ 948 public RangeSet queryDisc(Pointing ptg, double radius) 949 throws Exception 950 { return queryDiscInternal (ptg, radius, 0); } 951 /** Returns a range set of pixels which overlap with a given disk. <p> 952 This method is more efficient in the RING scheme. <p> 953 This method may return some pixels which don't overlap with 954 the polygon at all. The higher {@code fact} is chosen, the fewer false 955 positives are returned, at the cost of increased run time. 956 @param ptg the angular coordinates of the disk center 957 @param radius the radius (in radians) of the disk 958 @param fact The overlapping test will be done at the resolution 959 {@code fact*nside}. For NESTED ordering, {@code fact} must be a power 960 of 2, else it can be any positive integer. A typical choice would be 4. 961 @return the requested set of pixel number ranges */ 962 public RangeSet queryDiscInclusive (Pointing ptg, double radius, int fact) 963 throws Exception 964 { return queryDiscInternal (ptg, radius, fact); } 965 966 private RangeSet queryMultiDisc (Vec3[] norm, double[] rad, 967 int fact) throws Exception 968 { 969 boolean inclusive = (fact!=0); 970 int nv=norm.length; 971 HealpixUtils.check(nv==rad.length,"inconsistent input arrays"); 972 RangeSet res = new RangeSet(); 973 974 if (scheme==Scheme.RING) 975 { 976 int fct=1; 977 if (inclusive) 978 { 979 HealpixUtils.check (((1L<<order_max)/nside)>=fact, 980 "invalid oversampling factor"); 981 fct = fact; 982 } 983 HealpixBase b2=null; 984 double rpsmall,rpbig; 985 if (fct>1) 986 { 987 b2=new HealpixBase(fct*nside,Scheme.RING); 988 rpsmall = b2.maxPixrad(); 989 rpbig = maxPixrad(); 990 } 991 else 992 rpsmall = rpbig = inclusive ? maxPixrad() : 0; 993 994 long irmin=1, irmax=nl4-1; 995 int nd=0; 996 double[] z0=new double[nv]; 997 double[] xa=new double[nv]; 998 double[] cosrsmall=new double[nv]; 999 double[] cosrbig=new double[nv]; 1000 Pointing[] ptg=new Pointing[nv]; 1001 long[] cpix=new long[nv]; 1002 Zphi[] czphi=new Zphi[nv]; 1003 for (int i=0; i<nv; ++i) 1004 { 1005 double rsmall=rad[i]+rpsmall; 1006 if (rsmall<Math.PI) 1007 { 1008 double rbig = Math.min(Math.PI,rad[i]+rpbig); 1009 Pointing pnt= new Pointing(norm[i]); 1010 cosrsmall[nd]=FastMath.cos(rsmall); 1011 cosrbig[nd]=FastMath.cos(rbig); 1012 double cth=FastMath.cos(pnt.theta); 1013 z0[nd]=cth; 1014 if (fct>1) cpix[nd]=zphi2pix(new Zphi(cth,pnt.phi)); 1015 if (fct>1) czphi[nd]=new Zphi(cth,pnt.phi); 1016 xa[nd]=1./Math.sqrt((1-cth)*(1+cth)); 1017 ptg[nd]=pnt; 1018 1019 double rlat1 = pnt.theta - rsmall; 1020 double zmax = FastMath.cos(rlat1); 1021 long irmin_t = (rlat1<=0) ? 1 : ringAbove (zmax)+1; 1022 1023 if ((fct>1) && (rlat1>0)) irmin_t=Math.max(1L,irmin_t-1); 1024 1025 double rlat2 = pnt.theta + rsmall; 1026 double zmin = FastMath.cos(rlat2); 1027 long irmax_t = (rlat2>=Math.PI) ? nl4-1 : ringAbove (zmin); 1028 1029 if ((fct>1) && (rlat2<Math.PI)) 1030 irmax_t=Math.min(nl4-1,irmax_t+1); 1031 if (irmax_t < irmax) irmax=irmax_t; 1032 if (irmin_t > irmin) irmin=irmin_t; 1033 1034 ++nd; 1035 } 1036 } 1037 1038 for (long iz=irmin; iz<=irmax; ++iz) 1039 { 1040 double z=ring2z(iz); 1041 RingInfoSmall ris=get_ring_info_small(iz); 1042 long ipix1=ris.startpix, nr=ris.ringpix; 1043 long ipix2 = ipix1 + nr - 1; // highest pixel number in the ring 1044 double shift = ris.shifted ? 0.5 : 0.; 1045 RangeSet rstmp = new RangeSet(); 1046 rstmp.append(ipix1,ipix2+1); 1047 1048 for (int j=0; (j<nd)&&(rstmp.nranges()>0); ++j) 1049 { 1050 double x = (cosrbig[j]-z*z0[j])*xa[j]; 1051 double ysq = 1.-z*z-x*x; 1052 double dphi = (ysq<=0) ? Math.PI-1e-15 : 1053 FastMath.atan2(Math.sqrt(ysq),x); 1054 1055 long ip_lo = (long)Math.floor 1056 (nr*Constants.inv_twopi*(ptg[j].phi-dphi)-shift)+1; 1057 long ip_hi = (long)Math.floor 1058 (nr*Constants.inv_twopi*(ptg[j].phi+dphi)-shift); 1059 1060 if (fct>1) 1061 { 1062 while ((ip_lo<=ip_hi) && checkPixelRing 1063 (b2,ip_lo,nr,ipix1,fct,czphi[j],cosrsmall[j],cpix[j])) 1064 ++ip_lo; 1065 while ((ip_hi>ip_lo) && checkPixelRing 1066 (b2,ip_hi,nr,ipix1,fct,czphi[j],cosrsmall[j],cpix[j])) 1067 --ip_hi; 1068 } 1069 1070 if (ip_hi>=nr) 1071 { ip_lo-=nr; ip_hi-=nr;} 1072 1073 if (ip_lo<0) 1074 { 1075 if (ip_hi+1<=ip_lo+nr-1) 1076 rstmp.remove(ipix1+ip_hi+1,ipix1+ip_lo+nr); 1077 } 1078 else 1079 rstmp.intersect(ipix1+ip_lo,ipix1+ip_hi+1); 1080 } 1081 res.append(rstmp); 1082 } 1083 } 1084 else // scheme == NEST 1085 { 1086 int oplus = 0; 1087 if (inclusive) 1088 { 1089 HealpixUtils.check ((1L<<(order_max-order))>=fact, 1090 "invalid oversampling factor"); 1091 HealpixUtils.check ((fact&(fact-1))==0, 1092 "oversampling factor must be a power of 2"); 1093 oplus=HealpixUtils.ilog2(fact); 1094 } 1095 int omax=order+oplus; // the order up to which we test 1096 1097 // TODO: ignore all disks with radius>=pi 1098 1099 double[][][] crlimit=new double[omax+1][nv][3]; 1100 for (int o=0; o<=omax; ++o) // prepare data at the required orders 1101 { 1102 double dr=HealpixProc.bn[o].maxPixrad(); // safety distance 1103 for (int i=0; i<nv; ++i) 1104 { 1105 crlimit[o][i][0] = (rad[i]+dr>Math.PI) ? -1: FastMath.cos(rad[i]+dr); 1106 crlimit[o][i][1] = (o==0) ? FastMath.cos(rad[i]) : crlimit[0][i][1]; 1107 crlimit[o][i][2] = (rad[i]-dr<0.) ? 1. : FastMath.cos(rad[i]-dr); 1108 } 1109 } 1110 1111 pstack stk=new pstack(12+3*omax); 1112 for (int i=0; i<12; i++) // insert the 12 base pixels in reverse order 1113 stk.push(11-i,0); 1114 1115 while (stk.size()>0) { // as long as there are pixels on the stack 1116 // pop current pixel number and order from the stack 1117 long pix=stk.ptop(); 1118 int o=stk.otop(); 1119 stk.pop(); 1120 1121 Vec3 pv = HealpixProc.bn[o].pix2vec(pix); 1122 1123 int zone=3; 1124 for (int i=0; (i<nv)&&(zone>0); ++i) 1125 { 1126 double crad=pv.dot(norm[i]); 1127 for (int iz=0; iz<zone; ++iz) 1128 if (crad<crlimit[o][i][iz]) 1129 zone=iz; 1130 } 1131 1132 if (zone>0) check_pixel (o, omax, zone, res, pix, stk, inclusive); 1133 } 1134 } 1135 return res; 1136 } 1137 1138 private RangeSet queryPolygonInternal (Pointing[] vertex, int fact) 1139 throws Exception 1140 { 1141 boolean inclusive = (fact!=0); 1142 int nv=vertex.length; 1143 int ncirc = inclusive ? nv+1 : nv; 1144 HealpixUtils.check(nv>=3,"not enough vertices in polygon"); 1145 Vec3[] vv = new Vec3[nv]; 1146 for (int i=0; i<nv; ++i) 1147 vv[i] = new Vec3(vertex[i]); 1148 Vec3[] normal = new Vec3[ncirc]; 1149 int flip=0; 1150 for (int i=0; i<nv; ++i) 1151 { 1152 normal[i]=vv[i].cross(vv[(i+1)%nv]).norm(); 1153 double hnd=normal[i].dot(vv[(i+2)%nv]); 1154 HealpixUtils.check(Math.abs(hnd)>1e-10,"degenerate corner"); 1155 if (i==0) 1156 flip = (hnd<0.) ? -1 : 1; 1157 else 1158 HealpixUtils.check(flip*hnd>0,"polygon is not convex"); 1159 normal[i].scale(flip); 1160 } 1161 double[] rad = new double[ncirc]; 1162 Arrays.fill(rad,Constants.halfpi); 1163 if (inclusive) 1164 { 1165 CircleFinder cf = new CircleFinder(vv); 1166 normal[nv]=cf.getCenter(); 1167 rad[nv]=FastMath.acos(cf.getCosrad()); 1168 } 1169 return queryMultiDisc(normal,rad,fact); 1170 } 1171 /** Returns a range set of pixels whose centers lie within the convex 1172 polygon defined by the {@code vertex} array. <p> 1173 This method is more efficient in the RING scheme. 1174 @param vertex an array containing the vertices of the requested convex 1175 polygon. 1176 @return the requested set of pixel number ranges */ 1177 public RangeSet queryPolygon (Pointing[] vertex) throws Exception 1178 { return queryPolygonInternal (vertex, 0); } 1179 /** Returns a range set of pixels that overlap with the convex 1180 polygon defined by the {@code vertex} array. <p> 1181 This method is more efficient in the RING scheme.<p> 1182 This method may return some pixels which don't overlap with 1183 the polygon at all. The higher {@code fact} is chosen, the fewer false 1184 positives are returned, at the cost of increased run time. 1185 @param vertex an array containing the vertices of the requested convex 1186 polygon. 1187 @param fact The overlapping test will be done at the resolution 1188 {@code fact*nside}. For NESTED ordering, {@code fact} must be a power 1189 of 2, else it can be any positive integer. A typical choice would be 4. 1190 @return the requested set of pixel number ranges */ 1191 public RangeSet queryPolygonInclusive (Pointing[] vertex, int fact) 1192 throws Exception 1193 { return queryPolygonInternal (vertex, fact); } 1194 1195 private void check_pixel (int o, int omax, int zone, 1196 RangeSet pixset, long pix, pstack stk, boolean inclusive) 1197 { 1198 if (zone==0) return; 1199 1200 if (o<order) 1201 { 1202 if (zone>=3) // output all subpixels 1203 { 1204 int sdist=2*(order-o); // the "bit-shift distance" between map orders 1205 pixset.append(pix<<sdist,((pix+1)<<sdist)); 1206 } 1207 else // (zone>=1) 1208 for (int i=0; i<4; ++i) 1209 stk.push(4*pix+3-i,o+1); // add children 1210 } 1211 else if (o>order) // this implies that inclusive==true 1212 { 1213 if (zone>=2) // pixel center in shape 1214 { 1215 pixset.append(pix>>>(2*(o-order))); // output the parent pixel at order 1216 stk.popToMark(); // unwind the stack 1217 } 1218 else // (zone>=1): pixel center in safety range 1219 { 1220 if (o<omax) // check sublevels 1221 for (int i=0; i<4; ++i) // add children in reverse order 1222 stk.push(4*pix+3-i,o+1); // add children 1223 else // at resolution limit 1224 { 1225 pixset.append(pix>>>(2*(o-order)));// output the parent pixel at order 1226 stk.popToMark(); // unwind the stack 1227 } 1228 } 1229 } 1230 else // o==order 1231 { 1232 if (zone>=2) 1233 pixset.append(pix); 1234 else if (inclusive) // and (zone>=1) 1235 { 1236 if (order<omax) // check sublevels 1237 { 1238 stk.mark(); // remember current stack position 1239 for (int i=0; i<4; ++i) // add children in reverse order 1240 stk.push(4*pix+3-i,o+1); // add children 1241 } 1242 else // at resolution limit 1243 pixset.append(pix); // output the pixel 1244 } 1245 } 1246 } 1247 1248 /** Compute ring index from pixel number. 1249 Works in both RING and NESTED schemes 1250 @param pix pixel number 1251 @return ring index (1 to 4Nside-1) */ 1252 public long pix2ring (long pix) 1253 { 1254 if (scheme==Scheme.RING) 1255 { 1256 if (pix<ncap) // North Polar cap 1257 return (1+HealpixUtils.isqrt(1+2*pix))>>>1; // counted from North pole 1258 else if (pix<(npix-ncap)) // Equatorial region 1259 return (pix-ncap)/nl4 + nside; // counted from North pole 1260 else // South Polar cap 1261 return nl4-((1+HealpixUtils.isqrt(2*(npix-pix)-1))>>>1); 1262 } 1263 else 1264 { 1265 Xyf xyf = nest2xyf(pix); 1266 return ((long)(jrll[xyf.face])<<order) - xyf.ix - xyf.iy - 1; 1267 } 1268 } 1269 1270 /** Returns a set of points along the boundary of the given pixel. 1271 * Step 1 gives 4 points on the corners. The first point corresponds 1272 * to the northernmost corner, the subsequent points follow the pixel 1273 * boundary through west, south and east corners. 1274 * 1275 * @param pix pixel index number 1276 * @param step the number of returned points is 4*step 1277 * @return {@link Vec3} for each point 1278 * @throws Exception 1279 */ 1280 public Vec3[] boundaries(long pix, int step) throws Exception 1281 { 1282 HealpixUtils.check(step>0,"step must be positive"); 1283 Vec3[] points = new Vec3[4*step]; 1284 Xyf xyf = pix2xyf(pix); 1285 double dc=0.5/nside; 1286 double xc=(xyf.ix+0.5)/nside, yc=(xyf.iy+0.5)/nside; 1287 double d = 1./(step*nside); 1288 for (int i=0; i<step; ++i) 1289 { 1290 points[i ]=new Fxyf(xc+dc-i*d, yc+dc, xyf.face).toVec3(); 1291 points[i+ step]=new Fxyf(xc-dc, yc+dc-i*d, xyf.face).toVec3(); 1292 points[i+2*step]=new Fxyf(xc-dc+i*d, yc-dc, xyf.face).toVec3(); 1293 points[i+3*step]=new Fxyf(xc+dc, yc-dc+i*d, xyf.face).toVec3(); 1294 } 1295 return points; 1296 } 1297 } 1298