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

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

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more