Add benchmark + refer to Predicates directly

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

@@ -5,8 +5,9 @@ module GeometryOps
 using GeoInterface
 using GeometryBasics
 import Tables
-using ExactPredicates, LinearAlgebra, Statistics
+using LinearAlgebra, Statistics
 import GeometryBasics.StaticArrays
+import ExactPredicates
 import Base.@kwdef
 
 using GeoInterface.Extents: Extents

src/methods/clipping/intersection.jl

@@ -263,8 +263,8 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge) where T
     r_cross_s = rx * sy - ry * sx
     if Predicates.isparallel((rx, ry), (sx, sy)) != 0  # non-parallel lines
         # Calculate α ratio if lines cross or touch
-        a1_orient = ExactPredicates.orient(b1, b2, a1)
-        a2_orient = ExactPredicates.orient(b1, b2, a2)
+        a1_orient = Predicates.orient(b1, b2, a1)
+        a2_orient = Predicates.orient(b1, b2, a2)
         α = if a1_orient == 0  # α = 0
             zero(T)
         elseif a2_orient == 0  # α = 1
@@ -276,8 +276,8 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge) where T
             return line_orient, intr1, intr2
         end
         # Calculate β ratio if lines touch or cross
-        b1_orient = ExactPredicates.orient(a1, a2, b1)
-        b2_orient = ExactPredicates.orient(a1, a2, b2)
+        b1_orient = Predicates.orient(a1, a2, b1)
+        b2_orient = Predicates.orient(a1, a2, b2)
         β = if b1_orient == 0  # β = 0
             zero(T)
         elseif b2_orient == 0  # β = 1
@@ -295,10 +295,10 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge) where T
         line_orient = (α == 0 || α == 1 || β == 0 || β == 1) ? line_hinge : line_cross
     elseif Predicates.iscollinear((Δqp_x, Δqp_y), (sx, sy)) == 0 # collinear parallel lines
         # Determine if lines touch or overlap and with what α and β values
-        a1_side = ExactPredicates.sameside(a1, b1, b2)
-        a2_side = ExactPredicates.sameside(a2, b1, b2)
-        b1_side = ExactPredicates.sameside(b1, a1, a2)
-        b2_side = ExactPredicates.sameside(b2, a1, a2)
+        a1_side = Predicates.sameside(a1, b1, b2)
+        a2_side = Predicates.sameside(a2, b1, b2)
+        b1_side = Predicates.sameside(b1, a1, a2)
+        b2_side = Predicates.sameside(b2, a1, a2)
         # Lines touch or overlap if endpoints of line a are on/in line b and visa versa
         r_dot_s = rx * sx + ry * sy
         # Determine which endpoints start and end the overlapping region

src/methods/clipping/predicates.jl

@@ -6,6 +6,9 @@ module Predicates
 
     using ExactPredicates: ext, inp
 
+    orient(args...) = ExactPredicates.orient(args...)
+    sameside(args...) = ExactPredicates.sameside(args...)
+
     # This is the implementation of r_cross_s from 2 points in the `_intersection_points` file.
     # 0 == parallel
     function isparallel(a1, a2, b1, b2)