Tag Archives: isosurfaces

Distance Fields for Smoothed, Convex Polyhedra

Today was a very fun math day, and I want to share a bit of my experimentation.  Of all the things I did today, one of the most useful was to come up with a formula for the distance field of a convex polyhedron, defined by some set of planes.  I also implemented a smoothed version, which is the one I'll post, because it's quite interesting.

Disclaimer: I'm calling it a distance field because it looks like one and because I don't see many people properly differentiate (teehee?) between distance fields and things that look like distance fields but aren't actually.  So I'll join that club.  The 0-isosurface is the convex polyhedron (or smoothed version thereof), and it does exhibit a nice differentiable behavior, but I don't think that the scalar value actually represents a distance.  Certainly not a Euclidean one if it does.

Random, Smooth Convex Polyhedron


Suppose we have a set of k planes defined by normals {n_1, ..., n_k} and d-values (the dot product of some point on the plane with the plane normal) {d_1, ..., d_k}.  Define the shape factor, p, to represent how sharp the edges of the polyhedron are, where p = 2 corresponds to sphere-like edges, and p = \infty corresponds to hard edges.  This definition is a reference to the p-norm, to which the following formula is closely related.  The field value at a point x is given by:

f(x) = \left(\sum_i^k{\left({\frac{x \cdot n_i}{d_i}}\right)}^p \right)^\frac{1}{p} - 1

Just so we're clear, that inner bit is the dot product of the position in the field x and ith plane's normal, n_i. Intuitively, you can see that, with the \infty-norm, this gives the maximal (scaled) distance to any plane - 1, which will be 0 if we are on the polyhedron, and will fall off with the closest plane's distance.  This corresponds to our idea of what should happen for a hard-edged polyhedron.  For other norms, you can verify that it gives interesting, smooth results.

Now comes the real work.  This is great, but we also want analytic gradients for our fields, so that we can compute perfect normals for the surfaces rather than using ugly central differences.  Perfect normals go a long way towards making the polygonized result look good.  It's a bit painful, but if you just remember your multivariable calculus (I'll admit I struggled with trying to remember it for quite a while), it's not too terrible.  It's just application of a bunch of  rules.  I'm going to spare the derivation and give the result:

\nabla f(x) = \left(\sum_i^k{\left({\frac{x \cdot n_i}{d_i}}\right)}^p \right)^{\frac{1}{p}-1} \left(\sum_i^k n_i {\left({\frac{x \cdot n_i}{d_i}}\right)}^{p-1} \right)

The left side is a scalar factor, closely related to the field value.  The right side is the actual vector component, which is just a sum of scaled plane normals.  All of it can be easily implemented in a single loop over the planes!

With these formulas in hand, you can create lots of cool, smooth shapes. Boxes, discrete cylinders, prisms, etc...or just random convex shapes!

Smooth Convex Cylinder

Smoother Convex Cylinder


If you're wondering how I came up with these formulas, it was by studying and generalizing the field formula for a smooth box, which is pretty obvious. This polyhedron formula is less obvious, in my opinion, but is still basically just a straightforward generalization of the box formula.

PS ~ As usual, I invite any lurking mathies to speak up if I've made a blunder :)

PPS ~ I hope you like the Latex plugin that I installed, solely for the purpose of this post..


99% of the time, algorithms end up being way more frustrating and way less effective in code than they were in your head. Then, there's that 1% of the time wherein divine powder falls from the sky, striking your code at a secret, magical angle that induces ridiculous levels of wow-this-just-works-ness. And then you ask "where's the catch?" and hope that there isn't one.

Such has been the story of SurfaceWrap, my work-in-progress surface-based scalar field polygonization algorithm. It's remarkable how well it's working so far. Granted, it's still not adaptive nor self-stitching, so I'm hesitant to terminate my relationship with marching cubes just yet...but things are looking really, really good. I've run into a surprisingly small number of headaches, especially considering how hard I usually fail at writing topological algorithms.

Here's the current state of the algorithm:

A few things to notice:

  • The triangles cover the surface very smoothly and gradually without slivers or jaggies
  • The consistency of the triangles is really nice! Looks like this algorithm is going to dominate with symmetric surfaces
  • Approximate normals are calculated for free (no gradient evaluations!) during the process! BIG WIN!
  • Although stitching isn't yet implemented, no self-intersection has been observed (so long as the boundary isn't allowed to wrap back around on itself). Looks like the small-angle-elimination feature is keeping everything in check

Next step: boundary merging/splitting. Pray for me.

Surface-Based Polygonization

I've set aside the physics engine for a while to undertake another huge project - one that lies at the core of everything I do. My new goal: eliminate marching cubes. It's slow, it's memory intensive, it has no respect for surface topology, and it's not easy to make adaptive.

Instead, I'm building an algorithm that will polygonize a surface using....the surface! Just like my adapted version of marching cubes ('crawling cubes') crawls along the surface, so will my new algorithm. But, unlike MC, my algorithm will orient triangles such that they actually approximate the surface.

This is definitely going to be the most challenging algorithm I've ever designed, and I'm sure edge cases and slight things that I haven't considered yet are going to turn this into an absolute nightmare...but the reward for success is huge: polygonization quality and speed are absolutely essential to my entire engine, since all meshes will be converted at some point. MC is currently responsible for many of the annoyances that I see on a daily basis: poorly-calculated normals (due to sliver triangles), high poly count (due to the lack of adaptiveness), meshes with holes (due to lack of a bulletproof method for determining the isosurface bounds). A good surface-based algorithm will fix all of the above.

At any rate, I'm excited to have some initial results to share! It's not much, but you can see the formation of something that loosely resembles a surface! There's no adaptability yet, and the surface doesn't handle self-intersection (which I anticipate will be the most difficult part of the algorithm). Still, it's an encouraging start, and the surface quality clearly beats out anything MC could produce.

Procedural Trees as Isosurfaces

To demonstrate the capabilities of the advanced metaballs described in the last post, it seems fitting that I should try to model something not easily modeled with standard metaballs. Trees certainly come to mind. The asymmetric and non-axis-aligned qualities of tree models make them difficult to recreate with isosurfaces.

With the power of advanced metaballs and a few hours of work, I was able to come up with some shapes that resemble alien-like trees. No, they're not perfect evergreens, but they're a start. They're a heck of a lot more interesting than the cylinders with billboarded leaves that you see in many low-budget games. Most importantly, however, they show that complex mesh are possible with metaballs.

Advanced Metaballs

Ever since the successful implementation of Marching Cubes into my engine, I've been experimenting with ways to model things procedurally. Naturally, metaball modeling is the first thing that comes to mind. Unfortunately, metaballs, as described in literature, are somewhat restricted. They usually contain four pieces of information: x, y, z, and radius. The first three control the origin of the metaball's field of influence, and the radius describes how far that field extends.

In an attempt to allow more flexibility, I introduced several new parameters: scale and power. This adds an addition six pieces of information to the metaball (scale and power are each three-dimensional vectors). The scale controls how far an influence field extends in each dimension. The power then controls how quickly the field falls off in each direction. Interestingly, one can achieve completely different primitives just by changing the power: spheres, cylinders, capsules, and even cubes.

Even with these additions, metaballs did not have the power to recreate more complex topologies. For this reason, I began to explore rotation as another parameter of the metaball. After much trial-and-error, brushing up on matrix math, and poking around forums, I was able to get rotation right.

The advanced metaball, as I casually call it, has still four more pieces of information on top of the position, radius, scale, and power information contained in the original. It contains a 3-component vector that specifies a rotation axis, and a scalar rotation angle. Together, these components tell the metaball how to orient itself about an arbitrary axis of rotation.

With the power of rotation now at hand, it remains to be seen just how flexible advanced metaballs will prove to be as tools in implicit surface modeling.

Isosurface Extraction Optimization

I've spent an insane amount of time over the past few days trying to optimize the Marching Cubes algorithm, since it's going to be such a pivotal part of my project. Ideally, I'd like all meshes to be parametrized in terms of isosurfaces for memory efficiency as well as for the ability to quickly prototype models.

I experimented with both octree culling as well as optimizations to the original, brute-force algorithm. Surprisingly, the brute-force algorithm ended up outperforming the octree in all cases. Furthermore, I failed to find a reliable way of predicting whether or not an octree node would intersect the isosurface given just the density values at the corners, so reconstruction often produced meshes with holes in them.

Ultimately, I ended up cutting polygonization time down by a factor of four. I'm still using a brute-force technique. However, using some rather hacky methods to perform linear interpolation on selective parts of the density field, I've managed to avoid having to directly compute the entire field. This method basically involves computing a derivative of the field and then using a very reserved linear approximation of the density field to skip a few computations if the values are determined to be at a "safe" distance from the isolevel.

With this technique, the engine can build a 60k poly model (high quality) of a density field containing roughly 30 metaballs in less than 2 seconds. That's pretty fast! I'm going to continue to seek optimizations for this algorithm, but for now I'm much more pleased than I used to be with the speed. Isosurface extraction is starting to look very feasible for a game engine!

3D Perlin Noise

With the power of a working isosurface extraction algorithm, I decided, quite naturally, to test the capabilities of marching cubes on the 3-dimensional analogue of my favorite noise function - perlin noise.

Coding a perlin noise function to work in tandem with the marching cubes algorithm was simply a matter of transferring the 3D perlin noise function that I had already written in HLSL over to equivalent CPU instructions.

The results, though quite intriguing to look at, aren't particularly beautiful.

Unfortunately, 3D perlin noise does not seem to be feasible for real-time application. Generating one or two meshes that use 3D perlin noise may be acceptable, but making heavy use of the function simply isn't feasible if loading times are to be kept reasonable. It's slow. Really slow. Perlin noise is already a computationally-expensive function, and the fact that I'm now running a 3D analogue of it (which is, by nature, more expensive than 2D) on the CPU rather than the GPU makes the problem intractable. The CPU simply isn't capable of crunching the math fast enough.

In the long run, I'm going to look into alternatives to 3D perlin noise. The function would be a great asset for a new reality, but it's not practical for real-time application at the moment - at least on the CPU. In the future, I'll try to look into DirectX 10's stream out, as I know that it can be used to perform marching cubes on the GPU (probably in a ludicrously small fraction of the time that it takes the CPU).

Ellipsoid Modeling - Starship

Using an isosurface extraction algorithm, it is possible to model in far more intuitive ways than with extrusions, welds, and the other tedious tasks of typical 3D modeling software. One such method is via the use of ellipsoids that influence the density function from which an isosurface is extracted. Modeling using ellipsoids is fast, intuitive, and requires only a fraction of the space that a regular mesh would take to store.

The following simple starship mesh was created in about ten minutes using only 14 ellipsoids.

Note that ALL of the assets used in this demo, including my core engine, the precompiled drawing shaders, as well as the ship meta-mesh and extraction algorithm, fit into a single 33kb executable (after UPX compression). And no D3DX dependency, either! I'm quite proud of how space-efficient my project is turning out to be.

All hail the marching cubes!

Marching Cubes Algorithm Implemented

After a great deal of problems and long nights of debugging, I finally have the marching cubes polygonization algorithm up and working in my engine! I was hesitant to implement the algorithm in the first place, since it requires a lot of constant global memory (several KB for the lookup tables), and I'm trying to keep my executable size as low as possible. In the end, however, it was worth the few extra kilobytes. The results are incredible! Isosurface extraction may be the answer to all of A New Reality's problems!

Using isosurface extraction, it will be loads easier to create procedural models. My fear, however, is performance. The algorithm is sluggish at best, even on powerful CPUs. No wonder, considering it has to sample tens, possibly hundreds, of thousands of grid points in the density function. For high-quality procedural modeling to be feasible using isosurface extraction, I'm going to have to hope that CPU power increases a lot in the near future. Until then, I will continue attempting to optimize the algorithm as well as explore multithreaded options to prevent lockups during mesh creation.

I hope to post some results of the algorithm soon!