Text rendering is cursed. Anyone who has worked on text will tell you the same; whether it's about layout, bi-directional, shaping, Unicode, or the rendering itself, it's never a completely solved problem. In my personal case, I've been working on trying to render text in the context of a compositing engine for creative content. I needed crazy text effects, and I needed them to be reasonably fast, which implied working with the GPU as much as possible. The distance field was an obvious requirement because it unlocks anti-aliasing and the ability to make many great effects for basically free. In this article, we will see how to compute signed distance field on the GPU because it's much faster than doing it on the CPU, especially when targeting mobile devices. We will make the algorithm decently fast, then after lamenting about the limitations, we will see what kind of effects this opens up. Progressive build of the 'あ' glyph from the Mochiy Pop One font Extraction of glyph outlines Non-bitmap fonts contain glyphs defined by outlines made of closed sequences of lines and (quadratic or cubic) Bézier curves. Extracting them isn't exactly complicated: FreeType or ttf-parser typically expose a way to do that. For the purpose of this article, we're going to hard code the list of the Bézier curves inside the shader, but of course in a more serious setup those would be uploaded through storage buffers or similar. Using this tiny program, a glyph can be dumped as series of outlines into a fixed size array: struct Bezier { vec2 p0; // start point vec2 p1; // control point 1 vec2 p2; // control point 2 vec2 p3; // end point }; #define N 42 #define NC 2 const int glyph_A_count[2] = int[](33, 9); const Bezier glyph_A[42] = Bezier[]( Bezier(vec2( 0.365370, -0.570817), vec2( 0.374708, -0.631518), vec2( 0.332685, -0.687549), vec2( 0.339689, -0.748249)), Bezier(vec2( 0.339689, -0.748249), vec2( 0.339689, -0.748249), vec2( 0.344358, -0.764591), vec2( 0.351362, -0.771595)), Bezier(vec2( 0.351362, -0.771595), vec2( 0.384047, -0.822957), vec2( 0.402724, -0.855642), vec2( 0.442412, -0.885992)), // ... ); glyph_A_count contains how many Bézier curves there is for each sub-shape composing the glyph, and glyph_A contains that list of Bézier cubic curves. Even though glyphs are also composed of lines and quadratic curves, we expand them all into cubics: "who can do more can do less". We use these formulas to respectively expand lines and quadratics into cubics: Where P_n are the Bézier control points. For simplicity and because we want to make sure the most complex case is well tested, we will stick to this approach in this article. But it also means there is a lot of room for further optimizations. Since solving linear and quadratics is much simpler, this is left as an exercise for the reader. Warning You may be tempted to upload the polynomial form of these curves directly to save some computation in the shader. Don't. You will lose the exact stitching property because one evaluated polynomial end B_n(1) will not necessary match the next polynomial start B_{n+1}(0). This makes artificial "precision holes" that will break rendering in obscure ways. Signed distance to the shape In the previous article, we saw how to get the distance to a cubic Bézier curve. Each glyph being composed of multiple outlines, we can simply run over all of them and pick the shortest distance. float get_distance(vec2 p, Bezier buf[N], int counts[NC]) { int base = 0; float dist = 1e38; for (int j = 0; j < NC; j++) { int count = counts[j]; for (int i = 0; i < count; i++) { Bezier b = buf[base + i]; float d = bezier_sq(p, b.p0, b.p1, b.p2, b.p3); dist = min(dist, d); } base += count; } return sqrt(dist); } Where bezier_sq() is the distance to the Bézier curve, squared, as defined in the previous article. Distance to the 'A' glyph from the Virgil font This works just fine, but as you can imagine, it's not cheap to solve so many distances per pixel. A first straightforward optimization would be to ignore any curve with a bounding box further than our currently best distance, because none of them can give a shorter one: Box distance optimization Where each box encloses a Bézier curve like this: Most naive/conservative bounding box We could use a tighter bound but it would require more computation so this felt like a good trade-off. Implementing this in the inner loop is pretty simple: for (int i = 0; i < count; i++) { Bezier b = buf[base + i]; vec2 p0=b.p0, p1=b.p1, p2=b.p2, p3=b.p3; // Distance to box (0 if inside), squared vec2 q0 = min(p0, min(p1, min(p2, p3))); vec2 q1 = max(p0, max(p1, max(p2, p3))); vec2 v = max(abs(q0+q1-p-p)-q1+q0, 0.)*.5; float h = dot(v,v); // We can't get a shorter distance than h if we were to compute the // distance to that curve if (h > dist) continue; float d = bezier_sq(p, p0, p1, p2, p3); dist = min(dist, d); } The distance to the bounding box formula comes from this explicative video by Inigo Quilez (the basic one, without the inside distance), adapted to the Bézier control point coordinates. This saves a lot of computation in certain cases, but the worst case is still pretty terrible, as shown by the heat map of this 'C' glyph: Heat map of how many distances are evaluated Indeed sometimes, it takes a long time to reach a good Bézier curve that is small enough to disregard most of the others. We observe it the further we go away from the beginning of the shape. So the next step is to find a good initial candidate. One cheap way to do that is to first compute the distance to the center point of each curve, and pick the smallest: // Find a good initial guess int best = 0; float boxd = 1e38; for (int i = 0; i < count; i++) { Bezier b = buf[base + i]; vec2 p0=b.p0, p1=b.p1, p2=b.p2, p3=b.p3; vec2 q0 = min(p0, min(p1, min(p2, p3))); vec2 q1 = max(p0, max(p1, max(p2, p3))); vec2 v = (q0+q1)*.5 - p; float h = dot(v,v); if (h < boxd) best=i, boxd=h; } // Initial guess Bezier bb = buf[base + best]; dist = min(dist, bezier_sq(p, bb.p0, bb.p1, bb.p2, bb.p3)); for (int i = 0; i < count; i++) { if (i == best) // We already computed this one continue; Bezier b = buf[base + i]; vec2 p0=b.p0, p1=b.p1, p2=b.p2, p3=b.p3; // ... This optimization is immediately reflected on the heat map, where only the central point seems to become a critical point (this glyph is a pathological case as it forms a circle): Heat map with a rough initial guess Winding number The last step is to figure out whether we are inside or outside the shape. There are two schools here, the even-odd and the non-zero rules. We'll pick the latter because that's the expectation in the case of font rendering. In deconstructing Bézier curves, we explained the theory of that specific algorithm so we're not going to dive into the details again. The basic idea is to strike a ray in one direction from our current position, and get how many times we cross a given curve. Here we will arbitrarily choose a horizontal ray line y = P_y where P is our current coordinate. The topology of each curve can hint us on whether it's worth considering it or not. For example, if every control point is above or below our current position, it can be ignored. We can store all the signs in a mask and bail out as soon as the ray is either completely below or completely above the bounding box of the curve: int signs = int(p0.y < p.y) | int(p1.y < p.y) << 1 | int(p2.y < p.y) << 2 | int(p3.y < p.y) << 3; if (signs == 0 || signs == 15) // all signs are identical return 0; Each sign indicates the position of the control point with regard to the ray. We can use the relative position of the starting point as a reference for the overall orientation (if there is a crossing, we know it will come from below or above): int inc = (signs & 1) == 0 ? 1 : -1; We also need to convert the Bézier curves to the usual polynomial at^3+bt^2+ct+d: vec2 a = -p0 + 3.*(p1 - p2) + p3, b = 3. * (p0 - 2.*p1 + p2), c = 3. * (p1 - p0), d = p0 - p; Then we can find the y-roots and check every point on the x-axis. For every crossing point (at most 3), we switch the sign: float r[5]; int count = root_find3(r, a.y, b.y, c.y, d.y); vec3 t = vec3(r[0], r[1], r[2]); vec3 v = ((a.x*t + b.x)*t + c.x)*t + d.x; if (count > 0 && v.x >= 0.) w += inc; if (count > 1 && v.y >= 0.) w -= inc; if (count > 2 && v.z >= 0.) w += inc; Since we already have a 5th degree root finder from the previous article, we just have to build a tiny version for the 3rd degree: int root_find3(out float r[5], float a, float b, float c, float d) { float r2[5]; int n = root_find2(r2, 3.*a, b+b, c); return cy_find5(r, r2, n, 0., 0., a, b, c, d); } Note Our root finder doesn't return roots outside [0,1] so no filtering is required. To summarize: int bezier_winding(vec2 p, vec2 p0, vec2 p1, vec2 p2, vec2 p3) { int w = 0; int signs = int(p0.y < p.y) | int(p1.y < p.y) << 1 | int(p2.y < p.y) << 2 | int(p3.y < p.y) << 3; if (signs == 0 || signs == 15) return 0; int inc = (signs & 1) == 0 ? 1 : -1; vec2 a = -p0 + 3.*(p1 - p2) + p3, b = 3. * (p0 - 2.*p1 + p2), c = 3. * (p1 - p0), d = p0 - p; float r[5]; int count = root_find3(r, a.y, b.y, c.y, d.y); vec3 t = vec3(r[0], r[1], r[2]); vec3 v = ((a.x*t + b.x)*t + c.x)*t + d.x; if (count > 0 && v.x >= 0.) w += inc; if (count > 1 && v.y >= 0.) w -= inc; if (count > 2 && v.z >= 0.) w += inc; return w; } For every sub-shape, we can accumulate the winding number, and use it at the end to decide whether we're inside or outside: float get_distance(vec2 p, Bezier buf[N], int counts[NC]) { int w = 0; int base = 0; float dist = 1e38; for (int j = 0; j < NC; j++) { int count = counts[j]; // Get the sign of the distance for (int i = 0; i < count; i++) { Bezier b = buf[base + i]; w += bezier_winding(p, b.p0, b.p1, b.p2, b.p3); } // ... } // Positive outside, negative inside return (w != 0 ? -1. : 1.) * sqrt(dist); } And voilà: Signed distance to the 'A' glyph from the Virgil font Warning This winding number logic might be too fragile: it doesn't cover potential degenerate cases such as horizontal tangents / duplicated roots. But for some reason, while I fought these issues for years, none of the weird corner cases seemed to glitch in my extensive tests, probably because the root finder is more resilient than what I was using before. Limitations Wicked curves This may look satisfying, but it's only the beginning of the problems. For example, variadic fonts are typically following chaotic patterns: The glyph 'e' in the Quicksand font In addition to the self overlapping part, notice the reverse folding triangle on the right. This completely wreck the distance field: Glyph with a broken SDF due to overlaps Even with a simple character display (meaning something that doesn't exploit the wide range of effects available with an SDF), it starts to glitch: Glitching glyph due to broken SDF Little "cracks" should appears around the overlaps. This can be mitigated by lowering the distance by a tiny constant to avoid the zero-crossing, but it impacts the overall glyph (it gets more bold). And it's not just because of variadic problem, sometimes designers rely on overlaps for simplicity: The glyph 't' in the Quicksand font And sometimes... well let's say they have a legitimate reason to do it: A Bengali glyph This is not something that can be addressed easily. For example, take these two overlapping shapes: Distance inside two overlapping shapes We see that the actual distance (white circle) is not the smallest distance to either shape, and it's not even the smallest distance to any edge: it is at an intersection point between two curves, which we do not have. Here we're dealing with line segments, but with cubic curves, the problem explodes in complexity. At this point, we need another strategy, like feeding the GPU renderer with preprocessed outline-only curves. Many people rely on curves flattening to address this issue. This is unfortunately yet another field of research that we're not going to explore this time. Inigo talked about the combination of signed distance if you want some ideas, but aside from the first one (giving up), none seems particularly applicable here. Atlas and overlapping distances Some effects such as blur or glow expand beyond the boundaries of the characters, so the distance field needs to be larger than the glyph itself. This means when an effect spread too large, there will be an overlap. If we're making an effect on a word, the distance field must be the unified version of all the word glyphs (or sometimes even the sentence). The classic approach of an atlas of glyph distances will not work reliably. In the following illustration, a geometry per glyph is used, each geometry is enlarged to account for the larger distance field, and we end up with potential overlaps when applying effects. Overlapping character geometries due to larger distance Rounded corners Like all distance maps, it suffers from the same limitations. The most common one is the rounded corners problem. This is typically addressed using a multi-channel signed distance field generator, but it's hard for me to tell how accessible it is for a portage on the GPU. msdfgen demonstration of corners improvement Note This problem only appears with intermediate textures. When computing exact distances like here directly in the shaders, this is not an issue. Effects Despite all these limitations, we can already do so much, so let's close this article on a positive note. This is not done here, but all of these effects are free as soon as we have the distance field stored in an intermediate texture. First, we have anti-aliasing / blur: AA / blur effect I wrote a dedicated article on the subject of AA (and blur) on SDF if you want more information on how to achieve that. The shape can also be drastically altered with a simple operator such as "rounding": d -= rounding; Rounding effect This is the same technique we suggested to cover up for the overlap glitch earlier, just rebranded as an effect. In the same spirit we can also create an outline stroke (on the outer edge to preserve the original glyph design): Outline effect This is sooo useful because it makes it possible for our text to be visible no matter what the background is. So many editors don't have this feature because it's hard and expensive to do correctly. Given a distance field though, all we have to do is this (which also includes anti-aliasing on every border): float aa = fwidth(d); // pixel width estimates float w = aa * .5; // half diffuse width vec2 b = vec2(0,1)*outline - d; // inner and outer boundaries; vec2(-1,0) for inner, vec2(-.5,.5) for centered float inner_mask = smoothstep(-w, w, b.x); // cut-off between the outline and the outside (whole shape w/ outline) float outer_mask = smoothstep(-w, w, b.y); // cut-off between the fill color and the outline (whole shape w/o outline) float outline_mask = outer_mask - inner_mask; vec3 o = (inner_color*inner_mask + outline_color*outline_mask) * outer_mask; We can also dig into our character with d = abs(d)-ring: Ring effect And maybe apply some glow to create a neon effect: Ring combined with a neon/glow effect float glow_power = glow * exp(-max(d, 0.) * 10.); o += glow_color * glow_power; We could also do drop shadows, all sorts of distortions, or so many other creative way exploiting this distance. You get the idea: it is fundamental as soon as you want fast visual effects. Conclusion This article is the last of the series on 2D rendering for me. I've wanted to share this experience and knowledge after many years of struggling (mostly alone) on these issues. I wish I could have succeeded in providing a good free and open-source text effects rendering engine to compete with the industry standards. (Un)fortunately for me, the adventure stops here, but I hope this will benefit creators and future tinkerers interested in the subject.

Bézier curves are a core building block of text and 2D shapes rendering. There are several approaches to rendering them, but one especially challenging problem, both mathematically and technically, is computing the distance to a Bézier curve. For quadratic curves (one control point), this is fairly accessible, but for cubic (two control points) we're going to see why it is so hard. A glyph from the Virgil font, composed of multiple Bézier curves Having this distance field opens up many rendering possibilities. It's hard, but it's possible; here is a live proof: Distance to a cubic Bézier curve In this visualization, I'm borrowing your device resources to compute the distance to the curve for every single pixel. The yellow points are the control points of the curve (in white) and the blue zone is a representation of the distance field. Note All the demos and code in this article are self-contained GLSL fragment shaders. Most of the code can be found in the article, but feel free to inspect the source code of any of these WebGL demo for the complete code. They can be run verbatim using ShaderWorkshop. The basic maths In a previous article, we explained that a Bézier curve can be expressed as a polynomial. In our case, a cubic polynomial: Where a, b, c and d are the vector coefficients derived from the start (P_0), end (P_3), and control points (P_1, P_2) using the following formulas (you can refer to the previous article for details): For a given point p in 2D space, the distance to that Bézier curve can be expressed as a length between our curve and p: Our goal is to find the t value where d(t) is the smallest. The length formula has an annoying square root, so we start with the distance squared for simplicity, which we are going to unroll: The derivative of that function will allow us to identify critical points: that is, points where the distance starts growing or reducing. Said differently, solving D'(t)=0 will identify all the maximums and minimums (we're interested in the latter) of D(t) (and thus d(t) as well). It is a bit convoluted in our case but straightforward to compute: A polynomial, this time of degree 5, emerges here. For conciseness, we can express D'(t) polynomial coefficients as a bunch of dot products: Finally, we notice that solving D'(t)=0 is equivalent to solving D'(t)/2 = 0, so we simplify the expression: Assuming we are able to solve this equation, we will get at most 5 values of t, among which we should find the shortest distance from p to the curve. Since t is bound within 0 and 1 (start and end of the curve), we will also have to test the distance at these locations. Note We could also compute the 2nd derivative in order to differentiate minimums from maximums, but simply evaluating the 5(+2) potential t values and keeping the smallest works just fine. Distance from a random point to critical testing points of the curve The red dot in the blue field is a random point in space. The red lines show which distances are evaluated (at most 5+2) to find the smallest one. Translated to GLSL code Transposing these formulas into code gives us this base template code: float bezier_distance(vec2 p, vec2 p0, vec2 p1, vec2 p2, vec2 p3) { // Start by testing the distance to the boundary points at t=0 (p0) and t=1 (p3) vec2 dp0 = p0 - p, dp3 = p3 - p; float dist = min(dot(dp0, dp0), dot(dp3, dp3)); // Bezier cubic points to polynomial coefficients vec2 a = -p0 + 3.0*(p1 - p2) + p3, b = 3.0 * (p0 - 2.0*p1 + p2), c = 3.0 * (p1 - p0), d = p0; // Solve D'(t)=0 where D(t) is the distance squared vec2 dmp = d - p; float da = 3.0 * dot(a, a), db = 5.0 * dot(a, b), dc = 4.0 * dot(a, c) + 2.0 * dot(b, b), dd = 3.0 * (dot(a, dmp) + dot(b, c)), de = 2.0 * dot(b, dmp) + dot(c, c), df = dot(c, dmp); float roots[5]; int count = root_find5(roots, da, db, dc, dd, de, df); for (int i = 0; i < count; i++) { float t = roots[i]; // Evaluate the distance to our point p and keep the smallest vec2 dp = ((a * t + b) * t + c) * t + dmp; dist = min(dist, dot(dp, dp)); } // We've been working with the squared distance so far, it's time to get its // square root return sqrt(dist); } Note dot(dp,dp) is a shorthand for the length squared, of course cheaper than computing length() which contains a square root. Warning We assume here the root finder only returns the roots that are within [0,1]. root_find5() is our 5th degree root finder, that is the function that gives us all the t (at most 5) which satisfy: But before we are able to solve that, we need to study the simpler 2nd degree polynomial solving: Solving quadratic polynomial equations Diving into the rabbit hole of solving polynomial numerically will lead you to insanity. But we still have to scratch the surface because superior degree solvers usually rely on it. My favorite quadratic root finding formula is the super simple one introduced by 3Blue1Brown, which involves locating a mid point m from which you get the 2 surrounding roots r: In GLSL, a code to cover most common corner cases would look like this: // Return true if x is not a NaN nor an infinite // highp is probably mandatory to force IEEE 754 compliance bool isfinite(highp float x) { return (floatBitsToUint(x) & 0x7f800000u) != 0x7f800000u; } // Quadratic: solve ax²+bx+c=0 int root_find2(out float r[5], float a, float b, float c) { int count = 0; float m = -b / (2.*a); float d = m*m - c/a; if (!isfinite(m) || !isfinite(d)) { // a is (probably) too small // Linear: solve bx+c=0 float s = -c / b; if (isfinite(s)) r[count++] = s; return count; } if (d < 0.) // no root return count; if (d == 0.) { r[count++] = m; // single root return count; } float z = sqrt(d); r[count++] = m - z; r[count++] = m + z; return count; } Not quite as straightforward as the math formula, isn't it? We cannot know in advance whether the division is going to succeed, so we do run divisions and only then check if they failed (and assume a reason for the failing). This is much more reliable than an arbitrary epsilon value. We also try to avoid duplicated roots. Note The roots are automatically sorted because z is always positive. Warning isfinite() may not be as reliable because in GLSL "NaNs are not required to be generated", meaning some edge case may not be supported depending on the hardware, drivers, and the current weather in Yokohama. As much as I like it, this implementation might not be the most stable numerically (even though I don't have have strong data to back this claim). Instead, we may prefer the formula from Numerical Recipes: Leading to the following alternative implementation: int root_find2(out float r[5], float a, float b, float c) { int count = 0; float d = b*b - 4.*a*c; if (d < 0.) return count; if (d == 0.) { float s = -.5 * b / a; if (isfinite(s)) r[count++] = s; return count; } float h = sqrt(d); float q = -.5 * (b + (b > 0. ? h : -h)); float r0 = q/a, r1 = c/q; if (isfinite(r0)) r[count++] = r0; if (isfinite(r1)) r[count++] = r1; return count; } This is not perfect at all (especially with the b²-4ac part). There are actually many other possible implementations, and this HAL CNRS paper shows how near impossible it is to make a correct one. It is an interesting but depressing read, especially since it "only" covers IEEE 754 floats, and we have no such guarantee on GPUs. We also don't have fma() in WebGL, which greatly limits improvements. For now, it will have to do. Solving quintic polynomial equations: attempt 1 Solving polynomials of degree 5 cannot be solved analytically like quadratics. And even if they were, we probably wouldn't do it because of numerical instability. Typically, in my experience, analytical 3rd degree polynomials solver do not provide reliable results. The first iterative algorithm I picked was the Aberth–Ehrlich method. Nowadays, more appropriate algorithms exist, but at the time I started messing up with these problems (several years ago), it was a fairly good contender. This video explores how it works. The convergence to the roots is quick, and it's overall simple to implement. But it's not without flaws. The main problem is that it works in complex space. We can't ignore the complex roots because they all "respond" to each others. And filtering these roots out at the end implies some unreliable arbitrary threshold mechanism (we keep the root only when the imaginary part is close to 0). The initialization process also annoyingly requires you to come up with a guess at what the roots are, and doesn't provide anything relevant. Aberth-Ehrlich works by refining these initial roots, similar to a more elaborate Newton iterations algorithm. Choosing better initial estimates leads to a faster convergence (meaning less iterations). The Cauchy bound specifies a space by defining the radius of a disk (complex numbers are in 2D space) where all the roots of a polynomial should lie within. We are going to use it for the initial guess, and more specifically its "tight" version (which unfortunately relies on pow()). Since Aberth-Ehrlich is a refinement and not just a shrinking process, we define and use an inner disk that has half the area of Cauchy bound disk. That way, we're more likely to start with initial guesses spread in the "middle" of the roots; this is where the √2 comes from in the formula below. #define K5_0 vec2( 0.951056516295154, 0.309016994374947) #define K5_1 vec2( 0.000000000000000, 1.000000000000000) #define K5_2 vec2(-0.951056516295154, 0.309016994374948) #define K5_3 vec2(-0.587785252292473, -0.809016994374947) #define K5_4 vec2( 0.587785252292473, -0.809016994374948) int root_find5_aberth(out float roots[5], float a, float b, float c, float d, float e, float f) { // Initial candidates set mid-way of the tight Cauchy bound estimate float r = (1.0 + max_5( pow(abs(b/a), 1.0/5.0), pow(abs(c/a), 1.0/4.0), pow(abs(d/a), 1.0/3.0), pow(abs(e/a), 1.0/2.0), abs(f/a))) / sqrt(2.0); // Spread in a circle vec2 r0 = r * K5_0, r1 = r * K5_1, r2 = r * K5_2, r3 = r * K5_3, r4 = r * K5_4; The circle constants are generated with the following script: import math import sys n = int(sys.argv[1]) for k in range(n): angle = 2 * math.pi / n off = math.pi / (2 * n) z = angle * k + off c, s = math.cos(z), math.sin(z) print(f"#define K{n}_{k} vec2({c:18.15f}, {s:18.15f})") Next, it's basically a simple iterative process. Unrolling everything for degree 5 looks like this: #define close_to_zero(x) (abs(x) < eps) // This also filters out roots out of the [0,1] range #define ADD_ROOT_IF_REAL(r) if (close_to_zero(r.y) && r.x >= 0. && r.x <= 1.) roots[count++] = r.x #define SMALL_OFF(off) (dot(off, off) <= eps*eps) /* Complex multiply, divide, inverse */ vec2 c_mul(vec2 a, vec2 b) { return mat2(a, -a.y, a.x) * b; } vec2 c_div(vec2 a, vec2 b) { return mat2(a, a.y, -a.x) * b / dot(b, b); } vec2 c_inv(vec2 z) { return vec2(z.x, -z.y) / dot(z, z); } // Compute f(x)/f'(x): complex polynomial evaluation (y) divided by their // derivatives (q) using Horner's method in one pass vec2 c_poly5d4(float a, float b, float c, float d, float e, float f, vec2 x) { vec2 y = a*x + vec2(b, 0), q = a*x + y; y = c_mul(y,x) + vec2(c, 0); q = c_mul(q,x) + y; y = c_mul(y,x) + vec2(d, 0); q = c_mul(q,x) + y; y = c_mul(y,x) + vec2(e, 0); q = c_mul(q,x) + y; y = c_mul(y,x) + vec2(f, 0); return c_div(y, q); } vec2 sum_of_inv(vec2 z0, vec2 z1, vec2 z2, vec2 z3, vec2 z4) { return c_inv(z0 - z1) + c_inv(z0 - z2) + c_inv(z0 - z3) + c_inv(z0 - z4); } int root_find5_aberth(out float roots[5], float a, float b, float c, float d, float e, float f) { if (close_to_zero(a)) return root_find4_aberth(roots, b, c, d, e, f); // Code snip: see previous snippet // float r = ... // vec2 r0, r1, r2, ... for (int m = 0; m < 16; m++) { vec2 d0 = c_poly5d4(a, b, c, d, e, f, r0), d1 = c_poly5d4(a, b, c, d, e, f, r1), d2 = c_poly5d4(a, b, c, d, e, f, r2), d3 = c_poly5d4(a, b, c, d, e, f, r3), d4 = c_poly5d4(a, b, c, d, e, f, r4); vec2 off0 = c_div(d0, vec2(1,0) - c_mul(d0, sum_of_inv(r0, r1, r2, r3, r4))), off1 = c_div(d1, vec2(1,0) - c_mul(d1, sum_of_inv(r1, r0, r2, r3, r4))), off2 = c_div(d2, vec2(1,0) - c_mul(d2, sum_of_inv(r2, r0, r1, r3, r4))), off3 = c_div(d3, vec2(1,0) - c_mul(d3, sum_of_inv(r3, r0, r1, r2, r4))), off4 = c_div(d4, vec2(1,0) - c_mul(d4, sum_of_inv(r4, r0, r1, r2, r3))); r0 -= off0; r1 -= off1; r2 -= off2; r3 -= off3; r4 -= off4; if (SMALL_OFF(off0) && SMALL_OFF(off1) && SMALL_OFF(off2) && SMALL_OFF(off3) && SMALL_OFF(off4)) break; } int count = 0; ADD_ROOT_IF_REAL(r0); ADD_ROOT_IF_REAL(r1); ADD_ROOT_IF_REAL(r2); ADD_ROOT_IF_REAL(r3); ADD_ROOT_IF_REAL(r4); return count; } When the main coefficient is too small, we fall back on the 4th degree (and so on until we reach the analytic quadratic). The 4th and 3rd degree versions of this function are easy to guess (they're pretty much identical, just removing one coefficient at each degree). We're also hardcoding a maximum of 16 iterations here because it's usually enough. To have an idea of how many iterations are required in practice, following is a visualization of the heat map of the number of iterations for every pixel: Heat map of the iterations of the Aberth-Ehrlich algorithm The big picture and the weaknesses of the algorithm should be pretty obvious by now. Among all drawbacks of this approach, there are also surprising pathological cases where the algorithm is not performing well. Fortunately, there were some progress on the state of the art in recent years. Solving quintic polynomial equations: the state of the art In 2022, Cem Yuksel published a new algorithm for polynomial root solving. Initially I had my reservations because the official implementation had a few shortcomings on some edge cases, which made me question its reliability. It's also optimized for CPU computation and is, to my very personal taste, overly complex. Fortunately, Christoph Peters showed that it was possible on the GPU by implementing it for very large degrees, and without any recursion. Inspired by that, I decided to unroll it myself for degree 5. One core difference with Aberth approach is that it is designed for arbitrary ranges. In our case this is actually convenient because, due to how Bézier curves are defined, we are only interested in roots between 0 and 1. We will need to adjust the Quadratic function to work in this range, as well as keeping the roots ordered: } float h = sqrt(d); float q = -.5 * (b + (b > 0. ? h : -h)); - float r0 = q/a, r1 = c/q; - if (isfinite(r0)) r[count++] = r0; - if (isfinite(r1)) r[count++] = r1; + vec2 v = vec2(q/a, c/q); + if (v.x > v.y) v.xy = v.yx; // keep them ordered + if (isfinite(v.x) && v.x >= 0. && v.x <= 1.) r[count++] = v.x; + if (isfinite(v.y) && v.y >= 0. && v.y <= 1.) r[count++] = v.y; return r; } The core logic of the algorithm relies on a cascade of derivatives for every degree. Christoph Peters provides an analytic formula to obtain the derivative for any degree. This is a huge helper when we need to work for an arbitrary degree, but in our case we can just differentiate manually: Since we're only interested in the roots, similar to what we did to D(t), we can simplify some of these expressions: The purpose of that cascade of derivatives is to cut the curve into monotonic segments. In practice, the core function looks like this: int root_find5_cy(out float r[5], float a, float b, float c, float d, float e, float f) { float r2[5], r3[5], r4[5]; int n = root_find2(r2, 10.*a, 4.*b, c); // degree 2 n = cy_find5(r3, r2, n, 0., 0., 10.*a, 6.*b, 3.*c, d); // degree 3 n = cy_find5(r4, r3, n, 0., 5.*a, 4.*b, 3.*c, d+d, e); // degree 4 n = cy_find5(r, r4, n, a, b, c, d, e, f); // degree 5 reutnr n; } We could unroll cy_find3, cy_find4, and cy_find5, but to keep the code simple, the degree 3 to 5 will share the same function, with leading coefficients set to 0 (hopefully the compiler does its job properly). cy_find5 relies on roots found (at most 4) at previous stages to define intervals of search: Such an approach has the nice side effect of keeping the roots ordered. The solver itself is not that complex either: float poly5(float a, float b, float c, float d, float e, float f, float t) { return ((((a * t + b) * t + c) * t + d) * t + e) * t + f; } // Quintic: solve ax⁵+bx⁴+cx³+dx²+ex+f=0 iint cy_find5(out float r[5], float r4[5], int n, float a, float b, float c, float d, float e, float f) { int count = 0; vec2 p = vec2(0, poly5(a,b,c,d,e,f, 0.)); for (int i = 0; i <= n; i++) { float x = i == n ? 1. : r4[i], y = poly5(a,b,c,d,e,f, x); if (p.y * y > 0.) continue; float v = bisect5(a,b,c,d,e,f, vec2(p.x,x), vec2(p.y,y)); r[count++] = v; p = vec2(x, y); } return count; } The last brick of the algorithm is the Newton bisection, the slowest part: // Newton bisection // // a,b,c,d,e,f: 5th degree polynomial parameters // t: x-axis boundaries // v: respectively f(t.x) and f(t.y) float bisect5(float a, float b, float c, float d, float e, float f, vec2 t, vec2 v) { float x = (t.x+t.y) * .5; // mid point float s = v.x < v.y ? 1. : -1.; // sign flip for (int i = 0; i < 32; i++) { // Evaluate polynomial (y) and its derivative (q) using Horner's method in one pass float y = a*x + b, q = a*x + y; y = y*x + c; q = q*x + y; y = y*x + d; q = q*x + y; y = y*x + e; q = q*x + y; y = y*x + f; t = s*y < 0. ? vec2(x, t.y) : vec2(t.x, x); float next = x - y/q; // Newton iteration next = next >= t.x && next <= t.y ? next : (t.x+t.y) * .5; if (abs(next - x) < eps) return next; x = next; } return x; } And that's pretty much it. Looking at its heat map, it has a completely different look than Aberth: Heat map of the iterations of Cem Yuksel algorithm The number of iterations might be larger but it is much faster (I observed a factor 3 on my machine), the code is shorter, and actually more reliable. Note The scale used to represent the heat map is not the same as the one used in Aberth, but it is identical with the method presented in the next section. Exploring ITP convergence The bisection being the hot loop, it is interesting to ponder how to make this faster. A while back, Raph Levien hypothesized about how the ITP method could perform. Out of curiosity, I gave it a chance. The function is designed to work like a bisection, claiming to be as performant in the worst case. There isn't a lot of code, and the paper provides a pseudo-code. But implementing it was actually challenging in many ways. First of all, the authors didn't seem to find relevant to mention that it only works if f(a)<0<f(b). If f(a)>0>f(b), you're pretty much on your own. It requires just 2 lines of adjustments but figuring out this shortcoming of the algorithm was particularly unexpected. Another bothering aspect concerns the parameters: K_1, K_2 and n_0. The paper proposes those: I played with them for a while and couldn't find any set that would make a real difference, so I ended up with the following: For performance reasons, reducing K_2 to a value of 2 saves a call to pow(). For K_1, CRAN seems to suggest \frac{0.2}{b-a} so I went along with it And for n_0, well 1 or 2 seem to be the usual parameter. In the end, the function looks like this: // ITP algorithm (2020) by Oliveira & Takahashi // "An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality" // // a,b,c,d,e,f: 5th degree polynomial parameters // t: x-axis boundaries (a and b in the paper) // v: respectively f(a) and f(b) in the paper (evaluation of the function with t.x and t.y) float itp5(float a, float b, float c, float d, float e, float f, vec2 t, vec2 v) { float diff = t.y-t.x; // K1 and n0 suggested by CRAN float K1 = .2 / diff; int n0 = 1; // The paper has the assumption that f(a)<0<f(b) but we want to // support f(a)>0>f(b) too, so we keep a sign flip float s = v.x < v.y ? 1. : -1.; // Using log(ab)=log(a)+log(b): log2(x/(2ε)) <=> log2(x/ε)-1 int nh = int(ceil(log2(diff/eps)-1.)); // n_{1/2} (half point) int n_max = nh + n0; // ε 2^(n_max-k) = ε 2^n_max 2^-k = ε 2^n_max ½^k // ½^k is done iteratively in the loop, simplifying the arithmetic float q = eps * float(1<<n_max); while (diff > eps+eps) { // Interpolation float xf = (v.y*t.x - v.x*t.y) / (v.y-v.x); // Regula-Falsi // Truncation float xh = (t.x+t.y) * .5; // x half point float sigma = sign(xh - xf); float delta = K1*diff*diff; // save a pow() by forcing K2=2 float xt = delta <= abs(xh - xf) ? xf + sigma*delta : xh; // xt: truncation of xf // Projection float r = q - diff*.5; float x = abs(xt-xh) <= r ? xt : xh-sigma*r; // Updating float y = poly5(a,b,c,d,e,f, x); float side = s*y; if (side > 0.) t.y=x, v.y=y; else if (side < 0.) t.x=x, v.x=y; else return x; diff = t.y-t.x; q *= .5; } return (t.x+t.y) * .5; } This function can be used as a drop'in replacement for bisect5. I had a lot of expectations about it, but in the end it requires more iterations than the bisection we implemented. The paper claims to perform at least as good as a bisection, but our bisect5 is driven by the derivatives so it converges much faster. Here is the heat map with itp5 instead of bisect5: Heat map of the iterations of Cem Yuksel algorithm with ITP method Conclusion The naive unrolled version of Cem Yuksel paper definitely is, so far, the best choice for our problem. I have still concerns about how to implement a good quadratic formula, and I have my reservations about various edge cases. There is also still room for improvements in the cubic solver (degree 3) because it's still a special case where analytical formulas exist, but in general this implementation is satisfying. The next step is to work with chains of Bézier curves to make up complex shapes (such as font glyphs). It will lead us to build a signed distance field. This is not trivial at all and mandates one or several dedicated articles. We will hopefully study these subjects in the not-so-distant future.

A few weeks ago, I made a tiny demo that fits into 448 characters: Red Alp GLSL demo in 448 characters void main(){vec3 c,p,K=vec3(3,1,0);for(float z,i,a,g=1.,t,h,d,w,k=.15;i++<1e2;d=max(max(d-3. ,-d),a=z)*k,w=g-g/exp(h>.001?a++,d/.4:h*3e2),g-=a*=w,c+=a*d*4.5+(d>z?z:h/2e2)*K,a=min(p.y+2. ,1.),c.r+=w*a*a*.1,t+=min(h*.2,k/=.985))for(p=normalize(vec3(P+P-R,R.y))*t,p.xz*=mat2(cos( sin(T*.2)+K.zyxz*11.)),p.z+=T*.3,d=p.y,h=d+.5,a=.01;a<1.;a+=a)p.xz*=mat2(8,6,-6,8)*.1,d+=abs (dot(sin((p/a+T)*.3),p-p+a)),h+=abs(dot(sin(p.xz*.6/a),P-P+a));O=vec4(tanh(c),1);} Note The number of characters was 464 characters at first, but thanks to the community it got reduced further, and the article updated accordingly. There is no texture, no mesh, no 3D helper: it's simply a procedural mathematical formula evaluated at each pixel assigning them a color. Code golfing is about making it as short as possible, and thus is part of the art performance. To put things into perspective, the 853x480 JPEG thumbnail of this article is 167x larger than this code. You can watch a larger version on its main dedicated page, or a portage on Shadertoy (484 chars). If your device is not powerful enough (I'm sorry for the lag on this page) or doesn't support WebGL2, a short preview video can be seen on Mastodon. I'm guessing the wizardry of the code has confused many people so we're going to dive through the making-of together. Overall, this demo is a particularly dense and entangled compilation of different techniques, where each aspect could mandate a dedicated article. For that reason, some parts will prefer to link to external resources when the literacy is verbose on the subject. Warning Some demos in this article will start "decaying" over time due to floating point variables getting too large. Reloading the page should fix that. The base template The code is written in GLSL and is executed for each pixel (technically each fragment) on a simple quad geometry (to be accurate it's even a single big triangle). There is no geometry aside from that, it's basically just a fragment shader. The fragment receives 3 different inputs: the canvas resolution vec2 R the time float T the pixel position vec2 P (basically gl_FragCoord.xy) And it has to output a sRGB color in out vec3 O. The code has to be written in a void main() function, and that's pretty much all we need to start. If you're curious about the glue to setup WebGL2, just look at the source code on the dedicated page. There is no external dependency and the canvas setup code is pretty simple. Development setup For development, people usually directly use Shadertoy. I prefer to use my own local live coding environment: ShaderWorkshop. It can be run without setting up anything, just uv run --with shader-workshop sw-server (assuming the uv Python package manager is installed on the machine). Aside from the comfort of being able to use your favorite code editor, it allows instancing live controls for uniforms very easily, making it smooth to interact with any value and get an immediate feedback. Red Alp demo with user controls as seen from ShaderWorkshop Noise One of the core primitive we need is a noise function: it is required for the mountains, the fog, and the clouds. In a recent article, I talked about gradient noise. We could technically use that, but it will have a lot of drawbacks. First of all, it's super expensive. I know because I made a demo using it the other day, and it was awfully slow. Once per pixel would be fine, but in our case it will have to be evaluated a hundred times, so we need something faster. Secondly, we're trying to make it as short as possible, and the 2D gradient noise, even minified, is already twice as big as the size of the full demo. We will also need a 3D noise for the clouds and fog, which is even larger and more expensive. And that's not even accounting for the fbm signal combination code. Inigo Quilez, in his famous Rainforest, used value noise. It is faster, but it still won't do for us for the same reasons, just somehow mitigated. And since we're professionals, we're not going to cheat by sampling a noise texture. Fortunately, while reverse engineering some Shadertoy demos, in particular the ones from diatribes, I came across some code that made use of this incredible technique of accumulating sine waves. Combining sin waves Let's say we want to combine two sine waves in order to get a height map as a 3rd dimension. There are multiple ways of achieving that. For example, we can multiply them: But we could also add them together: The surprising take here is that... it's pretty much equivalent. It doesn't give the same result for sure, but visually it could be considered the same, just with a frequency and amplitude a bit different, and rotated on the z axis by 45°. Similarly, you may think using cosines instead of sinusoids would make a difference, but no, even when combined together, they always give the same base pattern we just saw. So let's pick one, let's say z=\sin x + \sin y. But this time, we're going to take the absolute value to transform the up and down pattern into bumps: These bumps are the perfect base for clouds, but not so much for spiky mountains going through aggressive erosion. But with the help of this weird little trick, we can just flip the shape upside down to get sharp edges: We now have the basis for both our clouds and mountains, but it's not yet convincing. The next step is to use the fbm loop as if we were dealing with Gaussian or value noise: we accumulate several frequencies of our signal together: S is the sign (-1 for spiky, 1 for bobby) i is the octave identifier going from 0 to N-1 (included). F(x,y) is usually the noise signal function, in our case it's the sinusoid combination function, we choose |\sin x + \sin y| here. l is the lacunarity factor, that is how frequency changes at each octaves; this is usually a multiply by 2 or a close value. g is the gain, that is by how the amplitude changes at each octaves; this is usually a multiply by 0.5 or a close value. Without surprise this is still very periodic, but we can see a glimpse of chaos emerging. The final touch does all the magic: all we have to do now is simply rotate each layer by like, 30° or something (I'll pick 0.5 radians here, or about 29°): The symmetry around the origin is still noticeable, but the illusion will work as we will move away from it. It's also possible to add some phase or offsetting (arbitrary addition within the sin or between each layer). I implemented this in a Desmos 3D scene with all the parameters if one wants to play with it. The formula there has a few more controls, for example the vertical location, an optional transition offset in addition to the rotation, and controls for the base frequency and amplitude. If this mathematical gibberish is above your head, a GLSL code for the 2D noise could look like this with a lacunarity of 2, a gain of 0.5 and 5 octaves: float noise(vec2 p) { float v = 0.0; float amplitude = 1.0; for (int i = 0; i < 5; i++) { p = rotate(0.5) * p; // rotate our space (more on this in the next section) v += abs(sin(p.x) + sin(p.y)) * amplitude; // accumulate noise p *= 2.0; // double the frequency at each octave amplitude *= 0.5; // half the amplitude at each octave } return v; } One cool trick here: abs(sin(p.x)+sin(p.y)) could also be written abs(dot(sin(p),vec2(1))). This is interesting because now we can operate on the two components of p, easing the possibility to modify them at once (for example doing p*A+B). The dot trick doesn't work with sin(p.x)*sin(p.y), but fortunately, as we saw before, multiply and addition are similar and could be swapped in various situations. Rotations We needed some rotation for the noise, and they will be required again soon, so we need to have a closer look to them. Let's start with the formula most people are familiar with: A matrix can be seen as a function, so mathematically writing p'=M \cdot p would be equivalent to the code p=rotate(angle)*p with: // Matrix for a counter-clockwise rotation mat2 rotate(float a) { return mat2( cos(a), sin(a), // column 1 -sin(a), cos(a) // column 2 ); } Doing p'=M \cdot p is rotating the space p lies into, which means it gives the illusion the object is rotating clockwise. Though, in the expression p=rotate(angle)*p, I can't help but be bothered by the redundancy of p, so I would prefer to write p*=rotate(angle) instead. Since matrices are not commutative, this will instead do a counter-clockwise rotation of the object. The inlined rotation ends up being: p *= mat2(cos(a),sin(a),-sin(a),cos(a)); // counter-clockwise rotation of object at point p Note To make the rotation clockwise, we can of course use -a, or we can transpose the matrix: mat2(cos(a),-sin(a),sin(a),cos(a)). This is problematic though: we need to repeat the angle 4 times, which can be particularly troublesome if we want to create a macro and/or don't want an intermediate variable for the angle. But I got you covered: trigonometry has a shitton of identities, and we can express every sin according to a cos (and the other way around). For example, here is another formulation of the same expression: p *= mat2(cos(a + vec4(0,3,1,0)*PI/2.0)); Now the angle appears only once, in a vectorized cosine call. GLSL has degrees() and radians() functions, but it doesn't expose anything for \pi nor \tau constants. And of course, it doesn't have sinpi and cospi implementation either. So it's obvious they want us to use \arccos(-1) for \pi and \arccos(0) for \pi/2: p *= mat2(cos(a + vec4(0,3,1,0)*acos(0.))); Note To specify a as a normalized value, we can use mat2(cos((a*4.+vec4(0,3,1,0))*acos(0.))). On his Unofficial Shadertoy blog, Fabrice Neyret goes further and provide us with a very cute approximation, which is the one we will use: p *= mat2(cos(a + vec4(0,11,33,0))); I checked for the best numbers in 2 digits, and I can confirm they are indeed the ones providing the best accuracy. Comparison of the 2 rotations matrices On this last figure, the slight red/green on the outline of the circle represents the loss of precision. Note With 3 digits, 344 and 699 can respectively be used instead of 11 and 33. This is good when we want a dynamic rotation angle (we will need that for the camera panning typically), but sometimes we just need a hardcoded value: for example in the rotate(0.5) of our combined noise function. mat2(cos(.5+vec4(0,11,33,0))) is fine but we can do better. Through Inigo's demos I found the following: mat2(.8,.6,-.6,.8). It makes a rotation angle of about 37° (around 0.64 radians) in a very tiny form. Since 0.5 was pretty much arbitrary, we can just use this matrix as well. And we can make it even smaller (thank you jolle): p *= mat2(8,6,-6,8)*.1; // rotate p counter-clockwise by about 37° without any trigo One last rotation tip from Fabrice's bag of tricks: rotating in 3D around an axis can be done with the help of GLSL swizzling: p.xz *= rotate(0.5); // 3D rotation around y-axis (the absent component) We will use this too. Note p.zy *= rotate(.5) is the same p.yz *= rotate(-.5), if we need to save one character and can't transpose the matrix. Camera (and axis) setup One last essential before going creative is the camera setup. We start with the 2D P pixel coordinates which we are going to make resolution independent by transforming them into a traditional mathematical coordinates system: // 1:1 ratio with [-1,1] along the shortest axis (horizontal or vertical) vec2 u = (2.0*P - R) / min(R.x, R.y); Since we know our demo will be rendered in landscape mode, dividing by R.y is enough. We can also save one character using P+P: // 1:1 ratio with [-1,1] along the vertical axis vec2 u = (P+P - R) / R.y; To enter 3D space, we append a third component, giving us either a right or a left-handed Y-up coordinates system. This choice is not completely random. Indeed, it's easier/shorter to add a 3rd dimension at the end compared to interleaving a middle component. Compare the length of vec3(P, z) to vec3(P.x, z, P.y) (Z-up convention). In the former case, picking just a plane remains short and easy thanks to swizzling: p.xz instead of p.xy. To work in 3D, we need an origin point (ro for ray origin) and a looking direction (rd for ray direction). ro is picked arbitrarily for the eye position, while rd is usually calculated thanks to a lookAt helper: // Right-hand with Y-up (like Godot) mat3 lookAt(vec3 origin /* where we are */, vec3 target /* where we look */) { vec3 w = normalize(target - origin); vec3 u = normalize(cross(w, vec3(0,1,0))); vec3 v = normalize(cross(u, w)); // Note: normalize() can be ditched here return mat3(u, v, w); } Right-hand Y-up 3D coordinates system Which is then used like that, for example: vec2 u = (P+P - R) / R.y; vec3 target = /* ... */; vec3 ro = /* ... */; vec3 rd = normalize(lookAt(ro, target) * vec3(u, 1)); Note I made a Shadertoy demo to experiment with different 3D coordinate spaces if you are interested in digging this further. All of this is perfectly fine because it is flexible, but it's also way too much unnecessary code for our needs, so we need to shrink it. One approach is to pick a simple origin and straight target point so that the matrix is as simple as possible. And then later on apply some transformations on the point. If we give ro=vec3(0) and target=vec3(0,0,1), we end up with an identity matrix, so we can ditch everything and just write: vec3 rd = normalize(vec3((P+P - R) / R.y, 1)); This can be shorten further: since the vector is normalized anyway, we can scale it at will, for example by a factor R.y, saving us precious characters: vec3 rd = normalize(vec3(P+P - R, R.y)); And just like that, we are located at the origin vec3(0), looking toward Z+, ready to render our scene. Mountain height map It's finally time to build our scene. We're going to start with our noise function previously defined, but we're going tweak it in various ways to craft a mountain height map function. Here is our first draft: const float mountain_y = -0.5; // mountain y-axis position const float mountain_f = 0.6; // mountain base frequency float mountain_height_map(vec2 p) { float h = mountain_y; for (float a = 1.0; a > 0.01; a /= 2.0) { p *= rotate(0.5); h += abs(dot(sin(p*mountain_f / a), vec2(1))) * a; // dot(sin(v),1) -> sin(v.x)+sin(v.y) } return -h; // minus for the spiky version of the noise } We're exploiting one important correlation of the noise function: at every octave, the amplitude is halving while the frequency is doubling. So instead of having 2 running variables, we just have an amplitude a getting halved every octave, and we divide our position p by a (which is the same as multiplying by a frequency that doubles itself). I actually like this way of writing the loop because we can stop the loop when the amplitude is meaningless (a>0.01 acts as a precision stopper). Unfortunately, we'll have to change it to save one character: a/=2. is too long for the iteration, we're going to double instead by using a+=a which saves one character. So instead the loop will be written the other way around: for (float a=.01; a<1.; a+=a). It's not exactly equivalent, but it's good enough (and we can still tweak the values if necessary). We're going to inline the constants and rotate, and use one more cool trick: vec2(1) can be shortened: we just need another vec2. Luckily we have p, so we can simply replace it with p/p. Finally, we can get rid of the braces of the for loop by using the , in its local scope: float mountain_height_map(vec2 p) { float h = -.5; for (float a=.01; a<1.; a+=a) p *= mat2(8,6,-6,8)*.1, h += abs(dot(sin(p*.6/a), p/p))*a; return -h; } p/p works fine as long as it's not zero. In this particular case, we can instead use vec2(0) (obtained with p-p) and then include the a amplitude multiplier within the expression: abs(dot(sin(p*.6/a), p-p+a)). (p-p+a is the same as vec2(a) when p is a vec2). We end up with the following safer version: float mountain_height_map(vec2 p) { float h = -.5; for (float a=.01; a<1.; a+=a) p *= mat2(8,6,-6,8)*.1, h += abs(dot(sin(p*.6/a), p-p+a)); return -h; } Mountain height map in 2D (rescaled for display) To render this in 3D, we are going to do some ray-marching. Solid ray-marching The main technique used in most Shadertoy demos is ray-marching. I will assume familiarity with the technique, but if that's not the case, An introduction to Raymarching (YouTube) by kishimisu and Painting with Math: A Gentle Study of Raymarching by Maxime Heckel were good resources for me. In short: we start from a position in space called the ray origin ro and we project it toward a ray direction rd. At every iteration we check the distance to the closest solid in our scene, and step toward that distance, hoping to converge closer and closer to the object boundary. We end up with this main loop template: float t = 0.0; vec3 ro = vec3(0); // ray origin vec3 rd = normalize(vec3(P+P - R, R.y)); // ray direction // 100 iterations should be enough to hit something if there is any for (int i = 0; i < 100; i++) { vec3 p = ro + rd*t; // t amount in rd direction from ro origin float h = distance_to_solid(p); // 3D distance function if (h < 0.001) { // we converged close enough to a solid // Here we assign a color according to where p is // [...] break; } t += h; // there is no solid closer than h so we step by that much } This works fine for solids expressed with 3D distance fields, that is functions that for a given point give the distance to the object. We will use it for our mountain, with one subtlety: the noise height map of the mountain is not exactly a distance (it is only the distance to what's below our current point p): float distance_to_solid(vec3 p) { // positive outside, negative inside return p.y - mountain_height_map(p.xz); } Because of this, we can't step by the distance directly, or we're likely to go through mountains during the stepping (t += h). A common workaround here is to step a certain percentage of that distance to play it safe. Technically we should figure out the theorical proper shrink factor, but we're going to take a shortcut today and just arbitrarily cut. Using trial and error I ended up with 20% of the distance. After a few simplifications, we end up with the following (complete) code: float mountain_height_map(vec2 p) { float h = .5; for (float a=.01; a<1.; a+=a) p *= mat2(8,6,-6,8)*.1, h += abs(dot(sin(p*.6/a), p-p+a)); return -h; } float distance_to_solid(vec3 p) { return p.y - mountain_height_map(p.xz); } void main() { vec3 rd = normalize(vec3(P+P - R, R.y)); float t = 0.0, color = 0.0; for (int i = 0; i < 100; i++) { vec3 p = rd*t; p.z += T*.2; // move forward float h = distance_to_solid(p); if (h < 0.001) { color = exp(-t*t*.01); // depth map like "coloring" break; } t += h * 0.2; } O = vec4(vec3(pow(color, 3.0/2.2)), 1); } Basic ray-marching of the mountain height map We start at ro=vec3(0) so I dropped the variable entirely. You may be curious about the power at the end; this is just a combination of luminance perception with gamma 2.2 (sRGB) transfer function. It only works well for grayscale; for more information, see my previous article on blending. Clouds and fog Compared to the mountain, the clouds and fog will need a 3 dimensional noise. Well, we don't need to be very original here; we simply extend the 2D noise to 3D: float noise3(vec3 p) { float v; for (float a=.01; a<1.; a+=a) p.xz *= mat2(8,6,-6,8)*.1, v += abs(dot(sin(p*.3/a + T*.3), p-p+a)); return v; } The base frequency is lowered to 0.3 to make it smoother, and the p goes from 2 to 3 dimensions. Notice how the rotation is only done on the y-axis, the one pointing up): don't worry, it's good enough for our purpose. We also add a phase (meaning we are offsetting the sinusoid) of T*0.3 (T is the time in seconds, slowed down by the multiply) to slowly morph it over time. The base frequency and time scale being identical is a happy "coincidence" to be factored out later (I actually forgot about it until jolle reminded me of it). You also most definitely noticed v isn't explicitly initialized: while only true WebGL, it guarantees zero initialization so we're saving a few characters here. Volumetric ray-marching For volumetric material (clouds and fog), the loop is a bit different: instead of calculating the distance to the solid for our current point p, we do compute the density of our target "object". Funny enough, it can be thought as a 3D SDF but with the sign flipped: positive inside (because the density increases as we go deeper) and negative outside (there is no density, we're not in it). const float clouds_y = 3.0; // vertical position float clouds_density(vec3 p) { float n = noise3(p); // random value associated with a 3D position in space float h = -clouds_y + n; // similar to mountain_height_map() but 3d and bobby float d = p.y - h; // similar to distance_to_solid() d = -d; // flip sign: distance to density // We are only interested in the density within the material, // the density will be considered 0 when outside of it. return max(d, 0.0); } For simplicity, we're going to rewrite the function like this: const float clouds_y = 3.0; float clouds_density(vec3 p) { float n = noise3(p); float d = -p.y - cloud_y + n; return max(d, 0.0); } Compared to the solid ray-marching loop, the volumetric one doesn't bail out when it reaches the target. Instead, it slowly steps into it, damping the light as the density increases: const float absorption = 0.15; const float radiance = 1.0; void main() { float step_len = 0.15; float t; vec3 rd = normalize(vec3(P+P-R,R.y)); vec3 color; float transmittance = 1.0; // remaining visibility for (int i = 0; i < 100; i++) { vec3 p = rd*t; // Move camera forward p.z += T * 1.5; // How many particules of the material we can find at that position // If negative, we're not in the element yet, otherwise it's the density // (getting higher as we go deeper into it typically). float d = clouds_density(p); // Integrate the density discretely: we assume the segment of length // we're walking has the same point density all along d *= step_len; // The fraction of light that survives through this segment (Beer-Lambert law) // The denser, the closer to 0 this gets float attenuation = exp(-d*absorption); float emission = d*radiance; // how much light is emitted along the segment (glow) float alpha = 1.0 - attenuation; // fraction of light removed for that given density segment float weight = alpha * transmittance; // Accumulate color emission color += weight * emission; transmittance -= weight; // could also be written transmittance *= attenuation // Normal volumetric marching (step_len) clamped to the distance to the // solid (mountain) t += step_len; // Larger volumetric steps as we go far step_len *= 1.015; } O = vec4(pow(color, vec3(3.0/2.2)), 1); } The core idea is that the volumetric material emit some radiance but also absorbs the atmospheric light. The deeper we get, the smaller the transmittance gets, til it converges to 0 and stops all light. All the threshold you see are chosen by tweaking them through trial and error, not any particular logic. It is also highly dependent on the total number of iterations. Note Steps get larger and larger as the distance increases; this is because we don't need as much precision per "slice", but we still want to reach a long distance. Basic volumetric ray-marching of the clouds density map We want to be positioned below the clouds, so we're going to need a simple sign flip in the function. The fog will take the place at the bottom, except upside down (the sharpness will give a mountain-hug feeling) and at a different position. clouds_density() becomes: const float clouds_y = 3.0; const float fog_y = 0.0; float clouds_fog_density(vec3 p) { float n = noise3(p); float clouds_d = p.y - clouds_y + n; float fog_d = p.y - fog_y + n; // Pick the element with the highest density (they don't overlap anyway) float d = max(clouds_d, -fog_d); return max(d, 0.0); } Both clouds and fog with volumetric ray-marching For more resources on volumetric rendering, following are the ones I studied the most: Volumetric Rendering in 2 parts, by Chris Real-time dreamy Cloudscapes with Volumetric Raymarching, by Maxime Heckel again Volumetric Raymarching, by Xor Combining ray-marching Having a single ray-marching loop combining the two methods (solid and volumetric) can be challenging. In theory, we should stop the marching when we hit a solid, bail out of the loop, do some fancy normal calculations along with light position. We can't afford any of that, so we're going to start doing art from now on. We start from the volumetric ray-marching loop, and add the distance to the mountain: for (int i = 0; i < 100; i++) { vec3 p = rd*t; // ... float d = clouds_fog_density(p); float h = distance_to_solid(p); // ... } If h gets small enough, we can assume we hit a solid: bool solid = h < 0.001; In volumetric, the attenuation is calculated with the Beer-Lambert law. For solid, we're simply going to make it fairly high: - float attenuation = exp(-d*absorption); + float attenuation = solid ? 0.95 : exp(-d*absorption); This has the effect of making the mountain like a very dense gas. We're also going to disable the light emission from the solid (it will be handled differently down the line): - float emission = d*radiance; + float emission = solid ? 0.0 : d*radiance; The transmittance is not going to be changed when we hit a solid as we just want to accumulate light onto it: - transmittance -= weight; + if (!solid) transmittance -= weight; Finally, we have to combine the volumetric stepping (t += step_len) with the solid stepping (t += h*0.2) by choosing the safest step length, that is the minimum: - t += step_len; + t += min(h*0.2, step_len); We end up with the following: Combination of volumetric and solid ray-marching We can notice the mountain from negative space and the discrete presence of the fog, but it's definitely way too dark. So the first thing we're going to do is boost the radiance, as well as the absorption for the contrast: -const float absorption = 0.15; -const float radiance = 1.0; +const float absorption = 2.5; +const float radiance = 4.5; This will make the light actually overshoot, so we also have to replace the current gamma 2.2 correction with a cheap and simple tone mapping hack: tanh(): - O = vec4(pow(color, vec3(3.0/2.2)), 1); + O = vec4(tanh(color), 1); Tonemapping the scene The clouds and fog are much better but the mountain is still trying to act cool. So we're going to tweak it in the loop: emission += 0.1; This boosts the overall emission. While at it, since the horizon is also sadly dark, we want to blast some light into it: color += d == 0.0 ? 0.005*h : 0.0; mkbosmans from HackerNews noticed that the opposite of d==0.0 is actually d>0.0 due to the max(...,0). So we could write it more simply: color += d > 0.0 ? 0.0 : 0.005*h; When the density is null (meaning we're outside clouds and fog), an additional light is added, proportional to how far we are from any solid (the sky gets the most boost basically). More atmospheric light The mountain looks fine but I wanted a more eerie atmosphere, so I changed the attenuation: - float attenuation = solid ? 0.95 : exp(-d*absorption); + float attenuation = exp(solid ? -h*300.0 : -d*absorption); Now instead of being a hard value, the attenuation is correlated with the proximity to the solid (when getting close to it). This has nothing to with any physics formula or anything, it's more of an implementation trick which relies on the ray-marching algorithm. The effect it creates is those crack-like polygon edges on the mountain. To add more to the effect, the emission boost is tweaked into: - emission += 0.1; + float e = min(p.y - mountain_y + 1.5, 1.0); + emission += e*e * 0.1; This makes the bottom of the mountain darker quadratically: only the tip of the mountain would have the glowing cracks. Making mountains eerie Color We've been working in grayscale so far, which is a usually a sound approach to visual art in general. But we can afford a few more characters to move the scene to a decent piece of art from the 21st century. Adding the color just requires very tiny changes. First, the emission boost is going to target only the red component of the color: - emission += e*e * 0.1; float alpha = 1. - attenuation; float weight = alpha * transmittance; color += weight * emission; + color.r += weight * e*e * 0.1; And similarly, the overall addition of light into the horizon/atmosphere is going to get a redish/orange tint: - color += d > 0.0 ? 0.0 : 0.005*h; + color += (d > 0.0 ? 0.0 : 0.005*h) * vec3(3,1,0); Add a red/orange tint Last tweaks We're almost done. For the last tweak, we're going to add a cyclic panning rotation of the camera, and adjust the moving speed: p.xz *= mat2(cos(sin(T*.2)+vec4(0,11,33,0))); p.z += T*.3; Note I'm currently satisfied with the "seed" of the scene, but otherwise it would have been possible to nudge the noise in different ways. For example, remember the sin can be replaced with cos in either or both volumetric and mountain related noises. Similarly, the offsetting +T could be changed into -T for a different morphing effect. And of course the rotations can be swapped (either by changing .xz into .zx or transposing the values). Code golfing At this point, our code went through early stages of code golfing, but it still needs some work to reach perfection. Stripped out of its comments, it looks like this: // Reference code: 1278 chars (unnecessary spaces and line breaks are not counted) const float fog_y = 0.0; const float clouds_y = 3.0; const float mountain_y = -0.5; const float absorption = 2.5; const float radiance = 4.5; float noise3(vec3 p) { float v; for(float a=.01; a<1.; a+=a) p.xz *= mat2(8,6,-6,8)*.1, v += abs(dot(sin(p*.3/a + T*.3), vec3(1)))*a; return v; } float clouds_fog_density(vec3 p) { float n = noise3(p); float clouds_d = p.y-clouds_y+n; float fog_d = p.y-fog_y+n; float d = max(clouds_d, -fog_d); return max(d, 0.0); } float mountain_height_map(vec2 p) { float h = -mountain_y; for (float a=.01; a<1.; a+=a) p *= mat2(8,6,-6,8)*.1, h += abs(dot(sin(p*.6/a), vec2(1)))*a; return -h; } float distance_to_solid(vec3 p) { return p.y - mountain_height_map(p.xz); } void main() { float step_len = 0.15; float t; vec3 color; float transmittance = 1.0; vec3 rd = normalize(vec3(P+P-R,R.y)); for (int i = 0; i < 100; i++) { vec3 p = rd*t; p.xz *= mat2(cos(sin(T*.2)+vec4(0,11,33,0))); p.z += T*.3; float d = clouds_fog_density(p); float h = distance_to_solid(p); bool solid = h < 0.001; d *= step_len; float attenuation = exp(solid ? -h*300.0 : -d*absorption); float emission = solid ? 0.0 : d*radiance; float e = min(p.y - mountain_y + 1.5, 1.0); float alpha = 1. - attenuation; float weight = alpha * transmittance; color += weight * emission; color.r += weight * e*e * 0.1; color += (d > 0.0 ? 0.0 : 0.005*h) * vec3(3,1,0); if (!solid) transmittance -= weight; t += min(h*0.2, step_len); step_len *= 1.015; } O = vec4(tanh(color), 1); } The first thing we're going to do is notice that both the mountain, clouds, and fog use the exact same loop. Factoring them out and inlining the whole thing in the main function is the obvious move: // 922 chars const float fog_y = 0.0; const float clouds_y = 3.0; const float mountain_y = -0.5; const float absorption = 2.5; const float radiance = 4.5; void main() { float step_len = 0.15; float t; vec3 color; float transmittance = 1.0; vec3 rd = normalize(vec3(P+P-R,R.y)); for (int i = 0; i < 100; i++) { vec3 p = rd*t; p.xz *= mat2(cos(sin(T*.2)+vec4(0,11,33,0))); p.z += T*.3; float d = p.y; float h = p.y-mountain_y; for (float a=.01; a<1.; a+=a) p.xz *= mat2(8,6,-6,8)*.1, d += abs(dot(sin(p*.3/a + T*.3), vec3(1)))*a, h += abs(dot(sin(p.xz*.6/a), vec2(1)))*a; d = max(max(d-clouds_y, -(d-fog_y)), 0.0); bool solid = h < 0.001; d *= step_len; float attenuation = exp(solid ? -h*300.0 : -d*absorption); float emission = solid ? 0.0 : d*radiance; float e = min(p.y - mountain_y + 1.5, 1.0); float alpha = 1. - attenuation; float weight = alpha * transmittance; color += weight * emission; color.r += weight * e*e * 0.1; color += (d > 0.0 ? 0.0 : 0.005*h) * vec3(3,1,0); if (!solid) transmittance -= weight; t += min(h*0.2, step_len); step_len *= 1.015; } O = vec4(tanh(color), 1); } Next, we are going to do the following changes: Rename every variable to single letter or inline them whenever possible Inline all constants Remove any explicit zero initialization Use float instead of int for the iterator and bool for the solid flag Pack all float and vec3 declarations together Simplify numbers: 1e2 instead of 100.0, 3. instead of 3.0, etc. vec*() constructor act like cast, so you can pass down integers Instead of *x, /(1/x) is sometimes shorter (for example /.4 instead of *2.5) (thanks coyote) // 491 chars void main() { vec3 c, p; for (float i, a, g=1., t, h, d, w, k=.15, x, e; i < 1e2; i++) { p = normalize(vec3(P+P-R,R.y))*t; p.xz *= mat2(cos(sin(T*.2)+vec4(0,11,33,0))); p.z += T*.3; d = p.y; h = p.y+.5; for (a=.01; a<1.; a+=a) p.xz *= mat2(8,6,-6,8)*.1, d += abs(dot(sin(p*.3/a + T*.3), vec3(1)))*a, h += abs(dot(sin(p.xz*.6/a), vec2(1)))*a; d = max(max(d-3., -d), 0.); x = h < .001 ? 0. : 1.; d *= k; e = min(p.y+2., 1.); w = g * (1. - exp(x==0. ? -h*3e2 : -d/.4)); c += w * x*d*4.5; c.r += w * e*e * .1; c += (d > 0. ? .0 : h/2e2) * vec3(3,1,0); g -= w * x; t += min(h*.2, k); k /= .985; } O = vec4(tanh(c), 1); } Last pass of tricks: Merge and unroll more expressions together Use alternative forms for vec*(1) Rely on mathematical equivalences such as e^{-x}=1/e^x Some symbol names can be reused (see a) Notice how the rotation matrix coefficients (0,11,33,0) are close to the red factors (3,1,0)? That's right, we can factor that out into a shared constant K. Iterate i within the condition We're going to inline k*=1.015 inside the min(): this is not equivalent, but in practice it makes no difference The first 5 instructions of the main loop go into the initialization placeholder of the inner for, and all the others go into the iteration placeholder of the outer for, so that we can remove all {} Declare a z to be used instead of 0. since we have a bunch of them (thanks coyote) The x = h < .001 ? 0. : 1. can also be obtained progressively through some increment trick (thanks coyote) I'm also reordering a bit some instructions for clarity 🙃 // 448 chars void main() { vec3 c,p,K=vec3(3,1,0); for(float z,i,a,g=1.,t,h,d,w,k=.15; i++<1e2; d = max(max(d-3.,-d),a=z)*k, w = g-g/exp(h>.001?a++,d/.4:h*3e2), g -= a*=w, c += a*d*4.5+(d>z?z:h/2e2)*K, a = min(p.y+2.,1.), c.r += w*a*a*.1, t += min(h*.2,k/=.985)) for(p=normalize(vec3(P+P-R,R.y))*t, p.xz*=mat2(cos(sin(T*.2)+K.zyxz*11.)), p.z+=T*.3, d=p.y,h=d+.5,a=.01;a<1.;a+=a) p.xz *= mat2(8,6,-6,8)*.1, d += abs(dot(sin((p/a+T)*.3),p-p+a)), h += abs(dot(sin(p.xz*.6/a),P-P+a)); O = vec4(tanh(c),1); } And here we are. All we have to do now is remove all unnecessary spaces and line breaks to obtain the final version. I'll leave you here with this readable version. Golfer by Courtney Cook (Unsplash) Forewords I'm definitely breaking the magic of that artwork by explaining everything in detail here. But it should be replaced with an appreciation for how much concepts, math, and art can be packed in so little space. Maybe this is possible because they fundamentally overlap? Nevertheless, writing such a piece was extremely refreshing and liberating. As a developer, we're so used to navigate through mountains of abstractions, dealing with interoperability issues, and pissing glue code like robots. Here, even though GLSL is a very crude language, I can't stop but being in awe by how much beauty we can produce with a standalone shader. It's just... Pure code and math, and I just love it.

The idea that we must always linearize sRGB gradients or work in a perceptually uniform colorspace is starting to be accepted universally. But is it that simple? When I learned about the subject, it felt like being handed a hammer and using it everywhere. The reality is a bit more nuanced. In this article we will see when to use which, how to use them, and we will then see why the situation is more dire than it looks. Code snippets Before we start, since we are going to use GLSL as language, following are the reference functions we will use for the rest of the article. vec3 s2l(vec3 c) { // sRGB to linear return mix(c/12.92, pow((max(c,0.0)+0.055)/1.055,vec3(2.4)), step(vec3(0.04045),c)); } vec3 l2s(vec3 c) { // linear to sRGB return mix(c*12.92, 1.055*pow(max(c,0.0),vec3(1./2.4))-0.055, step(vec3(0.0031308),c)); } vec3 l2oklab(vec3 rgb) { // linear to OkLab const mat3 rgb2lms = mat3( +0.4122214708, +0.2119034982, +0.0883024619, +0.5363325363, +0.6806995451, +0.2817188376, +0.0514459929, +0.1073969566, +0.6299787005); const mat3 lms2lab = mat3( +0.2104542553, +1.9779984951, +0.0259040371, +0.7936177850, -2.4285922050, +0.7827717662, -0.0040720468, +0.4505937099, -0.8086757660); vec3 lms = rgb2lms * rgb; return lms2lab * pow(lms, vec3(1.0/3.0)); } vec3 oklab2l(vec3 lab) { // OkLab to linear const mat3 lab2lms = mat3( +1.0000000000, +1.0000000000, +1.0000000000, +0.3963377774, -0.1055613458, -0.0894841775, +0.2158037573, -0.0638541728, -1.2914855480); const mat3 lms2rgb = mat3( +4.0767416621, -1.2684380046, -0.0041960863, -3.3077115913, +2.6097574011, -0.7034186147, +0.2309699292, -0.3413193965, +1.7076147010); vec3 lms = lab2lms * lab; return lms2rgb * (lms*lms*lms); } Also, the output of the pipeline will be expected to be sRGB all the time. Color gradients To illustrate how sRGB, linear RGB and OkLab respectively look like, let's interpolate between two colors with each one of them: Color gradients from top to bottom: sRGB, linear, OkLab The 3 stripes were generated like this: vec3 o_srgb = mix(c0, c1, v); vec3 o_linear = l2s(mix(s2l(c0), s2l(c1), v)); vec3 o_oklab = l2s(oklab2l(mix(l2oklab(s2l(c0)), l2oklab(s2l(c1)), v))); Where v is simply the x coordinate between 0 and 1, c0 the left color, and c1 the right one. Note The input colors are considered to be sRGB in input. Similarly, we always make sure to output sRGB at the end (with l2s()) because that's what the pipeline expects. Key takeaways: sRGB is not acceptable because of this grayish/brownish zone, which is also perceived darker. In various situations this creates undesirable muddy midtones. In general, it's wrong and broken to do that. Linear is better from a purely physical point of view as it models the mixing of light energy properly. But from a color perception point of view it's not ideal, for example here it has this transition into pinkish which might not be desirable. The last one is using OkLab for a perceptually uniform gradient; it is the one providing the best result for our human perception, at a certain performance cost. The general consensus is as follow: if you need a color transition within a shape or texture, or some sort of color map, OkLab is the best tool, while linear is cheap, physically correct, and usually acceptable visually. But what about monochrome gradients? Things are not as obvious as they seem when we work in monochrome. If instead of red and blue we pick black and white, this is what happens: Grayscale gradients from top to bottom: sRGB, linear, OkLab Suddenly this tells a whole different story. sRGB becomes perfectly acceptable, while linear favors way too much the lightness, and OkLab remains the best. The linear gradient felt acceptable before, but now it is highly questionable. Just to be clear, the linear strip is linear, you can see it as linear energy or casually said "wattage", to which our perception does respond non-linearly. At this point one may even argue that sRGB looks best. So what can we do about this? First of all, we always need to question what we are trying to achieve, and fortunately sometimes we can take a few shortcuts. For example, let's say we want to depict a heat map in black and white. In my previous article I had to display 2D noise, so I wanted the observer to experience a linear perception of the "height" of the noise. In this case, working in sRGB (that is, doing zero effort with regards to perception) is actually a better call than mixing between black and white in linear space: Noise 2D with height as sRGB (left) or linear (right) Here we are comparing these two: vec3 o_srgb = vec3(v); // equivalent to mix(black, white, v) vec3 o_linear = l2s(vec3(v)); // equivalent to l2s(mix(black, white, v)) Note We removed the mix out of the formulas because black=vec3(0) and white=vec3(1), which have the same value when uncompressed to linear space. To do things right we may want to use OkLab but this feels overkill since this is just a straightforward monochromatic signal. Fortunately, the perceptual lightness can be fairly simple to model. With monochromatic input, OkLab uses L=x³, which is basically equivalent to do a gamma correction with γ=3. This means that we can simplify the OkLab interpolation we used before to the very simple: vec3 o_oklab = l2s(vec3(v*v*v)); // equivalent to l2s(oklab2l(vec3(v,0,0))) Noise 2D with height remapped to human lightness perception Doing this simple operation is exactly equivalent to interpolating between black and white in OkLab space, except it's just 2 extra multiplications. We still need to be extra careful if we want to swap the black and white. v needs to be swapped before the gamma encoding, and that means before the sRGB gamma encoding as well: Top to bottom: srgb(1-v³) (incorrect), 1-srgb(v³) (incorrect), srgb((1-v)³) One extra trick: combining gammas sRGB has a curve that closely approximates a gamma correction γ=2.2. So sometimes, instead of using l2s(rgb), we may prefer to use the simpler pow(rgb,vec3(1.0/2.2)). It means we could replace l2s(vec3(v*v*v)) with the following to merge the two operations: vec3 o_oklab = vec3(pow(v, 3.0/2.2)); // combination of v³ and gamma 2.2 (sRGB-like) encoding And the white-to-black version: vec3 o_oklab = vec3(pow(1.0-v, 3.0/2.2)); Warning Whenever you use pow, make sure your input is positive. Adding a max(v,0.0) for safety might be reasonable in certain cases. The difference between a proper sRGB conversion and the combined gamma is pretty small: Top: srgb(v³), bottom: v^(3/2.2) Alpha blending and pre-multiplication Sometimes, instead of fading colors into each others, we need to compose shapes, textures, masks, ... This need for compositing, or blending, arises when the pipelines are separated, meaning we are not working in the same fragment shader for everything. For example, we could have a shape generated in a fragment, which we need to overlay onto a surface. That shape might have some non-binary transparency, either for anti-aliasing purposes, a blur, or similar. An example of a shape partially transparent If that shape were to be blend onto another colored surface, we would like to have the same effect as the gradient earlier. For well-known reasons, it is likely that this shape would end up as a pre-multiplied color, which would be blend onto one or more layers. If what I just said is confusing, I recommend checking out this good article on alpha compositing from Bartosz Ciechanowski. The literature is quite extensive on the subject so I will assume familiarity with it. Of course, if we are to do things right, the blending would have to happen in linear space. Do not consider sRGB blending in alpha blending, it's an even more terrible idea than before because of the bilinear filtering, transforms or mipmaps that can happen between the pre-multiplication and the blending itself. But that means we would end up with the linear gradient shortcomings from earlier, wouldn't we? And this is where things get ugly. Look at the difference between a linear and an OkLab blending, in black and white: Blending of a blurry white circle onto black, left is linear, right is OkLab If we invert the colors: Blending of a blurry black circle onto white, left is linear, right is OkLab We have the exact same problem as earlier, but seeing it with an actual blending of shapes makes the problem particularly striking. The white and black OkLab circles look the same size (because they are), and they don't have the unfortunate "bobbing" effect of the linear version (on the white onto black). Warning The OkLab blending is done with pre-multiplied Lab colors. It is important not to pre-multiply linear values which are then converted to OkLab, this will give very unexpected results. The problem is, it is very unlikely that your whole graphics pipeline would switch to OkLab for every textures and buffers. And since most of the time the pipelines are built for more than just black and white, the cube hack suggested earlier has a very limited scope. In the case of shape blending, it is almost certain that the whole pipeline would not be contained in a single shader where you can just mix in OkLab. You're probably thinking of using sRGB, but in a blending pipeline this really is a terrible idea. Final words In practice, neither sRGB nor pure linear blending give good results, and using OkLab is not always an option. And unfortunately, I don't have a good answer to this whole situation. My next article is about anti-aliasing where this problem also exists, and I must admit this whole ordeal puts me in quite some distress; I had to talk about this issue first.