3D Design in Mathematica: Generative Jewelry

Happy July everyone. This month I'll share a bit of Mathematica code that has helped me to add some ordered randomness to my art and jewelry. And we'll design a new piece of generative jewelry!
 
Before that, let me share something else I've been working on. This June I gave myself a goal to share with the public each day one piece of jewelry I designed. Designing, refining, and posting took some effort, and I have completed my goal! In the process I designed some new jewelry that I really like, and I will be working over the coming months to develop them into coherent jewelry collections. I have been sharing on various social media accounts: @hanusadesign on instagram, facebook, twitter, and pinterest. Be sure to follow me to see my latest work!
 
Now on to the Mathematica tutorial. What do I mean by ordered randomness? Compare the following sets of 100 random points:
 

While the points seem rather "random" in both sets, you can see that the points on the left often land rather close to each other while the points on the right are pretty regularly spaced. This last observation should tell you that the points on the right are not actually "random"! Indeed, the points on the left are chosen randomly using Mathematica's RandomReal command. (Even though the numbers are generated by a deterministic process in a computer program and thus inherently not completely random, we will say that pseudorandom is close enough for our purposes.) The command

ptlist = RandomReal[{-1, 1}, {100, 2}];

generates 100 pairs of real numbers that all fall between $-1$ and $1$. On the other hand, the points on the right are chosen by a selection process in which we grow a list of points, starting with any point in the square. Then when we want to add a new point to the list, we generate a random point in the square and see if it is greater than a specified distance to any of the points already in our list. If so, it is added to our list; otherwise it is thrown out! Here is the code:

ptlist = {RandomReal[{-1, 1}, 2]};
While[Length[ptlist] < 100,
    newpt = RandomReal[{-1, 1}, 2];
    If[Min[Map[Norm[newpt - #] &, ptlist]] > .1,
       ptlist = Append[ptlist, newpt]]]

Running the same code multiple times gives different sets of points.

Generative jewelry

Let's take this idea and turn it into a nice piece of generative jewelry, in that we don't actually know what it is going to look like before our program completes!
 
My goal is to create a piece of jewelry that is a collection of overlapping rings. We'll first generate a set of points and then construct rings centered at those points. We will work to complete three tasks to improve the aesthetic appeal of the final product:

1. The points should be lie in a circle instead of a square.
2. The points should be generously spaced.
3. The rings should be different sizes

For Task 1, we will specify a shape and ensure that each point that we consider is a member of that region. (Here we choose a circle, but you can easily modify it to be ANY region.) For Task 2, we will modify our selection process from above to add more distance between the points (radius .21 instead of .1) and have fewer points (50 instead of 100). By using our selection process, the rings won't overlap too awkwardly and it will be pleasing to the eye. The values I have chosen are purely by trial and error with a view toward making the proportions of the final product look good.

shape = Disk[{0, 0}, 1];
ptlist = {RandomReal[{-1, 1}, 2]};
While[! RegionMember[shape, ptlist[[1]]],
    ptlist = {RandomReal[{-1, 1}, 2]}]
While[Length[ptlist] < 50, newpt = RandomReal[{-1, 1}, 2];
    If[Min[Map[Norm[newpt - #] &, ptlist]] > .21
       && RegionMember[shape, newpt],
       ptlist = Append[ptlist, newpt]]]
Graphics[{
    {Lighter[Blue, .8], Rectangle[{-1, -1}, {1, 1}]},
    {EdgeForm[{Black, Thick}], White, shape},
    {Purple, Point[ptlist]}
}]

Notice that we have also ensured that the first point is also in the circular region. The result of this code looks like this:
 

Now we want to build a ring centered at each point. A basic torus that has an outer radius of .14 and tube radius of .03 looks like this:

torus[coords_] := Module[{thickness = .14, innerradius = .03},
    ParametricPlot3D[{
       ( thickness + innerradius Cos[v]) Cos[u],
       ( thickness + innerradius Cos[v]) Sin[u],
       innerradius Sin[v]} + Append[coords, 0],
    {u, 0, 2 Pi}, {v, 0, 2 Pi}, Mesh -> None, PlotPoints -> 50]];

Mapping this function to the points in ptlist gives the following.
 

Now we need to attack Task 3, making the rings have different sizes. To do this we modify our torus function:

torus[coords_, thickness_] := Module[{innerradius = .03},
    ParametricPlot3D[{
       ( thickness + innerradius Cos[v]) Cos[u],
       ( thickness + innerradius Cos[v]) Sin[u],
       innerradius Sin[v]} + Append[coords, 0],
    {u, 0, 2 Pi}, {v, 0, 2 Pi}, Mesh -> None, PlotPoints -> 50]];

and generate a set of random outer radii. I first chose the outer radii to be between .1 and .15, but when I displayed the rings, they were disconnected, which is not what we wanted

thicknesses = RandomReal[{.1, .15}, 50]
Show[MapThread[torus, {ptlist, thicknesses}], PlotRange -> All]

So I modified the the outer radii to be between .09 and .17. I definitely re-ran the code multiple times until I was happy with the final product. Here you go:

thicknesses = RandomReal[{.09, .17}, 50]
Show[MapThread[torus, {ptlist, thicknesses}], PlotRange -> All]

Pro Tip: Re-running the code always gives you a new arrangement. When you get a random arrangement that you like, you need to save the data that created it, so you can recreate it next time!
 
In essence, when working to create a piece of generative art, you impose randomness until you find something that you like, at which time you need to save this input so that this data will not change again — which seems to be the opposite of random! It was a revelation when I realized that adding randomness to my art involved saving the random data I generated.
 

The Final Result

After exporting our final work to an STL, here is a rendering of our new piece of generative jewelry on Sketchfab:
 

And Shapeways gives the following beautiful rendering of our pendant in Raw Bronze:
 

And now you have the tools to make your own! If you have suggestions for how to modify this in an interesting way or for something else I should tackle, let me know. Happy July everyone!

Christopher Hanusa's portfolio is available online at

Purchase a copy of this and other 3D printed artwork
at Math Art Shop by Hanusa ≀ Design on Shapeways.

3D Design in Mathematica: Custom Dice

Today we continue our series in 3D Design in Mathematica. I'll teach you how to use the RegionPlot3D command to create custom dice. Let's first start off with a story I teach in my combinatorics class. Most people are familiar with a cubical die having six sides labeled with pips representing the numbers 1 through 6, each of which has the same probability of coming up when the die is rolled.
 

What happens when two dice are rolled and we consider their sum? The possible values are between 2 and 12 with a predictable distribution.
 

	1	2	3	4	5	6
1	2	3	4	5	6	7
2	3	4	5	6	7	8
3	4	5	6	7	8	9
4	5	6	7	8	9	10
5	6	7	8	9	10	11
6	7	8	9	10	11	12

Question of the day: Are there two other dice (not labeled 1 through 6) that when rolled and summed have the same chance of giving the sums in this table?
 
Answer of the day: Yes! Consider the pair of Sicherman dice. The first die has sides 1, 2, 2, 3, 3, and 4, while the second die has sides 1, 3, 4, 5, 6, and 8. Here is the distribution of sums of these two dice.
 

	1	2	2	3	3	4
1	2	3	3	4	4	5
3	4	5	5	6	6	7
4	5	6	6	7	7	8
5	6	7	7	8	8	9
6	7	8	8	9	9	10
8	9	10	10	11	11	12

You can count how many times each number appears to verify that the number of appearances of 2, 3, ... 12 are the same for each pair of dice. Let's modify these dice even further by subtracting one from each number on the first die (0, 1, 1, 2, 2, 3) and adding one to each number on the second die (2, 4, 5, 6, 7, 9), which doesn't change any sum one bit!
 

	0	1	1	2	2	3
2	2	3	3	4	4	5
4	4	5	5	6	6	7
5	5	6	6	7	7	8
6	6	7	7	8	8	9
7	7	8	8	9	9	10
9	9	10	10	11	11	12

This example is a great application of generating functions, because the existence of these dice is a direct result of the polynomial identity: \[(x^1+x^2+x^3+x^4+x^5+x^6)^2=(x^0+2x^1+2x^2+x^3)(x^2+x^4+x^5+x^6+x^7+x^9).\] Since these modified Sicherman dice are not easy to find at your local store (or at ANY store), we'll design our own copies as our task for today. If you're interested in dice, I suggest checking out The Dice Lab by Henry Segerman and Robert Fathauer where they've taken a mathematician's view toward perfecting dice.
 
Also, if you've missed an installment in the series, you can learn about designing a name ring, a geometric terrarium, or a country keychain.
 
Now on to the tutorial!
 

Step One: Determine a plan for the dice

One of the things about standard dice that I consider crucial is that the opposite sides of each die add to the same number. The 1 is opposite the 6, the 2 is opposite the 5, and the 3 is opposite the four. If we were to build a die out of paper and unfold the sides, we would get the following two-dimensional visualization:

We need to now place the numbers 0, 1, 1, 2, 2, 3 on the first die; I'll choose to continue the convention by ensuring that the 0 is opposite the 3 and each 1 is opposite a 2. (If you prefer to place the 1s opposite each other and the 2s opposite each other, you can modify our code to make your own dice!) I also plan to organize it so that the diagonals formed by the pips of the 2s are oriented in the same way.

The decisions in placing the numbers 2, 4, 5, 6, 7, 9 on the second die include the slant of the pips on the two and the orientation of the pips on the 6 and the 7. I think I'll just wing it and if my work turns out differently from my sketch then so be it!
 

Step Two: Understand how RegionPlot3D works

Today we'll be using Mathematica's RegionPlot3D command, which is one of the few ways in which Mathematica allows you to intersect three-dimensional regions. (Disappointingly, you can't intersect 3D MeshRegion objects like we used last time!)
 
RegionPlot3D outputs the three-dimensional object defined by a set of restrictions you specify along with boolean operators that indicate how your restrictions interact.
 
For instance, if you want to find the intersection of two spheres, you want points that are in both the first sphere AND the second sphere, so you would use the AND operator && like this:

RegionPlot3D[ x^2+y^2+z^2 <= 1 && (x-1)^2+y^2+z^2 <= 1,
 {x, 0, 1}, {y, -1, 1}, {z, -1, 1}]

If on the other hand you would like to find the union of two spheres, you want points that are in either the first sphere OR the second sphere, so you would use the OR operator || like this:

RegionPlot3D[ x^2+y^2+z^2 <= 1 || (x-1)^2+y^2+z^2 <= 1,
 {x, -1, 2}, {y, -1.5, 1.5}, {z, -1.5, 1.5}]

You must specify the bounds on the $x$, $y$, and $z$ variables where Mathematica will be looking for the valid points.
 
When we are making our dice, we will start with a basic cube, defined by points where the $x$, $y$, and $z$ values are between $-1$ and $1$:

RegionPlot3D[Abs[x] <= 1 && Abs[y] <= 1 && Abs[z] <= 1, {x, -1.1, 1.1}, {y, -1.1, 1.1}, {z, -1.1, 1.1}]

so that the equations of the sides of the cube are $x=1$, $x=-1$, $y=1$, $y=-1$, $z=1$, $z=-1$.
 

Now suppose we want to create the pips of the dice. I will do so by removing spheres from each of the faces. For instance, if we want to create one pip on the $z=1$ face, which would be centered at coordinates $(0,0,1)$, we will want the points that are in the cube with the points AND outside the sphere \[x^2+y^2+(z-1)^2=r^2,\] which we do with the code

pip = .04; RegionPlot3D[
  Abs[x] <= 1 && Abs[y] <= 1 && Abs[z] <= 1 &&
  (x^2 + y^2 + (z - 1)^2 >= pip)
  , {x, -1.1, 1.1}, {y, -1.1, 1.1}, {z, -1.1, 1.1},
  PlotStyle -> White, Mesh -> None, Axes -> False, Boxed -> False]

We have chosen $r=.2$, so $r^2=.04$. Note that we have changed the color using PlotStyle -> White, removed the lines from the surface with the option Mesh -> None, and removed the bounding wire cube and tickmarks using Axes -> False and Boxed -> False.
 
One of the issues with RegionPlot3D is the resolution of the file you create. In order to have a better approximation of the three-dimensional object that is determined by the inequalities you provide, you can use the PlotPoints option, which tells how much sampling you want Mathematica to do in the range you specify. Here you can see the difference in the result when you specify PlotPoints -> 10, PlotPoints -> 50, PlotPoints -> 100, and PlotPoints -> 200:
 

During the testing phase I used PlotPoints -> 50, but when creating the final STL file I specified PlotPoints -> 250 which took quite a while to compute.
 

Step Three: Rounding the Edges

We'll get back to the pips in a little bit, but first let's perfect the cube. You'll notice that when PlotPoints -> 200 is specified, the die starts to look more and more like a perfect cube. If you look at a normal die, you'll see that it is not a perfect cube—the edges are rounded. What sort of inequality restriction works to make a rounded edge? We need some sort of cylindrical constraint like \[(x-h)^2+(y-k)^2 < r^2.\] To line up with the faces of the cube, we need the radius of the cylinder to be the same as the distance from the point (h,k) to both faces. So in the simplest case we need something that looks like \[(x - .8)^2 + (y - .8)^2 < .2^2\] since it says that the distance to the $x=1$ face and to the $y=1$ face equals the radius $0.2$. (I tried a couple different radii and this one looked the best to my eyes.)
 
But simply including this cylindrical restriction would only include the part of the cube that is in the cylinder, which defeats the purpose.

RegionPlot[((x - .8)^2 + (y - .8)^2 < .2^2), {x, 0, 1.1}, {y, 0, 1.1}]

We need to instead merge this with two additional restrictions that either you are in the cylinder OR you satisfy $x<.8$ OR you satisfy $y<.8$. By doing this, the restriction on our cube is exactly the rounded bezel we wanted, as you can see here.

RegionPlot[((x-.8)^2+(y-.8)^2<.2^2||x<=.8||y<=.8)
  {x, 0, 1.1}, {y, 0, 1.1}]

We'll have to include this restriction for all twelve edges of the cube, which is coded as follows:

RegionPlot3D[Abs[x]<=1&&Abs[y]<=1&&Abs[z]<=1&&
 ((x-.8)^2+(y-.8)^2<.2^2||x<=.8||y<=.8)&&
 ((x-.8)^2+(y+.8)^2<.2^2||x<=.8||y>=-.8)&&
 ((x+.8)^2+(y-.8)^2<.2^2||x>=-.8||y<=.8)&&
 ((x+.8)^2+(y+.8)^2<.2^2||x>=-.8||y>=-.8)&&
 ((x-.8)^2+(z-.8)^2<.2^2||x<=.8||z<=.8)&&
 ((x-.8)^2+(z+.8)^2<.2^2||x<=.8||z>=-.8)&&
 ((x+.8)^2+(z-.8)^2<.2^2||x>=-.8||z<=.8)&&
 ((x+.8)^2+(z+.8)^2<.2^2||x>=-.8||z>=-.8)&&
 ((y-.8)^2+(z-.8)^2<.2^2||y<=.8||z<=.8)&&
 ((y-.8)^2+(z+.8)^2<.2^2||y<=.8||z>=-.8)&&
 ((y+.8)^2+(z-.8)^2<.2^2||y>=-.8||z<=.8)&&
 ((y+.8)^2+(z+.8)^2<.2^2||y>=-.8||z>=-.8)&&
 (x^2+y^2+(z-1)^2>=pip)
 ,{x,-1.1,1.1},{y,-1.1,1.1},{z,-1.1,1.1},
 PlotPoints->100,PlotStyle->White,Mesh->None,Axes->False,Boxed->False]

and looks like this:
 

Note that we could have used a similar method to round the corners, but I preferred this styling.
 

Step Four: Lay out the pips

Now onto the task of placing the pips. Let's start with my Modified Sicherman Die I with sides (0,1,1,2,2,3), which will be easier for us to code. We'll put the one-pip sides on the $x=-1$ and $y=1$ faces, the two-pip sides on the $x=1$ and $y=-1$ faces, and the three-pip side on the $z=-1$ face. Here is the Mathematica code for this:

(* 1 side *)
  ((x + 1)^2 + y^2 + z^2 >= pip) &&
 
(* 1 side *)
  (x^2 + (y - 1)^2 + z^2 >= pip) &&
 
(* 2 side *)
  ((x - 1)^2 + (y + .5)^2 + (z + .5)^2 >= pip) &&
  ((x - 1)^2 + (y - .5)^2 + (z - .5)^2 >= pip) &&
 
(* 2 side *)
  ((x + .5)^2 + (y + 1)^2 + (z + .5)^2 >= pip) &&
  ((x - .5)^2 + (y + 1)^2 + (z - .5)^2 >= pip)
 
(* 3 side *)
  ( (x + .5)^2 + (y + .5)^2 + (z + 1)^2 >= pip) &&
  ( x^2 + y^2 + (z + 1)^2 >= pip) &&
  ( (x - .5)^2 + (y - .5)^2 + (z + 1)^2 >= pip) &&

As you can see, the three-pip side has three spheres, centered at $(-.5,-.5,-1)$, $(0,0,-1)$, and $(.5,.5,-1)$, all indeed lying on the face $z=-1$.
 
For my Modified Sicherman Die II with sides (2,4,5,6,7,9), we place the two- and nine-pip sides on the $z$ faces, the five- and six-pip sides on the $x$ faces, and the four- and seven-pip sides on the $y$ faces:

(* 2 side *)
  ((x + .5)^2 + (y + .5)^2 + (z - 1)^2 >= pip) &&
  ((x - .5)^2 + (y - .5)^2 + (z - 1)^2 >= pip) &&
 
(* 9 side *)
  ((x + .5)^2 + (y + .5)^2 + (z + 1)^2 >= pip) &&
  ((x + .5)^2 + y^2 + (z + 1)^2 >= pip) &&
  ((x + .5)^2 + (y - .5)^2 + (z + 1)^2 >= pip) &&
  (x^2 + (y + .5)^2 + (z + 1)^2 >= pip) &&
  (x^2 + y^2 + (z + 1)^2 >= pip) &&
  (x^2 + (y - .5)^2 + (z + 1)^2 >= pip) &&
  ((x - .5)^2 + (y + .5)^2 + (z + 1)^2 >= pip) &&
  ((x - .5)^2 + y^2 + (z + 1)^2 >= pip) &&
  ((x - .5)^2 + (y - .5)^2 + (z + 1)^2 >= pip) &&
 
(* 6 side *)
  ((x - 1)^2 + (y + .5)^2 + (z + .5)^2 >= pip) &&
  ((x - 1)^2 + (y)^2 + (z + .5)^2 >= pip) &&
  ((x - 1)^2 + (y - .5)^2 + (z - .5)^2 >= pip) &&
  ((x - 1)^2 + (y + .5)^2 + (z - .5)^2 >= pip) &&
  ((x - 1)^2 + (y)^2 + (z - .5)^2 >= pip) &&
  ((x - 1)^2 + (y - .5)^2 + (z + .5)^2 >= pip) &&
 
(* 5 side *)
  ((x + 1)^2 + y^2 + z^2 >= pip) &&
  ((x + 1)^2 + (y + .5)^2 + (z + .5)^2 >= pip) &&
  ((x + 1)^2 + (y + .5)^2 + (z - .5)^2 >= pip) &&
  ((x + 1)^2 + (y - .5)^2 + (z + .5)^2 >= pip) &&
  ((x + 1)^2 + (y - .5)^2 + (z - .5)^2 >= pip) &&
 
(* 4 side *)
  ((x + .5)^2 + (y - 1)^2 + (z + .5)^2 >= pip) &&
  ((x + .5)^2 + (y - 1)^2 + (z - .5)^2 >= pip) &&
  ((x - .5)^2 + (y - 1)^2 + (z + .5)^2 >= pip) &&
  ((x - .5)^2 + (y - 1)^2 + (z - .5)^2 >= pip) &&
 
(* 7 side *)
  ((x + .5)^2 + (y + 1)^2 + (z + .5)^2 >= pip) &&
  ((x)^2 + (y + 1)^2 + (z - .5)^2 >= pip) &&
  ((x + .5)^2 + (y + 1)^2 + (z - .5)^2 >= pip) &&
  ((x)^2 + (y + 1)^2 + (z)^2 >= pip) &&
  ((x - .5)^2 + (y + 1)^2 + (z + .5)^2 >= pip) &&
  ((x)^2 + (y + 1)^2 + (z + .5)^2 >= pip) &&
  ((x - .5)^2 + (y + 1)^2 + (z - .5)^2 >= pip)

All these spheres were placed by hand with coordinates in the set $\{-1,-.5,0,.5,1)$. Here's a fun picture of what it looks likes when all the pips' spheres are added to the cube!
 

Step Five: Get them printed!

To serve as a recap, here is the entirety of the code that generates my Modified Sicherman Die I.

pip = .04;
die = RegionPlot3D[
  Abs[x] <= 1 && Abs[y] <= 1 && Abs[z] <= 1 &&
  ((x-.8)^2+(y-.8)^2<.2^2||x<=.8||y<=.8)&&
  ((x-.8)^2+(y+.8)^2<.2^2||x<=.8||y>=-.8)&&
  ((x+.8)^2+(y-.8)^2<.2^2||x>=-.8||y<=.8)&&
  ((x+.8)^2+(y+.8)^2<.2^2||x>=-.8||y>=-.8)&&
  ((x-.8)^2+(z-.8)^2<.2^2||x<=.8||z<=.8)&&
  ((x-.8)^2+(z+.8)^2<.2^2||x<=.8||z>=-.8)&&
  ((x+.8)^2+(z-.8)^2<.2^2||x>=-.8||z<=.8)&&
  ((x+.8)^2+(z+.8)^2<.2^2||x>=-.8||z>=-.8)&&
  ((y-.8)^2+(z-.8)^2<.2^2||y<=.8||z<=.8)&&
  ((y-.8)^2+(z+.8)^2<.2^2||y<=.8||z>=-.8)&&
  ((y+.8)^2+(z-.8)^2<.2^2||y>=-.8||z<=.8)&&
  ((y+.8)^2+(z+.8)^2<.2^2||y>=-.8||z>=-.8)&&
  ((x + 1)^2 + y^2 + z^2 >= pip) &&
  (x^2 + (y - 1)^2 + z^2 >= pip) &&
  ((x - 1)^2 + (y + .5)^2 + (z + .5)^2 >= pip) &&
  ((x - 1)^2 + (y - .5)^2 + (z - .5)^2 >= pip) &&
  ((x + .5)^2 + (y + 1)^2 + (z + .5)^2 >= pip) &&
  ((x - .5)^2 + (y + 1)^2 + (z - .5)^2 >= pip)
  ( (x + .5)^2 + (y + .5)^2 + (z + 1)^2 >= pip) &&
  ( x^2 + y^2 + (z + 1)^2 >= pip) &&
  ( (x - .5)^2 + (y - .5)^2 + (z + 1)^2 >= pip) &&
  , {x, -1.1, 1.1}, {y, -1.1, 1.1}, {z, -1.1, 1.1}, PlotPoints -> 250,
  PlotStyle -> White, Mesh -> None, Axes -> False, Boxed -> False]

After exporting the STL files for my two dice (using PlotPoints -> 250), here is a rendering of the final version of my modified Sicherman dice on Sketchfab:
 

I printed them on the local Lulzbot Mini 3D printers that the Queens College Art department has so graciously loaned us. They look good too!
 

You can buy your own pair of modified Sicherman dice on Shapeways or print your own copy of the STL files I posted on Thingiverse. I don't guarantee that my modified Sicherman dice are unbiased because by removing material from the sides of the dice they no longer symmetrical.
 
The same technique above shows that a 4-sided die with sides (1,4,4,7) and a 9-sided die with (1,2,2,3,3,3,4,4,5) would also have the same distribution of sums—now THAT would be a unique pair of dice! What other dice would you create?

Christopher Hanusa's portfolio is available online at

Purchase a copy of this and other 3D printed artwork
at Hanusa ≀ Design on Shapeways.

3D Design in Mathematica: Costa Rica Keychain

Hi everyone! Welcome back to another 3D Design in Mathematica Tutorial. Last time I showed how to create a polyhedral terrarium. Today I will use Mathematica's MeshRegion command to design a keychain shaped like a country. I had the opportunity to visit Costa Rica in February and especially enjoyed my visit to the cloud forests of Monteverde, so let's choose Costa Rica for our example.
 

Step One: Import the shape of Costa Rica

We'll be relying on one of the strengths of Mathematica and the whole Wolfram suite—its vast amount of curated data. The command CountryData can access diverse information about countries such as their GDP, their population, their borders, their languages, and even their flags. For instance, we can access the shape of Costa Rica by using the following command.

CountryData["CostaRica", "Shape"]

What we would really like are the coordinates of the vertices of this shape, which we can acquire using the command

vertices =
     CountryData[Entity["Country", "CostaRica"], "Polygon"][[1, 1, 1]];

The first few elements of this list are

{{11.0785, -85.6911}, {11.2142, -85.6144}, {11.21, -85.5642},
  {11.1667, -85.5264}, ...

which are the longitude and latitude of each of the coordinates.
 

Step Two: Simplify the boundary

Keeping an eye toward printability, I realized that Costa Rica's actual land and sea borders are too complicated to be able to be printed in metal. (Shapeways has certain tolerances that must be adhered to for them to accept your file.) So I needed to find a way to smooth the border to reduce the complexity. After searching around for quite a while, I learned about a couple well-known simplification algorithms including Mike Bostock's awesome visualization implementing the Visvalingam-Whyatt algorithm.
 
I found this list of implementations of Mike's simplify code and was able to run the python implementation on repl.it. I needed to take the data from Mathematica and put it into the right python format to input into simplify, which would take a pair of coordinates {11.0785, -85.6911} and output them as the string {'x':11.0785, 'y':-85.6911}. So I used Mathematica's String operations.

StringJoin @@
  Map["{'x':" <> ToString[#[[1]]] <>
     ", 'y':" <> ToString[#[[2]]] <> "}," &, vertices]

I input the resulting code:

[{'x':11.0785, 'y':-85.6911}, {'x':11.2142, 'y':-85.6144}, ... ,
   {'x':11.0782, 'y':-85.6907}]

into simplify with the parameter 0.05. The output was again in python format, which I needed to reconvert into Mathematica format by deleting extraneous characters:

StringDelete[StringDelete[newpoints, "'y':"], "'x':"]

We can compare the before and after of this process:
 

Notice that there are still some jagged parts of the polygon, so I decided to make a few manual adjustments too. The best way I have found to do this is to use the ToolTip command

Graphics[{Polygon@simplified, Red, Map[Tooltip[Point@#, #] &, simplified]}]

which displays a vertex's coordinates when you hover over it so that I know I am removing the correct vertices.
 

By the end of this process, I had reduced the original polygon with over 1400 vertices to a polygon with 55 vertices that still contains the essence of the shape of Costa Rica.

simplified = {{ -85.6911, 11.0785}, { -85.6144, 11.2142}, { -84.9387, 10.9541}, { -84.7104, 11.091}, { -84.4801, 10.973}, { -84.3577, 10.9965}, { -84.2041, 10.7839}, { -84.0131, 10.7908}, { -83.9293, 10.7087}, { -83.6779, 10.7938}, { -83.6978, 10.8801}, { -83.6401, 10.9096}, { -83.3902, 10.3557}, { -82.7786, 9.66272}, { -82.62, 9.57017}, { -82.6712, 9.49412}, { -82.8549, 9.57112}, { -82.9361, 9.47312}, { -82.9359, 9.07847}, { -82.7118, 8.92455}, { -82.9192, 8.76513}, { -82.83, 8.63599}, { -82.8396, 8.48029}, { -83.0517, 8.33343}, { -83.1454, 8.36625}, { -83.1207, 8.60242}, { -83.3745, 8.74827}, { -83.4645, 8.72645}, { -83.2844, 8.54143}, { -83.2774, 8.3865}, { -83.5892, 8.45915}, { -83.7345, 8.59876}, { -83.5704, 8.84997}, { -83.6403, 9.05175}, { -84.2461, 9.49397}, { -84.6181, 9.58172}, { -84.703, 9.92881}, { -84.774, 9.99558}, { -85.2166, 10.1833}, { -85.1987, 10.0272}, { -84.9217, 9.9254}, { -84.8642, 9.82857}, { -85.1111, 9.559}, { -85.2539, 9.78983}, { -85.6684, 9.90357}, { -85.8752, 10.3566}, { -85.7732, 10.4469}, { -85.811, 10.5156}, { -85.6993, 10.6107}, { -85.6584, 10.766}, { -85.9497, 10.8896}, { -85.8843, 10.9498}, { -85.7141, 10.9201}, { -85.755, 11.0254}, { -85.6907, 11.0782}}

Step Three: Define an interior

I would like to make a raised lip around the outside of the keychain, so I will use some code I wrote last year that creates a boundary of a fixed thickness on the inside or outside of a given polygon.
 
Before I apply it, I'll remove some of the peninsulas from our polygon and some of the redundant vertices and visualize the difference:

minimal = { { -85.6144, 11.2142}, { -84.9387, 10.9541}, { -84.7104, 11.091}, { -84.4801, 10.973}, { -84.3577, 10.9965}, { -84.2041, 10.7839}, { -84.0131, 10.7908}, { -83.9293, 10.7087}, { -83.5979, 10.8238}, { -83.3902, 10.3557}, { -82.7786, 9.66272}, { -82.9361, 9.47312}, { -82.9359, 9.07847}, { -82.9192, 8.76513}, { -83.1207, 8.60242}, { -83.3745, 8.74827}, { -83.5704, 8.84997}, { -83.6403, 9.05175}, { -84.2461, 9.49397}, { -84.6181, 9.58172}, { -84.703, 9.92881}, { -84.774, 9.99558}, { -85.2166, 10.1833}, { -85.1987, 10.0272}, { -85.2539, 9.78983}, { -85.6684, 9.90357}, { -85.8752, 10.3566}, { -85.7732, 10.4469}, { -85.811, 10.5156}, { -85.6584, 10.766}, { -85.7141, 10.9201}};
Graphics[{Polygon@simplified, Gray, Polygon@minimal]}]

Now let's describe the algorithm to create a polygon nested inside our bounding shape—we'll be using geometry and trigonometry! We define the vertices of the interior polygon in relation to the vertices of the exterior polygon. The interior vertex is found along the angle bisector of the edges incident with the exterior vertex. We determine its final position by ensuring that the distance to each exterior edge is the prescribed distance ep.
 
Let me give you the code and then dissect it.

ep = .15;
boundarypoints = Reverse@minimal;

interiorhalfangles =
  Table[FullSimplify[
   Mod[
    Apply[ArcTan,
     (boundarypoints[[i]] -
      boundarypoints[[Mod[i + 1, Length[boundarypoints], 1]]])
    ] -
    Apply[ArcTan,
     (boundarypoints[[Mod[i + 2, Length[boundarypoints], 1]]] -
      boundarypoints[[Mod[i + 1, Length[boundarypoints], 1]]])
    ],
   2 Pi]],
  {i, Length[boundarypoints]}]/2;

edgedirections =
  Table[FullSimplify[
   Apply[ArcTan,
    boundarypoints[[Mod[i + 2, Length[boundarypoints], 1]]] -
    boundarypoints[[Mod[i + 1, Length[boundarypoints], 1]]]]],
  {i, Length[boundarypoints]}];

interiorvectors =
  Map[AngleVector, interiorhalfangles + edgedirections];

interiorvertices =
  MapThread[#1 + ep/Sin[#3] #2 &,
   {RotateLeft[boundarypoints], interiorvectors, interiorhalfangles}]

We choose our thickness (ep for epsilon) to be .15. Notice that the following command is to invert the order of the points on the bounding polygon. If we had left the list as is, the new polygon would have been constructed on the outside of the original polygon instead of on the inside.
 
Then I find the interior half angles by using a special overloading of the ArcTan command. The expected functionality of ArcTan is that ArcTan applied to a single number finds the angle whose tangent is that number. More useful in this context is the two argument version of ArcTan:
 
Useful Command: ArcTan[x,y] gives the angle that the line through (x,y) makes with the positive x-axis.
 
I have used ArcTan to find the angles that each of the directed edges make with the positive x-axis and taking their difference modulo 2π so that the number is between 0 and 2π. Another command that seemed useful when I started programming was VectorAngle, which finds the angle between two vectors, but that gives the unsigned angle instead of the signed angle between the two vectors.
 
If we add the interior half angle to the angle that the corresponding edge makes with the x-axis, we get the direction of the vector leaving the exterior vertex toward the interior vertex. Applying the command AngleVector gives the unit vector in that direction. (It gives me warm fuzzy feelings to know that Mathematica has both a command VectorAngle and a command AngleVector!)
 
The interior vertices are finally calculated using trigonometry to calculate the distance along this unit vector from the exterior vertex.
 
The final result is a polygon that nests perfectly inside the original Costa Rica shape.

Graphics[{Polygon@simplified, Red, Polygon@interiorvertices}]

Step Four: Create a Mesh Object

Now that we have our exterior polygon and our interior polygon, we will convert these Graphics objects into MeshRegion objects, since this seems to be the underlying structure necessary to export three dimensional objects to STL files. The first step is to use DiscretizeGraphics to convert the polygons into MeshRegion objects.

interior = DiscretizeGraphics@Graphics@Polygon@interiorvertices
exterior = DiscretizeGraphics@Graphics@Polygon@simplified

Then we use RegionDifference to excise the interior from the exterior.

lip = RegionDifference[exterior, interior]

The last pieces we need are the boundaries of these regions.

exteriorbdry = RegionBoundary@exterior
interiorbdry = RegionBoundary@interior

With these building blocks, we can create our three-dimensional model. We define four variables that are the heights of each of the layers, making sure that they are numerical / decimal instead of exact / infinitely precise numbers. There are errors if you use exact numbers like '0' instead of '0.'. I am not sure why.
 
Pro tip: Notice that I have defined these variables outside the subsequent code. By separating the definition of the variables, it makes it easier to modify the rest of the code later if you need to change a parameter for aesthetic or printing reasons.

extbottom = 0.;
intbottom = 0.1;
inttop = 0.2;
exttop = 0.3;

And then build the pieces that comprise the hull of the model by using a RegionProduct command.

Show[
  RegionProduct[interior, Point[{{intbottom}}]],
  RegionProduct[interior, Point[{{inttop}}]],
  RegionProduct[interiorbdry, Line[{{extbottom}, {intbottom}}]],
  RegionProduct[interiorbdry, Line[{{inttop}, {exttop}}]],
  RegionProduct[lip, Point[{{extbottom}}]],
  RegionProduct[lip, Point[{{exttop}}]],
  RegionProduct[RegionBoundary@exteriorbdry,
                 Line[{{extbottom}, {exttop}}]]]

You can see each of these objects individually and all assembled together.

Step Five: Add a loop

Since we want to use this object as a keychain, we need to add in a ring onto which we can clip a key ring. To do this we will be adding a torus into the object, which has this general form.

thickness = .15;
innerradius = .8;
outerradius = 2;
xcoord = -84;
ycoord = 10.75;
zcoord = Mean[{extbottom, exttop}];
loop = ParametricPlot3D[{
  (outerradius thickness + innerradius thickness Cos[v]) Cos[u] +
     xcoord,
  (outerradius thickness + innerradius thickness Cos[v]) Sin[u] +
     ycoord + 2 thickness,
  thickness Sin[v] + zcoord},
  {u, 0, 2 Pi}, {v, 0, 2 Pi}, Mesh -> None, PlotPoints -> 200];

I chose the ring to have the same thickness as the body and played with the other parameters, especially the x-, y-, and z-translations to figure out where best on the model to place it. One reason why I think my choices work well is because the circle is basically just North of the center of gravity—the map of Costa Rica will have North oriented upward when dangling from the keychain!
 

Step Six: Send to the printer!

Now let's print out the model! (Remember that there is more detail on the prototyping and printing processes at the bottom of the name ring design post.) One thing we have to take into consideration is that stainless steel requires a 1 mm thickness everywhere, which may mean that you will have to change the thickness parameters above if you want to shrink the keychain much smaller. Here is what Shapeways gives as a 3D rendering of the keychain
 

And here is a 3D rendering of the keychain by Sketchfab
 

Now you have the tools to make your own keychain for any country or administrative region. What countries are on your bucket keychain list? Until next time!

Christopher Hanusa's portfolio is available online at

Purchase a copy of this and other 3D printed artwork
at Hanusa ≀ Design on Shapeways.

3D Design in Mathematica: Polyhedral Terrarium

Today I continue my series of Mathematica Tutorials for 3D printing. Last time I showed how to create a custom-designed name ring. In this post I will discuss how to design a polyhedral terrarium. You may have noticed that geometric objects are very fashionable these days. A quick Google Image search for Geometric Terrarium gives a wide variety of styles and manufacturers.
 

We're going to make our own terrarium using Mathematica. I'm inspired to write about this because I'm giving a talk in May at Construct3D about 3D printing techniques in Mathematica that I teach my students, and this is a case in point—my student Isaac Deonarine used similar techniques to create a Dodecahedral Terrarium. Click here to interact with a 3D rendering of his work:

or see the Geometric Terrarium at Shapeways. Now, on to the tutorial!
 

Step One: Choose a Polyhedron

The Mathematica team curated a large amount of polyhedra that can be accessed using the command PolyhedronData. The full list of 195 entries can be accessed using the command PolyhedronData[All], and all 195 polyhedra can be visualized by typing

Map[PolyhedronData, PolyhedronData[All]]

Alternatively, Mathematica's thorough Documentation Center suggests using the code

Manipulate[
 Column[{PolyhedronData[g], PolyhedronData[g, p]}],
  {g, PolyhedronData[All]},
  {p, Complement @@ PolyhedronData /@ {"Properties", "Classes"}}
]

to generate an interactive interface for navigating the data:

I've chosen to work with the Snub Cube.

The snub cube is one of my favorite polyhedra because it has chirality—there are two distinct snub cubes based on the way it is built— (also an important concept in chemistry...) and because it cannot be built using the mathematical building blocks ZomeTools!
 

Step Two: Develop the schematics

Now let's gather some data about snub cubes. It is useful to know the coordinates of the vertices, the pairs of vertices that form the edges, and faces of the polyhedron. These are all simple to request.

vertices = N@PolyhedronData["SnubCube", "VertexCoordinates"]
edges = N@PolyhedronData["SnubCube", "Edges"]
faces = N@PolyhedronData["SnubCube", "Faces"]

The N@ prefix asks for numerical approximations of these coordinates so that we are not dragging around the exact coordinates which involve roots of the polynomial \[32 x^6-32 x^4-12 x^2-1.\]
Fun Fact: N@ is also a colloquialism from Pittsburgh, PA, where I grew up. It is a filler that means "and so on" and might be used in a sentence such as "Yinz goin' dahntahn n'at?".

What I would like to do is to remove some of the faces. If you look at the structure of the faces variable, you see that it is a GraphicsComplex object listing its coordinates and the list of vertices involved in each face. A very helpful thing to do is to put the labels of the vertices on the vertices to understand how the polyhedron is constructed. We can do that by putting the vertex index at the coordinates of the vertex.

coords = faces[[1]]
facelist = faces[[2, 1]]
Graphics3D[{
  Table[Text[Style[i, Large], verts[[i]]], {i, Length[verts]}],
  GraphicsComplex[coords, Polygon@facelist]
}, Boxed -> False]

which looks like this:

The variable facelist contains the list of vertices for each of the polygons:

{{3, 1, 17}, {3, 17, 9}, {3, 19, 2}, {3, 9, 19}, {1, 4, 20}, {1, 20, 11}, {1, 11, 17}, {2, 19, 12}, {2, 18, 4}, {2, 12, 18}, {4, 18, 10}, {4, 10, 20}, {17, 11, 13}, {19, 9, 15}, {18, 12, 14}, {20, 10, 16}, {9, 21, 15}, {11, 23, 13}, {12, 24, 14}, {10, 22, 16}, {13, 23, 7}, {13, 7, 21}, {15, 21, 5}, {15, 5, 24}, {16, 22, 6}, {16, 6, 23}, {14, 24, 8}, {14, 8, 22}, {21, 7, 5}, {23, 6, 7}, {24, 5, 8}, {22, 8, 6}, {1, 3, 2, 4}, {21, 9, 17, 13}, {24, 12, 19, 15}, {10, 18, 14, 22}, {11, 20, 16, 23}, {8, 5, 7, 6}}

For example, in the above figure you can see the triangular face with vertices {3, 17, 9} on the left and the square face with vertices {8, 5, 7, 6} on the right.

I am going to define some vertex sets in order to simplify some future commands. Let's define the sets of vertices on the top and bottom faces of the snub cube and the vertices on the upper and lower halves of the snub cube.

topvs = {12, 15, 19, 24};
bottomvs = {16, 20, 11, 23};
uppervs = Quiet@Flatten@Position[vertices, _?(#[[3]] > 0 &)];
lowervs = Quiet@Flatten@Position[vertices, _?(#[[3]] < 0 &)];

I'd like to only keep the faces at the bottom of the snub cube, which I could have done by hand, but because of the previous definitions, I can calculate algorithmically by asking which faces in facelist do not contain any of the vertices in the upper half:

newfacelist = Select[facelist, Intersection[uppervs, #] == {} &]
Graphics3D[{
  GraphicsComplex[coords, Polygon@newfacelist]
}, Boxed -> False]

which looks like this:

When we removed the faces, we also removed the incident edges, and I'd like to add many of them back in. In particular, I'd like to reinsert the lost edges that are not are not touching the top face and also not edges of the faces of the base.

lostedges = Select[edgelist,
  (Length@Intersection[lowervs, #] < 2 &&
   Length@Intersection[topvs, #] == 0) &]
Graphics3D[{
  GraphicsComplex[coords, Polygon@newfacelist],
  Thick, GraphicsComplex[coords, Line@lostedges]
}, Boxed -> False]

This code yields the following image. The basic structure of our terrarium is finally coming in view!

 

Step Three: Make it 3D printable

The structure of what we want to 3D print is there, but we can't print it because each face is only two-dimensional and each edge is only one-dimensional! So we'll need to add some thickness to the faces and the edges. In fact, we won't be adding thickness to a face as much as we will be constructing a new solid that has two layers of faces. We will be generating a new GraphicsComplex object, which means we need to mark down all vertices that form the corners of our object, and keep track of which of these vertices are connected to form the polygonal faces.

Just like blowing up a balloon appears to dilate the surface of the balloon by a constant factor, we can create an outer layer of faces 10% bigger than the original faces by multiplying all coordinates by a factor of 1.1; we save these coordinates as dilatedcoords. Thinking of the original vertices as vertices 1 through 24, we consider these dilated vertices to be vertices 25 through 48; furthermore, they have the exact same incidences as the vertices 1 through 24. So determining the vertex indices in the list of dilated faces is as simple as adding 24 to the index of each of the original vertices.

dilatedcoords = 1.1 coords;
dilatedfacelist = newfacelist + 24;
allcoords = Join[coords, dilatedcoords];
allfaces = Join[newfacelist, dilatedfacelist];
Graphics3D[{GraphicsComplex[allcoords, Polygon@allfaces]}, Boxed -> False]

Combining the original and dilated coordinates and faces gives the inner and outer shell of polygonal faces, but this is not a closed object. We also need to complete the base by adding in quadrilaterals that bridge the shells. To do this, we (manually) keep track of the indices of the original vertices that work around the boundary of the base. The quadrilaterals we create need to join two adjacent boundary vertices with their dilates (remember they are indexed by 24 more than their non-dilated counterparts.

boundaryedges = {{1, 17}, {17, 13}, {13, 7}, {7, 6}, {6, 22}, {22, 10}, {10, 4}, {4, 1}};
boundaryfaces = Map[
  {#[[1]], #[[2]], #[[2]] + 24, #[[1]] + 24} &,
  boundaryedges];
Graphics3D[{
  GraphicsComplex[allcoords, Polygon@allfaces],
  GraphicsComplex[allcoords, Polygon@boundaryfaces]
}, Boxed -> False]

Putting these all together gives the base; we visualize the result:

Now we need to construct the three dimensional version of the edges. I want to make thin cylindrical edges for the black lines on the sketch above, and I also want to curve the top of each of the boundary edges to soften the edges. I introduced some parameters to be able to play around with them to work out the best visualization. The value midpt is the coordinate multiplier. It is 1.05 since that is halfway between the original coordinate (multiplier 1.0) and the dilated coordinate (multiplier 1.1). The big radius is the radius of the cylinder along the `boundary edges'. I was originally surprised that .05 was not the correct value here. The small radius is the radius of the cylinder along the `lost edges' (the black lines above). Map then puts a Tube (cylinder) along each of those edges at the specified radii.

midpt = 1.05; bigr = .0615; smallr = .04;
tubes = Join[
  Map[Tube[#, bigr] &,
   Table[midpt {coords[[boundaryedges[[i, 1]]]],
    coords[[boundaryedges[[i, 2]]]]},
    {i, Length[boundaryedges]}]],
  Map[Tube[#, smallr] &,
   Table[midpt {coords[[lostedges[[i, 1]]]],
    coords[[lostedges[[i, 2]]]]},
   {i, Length[lostedges]}]]
];

Just as when we constructed the name ring, we need to finish the cylinders with spheres at each endpoint. We need spheres of the same radius as the radius of the adjacent cylinders.

spheres = {
  Map[Sphere[#, bigr] &,
   midpt coords[[Complement[lowervs, bottomvs]]]],
  Map[Sphere[#, smallr] &,
   midpt coords[[Complement[uppervs, topvs]]]]
}

Now put everything together and export the model to an STL file.

terrarium = Graphics3D[{
  GraphicsComplex[allcoords, Polygon@allfaces],
  GraphicsComplex[allcoords, Polygon@boundaryfaces],
  tubes, spheres
}, Boxed -> False]
Export[NotebookDirectory[] <> "terrarium.stl", terrarium]

After exporting, we can upload it where we need to. I've uploaded it to Sketchfab so you can play with a rendering of the final product!

Step Four: Send to the printer!

Now let's print out the model! (Note that I provide much more detail on the prototyping and printing processes at the bottom of the name ring design post.)

I highly suggest a final print of your terrarium in glazed ceramic because it will be able to withstand the elements and the cleaning required by dirt and plants. At the same time, I would also suggest prototyping on a local 3D printer instead of prototyping through Shapeways. At 4.5 inches by 4.5 inches by 3.5 inches, a basic print is already approaching $60! I printed it on the Lulzbot Mini printer that my colleague in the Queens College Art department is graciously loaning me, and, after 9 hours of printing, here is the stunning result: (click to zoom)

Pro Tip: The standard Cura software that comes with the Lulzbot was not able to correctly process my file—it either completely fills in the bottom half of the model or does not connect the edges together. So I needed to do some post-processing of the STL file. I uploaded it to Netfabb's online STL repair service and the Lulzbot was able to print the STL output from there with no issues.

Shapeways has produced a rendering of what it anticipates that a glazed ceramic version of my file will look like.

I can't wait to get one for myself and house some succulents inside. If you do your own modifications, do let me know how it turns out. Until next time!

Christopher Hanusa's portfolio is available online at

Purchase a copy of this and other 3D printed artwork
at Hanusa ≀ Design on Shapeways.