Ambient Occlusion in Raymarching
<< Home

Ambient Occlusion in Raymarching

Ambient Occlusion (AO) is another key to photorealism, and it adds another great deal to the scene. It mimicks diffuse reflection. I used to implement it in OpenGL, but I was less familiar back then. Well, am I any better now? We will see!

Introduction

Let’s take a look, again, at my recent raymarching lab rat, Fruxis:

Fruxis

See the darkened corers near the walls? Now that’s ambient occlusion. Removing it will make the scene looks a tad more artificial:

Withou AO

Even though it’s not too bad without AO, the melon already seems like it’s floating, and the walls seems weird. So, are you ready? Let’s venture into the brave world of Gloabl Illumination again! First, a basic scene:

const vec3 lightPos = vec3(1.0, 3.0, 1.0);
const vec3 lightDir = normalize(lightPos);

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

float ball(vec3 p, vec3 off) {
    p -= off;
    return length(p) - 0.5;
}

float wall(vec3 p) {
    float d = min(p.x + 2.0, p.z + 2.0);
    return d;
}

vec3 getWallColor(vec3 p) {
    p *= 2.0;
    vec3 baseColor = vec3(1.0);
    vec3 u = mod(floor(p), 2.0);
    return baseColor - vec3(abs(u.y - max(u.x, u.z))) * 0.095;
}

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

vec2 map(vec3 p) {
    float closest = 1000.0;
    float id = -1.0;

    float dist = sol(p);
    if (dist < closest) { closest = dist; id = 0.5; }

    dist = ball(p, vec3(0.0, 0.1, 0.0));
    if (dist < closest) { closest = dist; id = 1.0; }

    dist = wall(p);
    if (dist < closest) { closest = dist; id = 2.0; }

    return vec2(closest, id);
}

vec2 intersect(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 (info.x < 0.001) {
            id = info.y;
            break;
        }
        depth += info.x;
    }
    return vec2(depth, id);
}

vec3 getColor(float id, vec3 p) {
    if (id < -0.5) { return vec3(0.0); }
    if (id < 1.0) { return getFloorColor(p); }
    if (id < 2.0) { return vec3(0.4, 0.2, 0.1) * 2.0; }
    if (id < 3.0) { return getWallColor(p); }
    return vec3(1.0, 0.0, 0.0);
}

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

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    vec3 ro = vec3(3.0 + sin(time) * 0.5, 3.0, 2.0 + sin(time * 0.2) * 0.5);
    vec3 center = vec3(0.0, 0.0, 0.0);
    
    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)
    );
    vec3 rd = normalize(vec3(lookAt * vec4(xy, 2.0, 1.0)));
    vec2 info = intersect(ro, rd);
    vec3 pos = ro + rd * info.x;
    vec3 n = getNormal(pos);

    float ambient = 1.0;
    float diffuse = max(dot(n, lightDir), 0.0);
    float back = max(dot(n, vec3(-lightDir.x, 0.0, -lightDir.z)), 0.0);
    float dome = 0.2 + 0.8 * clamp(rd.y, 0.0, 1.0);
    float sol = 0.2 + 0.8 * clamp(-rd.y, 0.0, 1.0);
    
    vec3 light = vec3(0.0);
    light += ambient * vec3(0.1, 0.1, 0.2);
    light += diffuse * vec3(0.65, 0.6, 0.4);
    light += back * vec3(0.3, 0.2, 0.1);
    light += dome * vec3(0.5, 0.5, 0.7);
    light += sol * vec3(0.3, 0.3, 0.32);
    
    vec3 objColor = getColor(info.y, pos) * 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.

Implementation

Now that we have a pretty basic scene, which consists of walls, floors and a ball, it’s time to talk about ambient occlusion!

First, you should know about diffuse reflection, which is happening everywhere in our lives (as nothing could be completely smooth). Photons were shot from some sort of light source to a rough surface, scatters randomly (because the surface was rough), and bounces into other obstacles and got absorbed by them, never to be caught by the eye again.

Diffuse reflection

This is (nearly exactly) what ambient occlusion is simulating. We really don’t need a light source, though; here’s how it’s gonna work, Monte Carlo way.

  1. Get the sample point (a smidge off pos + n * depth)

Sample point

  1. Randomly shoots out rays at that point

Shoots rays

  1. Power those ray vectors by three so they appear to be ball-like pattern
  2. Reverse the directions of the rays in the opposite direction to the sampled normal; we do this by dotting the ray with normal and getting the sign

Reverse

  1. Sample those rays.

Sample

By generating randomly reflected rays and estimating whether a majority of those rays ends up in obstacles or not, we can almost correctly estimate the ambient occlusion value of said point:

const vec3 lightPos = vec3(1.0, 3.0, 1.0);
const vec3 lightDir = normalize(lightPos);

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

float ball(vec3 p, vec3 off) {
    p -= off;
    return length(p) - 0.5;
}

float wall(vec3 p) {
    float d = min(p.x + 2.0, p.z + 2.0);
    return d;
}

vec3 getWallColor(vec3 p) {
    p *= 2.0;
    vec3 baseColor = vec3(1.0);
    vec3 u = mod(floor(p), 2.0);
    return baseColor - vec3(abs(u.y - max(u.x, u.z))) * 0.095;
}

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

vec2 map(vec3 p) {
    float closest = 1000.0;
    float id = -1.0;

    float dist = sol(p);
    if (dist < closest) { closest = dist; id = 0.5; }

    dist = ball(p, vec3(0.0, 0.1, 0.0));
    if (dist < closest) { closest = dist; id = 1.0; }

    dist = wall(p);
    if (dist < closest) { closest = dist; id = 2.0; }

    return vec2(closest, id);
}

vec2 intersect(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 (info.x < 0.001) {
            id = info.y;
            break;
        }
        depth += info.x;
    }
    return vec2(depth, id);
}

vec3 getColor(float id, vec3 p) {
    if (id < -0.5) { return vec3(0.0); }
    if (id < 1.0) { return getFloorColor(p); }
    if (id < 2.0) { return vec3(0.4, 0.2, 0.1) * 2.0; }
    if (id < 3.0) { return getWallColor(p); }
    return vec3(1.0, 0.0, 0.0);
}

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

vec3 rand3(float i) {
    return fract(sin(vec3(
         i * 123.456,
         i * 61.8,
         i * 618.21
    )) * 41234.56) * 2.0 - 1.0;
}

float ambientOcclusion(vec3 p, vec3 n);

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    vec3 ro = vec3(3.0 + sin(time) * 0.5, 3.0, 2.0 + sin(time * 0.2) * 0.5);
    vec3 center = vec3(0.0, 0.0, 0.0);
    
    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)
    );
    vec3 rd = normalize(vec3(lookAt * vec4(xy, 2.0, 1.0)));
    vec2 info = intersect(ro, rd);
    vec3 pos = ro + rd * info.x;
    vec3 n = getNormal(pos);

    float ambient = 1.0;
    float diffuse = max(dot(n, lightDir), 0.0);
    float back = max(dot(n, vec3(-lightDir.x, 0.0, -lightDir.z)), 0.0);
    float dome = 0.2 + 0.8 * clamp(rd.y, 0.0, 1.0);
    float sol = 0.2 + 0.8 * clamp(-rd.y, 0.0, 1.0);
    float ao = ambientOcclusion(pos, n);
    
    vec3 light = vec3(0.0);
    light += ambient * vec3(0.1, 0.1, 0.2);
    light += diffuse * vec3(0.65, 0.6, 0.4) * ao;
    light += back * vec3(0.3, 0.2, 0.1);
    light += dome * vec3(0.5, 0.5, 0.7);
    light += sol * vec3(0.3, 0.3, 0.32);
    
    vec3 objColor = getColor(info.y, pos) * light;
    objColor = pow(objColor, vec3(0.4545));
    color = vec4(objColor, 1.0);
}
float ambientOcclusion(vec3 p, vec3 n) {
    float totalAO = 0.0;
    for (int i = 0; i < 20; i++) {
        vec3 ray = rand3(float(i * 92));
        ray = pow(ray, vec3(3.0));
        ray *= sign(dot(ray, n));
        totalAO += clamp(map(p + n * 0.001 + 0.15 * ray).x * 48.0, 0.0, 1.0);
    }
    totalAO /= 20.0; // normalize
    return clamp(totalAO * totalAO, 0.0, 1.0);
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Notice the very obvious lines along the walls now? It is now quite obvious that there’s a wall there. Also the ball gains some kind of shadow - but not too much. We still need a shadow calculating algorithm!

totalAO, the ambient occlusion-ness value was averaged in the end, so we know where did most of those rays end up. The 48.0 above is some kind of factor: the less, the laxer, which means less ambient occlusion value, which means generally darker scene and a fuzzier ambient occlusion result. Other than that, everything’s fine.

Also after applying ambient occlusion to the scene, it becomes noticably darker. That’s because those rays will occasionally sample itself, so it is unavoidable. So we need to adjust some colors to brighten it up again. Here’s an (extremely) crude way: multiplying the output lighting calculations.

const vec3 lightPos = vec3(1.0, 3.0, 1.0);
const vec3 lightDir = normalize(lightPos);

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

float ball(vec3 p, vec3 off) {
    p -= off;
    return length(p) - 0.5;
}

float wall(vec3 p) {
    float d = min(p.x + 2.0, p.z + 2.0);
    return d;
}

vec3 getWallColor(vec3 p) {
    p *= 2.0;
    vec3 baseColor = vec3(1.0);
    vec3 u = mod(floor(p), 2.0);
    return baseColor - vec3(abs(u.y - max(u.x, u.z))) * 0.095;
}

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

vec2 map(vec3 p) {
    float closest = 1000.0;
    float id = -1.0;

    float dist = sol(p);
    if (dist < closest) { closest = dist; id = 0.5; }

    dist = ball(p, vec3(0.0, 0.1, 0.0));
    if (dist < closest) { closest = dist; id = 1.0; }

    dist = wall(p);
    if (dist < closest) { closest = dist; id = 2.0; }

    return vec2(closest, id);
}

vec2 intersect(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 (info.x < 0.001) {
            id = info.y;
            break;
        }
        depth += info.x;
    }
    return vec2(depth, id);
}

vec3 getColor(float id, vec3 p) {
    if (id < -0.5) { return vec3(0.0); }
    if (id < 1.0) { return getFloorColor(p); }
    if (id < 2.0) { return vec3(0.4, 0.2, 0.1) * 2.0; }
    if (id < 3.0) { return getWallColor(p); }
    return vec3(1.0, 0.0, 0.0);
}

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

vec3 rand3(float i) {
    return fract(sin(vec3(
         i * 123.456,
         i * 61.8,
         i * 618.21
    )) * 41234.56) * 2.0 - 1.0;
}

float ambientOcclusion(vec3 p, vec3 n) {
    float totalAO = 0.0;
    for (int i = 0; i < 20; i++) {
        vec3 ray = rand3(float(i * 92));
        ray = pow(ray, vec3(3.0));
        ray *= sign(dot(ray, n));
        totalAO += clamp(map(p + n * 0.001 + 0.15 * ray).x * 48.0, 0.0, 1.0);
    }
    totalAO /= 20.0; // normalize
    return clamp(totalAO * totalAO, 0.0, 1.0);
}
void main() {
    vec2 xy = uv * 2.0 - 1.0;
    vec3 ro = vec3(3.0 + sin(time) * 0.5, 3.0, 2.0 + sin(time * 0.2) * 0.5);
    vec3 center = vec3(0.0, 0.0, 0.0);
    
    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)
    );
    vec3 rd = normalize(vec3(lookAt * vec4(xy, 2.0, 1.0)));
    vec2 info = intersect(ro, rd);
    vec3 pos = ro + rd * info.x;
    vec3 n = getNormal(pos);

    float ambient = 1.0;
    float diffuse = max(dot(n, lightDir), 0.0);
    float back = max(dot(n, vec3(-lightDir.x, 0.0, -lightDir.z)), 0.0);
    float dome = 0.2 + 0.8 * clamp(rd.y, 0.0, 1.0);
    float sol = 0.2 + 0.8 * clamp(-rd.y, 0.0, 1.0);
    float ao = ambientOcclusion(pos, n);
    
    // Make them brighter: just multiply them!
    vec3 light = vec3(0.0);
    light += ambient * vec3(0.3, 0.3, 0.3) * 1.1;
    light += diffuse * vec3(0.65, 0.6, 0.4) * ao * 1.2;
    light += back * vec3(0.3, 0.2, 0.1) * 1.1;
    light += dome * vec3(0.5, 0.5, 0.7) * 1.1;
    light += sol * vec3(0.3, 0.3, 0.32) * 1.1;
    
    vec3 objColor = getColor(info.y, pos) * 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.

Conclusion

And there we have it! Ambient occlusion, which is surprisingly intuitive. This guy, paired with shadow, and a bunch of other lighting calculations, could make a fairly good-looking scene. I am not saying it’s photorealistic, because it is not; but still!

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

float ball(vec3 p, vec3 off) {
    p -= off;
    return length(p) - 0.5;
}

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

float wall(vec3 p) {
    float d = min(p.x + 2.0, p.z + 2.0);
    return d;
}

vec3 getWallColor(vec3 p) {
    p *= 2.0;
    vec3 baseColor = vec3(1.0);
    vec3 u = mod(floor(p), 2.0);
    return baseColor - vec3(abs(u.y - max(u.x, u.z))) * 0.095;
}

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

vec2 map(vec3 p) {
    float closest = 1000.0;
    float id = -1.0;

    float dist = sol(p);
    if (dist < closest) { closest = dist; id = 0.5; }

    dist = ball(p, vec3(1.0, 1.4, -0.7));
    if (dist < closest) { closest = dist; id = 1.0; }
    
    dist = ball(p, vec3(0.0, 0.4, 0.3));
    if (dist < closest) { closest = dist; id = 1.0; }
    
    dist = cube(p, vec3(1.0, 0.5, -0.7));
    if (dist < closest) { closest = dist; id = 1.0; }

    dist = wall(p);
    if (dist < closest) { closest = dist; id = 2.0; }

    return vec2(closest, id);
}

vec2 intersect(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 (info.x < 0.001) {
            id = info.y;
            break;
        }
        depth += info.x;
    }
    return vec2(depth, id);
}

vec3 getColor(float id, vec3 p) {
    if (id < -0.5) { return vec3(0.0); }
    if (id < 1.0) { return getFloorColor(p); }
    if (id < 2.0) { return vec3(0.4, 0.2, 0.1) * 2.0; }
    if (id < 3.0) { return getWallColor(p); }
    return vec3(1.0, 0.0, 0.0);
}

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

vec3 rand3(float i) {
    return fract(sin(vec3(
         i * 123.456,
         i * 61.8,
         i * 618.21
    )) * 41234.56) * 2.0 - 1.0;
}

float ambientOcclusion(vec3 p, vec3 n) {
    float totalAO = 0.0;
    for (int i = 0; i < 20; i++) {
        vec3 ray = rand3(float(i * 92));
        ray = pow(ray, vec3(3.0));
        ray *= sign(dot(ray, n));
        totalAO += clamp(map(p + n * 0.001 + 0.15 * ray).x * 48.0, 0.0, 1.0);
    }
    totalAO /= 20.0;
    return clamp(totalAO * totalAO, 0.0, 1.0);
}

float getShadowIntensity(vec3 ro, vec3 rd) {
    float res = 1.0;
    float depth = 0.001;
    for (int i = 0; i < 35; i++) {
        vec2 info = map(ro + rd * depth);
        res = min(res, clamp(20.0 * info.x / depth, 0.0, 1.0));
        if (res < 0.001) { break; }
        depth += info.x;
    }
    return res;
}

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    vec3 ro = vec3(3.0 + sin(time) * 0.5, 3.0, 2.0 + sin(time * 0.2) * 0.5);
    vec3 center = vec3(0.0, 0.0, 0.0);
    
    vec3 lightPos = vec3((sin(time) * 1.0 + 2.0) * (sin(time) * 0.5 + 0.6), 3.0, cos(time * 0.5) * 1.0 + 5.0);
    vec3 lightDir = normalize(lightPos);

    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)
    );
    vec3 rd = normalize(vec3(lookAt * vec4(xy, 2.0, 1.0)));
    vec2 info = intersect(ro, rd);
    vec3 pos = ro + rd * info.x;
    vec3 n = getNormal(pos);

    float ambient = 1.0;
    float diffuse = max(dot(n, lightDir), 0.0);
    float back = max(dot(n, vec3(-lightDir.x, 0.0, -lightDir.z)), 0.0);
    float dome = 0.2 + 0.8 * clamp(rd.y, 0.0, 1.0);
    float sol = 0.2 + 0.8 * clamp(-rd.y, 0.0, 1.0);
    float ao = ambientOcclusion(pos, n);
    float shadowIntensity = getShadowIntensity(pos + n * 0.001, lightDir);
    
    vec3 light = vec3(0.0);
    light += ambient * vec3(0.3, 0.3, 0.3) * 1.1 * shadowIntensity;
    light += diffuse * vec3(0.65, 0.6, 0.4) * ao * 1.2 * shadowIntensity;
    light += back * vec3(0.3, 0.2, 0.1) * 1.1;
    light += dome * vec3(0.5, 0.5, 0.7) * 1.1 * shadowIntensity;
    light += sol * vec3(0.3, 0.3, 0.32) * 1.1;
    
    vec3 objColor = getColor(info.y, pos) * 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.

References

  1. Fruxis, iq, Shadertoy