geodesic_missing_cells.c

#include <stdio.h>
#include <stdlib.h>

#include "h3Index.h"
#include "iterators.h"
#include "polyfill.h"
#include "polygon.h"
#include "test.h"
#include "utility.h"

// Owning array of H3 cells.
typedef struct {
    H3Index *cells;
    int64_t count;
} CellSet;

static CellSet cellSetNew(int64_t capacity) {
    CellSet s;
    s.cells = calloc(capacity, sizeof(H3Index));
    s.count = capacity;
    return s;
}

static void cellSetFree(CellSet *s) {
    free(s->cells);
    s->cells = NULL;
    s->count = 0;
}

// Remove zero entries in place.
static void cellSetRemoveZeros(CellSet *s) {
    int64_t out = 0;
    for (int64_t i = 0; i < s->count; i++) {
        if (s->cells[i] != 0) s->cells[out++] = s->cells[i];
    }
    s->count = out;
}

static void cellSetPrint(CellSet s, const char *label) {
    printf("%s (%lld cells):\n", label, (long long)s.count);
    for (int64_t i = 0; i < s.count; i++) {
        printf("  0x%015llx\n", (unsigned long long)s.cells[i]);
    }
}

// Check whether `b` is a superset of `a`.
static int cellSetIsSuperset(CellSet a, CellSet b) {
    for (int64_t i = 0; i < a.count; i++) {
        int found = 0;
        for (int64_t j = 0; j < b.count; j++) {
            if (b.cells[j] == a.cells[i]) {
                found = 1;
                break;
            }
        }
        if (!found) return 0;
    }
    return 1;
}

// Print cells in `a` that are not in `b`.
static void cellSetPrintMissing(CellSet a, CellSet b) {
    for (int64_t i = 0; i < a.count; i++) {
        int found = 0;
        for (int64_t j = 0; j < b.count; j++) {
            if (b.cells[j] == a.cells[i]) {
                found = 1;
                break;
            }
        }
        if (!found) {
            printf("  MISSING: 0x%015llx\n", (unsigned long long)a.cells[i]);
        }
    }
}

static CellSet grid_disk(H3Index center, int k) {
    int64_t maxSize;
    t_assertSuccess(H3_EXPORT(maxGridDiskSize)(k, &maxSize));
    CellSet s = cellSetNew(maxSize);
    t_assertSuccess(H3_EXPORT(gridDisk)(center, k, s.cells));
    cellSetRemoveZeros(&s);
    return s;
}

static CellSet geodesicPolyfill(const GeoPolygon *polygon, int res) {
    uint32_t flags = CONTAINMENT_OVERLAPPING;
    FLAG_SET_GEODESIC(flags);

    int64_t maxSize = 0;
    t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
        polygon, res, flags, &maxSize));

    CellSet s = cellSetNew(maxSize);
    t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(polygon, res, flags,
                                                          maxSize, s.cells));
    cellSetRemoveZeros(&s);
    return s;
}

// Convert cells to a single polygon with no holes.
static GeoMultiPolygon cellSetToPolygon(CellSet s) {
    GeoMultiPolygon mpoly = {0};
    t_assertSuccess(H3_EXPORT(cellsToMultiPolygon)(s.cells, s.count, &mpoly));
    t_assert(mpoly.numPolygons == 1, "single polygon");
    t_assert(mpoly.polygons[0].numHoles == 0, "no holes");
    return mpoly;
}

static void geoPolygonPrintGeoJSON(const GeoPolygon *polygon) {
    printf("{\n  \"type\": \"Feature\",\n  \"properties\": {},\n");
    printf("  \"geometry\": {\n    \"type\": \"Polygon\",\n");
    printf("    \"coordinates\": [[\n");
    const GeoLoop *loop = &polygon->geoloop;
    for (int i = 0; i <= loop->numVerts; i++) {
        const LatLng *v = &loop->verts[i % loop->numVerts];
        printf("      [%.15f, %.15f]%s\n", H3_EXPORT(radsToDegs)(v->lng),
               H3_EXPORT(radsToDegs)(v->lat), (i < loop->numVerts) ? "," : "");
    }
    printf("    ]]\n  }\n}\n");
}

static void gridDiskPolyfillCheck(H3Index center, int k) {
    CellSet disk = grid_disk(center, k);
    GeoMultiPolygon mpoly = cellSetToPolygon(disk);
    int res = H3_EXPORT(getResolution)(center);
    CellSet out = geodesicPolyfill(&mpoly.polygons[0], res);

    if (!cellSetIsSuperset(disk, out)) {
        printf("\n=== h = 0x%015llx,  k=%d ===\n", (unsigned long long)center,
               k);
        cellSetPrint(disk, "Input");
        cellSetPrint(out, "Output");
        cellSetPrintMissing(disk, out);
        geoPolygonPrintGeoJSON(&mpoly.polygons[0]);
    }

    cellSetFree(&out);
    cellSetFree(&disk);
    H3_EXPORT(destroyGeoMultiPolygon)(&mpoly);
}

SUITE(GeodesicGridDiskRoundTrip) {
    // TEST(failingCells) {
    //     int k = 1;
    //     gridDiskPolyfillCheck(0x8027fffffffffff, k);
    //     gridDiskPolyfillCheck(0x804bfffffffffff, k);
    //     gridDiskPolyfillCheck(0x80d3fffffffffff, k);
    // }

    TEST(allRes0k1) {
        for (IterCellsResolution it = iterInitRes(0); it.h; iterStepRes(&it)) {
            gridDiskPolyfillCheck(it.h, 1);
        }
    }

    // TEST(allRes0k2) {
    //     int count = H3_EXPORT(res0CellCount)();
    //     H3Index *res0 = calloc(count, sizeof(H3Index));
    //     t_assert(res0 != NULL, "alloc res0");
    //     t_assertSuccess(H3_EXPORT(getRes0Cells)(res0));
    //     for (int i = 0; i < count; i++) {
    //         gridDiskPolyfillCheck(res0[i], 2);
    //     }
    //     free(res0);
    // }

    // TEST(allRes1k1) {
    //     for (IterCellsResolution it = iterInitRes(1); it.h; iterStepRes(&it))
    //     {
    //         gridDiskPolyfillCheck(it.h, 1);
    //     }
    // }

    // TEST(allRes1k2) {
    //     for (IterCellsResolution it = iterInitRes(1); it.h; iterStepRes(&it))
    //     {
    //         gridDiskPolyfillCheck(it.h, 2);
    //     }
    // }
}