Voronoi tessellations on meshed surfaces
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}
]
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"]
