Curvature of a Distance field / Implicit surface

Draft notes

Curvature signed distance field (SDF), implicit surface, level set

Mean curvature with cheap evaluation:

https://www.shadertoy.com/view/Xts3WM

//Cheap curvature by nimitz (twitter: @stormoid)

/*
	Not sure if this method is new, but I haven't seen it around.
	Using a single extra tap to return analytic curvature of an SDF.

	My first implementation required 16 taps, then 12, then thanks 
	to some math help from a friend (austere) I got it down to 7 taps.
	And this is the	final optimization which requires 5 taps or a
	single tap if you're already computing normals.


	Edit (April 2016):

	Coming back to this I know realize that what I am returning is an
	approximation of the Laplacian of the SDF and the Laplacian being 
	the divergence of the gradient of the field means that any point
	on the surface being a sink will return negative "curvature" and
	vice-versa (since the gradient of a SDF should point to the centroid
	at any point in space).

	So now, I include the 7-tap version which is a more accurate Laplacian,
	while the 5-tap version while cheaper, it is less accurate.

	N.B. Mean curvature is computed here, not Gaussian curvature. (thanks to Fabrice)
*/

#define ITR 80
#define FAR 10.
#define time iTime

mat2 mm2(in float a){float c = cos(a), s = sin(a);return mat2(c,s,-s,c);}

float map(vec3 p)
{
    p.x += sin(p.y*4.+time+sin(p.z))*0.15;
    float d = length(p)-1.;
    float st = sin(time*0.42)*.5+0.5; 
    const float frq = 10.;
    d += sin(p.x*frq + time*.3 + sin(p.z*frq+time*.5+sin(p.y*frq+time*.7)))*0.075*st;
    
    return d;
}

float march(in vec3 ro, in vec3 rd)
{
	float precis = 0.001;
    float h=precis*2.0;
    float d = 0.;
    for( int i=0; i<ITR; i++ )
    {
        if( abs( h ) < precis || d > FAR ) break;
        d += h;
	    float res = map(ro+rd*d);
        h = res;
    }
	return d;
}

//5 taps total, returns both normal and curvature
vec3 norcurv(in vec3 p, out float curv)
{
    vec2 e = vec2(-1., 1.)*0.01;   
    float t1 = map(p + e.yxx), t2 = map(p + e.xxy);
    float t3 = map(p + e.xyx), t4 = map(p + e.yyy);

    curv = .25/e.y*(t1 + t2 + t3 + t4 - 4.0*map(p));
    return normalize(e.yxx*t1 + e.xxy*t2 + e.xyx*t3 + e.yyy*t4);
}

//Curvature only, 5 taps, with epsilon width as input
float curv(in vec3 p, in float w)
{
    vec2 e = vec2(-1., 1.)*w;   
    
    float t1 = map(p + e.yxx), t2 = map(p + e.xxy);
    float t3 = map(p + e.xyx), t4 = map(p + e.yyy);
    
    return .25/e.y*(t1 + t2 + t3 + t4 - 4.0*map(p));
}

//Curvature in 7-tap (more accurate)
float curv2(in vec3 p, in float w)
{
    vec3 e = vec3(w, 0, 0);
    
    float t1 = map(p + e.xyy), t2 = map(p - e.xyy);
    float t3 = map(p + e.yxy), t4 = map(p - e.yxy);
    float t5 = map(p + e.yyx), t6 = map(p - e.yyx);
    /*
FabriceNeyret2, 2019-02-07
@nimitz :  Some remarks:

- you should state that what you propose is the mean curvature, not the Gaussian curvature.

- there is a scaling problem ( curvature on a radius-1 sphere is very wrong ), because there is a unit problem: as a second finite derivative, the denominator should be e.x*e.x , not e.x .
And for a mysterious reason the factor should be 0.5, not 0.25 .
With .5/(e.x*e.x) , I got curvature of 1-radius sphere = 1.
*/
    // Well with fabrice correction now curvature looks like it goes beyond the range [-1 1]
    // Is it normal?
    return 0.5/(e.x*e.x)*(t1 + t2 + t3 + t4 + t5 + t6 - 6.0*map(p)); //return .25/e.x*(t1 + t2 + t3 + t4 + t5 + t6 - 6.0*map(p)); } void mainImage( out vec4 fragColor, in vec2 fragCoord ) { vec2 p = fragCoord.xy/iResolution.xy-0.5; p.x*=iResolution.x/iResolution.y; vec2 mo = iMouse.xy / iResolution.xy-.5; mo.x *= iResolution.x/iResolution.y; vec3 ro = vec3(0.,0.,4.); ro.xz *= mm2(time*0.05+mo.x*3.); vec3 ta = vec3(0); vec3 ww = normalize(ta - ro); vec3 uu = normalize(cross(vec3(0.,1.,0.), ww)); vec3 vv = normalize(cross(ww, uu)); vec3 rd = normalize(p.x*uu + p.y*vv + 1.5*ww); float rz = march(ro,rd); vec3 col = texture(iChannel0, rd).rgb; if ( rz < FAR ) { vec3 pos = ro+rz*rd; float crv; vec3 nor = norcurv(pos, crv); crv = curv2(pos, 0.01); vec3 ligt = normalize( vec3(.0, 1., 0.) ); float dif = clamp(dot( nor, ligt ), 0., 1.); float bac = clamp( dot( nor, -ligt), 0.0, 1.0 ); float spe = pow(clamp( dot( reflect(rd,nor), ligt ), 0.0, 1.0 ),400.); float fre = pow( clamp(1.0+dot(nor,rd),0.0,1.0), 2.0 ); vec3 brdf = vec3(0.10,0.11,0.13); brdf += bac*vec3(0.1); brdf += dif*vec3(1.00,0.90,0.60); col = abs(sin(vec3(0.2,0.5,.9)+clamp(crv*80.,0.,1.)*1.2)); col = mix(col,texture(iChannel0,reflect(rd,nor)).rgb,.5); col = col*brdf + col*spe +.3*fre*mix(col,vec3(1),0.5); col *= smoothstep(-1.,-.9,sin(crv*200.))*0.15+0.85; } fragColor = vec4( col, 1.0 ); }

// My version with traditional coloring of the curvature:

https://www.shadertoy.com/view/wsGfRW

A version using the glsl instruction fwidth:
https://www.shadertoy.com/view/ldd3DH

//////////////////////
        // Getting the curvature from 2 derivates : position and normal.
        vec3 gw = fwidth(g);
        vec3 pw = fwidth(rp);
        
        float wfcurvature = length(gw) / length(pw);
        wfcurvature = smoothstep(0.0, 30., wfcurvature);
        color = vec4(wfcurvature);

SDF and hessian stuff:

https://www.shadertoy.com/view/4stSRf
https://www.shadertoy.com/view/Mt3SWX

Laplacian to detect discontinuities in the field!

https://www.shadertoy.com/view/4djyzK
https://www.shadertoy.com/view/MsSfD3

No comments

All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Anti-spam question: