Ll/polygon intersects (#53)

* Fix up intersection point base calculation

* Update intersects and add line tests

* Add more tests and debug intersects

* Add comments to point_in_poly

* adds mutable struct for node

* Remove CairoMakie

* implements 1st draft of intersection for polygons

* correctly inserts intersection points into b_list

* takes into account alpha vaks when adding intr pts

* returns intersection of polygons

* Update equals and overlaps

* simplifies loops, & tracking of unprocessed pts

* Remove use of findfirst for 1.6 compat

* Updated geom, multi-geom equality

* adds latest changes before pulling equals func

* adds compareGO_LG

* expands tests for intersections

* adds tolerance when comparing polygons

* Add files for new functions

* deals with merge on 10_27

* Start cleaning up geom on and in geom

* Uncomment polygon_in_polygon

* Add more geom in geom

* adds more polygon to tests & checks if inter empty

* adjusts intersects for new point_in_polygon

* Add more comments

* adds union and intersection

* fixes difference

* inter/diff/union working for polys with holes

* Moves polygon clipping functions into own folder

* cleans up difference and intersection tests

* fix to_points and add repeat_last argument

* Changes k from 30 to 4

* cleans up more functions

* Fix merge issues

* gets rid of edges_a and edges_a in build_a_list

* removes all edges_a and edges_b variables

* better spacing

* fixes out-of-memory error

* fixes/cleans up _flag_ent_exit helper func

* changes alpha values to tuples

* trying to fix _build_b_list2

* makes _build_b_list_2 works

* Remove inter_list from tracing functions

* Remove unneeded vectors

* Stylistic fixes

* Remove use of idx field for starting point

* Finish updating clipping_processor

* Cleanup exterior clipping functions

* Clean up trace intersection

* Finish swapping out all traces for helper functions

* Cleaned up intersection hole code

* Major clipping cleanup done

* Switch return type

* Finish documenting and changing typing

* Update doc tests

* Update clipping doc tests

* Update doc test output types

* Update doc test spacing

* Small cleanups

* Remove clipping test utils

* commit so that I can merge Skylar's changes

* wrote nested for-loop in one line in intersection.jl

* fixes "exists" typo

* changes assert to throw(ArgumentError)

* fixes doc string for intersection

* fixes docstring for difference

* Update src/methods/clipping/clipping_processor.jl

Co-authored-by: Rafael Schouten <[email protected]>

* uses continue to avoid deep nesting

* uses continue to avoid deep nesting

* adds type 'T' to union docstring

* fixes line 5 union docstring

* fixes and error with the continues

* makes PolyNode Immutable

* makes _tuple_points earlier

* Update types

* Remove extra .gitignore files

---------

Co-authored-by: Skylar A Gering <[email protected]>
Co-authored-by: Rafael Schouten <[email protected]>

src/GeometryOps.jl

@@ -28,6 +28,10 @@ include("methods/bools.jl")
 include("methods/centroid.jl")
 include("methods/distance.jl")
 include("methods/equals.jl")
+include("methods/clipping/clipping_processor.jl")
+include("methods/clipping/intersection.jl")
+include("methods/clipping/difference.jl")
+include("methods/clipping/union.jl")
 include("methods/geom_relations/contains.jl")
 include("methods/geom_relations/coveredby.jl")
 include("methods/geom_relations/covers.jl")

src/methods/clipping/clipping_processor.jl

@@ -0,0 +1,291 @@
+# # This file contains the shared helper functions forlyNode the polygon clipping functionalities.
+
+#= This is the struct that makes up a_list and b_list. Many values are only used if point is
+an intersection point (ipt). =#
+struct PolyNode{T <: AbstractFloat}
+    idx::Int           # If ipt, index of point in a_idx_list, else 0
+    point::Tuple{T,T}  # (x, y) values of given point
+    inter::Bool        # If ipt, true, else 0
+    neighbor::Int      # If ipt, index of equivalent point in a_list or b_list, else 0
+    ent_exit::Bool     # If ipt, true if enter and false if exit, else false
+    fracs::Tuple{T,T}  # If ipt, fractions along edges to ipt (a_frac, b_frac), else (0, 0)
+end
+
+#=
+    _build_ab_list(::Type{T}, poly_a, poly_b) -> (a_list, b_list, a_idx_list)
+
+This function takes in two polygon rings and calls '_build_a_list', '_build_b_list', and
+'_flag_ent_exit' in order to fully form a_list and b_list. The 'a_list' and 'b_list' that it
+returns are the fully updated vectors of PolyNodes that represent the rings 'poly_a' and
+'poly_b', respectively. This function also returns 'a_idx_list', which at its "ith" index
+stores the index in 'a_list' at which the "ith" intersection point lies.
+=#
+function _build_ab_list(::Type{T}, poly_a, poly_b) where T
+    # Make a list for nodes of each polygon
+    a_list, a_idx_list = _build_a_list(T, poly_a, poly_b)
+    b_list = _build_b_list(T, a_idx_list, a_list, poly_b)
+
+    # Flag the entry and exits
+    _flag_ent_exit!(poly_b, a_list)
+    _flag_ent_exit!(poly_a, b_list)
+
+    return a_list, b_list, a_idx_list
+end
+
+#=
+    _build_a_list(::Type{T}, poly_a, poly_b) -> (a_list, a_idx_list)
+
+This function take in two polygon rings and creates a vector of PolyNodes to represent
+poly_a, including its intersection points with poly_b. The information stored in each
+PolyNode is needed for clipping using the Greiner-Hormann clipping algorithm.
+    
+Note: After calling this function, a_list is not fully formed because the neighboring
+indicies of the intersection points in b_list still need to be updated. Also we still have
+not update the entry and exit flags for a_list.
+    
+The a_idx_list is a list of the indicies of intersection points in a_list. The value at
+index i of a_idx_list is the location in a_list where the ith intersection point lies.
+=#
+function _build_a_list(::Type{T}, poly_a, poly_b) where T
+    n_a_edges = _nedge(poly_a)
+    a_list = Vector{PolyNode{T}}(undef, n_a_edges)  # list of points in poly_a
+    a_idx_list = Vector{Int}()  # finds indices of intersection points in a_list
+    intr_count = 0  # number of intersection points found
+    a_count = 0  # number of points added to a_list
+    # Loop through points of poly_a
+    local a_pt1
+    for (i, a_p2) in enumerate(GI.getpoint(poly_a))
+        a_pt2 = (T(GI.x(a_p2)), T(GI.y(a_p2)))
+        if i <= 1
+            a_pt1 = a_pt2
+            continue
+        end
+        # Add the first point of the edge to the list of points in a_list
+        new_point = PolyNode(0, a_pt1, false, 0, false, (zero(T), zero(T)))
+        a_count += 1
+        _add!(a_list, a_count, new_point, n_a_edges)
+        # Find intersections with edges of poly_b
+        local b_pt1
+        prev_counter = intr_count
+        for (j, b_p2) in enumerate(GI.getpoint(poly_b))
+            b_pt2 = _tuple_point(b_p2)
+            if j <=1
+                b_pt1 = b_pt2
+                continue
+            end
+            int_pt, fracs = _intersection_point(T, (a_pt1, a_pt2), (b_pt1, b_pt2))
+            # if no intersection point, skip this edge
+            if !isnothing(int_pt) && all(0 .≤ fracs .≤ 1)
+                # Set neighbor field to b edge (j-1) to keep track of intersection
+                new_intr = PolyNode(intr_count, int_pt, true, j - 1, false, fracs)
+                a_count += 1
+                intr_count += 1
+                _add!(a_list, a_count, new_intr, n_a_edges)
+                push!(a_idx_list, a_count)
+            end
+            b_pt1 = b_pt2
+        end
+
+        # Order intersection points by placement along edge using fracs value
+        if prev_counter < intr_count
+            Δintrs = intr_count - prev_counter
+            inter_points = @view a_list[(a_count - Δintrs + 1):a_count]
+            sort!(inter_points, by = x -> x.fracs[1])
+            for (i, p) in enumerate(inter_points)
+                inter_points[i] = PolyNode(prev_counter + i, p.point, p.inter, p.neighbor, p.ent_exit, p.fracs)
+            end
+        end
+
+        a_pt1 = a_pt2
+    end
+    return a_list, a_idx_list
+end
+
+# Add value x at index i to given array - if list isn't long enough, push value to array
+function _add!(arr, i, x, l = length(arr))
+    if i <= l
+        arr[i] = x
+    else
+        push!(arr, x)
+    end
+    return
+end
+
+#=
+    _build_b_list(::Type{T}, a_idx_list, a_list, poly_b) -> b_list
+
+This function takes in the a_list and a_idx_list build in _build_a_list and poly_b and
+creates a vector of PolyNodes to represent poly_b. The information stored in each PolyNode
+is needed for clipping using the Greiner-Hormann clipping algorithm.
+    
+Note: after calling this function, b_list is not fully updated. The entry/exit flags still
+need to be updated. However, the neightbor value in a_list is now updated.
+=#
+function _build_b_list(::Type{T}, a_idx_list, a_list, poly_b) where T
+    # Sort intersection points by insertion order in b_list
+    sort!(a_idx_list, by = x-> a_list[x].neighbor + a_list[x].fracs[2])
+    # Initialize needed values and lists
+    n_b_edges = _nedge(poly_b)
+    n_intr_pts = length(a_idx_list)
+    b_list = Vector{PolyNode{T}}(undef, n_b_edges + n_intr_pts)
+    intr_curr = 1
+    b_count = 0
+    # Loop over points in poly_b and add each point and intersection point
+    for (i, p) in enumerate(GI.getpoint(poly_b))
+        (i == n_b_edges + 1) && break
+        b_count += 1
+        pt = (T(GI.x(p)), T(GI.y(p)))
+        b_list[b_count] = PolyNode(0, pt, false, 0, false, (zero(T), zero(T)))
+        if intr_curr ≤ n_intr_pts
+            curr_idx = a_idx_list[intr_curr]
+            curr_node = a_list[curr_idx]
+            while curr_node.neighbor == i  # Add all intersection points in current edge
+                b_count += 1
+                b_list[b_count] = PolyNode(curr_node.idx, curr_node.point, true, curr_idx, false, curr_node.fracs)
+                a_list[curr_idx] = PolyNode(curr_node.idx, curr_node.point, curr_node.inter, b_count, curr_node.ent_exit, curr_node.fracs)
+                curr_node = a_list[curr_idx]
+                intr_curr += 1
+                intr_curr > n_intr_pts && break
+                curr_idx = a_idx_list[intr_curr]
+                curr_node = a_list[curr_idx]
+            end
+        end
+    end
+    sort!(a_idx_list)  # return a_idx_list to order of points in a_list
+    return b_list
+end
+
+#=
+    _flag_ent_exit(poly_b, a_list)
+
+This function flags all the intersection points as either an 'entry' or 'exit' point in
+relation to the given polygon.
+=#
+function _flag_ent_exit!(poly, pt_list)
+    local status
+    for ii in eachindex(pt_list)
+        if ii == 1
+            status = !_point_filled_curve_orientation(
+                pt_list[ii].point, poly;
+                in = true, on = false, out = false
+            )
+        elseif pt_list[ii].inter
+            pt_list[ii] = PolyNode(pt_list[ii].idx, pt_list[ii].point, pt_list[ii].inter, pt_list[ii].neighbor, status, pt_list[ii].fracs)
+            status = !status
+        end
+    end
+    return
+end
+
+#=
+    _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step)::Vector{GI.Polygon}
+
+This function takes the outputs of _build_ab_list and traces the lists to determine which
+polygons are formed as described in Greiner and Hormann. The function f_step determines in
+which direction the lists are traced.  This function is different for intersection,
+difference, and union. f_step must take in two arguments: the most recent intersection
+node's entry/exit status and a boolean that is true if we are currently tracing a_list and
+false if we are tracing b_list. The functions used for each clipping operation are follows:
+    - Intersection: (x, y) -> x ? 1 : (-1)
+    - Difference: (x, y) -> (x ⊻ y) ? 1 : (-1)
+    - Union: (x, y) -> (x ⊻ y) ? 1 : (-1)
+
+A list of GeoInterface polygons is returned from this function. 
+=#
+function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step) where T
+    n_a_pts, n_b_pts = length(a_list), length(b_list)
+    n_intr_pts = length(a_idx_list)
+    return_polys = Vector{_get_poly_type(T)}(undef, 0)
+    # Keep track of number of processed intersection points
+    processed_pts = 0
+    while processed_pts < n_intr_pts
+        curr_list, curr_npoints = a_list, n_a_pts
+        on_a_list = true
+        # Find first unprocessed intersecting point in subject polygon
+        processed_pts += 1
+        first_idx = findnext(x -> x != 0, a_idx_list, processed_pts)
+        idx = a_idx_list[first_idx]
+        a_idx_list[first_idx] = 0
+        start_pt = a_list[idx]
+
+        # Set first point in polygon
+        curr = curr_list[idx]
+        pt_list = [curr.point]
+
+        curr_not_start = true
+        while curr_not_start
+            step = f_step(curr.ent_exit, on_a_list)
+            curr_not_intr = true
+            while curr_not_intr
+                # Traverse polygon either forwards or backwards
+                idx += step
+                idx = (idx > curr_npoints) ? mod(idx, curr_npoints) : idx
+                idx = (idx == 0) ? curr_npoints : idx
+
+                # Get current node and add to pt_list
+                curr = curr_list[idx]
+                push!(pt_list, curr.point)
+                if curr.inter 
+                    # Keep track of processed intersection points
+                    curr_not_start = curr != start_pt && curr != b_list[start_pt.neighbor]
+                    if curr_not_start
+                        processed_pts += 1
+                        a_idx_list[curr.idx] = 0
+                    end
+                    curr_not_intr = false
+                end
+            end
+
+            # Switch to next list and next point
+            curr_list, curr_npoints = on_a_list ? (b_list, n_b_pts) : (a_list, n_a_pts)
+            on_a_list = !on_a_list
+            idx = curr.neighbor
+            curr = curr_list[idx]
+        end
+        push!(return_polys, GI.Polygon([pt_list]))
+    end
+    return return_polys
+end
+
+# Get type of polygons that will be made
+# TODO: Increase type options
+_get_poly_type(::Type{T}) where T =
+    GI.Polygon{false, false, Vector{GI.LinearRing{false, false, Vector{Tuple{T, T}}, Nothing, Nothing}}, Nothing, Nothing}
+
+#=
+    _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator)
+
+The holes specified by the hole iterator are added to the polygons in the return_polys list.
+If this creates more polygon, they are added to the end of the list. If this removes
+polygons, they are removed from the list
+=#
+function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator) where T
+    n_polys = length(return_polys)
+    # Remove set of holes from all polygons
+    for i in 1:n_polys
+        n_new_per_poly = 0
+        for hole in hole_iterator  # loop through all holes
+            hole_poly = GI.Polygon([hole])
+            # loop through all pieces of original polygon (new pieces added to end of list)
+            for j in Iterators.flatten((i:i, (n_polys + 1):(n_polys + n_new_per_poly)))
+                if !isnothing(return_polys[j])
+                    new_polys = difference(return_polys[j], hole_poly, T; target = GI.PolygonTrait)
+                    n_new_polys = length(new_polys)
+                    if n_new_polys == 0  # hole covered whole polygon
+                        return_polys[j] = nothing
+                    else
+                        return_polys[j] = new_polys[1]  # replace original
+                        if n_new_polys > 1  # add any extra pieces
+                            append!(return_polys, @view new_polys[2:end])
+                            n_new_per_poly += n_new_polys - 1
+                        end
+                    end
+                end
+            end
+        end
+        n_polys += n_new_per_poly
+    end
+    # Remove all polygon that were marked for removal
+    filter!(!isnothing, return_polys)
+    return
+end

src/methods/clipping/difference.jl

@@ -0,0 +1,83 @@
+# #  Difference Polygon Clipping
+export difference
+
+
+"""
+    difference(geom_a, geom_b, [T::Type]; target::Type)
+
+Return the difference between two geometries as a list of geometries. Return an empty list
+if none are found. The type of the list will be constrained as much as possible given the
+input geometries. Furthermore, the user can provide a `taget` type as a keyword argument and
+a list of target geometries found in the difference will be returned. The user can also
+provide a float type that they would like the points of returned geometries to be. 
+
+## Example 
+
+```jldoctest
+import GeoInterface as GI, GeometryOps as GO
+
+poly1 = GI.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]])
+poly2 = GI.Polygon([[[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]])
+diff_poly = GO.difference(poly1, poly2; target = GI.PolygonTrait)
+GI.coordinates.(diff_poly)
+
+# output
+1-element Vector{Vector{Vector{Vector{Float64}}}}:
+ [[[6.5, 3.5], [5.0, 5.0], [0.0, 0.0], [5.0, -5.0], [6.5, -3.5], [3.0, 0.0], [6.5, 3.5]]]
+```
+"""
+function difference(
+    geom_a, geom_b, ::Type{T} = Float64; target::Type{Target} = Nothing,
+) where {T <: AbstractFloat, Target <: GI.AbstractTrait}
+    return _difference(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+end
+
+#= The 'difference' function returns the difference of two polygons as a list of polygons.
+The algorithm to determine the difference was adapted from "Efficient clipping of efficient
+polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
+function _difference(
+    ::Type{GI.PolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.PolygonTrait, poly_b,
+) where T
+    # Get the exterior of the polygons
+    ext_poly_a = GI.getexterior(poly_a)
+    ext_poly_b = GI.getexterior(poly_b)
+    # Find the difference of the exterior of the polygons
+    a_list, b_list, a_idx_list = _build_ab_list(T, ext_poly_a, ext_poly_b)
+    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> (x ⊻ y) ? 1 : (-1))
+    if isempty(polys)
+        if _point_filled_curve_orientation(b_list[1].point, ext_poly_a) == point_in
+            poly_a_b_hole = GI.Polygon([ext_poly_a, ext_poly_b])
+            push!(polys, poly_a_b_hole)
+        elseif _point_filled_curve_orientation(a_list[1].point, ext_poly_b) != point_in
+            # Two polygons don't intersect and are not contained in one another
+            push!(polys, GI.Polygon([ext_poly_a]))
+        end
+    end
+
+    # If the original polygons had holes, take that into account.
+    if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
+        _add_holes_to_polys!(T, polys, GI.gethole(poly_a))
+        for hole in GI.gethole(poly_b)
+            new_polys = intersection(GI.Polygon([hole]), poly_a, T; target = GI.PolygonTrait)
+            if length(new_polys) > 0
+                append!(polys, new_polys)
+            end
+        end
+    end
+    return polys
+end
+
+# Many type and target combos aren't implemented
+function _difference(
+    ::Type{Target}, ::Type{T},
+    trait_a::GI.AbstractTrait, geom_a,
+    trait_b::GI.AbstractTrait, geom_b,
+) where {Target, T}
+    @assert(
+        false,
+        "Intersection between $trait_a and $trait_b with target $Target isn't implemented yet.",
+    )
+    return nothing
+end