Debug and add to IntersectingPolygons (#107)

* Debug and add to IntersectingPolygons * Update src/transformations/correction/intersecting_polygons.jl Co-authored-by: Anshul Singhvi <[email protected]> --------- Co-authored-by: Anshul Singhvi <[email protected]>
JuliaGeo · Apr 14, 2024 · 1d65f97 · 1d65f97
1 parent 6df27de
commit 1d65f97
Show file tree

Hide file tree

Showing 3 changed files with 138 additions and 19 deletions.
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
@@ -382,6 +382,7 @@ function _flag_ent_exit!(::GI.LinearRingTrait, poly, pt_list, delay_cross_f, del
     npts = length(pt_list)
     next_idx = start_idx < npts ? (start_idx + 1) : 1
     start_val = (pt_list[start_idx].point .+ pt_list[next_idx].point) ./ 2
+    start_idx = next_idx - 1  # reset for iterating below
     status = !_point_filled_curve_orientation(start_val, poly; in = true, on = false, out = false)
     # Loop over points and mark entry and exit status
     start_chain_idx = 0

diff --git a/src/transformations/correction/intersecting_polygons.jl b/src/transformations/correction/intersecting_polygons.jl
@@ -59,13 +59,86 @@ struct UnionIntersectingPolygons <: GeometryCorrection end
 application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait
 
 function (::UnionIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
-    union_multipoly = if GI.npolygon(multipoly) > 1
+    union_multipoly = tuples(multipoly)
+    n_polys = GI.npolygon(multipoly)
+    if n_polys > 1
+        keep_idx = trues(n_polys)  # keep track of sub-polygons to remove
         # Combine any sub-polygons that intersect
-        first_poly = GI.getpolygon(multipoly, 1)
-        exclude_first_poly = GI.MultiPolygon(collect(Iterators.drop(GI.getpolygon(multipoly), 1)))
-        GI.MultiPolygon(union(first_poly, exclude_first_poly; target = GI.PolygonTrait(), fix_multipoly = nothing))
-    else
-        tuples(multipoly)
+        for (curr_idx, _) in Iterators.filter(last, Iterators.enumerate(keep_idx))
+            curr_poly = union_multipoly.geom[curr_idx]
+            poly_disjoint = false
+            while !poly_disjoint
+                poly_disjoint = true  # assume current polygon is disjoint from others
+                for (next_idx, _) in Iterators.filter(last, Iterators.drop(Iterators.enumerate(keep_idx), curr_idx))
+                    next_poly = union_multipoly.geom[next_idx]
+                    if intersects(curr_poly, next_poly)  # if two polygons intersect
+                        new_polys = union(curr_poly, next_poly; target = GI.PolygonTrait())
+                        n_new_polys = length(new_polys)
+                        if n_new_polys == 1  # if polygons combined
+                            poly_disjoint = false
+                            union_multipoly.geom[curr_idx] = new_polys[1]
+                            curr_poly = union_multipoly.geom[curr_idx]
+                            keep_idx[next_idx] = false
+                        end
+                    end
+                end
+            end
+        end
+        keepat!(union_multipoly.geom, keep_idx)
     end
     return union_multipoly
+end
+
+"""
+    DiffIntersectingPolygons() <: GeometryCorrection
+This correction ensures that the polygons included in a multipolygon aren't intersecting.
+If any polygon's are intersecting, they will be made nonintersecting through the [`difference`](@ref) 
+operation to create a unique set of disjoint (other than potentially connections by a single point)
+polygons covering the same area.
+See also [`GeometryCorrection`](@ref), [`UnionIntersectingPolygons`](@ref).
+"""
+struct DiffIntersectingPolygons <: GeometryCorrection end
+
+application_level(::DiffIntersectingPolygons) = GI.MultiPolygonTrait
+
+function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
+    diff_multipoly = tuples(multipoly)
+    n_starting_polys = GI.npolygon(multipoly)
+    n_polys = n_starting_polys
+    if n_polys > 1
+        keep_idx = trues(n_polys)  # keep track of sub-polygons to remove
+        # Break apart any sub-polygons that intersect
+        for curr_idx in 1:n_starting_polys
+            !keep_idx[curr_idx] && continue
+            for next_idx in (curr_idx + 1):n_starting_polys
+                !keep_idx[next_idx] && continue
+                next_poly = diff_multipoly.geom[next_idx]
+                n_new_polys = 0
+                curr_pieces_added = (n_polys + 1):(n_polys + n_new_polys)
+                for curr_piece_idx in Iterators.flatten((curr_idx:curr_idx, curr_pieces_added))
+                    !keep_idx[curr_piece_idx] && continue
+                    curr_poly = diff_multipoly.geom[curr_piece_idx]
+                    if intersects(curr_poly, next_poly)  # if two polygons intersect
+                        new_polys = difference(curr_poly, next_poly; target = GI.PolygonTrait())
+                        n_new_pieces = length(new_polys) - 1
+                        if n_new_pieces < 0  # current polygon is covered by next_polygon
+                            keep_idx[curr_piece_idx] = false
+                            break
+                        elseif n_new_pieces ≥ 0
+                            diff_multipoly.geom[curr_piece_idx] = new_polys[1]
+                            curr_poly = diff_multipoly.geom[curr_piece_idx]
+                            if n_new_pieces > 0 # current polygon breaks into several pieces
+                                append!(diff_multipoly.geom, @view new_polys[2:end])
+                                append!(keep_idx, trues(n_new_pieces))
+                                n_new_polys += n_new_pieces
+                            end
+                        end
+                    end
+                end
+                n_polys += n_new_polys
+            end
+        end
+        keepat!(diff_multipoly.geom, keep_idx)
+    end
+    return diff_multipoly
 end
diff --git a/test/transformations/correction/intersecting_polygons.jl b/test/transformations/correction/intersecting_polygons.jl
@@ -10,25 +10,70 @@ p2 = GI.Polygon([[(1.0, 0.0), (2.0, 0.0), (2.0, 1.0), (1.0, 1.0), (1.0, 0.0)]])
 p3 = GI.Polygon([[(1.0, 0.0), (2.0, 0.0), (2.0, -1.0), (1.0, -1.0), (1.0, 0.0)]])
 # (p1, p4) -> polygons are completly disjoint (no holes)
 p4 = GI.Polygon([[(1.0, -1.0), (2.0, -1.0), (2.0, -2.0), (1.0, -2.0), (1.0, -1.0)]])
+# (p1, p5) -> polygons cut through one another
+p5 = GI.Polygon([[(1.0, -1.0), (2.0, -1.0), (2.0, 4.0), (1.0, 4.0), (1.0, -1.0)]])
 
 mp1 = GI.MultiPolygon([p1])
 mp2 = GI.MultiPolygon([p1, p1])
 mp3 = GI.MultiPolygon([p1, p2, p3])
 mp4 = GI.MultiPolygon([p1, p4])
+mp5 = GI.MultiPolygon([p1, p5])
+mp6 = GI.MultiPolygon([  # four interlocking polygons forming a hole
+    [[(-5.0, 10.0), (-5.0, 15.0), (15.0, 15.0), (15.0, 10.0), (-5.0, 10.0)]],
+    [[(-5.0, -5.0), (-5.0, 0.0), (15.0, 0.0), (15.0, -5.0), (-5.0, -5.0)]],
+    [[(10.0, -5.0), (10.0, 15.0), (15.0, 15.0), (15.0, -5.0), (10.0, -5.0)]],
+    [[(-5.0, -5.0), (-5.0, 15.0), (0.0, 15.0), (0.0, -5.0), (-5.0, -5.0)]],
+])
 
-fixed_mp1 = GO.UnionIntersectingPolygons()(mp1)
-@test GI.npolygon(fixed_mp1) == 1
-@test GO.equals(GI.getpolygon(fixed_mp1, 1), p1)
+union_fixed_mp1 = GO.UnionIntersectingPolygons()(mp1)
+@test GI.npolygon(union_fixed_mp1) == 1
+@test GO.equals(GI.getpolygon(union_fixed_mp1, 1), p1)
 
-fixed_mp2 = GO.UnionIntersectingPolygons()(mp2)
-@test GI.npolygon(fixed_mp2) == 1
-@test GO.equals(GI.getpolygon(fixed_mp2, 1), p1)
+diff_fixed_mp1 = GO.DiffIntersectingPolygons()(mp1)
+@test GO.equals(diff_fixed_mp1, union_fixed_mp1)
 
-fixed_mp3 = GO.UnionIntersectingPolygons()(mp3)
-@test GI.npolygon(fixed_mp3) == 1
-@test GO.coveredby(p1, fixed_mp3) && GO.coveredby(p2, fixed_mp3) && GO.coveredby(p3, fixed_mp3)
+union_fixed_mp2 = GO.UnionIntersectingPolygons()(mp2)
+@test GI.npolygon(union_fixed_mp2) == 1
+@test GO.equals(GI.getpolygon(union_fixed_mp2, 1), p1)
 
-fixed_mp4 = GO.UnionIntersectingPolygons()(mp4)
-@test GI.npolygon(fixed_mp4) == 2
-@test (GO.equals(GI.getpolygon(fixed_mp4, 1), p1) && GO.equals(GI.getpolygon(fixed_mp4, 2), p4)) ||
-    (GO.equals(GI.getpolygon(fixed_mp4, 2), p1) && GO.equals(GI.getpolygon(fixed_mp4, 1), p4))
+diff_fixed_mp2 = GO.DiffIntersectingPolygons()(mp2)
+@test GO.equals(diff_fixed_mp2, union_fixed_mp2)
+
+union_fixed_mp3 = GO.UnionIntersectingPolygons()(mp3)
+@test GI.npolygon(union_fixed_mp3) == 1
+@test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(mp3)))
+diff_polys = GO.difference(union_fixed_mp3, mp3; target = GI.PolygonTrait(), fix_multipoly = nothing)
+@test sum(GO.area, diff_polys; init = 0.0) == 0
+
+diff_fixed_mp3 = GO.DiffIntersectingPolygons()(mp3)
+@test GI.npolygon(diff_fixed_mp3) == 3
+@test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(diff_fixed_mp3)))
+
+union_fixed_mp4 = GO.UnionIntersectingPolygons()(mp4)
+@test GI.npolygon(union_fixed_mp4) == 2
+@test (GO.equals(GI.getpolygon(union_fixed_mp4, 1), p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 2), p4)) ||
+    (GO.equals(GI.getpolygon(union_fixed_mp4, 2), p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 1), p4))
+
+diff_fixed_mp4 = GO.DiffIntersectingPolygons()(mp4)
+@test GO.equals(diff_fixed_mp4, union_fixed_mp4)
+
+union_fixed_mp5 = GO.UnionIntersectingPolygons()(mp5)
+@test GI.npolygon(union_fixed_mp5) == 1
+@test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(mp5)))
+diff_polys = GO.difference(union_fixed_mp5, mp5; target = GI.PolygonTrait(), fix_multipoly = nothing)
+@test sum(GO.area, diff_polys; init = 0.0) == 0
+
+diff_fixed_mp5 = GO.DiffIntersectingPolygons()(mp5)
+@test GI.npolygon(diff_fixed_mp5) == 3
+@test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(diff_fixed_mp5)))
+
+union_fixed_mp6 = GO.UnionIntersectingPolygons()(mp6)
+@test GI.npolygon(union_fixed_mp6) == 1
+@test GI.nhole(GI.getpolygon(union_fixed_mp6, 1)) == 1
+@test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(mp6)))
+diff_polys = GO.difference(union_fixed_mp6, mp6; target = GI.PolygonTrait(), fix_multipoly = nothing)
+@test sum(GO.area, diff_polys; init = 0.0) == 0
+
+diff_fixed_mp6 = GO.DiffIntersectingPolygons()(mp6)
+@test GI.npolygon(diff_fixed_mp6) == 4
+@test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(diff_fixed_mp6)))