boost/geometry/srs/projections/proj/isea.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 // This code was entirely written by Nathan Wagner
  23 // and is in the public domain.
  24
  25 // Permission is hereby granted, free of charge, to any person obtaining a
  26 // copy of this software and associated documentation files (the "Software"),
  27 // to deal in the Software without restriction, including without limitation
  28 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  29 // and/or sell copies of the Software, and to permit persons to whom the
  30 // Software is furnished to do so, subject to the following conditions:
  31
  32 // The above copyright notice and this permission notice shall be included
  33 // in all copies or substantial portions of the Software.
  34
  35 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  36 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  37 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  38 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  39 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  40 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  41 // DEALINGS IN THE SOFTWARE.
  42
  43 #ifndef BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP
  44 #define BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP
  45
  46 #include <sstream>
  47
  48 #include <boost/core/ignore_unused.hpp>
  49
  50 #include <boost/geometry/core/assert.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_param.hpp>
  56 #include <boost/geometry/srs/projections/impl/projects.hpp>
  57
  58 #include <boost/geometry/util/math.hpp>
  59
  60 namespace boost { namespace geometry
  61 {
  62
  63 namespace projections
  64 {
  65     #ifndef DOXYGEN_NO_DETAIL
  66     namespace detail { namespace isea
  67     {
  68             static const double epsilon = std::numeric_limits<double>::epsilon();
  69
  70             /* sqrt(5)/M_PI */
  71             static const double isea_scale = 0.8301572857837594396028083;
  72             /* 26.565051177 degrees */
  73             static const double v_lat = 0.46364760899944494524;
  74             /* 52.62263186 */
  75             static const double e_rad = 0.91843818702186776133;
  76             /* 10.81231696 */
  77             static const double f_rad = 0.18871053072122403508;
  78             /* R tan(g) sin(60) */
  79             static const double table_g = 0.6615845383;
  80             /* H = 0.25 R tan g = */
  81             static const double table_h = 0.1909830056;
  82             //static const double RPRIME = 0.91038328153090290025;
  83             static const double isea_std_lat = 1.01722196792335072101;
  84             static const double isea_std_lon = .19634954084936207740;
  85
  86             template <typename T>
  87             inline T deg30_rad() { return T(30) * geometry::math::d2r<T>(); }
  88             template <typename T>
  89             inline T deg120_rad() { return T(120) * geometry::math::d2r<T>(); }
  90             template <typename T>
  91             inline T deg72_rad() { return T(72) * geometry::math::d2r<T>(); }
  92             template <typename T>
  93             inline T deg90_rad() { return geometry::math::half_pi<T>(); }
  94             template <typename T>
  95             inline T deg144_rad() { return T(144) * geometry::math::d2r<T>(); }
  96             template <typename T>
  97             inline T deg36_rad() { return T(36) * geometry::math::d2r<T>(); }
  98             template <typename T>
  99             inline T deg108_rad() { return T(108) * geometry::math::d2r<T>(); }
 100             template <typename T>
 101             inline T deg180_rad() { return geometry::math::pi<T>(); }
 102
 103             inline bool downtri(int tri) { return (((tri - 1) / 5) % 2 == 1); }
 104
 105             /*
 106              * Proj 4 provides its own entry points into
 107              * the code, so none of the library functions
 108              * need to be global
 109              */
 110
 111             struct hex {
 112                     int iso;
 113                     int x, y, z;
 114             };
 115
 116             /* y *must* be positive down as the xy /iso conversion assumes this */
 117             inline
 118             int hex_xy(struct hex *h) {
 119                 if (!h->iso) return 1;
 120                 if (h->x >= 0) {
 121                     h->y = -h->y - (h->x+1)/2;
 122                 } else {
 123                     /* need to round toward -inf, not toward zero, so x-1 */
 124                     h->y = -h->y - h->x/2;
 125                 }
 126                 h->iso = 0;
 127
 128                 return 1;
 129             }
 130
 131             inline
 132             int hex_iso(struct hex *h) {
 133                 if (h->iso) return 1;
 134
 135                 if (h->x >= 0) {
 136                     h->y = (-h->y - (h->x+1)/2);
 137                 } else {
 138                     /* need to round toward -inf, not toward zero, so x-1 */
 139                     h->y = (-h->y - (h->x)/2);
 140                 }
 141
 142                 h->z = -h->x - h->y;
 143                 h->iso = 1;
 144                 return 1;
 145             }
 146
 147             template <typename T>
 148             inline
 149             int hexbin2(T const& width, T x, T y, int *i, int *j)
 150             {
 151                 T z, rx, ry, rz;
 152                 T abs_dx, abs_dy, abs_dz;
 153                 int ix, iy, iz, s;
 154                 struct hex h;
 155
 156                 static const T cos_deg30 = cos(deg30_rad<T>());
 157
 158                 x = x / cos_deg30; /* rotated X coord */
 159                 y = y - x / 2.0; /* adjustment for rotated X */
 160
 161                 /* adjust for actual hexwidth */
 162                 x /= width;
 163                 y /= width;
 164
 165                 z = -x - y;
 166
 167                 rx = floor(x + 0.5);
 168                 ix = (int)rx;
 169                 ry = floor(y + 0.5);
 170                 iy = (int)ry;
 171                 rz = floor(z + 0.5);
 172                 iz = (int)rz;
 173
 174                 s = ix + iy + iz;
 175
 176                 if (s) {
 177                     abs_dx = fabs(rx - x);
 178                     abs_dy = fabs(ry - y);
 179                     abs_dz = fabs(rz - z);
 180
 181                     if (abs_dx >= abs_dy && abs_dx >= abs_dz) {
 182                         ix -= s;
 183                     } else if (abs_dy >= abs_dx && abs_dy >= abs_dz) {
 184                         iy -= s;
 185                     } else {
 186                         iz -= s;
 187                     }
 188                 }
 189                 h.x = ix;
 190                 h.y = iy;
 191                 h.z = iz;
 192                 h.iso = 1;
 193
 194                 hex_xy(&h);
 195                 *i = h.x;
 196                 *j = h.y;
 197                 return ix * 100 + iy;
 198             }
 199
 200             //enum isea_poly { isea_none = 0, isea_icosahedron = 20 };
 201             //enum isea_topology { isea_hexagon=6, isea_triangle=3, isea_diamond=4 };
 202             enum isea_address_form {
 203                 /*isea_addr_geo,*/ isea_addr_q2di, isea_addr_seqnum,
 204                 /*isea_addr_interleave,*/ isea_addr_plane, isea_addr_q2dd,
 205                 isea_addr_projtri, isea_addr_vertex2dd, isea_addr_hex
 206             };
 207
 208             template <typename T>
 209             struct isea_dgg {
 210                 T                 o_lat, o_lon, o_az; /* orientation, radians */
 211                 T                 radius; /* radius of the earth in meters, ignored 1.0 */
 212                 unsigned long     serial;
 213                 //int               pole; /* true if standard snyder */
 214                 int               aperture; /* valid values depend on partitioning method */
 215                 int               resolution;
 216                 int               triangle; /* triangle of last transformed point */
 217                 int               quad; /* quad of last transformed point */
 218                 //isea_poly         polyhedron; /* ignored, icosahedron */
 219                 //isea_topology     topology; /* ignored, hexagon */
 220                 isea_address_form output; /* an isea_address_form */
 221             };
 222
 223             template <typename T>
 224             struct isea_pt {
 225                 T x, y;
 226             };
 227
 228             template <typename T>
 229             struct isea_geo {
 230                 T lon, lat;
 231             };
 232
 233             template <typename T>
 234             struct isea_address {
 235                 T      x,y; /* or i,j or lon,lat depending on type */
 236                 int    type; /* enum isea_address_form */
 237                 int    number;
 238             };
 239
 240             /* ENDINC */
 241
 242             enum snyder_polyhedron {
 243                 snyder_poly_hexagon = 0, snyder_poly_pentagon = 1,
 244                 snyder_poly_tetrahedron = 2, snyder_poly_cube = 3,
 245                 snyder_poly_octahedron = 4, snyder_poly_dodecahedron = 5,
 246                 snyder_poly_icosahedron = 6
 247             };
 248
 249             template <typename T>
 250             struct snyder_constants {
 251                 T          g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
 252             };
 253
 254             template <typename T>
 255             inline const snyder_constants<T> * constants()
 256             {
 257                 /* TODO put these in radians to avoid a later conversion */
 258                 static snyder_constants<T> result[] = {
 259                     {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
 260                     {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
 261                     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 262                     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 263                     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 264                     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 265                     {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0}
 266                 };
 267                 return result;
 268             }
 269
 270             template <typename T>
 271             inline const isea_geo<T> * vertex()
 272             {
 273                 static isea_geo<T> result[] = {
 274                     { 0.0,              deg90_rad<T>()},
 275                     { deg180_rad<T>(),  v_lat},
 276                     {-deg108_rad<T>(),  v_lat},
 277                     {-deg36_rad<T>(),   v_lat},
 278                     { deg36_rad<T>(),   v_lat},
 279                     { deg108_rad<T>(),  v_lat},
 280                     {-deg144_rad<T>(), -v_lat},
 281                     {-deg72_rad<T>(),  -v_lat},
 282                     { 0.0,             -v_lat},
 283                     { deg72_rad<T>(),  -v_lat},
 284                     { deg144_rad<T>(), -v_lat},
 285                     { 0.0,             -deg90_rad<T>()}
 286                 };
 287                 return result;
 288             }
 289
 290             /* TODO make an isea_pt array of the vertices as well */
 291
 292             static int      tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};
 293
 294             /* triangle Centers */
 295             template <typename T>
 296             inline const isea_geo<T> * icostriangles()
 297             {
 298                 static isea_geo<T> result[] = {
 299                     { 0.0,              0.0},
 300                     {-deg144_rad<T>(),  e_rad},
 301                     {-deg72_rad<T>(),   e_rad},
 302                     { 0.0,              e_rad},
 303                     { deg72_rad<T>(),   e_rad},
 304                     { deg144_rad<T>(),  e_rad},
 305                     {-deg144_rad<T>(),  f_rad},
 306                     {-deg72_rad<T>(),   f_rad},
 307                     { 0.0,              f_rad},
 308                     { deg72_rad<T>(),   f_rad},
 309                     { deg144_rad<T>(),  f_rad},
 310                     {-deg108_rad<T>(), -f_rad},
 311                     {-deg36_rad<T>(),  -f_rad},
 312                     { deg36_rad<T>(),  -f_rad},
 313                     { deg108_rad<T>(), -f_rad},
 314                     { deg180_rad<T>(), -f_rad},
 315                     {-deg108_rad<T>(), -e_rad},
 316                     {-deg36_rad<T>(),  -e_rad},
 317                     { deg36_rad<T>(),  -e_rad},
 318                     { deg108_rad<T>(), -e_rad},
 319                     { deg180_rad<T>(), -e_rad},
 320                 };
 321                 return result;
 322             }
 323
 324             template <typename T>
 325             inline T az_adjustment(int triangle)
 326             {
 327                 T          adj;
 328
 329                 isea_geo<T> v;
 330                 isea_geo<T> c;
 331
 332                 v = vertex<T>()[tri_v1[triangle]];
 333                 c = icostriangles<T>()[triangle];
 334
 335                 /* TODO looks like the adjustment is always either 0 or 180 */
 336                 /* at least if you pick your vertex carefully */
 337                 adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
 338                         cos(c.lat) * sin(v.lat)
 339                         - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
 340                 return adj;
 341             }
 342
 343             template <typename T>
 344             inline isea_pt<T> isea_triangle_xy(int triangle)
 345             {
 346                 isea_pt<T>  c;
 347                 T Rprime = 0.91038328153090290025;
 348
 349                 triangle = (triangle - 1) % 20;
 350
 351                 c.x = table_g * ((triangle % 5) - 2) * 2.0;
 352                 if (triangle > 9) {
 353                     c.x += table_g;
 354                 }
 355                 switch (triangle / 5) {
 356                 case 0:
 357                     c.y = 5.0 * table_h;
 358                     break;
 359                 case 1:
 360                     c.y = table_h;
 361                     break;
 362                 case 2:
 363                     c.y = -table_h;
 364                     break;
 365                 case 3:
 366                     c.y = -5.0 * table_h;
 367                     break;
 368                 default:
 369                     /* should be impossible */
 370                     BOOST_THROW_EXCEPTION( projection_exception() );
 371                 };
 372                 c.x *= Rprime;
 373                 c.y *= Rprime;
 374
 375                 return c;
 376             }
 377
 378             /* snyder eq 14 */
 379             template <typename T>
 380             inline T sph_azimuth(T const& f_lon, T const& f_lat, T const& t_lon, T const& t_lat)
 381             {
 382                 T          az;
 383
 384                 az = atan2(cos(t_lat) * sin(t_lon - f_lon),
 385                        cos(f_lat) * sin(t_lat)
 386                        - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
 387                     );
 388                 return az;
 389             }
 390
 391             /* coord needs to be in radians */
 392             template <typename T>
 393             inline int isea_snyder_forward(isea_geo<T> * ll, isea_pt<T> * out)
 394             {
 395                 static T const two_pi = detail::two_pi<T>();
 396                 static T const d2r = geometry::math::d2r<T>();
 397
 398                 int             i;
 399
 400                 /*
 401                  * spherical distance from center of polygon face to any of its
 402                  * vertexes on the globe
 403                  */
 404                 T          g;
 405
 406                 /*
 407                  * spherical angle between radius vector to center and adjacent edge
 408                  * of spherical polygon on the globe
 409                  */
 410                 T          G;
 411
 412                 /*
 413                  * plane angle between radius vector to center and adjacent edge of
 414                  * plane polygon
 415                  */
 416                 T          theta;
 417
 418                 /* additional variables from snyder */
 419                 T          q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
 420                                 x, y;
 421
 422                 /* variables used to store intermediate results */
 423                 T          cot_theta, tan_g, az_offset;
 424
 425                 /* how many multiples of 60 degrees we adjust the azimuth */
 426                 int             Az_adjust_multiples;
 427
 428                 snyder_constants<T> c;
 429
 430                 /*
 431                  * TODO by locality of reference, start by trying the same triangle
 432                  * as last time
 433                  */
 434
 435                 /* TODO put these constants in as radians to begin with */
 436                 c = constants<T>()[snyder_poly_icosahedron];
 437                 theta = c.theta * d2r;
 438                 g = c.g * d2r;
 439                 G = c.G * d2r;
 440
 441                 for (i = 1; i <= 20; i++) {
 442                     T          z;
 443                     isea_geo<T> center;
 444
 445                     center = icostriangles<T>()[i];
 446
 447                     /* step 1 */
 448                     z = acos(sin(center.lat) * sin(ll->lat)
 449                          + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
 450
 451                     /* not on this triangle */
 452                     if (z > g + 0.000005) { /* TODO DBL_EPSILON */
 453                         continue;
 454                     }
 455
 456                     Az = sph_azimuth(center.lon, center.lat, ll->lon, ll->lat);
 457
 458                     /* step 2 */
 459
 460                     /* This calculates "some" vertex coordinate */
 461                     az_offset = az_adjustment<T>(i);
 462
 463                     Az -= az_offset;
 464
 465                     /* TODO I don't know why we do this.  It's not in snyder */
 466                     /* maybe because we should have picked a better vertex */
 467                     if (Az < 0.0) {
 468                         Az += two_pi;
 469                     }
 470                     /*
 471                      * adjust Az for the point to fall within the range of 0 to
 472                      * 2(90 - theta) or 60 degrees for the hexagon, by
 473                      * and therefore 120 degrees for the triangle
 474                      * of the icosahedron
 475                      * subtracting or adding multiples of 60 degrees to Az and
 476                      * recording the amount of adjustment
 477                      */
 478
 479                     Az_adjust_multiples = 0;
 480                     while (Az < 0.0) {
 481                         Az += deg120_rad<T>();
 482                         Az_adjust_multiples--;
 483                     }
 484                     while (Az > deg120_rad<T>() + epsilon) {
 485                         Az -= deg120_rad<T>();
 486                         Az_adjust_multiples++;
 487                     }
 488
 489                     /* step 3 */
 490                     cot_theta = 1.0 / tan(theta);
 491                     tan_g = tan(g);    /* TODO this is a constant */
 492
 493                     /* Calculate q from eq 9. */
 494                     /* TODO cot_theta is cot(30) */
 495                     q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
 496
 497                     /* not in this triangle */
 498                     if (z > q + 0.000005) {
 499                         continue;
 500                     }
 501                     /* step 4 */
 502
 503                     /* Apply equations 5-8 and 10-12 in order */
 504
 505                     /* eq 5 */
 506                     /* Rprime = 0.9449322893 * R; */
 507                     /* R' in the paper is for the truncated */
 508                     Rprime = 0.91038328153090290025;
 509
 510                     /* eq 6 */
 511                     H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
 512
 513                     /* eq 7 */
 514                     /* Ag = (Az + G + H - deg180_rad) * M_PI * R * R / deg180_rad; */
 515                     Ag = Az + G + H - deg180_rad<T>();
 516
 517                     /* eq 8 */
 518                     Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
 519
 520                     /* eq 10 */
 521                     /* cot(theta) = 1.73205080756887729355 */
 522                     dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
 523
 524                     /* eq 11 */
 525                     f = dprime / (2.0 * Rprime * sin(q / 2.0));
 526
 527                     /* eq 12 */
 528                     rho = 2.0 * Rprime * f * sin(z / 2.0);
 529
 530                     /*
 531                      * add back the same 60 degree multiple adjustment from step
 532                      * 2 to Azprime
 533                      */
 534
 535                     Azprime += deg120_rad<T>() * Az_adjust_multiples;
 536
 537                     /* calculate rectangular coordinates */
 538
 539                     x = rho * sin(Azprime);
 540                     y = rho * cos(Azprime);
 541
 542                     /*
 543                      * TODO
 544                      * translate coordinates to the origin for the particular
 545                      * hexagon on the flattened polyhedral map plot
 546                      */
 547
 548                     out->x = x;
 549                     out->y = y;
 550
 551                     return i;
 552                 }
 553
 554                 /*
 555                  * should be impossible, this implies that the coordinate is not on
 556                  * any triangle
 557                  */
 558
 559                 //fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
 560                 //    ll->lon * geometry::math::r2d<double>(), ll->lat * geometry::math::r2d<double>());
 561                 std::stringstream ss;
 562                 ss << "impossible transform: " << ll->lon * geometry::math::r2d<T>()
 563                    << " " << ll->lat * geometry::math::r2d<T>() << " is not on any triangle.";
 564
 565                 BOOST_THROW_EXCEPTION( projection_exception(ss.str()) );
 566
 567                 /* not reached */
 568                 return 0;        /* supresses a warning */
 569             }
 570
 571             /*
 572              * return the new coordinates of any point in orginal coordinate system.
 573              * Define a point (newNPold) in orginal coordinate system as the North Pole in
 574              * new coordinate system, and the great circle connect the original and new
 575              * North Pole as the lon0 longitude in new coordinate system, given any point
 576              * in orginal coordinate system, this function return the new coordinates.
 577              */
 578
 579
 580             /* formula from Snyder, Map Projections: A working manual, p31 */
 581             /*
 582              * old north pole at np in new coordinates
 583              * could be simplified a bit with fewer intermediates
 584              *
 585              * TODO take a result pointer
 586              */
 587             template <typename T>
 588             inline isea_geo<T> snyder_ctran(isea_geo<T> * np, isea_geo<T> * pt)
 589             {
 590                 static T const pi = detail::pi<T>();
 591                 static T const two_pi = detail::two_pi<T>();
 592
 593                 isea_geo<T> npt;
 594                 T           alpha, phi, lambda, lambda0, beta, lambdap, phip;
 595                 T           sin_phip;
 596                 T           lp_b;    /* lambda prime minus beta */
 597                 T           cos_p, sin_a;
 598
 599                 phi = pt->lat;
 600                 lambda = pt->lon;
 601                 alpha = np->lat;
 602                 beta = np->lon;
 603                 lambda0 = beta;
 604
 605                 cos_p = cos(phi);
 606                 sin_a = sin(alpha);
 607
 608                 /* mpawm 5-7 */
 609                 sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
 610
 611                 /* mpawm 5-8b */
 612
 613                 /* use the two argument form so we end up in the right quadrant */
 614                 lp_b = atan2(cos_p * sin(lambda - lambda0),
 615                    (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
 616
 617                 lambdap = lp_b + beta;
 618
 619                 /* normalize longitude */
 620                 /* TODO can we just do a modulus ? */
 621                 lambdap = fmod(lambdap, two_pi);
 622                 while (lambdap > pi)
 623                     lambdap -= two_pi;
 624                 while (lambdap < -pi)
 625                     lambdap += two_pi;
 626
 627                 phip = asin(sin_phip);
 628
 629                 npt.lat = phip;
 630                 npt.lon = lambdap;
 631
 632                 return npt;
 633             }
 634
 635             template <typename T>
 636             inline isea_geo<T> isea_ctran(isea_geo<T> * np, isea_geo<T> * pt, T const& lon0)
 637             {
 638                 static T const pi = detail::pi<T>();
 639                 static T const two_pi = detail::two_pi<T>();
 640
 641                 isea_geo<T> npt;
 642
 643                 np->lon += pi;
 644                 npt = snyder_ctran(np, pt);
 645                 np->lon -= pi;
 646
 647                 npt.lon -= (pi - lon0 + np->lon);
 648
 649                 /*
 650                  * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
 651                  * vertex 1 these are 180 degrees apart
 652                  */
 653                 npt.lon += pi;
 654                 /* normalize longitude */
 655                 npt.lon = fmod(npt.lon, two_pi);
 656                 while (npt.lon > pi)
 657                     npt.lon -= two_pi;
 658                 while (npt.lon < -pi)
 659                     npt.lon += two_pi;
 660
 661                 return npt;
 662             }
 663
 664             /* in radians */
 665
 666             /* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
 667
 668             template <typename T>
 669             inline int isea_grid_init(isea_dgg<T> * g)
 670             {
 671                 if (!g)
 672                     return 0;
 673
 674                 //g->polyhedron = isea_icosahedron;
 675                 g->o_lat = isea_std_lat;
 676                 g->o_lon = isea_std_lon;
 677                 g->o_az = 0.0;
 678                 g->aperture = 4;
 679                 g->resolution = 6;
 680                 g->radius = 1.0;
 681                 //g->topology = isea_hexagon;
 682
 683                 return 1;
 684             }
 685
 686             template <typename T>
 687             inline int isea_orient_isea(isea_dgg<T> * g)
 688             {
 689                 if (!g)
 690                     return 0;
 691                 g->o_lat = isea_std_lat;
 692                 g->o_lon = isea_std_lon;
 693                 g->o_az = 0.0;
 694                 return 1;
 695             }
 696
 697             template <typename T>
 698             inline int isea_orient_pole(isea_dgg<T> * g)
 699             {
 700                 static T const half_pi = detail::half_pi<T>();
 701
 702                 if (!g)
 703                     return 0;
 704                 g->o_lat = half_pi;
 705                 g->o_lon = 0.0;
 706                 g->o_az = 0;
 707                 return 1;
 708             }
 709
 710             template <typename T>
 711             inline int isea_transform(isea_dgg<T> * g, isea_geo<T> * in,
 712                                       isea_pt<T> * out)
 713             {
 714                 isea_geo<T> i, pole;
 715                 int         tri;
 716
 717                 pole.lat = g->o_lat;
 718                 pole.lon = g->o_lon;
 719
 720                 i = isea_ctran(&pole, in, g->o_az);
 721
 722                 tri = isea_snyder_forward(&i, out);
 723                 out->x *= g->radius;
 724                 out->y *= g->radius;
 725                 g->triangle = tri;
 726
 727                 return tri;
 728             }
 729
 730
 731             template <typename T>
 732             inline void isea_rotate(isea_pt<T> * pt, T const& degrees)
 733             {
 734                 static T const d2r = geometry::math::d2r<T>();
 735                 static T const two_pi = detail::two_pi<T>();
 736
 737                 T          rad;
 738
 739                 T          x, y;
 740
 741                 rad = -degrees * d2r;
 742                 while (rad >= two_pi) rad -= two_pi;
 743                 while (rad <= -two_pi) rad += two_pi;
 744
 745                 x = pt->x * cos(rad) + pt->y * sin(rad);
 746                 y = -pt->x * sin(rad) + pt->y * cos(rad);
 747
 748                 pt->x = x;
 749                 pt->y = y;
 750             }
 751
 752             template <typename T>
 753             inline int isea_tri_plane(int tri, isea_pt<T> *pt, T const& radius)
 754             {
 755                 isea_pt<T> tc; /* center of triangle */
 756
 757                 if (downtri(tri)) {
 758                     isea_rotate(pt, 180.0);
 759                 }
 760                 tc = isea_triangle_xy<T>(tri);
 761                 tc.x *= radius;
 762                 tc.y *= radius;
 763                 pt->x += tc.x;
 764                 pt->y += tc.y;
 765
 766                 return tri;
 767             }
 768
 769             /* convert projected triangle coords to quad xy coords, return quad number */
 770             template <typename T>
 771             inline int isea_ptdd(int tri, isea_pt<T> *pt)
 772             {
 773                 int             downtri, quad;
 774
 775                 downtri = (((tri - 1) / 5) % 2 == 1);
 776                 quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
 777
 778                 isea_rotate(pt, downtri ? 240.0 : 60.0);
 779                 if (downtri) {
 780                     pt->x += 0.5;
 781                     /* pt->y += cos(30.0 * M_PI / 180.0); */
 782                     pt->y += .86602540378443864672;
 783                 }
 784                 return quad;
 785             }
 786
 787             template <typename T>
 788             inline int isea_dddi_ap3odd(isea_dgg<T> *g, int quad, isea_pt<T> *pt, isea_pt<T> *di)
 789             {
 790                 static T const pi = detail::pi<T>();
 791
 792                 isea_pt<T> v;
 793                 T          hexwidth;
 794                 T          sidelength;    /* in hexes */
 795                 int        d, i;
 796                 int        maxcoord;
 797                 hex        h;
 798
 799                 /* This is the number of hexes from apex to base of a triangle */
 800                 sidelength = (math::pow(T(2), g->resolution) + T(1)) / T(2);
 801
 802                 /* apex to base is cos(30deg) */
 803                 hexwidth = cos(pi / 6.0) / sidelength;
 804
 805                 /* TODO I think sidelength is always x.5, so
 806                  * (int)sidelength * 2 + 1 might be just as good
 807                  */
 808                 maxcoord = (int) (sidelength * 2.0 + 0.5);
 809
 810                 v = *pt;
 811                 hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
 812                 h.iso = 0;
 813                 hex_iso(&h);
 814
 815                 d = h.x - h.z;
 816                 i = h.x + h.y + h.y;
 817
 818                 /*
 819                  * you want to test for max coords for the next quad in the same
 820                  * "row" first to get the case where both are max
 821                  */
 822                 if (quad <= 5) {
 823                     if (d == 0 && i == maxcoord) {
 824                         /* north pole */
 825                         quad = 0;
 826                         d = 0;
 827                         i = 0;
 828                     } else if (i == maxcoord) {
 829                         /* upper right in next quad */
 830                         quad += 1;
 831                         if (quad == 6)
 832                             quad = 1;
 833                         i = maxcoord - d;
 834                         d = 0;
 835                     } else if (d == maxcoord) {
 836                         /* lower right in quad to lower right */
 837                         quad += 5;
 838                         d = 0;
 839                     }
 840                 } else if (quad >= 6) {
 841                     if (i == 0 && d == maxcoord) {
 842                         /* south pole */
 843                         quad = 11;
 844                         d = 0;
 845                         i = 0;
 846                     } else if (d == maxcoord) {
 847                         /* lower right in next quad */
 848                         quad += 1;
 849                         if (quad == 11)
 850                             quad = 6;
 851                         d = maxcoord - i;
 852                         i = 0;
 853                     } else if (i == maxcoord) {
 854                         /* upper right in quad to upper right */
 855                         quad = (quad - 4) % 5;
 856                         i = 0;
 857                     }
 858                 }
 859
 860                 di->x = d;
 861                 di->y = i;
 862
 863                 g->quad = quad;
 864                 return quad;
 865             }
 866
 867             template <typename T>
 868             inline int isea_dddi(isea_dgg<T> *g, int quad, isea_pt<T> *pt, isea_pt<T> *di)
 869             {
 870                 isea_pt<T> v;
 871                 T          hexwidth;
 872                 int        sidelength;    /* in hexes */
 873                 hex        h;
 874
 875                 if (g->aperture == 3 && g->resolution % 2 != 0) {
 876                     return isea_dddi_ap3odd(g, quad, pt, di);
 877                 }
 878                 /* todo might want to do this as an iterated loop */
 879                 if (g->aperture >0) {
 880                     sidelength = (int) (math::pow(T(g->aperture), T(g->resolution / T(2))) + T(0.5));
 881                 } else {
 882                     sidelength = g->resolution;
 883                 }
 884
 885                 hexwidth = 1.0 / sidelength;
 886
 887                 v = *pt;
 888                 isea_rotate(&v, -30.0);
 889                 hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
 890                 h.iso = 0;
 891                 hex_iso(&h);
 892
 893                 /* we may actually be on another quad */
 894                 if (quad <= 5) {
 895                     if (h.x == 0 && h.z == -sidelength) {
 896                         /* north pole */
 897                         quad = 0;
 898                         h.z = 0;
 899                         h.y = 0;
 900                         h.x = 0;
 901                     } else if (h.z == -sidelength) {
 902                         quad = quad + 1;
 903                         if (quad == 6)
 904                             quad = 1;
 905                         h.y = sidelength - h.x;
 906                         h.z = h.x - sidelength;
 907                         h.x = 0;
 908                     } else if (h.x == sidelength) {
 909                         quad += 5;
 910                         h.y = -h.z;
 911                         h.x = 0;
 912                     }
 913                 } else if (quad >= 6) {
 914                     if (h.z == 0 && h.x == sidelength) {
 915                         /* south pole */
 916                         quad = 11;
 917                         h.x = 0;
 918                         h.y = 0;
 919                         h.z = 0;
 920                     } else if (h.x == sidelength) {
 921                         quad = quad + 1;
 922                         if (quad == 11)
 923                             quad = 6;
 924                         h.x = h.y + sidelength;
 925                         h.y = 0;
 926                         h.z = -h.x;
 927                     } else if (h.y == -sidelength) {
 928                         quad -= 4;
 929                         h.y = 0;
 930                         h.z = -h.x;
 931                     }
 932                 }
 933                 di->x = h.x;
 934                 di->y = -h.z;
 935
 936                 g->quad = quad;
 937                 return quad;
 938             }
 939
 940             template <typename T>
 941             inline int isea_ptdi(isea_dgg<T> *g, int tri, isea_pt<T> *pt,
 942                                  isea_pt<T> *di)
 943             {
 944                 isea_pt<T> v;
 945                 int        quad;
 946
 947                 v = *pt;
 948                 quad = isea_ptdd(tri, &v);
 949                 quad = isea_dddi(g, quad, &v, di);
 950                 return quad;
 951             }
 952
 953             /* q2di to seqnum */
 954             template <typename T>
 955             inline int isea_disn(isea_dgg<T> *g, int quad, isea_pt<T> *di)
 956             {
 957                 int             sidelength;
 958                 int             sn, height;
 959                 int             hexes;
 960
 961                 if (quad == 0) {
 962                     g->serial = 1;
 963                     return g->serial;
 964                 }
 965                 /* hexes in a quad */
 966                 hexes = (int) (math::pow(T(g->aperture), T(g->resolution)) + T(0.5));
 967                 if (quad == 11) {
 968                     g->serial = 1 + 10 * hexes + 1;
 969                     return g->serial;
 970                 }
 971                 if (g->aperture == 3 && g->resolution % 2 == 1) {
 972                     height = (int) (math::pow(T(g->aperture), T((g->resolution - 1) / T(2))));
 973                     sn = ((int) di->x) * height;
 974                     sn += ((int) di->y) / height;
 975                     sn += (quad - 1) * hexes;
 976                     sn += 2;
 977                 } else {
 978                     sidelength = (int) (math::pow(T(g->aperture), T(g->resolution / T(2))) + T(0.5));
 979                     sn = (int) ((quad - 1) * hexes + sidelength * di->x + di->y + 2);
 980                 }
 981
 982                 g->serial = sn;
 983                 return sn;
 984             }
 985
 986             /* TODO just encode the quad in the d or i coordinate
 987              * quad is 0-11, which can be four bits.
 988              * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
 989              */
 990             /* convert a q2di to global hex coord */
 991             template <typename T>
 992             inline int isea_hex(isea_dgg<T> *g, int tri, isea_pt<T> *pt,
 993                                 isea_pt<T> *hex)
 994             {
 995                 isea_pt<T> v;
 996 #ifdef BOOST_GEOMETRY_PROJECTIONS_FIXME
 997                 int sidelength;
 998                 int d, i, x, y;
 999 #endif // BOOST_GEOMETRY_PROJECTIONS_FIXME
1000                 int quad;
1001
1002                 quad = isea_ptdi(g, tri, pt, &v);
1003
1004                 hex->x = ((int)v.x << 4) + quad;
1005                 hex->y = v.y;
1006
1007                 return 1;
1008 #ifdef BOOST_GEOMETRY_PROJECTIONS_FIXME
1009                 d = (int)v.x;
1010                 i = (int)v.y;
1011
1012                 /* Aperture 3 odd resolutions */
1013                 if (g->aperture == 3 && g->resolution % 2 != 0) {
1014                     int offset = (int)(pow(T(3.0), T(g->resolution - 1)) + 0.5);
1015
1016                     d += offset * ((g->quad-1) % 5);
1017                     i += offset * ((g->quad-1) % 5);
1018
1019                     if (quad == 0) {
1020                         d = 0;
1021                         i = offset;
1022                     } else if (quad == 11) {
1023                         d = 2 * offset;
1024                         i = 0;
1025                     } else if (quad > 5) {
1026                         d += offset;
1027                     }
1028
1029                     x = (2*d - i) /3;
1030                     y = (2*i - d) /3;
1031
1032                     hex->x = x + offset / 3;
1033                     hex->y = y + 2 * offset / 3;
1034                     return 1;
1035                 }
1036
1037                 /* aperture 3 even resolutions and aperture 4 */
1038                 sidelength = (int) (pow(T(g->aperture), T(g->resolution / 2.0)) + 0.5);
1039                 if (g->quad == 0) {
1040                     hex->x = 0;
1041                     hex->y = sidelength;
1042                 } else if (g->quad == 11) {
1043                     hex->x = sidelength * 2;
1044                     hex->y = 0;
1045                 } else {
1046                     hex->x = d + sidelength * ((g->quad-1) % 5);
1047                     if (g->quad > 5) hex->x += sidelength;
1048                     hex->y = i + sidelength * ((g->quad-1) % 5);
1049                 }
1050
1051                 return 1;
1052 #endif // BOOST_GEOMETRY_PROJECTIONS_FIXME
1053             }
1054
1055             template <typename T>
1056             inline isea_pt<T> isea_forward(isea_dgg<T> *g, isea_geo<T> *in)
1057             {
1058                 int        tri;
1059                 isea_pt<T> out, coord;
1060
1061                 tri = isea_transform(g, in, &out);
1062
1063                 if (g->output == isea_addr_plane) {
1064                     isea_tri_plane(tri, &out, g->radius);
1065                     return out;
1066                 }
1067
1068                 /* convert to isea standard triangle size */
1069                 out.x = out.x / g->radius * isea_scale;
1070                 out.y = out.y / g->radius * isea_scale;
1071                 out.x += 0.5;
1072                 out.y += 2.0 * .14433756729740644112;
1073
1074                 switch (g->output) {
1075                     case isea_addr_projtri:
1076                         /* nothing to do, already in projected triangle */
1077                         break;
1078                     case isea_addr_vertex2dd:
1079                         g->quad = isea_ptdd(tri, &out);
1080                         break;
1081                     case isea_addr_q2dd:
1082                         /* Same as above, we just don't print as much */
1083                         g->quad = isea_ptdd(tri, &out);
1084                         break;
1085                     case isea_addr_q2di:
1086                         g->quad = isea_ptdi(g, tri, &out, &coord);
1087                         return coord;
1088                         break;
1089                     case isea_addr_seqnum:
1090                         isea_ptdi(g, tri, &out, &coord);
1091                         /* disn will set g->serial */
1092                         isea_disn(g, g->quad, &coord);
1093                         return coord;
1094                         break;
1095                     case isea_addr_hex:
1096                         isea_hex(g, tri, &out, &coord);
1097                         return coord;
1098                         break;
1099                     default:
1100                         // isea_addr_plane handled above
1101                         BOOST_GEOMETRY_ASSERT(false);
1102                         break;
1103                 }
1104
1105                 return out;
1106             }
1107             /*
1108              * Proj 4 integration code follows
1109              */
1110
1111             template <typename T>
1112             struct par_isea
1113             {
1114                 isea_dgg<T> dgg;
1115             };
1116
1117             template <typename T, typename Parameters>
1118             struct base_isea_spheroid
1119             {
1120                 par_isea<T> m_proj_parm;
1121
1122                 // FORWARD(s_forward)
1123                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
1124                 inline void fwd(Parameters const& , T const& lp_lon, T const& lp_lat, T& xy_x, T& xy_y) const
1125                 {
1126                     isea_pt<T> out;
1127                     isea_geo<T> in;
1128
1129                     in.lon = lp_lon;
1130                     in.lat = lp_lat;
1131
1132                     isea_dgg<T> copy = this->m_proj_parm.dgg;
1133                     out = isea_forward(&copy, &in);
1134
1135                     xy_x = out.x;
1136                     xy_y = out.y;
1137                 }
1138
1139                 static inline std::string get_name()
1140                 {
1141                     return "isea_spheroid";
1142                 }
1143
1144             };
1145
1146             template <typename T>
1147             inline void isea_orient_init(srs::detail::proj4_parameters const& params,
1148                                          par_isea<T>& proj_parm)
1149             {
1150                 std::string opt = pj_get_param_s(params, "orient");
1151                 if (! opt.empty()) {
1152                     if (opt == std::string("isea")) {
1153                         isea_orient_isea(&proj_parm.dgg);
1154                     } else if (opt == std::string("pole")) {
1155                         isea_orient_pole(&proj_parm.dgg);
1156                     } else {
1157                         BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
1158                     }
1159                 }
1160             }
1161
1162             template <typename T>
1163             inline void isea_orient_init(srs::dpar::parameters<T> const& params,
1164                                          par_isea<T>& proj_parm)
1165             {
1166                 typename srs::dpar::parameters<T>::const_iterator
1167                     it = pj_param_find(params, srs::dpar::orient);
1168                 if (it != params.end()) {
1169                     srs::dpar::value_orient o = static_cast<srs::dpar::value_orient>(it->template get_value<int>());
1170                     if (o == srs::dpar::orient_isea) {
1171                         isea_orient_isea(&proj_parm.dgg);
1172                     } else if (o == srs::dpar::orient_pole) {
1173                         isea_orient_pole(&proj_parm.dgg);
1174                     } else {
1175                         BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
1176                     }
1177                 }
1178             }
1179
1180             template <typename T>
1181             inline void isea_mode_init(srs::detail::proj4_parameters const& params,
1182                                        par_isea<T>& proj_parm)
1183             {
1184                 std::string opt = pj_get_param_s(params, "mode");
1185                 if (! opt.empty()) {
1186                     if (opt == std::string("plane")) {
1187                         proj_parm.dgg.output = isea_addr_plane;
1188                     } else if (opt == std::string("di")) {
1189                         proj_parm.dgg.output = isea_addr_q2di;
1190                     } else if (opt == std::string("dd")) {
1191                         proj_parm.dgg.output = isea_addr_q2dd;
1192                     } else if (opt == std::string("hex")) {
1193                         proj_parm.dgg.output = isea_addr_hex;
1194                     } else {
1195                         BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
1196                     }
1197                 }
1198             }
1199
1200             template <typename T>
1201             inline void isea_mode_init(srs::dpar::parameters<T> const& params,
1202                                        par_isea<T>& proj_parm)
1203             {
1204                 typename srs::dpar::parameters<T>::const_iterator
1205                     it = pj_param_find(params, srs::dpar::mode);
1206                 if (it != params.end()) {
1207                     srs::dpar::value_mode m = static_cast<srs::dpar::value_mode>(it->template get_value<int>());
1208                     if (m == srs::dpar::mode_plane) {
1209                         proj_parm.dgg.output = isea_addr_plane;
1210                     } else if (m == srs::dpar::mode_di) {
1211                         proj_parm.dgg.output = isea_addr_q2di;
1212                     } else if (m == srs::dpar::mode_dd) {
1213                         proj_parm.dgg.output = isea_addr_q2dd;
1214                     } else if (m == srs::dpar::mode_hex) {
1215                         proj_parm.dgg.output = isea_addr_hex;
1216                     } else {
1217                         BOOST_THROW_EXCEPTION( projection_exception(error_ellipsoid_use_required) );
1218                     }
1219                 }
1220             }
1221
1222             // Icosahedral Snyder Equal Area
1223             template <typename Params, typename T>
1224             inline void setup_isea(Params const& params, par_isea<T>& proj_parm)
1225             {
1226                 std::string opt;
1227
1228                 isea_grid_init(&proj_parm.dgg);
1229
1230                 proj_parm.dgg.output = isea_addr_plane;
1231             /*        proj_parm.dgg.radius = par.a; / * otherwise defaults to 1 */
1232                 /* calling library will scale, I think */
1233
1234                 isea_orient_init(params, proj_parm);
1235
1236                 pj_param_r<srs::spar::azi>(params, "azi", srs::dpar::azi, proj_parm.dgg.o_az);
1237                 pj_param_r<srs::spar::lon_0>(params, "lon_0", srs::dpar::lon_0, proj_parm.dgg.o_lon);
1238                 pj_param_r<srs::spar::lat_0>(params, "lat_0", srs::dpar::lat_0, proj_parm.dgg.o_lat);
1239                 // TODO: this parameter is set below second time
1240                 pj_param_i<srs::spar::aperture>(params, "aperture", srs::dpar::aperture, proj_parm.dgg.aperture);
1241                 // TODO: this parameter is set below second time
1242                 pj_param_i<srs::spar::resolution>(params, "resolution", srs::dpar::resolution, proj_parm.dgg.resolution);
1243
1244                 isea_mode_init(params, proj_parm);
1245
1246                 // TODO: pj_param_exists -> pj_get_param_b ?
1247                 if (pj_param_exists<srs::spar::rescale>(params, "rescale", srs::dpar::rescale)) {
1248                     proj_parm.dgg.radius = isea_scale;
1249                 }
1250
1251                 if (pj_param_i<srs::spar::resolution>(params, "resolution", srs::dpar::resolution, proj_parm.dgg.resolution)) {
1252                     /* empty */
1253                 } else {
1254                     proj_parm.dgg.resolution = 4;
1255                 }
1256
1257                 if (pj_param_i<srs::spar::aperture>(params, "aperture", srs::dpar::aperture, proj_parm.dgg.aperture)) {
1258                     /* empty */
1259                 } else {
1260                     proj_parm.dgg.aperture = 3;
1261                 }
1262             }
1263
1264     }} // namespace detail::isea
1265     #endif // doxygen
1266
1267     /*!
1268         \brief Icosahedral Snyder Equal Area projection
1269         \ingroup projections
1270         \tparam Geographic latlong point type
1271         \tparam Cartesian xy point type
1272         \tparam Parameters parameter type
1273         \par Projection characteristics
1274          - Spheroid
1275         \par Projection parameters
1276          - orient (string)
1277          - azi: Azimuth (or Gamma) (degrees)
1278          - lon_0: Central meridian (degrees)
1279          - lat_0: Latitude of origin (degrees)
1280          - aperture (integer)
1281          - resolution (integer)
1282          - mode (string)
1283          - rescale
1284         \par Example
1285         \image html ex_isea.gif
1286     */
1287     template <typename T, typename Parameters>
1288     struct isea_spheroid : public detail::isea::base_isea_spheroid<T, Parameters>
1289     {
1290         template <typename Params>
1291         inline isea_spheroid(Params const& params, Parameters const& )
1292         {
1293             detail::isea::setup_isea(params, this->m_proj_parm);
1294         }
1295     };
1296
1297     #ifndef DOXYGEN_NO_DETAIL
1298     namespace detail
1299     {
1300
1301         // Static projection
1302         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_F(srs::spar::proj_isea, isea_spheroid)
1303
1304         // Factory entry(s)
1305         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_F(isea_entry, isea_spheroid)
1306
1307         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(isea_init)
1308         {
1309             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(isea, isea_entry)
1310         }
1311
1312     } // namespace detail
1313     #endif // doxygen
1314
1315 } // namespace projections
1316
1317 }} // namespace boost::geometry
1318
1319 #endif // BOOST_GEOMETRY_PROJECTIONS_ISEA_HPP
1320