diff --git a/scripts/make_countries.js b/scripts/make_countries.js
index 4e8a495d4..8e810dff3 100644
--- a/scripts/make_countries.js
+++ b/scripts/make_countries.js
@@ -131,6 +131,7 @@ async function makeCountries(sourceUrl, targetPath) {
 #include "h3api.h"
 #include "mathExtensions.h"
 #include "polyfill.h"
+#include "polygon.h"
 
 const GeoPolygon COUNTRIES[${polygons.length}] = {${
   polygons.map((poly, i) => formatGeoPolygon(poly, names[i])).join(',')
@@ -149,13 +150,25 @@ for (int res = 0; res < MAX_RES + 1; res++) {
 
   BENCHMARK(polygonToCells_AllCountries1, 5, {
     for (int index = 0; index < ${polygons.length}; index++) {
-      H3_EXPORT(polygonToCells)(&COUNTRIES[index], res, 0, hexagons);
+      H3_EXPORT(polygonToCells)(&COUNTRIES[index], res, CONTAINMENT_CENTER, hexagons);
     }
   });
 
   BENCHMARK(polygonToCells_AllCountries2, 5, {
     for (int index = 0; index < ${polygons.length}; index++) {
-      H3_EXPORT(polygonToCellsExperimental)(&COUNTRIES[index], res, 0, hexagons);
+      H3_EXPORT(polygonToCellsExperimental)(&COUNTRIES[index], res, CONTAINMENT_CENTER, hexagons);
+    }
+  });
+
+  BENCHMARK(polygonToCells_AllCountries3, 5, {
+    for (int index = 0; index < ${polygons.length}; index++) {
+      H3_EXPORT(polygonToCellsExperimental)(&COUNTRIES[index], res, CONTAINMENT_FULL, hexagons);
+    }
+  });
+
+  BENCHMARK(polygonToCells_AllCountries4, 5, {
+    for (int index = 0; index < ${polygons.length}; index++) {
+      H3_EXPORT(polygonToCellsExperimental)(&COUNTRIES[index], res, CONTAINMENT_OVERLAPPING, hexagons);
     }
   });
 
diff --git a/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c b/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c
index b0837c7fc..1034f3372 100644
--- a/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c
+++ b/src/apps/benchmarks/benchmarkPolygonToCellsExperimental.c
@@ -17,6 +17,7 @@
 #include "benchmark.h"
 #include "h3api.h"
 #include "polyfill.h"
+#include "polygon.h"
 
 // Fixtures
 LatLng sfVerts[] = {
@@ -122,24 +123,84 @@ southernGeoPolygon.geoloop = southernGeoLoop;
 int64_t numHexagons;
 H3Index *hexagons;
 
-BENCHMARK(polygonToCellsSF, 500, {
-    H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0, &numHexagons);
+BENCHMARK(polygonToCellsSF_Center, 500, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&sfGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons);
     hexagons = calloc(numHexagons, sizeof(H3Index));
-    H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, 0, hexagons);
+    H3_EXPORT(polygonToCellsExperimental)
+    (&sfGeoPolygon, 9, CONTAINMENT_CENTER, hexagons);
     free(hexagons);
 });
 
-BENCHMARK(polygonToCellsAlameda, 500, {
-    H3_EXPORT(maxPolygonToCellsSize)(&alamedaGeoPolygon, 9, 0, &numHexagons);
+BENCHMARK(polygonToCellsSF_Full, 500, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&sfGeoPolygon, 9, CONTAINMENT_FULL, &numHexagons);
     hexagons = calloc(numHexagons, sizeof(H3Index));
-    H3_EXPORT(polygonToCellsExperimental)(&alamedaGeoPolygon, 9, 0, hexagons);
+    H3_EXPORT(polygonToCellsExperimental)
+    (&sfGeoPolygon, 9, CONTAINMENT_FULL, hexagons);
     free(hexagons);
 });
 
-BENCHMARK(polygonToCellsSouthernExpansion, 10, {
-    H3_EXPORT(maxPolygonToCellsSize)(&southernGeoPolygon, 9, 0, &numHexagons);
+BENCHMARK(polygonToCellsSF_Overlapping, 500, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&sfGeoPolygon, 9, CONTAINMENT_OVERLAPPING, &numHexagons);
     hexagons = calloc(numHexagons, sizeof(H3Index));
-    H3_EXPORT(polygonToCellsExperimental)(&southernGeoPolygon, 9, 0, hexagons);
+    H3_EXPORT(polygonToCellsExperimental)
+    (&sfGeoPolygon, 9, CONTAINMENT_OVERLAPPING, hexagons);
+    free(hexagons);
+});
+
+BENCHMARK(polygonToCellsAlameda_Center, 500, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&alamedaGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons);
+    hexagons = calloc(numHexagons, sizeof(H3Index));
+    H3_EXPORT(polygonToCellsExperimental)
+    (&alamedaGeoPolygon, 9, CONTAINMENT_CENTER, hexagons);
+    free(hexagons);
+});
+
+BENCHMARK(polygonToCellsAlameda_Full, 500, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&alamedaGeoPolygon, 9, CONTAINMENT_FULL, &numHexagons);
+    hexagons = calloc(numHexagons, sizeof(H3Index));
+    H3_EXPORT(polygonToCellsExperimental)
+    (&alamedaGeoPolygon, 9, CONTAINMENT_FULL, hexagons);
+    free(hexagons);
+});
+
+BENCHMARK(polygonToCellsAlameda_Overlapping, 500, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&alamedaGeoPolygon, 9, CONTAINMENT_OVERLAPPING, &numHexagons);
+    hexagons = calloc(numHexagons, sizeof(H3Index));
+    H3_EXPORT(polygonToCellsExperimental)
+    (&alamedaGeoPolygon, 9, CONTAINMENT_OVERLAPPING, hexagons);
+    free(hexagons);
+});
+
+BENCHMARK(polygonToCellsSouthernExpansion_Center, 10, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&southernGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons);
+    hexagons = calloc(numHexagons, sizeof(H3Index));
+    H3_EXPORT(polygonToCellsExperimental)
+    (&southernGeoPolygon, 9, CONTAINMENT_CENTER, hexagons);
+    free(hexagons);
+});
+
+BENCHMARK(polygonToCellsSouthernExpansion_Full, 10, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&southernGeoPolygon, 9, CONTAINMENT_FULL, &numHexagons);
+    hexagons = calloc(numHexagons, sizeof(H3Index));
+    H3_EXPORT(polygonToCellsExperimental)
+    (&southernGeoPolygon, 9, CONTAINMENT_FULL, hexagons);
+    free(hexagons);
+});
+
+BENCHMARK(polygonToCellsSouthernExpansion_Overlapping, 10, {
+    H3_EXPORT(maxPolygonToCellsSize)
+    (&southernGeoPolygon, 9, CONTAINMENT_OVERLAPPING, &numHexagons);
+    hexagons = calloc(numHexagons, sizeof(H3Index));
+    H3_EXPORT(polygonToCellsExperimental)
+    (&southernGeoPolygon, 9, CONTAINMENT_OVERLAPPING, hexagons);
     free(hexagons);
 });
 
diff --git a/src/apps/testapps/testH3Memory.c b/src/apps/testapps/testH3Memory.c
index 96627b543..54a7c688e 100644
--- a/src/apps/testapps/testH3Memory.c
+++ b/src/apps/testapps/testH3Memory.c
@@ -26,6 +26,7 @@
 #include "h3api.h"
 #include "latLng.h"
 #include "polyfill.h"
+#include "polygon.h"
 #include "test.h"
 #include "utility.h"
 
@@ -230,22 +231,22 @@ SUITE(h3Memory) {
         sfGeoPolygon.numHoles = 0;
 
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &sfGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
         resetMemoryCounters(0);
         failAlloc = true;
-        H3Error err = H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, 0,
-                                                            hexagons);
+        H3Error err = H3_EXPORT(polygonToCellsExperimental)(
+            &sfGeoPolygon, 9, CONTAINMENT_CENTER, hexagons);
         t_assert(err == E_MEMORY_ALLOC,
                  "polygonToCellsExperimental failed (1)");
         t_assert(actualAllocCalls == 1, "alloc called once");
         t_assert(actualFreeCalls == 0, "free not called");
 
         resetMemoryCounters(1);
-        err = H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9, 0,
-                                                    hexagons);
+        err = H3_EXPORT(polygonToCellsExperimental)(
+            &sfGeoPolygon, 9, CONTAINMENT_CENTER, hexagons);
         t_assert(err == E_SUCCESS, "polygonToCellsExperimental succeeded (1)");
         t_assert(actualAllocCalls == 1, "alloc called one time");
         t_assert(actualFreeCalls == 1, "free called one time");
diff --git a/src/apps/testapps/testPolyfillInternal.c b/src/apps/testapps/testPolyfillInternal.c
index db393cc30..4e44255ad 100644
--- a/src/apps/testapps/testPolyfillInternal.c
+++ b/src/apps/testapps/testPolyfillInternal.c
@@ -14,11 +14,14 @@
  * limitations under the License.
  */
 
+#include <math.h>
+
 #include "bbox.h"
 #include "h3Index.h"
 #include "h3api.h"
 #include "latLng.h"
 #include "polyfill.h"
+#include "polygon.h"
 #include "test.h"
 #include "utility.h"
 
@@ -33,16 +36,24 @@ static GeoPolygon sfGeoPolygon = {
                                     {0.6599990002976, -2.1376771158464}}},
     .numHoles = 0};
 
+static GeoPolygon invalidGeoPolygon = {
+    .geoloop = {.numVerts = 4,
+                .verts = (LatLng[]){{NAN, -2.1364398519396},
+                                    {0.6595011102219, NAN},
+                                    {NAN, -2.1354884206045},
+                                    {0.6581220034068, NAN}}},
+    .numHoles = 0};
+
 SUITE(polyfillInternal) {
     TEST(iterInitPolygonCompact_errors) {
         IterCellsPolygonCompact iter;
 
-        iter = iterInitPolygonCompact(&sfGeoPolygon, -1, 0);
+        iter = iterInitPolygonCompact(&sfGeoPolygon, -1, CONTAINMENT_CENTER);
         t_assert(iter.error == E_RES_DOMAIN,
                  "Got expected error for invalid res");
         t_assert(iter.cell == H3_NULL, "Got null output for invalid res");
 
-        iter = iterInitPolygonCompact(&sfGeoPolygon, 16, 0);
+        iter = iterInitPolygonCompact(&sfGeoPolygon, 16, CONTAINMENT_CENTER);
         t_assert(iter.error == E_RES_DOMAIN,
                  "Got expected error for invalid res");
         t_assert(iter.cell == H3_NULL, "Got null output for invalid res");
@@ -53,11 +64,12 @@ SUITE(polyfillInternal) {
         t_assert(iter.cell == H3_NULL, "Got null output for invalid flags");
     }
 
-    TEST(iterStepPolygonCompact_errors) {
+    TEST(iterStepPolygonCompact_invalidCellErrors) {
         IterCellsPolygonCompact iter;
         H3Index cell;
 
-        iter = iterInitPolygonCompact(&sfGeoPolygon, 9, 0);
+        iter = iterInitPolygonCompact(&sfGeoPolygon, 9, CONTAINMENT_CENTER);
+        t_assertSuccess(iter.error);
 
         // Give the iterator a cell with a bad base cell
         cell = 0x85283473fffffff;
@@ -69,7 +81,37 @@ SUITE(polyfillInternal) {
                  "Got expected error for invalid cell");
         t_assert(iter.cell == H3_NULL, "Got null output for invalid cell");
 
-        iter = iterInitPolygonCompact(&sfGeoPolygon, 9, 0);
+        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;
+
+        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");
+
+        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;
+
+        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");
+
+        iter = iterInitPolygonCompact(&sfGeoPolygon, 9, CONTAINMENT_CENTER);
+        t_assertSuccess(iter.error);
 
         // 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
@@ -84,9 +126,28 @@ SUITE(polyfillInternal) {
         t_assert(iter.cell == H3_NULL, "Got null output for invalid cell");
     }
 
+    TEST(iterStepPolygonCompact_invalidPolygonErrors) {
+        IterCellsPolygonCompact iter;
+
+        // Start with a good polygon, otherwise we error out early
+        iter =
+            iterInitPolygonCompact(&sfGeoPolygon, 5, CONTAINMENT_OVERLAPPING);
+        t_assertSuccess(iter.error);
+
+        // Give the iterator a bad polygon and a cell at target res
+        iter._polygon = &invalidGeoPolygon;
+        iter.cell = 0x85283473fffffff;
+
+        iterStepPolygonCompact(&iter);
+        t_assert(iter.error == E_LATLNG_DOMAIN,
+                 "Got expected error for invalid polygon");
+        t_assert(iter.cell == H3_NULL, "Got null output for invalid cell");
+    }
+
     TEST(iterDestroyPolygonCompact) {
         IterCellsPolygonCompact iter =
-            iterInitPolygonCompact(&sfGeoPolygon, 9, 0);
+            iterInitPolygonCompact(&sfGeoPolygon, 9, CONTAINMENT_CENTER);
+        t_assertSuccess(iter.error);
 
         iterDestroyPolygonCompact(&iter);
         t_assert(iter.error == E_SUCCESS, "Got success for destroyed iterator");
@@ -101,7 +162,9 @@ SUITE(polyfillInternal) {
     }
 
     TEST(iterDestroyPolygon) {
-        IterCellsPolygon iter = iterInitPolygon(&sfGeoPolygon, 9, 0);
+        IterCellsPolygon iter =
+            iterInitPolygon(&sfGeoPolygon, 9, CONTAINMENT_CENTER);
+        t_assertSuccess(iter.error);
 
         iterDestroyPolygon(&iter);
         t_assert(iter.error == E_SUCCESS, "Got success for destroyed iterator");
diff --git a/src/apps/testapps/testPolygonInternal.c b/src/apps/testapps/testPolygonInternal.c
index 4b0f80286..42d19d14a 100644
--- a/src/apps/testapps/testPolygonInternal.c
+++ b/src/apps/testapps/testPolygonInternal.c
@@ -626,55 +626,55 @@ SUITE(polygonInternal) {
         H3_EXPORT(destroyLinkedMultiPolygon)(&polygon);
     }
 
-    TEST(lineIntersectsLine) {
+    TEST(lineCrossesLine) {
         LatLng lines1[4] = {{0, 0}, {1, 1}, {0, 1}, {1, 0}};
         t_assert(
-            lineIntersectsLine(&lines1[0], &lines1[1], &lines1[2], &lines1[3]),
+            lineCrossesLine(&lines1[0], &lines1[1], &lines1[2], &lines1[3]),
             "diagonal intersection");
 
         LatLng lines2[4] = {{1, 1}, {0, 0}, {1, 0}, {0, 1}};
         t_assert(
-            lineIntersectsLine(&lines2[0], &lines2[1], &lines2[2], &lines2[3]),
+            lineCrossesLine(&lines2[0], &lines2[1], &lines2[2], &lines2[3]),
             "diagonal intersection, reverse vertexes");
 
         LatLng lines3[4] = {{0.5, 0}, {0.5, 1}, {0, 0.5}, {1, 0.5}};
         t_assert(
-            lineIntersectsLine(&lines3[0], &lines3[1], &lines3[2], &lines3[3]),
+            lineCrossesLine(&lines3[0], &lines3[1], &lines3[2], &lines3[3]),
             "horizontal/vertical intersection");
 
         LatLng lines4[4] = {{0.5, 1}, {0.5, 0}, {1, 0.5}, {0, 0.5}};
         t_assert(
-            lineIntersectsLine(&lines4[0], &lines4[1], &lines4[2], &lines4[3]),
+            lineCrossesLine(&lines4[0], &lines4[1], &lines4[2], &lines4[3]),
             "horizontal/vertical intersection, reverse vertexes");
 
         LatLng lines5[4] = {{0, 0}, {0.4, 0.4}, {0, 1}, {1, 0}};
         t_assert(
-            !lineIntersectsLine(&lines5[0], &lines5[1], &lines5[2], &lines5[3]),
+            !lineCrossesLine(&lines5[0], &lines5[1], &lines5[2], &lines5[3]),
             "diagonal non-intersection, below");
 
         LatLng lines6[4] = {{0.6, 0.6}, {1, 1}, {0, 1}, {1, 0}};
         t_assert(
-            !lineIntersectsLine(&lines6[0], &lines6[1], &lines6[2], &lines6[3]),
+            !lineCrossesLine(&lines6[0], &lines6[1], &lines6[2], &lines6[3]),
             "diagonal non-intersection, above");
 
         LatLng lines7[4] = {{0.5, 0}, {0.5, 1}, {0, 0.5}, {0.4, 0.5}};
         t_assert(
-            !lineIntersectsLine(&lines7[0], &lines7[1], &lines7[2], &lines7[3]),
+            !lineCrossesLine(&lines7[0], &lines7[1], &lines7[2], &lines7[3]),
             "horizontal/vertical non-intersection, below");
 
         LatLng lines8[4] = {{0.5, 0}, {0.5, 1}, {0.6, 0.5}, {1, 0.5}};
         t_assert(
-            !lineIntersectsLine(&lines8[0], &lines8[1], &lines8[2], &lines8[3]),
+            !lineCrossesLine(&lines8[0], &lines8[1], &lines8[2], &lines8[3]),
             "horizontal/vertical non-intersection, above");
 
         LatLng lines9[4] = {{0.5, 0}, {0.5, 0.4}, {0, 0.5}, {1, 0.5}};
         t_assert(
-            !lineIntersectsLine(&lines9[0], &lines9[1], &lines9[2], &lines9[3]),
+            !lineCrossesLine(&lines9[0], &lines9[1], &lines9[2], &lines9[3]),
             "horizontal/vertical non-intersection, left");
 
         LatLng lines10[4] = {{0.5, 0.6}, {0.5, 1}, {0, 0.5}, {1, 0.5}};
-        t_assert(!lineIntersectsLine(&lines10[0], &lines10[1], &lines10[2],
-                                     &lines10[3]),
+        t_assert(!lineCrossesLine(&lines10[0], &lines10[1], &lines10[2],
+                                  &lines10[3]),
                  "horizontal/vertical non-intersection, right");
     }
 
diff --git a/src/apps/testapps/testPolygonToCells.c b/src/apps/testapps/testPolygonToCells.c
index 22a0cdecf..b35dedd58 100644
--- a/src/apps/testapps/testPolygonToCells.c
+++ b/src/apps/testapps/testPolygonToCells.c
@@ -461,19 +461,21 @@ SUITE(polygonToCells) {
 
     TEST(invalidFlags) {
         int64_t numHexagons;
-        for (uint32_t flags = 1; flags <= 32; flags++) {
+        for (uint32_t flags = 3; flags <= 32; flags++) {
             t_assert(
                 H3_EXPORT(maxPolygonToCellsSize)(
                     &sfGeoPolygon, 9, flags, &numHexagons) == E_OPTION_INVALID,
-                "Flags other than 0 are invalid for maxPolygonToCellsSize");
+                "Flags other than polyfill modes are invalid for "
+                "maxPolygonToCellsSize");
         }
         t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0,
                                                          &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
-        for (uint32_t flags = 1; flags <= 32; flags++) {
+        for (uint32_t flags = 3; flags <= 32; flags++) {
             t_assert(H3_EXPORT(polygonToCells)(&sfGeoPolygon, 9, flags,
                                                hexagons) == E_OPTION_INVALID,
-                     "Flags other than 0 are invalid for polygonToCells");
+                     "Flags other than polyfill modes are invalid for "
+                     "polygonToCells");
         }
         free(hexagons);
     }
diff --git a/src/apps/testapps/testPolygonToCellsExperimental.c b/src/apps/testapps/testPolygonToCellsExperimental.c
index f389456ab..24464597d 100644
--- a/src/apps/testapps/testPolygonToCellsExperimental.c
+++ b/src/apps/testapps/testPolygonToCellsExperimental.c
@@ -22,6 +22,7 @@
 #include "h3Index.h"
 #include "latLng.h"
 #include "polyfill.h"
+#include "polygon.h"
 #include "test.h"
 #include "utility.h"
 
@@ -99,7 +100,7 @@ static void fillIndex_assertions(H3Index h) {
         H3Index *polygonToCellsOut =
             calloc(polygonToCellsSize, sizeof(H3Index));
         t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
-            &polygon, nextRes, 0, polygonToCellsOut));
+            &polygon, nextRes, CONTAINMENT_CENTER, polygonToCellsOut));
 
         int64_t polygonToCellsCount =
             countNonNullIndexes(polygonToCellsOut, polygonToCellsSize);
@@ -160,26 +161,56 @@ SUITE(polygonToCells) {
 
     TEST(polygonToCells) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &sfGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&sfGeoPolygon, 9,
-                                                              0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &sfGeoPolygon, 9, CONTAINMENT_CENTER, hexagons));
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
         t_assert(actualNumIndexes == 1253, "got expected polygonToCells size");
         free(hexagons);
     }
 
+    TEST(polygonToCells_FullContainment) {
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &sfGeoPolygon, 9, CONTAINMENT_FULL, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &sfGeoPolygon, 9, CONTAINMENT_FULL, hexagons));
+        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+
+        t_assert(actualNumIndexes == 1175,
+                 "got expected polygonToCells size (full containment mode)");
+        free(hexagons);
+    }
+
+    TEST(polygonToCells_Overlapping) {
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &sfGeoPolygon, 9, CONTAINMENT_OVERLAPPING, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &sfGeoPolygon, 9, CONTAINMENT_OVERLAPPING, hexagons));
+        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+
+        t_assert(actualNumIndexes == 1334,
+                 "got expected polygonToCells size (overlapping mode)");
+        free(hexagons);
+    }
+
     TEST(polygonToCellsHole) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&holeGeoPolygon, 9, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &holeGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&holeGeoPolygon,
-                                                              9, 0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &holeGeoPolygon, 9, CONTAINMENT_CENTER, hexagons));
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
         t_assert(actualNumIndexes == 1214,
@@ -187,14 +218,45 @@ SUITE(polygonToCells) {
         free(hexagons);
     }
 
+    TEST(polygonToCellsHoleFullContainment) {
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &holeGeoPolygon, 9, CONTAINMENT_FULL, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &holeGeoPolygon, 9, CONTAINMENT_FULL, hexagons));
+        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+
+        t_assert(
+            actualNumIndexes == 1118,
+            "got expected polygonToCells size (hole, full containment mode)");
+        free(hexagons);
+    }
+
+    TEST(polygonToCellsHoleOverlapping) {
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &holeGeoPolygon, 9, CONTAINMENT_OVERLAPPING, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &holeGeoPolygon, 9, CONTAINMENT_OVERLAPPING, hexagons));
+        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+
+        t_assert(actualNumIndexes == 1311,
+                 "got expected polygonToCells size (hole, overlapping mode)");
+        free(hexagons);
+    }
+
     TEST(polygonToCellsEmpty) {
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&emptyGeoPolygon, 9, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &emptyGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&emptyGeoPolygon,
-                                                              9, 0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &emptyGeoPolygon, 9, CONTAINMENT_CENTER, hexagons));
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
         t_assert(actualNumIndexes == 0,
@@ -202,6 +264,77 @@ SUITE(polygonToCells) {
         free(hexagons);
     }
 
+    TEST(polygonToCellsContainsPolygon) {
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &sfGeoPolygon, 4, CONTAINMENT_CENTER, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &sfGeoPolygon, 4, CONTAINMENT_CENTER, hexagons));
+        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+
+        t_assert(actualNumIndexes == 0, "got expected polygonToCells size");
+        free(hexagons);
+    }
+
+    TEST(polygonToCellsContainsPolygon_CenterContained) {
+        // Contains the center point of a res 4 cell
+        LatLng centerVerts[] = {{0.6595645, -2.1353315},
+                                {0.6595645, -2.1353314},
+                                {0.6595644, -2.1353314},
+                                {0.6595644, -2.1353314265}};
+        GeoLoop centerGeoLoop = {.numVerts = 4, .verts = centerVerts};
+        GeoPolygon centerGeoPolygon;
+        centerGeoPolygon.geoloop = centerGeoLoop;
+        centerGeoPolygon.numHoles = 0;
+
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &centerGeoPolygon, 4, CONTAINMENT_CENTER, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &centerGeoPolygon, 4, CONTAINMENT_CENTER, hexagons));
+        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+
+        t_assert(actualNumIndexes == 1, "got expected polygonToCells size");
+        t_assert(hexagons[0] == 0x8428309ffffffff, "got expected hexagon");
+
+        free(hexagons);
+    }
+
+    TEST(polygonToCellsContainsPolygon_FullContainment) {
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &sfGeoPolygon, 4, CONTAINMENT_FULL, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &sfGeoPolygon, 4, CONTAINMENT_FULL, hexagons));
+        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+
+        t_assert(actualNumIndexes == 0,
+                 "got expected polygonToCells size (full containment mode)");
+        free(hexagons);
+    }
+
+    TEST(polygonToCellsContainsPolygon_Overlapping) {
+        int64_t numHexagons;
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &sfGeoPolygon, 4, CONTAINMENT_OVERLAPPING, &numHexagons));
+        H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &sfGeoPolygon, 4, CONTAINMENT_OVERLAPPING, hexagons));
+        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+
+        t_assert(actualNumIndexes == 1,
+                 "got expected polygonToCells size (overlapping mode)");
+        t_assert(hexagons[0] == 0x8428309ffffffff, "got expected hexagon");
+        free(hexagons);
+    }
+
     TEST(polygonToCellsExact) {
         LatLng somewhere = {1, 2};
         H3Index origin;
@@ -223,15 +356,27 @@ SUITE(polygonToCells) {
         someHexagon.numHoles = 0;
 
         int64_t numHexagons;
-        t_assertSuccess(
-            H3_EXPORT(maxPolygonToCellsSize)(&someHexagon, 9, 0, &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &someHexagon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&someHexagon, 9,
-                                                              0, hexagons));
-        int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+        int64_t actualNumIndexes;
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &someHexagon, 9, CONTAINMENT_CENTER, hexagons));
+        actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+        t_assert(actualNumIndexes == 1,
+                 "got expected polygonToCells size for center containment (1)");
+
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &someHexagon, 9, CONTAINMENT_FULL, hexagons));
+        actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
+        t_assert(actualNumIndexes == 1,
+                 "got expected polygonToCells size for full containment (1)");
+
+        // TODO: CONTAINMENT_OVERLAPPING yields 7 cells, presumably due to FPE
+        // in the various cell boundaries
 
-        t_assert(actualNumIndexes == 1, "got expected polygonToCells size (1)");
         free(hexagons);
         free(verts);
     }
@@ -272,11 +417,11 @@ SUITE(polygonToCells) {
         expectedSize = 4228;
         int64_t numHexagons;
         t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
-            &primeMeridianGeoPolygon, 7, 0, &numHexagons));
+            &primeMeridianGeoPolygon, 7, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
         t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
-            &primeMeridianGeoPolygon, 7, 0, hexagons));
+            &primeMeridianGeoPolygon, 7, CONTAINMENT_CENTER, hexagons));
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
         t_assert(actualNumIndexes == expectedSize,
@@ -287,11 +432,11 @@ SUITE(polygonToCells) {
         // differences in hex size and grid offset between the two cases
         expectedSize = 4238;
         t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
-            &transMeridianGeoPolygon, 7, 0, &numHexagons));
+            &transMeridianGeoPolygon, 7, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagonsTM = calloc(numHexagons, sizeof(H3Index));
 
         t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
-            &transMeridianGeoPolygon, 7, 0, hexagonsTM));
+            &transMeridianGeoPolygon, 7, CONTAINMENT_CENTER, hexagonsTM));
         actualNumIndexes = countNonNullIndexes(hexagonsTM, numHexagons);
 
         t_assert(actualNumIndexes == expectedSize,
@@ -300,21 +445,23 @@ SUITE(polygonToCells) {
         // Transmeridian filled hole case -- only needed for calculating hole
         // size
         t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
-            &transMeridianFilledHoleGeoPolygon, 7, 0, &numHexagons));
+            &transMeridianFilledHoleGeoPolygon, 7, CONTAINMENT_CENTER,
+            &numHexagons));
         H3Index *hexagonsTMFH = calloc(numHexagons, sizeof(H3Index));
 
         t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
-            &transMeridianFilledHoleGeoPolygon, 7, 0, hexagonsTMFH));
+            &transMeridianFilledHoleGeoPolygon, 7, CONTAINMENT_CENTER,
+            hexagonsTMFH));
         int64_t actualNumHoleIndexes =
             countNonNullIndexes(hexagonsTMFH, numHexagons);
 
         // Transmeridian hole case
         t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
-            &transMeridianHoleGeoPolygon, 7, 0, &numHexagons));
+            &transMeridianHoleGeoPolygon, 7, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagonsTMH = calloc(numHexagons, sizeof(H3Index));
 
         t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
-            &transMeridianHoleGeoPolygon, 7, 0, hexagonsTMH));
+            &transMeridianHoleGeoPolygon, 7, CONTAINMENT_CENTER, hexagonsTMH));
         actualNumIndexes = countNonNullIndexes(hexagonsTMH, numHexagons);
 
         t_assert(actualNumIndexes == expectedSize - actualNumHoleIndexes,
@@ -337,12 +484,12 @@ SUITE(polygonToCells) {
         GeoPolygon polygon = {.geoloop = geoloop, .numHoles = 0};
 
         int64_t numHexagons;
-        t_assertSuccess(
-            H3_EXPORT(maxPolygonToCellsSize)(&polygon, 4, 0, &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &polygon, 4, CONTAINMENT_CENTER, &numHexagons));
 
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
-        t_assertSuccess(
-            H3_EXPORT(polygonToCellsExperimental)(&polygon, 4, 0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &polygon, 4, CONTAINMENT_CENTER, hexagons));
 
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
@@ -389,12 +536,12 @@ SUITE(polygonToCells) {
         polygon.numHoles = 0;
 
         int64_t numHexagons;
-        t_assertSuccess(
-            H3_EXPORT(maxPolygonToCellsSize)(&polygon, 9, 0, &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &polygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(
-            H3_EXPORT(polygonToCellsExperimental)(&polygon, 9, 0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &polygon, 9, CONTAINMENT_CENTER, hexagons));
 
         int found = 0;
         int numPentagons = 0;
@@ -413,19 +560,21 @@ SUITE(polygonToCells) {
 
     TEST(invalidFlags) {
         int64_t numHexagons;
-        for (uint32_t flags = 1; flags <= 32; flags++) {
+        for (uint32_t flags = 3; flags <= 32; flags++) {
             t_assert(
                 H3_EXPORT(maxPolygonToCellsSize)(
                     &sfGeoPolygon, 9, flags, &numHexagons) == E_OPTION_INVALID,
-                "Flags other than 0 are invalid for maxPolygonToCellsSize");
+                "Flags other than polyfill modes are invalid for "
+                "maxPolygonToCellsSize");
         }
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&sfGeoPolygon, 9, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &sfGeoPolygon, 9, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
-        for (uint32_t flags = 1; flags <= 32; flags++) {
+        for (uint32_t flags = 3; flags <= 32; flags++) {
             t_assert(H3_EXPORT(polygonToCellsExperimental)(
                          &sfGeoPolygon, 9, flags, hexagons) == E_OPTION_INVALID,
-                     "Flags other than 0 are invalid for polygonToCells");
+                     "Flags other than polyfill modes are invalid for "
+                     "polygonToCells");
         }
         free(hexagons);
     }
diff --git a/src/apps/testapps/testPolygonToCellsReportedExperimental.c b/src/apps/testapps/testPolygonToCellsReportedExperimental.c
index 58e688504..b3fe8e5b8 100644
--- a/src/apps/testapps/testPolygonToCellsReportedExperimental.c
+++ b/src/apps/testapps/testPolygonToCellsReportedExperimental.c
@@ -21,6 +21,7 @@
 #include "h3Index.h"
 #include "latLng.h"
 #include "polyfill.h"
+#include "polygon.h"
 #include "test.h"
 #include "utility.h"
 
@@ -42,23 +43,26 @@ SUITE(polygonToCells_reported) {
         for (int res = 0; res < 3; res++) {
             int64_t polygonToCellsSize;
             t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
-                &worldGeoPolygon, res, 0, &polygonToCellsSize));
+                &worldGeoPolygon, res, CONTAINMENT_CENTER,
+                &polygonToCellsSize));
             H3Index *polygonToCellsOut =
                 calloc(polygonToCellsSize, sizeof(H3Index));
 
             t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
-                &worldGeoPolygon, res, 0, polygonToCellsOut));
+                &worldGeoPolygon, res, CONTAINMENT_CENTER, polygonToCellsOut));
             int64_t actualNumIndexes =
                 countNonNullIndexes(polygonToCellsOut, polygonToCellsSize);
 
             int64_t polygonToCellsSize2;
             t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
-                &worldGeoPolygon2, res, 0, &polygonToCellsSize2));
+                &worldGeoPolygon2, res, CONTAINMENT_CENTER,
+                &polygonToCellsSize2));
             H3Index *polygonToCellsOut2 =
                 calloc(polygonToCellsSize2, sizeof(H3Index));
 
             t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
-                &worldGeoPolygon2, res, 0, polygonToCellsOut2));
+                &worldGeoPolygon2, res, CONTAINMENT_CENTER,
+                polygonToCellsOut2));
             int64_t actualNumIndexes2 =
                 countNonNullIndexes(polygonToCellsOut2, polygonToCellsSize2);
 
@@ -105,12 +109,12 @@ SUITE(polygonToCells_reported) {
 
         int res = 7;
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&testPolygon, res, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &testPolygon, res, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&testPolygon, res,
-                                                              0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &testPolygon, res, CONTAINMENT_CENTER, hexagons));
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
         t_assert(actualNumIndexes == 4499,
@@ -134,12 +138,12 @@ SUITE(polygonToCells_reported) {
 
         int res = 7;
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&testPolygon, res, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &testPolygon, res, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&testPolygon, res,
-                                                              0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &testPolygon, res, CONTAINMENT_CENTER, hexagons));
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
         t_assert(actualNumIndexes == 4609,
@@ -160,12 +164,12 @@ SUITE(polygonToCells_reported) {
 
         int res = 13;
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&testPolygon, res, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &testPolygon, res, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&testPolygon, res,
-                                                              0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &testPolygon, res, CONTAINMENT_CENTER, hexagons));
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
         t_assert(actualNumIndexes == 4353, "got expected polygonToCells size");
@@ -193,12 +197,12 @@ SUITE(polygonToCells_reported) {
 
         int res = 5;
         int64_t numHexagons;
-        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(&testPolygon, res, 0,
-                                                         &numHexagons));
+        t_assertSuccess(H3_EXPORT(maxPolygonToCellsSize)(
+            &testPolygon, res, CONTAINMENT_CENTER, &numHexagons));
         H3Index *hexagons = calloc(numHexagons, sizeof(H3Index));
 
-        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(&testPolygon, res,
-                                                              0, hexagons));
+        t_assertSuccess(H3_EXPORT(polygonToCellsExperimental)(
+            &testPolygon, res, CONTAINMENT_CENTER, hexagons));
         int64_t actualNumIndexes = countNonNullIndexes(hexagons, numHexagons);
 
         t_assert(actualNumIndexes == 8, "got expected polygonToCells size");
diff --git a/src/h3lib/include/polygon.h b/src/h3lib/include/polygon.h
index 387ed034a..d7d24ffbc 100644
--- a/src/h3lib/include/polygon.h
+++ b/src/h3lib/include/polygon.h
@@ -40,18 +40,38 @@
 /** Macro: Whether a GeoLoop is empty */
 #define IS_EMPTY_GEOFENCE(geoloop) geoloop->numVerts == 0
 
+/**
+ * Values representing polyfill containment modes, to be used in
+ * the `flags` bit field.
+ */
+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
+} ContainmentMode;
+
+// 1s in the 4 bits defining the polyfill containment mode, 0s elsewhere
+#define FLAG_CONTAINMENT_MODE_MASK ((uint32_t)(15))
+#define FLAG_GET_CONTAINMENT_MODE(flags) (flags & FLAG_CONTAINMENT_MODE_MASK)
+
 // Defined directly in polygon.c:
+H3Error validatePolygonFlags(uint32_t flags);
 void bboxesFromGeoPolygon(const GeoPolygon *polygon, BBox *bboxes);
 bool pointInsidePolygon(const GeoPolygon *geoPolygon, const BBox *bboxes,
                         const LatLng *coord);
 bool cellBoundaryInsidePolygon(const GeoPolygon *geoPolygon, const BBox *bboxes,
                                const CellBoundary *boundary,
                                const BBox *boundaryBBox);
+bool cellBoundaryCrossesPolygon(const GeoPolygon *geoPolygon,
+                                const BBox *bboxes,
+                                const CellBoundary *boundary,
+                                const BBox *boundaryBBox);
 bool cellBoundaryCrossesGeoLoop(const GeoLoop *geoloop, const BBox *loopBBox,
                                 const CellBoundary *boundary,
                                 const BBox *boundaryBBox);
-bool lineIntersectsLine(const LatLng *a1, const LatLng *a2, const LatLng *b1,
-                        const LatLng *b2);
+bool lineCrossesLine(const LatLng *a1, const LatLng *a2, const LatLng *b1,
+                     const LatLng *b2);
 
 // The following functions are created via macro in polygonAlgos.h,
 // so their signatures are documented here:
diff --git a/src/h3lib/lib/algos.c b/src/h3lib/lib/algos.c
index 79666d843..949988a08 100644
--- a/src/h3lib/lib/algos.c
+++ b/src/h3lib/lib/algos.c
@@ -774,8 +774,9 @@ H3Error H3_EXPORT(gridRingUnsafe)(H3Index origin, int k, H3Index *out) {
  */
 H3Error H3_EXPORT(maxPolygonToCellsSize)(const GeoPolygon *geoPolygon, int res,
                                          uint32_t flags, int64_t *out) {
-    if (flags != 0) {
-        return E_OPTION_INVALID;
+    H3Error flagErr = validatePolygonFlags(flags);
+    if (flagErr) {
+        return flagErr;
     }
     // Get the bounding box for the GeoJSON-like struct
     BBox bbox;
@@ -890,8 +891,9 @@ H3Error _getEdgeHexagons(const GeoLoop *geoloop, int64_t numHexagons, int res,
  */
 H3Error H3_EXPORT(polygonToCells)(const GeoPolygon *geoPolygon, int res,
                                   uint32_t flags, H3Index *out) {
-    if (flags != 0) {
-        return E_OPTION_INVALID;
+    H3Error flagErr = validatePolygonFlags(flags);
+    if (flagErr) {
+        return flagErr;
     }
     // One of the goals of the polygonToCells algorithm is that two adjacent
     // polygons with zero overlap have zero overlapping hexagons. That the
diff --git a/src/h3lib/lib/polyfill.c b/src/h3lib/lib/polyfill.c
index 7ba5a28b5..bde5eedbc 100644
--- a/src/h3lib/lib/polyfill.c
+++ b/src/h3lib/lib/polyfill.c
@@ -354,8 +354,9 @@ IterCellsPolygonCompact iterInitPolygonCompact(const GeoPolygon *polygon,
         return iter;
     }
 
-    if (flags != 0) {
-        iterErrorPolygonCompact(&iter, E_OPTION_INVALID);
+    H3Error flagErr = validatePolygonFlags(flags);
+    if (flagErr) {
+        iterErrorPolygonCompact(&iter, flagErr);
         return iter;
     }
 
@@ -406,23 +407,80 @@ void iterStepPolygonCompact(IterCellsPolygonCompact *iter) {
         iter->_started = true;
     }
 
+    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) {
-            // Check if the cell is in the polygon
-            // TODO: Handle other polyfill modes here
-            LatLng center;
-            H3Error centerErr = H3_EXPORT(cellToLatLng)(cell, &center);
-            if (NEVER(centerErr != E_SUCCESS)) {
-                iterErrorPolygonCompact(iter, centerErr);
-                return;
+            if (mode == CONTAINMENT_CENTER || mode == CONTAINMENT_OVERLAPPING) {
+                // 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 (pointInsidePolygon(iter->_polygon, iter->_bboxes, &center)) {
-                // Set to next output
-                iter->cell = cell;
-                return;
+            if (mode == CONTAINMENT_OVERLAPPING) {
+                // 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.
+                H3Index polygonCell;
+                H3Error polygonCellErr = H3_EXPORT(latLngToCell)(
+                    &(iter->_polygon->geoloop.verts[0]), cellRes, &polygonCell);
+                if (polygonCellErr != E_SUCCESS) {
+                    iterErrorPolygonCompact(iter, polygonCellErr);
+                    return;
+                }
+                if (polygonCell == cell) {
+                    // Set to next output
+                    iter->cell = cell;
+                    return;
+                }
+            }
+            if (mode == CONTAINMENT_FULL || mode == CONTAINMENT_OVERLAPPING) {
+                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 &&
+                    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 &&
+                         cellBoundaryCrossesPolygon(
+                             iter->_polygon, iter->_bboxes, &boundary, &bbox)) {
+                    // Set to next output
+                    iter->cell = cell;
+                    return;
+                }
             }
         }
 
diff --git a/src/h3lib/lib/polygon.c b/src/h3lib/lib/polygon.c
index a92fecfba..1944e5bda 100644
--- a/src/h3lib/lib/polygon.c
+++ b/src/h3lib/lib/polygon.c
@@ -42,6 +42,20 @@
 #undef ITERATE
 #undef IS_EMPTY
 
+/**
+ * Whether the flags for the polyfill operation are valid
+ * TODO: Move to polyfill.c when the old algo is removed
+ * @param  flags Flags to validate
+ * @return       Whether the flags are valid
+ */
+H3Error validatePolygonFlags(uint32_t flags) {
+    if (flags & (~FLAG_CONTAINMENT_MODE_MASK) ||
+        FLAG_GET_CONTAINMENT_MODE(flags) >= CONTAINMENT_INVALID) {
+        return E_OPTION_INVALID;
+    }
+    return E_SUCCESS;
+}
+
 /**
  * Create a bounding box from a GeoPolygon
  * @param polygon Input GeoPolygon
@@ -115,6 +129,33 @@ bool cellBoundaryInsidePolygon(const GeoPolygon *geoPolygon, const BBox *bboxes,
     return true;
 }
 
+/**
+ * Whether any part of a cell boundary crosses a polygon. Crossing in this case
+ * means whether any line segments intersect; it does not include containment.
+ * @param  geoPolygon The polygon to test
+ * @param  bboxes     The bboxes for the main geoloop and each of its holes
+ * @param  boundary   The cell boundary to test
+ * @return            Whether the cell boundary is contained
+ */
+bool cellBoundaryCrossesPolygon(const GeoPolygon *geoPolygon,
+                                const BBox *bboxes,
+                                const CellBoundary *boundary,
+                                const BBox *boundaryBBox) {
+    // Check for line intersections with outer loop
+    if (cellBoundaryCrossesGeoLoop(&(geoPolygon->geoloop), &bboxes[0], boundary,
+                                   boundaryBBox)) {
+        return true;
+    }
+    // Check for line intersections with any hole
+    for (int i = 0; i < geoPolygon->numHoles; i++) {
+        if (cellBoundaryCrossesGeoLoop(&(geoPolygon->holes[i]), &bboxes[i + 1],
+                                       boundary, boundaryBBox)) {
+            return true;
+        }
+    }
+    return false;
+}
+
 /**
  * Whether a cell boundary crosses a geo loop. Crossing in this case means
  * whether any line segments intersect; it does not include containment.
@@ -167,7 +208,7 @@ bool cellBoundaryCrossesGeoLoop(const GeoLoop *geoloop, const BBox *loopBBox,
         }
 
         for (int j = 0; j < normalBoundary.numVerts; j++) {
-            if (lineIntersectsLine(
+            if (lineCrossesLine(
                     &loop1, &loop2, &normalBoundary.verts[j],
                     &normalBoundary.verts[(j + 1) % normalBoundary.numVerts])) {
                 return true;
@@ -187,8 +228,8 @@ bool cellBoundaryCrossesGeoLoop(const GeoLoop *geoloop, const BBox *loopBBox,
  * @param  b2 End of line B
  * @return    Whether the lines intersect
  */
-bool lineIntersectsLine(const LatLng *a1, const LatLng *a2, const LatLng *b1,
-                        const LatLng *b2) {
+bool lineCrossesLine(const LatLng *a1, const LatLng *a2, const LatLng *b1,
+                     const LatLng *b2) {
     double denom = ((b2->lng - b1->lng) * (a2->lat - a1->lat) -
                     (b2->lat - b1->lat) * (a2->lng - a1->lng));
     if (!denom) return false;