Merge pull request #4225 from ecere/master

ISEA: Fix bug introduced in 825682a
OSGeo · Aug 14, 2024 · 892a313 · 892a313
2 parents ae23ad5 + 724878b
commit 892a313
Show file tree

Hide file tree

Showing 2 changed files with 59 additions and 8 deletions.
diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp
@@ -342,6 +342,24 @@ struct pj_isea_data {
 #pragma warning(disable : 4702)
 #endif
 
+#define SAFE_ARC_EPSILON 1E-15
+
+static inline double safeArcSin(double t) {
+    return fabs(t) < SAFE_ARC_EPSILON         ? 0
+           : fabs(t - 1.0) < SAFE_ARC_EPSILON ? M_PI / 2
+           : fabs(t + 1.0) < SAFE_ARC_EPSILON ? -M_PI / 2
+                                              : asin(t);
+}
+
+static inline double safeArcCos(double t) {
+    return fabs(t) < SAFE_ARC_EPSILON       ? M_PI / 2
+           : fabs(t + 1) < SAFE_ARC_EPSILON ? M_PI
+           : fabs(t - 1) < SAFE_ARC_EPSILON ? 0
+                                            : acos(t);
+}
+
+#undef SAFE_ARC_EPSILON
+
 /* coord needs to be in radians */
 static int isea_snyder_forward(const struct pj_isea_data *data,
                                const struct GeoPoint *ll, struct isea_pt *out) {
@@ -371,7 +389,8 @@ static int isea_snyder_forward(const struct pj_isea_data *data,
         double sinAz, cosAz;
 
         /* step 1 */
-        double z = fabs(cosZ) < 1E-15 ? 0 : acos(cosZ);
+        double z = safeArcCos(cosZ);
+
         /* not on this triangle */
         if (z > sdc2vos /*g*/ + 0.000005) { /* TODO DBL_EPSILON */
             continue;
@@ -533,9 +552,7 @@ static struct GeoPoint snyder_ctran(const struct GeoPoint &np,
     while (lambdap < -M_PI)
         lambdap += 2 * M_PI;
 
-    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.lat = safeArcSin(sin_phip);
     result.lon = lambdap;
     return result;
 }
@@ -1242,10 +1259,7 @@ class ISEAPlanarProjection {
                         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 lat = safeArcSin(latSin);
                         double lon =
                             facesCenterDodecahedronVertices[c.face].lon +
                             atan2(sin(Az_earth) * cosLat0SinZ,

diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
@@ -2820,6 +2820,43 @@ accept 0 45
 expect -837334.699958428042009 8323409.759132191538811
 roundtrip 1
 
+-------------------------------------------------------------------------------
+operation +proj=isea +R=6371007.18091875 +orient=pole
+-------------------------------------------------------------------------------
+tolerance 0.2 mm
+
+accept -168.75 58.282525588539
+expect -16702163.549901897087693 6386395.630649688653648
+roundtrip 1
+
+accept 11.25 58.282525588539
+expect 619648.646531744743697 6212947.536539182066917
+roundtrip 1
+
+accept -110 54
+expect -13285649.857057726010680 6149501.348902118392289
+roundtrip 1
+
+accept -75 45
+expect -7921366.529368571005762 4728387.055336073972285
+roundtrip 1
+
+accept 2 49
+expect 152616.434999307675753 5152048.791301283054054
+roundtrip 1
+
+accept 0 0
+expect 0 -195097.133640714135254
+roundtrip 1
+
+accept 90 0
+expect 9593072.435467451811 0
+roundtrip 1
+
+accept 0 45
+expect 0 4726854.770339427515864
+roundtrip 1
+
 ===============================================================================
 # Kavrayskiy V
 # 	PCyl., Sph.