# 3D printing

## Curve Functions as Voxels

- Intro
- Cartesian Functions
- Spherical Functions
- Curve functions (this post)

Curve functions are defined as a circular co-ordinate function extruded along the length of a parametric curve. Hence, just like Cartesian and spherical functions, we have just two inputs (in this case the angle *θ* and parametric variable *s*) and three outputs (the *x*, *y*, *z* co-ordinates of the surface point):

(*x*, *y*, *z*) = *f*(*s*, *θ*).

This is constructed from the parametric curve

(*x*, *y*, *z*) = *p*(*s*)

and the radius functions

*r* = *q*(*s*, *θ*).

Practically speaking there’s some ambiguity here, because the angle *θ* needs to be given an orientation. One rotational axis is defined by the tangent to the curve, but where does that leave the other two so we know where to start our rotation from? To solve this we use the Frenet-frame of the curve to define an orthogonal co-ordinate space independently at each point along its length. It’s generally a good solution because it tries to maintain a consistent orientation along the length of the curve (although this can fail rather miserably if the curve becomes a straight line, which is a discussion for another time) while also being independently defined at each point on the curve. This last point is important in order to be able to calculate the curve surface points in parallel using the GPU.

This arrangement works well transforming from parameters to surface co-ordinates. But for a voxel map we need the reverse: take a voxel co-ordinate (*x*, *y*, *z*) and establish how far it is from the curve. We have the equation of our curve and can write an equation to capture the distance, but actually solving the equation is a different matter.

What we need is a solution to the following equation.

d*p*(*s*) / d*s* = 0

Calculating the algebraic form of d*p*(*s*) / d*s* isn’t hard (we do it already for the lighting calculations), but solving the equation can be. We’d have to rearrange the equation to establish it in terms of *s*. It could have multiple solutions, as the diagram above shows. In short, given *p*(*s*) is user defined and could be practically any function, we simply can’t rely on being able to solve the equation.

So is there nothing we can do? We can’t find the parametric variables that relate to a specific voxel algebraically. This is a real shame; it means we have to use a different approach entirely. The solution I’ve come up with is to effectively rasterise an approximation of the volume as a series of tetrahedrons that can make up the volume. Sadly this is less accurate and less efficient than the ideal algebraic approach, but at least it seems to work and reasonably well.

A tetrahedron is the 3D primitive counterpart to a triangle (the 2D primitive): we can represent any solid approximately as a series of tetrahedrons just as we can approximate any surface as a series of triangles. We can similarly render the solid into a 3D texture as a series of tetrahedrons in the same way we usually render a surface as triangles onto a 2D texture.

To understand this better, think of the tube defined by our curve as if it were partitioned into a series of cuboids, as in the diagram below. In the diagram the cuboids form a very crude approximation of the cross-section of the tube. This is crude partly because of the low number of cuboids used (in practice we’d usually want more) and partly because of my poor Blender skills (which you’ll just have to try to look past).

Here’s how one of those cuboids might look:

Our voxel space is made up of a series of slices (each slice has the same *z* co-ordinates but with varying *x*-*y* co-ordinates) so we need to be able to slice this cuboid. Performing slicing directly on the cuboid turns out to be hard (for me, at least), so the solution is to convert each cuboid into five tetrahedrons, like this:

This is actually really easy to do, because every vertex of each tetrahedron is just one of the vertices of the cuboid: we just need to pick them out in the right order.

Slicing a tetrahedron is a lot easier. Assuming (for the sake of simplicity) we don’t slice along one of the edges, a plane will always cut the edges of the tetrahedron in either three or four points. If it’s three points, we can just join them to create the triangle that represents the slice through the tetrahedron. If it’s four points, we have a quadrilateral, which we can turn into two triangles. We have to take a bit of care over this to avoid accidentally choosing the concave variant where two of the lines cross:

To avoid this, we check whether two of the lines cross, and if they do, reorder the vertices to fix it. We end up with two triangle which we can then rasterize onto the surface slicing the volume. Computers have been rasterizing triangles since the dawn of time (well, at least 1967, which is as good as).

So we have our sequence of simplification: curve volume, cuboids, tetrahedrons, triangles, voxels. It’s a bit long-winded, but we get there in the end. The inefficiency is compounded by the fact we end up performing this ‘volume rasterization’ repeatedly for each slice of the voxel space. We could do it all in one go, but that would require the entire voxel space to be held in memory while the render is performed. Given the amount of data involved (e.g. a 1024 × 1024 × 1024 voxel space would require a gigabyte of RAM) rendering each slice individually is considerably more memory-efficient, if not time-efficient.

The final results are generally good, although it’s still an approximation as compared to the results for Cartesian and spherical functions. As far as I can tell, there’s no straightforward way to avoid this based on Functy’s current design.

Here’s the images from the original post showing how all three curves - Cartesian, spherical and curve - are voxelised into a set of slices making up their combined volume.

## Spherical Functions as Voxels

- Intro
- Cartesian Functions
- Spherical Functions (this post)
- Curve functions

To convert spherical functions into a voxel volume is a bit, but not much, more coplicated than for Cartesian functions.

In the case of Cartesian functions the key criterion was the height above the *x*-*y* plane. In the case of a spherical function, it’s the distance from the centre point of the function that we’ll need to perform the comparison on.

We define our spherical functions using angles *θ* and *φ* circumlocating this point. The angle *θ* represents the angle around the *z*-axis through the centre point of the object. The angle *φ* represents the elevation above the *x*-*y* plane also through the centre.

Converting to a triangulation involves stepping through *θ* and *φ* with a given step distance and calculating the radius *r* = *f*(*θ*, *φ*) for the function. This radius acts like a spine emerging from the centre at the given angle and with length *r* as defined by the function. These spines are used to define the vertices of the triangles that form the surface of the object.

When voxelating the volume, the question becomes one of whether a given voxel, defined by its (*x*, *y*, *z*) position, is inside or outside the object. We do this by considering the spine that runs between the voxel co-ordinates and the centre of the object, and calculating the *θ* and *φ* values that would apply in this case.

The contrasting nature of the two approaches shows itself here. To create triangles you have to calculate the (*x*, *y*, *z*) position of a triangle vertex from the *θ* and *φ* values. In both cases the *θ* and *φ* are needed to calculate the radius.

Calculating *θ* (rotation) and *φ* (elevation) values isn’t hard with a bit of trigonometry.

*θ* = tan^{-1} (*y* / *x*)

and

*φ* = tan^{-1} (*z* / (*x*^{2} + *y*^{2})^{1/2})

Here’s the code we use to do this.

fZFunc = psFuncData->fZMin + (nSlice * fZStep); for (nX = 0; nX < nResolution; nX++) { fXFunc = psFuncData->fXMin + (nX * fXStep); for (nY = 0; nY < nResolution; nY++) { fYFunc = psFuncData->fYMin + (nY * fYStep); // Calculate the angles fOpposite = (fXFunc - psSphericalData->fXCentre); fAdjacent = (fYFunc - psSphericalData->fYCentre); fElevation = (fZFunc - psSphericalData->fZCentre); fAFunc = atan2 (fOpposite, fAdjacent); fDiag = sqrt ((fAdjacent * fAdjacent) + (fOpposite * fOpposite)); fPFunc = atan2 (fElevation, fDiag); if (psSphericalData->psVariableA) { SetVariable (psSphericalData->psVariableA, fAFunc); } if (psSphericalData->psVariableP) { SetVariable (psSphericalData->psVariableP, fPFunc); } fRFunc = ApproximateOperation (psFuncData->psFunction); fDistance = sqrt ((fAdjacent * fAdjacent) + (fOpposite * fOpposite) + (fElevation * fElevation)); if ((fDistance < fRFunc) && ((psFuncData->boMaterialFill == TRUE) || (fDistance > (fRFunc - psFuncData->fMaterialThickness)))) { ucFill = 255u; } else { ucFill = 0u; } } }

Once we have them, we can plug them into the function that defines our spherical object to calculate the radius at that point *r*‘ = *f*(*θ*, *φ*). We compare this with *r*, the length of the spine from the centre of the object to the position of our voxel.

*r* = (*x*^{2} + *y*^{2} + *z*^{2})^{1/2}.

Notice we’re pretending here the co-ordinates of our object are at the origin. If they aren’t, we need to subtract the centre co-ordinates component-wise from our voxel co-ordinates first.

The comparison is then simply *r* ≤ *r*‘ if we’re inside the object and *r* > *r*‘ if we’re outside.

We perform these calculations and comparison for each voxel in our voxel space. As before, if we’re inside we set the voxel to be filled, otherwise we leave it empty.

That deals with Cartesian and spherical functions, which are the easy ones, because the mapping between voxels and the co-ordinate system is one-to-one. Sadly this isn’t the case for curve functions, which I’ll tackle in the next post.

## Cartesian Functions as Voxels

- Intro
- Cartesian Functions (this post)
- Spherical Functions
- Curve functions

Of the three, the easiest class of functions to turn into voxels was the Cartesian functions. In order to understand the process the key thing I eventually realised was that voxelisation is almost the reverse of triangulation.

For the latter, we step through the *x* and *y* co-ordinates at a given resolution and calculate the consequent *z*-values to define the vertex co-ordinates of the triangle.

In contrast, for the former we step through not just the *x* and *y* co-ordinates, but also the *z* co-ordinate. The question we then have to answer is whether this (*x*, *y*, *z*) point is inside or outside the volume.

In practice, we can’t answer this question without defining the boundaries of the volume. Functy has an option to specify the thickness for a function, but the simplest case is to say anything on or below the function that defines the surface is inside the volume, and anything above is outside.

The *z*-value generated by the function for given *x* and *y* co-ordinates represents the height of the function above the *x*-*y* plane. For a given voxel position (*x*, *y*, *z*) we therefore calculate the height of the function *z*‘ = *f*(*x*, *y*) at the point *x*, *y* and compare this to our actual position (*x*, *y*, *z*). If *z* ≤ *z*‘ then we’re inside the volume, so should fill in the voxel. If *z* > *z*‘ on the other hand, we’re outside the volume, so should leave the voxel empty.

And that’s it. To voxelise the function we check each of the points that make up the voxel space, perform the comparison, and either fill in or leave empty the voxels as we go. There are a lot of voxels to go through because we’re working with cubic dimensions (so even a low resolution of 100 steps per dimension gives us a million separate voxels to consider), but the comparison to perform is pretty straightforward in itself, as is clear from the code that performs the check.

fZSlice = psFuncData->fZMin + (nSlice * fZStep); for (nX = 0; nX < nResolution; nX++) { fXFunc = psFuncData->fXMin + (nX * fXStep); for (nY = 0; nY < nResolution; nY++) { fYFunc = psFuncData->fYMin + (nY * fYStep); if (psCartesianData->psVariableX) { SetVariable (psCartesianData->psVariableX, fXFunc); } if (psCartesianData->psVariableY) { SetVariable (psCartesianData->psVariableY, fYFunc); } fZFunc = ApproximateOperation (psFuncData->psFunction); if ((fZFunc < fZSlice) && ((psFuncData->boMaterialFill == TRUE) || (fZFunc > (fZSlice - psFuncData->fMaterialThickness)))) { ucFill = 255u; } else { ucFill = 0u; } } }

Spherical functions are a little more complicated, but not much. I’ll write about them in the next post.

## Voxel Volumes

One of the main feature additions of the latests version of Functy has been the ability to export as SVX files. Functy could already export in PLY and STL, but both of these are triangle based. They represent the 3D functions as surfaces defined by carefully aligned triangle meshes. Rendering objects using a graphics card also uses the same triangulation process, so exporting as PLY or STL is a very natural extension of the existing rendering.

The SVX format is different though. It stores the models as a voxel image (a voxel being a three dimensional pixel, for those who didn’t grow up through the 90s demo scene). As a result, SVX doesn’t just store the surface, but also the volume of a function.

Turning a triangulated surface into a voxelated volume isn’t necessarily straightforward, but Functy has the advantage of having all its objects originate as purely mathematical forms. In theory, this means voxel rendering them as volumes should be quite easily.

What I found in practice is that for Cartesian functions and spherical functions this is true: they can be turned into voxel volumes in a very natural way. Curve functions are a different story though. In the next few posts I’ll go through each of the processes separately, to give an idea about how the solutions for each of the three function types were coded.

- Intro (this post)
- Cartesian Functions
- Spherical Functions
- Curve functions

## Functy 0.33 now available

This latest version has two main new features. First it allows animations to be exported as a series of PNG frames. Second you can also now export in Simple Voxels (SVX) format. This is a neat new format that represents a model as a volume rather than a surface as would normally be the case with STL or PLY files. Functy is ideally set up for this, since the objects are mathematically defined anyway. SVX files are also ideally suited to 3D printing, which ultimately is working with 3D volumes too.

Even though it’s new, Shapeways already allow upload of models in SVX format, and you can see an example - my first attempt - on the Shapeways site.

There have been some more minor improvements too. For the first time ever, Functy now has an About window. Only a small thing, but long overdue. It’s been on the back-burner for a while, since until recently Ubuntu’s version of Glade crashed when creating dialogues. Ironically I didn’t end up including it in the gtk-builder file anyway. Never mind. There’s also a new Progress window to make long exports bearable. SVX export can take some time, so this became a necessity. But it also helps with exporting of animations too, since these can also take a surprising length of time when every frame has to be generated as a separate file.

Please get yourself a copy and try out the new features. As with every release of Functy since the dawn of time, this is still a beta version, so please bear this in mind and let me know if things go wrong.

As always, the binary and source is available direct from Sourceforge or via the downloads page.

## Projections

The sun was out in Liverpool today, creating crisp and long evening shadows. So it seemed like a great opportunity to take photos of recent 3D printed Functy objects. The full images are rather large, but show the grain of the printing, which I think is rather interesting in itself. Click on the images for the full views.

The original Lissajous is up on deviantArt and Shapeways; the alien egg is also on deviantArt and Shapeways.

## Comparing renders with reality

Shapeways delivered a new batch of 3D prints recently. I was particularly pleased with the Alien Egg print of a spherical function with the formula:

radius = (3*(0.5+(sin(cos(a)+(p*((3+cos((pi*0.3)-1.5))/4))*10)**2)+((6/6.6)*((2+sin(8*a))/3)**4)))*sin(p)+(4.3*(1-(sin(p)**2))),

colour (R, G, B) = (r/8, (3+sin((a)*8))/5, (3+cos(a*8))/5).

The print is actually rather small (just 6 cm diameter) and the ridges of the shape are really quite delicate. In spite of this, the 3D print has come out really very similar to the original design. I guess you might expect it to be pretty similar, given the way it was produced directly from the model! However, if you look really closely at the original you can see the strata through the object created by the printing process. What I’m really impressed with, though, is the colour produced. I’d expected this to be a bit washed out, but in practice it’s a pretty impressive match.

Below is a comparison of (from top to bottom) the Functy render, the 3D print and a render done using Blender Cycles. In case you’re interested and your browser supports APNGs, there’s also a peculiar animated version!

## Lissajous Looping

Following on from my previous post, I thought it’d be interesting to make an animated render of the Lissajous figure. If you have an APNG-capable browser (e.g. Firefox) you can see the result on DeviantArt.

While it’s neat to be able to print static versions of these Lissajous figures, in the future I’m sure it’ll be possible to make the fully moving version as well. Now *that* would be really something!

## Lissajous Loops

Sines and Cosines have been responsible for some of the most elegant mathematical constructs. Lissajous curves are a particularly simple, yet elegant example. Put simply, a Lissajous is a parametric curve where each axis follows a sinusoidal path. By tweaking the amplitude and cycle length for each axis, a myriad of different patterns can be generated, from circles to intricately woven lattices.

The parametric curves in Functy are particularly suitable for generating nice Lissajous curves, and as usual, they can be output for 3D printing. The results of pumping them through a 3D printer, courtesy of Shapeways, can be seen in the photos below, along with a Blender Cycles render of one of the curves.

If you fancy getting really up-close-and-personal with them, you can order your own copies as unusual desk ornaments, from the Shapeways site.

## Network structures

As part of an experimental game project I’ve been trying to use the Functy rendering routines to visualise network structures. At the moment it’s at a very early stage, but has - I think - already generated some interesting results.

The screenshot below shows a network of 60 nodes, each one rendered as a spherical co-ordinate function, joined together using links rendered as curves. I just plucked some simple functions out of the air to see what the results would be like but am hoping to extend it with more interesting shapes as things progress.

The various parts of the network are a little hard to discern with a static image, but when I tried to capture a video the result was a mess of fuzzy artefacts (I think there must be something going wrong with my screen capture software), so I gave up on that.

The next step, after neatening up the code, is to arrange better animation of the nodes and links, with dynamic movement based on things like the forces between the nodes. I’m hoping this will produce some really nice effects, and if anything comes of it I’ll put a bit more effort into getting a successful video capture.