# Voronoi mesh from non-parametric surfaces for 3D printing

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,
divMat, divOp, Delta, t1, t2, t3, t4, t5, , Ac, ct, wc, deltacot,
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];
SparseArray[
Join[({#1, #2} -> 1) & @@@ edges, ({#2, #1} -> 1) & @@@
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}]];
arear2 =
SparseArray[Table[{i, i} -> (2*faceareas[[i]]), {i, nfaces}]];
divOp[q_] :=
divMat[[1]].q[[All, 1]] + divMat[[2]].q[[All, 2]] +
divMat[[3]].q[[All, 3]];
Delta =
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}];

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]];
divOp = prepvals[[4]];
delta = Table[If[i == i0, 1, 0], {i, nvertices}];
t = t0;
u = LinearSolve[(Ac + t*Delta), delta];
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]

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. *)

{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]][

(* 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[
(*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.

• So impressive and so quick--thanks for the reply. It's not quite running on my instance quite yet, so I'm curious about what version of Mathematica you are using and if any other packages are required. Every instance of MeshCells[R, 2, "Multicells" -> True] is shadowed (warning in a red font), and it produces a stream of errors, starting with those below. Thanks again for the reply. Transpose: The first two levels of {1, 2, 3} cannot be transposed. Set: Lists {i1, i2, i3} and Transpose[{1, 2, 3}] are not the same shape. Part: The expression i3 cannot be used as a part specification. – user4660484 Jul 31 '18 at 3:21
• I'm working with version 11.3. The option Multicells was added only recently (and the syntax highlighting is still broken). Apperently, you work with an older version. In that case, you can replace all occurences of MeshCells[R, 2, "Multicells" -> True][[1, 1]] etc. by DeveloperToPackedArray[MeshCells[R, 2][[All, 1]]]. (This will result in a slight performance loss.) Notice that you have to do that also in the linked code for the heat method. – Henrik Schumacher Jul 31 '18 at 5:48
• Brilliant--that fix worked with version 11.2; who would have imagined there would be such a difference. One nuance worth mentioning: It is implied in your comment, but edges = DeveloperToPackedArray[MeshCells[R, 1][[All, 1]]]; was required instead of DeveloperToPackedArray[MeshCells[R, 2][[All, 1]]]. Thanks again. – user4660484 Aug 1 '18 at 3:06