diff --git a/CHANGELOG.md b/CHANGELOG.md
index 547b3edd6..4314f0a7a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,10 @@ The format of this changelog is based on
 [Keep a Changelog](https://keepachangelog.com/), and this project adheres to
 [Semantic Versioning](https://semver.org/).
 
+## Unreleased
+
+  - Added `recover_curves` and the curve-preserving Boolean variants `union2d_curved`, `difference2d_curved`, `intersect2d_curved`, and `xor2d_curved`, which recover original curves (arcs, splines) from a clipped result wherever their discretized footprint survived the operation intact, returning a `Vector{CurvilinearRegion}`
+
 ## 1.14.0 (2026-05-28)
 
   - Added `ParameterSet`, a nested dictionary wrapper with dot-access for reading and
diff --git a/docs/src/concepts/polygons.md b/docs/src/concepts/polygons.md
index a42a3f671..39a8d5cea 100644
--- a/docs/src/concepts/polygons.md
+++ b/docs/src/concepts/polygons.md
@@ -21,7 +21,7 @@ Geometric Boolean operations on polygons are called "clipping" operations. For 2
 
 !!! info
 
-    Boolean operations in 3D with `SolidModel` are handled by the Open CASCADE Technology kernel, which works directly with rich geometry types rendered from our native `CoordinateSystem`. If you need boolean operations involving curved geometry whose results can't be achieved by clipping-then-rounding, then your 2D geometry should defer the boolean operation until `SolidModel` postrendering so that the result will still be represented with curves.
+    Boolean operations in 3D with `SolidModel` are handled by the Open CASCADE Technology kernel, which works directly with rich geometry types rendered from our native `CoordinateSystem`. If you need boolean operations involving curved geometry whose results can't be achieved by clipping-then-rounding, you have two options: keep curves in 2D using the curve-preserving Boolean variants (see [Recovering curves through clipping](@ref) below), or defer the boolean operation until `SolidModel` postrendering so that the result will still be represented with curves.
 
 For many use cases, `union2d`, `difference2d`, `intersect2d`, and `xor2d` behave as expected and are easiest to use.
 More general operations may be accomplished using the `clip` function.
@@ -45,6 +45,34 @@ A [`CurvilinearRegion`](@ref) pairs a `CurvilinearPolygon` exterior with zero or
 
 See [API Reference: Curvilinear geometry](@ref api-curvilinear).
 
+### Recovering curves through clipping
+
+Boolean operations (`union2d`, `difference2d`, etc.) discretize curved geometry to polygons
+before passing them to the Clipper library. Normally, the original curves (arcs, splines)
+are lost in this process. The [`recover_curves`](@ref) function and its convenience
+wrappers (`difference2d_curved`, `union2d_curved`, `intersect2d_curved`, `xor2d_curved`)
+track each input curve's discretized integer-grid footprint and substitute the original
+curve back into the result wherever that footprint survived the boolean operation intact.
+
+The curve-preserving variants return a `Vector{CurvilinearRegion}` rather than a single
+`ClippedPolygon`. Each region in the vector corresponds to one outer contour in the clipped
+result (the disjoint pieces each become a separate region). For example:
+
+```julia
+# Standard clipping discretizes curves to polygons:
+result = difference2d(a, b)  # ClippedPolygon
+
+# Curve-preserving variant recovers arcs where possible:
+regions = difference2d_curved(a, b)  # Vector{CurvilinearRegion}, arcs preserved
+```
+
+**Current limitations:** A curve is recovered only if its entire discretized run survives
+the boolean operation with exact integer equality. If the operation cuts through a curve
+(e.g., a straight edge crossing an arc's interior), that curve falls back to a polyline.
+Additionally, curves can only be recovered on CurvilinearRegion/CurvilinearPolygon,
+Path nodes, and `Rounded`-styled `Polygon`/`Rectangle` entities. `Rounded` applied to
+other entities, styled Curvilinear entities, and nested styles do not yet support curve recovery.
+
 ## Styles
 
 In addition to other generic [entity styles](./geometry.md#Entity-Styles) like `NoRender`, `AbstractPolygon`s can be paired with the `Rounded` style. `ClippedPolygon`s support `StyleDict`, which allows for different styles to be applied to different contours in its tree.
diff --git a/docs/src/reference/api.md b/docs/src/reference/api.md
index 3eda018c4..a5b9c51dc 100644
--- a/docs/src/reference/api.md
+++ b/docs/src/reference/api.md
@@ -194,6 +194,11 @@
     CurvilinearRegion
     Curvilinear.edge_type_at_vertex
     Curvilinear.line_arc_cornerindices
+    difference2d_curved
+    intersect2d_curved
+    recover_curves
+    union2d_curved
+    xor2d_curved
 ```
 
 ### [Shapes](@id api-shapes)
diff --git a/src/DeviceLayout.jl b/src/DeviceLayout.jl
index 2fef9d2f7..8f1697550 100644
--- a/src/DeviceLayout.jl
+++ b/src/DeviceLayout.jl
@@ -563,7 +563,15 @@ include("render/render.jl")
 # After render.jl to ensure access to `to_polygons`
 include("curvilinear.jl")
 using .Curvilinear
-export Curvilinear, CurvilinearPolygon, CurvilinearRegion, pathtopolys
+export Curvilinear,
+    CurvilinearPolygon,
+    CurvilinearRegion,
+    pathtopolys,
+    recover_curves,
+    difference2d_curved,
+    union2d_curved,
+    intersect2d_curved,
+    xor2d_curved
 
 include("simple_shapes.jl")
 import .SimpleShapes:
diff --git a/src/curve_recovery.jl b/src/curve_recovery.jl
new file mode 100644
index 000000000..325a2a1fe
--- /dev/null
+++ b/src/curve_recovery.jl
@@ -0,0 +1,598 @@
+# Provenance-based curve recovery.
+# Included in Curvilinear — needs CurvilinearRegion, discretize_curve, and the Polygons clip functions.
+
+using .Polygons: clipperize, ClippedPolygon, union2d, difference2d, intersect2d, xor2d
+
+# One discretized curve's footprint on the integer grid, paired with its source segment.
+# The integer-grid coordinate type P matches whatever Polygons.clipperize(::Point{R}) produces.
+struct ProvenanceRun{P}
+    run::Vector{P}
+    curve::Paths.Segment
+end
+
+# Returns (polys::Vector{Polygon{R}}, runs::Vector{ProvenanceRun}) for one operand.
+# `R` is the promoted coordinate type clip will use; discretization uses clip's default atol.
+function discretize_with_provenance(entities, ::Type{R}) where {R}
+    polys = Polygon{R}[]
+    runs = ProvenanceRun[]
+    atol = DeviceLayout.onenanometer(R)
+    for ent in entities
+        _collect_provenance!(polys, runs, ent, R, atol)
+    end
+    return polys, runs
+end
+
+# CurvilinearPolygon: walk exactly like to_polygons(::CurvilinearPolygon), but also
+# record each curve's full [start, interior…, end] run, snapped to the integer grid.
+function _collect_provenance!(polys, runs, e::CurvilinearPolygon, ::Type{R}, atol) where {R}
+    ec = convert(CurvilinearPolygon{R}, e)
+    i = 1
+    p = Point{R}[]
+    # Walk in ascending start-index order — must match to_polygons(::CurvilinearPolygon) exactly
+    # so the Clipper polygon and these provenance runs share the same point ordering. See the
+    # note there: an out-of-order (e.g. wrap-seam-first) curve list otherwise emits the wrong
+    # `ec.p[i:csi]` spans and a sub-µm near-pinch.
+    order =
+        issorted(ec.curve_start_idx) ? eachindex(ec.curve_start_idx) :
+        sortperm(ec.curve_start_idx)
+    for idx in order
+        csi = ec.curve_start_idx[idx]
+        c = ec.curves[idx]
+        append!(p, ec.p[i:csi])
+        wrapped_end = mod1(csi + 1, length(ec.p))
+        pp = DeviceLayout.discretize_curve(c, atol; rtol=nothing)
+        run_pts = vcat([ec.p[csi]], pp[2:(end - 1)], [ec.p[wrapped_end]])
+        push!(runs, ProvenanceRun(clipperize.(run_pts), c))
+        append!(p, pp[2:(end - 1)])
+        i = csi + 1
+    end
+    append!(p, ec.p[i:end])
+    push!(polys, Polygon{R}(p))
+    return nothing
+end
+
+# CurvilinearRegion: exterior + holes, each a CurvilinearPolygon.
+function _collect_provenance!(polys, runs, e::CurvilinearRegion, ::Type{R}, atol) where {R}
+    _collect_provenance!(polys, runs, e.exterior, R, atol)
+    for h in e.holes
+        _collect_provenance!(polys, runs, _reverse(h), R, atol)
+    end
+    return nothing
+end
+
+# Count, in already-clipperized integer space, runs of short edges that fit a circle — a cheap
+# proxy for "this discretized polygon CONTAINS curve geometry we are about to lose". Used only to
+# warn on the silent-discretization path below; not exact, just enough to flag a likely loss.
+function _count_arclike_runs(pts; short_nm=6000.0, min_run=4, resid_nm=2.0)
+    n = length(pts)
+    n < min_run + 1 && return 0
+    xy = [
+        (
+            Float64(DeviceLayout.ustrip(getx(p))) / 1000,
+            Float64(DeviceLayout.ustrip(gety(p))) / 1000
+        ) for p in pts
+    ]  # → ~nm scale-agnostic
+    # edge lengths in the same units as xy
+    el = [
+        hypot(xy[mod1(i + 1, n)][1] - xy[i][1], xy[mod1(i + 1, n)][2] - xy[i][2]) for
+        i = 1:n
+    ]
+    short = short_nm / 1000
+    cnt = 0
+    i = 1
+    while i <= n
+        if !(0 < el[i] < short)
+            i += 1
+            continue
+        end
+        j = i
+        while j < n && 0 < el[mod1(j + 1, n)] < short
+            j += 1
+        end
+        rl = j - i + 1
+        if rl >= min_run
+            vx = [xy[mod1(i + k, n)][1] for k = 0:rl]
+            vy = [xy[mod1(i + k, n)][2] for k = 0:rl]
+            m = length(vx)
+            sx = sum(vx)
+            sy = sum(vy)
+            sxx = sum(vx .^ 2)
+            syy = sum(vy .^ 2)
+            sxy = sum(vx .* vy)
+            sxxx = sum(vx .^ 3)
+            syyy = sum(vy .^ 3)
+            sxyy = sum(vx .* vy .^ 2)
+            sxxy = sum(vx .^ 2 .* vy)
+            A = m * sxx - sx^2
+            B = m * sxy - sx * sy
+            Cc = m * syy - sy^2
+            D = 0.5 * (m * sxxx + m * sxyy - sx * sxx - sx * syy)
+            E = 0.5 * (m * syyy + m * sxxy - sy * syy - sy * sxx)
+            den = A * Cc - B^2
+            if abs(den) > 1e-12
+                cx = (D * Cc - B * E) / den
+                cy = (A * E - B * D) / den
+                r = sqrt(max(0.0, (sxx + syy - 2cx * sx - 2cy * sy) / m + cx^2 + cy^2))
+                if 1e-3 < r < 1e5
+                    mr = maximum(abs(hypot(vx[k] - cx, vy[k] - cy) - r) for k = 1:m)
+                    mr < resid_nm / 1000 && (cnt += 1)
+                end
+            end
+        end
+        i = j + 1
+    end
+    return cnt
+end
+
+# Module-level switch for the silent-discretization warning (default ON; set false to silence).
+const warn_on_curve_loss = Ref(true)
+# Track per-(entity-type) loss counts so a run can report which entity classes lost curves.
+const _curve_loss_log = Dict{String, Int}()
+curve_loss_log() = _curve_loss_log
+reset_curve_loss_log!() = empty!(_curve_loss_log)
+
+# Any other entity: no curves to recover; discretize to polygons. THIS is the silent-loss path —
+# a curve-bearing entity whose (type, style) combination has no `_as_entities`/`_collect_provenance!`
+# method falls here and is discretized to polyline with NO provenance, so its curves can never be
+# recovered downstream. We detect the likely-loss case (the discretized polygon contains short-edge
+# runs that fit a circle) and warn once per entity type, and tally it in `_curve_loss_log`, so such
+# losses are observable instead of silent. (Plain rectilinear polygons fit nothing → no warning.)
+function _collect_provenance!(polys, runs, e, ::Type{R}, atol) where {R}
+    # Convert to polygons and append. to_polygons returns a polygon or array.
+    poly_result = DeviceLayout.to_polygons(e)
+    polylist = poly_result isa Polygon ? (poly_result,) : poly_result
+    arclike = 0
+    for poly in polylist
+        push!(polys, convert(Polygon{R}, poly))
+        warn_on_curve_loss[] && (arclike += _count_arclike_runs(points(poly)))
+    end
+    if warn_on_curve_loss[] && arclike > 0
+        key = string(nameof(typeof(e)))
+        e isa StyledEntity && (
+            key =
+                "Styled{" *
+                string(nameof(typeof(e.ent))) *
+                "," *
+                string(nameof(typeof(e.sty))) *
+                "}"
+        )
+        first_seen = !haskey(_curve_loss_log, key)
+        _curve_loss_log[key] = get(_curve_loss_log, key, 0) + arclike
+        first_seen &&
+            @warn "recover_curves: discretizing a curve-bearing entity with no provenance — its arcs are LOST. " *
+                  "Add an _as_entities method for this (type, style) to recover them." entity_type =
+                key arclike_runs = arclike
+    end
+    return nothing
+end
+
+# Search `contour` (treated as cyclic) for `run` as a contiguous block, forward or reversed.
+# Returns (start, reversed) of the first hit (1-based start index into contour), or nothing.
+# Exact integer equality — no tolerance.
+function match_run(contour::AbstractVector{P}, run::AbstractVector{P}) where {P}
+    n = length(contour)
+    m = length(run)
+    (m == 0 || m > n) && return nothing
+    rev = reverse(run)
+    for s = 1:n
+        fwd_ok = true
+        rev_ok = true
+        for k = 0:(m - 1)
+            cv = contour[mod1(s + k, n)]
+            fwd_ok &= (cv == run[k + 1])
+            rev_ok &= (cv == rev[k + 1])
+            (fwd_ok || rev_ok) || break
+        end
+        fwd_ok && return (start=s, reversed=false)
+        rev_ok && return (start=s, reversed=true)
+    end
+    return nothing
+end
+
+# Find ALL maximal contiguous sub-blocks of `run` that exist in `contour` (treated as cyclic),
+# considering forward and reversed orientations. Returns a Vector of named-tuples
+# `(contour_start, run_lo, run_hi, reversed)` where:
+#   - `contour_start` is the 1-based start position in the contour
+#   - `run_lo`, `run_hi` are 1-based inclusive indices into `run` (or `reverse(run)` if reversed)
+#   - `reversed::Bool` indicates orientation
+# Each returned block has length ≥ `min_len` (default 3 — meaningful sub-arc, not a single edge).
+# Sub-blocks that are subsumed by a larger block at the same contour position are dropped.
+#
+# Used by `substitute_curves` for partial recovery: when a curve is "genuinely cut" by a
+# boolean operation (e.g. a Turn arc bisected by another polygon's edge), the run survives
+# as multiple disjoint contiguous sub-blocks. Each sub-block becomes a sub-curve via the
+# original segment's `Paths.split` interface.
+function match_run_partial(
+    contour::AbstractVector{P},
+    run::AbstractVector{P};
+    min_len::Int=3
+) where {P}
+    n = length(contour)
+    m = length(run)
+    (m == 0 || n == 0) && return NamedTuple{
+        (:contour_start, :run_lo, :run_hi, :reversed),
+        Tuple{Int, Int, Int, Bool}
+    }[]
+    # Build a position lookup: for each run point value, what indices in run have that value?
+    # Run is short (typically 50-200 points for an arc), contour can be longer; this trades a
+    # bit of memory for fast search.
+    run_pos = Dict{P, Vector{Int}}()
+    for (i, v) in enumerate(run)
+        push!(get!(run_pos, v, Int[]), i)
+    end
+    rev_run = reverse(run)
+    rev_pos = Dict{P, Vector{Int}}()
+    for (i, v) in enumerate(rev_run)
+        push!(get!(rev_pos, v, Int[]), i)
+    end
+
+    matches = NamedTuple{
+        (:contour_start, :run_lo, :run_hi, :reversed),
+        Tuple{Int, Int, Int, Bool}
+    }[]
+
+    # Try to extend a contiguous run-match starting at contour position s, run position r0,
+    # in the given orientation. Returns the longest extension length k such that
+    # contour[s..s+k-1] == run[r0..r0+k-1] (orientation-respecting).
+    extend_match(s::Int, r0::Int, src::AbstractVector{P}) = begin
+        m_local = length(src)
+        k = 0
+        while r0 + k <= m_local && contour[mod1(s + k, n)] == src[r0 + k]
+            k += 1
+        end
+        k
+    end
+
+    for orient in (false, true)
+        src = orient ? rev_run : run
+        pos_lookup = orient ? rev_pos : run_pos
+        s = 1
+        while s <= n
+            cv = contour[s]
+            candidates = get(pos_lookup, cv, nothing)
+            if candidates === nothing
+                s += 1
+                continue
+            end
+            best_k = 0
+            best_r0 = 0
+            for r0 in candidates
+                k = extend_match(s, r0, src)
+                if k > best_k
+                    best_k = k
+                    best_r0 = r0
+                end
+            end
+            if best_k >= min_len
+                run_lo = orient ? (m + 1 - (best_r0 + best_k - 1)) : best_r0
+                run_hi = orient ? (m + 1 - best_r0) : (best_r0 + best_k - 1)
+                push!(
+                    matches,
+                    (contour_start=s, run_lo=run_lo, run_hi=run_hi, reversed=orient)
+                )
+                s += best_k                # skip the matched block; no overlap
+            else
+                s += 1
+            end
+        end
+    end
+
+    # Drop entries that are subsumed by a longer entry covering the same contour range.
+    # Sort by length descending; keep an entry only if no kept entry already covers its
+    # contour_start..contour_start+(hi-lo).
+    isempty(matches) && return matches
+    by_len = sort(matches; by=t -> -(t.run_hi - t.run_lo + 1))
+    kept = empty(by_len)
+    for t in by_len
+        len = t.run_hi - t.run_lo + 1
+        overlap = false
+        for kt in kept
+            klen = kt.run_hi - kt.run_lo + 1
+            # Check cyclic overlap of [t.contour_start, t.contour_start+len-1]
+            # with [kt.contour_start, kt.contour_start+klen-1].
+            for i = 0:(len - 1)
+                pos = mod1(t.contour_start + i, n)
+                # Is pos inside kt's span?
+                for j = 0:(klen - 1)
+                    if mod1(kt.contour_start + j, n) == pos
+                        overlap = true
+                        break
+                    end
+                end
+                overlap && break
+            end
+            overlap && break
+        end
+        overlap || push!(kept, t)
+    end
+    return kept
+end
+
+# Given a curve and a sub-run [run_lo, run_hi] of its full discretization (length m),
+# return the corresponding sub-segment (the same kind of curve restricted to the
+# parametric range `t ∈ [(run_lo-1)/(m-1), (run_hi-1)/(m-1)] * pathlength(curve)`).
+# Reversed sub-runs return a reversed sub-segment.
+function _sub_curve(curve::Paths.Segment, m::Int, run_lo::Int, run_hi::Int, reversed::Bool)
+    pl = Paths.pathlength(curve)
+    t_lo = (run_lo - 1) / (m - 1) * pl
+    t_hi = (run_hi - 1) / (m - 1) * pl
+    # Take three-way split: pre, middle (the surviving sub-arc), post
+    if t_lo > zero(pl)
+        _, rest = Paths.split(curve, t_lo)
+    else
+        rest = curve
+    end
+    mid_len = t_hi - t_lo
+    if mid_len < Paths.pathlength(rest)
+        mid, _ = Paths.split(rest, mid_len)
+    else
+        mid = rest
+    end
+    return reversed ? Paths.reverse(mid) : mid
+end
+
+# Walk a ClippedPolygon's PolyNode tree, substituting known curves back into each contour
+# wherever a ProvenanceRun's discretized integer-grid point-run survived a boolean op intact.
+# Returns Vector{CurvilinearRegion{T}}: one region per outer contour (its direct children are
+# holes; grandchildren start new regions), mirroring to_primitives(::SolidModel, ::ClippedPolygon).
+#
+# Each run matches at most once globally — `used` is shared across all contours.
+# `report`, if given, collects (status::Symbol, curve, contour_index) tuples with
+# status ∈ (:recovered, :clipped); :clipped runs (never matched) carry contour_index 0.
+# When true, `substitute_curves` falls back to longest-contiguous-substring matching
+# for curves whose full run was cut by another polygon's edge — recovering each surviving
+# sub-block as a sub-segment via `Paths.split`, instead of dropping the whole curve to polyline.
+const substitute_curves_partial = Ref(false)
+
+function substitute_curves(clipped::ClippedPolygon{T}, runs; report=nothing) where {T}
+    out = CurvilinearRegion{T}[]
+    contour_index = Ref(0)
+    used = falses(length(runs))                  # shared across ALL contours
+    function build_cpoly(node)
+        pts = collect(node.contour)              # Vector{Point{T}}
+        snapped = clipperize.(pts)
+        n = length(pts)
+        contour_index[] += 1
+        ci = contour_index[]
+        # Find each surviving run's span: (start, m, segment), start 1-based into `pts`,
+        # m = number of contour vertices the run occupies (cyclically start … start+m-1).
+        matched = Tuple{Int, Int, Paths.Segment}[]
+        # match_run is exact and translation-sensitive: two geometrically distinct curves
+        # produce distinct integer runs, so first-match-wins cannot misattribute. Only
+        # overlapping (degenerate) curves could share a run, which is not handled.
+        for (ri, pr) in enumerate(runs)
+            used[ri] && continue
+            hit = match_run(snapped, pr.run)
+            if !isnothing(hit)
+                seg = hit.reversed ? reverse(pr.curve) : pr.curve
+                push!(matched, (hit.start, length(pr.run), seg))
+                used[ri] = true
+                !isnothing(report) && push!(report, (:recovered, pr.curve, ci))
+                continue
+            end
+            # Partial recovery: the full run didn't match (curve was cut by a boolean edge).
+            # Try longest-contiguous-substring matching: each surviving sub-block of the
+            # original run becomes a sub-segment of the original curve via Paths.split.
+            substitute_curves_partial[] || continue
+            partial_matches = match_run_partial(snapped, pr.run; min_len=3)
+            isempty(partial_matches) && continue
+            m_full = length(pr.run)
+            for pm in partial_matches
+                sub_len = pm.run_hi - pm.run_lo + 1
+                # sub_len includes both endpoints; the contour occupies sub_len positions.
+                sub_seg = _sub_curve(pr.curve, m_full, pm.run_lo, pm.run_hi, pm.reversed)
+                push!(matched, (pm.contour_start, sub_len, sub_seg))
+            end
+            used[ri] = true
+            !isnothing(report) && push!(report, (:recovered, pr.curve, ci))
+        end
+        isempty(matched) && return CurvilinearPolygon{T}(pts, Paths.Segment[], Int[])
+
+        # A curve replaces its discretized run: keep only the run's two endpoints, drop the
+        # m-2 interior vertices, and point curve_start_idx at the (reduced-list) start vertex.
+        # Mark interior positions to drop, and starts to tag.
+        drop = falses(n)
+        start_seg = Dict{Int, Paths.Segment}()
+        for (s, m, seg) in matched
+            start_seg[s] = seg
+            for k = 1:(m - 2)                    # strictly-interior positions, cyclic
+                drop[mod1(s + k, n)] = true
+            end
+        end
+        # Single ordered pass: build reduced point list and record curve starts at their
+        # index in the reduced list. (Starts/ends are never dropped, so they survive.)
+        reduced = Point{T}[]
+        curves = Paths.Segment[]
+        csis = Int[]
+        for i = 1:n
+            drop[i] && continue
+            push!(reduced, pts[i])
+            if haskey(start_seg, i)
+                push!(curves, start_seg[i])
+                push!(csis, length(reduced))      # this vertex's index in `reduced`
+            end
+        end
+        return CurvilinearPolygon{T}(reduced, curves, csis)
+    end
+    function add_region(node)
+        ext = build_cpoly(node)
+        # Holes must be CCW, but come out CW from Clipper
+        holes =
+            CurvilinearPolygon{T}[_reverse(build_cpoly(child)) for child in node.children]
+        push!(out, CurvilinearRegion{T}(ext, holes))
+        for child in node.children
+            for gc in child.children
+                add_region(gc)
+            end
+        end
+    end
+    for child in clipped.tree.children
+        add_region(child)
+    end
+    if !isnothing(report)
+        for (ri, pr) in enumerate(runs)
+            used[ri] || push!(report, (:clipped, pr.curve, 0))
+        end
+    end
+    return out
+end
+
+# Normalize a clip operand into a flat Vector of entities, preserving curve-bearing
+# entities intact (unlike _normalize_clip_arg, which discretizes via to_polygons).
+_as_entities(p::DeviceLayout.GeometryEntity) = [p]
+_as_entities(p::AbstractArray) = collect(Iterators.flatten(_as_entities.(p)))
+_as_entities(p::Union{GeometryStructure, GeometryReference}) =
+    _as_entities(flat_elements(p))
+_as_entities(p::Pair{<:Union{GeometryStructure, GeometryReference}}) =
+    _as_entities(flat_elements(p))
+function _as_entities(p::Paths.Node)
+    # Paths.-qualified: ContinuousStyle is defined in the Paths submodule but is NOT in scope
+    # unqualified inside Curvilinear (a bare `ContinuousStyle` throws UndefVarError on the L1
+    # zero-length-node path). Qualifying is robust and needs no import.
+    iszero(pathlength(p.seg)) && p.sty isa Paths.ContinuousStyle && return []
+    return _as_entities(pathtopolys(p.seg, p.sty)) # Use `islinear` dispatch on segment and style
+end
+# A Rounded-styled straight Polygon recovers as its exact-arc CurvilinearPolygon, so
+# corners survive the clip. Mirrors to_primitives(::SolidModel, ::StyledEntity{Polygon,Rounded}).
+function _as_entities(
+    p::StyledEntity{T, <:Union{Polygon{T}, Polygons.Rectangle{T}}, <:Polygons.Rounded}
+) where {T}
+    return [
+        round_to_curvilinearpolygon(
+            p.ent,
+            radius(p.sty),
+            min_side_len=p.sty.min_side_len,
+            corner_indices=cornerindices(p.ent, p.sty),
+            min_angle=p.sty.min_angle
+        )
+    ]
+end
+# A Rounded- or StyleDict{Rounded}-styled ClippedPolygon (e.g. the output of a non-curved
+# `difference2d`/`union2d` followed by post-clip rounding) recovers as exact fillet arcs
+# per contour, so its corners survive the clip instead of discretizing. The render
+# path already does this conversion (`to_curvilinear_regions` walks the clipped tree, applying
+# `styled_loop`+`round_to_curvilinearpolygon` per contour → CurvilinearRegions with arcs); we
+# reuse that exact function so boolean recovery matches the render. `SolidModels` loads after
+# this file, but `_as_entities` only runs at boolean time (all modules loaded), so the qualified
+# name resolves at call time. Returns a Vector{CurvilinearRegion}; `_as_entities(::AbstractArray)`
+# flattens it, and `_collect_provenance!(::CurvilinearRegion)` records each contour's curve runs.
+function _as_entities(p::StyledEntity{T, ClippedPolygon{T}, <:Polygons.StyleDict}) where {T}
+    return DeviceLayout.SolidModels.to_curvilinear_regions(p.ent, p.sty)
+end
+function _as_entities(p::StyledEntity{T, ClippedPolygon{T}, <:Polygons.Rounded}) where {T}
+    return DeviceLayout.SolidModels.to_curvilinear_regions(p.ent, Polygons.StyleDict(p.sty))
+end
+# A geometrically-transparent style (e.g. MeshSized — a mesh-density hint applied via
+# `to_polygons(ent, ::MeshSized) = to_polygons(ent)`) must NOT block curve recovery: a
+# MeshSized-wrapped Node / Rounded-Polygon / CurvilinearPolygon carries the same curves as
+# the bare entity. Without this, such a wrapped entity falls to the generic
+# `_as_entities(::GeometryEntity) = [p]` and is discretized by `to_polygons`, silently losing
+# its arcs through the boolean. Recurse to the inner entity (mirrors the render-side default
+# `to_primitives(::SolidModel, ::StyledEntity) = to_primitives(sm, ent.ent)`). The specific
+# Rounded / ClippedPolygon methods above are more specific and still win.
+_as_entities(p::StyledEntity) = _as_entities(p.ent)
+
+# Promoted coordinate type matching what clip's promote_type would pick.
+function _recover_coordtype(plus, minus)
+    return promote_type(
+        DeviceLayout.coordinatetype(plus),
+        DeviceLayout.coordinatetype(minus)
+    )
+end
+
+"""
+    recover_curves(op, plus, minus; report=nothing)
+
+Run a boolean operation `op` on `plus` and `minus`, then recover original curves from the
+discretized result wherever they survived the operation intact. Returns
+`Vector{CurvilinearRegion}` — one region per outer contour in the clipped result (each
+disjoint piece becomes a separate region).
+
+## Positional arguments
+
+The first argument `op` must be one of the polygon clipping operations: `difference2d`,
+`union2d`, `intersect2d`, or `xor2d`.
+
+The second and third arguments (`plus` and `minus`) accept the same forms as the
+corresponding clipping operation: a `GeometryEntity` or array of `GeometryEntity`, a
+`GeometryStructure` or `GeometryReference` (whose flattened elements are used), or a pair
+`geom => layer` selecting only elements in those layers from the flattened structure.
+
+Curve-bearing entities have their curves tracked through discretization and recovered in
+the result: `CurvilinearPolygon`, `CurvilinearRegion`, `Path` nodes (e.g. `Turn`/`BSpline`
+segments rendered with a `Style`), and `Rounded`-styled `Polygon`/`Rectangle`. All other
+entities are discretized via `to_polygons`.
+
+## Keyword arguments
+
+`report` may be a `Vector` that will be filled with `(status, curve, contour_index)` tuples
+tracking recovery results. `status` is `:recovered` if the curve's entire discretized run
+survived the boolean operation intact, or `:clipped` if it was cut and fell back to a
+polyline. `contour_index` is the 1-based index of the output contour where the curve was
+recovered, or `0` for `:clipped` curves.
+
+## Return type
+
+Returns `Vector{CurvilinearRegion}` (one region per outer contour). This differs from
+`difference2d` and similar functions, which return a single `ClippedPolygon`. Callers
+migrating from `difference2d(a, b)` to `recover_curves(difference2d, a, b)` must handle a
+vector result.
+
+## Limitations
+
+**All-or-nothing recovery:** A curve is recovered only if its entire discretized run (the
+sequence of integer-grid vertices produced by `discretize_curve`) survives the boolean
+operation with exact integer equality. If the operation cuts through a curve, that curve is
+reported `:clipped` and falls back to a polyline. Partial-curve recovery is not supported.
+
+Additionally, curves can only be recovered on CurvilinearRegion/CurvilinearPolygon,
+Path nodes, and `Rounded`-styled `Polygon`/`Rectangle` entities. `Rounded` applied to
+other entities, styled Curvilinear entities, and nested styles do not yet support curve recovery.
+
+See also [`difference2d`](@ref), [`union2d`](@ref), [`intersect2d`](@ref), [`xor2d`](@ref),
+[`difference2d_curved`](@ref), [`union2d_curved`](@ref), [`intersect2d_curved`](@ref),
+[`xor2d_curved`](@ref).
+"""
+function recover_curves(op, plus, minus; report=nothing)
+    R = _recover_coordtype(plus, minus)
+    pp, runs_p = discretize_with_provenance(_as_entities(plus), R)
+    pm, runs_m = discretize_with_provenance(_as_entities(minus), R)
+    # Annotate the result so a wrong `op` (not one of difference2d/union2d/intersect2d/
+    # xor2d) fails here with a clear TypeError rather than deep inside substitute_curves.
+    clipped = op(pp, pm)::ClippedPolygon
+    return substitute_curves(clipped, vcat(runs_p, runs_m); report=report)
+end
+
+"""
+    difference2d_curved(plus, minus; report=nothing)
+
+Curve-preserving variant of [`difference2d`](@ref), returning `Vector{CurvilinearRegion}`.
+See [`recover_curves`](@ref).
+"""
+difference2d_curved(p, m; kwargs...) = recover_curves(difference2d, p, m; kwargs...)
+"""
+    union2d_curved(p1, p2; report=nothing)
+    union2d_curved(p; report=nothing)
+
+Curve-preserving variant of [`union2d`](@ref), returning `Vector{CurvilinearRegion}`.
+The single-argument form self-unions `p` (equivalent to `union2d_curved(p, [])`), which is
+useful for merging a collection of overlapping curved entities into one region per piece.
+See [`recover_curves`](@ref).
+"""
+union2d_curved(p, m; kwargs...) = recover_curves(union2d, p, m; kwargs...)
+union2d_curved(p; kwargs...) = recover_curves(union2d, p, []; kwargs...)
+
+"""
+    intersect2d_curved(p1, p2; report=nothing)
+
+Curve-preserving variant of [`intersect2d`](@ref), returning `Vector{CurvilinearRegion}`.
+See [`recover_curves`](@ref).
+"""
+intersect2d_curved(p, m; kwargs...) = recover_curves(intersect2d, p, m; kwargs...)
+"""
+    xor2d_curved(p1, p2; report=nothing)
+
+Curve-preserving variant of [`xor2d`](@ref), returning `Vector{CurvilinearRegion}`.
+See [`recover_curves`](@ref).
+"""
+xor2d_curved(p, m; kwargs...) = recover_curves(xor2d, p, m; kwargs...)
diff --git a/src/curvilinear.jl b/src/curvilinear.jl
index bc072cf74..b7b646a21 100644
--- a/src/curvilinear.jl
+++ b/src/curvilinear.jl
@@ -11,14 +11,27 @@ import DeviceLayout:
     AbstractPolygon,
     GeometryEntity,
     GeometryEntityStyle,
+    GeometryStructure,
+    GeometryReference,
     Paths,
     Reflection,
     Rotation,
     ScaledIsometry,
+    StyledEntity,
     Transformation,
     Translation
 import DeviceLayout:
-    to_polygons, points, rotation, origin, mag, xrefl, transform, perimeter, isapprox_angle
+    flat_elements,
+    to_polygons,
+    points,
+    rotation,
+    round_to_curvilinearpolygon,
+    origin,
+    mag,
+    xrefl,
+    transform,
+    perimeter,
+    isapprox_angle
 using DeviceLayout.Paths
 import DeviceLayout.Polygons: cornerindices
 
@@ -27,6 +40,7 @@ using ..Polygons
 using ..Paths
 
 export CurvilinearPolygon, CurvilinearRegion, pathtopolys, line_arc_cornerindices
+export recover_curves, difference2d_curved, intersect2d_curved, union2d_curved, xor2d_curved
 
 """
     struct CurvilinearPolygon{T} <: GeometryEntity{T}
@@ -85,6 +99,16 @@ struct CurvilinearPolygon{T} <: GeometryEntity{T}
         for (curve_idx, start_idx) in enumerate(csi)
             csi[curve_idx] = start_idx - count(dup_idx .< start_idx)
         end
+        # Consumers (`to_polygons`, `_collect_provenance!`) walk curves with a running
+        # cursor and slice `p[i:csi]`, which requires `csi` ascending. Producers such as
+        # `_reverse` and reflective `transform` can emit unsorted indices when a curve
+        # wraps around to the last vertex, so sort curves and indices jointly here to
+        # restore the invariant for every producer.
+        if length(csi) > 1
+            perm = sortperm(csi)
+            c = c[perm]
+            csi = csi[perm]
+        end
         return new{T}(p, c, csi)
     end
 end
@@ -94,7 +118,7 @@ function CurvilinearPolygon(points::Vector{Point{T}}) where {T}
     # Straight segments are implicit
     return CurvilinearPolygon{T}(points, Paths.Segment[], Int[])
 end
-CurvilinearPolygon(p::Polygon{T}) where {T} = CurvilinearPolygon{T}(points(p))
+CurvilinearPolygon(p::Polygon{T}) where {T} = CurvilinearPolygon(points(p))
 
 ### Conversion methods
 function to_polygons(
@@ -106,7 +130,20 @@ function to_polygons(
     i = 1
     p = Point{T}[]
 
-    for (idx, (csi, c)) ∈ enumerate(zip(e.curve_start_idx, e.curves))
+    # Walk curves in ascending start-index order. The running index `i` and `e.p[i:csi]` below
+    # assume `curve_start_idx` is monotonically increasing; a curve supplied out of order — in
+    # particular one at the wrap seam (csi == length(e.p)) listed first — would make `e.p[i:csi]`
+    # dump the whole point ring up front and append that curve's arc points at the tail, ending
+    # one step short of its endpoint and leaving a sub-µm near-pinch. Sorting locally (without
+    # mutating `e`) makes the walk robust to any curve ordering, including deserialized polygons
+    # that predate constructor-side normalization.
+    order =
+        issorted(e.curve_start_idx) ? eachindex(e.curve_start_idx) :
+        sortperm(e.curve_start_idx)
+
+    for idx in order
+        csi = e.curve_start_idx[idx]
+        c = e.curves[idx]
         # Add the points from current to start of curve
         append!(p, e.p[i:csi])
 
@@ -132,6 +169,15 @@ function to_polygons(
     return Polygon{T}(p)
 end
 
+function _reverse(e::CurvilinearPolygon)
+    csi_rev = (i, N) -> mod1(i + 1, N) - N - 1
+    return CurvilinearPolygon(
+        reverse(e.p),
+        reverse(e.curves),
+        reverse(csi_rev.(e.curve_start_idx, length(e.p)))
+    )
+end
+
 function transform(e::CurvilinearPolygon, f::Transformation)
     # If the transformation is a reflection, have to fix the winding.
     # curve_start_idx are shifted forward 1, reversed, then negated.
@@ -960,4 +1006,6 @@ function rounded_corner_line_arc(
     return (; points=fillet_pts, T_arc=T_arc)
 end
 
+include("curve_recovery.jl")
+
 end # module
diff --git a/src/render/render.jl b/src/render/render.jl
index 61758aeba..f5f2144fa 100644
--- a/src/render/render.jl
+++ b/src/render/render.jl
@@ -213,3 +213,91 @@ Additional keyword arguments are passed to [`to_polygons`](@ref) for each entity
 certain entity types to control how they are converted to polygons.
 """
 render!(c::Cell, s::GeometryStructure; kwargs...) = _render!(c, s; kwargs...)
+
+####### Rendering pathway through Curvilinear
+function round_to_curvilinearpolygon(
+    pol::GeometryEntity{T},
+    radius::S;
+    corner_indices=eachindex(points(pol)),
+    line_arc_corner_indices=nothing,
+    min_angle=1e-3,
+    relative::Bool=(T <: Length) && (S <: Real),
+    min_side_len=relative ? zero(T) : radius
+) where {T, S <: Coordinate}
+    # If radius is dimensional, non-relative rounding.
+    V = float(T)
+    # Tie break for Real, Real introduces a type instability for non-dimensional.
+    relative = ((T <: Length) && (S <: Real)) || (relative && T <: Real && S <: Real)
+
+    poly = points(pol)
+    len = length(poly)
+    new_points = Point{V}[]
+    new_curves = Paths.Turn{V}[]
+    new_curve_start_idx = Int[]
+
+    for i in eachindex(poly)
+        if !(i in corner_indices)
+            push!(new_points, poly[i])
+        else
+            p0 = poly[mod1(i - 1, len)] # handles the cyclic boundary condition
+            p1 = poly[i]
+            p2 = poly[mod1(i + 1, len)]
+            radius_dim = relative ? radius * min(norm(p0 - p1), norm(p1 - p2)) : radius
+            seg_or_p1 = rounded_corner_segment(
+                p0,
+                p1,
+                p2,
+                radius_dim,
+                min_side_len=min_side_len,
+                min_angle=min_angle
+            )
+            if seg_or_p1 isa Paths.Turn
+                push!(new_points, Paths.p0(seg_or_p1))
+                push!(new_curves, seg_or_p1)
+                push!(new_curve_start_idx, length(new_points))
+                push!(new_points, Paths.p1(seg_or_p1))
+            else
+                push!(new_points, seg_or_p1)
+            end
+        end
+    end
+
+    return CurvilinearPolygon(new_points, new_curves, new_curve_start_idx)
+end
+
+function rounded_corner_segment(
+    p0::Point{T},
+    p1::Point{T},
+    p2::Point{T},
+    radius::S;
+    min_side_len=radius,
+    min_angle=1e-3
+) where {T, S <: Coordinate}
+    V = float(T)
+    rad = convert(V, radius)
+
+    v1 = (p1 - p0) / norm(p1 - p0)
+    v2 = (p2 - p1) / norm(p2 - p1)
+    α1 = atan(v1.y, v1.x) # between -π and π
+    α2 = atan(v2.y, v2.x)
+
+    if min_side_len > norm(p1 - p0) || min_side_len > norm(p2 - p1) # checks that the side lengths against min_side_len
+        return p1
+    elseif isapprox(rem2pi(α1 - α2, RoundNearest), 0, atol=min_angle) # checks if the points are collinear, within tolerance
+        return p1
+    end
+
+    dir = orientation(p0, p1, p2) # checks the direction of the corner
+    dα = α2 - α1 # always between +/- 2π
+    if sign(dα) != dir # Make sure turn is in the correct direction
+        dα = dα + dir * 2π # Still between +/- 2π
+    end
+
+    # p0_seg is the start of the arc, determined by the intersection
+    # of lines parallel to v1, v2
+    k =
+        inv([v1.x -v2.x; v1.y -v2.y]) *
+        [p2.x - p0.x + dir * rad * (v1.y - v2.y), p2.y - p0.y + dir * rad * (v2.x - v1.x)]
+    p0_seg = p0 + k[1] * v1
+    return Paths.Turn(uconvert(°, dα), rad, p0_seg, uconvert(°, α1))
+end
diff --git a/src/solidmodels/render.jl b/src/solidmodels/render.jl
index 16987f05b..84cf6902f 100644
--- a/src/solidmodels/render.jl
+++ b/src/solidmodels/render.jl
@@ -188,56 +188,6 @@ function to_primitives(
     )
 end
 
-function round_to_curvilinearpolygon(
-    pol::GeometryEntity{T},
-    radius::S;
-    corner_indices=eachindex(points(pol)),
-    line_arc_corner_indices=nothing,
-    min_angle=1e-3,
-    relative::Bool=(T <: Length) && (S <: Real),
-    min_side_len=relative ? zero(T) : radius
-)::CurvilinearPolygon{T} where {T, S <: Coordinate}
-    # If radius is dimensional, non-relative rounding.
-    V = float(T)
-    # Tie break for Real, Real introduces a type instability for non-dimensional.
-    relative = ((T <: Length) && (S <: Real)) || (relative && T <: Real && S <: Real)
-
-    poly = points(pol)
-    len = length(poly)
-    new_points = Point{V}[]
-    new_curves = Paths.Turn{V}[]
-    new_curve_start_idx = Int[]
-
-    for i in eachindex(poly)
-        if !(i in corner_indices)
-            push!(new_points, poly[i])
-        else
-            p0 = poly[mod1(i - 1, len)] # handles the cyclic boundary condition
-            p1 = poly[i]
-            p2 = poly[mod1(i + 1, len)]
-            radius_dim = relative ? radius * min(norm(p0 - p1), norm(p1 - p2)) : radius
-            seg_or_p1 = rounded_corner_segment(
-                p0,
-                p1,
-                p2,
-                radius_dim,
-                min_side_len=min_side_len,
-                min_angle=min_angle
-            )
-            if seg_or_p1 isa Paths.Turn
-                push!(new_points, Paths.p0(seg_or_p1))
-                push!(new_curves, seg_or_p1)
-                push!(new_curve_start_idx, length(new_points))
-                push!(new_points, Paths.p1(seg_or_p1))
-            else
-                push!(new_points, seg_or_p1)
-            end
-        end
-    end
-
-    return CurvilinearPolygon(new_points, new_curves, new_curve_start_idx)
-end
-
 function round_to_curvilinearpolygon(
     pol::CurvilinearPolygon{T},
     radius::S;
@@ -363,54 +313,10 @@ function round_to_curvilinearpolygon(
             push!(new_curve_start_idx, csi)
         end
     end
-
-    # Sort curves by start index so to_polygons can iterate in vertex order
-    if length(new_curve_start_idx) > 1
-        perm = sortperm(new_curve_start_idx)
-        new_curves = new_curves[perm]
-        new_curve_start_idx = new_curve_start_idx[perm]
-    end
-
+    # Constructor will sort curves by start index
     return CurvilinearPolygon(new_points, new_curves, new_curve_start_idx)
 end
 
-function rounded_corner_segment(
-    p0::Point{T},
-    p1::Point{T},
-    p2::Point{T},
-    radius::S;
-    min_side_len=radius,
-    min_angle=1e-3
-) where {T, S <: Coordinate}
-    V = float(T)
-    rad = convert(V, radius)
-
-    v1 = (p1 - p0) / norm(p1 - p0)
-    v2 = (p2 - p1) / norm(p2 - p1)
-    α1 = atan(v1.y, v1.x) # between -π and π
-    α2 = atan(v2.y, v2.x)
-
-    if min_side_len > norm(p1 - p0) || min_side_len > norm(p2 - p1) # checks that the side lengths against min_side_len
-        return p1
-    elseif isapprox(rem2pi(α1 - α2, RoundNearest), 0, atol=min_angle) # checks if the points are collinear, within tolerance
-        return p1
-    end
-
-    dir = orientation(p0, p1, p2) # checks the direction of the corner
-    dα = α2 - α1 # always between +/- 2π
-    if sign(dα) != dir # Make sure turn is in the correct direction
-        dα = dα + dir * 2π # Still between +/- 2π
-    end
-
-    # p0_seg is the start of the arc, determined by the intersection
-    # of lines parallel to v1, v2
-    k =
-        inv([v1.x -v2.x; v1.y -v2.y]) *
-        [p2.x - p0.x + dir * rad * (v1.y - v2.y), p2.y - p0.y + dir * rad * (v2.x - v1.x)]
-    p0_seg = p0 + k[1] * v1
-    return Paths.Turn(uconvert(°, dα), rad, p0_seg, uconvert(°, α1))
-end
-
 # TODO: The fillet center geometry here duplicates `rounded_corner_line_arc` in
 # curvilinear.jl. Extract the shared solver into a common helper when consolidating
 # rounding logic through CurvilinearPolygon.
diff --git a/src/solidmodels/solidmodels.jl b/src/solidmodels/solidmodels.jl
index 38cc920a2..1acc806a0 100644
--- a/src/solidmodels/solidmodels.jl
+++ b/src/solidmodels/solidmodels.jl
@@ -66,6 +66,8 @@ import DeviceLayout:
     level,
     norm,
     render!,
+    round_to_curvilinearpolygon,
+    rounded_corner_segment,
     uniquename,
     isapprox_angle
 import DeviceLayout.Paths: trace, gap, offset, extent, pathlength, bspline_approximation
diff --git a/test/test_line_arc_rounding.jl b/test/test_line_arc_rounding.jl
index 468b5a1bb..5edcbfa25 100644
--- a/test/test_line_arc_rounding.jl
+++ b/test/test_line_arc_rounding.jl
@@ -739,3 +739,63 @@ end
     place!(cs_sm, tiny_poly, GDSMeta(0, 0))
     @test_nowarn render!(sm, cs_sm)
 end
+
+@testitem "CurvilinearPolygon walk is curve-order independent (wrap-seam regression)" setup =
+    [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, Rounded, to_polygons, points
+    import DeviceLayout.SolidModels: styled_loop
+    # Regression for the seam-ordering bug: to_polygons(::CurvilinearPolygon) and
+    # _collect_provenance! advance a running index `i` and emit `p[i:csi]` before each curve,
+    # which assumes `curve_start_idx` is ascending. A curve at the wrap seam (csi == length(p))
+    # listed out of order makes the walk dump the whole ring up front and append that curve's arc
+    # at the tail, ending one discretization step (~tens of nm) short of its endpoint → a sub-µm
+    # near-pinch that a downstream boolean closes into a zero-area sliver (a degenerate mesh
+    # element). The walk must be independent of the order curves are supplied in.
+
+    # styled_loop on a box produces a valid multi-arc CurvilinearPolygon (one fillet per corner),
+    # including a corner at the wrap seam — exactly the shape that triggered the bug.
+    box = CurvilinearPolygon(Point{Float64}[(0, 0), (10, 0), (10, 10), (0, 10)] .* μm)
+    cp = styled_loop(box, Rounded(1.0μm))
+    @test cp isa CurvilinearPolygon
+    @test !isempty(cp.curves)
+
+    # Antenna signature: two vertices far apart on the ring (gap > 5) but spatially < 0.5µm,
+    # both sitting on long (> 0.5µm) edges — i.e. two boundary strands that should have met.
+    function has_pinch(poly)
+        q = [
+            (
+                Float64(ustrip(μm, DeviceLayout.getx(z))),
+                Float64(ustrip(μm, DeviceLayout.gety(z)))
+            ) for z in points(poly)
+        ]
+        m = length(q)
+        for i = 1:m, j = (i + 5):m
+            (i <= 2 && j >= m - 1) && continue
+            d = hypot(q[i][1] - q[j][1], q[i][2] - q[j][2])
+            (0.001 < d < 0.5) || continue
+            ei = max(
+                hypot(q[mod1(i + 1, m)][1] - q[i][1], q[mod1(i + 1, m)][2] - q[i][2]),
+                hypot(q[i][1] - q[mod1(i - 1, m)][1], q[i][2] - q[mod1(i - 1, m)][2])
+            )
+            ej = max(
+                hypot(q[mod1(j + 1, m)][1] - q[j][1], q[mod1(j + 1, m)][2] - q[j][2]),
+                hypot(q[j][1] - q[mod1(j - 1, m)][1], q[j][2] - q[mod1(j - 1, m)][2])
+            )
+            (ei > 0.5 && ej > 0.5) && return true
+        end
+        return false
+    end
+
+    # Natural order (as styled_loop produced it) — must be clean.
+    poly_natural = to_polygons(cp)
+    @test !has_pinch(poly_natural)
+
+    # Reverse the curve order so a high-index (seam-adjacent) curve is processed first. This is
+    # the condition a deserialized polygon hits when its curves were stored out of canonical order. The constructor normalizes by sorting,
+    # and the walk sorts locally, so the result must match the natural order and stay pinch-free.
+    perm = reverse(eachindex(cp.curve_start_idx))
+    cp_perm = CurvilinearPolygon(copy(cp.p), cp.curves[perm], cp.curve_start_idx[perm])
+    poly_perm = to_polygons(cp_perm)
+    @test !has_pinch(poly_perm)
+    @test length(points(poly_perm)) == length(points(poly_natural))
+end
diff --git a/test/test_provenance_recovery.jl b/test/test_provenance_recovery.jl
new file mode 100644
index 000000000..1f128df03
--- /dev/null
+++ b/test/test_provenance_recovery.jl
@@ -0,0 +1,570 @@
+@testitem "Provenance — preserved vertices survive Clipper" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, Polygon, union2d
+    snap(p) = DeviceLayout.Polygons.clipperize(p)
+
+    R = 10.0
+    N = 64
+    circ = [Point(R * cos(t), R * sin(t)) for t in range(0, 2π, length=N + 1)[1:(end - 1)]]
+    poly = Polygon(circ)
+    in_int = snap.(circ)
+
+    sq = Polygon([
+        Point(1000.0, 1000.0),
+        Point(1001.0, 1000.0),
+        Point(1001.0, 1001.0),
+        Point(1000.0, 1001.0)
+    ])
+    res = union2d(poly, sq)
+
+    function allcontours(node, acc)
+        for c in node.children
+            push!(acc, c.contour)
+            allcontours(c, acc)
+        end
+        return acc
+    end
+    cons = allcontours(res.tree, Vector{Point{Float64}}[])
+    circ_con = argmax(length, cons)
+    out_set = Set(snap.(circ_con))
+
+    @test all(p -> p in out_set, in_int)
+    @test length(circ_con) == N
+end
+
+@testitem "Provenance — discretize_with_provenance captures arc runs" setup =
+    [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, Paths, Polygon
+    using DeviceLayout.Curvilinear: discretize_with_provenance
+    # Create a square with one curved edge: (0,0) -> (10,0) -> Turn -> (0,10)
+    # Turn from (10,0) with α0=90° (pointing up), radius 10, sweeping 90° ends at (0,10)
+    pts = Point{Float64}[(0, 0), (10, 0), (0, 10)]
+    turn = Paths.Turn(90.0°, 10.0; p0=Point(10.0, 0.0), α0=90.0°)
+    cpoly = CurvilinearPolygon(pts, [turn], [2])
+    polys, runs = discretize_with_provenance([cpoly], Float64)
+    @test polys isa Vector{<:Polygon}
+    @test length(runs) == 1
+    @test runs[1].curve.p0 === turn.p0 # Can't just check segment identity bc Julia 1.10 makes a copy in Curvilinear constructor
+    @test runs[1].curve.α0 === turn.α0
+    @test runs[1].curve.α === turn.α
+    @test runs[1].curve.r === turn.r
+    @test runs[1].run[1] == DeviceLayout.Polygons.clipperize(Point(10.0, 0.0))
+    @test length(runs[1].run) ≥ 3
+end
+
+@testitem "Provenance — match_run cyclic + bidirectional" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point
+    M = DeviceLayout.Curvilinear   # match_run is internal (not exported)
+    ip(x) = Point{Int64}(x, 0)
+    contour = ip.([1, 2, 3, 4, 5, 6])
+    @test M.match_run(contour, ip.([2, 3, 4])) == (start=2, reversed=false)
+    @test M.match_run(contour, ip.([5, 6, 1])) == (start=5, reversed=false)  # cyclic wrap
+    @test M.match_run(contour, ip.([4, 3, 2])) == (start=2, reversed=true)   # reversed
+    @test M.match_run(contour, ip.([2, 4, 3])) === nothing                   # no contiguous match
+end
+
+@testitem "Provenance — substitute_curves recovers one arc" setup = [CommonTestSetup] begin
+    using DeviceLayout:
+        Point, CurvilinearPolygon, CurvilinearRegion, Paths, union2d, Polygon
+    using DeviceLayout.Curvilinear: discretize_with_provenance, substitute_curves
+    # A CurvilinearPolygon whose exterior has a single recoverable arc.
+    pts = Point{Float64}[(0, 0), (10, 0), (0, 10)]
+    turn = Paths.Turn(90.0°, 10.0; p0=Point(10.0, 0.0), α0=90.0°)
+    region = CurvilinearRegion(CurvilinearPolygon(pts, [turn], [2]))
+
+    polys, runs = discretize_with_provenance([region], Float64)
+    sq = Polygon([
+        Point(1e4, 1e4),
+        Point(1e4 + 1, 1e4),
+        Point(1e4 + 1, 1e4 + 1),
+        Point(1e4, 1e4 + 1)
+    ])
+    clipped = union2d(polys, [sq])              # ClippedPolygon, arc untouched
+    report = Tuple[]
+    out = substitute_curves(clipped, runs; report=report)
+    @test out isa Vector{<:CurvilinearRegion}
+    @test sum(length(r.exterior.curves) for r in out) == 1
+    @test count(t -> t[1] == :recovered, report) == 1
+    @test count(t -> t[1] == :clipped, report) == 0
+end
+
+@testitem "Provenance — substitute_curves recovers two arcs in order" setup =
+    [CommonTestSetup] begin
+    using DeviceLayout:
+        Point, CurvilinearPolygon, CurvilinearRegion, Paths, union2d, Polygon
+    using DeviceLayout.Curvilinear: discretize_with_provenance, substitute_curves
+    # A CurvilinearPolygon with two recoverable arcs (two rounded corners) on its exterior.
+    # Exercises the multi-curve ordering path: matched runs must be sorted by start index
+    # so to_polygons' monotonic cursor doesn't drop/misorder vertices.
+    t1 = Paths.Turn(90.0°, 2.0; p0=Point(10.0, 0.0), α0=0.0°)   # (10,0) heading +x → (12,2)
+    t2 = Paths.Turn(90.0°, 2.0; p0=Point(12.0, 10.0), α0=90.0°) # (12,10) heading +y → (10,12)
+    pts = Point{Float64}[(0, 0), (10, 0), (12, 2), (12, 10), (10, 12), (0, 12)]
+    region = CurvilinearRegion(CurvilinearPolygon(pts, [t1, t2], [2, 4]))
+
+    polys, runs = discretize_with_provenance([region], Float64)
+    sq = Polygon([
+        Point(1e4, 1e4),
+        Point(1e4 + 1, 1e4),
+        Point(1e4 + 1, 1e4 + 1),
+        Point(1e4, 1e4 + 1)
+    ])
+    clipped = union2d(polys, [sq])
+    report = Tuple[]
+    out = substitute_curves(clipped, runs; report=report)
+    # Clipper may order the disjoint square's region first, so select the curved
+    # region by curve count rather than assuming an index.
+    curved = out[findfirst(r -> !isempty(r.exterior.curves), out)]
+    @test length(curved.exterior.curves) == 2
+    # curve_start_idx must be ascending after sorting.
+    @test issorted(curved.exterior.curve_start_idx)
+    # to_polygons must not error or drop everything when walking the two curves.
+    @test length(DeviceLayout.to_polygons(curved.exterior).p) > 0
+    @test count(t -> t[1] == :recovered, report) == 2
+    @test count(t -> t[1] == :clipped, report) == 0
+end
+
+@testitem "Provenance — recover_curves end to end" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, union2d, union2d_curved, difference2d_curved
+    pts = Point{Float64}[(0, 0), (10, 0), (0, 10)]
+    turn = Paths.Turn(90.0°, 10.0; p0=Point(10.0, 0.0), α0=90.0°)
+    region = CurvilinearRegion(CurvilinearPolygon(pts, [turn], [2]))
+    sq = Polygon([
+        Point(1e4, 1e4),
+        Point(1e4 + 1, 1e4),
+        Point(1e4 + 1, 1e4 + 1),
+        Point(1e4, 1e4 + 1)
+    ])
+
+    a = recover_curves(union2d, region, sq)
+    b = union2d_curved(region, sq)
+    @test a isa Vector{<:CurvilinearRegion}
+    @test b isa Vector{<:CurvilinearRegion}
+    @test sum(length(r.exterior.curves) for r in a) == 1
+    @test sum(length(r.exterior.curves) for r in b) == 1
+    # report kwarg flows through
+    report = Tuple[]
+    recover_curves(union2d, region, sq; report=report)
+    @test count(t -> t[1] == :recovered, report) == 1
+end
+
+@testitem "Provenance — interface methods" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, union2d, union2d_curved, difference2d_curved
+    pts = Point{Float64}[(0, 0), (10, 0), (0, 10)]
+    turn = Paths.Turn(90.0°, 10.0; p0=Point(10.0, 0.0), α0=90.0°)
+    region = CurvilinearRegion(CurvilinearPolygon(pts, [turn], [2]))
+    sq = Polygon([
+        Point(1e4, 1e4),
+        Point(1e4 + 1, 1e4),
+        Point(1e4 + 1, 1e4 + 1),
+        Point(1e4, 1e4 + 1)
+    ])
+    cs1 = CoordinateSystem{Float64}("test1")
+    cs2 = CoordinateSystem{Float64}("test2")
+    place!(cs1, region, :curved)
+    place!(cs1, sq, :square)
+    place!(cs2, sq, :square)
+    out = union2d_curved(region, sq)
+    a = union2d_curved(cs1)
+    b = union2d_curved(region, cs2)
+    c = union2d_curved(cs1 => :curved, cs1 => :square)
+    d = union2d_curved([sq, cs1 => :curved])
+    @test isempty(to_polygons(xor2d(out, a)))
+    @test isempty(to_polygons(xor2d(a, b)))
+    @test isempty(to_polygons(xor2d(b, c)))
+    @test isempty(to_polygons(xor2d(c, d)))
+    # Mutli-polygon ClippedPolygon input
+    multi_cp = union2d(sq, sq + Point(10, 10))
+    @test isempty(to_polygons(xor2d(multi_cp, union2d_curved(multi_cp))))
+end
+
+@testitem "Provenance — Path round trips" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, union2d, union2d_curved, difference2d_curved
+    # Trace
+    pa = Path()
+    turn!(pa, 90°, 50μm, Paths.Trace(5μm))
+    turn!(pa, 90°, 50μm, Paths.Trace(5μm))
+    out = union2d_curved(pa)
+    @test length(out) == 1 # Unioned into a single CurvilinearRegion
+    curved = out[1]
+    @test length(curved.exterior.curves) == 4 # original 4 curves
+    @test length(curved.exterior.p) == 6 # original 6 points
+    @test isempty(to_polygons(xor2d(curved, pathtopolys(pa))))
+
+    # CPW
+    pa = Path()
+    turn!(pa, 90°, 50μm, Paths.CPW(5μm, 5μm))
+    turn!(pa, 90°, 50μm, Paths.CPW(5μm, 5μm))
+    out = union2d_curved(pa)
+    @test length(out) == 2
+    curved = out[1]
+    @test length(curved.exterior.curves) == 4 # original 4 curves
+    @test length(curved.exterior.p) == 6 # original 6 points
+    @test isempty(to_polygons(xor2d(out, pathtopolys(pa))))
+
+    # Generic curves
+    pa = Path()
+    turn!(pa, 90°, 50μm, Paths.TaperTrace(5μm, 10μm))
+    out = union2d_curved(pa)
+    @test length(out) == 1 # Unioned into a single CurvilinearRegion
+    curved = out[1]
+    @test length(curved.exterior.curves) == 2 # Recognizes offset turn
+    @test isempty(to_polygons(xor2d(curved, pathtopolys(pa))))
+    # BSpline offset
+    pa = Path()
+    bspline!(pa, [Point(0μm, 1mm)], 180°, Paths.Trace(10μm))
+    out = union2d_curved(pa)
+    @test length(out) == 1 # Unioned into a single CurvilinearRegion
+    curved = out[1]
+    @test length(curved.exterior.curves) == 2
+    @test isempty(to_polygons(xor2d(curved, pathtopolys(pa))))
+end
+
+@testitem "Provenance — Rounded polygon round trips" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, union2d, union2d_curved, difference2d_curved
+    r = centered(Rectangle(10.0μm, 10.0μm))
+    rr = Rounded(r, 1μm)
+    out = intersect2d_curved(r, rr)
+    curved = out[1]
+    @test length(curved.exterior.curves) == 4
+    @test length(curved.exterior.p) == 8
+    @test isempty(to_polygons(xor2d(curved, Curvilinear._as_entities(rr))))
+    r2 = centered(Rectangle(100μm, 2μm))
+    out = union2d_curved(rr, r2)
+    curved = out[1]
+    @test length(curved.exterior.curves) == 4
+    @test length(curved.exterior.p) == 16
+    @test isempty(to_polygons(xor2d(curved, [r2, Curvilinear._as_entities(rr)])))
+    out = difference2d_curved(rr, r2)
+    @test length(out) == 2
+    @test length(out[1].exterior.curves) == 2
+    @test length(out[2].exterior.curves) == 2
+    @test isempty(to_polygons(xor2d(out, difference2d(Curvilinear._as_entities(rr), r2))))
+end
+
+@testitem "Provenance — single untouched arc, varied radius" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, union2d, to_polygons
+    # An arc-bearing region unioned with a DISJOINT square recovers its arc intact
+    # across a range of radii (the integer-grid run survives byte-identical).
+    for R in (1.0, 100.0, 1.0μm, 100.0μm, 1.0mm)
+        z = zero(R)
+        T = typeof(R)
+        um = DeviceLayout.onemicron(T)
+        turn = Paths.Turn(90.0°, R; p0=Point(R, z), α0=90.0°)  # (R,0) heading +y → (0,R)
+        region = CurvilinearRegion(
+            CurvilinearPolygon(Point{T}[(z, z), (R, z), (z, R)], [turn], [2])
+        )
+        sq = Polygon([
+            Point(1e4, 1e4)um,
+            Point(1e4 + 1, 1e4)um,
+            Point(1e4 + 1, 1e4 + 1)um,
+            Point(1e4, 1e4 + 1)um
+        ])
+        report = Tuple[]
+        out = recover_curves(union2d, region, sq; report=report)
+        @test count(t -> t[1] == :recovered, report) == 1
+        @test count(t -> t[1] == :clipped, report) == 0
+        curved = out[findfirst(r -> !isempty(r.exterior.curves), out)]
+        @test length(curved.exterior.curves) == 1
+        @test length(curved.exterior.p) == 3
+        @test length(to_polygons(curved.exterior).p) > 0
+        @test isempty(to_polygons(xor2d(region, curved)))
+    end
+end
+
+@testitem "Provenance — seam-rotation (union with self)" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, union2d
+    # A rectangle with one rounded corner unioned with a copy of itself. Clipper
+    # rotates the start vertex of the output contour so the seam lands on a kink.
+    # The cyclic match_run handles it.
+    t = Paths.Turn(90.0°, 2.0; p0=Point(18.0, 0.0), α0=0.0°)  # (18,0) heading +x → (20,2)
+    region = CurvilinearRegion(
+        CurvilinearPolygon(
+            Point{Float64}[(0, 0), (18, 0), (20, 2), (20, 10), (0, 10)],
+            [t],
+            [2]
+        )
+    )
+    report = Tuple[]
+    out = recover_curves(union2d, region, region; report=report)
+    # The arc survives the union-with-self and is recovered onto the merged contour.
+    @test count(t -> t[1] == :recovered, report) ≥ 1
+    @test count(t -> t[1] == :clipped, report) == 0
+    curved = out[findfirst(r -> !isempty(r.exterior.curves), out)]
+    @test length(curved.exterior.curves) == 1
+    @test length(curved.exterior.p) == 5
+    @test curved.exterior.curves[1] isa Paths.Turn
+    @test isempty(to_polygons(xor2d(region, out)))
+end
+
+@testitem "Provenance — reversed winding (arc on a hole)" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout.Curvilinear: recover_curves, difference2d, discretize_with_provenance
+    # Differencing a curved region out of a solid square puts the arc on the
+    # resulting hole boundary, whose winding is reversed relative to the input
+    # exterior. The reversed branch of match_run must fire to recover it.
+    plus = Polygon([Point(0.0, 0), Point(30, 0), Point(30, 30), Point(0, 30)])
+    tm = Paths.Turn(90.0°, 2.0; p0=Point(18.0, 10.0), α0=0.0°)  # (18,10) → (20,12)
+    minus = CurvilinearRegion(
+        CurvilinearPolygon(
+            Point{Float64}[(10, 10), (18, 10), (20, 12), (20, 20), (10, 20)],
+            [tm],
+            [2]
+        )
+    )
+    report = Tuple[]
+    out = recover_curves(difference2d, plus, minus; report=report)
+    @test count(t -> t[1] == :recovered, report) == 1
+    @test count(t -> t[1] == :clipped, report) == 0
+    # The recovered curve lands on a hole, not the exterior.
+    @test sum(length(r.exterior.curves) for r in out) == 0
+    @test sum(sum(length(h.curves) for h in r.holes; init=0) for r in out) == 1
+
+    # Confirm the reversed branch of match_run actually fired on the hole contour.
+    R = Float64
+    polys_plus, _ = discretize_with_provenance([plus], R)
+    polys_minus, runs_m = discretize_with_provenance([minus], R)
+    clipped = difference2d(polys_plus, polys_minus)
+    # Output is geometrically identical
+    @test isempty(to_polygons(xor2d(clipped, out)))
+    saw_reversed = Ref(false)
+    walk(node) =
+        for c in node.children
+            snapped = DeviceLayout.Polygons.clipperize.(collect(c.contour))
+            for pr in runs_m
+                hit = Curvilinear.match_run(snapped, pr.run)
+                (hit !== nothing && hit.reversed) && (saw_reversed[] = true)
+            end
+            walk(c)
+        end
+    walk(clipped.tree)
+    @test saw_reversed[]
+
+    # Now curve recovery with CR with a hole as input
+    self_union = union2d_curved(out)
+    @test sum(length(r.exterior.curves) for r in self_union) == 0
+    @test sum(sum(length(h.curves) for h in r.holes; init=0) for r in self_union) == 1
+    @test isempty(to_polygons(xor2d(self_union, out)))
+    @test isempty(to_polygons(xor2d(self_union, union2d(out))))
+
+    # Clipper op returns ClippedPolygon with grandchildren
+    outer = centered(Rectangle(100.0, 100.0))
+    curved_gc = difference2d_curved(outer, out)
+    @test length(curved_gc) == 2 # Inner is separate region
+    @test sum(length(r.exterior.curves) for r in curved_gc) == 1
+    @test sum(sum(length(h.curves) for h in r.holes; init=0) for r in curved_gc) == 0
+    @test isempty(to_polygons(xor2d(curved_gc, difference2d(outer, out))))
+    @test isempty(to_polygons(xor2d(curved_gc, difference2d(outer, out))))
+end
+
+@testitem "Provenance — arcs on exterior and a hole" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, difference2d
+    # A region with an arc on its exterior, differenced by a curved region that
+    # becomes a hole carrying its own arc. Both arcs must be recovered, one on
+    # the exterior contour and one on the hole contour.
+    te = Paths.Turn(90.0°, 2.0; p0=Point(28.0, 0.0), α0=0.0°)  # exterior arc (28,0) → (30,2)
+    plus = CurvilinearRegion(
+        CurvilinearPolygon(
+            Point{Float64}[(0, 0), (28, 0), (30, 2), (30, 30), (0, 30)],
+            [te],
+            [2]
+        )
+    )
+    th = Paths.Turn(90.0°, 2.0; p0=Point(18.0, 10.0), α0=0.0°)  # hole arc (18,10) → (20,12)
+    minus = CurvilinearRegion(
+        CurvilinearPolygon(
+            Point{Float64}[(10, 10), (18, 10), (20, 12), (20, 20), (10, 20)],
+            [th],
+            [2]
+        )
+    )
+    report = Tuple[]
+    out = recover_curves(difference2d, plus, minus; report=report)
+    @test count(t -> t[1] == :recovered, report) == 2
+    @test count(t -> t[1] == :clipped, report) == 0
+    n_ext = sum(length(r.exterior.curves) for r in out)
+    n_hole = sum(sum(length(h.curves) for h in r.holes; init=0) for r in out)
+    @test n_ext == 1   # one arc on an exterior contour
+    @test n_hole == 1  # one arc on a hole contour
+    @test isempty(to_polygons(xor2d(difference2d(plus, minus), out)))
+end
+
+@testitem "Provenance — annulus / collision probe" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout.Curvilinear: recover_curves, union2d, discretize_with_provenance
+    # Two arcs of different radii in one geometry. Each must map to its own curve.
+    # Different radii produce different integer runs, so no cross-match collision
+    # is possible — we assert the runs differ to make that intent explicit.
+    t1 = Paths.Turn(90.0°, 2.0; p0=Point(10.0, 0.0), α0=0.0°)   # r=2: (10,0) → (12,2)
+    t2 = Paths.Turn(90.0°, 5.0; p0=Point(12.0, 10.0), α0=90.0°) # r=5: (12,10) → (7,15)
+    pts = Point{Float64}[(0, 0), (10, 0), (12, 2), (12, 10), (7, 15), (0, 15)]
+    region = CurvilinearRegion(CurvilinearPolygon(pts, [t1, t2], [2, 4]))
+    sq = Polygon([
+        Point(1e4, 1e4),
+        Point(1e4 + 1, 1e4),
+        Point(1e4 + 1, 1e4 + 1),
+        Point(1e4, 1e4 + 1)
+    ])
+    report = Tuple[]
+    out = recover_curves(union2d, region, sq; report=report)
+    @test count(t -> t[1] == :recovered, report) == 2
+    @test count(t -> t[1] == :clipped, report) == 0
+    curved = out[findfirst(r -> !isempty(r.exterior.curves), out)]
+    @test length(curved.exterior.curves) == 2
+    radii = sort([c.r for c in curved.exterior.curves])
+    @test radii == [2.0, 5.0]   # each arc recovered as its own distinct curve
+    @test isempty(to_polygons(xor2d(region, curved)))
+
+    # The two provenance runs differ → no cross-match collision is possible.
+    _, runs = discretize_with_provenance([region], Float64)
+    @test runs[1].run != runs[2].run
+end
+
+@testitem "Provenance — straight edge cuts an arc (clipped, not recovered)" setup =
+    [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, difference2d
+    # A straight rectangle that crosses the arc interior. Clipper inserts a new
+    # intersection vertex on the arc, breaking its run. Recovery is all-or-nothing,
+    # so the cut arc is reported :clipped (not recovered) and falls back to polyline.
+    t = Paths.Turn(90.0°, 2.0; p0=Point(18.0, 10.0), α0=0.0°)  # (18,10) → (20,12)
+    region = CurvilinearRegion(
+        CurvilinearPolygon(
+            Point{Float64}[(10, 10), (18, 10), (20, 12), (20, 20), (10, 20)],
+            [t],
+            [2]
+        )
+    )
+    cutter = Polygon([
+        Point(18.5, 10.5),
+        Point(25.0, 10.5),
+        Point(25.0, 13.0),
+        Point(18.5, 13.0)
+    ])
+    report = Tuple[]
+    out = recover_curves(difference2d, region, cutter; report=report)
+    @test count(t -> t[1] == :recovered, report) == 0
+    @test count(t -> t[1] == :clipped, report) == 1
+    @test sum(length(r.exterior.curves) for r in out) == 0  # polyline fallback
+end
+
+@testitem "Provenance — full circle across the seam" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, union2d, to_polygons
+    # A closed-loop circle built from two half-circle Turns (a single 360° Turn is
+    # degenerate — start and end vertices coincide), unioned with a disjoint square.
+    # Both half-arcs must be recovered onto the single closed contour.
+    R = 10.0
+    t1 = Paths.Turn(180.0°, R; p0=Point(R, 0.0), α0=90.0°)    # top half: (R,0) → (-R,0)
+    t2 = Paths.Turn(180.0°, R; p0=Point(-R, 0.0), α0=-90.0°)  # bottom half: (-R,0) → (R,0)
+    region = CurvilinearRegion(
+        CurvilinearPolygon(Point{Float64}[(R, 0), (-R, 0)], [t1, t2], [1, 2])
+    )
+    sq = Polygon([
+        Point(1e4, 1e4),
+        Point(1e4 + 1, 1e4),
+        Point(1e4 + 1, 1e4 + 1),
+        Point(1e4, 1e4 + 1)
+    ])
+    report = Tuple[]
+    out = recover_curves(union2d, region, sq; report=report)
+    @test count(t -> t[1] == :recovered, report) == 2
+    @test count(t -> t[1] == :clipped, report) == 0
+    curved = out[findfirst(r -> !isempty(r.exterior.curves), out)]
+    @test length(curved.exterior.p) == 2
+    @test length(curved.exterior.curves) == 2
+    @test length(to_polygons(curved.exterior).p) > 0
+    @test isempty(to_polygons(xor2d(region, curved)))
+end
+
+@testitem "Provenance — BSpline segment (type-agnostic)" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, union2d
+    # A BSpline-bearing region unioned with a disjoint square. The recovery
+    # mechanism is curve-type-agnostic: it must recover the BSpline, not just arcs.
+    bpts = Point{Float64}[(10, 0), (13, 3), (16, 2), (18, 6)]
+    bs = Paths.BSpline(bpts, Point(1.0, 1.0), Point(1.0, 1.0))
+    pts = Point{Float64}[(0, 0), bs.p0, bs.p1, (0, 10)]
+    region = CurvilinearRegion(CurvilinearPolygon(pts, [bs], [2]))
+    sq = Polygon([
+        Point(1e4, 1e4),
+        Point(1e4 + 1, 1e4),
+        Point(1e4 + 1, 1e4 + 1),
+        Point(1e4, 1e4 + 1)
+    ])
+    report = Tuple[]
+    out = recover_curves(union2d, region, sq; report=report)
+    @test count(t -> t[1] == :recovered, report) == 1
+    @test count(t -> t[1] == :clipped, report) == 0
+    curved = out[findfirst(r -> !isempty(r.exterior.curves), out)]
+    @test length(curved.exterior.curves) == 1
+    @test curved.exterior.curves[1] isa Paths.BSpline
+    @test isempty(to_polygons(xor2d(region, curved)))
+end
+
+@testitem "Provenance — clipped curve not spuriously recovered" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, CurvilinearRegion, Paths, Polygon
+    using DeviceLayout: recover_curves, difference2d
+    # Safety property: a destroyed arc must NEVER emit a (wrong) curve. Uses the same
+    # cut-arc geometry as the "straight edge cuts an arc" case and asserts zero
+    # recovery, a :clipped report, and a polyline (zero-curve) exterior.
+    t = Paths.Turn(90.0°, 2.0; p0=Point(18.0, 10.0), α0=0.0°)
+    region = CurvilinearRegion(
+        CurvilinearPolygon(
+            Point{Float64}[(10, 10), (18, 10), (20, 12), (20, 20), (10, 20)],
+            [t],
+            [2]
+        )
+    )
+    cutter = Polygon([
+        Point(18.5, 10.5),
+        Point(25.0, 10.5),
+        Point(25.0, 13.0),
+        Point(18.5, 13.0)
+    ])
+    report = Tuple[]
+    out = recover_curves(difference2d, region, cutter; report=report)
+    @test count(t -> t[1] == :recovered, report) == 0
+    @test any(t -> t[1] == :clipped, report)
+    @test all(isempty(r.exterior.curves) for r in out)
+end
+
+@testitem "Provenance — curve_start_idx sorted invariant" setup = [CommonTestSetup] begin
+    using DeviceLayout: Point, CurvilinearPolygon, Paths, to_polygons, Reflection, °
+    using DeviceLayout.Curvilinear: _reverse
+    # to_polygons walks curves with a running cursor and slices p[i:csi], which requires
+    # curve_start_idx ascending. _reverse and reflective transform both flip a curve that
+    # wraps to the last vertex into the unsorted order [4, 2], which either throws an
+    # AssertionError or silently duplicates vertices (corrupt geometry). The constructor
+    # sorts curves and indices jointly so every producer preserves the invariant.
+
+    # Stadium: straight bottom/top, semicircle caps. The left cap (vertex 4 → vertex 1)
+    # is the wrap-around curve that triggers the bug under reversal/reflection.
+    pts = Point{Float64}[(0, 0), (10, 0), (10, 10), (0, 10)]
+    cap_r = Paths.Turn(180.0°, 5.0; p0=Point(10.0, 0.0), α0=0.0°)    # (10,0) → (10,10)
+    cap_l = Paths.Turn(180.0°, 5.0; p0=Point(0.0, 10.0), α0=180.0°)  # (0,10) → (0,0)
+    cp = CurvilinearPolygon(copy(pts), [cap_r, cap_l], [2, 4])
+    n_fwd = length(to_polygons(cp).p)
+
+    # Constructor reorders curves jointly with indices when given unsorted input.
+    unsorted = CurvilinearPolygon(copy(pts), [cap_l, cap_r], [4, 2])
+    @test issorted(unsorted.curve_start_idx)
+    @test unsorted.curves == cp.curves
+    @test to_polygons(unsorted).p == to_polygons(cp).p
+
+    # _reverse produces a wrap-around curve (csi would be [4, 2]) — must stay sorted and
+    # round-trip without the extra vertices the corrupt path introduced.
+    rev = _reverse(cp)
+    @test issorted(rev.curve_start_idx)
+    @test length(to_polygons(rev).p) == n_fwd
+
+    # Reflective transform uses the same csi-flip logic and had the identical bug.
+    refl = DeviceLayout.transform(cp, Reflection(0.0°))
+    @test issorted(refl.curve_start_idx)
+    @test length(to_polygons(refl).p) == n_fwd
+end