diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4cbe1ba1..fef98e63 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,12 @@ The format of this changelog is based on
 [Keep a Changelog](https://keepachangelog.com/), and this project adheres to
 [Semantic Versioning](https://semver.org/).
 
+## Upcoming
+
+  - Added layerwise Booleans `union2d_layerwise`, `difference2d_layerwise`, `intersect2d_layerwise`, and `xor2d_layerwise`
+  - Added `Polygons.area`
+  - Fixed overly-strict argument types for polygon clipping methods
+
 ## 1.12.0 (2026-04-13)
 
   - Added `auto_union` SolidModel rendering option; if `true`, self-unions every 2D group before any other postrendering (default `false`)
diff --git a/docs/make.jl b/docs/make.jl
index 443593ed..1975611a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -30,7 +30,7 @@ makedocs(
     format=Documenter.HTML(
         prettyurls=true,
         assets=["assets/favicon.ico"],
-        size_threshold=300_000,
+        size_threshold=400_000,
         collapselevel=1
     ),
     sitename="DeviceLayout.jl",
diff --git a/docs/src/reference/api.md b/docs/src/reference/api.md
index 8d74e06f..598eaf2f 100644
--- a/docs/src/reference/api.md
+++ b/docs/src/reference/api.md
@@ -158,6 +158,7 @@
     Polygon(::AbstractVector{Point{T}}) where {T}
     Polygon(::Point, ::Point, ::Point, ::Point...)
     Rectangle
+    Polygons.area
     bounds
     circle_polygon
     gridpoints_in_polygon
@@ -174,10 +175,15 @@
 ```@docs
     Polygons.ClippedPolygon
     difference2d
+    difference2d_layerwise
     intersect2d
+    intersect2d_layerwise
     union2d
+    union2d_layerwise
     xor2d
+    xor2d_layerwise
     clip
+    clip_tiled
     Polygons.StyleDict
 ```
 
diff --git a/src/DeviceLayout.jl b/src/DeviceLayout.jl
index 560eff75..b693eccd 100644
--- a/src/DeviceLayout.jl
+++ b/src/DeviceLayout.jl
@@ -19,6 +19,8 @@ import Unitful: Length, LengthUnits, DimensionlessQuantity, NoUnits, DimensionEr
 import Unitful: ustrip, unit, inch
 Unitful.@derived_dimension InverseLength inv(Unitful.𝐋)
 
+import SpatialIndexing
+
 function render! end
 export render!
 
@@ -195,6 +197,7 @@ coordinatetype(::AbstractArray{S}) where {T, S <: AbstractGeometry{T}} = T
 coordinatetype(iterable) = promote_type(coordinatetype.(iterable)...)
 coordinatetype(::Point{T}) where {T} = T
 coordinatetype(::Type{Point{T}}) where {T} = T
+coordinatetype(::Pair{<:AbstractGeometry{T}}) where {T} = T
 
 # Entity interface
 include("entities.jl")
@@ -378,10 +381,13 @@ import .Polygons:
     circle,
     circle_polygon,
     clip,
+    clip_tiled,
     cliptree,
     difference2d,
+    difference2d_layerwise,
     gridpoints_in_polygon,
     intersect2d,
+    intersect2d_layerwise,
     offset,
     perimeter,
     points,
@@ -390,7 +396,9 @@ import .Polygons:
     sweep_poly,
     unfold,
     union2d,
-    xor2d
+    union2d_layerwise,
+    xor2d,
+    xor2d_layerwise
 export Polygons,
     Polygon,
     Ellipse,
@@ -406,8 +414,10 @@ export Polygons,
     clip,
     cliptree,
     difference2d,
+    difference2d_layerwise,
     gridpoints_in_polygon,
     intersect2d,
+    intersect2d_layerwise,
     offset,
     perimeter,
     points,
@@ -416,7 +426,9 @@ export Polygons,
     sweep_poly,
     unfold,
     union2d,
-    xor2d
+    union2d_layerwise,
+    xor2d,
+    xor2d_layerwise
 
 include("align.jl")
 using .Align
diff --git a/src/curvilinear.jl b/src/curvilinear.jl
index ebb97aa0..f8c9d0ba 100644
--- a/src/curvilinear.jl
+++ b/src/curvilinear.jl
@@ -52,9 +52,16 @@ struct CurvilinearPolygon{T} <: GeometryEntity{T}
     p::Vector{Point{T}}
     curves::Vector{<:Paths.Segment} # Only need to store non-line-segment curves
     curve_start_idx::Vector{Int} # And the indices at which they start
-    # A negative start idx like -3 means that the corresponding curve
-    # between p[3] and p[4] is actually parameterized from p[4] to p[3]
+    # Backward-parameterized curves (negative start idx) are normalized in the constructor:
+    # the segment is reversed and the index flipped positive.
     function CurvilinearPolygon{T}(p, c, csi) where {T} # Make sure you don't have zero-length curves
+        # Normalize backward-parameterized curves: reverse segment, flip index positive.
+        for i in eachindex(csi)
+            if csi[i] < 0
+                c[i] = reverse(c[i])
+                csi[i] = -csi[i]
+            end
+        end
         # Don't treat duplicates in any different fashion -> view as user error
         # Some endpoint pairs may be identical; delete the duplicates
         # Maybe inefficient but least confusing to iterate to find them and then delete
@@ -100,24 +107,24 @@ function to_polygons(
 
     for (idx, (csi, c)) ∈ enumerate(zip(e.curve_start_idx, e.curves))
         # Add the points from current to start of curve
-        append!(p, e.p[i:abs(csi)])
+        append!(p, e.p[i:csi])
 
         # Discretize segment using tolerance-based adaptive grid.
-        wrapped_i = mod1(abs(csi) + 1, length(e.p))
+        wrapped_i = mod1(csi + 1, length(e.p))
         pp = DeviceLayout.discretize_curve(c, atol)
 
         # Remove the calculated points corresponding to start and end.
-        term_p = csi < 0 ? popfirst!(pp) : pop!(pp)
-        init_p = csi < 0 ? pop!(pp) : popfirst!(pp)
+        term_p = pop!(pp)
+        init_p = popfirst!(pp)
 
         # Add interior points and bump counter.
-        append!(p, csi < 0 ? reverse(pp) : pp)
-        i = abs(csi) + 1
+        append!(p, pp)
+        i = csi + 1
 
         # Ensure that the calculated start and end points match the non-calculated points.
         @assert !isapprox(init_p, term_p; atol=1e-3 * DeviceLayout.onenanometer(T)) "Curve $idx must have non-zero length!"
-        @assert isapprox(term_p, e.p[wrapped_i]; atol=1e-3 * DeviceLayout.onenanometer(T)) "Curve $idx must $(csi < 0 ? "start" : "end") at point $(wrapped_i)!"
-        @assert isapprox(init_p, e.p[abs(csi)]; atol=1e-3 * DeviceLayout.onenanometer(T)) "Curve $idx must $(csi < 0 ? "end" : "start") at point $(abs(csi))!"
+        @assert isapprox(term_p, e.p[wrapped_i]; atol=1e-3 * DeviceLayout.onenanometer(T)) "Curve $idx must end at point $(wrapped_i)!"
+        @assert isapprox(init_p, e.p[csi]; atol=1e-3 * DeviceLayout.onenanometer(T)) "Curve $idx must start at point $(csi)!"
     end
     append!(p, e.p[i:end])
 
@@ -463,8 +470,7 @@ end
 # Only indices that don't start or end a curve are available for rounding.
 # cornerindices(p::CurvilinearPolygon, s::GeometryEntityStyle) = cornerindices(p, p0(s))
 function cornerindices(p::CurvilinearPolygon{T}) where {T}
-    curve_bound_ind =
-        vcat((x -> [abs(x), (abs(x) % length(p.p)) + 1]).(p.curve_start_idx)...)
+    curve_bound_ind = vcat((x -> [x, (x % length(p.p)) + 1]).(p.curve_start_idx)...)
     valid_ind = setdiff(1:length(p.p), curve_bound_ind)
     return valid_ind
 end
@@ -542,9 +548,9 @@ Returns a NamedTuple `(incoming=..., outgoing=...)` where each field is either
 `:straight` or the `Paths.Segment` (e.g., `Paths.Turn`) for that edge.
 
   - **Outgoing edge** (from `p[i]` to `p[i+1]`): curved if any
-    `abs(curve_start_idx[k]) == i`
+    `curve_start_idx[k] == i`
   - **Incoming edge** (from `p[i-1]` to `p[i]`): curved if any
-    `abs(curve_start_idx[k]) == mod1(i-1, n)`
+    `curve_start_idx[k] == mod1(i-1, n)`
 """
 function edge_type_at_vertex(p::CurvilinearPolygon, i::Int)
     n = length(p.p)
@@ -552,16 +558,12 @@ function edge_type_at_vertex(p::CurvilinearPolygon, i::Int)
 
     incoming = :straight
     outgoing = :straight
-    # k = index into curves array, csi = vertex index where curve k starts.
-    # Negative csi means the curve is parameterized in reverse; use reverse(curve)
-    # so p0/p1 match the vertex order expected by downstream rounding code.
     for (k, csi) in enumerate(p.curve_start_idx)
-        curve = csi < 0 ? reverse(p.curves[k]) : p.curves[k]
-        if abs(csi) == prev_i
-            incoming = curve
+        if csi == prev_i
+            incoming = p.curves[k]
         end
-        if abs(csi) == i
-            outgoing = curve
+        if csi == i
+            outgoing = p.curves[k]
         end
     end
     return (; incoming=incoming, outgoing=outgoing)
@@ -596,7 +598,7 @@ function to_polygons(
     la_corners = Set(line_arc_cornerindices(ent, sty))
 
     # Precompute vertex for curve index lookup
-    vertex_to_curve = Dict(abs(csi) => k for (k, csi) in enumerate(ent.curve_start_idx))
+    vertex_to_curve = Dict(csi => k for (k, csi) in enumerate(ent.curve_start_idx))
 
     # Round corners, tracking which original vertex maps to which output range.
     # Also track T_arc for each line-arc corner so we can trim the curves.
@@ -696,19 +698,13 @@ function to_polygons(
             csi = ent.curve_start_idx[curve_k]
             arc_len = pathlength(c)
 
-            # Determine trim parameters: find t values for trim points on the arc.
-            # For negative csi, the curve is parameterized in reverse, so trim_start
-            # (at the forward-start vertex) maps to a high t value and trim_end to a
-            # low t value. Swap them so t_start < t_end for correct discretization.
             t_start = zero(S)
             t_end = arc_len
-            ts_dict = csi < 0 ? trim_end : trim_start
-            te_dict = csi < 0 ? trim_start : trim_end
-            if haskey(ts_dict, curve_k)
-                t_start = Paths.pathlength_nearest(c, ts_dict[curve_k])
+            if haskey(trim_start, curve_k)
+                t_start = Paths.pathlength_nearest(c, trim_start[curve_k])
             end
-            if haskey(te_dict, curve_k)
-                t_end = Paths.pathlength_nearest(c, te_dict[curve_k])
+            if haskey(trim_end, curve_k)
+                t_end = Paths.pathlength_nearest(c, trim_end[curve_k])
             end
 
             # Discretize the trimmed portion of the curve.
@@ -725,7 +721,7 @@ function to_polygons(
                 )
                 # Remove endpoints (already present as fillet tangent points)
                 inner = grid[(begin + 1):(end - 1)] .* l
-                pp = c.(csi < 0 ? reverse(inner) : inner)
+                pp = c.(inner)
                 append!(final_points, pp)
             end
         end
diff --git a/src/entities.jl b/src/entities.jl
index f252ab39..12432978 100644
--- a/src/entities.jl
+++ b/src/entities.jl
@@ -403,3 +403,53 @@ lowerleft(aent::ArrayEntity) = lowerleft(aent.a)
 upperright(aent::ArrayEntity) = upperright(aent.a)
 halo(aent::ArrayEntity, outer_delta, inner_delta=nothing) =
     ArrayEntity(halo(aent.a, outer_delta, inner_delta))
+
+######## Entity selection
+function SpatialIndexing.mbr(ent::AbstractGeometry{T}) where {T}
+    r = bounds(ent)
+    return SpatialIndexing.Rect(
+        (ustrip(unit(onemicron(T)), float.(r.ll))...,),
+        (ustrip(unit(onemicron(T)), float.(r.ur))...,)
+    )
+end
+
+"""
+    mbr_spatial_index(geoms)
+
+An `RTree` of minimum bounding rectangles for `geoms`, with indices in `geoms` as values.
+
+See also [`findbox`](@ref).
+"""
+function mbr_spatial_index(geoms)
+    tree = SpatialIndexing.RTree{Float64, 2}(Int)
+    function convertel(enum_ent)
+        idx, ent = enum_ent
+        return SpatialIndexing.SpatialElem(SpatialIndexing.mbr(ent), nothing, idx)
+    end
+    SpatialIndexing.load!(tree, enumerate(geoms), convertel=convertel)
+    return tree
+end
+
+"""
+    findbox(box, geoms; intersects=false)
+    findbox(box, tree::SpatialIndexing.RTree; intersects=false)
+
+Return `indices` such that `geoms[indices]` gives all elements of `geoms` with minimum bounding rectangle contained in `bounds(box)`.
+
+A spatial index created with [`mbr_spatial_index(geoms)`](@ref) can be supplied explicitly to avoid re-indexing for multiple `findbox` operations.
+
+By default, `findbox` will find only entities with bounds contained in `bounds(box)`. If `intersects` is `true`,
+it also includes entities with bounds intersecting `bounds(box)` (including those only touching edge-to-edge).
+"""
+function findbox(box, geoms; intersects=false)
+    tree = mbr_spatial_index(geoms)
+    return findbox(box, tree; intersects)
+end
+
+function findbox(box, tree::SpatialIndexing.RTree; intersects=false)
+    intersects && return map(
+        x -> x.val,
+        SpatialIndexing.intersects_with(tree, SpatialIndexing.mbr(box))
+    )
+    return map(x -> x.val, SpatialIndexing.contained_in(tree, SpatialIndexing.mbr(box)))
+end
diff --git a/src/polygons.jl b/src/polygons.jl
index 5e3b487f..d2ff5667 100644
--- a/src/polygons.jl
+++ b/src/polygons.jl
@@ -11,6 +11,8 @@ import CoordinateTransformations: AffineMap, LinearMap, Translation, Transformat
 import Clipper
 import Clipper: children, contour
 import StaticArrays
+import SpatialIndexing
+import SpatialIndexing: mbr
 
 import DeviceLayout
 import DeviceLayout:
@@ -83,6 +85,9 @@ export circle,
     union2d,
     xor2d
 
+export difference2d_layerwise,
+    intersect2d_layerwise, union2d_layerwise, xor2d_layerwise, clip_tiled
+
 const USCALE = 1.0 * Unitful.fm
 const SCALE  = 10.0^9
 
@@ -536,6 +541,20 @@ function perimeter(e::Ellipse)
     return π * ((a + b) + 3 * (a - b)^2 / (10 * (a + b) + sqrt(a^2 + 14 * a * b + b^2)))
 end
 
+function signed_area(p::Polygon) # Can be negative
+    return sum(
+        (gety.(p.p) + gety.(circshift(p.p, -1))) .*
+        (getx.(p.p) - getx.(circshift(p.p, -1)))
+    ) / 2
+end
+
+"""
+    area(poly::Polygon)
+
+Area of a non-self-intersecting polygon. Always positive.
+"""
+area(p::Polygon) = abs(signed_area(p))
+
 """
     circle_polygon(r, Δθ=10°)
 
@@ -942,8 +961,8 @@ end
         kwargs...) where {S, T, A<:Polygon{S}, B<:Polygon{T}}
     clip(op::Clipper.ClipType,
         s::AbstractVector{Polygon{T}}, c::AbstractVector{Polygon{T}};
-        pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd,
-        pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd) where {T}
+        pfs::Clipper.PolyFillType=Clipper.PolyFillTypePositive,
+        pfc::Clipper.PolyFillType=Clipper.PolyFillTypePositive) where {T}
 
 Return the `ClippedPolygon` resulting from a polygon clipping operation.
 
@@ -982,7 +1001,7 @@ These arguments may include:
 See the [`Clipper` docs](http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Types/PolyFillType.htm)
 for further information.
 
-See also [union2d](@ref), [difference2d](@ref), and [intersect2d](@ref).
+See also [union2d](@ref), [difference2d](@ref), [intersect2d](@ref), and [xor2d](@ref).
 """
 function clip(op::Clipper.ClipType, s, c; kwargs...)
     return clip(op, _normalize_clip_arg(s), _normalize_clip_arg(c); kwargs...)
@@ -1000,6 +1019,8 @@ _normalize_clip_arg(p::Union{GeometryStructure, GeometryReference}) =
     _normalize_clip_arg(flat_elements(p))
 _normalize_clip_arg(p::Pair{<:Union{GeometryStructure, GeometryReference}}) =
     _normalize_clip_arg(flat_elements(p))
+_normalize_clip_arg(p::AbstractArray) =
+    reduce(vcat, _normalize_clip_arg.(p); init=Polygon{DeviceLayout.coordinatetype(p)}[])
 
 # Clipping arrays of AbstractPolygons
 function clip(
@@ -1024,8 +1045,8 @@ function clip(
     op::Clipper.ClipType,
     s::AbstractVector{Polygon{T}},
     c::AbstractVector{Polygon{T}};
-    pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd,
-    pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd
+    pfs::Clipper.PolyFillType=Clipper.PolyFillTypePositive,
+    pfc::Clipper.PolyFillType=Clipper.PolyFillTypePositive
 ) where {T}
     sc, cc = clipperize(s), clipperize(c)
     polys = _clip(op, sc, cc; pfs, pfc)
@@ -1039,8 +1060,8 @@ end
         kwargs...) where {S, T, A<:AbstractPolygon{S}, B<:AbstractPolygon{T}}
     cliptree(op::Clipper.ClipType,
         s::AbstractVector{Polygon{T}}, c::AbstractVector{Polygon{T}};
-        pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd,
-        pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd) where {T}
+        pfs::Clipper.PolyFillType=Clipper.PolyFillTypePositive,
+        pfc::Clipper.PolyFillType=Clipper.PolyFillTypePositive) where {T}
 
 Return a `Clipper.PolyNode` representing parent-child relationships between polygons and
 interior holes. The units and number type may need to be converted.
@@ -1104,8 +1125,8 @@ function cliptree(
     op::Clipper.ClipType,
     s::AbstractVector{Polygon{T}},
     c::AbstractVector{Polygon{T}};
-    pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd,
-    pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd
+    pfs::Clipper.PolyFillType=Clipper.PolyFillTypePositive,
+    pfc::Clipper.PolyFillType=Clipper.PolyFillTypePositive
 ) where {T}
     sc, cc = clipperize(s), clipperize(c)
     cpoly = _clip(op, sc, cc; pfs, pfc)
@@ -1150,7 +1171,8 @@ end
 Return the geometric union of `p` or all entities in `p`.
 """
 union2d(p::AbstractGeometry{T}) where {T} = union2d(p, Polygon{T}[])
-union2d(p::AbstractArray{<:AbstractGeometry{T}}) where {T} = union2d(p, Polygon{T}[])
+union2d(p::AbstractArray) = union2d(p, Polygon{DeviceLayout.coordinatetype(p)}[])
+union2d(p::Pair{<:AbstractGeometry{T}}) where {T} = union2d(p, Polygon{T}[])
 
 """
     difference2d(p1, p2)
@@ -1254,8 +1276,8 @@ function _clip(
     op::Clipper.ClipType,
     s::AbstractVector{Polygon{T}},
     c::AbstractVector{Polygon{T}};
-    pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd,
-    pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd
+    pfs::Clipper.PolyFillType=Clipper.PolyFillTypePositive,
+    pfc::Clipper.PolyFillType=Clipper.PolyFillTypePositive
 ) where {T <: Union{Int64, Unitful.Quantity{Int64}}}
     clip = clipper()
     Clipper.clear!(clip)
@@ -2293,4 +2315,242 @@ function transform(d::StyleDict, f::Transformation)
     return newdict
 end
 
+### Layerwise
+function xor2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...)
+    return clip_layerwise(Clipper.ClipTypeXor, obj, tool; kwargs...)
+end
+
+function difference2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...)
+    return clip_layerwise(Clipper.ClipTypeDifference, obj, tool; kwargs...)
+end
+
+function intersect2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...)
+    return clip_layerwise(Clipper.ClipTypeIntersection, obj, tool; kwargs...)
+end
+
+function union2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...)
+    return clip_layerwise(Clipper.ClipTypeUnion, obj, tool; kwargs...)
+end
+
+for func in ("union2d", "difference2d", "intersect2d", "xor2d")
+    func_layerwise = Symbol(string(func) * "_layerwise")
+    doc = """
+        $(func)_layerwise(obj::GeometryStructure, tool::GeometryStructure;
+            only_layers=[],
+            ignore_layers=[],
+            depth=-1,
+            tiled=false,
+            tile_size=nothing
+        )
+
+    Return a `Dict` of `meta => [$(func)(obj => meta, tool => meta)]` for each unique element metadata `meta` in `obj` and `tool`.
+
+    Entities with metadata matching `only_layers` or `ignore_layers` are included or excluded based on [`layer_inclusion`](@ref).
+
+    Entities in references up to a depth of `depth` are included, where `depth=0` uses only top-level entities in `obj` and `tool`.
+    Depth is unlimited by default.
+
+    See also [`$(func)`](@ref).
+
+    # Tiling
+
+    Using `tiled=true` or manually setting a `tile_size` can significantly speed up operations and reduce maximum memory usage for large geometries.
+    It does this by breaking up the geometry into smaller portions ("tiles") and operating on them one at a time.
+
+    If a length is provided to `tile_size`, the bounds of the combined geometries are tiled with squares with that
+    edge length, starting from the lower left corner.
+
+    If `tiled` is `true` but `tile_size` is not specified, a tile size will be set automatically based on the total number of entities in the operation,
+    such that there is about one square tile per 100 entities. This is usually a reasonable choice, but you may want to benchmark your use case.
+
+    Entities crossing between tiles are split into their intersections with each tile before clipping.
+    For each tile, those intersection results and all entities inside that tile are selected.
+    The values for each layer in the returned `Dict` are then lazy iterators over clipping results for
+    selected entities in each tile.
+    """
+    eval(quote
+        @doc $doc $func_layerwise
+    end)
+end
+
+function clip_layerwise(
+    op::Clipper.ClipType,
+    obj::GeometryStructure,
+    tool::GeometryStructure;
+    only_layers=[],
+    ignore_layers=[],
+    tiled=false,
+    tile_size=nothing,
+    depth=-1,
+    pfs=Clipper.PolyFillTypePositive,
+    pfc=Clipper.PolyFillTypePositive
+)
+    metadata_filter = DeviceLayout.layer_inclusion(only_layers, ignore_layers)
+    if metadata_filter == DeviceLayout.trivial_inclusion
+        metadata_filter = nothing
+    end
+    obj_flat = DeviceLayout.flatten(obj; metadata_filter, depth)
+    tool_flat = DeviceLayout.flatten(tool; metadata_filter, depth)
+    obj_metas = unique(DeviceLayout.element_metadata(obj_flat))
+    tool_metas = unique(DeviceLayout.element_metadata(tool_flat))
+    all_metas = unique([obj_metas; tool_metas])
+    if !tiled && isnothing(tile_size)
+        res = Dict(
+            meta => [
+                clip(
+                    op,
+                    DeviceLayout.elements(obj_flat)[DeviceLayout.element_metadata(
+                        obj_flat
+                    ) .== meta],
+                    DeviceLayout.elements(tool_flat)[DeviceLayout.element_metadata(
+                        tool_flat
+                    ) .== meta],
+                    pfs=pfs,
+                    pfc=pfc
+                )
+            ] for meta in all_metas
+        )
+    else
+        res = Dict(
+            meta => clip_tiled(
+                op,
+                DeviceLayout.elements(obj_flat)[DeviceLayout.element_metadata(
+                    obj_flat
+                ) .== meta],
+                DeviceLayout.elements(tool_flat)[DeviceLayout.element_metadata(
+                    tool_flat
+                ) .== meta],
+                tile_size,
+                pfs=pfs,
+                pfc=pfc
+            ) for meta in all_metas
+        )
+    end
+    return res
+end
+
+### Tiling
+function tiles_and_edges(r::Rectangle, tile_size)
+    nx = ceil(width(r) / tile_size)
+    ny = ceil(height(r) / tile_size)
+    tile0 = r.ll + Rectangle(tile_size, tile_size)
+    tiles =
+        [tile0 + Point((i - 1) * tile_size, (j - 1) * tile_size) for i = 1:nx for j = 1:ny]
+    h_edges = [
+        Rectangle(
+            Point(r.ll.x, r.ll.y + (i - 1) * tile_size),
+            Point(r.ur.x, r.ll.y + (i - 1) * tile_size)
+        ) for i = 2:ny
+    ]
+    v_edges = [
+        Rectangle(
+            Point(r.ll.x + (i - 1) * tile_size, r.ll.y),
+            Point(r.ll.x + (i - 1) * tile_size, r.ur.y)
+        ) for i = 2:nx
+    ]
+    return tiles, vcat(h_edges, v_edges) # Edges could be used for healing
+end
+
+function _auto_tile_size(bnds::Rectangle, num_ents)
+    target_ents_per_tile = 100
+    target_num_tiles = max(1.0, round(num_ents / target_ents_per_tile))
+    target_size = sqrt(width(bnds) * height(bnds) / target_num_tiles)
+    # Return the closest tile size that evenly divides width
+    return width(bnds) / max(1.0, round(width(bnds) / target_size))
+end
+
+"""
+    function clip_tiled(
+        op,
+        ents1::AbstractArray{<:GeometryEntity{T}},
+        ents2::AbstractArray{<:GeometryEntity{T}},
+        tile_size=nothing;
+        pfs=Clipper.PolyFillTypePositive,
+        pfc=Clipper.PolyFillTypePositive
+    )
+
+Return a lazy iterator that applies `op(ents1, ents2)` tile by tile.
+
+The bounds of the combined geometries are tiled with squares with edge length `tile_size`, starting at the bottom left
+corner. Entities crossing between tiles are split into their intersections with each tile before clipping.
+For each tile, those intersection results and all entities inside that tile are selected.
+The return value is a lazy iterator over clipping results for selected entities per tile.
+
+A rough guideline for choosing tile size is to aim for 100 polygons per tile, but you may want to
+benchmark your use case. If an explicit tile size is not provided, then tile size will be set automatically
+based on the total number of entities in the operation, such that there is about one square tile per 100 entities.
+"""
+function clip_tiled(
+    op,
+    ents1::AbstractArray{<:GeometryEntity{T}},
+    ents2::AbstractArray{<:GeometryEntity{T}},
+    tile_size=nothing;
+    pfs=Clipper.PolyFillTypePositive,
+    pfc=Clipper.PolyFillTypePositive
+) where {T}
+    (isempty(ents1) && isempty(ents2)) &&
+        return Iterators.map(identity, ClippedPolygon{T}[])
+    # Create spatial index for each set of polygons
+    if isempty(ents2)
+        tree1 = DeviceLayout.mbr_spatial_index(ents1)
+        bnds = mbr(tree1)
+    elseif isempty(ents1)
+        tree2 = DeviceLayout.mbr_spatial_index(ents2)
+        bnds = mbr(tree2)
+    else
+        tree1 = DeviceLayout.mbr_spatial_index(ents1)
+        tree2 = DeviceLayout.mbr_spatial_index(ents2)
+        bnds = SpatialIndexing.combine(mbr(tree1), mbr(tree2))
+    end
+
+    # Get tiles and indices of polygons intersecting tiles
+    bnds_dl = Rectangle(
+        Point(bnds.low...) * DeviceLayout.onemicron(T),
+        Point(bnds.high...) * DeviceLayout.onemicron(T)
+    )
+
+    if isnothing(tile_size)
+        tile_size = _auto_tile_size(bnds_dl, length(ents1) + length(ents2))
+    end
+
+    tiles, edges = tiles_and_edges(bnds_dl, tile_size) # DeviceLayout Rectangles
+    # Get single vector of all entity indices touching any edge
+    edge_touching_idx1 = Set( # Use Set because we'll be testing membership
+        mapreduce(vcat, edges, init=Int[]) do edge
+            return isempty(ents1) ? Int[] : DeviceLayout.findbox(edge, tree1; intersects=true)
+        end
+    )
+    # Same for ents2
+    edge_touching_idx2 = Set(
+        mapreduce(vcat, edges, init=Int[]) do edge
+            return isempty(ents2) ? Int[] : DeviceLayout.findbox(edge, tree2; intersects=true)
+        end
+    )
+    # Get vector of (ents1 indices, ents2 indices) for each tile
+    tile_poly_indices = map(tiles) do tile
+        idx1 = isempty(ents1) ? Int[] : DeviceLayout.findbox(tile, tree1; intersects=true)
+        idx2 =
+            isempty(ents2) ? Int[] : DeviceLayout.findbox(tile, tree2; intersects=true)
+        return (idx1, idx2)
+    end
+    # Clip within each tile
+    res = Iterators.map(enumerate(tile_poly_indices)) do (tile_idx, poly_idxs)
+        idx1, idx2 = poly_idxs
+        # If an entity is touching a tile edge, clip it to the tile with intersect2d
+        idx1_on_edge = in.(idx1, Ref(edge_touching_idx1))
+        idx2_on_edge = in.(idx2, Ref(edge_touching_idx2))
+        edge_idx1 = idx1[idx1_on_edge]
+        bulk_idx1 = idx1[(!).(idx1_on_edge)]
+        edge_idx2 = idx2[idx2_on_edge]
+        bulk_idx2 = idx2[(!).(idx2_on_edge)]
+        tile = tiles[tile_idx]
+        ents1_clipped_to_tile =
+            vcat(ents1[bulk_idx1], intersect2d(tile, ents1[edge_idx1]))
+        ents2_clipped_to_tile =
+            vcat(ents2[bulk_idx2], intersect2d(tile, ents2[edge_idx2]))
+        return clip(op, ents1_clipped_to_tile, ents2_clipped_to_tile; pfs, pfc)
+    end
+    return res
+end
+
 end # module
diff --git a/src/schematics/ExamplePDK/components/ReadoutResonators/clawed_meander.jl b/src/schematics/ExamplePDK/components/ReadoutResonators/clawed_meander.jl
index 5a3980fb..01a57fa5 100644
--- a/src/schematics/ExamplePDK/components/ReadoutResonators/clawed_meander.jl
+++ b/src/schematics/ExamplePDK/components/ReadoutResonators/clawed_meander.jl
@@ -1,6 +1,6 @@
 import DeviceLayout: flushtop, flushleft, flushright, below, above
 
-import .ExamplePDK: filter_params, tap!
+import .ExamplePDK: tap!
 import .ExamplePDK.Transmons: ExampleRectangleTransmon
 
 """
diff --git a/src/solidmodels/render.jl b/src/solidmodels/render.jl
index 8f9d8011..6aa734f1 100644
--- a/src/solidmodels/render.jl
+++ b/src/solidmodels/render.jl
@@ -296,7 +296,7 @@ function round_to_curvilinearpolygon(
                 push!(new_points, Paths.p1(result.fillet))
                 # Record trim for the original arc
                 arc_start_vtx = arc_is_outgoing ? i : mod1(i - 1, len)
-                curve_k = findfirst(csi -> abs(csi) == arc_start_vtx, pol.curve_start_idx)
+                curve_k = findfirst(isequal(arc_start_vtx), pol.curve_start_idx)
                 if !isnothing(curve_k)
                     if arc_is_outgoing
                         trim_start_pts[curve_k] = result.T_arc
@@ -366,7 +366,7 @@ function round_to_curvilinearpolygon(
 
     # 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, by=abs)
+        perm = sortperm(new_curve_start_idx)
         new_curves = new_curves[perm]
         new_curve_start_idx = new_curve_start_idx[perm]
     end
@@ -1652,17 +1652,9 @@ function _add_loop!(
     endpoint_pairs = zip(pts, circshift(pts, -1))
     curves = map(enumerate(endpoint_pairs)) do (i, endpoints)
         curve_idx = findfirst(isequal(i), cl.curve_start_idx)
-        if isnothing(curve_idx) # not found
-            # see if the negative of the index is in there
-            curve_idx = findfirst(isequal(-i), cl.curve_start_idx)
-            if isnothing(curve_idx) # nope, just a line
-                return k.add_line(endpoints[1], endpoints[2])
-            else # negative index => reverse endpoints
-                c = _add_curve!(reverse(endpoints), cl.curves[curve_idx], k, z; atol)
-                !isempty(size(c)) && return reverse(c)
-                return c
-            end
-        else # add the curve whose start index is i
+        if isnothing(curve_idx)
+            return k.add_line(endpoints[1], endpoints[2])
+        else
             return _add_curve!(endpoints, cl.curves[curve_idx], k, z; atol)
         end
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index 11fdc5b2..7718b478 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,6 +23,18 @@ using TestItemRunner
     import Unitful: pm, nm, μm, mm, cm, m
 
     p(x, y) = Point(x, y)
+
+    """
+        is_sliver(p::Polygon{T}; atol=DeviceLayout.onenanometer(T))
+
+    Return `true` if `2 * area(p) / perimeter(p) < atol`, and `false` otherwise.
+
+    In other words, if `p` has an "average width" less than `atol`, it is counted as a sliver.
+    """
+    function is_sliver(p::Polygon{T}; atol=DeviceLayout.onenanometer(T)) where {T}
+        return 2 * Polygons.area(p) / Polygons.perimeter(p) < atol
+    end
+
     const tdir = mktempdir()
 
     # G1 continuity check: verify no angle jump exceeds the discretization step
diff --git a/test/test_bspline.jl b/test/test_bspline.jl
index 2663789f..1b856fe0 100644
--- a/test/test_bspline.jl
+++ b/test/test_bspline.jl
@@ -180,13 +180,8 @@ end
     approx = Paths.bspline_approximation(c, atol=100.0nm)
     pts = DeviceLayout.discretize_curve(c, 100.0nm)
     pts_approx = vcat(DeviceLayout.discretize_curve.(approx.segments, 100.0nm)...)
-    area(p) =
-        sum(
-            (gety.(p.p) + gety.(circshift(p.p, -1))) .*
-            (getx.(p.p) - getx.(circshift(p.p, -1)))
-        ) / 2
     poly = Polygon([pts; reverse(pts_approx)])
-    @test abs(area(poly) / perimeter(poly)) < 100nm # It's actually ~25nm but the guarantee is ~< tolerance
+    @test abs(Polygons.area(poly) / perimeter(poly)) < 100nm # It's actually ~25nm but the guarantee is ~< tolerance
     c = curv[8].curves[1] # ConstantOffset Turn
     @test Paths.curvatureradius(c, 10μm) == sign(c.seg.α) * c.seg.r - c.offset
     @test curvatureradius_fd(c, 10μm) ≈ Paths.curvatureradius(c, 10μm) atol = 1nm
diff --git a/test/test_clipping.jl b/test/test_clipping.jl
index 05f33034..ab2ec771 100644
--- a/test/test_clipping.jl
+++ b/test/test_clipping.jl
@@ -68,12 +68,28 @@
             small_square = Polygon(square_pts)
             big_square = Polygon(2 .* square_pts)
             mask = first(
-                to_polygons(clip(Clipper.ClipTypeDifference, big_square, small_square))
+                to_polygons(
+                    clip(
+                        Clipper.ClipTypeDifference,
+                        big_square,
+                        small_square,
+                        pfs=Clipper.PolyFillTypeEvenOdd,
+                        pfc=Clipper.PolyFillTypeEvenOdd
+                    )
+                )
             )
             c = Cell{Float64}("test")
             render!(c, mask)
             canvas = Polygon(5 .* square_pts)
-            x = to_polygons(clip(Clipper.ClipTypeDifference, canvas, mask))
+            x = to_polygons(
+                clip(
+                    Clipper.ClipTypeDifference,
+                    canvas,
+                    mask,
+                    pfs=Clipper.PolyFillTypeEvenOdd,
+                    pfc=Clipper.PolyFillTypeEvenOdd
+                )
+            )
             c2 = Cell{Float64}("squares")
             for z in x
                 render!(c2, z, GDSMeta(0))
@@ -85,11 +101,23 @@
             square_pts = Point.([(-1, -1), (-1, 1), (1, 1), (1, -1)])
             small_square = Polygon(square_pts)
             big_square = Polygon(2 .* square_pts)
-            mask = clip(Clipper.ClipTypeDifference, big_square, small_square)
+            mask = clip(
+                Clipper.ClipTypeDifference,
+                big_square,
+                small_square,
+                pfs=Clipper.PolyFillTypeEvenOdd,
+                pfc=Clipper.PolyFillTypeEvenOdd
+            )
             c = Cell{Float64}("test")
             render!(c, mask)
             canvas = Polygon(5 .* square_pts)
-            x = clip(Clipper.ClipTypeDifference, canvas, mask)
+            x = clip(
+                Clipper.ClipTypeDifference,
+                canvas,
+                mask,
+                pfs=Clipper.PolyFillTypeEvenOdd,
+                pfc=Clipper.PolyFillTypeEvenOdd
+            )
             c2 = Cell{Float64}("squares")
             render!(c2, x, GDSMeta(0))
 
@@ -479,6 +507,23 @@
         pp = Point.([(185.0, -100.0), (300.0, -215.0), (300.0, -185.0), (215.0, -100.0)])
         @test Polygons.orientation(Polygon(pp)) == 1
     end
+
+    @testset "> Mixed argument types" begin
+        p1 = Rectangle(1mm, 1mm)
+        p2 = centered(Rectangle(1µm, 1µm))
+        # Clipped polygons with mixed units
+        p1_clip = union2d(p1)
+        p2_clip = union2d(p2)
+        @test union2d([p1_clip, p2_clip]) == union2d([p1, p2])
+        # add in a GeometryStructure
+        cs = CoordinateSystem("test")
+        place!(cs, p2, :test)
+        @test coordinatetype(cs => :test) == coordinatetype(cs)
+        @test union2d(cs => :test) == p2_clip
+        @test union2d([p1, cs => :test]) == union2d([p1, p2])
+        @test union2d(p1, cs => :test) == union2d([p1, p2])
+        @test union2d(p1, [cs => :test]) == union2d([p1, p2])
+    end
 end
 
 @testitem "Polygon offsetting" setup = [CommonTestSetup] begin
@@ -923,13 +968,7 @@ end
 
 @testitem "Hole sorting" setup = [CommonTestSetup] begin
     # Issue #175
-    # Shoelace formula for signed polygon area
-    area(p) =
-        sum(
-            (gety.(p.p) + gety.(circshift(p.p, -1))) .*
-            (getx.(p.p) - getx.(circshift(p.p, -1)))
-        ) / 2
-
+    import .Polygons: area
     @testset "Hole sorting" begin
         c1 = rotate(centered(Rectangle(0.5mm, 0.5mm)), 45°) - Point(0mm, 0.5mm)
         c2 = c1 + Point(0mm, 1mm)
@@ -970,3 +1009,96 @@ end
         @test area(off[2]) ≈ (3.1mm)^2 - 2 * (0.4mm)^2
     end
 end
+
+@testitem "Layerwise clipping" setup = [CommonTestSetup] begin
+    c1 = CoordinateSystem("test1")
+    c2 = CoordinateSystem("test2")
+    r1 = Rectangle(10μm, 10μm)
+    r2 = r1 + Point(5μm, 5μm)
+    overlap = intersect2d(r1, r2)
+    x = xor2d(r1, r2)
+    d1 = difference2d(r1, r2)
+    d2 = difference2d(r2, r1)
+    uni = union2d(r1, r2)
+    lyr_a = SemanticMeta(:a)
+    lyr_b = SemanticMeta(:b)
+    place!(c1, r1, lyr_a)
+    place!(c1, r2, lyr_b)
+    place!(c2, r1, lyr_b)
+    place!(c2, r2, lyr_a)
+
+    @test isempty(to_polygons(xor2d_layerwise(c1, c1)[lyr_a][1]))
+    @test isempty(to_polygons(xor2d_layerwise(c1, c1)[lyr_b][1]))
+    @test xor2d_layerwise(c1, c2)[lyr_a][1] == x
+    @test xor2d_layerwise(c1, c2)[lyr_b][1] == x
+    @test difference2d_layerwise(c1, c2)[lyr_a][1] == d1
+    @test difference2d_layerwise(c1, c2)[lyr_b][1] == d2
+    @test union2d_layerwise(c1, c2)[lyr_a][1] == uni
+    @test intersect2d_layerwise(c1, c2)[lyr_b][1] == overlap
+
+    # Findbox
+    @test DeviceLayout.findbox(r1, [r1, r2]) == [1]
+    @test DeviceLayout.findbox(r1, [r1, r2]; intersects=true) == [1, 2]
+    @test DeviceLayout.findbox(r1, [c1]) == []
+    @test DeviceLayout.findbox(r1, [r1, c1], intersects=true) == [1, 2]
+
+    # Tiling
+    ca_1 = CoordinateSystem("array1")
+    ca_2 = CoordinateSystem("array2")
+    addref!(ca_1, aref(c1, 100μm * (-1:1), 100μm * (-1:1)))
+    addref!(ca_2, aref(c2, 100μm * (-1:1), 100μm * (-1:1)))
+    xa = xor2d_layerwise(ca_1, ca_2)
+    xa_tiled = xor2d_layerwise(ca_1, ca_2, tile_size=99μm2nm)
+    @test length(xa_tiled[lyr_a]) == 3 * 3
+    @test length(xa_tiled[lyr_b]) == 3 * 3
+    @test length(vcat(to_polygons.(xa_tiled[lyr_a])...)) == 3 * 3 * 2
+    @test isempty(to_polygons(xor2d(vcat(xa_tiled[lyr_a]...), xa[lyr_a])))
+
+    # Tiling with entities on edges
+    xa_tiled_edges = xor2d_layerwise(ca_1, ca_2, tile_size=106μm2nm)
+    @test length(xa_tiled_edges[lyr_a]) == 3 * 3
+    @test length(xa_tiled_edges[lyr_b]) == 3 * 3
+    all_polys = vcat(to_polygons.(xa_tiled_edges[lyr_a])...)
+    # EvenOdd union to remove regions where polygons overlap
+    all_no_overlap = clip(
+        Polygons.Clipper.ClipTypeUnion,
+        all_polys,
+        Polygon{typeof(1.0nm)}[],
+        pfs=Polygons.Clipper.PolyFillTypeEvenOdd,
+        pfc=Polygons.Clipper.PolyFillTypeEvenOdd
+    )
+    @test length(all_polys) > 3 * 3 * 2 # Some polygons were split
+    @test isempty(to_polygons(xor2d(all_polys, xa[lyr_a]))) # Split polygons are still correct
+    @test length(to_polygons(all_no_overlap)) == 3 * 3 * 2 # Split polygons are not overlapping
+    @test isempty(to_polygons(xor2d(all_no_overlap, xa[lyr_a]))) # Split polygons add up correctly
+
+    # Tiling with empty layers
+    uae = union2d_layerwise(ca_1, CoordinateSystem("empty"))
+    uae_tiled = union2d_layerwise(ca_1, CoordinateSystem("empty"), tile_size=99μm2nm)
+    @test length(vcat(to_polygons.(uae_tiled[lyr_a])...)) == 3 * 3
+    @test isempty(to_polygons(xor2d(vcat(uae_tiled[lyr_a]...), ca_1 => lyr_a)))
+    uae_tiled = union2d_layerwise(CoordinateSystem("empty"), ca_1, tile_size=99μm2nm)
+    @test length(vcat(to_polygons.(uae_tiled[lyr_a])...)) == 3 * 3
+    @test isempty(to_polygons(xor2d(vcat(uae_tiled[lyr_a]...), ca_1 => lyr_a)))
+
+    # Auto tile size
+    ca_3 = CoordinateSystem("array1")
+    ca_4 = CoordinateSystem("array2")
+    addref!(ca_3, aref(c1, 100μm * (-10:10), 100μm * (-10:10)))
+    addref!(ca_4, aref(c2, 100μm * (-10:10), 100μm * (-10:10)))
+    xa_auto_tiled = xor2d_layerwise(ca_3, ca_4, tiled=true)
+    @test all(length.(values(xa_auto_tiled)) .== 9) # 9 tiles for about 900 entities
+
+    # Auto tile size is decent for large n
+    for n in [12000, 20000, 33000]
+        w, l, n = (10mm, 4mm, n)
+        l_tile = Polygons._auto_tile_size(Rectangle(w, l), n)
+        num_tiles = ceil(w / l_tile) * ceil(l / l_tile)
+        @test abs(100 - n / num_tiles) <= 10
+    end
+    # Works for sub-single-tile
+    w, l, n = (400μm, 400μm, 40)
+    l_tile = Polygons._auto_tile_size(Rectangle(w, l), n)
+    num_tiles = ceil(w / l_tile) * ceil(l / l_tile)
+    @test num_tiles == 1
+end
diff --git a/test/test_line_arc_rounding.jl b/test/test_line_arc_rounding.jl
index 7e0b7ca0..6c805780 100644
--- a/test/test_line_arc_rounding.jl
+++ b/test/test_line_arc_rounding.jl
@@ -447,7 +447,7 @@
         [neg_arc],
         [-3]  # negative: curve between p[3] and p[4], parameterized from p[4] to p[3]
     )
-    @test neg_cp.curve_start_idx[1] < 0
+    @test neg_cp.curve_start_idx[1] > 0  # constructor normalizes to positive
 
     # Plain rendering must work (catches invalid vertex/curve mismatch)
     neg_plain = to_polygons(neg_cp)
diff --git a/test/test_render.jl b/test/test_render.jl
index fbf6b3f0..da29bba0 100644
--- a/test/test_render.jl
+++ b/test/test_render.jl
@@ -1005,6 +1005,11 @@ end
         @test Polygons.orientation(halo(dr1, 0.1μm)[2]) == -1 # still has a hole
         @test footprint(union2d(r1, r1 + Point(40, 0)μm)) isa Rectangle # multipolygon => use bounds
         @test halo(union2d(r3), 1μm, -0.5μm) == dr1 # ClippedPolygon halo with inner delta
+
+        @test Polygons.area(to_polygons(difference2d(r3, r4))[1]) ==
+              Polygons.area(to_polygons(r3)) - Polygons.area(to_polygons(r4))
+        @test is_sliver(to_polygons(difference2d(r3, r3 + Point(5, 5)nm))[1]; atol=10nm)
+        @test (!).(is_sliver(to_polygons(difference2d(r3, r3 + Point(5, 5)nm))[1])) # default 1nm
     end
 end
 
diff --git a/test/test_shapes.jl b/test/test_shapes.jl
index 6840e003..75fd48d9 100644
--- a/test/test_shapes.jl
+++ b/test/test_shapes.jl
@@ -429,6 +429,9 @@ end
     c = Cell("main", nm)
     @test_nowarn render!(c, cs)
 
+    # Constructor normalizes negative curve_start_idx to positive via reverse()
+    @test cp.curve_start_idx[1] == 2
+
     # Reverse parameterized and forward parameterized should produce same number of points
     @test length(points(to_polygons(cp))) == length(pgen)
     @test length(points(to_polygons(t(cp)))) == length(ptgen)
@@ -438,10 +441,22 @@ end
     c = Cell("main", nm)
     @test_nowarn render!(c, cs)
 
-    # Clipping the transformed inverse and forward should give negligible difference.
-    # Adaptive discretization may produce thin slivers rather than exactly empty.
-    diff_poly = difference2d(to_polygons(cpt), to_polygons(t(cp)))
-    @test perimeter(diff_poly) < 0.1μm
+    # Transformed forward and reverse should give identical geometry.
+    # (Constructor normalization makes them internally equivalent.)
+    cpt_pts = points(to_polygons(cpt))
+    tcp_pts = points(to_polygons(t(cp)))
+    @test length(cpt_pts) == length(tcp_pts)
+    @test all(isapprox.(cpt_pts, tcp_pts; atol=0.01nm))
+
+    # Negative csi must be normalized before duplicate-point cleanup, otherwise
+    # -3 ∉ [3] and the curve between duplicate points survives deletion.
+    dup_pts =
+        [Point(0.0μm, 0.0μm), Point(1.0μm, 0.0μm), Point(0.5μm, 1.0μm), Point(0.5μm, 1.0μm)]
+    dup_seg = Paths.Turn(90°, 1.0μm, α0=0°, p0=dup_pts[4])
+    dup_cp = CurvilinearPolygon(dup_pts, [dup_seg], [-3])
+    @test length(dup_cp.p) == 3
+    @test isempty(dup_cp.curves)
+    @test isempty(dup_cp.curve_start_idx)
 
     # Convert a SimpleTrace to a CurvilinearRegion
     pa = Path(0nm, 0nm)
@@ -609,14 +624,9 @@ end
     @test length(points(circ_default)) == length(circular_arc(2pi, 1.0μm, 1.0nm)) - 1 # last pt duplicated
 
     # Make sure error is as small as tolerance says
-    area(p) =
-        sum(
-            (gety.(p.p) + gety.(circshift(p.p, -1))) .*
-            (getx.(p.p) - getx.(circshift(p.p, -1)))
-        ) / 2
     e_fine = to_polygons(e; atol=0.1nm)
     poly = to_polygons(difference2d(e_fine, e_default))[1]
-    @test abs(area(poly) / perimeter(poly)) < 1nm # on average better than 1nm
+    @test is_sliver(poly; atol=1nm) # on average better than 1nm (area/perimeter)
 
     # Last two points are not too close together
     poly = points(to_polygons(e, atol=60nm))
diff --git a/test/tests.jl b/test/tests.jl
index ec4e09b4..b298c86d 100644
--- a/test/tests.jl
+++ b/test/tests.jl
@@ -200,6 +200,11 @@ end
         @test perimeter(pfloat) ≈ (2 + sqrt(2))m
         badpoly = Polygon(Point{Float64}[])
         @test perimeter(badpoly) == 0.0
+
+        @test Polygons.signed_area(pfloat) == 0.5m^2
+        @test Polygons.signed_area(Polygon(reverse(points(pfloat)))) == -0.5m^2
+        @test Polygons.area(Polygon(reverse(points(pfloat)))) == 0.5m^2
+        @test (!).(is_sliver(pfloat))
     end
 end