1 /* 2 * This file is part of Healpix_cxx. 3 * 4 * Healpix_cxx 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 * Healpix_cxx 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 Healpix_cxx; 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 /* 22 * Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik 23 * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt 24 * (DLR). 25 */ 26 27 /*! \file healpix_base.h 28 * Copyright (C) 2003-2017 Max-Planck-Society 29 * \author Martin Reinecke 30 */ 31 32 #ifndef HEALPIX_BASE_H 33 #define HEALPIX_BASE_H 34 35 #include <vector> 36 #include "healpix_tables.h" 37 #include "pointing.h" 38 #include "arr.h" 39 #include "rangeset.h" 40 41 template<typename I> struct Orderhelper__ {}; 42 template<> struct Orderhelper__<int> {enum{omax=13};}; 43 template<> struct Orderhelper__<int64> {enum{omax=29};}; 44 45 /*! Functionality related to the HEALPix pixelisation. */ 46 template<typename I> class T_Healpix_Base: public Healpix_Tables 47 { 48 protected: 49 /*! The order of the map; -1 for nonhierarchical map. */ 50 int order_; 51 /*! The N_side parameter of the map; 0 if not allocated. */ 52 I nside_; 53 I npface_, ncap_, npix_; 54 double fact1_, fact2_; 55 /*! The map's ordering scheme. */ 56 Healpix_Ordering_Scheme scheme_; 57 58 /*! Returns the number of the next ring to the north of \a z=cos(theta). 59 It may return 0; in this case \a z lies north of all rings. */ 60 inline I ring_above (double z) const; 61 void in_ring (I iz, double phi0, double dphi, rangeset<I> &pixset) const; 62 63 template<typename I2> void query_multidisc (const arr<vec3> &norm, 64 const arr<double> &rad, int fact, rangeset<I2> &pixset) const; 65 66 void query_multidisc_general (const arr<vec3> &norm, const arr<double> &rad, 67 bool inclusive, const std::vector<int> &cmds, rangeset<I> &pixset) const; 68 69 void query_strip_internal (double theta1, double theta2, bool inclusive, 70 rangeset<I> &pixset) const; 71 72 inline I spread_bits (int v) const; 73 inline int compress_bits (I v) const; 74 75 I xyf2nest(int ix, int iy, int face_num) const; 76 void nest2xyf(I pix, int &ix, int &iy, int &face_num) const; 77 I xyf2ring(int ix, int iy, int face_num) const; 78 void ring2xyf(I pix, int &ix, int &iy, int &face_num) const; 79 80 I loc2pix (double z, double phi, double sth, bool have_sth) const; 81 void pix2loc (I pix, double &z, double &phi, double &sth, bool &have_sth) 82 const; 83 84 void xyf2loc(double x, double y, int face, double &z, double &ph, 85 double &sth, bool &have_sth) const; 86 87 I nest_peano_helper (I pix, int dir) const; 88 89 typedef I (T_Healpix_Base::*swapfunc)(I pix) const; 90 91 public: 92 enum {order_max=Orderhelper__<I>::omax}; 93 94 /*! Calculates the map order from its \a N_side parameter. 95 Returns -1 if \a nside is not a power of 2. 96 \param nside the \a N_side parameter */ 97 static int nside2order (I nside); 98 /*! Calculates the \a N_side parameter from the number of pixels. 99 \param npix the number of pixels */ 100 static I npix2nside (I npix); 101 /*! Constructs an unallocated object. */ 102 T_Healpix_Base (); 103 /*! Constructs an object with a given \a order and the ordering 104 scheme \a scheme. */ 105 T_Healpix_Base (int order, Healpix_Ordering_Scheme scheme) 106 { Set (order, scheme); } 107 /*! Constructs an object with a given \a nside and the ordering 108 scheme \a scheme. The \a nside_dummy parameter must be set to 109 SET_NSIDE. */ 110 T_Healpix_Base (I nside, Healpix_Ordering_Scheme scheme, const nside_dummy) 111 { SetNside (nside, scheme); } 112 113 /*! Adjusts the object to \a order and \a scheme. */ 114 void Set (int order, Healpix_Ordering_Scheme scheme); 115 /*! Adjusts the object to \a nside and \a scheme. */ 116 void SetNside (I nside, Healpix_Ordering_Scheme scheme); 117 118 /*! Returns the z-coordinate of the ring \a ring. This also works 119 for the (not really existing) rings 0 and 4*nside. */ 120 double ring2z (I ring) const; 121 /*! Returns the number of the ring in which \a pix lies. */ 122 I pix2ring (I pix) const; 123 124 I xyf2pix(int ix, int iy, int face_num) const 125 { 126 return (scheme_==RING) ? 127 xyf2ring(ix,iy,face_num) : xyf2nest(ix,iy,face_num); 128 } 129 void pix2xyf(I pix, int &ix, int &iy, int &face_num) const 130 { 131 (scheme_==RING) ? 132 ring2xyf(pix,ix,iy,face_num) : nest2xyf(pix,ix,iy,face_num); 133 } 134 135 /*! Translates a pixel number from NEST to RING. */ 136 I nest2ring (I pix) const; 137 /*! Translates a pixel number from RING to NEST. */ 138 I ring2nest (I pix) const; 139 /*! Translates a pixel number from NEST to its Peano index. */ 140 I nest2peano (I pix) const; 141 /*! Translates a pixel number from its Peano index to NEST. */ 142 I peano2nest (I pix) const; 143 144 /*! Returns the number of the pixel which contains the angular coordinates 145 (\a z:=cos(theta), \a phi). 146 \note This method is inaccurate near the poles at high resolutions. */ 147 I zphi2pix (double z, double phi) const 148 { return loc2pix(z,phi,0.,false); } 149 150 /*! Returns the number of the pixel which contains the angular coordinates 151 \a ang. */ 152 I ang2pix (const pointing &ang) const 153 { 154 const double pi_=3.141592653589793238462643383279502884197; 155 planck_assert((ang.theta>=0)&&(ang.theta<=pi_),"invalid theta value"); 156 return ((ang.theta<0.01) || (ang.theta > 3.14159-0.01)) ? 157 loc2pix(cos(ang.theta),ang.phi,sin(ang.theta),true) : 158 loc2pix(cos(ang.theta),ang.phi,0.,false); 159 } 160 /*! Returns the number of the pixel which contains the vector \a vec 161 (\a vec is normalized if necessary). */ 162 I vec2pix (const vec3 &vec) const 163 { 164 double xl = 1./vec.Length(); 165 double phi = safe_atan2(vec.y,vec.x); 166 double nz = vec.z*xl; 167 if (std::abs(nz)>0.99) 168 return loc2pix (nz,phi,sqrt(vec.x*vec.x+vec.y*vec.y)*xl,true); 169 else 170 return loc2pix (nz,phi,0,false); 171 } 172 173 /*! Returns the angular coordinates (\a z:=cos(theta), \a phi) of the center 174 of the pixel with number \a pix. 175 \note This method is inaccurate near the poles at high resolutions. */ 176 void pix2zphi (I pix, double &z, double &phi) const 177 { 178 bool dum_b; 179 double dum_d; 180 pix2loc(pix,z,phi,dum_d,dum_b); 181 } 182 183 /*! Returns the angular coordinates of the center of the pixel with 184 number \a pix. */ 185 pointing pix2ang (I pix) const 186 { 187 double z, phi, sth; 188 bool have_sth; 189 pix2loc (pix,z,phi,sth,have_sth); 190 return have_sth ? pointing(atan2(sth,z),phi) : pointing(acos(z),phi); 191 } 192 /*! Returns the vector to the center of the pixel with number \a pix. */ 193 vec3 pix2vec (I pix) const 194 { 195 double z, phi, sth; 196 bool have_sth; 197 pix2loc (pix,z,phi,sth,have_sth); 198 if (have_sth) 199 return vec3(sth*cos(phi),sth*sin(phi),z); 200 else 201 { 202 vec3 res; 203 res.set_z_phi (z, phi); 204 return res; 205 } 206 } 207 /*! Returns the pixel number for this T_Healpix_Base corresponding to the 208 pixel number \a pix in \a b. 209 \note \a b.Nside()\%Nside() must be 0. */ 210 I pixel_import (I pix, const T_Healpix_Base &b) const 211 { 212 I ratio = b.nside_/nside_; 213 planck_assert(nside_*ratio==b.nside_,"bad nside ratio"); 214 int x, y, f; 215 b.pix2xyf(pix, x, y, f); 216 x/=ratio; y/=ratio; 217 return xyf2pix(x, y, f); 218 } 219 220 template<typename I2> void query_disc_internal (pointing ptg, double radius, 221 int fact, rangeset<I2> &pixset) const; 222 223 /*! Returns the range set of all pixels whose centers lie within the disk 224 defined by \a dir and \a radius. 225 \param dir the angular coordinates of the disk center 226 \param radius the radius (in radians) of the disk 227 \param pixset a \a rangeset object containing the indices of all pixels 228 whose centers lie inside the disk 229 \note This method is more efficient in the RING scheme. */ 230 void query_disc (pointing ptg, double radius, rangeset<I> &pixset) const; 231 /*! Returns the range set of all pixels whose centers lie within the disk 232 defined by \a dir and \a radius. 233 \param dir the angular coordinates of the disk center 234 \param radius the radius (in radians) of the disk 235 \note This method is more efficient in the RING scheme. */ 236 rangeset<I> query_disc (pointing ptg, double radius) const 237 { 238 rangeset<I> res; 239 query_disc(ptg, radius, res); 240 return res; 241 } 242 /*! Returns the range set of all pixels which overlap with the disk 243 defined by \a dir and \a radius. 244 \param dir the angular coordinates of the disk center 245 \param radius the radius (in radians) of the disk 246 \param pixset a \a rangeset object containing the indices of all pixels 247 overlapping with the disk. 248 \param fact The overlapping test will be done at the resolution 249 \a fact*nside. For NESTED ordering, \a fact must be a power of 2, 250 else it can be any positive integer. A typical choice would be 4. 251 \note This method may return some pixels which don't overlap with 252 the disk at all. The higher \a fact is chosen, the fewer false 253 positives are returned, at the cost of increased run time. 254 \note This method is more efficient in the RING scheme. */ 255 void query_disc_inclusive (pointing ptg, double radius, rangeset<I> &pixset, 256 int fact=1) const; 257 /*! Returns the range set of all pixels which overlap with the disk 258 defined by \a dir and \a radius. 259 \param dir the angular coordinates of the disk center 260 \param radius the radius (in radians) of the disk 261 \param fact The overlapping test will be done at the resolution 262 \a fact*nside. For NESTED ordering, \a fact must be a power of 2, 263 else it can be any positive integer. A typical choice would be 4. 264 \note This method may return some pixels which don't overlap with 265 the disk at all. The higher \a fact is chosen, the fewer false 266 positives are returned, at the cost of increased run time. 267 \note This method is more efficient in the RING scheme. */ 268 rangeset<I> query_disc_inclusive (pointing ptg, double radius, 269 int fact=1) const 270 { 271 rangeset<I> res; 272 query_disc_inclusive(ptg, radius, res, fact); 273 return res; 274 } 275 276 /*! \deprecated Please use the version based on \a rangeset */ 277 void query_disc (const pointing &dir, double radius, 278 std::vector<I> &listpix) const 279 { 280 rangeset<I> pixset; 281 query_disc(dir,radius,pixset); 282 pixset.toVector(listpix); 283 } 284 /*! \deprecated Please use the version based on \a rangeset */ 285 void query_disc_inclusive (const pointing &dir, double radius, 286 std::vector<I> &listpix, int fact=1) const 287 { 288 rangeset<I> pixset; 289 query_disc_inclusive(dir,radius,pixset,fact); 290 pixset.toVector(listpix); 291 } 292 293 template<typename I2> void query_polygon_internal 294 (const std::vector<pointing> &vertex, int fact, 295 rangeset<I2> &pixset) const; 296 297 /*! Returns a range set of pixels whose centers lie within the convex 298 polygon defined by the \a vertex array. 299 \param vertex array containing the vertices of the polygon. 300 \param pixset a \a rangeset object containing the indices of all pixels 301 whose centers lie inside the polygon 302 \note This method is more efficient in the RING scheme. */ 303 void query_polygon (const std::vector<pointing> &vertex, 304 rangeset<I> &pixset) const; 305 /*! Returns a range set of pixels whose centers lie within the convex 306 polygon defined by the \a vertex array. 307 \param vertex array containing the vertices of the polygon. 308 \note This method is more efficient in the RING scheme. */ 309 rangeset<I> query_polygon (const std::vector<pointing> &vertex) const 310 { 311 rangeset<I> res; 312 query_polygon(vertex, res); 313 return res; 314 } 315 316 /*! Returns a range set of pixels which overlap with the convex 317 polygon defined by the \a vertex array. 318 \param vertex array containing the vertices of the polygon. 319 \param pixset a \a rangeset object containing the indices of all pixels 320 overlapping with the polygon. 321 \param fact The overlapping test will be done at the resolution 322 \a fact*nside. For NESTED ordering, \a fact must be a power of 2, 323 else it can be any positive integer. A typical choice would be 4. 324 \note This method may return some pixels which don't overlap with 325 the polygon at all. The higher \a fact is chosen, the fewer false 326 positives are returned, at the cost of increased run time. 327 \note This method is more efficient in the RING scheme. */ 328 void query_polygon_inclusive (const std::vector<pointing> &vertex, 329 rangeset<I> &pixset, int fact=1) const; 330 /*! Returns a range set of pixels which overlap with the convex 331 polygon defined by the \a vertex array. 332 \param vertex array containing the vertices of the polygon. 333 \param fact The overlapping test will be done at the resolution 334 \a fact*nside. For NESTED ordering, \a fact must be a power of 2, 335 else it can be any positive integer. A typical choice would be 4. 336 \note This method may return some pixels which don't overlap with 337 the polygon at all. The higher \a fact is chosen, the fewer false 338 positives are returned, at the cost of increased run time. 339 \note This method is more efficient in the RING scheme. */ 340 rangeset<I> query_polygon_inclusive (const std::vector<pointing> &vertex, 341 int fact=1) const 342 { 343 rangeset<I> res; 344 query_polygon_inclusive(vertex, res, fact); 345 return res; 346 } 347 348 /*! Returns a range set of pixels whose centers lie within the colatitude 349 range defined by \a theta1 and \a theta2 (if \a inclusive==false), or 350 which overlap with this region (if \a inclusive==true). If 351 \a theta1<theta2, the region between both angles is considered, 352 otherwise the regions \a 0<theta<theta2 and \a theta1<theta<pi. 353 \param theta1 first colatitude 354 \param theta2 second colatitude 355 \param inclusive if \a false, return the exact set of pixels whose 356 pixels centers lie within the region; if \a true, return all pixels 357 that overlap with the region. */ 358 void query_strip (double theta1, double theta2, bool inclusive, 359 rangeset<I> &pixset) const; 360 /*! Returns a range set of pixels whose centers lie within the colatitude 361 range defined by \a theta1 and \a theta2 (if \a inclusive==false), or 362 which overlap with this region (if \a inclusive==true). If 363 \a theta1<theta2, the region between both angles is considered, 364 otherwise the regions \a 0<theta<theta2 and \a theta1<theta<pi. 365 \param theta1 first colatitude 366 \param theta2 second colatitude 367 \param inclusive if \a false, return the exact set of pixels whose 368 pixels centers lie within the region; if \a true, return all pixels 369 that overlap with the region. */ 370 rangeset<I> query_strip (double theta1, double theta2, bool inclusive) const 371 { 372 rangeset<I> res; 373 query_strip(theta1, theta2, inclusive, res); 374 return res; 375 } 376 377 /*! Returns useful information about a given ring of the map. 378 \param ring the ring number (the number of the first ring is 1) 379 \param startpix the number of the first pixel in the ring 380 (NOTE: this is always given in the RING numbering scheme!) 381 \param ringpix the number of pixels in the ring 382 \param costheta the cosine of the colatitude of the ring 383 \param sintheta the sine of the colatitude of the ring 384 \param shifted if \a true, the center of the first pixel is not at 385 \a phi=0 */ 386 void get_ring_info (I ring, I &startpix, I &ringpix, 387 double &costheta, double &sintheta, bool &shifted) const; 388 /*! Returns useful information about a given ring of the map. 389 \param ring the ring number (the number of the first ring is 1) 390 \param startpix the number of the first pixel in the ring 391 (NOTE: this is always given in the RING numbering scheme!) 392 \param ringpix the number of pixels in the ring 393 \param theta the colatitude (in radians) of the ring 394 \param shifted if \a true, the center of the first pixel is not at 395 \a phi=0 */ 396 void get_ring_info2 (I ring, I &startpix, I &ringpix, 397 double &theta, bool &shifted) const; 398 /*! Returns useful information about a given ring of the map. 399 \param ring the ring number (the number of the first ring is 1) 400 \param startpix the number of the first pixel in the ring 401 (NOTE: this is always given in the RING numbering scheme!) 402 \param ringpix the number of pixels in the ring 403 \param shifted if \a true, the center of the first pixel is not at 404 \a phi=0 */ 405 void get_ring_info_small (I ring, I &startpix, I &ringpix, 406 bool &shifted) const; 407 /*! Returns the neighboring pixels of \a pix in \a result. 408 On exit, \a result contains (in this order) 409 the pixel numbers of the SW, W, NW, N, NE, E, SE and S neighbor 410 of \a pix. If a neighbor does not exist (this can only be the case 411 for the W, N, E and S neighbors), its entry is set to -1. 412 413 \note This method works in both RING and NEST schemes, but is 414 considerably faster in the NEST scheme. */ 415 void neighbors (I pix, fix_arr<I,8> &result) const; 416 /*! Returns interpolation information for the direction \a ptg. 417 The surrounding pixels are returned in \a pix, their corresponding 418 weights in \a wgt. 419 \note This method works in both RING and NEST schemes, but is 420 considerably faster in the RING scheme. */ 421 void get_interpol (const pointing &ptg, fix_arr<I,4> &pix, 422 fix_arr<double,4> &wgt) const; 423 424 /*! Returns the order parameter of the object. */ 425 int Order() const { return order_; } 426 /*! Returns the \a N_side parameter of the object. */ 427 I Nside() const { return nside_; } 428 /*! Returns the number of pixels of the object. */ 429 I Npix() const { return npix_; } 430 /*! Returns the ordering scheme of the object. */ 431 Healpix_Ordering_Scheme Scheme() const { return scheme_; } 432 433 /*! Returns \a true, if both objects have the same nside and scheme, 434 else \a false. */ 435 bool conformable (const T_Healpix_Base &other) const 436 { return ((nside_==other.nside_) && (scheme_==other.scheme_)); } 437 438 /*! Swaps the contents of two Healpix_Base objects. */ 439 void swap (T_Healpix_Base &other); 440 441 /*! Returns the maximum angular distance (in radian) between any pixel 442 center and its corners. */ 443 double max_pixrad() const; 444 445 /*! Returns the maximum angular distance (in radian) between any pixel 446 center and its corners in a given ring. */ 447 double max_pixrad(I ring) const; 448 449 /*! Returns a set of points along the boundary of the given pixel. 450 \a step=1 gives 4 points on the corners. The first point corresponds 451 to the northernmost corner, the subsequent points follow the pixel 452 boundary through west, south and east corners. 453 \param pix pixel index number 454 \param step the number of returned points is 4*step. */ 455 void boundaries (I pix, tsize step, std::vector<vec3> &out) const; 456 457 arr<int> swap_cycles() const; 458 }; 459 460 /*! T_Healpix_Base for Nside up to 2^13. */ 461 typedef T_Healpix_Base<int> Healpix_Base; 462 /*! T_Healpix_Base for Nside up to 2^29. */ 463 typedef T_Healpix_Base<int64> Healpix_Base2; 464 465 #endif 466