This notebook walks through a strategy for adaptive contouring in a fragment shader, with particular application to domain coloring. Domain coloring is a simple and well-explored technique for visualizing functions in the complex plane using contours for magnitude and color for phase
This notebook walks through a strategy for adaptive contouring in a fragment shader, with particular application to domain coloring. Domain coloring is a simple and well-explored technique for visualizing functions in the complex plane using contours for magnitude and color for phase.
Over the past few years and with particular help and inspiration from Jacob Rus, I’ve slowly iterated on my approach to domain coloring. A couple of the core ideas have proved so useful and generally applicable that I’ve wanted to write them up carefully and clearly for a long time. In this notebook, rather than focusing on the math and interpretation of domain coloring, I’m going to focus on some of the core tricks I use for everything from domain coloring to contour plots, wireframe meshes, 3D surface visualization, gridlines, cartoon edges, and trippy animated gifs.
Even without working through the mathematical interpretation, the plot on the lefthand side above has two particularly obnoxious features (apart from the garrish colorscale). The first is that the lines have non-uniform width, and the second is that the contour spacing is also highly non-uniform, ranging from half the plot size in slowly varying regions to infinitely dense near zeros and poles. Most regions, if zoomed in, would appear as nothing but a solid color. This technique from Inigo Quilez is a wonderful way to resolve the aliasing problems which result from dense spacing, but ideally we’d like to control the uniformity of information a bit more thoughtfully.
Remedying these problems is subject to one primary constraint: computing plots in a GPU shader is highly performant but means that—short of more advanced parallelization tricks—each pixel needs to be evaluated independently from all others. As a result, our ideal approach uses strictly local information to display lines and select contour spacing.
In this notebook, the first in a couple-part series, we’ll walk through the core techniques by which we can resolve the above challenges. We’ll use screen-space derivatives to achieve uniform-width features and the graphical equivalent of a Shepard tone to achieve ideally-spaced contours.
This notebook is presented in the following sections:
- Screen-space-scaled raster contours: Using screen-space derivatives to display uniform-width lines
- Local contour spacing: Using the local gradient to select nice contour spacing
- Blending over octaves: Smoothly blending between locally scaled contouring
- Shaded contouring: Repurposing line contour techniques to achieve nicely antialiased shading
The contour plots in this notebook are computed and displayed in WebGL via the provided GLSL shaders. And although I’ve used the regl library to create them, the techniques should be applicable to any sort of language, whether on the GPU or otherwise. To make understanding easier, I’ve tried to keep the plotting very simple at the cost of luxuries like interactions and axes. See the Definitions section for implementation.
Vector graphics have many desirable properties, but in this notebook we explore lines as a raster image effect. Particular benefits are that the workload is independent of the complexity of the contours, and that we may avoid computing, storing, updating, and rendering explicit contour line geometry, commonly computed via marching squares.
We start by plotting the log of the magnitude of a complex function using a repeating triangle-wave colormap.
Adjust the threshold below and observe that although we can draw bands, we can’t yet draw uniform-width lines.
Applying a threshold to the plot above produces bands, but of course the width is highly non-uniform. To achieve uniform pixel-width lines, we must select just the right threshold at every point.
To select this threshold, we use screen-space derivatives. As the name suggests, screen-space derivatives represent derivatives with respect to pixels on the screen. Computed as the difference between adjacent pixel values—in the context of GPUs called fragments—they’re built into the GPU and are commonly used in mipmapping to avoid aliasing by gauging how fast a texture changes from fragment to fragment.
GPUs evaluate fragments in blocks and offer horizontal and vertical differences between fragments as built-in functions. These differences are equivalent to first order finite difference approximations of the derivative with units of units per pixel. Horizontal derivatives within a block are defined as
where and are pixel coordinates. Vertical derivatives are defined similarly. The screen-space gradient magnitude is then the sum of the squares of the and derivatives.
The key insight is dimensional. Division of a function (in units) by its screen-space gradient (in units per pixel) yields a result in pixels.
Another way to think about it: at pixel scale, the function is approximately linear. Near each contour line (where the triangle wave crosses zero), dividing by the local slope normalizes the function so that it increases by exactly one unit per pixel. The result is the signed distance, in pixels, to the nearest contour.
So if we divide the triangle wave by the screen-space gradient magnitude, we get a value that represents, in pixels, how far each point is from the nearest contour. Thresholding this at, say, 1.5 pixels gives us uniform-width antialiased contour lines everywhere on screen.
Numerical derivatives, even though only accurate to the first order, will be adequate for us. We could of course compute derivatives by analytical methods like automatic differentiation, but in practice, for the accuracy required to achieve the effects in this notebook, that proves to be massive overkill.
Note finally that this doesn’t have to be a two-dimensional function. Since all other transforms, whether 2D, 3D, or otherwise happen before screen space differentiation, all of the techniques in this notebook can be used without modification, for example, to draw wireframes or other function contours on 3D models.
As for implementation, WebGL offers the OES_standard_derivatives extension with the following functions:
dFdx(f): differences between horizontally adjacent values offdFdy(f): differences between vertically adjacent values offfwidth(f): sum of absolute values,abs(dFdx(f)) + abs(dFdy(f))
The function fwidth is often used as a square-root-free proxy for the gradient magnitude but results in about a 40% anisotropy in diagonal directions (i.e. ), so that length(vec2(dFdx(f), dFdy(f))) is preferable for the best quality.
The GLSL function below implements a contouring function that receives input f and returns 0.0 in the background and 1.0 on a contour, with a smooth fade between for antialiasing.
The figures below illustrate usage of this function, first for an increasing function of one variable and then for a two-dimensional function. Modify the slider values and observe how the function is transformed into equally-spaced, antialiased contours.
We now turn our attention to contour spacing. As successful as the above plot is, the spacing for this very tame function still varies from less than a pixel to more than a hundred pixels.
We’ve already established that dividing by the screen-space gradient tells us how many pixels away we are from a contour line. That same gradient also tells us how densely contours are packed at any given point. If the gradient is large, contours are dense; if small, contours are sparse. So we can use this information to choose which contours to draw in any given region.
The approach is to partition the plane into regions we’ll call octaves, where each octave draws contours at a different spacing. Within each octave, we draw evenly spaced contours. When the contours in one octave get too close together, we switch to the next octave and draw times fewer contours.
The key trick is to use a log-ceil-exp transform. Taking the logarithm of the gradient, rounding up to an integer, and exponentiating back produces discrete steps. Regions sharing the same integer octave value will draw contours at the same spacing. And because logarithms convert multiplication to addition, zooming in (which scales the gradient) simply shifts us between octaves rather than continuously changing the contour density.
To select the locally appropriate octave, we use the equation
where is an integer indicating which octave a given region belongs to, is the number of divisions per octave, and is the minimum spacing, in pixels, between divisions.
Why and not just ? The answer is that we want to place contours at evenly spaced values of , not itself. Logarithmic spacing places more contours near zeros and poles, which helps reveal the structure of the function. Since we’re effectively plotting , we need the gradient of By the logarithmic derivative,
This is convenient because it means we don’t have to compute the logarithm to get its gradient. We just divide the gradient of by .
Once we’ve determined the octave, we compute the contour spacing as
The toggle below lets you switch between logarithmic and linear contour placement. For linear spacing, we remove the division by abs(f) and plot abs(f) directly instead of log2(abs(f)).
Whether to place contours linearly or logarithmically is a bit of a subtle decision. I think linearly spaced contours often make more sense, but when shading is applied, it’s easier to tell the difference between roots and poles when contours are placed logarithmically. Fortunately both are valid and we don’t have to commit to just one.
The shader below implements locally scaled contouring. Try enabling the octave visualization to see how the plane is partitioned into regions of equal contour spacing.
Change the parameters below and observe the effect on the contours. Note in particular that changing the minimum contour spacing does not affect the position of the contours. Instead, it only affects the location of the thresholds between the octaves.
At this point, what we’ve created is sort of like a fractal pattern, carefully band-limited so that the smallest-scale pattern is pinned to the scale of just a few pixels. The number of octaves may be chosen by the user, but in fact there is a fairly concrete answer to the question of how many octaves should be used. The number of octaves should be chosen so that the largest scale pattern precisely covers the entire image.
This is where we’ll leave it for now! These are the core components of my currently-favorite approach to domain coloring, and I’ve found many other uses for them. In a followup up post, I’ll work through coloring by phase and applying these contours to the phase, finally arriving at a flexible, configurable shader for domain coloring.
In the mean time, though the code is messier than I’d like, you can play with domain coloring using my Complex Function Plotter notebook.
And finally, I’d especially like to thank Observable for making this notebook possible, Mikola Lysenko for his excellent regl library, and Jacob Rus for help, thought and feedback along the way.
Feel free to leave comments on this notebook or contact me on Bluesky or on Mastodon.