Improve H3 Vec2d#v2dIntersect method (#92433)

libs/h3/src/main/java/org/elasticsearch/h3/FaceIJK.java

@@ -628,7 +628,7 @@ final class FaceIJK {
                 adjacent hexagon edge will lie completely on a single icosahedron
                 face, and no additional vertex is required.
                 */
-                final boolean isIntersectionAtVertex = orig2d0.equals(inter) || orig2d1.equals(inter);
+                final boolean isIntersectionAtVertex = orig2d0.numericallyIdentical(inter) || orig2d1.numericallyIdentical(inter);
                 if (isIntersectionAtVertex == false) {
                     final LatLng point = inter.hex2dToGeo(this.face, adjRes, true);
                     boundary.add(point);

libs/h3/src/main/java/org/elasticsearch/h3/Vec2d.java

@@ -32,6 +32,8 @@ final class Vec2d {
     /** sin(60') */
     private static final double M_SIN60 = Constants.M_SQRT3_2;
 
+    private static final double VEC2D_RESOLUTION = 1e-7;
+
     /**
      * icosahedron face centers in lat/lng radians
      */
@@ -257,6 +259,10 @@ final class Vec2d {
         return coordIJK;
     }
 
+    public boolean numericallyIdentical(Vec2d vec2d) {
+        return Math.abs(vec2d.x - x) < VEC2D_RESOLUTION && Math.abs(vec2d.y - y) < VEC2D_RESOLUTION;
+    }
+
     @Override
     public boolean equals(Object o) {
         if (this == o) return true;
@@ -280,16 +286,14 @@ final class Vec2d {
      * @param p3 The second endpoint of the second line.
      */
     public static Vec2d v2dIntersect(Vec2d p0, Vec2d p1, Vec2d p2, Vec2d p3) {
-        double[] s1 = new double[2], s2 = new double[2];
-        s1[0] = p1.x - p0.x;
-        s1[1] = p1.y - p0.y;
-        s2[0] = p3.x - p2.x;
-        s2[1] = p3.y - p2.y;
+        final double s1x = p1.x - p0.x;
+        final double s1y = p1.y - p0.y;
+        final double s2x = p3.x - p2.x;
+        final double s2y = p3.y - p2.y;
 
-        float t;
-        t = (float) ((s2[0] * (p0.y - p2.y) - s2[1] * (p0.x - p2.x)) / (-s2[0] * s1[1] + s1[0] * s2[1]));
+        final double t = ((s2x * (p0.y - p2.y) - s2y * (p0.x - p2.x)) / (-s2x * s1y + s1x * s2y));
 
-        return new Vec2d(p0.x + (t * s1[0]), p0.y + (t * s1[1]));
+        return new Vec2d(p0.x + (t * s1x), p0.y + (t * s1y));
     }
 
     /**

libs/h3/src/test/java/org/elasticsearch/h3/CellBoundaryTests.java

@@ -18,6 +18,8 @@
  */
 package org.elasticsearch.h3;
 
+import org.apache.lucene.geo.GeoEncodingUtils;
+import org.apache.lucene.tests.geo.GeoTestUtil;
 import org.elasticsearch.test.ESTestCase;
 
 import java.io.BufferedReader;
@@ -30,6 +32,9 @@ import java.util.List;
 import java.util.StringTokenizer;
 import java.util.zip.GZIPInputStream;
 
+import static org.hamcrest.Matchers.either;
+import static org.hamcrest.Matchers.equalTo;
+
 public class CellBoundaryTests extends ESTestCase {
 
     public void testRes0() throws Exception {
@@ -171,8 +176,46 @@ public class CellBoundaryTests extends ESTestCase {
         CellBoundary boundary = H3.h3ToGeoBoundary(h3Address);
         assert boundary.numPoints() == points.size();
         for (int i = 0; i < boundary.numPoints(); i++) {
-            assertEquals(h3Address, points.get(i)[0], boundary.getLatLon(i).getLatDeg(), 1e-8);
-            assertEquals(h3Address, points.get(i)[1], boundary.getLatLon(i).getLonDeg(), 1e-8);
+            assertEquals(h3Address, points.get(i)[0], boundary.getLatLon(i).getLatDeg(), 5e-7);
+            assertEquals(h3Address, points.get(i)[1], boundary.getLatLon(i).getLonDeg(), 5e-7);
+        }
+    }
+
+    public void testNumericEquivalentSharedBoundary() {
+        // we consider boundaries numerical equivalent if after encoded them using lucene, they resolve to the same number.
+        long h3 = H3.geoToH3(GeoTestUtil.nextLatitude(), GeoTestUtil.nextLongitude(), randomIntBetween(0, 15));
+        CellBoundary boundary = H3.h3ToGeoBoundary(h3);
+        for (long r : H3.hexRing(h3)) {
+            int count = 0;
+            CellBoundary ringBoundary = H3.h3ToGeoBoundary(r);
+            for (int i = 0; i < boundary.numPoints(); i++) {
+                LatLng latLng1 = boundary.getLatLon(i % boundary.numPoints());
+                LatLng latLng2 = boundary.getLatLon((i + 1) % boundary.numPoints());
+                int lon1 = GeoEncodingUtils.encodeLongitude(latLng1.getLonDeg());
+                int lat1 = GeoEncodingUtils.encodeLatitude(latLng1.getLatDeg());
+                int lon2 = GeoEncodingUtils.encodeLongitude(latLng2.getLonDeg());
+                int lat2 = GeoEncodingUtils.encodeLatitude(latLng2.getLatDeg());
+                if (isSharedBoundary(lon1, lat1, lon2, lat2, ringBoundary)) {
+                    count++;
+                }
+            }
+            assertThat("For cell " + H3.h3ToString(h3), count, either(equalTo(1)).or(equalTo(2)));
+        }
+    }
+
+    private boolean isSharedBoundary(int clon1, int clat1, int clon2, int clat2, CellBoundary boundary) {
+        for (int i = 0; i < boundary.numPoints(); i++) {
+            LatLng latLng1 = boundary.getLatLon(i % boundary.numPoints());
+            LatLng latLng2 = boundary.getLatLon((i + 1) % boundary.numPoints());
+            int lon1 = GeoEncodingUtils.encodeLongitude(latLng1.getLonDeg());
+            int lat1 = GeoEncodingUtils.encodeLatitude(latLng1.getLatDeg());
+            int lon2 = GeoEncodingUtils.encodeLongitude(latLng2.getLonDeg());
+            int lat2 = GeoEncodingUtils.encodeLatitude(latLng2.getLatDeg());
+            // edges are in opposite directions.
+            if (clon1 == lon2 & clat1 == lat2 && clon2 == lon1 && clat2 == lat1) {
+                return true;
+            }
         }
+        return false;
     }
 }