- 0.3.5
- 0.23
- 0.24
- 0.32
- 3D
- 3D print
- 3D printing
- 033
- ambient occlusion
- Anderson
- animate
- animated
- animation
- announcement
- apt-get
- ARM
- binary
- blender
- blender cycles
- blog
- blur
- calculators
- Cartesian
- co-ordinates
- Coding
- comments
- comparison
- contribute
- contributors
- cosine
- curvature
- curve
- Curves
- cylindrical
- Dandelion
- deb
- Debian
- demo
- depth
- description
- Development
- deviantart
- download
- downtime
- Enzyme
- failure
- focus
- framebuffer
- Frenet Frame
- freshmeat
- function
- functions
- Functy
- gnome
- gnomefiles
- Goldman
- graphics
- GTK+2
- hot picks
- HTML5
- image
- impicit representation
- implicit curve
- install
- island
- jewellery
- jewelry
- knot
- landscape
- Launchpad
- lava lamp
- library
- link
- linux
- linux format
- linuxlinks
- lissajous
- literuture
- lxf
- magazine
- management
- maths
- MIT licence
- Model
- negative
- network
- network visualisation
- O3D
- object oriented
- OpenGL
- open source
- openvdb
- outage
- packaging
- paper
- parallelisation
- parametric
- parametric curve
- parametric curvec
- parametric surface
- ppa
- publication
- Raspberry Pi
- Release
- Render
- rendering
- repository
- research
- resolution independent
- review
- rings
- rpm
- scenary
- Screenshot
- Screenshots
- Sederberg
- Shader
- shaders
- shades
- shadows
- Shapeways
- sine
- slice
- slicing
- softpedia
- software
- source
- sourceforge
- spam
- sphere
- spherical
- string
- subtraction
- surface
- svn
- svx
- Symbolic
- Tony Ralano
- torsion
- Trac
- tubes
- tuxradar
- ubuntu
- update
- user interface
- version 0.1
- version 0.2
- version 0.21
- version 0.22
- vertex shader
- Video
- visualisation
- visualization
- volume
- voxel
- voxels
- Web
- website
- windows

# Posts tagged with 'voxel'

## 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.