Tracing metaballs in two dimensions
This article describes an algorithm for efficiently tracing the outline of two-dimensional metaballs. A metaball (not to be confused with meatballs, which, while delicious, are not the topic in which we are presently interested) is a computer graphics trick developed in the 1980s that creates globule-like shapes not unlike the blobs in a lava lamp. Despite mimicking the appearance of a fluid, as the #metaball hashtag on Twitter succintly demonstrates, the actual math behind metaballs is surprisingly simple.
In two dimensions, a metaball is a type of isoline: a line that connects points of equal value on a plane, where the value of each point is determined by some function. Chances are you've seen isolines before. The elevation lines on a topographic map, for example, are isolines that connect points of equal altitude. Isobars on a weather map are isolines that connect points of equal atmospheric pressure. The same concept that creates a two-dimensional metaball can be (and often is) extended into three dimensions to create an isosurface, which is, as you guessed, a three-dimensional surface connecting points of equal value in space.
Should you need to actually trace and then fill the outline of a two-dimensional metaball, however, neither of the two aforelinked techniques quite fit the bill. The former method leaves us with a jumble of tiny line segments and no way to stitch them together into subpaths, rather like a game of virtual pick-up sticks. The latter method, despite its spectacular appearance, leaves us with no lines at all, and thus no way to (say) outline the shape with one color and fill the shape with another. To get where we want to go, we need something slightly different, so we'll borrow the sampling technique described in Jamie Wong's article and then strike out in search of a way to derive the actual metaball path. Our destination is an algorithm fast enough to draw metaballs composed of several tens of circles in real time on a modern PC.
Metaballs begin with circles, so let's review the math behind circles. If you and circles are best buddies, feel free to skip this section. If you need a refresher, fear not: we'll use just a smidge of algebra and a dash of geometry ahead.
We know that when we talk about a circle, we're really talking about a set of points: pairs of x and y values that identify locations in two-dimensional space. More specifically, we're talking about all points that are all the same distance from the center point. We also know we can't possibly list out every one of these points, what with there being an infinite number of points on any circle. Fortunately, we also know that with a single equation, we can describe all the points that lie precisely on a circle with radius r:
x2 + y2 = r2
But wait; how do we actually know that this equation is talking about a circle? Being the clever reader, you no doubt recognize this equation as none other than that old classic, the Pythagorean theorem. In fact, if we imagine a circle centered on (0, 0) and we draw any path from the center to the edge, then straight down (or up) to the x-axis, and finally back to the origin, we've drawn ourselves a nifty right triangle where the hypoteneuse is length r, the vertical side is length y, and the horizontal side is length x. So for any such triangle we could possibly draw (there being an infinite number of those, too), we know there is some x and y that satisfies the above equation.
But for the sake of illustration, let's suppose we didn't know about the Pythagorean theorem because we live in an alternate reality where Pythagoras spent most of his time playing Greek checkers and the Babylonians, who may very well have been the first to happen upon the equation above, never existed. If you and I were toiling in the mud (as we would undoubtedly be in a world lacking such essential mathematical knowledge) and we unearthed (unmudded?) an obelisk bearing the inscription x2 + y2 = r2, how could we know what it means?
One thing we could try is to sample the equation for various values of x, y, and r. We can't try every value (infinity and all that), but we can try a large enough number to get a rough idea of what the equation is trying to say.
Before we start, we'll make one small adjustment to the formula to get everything on one side.
r2 ÷ (x2 + y2) =1
Arranging the equation this way tells us that the number 1.0 is probably important, perhaps as a natural dividing line between one set of values from the left side of the equation and another. Let's represent this visually by plugging in a range of values for all three variables. Each pair of x and y will represent a square in an imaginary grid, and we'll color in that square only if the resulting value is greater than or equal to 1.0.
As long as we take enough samples, we can experimentally verify that the equation does seem to describe a circle. And if we take the same samples for different values of r, we see that the point where the samples cross the x- or y-axis is roughly a distance of r from the origin. Naturally, given our current hypothetical employment as mud-toilers, we would likely not achieve this level of insight and instead see the mysterious obelisk as a fine platform for forming mud-cakes and thereby becoming the envy of our fellow mud-ensconced workers.
If, however, we were not quite so intellectually impoverished and were to think more abstractly about the concept of sampling, we would realize that all sampling involves a tradeoff between quality (how accurately the samples represent the underlying subject) and efficiency (how much time or space we need to compute or store all the samples). The quality depends directly on the resolution—the number of samples we take in a given area. If the resolution is too low, the image devolves into a collection of unrecognizable blocks. As we increase the resolution (and thus the number of samples), each sample represents a smaller and smaller portion of the subject, and fine details previously smudged out by low-resolution samples appear. A digital camera operates on the same principle: picture quality improves as the sensor resolution increases and the amount of light represented by each “sample” decreases.
If we were mathematicians enamored with the concept of infinity, we would jump to the conclusion that if our resolution were infinitely high and our samples infinitely small, we would be able to represent the subject with infinite precision. True though this may be, we are programmers, and as programmers we tend to think much about the other side of that tradeoff: efficiency. If we want to sample each pixel of a 1920- by 1080-pixel image, that's more than two million calculations, not to mention the memory required to store those samples. The resulting image would be quite faithful to the original, but depending on the complexity of each calculation, the time we can afford to spend on sampling, and how much processing power we have at our disposal, the task may be simply too big.
Our job, then, is to find a way to compromise between quality and efficiency that gets us to the point of “good enough to fool the user”.
From one to many
With a firm grasp on the concept of sampling, we can now experiment with multiple circles. The usual circle equation assumes the circle is centered on the origin. We need only make a minor change to allow for circles centered on any point (x0, y0).
r2 ÷ ((x − x0)2 + (y − y0)2) =1
(Subtracting x0 and y0 has the effect of temporarily shifting (or translating) the circle back to the origin so that the original equation will work as we expect.)
Let's supply ourselves with an assortment of circles of various sizes:
With more than one circle in play, what does this mean for our sampling? Previously we had only one equation to evaluate, so we had only one value for each sample. With ten circles, each sample could yield a different value depending on which circle we look at.
What we need here is some way to aggregate (or reduce) the values from each circle into a single sample value. Let's see what happens if we take the largest value at each sample point.
Taking the largest value reveals the locations and bounds of the original circles, which is not all that exciting. What if we instead add up the values from each circle and use the sum as the sample value?
Now we're on to something. When we sum the values, the boundaries of nearby circles “reach out” toward nearby circles, creating that liquid-like congealing effect that is the hallmark of metaballs.
To understand why the “congealing” effect is limited to only a circle's nearest neighbors, we can make a minor change to the way we draw the sample squares. This time, rather than turning each sample square “on” or “off”, we'll set the opacity of each square to the sample value, clamped to the range [0, 1]. This will make larger sample values appear more opaque and smaller values appear more transparent.
With this visualization, we see that the locations with the greatest chance for the sample values to add up to at least 1.0 are either inside a circle (naturally) or in between two nearby circles. Beyond a circle's boundary, the sample values quickly fall toward zero, so each circle “contributes” little to the total sample values in distant locations.
This makes sense if we think about the equation. For a given circle, the numerator (r2) remains constant, but the denominator (x2 + y2) grows rapidly the further we go from the circle's boundary. A rapidly increasing denominator makes for a rapidly decreasing quotient.
A walk through a field
Unlocking the math behind metaballs gets us one step closer to our final goal. All that remains is to turn these piles of chunky pixels into an actual path.
Recall that a metaball is an isoline: a path connecting points of equal value according to some function. We have our function (the sum of the left side of all circle equations) and a grid of samples of that function's outputs. Somewhere in that scramble of numbers is a path that neatly connects the points where the sample value is exactly 1.0. We want to reproduce that path as accurately as possible given the sample resolution.
We can start with the important observation that the path must go between (or precisely through, but more likely between) pairs of sample values. Say we have two neighboring sample squares. The value of one is 0.88, while the value of the other is 1.42. We know that the former lies outside the path because it is less than the threshold (1.0), but the latter lies inside the path. Ergo, the path must travel between these two samples.
If we start from the top-left corner and work left-to-right, top-to-bottom, we'll eventually reach a sample value greater than or equal to 1.0, and this sample will have a neighbor to the left with a value less than 1.0. From that point we should be able to trace out the rest of the path. Let's test this with a Z-tetromino shape.
The first square we encounter with a sample value greater than 1.0 is at (2, 0). The square to the left has a sample value less than 1.0, so we add the point to the left of (2, 0) to the path. We then make one quarter-turn clockwise to examine the sample above (2, 0).
The coordinate (2, -1) happens to lie beyond the edge of the sample field, so we don't have a sample value. For the sake of convenience, we'll consider all out-of-bounds coordinates to have a sample value of zero. This means that the path travels between (2, 0) and (2, -1), so we add another point to the path and turn clockwise.
(3, 0) also lies outside the path, so after adding the point to the right of (2, 0) and turning clockwise we are now looking at (2, 1). (2, 1) is inside the path. We handle this case by moving the “cursor” (or “turtle” if you come from a Logo background) into that square and making one quarter-turn anticlockwise.
This leaves us facing (3, 1), which lies outside the path and thus leaves us with another point and another quarter-turn to the right.
Repeating these steps (adding a point and turning clockwise if the adjacent square is outside the path, advancing and turning anticlockwise if the adjacent square is inside the path) eventually lands us back where we started. At this point, we know the path is complete and we can examine the remaining samples.
As we encounter other samples inside the path (such as at (1, 1) and (1, 2)), we check to see whether we have already placed a point to the left. If we have, it means we've already discovered this path and can skip the process above.
Applying this algorithm to our friendly set of circles yields the following metaball-like blob.
Polish and shine
Cranking up the resolution with the above algorithm helps, but something's still not quite right. Because we always take the midpoint between two samples, every metaball ends up being a jagged path made of exclusively vertical, horizontal, or diagonal line segments. Can we do better?
Let's put on our thinking caps and ponder the case where we've determined that the path passes between two sample points, one on the left and one on the right. Since the sample points are horizontal neighbors and therefore share a y-value, the point we're looking for must have the same y-value. What we don't know is the x-value. Ideally we want to choose a value for x that's as close as possible to the point where the sample value would equal 1.0.
One strategy would be to interpolate between the two sample values. “Interpolation” essentially means we make guesses about the data we don't have based on the data we do have. There are a variety of fancy ways we can guess about missing data, but the simplest and least costly method is linear interpolation, which assumes that the transition between two values follows a straight line.
As an example, suppose we're tracking the population of a city. On January 1st of one year, we take a census and discover that city's population is 96,400, plus or minus a few dogs and cats. On January 1st of the next year, we take the same census and find that the population has increased to 101,200. When did the population reach 100,000?
We can't answer this accurately with just two points, but we can assume that the rate at which the population increased remained steady over the course of the year. Based on that assumption, we would expect to find that the city's population crossed the 100,000-mark sometime in the second half of the year, since 101,200 is closer to 100,000 than is 96,400. Let's do the math for a 365-day year:
(100,000 − 96,400) ÷ (101,200 − 96,400) ×365=273.5
So based on the best data we have, we guess that September 30th, being the 273rd day of a 365-day year, was the day the population reached 100,000. This may be inaccurate for a variety of reasons, of course. Perhaps the population actually decreased during the first half of the year before booming in the second half. Having higher-resolution data (such as census results from every month) would improve the accuracy of our interpolation. But with only two data points at our disposal, this is the best we can do.
Since our path-finding algorithm always compares two sample values at a time, it aligns perfectly with the interpolation example above: given two sample values on opposite sides of 1.0, where did the metaball function cross 1.0? In reality, we know the metaball function doesn't produce a straight line of sample values. But we anticipate that if our resolution is sufficiently high, the distance between samples is small enough that we can treat the change between neighboring samples as linear. Let's see if this holds true.
Sure enough, interpolating between sample values improves the appearance greatly without changing the sample resolution. But linear interpolation tends to produce odd points and protrusions where two metaballs come close to each other without quite touching. These points are where our assumption about the linearity of the sample function fits the underlying data most poorly. Increasing the sample resolution helps, but every doubling of the resolution quadruples the number of samples we have to take and costs us more time and memory.
Rather than increasing the resolution globally, what if we increase it locally? Given a sample value on the left of 0.9 and a sample value on the right of 1.1, we can calculate the sample value for the point exactly halfway between the left and right samples. Suppose this middle sample value is 1.05. Now we know that 1.0 lies somewhere between the left point and the midpoint. We can repeat this process (calculating the sample value at the midpoint and narrowing down the range) either a fixed number of times or until we get within some small distance of 1.0. Let's call this process “midpoint with iterative refinement” and see if it helps any.
Indeed it does! We were able to eliminate many of the unsightly jagged edges without paying for increased resolution.
There are both obvious and non-obvious real-world applications for the path-tracing algorithm described above.
Perhaps the most obvious is the rendering of elevation lines in maps as mentioned in the introduction. Given a grid of elevation samples, repeating the metaball algorithm for various thresholds (such as every 100-meter increment) and overlaying the outlines gives you a topographic map.
Metaballs can also be applied to scatter plots similar to those produced by the Gapminder tool popularized by the late Hans Rosling. Rather than plotting each point individually, a metaball can better show the overall range of a group of points and how that group compares to other groups.
A less obvious application is in identicons, particularly the kind used by Move by Numbers, a multiplayer chess game developed by Eight Squared Software.
At first blush these identicons appear to be simply made of light and dark squares arranged on a 64-square grid. A closer inspection, however, reveals that the light foreground area is actually a single path drawn above a background rectangle.
As it so happens, the same path-tracing algorithm that works for metaballs also works for tracing the above identicons' foreground shapes. But why go to such lengths? Surely a collection of light and dark squares would produce the same result, would it not?
The answer to this question is that resizing that collection of squares may produce artifacts—graphical undesirables—depending on the ratio of the original image to the resized image. When drawn as individual squares, despite sizing each square so that it extends completely to the edges of neighboring squares, resizing the image sometimes reveals hairline gaps between squares and allows the background color to show through. When the contiguous foreground regions are drawn as a complete path, however, these gaps disappear completely at any image size.
Depending on the slice of the programming universe you find yourself in, you may never need to actually code up the algorithm above. But programming is often just as much about having a well-stocked toolbox of ideas and principles as it is about raw code-slinging prowess. Sampling and interpolation are useful techniques that have far-reaching applications beyond metaballs and beyond even the realm of computer graphics, and will undoubtedly serve you well.