diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index ac3d3e06d..d66d80fff 100644 --- 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 index dc5888a7a..e8858d45f 100644 --- 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 \ No newline at end of file diff --git a/test/transformations/correction/intersecting_polygons.jl b/test/transformations/correction/intersecting_polygons.jl index 24b3d97dc..c55618b72 100644 --- 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)) \ No newline at end of file +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)))