In developing virtual globes like Cesium, we need to be able to quickly determine when objects in the scene, like terrain tiles, satellites, buildings, vehicles, etc. are invisible and therefore do not need to be rendered. We do view-frustum culling, of course. But another important type of culling is horizon culling.

In the figure above, the green points are visible to the viewer. The red points are not visible because they are outside the view frustum, represented as heavy white lines. The blue point is inside the view frustum, but it is not visible to the viewer because it is occluded by the Earth. In other words, it is below the horizon. Horizon culling is the straightforward idea that you do not need to render objects that lie below the horizon as viewed from the current viewer position. As straightforward as this sounds, the details get tricky, especially because it needs to be very fast. Cesium will do this test hundreds of times per render frame in order to test the visibility of terrain tiles. It’s an important test, though. In the configuration in the figure above, terrain tiles covering the entire Earth lie inside the view frustum. Over half of them, however, are below the horizon and do not need to be rendered.

A few years back, Deron Ohlarik wrote two excellent articles on horizon culling. They are here and here. We’ve since developed an extension to his techniques that I would like to share here. While it only applies to static data like terrain tiles, we’ve found it to be very useful because it is both faster and more accurate than previous techniques. The accuracy improvement comes from horizon culling against an ellipsoidal model of the Earth rather than a spherical approximation.

I should mention up front that credit for this technique is entirely due to my colleague, Frank Stoner. My only contributions are implementing it in Cesium and writing about it here, after he did the hard work of deriving it.

Horizon culling a point against a sphere

As described by Ohlarik, for horizon-culling purposes, we can compute a bounding sphere for a static object like a terrain tile that is so tight that it is just a single point. If that point is below the horizon, then we can be sure that the entire tile is below the horizon, too. Our new technique is limited to culling a single point against an ellipsoid, so we start off assuming that this “occludee point” has already been computed. For details about how that is done, see the followup blog post.

I promised we would implement horizon culling against a general ellipsoid, and I’ll make good on that promise, but let’s start by using a simple unit sphere for horizon culling. Then, I’ll show that we can easily generalize this to an arbitrary ellipsoid. Consider the figure below:

In this figure, the blue circle is our unit sphere. The lines extending from the camera position and tangent to the sphere represent the horizon. The black, vertical line represents all of the horizon points. On our unit sphere, the horizon points lie on a plane and form a circle. The vectors from the camera position to all horizon points form an infinite cone.

The portion of the sphere and the space around it that is shaded grey represents the region “below” the horizon. Any point in the shaded region is not visible from the camera position. Intuitively, the point is below the horizon if it is inside the infinite cone formed by the tangent vectors, and it is behind the plane containing all the horizon points.

The plane test

First, let’s develop an inexpensive test to determine which side of the plane a point is on. Consider the figure below:

We know \(\vec{VC}\) and \(\vec{VT}\): they’re just the vectors from the viewer to the target point and to the center of the ellipsoid, respectively. We also know that \(\vec{HC}\) is a unit vector, because we’re working with a unit sphere. By the Pythagorean theorem:

Next, we note that triangles \(\triangle VCH\) and \(\triangle HCP\) are similar triangles. They share an angle at point \(C\) and each have another angle that is 90 degrees. So:

So the distance from the viewer to the plane is:

If the projection of \(\vec{VT}\) onto \(\vec{VC}\) is less than \(|\vec{VP}|\), the target point is in front of the plane. In other words, the target point is behind the horizon plane when:

We can simplify this by multiplying both sides by \(|\vec{VC}|\), yielding:

So there it is. To determine if the target point is behind the horizon plane, take the dot product of the vector from the viewer to the target with the vector from the viewer to the center of the ellipsoid. If that’s larger than the magnitude squared of the vector from the viewer to the center of the ellipsoid, minus one, then the target is behind the plane. No square roots or trigonometric functions required.

The cone test

If the target point is in front of the horizon plane, then the target point is definitely not occluded by the sphere and we’re done. If it’s behind the plane, however, the target point may or may not be occluded. If the target point is also inside the infinite cone formed by connecting the viewer to all of the horizon points, then it is occluded. If it is outside that cone, then it is not occluded. So how do we test the point against the cone?

Let’s look at our diagram again, this time with angle \(\angle HVC\) labelled \(\alpha\) and \(\angle TVC\) labelled \(\beta\).

We can see that point \(T\) is inside the cone when:

Or, for angles \(0 <= \theta <= \pi\):

Angle \(\alpha\) is part of right-triangle \(\triangle VCH\), so let’s rewrite the right-hand side of the inequality using the trig identity \(cos(\theta) = \frac{adjacent}{hypotenuse}\):

And let’s rewrite the left-hand side in terms of the dot product and multiply both sides by \(|\vec{VC}|\):

Now we’ll square both sides to eliminate the magnitudes, which otherwise involve a square root:

By squaring both sides, we’re effectively testing the target point against a double cone, where the second cone points from the viewer away from the ellipsoid. This does not impact our results, however, because a target point behind the viewer is guaranteed to be in front of the horizon plane. Any point in front of the horizon plane cannot be horizon culled, so this second cone test is not required.

Now where do we stand? \(\vec{VC}\) and \(\vec{VT}\) are easy to compute from our known ellipsoid center, target point, and viewer positions. \(\vec{VH}\) is not so obvious. But remember way back in the section on testing against the plane? We discovered that:

Not only is this easy to compute, but we already did so in the course of determining which side of the plane the point was on. Similarly, we already computed \(\vec{VT} \cdot \vec{VC}\).

So our final inequality, which requires only a few more arithmetic operations to evaluate, is:

If this inequality is true, then the target point is inside the cone. If it was also behind the horizon plane, then the target point is occluded.

Generalizing to an ellipsoid

This is all very elegant in our nice little unit-sphere world. How do we generalize it to an arbitrary ellipsoid?

The equation for our unit sphere is:

And the equation for an ellipsoid is:

where \(a\), \(b\), and \(c\) are the radii of the ellipsoid along the \(x\), \(y\), and \(z\) axes, respectively.

Given an ellipsoid centered at the origin, a viewer position, and a target position, we can apply a scaling transformation to all coordinates in order to create an equivalent problem in which the ellipsoid is actually a unit sphere. A matrix to perform the scaling operation looks like this:

We call this scaled coordinate system ellipsoid-scaled space, and have found it useful for solving a variety of problems on an ellipsoid.

A more rigorous treatment of this subject can be found in section 2 of the supplement to GPU Ray Casting of Virtual Globes, a poster presented at SIGGRAPH 2010.

Code

I thought it was important to write out all the math, but it all boils down to some simple code. Each time the camera position changes, we execute:

// Ellipsoid radii - WGS84 shown here
var rX = 6378137.0;
var rY = 6378137.0;
var rZ = 6356752.3142451793;

// Vector CV
var cvX = cameraPosition.x / rX;
var cvY = cameraPosition.y / rY;
var cvZ = cameraPosition.z / rZ;

var vhMagnitudeSquared = cvX * cvX + cvY * cvY + cvZ * cvZ - 1.0;

Then, for each point we wish to test for occlusion culling:

// Target position, transformed to scaled space
var tX = position.x / rX;
var tY = position.y / rY;
var tZ = position.z / rZ;

// Vector VT
var vtX = tX - cvX;
var vtY = tY - cvY;
var vtZ = tZ - cvZ;
var vtMagnitudeSquared = vtX * vtX + vtY * vtY + vtZ * vtZ;

// VT dot VC is the inverse of VT dot CV
var vtDotVc = -(vtX * cvX + vtY * cvY + vtZ * cvZ);

var isOccluded = vtDotVc > vhMagnitudeSquared &&
                 vtDotVc * vtDotVc / vtMagnitudeSquared > vhMagnitudeSquared;

In Cesium, we pre-compute the scaled-space position instead of doing it before each test as shown above.

Future

Using this technique for terrain culling in Cesium allowed us to avoid drawing about 15% of the tiles we otherwise would have drawn in common scenes, versus our previous technique of culling with a minimum-radius bounding sphere. Happily, the new test is faster to execute for each tile as well!

One detail that we’ve skirted around so far is how we’re generating the “occludee” test points from our terrain tiles and other static geometry. Currently, we’re computing each tile’s occludee point based on the (false, but conservative) assumption that occlusion will be performed using a sphere formed from the minimum radius of the ellipsoid. By using a more accurate computation for the occludee point, we should be able to cull more tiles.

Update: This is covered in more detail in a followup post.

However, while the ellipsoid is a convenient and reasonably accurate surface for horizon culling, we must always keep in mind that real terrain is often below the ellipsoid. If we improve the computation of the occludee point, we must take care that the more accurate, relative to the ellipsoid, horizon culling does not end up culling tiles that are actually still visible, relative to the real terrain. This is especially likely to be a concern when rendering underwater terrain.