# Mesh/grid over convex polytope
**URL:** https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047
**Category:** Numerics
**Tags:** question, polygons, convex-hull
**Created:** [March 16, 2020, 4:30pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047 "2020-03-16T16:30:24Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 16, 2020, 4:30pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/1 "2020-03-16T16:30:24Z")
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
1. make a _more or less equispaced_ (adjusted when necessary) grid \{ z\_j \}\_j over this polytope,
2. 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,
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](https://github.com/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.)
---
### Author:  [@jlchan](https://discourse.julialang.org/u/jlchan)
#### Post date: [March 16, 2020, 4:40pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/2 "2020-03-16T16:40:38Z")
For clarification, are the points z\_j meant to be only on the boundary or distributed throughout the interior of the polytope?
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 16, 2020, 5:05pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/3 "2020-03-16T17:05:41Z")
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.
Illustration (just the concept):
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [March 16, 2020, 5:21pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/4 "2020-03-16T17:21:39Z")
Why do you need it? For interpolation, or integration?
How about tiling the convex hull with tetrahedra (Delaunay) and then distributing the job across these tetrahedra?
---
### Author:  [@nboyd](https://discourse.julialang.org/u/nboyd)
#### Post date: [March 16, 2020, 6:03pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/5 "2020-03-16T18:03:12Z")
If speed isn’t important one possibility (almost certainly not the best idea…) for 2 is a small LP:
\min\_{w \ge 0, \sum w = 1, Zw = y} c^Tw
where c could be, for instance, c\_i = \| z\_i - y \|^2.
```julia
using Convex, SCS
using LinearAlgebra
d = 2
n = 100
Z = randn(d,n)
function interpolation_weights(points, point)
(d,n) = size(points)
w = Variable(n)
cost = map(p -> norm(p - point)^2, eachcol(points))
problem = minimize(dot(cost, w), [w >= 0, sum(w)==1, points*w == point])
solve!(problem, SCS.Optimizer)
vec(w.value)
end
y = randn(2)
@time w = interpolation_weights(Z, y)
using Makie
scatter(map(Point2, eachcol(Z)))
scatter!([Point2(y)], marker=:x, color=:red)
scatter!(map(Point2, eachcol(Z)), color=(:blue,0.5), markersize=w, transparency=true)
scatter!(map(Point2, eachcol(Z[:,findall(w .> 1E-5)])), color=:red)
```
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 17, 2020, 8:10am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/6 "2020-03-17T08:10:54Z")
> [@PetrKryslUCSD](#):
>
> Why do you need it? For interpolation, or integration?
For simplifying a simulation, as described [here](https://discourse.julialang.org/t/simplifying-simulation-by-clusters/35870).
> [@PetrKryslUCSD](#):
>
> How about tiling the convex hull with tetrahedra (Delaunay) and then distributing the job across these tetrahedra?
Tetrahedra could be a good idea, but maybe V-D would not be a good match: I want to keep it as regular as feasible.
Here is what I came up with (2D illustration, but remember the problem is 3D):
1. Find the facets [F, blue]
2. 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)
3. Partition the convex hull into simplexes (tetrahedra), [red lines]
4. 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).
---
### Author:  [@leethargo](https://discourse.julialang.org/u/leethargo)
#### Post date: [March 17, 2020, 8:23am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/7 "2020-03-17T08:23:56Z")
> [@Tamas\_Papp](#):
>
> 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(n2)
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](https://github.com/JuliaPolyhedra/CDDLib.jl).
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 17, 2020, 8:39am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/8 "2020-03-17T08:39:10Z")
> [@leethargo](#):
>
> I guess you need the combinatorial structure of the faces, that is, sets of vertices that span a face?
Yes, ideally organized into a graph by adjacency of vertices, but that can be figured out.
> [@leethargo](#):
>
> very high complexity (doubly exponential?)
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.
---
### Author:  [@leethargo](https://discourse.julialang.org/u/leethargo)
#### Post date: [March 17, 2020, 8:50am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/9 "2020-03-17T08:50:37Z")
Naive algorithm:
1. Compute outer description from points (with CDDLib.jl) → list of inequalities (one for each facet).
2. For each inequality, evaluate on every point. The points that satisfy with equality (no slack) are part of that facet.
I guess the tricky part is defining the tolerance for evaluating the inequality, if you don’t with integer coefficients or similar.
---
### Author:  [@leethargo](https://discourse.julialang.org/u/leethargo)
#### Post date: [March 17, 2020, 8:53am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/10 "2020-03-17T08:53:16Z")
I just saw that `Polyhedra.jl` already defines some [methods for incidence](https://juliapolyhedra.github.io/Polyhedra.jl/stable/polyhedron/#Incidence-1), but I’m not familiar with them.
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [March 17, 2020, 11:08am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/11 "2020-03-17T11:08:53Z")
> [@Tamas\_Papp](#):
>
> I am totally fine with whatever approximation/definition of “equispaced” and "nearest’ that make this problem easy.
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](https://github.com/JuliaReach/LazySets.jl/). I used 2D to simplify, but it can be generalized straightforwardly.
```julia
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:
```julia
V = rand(VPolygon, num_vertices=5)
plot(V, color=:orange)
U = UnionSetArray(mesh(V))
plot!(U, color=:red)
```
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [March 17, 2020, 11:39am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/12 "2020-03-17T11:39:35Z")
The version below is for arbitrary dimension (`n >= 1` for any `n`). Now `V` is annotated as an [abstract polytpe](https://juliareach.github.io/LazySets.jl/dev/lib/interfaces/#LazySets.AbstractPolytope) 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.
```julia
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](https://juliareach.github.io/LazySets.jl/dev/lib/sets/VPolygon/) (resp. [VPolytope](https://juliareach.github.io/LazySets.jl/dev/lib/sets/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)`.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 17, 2020, 12:16pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/13 "2020-03-17T12:16:24Z")
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.
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [March 17, 2020, 12:27pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/14 "2020-03-17T12:27:22Z")
> [@Tamas\_Papp](#):
>
> vertices should be there,
The vertices can be added with `push!(inner_points, vertices_list(V)...)`.
> [@Tamas\_Papp](#):
>
> and some points on the edges
If `V` is a `VPolygon` (ie. 2D), the 1dim “facets” or edges can be obtained like so
```julia
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.
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [March 17, 2020, 8:07pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/15 "2020-03-17T20:07:19Z")
Here is a proof of principle implementation of the idea in my last post. The code works for any dimension `n >= 2`.
```julia
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:
```julia
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)
```
Some other comments:
- 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.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 18, 2020, 8:21am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/16 "2020-03-18T08:21:09Z")
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).
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [March 18, 2020, 11:35am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/17 "2020-03-18T11:35:59Z")
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:
```julia
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:
```julia
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)
```
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 18, 2020, 11:52am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/18 "2020-03-18T11:52:59Z")
I just realized that my [naive algorithm](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/6) 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.
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [March 18, 2020, 2:05pm UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/19 "2020-03-18T14:05:25Z")
If convex, I would use one of the Delaunay triangulation packages from [https://github.com/JuliaPDE/SurveyofPDEPackages](https://github.com/JuliaPDE/SurveyofPDEPackages).
---
### Author:  [@mforets](https://discourse.julialang.org/u/mforets)
#### Post date: [March 19, 2020, 10:47am UTC](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/20 "2020-03-19T10:47:31Z")
Interesting! We’ve tagged a new LazySets release adding MiniQhull as an optional dependency, so one can triangulate with:
```julia
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
```
