Change float to double in vec2d

CHANGELOG.md

@@ -6,6 +6,8 @@ The public API of this library consists of the functions declared in file
 [h3api.h.in](./src/h3lib/include/h3api.h.in).
 
 ## [Unreleased]
+### Fixed
+- Changing an internal `float` to `double` improves the precision of geographic coordinate output (#652)
 
 ## [4.0.0] - 2022-08-23
 ### Breaking changes

CMakeLists.txt

@@ -190,6 +190,7 @@ set(OTHER_SOURCE_FILES
     src/apps/testapps/testPentagonIndexes.c
     src/apps/testapps/testGridDisk.c
     src/apps/testapps/testCellToBoundary.c
+    src/apps/testapps/testCellToBoundaryEdgeCases.c
     src/apps/testapps/testCellToParent.c
     src/apps/testapps/testH3Index.c
     src/apps/testapps/mkRandGeoBoundary.c

CMakeTests.cmake

@@ -167,6 +167,7 @@ foreach(file ${all_cells})
     add_h3_test_with_file(testCellToBoundary src/apps/testapps/testCellToBoundary.c ${file})
 endforeach()
 
+add_h3_test(testCellToBoundaryEdgeCases src/apps/testapps/testCellToBoundaryEdgeCases.c)
 add_h3_test(testCompactCells src/apps/testapps/testCompactCells.c)
 add_h3_test(testGridDisk src/apps/testapps/testGridDisk.c)
 add_h3_test(testGridRingUnsafe src/apps/testapps/testGridRingUnsafe.c)

src/apps/miscapps/README.md

@@ -0,0 +1,18 @@
+# miscapps
+
+These apps generate data, either static tables used in the library or test data used in regression tests.
+
+## Generate Test Data
+
+For a specific base cell at a specific res:
+
+    ./bin/cellToLatLngHier -p 800bfffffffffff -r 8 > tests/inputfiles/bc05r08centers.txt
+    ./bin/cellToBoundaryHier -p 800bfffffffffff -r 8 > tests/inputfiles/bc05r08cells.txt
+
+For all cells at a given res:
+
+    ./bin/h3ToHier -r 0 | xargs -I {} ./bin/cellToBoundaryHier -p {} -r 3 > tests/inputfiles/res03cells.txt
+
+For a random sample of cells at a given res:
+
+    ./bin/mkRandGeoBoundary -n 5000 -r 5 > tests/inputfiles/rand05cells.txt

src/apps/testapps/testCellToBoundaryEdgeCases.c

@@ -0,0 +1,60 @@
+/*
+ * Copyright 2017-2021 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <float.h>
+#include <stdlib.h>
+
+#include "bbox.h"
+#include "h3api.h"
+#include "latLng.h"
+#include "polygon.h"
+#include "test.h"
+#include "utility.h"
+
+SUITE(cellToBoundaryEdgeCases) {
+    TEST(doublePrecisionVertex) {
+        // The carefully constructed case here:
+        // - A res 1 pentagon cell with distortion vertexes that change
+        //   when we use a double instead of a float in _v2dIntersect
+        // - One of the previous (float-based) distortion vertexes
+        // This is the only case yet found where a point indexed to the
+        // cell is shown to be incorrectly outside of the geo boundary
+        // when we use the float version. Presumably more could be found.
+        H3Index cell = 0x81083ffffffffff;
+        LatLng point = {.lat = H3_EXPORT(degsToRads)(61.890838431),
+                        .lng = H3_EXPORT(degsToRads)(8.644221328)};
+
+        CellBoundary boundary;
+        t_assertSuccess(H3_EXPORT(cellToBoundary)(cell, &boundary));
+
+        LatLng *verts = boundary.verts;
+        GeoLoop geoloop = {.numVerts = boundary.numVerts, .verts = verts};
+
+        BBox bbox;
+        bboxFromGeoLoop(&geoloop, &bbox);
+
+        H3Index cell2;
+        t_assertSuccess(H3_EXPORT(latLngToCell)(&point, 1, &cell2));
+        // Check whether the point is physically inside the geo boundary
+        if (cell2 == cell) {
+            t_assert(pointInsideGeoLoop(&geoloop, &bbox, &point),
+                     "Boundary contains input point");
+        } else {
+            t_assert(!pointInsideGeoLoop(&geoloop, &bbox, &point),
+                     "Boundary does not contain input point");
+        }
+    }
+}

src/apps/testapps/testCellsToLinkedMultiPolygon.c

@@ -308,4 +308,69 @@ SUITE(cellsToLinkedMultiPolygon) {
                                                       &polygon) == E_FAILED,
                  "invalid cells fail");
     }
+
+    TEST(gridDiskResolutions) {
+        // This is a center-face base cell, no pentagon siblings
+        H3Index baseCell = 0x8073fffffffffff;
+        H3Index origin = baseCell;
+
+        H3Index indexes[19] = {0};
+        int numHexes = 7;
+
+        for (int res = 1; res < 15; res++) {
+            // Take the 2-disk of the center child at res
+            t_assertSuccess(
+                H3_EXPORT(cellToCenterChild)(baseCell, res, &origin));
+            t_assertSuccess(H3_EXPORT(gridDisk)(origin, 2, indexes));
+
+            // Test the polygon output
+            LinkedGeoPolygon polygon;
+            t_assertSuccess(H3_EXPORT(cellsToLinkedMultiPolygon)(
+                indexes, numHexes, &polygon));
+            t_assert(countLinkedPolygons(&polygon) == 1, "1 polygon added");
+            t_assert(countLinkedLoops(&polygon) == 1,
+                     "1 loop on the first polygon");
+            t_assert(countLinkedCoords(polygon.first) == 18,
+                     "All coords for all hexes added to first loop");
+
+            H3_EXPORT(destroyLinkedMultiPolygon)(&polygon);
+        }
+    }
+
+    TEST(gridDiskResolutionsPentagon) {
+        // This is a pentagon base cell
+        H3Index baseCell = 0x8031fffffffffff;
+        H3Index origin = baseCell;
+
+        H3Index diskIndexes[7] = {0};
+        H3Index indexes[6] = {0};
+
+        for (int res = 1; res < 15; res++) {
+            // Take the 1-disk of the center child at res. Note: We can't take
+            // the 2-disk here, as increased distortion around the pentagon will
+            // still fail at res 1. TODO: Use a 2-ring, start at res 0
+            // when output is correct.
+            t_assertSuccess(
+                H3_EXPORT(cellToCenterChild)(baseCell, res, &origin));
+            t_assertSuccess(H3_EXPORT(gridDisk)(origin, 1, diskIndexes));
+
+            int j = 0;
+            for (int i = 0; i < 7; i++) {
+                if (diskIndexes[i]) indexes[j++] = diskIndexes[i];
+            }
+            t_assert(j == 6, "Filled all 6 indexes");
+
+            // Test the polygon output
+            LinkedGeoPolygon polygon;
+            t_assertSuccess(
+                H3_EXPORT(cellsToLinkedMultiPolygon)(indexes, 6, &polygon));
+            t_assert(countLinkedPolygons(&polygon) == 1, "1 polygon added");
+            t_assert(countLinkedLoops(&polygon) == 1,
+                     "1 loop on the first polygon");
+            t_assert(countLinkedCoords(polygon.first) == 15,
+                     "All coords for all hexes added to first loop");
+
+            H3_EXPORT(destroyLinkedMultiPolygon)(&polygon);
+        }
+    }
 }

src/apps/testapps/testH3Api.c

@@ -55,16 +55,14 @@ SUITE(h3Api) {
 
     TEST(cellToBoundary_classIIIEdgeVertex) {
         // Bug test for https://github.com/uber/h3/issues/45
-        char *hexes[] = {
-            "894cc5349b7ffff", "894cc534d97ffff", "894cc53682bffff",
-            "894cc536b17ffff", "894cc53688bffff", "894cead92cbffff",
-            "894cc536537ffff", "894cc5acbabffff", "894cc536597ffff"};
+        H3Index hexes[] = {
+            0x894cc5349b7ffff, 0x894cc534d97ffff, 0x894cc53682bffff,
+            0x894cc536b17ffff, 0x894cc53688bffff, 0x894cead92cbffff,
+            0x894cc536537ffff, 0x894cc5acbabffff, 0x894cc536597ffff};
         int numHexes = sizeof(hexes) / sizeof(hexes[0]);
-        H3Index h3;
         CellBoundary b;
         for (int i = 0; i < numHexes; i++) {
-            t_assertSuccess(H3_EXPORT(stringToH3)(hexes[i], &h3));
-            H3_EXPORT(cellToBoundary)(h3, &b);
+            H3_EXPORT(cellToBoundary)(hexes[i], &b);
             t_assert(b.numVerts == 7, "got expected vertex count");
         }
     }

src/apps/testapps/testH3CellArea.c

@@ -26,7 +26,7 @@
 #include "test.h"
 
 static const double areasKm2[] = {
-    2.562182162955496e+06, 4.476842018179411e+05, 6.596162242711056e+04,
+    2.562182162955496e+06, 4.476842017201860e+05, 6.596162242711056e+04,
     9.228872919002590e+03, 1.318694490797110e+03, 1.879593512281298e+02,
     2.687164354763186e+01, 3.840848847060638e+00, 5.486939641329893e-01,
     7.838600808637444e-02, 1.119834221989390e-02, 1.599777169186614e-03,

src/apps/testapps/testVec2d.c

@@ -47,14 +47,16 @@ SUITE(Vec2d) {
                  "Y coord as expected");
     }
 
-    TEST(_v2dEquals) {
+    TEST(_v2dAlmostEquals) {
         Vec2d v1 = {3.0, 4.0};
         Vec2d v2 = {3.0, 4.0};
         Vec2d v3 = {3.5, 4.0};
         Vec2d v4 = {3.0, 4.5};
+        Vec2d v5 = {3.0 + DBL_EPSILON, 4.0 - DBL_EPSILON};
 
-        t_assert(_v2dEquals(&v1, &v2), "true for equal vectors");
-        t_assert(!_v2dEquals(&v1, &v3), "false for different x");
-        t_assert(!_v2dEquals(&v1, &v4), "false for different y");
+        t_assert(_v2dAlmostEquals(&v1, &v2), "true for equal vectors");
+        t_assert(!_v2dAlmostEquals(&v1, &v3), "false for different x");
+        t_assert(!_v2dAlmostEquals(&v1, &v4), "false for different y");
+        t_assert(_v2dAlmostEquals(&v1, &v5), "true for almost equal");
     }
 }

src/h3lib/include/vec2d.h

@@ -35,6 +35,6 @@ typedef struct {
 double _v2dMag(const Vec2d *v);
 void _v2dIntersect(const Vec2d *p0, const Vec2d *p1, const Vec2d *p2,
                    const Vec2d *p3, Vec2d *inter);
-bool _v2dEquals(const Vec2d *p0, const Vec2d *p1);
+bool _v2dAlmostEquals(const Vec2d *p0, const Vec2d *p1);
 
 #endif

src/h3lib/lib/faceijk.c

@@ -747,8 +747,8 @@ void _faceIjkToCellBoundary(const FaceIJK *h, int res, int start, int length,
             adjacent hexagon edge will lie completely on a single icosahedron
             face, and no additional vertex is required.
             */
-            bool isIntersectionAtVertex =
-                _v2dEquals(&orig2d0, &inter) || _v2dEquals(&orig2d1, &inter);
+            bool isIntersectionAtVertex = _v2dAlmostEquals(&orig2d0, &inter) ||
+                                          _v2dAlmostEquals(&orig2d1, &inter);
             if (!isIntersectionAtVertex) {
                 _hex2dToGeo(&inter, centerIJK.face, adjRes, 1,
                             &g->verts[g->numVerts]);

src/h3lib/lib/vec2d.c

@@ -19,6 +19,7 @@
 
 #include "vec2d.h"
 
+#include <float.h>
 #include <math.h>
 #include <stdio.h>
 
@@ -46,7 +47,7 @@ void _v2dIntersect(const Vec2d *p0, const Vec2d *p1, const Vec2d *p2,
     s2.x = p3->x - p2->x;
     s2.y = p3->y - p2->y;
 
-    float t;
+    double t;
     t = (s2.x * (p0->y - p2->y) - s2.y * (p0->x - p2->x)) /
         (-s2.x * s1.y + s1.x * s2.y);
 
@@ -55,12 +56,12 @@ void _v2dIntersect(const Vec2d *p0, const Vec2d *p1, const Vec2d *p2,
 }
 
 /**
- * Whether two 2D vectors are equal. Does not consider possible false
- * negatives due to floating-point errors.
+ * Whether two 2D vectors are almost equal, within some threshold
  * @param v1 First vector to compare
  * @param v2 Second vector to compare
- * @return Whether the vectors are equal
+ * @return Whether the vectors are almost equal
  */
-bool _v2dEquals(const Vec2d *v1, const Vec2d *v2) {
-    return v1->x == v2->x && v1->y == v2->y;
+bool _v2dAlmostEquals(const Vec2d *v1, const Vec2d *v2) {
+    return fabs(v1->x - v2->x) < FLT_EPSILON &&
+           fabs(v1->y - v2->y) < FLT_EPSILON;
 }