Add a convex hull algorithm

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.10"
 
 [deps]
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
+DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
@@ -26,6 +27,7 @@ GeometryOpsProjExt = "Proj"
 
 [compat]
 CoordinateTransformations = "0.5, 0.6"
+DelaunayTriangulation = "1.0.4"
 ExactPredicates = "2.2.8"
 FlexiJoins = "0.1.30"
 GeoInterface = "1.2"

ext/GeometryOpsLibGEOSExt/simple_overrides.jl

@@ -64,3 +64,17 @@ function GO.intersects(::GEOS, geom_a, geom_b)
     return LG.intersects(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
 end
 
+# ## Convex hull
+function GO.convex_hull(::GEOS, geoms)
+    return LG.convexhull(
+        LG.MultiPoint(
+            collect(
+                GO.flatten(
+                    x -> GI.convert(LG, x), 
+                    GI.PointTrait, 
+                    geoms
+                )
+            )
+        )
+    )
+end

src/GeometryOps.jl

@@ -4,13 +4,14 @@ module GeometryOps
 
 using GeoInterface
 using GeometryBasics
-import Tables
 using LinearAlgebra, Statistics
+
+import Tables
 import GeometryBasics.StaticArrays
+import DelaunayTriangulation # for convex hull and triangulation
 import ExactPredicates
 import Base.@kwdef
-
-using GeoInterface.Extents: Extents
+import GeoInterface.Extents: Extents
 
 const GI = GeoInterface
 const GB = GeometryBasics
@@ -28,6 +29,7 @@ include("methods/area.jl")
 include("methods/barycentric.jl")
 include("methods/buffer.jl")
 include("methods/centroid.jl")
+include("methods/convex_hull.jl")
 include("methods/distance.jl")
 include("methods/equals.jl")
 include("methods/clipping/predicates.jl")

src/methods/convex_hull.jl

@@ -0,0 +1,131 @@
+#=
+# Convex hull
+
+The [_convex hull_](https://en.wikipedia.org/wiki/Convex_hull) of a set of points is the smallest [**convex**](https://en.wikipedia.org/wiki/Convex_set) polygon that contains all the points.
+
+GeometryOps.jl provides a number of methods for computing the convex hull of a set of points, usually
+linked to other Julia packages.  
+
+For now, we expose one algorithm, [MonotoneChainMethod](@ref), which uses the [DelaunayTriangulation.jl](https://github.com/JuliaGeometry/DelaunayTriangulation.jl) 
+package.  The `GEOS()` interface also supports convex hulls.  
+
+Future work could include other algorithms, such as [Quickhull.jl](https://github.com/augustt198/Quickhull.jl), or similar, via package extensions.
+
+
+## Example
+
+### Simple hull
+```@example simple
+import GeometryOps as GO, GeoInterface as GI
+using CairoMakie # to plot
+
+points = randn(GO.Point2f, 100)
+f, a, p = plot(points; label = "Points")
+hull_poly = GO.convex_hull(points)
+lines!(a, hull_poly; label = "Convex hull", color = Makie.wong_colors()[2])
+axislegend(a)
+f
+```
+
+## Convex hull of the USA
+```@example usa
+import GeometryOps as GO, GeoInterface as GI
+using CairoMakie # to plot
+using NaturalEarth # for data
+
+all_adm0 = naturalearth("admin_0_countries", 110)
+usa = all_adm0.geometry[findfirst(==("USA"), all_adm0.ADM0_A3)]
+f, a, p = lines(usa)
+lines!(a, GO.convex_hull(usa); color = Makie.wong_colors()[2])
+f
+```
+
+## Investigating the winding order
+
+The winding order of the monotone chain method is counterclockwise,
+while the winding order of the GEOS method is clockwise.
+
+GeometryOps' convexity detection says that the GEOS hull is convex,
+while the monotone chain method hull is not.  However, they are both going
+over the same points (we checked), it's just that the winding order is different.
+
+In reality, both sets are convex, but we need to fix the GeometryOps convexity detector 
+([`isconcave`](@ref))!
+
+We may also decide at a later date to change the returned winding order of the polygon, but
+most algorithms are robust to that, and you can always [`fix`](@ref) it...
+
+```@example windingorder
+import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+using CairoMakie # to plot
+
+points = rand(Point2{Float64}, 100)
+go_hull = GO.convex_hull(GO.MonotoneChainMethod(), points)
+lg_hull = GO.convex_hull(GO.GEOS(), points)
+
+fig = Figure()
+a1, p1 = lines(fig[1, 1], go_hull; color = 1:GI.npoint(go_hull), axis = (; title = "MonotoneChainMethod()"))
+a2, p2 = lines(fig[2, 1], lg_hull; color = 1:GI.npoint(lg_hull), axis = (; title = "GEOS()"))
+cb = Colorbar(fig[1:2, 2], p1; label = "Vertex number")
+fig
+```
+
+## Implementation
+
+=#
+
+"""
+    convex_hull([method], geometries)
+
+Compute the convex hull of the points in `geometries`.  
+Returns a `GI.Polygon` representing the convex hull.
+
+Note that the polygon returned is wound counterclockwise
+as in the Simple Features standard by default.  If you 
+choose GEOS, the winding order will be inverted.
+
+!!! warning
+    This interface only computes the 2-dimensional convex hull!
+
+    For higher dimensional hulls, use the relevant package (Qhull.jl, Quickhull.jl, or similar).
+"""
+function convex_hull end
+
+"""
+    MonotoneChainMethod()
+
+This is an algorithm for the [`convex_hull`](@ref) function.
+
+Uses [`DelaunayTriangulation.jl`](https://github.com/JuliaGeometry/DelaunayTriangulation.jl) to compute the convex hull.
+This is a pure Julia algorithm which provides an optimal Delaunay triangulation.
+
+See also [`convex_hull`](@ref)
+"""
+struct MonotoneChainMethod end
+
+# GrahamScanMethod, etc. can be implemented in GO as well, if someone wants to.
+# If we add an extension on Quickhull.jl, then that would be another algorithm.
+
+convex_hull(geometries) = convex_hull(MonotoneChainMethod(), geometries)
+
+# TODO: have this respect the CRS by pulling it out of `geometries`.
+function convex_hull(::MonotoneChainMethod, geometries)
+    # Extract all points as tuples.  We have to collect and allocate
+    # here, because DelaunayTriangulation only accepts vectors of 
+    # point-like geoms.
+
+    # Cleanest would be to use the iterable from GO.flatten directly,
+    # but that would require us to implement the convex hull algorithm
+    # directly.
+
+    # TODO: create a specialized method that extracts only the information
+    # required, GeometryBasics points can be passed through directly.
+    points = collect(flatten(tuples, GI.PointTrait, geometries))
+    # Compute the convex hull using DelTri (shorthand for DelaunayTriangulation.jl).
+    hull = DelaunayTriangulation.convex_hull(points)
+    # Convert the result to a `GI.Polygon` and return it.
+    # View would be more efficient here, but re-allocating
+    # is cleaner.
+    point_vec = DelaunayTriangulation.get_points(hull)[DelaunayTriangulation.get_vertices(hull)]
+    return GI.Polygon([GI.LinearRing(point_vec)])
+end

test/methods/convex_hull.jl

@@ -0,0 +1,42 @@
+using Test
+import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+
+@testset "Basic example" begin
+    points = tuple.(rand(100), rand(100))
+    hull = GO.convex_hull(points)
+    @test !GO.isconcave(hull) skip=true # TODO: fix
+    @test !GO.isclockwise(GI.getexterior(hull)) # exterior should be ccw
+    @test all(x -> GO.covers(hull, x), points) # but the orientation is right
+    # Test against LibGEOS, by testing that all the points are the same
+    # This is robust against winding order and starting at a different point, etc.
+    @test isempty(
+        setdiff(
+            collect(GO.flatten(GO.tuples, GI.PointTrait, hull)), 
+            collect(GO.flatten(GO.tuples, GI.PointTrait, GO.convex_hull(GEOS(), points)))
+        )
+    )
+end
+
+@testset "Duplicated points" begin
+    points = tuple.(rand(100), rand(100))
+    @test_nowarn hull = GO.convex_hull(vcat(points, points))
+    single_hull = GO.convex_hull(points)
+    double_hull = GO.convex_hull(vcat(points, points))
+
+    @test GO.equals(GI.getexterior(single_hull), GI.getexterior(double_hull))
+    @test !GO.isconcave(double_hull) skip=true # TODO: fix
+    @test !GO.isclockwise(GI.getexterior(double_hull)) # exterior should be ccw
+    @test all(x -> GO.covers(single_hull, x), points)
+    @test all(x -> GO.covers(double_hull, x), points)
+    # Test against LibGEOS, by testing that all the points are the same
+    # This is robust against winding order and starting at a different point, etc.
+    @test isempty(
+        setdiff(
+            collect(GO.flatten(GO.tuples, GI.PointTrait, double_hull)), 
+            collect(GO.flatten(GO.tuples, GI.PointTrait, GO.convex_hull(GEOS(), points)))
+        )
+    )
+end
+
+# The rest of the tests are handled in DelaunayTriangulation.jl, this is simply testing
+# that the methods work as expected.

test/runtests.jl

@@ -22,6 +22,7 @@ const GO = GeometryOps
     @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
     @testset "Orientation" begin include("methods/orientation.jl") end
     @testset "Centroid" begin include("methods/centroid.jl") end
+    @testset "Convex Hull" begin include("methods/convex_hull.jl") end
     @testset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end
     @testset "Distance" begin include("methods/distance.jl") end
     @testset "Equals" begin include("methods/equals.jl") end