I have x_1, \dots, x_m points in \mathbb{R}^n defining a polytope as their convex hull, m > n. Typically n = 2 or 3 and m = 4 or 5.
I would like to
make a more or less equispaced (adjusted when necessary) grid \{ z_j \}_j over this polytope,
for an arbitrary point y inside the convex hull, compute positive weights \sum_i w_i = 1 such that for some set of “nearest” points z_1, \dots, z_k of the grid,
y = \sum_i w_i z_i
I am totally fine with whatever approximation/definition of “equispaced” and "nearest’ that make this problem easy. I am solving this problem a few times for a very long computation, on about 50-200 gridpoints, so it does not have to be super-fast.
I have been looking at packages in JuliaGeometry but this is outside my area of expertise, so I don’t know where to start. Maybe FEM uses methods like this?
(If my domain was a generalized hypercube, I would just cut it up to little hypercubes and use the tensored interpolating weights. But unfortunately it isn’t.)
I need some one the boundary to be able to approximate all interior points, but I would prefer some in the interior too so that I always have something close to all interior points.
Using the “center of mass” [C, red], (could be weighted, all that matters is that it is a nice interior point that gives more or less balanced partition below)
Partition the convex hull into simplexes (tetrahedra), [red lines]
Cut each simplex into self-similar pieces, effectively getting d^n subcells.
This would give me a very fast and simple way to do everything using barycentric coordinates once I determine the relevant simplex.
I can implement most of this, except finding the facets of the convex hull. Help with this would be appreciated — I don’t need the most efficient algorithm, in fact I am happy to brute force this with something O(n^2). I am just hoping that building blocks exist within Julia so that I don’t have to figure this out from first principles or learn the relevant geometry (again, it is fascinating, but not my field).
I guess you need the combinatorial structure of the faces, that is, sets of vertices that span a face?
Not sure how to compute that, but going from the inner description (convex hull of points) to the outer description (intersection of halfspaces, defined by linear inequalities) has a very high complexity (doubly exponential?). See, e.g., CDDLib.jl.
Yes, ideally organized into a graph by adjacency of vertices, but that can be figured out.
Again, this is not a concern for me as I have 4-6 vertices. Basically anything goes, the simpler the better.
Naively, I could even form all the \binom{m}{3} vertices and eliminate the ones I don’t need (just don’t know how). Or form all the \binom{m}{2} edge vectors.
I am just not used to doing these things, so if someone could help me with the algorithm that is so embarrassingly naive and costly that it is only taught to people in the field so that they can get a good laugh about how people thought about this in the stone age/ancient Greece/the middle ages, even that could help.
Since the dimension (n) is low, I think that rejection sampling should work pretty well for subproblem 1: start with an equally-spaced grid over the bounding box of the polytope, then do membership tests to reject those points that are outside. The following picture is an example in 2D:
Below is an implementation using LazySets.jl. I used 2D to simplify, but it can be generalized straightforwardly.
using LazySets, Plots
using LazySets: center
function mesh(V::VPolygon{N, VN}, Δ=0.2) where {N, VN<:AbstractVector{N}}
B = overapproximate(V, BallInf) # bounding box
c = center(B)
r = radius(B)
# equally spaced mesh
mesh_x = range(c[1] - r, c[1] + r, step=Δ)
mesh_y = range(c[2] - r, c[2] + r, step=Δ)
# set union of singletons
U = UnionSetArray([Singleton([mx, my]) for mx in mesh_x for my in mesh_y])
inner_points = Vector{Singleton{N, VN}}()
for s in array(U)
p = element(s)
if p ∈ V
push!(inner_points, s)
end
end
return inner_points
end
Example:
V = rand(VPolygon, num_vertices=5)
plot(V, color=:orange)
U = UnionSetArray(mesh(V))
plot!(U, color=:red)
The version below is for arbitrary dimension (n >= 1 for any n). Now V is annotated as an abstract polytpe so it can use VPolytope as a special case. I also changed the bounding box to be a bounding hyperrectangle (Hyperrectangle) so that the outer box is tight in every dimension.
function mesh(V::AbstractPolytope{N}, Δ=0.2) where {N}
n = dim(V)
B = overapproximate(V, Hyperrectangle)
c = center(B)
r = radius_hyperrectangle(B)
# equally spaced mesh
aux = [range(c[i] - r[i], c[i] + r[i], step=Δ) for i in 1:n]
mesh_points = Base.Iterators.product(aux...)
vlist = vertices_list(V)
VN = eltype(vlist)
inner_points = Vector{VN}()
for tup in mesh_points
p = collect(tup)
if p ∈ V
push!(inner_points, p)
end
end
return inner_points
end
For plotting, you have to use the plot recipe of a LazySet, to plot U = UnionSetArray([Singleton(p) for p in mesh(V)]). Here U represents the set union of singleton sets (points).
(For 3D plotting i have used Makie in the past, and we have a plot3d plot recipe in LazySets, but i haven’t used it since a long time for compilation issues with Makie’s dependencies on my machine.)
One note about the input sets: in LazySets, the VPolygon (resp. VPolytope) types have constructor flags apply_convex_hull which default to true, meaning that if you pass VPolytope(vertices) it will run a convex hull unless you pass apply_convex_hull=false to the constructor. On the other hand, the (concrete) convex convex hull of a set of points can be computed with convex_hull(points).
This is neat, but I need to be able to approximate all possible points in the polytope as convex combinations, so eg vertices should be there, and some points on the edges.
The vertices can be added with push!(inner_points, vertices_list(V)...).
If V is a VPolygon (ie. 2D), the 1dim “facets” or edges can be obtained like so
v = vertices_list(V)
m = length(v)
segments = Vector{LineSegment{Float64}}()
push!(segments, LineSegment(v[m], v[1]))
for i in 1:m-1
l = LineSegment(v[i], v[i+1])
push!(segments, l)
end
(see with plot([segments[i] for i in 1:length(segments)], color=:black)). Then, for each line segment you can sample it as desired. The implementation above works because in LazySets, VPolygon have their vertices sorted in counter-clockwise fashion.
In higher dimension if you want to iterate over the facets then this is not so immediate, what i would try as a workaround is to intersect each defining half-space with the bounding box, then use the n-1-dimensional sampling idea from my previous post.
Here is a proof of principle implementation of the idea in my last post. The code works for any dimension n >= 2.
function _convex_combination(vertices)
m = length(vertices)
weights = normalize(abs.(rand(m)), 1)
return sum(weights .* vertices)
end
# for each n-1 dimensional facet we store npoints random points
function mesh_facets(V::AbstractPolytope{N}, npoints=10) where {N}
# compute bounding box
B = overapproximate(V, Hyperrectangle)
# preallocate points in the "border"
vertices = vertices_list(V)
VN = eltype(vertices)
border_points = Vector{VN}()
# add vertices
push!(border_points, vertices...)
# add samples for each n-1 dimensional facet
H = convert(HPolytope, V)
constraints = constraints_list(H)
facets = [intersection(Hyperplane(c.a, c.b), H) for c in constraints]
# random sampling on each facet
for F in facets
for k in 1:npoints
p = _convex_combination(vertices_list(F))
push!(border_points, p)
end
end
return border_points
end
Example:
V = rand(VPolygon, num_vertices=5)
points_inner = mesh(V)
points_border = mesh_facets(V)
plot(V)
plot!(UnionSetArray([Singleton(p) for p in mesh_inner]), color=:red)
plot!(UnionSetArray([Singleton(p) for p in mesh_border]), color=:blue)
This code can be optimized rather easily for 2D, just try to use HPolygon/VPolygon and dispatch will take care.
I used the default polyhedra backend but for efficiency you can try CDDLib, just pas it as backend to the appropriate conversion/intersection functions, or ask back if you have doubts.
Note that “border” here is only the n-1 dimensional facets (and of course the 0 dimensional facets which are the vertices). So in 3d, there’ll be low probability to hit points on the edges (but there will be points on the faces of the polytope). I think that the same idea can be extrapolated to include n-2 dimensional facets and so on.
Thanks for investing in a solution. What I do not yet see how to do with this scheme is construct the interpolations, for which I would need to find the nearest points (see the original question).
For subproblem 2 as suggested above by @nboyd one possibility is to use linear programming. An alternative to that code is to pre-select the z_1, ..., z_k points (either take all m mesh points or keep k < m closest to y for a given norm) and then solve a feasibility LP for such k points to find the weights.
Solution:
function preselect(points::Vector{VN}, y::VN, tol=0.2, pnorm=Inf) where {N, VN<:AbstractVector{N}}
closest_points = Vector{VN}()
point_ids = Vector{Int}()
for (i, p) in enumerate(points)
δ = norm(p - y, pnorm)
if δ < tol
push!(closest_points, p)
push!(point_ids, i)
end
end
return closest_points, point_ids
end
using JuMP, GLPK
function interpolation_weights(points, y)
model = Model(with_optimizer(GLPK.Optimizer))
m = length(points)
@variable(model, w[1:m] >= 0)
@constraint(model, sum(w) == 1)
Z = hcat(points...)
@constraint(model, Z * w .== y)
optimize!(model)
return value.(w)
end
Usage:
V = rand(VPolygon, num_vertices=5)
mesh_inner = mesh(V)
mesh_border = mesh_facets(V);
mesh_points = vcat(mesh_inner, mesh_border);
plot(V)
plot!(UnionSetArray([Singleton(p) for p in mesh_points]), color=:red)
v = vertices_list(V)
y = sum(v) / length(v) # "test" point
plot!(Singleton(y), color=:yellow)
zk, ids = preselect(mesh_points, y);
weights = interpolation_weights(zk, y)
@assert sum(weights) == 1
# check that we get y
plot!(UnionSetArray([Singleton(p) for p in zk]), color=:cyan, marker=:x)
plot!(Singleton(sum(zk .* weights)), marker=:x)
I just realized that my naive algorithm for partitioning won’t work because the faces may not be triangles.
Is there a general algorithm for partitioning a polytope, specified by vertices, into simplexes? Possibly available in Julia, or easy to implement? Again, does not need to be efficient.
Interesting! We’ve tagged a new LazySets release adding MiniQhull as an optional dependency, so one can triangulate with:
using LazySets, MiniQhull
V = rand(VPolygon, num_vertices=5)
plot(delaunay(V)) # here delaunay(V) returns the set union of VPolytopes
# + rejection-sampled based grid from previous comment
