Curvature of a Distance field / Implicit surface
- 11/2020
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