::     ::

### Intro

We know that proper 2D and 3D SDFs have a gradient of length 1.0 everywhere, since SDFs measure distances, and the rate of change of distance with respect to distance will always be the identity. However, knowing the exact direction of the gradient is useful for computer graphics (for lighting, aligning objects to surfaces, etc etc). Of course, one can use central differences to compute gradients, but that has a high computational cost and can get tricky to tweak. Automatic differentiation works well with SDFs of course, but comes with a cost.

In this page you'll find a list of primitives and operators and their analytic gradients that I developed around 2020, where lots of terms are reused between the SDF computation and its gradient computation, making them very cheap. The list is incomplete, so feel free to carry on all primitives and operations that I haven't defined with automatic differentiation/duals. Also, this list just contains 2D SDFs, but extending them to 3D through revolution or extrusion is pretty easy.

Lastly, note that all the primitives listed here come with a link to a realtime online demo in Shadertoy. In fact, this public playlist contains all the shadertoy examples: https://www.shadertoy.com/playlist/M3dSRf, so don't miss that out.

### Primitives

All primitives are centered at the origin; you will have to transform the point to get arbitrarily rotated, translated and scaled objects, and transform the gradient back accordingly (remember not to use non uniform scales since such transforms do not preserve distances). All SDFs below return the actual distance in the .x component, and its partial derivatives with respect to x and y in the .y and .z components. In other words,

.x = f(p)
.y = ∂f(p)/∂x
.z = ∂f(p)/∂y
.yz = ∇f(p) with ‖∇f(p)‖ = 1 vec3 sdgCircle( in vec2 p, in float r ) { float d = length(p); return vec3( d-r, p/d ); } // sc = sin/cos of aperture vec3 sdgPie( in vec2 p, in vec2 sc, in float r ) { float s = sign(p.x); p.x = abs(p.x); float l = length(p); float n = l - r; vec2 q = p - sc*clamp(dot(p,sc),0.0,r); float m = length(q) * sign(sc.y*p.x-sc.x*p.y); vec3 res = (n>m) ? vec3(n,p/l) : vec3(m,q/m); return vec3(res.x,s*res.y,res.z); } // sc = sin/cos of aperture vec3 sdgArc( in vec2 p, in vec2 sc, in float ra, in float rb ) { vec2 q = p; float s = sign(p.x); p.x = abs(p.x); if( sc.y*p.x > sc.x*p.y ) { vec2 w = p - ra*sc; float d = length(w); return vec3( d-rb, vec2(s*w.x,w.y)/d ); } else { float l = length(q); float w = l - ra; return vec3( abs(w)-rb, sign(w)*q/l ); } } vec3 sdgSegment( in vec2 p, in vec2 a, in vec2 b, in float r ) { vec2 ba = b-a, pa = p-a; float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 ); vec2 q = pa-h*ba; float d = length(q); return vec3(d-r,q/d); } vec3 sdgVesica(vec2 p, float r, float d) { vec2 s = sign(p); p = abs(p); float b = sqrt(r*r-d*d); if( (p.y-b)*d>p.x*b ) { vec2 q = vec2(p.x,p.y-b); float l = length(q)*sign(d); return vec3( l, s*q/l ); } else { vec2 q = vec2(p.x+d,p.y); float l = length(q); return vec3( l-r, s*q/l ); } } vec3 sdgBox( in vec2 p, in vec2 b ) { vec2 w = abs(p)-b; vec2 s = vec2(p.x<0.0?-1:1,p.y<0.0?-1:1); float g = max(w.x,w.y); vec2 q = max(w,0.0); float l = length(q); return vec3( (g>0.0)?l :g, s*((g>0.0)?q/l:((w.x>w.y)?vec2(1,0):vec2(0,1)))); } vec3 sdgCross( in vec2 p, in vec2 b ) { vec2 s = sign(p); p = abs(p); vec2 q = ((p.y>p.x)?p.yx:p.xy) - b; float h = max( q.x, q.y ); vec2 o = max( (h<0.0)?vec2(b.y-b.x,0.0)-q:q, 0.0 ); float l = length(o); vec3 r = (h<0.0 && -q.x<l)?vec3(-q.x,1.0,0.0):vec3(l,o/l); return vec3( sign(h)*r.x, s*((p.y>p.x)?r.zy:r.yz) ); } vec3 sdgHexagon( in vec2 p, in float r ) { const vec3 k = vec3(-0.866025404,0.5,0.577350269); vec2 s = sign(p); p = abs(p); float w = dot(k.xy,p); p -= 2.0*min(w,0.0)*k.xy; p -= vec2(clamp(p.x, -k.z*r, k.z*r), r); float d = length(p)*sign(p.y); vec2 g = (w<0.0) ? mat2(-k.y,-k.x,-k.x,k.y)*p : p; return vec3( d, s*g/d ); } vec3 sdgTriangleIsosceles( in vec2 p, in vec2 q ) { float w = sign(p.x); p.x = abs(p.x); vec2 a = p - q*clamp( dot(p,q)/dot(q,q), 0.0, 1.0 ); vec2 b = p - q*vec2( clamp( p.x/q.x, 0.0, 1.0 ), 1.0 ); float k = sign(q.y); float l1 = dot(a,a); float l2 = dot(b,b); float d = sqrt((l1<l2)?l1:l2); vec2 g = (l1<l2)? a: b; float s = max( k*(p.x*q.y-p.y*q.x),k*(p.y-q.y) ); return vec3(d,vec2(w*g.x,g.y)/d)*sign(s); } float cro( in vec2 a, in vec2 b ) { return a.x*b.y - a.y*b.x; } vec3 sdgTriangle( in vec2 p, in vec2 v ) { float gs = cro(v-v,v-v); vec4 res; { vec2 e = v-v, w = p-v; vec2 q = w-e*clamp(dot(w,e)/dot(e,e),0.0,1.0); float d = dot(q,q), s = gs*cro(w,e); res = vec4(d,q,s); } { vec2 e = v-v, w = p-v; vec2 q = w-e*clamp(dot(w,e)/dot(e,e),0.0,1.0); float d = dot(q,q), s = gs*cro(w,e); res = vec4( (d<res.x) ? vec3(d,q) : res.xyz, (s>res.w) ? s : res.w ); } { vec2 e = v-v, w = p-v; vec2 q = w-e*clamp(dot(w,e)/dot(e,e),0.0,1.0); float d = dot(q,q), s = gs*cro(w,e); res = vec4( (d<res.x) ? vec3(d,q) : res.xyz, (s>res.w) ? s : res.w ); } float d = sqrt(res.x)*sign(res.w); return vec3(d,res.yz/d); } float cro( in vec2 a, in vec2 b ) { return a.x*b.y - a.y*b.x; } vec3 sdgQuad( in vec2 p, in vec2 v ) { float gs = cro(v-v,v-v); vec4 res; { vec2 e = v-v, w = p-v; vec2 q = w-e*clamp(dot(w,e)/dot(e,e),0.0,1.0); float d = dot(q,q), s = gs*cro(w,e); res = vec4(d,q,s); } { vec2 e = v-v, w = p-v; vec2 q = w-e*clamp(dot(w,e)/dot(e,e),0.0,1.0); float d = dot(q,q), s = gs*cro(w,e); res = vec4( (d<res.x) ? vec3(d,q) : res.xyz, (s>res.w) ? s : res.w ); } { vec2 e = v-v, w = p-v; vec2 q = w-e*clamp(dot(w,e)/dot(e,e),0.0,1.0); float d = dot(q,q), s = gs*cro(w,e); res = vec4( (d<res.x) ? vec3(d,q) : res.xyz, (s>res.w) ? s : res.w ); } { vec2 e = v-v, w = p-v; vec2 q = w-e*clamp(dot(w,e)/dot(e,e),0.0,1.0); float d = dot(q,q), s = gs*cro(w,e); res = vec4( (d<res.x) ? vec3(d,q) : res.xyz, (s>res.w) ? s : res.w ); } float d = sqrt(res.x)*sign(res.w); return vec3(d,res.yz/d); } vec3 sdgEllipse( vec2 p, in vec2 ab ) { vec2 sp = sign(p); p = abs( p ); bool s = dot(p/ab,p/ab)>1.0; float w = atan(p.y*ab.x, p.x*ab.y); if(!s) w=(ab.x*(p.x-ab.x)<ab.y*(p.y-ab.y))? 1.570796327 : 0.0; for( int i=0; i<4; i++ ) { vec2 cs = vec2(cos(w),sin(w)); vec2 u = ab*vec2( cs.x,cs.y); vec2 v = ab*vec2(-cs.y,cs.x); w = w + dot(p-u,v)/(dot(p-u,u)+dot(v,v)); } vec2 q = ab*vec2(cos(w),sin(w)); float d = length(p-q); return vec3( d, sp*(p-q)/d ) * (s?1.0:-1.0); }

### Operators/modifiers

As said, once the gradients of the basic primitives are hardcoded, automatic differentiation can take over for most of the remaining composition work. However, if you really want to make it all crystal clear and remove all opaque aspects of your code, you might want to implement analytic gradient computation for the composition operators too. Here go a few:

Rounding

As we know, inflating any shape into a rounded self is pretty easy with the single subtraction of a quantity r. If such r is a constant, the derivatives of the SDF won't change, so we don't really need to do anything in order to continue getting correct gradients.

vec3 sdgRound( in vec2 p, in float r ) { vec3 dis_gra = sdgShape(p); return vec3( dis_gra.x - r, dis_gra.yz ); }

These are a few examples:    Onion - annular shapes

Similarly, shapes can be made annular (like a ring or the layers of an onion) by taking their absolute value and then subtracting a constant from their field. This operation is also very easy to carry over for the gradients:

vec3 sdgOnion( in vec2 p, in float r ) { vec3 dis_gra = sdgShape(p); return vec3( abs(dis_gra.x) - r, sign(dis_gra.x)*dis_gra.yz ); }

These are a few examples:    Minimum and Smooth Minimum   (https://www.shadertoy.com/view/tdGBDt)

One way to combine shapes is to take the minimum of two SDFs. Naturally, the same conditional branching that picks the smallest of the SDFs should be used to select and propagate the gradient. However, when a smooth minimum is used to smoothly connect shapes, the derivative selection is more involved. Luckily, using polynomial smooth-minimums have simple derivatives and accept an analytical propagation. Of course, bear in mind that the smooth-minimum does not respect the constraints of having unit length everywhere, so its output is only a true SDF (or almost) far away from the surface of the resulting shape. Anyways, here's the code for the analytic derivatives/gradient of a quadratic smooth-minimum.

vec3 sdgMin( in vec3 a, in vec3 b ) { return (a.x<b.x)?a:b; }

vec3 sdgSmoothMin( in vec3 a, in vec3 b, in float k ) { float h = max(k-abs(a.x-b.x),0.0); float m = 0.25*h*h/k; float n = 0.50* h/k; return vec3( min(a.x, b.x) - m, mix(a.yz, b.yz, (a.x<b.x)?n:1.0-n) ); }

Here's on the left the minimum, and on the right the smooth-minimum:  