diff --git a/h3.go b/h3.go
index b04e4cd..e684585 100644
--- a/h3.go
+++ b/h3.go
@@ -25,6 +25,8 @@ package h3
 #include <stdlib.h>
 #include <h3_h3api.h>
 #include <h3_h3Index.h>
+#include <h3_polygon.h>
+#include <h3_polyfill.h>
 */
 import "C"
 
@@ -66,6 +68,13 @@ const (
 
 	DegsToRads = math.Pi / 180.0
 	RadsToDegs = 180.0 / math.Pi
+
+	// PolygonToCells containment modes
+	ContainmentCenter          ContainmentMode = C.CONTAINMENT_CENTER           // Cell center is contained in the shape
+	ContainmentFull            ContainmentMode = C.CONTAINMENT_FULL             // Cell is fully contained in the shape
+	ContainmentOverlapping     ContainmentMode = C.CONTAINMENT_OVERLAPPING      // Cell overlaps the shape at any point
+	ContainmentOverlappingBbox ContainmentMode = C.CONTAINMENT_OVERLAPPING_BBOX // Cell bounding box overlaps shape
+	ContainmentInvalid         ContainmentMode = C.CONTAINMENT_INVALID          // This mode is invalid and should not be used
 )
 
 // Error codes.
@@ -137,6 +146,9 @@ type (
 		GeoLoop GeoLoop
 		Holes   []GeoLoop
 	}
+
+	// ContainmentMode is an int for specifying PolygonToCell containment behavior.
+	ContainmentMode C.uint32_t
 )
 
 func NewLatLng(lat, lng float64) LatLng {
@@ -222,6 +234,7 @@ func GridDiskDistances(origin Cell, k int) ([][]Cell, error) {
 	rsz := maxGridDiskSize(k)
 	outHexes := make([]C.H3Index, rsz)
 	outDists := make([]C.int, rsz)
+
 	if err := toErr(C.gridDiskDistances(C.H3Index(origin), C.int(k), &outHexes[0], &outDists[0])); err != nil {
 		return nil, err
 	}
@@ -276,6 +289,36 @@ func PolygonToCells(polygon GeoPolygon, resolution int) ([]Cell, error) {
 	return cellsFromC(out, true, false), toErr(errC)
 }
 
+// PolygonToCells takes a given GeoJSON-like data structure fills it with the
+// hexagon cells that are contained by the GeoJSON-like data structure.
+//
+// This implementation traces the GeoJSON geoloop(s) in cartesian space with
+// hexagons, tests them and their neighbors to be contained by the geoloop(s),
+// and then any newly found hexagons are used to test again until no new
+// hexagons are found.
+func PolygonToCellsExperimental(polygon GeoPolygon, resolution int, mode ContainmentMode, maxNumCellsReturn ...int64) ([]Cell, error) {
+	var maxNumCells int64 = math.MaxInt64
+	if len(maxNumCellsReturn) > 0 {
+		maxNumCells = maxNumCellsReturn[0]
+	}
+	if len(polygon.GeoLoop) == 0 {
+		return nil, nil
+	}
+	cpoly := allocCGeoPolygon(polygon)
+
+	defer freeCGeoPolygon(&cpoly)
+
+	maxLen := new(C.int64_t)
+	if err := toErr(C.maxPolygonToCellsSizeExperimental(&cpoly, C.int(resolution), C.uint32_t(mode), maxLen)); err != nil {
+		return nil, err
+	}
+
+	out := make([]C.H3Index, *maxLen)
+	errC := C.polygonToCellsExperimental(&cpoly, C.int(resolution), C.uint32_t(mode), C.int64_t(maxNumCells), &out[0])
+
+	return cellsFromC(out, true, false), toErr(errC)
+}
+
 // PolygonToCells takes a given GeoJSON-like data structure fills it with the
 // hexagon cells that are contained by the GeoJSON-like data structure.
 //
diff --git a/h3_test.go b/h3_test.go
index f269ba4..39d31ef 100644
--- a/h3_test.go
+++ b/h3_test.go
@@ -637,6 +637,90 @@ func TestCellsToMultiPolygon(t *testing.T) {
 	assertNil(t, res)
 }
 
+func TestPolygonToCellsExperimental(t *testing.T) {
+	t.Parallel()
+
+	t.Run("empty", func(t *testing.T) {
+		t.Parallel()
+
+		for _, flag := range []ContainmentMode{
+			ContainmentCenter,
+			ContainmentFull,
+			ContainmentOverlapping,
+			ContainmentOverlappingBbox,
+		} {
+			cells, err := PolygonToCellsExperimental(GeoPolygon{}, 6, flag)
+			if err != nil {
+				t.Error(t)
+			}
+
+			assertEqual(t, 0, len(cells))
+		}
+	})
+
+	t.Run("without holes", func(t *testing.T) {
+		t.Parallel()
+
+		for _, flag := range []ContainmentMode{
+			ContainmentCenter,
+			ContainmentFull,
+			ContainmentOverlapping,
+			ContainmentOverlappingBbox,
+		} {
+			cells, err := PolygonToCellsExperimental(validGeoPolygonNoHoles, 6, flag)
+			if err != nil {
+				t.Error(t)
+			}
+			expectedCellCounts := map[ContainmentMode]int{
+				ContainmentCenter:          7,
+				ContainmentFull:            1,
+				ContainmentOverlapping:     14,
+				ContainmentOverlappingBbox: 21,
+			}
+			assertEqual(t, expectedCellCounts[flag], len(cells))
+		}
+	})
+
+	t.Run("with holes", func(t *testing.T) {
+		t.Parallel()
+
+		for _, flag := range []ContainmentMode{
+			ContainmentCenter,
+			ContainmentFull,
+			ContainmentOverlapping,
+			ContainmentOverlappingBbox,
+		} {
+			cells, err := PolygonToCellsExperimental(validGeoPolygonHoles, 6, flag)
+			if err != nil {
+				t.Error(t)
+			}
+			expectedCellCounts := map[ContainmentMode]int{
+				ContainmentCenter:          6,
+				ContainmentFull:            0,
+				ContainmentOverlapping:     14,
+				ContainmentOverlappingBbox: 21,
+			}
+
+			assertEqual(t, expectedCellCounts[flag], len(cells))
+		}
+	})
+
+	t.Run("busting memory", func(t *testing.T) {
+		t.Parallel()
+
+		for _, flag := range []ContainmentMode{
+			ContainmentCenter,
+			ContainmentOverlapping,
+			ContainmentOverlappingBbox,
+		} {
+			_, err := PolygonToCellsExperimental(validGeoPolygonHoles, 6, flag, 3)
+			if err != ErrMemoryBounds {
+				t.Error(t)
+			}
+		}
+	})
+}
+
 func TestGridPath(t *testing.T) {
 	t.Parallel()
 	path, err := lineStartCell.GridPath(lineEndCell)