Are you sick of horizon culling yet? Great, me neither!

Last time, we explained what horizon culling is all about, and showed a very efficient way to test a point for occlusion by an ellipsoid. The objects we want to test for occlusion are rarely simple points, however. In particular, we’d like to be able to test terrain tiles for occlusion by the ellipsoid. But terrain tiles are complicated objects composed of thousands of vertices.

Deron Ohlarik addressed this in a previous blog post by explaining that, for any arbitrary geometry, we can compute the position of a point, which we call the horizon occlusion point, with a special relationship to the geometry. No matter what direction a viewer approaches the geometry from, the point will become visible to the viewer at the same time or before any part of the geometry becomes visible. This is exactly what we need! Many details of how to compute such a point were left as an exercise for the reader, however. Furthermore, it was unclear if this approach could be generalized to an ellipsoid rather than a sphere. This blog aims to fill in both of these gaps.

Once again, credit for the technique presented here is due entirely to Frank Stoner.

Let’s take a look at our situation. As before, we transform all of our coordinates to the ellipsoid-scaled space by multiplying each component, X, Y, and Z, by the inverse of the ellipsoid’s radii along that axis.

In this figure, the Earth is shown in blue and a terrain tile is shown in brown. In the scaled-space, the Earth is a unit sphere. The center of the bounding sphere surrounding the terrain tile is shown as point \(C\). The bounding sphere is not a sphere in the scaled space, but that does not concern us because we are only going to use its center.

To start, we make the arbitrary decision that our horizon occlusion point will lie somewhere along this center line, \(\vec{OC}\), the vector from the center of the Earth to the center of the terrain tile’s bounding sphere. We only need to compute its distance along that vector.

Point \(V\) is a vertex in the terrain tile. Point \(H\) is a point on the horizon from the perspective of \(V\). There are an infinite number of horizon points from the perspective of \(V\), forming a circle on the unit sphere, but only two of these horizon points form a vector through \(V\) that intersects with the center line. One is shown as a solid line, \(\vec{HP}\). The other is shown connected to \(V\) as a dashed line. On the dashed line, the intersection with the center line occurs before point \(V\), so it will be closer to the center of the ellipsoid than the other intersection point and we do not need to be concerned with it. If point \(V\) were the only vertex in the terrain tile, then point \(P\) in this diagram would be our horizon occlusion point. With multiple vertices, we repeat the computation of \(P\) for each vertex, and then choose the one that is farthest from the ellipsoid.

So how do we compute point \(P\) for a given terrain tile vertex? Let’s label the various angles in the diagram.

After labeling angles \(\alpha\) and \(\beta\), we can express the other angles in terms of them, simply using the knowledge that the angles in a triangle add up to 180 degrees.

Then, by the law of sines:

Of course, \(sin(90 + \theta) = cos(\theta)\), so we can rewrite that as:

\(\beta\) is an angle in a right-triangle, so:

Then, we apply the compound angle formula:

We already know how to compute \(cos(\beta)\). We can compute \(sin(\beta)\) in the same way after applying the Pythagorean theorem:

We can rewrite \(cos(\alpha)\) in terms of the dot product:

Finally, we can express \(sin(\alpha)\) in terms of the magnitude of the cross product:

We now have everything we need to compute the magnitude of \(\vec{OP}\). Let’s pull it together:

Remember that we know \(\hat{OP}\) by construction; we chose it to point from the center of the ellipsoid to the center of the terrain tile’s bounding sphere. To compute the position of point \(P\) in the ellipsoid-scaled space, we simply multiply the direction by the magnitude computed above. Since our occlusion test works with a point expressed in the scaled space, we’re done. If we also want to know the position in real, unscaled coordinates, we only need to multiply each component of the position by the ellipsoid radius along the corresponding axis.

Here’s what the code looks like in Cesium, tweaked slightly for clarity:

```
function computeMagnitude(ellipsoid, position, scaledSpaceDirectionToPoint) {
var scaledSpacePosition = ellipsoid.transformPositionToScaledSpace(position);
var magnitudeSquared = scaledSpacePosition.magnitudeSquared();
var magnitude = Math.sqrt(magnitudeSquared);
var direction = scaledSpacePosition.divideByScalar(magnitude);
// For the purpose of this computation, points below the ellipsoid
// are considered to be on it instead.
magnitudeSquared = Math.max(1.0, magnitudeSquared);
magnitude = Math.max(1.0, magnitude);
var cosAlpha = direction.dot(scaledSpaceDirectionToPoint);
var sinAlpha = direction.cross(scaledSpaceDirectionToPoint).magnitude();
var cosBeta = 1.0 / magnitude;
var sinBeta = Math.sqrt(magnitudeSquared - 1.0) * cosBeta;
return 1.0 / (cosAlpha * cosBeta - sinAlpha * sinBeta);
}
```

As you can see, this computation is more expensive than the computation we described last time to test this point against the horizon. It’s probably possible to optimize it by testing each vertex using the cone test described previously and only computing the precise horizon occlusion point for the vertex if the vertex is found to be outside the cone. I’ll leave that as an exercise for the reader.

In any case, the cost of this computation is the main reason that it is primarily only appropriate for static geometry. If the geometry changes with respect to the ellipsoid, then this computation will need to be repeated on each change. That’s likely to get expensive.

Also, please keep in mind an important caveat when using this approach. In the real world, objects occluded by the WGS84 ellipsoid are not necessarily occluded by the real surface of the Earth. This is because the Earth’s surface is actually slightly below the ellipsoid in parts of the world. Depending on your application, using WGS84 as the occluding volume may be acceptable, or you may need to use a more conservative ellipsoid.