December 2014 Archives
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.
dp(s) / ds = 0
Calculating the algebraic form of dp(s) / ds 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.