Skip to content

Commit

Permalink
Improve bools (#16)
Browse files Browse the repository at this point in the history 
* start improving and testing bools

* add turf boolean tests

* fix bool methods

* more fixes

* tests passing

* cleanup

* fix doctests

* dont check poly2 intersections check that

* fix doctests

* fix doctest
  • Loading branch information
@rafaqz
rafaqz authored Sep 1, 2023
1 parent 8375cb6 commit 98bfd3b
Show file tree
Hide file tree
Showing 14 changed files with 566 additions and 316 deletions.
1 change: 1 addition & 0 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ authors = ["Anshul Singhvi <[email protected]> and contributors"]
version = "0.0.1-DEV"

[deps]
ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
5 changes: 5 additions & 0 deletions src/GeometryOps.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ using GeoInterface
using GeometryBasics
import Proj
using LinearAlgebra
import ExactPredicates

using GeoInterface.Extents: Extents

Expand All @@ -21,6 +22,10 @@ include("methods/signed_area.jl")
include("methods/centroid.jl")
include("methods/intersects.jl")
include("methods/contains.jl")
include("methods/crosses.jl")
include("methods/disjoint.jl")
include("methods/overlaps.jl")
include("methods/within.jl")
include("methods/polygonize.jl")
include("methods/barycentric.jl")

Expand Down
239 changes: 114 additions & 125 deletions src/methods/bools.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,11 +13,14 @@ export line_on_line, line_in_polygon, polygon_in_polygon
Take a ring and return true or false whether or not the ring is clockwise or counter-clockwise.
## Examples
## Example
```jldoctest
import GeoInterface as GI, GeometryOps as GO
line = GI.LineString([(0, 0), (1, 1), (1, 0), (0, 0)])
GO.isclockwise(line)
ring = GI.LinearRing([(0, 0), (1, 1), (1, 0), (0, 0)])
GO.isclockwise(ring)
# output
true
```
Expand All @@ -43,8 +46,10 @@ Take a polygon and return true or false as to whether it is concave or not.
## Examples
```jldoctest
import GeoInterface as GI, GeometryOps as GO
poly = GI.Polygon([[(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]])
GO.isconcave(poly)
# output
false
```
Expand All @@ -53,8 +58,9 @@ function isconcave(poly)::Bool
sign = false

exterior = GI.getexterior(poly)
GI.npoint(exterior) <= 4 && return false

# FIXME handle not closed polygons
GI.npoint(exterior) <= 4 && return false
n = GI.npoint(exterior) - 1

for i in 1:n
Expand Down Expand Up @@ -82,33 +88,32 @@ function isconcave(poly)::Bool
return false
end

equals(geo1, geo2) = _equals(trait(geo1), geo1, trait(geo2), geo2)

function equals(geo1, geo2)
GI.geomtrait(geo1) !== GI.geomtrait(geo2) && return false

GI.geomtrait(geo1) isa PointTrait && return compare_points(geo1, geo2)
GI.geomtrait(geo1) isa LineStringTrait && return compare_lines(geo1, geo2)

error("Cant compare $(GI.trait(geo1)) and $(GI.trait(geo2)) yet")
end

function compare_points(p1, p2)
length(p1) !== length(p2) && return false

for i in eachindex(p1)
round(p1[i]; digits=10) !== round(p2[i]; digits=10) && return false
_equals(::T, geo1, ::T, geo2) where T = error("Cant compare $T yet")
function _equals(::T, p1, ::T, p2) where {T<:PointTrait}
GI.ncoord(p1) == GI.ncoord(p2) || return false
GI.x(p1) == GI.x(p2) || return false
GI.y(p1) == GI.y(p2) || return false
if GI.is3d(p1)
GI.z(p1) == GI.z(p2) || return false
end

return true
end
function _equals(::T, l1, ::T, l2) where {T<:AbstractCurveTrait}
# Check line lengths match
GI.npoint(l1) == GI.npoint(l2) || return false

function compare_lines(p1::Vector, p2::Vector)
# TODO: complete this
length(p1[1]) !== length(p2[1]) && return false
# Then check all points are the same
for (p1, p2) in zip(GI.getpoint(l1), GI.getpoint(l2))
equals(p1, p2) || return false
end
return true
end
_equals(t1, geo1, t2, geo2) = false

# """
# parallel(line1::LineString, line2::LineString)::Bool
# isparallel(line1::LineString, line2::LineString)::Bool

# Return `true` if each segment of `line1` is parallel to the correspondent segment of `line2`

Expand Down Expand Up @@ -148,17 +153,20 @@ end


"""
point_on_line(point::Point, line::LineString, ignoreEndVertices::Bool=false)::Bool
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 = GI.Point(1, 1)
point = (1, 1)
line = GI.LineString([(0, 0), (3, 3), (4, 4)])
GO.point_on_line(point, line)
# output
true
```
Expand All @@ -167,25 +175,25 @@ function point_on_line(point, line; ignore_end_vertices::Bool=false)::Bool
line_points = tuple_points(line)
n = length(line_points)

ignore = :none
exclude_boundary = :none
for i in 1:n - 1
if ignore_end_vertices == true
if ignore_end_vertices
if i === 1
ignore = :start
exclude_boundary = :start
elseif i === n - 2
ignore = :end
exclude_boundary = :end
elseif (i === 1 && i + 1 === n - 1)
ignore = :both
exclude_boundary = :both
end
end
if point_on_segment(line_points[i], line_points[i + 1], point, ignore)
if point_on_segment(point, (line_points[i], line_points[i + 1]); exclude_boundary)
return true
end
end
return false
end

function point_on_segment(start, stop, point, exclude_boundary::Symbol=:none)::Bool
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)
Expand Down Expand Up @@ -225,46 +233,72 @@ function point_on_segment(start, stop, point, exclude_boundary::Symbol=:none)::B
end

"""
point_in_polygon(point::Point, polygon::Union{Polygon, MultiPolygon}, ignoreBoundary::Bool=false)::Bool
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]]])
poly = GI.Polygon([[(-81, 41), (-81, 47), (-72, 47), (-72, 41), (-81, 41)]])
GO.point_in_polygon(point, poly)
# output
true
```
"""
function point_in_polygon(p, polygon, ignore_boundary::Bool=false)::Bool
GI.trait(polygon) isa PolygonTrait || throw(ArgumentError("Not a polygon"))

point_in_extent(p, GI.extent(polygon)) || return false
point_in_ring(p, GI.getexterior(polygon), ignore_boundary) || return false

for ring in GI.gethole(polygon)
point_in_ring(pt, ring, !ignore_boundary) && return false
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
# Cheaply check that the point is inside the ring extent
if check_extent
point_in_extent(point, GI.extent(ring)) || return false
end

function point_in_ring(pt, ring, ignore_boundary::Bool=false)
GI.trait(ring) isa Union{LineStringTrait,LinearRingTrait} || throw(ArgumentError("Not a ring"))
# Then check the point is inside the ring
inside = false
n = GI.npoint(ring)
p1 = first(GI.getpoint(ring))
p_start = GI.getpoint(ring, 1)
p_end = GI.getpoint(ring, n)

l = if GI.x(p1) == GI.x(p_end) && GI.y(p1) == GI.y(p_end)
l = n -1
# Handle closed on non-closed rings
l = if GI.x(p_start) == GI.x(p_end) && GI.y(p_start) == GI.y(p_end)
l = n - 1
else
n
end

# Loop over all points in the ring
for i in 1:l - 1
j = i + 1

Expand Down Expand Up @@ -292,103 +326,58 @@ function point_in_ring(pt, ring, ignore_boundary::Bool=false)
end

function point_in_extent(p, extent::Extents.Extent)
extent.X[1] <= GI.x(p) && extent.Y[1] <= GI.y(p) &&
extent.X[2] >= GI.x(p) && extent.Y[2] >= GI.y(p)
end

function line_in_polygon(poly, line)
out = false

polybox = bbox(poly)
linebox = bbox(line)

!(bboxOverlap(polybox, linebox)) && return false

coords = line.coordinates

for i in 1:length(coords) - 1
mid = [(coords[i][1] + coords[i + 1][1]) / 2, (coords[i][2] + coords[i + 1][2]) / 2]
if point_in_polygon(Point(mid), poly, true)
out = true
break
end
end
return out
(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(::LineStringTrait, line, ::PolygonTrait, poly)
polybox = bbox(poly)
linebox = bbox(line)

!(bboxOverlap(polybox, linebox)) && return false
function line_in_polygon(
::AbstractCurveTrait, line,
::Union{AbstractPolygonTrait,LinearRingTrait}, poly
)
Extents.intersects(GI.extent(poly), GI.extent(line)) || return false

coords = line.coordinates
inside = false

for i in 1:length(coords) - 1
!(point_in_polygon(Point(coords[i]), poly)) && return false
!inside && (inside = point_in_polygon(Point(coords[i]), poly, true))
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 = [(coords[i][1] + coords[i + 1][1]) / 2, (coords[i][2] + coords[i + 1][2]) / 2]
inside = point_in_polygon(Point(mid), poly, true)
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

# TODO - why were there two methods for this in Turf.jl?
function polygon_in_polygon(ft1, ft2, reverse::Bool=false)
polybox1 = bbox(ft1)
polybox2 = bbox(ft2)
coords = []

if reverse
!(bbox_overlap(polybox2, polybox1)) && return false

for point in GI.getpoint(ft1)
!(point_in_polygon(point, ft2)) && return false
end
else
!(bbox_overlap(polybox1, polybox2)) && return false

for point in GI.getpoint(ft2)
!(point_in_polygon(point, ft1)) && return false
end
end

return true
end
function poly_in_poly(poly1, poly2)

for point in GI.getpoint(poly1)
(point_in_polygon(point, poly2)) && return true
end

for point in GI.getpoint(poly2)
(point_in_polygon(point, poly1)) && return true
end

inter = line_intersects(polygon_to_line(poly1), polygon_to_line(poly2))
inter != nothing && return true
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

return false
# Check all points in poly1 are in poly2
for point in GI.getpoint(poly1)
point_in_polygon(point, poly2) || return false
end

end

function bbox_overlap(box1::Vector{T}, box2::Vector{T}) where {T <: Real}
box1[1] > box2[1] && return false
box1[3] < box2[3] && return false
box1[2] > box2[2] && return false
box1[4] < box2[4] && return false
return true
end
# Check the line of poly1 does not intersect the line of poly2
line_intersects(poly1, poly2) && return false

# poly1 must be in poly2
return true
end
Loading

0 comments on commit 98bfd3b

Please sign in to comment.