From a31e5db55423bbaf3ea671549cfba9736379912f Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Fri, 14 Aug 2026 09:06:29 +0200 Subject: [PATCH] Check all clip vertices in the spherical Sutherland-Hodgman containment fallback MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When `_intersection_sutherland_hodgman` clips the subject down to nothing, it falls back to asking whether the clip polygon was inside the subject all along. That check only tested `clip_points[1]`, and `spherical_orient >= 0` counts the boundary as inside — so any subject merely touching that one vertex was credited with the clip polygon's whole area. Test every clip vertex instead. Under the documented convex-convex precondition this changes no answer: a sweep of 159600 convex pairs hits the empty-output branch 158732 times and never once finds `clip_points[1]` inside the subject, and genuine containment does not empty the output in the first place. It only costs O(m) point tests on a branch that is already the slow path. Fixes #467 Co-Authored-By: Claude Opus 5 --- src/methods/clipping/sutherland_hodgman.jl | 9 ++++-- test/methods/clipping/sutherland_hodgman.jl | 33 +++++++++++++++++++++ 2 files changed, 40 insertions(+), 2 deletions(-) diff --git a/src/methods/clipping/sutherland_hodgman.jl b/src/methods/clipping/sutherland_hodgman.jl index a4649074c0..e81869808c 100644 --- a/src/methods/clipping/sutherland_hodgman.jl +++ b/src/methods/clipping/sutherland_hodgman.jl @@ -314,9 +314,14 @@ function _intersection_sutherland_hodgman( output = _sh_clip_to_edge_spherical(output, edge_start, edge_end, T) end - # Handle empty result - check if clip polygon is fully inside the original subject + # Handle empty result - check if clip polygon is fully inside the original subject. + # "Fully inside" needs EVERY clip vertex inside the subject, not just the first: + # one vertex lying on the subject's boundary is a graze, and `spherical_orient >= 0` + # counts the boundary as inside, so testing `clip_points[1]` alone credits every + # subject that merely touches that one vertex with the clip polygon's whole area. if isempty(output) - if !isempty(clip_points) && _point_in_convex_spherical_polygon(clip_points[1], original_subject) + if !isempty(clip_points) && + all(p -> _point_in_convex_spherical_polygon(p, original_subject), clip_points) # Subject contains clip - return clip polygon result = copy(clip_points) push!(result, result[1]) diff --git a/test/methods/clipping/sutherland_hodgman.jl b/test/methods/clipping/sutherland_hodgman.jl index 0f5f13cfec..cb2a492f31 100644 --- a/test/methods/clipping/sutherland_hodgman.jl +++ b/test/methods/clipping/sutherland_hodgman.jl @@ -379,6 +379,39 @@ import GeoInterface as GI @test a_oh ≈ a_ho rtol = 1e-6 @test 0 < a_oh < 0.1 * spherical_area(healpix) # sliver, not the whole cell end + + @testset "Grazing tiles are not handed the whole clip" begin + # The empty-output fallback returns the clip polygon when the subject + # contains it. "Contains" needs EVERY clip vertex inside the subject: + # `spherical_orient >= 0` counts the boundary as inside, so testing only + # `clip_points[1]` credits every subject that merely touches that one + # vertex with the clip's whole area. + # + # Reaching the fallback at all takes a non-convex clip -- with a convex + # clip a grazing subject survives every half-space test and leaves a + # 1-2 point output, which the degenerate branch already zeroes. So this + # is a graceful-degradation test for out-of-contract input (a cell that + # goes non-convex numerically, see #466): the answer may be wrong, but it + # must be wrong towards zero, never towards a whole fabricated cell. + clip = spherical_polygon([(0.0, 0.0), (20.0, 0.0), (20.0, 20.0), (10.0, 6.0), (0.0, 20.0)]) + + # Three convex, pairwise-disjoint tiles of the quadrants meeting at the + # clip ring's FIRST vertex (0, 0). None overlaps the clip's interior. + tiles = [ + spherical_polygon([(-5.0, -5.0), (0.0, -5.0), (0.0, 0.0), (-5.0, 0.0)]), # SW + spherical_polygon([(0.0, -5.0), (5.0, -5.0), (5.0, 0.0), (0.0, 0.0)]), # SE, shares the clip's south edge + spherical_polygon([(-5.0, 0.0), (0.0, 0.0), (0.0, 5.0), (-5.0, 5.0)]), # NW, shares the clip's west edge + ] + + alg = GO.ConvexConvexSutherlandHodgman(GO.Spherical()) + clip_a = spherical_area(clip) + for tile in tiles + @test spherical_area(GO.intersection(alg, tile, clip)) < 1e-9 * clip_a + end + # Summed over the disjoint tiles this is what breaks conservative + # regridding: the total was 3x the clip's own area. + @test sum(spherical_area(GO.intersection(alg, tile, clip)) for tile in tiles) < 1e-9 * clip_a + end end end