Latitude precision issue in Globe Projection with Mali GPUs

[tavimori]

maplibre-gl-js version: 5.22.0

browser: any phone with Mali GPU?

Steps to Trigger Behavior

1. Mali GPU, zoom around 11, globe projection, latitude around 30.9 N.
2. Open the demo links on phones with Mali GPU.
3. The y-direction projection has a precision issue at zoom around 11 under the globe projection. Some distinct latitude values will overlap (see the image below); this applies to both custom layers and map tiles.

BTW, I noticed that there is a discussion on GPU "atan" error correction in https://github.com/maplibre/maplibre-gl-js/blob/main/developer-guides/globe.md, not sure if related.

Link to Demonstration

https://jsbin.com/zenatuk/edit?html,output
(modified based on the official example https://maplibre.org/maplibre-gl-js/docs/examples/add-a-custom-layer-with-tiles-to-a-globe/ )

Expected Behavior

(Screenshot on Xiaomi 13, Adreno 740)

(Screenshot on iPhone 16e, Apple A18)

Actual Behavior

(Screenshot on Huawei Nova 13, Mali-G610?)

We also tested on the Oppo Find N2 Flip (PGT110) with Mali-G710, which has the same issue.

[tavimori]

The root cause: the precision of sine and cosine.

I can confirm this is an issue with sin precision on certain devices.
A device test page has been created: https://jsbin.com/yamahed
Most GPUs handle sin and cos with highp precision correctly, which uses 32 bits (the Earth's surface-level resolution is just a few meters).

However, for Mali-G710 and Mali-G610, the error can vary by up to 200m. Which is significantly above the 32bits precision.

Potential solution

Having said that above, I also want to discuss a potential fix that can skip all trigonometric computation when handling latitudes. (and I noticed that the existing atan fix may also not be necessary anymore after the fix).

Currently, the GPU side trigonometric computation is mainly for converting the Mercator coordinates to Cartesian in 3D space.
The pipeline looks like this:

spherical_y = 2 * atan(exp(PI - mercator_y * 2*PI)) - PI/2   ... (2) inverse Mercator
  |
  v
sphere = ( sin(spherical_x) * cos(spherical_y),        ... (3) spherical -> Cartesian
           sin(spherical_y),
           cos(spherical_x) * cos(spherical_y) )

In this following, instead of calculating spherical_y as an intermediate step, we can calculate cos(spherical_y) and sin(spherical_y) directly, due to the property of Gudermannian.

The fix: algebraic identities that bypass sin/cos

The original code computes an intermediate angle and then calls sin/cos on it:

// Original (3 trig calls for Y: exp, atan, sin+cos)
float spherical_y = 2.0 * atan(exp(PI - mercator_y * 2.0 * PI)) - PI / 2.0;
float sin_sy = sin(spherical_y);
float cos_sy = cos(spherical_y);

We can eliminate atan, sin, and cos entirely using closed-form identities.

Derivation

Let t = exp(PI - mercator_y * 2 * PI). The original code computes:

spherical_y = 2 * atan(t) - PI/2

We need sin(spherical_y) and cos(spherical_y). Using the double-angle formulas for atan:

sin(2 * atan(t)) = 2t / (1 + t^2)
cos(2 * atan(t)) = (1 - t^2) / (1 + t^2)

These come from the substitution atan(t) = theta, so tan(theta) = t, which gives sin(theta) = t / sqrt(1 + t^2) and cos(theta) = 1 / sqrt(1 + t^2), then applying the double-angle identities sin(2*theta) and cos(2*theta).

Now applying the phase shift (- PI/2):

sin(alpha - PI/2) = -cos(alpha)
cos(alpha - PI/2) =  sin(alpha)

Where alpha = 2 * atan(t). Substituting:

sin(2*atan(t) - PI/2) = -cos(2*atan(t)) = -(1 - t^2) / (1 + t^2) = (t^2 - 1) / (t^2 + 1)
cos(2*atan(t) - PI/2) =  sin(2*atan(t)) =  2t / (t^2 + 1)

Result

// New (1 trig call: exp only, rest is arithmetic)
float t  = exp(PI - mercator_y * 2.0 * PI);
float t2 = t * t;
float denom = t2 + 1.0;
float sin_sy = (t2 - 1.0) / denom;    // exact equivalent of sin(spherical_y)
float cos_sy = (2.0 * t)  / denom;    // exact equivalent of cos(spherical_y)

This produces mathematically identical results to the original formula, but using only exp(), multiplications, additions, and divisions -- no atan, sin, or cos.

Performance impact

The fix also reduces computation cost:

	Original	New
Transcendental ops	exp + atan + sin + cos = 4	exp = 1
ALU ops (mul/add/div)	~3	~6 (mul, add, sub, 2x div)
Net	4 transcendental + 3 ALU	1 transcendental + 6 ALU

On GPUs, transcendental functions run on Special Function Units (SFUs) which are shared/limited resources and take multiple cycles each. Basic ALU operations are significantly cheaper. The net effect is a small performance improvement on all GPUs.

Lastly, there may also have GPU side sin call when dealing with longitude. These may still have issue.

[HarelM]

This is a very nice write up, thanks for looking into it!
I haven't verified the math, but I trust you did.
Feel free to open a PR, if it passes all the existing render tests it means it's an in-place replacement and is a good workaround to this issue.
It's worth filling issues for the GPU manufacturer if you know how to do it as I believe we will end up using trigonometric functions in the future and this workaround will not hold...

[kubapelc]

Very nice work!

If the remaining sin call is problematic, you can probably correct its error with similar periodic measurement trick as what is currently used for atan error, as long as the error is continuous and changes gradually.

[tavimori]

Very nice work!

If the remaining sin call is problematic, you can probably correct its error with similar periodic measurement trick as what is currently used for atan error, as long as the error is continuous and changes gradually.

Thank you for your reply. I noticed that you may be the author of the previous shader version, so I wanted to ask for your opinion as well. I currently believe the algorithmic changes mentioned above can fully resolve the precision issues with the atan function; therefore, I plan to implement this improvement and remove the specific atan precision correction (assuming this helps reduce the overall complexity of the project).

However, we do indeed still face a lingering issue. I have conducted some tests regarding the behavior of trigonometric functions—as well as functions like exp—within WebGL across various GPU architectures: https://webgl-test.lcblue.net/dashboard.
It appears that, currently, WebGL's trigonometric functions do suffer from certain issues; specifically on particular devices (primarily those featuring Intel UHD or ARM Mali GPUs), the precision error can be as large as approximately 300 meters. Furthermore, this precision issue is not consistent; at least on Mali GPUs, the sin function is not strictly monotonic—at the scale of the Earth's surface, the calculated values ​​can sometimes "reverse direction," deviating by several hundred meters.
This issue seems to have existed for many years: CesiumGS/cesium#2110.

Consequently, the "periodic measurement trick" you mentioned may not be sufficient to resolve this specific problem. However, given that this issue is not visually prominent—and considering that GPU technology is constantly evolving, which may lead to improvements over time—I suggest that a potentially better solution would be to trigger the transition to the Mercator projection at an earlier stage. For instance, to accommodate an error margin of around 300 meters, we could switch to the Mercator projection at a Z-level of approximately 8 (thereby, I assume, offloading the trigonometric calculations from the GPU to the CPU). Additionally, we could perhaps add a warning regarding GPU precision limitations on specific devices to the documentation for the vertical-perspective projection.

What does everyone think of this approach?

[HarelM]

If the current code fixes the current issue and simplifies the problem let's merge that.
If we have additional issues, let's open an issue and discuss the possible solutions and how they affect and if it's worth pursuing in terms of edge cases etc.

[kubapelc]

Go for it, I think your method of computing spherical coordinates from mercator is a big improvement.

As for the sin error: The atan error is also not monotonic and can also be in the range of hundreds of meters, but it changes very gradually. So the current correction cannot fix the error at the scale of the planet, but it is not noticeable at those scales anyway. Once the map is zoomed in enough for the error to be relevant, the entire screen usually covers a small enough region of the planet for the error to be roughly the same everywhere the user can see, so a single correction factor is good enough. Maybe there is still hope that this approach could work for sin as well?

I found a plot of the atan error once on some blog post, but it has since been deleted, so I reposting it here.

(note: the image seems to have black text on transparent background, so it is hard to see the label on Github'S dark mode)