Raymarching Cheap Subsurface Scattering
<< Home

Raymarching Cheap Subsurface Scattering

It’s been a while, school’s been busy, yada yada yada. I’ve got a new phone! I am considering about migrating this site to 42yeah.me or 42yeah.com, which miraculously isn’t occupied by anyone! Long story short, I’m back! At least for a while. And as I am back, let’s do some Subsurface Scattering!

Introduction

Subsurface Scattering (SSS) itself in nature is a very complex aspect, as light photons shoots through surface and bounce inside it for a little bit, before shooting out again. This can make things around the outer edge much brighter, as if it’s partially transparent:

Skin SSS

And I am gonna go ahead and make just a scene like this:

mat2 rot2d(float rad) {
    float a = sin(rad), b = cos(rad);
    return mat2(
        vec2(b, a),
        vec2(-a, b)
    );
}

float sol(vec3 p) {
    return p.y;
}

float cube(vec3 p, vec3 off, float rot) {
    p -= off;
    p.xz = rot2d(rot) * p.xz;
    vec3 d = abs(p) - vec3(0.5, 0.5, 0.5);
    return length(max(d, 0.0)) + min(max(d.x, max(d.y, d.z)), 0.0);
}

vec2 map(vec3 p) {
    float shortest = 1000.0;
    float id = -1.0;
    
    float dist = sol(p);
    if (dist < shortest) { shortest = dist; id = 0.5; }
    
    dist = cube(p, vec3(0.0, 0.48, 0.0), radians(45.0));
    if (dist < shortest) { shortest = dist; id = 1.5; }
    
    dist = cube(p, vec3(1.45, 0.48, 0.0), radians(45.0));
    if (dist < shortest) { shortest = dist; id = 1.5; }
    
    return vec2(shortest, id);
}

vec2 trace(vec3 ro, vec3 rd) {
    float depth = 0.0;
    float id = -1.0;
    for (int i = 0; i < 200; i++) {
        vec2 info = map(ro + rd * depth);
        if (abs(info.x) < 1e-3) {
            id = info.y;
            break;
        }
        depth += info.x;
    }
    return vec2(depth, id);
}

vec3 getRd(vec3 ro, vec3 center, vec2 uv) {
    vec3 front = normalize(center - ro);
    vec3 right = normalize(cross(front, vec3(0.0, 1.0, 0.0)));
    vec3 up = normalize(cross(right, front));
    mat4 lookAt = mat4(
        vec4(right, 0.0),
        vec4(up, 0.0),
        vec4(front, 0.0),
        vec4(0.0, 0.0, 0.0, 1.0)
    );
    return normalize(vec3(lookAt * vec4(uv, 2.0, 1.0)));
}

vec3 getFloorColor(vec3 p) {
    p *= 2.0;
    vec3 baseColor = vec3(1.0, 1.0, 1.0);
    vec3 u = mod(floor(p), 2.0);
    return baseColor - clamp(abs(u.x - u.z), 0.0, 0.1);
}

vec3 getColor(float id, vec3 p, vec3 rd) {
    if (id < -0.5) {
        vec3 skyColor = vec3(0.67, 0.86, 0.98);
        return mix(skyColor, vec3(1.0), clamp(rd.y * 5.0, 0.0, 1.0));
    }
    if (id < 1.0) { return getFloorColor(p); }
    if (id < 2.0) { return vec3(1.0, 0.5, 0.0); }
    return vec3(1.0, 0.0, 0.0);
}

vec3 getNormal(vec3 p) {
    const float epsilon = 1e-3;
    float m = map(p).x;
    return normalize(vec3(
        m - map(vec3(p.x - epsilon, p.yz)).x,
        m - map(vec3(p.x, p.y - epsilon, p.z)).x,
        m - map(vec3(p.x, p.y, p.z - epsilon)).x
    ));
}

float getShadowIntensity(vec3 ro, vec3 rd) {
    float depth = 0.001;
    float res = 1.0;
    for (int i = 0; i < 25; i++) {
        float dist = map(ro + rd * depth).x;
        res = min(res, 20.0 * dist / depth);
        if (res < 1e-6) { break; }
        depth += clamp(dist, 0.001, 2.0);
    }
    return clamp(res, 0.0, 1.0);
}

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    
    vec3 ro = vec3(0.65, 2.5, 3.0);
    vec3 center = vec3(0.65, 0.5, 0.0);
    vec3 rd = getRd(ro, center, xy);
    
    vec2 info = trace(ro, rd);
    vec3 p = ro + rd * info.x;
    vec3 n = getNormal(p);
    
    vec3 lightDir = normalize(vec3(1.0, 1.0, 1.0));
     float ambient = 1.0;
    float diffuse = max(dot(n, lightDir), 0.0);
    float dome = 0.2 + 0.8 * clamp(n.y, 0.0, 1.0);
    float sol = 0.2 + 0.8 * clamp(-n.y, 0.0, 1.0);
    float back = max(dot(n, vec3(-lightDir.x, 0.0, -lightDir.z)), 0.0);
    float shadow = getShadowIntensity(p + n * 1e-3, lightDir);

    vec3 light = vec3(1.0);
    if (info.y > 0.0) {
        light = vec3(0.0);
        light += ambient * vec3(0.2, 0.2, 0.2) * shadow;
        light += diffuse * vec3(0.82, 0.80, 0.82) * shadow;
        light += dome * vec3(0.26, 0.32, 0.334);
        light += sol * vec3(0.3, 0.21, 0.23);
        light += back * vec3(0.2, 0.21, 0.23);
    }
    
    vec3 objColor = getColor(info.y, p, rd) * light;
    objColor = pow(objColor, vec3(0.4545));
    color = vec4(objColor, 1.0);
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Well, it looks nothing like a hand. But whatever! I just feel like adding two cubes right now. Let’s go!

Implementation

Nature’s way of implementing SSS is certainly not feasible in a tiny GPU. The light rays went through complex reflections to go to the other side. We have two ways to achieve sounding effects: shooting random rays (which kinda looks like reverse AO):

  • You get the shading point
  • You shoot random rays in the direction of the negative normal’s hemisphere
  • The more rays it is able to shoot out of a geometry, the more SSS it gets
  • If there are little to none rays that are able to shoot out of a geometry, the shading point should be deep in the current geometry, and thus gets no SSS

Sketch

and estimating thickness, which is what Shane’s shader does, and thus is going to be covered here. It’s kinda cheaper & quicker than shooting random rays. So how are we gonna do that?

We estimate thickness by rapidly shooting rays in the reverse direction of the normal. At first, the ray is quite short, so it is quite likely that it will end up in the geometry; but as the ray gets longer and longer, it becomes more and more possible that it will shoot out of the current geometry and end up in open space. And that’s exactly what this piece of code do!

mat2 rot2d(float rad) {
    float a = sin(rad), b = cos(rad);
    return mat2(
        vec2(b, a),
        vec2(-a, b)
    );
}

float sol(vec3 p) {
    return p.y;
}

float cube(vec3 p, vec3 off, float rot) {
    p -= off;
    p.xz = rot2d(rot) * p.xz;
    vec3 d = abs(p) - vec3(0.5, 0.5, 0.5);
    return length(max(d, 0.0)) + min(max(d.x, max(d.y, d.z)), 0.0);
}

vec2 map(vec3 p) {
    float shortest = 1000.0;
    float id = -1.0;
    
    float dist = sol(p);
    if (dist < shortest) { shortest = dist; id = 0.5; }
    
    dist = cube(p, vec3(0.0, 0.48, 0.0), radians(45.0));
    if (dist < shortest) { shortest = dist; id = 1.5; }
    
    dist = cube(p, vec3(1.45, 0.48, 0.0), radians(45.0));
    if (dist < shortest) { shortest = dist; id = 1.5; }
    
    return vec2(shortest, id);
}

vec2 trace(vec3 ro, vec3 rd) {
    float depth = 0.0;
    float id = -1.0;
    for (int i = 0; i < 200; i++) {
        vec2 info = map(ro + rd * depth);
        if (abs(info.x) < 1e-3) {
            id = info.y;
            break;
        }
        depth += info.x;
    }
    return vec2(depth, id);
}

vec3 getRd(vec3 ro, vec3 center, vec2 uv) {
    vec3 front = normalize(center - ro);
    vec3 right = normalize(cross(front, vec3(0.0, 1.0, 0.0)));
    vec3 up = normalize(cross(right, front));
    mat4 lookAt = mat4(
        vec4(right, 0.0),
        vec4(up, 0.0),
        vec4(front, 0.0),
        vec4(0.0, 0.0, 0.0, 1.0)
    );
    return normalize(vec3(lookAt * vec4(uv, 2.0, 1.0)));
}

vec3 getFloorColor(vec3 p) {
    p *= 2.0;
    vec3 baseColor = vec3(1.0, 1.0, 1.0);
    vec3 u = mod(floor(p), 2.0);
    return baseColor - clamp(abs(u.x - u.z), 0.0, 0.1);
}

vec3 getSkyColor(vec3 rd) {
    vec3 skyColor = vec3(0.67, 0.86, 0.98);
    return mix(skyColor, vec3(1.0), clamp(rd.y * 5.0, 0.0, 1.0));
}

vec3 getColor(float id, vec3 p, vec3 rd) {
    if (id < -0.5) {
        return getSkyColor(rd);
    }
    if (id < 1.0) { return getFloorColor(p); }
    if (id < 2.0) { return vec3(1.0, 0.5, 0.0); }
    return vec3(1.0, 0.0, 0.0);
}

vec3 rand3d(float i) {
    return fract(sin(vec3(
        i * 12.345,
        i * 61.12,
        i * 512.3
    )) * 4123.123) * 2.0 - 1.0;
}

vec3 getNormal(vec3 p) {
    const float epsilon = 1e-3;
    float m = map(p).x;
    return normalize(vec3(
        m - map(vec3(p.x - epsilon, p.yz)).x,
        m - map(vec3(p.x, p.y - epsilon, p.z)).x,
        m - map(vec3(p.x, p.y, p.z - epsilon)).x
    ));
}

float getShadowIntensity(vec3 ro, vec3 rd) {
    float depth = 0.001;
    float res = 1.0;
    for (int i = 0; i < 25; i++) {
        float dist = map(ro + rd * depth).x;
        res = min(res, 20.0 * dist / depth);
        if (res < 1e-6) { break; }
        depth += clamp(dist, 0.001, 2.0);
    }
    return clamp(res, 0.0, 1.0);
}

float thinness(vec3 p, vec3 n);

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    
    vec3 ro = vec3(cos(time) * 3.0, 2.5, sin(time) * 3.0);
    vec3 center = vec3(0.65, 0.5, 0.0);
    vec3 rd = getRd(ro, center, xy);
    
    vec2 info = trace(ro, rd);
    vec3 p = ro + rd * info.x;
    vec3 n = getNormal(p);
    
    vec3 lightDir = normalize(vec3(1.0, 1.0, 1.0));
    float ambient = 1.0;
    float diffuse = max(dot(n, lightDir), 0.0);
    float dome = 0.2 + 0.8 * clamp(n.y, 0.0, 1.0);
    float sol = 0.2 + 0.8 * clamp(-n.y, 0.0, 1.0);
    float back = max(dot(n, vec3(-lightDir.x, 0.0, -lightDir.z)), 0.0);
    float shadow = getShadowIntensity(p + n * 1e-3, lightDir);
    
    float th = thinness(p, n);
    color = vec4(th, th < 0.0, th == 0.0, 1.0);
    return;
}
float thinness(vec3 p, vec3 n) {
    float sub = 0.0;
    for (int i = 0; i < 30; i++) {
        float dist = float(i) * 0.4 / 30.0;
        sub += dist + map(p - n * dist).x;
    }
    sub /= 30.0;
    return sub;
}
Recompile | Click on the canvas to toggle running. It is paused by default.

The thinner a part is, the more it is possible to end up in open space. Thus we add those results up, and average it. We can see that indeed the thinnest parts are reddest. Which is good!

After getting the thinness value, we can utilize it to create the sss effect. It is straightforward:

light += th * getSkyColor(vec3(rd.x, -rd.y, rd.z)) * 0.9;

As a matter of fact, the color is really your pick. It could also be the color of the stuff behind it, or in my case, I use the sky’s color. It does not matter. Just test until you find a suitable color! And this is how it looks like when we put all stuffs we know for now together:

mat2 rot2d(float rad) {
    float a = sin(rad), b = cos(rad);
    return mat2(
        vec2(b, a),
        vec2(-a, b)
    );
}

float sol(vec3 p) {
    return p.y;
}

float cube(vec3 p, vec3 off, float rot) {
    p -= off;
    p.xz = rot2d(rot) * p.xz;
    vec3 d = abs(p) - vec3(0.5, 0.5, 0.5);
    return length(max(d, 0.0)) + min(max(d.x, max(d.y, d.z)), 0.0);
}

vec2 map(vec3 p) {
    float shortest = 1000.0;
    float id = -1.0;
    
    float dist = sol(p);
    if (dist < shortest) { shortest = dist; id = 0.5; }
    
    dist = cube(p, vec3(0.0, 0.48, 0.0), radians(45.0));
    if (dist < shortest) { shortest = dist; id = 1.5; }
    
    dist = cube(p, vec3(1.45, 0.48, 0.0), radians(45.0));
    if (dist < shortest) { shortest = dist; id = 1.5; }
    
    return vec2(shortest, id);
}

vec2 trace(vec3 ro, vec3 rd) {
    float depth = 0.0;
    float id = -1.0;
    for (int i = 0; i < 200; i++) {
        vec2 info = map(ro + rd * depth);
        if (abs(info.x) < 1e-3) {
            id = info.y;
            break;
        }
        depth += info.x;
    }
    return vec2(depth, id);
}

vec3 getRd(vec3 ro, vec3 center, vec2 uv) {
    vec3 front = normalize(center - ro);
    vec3 right = normalize(cross(front, vec3(0.0, 1.0, 0.0)));
    vec3 up = normalize(cross(right, front));
    mat4 lookAt = mat4(
        vec4(right, 0.0),
        vec4(up, 0.0),
        vec4(front, 0.0),
        vec4(0.0, 0.0, 0.0, 1.0)
    );
    return normalize(vec3(lookAt * vec4(uv, 2.0, 1.0)));
}

vec3 getFloorColor(vec3 p) {
    p *= 2.0;
    vec3 baseColor = vec3(1.0, 1.0, 1.0);
    vec3 u = mod(floor(p), 2.0);
    return baseColor - clamp(abs(u.x - u.z), 0.0, 0.1);
}

vec3 getSkyColor(vec3 rd) {
    vec3 skyColor = vec3(0.67, 0.86, 0.98);
    return mix(skyColor, vec3(1.0), clamp(rd.y * 5.0, 0.0, 1.0));
}

vec3 getColor(float id, vec3 p, vec3 rd) {
    if (id < -0.5) {
        return getSkyColor(rd);
    }
    if (id < 1.0) { return getFloorColor(p); }
    if (id < 2.0) { return vec3(1.0, 0.5, 0.0); }
    return vec3(1.0, 0.0, 0.0);
}

vec3 rand3d(float i) {
    return fract(sin(vec3(
        i * 12.345,
        i * 61.12,
        i * 512.3
    )) * 4123.123) * 2.0 - 1.0;
}

vec3 getNormal(vec3 p) {
    const float epsilon = 1e-3;
    float m = map(p).x;
    return normalize(vec3(
        m - map(vec3(p.x - epsilon, p.yz)).x,
        m - map(vec3(p.x, p.y - epsilon, p.z)).x,
        m - map(vec3(p.x, p.y, p.z - epsilon)).x
    ));
}

float getShadowIntensity(vec3 ro, vec3 rd) {
    float depth = 0.001;
    float res = 1.0;
    for (int i = 0; i < 25; i++) {
        float dist = map(ro + rd * depth).x;
        res = min(res, 20.0 * dist / depth);
        if (res < 1e-6) { break; }
        depth += clamp(dist, 0.001, 2.0);
    }
    return clamp(res, 0.0, 1.0);
}
float thinness(vec3 p, vec3 n) {
    float sub = 0.0;
    for (int i = 0; i < 30; i++) {
        float dist = float(i) * 0.4 / 30.0;
        sub += dist + map(p - n * dist).x;
    }
    sub /= 30.0;
    return sub;
}

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    
    vec3 ro = vec3(cos(time) * 3.0, 2.5, sin(time) * 3.0);
    vec3 center = vec3(0.65, 0.5, 0.0);
    vec3 rd = getRd(ro, center, xy);
    
    vec2 info = trace(ro, rd);
    vec3 p = ro + rd * info.x;
    vec3 n = getNormal(p);
    
    vec3 lightDir = normalize(vec3(1.0, 1.0, 1.0));
    float ambient = 1.0;
    float diffuse = max(dot(n, lightDir), 0.0);
    float dome = 0.2 + 0.8 * clamp(n.y, 0.0, 1.0);
    float sol = 0.2 + 0.8 * clamp(-n.y, 0.0, 1.0);
    float back = max(dot(n, vec3(-lightDir.x, 0.0, -lightDir.z)), 0.0);
    float shadow = getShadowIntensity(p + n * 1e-3, lightDir);
    float th = thinness(p, n);

    vec3 light = vec3(1.0);
    if (info.y > 0.0) {
        light = vec3(0.0);
        light += ambient * vec3(0.2, 0.2, 0.2) * shadow;
        light += diffuse * vec3(0.82, 0.80, 0.82) * shadow;
        light += dome * vec3(0.26, 0.32, 0.334);
        light += sol * vec3(0.3, 0.21, 0.23);
        light += back * vec3(0.2, 0.21, 0.23);
        light += th * getSkyColor(vec3(rd.x, -rd.y, rd.z)) * 0.9;
    }
    
    vec3 objColor = getColor(info.y, p, rd) * light;
    objColor = pow(objColor, vec3(0.4545));
    color = vec4(objColor, 1.0);
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Not bad! Notice the edges became brighter and the cube feels jellier now? However, there are artifacts. They are quite obvious when facing back: it looks like a dark “X”. There are ways to weaken it, though (as far as I know) it will also weaken the sss effect. One way is to scale the result, so the further the ray, the weaker its influence is:

mat2 rot2d(float rad) {
    float a = sin(rad), b = cos(rad);
    return mat2(
        vec2(b, a),
        vec2(-a, b)
    );
}

float sol(vec3 p) {
    return p.y;
}

float cube(vec3 p, vec3 off, float rot) {
    p -= off;
    p.xz = rot2d(rot) * p.xz;
    vec3 d = abs(p) - vec3(0.5, 0.5, 0.5);
    return length(max(d, 0.0)) + min(max(d.x, max(d.y, d.z)), 0.0);
}

vec2 map(vec3 p) {
    float shortest = 1000.0;
    float id = -1.0;
    
    float dist = sol(p);
    if (dist < shortest) { shortest = dist; id = 0.5; }
    
    dist = cube(p, vec3(0.0, 0.48, 0.0), radians(45.0));
    if (dist < shortest) { shortest = dist; id = 1.5; }
    
    dist = cube(p, vec3(1.45, 0.48, 0.0), radians(45.0));
    if (dist < shortest) { shortest = dist; id = 1.5; }
    
    return vec2(shortest, id);
}

vec2 trace(vec3 ro, vec3 rd) {
    float depth = 0.0;
    float id = -1.0;
    for (int i = 0; i < 200; i++) {
        vec2 info = map(ro + rd * depth);
        if (abs(info.x) < 1e-3) {
            id = info.y;
            break;
        }
        depth += info.x;
    }
    return vec2(depth, id);
}

vec3 getRd(vec3 ro, vec3 center, vec2 uv) {
    vec3 front = normalize(center - ro);
    vec3 right = normalize(cross(front, vec3(0.0, 1.0, 0.0)));
    vec3 up = normalize(cross(right, front));
    mat4 lookAt = mat4(
        vec4(right, 0.0),
        vec4(up, 0.0),
        vec4(front, 0.0),
        vec4(0.0, 0.0, 0.0, 1.0)
    );
    return normalize(vec3(lookAt * vec4(uv, 2.0, 1.0)));
}

vec3 getFloorColor(vec3 p) {
    p *= 2.0;
    vec3 baseColor = vec3(1.0, 1.0, 1.0);
    vec3 u = mod(floor(p), 2.0);
    return baseColor - clamp(abs(u.x - u.z), 0.0, 0.1);
}

vec3 getSkyColor(vec3 rd) {
    vec3 skyColor = vec3(0.67, 0.86, 0.98);
    return mix(skyColor, vec3(1.0), clamp(rd.y * 5.0, 0.0, 1.0));
}

vec3 getColor(float id, vec3 p, vec3 rd) {
    if (id < -0.5) {
        return getSkyColor(rd);
    }
    if (id < 1.0) { return getFloorColor(p); }
    if (id < 2.0) { return vec3(1.0, 0.5, 0.0); }
    return vec3(1.0, 0.0, 0.0);
}

vec3 rand3d(float i) {
    return fract(sin(vec3(
        i * 12.345,
        i * 61.12,
        i * 512.3
    )) * 4123.123) * 2.0 - 1.0;
}

vec3 getNormal(vec3 p) {
    const float epsilon = 1e-3;
    float m = map(p).x;
    return normalize(vec3(
        m - map(vec3(p.x - epsilon, p.yz)).x,
        m - map(vec3(p.x, p.y - epsilon, p.z)).x,
        m - map(vec3(p.x, p.y, p.z - epsilon)).x
    ));
}

float getShadowIntensity(vec3 ro, vec3 rd) {
    float depth = 0.001;
    float res = 1.0;
    for (int i = 0; i < 25; i++) {
        float dist = map(ro + rd * depth).x;
        res = min(res, 20.0 * dist / depth);
        if (res < 1e-6) { break; }
        depth += clamp(dist, 0.001, 2.0);
    }
    return clamp(res, 0.0, 1.0);
}

float thinness(vec3 p, vec3 n);

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    
    vec3 ro = vec3(cos(time) * 3.0, 2.5, sin(time) * 3.0);
    vec3 center = vec3(0.65, 0.5, 0.0);
    vec3 rd = getRd(ro, center, xy);
    
    vec2 info = trace(ro, rd);
    vec3 p = ro + rd * info.x;
    vec3 n = getNormal(p);
    
    vec3 lightDir = normalize(vec3(1.0, 1.0, 1.0));
    float ambient = 1.0;
    float diffuse = max(dot(n, lightDir), 0.0);
    float dome = 0.2 + 0.8 * clamp(n.y, 0.0, 1.0);
    float sol = 0.2 + 0.8 * clamp(-n.y, 0.0, 1.0);
    float back = max(dot(n, vec3(-lightDir.x, 0.0, -lightDir.z)), 0.0);
    float shadow = getShadowIntensity(p + n * 1e-3, lightDir);
    float th = thinness(p, n);

    color = vec4(th, th < 0.0, th == 0.0, 1.0);
    return;
}
float thinness(vec3 p, vec3 n) {
    float sub = 0.0;
    float scl = 1.0;
    for (int i = 0; i < 30; i++) {
        float dist = float(i) * 0.4 / 30.0;
        sub += dist + map(p - n * dist).x * scl;
        scl *= 0.9;
    }
    sub /= 30.0;
    return sub;
}
Recompile | Click on the canvas to toggle running. It is paused by default.

However as it becomes more uniform, the sss effect becomes more subtle, too. Well, I still haven’t figure this out yet, so I am open to suggestions!

Conclusion

And that concludes the SSS. It is cheap, and I know this is not the best practice of SSS, and I am still learning. Hopefully one day I can come back and optimize it! And for now, cheerio!