Raymarching Water Body - Part I
<< Home

Raymarching Water Body - Part I

WARNING: WebGL2 required! So if you are not using a PC with a better-than-below-average GPU, stuffs below might not work.

Introduction

Water is beloved by all humankind (I think). And who doesn’t like water generated by pure math? Well, not me! So today, I am going to break down Seascape (finally not fruxis anymore!), and let’s see how that works out! Anyway first, a really basic scene:

vec3 getSkyColor(vec3 rd) {
    vec3 sky = vec3(0.58, 0.77, 1.0);
    return mix(sky, vec3(1.0), clamp(rd.y * 1.5, 0.0, 1.0));
}

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

vec2 map(vec3 p) {
    float id = -1.0;
    float closest = 1000.0;
    
    float dist = sol(p);
    if (dist < closest) { closest = dist; id = 0.5; }
    
    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.0001) {
            id = info.y;
            break;
        }
        depth += clamp(info.x, 0.01, 2.0);
    }
    return vec2(depth, id);
}

vec3 getColor(float id, vec3 p, vec3 rd) {
    if (id < -0.5) { return getSkyColor(rd); }
    if (id < 1.0) { return vec3(0.25, 0.55, 1.0); }
    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(0.0, 1.0, 5.0);
    vec3 center = vec3(0.0, 0.0, 0.0);
    
    vec3 lightDir = normalize(vec3(1.0, 1.0, 1.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);
    
    vec3 light = vec3(0.0);
    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(n.y, 0.0, 1.0);
    
    light += ambient * vec3(0.0, 0.0, 0.0);
    light += diffuse * vec3(1.2, 1.2, 1.2);
    light += back * vec3(0.52, 0.42, 0.34);
    light += dome * vec3(0.1, 0.1, 0.1);
    if (info.y < -0.5) {
        light = vec3(1.0);
    }

    vec3 objColor = getColor(info.y, pos, 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.

Nothing else beside a bright blue sky and an… azure… plane. Well, not to worry. It’s time to dive into WATER (pun intended)!

Implementation

First, let’s go back and take a good look at the seascape:

Seascape

What makes it look real? It’s the water details. It’s basically water bumping around, with nice lighting. And that’s why at first, we need to implement its basic water shape, which as we can see here, is clearly the work of fBm, but the variant is not known. Well as a matter of fact, this isn’t any variant we’ve see before all right; it’s the result of sine and cosine. First, let’s make wavey patterns in 1D. Take a look at this graph:

float absSin(vec2 uv);
float absCos(vec2 uv);
float amalgamation(vec2 uv);

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    xy *= 2.0;
    xy.x += time * 5.0;
    float dist = pow(
        max(1.0 - abs(xy.y - absSin(xy) * 1.5 + 0.5), 0.0),
        9.0);
    float v = pow(
        max(1.0 - abs(xy.y - absCos(xy) * 1.5 + 0.5), 0.0),
        9.0);
    float a = pow(
        max(1.0 - abs(xy.y - amalgamation(xy) * 1.5 + 0.5), 0.0),
        9.0);
    if (mod(xy, 0.2).x <= 0.1) {
        dist = 0.0;
        v = 0.0;
    }
    color = vec4(dist, v, a, 1.0);
}
// red line
float absSin(vec2 uv) {
    return 1.0 - abs(sin(uv.x));
}

// green line
float absCos(vec2 uv) {
    return abs(cos(uv.x));
}

// blue line
float amalgamation(vec2 uv) {
    return mix(absSin(uv), absCos(uv), absSin(uv));
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Look at the final amalgamation line. Sharp, but not that sharp. Smooth, but not that smooth. The best of both worlds. Powering it down a little will make it look a tad sharper:

float wave(vec2 uv);

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    xy *= 2.0;
    xy.x += time * 5.0;
    float dist = pow(
        max(1.0 - abs(xy.y - wave(xy) * 1.5 + 0.5), 0.0),
        9.0);
    color = vec4(dist, dist, dist, 1.0);
}
float wave(vec2 uv) {
    float a = 1.0 - abs(sin(uv.x));
    float b = abs(cos(uv.x));
    return pow(mix(a, b, a), 0.5);
}
Recompile | Click on the canvas to toggle running. It is paused by default.

And this is the wavey pattern we are using! Together with perlin noise, it will look a lot like natural, random waves.

float wave(vec2 uv);

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    xy *= 2.0;
    xy.x += time * 5.0;
    float dist = pow(
        max(1.0 - abs(xy.y - wave(xy) * 1.5 + 0.5), 0.0),
        9.0);
    color = vec4(dist, dist, dist, 1.0);
}
float rand(float i) {
    return fract(sin(i * 41245.123)) * 2.0 - 1.0;
}

float perlin(float v) {
    float u = floor(v);
    float f = smoothstep(0.0, 1.0, fract(v));
    return mix(rand(u), rand(u + 1.0), f);
}

float pattern(vec2 uv) {
    float choppiness = 2.0;
    float a = 1.0 - abs(sin(uv.x));
    float b = abs(cos(uv.x));
    return pow(1.0 - pow(mix(a, b, a), 0.5), choppiness);
}

float wave(vec2 uv) {
    uv.x += perlin(uv.x); // offset by itself
    return pattern(uv);
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Introducing wave choppiness: the choppier it is, the sharper those waves get. And finally, wrapping wave() in an fBm function to add details.

float rand(float i) {
    return fract(sin(i * 41245.123)) * 2.0 - 1.0;
}

float perlin(float v) {
    float u = floor(v);
    float f = smoothstep(0.0, 1.0, fract(v));
    return mix(rand(u), rand(u + 1.0), f);
}

float wave(vec2 uv, float choppiness) {
    uv.x += perlin(uv.x); // offset by itself
    float a = 1.0 - abs(sin(uv.x));
    float b = abs(cos(uv.x));
    return pow(1.0 - pow(mix(a, b, a), 0.5), choppiness);
}

float fbm(vec2 uv);

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    xy *= 2.0;
    float dist = pow(
        max(1.0 - abs(xy.y - fbm(xy) * 1.5 + 0.5), 0.0),
        9.0);
    color = vec4(dist, dist, dist, 1.0);
}
float fbm(vec2 uv) {
    float height = 0.6;
    float amplitude = 0.9;
    const int octaves = 5;
    float frequency = 0.16;
    float value = 0.0;
    float choppiness = 4.0;
    for (int i = 0; i < octaves; i++) {
        float d = wave((uv + time) * frequency, choppiness);
        uv *= 2.0;
        value += height * d;
        frequency *= 2.0;
        height *= 0.22;
        choppiness = mix(choppiness, 1.0, 0.2);
    }
    return value;
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Finally, we have this irregular moving, cool water shaped pattern. Changing the palette and you will see how closely this resembles water body:

float rand(float i) {
    return fract(sin(i * 41245.123)) * 2.0 - 1.0;
}

float perlin(float v) {
    float u = floor(v);
    float f = smoothstep(0.0, 1.0, fract(v));
    return mix(rand(u), rand(u + 1.0), f);
}

float wave(vec2 uv, float choppiness) {
    uv.x += perlin(uv.x); // offset by itself
    float a = 1.0 - abs(sin(uv.x));
    float b = abs(cos(uv.x));
    return pow(1.0 - pow(mix(a, b, a), 0.5), choppiness);
}

float fbm(vec2 uv) {
    float height = 0.6;
    float amplitude = 0.9;
    const int octaves = 5;
    float frequency = 0.16;
    float value = 0.0;
    float choppiness = 4.0;
    for (int i = 0; i < octaves; i++) {
        float d = wave((uv + time) * frequency, choppiness);
        uv *= 2.0;
        value += height * d;
        frequency *= 2.0;
        height *= 0.22;
        choppiness = mix(choppiness, 1.0, 0.2);
    }
    return value;
}
void main() {
    vec2 xy = uv * 2.0 - 1.0;
    xy *= 2.0;
    float dist = fbm(xy) * 1.5 - 0.5;
    
    if (xy.y - dist < 0.0) {
        // in water
        color = vec4(0.3, 0.52, 0.67, 1.0);
    } else {
        // sky
        color = vec4(0.6 + xy.y * 0.1, 0.8 + xy.y * 0.1, 1.0, 1.0);
    }
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Another palette, with a non-changin uv in fBm, will also make this a perfect desert:

float rand(float i) {
    return fract(sin(i * 41245.123)) * 2.0 - 1.0;
}

float perlin(float v) {
    float u = floor(v);
    float f = smoothstep(0.0, 1.0, fract(v));
    return mix(rand(u), rand(u + 1.0), f);
}

float wave(vec2 uv, float choppiness) {
    uv.x += perlin(uv.x); // offset by itself
    float a = 1.0 - abs(sin(uv.x));
    float b = abs(cos(uv.x));
    return pow(1.0 - pow(mix(a, b, a), 0.5), choppiness);
}

float fbm(vec2 uv) {
    float height = 0.6;
    float amplitude = 0.9;
    const int octaves = 5;
    float frequency = 0.16;
    float value = 0.0;
    float choppiness = 4.0;
    for (int i = 0; i < octaves; i++) {
        float d = wave((uv + time) * frequency, choppiness);
        value += height * d;
        frequency *= 2.0;
        height *= 0.22;
        choppiness = mix(choppiness, 1.0, 0.2);
    }
    return value;
}
void main() {
    vec2 xy = uv * 2.0 - 1.0;
    xy *= 2.0;
    float dist = fbm(xy) * 1.5 - 0.5;
    
    if (xy.y - dist < 0.0) {
        // in desert
        color = vec4(0.95, 0.92, 0.6, 1.0);
    } else {
        // sky
        color = vec4(0.92 + xy.y * 0.2, 0.8 + xy.y * 0.2, 0.3 + xy.y * 0.4, 1.0);
    }
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Now even though we have water (yay), we still have one problem (aww): the water is static, and moving in only one direction. In the ocean, waves in all sorts of direction. How to move the water around in this 1D world? Well, it’s really, really straightforward. We can do this by sampling another wave in the fbm function!

float rand(float i) {
    return fract(sin(i * 41245.123)) * 2.0 - 1.0;
}

float perlin(float v) {
    float u = floor(v);
    float f = smoothstep(0.0, 1.0, fract(v));
    return mix(rand(u), rand(u + 1.0), f);
}

float wave(vec2 uv, float choppiness) {
    uv.x += perlin(uv.x); // offset by itself
    float a = 1.0 - abs(sin(uv.x));
    float b = abs(cos(uv.x));
    return pow(1.0 - pow(mix(a, b, a), 0.5), choppiness);
}

float fbm(vec2 uv);

void main() {
    vec2 xy = uv * 2.0 - 1.0;
    xy *= 2.0;
    float dist = fbm(xy) * 1.5 - 0.5;
    
    if (xy.y - dist < 0.0) {
        // in water
        color = vec4(0.3, 0.52, 0.67, 1.0);
    } else {
        // sky
        color = vec4(0.6 + xy.y * 0.1, 0.8 + xy.y * 0.1, 1.0, 1.0);
    }
}
float fbm(vec2 uv) {
    float height = 0.6;
    float amplitude = 0.9;
    const int octaves = 5;
    float frequency = 0.16;
    float value = 0.0;
    float choppiness = 4.0;
    for (int i = 0; i < octaves; i++) {
        float d = wave((uv + time) * frequency, choppiness);
        d += wave((uv - time * 0.8) * frequency, choppiness);
        uv *= 2.0;
        value += height * d;
        frequency *= 2.0;
        height *= 0.22;
        choppiness = mix(choppiness, 1.0, 0.2);
    }
    return value;
}
Recompile | Click on the canvas to toggle running. It is paused by default.

Expanding into 2D

After knowing how water work in 1D, expanding it into 2D is incredibly easy:

  1. Upgrade our random function into a vec2 → float noise function.
  2. Upgrade our perlin noise function to a 2D perlin noise function.
  3. Changing the wave function’s product to the following:
float wave(vec2 uv, float choppiness);

vec3 getSkyColor(vec3 rd) {
    vec3 sky = vec3(0.58, 0.77, 1.0);
    return mix(sky, vec3(1.0), clamp(rd.y * 1.5, 0.0, 1.0));
}

// === WAVE CODE START === //
float rand2d(vec2 p) {
    return fract(sin(
        dot(p, vec2(12.345, 67.890))
    ) * 41234.45) * 2.0 - 1.0;
}

float perlin(vec2 p) {
    vec2 u = floor(p);
    vec2 f = fract(p);
    vec2 s = smoothstep(0.0, 1.0, f);
    
    float a = rand2d(u);
    float b = rand2d(u + vec2(1.0, 0.0));
    float c = rand2d(u + vec2(0.0, 1.0));
    float d = rand2d(u + vec2(1.0, 1.0));
    
    return mix(mix(a, b, s.x), mix(c, d, s.x), s.y);
}

float fbm(vec2 uv) {
    float height = 0.6;
    float amplitude = 0.9;
    const int octaves = 5;
    float frequency = 0.16;
    float value = 0.0;
    float choppiness = 4.0;
    for (int i = 0; i < octaves; i++) {
        float d = wave((uv + time) * frequency, choppiness);
        d += wave((uv - time * 0.8) * frequency, choppiness);
        uv *= 2.0;
        value += height * d;
        frequency *= 2.0;
        height *= 0.22;
        choppiness = mix(choppiness, 1.0, 0.2);
    }
    return value;
}
// === WAVE CODE END === //

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

vec2 map(vec3 p) {
    float id = -1.0;
    float closest = 1000.0;
    
    float dist = sol(p);
    if (dist < closest) { closest = dist; id = 0.5; }
    
    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.0001) {
            id = info.y;
            break;
        }
        depth += clamp(info.x, 0.01, 2.0);
    }
    return vec2(depth, id);
}

vec3 getColor(float id, vec3 p, vec3 rd) {
    if (id < -0.5) { return getSkyColor(rd); }
    if (id < 1.0) { return vec3(0.25, 0.55, 1.0); }
    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(0.0, 5.0, 5.0);
    vec3 center = vec3(0.0, 0.0, 0.0);
    
    vec3 lightDir = normalize(vec3(1.0, 1.0, 1.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);
    
    vec3 light = vec3(0.0);
    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(n.y, 0.0, 1.0);
    
    light += ambient * vec3(0.0, 0.0, 0.0);
    light += diffuse * vec3(1.2, 1.2, 1.2);
    light += back * vec3(0.52, 0.42, 0.34);
    light += dome * vec3(0.1, 0.1, 0.1);
    if (info.y < -0.5) {
        light = vec3(1.0);
    }

    vec3 objColor = getColor(info.y, pos, rd) * light;
    objColor = pow(objColor, vec3(0.4545));

    color = vec4(objColor, 1.0);
}
float wave(vec2 uv, float choppiness) {
    uv += perlin(uv); // offset by itself
    vec2 a = 1.0 - abs(sin(uv));
    vec2 b = abs(cos(uv));
    vec2 smoothed = mix(a, b, a);
    return pow(1.0 - pow(smoothed.x * smoothed.y, 0.5), choppiness);
}
Recompile | Click on the canvas to toggle running. It is paused by default.

By multiplying the product, we can made the 2 dimensional result back to one dimension again, which is good. The water’s quality could be further increased with a high-quality fbm, including rotation and shifting each octave:

vec3 getSkyColor(vec3 rd) {
    vec3 sky = vec3(0.58, 0.77, 1.0);
    return mix(sky, vec3(1.0), clamp(rd.y * 1.5, 0.0, 1.0));
}

// === WAVE CODE START === //
float rand2d(vec2 p) {
    return fract(sin(
        dot(p, vec2(12.345, 67.890))
    ) * 41234.45) * 2.0 - 1.0;
}

float perlin(vec2 p) {
    vec2 u = floor(p);
    vec2 f = fract(p);
    vec2 s = smoothstep(0.0, 1.0, f);
    
    float a = rand2d(u);
    float b = rand2d(u + vec2(1.0, 0.0));
    float c = rand2d(u + vec2(0.0, 1.0));
    float d = rand2d(u + vec2(1.0, 1.0));
    
    return mix(mix(a, b, s.x), mix(c, d, s.x), s.y);
}

float fbm(vec2 uv);

float wave(vec2 uv, float choppiness) {
    uv += perlin(uv); // offset by itself
    vec2 a = 1.0 - abs(sin(uv));
    vec2 b = abs(cos(uv));
    vec2 smoothed = mix(a, b, a);
    return pow(1.0 - pow(smoothed.x * smoothed.y, 0.5), choppiness);
}
// === WAVE CODE END === //

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

vec2 map(vec3 p) {
    float id = -1.0;
    float closest = 1000.0;
    
    float dist = sol(p);
    if (dist < closest) { closest = dist; id = 0.5; }
    
    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.0001) {
            id = info.y;
            break;
        }
        depth += clamp(info.x, 0.01, 2.0);
    }
    return vec2(depth, id);
}

vec3 getColor(float id, vec3 p, vec3 rd) {
    if (id < -0.5) { return getSkyColor(rd); }
    if (id < 1.0) { return vec3(0.25, 0.55, 1.0); }
    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(0.0, 5.0, 5.0);
    vec3 center = vec3(0.0, 0.0, 0.0);
    
    vec3 lightDir = normalize(vec3(1.0, 1.0, 1.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);
    
    vec3 light = vec3(0.0);
    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(n.y, 0.0, 1.0);
    
    light += ambient * vec3(0.0, 0.0, 0.0);
    light += diffuse * vec3(1.2, 1.2, 1.2);
    light += back * vec3(0.52, 0.42, 0.34);
    light += dome * vec3(0.1, 0.1, 0.1);
    if (info.y < -0.5) {
        light = vec3(1.0);
    }

    vec3 objColor = getColor(info.y, pos, rd) * light;
    objColor = pow(objColor, vec3(0.4545));

    color = vec4(objColor, 1.0);
}
mat2 rot2d(float deg) {
    float rad = radians(deg);
    float a = sin(rad), b = cos(rad);
    return mat2(b, a, -a, b);
}

float fbm(vec2 uv) {
    float height = 0.6;
    float amplitude = 0.9;
    const int octaves = 5;
    float frequency = 0.16;
    float value = 0.0;
    float choppiness = 4.0;
    vec2 shift = vec2(100.0, 0.0);
    for (int i = 0; i < octaves; i++) {
        float d = wave((uv + time) * frequency, choppiness);
        d += wave((uv - time * 0.8) * frequency, choppiness);
        uv = 2.0 * rot2d(45.0) * uv + shift;
        value += height * d;
        frequency *= 2.0;
        height *= 0.22;
        choppiness = mix(choppiness, 1.0, 0.2);
    }
    return value;
}
Recompile | Click on the canvas to toggle running. It is paused by default.

End of Part I

Well, that’s the end of part one. In this part, we only covered how to render the basic shape of 1D and 2D water. Nevertheless, here we are! (Really basic) Water with only basic lighting (not even specular!). And in next part, I am going to explore how to add more details to the water body. Stay tuned & cheerio!