Fix T78324: Different Sky Texture results between CPU and GPU

The problem here was numerical precision: The code calculates the angle between
sun and view direction, and the usual acos(dot(a, b)) approach for that has
poor numerical performance for almost parallel angles.

As a result, the generally tiny difference between floating point computation
between CPU and GPU was enough to make the sun vanish at different radii,
causing different results.

The new version fixes the difference by making the computation much more robust
on both platforms.

3 files changed, 17 insertions, 2 deletions

diff --git a/intern/cycles/kernel/shaders/node_sky_texture.osl b/intern/cycles/kernel/shaders/node_sky_texture.osl
index 08bc8f85120..20d379939ab 100644
--- a/intern/cycles/kernel/shaders/node_sky_texture.osl
+++ b/intern/cycles/kernel/shaders/node_sky_texture.osl
@@ -122,6 +122,11 @@ vector geographical_to_direction(float lat, float lon)
   return vector(cos(lat) * cos(lon), cos(lat) * sin(lon), sin(lat));
 }
 
+float precise_angle(vector a, vector b)
+{
+  return 2.0 * atan2(length(a - b), length(a + b));
+}
+
 color sky_radiance_nishita(vector dir, float nishita_data[9], string filename)
 {
   /* definitions */
@@ -138,7 +143,7 @@ color sky_radiance_nishita(vector dir, float nishita_data[9], string filename)
   if (dir[2] >= 0.0) {
     /* definitions */
     vector sun_dir = geographical_to_direction(sun_elevation, sun_rotation + M_PI_2);
-    float sun_dir_angle = acos(dot(dir, sun_dir));
+    float sun_dir_angle = precise_angle(dir, sun_dir);
     float half_angular = angular_diameter / 2.0;
     float dir_elevation = M_PI_2 - direction[0];
 
diff --git a/intern/cycles/kernel/svm/svm_sky.h b/intern/cycles/kernel/svm/svm_sky.h
index e877bd9a5c8..214c0cd1a9a 100644
--- a/intern/cycles/kernel/svm/svm_sky.h
+++ b/intern/cycles/kernel/svm/svm_sky.h
@@ -145,7 +145,7 @@ ccl_device float3 sky_radiance_nishita(KernelGlobals *kg,
   if (dir.z >= 0.0f) {
     /* definitions */
     float3 sun_dir = geographical_to_direction(sun_elevation, sun_rotation + M_PI_2_F);
-    float sun_dir_angle = acos(dot(dir, sun_dir));
+    float sun_dir_angle = precise_angle(dir, sun_dir);
     float half_angular = angular_diameter / 2.0f;
     float dir_elevation = M_PI_2_F - direction.x;
 
diff --git a/intern/cycles/util/util_math.h b/intern/cycles/util/util_math.h
index 737c834e073..8caabf6eac3 100644
--- a/intern/cycles/util/util_math.h
+++ b/intern/cycles/util/util_math.h
@@ -787,6 +787,16 @@ ccl_device_inline float compare_floats(float a, float b, float abs_diff, int ulp
   return (abs(__float_as_int(a) - __float_as_int(b)) < ulp_diff);
 }
 
+/* Calculate the angle between the two vectors a and b.
+ * The usual approach acos(dot(a, b)) has severe precision issues for small angles,
+ * which are avoided by this method.
+ * Based on "Mangled Angles" from https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf
+ */
+ccl_device_inline float precise_angle(float3 a, float3 b)
+{
+  return 2.0f * atan2f(len(a - b), len(a + b));
+}
+
 CCL_NAMESPACE_END
 
 #endif /* __UTIL_MATH_H__ */