Skip to content

Commit 1d65f97

Browse files
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]>
1 parent 6df27de commit 1d65f97

File tree

3 files changed

+138
-19
lines changed

3 files changed

+138
-19
lines changed

src/methods/clipping/clipping_processor.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -382,6 +382,7 @@ function _flag_ent_exit!(::GI.LinearRingTrait, poly, pt_list, delay_cross_f, del
382382
npts = length(pt_list)
383383
next_idx = start_idx < npts ? (start_idx + 1) : 1
384384
start_val = (pt_list[start_idx].point .+ pt_list[next_idx].point) ./ 2
385+
start_idx = next_idx - 1 # reset for iterating below
385386
status = !_point_filled_curve_orientation(start_val, poly; in = true, on = false, out = false)
386387
# Loop over points and mark entry and exit status
387388
start_chain_idx = 0

src/transformations/correction/intersecting_polygons.jl

Lines changed: 79 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -59,13 +59,86 @@ struct UnionIntersectingPolygons <: GeometryCorrection end
5959
application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait
6060

6161
function (::UnionIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
62-
union_multipoly = if GI.npolygon(multipoly) > 1
62+
union_multipoly = tuples(multipoly)
63+
n_polys = GI.npolygon(multipoly)
64+
if n_polys > 1
65+
keep_idx = trues(n_polys) # keep track of sub-polygons to remove
6366
# Combine any sub-polygons that intersect
64-
first_poly = GI.getpolygon(multipoly, 1)
65-
exclude_first_poly = GI.MultiPolygon(collect(Iterators.drop(GI.getpolygon(multipoly), 1)))
66-
GI.MultiPolygon(union(first_poly, exclude_first_poly; target = GI.PolygonTrait(), fix_multipoly = nothing))
67-
else
68-
tuples(multipoly)
67+
for (curr_idx, _) in Iterators.filter(last, Iterators.enumerate(keep_idx))
68+
curr_poly = union_multipoly.geom[curr_idx]
69+
poly_disjoint = false
70+
while !poly_disjoint
71+
poly_disjoint = true # assume current polygon is disjoint from others
72+
for (next_idx, _) in Iterators.filter(last, Iterators.drop(Iterators.enumerate(keep_idx), curr_idx))
73+
next_poly = union_multipoly.geom[next_idx]
74+
if intersects(curr_poly, next_poly) # if two polygons intersect
75+
new_polys = union(curr_poly, next_poly; target = GI.PolygonTrait())
76+
n_new_polys = length(new_polys)
77+
if n_new_polys == 1 # if polygons combined
78+
poly_disjoint = false
79+
union_multipoly.geom[curr_idx] = new_polys[1]
80+
curr_poly = union_multipoly.geom[curr_idx]
81+
keep_idx[next_idx] = false
82+
end
83+
end
84+
end
85+
end
86+
end
87+
keepat!(union_multipoly.geom, keep_idx)
6988
end
7089
return union_multipoly
90+
end
91+
92+
"""
93+
DiffIntersectingPolygons() <: GeometryCorrection
94+
This correction ensures that the polygons included in a multipolygon aren't intersecting.
95+
If any polygon's are intersecting, they will be made nonintersecting through the [`difference`](@ref)
96+
operation to create a unique set of disjoint (other than potentially connections by a single point)
97+
polygons covering the same area.
98+
See also [`GeometryCorrection`](@ref), [`UnionIntersectingPolygons`](@ref).
99+
"""
100+
struct DiffIntersectingPolygons <: GeometryCorrection end
101+
102+
application_level(::DiffIntersectingPolygons) = GI.MultiPolygonTrait
103+
104+
function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
105+
diff_multipoly = tuples(multipoly)
106+
n_starting_polys = GI.npolygon(multipoly)
107+
n_polys = n_starting_polys
108+
if n_polys > 1
109+
keep_idx = trues(n_polys) # keep track of sub-polygons to remove
110+
# Break apart any sub-polygons that intersect
111+
for curr_idx in 1:n_starting_polys
112+
!keep_idx[curr_idx] && continue
113+
for next_idx in (curr_idx + 1):n_starting_polys
114+
!keep_idx[next_idx] && continue
115+
next_poly = diff_multipoly.geom[next_idx]
116+
n_new_polys = 0
117+
curr_pieces_added = (n_polys + 1):(n_polys + n_new_polys)
118+
for curr_piece_idx in Iterators.flatten((curr_idx:curr_idx, curr_pieces_added))
119+
!keep_idx[curr_piece_idx] && continue
120+
curr_poly = diff_multipoly.geom[curr_piece_idx]
121+
if intersects(curr_poly, next_poly) # if two polygons intersect
122+
new_polys = difference(curr_poly, next_poly; target = GI.PolygonTrait())
123+
n_new_pieces = length(new_polys) - 1
124+
if n_new_pieces < 0 # current polygon is covered by next_polygon
125+
keep_idx[curr_piece_idx] = false
126+
break
127+
elseif n_new_pieces 0
128+
diff_multipoly.geom[curr_piece_idx] = new_polys[1]
129+
curr_poly = diff_multipoly.geom[curr_piece_idx]
130+
if n_new_pieces > 0 # current polygon breaks into several pieces
131+
append!(diff_multipoly.geom, @view new_polys[2:end])
132+
append!(keep_idx, trues(n_new_pieces))
133+
n_new_polys += n_new_pieces
134+
end
135+
end
136+
end
137+
end
138+
n_polys += n_new_polys
139+
end
140+
end
141+
keepat!(diff_multipoly.geom, keep_idx)
142+
end
143+
return diff_multipoly
71144
end

test/transformations/correction/intersecting_polygons.jl

Lines changed: 58 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -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)]])
1010
p3 = GI.Polygon([[(1.0, 0.0), (2.0, 0.0), (2.0, -1.0), (1.0, -1.0), (1.0, 0.0)]])
1111
# (p1, p4) -> polygons are completly disjoint (no holes)
1212
p4 = GI.Polygon([[(1.0, -1.0), (2.0, -1.0), (2.0, -2.0), (1.0, -2.0), (1.0, -1.0)]])
13+
# (p1, p5) -> polygons cut through one another
14+
p5 = GI.Polygon([[(1.0, -1.0), (2.0, -1.0), (2.0, 4.0), (1.0, 4.0), (1.0, -1.0)]])
1315

1416
mp1 = GI.MultiPolygon([p1])
1517
mp2 = GI.MultiPolygon([p1, p1])
1618
mp3 = GI.MultiPolygon([p1, p2, p3])
1719
mp4 = GI.MultiPolygon([p1, p4])
20+
mp5 = GI.MultiPolygon([p1, p5])
21+
mp6 = GI.MultiPolygon([ # four interlocking polygons forming a hole
22+
[[(-5.0, 10.0), (-5.0, 15.0), (15.0, 15.0), (15.0, 10.0), (-5.0, 10.0)]],
23+
[[(-5.0, -5.0), (-5.0, 0.0), (15.0, 0.0), (15.0, -5.0), (-5.0, -5.0)]],
24+
[[(10.0, -5.0), (10.0, 15.0), (15.0, 15.0), (15.0, -5.0), (10.0, -5.0)]],
25+
[[(-5.0, -5.0), (-5.0, 15.0), (0.0, 15.0), (0.0, -5.0), (-5.0, -5.0)]],
26+
])
1827

19-
fixed_mp1 = GO.UnionIntersectingPolygons()(mp1)
20-
@test GI.npolygon(fixed_mp1) == 1
21-
@test GO.equals(GI.getpolygon(fixed_mp1, 1), p1)
28+
union_fixed_mp1 = GO.UnionIntersectingPolygons()(mp1)
29+
@test GI.npolygon(union_fixed_mp1) == 1
30+
@test GO.equals(GI.getpolygon(union_fixed_mp1, 1), p1)
2231

23-
fixed_mp2 = GO.UnionIntersectingPolygons()(mp2)
24-
@test GI.npolygon(fixed_mp2) == 1
25-
@test GO.equals(GI.getpolygon(fixed_mp2, 1), p1)
32+
diff_fixed_mp1 = GO.DiffIntersectingPolygons()(mp1)
33+
@test GO.equals(diff_fixed_mp1, union_fixed_mp1)
2634

27-
fixed_mp3 = GO.UnionIntersectingPolygons()(mp3)
28-
@test GI.npolygon(fixed_mp3) == 1
29-
@test GO.coveredby(p1, fixed_mp3) && GO.coveredby(p2, fixed_mp3) && GO.coveredby(p3, fixed_mp3)
35+
union_fixed_mp2 = GO.UnionIntersectingPolygons()(mp2)
36+
@test GI.npolygon(union_fixed_mp2) == 1
37+
@test GO.equals(GI.getpolygon(union_fixed_mp2, 1), p1)
3038

31-
fixed_mp4 = GO.UnionIntersectingPolygons()(mp4)
32-
@test GI.npolygon(fixed_mp4) == 2
33-
@test (GO.equals(GI.getpolygon(fixed_mp4, 1), p1) && GO.equals(GI.getpolygon(fixed_mp4, 2), p4)) ||
34-
(GO.equals(GI.getpolygon(fixed_mp4, 2), p1) && GO.equals(GI.getpolygon(fixed_mp4, 1), p4))
39+
diff_fixed_mp2 = GO.DiffIntersectingPolygons()(mp2)
40+
@test GO.equals(diff_fixed_mp2, union_fixed_mp2)
41+
42+
union_fixed_mp3 = GO.UnionIntersectingPolygons()(mp3)
43+
@test GI.npolygon(union_fixed_mp3) == 1
44+
@test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(mp3)))
45+
diff_polys = GO.difference(union_fixed_mp3, mp3; target = GI.PolygonTrait(), fix_multipoly = nothing)
46+
@test sum(GO.area, diff_polys; init = 0.0) == 0
47+
48+
diff_fixed_mp3 = GO.DiffIntersectingPolygons()(mp3)
49+
@test GI.npolygon(diff_fixed_mp3) == 3
50+
@test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(diff_fixed_mp3)))
51+
52+
union_fixed_mp4 = GO.UnionIntersectingPolygons()(mp4)
53+
@test GI.npolygon(union_fixed_mp4) == 2
54+
@test (GO.equals(GI.getpolygon(union_fixed_mp4, 1), p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 2), p4)) ||
55+
(GO.equals(GI.getpolygon(union_fixed_mp4, 2), p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 1), p4))
56+
57+
diff_fixed_mp4 = GO.DiffIntersectingPolygons()(mp4)
58+
@test GO.equals(diff_fixed_mp4, union_fixed_mp4)
59+
60+
union_fixed_mp5 = GO.UnionIntersectingPolygons()(mp5)
61+
@test GI.npolygon(union_fixed_mp5) == 1
62+
@test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(mp5)))
63+
diff_polys = GO.difference(union_fixed_mp5, mp5; target = GI.PolygonTrait(), fix_multipoly = nothing)
64+
@test sum(GO.area, diff_polys; init = 0.0) == 0
65+
66+
diff_fixed_mp5 = GO.DiffIntersectingPolygons()(mp5)
67+
@test GI.npolygon(diff_fixed_mp5) == 3
68+
@test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(diff_fixed_mp5)))
69+
70+
union_fixed_mp6 = GO.UnionIntersectingPolygons()(mp6)
71+
@test GI.npolygon(union_fixed_mp6) == 1
72+
@test GI.nhole(GI.getpolygon(union_fixed_mp6, 1)) == 1
73+
@test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(mp6)))
74+
diff_polys = GO.difference(union_fixed_mp6, mp6; target = GI.PolygonTrait(), fix_multipoly = nothing)
75+
@test sum(GO.area, diff_polys; init = 0.0) == 0
76+
77+
diff_fixed_mp6 = GO.DiffIntersectingPolygons()(mp6)
78+
@test GI.npolygon(diff_fixed_mp6) == 4
79+
@test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(diff_fixed_mp6)))

0 commit comments

Comments
 (0)