diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index f55a042f1..6d9a6dc84 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -131,6 +131,20 @@ include("methods/clipping/overlayng/noding/node_identity.jl") include("methods/clipping/overlayng/noding/split.jl") include("methods/clipping/overlayng/noding/emit.jl") +# OverlayNG engine core (phase 2a): the half-edge graph over the arrangement. +include("methods/clipping/overlayng/overlay_label.jl") +include("methods/clipping/overlayng/edge_source.jl") +include("methods/clipping/overlayng/half_edge.jl") +include("methods/clipping/overlayng/overlay_graph.jl") + +# OverlayNG engine core (phase 2b): labeller, result builders, and the driver. +include("methods/clipping/overlayng/overlay_labeller.jl") +include("methods/clipping/overlayng/maximal_edge_ring.jl") +include("methods/clipping/overlayng/polygon_builder.jl") +include("methods/clipping/overlayng/line_builder.jl") +include("methods/clipping/overlayng/intersection_point_builder.jl") +include("methods/clipping/overlayng/overlay_ng.jl") + include("methods/orientation.jl") include("methods/polygonize.jl") include("methods/minimum_bounding_circle.jl") diff --git a/src/methods/clipping/overlayng/edge_source.jl b/src/methods/clipping/overlayng/edge_source.jl new file mode 100644 index 000000000..4208ff257 --- /dev/null +++ b/src/methods/clipping/overlayng/edge_source.jl @@ -0,0 +1,59 @@ +# # EdgeSourceInfo — per-ring/line source topology (port of JTS `EdgeSourceInfo`) +# +# Phase 2a of the OverlayNG port (design doc §2.7, §3 amendment 3). Records the +# topological information carried by one source linework (a polygon ring or a +# line) from a single input geometry: which input it came from, its dimension, +# its ring role, and — for area rings — the signed `depth_delta`. +# +# One `EdgeSourceInfo` exists per `RelateSegmentString`, indexed by the segment +# string's position in the arrangement's `segstrings`. A `NodedEdge` reaches its +# source info through `NodedEdge.string_idx`; nothing is copied per noded edge +# (design §2.1 invariant 5), so nothing can desynchronize. +# +# Everything here is internal to GeometryOps — nothing is exported. + +""" + EdgeSourceInfo + +Source topology of one input linework (port of JTS `EdgeSourceInfo`): `index` +(`0` = A, `1` = B), `dim` (`DIM_A` for a ring, `DIM_L` for a line), `is_hole` +(ring role), and `depth_delta` (the signed side-labelling delta of an area ring, +`0` for a line). Consumed by the edge merger to build `OverlayLabel`s. +""" +struct EdgeSourceInfo + index :: Int8 + dim :: Int8 + is_hole :: Bool + depth_delta :: Int8 +end + +# The depth delta of one area ring, derived ONCE from the material-interior +# authority (design §2.7). JTS `EdgeNodingBuilder.computeDepthDelta` assigns the +# canonical delta `+1` to a ring in "Exterior-on-Left" orientation and flips to +# `-1` otherwise; equivalently (`Edge.locationLeft/Right`) a positive delta means +# Left = EXTERIOR, Right = INTERIOR — i.e. the material interior lies on the +# RIGHT. So `depth_delta = material_interior_on_left ? -1 : +1`, computed on the +# ring's stored order and matching that convention exactly. Folding the hole flip +# into `_ring_material_interior_on_left` keeps relate / overlay / extents in +# agreement by construction. +function _ring_depth_delta(m::Manifold, pts::Vector, is_hole::Bool; exact) + return _ring_material_interior_on_left(m, pts, is_hole; exact) ? Int8(-1) : Int8(1) +end + +# Build the source info for one segment string. Area rings carry a depth delta; +# lines carry none (`dim == DIM_L`, `depth_delta == 0`). +function _edge_source_info(m::Manifold, ss::RelateSegmentString; exact = True()) + index = ss.is_a ? Int8(0) : Int8(1) + if ss.dim == DIM_A + is_hole = _ss_is_hole(ss) + return EdgeSourceInfo(index, DIM_A, is_hole, _ring_depth_delta(m, ss.pts, is_hole; exact)) + end + #-- line (or any non-area edge): no side labelling + return EdgeSourceInfo(index, DIM_L, false, Int8(0)) +end + +# The per-string source-info table of an arrangement, aligned to `segstrings` +# (so `sources[NodedEdge.string_idx]` is the noded edge's source info). +function _edge_source_infos(m::Manifold, arr::NodedArrangement; exact = True()) + return [_edge_source_info(m, ss; exact) for ss in arr.segstrings] +end diff --git a/src/methods/clipping/overlayng/half_edge.jl b/src/methods/clipping/overlayng/half_edge.jl new file mode 100644 index 000000000..91e08d918 --- /dev/null +++ b/src/methods/clipping/overlayng/half_edge.jl @@ -0,0 +1,139 @@ +# # Half-edge substrate — the winged-edge algebra (port of JTS `edgegraph/HalfEdge`) +# +# Phase 2a of the OverlayNG port (design doc §3). A port of +# `edgegraph/HalfEdge.java`, adapted to integer-indexed storage: half-edges live +# in a `Vector` and refer to each other by 1-based index (`0` = null), not by +# object reference — the spike-validated representation ("indices, not object +# references, proved painless"). Directed half-edges occur in symmetric pairs; +# each stores the index of its symmetric partner (`sym`) and of the next +# half-edge CCW around its origin's star. +# +# Field contract (every element of the edge store must provide these mutable +# fields; `OverlayEdge` in overlay_graph.jl is the sole concrete implementation): +# - `origin :: Int32` — the origin node id (into the arrangement's node table). +# - `sym :: Int32` — index of the symmetric (oppositely-oriented) half-edge. +# - `o_next :: Int32` — index of the next half-edge CCW around the origin, with +# the same origin (JTS `oNext`); a degree-1 origin points +# it at the half-edge itself. +# - `dir_pt :: P` — the direction point for angular ordering: the parent +# segment's far endpoint, an original kernel vertex +# (design §3 amendment 1) — never a constructed coordinate. +# +# Adaptation note (winged-edge field choice): JTS stores `next` (the next edge CCW +# around the DESTINATION) and derives `oNext() == sym.next`. We store `o_next` +# (around the ORIGIN) directly, because it is the primitive every overlay +# traversal uses and it makes `oNext`/`degree`/`prev` index-only with no `sym` +# indirection; JTS's `next` is recoverable as `next(e) == o_next(sym(e))`. We +# bulk-build each node's star and sort it once with the exact angular comparator +# (the spike's `_order_stars!`), in place of JTS's incremental `insert`: the exact +# comparator is costly, so one `O(d log d)` sort per node beats `insert`'s repeated +# rescans, and integer storage makes a bulk sort trivial. +# +# Everything here is internal to GeometryOps — nothing is exported. + +# ## Core edge algebra (integer-indexed ports of the HalfEdge methods) + +# The symmetric (oppositely-oriented) partner (JTS `sym`). +@inline he_sym(edges, i::Integer) = @inbounds edges[i].sym + +# The next half-edge CCW around the origin, same origin (JTS `oNext`). +@inline he_onext(edges, i::Integer) = @inbounds edges[i].o_next + +# The origin node id. +@inline he_origin(edges, i::Integer) = @inbounds edges[i].origin + +# The destination node id: the origin of the sym edge (JTS `dest`). +@inline he_dest(edges, i::Integer) = @inbounds edges[he_sym(edges, i)].origin + +# The next half-edge CCW around the destination (JTS `next`), recovered from the +# stored `o_next`: `next(e) == oNext(sym(e))`. +@inline he_next(edges, i::Integer) = he_onext(edges, he_sym(edges, i)) + +# The previous half-edge CW around the origin, with that vertex as its +# destination. Always `he_next(he_prev(e)) == e`. Port of JTS `prev` (scan the +# origin star, return the last edge's sym). +function he_prev(edges, i::Integer) + curr = Int32(i); prev = Int32(i) + while true + prev = curr + curr = he_onext(edges, curr) + curr == i && break + end + return he_sym(edges, prev) +end + +# The degree of the origin vertex: the number of half-edges originating there +# (port of JTS `degree`). +function he_degree(edges, i::Integer) + d = 0 + e = Int32(i) + while true + d += 1 + e = he_onext(edges, e) + e == i && break + end + return d +end + +# The first node (degree != 2) reached walking `prev` from this edge, or `0` if +# the edge is part of a ring with no such node (port of JTS `prevNode`). +function he_prev_node(edges, i::Integer) + e = Int32(i) + while he_degree(edges, e) == 2 + e = he_prev(edges, e) + e == i && return Int32(0) + end + return e +end + +# The half-edge originating at this edge's origin with the given destination node +# id, or `0` if none (port of JTS `find`, on node ids rather than coordinates). +function he_find(edges, i::Integer, dest::Integer) + o = Int32(i) + while true + he_dest(edges, o) == dest && return o + o = he_onext(edges, o) + o == i && break + end + return Int32(0) +end + +# ## Symmetric-pair linkage and star ordering + +# Link a freshly created symmetric pair `(i, j)`: set the sym pointers and the +# single-segment star (each edge's `o_next` is itself, so a degree-1 origin's +# `oNext` is the edge itself — JTS `link`, restated for the origin star). Both are +# overwritten for interior-star edges once their origin's star is ordered. +@inline function he_link!(edges, i::Integer, j::Integer) + @inbounds edges[i].sym = Int32(j) + @inbounds edges[j].sym = Int32(i) + @inbounds edges[i].o_next = Int32(i) + @inbounds edges[j].o_next = Int32(j) + return nothing +end + +# Angular comparison of two half-edges that share an origin, about that origin's +# symbolic apex (the node's `NodeKey`), via the exact kernel comparator +# (`rk_compare_edge_dir`, design §3 amendment 1). Foreign directions at a +# coincidence-merged apex take the kernel's exact-rational slow path — no special +# casing here. +@inline function he_compare_angular(m::Manifold, edges, keys, i::Integer, j::Integer; exact) + apex = @inbounds keys[he_origin(edges, i)] + return rk_compare_edge_dir(m, apex, (@inbounds edges[i].dir_pt), (@inbounds edges[j].dir_pt); exact) +end + +# Order one node's star of outgoing half-edge indices CCW and wire the origin ring +# so that `he_onext` walks it (the spike's `_order_stars!`). `star` is mutated +# into CCW order. All members must share the origin whose apex is `keys[origin]`. +function he_order_star!(m::Manifold, edges, keys, star::Vector{Int32}; exact) + n = length(star) + if n == 1 + @inbounds edges[star[1]].o_next = star[1] + return nothing + end + sort!(star; lt = (i, j) -> he_compare_angular(m, edges, keys, i, j; exact) < 0) + @inbounds for t in 1:n + edges[star[t]].o_next = star[t == n ? 1 : t + 1] + end + return nothing +end diff --git a/src/methods/clipping/overlayng/intersection_point_builder.jl b/src/methods/clipping/overlayng/intersection_point_builder.jl new file mode 100644 index 000000000..97907e980 --- /dev/null +++ b/src/methods/clipping/overlayng/intersection_point_builder.jl @@ -0,0 +1,41 @@ +# # IntersectionPointBuilder — result points (port of JTS `IntersectionPointBuilder`) +# +# Phase 2b of the OverlayNG port (design doc §3). A port of +# `operation/overlayng/IntersectionPointBuilder.java`: extracts `Point` results +# from an Intersection of non-point inputs, at nodes incident on edges of BOTH +# inputs where none of the incident edges is itself in the result (i.e. the +# inputs meet at an isolated point). `isAllowCollapseLines` is always true here +# (non-strict), so the boundary-collapse guard in `isEdgeOf` is elided. +# +# Returns the result points as coordinate tuples; the driver wraps them. Nothing +# here is exported. + +function _build_points(g::OverlayGraph) + edges = g.edges + points = Tuple{Float64, Float64}[] + for ne in graph_node_edges(g) + if _is_result_point(edges, ne) + push!(points, node_point(g.arr, he_origin(edges, ne))) + end + end + return points +end + +# Port of `isResultPoint`: a node incident on edges of both inputs, none of them +# in the result. +function _is_result_point(edges, node_edge::Integer) + is_edge_of_a = false + is_edge_of_b = false + e = Int32(node_edge) + while true + oe_in_result(edges, e) && return false + lbl = oe_label(edges, e) + is_edge_of_a |= _is_edge_of(lbl, 0) + is_edge_of_b |= _is_edge_of(lbl, 1) + e = he_onext(edges, e) + e == node_edge && break + end + return is_edge_of_a && is_edge_of_b +end + +@inline _is_edge_of(lbl::OverlayLabel, i::Integer) = is_boundary(lbl, i) || is_line(lbl, i) diff --git a/src/methods/clipping/overlayng/line_builder.jl b/src/methods/clipping/overlayng/line_builder.jl new file mode 100644 index 000000000..e98cf5f67 --- /dev/null +++ b/src/methods/clipping/overlayng/line_builder.jl @@ -0,0 +1,77 @@ +# # LineBuilder — result lines from the overlay graph (port of JTS `LineBuilder`) +# +# Phase 2b of the OverlayNG port (design doc §3). A port of +# `operation/overlayng/LineBuilder.java`, restricted to the raw-edge extraction +# JTS actually uses (`addResultLines`); the merged path (`addResultLinesMerged`) +# is marked NOT USED in JTS and is skipped. Strict-mode branches are dropped — +# this engine runs the original (non-strict) JTS semantics, so +# `isAllowMixedResult` and `isAllowCollapseLines` are always true. +# +# Each `OverlayEdge` is a single segment, so every result line is one two-point +# `LineString` (the raw noded-edge output; JTS emits the same edges, just as +# possibly-longer chains). Everything here is internal — nothing is exported. + +# Build the result lines from the graph (port of `getLines`: `markResultLines` +# then `addResultLines`). +function _build_lines(m::Manifold, g::OverlayGraph, input, has_result_area::Bool, + op::_OverlayOpCode; exact) + edges = g.edges + area_index = _input_area_index(input) + #-- markResultLines + for i in eachindex(edges) + oe_in_result_either(edges, i) && continue + if _is_result_line(oe_label(edges, i), op, has_result_area, area_index) + oe_mark_in_result_line!(edges, i) + end + end + #-- addResultLines (raw noded edges) + lines = Vector{typeof(_edge_line(g, 1))}() + for i in eachindex(edges) + oe_in_result_line(edges, i) || continue + oe_is_visited(edges, i) && continue + push!(lines, _edge_line(g, i)) + oe_mark_visited_both!(edges, i) + end + return lines +end + +# The `LineString` of one result edge (a single segment, in the half-edge's +# direction). +_edge_line(g::OverlayGraph, i::Integer) = + GI.LineString([node_point(g.arr, he_origin(g.edges, i)), + node_point(g.arr, he_dest(g.edges, i))]) + +# Port of `isResultLine`. `is_allow_collapse_lines` / `is_allow_mixed_result` +# are always true here (non-strict), so those guards are elided. +function _is_result_line(lbl::OverlayLabel, op::_OverlayOpCode, has_result_area::Bool, + area_index::Integer) + #-- a boundary of a single geometry is only in the result as part of an area + is_boundary_singleton(lbl) && return false + #-- a collapse interior to its parent area (narrow gore / hole spike) + is_interior_collapse(lbl) && return false + + if op != OVERLAY_INTERSECTION + #-- collapsed edge in the other area's interior + is_collapse_and_not_part_interior(lbl) && return false + #-- a line inside the result area is subsumed by it + if has_result_area && area_index >= 0 && is_line_in_area(lbl, area_index) + return false + end + end + + #-- touching area boundaries produce a line for Intersection (mixed result) + op == OVERLAY_INTERSECTION && is_boundary_touch(lbl) && return true + + #-- otherwise, the op boolean logic over the effective line locations + a_loc = _effective_location(lbl, 0) + b_loc = _effective_location(lbl, 1) + return _is_result_of_op(op, a_loc, b_loc) +end + +# Port of `effectiveLocation`: line and collapse edges report INTERIOR so the op +# logic can include them where warranted. +@inline function _effective_location(lbl::OverlayLabel, gi::Integer) + is_collapse(lbl, gi) && return LOC_INTERIOR + is_line(lbl, gi) && return LOC_INTERIOR + return get_line_location(lbl, gi) +end diff --git a/src/methods/clipping/overlayng/maximal_edge_ring.jl b/src/methods/clipping/overlayng/maximal_edge_ring.jl new file mode 100644 index 000000000..c0db22e51 --- /dev/null +++ b/src/methods/clipping/overlayng/maximal_edge_ring.jl @@ -0,0 +1,299 @@ +# # Result ring linking — maximal rings, minimal-ring split, edge rings +# +# Phase 2b of the OverlayNG port (design doc §3). Ports two JTS files kept +# together because they form one pipeline over the marked result-area edges: +# - `MaximalEdgeRing.java` → maximal-ring linking + the minimal-ring split +# at self-touching nodes (`_MaxEdgeRing` + `_link_result_area_max_ring_at_node!`). +# - `OverlayEdgeRing.java` → one minimal result ring: its coordinates, its +# shell/hole role (via `_ring_is_ccw`), and hole containment (design §3 +# amendment 5 — never planar even-odd on emitted coordinates). +# +# Ring linkage lives in the phase-2a `OverlayEdge` handles: `next_result_max` +# links maximal rings, `next_result` links minimal rings, and `max_edge_ring` / +# `edge_ring` are integer handles (`0` = null) into a builder's ring vectors. +# Each `NodedEdge`/`OverlayEdge` is a single straight segment between two nodes, +# so a ring's coordinates are exactly the sequence of its edges' node points — +# no intermediate `addCoordinates` bookkeeping is needed. +# +# Everything here is internal to GeometryOps — nothing is exported. + +# A maximal edge ring: a cycle of result-area half-edges linked by +# `next_result_max`. Identified by `id` (stored in each member's `max_edge_ring` +# handle) so `_attach_max_edges!` can detect revisits (port of JTS +# `MaximalEdgeRing`). +mutable struct _MaxEdgeRing + id::Int32 + start_edge::Int32 +end + +# One minimal result ring (port of JTS `OverlayEdgeRing`): the emitted output +# coordinates (`ring_pts`, closed), the manifold kernel points backing the +# orientation and PIP predicates (`kernel_pts`; identical to `ring_pts` on the +# plane), its shell/hole role, a bounding box for extent pruning, its assigned +# shell / contained holes (handles, `0` = null), and a lazily-built indexed +# point-in-area locator over its own ring. +mutable struct _OverlayEdgeRing{P} + id::Int32 + start_edge::Int32 + ring_pts::Vector{Tuple{Float64, Float64}} + kernel_pts::Vector{P} + is_hole::Bool + bbox::NTuple{4, Float64} # (xmin, xmax, ymin, ymax) of ring_pts + shell::Int32 + holes::Vector{Int32} + locator::Any # Union{Nothing, IndexedPointInAreaLocator}, lazy +end + +# The polygon-builder working context (JTS `PolygonBuilder`'s mutable state). Held +# together so the ring-linking and placement functions share the graph edge store, +# the arrangement (for `node_point`), the manifold/exact predicate context, and the +# growing ring collections. Parameterized on `M`/`P`/`E` so `m`/`exact` stay +# concrete and the `Planar`/`Spherical` methods dispatch. +mutable struct _PolyBuilderCtx{M <: Manifold, P, E} + m::M + edges::Vector{OverlayEdge{P}} + arr::NodedArrangement{P} + exact::E + max_rings::Vector{_MaxEdgeRing} + edge_rings::Vector{_OverlayEdgeRing{P}} + shell_list::Vector{Int32} # handles into edge_rings + free_hole_list::Vector{Int32} +end + +@inline _ctx_point_type(::_PolyBuilderCtx{M, P}) where {M, P} = P + +# ## Maximal-ring linking at a node (port of `linkResultAreaMaxRingAtNode`) +# +# Design §3 amendment 4 (the known trap): this is called UNGATED for every +# in-result edge (JTS's `// TODO: skip already-linked` is deliberately +# unfulfilled — gating on already-linked loses degree-2 nodes whose lone +# out-edge was pre-linked as an in-edge). The per-node scan's own early return +# on an already-linked in-edge provides the necessary idempotency. +function _link_result_area_max_ring_at_node!(edges, node_edge::Integer) + #-- precondition: node_edge is in the result area + end_out = he_onext(edges, node_edge) + curr_out = end_out + #-- state machine: 1 = find an incoming result edge, 2 = link to an outgoing + state = 1 + curr_result_in = Int32(0) + while true + #-- if the found in-edge is already linked, this node is done + (curr_result_in != 0 && oe_is_result_max_linked(edges, curr_result_in)) && return nothing + if state == 1 + curr_in = he_sym(edges, curr_out) + if oe_in_result_area(edges, curr_in) + curr_result_in = curr_in + state = 2 + end + else # state == 2 + if oe_in_result_area(edges, curr_out) + oe_set_next_result_max!(edges, curr_result_in, curr_out) + state = 1 + end + end + curr_out = he_onext(edges, curr_out) + curr_out == end_out && break + end + state == 2 && throw(_OverlayTopologyError("no outgoing edge found")) + return nothing +end + +# Attach the edges of a maximal ring, tagging each with `mr.id` (port of the +# `MaximalEdgeRing` constructor / `attachEdges`). +function _attach_max_edges!(ctx, mr::_MaxEdgeRing) + edges = ctx.edges + edge = mr.start_edge + while true + edge == 0 && throw(_OverlayTopologyError("Ring edge is null")) + edges[edge].max_edge_ring == mr.id && + throw(_OverlayTopologyError("Ring edge visited twice at max-ring build")) + oe_next_result_max(edges, edge) == 0 && + throw(_OverlayTopologyError("Ring edge missing at max-ring build")) + edges[edge].max_edge_ring = mr.id + edge = oe_next_result_max(edges, edge) + edge == mr.start_edge && break + end + return nothing +end + +# ## Minimal-ring split (ports of `buildMinimalRings` / `linkMinimalRings`) +# +# Splits a self-touching maximal ring into OGC-valid minimal rings by relinking +# the max-ring edges in the OPPOSITE (CW) orientation via `next_result`. This is +# exactly the piece the spike prototypes faked (which produced invalid unions on +# many-island geometries). + +# Build the minimal rings of a maximal ring, returning their handles. +function _build_minimal_rings!(ctx, mr::_MaxEdgeRing) + _link_minimal_rings!(ctx, mr) + min_rings = Int32[] + edges = ctx.edges + e = mr.start_edge + while true + edges[e].edge_ring == 0 && push!(min_rings, _new_edge_ring!(ctx, e)) + e = oe_next_result_max(edges, e) + e == mr.start_edge && break + end + return min_rings +end + +function _link_minimal_rings!(ctx, mr::_MaxEdgeRing) + edges = ctx.edges + e = mr.start_edge + while true + _link_min_ring_edges_at_node!(ctx, e, mr) + e = oe_next_result_max(edges, e) + e == mr.start_edge && break + end + return nothing +end + +# Port of `linkMinRingEdgesAtNode`: relink this max ring's edges around one node +# into minimal rings (CW orientation, via `next_result`). +function _link_min_ring_edges_at_node!(ctx, node_edge::Integer, mr::_MaxEdgeRing) + edges = ctx.edges + end_out = Int32(node_edge) + curr_max_ring_out = Int32(node_edge) + curr_out = he_onext(edges, node_edge) + while true + _is_already_linked_min(edges, he_sym(edges, curr_out), mr) && return nothing + if curr_max_ring_out == 0 + curr_max_ring_out = _select_max_out_edge(edges, curr_out, mr) + else + curr_max_ring_out = _link_max_in_edge!(edges, curr_out, curr_max_ring_out, mr) + end + curr_out = he_onext(edges, curr_out) + curr_out == end_out && break + end + curr_max_ring_out != 0 && + throw(_OverlayTopologyError("Unmatched edge found during min-ring linking")) + return nothing +end + +@inline _is_already_linked_min(edges, edge::Integer, mr::_MaxEdgeRing) = + edges[edge].max_edge_ring == mr.id && oe_is_result_linked(edges, edge) + +@inline _select_max_out_edge(edges, curr_out::Integer, mr::_MaxEdgeRing) = + edges[curr_out].max_edge_ring == mr.id ? Int32(curr_out) : Int32(0) + +@inline function _link_max_in_edge!(edges, curr_out::Integer, curr_max_ring_out::Integer, + mr::_MaxEdgeRing) + curr_in = he_sym(edges, curr_out) + edges[curr_in].max_edge_ring != mr.id && return Int32(curr_max_ring_out) + oe_set_next_result!(edges, curr_in, curr_max_ring_out) + return Int32(0) +end + +# ## Minimal ring construction (port of the `OverlayEdgeRing` constructor) + +function _new_edge_ring!(ctx, start::Integer) + P = _ctx_point_type(ctx) + id = Int32(length(ctx.edge_rings) + 1) + ring = _OverlayEdgeRing{P}(id, Int32(start), Tuple{Float64, Float64}[], P[], + false, (0.0, 0.0, 0.0, 0.0), Int32(0), Int32[], nothing) + push!(ctx.edge_rings, ring) + _compute_ring!(ctx, ring) + return id +end + +# Port of `computeRingPts` + `computeRing`: walk the minimal ring via +# `next_result`, collecting node points; then derive the shell/hole role and the +# bounding box. +function _compute_ring!(ctx, ring::_OverlayEdgeRing) + edges = ctx.edges + pts = Tuple{Float64, Float64}[] + push!(pts, node_point(ctx.arr, he_origin(edges, ring.start_edge))) + edge = ring.start_edge + while true + edges[edge].edge_ring == ring.id && + throw(_OverlayTopologyError("Edge visited twice during ring-building")) + push!(pts, node_point(ctx.arr, he_dest(edges, edge))) + edges[edge].edge_ring = ring.id + ne = oe_next_result(edges, edge) + ne == 0 && throw(_OverlayTopologyError("Found null edge in ring")) + edge = ne + edge == ring.start_edge && break + end + #-- the last dest is the start origin, so pts is already closed; be defensive + pts[end] == pts[1] || push!(pts, pts[1]) + + ring.ring_pts = pts + ring.kernel_pts = _ring_kernel_pts(ctx.m, pts) + ring.is_hole = _ring_is_ccw(ctx.m, ring.kernel_pts; exact = ctx.exact) + ring.bbox = _ring_bbox(pts) + return nothing +end + +# Emitted output coordinates ARE the planar kernel points; the sphere converts +# each realized (lon, lat) back to a unit vector for the kernel predicates. +_ring_kernel_pts(::Planar, pts::Vector{Tuple{Float64, Float64}}) = pts +_ring_kernel_pts(m::Spherical, pts::Vector{Tuple{Float64, Float64}}) = + [_to_kernel_point(m, p) for p in pts] + +function _ring_bbox(pts::Vector{Tuple{Float64, Float64}}) + xmin = xmax = pts[1][1] + ymin = ymax = pts[1][2] + for p in pts + xmin = min(xmin, p[1]); xmax = max(xmax, p[1]) + ymin = min(ymin, p[2]); ymax = max(ymax, p[2]) + end + return (xmin, xmax, ymin, ymax) +end + +# ## Hole containment (ports of `OverlayEdgeRing.locate` / `contains` / …) + +# Lazily builds an indexed point-in-area locator over this ring and locates `p` +# (design §3 amendment 5: robust ray crossing via `rk_orient`, over the ring's +# emitted coordinates — never naive even-odd). +function _ring_locate(ctx, ring::_OverlayEdgeRing, p) + if ring.locator === nothing + ring.locator = IndexedPointInAreaLocator(ctx.m, GI.Polygon([ring.ring_pts]); + exact = ctx.exact) + end + return locate(ring.locator, p) +end + +# Whether `shell` contains `hole` (port of `contains` + `isPointInOrOut`). On the +# plane a bounding-box reject prefilters before the point tests; on the sphere the +# lon/lat box is unreliable near the poles/antimeridian, so containment is decided +# purely by the point tests (free holes are rare, this path is cold). +# +# Adaptation to the non-self-noding substrate: JTS uses `containsProperly` here, +# because a hole touching its own shell would have been connected into one maximal +# ring (JTS nodes A against itself). This substrate does NOT self-node a single +# input (design §2.2), so such a hole surfaces as a disconnected free hole whose +# bbox touches its shell's; the prefilter must therefore be non-strict +# (`_bbox_contains`), letting the point-in-area test — the real decision — run. +_ring_contains(ctx::_PolyBuilderCtx{<:Planar}, shell::_OverlayEdgeRing, hole::_OverlayEdgeRing) = + _bbox_contains(shell.bbox, hole.bbox) && _is_point_in_or_out(ctx, shell, hole) +_ring_contains(ctx::_PolyBuilderCtx{<:Spherical}, shell::_OverlayEdgeRing, hole::_OverlayEdgeRing) = + _is_point_in_or_out(ctx, shell, hole) + +function _is_point_in_or_out(ctx, shell::_OverlayEdgeRing, hole::_OverlayEdgeRing) + for p in hole.ring_pts + loc = _ring_locate(ctx, shell, p) + loc == LOC_INTERIOR && return true + loc == LOC_EXTERIOR && return false + #-- LOC_BOUNDARY: inconclusive, keep checking + end + return false +end + +@inline _bbox_contains(outer::NTuple{4, Float64}, inner::NTuple{4, Float64}) = + outer[1] <= inner[1] && outer[2] >= inner[2] && outer[3] <= inner[3] && outer[4] >= inner[4] + +# Port of `findEdgeRingContaining`: the innermost (smallest-envelope) shell in +# `candidates` that contains `hole`, or `0`. +function _find_edge_ring_containing(ctx, hole::_OverlayEdgeRing, candidates) + min_containing = Int32(0) + for sh in candidates + shell = ctx.edge_rings[sh] + if _ring_contains(ctx, shell, hole) + if min_containing == 0 || + _bbox_contains(ctx.edge_rings[min_containing].bbox, shell.bbox) + min_containing = Int32(sh) + end + end + end + return min_containing +end diff --git a/src/methods/clipping/overlayng/overlay_graph.jl b/src/methods/clipping/overlayng/overlay_graph.jl new file mode 100644 index 000000000..84cfa7fc8 --- /dev/null +++ b/src/methods/clipping/overlayng/overlay_graph.jl @@ -0,0 +1,341 @@ +# # Overlay graph — edge merge, half-edge pairs, and the topology graph +# +# Phase 2a of the OverlayNG port (design doc §3). Consumes the phase-1 +# `NodedArrangement` and the per-string `EdgeSourceInfo` table and produces the +# `OverlayGraph`: for every distinct noded edge a symmetric pair of `OverlayEdge` +# half-edges sharing one `OverlayLabel`, with each node's star ordered CCW. +# +# Three JTS files are ported here, kept together because they form one pipeline: +# - `Edge.java` → `MergeEdge` + `merge!` (depth-delta summing, collapse +# detection, hole-role merge, label creation). +# - `EdgeMerger.java` → `_merge_noded_edges` (coincident edges grouped by the +# unordered node-id pair — design §3 amendment 2 — +# never by coordinates). +# - `OverlayEdge.java` / `OverlayGraph.java` → `OverlayEdge` + `OverlayGraph`. +# +# Everything here is internal to GeometryOps — nothing is exported. + +# ========================================================================== +# Edge (JTS `Edge`): the merge intermediate + label creation +# ========================================================================== +# +# Renamed `MergeEdge` to avoid collision with `OverlayEdge` and RelateNG's +# `RelateEdge` in the single GeometryOps module. It carries no coordinates: its +# geometry is the node pair `(node_lo, node_hi)` (canonical direction = the base +# contributor's parent traversal order) plus the base contributor's parent +# `(string_idx, seg_idx)`, from which the half-edge direction points (the parent +# segment's original endpoints) are read at graph-build time. Topology info is the +# JTS `Edge` a/b `dim`/`depth_delta`/`is_hole` fields. + +mutable struct MergeEdge + node_lo :: Int32 # canonical origin node id (base contributor's sub-edge start) + node_hi :: Int32 # canonical dest node id (base contributor's sub-edge end) + string_idx :: Int32 # base contributor's parent string (for direction points) + seg_idx :: Int32 # base contributor's parent segment + + a_dim :: Int8 + a_depth_delta :: Int32 # summed across merges — widened from the Int8 source delta + a_is_hole :: Bool + + b_dim :: Int8 + b_depth_delta :: Int32 + b_is_hole :: Bool +end + +# Build the base `MergeEdge` for a noded edge, carrying its single source's info +# on the matching input index (port of JTS `Edge(pts, info)` + `copyInfo`). +function _merge_edge(ne::NodedEdge, src::EdgeSourceInfo) + if src.index == 0 + return MergeEdge(ne.node_lo, ne.node_hi, ne.string_idx, ne.seg_idx, + src.dim, Int32(src.depth_delta), src.is_hole, + DIM_NOT_PART, Int32(0), false) + else + return MergeEdge(ne.node_lo, ne.node_hi, ne.string_idx, ne.seg_idx, + DIM_NOT_PART, Int32(0), false, + src.dim, Int32(src.depth_delta), src.is_hole) + end +end + +# Whether this edge is part of a shell for input `gi` (a non-hole boundary). +_me_is_shell(e::MergeEdge, gi::Integer) = gi == 0 ? (e.a_dim == DIM_A && !e.a_is_hole) : + (e.b_dim == DIM_A && !e.b_is_hole) + +# The merged hole role for input `gi`: a shell if any contributor is a shell +# (port of JTS `isHoleMerged`; `is_hole` is stored, hence the flip). +_me_is_hole_merged(gi::Integer, e1::MergeEdge, e2::MergeEdge) = + !(_me_is_shell(e1, gi) || _me_is_shell(e2, gi)) + +# Merge a coincident contributor `inc` into `base` (port of JTS `Edge.merge`). +# Hole status is updated first (it reads the pre-update dims); dims take the max; +# depth deltas sum with a direction flip. The flip is exact from node ids: the +# unordered pair is unique, so `inc` runs the same direction as `base` iff their +# `node_lo` agree (design §3 amendment 2) — no coordinate comparison. +function _merge!(base::MergeEdge, inc::MergeEdge) + base.a_is_hole = _me_is_hole_merged(0, base, inc) + base.b_is_hole = _me_is_hole_merged(1, base, inc) + inc.a_dim > base.a_dim && (base.a_dim = inc.a_dim) + inc.b_dim > base.b_dim && (base.b_dim = inc.b_dim) + flip = inc.node_lo == base.node_lo ? Int32(1) : Int32(-1) + base.a_depth_delta += flip * inc.a_depth_delta + base.b_depth_delta += flip * inc.b_depth_delta + return base +end + +# ## Label creation from a merged edge (port of JTS `Edge.createLabel` + statics) + +@inline _del_sign(d::Integer) = d > 0 ? 1 : (d < 0 ? -1 : 0) + +# Positive delta ⇒ Left = EXTERIOR, Right = INTERIOR (JTS `locationLeft/Right`). +@inline function _location_left(depth_delta::Integer) + s = _del_sign(depth_delta) + return s == 1 ? LOC_EXTERIOR : (s == -1 ? LOC_INTERIOR : LOC_NONE) +end +@inline function _location_right(depth_delta::Integer) + s = _del_sign(depth_delta) + return s == 1 ? LOC_INTERIOR : (s == -1 ? LOC_EXTERIOR : LOC_NONE) +end + +# The effective edge-state dimension of a source (port of JTS `Edge.labelDim`): +# an area edge with zero summed depth delta is a collapse. +@inline function _label_dim(dim::Int8, depth_delta::Integer) + dim == DIM_NOT_PART && return DIM_NOT_PART + dim == DIM_L && return DIM_L + #-- dim == DIM_A (area) + return depth_delta == 0 ? DIM_COLLAPSE : DIM_A +end + +# Initialize one source's slot of a label (port of JTS `Edge.initLabel`). +function _init_label!(l::OverlayLabel, gi::Integer, dim::Int8, depth_delta::Integer, is_hole::Bool) + dl = _label_dim(dim, depth_delta) + if dl == DIM_NOT_PART + init_not_part!(l, gi) + elseif dl == DIM_A # boundary + init_boundary!(l, gi, _location_left(depth_delta), _location_right(depth_delta), is_hole) + elseif dl == DIM_COLLAPSE + init_collapse!(l, gi, is_hole) + elseif dl == DIM_L # line + init_line!(l, gi) + end + return l +end + +# The shared `OverlayLabel` for a merged edge (port of JTS `Edge.createLabel`). +function _create_label(e::MergeEdge) + l = OverlayLabel() + _init_label!(l, 0, e.a_dim, e.a_depth_delta, e.a_is_hole) + _init_label!(l, 1, e.b_dim, e.b_depth_delta, e.b_is_hole) + return l +end + +# ## EdgeMerger (port of JTS `EdgeMerger.merge`, keyed by unordered node pair) + +@inline _edge_key(ne::NodedEdge) = + ne.node_lo < ne.node_hi ? (ne.node_lo, ne.node_hi) : (ne.node_hi, ne.node_lo) + +# Merge all noded edges of the arrangement. Coincident edges (same unordered node +# pair; the segment/minor-arc between two nodes is unique, antipodal edges +# excluded at ingest) collapse to one `MergeEdge`, the first seen setting the +# canonical direction and later ones merged into it. +function _merge_noded_edges(arr::NodedArrangement, sources::Vector{EdgeSourceInfo}) + edgemap = Dict{Tuple{Int32, Int32}, Int}() + merged = MergeEdge[] + for ne in arr.edges + src = sources[ne.string_idx] + key = _edge_key(ne) + idx = get(edgemap, key, 0) + if idx == 0 + push!(merged, _merge_edge(ne, src)) + edgemap[key] = length(merged) + else + _merge!(merged[idx], _merge_edge(ne, src)) + end + end + return merged +end + +# ========================================================================== +# OverlayEdge (JTS `OverlayEdge`): the graph half-edge +# ========================================================================== +# +# A directed half-edge in the graph. Carries the winged-edge fields required by +# half_edge.jl (`origin`, `sym`, `o_next`, `dir_pt`), the shared `label`, the +# `is_forward` direction (which reinterprets the label's Left/Right), the result- +# marking flags, and the ring-linkage pointers phase 2b fills in. Ring pointers +# are integer half-edge indices (`0` = null); the `edge_ring`/`max_edge_ring` +# fields are integer handles into phase 2b's ring collections (`0` = null). + +mutable struct OverlayEdge{P} + origin :: Int32 # origin node id + sym :: Int32 # index of the symmetric half-edge + o_next :: Int32 # next half-edge CCW around the origin (same origin) + dir_pt :: P # direction point: parent segment's far endpoint (original vertex) + + is_forward :: Bool # direction relative to the shared label + label :: OverlayLabel + + in_result_area :: Bool + in_result_line :: Bool + visited :: Bool + + next_result :: Int32 # next half-edge in the result ring (0 = null) + next_result_max :: Int32 # next half-edge in the result maximal ring (0 = null) + edge_ring :: Int32 # phase 2b OverlayEdgeRing handle (0 = null) + max_edge_ring :: Int32 # phase 2b MaximalEdgeRing handle (0 = null) +end + +_overlay_edge(origin::Int32, dir_pt::P, is_forward::Bool, label::OverlayLabel) where {P} = + OverlayEdge{P}(origin, Int32(0), Int32(0), dir_pt, is_forward, label, + false, false, false, Int32(0), Int32(0), Int32(0), Int32(0)) + +# ## OverlayEdge accessors (ports of the JTS `OverlayEdge` surface, index-based) + +@inline oe_is_forward(edges, i::Integer) = @inbounds edges[i].is_forward +@inline oe_label(edges, i::Integer) = @inbounds edges[i].label + +# Location of `position` for input `index`, resolved for this edge's orientation +# (port of JTS `OverlayEdge.getLocation`). +@inline oe_get_location(edges, i::Integer, index::Integer, position::Integer) = + get_location((@inbounds edges[i].label), index, position, (@inbounds edges[i].is_forward)) + +@inline oe_get_location_boundary_or_line(edges, i::Integer, index::Integer, position::Integer) = + get_location_boundary_or_line((@inbounds edges[i].label), index, position, (@inbounds edges[i].is_forward)) + +# Result-area marking (ports of the `markInResultArea*` / `unmark*` family). +@inline oe_in_result_area(edges, i::Integer) = @inbounds edges[i].in_result_area +@inline oe_in_result_area_both(edges, i::Integer) = + (@inbounds edges[i].in_result_area) && (@inbounds edges[he_sym(edges, i)].in_result_area) +@inline function oe_mark_in_result_area!(edges, i::Integer) + @inbounds edges[i].in_result_area = true; return nothing +end +@inline function oe_mark_in_result_area_both!(edges, i::Integer) + @inbounds edges[i].in_result_area = true + @inbounds edges[he_sym(edges, i)].in_result_area = true + return nothing +end +@inline function oe_unmark_from_result_area_both!(edges, i::Integer) + @inbounds edges[i].in_result_area = false + @inbounds edges[he_sym(edges, i)].in_result_area = false + return nothing +end + +# Result-line marking (marks both members of the pair, per JTS). +@inline oe_in_result_line(edges, i::Integer) = @inbounds edges[i].in_result_line +@inline function oe_mark_in_result_line!(edges, i::Integer) + @inbounds edges[i].in_result_line = true + @inbounds edges[he_sym(edges, i)].in_result_line = true + return nothing +end +@inline oe_in_result(edges, i::Integer) = + (@inbounds edges[i].in_result_area) || (@inbounds edges[i].in_result_line) +@inline oe_in_result_either(edges, i::Integer) = + oe_in_result(edges, i) || oe_in_result(edges, he_sym(edges, i)) + +# Visited flag (marks both members, per JTS `markVisitedBoth`). +@inline oe_is_visited(edges, i::Integer) = @inbounds edges[i].visited +@inline function oe_mark_visited_both!(edges, i::Integer) + @inbounds edges[i].visited = true + @inbounds edges[he_sym(edges, i)].visited = true + return nothing +end + +# Result-ring linkage getters/setters (`0` = null). +@inline oe_next_result(edges, i::Integer) = @inbounds edges[i].next_result +@inline function oe_set_next_result!(edges, i::Integer, e::Integer) + @inbounds edges[i].next_result = Int32(e); return nothing +end +@inline oe_is_result_linked(edges, i::Integer) = (@inbounds edges[i].next_result) != 0 +@inline oe_next_result_max(edges, i::Integer) = @inbounds edges[i].next_result_max +@inline function oe_set_next_result_max!(edges, i::Integer, e::Integer) + @inbounds edges[i].next_result_max = Int32(e); return nothing +end +@inline oe_is_result_max_linked(edges, i::Integer) = (@inbounds edges[i].next_result_max) != 0 + +# ========================================================================== +# OverlayGraph (JTS `OverlayGraph`): build from arrangement + sources +# ========================================================================== + +""" + OverlayGraph{P} + +The topology graph of an overlay operation (port of JTS `OverlayGraph`). Holds +the arrangement it was built from, the vector of half-edges (both orientations of +every merged edge), and one representative outgoing half-edge index per node id +(`node_edges[nid]`, `0` if the node has no incident edges). `P` is the manifold +kernel point type, so the graph is type-erased over input geometry types. +""" +struct OverlayGraph{P} + arr :: NodedArrangement{P} + edges :: Vector{OverlayEdge{P}} + node_edges :: Vector{Int32} +end + +""" + OverlayGraph(m, arr::NodedArrangement, sources) -> OverlayGraph + +Build the overlay graph from a noded arrangement and its `EdgeSourceInfo` table. +Coincident noded edges are merged (JTS `Edge.merge` semantics), each merged edge +becomes a symmetric `OverlayEdge` pair sharing one label, and every node's star +is ordered CCW about its symbolic apex via the exact kernel comparator. +""" +function OverlayGraph(m::Manifold, arr::NodedArrangement{P}, sources::Vector{EdgeSourceInfo}; + exact = True()) where {P} + merged = _merge_noded_edges(arr, sources) + nnodes = num_nodes(arr) + edges = Vector{OverlayEdge{P}}() + sizehint!(edges, 2 * length(merged)) + stars = [Int32[] for _ in 1:nnodes] + for me in merged + label = _create_label(me) + ss = arr.segstrings[me.string_idx] + #-- direction points are the parent segment's original endpoints, so the + #-- forward half-edge (origin node_lo) heads toward the node_hi side and + #-- the reverse (origin node_hi) toward the node_lo side (design §3.1). + fwd_dir = ss.pts[me.seg_idx + 1] + bwd_dir = ss.pts[me.seg_idx] + push!(edges, _overlay_edge(me.node_lo, fwd_dir, true, label)) + i_fwd = Int32(length(edges)) + push!(edges, _overlay_edge(me.node_hi, bwd_dir, false, label)) + i_rev = Int32(length(edges)) + he_link!(edges, i_fwd, i_rev) + push!(stars[me.node_lo], i_fwd) + push!(stars[me.node_hi], i_rev) + end + node_edges = zeros(Int32, nnodes) + _order_all_stars!(m, edges, arr.nodes.keys, stars, node_edges; exact) + return OverlayGraph{P}(arr, edges, node_edges) +end + +# Convenience: build the sources and the graph directly from an arrangement. +function OverlayGraph(m::Manifold, arr::NodedArrangement; exact = True()) + return OverlayGraph(m, arr, _edge_source_infos(m, arr; exact); exact) +end + +# Order every node's star once (function barrier: the abstract `m` dispatches into +# the exact comparator here, off the type-stable build loop). +function _order_all_stars!(m::Manifold, edges, keys, stars, node_edges; exact) + for nid in eachindex(stars) + star = stars[nid] + isempty(star) && continue + he_order_star!(m, edges, keys, star; exact) + @inbounds node_edges[nid] = star[1] + end + return nothing +end + +# ## Graph queries (ports of the `OverlayGraph` accessor surface) + +# All half-edges (both orientations), matching JTS `getEdges()`. +graph_edges(g::OverlayGraph) = g.edges + +# A representative outgoing half-edge index for node `nid`, or `0` (JTS +# `getNodeEdge`). +graph_node_edge(g::OverlayGraph, nid::Integer) = @inbounds g.node_edges[nid] + +# The representative node edges (one outgoing half-edge per non-empty node), for +# node-star iteration (JTS `getNodeEdges()`). +graph_node_edges(g::OverlayGraph) = Int32[e for e in g.node_edges if e != 0] + +# The half-edge indices marked as being in the result area (JTS +# `getResultAreaEdges`). +graph_result_area_edges(g::OverlayGraph) = + Int32[i for i in eachindex(g.edges) if g.edges[i].in_result_area] diff --git a/src/methods/clipping/overlayng/overlay_label.jl b/src/methods/clipping/overlayng/overlay_label.jl new file mode 100644 index 000000000..df2d008d7 --- /dev/null +++ b/src/methods/clipping/overlayng/overlay_label.jl @@ -0,0 +1,211 @@ +# # OverlayLabel — the per-edge topological label (port of JTS `OverlayLabel`) +# +# Phase 2a of the OverlayNG port (design doc `2026-07-16-overlayng-noding-substrate.md`, +# §3). A faithful port of `operation/overlayng/OverlayLabel.java`: a mutable plain +# struct (NOT bit-packed), one instance shared between the two oppositely-oriented +# `OverlayEdge`s of a symmetric pair. Orientation-sensitive accessors take the +# containing edge's `is_forward` flag and swap Left/Right when it is false. +# +# A label records, for each of the (up to) two input geometries (index `0` = A, +# `1` = B, matching JTS), one of four states via its `dim` field: +# - Boundary (`DIM_A`) : an area-boundary edge; `loc_left`/`loc_right` set, +# `loc_line = LOC_INTERIOR`, `is_hole` = ring role. +# - Collapse (`DIM_COLLAPSE`): two or more coincident area edges summed to a +# zero depth delta; `loc_line` filled later. +# - Line (`DIM_L`) : a linear edge; only `loc_line` is meaningful. +# - NotPart (`DIM_NOT_PART`): the edge is not part of this input. +# +# Everything here is internal to GeometryOps — nothing is exported. + +# ## Edge-state dimension codes (port of `OverlayLabel.DIM_*`) +# +# `DIM_BOUNDARY` / `DIM_LINE` reuse the DE9IM `DIM_A` (2) / `DIM_L` (1) values, +# which JTS's `OverlayLabel` constants numerically coincide with (JTS +# `DIM_BOUNDARY == Dimension.A`, `DIM_LINE == Dimension.L`); only the two overlay- +# specific states need fresh values. +const DIM_NOT_PART = Int8(-1) # JTS OverlayLabel.DIM_NOT_PART == DIM_UNKNOWN == Dimension.FALSE +const DIM_COLLAPSE = Int8(3) # JTS OverlayLabel.DIM_COLLAPSE + +# `Position` codes (`POS_ON`/`POS_LEFT`/`POS_RIGHT`) and location codes +# (`LOC_INTERIOR`/`LOC_BOUNDARY`/`LOC_EXTERIOR`/`LOC_NONE`) are the RelateNG ports +# already in scope (relate_node.jl / de9im.jl); `LOC_NONE` is JTS `LOC_UNKNOWN`. + +""" + OverlayLabel + +The topological label of one edge of the overlay graph (port of JTS +`OverlayLabel`). Mutable plain fields, one instance shared between a symmetric +`OverlayEdge` pair; the `is_forward`-parameterized accessors swap Left/Right for +the reverse half-edge. Index `0` selects input A, index `1` input B. +""" +mutable struct OverlayLabel + a_dim :: Int8 + a_is_hole :: Bool + a_loc_left :: Int8 + a_loc_right:: Int8 + a_loc_line :: Int8 + + b_dim :: Int8 + b_is_hole :: Bool + b_loc_left :: Int8 + b_loc_right:: Int8 + b_loc_line :: Int8 +end + +# Uninitialized label: both inputs NotPart, all locations unknown. +OverlayLabel() = OverlayLabel(DIM_NOT_PART, false, LOC_NONE, LOC_NONE, LOC_NONE, + DIM_NOT_PART, false, LOC_NONE, LOC_NONE, LOC_NONE) + +Base.copy(l::OverlayLabel) = OverlayLabel(l.a_dim, l.a_is_hole, l.a_loc_left, l.a_loc_right, l.a_loc_line, + l.b_dim, l.b_is_hole, l.b_loc_left, l.b_loc_right, l.b_loc_line) + +# ## Initializers (port of `initBoundary`/`initCollapse`/`initLine`/`initNotPart`) + +# Area-boundary edge: side locations known, `loc_line = LOC_INTERIOR`. +function init_boundary!(l::OverlayLabel, index::Integer, loc_left, loc_right, is_hole::Bool) + if index == 0 + l.a_dim = DIM_A; l.a_is_hole = is_hole + l.a_loc_left = loc_left; l.a_loc_right = loc_right; l.a_loc_line = LOC_INTERIOR + else + l.b_dim = DIM_A; l.b_is_hole = is_hole + l.b_loc_left = loc_left; l.b_loc_right = loc_right; l.b_loc_line = LOC_INTERIOR + end + return l +end + +# Collapsed area edge: location unknown, resolved later from the graph topology. +function init_collapse!(l::OverlayLabel, index::Integer, is_hole::Bool) + if index == 0 + l.a_dim = DIM_COLLAPSE; l.a_is_hole = is_hole + else + l.b_dim = DIM_COLLAPSE; l.b_is_hole = is_hole + end + return l +end + +# Line edge: only the line location is meaningful, initialized unknown. +function init_line!(l::OverlayLabel, index::Integer) + if index == 0 + l.a_dim = DIM_L; l.a_loc_line = LOC_NONE + else + l.b_dim = DIM_L; l.b_loc_line = LOC_NONE + end + return l +end + +# Not part of this input (locations assumed already unknown). +function init_not_part!(l::OverlayLabel, index::Integer) + index == 0 ? (l.a_dim = DIM_NOT_PART) : (l.b_dim = DIM_NOT_PART) + return l +end + +# ## Location setters (used during label propagation, phase 2b) + +function set_location_line!(l::OverlayLabel, index::Integer, loc) + index == 0 ? (l.a_loc_line = loc) : (l.b_loc_line = loc) + return l +end + +function set_location_all!(l::OverlayLabel, index::Integer, loc) + if index == 0 + l.a_loc_line = loc; l.a_loc_left = loc; l.a_loc_right = loc + else + l.b_loc_line = loc; l.b_loc_left = loc; l.b_loc_right = loc + end + return l +end + +# A collapsed edge with no boundary information takes its parent ring role: +# a hole collapse is INTERIOR, a shell collapse is EXTERIOR. +function set_location_collapse!(l::OverlayLabel, index::Integer) + loc = is_hole(l, index) ? LOC_INTERIOR : LOC_EXTERIOR + index == 0 ? (l.a_loc_line = loc) : (l.b_loc_line = loc) + return l +end + +# ## State predicates (port of the `OverlayLabel` boolean surface) + +dimension(l::OverlayLabel, index::Integer) = index == 0 ? l.a_dim : l.b_dim + +is_line(l::OverlayLabel) = l.a_dim == DIM_L || l.b_dim == DIM_L +is_line(l::OverlayLabel, index::Integer) = (index == 0 ? l.a_dim : l.b_dim) == DIM_L +is_linear(l::OverlayLabel, index::Integer) = + (d = index == 0 ? l.a_dim : l.b_dim; d == DIM_L || d == DIM_COLLAPSE) +is_known(l::OverlayLabel, index::Integer) = (index == 0 ? l.a_dim : l.b_dim) != DIM_NOT_PART +is_not_part(l::OverlayLabel, index::Integer) = (index == 0 ? l.a_dim : l.b_dim) == DIM_NOT_PART + +is_boundary_either(l::OverlayLabel) = l.a_dim == DIM_A || l.b_dim == DIM_A +is_boundary_both(l::OverlayLabel) = l.a_dim == DIM_A && l.b_dim == DIM_A +is_boundary(l::OverlayLabel, index::Integer) = (index == 0 ? l.a_dim : l.b_dim) == DIM_A + +# A collapse coincident with the other input's (non-collapsed) boundary. +is_boundary_collapse(l::OverlayLabel) = is_line(l) ? false : !is_boundary_both(l) + +# Two areas touching along their common boundary (opposite Right locations). +is_boundary_touch(l::OverlayLabel) = + is_boundary_both(l) && + get_location(l, 0, POS_RIGHT, true) != get_location(l, 1, POS_RIGHT, true) + +is_boundary_singleton(l::OverlayLabel) = + (l.a_dim == DIM_A && l.b_dim == DIM_NOT_PART) || + (l.b_dim == DIM_A && l.a_dim == DIM_NOT_PART) + +is_line_location_unknown(l::OverlayLabel, index::Integer) = + (index == 0 ? l.a_loc_line : l.b_loc_line) == LOC_NONE + +is_line_in_area(l::OverlayLabel, index::Integer) = + (index == 0 ? l.a_loc_line : l.b_loc_line) == LOC_INTERIOR + +is_hole(l::OverlayLabel, index::Integer) = index == 0 ? l.a_is_hole : l.b_is_hole + +is_collapse(l::OverlayLabel, index::Integer) = dimension(l, index) == DIM_COLLAPSE + +is_interior_collapse(l::OverlayLabel) = + (l.a_dim == DIM_COLLAPSE && l.a_loc_line == LOC_INTERIOR) || + (l.b_dim == DIM_COLLAPSE && l.b_loc_line == LOC_INTERIOR) + +is_collapse_and_not_part_interior(l::OverlayLabel) = + (l.a_dim == DIM_COLLAPSE && l.b_dim == DIM_NOT_PART && l.b_loc_line == LOC_INTERIOR) || + (l.b_dim == DIM_COLLAPSE && l.a_dim == DIM_NOT_PART && l.a_loc_line == LOC_INTERIOR) + +is_line_interior(l::OverlayLabel, index::Integer) = + (index == 0 ? l.a_loc_line : l.b_loc_line) == LOC_INTERIOR + +get_line_location(l::OverlayLabel, index::Integer) = index == 0 ? l.a_loc_line : l.b_loc_line + +has_sides(l::OverlayLabel, index::Integer) = + index == 0 ? (l.a_loc_left != LOC_NONE || l.a_loc_right != LOC_NONE) : + (l.b_loc_left != LOC_NONE || l.b_loc_right != LOC_NONE) + +# ## Orientation-sensitive location accessor (the shared-label L/R swap) + +""" + get_location(label, index, position, is_forward) -> Int8 + +The location of `position` (`POS_LEFT` / `POS_RIGHT` / `POS_ON`) of input `index`, +for a containing half-edge whose orientation is `is_forward`. When the edge is the +reverse half-edge (`is_forward == false`) the Left/Right stored sides are swapped, +so a single label serves both members of a symmetric pair (port of JTS +`OverlayLabel.getLocation(index, position, isForward)`). +""" +function get_location(l::OverlayLabel, index::Integer, position::Integer, is_forward::Bool) + if index == 0 + position == POS_LEFT && return is_forward ? l.a_loc_left : l.a_loc_right + position == POS_RIGHT && return is_forward ? l.a_loc_right : l.a_loc_left + position == POS_ON && return l.a_loc_line + else + position == POS_LEFT && return is_forward ? l.b_loc_left : l.b_loc_right + position == POS_RIGHT && return is_forward ? l.b_loc_right : l.b_loc_left + position == POS_ON && return l.b_loc_line + end + return LOC_NONE +end + +# The linear (ON) location of a source (port of the 2-arg `getLocation`). +get_location(l::OverlayLabel, index::Integer) = index == 0 ? l.a_loc_line : l.b_loc_line + +# For a boundary edge the side location; otherwise the line location — the +# quantity `markInResultArea` tests (port of `getLocationBoundaryOrLine`). +get_location_boundary_or_line(l::OverlayLabel, index::Integer, position::Integer, is_forward::Bool) = + is_boundary(l, index) ? get_location(l, index, position, is_forward) : + get_line_location(l, index) diff --git a/src/methods/clipping/overlayng/overlay_labeller.jl b/src/methods/clipping/overlayng/overlay_labeller.jl new file mode 100644 index 000000000..123870cd7 --- /dev/null +++ b/src/methods/clipping/overlayng/overlay_labeller.jl @@ -0,0 +1,242 @@ +# # OverlayLabeller — the five-pass edge labelling (port of JTS `OverlayLabeller`) +# +# Phase 2b of the OverlayNG port (design doc `2026-07-16-overlayng-noding-substrate.md`, +# §3). A faithful port of `operation/overlayng/OverlayLabeller.java`, operating on +# the phase-2a index-based `OverlayGraph`. The single shared `OverlayLabel` between a +# symmetric pair means every `set_location_*!` here is seen by both half-edges — +# JTS relies on this, and so does this port. +# +# Also hosts the overlay op-code enum and `_is_result_of_op` (JTS +# `OverlayNG.isResultOfOp`, the boolean core used by both the labeller and the line +# builder), and the topology-error type (JTS `TopologyException`). +# +# Everything here is internal to GeometryOps — nothing is exported. + +# ## Overlay operation codes and the op boolean core + +# The four set-theoretic overlay operations (port of the `OverlayNG.INTERSECTION` +# / `UNION` / `DIFFERENCE` / `SYMDIFFERENCE` op codes). +@enum _OverlayOpCode::UInt8 begin + OVERLAY_INTERSECTION + OVERLAY_UNION + OVERLAY_DIFFERENCE + OVERLAY_SYMDIFFERENCE +end + +""" + _OverlayTopologyError(msg) + +A robustness/topology error raised by the overlay engine (port of JTS +`TopologyException`). Signals an inconsistency the builder could not resolve +(e.g. a side-location conflict during area propagation, or a ring that cannot +be closed). +""" +struct _OverlayTopologyError <: Exception + msg::String +end +Base.showerror(io::IO, e::_OverlayTopologyError) = print(io, "OverlayTopologyError: ", e.msg) + +# Port of `OverlayNG.isResultOfOp`: whether a point with the given per-input +# locations lies in the result of `op`. `LOC_BOUNDARY` counts as `LOC_INTERIOR`. +@inline function _is_result_of_op(op::_OverlayOpCode, loc0::Integer, loc1::Integer) + loc0 == LOC_BOUNDARY && (loc0 = LOC_INTERIOR) + loc1 == LOC_BOUNDARY && (loc1 = LOC_INTERIOR) + if op == OVERLAY_INTERSECTION + return loc0 == LOC_INTERIOR && loc1 == LOC_INTERIOR + elseif op == OVERLAY_UNION + return loc0 == LOC_INTERIOR || loc1 == LOC_INTERIOR + elseif op == OVERLAY_DIFFERENCE + return loc0 == LOC_INTERIOR && loc1 != LOC_INTERIOR + else # OVERLAY_SYMDIFFERENCE + return (loc0 == LOC_INTERIOR && loc1 != LOC_INTERIOR) || + (loc0 != LOC_INTERIOR && loc1 == LOC_INTERIOR) + end +end + +# ## computeLabelling — the five passes, in order (port of `computeLabelling`) + +function _compute_labelling!(g::OverlayGraph, input) + edges = g.edges + for ne in graph_node_edges(g) + _propagate_area_locations!(edges, input, ne, 0) + _input_has_edges(input, 1) && _propagate_area_locations!(edges, input, ne, 1) + end + _label_connected_linear_edges!(edges, input) + _label_collapsed_edges!(edges) + _label_connected_linear_edges!(edges, input) + _label_disconnected_edges!(g, input) + return nothing +end + +# ### Pass 1: area-node propagation (port of `propagateAreaLocations`) + +# Scans a node's star CCW, propagating side labels for one area input to every +# edge whose location for that input is still unknown. A side-location conflict +# between two boundary edges is a topology error (port of the JTS check). +function _propagate_area_locations!(edges, input, node_edge::Integer, gi::Integer) + _input_is_area(input, gi) || return nothing + #-- one-edge node: nothing to propagate (dangling edge) + he_degree(edges, node_edge) == 1 && return nothing + + e_start = _find_propagation_start_edge(edges, node_edge, gi) + e_start == 0 && return nothing + + curr_loc = oe_get_location(edges, e_start, gi, POS_LEFT) + e = he_onext(edges, e_start) + while e != e_start + label = oe_label(edges, e) + if !is_boundary(label, gi) + #-- non-boundary edge: its location relative to this area is now known + set_location_line!(label, gi, curr_loc) + else + loc_right = oe_get_location(edges, e, gi, POS_RIGHT) + loc_right == curr_loc || + throw(_OverlayTopologyError("side location conflict: arg $gi")) + loc_left = oe_get_location(edges, e, gi, POS_LEFT) + #-- loc_left == LOC_NONE should never happen for a boundary edge + curr_loc = loc_left + end + e = he_onext(edges, e) + end + return nothing +end + +# Port of `findPropagationStartEdge`: a boundary edge for `gi` in the node's +# star, or `0` if none. +function _find_propagation_start_edge(edges, node_edge::Integer, gi::Integer) + e = Int32(node_edge) + while true + is_boundary(oe_label(edges, e), gi) && return e + e = he_onext(edges, e) + e == node_edge && break + end + return Int32(0) +end + +# ### Pass 2/4: connected linear propagation (port of `labelConnectedLinearEdges`) + +function _label_connected_linear_edges!(edges, input) + _propagate_linear_locations!(edges, input, 0) + _input_has_edges(input, 1) && _propagate_linear_locations!(edges, input, 1) + return nothing +end + +# BFS over linear (line or collapse) edges with a known location, propagating +# that location to connected unknown edges (port of `propagateLinearLocations`). +function _propagate_linear_locations!(edges, input, gi::Integer) + stack = Int32[] + for i in eachindex(edges) + lbl = oe_label(edges, i) + if is_linear(lbl, gi) && !is_line_location_unknown(lbl, gi) + push!(stack, Int32(i)) + end + end + isempty(stack) && return nothing + + is_input_line = _input_is_line(input, gi) + while !isempty(stack) + e_node = pop!(stack) + _propagate_linear_at_node!(edges, e_node, gi, is_input_line, stack) + end + return nothing +end + +# Port of `propagateLinearLocationAtNode`. Line parents propagate EXTERIOR only. +function _propagate_linear_at_node!(edges, e_node::Integer, gi::Integer, + is_input_line::Bool, stack::Vector{Int32}) + line_loc = get_line_location(oe_label(edges, e_node), gi) + #-- a Line parent only propagates EXTERIOR locations + is_input_line && line_loc != LOC_EXTERIOR && return nothing + + e = he_onext(edges, e_node) + while e != e_node + label = oe_label(edges, e) + if is_line_location_unknown(label, gi) + set_location_line!(label, gi, line_loc) + #-- continue the traversal from the far node (don't re-add e itself) + push!(stack, he_sym(edges, e)) + end + e = he_onext(edges, e) + end + return nothing +end + +# ### Pass 3: collapsed-edge ring-role labelling (port of `labelCollapsedEdges`) + +function _label_collapsed_edges!(edges) + for i in eachindex(edges) + label = oe_label(edges, i) + is_line_location_unknown(label, 0) && _label_collapsed_edge!(label, 0) + is_line_location_unknown(label, 1) && _label_collapsed_edge!(label, 1) + end + return nothing +end + +function _label_collapsed_edge!(label::OverlayLabel, gi::Integer) + is_collapse(label, gi) || return nothing + #-- disconnected collapsed edge: label from its parent ring role (shell/hole) + set_location_collapse!(label, gi) + return nothing +end + +# ### Pass 5: disconnected-edge PIP labelling (port of `labelDisconnectedEdges`) + +function _label_disconnected_edges!(g::OverlayGraph, input) + edges = g.edges + for i in eachindex(edges) + label = oe_label(edges, i) + is_line_location_unknown(label, 0) && _label_disconnected_edge!(g, input, i, 0) + is_line_location_unknown(label, 1) && _label_disconnected_edge!(g, input, i, 1) + end + return nothing +end + +# Locates a disconnected edge against the ORIGINAL input area (design §3 +# amendment 7: never against the reduced/collapsed linework), using both +# endpoints for robustness (port of `labelDisconnectedEdge` + +# `locateEdgeBothEnds`). +function _label_disconnected_edge!(g::OverlayGraph, input, i::Integer, gi::Integer) + label = oe_label(g.edges, i) + if !_input_is_area(input, gi) + #-- non-area target: a disconnected edge must be EXTERIOR + set_location_all!(label, gi, LOC_EXTERIOR) + return nothing + end + loc_orig = _input_locate_in_area(input, gi, node_point(g.arr, he_origin(g.edges, i))) + loc_dest = _input_locate_in_area(input, gi, node_point(g.arr, he_dest(g.edges, i))) + is_int = loc_orig != LOC_EXTERIOR && loc_dest != LOC_EXTERIOR + set_location_all!(label, gi, is_int ? LOC_INTERIOR : LOC_EXTERIOR) + return nothing +end + +# ## Result-area marking (ports of `markResultAreaEdges` / `unmarkDuplicate…`) + +function _mark_result_area_edges!(g::OverlayGraph, op::_OverlayOpCode) + edges = g.edges + for i in eachindex(edges) + _mark_in_result_area!(edges, i, op) + end + return nothing +end + +# Port of `markInResultArea`: mark an edge whose right-side (boundary) or line +# location makes it part of the result-area boundary under `op`. +@inline function _mark_in_result_area!(edges, i::Integer, op::_OverlayOpCode) + label = oe_label(edges, i) + if is_boundary_either(label) && _is_result_of_op(op, + oe_get_location_boundary_or_line(edges, i, 0, POS_RIGHT), + oe_get_location_boundary_or_line(edges, i, 1, POS_RIGHT)) + oe_mark_in_result_area!(edges, i) + end + return nothing +end + +# Port of `unmarkDuplicateEdgesFromResultArea`: an edge whose sym is also in the +# result area cancels (merges edge-adjacent result areas per polygon validity). +function _unmark_duplicate_edges_from_result_area!(g::OverlayGraph) + edges = g.edges + for i in eachindex(edges) + oe_in_result_area_both(edges, i) && oe_unmark_from_result_area_both!(edges, i) + end + return nothing +end diff --git a/src/methods/clipping/overlayng/overlay_ng.jl b/src/methods/clipping/overlayng/overlay_ng.jl new file mode 100644 index 000000000..8abb7cd24 --- /dev/null +++ b/src/methods/clipping/overlayng/overlay_ng.jl @@ -0,0 +1,241 @@ +# # OverlayNG driver — the internal end-to-end overlay engine +# +# Phase 2b of the OverlayNG port (design doc §3, §4 preview). Ties the phase-1 +# noding substrate, the phase-2a graph, and the phase-2b labeller/builders into +# the internal driver `_overlay_ng`. Ports the engine core of `OverlayNG.java` +# (`computeEdgeOverlay` phase order, `extractResult` dimensional priority) and +# `OverlayUtil.java` (`resultDimension`, result assembly), plus the input model +# of `InputGeometry.java`. Skipped per the design: ElevationModel, +# FastOverlayFilter, strict mode, precision/PM, RingClipper/LineLimiter, +# OverlayPoints/OverlayMixedPoints (phase 3). +# +# This is NOT public: no exports, no `@ref` docstrings. The public opt-in +# `OverlayNG{M}` algorithm and the differential-validation harnesses are phase 3; +# the existing `intersection`/`union`/`difference` are untouched. +# +# Accepted inputs: line (`LineString`/`LinearRing`/`MultiLineString`) and area +# (`Polygon`/`MultiPolygon`) geometries, in any A×B combination. Point inputs +# and geometry collections are rejected with a clear error (phase 3). + +# ## Input model (port of `InputGeometry`) + +# The two overlay operands with their dimensions, lazily-built area locators, and +# empty flags. `a`/`b` are boxed (`Any`) — only the cold locator-build path reads +# them; the hot pipeline runs on the type-erased graph. +mutable struct _OverlayInput{M <: Manifold, E} + m::M + a::Any + b::Any + dim_a::Int + dim_b::Int + exact::E + empty_a::Bool + empty_b::Bool + loc_a::Any # Union{Nothing, IndexedPointInAreaLocator}, lazy + loc_b::Any +end + +@inline _input_dim(input::_OverlayInput, gi::Integer) = gi == 0 ? input.dim_a : input.dim_b +@inline _input_is_area(input::_OverlayInput, gi::Integer) = _input_dim(input, gi) == 2 +@inline _input_is_line(input::_OverlayInput, gi::Integer) = _input_dim(input, gi) == 1 +@inline _input_has_edges(input::_OverlayInput, gi::Integer) = _input_dim(input, gi) > 0 + +# Port of `InputGeometry.getAreaIndex`: the index of an area input, or `-1`. +@inline _input_area_index(input::_OverlayInput) = + input.dim_a == 2 ? 0 : (input.dim_b == 2 ? 1 : -1) + +# Port of `InputGeometry.locatePointInArea`: locate an emitted point against the +# ORIGINAL input area (design §3 amendment 7). Empty geometries locate EXTERIOR; +# the indexed locator is built once per side on first use. +function _input_locate_in_area(input::_OverlayInput, gi::Integer, pt) + if gi == 0 + input.empty_a && return LOC_EXTERIOR + input.loc_a === nothing && + (input.loc_a = IndexedPointInAreaLocator(input.m, input.a; exact = input.exact)) + return locate(input.loc_a, pt) + else + input.empty_b && return LOC_EXTERIOR + input.loc_b === nothing && + (input.loc_b = IndexedPointInAreaLocator(input.m, input.b; exact = input.exact)) + return locate(input.loc_b, pt) + end +end + +# Dimension of an overlay operand: 2 (area), 1 (line), 0 (point). Geometry +# collections and other traits are unsupported here (phase 3). +function _overlay_dimension(geom) + t = GI.trait(geom) + if t isa GI.PolygonTrait || t isa GI.MultiPolygonTrait + return 2 + elseif t isa GI.LineStringTrait || t isa GI.LinearRingTrait || t isa GI.MultiLineStringTrait + return 1 + elseif t isa GI.PointTrait || t isa GI.MultiPointTrait + return 0 + end + throw(ArgumentError("_overlay_ng: unsupported input geometry trait $(typeof(t))")) +end + +# ## The driver (port of `getResult` / `computeEdgeOverlay`) + +""" + _overlay_ng(m, op::_OverlayOpCode, a, b; exact=True(), tree_a=nothing, tree_b=nothing) + +Compute the overlay of `a` and `b` under `op` on manifold `m`, returning a +GeoInterface geometry. Internal engine entry point for OverlayNG phase 2b — +line and area inputs are supported (any A×B combination); point inputs raise an +error (phase 3). `tree_a`/`tree_b` accept caller-prebuilt segment indices +(threaded to the noding substrate). +""" +function _overlay_ng(m::Manifold, op::_OverlayOpCode, a, b; + exact = True(), tree_a = nothing, tree_b = nothing) + dim_a = _overlay_dimension(a) + dim_b = _overlay_dimension(b) + (dim_a < 1 || dim_b < 1) && throw(ArgumentError( + "_overlay_ng supports line and area inputs only (got dimensions " * + "$dim_a and $dim_b); point inputs are handled in phase 3")) + + input = _OverlayInput(m, a, b, dim_a, dim_b, exact, GI.isempty(a), GI.isempty(b), + nothing, nothing) + + #-- empty-input / disjoint-envelope short circuits (port of isEmptyResult) + er = _empty_result_short_circuit(m, op, input) + er === nothing || return er + + arr = NodedArrangement(m, a, b; exact, tree_a, tree_b) + g = OverlayGraph(m, arr; exact) + + _compute_labelling!(g, input) + _mark_result_area_edges!(g, op) + _unmark_duplicate_edges_from_result_area!(g) + + return _extract_result(m, op, g, input; exact) +end + +# ## Result extraction (port of `extractResult`) + +function _extract_result(m::Manifold, op::_OverlayOpCode, g::OverlayGraph, input; exact) + result_area_edges = graph_result_area_edges(g) + polys = _build_polygons(m, g, result_area_edges; exact) + has_result_area = !isempty(polys) + + #-- non-strict semantics always allow result lines + lines = _build_lines(m, g, input, has_result_area, op; exact) + + #-- only Intersection produces points from non-point inputs + points = op == OVERLAY_INTERSECTION ? _build_points(g) : Tuple{Float64, Float64}[] + + if isempty(polys) && isempty(lines) && isempty(points) + return _resolve_empty_result(m, op, input) + end + return _create_result_geometry(polys, lines, points) +end + +# Port of `OverlayUtil.createResultGeometry` + `GeometryFactory.buildGeometry`: +# the most specific geometry over the A, L, P components. +function _create_result_geometry(polys, lines, points) + has_p = !isempty(polys) + has_l = !isempty(lines) + has_pt = !isempty(points) + if has_p && !has_l && !has_pt + return length(polys) == 1 ? polys[1] : GI.MultiPolygon(polys) + elseif !has_p && has_l && !has_pt + return length(lines) == 1 ? lines[1] : GI.MultiLineString(lines) + elseif !has_p && !has_l && has_pt + return length(points) == 1 ? GI.Point(points[1]) : GI.MultiPoint(points) + end + #-- mixed dimensions: a geometry collection in A, L, P order + comps = Any[] + append!(comps, polys) + append!(comps, lines) + for p in points + push!(comps, GI.Point(p)) + end + return GI.GeometryCollection(comps) +end + +# ## Empty / full-sphere handling + +# Port of `OverlayUtil.isEmptyResult` (the input-driven short circuit). Returns +# an empty result geometry, or `nothing` if the pipeline must run. +function _empty_result_short_circuit(m::Manifold, op::_OverlayOpCode, input::_OverlayInput) + if op == OVERLAY_INTERSECTION + (input.empty_a || input.empty_b) && return _empty_result(op, input) + #-- disjoint-envelope reject (planar only; the spherical box is unreliable) + m isa Planar && _env_disjoint(input.a, input.b) && return _empty_result(op, input) + elseif op == OVERLAY_DIFFERENCE + input.empty_a && return _empty_result(op, input) + else # UNION / SYMDIFFERENCE + input.empty_a && input.empty_b && return _empty_result(op, input) + end + return nothing +end + +@inline function _env_disjoint(a, b) + ea = GI.extent(a); eb = GI.extent(b) + (ea === nothing || eb === nothing) && return false + return !Extents.intersects(ea, eb) +end + +# Port of `OverlayUtil.resultDimension`. +function _result_dimension(op::_OverlayOpCode, d0::Integer, d1::Integer) + op == OVERLAY_INTERSECTION && return min(d0, d1) + op == OVERLAY_UNION && return max(d0, d1) + op == OVERLAY_DIFFERENCE && return d0 + return max(d0, d1) # SYMDIFFERENCE +end + +# Resolve a pipeline that produced no components. On the plane an empty result is +# always the empty geometry. On the sphere (design §3 amendment 6) a boundaryless +# area result is ambiguous between empty and the whole sphere; disambiguate by +# locating one input vertex under the op semantics, and reject a full-sphere +# result as unrepresentable under enclosed-region semantics. +function _resolve_empty_result(m::Manifold, op::_OverlayOpCode, input::_OverlayInput) + if m isa Spherical && + _result_dimension(op, input.dim_a, input.dim_b) == 2 && + _covers_everything(m, op, input) + throw(ArgumentError( + "OverlayNG: the overlay result covers the whole sphere, which is not " * + "representable as an enclosed-region polygon under GeometryOps' spherical " * + "overlay semantics (a documented phase-3 refinement point)")) + end + return _empty_result(op, input) +end + +# Whether a boundaryless result covers the entire manifold: since there is no +# result boundary the result is uniform, so evaluating the op at any single point +# decides it. Uses a vertex of an area input (boundary counts as interior). +function _covers_everything(m::Manifold, op::_OverlayOpCode, input::_OverlayInput) + p = _first_area_vertex(input) + loc0 = _input_is_area(input, 0) ? _input_locate_in_area(input, 0, p) : LOC_EXTERIOR + loc1 = _input_is_area(input, 1) ? _input_locate_in_area(input, 1, p) : LOC_EXTERIOR + return _is_result_of_op(op, loc0, loc1) +end + +function _first_area_vertex(input::_OverlayInput) + geom = _input_is_area(input, 0) ? input.a : input.b + p = first(GI.getpoint(geom)) + return (Float64(GI.x(p)), Float64(GI.y(p))) +end + +function _empty_result(op::_OverlayOpCode, input::_OverlayInput) + dim = _result_dimension(op, input.dim_a, input.dim_b) + return _empty_geom(dim) +end + +# Empty geometry of the given dimension (2 → area, 1 → line). Dimension 0 is not +# reachable here (point inputs are rejected up front). GeoInterface's +# auto-detecting wrapper constructors inspect `first(geom)`, so an empty geometry +# must be built through the raw typed (`{Z,M,T,E,C}`) constructor. +const _EmptyRing = Vector{Tuple{Float64, Float64}} +const _EmptyPoly = GI.Polygon{false, false, Vector{_EmptyRing}, Nothing, Nothing} +const _EmptyLine = GI.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing} + +function _empty_geom(dim::Integer) + if dim == 2 + return GI.MultiPolygon{false, false, Vector{_EmptyPoly}, Nothing, Nothing}( + _EmptyPoly[], nothing, nothing) + end + #-- dim == 1 (line); dim 0 unreachable + return GI.MultiLineString{false, false, Vector{_EmptyLine}, Nothing, Nothing}( + _EmptyLine[], nothing, nothing) +end diff --git a/src/methods/clipping/overlayng/polygon_builder.jl b/src/methods/clipping/overlayng/polygon_builder.jl new file mode 100644 index 000000000..fada27966 --- /dev/null +++ b/src/methods/clipping/overlayng/polygon_builder.jl @@ -0,0 +1,150 @@ +# # PolygonBuilder — result polygons from the marked result-area edges +# +# Phase 2b of the OverlayNG port (design doc §3). A port of +# `operation/overlayng/PolygonBuilder.java`: link the result-area edges into +# maximal rings, split those into OGC-valid minimal rings (the ring types and +# linking live in maximal_edge_ring.jl), classify shells vs holes, and assign +# each free hole to its containing shell (design §3 amendment 5: containment via +# the indexed point-in-area locators over kernel points, with an `RTree(STR())` +# over shell extents for pruning, mirroring JTS's HPRtree). +# +# Each `OverlayEdge` is a single segment, so a ring's coordinates are the node +# points of its edges — no `RingClipper`/`LineLimiter` (design §3: replaced by +# whole-ring extent pruning + one PIP per pruned ring). +# +# Everything here is internal to GeometryOps — nothing is exported. + +# Build the result polygons from the graph's result-area edges (port of the +# `PolygonBuilder` constructor + `getPolygons`). +function _build_polygons(m::Manifold, g::OverlayGraph{P}, result_area_edges; exact) where {P} + ctx = _PolyBuilderCtx(m, g.edges, g.arr, exact, _MaxEdgeRing[], _OverlayEdgeRing{P}[], + Int32[], Int32[]) + _build_rings!(ctx, result_area_edges) + return [_ring_to_polygon(ctx, sh) for sh in ctx.shell_list] +end + +# Port of `buildRings`. +function _build_rings!(ctx, result_edges) + #-- design §3 amendment 4: link UNGATED for every in-result edge + for e in result_edges + _link_result_area_max_ring_at_node!(ctx.edges, e) + end + max_rings = _build_maximal_rings!(ctx, result_edges) + for mrid in max_rings + min_rings = _build_minimal_rings!(ctx, ctx.max_rings[mrid]) + _assign_shells_and_holes!(ctx, min_rings) + end + _place_free_holes!(ctx) + return nothing +end + +# Port of `buildMaximalRings`: one `_MaxEdgeRing` per unprocessed in-result +# boundary edge. Returns the max-ring handles. +function _build_maximal_rings!(ctx, edges_iter) + max_rings = Int32[] + for e in edges_iter + if oe_in_result_area(ctx.edges, e) && is_boundary_either(oe_label(ctx.edges, e)) && + ctx.edges[e].max_edge_ring == 0 + id = Int32(length(ctx.max_rings) + 1) + mr = _MaxEdgeRing(id, Int32(e)) + push!(ctx.max_rings, mr) + _attach_max_edges!(ctx, mr) + push!(max_rings, id) + end + end + return max_rings +end + +# Port of `assignShellsAndHoles`: the minimal rings of one maximal ring are +# either a shell + its holes, or a set of (connected) holes whose shell is found +# later (free holes). +function _assign_shells_and_holes!(ctx, min_rings) + shell = _find_single_shell(ctx, min_rings) + if shell != 0 + for er in min_rings + ctx.edge_rings[er].is_hole && _set_shell!(ctx, er, shell) + end + push!(ctx.shell_list, shell) + else + append!(ctx.free_hole_list, min_rings) + end + return nothing +end + +# Port of `findSingleShell`: the single non-hole ring, or `0` (all holes). +function _find_single_shell(ctx, min_rings) + shell = Int32(0) + shell_count = 0 + for er in min_rings + if !ctx.edge_rings[er].is_hole + shell = Int32(er) + shell_count += 1 + end + end + shell_count <= 1 || throw(_OverlayTopologyError("found two shells in EdgeRing list")) + return shell +end + +# Port of `OverlayEdgeRing.setShell` (+ `addHole`). +function _set_shell!(ctx, hole_er::Integer, shell::Integer) + ctx.edge_rings[hole_er].shell = Int32(shell) + shell != 0 && push!(ctx.edge_rings[shell].holes, Int32(hole_er)) + return nothing +end + +# ## Free-hole placement (port of `placeFreeHoles`) + +# Planar: prune candidate shells with an `RTree(STR())` over shell extents +# (design §3 amendment 5, the HPRtree analogue). +function _place_free_holes!(ctx::_PolyBuilderCtx{<:Planar}) + isempty(ctx.free_hole_list) && return nothing + shells = ctx.shell_list + if isempty(shells) + throw(_OverlayTopologyError("unable to assign free hole to a shell")) + end + exts = [_ext_of(ctx.edge_rings[s].bbox) for s in shells] + index = RTree(STR(), collect(shells); extents = exts) + for hole_er in ctx.free_hole_list + ctx.edge_rings[hole_er].shell == 0 || continue + hole_ext = _ext_of(ctx.edge_rings[hole_er].bbox) + cand = Int32[] + SpatialTreeInterface.depth_first_search(Base.Fix1(Extents.intersects, hole_ext), index) do i + push!(cand, index.data[i]) + end + shell = _find_edge_ring_containing(ctx, ctx.edge_rings[hole_er], cand) + shell == 0 && throw(_OverlayTopologyError("unable to assign free hole to a shell")) + _set_shell!(ctx, hole_er, shell) + end + return nothing +end + +# Spherical: the lon/lat extent prune is unreliable near the poles/antimeridian, +# so test every shell (free holes are rare — this path is cold). +function _place_free_holes!(ctx::_PolyBuilderCtx{<:Spherical}) + isempty(ctx.free_hole_list) && return nothing + for hole_er in ctx.free_hole_list + ctx.edge_rings[hole_er].shell == 0 || continue + shell = _find_edge_ring_containing(ctx, ctx.edge_rings[hole_er], ctx.shell_list) + shell == 0 && throw(_OverlayTopologyError("unable to assign free hole to a shell")) + _set_shell!(ctx, hole_er, shell) + end + return nothing +end + +@inline _ext_of(bbox::NTuple{4, Float64}) = + Extents.Extent(X = (bbox[1], bbox[2]), Y = (bbox[3], bbox[4])) + +# ## Polygon assembly (port of `OverlayEdgeRing.toPolygon`) + +# Emit the polygon of one shell ring and its assigned holes. Ring windings are +# left as the graph produced them (JTS does the same); `GO.area` on either +# manifold is orientation-independent, and validity does not depend on winding. +function _ring_to_polygon(ctx, shell_handle::Integer) + sh = ctx.edge_rings[shell_handle] + rings = Vector{Vector{Tuple{Float64, Float64}}}() + push!(rings, sh.ring_pts) + for h in sh.holes + push!(rings, ctx.edge_rings[h].ring_pts) + end + return GI.Polygon(rings) +end diff --git a/test/methods/clipping/overlayng/overlay_graph.jl b/test/methods/clipping/overlayng/overlay_graph.jl new file mode 100644 index 000000000..a5e635ae8 --- /dev/null +++ b/test/methods/clipping/overlayng/overlay_graph.jl @@ -0,0 +1,329 @@ +# Tests for the OverlayNG phase-2a engine core (design §3): the half-edge graph +# over the phase-1 `NodedArrangement` — edge sources / depth deltas, the shared +# `OverlayLabel`, the node-pair edge merger, and the CCW-ordered half-edge stars. + +using Test +import GeometryOps as GO +import GeometryOps: Planar, Spherical, True, False +import GeometryOps: POS_LEFT, POS_RIGHT, POS_ON +import GeometryOps: LOC_INTERIOR, LOC_EXTERIOR, LOC_BOUNDARY, LOC_NONE +import GeometryOps: DIM_A, DIM_L, DIM_COLLAPSE, DIM_NOT_PART +import GeoInterface as GI + +const EX = True() + +# --------------------------------------------------------------------------- +# helpers +# --------------------------------------------------------------------------- + +_crossing(arr) = [i for i in 1:GO.num_nodes(arr) if arr.nodes.keys[i].is_crossing] + +# The onext cycle (outgoing half-edge indices, CCW) at node `nid`. +function star_of(g, nid) + ne = GO.graph_node_edge(g, nid) + ne == 0 && return Int32[] + star = Int32[] + e = ne + while true + push!(star, e) + e = GO.he_onext(g.edges, e) + e == ne && break + end + return star +end + +# Every half-edge has a well-formed sym, every node's star is a CCW-ordered +# permutation of its outgoing edges, and degrees match the arrangement incidence. +function check_graph(m, g) + for i in eachindex(g.edges) + s = GO.he_sym(g.edges, i) + @test GO.he_sym(g.edges, s) == i # sym is involutive + @test s != i # a half-edge is never its own sym + @test g.edges[s].origin == GO.he_dest(g.edges, i) # dest == sym.origin + end + for nid in 1:GO.num_nodes(g.arr) + star = star_of(g, nid) + isempty(star) && continue + incidence = count(i -> g.edges[i].origin == nid, eachindex(g.edges)) + @test length(star) == incidence # degree == incidence + @test GO.he_degree(g.edges, star[1]) == incidence + @test all(i -> g.edges[i].origin == nid, star) # all share the origin + #-- strictly CCW-increasing from the representative (lowest) edge + for t in 1:(length(star) - 1) + @test GO.he_compare_angular(m, g.edges, g.arr.nodes.keys, star[t], star[t + 1]; exact = EX) < 0 + end + end +end + +# --------------------------------------------------------------------------- +# 1. Angular star ordering torture set (§3.1) +# --------------------------------------------------------------------------- + +@testset "degree-6 coincidence-merged node (foreign-direction slow path)" begin + # two A lines and one B line through the origin -> one merged node of degree 6 + # (each line contributes an in- and an out-going half-edge), exercising + # rk_compare_edge_dir's exact-rational foreign-direction path at a crossing apex. + A = GI.MultiLineString([[(-1.0, -1.0), (1.0, 1.0)], [(-1.0, 1.0), (1.0, -1.0)]]) + B = GI.LineString([(-1.0, 0.0), (1.0, 0.0)]) + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = EX) + g = GO.OverlayGraph(m, arr; exact = EX) + cid = findfirst(i -> arr.nodes.keys[i].is_crossing, 1:GO.num_nodes(arr)) + @test cid !== nothing + ne = GO.graph_node_edge(g, cid) + @test GO.he_degree(g.edges, ne) == 6 + check_graph(m, g) + end +end + +@testset "collinear shared boundary -> one merged edge, both labels" begin + # edge-adjacent unit squares sharing the vertex-identical edge x = 2 + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(2.0, 0.0), (4.0, 0.0), (4.0, 2.0), (2.0, 2.0), (2.0, 0.0)]]) + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = EX) + g = GO.OverlayGraph(m, arr; exact = EX) + @test length(_crossing(arr)) == 0 + #-- exactly one merged edge is a boundary of BOTH inputs (its sym pair) + both = [i for i in eachindex(g.edges) if GO.is_boundary_both(g.edges[i].label)] + @test length(both) == 2 + @test GO.he_sym(g.edges, both[1]) == both[2] + check_graph(m, g) + end + #-- per-side locations on the plane: on the shared x=2 edge the WEST side + #-- (x<2) is A-interior/B-exterior and the EAST side is B-interior/A-exterior, + #-- regardless of the canonical ring reorientation applied at ingest — derive + #-- west/east from the forward half-edge's actual heading. + m = Planar() + arr = GO.NodedArrangement(m, A, B; exact = EX) + g = GO.OverlayGraph(m, arr; exact = EX) + bf = first(i for i in eachindex(g.edges) + if GO.is_boundary_both(g.edges[i].label) && g.edges[i].is_forward) + lbl = g.edges[bf].label + o = GO.node_point(g.arr, g.edges[bf].origin) + d = GO.node_point(g.arr, GO.he_dest(g.edges, bf)) + heading_up = d[2] > o[2] # both on x=2 + aL = GO.get_location(lbl, 0, POS_LEFT, true); aR = GO.get_location(lbl, 0, POS_RIGHT, true) + bL = GO.get_location(lbl, 1, POS_LEFT, true); bR = GO.get_location(lbl, 1, POS_RIGHT, true) + a_west, a_east = heading_up ? (aL, aR) : (aR, aL) + b_west, b_east = heading_up ? (bL, bR) : (bR, bL) + @test a_west == LOC_INTERIOR && a_east == LOC_EXTERIOR # A is west + @test b_west == LOC_EXTERIOR && b_east == LOC_INTERIOR # B is east + @test GO.is_boundary_touch(lbl) # two areas meeting +end + +@testset "spherical crossing node tangent ordering (both manifolds)" begin + # a near-equatorial line and a meridian crossing at (0,0): the crossing node + # has degree 4 and its star is CCW-consistent on both manifolds. + A = GI.LineString([(-10.0, 0.0), (10.0, 0.0)]) + B = GI.LineString([(0.0, -10.0), (0.0, 10.0)]) + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = EX) + g = GO.OverlayGraph(m, arr; exact = EX) + cid = findfirst(i -> arr.nodes.keys[i].is_crossing, 1:GO.num_nodes(arr)) + @test cid !== nothing + @test GO.he_degree(g.edges, GO.graph_node_edge(g, cid)) == 4 + check_graph(m, g) + end +end + +# --------------------------------------------------------------------------- +# 2. Depth-delta algebra and Edge.merge semantics (§2.7, §3.3) +# --------------------------------------------------------------------------- + +@testset "depth delta cross-check vs JTS convention" begin + # JTS: canonical depth delta +1 == Exterior-on-Left (material interior RIGHT), + # flipped to -1 otherwise. So material_interior_on_left => -1. + for m in (Planar(), Spherical()) + ccw = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)] + kccw = GO._to_kernel_points(m, GI.LinearRing(ccw)) + d = GO._ring_depth_delta(m, kccw, false; exact = EX) + # CCW shell has material interior on the left -> delta -1 -> L=INT, R=EXT + moil = GO._ring_material_interior_on_left(m, kccw, false; exact = EX) + @test d == (moil ? -1 : 1) + @test GO._location_left(d) == (moil ? LOC_INTERIOR : LOC_EXTERIOR) + @test GO._location_right(d) == (moil ? LOC_EXTERIOR : LOC_INTERIOR) + # reversing the ring flips the delta and the side locations + kcw = GO._to_kernel_points(m, GI.LinearRing(reverse(ccw))) + dcw = GO._ring_depth_delta(m, kcw, false; exact = EX) + @test dcw == -d + @test GO._location_left(dcw) == GO._location_right(d) + end +end + +@testset "EdgeSourceInfo: shell/hole deltas and line" begin + m = Planar() + shell = GO._to_kernel_points(m, GI.LinearRing([(0.0,0.0),(4.0,0.0),(4.0,4.0),(0.0,4.0),(0.0,0.0)])) + # a hole is labelled opposite to a shell of the same stored orientation + d_shell = GO._ring_depth_delta(m, shell, false; exact = EX) + d_hole = GO._ring_depth_delta(m, shell, true; exact = EX) + @test d_hole == -d_shell + # line source carries no side labelling + li = GO.EdgeSourceInfo(Int8(0), DIM_L, false, Int8(0)) + @test li.dim == DIM_L && li.depth_delta == 0 +end + +@testset "Edge.merge: same-dir sum, opposite-dir collapse, hole role" begin + srcA(dd, hole) = GO.EdgeSourceInfo(Int8(0), DIM_A, hole, Int8(dd)) + ne(lo, hi) = GO.NodedEdge(Int32(1), Int32(1), Int32(lo), Int32(hi)) + # same direction: deltas add + base = GO._merge_edge(ne(3, 7), srcA(1, false)) + GO._merge!(base, GO._merge_edge(ne(3, 7), srcA(1, false))) + @test base.a_depth_delta == 2 + # opposite direction (a-b-a spike): deltas cancel -> DIM_COLLAPSE + base = GO._merge_edge(ne(3, 7), srcA(-1, false)) + GO._merge!(base, GO._merge_edge(ne(7, 3), srcA(-1, false))) # reversed + @test base.a_depth_delta == 0 + lbl = GO._create_label(base) + @test GO.is_collapse(lbl, 0) + @test lbl.a_dim == DIM_COLLAPSE + # hole-role merge: a shell contributor makes the merged edge a shell + base = GO._merge_edge(ne(3, 7), srcA(-1, true)) # hole + GO._merge!(base, GO._merge_edge(ne(3, 7), srcA(-1, false))) # shell + @test base.a_is_hole == false + # two holes stay a hole + base = GO._merge_edge(ne(3, 7), srcA(-1, true)) + GO._merge!(base, GO._merge_edge(ne(3, 7), srcA(-1, true))) + @test base.a_is_hole == true +end + +@testset "geometric a-b-a spike ring produces a collapse edge" begin + # a self-touching ring with an out-and-back spike (2,4)->(2,6)->(2,4); + # the noder pairs the two coincident segments and the merger collapses them. + A = GI.Polygon([[(0.0, 0.0), (4.0, 0.0), (4.0, 4.0), (2.0, 4.0), (2.0, 6.0), + (2.0, 4.0), (0.0, 4.0), (0.0, 0.0)]]) + B = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]]) + m = Planar() + arr = GO.NodedArrangement(m, A, B; exact = EX) + g = GO.OverlayGraph(m, arr; exact = EX) + @test any(i -> GO.is_collapse(g.edges[i].label, 0), eachindex(g.edges)) + check_graph(m, g) +end + +# --------------------------------------------------------------------------- +# 3. Label semantics — the four states, predicates, and the isForward swap +# --------------------------------------------------------------------------- + +@testset "the four label states from Edge sources" begin + # Boundary (area, nonzero delta) for A, NotPart for B + lbl = GO._create_label(GO.MergeEdge(Int32(1), Int32(2), Int32(1), Int32(1), + DIM_A, Int32(-1), false, DIM_NOT_PART, Int32(0), false)) + @test GO.is_boundary(lbl, 0) && GO.is_not_part(lbl, 1) + @test GO.is_boundary_either(lbl) && !GO.is_boundary_both(lbl) + @test GO.is_boundary_singleton(lbl) + @test GO.has_sides(lbl, 0) && !GO.has_sides(lbl, 1) + @test GO.get_line_location(lbl, 0) == LOC_INTERIOR # boundary line loc is INTERIOR + @test GO.is_known(lbl, 0) && !GO.is_known(lbl, 1) + + # Collapse (area, zero delta) + lbl = GO._create_label(GO.MergeEdge(Int32(1), Int32(2), Int32(1), Int32(1), + DIM_A, Int32(0), true, DIM_NOT_PART, Int32(0), false)) + @test GO.is_collapse(lbl, 0) && GO.is_linear(lbl, 0) + @test GO.is_line_location_unknown(lbl, 0) + GO.set_location_collapse!(lbl, 0) + @test GO.get_line_location(lbl, 0) == LOC_INTERIOR # hole collapse -> INTERIOR + + # Line + lbl = GO._create_label(GO.MergeEdge(Int32(1), Int32(2), Int32(1), Int32(1), + DIM_L, Int32(0), false, DIM_NOT_PART, Int32(0), false)) + @test GO.is_line(lbl) && GO.is_line(lbl, 0) && GO.is_linear(lbl, 0) + @test !GO.is_boundary(lbl, 0) + + # NotPart both + lbl = GO.OverlayLabel() + @test GO.is_not_part(lbl, 0) && GO.is_not_part(lbl, 1) + @test !GO.is_known(lbl, 0) && !GO.is_boundary_either(lbl) +end + +@testset "getLocation isForward L/R swap consistency" begin + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(2.0, 0.0), (4.0, 0.0), (4.0, 2.0), (2.0, 2.0), (2.0, 0.0)]]) + g = GO.OverlayGraph(Planar(), GO.NodedArrangement(Planar(), A, B; exact = EX); exact = EX) + for i in eachindex(g.edges) + s = GO.he_sym(g.edges, i) + lbl = g.edges[i].label + fi = g.edges[i].is_forward; fs = g.edges[s].is_forward + @test fi != fs # a pair has opposite directions + for idx in (0, 1) + #-- forward.L == reverse.R and forward.R == reverse.L, from ONE shared label + @test GO.get_location(lbl, idx, POS_LEFT, fi) == GO.get_location(lbl, idx, POS_RIGHT, fs) + @test GO.get_location(lbl, idx, POS_RIGHT, fi) == GO.get_location(lbl, idx, POS_LEFT, fs) + @test GO.get_location(lbl, idx, POS_ON, fi) == GO.get_location(lbl, idx, POS_ON, fs) + end + #-- the OverlayEdge accessor agrees with the label + orientation + @test GO.oe_get_location(g.edges, i, 0, POS_LEFT) == GO.get_location(lbl, 0, POS_LEFT, fi) + end +end + +@testset "result-marking flags and accessors" begin + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]]) + g = GO.OverlayGraph(Planar(), GO.NodedArrangement(Planar(), A, B; exact = EX); exact = EX) + i = 1; s = GO.he_sym(g.edges, i) + @test !GO.oe_in_result_area(g.edges, i) + GO.oe_mark_in_result_area!(g.edges, i) + @test GO.oe_in_result_area(g.edges, i) && !GO.oe_in_result_area(g.edges, s) + GO.oe_mark_in_result_area_both!(g.edges, i) + @test GO.oe_in_result_area_both(g.edges, i) + GO.oe_unmark_from_result_area_both!(g.edges, i) + @test !GO.oe_in_result_area(g.edges, i) && !GO.oe_in_result_area(g.edges, s) + GO.oe_set_next_result_max!(g.edges, i, s) + @test GO.oe_next_result_max(g.edges, i) == s && GO.oe_is_result_max_linked(g.edges, i) + GO.oe_mark_in_result_line!(g.edges, i) + @test GO.oe_in_result_line(g.edges, i) && GO.oe_in_result_line(g.edges, s) +end + +# --------------------------------------------------------------------------- +# 4. Graph invariants on small real inputs (env-gated, phase-1 smoke pattern) +# --------------------------------------------------------------------------- + +ne_ok = false +ne_names = String[]; ne_geoms = Any[] +try + import NaturalEarth, GeoJSON + fc = NaturalEarth.naturalearth("admin_0_countries", 110) + for f in fc + gg = GeoJSON.geometry(f) + (gg === nothing || GI.npoint(gg) == 0) && continue + nm = try; string(f.NAME); catch; "?"; end + push!(ne_names, nm); push!(ne_geoms, GO.tuples(gg)) + end + global ne_ok = length(ne_geoms) > 0 +catch err + @info "Natural Earth subset skipped (data unavailable)" err +end + +@testset "Natural Earth country-pair graph invariants" begin + if !ne_ok + @test_skip "Natural Earth data unavailable" + else + import LibGEOS as LG + picks = String["Brazil", "France", "Egypt", "India", "Australia"] + tested = 0 + for nm in picks + idx = findfirst(==(nm), ne_names) + idx === nothing && continue + A = ne_geoms[idx] + LG.isValid(GI.convert(LG, A)) || continue + B = GO.apply(GI.PointTrait(), A) do p + (GI.x(p) + 0.5, GI.y(p)) + end + tested += 1 + for m in (Planar(), Spherical()) + arr = GO.NodedArrangement(m, A, B; exact = EX) + g = GO.OverlayGraph(m, arr; exact = EX) + #-- every half-edge has a sym; stars are CCW-consistent; degrees + #-- match arrangement incidence (checks all invariants at once) + check_graph(m, g) + #-- each merged edge yields exactly one symmetric pair + @test iseven(length(g.edges)) + #-- node_edges representatives all originate at their node + for nid in 1:GO.num_nodes(arr) + e = GO.graph_node_edge(g, nid) + e == 0 || @test g.edges[e].origin == nid + end + end + end + @test tested >= 2 + end +end diff --git a/test/methods/clipping/overlayng/overlay_ng.jl b/test/methods/clipping/overlayng/overlay_ng.jl new file mode 100644 index 000000000..515f1b4f1 --- /dev/null +++ b/test/methods/clipping/overlayng/overlay_ng.jl @@ -0,0 +1,408 @@ +# Tests for the OverlayNG phase-2b engine core (design §3): the labeller, the +# result builders (polygons with the real minimal-ring split + hole nesting, +# lines, points), and the internal `_overlay_ng` driver — end to end, over the +# phase-1 arrangement and the phase-2a graph. +# +# Equality strategy: planar results are checked against LibGEOS's own overlay +# (high-level API only) with GEOS topological `equals` (order/orientation/merge- +# granularity independent), `isValid`, and area agreement at rtol 1e-12 (compared +# through the same `GO.area`, so it is a machine-precision gate). Spherical is +# gated on area conservation. Ported JTS cases keep JTS's expected WKT answers, +# compared via GEOS `equals`. + +using Test +import GeometryOps as GO +import GeoInterface as GI +import LibGEOS as LG +import GeometryOps: Planar, Spherical, True + +const EX = True() + +lgc(g) = GI.convert(LG, g) +giwkt(wkt) = GO.tuples(LG.readgeom(wkt)) + +const OPS = (GO.OVERLAY_INTERSECTION, GO.OVERLAY_UNION, + GO.OVERLAY_DIFFERENCE, GO.OVERLAY_SYMDIFFERENCE) +opname(op) = op == GO.OVERLAY_INTERSECTION ? "intersection" : + op == GO.OVERLAY_UNION ? "union" : + op == GO.OVERLAY_DIFFERENCE ? "difference" : "symdifference" + +geos_op(op, A, B) = + op == GO.OVERLAY_INTERSECTION ? LG.intersection(lgc(A), lgc(B)) : + op == GO.OVERLAY_UNION ? LG.union(lgc(A), lgc(B)) : + op == GO.OVERLAY_DIFFERENCE ? LG.difference(lgc(A), lgc(B)) : + LG.symmetricDifference(lgc(A), lgc(B)) + +# Planar: check `_overlay_ng` against LibGEOS for one op. +function check_planar(op, A, B; areatol = 1e-12) + r = GO._overlay_ng(Planar(), op, A, B; exact = EX) + geos = geos_op(op, A, B) + @test LG.isValid(lgc(r)) + @test LG.equals(lgc(r), geos) + @test isapprox(GO.area(Planar(), r), GO.area(Planar(), GO.tuples(geos)); + rtol = areatol, atol = 1e-12) + return r +end + +check_all_ops(A, B; areatol = 1e-12) = + for op in OPS + @testset "$(opname(op))" begin check_planar(op, A, B; areatol) end + end + +# The rounded coordinate multiset of a geometry (for vertex-set equality). +function coordset(g; digits = 9) + s = Set{Tuple{Float64, Float64}}() + for p in GI.getpoint(g) + push!(s, (round(GI.x(p); digits), round(GI.y(p); digits))) + end + return s +end + +# --------------------------------------------------------------------------- +# 1. S2/S3 case suite — all four ops vs LibGEOS (planar) +# --------------------------------------------------------------------------- + +@testset "overlapping squares (all ops)" begin + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]]) + #-- analytic areas: overlap [1,2]² = 1 + @test isapprox(GO.area(GO._overlay_ng(Planar(), GO.OVERLAY_INTERSECTION, A, B; exact = EX)), 1.0; rtol = 1e-12) + @test isapprox(GO.area(GO._overlay_ng(Planar(), GO.OVERLAY_UNION, A, B; exact = EX)), 7.0; rtol = 1e-12) + @test isapprox(GO.area(GO._overlay_ng(Planar(), GO.OVERLAY_DIFFERENCE, A, B; exact = EX)), 3.0; rtol = 1e-12) + @test isapprox(GO.area(GO._overlay_ng(Planar(), GO.OVERLAY_SYMDIFFERENCE, A, B; exact = EX)), 6.0; rtol = 1e-12) + check_all_ops(A, B) + #-- vertex-set equality against GEOS (intersection = the overlap square) + ri = GO._overlay_ng(Planar(), GO.OVERLAY_INTERSECTION, A, B; exact = EX) + @test coordset(ri) == coordset(GO.tuples(geos_op(GO.OVERLAY_INTERSECTION, A, B))) +end + +@testset "polygon-with-hole, B overlaps into the hole (§2.7 regression)" begin + #-- the wrong-area-hole case: B reaches into A's hole. The material-interior + #-- authority (§2.7) must give the hole the right side, or the areas invert. + A = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)], + [(3.0, 3.0), (7.0, 3.0), (7.0, 7.0), (3.0, 7.0), (3.0, 3.0)]]) + B = GI.Polygon([[(5.0, 5.0), (12.0, 5.0), (12.0, 12.0), (5.0, 12.0), (5.0, 5.0)]]) + #-- A area = 100 - 16 = 84; overlap of B with A-material = B∩A = 21 + @test isapprox(GO.area(GO._overlay_ng(Planar(), GO.OVERLAY_INTERSECTION, A, B; exact = EX)), 21.0; rtol = 1e-12) + check_all_ops(A, B) +end + +@testset "collinear shared boundary (degenerate intersection, merged union)" begin + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + B = GI.Polygon([[(2.0, 0.0), (4.0, 0.0), (4.0, 2.0), (2.0, 2.0), (2.0, 0.0)]]) + #-- intersection is the shared boundary line (1-D), not an area + ri = GO._overlay_ng(Planar(), GO.OVERLAY_INTERSECTION, A, B; exact = EX) + @test GI.trait(ri) isa GI.LineStringTrait + @test LG.equals(lgc(ri), geos_op(GO.OVERLAY_INTERSECTION, A, B)) + #-- union merges into one 2×4 box + ru = GO._overlay_ng(Planar(), GO.OVERLAY_UNION, A, B; exact = EX) + @test GI.trait(ru) isa GI.PolygonTrait + @test isapprox(GO.area(ru), 8.0; rtol = 1e-12) + @test LG.isValid(lgc(ru)) && LG.equals(lgc(ru), geos_op(GO.OVERLAY_UNION, A, B)) +end + +@testset "degree-6 coincident crossing (all ops)" begin + #-- two A squares touching at (1,1) + a B triangle edge through (1,1) + A = GI.MultiPolygon([[[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]], + [[(1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1.0, 2.0), (1.0, 1.0)]]]) + B = GI.Polygon([[(0.0, 0.0), (2.0, 2.0), (2.0, 0.0), (0.0, 0.0)]]) + arr = GO.NodedArrangement(Planar(), A, B; exact = EX) + g = GO.OverlayGraph(Planar(), arr; exact = EX) + nid = findfirst(1:GO.num_nodes(arr)) do i + p = GO.node_point(arr, i) + isapprox(p[1], 1.0; atol = 1e-9) && isapprox(p[2], 1.0; atol = 1e-9) + end + @test nid !== nothing + e = GO.graph_node_edge(g, nid) + @test e != 0 + @test GO.he_degree(g.edges, e) == 6 # 2 (square 1) + 2 (square 2) + 2 (B diagonal) + check_all_ops(A, B) +end + +@testset "concave L-shapes (all ops)" begin + A = GI.Polygon([[(0.0, 0.0), (3.0, 0.0), (3.0, 1.0), (1.0, 1.0), (1.0, 3.0), (0.0, 3.0), (0.0, 0.0)]]) + B = GI.Polygon([[(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (2.0, 3.0), (2.0, 1.0), (0.0, 1.0), (0.0, 0.0)]]) + check_all_ops(A, B) +end + +@testset "MultiPolygon input (all ops)" begin + A = GI.MultiPolygon([[[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]], + [[(5.0, 5.0), (7.0, 5.0), (7.0, 7.0), (5.0, 7.0), (5.0, 5.0)]]]) + B = GI.Polygon([[(1.0, 1.0), (6.0, 1.0), (6.0, 6.0), (1.0, 6.0), (1.0, 1.0)]]) + check_all_ops(A, B) +end + +# --------------------------------------------------------------------------- +# 2. Ported JTS OverlayNGTest subset (floating-safe area/line cases) +# --------------------------------------------------------------------------- +# +# Equality: GEOS `equals` against JTS's expected WKT (topological — handles +# JTS's coordinate ordering and OverlayNG's line-merging, which this engine +# emits as raw noded segments). +# +# SKIPPED, with reasons (all fixed-precision / topology-collapse tests that do +# not apply to a floating-only, non-snapping engine — design §0): +# testTriangleFillingHoleUnion(Prec10), testBoxTri{Intersection,Union}, +# test2spikes{Intersection,Union}, testTriBoxIntersection, +# testCollapse* / testSnapBoxGore* (topology collapse from precision rounding), +# testVerySmallBIntersection (scale 1e8), testEdgeDisappears (scale 1e6), +# testBcollapse* / testBNearVertexSnappingCausesInversion / +# testBCollapsedHoleEdgeLabelledExterior (snap-rounding collapse), +# testDisjointLinesRoundedIntersection (coordinate rounding to a point). +# SKIPPED (substrate limitation, not precision): testTouchingPolyDifference — a +# single input whose hole touches its own shell, where the DIFFERENCE splits the +# result into point-touching polygons; the substrate does not self-node one +# input (design §2.2), so it yields a correct-area but non-simple result. + +const JTS_CASES = [ + # (name, op, A_wkt, B_wkt, expected_wkt) + ("NestedShellsIntersection", GO.OVERLAY_INTERSECTION, + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", + "POLYGON ((120 180, 180 180, 180 120, 120 120, 120 180))", + "POLYGON ((120 180, 180 180, 180 120, 120 120, 120 180))"), + ("NestedShellsUnion", GO.OVERLAY_UNION, + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", + "POLYGON ((120 180, 180 180, 180 120, 120 120, 120 180))", + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))"), + ("AdjacentBoxesIntersection", GO.OVERLAY_INTERSECTION, + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", + "POLYGON ((300 200, 300 100, 200 100, 200 200, 300 200))", + "LINESTRING (200 100, 200 200)"), + ("AdjacentBoxesUnion", GO.OVERLAY_UNION, + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", + "POLYGON ((300 200, 300 100, 200 100, 200 200, 300 200))", + "POLYGON ((100 100, 100 200, 200 200, 300 200, 300 100, 200 100, 100 100))"), + ("TouchingHoleUnion", GO.OVERLAY_UNION, + "POLYGON ((100 300, 300 300, 300 100, 100 100, 100 300), (200 200, 150 200, 200 300, 200 200))", + "POLYGON ((130 160, 260 160, 260 120, 130 120, 130 160))", + "POLYGON ((100 100, 100 300, 200 300, 300 300, 300 100, 100 100), (150 200, 200 200, 200 300, 150 200))"), + ("TouchingMultiHoleUnion", GO.OVERLAY_UNION, + "POLYGON ((100 300, 300 300, 300 100, 100 100, 100 300), (200 200, 150 200, 200 300, 200 200), (250 230, 216 236, 250 300, 250 230), (235 198, 300 200, 237 175, 235 198))", + "POLYGON ((130 160, 260 160, 260 120, 130 120, 130 160))", + "POLYGON ((100 300, 200 300, 250 300, 300 300, 300 200, 300 100, 100 100, 100 300), (200 300, 150 200, 200 200, 200 300), (250 300, 216 236, 250 230, 250 300), (300 200, 235 198, 237 175, 300 200))"), + ("ATouchingNestedPolyUnion", GO.OVERLAY_UNION, + "MULTIPOLYGON (((0 200, 200 200, 200 0, 0 0, 0 200), (50 50, 190 50, 50 200, 50 50)), ((60 100, 100 60, 50 50, 60 100)))", + "POLYGON ((135 176, 180 176, 180 130, 135 130, 135 176))", + "MULTIPOLYGON (((0 0, 0 200, 50 200, 200 200, 200 0, 0 0), (50 50, 190 50, 50 200, 50 50)), ((50 50, 60 100, 100 60, 50 50)))"), + ("BoxLineIntersection", GO.OVERLAY_INTERSECTION, + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", + "LINESTRING (50 150, 150 150)", + "LINESTRING (100 150, 150 150)"), + ("BoxLineUnion", GO.OVERLAY_UNION, + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", + "LINESTRING (50 150, 150 150)", + "GEOMETRYCOLLECTION (POLYGON ((200 200, 200 100, 100 100, 100 150, 100 200, 200 200)), LINESTRING (50 150, 100 150))"), + ("LinePolygonUnion", GO.OVERLAY_UNION, + "LINESTRING (50 150, 150 150)", + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", + "GEOMETRYCOLLECTION (LINESTRING (50 150, 100 150), POLYGON ((100 200, 200 200, 200 100, 100 100, 100 150, 100 200)))"), + ("LinePolygonUnionAlongPolyBoundary", GO.OVERLAY_UNION, + "LINESTRING (150 300, 250 300)", + "POLYGON ((100 400, 200 400, 200 300, 100 300, 100 400))", + "GEOMETRYCOLLECTION (LINESTRING (200 300, 250 300), POLYGON ((200 300, 150 300, 100 300, 100 400, 200 400, 200 300)))"), + ("LinePolygonIntersectionAlongPolyBoundary", GO.OVERLAY_INTERSECTION, + "LINESTRING (150 300, 250 300)", + "POLYGON ((100 400, 200 400, 200 300, 100 300, 100 400))", + "LINESTRING (200 300, 150 300)"), + ("PolygonLineVerticalIntersection", GO.OVERLAY_INTERSECTION, + "POLYGON ((-200 -200, 200 -200, 200 200, -200 200, -200 -200))", + "LINESTRING (-100 100, -100 -100)", + "LINESTRING (-100 100, -100 -100)"), + ("PolygonLineHorizontalIntersection", GO.OVERLAY_INTERSECTION, + "POLYGON ((10 90, 90 90, 90 10, 10 10, 10 90))", + "LINESTRING (20 50, 80 50)", + "LINESTRING (20 50, 80 50)"), + ("PolygonMultiLineUnion", GO.OVERLAY_UNION, + "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", + "MULTILINESTRING ((150 250, 150 50), (250 250, 250 50))", + "GEOMETRYCOLLECTION (LINESTRING (150 50, 150 100), LINESTRING (150 200, 150 250), LINESTRING (250 50, 250 250), POLYGON ((100 100, 100 200, 150 200, 200 200, 200 100, 150 100, 100 100)))"), + ("PolygonLineIntersectionOrder", GO.OVERLAY_INTERSECTION, + "POLYGON ((1 1, 1 9, 9 9, 9 7, 3 7, 3 3, 9 3, 9 1, 1 1))", + "MULTILINESTRING ((2 10, 2 0), (4 10, 4 0))", + "MULTILINESTRING ((2 9, 2 1), (4 9, 4 7), (4 3, 4 1))"), + ("AreaLineIntersection", GO.OVERLAY_INTERSECTION, + "POLYGON ((360 200, 220 200, 220 180, 300 180, 300 160, 300 140, 360 200))", + "MULTIPOLYGON (((280 180, 280 160, 300 160, 300 180, 280 180)), ((220 230, 240 230, 240 180, 220 180, 220 230)))", + "GEOMETRYCOLLECTION (LINESTRING (280 180, 300 180), LINESTRING (300 160, 300 180), POLYGON ((220 180, 220 200, 240 200, 240 180, 220 180)))"), + ("LineUnion", GO.OVERLAY_UNION, + "LINESTRING (0 0, 1 1)", "LINESTRING (1 1, 2 2)", + "MULTILINESTRING ((0 0, 1 1), (1 1, 2 2))"), + ("Line2Union", GO.OVERLAY_UNION, + "LINESTRING (0 0, 1 1, 0 1)", "LINESTRING (1 1, 2 2, 3 3)", + "MULTILINESTRING ((0 0, 1 1), (0 1, 1 1), (1 1, 2 2, 3 3))"), + ("Line3Union", GO.OVERLAY_UNION, + "MULTILINESTRING ((0 1, 1 1), (2 2, 2 0))", "LINESTRING (0 0, 1 1, 2 2, 3 3)", + "MULTILINESTRING ((0 0, 1 1), (0 1, 1 1), (1 1, 2 2), (2 0, 2 2), (2 2, 3 3))"), + ("Line4Union", GO.OVERLAY_UNION, + "LINESTRING (100 300, 200 300, 200 100, 100 100)", + "LINESTRING (300 300, 200 300, 200 300, 200 100, 300 100)", + "MULTILINESTRING ((200 100, 100 100), (300 300, 200 300), (200 300, 200 100), (200 100, 300 100), (100 300, 200 300))"), + ("LineFigure8Union", GO.OVERLAY_UNION, + "LINESTRING (5 1, 2 2, 5 3, 2 4, 5 5)", "LINESTRING (5 1, 8 2, 5 3, 8 4, 5 5)", + "MULTILINESTRING ((5 1, 2 2, 5 3), (5 1, 8 2, 5 3), (5 3, 2 4, 5 5), (5 3, 8 4, 5 5))"), + ("LineRingUnion", GO.OVERLAY_UNION, + "LINESTRING (1 1, 5 5, 9 1)", "LINESTRING (1 1, 9 1)", + "MULTILINESTRING ((1 1, 5 5, 9 1), (1 1, 9 1))"), + ("PolygonFlatCollapseIntersection", GO.OVERLAY_INTERSECTION, + "POLYGON ((200 100, 150 200, 250 200, 150 200, 100 100, 200 100))", + "POLYGON ((50 150, 250 150, 250 50, 50 50, 50 150))", + "POLYGON ((175 150, 200 100, 100 100, 125 150, 175 150))"), +] + +@testset "ported JTS OverlayNGTest subset ($(length(JTS_CASES)) cases)" begin + for (name, op, awkt, bwkt, ewkt) in JTS_CASES + @testset "$name" begin + r = GO._overlay_ng(Planar(), op, giwkt(awkt), giwkt(bwkt); exact = EX) + @test LG.isValid(lgc(r)) + @test LG.equals(lgc(r), LG.readgeom(ewkt)) + end + end +end + +# --------------------------------------------------------------------------- +# 3. Ring-builder specifics +# --------------------------------------------------------------------------- + +@testset "self-touching result ring (minimal-ring split)" begin + #-- A minus a triangle B that touches A's boundary at a single point (0,3): + #-- the result ring self-touches there and must split into shell + hole. + A = GI.Polygon([[(0.0, 0.0), (6.0, 0.0), (6.0, 6.0), (0.0, 6.0), (0.0, 0.0)]]) + B = GI.Polygon([[(0.0, 3.0), (4.0, 1.0), (4.0, 5.0), (0.0, 3.0)]]) + r = GO._overlay_ng(Planar(), GO.OVERLAY_DIFFERENCE, A, B; exact = EX) + @test GI.trait(r) isa GI.PolygonTrait + @test GI.nring(r) == 2 # shell + one split-out hole + @test LG.isValid(lgc(r)) + @test LG.equals(lgc(r), geos_op(GO.OVERLAY_DIFFERENCE, A, B)) + @test isapprox(GO.area(r), 28.0; rtol = 1e-12) # 36 - 8 +end + +@testset "free-hole assignment (strictly interior hole)" begin + #-- difference of a strictly-interior square from a shell → one free hole + Big = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]]) + Inner = GI.Polygon([[(3.0, 3.0), (7.0, 3.0), (7.0, 7.0), (3.0, 7.0), (3.0, 3.0)]]) + r = GO._overlay_ng(Planar(), GO.OVERLAY_DIFFERENCE, Big, Inner; exact = EX) + @test GI.trait(r) isa GI.PolygonTrait + @test GI.nring(r) == 2 + @test LG.isValid(lgc(r)) + @test isapprox(GO.area(r), 84.0; rtol = 1e-12) +end + +@testset "union of a multi-island geometry (France-class nesting)" begin + #-- 20 disjoint island squares, unioned with a shifted copy (overlaps) — the + #-- many-component case the spike prototypes faked, exercising minimal-ring + #-- split + hole nesting across many shells. + isl = Vector{Vector{Vector{Tuple{Float64, Float64}}}}() + for i in 0:19 + x = (i % 5) * 3.0; y = (i ÷ 5) * 3.0 + push!(isl, [[(x, y), (x + 2, y), (x + 2, y + 2), (x, y + 2), (x, y)]]) + end + MI = GI.MultiPolygon(isl) + #-- shift by (1,1) so islands overlap their diagonal neighbours + MI2 = GO.apply(GI.PointTrait(), MI) do p + (GI.x(p) + 1.0, GI.y(p) + 1.0) + end + r = GO._overlay_ng(Planar(), GO.OVERLAY_UNION, MI, MI2; exact = EX) + geos = geos_op(GO.OVERLAY_UNION, MI, MI2) + @test LG.isValid(lgc(r)) + @test LG.equals(lgc(r), geos) + @test isapprox(GO.area(r), GO.area(GO.tuples(geos)); rtol = 1e-12) +end + +# --------------------------------------------------------------------------- +# 4. Spherical end-to-end +# --------------------------------------------------------------------------- + +@testset "spherical overlapping quads — area conservation" begin + A = GI.Polygon([[(0.0, 0.0), (20.0, 0.0), (20.0, 20.0), (0.0, 20.0), (0.0, 0.0)]]) + B = GI.Polygon([[(10.0, 10.0), (30.0, 10.0), (30.0, 30.0), (10.0, 30.0), (10.0, 10.0)]]) + ai = GO.area(Spherical(), GO._overlay_ng(Spherical(), GO.OVERLAY_INTERSECTION, A, B; exact = EX)) + au = GO.area(Spherical(), GO._overlay_ng(Spherical(), GO.OVERLAY_UNION, A, B; exact = EX)) + aA = GO.area(Spherical(), A); aB = GO.area(Spherical(), B) + @test isapprox(au + ai, aA + aB; rtol = 1e-12) + #-- difference + intersection reconstruct A + ad = GO.area(Spherical(), GO._overlay_ng(Spherical(), GO.OVERLAY_DIFFERENCE, A, B; exact = EX)) + @test isapprox(ad + ai, aA; rtol = 1e-12) + #-- symdifference = union - intersection + asd = GO.area(Spherical(), GO._overlay_ng(Spherical(), GO.OVERLAY_SYMDIFFERENCE, A, B; exact = EX)) + @test isapprox(asd, au - ai; rtol = 1e-12) +end + +@testset "spherical empty-vs-full disambiguation (§3 amendment 6)" begin + A = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]]) + Bdisjoint = GI.Polygon([[(40.0, 0.0), (50.0, 0.0), (50.0, 10.0), (40.0, 10.0), (40.0, 0.0)]]) + #-- disjoint intersection → empty (NOT full-sphere) + ri = GO._overlay_ng(Spherical(), GO.OVERLAY_INTERSECTION, A, Bdisjoint; exact = EX) + @test GI.npoint(ri) == 0 + @test isapprox(GO.area(Spherical(), ri), 0.0; atol = 1e-9) + + #-- the disambiguation function directly: a boundaryless union that covers + #-- everything is the full sphere and must throw; an empty intersection must not. + inp_union = GO._OverlayInput(Spherical(), A, A, 2, 2, EX, false, false, nothing, nothing) + @test GO._covers_everything(Spherical(), GO.OVERLAY_UNION, inp_union) + @test_throws ArgumentError GO._resolve_empty_result(Spherical(), GO.OVERLAY_UNION, inp_union) + + inp_int = GO._OverlayInput(Spherical(), A, Bdisjoint, 2, 2, EX, false, false, nothing, nothing) + @test !GO._covers_everything(Spherical(), GO.OVERLAY_INTERSECTION, inp_int) + rr = GO._resolve_empty_result(Spherical(), GO.OVERLAY_INTERSECTION, inp_int) + @test GI.npoint(rr) == 0 +end + +# --------------------------------------------------------------------------- +# 5. Input validation + empty inputs +# --------------------------------------------------------------------------- + +@testset "input validation and empty short-circuits" begin + A = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]]) + #-- point inputs are rejected (phase 3) + @test_throws ArgumentError GO._overlay_ng(Planar(), GO.OVERLAY_INTERSECTION, GI.Point((1.0, 1.0)), A; exact = EX) + #-- disjoint intersection short-circuits to empty (planar envelope) + Far = GI.Polygon([[(100.0, 100.0), (102.0, 100.0), (102.0, 102.0), (100.0, 102.0), (100.0, 100.0)]]) + r = GO._overlay_ng(Planar(), GO.OVERLAY_INTERSECTION, A, Far; exact = EX) + @test GI.npoint(r) == 0 +end + +# --------------------------------------------------------------------------- +# 6. Spherical NE shifted-self smoke (env-gated, phase-1 smoke pattern) +# --------------------------------------------------------------------------- + +ne_ok = false +ne_names = String[]; ne_geoms = Any[] +try + import NaturalEarth, GeoJSON + fc = NaturalEarth.naturalearth("admin_0_countries", 110) + for f in fc + gg = GeoJSON.geometry(f) + (gg === nothing || GI.npoint(gg) == 0) && continue + nm = try; string(f.NAME); catch; "?"; end + push!(ne_names, nm); push!(ne_geoms, GO.tuples(gg)) + end + global ne_ok = length(ne_geoms) > 0 +catch err + @info "Natural Earth subset skipped (data unavailable)" err +end + +@testset "Natural Earth shifted-self area conservation (spherical + planar)" begin + if !ne_ok + @test_skip "Natural Earth data unavailable" + else + picks = String["Brazil", "France", "Egypt", "Australia"] + tested = 0 + for nm in picks + idx = findfirst(==(nm), ne_names) + idx === nothing && continue + A = ne_geoms[idx] + LG.isValid(lgc(A)) || continue + B = GO.apply(GI.PointTrait(), A) do p + (GI.x(p) + 0.5, GI.y(p)) + end + tested += 1 + for m in (Planar(), Spherical()) + ri = GO._overlay_ng(m, GO.OVERLAY_INTERSECTION, A, B; exact = EX) + ru = GO._overlay_ng(m, GO.OVERLAY_UNION, A, B; exact = EX) + aA = GO.area(m, A); aB = GO.area(m, B) + @test isapprox(GO.area(m, ru) + GO.area(m, ri), aA + aB; rtol = 1e-9) + end + end + @test tested >= 2 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index ec62ff8c8..e1c4f1b84 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,6 +45,8 @@ end @safetestset "Polygon Clipping" begin include("methods/clipping/polygon_clipping.jl") end @safetestset "Sutherland-Hodgman" begin include("methods/clipping/sutherland_hodgman.jl") end @safetestset "OverlayNG noding" begin include("methods/clipping/overlayng/noding.jl") end +@safetestset "OverlayNG graph" begin include("methods/clipping/overlayng/overlay_graph.jl") end +@safetestset "OverlayNG engine" begin include("methods/clipping/overlayng/overlay_ng.jl") end # Transformations @safetestset "Embed Extent" begin include("transformations/extent.jl") end @safetestset "Reproject" begin include("transformations/reproject.jl") end