Add allocation cache for ConvexConvexSutherlandHodgman clipping

src/methods/clipping/sutherland_hodgman.jl

@@ -1,5 +1,5 @@
 # # Sutherland-Hodgman Convex-Convex Clipping
-export ConvexConvexSutherlandHodgman
+export ConvexConvexSutherlandHodgman, SutherlandHodgmanCache
 
 """
     ConvexConvexSutherlandHodgman{M <: Manifold} <: GeometryOpsCore.Algorithm{M}
@@ -11,6 +11,9 @@ Both input polygons MUST be convex. If either polygon is non-convex, results are
 This is simpler and faster than Foster-Hormann for small convex polygons, with O(n*m)
 complexity where n and m are vertex counts.
 
+When intersecting many polygon pairs in a hot loop, pass a [`SutherlandHodgmanCache`](@ref)
+via the `cache` keyword argument to `intersection` to avoid intermediate allocations.
+
 ## Spherical manifold
 
 For `Spherical()` manifold, input polygons must have **counter-clockwise winding** when
@@ -37,18 +40,78 @@ end
 # Default constructor uses Planar
 ConvexConvexSutherlandHodgman() = ConvexConvexSutherlandHodgman(Planar())
 
+"""
+    SutherlandHodgmanCache{P}()
+    SutherlandHodgmanCache(manifold::Manifold, [T = Float64])
+    SutherlandHodgmanCache(alg::ConvexConvexSutherlandHodgman, [T = Float64])
+
+Preallocated buffers for [`ConvexConvexSutherlandHodgman`](@ref) clipping.
+
+Pass this as the `cache` keyword argument to `intersection` to avoid allocating
+intermediate vectors on every call - useful when intersecting many polygon pairs
+in a hot loop. The returned polygon never aliases the cache, so results remain
+valid after the cache is reused.
+
+The point type `P` must match the algorithm's manifold and float type:
+`Tuple{T,T}` for `Planar()`, `UnitSpherical.UnitSphericalPoint{T}` for
+`Spherical()`. The manifold/algorithm constructors take care of this.
+
+!!! warning "Thread safety"
+    A cache must not be shared across concurrent tasks or threads. Create one
+    cache per task. The default (`cache = nothing`) allocates fresh buffers on
+    each call and is always safe.
+
+# Example
+
+```julia
+import GeometryOps as GO
+
+alg = GO.ConvexConvexSutherlandHodgman()
+cache = GO.SutherlandHodgmanCache(alg)
+for (a, b) in polygon_pairs
+    result = GO.intersection(alg, a, b; cache)
+end
+```
+"""
+struct SutherlandHodgmanCache{P}
+    input::Vector{P}    # ping buffer
+    output::Vector{P}   # pong buffer
+    clip::Vector{P}     # spherical only: clip polygon vertices
+    subject::Vector{P}  # spherical only: copy of original subject for containment check
+end
+SutherlandHodgmanCache{P}() where P = SutherlandHodgmanCache{P}(P[], P[], P[], P[])
+
+SutherlandHodgmanCache(::Planar, ::Type{T} = Float64) where {T} =
+    SutherlandHodgmanCache{Tuple{T,T}}()
+SutherlandHodgmanCache(::Spherical, ::Type{T} = Float64) where {T} =
+    SutherlandHodgmanCache{UnitSpherical.UnitSphericalPoint{T}}()
+SutherlandHodgmanCache(alg::ConvexConvexSutherlandHodgman, ::Type{T} = Float64) where {T} =
+    SutherlandHodgmanCache(alg.manifold, T)
+
+# Validate that a user-supplied cache has the point type the manifold/float combination needs
+function _sh_check_cache(cache::SutherlandHodgmanCache{P}, ::Type{PT}) where {P, PT}
+    P === PT || throw(ArgumentError(
+        "SutherlandHodgmanCache point type mismatch: this intersection requires " *
+        "SutherlandHodgmanCache{$PT}, got SutherlandHodgmanCache{$P}. " *
+        "Construct the cache with `SutherlandHodgmanCache(alg, T)` to match the algorithm."
+    ))
+    return cache
+end
+
 # Main entry point - algorithm dispatch
 function intersection(
     alg::ConvexConvexSutherlandHodgman,
     geom_a,
     geom_b,
     ::Type{T}=Float64;
+    cache::Union{Nothing, SutherlandHodgmanCache}=nothing,
     kwargs...
 ) where {T<:AbstractFloat}
     return _intersection_sutherland_hodgman(
         alg, T,
         GI.trait(geom_a), geom_a,
-        GI.trait(geom_b), geom_b
+        GI.trait(geom_b), geom_b;
+        cache
     )
 end
 
@@ -57,51 +120,60 @@ function _intersection_sutherland_hodgman(
     alg::ConvexConvexSutherlandHodgman{Planar},
     ::Type{T},
     ::GI.PolygonTrait, poly_a,
-    ::GI.PolygonTrait, poly_b
+    ::GI.PolygonTrait, poly_b;
+    cache::Union{Nothing, SutherlandHodgmanCache}=nothing
 ) where {T}
+    cache = isnothing(cache) ? SutherlandHodgmanCache(Planar(), T) : _sh_check_cache(cache, Tuple{T,T})
+
     # Get exterior rings (convex polygons have no holes)
     ring_a = GI.getexterior(poly_a)
     ring_b = GI.getexterior(poly_b)
 
-    # Start with vertices of poly_a as the output list (excluding closing point)
-    output = Tuple{T,T}[]
+    # Start with vertices of poly_a as the subject list (excluding closing point),
+    # ping-ponging between the two cache buffers as we clip
+    buf_in, buf_out = cache.input, cache.output
+    empty!(buf_in)
     for point in GI.getpoint(ring_a)
         pt = _tuple_point(point, T)
         # Skip the closing point (same as first)
-        if !isempty(output) && pt == output[1]
+        if !isempty(buf_in) && pt == buf_in[1]
             continue
         end
-        push!(output, pt)
+        push!(buf_in, pt)
     end
 
     # Clip against each edge of poly_b
     for (edge_start, edge_end) in eachedge(ring_b, T)
-        isempty(output) && break
-        output = _sh_clip_to_edge(output, edge_start, edge_end, T)
+        isempty(buf_in) && break
+        _sh_clip_to_edge!(buf_out, buf_in, edge_start, edge_end, T)
+        buf_in, buf_out = buf_out, buf_in
     end
 
     # Handle empty result (no intersection) - return degenerate polygon with zero area
-    if isempty(output)
+    if isempty(buf_in)
         zero_pt = (zero(T), zero(T))
         return GI.Polygon([[zero_pt, zero_pt, zero_pt]])
     end
 
-    # Close the ring
-    push!(output, output[1])
+    # Copy into a fresh closed ring so the result doesn't alias the cache
+    result = Vector{Tuple{T,T}}(undef, length(buf_in) + 1)
+    copyto!(result, buf_in)
+    result[end] = buf_in[1]
 
     # Return polygon
-    return GI.Polygon([output])
+    return GI.Polygon([result])
 end
 
-# Clip polygon against a single edge using Sutherland-Hodgman rules
-function _sh_clip_to_edge(polygon_points::Vector{Tuple{T,T}}, edge_start, edge_end, ::Type{T}) where T
-    output = Tuple{T,T}[]
-    n = length(polygon_points)
+# Clip polygon (read from `input`) against a single edge using Sutherland-Hodgman
+# rules, writing the surviving points into `output`
+function _sh_clip_to_edge!(output::Vector{Tuple{T,T}}, input::Vector{Tuple{T,T}}, edge_start, edge_end, ::Type{T}) where T
+    empty!(output)
+    n = length(input)
     n == 0 && return output
 
     for i in 1:n
-        current = polygon_points[i]
-        next_pt = polygon_points[mod1(i + 1, n)]
+        current = input[i]
+        next_pt = input[mod1(i + 1, n)]
 
         # Determine if points are inside (left of or on the edge)
         # orient > 0 means left (inside for CCW polygon), == 0 means on edge, < 0 means right (outside)
@@ -221,24 +293,26 @@ function _sh_spherical_intersection(
     return UnitSpherical.UnitSphericalPoint{T}(p1)
 end
 
-# Clip polygon against a single edge using Sutherland-Hodgman rules (spherical version)
-function _sh_clip_to_edge_spherical(
-    polygon_points::Vector{UnitSpherical.UnitSphericalPoint{T}},
+# Clip polygon (read from `input`) against a single edge using Sutherland-Hodgman
+# rules, writing the surviving points into `output` (spherical version)
+function _sh_clip_to_edge_spherical!(
+    output::Vector{UnitSpherical.UnitSphericalPoint{T}},
+    input::Vector{UnitSpherical.UnitSphericalPoint{T}},
     edge_start::UnitSpherical.UnitSphericalPoint,
     edge_end::UnitSpherical.UnitSphericalPoint,
     ::Type{T}
 ) where T
-    output = UnitSpherical.UnitSphericalPoint{T}[]
-    n = length(polygon_points)
+    empty!(output)
+    n = length(input)
     n == 0 && return output
 
     # Track actual orient values to handle edge cases (orient=0 means exactly on the edge)
-    current_orient = UnitSpherical.spherical_orient(edge_start, edge_end, polygon_points[1])
+    current_orient = UnitSpherical.spherical_orient(edge_start, edge_end, input[1])
 
     for i in 1:n
-        current = polygon_points[i]
+        current = input[i]
         next_idx = mod1(i + 1, n)
-        next_pt = polygon_points[next_idx]
+        next_pt = input[next_idx]
 
         next_orient = UnitSpherical.spherical_orient(edge_start, edge_end, next_pt)
         current_inside = current_orient >= 0
@@ -276,8 +350,11 @@ function _intersection_sutherland_hodgman(
     alg::ConvexConvexSutherlandHodgman{Spherical{F}},
     ::Type{T},
     ::GI.PolygonTrait, poly_a,
-    ::GI.PolygonTrait, poly_b
+    ::GI.PolygonTrait, poly_b;
+    cache::Union{Nothing, SutherlandHodgmanCache}=nothing
 ) where {F, T}
+    cache = isnothing(cache) ? SutherlandHodgmanCache(alg.manifold, T) : _sh_check_cache(cache, UnitSpherical.UnitSphericalPoint{T})
+
     ring_a = GI.getexterior(poly_a)
     ring_b = GI.getexterior(poly_b)
 
@@ -291,41 +368,47 @@ function _intersection_sutherland_hodgman(
     end
 
     # Collect clip polygon points (excluding closing point)
-    clip_points = UnitSpherical.UnitSphericalPoint{T}[]
+    clip_points = cache.clip
+    empty!(clip_points)
     for point in GI.getpoint(ring_b)
         if !isempty(clip_points) && point ≈ clip_points[1]
             continue
         end
         push!(clip_points, UnitSpherical.UnitSphericalPoint{T}(point))
     end
 
-    # Build initial output list from poly_a (excluding closing point)
-    output = UnitSpherical.UnitSphericalPoint{T}[]
+    # Build initial subject list from poly_a (excluding closing point),
+    # ping-ponging between the two cache buffers as we clip
+    buf_in, buf_out = cache.input, cache.output
+    empty!(buf_in)
     for point in GI.getpoint(ring_a)
-        if !isempty(output) && point == output[1]
+        if !isempty(buf_in) && point == buf_in[1]
             continue
         end
-        push!(output, UnitSpherical.UnitSphericalPoint{T}(point))
+        push!(buf_in, UnitSpherical.UnitSphericalPoint{T}(point))
     end
 
     # Save original subject for containment check
-    original_subject = copy(output)
+    original_subject = cache.subject
+    empty!(original_subject)
+    append!(original_subject, buf_in)
 
     # Clip against each edge of poly_b
     n_clip = length(clip_points)
     for i in 1:n_clip
-        isempty(output) && break
+        isempty(buf_in) && break
         edge_start = clip_points[i]
         edge_end = clip_points[mod1(i + 1, n_clip)]
         # Skip degenerate edges (duplicate vertices)
         edge_start == edge_end && continue
-        output = _sh_clip_to_edge_spherical(output, edge_start, edge_end, T)
+        _sh_clip_to_edge_spherical!(buf_out, buf_in, edge_start, edge_end, T)
+        buf_in, buf_out = buf_out, buf_in
     end
 
     # Handle empty result - check if clip polygon is fully inside the original subject
-    if isempty(output)
+    if isempty(buf_in)
         if !isempty(clip_points) && _point_in_convex_spherical_polygon(clip_points[1], original_subject)
-            # Subject contains clip - return clip polygon
+            # Subject contains clip - return clip polygon (copied, so it doesn't alias the cache)
             result = copy(clip_points)
             push!(result, result[1])
             return GI.Polygon([result])
@@ -337,22 +420,25 @@ function _intersection_sutherland_hodgman(
 
     # Handle degenerate result (1-2 points can't form a valid ring)
     # This happens for adjacent polygons that share only an edge or corner
-    if length(output) < 3
+    if length(buf_in) < 3
         north_pole = UnitSpherical.UnitSphericalPoint{T}(0, 0, 1)
         return GI.Polygon([[north_pole, north_pole, north_pole]])
     end
 
-    # Close the ring and return
-    push!(output, output[1])
-    return GI.Polygon([output])
+    # Copy into a fresh closed ring so the result doesn't alias the cache
+    result = Vector{UnitSpherical.UnitSphericalPoint{T}}(undef, length(buf_in) + 1)
+    copyto!(result, buf_in)
+    result[end] = buf_in[1]
+    return GI.Polygon([result])
 end
 
 # Fallback for unsupported geometry combinations
 function _intersection_sutherland_hodgman(
     alg::ConvexConvexSutherlandHodgman,
     ::Type{T},
     trait_a, geom_a,
-    trait_b, geom_b
+    trait_b, geom_b;
+    kwargs...
 ) where {T}
     throw(ArgumentError(
         "ConvexConvexSutherlandHodgman only supports Polygon-Polygon intersection, " *

test/methods/clipping/sutherland_hodgman.jl

@@ -181,6 +181,78 @@ import GeoInterface as GI
         @test GO.area(result3) ≈ 0.0 atol=1e-10
     end
 
+    @testset "Cache" begin
+        alg = GO.ConvexConvexSutherlandHodgman()
+        square1 = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
+        square2 = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]])
+        tri1 = GI.Polygon([[(0.0, 0.0), (4.0, 0.0), (2.0, 4.0), (0.0, 0.0)]])
+        tri2 = GI.Polygon([[(0.0, 2.0), (2.0, -2.0), (4.0, 2.0), (0.0, 2.0)]])
+
+        @testset "Same results with and without cache" begin
+            cache = GO.SutherlandHodgmanCache(alg)
+            for (a, b) in [(square1, square2), (tri1, tri2), (square1, tri1)]
+                uncached = GO.intersection(alg, a, b)
+                cached = GO.intersection(alg, a, b; cache)
+                @test collect(GI.getpoint(GI.getexterior(cached))) ==
+                      collect(GI.getpoint(GI.getexterior(uncached)))
+            end
+        end
+
+        @testset "Result does not alias the cache" begin
+            cache = GO.SutherlandHodgmanCache(alg)
+            first_result = GO.intersection(alg, square1, square2; cache)
+            first_points = collect(GI.getpoint(GI.getexterior(first_result)))
+            # Reuse the cache - this must not corrupt the previously returned polygon
+            GO.intersection(alg, tri1, tri2; cache)
+            @test collect(GI.getpoint(GI.getexterior(first_result))) == first_points
+        end
+
+        @testset "Wrong cache point type throws" begin
+            f32_cache = GO.SutherlandHodgmanCache(GO.Planar(), Float32)
+            @test_throws ArgumentError GO.intersection(alg, square1, square2; cache=f32_cache)
+
+            spherical_cache = GO.SutherlandHodgmanCache(GO.Spherical())
+            @test_throws ArgumentError GO.intersection(alg, square1, square2; cache=spherical_cache)
+        end
+
+        @testset "No intermediate allocations after warm-up" begin
+            cache = GO.SutherlandHodgmanCache(alg)
+            GO.intersection(alg, square1, square2; cache)  # warm up buffers and compilation
+            allocated = @allocated GO.intersection(alg, square1, square2; cache)
+            # Only the returned polygon should allocate; loose bound to avoid version flakiness
+            @test allocated <= 512
+        end
+
+        @testset "Spherical cache" begin
+            using GeometryOps.UnitSpherical: UnitSphereFromGeographic
+
+            function spherical_polygon(coords)
+                transform = UnitSphereFromGeographic()
+                points = [transform((lon, lat)) for (lon, lat) in coords]
+                push!(points, points[1])
+                return GI.Polygon([points])
+            end
+
+            salg = GO.ConvexConvexSutherlandHodgman(GO.Spherical())
+            scache = GO.SutherlandHodgmanCache(salg)
+            sq1 = spherical_polygon([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0)])
+            sq2 = spherical_polygon([(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0)])
+            large = spherical_polygon([(-5.0, -5.0), (5.0, -5.0), (5.0, 5.0), (-5.0, 5.0)])
+
+            for (a, b) in [(sq1, sq2), (large, sq1), (sq1, large)]
+                uncached = GO.intersection(salg, a, b)
+                cached = GO.intersection(salg, a, b; cache=scache)
+                @test GO.area(GO.Spherical(), cached) ≈ GO.area(GO.Spherical(), uncached)
+            end
+
+            # Aliasing check on the spherical path too
+            first_result = GO.intersection(salg, sq1, sq2; cache=scache)
+            first_points = collect(GI.getpoint(GI.getexterior(first_result)))
+            GO.intersection(salg, large, sq1; cache=scache)
+            @test collect(GI.getpoint(GI.getexterior(first_result))) == first_points
+        end
+    end
+
     @testset "Spherical helpers" begin
         using GeometryOps.UnitSpherical: UnitSphericalPoint, UnitSphereFromGeographic