boolean.go

package polyclip

import (
	"errors"
	"fmt"
	"math"
	"sort"

	"github.com/lestrrat-go/polyclip/geom"
	"github.com/lestrrat-go/polyclip/internal/clip"
	"github.com/lestrrat-go/polyclip/internal/fixed"
)

// ErrHorizontalNotSupported is returned by the boolean ops when the
// bound-model pre-pass [clip.BuildLocalMinima] fails (typically because
// shared vertices between rings broke topology reconstruction) AND the
// legacy per-edge fallback's [clip.ClassifyHorizontals] then encounters a
// mid-bound horizontal it can't classify. The bound model itself handles
// mid-bound horizontals natively (DESIGN.md §12.10), so well-formed input
// rings without shared vertices never hit this path.
var ErrHorizontalNotSupported = errors.New("polyclip: input contains a horizontal edge that is neither a local minimum nor a local maximum of its ring")

// Union returns a ∪ b.
//
// Handled cases:
//
//   - Empty inputs: Union(empty, b) returns b unchanged, Union(a, empty)
//     returns a.
//   - Strictly disjoint bounding boxes: equivalent to concatenation. The
//     two MultiPolygons are returned spliced together with no engine work.
//   - Inputs with non-horizontal edges or with horizontal edges that are
//     each a local minimum (polygon bottom) or local maximum (polygon
//     top) of their ring: the Vatti engine in
//     [github.com/lestrrat-go/polyclip/internal/clip] runs over the snapped
//     segments. Output rings are converted back to a float64
//     MultiPolygon. Hole assignment uses signed-area sign and bbox-prefilter
//     point-in-polygon (DESIGN.md §11.9).
//
// Inputs containing a mid-bound horizontal (a staircase step) return
// [ErrHorizontalNotSupported] when the bound-model pre-pass fails on
// shared-vertex inputs that fall back to the per-edge path.
func Union(a, b geom.MultiPolygon) (geom.MultiPolygon, error) {
	return execOp(a, b, OpUnion, FillNonZero, nil)
}

// UnionAll returns the union of all inputs. It is functionally equivalent
// to repeated [Union], but pairs inputs in a tournament so the total work
// is O(n) Union calls of roughly balanced size rather than the O(n²)
// cumulative reduction `Union(Union(Union(p0, p1), p2), p3)…`.
//
// Empty input slice returns an empty [MultiPolygon]; a single-element
// slice returns that element unchanged.
func UnionAll(polys ...geom.MultiPolygon) (geom.MultiPolygon, error) {
	if len(polys) == 0 {
		return geom.MultiPolygon{}, nil
	}
	if len(polys) == 1 {
		return polys[0], nil
	}
	// Work on a local copy so the caller's slice isn't mutated when we
	// overwrite entries between rounds.
	current := make([]geom.MultiPolygon, len(polys))
	copy(current, polys)
	for len(current) > 1 {
		n := len(current)
		// Pair-merge in place: result of i and i+1 goes into slot i/2.
		// An unpaired trailing element survives to the next round.
		write := 0
		for i := 0; i+1 < n; i += 2 {
			merged, err := Union(current[i], current[i+1])
			if err != nil {
				return nil, err
			}
			current[write] = merged
			write++
		}
		if n%2 == 1 {
			current[write] = current[n-1]
			write++
		}
		current = current[:write]
	}
	return current[0], nil
}

// Intersect returns a ∩ b.
//
// Empty input or disjoint bounding boxes short-circuit to the empty
// MultiPolygon. Otherwise the Vatti engine runs with [clip.OpIntersect]
// and the §11.4 / §12.5 classification rules emit exactly the region
// covered by BOTH inputs.
func Intersect(a, b geom.MultiPolygon) (geom.MultiPolygon, error) {
	return execOp(a, b, OpIntersect, FillNonZero, nil)
}

// Difference returns a ∖ b — the region covered by a but not by b.
//
// Empty subject (a) short-circuits to empty; empty clip (b) returns a
// unchanged. Disjoint bounding boxes return a unchanged. Otherwise the
// Vatti engine runs with [clip.OpDifference].
func Difference(a, b geom.MultiPolygon) (geom.MultiPolygon, error) {
	return execOp(a, b, OpDifference, FillNonZero, nil)
}

// Xor returns the symmetric difference (a ∪ b) ∖ (a ∩ b) — the region
// covered by exactly one of the inputs.
//
// Empty operands short-circuit to the other input (or empty if both are
// empty). Disjoint bounding boxes return the concatenation, equivalent to
// Union. Otherwise the Vatti engine runs with [clip.OpXor].
func Xor(a, b geom.MultiPolygon) (geom.MultiPolygon, error) {
	return execOp(a, b, OpXor, FillNonZero, nil)
}

// Simplify resolves self-intersections and self-overlaps in m, returning an
// equivalent MultiPolygon whose rings are simple (non-self-intersecting) under
// the non-zero winding rule — the same convention the boolean ops use (CCW
// outer, CW hole). It runs the Vatti engine over m as a single source: a
// figure-eight splits into its two oppositely-wound loops, a ring traced twice
// collapses to one, and a doubled-back spur cancels.
//
// Simplify is the correct way to clean a self-intersecting polygon — for
// example the outward offset of a sharply concave ring. Running [Union] of m
// with itself does NOT work: identical inputs hit an idempotency
// short-circuit and are returned unchanged, leaving the self-intersection in
// place. [MultiPolygon.Clean] is purely geometric and also cannot resolve it.
//
// Empty input returns an empty MultiPolygon. Transversal (general-position)
// self-intersections resolve exactly. A snapped collinear degeneracy — e.g.
// parallel coincident walls landing on a single integer grid line — is subject
// to the same fixed-point limitation as the rest of the engine (DESIGN.md
// §7.2); callers needing robustness there should use [Offset], which adds its
// own multi-frame resolution.
func Simplify(m geom.MultiPolygon) (geom.MultiPolygon, error) {
	if len(m) == 0 {
		return geom.MultiPolygon{}, nil
	}
	bbox := m.BoundingBox()
	scale := fixed.ScaleFromBBox(bbox.Min.X, bbox.Min.Y, bbox.Max.X, bbox.Max.Y)

	segs := collectSegments(m, clip.Subject, scale, nil)
	return sweepSegments(segs, clip.OpUnion, clip.FillNonZero, scale, nil)
}

// mpolyEqual reports whether two MultiPolygons are deeply equal — same
// piece count, same outer ring vertices in the same order, same holes.
// Used by the idempotency short-circuits in [Union] / [Intersect] /
// [Difference] / [Xor] to bypass the engine when inputs are identical.
// Non-identical diff-src coincident inputs are resolved by the sweep's
// winding classification over first-class horizontal AEL edges (DESIGN.md
// §12.6.1); identical inputs are a degenerate case where every edge collapses
// to a same-vertex local minimum, which the bound-model pre-pass can't
// disambiguate.
func mpolyEqual(a, b geom.MultiPolygon) bool {
	if len(a) != len(b) {
		return false
	}
	for i := range a {
		if !polyEqual(a[i].Outer, b[i].Outer) {
			return false
		}
		if len(a[i].Holes) != len(b[i].Holes) {
			return false
		}
		for j := range a[i].Holes {
			if !polyEqual(a[i].Holes[j], b[i].Holes[j]) {
				return false
			}
		}
	}
	return true
}

func polyEqual(a, b geom.Polygon) bool {
	if len(a) != len(b) {
		return false
	}
	for i := range a {
		if a[i] != b[i] {
			return false
		}
	}
	return true
}

// sweepSegments runs the shared engine pipeline over source-tagged segments:
// split overlaps and T-junctions, dedup coincident pairs, sweep under op, then
// assemble the output rings back to a user-space MultiPolygon at scale (the same
// fixed-point scale used to build segs). It is the single seam shared by
// [runBooleanOp] and [Simplify]; op- and operand-specific post-filtering (the
// subset invariant) stays in the caller.
func sweepSegments(segs []clip.Segment, op clip.Operation, fill clip.FillRule, scale fixed.Scale, zt *zTracker) (geom.MultiPolygon, error) {
	segs = clip.Preprocess(segs)
	// SweepFillZ additionally records crossings for Z assignment; it is X/Y-
	// identical to SweepFill, so the standard (zt == nil) path is unchanged.
	var sw *clip.SweepResult
	if zt != nil {
		sw = clip.SweepFillZ(segs, op, fill)
	} else {
		sw = clip.SweepFill(segs, op, fill)
	}
	if sw.Err != nil {
		if errors.Is(sw.Err, clip.ErrUnsupportedHorizontal) {
			return nil, fmt.Errorf("%w: %v", ErrHorizontalNotSupported, sw.Err)
		}
		return nil, sw.Err
	}
	zt.applyCrossings(sw.Crossings, scale)
	return assembleResult(sw.Rings, scale, zt), nil
}

// runBooleanOp is the engine path: snap inputs to a fixed-point grid, feed
// segments through the sweep, and convert rings back to a user-space
// MultiPolygon.
func runBooleanOp(a, b geom.MultiPolygon, op clip.Operation, fill clip.FillRule, za ZAssigner) (geom.MultiPolygon, error) {
	bbox := a.BoundingBox().Union(b.BoundingBox())
	scale := fixed.ScaleFromBBox(bbox.Min.X, bbox.Min.Y, bbox.Max.X, bbox.Max.Y)

	zt := newZTracker(za)
	segs := collectSegments(a, clip.Subject, scale, zt)
	segs = append(segs, collectSegments(b, clip.Clip, scale, zt)...)

	res, err := sweepSegments(segs, op, fill, scale, zt)
	if err != nil {
		return nil, err
	}
	// Subset invariant: A∖B ⊆ A and A∩B ⊆ A∩B. The sweep can over-trace a
	// cross-source bound past where it exits the op-region (DESIGN.md §7.6),
	// emitting a spurious piece that lies OUTSIDE the required superset. Drop any
	// result piece whose interior point violates the op's subset invariant — a
	// correct piece always satisfies it, so this never drops valid output. (Only
	// whole entirely-spurious pieces are caught; a diagonal on an otherwise-valid
	// piece is not — those are handled in the sweep.)
	if op == clip.OpDifference || op == clip.OpIntersect {
		kept := res[:0]
		for _, ex := range res {
			pt, ok := interiorPoint(ex.Outer)
			// Only test hole-free pieces: interiorPoint(Outer) ignores holes and may
			// land inside one of ex's own holes, where a.Contains is false, wrongly
			// dropping a valid holed piece. Spurious over-traced fragments are simple
			// (hole-free), so this still catches them.
			if !ok || len(ex.Holes) > 0 {
				kept = append(kept, ex)
				continue
			}
			keep := a.Contains(pt) // Difference: ⊆ A
			if op == clip.OpIntersect {
				keep = a.Contains(pt) && b.Contains(pt)
			}
			if keep {
				kept = append(kept, ex)
			}
		}
		res = kept
	}
	return res, nil
}

// collectSegments converts every input edge into a fixed-point Segment and
// returns the slice. Horizontal segments are kept; the engine classifies
// them in a pre-pass.
func collectSegments(m geom.MultiPolygon, src clip.Source, scale fixed.Scale, zt *zTracker) []clip.Segment {
	var out []clip.Segment
	for _, ex := range m {
		appendRing(&out, ex.Outer, src, scale, true, zt)
		for _, h := range ex.Holes {
			appendRing(&out, h, src, scale, false, zt)
		}
	}
	return out
}

// appendRing emits ring's edges as fixed-point segments. The sweep's winding
// model derives each edge's WindDx from its direction (clip.signedContribution),
// so it assumes a canonical input orientation: CCW outers, CW holes. Inputs
// with the opposite winding invert WindDx for that source and misclassify
// (e.g. a CW subject made Union/Xor under-count). Normalize here by walking the
// ring in reverse when its signed area disagrees with wantCCW, so callers may
// pass either orientation.
//
// Collinear-through vertices (a vertex exactly on the straight line between its
// neighbours) are removed before emitting segments. They are geometrically
// redundant, but the bound model treats the extra segment's shared endpoint as
// a turn/maximum of its bound, so a flat-topped ring with a collinear vertex on
// the top edge mis-builds its rings (the collinear-mid-vertex degeneracy,
// DESIGN.md §12.11). Removal happens here, on the INPUT ring, before
// [clip.SplitOverlaps] / [clip.SplitTJunctions] introduce their own (intended)
// collinear split vertices at crossings.
func appendRing(dst *[]clip.Segment, ring geom.Polygon, src clip.Source, scale fixed.Scale, wantCCW bool, zt *zTracker) {
	n := len(ring)
	if n < 3 {
		return
	}
	reverse := (ring.SignedArea() < 0) == wantCCW
	pts := make([]fixed.Point, 0, n)
	for i := range n {
		k := i
		if reverse {
			k = n - 1 - i
		}
		p := scale.Snap(ring[k].X, ring[k].Y)
		// Record this input vertex's Z under its grid point so the output can
		// carry it back (no-op when Z tracking is off). Done before the dedup
		// skip so every input vertex contributes.
		zt.recordInput(p, ring[k].Z)
		if len(pts) > 0 && pts[len(pts)-1] == p {
			continue
		}
		pts = append(pts, p)
	}
	// Drop a wrap-around duplicate (first == last) so simplifyCollinearRing never
	// sees a zero-length neighbour, which Orient2D reads as collinear and would
	// wrongly delete a real corner of a ring with a repeated vertex.
	for len(pts) >= 2 && pts[0] == pts[len(pts)-1] {
		pts = pts[:len(pts)-1]
	}
	pts = simplifyCollinearRing(pts)
	m := len(pts)
	if m < 3 {
		return
	}
	for i := range m {
		j := i + 1
		if j == m {
			j = 0
		}
		seg := clip.NewSegment(pts[i], pts[j], src)
		if seg.Degenerate() {
			continue
		}
		*dst = append(*dst, seg)
	}
}

// simplifyCollinearRing removes vertices that lie exactly on the straight line
// between their cyclic neighbours from a closed ring of grid points. For a
// simple (non-self-intersecting) polygon every collinear vertex is a redundant
// through-vertex, so removal is an exact geometric no-op. The pass repeats to
// collapse runs of three or more collinear vertices down to the two endpoints.
// Consecutive duplicate points (including the wrap) are assumed already removed
// by the caller.
func simplifyCollinearRing(pts []fixed.Point) []fixed.Point {
	for len(pts) >= 3 {
		n := len(pts)
		kept := make([]fixed.Point, 0, n)
		for i := range n {
			prev := pts[(i-1+n)%n]
			next := pts[(i+1)%n]
			if fixed.Orient2D(prev, pts[i], next) == 0 {
				continue
			}
			kept = append(kept, pts[i])
		}
		if len(kept) == n {
			return pts
		}
		pts = kept
	}
	return pts
}

// classifiedRing is an output ring prepared for containment classification: its
// polygon, bounding box, absolute area, and a precomputed interior point.
type classifiedRing struct {
	poly geom.Polygon
	bbox geom.BBox
	area float64 // absolute (unsigned) area
	// interior is a point of poly's OPEN interior, precomputed once per ring
	// (it depends only on poly). hasInterior is false for a degenerate ring
	// with no interior. Hoisted out of ringContains so the O(n^2) containment
	// scan does not recompute it on every candidate container.
	interior    geom.Point
	hasInterior bool
}

// newClassifiedRing precomputes the classification fields for poly.
func newClassifiedRing(poly geom.Polygon) classifiedRing {
	interior, hasInterior := interiorPoint(poly)
	return classifiedRing{
		poly: poly, bbox: poly.BoundingBox(),
		area:        math.Abs(poly.SignedArea()),
		interior:    interior,
		hasInterior: hasInterior,
	}
}

// ringContains reports whether inner is nested within outer. The sweep's
// output rings have pairwise-disjoint interiors (they partition the plane
// into in/out), so two rings are either disjoint or one strictly contains
// the other — never partially overlapping. Under that invariant, inner is
// nested in outer iff a point of inner's OPEN interior lies inside outer.
//
// The sample must be a genuine interior point of inner, not a vertex or the
// vertex centroid. When two rings merely touch, their shared vertices — and,
// for a collinear shared edge, even the vertex centroid — land ON the other
// ring's boundary, which Polygon.Contains counts as inside, wrongly nesting
// polygons that only touch (the shared-vertex bug, DESIGN.md §12.11).
// Conversely a hole emitted by the sweep can have ALL its vertices on the
// enclosing outer's boundary (e.g. the Xor overlap rectangle whose corners
// sit on the union outline), so a vertex-based test gives the opposite false
// negative. An interior point of inner avoids both: if inner is nested it is
// strictly inside outer; if the rings only touch it is strictly outside.
func ringContains(inner, outer classifiedRing) bool {
	if !inner.hasInterior {
		return false
	}
	return outer.bbox.Contains(inner.interior) && outer.poly.Contains(inner.interior)
}

// canContainRing reports whether outer may be inner's container. A strictly
// larger ring can. Two coincident rings have EQUAL area and mutually contain
// each other (same boundary; e.g. Difference/Xor of identical inputs emits the
// region once CCW and once CW, which must cancel to zero area): the tie is
// broken by orientation — the filled (CCW) ring is the parent, the hole (CW)
// ring the child — so the pair nests as outer+hole and cancels.
func canContainRing(outer, inner classifiedRing) bool {
	if outer.area > inner.area {
		return true
	}
	return outer.area == inner.area &&
		outer.poly.SignedArea() > 0 && inner.poly.SignedArea() < 0
}

// buildContainmentForest computes, for each ring, the index of its immediate
// (smallest) container, or -1 if it is top-level. The rings' interiors are
// pairwise disjoint, so any two are either disjoint or one strictly contains
// the other. The sweep's own CCW/CW orientation is NOT used to classify
// outer-vs-hole: it is locally meaningful but does not encode nesting depth
// (an island inside a hole is CCW yet must become a filled top-level piece;
// an Intersect cycle can emit a sole region CW). Depth parity over this forest
// is the only reliable signal — see DESIGN.md §11.9 / §12.10.
//
// This is a stabbing query: a true container's bbox must contain the query
// ring's interior point (ringContains's first conjunct), so its X-interval
// [Min.X, Max.X] covers the query's interior.X. Candidates are therefore the
// rings whose Min.X <= query.X (a prefix of an order sorted by Min.X, found by
// binary search) AND whose Max.X >= query.X. This prune is exact — it never
// drops a real container — so it shrinks the candidate set without altering the
// smallest-container selection (DESIGN.md §7.10). Worst case (every X-interval
// overlaps) is still O(n²); the common disjoint slicer output collapses to
// near-linear.
func buildContainmentForest(rings []classifiedRing) []int {
	n := len(rings)
	parent := make([]int, n)

	// order: ring indices sorted ascending by bbox.Min.X; minXs is the parallel
	// Min.X array so the binary search avoids an indirection per probe.
	// Degenerate rings stay in order — they are filtered in the scan, not here,
	// so the prefix logic stays simple.
	order := make([]int, n)
	for i := range order {
		order[i] = i
	}
	sort.Slice(order, func(a, b int) bool {
		return rings[order[a]].bbox.Min.X < rings[order[b]].bbox.Min.X
	})
	minXs := make([]float64, n)
	for k, idx := range order {
		minXs[k] = rings[idx].bbox.Min.X
	}

	for i := range rings {
		parent[i] = -1
		ri := rings[i]
		// No interior point to stab with: ringContains would reject it as inner
		// anyway, so it can have no container.
		if len(ri.poly) == 0 || !ri.hasInterior {
			continue
		}
		qx := ri.interior.X

		// Candidate prefix order[0:hi) holds every ring with Min.X <= qx; one
		// with Min.X > qx cannot contain ri.interior on X. sort.Search returns
		// the first index whose minXs > qx, i.e. the exclusive prefix end.
		hi := sort.Search(n, func(k int) bool { return minXs[k] > qx })

		for k := range hi {
			j := order[k]
			if j == i || len(rings[j].poly) == 0 {
				continue
			}
			// Complete the X-interval stab (Min.X <= qx holds by the prefix).
			if rings[j].bbox.Max.X < qx {
				continue
			}
			if !canContainRing(rings[j], ri) {
				continue
			}
			if !ringContains(ri, rings[j]) {
				continue
			}
			// Prefer the smallest container; among equal-area coincident
			// containers any is fine (there is normally just one).
			if parent[i] == -1 || rings[j].area < rings[parent[i]].area {
				parent[i] = j
			}
		}
	}
	return parent
}

// ringDepth returns the number of rings strictly containing ring i, walking the
// parent chain of a forest produced by [buildContainmentForest].
func ringDepth(parent []int, i int) int {
	d := 0
	for parent[i] != -1 {
		i = parent[i]
		d++
	}
	return d
}

// classifyOutRecs turns the sweep's closed output rings into classifiedRings:
// dedup consecutive points, split self-touching cycles into simple loops, and
// unsnap to user space at scale.
func classifyOutRecs(rings []*clip.OutRec, scale fixed.Scale, zt *zTracker) []classifiedRing {
	var rings2 []classifiedRing
	for _, r := range rings {
		if r.Pts == nil {
			continue
		}
		base := dedupConsecutive(r.Points())
		if len(base) < 3 {
			continue
		}
		// A ring that revisits a vertex is self-touching — two boundary loops
		// meeting at a shared point, produced when an input vertex lies on the
		// other source's edge (vertex-on-edge degeneracy) or when the sweep
		// merges rings at a same-side maximum (AddLocalMaxPoly's figure-8 pinch).
		// Decompose into simple loops so each is classified independently; a
		// single self-touching cycle would yield a wrong net shoelace area.
		// Simple rings have no repeats and pass through unchanged.
		for _, fixedPts := range splitSelfTouchingRings(base) {
			if len(fixedPts) < 3 {
				continue
			}
			poly := make(geom.Polygon, len(fixedPts))
			for i, fp := range fixedPts {
				poly[i].X, poly[i].Y = scale.Unsnap(fp)
				poly[i].Z = zt.lookup(fp)
			}
			rings2 = append(rings2, newClassifiedRing(poly))
		}
	}
	return rings2
}

// assembleResult converts the sweep's closed output rings into a user-space
// MultiPolygon, classifying each ring as outer or hole by its signed area
// and grouping holes into their containing outer.
func assembleResult(rings []*clip.OutRec, scale fixed.Scale, zt *zTracker) geom.MultiPolygon {
	rings2 := classifyOutRecs(rings, scale, zt)
	parent := buildContainmentForest(rings2)

	// depth(i) = number of rings strictly containing i. Even depth = a filled
	// region (the outer of an ExPolygon); odd depth = a hole. A hole's parent
	// has depth one less, so it is always an even-depth filled ring — attach
	// the hole directly to it. A filled ring at depth ≥ 2 (an island inside a
	// hole) is a top-level ExPolygon of its own (MultiPolygon is flat).
	depth := func(i int) int {
		return ringDepth(parent, i)
	}

	idxMap := make(map[int]int, len(rings2))
	var result geom.MultiPolygon
	for i := range rings2 {
		if len(rings2[i].poly) == 0 || depth(i)%2 != 0 {
			continue
		}
		if rings2[i].poly.SignedArea() < 0 { // a filled ring must be CCW
			rings2[i].poly.Reverse()
		}
		idxMap[i] = len(result)
		result = append(result, geom.ExPolygon{Outer: rings2[i].poly})
	}
	for i := range rings2 {
		if len(rings2[i].poly) == 0 || depth(i)%2 == 0 {
			continue
		}
		if rings2[i].poly.SignedArea() > 0 { // a hole must be CW
			rings2[i].poly.Reverse()
		}
		if k, ok := idxMap[parent[i]]; ok {
			result[k].Holes = append(result[k].Holes, rings2[i].poly)
		}
	}

	return result
}

// buildPolyTree reconstructs the full nested containment hierarchy of m as a
// [PolyTree]. A MultiPolygon is flat — an island inside a hole is a top-level
// ExPolygon — so the depth-≥2 nesting is recovered by rebuilding the
// containment forest (§11.9) over every ring of m (each outer and each hole).
// Operating on the finished MultiPolygon, rather than the sweep's raw rings,
// lets ExecuteTree reuse every execOp path (short-circuits, Xor composition,
// per-piece Difference, alternate fills) unchanged: each yields a well-formed
// MultiPolygon whose rings have pairwise-disjoint interiors, exactly the
// invariant the forest needs.
//
// Each ring becomes a node; a node's Children are the rings whose immediate
// container is it. Even depth = a filled region, odd depth = a hole; winding is
// normalized as elsewhere (filled CCW, hole CW). Top-level nodes (no container)
// are the tree's roots.
func buildPolyTree(m geom.MultiPolygon) *PolyTree {
	var rings []classifiedRing
	for _, ex := range m {
		if len(ex.Outer) >= 3 {
			rings = append(rings, newClassifiedRing(ex.Outer))
		}
		for _, h := range ex.Holes {
			if len(h) >= 3 {
				rings = append(rings, newClassifiedRing(h))
			}
		}
	}

	parent := buildContainmentForest(rings)

	nodes := make([]*PolyTreeNode, len(rings))
	for i := range rings {
		isHole := ringDepth(parent, i)%2 == 1
		// Copy so the tree owns its rings: execOp short-circuits (e.g.
		// Union(s,∅)→s) can return the caller's own polygons, which Reverse must
		// not mutate. For a well-formed MultiPolygon the orientation already
		// matches depth parity, so the reversal below is normally a no-op.
		poly := append(geom.Polygon(nil), rings[i].poly...)
		switch {
		case isHole && poly.SignedArea() > 0: // a hole must be CW
			poly.Reverse()
		case !isHole && poly.SignedArea() < 0: // a filled ring must be CCW
			poly.Reverse()
		}
		nodes[i] = &PolyTreeNode{Polygon: poly, IsHole: isHole}
	}

	tree := &PolyTree{}
	for i := range rings {
		if parent[i] == -1 {
			tree.Children = append(tree.Children, nodes[i])
			continue
		}
		p := parent[i]
		nodes[p].Children = append(nodes[p].Children, nodes[i])
	}
	return tree
}

// dedupConsecutive removes consecutive identical points from a closed ring,
// including the wrap-around (last == first). The sweep can emit a vertex
// twice at a maxima confluence — once when a hot maxima edge is crossed past
// a cold co-maximum edge (the one-hot IntersectEdges branch adds the apex on
// one ring side) and again when AddLocalMaxPoly closes the ring on the other
// side — leaving a zero-length edge. Clipper2 strips these in its output
// stage (BuildPath); this is the equivalent cleanup.
// splitSelfTouchingRings decomposes a closed ring whose vertex list repeats a
// vertex into simple loops. Walking the vertices, whenever the walk returns to
// a vertex already on the open path, the run since that vertex forms a closed
// sub-loop and is split off; the shared vertex stays on the continuing path.
// The result loops are each simple (no internal repeats). A ring with no
// repeated vertex returns as a single loop unchanged.
func splitSelfTouchingRings(pts []fixed.Point) [][]fixed.Point {
	var loops [][]fixed.Point
	stack := make([]fixed.Point, 0, len(pts))
	at := make(map[fixed.Point]int, len(pts))
	for _, p := range pts {
		if j, ok := at[p]; ok {
			loop := make([]fixed.Point, len(stack)-j)
			copy(loop, stack[j:])
			loops = append(loops, loop)
			for k := j; k < len(stack); k++ {
				delete(at, stack[k])
			}
			stack = stack[:j]
		}
		stack = append(stack, p)
		at[p] = len(stack) - 1
	}
	if len(stack) > 0 {
		loops = append(loops, stack)
	}
	return loops
}

func dedupConsecutive(pts []fixed.Point) []fixed.Point {
	if len(pts) < 2 {
		return pts
	}
	out := make([]fixed.Point, 0, len(pts))
	out = append(out, pts[0])
	for _, p := range pts[1:] {
		if p == out[len(out)-1] {
			continue
		}
		out = append(out, p)
	}
	for len(out) >= 2 && out[len(out)-1] == out[0] {
		out = out[:len(out)-1]
	}
	return out
}

// interiorPoint returns a point strictly inside the simple polygon p, and a
// bool reporting success. It casts a horizontal ray through the polygon's
// vertex-Y centroid, collects the edge crossings, and returns the midpoint of
// the widest interior span. That point is guaranteed strictly inside p (it sits
// between an entering and leaving crossing of a well-formed ring), independent
// of whether p is convex — unlike the vertex centroid, which can fall outside a
// concave ring or land on a neighbour's boundary. Used by nesting detection so
// that two rings which merely touch are never reported as nested (DESIGN.md
// §12.11). Returns false for degenerate rings (<3 vertices, or no interior span
// found, e.g. zero area).
func interiorPoint(p geom.Polygon) (geom.Point, bool) {
	n := len(p)
	if n < 3 {
		return geom.Point{}, false
	}
	// Choose the scanline Y strictly between two adjacent distinct vertex Ys
	// (the widest gap), so it grazes no vertex and runs along no horizontal
	// edge. A Y equal to a vertex — e.g. the mean coinciding with a horizontal
	// edge the ring shares with another — makes the "interior" span run along
	// the ring's own boundary, returning a boundary point that Polygon.Contains
	// treats as inside the other ring and wrongly nests touching polygons
	// (DESIGN.md §12.11).
	ys := make([]float64, n)
	for i, v := range p {
		ys[i] = v.Y
	}
	sort.Float64s(ys)
	gapLo, gap := 0.0, 0.0
	for i := 0; i+1 < n; i++ {
		if g := ys[i+1] - ys[i]; g > gap {
			gap, gapLo = g, ys[i]
		}
	}
	if gap <= 0 {
		return geom.Point{}, false // degenerate: all vertices share one Y
	}
	y := gapLo + gap/2

	var xs []float64
	for i := range n {
		a := p[i]
		b := p[(i+1)%n]
		// Half-open [min.Y, max.Y) crossing test: counts each edge once and
		// avoids double-counting at shared vertices.
		if (a.Y <= y) == (b.Y <= y) {
			continue
		}
		t := (y - a.Y) / (b.Y - a.Y)
		xs = append(xs, a.X+t*(b.X-a.X))
	}
	if len(xs) < 2 {
		return geom.Point{}, false
	}
	sort.Float64s(xs)

	// Interior spans are the (0,1),(2,3),… pairs. Pick the widest so the
	// midpoint sits well clear of both boundaries.
	bestLo, bestHi, bestW := 0.0, 0.0, -1.0
	for i := 0; i+1 < len(xs); i += 2 {
		if w := xs[i+1] - xs[i]; w > bestW {
			bestW, bestLo, bestHi = w, xs[i], xs[i+1]
		}
	}
	if bestW <= 0 {
		return geom.Point{}, false
	}
	return geom.Point{X: (bestLo + bestHi) / 2, Y: y}, true
}