The problem
Say you are using some random implicit function, and imagine we are using this function to draw some shape (defined by it's f = 0 isosurface). Let's take one simple formula as an example:
where, as usual, r = sqrt(x˛+y˛) and a = atan(y,x), meaning f is in the end of the day a function of regular x and y cartesian coordiantes. This formula produces a simple 3-lobe fan shape, and it could easily be one of the many elements making a bigger procedural image/texture or something. Now, if you were to do nothing special but just assign a color to each pixel based on the value of f, in a way similar to the code below does, then you would be producing an image like the one here in the left
float color( in vec2 x )
{
float v = f( x );
return smoothstep( 0.19, 0.20, abs(v) );
}
|
 |
 |
||
| code for color(x) = ramp( f(x) ) | color(x) = ramp( f(x) ) | color(x) = ramp( distanceEstimation( f(x) ) |
Interesting as it is, imagine that it didn't look as you wanted, and that in particular the non constant thicknes of the shape was a problem. Perhaps you would be more interested in a constant thickness outline image, probably based on a more homogeneous scalar field. In fact, if you could compute some sort of distance to the zero-isoline of f, that would be great, wouldn't it? (now you probably see how this is related to doing distance based raymarching and you try to find intersections with the zero-isosurface). Perhaps, you would be interested in being able to compute an image like the one on the right. Let's see how we can achieve this!
The maths
|
Let's do some maths in order to understand the technique. The implementation is going to be terribly simple and you can skip this paragraph and go straigh there,
although I always think it's a good idea to understand what's going on behind the code. So, we start by having a look to the distance field of our example function f. The image to the right of this text shows it in grey scale, on which I have mapped the values of f to a shade of grey such that f = 0 is represented in black and |f| > 0 shows a brigher color, through a pow() function for the purpose of making the f = 0 isoline more pronounced. Of course, in the interior of the shape f < 0 and in the exterior f > 0 Now, we are interested to compute the distance to our f = 0 shape at a given point in space x, which is represented by the yellow dot. The distance to the isoline is the distance to the red point, which is the closest point in the isoline to the yellow point. Let's call that closest red point x+e, such that e is the vector going from the yellow to the red one. What we are after here is |e|, the length of e, which is the distance from the x to the isoline/isosurface. In this setup, since x+e is in the 0-isoline, we clearly have  Let's assume we are pretty close to the shape, meaning, that |e| is small. Then, we can expand f(x+e) in its Taylor's decomposition  where the dot (ˇ) means dot product, as usual (remember that both the gradient and e are vectors here). If we where close enough, indeed, then we could just use this linear approximation of f and say that  meaning that  the first inequality being due to the triangle inequality, and the second one to the basic properties of dot product. From there, we could isolate  which gives us an upper bound to the estimated distance from x to the zero-isosurface of f. Another way to arrive to the same results would be to construct a plane tangent to the surface at x+e, defined by the gradient/normal of the surface, and intersect a line going though x that is parallel to the gradient/normal with the plane. |
 distance field for note the dark spiral that leaves from the 3 spikes of the main shape, that explains their elongation when f was rendered by direct color mapping in the first image in the header of the articles. |
The implementation
Of course, most of the times we won't have an analytical gradient vector for our distance field, so we are forced to do a numerical approximation with the regular central differences method
vec2 grad( in vec2 x )
{
vec2 h = vec2( 0.01, 0.0 );
return vec2( f(p+h.xy) - f(p-h.xy),
f(p+h.yx) - f(p-h.yx) )/(2.0*h.x);
}
and then use the formula above to compute the distance estimation (de in the code) to compute the color as before:
float color( in vec2 x )
{
float v = f( x );
float g = grad( x );
float de = abs(f)/length(g)
return smoothstep( 0.19, 0.20, de );
}
This method evaluates f four extra times, which makes it more expensive. Besides, choosing the right value for h might be difficult sometimes.
But more often than not, it simply works pretty well. When analytical gradients can be computed, it is strongly recommended to use them instead, like when rasterizing
2D fractals or raymarching mandelbulbs or
julia sets.This method shines for computing distances to 1D functions. In fact, the gradient of f is perpendicular to it's tangent/derivative, so it simplifies to (1, f'), making the distance estimation even simpler:
This is particulary nice to graph 1D functions in tools which scan the xy plane and compute a color from x and y, like Shadertoy.
 using |f(x)-y|. note how some parts of the graph go so thin that become hardly visible |
 using a scaled |f(x)-y|. some parts are still too thin, and others got too thick |
 using |f(x)-y| / sqrt(1+f'(x)˛) solves the problem and produces a much more uniform curve |