From 825682a2ea70fb923652fb921e3f9122adc0e22f Mon Sep 17 00:00:00 2001 From: Jerome St-Louis Date: Wed, 31 Jul 2024 11:34:41 -0400 Subject: [PATCH 1/4] ISEA Inverse projection (#3047) - Integrating code ported from Java to eC, then eC to C++ implementing the inverse ISEA projection - Originally from Franz-Benjamin Mocnik's ISEA implementation found at https://github.com/mocnik-science/geogrid/blob/master/src/main/java/org/giscience/utils/geogrid/projections/ISEAProjection.java (MIT License) - NOTE: The inverse only supports the default planar options --- docs/source/operations/projections/isea.rst | 34 +- src/projections/isea.cpp | 672 ++++++++++++++++---- test/gie/builtins.gie | 56 +- 3 files changed, 624 insertions(+), 138 deletions(-) diff --git a/docs/source/operations/projections/isea.rst b/docs/source/operations/projections/isea.rst index 73f0af7aca..686fe8fe73 100644 --- a/docs/source/operations/projections/isea.rst +++ b/docs/source/operations/projections/isea.rst @@ -4,18 +4,16 @@ Icosahedral Snyder Equal Area ******************************************************************************** -Snyder's Icosahedral Equal Area map projections on polyhedral globes for the -dodecahedron and truncated icosahedron offer relatively low scale and -angular distortion. The equations involved are relatively straight-forward, -and for certain instructional tools and databases, the projections are useful -for world maps. The interruptions remain a disadvantage, as with any low-error -projection system applied to the entire globe :cite:`Snyder1992`. - +Snyder's Icosahedral Equal Area map projection on an icosahedron polyhedral globe +offers relatively low scale and angular distortion. The equations involved are +relatively straight-forward. The interruptions remain a disadvantage, as with +any low-error projection system applied to the entire globe :cite:`Snyder1992`. +This projection is used as a basis for defining discrete global grid hierarchies. +---------------------+----------------------------------------------------------+ | **Classification** | Polyhedral, equal area | +---------------------+----------------------------------------------------------+ -| **Available forms** | Forward, spherical | +| **Available forms** | Forward and inverse, spherical | +---------------------+----------------------------------------------------------+ | **Defined area** | Global | +---------------------+----------------------------------------------------------+ @@ -36,6 +34,14 @@ projection system applied to the entire globe :cite:`Snyder1992`. proj-string: ``+proj=isea`` +.. note:: + As the projection is only defined on a sphere, it should only be used + with a spherical ellipsoid e.g., ``+R=6371007.18091875`` for a sphere with the + authalic radius of the WGS84 ellipsoid. For mapping coordinates on the WGS84 + ellipsoid to the authalic sphere, the input latitude should be converted + from geodetic latitude to authalic latitude. A future version may + automatically perform this conversion when using a non-spherical ellipsoid. + Parameters ################################################################################ @@ -44,27 +50,35 @@ Parameters .. option:: +orient= Can be set to either ``isea`` or ``pole``. See Snyder's Figure 12 for pole orientation :cite:`Snyder1992`. - + *Defaults to isea* .. option:: +azi= Azimuth. + Not supported by the inverse. + *Defaults to 0.0* .. option:: +aperture= + Not supported by the inverse. + *Defaults to 3.0* .. option:: +resolution= + Not supported by the inverse. + *Defaults to 4.0* .. option:: +mode= Can be either ``plane``, ``di``, ``dd`` or ``hex``. - + + Only ``plane`` supported by the inverse. + *Defaults to plane* .. include:: ../options/lon_0.rst diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp index 8ff5009edc..6f8ba66e49 100644 --- a/src/projections/isea.cpp +++ b/src/projections/isea.cpp @@ -1,6 +1,57 @@ /* - * This code was entirely written by Nathan Wagner - * and is in the public domain. + The public domain code for the forward direction was initially + written by Nathan Wagner. + + The inverse projection was adapted from Java and eC by + Jérôme Jacovella-St-Louis, originally from the Franz-Benjamin Mocnik's ISEA + implementation found at + https://github.com/mocnik-science/geogrid/blob/master/ + src/main/java/org/giscience/utils/geogrid/projections/ISEAProjection.java + with the following license: + -------------------------------------------------------------------------- + MIT License + + Copyright (c) 2017-2019 Heidelberg University + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in + all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ + +/* The ISEA projection a projects a sphere on the icosahedron. Thereby the size + * of areas mapped to the icosahedron are preserved. Angles and distances are + * however slightly distorted. The angular distortion is below 17.27 degree, and + * the scale variation is less than 16.3 per cent. + * + * The projection has been proposed and has been described in detail by: + * + * John P. Snyder: An equal-area map projection for polyhedral globes. + * Cartographica, 29(1), 10–21, 1992. doi:10.3138/27H7-8K88-4882-1752 + * + * Another description and improvements can be found in: + * + * Erika Harrison, Ali Mahdavi-Amiri, and Faramarz Samavati: Optimization of + * inverse Snyder polyhedral projection. International Conference on Cyberworlds + * 2011. doi:10.1109/CW.2011.36 + * + * Erika Harrison, Ali Mahdavi-Amiri, and Faramarz Samavati: Analysis of inverse + * Snyder optimizations. In: Marina L. Gavrilova, and C. J. Kenneth Tan (Eds): + * Transactions on Computational Science XVI. Heidelberg, Springer, 2012. pp. + * 134–148. doi:10.1007/978-3-642-32663-9_8 */ #include @@ -10,6 +61,7 @@ #include #include +#include #include #include "proj.h" @@ -30,17 +82,52 @@ /* 26.565051177 degrees */ #define V_LAT 0.46364760899944494524 -/* 52.62263186 */ -#define E_RAD 0.91843818702186776133 +// Latitude of center of top icosahedron faces +// atan((3 + sqrt(5)) / 4) = 52.6226318593487 degrees +#define E_RAD 0.91843818701052843323 + +// Latitude of center of faces mirroring top icosahedron face +// atan((3 - sqrt(5)) / 4) = 10.8123169635739 degrees +#define F_RAD 0.18871053078356206978 + +// #define phi ((1 + sqrt(5)) / 2) +// #define atanphi 1.01722196789785136772 + +// g: Spherical distance from center of polygon face to +// any of its vertices on the sphere +// g = F + 2 * atan(phi) - 90 deg -- sdc2vos +#define sdc2vos 0.6523581397843681859886783 -/* 10.81231696 */ -#define F_RAD 0.18871053072122403508 +#define tang 0.76393202250021030358019673567 // tan(sdc2vos) -/* R tan(g) sin(60) */ -#define TABLE_G 0.6615845383 +// theta (30 degrees) is plane angle between radius +// vector to center and adjacent edge of plane polygon + +#define tan30 0.57735026918962576450914878 // tan(DEG_TO_RAD * 30) +#define cotTheta (1.0 / tan30) + +// G: spherical angle between radius vector to center and adjacent edge +// of spherical polygon on the globe (36 degrees) +// cos(DEG_TO_RAD * 36) +#define cosG 0.80901699437494742410229341718281905886 +// sin(DEG_TO_RAD * 36) +#define sinG 0.587785252292473129168705954639072768597652 +// cos(g) +#define cosSDC2VoS 0.7946544722917661229596057297879189448539 + +#define sinGcosSDC2VoS (sinG * cosSDC2VoS) // sin G cos g + +#define SQRT3 1.73205080756887729352744634150587236694280525381038 +#define sin60 (SQRT3 / 2.0) + +// tang * sin(60 deg) +#define TABLE_G (tang * sin60) + +// (1 / (2 * sqrt(5)) + 1 / 6.0) * sqrt(M_PI * sqrt(3)) +#define RprimeOverR 0.9103832815095032 // R' / R /* H = 0.25 R tan g = */ -#define TABLE_H 0.1909830056 +#define TABLE_H (0.25 * tang) /* in radians */ #define ISEA_STD_LAT 1.01722196792335072101 @@ -135,9 +222,9 @@ static void hexbin2(double width, double x, double y, long *i, long *j) { *j = h.y; } +#define numIcosahedronFaces 20 + namespace { // anonymous namespace -enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 }; -enum isea_topology { ISEA_HEXAGON = 6, ISEA_TRIANGLE = 3, ISEA_DIAMOND = 4 }; enum isea_address_form { ISEA_GEO, ISEA_Q2DI, @@ -149,9 +236,7 @@ enum isea_address_form { ISEA_VERTEX2DD, ISEA_HEX }; -} // anonymous namespace -namespace { // anonymous namespace struct isea_dgg { int polyhedron; /* ignored, icosahedron */ double o_lat, o_lon, o_az; /* orientation, radians */ @@ -164,23 +249,22 @@ struct isea_dgg { int quad; /* quad of last transformed point */ unsigned long serial; }; -} // anonymous namespace -namespace { // anonymous namespace struct isea_pt { double x, y; }; -} // anonymous namespace -namespace { // anonymous namespace struct isea_geo { double longitude, lat; }; + } // anonymous namespace -/* ENDINC */ +/* No longer used + +enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = numIcosahedronFaces }; +enum isea_topology { ISEA_HEXAGON = 6, ISEA_TRIANGLE = 3, ISEA_DIAMOND = 4 }; -namespace { // anonymous namespace enum snyder_polyhedron { SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON, @@ -190,17 +274,14 @@ enum snyder_polyhedron { SNYDER_POLY_DODECAHEDRON, SNYDER_POLY_ICOSAHEDRON }; -} // anonymous namespace -namespace { // anonymous namespace struct snyder_constants { double g, G, theta; - /* cppcheck-suppress unusedStructMember */ + // cppcheck-suppress unusedStructMember double ea_w, ea_a, ea_b, g_w, g_a, g_b; }; -} // anonymous namespace -/* TODO put these in radians to avoid a later conversion */ +// TODO put these in radians to avoid a later conversion static const struct snyder_constants constants[] = { {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0}, {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027}, @@ -210,6 +291,7 @@ static const struct snyder_constants constants[] = { {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0}, }; +*/ static struct isea_geo vertex[] = { {0.0, DEG90}, {DEG180, V_LAT}, {-DEG108, V_LAT}, {-DEG36, V_LAT}, @@ -250,9 +332,8 @@ static double az_adjustment(int triangle) { static struct isea_pt isea_triangle_xy(int triangle) { struct isea_pt c; - const double Rprime = 0.91038328153090290025; - triangle = (triangle - 1) % 20; + triangle = (triangle - 1) % numIcosahedronFaces; c.x = TABLE_G * ((triangle % 5) - 2) * 2.0; if (triangle > 9) { @@ -275,8 +356,8 @@ static struct isea_pt isea_triangle_xy(int triangle) { /* should be impossible */ exit(EXIT_FAILURE); } - c.x *= Rprime; - c.y *= Rprime; + c.x *= RprimeOverR; + c.y *= RprimeOverR; return c; } @@ -302,58 +383,29 @@ static double sph_azimuth(double f_lon, double f_lat, double t_lon, static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { int i; - /* - * spherical distance from center of polygon face to any of its - * vertices on the globe - */ - double g; - - /* - * spherical angle between radius vector to center and adjacent edge - * of spherical polygon on the globe - */ - double G; - - /* - * plane angle between radius vector to center and adjacent edge of - * plane polygon - */ - double theta; - /* additional variables from snyder */ double q, H, Ag, Azprime, Az, dprime, f, rho, x, y; /* variables used to store intermediate results */ - double cot_theta, tan_g, az_offset; + double az_offset; /* how many multiples of 60 degrees we adjust the azimuth */ int Az_adjust_multiples; - struct snyder_constants c; - /* * TODO by locality of reference, start by trying the same triangle * as last time */ - - /* TODO put these constants in as radians to begin with */ - c = constants[SNYDER_POLY_ICOSAHEDRON]; - theta = PJ_TORAD(c.theta); - g = PJ_TORAD(c.g); - G = PJ_TORAD(c.G); - - for (i = 1; i <= 20; i++) { - double z; - struct isea_geo center; - - center = icostriangles[i]; + for (i = 1; i <= numIcosahedronFaces; i++) { + struct isea_geo center = icostriangles[i]; + double cosZ = sin(center.lat) * sin(ll->lat) + + cos(center.lat) * cos(ll->lat) * + cos(ll->longitude - center.longitude); /* step 1 */ - z = acos(sin(center.lat) * sin(ll->lat) + - cos(center.lat) * cos(ll->lat) * - cos(ll->longitude - center.longitude)); + double z = fabs(cosZ) < 1E-15 ? 0 : acos(cosZ); /* not on this triangle */ - if (z > g + 0.000005) { /* TODO DBL_EPSILON */ + if (z > sdc2vos /*g*/ + 0.000005) { /* TODO DBL_EPSILON */ continue; } @@ -391,12 +443,9 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { } /* step 3 */ - cot_theta = 1.0 / tan(theta); - tan_g = tan(g); /* TODO this is a constant */ /* Calculate q from eq 9. */ - /* TODO cot_theta is cot(30) */ - q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta); + q = atan2(tang, cos(Az) + sin(Az) * cotTheta); /* not in this triangle */ if (z > q + 0.000005) { @@ -407,30 +456,29 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { /* Apply equations 5-8 and 10-12 in order */ /* eq 5 */ - /* Rprime = 0.9449322893 * R; */ - /* R' in the paper is for the truncated */ - const double Rprime = 0.91038328153090290025; + /* R' in the paper is for the truncated (icosahedron?) */ /* eq 6 */ - H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G)); + H = acos(sin(Az) * sinGcosSDC2VoS /* sin(G) * cos(g) */ - + cos(Az) * cosG); /* eq 7 */ /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */ - Ag = Az + G + H - DEG180; + Ag = Az + DEG_TO_RAD * 36 /* G */ + H - DEG180; /* eq 8 */ - Azprime = atan2(2.0 * Ag, - Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta); + Azprime = atan2(2.0 * Ag, RprimeOverR * RprimeOverR * tang * tang - + 2.0 * Ag * cotTheta); /* eq 10 */ /* cot(theta) = 1.73205080756887729355 */ - dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta); + dprime = RprimeOverR * tang / (cos(Azprime) + sin(Azprime) * cotTheta); /* eq 11 */ - f = dprime / (2.0 * Rprime * sin(q / 2.0)); + f = dprime / (2.0 * RprimeOverR * sin(q / 2.0)); /* eq 12 */ - rho = 2.0 * Rprime * f * sin(z / 2.0); + rho = 2.0 * RprimeOverR * f * sin(z / 2.0); /* * add back the same 60 degree multiple adjustment from step @@ -486,33 +534,25 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { * * TODO take a result pointer */ -static struct isea_geo snyder_ctran(struct isea_geo *np, struct isea_geo *pt) { - struct isea_geo npt; - double alpha, phi, lambda, lambda0, beta, lambdap, phip; - double sin_phip; - double lp_b; /* lambda prime minus beta */ - double cos_p, sin_a; - - phi = pt->lat; - lambda = pt->longitude; - alpha = np->lat; - beta = np->longitude; - lambda0 = beta; - - cos_p = cos(phi); - sin_a = sin(alpha); +static struct isea_geo snyder_ctran(const struct isea_geo &np, + const struct isea_geo &pt) { + struct isea_geo result; + double phi = pt.lat, lambda = pt.longitude; + double alpha = np.lat, beta = np.longitude; + double dlambda = lambda - beta /* lambda0 */; + double cos_p = cos(phi), sin_p = sin(phi); + double cos_a = cos(alpha), sin_a = sin(alpha); + double cos_dlambda = cos(dlambda), sin_dlambda = sin(dlambda); /* mpawm 5-7 */ - sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0); + double sin_phip = sin_a * sin_p - cos_a * cos_p * cos_dlambda; /* mpawm 5-8b */ /* use the two argument form so we end up in the right quadrant */ - lp_b = - atan2(cos_p * sin(lambda - lambda0), - (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi))); - - lambdap = lp_b + beta; + double lp_b = /* lambda prime minus beta */ + atan2(cos_p * sin_dlambda, sin_a * cos_p * cos_dlambda + cos_a * sin_p); + double lambdap = lp_b + beta; /* normalize longitude */ /* TODO can we just do a modulus ? */ @@ -522,29 +562,25 @@ static struct isea_geo snyder_ctran(struct isea_geo *np, struct isea_geo *pt) { while (lambdap < -M_PI) lambdap += 2 * M_PI; - phip = asin(sin_phip); - - npt.lat = phip; - npt.longitude = lambdap; - - return npt; + result.lat = fabs(sin_phip - 1.0) < 1E-15 ? M_PI / 2 + : fabs(sin_phip + 1.0) < 1E-15 ? -M_PI / 2 + : asin(sin_phip); + result.longitude = lambdap; + return result; } -static struct isea_geo isea_ctran(struct isea_geo *np, struct isea_geo *pt, - double lon0) { - struct isea_geo npt; - - np->longitude += M_PI; - npt = snyder_ctran(np, pt); - np->longitude -= M_PI; - - npt.longitude -= (M_PI - lon0 + np->longitude); +static struct isea_geo isea_ctran(const struct isea_geo *np, + const struct isea_geo *pt, double lon0) { + struct isea_geo cnp = {np->longitude + M_PI, np->lat}; + struct isea_geo npt = snyder_ctran(cnp, *pt); + npt.longitude -= (/* M_PI */ -lon0 + np->longitude); /* * snyder is down tri 3, isea is along side of tri1 from vertex 0 to * vertex 1 these are 180 degrees apart */ - npt.longitude += M_PI; + // npt.longitude += M_PI; + /* normalize longitude */ npt.longitude = fmod(npt.longitude, 2 * M_PI); while (npt.longitude > M_PI) @@ -561,7 +597,7 @@ static int isea_grid_init(struct isea_dgg *g) { if (!g) return 0; - g->polyhedron = 20; + g->polyhedron = numIcosahedronFaces; g->o_lat = ISEA_STD_LAT; g->o_lon = ISEA_STD_LONG; g->o_az = 0.0; @@ -984,8 +1020,33 @@ static struct isea_pt isea_forward(struct isea_dgg *g, struct isea_geo *in) { PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph"; namespace { // anonymous namespace + +struct GeoPoint { + double lat, lon; +}; // In radians + +class ISEAPlanarProjection; + +struct isea_inverse_params { + double R2; + double Rprime; + double Rprime2X; + double RprimeTang; + double Rprime2Tan2g; + double triTang; + double centerToBase; + double triWidth; + double yOffsets[4]; + double xo, yo; + double sx, sy; + ISEAPlanarProjection *p; + + void initialize(const PJ *P); +}; + struct pj_isea_data { struct isea_dgg dgg; + struct isea_inverse_params inv; }; } // anonymous namespace @@ -999,6 +1060,8 @@ static PJ_XY isea_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ in.lat = lp.phi; try { + // TODO: Convert geodetic latitude to authalic latitude if not + // spherical as in eqearth, healpix, laea, etc. out = isea_forward(&Q->dgg, &in); } catch (const char *) { proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); @@ -1011,6 +1074,8 @@ static PJ_XY isea_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ return xy; } +static PJ_LP isea_s_inverse(PJ_XY xy, PJ *P); + PJ *PJ_PROJECTION(isea) { char *opt; struct pj_isea_data *Q = static_cast( @@ -1022,9 +1087,10 @@ PJ *PJ_PROJECTION(isea) { // NOTE: if a inverse was needed, there is some material at // https://brsr.github.io/2021/08/31/snyder-equal-area.html P->fwd = isea_s_forward; + P->inv = isea_s_inverse; isea_grid_init(&Q->dgg); - Q->dgg.output = ISEA_PLANE; + /* P->dgg.radius = P->a; / * otherwise defaults to 1 */ /* calling library will scale, I think */ @@ -1089,9 +1155,375 @@ PJ *PJ_PROJECTION(isea) { Q->dgg.aperture = 3; } + Q->inv.initialize(P); + return P; } +#define Min std::min +#define Max std::max + +#define inf std::numeric_limits::infinity() + +// distortion +// static double maximumAngularDistortion = 17.27; +// static double maximumScaleVariation = 1.163; +// static double minimumScaleVariation = .860; + +// Vertices of dodecahedron centered in icosahedron faces +static const GeoPoint facesCenterDodecahedronVertices[numIcosahedronFaces] = { + {E_RAD, DEG_TO_RAD * -144}, {E_RAD, DEG_TO_RAD * -72}, + {E_RAD, DEG_TO_RAD * 0}, {E_RAD, DEG_TO_RAD * 72}, + {E_RAD, DEG_TO_RAD * 144}, {F_RAD, DEG_TO_RAD * -144}, + {F_RAD, DEG_TO_RAD * -72}, {F_RAD, DEG_TO_RAD * 0}, + {F_RAD, DEG_TO_RAD * 72}, {F_RAD, DEG_TO_RAD * 144}, + {-F_RAD, DEG_TO_RAD * -108}, {-F_RAD, DEG_TO_RAD * -36}, + {-F_RAD, DEG_TO_RAD * 36}, {-F_RAD, DEG_TO_RAD * 108}, + {-F_RAD, DEG_TO_RAD * 180}, {-E_RAD, DEG_TO_RAD * -108}, + {-E_RAD, DEG_TO_RAD * -36}, {-E_RAD, DEG_TO_RAD * 36}, + {-E_RAD, DEG_TO_RAD * 108}, {-E_RAD, DEG_TO_RAD * 180}}; + +// static define precision = DEG_TO_RAD * 1e-9; +#define precision (DEG_TO_RAD * 1e-11) +#define precisionPerDefinition (DEG_TO_RAD * 1e-5) + +#define AzMax (DEG_TO_RAD * 120) + +#define westVertexLon (DEG_TO_RAD * -144) + +namespace { // anonymous namespace + +struct ISEAFacePoint { + int face; + double x, y; +}; + +class ISEAPlanarProjection { + public: + explicit ISEAPlanarProjection(const GeoPoint &value) + : orientation(value), cosOrientationLat(cos(value.lat)), + sinOrientationLat(sin(value.lat)) {} + + bool cartesianToGeo(const PJ_XY &inPosition, + const isea_inverse_params ¶ms, GeoPoint &result) { + bool r = false; + static const double epsilon = + 2E-8; // NOTE: 1E-11 seems too small for forward projection + // precision at boundaries + int face = 0; + PJ_XY position = inPosition; + +#define sr -sin60 // sin(-60) +#define cr 0.5 // cos(-60) + if (position.x < 0 || + (position.x < params.triWidth / 2 && position.y < 0 && + position.y * cr < position.x * sr)) + position.x += 5 * params.triWidth; // Wrap around +// Rotate and shear to determine face if not stored in position.z +#define shearX (1.0 / SQRT3) + double yp = -(position.x * sr + position.y * cr); + double x = + (position.x * cr - position.y * sr + yp * shearX) * params.sx; + double y = yp * params.sy; +#undef shearX +#undef sr +#undef cr + + if (x < 0 || (y > x && x < 5 - epsilon)) + x += epsilon; + else if (x > 5 || (y < x && x > 0 + epsilon)) + x -= epsilon; + + if (y < 0 || (x > y && y < 6 - epsilon)) + y += epsilon; + else if (y > 6 || (x < y && y > 0 + epsilon)) + y -= epsilon; + + if (x >= 0 && x <= 5 && y >= 0 && y <= 6) { + int ix = Max(0, Min(4, (int)x)); + int iy = Max(0, Min(5, (int)y)); + + if (iy == ix || iy == ix + 1) { + int rhombus = ix + iy; + bool top = x - ix > y - iy; + face = -1; + + switch (rhombus) { + case 0: + face = top ? 0 : 5; + break; + case 2: + face = top ? 1 : 6; + break; + case 4: + face = top ? 2 : 7; + break; + case 6: + face = top ? 3 : 8; + break; + case 8: + face = top ? 4 : 9; + break; + + case 1: + face = top ? 10 : 15; + break; + case 3: + face = top ? 11 : 16; + break; + case 5: + face = top ? 12 : 17; + break; + case 7: + face = top ? 13 : 18; + break; + case 9: + face = top ? 14 : 19; + break; + } + face++; + } + } + + if (face) { + int fy = (face - 1) / 5, fx = (face - 1) - 5 * fy; + double rx = + position.x - (2 * fx + fy / 2 + 1) * params.triWidth / 2; + double ry = + position.y - (params.yOffsets[fy] + 3 * params.centerToBase); + GeoPoint dst; + + r = icosahedronToSphere({face - 1, rx, ry}, params, dst); + + if (dst.lon < -M_PI - epsilon) + dst.lon += 2 * M_PI; + else if (dst.lon > M_PI + epsilon) + dst.lon -= 2 * M_PI; + + result = {dst.lat, dst.lon}; + } + return r; + } + + // Converts coordinates on the icosahedron to geographic coordinates + // (inverse projection) + bool icosahedronToSphere(const ISEAFacePoint &c, + const isea_inverse_params ¶ms, GeoPoint &r) { + if (c.face >= 0 && c.face < numIcosahedronFaces) { + double Az = atan2(c.x, c.y); // Az' + double rho = sqrt(c.x * c.x + c.y * c.y); + double AzAdjustment = faceOrientation(c.face); + + Az += AzAdjustment; + while (Az < 0) { + AzAdjustment += AzMax; + Az += AzMax; + } + while (Az > AzMax) { + AzAdjustment -= AzMax; + Az -= AzMax; + } + { + double sinAz = sin(Az), cosAz = cos(Az); + double cotAz = cosAz / sinAz; + double area = params.Rprime2Tan2g / + (2 * (cotAz + cotTheta)); // A_G or A_{ABD} + double deltaAz = 10 * precision; + double degAreaOverR2Plus180Minus36 = + area / params.R2 - westVertexLon; + double Az_earth = Az; + + while (fabs(deltaAz) > precision) { + double sinAzEarth = sin(Az_earth), + cosAzEarth = cos(Az_earth); + double H = + acos(sinAzEarth * sinGcosSDC2VoS - cosAzEarth * cosG); + double FAz_earth = degAreaOverR2Plus180Minus36 - H - + Az_earth; // F(Az) or g(Az) + double F2Az_earth = + (cosAzEarth * sinGcosSDC2VoS + sinAzEarth * cosG) / + sin(H) - + 1; // F'(Az) or g'(Az) + deltaAz = -FAz_earth / F2Az_earth; // Delta Az^0 or Delta Az + Az_earth += deltaAz; + } + { + double sinAz_earth = sin(Az_earth), + cosAz_earth = cos(Az_earth); + double q = + atan2(tang, (cosAz_earth + sinAz_earth * cotTheta)); + double d = + params.RprimeTang / (cosAz + sinAz * cotTheta); // d' + double f = d / (params.Rprime2X * sin(q / 2)); // f + double z = 2 * asin(rho / (params.Rprime2X * f)); + + Az_earth -= AzAdjustment; + { + double lat0 = + facesCenterDodecahedronVertices[c.face].lat; + double sinLat0 = sin(lat0), cosLat0 = cos(lat0); + double sinZ = sin(z), cosZ = cos(z); + double cosLat0SinZ = cosLat0 * sinZ; + double latSin = + sinLat0 * cosZ + cosLat0SinZ * cos(Az_earth); + double lat = fabs(latSin - 1.0) < 1E-15 ? M_PI / 2 + : fabs(latSin + 1.0) < 1E-15 + ? -M_PI / 2 + : asin(latSin); + double lon = + facesCenterDodecahedronVertices[c.face].lon + + atan2(sin(Az_earth) * cosLat0SinZ, + cosZ - sinLat0 * sin(lat)); + + revertOrientation({lat, lon}, r); + } + } + } + return true; + } + r = {inf, inf}; + return false; + } + + private: + GeoPoint orientation; + double cosOrientationLat, sinOrientationLat; + + inline void revertOrientation(const GeoPoint &c, GeoPoint &r) { + double lon = (c.lat < DEG_TO_RAD * -90 + precisionPerDefinition || + c.lat > DEG_TO_RAD * 90 - precisionPerDefinition) + ? 0 + : c.lon; + if (orientation.lat != 0.0 || orientation.lon != 0.0) { + double sinLat = sin(c.lat), cosLat = cos(c.lat); + double sinLon = sin(lon), cosLon = cos(lon); + double cosLonCosLat = cosLon * cosLat; + r = {asin(sinLat * cosOrientationLat - + cosLonCosLat * sinOrientationLat), + atan2(sinLon * cosLat, cosLonCosLat * cosOrientationLat + + sinLat * sinOrientationLat) - + orientation.lon}; + } else + r = {c.lat, lon}; + } + + static inline double faceOrientation(int face) { + return (face <= 4 || (10 <= face && face <= 14)) ? 0 : DEG_TO_RAD * 180; + } +}; + +// Orientation symmetric to equator (+proj=isea) +/* Sets the orientation of the icosahedron such that the north and the south + * poles are mapped to the edge midpoints of the icosahedron. The equator is + * thus mapped symmetrically. + */ +static ISEAPlanarProjection standardISEA( + /* DEG_TO_RAD * (90 - 58.282525589) = 31.7174744114613 */ + {(E_RAD + F_RAD) / 2, DEG_TO_RAD * -11.25}); + +// Polar orientation (+proj=isea +orient=pole) +/* + * One corner of the icosahedron is, by default, facing to the north pole, and + * one to the south pole. The provided orientation is relative to the default + * orientation. + * + * The orientation shifts every location by the angle orientation.lon in + * direction of positive longitude, and thereafter by the angle orientation.lat + * in direction of positive latitude. + */ +static ISEAPlanarProjection polarISEA({0, 0}); + +void isea_inverse_params::initialize(const PJ *P) { + struct pj_isea_data *Q = static_cast(P->opaque); + // Only supporting default planar options for now + if (Q->dgg.output == ISEA_PLANE && Q->dgg.o_az == 0.0 && + Q->dgg.aperture == 3.0 && Q->dgg.resolution == 4.) { + // Only supporting +orient=isea and +orient=pole for now + if (Q->dgg.o_lat == ISEA_STD_LAT && Q->dgg.o_lon == ISEA_STD_LONG) + p = &standardISEA; + else if (Q->dgg.o_lat == M_PI / 2.0 && Q->dgg.o_lon == 0) + p = &polarISEA; + else + p = nullptr; + } + + if (p != nullptr) { + if (P->e > 0) { + double a2 = P->a * P->a, c2 = P->b * P->b; + double log1pe_1me = log((1 + P->e) / (1 - P->e)); + double S = M_PI * (2 * a2 + c2 / P->e * log1pe_1me); + R2 = S / (4 * M_PI); // [WGS84] R = 6371007.1809184747 m + Rprime = RprimeOverR * sqrt(R2); // R' + } else { + R2 = P->a * P->a; // R^2 + Rprime = RprimeOverR * P->a; // R' + } + + Rprime2X = 2 * Rprime; + RprimeTang = Rprime * tang; // twice the center-to-base distance + centerToBase = RprimeTang / 2; + triWidth = RprimeTang * SQRT3; + Rprime2Tan2g = RprimeTang * RprimeTang; + + yOffsets[0] = -2 * centerToBase; + yOffsets[1] = -4 * centerToBase; + yOffsets[2] = -5 * centerToBase; + yOffsets[3] = -7 * centerToBase; + + xo = 2.5 * triWidth; + yo = -1.5 * centerToBase; + sx = 1.0 / triWidth; + sy = 1.0 / (3 * centerToBase); + } +} + +} // anonymous namespace + +static PJ_LP isea_s_inverse(PJ_XY xy, PJ *P) { + struct pj_isea_data *Q = static_cast(P->opaque); + ISEAPlanarProjection *p = Q->inv.p; + + if (p) { + // Default origin of +proj=isea is different (OGC:1534 is + // +x_0=19186144.870934911 +y_0=-3323137.7717836285) + PJ_XY input{xy.x * P->a + Q->inv.xo, xy.y * P->a + Q->inv.yo}; + GeoPoint result; + + if (p->cartesianToGeo(input, Q->inv, result)) + // TODO: Convert authalic latitude to geodetic latitude if not + // spherical as in eqearth, healpix, laea, etc. + return {result.lon, result.lat}; + else + return {inf, inf}; + } else + return {inf, inf}; +} + +#undef ISEA_STD_LAT +#undef ISEA_STD_LONG + +#undef numIcosahedronFaces +#undef precision +#undef precisionPerDefinition + +#undef AzMax +#undef sdc2vos +#undef tang +#undef cotTheta +#undef cosG +#undef sinGcosSDC2VoS +#undef westVertexLon + +#undef RprimeOverR + +#undef Min +#undef Max + +#undef inf + +#undef E_RAD +#undef F_RAD + #undef DEG36 #undef DEG72 #undef DEG90 @@ -1101,9 +1533,5 @@ PJ *PJ_PROJECTION(isea) { #undef DEG180 #undef ISEA_SCALE #undef V_LAT -#undef E_RAD -#undef F_RAD #undef TABLE_G #undef TABLE_H -#undef ISEA_STD_LAT -#undef ISEA_STD_LONG diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 98ea577449..c500c3fb04 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -2762,20 +2762,64 @@ expect 0.000898315284 0.000000000000 ------------------------------------------------------------------------------- operation +proj=isea +a=6400000 ------------------------------------------------------------------------------- -tolerance 0.1 mm +tolerance 0.2 mm accept 2 1 -expect -1097074.948022474 3442909.309037183 +expect -1097074.948153475765139 3442909.309747453313321 +roundtrip 1 + accept 2 -1 -expect -1097074.948264795 3233611.728585708 +expect -1097074.948149705072865 3233611.728292400948703 +roundtrip 1 + accept -2 1 -expect -1575486.353641554 3442168.342028188 +expect -1575486.353775786235929 3442168.342736063525081 +roundtrip 1 + accept -2 -1 -expect -1575486.353880283 3234352.695594706 +expect -1575486.353772019501776 3234352.695310209877789 +roundtrip 1 operation +proj=isea +mode=hex +resolution=31 accept 0 0 expect failure +------------------------------------------------------------------------------- +operation +proj=isea +R=6371007.18091875 +------------------------------------------------------------------------------- +tolerance 0.2 mm + +accept -168.75 58.282525588539 +expect -19186144.870842020958662 3323137.771944524254650 +roundtrip 1 + +accept 11.25 58.282525588539 +expect -15348915.896747918799520 9969413.315350906923413 +roundtrip 1 + +accept -110 54 +expect -15321401.505530973896384 3338358.859094056300819 +roundtrip 1 + +accept -75 45 +expect -12774358.709073608741164 4373188.646695702336729 +roundtrip 1 + +accept 2 49 +expect -642252.939347098814324 8796229.009143760427833 +roundtrip 1 + +accept 0 0 +expect -1331454.074623266700655 3323137.771634854841977 +roundtrip 1 + +accept 90 0 +expect 8564460.639100870117545 593869.297485541785136 +roundtrip 1 + +accept 0 45 +expect -837334.699958428042009 8323409.759132191538811 +roundtrip 1 + =============================================================================== # Kavrayskiy V # PCyl., Sph. @@ -5988,7 +6032,7 @@ accept 1 0.5 expect 45 0 accept 0 0 expect -45 -35.446011426401625 -accept 1 0 +accept 1 0 expect 45 -35.446011426401625 accept 1 1 expect 45 35.446011426401625 From 5f0aedf1e80a1a979649f9a24b809773513b9b36 Mon Sep 17 00:00:00 2001 From: Jerome St-Louis Date: Tue, 13 Aug 2024 21:50:29 -0400 Subject: [PATCH 2/4] ISEA: Optimizations, unification and simplification - Eliminating several redundant trigonometric function calls - Single shared definition of dodecahedron vertices (centers of icosahedral triangular faces) for forward and inverse - Shared data structure for forward and inverse - Consistent use of 0-base index for triangular faces --- src/projections/isea.cpp | 438 +++++++++++++++++---------------------- 1 file changed, 189 insertions(+), 249 deletions(-) diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp index 6f8ba66e49..1145a4ac49 100644 --- a/src/projections/isea.cpp +++ b/src/projections/isea.cpp @@ -119,6 +119,7 @@ #define SQRT3 1.73205080756887729352744634150587236694280525381038 #define sin60 (SQRT3 / 2.0) +#define cos30 (SQRT3 / 2.0) // tang * sin(60 deg) #define TABLE_G (tang * sin60) @@ -134,6 +135,11 @@ #define ISEA_STD_LONG .19634954084936207740 namespace { // anonymous namespace + +struct GeoPoint { + double lat, lon; +}; // In radians + struct hex { int iso; long x, y, z; @@ -237,103 +243,48 @@ enum isea_address_form { ISEA_HEX }; -struct isea_dgg { - int polyhedron; /* ignored, icosahedron */ - double o_lat, o_lon, o_az; /* orientation, radians */ - int topology; /* ignored, hexagon */ - int aperture; /* valid values depend on partitioning method */ - int resolution; - double radius; /* radius of the earth in meters, ignored 1.0 */ - int output; /* an isea_address_form */ - int triangle; /* triangle of last transformed point */ - int quad; /* quad of last transformed point */ - unsigned long serial; +struct isea_sincos { + double s, c; }; struct isea_pt { double x, y; }; -struct isea_geo { - double longitude, lat; -}; - } // anonymous namespace -/* No longer used - -enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = numIcosahedronFaces }; -enum isea_topology { ISEA_HEXAGON = 6, ISEA_TRIANGLE = 3, ISEA_DIAMOND = 4 }; - -enum snyder_polyhedron { - SNYDER_POLY_HEXAGON, - SNYDER_POLY_PENTAGON, - SNYDER_POLY_TETRAHEDRON, - SNYDER_POLY_CUBE, - SNYDER_POLY_OCTAHEDRON, - SNYDER_POLY_DODECAHEDRON, - SNYDER_POLY_ICOSAHEDRON -}; - -struct snyder_constants { - double g, G, theta; - // cppcheck-suppress unusedStructMember - double ea_w, ea_a, ea_b, g_w, g_a, g_b; -}; - -// TODO put these in radians to avoid a later conversion -static const struct snyder_constants constants[] = { - {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0}, - {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, - {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0}, -}; -*/ - -static struct isea_geo vertex[] = { - {0.0, DEG90}, {DEG180, V_LAT}, {-DEG108, V_LAT}, {-DEG36, V_LAT}, - {DEG36, V_LAT}, {DEG108, V_LAT}, {-DEG144, -V_LAT}, {-DEG72, -V_LAT}, - {0.0, -V_LAT}, {DEG72, -V_LAT}, {DEG144, -V_LAT}, {0.0, -DEG90}}; - -/* TODO make an isea_pt array of the vertices as well */ - -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}; - -/* triangle Centers */ -static const struct isea_geo icostriangles[] = { - {0.0, 0.0}, {-DEG144, E_RAD}, {-DEG72, E_RAD}, {0.0, E_RAD}, - {DEG72, E_RAD}, {DEG144, E_RAD}, {-DEG144, F_RAD}, {-DEG72, F_RAD}, - {0.0, F_RAD}, {DEG72, F_RAD}, {DEG144, F_RAD}, {-DEG108, -F_RAD}, - {-DEG36, -F_RAD}, {DEG36, -F_RAD}, {DEG108, -F_RAD}, {DEG180, -F_RAD}, - {-DEG108, -E_RAD}, {-DEG36, -E_RAD}, {DEG36, -E_RAD}, {DEG108, -E_RAD}, - {DEG180, -E_RAD}, -}; - -static double az_adjustment(int triangle) { - double adj; - - struct isea_geo v; - struct isea_geo c; +// distortion +// static double maximumAngularDistortion = 17.27; +// static double maximumScaleVariation = 1.163; +// static double minimumScaleVariation = .860; - v = vertex[tri_v1[triangle]]; - c = icostriangles[triangle]; +// Vertices of dodecahedron centered in icosahedron triangular faces +static const GeoPoint facesCenterDodecahedronVertices[numIcosahedronFaces] = { + {E_RAD, DEG_TO_RAD * -144}, {E_RAD, DEG_TO_RAD * -72}, + {E_RAD, DEG_TO_RAD * 0}, {E_RAD, DEG_TO_RAD * 72}, + {E_RAD, DEG_TO_RAD * 144}, {F_RAD, DEG_TO_RAD * -144}, + {F_RAD, DEG_TO_RAD * -72}, {F_RAD, DEG_TO_RAD * 0}, + {F_RAD, DEG_TO_RAD * 72}, {F_RAD, DEG_TO_RAD * 144}, + {-F_RAD, DEG_TO_RAD * -108}, {-F_RAD, DEG_TO_RAD * -36}, + {-F_RAD, DEG_TO_RAD * 36}, {-F_RAD, DEG_TO_RAD * 108}, + {-F_RAD, DEG_TO_RAD * 180}, {-E_RAD, DEG_TO_RAD * -108}, + {-E_RAD, DEG_TO_RAD * -36}, {-E_RAD, DEG_TO_RAD * 36}, + {-E_RAD, DEG_TO_RAD * 108}, {-E_RAD, DEG_TO_RAD * 180}}; - /* TODO looks like the adjustment is always either 0 or 180 */ - /* at least if you pick your vertex carefully */ - adj = atan2(cos(v.lat) * sin(v.longitude - c.longitude), - cos(c.lat) * sin(v.lat) - - sin(c.lat) * cos(v.lat) * cos(v.longitude - c.longitude)); - return adj; +// NOTE: Very similar to ISEAPlanarProjection::faceOrientation(), +// but the forward projection sometimes is returning a negative M_PI +static inline double az_adjustment(int triangle) { + if ((triangle >= 5 && triangle <= 9) || triangle == 15 || triangle == 16) + return M_PI; + else if (triangle >= 17) + return -M_PI; + return 0; } static struct isea_pt isea_triangle_xy(int triangle) { struct isea_pt c; - triangle = (triangle - 1) % numIcosahedronFaces; + triangle %= numIcosahedronFaces; c.x = TABLE_G * ((triangle % 5) - 2) * 2.0; if (triangle > 9) { @@ -362,16 +313,39 @@ static struct isea_pt isea_triangle_xy(int triangle) { return c; } -/* snyder eq 14 */ -static double sph_azimuth(double f_lon, double f_lat, double t_lon, - double t_lat) { - double az; +namespace { // anonymous namespace - az = atan2(cos(t_lat) * sin(t_lon - f_lon), - cos(f_lat) * sin(t_lat) - - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)); - return az; -} +class ISEAPlanarProjection; + +struct pj_isea_data { + int polyhedron; /* ignored, icosahedron */ + double o_lat, o_lon, o_az; /* orientation, radians */ + int topology; /* ignored, hexagon */ + int aperture; /* valid values depend on partitioning method */ + int resolution; + double radius; /* radius of the earth in meters, ignored 1.0 */ + int output; /* an isea_address_form */ + int triangle; /* triangle of last transformed point */ + int quad; /* quad of last transformed point */ + isea_sincos vertexLatSinCos[numIcosahedronFaces]; + unsigned long serial; + + double R2; + double Rprime; + double Rprime2X; + double RprimeTang; + double Rprime2Tan2g; + double triTang; + double centerToBase; + double triWidth; + double yOffsets[4]; + double xo, yo; + double sx, sy; + ISEAPlanarProjection *p; + + void initialize(const PJ *P); +}; +} // anonymous namespace #ifdef _MSC_VER #pragma warning(push) @@ -380,27 +354,32 @@ static double sph_azimuth(double f_lon, double f_lat, double t_lon, #endif /* coord needs to be in radians */ -static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { +static int isea_snyder_forward(const struct pj_isea_data *data, + const struct GeoPoint *ll, struct isea_pt *out) { int i; - - /* additional variables from snyder */ - double q, H, Ag, Azprime, Az, dprime, f, rho, x, y; - - /* variables used to store intermediate results */ - double az_offset; - - /* how many multiples of 60 degrees we adjust the azimuth */ - int Az_adjust_multiples; + double sinLat = sin(ll->lat), cosLat = cos(ll->lat); /* * TODO by locality of reference, start by trying the same triangle * as last time */ - for (i = 1; i <= numIcosahedronFaces; i++) { - struct isea_geo center = icostriangles[i]; - double cosZ = sin(center.lat) * sin(ll->lat) + - cos(center.lat) * cos(ll->lat) * - cos(ll->longitude - center.longitude); + for (i = 0; i < numIcosahedronFaces; i++) { + /* additional variables from snyder */ + double q, H, Ag, Azprime, Az, dprime, f, rho, x, y; + + /* variables used to store intermediate results */ + double az_offset; + + /* how many multiples of 60 degrees we adjust the azimuth */ + int Az_adjust_multiples; + + const struct GeoPoint *center = &facesCenterDodecahedronVertices[i]; + const struct isea_sincos *centerLatSinCos = &data->vertexLatSinCos[i]; + double dLon = ll->lon - center->lon; + double cosLat_cosLon = cosLat * cos(dLon); + double cosZ = + centerLatSinCos->s * sinLat + centerLatSinCos->c * cosLat_cosLon; + double sinAz, cosAz; /* step 1 */ double z = fabs(cosZ) < 1E-15 ? 0 : acos(cosZ); @@ -409,7 +388,9 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { continue; } - Az = sph_azimuth(center.longitude, center.lat, ll->longitude, ll->lat); + /* snyder eq 14 */ + Az = atan2(cosLat * sin(dLon), centerLatSinCos->c * sinLat - + centerLatSinCos->s * cosLat_cosLon); /* step 2 */ @@ -445,7 +426,9 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { /* step 3 */ /* Calculate q from eq 9. */ - q = atan2(tang, cos(Az) + sin(Az) * cotTheta); + cosAz = cos(Az); + sinAz = sin(Az); + q = atan2(tang, cosAz + sinAz * cotTheta); /* not in this triangle */ if (z > q + 0.000005) { @@ -459,8 +442,7 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { /* R' in the paper is for the truncated (icosahedron?) */ /* eq 6 */ - H = acos(sin(Az) * sinGcosSDC2VoS /* sin(G) * cos(g) */ - - cos(Az) * cosG); + H = acos(sinAz * sinGcosSDC2VoS /* sin(G) * cos(g) */ - cosAz * cosG); /* eq 7 */ /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */ @@ -510,7 +492,7 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { */ fprintf(stderr, "impossible transform: %f %f is not on any triangle\n", - PJ_TODEG(ll->longitude), PJ_TODEG(ll->lat)); + PJ_TODEG(ll->lon), PJ_TODEG(ll->lat)); exit(EXIT_FAILURE); } @@ -534,11 +516,11 @@ static int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out) { * * TODO take a result pointer */ -static struct isea_geo snyder_ctran(const struct isea_geo &np, - const struct isea_geo &pt) { - struct isea_geo result; - double phi = pt.lat, lambda = pt.longitude; - double alpha = np.lat, beta = np.longitude; +static struct GeoPoint snyder_ctran(const struct GeoPoint &np, + const struct GeoPoint &pt) { + struct GeoPoint result; + double phi = pt.lat, lambda = pt.lon; + double alpha = np.lat, beta = np.lon; double dlambda = lambda - beta /* lambda0 */; double cos_p = cos(phi), sin_p = sin(phi); double cos_a = cos(alpha), sin_a = sin(alpha); @@ -565,35 +547,37 @@ static struct isea_geo snyder_ctran(const struct isea_geo &np, result.lat = fabs(sin_phip - 1.0) < 1E-15 ? M_PI / 2 : fabs(sin_phip + 1.0) < 1E-15 ? -M_PI / 2 : asin(sin_phip); - result.longitude = lambdap; + result.lon = lambdap; return result; } -static struct isea_geo isea_ctran(const struct isea_geo *np, - const struct isea_geo *pt, double lon0) { - struct isea_geo cnp = {np->longitude + M_PI, np->lat}; - struct isea_geo npt = snyder_ctran(cnp, *pt); +static struct GeoPoint isea_ctran(const struct GeoPoint *np, + const struct GeoPoint *pt, double lon0) { + struct GeoPoint cnp = {np->lat, np->lon + M_PI}; + struct GeoPoint npt = snyder_ctran(cnp, *pt); - npt.longitude -= (/* M_PI */ -lon0 + np->longitude); + npt.lon -= (/* M_PI */ -lon0 + np->lon); /* * snyder is down tri 3, isea is along side of tri1 from vertex 0 to * vertex 1 these are 180 degrees apart */ - // npt.longitude += M_PI; + // npt.lon += M_PI; - /* normalize longitude */ - npt.longitude = fmod(npt.longitude, 2 * M_PI); - while (npt.longitude > M_PI) - npt.longitude -= 2 * M_PI; - while (npt.longitude < -M_PI) - npt.longitude += 2 * M_PI; + /* normalize lon */ + npt.lon = fmod(npt.lon, 2 * M_PI); + while (npt.lon > M_PI) + npt.lon -= 2 * M_PI; + while (npt.lon < -M_PI) + npt.lon += 2 * M_PI; return npt; } /* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */ -static int isea_grid_init(struct isea_dgg *g) { +static int isea_grid_init(struct pj_isea_data *g) { + int i; + if (!g) return 0; @@ -606,10 +590,15 @@ static int isea_grid_init(struct isea_dgg *g) { g->radius = 1.0; g->topology = 6; + for (i = 0; i < numIcosahedronFaces; i++) { + const GeoPoint *c = &facesCenterDodecahedronVertices[i]; + g->vertexLatSinCos[i].s = sin(c->lat); + g->vertexLatSinCos[i].c = cos(c->lat); + } return 1; } -static void isea_orient_isea(struct isea_dgg *g) { +static void isea_orient_isea(struct pj_isea_data *g) { if (!g) return; g->o_lat = ISEA_STD_LAT; @@ -617,7 +606,7 @@ static void isea_orient_isea(struct isea_dgg *g) { g->o_az = 0.0; } -static void isea_orient_pole(struct isea_dgg *g) { +static void isea_orient_pole(struct pj_isea_data *g) { if (!g) return; g->o_lat = M_PI / 2.0; @@ -625,17 +614,17 @@ static void isea_orient_pole(struct isea_dgg *g) { g->o_az = 0; } -static int isea_transform(struct isea_dgg *g, struct isea_geo *in, +static int isea_transform(struct pj_isea_data *g, struct GeoPoint *in, struct isea_pt *out) { - struct isea_geo i, pole; + struct GeoPoint i, pole; int tri; pole.lat = g->o_lat; - pole.longitude = g->o_lon; + pole.lon = g->o_lon; i = isea_ctran(&pole, in, g->o_az); - tri = isea_snyder_forward(&i, out); + tri = isea_snyder_forward(g, &i, out); out->x *= g->radius; out->y *= g->radius; g->triangle = tri; @@ -643,7 +632,7 @@ static int isea_transform(struct isea_dgg *g, struct isea_geo *in, return tri; } -#define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1) +#define DOWNTRI(tri) ((tri / 5) % 2 == 1) static void isea_rotate(struct isea_pt *pt, double degrees) { double rad; @@ -682,20 +671,20 @@ static int isea_tri_plane(int tri, struct isea_pt *pt, double radius) { static int isea_ptdd(int tri, struct isea_pt *pt) { int downtri, quadz; - downtri = (((tri - 1) / 5) % 2 == 1); - quadz = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1; + downtri = ((tri / 5) % 2 == 1); + quadz = (tri % 5) + (tri / 10) * 5 + 1; isea_rotate(pt, downtri ? 240.0 : 60.0); if (downtri) { pt->x += 0.5; /* pt->y += cos(30.0 * M_PI / 180.0); */ - pt->y += .86602540378443864672; + pt->y += cos30; } return quadz; } -static int isea_dddi_ap3odd(struct isea_dgg *g, int quadz, struct isea_pt *pt, - struct isea_pt *di) { +static int isea_dddi_ap3odd(struct pj_isea_data *g, int quadz, + struct isea_pt *pt, struct isea_pt *di) { struct isea_pt v; double hexwidth; double sidelength; /* in hexes */ @@ -771,7 +760,7 @@ static int isea_dddi_ap3odd(struct isea_dgg *g, int quadz, struct isea_pt *pt, return quadz; } -static int isea_dddi(struct isea_dgg *g, int quadz, struct isea_pt *pt, +static int isea_dddi(struct pj_isea_data *g, int quadz, struct isea_pt *pt, struct isea_pt *di) { struct isea_pt v; double hexwidth; @@ -850,7 +839,7 @@ static int isea_dddi(struct isea_dgg *g, int quadz, struct isea_pt *pt, return quadz; } -static int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt, +static int isea_ptdi(struct pj_isea_data *g, int tri, struct isea_pt *pt, struct isea_pt *di) { struct isea_pt v; int quadz; @@ -863,7 +852,7 @@ static int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt, /* q2di to seqnum */ -static long isea_disn(struct isea_dgg *g, int quadz, struct isea_pt *di) { +static long isea_disn(struct pj_isea_data *g, int quadz, struct isea_pt *di) { long sidelength; long sn, height; long hexes; @@ -900,7 +889,7 @@ static long isea_disn(struct isea_dgg *g, int quadz, struct isea_pt *di) { * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf */ /* convert a q2di to global hex coord */ -static int isea_hex(struct isea_dgg *g, int tri, struct isea_pt *pt, +static int isea_hex(struct pj_isea_data *g, int tri, struct isea_pt *pt, struct isea_pt *hex) { struct isea_pt v; #ifdef FIXME @@ -966,7 +955,8 @@ static int isea_hex(struct isea_dgg *g, int tri, struct isea_pt *pt, #endif } -static struct isea_pt isea_forward(struct isea_dgg *g, struct isea_geo *in) { +static struct isea_pt isea_forward(struct pj_isea_data *g, + struct GeoPoint *in) { int tri; struct isea_pt out, coord; @@ -1019,50 +1009,19 @@ static struct isea_pt isea_forward(struct isea_dgg *g, struct isea_geo *in) { PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph"; -namespace { // anonymous namespace - -struct GeoPoint { - double lat, lon; -}; // In radians - -class ISEAPlanarProjection; - -struct isea_inverse_params { - double R2; - double Rprime; - double Rprime2X; - double RprimeTang; - double Rprime2Tan2g; - double triTang; - double centerToBase; - double triWidth; - double yOffsets[4]; - double xo, yo; - double sx, sy; - ISEAPlanarProjection *p; - - void initialize(const PJ *P); -}; - -struct pj_isea_data { - struct isea_dgg dgg; - struct isea_inverse_params inv; -}; -} // anonymous namespace - static PJ_XY isea_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */ PJ_XY xy = {0.0, 0.0}; struct pj_isea_data *Q = static_cast(P->opaque); struct isea_pt out; - struct isea_geo in; + struct GeoPoint in; - in.longitude = lp.lam; + // TODO: Convert geodetic latitude to authalic latitude if not + // spherical as in eqearth, healpix, laea, etc. in.lat = lp.phi; + in.lon = lp.lam; try { - // TODO: Convert geodetic latitude to authalic latitude if not - // spherical as in eqearth, healpix, laea, etc. - out = isea_forward(&Q->dgg, &in); + out = isea_forward(Q, &in); } catch (const char *) { proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); return proj_coord_error().xy; @@ -1088,18 +1047,18 @@ PJ *PJ_PROJECTION(isea) { // https://brsr.github.io/2021/08/31/snyder-equal-area.html P->fwd = isea_s_forward; P->inv = isea_s_inverse; - isea_grid_init(&Q->dgg); - Q->dgg.output = ISEA_PLANE; + isea_grid_init(Q); + Q->output = ISEA_PLANE; - /* P->dgg.radius = P->a; / * otherwise defaults to 1 */ + /* P->radius = P->a; / * otherwise defaults to 1 */ /* calling library will scale, I think */ opt = pj_param(P->ctx, P->params, "sorient").s; if (opt) { if (!strcmp(opt, "isea")) { - isea_orient_isea(&Q->dgg); + isea_orient_isea(Q); } else if (!strcmp(opt, "pole")) { - isea_orient_pole(&Q->dgg); + isea_orient_pole(Q); } else { proj_log_error( P, @@ -1110,27 +1069,27 @@ PJ *PJ_PROJECTION(isea) { } if (pj_param(P->ctx, P->params, "tazi").i) { - Q->dgg.o_az = pj_param(P->ctx, P->params, "razi").f; + Q->o_az = pj_param(P->ctx, P->params, "razi").f; } if (pj_param(P->ctx, P->params, "tlon_0").i) { - Q->dgg.o_lon = pj_param(P->ctx, P->params, "rlon_0").f; + Q->o_lon = pj_param(P->ctx, P->params, "rlon_0").f; } if (pj_param(P->ctx, P->params, "tlat_0").i) { - Q->dgg.o_lat = pj_param(P->ctx, P->params, "rlat_0").f; + Q->o_lat = pj_param(P->ctx, P->params, "rlat_0").f; } opt = pj_param(P->ctx, P->params, "smode").s; if (opt) { if (!strcmp(opt, "plane")) { - Q->dgg.output = ISEA_PLANE; + Q->output = ISEA_PLANE; } else if (!strcmp(opt, "di")) { - Q->dgg.output = ISEA_Q2DI; + Q->output = ISEA_Q2DI; } else if (!strcmp(opt, "dd")) { - Q->dgg.output = ISEA_Q2DD; + Q->output = ISEA_Q2DD; } else if (!strcmp(opt, "hex")) { - Q->dgg.output = ISEA_HEX; + Q->output = ISEA_HEX; } else { proj_log_error(P, _("Invalid value for mode: only plane, di, dd or " "hex are supported")); @@ -1140,22 +1099,22 @@ PJ *PJ_PROJECTION(isea) { } if (pj_param(P->ctx, P->params, "trescale").i) { - Q->dgg.radius = ISEA_SCALE; + Q->radius = ISEA_SCALE; } if (pj_param(P->ctx, P->params, "tresolution").i) { - Q->dgg.resolution = pj_param(P->ctx, P->params, "iresolution").i; + Q->resolution = pj_param(P->ctx, P->params, "iresolution").i; } else { - Q->dgg.resolution = 4; + Q->resolution = 4; } if (pj_param(P->ctx, P->params, "taperture").i) { - Q->dgg.aperture = pj_param(P->ctx, P->params, "iaperture").i; + Q->aperture = pj_param(P->ctx, P->params, "iaperture").i; } else { - Q->dgg.aperture = 3; + Q->aperture = 3; } - Q->inv.initialize(P); + Q->initialize(P); return P; } @@ -1165,24 +1124,6 @@ PJ *PJ_PROJECTION(isea) { #define inf std::numeric_limits::infinity() -// distortion -// static double maximumAngularDistortion = 17.27; -// static double maximumScaleVariation = 1.163; -// static double minimumScaleVariation = .860; - -// Vertices of dodecahedron centered in icosahedron faces -static const GeoPoint facesCenterDodecahedronVertices[numIcosahedronFaces] = { - {E_RAD, DEG_TO_RAD * -144}, {E_RAD, DEG_TO_RAD * -72}, - {E_RAD, DEG_TO_RAD * 0}, {E_RAD, DEG_TO_RAD * 72}, - {E_RAD, DEG_TO_RAD * 144}, {F_RAD, DEG_TO_RAD * -144}, - {F_RAD, DEG_TO_RAD * -72}, {F_RAD, DEG_TO_RAD * 0}, - {F_RAD, DEG_TO_RAD * 72}, {F_RAD, DEG_TO_RAD * 144}, - {-F_RAD, DEG_TO_RAD * -108}, {-F_RAD, DEG_TO_RAD * -36}, - {-F_RAD, DEG_TO_RAD * 36}, {-F_RAD, DEG_TO_RAD * 108}, - {-F_RAD, DEG_TO_RAD * 180}, {-E_RAD, DEG_TO_RAD * -108}, - {-E_RAD, DEG_TO_RAD * -36}, {-E_RAD, DEG_TO_RAD * 36}, - {-E_RAD, DEG_TO_RAD * 108}, {-E_RAD, DEG_TO_RAD * 180}}; - // static define precision = DEG_TO_RAD * 1e-9; #define precision (DEG_TO_RAD * 1e-11) #define precisionPerDefinition (DEG_TO_RAD * 1e-5) @@ -1204,27 +1145,25 @@ class ISEAPlanarProjection { : orientation(value), cosOrientationLat(cos(value.lat)), sinOrientationLat(sin(value.lat)) {} - bool cartesianToGeo(const PJ_XY &inPosition, - const isea_inverse_params ¶ms, GeoPoint &result) { + bool cartesianToGeo(const PJ_XY &inPosition, const pj_isea_data *params, + GeoPoint &result) { bool r = false; - static const double epsilon = - 2E-8; // NOTE: 1E-11 seems too small for forward projection - // precision at boundaries + static const double epsilon = 1E-11; int face = 0; PJ_XY position = inPosition; #define sr -sin60 // sin(-60) #define cr 0.5 // cos(-60) if (position.x < 0 || - (position.x < params.triWidth / 2 && position.y < 0 && + (position.x < params->triWidth / 2 && position.y < 0 && position.y * cr < position.x * sr)) - position.x += 5 * params.triWidth; // Wrap around + position.x += 5 * params->triWidth; // Wrap around // Rotate and shear to determine face if not stored in position.z #define shearX (1.0 / SQRT3) double yp = -(position.x * sr + position.y * cr); double x = - (position.x * cr - position.y * sr + yp * shearX) * params.sx; - double y = yp * params.sy; + (position.x * cr - position.y * sr + yp * shearX) * params->sx; + double y = yp * params->sy; #undef shearX #undef sr #undef cr @@ -1288,9 +1227,9 @@ class ISEAPlanarProjection { if (face) { int fy = (face - 1) / 5, fx = (face - 1) - 5 * fy; double rx = - position.x - (2 * fx + fy / 2 + 1) * params.triWidth / 2; + position.x - (2 * fx + fy / 2 + 1) * params->triWidth / 2; double ry = - position.y - (params.yOffsets[fy] + 3 * params.centerToBase); + position.y - (params->yOffsets[fy] + 3 * params->centerToBase); GeoPoint dst; r = icosahedronToSphere({face - 1, rx, ry}, params, dst); @@ -1307,8 +1246,8 @@ class ISEAPlanarProjection { // Converts coordinates on the icosahedron to geographic coordinates // (inverse projection) - bool icosahedronToSphere(const ISEAFacePoint &c, - const isea_inverse_params ¶ms, GeoPoint &r) { + bool icosahedronToSphere(const ISEAFacePoint &c, const pj_isea_data *params, + GeoPoint &r) { if (c.face >= 0 && c.face < numIcosahedronFaces) { double Az = atan2(c.x, c.y); // Az' double rho = sqrt(c.x * c.x + c.y * c.y); @@ -1326,11 +1265,11 @@ class ISEAPlanarProjection { { double sinAz = sin(Az), cosAz = cos(Az); double cotAz = cosAz / sinAz; - double area = params.Rprime2Tan2g / + double area = params->Rprime2Tan2g / (2 * (cotAz + cotTheta)); // A_G or A_{ABD} double deltaAz = 10 * precision; double degAreaOverR2Plus180Minus36 = - area / params.R2 - westVertexLon; + area / params->R2 - westVertexLon; double Az_earth = Az; while (fabs(deltaAz) > precision) { @@ -1353,15 +1292,15 @@ class ISEAPlanarProjection { double q = atan2(tang, (cosAz_earth + sinAz_earth * cotTheta)); double d = - params.RprimeTang / (cosAz + sinAz * cotTheta); // d' - double f = d / (params.Rprime2X * sin(q / 2)); // f - double z = 2 * asin(rho / (params.Rprime2X * f)); + params->RprimeTang / (cosAz + sinAz * cotTheta); // d' + double f = d / (params->Rprime2X * sin(q / 2)); // f + double z = 2 * asin(rho / (params->Rprime2X * f)); Az_earth -= AzAdjustment; { - double lat0 = - facesCenterDodecahedronVertices[c.face].lat; - double sinLat0 = sin(lat0), cosLat0 = cos(lat0); + const isea_sincos *latSinCos = + ¶ms->vertexLatSinCos[c.face]; + double sinLat0 = latSinCos->s, cosLat0 = latSinCos->c; double sinZ = sin(z), cosZ = cos(z); double cosLat0SinZ = cosLat0 * sinZ; double latSin = @@ -1433,15 +1372,15 @@ static ISEAPlanarProjection standardISEA( */ static ISEAPlanarProjection polarISEA({0, 0}); -void isea_inverse_params::initialize(const PJ *P) { +void pj_isea_data::initialize(const PJ *P) { struct pj_isea_data *Q = static_cast(P->opaque); // Only supporting default planar options for now - if (Q->dgg.output == ISEA_PLANE && Q->dgg.o_az == 0.0 && - Q->dgg.aperture == 3.0 && Q->dgg.resolution == 4.) { + if (Q->output == ISEA_PLANE && Q->o_az == 0.0 && Q->aperture == 3.0 && + Q->resolution == 4.) { // Only supporting +orient=isea and +orient=pole for now - if (Q->dgg.o_lat == ISEA_STD_LAT && Q->dgg.o_lon == ISEA_STD_LONG) + if (Q->o_lat == ISEA_STD_LAT && Q->o_lon == ISEA_STD_LONG) p = &standardISEA; - else if (Q->dgg.o_lat == M_PI / 2.0 && Q->dgg.o_lon == 0) + else if (Q->o_lat == M_PI / 2.0 && Q->o_lon == 0) p = &polarISEA; else p = nullptr; @@ -1480,16 +1419,17 @@ void isea_inverse_params::initialize(const PJ *P) { } // anonymous namespace static PJ_LP isea_s_inverse(PJ_XY xy, PJ *P) { - struct pj_isea_data *Q = static_cast(P->opaque); - ISEAPlanarProjection *p = Q->inv.p; + const struct pj_isea_data *Q = + static_cast(P->opaque); + ISEAPlanarProjection *p = Q->p; if (p) { // Default origin of +proj=isea is different (OGC:1534 is // +x_0=19186144.870934911 +y_0=-3323137.7717836285) - PJ_XY input{xy.x * P->a + Q->inv.xo, xy.y * P->a + Q->inv.yo}; + PJ_XY input{xy.x * P->a + Q->xo, xy.y * P->a + Q->yo}; GeoPoint result; - if (p->cartesianToGeo(input, Q->inv, result)) + if (p->cartesianToGeo(input, Q, result)) // TODO: Convert authalic latitude to geodetic latitude if not // spherical as in eqearth, healpix, laea, etc. return {result.lon, result.lat}; From 7874a3aa1b2368d58e4e38db2382649d6cc2e6c3 Mon Sep 17 00:00:00 2001 From: Jerome St-Louis Date: Tue, 13 Aug 2024 23:42:41 -0400 Subject: [PATCH 3/4] ISEA: Removing extraneous multiplications, unused code - NOTE: There may have been an undocumented +rescale= parameter which will no longer work --- src/projections/isea.cpp | 137 ++++++++++----------------------------- 1 file changed, 35 insertions(+), 102 deletions(-) diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp index 1145a4ac49..b55dea28dc 100644 --- a/src/projections/isea.cpp +++ b/src/projections/isea.cpp @@ -231,17 +231,7 @@ static void hexbin2(double width, double x, double y, long *i, long *j) { #define numIcosahedronFaces 20 namespace { // anonymous namespace -enum isea_address_form { - ISEA_GEO, - ISEA_Q2DI, - ISEA_SEQNUM, - ISEA_INTERLEAVE, - ISEA_PLANE, - ISEA_Q2DD, - ISEA_PROJTRI, - ISEA_VERTEX2DD, - ISEA_HEX -}; +enum isea_address_form { ISEA_PLANE, ISEA_Q2DI, ISEA_Q2DD, ISEA_HEX }; struct isea_sincos { double s, c; @@ -318,17 +308,13 @@ namespace { // anonymous namespace class ISEAPlanarProjection; struct pj_isea_data { - int polyhedron; /* ignored, icosahedron */ double o_lat, o_lon, o_az; /* orientation, radians */ - int topology; /* ignored, hexagon */ int aperture; /* valid values depend on partitioning method */ int resolution; - double radius; /* radius of the earth in meters, ignored 1.0 */ - int output; /* an isea_address_form */ - int triangle; /* triangle of last transformed point */ - int quad; /* quad of last transformed point */ + isea_address_form output; /* an isea_address_form */ + int triangle; /* triangle of last transformed point */ + int quad; /* quad of last transformed point */ isea_sincos vertexLatSinCos[numIcosahedronFaces]; - unsigned long serial; double R2; double Rprime; @@ -581,14 +567,11 @@ static int isea_grid_init(struct pj_isea_data *g) { if (!g) return 0; - g->polyhedron = numIcosahedronFaces; g->o_lat = ISEA_STD_LAT; g->o_lon = ISEA_STD_LONG; g->o_az = 0.0; g->aperture = 4; g->resolution = 6; - g->radius = 1.0; - g->topology = 6; for (i = 0; i < numIcosahedronFaces; i++) { const GeoPoint *c = &facesCenterDodecahedronVertices[i]; @@ -625,8 +608,6 @@ static int isea_transform(struct pj_isea_data *g, struct GeoPoint *in, i = isea_ctran(&pole, in, g->o_az); tri = isea_snyder_forward(g, &i, out); - out->x *= g->radius; - out->y *= g->radius; g->triangle = tri; return tri; @@ -652,15 +633,13 @@ static void isea_rotate(struct isea_pt *pt, double degrees) { pt->y = y; } -static int isea_tri_plane(int tri, struct isea_pt *pt, double radius) { +static int isea_tri_plane(int tri, struct isea_pt *pt) { struct isea_pt tc; /* center of triangle */ if (DOWNTRI(tri)) { isea_rotate(pt, 180.0); } tc = isea_triangle_xy(tri); - tc.x *= radius; - tc.y *= radius; pt->x += tc.x; pt->y += tc.y; @@ -850,40 +829,6 @@ static int isea_ptdi(struct pj_isea_data *g, int tri, struct isea_pt *pt, return quadz; } -/* q2di to seqnum */ - -static long isea_disn(struct pj_isea_data *g, int quadz, struct isea_pt *di) { - long sidelength; - long sn, height; - long hexes; - - if (quadz == 0) { - g->serial = 1; - return g->serial; - } - /* hexes in a quad */ - hexes = lround(pow(static_cast(g->aperture), - static_cast(g->resolution))); - if (quadz == 11) { - g->serial = 1 + 10 * hexes + 1; - return g->serial; - } - if (g->aperture == 3 && g->resolution % 2 == 1) { - height = lround(floor((pow(g->aperture, (g->resolution - 1) / 2.0)))); - sn = ((long)di->x) * height; - sn += ((long)di->y) / height; - sn += (quadz - 1) * hexes; - sn += 2; - } else { - sidelength = lround((pow(g->aperture, g->resolution / 2.0))); - sn = lround( - floor(((quadz - 1) * hexes + sidelength * di->x + di->y + 2))); - } - - g->serial = sn; - return sn; -} - /* TODO just encode the quad in the d or i coordinate * quad is 0-11, which can be four bits. * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf @@ -957,49 +902,35 @@ static int isea_hex(struct pj_isea_data *g, int tri, struct isea_pt *pt, static struct isea_pt isea_forward(struct pj_isea_data *g, struct GeoPoint *in) { - int tri; - struct isea_pt out, coord; - - tri = isea_transform(g, in, &out); - - if (g->output == ISEA_PLANE) { - isea_tri_plane(tri, &out, g->radius); - return out; - } - - /* convert to isea standard triangle size */ - out.x = out.x / g->radius * ISEA_SCALE; - out.y = out.y / g->radius * ISEA_SCALE; - out.x += 0.5; - out.y += 2.0 * .14433756729740644112; - - switch (g->output) { - case ISEA_PROJTRI: - /* nothing to do, already in projected triangle */ - break; - case ISEA_VERTEX2DD: - g->quad = isea_ptdd(tri, &out); - break; - case ISEA_Q2DD: - /* Same as above, we just don't print as much */ - g->quad = isea_ptdd(tri, &out); - break; - case ISEA_Q2DI: - g->quad = isea_ptdi(g, tri, &out, &coord); - return coord; - break; - case ISEA_SEQNUM: - isea_ptdi(g, tri, &out, &coord); - /* disn will set g->serial */ - isea_disn(g, g->quad, &coord); - return coord; - break; - case ISEA_HEX: - isea_hex(g, tri, &out, &coord); - return coord; - break; + isea_pt out; + int tri = isea_transform(g, in, &out); + + if (g->output == ISEA_PLANE) + isea_tri_plane(tri, &out); + else { + isea_pt coord; + + /* convert to isea standard triangle size */ + out.x *= ISEA_SCALE; // / g->radius; + out.y *= ISEA_SCALE; // / g->radius; + out.x += 0.5; + out.y += 2.0 * .14433756729740644112; + + switch (g->output) { + case ISEA_PLANE: + /* already handled above -- GCC should not be complaining */ + case ISEA_Q2DD: + /* Same as above, we just don't print as much */ + g->quad = isea_ptdd(tri, &out); + break; + case ISEA_Q2DI: + g->quad = isea_ptdi(g, tri, &out, &coord); + return coord; + case ISEA_HEX: + isea_hex(g, tri, &out, &coord); + return coord; + } } - return out; } @@ -1098,9 +1029,11 @@ PJ *PJ_PROJECTION(isea) { } } + /* REVIEW: Was this an undocumented +rescale= parameter? if (pj_param(P->ctx, P->params, "trescale").i) { Q->radius = ISEA_SCALE; } + */ if (pj_param(P->ctx, P->params, "tresolution").i) { Q->resolution = pj_param(P->ctx, P->params, "iresolution").i; From d1d8a5a75e33f03d40ccf4dbe0cba432695ddd91 Mon Sep 17 00:00:00 2001 From: Jerome St-Louis Date: Wed, 14 Aug 2024 05:05:06 -0400 Subject: [PATCH 4/4] ISEA: Save four trigonometric calls - Avoid calling isea_rotate() in isea_tri_plane() to rotate coordinates in bottom triangles by 180 degrees. --- src/projections/isea.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp index b55dea28dc..fb5bb696c5 100644 --- a/src/projections/isea.cpp +++ b/src/projections/isea.cpp @@ -280,6 +280,9 @@ static struct isea_pt isea_triangle_xy(int triangle) { if (triangle > 9) { c.x += TABLE_G; } + + // REVIEW: This is likely related to + // pj_isea_data::yOffsets switch (triangle / 5) { case 0: c.y = 5.0 * TABLE_H; @@ -633,17 +636,16 @@ static void isea_rotate(struct isea_pt *pt, double degrees) { pt->y = y; } -static int isea_tri_plane(int tri, struct isea_pt *pt) { +static void isea_tri_plane(int tri, struct isea_pt *pt) { struct isea_pt tc; /* center of triangle */ if (DOWNTRI(tri)) { - isea_rotate(pt, 180.0); + pt->x *= -1; + pt->y *= -1; } tc = isea_triangle_xy(tri); pt->x += tc.x; pt->y += tc.y; - - return tri; } /* convert projected triangle coords to quad xy coords, return quad number */ @@ -653,6 +655,8 @@ static int isea_ptdd(int tri, struct isea_pt *pt) { downtri = ((tri / 5) % 2 == 1); quadz = (tri % 5) + (tri / 10) * 5 + 1; + // NOTE: This would always be a 60 degrees rotation if the flip were + // already done as in isea_tri_plane() isea_rotate(pt, downtri ? 240.0 : 60.0); if (downtri) { pt->x += 0.5;