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] :=
Module[{acm, ada, adi, adjMat, areas, del, divMat, edges, faces, vertices,
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];
adjMat = AdjacencyMatrix[UndirectedEdge @@@ edges];
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]]}];
sa1 = SparseArray[Flatten[MapThread[Rule,
{MapIndexed[Transpose[PadLeft[{#}, {2, 3}, #2]] &, faces],
Transpose[{wi1[[All, 1]], wi2[[All, 1]], wi3[[All, 1]]}]},
2]]];
sa2 = SparseArray[Flatten[MapThread[Rule,
{MapIndexed[Transpose[PadLeft[{#}, {2, 3}, #2]] &, faces],
Transpose[{wi1[[All, 2]], wi2[[All, 2]], wi3[[All, 2]]}]},
2]]];
sa3 = SparseArray[Flatten[MapThread[Rule,
{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];
gm1 = adi.sa1; gm2 = adi.sa2; gm3 = adi.sa3;
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"]},
Internal`WithLocalSettings[
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}],
nlen /= 2;
acm = SparseArray[MapThread[{#1, #1} -> #2 &,
{Flatten[Transpose[faces]],
Flatten[ConstantArray[nlen, 3]]}]],
SetSystemOptions[spopt]]];
{acm, del, gradOp, divMat, adjMat}]
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]];
h = Transpose[-Normalize /@ Normal[gradOp.u]];
LinearSolve[del, Total[MapThread[Dot, {divMat, h}]]]]
heatSolve[mesh_MeshRegion, acm_, del_,
gradOp_, divMat_][idx_ /; VectorQ[idx, IntegerQ],
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"]
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.
$\endgroup$
Commented
Apr 5, 2017 at 22:11
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}]];
gradmat1 = areaar.sumAr1;
gradmat2 = areaar.sumAr2;
gradmat3 = areaar.sumAr3;
gradOp[u_] := Transpose[{gradmat1.u, gradmat2.u, gradmat3.u}];
arear2 = SparseArray[Table[{i, i} -> (2*faceareas[[i]]), {i, nfaces}]];
divMat = {Transpose[gradmat1].arear2, Transpose[gradmat2].arear2,Transpose[gradmat3].arear2};
divOp[q_] := divMat[[1]].q[[All, 1]] + divMat[[2]].q[[All, 2]] + divMat[[3]].q[[All, 3]];
Delta = divMat[[1]].gradmat1 + divMat[[2]].gradmat2 + divMat[[3]].gradmat3;
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}];
{Ac, Delta, gradOp, divOp, nvertices, vertices, adjMat}
]
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]];
gradOp = prepvals[[3]];
divOp = prepvals[[4]];
delta = Table[If[i == i0, 1, 0], {i, nvertices}];
t = t0;
u = LinearSolve[(Ac + t*Delta), delta];
g = gradOp[u];
h = -Normalize /@ g;
phi = LinearSolve[Delta, divOp[h]];
plotdata = Map[Join[vertices[[#]], {phi[[#]]}] &, Range[Length[vertices]]];
{phi}
]
