projections/proj/healpix.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018, 2019.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // Purpose: Implementation of the HEALPix and rHEALPix projections.
0023 //          For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>.
0024 // Authors: Alex Raichev (raichev@cs.auckland.ac.nz)
0025 //          Michael Speth (spethm@landcareresearch.co.nz)
0026 // Notes:   Raichev implemented these projections in Python and
0027 //          Speth translated them into C here.
0028
0029 // Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
0030
0031 // Permission is hereby granted, free of charge, to any person obtaining a
0032 // copy of this software and associated documentation files (the "Software"),
0033 // to deal in the Software without restriction, including without limitation
0034 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0035 // and/or sell copies of the Software, and to permit persons to whom the
0036 // Software is furnished to do so, subject to the following conditions:
0037
0038 // The above copyright notice and this permission notice shall be included
0039 // in all copies or substantial portions of the Software.
0040
0041 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0042 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0043 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0044 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0045 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0046 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0047 // DEALINGS IN THE SOFTWARE.
0048
0049 #ifndef BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP
0050 #define BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP
0051
0052 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0053 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0054 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0055 #include <boost/geometry/srs/projections/impl/pj_auth.hpp>
0056 #include <boost/geometry/srs/projections/impl/pj_param.hpp>
0057 #include <boost/geometry/srs/projections/impl/pj_qsfn.hpp>
0058 #include <boost/geometry/srs/projections/impl/projects.hpp>
0059
0060 #include <boost/geometry/util/math.hpp>
0061
0062 namespace boost { namespace geometry
0063 {
0064
0065 namespace projections
0066 {
0067     #ifndef DOXYGEN_NO_DETAIL
0068     namespace detail { namespace healpix
0069     {
0070
0071             /* Fuzz to handle rounding errors: */
0072             static const double epsilon = 1e-15;
0073
0074             template <typename T>
0075             struct par_healpix
0076             {
0077                 T qp;
0078                 detail::apa<T> apa;
0079                 int north_square;
0080                 int south_square;
0081             };
0082
0083             template <typename T>
0084             struct cap_map
0085             {
0086                 T x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */
0087                 int cn; /* An integer 0--3 indicating the position of the polar cap. */
0088                 enum region_type {north, south, equatorial} region;
0089             };
0090             template <typename T>
0091             struct point_xy
0092             {
0093                 T x, y;
0094             };
0095
0096             /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/
0097             static double rot[7][2][2] = {
0098                 /* Identity matrix */
0099                 {{1, 0},{0, 1}},
0100                 /* Matrix for counterclockwise rotation by pi/2: */
0101                 {{ 0,-1},{ 1, 0}},
0102                 /* Matrix for counterclockwise rotation by pi: */
0103                 {{-1, 0},{ 0,-1}},
0104                 /* Matrix for counterclockwise rotation by 3*pi/2:  */
0105                 {{ 0, 1},{-1, 0}},
0106                 {{ 0, 1},{-1, 0}}, // 3*pi/2
0107                 {{-1, 0},{ 0,-1}}, // pi
0108                 {{ 0,-1},{ 1, 0}}  // pi/2
0109             };
0110
0111             /**
0112              * Returns the sign of the double.
0113              * @param v the parameter whose sign is returned.
0114              * @return 1 for positive number, -1 for negative, and 0 for zero.
0115              **/
0116             template <typename T>
0117             inline T pj_sign (T const& v)
0118             {
0119                 return v > 0 ? 1 : (v < 0 ? -1 : 0);
0120             }
0121             /**
0122              * 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}}}.
0123              * @param index ranges from -3 to 3.
0124              */
0125             inline int get_rotate_index(int index)
0126             {
0127                 switch(index) {
0128                 case 0:
0129                     return 0;
0130                 case 1:
0131                     return 1;
0132                 case 2:
0133                     return 2;
0134                 case 3:
0135                     return 3;
0136                 case -1:
0137                     return 4;
0138                 case -2:
0139                     return 5;
0140                 case -3:
0141                     return 6;
0142                 }
0143                 return 0;
0144             }
0145             /**
0146              * Return 1 if point (testx, testy) lies in the interior of the polygon
0147              * determined by the vertices in vert, and return 0 otherwise.
0148              * See http://paulbourke.net/geometry/polygonmesh/ for more details.
0149              * @param nvert the number of vertices in the polygon.
0150              * @param vert the (x, y)-coordinates of the polygon's vertices
0151              **/
0152             template <typename T>
0153             inline int pnpoly(int nvert, T vert[][2], T const& testx, T const& testy)
0154             {
0155                 int i;
0156                 int counter = 0;
0157                 T xinters;
0158                 point_xy<T> p1, p2;
0159
0160                 /* Check for boundrary cases */
0161                 for (i = 0; i < nvert; i++) {
0162                     if (testx == vert[i][0] && testy == vert[i][1]) {
0163                         return 1;
0164                     }
0165                 }
0166
0167                 p1.x = vert[0][0];
0168                 p1.y = vert[0][1];
0169
0170                 for (i = 1; i < nvert; i++) {
0171                     p2.x = vert[i % nvert][0];
0172                     p2.y = vert[i % nvert][1];
0173                     if (testy > (std::min)(p1.y, p2.y)  &&
0174                         testy <= (std::max)(p1.y, p2.y) &&
0175                         testx <= (std::max)(p1.x, p2.x) &&
0176                         p1.y != p2.y)
0177                     {
0178                         xinters = (testy-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x;
0179                         if (p1.x == p2.x || testx <= xinters)
0180                             counter++;
0181                     }
0182                     p1 = p2;
0183                 }
0184
0185                 if (counter % 2 == 0) {
0186                     return 0;
0187                 } else {
0188                     return 1;
0189                 }
0190             }
0191             /**
0192              * Return 1 if (x, y) lies in (the interior or boundary of) the image of the
0193              * HEALPix projection (in case proj=0) or in the image the rHEALPix projection
0194              * (in case proj=1), and return 0 otherwise.
0195              * @param north_square the position of the north polar square (rHEALPix only)
0196              * @param south_square the position of the south polar square (rHEALPix only)
0197              **/
0198             template <typename T>
0199             inline int in_image(T const& x, T const& y, int proj, int north_square, int south_square)
0200             {
0201                 static const T pi = detail::pi<T>();
0202                 static const T half_pi = detail::half_pi<T>();
0203                 static const T fourth_pi = detail::fourth_pi<T>();
0204
0205                 if (proj == 0) {
0206                     T healpixVertsJit[][2] = {
0207                         {-pi - epsilon,   fourth_pi},
0208                         {-3.0*fourth_pi,  half_pi + epsilon},
0209                         {-half_pi,        fourth_pi + epsilon},
0210                         {-fourth_pi,      half_pi + epsilon},
0211                         {0.0,             fourth_pi + epsilon},
0212                         {fourth_pi,       half_pi + epsilon},
0213                         {half_pi,         fourth_pi + epsilon},
0214                         {3.0*fourth_pi,   half_pi + epsilon},
0215                         {pi + epsilon,    fourth_pi},
0216                         {pi + epsilon,   -fourth_pi},
0217                         {3.0*fourth_pi,  -half_pi - epsilon},
0218                         {half_pi,        -fourth_pi - epsilon},
0219                         {fourth_pi,      -half_pi - epsilon},
0220                         {0.0,            -fourth_pi - epsilon},
0221                         {-fourth_pi,     -half_pi - epsilon},
0222                         {-half_pi,       -fourth_pi - epsilon},
0223                         {-3.0*fourth_pi, -half_pi - epsilon},
0224                         {-pi - epsilon,  -fourth_pi}
0225                     };
0226                     return pnpoly((int)sizeof(healpixVertsJit)/
0227                                   sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
0228                 } else {
0229                     T rhealpixVertsJit[][2] = {
0230                         {-pi - epsilon,                                 fourth_pi + epsilon},
0231                         {-pi + north_square*half_pi - epsilon,          fourth_pi + epsilon},
0232                         {-pi + north_square*half_pi - epsilon,          3.0*fourth_pi + epsilon},
0233                         {-pi + (north_square + 1.0)*half_pi + epsilon,  3.0*fourth_pi + epsilon},
0234                         {-pi + (north_square + 1.0)*half_pi + epsilon,  fourth_pi + epsilon},
0235                         {pi + epsilon,                                  fourth_pi + epsilon},
0236                         {pi + epsilon,                                 -fourth_pi - epsilon},
0237                         {-pi + (south_square + 1.0)*half_pi + epsilon, -fourth_pi - epsilon},
0238                         {-pi + (south_square + 1.0)*half_pi + epsilon, -3.0*fourth_pi - epsilon},
0239                         {-pi + south_square*half_pi - epsilon,         -3.0*fourth_pi - epsilon},
0240                         {-pi + south_square*half_pi - epsilon,         -fourth_pi - epsilon},
0241                         {-pi - epsilon,                                -fourth_pi - epsilon}
0242                     };
0243
0244                     return pnpoly((int)sizeof(rhealpixVertsJit)/
0245                                   sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y);
0246                 }
0247             }
0248             /**
0249              * Return the authalic latitude of latitude alpha (if inverse=0) or
0250              * return the approximate latitude of authalic latitude alpha (if inverse=1).
0251              * P contains the relavent ellipsoid parameters.
0252              **/
0253             template <typename Parameters, typename T>
0254             inline T auth_lat(const Parameters& par, const par_healpix<T>& proj_parm, T const& alpha, int inverse)
0255             {
0256                 if (inverse == 0) {
0257                     /* Authalic latitude. */
0258                     T q = pj_qsfn(sin(alpha), par.e, 1.0 - par.es);
0259                     T qp = proj_parm.qp;
0260                     T ratio = q/qp;
0261
0262                     if (math::abs(ratio) > 1) {
0263                         /* Rounding error. */
0264                         ratio = pj_sign(ratio);
0265                     }
0266
0267                     return asin(ratio);
0268                 } else {
0269                     /* Approximation to inverse authalic latitude. */
0270                     return pj_authlat(alpha, proj_parm.apa);
0271                 }
0272             }
0273             /**
0274              * Return the HEALPix projection of the longitude-latitude point lp on
0275              * the unit sphere.
0276             **/
0277             template <typename T>
0278             inline void healpix_sphere(T const& lp_lam, T const& lp_phi, T& xy_x, T& xy_y)
0279             {
0280                 static const T pi = detail::pi<T>();
0281                 static const T half_pi = detail::half_pi<T>();
0282                 static const T fourth_pi = detail::fourth_pi<T>();
0283
0284                 T lam = lp_lam;
0285                 T phi = lp_phi;
0286                 T phi0 = asin(T(2.0)/T(3.0));
0287
0288                 /* equatorial region */
0289                 if ( fabsl(phi) <= phi0) {
0290                     xy_x = lam;
0291                     xy_y = 3.0*pi/8.0*sin(phi);
0292                 } else {
0293                     T lamc;
0294                     T sigma = sqrt(3.0*(1 - math::abs(sin(phi))));
0295                     T cn = floor(2*lam / pi + 2);
0296                     if (cn >= 4) {
0297                         cn = 3;
0298                     }
0299                     lamc = -3*fourth_pi + half_pi*cn;
0300                     xy_x = lamc + (lam - lamc)*sigma;
0301                     xy_y = pj_sign(phi)*fourth_pi*(2 - sigma);
0302                 }
0303                 return;
0304             }
0305             /**
0306              * Return the inverse of healpix_sphere().
0307             **/
0308             template <typename T>
0309             inline void healpix_sphere_inverse(T const& xy_x, T const& xy_y, T& lp_lam, T& lp_phi)
0310             {
0311                 static const T pi = detail::pi<T>();
0312                 static const T half_pi = detail::half_pi<T>();
0313                 static const T fourth_pi = detail::fourth_pi<T>();
0314
0315                 T x = xy_x;
0316                 T y = xy_y;
0317                 T y0 = fourth_pi;
0318
0319                 /* Equatorial region. */
0320                 if (math::abs(y) <= y0) {
0321                     lp_lam = x;
0322                     lp_phi = asin(8.0*y/(3.0*pi));
0323                 } else if (fabsl(y) < half_pi) {
0324                     T cn = floor(2.0*x/pi + 2.0);
0325                     T xc, tau;
0326                     if (cn >= 4) {
0327                         cn = 3;
0328                     }
0329                     xc = -3.0*fourth_pi + (half_pi)*cn;
0330                     tau = 2.0 - 4.0*fabsl(y)/pi;
0331                     lp_lam = xc + (x - xc)/tau;
0332                     lp_phi = pj_sign(y)*asin(1.0 - math::pow(tau, 2)/3.0);
0333                 } else {
0334                     lp_lam = -1.0*pi;
0335                     lp_phi = pj_sign(y)*half_pi;
0336                 }
0337                 return;
0338             }
0339             /**
0340              * Return the vector sum a + b, where a and b are 2-dimensional vectors.
0341              * @param ret holds a + b.
0342              **/
0343             template <typename T>
0344             inline void vector_add(const T a[2], const T b[2], T ret[2])
0345             {
0346                 int i;
0347                 for(i = 0; i < 2; i++) {
0348                     ret[i] = a[i] + b[i];
0349                 }
0350             }
0351             /**
0352              * Return the vector difference a - b, where a and b are 2-dimensional vectors.
0353              * @param ret holds a - b.
0354              **/
0355             template <typename T>
0356             inline void vector_sub(const T a[2], const T b[2], T ret[2])
0357             {
0358                 int i;
0359                 for(i = 0; i < 2; i++) {
0360                     ret[i] = a[i] - b[i];
0361                 }
0362             }
0363             /**
0364              * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and
0365              * b is a 2 x 1 matrix.
0366              * @param ret holds a*b.
0367              **/
0368             template <typename T1, typename T2>
0369             inline void dot_product(const T1 a[2][2], const T2 b[2], T2 ret[2])
0370             {
0371                 int i, j;
0372                 int length = 2;
0373                 for(i = 0; i < length; i++) {
0374                     ret[i] = 0;
0375                     for(j = 0; j < length; j++) {
0376                         ret[i] += a[i][j]*b[j];
0377                     }
0378                 }
0379             }
0380             /**
0381              * Return the number of the polar cap, the pole point coordinates, and
0382              * the region that (x, y) lies in.
0383              * If inverse=0, then assume (x,y) lies in the image of the HEALPix
0384              * projection of the unit sphere.
0385              * If inverse=1, then assume (x,y) lies in the image of the
0386              * (north_square, south_square)-rHEALPix projection of the unit sphere.
0387              **/
0388             template <typename T>
0389             inline cap_map<T> get_cap(T x, T const& y, int north_square, int south_square,
0390                                      int inverse)
0391             {
0392                 static const T pi = detail::pi<T>();
0393                 static const T half_pi = detail::half_pi<T>();
0394                 static const T fourth_pi = detail::fourth_pi<T>();
0395
0396                 cap_map<T> capmap;
0397                 T c;
0398                 capmap.x = x;
0399                 capmap.y = y;
0400                 if (inverse == 0) {
0401                     if (y > fourth_pi) {
0402                         capmap.region = cap_map<T>::north;
0403                         c = half_pi;
0404                     } else if (y < -fourth_pi) {
0405                         capmap.region = cap_map<T>::south;
0406                         c = -half_pi;
0407                     } else {
0408                         capmap.region = cap_map<T>::equatorial;
0409                         capmap.cn = 0;
0410                         return capmap;
0411                     }
0412                     /* polar region */
0413                     if (x < -half_pi) {
0414                         capmap.cn = 0;
0415                         capmap.x = (-3.0*fourth_pi);
0416                         capmap.y = c;
0417                     } else if (x >= -half_pi && x < 0) {
0418                         capmap.cn = 1;
0419                         capmap.x = -fourth_pi;
0420                         capmap.y = c;
0421                     } else if (x >= 0 && x < half_pi) {
0422                         capmap.cn = 2;
0423                         capmap.x = fourth_pi;
0424                         capmap.y = c;
0425                     } else {
0426                         capmap.cn = 3;
0427                         capmap.x = 3.0*fourth_pi;
0428                         capmap.y = c;
0429                     }
0430                 } else {
0431                     if (y > fourth_pi) {
0432                         capmap.region = cap_map<T>::north;
0433                         capmap.x = (-3.0*fourth_pi + north_square*half_pi);
0434                         capmap.y = half_pi;
0435                         x = x - north_square*half_pi;
0436                     } else if (y < -fourth_pi) {
0437                         capmap.region = cap_map<T>::south;
0438                         capmap.x = (-3.0*fourth_pi + south_square*pi/2);
0439                         capmap.y = -half_pi;
0440                         x = x - south_square*half_pi;
0441                     } else {
0442                         capmap.region = cap_map<T>::equatorial;
0443                         capmap.cn = 0;
0444                         return capmap;
0445                     }
0446                     /* Polar Region, find the HEALPix polar cap number that
0447                        x, y moves to when rHEALPix polar square is disassembled. */
0448                     if (capmap.region == cap_map<T>::north) {
0449                         if (y >= -x - fourth_pi - epsilon && y < x + 5.0*fourth_pi - epsilon) {
0450                             capmap.cn = (north_square + 1) % 4;
0451                         } else if (y > -x -fourth_pi + epsilon && y >= x + 5.0*fourth_pi - epsilon) {
0452                             capmap.cn = (north_square + 2) % 4;
0453                         } else if (y <= -x -fourth_pi + epsilon && y > x + 5.0*fourth_pi + epsilon) {
0454                             capmap.cn = (north_square + 3) % 4;
0455                         } else {
0456                             capmap.cn = north_square;
0457                         }
0458                     } else if (capmap.region == cap_map<T>::south) {
0459                         if (y <= x + fourth_pi + epsilon && y > -x - 5.0*fourth_pi + epsilon) {
0460                             capmap.cn = (south_square + 1) % 4;
0461                         } else if (y < x + fourth_pi - epsilon && y <= -x - 5.0*fourth_pi + epsilon) {
0462                             capmap.cn = (south_square + 2) % 4;
0463                         } else if (y >= x + fourth_pi - epsilon && y < -x - 5.0*fourth_pi - epsilon) {
0464                             capmap.cn = (south_square + 3) % 4;
0465                         } else {
0466                             capmap.cn = south_square;
0467                         }
0468                     }
0469                 }
0470                 return capmap;
0471             }
0472             /**
0473              * Rearrange point (x, y) in the HEALPix projection by
0474              * combining the polar caps into two polar squares.
0475              * Put the north polar square in position north_square and
0476              * the south polar square in position south_square.
0477              * If inverse=1, then uncombine the polar caps.
0478              * @param north_square integer between 0 and 3.
0479              * @param south_square integer between 0 and 3.
0480              **/
0481             template <typename T>
0482             inline void combine_caps(T& xy_x, T& xy_y, int north_square, int south_square,
0483                                      int inverse)
0484             {
0485                 static const T half_pi = detail::half_pi<T>();
0486                 static const T fourth_pi = detail::fourth_pi<T>();
0487
0488                 T v[2];
0489                 T c[2];
0490                 T vector[2];
0491                 T v_min_c[2];
0492                 T ret_dot[2];
0493                 const double (*tmpRot)[2];
0494                 int pole = 0;
0495
0496                 cap_map<T> capmap = get_cap(xy_x, xy_y, north_square, south_square, inverse);
0497                 if (capmap.region == cap_map<T>::equatorial) {
0498                     xy_x = capmap.x;
0499                     xy_y = capmap.y;
0500                     return;
0501                 }
0502
0503                 v[0] = xy_x; v[1] = xy_y;
0504                 c[0] = capmap.x; c[1] = capmap.y;
0505
0506                 if (inverse == 0) {
0507                     /* Rotate (xy_x, xy_y) about its polar cap tip and then translate it to
0508                        north_square or south_square. */
0509
0510                     if (capmap.region == cap_map<T>::north) {
0511                         pole = north_square;
0512                         tmpRot = rot[get_rotate_index(capmap.cn - pole)];
0513                     } else {
0514                         pole = south_square;
0515                         tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
0516                     }
0517                 } else {
0518                     /* Inverse function.
0519                      Unrotate (xy_x, xy_y) and then translate it back. */
0520
0521                     /* disassemble */
0522                     if (capmap.region == cap_map<T>::north) {
0523                         pole = north_square;
0524                         tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
0525                     } else {
0526                         pole = south_square;
0527                         tmpRot = rot[get_rotate_index(capmap.cn - pole)];
0528                     }
0529                 }
0530
0531                 vector_sub(v, c, v_min_c);
0532                 dot_product(tmpRot, v_min_c, ret_dot);
0533
0534                 {
0535                     T a[2];
0536                     /* Workaround cppcheck git issue */
0537                     T* pa = a;
0538                     // TODO: in proj4 5.0.0 this line is used instead
0539                     //pa[0] = -3.0*fourth_pi + ((inverse == 0) ? 0 : capmap.cn) *half_pi;
0540                     pa[0] = -3.0*fourth_pi + ((inverse == 0) ? pole : capmap.cn) *half_pi;
0541                     pa[1] = half_pi;
0542                     vector_add(ret_dot, a, vector);
0543                 }
0544
0545                 xy_x = vector[0];
0546                 xy_y = vector[1];
0547             }
0548
0549             template <typename T, typename Parameters>
0550             struct base_healpix_ellipsoid
0551             {
0552                 par_healpix<T> m_proj_parm;
0553
0554                 // FORWARD(e_healpix_forward)  ellipsoid
0555                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0556                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
0557                 {
0558                     lp_lat = auth_lat(par, m_proj_parm, lp_lat, 0);
0559                     return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y);
0560                 }
0561
0562                 // INVERSE(e_healpix_inverse)  ellipsoid
0563                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0564                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0565                 {
0566                     /* Check whether (x, y) lies in the HEALPix image. */
0567                     if (in_image(xy_x, xy_y, 0, 0, 0) == 0) {
0568                         lp_lon = HUGE_VAL;
0569                         lp_lat = HUGE_VAL;
0570                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) );
0571                     }
0572                     healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat);
0573                     lp_lat = auth_lat(par, m_proj_parm, lp_lat, 1);
0574                 }
0575
0576                 static inline std::string get_name()
0577                 {
0578                     return "healpix_ellipsoid";
0579                 }
0580
0581             };
0582
0583             template <typename T, typename Parameters>
0584             struct base_healpix_spheroid
0585             {
0586                 par_healpix<T> m_proj_parm;
0587
0588                 // FORWARD(s_healpix_forward)  sphere
0589                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0590                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0591                 {
0592                     return healpix_sphere(lp_lon, lp_lat, xy_x, xy_y);
0593                 }
0594
0595                 // INVERSE(s_healpix_inverse)  sphere
0596                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0597                 inline void inv(Parameters const& , T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0598                 {
0599                     /* Check whether (x, y) lies in the HEALPix image */
0600                     if (in_image(xy_x, xy_y, 0, 0, 0) == 0) {
0601                         lp_lon = HUGE_VAL;
0602                         lp_lat = HUGE_VAL;
0603                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) );
0604                     }
0605                     return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat);
0606                 }
0607
0608                 static inline std::string get_name()
0609                 {
0610                     return "healpix_spheroid";
0611                 }
0612
0613             };
0614
0615             template <typename T, typename Parameters>
0616             struct base_rhealpix_ellipsoid
0617             {
0618                 par_healpix<T> m_proj_parm;
0619
0620                 // FORWARD(e_rhealpix_forward)  ellipsoid
0621                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0622                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
0623                 {
0624                     lp_lat = auth_lat(par, m_proj_parm, lp_lat, 0);
0625                     healpix_sphere(lp_lon, lp_lat, xy_x, xy_y);
0626                     combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0);
0627                 }
0628
0629                 // INVERSE(e_rhealpix_inverse)  ellipsoid
0630                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0631                 inline void inv(Parameters const& par, T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0632                 {
0633                     /* Check whether (x, y) lies in the rHEALPix image. */
0634                     if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) {
0635                         lp_lon = HUGE_VAL;
0636                         lp_lat = HUGE_VAL;
0637                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) );
0638                     }
0639                     combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1);
0640                     healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat);
0641                     lp_lat = auth_lat(par, m_proj_parm, lp_lat, 1);
0642                 }
0643
0644                 static inline std::string get_name()
0645                 {
0646                     return "rhealpix_ellipsoid";
0647                 }
0648
0649             };
0650
0651             template <typename T, typename Parameters>
0652             struct base_rhealpix_spheroid
0653             {
0654                 par_healpix<T> m_proj_parm;
0655
0656                 // FORWARD(s_rhealpix_forward)  sphere
0657                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0658                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
0659                 {
0660                     healpix_sphere(lp_lon, lp_lat, xy_x, xy_y);
0661                     combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 0);
0662                 }
0663
0664                 // INVERSE(s_rhealpix_inverse)  sphere
0665                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0666                 inline void inv(Parameters const& , T xy_x, T xy_y, T& lp_lon, T& lp_lat) const
0667                 {
0668                     /* Check whether (x, y) lies in the rHEALPix image. */
0669                     if (in_image(xy_x, xy_y, 1, this->m_proj_parm.north_square, this->m_proj_parm.south_square) == 0) {
0670                         lp_lon = HUGE_VAL;
0671                         lp_lat = HUGE_VAL;
0672                         BOOST_THROW_EXCEPTION( projection_exception(error_invalid_x_or_y) );
0673                     }
0674                     combine_caps(xy_x, xy_y, this->m_proj_parm.north_square, this->m_proj_parm.south_square, 1);
0675                     return healpix_sphere_inverse(xy_x, xy_y, lp_lon, lp_lat);
0676                 }
0677
0678                 static inline std::string get_name()
0679                 {
0680                     return "rhealpix_spheroid";
0681                 }
0682
0683             };
0684
0685             // HEALPix
0686             template <typename Parameters, typename T>
0687             inline void setup_healpix(Parameters& par, par_healpix<T>& proj_parm)
0688             {
0689                 if (par.es != 0.0) {
0690                     proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */
0691                     proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */
0692                     par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */
0693                     pj_calc_ellipsoid_params(par, par.a, par.es); /* Ensure we have a consistent parameter set */
0694                 } else {
0695                 }
0696             }
0697
0698             // rHEALPix
0699             template <typename Params, typename Parameters, typename T>
0700             inline void setup_rhealpix(Params const& params, Parameters& par, par_healpix<T>& proj_parm)
0701             {
0702                 proj_parm.north_square = pj_get_param_i<srs::spar::north_square>(params, "north_square", srs::dpar::north_square);
0703                 proj_parm.south_square = pj_get_param_i<srs::spar::south_square>(params, "south_square", srs::dpar::south_square);
0704                 /* Check for valid north_square and south_square inputs. */
0705                 if ((proj_parm.north_square < 0) || (proj_parm.north_square > 3)) {
0706                     BOOST_THROW_EXCEPTION( projection_exception(error_axis) );
0707                 }
0708                 if ((proj_parm.south_square < 0) || (proj_parm.south_square > 3)) {
0709                     BOOST_THROW_EXCEPTION( projection_exception(error_axis) );
0710                 }
0711                 if (par.es != 0.0) {
0712                     proj_parm.apa = pj_authset<T>(par.es); /* For auth_lat(). */
0713                     proj_parm.qp = pj_qsfn(1.0, par.e, par.one_es); /* For auth_lat(). */
0714                     par.a = par.a*sqrt(0.5*proj_parm.qp); /* Set par.a to authalic radius. */
0715                     // TODO: why not the same as in healpix?
0716                     //pj_calc_ellipsoid_params(par, par.a, par.es);
0717                     par.ra = 1.0/par.a;
0718                 } else {
0719                 }
0720             }
0721
0722     }} // namespace detail::healpix
0723     #endif // doxygen
0724
0725     /*!
0726         \brief HEALPix projection
0727         \ingroup projections
0728         \tparam Geographic latlong point type
0729         \tparam Cartesian xy point type
0730         \tparam Parameters parameter type
0731         \par Projection characteristics
0732          - Spheroid
0733          - Ellipsoid
0734         \par Example
0735         \image html ex_healpix.gif
0736     */
0737     template <typename T, typename Parameters>
0738     struct healpix_ellipsoid : public detail::healpix::base_healpix_ellipsoid<T, Parameters>
0739     {
0740         template <typename Params>
0741         inline healpix_ellipsoid(Params const& , Parameters & par)
0742         {
0743             detail::healpix::setup_healpix(par, this->m_proj_parm);
0744         }
0745     };
0746
0747     /*!
0748         \brief HEALPix projection
0749         \ingroup projections
0750         \tparam Geographic latlong point type
0751         \tparam Cartesian xy point type
0752         \tparam Parameters parameter type
0753         \par Projection characteristics
0754          - Spheroid
0755          - Ellipsoid
0756         \par Example
0757         \image html ex_healpix.gif
0758     */
0759     template <typename T, typename Parameters>
0760     struct healpix_spheroid : public detail::healpix::base_healpix_spheroid<T, Parameters>
0761     {
0762         template <typename Params>
0763         inline healpix_spheroid(Params const& , Parameters & par)
0764         {
0765             detail::healpix::setup_healpix(par, this->m_proj_parm);
0766         }
0767     };
0768
0769     /*!
0770         \brief rHEALPix projection
0771         \ingroup projections
0772         \tparam Geographic latlong point type
0773         \tparam Cartesian xy point type
0774         \tparam Parameters parameter type
0775         \par Projection characteristics
0776          - Spheroid
0777          - Ellipsoid
0778         \par Projection parameters
0779          - north_square (integer)
0780          - south_square (integer)
0781         \par Example
0782         \image html ex_rhealpix.gif
0783     */
0784     template <typename T, typename Parameters>
0785     struct rhealpix_ellipsoid : public detail::healpix::base_rhealpix_ellipsoid<T, Parameters>
0786     {
0787         template <typename Params>
0788         inline rhealpix_ellipsoid(Params const& params, Parameters & par)
0789         {
0790             detail::healpix::setup_rhealpix(params, par, this->m_proj_parm);
0791         }
0792     };
0793
0794     /*!
0795         \brief rHEALPix projection
0796         \ingroup projections
0797         \tparam Geographic latlong point type
0798         \tparam Cartesian xy point type
0799         \tparam Parameters parameter type
0800         \par Projection characteristics
0801          - Spheroid
0802          - Ellipsoid
0803         \par Projection parameters
0804          - north_square (integer)
0805          - south_square (integer)
0806         \par Example
0807         \image html ex_rhealpix.gif
0808     */
0809     template <typename T, typename Parameters>
0810     struct rhealpix_spheroid : public detail::healpix::base_rhealpix_spheroid<T, Parameters>
0811     {
0812         template <typename Params>
0813         inline rhealpix_spheroid(Params const& params, Parameters & par)
0814         {
0815             detail::healpix::setup_rhealpix(params, par, this->m_proj_parm);
0816         }
0817     };
0818
0819     #ifndef DOXYGEN_NO_DETAIL
0820     namespace detail
0821     {
0822
0823         // Static projection
0824         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_healpix, healpix_spheroid, healpix_ellipsoid)
0825         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI2(srs::spar::proj_rhealpix, rhealpix_spheroid, rhealpix_ellipsoid)
0826
0827         // Factory entry(s)
0828         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(healpix_entry, healpix_spheroid, healpix_ellipsoid)
0829         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI2(rhealpix_entry, rhealpix_spheroid, rhealpix_ellipsoid)
0830
0831         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(healpix_init)
0832         {
0833             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(healpix, healpix_entry)
0834             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(rhealpix, rhealpix_entry)
0835         }
0836
0837     } // namespace detail
0838     #endif // doxygen
0839
0840 } // namespace projections
0841
0842 }} // namespace boost::geometry
0843
0844 #endif // BOOST_GEOMETRY_PROJECTIONS_HEALPIX_HPP
0845