Generate an STL 3D model with a surface described by an equation

I recently needed to generate a 3D model for printing based on an equation and couldn't find software to do the job.  In the end the easiest thing to do was to write a quick Python script called "stl-surface.py" to do the job.

At first the task seemed daunting, but generating the STL file is rather easy with the use of the numpy-stl library. The equation for the surface is calculated over a structured square/rectangular grid and each little cell of the grid is split in two to form triangles. You use many many of these small triangular faces to construct the model. All you have to do is to generate the coordinates for the 3 vertices that make up each triangle and numpy-stl does the rest. The bottom and sides are flat which makes things rather easy, and the top is defined by the equation.

When constructing the model it's important to list the coordinates of each face in a counter clockwise direction when looking at the model from the outside. This allows numpy-stl to later calculate things like volume of the model and centre of gravity and also test if the surface is closed (although they use a slightly buggy way to test this).

I've shown a few examples of the generated models below.

10 - 2*(1 - math.cos(2 * x * math.pi / 20)) * (1 - math.cos(2 * y * math.pi / 20))

Dimpled 3D model

10 + 2 * math.cos(math.sqrt((x - 50)**2 + (y - 50)**2)/2)

Rippled 3D model

Here you can see how the triangles are assembled to construct the model.

The triangular faces that make up the model

The triangles that form the edges of the model

Edge of the model

From the top you can see the grid points where the surface function is evaluated.

Top of the model

I hope this code can help someone else. I may not be the easiest thing to use and if you need help getting started get in touch.

Get the code!

Calculating Hilbert Curve Coordinates

You may have noticed from some of my previous blog posts that I've fallen down a rabbit hole learning about Gray codes and Hilbert curves.  They weren't just random topics.  I wrote those articles in order to understand the topic of this one, how to easily convert a Hilbert Index to Hilbert Coordinates via a method described in a paper by John Skilling.

Skilling, J. (2004), "Programming the Hilbert Curve", AIP Conference Proceedings.

Regrettably I don't yet understand why it works to a level that I'm satisfied with.  The paper describes the process and gives a summary of previous work by Butz and Lawder but doesn't really give a clear explanation of why the transform works.  It may be perfectly clear to someone who works in the field, but I don't get it... yet.

What I am capable of doing though is summarising the information here with some graphics to explain the concept to the next person that gets nerd sniped by this.  I'll start by recommending my last few blog posts as they will give you an idea as to my train of thought.

Entry and Exit Points for Space Filling Paths on a Grid

Hilbert Curve Generation With Lookup Tables

Converting Binary to Gray Code with XOR

A basic concept to understand is how binary numbers are treated as coordinates.  In the animation below, the index is 6 bits long and counts up from 000000 to 111111.  This can be used to cover every coordinate on an 8 x 8 two dimensional grid if every second bit belongs to the coordinates of one dimension while the others belong to the other dimension.  In this case the red bits belong to the x axis while the blue bits belong to the y axis.  Everything I'll describe in this post is for two dimensions, but it works for any number of dimensions.  For example an 8x8x8 grid could be covered with a 9 bit index where bits 0, 3, 6 belong to one axis, 1, 4, 7 to another, and 2, 5 ,8 to another.

Binary indexing in 2 dimensions

In the binary indexing animation you can see that there are diagonal lines.  These occur because more than one bit of the index can change at once affecting both dimensions.

The first step in the Skilling Transform is to take the binary reflected Gray code of the index.  This means that only one bit at a time can change.  Therefore the value of only one axis can change at a time leaving only vertical and horizontal movements.

Gray code indexing in 2 dimensions

The addressing scheme detailed above is formally described in this section of the paper.  I put this here just so we're all on the same page about the variable names.

$$$p$$$ is the number of bits in each axis
$$$n$$$ is the number of dimensions

Addressing description

Now onto the actual transform.  It's minimal and easy to compute.  It's not a recursive algorithm and takes a bit string that specifies the Hilbert index as an input and performs a series of swaps and inversions to return a bit string of the same length that describes the Hilbert coordinates.

Skilling transform

To see if I could understand the process better I wanted to visualise the transform.  The animation below show the transformation from Gray code indexing to the Hilbert curve.  You can see how the algorithm starts the transformation on the smaller features first and works up to larger features at each step.

Gray code to Hilbert curve via the Skilling transform in 2 dimensions.

Hopefully the worked example below will help people understand the process a little better.  The index is converted to its Gray code and the algorithm processes bits in a reverse order.  The grey boxes refer to the bit $$$r$$$ shown in the algorithm above.  The bits that are classified as low bits at each stage are underlined.

Worked example

I can't really offer much of a conclusion here.  The understanding I've gleaned so far is minimal and superficial at best.  I 'll leave this one for a while and come back to it in the future when I have more time.

I have a love hate relationship with problems like this.  I find them fascinating, but at the same time they remind me how much I miss formal study, and I get bummed out.  Kinda sucks having no one to bounce ideas off either.  Having said that though, I've finally generated some animations that I'm proud of.  You can find the undocumented scratch code here.

If you come across this page and have a better understanding of why this works I'd be glad to hear from you.

Calculating Dihedral Angles

Just a quick one today.  Up until recently I had never heard of the term dihedral angle.  It sounds complicated but it's something you already know.  If you have two surfaces in spaces that meet, they will form a line. The angle between the two surfaces at this intersection is called the dihedral angle.

I had recently been wondering about how to calculate the angle between two surfaces and decided to work through the problem myself.  First of all you need to describe the surfaces and I thought that the best way to do this was by supplying vectors normal to each surface, v1 and v2.

Angle between two surfaces

By simplifying the geometry and extending the normal vectors until they intersect you can start to see how to solve the problem.

Extended surface normal vectors

You can also derive the following formula from the geometry.

Relation between alpha and beta

From this you can see that the cosine of beta is equal to the negative cosine of alpha.

Simplified cosine term

By taking the dot product of the vectors v1 and v2 you can quickly derive an expression for the angle alpha.

Expression for alpha

Great!  Problem solved.  Or is it?  The inverse cosine function only returns values between 0 and 180. This means you will never get an answer greater than 180 degrees, meaning you can't have convex angles.  The image below demonstrates why.  Given the vectors v1 and v2 to define the surfaces in the first part of the image below you would quickly run into trouble.  You can see in the lower part of the image that the same two vectors can be used to define a different geometry with a concave angle.

Ambiguity between vectors

As the equation for alpha above doesn't change if the vectors v1 and v2 are swapped the two geometries above will give the same angle.  There are ways to work around this problem, you just need to be aware of the geometry you're working with.  I'm still having trouble getting my head around this.  The order the angles are given is important and would allow you to work out if the angle is concave or convex, but when generalising this to 3 dimensions you have trouble depending on what angle you look at joint.  I'm falling asleep typing this with pictures of vectors dancing around my head.  Maybe a fresh set of eyes will help.