Line-Sweep Ambient Occlusion
Generating ambient occlusion from heightmaps can be done with a line-sweep technique, where we sweep over a heightfield in several directions to gather an approximation of AO. In this post, we compute a sky visibility factor from each point to calculate the AO. Visibility is calculated using the iterative horizon-based approach described by Sean T. Barrett [1].
From a particular point on the heightmap, the horizon is the boundary between the sky and the heightmap as seen from that point. It represents the maximum elevation angle the horizon reaches in one azimuthal direction. For each azimuthal direction we sweep the heightfield, find the horizon for each point to calculate an occlusion factor and finally combine all the sweeps as final result. A sweep gives the contribution of sky visibility for a given direction.
Horizon
The first thing we need is getting the highest visible point for a heightmap location given a sweep direction. This horizon point is part of the convex hull of the heightmap. In this algorithm, it is situated before the sweep and is found by maximizing the $slope = \frac{\Delta y}{\Delta x}$ between a point and the horizon.
For a sweep, it is only necessary to look at the horizon located before the current iteration, as the sweep in the opposite direction will gather the AO from the horizon situated on the other side.
Line-Sweeps
The line sweeps over the heightmap can be done using the technique described by "A Fast Voxel Traversal Algorithm for Ray Tracing" [4] adjusted for 2D. This algorithm allows to iterate over the image in any given direction without missing pixels.
Once we have the horizon for a given point, we use the slope between this point and horizon as input value to an occlusion function. Choosing this function carefully is important and may result in a wide range of results, I recommend experimenting with it depending on your visual needs or artistic style.
In the demo above, $max(1 - e^{slope}, 0)$ is used as the occlusion function.
After getting all the single line-sweeps, we combine them together for the final result. Adding more sweeps generally gives better results but is more costly, rotating the sweep directions by 45 degrees increments is usually sufficient.
Falloff
In particular situations, the AO function used in the demo above might fail to give enough control over the final result, especially when the slope is too steep. To give more artistic control over this function, we can modulate the final AO value by a distance-based falloff. For example $e^{\frac{-x}{d}}$ or $\frac{d}{d + x^2}$ where $d$ is the distance point-horizon.
This can also be an optimization when tiling heightmaps, as it prevents having to look-up neighborhing tiles too far, since we know how far the horizon will stop contributing.
To expand control of depth over the heightmap features, it can be interesting to combine different falloff AO results. For example, combining a stronger narrow-scale with a lighter large-scale AO, to emphasize narrow terrain features while still retaining larger ones.
Limitation
The main limitation from this technique is the potential bands that may appear when the geometry is aligned with the line-sweep direction. This artifact may not be as common with terrain, but can be visible in heightmap of city landscapes or higher frequency heightmaps as shown below.
If banding is too much of an issue for the heightmap, it can be resolved by blurring the AO to soften the final result, or sweeping over more directions.
A simple example source code using this technique is available here.