h3_polyfill.c

/*
 * Copyright 2023 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.
 */
/** @file polyfill.c
 * @brief   Functions relating to the cell-to-polygon algorithm
 */

#include "h3_polyfill.h"

#include <math.h>
#include <string.h>

#include "h3_alloc.h"
#include "h3_baseCells.h"
#include "h3_coordijk.h"
#include "h3_h3Assert.h"
#include "h3_h3Index.h"
#include "h3_polygon.h"

// Factor by which to scale the cell bounding box to include all cells.
// This was determined empirically by finding the smallest factor that
// passed exhaustive tests.
#define CELL_SCALE_FACTOR 1.1

// Factor by which to scale the cell bounding box to include all children.
// This was determined empirically by finding the smallest factor that
// passed exhaustive tests.
#define CHILD_SCALE_FACTOR 1.4

/**
 * Max cell edge length, in radians, for each resolution. This was computed
 * by taking the max exact edge length for cells at the center of each base
 * cell at that resolution.
 */
static double MAX_EDGE_LENGTH_RADS[MAX_H3_RES + 1] = {
    0.21577206265130, 0.08308767068495, 0.03148970436439, 0.01190662871439,
    0.00450053330908, 0.00170105523619, 0.00064293917678, 0.00024300820659,
    0.00009184847087, 0.00003471545901, 0.00001312121017, 0.00000495935129,
    0.00000187445860, 0.00000070847876, 0.00000026777980, 0.00000010121125};

/** All cells that contain the north pole, by res */
static H3Index NORTH_POLE_CELLS[MAX_H3_RES + 1] = {
    0x8001fffffffffff, 0x81033ffffffffff, 0x820327fffffffff, 0x830326fffffffff,
    0x8403263ffffffff, 0x85032623fffffff, 0x860326237ffffff, 0x870326233ffffff,
    0x880326233bfffff, 0x890326233abffff, 0x8a0326233ab7fff, 0x8b0326233ab0fff,
    0x8c0326233ab03ff, 0x8d0326233ab03bf, 0x8e0326233ab039f, 0x8f0326233ab0399};

/** All cells that contain the south pole, by res */
static H3Index SOUTH_POLE_CELLS[MAX_H3_RES + 1] = {
    0x80f3fffffffffff, 0x81f2bffffffffff, 0x82f297fffffffff, 0x83f293fffffffff,
    0x84f2939ffffffff, 0x85f29383fffffff, 0x86f29380fffffff, 0x87f29380effffff,
    0x88f29380e1fffff, 0x89f29380e0fffff, 0x8af29380e0d7fff, 0x8bf29380e0d0fff,
    0x8cf29380e0d0dff, 0x8df29380e0d0cff, 0x8ef29380e0d0cc7, 0x8ff29380e0d0cc4};

/** Pre-calculated bounding boxes for all res 0 cells */
static BBox RES0_BBOXES[NUM_BASE_CELLS] = {
    {1.52480158339146, 1.20305471830087, -0.60664883654036, 0.00568297271999},
    {1.52480158339146, 1.17872424267511, -0.60664883654036, 2.54046980298264},
    {1.52480158339146, 1.09069387298096, -2.85286053297673, 1.64310689027893},
    {1.41845302535151, 1.01285145697208, 0.00568297271999, -1.16770379632602},
    {1.27950477868453, 0.97226652536306, 0.55556064983494, -0.18229924845326},
    {1.32929586572429, 0.91898920750071, 2.05622344943192, 1.08813154278274},
    {1.32899086063916, 0.94271815376360, -2.29875289606378, 3.01700008041993},
    {1.26020983864103, 0.84291228415618, -0.89971867664861, -1.75967359310997},
    {1.21114673854945, 0.86170600921069, 1.19129757609455, 0.43777608996454},
    {1.21075831414294, 0.83795331049498, -1.72022875779891, -2.43793861727138},
    {1.15546530929588, 0.78982455384253, 2.53659412229266, 1.85709133451243},
    {1.15528445067052, 0.76641428724335, -3.06738507202411, 2.53646110244042},
    {1.10121643537669, 0.71330093663066, 0.09640581900154, -0.52154514518248},
    {1.07042472765165, 0.67603948819406, -0.47984202840088, -1.10306159603090},
    {1.03270228748960, 0.72356358827215, -2.24990138725146, -2.74510220919157},
    {1.01929924623886, 0.65491232835426, 0.63035574240731, 0.03537030096470},
    {1.01786037568858, 0.58827636737638, 1.53192721817065, 0.93672682511233},
    {0.98081434136020, 0.61076063532947, -2.67100636598529, 3.06516463008733},
    {0.98106023192774, 0.58679836571570, 2.02829766214461, 1.51334374970280},
    {0.96374551790056, 0.55186491737474, -1.42976721313659, -1.96852202530104},
    {0.87536136210723, 0.50008952762292, -1.92435613571430, -2.41641343219793},
    {0.88611243445554, 0.52742963716774, -0.95781946324194, -1.47628966305930},
    {0.86881343251986, 0.50770567021439, 1.03236795495839, 0.50347284027426},
    {0.89235638181782, 0.48781264892508, 2.76430302119150, 2.29989716697031},
    {0.82570569254601, 0.52173101741059, 2.30921681461428, 1.93198541828980},
    {0.80599330438546, 0.40150819579319, -3.06417559403240, 2.70079300784409},
    {0.81612079704781, 0.38396800633226, -0.21614378891839, -0.70420149722178},
    {0.75822779851431, 0.39943555383751, -2.34059978084699, -2.82127373822444},
    {0.78861390967531, 0.38742018303868, 0.23115687731652, -0.22599491086066},
    {0.71515840341957, 0.33012478438475, -0.64847976163163, -1.08249728121219},
    {0.70359051048414, 0.29148673180722, 1.71441081857246, 1.28443348381696},
    {0.69190629544818, 0.28808313184381, 0.64863909244647, 0.16372369282557},
    {0.64863235654749, 0.26290420067147, 2.10318098268379, 1.69556122548344},
    {0.65722892279906, 0.28222653310929, 1.30918693285466, 0.87594416271685},
    {0.64750997738584, 0.24149865709850, -1.30272192474556, -1.68708570163242},
    {0.62380174028378, 0.25522080363509, -2.72428423026826, 3.10401473237630},
    {0.64228460410023, 0.21206753429148, -1.67639240992071, -2.11772366767341},
    {0.59919175361146, 0.21620460836570, 2.48592868387690, 2.07350353893591},
    {0.55637406851384, 0.25276557437230, -0.99885388505694, -1.32642489358939},
    {0.55648013300665, 0.15187401321019, 2.87032088421324, 2.44642320475367},
    {0.54603687970450, 0.15589091511369, -2.06789866067060, -2.49091419631961},
    {0.51206347752697, 0.15522020377124, 0.95446767315996, 0.54443262110414},
    {0.49767951537101, 0.10944898890579, -0.04335162263358, -0.42900268178569},
    {0.46538045483671, 0.06029968637720, -0.41240613713421, -0.80603623808166},
    {0.44686891066946, 0.06926857458503, 0.32053284794952, -0.07005748900849},
    {0.43208958202064, 0.07796440938140, -3.06232453079660, 2.80602499990282},
    {0.43103892586713, 0.02927431919853, -2.41589238618422, -2.85735809951951},
    {0.38073727558986, -0.00297016159959, -0.77039553861218, -1.14788248745028},
    {0.39113816687141, -0.01518764903038, 1.49130246958290, 1.14714731736311},
    {0.33421063142418, 0.02526613430348, 1.15141032578749, 0.85000706261644},
    {0.38915669778582, -0.04371359825454, 1.88046353933242, 1.48230231380717},
    {0.33787520825987, -0.04835090128296, -1.12274014380603, -1.49454408844749},
    {0.33601418932337, -0.06675068178541, 2.23792354204464, 1.85723423013211},
    {0.31838318078049, -0.05821955623722, 0.66058854060373, 0.25452572938783},
    {0.33630761471457, -0.07589541031521, -1.47957331741818, -1.85981735718264},
    {0.28924817322870, -0.09150638064667, -1.83561930288569, -2.21855897384292},
    {0.26678632252475, -0.10058088990867, -2.76808651991421, 3.12792953247061},
    {0.29285254112587, -0.13483165093783, 2.61406468380434, 2.20466422911705},
    {0.20150342788824, -0.10279852729762, 0.06881896344365, -0.23925229432978},
    {0.21283813275258, -0.18626835417891, 2.93800440256577, 2.57470747655623},
    {0.19587614179884, -0.17237030304155, -2.16941795427335, -2.55405165906601},
    {0.17237030304155, -0.19587614179884, 0.97217469931645, 0.58754099452378},
    {0.18626835417891, -0.21283813275258, -0.20358825102402, -0.56688517703356},
    {0.10279852729762, -0.20150342788824, -3.07277369014614, 2.90234035926002},
    {0.13483165093783, -0.29285254112587, -0.52752796978545, -0.93692842447275},
    {0.10058088990867, -0.26678632252475, 0.37350613367558, -0.01366312111919},
    {0.09150638064667, -0.28924817322870, 1.30597335070410, 0.92303367974687},
    {0.07589541031521, -0.33630761471457, 1.66201933617161, 1.28177529640715},
    {0.05821955623722, -0.31838318078049, -2.48100411298606, -2.88706692420196},
    {0.06675068178541, -0.33601418932337, -0.90366911154516, -1.28435842345769},
    {0.04835090128296, -0.33787520825987, 2.01885250978376, 1.64704856514230},
    {0.04371359825454, -0.38915669778582, -1.26112911425737, -1.65929033978262},
    {-0.02526613430348, -0.33421063142418, -1.99018232780231,
     -2.29158559097336},
    {0.01518764903038, -0.39113816687140, -1.65029018400690, -1.99444533622668},
    {0.00297016159959, -0.38073727558986, 2.37119711497761, 1.99371016613951},
    {-0.02927431919853, -0.43103892586713, 0.72570026740558, 0.28423455407029},
    {-0.07796440938140, -0.43208958202064, 0.07926812279319, -0.33556765368697},
    {-0.06926857458503, -0.44686891066946, -2.82105980564027, 3.07153516458131},
    {-0.06029968637720, -0.46538045483671, 2.72918651645558, 2.33555641550814},
    {-0.10944898890579, -0.49767951537101, 3.09824103095621, 2.71258997180410},
    {-0.15522020377124, -0.51206347752697, -2.18712498042983,
     -2.59716003248565},
    {-0.15589091511369, -0.54603687970450, 1.07369399291919, 0.65067845727018},
    {-0.15187401321019, -0.55648013300665, -0.27127176937655,
     -0.69516944883612},
    {-0.25276557437230, -0.55637406851385, 2.14273876853285, 1.81516776000041},
    {-0.21620460836570, -0.59919175361146, -0.65566396971290,
     -1.06808911465388},
    {-0.21206753429148, -0.64228460410023, 1.46520024366909, 1.02386898591638},
    {-0.25522080363509, -0.62380174028378, 0.41730842332153, -0.03757792121350},
    {-0.24149865709850, -0.64750997738584, 1.83887072884423, 1.45450695195737},
    {-0.28222653310929, -0.65722892279906, -1.83240572073513,
     -2.26564849087294},
    {-0.26290420067147, -0.64863235654749, -1.03841167090601,
     -1.44603142810635},
    {-0.28808313184381, -0.69190629544818, -2.49295356114332,
     -2.97786896076422},
    {-0.29148673180722, -0.70359051048414, -1.42718183501734,
     -1.85715916977284},
    {-0.33012478438475, -0.71515840341957, 2.49311289195816, 2.05909537237761},
    {-0.38742018303868, -0.78861390967531, -2.91043577627328, 2.91559774272914},
    {-0.39943555383751, -0.75822779851431, 0.80099287274280, 0.32031891536535},
    {-0.38396800633226, -0.81612079704781, 2.92544886467140, 2.43739115636801},
    {-0.40150819579319, -0.80599330438546, 0.07741705955739, -0.44079964574570},
    {-0.52173101741059, -0.82570569254601, -0.83237583897551,
     -1.20960723529999},
    {-0.48781264892508, -0.89235638181782, -0.37728963239830,
     -0.84169548661948},
    {-0.50770567021439, -0.86881343251986, -2.10922469863141,
     -2.63811981331554},
    {-0.52742963716774, -0.88611243445554, 2.18377319034785, 1.66530299053050},
    {-0.50008952762292, -0.87536136210723, 1.21723651787549, 0.72517922139186},
    {-0.55186491737474, -0.96374551790056, 1.71182544045320, 1.17307062828876},
    {-0.58679836571570, -0.98106023192774, -1.11329499144518,
     -1.62824890388699},
    {-0.61076063532947, -0.98081434136020, 0.47058628760450, -0.07642802350246},
    {-0.58827636737638, -1.01786037568858, -1.60966543541914,
     -2.20486582847747},
    {-0.65491232835426, -1.01929924623886, -2.51123691118248,
     -3.10622235262510},
    {-0.72356358827215, -1.03270228748960, 0.89169126633833, 0.39649044439822},
    {-0.67603948819406, -1.07042472765165, 2.66175062518892, 2.03853105755889},
    {-0.71330093663066, -1.10121643537669, -3.04518683458825, 2.62004750840731},
    {-0.76641428724335, -1.15528445067052, 0.07420758156568, -0.60513155114938},
    {-0.78982455384253, -1.15546530929588, -0.60499853129713,
     -1.28450131907736},
    {-0.83795331049498, -1.21075831414294, 1.42136389579088, 0.70365403631841},
    {-0.86170600921069, -1.21114673854945, -1.95029507749525,
     -2.70381656362525},
    {-0.84291228415618, -1.26020983864103, 2.24187397694118, 1.38191906047983},
    {-0.94271815376360, -1.32899086063916, 0.84283975752601, -0.12459257316986},
    {-0.91898920750071, -1.32929586572429, -1.08536920415787,
     -2.05346111080706},
    {-0.97226652536306, -1.27950477868453, -2.58603200375485, 2.95929340513654},
    {-1.01285145697208, -1.41845302535151, -3.13590968086981, 1.97388885726377},
    {-1.09069387298096, -1.52480158339146, 0.28873212061306, -1.49848576331087},
    {-1.17872424267511, -1.52480158339146, 2.53494381704943, -0.60112285060716},
    {-1.20305471830087, -1.52480158339146, -0.60112285060716,
     2.53494381704943}};

static BBox VALID_RANGE_BBOX = {M_PI_2, -M_PI_2, M_PI, -M_PI};

/**
 * For a given cell, return its bounding box. If coverChildren is true, the bbox
 * will be guaranteed to contain its children at any finer resolution. Note that
 * no guarantee is provided as to the level of accuracy, and the bounding box
 * may have a significant margin of error.
 * @param cell Cell to calculate bbox for
 * @param out  BBox to hold output
 * @param coverChildren Whether the bounding box should cover all children
 */
H3Error cellToBBox(H3Index cell, BBox *out, bool coverChildren) {
    // Adjust the BBox to handle poles, if needed
    int res = H3_GET_RESOLUTION(cell);

    if (res == 0) {
        int baseCell = H3_GET_BASE_CELL(cell);
        if (NEVER(baseCell < 0) || baseCell >= NUM_BASE_CELLS) {
            return E_CELL_INVALID;
        }
        *out = RES0_BBOXES[baseCell];
    } else {
        LatLng center;
        H3Error centerErr = H3_EXPORT(cellToLatLng)(cell, &center);
        if (centerErr != E_SUCCESS) {
            return centerErr;
        }
        double lngRatio = 1 / cos(center.lat);
        out->north = center.lat + MAX_EDGE_LENGTH_RADS[res];
        out->south = center.lat - MAX_EDGE_LENGTH_RADS[res];
        out->east = center.lng + MAX_EDGE_LENGTH_RADS[res] * lngRatio;
        out->west = center.lng - MAX_EDGE_LENGTH_RADS[res] * lngRatio;
    }

    // Buffer the bounding box to cover children. Call this even if no buffering
    // is required in order to normalize the bbox to lat/lng bounds
    scaleBBox(out, coverChildren ? CHILD_SCALE_FACTOR : CELL_SCALE_FACTOR);

    // Cell that contains the north pole
    if (cell == NORTH_POLE_CELLS[res]) {
        out->north = M_PI_2;
    }

    // Cell that contains the south pole
    if (cell == SOUTH_POLE_CELLS[res]) {
        out->south = -M_PI_2;
    }

    // If we contain a pole, expand the longitude to include the full domain,
    // effectively making the bbox a circle around the pole.
    if (out->north == M_PI_2 || out->south == -M_PI_2) {
        out->east = M_PI;
        out->west = -M_PI;
    }

    return E_SUCCESS;
}

/**
 * Get a base cell by number, or H3_NULL if out of bounds
 */
H3Index baseCellNumToCell(int baseCellNum) {
    if (baseCellNum < 0 || baseCellNum >= NUM_BASE_CELLS) {
        return H3_NULL;
    }
    H3Index baseCell;
    setH3Index(&baseCell, 0, baseCellNum, 0);
    return baseCell;
}

static void iterErrorPolygonCompact(IterCellsPolygonCompact *iter,
                                    H3Error error) {
    iterDestroyPolygonCompact(iter);
    iter->error = error;
}

/**
 * Given a cell, find the next cell in the sequence of all cells
 * to check in the iteration.
 */
static H3Index nextCell(H3Index cell) {
    int res = H3_GET_RESOLUTION(cell);
    while (true) {
        // If this is a base cell, set to next base cell (or H3_NULL if done)
        if (res == 0) {
            return baseCellNumToCell(H3_GET_BASE_CELL(cell) + 1);
        }

        // Faster cellToParent when we know the resolution is valid
        // and we're only moving up one level
        H3Index parent = cell;
        H3_SET_RESOLUTION(parent, res - 1);
        H3_SET_INDEX_DIGIT(parent, res, H3_DIGIT_MASK);

        // If not the last sibling of parent, return next sibling
        Direction digit = H3_GET_INDEX_DIGIT(cell, res);
        if (digit < INVALID_DIGIT - 1) {
            H3_SET_INDEX_DIGIT(cell, res,
                               digit + ((H3_EXPORT(isPentagon)(parent) &&
                                         digit == CENTER_DIGIT)
                                            ? 2  // Skip missing pentagon child
                                            : 1));
            return cell;
        }
        // Move up to the parent for the next loop iteration
        res--;
        cell = parent;
    }
}

/**
 * Internal function - initialize the iterator without stepping to the first
 * value
 */
static IterCellsPolygonCompact _iterInitPolygonCompact(
    const GeoPolygon *polygon, int res, uint32_t flags) {
    IterCellsPolygonCompact iter = {// Initialize output properties. The first
                                    // valid cell will be set in iterStep
                                    .cell = baseCellNumToCell(0),
                                    .error = E_SUCCESS,
                                    // Save input arguments
                                    ._polygon = polygon,
                                    ._res = res,
                                    ._flags = flags,
                                    ._bboxes = NULL,
                                    ._started = false};

    if (res < 0 || res > MAX_H3_RES) {
        iterErrorPolygonCompact(&iter, E_RES_DOMAIN);
        return iter;
    }

    H3Error flagErr = validatePolygonFlags(flags);
    if (flagErr) {
        iterErrorPolygonCompact(&iter, flagErr);
        return iter;
    }

    // Initialize bounding boxes for polygon and any holes. Memory allocated
    // here must be released through iterDestroyPolygonCompact
    iter._bboxes = H3_MEMORY(calloc)((polygon->numHoles + 1), sizeof(BBox));
    if (!iter._bboxes) {
        iterErrorPolygonCompact(&iter, E_MEMORY_ALLOC);
        return iter;
    }
    bboxesFromGeoPolygon(polygon, iter._bboxes);

    return iter;
}

/**
 * Initialize a IterCellsPolygonCompact struct representing the sequence of
 * compact cells within the target polygon. The test for including edge cells is
 * defined by the polyfill mode passed in the `flags` argument.
 *
 * Initialization of this object may fail, in which case the `error` property
 * will be set and all iteration will return H3_NULL. It is the responsibility
 * of the caller to check the error property after initialization.
 *
 * At any point in the iteration, starting once the struct is initialized, the
 * output value can be accessed through the `cell` property.
 *
 * Note that initializing the iterator allocates memory. If an iterator is
 * exhausted or returns an error that memory is released; otherwise it must be
 * released manually with iterDestroyPolygonCompact.
 *
 * @param  polygon Polygon to fill with compact cells
 * @param  res     Finest resolution for output cells
 * @param  flags   Bit mask of option flags
 * @return         Initialized iterator, with the first value available
 */
IterCellsPolygonCompact iterInitPolygonCompact(const GeoPolygon *polygon,
                                               int res, uint32_t flags) {
    IterCellsPolygonCompact iter = _iterInitPolygonCompact(polygon, res, flags);

    // Start the iterator by taking the first step.
    // This is necessary to have a valid value after initialization.
    iterStepPolygonCompact(&iter);

    return iter;
}

/**
 * Increment the polyfill iterator, running the polygon to cells algorithm.
 *
 * Briefly, the algorithm checks every cell in the global grid hierarchically,
 * starting with the base cells. Cells coarser than the target resolution are
 * checked for complete child inclusion using a bounding box guaranteed to
 * contain all children.
 * - If the bounding box is contained by the polygon, output is set to the cell
 * - If the bounding box intersects, recurse into the first child
 * - Otherwise, continue with the next cell in sequence
 *
 * For cells at the target resolution, a finer-grained check is used according
 * to the inclusion criteria set in flags.
 *
 * @param  iter Iterator to increment
 */
void iterStepPolygonCompact(IterCellsPolygonCompact *iter) {
    H3Index cell = iter->cell;

    // once the cell is H3_NULL, the iterator returns an infinite sequence of
    // H3_NULL
    if (cell == H3_NULL) return;

    // For the first step, we need to evaluate the current cell; after that, we
    // should start with the next cell.
    if (iter->_started) {
        cell = nextCell(cell);
    } else {
        iter->_started = true;
    }

    // Short-circuit iteration for 0-vert polygon
    if (iter->_polygon->geoloop.numVerts == 0) {
        iterDestroyPolygonCompact(iter);
        return;
    }

    ContainmentMode mode = FLAG_GET_CONTAINMENT_MODE(iter->_flags);

    while (cell) {
        int cellRes = H3_GET_RESOLUTION(cell);

        // Target res: Do a fine-grained check
        if (cellRes == iter->_res) {
            if (mode == CONTAINMENT_CENTER || mode == CONTAINMENT_OVERLAPPING ||
                mode == CONTAINMENT_OVERLAPPING_BBOX) {
                // Check if the cell center is inside the polygon
                LatLng center;
                H3Error centerErr = H3_EXPORT(cellToLatLng)(cell, &center);
                if (centerErr != E_SUCCESS) {
                    iterErrorPolygonCompact(iter, centerErr);
                    return;
                }
                if (pointInsidePolygon(iter->_polygon, iter->_bboxes,
                                       &center)) {
                    // Set to next output
                    iter->cell = cell;
                    return;
                }
            }
            if (mode == CONTAINMENT_OVERLAPPING ||
                mode == CONTAINMENT_OVERLAPPING_BBOX) {
                // For overlapping, we need to do a quick check to determine
                // whether the polygon is wholly contained by the cell. We
                // check the first polygon vertex, which if it is contained
                // could also mean we simply intersect.

                // Deferencing verts[0] is safe because we check numVerts above
                LatLng firstVertex = iter->_polygon->geoloop.verts[0];

                // We have to check whether the point is in the expected range
                // first, because out-of-bounds values will yield false
                // positives with latLngToCell
                if (bboxContains(&VALID_RANGE_BBOX, &firstVertex)) {
                    H3Index polygonCell;
                    H3Error polygonCellErr = H3_EXPORT(latLngToCell)(
                        &(iter->_polygon->geoloop.verts[0]), cellRes,
                        &polygonCell);
                    if (NEVER(polygonCellErr != E_SUCCESS)) {
                        // This should be unreachable with the bbox check
                        iterErrorPolygonCompact(iter, polygonCellErr);
                        return;
                    }
                    if (polygonCell == cell) {
                        // Set to next output
                        iter->cell = cell;
                        return;
                    }
                }
            }
            if (mode == CONTAINMENT_FULL || mode == CONTAINMENT_OVERLAPPING ||
                mode == CONTAINMENT_OVERLAPPING_BBOX) {
                CellBoundary boundary;
                H3Error boundaryErr =
                    H3_EXPORT(cellToBoundary)(cell, &boundary);
                if (boundaryErr != E_SUCCESS) {
                    iterErrorPolygonCompact(iter, boundaryErr);
                    return;
                }
                BBox bbox;
                H3Error bboxErr = cellToBBox(cell, &bbox, false);
                if (NEVER(bboxErr != E_SUCCESS)) {
                    // Should be unreachable - invalid cells would be caught in
                    // the previous boundaryErr
                    iterErrorPolygonCompact(iter, bboxErr);
                    return;
                }
                // Check if the cell is fully contained by the polygon
                if ((mode == CONTAINMENT_FULL ||
                     mode == CONTAINMENT_OVERLAPPING_BBOX) &&
                    cellBoundaryInsidePolygon(iter->_polygon, iter->_bboxes,
                                              &boundary, &bbox)) {
                    // Set to next output
                    iter->cell = cell;
                    return;
                }
                // For overlap, we've already checked for center point inclusion
                // above; if that failed, we only need to check for line
                // intersection
                else if ((mode == CONTAINMENT_OVERLAPPING ||
                          mode == CONTAINMENT_OVERLAPPING_BBOX) &&
                         cellBoundaryCrossesPolygon(
                             iter->_polygon, iter->_bboxes, &boundary, &bbox)) {
                    // Set to next output
                    iter->cell = cell;
                    return;
                }
            }
            if (mode == CONTAINMENT_OVERLAPPING_BBOX) {
                // Get a bounding box containing all the cell's children, so
                // this can work for the max size calculation
                BBox bbox;
                H3Error bboxErr = cellToBBox(cell, &bbox, true);
                if (bboxErr) {
                    iterErrorPolygonCompact(iter, bboxErr);
                    return;
                }
                if (bboxOverlapsBBox(&iter->_bboxes[0], &bbox)) {
                    CellBoundary bboxBoundary = bboxToCellBoundary(&bbox);
                    if (
                        // cell bbox contains the polygon
                        bboxContainsBBox(&bbox, &iter->_bboxes[0]) ||
                        // polygon contains cell bbox
                        pointInsidePolygon(iter->_polygon, iter->_bboxes,
                                           &bboxBoundary.verts[0]) ||
                        // polygon crosses cell bbox
                        cellBoundaryCrossesPolygon(iter->_polygon,
                                                   iter->_bboxes, &bboxBoundary,
                                                   &bbox)) {
                        iter->cell = cell;
                        return;
                    }
                }
            }
        }

        // Coarser cell: Check the bounding box
        if (cellRes < iter->_res) {
            // Get a bounding box for all of the cell's children
            BBox bbox;
            H3Error bboxErr = cellToBBox(cell, &bbox, true);
            if (bboxErr) {
                iterErrorPolygonCompact(iter, bboxErr);
                return;
            }
            if (bboxOverlapsBBox(&iter->_bboxes[0], &bbox)) {
                // Quick check for possible containment
                if (bboxContainsBBox(&iter->_bboxes[0], &bbox)) {
                    CellBoundary bboxBoundary = bboxToCellBoundary(&bbox);
                    // Do a fine-grained, more expensive check on the polygon
                    if (cellBoundaryInsidePolygon(iter->_polygon, iter->_bboxes,
                                                  &bboxBoundary, &bbox)) {
                        // Bounding box is fully contained, so all children are
                        // included. Set to next output.
                        iter->cell = cell;
                        return;
                    }
                }
                // Otherwise, the intersecting bbox means we need to test all
                // children, starting with the first child
                H3Index child;
                H3Error childErr =
                    H3_EXPORT(cellToCenterChild)(cell, cellRes + 1, &child);
                if (childErr) {
                    iterErrorPolygonCompact(iter, childErr);
                    return;
                }
                // Restart the loop with the child cell
                cell = child;
                continue;
            }
        }

        // Find the next cell in the sequence of all cells and continue
        cell = nextCell(cell);
    }
    // If we make it out of the loop, we're done
    iterDestroyPolygonCompact(iter);
}

/**
 * Destroy an iterator, releasing any allocated memory. Iterators destroyed in
 * this manner are safe to use but will always return H3_NULL.
 * @param  iter Iterator to destroy
 */
void iterDestroyPolygonCompact(IterCellsPolygonCompact *iter) {
    if (iter->_bboxes) {
        H3_MEMORY(free)(iter->_bboxes);
    }
    iter->cell = H3_NULL;
    iter->error = E_SUCCESS;
    iter->_polygon = NULL;
    iter->_res = -1;
    iter->_flags = 0;
    iter->_bboxes = NULL;
}

/**
 * Initialize a IterCellsPolygon struct representing the sequence of
 * cells within the target polygon. The test for including edge cells is defined
 * by the polyfill mode passed in the `flags` argument.
 *
 * Initialization of this object may fail, in which case the `error` property
 * will be set and all iteration will return H3_NULL. It is the responsibility
 * of the caller to check the error property after initialization.
 *
 * At any point in the iteration, starting once the struct is initialized, the
 * output value can be accessed through the `cell` property.
 *
 * Note that initializing the iterator allocates memory. If an iterator is
 * exhausted or returns an error that memory is released; otherwise it must be
 * released manually with iterDestroyPolygon.
 *
 * @param  polygon Polygon to fill with cells
 * @param  res     Resolution for output cells
 * @param  flags   Bit mask of option flags
 * @return         Initialized iterator, with the first value available
 */
IterCellsPolygon iterInitPolygon(const GeoPolygon *polygon, int res,
                                 uint32_t flags) {
    // Create the sub-iterator for compact cells
    IterCellsPolygonCompact cellIter =
        iterInitPolygonCompact(polygon, res, flags);
    // Create the sub-iterator for children
    IterCellsChildren childIter = iterInitParent(cellIter.cell, res);

    IterCellsPolygon iter = {.cell = childIter.h,
                             .error = cellIter.error,
                             ._cellIter = cellIter,
                             ._childIter = childIter};
    return iter;
}

/**
 * Increment the polyfill iterator, outputting the latest cell at the
 * desired resolution.
 *
 * @param  iter Iterator to increment
 */
void iterStepPolygon(IterCellsPolygon *iter) {
    if (iter->cell == H3_NULL) return;

    // See if there are more children to output
    iterStepChild(&(iter->_childIter));
    if (iter->_childIter.h) {
        iter->cell = iter->_childIter.h;
        return;
    }

    // Otherwise, increment the polyfill iterator
    iterStepPolygonCompact(&(iter->_cellIter));
    if (iter->_cellIter.cell) {
        _iterInitParent(iter->_cellIter.cell, iter->_cellIter._res,
                        &(iter->_childIter));
        iter->cell = iter->_childIter.h;
        return;
    }

    // All done, set to null and report errors if any
    iter->cell = H3_NULL;
    iter->error = iter->_cellIter.error;
}

/**
 * Destroy an iterator, releasing any allocated memory. Iterators destroyed in
 * this manner are safe to use but will always return H3_NULL.
 * @param  iter Iterator to destroy
 */
void iterDestroyPolygon(IterCellsPolygon *iter) {
    iterDestroyPolygonCompact(&(iter->_cellIter));
    // null out the child iterator by passing H3_NULL
    _iterInitParent(H3_NULL, 0, &(iter->_childIter));
    iter->cell = H3_NULL;
    iter->error = E_SUCCESS;
}

/**
 * polygonToCells takes a given GeoJSON-like data structure and preallocated,
 * zeroed memory, and fills it with the hexagons that are contained by
 * the GeoJSON-like data structure. Polygons are considered in Cartesian space.
 *
 * @param polygon The geoloop and holes defining the relevant area
 * @param res The Hexagon resolution (0-15)
 * @param flags Algorithm flags such as containment mode
 * @param size Maximum number of indexes to write to `out`.
 * @param out The slab of zeroed memory to write to. Must be at least of size
 * `size`.
 */
H3Error H3_EXPORT(polygonToCellsExperimental)(const GeoPolygon *polygon,
                                              int res, uint32_t flags,
                                              int64_t size, H3Index *out) {
    IterCellsPolygon iter = iterInitPolygon(polygon, res, flags);
    int64_t i = 0;
    for (; iter.cell; iterStepPolygon(&iter)) {
        if (i >= size) {
            iterDestroyPolygon(&iter);
            return E_MEMORY_BOUNDS;
        }
        out[i++] = iter.cell;
    }
    return iter.error;
}

static int MAX_SIZE_CELL_THRESHOLD = 10;

static double getAverageCellArea(int res) {
    double hexAreaKm2;
    H3_EXPORT(getHexagonAreaAvgKm2)(res, &hexAreaKm2);
    return hexAreaKm2;
}

/**
 * maxPolygonToCellsSize returns the number of cells to allocate space for
 * when performing a polygonToCells on the given GeoJSON-like data structure.
 * @param geoPolygon A GeoJSON-like data structure indicating the poly to fill
 * @param res Hexagon resolution (0-15)
 * @param out number of cells to allocate for
 * @return 0 (E_SUCCESS) on success.
 */
H3Error H3_EXPORT(maxPolygonToCellsSizeExperimental)(const GeoPolygon *polygon,
                                                     int res, uint32_t flags,
                                                     int64_t *out) {
    // Special case: 0-vertex polygon
    if (polygon->geoloop.numVerts == 0) {
        *out = 0;
        return E_SUCCESS;
    }

    // Initialize the iterator without stepping, so we can adjust the res and
    // flags (after they are validated by the initialization) before we start
    IterCellsPolygonCompact iter = _iterInitPolygonCompact(polygon, res, flags);

    if (iter.error) {
        return iter.error;
    }

    // Ignore the requested flags and use the faster overlapping-bbox mode
    iter._flags = CONTAINMENT_OVERLAPPING_BBOX;

    // Get a (very) rough area of the polygon bounding box
    BBox *polygonBBox = &iter._bboxes[0];
    double polygonBBoxAreaKm2 =
        bboxHeightRads(polygonBBox) * bboxWidthRads(polygonBBox) /
        cos(fmin(fabs(polygonBBox->north), fabs(polygonBBox->south))) *
        EARTH_RADIUS_KM * EARTH_RADIUS_KM;

    // Determine the res for the size estimate, based on a (very) rough estimate
    // of the number of cells at various resolutions that would fit in the
    // polygon. All we need here is a general order of magnitude.
    while (iter._res > 0 &&
           polygonBBoxAreaKm2 / getAverageCellArea(iter._res - 1) >
               MAX_SIZE_CELL_THRESHOLD) {
        iter._res--;
    }

    // Now run the polyfill, counting the output in the target res.
    // We have to take the first step outside the loop, to get the first
    // valid output cell
    iterStepPolygonCompact(&iter);

    *out = 0;
    int64_t childrenSize;
    for (; iter.cell; iterStepPolygonCompact(&iter)) {
        H3_EXPORT(cellToChildrenSize)(iter.cell, res, &childrenSize);
        *out += childrenSize;
    }

    return iter.error;
}