Sg/clipping errors

Project.toml

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

benchmarks/benchmark_plots.jl

@@ -89,7 +89,6 @@ function plot_trials(
     end
 
     tag_attrs = capture_tag_attrs(results.tags)
-    tag_theme = merge(theme, tag_attrs)
 
     lp = if legend_position isa Makie.Automatic 
         gp.layout[gp.span.rows, gp.span.cols, TopRight()] 
@@ -101,12 +100,15 @@ function plot_trials(
         error()
     end
 
-    return Makie.with_theme(tag_theme) do
+    ax = Makie.with_theme(theme) do
         ax = Axis(
             gp;
             tag_attrs.Axis...,
             xlabel = "Number of points", ylabel = "Time to calculate",
             xscale = log10, yscale = log10, ytickformat = _prettytime,
+            xticksvisible = true, xticklabelsvisible = true,
+            yticks = Makie.LogTicks(Makie.WilkinsonTicks(7; k_min = 4)),
+            ygridwidth = 0.75,
         ) 
         plots = [scatterlines!(ax, x, y; label = label) for (x, y, label) in zip(xs, ys, labels)]
         setproperty!.(getindex.(getproperty.(plots, :plots), 1), :alpha, 0.1)
@@ -118,13 +120,11 @@ function plot_trials(
             valign = legend_valign, 
             orientation = legend_orientation
         )
-        ax.xticksvisible[] = true
         ax.xtickcolor[] = ax.xgridcolor[]
-        ax.xticklabelsvisible[] = true
-        ax.yticks[] = Makie.LogTicks(Makie.WilkinsonTicks(7; k_min = 4))
-        ax.ygridwidth[] = 0.75
-        return ax
+        ax
     end
+
+    return ax
 end
 
 const _tag_includelist = ["title", "subtitle"]

benchmarks/polygon_scaling.jl

@@ -0,0 +1,19 @@
+using CairoMakie, Chairmarks, BenchmarkTools
+import GeometryOps as GO, GeoInterface as GI, LibGEOS as LG
+include("benchmark_plots.jl")
+
+p25 = GI.Polygon([[(0.0, 0.0), (5.0, 0.0), (5.0, 8.0), (0.0, 8.0), (0.0, 0.0)], [(4.0, 0.5), (4.5, 0.5), (4.5, 3.5), (4.0, 3.5), (4.0, 0.5)], [(2.0, 4.0), (4.0, 4.0), (4.0, 6.0), (2.0, 6.0), (2.0, 4.0)]])
+p26 = GI.Polygon([[(3.0, 1.0), (8.0, 1.0), (8.0, 7.0), (3.0, 7.0), (3.0, 5.0), (6.0, 5.0), (6.0, 3.0), (3.0, 3.0), (3.0, 1.0)], [(3.5, 5.5), (6.0, 5.5), (6.0, 6.5), (3.5, 6.5), (3.5, 5.5)],  [(5.5, 1.5), (5.5, 2.5), (3.5, 2.5), (3.5, 1.5), (5.5, 1.5)]])
+
+suite = BenchmarkGroup(["title:Polygon intersection timing","subtitle:Single polygon, densified"])
+
+for max_distance in exp10.(LinRange(-1, 1.5, 10))
+    p25s = GO.segmentize(p25; max_distance)
+    p26s = GO.segmentize(p26; max_distance)
+    n_verts = GI.npoint(p25s)
+    suite["GeometryOps"][n_verts] = @be GO.intersection($p25s, $p26s; target = $(GI.PolygonTrait()), fix_multipoly = $nothing)
+    suite["LibGEOS"][n_verts] = @be LG.intersection($(GI.convert(LG, p25s)), $(GI.convert(LG, p26s)))
+end
+
+
+plot_trials(suite)

src/GeometryOps.jl

@@ -7,6 +7,7 @@ using GeometryBasics
 import Tables
 using LinearAlgebra, Statistics
 import GeometryBasics.StaticArrays
+import ExactPredicates
 import Base.@kwdef
 
 using GeoInterface.Extents: Extents
@@ -26,6 +27,7 @@ include("methods/barycentric.jl")
 include("methods/centroid.jl")
 include("methods/distance.jl")
 include("methods/equals.jl")
+include("methods/clipping/predicates.jl")
 include("methods/clipping/clipping_processor.jl")
 include("methods/clipping/coverage.jl")
 include("methods/clipping/cut.jl")

src/methods/clipping/coverage.jl

@@ -70,23 +70,23 @@ function coverage(geom, cell_ext::Extents.Extent, ::Type{T} = Float64; threaded=
 end
 
 # Points, MultiPoints, Curves, MultiCurves
-_coverage(::Type{T}, ::GI.AbstractGeometryTrait, geom, xmin, xmax, ymin, ymax) where T = zero(T)
+_coverage(::Type{T}, ::GI.AbstractGeometryTrait, geom, xmin, xmax, ymin, ymax; kwargs...) where T = zero(T)
 
 # Polygons
-function _coverage(::Type{T}, ::GI.PolygonTrait, poly, xmin, xmax, ymin, ymax) where T
+function _coverage(::Type{T}, ::GI.PolygonTrait, poly, xmin, xmax, ymin, ymax; exact = _False()) where T
     GI.isempty(poly) && return zero(T)
-    cov_area = _coverage(T, GI.getexterior(poly), xmin, xmax, ymin, ymax)
+    cov_area = _coverage(T, GI.getexterior(poly), xmin, xmax, ymin, ymax; exact)
     cov_area == 0 && return cov_area
     # Remove hole coverage from total
     for hole in GI.gethole(poly)
-        cov_area -= _coverage(T, hole, xmin, xmax, ymin, ymax)
+        cov_area -= _coverage(T, hole, xmin, xmax, ymin, ymax; exact)
     end
     return cov_area
 end
 
 #= Calculates the area of the filled ring within the cell defined by corners with (xmin, ymin),
 (xmin, ymax), (xmax, ymax), and (xmax, ymin). =#
-function _coverage(::Type{T}, ring, xmin, xmax, ymin, ymax) where T
+function _coverage(::Type{T}, ring, xmin, xmax, ymin, ymax; exact) where T
     cov_area = zero(T)
     unmatched_out_wall, unmatched_out_point = UNKNOWN, (zero(T), zero(T))
     unmatched_in_wall, unmatched_in_point = unmatched_out_wall, unmatched_out_point
@@ -137,7 +137,7 @@ function _coverage(::Type{T}, ring, xmin, xmax, ymin, ymax) where T
             else
                 check_point = find_point_on_cell(unmatched_out_point, start_point,
                     unmatched_out_wall, start_wall,xmin, xmax, ymin, ymax)
-                if _point_filled_curve_orientation(check_point, ring; in = true, on = false, out = false)
+                if _point_filled_curve_orientation(check_point, ring; in = true, on = false, out = false, exact)
                     cov_area += connect_edges(T, unmatched_out_point, start_point,
                         unmatched_out_wall, start_wall,xmin, xmax, ymin, ymax)
                 else
@@ -160,7 +160,7 @@ function _coverage(::Type{T}, ring, xmin, xmax, ymin, ymax) where T
     cov_area = abs(cov_area) / 2
     #  if grid cell is within polygon then the area is grid cell area
     if cov_area == 0
-        if _point_filled_curve_orientation((xmin, ymin), ring; in = true, on = true, out = false)
+        if _point_filled_curve_orientation((xmin, ymin), ring; in = true, on = true, out = false, exact)
             cov_area = abs((xmax - xmin) * (ymax - ymin))
         end
     end

src/methods/clipping/cut.jl

@@ -59,33 +59,33 @@ GI.coordinates.(cut_polys)
 ```
 """
 cut(geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} =
-    _cut(T, GI.trait(geom), geom, GI.trait(line), line)
+    _cut(T, GI.trait(geom), geom, GI.trait(line), line; exact = _True())
 
 #= Cut a given polygon by given line. Add polygon holes back into resulting pieces if there
 are any holes. =#
-function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line) where T
+function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) where T
     ext_poly = GI.getexterior(poly)
-    poly_list, intr_list = _build_a_list(T, ext_poly, line)
+    poly_list, intr_list = _build_a_list(T, ext_poly, line; exact)
     n_intr_pts = length(intr_list)
     # If an impossible number of intersection points, return original polygon
     if n_intr_pts < 2 || isodd(n_intr_pts)
         return [tuples(poly)]
     end
     # Cut polygon by line
-    cut_coords = _cut(T, ext_poly, line, poly_list, intr_list, n_intr_pts)
+    cut_coords = _cut(T, ext_poly, line, poly_list, intr_list, n_intr_pts; exact)
     # Close coords and create polygons
     for c in cut_coords
         push!(c, c[1])
     end
     cut_polys = [GI.Polygon([c]) for c in cut_coords]
     # Add original polygon holes back in
     remove_idx = falses(length(cut_polys))
-    _add_holes_to_polys!(T, cut_polys, GI.gethole(poly), remove_idx)
+    _add_holes_to_polys!(T, cut_polys, GI.gethole(poly), remove_idx; exact)
     return cut_polys
 end
 
 # Many types aren't implemented
-function _cut(_, trait::GI.AbstractTrait, geom, line)
+function _cut(::Type{T}, trait::GI.AbstractTrait, geom, line; kwargs...) where T
     @assert(
         false,
         "Cutting of $trait isn't implemented yet.",
@@ -97,10 +97,10 @@ end
 of cut geometry in Vector{Vector{Tuple}} format. 
 
 Note: degenerate cases where intersection points are vertices do not work right now. =#
-function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts) where T
+function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where T
     # Sort and catagorize the intersection points
     sort!(intr_list, by = x -> geom_list[x].fracs[2])
-    _flag_ent_exit!(GI.LineTrait(), line, geom_list)
+    _flag_ent_exit!(GI.LineTrait(), line, geom_list; exact)
     # Add first point to output list
     return_coords = [[geom_list[1].point]]
     cross_backs = [(T(Inf),T(Inf))]

src/methods/clipping/difference.jl

@@ -35,7 +35,10 @@ GI.coordinates.(diff_poly)
 function difference(
     geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...,
 ) where {T<:AbstractFloat}
-    return _difference(TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; kwargs...)
+    return _difference(
+        TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b;
+        exact = _True(), kwargs...,
+    )
 end
 
 #= The 'difference' function returns the difference of two polygons as a list of polygons.
@@ -45,17 +48,17 @@ function _difference(
     ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
     ::GI.PolygonTrait, poly_b;
-    kwargs...
+    exact, kwargs...
 ) where T
     # Get the exterior of the polygons
     ext_a = GI.getexterior(poly_a)
     ext_b = GI.getexterior(poly_b)
     # Find the difference of the exterior of the polygons
-    a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f)
-    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _diff_step)
+    a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f; exact)
+    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _diff_step, poly_a, poly_b)
     # if no crossing points, determine if either poly is inside of the other
     if isempty(polys)
-        a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b)
+        a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact)
         # add case for if they polygons are the same (all intersection points!)
         # add a find_first check to find first non-inter poly!
         if b_in_a && !a_in_b  # b in a and can't be the same polygon
@@ -69,7 +72,7 @@ function _difference(
     remove_idx = falses(length(polys))
     # If the original polygons had holes, take that into account.
     if GI.nhole(poly_a) != 0
-        _add_holes_to_polys!(T, polys, GI.gethole(poly_a), remove_idx)
+        _add_holes_to_polys!(T, polys, GI.gethole(poly_a), remove_idx; exact)
     end
     if GI.nhole(poly_b) != 0
         for hole in GI.gethole(poly_b)
@@ -81,9 +84,7 @@ function _difference(
         end
     end
     # Remove uneeded collinear points on same edge
-    for p in polys
-        _remove_collinear_points!(p, remove_idx)
-    end
+    _remove_collinear_points!(polys, remove_idx, poly_a, poly_b)
     return polys
 end