boost/geometry/srs/projections/proj/healpix.hpp

   1 // Boost.Geometry - gis-projections (based on PROJ4)
   2
   3 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
   4
   5 // This file was modified by Oracle on 2017, 2018, 2019.
   6 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
   7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
   8
   9 // Use, modification and distribution is subject to the Boost Software License,
  10 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  11 // http://www.boost.org/LICENSE_1_0.txt)
  12
  13 // This file is converted from PROJ4, http://trac.osgeo.org/proj
  14 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  15 // PROJ4 is maintained by Frank Warmerdam
  16 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
  17
  18 // Last updated version of proj: 5.0.0
  19
  20 // Original copyright notice:
  21
  22 // Purpose: Implementation of the HEALPix and rHEALPix projections.
  23 //          For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>.
  24 // Authors: Alex Raichev (raichev@cs.auckland.ac.nz)
  25 //          Michael Speth (spethm@landcareresearch.co.nz)
  26 // Notes:   Raichev implemented these projections in Python and
  27 //          Speth translated them into C here.
  28
  29 // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
  30
  31 // Permission is hereby granted, free of charge, to any person obtaining a
  32 // copy of this software and associated documentation files (the "Software"),
  33 // to deal in the Software without restriction, including without limitation
  34 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  35 // and/or sell copies of the Software, and to permit persons to whom the
  36 // Software is furnished to do so, subject to the following conditions:
  37
  38 // The above copyright notice and this permission notice shall be included
  39 // in all copies or substantial portions of the Software.
  40
  41 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  42 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  43 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  44 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  45 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  46 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  47 // DEALINGS IN THE SOFTWARE.
  48
  49 #ifndef BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP
  50 #define BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP
  51
  52 #include <boost/geometry/srs/projections/impl/base_static.hpp>
  53 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
  54 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
  55 #include <boost/geometry/srs/projections/impl/pj_auth.hpp>
  56 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
  57 #include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>
  58 #include <boost/geometry/srs/projections/impl/projects.hpp>
  59
  60 #include <boost/geometry/util/math.hpp>
  61
  62 namespace boost { namespace geometry
  63 {
  64
  65 namespace projections
  66 {
  67     #ifndef DOXYGEN_NO_DETAIL
  68     namespace detail { namespace healpix
  69     {
  70
  71             /* Fuzz to handle rounding errors: */
  72             static const double epsilon = 1e-15;
  73
  74             template <typename T>
  75             struct par_healpix
  76             {
  77                 T qp;
  78                 detail::apa<T> apa;
  79                 int north_square;
  80                 int south_square;
  81             };
  82
  83             template <typename T>
  84             struct cap_map
  85             {
  86                 T x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */
  87                 int cn; /* An integer 0--3 indicating the position of the polar cap. */
  88                 enum region_type {north, south, equatorial} region;
  89             };
  90             template <typename T>
  91             struct point_xy
  92             {
  93                 T x, y;
  94             };
  95
  96             /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/
  97             static double rot[7][2][2] = {
  98                 /* Identity matrix */
  99                 {{1, 0},{0, 1}},
 100                 /* Matrix for counterclockwise rotation by pi/2: */
 101                 {{ 0,-1},{ 1, 0}},
 102                 /* Matrix for counterclockwise rotation by pi: */
 103                 {{-1, 0},{ 0,-1}},
 104                 /* Matrix for counterclockwise rotation by 3*pi/2:  */
 105                 {{ 0, 1},{-1, 0}},
 106                 {{ 0, 1},{-1, 0}}, // 3*pi/2
 107                 {{-1, 0},{ 0,-1}}, // pi
 108                 {{ 0,-1},{ 1, 0}}  // pi/2
 109             };
 110
 111             /**
 112              * Returns the sign of the double.
 113              * @param v the parameter whose sign is returned.
 114              * @return 1 for positive number, -1 for negative, and 0 for zero.
 115              **/
 116             template <typename T>
 117             inline T pj_sign (T const& v)
 118             {
 119                 return v > 0 ? 1 : (v < 0 ? -1 : 0);
 120             }
 121             /**
 122              * Return the index of the matrix in {{{1, 0},{0, 1}}, {{ 0,-1},{ 1, 0}}, {{-1, 0},{ 0,-1}}, {{ 0, 1},{-1, 0}}, {{ 0, 1},{-1, 0}}, {{-1, 0},{ 0,-1}}, {{ 0,-1},{ 1, 0}}}.
 123              * @param index ranges from -3 to 3.
 124              */
 125             inline int get_rotate_index(int index)
 126             {
 127                 switch(index) {
 128                 case 0:
 129                     return 0;
 130                 case 1:
 131                     return 1;
 132                 case 2:
 133                     return 2;
 134                 case 3:
 135                     return 3;
 136                 case -1:
 137                     return 4;
 138                 case -2:
 139                     return 5;
 140                 case -3:
 141                     return 6;
 142                 }
 143                 return 0;
 144             }
 145             /**
 146              * Return 1 if point (testx, testy) lies in the interior of the polygon
 147              * determined by the vertices in vert, and return 0 otherwise.
 148              * See http://paulbourke.net/geometry/polygonmesh/ for more details.
 149              * @param nvert the number of vertices in the polygon.
 150              * @param vert the (x, y)-coordinates of the polygon's vertices
 151              **/
 152             template <typename T>
 153             inline int pnpoly(int nvert, T vert[][2], T const& testx, T const& testy)
 154             {
 155                 int i;
 156                 int counter = 0;
 157                 T xinters;
 158                 point_xy<T> p1, p2;
 159
 160                 /* Check for boundrary cases */
 161                 for (i = 0; i < nvert; i++) {
 162                     if (testx == vert[i][0] && testy == vert[i][1]) {
 163                         return 1;
 164                     }
 165                 }
 166
 167                 p1.x = vert[0][0];
 168                 p1.y = vert[0][1];
 169
 170                 for (i = 1; i < nvert; i++) {
 171                     p2.x = vert[i % nvert][0];
 172                     p2.y = vert[i % nvert][1];
 173                     if (testy > (std::min)(p1.y, p2.y)  &&
 174                         testy <= (std::max)(p1.y, p2.y) &&
 175                         testx <= (std::max)(p1.x, p2.x) &&
 176                         p1.y != p2.y)
 177                     {
 178                         xinters = (testy-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x;
 179                         if (p1.x == p2.x || testx <= xinters)
 180                             counter++;
 181                     }
 182                     p1 = p2;
 183                 }
 184
 185                 if (counter % 2 == 0) {
 186                     return 0;
 187                 } else {
 188                     return 1;
 189                 }
 190             }
 191             /**
 192              * Return 1 if (x, y) lies in (the interior or boundary of) the image of the
 193              * HEALPix projection (in case proj=0) or in the image the rHEALPix projection
 194              * (in case proj=1), and return 0 otherwise.
 195              * @param north_square the position of the north polar square (rHEALPix only)
 196              * @param south_square the position of the south polar square (rHEALPix only)
 197              **/
 198             template <typename T>
 199             inline int in_image(T const& x, T const& y, int proj, int north_square, int south_square)
 200             {
 201                 static const T pi = detail::pi<T>();
 202                 static const T half_pi = detail::half_pi<T>();
 203                 static const T fourth_pi = detail::fourth_pi<T>();
 204
 205                 if (proj == 0) {
 206                     T healpixVertsJit[][2] = {
 207                         {-pi - epsilon,   fourth_pi},
 208                         {-3.0*fourth_pi,  half_pi + epsilon},
 209                         {-half_pi,        fourth_pi + epsilon},
 210                         {-fourth_pi,      half_pi + epsilon},
 211                         {0.0,             fourth_pi + epsilon},
 212                         {fourth_pi,       half_pi + epsilon},
 213                         {half_pi,         fourth_pi + epsilon},
 214                         {3.0*fourth_pi,   half_pi + epsilon},
 215                         {pi + epsilon,    fourth_pi},
 216                         {pi + epsilon,   -fourth_pi},
 217                         {3.0*fourth_pi,  -half_pi - epsilon},
 218                         {half_pi,        -fourth_pi - epsilon},
 219                         {fourth_pi,      -half_pi - epsilon},
 220                         {0.0,            -fourth_pi - epsilon},
 221                         {-fourth_pi,     -half_pi - epsilon},
 222                         {-half_pi,       -fourth_pi - epsilon},
 223                         {-3.0*fourth_pi, -half_pi - epsilon},
 224                         {-pi - epsilon,  -fourth_pi}
 225                     };
 226                     return pnpoly((int)sizeof(healpixVertsJit)/
 227                                   sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
 228                 } else {
 229                     T rhealpixVertsJit[][2] = {
 230                         {-pi - epsilon,                                 fourth_pi + epsilon},
 231                         {-pi + north_square*half_pi - epsilon,          fourth_pi + epsilon},
 232                         {-pi + north_square*half_pi - epsilon,          3.0*fourth_pi + epsilon},
 233                         {-pi + (north_square + 1.0)*half_pi + epsilon,  3.0*fourth_pi + epsilon},
 234                         {-pi + (north_square + 1.0)*half_pi + epsilon,  fourth_pi + epsilon},
 235                         {pi + epsilon,                                  fourth_pi + epsilon},
 236                         {pi + epsilon,                                 -fourth_pi - epsilon},
 237                         {-pi + (south_square + 1.0)*half_pi + epsilon, -fourth_pi - epsilon},
 238                         {-pi + (south_square + 1.0)*half_pi + epsilon, -3.0*fourth_pi - epsilon},
 239                         {-pi + south_square*half_pi - epsilon,         -3.0*fourth_pi - epsilon},
 240                         {-pi + south_square*half_pi - epsilon,         -fourth_pi - epsilon},
 241                         {-pi - epsilon,                                -fourth_pi - epsilon}
 242                     };
 243
 244                     return pnpoly((int)sizeof(rhealpixVertsJit)/
 245                                   sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y);
 246                 }
 247             }
 248             /**
 249              * Return the authalic latitude of latitude alpha (if inverse=0) or
 250              * return the approximate latitude of authalic latitude alpha (if inverse=1).
 251              * P contains the relavent ellipsoid parameters.
 252              **/
 253             template <typename Parameters, typename T>
 254             inline T auth_lat(const Parameters& par, const par_healpix<T>& proj_parm, T const& alpha, int inverse)
 255             {
 256                 if (inverse == 0) {
 257                     /* Authalic latitude. */
 258                     T q = pj_qsfn(sin(alpha), par.e, 1.0 - par.es);
 259                     T qp = proj_parm.qp;
 260                     T ratio = q/qp;
 261
 262                     if (math::abs(ratio) > 1) {
 263                         /* Rounding error. */
 264                         ratio = pj_sign(ratio);
 265                     }
 266
 267                     return asin(ratio);
 268                 } else {
 269                     /* Approximation to inverse authalic latitude. */
 270                     return pj_authlat(alpha, proj_parm.apa);
 271                 }
 272             }
 273             /**
 274              * Return the HEALPix projection of the longitude-latitude point lp on
 275              * the unit sphere.
 276             **/
 277             template <typename T>
 278             inline void healpix_sphere(T const& lp_lam, T const& lp_phi, T& xy_x, T& xy_y)
 279             {
 280                 static const T pi = detail::pi<T>();
 281                 static const T half_pi = detail::half_pi<T>();
 282                 static const T fourth_pi = detail::fourth_pi<T>();
 283
 284                 T lam = lp_lam;
 285                 T phi = lp_phi;
 286                 T phi0 = asin(T(2.0)/T(3.0));
 287
 288                 /* equatorial region */
 289                 if ( fabsl(phi) <= phi0) {
 290                     xy_x = lam;
 291                     xy_y = 3.0*pi/8.0*sin(phi);
 292                 } else {
 293                     T lamc;
 294                     T sigma = sqrt(3.0*(1 - math::abs(sin(phi))));
 295                     T cn = floor(2*lam / pi + 2);
 296                     if (cn >= 4) {
 297                         cn = 3;
 298                     }
 299                     lamc = -3*fourth_pi + half_pi*cn;
 300                     xy_x = lamc + (lam - lamc)*sigma;
 301                     xy_y = pj_sign(phi)*fourth_pi*(2 - sigma);
 302                 }
 303                 return;
 304             }
 305             /**
 306              * Return the inverse of healpix_sphere().
 307             **/
 308             template <typename T>
 309             inline void healpix_sphere_inverse(T const& xy_x, T const& xy_y, T& lp_lam, T& lp_phi)
 310             {
 311                 static const T pi = detail::pi<T>();
 312                 static const T half_pi = detail::half_pi<T>();
 313                 static const T fourth_pi = detail::fourth_pi<T>();
 314
 315                 T x = xy_x;
 316                 T y = xy_y;
 317                 T y0 = fourth_pi;
 318
 319                 /* Equatorial region. */
 320                 if (math::abs(y) <= y0) {
 321                     lp_lam = x;
 322                     lp_phi = asin(8.0*y/(3.0*pi));
 323                 } else if (fabsl(y) < half_pi) {
 324                     T cn = floor(2.0*x/pi + 2.0);
 325                     T xc, tau;
 326                     if (cn >= 4) {
 327                         cn = 3;
 328                     }
 329                     xc = -3.0*fourth_pi + (half_pi)*cn;
 330                     tau = 2.0 - 4.0*fabsl(y)/pi;
 331                     lp_lam = xc + (x - xc)/tau;
 332                     lp_phi = pj_sign(y)*asin(1.0 - math::pow(tau, 2)/3.0);
 333                 } else {
 334                     lp_lam = -1.0*pi;
 335                     lp_phi = pj_sign(y)*half_pi;
 336                 }
 337                 return;
 338             }
 339             /**
 340              * Return the vector sum a + b, where a and b are 2-dimensional vectors.
 341              * @param ret holds a + b.
 342              **/
 343             template <typename T>
 344             inline void vector_add(const T a[2], const T b[2], T ret[2])
 345             {
 346                 int i;
 347                 for(i = 0; i < 2; i++) {
 348                     ret[i] = a[i] + b[i];
 349                 }
 350             }
 351             /**
 352              * Return the vector difference a - b, where a and b are 2-dimensional vectors.
 353              * @param ret holds a - b.
 354              **/
 355             template <typename T>
 356             inline void vector_sub(const T a[2], const T b[2], T ret[2])
 357             {
 358                 int i;
 359                 for(i = 0; i < 2; i++) {
 360                     ret[i] = a[i] - b[i];
 361                 }
 362             }
 363             /**
 364              * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and
 365              * b is a 2 x 1 matrix.
 366              * @param ret holds a*b.
 367              **/
 368             template <typename T1, typename T2>
 369             inline void dot_product(const T1 a[2][2], const T2 b[2], T2 ret[2])
 370             {
 371                 int i, j;
 372                 int length = 2;
 373                 for(i = 0; i < length; i++) {
 374                     ret[i] = 0;
 375                     for(j = 0; j < length; j++) {
 376                         ret[i] += a[i][j]*b[j];
 377                     }
 378                 }
 379             }
 380             /**
 381              * Return the number of the polar cap, the pole point coordinates, and
 382              * the region that (x, y) lies in.
 383              * If inverse=0, then assume (x,y) lies in the image of the HEALPix
 384              * projection of the unit sphere.
 385              * If inverse=1, then assume (x,y) lies in the image of the
 386              * (north_square, south_square)-rHEALPix projection of the unit sphere.
 387              **/
 388             template <typename T>
 389             inline cap_map<T> get_cap(T x, T const& y, int north_square, int south_square,
 390                                      int inverse)
 391             {
 392                 static const T pi = detail::pi<T>();
 393                 static const T half_pi = detail::half_pi<T>();
 394                 static const T fourth_pi = detail::fourth_pi<T>();
 395
 396                 cap_map<T> capmap;
 397                 T c;
 398                 capmap.x = x;
 399                 capmap.y = y;
 400                 if (inverse == 0) {
 401                     if (y > fourth_pi) {
 402                         capmap.region = cap_map<T>::north;
 403                         c = half_pi;
 404                     } else if (y < -fourth_pi) {
 405                         capmap.region = cap_map<T>::south;
 406                         c = -half_pi;
 407                     } else {
 408                         capmap.region = cap_map<T>::equatorial;
 409                         capmap.cn = 0;
 410                         return capmap;
 411                     }
 412                     /* polar region */
 413                     if (x < -half_pi) {
 414                         capmap.cn = 0;
 415                         capmap.x = (-3.0*fourth_pi);
 416                         capmap.y = c;
 417                     } else if (x >= -half_pi && x < 0) {
 418                         capmap.cn = 1;
 419                         capmap.x = -fourth_pi;
 420                         capmap.y = c;
 421                     } else if (x >= 0 && x < half_pi) {
 422                         capmap.cn = 2;
 423                         capmap.x = fourth_pi;
 424                         capmap.y = c;
 425                     } else {
 426                         capmap.cn = 3;
 427                         capmap.x = 3.0*fourth_pi;
 428                         capmap.y = c;
 429                     }
 430                 } else {
 431                     if (y > fourth_pi) {
 432                         capmap.region = cap_map<T>::north;
 433                         capmap.x = (-3.0*fourth_pi + north_square*half_pi);
 434                         capmap.y = half_pi;
 435                         x = x - north_square*half_pi;
 436                     } else if (y < -fourth_pi) {
 437                         capmap.region = cap_map<T>::south;
 438                         capmap.x = (-3.0*fourth_pi + south_square*pi/2);
 439                         capmap.y = -half_pi;
 440                         x = x - south_square*half_pi;
 441                     } else {
 442                         capmap.region = cap_map<T>::equatorial;
 443                         capmap.cn = 0;
 444                         return capmap;
 445                     }
 446                     /* Polar Region, find the HEALPix polar cap number that
 447                        x, y moves to when rHEALPix polar square is disassembled. */
 448                     if (capmap.region == cap_map<T>::north) {
 449                         if (y >= -x - fourth_pi - epsilon && y < x + 5.0*fourth_pi - epsilon) {
 450                             capmap.cn = (north_square + 1) % 4;
 451                         } else if (y > -x -fourth_pi + epsilon && y >= x + 5.0*fourth_pi - epsilon) {
 452                             capmap.cn = (north_square + 2) % 4;
 453                         } else if (y <= -x -fourth_pi + epsilon && y > x + 5.0*fourth_pi + epsilon) {
 454                             capmap.cn = (north_square + 3) % 4;
 455                         } else {
 456                             capmap.cn = north_square;
 457                         }
 458                     } else if (capmap.region == cap_map<T>::south) {
 459                         if (y <= x + fourth_pi + epsilon && y > -x - 5.0*fourth_pi + epsilon) {
 460                             capmap.cn = (south_square + 1) % 4;
 461                         } else if (y < x + fourth_pi - epsilon && y <= -x - 5.0*fourth_pi + epsilon) {
 462                             capmap.cn = (south_square + 2) % 4;
 463                         } else if (y >= x + fourth_pi - epsilon && y < -x - 5.0*fourth_pi - epsilon) {
 464                             capmap.cn = (south_square + 3) % 4;
 465                         } else {
 466                             capmap.cn = south_square;
 467                         }
 468                     }
 469                 }
 470                 return capmap;
 471             }
 472             /**
 473              * Rearrange point (x, y) in the HEALPix projection by
 474              * combining the polar caps into two polar squares.
 475              * Put the north polar square in position north_square and
 476              * the south polar square in position south_square.
 477              * If inverse=1, then uncombine the polar caps.
 478              * @param north_square integer between 0 and 3.
 479              * @param south_square integer between 0 and 3.
 480              **/
 481             template <typename T>
 482             inline void combine_caps(T& xy_x, T& xy_y, int north_square, int south_square,
 483                                      int inverse)
 484             {
 485                 static const T half_pi = detail::half_pi<T>();
 486                 static const T fourth_pi = detail::fourth_pi<T>();
 487
 488                 T v[2];
 489                 T c[2];
 490                 T vector[2];
 491                 T v_min_c[2];
 492                 T ret_dot[2];
 493                 const double (*tmpRot)[2];
 494                 int pole = 0;
 495
 496                 cap_map<T> capmap = get_cap(xy_x, xy_y, north_square, south_square, inverse);
 497                 if (capmap.region == cap_map<T>::equatorial) {
 498                     xy_x = capmap.x;
 499                     xy_y = capmap.y;
 500                     return;
 501                 }
 502
 503                 v[0] = xy_x; v[1] = xy_y;
 504                 c[0] = capmap.x; c[1] = capmap.y;
 505
 506                 if (inverse == 0) {
 507                     /* Rotate (xy_x, xy_y) about its polar cap tip and then translate it to
 508                        north_square or south_square. */
 509
 510                     if (capmap.region == cap_map<T>::north) {
 511                         pole = north_square;
 512                         tmpRot = rot[get_rotate_index(capmap.cn - pole)];
 513                     } else {
 514                         pole = south_square;
 515                         tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
 516                     }
 517                 } else {
 518                     /* Inverse function.
 519                      Unrotate (xy_x, xy_y) and then translate it back. */
 520
 521                     /* disassemble */
 522                     if (capmap.region == cap_map<T>::north) {
 523                         pole = north_square;
 524                         tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
 525                     } else {
 526                         pole = south_square;
 527                         tmpRot = rot[get_rotate_index(capmap.cn - pole)];
 528                     }
 529                 }
 530
 531                 vector_sub(v, c, v_min_c);
 532                 dot_product(tmpRot, v_min_c, ret_dot);
 533
 534                 {
 535                     T a[2];
 536                     /* Workaround cppcheck git issue */
 537                     T* pa = a;
 538                     // TODO: in proj4 5.0.0 this line is used instead
 539                     //pa[0] = -3.0*fourth_pi + ((inverse == 0) ? 0 : capmap.cn) *half_pi;
 540                     pa[0] = -3.0*fourth_pi + ((inverse == 0) ? pole : capmap.cn) *half_pi;
 541                     pa[1] = half_pi;
 542                     vector_add(ret_dot, a, vector);
 543                 }
 544
 545                 xy_x = vector[0];
 546                 xy_y = vector[1];
 547             }
 548
 549             template <typename T, typename Parameters>
 550             struct base_healpix_ellipsoid
 551             {
 552                 par_healpix<T> m_proj_parm;
 553
 554                 // FORWARD(e_healpix_forward)  ellipsoid
 555                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 556                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
 557                 {
 558                     lp_lat = auth_lat(par, m_proj_parm, lp_lat, 0);
 559                     return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y);
 560                 }
 561
 562                 // INVERSE(e_healpix_inverse)  ellipsoid
 563                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 564                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
 565                 {
 566                     /* Check whether (x, y) lies in the HEALPix image. */
 567                     if (in_image(xy_x, xy_y, 0, 0, 0) == 0) {
 568                         lp_lon = HUGE_VAL;
 569                         lp_lat = HUGE_VAL;
 570                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) );
 571                     }
 572                     healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat);
 573                     lp_lat = auth_lat(par, m_proj_parm, lp_lat, 1);
 574                 }
 575
 576                 static inline std::string get_name()
 577                 {
 578                     return "healpix_ellipsoid";
 579                 }
 580
 581             };
 582
 583             template <typename T, typename Parameters>
 584             struct base_healpix_spheroid
 585             {
 586                 par_healpix<T> m_proj_parm;
 587
 588                 // FORWARD(s_healpix_forward)  sphere
 589                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 590                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 591                 {
 592                     return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y);
 593                 }
 594
 595                 // INVERSE(s_healpix_inverse)  sphere
 596                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 597                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
 598                 {
 599                     /* Check whether (x, y) lies in the HEALPix image */
 600                     if (in_image(xy_x, xy_y, 0, 0, 0) == 0) {
 601                         lp_lon = HUGE_VAL;
 602                         lp_lat = HUGE_VAL;
 603                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) );
 604                     }
 605                     return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat);
 606                 }
 607
 608                 static inline std::string get_name()
 609                 {
 610                     return "healpix_spheroid";
 611                 }
 612
 613             };
 614
 615             template <typename T, typename Parameters>
 616             struct base_rhealpix_ellipsoid
 617             {
 618                 par_healpix<T> m_proj_parm;
 619
 620                 // FORWARD(e_rhealpix_forward)  ellipsoid
 621                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 622                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
 623                 {
 624                     lp_lat = auth_lat(par, m_proj_parm, lp_lat, 0);
 625                     healpix_sphere(lp_lon, lp_lat, xy_x, xy_y);
 626                     combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0);
 627                 }
 628
 629                 // INVERSE(e_rhealpix_inverse)  ellipsoid
 630                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 631                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 632                 {
 633                     /* Check whether (x, y) lies in the rHEALPix image. */
 634                     if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) {
 635                         lp_lon = HUGE_VAL;
 636                         lp_lat = HUGE_VAL;
 637                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) );
 638                     }
 639                     combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1);
 640                     healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat);
 641                     lp_lat = auth_lat(par, m_proj_parm, lp_lat, 1);
 642                 }
 643
 644                 static inline std::string get_name()
 645                 {
 646                     return "rhealpix_ellipsoid";
 647                 }
 648
 649             };
 650
 651             template <typename T, typename Parameters>
 652             struct base_rhealpix_spheroid
 653             {
 654                 par_healpix<T> m_proj_parm;
 655
 656                 // FORWARD(s_rhealpix_forward)  sphere
 657                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
 658                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
 659                 {
 660                     healpix_sphere(lp_lon, lp_lat, xy_x, xy_y);
 661                     combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0);
 662                 }
 663
 664                 // INVERSE(s_rhealpix_inverse)  sphere
 665                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
 666                 inline void inv(Parameters const& , T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
 667                 {
 668                     /* Check whether (x, y) lies in the rHEALPix image. */
 669                     if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) {
 670                         lp_lon = HUGE_VAL;
 671                         lp_lat = HUGE_VAL;
 672                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) );
 673                     }
 674                     combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1);
 675                     return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat);
 676                 }
 677
 678                 static inline std::string get_name()
 679                 {
 680                     return "rhealpix_spheroid";
 681                 }
 682
 683             };
 684
 685             // HEALPix
 686             template <typename Parameters, typename T>
 687             inline void setup_healpix(Parameters& par, par_healpix<T>& proj_parm)
 688             {
 689                 if (par.es != 0.0) {
 690                     proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */
 691                     proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */
 692                     par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */
 693                     pj_calc_ellipsoid_params(par, par.a, par.es); /* Ensure we have a consistent parameter set */
 694                 } else {
 695                 }
 696             }
 697
 698             // rHEALPix
 699             template <typename Params, typename Parameters, typename T>
 700             inline void setup_rhealpix(Params const& params, Parameters& par, par_healpix<T>& proj_parm)
 701             {
 702                 proj_parm.north_square = pj_get_param_i<srs::spar::north_square>(params, "north_square", srs::dpar::north_square);
 703                 proj_parm.south_square = pj_get_param_i<srs::spar::south_square>(params, "south_square", srs::dpar::south_square);
 704                 /* Check for valid north_square and south_square inputs. */
 705                 if ((proj_parm.north_square < 0) || (proj_parm.north_square > 3)) {
 706                     BOOST_THROW_EXCEPTION( projection_exception(error_axis) );
 707                 }
 708                 if ((proj_parm.south_square < 0) || (proj_parm.south_square > 3)) {
 709                     BOOST_THROW_EXCEPTION( projection_exception(error_axis) );
 710                 }
 711                 if (par.es != 0.0) {
 712                     proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */
 713                     proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */
 714                     par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */
 715                     // TODO: why not the same as in healpix?
 716                     //pj_calc_ellipsoid_params(par, par.a, par.es);
 717                     par.ra = 1.0/par.a;
 718                 } else {
 719                 }
 720             }
 721
 722     }} // namespace detail::healpix
 723     #endif // doxygen
 724
 725     /*!
 726         \brief HEALPix projection
 727         \ingroup projections
 728         \tparam Geographic latlong point type
 729         \tparam Cartesian xy point type
 730         \tparam Parameters parameter type
 731         \par Projection characteristics
 732          - Spheroid
 733          - Ellipsoid
 734         \par Example
 735         \image html ex_healpix.gif
 736     */
 737     template <typename T, typename Parameters>
 738     struct healpix_ellipsoid : public detail::healpix::base_healpix_ellipsoid<T, Parameters>
 739     {
 740         template <typename Params>
 741         inline healpix_ellipsoid(Params const& , Parameters & par)
 742         {
 743             detail::healpix::setup_healpix(par, this->m_proj_parm);
 744         }
 745     };
 746
 747     /*!
 748         \brief HEALPix projection
 749         \ingroup projections
 750         \tparam Geographic latlong point type
 751         \tparam Cartesian xy point type
 752         \tparam Parameters parameter type
 753         \par Projection characteristics
 754          - Spheroid
 755          - Ellipsoid
 756         \par Example
 757         \image html ex_healpix.gif
 758     */
 759     template <typename T, typename Parameters>
 760     struct healpix_spheroid : public detail::healpix::base_healpix_spheroid<T, Parameters>
 761     {
 762         template <typename Params>
 763         inline healpix_spheroid(Params const& , Parameters & par)
 764         {
 765             detail::healpix::setup_healpix(par, this->m_proj_parm);
 766         }
 767     };
 768
 769     /*!
 770         \brief rHEALPix projection
 771         \ingroup projections
 772         \tparam Geographic latlong point type
 773         \tparam Cartesian xy point type
 774         \tparam Parameters parameter type
 775         \par Projection characteristics
 776          - Spheroid
 777          - Ellipsoid
 778         \par Projection parameters
 779          - north_square (integer)
 780          - south_square (integer)
 781         \par Example
 782         \image html ex_rhealpix.gif
 783     */
 784     template <typename T, typename Parameters>
 785     struct rhealpix_ellipsoid : public detail::healpix::base_rhealpix_ellipsoid<T, Parameters>
 786     {
 787         template <typename Params>
 788         inline rhealpix_ellipsoid(Params const& params, Parameters & par)
 789         {
 790             detail::healpix::setup_rhealpix(params, par, this->m_proj_parm);
 791         }
 792     };
 793
 794     /*!
 795         \brief rHEALPix projection
 796         \ingroup projections
 797         \tparam Geographic latlong point type
 798         \tparam Cartesian xy point type
 799         \tparam Parameters parameter type
 800         \par Projection characteristics
 801          - Spheroid
 802          - Ellipsoid
 803         \par Projection parameters
 804          - north_square (integer)
 805          - south_square (integer)
 806         \par Example
 807         \image html ex_rhealpix.gif
 808     */
 809     template <typename T, typename Parameters>
 810     struct rhealpix_spheroid : public detail::healpix::base_rhealpix_spheroid<T, Parameters>
 811     {
 812         template <typename Params>
 813         inline rhealpix_spheroid(Params const& params, Parameters & par)
 814         {
 815             detail::healpix::setup_rhealpix(params, par, this->m_proj_parm);
 816         }
 817     };
 818
 819     #ifndef DOXYGEN_NO_DETAIL
 820     namespace detail
 821     {
 822
 823         // Static projection
 824         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_healpix, healpix_spheroid, healpix_ellipsoid)
 825         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_rhealpix, rhealpix_spheroid, rhealpix_ellipsoid)
 826
 827         // Factory entry(s)
 828         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(healpix_entry, healpix_spheroid, healpix_ellipsoid)
 829         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(rhealpix_entry, rhealpix_spheroid, rhealpix_ellipsoid)
 830
 831         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(healpix_init)
 832         {
 833             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(healpix, healpix_entry)
 834             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(rhealpix, rhealpix_entry)
 835         }
 836
 837     } // namespace detail
 838     #endif // doxygen
 839
 840 } // namespace projections
 841
 842 }} // namespace boost::geometry
 843
 844 #endif // BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP
 845