14
$\begingroup$

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.

$\endgroup$
9
$\begingroup$

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"]

Voronoi bunny

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Dunlop Apr 5 '17 at 20:18
  • 1
    $\begingroup$ 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. $\endgroup$ – J. M. is away Apr 5 '17 at 22:11
  • $\begingroup$ 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. $\endgroup$ – Dunlop Apr 6 '17 at 3:10
17
$\begingroup$

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}]

Bunny

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}
] 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.