Intro
One of the basic building blocks of implicit procedural modeling (such as when building a distance field for raymarching based on basic primitives) is the union operator.
float opU( float d1, float d2 )
{
return min( d1, d2 );
}
This operator works great, but has the problem that the resulting shape has discontinuities in its derivatives. Or in other words, the resulting surface of unifying two smooth objects is not a smooth surface anymore. This is often inconvenient from a looks perspective, such as when trying to model organic shapes.
smooth-min() based primitive union
Regular min() based primitive union
Several implementations
The way to smoothly blend the shapes is to get rid of the discontinuity of the min() function, of course. But we want our smooth-min function to behave quite like min() when one of the two primitives is way further that the other. It's only in the area where the two values get similar that we want to apply the smoothness.
// exponential smooth min (k=32)
float smin( float a, float b, float k )
{
float res = exp2( -k*a ) + exp2( -k*b );
return -log2( res )/k;
}
// power smooth min (k=8)
float smin( float a, float b, float k )
{
a = pow( a, k ); b = pow( b, k );
return pow( (a*b)/(a+b), 1.0/k );
}
// root smooth min (k=0.01)
float smin( float a, float b, float k )
{
float h = a-b;
return 0.5*( (a+b) - sqrt(h*h+k) );
}
// polynomial smooth min 1 (k=0.1)
float smin( float a, float b, float k )
{
float h = clamp( 0.5+0.5*(b-a)/k, 0.0, 1.0 );
return mix( b, a, h ) - k*h*(1.0-h);
}
// polynomial smooth min 2 (k=0.1)
float smin( float a, float b, float k )
{
float h = max( k-abs(a-b), 0.0 )/k;
return min( a, b ) - h*h*k*(1.0/4.0);
}
It worth noting that the exponential and power based smooth-min functions, unlike the polynomial or the root smooth-minimum, both generalize to more than two distances, so they are probably better suited for computing minimun distances to big sets of points beyond 2, for example when you want to compute smooth voronoi patterns or interpolate pointclouds. In the case of the power based smooth-min function, the expression a*b/(a+b) generalizes with the same formula as when computing the global resistance of N parallel resistors: 1/ ( 1/a + 1/b + 1/c + ... ). For example, for three distances, you get a*b*c / (a*b + b*c + c*a).
Besides accepting an arbitrary number of points, another advantage of the exponential smin() over the polynomial or root smin() is that when called multiple times with two arguments at a time, the exponential smin() produces the same result regardless of the order of the operations. The polynomial and root smin() functions however are not order independent. To make it more explicit, smin(a,smin(b,c)) is equal to smin(b,smin(a,c)) for the exponential smin(), but not for the polynomial or root. That means that the exponential smin() allows one to process long lists of distances in any arbitrary order and slowly compute the smin, while the polynomial is ordering dependent. This can be useful to compute smooth voronoi patterns for example.
About 5 years later after writing about the first polynomial smooth minimum here, I learnt that Media Molecule used the same polynomial smin() for their game "Dreams" (credited to Dave Smith),although their rewrote it in the form showed at the top and labeled as "smooth min 2" above. This form is mathematically equivalent but more efficient from a computational standpoint, and is the form I use these days the most:
// polynomial smooth min
float smin( float a, float b, float k )
{
float h = max( k-abs(a-b), 0.0 )/k;
return min( a, b ) - h*h*k*(1.0/4.0);
}
As noted by Shadertoy user TinyTexel, this can be generalized to higher levels of continuity than the quadratic polynomial offers (C1), which might be important for preventing lighting artifacts. Moving on to a cubic curve gives us C2 continuity, and doesn't get a lot more expensive than the quadratic one anyways:
// polynomial smooth min
float sminCubic( float a, float b, float k )
{
float h = max( k-abs(a-b), 0.0 )/k;
return min( a, b ) - h*h*h*k*(1.0/6.0);
}
Lastly, it's worth mentioning that since the polynomial smin() is always smaller or equal to regular min() by design, it is well suited for raymarching algorithms since the ray will never overshoot past the original sharp intersection.
Two functions, f(x)=e-x and g(x)=sin(8x)
Their min() in blue min(f(x),g(x)), and polynomial smin() in purple smin(f(x),g(x))
Derivation
Deriving the polynomial smin() function is not difficult. The easiest is probably to start with the simplified version:
smin( f(x), g(x), k ) = min( f(x), g(x) ) - w(x,k)
where w(x,k) is the correction we do to the regular min(), that is, the little polynomial arc between the points at which f(x) and g(x) differ by k. That is, if x=a is the point at which f(x)-g(x) = -k, and x=b is the point where f(x)-g(x) = k, and a middle point x=c is where f(x)=g(x), then we know that
w(a)=0, w(b)=0, w(c)=s
with s being the maximum value. Since we want w(x) to be a smooth function that connects nicely to the curves f(x) and g(x) at the points x=a and x=b, we can choose
w(x) = s·hn(x)
with h(x) = 1 ± [f(x)-g(x)]/k
meaning h(x) is a power curve of degree n, with range between 0 and 1, and with the signs being positive if x < c and negative if x > c. This gives two versions of w(x), which we can call wl(x) and wr(x) for "left" and "right".
Since we want continuity also at x=c, we need to make sure that the derivatives of smin match when coming from either left or right of x=c. Therefore,
f'(c) + wl'(c) = g'(c) + wr'(c)
which means that
f'(c) + [f'(c)-g'(c)]nS/k = g'(c) - [f'(c)-g'(c)]·n·s/k
This can only be solved if 1+2·n·s/k = 0, which gives
s = -k/2n
which is what we used in the quadratic and cubic implementations above.
The quadratic smin() though doesn't have second derivatives for w(x), meaning that the only way the condition
f''(c) + wl''(c) = g''(c) + wr''(c)
can be met is by making w(x) a cubic. Fortunately, since h(c)=1, the same condition 1+2·n·s/k = 0 needs to be met for ensuring continuity of the second (or any higher order) derivative.
This method can be used to generate smooth minimums of higher degrees easily.
Mix factor
Besides smoothly blending values, it might be useful to compute also a blending factor that can be used for shading. For example, if the smooth-minimum is being used to blend SDF shapes, having a blend factor could be useful to blend the material properties of the two shapes during the transition area. For example, the image below shows the mix of a red and a blue materials based on this mix factor as computed by the code below, which returns the smooth-minimum in .x and the blend factor in .y:
vec2 smin( float a, float b, float k )
{
float h = max( k-abs(a-b), 0.0 )/k;
float m = h*h*0.5;
float s = m*k*(1.0/2.0);
return (a<b) ? vec2(a-s,m) : vec2(b-s,1.0-m);
}
vec2 sminCubic( float a, float b, float k )
{
float h = max( k-abs(a-b), 0.0 )/k;
float m = h*h*h*0.5;
float s = m*k*(1.0/3.0);
return (a<b) ? vec2(a-s,m) : vec2(b-s,1.0-m);
}
The generalization to any power n being:
vec2 sminN( float a, float b, float k, float n )
{
float h = max( k-abs(a-b), 0.0 )/k;
float m = pow(h, n)*0.5;
float s = m*k/n;
return (a<b) ? vec2(a-s,m) : vec2(b-s,1.0-m);
}
Mixing between materials together with the geometry
Analytical Gradient
Oftentimes we want to compute the gradient of an SDF without resorting to numerical methods. There are multiple reasons to do so, including placement and alignment of objects in 3D space. There's also theoretical interest in seeing what form it takes and therefore how the smooth minimum operator will behave within a raymarcher. We'll compute it for the quadratic polynomial smooths step, but computing it for the others is also easy. The result is that for smin(f,g,k), we have that
∇smin(f,g,k) = ∇f⋅h + (1-h)⋅∇g for f>g
∇smin(f,g,k) = ∇f⋅(1-h) + h⋅∇g for f<g
or in other words, the gradient of the smooth minimum is the linear blend or mix of the gradients of the input SDFs, ∇f and ∇g. These of course have length one, since they are SDFs:
∥∇f∥ = 1
∥∇g∥ = 1
However, as we know from lerping vs slerping vectors, even though both ∇f and ∇g have length of one the linear interpolation of two unit length vectors won't have unit length, but it will be shorter. We'll be cutting through the unit sphere on our way from one vector to the other, rathr than going through its surface. We can checkthis numerically too; indeed, if we do the math we get
∥∇f∥2 = 1 - 2h(1-h)(1-<∇f,∇g>)
which is clearly less than one. This means the smooth minimum of two SDFs is not an SDF, and that it's not the optimal way to raymarch things. That's bad, but not disastrous. After all, the gradient is always less than 1, which means that at least out raymarcher won't step over the geometry during our raymarching loop. You can see this clearly in the image below in the areas where the capsule and the rectangle meet and the field lines separate and no longer are equidistant.
Another way to think about this is that the smooth minimum is, by construction, always less than the minimum of a and b, which means that, again, it's a conservative computation that will prevent the raymarcher from tunneling through the geometry.
The code to compute the gradient is pretty simple, and can be done part of the smooth minimum function itself like this, assuming now a and b carry both the value of the SDF and its gradient together in a single vector:
// .x = f(p)
// .y = ∂s(p)/∂x
// .z = ∂s(p)/∂y
// .w = ∂s(p)/∂z
// .yzw = ∇s(p) with ∥∇s(p)∥<1 sadly
vec4 smin( in vec4 a, in vec4 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.yzw,b.yzw,(a.x<b.x)?n:1.0-n));
}
The gradient of smin() in 2D
You can see this code in action here: https://www.shadertoy.com/view/tdGBDt
Results
In general, the polynomial smooth-min function works very well, predictably and fast. It can be used to connect surfaces, such as snow and bridge in the images below.
Regular min()
Polynomial smooth-min()
Regular min()
Polynomial smooth-min().