(Updated on 14 April 2012: added a paragraph about step 2 and fixed the example code)

Let's say it's Easter and the weather sucks.
You'd like to be beside the seaside, watching a beautiful sunset but...
The wind is blowing, the sea is quite far and you are lazy.

If you are like me you start coding a digital version of it. I know, there's not much people like me :)

Jokes aside, I was thinking about terrain rendering through raymarching with distance fields (also known as sphere tracing). The nice thing about this technique is that you can render almost any function and shade it just like you would do with raytracing. In other words, it's fast and easy to code (at least if you stick to the basics) and it's fun to play with it.
The results can be absolutely stunning, as the famous 4KB intro Elevated taught us in 2009.

If your browser supports WebGL you can watch my "Sunset on the sea" experiment
and play with the code directly in your browser on GLSL Sandbox:
Sunset on the sea v.1.0.1 (original version here)



And here is how it works...

If you need an introduction on raymarching, the best site where you can go is the one by I. Quilez, where you can find a useful presentation on this technique and some ideas about terrain rendering.

Rendering a terrain (or the sea in this case) by using distance fields is, computationally wise, almost the worst-case scenario, as the primary rays have to travel near the surface, producing a very high number of iterations, as you can see in the image below.



I decided to avoid recurring to existing solutions and think by myself on a way to solve this problem.
After trying some different ideas I settled for the following algorithm, that you must apply to each ray casted from the camera:

1. Intersect the ray with a imaginary plane placed just above the sea (simple and fast raytracing technique).
2. Raymarch by the maximum between the distance function and a fixed step size.
3. As soon as you detect that you have stepped inside the sea bisect in this last interval until you are within a desired error threshold.

As you can see in the picture above, this algorithm is a lot faster and, if well tuned, can be quite accurate.

The purpose of step 2 is to find a rough interval where lies the first intersection between the ray and the sea. Note that a ray may intersect the sea in many different places (the front and the back of the same wave for example) but we are interested only in the nearest intersection.
In order to detect this first intersection we may use regular raymarching. This works well when a ray is distant from waves and when it's almost perpendicular to the surface, but becomes very slow when rays are almost parallel to the sea surface (as shown in the first graph).
Proceeding by fixed steps avoids this problem, but is less optimal than regular raymarching in some cases.
We can combine these two approches and get the best of each one. We can simply proceed by the maximum between the fixed step and the distance from the objects.
Additionally, I found that slightly increasing at each iteration the size of the fixed step works even better.

In step 3 of the algorithm I mentioned the bisection method. The idea is that you have two points: one over the surface and more near to the camera that we call "distMin" and one below the surface and more distant that we call "distMax".
Remember that determining if a point is over or below the surface in raymarching is as simple as determining if the distance function is greater or lower than zero.
The bisection algorithm that you would normally apply is something like this:

for (int i = 0; i < MAX_STEPS; i++) {
    distMid = distMin + (distMax - distMin) / 2.0;
    distFromObjects = map(ray_start + ray_dir * distMid);
    if (abs(distFromObjects) < PRECISION_LIMIT) return distMid; // Enough precision, stop
    if (distFromObjects > 0.0) {
        // distMid is outside the sea, then we can increase distMin
        distMin = distMid;
    } else {
        // distMid is inside the sea, then we can decrease distMax
        distMax = distMid;
    }
}
return distMid;

To this algorithm I applied a nice little optimization: as you are already computing the distance function for distMid at each step, you can use this information to further reduce the interval you are considering, then you can change the code above to the following one (differences in red):

for (int i = 0; i < MAX_STEPS; i++) {
    distMid = distMin + (distMax - distMin) / 2.0;
    distFromObjects = map(ray_start + ray_dir * distMid);
    if (abs(distFromObjects) < PRECISION_LIMIT) return distMid; // Enough precision, stop
    if (distFromObjects > 0.0) {
        // distMid is outside the sea, then we can increase distMin
        distMin = distMid + distFromObjects;
    } else {
        // distMid is inside the sea, then we can decrease distMax
        distMax = distMid + distFromObjects; // Note: distFromObjects is negative here
    }
}
return distMid;

A graphical representation of this idea (the numbers 1, 2 and 3 represent the order of the bisection steps):



Another little optimization i could think of is in the "gradientNormalFast" function.
In raymarching you generally retrieve normals by computing the gradient from six different points near the point you are considering. In this case I could use a less accurate version where I only compute three of these values.