Fix cut

Author: asinghvi17

src/methods/clipping/clipping_processor.jl

@@ -892,7 +892,7 @@ or they are separate polygons with no intersection (other than an edge or point)
 Return two booleans that represent if a is inside b (potentially with shared edges / points)
 and visa versa if b is inside of a.
 =#
-function _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact)
+function _find_non_cross_orientation(m::M, a_list, b_list, a_poly, b_poly; exact) where {M <: Manifold}
     non_intr_a_idx = findfirst(x -> !x.inter, a_list)
     non_intr_b_idx = findfirst(x -> !x.inter, b_list)
     #= Determine if non-intersection point is in or outside of polygon - if there isn't A
@@ -906,6 +906,9 @@ function _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact)
     return a_in_b, b_in_a
 end
 
+_find_non_cross_orientation(alg::FosterHormannClipping{M}, a_list, b_list, a_poly, b_poly; exact) where {M <: Manifold} =
+    _find_non_cross_orientation(alg.manifold, a_list, b_list, a_poly, b_poly; exact)
+
 #=
     _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_poly_idx; exact)
 

src/methods/clipping/cut.jl

@@ -58,29 +58,32 @@ GI.coordinates.(cut_polys)
  [[[5.0, 0.0], [10.0, 0.0], [10.0, 10.0], [5.0, 10.0], [5.0, 0.0]]]
 ```
 """
-cut(geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} =
-    _cut(T, GI.trait(geom), geom, GI.trait(line), line; exact = _True())
+cut(geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = cut(FosterHormannClipping(), geom, line, T)
+cut(m::Manifold, geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = cut(FosterHormannClipping(m), geom, line, T)
+
+cut(alg::FosterHormannClipping{M, A}, geom, line, ::Type{T} = Float64) where {T <: AbstractFloat, M, A} =
+    _cut(alg, T, GI.trait(geom), geom, GI.trait(line), line; exact = _True())
 
 #= Cut a given polygon by given line. Add polygon holes back into resulting pieces if there
 are any holes. =#
-function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) where T
+function _cut(alg::FosterHormannClipping{M, A}, ::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) where {T, M, A}
     ext_poly = GI.getexterior(poly)
-    poly_list, intr_list = _build_a_list(T, ext_poly, line; exact)
+    poly_list, intr_list = _build_a_list(alg, T, ext_poly, line; exact)
     n_intr_pts = length(intr_list)
     # If an impossible number of intersection points, return original polygon
     if n_intr_pts < 2 || isodd(n_intr_pts)
         return [tuples(poly)]
     end
     # Cut polygon by line
-    cut_coords = _cut(T, ext_poly, line, poly_list, intr_list, n_intr_pts; exact)
+    cut_coords = _cut(alg, T, ext_poly, line, poly_list, intr_list, n_intr_pts; exact)
     # Close coords and create polygons
     for c in cut_coords
         push!(c, c[1])
     end
     cut_polys = [GI.Polygon([c]) for c in cut_coords]
     # Add original polygon holes back in
     remove_idx = falses(length(cut_polys))
-    _add_holes_to_polys!(T, cut_polys, GI.gethole(poly), remove_idx; exact)
+    _add_holes_to_polys!(alg, T, cut_polys, GI.gethole(poly), remove_idx; exact)
     return cut_polys
 end
 
@@ -97,10 +100,10 @@ end
 of cut geometry in Vector{Vector{Tuple}} format. 
 
 Note: degenerate cases where intersection points are vertices do not work right now. =#
-function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where T
+function _cut(alg::FosterHormannClipping{M, A}, ::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where {T, M, A}
     # Sort and categorize the intersection points
     sort!(intr_list, by = x -> geom_list[x].fracs[2])
-    _flag_ent_exit!(GI.LineTrait(), line, geom_list; exact)
+    _flag_ent_exit!(alg, GI.LineTrait(), line, geom_list; exact)
     # Add first point to output list
     return_coords = [[geom_list[1].point]]
     cross_backs = [(T(Inf),T(Inf))]

test/methods/clipping/polygon_clipping.jl

@@ -166,7 +166,7 @@ test_pairs = [
 const ϵ = 1e-10
 # Compare clipping results from GeometryOps and LibGEOS
 function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
-    GO_result_list = GO_f(p1, p2; target = GI.PolygonTrait())
+    
     LG_result_geom = LG_f(p1, p2)
     if LG_result_geom isa LG.GeometryCollection
         poly_list = LG.Polygon[]
@@ -175,38 +175,43 @@ function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
         end
         LG_result_geom = LG.MultiPolygon(poly_list)
     end
-    # Check if nothing is returned
-    if isempty(GO_result_list) && (LG.isEmpty(LG_result_geom) || LG.area(LG_result_geom) == 0)
-        return true
-    end
-    # Check for unnecessary points
-    if sum(GI.npoint, GO_result_list; init = 0.0) > GI.npoint(LG_result_geom)
-        return false
-    end
-    # Make sure last point is repeated
-    for poly in GO_result_list
-        for ring in GI.getring(poly)
-            GI.getpoint(ring, 1) != GI.getpoint(ring, GI.npoint(ring)) && return false
+
+    for _accelerator in (GO.AutoAccelerator(), GO.NestedLoop(), GO.SingleSTRtree(), GO.DoubleSTRtree())
+    @testset let accelerator = _accelerator # this is a ContextTestSet that is otherwise invisible but adds context to the testset
+        GO_result_list = GO_f(GO.FosterHormannClipping(accelerator), p1, p2; target = GI.PolygonTrait())
+        # Check if nothing is returned
+        if isempty(GO_result_list) && (LG.isEmpty(LG_result_geom) || LG.area(LG_result_geom) == 0)
+            @test true
+            continue
+        end
+        # Check for unnecessary points
+        @test !(sum(GI.npoint, GO_result_list; init = 0.0) > GI.npoint(LG_result_geom))
+        # Make sure last point is repeated
+        for poly in GO_result_list
+            for ring in GI.getring(poly)
+                @test !(GI.getpoint(ring, 1) != GI.getpoint(ring, GI.npoint(ring)))
+            end
         end
-    end
 
-    # Check if polygons cover the same area
-    local GO_result_geom
-    if length(GO_result_list) == 1
-        GO_result_geom = GO_result_list[1]
-    else
-        GO_result_geom = GI.MultiPolygon(GO_result_list)
-    end
-    diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
-    diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
-    return diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+        # Check if polygons cover the same area
+        local GO_result_geom
+        if length(GO_result_list) == 1
+            GO_result_geom = GO_result_list[1]
+        else
+            GO_result_geom = GI.MultiPolygon(GO_result_list)
+        end
+        diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
+        diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
+        @test diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+    end # testset
+    end # loop
 end
 
 # Test clipping functions and print error message if tests fail
 function test_clipping(GO_f, LG_f, f_name)
     for (p1, p2, sg1, sg2, sdesc) in test_pairs
         @testset_implementations "$sg1 $f_name $sg2 - $sdesc" begin
-            @test compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) 
+            compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) # this executes tests internally
         end
     end
     return