Skip to content

Commit 2eb7328

Browse files
authored
Improvements from the Foster Extension (#78)
* Glued edges for intersection * Glued edges finished * Add new tests * Add tests back in * Allow approx intersection comparison * Update finding trace starting point * Don't add repeated points * Clean up and comment clipping processor * Clean up clipping tests * Add comments to clipping helper functions * Small changes to reduce memory * Add back idx field * Remove ExactPredicates * Add clipping file suggestions * Make hole points tuples * Improve intersection comments * Start code to remove collinear points * Debugged intersection failure * Debug union hole errors * Test collinear and finish comments * Code cleanup
1 parent e834037 commit 2eb7328

File tree

9 files changed

+531
-197
lines changed

9 files changed

+531
-197
lines changed

Project.toml

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,6 @@ version = "0.1.1"
55

66
[deps]
77
CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
8-
ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
98
GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
109
GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
1110
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -19,7 +18,6 @@ GeometryOpsProjExt = "Proj"
1918

2019
[compat]
2120
CoordinateTransformations = "0.5, 0.6"
22-
ExactPredicates = "2"
2321
GeoInterface = "1.2"
2422
GeometryBasics = "0.4.7"
2523
LinearAlgebra = "1"
@@ -33,9 +31,9 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
3331
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
3432
GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
3533
GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
36-
Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
3734
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
3835
LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
36+
Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
3937
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
4038
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
4139

src/methods/clipping/clipping_processor.jl

Lines changed: 246 additions & 139 deletions
Large diffs are not rendered by default.

src/methods/clipping/cut.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,8 @@ function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line) where T
7979
end
8080
cut_polys = [GI.Polygon([c]) for c in cut_coords]
8181
# Add original polygon holes back in
82-
_add_holes_to_polys!(T, cut_polys, GI.gethole(poly))
82+
remove_idx = falses(length(cut_polys))
83+
_add_holes_to_polys!(T, cut_polys, GI.gethole(poly), remove_idx)
8384
return cut_polys
8485
end
8586

src/methods/clipping/difference.jl

Lines changed: 32 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -44,36 +44,61 @@ function _difference(
4444
ext_a = GI.getexterior(poly_a)
4545
ext_b = GI.getexterior(poly_b)
4646
# Find the difference of the exterior of the polygons
47-
a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b)
48-
polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> (x y) ? 1 : (-1))
47+
a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f)
48+
polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _diff_step)
4949
# if no crossing points, determine if either poly is inside of the other
5050
if isempty(polys)
5151
a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
5252
# add case for if they polygons are the same (all intersection points!)
5353
# add a find_first check to find first non-inter poly!
5454
if b_in_a && !a_in_b # b in a and can't be the same polygon
55-
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"
5655
poly_a_b_hole = GI.Polygon([tuples(ext_a), tuples(ext_b)])
5756
push!(polys, poly_a_b_hole)
5857
elseif !b_in_a && !a_in_b # polygons don't intersect
5958
push!(polys, tuples(poly_a))
6059
return polys
6160
end
6261
end
63-
62+
remove_idx = falses(length(polys))
6463
# If the original polygons had holes, take that into account.
65-
if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
66-
_add_holes_to_polys!(T, polys, GI.gethole(poly_a))
64+
if GI.nhole(poly_a) != 0
65+
_add_holes_to_polys!(T, polys, GI.gethole(poly_a), remove_idx)
66+
end
67+
if GI.nhole(poly_b) != 0
6768
for hole in GI.gethole(poly_b)
68-
new_polys = intersection(GI.Polygon([hole]), poly_a, T; target = GI.PolygonTrait)
69+
hole_poly = GI.Polygon(StaticArrays.SVector(hole))
70+
new_polys = intersection(hole_poly, poly_a, T; target = GI.PolygonTrait)
6971
if length(new_polys) > 0
7072
append!(polys, new_polys)
7173
end
7274
end
7375
end
76+
# Remove uneeded collinear points on same edge
77+
for p in polys
78+
_remove_collinear_points!(p, remove_idx)
79+
end
7480
return polys
7581
end
7682

83+
# # Helper functions for Differences with Greiner and Hormann Polygon Clipping
84+
85+
#= When marking the crossing status of a delayed crossing, the chain start point is crossing
86+
when the start point is a entry point and is a bouncing point when the start point is an
87+
exit point. The end of the chain has the opposite crossing / bouncing status. =#
88+
_diff_delay_cross_f(x) = (x, !x)
89+
#= When marking the crossing status of a delayed bouncing, the chain start and end points
90+
are crossing if the current polygon's adjacent edges are within the non-tracing polygon and
91+
we are tracing b_list or if the edges are outside and we are on a_list. Otherwise the
92+
endpoints are marked as crossing. x is a boolean representing if the edges are inside or
93+
outside of the polygon and y is a variable that is true if we are on a_list and false if we
94+
are on b_list. =#
95+
_diff_delay_bounce_f(x, y) = x y
96+
#= When tracing polygons, step forwards if the most recent intersection point was an entry
97+
point and we are currently tracing b_list or if it was an exit point and we are currently
98+
tracing a_list, else step backwards, where x is the entry/exit status and y is a variable
99+
that is true if we are on a_list and false if we are on b_list. =#
100+
_diff_step(x, y) = (x y) ? 1 : (-1)
101+
77102
# Many type and target combos aren't implemented
78103
function _difference(
79104
::TraitTarget{Target}, ::Type{T},

src/methods/clipping/intersection.jl

Lines changed: 53 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,6 @@ _intersection(
4646
trait_b::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_b,
4747
) where T = _intersection_points(T, trait_a, geom_a, trait_b, geom_b)
4848

49-
5049
#= Polygon-Polygon Intersections with target Polygon
5150
The algorithm to determine the intersection was adapted from "Efficient clipping
5251
of efficient polygons," by Greiner and Hormann (1998).
@@ -60,8 +59,8 @@ function _intersection(
6059
ext_a = GI.getexterior(poly_a)
6160
ext_b = GI.getexterior(poly_b)
6261
# Then we find the intersection of the exteriors
63-
a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b)
64-
polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> x ? 1 : (-1))
62+
a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _inter_delay_cross_f, _inter_delay_bounce_f)
63+
polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _inter_step)
6564
if isempty(polys) # no crossing points, determine if either poly is inside the other
6665
a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
6766
if a_in_b
@@ -70,14 +69,35 @@ function _intersection(
7069
push!(polys, GI.Polygon([tuples(ext_b)]))
7170
end
7271
end
72+
remove_idx = falses(length(polys))
7373
# If the original polygons had holes, take that into account.
7474
if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
7575
hole_iterator = Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b)))
76-
_add_holes_to_polys!(T, polys, hole_iterator)
77-
end
76+
_add_holes_to_polys!(T, polys, hole_iterator, remove_idx)
77+
end
78+
# Remove uneeded collinear points on same edge
79+
for p in polys
80+
_remove_collinear_points!(p, remove_idx)
81+
end
7882
return polys
7983
end
8084

85+
# # Helper functions for Intersections with Greiner and Hormann Polygon Clipping
86+
87+
#= When marking the crossing status of a delayed crossing, the chain start point is bouncing
88+
when the start point is a entry point and is a crossing point when the start point is an
89+
exit point. The end of the chain has the opposite crossing / bouncing status. x is the
90+
entry/exit status. =#
91+
_inter_delay_cross_f(x) = (!x, x)
92+
#= When marking the crossing status of a delayed bouncing, the chain start and end points
93+
are crossing if the current polygon's adjacent edges are within the non-tracing polygon. If
94+
the edges are outside then the chain endpoints are marked as bouncing. x is a boolean
95+
representing if the edges are inside or outside of the polygon. =#
96+
_inter_delay_bounce_f(x, _) = x
97+
#= When tracing polygons, step forward if the most recent intersection point was an entry
98+
point, else step backwards where x is the entry/exit status. =#
99+
_inter_step(x, _) = x ? 1 : (-1)
100+
81101
# Many type and target combos aren't implemented
82102
function _intersection(
83103
::TraitTarget{Target}, ::Type{T},
@@ -175,17 +195,38 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge) where T
175195
# Second line runs from q to q + s
176196
qx, qy = GI.x(b1), GI.y(b1)
177197
sx, sy = GI.x(b2) - qx, GI.y(b2) - qy
178-
# Intersections will be where p + αr = q + βs where 0 < α, β < 1 and
198+
# Intersections will be where p + αr = q + βs where 0 < α, β < 1
179199
r_cross_s = rx * sy - ry * sx
180200
Δqp_x = qx - px
181201
Δqp_y = qy - py
182202
if r_cross_s != 0 # if lines aren't parallel
183-
α = (Δqp_x * sy - Δqp_y * sx) / r_cross_s
184-
β = (Δqp_x * ry - Δqp_y * rx) / r_cross_s
185-
x = px + α * rx
186-
y = py + α * ry
187-
if 0 α 1 && 0 β 1
188-
intr1 = (T(x), T(y)), (T(α), T(β))
203+
#= Calculate α = (Δqp_x * sy - Δqp_y * sx) / r_cross_s where we use approx
204+
comparisons due to inexact calculations =#
205+
α_num_t1, α_num_t2 = Δqp_x * sy, Δqp_y * sx # α numerator terms
206+
α_num = α_num_t1 - α_num_t2
207+
α = if α_num_t1 α_num_t2 # α = 0
208+
zero(T)
209+
elseif α_num r_cross_s # α = 1
210+
one(T)
211+
else # α != 0, 1
212+
T(α_num / r_cross_s)
213+
end
214+
#= Calculate β = (Δqp_x * ry - Δqp_y * rx) / r_cross_s where we use approx
215+
comparisons due to inexact calculations =#
216+
β_num_t1, β_num_t2 = Δqp_x * ry, Δqp_y * rx
217+
β_num = β_num_t1 - β_num_t2
218+
β = if β_num_t1 β_num_t2 # β = 0
219+
zero(T)
220+
elseif β_num r_cross_s # β = 1
221+
one(T)
222+
else # β != 0, 1
223+
T(β_num / r_cross_s)
224+
end
225+
# Calculate intersection point using α and β
226+
x = T(px + α * rx)
227+
y = T(py + α * ry)
228+
if 0 α 1 && 0 β 1 # only a valid intersection is 0 ≤ α, β ≤ 1
229+
intr1 = (x, y), (α, β)
189230
line_orient === 0 || α == 1 || β == 0 || β == 1) ? line_hinge : line_cross
190231
end
191232
elseif sx * Δqp_y == sy * Δqp_x # if parallel lines are collinear

src/methods/clipping/union.jl

Lines changed: 135 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -44,36 +44,158 @@ function _union(
4444
ext_a = GI.getexterior(poly_a)
4545
ext_b = GI.getexterior(poly_b)
4646
# Then, I get the union of the exteriors
47-
a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b)
48-
polys = _trace_polynodes(T, a_list, b_list, a_idx_list, (x, y) -> x ? (-1) : 1)
47+
a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _union_delay_cross_f, _union_delay_bounce_f)
48+
polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _union_step)
4949
n_pieces = length(polys)
50-
# Check if one polygon totally within other and if so, return the larger polygon.
50+
# Check if one polygon totally within other and if so, return the larger polygon
51+
a_in_b, b_in_a = false, false
5152
if n_pieces == 0 # no crossing points, determine if either poly is inside the other
5253
a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
5354
if a_in_b
5455
push!(polys, GI.Polygon([tuples(ext_b)]))
5556
elseif b_in_a
5657
push!(polys, GI.Polygon([tuples(ext_a)]))
5758
else
58-
share_edge_warn(a_list, "Edge case: polygons share edge but can't be combined.") # will get taken care of with "glued edges"
5959
push!(polys, tuples(poly_a))
6060
push!(polys, tuples(poly_b))
6161
return polys
6262
end
63-
elseif n_pieces > 1 # extra polygons are holes (n_pieces == 1 is the desired state)
64-
sort!(polys, by = area, rev = true) # sort so first element is the exterior
63+
elseif n_pieces > 1
64+
#= extra polygons are holes (n_pieces == 1 is the desired state) and since
65+
holes are formed by regions exterior to both poly_a and poly_b, they can't interact
66+
with pre-existing holes =#
67+
sort!(polys, by = area, rev = true) # sort by area so first element is the exterior
68+
# the first element is the exterior, the rest are holes
69+
@views append!(polys[1].geom, (GI.getexterior(p) for p in polys[2:end]))
70+
keepat!(polys, 1)
6571
end
66-
# the first element is the exterior, the rest are holes
67-
new_holes = @views (GI.getexterior(p) for p in polys[2:end])
68-
polys = n_pieces > 1 ? polys[1:1] : polys
69-
# Add holes back in for there are any
70-
if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0 || n_pieces > 1
71-
hole_iterator = Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b), new_holes))
72-
_add_holes_to_polys!(T, polys, hole_iterator)
72+
# Add in holes
73+
if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
74+
_add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b)
75+
end
76+
# Remove uneeded collinear points on same edge
77+
for p in polys
78+
_remove_collinear_points!(p, [false])
7379
end
7480
return polys
7581
end
7682

83+
# # Helper functions for Unions with Greiner and Hormann Polygon Clipping
84+
85+
#= When marking the crossing status of a delayed crossing, the chain start point is crossing
86+
when the start point is a entry point and is a bouncing point when the start point is an
87+
exit point. The end of the chain has the opposite crossing / bouncing status. =#
88+
_union_delay_cross_f(x) = (x, !x)
89+
90+
#= When marking the crossing status of a delayed bouncing, the chain start and end points
91+
are bouncing if the current polygon's adjacent edges are within the non-tracing polygon. If
92+
the edges are outside then the chain endpoints are marked as crossing. x is a boolean
93+
representing if the edges are inside or outside of the polygon. =#
94+
_union_delay_bounce_f(x, _) = !x
95+
96+
#= When tracing polygons, step backwards if the most recent intersection point was an entry
97+
point, else step forwards where x is the entry/exit status. =#
98+
_union_step(x, _) = x ? (-1) : 1
99+
100+
#= Add holes from two polygons to the exterior polygon formed by their union. If adding the
101+
the holes reveals that the polygons aren't actually intersecting, return the original
102+
polygons. =#
103+
function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b)
104+
if a_in_b
105+
_add_union_holes_contained_polys!(polys, poly_a, poly_b)
106+
elseif b_in_a
107+
_add_union_holes_contained_polys!(polys, poly_b, poly_a)
108+
else # Polygons intersect, but neither is contained in the other
109+
n_a_holes = GI.nhole(poly_a)
110+
ext_poly_a = GI.Polygon(StaticArrays.SVector(GI.getexterior(poly_a)))
111+
ext_poly_b = GI.Polygon(StaticArrays.SVector(GI.getexterior(poly_b)))
112+
#= Start with poly_b when comparing with holes from poly_a and then switch to poly_a
113+
to compare with holes from poly_b. For current_poly, use ext_poly_b to avoid
114+
repeating overlapping holes in poly_a and poly_b =#
115+
curr_exterior_poly = n_a_holes > 0 ? ext_poly_b : ext_poly_a
116+
current_poly = n_a_holes > 0 ? ext_poly_b : poly_a
117+
# Loop over all holes in both original polygons
118+
for (i, ih) in enumerate(Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b))))
119+
in_ext, _, _ = _line_polygon_interactions(ih, curr_exterior_poly; closed_line = true)
120+
if !in_ext
121+
#= if the hole isn't in the overlapping region between the two polygons, add
122+
the hole to the resulting polygon as we know it can't interact with any
123+
other holes =#
124+
push!(polys[1].geom, ih)
125+
else
126+
#= if the hole is at least partially in the overlapping region, take the
127+
difference of the hole from the polygon it didn't originate from - note that
128+
when current_poly is poly_a this includes poly_a holes so overlapping holes
129+
between poly_a and poly_b within the overlap are added, in addition to all
130+
holes in non-overlapping regions =#
131+
h_poly = GI.Polygon(StaticArrays.SVector(ih))
132+
new_holes = difference(h_poly, current_poly; target = GI.PolygonTrait())
133+
append!(polys[1].geom, (GI.getexterior(new_h) for new_h in new_holes))
134+
end
135+
if i == n_a_holes
136+
curr_exterior_poly = ext_poly_a
137+
current_poly = poly_a
138+
end
139+
end
140+
end
141+
return
142+
end
143+
144+
#= Add holes holes to the union of two polygons where one of the original polygons was
145+
inside of the other. If adding the the holes reveal that the polygons aren't actually
146+
intersecting, return the original polygons.=#
147+
function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly)
148+
union_poly = polys[1]
149+
ext_int_ring = GI.getexterior(interior_poly)
150+
for (i, ih) in enumerate(GI.gethole(exterior_poly))
151+
poly_ih = GI.Polygon(StaticArrays.SVector(ih))
152+
in_ih, on_ih, out_ih = _line_polygon_interactions(ext_int_ring, poly_ih; closed_line = true)
153+
if in_ih # at least part of interior polygon exterior is within the ith hole
154+
if !on_ih && !out_ih
155+
#= interior polygon is completly within the ith hole - polygons aren't
156+
touching and do not actually form a union =#
157+
polys[1] = tuples(interior_poly)
158+
push!(polys, tuples(exterior_poly))
159+
return polys
160+
else
161+
#= interior polygon is partially within the ith hole - area of interior
162+
polygon reduces the size of the hole =#
163+
new_holes = difference(poly_ih, interior_poly; target = GI.PolygonTrait())
164+
append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes))
165+
end
166+
else # none of interior polygon exterior is within the ith hole
167+
if !out_ih
168+
#= interior polygon's exterior is the same as the ith hole - polygons do
169+
form a union, but do not overlap so all holes stay in final polygon =#
170+
append!(union_poly.geom, Iterators.drop(GI.gethole(exterior_poly), i))
171+
append!(union_poly.geom, GI.gethole(interior_poly))
172+
return polys
173+
else
174+
#= interior polygon's exterior is outside of the ith hole - the interior
175+
polygon could either be disjoint from the hole, or contain the hole =#
176+
ext_int_poly = GI.Polygon(StaticArrays.SVector(ext_int_ring))
177+
in_int, _, _ = _line_polygon_interactions(ih, ext_int_poly; closed_line = true)
178+
if in_int
179+
#= interior polygon contains the hole - overlapping holes between the
180+
interior and exterior polygons will be added =#
181+
for jh in GI.gethole(interior_poly)
182+
poly_jh = GI.Polygon(StaticArrays.SVector(jh))
183+
if intersects(poly_ih, poly_jh)
184+
new_holes = intersection(poly_ih, poly_jh; target = GI.PolygonTrait())
185+
append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes))
186+
end
187+
end
188+
else
189+
#= interior polygon and the exterior polygon are disjoint - add the ith
190+
hole as it is not covered by the interior polygon =#
191+
push!(union_poly.geom, ih)
192+
end
193+
end
194+
end
195+
end
196+
return
197+
end
198+
77199

78200
# Many type and target combos aren't implemented
79201
function _union(

src/transformations/tuples.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -12,14 +12,14 @@ geometries wrapping `Tuple` points.
1212
1313
$APPLY_KEYWORDS
1414
"""
15-
function tuples(geom; kw...)
15+
function tuples(geom, ::Type{T} = Float64; kw...) where T
1616
if _is3d(geom)
1717
return apply(PointTrait(), geom; kw...) do p
18-
(Float64(GI.x(p)), Float64(GI.y(p)), Float64(GI.z(p)))
18+
(T(GI.x(p)), T(GI.y(p)), T(GI.z(p)))
1919
end
2020
else
2121
return apply(PointTrait(), geom; kw...) do p
22-
(Float64(GI.x(p)), Float64(GI.y(p)))
22+
(T(GI.x(p)), T(GI.y(p)))
2323
end
2424
end
2525
end

0 commit comments

Comments
 (0)