Sg/cut poly (#56)

* Cup polygon with line code * Change target type * Add cut tests that throw warnings * replace str name * Remove warning * Add more comments Co-authored-by: Anshul Singhvi <[email protected]> * Fix comments * Add line type dispatch --------- Co-authored-by: Anshul Singhvi <[email protected]>
JuliaGeo · Feb 12, 2024 · 854b036 · 854b036
1 parent 84e9d53
commit 854b036
Show file tree

Hide file tree

Showing 7 changed files with 207 additions and 6 deletions.
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
@@ -29,6 +29,7 @@ include("methods/centroid.jl")
 include("methods/distance.jl")
 include("methods/equals.jl")
 include("methods/clipping/clipping_processor.jl")
+include("methods/clipping/cut.jl")
 include("methods/clipping/intersection.jl")
 include("methods/clipping/difference.jl")
 include("methods/clipping/union.jl")

diff --git a/src/methods/clipping/cut.jl b/src/methods/clipping/cut.jl
@@ -0,0 +1,132 @@
+# # Polygon cutting
+
+export cut
+
+#=
+## What is cut?
+
+The cut function cuts a polygon through a line segment. This is inspired by functions such
+as Matlab's [`cutpolygon`](https://www.mathworks.com/matlabcentral/fileexchange/24449-cutpolygon)
+function.
+
+To provide an example, consider the following polygon and line:
+```julia
+import GeoInterface as GI, GeometryOps as GO
+using CairoMakie
+using Makie
+
+poly = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]])
+line = GI.Line([(5.0, -5.0), (5.0, 15.0)])
+cut_polys = GO.cut(poly, line)
+
+f, a, p1 = Makie.poly(collect(GI.getpoint(cut_polys[1])); color = :blue)
+Makie.poly!(collect(GI.getpoint(cut_polys[2])); color = :orange)
+Makie.lines!(GI.getpoint(line); color = :black)
+f
+```
+
+## Implementation
+
+This function depends on polygon clipping helper function and is inspired by the
+Greiner-Hormann clipping algorithm used elsewhere in this library. The inspiration came from
+[this](https://stackoverflow.com/questions/3623703/how-can-i-split-a-polygon-by-a-line)
+Stack Overflow discussion. 
+=#
+
+"""
+    cut(geom, line, [T::Type])
+
+Return given geom cut by given line as a list of geometries of the same type as the input
+geom. Return the original geometry as only list element if none are found. Line must cut
+fully through given geometry or the original geometry will be returned.
+
+## Example 
+
+```jldoctest
+import GeoInterface as GI, GeometryOps as GO
+
+poly = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]])
+line = GI.Line([(5.0, -5.0), (5.0, 15.0)])
+cut_polys = GO.cut(poly, line)
+GI.coordinates.(cut_polys)
+
+# output
+2-element Vector{Vector{Vector{Vector{Float64}}}}:
+ [[[0.0, 0.0], [5.0, 0.0], [5.0, 10.0], [0.0, 10.0], [0.0, 0.0]]]
+ [[[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)
+
+#= 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) where T
+    ext_poly = GI.getexterior(poly)
+    poly_list, intr_list = _build_a_list(T, ext_poly, line)
+    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, poly_list, intr_list, n_intr_pts)
+    # 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
+    _add_holes_to_polys!(T, cut_polys, GI.gethole(poly))
+    return cut_polys
+end
+
+# Many types aren't implemented
+function _cut(_, trait::GI.AbstractTrait, geom, line)
+    @assert(
+        false,
+        "Cutting of $trait isn't implemented yet.",
+    )
+    return nothing
+end
+
+#= Cutting algorithm inspired by Greiner and Hormann clipping algorithm. Returns coordinates
+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, geom_list, intr_list, n_intr_pts) where T
+    # Sort and catagorize the intersection points
+    sort!(intr_list, by = x -> geom_list[x].fracs[2])
+    _flag_ent_exit!(geom, geom_list)
+    # Add first point to output list
+    return_coords = [[geom_list[1].point]]
+    cross_backs = [(T(Inf),T(Inf))]
+    poly_idx = 1
+    n_polys = 1
+    # Walk around original polygon to find split polygons
+    for (pt_idx, curr) in enumerate(geom_list)
+        if pt_idx > 1
+            push!(return_coords[poly_idx], curr.point)
+        end
+        if curr.inter
+            # Find cross back point for current polygon
+            intr_idx = findfirst(x -> equals(curr.point, geom_list[x].point), intr_list)
+            cross_idx = intr_idx + (curr.ent_exit ? 1 : -1)
+            cross_idx = cross_idx < 1 ? n_intr_pts : cross_idx
+            cross_idx = cross_idx > n_intr_pts ? 1 : cross_idx
+            cross_backs[poly_idx] = geom_list[intr_list[cross_idx]].point
+            # Check if current point is a cross back point
+            next_poly_idx = findfirst(x -> equals(x, curr.point), cross_backs)
+            if isnothing(next_poly_idx)
+                push!(return_coords, [curr.point])
+                push!(cross_backs, curr.point)
+                n_polys += 1
+                poly_idx = n_polys
+            else
+                push!(return_coords[next_poly_idx], curr.point)
+                poly_idx = next_poly_idx
+            end
+        end
+    end
+    return return_coords
+end
diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl
@@ -28,7 +28,7 @@ GI.coordinates.(diff_poly)
 """
 function difference(
     geom_a, geom_b, ::Type{T} = Float64; target::Type{Target} = Nothing,
-) where {T <: AbstractFloat, Target <: GI.AbstractTrait}
+) where {T <: AbstractFloat, Target <: Union{Nothing, GI.AbstractTrait}}
     return _difference(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
 end
 
@@ -77,7 +77,8 @@ function _difference(
 ) where {Target, T}
     @assert(
         false,
-        "Intersection between $trait_a and $trait_b with target $Target isn't implemented yet.",
+        "Difference between $trait_a and $trait_b with target $Target isn't implemented yet.",
     )
     return nothing
-end
+end
+
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
@@ -28,7 +28,7 @@ GI.coordinates.(inter_points)
 """
 function intersection(
     geom_a, geom_b, ::Type{T} = Float64; target::Type{Target} = Nothing,
-) where {T <: AbstractFloat, Target <: GI.AbstractTrait}
+) where {T <: AbstractFloat, Target <: Union{Nothing, GI.AbstractTrait}}
     return _intersection(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
 end
 

diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl
@@ -28,7 +28,7 @@ GI.coordinates.(union_poly)
 """
 function union(
     geom_a, geom_b, ::Type{T} = Float64; target::Type{Target} = Nothing,
-) where {T <: AbstractFloat, Target <: GI.AbstractTrait}
+) where {T <: AbstractFloat, Target <: Union{Nothing, GI.AbstractTrait}}
     _union(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
 end
 
@@ -87,6 +87,6 @@ function _union(
     trait_a::GI.AbstractTrait, geom_a,
     trait_b::GI.AbstractTrait, geom_b,
 ) where {Target, T}
-    throw(ArgumentError("Intersection between $trait_a and $trait_b with target $Target isn't implemented yet."))
+    throw(ArgumentError("Union between $trait_a and $trait_b with target $Target isn't implemented yet."))
     return nothing
 end
diff --git a/test/methods/clipping/cut.jl b/test/methods/clipping/cut.jl
@@ -0,0 +1,66 @@
+l1 = GI.Line([(5.0, -5.0), (5.0, 15.0)])
+l2 = GI.Line([(-1.0, 8.0), (2.0, 11.0)])
+l3 = GI.Line([(-1.0, 6.0), (11.0, 6.0)])
+l4 = GI.Line([(-10.0, -10.0), (-1.0, -1.0)])
+l5 = GI.Line([(-10.0, 5.0), (5.0, 5.0)])
+r1 = GI.LinearRing([(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)])
+r2 = GI.LinearRing([(2.0, 2.0), (8.0, 2.0), (8.0, 4.0), (2.0, 4.0), (2.0, 2.0)])
+p1 = GI.Polygon([r1])
+p2 = GI.Polygon([r1, r2])
+p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (10.0, 5.0), (10.0, 0.0), (0.0, 0.0)]])
+
+@testset "Cut Polygon" begin
+    # Cut convex polygon
+    cut_polys = GO.cut(p1, l1)
+    @test all(GO.equals.(
+        cut_polys,
+        [
+            GI.Polygon([[(0.0, 0.0), (5.0, 0.0), (5.0, 10.0), (0.0, 10.0), (0.0, 0.0)]]),
+            GI.Polygon([[(5.0, 0.0), (10.0, 0.0), (10.0, 10.0), (5.0, 10.0), (5.0, 0.0)]]),
+        ],
+    ))
+
+    cut_polys = GO.cut(p1, l2)
+    @test all(GO.equals.(
+        cut_polys,
+        [
+            GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (1.0, 10.0), (0.0, 9.0), (0.0, 0.0)]]),
+            GI.Polygon([[(0.0, 9.0), (0.0, 10.0), (1.0, 10.0), (0.0, 9.0)]]),
+        ],
+    ))
+
+    # Cut convex polygon with hole
+    cut_polys = GO.cut(p2, l1)
+    @test all(GO.equals.(
+        cut_polys,
+        [
+            GI.Polygon([[(0.0, 0.0), (5.0, 0.0), (5.0, 2.0), (2.0, 2.0), (2.0, 4.0), (5.0, 4.0), (5.0, 10.0), (0.0, 10.0), (0.0, 0.0)]]),
+            GI.Polygon([[(5.0, 0.0), (10.0, 0.0), (10.0, 10.0), (5.0, 10.0), (5.0, 4.0), (8.0, 4.0), (8.0, 2.0), (5.0, 2.0), (5.0, 0.0)]]),
+        ],
+    ))
+    cut_polys = GO.cut(p2, l3)
+    @test all(GO.equals.(
+        cut_polys,
+        [
+            GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 6.0), (0.0, 6.0), (0.0, 0.0)], r2]),
+            GI.Polygon([[(0.0, 6.0), (10.0, 6.0), (10.0, 10.0), (0.0, 10.0), (0.0, 6.0)]]),
+        ],
+    ))
+
+    # Cut concave polygon into three pieces
+    cut_polys = GO.cut(p3, l3)
+    @test all(GO.equals.(
+        cut_polys,
+        [
+            GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (1.0, 6.0), (4.0, 6.0), (5.0, 5.0), (6.0, 6.0), (9.0, 6.0), (10.0, 5.0), (10.0, 0.0), (0.0, 0.0)]]),
+            GI.Polygon([[(1.0, 6.0), (2.5, 7.5), (4.0, 6.0), (1.0, 6.0)]]),
+            GI.Polygon([[(6.0, 6.0), (7.5, 7.5), (9.0, 6.0), (6.0, 6.0)]]),
+        ],
+    ))
+
+    # Line doesn't cut through polygon
+    cut_polys = GO.cut(p1, l4)
+    @test all(GO.equals.(cut_polys, [p1]))
+    cut_polys = GO.cut(p1, l5)
+    @test all(GO.equals.(cut_polys, [p1]))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -24,6 +24,7 @@ const GO = GeometryOps
     @testset "Distance" begin include("methods/distance.jl") end
     @testset "Equals" begin include("methods/equals.jl") end
     # Clipping
+    @testset "Cut" begin include("methods/clipping/cut.jl") end
     @testset "Difference" begin include("methods/clipping/difference.jl") end
     @testset "Intersection" begin include("methods/clipping/intersection.jl") end
     @testset "Union" begin include("methods/clipping/union.jl") end