Convex Hulls
In this section we illustrate the convex hull operation. We give examples of the symbolic implementation, and the concrete convex hull in low dimensions. We show how to test if a point lies in the convex hull of a set of points in the plane using LazySets
. Moreover, we give examples of creating the convex hull of sets whose vertices are represented as static vectors, which can dramatically improve performance in many use cases. Finally, we give an example of creating the convex hull of points in higher dimensions.
Symbolic convex hull
The lazy convex hull, ConvexHull
, is the binary operator that implements the convex hull of the union between two convex sets.
using Plots, LazySets
A = 1/sqrt(2.) * [1 -1; 1 1]
Bn = n -> BallInf(ones(n), 0.2)
X = Bn(2)
Y = CH(X, exp(A) * X)
ConvexHull{Float64, BallInf{Float64, Vector{Float64}}, LinearMap{Float64, BallInf{Float64, Vector{Float64}}, Float64, Matrix{Float64}}}(BallInf{Float64, Vector{Float64}}([1.0, 1.0], 0.2), LinearMap{Float64, BallInf{Float64, Vector{Float64}}, Float64, Matrix{Float64}}([1.5418634570456318 -1.3175384087798812; 1.3175384087798812 1.5418634570456318], BallInf{Float64, Vector{Float64}}([1.0, 1.0], 0.2)))
The name CH
is an alias for ConvexHull
, so you can use both interchangeably. This type is parametric in the operands's types.
p = plot(X, 1e-2, color="blue")
plot!(p, exp(A) * X, color="green")
plot!(p, Y, color="red", alpha=0.2)
We can as well work with a 100-dimensional set:
using SparseArrays
X = Bn(100)
A = blockdiag([sparse(A) for i in 1:50]...)
Y = CH(X, exp(Matrix(A)) * X)
dim(Y)
100
To take the convex hull of a large number of sets, there is the n
-ary type ConvexHullArray
. For instance, below we create a collection of balls b
via list comprehension, and pass them to create a new ConvexHullArray
instance.
b = [Ball2([2*pi*i/100, sin(2*pi*i/100)], 0.05) for i in 1:100];
c = ConvexHullArray(b);
plot(c, alpha=0.1, color="blue")
plot!(b, alpha=0.5, color="red")
2D convex hull
In two dimensions the convex_hull
function computes the concrete convex hull of a set of points.
points = N -> [randn(2) for i in 1:N]
v = points(30)
hull = convex_hull(v)
typeof(hull), length(v), length(hull)
(Vector{Vector{Float64}}, 30, 8)
Notice that the output is a vector of floating point numbers representing the coordinates of the points, and that the number of points in the convex hull has decreased.
We can plot both the random points and the polygon generated by the convex hull with the plot
function:
p = plot([Singleton(vi) for vi in v])
plot!(p, VPolygon(hull), alpha=0.2)
Test point in convex hull
One can check whether a point lies inside or outside of a convex hull efficiently in two dimensions, using the fact that the output of convex_hull
returns the points ordered in counter-clockwise fashion.
To check if a point p::AbstractVector
is in another set, e.g. a polygon in vertex representation V
, use p ∈ V
. However, if you are working with a Singleton
, which is a set with one element, use set inclusion ⊆
. The following example illustrates this difference.
julia> Singleton(v[1]) ∈ VPolygon(hull)
ERROR: cannot make a point-in-set check if the left-hand side is a set; either
check for set inclusion, as in `S ⊆ X`, or check for membership, as in
`element(S) ∈ X` (the results are equivalent but the implementations may differ)
As the error suggests, either use element
to access the element of the singleton and check if it belongs to the right-hand side set:
element(Singleton(v[1])) ∈ VPolygon(hull)
true
Or use set inclusion between the singleton and the right-hand side set:
Singleton(v[1]) ⊆ VPolygon(hull)
true
Let us note that one can also make the point-in-convex-hull test by solving a feasibility problem; actually, this is the fallback implementation used for in any dimension. However, the specialized approach in 2D is more efficient.
Using static vectors
Usual vectors are such that you can push!
and pop!
without changing its type: the size is not a static property. Vectors of fixed size, among other types, are provided by the StaticArrays.jl package from the JuliaArrays ecosystem. Using static arrays for vectors of "small" dimension can dramatically improve performance.
Since the convex hull algorithm supports any AbstractVector
, it can be applied with static vectors. The following example illustrates this fact.
v = points(1000)
convex_hull(points(3)) # warm-up
@time convex_hull(v)
13-element Vector{Vector{Float64}}:
[-3.1018942487027936, -0.5138625366055106]
[-2.2835214719541685, -1.7520326334993468]
[-1.7121083164595972, -2.5655428095216295]
[-0.5674027525415846, -2.992526926274017]
[0.8313656167997135, -2.5193347256443763]
[2.085945193325355, -2.0573707432814796]
[2.5869306117994584, -1.2549726253197953]
[2.9094031573517074, -0.3070263535954405]
[2.9684509920968605, 0.5953353058671058]
[2.250475430291036, 2.351070953959735]
[1.569332645014958, 2.676139601555899]
[-0.23562340136345455, 3.37989598199693]
[-3.048159617421087, 1.6798228422033785]
Now working with static vectors:
using StaticArrays
convex_hull([@SVector(rand(2)) for i in 1:3]) # warm-up
v_static = [SVector{2, Float64}(vi) for vi in v]
@time convex_hull(v_static)
13-element Vector{StaticArraysCore.SVector{2, Float64}}:
[-3.1018942487027936, -0.5138625366055106]
[-2.2835214719541685, -1.7520326334993468]
[-1.7121083164595972, -2.5655428095216295]
[-0.5674027525415846, -2.992526926274017]
[0.8313656167997135, -2.5193347256443763]
[2.085945193325355, -2.0573707432814796]
[2.5869306117994584, -1.2549726253197953]
[2.9094031573517074, -0.3070263535954405]
[2.9684509920968605, 0.5953353058671058]
[2.250475430291036, 2.351070953959735]
[1.569332645014958, 2.676139601555899]
[-0.23562340136345455, 3.37989598199693]
[-3.048159617421087, 1.6798228422033785]
Higher-dimensional convex hull
One can compute the convex hull of points in higher dimensions using convex_hull
. The appropriate algorithm is decided based on the dimensionality of the given points.
using Polyhedra
v = [randn(3) for _ in 1:30]
hull = convex_hull(v)
typeof(hull), length(v), length(hull)
(Vector{Vector{Float64}}, 30, 15)
Here, convex_hull
is now using the concrete polyhedra library Polyhedra, hence it needs to be loaded beforehand.
One can check whether a point belongs to the convex hull using ∈
as follows:
P = VPolytope(hull)
x = sum(hull)/length(hull)
x ∈ P
true
Here x ∈ P
solves a feasibility problem; see the docs of ?∈
for details. Equivalently, using set inclusion:
Singleton(x) ⊆ P
true
If no additional arguments are passed, convex_hull
uses the default polyhedra library from default_polyhedra_backend
for the given input; different options can be passed through the backend
keyword; see the Julia polyhedra website for all the available backends.