1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <math.h>
7 #include <assert.h>
8 #include <stdio.h>
9 #include <string.h>
10
11 #include "os-features.h"
12 #include "healpix.h"
13 #include "mathutil.h"
14 #include "starutil.h"
15 #include "keywords.h"
16 #include "permutedsort.h"
17 #include "log.h"
18
19 // Internal type
20 struct hp_s {
21 int bighp;
22 int x;
23 int y;
24 };
25 typedef struct hp_s hp_t;
26
hptointl(hp_t hp,int Nside)27 static int64_t hptointl(hp_t hp, int Nside) {
28 return healpix_compose_xyl(hp.bighp, hp.x, hp.y, Nside);
29 }
30
hptoint(hp_t hp,int Nside)31 static int hptoint(hp_t hp, int Nside) {
32 return healpix_compose_xy(hp.bighp, hp.x, hp.y, Nside);
33 }
34
intltohp(int64_t pix,hp_t * hp,int Nside)35 static void intltohp(int64_t pix, hp_t* hp, int Nside) {
36 healpix_decompose_xyl(pix, &hp->bighp, &hp->x, &hp->y, Nside);
37 }
38
inttohp(int pix,hp_t * hp,int Nside)39 static void inttohp(int pix, hp_t* hp, int Nside) {
40 healpix_decompose_xy(pix, &hp->bighp, &hp->x, &hp->y, Nside);
41 }
42
hp_decompose(hp_t * hp,int * php,int * px,int * py)43 static void hp_decompose(hp_t* hp, int* php, int* px, int* py) {
44 if (php)
45 *php = hp->bighp;
46 if (px)
47 *px = hp->x;
48 if (py)
49 *py = hp->y;
50 }
51
52 // I've had troubles with rounding functions being declared properly
53 // in other contexts... Declare it here so the compiler complains if
54 // something is wrong.
55 double round(double x);
56
mysquare(double d)57 Const static Inline double mysquare(double d) {
58 return d*d;
59 }
60
healpix_xy_to_nested(int hp,int Nside)61 Const int healpix_xy_to_nested(int hp, int Nside) {
62 int bighp,x,y;
63 int index;
64 int i;
65
66 healpix_decompose_xy(hp, &bighp, &x, &y, Nside);
67 if (!is_power_of_two(Nside)) {
68 fprintf(stderr, "healpix_xy_to_nested: Nside must be a power of two.\n");
69 return -1;
70 }
71
72 // We construct the index called p_n' in the healpix paper, whose bits
73 // are taken from the bits of x and y:
74 // x = ... b4 b2 b0
75 // y = ... b5 b3 b1
76 // We go through the bits of x,y, building up "index":
77 index = 0;
78 for (i=0; i<(8*sizeof(int)/2); i++) {
79 index |= (((y & 1) << 1) | (x & 1)) << (i*2);
80 y >>= 1;
81 x >>= 1;
82 if (!x && !y) break;
83 }
84
85 return index + bighp * Nside * Nside;
86 }
87
healpix_nested_to_xy(int hp,int Nside)88 Const int healpix_nested_to_xy(int hp, int Nside) {
89 int bighp, x, y;
90 int index;
91 int i;
92 if (!is_power_of_two(Nside)) {
93 fprintf(stderr, "healpix_xy_to_nested: Nside must be a power of two.\n");
94 return -1;
95 }
96 bighp = hp / (Nside*Nside);
97 index = hp % (Nside*Nside);
98 x = y = 0;
99 for (i=0; i<(8*sizeof(int)/2); i++) {
100 x |= (index & 0x1) << i;
101 index >>= 1;
102 y |= (index & 0x1) << i;
103 index >>= 1;
104 if (!index) break;
105 }
106 return healpix_compose_xy(bighp, x, y, Nside);
107 }
108
healpix_compose_ring(int ring,int longind,int Nside)109 Const int healpix_compose_ring(int ring, int longind, int Nside) {
110 if (ring <= Nside)
111 // north polar
112 return ring * (ring-1) * 2 + longind;
113 if (ring < 3*Nside)
114 // equatorial
115 return Nside*(Nside-1)*2 + Nside*4*(ring-Nside) + longind;
116 {
117 int ri;
118 ri = 4*Nside - ring;
119 return 12*Nside*Nside-1 - ( ri*(ri-1)*2 + (ri*4 - 1 - longind) );
120 }
121 }
122
healpix_decompose_ring(int hp,int Nside,int * p_ring,int * p_longind)123 void healpix_decompose_ring(int hp, int Nside, int* p_ring, int* p_longind) {
124 // FIXME: this could be written in closed form...
125 int longind;
126 int ring;
127 int offset = 0;
128 for (ring=1; ring<=Nside; ring++) {
129 if (offset + ring*4 > hp) {
130 longind = hp - offset;
131 goto gotit;
132 }
133 offset += ring*4;
134 }
135 for (; ring<(3*Nside); ring++) {
136 if (offset + Nside*4 > hp) {
137 longind = hp - offset;
138 goto gotit;
139 }
140 offset += Nside*4;
141 }
142 for (; ring<(4*Nside); ring++) {
143 if (offset + (Nside*4 - ring)*4 > hp) {
144 longind = hp - offset;
145 goto gotit;
146 }
147 offset += (Nside*4 - ring)*4;
148 }
149 fprintf(stderr, "healpix_decompose_ring: shouldn't get here!\n");
150 if (p_ring) *p_ring = -1;
151 if (p_longind) *p_longind = -1;
152 return;
153 gotit:
154 if (p_ring)
155 *p_ring = ring;
156 if (p_longind)
157 *p_longind = longind;
158 }
159
healpix_ring_to_xy(int ring,int Nside)160 Const int healpix_ring_to_xy(int ring, int Nside) {
161 int bighp, x, y;
162 int ringind, longind;
163 healpix_decompose_ring(ring, Nside, &ringind, &longind);
164 if (ringind <= Nside) {
165 int ind;
166 int v;
167 int F1;
168 int frow;
169 bighp = longind / ringind;
170 ind = longind - bighp * ringind;
171 y = (Nside - 1 - ind);
172 frow = bighp / 4;
173 F1 = frow + 2;
174 v = F1*Nside - ringind - 1;
175 x = v - y;
176 return healpix_compose_xy(bighp, x, y, Nside);
177 } else if (ringind < 3*Nside) {
178 int panel;
179 int ind;
180 int bottomleft;
181 int topleft;
182 int frow, F1, F2, s, v, h;
183 int bighp = -1;
184 int x, y;
185 int hp;
186 int R = 0;
187
188 panel = longind / Nside;
189 ind = longind % Nside;
190 bottomleft = ind < (ringind - Nside + 1) / 2;
191 topleft = ind < (3*Nside - ringind + 1)/2;
192
193 if (!bottomleft && topleft) {
194 // top row.
195 bighp = panel;
196 } else if (bottomleft && !topleft) {
197 // bottom row.
198 bighp = 8 + panel;
199 } else if (bottomleft && topleft) {
200 // left side.
201 bighp = 4 + panel;
202 } else if (!bottomleft && !topleft) {
203 // right side.
204 bighp = 4 + (panel + 1) % 4;
205 if (bighp == 4) {
206 longind -= (4*Nside - 1);
207 // Gah! Wacky hack - it seems that since
208 // "longind" is negative in this case, the
209 // rounding behaves differently, so we end up
210 // computing the wrong "h" and have to correct
211 // for it.
212 R = 1;
213 }
214 }
215
216 frow = bighp / 4;
217 F1 = frow + 2;
218 F2 = 2*(bighp % 4) - (frow % 2) + 1;
219 s = (ringind - Nside) % 2;
220 v = F1*Nside - ringind - 1;
221 h = 2*longind - s - F2*Nside;
222 if (R)
223 h--;
224 x = (v + h) / 2;
225 y = (v - h) / 2;
226 //fprintf(stderr, "bighp=%i, frow=%i, F1=%i, F2=%i, s=%i, v=%i, h=%i, x=%i, y=%i.\n", bighp, frow, F1, F2, s, v, h, x, y);
227
228 if ((v != (x+y)) || (h != (x-y))) {
229 h++;
230 x = (v + h) / 2;
231 y = (v - h) / 2;
232 //fprintf(stderr, "tweak h=%i, x=%i, y=%i\n", h, x, y);
233
234 if ((v != (x+y)) || (h != (x-y))) {
235 //fprintf(stderr, "still not right.\n");
236 }
237 }
238 hp = healpix_compose_xy(bighp, x, y, Nside);
239 //fprintf(stderr, "hp %i\n", hp);
240 return hp;
241 } else {
242 int ind;
243 int v;
244 int F1;
245 int frow;
246 int ri;
247 ri = 4*Nside - ringind;
248 bighp = 8 + longind / ri;
249 ind = longind - (bighp%4) * ri;
250 y = (ri-1) - ind;
251 frow = bighp / 4;
252 F1 = frow + 2;
253 v = F1*Nside - ringind - 1;
254 x = v - y;
255 return healpix_compose_xy(bighp, x, y, Nside);
256 }
257 }
258
healpix_xy_to_ring(int hp,int Nside)259 Const int healpix_xy_to_ring(int hp, int Nside) {
260 int bighp,x,y;
261 int frow;
262 int F1;
263 int v;
264 int ring;
265 int index;
266
267 healpix_decompose_xy(hp, &bighp, &x, &y, Nside);
268 frow = bighp / 4;
269 F1 = frow + 2;
270 v = x + y;
271 // "ring" starts from 1 at the north pole and goes to 4Nside-1 at
272 // the south pole; the pixels in each ring have the same latitude.
273 ring = F1*Nside - v - 1;
274 /*
275 ring:
276 [1, Nside] : n pole
277 (Nside, 2Nside] : n equatorial
278 (2Nside+1, 3Nside) : s equat
279 [3Nside, 4Nside-1] : s pole
280 */
281 // this probably can't happen...
282 if ((ring < 1) || (ring >= 4*Nside)) {
283 fprintf(stderr, "Invalid ring index: %i\n", ring);
284 return -1;
285 }
286 if (ring <= Nside) {
287 // north polar.
288 // left-to-right coordinate within this healpix
289 index = (Nside - 1 - y);
290 // offset from the other big healpixes
291 index += ((bighp % 4) * ring);
292 // offset from the other rings
293 index += ring*(ring-1)*2;
294 } else if (ring >= 3*Nside) {
295 // south polar.
296 // Here I first flip everything so that we label the pixels
297 // at zero starting in the southeast corner, increasing to the
298 // west and north, then subtract that from the total number of
299 // healpixels.
300 int ri = 4*Nside - ring;
301 // index within this healpix
302 index = (ri-1) - x;
303 // big healpixes
304 index += ((3-(bighp % 4)) * ri);
305 // other rings
306 index += ri*(ri-1)*2;
307 // flip!
308 index = 12*Nside*Nside - 1 - index;
309 } else {
310 // equatorial.
311 int s, F2, h;
312 s = (ring - Nside) % 2;
313 F2 = 2*((int)bighp % 4) - (frow % 2) + 1;
314 h = x - y;
315 index = (F2 * (int)Nside + h + s) / 2;
316 // offset from the north polar region:
317 index += Nside*(Nside-1)*2;
318 // offset within the equatorial region:
319 index += Nside * 4 * (ring - Nside);
320 // handle healpix #4 wrap-around
321 if ((bighp == 4) && (y > x))
322 index += (4 * Nside - 1);
323 //fprintf(stderr, "frow=%i, F1=%i, v=%i, ringind=%i, s=%i, F2=%i, h=%i, longind=%i.\n", frow, F1, v, ring, s, F2, h, (F2*(int)Nside+h+s)/2);
324 }
325 return index;
326 }
327
healpix_side_length_arcmin(int Nside)328 Const double healpix_side_length_arcmin(int Nside) {
329 return sqrt((4.0 * M_PI * mysquare(180.0 * 60.0 / M_PI)) /
330 (12.0 * Nside * Nside));
331 }
332
healpix_nside_for_side_length_arcmin(double arcmin)333 double healpix_nside_for_side_length_arcmin(double arcmin) {
334 return sqrt(4.0*M_PI / (mysquare(arcmin2rad(arcmin)) * 12.0));
335 }
336
swap(int * i1,int * i2)337 static Inline void swap(int* i1, int* i2) {
338 int tmp;
339 tmp = *i1;
340 *i1 = *i2;
341 *i2 = tmp;
342 }
343
swap_double(double * i1,double * i2)344 static Inline void swap_double(double* i1, double* i2) {
345 double tmp;
346 tmp = *i1;
347 *i1 = *i2;
348 *i2 = tmp;
349 }
350
ispolar(int healpix)351 static Inline anbool ispolar(int healpix)
352 {
353 // the north polar healpixes are 0,1,2,3
354 // the south polar healpixes are 8,9,10,11
355 return (healpix <= 3) || (healpix >= 8);
356 }
357
isequatorial(int healpix)358 static Inline anbool isequatorial(int healpix)
359 {
360 // the north polar healpixes are 0,1,2,3
361 // the south polar healpixes are 8,9,10,11
362 return (healpix >= 4) && (healpix <= 7);
363 }
364
isnorthpolar(int healpix)365 static Inline anbool isnorthpolar(int healpix)
366 {
367 return (healpix <= 3);
368 }
369
issouthpolar(int healpix)370 static Inline anbool issouthpolar(int healpix)
371 {
372 return (healpix >= 8);
373 }
374
compose_xy(int x,int y,int Nside)375 static int compose_xy(int x, int y, int Nside) {
376 assert(Nside > 0);
377 assert(x >= 0);
378 assert(x < Nside);
379 assert(y >= 0);
380 assert(y < Nside);
381 return (x * Nside) + y;
382 }
383
healpix_compose_xy(int bighp,int x,int y,int Nside)384 int healpix_compose_xy(int bighp, int x, int y, int Nside) {
385 assert(bighp >= 0);
386 assert(bighp < 12);
387 return (bighp * Nside * Nside) + compose_xy(x, y, Nside);
388 }
389
healpix_compose_xyl(int bighp,int x,int y,int Nside)390 int64_t healpix_compose_xyl(int bighp, int x, int y, int Nside) {
391 int64_t ns = Nside;
392 assert(Nside > 0);
393 assert(bighp >= 0);
394 assert(bighp < 12);
395 assert(x >= 0);
396 assert(x < Nside);
397 assert(y >= 0);
398 assert(y < Nside);
399 return ((((int64_t)bighp * ns) + x) * ns) + y;
400 }
401
healpix_convert_nside(int hp,int nside,int outnside,int * outhp)402 void healpix_convert_nside(int hp, int nside, int outnside, int* outhp) {
403 int basehp, x, y;
404 int ox, oy;
405 healpix_decompose_xy(hp, &basehp, &x, &y, nside);
406 healpix_convert_xy_nside(x, y, nside, outnside, &ox, &oy);
407 *outhp = healpix_compose_xy(basehp, ox, oy, outnside);
408 }
409
healpix_convert_xy_nside(int x,int y,int nside,int outnside,int * outx,int * outy)410 void healpix_convert_xy_nside(int x, int y, int nside, int outnside,
411 int* outx, int* outy) {
412 double fx, fy;
413 int ox, oy;
414 assert(x >= 0);
415 assert(x < nside);
416 assert(y >= 0);
417 assert(y < nside);
418
419 // MAGIC 0.5: assume center of pixel...
420 fx = (x + 0.5) / (double)nside;
421 fy = (y + 0.5) / (double)nside;
422
423 ox = floor(fx * outnside);
424 oy = floor(fy * outnside);
425
426 if (outx)
427 *outx = ox;
428 if (outy)
429 *outy = oy;
430 }
431
healpix_decompose_xy(int finehp,int * pbighp,int * px,int * py,int Nside)432 void healpix_decompose_xy(int finehp, int* pbighp, int* px, int* py, int Nside) {
433 int hp;
434 assert(Nside > 0);
435 assert(finehp < (12 * Nside * Nside));
436 assert(finehp >= 0);
437 if (pbighp) {
438 int bighp = finehp / (Nside * Nside);
439 assert(bighp >= 0);
440 assert(bighp < 12);
441 *pbighp = bighp;
442 }
443 hp = finehp % (Nside * Nside);
444 if (px) {
445 *px = hp / Nside;
446 assert(*px >= 0);
447 assert(*px < Nside);
448 }
449 if (py) {
450 *py = hp % Nside;
451 assert(*py >= 0);
452 assert(*py < Nside);
453 }
454 }
455
healpix_decompose_xyl(int64_t finehp,int * pbighp,int * px,int * py,int Nside)456 void healpix_decompose_xyl(int64_t finehp,
457 int* pbighp, int* px, int* py,
458 int Nside) {
459 int64_t hp;
460 int64_t ns2 = (int64_t)Nside * (int64_t)Nside;
461 assert(Nside > 0);
462 assert(finehp < (12L * ns2));
463 assert(finehp >= 0);
464 if (pbighp) {
465 int bighp = finehp / ns2;
466 assert(bighp >= 0);
467 assert(bighp < 12);
468 *pbighp = bighp;
469 }
470 hp = finehp % ns2;
471 if (px) {
472 *px = hp / Nside;
473 assert(*px >= 0);
474 assert(*px < Nside);
475 }
476 if (py) {
477 *py = hp % Nside;
478 assert(*py >= 0);
479 assert(*py < Nside);
480 }
481 }
482
483 /**
484 Given a large-scale healpix number, computes its neighbour in the
485 direction (dx,dy). Returns -1 if there is no such neighbour.
486 */
healpix_get_neighbour(int hp,int dx,int dy)487 static int healpix_get_neighbour(int hp, int dx, int dy)
488 {
489 if (isnorthpolar(hp)) {
490 if ((dx == 1) && (dy == 0)) return (hp + 1) % 4;
491 if ((dx == 0) && (dy == 1)) return (hp + 3) % 4;
492 if ((dx == 1) && (dy == 1)) return (hp + 2) % 4;
493 if ((dx == -1) && (dy == 0)) return (hp + 4);
494 if ((dx == 0) && (dy == -1)) return 4 + ((hp + 1) % 4);
495 if ((dx == -1) && (dy == -1)) return hp + 8;
496 return -1;
497 } else if (issouthpolar(hp)) {
498 if ((dx == 1) && (dy == 0)) return 4 + ((hp + 1) % 4);
499 if ((dx == 0) && (dy == 1)) return hp - 4;
500 if ((dx == -1) && (dy == 0)) return 8 + ((hp + 3) % 4);
501 if ((dx == 0) && (dy == -1)) return 8 + ((hp + 1) % 4);
502 if ((dx == -1) && (dy == -1)) return 8 + ((hp + 2) % 4);
503 if ((dx == 1) && (dy == 1)) return hp - 8;
504 return -1;
505 } else {
506 if ((dx == 1) && (dy == 0)) return hp - 4;
507 if ((dx == 0) && (dy == 1)) return (hp + 3) % 4;
508 if ((dx == -1) && (dy == 0)) return 8 + ((hp + 3) % 4);
509 if ((dx == 0) && (dy == -1)) return hp + 4;
510 if ((dx == 1) && (dy == -1)) return 4 + ((hp + 1) % 4);
511 if ((dx == -1) && (dy == 1)) return 4 + ((hp - 1) % 4);
512 return -1;
513 }
514 return -1;
515 }
516
get_neighbours(hp_t hp,hp_t * neighbour,int Nside)517 static int get_neighbours(hp_t hp, hp_t* neighbour, int Nside) {
518 int base;
519 int x, y;
520 int nn = 0;
521 int nbase;
522 int nx, ny;
523
524 base = hp.bighp;
525 x = hp.x;
526 y = hp.y;
527
528 // ( + , 0 )
529 nx = (x + 1) % Nside;
530 ny = y;
531 if (x == (Nside - 1)) {
532 nbase = healpix_get_neighbour(base, 1, 0);
533 if (isnorthpolar(base)) {
534 nx = x;
535 swap(&nx, &ny);
536 }
537 } else
538 nbase = base;
539
540 neighbour[nn].bighp = nbase;
541 neighbour[nn].x = nx;
542 neighbour[nn].y = ny;
543 nn++;
544
545 // ( + , + )
546 nx = (x + 1) % Nside;
547 ny = (y + 1) % Nside;
548 if ((x == Nside - 1) && (y == Nside - 1)) {
549 if (ispolar(base))
550 nbase = healpix_get_neighbour(base, 1, 1);
551 else
552 nbase = -1;
553 } else if (x == (Nside - 1))
554 nbase = healpix_get_neighbour(base, 1, 0);
555 else if (y == (Nside - 1))
556 nbase = healpix_get_neighbour(base, 0, 1);
557 else
558 nbase = base;
559
560 if (isnorthpolar(base)) {
561 if (x == (Nside - 1))
562 nx = Nside - 1;
563 if (y == (Nside - 1))
564 ny = Nside - 1;
565 if ((x == (Nside - 1)) || (y == (Nside - 1)))
566 swap(&nx, &ny);
567 }
568
569 //printf("(+ +): nbase=%i, nx=%i, ny=%i, pix=%i\n", nbase, nx, ny, nbase*Ns2+xy_to_pnprime(nx,ny,Nside));
570
571 if (nbase != -1) {
572 neighbour[nn].bighp = nbase;
573 neighbour[nn].x = nx;
574 neighbour[nn].y = ny;
575 nn++;
576 }
577
578 // ( 0 , + )
579 nx = x;
580 ny = (y + 1) % Nside;
581 if (y == (Nside - 1)) {
582 nbase = healpix_get_neighbour(base, 0, 1);
583 if (isnorthpolar(base)) {
584 ny = y;
585 swap(&nx, &ny);
586 }
587 } else
588 nbase = base;
589
590 //printf("(0 +): nbase=%i, nx=%i, ny=%i, pix=%i\n", nbase, nx, ny, nbase*Ns2+xy_to_pnprime(nx,ny,Nside));
591
592 neighbour[nn].bighp = nbase;
593 neighbour[nn].x = nx;
594 neighbour[nn].y = ny;
595 nn++;
596
597 // ( - , + )
598 nx = (x + Nside - 1) % Nside;
599 ny = (y + 1) % Nside;
600 if ((x == 0) && (y == (Nside - 1))) {
601 if (isequatorial(base))
602 nbase = healpix_get_neighbour(base, -1, 1);
603 else
604 nbase = -1;
605 } else if (x == 0) {
606 nbase = healpix_get_neighbour(base, -1, 0);
607 if (issouthpolar(base)) {
608 nx = 0;
609 swap(&nx, &ny);
610 }
611 } else if (y == (Nside - 1)) {
612 nbase = healpix_get_neighbour(base, 0, 1);
613 if (isnorthpolar(base)) {
614 ny = y;
615 swap(&nx, &ny);
616 }
617 } else
618 nbase = base;
619
620 //printf("(- +): nbase=%i, nx=%i, ny=%i, pix=%i\n", nbase, nx, ny, nbase*Ns2+xy_to_pnprime(nx,ny,Nside));
621
622 if (nbase != -1) {
623 neighbour[nn].bighp = nbase;
624 neighbour[nn].x = nx;
625 neighbour[nn].y = ny;
626 nn++;
627 }
628
629 // ( - , 0 )
630 nx = (x + Nside - 1) % Nside;
631 ny = y;
632 if (x == 0) {
633 nbase = healpix_get_neighbour(base, -1, 0);
634 if (issouthpolar(base)) {
635 nx = 0;
636 swap(&nx, &ny);
637 }
638 } else
639 nbase = base;
640
641 //printf("(- 0): nbase=%i, nx=%i, ny=%i, pix=%i\n", nbase, nx, ny, nbase*Ns2+xy_to_pnprime(nx,ny,Nside));
642
643 neighbour[nn].bighp = nbase;
644 neighbour[nn].x = nx;
645 neighbour[nn].y = ny;
646 nn++;
647
648 // ( - , - )
649 nx = (x + Nside - 1) % Nside;
650 ny = (y + Nside - 1) % Nside;
651 if ((x == 0) && (y == 0)) {
652 if (ispolar(base))
653 nbase = healpix_get_neighbour(base, -1, -1);
654 else
655 nbase = -1;
656 } else if (x == 0)
657 nbase = healpix_get_neighbour(base, -1, 0);
658 else if (y == 0)
659 nbase = healpix_get_neighbour(base, 0, -1);
660 else
661 nbase = base;
662
663 if (issouthpolar(base)) {
664 if (x == 0)
665 nx = 0;
666 if (y == 0)
667 ny = 0;
668 if ((x == 0) || (y == 0))
669 swap(&nx, &ny);
670 }
671
672 //printf("(- -): nbase=%i, nx=%i, ny=%i, pix=%i\n", nbase, nx, ny, nbase*Ns2+xy_to_pnprime(nx,ny,Nside));
673
674 if (nbase != -1) {
675 neighbour[nn].bighp = nbase;
676 neighbour[nn].x = nx;
677 neighbour[nn].y = ny;
678 nn++;
679 }
680
681 // ( 0 , - )
682 ny = (y + Nside - 1) % Nside;
683 nx = x;
684 if (y == 0) {
685 nbase = healpix_get_neighbour(base, 0, -1);
686 if (issouthpolar(base)) {
687 ny = y;
688 swap(&nx, &ny);
689 }
690 } else
691 nbase = base;
692
693 //printf("(0 -): nbase=%i, nx=%i, ny=%i, pix=%i\n", nbase, nx, ny, nbase*Ns2+xy_to_pnprime(nx,ny,Nside));
694
695 neighbour[nn].bighp = nbase;
696 neighbour[nn].x = nx;
697 neighbour[nn].y = ny;
698 nn++;
699
700 // ( + , - )
701 nx = (x + 1) % Nside;
702 ny = (y + Nside - 1) % Nside;
703 if ((x == (Nside - 1)) && (y == 0)) {
704 if (isequatorial(base)) {
705 nbase = healpix_get_neighbour(base, 1, -1);
706 } else
707 nbase = -1;
708
709 } else if (x == (Nside - 1)) {
710 nbase = healpix_get_neighbour(base, 1, 0);
711 if (isnorthpolar(base)) {
712 nx = x;
713 swap(&nx, &ny);
714 }
715 } else if (y == 0) {
716 nbase = healpix_get_neighbour(base, 0, -1);
717 if (issouthpolar(base)) {
718 ny = y;
719 swap(&nx, &ny);
720 }
721 } else
722 nbase = base;
723
724 //printf("(+ -): nbase=%i, nx=%i, ny=%i, pix=%i\n", nbase, nx, ny, nbase*Ns2+xy_to_pnprime(nx,ny,Nside));
725
726 if (nbase != -1) {
727 neighbour[nn].bighp = nbase;
728 neighbour[nn].x = nx;
729 neighbour[nn].y = ny;
730 nn++;
731 }
732
733 return nn;
734 }
735
healpix_get_neighbours(int pix,int * neighbour,int Nside)736 int healpix_get_neighbours(int pix, int* neighbour, int Nside) {
737 hp_t neigh[8];
738 hp_t hp;
739 int nn;
740 int i;
741 inttohp(pix, &hp, Nside);
742 nn = get_neighbours(hp, neigh, Nside);
743 for (i=0; i<nn; i++)
744 neighbour[i] = hptoint(neigh[i], Nside);
745 return nn;
746 }
747
healpix_get_neighboursl(int64_t pix,int64_t * neighbour,int Nside)748 int healpix_get_neighboursl(int64_t pix, int64_t* neighbour, int Nside) {
749 hp_t neigh[8];
750 hp_t hp;
751 int nn;
752 int i;
753 intltohp(pix, &hp, Nside);
754 nn = get_neighbours(hp, neigh, Nside);
755 for (i=0; i<nn; i++)
756 neighbour[i] = hptointl(neigh[i], Nside);
757 return nn;
758 }
759
xyztohp(double vx,double vy,double vz,int Nside,double * p_dx,double * p_dy)760 static hp_t xyztohp(double vx, double vy, double vz, int Nside,
761 double* p_dx, double* p_dy) {
762 double phi;
763 double twothirds = 2.0 / 3.0;
764 double pi = M_PI;
765 double twopi = 2.0 * M_PI;
766 double halfpi = 0.5 * M_PI;
767 double dx, dy;
768 int basehp;
769 int x, y;
770 double sector;
771 int offset;
772 double phi_t;
773 hp_t hp;
774
775 // only used in asserts()
776 VarUnused double EPS = 1e-8;
777
778 assert(Nside > 0);
779
780 /* Convert our point into cylindrical coordinates for middle ring */
781 phi = atan2(vy, vx);
782 if (phi < 0.0)
783 phi += twopi;
784 phi_t = fmod(phi, halfpi);
785 assert (phi_t >= 0.0);
786
787 // North or south polar cap.
788 if ((vz >= twothirds) || (vz <= -twothirds)) {
789 double zfactor;
790 anbool north;
791 int column;
792 double root;
793 double xx, yy, kx, ky;
794
795 // Which pole?
796 if (vz >= twothirds) {
797 north = TRUE;
798 zfactor = 1.0;
799 } else {
800 north = FALSE;
801 zfactor = -1.0;
802 }
803
804 // solve eqn 20: k = Ns - xx (in the northern hemi)
805 root = (1.0 - vz*zfactor) * 3.0 * mysquare(Nside * (2.0 * phi_t - pi) / pi);
806 kx = (root <= 0.0) ? 0.0 : sqrt(root);
807
808 // solve eqn 19 for k = Ns - yy
809 root = (1.0 - vz*zfactor) * 3.0 * mysquare(Nside * 2.0 * phi_t / pi);
810 ky = (root <= 0.0) ? 0.0 : sqrt(root);
811
812 if (north) {
813 xx = Nside - kx;
814 yy = Nside - ky;
815 } else {
816 xx = ky;
817 yy = kx;
818 }
819
820 // xx, yy should be in [0, Nside].
821 x = MIN(Nside-1, floor(xx));
822 assert(x >= 0);
823 assert(x < Nside);
824
825 y = MIN(Nside-1, floor(yy));
826 assert(y >= 0);
827 assert(y < Nside);
828
829 dx = xx - x;
830 dy = yy - y;
831
832 sector = (phi - phi_t) / (halfpi);
833 offset = (int)round(sector);
834 assert(fabs(sector - offset) < EPS);
835 offset = ((offset % 4) + 4) % 4;
836 assert(offset >= 0);
837 assert(offset <= 3);
838 column = offset;
839
840 if (north)
841 basehp = column;
842 else
843 basehp = 8 + column;
844
845 } else {
846 // could be polar or equatorial.
847 double sector;
848 int offset;
849 double u1, u2;
850 double zunits, phiunits;
851 double xx, yy;
852
853 // project into the unit square z=[-2/3, 2/3], phi=[0, pi/2]
854 zunits = (vz + twothirds) / (4.0 / 3.0);
855 phiunits = phi_t / halfpi;
856 // convert into diagonal units
857 // (add 1 to u2 so that they both cover the range [0,2].
858 u1 = zunits + phiunits;
859 u2 = zunits - phiunits + 1.0;
860 assert(u1 >= 0.);
861 assert(u1 <= 2.);
862 assert(u2 >= 0.);
863 assert(u2 <= 2.);
864 // x is the northeast direction, y is the northwest.
865 xx = u1 * Nside;
866 yy = u2 * Nside;
867
868 // now compute which big healpix it's in.
869 // (note that we subtract off the modded portion used to
870 // compute the position within the healpix, so this should be
871 // very close to one of the boundaries.)
872 sector = (phi - phi_t) / (halfpi);
873 offset = (int)round(sector);
874 assert(fabs(sector - offset) < EPS);
875 offset = ((offset % 4) + 4) % 4;
876 assert(offset >= 0);
877 assert(offset <= 3);
878
879 // we're looking at a square in z,phi space with an X dividing it.
880 // we want to know which section we're in.
881 // xx ranges from 0 in the bottom-left to 2Nside in the top-right.
882 // yy ranges from 0 in the bottom-right to 2Nside in the top-left.
883 // (of the phi,z unit box)
884 if (xx >= Nside) {
885 xx -= Nside;
886 if (yy >= Nside) {
887 // north polar.
888 yy -= Nside;
889 basehp = offset;
890 } else {
891 // right equatorial.
892 basehp = ((offset + 1) % 4) + 4;
893 }
894 } else {
895 if (yy >= Nside) {
896 // left equatorial.
897 yy -= Nside;
898 basehp = offset + 4;
899 } else {
900 // south polar.
901 basehp = 8 + offset;
902 }
903 }
904
905 assert(xx >= -EPS);
906 assert(xx < (Nside+EPS));
907 x = MAX(0, MIN(Nside-1, floor(xx)));
908 assert(x >= 0);
909 assert(x < Nside);
910 dx = xx - x;
911
912 assert(yy >= -EPS);
913 assert(yy < (Nside+EPS));
914 y = MAX(0, MIN(Nside-1, floor(yy)));
915 assert(y >= 0);
916 assert(y < Nside);
917 dy = yy - y;
918 }
919
920 hp.bighp = basehp;
921 hp.x = x;
922 hp.y = y;
923
924 if (p_dx) *p_dx = dx;
925 if (p_dy) *p_dy = dy;
926
927 return hp;
928 }
929
xyztohealpix(double x,double y,double z,int Nside)930 int xyztohealpix(double x, double y, double z, int Nside) {
931 return xyztohealpixf(x, y, z, Nside, NULL, NULL);
932 }
933
xyztohealpixl(double x,double y,double z,int Nside)934 int64_t xyztohealpixl(double x, double y, double z, int Nside) {
935 return xyztohealpixlf(x, y, z, Nside, NULL, NULL);
936 }
937
xyztohealpixlf(double x,double y,double z,int Nside,double * p_dx,double * p_dy)938 int64_t xyztohealpixlf(double x, double y, double z, int Nside,
939 double* p_dx, double* p_dy) {
940 hp_t hp = xyztohp(x,y,z, Nside, p_dx,p_dy);
941 return hptointl(hp, Nside);
942 }
943
xyztohealpixf(double x,double y,double z,int Nside,double * p_dx,double * p_dy)944 int xyztohealpixf(double x, double y, double z, int Nside,
945 double* p_dx, double* p_dy) {
946 hp_t hp = xyztohp(x,y,z, Nside, p_dx,p_dy);
947 return hptoint(hp, Nside);
948 }
949
radectohealpix(double ra,double dec,int Nside)950 int radectohealpix(double ra, double dec, int Nside) {
951 return xyztohealpix(radec2x(ra,dec), radec2y(ra,dec), radec2z(ra,dec), Nside);
952 }
953
radectohealpixlf(double ra,double dec,int Nside,double * dx,double * dy)954 int64_t radectohealpixlf(double ra, double dec, int Nside, double* dx, double* dy) {
955 return xyztohealpixlf(radec2x(ra,dec), radec2y(ra,dec), radec2z(ra,dec), Nside, dx, dy);
956 }
957
radectohealpixl(double ra,double dec,int Nside)958 int64_t radectohealpixl(double ra, double dec, int Nside) {
959 return xyztohealpixl(radec2x(ra,dec), radec2y(ra,dec), radec2z(ra,dec), Nside);
960 }
961
radectohealpixf(double ra,double dec,int Nside,double * dx,double * dy)962 int radectohealpixf(double ra, double dec, int Nside, double* dx, double* dy) {
963 return xyztohealpixf(radec2x(ra,dec), radec2y(ra,dec), radec2z(ra,dec),
964 Nside, dx, dy);
965 }
966
radecdegtohealpix(double ra,double dec,int Nside)967 Const int radecdegtohealpix(double ra, double dec, int Nside) {
968 return radectohealpix(deg2rad(ra), deg2rad(dec), Nside);
969 }
970
radecdegtohealpixl(double ra,double dec,int Nside)971 Const int64_t radecdegtohealpixl(double ra, double dec, int Nside) {
972 return radectohealpixl(deg2rad(ra), deg2rad(dec), Nside);
973 }
974
radecdegtohealpixlf(double ra,double dec,int Nside,double * dx,double * dy)975 int64_t radecdegtohealpixlf(double ra, double dec, int Nside, double* dx, double* dy) {
976 return radectohealpixlf(deg2rad(ra), deg2rad(dec), Nside, dx, dy);
977 }
978
radecdegtohealpixf(double ra,double dec,int Nside,double * dx,double * dy)979 int radecdegtohealpixf(double ra, double dec, int Nside, double* dx, double* dy) {
980 return radectohealpixf(deg2rad(ra), deg2rad(dec), Nside, dx, dy);
981 }
982
xyzarrtohealpix(const double * xyz,int Nside)983 int xyzarrtohealpix(const double* xyz, int Nside) {
984 return xyztohealpix(xyz[0], xyz[1], xyz[2], Nside);
985 }
986
xyzarrtohealpixl(const double * xyz,int Nside)987 int64_t xyzarrtohealpixl(const double* xyz, int Nside) {
988 return xyztohealpixl(xyz[0], xyz[1], xyz[2], Nside);
989 }
990
xyzarrtohealpixf(const double * xyz,int Nside,double * p_dx,double * p_dy)991 int xyzarrtohealpixf(const double* xyz,int Nside, double* p_dx, double* p_dy) {
992 return xyztohealpixf(xyz[0], xyz[1], xyz[2], Nside, p_dx, p_dy);
993 }
994
hp_to_xyz(hp_t * hp,int Nside,double dx,double dy,double * rx,double * ry,double * rz)995 static void hp_to_xyz(hp_t* hp, int Nside,
996 double dx, double dy,
997 double* rx, double *ry, double *rz) {
998 int chp;
999 anbool equatorial = TRUE;
1000 double zfactor = 1.0;
1001 int xp, yp;
1002 double x, y, z;
1003 double pi = M_PI, phi;
1004 double rad;
1005
1006 hp_decompose(hp, &chp, &xp, &yp);
1007
1008 // this is x,y position in the healpix reference frame
1009 x = xp + dx;
1010 y = yp + dy;
1011
1012 if (isnorthpolar(chp)) {
1013 if ((x + y) > Nside) {
1014 equatorial = FALSE;
1015 zfactor = 1.0;
1016 }
1017 }
1018 if (issouthpolar(chp)) {
1019 if ((x + y) < Nside) {
1020 equatorial = FALSE;
1021 zfactor = -1.0;
1022 }
1023 }
1024
1025 if (equatorial) {
1026 double zoff=0;
1027 double phioff=0;
1028 x /= (double)Nside;
1029 y /= (double)Nside;
1030
1031 if (chp <= 3) {
1032 // north
1033 phioff = 1.0;
1034 } else if (chp <= 7) {
1035 // equator
1036 zoff = -1.0;
1037 chp -= 4;
1038 } else if (chp <= 11) {
1039 // south
1040 phioff = 1.0;
1041 zoff = -2.0;
1042 chp -= 8;
1043 } else {
1044 // should never get here
1045 assert(0);
1046 }
1047
1048 z = 2.0/3.0*(x + y + zoff);
1049 phi = pi/4*(x - y + phioff + 2*chp);
1050
1051 } else {
1052 /*
1053 Rearrange eqns (19) and (20) to find phi_t in terms of x,y.
1054
1055 y = Ns - k in eq(19)
1056 x - Ns - k in eq(20)
1057
1058 (Ns - y)^2 / (Ns - x)^2 = (2 phi_t)^2 / (2 phi_t - pi)^2
1059
1060 Recall than y<=Ns, x<=Ns and 0<=phi_t<pi/2, so we can choose the
1061 root we want by taking square roots:
1062
1063 (Ns - y) (pi - 2 phi_t) = 2 phi_t (Ns - x)
1064 (Ns - y) pi = 2 phi_t (Ns - x + Ns - y)
1065 phi_t = pi (Ns-y) / (2 (Ns - x) + (Ns - y))
1066 */
1067 double phi_t;
1068
1069 if (zfactor == -1.0) {
1070 swap_double(&x, &y);
1071 x = (Nside - x);
1072 y = (Nside - y);
1073 }
1074
1075 if (y == Nside && x == Nside)
1076 phi_t = 0.0;
1077 else
1078 phi_t = pi * (Nside-y) / (2.0 * ((Nside-x) + (Nside-y)));
1079
1080 if (phi_t < pi/4.) {
1081 z = 1.0 - mysquare(pi * (Nside - x) / ((2.0 * phi_t - pi) * Nside)) / 3.0;
1082 } else {
1083 z = 1.0 - mysquare(pi * (Nside - y) / (2.0 * phi_t * Nside)) / 3.0;
1084 }
1085 assert(0.0 <= fabs(z) && fabs(z) <= 1.0);
1086 z *= zfactor;
1087 assert(0.0 <= fabs(z) && fabs(z) <= 1.0);
1088
1089 // The big healpix determines the phi offset
1090 if (issouthpolar(chp))
1091 phi = pi/2.0* (chp-8) + phi_t;
1092 else
1093 phi = pi/2.0 * chp + phi_t;
1094 }
1095
1096 if (phi < 0.0)
1097 phi += 2*pi;
1098
1099 rad = sqrt(1.0 - z*z);
1100 *rx = rad * cos(phi);
1101 *ry = rad * sin(phi);
1102 *rz = z;
1103 }
1104
healpix_to_xyz(int ihp,int Nside,double dx,double dy,double * px,double * py,double * pz)1105 void healpix_to_xyz(int ihp, int Nside,
1106 double dx, double dy,
1107 double* px, double *py, double *pz) {
1108 hp_t hp;
1109 inttohp(ihp, &hp, Nside);
1110 hp_to_xyz(&hp, Nside, dx, dy, px, py, pz);
1111 }
1112
healpixl_to_radecdeg(int64_t ihp,int Nside,double dx,double dy,double * ra,double * dec)1113 void healpixl_to_radecdeg(int64_t ihp, int Nside, double dx, double dy,
1114 double* ra, double* dec) {
1115 hp_t hp;
1116 double xyz[3];
1117 intltohp(ihp, &hp, Nside);
1118 hp_to_xyz(&hp, Nside, dx, dy, xyz, xyz+1, xyz+2);
1119 xyzarr2radecdeg(xyz, ra, dec);
1120 }
1121
healpix_to_xyzarr(int ihp,int Nside,double dx,double dy,double * xyz)1122 void healpix_to_xyzarr(int ihp, int Nside,
1123 double dx, double dy,
1124 double* xyz) {
1125 hp_t hp;
1126 inttohp(ihp, &hp, Nside);
1127 hp_to_xyz(&hp, Nside, dx, dy, xyz, xyz+1, xyz+2);
1128 }
1129
healpix_to_radec(int hp,int Nside,double dx,double dy,double * ra,double * dec)1130 void healpix_to_radec(int hp, int Nside,
1131 double dx, double dy,
1132 double* ra, double* dec) {
1133 double xyz[3];
1134 healpix_to_xyzarr(hp, Nside, dx, dy, xyz);
1135 xyzarr2radec(xyz, ra, dec);
1136 }
1137
healpix_to_radecdeg(int hp,int Nside,double dx,double dy,double * ra,double * dec)1138 void healpix_to_radecdeg(int hp, int Nside,
1139 double dx, double dy,
1140 double* ra, double* dec) {
1141 double xyz[3];
1142 healpix_to_xyzarr(hp, Nside, dx, dy, xyz);
1143 xyzarr2radecdeg(xyz, ra, dec);
1144 }
1145
healpix_to_radecarr(int hp,int Nside,double dx,double dy,double * radec)1146 void healpix_to_radecarr(int hp, int Nside,
1147 double dx, double dy,
1148 double* radec) {
1149 double xyz[3];
1150 healpix_to_xyzarr(hp, Nside, dx, dy, xyz);
1151 xyzarr2radec(xyz, radec, radec+1);
1152 }
1153
healpix_to_radecdegarr(int hp,int Nside,double dx,double dy,double * radec)1154 void healpix_to_radecdegarr(int hp, int Nside,
1155 double dx, double dy,
1156 double* radec) {
1157 double xyz[3];
1158 healpix_to_xyzarr(hp, Nside, dx, dy, xyz);
1159 xyzarr2radecdeg(xyz, radec, radec+1);
1160 }
1161
1162 struct neighbour_dirn {
1163 double x, y;
1164 double dx, dy;
1165 };
1166
healpix_get_neighbours_within_range_radec(double ra,double dec,double radius,int * healpixes,int Nside)1167 int healpix_get_neighbours_within_range_radec(double ra, double dec, double radius,
1168 int* healpixes, int Nside) {
1169 double xyz[3];
1170 double r;
1171 radecdeg2xyzarr(ra, dec, xyz);
1172 r = deg2dist(radius);
1173 return healpix_get_neighbours_within_range(xyz, r, healpixes, Nside);
1174 }
1175
healpix_get_neighbours_within_range(double * xyz,double range,int * out_healpixes,int Nside)1176 int healpix_get_neighbours_within_range(double* xyz, double range, int* out_healpixes,
1177 int Nside) {
1178 int hp;
1179 int i,j;
1180 double fx, fy;
1181 int nhp = 0;
1182
1183 // HACK -- temp array to avoid cleverly avoiding duplicates
1184 int healpixes[100];
1185
1186 //assert(Nside > 0);
1187 if (Nside <= 0) {
1188 logerr("healpix_get_neighbours_within_range: Nside must be > 0.\n");
1189 return -1;
1190 }
1191
1192 hp = xyzarrtohealpixf(xyz, Nside, &fx, &fy);
1193 healpixes[nhp] = hp;
1194 nhp++;
1195
1196 {
1197 struct neighbour_dirn dirs[] = {
1198 // edges
1199 { fx, 0, 0, -1 },
1200 { fx, 1, 0, 1 },
1201 { 0 , fy,-1, 0 },
1202 { 1 , fy, 1, 0 },
1203 // bottom corner
1204 { 0, 0, -1, 1 },
1205 { 0, 0, -1, 0 },
1206 { 0, 0, -1, -1 },
1207 { 0, 0, 0, -1 },
1208 { 0, 0, 1, -1 },
1209 // right corner
1210 { 1, 0, 1, 1 },
1211 { 1, 0, 1, 0 },
1212 { 1, 0, 1, -1 },
1213 { 1, 0, 0, -1 },
1214 { 1, 0, -1, -1 },
1215 // left corner
1216 { 0, 1, 1, 1 },
1217 { 0, 1, 0, 1 },
1218 { 0, 1, -1, 1 },
1219 { 0, 1, -1, 0 },
1220 { 0, 1, -1, -1 },
1221 // top corner
1222 { 1, 1, -1, 1 },
1223 { 1, 1, 0, 1 },
1224 { 1, 1, 1, 1 },
1225 { 1, 1, 1, 0 },
1226 { 1, 1, 1, -1 },
1227 };
1228 int ndirs = sizeof(dirs) / sizeof(struct neighbour_dirn);
1229
1230 double ptx, pty, ptdx, ptdy;
1231 int pthp;
1232
1233 for (i=0; i<ndirs; i++) {
1234 double pt[3];
1235 double ptstepx[3];
1236 double ptstepy[3];
1237 double across[3];
1238 double step = 1e-3;
1239 double d2;
1240 double stepdirx, stepdiry;
1241 struct neighbour_dirn* dir = dirs+i;
1242 ptx = dir->x;
1243 pty = dir->y;
1244 ptdx = dir->dx;
1245 ptdy = dir->dy;
1246
1247 // pt = point on the edge nearest to the query point.
1248 // FIXME -- check that this is true, esp in the polar regions!
1249 healpix_to_xyzarr(hp, Nside, ptx, pty, pt);
1250 d2 = distsq(pt, xyz, 3);
1251
1252 // delta vector should be outside the healpix
1253 assert((ptx+step*ptdx < 0) ||
1254 (ptx+step*ptdx > 1) ||
1255 (pty+step*ptdy < 0) ||
1256 (pty+step*ptdy > 1));
1257
1258 if (d2 > range*range)
1259 continue;
1260
1261 // compute dx and dy directions that are toward the interior of
1262 // the healpix.
1263 stepdirx = (ptx < step) ? 1 : -1;
1264 stepdiry = (pty < step) ? 1 : -1;
1265
1266 // take steps in those directions.
1267 healpix_to_xyzarr(hp, Nside, ptx + stepdirx * step, pty, ptstepx);
1268 healpix_to_xyzarr(hp, Nside, ptx, pty + stepdiry * step, ptstepy);
1269
1270 // convert the steps into dx,dy vectors.
1271 for (j=0; j<3; j++) {
1272 ptstepx[j] = stepdirx * (ptstepx[j] - pt[j]);
1273 ptstepy[j] = stepdiry * (ptstepy[j] - pt[j]);
1274 }
1275
1276 // take a small step in the specified direction.
1277 for (j=0; j<3; j++)
1278 across[j] = pt[j] + ptdx * ptstepx[j] + ptdy * ptstepy[j];
1279
1280 // see which healpix is at the end of the step.
1281 normalize_3(across);
1282 pthp = xyzarrtohealpix(across, Nside);
1283
1284 healpixes[nhp] = pthp;
1285 nhp++;
1286 }
1287 }
1288
1289 // Remove duplicates...
1290 for (i=0; i<nhp; i++) {
1291 for (j=i+1; j<nhp; j++) {
1292 if (healpixes[i] == healpixes[j]) {
1293 int k;
1294 for (k=j+1; k<nhp; k++)
1295 healpixes[k-1] = healpixes[k];
1296 nhp--;
1297 i=-1;
1298 break;
1299 }
1300 }
1301 }
1302
1303 for (i=0; i<nhp; i++)
1304 out_healpixes[i] = healpixes[i];
1305
1306 return nhp;
1307 }
1308
healpix_distance_to_xyz(int hp,int Nside,const double * xyz,double * closestxyz)1309 double healpix_distance_to_xyz(int hp, int Nside, const double* xyz,
1310 double* closestxyz) {
1311 int thehp;
1312 // corners
1313 double cdx[4], cdy[4];
1314 double cdists[4];
1315 int corder[4];
1316 int i;
1317 double dxA,dxB,dyA,dyB;
1318 double dxmid, dymid;
1319 double dist2A, dist2B;
1320 double midxyz[3];
1321 double dist2mid = 0.0;
1322
1323 double EPS = 1e-16;
1324
1325 // If the point is actually inside the healpix, dist = 0.
1326 thehp = xyzarrtohealpix(xyz, Nside);
1327 if (thehp == hp) {
1328 if (closestxyz)
1329 memcpy(closestxyz, xyz, 3*sizeof(double));
1330 return 0;
1331 }
1332
1333 // Find two nearest corners.
1334 for (i=0; i<4; i++) {
1335 double cxyz[3];
1336 cdx[i] = i/2;
1337 cdy[i] = i%2;
1338 healpix_to_xyzarr(hp, Nside, cdx[i], cdy[i], cxyz);
1339 cdists[i] = distsq(xyz, cxyz, 3);
1340 }
1341 permutation_init(corder, 4);
1342 permuted_sort(cdists, sizeof(double), compare_doubles_asc, corder, 4);
1343 // now "corder" [0] and [1] are the indices of the two nearest corners.
1344
1345 // Binary search for closest dx,dy.
1346 dxA = cdx[corder[0]];
1347 dyA = cdy[corder[0]];
1348 dist2A = cdists[corder[0]];
1349 dxB = cdx[corder[1]];
1350 dyB = cdy[corder[1]];
1351 dist2B = cdists[corder[1]];
1352 // the closest two points should share an edge... unless we're in some
1353 // weird configuration like the opposite side of the sphere.
1354 if (!((dxA == dxB) || (dyA == dyB))) {
1355 // weird configuration -- return the closest corner.
1356 if (closestxyz)
1357 healpix_to_xyzarr(hp, Nside, cdx[corder[0]], cdy[corder[0]], closestxyz);
1358 return distsq2deg(cdists[corder[0]]);
1359 }
1360 assert(dxA == dxB || dyA == dyB);
1361 assert(dist2A <= dist2B);
1362
1363 while (1) {
1364 dxmid = (dxA + dxB) / 2.0;
1365 dymid = (dyA + dyB) / 2.0;
1366 // converged to EPS?
1367 if ((dxA != dxB && (fabs(dxmid - dxA) < EPS || fabs(dxmid - dxB) < EPS)) ||
1368 (dyA != dyB && (fabs(dymid - dyA) < EPS || fabs(dymid - dyB) < EPS)))
1369 break;
1370 healpix_to_xyzarr(hp, Nside, dxmid, dymid, midxyz);
1371 dist2mid = distsq(xyz, midxyz, 3);
1372 //printf(" dx,dy (%g,%g) %g (%g,%g) %g (%g,%g) %g\n", dxA, dyA, dist2A, dxmid, dymid, dist2mid, dxB, dyB, dist2B);
1373 if ((dist2mid >= dist2A) && (dist2mid >= dist2B))
1374 break;
1375 if (dist2A < dist2B) {
1376 dist2B = dist2mid;
1377 dxB = dxmid;
1378 dyB = dymid;
1379 } else {
1380 dist2A = dist2mid;
1381 dxA = dxmid;
1382 dyA = dymid;
1383 }
1384 }
1385
1386 // Check whether endpoint A is actually closer.
1387 dist2A = cdists[corder[0]];
1388 if (dist2A < dist2mid) {
1389 dxA = cdx[corder[0]];
1390 dyA = cdy[corder[0]];
1391 healpix_to_xyzarr(hp, Nside, dxA, dyA, midxyz);
1392 dist2mid = dist2A;
1393 }
1394
1395 if (closestxyz)
1396 memcpy(closestxyz, midxyz, 3*sizeof(double));
1397 /*{
1398 double ra,dec;
1399 double close[2];
1400 xyzarr2radecdeg(xyz, &ra,&dec);
1401 xyzarr2radecdegarr(midxyz, close);
1402 printf("healpix_distance_to_xyz: %.4f,%.4f, hp %i -> closest %.4f, %.4f -> dist %g deg\n",
1403 ra, dec, hp, close[0],close[1], distsq2deg(dist2mid));
1404 }*/
1405 return distsq2deg(dist2mid);
1406 }
1407
healpix_distance_to_radec(int hp,int Nside,double ra,double dec,double * closestradec)1408 double healpix_distance_to_radec(int hp, int Nside, double ra, double dec,
1409 double* closestradec) {
1410 double xyz[3];
1411 double closestxyz[3];
1412 double dist;
1413 radecdeg2xyzarr(ra, dec, xyz);
1414 dist = healpix_distance_to_xyz(hp, Nside, xyz, closestxyz);
1415 if (closestradec)
1416 xyzarr2radecdegarr(closestxyz, closestradec);
1417 return dist;
1418 }
1419
healpix_within_range_of_radec(int hp,int Nside,double ra,double dec,double radius)1420 int healpix_within_range_of_radec(int hp, int Nside, double ra, double dec,
1421 double radius) {
1422 // the trivial implementation...
1423 return (healpix_distance_to_radec(hp, Nside, ra, dec, NULL) <= radius);
1424 }
healpix_within_range_of_xyz(int hp,int Nside,const double * xyz,double radius)1425 int healpix_within_range_of_xyz(int hp, int Nside, const double* xyz,
1426 double radius) {
1427 return (healpix_distance_to_xyz(hp, Nside, xyz, NULL) <= radius);
1428 }
1429
1430
healpix_radec_bounds(int hp,int nside,double * p_ralo,double * p_rahi,double * p_declo,double * p_dechi)1431 void healpix_radec_bounds(int hp, int nside,
1432 double* p_ralo, double* p_rahi,
1433 double* p_declo, double* p_dechi) {
1434 // corners!
1435 double ralo,rahi,declo,dechi;
1436 double ra,dec;
1437 double dx, dy;
1438 ralo = declo = HUGE_VAL;
1439 rahi = dechi = -HUGE_VAL;
1440 for (dy=0; dy<2; dy+=1.0) {
1441 for (dx=0; dx<2; dx+=1.0) {
1442 healpix_to_radecdeg(hp, nside, dx, dy, &ra, &dec);
1443 // FIXME -- wrap-around.
1444 ralo = MIN(ra, ralo);
1445 rahi = MAX(ra, rahi);
1446 declo = MIN(dec, declo);
1447 dechi = MAX(dec, dechi);
1448 }
1449 }
1450 if (p_ralo) *p_ralo = ralo;
1451 if (p_rahi) *p_rahi = rahi;
1452 if (p_declo) *p_declo = declo;
1453 if (p_dechi) *p_dechi = dechi;
1454 }
1455