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

**Sorry, the comment form is now closed.**