# Voronoi tessellations on meshed surfaces

Given a meshed surface, on which we can calculate geodesic distances between vertices, how can one calculate the Voronoi tessellation of a set of points located on this surface?

This is somewhat related to the question here, although that question is restricted to a Voronoi tessellation on a unit sphere. See also here for other approaches.

For this answer, I've slightly streamlined Dunlop's code. As with his routines, the initialization and solving steps are separate; one particular wrinkle in mine is that I wrote special routines for solving the heat equation for the case of multiple points (represented as indices of the associated mesh's vertices), as well as for a single point. The multiple point solver is more efficient than mapping the single point solver across the multiple points.

heatMethodInitialize[mesh_MeshRegion] :=
gm1, gm2, gm3, gradOp, nlen, nrms, oped, polys, sa1, sa2, sa3,
tmp, wi1, wi2, wi3},

vertices = MeshCoordinates[mesh];
faces = First /@ MeshCells[mesh, 2];
polys = Map[vertices[[#]] &, faces];

edges = First /@ MeshCells[mesh, 1];

tmp = Transpose[polys, {1, 3, 2}];
nrms = MapThread[Dot, {ListConvolve[{{-1, 1}}, #, {{2, -1}}] & /@ tmp,
ListConvolve[{{1, 1}}, #, {{-2, 2}}] & /@ tmp}, 2];
nlen = Norm /@ nrms; nrms /= nlen;

oped = ListCorrelate[{{1}, {-1}}, #, {{3, 1}}] & /@ polys;

wi1 = MapThread[Cross, {nrms, oped[[All, 1]]}];
wi2 = MapThread[Cross, {nrms, oped[[All, 2]]}];
wi3 = MapThread[Cross, {nrms, oped[[All, 3]]}];

{MapIndexed[Transpose[PadLeft[{#}, {2, 3}, #2]] &, faces],
Transpose[{wi1[[All, 1]], wi2[[All, 1]], wi3[[All, 1]]}]},
2]]];
{MapIndexed[Transpose[PadLeft[{#}, {2, 3}, #2]] &, faces],
Transpose[{wi1[[All, 2]], wi2[[All, 2]], wi3[[All, 2]]}]},
2]]];
{MapIndexed[Transpose[PadLeft[{#}, {2, 3}, #2]] &, faces],
Transpose[{wi1[[All, 3]], wi2[[All, 3]], wi3[[All, 3]]}]},
2]]];
adi = SparseArray[Band[{1, 1}] -> 1/nlen];

gradOp = Transpose[SparseArray[{gm1, gm2, gm3}], {2, 1, 3}];

areas = PropertyValue[{mesh, 2}, MeshCellMeasure];
ada = SparseArray[Band[{1, 1}] -> 2 areas];
divMat = Transpose[#].ada & /@ {gm1, gm2, gm3};
del = divMat[[1]].gm1 + divMat[[2]].gm2 + divMat[[3]].gm3;

With[{spopt = SystemOptions["SparseArrayOptions"]},
InternalWithLocalSettings[
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}],

nlen /= 2;
acm = SparseArray[MapThread[{#1, #1} -> #2 &,
{Flatten[Transpose[faces]],
Flatten[ConstantArray[nlen, 3]]}]],
SetSystemOptions[spopt]]];

heatSolve[mesh_MeshRegion, acm_, del_,
gradOp_, divMat_][idx_Integer, t : (_?NumericQ | Automatic) : Automatic] :=
Module[{h, tm, u},
tm = If[t === Automatic,
Max[PropertyValue[{mesh, 1}, MeshCellMeasure]]^2, t];
u = LinearSolve[acm + tm del, UnitVector[MeshCellCount[mesh, 0], idx]];

heatSolve[mesh_MeshRegion, acm_, del_,
t : (_?NumericQ | Automatic) : Automatic] :=
Module[{h, tm, u},
tm = If[t === Automatic,
Max[PropertyValue[{mesh, 1}, MeshCellMeasure]]^2, t];
u = Transpose[LinearSolve[acm + tm del, Normal[SparseArray[
MapIndexed[Prepend[#2, #1] &, idx] -> 1,
{MeshCellCount[mesh, 0], Length[idx]}]]]];
h = Transpose[-Normalize /@ Normal[gradOp.#]] & /@ u;
h = Transpose[Total[MapThread[Dot, {divMat, #}]] & /@ h];
LinearSolve[del, h]]


With these routines, here's how to generate a(n approximate) Voronoi diagram on the Stanford bunny:

bunny = ExampleData[{"Geometry3D", "StanfordBunny"}, "MeshRegion"];
vertices = MeshCoordinates[bunny]; faces = First /@ MeshCells[bunny, 2];

npoints = 9;
randvertlist = BlockRandom[SeedRandom[42, Method -> "Legacy"]; (* for reproducibility *)
RandomSample[Range[MeshCellCount[bunny, 0]], npoints]];

{am, Δ, gr, dv} = Most @ heatMethodInitialize[bunny];
Φ = heatSolve[bunny, am, Δ, gr, dv][randvertlist, 0.5];

cols = Table[ColorData[61] @ Ordering[v, 1][[1]], {v, Φ}];

Graphics3D[{{Green, Sphere[vertices[[randvertlist]], 0.003]},
GraphicsComplex[vertices, {EdgeForm[], Polygon[faces]},
VertexColors -> cols]},
Boxed -> False, Lighting -> "Neutral"]


• Looks great, will have to try it out when I get a bit more time. I have been struggling the last weeks to get the boundaries nicely highlighted and hope to get that implemented soon. Great work! Is it possible to win much speed by compiling, or is this problematic as one loses the use of the built in mesh functions that would have to be rewritten? Commented Apr 5, 2017 at 20:18
• The SparseArray[] stuff definitely isn't compilable, and in this case, they're the least memory-intensive storage format. The fuzziness in this case is expected, since at the moment, one is only classifying mesh vertices according to which feature point is nearest. For a sharp boundary, one would need to be able to draw across mesh polygons that might get split. Commented Apr 5, 2017 at 22:11
• That is exactly what is beating me, trying to split the polygons to get a sharper mesh. On smaller meshes I managed but larger ones (like the bunny) slowed things down enormously that one loses all the advantage of the speed of this algorithm. Commented Apr 6, 2017 at 3:10

Using the Geodesics in Heat Algorithm implemented here, we can calculate the distances of all vertices on the surface to a given vertex. By repeating this algorithm on a selected subset of vertices on the surface, we can calculate readily how close all the other vertices on these surfaces are, and label them according to which selected vertex they are closest to.

This can be done with the following code, with the Stanford bunny example again:

a = DiscretizeGraphics[ExampleData[{"Geometry3D", "StanfordBunny"}]]
prep = heatDistprep[a];
npoints = 10;
nvertices = prep[[5]];
vertices = prep[[6]];
faces = MeshCells[a, 2] /. Polygon[p_] :> p;
phiall = {};
randvertlist =
DeleteDuplicates[RandomInteger[{1, nvertices}, npoints]];
npoints = Length[randvertlist];
i = 1;
While[i < npoints + 1,
phi = solveHeat[a, prep, randvertlist[[i]], 0.5];
AppendTo[phiall, phi[[1]]];
i++];
labels = Map[Ordering[phiall[[All, #]]][[1]] &, Range[nvertices]]/npoints;
plotdata =
Map[Join[vertices[[#]], {labels[[#]]}] &, Range[Length[vertices]]];
labelplot =
Graphics3D[{EdgeForm[],
GraphicsComplex[vertices, Map[Polygon, faces],
VertexColors ->
Table[ColorData["BrightBands"][labels[[i]]], {i, 1,
nvertices}]]}, Boxed -> False, Lighting -> "Neutral"];
pointplot =
Graphics3D[{Black, Ball[Map[vertices[[#]] &, randvertlist], 0.003]},
Boxed -> False];
Show[{pointplot, labelplot}]


An issue with this approach is that the boundaries in the visualisation are somewhat "rough", and we don't get directly the edges of the Voronoi cells. Any hints on how to do this would be great.

I hope someone finds this useful.

To make the answer self-contained, the Geodesics code is given here:

heatDistprep[mesh0_] := Module[{a = mesh0, vertices, nvertices, edges, edgelengths, nedges, faces, faceareas, unnormfacenormals, acalc, facesnormals, facecenters, nfaces, oppedgevect, wi1, wi2, wi3, sumAr1, sumAr2, sumAr3, areaar, gradmat1, gradmat2, gradmat3, gradOp, arear2, divMat, divOp, Delta, t1, t2, t3, t4, t5, , Ac, ct, wc, deltacot, vertexcoordtrips, adjMat},
vertices = MeshCoordinates[a]; (*List of vertices*)
edges = MeshCells[a, 1] /. Line[p_] :> p; (*List of edges*)
faces = MeshCells[a, 2] /. Polygon[p_] :> p; (*List of faces*)
nvertices = Length[vertices];
nedges = Length[edges];
nfaces = Length[faces];
adjMat = SparseArray[Join[({#1, #2} -> 1) & @@@ edges, ({#2, #1} -> 1) & @@@edges]]; (*Adjacency Matrix for vertices*)
edgelengths = PropertyValue[{a, 1}, MeshCellMeasure];
faceareas = PropertyValue[{a, 2}, MeshCellMeasure];
vertexcoordtrips = Map[vertices[[#]] &, faces];
unnormfacenormals = Cross[#3 - #2, #1 - #2] & @@@ vertexcoordtrips;
acalc = (Norm /@ unnormfacenormals)/2;
facesnormals = Normalize /@ unnormfacenormals;
facecenters = Total[{#1, #2, #3}]/3 & @@@ vertexcoordtrips;
oppedgevect = (#1 - #2) & @@@ Partition[#, 2, 1, 3] & /@vertexcoordtrips;
wi1 = -Cross[oppedgevect[[#, 1]], facesnormals[[#]]] & /@Range[nfaces];
wi2 = -Cross[oppedgevect[[#, 2]], facesnormals[[#]]] & /@Range[nfaces];
wi3 = -Cross[oppedgevect[[#, 3]], facesnormals[[#]]] & /@Range[nfaces];
sumAr1 = SparseArray[Join[Map[{#, faces[[#, 1]]} -> wi1[[#, 1]] &, Range[nfaces]],Map[{#, faces[[#, 2]]} -> wi2[[#, 1]] &, Range[nfaces]],Map[{#, faces[[#, 3]]} -> wi3[[#, 1]] &, Range[nfaces]]]];
sumAr2 = SparseArray[Join[Map[{#, faces[[#, 1]]} -> wi1[[#, 2]] &, Range[nfaces]], Map[{#, faces[[#, 2]]} -> wi2[[#, 2]] &, Range[nfaces]],Map[{#, faces[[#, 3]]} -> wi3[[#, 2]] &, Range[nfaces]]]];
sumAr3 =SparseArray[Join[Map[{#, faces[[#, 1]]} -> wi1[[#, 3]] &, Range[nfaces]], Map[{#, faces[[#, 2]]} -> wi2[[#, 3]] &, Range[nfaces]], Map[{#, faces[[#, 3]]} -> wi3[[#, 3]] &, Range[nfaces]]]];
areaar = SparseArray[Table[{i, i} -> 1/(2*acalc[[i]]), {i, nfaces}]];
arear2 = SparseArray[Table[{i, i} -> (2*faceareas[[i]]), {i, nfaces}]];
divOp[q_] := divMat[[1]].q[[All, 1]] + divMat[[2]].q[[All, 2]] + divMat[[3]].q[[All, 3]];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}]; (*Required to allow addition of value assignment to Sparse Array*)
t1 = Join[faces[[All, 1]], faces[[All, 2]], faces[[All, 3]]];
t2 = Join[acalc, acalc, acalc];
Ac = SparseArray[Table[{t1[[i]], t1[[i]]} -> t2[[i]], {i, nfaces*3}]];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
]

solveHeat[mesh0_, prepvals_, i0_, t0_] := Module[{nvertices, delta, t, u, Ac, Delta, g, h, phi, gradOp, divOp, vertices, plotdata},
vertices = prepvals[[6]];
nvertices = prepvals[[5]];
Ac = prepvals[[1]];
Delta = prepvals[[2]];
`