From 16a00ba76e05eaab7c78ed6cc549e645d1c99efd Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 11:47:36 -0800 Subject: [PATCH 01/18] Added generalized code --- src/GeometryOps.jl | 1 + .../geom_relations/geom_geom_processors.jl | 739 ++++++++++++++++++ 2 files changed, 740 insertions(+) create mode 100644 src/methods/geom_relations/geom_geom_processors.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 67bac7453..7feaf6474 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -30,6 +30,7 @@ include("methods/geom_relations/intersects.jl") include("methods/geom_relations/contains.jl") include("methods/geom_relations/crosses.jl") include("methods/geom_relations/disjoint.jl") +include("methods/geom_relations/geom_geom_processors.jl") include("methods/geom_relations/overlaps.jl") include("methods/geom_relations/within.jl") include("methods/polygonize.jl") diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl new file mode 100644 index 000000000..e1fee57db --- /dev/null +++ b/src/methods/geom_relations/geom_geom_processors.jl @@ -0,0 +1,739 @@ +#= Code is based off of DE-9IM Standards (https://en.wikipedia.org/wiki/DE-9IM) +and attempts a standardized solution for most of the functions. +=# + +@enum PointOrientation point_in=1 point_on=2 point_out=3 + +@enum LineOrientation line_cross=1 line_hinge=2 line_over=3 line_out=4 + + +#= +Determines if a point meets the given checks with respect to a curve. + +If in_allow is true, the point can be on the curve interior. +If on_allow is true, the point can be on the curve boundary. +If out_allow is true, the point can be disjoint from the curve. + +If the point is in an "allowed" location, return true. Else, return false. + +If closed_curve is true, curve is treated as a closed curve where the first and +last point are connected by a segment. +=# +function _point_curve_process( + point, curve; + in_allow, on_allow, out_allow, + closed_curve = false, +) + # Determine if curve is closed + n = GI.npoint(curve) + first_last_equal = equals(GI.getpoint(curve, 1), GI.getpoint(curve, n)) + closed_curve |= first_last_equal + n -= first_last_equal ? 1 : 0 + # Loop through all curve segments + p_start = GI.getpoint(curve, closed_curve ? n : 1) + @inbounds for i in (closed_curve ? 1 : 2):n + p_end = GI.getpoint(curve, i) + seg_val = _point_segment_orientation(point, p_start, p_end) + seg_val == point_in && return in_allow + if seg_val == point_on + if !closed_curve # if point is on curve endpoints, it is "on" + i == 2 && equals(point, p_start) && return on_allow + i == n && equals(point, p_end) && return on_allow + end + return in_allow + end + p_start = p_end + end + return out_allow +end + +#= +Determines if a point meets the given checks with respect to a polygon. + +If in_allow is true, the point can be within the polygon interior +If on_allow is true, the point can be on the polygon boundary. +If out_allow is true, the point can be disjoint from the polygon. + +If the point is in an "allowed" location, return true. Else, return false. +=# +function _point_polygon_process( + point, polygon; + in_allow, on_allow, out_allow, +) + # Check interaction of geom with polygon's exterior boundary + ext_val = _point_filled_curve_orientation(point, GI.getexterior(polygon)) + # If a point is outside, it isn't interacting with any holes + ext_val == point_out && return out_allow + # if a point is on an external boundary, it isn't interacting with any holes + ext_val == point_on && return on_allow + + # If geom is within the polygon, need to check interactions with holes + for hole in GI.gethole(polygon) + hole_val = _point_filled_curve_orientation(point, hole) + # If a point in in a hole, it is outside of the polygon + hole_val == point_in && return out_allow + # If a point in on a hole edge, it is on the edge of the polygon + hole_val == point_on && return on_allow + end + + # Point is within external boundary and on in/on any holes + return in_allow +end + +#= +Determines if a line meets the given checks with respect to a curve. + +If over_allow is true, segments of the line and curve can be co-linear. +If cross_allow is true, segments of the line and curve can cross. +If on_allow is true, endpoints of either the line or curve can intersect a + segment of the other geometry. +If cross_allow is true, segments of the line and curve can be disjoint. + +If in_require is true, the interiors of the line and curve must meet in at least + one point. +If on_require is true, the bounday of one of the two geometries can meet the + interior or boundary of the other geometry in at least one point. +If out_require is true, there must be at least one point of the given line that + is exterior of the curve. + +If the point is in an "allowed" location and meets all requirments, return true. +Else, return false. + +If closed_line is true, line is treated as a closed line where the first and +last point are connected by a segment. Same with closed_curve. +=# +function _line_curve_process( + line, curve; + over_allow, cross_allow, on_allow, out_allow, + in_require, on_require, out_require, + closed_line = false, + closed_curve = false, +) + # Set up requirments + in_req_met = !in_require + on_req_met = !on_require + out_req_met = !out_require + # Determine curve endpoints + nl = GI.npoint(line) + nc = GI.npoint(curve) + first_last_equal_line = equals(GI.getpoint(line, 1), GI.getpoint(line, nl)) + first_last_equal_curve = equals(GI.getpoint(curve, 1), GI.getpoint(curve, nc)) + nl -= first_last_equal_line ? 1 : 0 + nc -= first_last_equal_curve ? 1 : 0 + closed_line |= first_last_equal_line + closed_curve |= first_last_equal_curve + # Loop over each line segment + l_start = GI.getpoint(line, closed_line ? nl : 1) + i = closed_line ? 1 : 2 + while i ≤ nl + l_end = GI.getpoint(line, i) + c_start = GI.getpoint(curve, closed_curve ? nc : 1) + # Loop over each curve segment + for j in (closed_curve ? 1 : 2):nc + c_end = GI.getpoint(curve, j) + # Check if line and curve segments meet + seg_val = _segment_segment_orientation( + (l_start, l_end), + (c_start, c_end), + ) + # If segments are co-linear + if seg_val == line_over + !over_allow && return false + # at least one point in, meets requirments + in_req_met = true + point_val = _point_segment_orientation( + l_start, + c_start, c_end, + ) + # If entire segment isn't covered, consider remaining section + if point_val != point_out + if _point_segment_orientation( + l_end, + c_start, c_end, + ) != point_out + l_start = l_end + i += 1 + break + elseif _point_segment_orientation( + c_start, + l_start, l_end, + ) != point_out && !equals(l_start, c_start) + l_start = c_start + break + elseif _point_segment_orientation( + c_end, + l_start, l_end, + ) != point_out && !equals(l_start, c_end) + l_start = c_end + break + end + end + else + if seg_val == line_cross + !cross_allow && return false + in_req_met = true + elseif seg_val == line_hinge # could cross or overlap + # Determine location of intersection point on each segment + _, fracs = _intersection_point( + (_tuple_point(l_start), _tuple_point(l_end)), + (_tuple_point(c_start), _tuple_point(c_end)) + ) + (α, β) = + if !isnothing(fracs) + fracs + else # line and curve segments are parallel + if equals(l_start, c_start) + (0, 0) + elseif equals(l_start, c_end) + (0, 1) + elseif equals(l_end, c_start) + (1, 0) + else # equals(l_end, c_end) + (1, 1) + end + end + + if ( # don't consider edges of curves as they can't cross + (β == 0 && !closed_curve && j == 2) || + (β == 1 && !closed_curve && j == nc) || + (α == 0 && !closed_line && i == 2) || + (α == 1 && !closed_line && i == nl) + ) + !on_allow && return false + on_req_met = true + else + in_req_met = true + # If needed, determine if hinge actually crosses + if (!cross_allow || !over_allow) && α != 0 && β != 0 + l, c = if β == 1 + if α == 1 + ( + (l_end, GI.getpoint(line, i + 1)), + (c_end, GI.getpoint(curve, j + 1)) + ) + else + ( + (l_start, l_end), + (c_end, GI.getpoint(curve, j + 1)) + ) + end + else # β ≠ 1 and α == 1 + ( + (l_end, GI.getpoint(line, i + 1)), + (c_start, c_end) + ) + end + if _segment_segment_orientation(l, c) == line_hinge + !cross_allow && return false + else + !over_allow && return false + end + end + end + end + # no overlap for a give segment + if j == nc + !out_allow && return false + out_req_met = true + end + end + c_start = c_end + if j == nc + i += 1 + l_start = l_end + end + end + end + return in_req_met && on_req_met && out_req_met +end + +#= +Determines if a line meets the given checks with respect to a polygon. + +If in_allow is true, segments of the line can be in the polygon interior. +If on_allow is true, segments of the line can be on the polygon's boundary. +If out_allow is true, segments of the line can be outside of the polygon. + +If in_require is true, the interiors of the line and polygon must meet in at + least one point. +If on_require is true, the line must have at least one point on the polygon'same + boundary. +If out_require is true, the line must have at least one point outside of the + polygon. + +If the point is in an "allowed" location and meets all requirments, return true. +Else, return false. + +If closed_line is true, line is treated as a closed line where the first and +last point are connected by a segment. +=# +function _line_polygon_process( + line, polygon; + in_allow, on_allow, out_allow, + in_require, on_require, out_require, + closed_line = false, +) + in_req_met = !in_require + on_req_met = !on_require + out_req_met = !out_require + # Check interaction of line with polygon's exterior boundary + in_curve, on_curve, out_curve = _line_filled_curve_interactions( + line, GI.getexterior(polygon); + closed_line = closed_line, + ) + if on_curve + !on_allow && return false + on_req_met = true + end + if out_curve + !out_allow && return false + out_req_met = true + end + # If no points within the polygon, the line is disjoint and we are done + !in_curve && return in_req_met && on_req_met && out_req_met + + # Loop over polygon holes + for hole in GI.gethole(polygon) + in_hole, on_hole, out_hole =_line_filled_curve_interactions( + line, hole; + closed_line = closed_line, + ) + if in_hole # line in hole is equivalent to being out of polygon + !out_allow && return false + out_req_met = true + end + if on_hole # hole bounday is polygon boundary + !on_allow && return false + on_req_met = true + end + if !out_hole # entire line is in/on hole, can't be in/on other holes + in_curve = false + break + end + end + if in_curve # entirely of curve isn't within a hole + !in_allow && return false + in_req_met = true + end + return in_req_met && on_req_met && out_req_met +end + +#= +Determines if a polygon meets the given checks with respect to a polygon. + +If in_allow is true, the polygon's interiors must intersect. +If on_allow is true, the one of the polygon's boundaries must either interact + with the other polygon's boundary or interior. +If out_allow is true, the first polygon must have interior regions outside of + the second polygon. + +If in_require is true, the polygon interiors must meet in at least one point. +If on_require is true, one of the polygon's must have at least one boundary + point in or on the other polygon. +If out_require is true, the first polygon must have at least one interior point + outside of the second polygon. + +If the point is in an "allowed" location and meets all requirments, return true. +Else, return false. +=# +function _polygon_polygon_process( + poly1, poly2; + in_allow, on_allow, out_allow, + in_require, on_require, out_require, +) + in_req_met = !in_require + on_req_met = !on_require + out_req_met = !out_require + # Check if exterior of poly1 is within poly2 + ext1 = GI.getexterior(poly1) + ext2 = GI.getexterior(poly2) + # Check if exterior of poly1 is in polygon 2 + e1_in_p2, e1_on_p2, e1_out_p2 = _line_polygon_interactions( + ext1, poly2; + closed_line = true, + ) + if e1_on_p2 + !on_allow && return false + on_req_met = true + end + if e1_out_p2 + !out_allow && return false + out_req_met = true + end + + if !e1_in_p2 + # if exterior ring isn't in poly2, check if it surrounds poly2 + _, _, e2_out_e1 = _line_filled_curve_interactions( + ext2, ext1; + closed_line = true, + ) # if they really are disjoint, we are done + e2_out_e1 && return in_req_met && on_req_met && out_req_met + end + # If interiors interact, check if poly2 interacts with any of poly1's holes + for h1 in GI.gethole(poly1) + h1_in_p2, h1_on_p2, h1_out_p2 = _line_polygon_interactions( + h1, poly2; + closed_line = true, + ) + if h1_on_p2 + !on_allow && return false + on_req_met = true + end + if h1_out_p2 + !out_allow && return false + out_req_met = true + end + if !h1_in_p2 + # If hole isn't in poly2, see if poly2 is in hole + _, _, e2_out_h1 = _line_filled_curve_interactions( + ext2, h1; + closed_line = true, + ) + # hole encompasses all of poly2 + !e2_out_h1 && return in_req_met && on_req_met && out_req_met + break + end + end + #= + poly2 isn't outside of poly1 and isn't in a hole, poly1 interior must + interact with poly2 interior + =# + !in_allow && return false + in_req_met = true + + # If any of poly2 holes are within poly1, part of poly1 is exterior to poly2 + for h2 in GI.gethole(poly2) + h2_in_p1, h2_on_p1, _ = _line_polygon_interactions( + h2, poly1; + closed_line = true, + ) + if h2_on_p1 + !on_allow && return false + on_req_met = true + end + if h2_in_p1 + !out_allow && return false + out_req_met = true + end + end + return in_req_met && on_req_met && out_req_met +end + +#= +Determines if a point is in, on, or out of a segment. If the point is `on` the +segment it is on one of the segments endpoints. If it is `in`, it is on any +other point of the segment. If the point is not on any part of the segment, it +is `out` of the segment. + +Point should be an object of point trait and curve should be an object with a +linestring or linearring trait. + +Can provide values of in, on, and out keywords, which determines return values +for each scenario. +=# +function _point_segment_orientation( + point, start, stop; + in::T = point_in, on::T = point_on, out::T = point_out, +) where {T} + # Parse out points + x, y = GI.x(point), GI.y(point) + x1, y1 = GI.x(start), GI.y(start) + x2, y2 = GI.x(stop), GI.y(stop) + Δx_seg = x2 - x1 + Δy_seg = y2 - y1 + Δx_pt = x - x1 + Δy_pt = y - y1 + if (Δx_pt == 0 && Δy_pt == 0) || (Δx_pt == Δx_seg && Δy_pt == Δy_seg) + # If point is equal to the segment start or end points + return on + else + #= + Determine if the point is on the segment -> see if vector from segment + start to point is parallel to segment and if point is between the + segment endpoints + =# + on_line = _isparallel(Δx_seg, Δy_seg, Δx_pt, Δy_pt) + !on_line && return out + between_endpoints = + (x2 > x1 ? x1 <= x <= x2 : x2 <= x <= x1) && + (y2 > y1 ? y1 <= y <= y2 : y2 <= y <= y1) + !between_endpoints && return out + end + return in +end + +#= +Determine if point is in, on, or out of a closed curve, which includes the space +enclosed by the closed curve. + +`In` means the point is within the closed curve (excluding edges and vertices). +`On` means the point is on an edge or a vertex of the closed curve. +`Out` means the point is outside of the closed curve. + +Point should be an object of point trait and curve should be an object with a +linestring or linearring trait, that is assumed to be closed, regardless of +repeated last point. + +Can provide values of in, on, and out keywords, which determines return values +for each scenario. + +Note that this uses the Algorithm by Hao and Sun (2018): +https://doi.org/10.3390/sym10100477 +Paper seperates orientation of point and edge into 26 cases. For each case, it +is either a case where the point is on the edge (returns on), where a ray from +the point (x, y) to infinity along the line y = y cut through the edge (k += 1), +or the ray does not pass through the edge (do nothing and continue). If the ray +passes through an odd number of edges, it is within the curve, else outside of +of the curve if it didn't return 'on'. +See paper for more information on cases denoted in comments. +=# +function _point_filled_curve_orientation( + point, curve; + in::T = point_in, on::T = point_on, out::T = point_out, +) where {T} + x, y = GI.x(point), GI.y(point) + n = GI.npoint(curve) + 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) + 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 + u1 = GI.x(p_start) - x + u2 = GI.x(p_end) - x + c1 = u1 * v2 # first element of cross product summation + c2 = u2 * v1 # second element of cross product summation + f = c1 - c2 + if v2 > 0 && v1 ≤ 0 # Case 3, 9, 16, 21, 13, or 24 + (c1 ≈ c2) && return on # Case 16 or 21 + f > 0 && (k += 1) # Case 3 or 9 + elseif v1 > 0 && v2 ≤ 0 # Case 4, 10, 19, 20, 12, or 25 + (c1 ≈ c2) && return on # Case 19 or 20 + f < 0 && (k += 1) # Case 4 or 10 + elseif v2 == 0 && v1 < 0 # Case 7, 14, or 17 + (c1 ≈ c2) && return on # Case 17 + elseif v1 == 0 && v2 < 0 # Case 8, 15, or 18 + (c1 ≈ c2) && return on # Case 18 + elseif v1 == 0 && v2 == 0 # Case 1, 2, 5, 6, 22, or 23 + u2 ≤ 0 && u1 ≥ 0 && return on # Case 1 + u1 ≤ 0 && u2 ≥ 0 && return on # Case 2 + end + end + p_start = p_end + end + return iseven(k) ? out : in +end + +#= +Determines the type of interaction between two line segments. If the segments +`cross`, this means that they have a single intersection point that isn't on +either of their enpoints. If they form a `hinge`, they meet at one of the +segments endpoints. If they are `over`, then they are co-linear for at least +some of the length of the segments. Finally, if they are `out`, then the +segments are disjoint. + +Point should be an object of point trait and curve should be an object with a +linestring or linearring trait. + +Can provide values of in, on, and out keywords, which determines return values +for each scenario. +=# +function _segment_segment_orientation( + (a_point, b_point), (c_point, d_point); + cross::T = line_cross, hinge::T = line_hinge, + over::T = line_over, out::T = line_out, +) where T + (ax, ay) = _tuple_point(a_point) + (bx, by) = _tuple_point(b_point) + (cx, cy) = _tuple_point(c_point) + (dx, dy) = _tuple_point(d_point) + meet_type = ExactPredicates.meet((ax, ay), (bx, by), (cx, cy), (dx, dy)) + # Lines meet at one point within open segments + meet_type == 1 && return cross + # Lines don't meet at any points + meet_type == -1 && return out + # Lines meet at one or more points within closed segments + if _isparallel(((ax, ay), (bx, by)), ((cx, cy), (dx, dy))) + min_x, max_x = cx < dx ? (cx, dx) : (dx, cx) + min_y, max_y = cy < dy ? (cy, dy) : (dy, cy) + if ( + ((ax ≤ min_x && bx ≤ min_x) || (ax ≥ max_x && bx ≥ max_x)) && + ((ay ≤ min_y && by ≤ min_y) || (ay ≥ max_y && by ≥ max_y)) + ) + # a_point and b_point are on the same side of segment, don't overlap + return hinge + else + return over + end + end + # if lines aren't parallel then they must hinge + return hinge +end + +#= +Determines the types of interactions of a line with a filled-in curve. By +filled-in curve, I am referring to the exterior ring of a poylgon, for example. + +Returns a tuple of booleans: (in_curve, on_curve, out_curve). + +If in_curve is true, some of the lines interior points interact with the curve's + interior points. +If on_curve is true, endpoints of either the line intersect with the curve or + the line interacts with the polygon boundary. +If out_curve is true, at least one segments of the line is outside the curve. + +If closed_line is true, line is treated as a closed line where the first and +last point are connected by a segment. +=# +function _line_filled_curve_interactions( + line, curve; + closed_line = false, +) + in_curve = false + on_curve = false + out_curve = false + + # Determine number of points in curve and line + nl = GI.npoint(line) + nc = GI.npoint(curve) + first_last_equal_line = equals(GI.getpoint(line, 1), GI.getpoint(line, nl)) + first_last_equal_curve = equals(GI.getpoint(curve, 1), GI.getpoint(curve, nc)) + nl -= first_last_equal_line ? 1 : 0 + nc -= first_last_equal_curve ? 1 : 0 + closed_line |= first_last_equal_line + + # See if first point is in an acceptable orientation + l_start = GI.getpoint(line, closed_line ? nl : 1) + point_val = _point_filled_curve_orientation(l_start, curve) + if point_val == point_in + in_curve = true + elseif point_val == point_on + on_curve = true + else # point_val == point_out + out_curve = true + end + + # Check for any intersections between line and curve + for i in (closed_line ? 1 : 2):nl + l_end = GI.getpoint(line, i) + c_start = GI.getpoint(curve, nc) + # If already interacted with all regions of curve, can stop + in_curve && on_curve && out_curve && break + # Check next segment of line against curve + for j in 1:nc + c_end = GI.getpoint(curve, j) + # Check if two line and curve segments meet + seg_val = _segment_segment_orientation( + (l_start, l_end), + (c_start, c_end), + ) + if seg_val != line_out + # If line and curve meet, then at least one point is on boundary + on_curve = true + if seg_val == line_cross + # When crossing boundary, line is both in and out of curve + in_curve = true + out_curve = true + else + if seg_val == line_over + sp = _point_segment_orientation(l_start, c_start, c_end) + lp = _point_segment_orientation(l_end, c_start, c_end) + if sp != point_in || lp != point_in + #= + Line crosses over segment endpoint, creating a hinge + with another segment. + =# + seg_val = line_hinge + end + end + if seg_val == line_hinge + #= + Can't determine all types of interactions (in, out) with + hinge as it could pass through multiple other segments + so calculate if segment endpoints and intersections are + in/out of filled curve + =# + ipoints = intersection_points( + GI.Line([l_start, l_end]), + curve + ) + npoints = length(ipoints) # since hinge, at least one + sort!(ipoints, by = p -> _euclid_distance(Float64, p, l_start)) + p_start = _tuple_point(l_start) + for i in 1:(npoints + 1) + p_end = i ≤ npoints ? + ipoints[i] : + _tuple_point(l_end) + mid_val = _point_filled_curve_orientation( + (p_start .+ p_end) ./ 2, + curve, + ) + if mid_val == point_in + in_curve = true + elseif mid_val == point_out + out_curve = true + end + end + # already checked segment against whole filled curve + l_start = l_end + break + end + end + end + c_start = c_end + end + l_start = l_end + end + return in_curve, on_curve, out_curve +end + +#= +Determines the types of interactions of a line with a polygon. + +Returns a tuple of booleans: (in_poly, on_poly, out_poly). + +If in_poly is true, some of the lines interior points interact with the polygon + interior points. +If in_poly is true, endpoints of either the line intersect with the polygon or + the line interacts with the polygon boundary, including hole bounaries. +If out_curve is true, at least one segments of the line is outside the polygon, + including inside of holes. + +If closed_line is true, line is treated as a closed line where the first and +last point are connected by a segment. +=# +function _line_polygon_interactions( + line, polygon; + closed_line = false, +) + in_poly, on_poly, out_poly = _line_filled_curve_interactions( + line, GI.getexterior(polygon); + closed_line = closed_line, + ) + !in_poly && return (in_poly, on_poly, out_poly) + # Loop over polygon holes + for hole in GI.gethole(polygon) + in_hole, on_hole, out_hole =_line_filled_curve_interactions( + line, hole; + closed_line = closed_line, + ) + if in_hole + out_poly = true + end + if on_hole + on_poly = true + end + if !out_hole # entire line is in/on hole, can't be in/on other holes + in_poly = false + return (in_poly, on_poly, out_poly) + end + end + return in_poly, on_poly, out_poly +end + +function _point_in_extent(p, extent::Extents.Extent) + (x1, x2), (y1, y2) = extent.X, extent.Y + return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2 +end \ No newline at end of file From 3d44ea3518d39af98eb935c882955850499af53a Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 11:48:58 -0800 Subject: [PATCH 02/18] Update contains --- src/methods/geom_relations/contains.jl | 53 ++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 7 deletions(-) diff --git a/src/methods/geom_relations/contains.jl b/src/methods/geom_relations/contains.jl index a58e2a48b..1f7b2dd8d 100644 --- a/src/methods/geom_relations/contains.jl +++ b/src/methods/geom_relations/contains.jl @@ -1,13 +1,52 @@ -# # Containment +# # Contains export contains +#= +## What is contains? + +The contains function checks if a given geometry completly contains another +geometry, or in other words, that the second geometry is completly within the +first. This requires that the two interiors intersect and that the interior and +boundary of the second geometry is not in the exterior of the first geometry. + +To provide an example, consider these two lines: +```@example contains +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = GI.LineString([(0.0, 0.0), (1.0, 0.0), (0.0, 0.1)]) +l2 = GI.LineString([(0.25, 0.0), (0.75, 0.0)]) +f, a, p = lines(GI.getpoint(l1), color = :blue) +scatter!(GI.getpoint(l1), color = :blue) +lines!(GI.getpoint(l2), color = :orange) +scatter!(GI.getpoint(l2), color = :orange) +``` +We can see that all of the points and edges of l2 are within l1, so l1 contains +l2. However, l2 does not contain l1. +```@example contains +contains(l1, l2) # returns true +contains(l2, l1) # returns false +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +Given that contains is the exact opposite of within, we simply pass the two +inputs variables, swapped in order, to within. +=# + """ - contains(ft1::AbstractGeometry, ft2::AbstractGeometry)::Bool + contains(g1::AbstractGeometry, g2::AbstractGeometry)::Bool + +Return true if the second geometry is completely contained by the first +geometry. The interiors of both geometries must intersect and the interior and +boundary of the secondary (g2) must not intersect the exterior of the first +(g1). -Return true if the second geometry is completely contained by the first geometry. -The interiors of both geometries must intersect and, the interior and boundary of the secondary (geometry b) -must not intersect the exterior of the primary (geometry a). `contains` returns the exact opposite result of `within`. ## Examples @@ -15,11 +54,11 @@ must not intersect the exterior of the primary (geometry a). ```jldoctest import GeometryOps as GO, GeoInterface as GI line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) -point = (1, 2) +point = GI.Point((1, 2)) GO.contains(line, point) # output true ``` """ -contains(g1, g2)::Bool = within(g2, g1) +contains(g1, g2) = GeometryOps.within(g2, g1) \ No newline at end of file From c53144beb74fff16e3c837fc8d6211aa9be875b4 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 11:50:25 -0800 Subject: [PATCH 03/18] Add coveredby --- src/GeometryOps.jl | 3 +- src/methods/geom_relations/coveredby.jl | 265 ++++++++++++++++++++++++ 2 files changed, 267 insertions(+), 1 deletion(-) create mode 100644 src/methods/geom_relations/coveredby.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 7feaf6474..729ca7e08 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -26,11 +26,12 @@ include("methods/bools.jl") include("methods/centroid.jl") include("methods/distance.jl") include("methods/equals.jl") -include("methods/geom_relations/intersects.jl") include("methods/geom_relations/contains.jl") +include("methods/geom_relations/coveredby.jl") include("methods/geom_relations/crosses.jl") include("methods/geom_relations/disjoint.jl") include("methods/geom_relations/geom_geom_processors.jl") +include("methods/geom_relations/intersects.jl") include("methods/geom_relations/overlaps.jl") include("methods/geom_relations/within.jl") include("methods/polygonize.jl") diff --git a/src/methods/geom_relations/coveredby.jl b/src/methods/geom_relations/coveredby.jl new file mode 100644 index 000000000..e6ee52b15 --- /dev/null +++ b/src/methods/geom_relations/coveredby.jl @@ -0,0 +1,265 @@ +# # CoveredBy + +export coveredby + +#= +## What is coveredby? + +The coveredby function checks if one geometry is covered by another geometry. +This is an extension of within that does not require the interiors of the two +geometries to intersect, but still does require that the interior and boundary +of the first geometry isn't outside of the second geometry. + +To provide an example, consider this point and line: +```@example coveredby +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +p1 = Point(0.0, 0.0) +p2 = Point(1.0, 1.0) +l1 = Line(p1, p2) + +f, a, p = lines([p1, p2]) +scatter!(p1, color = :red) +``` +As we can see, `p1` is on the endpoint of l1. This means it is not `within`, but +it does meet the definition of `coveredby`. +```@example coveredby +coveredby(p1, l1) # true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`coveredby` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are allowed to be in the interior of g2 (either through + overlap or crossing for lines) + - points of g1 are allowed to be on the boundary of g2 + - points of g1 are not allowed to be in the exterior of g2 + - no points of g1 are required to be in the interior of g2 + - no points of g1 are required to be on the boundary of g2 + - no points of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# +const COVEREDBY_ALLOWS = (in_allow = true, on_allow = true, out_allow = false) +const COVEREDBY_CURVE_ALLOWS = (over_allow = true, cross_allow = true, on_allow = true, out_allow = false) +const COVEREDBY_CURVE_REQUIRES = (in_require = false, on_require = false, out_require = false) +const COVEREDBY_POLYGON_REQUIRES = (in_require = true, on_require = false, out_require = false,) + +""" + coveredby(g1, g2)::Bool + +Return `true` if the first geometry is completely covered by the second +geometry. The interior and boundary of the primary geometry (g1) must not +intersect the exterior of the secondary geometry (g2). + +Furthermore, `coveredby` returns the exact opposite result of `covers`. They are +equivalent with the order of the arguments swapped. + +## Examples +```jldoctest setup=:(using GeometryOps, GeometryBasics) +import GeometryOps as GO, GeoInterface as GI +p1 = GI.Point(0.0, 0.0) +p2 = GI.Point(1.0, 1.0) +l1 = GI.Line([p1, p2]) + +GO.coveredby(p1, l1) +# output +true +``` +""" +coveredby(g1, g2) = _coveredby(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries +_coveredby(::GI.FeatureTrait, g1, ::Any, g2) = coveredby(GI.geometry(g1), g2) +_coveredby(::Any, g1, t2::GI.FeatureTrait, g2) = coveredby(g1, GI.geometry(g2)) + + +# # Points coveredby geometries + +# Point is coveredby another point if those points are equal +_coveredby( + ::GI.PointTrait, g1, + ::GI.PointTrait, g2, +) = equals(g1, g2) + +# Point is coveredby a line/linestring if it is on a line vertex or an edge +_coveredby( + ::GI.PointTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _point_curve_process( + g1, g2; + COVEREDBY_ALLOWS..., + closed_curve = false, +) + +# Point is coveredby a linearring if it is on a vertex or an edge of ring +_coveredby( + ::GI.PointTrait, g1, + ::GI.LinearRingTrait, g2, +) = _point_curve_process( + g1, g2; + COVEREDBY_ALLOWS..., + closed_curve = true, +) + +# Point is coveredby a polygon if it is inside polygon, including edges/vertices +_coveredby( + ::GI.PointTrait, g1, + ::GI.PolygonTrait, g2, +) = _point_polygon_process( + g1, g2; + COVEREDBY_ALLOWS..., +) + +# Points cannot cover any geometry other than points +_coveredby( + ::Union{GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::GI.PointTrait, g2, +) = false + + +# # Lines coveredby geometries + +#= Linestring is coveredby a line if all interior and boundary points of the +first line are on the interior/boundary points of the second line. =# +_coveredby( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + COVEREDBY_CURVE_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = false, + closed_curve = false, +) + +#= Linestring is coveredby a ring if all interior and boundary points of the +line are on the edges of the ring. =# +_coveredby( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + COVEREDBY_CURVE_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring is coveredby a polygon if all interior and boundary points of the +line are in the polygon interior or on its edges, inlcuding hole edges. =# +_coveredby( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + COVEREDBY_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = false, +) + +# # Rings covered by geometries + +#= Linearring is covered by a line if all vertices and edges of the ring are on +the edges and vertices of the line. =# +_coveredby( + ::GI.LinearRingTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + COVEREDBY_CURVE_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = true, + closed_curve = false, +) + +#= Linearring is covered by another linear ring if all vertices and edges of the +first ring are on the edges/vertices of the second ring. =# +_coveredby( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + COVEREDBY_CURVE_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = true, + closed_curve = true, +) + +#= Linearring is coveredby a polygon if all vertices and edges of the ring are +in the polygon interior or on the polygon edges, inlcuding hole edges. =# +_coveredby( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + COVEREDBY_ALLOWS..., + COVEREDBY_CURVE_REQUIRES..., + closed_line = true, +) + + +# # Polygons covered by geometries + +#= Polygon is covered by another polygon if if the interior and edges of the +first polygon are in the second polygon interior or on polygon edges, including +hole edges.=# +_coveredby( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + COVEREDBY_ALLOWS..., + COVEREDBY_POLYGON_REQUIRES..., +) + +# Polygons cannot covered by any curves +_coveredby( + ::GI.PolygonTrait, g1, + ::GI.AbstractCurveTrait, g2, +) = false + + +# # Geometries coveredby multi-geometry/geometry collections + +#= Geometry is covered by a multi-geometry or a collection if one of the elements +of the collection cover the geometry. =# +function _coveredby( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + coveredby(g1, sub_g2) && return true + end + return false +end + +# # Multi-geometry/geometry collections coveredby geometries + +#= Multi-geometry or a geometry collection is covered by a geometry if all +elements of the collection are covered by the geometry. =# +function _coveredby( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + !coveredby(sub_g1, g2) && return false + end + return true +end \ No newline at end of file From dbc93acc54a6de48a6fa3cbdbf89b41a703e7131 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 11:51:25 -0800 Subject: [PATCH 04/18] Add covers --- src/GeometryOps.jl | 1 + src/methods/geom_relations/covers.jl | 63 ++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+) create mode 100644 src/methods/geom_relations/covers.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 729ca7e08..e6ca32e89 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -28,6 +28,7 @@ include("methods/distance.jl") include("methods/equals.jl") include("methods/geom_relations/contains.jl") include("methods/geom_relations/coveredby.jl") +include("methods/geom_relations/covers.jl") include("methods/geom_relations/crosses.jl") include("methods/geom_relations/disjoint.jl") include("methods/geom_relations/geom_geom_processors.jl") diff --git a/src/methods/geom_relations/covers.jl b/src/methods/geom_relations/covers.jl new file mode 100644 index 000000000..714d14643 --- /dev/null +++ b/src/methods/geom_relations/covers.jl @@ -0,0 +1,63 @@ +# # Covers + +export covers + +#= +## What is covers? + +The covers function checks if a given geometry completly covers another +geometry. For this to be true, the "contained" geometry's interior and +boundaries must be covered by the "covering" geometry's interior and boundaries. +The interiors do not need to overlap. + +To provide an example, consider these two lines: +```@example covers +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +p1 = Point(0.0, 0.0) +p2 = Point(1.0, 1.0) +l1 = Line(p1, p2) + +f, a, p = lines([p1, p2]) +scatter!(p1, color = :red) +``` + +```@example covers +covers(l1, p1) # returns true +covers(p1, l1) # returns false +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +Given that covers is the exact opposite of coveredby, we simply pass the two +inputs variables, swapped in order, to coveredby. +=# + +""" + covers(g1::AbstractGeometry, g2::AbstractGeometry)::Bool + +Return true if the first geometry is completely covers the second geometry, +The exterior and boundary of the second geometry must not be outside of the +interior and boundary of the first geometry. However, the interiors need not +intersect. + +`covers` returns the exact opposite result of `coveredby`. + +## Examples + +```jldoctest +import GeometryOps as GO, GeoInterface as GI +l1 = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) +l2 = GI.LineString([(1, 1), (1, 2)]) + +GO.covers(l1, l2) +# output +true +``` +""" +covers(g1, g2)::Bool = GeometryOps.coveredby(g2, g1) From 4620374c92a142d64be532e50076a739d47f081e Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 11:52:16 -0800 Subject: [PATCH 05/18] Update disjoint --- src/methods/geom_relations/disjoint.jl | 253 ++++++++++++++++++++++--- 1 file changed, 229 insertions(+), 24 deletions(-) diff --git a/src/methods/geom_relations/disjoint.jl b/src/methods/geom_relations/disjoint.jl index 02b7d4e46..f3e051c05 100644 --- a/src/methods/geom_relations/disjoint.jl +++ b/src/methods/geom_relations/disjoint.jl @@ -1,42 +1,247 @@ -# # Disjointness checks +# # Disjoint +export disjoint +#= +## What is disjoint? + +The disjoint function checks if one geometry is outside of another geometry, +without sharing any boundaries or interiors. + +To provide an example, consider these two lines: +```@example disjoint +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = GI.LineString([(0.0, 0.0), (1.0, 0.0), (0.0, 0.1)]) +l2 = GI.LineString([(2.0, 0.0), (2.75, 0.0)]) +f, a, p = lines(GI.getpoint(l1), color = :blue) +scatter!(GI.getpoint(l1), color = :blue) +lines!(GI.getpoint(l2), color = :orange) +scatter!(GI.getpoint(l2), color = :orange) +``` +We can see that none of the edges or vertices of l1 interact with l2 so they are +disjoint. +```@example disjoint +disjoint(l1, l2) # returns true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`disjoint` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are not allowed to be in the interior of g2 + - points of g1 are not allowed to be on the boundary of g2 + - points of g1 are allowed to be in the exterior of g2 + - no points required to be in the interior of g2 + - no points of g1 are required to be on the boundary of g2 + - no points of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# + +const DISJOINT_ALLOWS = (in_allow = false, on_allow = false, out_allow = true) +const DISJOINT_CURVE_ALLOWS = (over_allow = false, cross_allow = false, on_allow = false, out_allow = true) +const DISJOINT_REQUIRES = (in_require = false, on_require = false, out_require = false) """ disjoint(geom1, geom2)::Bool -Return `true` if the intersection of the two geometries is an empty set. +Return `true` if the first geometry is disjoint from the second geometry. -# Examples +Return `true` if the first geometry is disjoint from the second geometry. The +interiors and boundaries of both geometries must not intersect. -```jldoctest +## Examples +```jldoctest setup=:(using GeometryOps, GeometryBasics) import GeometryOps as GO, GeoInterface as GI -poly = GI.Polygon([[(-1, 2), (3, 2), (3, 3), (-1, 3), (-1, 2)]]) -point = (1, 1) -GO.disjoint(poly, point) +line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)]) +point = (2, 2) +GO.disjoint(point, line) # output true ``` """ -disjoint(g1, g2)::Bool = disjoint(trait(g1), g1, trait(g2), g2) +disjoint(g1, g2)::Bool = _disjoint(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries disjoint(::FeatureTrait, g1, ::Any, g2)::Bool = disjoint(GI.geometry(g1), g2) disjoint(::Any, g1, t2::FeatureTrait, g2)::Bool = disjoint(g1, geometry(g2)) -disjoint(::PointTrait, g1, ::PointTrait, g2)::Bool = !point_equals_point(g1, g2) -disjoint(::PointTrait, g1, ::LineStringTrait, g2)::Bool = !point_on_line(g1, g2) -disjoint(::PointTrait, g1, ::PolygonTrait, g2)::Bool = !point_in_polygon(g1, g2) -disjoint(::LineStringTrait, g1, ::PointTrait, g2)::Bool = !point_on_line(g2, g1) -disjoint(::LineStringTrait, g1, ::LineStringTrait, g2)::Bool = !line_on_line(g1, g2) -disjoint(::LineStringTrait, g1, ::PolygonTrait, g2)::Bool = !line_in_polygon(g2, g1) -disjoint(::PolygonTrait, g1, ::PointTrait, g2)::Bool = !point_in_polygon(g2, g1) -disjoint(::PolygonTrait, g1, ::LineStringTrait, g2)::Bool = !line_in_polygon(g2, g1) -disjoint(::PolygonTrait, g1, ::PolygonTrait, g2)::Bool = polygon_disjoint(g2, g1) - -function polygon_disjoint(poly1, poly2) - for point in GI.getpoint(poly1) - point_in_polygon(point, poly2) && return false + + +# # Point disjoint geometries + +# Point is disjoint from another point if the points are not equal. +_disjoint( + ::GI.PointTrait, g1, + ::GI.PointTrait, g2, +) = !equals(g1, g2) + +# Point is disjoint from a linestring if it is not on the line's edges/vertices. +_disjoint( + ::GI.PointTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _point_curve_process( + g1, g2; + DISJOINT_ALLOWS..., + closed_curve = false, +) + +# Point is disjoint from a linearring if it is not on the ring's edges/vertices. +_disjoint( + ::GI.PointTrait, g1, + ::GI.LinearRingTrait, g2, +) = _point_curve_process( + g1, g2; + DISJOINT_ALLOWS..., + closed_curve = true, +) + +#= Point is disjoint from a polygon if it is not on any edges, vertices, or +within the polygon's interior. =# +_disjoint( + ::GI.PointTrait, g1, + ::GI.PolygonTrait, g2, +) = _point_polygon_process( + g1, g2; + DISJOINT_ALLOWS..., +) + +#= Geometry is disjoint from a point if the point is not in the interior or on +the boundary of the geometry. =# +_disjoint( + trait1::Union{GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + trait2::GI.PointTrait, g2, +) = _disjoint(trait2, g2, trait1, g1) + + +# # Lines disjoint geometries + +#= Linestring is disjoint from another line if they do not share any interior +edge/vertex points or boundary points. =# +_disjoint( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + DISJOINT_CURVE_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = false, + closed_curve = false, +) + +#= Linestring is disjoint from a linearring if they do not share any interior +edge/vertex points or boundary points. =# +_disjoint( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + DISJOINT_CURVE_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring is disjoint from a polygon if the interior and boundary points of +the line are not in the polygon's interior or on the polygon's boundary. =# +_disjoint( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + DISJOINT_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = false, +) + +#= Geometry is disjoint from a linestring if the line's interior and boundary +points don't intersect with the geometrie's interior and boundary points. =# +_disjoint( + trait1::Union{GI.LinearRingTrait, GI.PolygonTrait}, g1, + trait2::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _disjoint(trait2, g2, trait1, g1) + + +# # Rings disjoint geometries + +#= Linearrings is disjoint from another linearring if they do not share any +interior edge/vertex points or boundary points.=# +_disjoint( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + DISJOINT_CURVE_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = true, + closed_curve = true, +) + +#= Linearring is disjoint from a polygon if the interior and boundary points of +the ring are not in the polygon's interior or on the polygon's boundary. =# +_disjoint( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + DISJOINT_ALLOWS..., + DISJOINT_REQUIRES..., + closed_line = true, +) + +# # Polygon disjoint geometries + +#= Polygon is disjoint from another polygon if they do not share any edges or +vertices and if their interiors do not intersect, excluding any holes. =# +_disjoint( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + DISJOINT_ALLOWS..., + DISJOINT_REQUIRES..., +) + + +# # Geometries disjoint multi-geometry/geometry collections + +#= Geometry is disjoint from a multi-geometry or a collection if all of the +elements of the collection are disjoint from the geometry. =# +function _disjoint( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + !disjoint(g1, sub_g2) && return false end - for point in GI.getpoint(poly2) - point_in_polygon(point, poly1) && return false + return true +end + +# # Multi-geometry/geometry collections coveredby geometries + +#= Multi-geometry or a geometry collection is covered by a geometry if all +elements of the collection are covered by the geometry. =# +function _disjoint( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + !disjoint(sub_g1, g2) && return false end - return !intersects(poly1, poly2) + return true end From f4d0faf5035785d79e98acd75a6210268c49b840 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 11:52:55 -0800 Subject: [PATCH 06/18] Update intersects --- src/methods/geom_relations/intersects.jl | 94 ++++-------------------- 1 file changed, 13 insertions(+), 81 deletions(-) diff --git a/src/methods/geom_relations/intersects.jl b/src/methods/geom_relations/intersects.jl index 2efaf1b78..834745a33 100644 --- a/src/methods/geom_relations/intersects.jl +++ b/src/methods/geom_relations/intersects.jl @@ -5,7 +5,10 @@ export intersects, intersection, intersection_points #= ## What is `intersects` vs `intersection` vs `intersection_points`? -The `intersects` methods check whether two geometries intersect with each other. +The intersects function checks if a given geometry intersects with another +geometry, or in other words, the either the interiors or boundaries of the two +geometries intersect. + The `intersection` methods return the geometry intersection between the two input geometries. The `intersection_points` method returns a list of intersection points between two geometries. @@ -44,18 +47,16 @@ f This is the GeoInterface-compatible implementation. -First, we implement a wrapper method for intersects, intersection, and -intersection_points that dispatches to the correct implementation based on the -geometry trait. The two underlying helper functions that are widely used in all -geometry dispatches are _line_intersects, which determines if two line segments -intersect and _intersection_point which determines the intersection point -between two line segments. +Given that intersects is the exact opposite of disjoint, we simply pass the two +inputs variables, swapped in order, to disjoint. =# """ intersects(geom1, geom2)::Bool -Check if two geometries intersect, returning true if so and false otherwise. +Return true if the interiors or boundaries of the two geometries interact. + +`intersects` returns the exact opposite result of `disjoint`. ## Example @@ -70,77 +71,8 @@ GO.intersects(line1, line2) true ``` """ -intersects(geom1, geom2) = intersects( - GI.trait(geom1), - geom1, - GI.trait(geom2), - geom2 -) - -""" - intersects(::GI.LineTrait, a, ::GI.LineTrait, b)::Bool - -Returns true if two line segments intersect and false otherwise. -""" -function intersects(::GI.LineTrait, a, ::GI.LineTrait, b) - a1 = _tuple_point(GI.getpoint(a, 1)) - a2 = _tuple_point(GI.getpoint(a, 2)) - b1 = _tuple_point(GI.getpoint(b, 1)) - b2 = _tuple_point(GI.getpoint(b, 2)) - meet_type = ExactPredicates.meet(a1, a2, b1, b2) - return meet_type == 0 || meet_type == 1 -end - -""" - intersects(::GI.AbstractTrait, a, ::GI.AbstractTrait, b)::Bool - -Returns true if two geometries intersect with one another and false -otherwise. For all geometries but lines, convert the geometry to a list of edges -and cross compare the edges for intersections. -""" -function intersects( - trait_a::GI.AbstractTrait, a_geom, - trait_b::GI.AbstractTrait, b_geom, -) edges_a, edges_b = map(sort! ∘ to_edges, (a_geom, b_geom)) - return _line_intersects(edges_a, edges_b) || - within(trait_a, a_geom, trait_b, b_geom) || - within(trait_b, b_geom, trait_a, a_geom) -end +intersects(geom1, geom2) = !disjoint(geom1, geom2) -""" - _line_intersects( - edges_a::Vector{Edge}, - edges_b::Vector{Edge} - )::Bool - -Returns true if there is at least one intersection between edges within the -two lists of edges. -""" -function _line_intersects( - edges_a::Vector{Edge}, - edges_b::Vector{Edge} -) - # Extents.intersects(to_extent(edges_a), to_extent(edges_b)) || return false - for edge_a in edges_a - for edge_b in edges_b - _line_intersects(edge_a, edge_b) && return true - end - end - return false -end - -""" - _line_intersects( - edge_a::Edge, - edge_b::Edge, - )::Bool - -Returns true if there is at least one intersection between two edges. -""" -function _line_intersects(edge_a::Edge, edge_b::Edge) - meet_type = ExactPredicates.meet(edge_a..., edge_b...) - return meet_type == 0 || meet_type == 1 -end """ intersection(geom_a, geom_b)::Union{Tuple{::Real, ::Real}, ::Nothing} @@ -229,8 +161,8 @@ Calculates the intersection between two line segments. Return nothing if there isn't one. """ function intersection( - trait_a::GI.AbstractTrait, geom_a, - trait_b::GI.AbstractTrait, geom_b, + trait_a::GI.AbstractGeometryTrait, geom_a, + trait_b::GI.AbstractGeometryTrait, geom_b, ) @assert( false, @@ -269,7 +201,7 @@ line segments, line strings, linear rings, polygons, and multipolygons. If no intersection points were possible given geometry extents, return nothing. If none are found, return an empty list. """ -function intersection_points(::GI.AbstractTrait, a, ::GI.AbstractTrait, b) +function intersection_points(::GI.AbstractGeometryTrait, a, ::GI.AbstractGeometryTrait, b) # Check if the geometries extents even overlap Extents.intersects(GI.extent(a), GI.extent(b)) || return nothing # Create a list of edges from the two input geometries From cca11bfab50cda40d53a5b4aac69e9016d340f8f Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 11:54:11 -0800 Subject: [PATCH 07/18] Add touches --- src/GeometryOps.jl | 1 + src/methods/geom_relations/touches.jl | 250 ++++++++++++++++++++++++++ 2 files changed, 251 insertions(+) create mode 100644 src/methods/geom_relations/touches.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index e6ca32e89..c94ff82d0 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -34,6 +34,7 @@ include("methods/geom_relations/disjoint.jl") include("methods/geom_relations/geom_geom_processors.jl") include("methods/geom_relations/intersects.jl") include("methods/geom_relations/overlaps.jl") +include("methods/geom_relations/touches.jl") include("methods/geom_relations/within.jl") include("methods/polygonize.jl") diff --git a/src/methods/geom_relations/touches.jl b/src/methods/geom_relations/touches.jl new file mode 100644 index 000000000..24bf4b2ea --- /dev/null +++ b/src/methods/geom_relations/touches.jl @@ -0,0 +1,250 @@ +# # Touches + +export touches + +#= +## What is touches? + +The touches function checks if one geometry touches another geometry. In other +words, the interiors of the two geometries don't interact, but one of the +geometries must have a boundary point that interacts with either the other +geometies interior or boundary. + + +To provide an example, consider these two lines: +```@example touches +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = Line([Point(0.0, 0.0), Point(1.0, 0.0)]) +l2 = Line([Point(1.0, 0.0), Point(1.0, -1.0)]) + +f, a, p = lines(l1) +lines!(l2) +``` +We can see that these two lines touch only at their endpoints. +```@example touches +touches(l1, l2) # true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`touches` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are not allowed to be in the interior of g2 + - points of g1 are allowed to be on the boundary of g2 + - points of g1 are allowed to be in the exterior of g2 + - no points of g1 are required to be in the interior of g2 + - at least one point of g1 is required to be on the boundary of g2 + - no points of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# + +const TOUCHES_POINT_ALLOWED = (in_allow = false, on_allow = true, out_allow = false) +const TOUCHES_CURVE_ALLOWED = (over_allow = false, cross_allow = false, on_allow = true, out_allow = true) +const TOUCHES_POLYGON_ALLOWS = (in_allow = false, on_allow = true, out_allow = true) +const TOUCHES_REQUIRES = (in_require = false, on_require = true, out_require = false) + +""" + touches(geom1, geom2)::Bool + +Return `true` if the first geometry touches the second geometry. In other words, +the two interiors cannot interact, but one of the geometries must have a +boundary point that interacts with either the other geometies interior or +boundary. + +## Examples +```jldoctest setup=:(using GeometryOps, GeometryBasics) +import GeometryOps as GO, GeoInterface as GI + +l1 = GI.Line([(0.0, 0.0), (1.0, 0.0)]) +l2 = GI.Line([(0.5, 1.0), (0.5, -1.0)]) + +GO.crosses(l1, l2) +# output +true +``` +""" +touches(g1, g2)::Bool = _touches(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries +_touches(::GI.FeatureTrait, g1, ::Any, g2) = touches(GI.geometry(g1), g2) +_touches(::Any, g1, t2::GI.FeatureTrait, g2) = touches(g1, GI.geometry(g2)) + + +# # Point touches geometries + +# Point cannot touch another point as if they are equal, interiors interact +_touches( + ::GI.PointTrait, g1, + ::GI.PointTrait, g2, +) = false + +# Point touches a linestring if it equal to the first of last point of the line +function _touches( + ::GI.PointTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) + n = GI.npoint(g2) + p1 = GI.getpoint(g2, 1) + pn = GI.getpoint(g2, n) + equals(p1, pn) && return false + return equals(g1, p1) || equals(g1, pn) +end + +# Point cannot 'touch' a linearring given that the ring has no boundary points +_touches( + ::GI.PointTrait, g1, + ::GI.LinearRingTrait, g2, +) = false + +# Point touches a polygon if it is on the boundary of that polygon +_touches( + ::GI.PointTrait, g1, + ::GI.PolygonTrait, g2, +) = _point_polygon_process( + g1, g2; + TOUCHES_POINT_ALLOWED..., +) + +#= Geometry touches a point if the point is on the geometry boundary. =# +_touches( + trait1::Union{GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + trait2::GI.PointTrait, g2, +) = _touches(trait2, g2, trait1, g1) + + +# # Lines touching geometries + +#= Linestring touches another line if at least one bounday point interacts with +the bounday of interior of the other line, but the interiors don't interact. =# +_touches( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + TOUCHES_CURVE_ALLOWED..., + TOUCHES_REQUIRES..., + closed_line = false, + closed_curve = false, +) + + +#= Linestring touches a linearring if at least one of the boundary points of the +line interacts with the linear ring, but their interiors can't interact. =# +_touches( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + TOUCHES_CURVE_ALLOWED..., + TOUCHES_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring touches a polygon if at least one of the boundary points of the +line interacts with the boundary of the polygon. =# +_touches( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + TOUCHES_POLYGON_ALLOWS..., + TOUCHES_REQUIRES..., + closed_line = false, +) + + +# # Rings touch geometries + +#= Linearring touches a linestring if at least one of the boundary points of the +line interacts with the linear ring, but their interiors can't interact. =# +_touches( + trait1::GI.LinearRingTrait, g1, + trait2::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _touches(trait2, g2, trait1, g1) + +#= Linearring cannot touch another linear ring since they are both exclusively +made up of interior points and no bounday points =# +_touches( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = false + +#= Linearring touches a polygon if at least one of the points of the ring +interact with the polygon bounday and non are in the polygon interior. =# +_touches( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + TOUCHES_POLYGON_ALLOWS..., + TOUCHES_REQUIRES..., + closed_line = true, +) + + +# # Polygons touch geometries + +#= Polygon touches a curve if at least one of the curve bounday points interacts +with the polygon's bounday and no curve points interact with the interior.=# +_touches( + trait1::GI.PolygonTrait, g1, + trait2::GI.AbstractCurveTrait, g2 +) = _touches(trait2, g2, trait1, g1) + + +#= Polygon touches another polygon if they share at least one boundary point and +no interior points. =# +_touches( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + TOUCHES_POLYGON_ALLOWS..., + TOUCHES_REQUIRES..., +) + +# # Geometries touch multi-geometry/geometry collections + +#= Geometry touch a multi-geometry or a collection if the geometry touches at +least one of the elements of the collection. =# +function _touches( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + !touches(g1, sub_g2) && return false + end + return true +end + +# # Multi-geometry/geometry collections cross geometries + +#= Multi-geometry or a geometry collection touches a geometry if at least one +elements of the collection touches the geometry. =# +function _touches( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + !touches(sub_g1, g2) && return false + end + return true +end \ No newline at end of file From 3bef37f123291478ef8095dc2b4ae6f0138c34e3 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 11:54:54 -0800 Subject: [PATCH 08/18] Update within --- src/methods/geom_relations/within.jl | 272 ++++++++++++++++++++++++--- 1 file changed, 249 insertions(+), 23 deletions(-) diff --git a/src/methods/geom_relations/within.jl b/src/methods/geom_relations/within.jl index 16366f944..9b28c0ee7 100644 --- a/src/methods/geom_relations/within.jl +++ b/src/methods/geom_relations/within.jl @@ -1,15 +1,70 @@ -# # Containment/withinness +# # Within export within +#= +## What is within? + +The within function checks if one geometry is inside another geometry. This +requires that the two interiors intersect and that the interior and +boundary of the first geometry is not in the exterior of the second geometry. + +To provide an example, consider these two lines: +```@example cshape +using GeometryOps +using GeometryOps.GeometryBasics +using Makie +using CairoMakie + +l1 = GI.LineString([(0.0, 0.0), (1.0, 0.0), (0.0, 0.1)]) +l2 = GI.LineString([(0.25, 0.0), (0.75, 0.0)]) +f, a, p = lines(GI.getpoint(l1), color = :blue) +scatter!(GI.getpoint(l1), color = :blue) +lines!(GI.getpoint(l2), color = :orange) +scatter!(GI.getpoint(l2), color = :orange) +``` +We can see that all of the points and edges of l2 are within l1, so l2 is +within l1, but l1 is not within l2 +```@example cshape +within(l1, l2) # returns false +within(l2, l1) # returns true +``` + +## Implementation + +This is the GeoInterface-compatible implementation. + +First, we implement a wrapper method that dispatches to the correct +implementation based on the geometry trait. + +Each of these calls a method in the geom_geom_processors file. The methods in +this file determine if the given geometries meet a set of criteria. For the +`within` function and arguments g1 and g2, this criteria is as follows: + - points of g1 are allowed to be in the interior of g2 (either through + overlap or crossing for lines) + - points of g1 are allowed to be on the boundary of g2 + - points of g1 are not allowed to be in the exterior of g2 + - at least one point of g1 is required to be in the interior of g2 + - no points of g1 are required to be on the boundary of g2 + - no points of g1 are required to be in the exterior of g2 + +The code for the specific implementations is in the geom_geom_processors file. +=# + +const WITHIN_POINT_ALLOWS = (in_allow = true, on_allow = false, out_allow = false) +const WITHIN_CURVE_ALLOWS = (over_allow = true, cross_allow = true, on_allow = true, out_allow = false) +const WITHIN_POLYGON_ALLOWS = (in_allow = true, on_allow = true, out_allow = false) +const WITHIN_REQUIRES = (in_require = true, on_require = false, out_require = false) """ - within(geom1, geom)::Bool + within(geom1, geom2)::Bool Return `true` if the first geometry is completely within the second geometry. -The interiors of both geometries must intersect and, the interior and boundary of the primary (geometry a) -must not intersect the exterior of the secondary (geometry b). -`within` returns the exact opposite result of `contains`. +The interiors of both geometries must intersect and the interior and boundary of +the primary geometry (geom1) must not intersect the exterior of the secondary +geometry (geom2). + +Furthermore, `within` returns the exact opposite result of `contains`. ## Examples ```jldoctest setup=:(using GeometryOps, GeometryBasics) @@ -23,21 +78,192 @@ GO.within(point, line) true ``` """ -# Syntactic sugar -within(g1, g2)::Bool = within(trait(g1), g1, trait(g2), g2)::Bool -within(::GI.FeatureTrait, g1, ::Any, g2)::Bool = within(GI.geometry(g1), g2) -within(::Any, g1, t2::GI.FeatureTrait, g2)::Bool = within(g1, GI.geometry(g2)) -# Points in geometries -within(::GI.PointTrait, g1, ::GI.LineStringTrait, g2)::Bool = point_on_line(g1, g2; ignore_end_vertices=true) -within(::GI.PointTrait, g1, ::GI.LinearRingTrait, g2)::Bool = point_on_line(g1, g2; ignore_end_vertices=true) -within(::GI.PointTrait, g1, ::GI.PolygonTrait, g2)::Bool = point_in_polygon(g1, g2; ignore_boundary=true) -# Lines in geometries -within(::GI.LineStringTrait, g1, ::GI.LineStringTrait, g2)::Bool = line_on_line(g1, g2) -within(::GI.LineStringTrait, g1, ::GI.LinearRingTrait, g2)::Bool = line_on_line(g1, g2) -within(::GI.LineStringTrait, g1, ::GI.PolygonTrait, g2)::Bool = line_in_polygon(g1, g2) -# Polygons within geometries -within(::GI.PolygonTrait, g1, ::GI.PolygonTrait, g2)::Bool = polygon_in_polygon(g1, g2) - -# Everything not specified -# TODO: Add multipolygons -within(::GI.AbstractTrait, g1, ::GI.AbstractCurveTrait, g2)::Bool = false \ No newline at end of file +within(g1, g2) = _within(trait(g1), g1, trait(g2), g2) + +# # Convert features to geometries +_within(::GI.FeatureTrait, g1, ::Any, g2) = within(GI.geometry(g1), g2) +_within(::Any, g1, t2::GI.FeatureTrait, g2) = within(g1, GI.geometry(g2)) + + +# # Points within geometries + +# Point is within another point if those points are equal. +_within( + ::GI.PointTrait, g1, + ::GI.PointTrait, g2, +) = equals(g1, g2) + +#= Point is within a linestring if it is on a vertex or an edge of that line, +excluding the start and end vertex if the line is not closed. =# +_within( + ::GI.PointTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _point_curve_process( + g1, g2; + WITHIN_POINT_ALLOWS..., + closed_curve = false, +) + +# Point is within a linearring if it is on a vertex or an edge of that ring. +_within( + ::GI.PointTrait, g1, + ::GI.LinearRingTrait, g2, +) = _point_curve_process( + g1, g2; + WITHIN_POINT_ALLOWS..., + closed_curve = true, +) + +#= Point is within a polygon if it is inside of that polygon, excluding edges, +vertices, and holes. =# +_within( + ::GI.PointTrait, g1, + ::GI.PolygonTrait, g2, +) = _point_polygon_process( + g1, g2; + WITHIN_POINT_ALLOWS..., +) + +# No geometries other than points can be within points +_within( + ::Union{GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::GI.PointTrait, g2, +) = false + + +# # Lines within geometries + +#= Linestring is within another linestring if their interiors intersect and no +points of the first line are in the exterior of the second line. =# +_within( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + WITHIN_CURVE_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = false, + closed_curve = false, +) + +#= Linestring is within a linear ring if their interiors intersect and no points +of the line are in the exterior of the ring. =# +_within( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + WITHIN_CURVE_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = false, + closed_curve = true, +) + +#= Linestring is within a polygon if their interiors intersect and no points of +the line are in the exterior of the polygon, although they can be on an edge. =# +_within( + ::Union{GI.LineTrait, GI.LineStringTrait}, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + WITHIN_POLYGON_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = false, +) + + +# # Rings covered by geometries + +#= Linearring is within a linestring if their interiors intersect and no points +of the ring are in the exterior of the line. =# +_within( + ::GI.LinearRingTrait, g1, + ::Union{GI.LineTrait, GI.LineStringTrait}, g2, +) = _line_curve_process( + g1, g2; + WITHIN_CURVE_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = true, + closed_curve = false, +) + +#= Linearring is within another linearring if their interiors intersect and no +points of the first ring are in the exterior of the second ring. =# +_within( + ::GI.LinearRingTrait, g1, + ::GI.LinearRingTrait, g2, +) = _line_curve_process( + g1, g2; + WITHIN_CURVE_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = true, + closed_curve = true, +) + +#= Linearring is within a polygon if their interiors intersect and no points of +the ring are in the exterior of the polygon, although they can be on an edge. =# +_within( + ::GI.LinearRingTrait, g1, + ::GI.PolygonTrait, g2, +) = _line_polygon_process( + g1, g2; + WITHIN_POLYGON_ALLOWS..., + WITHIN_REQUIRES..., + closed_line = true, +) + + +# # Polygons within geometries + +#= Polygon is within another polygon if the interior of the first polygon +intersects with the interior of the second and no points of the first polygon +are outside of the second polygon. =# +_within( + ::GI.PolygonTrait, g1, + ::GI.PolygonTrait, g2, +) = _polygon_polygon_process( + g1, g2; + WITHIN_POLYGON_ALLOWS..., + WITHIN_REQUIRES..., +) + +# Polygons cannot be within any curves +_within( + ::GI.PolygonTrait, g1, + ::GI.AbstractCurveTrait, g2, +) = false + + +# # Geometries within multi-geometry/geometry collections + +#= Geometry is within a multi-geometry or a collection if the geometry is within +at least one of the collection elements. =# +function _within( + ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g2, +) + for sub_g2 in GI.getgeom(g2) + within(g1, sub_g2) && return true + end + return false +end + +# # Multi-geometry/geometry collections within geometries + +#= Multi-geometry or a geometry collection is within a geometry if all +elements of the collection are within the geometry. =# +function _within( + ::Union{ + GI.MultiPointTrait, GI.AbstractMultiCurveTrait, + GI.MultiPolygonTrait, GI.GeometryCollectionTrait, + }, g1, + ::GI.AbstractGeometryTrait, g2, +) + for sub_g1 in GI.getgeom(g1) + !within(sub_g1, g2) && return false + end + return true +end From 63ebb0a26dc44a9d47838531d0ff88edfcfd683f Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 12:03:06 -0800 Subject: [PATCH 09/18] Update test file --- test/methods/geom_relations.jl | 362 +++++++++++++++------------------ 1 file changed, 169 insertions(+), 193 deletions(-) diff --git a/test/methods/geom_relations.jl b/test/methods/geom_relations.jl index 88e977866..815cd128d 100644 --- a/test/methods/geom_relations.jl +++ b/test/methods/geom_relations.jl @@ -1,151 +1,182 @@ -@testset "intersects" begin - @testset "Lines/Rings" begin - # Line test intersects ----------------------------------------------------- - - # Test for parallel lines - l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) - l2 = GI.Line([(0.0, 1.0), (2.5, 1.0)]) - @test !GO.intersects(l1, l2) - @test isnothing(GO.intersection(l1, l2)) - - # Test for non-parallel lines that don't intersect - l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) - l2 = GI.Line([(2.0, -3.0), (3.0, 0.0)]) - @test !GO.intersects(l1, l2) - @test isnothing(GO.intersection(l1, l2)) - - # Test for lines only touching at endpoint - l1 = GI.Line([(0.0, 0.0), (2.5, 0.0)]) - l2 = GI.Line([(2.0, -3.0), (2.5, 0.0)]) - @test GO.intersects(l1, l2) - @test all(GO.intersection(l1, l2) .≈ (2.5, 0.0)) +# Tests of DE-9IM Methods - # Test for lines that intersect in the middle - l1 = GI.Line([(0.0, 0.0), (5.0, 5.0)]) - l2 = GI.Line([(0.0, 5.0), (5.0, 0.0)]) - @test GO.intersects(l1, l2) - @test all(GO.intersection(l1, l2) .≈ (2.5, 2.5)) +pt1 = LG.Point([0.0, 0.0]) +pt2 = LG.Point([5.0, 5.0]) +pt3 = LG.Point([1.0, 0.0]) +pt4 = LG.Point([0.5, 0.0]) +pt5 = LG.Point([0.5, 0.25]) +pt6 = LG.Point([0.6, 0.4]) +pt7 = LG.Point([0.4, 0.8]) - # Line string test intersects ---------------------------------------------- +l1 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]]) +l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0]]) +l3 = LG.LineString([[-1.0, 0.0], [0.0, 0.0], [1.0, 0.0]]) +l4 = LG.LineString([[0.5, 0.0], [1.0, 0.0], [1.0, 0.5]]) +l5 = LG.LineString([[0.0, 0.0], [-1.0, -1.0]]) +l6 = LG.LineString([[2.0, 2.0], [0.0, 1.0]]) +l7 = LG.LineString([[0.5, 1.0], [0.5, -1.0]]) +l8 = LG.LineString([[0.0, 0.0], [0.5, 0.0], [0.5, 0.5], [1.0, -0.5]]) +l9 = LG.LineString([[0.0, 1.0], [0.0, -1.0], [1.0, 1.0]]) +l10 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 0.0]]) +l11 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 0.0], [-1.0, 0.0]]) +l12 = LG.LineString([[0.6, 0.5], [0.6, 0.9]]) +l13 = LG.LineString([[2.0, 2.0], [3.0, 3.0]]) +l14 = LG.LineString([[0.6, 0.25], [0.6, 0.35]]) +l15 = LG.LineString([[0.0, 3.0], [4.0, 3.0]]) +l16 = LG.LineString([[0.3, -0.7], [1.0, 0.0], [3.0, 0.6]]) - # Single element line strings crossing over each other - l1 = LG.LineString([[5.5, 7.2], [11.2, 12.7]]) - l2 = LG.LineString([[4.3, 13.3], [9.6, 8.1]]) - @test GO.intersects(l1, l2) - go_inter = GO.intersection(l1, l2) - lg_inter = LG.intersection(l1, l2) - @test go_inter[1][1] .≈ GI.x(lg_inter) - @test go_inter[1][2] .≈ GI.y(lg_inter) +r1 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 0.0]]) +r2 = LG.LinearRing([[0.5, 0.2], [0.6, 0.4], [0.7, 0.2], [0.5, 0.2]]) +r3 = LG.LinearRing([[0.2, 0.7], [0.4, 0.9], [0.5, 0.6], [0.2, 0.7]]) +r4 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]]) +r5 = LG.LinearRing([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [2.0, 1.0], [1.0, 2.0], [1.0, 1.0], [0.0, 0.0]]) +r6 = LG.LinearRing([[0.0, 0.0], [-1.0, 0.0], [-1.0, 1.0], [0.0, 1.0], [0.0, 0.0]]) +r7 = LG.LinearRing([[0.5, 0.5], [1.5, 0.5], [1.5, 1.5], [0.5, 1.5], [0.5, 0.5]]) +r8 = LG.LinearRing([[0.1, 0.1], [0.2, 0.1], [0.2, 0.2], [0.1, 0.2], [0.1, 0.1]]) +r9 = LG.LinearRing([[1.0, -0.5], [0.5, -1.0], [1.5, -1.0], [1.0, -0.5]]) +r10 = LG.LinearRing([[0.5, 0.2], [0.6, 0.4], [0.3, 0.3], [0.5, 0.2]]) +r11 = LG.LinearRing([[0.55, 0.21], [0.55, 0.23], [0.65, 0.23], [0.66, 0.21], [0.55, 0.21]]) - # Multi-element line strings crossing over on vertex - l1 = LG.LineString([[0.0, 0.0], [2.5, 0.0], [5.0, 0.0]]) - l2 = LG.LineString([[2.0, -3.0], [3.0, 0.0], [4.0, 3.0]]) - @test GO.intersects(l1, l2) - go_inter = GO.intersection(l1, l2) - @test length(go_inter) == 1 - lg_inter = LG.intersection(l1, l2) - @test go_inter[1][1] .≈ GI.x(lg_inter) - @test go_inter[1][2] .≈ GI.y(lg_inter) +p1 = LG.Polygon(r4, [r2, r3]) +p2 = LG.Polygon(r4) +p3 = LG.Polygon(r2) +p4 = LG.Polygon(r6) +p5 = LG.Polygon(r9) +p6 = LG.Polygon(r11) +p7 = LG.Polygon(r7) +p8 = LG.Polygon([[[0.0, 0.0], [0.0, 3.0], [1.0, 2.0], [2.0, 3.0], [3.0, 2.0], [4.0, 3.0], [4.0, 0.0], [0.0, 0.0]]]) +p9 = LG.Polygon([[[0.0, 0.0], [0.0, 3.0], [1.0, 4.0], [2.0, 3.0], [3.0, 4.0], [4.0, 3.0], [4.0, 0.0], [0.0, 0.0]]]) +p10 = LG.Polygon([ + [[0.1, 0.5], [0.1, 0.99], [0.6, 0.99], [0.6, 0.5], [0.1, 0.5]], + [[0.15, 0.55], [0.15, 0.95], [0.55, 0.95], [0.55, 0.55], [0.15, 0.55]] +]) +p11 = LG.Polygon(r3) - # Multi-element line strings crossing over with multiple intersections - l1 = LG.LineString([[0.0, -1.0], [1.0, 1.0], [2.0, -1.0], [3.0, 1.0]]) - l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [3.0, 0.0]]) - @test GO.intersects(l1, l2) - go_inter = GO.intersection(l1, l2) - @test length(go_inter) == 3 - lg_inter = LG.intersection(l1, l2) - @test issetequal( - Set(go_inter), - Set(GO._tuple_point.(GI.getpoint(lg_inter))) - ) - - # Line strings far apart so extents don't overlap - l1 = LG.LineString([[100.0, 0.0], [101.0, 0.0], [103.0, 0.0]]) - l2 = LG.LineString([[0.0, 0.0], [1.0, 0.0], [3.0, 0.0]]) - @test !GO.intersects(l1, l2) - @test isnothing(GO.intersection(l1, l2)) +mpt1 = LG.MultiPoint([pt1, pt2]) +mpt2 = LG.MultiPoint([pt2, pt3]) +mpt3 = LG.MultiPoint([pt4, pt5]) +ml1 = LG.MultiLineString([l5, l6, l7]) +ml2 = LG.MultiLineString([l1]) +mp1 = LG.MultiPolygon([p1]) +mp2 = LG.MultiPolygon([p6, p7]) +gc1 = LG.GeometryCollection([pt1, l5, p6]) - # Line strings close together that don't overlap - l1 = LG.LineString([[3.0, 0.25], [5.0, 0.25], [7.0, 0.25]]) - l2 = LG.LineString([[0.0, 0.0], [5.0, 10.0], [10.0, 0.0]]) - @test !GO.intersects(l1, l2) - @test isempty(GO.intersection(l1, l2)) - - # Closed linear ring with open line string - r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) - l2 = LG.LineString([[0.0, -2.0], [12.0, 10.0],]) - @test GO.intersects(r1, l2) - go_inter = GO.intersection(r1, l2) - @test length(go_inter) == 2 - lg_inter = LG.intersection(r1, l2) - @test issetequal( - Set(go_inter), - Set(GO._tuple_point.(GI.getpoint(lg_inter))) - ) +test_pairs = [ + # Points and geometries + (pt1, pt1, "pt1", "pt1"), # Same point + (pt1, pt2, "pt1", "pt2"), # Different point + (pt1, l1, "pt1", "l1"), # Point on line endpoint + (pt2, l1, "pt2", "l1"), # Point outside line + (pt3, l1, "pt3", "l1"), # Point on line segment + (pt4, l1, "pt4", "l1"), # Point on line vertex between segments + (l1, pt3, "l1", "pt3"), # oint on line segment (order swapped) + (pt1, r1, "pt1", "r1"), # Point on ring "endpoint" + (pt2, r1, "pt2", "r1"), # Point outside ring + (pt3, r1, "pt3", "r1"), # Point on ring segment + (pt4, r1, "pt4", "r1"), # Point on ring vertex between segments + (r1, pt3, "r1", "pt3"), # Point on ring segment (order swapped) + (p1, pt1, "p1", "pt1"), # Point on vertex of polygon + (pt2, p1, "pt2", "p1"), # Point outside of polygon's external ring + (pt4, p1, "pt4", "p1"), # Point on polygon's edge + (pt5, p1, "pt5", "p1"), # Point inside of polygon + (pt6, p1, "pt6", "p1"), # Point on hole edge + (pt7, p1, "pt7", "p1"), # Point inside of polygon hole + (p1, pt5, "p1", "pt5"), # Point inside of polygon (order swapped) + # # Lines and geometries + (l1, l1, "l1", "l1"), # Same line + (l2, l1, "l2", "l1"), # L2 is one segment of l1 + (l3, l1, "l3", "l1"), # L3 shares one segment with l1 and has one segment outside + (l4, l1, "l4", "l1"), # L4 shares half of each of l1 segments + (l5, l1, "l5", "l1"), # L5 shares one endpoint with l1 but not segments + (l6, l1, "l6", "l1"), # Lines are disjoint + (l7, l1, "l7", "l1"), # L7 crosses through one of l1's segments + (l8, l1, "l8", "l1"), # Overlaps one segment and crosses through another segment + (l9, l1, "l9", "l1"), # Two segments touch and two segments crosses + (l16, l1, "l16", "l1"), # L16 bounces off of l1's corner + (l1, r1, "l1", "r1"), # Line inside of ring + (l3, r1, "l3", "r1"), # Line covers one edge of linear ring and has segment outside + (l5, r1, "l5", "r1"), # Line and linear ring are only covered by vertex + (l6, r1, "l6", "r1"), # Line and ring are disjoint + (l7, r1, "l7", "r1"), # Line crosses through two ring edges + (l8, r1, "l8", "r1"), # Line crosses through two ring edges and touches third edge + (l10, r1, "l10", "r1"), # Line is equal to linear ring + (l11, r1, "l11", "r1"), # Line covers linear ring and then has extra segment + (l1, p1, "l1", "p1"), # Line on polygon edge + (l3, p1, "l3", "p1"), # Line on polygon edge and extending beyond polygon edge + (l5, p1, "l5", "p1"), # Line outside polygon connected by a vertex + (l7, p1, "l7", "p1"), # Line through polygon cutting to the outside + (l12, p1, "l12", "p1"), # Line inside polygon + (l13, p1, "l13", "p1"), # Line outside of polygon + (l14, p1, "l14", "p1"), # Line in polygon hole + (l15, p8, "l15", "p8"), # Line outside crown-shaped polygon but touching edges + (l15, p9, "l15", "p9"), # Line within crown-shaped polygon but touching edges + # Ring and geometries + (r1, l1, "r1", "l1"), # Line is within linear ring + (r1, l3, "r1", "l3"), # Line covers one edge of linear ring and has segment outside + (r1, l5, "r1", "l5"), # Line and linear ring are only connected at vertex + (r1, l6, "r1", "l6"), # Line and linear ring are disjoint + (r1, l7, "r1", "l7"), # Line crosses though two ring edges + (r1, l8, "r1", "l8"), # Line crosses through two ring edges and touches third edge + (r1, l10, "r1", "l10"), # Line is equal to linear ring + (r1, l11, "r1", "l11"), # Line covers linear ring and then has extra segment + (r1, r1, "r1", "r1"), # Same rings + (r2, r1, "r2", "r1"), # Disjoint ring with one "inside" of hole created + (r3, r1, "r3", "r1"), # Disjoint ring with one "outside" of hole created + (r4, r1, "r4", "r1"), # Rings share two sides and rest of sides dont touch + (r1, r5, "r1", "r5"), # Ring shares all edges with other ring, plus an extra loop + (r1, r6, "r1", "r6"), # Rings share just one vertex + (r1, r7, "r1", "r7"), # Rings cross one another + (r4, p1, "r4", "p1"), # Ring on boundary of polygon + (r1, p1, "r1", "p1"), # Ring on boundary and cutting through polygon + (r2, p1, "r2", "p1"), # Ring on hole bounday + (r6, p1, "r6", "p1"), # Ring touches polygon at one vertex + (r7, p1, "r7", "p1"), # Ring crosses through polygon + (r8, p1, "r8", "p1"), # Ring inside of polygon + (r9, p1, "r9", "p1"), # Ring outside of polygon + (r10, p1, "r10", "p1"), # Ring inside of polygon and shares hole's edge + (r11, p1, "r11", "p1"), # Ring inside of polygon hole + # Polygon and geometries + (p1, p1, "p1", "p1"), # Same polygons + (p1, p2, "p1", "p2"), # P1 and p2 are the same but p1 has holes + (p2, p1, "p2", "p1"), # P1 and p2 are the same but p1 has holes (order swapped) + (p3, p1, "p3", "p1"), # P3 is equal to one of p1's holes + (p4, p1, "p4", "p1"), # Polygon's share just one vertex + (p5, p1, "p5", "p1"), # Polygon outside of other polygon + (p6, p1, "p6", "p1"), # Polygon inside of other polygon's hole + (p7, p1, "p7", "p1"), # Polygons overlap + (p10, p1, "p10", "p1"), # Polygon's with nested holes + # Multigeometries + (mpt1, mpt1, "mpt1", "mpt1"), # Same set of points for multipoints + (mpt1, mpt2, "mpt1", "mpt2"), # Some point matches, others are different + (mpt1, mpt3, "mpt1", "mpt3"), # No shared points + (ml1, ml2, "ml1", "ml2"), # Lines in ml1 cross and touch ml2 + (mp1, mp2, "mp1", "mp2"), # Polygons in mp1 are inside hole and overlap + (gc1, ml1, "gc1", "ml1"), # Make sure collection works with multi-geom +] - # Closed linear ring with closed linear ring - r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]) - r2 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]) - @test GO.intersects(r1, r2) - go_inter = GO.intersection(r1, r2) - @test length(go_inter) == 2 - lg_inter = LG.intersection(r1, r2) - @test issetequal( - Set(go_inter), - Set(GO._tuple_point.(GI.getpoint(lg_inter))) - ) - end - - @testset "Polygons" begin - # Two polygons that intersect - p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) - p2 = LG.Polygon([[[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]]) - @test GO.intersects(p1, p2) - @test all(GO.intersection_points(p1, p2) .== [(6.5, 3.5), (6.5, -3.5)]) - - # Two polygons that don't intersect - p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) - p2 = LG.Polygon([[[13.0, 0.0], [18.0, 5.0], [23.0, 0.0], [18.0, -5.0], [13.0, 0.0]]]) - @test !GO.intersects(p1, p2) - @test isnothing(GO.intersection_points(p1, p2)) - - # Polygon that intersects with linestring - p1 = LG.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]]) - l2 = LG.LineString([[0.0, 0.0], [10.0, 0.0]]) - @test GO.intersects(p1, l2) - GO.intersection_points(p1, l2) - @test all(GO.intersection_points(p1, l2) .== [(0.0, 0.0), (10.0, 0.0)]) - - # Polygon with a hole, line through polygon and hole - p1 = LG.Polygon([ - [[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]], - [[2.0, -1.0], [2.0, 1.0], [3.0, 1.0], [3.0, -1.0], [2.0, -1.0]] - ]) - l2 = LG.LineString([[0.0, 0.0], [10.0, 0.0]]) - @test GO.intersects(p1, l2) - @test all(GO.intersection_points(p1, l2) .== [(0.0, 0.0), (2.0, 0.0), (3.0, 0.0), (10.0, 0.0)]) - - # Polygon with a hole, line only within the hole - p1 = LG.Polygon([ - [[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]], - [[2.0, -1.0], [2.0, 1.0], [3.0, 1.0], [3.0, -1.0], [2.0, -1.0]] - ]) - l2 = LG.LineString([[2.25, 0.0], [2.75, 0.0]]) - @test !GO.intersects(p1, l2) - @test isempty(GO.intersection_points(p1, l2)) +function test_geom_relation(GO_f, LG_f, f_name; swap_points = false) + for (g1, g2, sg1, sg2) in test_pairs + if swap_points + g1, g2 = g2, g1 + sg1, sg2 = sg2, sg1 + end + go_val = GO_f(g1, g2) + lg_val = LG_f(g1, g2) + @test go_val == lg_val + go_val != lg_val && println("\n" * "↑ TEST INFO: " * sg1 * " " * f_name * " " * sg2 * "\n\n") end +end - @testset "MultiPolygons" begin - # TODO: Add these tests - # Multi-polygon and polygon that intersect +@testset "Contains" begin test_geom_relation(GO.contains, LG.contains, "contains"; swap_points = true) end +@testset "Covered By" begin test_geom_relation(GO.coveredby, LG.coveredby, "coveredby") end +@testset "Covers" begin test_geom_relation(GO.covers, LG.covers, "covers"; swap_points = true) end +@testset "Disjoint" begin test_geom_relation(GO.disjoint, LG.disjoint, "disjoint")end +@testset "Intersect" begin test_geom_relation(GO.intersects, LG.intersects, "intersects") end +@testset "Touches" begin test_geom_relation(GO.touches, LG.touches, "touches") end +@testset "Within" begin test_geom_relation(GO.within, LG.within, "within") end - # Multi-polygon and polygon that don't intersect - # Multi-polygon that intersects with linestring - - end -end -@testset "overlaps" begin +@testset "Overlaps" begin @testset "Points/MultiPoints" begin p1 = LG.Point([0.0, 0.0]) p2 = LG.Point([0.0, 1.0]) @@ -251,64 +282,9 @@ end @test GO.overlaps(m1, p4) == LG.overlaps(m1, p4) end end -@testset "Contains, within, disjoint" begin - poly3 = GI.Polygon([[(1, 1), (1, 10), (10, 10), (10, 1), (1, 1)]]) - poly4 = GI.Polygon([[(1, 1), (2, 2), (3, 2), (1, 1)]]) - line5 = GI.LineString([(1.0, 1.0), (2.0, 3.0), (2.0, 3.5)]) - +@testset "Crosses" begin line6 = GI.LineString([(1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0)]) - poly5 = GI.Polygon([[(1.0, 1.0), (1.0, 20.0), (1.0, 3.0), (1.0, 4.0), (1.0, 1.0)]]) - line7 = GI.LineString([(1.0, 2.0), (1.0, 3.0), (1.0, 3.5)]) - - @test GO.contains(poly3, poly4) == true - @test GO.contains(poly3, line5) == true - @test GO.contains(line6, (1, 2)) == true - @test GO.contains(poly3, poly5) == false - @test GO.contains(poly3 , line7) == false - - @test GO.within(poly4, poly3) == true - @test GO.within(line5, poly3) == true - @test GO.within(poly5, poly3) == false - @test GO.within((1, 2), line6) == true - @test GO.within(line7, poly3) == false - - poly6 = GI.Polygon([[(-11, -12), (-13, -12), (-13, -13), (-11, -13), (-11, -12)]]) poly7 = GI.Polygon([[(-1, 2), (3, 2), (3, 3), (-1, 3), (-1, 2)]]) - poly8 = GI.Polygon([[(-1, 2), (-13, -12), (-13, -13), (-11, -13), (-1, 2)]]) - - @test GO.disjoint(poly7, poly6) == true - @test GO.disjoint(poly7, (1, 1)) == true - @test GO.disjoint(poly7, GI.LineString([(0, 0), (12, 2), (12, 3), (12, 4)])) == true - @test GO.disjoint(poly8, poly7) == false - - line8 = GI.LineString([(124.584961, -12.768946), (126.738281, -17.224758)]) - line9 = GI.LineString([(123.354492, -15.961329), (127.22168, -14.008696)]) - - @test all(GO.intersection(line8, line9)[1] .≈ (125.583754, -14.835723)) - - line10 = GI.LineString([ - (142.03125, -11.695273), - (138.691406, -16.804541), - (136.40625, -14.604847), - (135.966797, -12.039321), - (131.308594, -11.436955), - (128.232422, -15.36895), - (125.947266, -13.581921), - (121.816406, -18.729502), - (117.421875, -20.632784), - (113.378906, -23.402765), - (114.169922, -26.667096), - ]) - line11 = GI.LineString([ - (117.861328, -15.029686), - (122.124023, -24.886436), - (132.583008, -22.309426), - (132.890625, -7.754537), - ]) - - points = GO.intersection(line10, line11) - @test all(points[1] .≈ (119.832884, -19.58857)) - @test all(points[2] .≈ (132.808697, -11.6309378)) @test GO.crosses(GI.LineString([(-2, 2), (4, 2)]), line6) == true @test GO.crosses(GI.LineString([(0.5, 2.5), (1.0, 1.0)]), poly7) == true From 1164784f2538a663fe3232a9b6593c549ad83d37 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 12:08:33 -0800 Subject: [PATCH 10/18] Add _isparallel --- src/methods/bools.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/methods/bools.jl b/src/methods/bools.jl index 30b8716e1..0f32d1c0e 100644 --- a/src/methods/bools.jl +++ b/src/methods/bools.jl @@ -128,7 +128,10 @@ end # return slope1 === slope2 # end +_isparallel(((ax, ay), (bx, by)), ((cx, cy), (dx, dy))) = + _isparallel(bx - ax, by - ay, dx - cx, dy - cy) +_isparallel(Δx1, Δy1, Δx2, Δy2) = (Δx1 * Δy2 == Δy1 * Δx2) """ point_on_line(point::Point, line::LineString; ignore_end_vertices::Bool=false)::Bool From 87268273db08933a2c96bd17480338dbf7f4e354 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 12:32:59 -0800 Subject: [PATCH 11/18] Fix overlaps and crosses tests --- src/methods/geom_relations/crosses.jl | 2 +- .../geom_relations/geom_geom_processors.jl | 8 ++--- src/methods/geom_relations/overlaps.jl | 34 +++++++++++++------ test/methods/geom_relations.jl | 8 ++--- 4 files changed, 33 insertions(+), 19 deletions(-) diff --git a/src/methods/geom_relations/crosses.jl b/src/methods/geom_relations/crosses.jl index f8a580db0..00a5ecba5 100644 --- a/src/methods/geom_relations/crosses.jl +++ b/src/methods/geom_relations/crosses.jl @@ -55,7 +55,7 @@ end function line_crosses_line(line1, line2) np2 = GI.npoint(line2) - if intersects(line1, line2) + if GeometryOps.intersects(line1, line2) for i in 1:GI.npoint(line1) - 1 for j in 1:GI.npoint(line2) - 1 exclude_boundary = (j === 1 || j === np2 - 2) ? :none : :both diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl index e1fee57db..f5b951f57 100644 --- a/src/methods/geom_relations/geom_geom_processors.jl +++ b/src/methods/geom_relations/geom_geom_processors.jl @@ -545,10 +545,10 @@ function _segment_segment_orientation( cross::T = line_cross, hinge::T = line_hinge, over::T = line_over, out::T = line_out, ) where T - (ax, ay) = _tuple_point(a_point) - (bx, by) = _tuple_point(b_point) - (cx, cy) = _tuple_point(c_point) - (dx, dy) = _tuple_point(d_point) + (ax, ay) = Float64.(_tuple_point(a_point)) + (bx, by) = Float64.(_tuple_point(b_point)) + (cx, cy) = Float64.(_tuple_point(c_point)) + (dx, dy) = Float64.(_tuple_point(d_point)) meet_type = ExactPredicates.meet((ax, ay), (bx, by), (cx, cy), (dx, dy)) # Lines meet at one point within open segments meet_type == 1 && return cross diff --git a/src/methods/geom_relations/overlaps.jl b/src/methods/geom_relations/overlaps.jl index f99b75c9d..3cbaa380f 100644 --- a/src/methods/geom_relations/overlaps.jl +++ b/src/methods/geom_relations/overlaps.jl @@ -157,7 +157,7 @@ function overlaps( trait_b::GI.PolygonTrait, poly_b, ) edges_a, edges_b = map(sort! ∘ to_edges, (poly_a, poly_b)) - return _line_intersects(edges_a, edges_b) && + return _edge_intersects(edges_a, edges_b) && !equals(trait_a, poly_a, trait_b, poly_b) end @@ -211,15 +211,8 @@ function overlaps( return false end -""" - _overlaps( - (a1, a2)::Edge, - (b1, b2)::Edge - )::Bool - -If the edges overlap, meaning that they are colinear but each have one endpoint -outside of the other edge, return true. Else false. -""" +#= If the edges overlap, meaning that they are colinear but each have one endpoint +outside of the other edge, return true. Else false. =# function _overlaps( (a1, a2)::Edge, (b1, b2)::Edge @@ -231,3 +224,24 @@ function _overlaps( b_fully_within = point_on_seg(b1, a1, a2) && point_on_seg(b2, a1, a2) return on_top && (!a_fully_within && !b_fully_within) end + +# Checks it vectors of edges intersect +function _edge_intersects( + edges_a::Vector{Edge}, + edges_b::Vector{Edge} +) + # Extents.intersects(to_extent(edges_a), to_extent(edges_b)) || return false + for edge_a in edges_a + for edge_b in edges_b + _edge_intersects(edge_a, edge_b) && return true + end + end + return false +end + + +# Checks if two edges intersect +function _edge_intersects(edge_a::Edge, edge_b::Edge) + meet_type = ExactPredicates.meet(edge_a..., edge_b...) + return meet_type == 0 || meet_type == 1 +end diff --git a/test/methods/geom_relations.jl b/test/methods/geom_relations.jl index 815cd128d..44b022670 100644 --- a/test/methods/geom_relations.jl +++ b/test/methods/geom_relations.jl @@ -284,11 +284,11 @@ end end @testset "Crosses" begin line6 = GI.LineString([(1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0)]) - poly7 = GI.Polygon([[(-1, 2), (3, 2), (3, 3), (-1, 3), (-1, 2)]]) + poly7 = GI.Polygon([[(-1.0, 2.0), (3.0, 2.0), (3.0, 3.0), (-1.0, 3.0), (-1.0, 2.0)]]) @test GO.crosses(GI.LineString([(-2, 2), (4, 2)]), line6) == true @test GO.crosses(GI.LineString([(0.5, 2.5), (1.0, 1.0)]), poly7) == true - @test GO.crosses(GI.MultiPoint([(1, 2), (12, 12)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == true - @test GO.crosses(GI.MultiPoint([(1, 0), (12, 12)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == false - @test GO.crosses(GI.LineString([(-2, 2), (-4, 2)]), poly7) == false + @test GO.crosses(GI.MultiPoint([(1.0, 2.0), (12.0, 12.0)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == true + @test GO.crosses(GI.MultiPoint([(1.0, 0.0), (12.0, 12.0)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == false + @test GO.crosses(GI.LineString([(-2.0, 2.0), (-4.0, 2.0)]), poly7) == false end \ No newline at end of file From 6d2b6f3c64827788e2e5dfe6b7babc4be0e1d978 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 12:53:55 -0800 Subject: [PATCH 12/18] clean up bools --- src/methods/bools.jl | 248 ------------------------- src/methods/equals.jl | 4 +- src/methods/geom_relations/crosses.jl | 51 ++++- src/methods/geom_relations/overlaps.jl | 33 +++- src/methods/geom_relations/within.jl | 4 +- test/methods/bools.jl | 6 +- 6 files changed, 81 insertions(+), 265 deletions(-) diff --git a/src/methods/bools.jl b/src/methods/bools.jl index 0f32d1c0e..b1cb9e878 100644 --- a/src/methods/bools.jl +++ b/src/methods/bools.jl @@ -1,8 +1,6 @@ # # Boolean conditions export isclockwise, isconcave -export point_on_line, point_in_polygon, point_in_ring -export line_on_line, line_in_polygon, polygon_in_polygon # These are all adapted from Turf.jl. @@ -133,256 +131,10 @@ _isparallel(((ax, ay), (bx, by)), ((cx, cy), (dx, dy))) = _isparallel(Δx1, Δy1, Δx2, Δy2) = (Δx1 * Δy2 == Δy1 * Δx2) -""" - point_on_line(point::Point, line::LineString; ignore_end_vertices::Bool=false)::Bool - -Return true if a point is on a line. Accept a optional parameter to ignore the -start and end vertices of the linestring. - -## Examples - -```jldoctest -import GeoInterface as GI, GeometryOps as GO - -point = (1, 1) -line = GI.LineString([(0, 0), (3, 3), (4, 4)]) -GO.point_on_line(point, line) - -# output -true -``` -""" -function point_on_line(point, line; ignore_end_vertices::Bool=false)::Bool - line_points = tuple_points(line) - n = length(line_points) - - exclude_boundary = :none - for i in 1:n - 1 - if ignore_end_vertices - if i === 1 - exclude_boundary = :start - elseif i === n - 2 - exclude_boundary = :end - elseif (i === 1 && i + 1 === n - 1) - exclude_boundary = :both - end - end - if point_on_segment(point, (line_points[i], line_points[i + 1]); exclude_boundary) - return true - end - end - return false -end - -function point_on_seg(point, start, stop) - # Parse out points - x, y = GI.x(point), GI.y(point) - x1, y1 = GI.x(start), GI.y(start) - x2, y2 = GI.x(stop), GI.y(stop) - Δxl = x2 - x1 - Δyl = y2 - y1 - # Determine if point is on segment - cross = (x - x1) * Δyl - (y - y1) * Δxl - if cross == 0 # point is on line extending to infinity - # is line between endpoints - if abs(Δxl) >= abs(Δyl) # is line between endpoints - return Δxl > 0 ? x1 <= x <= x2 : x2 <= x <= x1 - else - return Δyl > 0 ? y1 <= y <= y2 : y2 <= y <= y1 - end - end - return false -end - -function point_on_segment(point, (start, stop); exclude_boundary::Symbol=:none)::Bool - x, y = GI.x(point), GI.y(point) - x1, y1 = GI.x(start), GI.y(start) - x2, y2 = GI.x(stop), GI.y(stop) - - dxc = x - x1 - dyc = y - y1 - dx1 = x2 - x1 - dy1 = y2 - y1 - - # TODO use better predicate for crossing here - cross = dxc * dy1 - dyc * dx1 - cross != 0 && return false - - # Will constprop optimise these away? - if exclude_boundary === :none - if abs(dx1) >= abs(dy1) - return dx1 > 0 ? x1 <= x && x <= x2 : x2 <= x && x <= x1 - end - return dy1 > 0 ? y1 <= y && y <= y2 : y2 <= y && y <= y1 - elseif exclude_boundary === :start - if abs(dx1) >= abs(dy1) - return dx1 > 0 ? x1 < x && x <= x2 : x2 <= x && x < x1 - end - return dy1 > 0 ? y1 < y && y <= y2 : y2 <= y && y < y1 - elseif exclude_boundary === :end - if abs(dx1) >= abs(dy1) - return dx1 > 0 ? x1 <= x && x < x2 : x2 < x && x <= x1 - end - return dy1 > 0 ? y1 <= y && y < y2 : y2 < y && y <= y1 - elseif exclude_boundary === :both - if abs(dx1) >= abs(dy1) - return dx1 > 0 ? x1 < x && x < x2 : x2 < x && x < x1 - end - return dy1 > 0 ? y1 < y && y < y2 : y2 < y && y < y1 - end - return false -end - -""" - point_in_polygon(point::Point, polygon::Union{Polygon, MultiPolygon}, ignore_boundary::Bool=false)::Bool - -Take a Point and a Polygon and determine if the point -resides inside the polygon. The polygon can be convex or concave. The function accounts for holes. - -## Examples - -```jldoctest -import GeoInterface as GI, GeometryOps as GO - -point = (-77.0, 44.0) -poly = GI.Polygon([[(-81, 41), (-81, 47), (-72, 47), (-72, 41), (-81, 41)]]) -GO.point_in_polygon(point, poly) - -# output -true -``` -""" -point_in_polygon(point, polygon; kw...)::Bool = - point_in_polygon(GI.trait(point), point, GI.trait(polygon), polygon; kw...) -function point_in_polygon( - ::PointTrait, point, - ::PolygonTrait, poly; - ignore_boundary::Bool=false, - check_extent::Bool=false, -)::Bool - # Cheaply check that the point is inside the polygon extent - if check_extent - point_in_extent(point, GI.extent(poly)) || return false - end - - # Then check the point is inside the exterior ring - point_in_polygon( - point,GI.getexterior(poly); - ignore_boundary, check_extent=false, - ) || return false - - # Finally make sure the point is not in any of the holes, - # flipping the boundary condition - for ring in GI.gethole(poly) - point_in_polygon( - point, ring; - ignore_boundary=!ignore_boundary, - ) && return false - end - return true -end - -function point_in_polygon( - ::PointTrait, pt, - ::Union{LineStringTrait,LinearRingTrait}, ring; - ignore_boundary::Bool=false, - check_extent::Bool=false, -)::Bool - x, y = GI.x(pt), GI.y(pt) - # Cheaply check that the point is inside the ring extent - if check_extent - point_in_extent(point, GI.extent(ring)) || return false - end - # Then check the point is inside the ring - inside = false - n = GI.npoint(ring) - p_start = GI.getpoint(ring, 1) - p_end = GI.getpoint(ring, n) - # Handle closed vs opne rings - if GI.x(p_start) == GI.x(p_end) && GI.y(p_start) == GI.y(p_end) - n -= 1 - end - # Loop over all points in the ring - for i in 1:(n - 1) - # First point on edge - p_i = GI.getpoint(ring, i) - xi, yi = GI.x(p_i), GI.y(p_i) - # Second point on edge (j = i + 1) - p_j = GI.getpoint(ring, i + 1) - xj, yj = GI.x(p_j), GI.y(p_j) - # Check if point is on the ring boundary - on_boundary = ( # vertex to point has same slope as edge - yi * (xj - x) + yj * (x - xi) == y * (xj - xi) && - (xi - x) * (xj - x) <= 0 && # x is between xi and xj - (yi - y) * (yj - y) <= 0 # y is between yi and yj - ) - on_boundary && return !ignore_boundary - # Check if ray from point passes through edge - intersects = ( - (yi > y) !== (yj > y) && - (x < (xj - xi) * (y - yi) / (yj - yi) + xi) - ) - if intersects - inside = !inside - end - end - return inside -end function point_in_extent(p, extent::Extents.Extent) (x1, x2), (y1, y1) = extent.X, extent.Y return x1 <= GI.x(p) && y1 <= GI.y(p) && x2 >= GI.x(p) && y2 >= GI.y(p) end -line_on_line(line1, line2) = line_on_line(trait(line1), line1, trait(line2), line2) -function line_on_line(t1::GI.AbstractCurveTrait, line1, t2::AbstractCurveTrait, line2) - for p in GI.getpoint(line1) - # FIXME: all points being on the line doesn't - # actually mean the whole line is on the line... - point_on_line(p, line2) || return false - end - return true -end - -line_in_polygon(line, poly) = line_in_polygon(trait(line), line, trait(poly), poly) - -function line_in_polygon( - ::AbstractCurveTrait, line, - ::Union{AbstractPolygonTrait,LinearRingTrait}, poly -) - Extents.intersects(GI.extent(poly), GI.extent(line)) || return false - - inside = false - for i in 1:GI.npoint(line) - 1 - p = GI.getpoint(line, i) - p2 = GI.getpoint(line, i + 1) - point_in_polygon(p, poly) || return false - if !inside - inside = point_in_polygon(p, poly; ignore_boundary=true) - end - # FIXME This seems like a hack, we should check for intersections rather than midpoint?? - if !inside - mid = ((GI.x(p) + GI.x(p2)) / 2, (GI.y(p) + GI.y(p2)) / 2) - inside = point_in_polygon(mid, poly; ignore_boundary=true) - end - end - return inside -end - -function polygon_in_polygon(poly1, poly2) - # edges1, edges2 = to_edges(poly1), to_edges(poly2) - # extent1, extent2 = to_extent(edges1), to_extent(edges2) - # Check the extents intersect - Extents.intersects(GI.extent(poly1), GI.extent(poly2)) || return false - - # Check all points in poly1 are in poly2 - for point in GI.getpoint(poly1) - point_in_polygon(point, poly2) || return false - end - - # Check the line of poly1 does not intersect the line of poly2 - #intersects(poly1, poly2) && return false - # poly1 must be in poly2 - return true - end diff --git a/src/methods/equals.jl b/src/methods/equals.jl index f4b710b81..c797432d5 100644 --- a/src/methods/equals.jl +++ b/src/methods/equals.jl @@ -9,7 +9,7 @@ The equals function checks if two geometries are equal. They are equal if they share the same set of points and edges to define the same shape. To provide an example, consider these two lines: -```@example cshape +```@example equals using GeometryOps using GeometryOps.GeometryBasics using Makie @@ -24,7 +24,7 @@ scatter!(GI.getpoint(l2), color = :orange) ``` We can see that the two lines do not share a commen set of points and edges in the plot, so they are not equal: -```@example cshape +```@example equals equals(l1, l2) # returns false ``` diff --git a/src/methods/geom_relations/crosses.jl b/src/methods/geom_relations/crosses.jl index 00a5ecba5..fc5f7a843 100644 --- a/src/methods/geom_relations/crosses.jl +++ b/src/methods/geom_relations/crosses.jl @@ -41,7 +41,7 @@ function multipoint_crosses_line(geom1, geom2) while i < GI.npoint(geom1) && !int_point && !ext_point for j in 1:GI.npoint(geom2) - 1 exclude_boundary = (j === 1 || j === np2 - 2) ? :none : :both - if point_on_segment(GI.getpoint(geom1, i), (GI.getpoint(geom2, j), GI.getpoint(geom2, j + 1)); exclude_boundary) + if _point_on_segment(GI.getpoint(geom1, i), (GI.getpoint(geom2, j), GI.getpoint(geom2, j + 1)); exclude_boundary) int_point = true else ext_point = true @@ -62,7 +62,7 @@ function line_crosses_line(line1, line2) pa = GI.getpoint(line1, i) pb = GI.getpoint(line1, i + 1) p = GI.getpoint(line2, j) - point_on_segment(p, (pa, pb); exclude_boundary) && return true + _point_on_segment(p, (pa, pb); exclude_boundary) && return true end end end @@ -81,7 +81,10 @@ function multipoint_crosses_poly(mp, poly) ext_point = false for p in GI.getpoint(mp) - if point_in_polygon(p, poly) + if _point_polygon_process( + p, poly; + in_allow = true, on_allow = true, out_allow = false, + ) int_point = true else ext_point = true @@ -90,3 +93,45 @@ function multipoint_crosses_poly(mp, poly) end return false end + +#= TODO: Once crosses is swapped over to use the geom relations workflow, can +delete these helpers. =# + +function _point_on_segment(point, (start, stop); exclude_boundary::Symbol=:none)::Bool + x, y = GI.x(point), GI.y(point) + x1, y1 = GI.x(start), GI.y(start) + x2, y2 = GI.x(stop), GI.y(stop) + + dxc = x - x1 + dyc = y - y1 + dx1 = x2 - x1 + dy1 = y2 - y1 + + # TODO use better predicate for crossing here + cross = dxc * dy1 - dyc * dx1 + cross != 0 && return false + + # Will constprop optimise these away? + if exclude_boundary === :none + if abs(dx1) >= abs(dy1) + return dx1 > 0 ? x1 <= x && x <= x2 : x2 <= x && x <= x1 + end + return dy1 > 0 ? y1 <= y && y <= y2 : y2 <= y && y <= y1 + elseif exclude_boundary === :start + if abs(dx1) >= abs(dy1) + return dx1 > 0 ? x1 < x && x <= x2 : x2 <= x && x < x1 + end + return dy1 > 0 ? y1 < y && y <= y2 : y2 <= y && y < y1 + elseif exclude_boundary === :end + if abs(dx1) >= abs(dy1) + return dx1 > 0 ? x1 <= x && x < x2 : x2 < x && x <= x1 + end + return dy1 > 0 ? y1 <= y && y < y2 : y2 < y && y <= y1 + elseif exclude_boundary === :both + if abs(dx1) >= abs(dy1) + return dx1 > 0 ? x1 < x && x < x2 : x2 < x && x < x1 + end + return dy1 > 0 ? y1 < y && y < y2 : y2 < y && y < y1 + end + return false +end diff --git a/src/methods/geom_relations/overlaps.jl b/src/methods/geom_relations/overlaps.jl index 3cbaa380f..32b9b6583 100644 --- a/src/methods/geom_relations/overlaps.jl +++ b/src/methods/geom_relations/overlaps.jl @@ -14,7 +14,7 @@ single point and a line only overlaps with another line if only a section of each line is colinear. To provide an example, consider these two lines: -```@example cshape +```@example overlaps using GeometryOps using GeometryOps.GeometryBasics using Makie @@ -28,7 +28,7 @@ lines!(GI.getpoint(l2), color = :orange) scatter!(GI.getpoint(l2), color = :orange) ``` We can see that the two lines overlap in the plot: -```@example cshape +```@example overlaps overlap(l1, l2) ``` @@ -220,11 +220,14 @@ function _overlaps( # meets in more than one point on_top = ExactPredicates.meet(a1, a2, b1, b2) == 0 # one end point is outside of other segment - a_fully_within = point_on_seg(a1, b1, b2) && point_on_seg(a2, b1, b2) - b_fully_within = point_on_seg(b1, a1, a2) && point_on_seg(b2, a1, a2) + a_fully_within = _point_on_seg(a1, b1, b2) && _point_on_seg(a2, b1, b2) + b_fully_within = _point_on_seg(b1, a1, a2) && _point_on_seg(b2, a1, a2) return on_top && (!a_fully_within && !b_fully_within) end +#= TODO: Once overlaps is swapped over to use the geom relations workflow, can +delete these helpers. =# + # Checks it vectors of edges intersect function _edge_intersects( edges_a::Vector{Edge}, @@ -239,9 +242,29 @@ function _edge_intersects( return false end - # Checks if two edges intersect function _edge_intersects(edge_a::Edge, edge_b::Edge) meet_type = ExactPredicates.meet(edge_a..., edge_b...) return meet_type == 0 || meet_type == 1 end + +# Checks if point is on a segment +function _point_on_seg(point, start, stop) + # Parse out points + x, y = GI.x(point), GI.y(point) + x1, y1 = GI.x(start), GI.y(start) + x2, y2 = GI.x(stop), GI.y(stop) + Δxl = x2 - x1 + Δyl = y2 - y1 + # Determine if point is on segment + cross = (x - x1) * Δyl - (y - y1) * Δxl + if cross == 0 # point is on line extending to infinity + # is line between endpoints + if abs(Δxl) >= abs(Δyl) # is line between endpoints + return Δxl > 0 ? x1 <= x <= x2 : x2 <= x <= x1 + else + return Δyl > 0 ? y1 <= y <= y2 : y2 <= y <= y1 + end + end + return false +end diff --git a/src/methods/geom_relations/within.jl b/src/methods/geom_relations/within.jl index 9b28c0ee7..cab30c80b 100644 --- a/src/methods/geom_relations/within.jl +++ b/src/methods/geom_relations/within.jl @@ -10,7 +10,7 @@ requires that the two interiors intersect and that the interior and boundary of the first geometry is not in the exterior of the second geometry. To provide an example, consider these two lines: -```@example cshape +```@example within using GeometryOps using GeometryOps.GeometryBasics using Makie @@ -25,7 +25,7 @@ scatter!(GI.getpoint(l2), color = :orange) ``` We can see that all of the points and edges of l2 are within l1, so l2 is within l1, but l1 is not within l2 -```@example cshape +```@example within within(l1, l2) # returns false within(l2, l1) # returns true ``` diff --git a/test/methods/bools.jl b/test/methods/bools.jl index 8f621ff07..35706e2bf 100644 --- a/test/methods/bools.jl +++ b/test/methods/bools.jl @@ -42,12 +42,8 @@ import GeometryOps as GO p3 = GI.Point([20, 20]) l5 = GI.LineString([[0, 0], [3, 3], [38.32, 5.96]]) - @test GO.point_on_line(p2, l4; ignore_end_vertices=true) == false - @test GO.point_on_line(p3, l5; ignore_end_vertices=true) == false - @test GO.point_on_line(p1, l3) == true - pt = (-77, 44) poly = GI.Polygon([[[-81, 41], [-81, 47], [-72, 47], [-72, 41], [-81, 41]]]) - @test point_in_polygon(pt, poly) == true + end From 790bd6a33429b0b03c8b049d2f88acf078713d48 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 13:13:41 -0800 Subject: [PATCH 13/18] Add more info to test printouts --- test/methods/geom_relations.jl | 174 ++++++++++++++++----------------- 1 file changed, 87 insertions(+), 87 deletions(-) diff --git a/test/methods/geom_relations.jl b/test/methods/geom_relations.jl index 44b022670..862cb5cb6 100644 --- a/test/methods/geom_relations.jl +++ b/test/methods/geom_relations.jl @@ -63,99 +63,99 @@ gc1 = LG.GeometryCollection([pt1, l5, p6]) test_pairs = [ # Points and geometries - (pt1, pt1, "pt1", "pt1"), # Same point - (pt1, pt2, "pt1", "pt2"), # Different point - (pt1, l1, "pt1", "l1"), # Point on line endpoint - (pt2, l1, "pt2", "l1"), # Point outside line - (pt3, l1, "pt3", "l1"), # Point on line segment - (pt4, l1, "pt4", "l1"), # Point on line vertex between segments - (l1, pt3, "l1", "pt3"), # oint on line segment (order swapped) - (pt1, r1, "pt1", "r1"), # Point on ring "endpoint" - (pt2, r1, "pt2", "r1"), # Point outside ring - (pt3, r1, "pt3", "r1"), # Point on ring segment - (pt4, r1, "pt4", "r1"), # Point on ring vertex between segments - (r1, pt3, "r1", "pt3"), # Point on ring segment (order swapped) - (p1, pt1, "p1", "pt1"), # Point on vertex of polygon - (pt2, p1, "pt2", "p1"), # Point outside of polygon's external ring - (pt4, p1, "pt4", "p1"), # Point on polygon's edge - (pt5, p1, "pt5", "p1"), # Point inside of polygon - (pt6, p1, "pt6", "p1"), # Point on hole edge - (pt7, p1, "pt7", "p1"), # Point inside of polygon hole - (p1, pt5, "p1", "pt5"), # Point inside of polygon (order swapped) + (pt1, pt1, "pt1", "pt1", "Same point"), + (pt1, pt2, "pt1", "pt2", "Different point"), + (pt1, l1, "pt1", "l1", "Point on line endpoint"), + (pt2, l1, "pt2", "l1", "Point outside line"), + (pt3, l1, "pt3", "l1", "Point on line segment"), + (pt4, l1, "pt4", "l1", "Point on line vertex between segments"), + (l1, pt3, "l1", "pt3", "Point on line segment (order swapped)"), + (pt1, r1, "pt1", "r1", "Point on ring 'endpoint'"), + (pt2, r1, "pt2", "r1", "Point outside ring"), + (pt3, r1, "pt3", "r1", "Point on ring segment"), + (pt4, r1, "pt4", "r1", "Point on ring vertex between segments"), + (r1, pt3, "r1", "pt3", "Point on ring segment (order swapped)"), + (p1, pt1, "p1", "pt1", "Point on vertex of polygon"), + (pt2, p1, "pt2", "p1", "Point outside of polygon's external ring"), + (pt4, p1, "pt4", "p1", "Point on polygon's edge"), + (pt5, p1, "pt5", "p1", "Point inside of polygon"), + (pt6, p1, "pt6", "p1", "Point on hole edge"), + (pt7, p1, "pt7", "p1", "Point inside of polygon hole"), + (p1, pt5, "p1", "pt5", "Point inside of polygon (order swapped)"), # # Lines and geometries - (l1, l1, "l1", "l1"), # Same line - (l2, l1, "l2", "l1"), # L2 is one segment of l1 - (l3, l1, "l3", "l1"), # L3 shares one segment with l1 and has one segment outside - (l4, l1, "l4", "l1"), # L4 shares half of each of l1 segments - (l5, l1, "l5", "l1"), # L5 shares one endpoint with l1 but not segments - (l6, l1, "l6", "l1"), # Lines are disjoint - (l7, l1, "l7", "l1"), # L7 crosses through one of l1's segments - (l8, l1, "l8", "l1"), # Overlaps one segment and crosses through another segment - (l9, l1, "l9", "l1"), # Two segments touch and two segments crosses - (l16, l1, "l16", "l1"), # L16 bounces off of l1's corner - (l1, r1, "l1", "r1"), # Line inside of ring - (l3, r1, "l3", "r1"), # Line covers one edge of linear ring and has segment outside - (l5, r1, "l5", "r1"), # Line and linear ring are only covered by vertex - (l6, r1, "l6", "r1"), # Line and ring are disjoint - (l7, r1, "l7", "r1"), # Line crosses through two ring edges - (l8, r1, "l8", "r1"), # Line crosses through two ring edges and touches third edge - (l10, r1, "l10", "r1"), # Line is equal to linear ring - (l11, r1, "l11", "r1"), # Line covers linear ring and then has extra segment - (l1, p1, "l1", "p1"), # Line on polygon edge - (l3, p1, "l3", "p1"), # Line on polygon edge and extending beyond polygon edge - (l5, p1, "l5", "p1"), # Line outside polygon connected by a vertex - (l7, p1, "l7", "p1"), # Line through polygon cutting to the outside - (l12, p1, "l12", "p1"), # Line inside polygon - (l13, p1, "l13", "p1"), # Line outside of polygon - (l14, p1, "l14", "p1"), # Line in polygon hole - (l15, p8, "l15", "p8"), # Line outside crown-shaped polygon but touching edges - (l15, p9, "l15", "p9"), # Line within crown-shaped polygon but touching edges + (l1, l1, "l1", "l1", "Same line"), + (l2, l1, "l2", "l1", "L2 is one segment of l1"), + (l3, l1, "l3", "l1", "L3 shares one segment with l1 and has one segment outside"), + (l4, l1, "l4", "l1", "L4 shares half of each of l1 segments"), + (l5, l1, "l5", "l1", "L5 shares one endpoint with l1 but not segments"), + (l6, l1, "l6", "l1", "Lines are disjoint"), + (l7, l1, "l7", "l1", "L7 crosses through one of l1's segments"), + (l8, l1, "l8", "l1", "Overlaps one segment and crosses through another segment"), + (l9, l1, "l9", "l1", "Two segments touch and two segments crosses"), + (l16, l1, "l16", "l1", "L16 bounces off of l1's corner"), + (l1, r1, "l1", "r1", "Line inside of ring"), + (l3, r1, "l3", "r1", "Line covers one edge of linear ring and has segment outside"), + (l5, r1, "l5", "r1", "Line and linear ring are only covered by vertex"), + (l6, r1, "l6", "r1", "Line and ring are disjoint"), + (l7, r1, "l7", "r1", "Line crosses through two ring edges"), + (l8, r1, "l8", "r1", "Line crosses through two ring edges and touches third edge"), + (l10, r1, "l10", "r1", "Line is equal to linear ring"), + (l11, r1, "l11", "r1", "Line covers linear ring and then has extra segment"), + (l1, p1, "l1", "p1", "Line on polygon edge"), + (l3, p1, "l3", "p1", "Line on polygon edge and extending beyond polygon edge"), + (l5, p1, "l5", "p1", "Line outside polygon connected by a vertex"), + (l7, p1, "l7", "p1", "Line through polygon cutting to the outside"), + (l12, p1, "l12", "p1", "Line inside polygon"), + (l13, p1, "l13", "p1", "Line outside of polygon"), + (l14, p1, "l14", "p1", "Line in polygon hole"), + (l15, p8, "l15", "p8", "Line outside crown-shaped polygon but touching edges"), + (l15, p9, "l15", "p9", "Line within crown-shaped polygon but touching edges"), # Ring and geometries - (r1, l1, "r1", "l1"), # Line is within linear ring - (r1, l3, "r1", "l3"), # Line covers one edge of linear ring and has segment outside - (r1, l5, "r1", "l5"), # Line and linear ring are only connected at vertex - (r1, l6, "r1", "l6"), # Line and linear ring are disjoint - (r1, l7, "r1", "l7"), # Line crosses though two ring edges - (r1, l8, "r1", "l8"), # Line crosses through two ring edges and touches third edge - (r1, l10, "r1", "l10"), # Line is equal to linear ring - (r1, l11, "r1", "l11"), # Line covers linear ring and then has extra segment - (r1, r1, "r1", "r1"), # Same rings - (r2, r1, "r2", "r1"), # Disjoint ring with one "inside" of hole created - (r3, r1, "r3", "r1"), # Disjoint ring with one "outside" of hole created - (r4, r1, "r4", "r1"), # Rings share two sides and rest of sides dont touch - (r1, r5, "r1", "r5"), # Ring shares all edges with other ring, plus an extra loop - (r1, r6, "r1", "r6"), # Rings share just one vertex - (r1, r7, "r1", "r7"), # Rings cross one another - (r4, p1, "r4", "p1"), # Ring on boundary of polygon - (r1, p1, "r1", "p1"), # Ring on boundary and cutting through polygon - (r2, p1, "r2", "p1"), # Ring on hole bounday - (r6, p1, "r6", "p1"), # Ring touches polygon at one vertex - (r7, p1, "r7", "p1"), # Ring crosses through polygon - (r8, p1, "r8", "p1"), # Ring inside of polygon - (r9, p1, "r9", "p1"), # Ring outside of polygon - (r10, p1, "r10", "p1"), # Ring inside of polygon and shares hole's edge - (r11, p1, "r11", "p1"), # Ring inside of polygon hole + (r1, l1, "r1", "l1", "Line is within linear ring"), + (r1, l3, "r1", "l3", "Line covers one edge of linear ring and has segment outside"), + (r1, l5, "r1", "l5", "Line and linear ring are only connected at vertex"), + (r1, l6, "r1", "l6", "Line and linear ring are disjoint"), + (r1, l7, "r1", "l7", "Line crosses though two ring edges"), + (r1, l8, "r1", "l8", "Line crosses through two ring edges and touches third edge"), + (r1, l10, "r1", "l10", "Line is equal to linear ring"), + (r1, l11, "r1", "l11", "Line covers linear ring and then has extra segment"), + (r1, r1, "r1", "r1", "Same rings"), + (r2, r1, "r2", "r1", "Disjoint ring with one 'inside' of hole created"), + (r3, r1, "r3", "r1", "Disjoint ring with one 'outside' of hole created"), + (r4, r1, "r4", "r1", "Rings share two sides and rest of sides dont touch"), + (r1, r5, "r1", "r5", "Ring shares all edges with other ring, plus an extra loop"), + (r1, r6, "r1", "r6", "Rings share just one vertex"), + (r1, r7, "r1", "r7", "Rings cross one another"), + (r4, p1, "r4", "p1", "Ring on boundary of polygon"), + (r1, p1, "r1", "p1", "Ring on boundary and cutting through polygon"), + (r2, p1, "r2", "p1", "Ring on hole bounday"), + (r6, p1, "r6", "p1", "Ring touches polygon at one vertex"), + (r7, p1, "r7", "p1", "Ring crosses through polygon"), + (r8, p1, "r8", "p1", "Ring inside of polygon"), + (r9, p1, "r9", "p1", "Ring outside of polygon"), + (r10, p1, "r10", "p1", "Ring inside of polygon and shares hole's edge"), + (r11, p1, "r11", "p1", "Ring inside of polygon hole"), # Polygon and geometries - (p1, p1, "p1", "p1"), # Same polygons - (p1, p2, "p1", "p2"), # P1 and p2 are the same but p1 has holes - (p2, p1, "p2", "p1"), # P1 and p2 are the same but p1 has holes (order swapped) - (p3, p1, "p3", "p1"), # P3 is equal to one of p1's holes - (p4, p1, "p4", "p1"), # Polygon's share just one vertex - (p5, p1, "p5", "p1"), # Polygon outside of other polygon - (p6, p1, "p6", "p1"), # Polygon inside of other polygon's hole - (p7, p1, "p7", "p1"), # Polygons overlap - (p10, p1, "p10", "p1"), # Polygon's with nested holes + (p1, p1, "p1", "p1", "Same polygons"), + (p1, p2, "p1", "p2", "P1 and p2 are the same but p1 has holes"), + (p2, p1, "p2", "p1", "P1 and p2 are the same but p1 has holes (order swapped)"), + (p3, p1, "p3", "p1", "P3 is equal to one of p1's holes"), + (p4, p1, "p4", "p1", "Polygon's share just one vertex"), + (p5, p1, "p5", "p1", "Polygon outside of other polygon"), + (p6, p1, "p6", "p1", "Polygon inside of other polygon's hole"), + (p7, p1, "p7", "p1", "Polygons overlap"), + (p10, p1, "p10", "p1", "Polygon's with nested holes"), # Multigeometries - (mpt1, mpt1, "mpt1", "mpt1"), # Same set of points for multipoints - (mpt1, mpt2, "mpt1", "mpt2"), # Some point matches, others are different - (mpt1, mpt3, "mpt1", "mpt3"), # No shared points - (ml1, ml2, "ml1", "ml2"), # Lines in ml1 cross and touch ml2 - (mp1, mp2, "mp1", "mp2"), # Polygons in mp1 are inside hole and overlap - (gc1, ml1, "gc1", "ml1"), # Make sure collection works with multi-geom + (mpt1, mpt1, "mpt1", "mpt1", "Same set of points for multipoints"), + (mpt1, mpt2, "mpt1", "mpt2", "Some point matches, others are different"), + (mpt1, mpt3, "mpt1", "mpt3", "No shared points"), + (ml1, ml2, "ml1", "ml2", "Lines in ml1 cross and touch ml2"), + (mp1, mp2, "mp1", "mp2", "Polygons in mp1 are inside hole and overlap"), + (gc1, ml1, "gc1", "ml1", "Make sure collection works with multi-geom"), ] function test_geom_relation(GO_f, LG_f, f_name; swap_points = false) - for (g1, g2, sg1, sg2) in test_pairs + for (g1, g2, sg1, sg2, sdesc) in test_pairs if swap_points g1, g2 = g2, g1 sg1, sg2 = sg2, sg1 @@ -163,7 +163,7 @@ function test_geom_relation(GO_f, LG_f, f_name; swap_points = false) go_val = GO_f(g1, g2) lg_val = LG_f(g1, g2) @test go_val == lg_val - go_val != lg_val && println("\n" * "↑ TEST INFO: " * sg1 * " " * f_name * " " * sg2 * "\n\n") + go_val != lg_val && println("\n↑ TEST INFO: - $sg1 $f_name $sg2: $sdesc \n\n") end end From 431273cbdc9d8d9b347775b2bf455b6f990b6771 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 14:01:01 -0800 Subject: [PATCH 14/18] Update touches doc fail --- src/methods/geom_relations/touches.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/methods/geom_relations/touches.jl b/src/methods/geom_relations/touches.jl index 24bf4b2ea..1991e9012 100644 --- a/src/methods/geom_relations/touches.jl +++ b/src/methods/geom_relations/touches.jl @@ -69,7 +69,7 @@ import GeometryOps as GO, GeoInterface as GI l1 = GI.Line([(0.0, 0.0), (1.0, 0.0)]) l2 = GI.Line([(0.5, 1.0), (0.5, -1.0)]) -GO.crosses(l1, l2) +GO.touches(l1, l2) # output true ``` From 74985387ead5684764be95bba05e5d0c5dd563e2 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Tue, 2 Jan 2024 14:18:39 -0800 Subject: [PATCH 15/18] Actually fix touches docs --- src/methods/geom_relations/touches.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/methods/geom_relations/touches.jl b/src/methods/geom_relations/touches.jl index 1991e9012..48eac4dc8 100644 --- a/src/methods/geom_relations/touches.jl +++ b/src/methods/geom_relations/touches.jl @@ -67,7 +67,7 @@ boundary. import GeometryOps as GO, GeoInterface as GI l1 = GI.Line([(0.0, 0.0), (1.0, 0.0)]) -l2 = GI.Line([(0.5, 1.0), (0.5, -1.0)]) +l2 = GI.Line([(1.0, 1.0), (1.0, -1.0)]) GO.touches(l1, l2) # output From 13125da21c767fafedbe4658080345b3b2080224 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Wed, 3 Jan 2024 14:36:47 -0800 Subject: [PATCH 16/18] Add double geom feature function --- src/methods/geom_relations/coveredby.jl | 2 +- src/methods/geom_relations/disjoint.jl | 8 ++++---- src/methods/geom_relations/intersects.jl | 2 +- src/methods/geom_relations/touches.jl | 2 +- src/methods/geom_relations/within.jl | 1 + 5 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/methods/geom_relations/coveredby.jl b/src/methods/geom_relations/coveredby.jl index e6ee52b15..0f71317b9 100644 --- a/src/methods/geom_relations/coveredby.jl +++ b/src/methods/geom_relations/coveredby.jl @@ -82,7 +82,7 @@ coveredby(g1, g2) = _coveredby(trait(g1), g1, trait(g2), g2) # # Convert features to geometries _coveredby(::GI.FeatureTrait, g1, ::Any, g2) = coveredby(GI.geometry(g1), g2) _coveredby(::Any, g1, t2::GI.FeatureTrait, g2) = coveredby(g1, GI.geometry(g2)) - +_coveredby(::FeatureTrait, g1, ::FeatureTrait, g2) = coveredby(GI.geometry(g1), GI.geometry(g2)) # # Points coveredby geometries diff --git a/src/methods/geom_relations/disjoint.jl b/src/methods/geom_relations/disjoint.jl index f3e051c05..e65bc70fb 100644 --- a/src/methods/geom_relations/disjoint.jl +++ b/src/methods/geom_relations/disjoint.jl @@ -70,12 +70,12 @@ GO.disjoint(point, line) true ``` """ -disjoint(g1, g2)::Bool = _disjoint(trait(g1), g1, trait(g2), g2) +disjoint(g1, g2) = _disjoint(trait(g1), g1, trait(g2), g2) # # Convert features to geometries -disjoint(::FeatureTrait, g1, ::Any, g2)::Bool = disjoint(GI.geometry(g1), g2) -disjoint(::Any, g1, t2::FeatureTrait, g2)::Bool = disjoint(g1, geometry(g2)) - +_disjoint(::FeatureTrait, g1, ::Any, g2) = disjoint(GI.geometry(g1), g2) +_disjoint(::Any, g1, ::FeatureTrait, g2) = disjoint(g1, geometry(g2)) +_disjoint(::FeatureTrait, g1, ::FeatureTrait, g2) = disjoint(GI.geometry(g1), GI.geometry(g2)) # # Point disjoint geometries diff --git a/src/methods/geom_relations/intersects.jl b/src/methods/geom_relations/intersects.jl index 834745a33..524746181 100644 --- a/src/methods/geom_relations/intersects.jl +++ b/src/methods/geom_relations/intersects.jl @@ -201,7 +201,7 @@ line segments, line strings, linear rings, polygons, and multipolygons. If no intersection points were possible given geometry extents, return nothing. If none are found, return an empty list. """ -function intersection_points(::GI.AbstractGeometryTrait, a, ::GI.AbstractGeometryTrait, b) +function intersection_points(::GI.AbstractTrait, a, ::GI.AbstractTrait, b) # Check if the geometries extents even overlap Extents.intersects(GI.extent(a), GI.extent(b)) || return nothing # Create a list of edges from the two input geometries diff --git a/src/methods/geom_relations/touches.jl b/src/methods/geom_relations/touches.jl index 48eac4dc8..21f9259d7 100644 --- a/src/methods/geom_relations/touches.jl +++ b/src/methods/geom_relations/touches.jl @@ -79,7 +79,7 @@ touches(g1, g2)::Bool = _touches(trait(g1), g1, trait(g2), g2) # # Convert features to geometries _touches(::GI.FeatureTrait, g1, ::Any, g2) = touches(GI.geometry(g1), g2) _touches(::Any, g1, t2::GI.FeatureTrait, g2) = touches(g1, GI.geometry(g2)) - +_touches(::FeatureTrait, g1, ::FeatureTrait, g2) = touches(GI.geometry(g1), GI.geometry(g2)) # # Point touches geometries diff --git a/src/methods/geom_relations/within.jl b/src/methods/geom_relations/within.jl index cab30c80b..f2833a11b 100644 --- a/src/methods/geom_relations/within.jl +++ b/src/methods/geom_relations/within.jl @@ -83,6 +83,7 @@ within(g1, g2) = _within(trait(g1), g1, trait(g2), g2) # # Convert features to geometries _within(::GI.FeatureTrait, g1, ::Any, g2) = within(GI.geometry(g1), g2) _within(::Any, g1, t2::GI.FeatureTrait, g2) = within(g1, GI.geometry(g2)) +_within(::FeatureTrait, g1, ::FeatureTrait, g2) = within(GI.geometry(g1), GI.geometry(g2)) # # Points within geometries From 3fd9637e62152b8cb3f63246015eaa58ee5051b2 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Wed, 3 Jan 2024 14:37:08 -0800 Subject: [PATCH 17/18] Update test output text --- test/methods/geom_relations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/methods/geom_relations.jl b/test/methods/geom_relations.jl index 862cb5cb6..50b0dfffd 100644 --- a/test/methods/geom_relations.jl +++ b/test/methods/geom_relations.jl @@ -163,7 +163,7 @@ function test_geom_relation(GO_f, LG_f, f_name; swap_points = false) go_val = GO_f(g1, g2) lg_val = LG_f(g1, g2) @test go_val == lg_val - go_val != lg_val && println("\n↑ TEST INFO: - $sg1 $f_name $sg2: $sdesc \n\n") + go_val != lg_val && println("\n↑ TEST INFO: $sg1 $f_name $sg2 - $sdesc \n\n") end end From d0c7f6110b361639eb26f9e6ecc6829ae841f9d6 Mon Sep 17 00:00:00 2001 From: Skylar Gering Date: Wed, 3 Jan 2024 14:37:28 -0800 Subject: [PATCH 18/18] Split while loop into helper functions --- .../geom_relations/geom_geom_processors.jl | 141 +++++++++--------- 1 file changed, 68 insertions(+), 73 deletions(-) diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl index f5b951f57..d26818083 100644 --- a/src/methods/geom_relations/geom_geom_processors.jl +++ b/src/methods/geom_relations/geom_geom_processors.jl @@ -132,41 +132,17 @@ function _line_curve_process( for j in (closed_curve ? 1 : 2):nc c_end = GI.getpoint(curve, j) # Check if line and curve segments meet - seg_val = _segment_segment_orientation( - (l_start, l_end), - (c_start, c_end), - ) + seg_val = _segment_segment_orientation((l_start, l_end), (c_start, c_end)) # If segments are co-linear if seg_val == line_over !over_allow && return false # at least one point in, meets requirments in_req_met = true - point_val = _point_segment_orientation( - l_start, - c_start, c_end, - ) + point_val = _point_segment_orientation(l_start, c_start, c_end) # If entire segment isn't covered, consider remaining section if point_val != point_out - if _point_segment_orientation( - l_end, - c_start, c_end, - ) != point_out - l_start = l_end - i += 1 - break - elseif _point_segment_orientation( - c_start, - l_start, l_end, - ) != point_out && !equals(l_start, c_start) - l_start = c_start - break - elseif _point_segment_orientation( - c_end, - l_start, l_end, - ) != point_out && !equals(l_start, c_end) - l_start = c_end - break - end + i, l_start, break_off = _find_new_seg(i, l_start, l_end, c_start, c_end) + break_off && break end else if seg_val == line_cross @@ -174,30 +150,10 @@ function _line_curve_process( in_req_met = true elseif seg_val == line_hinge # could cross or overlap # Determine location of intersection point on each segment - _, fracs = _intersection_point( - (_tuple_point(l_start), _tuple_point(l_end)), - (_tuple_point(c_start), _tuple_point(c_end)) - ) - (α, β) = - if !isnothing(fracs) - fracs - else # line and curve segments are parallel - if equals(l_start, c_start) - (0, 0) - elseif equals(l_start, c_end) - (0, 1) - elseif equals(l_end, c_start) - (1, 0) - else # equals(l_end, c_end) - (1, 1) - end - end - - if ( # don't consider edges of curves as they can't cross - (β == 0 && !closed_curve && j == 2) || - (β == 1 && !closed_curve && j == nc) || - (α == 0 && !closed_line && i == 2) || - (α == 1 && !closed_line && i == nl) + (α, β) = _find_intersect_fracs(l_start, l_end, c_start, c_end) + if ( # Don't consider edges of curves as they can't cross + (!closed_line && ((α == 0 && i == 2) || (α == 1 && i == nl))) || + (!closed_curve && ((β == 0 && j == 2) || (β == 1 && j == nc))) ) !on_allow && return false on_req_met = true @@ -205,24 +161,11 @@ function _line_curve_process( in_req_met = true # If needed, determine if hinge actually crosses if (!cross_allow || !over_allow) && α != 0 && β != 0 - l, c = if β == 1 - if α == 1 - ( - (l_end, GI.getpoint(line, i + 1)), - (c_end, GI.getpoint(curve, j + 1)) - ) - else - ( - (l_start, l_end), - (c_end, GI.getpoint(curve, j + 1)) - ) - end - else # β ≠ 1 and α == 1 - ( - (l_end, GI.getpoint(line, i + 1)), - (c_start, c_end) - ) - end + # Find next pieces of hinge to see if line and curve cross + l, c = _find_hinge_next_segments( + α, β, l_start, l_end, c_start, c_end, + i, line, j, curve, + ) if _segment_segment_orientation(l, c) == line_hinge !cross_allow && return false else @@ -231,14 +174,14 @@ function _line_curve_process( end end end - # no overlap for a give segment + # no overlap for a give segment, some of segment must be out of curve if j == nc !out_allow && return false out_req_met = true end end - c_start = c_end - if j == nc + c_start = c_end # consider next segment of curve + if j == nc # move on to next line segment i += 1 l_start = l_end end @@ -247,6 +190,58 @@ function _line_curve_process( return in_req_met && on_req_met && out_req_met end +#= If entire segment (le to ls) isn't covered by segment (cs to ce), find remaining section +part of section outside of cs to ce. If completly covered, increase segment index i. =# +function _find_new_seg(i, ls, le, cs, ce) + break_off = true + if _point_segment_orientation(le, cs, ce) != point_out + ls = le + i += 1 + elseif !equals(ls, cs) && _point_segment_orientation(cs, ls, le) != point_out + ls = cs + elseif !equals(ls, ce) && _point_segment_orientation(ce, ls, le) != point_out + ls = ce + else + break_off = false + end + return i, ls, break_off +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). =# +function _find_intersect_fracs(ls, le, cs, ce) + _, fracs = _intersection_point( + (_tuple_point(ls), _tuple_point(le)), + (_tuple_point(cs), _tuple_point(ce)) + ) + (α, β) = if !isnothing(fracs) + fracs + else # line and curve segments are parallel + if equals(ls, cs) + (0, 0) + elseif equals(ls, ce) + (0, 1) + elseif equals(le, cs) + (1, 0) + else # equals(l_end, c_end) + (1, 1) + end + end + return α, β +end + +#= Find next set of segments needed to determine if given hinge segments cross or not.=# +_find_hinge_next_segments(α, β, ls, le, cs, ce, i, line, j, curve) = + if β == 1 + if α == 1 # hinge at endpoints, so next segment of both is needed + ((le, GI.getpoint(line, i + 1)), (ce, GI.getpoint(curve, j + 1))) + else # hinge at curve endpoint and line interior point, curve next segment needed + ((ls, le), (ce, GI.getpoint(curve, j + 1))) + end + else # hinge at curve interior point and line endpoint, line next segment needed + ((le, GI.getpoint(line, i + 1)), (cs, ce)) + end + #= Determines if a line meets the given checks with respect to a polygon.