unity - Top-Down Octree Generation of Procedural Terrain

I'm trying to implement a voxel-based terrain generation system in Unity3d (C#). I have successfully implemented a uniform 3d grid system, and have extracted the isosurface out using Marching Cubes and Surface Nets algorithms.

I quickly bumped up against the inherent issues in representing all space with a 3d grid. Much of this space was either above or below the surface, and I turned to space partitioning with octrees, as it didn't seem too difficult. One of the great things about octrees is that octree nodes that that don't have edges intersected by the surface don't need to be split again.

I researched, and found a couple resources to construct my code. One is Volume GFX, and another is the original Dual Contouring paper by Ju et al. My edge-checking is done by the Marching Cube's checking from Paul Bourke's code. If the "cubeindex" is either 0 or 255, then no edge is intersected, and the octree node does not need to be split.

My Octree generation code works for something like a quarter-sphere, where the surface features are extremely normal:

As seen in the picture, only octree nodes containing the surface are subdivided. Working great.

Now, however, when we move to something more complex, like Unity3d's PerlinNoise function: GASP! What's that hole doing in the mesh at the lower-right corner? Upon closer inspection, we see the octree didn't subdivide properly (red lines highlighting the offending node's dimensions): This turns out to be the same issue Jules Bloomenthal highlights in her paper on polygonization, page 10, Figure 10. Traditional methods at generating top-down octrees ("Adaptive Subdivision"), are sensitive to fine surface features relative to the size of the octree node.

The gist:

X-------X
|       |
|       |  <- Node
X--/\---X     X's - tested values, all outside of the surface!
  /  \  <- surface

The surface breaks the surface, but comes back down before it goes over a vertex. Because we calculate edge intersections by looking at the signs at the two vertices, this does not flag the edge as crossed.

Is there any method to determine if these anomalies exist? Preferably, the solution would work for more than just the Unity3d noise function, but for any 3d noise function (for cliffs/overhangs/floating islands, etc).

UPDATE

Thanks to Jason's awesome answer, I was able to answer my own questions I asked (in the comments under his answer). The issue I had was that I didn't understand how I could create a bounding function for trigonometric functions (sine,cosine,etc) due to their periodic nature.

For these guys, the periodicity is the key to bounds functions. We know that sin/cos reach their extrema values at a set interval, specifically every π/2. So, if the interval we are checking contains any multiple (cos) or half-multiple (sin) than 1/-1 is the certain extrema. If it doesnt contain a multiple (i.e. the interval [0.1,0.2]), then the range is simply the values of the function evaluated at its endpoints.

UPDATE 2:

If the above doesn't make sense, check Jason's answer to my comments.

Answer

I will assume that you have a function f such that the surface is defined by f(x,y,z)==0. Right now, all you have is the function f itself, and so with no further information, it is impossible to know when to subdivide. It is always possible for the surface to do what you describe, i.e. the function f can have arbitrarily thin "fins" which can require you do go arbitrarily deep in the octree.

The only solution that allows top down generation is to have additional information on f. Right now you have a point query: you can determine what the value of f is at any particular point, but you have no guarantees about what the value will be at any other point, even ones right next to it. Instead, you need a bounds query: given a AABB, it determines what the maximum and minimum value of f is over that range.

Now it is simple: for each octree node, perform a bounds query. If the minimum is greater than zero or the maximum is less than zero, do not subdivide. Otherwise, subdivide unless you are at the smallest size, in which case you can use the point query on the corners as usual to construct the surface.

It is easy to write the point query for interesting functions, but it is not clear how to write the corresponding bounds function. However, our job is made easier by the fact that we don't need to be precise: calculating the exact min/max is not going to be particularly useful anyway, because we are just going to compare them to zero. So really we just need to calculate and interval [a,b] which completely contains [true_min, true_max] over the given AABB bounds. If we are close to the true range, we will subdivide mostly things which actually contain the surface, and if we are loose and give much bigger ranges we will end up dividing more nodes than is necessary. However, as long as we make sure to contain the true range we will never miss part of the surface.

As an example, we can construct the bounds function for the sphere function. (To see why your technique doesn't work even in this case, start with a sphere which is entirely contained in the starting cube. The algorithm will terminate immediately!). The point function is f(x,y,z) = dist(, center) - radius. We see that this function depends only on the distance from the point center, so is really at its core a 1D function. First, we determine which range of radiuses our AABB covers. A simple but slow way to do this is to determine the distances for the 8 corners and take the minimum and maximum distances. The only special case is that if the center is contained in the AABB, then the minimum is zero. Now we have a range of radiuses, and by plugging in the function f(r) = r - radius will give the desired range for values of f.

Of course, that was a simple example; you probably want something more complicated than spheres. The nice thing about this approach is that you can combine such bounds functions. For example, if you add two functions together, say f + g, the bounds function for the sum can be easily constructed in a conservative manner as: the min and max are just the sum of the mins and maxes respectively of the two functions. To see this is conservative, consider the functions sin(x) and -sin(x) in 1D. Clearly over any range greater than 2pi these each have a range [-1,1], and so the calculated range is [-2,2], but the real range is the interval [0,0] containing just 0! All hope is not lost, however: if we zoom in on a small range for x, then the interval gets smaller (closer to [0,0]) as well.

This means if you can bound one octave of Perlin noise, then you can bound a full fractal noise, which is just a sum of many individual octaves. Actually, if we can bound individual functions, we can use the rules of interval arithmetic to construct more interesting functions. In particular for procedural synthesis, summation, minimum/maximum, domain distortions, and scaling by constants can all be handled by interval arithmetic, which let us add surface texture, place objects, and distort & scale without needing to write complicated bounds functions. So we can write f(x) + min(g(x), h(x)) * 4 + h(2*x) for the point query, and automatically (but imprecisely) construct the corresponding bounds function, assuming we have the bounds functions for f, g, and h.

Finally, some functions are not easily expressed as a simple combination of others. In that case, you can often bound the derivative and use that to bound the value of the function. This is somewhat more math intensive though, so I won't cover it here.

Implementation note

It is probably easiest to represent a function f by a pair of functions:

1. a point query f_point which is a function taking a point to a float.


2. a bounds query f_bounds which takes a AABB to an interval [min,max].

Then you can write combinators like Sum, Scale, etc. which take these pairs and return new ones which are the combination given. You can represent these as tuples, objects, etc. depending on the language.