Add layerwise and lazy tiled clipping

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`)

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",

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
 ```
 

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

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,26 +548,22 @@ 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)
     prev_i = mod1(i - 1, n)
 
     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

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