Materials for Polygons on Terrain and the Evils of Inverse Trigonometry

by

Cesium 1.46 added support for materials on terrain polygons, including textured materials for notational patterns or decals.

Materials

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.

Columbus View

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 a approximated substitute for a latitude/longitude coordinate.

Spherical coordinate system

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:

float altitude = asin(vector.z);
float azimuth = atan2(vector.y, vector.z);

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

Spherical Rectangle

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

floar u = (altitude - south) / (north - south);
floar v = (aziumth - west) / (east - west);

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.

Misalignment

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:

vec3 sphereNormal = normalize(worldCoord);
float latitude = asin(sphereNormal.z);
float longitude = atan2(sphereNormal.y, sphereNormal.x);

To my horror, it turn out my particular instance of WebGL’s asin and atan2 were 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.

Misalignment

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.

Platform discrepancies Left: Equatorial Discontinuity on Chrome, Windows 10, Nvidia. Right: Chrome, Linux, Intel (smaller area to show discrepancy)

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.

Desmos

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 input 0.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:

Desmos

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:

Mask Texture Left: reference use of a water mask texture. Right: using a cubic curve. Note: This use of materials for polygons on terrain is NOT high precision.

Eventually I realized we can use atan2 to 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.

float altitude = approximateAtan2(sqrt(vector.x * vector.x + vector.y * vector.y), vector.z);

atan and atan2 approximations often use quadratic functions, which have much better continuity and are possibly also faster on some GPUs than square roots. Here’s that Sandcastle toy modified to use our atan2 approximation instead. For the curious, there are a few more details on how we built our atan2 approximation right in the code documentation here.

So some takeaways:

  • If consistency is important, do not trust native GPU inverse trigonometry
  • When you have to use inverse trig, don’t get stuck with one method like I did with asin and spherical altitude/latitude
    • using atan2 was much better
    • insert joke about more than one way to skin a cat
    • kidding, I adore cats so unnecessary inverse trig makes me sad :(
  • Do not use materials for polygons on terrain to precisely map imagery features to terrain
    • spherical coordinates are not the same as latitude/longitude
    • for many x, approximateAtan2(x) !== atan2(x)!
    • use SingleTileImageryProvider for that

We put a lot of work into this feature, and we’re excited to see what users will build with it. In the meantime, we’re not fully done either - expect materials for polygons on 3D Tilesets soon!