Distance field composition and blending
How to blend or compose implicit surfaces - 09/2023 - #Jumble
To begin
This state of the art provides a nice start for simple operators:
https://inria.hal.science/hal-00753246/document
Even better Florian Canezin thesis:
https://theses.hal.science/tel-01585541/document
"Study of the composition models of field functions in computer graphics"
Inigo Quilez impressive list of operators (global support):
https://iquilezles.org/articles/distfunctions/
Ronja:
https://www.ronja-tutorials.com/post/035-2d-sdf-combination/
Other
French slides of my state of the art (soa) about implicit surface (SDF) blending operators.
sdf_blending_soa_master_presentation_fr.pdf
Note: see what's the difference with compact and global support on my blog post here.
Global support
For primitives like the sphere equation for instance $f(x,y,z) = \|c -x\|^2 $
R-function (Rvachech)
$C^1$ continuity, and $C^0$ when $f_1 = f_2$
$$G_\cup(f_1, f_2) = f_1 + f_2 + \sqrt{f_1^2 + f_2^2}$$ $$G_\cap(f_1, f_2) = -G_\cup(-f_1, -f_2)$$ $$G_/(f_1, f_2) = -G_\cap(f_1, -f_2)$$$C^m$ continuity
$$ G_\cup(f_1, f_2) = f_1 + f_2 + \sqrt{f_1^2 + f_2^2} (f_1^2 + f_2^2)^{m/2} $$Rounded blend
Pasko Blend
$$G_\cup(f_1, f_2) = f_1 + f_2 + g(f_1, f_2)$$ $$ g(f_1, f_2) = \frac{a_0}{1 + \frac{f_a}{a_1}^2 + \frac{f_2}{a_2}^2 } $$Tweaking the parameters $a_0$, $a_1$, $a2$ defines how sharp or rounded the blend will be.
Smooth min
Popular smin variant
\[ \text{smoothMin}(a, b, k) = -\frac{\log\left( e^{-ka} + e^{-kb} \right)}{k} \]
float smoothMax(float a, float b, float k) {
return log(exp(k*a) + exp(k*b)) / k;
}
float smoothMin(float a, float b, float k){
return -smoothMax(-a, -b, k);
}
Gradient Vector
\[ \nabla \text{smoothMin}(a, b, k) = \left( \frac{e^{-ka}}{e^{-ka} + e^{-kb}}, \quad \frac{e^{-kb}}{e^{-ka} + e^{-kb}} \right) \]
// Analatycal gradient
vec2 smoothMinGradient(float a, float b, float k) {
float expKa = exp(-k * a);
float expKb = exp(-k * b);
float denom = expKa + expKb;
float da = expKa / denom;
float db = expKb / denom;
return vec2(da, db);
}
// Central finite difference approximation of gradient
vec2 smoothMinGradientFiniteDiff(float a, float b, float k, float h) {
float dA = (smoothMin(a + h, b, k) - smoothMin(a - h, b, k)) / (2.0 * h);
float dB = (smoothMin(a, b + h, k) - smoothMin(a, b - h, k)) / (2.0 * h);
return vec2(dA, dB);
}
Compact support
Will use squared euclidean distance, for instance $d = x^2 + y^2 + z^2$. Typically applied compactly supported primitives such as the blob function $f(x,y,z) = b e^{-a \ d^2(x,y,z)}$
B.Wyvill
$f(d) = \frac{-4}{9R^6} (d^2 - R^2)^2 (d^2 - (\frac{3}{2}R) \quad \text{ if } \ d < R$
$f(d) = 0 \quad \text{ otherwise}$
N.Stolte
$f(d) = \frac{1}{R^8}(d^2-R^2)^4 \quad \text{ if } d < R$
$f(d) = 0 \quad \text{ otherwise}$
Ricci
$$G_\cup(f_1, f_2) = (f_1^n + f_2^n)^{1/n}$$ $$G_\cap(f_1, f_2) = ((-f_1)^n + (-f_2)^n)^{1/n}$$ $$G_/(f_1, f_2) = G_\cap(f_1, 1-f_2)$$As $n \rightarrow \infty$ the blending gets sharper
Note: Ricci min/max works as well (<- I don't remember what I meant by that...)
Next
Things to improve in the future:
- provide shader toy live examples.
- Give the research paper sources.
Visualize blending operators with colors
/*
The goal is to show how one can vizualize blending operators used to combine SDFs
This makes it easier to interpret the behavior of the 2D operator.
For instance we can see that smoothMin() is very similar to the min() function,
except it produces smoothed isolines around the diagonal y = x
which is where the two sdf values x and y are equal, that is to say where the shapes of the
sdf intersect.
When getting used to this way of graphing and interpreting the operators
and how it relates with the input SDFs values then you can start imagining
drawing operators by hand and find ways to express it mathematically in order to get
the desired blending between two SDFs.
See fig 16 of https://inria.hal.science/hal-00753246/document to see an operator
that produce a bulge when two sdf "collide"
(material seems to be pushed around the line of intersection).
This was crafted in part thanks to the graphical insight
on how the "shape" of the blending operator effects the final result.
*/
float level_curves_greyscale(float f, float scale)
{
float opacity = cos(f * 3.141592 * scale + 0.5);
opacity *= opacity;
opacity = 1.0 - opacity;
opacity *= opacity;
opacity *= opacity;
opacity = 1.0 - opacity;
return opacity;
}
//------------------------------------------------------------------------
vec4 level_curves_color(float f, float iso_level, float scale)
{
float opacity = level_curves_greyscale(f, scale);
vec4 color;
if(f > iso_level) color = vec4(0.9, 0.2, 0.1, 1.0); // Inside red
else if(f <= iso_level) color = vec4(0.1, 0.2, 0.9, 1.0); // Outside blue
if(true)
{
if(f >= 1.0 ) color = vec4(1.0, 0.75, 0.5, 1.0);
else if(f < -1.0 ) color = vec4(0.9, 0.90, 0.2, 1.0);
else if( isnan(f) ) color = vec4(1.0, 1.0, 0.0, 1.0); // yellow
// Color strip that belong to iso-surface
if( abs(f-iso_level) < (1./scale) * 0.5 )
color = vec4(0.2, 0.9, 0.1, 1.0); // greenish
//color = Vec3(0.5); // grey
}
/*
float eps = 1e-4f;
if( f < (1.f + eps) && f > (1.f - eps))
return vec4(0.0,0.0,0.0,1.0); // black
if( f < (eps*10.0) && f > 0.0 )
return vec4(0.0,0.0,0.0,1.0); // black
*/
// glsl equivalent of lerp:
color = mix(color, vec4(1.0), opacity);
color.a = 1.0 - opacity;
return color;
}
//------------------------------------------------------------------------
float smoothMax(float a, float b, float k) {
return log(exp(k*a) + exp(k*b)) / k;
}
float smoothMin(float a, float b, float k){
return -smoothMax(-a, -b, k);
}
//------------------------------------------------------------------------
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
// Normalized pixel coordinates (from 0 to 1 on the y axis)
vec2 p = fragCoord / iResolution.y;
// place 0 at iResolution.y/2 and scale 10 times
// so that y ranges from [-10(bottom); 10(top)]
p -= 0.5;
p *= 10.0;
// Switch every 3 seconds
int mode = int(iTime / 3.0) % 3;
// Output value of our 2D function
float value = 0.0;
if( mode == 1 )
// Intersection of SDFs
value = max(p.x, p.y);
else if( mode == 2 )
// Union of SDFs
value = min(p.x, p.y);
else
// Smooth union of SDFs
value = smoothMin(p.x, p.y, 1.0);
// Output color to screen
fragColor = level_curves_color(value,
0.0/*central isoline (green)*/,
6.5/*spacing between isos*/ );
}
/*
The goal is to show how one can vizualize blending operators float g( float f1, float f2) used to combine SDFs
This makes it easier to interpret the behavior of the 2D operator.
This is an extension to https://www.shadertoy.com/view/DtlyWB
with visualization of the 2d gradient g_x (red) and g_y (blue)
*/
#define M_PI 3.141592653589793238462643383279502
#define SHOW_GRAD_COLOR 1 // set to zero to only see isolines.
#if 1
float raise_power(float a)
{
float g = a;
g*= g;
g*= g;
g*= g;
g*= g;
g*= g;
return g;
}
// -----------------------------------------------------------------------------
vec2 grad_color(vec2 gf)
{
float s = gf.x + gf.y;
s = 1.0 - s;
if( s < 0.0 )
return gf;
else
return gf + s * 0.5;
}
// -----------------------------------------------------------------------------
vec4 color_bin_op_core(
float ddx,
float ddy,
float dtan,
vec2 gdres,
float ddres
)
{
const float line_min = 0.03;
float r, g, b;
if(false)
gdres = grad_color(gdres);
float spacing = 7.0 /* * 2.0 * 50.*/; //< spacing between iso lines
g = cos( ddres * M_PI * spacing);
g = raise_power( g );
g*= g;
#if SHOW_GRAD_COLOR
r = gdres.x;
b = gdres.y;
if(abs(ddres-0.5f)<line_min)
{
r = max(g, r);
b = max(g, b);
float u = g;
g = cos( (abs(dtan) < line_min) ? abs(dtan) * 50.0 : 1.0 );
g = raise_power(g);
r -= g;
b -= g;
g = max(g,u);
}
else
{
r *= (1.0 - g);
b *= (1.0 - g);
g = (abs(dtan) < 0.003f) ? 1.0 : 0.0;
g = cos( (abs(dtan) < line_min) ? abs(dtan) * 50.0 : 1.0 );
g = raise_power(g);
r -= g;
b -= g;
}
// Clamp between [0 1]
r = max(0.0, min(1.0, r));
g = max(0.0, min(1.0, g));
b = max(0.0, min(1.0, b));
if( ddres < 0.0 ){ // if g() is negative
// r = b = 0.0; // set to black
g = 0.3; // greenish tint
}
return vec4(r, g, b, 1.0);
#else
float a = g;
r = (1.0 - a);
g = (1.0 - a);
b = (1.0 - a);
a = (abs(dtan) < 0.003f) ? 1.0 : 0.0;
a = cos( (abs(dtan) < line_min) ? abs(dtan) * 50.0 : 1.0 );
a = raise_power(a);
r -= a;
g -= a;
b -= a;
if (abs(dtan) > 0.003){
if (ddres > 1.0+line_min) {// higher than 1 => green(0, 1, 0)
g = 1.0-a;
} else if (ddres > 1.0-line_min){// equal 1 => black(0, 0, 0)
;
} else if (ddres > 0.5+line_min) {// inside => red(0.9f, 0.2f, 0.1f)
r = 1.0-a;
} else if (ddres < 0.5-line_min) {// outside => blue(0.1f, 0.2f, 0.9f)
b = 1.0-a;
} else { // surface => black (0, 0, 0)
r = 1.0-a;
b = 1.0-a;
}
} else {
g = a;
r = a;
}
// Clamp between [0 1]
r = max(0.0, min(1.0, r));
g = max(0.0, min(1.0, g));
b = max(0.0, min(1.0, b));
return vec4(r, g, b, 1.0);
#endif
}
void color_bin_op(float ddx, //< coords where the op was evaluated
float ddy, //< coords where the op was evaluated
vec2 gdres, //< operator's gradient
float ddres, //< operator's value
bool gradient_based_operator)
{
float angle = 20.0; /* user value */
bool opening = false; /* user value */
float tant = angle;
float dtan = 100.0;
if( gradient_based_operator && opening){
// todo specific to gradient based blending operators (https://inria.hal.science/hal-00753246/document)
dtan = 100.0; // eval_opening(gb_operator->get_opening_kind(), ddx, ddy, tant);
}
color_bin_op_core(ddx, ddy, dtan, gdres, ddres);
}
#endif
float level_curves_greyscale(float f, float scale)
{
float opacity = cos(f * 3.141592 * scale + 0.5);
opacity *= opacity;
opacity = 1.0 - opacity;
opacity *= opacity;
opacity *= opacity;
opacity = 1.0 - opacity;
return opacity;
}
//------------------------------------------------------------------------
vec4 level_curves_color(float f, float iso_level, float scale)
{
float opacity = level_curves_greyscale(f, scale);
vec4 color;
if(f > iso_level) color = vec4(0.9, 0.2, 0.1, 1.0); // Inside red
else if(f <= iso_level) color = vec4(0.1, 0.2, 0.9, 1.0); // Outside blue
if(true)
{
if(f >= 1.0 ) color = vec4(1.0, 0.75, 0.5, 1.0);
else if(f < -1.0 ) color = vec4(0.9, 0.90, 0.2, 1.0);
else if( isnan(f) ) color = vec4(1.0, 1.0, 0.0, 1.0); // yellow
// Color strip that belong to iso-surface
if( abs(f-iso_level) < (1./scale) * 0.5 )
color = vec4(0.2, 0.9, 0.1, 1.0); // greenish
//color = Vec3(0.5); // grey
}
/*
float eps = 1e-4f;
if( f < (1.f + eps) && f > (1.f - eps))
return vec4(0.0,0.0,0.0,1.0); // black
if( f < (eps*10.0) && f > 0.0 )
return vec4(0.0,0.0,0.0,1.0); // black
*/
// glsl equivalent of lerp:
color = mix(color, vec4(1.0), opacity);
color.a = 1.0 - opacity;
return color;
}
//------------------------------------------------------------------------
float smoothMax(float a, float b, float k) {
return log(exp(k*a) + exp(k*b)) / k;
}
float smoothMin(float a, float b, float k){
return -smoothMax(-a, -b, k);
}
vec2 smoothMinGradient(float a, float b, float k) {
float expKa = exp(-k * a);
float expKb = exp(-k * b);
float denom = expKa + expKb;
float da = expKa / denom;
float db = expKb / denom;
return vec2(da, db);
}
// Central finite difference approximation of gradient
vec2 smoothMinGradientFiniteDiff(float a, float b, float k, float h) {
float dA = (smoothMin(a + h, b, k) - smoothMin(a - h, b, k)) / (2.0 * h);
float dB = (smoothMin(a, b + h, k) - smoothMin(a, b - h, k)) / (2.0 * h);
return vec2(dA, dB);
}
//------------------------------------------------------------------------
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
// Normalized pixel coordinates (from 0 to 1 on the y axis)
vec2 p = fragCoord / iResolution.y;
// place 0 at iResolution.y/2 and scale 10 times
// so that y ranges from [-10(bottom); 10(top)]
p -= 0.5;
p *= 10.0;
//ddx = ddx * 0.1 + 0.45;/////////////////DEBUG zoom in
//ddy = ddy * 0.1 + 0.45;/////////////////DEBUG zoom in
// Switch every 3 seconds
int mode = int(iTime / 3.0) % 3;
// mode = 3; // DEBUG
// Output value of our 2D function:
float value = 0.0;
//partial derivative in x and y of our 2D function:
vec2 op_grad = vec2(0.0, 0.0);
if( mode == 1 ){
// Intersection of SDFs
value = max(p.x, p.y);
op_grad = p.x > p.y ? vec2(1.0, 0.0) : vec2(0.0, 1.0);
}
else if( mode == 2 ) {
// Union of SDFs
value = min(p.x, p.y);
op_grad = p.x > p.y ? vec2(1.0, 0.0) : vec2(0.0, 1.0);
}
else {
// Smooth union of SDFs
value = smoothMin(p.x, p.y, 1.0);
//op_grad = smoothMinGradient(p.x, p.y, 1.0);
op_grad = smoothMinGradientFiniteDiff(p.x, p.y, 1.0, 0.0001);
}
// Output color to screen
fragColor = level_curves_color(value,
0.0/*central isoline (green)*/,
6.5/*spacing between isos*/ );
fragColor = color_bin_op_core(p.x, p.y, /*dummy value for dtan:*/100.0, op_grad, value);
}
No comments