Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix type stability of clipping methods #74

Merged
merged 8 commits into from
Mar 10, 2024
Merged

fix type stability of clipping methods #74

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
1aff5d1
fix difference
rafaqz Mar 7, 2024
a5374c3
switch all Targets to constructed objects
rafaqz Mar 7, 2024
5ea2972
fix tests
rafaqz Mar 7, 2024
a997fed
fix docs
rafaqz Mar 7, 2024
b95f1ca
TraitTarget
rafaqz Mar 10, 2024
5060a9c
fix intersection trait
rafaqz Mar 10, 2024
d200a02
Update src/primitives.jl
rafaqz Mar 10, 2024
70685aa
Update src/primitives.jl
rafaqz Mar 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/GeometryOps.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ module GeometryOps

using GeoInterface
using GeometryBasics
using GeometryBasics.StaticArrays
import Proj
using LinearAlgebra
import ExactPredicates
Expand Down
4 changes: 2 additions & 2 deletions src/methods/angles.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ geometry trait. This is also used in the implementation, since it's a lot less
work!
=#

const _ANGLE_TARGETS = Union{GI.PolygonTrait,GI.AbstractCurveTrait,GI.MultiPointTrait,GI.PointTrait}
const _ANGLE_TARGETS = TraitTarget{Union{GI.PolygonTrait,GI.AbstractCurveTrait,GI.MultiPointTrait,GI.PointTrait}}()

"""
angles(geom, ::Type{T} = Float64)
Expand Down Expand Up @@ -169,4 +169,4 @@ function _diffs_calc_angle(::Type{T}, (Δx_prev, Δy_prev), (Δx_curr, Δy_curr)
val = clamp(dot_prod / (prev_mag * curr_mag), -one(T), one(T))
angle = real(acos(val) * 180 / π)
return angle * (cross_prod < 0 ? -1 : 1)
end
end
2 changes: 1 addition & 1 deletion src/methods/area.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ This is why signed area is only implemented for polygons.
=#

# Targets for applys functions
const _AREA_TARGETS = Union{GI.PolygonTrait,GI.AbstractCurveTrait,GI.MultiPointTrait,GI.PointTrait}
const _AREA_TARGETS = TraitTarget{Union{GI.PolygonTrait,GI.AbstractCurveTrait,GI.MultiPointTrait,GI.PointTrait}}()

"""
area(geom, [T = Float64])::T
Expand Down
2 changes: 1 addition & 1 deletion src/methods/centroid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ end
Returns the centroid and area of a given geometry.
"""
function centroid_and_area(geom, ::Type{T}=Float64; threaded=false) where T
target = Union{GI.PolygonTrait,GI.LineStringTrait,GI.LinearRingTrait}
target = TraitTarget{Union{GI.PolygonTrait,GI.LineStringTrait,GI.LinearRingTrait}}()
init = (zero(T), zero(T)), zero(T)
applyreduce(_combine_centroid_and_area, target, geom; threaded, init) do g
_centroid_and_area(GI.trait(g), g, T)
Expand Down
4 changes: 2 additions & 2 deletions src/methods/clipping/coverage.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ with a linear ring.
=#

# Targets for applys functions
const _COVERAGE_TARGETS = Union{GI.PolygonTrait,GI.AbstractCurveTrait,GI.MultiPointTrait,GI.PointTrait}
const _COVERAGE_TARGETS = TraitTarget{Union{GI.PolygonTrait,GI.AbstractCurveTrait,GI.MultiPointTrait,GI.PointTrait}}()

# Wall types for coverage
const UNKNOWN, NORTH, EAST, SOUTH, WEST = 0:4
Expand Down Expand Up @@ -298,4 +298,4 @@ function _partial_edge_out_area((x1, y1), xmin, xmax, ymin, ymax, wall)
x_wall = (wall == NORTH || wall == EAST) ? xmax : xmin
y_wall = (wall == NORTH || wall == WEST) ? ymax : ymin
return x1 * y_wall - x_wall * y1
end
end
12 changes: 6 additions & 6 deletions src/methods/clipping/difference.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,16 +27,16 @@ GI.coordinates.(diff_poly)
```
"""
function difference(
geom_a, geom_b, ::Type{T} = Float64; target::Target = nothing,
) where {T <: AbstractFloat, Target <: Union{Nothing, GI.AbstractTrait}}
return _difference(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
geom_a, geom_b, ::Type{T} = Float64; target=nothing,
) where {T<:AbstractFloat}
return _difference(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
end
rafaqz marked this conversation as resolved.
Show resolved Hide resolved

#= The 'difference' function returns the difference of two polygons as a list of polygons.
The algorithm to determine the difference was adapted from "Efficient clipping of efficient
polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
function _difference(
::Type{GI.PolygonTrait}, ::Type{T},
::TraitTarget{GI.PolygonTrait}, ::Type{T},
::GI.PolygonTrait, poly_a,
::GI.PolygonTrait, poly_b,
) where T
Expand Down Expand Up @@ -65,7 +65,7 @@ function _difference(
if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
_add_holes_to_polys!(T, polys, GI.gethole(poly_a))
for hole in GI.gethole(poly_b)
new_polys = intersection(GI.Polygon([hole]), poly_a, T; target = GI.PolygonTrait())
new_polys = intersection(GI.Polygon([hole]), poly_a, T; target = GI.PolygonTrait)
if length(new_polys) > 0
append!(polys, new_polys)
end
Expand All @@ -76,7 +76,7 @@ end

# Many type and target combos aren't implemented
function _difference(
::Type{Target}, ::Type{T},
::TraitTarget{Target}, ::Type{T},
trait_a::GI.AbstractTrait, geom_a,
trait_b::GI.AbstractTrait, geom_b,
) where {Target, T}
Expand Down
10 changes: 5 additions & 5 deletions src/methods/clipping/intersection.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@

## Example

```jldoctest

Check failure on line 16 in src/methods/clipping/intersection.jl

View workflow job for this annotation

GitHub Actions / Documentation 

doctest failure in ~/work/GeometryOps.jl/GeometryOps.jl/src/methods/clipping/intersection.jl:16-27

```jldoctest
import GeoInterface as GI, GeometryOps as GO

line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
inter_points = GO.intersection(line1, line2; target = GI.PointTrait())
GI.coordinates.(inter_points)

# output
1-element Vector{Vector{Float64}}:
 [125.58375366067547, -14.83572303404496]
```

Subexpression:

import GeoInterface as GI, GeometryOps as GO

line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
inter_points = GO.intersection(line1, line2; target = GI.PointTrait())
GI.coordinates.(inter_points)

Evaluated output:

ERROR: AssertionError: Intersection between GeoInterface.LineTrait() and GeoInterface.LineTrait() with target GeoInterface.PointTrait isn't implemented yet.
Stacktrace:
 [1] _intersection(::GeometryOps.TraitTarget{GeoInterface.PointTrait}, ::Type{Float64}, trait_a::GeoInterface.LineTrait, geom_a::GeoInterface.Wrappers.Line{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}, trait_b::GeoInterface.LineTrait, geom_b::GeoInterface.Wrappers.Line{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing})
   @ GeometryOps ~/work/GeometryOps.jl/GeometryOps.jl/src/methods/clipping/intersection.jl:80
 [2] intersection(geom_a::GeoInterface.Wrappers.Line{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}, geom_b::GeoInterface.Wrappers.Line{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}, ::Type{Float64}; target::GeoInterface.PointTrait)
   @ GeometryOps ~/work/GeometryOps.jl/GeometryOps.jl/src/methods/clipping/intersection.jl:32
 [3] top-level scope
   @ none:1

Expected output:

1-element Vector{Vector{Float64}}:
 [125.58375366067547, -14.83572303404496]


  diff = 
    Warning: Diff output requires color.
    1-element Vector{Vector{Float64}}:
     [125.58375366067547, -14.83572303404496]ERROR: AssertionError: Intersection between GeoInterface.LineTrait() and GeoInterface.LineTrait() with target GeoInterface.PointTrait isn't implemented yet.
    Stacktrace:
     [1] _intersection(::GeometryOps.TraitTarget{GeoInterface.PointTrait}, ::Type{Float64}, trait_a::GeoInterface.LineTrait, geom_a::GeoInterface.Wrappers.Line{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}, trait_b::GeoInterface.LineTrait, geom_b::GeoInterface.Wrappers.Line{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing})
       @ GeometryOps ~/work/GeometryOps.jl/GeometryOps.jl/src/methods/clipping/intersection.jl:80
     [2] intersection(geom_a::GeoInterface.Wrappers.Line{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}, geom_b::GeoInterface.Wrappers.Line{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}, ::Type{Float64}; target::GeoInterface.PointTrait)
       @ GeometryOps ~/work/GeometryOps.jl/GeometryOps.jl/src/methods/clipping/intersection.jl:32
     [3] top-level scope
       @ none:1
import GeoInterface as GI, GeometryOps as GO

line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
Expand All @@ -27,9 +27,9 @@
```
"""
function intersection(
geom_a, geom_b, ::Type{T} = Float64; target::Target = nothing,
) where {T <: AbstractFloat, Target <: Union{Nothing, GI.AbstractTrait}}
return _intersection(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
geom_a, geom_b, ::Type{T}=Float64; target=nothing,
) where {T<:AbstractFloat}
return _intersection(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
end

# Curve-Curve Intersections with target Point
Expand All @@ -45,7 +45,7 @@
of efficient polygons," by Greiner and Hormann (1998).
DOI: https://doi.org/10.1145/274363.274364 =#
function _intersection(
::Type{GI.PolygonTrait}, ::Type{T},
::TraitTarget{GI.PolygonTrait}, ::Type{T},
::GI.PolygonTrait, poly_a,
::GI.PolygonTrait, poly_b,
) where {T}
Expand Down Expand Up @@ -73,7 +73,7 @@

# Many type and target combos aren't implemented
function _intersection(
::Type{Target}, ::Type{T},
::TraitTarget{Target}, ::Type{T},
trait_a::GI.AbstractTrait, geom_a,
trait_b::GI.AbstractTrait, geom_b,
) where {Target, T}
Expand Down
13 changes: 7 additions & 6 deletions src/methods/clipping/union.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,16 +27,16 @@ GI.coordinates.(union_poly)
```
"""
function union(
geom_a, geom_b, ::Type{T} = Float64; target::Target = nothing,
) where {T <: AbstractFloat, Target <: Union{Nothing, GI.AbstractTrait}}
_union(Target, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
geom_a, geom_b, ::Type{T}=Float64; target=nothing,
) where {T<:AbstractFloat}
_union(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
end

#= This 'union' implementation returns the union of two polygons. The algorithm to determine
the union was adapted from "Efficient clipping of efficient polygons," by Greiner and
Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
function _union(
::Type{GI.PolygonTrait}, ::Type{T},
::TraitTarget{GI.PolygonTrait}, ::Type{T},
::GI.PolygonTrait, poly_a,
::GI.PolygonTrait, poly_b,
) where T
Expand Down Expand Up @@ -74,12 +74,13 @@ function _union(
return polys
end


# Many type and target combos aren't implemented
function _union(
::Type{Target}, ::Type{T},
::TraitTarget{Target}, ::Type{T},
trait_a::GI.AbstractTrait, geom_a,
trait_b::GI.AbstractTrait, geom_b,
) where {Target, T}
) where {Target,T}
throw(ArgumentError("Union between $trait_a and $trait_b with target $Target isn't implemented yet."))
return nothing
end
2 changes: 1 addition & 1 deletion src/methods/distance.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ Also note that singed_distance only makes sense for "filled-in" shapes, like
polygons, so it isn't implemented for curves.
=#

const _DISTANCE_TARGETS = Union{GI.AbstractPolygonTrait,GI.LineStringTrait,GI.LinearRingTrait,GI.LineTrait,GI.PointTrait}
const _DISTANCE_TARGETS = TraitTarget{Union{GI.AbstractPolygonTrait,GI.LineStringTrait,GI.LinearRingTrait,GI.LineTrait,GI.PointTrait}}()

"""
distance(point, geom, ::Type{T} = Float64)::T
Expand Down
Loading
Loading