diff --git a/s2/polygon.go b/s2/polygon.go
index c14ccdf..7a1735f 100644
--- a/s2/polygon.go
+++ b/s2/polygon.go
@@ -184,6 +184,175 @@ func PolygonFromOrientedLoops(loops []*Loop) *Polygon {
 	return p
 }
 
+// PolygonOfCellUnionBorder creates a Polygon representing the border of the
+// given CellUnion. The CellUnion must be normalized, and should represent a
+// single connected region without discontiguous cells.
+//
+// This function extracts the boundary edges of the CellUnion (edges not shared
+// between adjacent cells) and stitches them into closed loops to form the
+// polygon boundary.
+//
+// Note: The C++ library implements InitToCellUnionBorder using S2Builder,
+// which handles edge snapping for numerical precision. This Go implementation
+// uses edge-stitching since S2Builder is not yet ported.
+func PolygonOfCellUnionBorder(cu CellUnion) (*Polygon, error) {
+	if len(cu) == 0 {
+		return PolygonFromOrientedLoops(nil), nil
+	}
+
+	level := 0
+	for _, cellID := range cu {
+		if l := cellID.Level(); l > level {
+			level = l
+		}
+	}
+
+	// Build a map of directed edges for border detection. Each cell contributes
+	// its 4 edges in counterclockwise order. Internal edges appear as both (A,B)
+	// and (B,A) from adjacent cells and cancel out; border edges appear only
+	// once. For mixed-level cells, edges are expanded for matching granularity.
+	directedEdges := make(map[Edge][]Edge)
+	for _, cellID := range cu {
+		for edgeIdx := range 4 {
+			cells := []CellID{cellID}
+			edges := []int{edgeIdx}
+			for i := 0; i < len(cells); i++ {
+				cell := CellFromCellID(cells[i])
+				vertexA := cell.Vertex(edges[i])
+				vertexB := cell.Vertex((edges[i] + 1) % 4)
+
+				if cells[i].Level() == level {
+					edge := Edge{V0: vertexA, V1: vertexB}
+					key := edge
+					if edge.V0.Vector.Cmp(edge.V1.Vector) > 0 {
+						key = edge.Reversed()
+					}
+					directedEdges[key] = append(directedEdges[key], edge)
+					continue
+				}
+
+				for _, child := range cells[i].Children() {
+					childCell := CellFromCellID(child)
+
+					// Find the child edge on the parent edge. Each child has at most
+					// one edge on any given parent edge.
+					for j := range 4 {
+						a := childCell.Vertex(j)
+						b := childCell.Vertex((j + 1) % 4)
+						if Project(a, vertexA, vertexB).ApproxEqual(a) &&
+							Project(b, vertexA, vertexB).ApproxEqual(b) {
+							cells = append(cells, child)
+							edges = append(edges, j)
+							break
+						}
+					}
+				}
+			}
+		}
+	}
+
+	// Edges appearing exactly once are border edges.
+	borderEdges := make(map[Point][]Point)
+	totalBorderEdges := 0
+	for _, directed := range directedEdges {
+		if len(directed) == 1 {
+			borderEdges[directed[0].V0] = append(borderEdges[directed[0].V0], directed[0].V1)
+			totalBorderEdges++
+		}
+	}
+
+	if len(borderEdges) == 0 {
+		return FullPolygon(), nil
+	}
+
+	var loops []*Loop
+	used := make(map[Edge]struct{})
+	for origin := range borderEdges {
+		for _, destination := range borderEdges[origin] {
+			edgeFirst := Edge{V0: origin, V1: destination}
+			if _, ok := used[edgeFirst]; ok {
+				continue
+			}
+			used[edgeFirst] = struct{}{}
+			vertex := destination
+			vertices := []Point{origin}
+
+			// A loop cannot have more edges than total border edges.
+			for range totalBorderEdges {
+				if vertex == origin {
+					break
+				}
+
+				vertices = append(vertices, vertex)
+
+				var vertexNext Point
+				found := false
+				for _, candidate := range borderEdges[vertex] {
+					edge := Edge{V0: vertex, V1: candidate}
+					if _, ok := used[edge]; !ok {
+						found = true
+						used[edge] = struct{}{}
+						vertexNext = candidate
+						break
+					}
+				}
+
+				if !found {
+					return nil, fmt.Errorf(
+						"broken edge chain at vertex %v: no unused outgoing edge",
+						vertex,
+					)
+				}
+
+				vertex = vertexNext
+			}
+
+			if vertex != origin {
+				return nil, fmt.Errorf(
+					"loop did not close: started at %v, ended at %v",
+					origin,
+					vertex,
+				)
+			}
+
+			if len(vertices) < 3 {
+				return nil, fmt.Errorf(
+					"degenerate loop with %d vertices at %v",
+					len(vertices),
+					origin,
+				)
+			}
+			loops = append(loops, LoopFromPoints(vertices))
+		}
+	}
+
+	if len(loops) == 0 {
+		return &Polygon{}, nil
+	}
+
+	for _, loop := range loops {
+		if loop.TurningAngle() < 0 {
+			loop.Invert()
+		}
+	}
+
+	contained := false
+	sampleCenter := CellFromCellID(cu[0]).Center()
+	for _, loop := range loops {
+		if loop.ContainsPoint(sampleCenter) {
+			contained = true
+			break
+		}
+	}
+	if !contained {
+		loops = append([]*Loop{FullLoop()}, loops...)
+	}
+
+	p := &Polygon{loops: loops}
+	p.initNested()
+	return p, nil
+}
+
 // Invert inverts the polygon (replaces it by its complement).
 func (p *Polygon) Invert() {
 	// Inverting any one loop will invert the polygon.  The best loop to invert
@@ -1217,7 +1386,6 @@ func (p *Polygon) decodeCompressed(d *decoder) {
 // ApproxSubtractFromPolyline
 // DestructiveUnion
 // DestructiveApproxUnion
-// InitToCellUnionBorder
 // IsNormalized
 // Equal/BoundaryEqual/BoundaryApproxEqual/BoundaryNear Polygons
 // BreakEdgesAndAddToBuilder
diff --git a/s2/polygon_test.go b/s2/polygon_test.go
index 52511a8..2cbdb86 100644
--- a/s2/polygon_test.go
+++ b/s2/polygon_test.go
@@ -1157,6 +1157,324 @@ func TestPolygonInvert(t *testing.T) {
 	}
 }
 
+func TestPolygonOfCellUnionBorder(t *testing.T) {
+	t.Run("empty cell union", func(t *testing.T) {
+		p, err := PolygonOfCellUnionBorder(CellUnion{})
+		if err != nil {
+			t.Fatalf("unexpected error: %v", err)
+		}
+		if !p.IsEmpty() {
+			t.Errorf("expected empty polygon, got %d loops", p.NumLoops())
+		}
+	})
+
+	t.Run("single cell", func(t *testing.T) {
+		cellID := CellIDFromFace(0).ChildBeginAtLevel(5)
+		cells := CellUnion{cellID}
+
+		p, err := PolygonOfCellUnionBorder(cells)
+		if err != nil {
+			t.Fatalf("unexpected error: %v", err)
+		}
+
+		if p.NumLoops() != 1 {
+			t.Errorf("expected 1 loop, got %d", p.NumLoops())
+		}
+		if p.Loop(0).NumVertices() != 4 {
+			t.Errorf("expected 4 vertices, got %d", p.Loop(0).NumVertices())
+		}
+		if !p.ContainsPoint(CellFromCellID(cellID).Center()) {
+			t.Error("polygon should contain cell center")
+		}
+	})
+
+	t.Run("two adjacent cells", func(t *testing.T) {
+		cell := CellIDFromFace(0).ChildBeginAtLevel(5)
+		neighbor := cell.EdgeNeighbors()[0]
+		cells := CellUnion{cell, neighbor}
+		cells.Normalize()
+
+		p, err := PolygonOfCellUnionBorder(cells)
+		if err != nil {
+			t.Fatalf("unexpected error: %v", err)
+		}
+
+		if p.NumLoops() != 1 {
+			t.Errorf("expected 1 loop, got %d", p.NumLoops())
+		}
+		// Two adjacent cells share one edge, so 4+4-2 = 6 vertices
+		if p.Loop(0).NumVertices() != 6 {
+			t.Errorf("expected 6 vertices, got %d", p.Loop(0).NumVertices())
+		}
+	})
+
+	t.Run("merged siblings (4 cells = square)", func(t *testing.T) {
+		parent := CellIDFromFace(1).ChildBeginAtLevel(8)
+		children := parent.Children()
+		cells := CellUnion{children[0], children[1], children[2], children[3]}
+		cells.Normalize()
+
+		p, err := PolygonOfCellUnionBorder(cells)
+		if err != nil {
+			t.Fatalf("unexpected error: %v", err)
+		}
+
+		if p.NumLoops() != 1 {
+			t.Fatalf("NumLoops() = %d, want 1", p.NumLoops())
+		}
+		// 4 sibling cells form a square with 4 vertices (internal edges cancel)
+		if p.Loop(0).NumVertices() != 4 {
+			t.Errorf("expected 4 vertices for 2x2 cell block, got %d", p.Loop(0).NumVertices())
+		}
+	})
+
+	t.Run("full sphere", func(t *testing.T) {
+		cells := make(CellUnion, NumFaces)
+		for face := 0; face < NumFaces; face++ {
+			cells[face] = CellIDFromFace(face)
+		}
+		cells.Normalize()
+
+		p, err := PolygonOfCellUnionBorder(cells)
+		if err != nil {
+			t.Fatalf("unexpected error: %v", err)
+		}
+
+		if !p.IsFull() {
+			t.Error("expected full polygon for complete sphere coverage")
+		}
+	})
+}
+
+func TestPolygonOfCellUnionBorder_Hole(t *testing.T) {
+	// Create a shape with a hole: a ring of cells surrounding an empty center.
+	center := CellIDFromFace(0).ChildBeginAtLevel(5)
+	neighbors := center.AllNeighbors(5)
+
+	// Use only the neighbors, not the center - creates a ring with a hole.
+	cells := CellUnion(neighbors)
+	cells.Normalize()
+
+	p, err := PolygonOfCellUnionBorder(cells)
+	if err != nil {
+		t.Fatalf("unexpected error: %v", err)
+	}
+
+	if p.NumLoops() != 2 {
+		t.Fatalf("NumLoops() = %d, want 2 (outer + hole)", p.NumLoops())
+	}
+
+	// The hole center should NOT be contained by the polygon.
+	centerPoint := CellFromCellID(center).Center()
+	if p.ContainsPoint(centerPoint) {
+		t.Error("polygon contains the hole center, should not")
+	}
+
+	// Each neighbor cell center should be contained by the polygon.
+	for _, neighborID := range neighbors {
+		neighborCenter := CellFromCellID(neighborID).Center()
+		if !p.ContainsPoint(neighborCenter) {
+			t.Errorf("polygon does not contain neighbor %v", neighborID)
+		}
+	}
+}
+
+func TestPolygonOfCellUnionBorder_FullSphereWithHole(t *testing.T) {
+	// Full sphere = 24 cells at level 1 (6 faces * 4 children)
+	var cells CellUnion
+	for face := 0; face < NumFaces; face++ {
+		for _, child := range CellIDFromFace(face).Children() {
+			cells = append(cells, child)
+		}
+	}
+
+	// Remove one cell to create a hole
+	hole := cells[0]
+	cells = cells[1:]
+	cells.Normalize()
+
+	p, err := PolygonOfCellUnionBorder(cells)
+	if err != nil {
+		t.Fatalf("unexpected error: %v", err)
+	}
+
+	if p.NumLoops() != 2 {
+		t.Fatalf("NumLoops() = %d, want 2", p.NumLoops())
+	}
+
+	// The hole center should NOT be contained
+	holeCenter := CellFromCellID(hole).Center()
+	if p.ContainsPoint(holeCenter) {
+		t.Error("polygon contains the hole center, should not")
+	}
+
+	// A cell in the union should be contained
+	testCell := cells[(len(cells) / 2)]
+	testCenter := CellFromCellID(testCell).Center()
+	if !p.ContainsPoint(testCenter) {
+		t.Error("polygon does not contain a cell that should be in the union")
+	}
+
+	// Area should be close to full sphere minus one cell
+	expectedArea := FullPolygon().Area() - CellFromCellID(hole).ExactArea()
+	if math.Abs(p.Area()-expectedArea) > 0.001 {
+		t.Errorf("Area() = %.6f, want %.6f", p.Area(), expectedArea)
+	}
+}
+
+func TestPolygonOfCellUnionBorder_MixedLevels(t *testing.T) {
+	// Test that cells at different levels properly merge without internal
+	// boundary artifacts. This creates a region where a large cell (level 3)
+	// is adjacent to smaller cells (level 5).
+
+	// Start with a level 3 cell
+	largeCell := CellIDFromFace(0).ChildBeginAtLevel(3)
+
+	// Get an adjacent cell at level 3, then break it into level 5 children,
+	// but OMIT ONE to prevent Normalize() from merging them back.
+	adjacentLarge := largeCell.Next()
+	var smallCells CellUnion
+	for i, child := range adjacentLarge.Children() {
+		for j, grandchild := range child.Children() {
+			// Skip one grandchild to prevent normalization merging
+			if i == 0 && j == 0 {
+				continue
+			}
+			smallCells = append(smallCells, grandchild)
+		}
+	}
+
+	// Combine the large cell with the small cells
+	cells := CellUnion{largeCell}
+	cells = append(cells, smallCells...)
+	cells.Normalize()
+
+	// Verify we actually have mixed levels
+	minLevel, maxLevel := 30, 0
+	for _, cell := range cells {
+		level := cell.Level()
+		if level < minLevel {
+			minLevel = level
+		}
+		if level > maxLevel {
+			maxLevel = level
+		}
+	}
+	if minLevel == maxLevel {
+		t.Fatalf(
+			"test setup error: all cells at level %d, expected mixed levels",
+			minLevel,
+		)
+	}
+
+	p, err := PolygonOfCellUnionBorder(cells)
+	if err != nil {
+		t.Fatalf("unexpected error: %v", err)
+	}
+
+	if p.NumLoops() < 1 {
+		t.Fatal("expected at least 1 loop")
+	}
+
+	// All cell centers should be contained
+	for _, cellID := range cells {
+		center := CellFromCellID(cellID).Center()
+		if !p.ContainsPoint(center) {
+			t.Errorf("polygon does not contain cell %v center", cellID)
+		}
+	}
+}
+
+func TestPolygonOfCellUnionBorder_CrossFaceBoundary(t *testing.T) {
+	// Test that adjacent cells crossing a face boundary form a single polygon.
+	// Face 0 and Face 5 share an edge, so cells at their shared boundary
+	// should produce a merged polygon with 6 vertices (two squares sharing an edge).
+	cellA := CellIDFromFace(0).ChildBeginAtLevel(2)
+	cellB := CellIDFromFace(5).ChildEndAtLevel(2).Prev()
+
+	cells := CellUnion{cellA, cellB}
+	cells.Normalize()
+
+	p, err := PolygonOfCellUnionBorder(cells)
+	if err != nil {
+		t.Fatalf("unexpected error: %v", err)
+	}
+
+	if p.NumLoops() != 1 {
+		t.Fatalf("NumLoops() = %d, want 1", p.NumLoops())
+	}
+
+	// Two adjacent cells share one edge, so 4+4-2 = 6 vertices
+	if p.Loop(0).NumVertices() != 6 {
+		t.Errorf("NumVertices() = %d, want 6", p.Loop(0).NumVertices())
+	}
+}
+
+func TestPolygonOfCellUnionBorder_LTJunction(t *testing.T) {
+	// Test L-junction: 3 cells where two coarse cells meet one fine cell.
+	// Cell A (level 3), Cell B (level 3) adjacent to A, Cell C (level 9) also
+	// adjacent to A on a different edge. Internal edges should not leak.
+	cellA := CellIDFromFace(0).ChildBeginAtLevel(3)
+	cellB := cellA.EdgeNeighbors()[0]
+	cellC := cellA.EdgeNeighbors()[1]
+	for cellC.Level() < 9 {
+		cellC = cellC.ChildBegin()
+	}
+
+	cells := CellUnion{cellA, cellB, cellC}
+	cells.Normalize()
+
+	p, err := PolygonOfCellUnionBorder(cells)
+	if err != nil {
+		t.Fatalf("unexpected error: %v", err)
+	}
+
+	if p.NumLoops() > 2 {
+		t.Errorf("NumLoops() = %d, want <= 2", p.NumLoops())
+	}
+
+	for _, cellID := range cells {
+		center := CellFromCellID(cellID).Center()
+		if !p.ContainsPoint(center) {
+			t.Errorf("polygon does not contain cell %v center", cellID)
+		}
+	}
+}
+
+func TestPolygonOfCellUnionBorder_ManyLevelTransitions(t *testing.T) {
+	// Test cells at multiple levels (11-15) adjacent to each other, simulating
+	// real-world coastline data. The bug produced many spurious loops due to
+	// edge level mismatches.
+	var cells CellUnion
+
+	seed := CellIDFromFace(0).ChildBeginAtLevel(11)
+	cells = append(cells, seed)
+	for _, n := range seed.AllNeighbors(11) {
+		cells = append(cells, n)
+	}
+
+	// Add finer cells to create level transitions
+	fineArea := seed.EdgeNeighbors()[0]
+	for fineArea.Level() < 15 {
+		fineArea = fineArea.ChildBegin()
+	}
+	cells = append(cells, fineArea)
+	for _, n := range fineArea.AllNeighbors(15) {
+		cells = append(cells, n)
+	}
+
+	cells.Normalize()
+
+	p, err := PolygonOfCellUnionBorder(cells)
+	if err != nil {
+		t.Fatalf("unexpected error: %v", err)
+	}
+
+	if p.NumLoops() > 15 {
+		t.Errorf("NumLoops() = %d, want <= 15", p.NumLoops())
+	}
+}
+
 // TODO(roberts): Remaining Tests
 // TestInit
 // TestMultipleInit
@@ -1176,7 +1494,6 @@ func TestPolygonInvert(t *testing.T) {
 // TestBug1 - Bug14
 // TestPolylineIntersection
 // TestSplitting
-// TestInitToCellUnionBorder
 // TestUnionWithAmbgiuousCrossings
 // TestInitToSloppySupportsEmptyPolygons
 // TestInitToSnappedDoesNotRotateVertices