Materials for Polygons on Terrain and the Evils of Inverse Trigonometry

One of the coolest details behind this feature is how texture coordinates are computed for polygons on terrain. Rather than create triangles that exactly match the underlying terrain geometry, Cesium implements polygons on terrain as image-space decals. One major benefit of this approach is that Cesium doesn’t have to rebuild geometry as the terrain level-of-detail changes, but the crucial difference for the purpose of this blog post is that we can’t use per-vertex texture coordinates.

We decided to compute texture coordinates per-fragment using the terrain depth texture and some per-geometry reference points. The basic idea is that you compute a terrain position for each fragment in the decal using the window coordinates and the terrain depth, then compare that position to some kind of “bounding box” to get a classic ([0, 1], [0, 1]) coordinate for texturing.

This is pretty straightforward in Columbus View, and also for polygons so small that they might as well be on a flat Earth - our “bounding box” is represented by a set of planes, so for each fragment we can compute the underlying terrain position’s distance to the planes.

But what about for polygons that span larger areas of the globe, or even cover more than 180 degrees of longitude? One possibility is to compute a spherical coordinate for each terrain position in the decal. This spherical coordinate is kind of like an approximated substitute for a latitude/longitude coordinate.

Note that to help support the analogy, we’re talking about altitude-azimuthal spherical coordinates instead of polar-azimuthal. Given a vector pointing from the center of the Earth to a position on terrain, the most straightfoward way to compute altitude and azimuth is:

The “bounding box” here then is similar to a cartographic Rectangle, with north, south, east, and west as spherical altitudes and azimuths.

Given the altitude and azimuth of a point on terrain inside the “bounding box,” you can compute a normalized texture coordinate with something like:

Elegant! Simple!

But then you go ahead and start proving out the concept by coloring a Rectangle on Cesium with the computed texture coordinates, red/green for u/v in [0, 1] range, blue for out-of-bounds… and you naively use glsl’s built-in inverse trigonometry.

What happened here was that we were computing the spherical rectangle using CPU inverse trigonometry and the per-fragment spherical coordinates using GPU-native inverse trig:

To my horror, it turn out my particular instance of WebGL’sasinandatan2were giving back significantly different values than what was being computed CPU side! In more detail, the altitude/azimuth computed at a corner of the rectangle on the CPU wasn’t the same as the value that the GPU was computing for the same spot, so the rectangle used for getting normalized texture coordinates wasn’t lining up. Here’s what the CPU thought the rectangle was in orange and what the GPU thought the rectangle was in pink. This wasn’t a constant offset either, it would shift depending on where the Rectangle Primitive was on the globe.

Using native asin also caused problems with polygons crossing the equator on my system. Even more disturbing, I only noticed this discontinuity on some platforms.

Here’s a demo in Sandcastle that reproduces the problem for me in Chrome on Windows 10 (Nvidia GT 750M) and on Android (Adreno 308) but not on Linux (Bay Trail Intel Graphics). However, on all platforms I still see the CPU/GPU alignment problem - click around your favorite places on the globe to move the rectangle, and marvel as the discrepancy changes. On some platforms like Intel Graphics on Linux the discrepancy might be less noticeable unless the rectangle covers a smaller area, so feel free to tweak the rectangleWidth variable.

I had already read in multiple places that inverse trigonometry in shaders is usually bad for performance, but discovering it also isn’t reliable and can differ across platforms was like learning that for some weird reason 2 + 3 = 5.1289 in Idaho and Nebraska. Built-in functions are arbitrary! Math can’t be trusted! Programming is anarchy!

(While we’re talking about anarchy and inverse trig, I found multiple articles stating that use of inverse trig in programming leads to ritual sacrifice of kittens. They can’t all be wrong right?!)

After I calmed down I started looking into fitting curves as inverse trig approximations. The basic idea is that with enough fiddling of constants and terms you can sometimes make an inexpensive function look like a more expensive one over the input range that you need. Sébastien Lagarde has an excellent blog post outlining how and why this should be done on the Playstation 4, including accuracy measures and tradeoffs of various approximations when applied to lighting computation. He even has GPU assembly instruction counts to prove a performance improvement!

Although better performance is always a nice plus, we looked to approximation more consistency between the CPU and GPU, as well as across platforms. This is almost the opposite of Lagarde’s problem - console developers should at least be able to count on their users having the same chips, drivers, and inverse trig functions under the hood, even if that inverse trig is slow. And so, Desmos plotter in hand, I embarked to find the perfect curve.

Unfortunately, the best approximations I found for asin still had discontinuity problems at the equator, leading me to suspect that my system’s native asin is a similar approximation buried somewhere in the graphics stack. The approximations looked spectacularly accurate when plotted out, but unfortunately I think the use of a square root and mirroring to cover the [-1, 1] range causes some kind of higher-degree discontinuity at input0.0, hence the pinching at the equator. There was actually a brief time when I thought a simple cubic function was the best curve for this case since there wasn’t a discontinuity, although you can probably tell from just the graph that this is pretty terrible for wide bands of latitude:

And here’s confirmation, using a texture on a ground primitive as a mask for continents and water. Materials on terrain polygons really shouldn’t be used for precisely mapping image features to latitude/longitude coordinates, but that was even more true when this was our best shot at an asin approximation:

Eventually I realized we can useatan2to compute the spherical altitude/latitude approximation as well - we just needed the magnitude of the vector when projected into the azimuthal plane, along with the vector’s vertical component.