Skip to main content

computational geometry - 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.



Answer



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


Comments

Post a Comment

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more