We wish to compute the amount of light hitting a piece of ground terrain from the sky. We can compute this by measuring the vertical angle to the horizon in a 360-degree angle around each point. We can sample this horizon in multiple directions independently, reducing this to a problem of only computing the horizon in a single direction. (A more detailed sketch of this is given later.) In the past a common recommendation was to cast short 2D rays in each direction and measure the horizon along each ray. Using short rays reduced the high cost of raycasting over the entire surface, although it causes large features to fail to cast a sufficiently large skylight shadow. Instead we suggest a solution that makes sure features cast shadows over the entire surface. Rather than view the problem as raycasting, we consider the problem of computing the horizon along a single direction, through a single series of points along a line in the same direction.
Given a single-valued function characterized as a sequence of connected line-segments, we wish to compute for each segment centerpoint or endpoint the arc above it which is not shadowed by other parts of the function. Here we will solve for endpoints; centerpoints can be reasonably approximated using data from both adjacent endpoints.
As shown in the second illustration, we can separate this into left and right arc regions, since directly above a point must always be visible. We compute the left and right separately. Without loss of generality, we describe only computing the left amount.
To compute the amount on the left we find a line below which is the function and above which is open area; that is, we want to compute a value which is the slope of a line which is "tangent" to the line segments with lower x-values. Formally, the line must pass through at least one point on the function, and all points on the function must be on or below the line. Restricting this to series of connected line segments means that the line must pass through one line segment vertex, and all line segment vertices must be below or on the line.
In practice, we are given an array of numbers; we interpret the numbers as y-coordinates and use the array indices as x-coordinates. This produces a series of vertices connected by line segments in order. (To control effects at the edges we may have to introduce an extra point which will cast an initial shadow.) Rather than explicitly computing a slope, we will just compute a horizon vertex, since the horizon slope line will always pass through a vertex; the slope can be computed from the vertex.
Ok, I dug up this: Fast Horizon Computation at All Points of a Terrain with Visibility and Shading Applications (1998). I didn't read the details, but it's O(N log N), but computes many directions at once, and avoids the undersampling artifacts the algorithm described here will have. Somebody should compare them.
Edit Feb 2013: Found out that someone else published my algorithm in 2010: Scalable Height-Field Self-Shadowing.
One closely-related and well-known idea is computing shadows on 1D height fields from an infinitely-far-away light source. There is a simple algorithm which sequentially propogates a current shadow height, lowering it at each step until it hits the ground at which point the shadow is reset to the ground height. I also don't have any references for that.
At each step, we determine the horizon vertex/slope based on the candidate previous horizon vertices, then update the list of candidate horizon vertices to include the new point.
A naive version of the above algorithm uses all of the vertices as candidate horizon vertices, and requires O(n2) time. We can improve the performance by keeping a smaller list. To do so, we need to determine which vertices can shadow the vertices to their right.
Given a sequence of two line segments (three vertices), if the middle vertex is below the line formed by the left and right vertices, it cannot be the horizon for any vertices further to the right.
Generalizing this, it turns out here can definitely be more than two points that can be the horizon points for vertices to the right. But the only vertices that can contribute to the left horizon of a given vertex are the points that lie on the convex hull of the left line segments and a "baseline" below the heightfield (just to establish which side is a concavity). In the following images, vertices that could be a horizon depending on the location of the purple dot are marked in green; those that cannot have any effect are marked in red.
If a given vertex becomes ineligible to be a horizon vertex for a vertex at x=k, it is also ineligble to be a horizon vertex for all x > k. Thus the algorithm, finding horizons for increasing x, becomes: determine the horizon vertex based on the current convex hull vertices, then update the convex hull by deleting some existing candidates and possibly adding the new vertex.
To understand how we can efficiently determine the horizon and update the convex hull, we can return to the definition of the slope line. As noted previously, the horizon slope line is defined by the new vertex and the old horizon vertex. All the other vertices are below the horizon slope line.
The old convex hull vertices that are right of the horizon vertex are marked in red, and fit the previous characterization of vertices below-a-line-between-surrounding-vertices which means the vertices can't cast any shadow further to the right. The old convex hull vertices that are left of the horizon vertex, marked in green, can still potentially shadow vertices further to the right.
Therefore we can classify the existing candidate vertices at each step into two sets: those to the left of the horizon vertex are kept as candidates for later points (as is the horizon vertex itself); those to the right of the horizon vertex are discarded.
Keep a stack of candidate points (which define a convex hull), initialized with a dummy value.
Iterate sequentially through the (x,y) points.
For each point in the list, if the stack has only one point on it, use that point as the horizon point. Otherwise, compare the slope formed by the current point to the top point on the stack, and the slope from the current point to the second-topmost point on the stack. If the slope to the second point is smaller, then pop the stack and repeat; otherwise use the topmost point on the stack as the horizon for the current point, push the current point on the stack, and advance to the next point.
This algorithm takes O(N) time to process N vertices. Each vertex is iterated once, pushed once, and popped at most once. Since the inner loop always pops a vertex or terminates, it can only iterate at most N times, thus at most 2N slope comparisons are performed.
Given an input array H[0..N-1] of N float point values
Output an array of V[1..N-1] of indices of horizon vertices
Using temp array S[0..N-1] to store a temporary stack
Set S[0]=0, k=0, i=1
while i < N do
while k > 0 AND (H[i]-H[S[k]])/(i-S[k]) > (H[i]-H[S[k-1]])/(i-S[k-1])
k -= 1
endwhile
V[i] = S[k]
k += 1
S[k] = i
i += 1
endwhile
We can compute the negative x-axis horizon on a single row of a 2D height field using the above algorithm. Applying the algorithm right-to-left will compute the positive x-axis horizon. Applying the algorithm with a stride to step along the y-axis will allow computing the negative and positive y-axis horizons.
Stepping through the heightfield diagonally allows computing the "exact" skylight horizon along all four diagonals (note that because the samples are more widely spaced, the conversion from horizon point to horizon arc must account for the sample spacing, but if the calculation uses the actual (x,y) coordinates for the computation it will be correct).
Stepping by alternate amounts of (1,0) and (1,1) allows computing an approximation to the horizon along fractional lines. Alternatively, the algorithm can step along knight's move paths, stepping by (2,1) at a time. (This will cause it to miss certain fine details, however.)
Without even going to totally arbitrary bresenhm/DDA paths, the above techniques allow computing 16 horizon sample points by using 16 passes.
(Another use of the horizon value is for precomputing shadows for a moving sun, especially if the sun moves through an overhead position so all light is always cast along positive or negative x. However, the stored values don't help compute sun shadows on other objects in the scene, and so this is more plausibly solved using other shadow methods anyway, so it's not important in practice.)