Skip to content

Commit

Permalink
Sg/de9im functions (#42)
Browse files Browse the repository at this point in the history 
* Added generalized code

* Update contains

* Add coveredby

* Add covers

* Update disjoint

* Update intersects

* Add touches

* Update within

* Update test file

* Add _isparallel

* Fix overlaps and crosses tests

* clean up bools

* Add more info to test printouts

* Update touches doc fail

* Actually fix touches docs

* Add double geom feature function

* Update test output text

* Split while loop into helper functions
  • Loading branch information
@skygering
skygering authored Jan 3, 2024
1 parent 35a3f43 commit 4918927
Show file tree
Hide file tree
Showing 15 changed files with 2,135 additions and 607 deletions.
6 changes: 5 additions & 1 deletion src/GeometryOps.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,15 @@ 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/covers.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/touches.jl")
include("methods/geom_relations/within.jl")
include("methods/polygonize.jl")

Expand Down
251 changes: 3 additions & 248 deletions src/methods/bools.jl
View file
Original file line number Diff line number Diff line change
@@ -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.

Expand Down Expand Up @@ -128,258 +126,15 @@ 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
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
4 changes: 2 additions & 2 deletions src/methods/equals.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
```
Expand Down
53 changes: 46 additions & 7 deletions src/methods/geom_relations/contains.jl
View file
Original file line number Diff line number Diff line change
@@ -1,25 +1,64 @@
# # 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
```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)
Loading

0 comments on commit 4918927

Please sign in to comment.