Voronoi mesh from non-parametric surfaces for 3D printing

Question

I am trying to make a Voronoi mesh from the Stanford Bunny surface for 3D printing, intended as a test case to make a Voronoi mesh from any arbitrary, non-parametric 3D surface. Related techniques have been published for exoskeletons and ornamental curve networks and is seen in the functionality of Meshmixer.

Several posted scripts come close, but so far I have been unable to produce a complete solution from start to finish. Question 141348 goes far, but I do not know how to convert those colored Voronoi surfaces to a useful Voronoi mesh. Mathematica has a similar demo for parametric surfaces.

Question 135094 detects edges on Image3D objects; question 3327 converts Graphics3D to Image3D. However, hooking these scripts together appears to leave a disconnect for non-parametric surfaces (and reduces to two dimensions, but I think that's a bug on my part).

Question 39879 has useful tips for converting a tiled surface on a parametric surface to a wire frame and then to a tube frame for printing, but it is not clear how to apply to an arbitrary surface.

As you can see, Mathematica is a little foreign to me--surely there is a more coherent way to do this than the mixed blocks of code below taken from the references above?

I'd appreciate any help or insight. Thanks!

(*From "Voronoi tessellations on meshed surfaces," 141348*)
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}]

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}]
    ]
   },
  Axes -> True,
  Boxed -> True,
  Lighting -> "Neutral"]

(* From "How to convert Graphics3D object into an Image3D object?", 33274*)
obj = labelplot;
slice[obj_, x_, dx_] := 
  Show[obj, ViewPoint -> {\[Infinity], 0, 0}, 
   PlotRange -> {{x, x + dx}, All, All}, Axes -> False, 
   Boxed -> False];
slice[obj, -0.1, 0.2]
frames = Table[
   ImageData@Thinning@ColorNegate@ColorConvert[#, "Grayscale"] &@
    Rasterize@slice[obj, x, 0.01], {x, -0.1, 0.1, 0.01}];
Image3D[frames]

Answer 1

Disclaimer: This is currently an imperfect solution.

Implementation

I use the optimized implementation of the heat method from 175570 to approximate geodesic distance functions. This is also the most time consuming part.

(* precompute data for heat method *)

R = ExampleData[{"Geometry3D", "StanfordBunny"}, "MeshRegion"];
data = heatDistprep2[R, 0.001];
(* pick some points from the triangulation as centers of the Voronoi cells *)
idx = RandomSample[Range[1, MeshCellCount[R, 0]], 400];
(* compute distance functions *)
ϕ = Transpose[Map[solveHeat2[R, data, #] &, idx]]; // AbsoluteTiming // First

14.9839

Moreover, I require the function IGMeshCellAdjacencyMatrix from Szabolcs' IGraphM package.

Needs["IGraphM`"]

nearest = Flatten[Ordering[#, 1] & /@ ϕ];
mindist = Min /@ ϕ;

pts = MeshCoordinates[R];

edges = MeshCells[R, 1, "Multicells" -> True][[1, 1]];

(* For each edge, determine whether its endpoints lie in different Voronoi cells. *)

edgenearestlists = Partition[nearest[[Flatten[edges]]], 2];
cutedgelist = 
  Flatten[SparseArray[
     Flatten[Differences@Transpose[edgenearestlists]]][
    "NonzeroPositions"]];
cutedges = edges[[cutedgelist]];
edgecutnumbers = 
  Normal[SparseArray[Partition[cutedgelist, 1] -> 1, 
    MeshCellCount[R, 1]]];

(* For each edge with endpoints in different Voronoi cells, determine the intersection point of the edge and the boundary of these cells. *)

(* Warning: Some (wrong) heuristic assumption go into this. *)

{i1, i2} = Transpose[cutedges];
Module[{t},
  f = {u1, u2, v1, v2} \[Function] Evaluate[
     t /. First@Solve[
        u1 (1 - t) + t u2 == v1 (1 - t) + t v2,
        t
        ]
     ]
  ];
t = f @@ Flatten[Outer[Extract[ϕ, Transpose[{#2, #1}]] &,
     {nearest[[i1]], nearest[[i2]]},
     {i1, i2},
     1], 1];
edgecutpts = pts[[i1]] (1 - t) + t pts[[i2]];

(* For each triangle, count the number of intersection points of its  edges with boundaries of Voronoi cells. *)

A = IGMeshCellAdjacencyMatrix[R, 2, 1];
{i1, i2, i3} = Transpose[MeshCells[R, 2, "Multicells" -> True][[1, 1]]];
c1 = nearest[[i1]];
c2 = nearest[[i2]];
c3 = nearest[[i3]];
trianglecutnumbers = Unitize[Transpose[{c2 - c1, c3 - c2, c1 - c3}]].ConstantArray[1, 3];

(* For each triangle that is intersected exactly twice, make an edge for the new line complex. *)

edges2 = A[[Flatten[Position[trianglecutnumbers, 2]], cutedgelist]][
   "AdjacencyLists"];

(* For each triangle that is intersected exactly thrice, find the intersection point p of Voronoi cell boundaries.  *)

cuttrianglelist = Flatten[Position[trianglecutnumbers, 3]];
Module[{x, y},
  g = {u1, u2, u3, v1, v2, v3, w1, w2, w3} \[Function]
    Evaluate[
     {x, y} /. First@Solve[{
         u1 + x (u2 - u1) + y (u3 - u1) == 
          v1 + x (v2 - v1) + y (v3 - v1),
         u1 + x (u2 - u1) + y (u3 - u1) == 
          w1 + x (w2 - w1) + y (w3 - w1)
         },
        {x, y}
        ]
     ]
  ];

{i1, i2, i3} = Transpose[MeshCells[R, 2, "Multicells" -> True][[1, 1]][[cuttrianglelist]]];
{x, y} = g @@ Flatten[Outer[Extract[ϕ, Transpose[{#2, #1}]] &,
     {nearest[[i1]], nearest[[i2]], nearest[[i3]]},
     {i1, i2, i3},
     1], 1];
trianglecutpts = pts[[i1]] (1 - x - y) + x pts[[i2]] + y pts[[i3]];

(* For each triangle that is intersected exactly thrice, and create \
three edges for the new line complex, each emanating from an \
intersection point on an edge and ending in p.  *)

edges3 = With[{cut3 = A[[cuttrianglelist, cutedgelist]]["AdjacencyLists"]},
   Transpose[{
     Flatten[Transpose[
       ConstantArray[
        Range[Length[edgecutpts] + 1, 
         Length[edgecutpts] + Length[cut3]], 3]
       ]],
     Flatten[cut3]
     }
    ]
   ];

The result looks like this:

θ = 0.0005;
Show[
 (*HighlightMesh[R,Thread[{2,cuttrianglelist}]],*)
 (*plot[R,mindist],*)
 R,
 Graphics3D[{
   Darker@Darker@Green, Specularity[White, 30], 
   Sphere[#, 2 θ] & /@ pts[[idx]],
   GraphicsComplex[
    Join[edgecutpts, trianglecutpts],
    {Gray, Specularity[White, 30], 
     Tube[Join[edges2, edges3], θ]}
    ],
   Darker@Darker@Red, Specularity[White, 30], 
   Sphere[#, 2 θ] & /@ trianglecutpts
   }],
 Lighting -> "Neutral"
 ]

The MeshRegion representing this line complex is

MeshRegion[
 Join[edgecutpts, trianglecutpts], 
 Line[Join[edges2, edges3]]
 ]

Remark

This is far from perfect. Right after I implemented this, I realized that I made the false assumption that each edge of the original triangle complex can intersect at most one boundary edge of a Voronoi cell. So my strategy to compute the points where three Voronoi cells meet (the red point in the plot; let's call them triple points) does not always produce correct results. Accordingly, also the edges of the final line complex that contain triple points are incorrect.

If the triangles of the original complex are sufficiently small compared to the Voronoi cells, this might be negligible, in particular if the line complex has to be thickened for printing, anyways.

Also, please notice that the Mathematica's bunny mesh is not an intrinsic Delaunay mesh. This may cause problems for computing the distance function with the heat method; for example, negative distances my arrise. So one should apply a remeshing step before applying the heat method.