Principles and Practice

A sample implementation is provided for two contouring algorithms: the classic
*Marching Cubes* and the newer *Dual Contouring*.

Download the zip file containing this page and the sample implementation.

Download a zip file containing pre-compiled binaries.

Consider the field function for a sphere. It is
`f(x,y,z) = x ^{2} + y^{2} + z^{2} - R^{2}`
where

This representation of volumes in three-dimensional space is very compact and precise.
To explicitly describe the sphere from the example above would require specifying a
list of points and polygons that make up an approximation of the contour of a sphere.
An implicit representation also allows for precise CSG (*constructive solid geometry*)
operations to occur in the domain of the function. To perform such an operation on an
explicitly-represented contour would first require figuring out what the actual volume is.

Explicitly-represented contours can be rendered much faster than implicit volumes, as the polygons are already known and can be easily drawn on the screen, often by sending them to a graphics card that is designed to render polygons. To render an implicit surface this way, however, it is first required to compute the polygonal approximation of its contour.

An implicit surface is defined by a derivated of the class `Isosurface`
(in files `threed/isosurface.*`). This class is the interface for implementing
the field function `d = f(x,y,z)` discussed above, by providing a pure virtual
function `Isosurface::fDensity` that computes the field function. To gain better
performance, `fDensity` can compute any number of points, starting at a point
specified by three parameters `(x0,y0,z0)` and progressing towards the +z axis by
a delta specified by the parameter `dz`.

As an example, an implementation of the field function of a sphere using `Isosurface`
would look like this:

class SphereIsosurface : public Isosurface { public: SphereIsosurface(float rad) : _radius(rad) { Vector center; // (0,0,0) Vector v1 = center - _radius; Vector v2 = center + _radius; Isosurface::addBoundingBox(BoundingBox(v1, v2)); } void SphereIsosurface::fDensity( float x0, float y0, float z0, float dz, int num_points, float *densities) { for (int i = 0; i < num_points; ++i) { float sqr_dist = x0 * x0 + y0 * y0 + z0 * z0; float sqr_rad = _radius * _radius; float d = sqr_dist - sqr_rad; densities[i] = d; z0 += dz; } } protected: float _radius; };

Note how the constructor of `SphereIsosurface` computes an axis-aligned
bounding box containing the sphere and reports it using `Isosurface::addBoundingBox`.
Generally speaking, it is not easy to figure out the bounding box just by looking at
the field function. A few algorithms have been proposed that address this problem,
but they are not perfect: for example, they may not work correctly when the
field function describes volumes that are not connected. (But see [7],
[8] for more information.) In many cases `Isosurface` is
derived to implement solid primitives, such as a sphere, a box, or a torus, whose bouding
box can be easily computed. Therefore I find it reasonable to expect the implementation
of an `Isosurface` to report its own bounding box.

The implementation of `Isosurface` also permits transformation of the volume.
The class `Matrix` (in files `threed/matrix.*`) implements an 4x4
transformation matrix. Such a matrix can cause any combination of translation,
rotation, scaling and shearing operations to occur on a point (or an object made
up of a collection of points) in 3D space. The `Matrix` class is also
compatible with the way OpenGL deals with matrices. The class `Transform`
(in files `threed/transform.*`), that inherits `Matrix` without adding
any new data members, has a function `Transform::glMultMatrix` which is
simply a wrapper around the OpenGL function `glMultMatrixf`.

A transformation matrix can be applied to the implicit surface using
`Isosurface::setTransform()`. Now consider a sphere of radius 1, that has been
translated on the x axis and is now centered at (‑5,0,0). Evaluating the field
function for the point at (0,0,0) should return a result indicating the point is outside
the volume. In this case `SphereIsosurface::fDensity()` as presented above will
evaluate to an incorrect result, as it does not take the translation into account. But
the definition can be easily revised to take any kind of transformation into account:

void SphereIsosurface::fDensity( float x0, float y0, float z0, float dz, int num_points, float *densities) { for (int i = 0; i < num_points; ++i) {float xt, yt, zt; inverted_matrix.transform(x0, y0, z0, &xt, &yt, &zt);float sqr_dist = xt * xt + yt * yt + zt * zt; float sqr_rad = _radius * _radius; float d = sqr_dist - sqr_rad; densities[i] = d; z0 += dz; } }

The revised version first transforms the specified point given by `(x0,y0,z0)`
by the inverted matrix of the transformation applied to the implicit surface. When
no transformation has been applied, the transformation matrix is the identity matrix,
and so is its inverse, in effect resulting in `xt=x0, yt=y0, zt=z0`. But to
revisit the last example where a sphere is translated by (‑5,0,0), the inverted
matrix would translate the point by (+5,0,0), resulting in `xt=(x0+5), yt=y0, zt=z0`.
Thus evaluating the field function at (0,0,0) would correctly return a result indicating
the point is outside the volume, while evaluating the point at (-5,0,0) would correctly
return a result indicating the point is inside the volume.

The following illustrations show a grid of size 5x5x3, and an outline of a single grid cube.

This is only an overview of Marching Cubes and far from a comprehensive description. See [3] for a much more in-depth discussion, including a sample implementation of the algorithm. My own implemention is described in the following section.

The algorithm is implemented in function `IsoMesher_MC::createMesh()`, where first
an axis-aligned bounding box is computed. The bounding box is computed in the
`Isosurface` itself by the function `Isosurface::getBoundingBox()` which
takes into account any transformations applied on the implicit surface, and also computes
the inverted transformation matrix for use by `Isosurface::fDensity()`. Once the
bounding box is known, the number of cubes it contains in each of the three axes x, y and
z is computed by a simple division by the size of the cube:

xsize = ( boundingBox.vmax().x() - boundingBox.vmin().x() ) / cubeSize.x(); ysize = ( boundingBox.vmax().y() - boundingBox.vmin().y() ) / cubeSize.y(); zsize = ( boundingBox.vmax().z() - boundingBox.vmin().z() ) / cubeSize.z();

`IsoMesher_MC::createMesh()` examines the grid in rows. A row is a slice of grid
points (not voxels) all having the same y coordinate, and is made up of two two-dimensional
arrays: `points[xsize][zsize]` which holds the (x,y,z) location of the grid point,
and `densities[xsize][zsize]` which holds the result of invoking
`Isosurface::fDensity()` on each of the points. Two rows of grid points describe
a single row of grid cubes, as each cube is bound by eight grid points.

Computing a row (as performed by `IsoMesher_MC::computeRow()`) means initializing
the `points[][]` array and evaluating the field function for each point. To
accomplish this, `computeRow()` has to know the bottom left near coordinate
of the row: this is essentially the minimum coordinate of the bounding box,
incremented along the y axis as the algorithm progresses throughout the rows.

void IsoMesher_MC::computeRow(const Vector &vRow) { float x0 = vRow.x(); // get the bottom left near float y0 = vRow.y(); // coordinate for the float z0 = vRow.z(); // current row of points float dx = cubeSize.x(); float dz = cubeSize.z(); for (x = 0; x < xsize; ++x) { // initialize the points structure for (z = 0; z < zsize; ++z) points[x][z] = Vector(x0 + x * dx, y0, z0 + z * dz); // evaluate the points (x,y0,zmin) through (x,y0,zmax) isosurface->fDensity(x0 + x * dx, y0, z0, dz, zsize, densities[x]); } }

The function `createMesh()` uses `computeRow()` to compute the first grid
row at the bottom of the grid, then starts a loop to process each of the subsequent rows
from the next-to-bottom row towards the topmost row. For each such subsequent row,
it invokes `computeRow()` to compute it, then calls `IsoMesher_MC::marchCubes()`
to examine the row of grid cubes bound by the two adjacent rows of grid points.

`IsoMesher_MC::marchCubes()` processes each cube in an x-major, z-minor order.
That is, it marches over the x axis of the grid from left to right, and for each x
coordinate it marches over the z axis from near to far. Each cube is processed separately:
no information that came to be known by processing one cube is needed for processing
any other cube. This means than an implementation of the Marching Cubes algorithm can be
easily parallelized for multiprocessor environments.

For each cube, the results of the evaluation of the field function for the eight grid points at the corners of the cube are examined. An 8-bit mask is initialized to zero, and each corner is associated with a specific bit position. If the result of the field function for a corner is negative, that is, if the grid point is inside the volume, then the associated bit in the mask is set to 1. This mask describes the relationship of the grid cube with the contour of the implicit surface. The mask is all zeros (0x00) or all ones (0xFF) if the cube is entirely outside or entirely inside the volume, respectively. Any other mask value indicates that some of the corners are inside the volume while others are outside it, and this means the cube intersects the contour in some way and polygons need to be generated.

A cube made of eight corners has 12-edges, any of which can interest the contour. A special
table, the *edge table*, which has 256 entries and is indexed by the bitmask,
indicates, using a 12-bit mask, exactly which edges interesect the contour. As each corner
was associated with a bit position in the 8-bit mask, so each edge is associated with
a bit position in the 12-bit mask. Therefore the bit positions associated with each corner
and each edge must be consistent with the values in the edge table. My implementation uses
an order similar (but not identical) to the one in [3].

x,y,z denote the bottom left near corner of the cube being processed. corner 0 is at (x, y, z ) corner 1 is at (x+1, y, z ) corner 2 is at (x+1, y+1, z ) corner 3 is at (x, y+1, z ) corner 4 is at (x, y, z+1) corner 5 is at (x+1 y, z+1) corner 6 is at (x+1, y+1, z+1) corner 7 is at (x, y+1, z+1) edge 0 is between corners 0 and 1 edge 1 is between corners 1 and 2 edge 2 is between corners 2 and 3 edge 3 is between corners 3 and 0 edge 4 is between corners 4 and 5 edge 5 is between corners 5 and 6 edge 6 is between corners 6 and 7 edge 7 is between corners 7 and 4 edge 8 is between corners 0 and 4 edge 9 is between corners 1 and 5 edge 10 is between corners 2 and 6 edge 11 is between corners 3 and 7 |

The edge table indicates which edges intersect the contour, but the exact intersection point can lie anywhere on that edge. For instance, if edge table indicates an intersection along edge 2, the exact intersection point lies somewhere on the horizontal line segment that has corners 2 and 3 as its endpoints. There are numerous ways to determine the exact intersection point. The most simplistic approach is to linearly interpolate the position using the results of the field function. (Again, see [3].) While very easy to implement, linear interpolation assumes the field function is continuous. While this may be true for a sphere, it is hardly the general case. Consider the field function for a box:

f(x,y,z) = max(x^{2}- Rx^{2}, max(y^{2}- Ry^{2}, z^{2}- Rz^{2}))

Iterative evaluation of this function for points that are nearby each other generally returns
a smooth progression of values, except at the corners of the box where the use of `max`
causes an abrupt switch from one axis to another. Another case is a volume that is the
result of CSG operations on other volumes, which means it is composed of any number of
independant field functions. In these cases, employing linearly interpolation to compute
the exact intersection point will yield incorrect results.

A far better choice for finding the exact intersection point is the midpoint (sometimes called bisection) algorithm for root solving. [9] According to the midpoint algorithm, the field function is evaluated for the middle point between the two endpoints, and the result is examined. If the result is zero, the exact intersection point was found. Otherwise, the midpoint replaces the endpoint for which the evaluation of the field function has the same sign, and the algorithm re-iterates. In this way the midpoint algorithm iteratively halves the line segment until it finds the line segment whose middle point intersects the contour.

Another choice is the false position algorithm for root solving [9]. It is a combination of both linear interpolation and the midpoint algorithm. Unlike simple linear interpolation, false position finds the exact intersection point even if the field function is not continuous. In terms of performance--the number of iterations made before finding the intersection point--it is no worse than the midpoint algorithm, and usually better. Thus it makes sense to use it as the algorithm of choice for finding the exact intersection point.

Both midpoint and false position algorithms are implemented by class `IsoMesher`
(in files `threed/isomesher.*`) in the functions `IsoMesher::intersect_xaxis()`,
`IsoMesher::intersect_yaxis()` and `IsoMesher::intersect_zaxis()`. The
`_xaxis` version should be used for finding the intersection point on edges 0, 2, 4
and 6 that have fixed y,z coordinates. The `_yaxis` version should be used on edges
1, 3, 5, 7 that have fixed x,z coordinates. And the `_zaxis` version should be used
on edges 8, 9, 10, 11 that have fixed x,y coordinates. The implementation is of the midpoint
algorithm if the preprocessor symbol `BISECTION` is defined, or of the false position
algorithm if the symbol `FALSE_POSITION` is defined.

An different approach to finding exact intersection points is presented in the Extending Marching Cubes algorithm [12]. EMC modifies the field function to return not one but three values, where each value measures the distance of the point from the contour of the implicit surface on a distinct axis. It employs linear interpolation to compute the x coordinate of the exact intersection point using only the x component of the result of the field function. Similarly it computes the y coordinate using only the y component, and the z coordinate using only the z component.

Once the intersection points are known, they can be connected to form polygons that
approximate the contour. The *triangle table*, which also has 256 entries and is
indexed by the 8-bit mask similarly to how the edge table is indexed, specifies which
intersection points to connect, and in what order, so that triangles may be formed.

Finding exact intersections points and connecting them to form triangles is accomplished
by the function `IsoMesher_MC::generateFaces()`, which is invoked by
`IsoMesher_MC::marchCubes()` for each cube that was found to intersect the
contour in some way. By the time all cubes have been marched through, a complete
polygonal approximation of the contour has been generated, and the process is
complete.

The example program `demo1` (in directory `demos/`) shows how to derive
the `Isosurface` class to define implicit surfaces for a sphere and a box, and
how to create a polygonal approximation of an implicit surface using `IsoMesher_MC`.
This is a screenshot of `demo1` taken just after invoking it.

When the demo is running, the keyboard and mouse may used to examine the objects in the scene. The gray cursor keys move the viewer around, as well as the gray plus and minus keys. The numeric pad keys 4 and 6 rotate the viewer around the y axis. Clicking the mouse on the object selects it for manipulation. The selected object may rotated using the keys a,d,w,x,q,e, and may be scaled larger or smaller using the keys 1 and 2.

Another problem of Marching Cubes--and generally of any implicit surface approximation method that uses a uniform sampling resolution--is that details finer than the resolution of the grid are missed. For instance, consider a box with a small hole carved in it. The polygonal approximation will not exhibit this hole if the hole falls entirely within one grid cube. Approximation errors of this sort may typically be corrected by increasing the resolution of the sampling grid, at the cost of generating far more triangles. This cost can be lowered by using an adaptively-refined grid (or an octree, dividing 3D space into a number of grids of varying resolution) in place of a uniform grid. Some work in this area is presented in [5], [6], [10], [11].

But as noted in [12], refining the sampling resolution (wheather uniformly or adaptively) cannot solve the problem of missing sharp features: "In principle, we could reduce the approximation error of the surface S* extracted from the discretizied field f* by excessively refining the grid cells in the vicinity of the feature. However, the normals of the extracted surface S* will never converge to the normals of S."

The Extended Marching Cubes algorithm presented in [12] attempts to solve two problems of the classic algorithm. The first is that of computing exact intersection points using linear interpolation, as described in the previous section. The second is that of approximating sharp edges and corners, by examining not only the exact intersection points, but also the normal vector at each point, and using the combination of both to estimate the location of the sharp corner or edge.

The normal vector is a vector perpendicular to some plane. In the context of an implicit surface, the normal vector of some point (x,y,z) in the volume can be approximated as the gradient of the field function in each of the three axes x, y, z.

nx = f(x + E, y, z) - f(x, y, z) ny = f(x, y + E, z) - f(x, y, z) nz = f(x, y, z + E) - f(x, y, z)

where E is some small epsilon (eg, 0.001). A plane is defined by a point on the plane and the normal to that point:

n . p = 0(The dot

Consider the case of a grid cube containing the corner point of an implicit box surface. There are three intersection points of the grid cube with the contour. The classic Marching Cubes algorithm connects these points to form the triangle approximating the corner. When the normal vector for each point is known, three planes can be computed. These three planes together form the corner point, and indeed the intersection of three planes is a point.

Given three points P1, P2, P3 and their respective normal vectors N1, N2, N3, solving the following system of linear equations finds the intersection point x.

N1 . (x - P1) = 0 N2 . (x - P2) = 0 N3 . (x - P3) = 0This system can be easily solved in matrix form as

[ P1 P1 P1 ] [ N1 . P1 ] [ P2 P2 P2 ] [ N2 . P2 ] [ P3 P3 P3 ] [ N3 . P3 ] matrix A vector B

Now instead of a corner, consider the case of a sharp edge of an implicit box surface. In this case the grid cube will intersect the contour at four points, but these define only two planes. These are the two planes that are on either side of the sharp edge. The intersection of two planes is a line, not a point. In this case the linear system is underdetermined and cannot be solved.

To get around this problem, EMC approximates an exact intersection point using the quadric error function

E[x] = x - Ni . Pi

This is the least squares solution to the problem. Put simply, the solution of the
quadric error function `E[x]` is the point that is at the same distance from all
planes involved. The QEF is best solved using the
singular value decomposition (also [13]) of the matrix A
from above into the matrices U, D and V, and solving

U D V^T x = b

void SphereIsosurface::fNormal( float x0, float y0, float z0, float *nx, float *ny, float *nz) { float xt, yt, zt; inverted_matrix.transform(x0, y0, z0, &xt, &yt, &zt); float sqr_dist = xt * xt + yt * yt + zt * zt; float sqr_rad = _radius * _radius; float d = sqr_dist - sqr_rad; sqr_dist = (xt + 0.001f) * (xt + 0.001f) + yt * yt + zt * zt; float nx1 = sqr_dist - sqr_rad; sqr_dist = xt * xt + (yt + 0.001f) * (yt + 0.001f) + zt * zt; float ny1 = sqr_dist - sqr_rad; sqr_dist = xt * xt + yt * yt + (zt + 0.001f) * (zt + 0.001f); float nz1 = sqr_dist - sqr_rad; inverted_matrix.transformNormal(nx1, ny1, nz1, nx, ny, nz); normalize_vector(nx, ny, nz); }Conceptually, the computation of a normal vector involves evaluating the field function four times. Therefore the polygonizer can compute the normal vector without invoking

To compute the quadric error function, the polygonizer populates a matrix and a vector
similar to the way for solving a linear system of equations as `Ax = B`. It then
invokes the static member function `QEF::evaluate()` to solve the system using
singular value decomposition.

double matrix_A[12][3]; double vector_B[12]; int rows = 0; for each intersection point { // px,py,pz is the intersection point // nx,ny,nz is the normal vector at that point matrix_A[rows][0] = nx; matrix_A[rows][1] = ny; matrix_A[rows][2] = nz; // compute dot product vector_B[rows] = (double)(px * nx + py * ny + pz * nz); ++rows; } // solve Ax=B using singular value decomposition Vector newPoint; QEF::evaluate(matrix_A, vector_B, rows, &newPoint);

The class `QEF` (in files `threed/qef.*`) implements the entire process
needed to solve a QEF as a black box usable through `QEF::evalute()`. The
source is documented by there is no practical need to describe the mathematics
involved. (The curious reader may refer to [13] for detailed
information.)

The second drawback has to do with how EMC actually approximates the contour in grid
cubes that contain sharp features. The EMC solves the QEF to find point `P`. It
then creates a triangle fan connecting each pair of exact intersection points to the point
`P`. This works well for corners, but for sharp edges the point `P` is
roughly the midpoint of the segment of the edge within the grid cube. (See the
illustration below.) The result of this is what appears to be a series of small pyramids
along the sharp edges. Therefore EMC requires an additional step of edge flipping, to
complete the approximation of the sharp feature.

The Dual Contouring algorithm presented in [14],
[15] takes a different approach to using the QEF. It solves the
QEF for every grid cube that intersects the contour of the implicit surface, wheather
that cube conains sharp features or not. That new point `P` is assocated with
the grid cubes. It approximates the contour by quadrilaterals connecting the `P`
points of four adjacent grid cubes sharing the grid cube edge where an intersection
with the contour occurs.

The grid is examined in a bottom-to-top, left-to-right, near-to-far sweeping motion. Therefore for each grid cube, only edges 5, 6, and 10 need to be examined to decide which four adjacent grid cubes participate in the generation of the quad polygon. (The number of edges is the same as with Marching Cubes.) If the reason for only these three edges being important is not immediately clear, the following example should clarify the issue.

If the bottom, left, near grid cube, call it cube A at position (x,y,z), intersects the
contour, it can only intersect it on edges 5, 6, or 10. All other edges would definately
be on the outside of the implicit surface. If the intersetion is on edge 6, then the
three other grid cubes sharing the edge are:

- cube B, immediately to the far side of cube A, at position (x,y,z+1)

- cube C, immediately to the above of cube A, at position (x,y+1,z)

- cube D, immediately to the above and far side of cube A, at position (x,y+1,z+1)

A quad is generated connecting the new points of each of the four adjacent cubes
A, B, C, D.

As the algorithm sweeps near-to-far along the grid, the next grid cube to be processed is cube B. It intersects the contour at least on its edge 2, since the grid cube to its near side (that is, cube A) intersects the contour on its edge 6. Cube A's edge 6 is the same as cube B's edge 2. But the quad connecting these two (and two other) grid cubes has already been generated during the processing of cube A.

This logic works the same way if cube A was interesecting the contour on edges 5 or 10. Therefore for any grid cube, examining any edge other than 5, 6, or 10, means re-doing work that was already done for a previous grid cube.

The first difference between the two polygonizers is that `IsoMesher_DC` works
with three rows of grid points instead of only two as in `IsoMesher_MC`. This
is due to the fact that Dual Contouring connects four adjacent grid cubes into one
quad, and if the intersection is on edges 6 or 10, two of the cubes will be from the
row of grid cubes (not grid points) immediately above the row of grid cubes being
processed.

`IsoMesher_DC` also has a notion of a cube structure, where `IsoMesher_MC`
only knows about rows of grid points. For a row of `xsize` by `zsize`
points, there are `(xsize - 1)` by `(zsize - 1)` cubes. Each cube
structure keeps its 8-bit mask (same 8-bit mask as in Marching Cubes) and the coordinates
computed for its new point.

The process of evaluating the field function (in `IsoMesher_DC::computePoints()`)
is identical to that in `IsoMesher_MC`. `IsoMesher_DC::computeCubes()`
is similar to `IsoMesher_MC::marchCubes()` in that is examines the grid cube
to see if it intersects the contour, but unlike the Marching Cubes version, it does not
generate any polygons. Instead, it calls `IsoMesher_DC::generateVertex()`
to compute the new point, and stores it along with the 8-bit mask in the cube structure.

The process of computing the new point in `IsoMesher_DC::generateVertex()`
is simply computing the exact intersection points as was done in Marching Cubes,
and storing these points along with their normals in the matrix A and vector B described
above for use with the QEF solver. It then invokes `QEF::evaluate()` to compute
the new point. One thing to note is that it computes the *mass point*, the average
point of all intersection points. The point stored in the vector is really the
intersection point minus the mass point. And the mass point is added back into the
result of the QEF to produce the new point that is stored in the cube structure.
The use of the mass point is needed to guarantee that the computed point will lie
within the grid cube even when the system of linear equations is underdetermined.

Once three rows of grid points and two rows of grid cubes have been computed, the
next step is to generate quads, accomplished by `IsoMesher_DC::generateQuads()`.
This function examines `(xsize - 2)` by `(zsize - 2)` grid cubes of
the row of grid cubes. Note that there is no need to examine the very last grid cube
because it does not have any neighboring cubes that share its edges 5, 6 or 10.
`IsoMesher_DC::generateQuads()` examines edges 5, 6, and 10 for each cube
to see which three adjacent grid cubes to use for the generation of the quad.

Edges 5, 6 and 10 are used to pick the three adjacent grid cubes, but this is not
enough to tell the ordering of vertices in the quad. The ordering of vertices is
reversed in these three cases:

- for edge 5, if corner 5 is inside the volume (as opposed to corner 6 being inside)

- for edge 6, if corner 7 is inside the volume (as opposed to corner 6 being inside)

- for edge 10, if corner 6 is inside the volume (as oppposed to corner 2 being inside)

This is only important when the graphics subsystem performs *backface culling* and
skips the drawing of polygons that are not facing the viewer. In this case examining the
corners guarantees a closed polygonal approximation of the implicit surface.

The example program `demo2` (in directory `demos/`) renders the same
image as `demo1` using `IsoMesher_DC`. This is a screenshot of
`demo2` taken just after invoking it, showing how well Dual Contouring performs
for curved surfaces as well as those having sharp features.

In CSG there are generally three ways to combine solid object. The *union*
operation adds or merges any number of solid objects. The *intersection* operation
keeps only those parts of the objects that are overlapping. The *difference*
operation subtracts from the first object the union of all other objects.

The class `CsgIsosurface` (in files `threed/csgisosurface.*`) derives
`Isosurface` to implement CSG operations. An instance of CsgIsosurface performs
a single CSG operation on its children. The specific operation is set with a call
to `CsgIsosurface::setCsgMode()` with a parameter of `CSG_UNION`,
`CSG_INTERSECTION`, or `CSG_DIFFERENCE`. Children are added one at a
time using `CsgIsosurface::addChild()`. The children are any derivative of
`Isosurface` such as the `SphereIsosurface` or `BoxIsosurface`
or even another instance of `CsgIsosurface`.

`CsgIsosurface::fDensity()` invokes `Isosurface::fDensity()` on each
of the children of the `CsgIsosurface`. But it has to return only one
value from this set of values. (Actually it has to store that value in the
`densities` parameter.) For a union operation, if a point is inside at
least one implicit surface, it considered inside the combined volume. Therefore
the minimal value from set of values is selected. For an intersection operation, all
points must be inside, so the maximal value is selected. For a difference operation,
which is basically an intersection of the first object with the inverse of all
other objects, all but the first value in the set are negated, and the maximal value
of this new set is selected.

The example program `demo3` (in directory `demos/`) renders a sphere
subtracted from a box using `IsoMesher_DC`. This is a screenshot of
`demo3` taken just after invoking it.

Notice how the different materials merge at the edges of where the CSG operation occurs. The polygonizer does not take materials into consideration, and this can lead to generation of polygons where the color drastically changes between vertices. To separate different colors (or materials) into distinct colors, it is necessary for the polygonizer to treat material changes in the implicit surface the same way as it does the change between inside and outside.

To compile and run the samples, you need to have cygwin installed. Make sure your installation of cygwin contains the OpenGL headers and libraries. You additionally need the SDL library.

You may not reprint this article in any form, electronic or otherwise, without expressed consent from the author.

The source code of the sample implementation--that is, all files in the
`threed` and `demos` directories contained in the zip archive--is hereby
put in the public domain. Do with it what you will.

THE AUTHOR MAKES NO REPRESENTATION ABOUT THE SUITABILITY OF THIS SOURCE CODE FOR ANY PURPOSE. IT IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF ANY KIND. THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOURCE CODE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOURCE CODE.

Have a nice day.

[2] Bloomenthal J.,
*Polygonization of Implicit Surfaces*,
Computer Aided Geometric Design 5, pp 341-356, 1988

[3] Burke P.,
*Polygonising a scalar field*
(http://astronomy.swin.edu.au/~pbourke/modelling/polygonise/),
May 1997

[5] Hall, M. and Warren, J.,
*Adaptive Polygonalization of Implicitly Defined Surfaces*,
IEEE Computer Graphics and Applications, pp 33-42, 1990

[6] Cignoni, P. et al,
*Reconstruction of Topologically Correct and Adaptive Trilinear Isosurfaces*,
Computer & Graphics, Vol 24 no 3, pp 399-418, 2000

[7] Witkin A. and Heckbert, P.,
*Using Particles to Sample and Control Implicit Surfaces*,
Computer Graphics (SIGGRAPH 94 Proceedings), pp 269-278, July 1994

[8] Akkouche S. and Galin E.,
*Adaptive Implicit Surface Polygonization Using Marching Triangles*,
Computer Graphics Forum, Vol 20 no 2, pp 67-80, June 2001.

[9] Pitman, B.
*Introduction to Numerical Analysis*
(http://www.math.buffalo.edu/~pitman/courses/mth437/na/na.html)

[10] Velho, L.
*Simple and Efficient Polygonization of Implicit Surfaces*,
Journal of Graphics Tools, Vol 1 no 2, pp 5-24, February 1996

[11] Livnat, Y. et al,
*A Near Optimal Isosurface Extraction Algorithm Using The Span Space*,
IEEE Transactions on Visualization and Computer Graphics, 1996

[12] Kobbelt, L. et al,
*Feature Sensitive Surface Extraction from Volume Data*,
SIGGRAPH 2001, Computer Graphics Proceedings, Annual Conference Series, pp 57-66, 2001

[14] Ju, T. et al,
*Dual Contouring of Hermite Data*,
ACM Transactions on Computer Graphics (SIGGRAPH 2002 Proceedings), pp 339-346, 2002

[15] Schaefer S. and Warren J.
*Dual Contouring: The Secret Sauce*