Clipping degeneracies (#58)

* Fix up intersection point base calculation

* Update intersects and add line tests

* Add more tests and debug intersects

* Add comments to point_in_poly

* Remove CairoMakie

* Update equals and overlaps

* Remove use of findfirst for 1.6 compat

* Updated geom, multi-geom equality

* Update frac calculations

* Fix non collinear intersection points

* Fix frac labeling

* Remove idx field

* Update intersection phase

* Add comments

* add helper functions for labelling degeneracies

* build_ab_list

* Update crossings and entry exit status

* adds tracining for simple intersections

* Remove unneeded list

* Fix entry exit tagging bug

* Update clipping edge cases

* Add entry/exit flagging for cutting

* Reorganize tests

* Fix test file

* Still debugging difference overflow

* Add comments for fix

* small adjustments in processing holes for union

* Fix add holes logic

* Consolidate clipping tests

* Turn all tests back on

* Fix merge remenants

* Fix inconsistent commenting

* add more test cases

* Remove try file

* Remove box

---------

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

src/methods/clipping/cut.jl

@@ -40,6 +40,8 @@ Return given geom cut by given line as a list of geometries of the same type as
 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.
 
+Note: This currently doesn't work for degenerate cases there line crosses through vertices.
+
 ## Example 
 
 ```jldoctest
@@ -70,7 +72,7 @@ function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line) where T
         return [tuples(poly)]
     end
     # Cut polygon by line
-    cut_coords = _cut(T, ext_poly, poly_list, intr_list, n_intr_pts)
+    cut_coords = _cut(T, ext_poly, line, poly_list, intr_list, n_intr_pts)
     # Close coords and create polygons
     for c in cut_coords
         push!(c, c[1])
@@ -94,10 +96,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, geom_list, intr_list, n_intr_pts) where T
+function _cut(::Type{T}, geom, line, 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)
+    _flag_ent_exit!(GI.LineTrait(), line, geom_list)
     # Add first point to output list
     return_coords = [[geom_list[1].point]]
     cross_backs = [(T(Inf),T(Inf))]

src/methods/clipping/difference.jl

@@ -41,18 +41,23 @@ function _difference(
     ::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)
+    ext_a = GI.getexterior(poly_a)
+    ext_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)
+    a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b)
     polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> (x ⊻ y) ? 1 : (-1))
+    # if no crossing points, determine if either poly is inside of the other
     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])
+        a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
+        # add case for if they polygons are the same (all intersection points!)
+        # add a find_first check to find first non-inter poly!
+        if b_in_a && !a_in_b  # b in a and can't be the same polygon
+            share_edge_warn(a_list, "Edge case: polygons share edge but one is hole of the other.")  # will get taken care of with "glued edges"
+            poly_a_b_hole = GI.Polygon([tuples(ext_a), tuples(ext_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]))
+        elseif !b_in_a && !a_in_b # polygons don't intersect
+            push!(polys, tuples(poly_a))
+            return polys
         end
     end
 

src/methods/clipping/intersection.jl

@@ -50,17 +50,17 @@ function _intersection(
     ::GI.PolygonTrait, poly_b,
 ) where {T}
     # First we get the exteriors of 'poly_a' and 'poly_b'
-    ext_poly_a = GI.getexterior(poly_a)
-    ext_poly_b = GI.getexterior(poly_b)
+    ext_a = GI.getexterior(poly_a)
+    ext_b = GI.getexterior(poly_b)
     # Then we find the intersection of the exteriors
-    a_list, b_list, a_idx_list = _build_ab_list(T, ext_poly_a, ext_poly_b)
+    a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b)
     polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> x ? 1 : (-1))
-
-    if isempty(polys)
-        if _point_filled_curve_orientation(a_list[1].point, ext_poly_b) == point_in
-            push!(polys, GI.Polygon([ext_poly_a]))
-        elseif _point_filled_curve_orientation(b_list[1].point, ext_poly_a) == point_in
-            push!(polys, GI.Polygon([ext_poly_b]))
+    if isempty(polys) # no crossing points, determine if either poly is inside the other
+        a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
+        if a_in_b
+            push!(polys, GI.Polygon([tuples(ext_a)]))
+        elseif b_in_a
+            push!(polys, GI.Polygon([tuples(ext_b)]))
         end
     end
     # If the original polygons had holes, take that into account.
@@ -161,14 +161,20 @@ function _intersection_point(::Type{T}, (a1, a2)::Tuple, (b1, b2)::Tuple) where
     sx, sy = GI.x(b2) - qx, GI.y(b2) - qy
     # Intersection will be where p + tr = q + us where 0 < t, u < 1 and
     r_cross_s = rx * sy - ry * sx
-    if r_cross_s != 0
-        Δqp_x = qx - px
-        Δqp_y = qy - py
+    Δqp_x = qx - px
+    Δqp_y = qy - py
+    point, fracs = if r_cross_s != 0  # if lines aren't parallel
         t = (Δqp_x * sy - Δqp_y * sx) / r_cross_s
         u = (Δqp_x * ry - Δqp_y * rx) / r_cross_s
         x = px + t * rx
         y = py + t * ry
-        return (T(x), T(y)), (T(t), T(u))
+        (T(x), T(y)), (T(t), T(u))
+    elseif  sx * Δqp_y == sy * Δqp_x  # if parallel lines are collinear
+        t = (Δqp_x * rx + Δqp_y * ry) / (rx^2 + ry^2)
+        u = -(Δqp_x * sx + Δqp_y * sy) / (sx^2 + sy^2)
+        nothing, (T(t), T(u))
+    else
+        nothing, nothing
     end
-    return nothing, nothing
+    return point, fracs
 end

src/methods/clipping/union.jl

@@ -41,42 +41,35 @@ function _union(
     ::GI.PolygonTrait, poly_b,
 ) where T
     # First, I get the exteriors of the two polygons
-    ext_poly_a = GI.getexterior(poly_a)
-    ext_poly_b = GI.getexterior(poly_b)
+    ext_a = GI.getexterior(poly_a)
+    ext_b = GI.getexterior(poly_b)
     # Then, I get the union of the exteriors
-    a_list, b_list, a_idx_list = _build_ab_list(T, ext_poly_a, ext_poly_b)
+    a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b)
     polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> x ? (-1) : 1)
+    n_pieces = length(polys)
     # Check if one polygon totally within other and if so, return the larger polygon.
-    if isempty(polys)
-        if _point_filled_curve_orientation(a_list[1].point, ext_poly_b) == point_in
-            push!(polys, GI.Polygon([ext_poly_b]))
-        elseif _point_filled_curve_orientation(b_list[1].point, ext_poly_a) == point_in
-            push!(polys,  GI.Polygon([ext_poly_a]))
+    if n_pieces == 0 # no crossing points, determine if either poly is inside the other
+        a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
+        if a_in_b
+            push!(polys, GI.Polygon([tuples(ext_b)]))
+        elseif b_in_a
+            push!(polys,  GI.Polygon([tuples(ext_a)]))
         else
-            push!(polys, poly_a)
-            push!(polys, poly_b)
+            share_edge_warn(a_list, "Edge case: polygons share edge but can't be combined.") # will get taken care of with "glued edges"
+            push!(polys, tuples(poly_a))
+            push!(polys, tuples(poly_b))
             return polys
         end
-    else
-        sort!(polys, by = area, rev = true)
-        polys = [GI.Polygon([GI.getexterior(p) for p in polys])]
+    elseif n_pieces > 1  # extra polygons are holes (n_pieces == 1 is the desired state)
+        sort!(polys, by = area, rev = true)  # sort so first element is the exterior
     end
-
-    n_b_holes = GI.nhole(poly_b)
-    if GI.nhole(poly_a) != 0 || n_b_holes != 0
-        new_poly = [GI.getexterior(polys[1]); collect(GI.gethole(polys[1]))]
-        current_poly = GI.Polygon([ext_poly_b])
-        for (i, hole) in enumerate(Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b))))
-            # Use ext_poly_b to not overcount overlapping holes in poly_a and in poly_b
-            new_hole = difference(GI.Polygon([hole]), current_poly, T; target = GI.PolygonTrait)
-            for h in new_hole
-                push!(new_poly, GI.getexterior(h))
-            end
-            if i == n_b_holes
-                current_poly = poly_a
-            end
-        end
-        polys[1] = GI.Polygon(new_poly)
+    # the first element is the exterior, the rest are holes
+    new_holes = @views (GI.getexterior(p) for p in polys[2:end])
+    polys = n_pieces > 1 ? polys[1:1] : polys
+    # Add holes back in for there are any
+    if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0 || n_pieces > 1
+        hole_iterator = Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b), new_holes))
+        _add_holes_to_polys!(T, polys, hole_iterator)
     end
     return polys
 end

src/methods/geom_relations/geom_geom_processors.jl

@@ -224,12 +224,12 @@ end
 #= Find where line and curve segments intersect by fraction of length. α is the fraction of
 the line (ls to le) and β is the traction of the curve (cs to ce). All inputs are tuples. =#
 function _find_intersect_fracs(ls, le, cs, ce)
-    _, fracs = _intersection_point(
+    point, fracs = _intersection_point(
         Float64,
         (ls, le),
         (cs, ce)
     )
-    (α, β) = if !isnothing(fracs)
+    (α, β) = if !isnothing(point)
         fracs
     else  # line and curve segments are parallel
         if equals(ls, cs)
@@ -507,8 +507,8 @@ function _point_filled_curve_orientation(
     n -= equals(GI.getpoint(curve, 1), GI.getpoint(curve, n)) ? 1 : 0
     k = 0  # counter for ray crossings
     p_start = GI.getpoint(curve, n)
-    @inbounds for i in 1:n
-        p_end = GI.getpoint(curve, i)
+    for (i, p_end) in enumerate(GI.getpoint(curve))
+        i > n && break
         v1 = GI.y(p_start) - y
         v2 = GI.y(p_end) - y
         if !((v1 < 0 && v2 < 0) || (v1 > 0 && v2 > 0)) # if not cases 11 or 26
@@ -671,7 +671,10 @@ function _line_filled_curve_interactions(
                             curve
                         )
                         npoints = length(ipoints)  # since hinge, at least one
-                        sort!(ipoints, by = p -> _euclid_distance(Float64, p, l_start))
+                        dist_from_lstart = let l_start = l_start
+                            x -> _euclid_distance(Float64, x, l_start)
+                        end
+                        sort!(ipoints, by = dist_from_lstart)
                         p_start = _tuple_point(l_start)
                         for i in 1:(npoints + 1)
                             p_end = i ≤ npoints ?