diff --git a/src/methods/geom_relations/relateng/edge_intersector.jl b/src/methods/geom_relations/relateng/edge_intersector.jl index 03ed7a644..96b0f8106 100644 --- a/src/methods/geom_relations/relateng/edge_intersector.jl +++ b/src/methods/geom_relations/relateng/edge_intersector.jl @@ -180,15 +180,15 @@ recorded on `computer`. - `NestedLoop`: a plain double loop over string pairs and segment pairs, with a per-pair segment-extent disjointness skip (on `Planar`). -- Any tree-backed accelerator (e.g. `DoubleSTRtree`): a spatial index - (`_relate_edge_index`, currently a `NaturalIndex`) is built over the - per-segment extents of each side and traversed with +- Any tree-backed accelerator (canonically `DoubleNaturalTree`): a spatial + index (`_relate_edge_index`, currently a `NaturalIndex`) is built over + the per-segment extents of each side and traversed with `SpatialTreeInterface.dual_depth_first_search` under the `Extents.intersects` predicate. - [`AutoAccelerator`](@ref): picks `NestedLoop` below the clipping size threshold (`GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS`) and on - non-`Planar` manifolds (planar extent trees are not valid there), and the - tree path otherwise. + manifolds without a segment-extent kernel (neither `Planar` nor + `Spherical`), and the tree path otherwise. After each processed pair `is_result_known(computer)` is consulted and the enumeration stops early once the predicate value is determined (the port of @@ -220,29 +220,20 @@ function process_edge_intersections!(tc::TopologyComputer, _select_edge_set_accelerator(m, ssa_list, ssb_list); m, exact) end -# STRtrees over planar extents are only valid on the Planar manifold, and -# below the clipping threshold the nested loop wins anyway. -function _select_edge_set_accelerator(::Planar, ssa_list, ssb_list) +# Below the clipping threshold the nested loop wins; above it, the tree +# path. Valid on both manifolds with a segment-extent kernel: planar boxes +# on `Planar`, 3D great-circle arc extents (`_segment_extent`) on +# `Spherical` — `NaturalIndex`, the dual DFS, and `Extents.intersects` are +# dimension-generic. Other manifolds have no extent kernel and always take +# the nested loop. +function _select_edge_set_accelerator(::Union{Planar, Spherical}, ssa_list, ssb_list) na = _total_segment_count(ssa_list) nb = _total_segment_count(ssb_list) if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS && nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS return NestedLoop() else - return DoubleSTRtree() - end -end -# The same threshold heuristic on the sphere: the tree path is valid because -# the segment extents are 3D great-circle arc extents (`_segment_extent`), and -# `NaturalIndex` / the dual DFS / `Extents.intersects` are dimension-generic. -function _select_edge_set_accelerator(::Spherical, ssa_list, ssb_list) - na = _total_segment_count(ssa_list) - nb = _total_segment_count(ssb_list) - if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS && - nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS - return NestedLoop() - else - return DoubleSTRtree() + return DoubleNaturalTree() end end _select_edge_set_accelerator(::Manifold, ssa_list, ssb_list) = NestedLoop() @@ -300,7 +291,7 @@ _segment_envs_disjoint(::Manifold, a0, a1, b0, b1) = false _relate_edge_index(extents::Vector{<:Extents.Extent}) = NaturalIndex(extents; nodecapacity = 16) -# Tree path (any other accelerator, canonically DoubleSTRtree): a spatial +# Tree path (any other accelerator, canonically DoubleNaturalTree): a spatial # index over the per-segment extents of each side, traversed simultaneously. function process_edge_intersections!(tc::TopologyComputer, ssa_list::AbstractVector{<:RelateSegmentString}, diff --git a/src/methods/geom_relations/relateng/indexed_point_in_area.jl b/src/methods/geom_relations/relateng/indexed_point_in_area.jl index b135c6e8d..a8c6d45f3 100644 --- a/src/methods/geom_relations/relateng/indexed_point_in_area.jl +++ b/src/methods/geom_relations/relateng/indexed_point_in_area.jl @@ -48,20 +48,15 @@ Port of JTS `SortedPackedIntervalRTree`, with two representation changes of level `k` (an unpaired trailing node is carried up unchanged, as in `buildLevel`). The last level is the root. - JTS always sorts the leaves by interval midpoint (`NodeComparator`) - before packing. Here `sort_leaves = false` skips that and keeps insertion - (ring) order — the `NaturalIndexing` observation. Query results are - order-independent (every visit is extent-checked), but the layouts trade - off differently: midpoint order groups same-`y` segments so a point query - descends few subtrees, while a long coastline ring in natural order - recrosses the query `y` in many separated runs. Measured on Natural Earth - 10m Canada: the sort is ~4× the rest of the build, and natural-order - queries are ~3× slower. So prepared mode sorts (build once, query - forever) and the lazily indexed unprepared path doesn't (its query count - is at most a few hundred, far below the ~1000-query crossover; see - `locate_on_polygonal`). + before packing; so does this port. The sort earns its cost in the + index's only (prepared, build-once-query-forever) use: midpoint order + groups same-`y` segments so a point query descends few subtrees, where + insertion (ring) order recrosses the query `y` in many separated runs — + measured on Natural Earth 10m Canada, the sort is ~4× the rest of the + build and ring-order queries are ~3× slower. """ struct SortedPackedIntervalRTree{I} - # leaf items: midpoint-sorted (`sort_leaves = true`) or insertion order + # leaf items, midpoint-sorted items::Vector{I} # level_min[1][i] / level_max[1][i] is the interval of leaf item i; # level k > 1 holds the pairwise-combined extents of level k - 1 @@ -70,19 +65,15 @@ struct SortedPackedIntervalRTree{I} end # Port of insert + init/buildRoot/buildTree/buildLevel, packed eagerly. -# Without `sort_leaves` the leaf arrays are taken over by the tree, not -# copied. function SortedPackedIntervalRTree(mins::Vector{Float64}, maxs::Vector{Float64}, - items::Vector{I}; sort_leaves::Bool = true) where {I} - if sort_leaves - #-- sort the leaf nodes (IntervalRTreeNode.NodeComparator: by - #-- midpoint; sortperm is stable, matching Collections.sort) - n = length(items) - perm = sortperm(Float64[(mins[i] + maxs[i]) / 2 for i in 1:n]) - mins = mins[perm] - maxs = maxs[perm] - items = items[perm] - end + items::Vector{I}) where {I} + #-- sort the leaf nodes (IntervalRTreeNode.NodeComparator: by + #-- midpoint; sortperm is stable, matching Collections.sort) + n = length(items) + perm = sortperm(Float64[(mins[i] + maxs[i]) / 2 for i in 1:n]) + mins = mins[perm] + maxs = maxs[perm] + items = items[perm] level_min = [mins] level_max = [maxs] #-- now group nodes into blocks of two and build tree up recursively @@ -285,14 +276,14 @@ struct IndexedPointInAreaLocator{M <: Manifold, E} is_empty::Bool end -function IndexedPointInAreaLocator(m::Manifold, geom; exact, sort_leaves::Bool = true) +function IndexedPointInAreaLocator(m::Manifold, geom; exact) mins = Float64[] maxs = Float64[] segs = _PIASegment[] n = GI.npoint(geom) sizehint!(mins, n); sizehint!(maxs, n); sizehint!(segs, n) _interval_index_add_geom!(mins, maxs, segs, GI.trait(geom), geom) - index = SortedPackedIntervalRTree(mins, maxs, segs; sort_leaves) + index = SortedPackedIntervalRTree(mins, maxs, segs) #-- IntervalIndexedGeometry.isEmpty: a (recursively) empty polygonal #-- geometry contributes no rings, hence no segments return IndexedPointInAreaLocator(m, exact, index, isempty(segs)) diff --git a/src/methods/geom_relations/relateng/point_locator.jl b/src/methods/geom_relations/relateng/point_locator.jl index 7c1ef96ed..b2b6a3095 100644 --- a/src/methods/geom_relations/relateng/point_locator.jl +++ b/src/methods/geom_relations/relateng/point_locator.jl @@ -235,12 +235,9 @@ bnRule)`; the manifold/`exact` parameters are the only additions (consistent with [`AdjacentEdgeLocator`](@ref)). As in JTS, prepared mode swaps the per-polygon `SimplePointInAreaLocator` ring loop for a cached [`IndexedPointInAreaLocator`](@ref) (indexed_point_in_area.jl), created -lazily on the first use per polygonal element (Task 22). Unprepared mode -deviates from Java (which keys indexing on `isPrepared` alone): the first -query on a polygonal element uses the direct ring loop, but repeat queries -build and reuse the indexed locator — one O(n) scan beats an O(n) index -build, while the many area-vertex locations of a multi-element relate -amortize the index (see `locate_on_polygonal`). +lazily on the first use per polygonal element (Task 22); unprepared mode +scans the rings directly on every query. Repeated point location against +one geometry is what [`prepare`](@ref) is for. """ mutable struct RelatePointLocator{M <: Manifold, E, G, BR <: BoundaryNodeRule, P} const m::M @@ -257,14 +254,10 @@ mutable struct RelatePointLocator{M <: Manifold, E, G, BR <: BoundaryNodeRule, P const polygons::Vector{Any} const line_boundary::LinearBoundary{BR, P} const is_empty::Bool - # per-polygonal-element indexed locators, created lazily by - # `_get_poly_locator` (Java: polyLocator, filled by getLocator). - # Prepared mode fills an entry on its first query; unprepared mode on - # its second (see `locate_on_polygonal`). + # per-polygonal-element indexed locators (prepared mode only), created + # lazily by `_get_poly_locator` on the element's first query + # (Java: polyLocator, filled by getLocator) const poly_locator::Vector{Union{Nothing, IndexedPointInAreaLocator{M, E}}} - # unprepared mode: queries seen per polygonal element, driving the lazy - # index heuristic above - const poly_query_count::Vector{Int32} # lazily built on the first multi-boundary point (Java: adjEdgeLocator) adj_edge_locator::Union{Nothing, AdjacentEdgeLocator{M, E, P}} end @@ -285,19 +278,18 @@ function RelatePointLocator(m::Manifold, geom; exact, # LinearBoundary behaves identically (no boundary, no boundary points), # so it is built unconditionally here. line_boundary = LinearBoundary(m, lines, boundary_rule) - # Java allocates `polyLocator` for both modes (Simple/Indexed); here both - # modes may cache indexed locator objects (unprepared lazily, on repeat - # queries), so it is allocated unconditionally. + # Java allocates `polyLocator` for both modes (its unprepared arm caches + # SimplePointInAreaLocator objects); the direct ring scan here is + # stateless, so only prepared mode fills it. poly_locator = Vector{Union{Nothing, IndexedPointInAreaLocator{typeof(m), typeof(exact)}}}( nothing, length(polygons)) - poly_query_count = zeros(Int32, length(polygons)) #-- P cannot be inferred from the `nothing` adj_edge_locator, so spell out #-- every type parameter return RelatePointLocator{typeof(m), typeof(exact), typeof(geom), typeof(boundary_rule), P}( m, exact, geom, is_prepared, boundary_rule, points, lines, polygons, line_boundary, is_empty, poly_locator, - poly_query_count, nothing) + nothing) end has_boundary(loc::RelatePointLocator) = has_boundary(loc.line_boundary) @@ -501,24 +493,9 @@ function locate_on_polygons(loc::RelatePointLocator, p, is_node::Bool, parent_po return LOC_EXTERIOR end -# Queries a polygonal element absorbs via the direct ring loop before its -# IndexedPointInAreaLocator is built. Both costs scale with the element's -# segment count, so one threshold fits all sizes: an unsorted index build -# costs ~10-13 ring scans (measured on Natural Earth coastlines), making -# the worst-case regret of switching at 8 about one build. Real relates are -# bimodal — a handful of queries (barely-touching neighbors, where indexing -# never pays) or hundreds (one area-vertex location per polygon element of -# the other geometry), so the threshold rarely sits near the break-even. -const _LAZY_INDEX_QUERY_THRESHOLD = Int32(8) - # Port of RelatePointLocator.locateOnPolygonal: Java dispatches to a # per-polygonal PointOnGeometryLocator — a cached IndexedPointInAreaLocator -# when prepared, a SimplePointInAreaLocator otherwise. Prepared mode does -# the same here (Task 22). Unprepared mode deviates from Java: the first -# query on an element uses the direct SimplePointInAreaLocator ring loop -# (one O(n) scan beats an O(n) index build + query), but repeat queries — -# e.g. one area-vertex location per polygon element of the other geometry -# in a multipolygon/multipolygon relate — build and amortize the index. +# when prepared, a SimplePointInAreaLocator otherwise. Same here (Task 22). function locate_on_polygonal(loc::RelatePointLocator, p, is_node::Bool, parent_polygonal, index::Int) polygonal = loc.polygons[index] if is_node && parent_polygonal === polygonal @@ -527,29 +504,18 @@ function locate_on_polygonal(loc::RelatePointLocator, p, is_node::Bool, parent_p #-- the RayCrossingCounter horizontal-ray sweep is coordinate-plane #-- logic (as is all of JTS), so a future non-planar kernel falls #-- through to its own rk_point_in_ring even when prepared - if loc.m isa Planar - use_index = loc.is_prepared - if !use_index - count = (loc.poly_query_count[index] += Int32(1)) - use_index = count > _LAZY_INDEX_QUERY_THRESHOLD - end - if use_index - return locate(_get_poly_locator(loc, index), p) - end + if loc.is_prepared && loc.m isa Planar + return locate(_get_poly_locator(loc, index), p) end return _locate_point_in_polygonal(loc.m, p, GI.trait(polygonal), polygonal; exact = loc.exact) end # Port of RelatePointLocator.getLocator (indexed arm): lazily create and -# cache the indexed locator for polygonal element `index`. Prepared mode -# pays for the midpoint-sorted layout (build once, query forever); the -# unprepared lazy index skips the sort, which dominates the build cost -# (see `SortedPackedIntervalRTree`). +# cache the indexed locator for polygonal element `index`. function _get_poly_locator(loc::RelatePointLocator, index::Int) locator = loc.poly_locator[index] if locator === nothing - locator = IndexedPointInAreaLocator(loc.m, loc.polygons[index]; - exact = loc.exact, sort_leaves = loc.is_prepared) + locator = IndexedPointInAreaLocator(loc.m, loc.polygons[index]; exact = loc.exact) loc.poly_locator[index] = locator end return locator diff --git a/src/methods/geom_relations/relateng/relate_ng.jl b/src/methods/geom_relations/relateng/relate_ng.jl index 7f3de205c..6034d28b7 100644 --- a/src/methods/geom_relations/relateng/relate_ng.jl +++ b/src/methods/geom_relations/relateng/relate_ng.jl @@ -709,21 +709,15 @@ relate_predicate(p::PreparedRelate, predicate::TopologyPredicate, b) = # Whether to prebuild the A-side segment tree, mirroring the dispatch of # `process_edge_intersections!` + `_select_edge_set_accelerator`: an explicit # `NestedLoop` accelerator never uses a tree; `AutoAccelerator` uses one on -# `Planar` above the clipping size threshold only (B is unknown at prepare -# time, so the decision is made on A's segment count alone); any other -# explicit accelerator always takes the tree path. +# the manifolds with a segment-extent kernel (`Planar`/`Spherical`) above +# the clipping size threshold only (B is unknown at prepare time, so the +# decision is made on A's segment count alone); any other explicit +# accelerator always takes the tree path. _build_prepared_edge_index(m::Manifold, ::IntersectionAccelerator, segs_a) = _make_prepared_edge_index(m, segs_a) _build_prepared_edge_index(::Manifold, ::NestedLoop, segs_a) = nothing _build_prepared_edge_index(::Manifold, ::AutoAccelerator, segs_a) = nothing -function _build_prepared_edge_index(m::Planar, ::AutoAccelerator, segs_a) - _total_segment_count(segs_a) >= GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS || - return nothing - return _make_prepared_edge_index(m, segs_a) -end -#-- the Spherical tree path is valid too (3D arc extents); above the threshold -#-- prepared mode indexes A just like Planar -function _build_prepared_edge_index(m::Spherical, ::AutoAccelerator, segs_a) +function _build_prepared_edge_index(m::Union{Planar, Spherical}, ::AutoAccelerator, segs_a) _total_segment_count(segs_a) >= GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS || return nothing return _make_prepared_edge_index(m, segs_a) diff --git a/src/methods/geom_relations/relateng/topology_predicate.jl b/src/methods/geom_relations/relateng/topology_predicate.jl index 35e99859e..ca3246722 100644 --- a/src/methods/geom_relations/relateng/topology_predicate.jl +++ b/src/methods/geom_relations/relateng/topology_predicate.jl @@ -211,10 +211,9 @@ intersects_exterior_of(p::IMPredicate, is_a::Bool) = is_a ? (is_intersects_entry(p, LOC_INTERIOR, LOC_EXTERIOR) || is_intersects_entry(p, LOC_BOUNDARY, LOC_EXTERIOR)) is_intersects_entry(p::IMPredicate, locA, locB) = p.im[locA, locB] >= DIM_P -# NOTE: unused; kept for JTS IMPredicate API parity. As ported it can never -# return `false`: matrix entries are initialized to DIM_FALSE and only ever -# increase, so they never hold DIM_UNKNOWN (-3). Do not use as a real check. -is_known_entry(p::IMPredicate, locA, locB) = p.im[locA, locB] != DIM_UNKNOWN +# JTS's isKnownEntry is not ported: entries here are initialized to +# DIM_FALSE and only ever increase, so they never hold DIM_UNKNOWN and the +# check could never return false. is_dimension_entry(p::IMPredicate, locA, locB, dim) = p.im[locA, locB] == dim get_dimension(p::IMPredicate, locA, locB) = p.im[locA, locB] diff --git a/test/methods/relateng/edge_intersector.jl b/test/methods/relateng/edge_intersector.jl index f9308a58c..36adb3a33 100644 --- a/test/methods/relateng/edge_intersector.jl +++ b/test/methods/relateng/edge_intersector.jl @@ -346,7 +346,7 @@ ngon(cx, cy, r, n) = GI.Polygon([[ ] for (ga, gb) in fixtures tc_truth, pred_truth = run_all_pairs(ga, gb) - for acc in (GO.NestedLoop(), GO.DoubleSTRtree(), GO.AutoAccelerator()) + for acc in (GO.NestedLoop(), GO.DoubleNaturalTree(), GO.AutoAccelerator()) tc, pred = run_enum(ga, gb, acc) @test node_counts(tc) == node_counts(tc_truth) @test imstr(pred) == imstr(pred_truth) @@ -395,7 +395,7 @@ end # stop long before enumerating every extent-overlapping pair ga = ngon(0.0, 0.0, 1.0, 64) gb = ngon(0.1, 0.0, 1.0, 64) - for acc in (GO.NestedLoop(), GO.DoubleSTRtree()) + for acc in (GO.NestedLoop(), GO.DoubleNaturalTree()) # baseline: the full-matrix predicate is never known early, so its # count is the total number of pairs the enumeration would process full = run_counted(ga, gb, acc, GO.RelateMatrixPredicate()) @@ -436,7 +436,7 @@ end end counts_loop, im_loop = run_self(GO.NestedLoop()) - counts_tree, im_tree = run_self(GO.DoubleSTRtree()) + counts_tree, im_tree = run_self(GO.DoubleNaturalTree()) @test counts_loop == counts_tree @test im_loop == im_tree # sanity: the return segment properly crosses all 40 zigzag segments @@ -461,7 +461,7 @@ end GO.process_self_intersections!(tc, ss_list, acc) return pred end - for acc in (GO.NestedLoop(), GO.DoubleSTRtree()) + for acc in (GO.NestedLoop(), GO.DoubleNaturalTree()) # baseline: the full-matrix predicate is never known early, so its # count is the total number of pairs the enumeration would process full = run_self_counted(GO.RelateMatrixPredicate(), acc) diff --git a/test/methods/relateng/indexed_point_in_area.jl b/test/methods/relateng/indexed_point_in_area.jl index f3ab2e260..7dd10b366 100644 --- a/test/methods/relateng/indexed_point_in_area.jl +++ b/test/methods/relateng/indexed_point_in_area.jl @@ -95,28 +95,19 @@ mp = GI.MultiPolygon([ function check_prepared_agreement(geom, pts) m = Planar() loc_prep = GO.RelatePointLocator(m, geom; exact = True(), is_prepared = true) - # long-lived unprepared locator: repeat queries flip each polygonal - # element to the lazily built index (deviation from Java; see - # `locate_on_polygonal`) - loc_lazy = GO.RelatePointLocator(m, geom; exact = True(), is_prepared = false) - fresh() = GO.RelatePointLocator(m, geom; exact = True(), is_prepared = false) + loc_unprep = GO.RelatePointLocator(m, geom; exact = True(), is_prepared = false) n_mismatch = 0 for pt in pts - # fresh unprepared locators so each query is the element's first and - # takes the direct ring loop — a true indexed-vs-simple differential - GO.locate(loc_prep, pt) == GO.locate(fresh(), pt) || (n_mismatch += 1) - GO.locate_with_dim(loc_prep, pt) == GO.locate_with_dim(fresh(), pt) || (n_mismatch += 1) - # the lazy locator must agree whichever path it is on - GO.locate(loc_lazy, pt) == GO.locate(loc_prep, pt) || (n_mismatch += 1) - GO.locate_with_dim(loc_lazy, pt) == GO.locate_with_dim(loc_prep, pt) || (n_mismatch += 1) + # unprepared = direct ring loop, prepared = indexed locator — a true + # indexed-vs-simple differential on every query + GO.locate(loc_prep, pt) == GO.locate(loc_unprep, pt) || (n_mismatch += 1) + GO.locate_with_dim(loc_prep, pt) == GO.locate_with_dim(loc_unprep, pt) || (n_mismatch += 1) end @test n_mismatch == 0 - # cache sanity: prepared mode built (at most) one locator per polygonal element + # cache sanity: prepared mode built (at most) one locator per polygonal + # element; unprepared mode never builds one @test length(loc_prep.poly_locator) == length(loc_prep.polygons) - # the long-lived unprepared locator saw enough repeat queries to switch - # to the index on every polygonal element - @test all(>(GO._LAZY_INDEX_QUERY_THRESHOLD), loc_lazy.poly_query_count) - @test all(!isnothing, loc_lazy.poly_locator) + @test all(isnothing, loc_unprep.poly_locator) end function check_prepared_relate(geom, pts) diff --git a/test/methods/relateng/spherical_end_to_end.jl b/test/methods/relateng/spherical_end_to_end.jl index 6e59edc61..618ad8906 100644 --- a/test/methods/relateng/spherical_end_to_end.jl +++ b/test/methods/relateng/spherical_end_to_end.jl @@ -36,7 +36,7 @@ end @testset "spherical tree accelerator agrees with NestedLoop (arc bulge)" begin A = GI.Polygon([GI.LinearRing([(0., 0.), (170., 0.), (85., 40.), (0., 0.)])]) B = GI.Polygon([GI.LinearRing([(88., -2.), (92., -2.), (92., 2.), (88., 2.), (88., -2.)])]) - tree = RelateNG(; manifold = Spherical(), accelerator = GO.DoubleSTRtree()) + tree = RelateNG(; manifold = Spherical(), accelerator = GO.DoubleNaturalTree()) loop = RelateNG(; manifold = Spherical(), accelerator = GO.NestedLoop()) @test GO.relate(tree, A, B) == GO.relate(loop, A, B) end