From 4636427014ff46c8b7c7e9e90b903778826edafc Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Mon, 13 Apr 2026 17:27:06 +0200 Subject: [PATCH 01/18] Normalize curvilinear curve_idx to be positive --- src/curvilinear.jl | 7 +++++++ test/test_line_arc_rounding.jl | 2 +- test/test_shapes.jl | 13 +++++++++---- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/src/curvilinear.jl b/src/curvilinear.jl index ebb97aa0..7d53f36a 100644 --- a/src/curvilinear.jl +++ b/src/curvilinear.jl @@ -55,6 +55,13 @@ struct CurvilinearPolygon{T} <: GeometryEntity{T} # 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] 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 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_shapes.jl b/test/test_shapes.jl index 6840e003..f98e8a9a 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,12 @@ 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)) # Convert a SimpleTrace to a CurvilinearRegion pa = Path(0nm, 0nm) From 19aae5b9ab3754be22655df482c6b48c09db4fc4 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Mon, 13 Apr 2026 18:59:33 +0200 Subject: [PATCH 02/18] Remove dead curve start index sign code --- src/curvilinear.jl | 56 ++++++++++++++++----------------------- src/solidmodels/render.jl | 18 ++++--------- test/test_shapes.jl | 10 +++++++ 3 files changed, 38 insertions(+), 46 deletions(-) diff --git a/src/curvilinear.jl b/src/curvilinear.jl index 7d53f36a..2d3268bc 100644 --- a/src/curvilinear.jl +++ b/src/curvilinear.jl @@ -52,8 +52,8 @@ 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) @@ -107,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]) @@ -471,7 +471,7 @@ end # 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)...) + 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 @@ -549,9 +549,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) @@ -559,16 +559,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) @@ -603,7 +599,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. @@ -703,19 +699,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. @@ -732,7 +722,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/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/test_shapes.jl b/test/test_shapes.jl index f98e8a9a..1f19187d 100644 --- a/test/test_shapes.jl +++ b/test/test_shapes.jl @@ -448,6 +448,16 @@ end @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) straight!(pa, 100μm, Paths.SimpleTrace(10.0μm)) From 2c7a9ee0e891c229ddd0fd496f1419f0f9609a90 Mon Sep 17 00:00:00 2001 From: Greg Peairs Date: Tue, 28 Oct 2025 19:11:52 +0000 Subject: [PATCH 03/18] Refactor clipping into separate file Move clipping to its own file --- src/clipping.jl | 1203 +++++++++++++++++++++++++++++++++++++++++++++++ src/polygons.jl | 4 + 2 files changed, 1207 insertions(+) create mode 100644 src/clipping.jl diff --git a/src/clipping.jl b/src/clipping.jl new file mode 100644 index 00000000..af2276cd --- /dev/null +++ b/src/clipping.jl @@ -0,0 +1,1203 @@ +# Clipping polygons one at a time +""" + clip(op::Clipper.ClipType, s, c; kwargs...) where {S<:Coordinate, T<:Coordinate} + clip(op::Clipper.ClipType, s::AbstractVector{A}, c::AbstractVector{B}; + 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} + +Return the `ClippedPolygon` resulting from a polygon clipping operation. + +Uses the [`Clipper`](http://www.angusj.com/delphi/clipper.php) library and the +[`Clipper.jl`](https://github.com/Voxel8/Clipper.jl) wrapper to perform polygon clipping. + +## Positional arguments + +The first argument must be one of the following types to specify a clipping operation: + + - `Clipper.ClipTypeDifference` + - `Clipper.ClipTypeIntersection` + - `Clipper.ClipTypeUnion` + - `Clipper.ClipTypeXor` + +Note that these are types; you should not follow them with `()`. + +The second and third argument may be a `GeometryEntity` or array of `GeometryEntity`. All entities +are first converted to polygons using [`to_polygons`](@ref). +Each can also be a `GeometryStructure` or `GeometryReference`, in which case +`elements(flatten(p))` will be converted to polygons. +Each can also be a pair `geom => layer`, where `geom` is a +`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection +of either, in which case only the elements in those layers will be taken from the flattened structure. + +## Keyword arguments + +`pfs` and `pfc` specify polygon fill rules for the `s` and `c` arguments, respectively. +These arguments may include: + + - `Clipper.PolyFillTypeNegative` + - `Clipper.PolyFillTypePositive` + - `Clipper.PolyFillTypeEvenOdd` + - `Clipper.PolyFillTypeNonZero` + +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). +""" +function clip(op::Clipper.ClipType, s, c; kwargs...) + return clip(op, _normalize_clip_arg(s), _normalize_clip_arg(c); kwargs...) +end + +# Clipping requires an AbstractVector{Polygon{T}} +_normalize_clip_arg(p::Polygon) = [p] +_normalize_clip_arg(p::GeometryEntity) = _normalize_clip_arg(to_polygons(p)) +_normalize_clip_arg(p::AbstractArray{Polygon{T}}) where {T} = p +_normalize_clip_arg(p::AbstractArray{<:GeometryEntity{T}}) where {T} = + reduce(vcat, to_polygons.(p); init=Polygon{T}[]) +_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)) + +# Clipping arrays of AbstractPolygons +function clip( + op::Clipper.ClipType, + s::AbstractVector{A}, + c::AbstractVector{B}; + kwargs... +) where {S, T, A <: Polygon{S}, B <: Polygon{T}} + dimension(S) != dimension(T) && throw(Unitful.DimensionError(oneunit(S), oneunit(T))) + R = promote_type(S, T) + + return clip( + op, + convert(Vector{Polygon{R}}, s), + convert(Vector{Polygon{R}}, c); + kwargs... + ) +end + +# Clipping two identically-typed arrays of <: Polygon +function clip( + op::Clipper.ClipType, + s::AbstractVector{Polygon{T}}, + c::AbstractVector{Polygon{T}}; + pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, + pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, + tiled=false +) where {T} + sc, cc = clipperize(s), clipperize(c) + polys = if !tiled + _clip(op, sc, cc; pfs, pfc) + else + _clip_tiled(op, sc, cc; pfs, pfc) + end + return declipperize(polys, T) +end + +""" + cliptree(op::Clipper.ClipType, s::AbstractPolygon{S}, c::AbstractPolygon{T}; + kwargs...) where {S<:Coordinate, T<:Coordinate} + cliptree(op::Clipper.ClipType, s::AbstractVector{A}, c::AbstractVector{B}; + 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} + +Return a `Clipper.PolyNode` representing parent-child relationships between polygons and +interior holes. The units and number type may need to be converted. + +Uses the [`Clipper`](http://www.angusj.com/delphi/clipper.php) library and the +[`Clipper.jl`](https://github.com/Voxel8/Clipper.jl) wrapper to perform polygon clipping. + +## Positional arguments + +The first argument must be one of the following types to specify a clipping operation: + + - `Clipper.ClipTypeDifference` + - `Clipper.ClipTypeIntersection` + - `Clipper.ClipTypeUnion` + - `Clipper.ClipTypeXor` + +Note that these are types; you should not follow them with `()`. The second and third +arguments are `AbstractPolygon`s or vectors thereof. + +## Keyword arguments + +`pfs` and `pfc` specify polygon fill rules for the `s` and `c` arguments, respectively. +These arguments may include: + + - `Clipper.PolyFillTypeNegative` + - `Clipper.PolyFillTypePositive` + - `Clipper.PolyFillTypeEvenOdd` + - `Clipper.PolyFillTypeNonZero` + +See the [`Clipper` docs](http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Types/PolyFillType.htm) +for further information. +""" +function cliptree( + op::Clipper.ClipType, + s::AbstractPolygon{S}, + c::AbstractPolygon{T}; + kwargs... +) where {S <: Coordinate, T <: Coordinate} + dimension(S) != dimension(T) && throw(Unitful.DimensionError(oneunit(S), oneunit(T))) + R = promote_type(S, T) + return cliptree(op, Polygon{R}[s], Polygon{R}[c]; kwargs...)::Vector{Polygon{R}} +end + +function cliptree( + op::Clipper.ClipType, + s::AbstractVector{A}, + c::AbstractVector{B}; + kwargs... +) where {S, T, A <: AbstractPolygon{S}, B <: AbstractPolygon{T}} + dimension(S) != dimension(T) && throw(Unitful.DimensionError(oneunit(S), oneunit(T))) + R = promote_type(S, T) + return cliptree( + op, + convert(Vector{Polygon{R}}, s), + convert(Vector{Polygon{R}}, c); + kwargs... + )::Vector{Polygon{R}} +end + +function cliptree( + op::Clipper.ClipType, + s::AbstractVector{Polygon{T}}, + c::AbstractVector{Polygon{T}}; + pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, + pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd +) where {T} + sc, cc = clipperize(s), clipperize(c) + cpoly = _clip(op, sc, cc; pfs, pfc) + return declipperize(cpoly, T).tree +end + +""" + union2d(p1, p2) + +Return the geometric union of p1 and p2 as a `ClippedPolygon`. + +Each of `p1` and `p2` may be a `GeometryEntity` or array of `GeometryEntity`. All entities +are first converted to polygons using [`to_polygons`](@ref). + +Each of `p1` and `p2` can also be a `GeometryStructure` or `GeometryReference`, in which case +`elements(flatten(p))` will be converted to polygons. + +Each can also be a pair `geom => layer`, where `geom` is a +`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection +of either, in which case only the elements in those layers will used. + +This is not implemented as a method of `union` because you can have a set union of arrays of +polygons, which is a distinct operation. + +The Clipper polyfill rule is PolyFillTypePositive, meaning as long as a +region lies within more non-hole (by orientation) than hole polygons, it lies +in the union. +""" +function union2d(p1, p2) + return clip( + Clipper.ClipTypeUnion, + p1, + p2, + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive + ) +end + +""" + union2d(p) + +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}[]) + +""" + difference2d(p1, p2) + +Return the geometric union of `p1` minus the geometric union of `p2` as a `ClippedPolygon`. + +Each of `p1` and `p2` may be a `GeometryEntity` or array of `GeometryEntity`. All entities +are first converted to polygons using [`to_polygons`](@ref). + +Each of `p1` and `p2` can also be a `GeometryStructure` or `GeometryReference`, in which case +`elements(flatten(p))` will be converted to polygons. + +Each can also be a pair `geom => layer`, where `geom` is a +`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection +of either, in which case only the elements in those layers will be used. +""" +function difference2d(plus, minus) + return clip( + Clipper.ClipTypeDifference, + plus, + minus, + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive + ) +end + +""" + intersect2d(p1, p2) + +Return the geometric union of `p1` intersected with the geometric union of `p2` as a `ClippedPolygon`. + +Each of `p1` and `p2` may be a `GeometryEntity` or array of `GeometryEntity`. All entities +are first converted to polygons using [`to_polygons`](@ref). + +Each of `p1` and `p2` can also be a `GeometryStructure` or `GeometryReference`, in which case +`elements(flatten(p))` will be converted to polygons. + +Each can also be a pair `geom => layer`, where `geom` is a +`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection +of either, in which case only the elements in those layers will be used. +""" +function intersect2d(plus, minus) + return clip( + Clipper.ClipTypeIntersection, + plus, + minus, + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive + ) +end + +""" + xor2d(p1, p2) + +Return the symmetric difference (XOR) of `p1` and `p2` as a `ClippedPolygon`. + +The XOR operation returns regions that are in either `p1` or `p2`, but not in both. +This is useful for finding non-overlapping regions between two sets of polygons. + +Each of `p1` and `p2` may be a `GeometryEntity` or array of `GeometryEntity`. All entities +are first converted to polygons using [`to_polygons`](@ref). + +Each of `p1` and `p2` can also be a `GeometryStructure` or `GeometryReference`, in which case +`elements(flatten(p))` will be converted to polygons. + +Each can also be a pair `geom => layer`, where `geom` is a +`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection +of either, in which case only the elements in those layers will be used. +""" +function xor2d(p1, p2) + return clip( + Clipper.ClipTypeXor, + p1, + p2, + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive + ) +end + +function add_path!( + c::Clipper.Clip, + path::Vector{Point{T}}, + polyType::Clipper.PolyType, + closed::Bool +) where {T <: Union{Int64, Unitful.Quantity{Int64}}} + return ccall( + (:add_path, libcclipper), + Cuchar, + (Ptr{Cvoid}, Ptr{Clipper.IntPoint}, Csize_t, Cint, Cuchar), + c.clipper_ptr, + path, + length(path), + Int(polyType), + closed + ) == 1 ? true : false +end + +# Clipping two identically-typed arrays of "Int64-based" Polygons. +# Internal method which should not be called by user (but does the heavy lifting) +function _clip( + op::Clipper.ClipType, + s::AbstractVector{Polygon{T}}, + c::AbstractVector{Polygon{T}}; + pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, + pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd +) where {T <: Union{Int64, Unitful.Quantity{Int64}}} + clip = clipper() + Clipper.clear!(clip) + for s0 in s + add_path!(clip, s0.p, Clipper.PolyTypeSubject, true) + end + for c0 in c + add_path!(clip, c0.p, Clipper.PolyTypeClip, true) + end + result = + convert(Clipper.PolyNode{Point{Int64}}, Clipper.execute_pt(clip, op, pfs, pfc)[2]) + + return ClippedPolygon(recast(Clipper.PolyNode{Point{T}}, result)) +end + +# recast(::Type{Clipper.PolyNode{T}}, x::Clipper.PolyNode}) where {T} +# Creates a `Clipper.PolyNode{T}` by reinterpreting vectors of points in `x`. +recast(::Type{Clipper.PolyNode{T}}, x::Clipper.PolyNode{T}) where {T} = x +function recast(::Type{Clipper.PolyNode{S}}, x::Clipper.PolyNode{T}) where {S, T} + pn = Clipper.PolyNode{S}( + reinterpret(S, Clipper.contour(x)), + Clipper.ishole(x), + Clipper.isopen(x) + ) + pn.children = [recast(y, pn) for y in Clipper.children(x)] + return pn.parent = pn +end +# recast(x::Clipper.PolyNode, parent::Clipper.PolyNode{S}) where {S} +# Creates a `Clipper.PolyNode{S}` from `x` given a new `parent` node. +function recast(x::Clipper.PolyNode, parent::Clipper.PolyNode{S}) where {S} + pn = Clipper.PolyNode{S}( + reinterpret(S, Clipper.contour(x)), + Clipper.ishole(x), + Clipper.isopen(x) + ) + pn.children = [recast(y, pn) for y in Clipper.children(x)] + pn.parent = parent + return pn +end + +# Int64like(x::Point{T}) where {T} +# Int64like(x::Polygon{T}) where {T} +# Converts Points or Polygons to an Int64-based representation (possibly with units). +Int64like(x::Point{T}) where {T} = convert(Point{typeof(Int64(1) * unit(T))}, x) +Int64like(x::Polygon{T}) where {T} = convert(Polygon{typeof(Int64(1) * unit(T))}, x) + +# prescale(x::Point{<:Real}) +# Since the Clipper library works on Int64-based points, we multiply floating-point-based +# `x` by `10.0^9` before rounding to retain high resolution. Since` 1.0` is interpreted +# to mean `1.0 um`, this yields `fm` resolution, which is more than sufficient for most uses. +prescale(x::Point{<:Real}) = x * SCALE # 2^29.897... + +# prescale(x::Point{<:Quantity}) +# Since the Clipper library works on Int64-based points, we unit-convert `x` to `fm` before +# rounding to retain high resolution, which is more than sufficient for most uses. +prescale(x::Point{<:Quantity}) = convert(Point{typeof(USCALE)}, x) + +# clipperize(A::AbstractVector{Polygon{T}}) where {T} +# clipperize(A::AbstractVector{Polygon{T}}) where {S<:Integer, T<:Union{S, Unitful.Quantity{S}}} +# clipperize(A::AbstractVector{Polygon{T}}) where {T <: Union{Int64, Unitful.Quantity{Int64}}} +# Prepare a vector of Polygons for being operated upon by the Clipper library, +# which expects Int64-based points (Quantity{Int64} is okay after using `reinterpret`). +function clipperize(A::AbstractVector{Polygon{T}}) where {T} + return [Polygon(clipperize.(points(x))) for x in A] +end + +# Already Integer-based, so no need to do rounding or scaling. Just convert to Int64-like. +function clipperize( + A::AbstractVector{Polygon{T}} +) where {S <: Integer, T <: Union{S, Unitful.Quantity{S}}} + return Int64like.(A) +end + +# Already Int64-based, so just pass through, nothing to do here. +function clipperize( + A::AbstractVector{Polygon{T}} +) where {T <: Union{Int64, Unitful.Quantity{Int64}}} + return A +end + +function clipperize(x::Point{T}) where {S <: Real, T <: Union{S, Unitful.Quantity{S}}} + return Int64like(unsafe_round(prescale(x))) +end +function clipperize( + x::Point{T} +) where {S <: Integer, D, U, T <: Union{S, Unitful.Quantity{S, D, U}}} + return Int64like(x) +end + +unscale(p::Point, ::Type{T}) where {T <: Quantity} = convert(Point{T}, p) +unscale(p::Point, ::Type{T}) where {T} = convert(Point{T}, p ./ SCALE) + +# Declipperize methods are used to get back to the original type. +declipperize(p, ::Type{T}) where {T} = Polygon{T}((x -> unscale(x, T)).(points(p))) +declipperize(p, ::Type{T}) where {T <: Union{Int64, Unitful.Quantity{Int64}}} = + Polygon{T}(reinterpret(Point{T}, points(p))) + +# Prepare a ClippedPolygon for use with Clipper. +function clipperize(p::ClippedPolygon) + R = typeof(clipperize(p.tree.children[1].contour[1])) + t = deepcopy(p.tree) + function prescale(p::Clipper.PolyNode) + Clipper.contour(p) .= (x -> unsafe_round(x * SCALE)).(Clipper.contour(p)) + for x in p.children + prescale(x) + end + end + prescale(t) + x = ClippedPolygon(convert(Clipper.PolyNode{R}, t)) + return x +end +function clipperize(p::ClippedPolygon{T}) where {T <: Quantity} + return ClippedPolygon(clipperize(p.tree)) +end + +# Prepare the data within a Clipper.PolyNode for use with Clipper. +function clipperize(p::Clipper.PolyNode) + # Create a tree by clipperizing contours recursively. + function buildtree(p) + T = typeof(clipperize.(p.contour)).parameters[1] + return Clipper.PolyNode{T}( + clipperize.(p.contour), + p.hole, + p.open, + buildtree.(p.children) + ) + end + t = buildtree(p) + + # Inform children of their heritage. + function labelchildren(node, parent) + for c ∈ node.children + c.parent = parent + labelchildren(c, node) + end + end + labelchildren(t, t) + return t +end + +# Convert a "clipperized" ClippedPolygon to a given type. +# Real valued clipping: convert the integer value back to float by dividing. +function declipperize(p::ClippedPolygon, ::Type{T}) where {T} + x = ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) + function unscale(p::Clipper.PolyNode) + Clipper.contour(p) .= (x -> x / SCALE).(Clipper.contour(p)) + for x in p.children + unscale(x) + end + end + unscale(x.tree) + return x +end +# Unitful quantities and integers use conversion directly. Extra methods resolve type +# ambiguities for aqua. +function declipperize( + p::ClippedPolygon{T}, + ::Type{T} +) where {T <: Union{Int, Quantity{Int}}} + return ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) +end +function declipperize(p::ClippedPolygon, ::Type{T}) where {T <: Union{Int, Quantity{Int}}} + return ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) +end +function declipperize(p::ClippedPolygon{<:Quantity}, ::Type{T}) where {T} + return ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) +end +function declipperize( + p::ClippedPolygon{<:Quantity}, + ::Type{T} +) where {T <: Union{Int, Quantity{Int}}} + return ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) +end + +""" + offset{S<:Coordinate}(s::AbstractPolygon{S}, delta::Coordinate; + j::Clipper.JoinType=Clipper.JoinTypeMiter, + e::Clipper.EndType=Clipper.EndTypeClosedPolygon) + offset{S<:AbstractPolygon}(subject::AbstractVector{S}, delta::Coordinate; + j::Clipper.JoinType=Clipper.JoinTypeMiter, + e::Clipper.EndType=Clipper.EndTypeClosedPolygon) + offset{S<:Polygon}(s::AbstractVector{S}, delta::Coordinate; + j::Clipper.JoinType=Clipper.JoinTypeMiter, + e::Clipper.EndType=Clipper.EndTypeClosedPolygon) + +Using the [`Clipper`](http://www.angusj.com/delphi/clipper.php) library and +the [`Clipper.jl`](https://github.com/Voxel8/Clipper.jl) wrapper, perform +polygon offsetting. + +The orientations of polygons must be consistent, such that outer polygons share the same +orientation, and any holes have the opposite orientation. Additionally, any holes should be +contained within outer polygons; offsetting hole edges may create positive artifacts at +corners. + +The first argument should be an [`AbstractPolygon`](@ref). The second argument +is how much to offset the polygon. Keyword arguments include a +[join type](http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Types/JoinType.htm): + + - `Clipper.JoinTypeMiter` + - `Clipper.JoinTypeRound` + - `Clipper.JoinTypeSquare` + +and also an +[end type](http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Types/EndType.htm): + + - `Clipper.EndTypeClosedPolygon` + - `Clipper.EndTypeClosedLine` + - `Clipper.EndTypeOpenSquare` + - `Clipper.EndTypeOpenRound` + - `Clipper.EndTypeOpenButt` +""" +function offset end + +function offset( + s::AbstractPolygon{T}, + delta::Coordinate; + j::Clipper.JoinType=Clipper.JoinTypeMiter, + e::Clipper.EndType=Clipper.EndTypeClosedPolygon +) where {T <: Coordinate} + dimension(T) != dimension(delta) && throw(Unitful.DimensionError(oneunit(T), delta)) + S = promote_type(T, typeof(delta)) + return offset(Polygon{S}[s], convert(S, delta); j=j, e=e) +end + +function offset( + s::AbstractVector{A}, + delta::Coordinate; + j::Clipper.JoinType=Clipper.JoinTypeMiter, + e::Clipper.EndType=Clipper.EndTypeClosedPolygon +) where {T, A <: AbstractPolygon{T}} + dimension(T) != dimension(delta) && throw(Unitful.DimensionError(oneunit(T), delta)) + S = promote_type(T, typeof(delta)) + + mask = typeof.(s) .<: ClippedPolygon + return offset( + convert( + Vector{Polygon{S}}, + [s[.!mask]..., reduce(vcat, to_polygons.(s[mask]); init=Polygon{S}[])...] + ), + convert(S, delta); + j=j, + e=e + ) +end + +prescaledelta(x::Real) = x * SCALE +prescaledelta(x::Integer) = x +prescaledelta(x::Length{<:Real}) = convert(typeof(USCALE), x) +prescaledelta(x::Length{<:Integer}) = x + +function offset( + s::AbstractVector{Polygon{T}}, + delta::T; + j::Clipper.JoinType=Clipper.JoinTypeMiter, + e::Clipper.EndType=Clipper.EndTypeClosedPolygon +) where {T <: Coordinate} + sc = clipperize(s) + d = prescaledelta(delta) + polys = _offset(sc, d, j=j, e=e) + return declipperize.(polys, T) +end + +function offset( + s::ClippedPolygon, + delta::T; + j::Clipper.JoinType=Clipper.JoinTypeMiter, + e::Clipper.EndType=Clipper.EndTypeClosedPolygon +) where {T <: Coordinate} + return offset(to_polygons(s), delta, j=j, e=e) +end + +function add_path!( + c::Clipper.ClipperOffset, + path::Vector{Point{T}}, + joinType::Clipper.JoinType, + endType::Clipper.EndType +) where {T <: Union{Int64, Unitful.Quantity{Int64}}} + return ccall( + (:add_offset_path, libcclipper), + Cvoid, + (Ptr{Cvoid}, Ptr{Clipper.IntPoint}, Csize_t, Cint, Cint), + c.clipper_ptr, + path, + length(path), + Int(joinType), + Int(endType) + ) +end + +function _offset( + s::AbstractVector{Polygon{T}}, + delta; + j::Clipper.JoinType=Clipper.JoinTypeMiter, + e::Clipper.EndType=Clipper.EndTypeClosedPolygon +) where {T <: Union{Int64, Unitful.Quantity{Int64}}} + c = coffset() + Clipper.clear!(c) + for s0 in s + add_path!(c, s0.p, j, e) + end + result = Clipper.execute(c, Float64(ustrip(delta))) #TODO: fix in clipper + return [Polygon(reinterpret(Point{T}, p)) for p in result] +end + +""" + orientation(p::Polygon) + +Return 1 if the points in the polygon contour are going counter-clockwise, -1 if clockwise. +Clipper considers clockwise-oriented polygons to be holes for some polygon fill types. +""" +function orientation(p::Polygon) + return ccall( + (:orientation, libcclipper), + Cuchar, + (Ptr{Clipper.IntPoint}, Csize_t), + reinterpret(Clipper.IntPoint, clipperize.(p.p)), + length(p.p) + ) == 1 ? 1 : -1 +end + +""" + ishole(p::Polygon) + +Return `true` if Clipper would consider this polygon to be a hole, for applicable +polygon fill rules. +""" +ishole(p::Polygon) = orientation(p) == -1 + +""" + orientation(p1::Point, p2::Point, p3::Point) + +Return 1 if the path `p1`--`p2`--`p3` is going counter-clockwise (increasing angle), +-1 if the path is going clockwise (decreasing angle), 0 if `p1`, `p2`, `p3` are colinear. +""" +function orientation(p1::Point, p2::Point, p3::Point) + return sign((p3.y - p2.y) * (p2.x - p1.x) - (p2.y - p1.y) * (p3.x - p2.x)) +end + +### cutting algorithm + +abstract type D1{T} end +Δy(d1::D1) = d1.p1.y - d1.p0.y +Δx(d1::D1) = d1.p1.x - d1.p0.x + +ab(p0, p1) = Point(gety(p1) - gety(p0), getx(p0) - getx(p1)) + +""" + LineSegment{T} <: D1{T} + +Represents a line segment. By construction, `p0.x <= p1.x`. +""" +struct LineSegment{T} <: D1{T} + p0::Point{T} + p1::Point{T} + function LineSegment(p0::Point{T}, p1::Point{T}) where {T} + if p1.x < p0.x + return new{T}(p1, p0) + else + return new{T}(p0, p1) + end + end +end +LineSegment(p0::Point{S}, p1::Point{T}) where {S, T} = LineSegment(promote(p0, p1)...) + +struct LineSegmentView{T} <: AbstractVector{T} + v::Vector{Point{T}} +end +Base.size(v::LineSegmentView) = size(v.v) +Base.length(v::LineSegmentView) = length(v.v) +Base.firstindex(v::LineSegmentView) = firstindex(v.v) +Base.lastindex(v::LineSegmentView) = lastindex(v.v) +function Base.getindex(v::LineSegmentView, i) + @boundscheck checkbounds(v.v, i) + return LineSegment(v.v[i], v.v[ifelse(i == length(v), 1, i + 1)]) +end + +""" + Ray{T} <: D1{T} + +Represents a ray. The ray starts at `p0` and goes toward `p1`. +""" +struct Ray{T} <: D1{T} + p0::Point{T} + p1::Point{T} +end +Ray(p0::Point{S}, p1::Point{T}) where {S, T} = Ray(promote(p0, p1)...) + +struct Line{T} <: D1{T} + p0::Point{T} + p1::Point{T} +end +Line(p0::Point{S}, p1::Point{T}) where {S, T} = Line(promote(p0, p1)...) +Line(seg::LineSegment) = Line(seg.p0, seg.p1) + +Base.promote_rule(::Type{Line{S}}, ::Type{Line{T}}) where {S, T} = Line{promote_type(S, T)} +Base.convert(::Type{Line{S}}, L::Line) where {S} = Line{S}(L.p0, L.p1) + +""" + segmentize(vertices, closed=true) + +Make an array of `LineSegment` out of an array of points. If `closed`, a segment should go +between the first and last point, otherwise nah. +""" +function segmentize(vertices, closed=true) + l = length(vertices) + if closed + return [LineSegment(vertices[i], vertices[i == l ? 1 : i + 1]) for i = 1:l] + else + return [LineSegment(vertices[i], vertices[i + 1]) for i = 1:(l - 1)] + end +end + +""" + uniqueray(v::Vector{Point{T}}) where {T <: Real} + +Given an array of points (thought to indicate a polygon or a hole in a polygon), +find the lowest / most negative y-coordinate[s] `miny`, then the lowest / most negative +x-coordinate `minx` of the points having that y-coordinate. This `Point(minx,miny)` ∈ `v`. +Return a ray pointing in -ŷ direction from that point. +""" +function uniqueray(v::Vector{Point{T}}) where {T <: Real} + nopts = reinterpret(T, v) + yarr = view(nopts, 2:2:length(nopts)) + miny, indy = findmin(yarr) + xarr = view(nopts, (findall(x -> x == miny, yarr) .* 2) .- 1) + minx, indx = findmin(xarr) + indv = findall(x -> x == Point(minx, miny), v)[1] + return Ray(Point(minx, miny), Point(minx, miny - 1)), indv +end + +isparallel(A::D1, B::D1) = Δy(A) * Δx(B) == Δy(B) * Δx(A) +isdegenerate(A::D1, B::D1) = + orientation(A.p0, A.p1, B.p0) == orientation(A.p0, A.p1, B.p1) == 0 +iscolinear(A::D1, B::Point) = orientation(A.p0, A.p1, B) == orientation(B, A.p1, A.p0) == 0 +iscolinear(A::Point, B::D1) = iscolinear(B, A) + +""" + intersects(A::LineSegment, B::LineSegment) + +Return two `Bool`s: + + 1. Does `A` intersect `B`? + 2. Did an intersection happen at a single point? (`false` if no intersection) +""" +function intersects(A::LineSegment, B::LineSegment) + sb0 = orientation(A.p0, A.p1, B.p0) + sb1 = orientation(A.p0, A.p1, B.p1) + sb = sb0 == sb1 + + sa0 = orientation(B.p0, B.p1, A.p0) + sa1 = orientation(B.p0, B.p1, A.p1) + sa = sa0 == sa1 + + if sa == false && sb == false + return true, true + else + # Test for special case of colinearity + if sb0 == sb1 == sa0 == sa1 == 0 + y0, y1 = minmax(A.p0.y, A.p1.y) + xinter = intersect(A.p0.x .. A.p1.x, B.p0.x .. B.p1.x) + yinter = intersect(A.p0.y .. A.p1.y, B.p0.y .. B.p1.y) + if !isempty(xinter) && !isempty(yinter) + if reduce(==, endpoints(xinter)) && reduce(==, endpoints(yinter)) + return true, true + else + return true, false + end + else + return false, false + end + else + return false, false + end + end +end + +""" + intersects_at_endpoint(A::LineSegment, B::LineSegment) + +Return three `Bool`s: + + 1. Does `A` intersect `B`? + 2. Did an intersection happen at a single point? (`false` if no intersection) + 3. Did an endpoint of `A` intersect an endpoint of `B`? +""" +function intersects_at_endpoint(A::LineSegment, B::LineSegment) + A_intersects_B, atapoint = intersects(A, B) + if A_intersects_B + if atapoint + if (A.p1 == B.p0) || (A.p1 == B.p1) || (A.p0 == B.p0) || (A.p0 == B.p1) + return A_intersects_B, atapoint, true + else + return A_intersects_B, atapoint, false + end + else + return A_intersects_B, atapoint, false + end + else + return A_intersects_B, atapoint, false + end +end + +""" + intersects(p::Point, A::Ray) + +Does `p` intersect `A`? +""" +function intersects(p::Point, A::Ray) + correctdir = sign(dot(A.p1 - A.p0, p - A.p0)) >= 0 + return iscolinear(p, A) && correctdir +end + +""" + in_bounds(p::Point, A::Ray) + +Is `p` in the halfspace defined by `A`? +""" +function in_bounds(p::Point, A::Ray) + return sign(dot(A.p1 - A.p0, p - A.p0)) >= 0 +end + +""" + intersects(p::Point, A::LineSegment) + +Does `p` intersect `A`? +""" +function intersects(p::Point, A::LineSegment) + if iscolinear(p, A) + y0, y1 = minmax(A.p0.y, A.p1.y) + xinter = intersect(A.p0.x .. A.p1.x, p.x .. p.x) + yinter = intersect(y0 .. y1, p.y .. p.y) + if !isempty(xinter) && !isempty(yinter) + return true + else + return false + end + else + return false + end +end + +""" + in_bounds(p::Point, A::LineSegment) + +Is `p` in the rectangle defined by the endpoints of `A`? +""" +function in_bounds(p::Point, A::LineSegment) + y0, y1 = minmax(A.p0.y, A.p1.y) + xinter = intersect(A.p0.x .. A.p1.x, p.x .. p.x) + yinter = intersect(y0 .. y1, p.y .. p.y) + return !isempty(xinter) && !isempty(yinter) +end + +function intersection(A::Ray{T}, B::LineSegment{T}) where {T} + fT = float(T) + if isparallel(A, B) + if isdegenerate(A, B) + # correct direction? + dist0 = dot(A.p1 - A.p0, B.p0 - A.p0) + dist1 = dot(A.p1 - A.p0, B.p1 - A.p0) + if sign(dist0) >= 0 + if sign(dist1) >= 0 + # Both in correct direction + return true, Point{fT}(min(dist0, dist1) == dist0 ? B.p0 : B.p1) + else + return true, Point{fT}(B.p0) + end + else + if sign(dist1) >= 0 + return true, Point{fT}(B.p1) + else + # Neither in correct direction + return false, zero(Point{fT}) + end + end + else + # no intersection + return false, zero(Point{fT}) + end + else + tf, w = intersection(Line(A.p0, A.p1), Line(B.p0, B.p1), false) + if tf && in_bounds(w, A) && in_bounds(w, B) + return true, w + else + return false, zero(Point{fT}) + end + end +end + +function intersection(A::Line{T}, B::Line{T}, checkparallel=true) where {T} + if checkparallel + # parallel checking goes here! + else + u = A.p1 - A.p0 + v = B.p1 - B.p0 + w = A.p0 - B.p0 + vp = Point{float(T)}(-v.y, v.x) # need float or hit overflow + + vp = vp / max(abs(vp.x), abs(vp.y)) # scale this, since its magnitude cancels out + # dot products will be smaller than maxintfloat(Float64) (assuming |w| and |u| are) + i = dot(-vp, w) / dot(vp, u) + return true, A.p0 + i * u + end +end + +mutable struct InteriorCutNode{T} + point::T + prev::InteriorCutNode{T} + next::InteriorCutNode{T} + + InteriorCutNode{T}(point, prev, next) where {T} = new{T}(point, prev, next) + function InteriorCutNode{T}(point) where {T} + node = new{T}(point) + node.prev = node + node.next = node + return node + end +end +segment(n::InteriorCutNode) = LineSegment(n.point, n.next.point) + +InteriorCutNode(val::T) where {T} = InteriorCutNode{T}(val) + +""" + interiorcuts(nodeortree::Clipper.PolyNode, outpolys::Vector{Polygon{T}}) where {T} + +Clipper gives polygons with holes as separate contours. The GDSII format doesn't support +this. This function makes cuts between the inner/outer contours so that ultimately there +is just one contour with one or more overlapping edges. + +Example: +┌────────────┐ ┌────────────┐ +│ ┌──┐ │ becomes... │ ┌──┐ │ +│ └──┘ ┌──┐ │ │ ├──┘ ┌──┐ │ +│ └──┘ │ │ │ ├──┘ │ +└────────────┘ └─┴─────┴────┘ +""" +function interiorcuts(nodeortree::Clipper.PolyNode, outpolys::Vector{Polygon{T}}) where {T} + # Assumes we have first element an enclosing polygon with the rest being holes. + # We also assume no hole collision. + + minpt = Point(-Inf, -Inf) + for enclosing in children(nodeortree) + enclosing_contour = contour(enclosing) + + # If a contour is empty, the PolyNode is effectively removed. This also effectively + # removes any further nodes, as they are no longer well defined. + isempty(enclosing_contour) && continue + + # No need to copy a large array of points, make a view giving line segments. + segs = LineSegmentView(enclosing_contour) + + # note to self: the problem has to do with segments reordering points... + + # Construct an interval tree of the x-extents of each line segment. + arr = reshape(reinterpret(Int, xinterval.(segs)), 2, :) + nodes = map(InteriorCutNode, enclosing_contour) + node1 = first(nodes) + for i in eachindex(nodes) + i == firstindex(nodes) || (nodes[i].prev = nodes[i - 1]) + i == lastindex(nodes) || (nodes[i].next = nodes[i + 1]) + end + IVT = IntervalValue{Int, InteriorCutNode{Point{Int}}} + iv = sort!(IVT.(view(arr, 1, :), view(arr, 2, :), nodes)) + itree = IntervalTree{Int, IVT}() + for v in iv + # We should be able to to bulk insertion, but it appears like this + # results in some broken trees for large enough initial insertion. + # see comments in merge request 21. + push!(itree, v) + end + loop_node = InteriorCutNode(enclosing_contour[1]) + loop_node.prev = last(nodes) + last(nodes).next = loop_node + + for hole in children(enclosing) + # process all the holes. + interiorcuts(hole, outpolys) + + # Intersect the unique ray with the line segments of the polygon. + hole_contour = contour(hole) + ray, m = uniqueray(hole_contour) + x0 = ray.p0.x + + # Find nearest intersection of the ray with the enclosing polygon. + best_intersection_point = minpt + local best_node + + # See which segments could possibly intersect with a line defined by `x = x0` + for interval in IntervalTrees.intersect(itree, (x0, x0)) + # Retrieve the segment index from the node. + node = IntervalTrees.value(interval) + seg = segment(node) + + # this is how we'll mark a "deleted" segment even though we don't + # actually remove it from the interval tree + (node.prev == node) && (node.next == node) && continue + + # See if it actually intersected with the segment + intersected, intersection_point = intersection(ray, seg) + if intersected + if gety(intersection_point) > gety(best_intersection_point) + best_intersection_point = intersection_point + best_node = node + end + end + end + + # Since the polygon was enclosing, an intersection had to happen *somewhere*. + if best_intersection_point != minpt + w = Point{Int64}( + round(getx(best_intersection_point)), + round(gety(best_intersection_point)) + ) + + # We are going to replace `best_node` + # need to do all of the following... + last_node = best_node.next + n0 = best_node.prev + + first_node = InteriorCutNode(best_node.point) + first_node.prev = n0 + n0.next = first_node + n0, p0 = first_node, w + + for r in (m:length(hole_contour), 1:m) + for i in r + n = InteriorCutNode(p0) + n.prev = n0 + n0.next = n + push!(itree, IntervalValue(xinterval(segment(n0))..., n0)) + n0, p0 = n, hole_contour[i] + end + end + + n = InteriorCutNode(p0) + n.prev = n0 + n0.next = n + push!(itree, IntervalValue(xinterval(segment(n0))..., n0)) + n0, p0 = n, w + + n = InteriorCutNode(p0) + n.prev = n0 + n0.next = n + push!(itree, IntervalValue(xinterval(segment(n0))..., n0)) + + n.next = last_node + last_node.prev = n + push!(itree, IntervalValue(xinterval(segment(n))..., n)) + + # serving the purpose of delete!(itree, best_node) + best_node.prev = best_node + best_node.next = best_node + + # in case we deleted node1... + if best_node === node1 + node1 = first_node + end + end + end + n = node1 + p = Point{Int}[] + while n.next != n + push!(p, n.point) + n = n.next + end + push!(outpolys, Polygon(reinterpret(Point{T}, p))) + end + return outpolys +end + +xinterval(l::LineSegment) = (l.p0.x, l.p1.x) +yinterval(l::LineSegment) = swap((l.p0.y, l.p1.y)) +swap(x) = x[1] > x[2] ? (x[2], x[1]) : x + +### Tiling +function SpatialIndexing.mbr(ent::GeometryEntity) + r = bounds(ent) + return SpatialIndexing.Rect((ustrip(r.ll)...,), (ustrip(r.ur)...,)) +end + +function spatial_index(ents::Vector{T}) where {T <: GeometryEntity} + tree = RTree{Float64, 2}(Int) + function convertel(enum_ent) + idx, ent = enum_ent + return SpatialIndexing.SpatialElem(mbr(ent), nothing, idx) + end + SpatialIndexing.load!(tree, enumerate(ents), convertel=convertel) + return tree +end + +# Goal is to have around 1000 polygons per tile +function intersect2d_tiled(poly1::Vector{Polygon{T}}, + poly2::Vector{Polygon{T}}, + max_tile_size=1mm; + heal=false) where T + # Create spatial index for each set of polygons + tree1 = spatial_index(poly1) + tree2 = spatial_index(poly2) + + # Get tiles and indices of polygons intersecting tiles + bnds = SpatialIndexing.combine(mbr(tree1), mbr(tree2)) + bnds_dl = Rectangle(Point(bnds.low...)*unit(T), Point(bnds.high...)*unit(T)) + tiles, edges = tiles_edges(bnds_dl, max_tile_size) # DeviceLayout Rectangles + tile_poly_indices = map(tiles) do tile + idx1 = intersecting_idx(tree1, tile) + idx2 = intersecting_idx(tree2, tile) + return (idx1, idx2) + end + # Intersect within each tile + output_poly = mapreduce(vcat, tile_poly_indices; init=Polygon{T}[]) do (idx1, idx2) + obj = @view poly1[idx1] + tool = @view poly2[idx2] + return to_polygons(intersect2d(obj, tool)) + end + + # Output from polygons touching edges may be duplicated + heal && heal_edges!(output_poly, edges) + # Output that does not itself touch an edge will still be duplicated + return output_poly +end + +function intersecting_idx(tree, tile) + return map(x -> x.val, intersects_with(tree, mbr(tile))) +end + +function heal_edges!(polygons::Vector{Polygon{T}}, edges) where T + tree = spatial_index(polygons) + touching_edge_idx = map(edges) do edge + intersects_with(tree, edge) + end + healed = reduce(vcat, + to_polygons(union2d(@view polygons[touching_edge_idx])), + init=Polygon{T}[]) + delete_at!(polygons, touching_edge_idx) + append!(polygons, healed) +end + +function tiles_edges(r::Rectangle, max_tile_size) + d = min(width(r), height(r)) + tile_size = d / ceil(d / max_tile_size) + nx = width(r) / tile_size + ny = 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 in 1:nx for j in 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:nx + ] + 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) +end + +import DeviceLayout: Cell, render!, aref, elements, flatten +using DeviceLayout.PreferredUnits +function benchmark_clip(ntot; tiled=false) + n = Int(round(sqrt(ntot))) + + circ1 = Cell("circ1", nm) + render!(circ1, Circle(10μm), GDSMeta(1)) + circ2 = Cell("circ2", nm) + render!(circ2, Circle(10μm), GDSMeta(2)) + + arr1 = aref(circ1, dc=Point(100μm, 0μm), dr=Point(0μm, 100μm), nc=n, nr=n) + arr2 = aref(circ2, Point(5μm, 5μm), dc=Point(100μm, 0μm), dr=Point(0μm, 100μm), nc=n, nr=n) + + poly1 = elements(flatten(arr1)) + poly2 = elements(flatten(arr2)) + + if tiled + @time "Tiled (total)" res = intersect2d_tiled(poly1, poly2) + else + @time "Direct (total)" res = to_polygons(intersect2d(poly1, poly2)) + end + return res +end \ No newline at end of file diff --git a/src/polygons.jl b/src/polygons.jl index 5e3b487f..714875bb 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 +using SpatialIndexing +import SpatialIndexing: mbr import DeviceLayout import DeviceLayout: @@ -2293,4 +2295,6 @@ function transform(d::StyleDict, f::Transformation) return newdict end +include("clipping.jl") + end # module From 36c7fb199d2d582db159fa23e359f7d91f095306 Mon Sep 17 00:00:00 2001 From: Greg Peairs Date: Fri, 6 Feb 2026 13:03:07 +0000 Subject: [PATCH 04/18] Add layerwise booleans Layerwise clipping Fix accidental reversion of 1D entity change Add tests with area, is_sliver Add more docstrings and CHANGELOG Add findbox and boolean interface tests Use exact tile size, wrap non-tiled result in vector Make sure mbr always uses float coordinates Fix bugged format Fix h_edges, remove unused code Re-format and add docstrings after rebase Add factor of 2 to sliver formula Move is_sliver definition to test setup --- CHANGELOG.md | 6 + docs/src/reference/api.md | 8 + src/DeviceLayout.jl | 18 +- src/clipping.jl | 311 ++++++++++++------ src/entities.jl | 50 +++ src/polygons.jl | 19 +- .../ReadoutResonators/clawed_meander.jl | 2 +- test/runtests.jl | 12 + test/test_bspline.jl | 7 +- test/test_clipping.jl | 62 ++++ test/test_render.jl | 5 + test/test_shapes.jl | 7 +- test/tests.jl | 5 + 13 files changed, 396 insertions(+), 116 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4cbe1ba1..f61fb81b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -64,8 +64,14 @@ There are also several minor features and fixes: - Added `SchematicDrivenLayout.filter_parameters` for sharing parameters between composite components and subcomponents - Added `rename_duplicates` option to `GDSWriterOptions` - Added experimental Text entity support to graphics backend + - Added layerwise Booleans `union2d_layerwise`, `difference2d_layerwise`, `intersect2d_layerwise`, and `xor2d_layerwise` + - Added `Polygons.area` + - Added `findbox(box, geoms; intersects=false)` for finding all elements of `geoms` whose bounding box is contained in or intersects `bounds(box)` + - Added `mbr_spatial_index` for creating an R-tree of minimum bounding rectangles + associated with an array of geometries, which can be used directly with `findbox` to avoid re-indexing for multiple `findbox` calls - Fixed bug where `map_metadata!` would map multiply-referenced structures multiple times - Fixed bug where `@composite_variant` would not forward `map_hooks` to base variant when defined with component instance rather than type + - Fixed overly-strict argument types for polygon clipping methods The documentation has also been reorganized and improved: diff --git a/docs/src/reference/api.md b/docs/src/reference/api.md index 8d74e06f..35f9aeb6 100644 --- a/docs/src/reference/api.md +++ b/docs/src/reference/api.md @@ -33,6 +33,8 @@ lowerleft(::DeviceLayout.GeometryEntity) upperright(::DeviceLayout.GeometryEntity) transform + findbox + mbr_spatial_index ``` ### [GeometryEntity](@id api-geometryentity) @@ -158,6 +160,7 @@ Polygon(::AbstractVector{Point{T}}) where {T} Polygon(::Point, ::Point, ::Point, ::Point...) Rectangle + Polygons.area bounds circle_polygon gridpoints_in_polygon @@ -174,10 +177,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..d53d4599 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") @@ -204,8 +207,10 @@ export GeometryEntity, center, centered, coordinatetype, + findbox, footprint, halo, + mbr_spatial_index, offset, to_polygons, lowerleft, @@ -378,10 +383,13 @@ import .Polygons: circle, circle_polygon, clip, + clip_tiled, cliptree, difference2d, + difference2d_layerwise, gridpoints_in_polygon, intersect2d, + intersect2d_layerwise, offset, perimeter, points, @@ -390,7 +398,9 @@ import .Polygons: sweep_poly, unfold, union2d, - xor2d + union2d_layerwise, + xor2d, + xor2d_layerwise export Polygons, Polygon, Ellipse, @@ -406,8 +416,10 @@ export Polygons, clip, cliptree, difference2d, + difference2d_layerwise, gridpoints_in_polygon, intersect2d, + intersect2d_layerwise, offset, perimeter, points, @@ -416,7 +428,9 @@ export Polygons, sweep_poly, unfold, union2d, - xor2d + union2d_layerwise, + xor2d, + xor2d_layerwise include("align.jl") using .Align diff --git a/src/clipping.jl b/src/clipping.jl index af2276cd..70019cd6 100644 --- a/src/clipping.jl +++ b/src/clipping.jl @@ -1,6 +1,6 @@ # Clipping polygons one at a time """ - clip(op::Clipper.ClipType, s, c; kwargs...) where {S<:Coordinate, T<:Coordinate} + clip(op::Clipper.ClipType, s, c; kwargs...) clip(op::Clipper.ClipType, s::AbstractVector{A}, c::AbstractVector{B}; kwargs...) where {S, T, A<:Polygon{S}, B<:Polygon{T}} clip(op::Clipper.ClipType, @@ -45,7 +45,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...) @@ -61,6 +61,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( @@ -86,15 +88,10 @@ function clip( s::AbstractVector{Polygon{T}}, c::AbstractVector{Polygon{T}}; pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, - pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, - tiled=false + pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd ) where {T} sc, cc = clipperize(s), clipperize(c) - polys = if !tiled - _clip(op, sc, cc; pfs, pfc) - else - _clip_tiled(op, sc, cc; pfs, pfc) - end + polys = _clip(op, sc, cc; pfs, pfc) return declipperize(polys, T) end @@ -216,7 +213,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) @@ -497,15 +495,9 @@ function declipperize( end """ - offset{S<:Coordinate}(s::AbstractPolygon{S}, delta::Coordinate; - j::Clipper.JoinType=Clipper.JoinTypeMiter, - e::Clipper.EndType=Clipper.EndTypeClosedPolygon) - offset{S<:AbstractPolygon}(subject::AbstractVector{S}, delta::Coordinate; - j::Clipper.JoinType=Clipper.JoinTypeMiter, - e::Clipper.EndType=Clipper.EndTypeClosedPolygon) - offset{S<:Polygon}(s::AbstractVector{S}, delta::Coordinate; - j::Clipper.JoinType=Clipper.JoinTypeMiter, - e::Clipper.EndType=Clipper.EndTypeClosedPolygon) + offset(s::AbstractPolygon{T}, delta::Coordinate; ...) where {T <: Coordinate} + offset(s::AbstractVector{A}, delta::Coordinate; ...) where {T, A <: AbstractPolygon{T}} + offset(s::AbstractVector{Polygon{T}}, delta::T; ...) where {T <: Coordinate} Using the [`Clipper`](http://www.angusj.com/delphi/clipper.php) library and the [`Clipper.jl`](https://github.com/Voxel8/Clipper.jl) wrapper, perform @@ -662,12 +654,16 @@ end ### cutting algorithm -abstract type D1{T} end +abstract type D1{T} <: GeometryEntity{T} end Δy(d1::D1) = d1.p1.y - d1.p0.y Δx(d1::D1) = d1.p1.x - d1.p0.x ab(p0, p1) = Point(gety(p1) - gety(p0), getx(p0) - getx(p1)) +to_polygons(::D1{T}) where {T} = Polygon{T}[] + +transform(d1::T, f::Transformation) where {T <: D1} = T(f(d1.p0), f(d1.p1)) + """ LineSegment{T} <: D1{T} @@ -1096,108 +1092,223 @@ xinterval(l::LineSegment) = (l.p0.x, l.p1.x) yinterval(l::LineSegment) = swap((l.p0.y, l.p1.y)) swap(x) = x[1] > x[2] ? (x[2], x[1]) : x -### Tiling -function SpatialIndexing.mbr(ent::GeometryEntity) - r = bounds(ent) - return SpatialIndexing.Rect((ustrip(r.ll)...,), (ustrip(r.ur)...,)) +### Layerwise +function xor2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...) + return clip_layerwise( + Clipper.ClipTypeXor, + obj, + tool; + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive, + kwargs... + ) end -function spatial_index(ents::Vector{T}) where {T <: GeometryEntity} - tree = RTree{Float64, 2}(Int) - function convertel(enum_ent) - idx, ent = enum_ent - return SpatialIndexing.SpatialElem(mbr(ent), nothing, idx) - end - SpatialIndexing.load!(tree, enumerate(ents), convertel=convertel) - return tree +function difference2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...) + return clip_layerwise( + Clipper.ClipTypeDifference, + obj, + tool; + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive, + kwargs... + ) end -# Goal is to have around 1000 polygons per tile -function intersect2d_tiled(poly1::Vector{Polygon{T}}, - poly2::Vector{Polygon{T}}, - max_tile_size=1mm; - heal=false) where T - # Create spatial index for each set of polygons - tree1 = spatial_index(poly1) - tree2 = spatial_index(poly2) - - # Get tiles and indices of polygons intersecting tiles - bnds = SpatialIndexing.combine(mbr(tree1), mbr(tree2)) - bnds_dl = Rectangle(Point(bnds.low...)*unit(T), Point(bnds.high...)*unit(T)) - tiles, edges = tiles_edges(bnds_dl, max_tile_size) # DeviceLayout Rectangles - tile_poly_indices = map(tiles) do tile - idx1 = intersecting_idx(tree1, tile) - idx2 = intersecting_idx(tree2, tile) - return (idx1, idx2) - end - # Intersect within each tile - output_poly = mapreduce(vcat, tile_poly_indices; init=Polygon{T}[]) do (idx1, idx2) - obj = @view poly1[idx1] - tool = @view poly2[idx2] - return to_polygons(intersect2d(obj, tool)) - end +function intersect2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...) + return clip_layerwise( + Clipper.ClipTypeIntersection, + obj, + tool; + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive, + kwargs... + ) +end - # Output from polygons touching edges may be duplicated - heal && heal_edges!(output_poly, edges) - # Output that does not itself touch an edge will still be duplicated - return output_poly +function union2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...) + return clip_layerwise( + Clipper.ClipTypeUnion, + obj, + tool; + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive, + kwargs... + ) end -function intersecting_idx(tree, tile) - return map(x -> x.val, intersects_with(tree, mbr(tile))) +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, + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive, + 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 + + Providing 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. For each tile, all entities with bounding box touching that tile + (including those touching edge-to-edge) are selected. + The values in the returned `Dict` are then lazy iterators over the results for each tile. + + Because entities touching more than one tile will be included in multiple operations, + resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal + is to find small differences between layouts: layouts are never incorrectly identified as identical, + and false positives are rare. In other cases, you may need to do some postprocessing or of the results. + + A rough guideline for choosing tile size is to aim for ~100 polygons per tile, but you may want to + benchmark your use case. + """ + eval(quote + @doc $doc $func_layerwise + end) end -function heal_edges!(polygons::Vector{Polygon{T}}, edges) where T - tree = spatial_index(polygons) - touching_edge_idx = map(edges) do edge - intersects_with(tree, edge) +function clip_layerwise( + op::Clipper.ClipType, + obj::GeometryStructure, + tool::GeometryStructure; + only_layers=[], + ignore_layers=[], + tile_size=nothing, + depth=-1, + pfs=Clipper.PolyFillTypeEvenOdd, + pfc=Clipper.PolyFillTypeEvenOdd +) + 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 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 - healed = reduce(vcat, - to_polygons(union2d(@view polygons[touching_edge_idx])), - init=Polygon{T}[]) - delete_at!(polygons, touching_edge_idx) - append!(polygons, healed) + return res end -function tiles_edges(r::Rectangle, max_tile_size) - d = min(width(r), height(r)) - tile_size = d / ceil(d / max_tile_size) - nx = width(r) / tile_size - ny = height(r) / tile_size +### 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 in 1:nx for j in 1:ny] + 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:nx + 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 + 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) + return tiles, vcat(h_edges, v_edges) # Edges could be used for healing end -import DeviceLayout: Cell, render!, aref, elements, flatten -using DeviceLayout.PreferredUnits -function benchmark_clip(ntot; tiled=false) - n = Int(round(sqrt(ntot))) +""" + function clip_tiled( + op, + ents1::AbstractArray{<:GeometryEntity{T}}, + ents2::AbstractArray{<:GeometryEntity{T}}, + tile_size=1000 * DeviceLayout.onemicron(T); + pfs=Clipper.PolyFillTypeEvenOdd, + pfc=Clipper.PolyFillTypeEvenOdd + ) - circ1 = Cell("circ1", nm) - render!(circ1, Circle(10μm), GDSMeta(1)) - circ2 = Cell("circ2", nm) - render!(circ2, Circle(10μm), GDSMeta(2)) +Return a lazy iterator that applies `op(ents1, ents2)` tile by tile. - arr1 = aref(circ1, dc=Point(100μm, 0μm), dr=Point(0μm, 100μm), nc=n, nr=n) - arr2 = aref(circ2, Point(5μm, 5μm), dc=Point(100μm, 0μm), dr=Point(0μm, 100μm), nc=n, nr=n) +The bounds of the combined geometries are tiled with squares with edge length `tile_size`, starting at the bottom left +corner. For each tile, all entities with bounding box touching that tile (including those touching edge-to-edge) are selected. +The return value is the lazy iterator over selected entities per tile. - poly1 = elements(flatten(arr1)) - poly2 = elements(flatten(arr2)) +Because entities touching more than one tile will be included in multiple operations, +resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal +is to find small differences between layouts: layouts are never incorrectly identified as identical, +and false positives are rare. In other cases, you may need to do some postprocessing or of the results. - if tiled - @time "Tiled (total)" res = intersect2d_tiled(poly1, poly2) - else - @time "Direct (total)" res = to_polygons(intersect2d(poly1, poly2)) +A rough guideline for choosing tile size is to aim for 100 polygons per tile, but you may want to +benchmark your use case. +""" +function clip_tiled( + op, + ents1::AbstractArray{<:GeometryEntity{T}}, + ents2::AbstractArray{<:GeometryEntity{T}}, + tile_size=1000 * DeviceLayout.onemicron(T); + pfs=Clipper.PolyFillTypeEvenOdd, + pfc=Clipper.PolyFillTypeEvenOdd +) where {T} + # Create spatial index for each set of polygons + tree1 = DeviceLayout.mbr_spatial_index(ents1) + tree2 = DeviceLayout.mbr_spatial_index(ents2) + + # Get tiles and indices of polygons intersecting tiles + bnds = SpatialIndexing.combine(mbr(tree1), mbr(tree2)) + bnds_dl = Rectangle( + Point(bnds.low...) * DeviceLayout.onemicron(T), + Point(bnds.high...) * DeviceLayout.onemicron(T) + ) + tiles, edges = tiles_and_edges(bnds_dl, tile_size) # DeviceLayout Rectangles + tile_poly_indices = map(tiles) do tile + idx1 = DeviceLayout.findbox(tile, tree1; intersects=true) + idx2 = DeviceLayout.findbox(tile, tree2; intersects=true) + return (idx1, idx2) + end + # Clip within each tile + res = Iterators.map(tile_poly_indices) do (idx1, idx2) + return clip(op, ents1[idx1], ents2[idx2]; pfs, pfc) end return res -end \ No newline at end of file +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 714875bb..674e9ce8 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -11,7 +11,7 @@ import CoordinateTransformations: AffineMap, LinearMap, Translation, Transformat import Clipper import Clipper: children, contour import StaticArrays -using SpatialIndexing +import SpatialIndexing import SpatialIndexing: mbr import DeviceLayout @@ -85,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 @@ -538,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°) 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/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..ae554ae9 100644 --- a/test/test_clipping.jl +++ b/test/test_clipping.jl @@ -479,6 +479,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 @@ -970,3 +987,48 @@ 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 findbox(r1, [r1, r2]) == [1] + @test findbox(r1, [r1, r2]; intersects=true) == [1, 2] + @test findbox(r1, [c1]) == [] + @test 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μm) + @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]))) +end 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 1f19187d..255512e5 100644 --- a/test/test_shapes.jl +++ b/test/test_shapes.jl @@ -624,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 From 555fd14675e6ce49cf112ac7861196298e7ee8f9 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Thu, 2 Apr 2026 12:29:54 +0200 Subject: [PATCH 05/18] Use PolyFillTypePositive as default in bare clip_layerwise and _tiled --- src/clipping.jl | 50 ++++++++++--------------------------------------- 1 file changed, 10 insertions(+), 40 deletions(-) diff --git a/src/clipping.jl b/src/clipping.jl index 70019cd6..4c76db8c 100644 --- a/src/clipping.jl +++ b/src/clipping.jl @@ -1094,47 +1094,19 @@ swap(x) = x[1] > x[2] ? (x[2], x[1]) : x ### Layerwise function xor2d_layerwise(obj::GeometryStructure, tool::GeometryStructure; kwargs...) - return clip_layerwise( - Clipper.ClipTypeXor, - obj, - tool; - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive, - 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; - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive, - 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; - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive, - 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; - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive, - kwargs... - ) + return clip_layerwise(Clipper.ClipTypeUnion, obj, tool; kwargs...) end for func in ("union2d", "difference2d", "intersect2d", "xor2d") @@ -1144,8 +1116,6 @@ for func in ("union2d", "difference2d", "intersect2d", "xor2d") only_layers=[], ignore_layers=[], depth=-1, - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive, tile_size=nothing ) @@ -1189,8 +1159,8 @@ function clip_layerwise( ignore_layers=[], tile_size=nothing, depth=-1, - pfs=Clipper.PolyFillTypeEvenOdd, - pfc=Clipper.PolyFillTypeEvenOdd + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive ) metadata_filter = DeviceLayout.layer_inclusion(only_layers, ignore_layers) if metadata_filter == DeviceLayout.trivial_inclusion @@ -1264,8 +1234,8 @@ end ents1::AbstractArray{<:GeometryEntity{T}}, ents2::AbstractArray{<:GeometryEntity{T}}, tile_size=1000 * DeviceLayout.onemicron(T); - pfs=Clipper.PolyFillTypeEvenOdd, - pfc=Clipper.PolyFillTypeEvenOdd + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive ) Return a lazy iterator that applies `op(ents1, ents2)` tile by tile. @@ -1287,8 +1257,8 @@ function clip_tiled( ents1::AbstractArray{<:GeometryEntity{T}}, ents2::AbstractArray{<:GeometryEntity{T}}, tile_size=1000 * DeviceLayout.onemicron(T); - pfs=Clipper.PolyFillTypeEvenOdd, - pfc=Clipper.PolyFillTypeEvenOdd + pfs=Clipper.PolyFillTypePositive, + pfc=Clipper.PolyFillTypePositive ) where {T} # Create spatial index for each set of polygons tree1 = DeviceLayout.mbr_spatial_index(ents1) From fa05043cb1baebb418d8191f336756ff2c007ccb Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Thu, 2 Apr 2026 12:58:26 +0200 Subject: [PATCH 06/18] Handle empty ents in clip_tiled --- src/clipping.jl | 21 ++++++++++++++++----- test/test_clipping.jl | 6 ++++++ 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/src/clipping.jl b/src/clipping.jl index 4c76db8c..f5ea64d1 100644 --- a/src/clipping.jl +++ b/src/clipping.jl @@ -1260,20 +1260,31 @@ function clip_tiled( 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 - tree1 = DeviceLayout.mbr_spatial_index(ents1) - tree2 = DeviceLayout.mbr_spatial_index(ents2) + if isempty(ents2) + tree1 = DeviceLayout.mbr_spatial_index(ents1) + bnds = mbr(tree1) + elseif isempty(ents2) + 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 = SpatialIndexing.combine(mbr(tree1), mbr(tree2)) bnds_dl = Rectangle( Point(bnds.low...) * DeviceLayout.onemicron(T), Point(bnds.high...) * DeviceLayout.onemicron(T) ) tiles, edges = tiles_and_edges(bnds_dl, tile_size) # DeviceLayout Rectangles tile_poly_indices = map(tiles) do tile - idx1 = DeviceLayout.findbox(tile, tree1; intersects=true) - idx2 = DeviceLayout.findbox(tile, tree2; intersects=true) + 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 diff --git a/test/test_clipping.jl b/test/test_clipping.jl index ae554ae9..1e4fee22 100644 --- a/test/test_clipping.jl +++ b/test/test_clipping.jl @@ -1031,4 +1031,10 @@ end @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 empty layers + uae = union2d_layerwise(ca_1, CoordinateSystem("empty")) + uae_tiled = union2d_layerwise(ca_1, CoordinateSystem("empty"), tile_size=99μm) + @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))) end From ad659fe13c0714d809539ad2cd4f22e5bdab661d Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Thu, 2 Apr 2026 13:20:32 +0200 Subject: [PATCH 07/18] Bump docs size threshold --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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", From b8cadfa63de7a49bf3f852a330896be40478f87a Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Fri, 3 Apr 2026 20:51:06 +0200 Subject: [PATCH 08/18] Move back to polygons.jl for a clean diff --- src/clipping.jl | 1295 ----------------------------------------------- src/polygons.jl | 209 +++++++- 2 files changed, 206 insertions(+), 1298 deletions(-) delete mode 100644 src/clipping.jl diff --git a/src/clipping.jl b/src/clipping.jl deleted file mode 100644 index f5ea64d1..00000000 --- a/src/clipping.jl +++ /dev/null @@ -1,1295 +0,0 @@ -# Clipping polygons one at a time -""" - clip(op::Clipper.ClipType, s, c; kwargs...) - clip(op::Clipper.ClipType, s::AbstractVector{A}, c::AbstractVector{B}; - 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} - -Return the `ClippedPolygon` resulting from a polygon clipping operation. - -Uses the [`Clipper`](http://www.angusj.com/delphi/clipper.php) library and the -[`Clipper.jl`](https://github.com/Voxel8/Clipper.jl) wrapper to perform polygon clipping. - -## Positional arguments - -The first argument must be one of the following types to specify a clipping operation: - - - `Clipper.ClipTypeDifference` - - `Clipper.ClipTypeIntersection` - - `Clipper.ClipTypeUnion` - - `Clipper.ClipTypeXor` - -Note that these are types; you should not follow them with `()`. - -The second and third argument may be a `GeometryEntity` or array of `GeometryEntity`. All entities -are first converted to polygons using [`to_polygons`](@ref). -Each can also be a `GeometryStructure` or `GeometryReference`, in which case -`elements(flatten(p))` will be converted to polygons. -Each can also be a pair `geom => layer`, where `geom` is a -`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection -of either, in which case only the elements in those layers will be taken from the flattened structure. - -## Keyword arguments - -`pfs` and `pfc` specify polygon fill rules for the `s` and `c` arguments, respectively. -These arguments may include: - - - `Clipper.PolyFillTypeNegative` - - `Clipper.PolyFillTypePositive` - - `Clipper.PolyFillTypeEvenOdd` - - `Clipper.PolyFillTypeNonZero` - -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), [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...) -end - -# Clipping requires an AbstractVector{Polygon{T}} -_normalize_clip_arg(p::Polygon) = [p] -_normalize_clip_arg(p::GeometryEntity) = _normalize_clip_arg(to_polygons(p)) -_normalize_clip_arg(p::AbstractArray{Polygon{T}}) where {T} = p -_normalize_clip_arg(p::AbstractArray{<:GeometryEntity{T}}) where {T} = - reduce(vcat, to_polygons.(p); init=Polygon{T}[]) -_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( - op::Clipper.ClipType, - s::AbstractVector{A}, - c::AbstractVector{B}; - kwargs... -) where {S, T, A <: Polygon{S}, B <: Polygon{T}} - dimension(S) != dimension(T) && throw(Unitful.DimensionError(oneunit(S), oneunit(T))) - R = promote_type(S, T) - - return clip( - op, - convert(Vector{Polygon{R}}, s), - convert(Vector{Polygon{R}}, c); - kwargs... - ) -end - -# Clipping two identically-typed arrays of <: Polygon -function clip( - op::Clipper.ClipType, - s::AbstractVector{Polygon{T}}, - c::AbstractVector{Polygon{T}}; - pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, - pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd -) where {T} - sc, cc = clipperize(s), clipperize(c) - polys = _clip(op, sc, cc; pfs, pfc) - return declipperize(polys, T) -end - -""" - cliptree(op::Clipper.ClipType, s::AbstractPolygon{S}, c::AbstractPolygon{T}; - kwargs...) where {S<:Coordinate, T<:Coordinate} - cliptree(op::Clipper.ClipType, s::AbstractVector{A}, c::AbstractVector{B}; - 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} - -Return a `Clipper.PolyNode` representing parent-child relationships between polygons and -interior holes. The units and number type may need to be converted. - -Uses the [`Clipper`](http://www.angusj.com/delphi/clipper.php) library and the -[`Clipper.jl`](https://github.com/Voxel8/Clipper.jl) wrapper to perform polygon clipping. - -## Positional arguments - -The first argument must be one of the following types to specify a clipping operation: - - - `Clipper.ClipTypeDifference` - - `Clipper.ClipTypeIntersection` - - `Clipper.ClipTypeUnion` - - `Clipper.ClipTypeXor` - -Note that these are types; you should not follow them with `()`. The second and third -arguments are `AbstractPolygon`s or vectors thereof. - -## Keyword arguments - -`pfs` and `pfc` specify polygon fill rules for the `s` and `c` arguments, respectively. -These arguments may include: - - - `Clipper.PolyFillTypeNegative` - - `Clipper.PolyFillTypePositive` - - `Clipper.PolyFillTypeEvenOdd` - - `Clipper.PolyFillTypeNonZero` - -See the [`Clipper` docs](http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Types/PolyFillType.htm) -for further information. -""" -function cliptree( - op::Clipper.ClipType, - s::AbstractPolygon{S}, - c::AbstractPolygon{T}; - kwargs... -) where {S <: Coordinate, T <: Coordinate} - dimension(S) != dimension(T) && throw(Unitful.DimensionError(oneunit(S), oneunit(T))) - R = promote_type(S, T) - return cliptree(op, Polygon{R}[s], Polygon{R}[c]; kwargs...)::Vector{Polygon{R}} -end - -function cliptree( - op::Clipper.ClipType, - s::AbstractVector{A}, - c::AbstractVector{B}; - kwargs... -) where {S, T, A <: AbstractPolygon{S}, B <: AbstractPolygon{T}} - dimension(S) != dimension(T) && throw(Unitful.DimensionError(oneunit(S), oneunit(T))) - R = promote_type(S, T) - return cliptree( - op, - convert(Vector{Polygon{R}}, s), - convert(Vector{Polygon{R}}, c); - kwargs... - )::Vector{Polygon{R}} -end - -function cliptree( - op::Clipper.ClipType, - s::AbstractVector{Polygon{T}}, - c::AbstractVector{Polygon{T}}; - pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, - pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd -) where {T} - sc, cc = clipperize(s), clipperize(c) - cpoly = _clip(op, sc, cc; pfs, pfc) - return declipperize(cpoly, T).tree -end - -""" - union2d(p1, p2) - -Return the geometric union of p1 and p2 as a `ClippedPolygon`. - -Each of `p1` and `p2` may be a `GeometryEntity` or array of `GeometryEntity`. All entities -are first converted to polygons using [`to_polygons`](@ref). - -Each of `p1` and `p2` can also be a `GeometryStructure` or `GeometryReference`, in which case -`elements(flatten(p))` will be converted to polygons. - -Each can also be a pair `geom => layer`, where `geom` is a -`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection -of either, in which case only the elements in those layers will used. - -This is not implemented as a method of `union` because you can have a set union of arrays of -polygons, which is a distinct operation. - -The Clipper polyfill rule is PolyFillTypePositive, meaning as long as a -region lies within more non-hole (by orientation) than hole polygons, it lies -in the union. -""" -function union2d(p1, p2) - return clip( - Clipper.ClipTypeUnion, - p1, - p2, - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive - ) -end - -""" - union2d(p) - -Return the geometric union of `p` or all entities in `p`. -""" -union2d(p::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) - -Return the geometric union of `p1` minus the geometric union of `p2` as a `ClippedPolygon`. - -Each of `p1` and `p2` may be a `GeometryEntity` or array of `GeometryEntity`. All entities -are first converted to polygons using [`to_polygons`](@ref). - -Each of `p1` and `p2` can also be a `GeometryStructure` or `GeometryReference`, in which case -`elements(flatten(p))` will be converted to polygons. - -Each can also be a pair `geom => layer`, where `geom` is a -`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection -of either, in which case only the elements in those layers will be used. -""" -function difference2d(plus, minus) - return clip( - Clipper.ClipTypeDifference, - plus, - minus, - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive - ) -end - -""" - intersect2d(p1, p2) - -Return the geometric union of `p1` intersected with the geometric union of `p2` as a `ClippedPolygon`. - -Each of `p1` and `p2` may be a `GeometryEntity` or array of `GeometryEntity`. All entities -are first converted to polygons using [`to_polygons`](@ref). - -Each of `p1` and `p2` can also be a `GeometryStructure` or `GeometryReference`, in which case -`elements(flatten(p))` will be converted to polygons. - -Each can also be a pair `geom => layer`, where `geom` is a -`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection -of either, in which case only the elements in those layers will be used. -""" -function intersect2d(plus, minus) - return clip( - Clipper.ClipTypeIntersection, - plus, - minus, - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive - ) -end - -""" - xor2d(p1, p2) - -Return the symmetric difference (XOR) of `p1` and `p2` as a `ClippedPolygon`. - -The XOR operation returns regions that are in either `p1` or `p2`, but not in both. -This is useful for finding non-overlapping regions between two sets of polygons. - -Each of `p1` and `p2` may be a `GeometryEntity` or array of `GeometryEntity`. All entities -are first converted to polygons using [`to_polygons`](@ref). - -Each of `p1` and `p2` can also be a `GeometryStructure` or `GeometryReference`, in which case -`elements(flatten(p))` will be converted to polygons. - -Each can also be a pair `geom => layer`, where `geom` is a -`GeometryStructure` or `GeometryReference`, while `layer` is a `DeviceLayout.Meta`, a layer name `Symbol`, and/or a collection -of either, in which case only the elements in those layers will be used. -""" -function xor2d(p1, p2) - return clip( - Clipper.ClipTypeXor, - p1, - p2, - pfs=Clipper.PolyFillTypePositive, - pfc=Clipper.PolyFillTypePositive - ) -end - -function add_path!( - c::Clipper.Clip, - path::Vector{Point{T}}, - polyType::Clipper.PolyType, - closed::Bool -) where {T <: Union{Int64, Unitful.Quantity{Int64}}} - return ccall( - (:add_path, libcclipper), - Cuchar, - (Ptr{Cvoid}, Ptr{Clipper.IntPoint}, Csize_t, Cint, Cuchar), - c.clipper_ptr, - path, - length(path), - Int(polyType), - closed - ) == 1 ? true : false -end - -# Clipping two identically-typed arrays of "Int64-based" Polygons. -# Internal method which should not be called by user (but does the heavy lifting) -function _clip( - op::Clipper.ClipType, - s::AbstractVector{Polygon{T}}, - c::AbstractVector{Polygon{T}}; - pfs::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd, - pfc::Clipper.PolyFillType=Clipper.PolyFillTypeEvenOdd -) where {T <: Union{Int64, Unitful.Quantity{Int64}}} - clip = clipper() - Clipper.clear!(clip) - for s0 in s - add_path!(clip, s0.p, Clipper.PolyTypeSubject, true) - end - for c0 in c - add_path!(clip, c0.p, Clipper.PolyTypeClip, true) - end - result = - convert(Clipper.PolyNode{Point{Int64}}, Clipper.execute_pt(clip, op, pfs, pfc)[2]) - - return ClippedPolygon(recast(Clipper.PolyNode{Point{T}}, result)) -end - -# recast(::Type{Clipper.PolyNode{T}}, x::Clipper.PolyNode}) where {T} -# Creates a `Clipper.PolyNode{T}` by reinterpreting vectors of points in `x`. -recast(::Type{Clipper.PolyNode{T}}, x::Clipper.PolyNode{T}) where {T} = x -function recast(::Type{Clipper.PolyNode{S}}, x::Clipper.PolyNode{T}) where {S, T} - pn = Clipper.PolyNode{S}( - reinterpret(S, Clipper.contour(x)), - Clipper.ishole(x), - Clipper.isopen(x) - ) - pn.children = [recast(y, pn) for y in Clipper.children(x)] - return pn.parent = pn -end -# recast(x::Clipper.PolyNode, parent::Clipper.PolyNode{S}) where {S} -# Creates a `Clipper.PolyNode{S}` from `x` given a new `parent` node. -function recast(x::Clipper.PolyNode, parent::Clipper.PolyNode{S}) where {S} - pn = Clipper.PolyNode{S}( - reinterpret(S, Clipper.contour(x)), - Clipper.ishole(x), - Clipper.isopen(x) - ) - pn.children = [recast(y, pn) for y in Clipper.children(x)] - pn.parent = parent - return pn -end - -# Int64like(x::Point{T}) where {T} -# Int64like(x::Polygon{T}) where {T} -# Converts Points or Polygons to an Int64-based representation (possibly with units). -Int64like(x::Point{T}) where {T} = convert(Point{typeof(Int64(1) * unit(T))}, x) -Int64like(x::Polygon{T}) where {T} = convert(Polygon{typeof(Int64(1) * unit(T))}, x) - -# prescale(x::Point{<:Real}) -# Since the Clipper library works on Int64-based points, we multiply floating-point-based -# `x` by `10.0^9` before rounding to retain high resolution. Since` 1.0` is interpreted -# to mean `1.0 um`, this yields `fm` resolution, which is more than sufficient for most uses. -prescale(x::Point{<:Real}) = x * SCALE # 2^29.897... - -# prescale(x::Point{<:Quantity}) -# Since the Clipper library works on Int64-based points, we unit-convert `x` to `fm` before -# rounding to retain high resolution, which is more than sufficient for most uses. -prescale(x::Point{<:Quantity}) = convert(Point{typeof(USCALE)}, x) - -# clipperize(A::AbstractVector{Polygon{T}}) where {T} -# clipperize(A::AbstractVector{Polygon{T}}) where {S<:Integer, T<:Union{S, Unitful.Quantity{S}}} -# clipperize(A::AbstractVector{Polygon{T}}) where {T <: Union{Int64, Unitful.Quantity{Int64}}} -# Prepare a vector of Polygons for being operated upon by the Clipper library, -# which expects Int64-based points (Quantity{Int64} is okay after using `reinterpret`). -function clipperize(A::AbstractVector{Polygon{T}}) where {T} - return [Polygon(clipperize.(points(x))) for x in A] -end - -# Already Integer-based, so no need to do rounding or scaling. Just convert to Int64-like. -function clipperize( - A::AbstractVector{Polygon{T}} -) where {S <: Integer, T <: Union{S, Unitful.Quantity{S}}} - return Int64like.(A) -end - -# Already Int64-based, so just pass through, nothing to do here. -function clipperize( - A::AbstractVector{Polygon{T}} -) where {T <: Union{Int64, Unitful.Quantity{Int64}}} - return A -end - -function clipperize(x::Point{T}) where {S <: Real, T <: Union{S, Unitful.Quantity{S}}} - return Int64like(unsafe_round(prescale(x))) -end -function clipperize( - x::Point{T} -) where {S <: Integer, D, U, T <: Union{S, Unitful.Quantity{S, D, U}}} - return Int64like(x) -end - -unscale(p::Point, ::Type{T}) where {T <: Quantity} = convert(Point{T}, p) -unscale(p::Point, ::Type{T}) where {T} = convert(Point{T}, p ./ SCALE) - -# Declipperize methods are used to get back to the original type. -declipperize(p, ::Type{T}) where {T} = Polygon{T}((x -> unscale(x, T)).(points(p))) -declipperize(p, ::Type{T}) where {T <: Union{Int64, Unitful.Quantity{Int64}}} = - Polygon{T}(reinterpret(Point{T}, points(p))) - -# Prepare a ClippedPolygon for use with Clipper. -function clipperize(p::ClippedPolygon) - R = typeof(clipperize(p.tree.children[1].contour[1])) - t = deepcopy(p.tree) - function prescale(p::Clipper.PolyNode) - Clipper.contour(p) .= (x -> unsafe_round(x * SCALE)).(Clipper.contour(p)) - for x in p.children - prescale(x) - end - end - prescale(t) - x = ClippedPolygon(convert(Clipper.PolyNode{R}, t)) - return x -end -function clipperize(p::ClippedPolygon{T}) where {T <: Quantity} - return ClippedPolygon(clipperize(p.tree)) -end - -# Prepare the data within a Clipper.PolyNode for use with Clipper. -function clipperize(p::Clipper.PolyNode) - # Create a tree by clipperizing contours recursively. - function buildtree(p) - T = typeof(clipperize.(p.contour)).parameters[1] - return Clipper.PolyNode{T}( - clipperize.(p.contour), - p.hole, - p.open, - buildtree.(p.children) - ) - end - t = buildtree(p) - - # Inform children of their heritage. - function labelchildren(node, parent) - for c ∈ node.children - c.parent = parent - labelchildren(c, node) - end - end - labelchildren(t, t) - return t -end - -# Convert a "clipperized" ClippedPolygon to a given type. -# Real valued clipping: convert the integer value back to float by dividing. -function declipperize(p::ClippedPolygon, ::Type{T}) where {T} - x = ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) - function unscale(p::Clipper.PolyNode) - Clipper.contour(p) .= (x -> x / SCALE).(Clipper.contour(p)) - for x in p.children - unscale(x) - end - end - unscale(x.tree) - return x -end -# Unitful quantities and integers use conversion directly. Extra methods resolve type -# ambiguities for aqua. -function declipperize( - p::ClippedPolygon{T}, - ::Type{T} -) where {T <: Union{Int, Quantity{Int}}} - return ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) -end -function declipperize(p::ClippedPolygon, ::Type{T}) where {T <: Union{Int, Quantity{Int}}} - return ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) -end -function declipperize(p::ClippedPolygon{<:Quantity}, ::Type{T}) where {T} - return ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) -end -function declipperize( - p::ClippedPolygon{<:Quantity}, - ::Type{T} -) where {T <: Union{Int, Quantity{Int}}} - return ClippedPolygon(convert(Clipper.PolyNode{Point{T}}, p.tree)) -end - -""" - offset(s::AbstractPolygon{T}, delta::Coordinate; ...) where {T <: Coordinate} - offset(s::AbstractVector{A}, delta::Coordinate; ...) where {T, A <: AbstractPolygon{T}} - offset(s::AbstractVector{Polygon{T}}, delta::T; ...) where {T <: Coordinate} - -Using the [`Clipper`](http://www.angusj.com/delphi/clipper.php) library and -the [`Clipper.jl`](https://github.com/Voxel8/Clipper.jl) wrapper, perform -polygon offsetting. - -The orientations of polygons must be consistent, such that outer polygons share the same -orientation, and any holes have the opposite orientation. Additionally, any holes should be -contained within outer polygons; offsetting hole edges may create positive artifacts at -corners. - -The first argument should be an [`AbstractPolygon`](@ref). The second argument -is how much to offset the polygon. Keyword arguments include a -[join type](http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Types/JoinType.htm): - - - `Clipper.JoinTypeMiter` - - `Clipper.JoinTypeRound` - - `Clipper.JoinTypeSquare` - -and also an -[end type](http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Types/EndType.htm): - - - `Clipper.EndTypeClosedPolygon` - - `Clipper.EndTypeClosedLine` - - `Clipper.EndTypeOpenSquare` - - `Clipper.EndTypeOpenRound` - - `Clipper.EndTypeOpenButt` -""" -function offset end - -function offset( - s::AbstractPolygon{T}, - delta::Coordinate; - j::Clipper.JoinType=Clipper.JoinTypeMiter, - e::Clipper.EndType=Clipper.EndTypeClosedPolygon -) where {T <: Coordinate} - dimension(T) != dimension(delta) && throw(Unitful.DimensionError(oneunit(T), delta)) - S = promote_type(T, typeof(delta)) - return offset(Polygon{S}[s], convert(S, delta); j=j, e=e) -end - -function offset( - s::AbstractVector{A}, - delta::Coordinate; - j::Clipper.JoinType=Clipper.JoinTypeMiter, - e::Clipper.EndType=Clipper.EndTypeClosedPolygon -) where {T, A <: AbstractPolygon{T}} - dimension(T) != dimension(delta) && throw(Unitful.DimensionError(oneunit(T), delta)) - S = promote_type(T, typeof(delta)) - - mask = typeof.(s) .<: ClippedPolygon - return offset( - convert( - Vector{Polygon{S}}, - [s[.!mask]..., reduce(vcat, to_polygons.(s[mask]); init=Polygon{S}[])...] - ), - convert(S, delta); - j=j, - e=e - ) -end - -prescaledelta(x::Real) = x * SCALE -prescaledelta(x::Integer) = x -prescaledelta(x::Length{<:Real}) = convert(typeof(USCALE), x) -prescaledelta(x::Length{<:Integer}) = x - -function offset( - s::AbstractVector{Polygon{T}}, - delta::T; - j::Clipper.JoinType=Clipper.JoinTypeMiter, - e::Clipper.EndType=Clipper.EndTypeClosedPolygon -) where {T <: Coordinate} - sc = clipperize(s) - d = prescaledelta(delta) - polys = _offset(sc, d, j=j, e=e) - return declipperize.(polys, T) -end - -function offset( - s::ClippedPolygon, - delta::T; - j::Clipper.JoinType=Clipper.JoinTypeMiter, - e::Clipper.EndType=Clipper.EndTypeClosedPolygon -) where {T <: Coordinate} - return offset(to_polygons(s), delta, j=j, e=e) -end - -function add_path!( - c::Clipper.ClipperOffset, - path::Vector{Point{T}}, - joinType::Clipper.JoinType, - endType::Clipper.EndType -) where {T <: Union{Int64, Unitful.Quantity{Int64}}} - return ccall( - (:add_offset_path, libcclipper), - Cvoid, - (Ptr{Cvoid}, Ptr{Clipper.IntPoint}, Csize_t, Cint, Cint), - c.clipper_ptr, - path, - length(path), - Int(joinType), - Int(endType) - ) -end - -function _offset( - s::AbstractVector{Polygon{T}}, - delta; - j::Clipper.JoinType=Clipper.JoinTypeMiter, - e::Clipper.EndType=Clipper.EndTypeClosedPolygon -) where {T <: Union{Int64, Unitful.Quantity{Int64}}} - c = coffset() - Clipper.clear!(c) - for s0 in s - add_path!(c, s0.p, j, e) - end - result = Clipper.execute(c, Float64(ustrip(delta))) #TODO: fix in clipper - return [Polygon(reinterpret(Point{T}, p)) for p in result] -end - -""" - orientation(p::Polygon) - -Return 1 if the points in the polygon contour are going counter-clockwise, -1 if clockwise. -Clipper considers clockwise-oriented polygons to be holes for some polygon fill types. -""" -function orientation(p::Polygon) - return ccall( - (:orientation, libcclipper), - Cuchar, - (Ptr{Clipper.IntPoint}, Csize_t), - reinterpret(Clipper.IntPoint, clipperize.(p.p)), - length(p.p) - ) == 1 ? 1 : -1 -end - -""" - ishole(p::Polygon) - -Return `true` if Clipper would consider this polygon to be a hole, for applicable -polygon fill rules. -""" -ishole(p::Polygon) = orientation(p) == -1 - -""" - orientation(p1::Point, p2::Point, p3::Point) - -Return 1 if the path `p1`--`p2`--`p3` is going counter-clockwise (increasing angle), --1 if the path is going clockwise (decreasing angle), 0 if `p1`, `p2`, `p3` are colinear. -""" -function orientation(p1::Point, p2::Point, p3::Point) - return sign((p3.y - p2.y) * (p2.x - p1.x) - (p2.y - p1.y) * (p3.x - p2.x)) -end - -### cutting algorithm - -abstract type D1{T} <: GeometryEntity{T} end -Δy(d1::D1) = d1.p1.y - d1.p0.y -Δx(d1::D1) = d1.p1.x - d1.p0.x - -ab(p0, p1) = Point(gety(p1) - gety(p0), getx(p0) - getx(p1)) - -to_polygons(::D1{T}) where {T} = Polygon{T}[] - -transform(d1::T, f::Transformation) where {T <: D1} = T(f(d1.p0), f(d1.p1)) - -""" - LineSegment{T} <: D1{T} - -Represents a line segment. By construction, `p0.x <= p1.x`. -""" -struct LineSegment{T} <: D1{T} - p0::Point{T} - p1::Point{T} - function LineSegment(p0::Point{T}, p1::Point{T}) where {T} - if p1.x < p0.x - return new{T}(p1, p0) - else - return new{T}(p0, p1) - end - end -end -LineSegment(p0::Point{S}, p1::Point{T}) where {S, T} = LineSegment(promote(p0, p1)...) - -struct LineSegmentView{T} <: AbstractVector{T} - v::Vector{Point{T}} -end -Base.size(v::LineSegmentView) = size(v.v) -Base.length(v::LineSegmentView) = length(v.v) -Base.firstindex(v::LineSegmentView) = firstindex(v.v) -Base.lastindex(v::LineSegmentView) = lastindex(v.v) -function Base.getindex(v::LineSegmentView, i) - @boundscheck checkbounds(v.v, i) - return LineSegment(v.v[i], v.v[ifelse(i == length(v), 1, i + 1)]) -end - -""" - Ray{T} <: D1{T} - -Represents a ray. The ray starts at `p0` and goes toward `p1`. -""" -struct Ray{T} <: D1{T} - p0::Point{T} - p1::Point{T} -end -Ray(p0::Point{S}, p1::Point{T}) where {S, T} = Ray(promote(p0, p1)...) - -struct Line{T} <: D1{T} - p0::Point{T} - p1::Point{T} -end -Line(p0::Point{S}, p1::Point{T}) where {S, T} = Line(promote(p0, p1)...) -Line(seg::LineSegment) = Line(seg.p0, seg.p1) - -Base.promote_rule(::Type{Line{S}}, ::Type{Line{T}}) where {S, T} = Line{promote_type(S, T)} -Base.convert(::Type{Line{S}}, L::Line) where {S} = Line{S}(L.p0, L.p1) - -""" - segmentize(vertices, closed=true) - -Make an array of `LineSegment` out of an array of points. If `closed`, a segment should go -between the first and last point, otherwise nah. -""" -function segmentize(vertices, closed=true) - l = length(vertices) - if closed - return [LineSegment(vertices[i], vertices[i == l ? 1 : i + 1]) for i = 1:l] - else - return [LineSegment(vertices[i], vertices[i + 1]) for i = 1:(l - 1)] - end -end - -""" - uniqueray(v::Vector{Point{T}}) where {T <: Real} - -Given an array of points (thought to indicate a polygon or a hole in a polygon), -find the lowest / most negative y-coordinate[s] `miny`, then the lowest / most negative -x-coordinate `minx` of the points having that y-coordinate. This `Point(minx,miny)` ∈ `v`. -Return a ray pointing in -ŷ direction from that point. -""" -function uniqueray(v::Vector{Point{T}}) where {T <: Real} - nopts = reinterpret(T, v) - yarr = view(nopts, 2:2:length(nopts)) - miny, indy = findmin(yarr) - xarr = view(nopts, (findall(x -> x == miny, yarr) .* 2) .- 1) - minx, indx = findmin(xarr) - indv = findall(x -> x == Point(minx, miny), v)[1] - return Ray(Point(minx, miny), Point(minx, miny - 1)), indv -end - -isparallel(A::D1, B::D1) = Δy(A) * Δx(B) == Δy(B) * Δx(A) -isdegenerate(A::D1, B::D1) = - orientation(A.p0, A.p1, B.p0) == orientation(A.p0, A.p1, B.p1) == 0 -iscolinear(A::D1, B::Point) = orientation(A.p0, A.p1, B) == orientation(B, A.p1, A.p0) == 0 -iscolinear(A::Point, B::D1) = iscolinear(B, A) - -""" - intersects(A::LineSegment, B::LineSegment) - -Return two `Bool`s: - - 1. Does `A` intersect `B`? - 2. Did an intersection happen at a single point? (`false` if no intersection) -""" -function intersects(A::LineSegment, B::LineSegment) - sb0 = orientation(A.p0, A.p1, B.p0) - sb1 = orientation(A.p0, A.p1, B.p1) - sb = sb0 == sb1 - - sa0 = orientation(B.p0, B.p1, A.p0) - sa1 = orientation(B.p0, B.p1, A.p1) - sa = sa0 == sa1 - - if sa == false && sb == false - return true, true - else - # Test for special case of colinearity - if sb0 == sb1 == sa0 == sa1 == 0 - y0, y1 = minmax(A.p0.y, A.p1.y) - xinter = intersect(A.p0.x .. A.p1.x, B.p0.x .. B.p1.x) - yinter = intersect(A.p0.y .. A.p1.y, B.p0.y .. B.p1.y) - if !isempty(xinter) && !isempty(yinter) - if reduce(==, endpoints(xinter)) && reduce(==, endpoints(yinter)) - return true, true - else - return true, false - end - else - return false, false - end - else - return false, false - end - end -end - -""" - intersects_at_endpoint(A::LineSegment, B::LineSegment) - -Return three `Bool`s: - - 1. Does `A` intersect `B`? - 2. Did an intersection happen at a single point? (`false` if no intersection) - 3. Did an endpoint of `A` intersect an endpoint of `B`? -""" -function intersects_at_endpoint(A::LineSegment, B::LineSegment) - A_intersects_B, atapoint = intersects(A, B) - if A_intersects_B - if atapoint - if (A.p1 == B.p0) || (A.p1 == B.p1) || (A.p0 == B.p0) || (A.p0 == B.p1) - return A_intersects_B, atapoint, true - else - return A_intersects_B, atapoint, false - end - else - return A_intersects_B, atapoint, false - end - else - return A_intersects_B, atapoint, false - end -end - -""" - intersects(p::Point, A::Ray) - -Does `p` intersect `A`? -""" -function intersects(p::Point, A::Ray) - correctdir = sign(dot(A.p1 - A.p0, p - A.p0)) >= 0 - return iscolinear(p, A) && correctdir -end - -""" - in_bounds(p::Point, A::Ray) - -Is `p` in the halfspace defined by `A`? -""" -function in_bounds(p::Point, A::Ray) - return sign(dot(A.p1 - A.p0, p - A.p0)) >= 0 -end - -""" - intersects(p::Point, A::LineSegment) - -Does `p` intersect `A`? -""" -function intersects(p::Point, A::LineSegment) - if iscolinear(p, A) - y0, y1 = minmax(A.p0.y, A.p1.y) - xinter = intersect(A.p0.x .. A.p1.x, p.x .. p.x) - yinter = intersect(y0 .. y1, p.y .. p.y) - if !isempty(xinter) && !isempty(yinter) - return true - else - return false - end - else - return false - end -end - -""" - in_bounds(p::Point, A::LineSegment) - -Is `p` in the rectangle defined by the endpoints of `A`? -""" -function in_bounds(p::Point, A::LineSegment) - y0, y1 = minmax(A.p0.y, A.p1.y) - xinter = intersect(A.p0.x .. A.p1.x, p.x .. p.x) - yinter = intersect(y0 .. y1, p.y .. p.y) - return !isempty(xinter) && !isempty(yinter) -end - -function intersection(A::Ray{T}, B::LineSegment{T}) where {T} - fT = float(T) - if isparallel(A, B) - if isdegenerate(A, B) - # correct direction? - dist0 = dot(A.p1 - A.p0, B.p0 - A.p0) - dist1 = dot(A.p1 - A.p0, B.p1 - A.p0) - if sign(dist0) >= 0 - if sign(dist1) >= 0 - # Both in correct direction - return true, Point{fT}(min(dist0, dist1) == dist0 ? B.p0 : B.p1) - else - return true, Point{fT}(B.p0) - end - else - if sign(dist1) >= 0 - return true, Point{fT}(B.p1) - else - # Neither in correct direction - return false, zero(Point{fT}) - end - end - else - # no intersection - return false, zero(Point{fT}) - end - else - tf, w = intersection(Line(A.p0, A.p1), Line(B.p0, B.p1), false) - if tf && in_bounds(w, A) && in_bounds(w, B) - return true, w - else - return false, zero(Point{fT}) - end - end -end - -function intersection(A::Line{T}, B::Line{T}, checkparallel=true) where {T} - if checkparallel - # parallel checking goes here! - else - u = A.p1 - A.p0 - v = B.p1 - B.p0 - w = A.p0 - B.p0 - vp = Point{float(T)}(-v.y, v.x) # need float or hit overflow - - vp = vp / max(abs(vp.x), abs(vp.y)) # scale this, since its magnitude cancels out - # dot products will be smaller than maxintfloat(Float64) (assuming |w| and |u| are) - i = dot(-vp, w) / dot(vp, u) - return true, A.p0 + i * u - end -end - -mutable struct InteriorCutNode{T} - point::T - prev::InteriorCutNode{T} - next::InteriorCutNode{T} - - InteriorCutNode{T}(point, prev, next) where {T} = new{T}(point, prev, next) - function InteriorCutNode{T}(point) where {T} - node = new{T}(point) - node.prev = node - node.next = node - return node - end -end -segment(n::InteriorCutNode) = LineSegment(n.point, n.next.point) - -InteriorCutNode(val::T) where {T} = InteriorCutNode{T}(val) - -""" - interiorcuts(nodeortree::Clipper.PolyNode, outpolys::Vector{Polygon{T}}) where {T} - -Clipper gives polygons with holes as separate contours. The GDSII format doesn't support -this. This function makes cuts between the inner/outer contours so that ultimately there -is just one contour with one or more overlapping edges. - -Example: -┌────────────┐ ┌────────────┐ -│ ┌──┐ │ becomes... │ ┌──┐ │ -│ └──┘ ┌──┐ │ │ ├──┘ ┌──┐ │ -│ └──┘ │ │ │ ├──┘ │ -└────────────┘ └─┴─────┴────┘ -""" -function interiorcuts(nodeortree::Clipper.PolyNode, outpolys::Vector{Polygon{T}}) where {T} - # Assumes we have first element an enclosing polygon with the rest being holes. - # We also assume no hole collision. - - minpt = Point(-Inf, -Inf) - for enclosing in children(nodeortree) - enclosing_contour = contour(enclosing) - - # If a contour is empty, the PolyNode is effectively removed. This also effectively - # removes any further nodes, as they are no longer well defined. - isempty(enclosing_contour) && continue - - # No need to copy a large array of points, make a view giving line segments. - segs = LineSegmentView(enclosing_contour) - - # note to self: the problem has to do with segments reordering points... - - # Construct an interval tree of the x-extents of each line segment. - arr = reshape(reinterpret(Int, xinterval.(segs)), 2, :) - nodes = map(InteriorCutNode, enclosing_contour) - node1 = first(nodes) - for i in eachindex(nodes) - i == firstindex(nodes) || (nodes[i].prev = nodes[i - 1]) - i == lastindex(nodes) || (nodes[i].next = nodes[i + 1]) - end - IVT = IntervalValue{Int, InteriorCutNode{Point{Int}}} - iv = sort!(IVT.(view(arr, 1, :), view(arr, 2, :), nodes)) - itree = IntervalTree{Int, IVT}() - for v in iv - # We should be able to to bulk insertion, but it appears like this - # results in some broken trees for large enough initial insertion. - # see comments in merge request 21. - push!(itree, v) - end - loop_node = InteriorCutNode(enclosing_contour[1]) - loop_node.prev = last(nodes) - last(nodes).next = loop_node - - for hole in children(enclosing) - # process all the holes. - interiorcuts(hole, outpolys) - - # Intersect the unique ray with the line segments of the polygon. - hole_contour = contour(hole) - ray, m = uniqueray(hole_contour) - x0 = ray.p0.x - - # Find nearest intersection of the ray with the enclosing polygon. - best_intersection_point = minpt - local best_node - - # See which segments could possibly intersect with a line defined by `x = x0` - for interval in IntervalTrees.intersect(itree, (x0, x0)) - # Retrieve the segment index from the node. - node = IntervalTrees.value(interval) - seg = segment(node) - - # this is how we'll mark a "deleted" segment even though we don't - # actually remove it from the interval tree - (node.prev == node) && (node.next == node) && continue - - # See if it actually intersected with the segment - intersected, intersection_point = intersection(ray, seg) - if intersected - if gety(intersection_point) > gety(best_intersection_point) - best_intersection_point = intersection_point - best_node = node - end - end - end - - # Since the polygon was enclosing, an intersection had to happen *somewhere*. - if best_intersection_point != minpt - w = Point{Int64}( - round(getx(best_intersection_point)), - round(gety(best_intersection_point)) - ) - - # We are going to replace `best_node` - # need to do all of the following... - last_node = best_node.next - n0 = best_node.prev - - first_node = InteriorCutNode(best_node.point) - first_node.prev = n0 - n0.next = first_node - n0, p0 = first_node, w - - for r in (m:length(hole_contour), 1:m) - for i in r - n = InteriorCutNode(p0) - n.prev = n0 - n0.next = n - push!(itree, IntervalValue(xinterval(segment(n0))..., n0)) - n0, p0 = n, hole_contour[i] - end - end - - n = InteriorCutNode(p0) - n.prev = n0 - n0.next = n - push!(itree, IntervalValue(xinterval(segment(n0))..., n0)) - n0, p0 = n, w - - n = InteriorCutNode(p0) - n.prev = n0 - n0.next = n - push!(itree, IntervalValue(xinterval(segment(n0))..., n0)) - - n.next = last_node - last_node.prev = n - push!(itree, IntervalValue(xinterval(segment(n))..., n)) - - # serving the purpose of delete!(itree, best_node) - best_node.prev = best_node - best_node.next = best_node - - # in case we deleted node1... - if best_node === node1 - node1 = first_node - end - end - end - n = node1 - p = Point{Int}[] - while n.next != n - push!(p, n.point) - n = n.next - end - push!(outpolys, Polygon(reinterpret(Point{T}, p))) - end - return outpolys -end - -xinterval(l::LineSegment) = (l.p0.x, l.p1.x) -yinterval(l::LineSegment) = swap((l.p0.y, l.p1.y)) -swap(x) = x[1] > x[2] ? (x[2], x[1]) : x - -### 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, - 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 - - Providing 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. For each tile, all entities with bounding box touching that tile - (including those touching edge-to-edge) are selected. - The values in the returned `Dict` are then lazy iterators over the results for each tile. - - Because entities touching more than one tile will be included in multiple operations, - resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal - is to find small differences between layouts: layouts are never incorrectly identified as identical, - and false positives are rare. In other cases, you may need to do some postprocessing or of the results. - - A rough guideline for choosing tile size is to aim for ~100 polygons per tile, but you may want to - benchmark your use case. - """ - eval(quote - @doc $doc $func_layerwise - end) -end - -function clip_layerwise( - op::Clipper.ClipType, - obj::GeometryStructure, - tool::GeometryStructure; - only_layers=[], - ignore_layers=[], - 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 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 clip_tiled( - op, - ents1::AbstractArray{<:GeometryEntity{T}}, - ents2::AbstractArray{<:GeometryEntity{T}}, - tile_size=1000 * DeviceLayout.onemicron(T); - 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. For each tile, all entities with bounding box touching that tile (including those touching edge-to-edge) are selected. -The return value is the lazy iterator over selected entities per tile. - -Because entities touching more than one tile will be included in multiple operations, -resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal -is to find small differences between layouts: layouts are never incorrectly identified as identical, -and false positives are rare. In other cases, you may need to do some postprocessing or of the results. - -A rough guideline for choosing tile size is to aim for 100 polygons per tile, but you may want to -benchmark your use case. -""" -function clip_tiled( - op, - ents1::AbstractArray{<:GeometryEntity{T}}, - ents2::AbstractArray{<:GeometryEntity{T}}, - tile_size=1000 * DeviceLayout.onemicron(T); - 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(ents2) - 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) - ) - tiles, edges = tiles_and_edges(bnds_dl, tile_size) # DeviceLayout Rectangles - 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(tile_poly_indices) do (idx1, idx2) - return clip(op, ents1[idx1], ents2[idx2]; pfs, pfc) - end - return res -end diff --git a/src/polygons.jl b/src/polygons.jl index 674e9ce8..c263e1ad 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -1001,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...) @@ -1019,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( @@ -1169,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) @@ -2312,6 +2315,206 @@ function transform(d::StyleDict, f::Transformation) return newdict end -include("clipping.jl") +### 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, + 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 + + Providing 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. For each tile, all entities with bounding box touching that tile + (including those touching edge-to-edge) are selected. + The values in the returned `Dict` are then lazy iterators over the results for each tile. + + Because entities touching more than one tile will be included in multiple operations, + resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal + is to find small differences between layouts: layouts are never incorrectly identified as identical, + and false positives are rare. In other cases, you may need to do some postprocessing or of the results. + + A rough guideline for choosing tile size is to aim for ~100 polygons per tile, but you may want to + benchmark your use case. + """ + eval(quote + @doc $doc $func_layerwise + end) +end + +function clip_layerwise( + op::Clipper.ClipType, + obj::GeometryStructure, + tool::GeometryStructure; + only_layers=[], + ignore_layers=[], + 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 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 clip_tiled( + op, + ents1::AbstractArray{<:GeometryEntity{T}}, + ents2::AbstractArray{<:GeometryEntity{T}}, + tile_size=1000 * DeviceLayout.onemicron(T); + 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. For each tile, all entities with bounding box touching that tile (including those touching edge-to-edge) are selected. +The return value is the lazy iterator over selected entities per tile. + +Because entities touching more than one tile will be included in multiple operations, +resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal +is to find small differences between layouts: layouts are never incorrectly identified as identical, +and false positives are rare. In other cases, you may need to do some postprocessing or of the results. + +A rough guideline for choosing tile size is to aim for 100 polygons per tile, but you may want to +benchmark your use case. +""" +function clip_tiled( + op, + ents1::AbstractArray{<:GeometryEntity{T}}, + ents2::AbstractArray{<:GeometryEntity{T}}, + tile_size=1000 * DeviceLayout.onemicron(T); + 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(ents2) + 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) + ) + tiles, edges = tiles_and_edges(bnds_dl, tile_size) # DeviceLayout Rectangles + 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(tile_poly_indices) do (idx1, idx2) + return clip(op, ents1[idx1], ents2[idx2]; pfs, pfc) + end + return res +end end # module From 1966d87a849f31d186fe0554a890713d6a8213e0 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Tue, 7 Apr 2026 15:13:56 +0200 Subject: [PATCH 09/18] Add auto tile size, make findbox/mbr private --- docs/src/reference/api.md | 2 -- src/DeviceLayout.jl | 2 -- src/polygons.jl | 36 +++++++++++++++++++++++++++--------- test/test_clipping.jl | 29 +++++++++++++++++++++++++---- 4 files changed, 52 insertions(+), 17 deletions(-) diff --git a/docs/src/reference/api.md b/docs/src/reference/api.md index 35f9aeb6..598eaf2f 100644 --- a/docs/src/reference/api.md +++ b/docs/src/reference/api.md @@ -33,8 +33,6 @@ lowerleft(::DeviceLayout.GeometryEntity) upperright(::DeviceLayout.GeometryEntity) transform - findbox - mbr_spatial_index ``` ### [GeometryEntity](@id api-geometryentity) diff --git a/src/DeviceLayout.jl b/src/DeviceLayout.jl index d53d4599..b693eccd 100644 --- a/src/DeviceLayout.jl +++ b/src/DeviceLayout.jl @@ -207,10 +207,8 @@ export GeometryEntity, center, centered, coordinatetype, - findbox, footprint, halo, - mbr_spatial_index, offset, to_polygons, lowerleft, diff --git a/src/polygons.jl b/src/polygons.jl index c263e1ad..d7328c67 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -2339,6 +2339,7 @@ for func in ("union2d", "difference2d", "intersect2d", "xor2d") only_layers=[], ignore_layers=[], depth=-1, + tiled=false, tile_size=nothing ) @@ -2353,11 +2354,16 @@ for func in ("union2d", "difference2d", "intersect2d", "xor2d") # Tiling - Providing a `tile_size` can significantly speed up operations and reduce maximum memory usage for large geometries. + 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. For each tile, all entities with bounding box touching that tile + 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. + + For each tile, all entities with bounding box touching that tile (including those touching edge-to-edge) are selected. The values in the returned `Dict` are then lazy iterators over the results for each tile. @@ -2365,9 +2371,6 @@ for func in ("union2d", "difference2d", "intersect2d", "xor2d") resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal is to find small differences between layouts: layouts are never incorrectly identified as identical, and false positives are rare. In other cases, you may need to do some postprocessing or of the results. - - A rough guideline for choosing tile size is to aim for ~100 polygons per tile, but you may want to - benchmark your use case. """ eval(quote @doc $doc $func_layerwise @@ -2380,6 +2383,7 @@ function clip_layerwise( tool::GeometryStructure; only_layers=[], ignore_layers=[], + tiled=false, tile_size=nothing, depth=-1, pfs=Clipper.PolyFillTypePositive, @@ -2394,7 +2398,7 @@ function clip_layerwise( obj_metas = unique(DeviceLayout.element_metadata(obj_flat)) tool_metas = unique(DeviceLayout.element_metadata(tool_flat)) all_metas = unique([obj_metas; tool_metas]) - if isnothing(tile_size) + if !tiled && isnothing(tile_size) res = Dict( meta => [ clip( @@ -2451,12 +2455,20 @@ function tiles_and_edges(r::Rectangle, tile_size) 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=1000 * DeviceLayout.onemicron(T); + tile_size=nothing; pfs=Clipper.PolyFillTypePositive, pfc=Clipper.PolyFillTypePositive ) @@ -2473,13 +2485,14 @@ is to find small differences between layouts: layouts are never incorrectly iden and false positives are rare. In other cases, you may need to do some postprocessing or of the results. A rough guideline for choosing tile size is to aim for 100 polygons per tile, but you may want to -benchmark your use case. +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=1000 * DeviceLayout.onemicron(T); + tile_size=nothing; pfs=Clipper.PolyFillTypePositive, pfc=Clipper.PolyFillTypePositive ) where {T} @@ -2503,6 +2516,11 @@ function clip_tiled( 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 tile_poly_indices = map(tiles) do tile idx1 = isempty(ents1) ? Int[] : DeviceLayout.findbox(tile, tree1; intersects=true) diff --git a/test/test_clipping.jl b/test/test_clipping.jl index 1e4fee22..aeba7102 100644 --- a/test/test_clipping.jl +++ b/test/test_clipping.jl @@ -1015,10 +1015,10 @@ end @test intersect2d_layerwise(c1, c2)[lyr_b][1] == overlap # Findbox - @test findbox(r1, [r1, r2]) == [1] - @test findbox(r1, [r1, r2]; intersects=true) == [1, 2] - @test findbox(r1, [c1]) == [] - @test findbox(r1, [r1, c1], intersects=true) == [1, 2] + @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") @@ -1037,4 +1037,25 @@ end uae_tiled = union2d_layerwise(ca_1, CoordinateSystem("empty"), tile_size=99μm) @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 From a92f02e4ddececbae26e8e7c1d99589ccbef47d6 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Tue, 7 Apr 2026 15:40:09 +0200 Subject: [PATCH 10/18] Fix empty first argument case --- src/polygons.jl | 2 +- test/test_clipping.jl | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/polygons.jl b/src/polygons.jl index d7328c67..37209d3f 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -2502,7 +2502,7 @@ function clip_tiled( if isempty(ents2) tree1 = DeviceLayout.mbr_spatial_index(ents1) bnds = mbr(tree1) - elseif isempty(ents2) + elseif isempty(ents1) tree2 = DeviceLayout.mbr_spatial_index(ents2) bnds = mbr(tree2) else diff --git a/test/test_clipping.jl b/test/test_clipping.jl index aeba7102..50b6cc77 100644 --- a/test/test_clipping.jl +++ b/test/test_clipping.jl @@ -1037,6 +1037,9 @@ end uae_tiled = union2d_layerwise(ca_1, CoordinateSystem("empty"), tile_size=99μm) @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μm) + @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") From c924393c50cf33dfa21d1414c481ab05fa322211 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Thu, 9 Apr 2026 12:50:28 +0200 Subject: [PATCH 11/18] Fix changelog after rebase --- CHANGELOG.md | 12 ++++++------ test/test_clipping.jl | 8 +------- 2 files changed, 7 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f61fb81b..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`) @@ -64,14 +70,8 @@ There are also several minor features and fixes: - Added `SchematicDrivenLayout.filter_parameters` for sharing parameters between composite components and subcomponents - Added `rename_duplicates` option to `GDSWriterOptions` - Added experimental Text entity support to graphics backend - - Added layerwise Booleans `union2d_layerwise`, `difference2d_layerwise`, `intersect2d_layerwise`, and `xor2d_layerwise` - - Added `Polygons.area` - - Added `findbox(box, geoms; intersects=false)` for finding all elements of `geoms` whose bounding box is contained in or intersects `bounds(box)` - - Added `mbr_spatial_index` for creating an R-tree of minimum bounding rectangles - associated with an array of geometries, which can be used directly with `findbox` to avoid re-indexing for multiple `findbox` calls - Fixed bug where `map_metadata!` would map multiply-referenced structures multiple times - Fixed bug where `@composite_variant` would not forward `map_hooks` to base variant when defined with component instance rather than type - - Fixed overly-strict argument types for polygon clipping methods The documentation has also been reorganized and improved: diff --git a/test/test_clipping.jl b/test/test_clipping.jl index 50b6cc77..ef749203 100644 --- a/test/test_clipping.jl +++ b/test/test_clipping.jl @@ -940,13 +940,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) From 0485e081b72c409765ae5489d915ca085df938a9 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Tue, 14 Apr 2026 15:55:28 +0200 Subject: [PATCH 12/18] Clip edge-touching polygons to tile --- src/polygons.jl | 24 ++++++++++++++++++++++-- test/test_clipping.jl | 13 +++++++++++++ 2 files changed, 35 insertions(+), 2 deletions(-) diff --git a/src/polygons.jl b/src/polygons.jl index 37209d3f..880b2872 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -2522,6 +2522,15 @@ function clip_tiled( 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 = mapreduce(vcat, edges, init=Int[]) do edge + return isempty(ents1) ? Int[] : DeviceLayout.findbox(edge, tree1; intersects=true) + end + # Same for ents2 + edge_touching_idx2 = 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 = @@ -2529,8 +2538,19 @@ function clip_tiled( return (idx1, idx2) end # Clip within each tile - res = Iterators.map(tile_poly_indices) do (idx1, idx2) - return clip(op, ents1[idx1], ents2[idx2]; pfs, pfc) + 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 diff --git a/test/test_clipping.jl b/test/test_clipping.jl index ef749203..ffe018b6 100644 --- a/test/test_clipping.jl +++ b/test/test_clipping.jl @@ -1026,6 +1026,19 @@ end @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μm) + @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μm) From e323952634b205b9fe9f57d93d16e43940f9b567 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Tue, 14 Apr 2026 15:59:04 +0200 Subject: [PATCH 13/18] Run formatter --- src/curvilinear.jl | 3 +-- src/polygons.jl | 6 ++++-- test/test_clipping.jl | 9 +++++++-- test/test_shapes.jl | 4 ++-- 4 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/curvilinear.jl b/src/curvilinear.jl index 2d3268bc..f8c9d0ba 100644 --- a/src/curvilinear.jl +++ b/src/curvilinear.jl @@ -470,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 -> [x, (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 diff --git a/src/polygons.jl b/src/polygons.jl index 880b2872..5651d659 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -2548,8 +2548,10 @@ function clip_tiled( 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])) + 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 diff --git a/test/test_clipping.jl b/test/test_clipping.jl index ffe018b6..8d8dec10 100644 --- a/test/test_clipping.jl +++ b/test/test_clipping.jl @@ -1032,8 +1032,13 @@ end @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) + 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 diff --git a/test/test_shapes.jl b/test/test_shapes.jl index 255512e5..75fd48d9 100644 --- a/test/test_shapes.jl +++ b/test/test_shapes.jl @@ -450,8 +450,8 @@ end # 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_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 From 1428d1f539944d81c34ab777904e50c4c22fb434 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Tue, 14 Apr 2026 17:38:58 +0200 Subject: [PATCH 14/18] Work around Unitful bug --- test/test_clipping.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_clipping.jl b/test/test_clipping.jl index 8d8dec10..86ca070a 100644 --- a/test/test_clipping.jl +++ b/test/test_clipping.jl @@ -1020,14 +1020,14 @@ end 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μm) + 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μm) + 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])...) @@ -1046,10 +1046,10 @@ end # Tiling with empty layers uae = union2d_layerwise(ca_1, CoordinateSystem("empty")) - uae_tiled = union2d_layerwise(ca_1, CoordinateSystem("empty"), tile_size=99μm) + 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μm) + 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))) From 8f604551ec4242bac0acfc052798c69e06fcd4e2 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Tue, 14 Apr 2026 19:18:38 +0200 Subject: [PATCH 15/18] Update docstrings for edge handling --- src/polygons.jl | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/src/polygons.jl b/src/polygons.jl index 5651d659..2c2702c6 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -2363,14 +2363,10 @@ for func in ("union2d", "difference2d", "intersect2d", "xor2d") 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. - For each tile, all entities with bounding box touching that tile - (including those touching edge-to-edge) are selected. - The values in the returned `Dict` are then lazy iterators over the results for each tile. - - Because entities touching more than one tile will be included in multiple operations, - resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal - is to find small differences between layouts: layouts are never incorrectly identified as identical, - and false positives are rare. In other cases, you may need to do some postprocessing or of the results. + 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 @@ -2476,13 +2472,9 @@ end 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. For each tile, all entities with bounding box touching that tile (including those touching edge-to-edge) are selected. -The return value is the lazy iterator over selected entities per tile. - -Because entities touching more than one tile will be included in multiple operations, -resulting polygons may be duplicated or incorrect. This is often acceptable, as in the case of `xor2d_layerwise` when the goal -is to find small differences between layouts: layouts are never incorrectly identified as identical, -and false positives are rare. In other cases, you may need to do some postprocessing or of the results. +corner. Entities crossing between tiles are split into their intersections with each tile before clipping. +For each tile, those intersection results and all entities and 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 From aefd43cf85b5298ca48091428dbca47f94758b8e Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Wed, 15 Apr 2026 12:13:09 +0200 Subject: [PATCH 16/18] Use Set for edge-touching-idx membership testing --- src/polygons.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/polygons.jl b/src/polygons.jl index 2c2702c6..5edf7f97 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -2473,7 +2473,7 @@ 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 and inside that tile are selected. +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 @@ -2515,13 +2515,17 @@ function clip_tiled( tiles, edges = tiles_and_edges(bnds_dl, tile_size) # DeviceLayout Rectangles # Get single vector of all entity indices touching any edge - edge_touching_idx1 = mapreduce(vcat, edges, init=Int[]) do edge - return isempty(ents1) ? Int[] : DeviceLayout.findbox(edge, tree1; intersects=true) - end + 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 = mapreduce(vcat, edges, init=Int[]) do edge - return isempty(ents2) ? Int[] : DeviceLayout.findbox(edge, tree2; intersects=true) - end + 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) From 64cace79ed6f199b5d6d015134056039f9f12a85 Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Wed, 15 Apr 2026 12:53:56 +0200 Subject: [PATCH 17/18] Change clip/cliptree default from EvenOdd to PolyFillTypePositive --- src/polygons.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/polygons.jl b/src/polygons.jl index 5edf7f97..d2ff5667 100644 --- a/src/polygons.jl +++ b/src/polygons.jl @@ -961,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. @@ -1045,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) @@ -1060,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. @@ -1125,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) @@ -1276,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) From 9d88fa4f8817015f05540b754ec621a602bedb2d Mon Sep 17 00:00:00 2001 From: Gregory Peairs Date: Wed, 15 Apr 2026 13:40:01 +0200 Subject: [PATCH 18/18] Fix clipping tests with clockwise polygons --- test/test_clipping.jl | 36 ++++++++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/test/test_clipping.jl b/test/test_clipping.jl index 86ca070a..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))