diff --git a/scripts/make_countries.js b/scripts/make_countries.js
index 8e810dff3..8741de46e 100644
--- a/scripts/make_countries.js
+++ b/scripts/make_countries.js
@@ -140,9 +140,8 @@ const GeoPolygon COUNTRIES[${polygons.length}] = {${
 BEGIN_BENCHMARKS();
 
 int MAX_RES = 5;
-// Max size of largest polygon at res 5, rounded up with padding
-int64_t numHexagons = 400000;
-H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+int64_t numHexagons;
+H3Index *hexagons;
 
 for (int res = 0; res < MAX_RES + 1; res++) {
 
@@ -150,31 +149,41 @@ for (int res = 0; res < MAX_RES + 1; res++) {
 
   BENCHMARK(polygonToCells_AllCountries1, 5, {
     for (int index = 0; index < ${polygons.length}; index++) {
+      H3_EXPORT(maxPolygonToCellsSize)(&COUNTRIES[index], res, CONTAINMENT_CENTER, &numHexagons);
+      hexagons = calloc(numHexagons, sizeof(H3Index));
       H3_EXPORT(polygonToCells)(&COUNTRIES[index], res, CONTAINMENT_CENTER, hexagons);
+      free(hexagons);
     }
   });
 
   BENCHMARK(polygonToCells_AllCountries2, 5, {
     for (int index = 0; index < ${polygons.length}; index++) {
+      H3_EXPORT(maxPolygonToCellsSizeExperimental)(&COUNTRIES[index], res, CONTAINMENT_CENTER, &numHexagons);
+      hexagons = calloc(numHexagons, sizeof(H3Index));
       H3_EXPORT(polygonToCellsExperimental)(&COUNTRIES[index], res, CONTAINMENT_CENTER, hexagons);
+      free(hexagons);
     }
   });
 
   BENCHMARK(polygonToCells_AllCountries3, 5, {
     for (int index = 0; index < ${polygons.length}; index++) {
+      H3_EXPORT(maxPolygonToCellsSizeExperimental)(&COUNTRIES[index], res, CONTAINMENT_CENTER, &numHexagons);
+      hexagons = calloc(numHexagons, sizeof(H3Index));
       H3_EXPORT(polygonToCellsExperimental)(&COUNTRIES[index], res, CONTAINMENT_FULL, hexagons);
     }
   });
 
   BENCHMARK(polygonToCells_AllCountries4, 5, {
     for (int index = 0; index < ${polygons.length}; index++) {
+      H3_EXPORT(maxPolygonToCellsSizeExperimental)(&COUNTRIES[index], res, CONTAINMENT_CENTER, &numHexagons);
+      hexagons = calloc(numHexagons, sizeof(H3Index));
       H3_EXPORT(polygonToCellsExperimental)(&COUNTRIES[index], res, CONTAINMENT_OVERLAPPING, hexagons);
+      free(hexagons);
     }
   });
 
 }
 
-free(hexagons);
 
 END_BENCHMARKS();
   `
diff --git a/src/apps/testapps/testPolyfillInternal.c b/src/apps/testapps/testPolyfillInternal.c
index 4e44255ad..7efeebf79 100644
--- a/src/apps/testapps/testPolyfillInternal.c
+++ b/src/apps/testapps/testPolyfillInternal.c
@@ -68,10 +68,10 @@ SUITE(polyfillInternal) {
         IterCellsPolygonCompact iter;
         H3Index cell;
 
+        // Give the iterator a cell with a bad base cell
         iter = iterInitPolygonCompact(&sfGeoPolygon, 9, CONTAINMENT_CENTER);
         t_assertSuccess(iter.error);
 
-        // Give the iterator a cell with a bad base cell
         cell = 0x85283473fffffff;
         H3_SET_BASE_CELL(cell, 123);
         iter.cell = cell;
@@ -81,10 +81,10 @@ SUITE(polyfillInternal) {
                  "Got expected error for invalid cell");
         t_assert(iter.cell == H3_NULL, "Got null output for invalid cell");
 
+        // Give the iterator a cell with a bad base cell, at the target res
         iter = iterInitPolygonCompact(&sfGeoPolygon, 9, CONTAINMENT_CENTER);
         t_assertSuccess(iter.error);
 
-        // Give the iterator a cell with a bad base cell, at the target res
         cell = 0x89283470003ffff;
         H3_SET_BASE_CELL(cell, 123);
         iter.cell = cell;
@@ -95,11 +95,11 @@ SUITE(polyfillInternal) {
         t_assert(iter.cell == H3_NULL,
                  "Got null output for invalid cell at res");
 
+        // Give the iterator a cell with a bad base cell, at the target res
+        // (full containment)
         iter = iterInitPolygonCompact(&sfGeoPolygon, 9, CONTAINMENT_FULL);
         t_assertSuccess(iter.error);
 
-        // Give the iterator a cell with a bad base cell, at the target res
-        // (full containment)
         cell = 0x89283470003ffff;
         H3_SET_BASE_CELL(cell, 123);
         iter.cell = cell;
@@ -110,12 +110,28 @@ SUITE(polyfillInternal) {
         t_assert(iter.cell == H3_NULL,
                  "Got null output for invalid cell at res");
 
-        iter = iterInitPolygonCompact(&sfGeoPolygon, 9, CONTAINMENT_CENTER);
+        // Give the iterator a cell with a bad base cell, at the target res
+        // (overlapping bounding box)
+        iter = iterInitPolygonCompact(&sfGeoPolygon, 9,
+                                      CONTAINMENT_OVERLAPPING_BBOX);
         t_assertSuccess(iter.error);
 
+        cell = 0x89283470003ffff;
+        H3_SET_BASE_CELL(cell, 123);
+        iter.cell = cell;
+
+        iterStepPolygonCompact(&iter);
+        t_assert(iter.error == E_CELL_INVALID,
+                 "Got expected error for invalid cell");
+        t_assert(iter.cell == H3_NULL,
+                 "Got null output for invalid cell at res");
+
         // Give the iterator a cell that's too fine for a child check,
         // and a target resolution that allows this to run. This cell has
         // to be inside the polygon to reach the error.
+        iter = iterInitPolygonCompact(&sfGeoPolygon, 9, CONTAINMENT_CENTER);
+        t_assertSuccess(iter.error);
+
         cell = 0x8f283080dcb019a;
         iter.cell = cell;
         iter._res = 42;
diff --git a/src/apps/testapps/testPolygonToCells.c b/src/apps/testapps/testPolygonToCells.c
index b35dedd58..f8724bac9 100644
--- a/src/apps/testapps/testPolygonToCells.c
+++ b/src/apps/testapps/testPolygonToCells.c
@@ -22,6 +22,7 @@
 #include "constants.h"
 #include "h3Index.h"
 #include "latLng.h"
+#include "polygon.h"
 #include "test.h"
 #include "utility.h"
 
@@ -461,7 +462,7 @@ SUITE(polygonToCells) {
 
     TEST(invalidFlags) {
         int64_t numHexagons;
-        for (uint32_t flags = 3; flags <= 32; flags++) {
+        for (uint32_t flags = CONTAINMENT_INVALID; flags <= 32; flags++) {
             t_assert(
                 H3_EXPORT(maxPolygonToCellsSize)(
                     &sfGeoPolygon, 9, flags, &numHexagons) == E_OPTION_INVALID,
@@ -471,7 +472,7 @@ SUITE(polygonToCells) {
         t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0,
                                                          &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
-        for (uint32_t flags = 3; flags <= 32; flags++) {
+        for (uint32_t flags = CONTAINMENT_INVALID; flags <= 32; flags++) {
             t_assert(H3_EXPORT(polygonToCells)(&sfGeoPolygon, 9, flags,
                                                hexagons) == E_OPTION_INVALID,
                      "Flags other than polyfill modes are invalid for "
diff --git a/src/apps/testapps/testPolygonToCellsExperimental.c b/src/apps/testapps/testPolygonToCellsExperimental.c
index 24464597d..045bf7261 100644
--- a/src/apps/testapps/testPolygonToCellsExperimental.c
+++ b/src/apps/testapps/testPolygonToCellsExperimental.c
@@ -95,8 +95,8 @@ static void fillIndex_assertions(H3Index h) {
             .holes = 0};
 
         int64_t polygonToCellsSize;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&polygon, nextRes, 0,
-                                                         &polygonToCellsSize));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
+            &polygon, nextRes, 0, &polygonToCellsSize));
         H3Index *polygonToCellsOut =
             calloc(polygonToCellsSize, sizeof(H3Index));
         t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
@@ -161,7 +161,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCells) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &sfGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -175,7 +175,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCells_FullContainment) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &sfGeoPolygon, 9, CONTAINMENT_FULL, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -190,7 +190,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCells_Overlapping) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &sfGeoPolygon, 9, CONTAINMENT_OVERLAPPING, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -205,7 +205,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCellsHole) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &holeGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -220,7 +220,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCellsHoleFullContainment) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &holeGeoPolygon, 9, CONTAINMENT_FULL, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -236,7 +236,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCellsHoleOverlapping) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &holeGeoPolygon, 9, CONTAINMENT_OVERLAPPING, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -251,7 +251,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCellsEmpty) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &emptyGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -266,7 +266,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCellsContainsPolygon) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &sfGeoPolygon, 4, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -290,7 +290,7 @@ SUITE(polygonToCells) {
         centerGeoPolygon.numHoles = 0;
 
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &centerGeoPolygon, 4, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -306,7 +306,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCellsContainsPolygon_FullContainment) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &sfGeoPolygon, 4, CONTAINMENT_FULL, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -321,7 +321,7 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCellsContainsPolygon_Overlapping) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &sfGeoPolygon, 4, CONTAINMENT_OVERLAPPING, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -356,7 +356,7 @@ SUITE(polygonToCells) {
         someHexagon.numHoles = 0;
 
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &someHexagon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -416,7 +416,7 @@ SUITE(polygonToCells) {
         // Prime meridian case
         expectedSize = 4228;
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &primeMeridianGeoPolygon, 7, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -431,7 +431,7 @@ SUITE(polygonToCells) {
         // This doesn't exactly match the prime meridian count because of slight
         // differences in hex size and grid offset between the two cases
         expectedSize = 4238;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &transMeridianGeoPolygon, 7, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagonsTM = calloc(numHexagons, sizeof(H3Index));
 
@@ -444,7 +444,7 @@ SUITE(polygonToCells) {
 
         // Transmeridian filled hole case -- only needed for calculating hole
         // size
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &transMeridianFilledHoleGeoPolygon, 7, CONTAINMENT_CENTER,
             &numHexagons));
         H3Index *hexagonsTMFH = calloc(numHexagons, sizeof(H3Index));
@@ -456,7 +456,7 @@ SUITE(polygonToCells) {
             countNonNullIndexes(hexagonsTMFH, numHexagons);
 
         // Transmeridian hole case
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &transMeridianHoleGeoPolygon, 7, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagonsTMH = calloc(numHexagons, sizeof(H3Index));
 
@@ -484,7 +484,7 @@ SUITE(polygonToCells) {
         GeoPolygon polygon = {.geoloop = geoloop, .numHoles = 0};
 
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &polygon, 4, CONTAINMENT_CENTER, &numHexagons));
 
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
@@ -536,7 +536,7 @@ SUITE(polygonToCells) {
         polygon.numHoles = 0;
 
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &polygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
@@ -560,21 +560,21 @@ SUITE(polygonToCells) {
 
     TEST(invalidFlags) {
         int64_t numHexagons;
-        for (uint32_t flags = 3; flags <= 32; flags++) {
+        for (uint32_t flags = CONTAINMENT_INVALID; flags <= 32; flags++) {
             t_assert(
-                H3_EXPORT(maxPolygonToCellsSize)(
+                H3_EXPORT(maxPolygonToCellsSizeExperimental)(
                     &sfGeoPolygon, 9, flags, &numHexagons) == E_OPTION_INVALID,
                 "Flags other than polyfill modes are invalid for "
-                "maxPolygonToCellsSize");
+                "maxPolygonToCellsSizeExperimental");
         }
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSizeExperimental)(
             &sfGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
-        for (uint32_t flags = 3; flags <= 32; flags++) {
+        for (uint32_t flags = CONTAINMENT_INVALID; flags <= 32; flags++) {
             t_assert(H3_EXPORT(polygonToCellsExperimental)(
                          &sfGeoPolygon, 9, flags, hexagons) == E_OPTION_INVALID,
                      "Flags other than polyfill modes are invalid for "
-                     "polygonToCells");
+                     "polygonToCellsExperimental");
         }
         free(hexagons);
     }
diff --git a/src/h3lib/include/bbox.h b/src/h3lib/include/bbox.h
index 3ae68cbf8..1c56982f4 100644
--- a/src/h3lib/include/bbox.h
+++ b/src/h3lib/include/bbox.h
@@ -43,6 +43,7 @@ bool bboxContains(const BBox *bbox, const LatLng *point);
 bool bboxContainsBBox(const BBox *a, const BBox *b);
 bool bboxOverlapsBBox(const BBox *a, const BBox *b);
 bool bboxEquals(const BBox *b1, const BBox *b2);
+CellBoundary bboxToCellBoundary(const BBox *bbox);
 H3Error bboxHexEstimate(const BBox *bbox, int res, int64_t *out);
 H3Error lineHexEstimate(const LatLng *origin, const LatLng *destination,
                         int res, int64_t *out);
diff --git a/src/h3lib/include/polyfill.h b/src/h3lib/include/polyfill.h
index 3bed9bc67..02e0632d9 100644
--- a/src/h3lib/include/polyfill.h
+++ b/src/h3lib/include/polyfill.h
@@ -81,6 +81,11 @@ DECLSPEC void iterDestroyPolygon(IterCellsPolygon *iter);
 DECLSPEC H3Error H3_EXPORT(polygonToCellsExperimental)(
     const GeoPolygon *polygon, int res, uint32_t flags, H3Index *out);
 
+DECLSPEC H3Error H3_EXPORT(polygonToCellsExperimental)(
+    const GeoPolygon *polygon, int res, uint32_t flags, H3Index *out);
+DECLSPEC H3Error H3_EXPORT(maxPolygonToCellsSizeExperimental)(
+    const GeoPolygon *polygon, int res, uint32_t flags, int64_t *out);
+
 H3Error cellToBBox(H3Index cell, BBox *out, bool coverChildren);
 H3Index baseCellNumToCell(int baseCellNum);
 
diff --git a/src/h3lib/include/polygon.h b/src/h3lib/include/polygon.h
index d7d24ffbc..b6fb37889 100644
--- a/src/h3lib/include/polygon.h
+++ b/src/h3lib/include/polygon.h
@@ -48,7 +48,8 @@ typedef enum {
     CONTAINMENT_CENTER = 0,       ///< Cell center is contained in the shape
     CONTAINMENT_FULL = 1,         ///< Cell is fully contained in the shape
     CONTAINMENT_OVERLAPPING = 2,  ///< Cell overlaps the shape at any point
-    CONTAINMENT_INVALID = 3  ///< This mode is invalid and should not be used
+    CONTAINMENT_OVERLAPPING_BBOX = 3,  ///< Cell bounding box overlaps shape
+    CONTAINMENT_INVALID = 4  ///< This mode is invalid and should not be used
 } ContainmentMode;
 
 // 1s in the 4 bits defining the polyfill containment mode, 0s elsewhere
diff --git a/src/h3lib/lib/bbox.c b/src/h3lib/lib/bbox.c
index b86d88000..70d56ffe2 100644
--- a/src/h3lib/lib/bbox.c
+++ b/src/h3lib/lib/bbox.c
@@ -140,6 +140,16 @@ bool bboxEquals(const BBox *b1, const BBox *b2) {
            b1->east == b2->east && b1->west == b2->west;
 }
 
+CellBoundary bboxToCellBoundary(const BBox *bbox) {
+    // Convert bbox to cell boundary, CCW vertex order
+    CellBoundary bboxBoundary = {.numVerts = 4,
+                                 .verts = {{bbox->north, bbox->east},
+                                           {bbox->north, bbox->west},
+                                           {bbox->south, bbox->west},
+                                           {bbox->south, bbox->east}}};
+    return bboxBoundary;
+}
+
 /**
  * _hexRadiusKm returns the radius of a given hexagon in Km
  *
diff --git a/src/h3lib/lib/polyfill.c b/src/h3lib/lib/polyfill.c
index bde5eedbc..688cb1529 100644
--- a/src/h3lib/lib/polyfill.c
+++ b/src/h3lib/lib/polyfill.c
@@ -316,28 +316,11 @@ static H3Index nextCell(H3Index cell) {
 }
 
 /**
- * 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
+ * Internal function - initialize the iterator without stepping to the first
+ * value
  */
-IterCellsPolygonCompact iterInitPolygonCompact(const GeoPolygon *polygon,
-                                               int res, uint32_t flags) {
+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),
@@ -369,6 +352,34 @@ IterCellsPolygonCompact iterInitPolygonCompact(const GeoPolygon *polygon,
     }
     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);
@@ -482,6 +493,32 @@ void iterStepPolygonCompact(IterCellsPolygonCompact *iter) {
                     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
@@ -496,13 +533,7 @@ void iterStepPolygonCompact(IterCellsPolygonCompact *iter) {
             if (bboxOverlapsBBox(&iter->_bboxes[0], &bbox)) {
                 // Quick check for possible containment
                 if (bboxContainsBBox(&iter->_bboxes[0], &bbox)) {
-                    // Convert bbox to cell boundary, CCW vertex order
-                    CellBoundary bboxBoundary = {
-                        .numVerts = 4,
-                        .verts = {{bbox.north, bbox.east},
-                                  {bbox.north, bbox.west},
-                                  {bbox.south, bbox.west},
-                                  {bbox.south, bbox.east}}};
+                    CellBoundary bboxBoundary = bboxToCellBoundary(&bbox);
                     // Do a fine-grained, more expensive check on the polygon
                     if (cellBoundaryInsidePolygon(iter->_polygon, iter->_bboxes,
                                                   &bboxBoundary, &bbox)) {
@@ -647,5 +678,66 @@ H3Error H3_EXPORT(polygonToCellsExperimental)(const GeoPolygon *polygon,
     for (; iter.cell; iterStepPolygon(&iter)) {
         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) {
+    // 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;
 }
\ No newline at end of file