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 #ifndef HEALPIX_H 7 #define HEALPIX_H 8 9 #include <sys/types.h> 10 #include <stdint.h> 11 12 #include "astrometry/keywords.h" 13 14 //#undef Const 15 //#define Const 16 17 /** 18 The HEALPix paper is here: 19 http://iopscience.iop.org/0004-637X/622/2/759/pdf/0004-637X_622_2_759.pdf 20 See: 21 http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=2005ApJ...622..759G&db_key=AST&high=41069202cf02947 22 */ 23 24 /** 25 In this documentation we talk about "base healpixes": these are the big, 26 top-level healpixes. There are 12 of these, with indices [0, 11]. 27 28 We say "fine healpixes" or "healpixes" or "pixels" when we mean the fine- 29 scale healpixes; there are Nside^2 of these in each base healpix, 30 for a total of 12*Nside^2, indexed from zero. 31 */ 32 33 /** 34 Some notes about the different indexing schemes: 35 36 The healpix paper discusses two different ways to number healpixes, and 37 there is a third way, which we prefer, which is (in my opinion) more 38 sensible and easy. 39 40 41 -RING indexing. Healpixes are numbered first in order of decreasing DEC, 42 then in order of increasing RA of the center of the pixel, ie: 43 44 . 0 1 2 3 45 . 4 5 6 7 8 9 10 11 46 . 12 13 14 15 16 17 18 19 47 . 20 21 22 23 24 25 26 27 48 . 28 29 30 31 32 33 34 35 49 . 36 37 38 39 40 41 42 43 50 . 44 45 46 47 51 52 Note that 12, 20 and 28 are part of base healpix 4, as is 27; it "wraps 53 around". 54 55 The RING index can be decomposed into the "ring number" and the index 56 within the ring (called "longitude index"). Note that different rings 57 contain different numbers of healpixes. Also note that the ring number 58 starts from 1, but the longitude index starts from zero. 59 60 61 -NESTED indexing. This only works for Nside parameters that are powers of 62 two. This scheme is hierarchical in the sense that each pair of bits of 63 the index tells you where the pixel center is to finer and finer 64 resolution. It doesn't really show with Nside=2, but here it is anyway: 65 66 . 3 7 11 15 67 . 2 1 6 5 10 9 14 13 68 . 19 0 23 4 27 8 31 12 69 . 17 22 21 26 25 30 29 18 70 . 16 35 20 39 24 43 28 47 71 . 34 33 38 37 42 41 46 45 72 . 32 36 40 44 73 74 Note that all the base healpixes have the same pattern; they're just 75 offset by factors of Nside^2. 76 77 Here's a zoom-in of the first base healpix, turned 45 degrees to the 78 right, for Nside=4: 79 80 . 10 11 14 15 81 . 8 9 12 13 82 . 2 3 6 7 83 . 0 1 4 5 84 85 Note that the bottom-left block of 4 have the smallest values, and within 86 that the bottom-left corner has the smallest value, followed by the 87 bottom-right, top-left, then top-right. 88 89 The NESTED index can't be decomposed into 'orthogonal' directions. 90 91 92 -XY indexing. This is arguably the most natural, at least for the 93 internal usage of the healpix code. Within each base healpix, the 94 healpixes are numbered starting with 0 for the southmost pixel, then 95 increasing first in the "y" (north-west), then in the "x" (north-east) 96 direction. In other words, within each base healpix there is a grid 97 and we number the pixels "lexicographically" (mod a 135 degree turn). 98 99 . 3 7 11 15 100 . 1 2 5 6 9 10 13 14 101 . 19 0 23 4 27 8 31 12 102 . 18 21 22 25 26 29 30 17 103 . 16 35 20 39 24 43 28 47 104 . 33 34 37 38 41 42 45 46 105 . 32 36 40 44 106 107 Zooming in on the first base healpix, turning 45 degrees to the right, 108 for Nside=4 we get: 109 110 . 3 7 11 15 111 . 2 6 10 14 112 . 1 5 9 13 113 . 0 4 8 12 114 115 Notice that the numbers first increase from bottom to top (y), then left to 116 right (x). 117 118 The XY indexing can be decomposed into 'x' and 'y' coordinates 119 (in case that wasn't obvious), where the above figure becomes (x,y): 120 121 . (0,3) (1,3) (2,3) (3,3) 122 . (0,2) (1,2) (2,2) (3,2) 123 . (0,1) (1,1) (2,1) (3,1) 124 . (0,0) (1,0) (2,0) (3,0) 125 126 Note that "x" increases in the north-east direction, and "y" increases in 127 the north-west direction. 128 129 The major advantage to this indexing scheme is that it extends to 130 fractional coordinates in a natural way: it is meaningful to talk about 131 the position (x,y) = (0.25, 0.6) and you can compute its position. 132 133 134 135 136 In this code, all healpix indexing uses the XY scheme. If you want to 137 use the other schemes you will have to use the conversion routines: 138 . healpix_xy_to_ring 139 . healpix_ring_to_xy 140 . healpix_xy_to_nested 141 . healpix_nested_to_xy 142 */ 143 144 // The maximum healpix Nside that leads to int-sized healpix indices. 145 // 12 * (13377+1)^2 > 2^31 (since we use signed ints) 146 // This corresponds to about 16 arcsec side length. 147 #define HP_MAX_INT_NSIDE 13377 148 149 /** 150 Converts a healpix index from the XY scheme to the RING scheme. 151 */ 152 Const int healpix_xy_to_ring(int hp, int Nside); 153 154 /** 155 Converts a healpix index from the RING scheme to the XY scheme. 156 */ 157 Const int healpix_ring_to_xy(int ring_index, int Nside); 158 159 /** 160 Converts a healpix index from the XY scheme to the NESTED scheme. 161 */ 162 Const int healpix_xy_to_nested(int hp, int Nside); 163 164 /** 165 Converts a healpix index from the NESTED scheme to the XY scheme. 166 */ 167 Const int healpix_nested_to_xy(int nested_index, int Nside); 168 169 /** 170 Decomposes a RING index into the "ring number" (each ring contain 171 healpixels of equal latitude) and "longitude index". Pixels within a 172 ring have longitude index starting at zero for the first pixel with 173 RA >= 0. Different rings contain different numbers of healpixels. 174 */ 175 void healpix_decompose_ring(int ring_index, int Nside, 176 int* p_ring_number, int* p_longitude_index); 177 178 /** 179 Composes a RING index given the "ring number" and "longitude index". 180 181 Does NOT check that the values are legal! Garbage in, garbage out. 182 */ 183 Const int healpix_compose_ring(int ring, int longind, int Nside); 184 185 /** 186 Decomposes an XY index into the "base healpix" and "x" and "y" coordinates 187 within that healpix. 188 */ 189 void healpix_decompose_xy(int finehp, int* bighp, int* x, int* y, int Nside); 190 191 void healpix_decompose_xyl(int64_t finehp, int* bighp, int* x, int* y, int Nside); 192 193 /** 194 Composes an XY index given the "base healpix" and "x" and "y" coordinates 195 within that healpix. 196 */ 197 Const int healpix_compose_xy(int bighp, int x, int y, int Nside); 198 199 Const int64_t healpix_compose_xyl(int bighp, int x, int y, int Nside); 200 201 /** 202 Given (x,y) coordinates of resolution "nside" within a base-level 203 healpixel, and an output resolution "outnside", returns the output 204 (x,y) coordinates at the output resolution. 205 */ 206 void healpix_convert_xy_nside(int x, int y, int nside, int outnside, 207 int* outx, int* outy); 208 209 /** 210 Given a healpix index (in the XY scheme) of resolution "nside", and 211 an output resolution "outnside", returns the healpix index at the 212 output resolution. 213 */ 214 void healpix_convert_nside(int hp, int nside, int outnside, int* outhp); 215 216 /** 217 Converts (RA, DEC) coordinates (in radians) to healpix index. 218 */ 219 Const int radectohealpix(double ra, double dec, int Nside); 220 221 int radectohealpixf(double ra, double dec, int Nside, double* dx, double* dy); 222 223 Const int64_t radectohealpixl(double ra, double dec, int Nside); 224 225 int64_t radectohealpixlf(double ra, double dec, int Nside, double* dx, double* dy); 226 227 /** 228 Converts (RA, DEC) coordinates (in degrees) to healpix index. 229 */ 230 Const int radecdegtohealpix(double ra, double dec, int Nside); 231 232 int radecdegtohealpixf(double ra, double dec, int Nside, double* dx, double* dy); 233 234 Const int64_t radecdegtohealpixl(double ra, double dec, int Nside); 235 236 int64_t radecdegtohealpixlf(double ra, double dec, int Nside, double* dx, double* dy); 237 238 /** 239 Converts (x,y,z) coordinates on the unit sphere into a healpix index. 240 */ 241 Const int xyztohealpix(double x, double y, double z, int Nside); 242 243 Const int64_t xyztohealpixl(double x, double y, double z, int Nside); 244 245 int xyztohealpixf(double x, double y, double z, int Nside, 246 double* p_dx, double* p_dy); 247 248 int64_t xyztohealpixlf(double x, double y, double z, int Nside, 249 double* p_dx, double* p_dy); 250 251 /** 252 Converts (x,y,z) coordinates (stored in an array) on the unit sphere into 253 a healpix index. 254 */ 255 int xyzarrtohealpix(const double* xyz, int Nside); 256 257 int64_t xyzarrtohealpixl(const double* xyz, int Nside); 258 259 int xyzarrtohealpixf(const double* xyz,int Nside, double* p_dx, double* p_dy); 260 261 /** 262 Converts a healpix index, plus fractional offsets (dx,dy), into (x,y,z) 263 coordinates on the unit sphere. (dx,dy) must be in [0, 1]. (0.5, 0.5) 264 is the center of the healpix. (0,0) is the southernmost corner, (1,1) is 265 the northernmost corner, (1,0) is the easternmost, and (0,1) the 266 westernmost. 267 */ 268 void healpix_to_xyz(int hp, int Nside, double dx, double dy, 269 double* p_x, double *p_y, double *p_z); 270 271 /** 272 Same as healpix_to_xyz, but (x,y,z) are stored in an array. 273 */ 274 void healpix_to_xyzarr(int hp, int Nside, double dx, double dy, 275 double* xyz); 276 277 /** 278 Same as healpix_to_xyz, but returns (RA,DEC) in radians. 279 */ 280 void healpix_to_radec(int hp, int Nside, double dx, double dy, 281 double* ra, double* dec); 282 283 void healpix_to_radecdeg(int hp, int Nside, double dx, double dy, 284 double* ra, double* dec); 285 286 void healpixl_to_radecdeg(int64_t hp, int Nside, double dx, double dy, 287 double* ra, double* dec); 288 289 /** 290 Same as healpix_to_radec, but (RA,DEC) are stored in an array. 291 */ 292 void healpix_to_radecarr(int hp, int Nside, double dx, double dy, 293 double* radec); 294 295 void healpix_to_radecdegarr(int hp, int Nside, double dx, double dy, 296 double* radec); 297 298 /** 299 Computes the approximate side length of a healpix, in arcminutes. 300 */ 301 Const double healpix_side_length_arcmin(int Nside); 302 303 /** 304 Computes the approximate Nside you need to get healpixes with side 305 length about "arcmin" arcminutes. (inverse of 306 healpix_side_length_arcmin) 307 */ 308 double healpix_nside_for_side_length_arcmin(double arcmin); 309 310 /** 311 Finds the healpixes neighbouring the given healpix, placing them in the 312 array "neighbour". Returns the number of neighbours. You must ensure 313 that "neighbour" has 8 elements. 314 315 Healpixes in the interior of a large healpix will have eight neighbours; 316 pixels near the edges can have fewer. 317 */ 318 int healpix_get_neighbours(int hp, int* neighbours, int Nside); 319 320 /** 321 Same as above, but for Nsides big enough that it overflows 32-bit int. 322 */ 323 int healpix_get_neighboursl(int64_t pix, int64_t* neighbour, int Nside); 324 325 /** 326 Finds the healpixes containing and neighbouring the given xyz 327 position which are within distance 'range' (in units of distance of 328 the unit sphere). Places the results in 'healpixes', which must have 329 at least 9 elements. Returns the number of 'healpixes' set. 330 331 Returns -1 if "Nside" < 0. 332 */ 333 int healpix_get_neighbours_within_range(double* xyz, double range, int* healpixes, 334 int Nside); 335 336 /** 337 Same as above, but RA,Dec,radius in degrees. 338 */ 339 int healpix_get_neighbours_within_range_radec(double ra, double dec, double radius, 340 int* healpixes, int Nside); 341 342 /** 343 Returns the minimum distance (in degrees) between the given healpix 344 and the given RA,Dec (in degrees). 345 */ 346 double healpix_distance_to_radec(int hp, int Nside, double ra, double dec, 347 double* closestradec); 348 349 /** 350 Returns the minimum distance (in degrees) between the given healpix 351 and the given xyz (point on unit sphere). 352 */ 353 double healpix_distance_to_xyz(int hp, int Nside, const double* xyz, 354 double* closestxyz); 355 356 /** 357 Returns true if the closest distance between the given healpix and 358 the given RA,Dec (in degrees) is less than then given radius (in degrees). 359 */ 360 int healpix_within_range_of_radec(int hp, int Nside, double ra, double dec, 361 double radius); 362 int healpix_within_range_of_xyz(int hp, int Nside, const double* xyz, 363 double radius); 364 365 366 /** 367 Computes the RA,Dec bounding-box of the given healpix. Results are 368 in degrees. RA may be wacky for healpixes spanning RA=0. 369 */ 370 void healpix_radec_bounds(int hp, int nside, 371 double* ralo, double* rahi, 372 double* declo, double* dechi); 373 374 #endif 375