Procedural Generation of Sculptural Forms

George W. Hart, Stony Brook University, http://www.georgehart.com, Jan 2008 
This notebook contains all the examples from my Bridges 2008 paper.

1. Built-in polyhedra, face representation, display, and stl output

PolyhedronData["Dodecahedron"]

Length[PolyhedronData[All]]

vertices[gc_]:=N[gc[[1,1]]]     (* vertices from graphics complex *)
vertices[PolyhedronData["Cube"]]

faces[gc_] := Map[vertices[gc][[#]]&, gc[[1,2,1]],{2}];   (* faces from gc *)
cube=faces[PolyhedronData["Cube"]]

view[obj_]:=Graphics3D[Map[Polygon,obj], Boxed->False];
view[cube]

stl[filename_,obj_]:=Export[filename,Graphics3D[Map[Polygon,obj]]];
stl["cube.stl", cube]

2. Tetrahedron, compounds

(* Test if given point p is above given triangle *)
above[p_,{v0_,v1_,v2_}]:=Module[{normal=Cross[v1-v0,v2-v1]},(p-v0).normal>0 ]

(* construct tetrahedron from four points*)
tetra[{p1_,p2_,p3_,p4_}]:=
If[above[p1,{p2,p3,p4}],{{p1,p2,p3},{p1,p3,p4},{p1,p4,p2},{p3,p2,p4}},
{{p1,p3,p2},{p1,p4,p3},{p1,p2,p4},{p3,p4,p2}}]

v4={{1,1,1},{1,-1,-1},{-1,1,-1},{-1,-1,1}};
view[tetra[v4]]

The following makes a Stella Octangula in which the large triangular faces intersect others, and there are internal walls. On some older SFF machines this "non-manifold" set of triangles confuses simple slicing algorithms and so will not build properly. But with recent machines the slicing algorithms are smarter and there is no problem. In other examples below, we will also use this technique of combining several overlapping objects, without further comment.

view[Join[tetra[v4],tetra[-v4]]]

To see which vertices go together in the following set of five tetrahedra, I separately printed out faces[PolyhedronData["Cube"]] and manually made a Schlegel diagram labeled with the particular vertex indices that Mathematica happens to use. In a future version of Mathematica, the indices might be permuted arbitrarily, requiring this code to be updated:

v=vertices[PolyhedronData["Dodecahedron"]];
view[Join[tetra[{v[[2]],v[[4]],v[[7]],v[[9]]}],
tetra[{v[[1]],v[[11]],v[[17]],v[[20]]}],
tetra[{v[[3]],v[[12]],v[[13]],v[[15]]}],
tetra[{v[[5]],v[[8]],v[[14]],v[[18]]}],
tetra[{v[[6]],v[[10]],v[[16]],v[[19]]}]]]

3. Rigid motions: translate, scale, and rotate

To translate an object, add a fixed offset to each XYZ coordinate in it. Map does this.  Our XYX's are down at "level 2" (inside faces which are inside objects).  Mathematical allows a third argument to Map of the form "{2}" which is a level specifier, just for this sort of thing.

translate[obj_,offset_]:=Map[(#+offset)&, obj,{2}];  (* offset should be {x,y,z} *)
view[Join[cube,translate[cube,{1.1,0,0}]]]

To scale an object, multiply each XYZ coordinate by a scale factor. We get this for free if we multiply the entire object by a scale factor, because in Mathematica  multiplication distributes across lists. Demonstrate this by joining two cubes together, one of which is scaled.

view[Join[cube,translate[0.5 cube,{.75,0,0}]]]

view[Apply[Join,Table[translate[i cube,{i^2/2,0,0}],{i,1,10}]]]

In a Linear Algebra class, they explain how a rotation matrix works. The built-in RotationMatrix function provides a matrix for rotating a given angle about a given axis. Multiply the matrix times a point {x,y,z} to get a new, rotated point. In Mathematica, the dot is used for this type of matrix-vector multiplication

rotate[obj_,angle_,axis_]:=Map[(RotationMatrix[angle,axis].#)&, obj,{2}];
view[Join[cube,rotate[cube,Pi/4,{0,0,1}]]]

Manipulate[view[Join[cube,rotate[cube,angle,{0,0,1}]]],{angle,0,Pi}]

view[Join[rotate[cube,Pi/4,{1,0,0}],
                   rotate[cube,Pi/4,{0,1,0}],
                   rotate[cube,Pi/4,{0,0,1}]]]

4. Poke

This method will create a unit-length vector in the direction of a given vector, by dividing out its length:

unit[vec_]:=vec/Sqrt[vec.vec];    (* will fail if given {0,0,0} *)

This method returns the mean of the elements in the list:

average[L_]:=Apply[Plus,L]/Length[L];

Our plan is to first make a function to poke one face, then map it on an object. The desired apex of a given poking can be found by starting at the center of the face and going in the direction of the normal a distance equal to the desired height. There will be n triangles generated from a given face. Each comes from two adjacent vertices and the apex. We will take position [[i]] and [[i+1]] of the face to get adjacent pairs, but the last pair wraps around as [[n]] and [[1]].  So plan ahead for this by appending the first vertex to the end of the list of vertices, to make a local variable, face1, that includes the wrap-around. Then a table of triples is easily formed:

pokeFace[face_,height_]:=Module[{apex,face1},
  apex=average[face]+height*unit[Cross[face[[1]]-face[[2]],face[[2]]-face[[3]]]];
  face1=Append[face,face[[1]]];
  Table[{face1[[i]], face1[[i+1]],apex}, {i,1,Length[face]}]
];

An object is a list of faces. So to construct the list of faces for a poked object, just poke each face and join the separate lists of triangles into a single list of triangles:

poke[obj_,height_]:=Apply[Join,Map[(pokeFace[#,height])&,obj]]

view[poke[cube,.7]]

Thgis makes a Stella Octangula without any faces intersecting, a single triangulated manifold:

view[poke[faces[PolyhedronData["Octahedron"]], .8]]

Manipulate[view[poke[faces[PolyhedronData["Icosahedron"]], height]],{height,-1,3}]

poke[obj_]:=poke[obj,1];
poke[obj_,height_,n_]:=Apply[Join,Map[(If[Length[#]==n,pokeFace[#,height],{#}])&,obj]]

view[poke[poke[cube], -.15]]

view[poke[faces[PolyhedronData["RhombicDodecahedron"]],.8]]

PolyhedronData["SmallRhombicuboctahedron"]

view[poke[faces[PolyhedronData["SmallRhombicuboctahedron"]],2,4]]

view[poke[poke[faces[PolyhedronData["Tetrahedron"]]]]]

5. Convex hull

This simple, but short and reasonably fast (for a quadratic time algorithm). It will fail if given five or more coplanar points to start with. But is is usually OK with coplanar points after it reaches a stage with some internal volume. There can also be occasional quirks from floating point roundoff. Note the search for unmatched edges is fast and easy if we use Mathematica's function definition mechanism in effect as a hash table.

(* To find unmateched edges, note that all matched edges show up twice, once in each order. This method, given a list of edges, returns edges not matched with their reverse *)
unmatchedEdges[List_]:=Module[{},
Clear[unmatched];
unmatched[x_]:=True;
Scan[(unmatched[Reverse[#]]=False)&,List];
Select[List,unmatched]
]

(* Utility method to delete 1st occurrence of x from a given List *)
remove[{},x_]:={};             
remove[L_,x_]:=If[First[L]==x,Rest[L],Prepend[remove[Rest[L],x],First[L]]]

(* Construct the three edges of given triangle *)
edgesOfTriangle[{v0_,v1_,v2_}]:={{v0,v1},{v1,v2},{v2,v0}}

(* 3D convex hull of list of points. *)
cHull[L_]:= Module[{faces,keepers,edges,newTriangles},        
If[Length[L]<4,Print["cHull needs at least 4 points"]];
faces=tetra[Take[L,4]];
For[i=5,i<=Length[L],i++,(
keepers=Select[faces,(!above[L[[i]],#])&];
edges=Apply[Join,Map[edgesOfTriangle,keepers]];
newTriangles=Map[({#[[2]],#[[1]],L[[i]]})&,unmatchedEdges[edges]];
faces=Join[keepers,newTriangles];
)];
faces
]

view[cHull[vertices[PolyhedronData["Icosahedron"]]]]

The bottom of this pyramid will be an n-gon divided into n-2 triangles.  But the division is irrelevant because they are coplanar, so does not affect the SFF result:

(* n-sided prism with base inscribed in unit circle in XY plane and apex at {0,0,1} *)
pyramid[n_]:=cHull[Prepend[Table[{Sin[2 Pi i/n],Cos[2 Pi i/n],0},{i,1,n}], {0,0,1}]//N];
view[pyramid[15]]

(* Prism of height 2 centered on origin.  Note the use of RandomSample to shuffle the points so it doesn't start with many in one plane *)
prism[n_]:=cHull[RandomSample[Join[Table[{Sin[2 Pi i/n],Cos[2 Pi i/n],1},{i,0,n-1}], 
Table[{Sin[2 Pi i/n],Cos[2 Pi i/n],-1},{i,0,n-1}]]//N]];
view[prism[17]]

The following coordinates were extracted from the 3D model at Mathworld, by Eric Weisstein:

durerPoints={{-0.70711`,-0.40825`,-0.44174`},{-0.70711`,0.40825`,0.44174`},{-0.43702`,-0.25231`,-0.7792`},{-0.43702`,0.25231`,0.7792`},{0.`,-0.8165`,0.44174`},{0.`,-0.50462`,0.7792`},{0.`,0.50462`,-0.7792`},{0.`,0.8165`,-0.44174`},{0.43702`,-0.25231`,-0.7792`},{0.43702`,0.25231`,0.7792`},{0.70711`,-0.40825`,-0.44174`},{0.70711`,0.40825`,0.44174`}};
view[cHull[durerPoints]]

6. Edge models

Note : The volume around each vertex is overlapped. Unfortunately, the orientation for different polyhedra in Mathematica's data base is not consistent, e.g., a cube and an octahedron do not have thewir symmetry axes aligned. So this code does not make symmetric models when two different polyhedra are chosen, even if they have identical symmetry. Hopefully a future version of Mathematica will be more sensitive to such uses and aligh all the polyhedra in a symmetry-consistent manner. 

Also, note that this code is based on Mathematica's GC format, rather than our list-of-faces data format.

edgeIndicesOfFace[face_]:=Select[Transpose[{face,RotateLeft[face]}],(#[[1]]<#[[2]])&];edgeIndices[gc_]:=Apply[Join, Map[edgeIndicesOfFace, gc[[1,2,1]]]];
edges[gc_]:= Map[vertices[gc][[#]]&, edgeIndices[gc],{2}]; 
translateVertices[name_,scale_,xyz_]:=Map[(xyz+#)&,scale vertices[PolyhedronData[name]]];
strut[name_,scale_,edge_]:= cHull[Join[translateVertices[name,scale,edge[[1]]],translateVertices[name,scale,edge[[2]]]]]
edgeModel[name_,nameVertex_,scale_]:=Apply[Join,Map[strut[nameVertex,scale,#]&,edges[PolyhedronData[name]]]]

view[edgeModel["Cube","Cube",.2]]

Manipulate[view[edgeModel[big,small,scale]], {big,PolyhedronData[All]},{small,PolyhedronData[All]},{scale,0.001,1}]

7. Fractal Tree

Each branch is a square prismatoid strut which surrounds the line connecting the two given points. Taking the Cross product with a "random" vector gives one direction orthogonal to the line, v1. Then v2 is constructed orthogonal to v1 and the line, and they provide a coordinate system for creating the two square bases of the strut.

branch[p1_,p2_]:=Module[{v1,v2,size=0.08Sqrt[(p1-p2).(p1-p2)]},
v1=size unit[Cross[(p2-p1),{1,2,3.4567}]];
v2=size unit[Cross[(p2-p1),v1]];
cHull[{p1+v1,p1-v1,p1+v2,p1-v2,p2+v1,p2-v1,p2+v2,p2-v2}]]

Standard recursive tree. Level 1 means just a single branch. The tree extendt starting at the given origin, in the "Z" direction of the Frenet frame, and the trunk has the given size. The three subtrees are rotated relative to the trunk by some fairly random angles I chose. Sow[] is used to collect all the branches into a set, returned by Reap[] in the main call, below.

tree[level_,origin_,frenet_,size_]:=Module[{top,M},
If[level==0,Return[]];
top=origin+size  frenet[[3]];
Sow[branch[origin,top]];
tree[level-1,top,RotationMatrix[Pi/5,{1,0,0}].frenet,.6 size];
tree[level-1,top,RotationMatrix[Pi/6,{-.5,.6,0}].frenet,.7 size];
tree[level-1,top,RotationMatrix[Pi/7,{-.5,-.6,0}].frenet,.5 size];
]

makeTree[level_]:=Reap[tree[level,{0,0,0},IdentityMatrix[3],1]][[2,1]]

view[makeTree[7]]

8. Fractal Polyhedron

This is kind of fun, once you understand how the parameters work. It would be much better if Mathematica's polyhedra were oriented with a consistent view to symmetry. The points[] function assembles together {vertices, edge midpoints, and face centers} so it is easy to choose one of the three options with a slider. The allSums[] function gives all nm possible sums obtained by adding one of the n elements of the first list to one of the m elements of the second list. The centers[] function applies this rrecursively to give a fractal pattern of points, uses as the centers of all the small polyhedra. The translates[] function gives many copies of a given object, one centered at each of the points in the given list of xyzs. The makeCell[] function optionally pokes (stellates) the unit cell.

points[name_] := With[{gc=PolyhedronData[name]}, {vertices[gc], Map[average,edges[gc]], Map[average,faces[gc]]}];
allSums[L1_,L2_]:=Apply[Join,Outer[Plus,L1,L2,1]] (* one from column A + one from column B *)
centers[depth_,type_,alpha_,points_]:=  (* recursive fn constructs pattern of ceners *)If[depth==0,{{0,0,0}},allSums[points[[type+1]],alpha centers[depth-1,type,alpha,points]]];
translates[obj_,xyzs_]:=Apply[Join,Map[translate[obj,#]&,xyzs]]; (* many copies of obj *)
makeCell[stell_,name_]:=If[stell, poke[faces[PolyhedronData[name]],2.83],faces[PolyhedronData[name]]];

To use this, first set the depth to 1. Set the two polyhedra and the VEF option. Optionally check the Stellated box f desired Next, adjust beta so the cells just barely overlap. Then set the depth to 2 and adjust alpha so te groups just barely contact. Then you can set the depth to 3 if desired.

Manipulate[Graphics3D[Map[Polygon,translates[(beta alpha^(depth-1)) makeCell[stellated,cell],centers[depth,VEF,alpha,points[name]]]] ],
{depth,0,3,1},{VEF,0,2,1},{{alpha,.1},0,1},{{beta,.1},0,1},{name,PolyhedronData[All]}, {cell,PolyhedronData[All]},{stellated,{False,True}}]

9. Marching Cubes

This is a very simple version, in which all the boundaries are voxel-separating squares. (No triangular corner clipping is used.) So it can result in non-manifold boundaries, with four faces meeting at an edge, or two sheets meeting at a vertex. But I have not observed any problems with the output being sliced on several different SFF machines. Note that the function g is defined internally as the intersection of f with the given cuboid, which is an indirect example of CSG, as made more explicit below.

boundary[f_,xmin_,xmax_,ymin_,ymax_,zmin_,zmax_,d_]:=Module[{g},
g[x_,y_,z_]:=If[xmin ≤ x≤ xmax &&ymin≤ y≤ ymax &&zmin≤ z≤ zmax,f[x,y,z],False];
Reap[
For[x=xmin,x≤xmax,x+=d,
For[y=ymin,y≤ymax,y+=d,
For[z=zmin,z≤zmax,z+=d,
If[g[x,y,z],( 
If[!g[x+d,y,z],Sow[{{x+d,y,z},{x+d,y+d,z},{x+d,y+d,z+d},{x+d,y,z+d}}]];
If[!g[x-d,y,z],Sow[{{x,y,z},{x,y,z+d},{x,y+d,z+d},{x,y+d,z}}]];
If[!g[x,y+d,z],Sow[{{x,y+d,z},{x,y+d,z+d},{x+d,y+d,z+d},{x+d,y+d,z}}]];
If[!g[x,y-d,z],Sow[{{x,y,z},{x+d,y,z},{x+d,y,z+d},{x,y,z+d}}]];
If[!g[x,y,z+d],Sow[{{x,y,z+d},{x+d,y,z+d},{x+d,y+d,z+d},{x,y+d,z+d}}]];
If[!g[x,y,z-d],Sow[{{x,y,z},{x,y+d,z},{x+d,y+d,z},{x+d,y,z}}]];
)]
]]]][[2,1]]]

fSphere[x_,y_,z_]:=(x^2+y^2+z^2)≤15^2
view[boundary[fSphere, 0,20,0,20,0,20,1]]

fWave[x_,y_,z_]:=Abs[Sin[x]+Sin[y]-z]<1;
view[boundary[fWave, 0,20,0,20,-3,3,.5]]

10. Menger Sponge

If you don't understand this, try cutting off two of your fingers on each hand:

fMenger[i_,j_,k_]:=Module[{id=IntegerDigits[i,3,7],jd=IntegerDigits[j,3,7],kd=IntegerDigits[k,3,7]},
For[n=1, n≤7, n++, 
If[(id[[n]]==jd[[n]]==1) || (id[[n]]==kd[[n]]==1) ||(jd[[n]]==kd[[n]]==1),
Return [False]]];
Return[True]]

Table[(max=3^depth-1;
view[boundary[fMenger,0,max,0,max,0,max,1]]), {depth,0,3}]

11. CSG

Any solid modeling textbook will explain about constructive solid geometry:

fDiff[x_,y_,z_]:=And[fMenger[x,y,z],Not[fSphere[x,y,z]]];
view[boundary[fDiff,0,26,0,26,0,26,1]]

12. 3D Hilbert curve

This is fairly routine, but I didn't look up a good source to cite. As implemented here, the makeHilbert[] function is used to set the definition of a voxel function, f, (defined only on integers) and then we apply Marching Cubes on f to get the boundary.

hilbert[level_,u_,v_,w_]:=Module[{},
If[level==0,Return[]];
hilbert[level-1,v,u,w];move[v];
hilbert[level-1,w,v,u];move[w];
hilbert[level-1,w,v,u];move[-v];
hilbert[level-1,u,-v,-w];move[u];
hilbert[level-1,u,-v,-w];move[v];
hilbert[level-1,-w,v,-u];move[-w];
hilbert[level-1,-w,v,-u];move[-v];
hilbert[level-1,-v,-u,w];
];
move[delta_]:=Module[{x,y,z, steps=2},
For[i=0,i<steps,i++,
turtle+=delta;
{x,y,z}=turtle;
fHilbert[x,y,z]=True;
]];
makeHilbert[level_]:=Module[{},
Clear[fHilbert];
fHilbert[x_,y_,z_]:=False;
fHilbert[0,0,0]:=True;
turtle={0,0,0};
hilbert[level,{1,0,0},{0,1,0},{0,0,1}];
]

Table[makeHilbert[i];
view[boundary[fHilbert,0,35,0,35,0,35,1]],{i,1,4}]