ODOS丨走向珊迪！RAY-MARCHING

by ERIN.Z, 2022-02-20

说点儿题外话，Unity3d的中文注册名是优美缔和优三缔....... 也太low了真的........

RayMarching

The term ray marching is more broad and refers to methods in which simulated rays are traversed iteratively, effectively dividing each ray into smaller ray segments, sampling some function at each step.

In SDF ray marching, or sphere tracing, an intersection point is approximated between the ray and a surface defined by a signed distance function (SDF). The SDF is evaluated for each iteration in order to be able take as large steps as possible without missing any part of the surface. A threshold is used to cancel further iteration when a point has reached that is close enough to the surface.

在shadertoy中的raymarching技术是基于SDF(Signed Distance Function)函数的: 从观察点发射一条射线，顺着射线的方向每次步进SDF的长度，直到SDF的长度小于阈值（击中）或步进的次数超过最大值（未击中）。 the art of code做了非常详尽的讲解，可以参考这两个视频：Ray Marching for Dummies!，Ray Marching Simple Shapes。

3维空间坐标

我们基于屏幕的uv构造出一个三维的世界坐标系，以图像平面中心为原点，向屏幕内为z方向。 float f(vec3 p){ vec4 sph = vec4(0,1,3,1);//shphe xyz-O w-r float ds = length(p-sph.xyz)-sph.w;//sdf sphere float dp = p.y;//sdf plane y=0 return min(ds,dp); } 这是一个非常简单的SDF函数，场景由一个球体和一个平面构成，函数返回输入点到这两个物体的距离最小值。 float rayMarch(in vec3 ro, in vec3 rd) { float t = TMIN; for(int i = 0; i < RAYMARCH_TIME ; i++) { vec3 p = ro + t * rd; float d = f(p); t += d; if(d < PRECISION || t > TMAX) break; } return t; } 我们将摄像机(ro)置于(0,1,0)，rd设置为normalize(vec3(uv,1))，获得的t值可视化效果如图：

Get Normal

为了使用我们熟悉的光照模型，首先我们要获取法线。对于某表示空间曲面的二元函数z=f(x,y)：

偏导数(Partial derivative)表示固定面上一点的切线斜率。偏导数f'x(x0,y0) 表示固定面上一点对x 轴的切线斜率；偏导数f'y(x0,y0) 表示固定面上一点对y 轴的切线斜率。

但我们这里求偏导的对象是SDF函数，所以对xyz分别求偏导的结果其实就是曲面的法线方向。 vec3 calcNormal( in vec3 p ) // for function f(p) { const float eps = 0.0001; // or some other value const vec2 h = vec2(eps,0); return normalize( vec3(f(p+h.xyy) - f(p-h.xyy), f(p+h.yxy) - f(p-h.yxy), f(p+h.yyx) - f(p-h.yyx) ) ); } 这个算法以输入点为中心进行计算，即求(f(p+h)-f(p-h))/2h，因为会进行归一化，所以"/2h"的部分可以省略。但这样就需要计算6次f(x)，我们也可以放弃部分精度，使用前向计算，这样就只需要4次计算。 return normalize( vec3(f(p+h.xyy) - f(p), f(p+h.yxy) - f(p), f(p+h.yyx) - f(p) ) ); iq老师提出了使用四面体顶点取样的方式，使用4次f(x)同样可以取得中心差异的法线信息，详细解读可以参考normalsSDF。 // https://www.iquilezles.org/www/articles/normalsSDF/normalsSDF.htm vec3 calcNormal( in vec3 & p ) // for function f(p) { const float h = 0.0001; // replace by an appropriate value const vec2 k = vec2(1,-1); return normalize( k.xyyf( p + k.xyyh ) + k.yyxf( p + k.yyxh ) + k.yxyf( p + k.yxyh ) + k.xxxf( p + k.xxxh ) ); }

Get Shadow

在raymarching中获取硬边缘阴影区域非常简单，只要从hit的位置向光源方向发射一条raymarch射线，如果距离小于点到光源的距离，说明该位置在阴影中。

float shadow = rayMarch(p,l);
if(shadow<length(lightPos-p)) dif*=.1;

注意，我这里raymarching的函数已经有了一个初始TMIN的offset，可以防止这里计算时第一步就小于阈值。 (远处白色应该是因为超出最大步进值了，不过因为我们是对diff做乘法所以不影响~）

Shading

Lambert模型，来：

vec3 render(vec2 uv){
    vec3 lightPos = vec3(5.*cos(iTime), 5., 5.*sin(iTime)+3.);//light
    vec3 ro = vec3(0,1,0);//camera
    vec3 rd = normalize(vec3(uv,1));

    vec3 color = vec3(0);//background
    vec3 amb = vec3(.1);//ambient color
    vec3 diffusecol = vec3(1.,1.,1.);//diffuse color

    float t = rayMarch(ro,rd);//raymarching

    vec3 p = ro+t*rd;
    vec3 n = calcNormal(p);
    vec3 l = normalize(lightPos-p);//lightDir
    vec3 dif = clamp(dot(l,n),0.,1.)*diffusecol;//lambert

    float shadow = rayMarch(p,l);
    if(shadow<length(lightPos-p)) dif*=.1;

    color = amb + dif;

    return color;
}

主函数里再进行一下AA，AA的原理与上一篇是一样的： for(int m = 0; m < AA; m++) { for(int n = 0; n < AA; n++) { vec2 offset = 2. (vec2(float(m), float(n)) / float(AA) - .5); vec2 uv = getuv(fragCoord + offset); color += render(uv); } } color /= float(AAAA);

More Primetives

float sdCapsule(vec3 p,vec3 a, vec3 b,float r){
    vec3 ab = b - a;
    vec3 ap = p - a;
    float t = dot(ab, ap)/dot(ab,ab);
    t = clamp(t,0.,1.);
    vec3 c = a +t*ab;
    return length(p-c)-r;
}
float sdTorus(vec3 p,vec3 o, float r1,float r2){
    p -= o;
    return length( vec2(length(vec2(p.x,p.z))-r1,p.y) ) - r2;
}
float sdBox(vec3 p, vec3 o, vec3 s){
    p-=o;
    return length(max(abs(p)-s,0.));
}
float sdCylinder(vec3 p, vec3 a, vec3 b, float r) {
vec3 ab = b-a;
    vec3 ap = p-a;
    float t = dot(ab, ap) / dot(ab, ab);
    vec3 c = a + t*ab;

    float x = length(p-c)-r;
    float y = (abs(t-.5)-.5)*length(ab);
    float e = length(max(vec2(x, y), 0.));
    float i = min(max(x, y), 0.);

    return e+i;
}

详细的解释可以参考https://www.youtube.com/watch?v=Ff0jJyyiVyw

参考资料

by ERIN.Z