Cycles: Add new Sky Texture method including direct sunlight

This commit adds a new model to the Sky Texture node, which is based on a
method by Nishita et al. and works by basically simulating volumetric
scattering in the atmosphere.

By making some approximations (such as only considering single scattering),
we get a fairly simple and fast simulation code that takes into account
Rayleigh and Mie scattering as well as Ozone absorption.

This code is used to precompute a 512x128 texture which is then looked up
during render time, and is fast enough to allow real-time tweaking in the
viewport.

Due to the nature of the simulation, it exposes several parameters that
allow for lots of flexibility in choosing the look and matching real-world
conditions (such as Air/Dust/Ozone density and altitude).

Additionally, the same volumetric approach can be used to compute absorption
of the direct sunlight, so the model also supports adding direct sunlight.
This makes it significantly easier to set up Sun+Sky illumination where
the direction, intensity and color of the sun actually matches the sky.

In order to support properly sampling the direct sun component, the commit
also adds logic for sampling a specific area to the kernel light sampling
code. This is combined with portal and background map sampling using MIS.

This sampling logic works for the common case of having one Sky texture
going into the Background shader, but if a custom input to the Vector
node is used or if there are multiple Sky textures, it falls back to using
only background map sampling (while automatically setting the resolution to
4096x2048 if auto resolution is used).

More infos and preview can be found here:
https://docs.google.com/document/d/1gQta0ygFWXTrl5Pmvl_nZRgUw0mWg0FJeRuNKS36m08/view

Underlying model, implementation and documentation by Marco (@nacioss).
Improvements, cleanup and sun sampling by @lukasstockner.

Differential Revision: https://developer.blender.org/D7896

2 files changed, 208 insertions, 95 deletions

diff --git a/intern/cycles/kernel/svm/svm_sky.h b/intern/cycles/kernel/svm/svm_sky.h
index 50fe0c8232f..e877bd9a5c8 100644
--- a/intern/cycles/kernel/svm/svm_sky.h
+++ b/intern/cycles/kernel/svm/svm_sky.h
@@ -37,16 +37,16 @@ ccl_device float sky_perez_function(float *lam, float theta, float gamma)
          (1.0f + lam[2] * expf(lam[3] * gamma) + lam[4] * cgamma * cgamma);
 }
 
-ccl_device float3 sky_radiance_old(KernelGlobals *kg,
-                                   float3 dir,
-                                   float sunphi,
-                                   float suntheta,
-                                   float radiance_x,
-                                   float radiance_y,
-                                   float radiance_z,
-                                   float *config_x,
-                                   float *config_y,
-                                   float *config_z)
+ccl_device float3 sky_radiance_preetham(KernelGlobals *kg,
+                                        float3 dir,
+                                        float sunphi,
+                                        float suntheta,
+                                        float radiance_x,
+                                        float radiance_y,
+                                        float radiance_z,
+                                        float *config_x,
+                                        float *config_y,
+                                        float *config_z)
 {
   /* convert vector to spherical coordinates */
   float2 spherical = direction_to_spherical(dir);
@@ -90,16 +90,16 @@ ccl_device float sky_radiance_internal(float *configuration, float theta, float
           configuration[6] * mieM + configuration[7] * zenith);
 }
 
-ccl_device float3 sky_radiance_new(KernelGlobals *kg,
-                                   float3 dir,
-                                   float sunphi,
-                                   float suntheta,
-                                   float radiance_x,
-                                   float radiance_y,
-                                   float radiance_z,
-                                   float *config_x,
-                                   float *config_y,
-                                   float *config_z)
+ccl_device float3 sky_radiance_hosek(KernelGlobals *kg,
+                                     float3 dir,
+                                     float sunphi,
+                                     float suntheta,
+                                     float radiance_x,
+                                     float radiance_y,
+                                     float radiance_z,
+                                     float *config_x,
+                                     float *config_y,
+                                     float *config_z)
 {
   /* convert vector to spherical coordinates */
   float2 spherical = direction_to_spherical(dir);
@@ -121,93 +121,206 @@ ccl_device float3 sky_radiance_new(KernelGlobals *kg,
   return xyz_to_rgb(kg, make_float3(x, y, z)) * (M_2PI_F / 683);
 }
 
+/* Nishita improved sky model */
+ccl_device float3 geographical_to_direction(float lat, float lon)
+{
+  return make_float3(cos(lat) * cos(lon), cos(lat) * sin(lon), sin(lat));
+}
+
+ccl_device float3 sky_radiance_nishita(KernelGlobals *kg,
+                                       float3 dir,
+                                       float *nishita_data,
+                                       uint texture_id)
+{
+  /* definitions */
+  float sun_elevation = nishita_data[6];
+  float sun_rotation = nishita_data[7];
+  float angular_diameter = nishita_data[8];
+  bool sun_disc = (angular_diameter > 0.0f);
+  float3 xyz;
+  /* convert dir to spherical coordinates */
+  float2 direction = direction_to_spherical(dir);
+
+  /* render above the horizon */
+  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 half_angular = angular_diameter / 2.0f;
+    float dir_elevation = M_PI_2_F - direction.x;
+
+    /* if ray inside sun disc render it, otherwise render sky */
+    if (sun_disc && sun_dir_angle < half_angular) {
+      /* get 3 pixels data */
+      float3 pixel_bottom = make_float3(nishita_data[0], nishita_data[1], nishita_data[2]);
+      float3 pixel_top = make_float3(nishita_data[3], nishita_data[4], nishita_data[5]);
+      float y;
+
+      /* sun interpolation */
+      if (sun_elevation - half_angular > 0.0f) {
+        if (sun_elevation + half_angular > 0.0f) {
+          y = ((dir_elevation - sun_elevation) / angular_diameter) + 0.5f;
+          xyz = interp(pixel_bottom, pixel_top, y);
+        }
+      }
+      else {
+        if (sun_elevation + half_angular > 0.0f) {
+          y = dir_elevation / (sun_elevation + half_angular);
+          xyz = interp(pixel_bottom, pixel_top, y);
+        }
+      }
+      /* limb darkening, coefficient is 0.6f */
+      float limb_darkening = (1.0f -
+                              0.6f * (1.0f - sqrtf(1.0f - sqr(sun_dir_angle / half_angular))));
+      xyz *= limb_darkening;
+    }
+    /* sky */
+    else {
+      /* sky interpolation */
+      float x = (direction.y + M_PI_F + sun_rotation) / M_2PI_F;
+      float y = dir_elevation / M_PI_2_F;
+      if (x > 1.0f) {
+        x -= 1.0f;
+      }
+      xyz = float4_to_float3(kernel_tex_image_interp(kg, texture_id, x, y));
+    }
+  }
+  /* ground */
+  else {
+    if (dir.z < -0.4f) {
+      xyz = make_float3(0.0f, 0.0f, 0.0f);
+    }
+    else {
+      /* black ground fade */
+      float fade = 1.0f + dir.z * 2.5f;
+      fade = sqr(fade) * fade;
+      /* interpolation */
+      float x = (direction.y + M_PI_F + sun_rotation) / M_2PI_F;
+      if (x > 1.0f) {
+        x -= 1.0f;
+      }
+      xyz = float4_to_float3(kernel_tex_image_interp(kg, texture_id, x, -0.5)) * fade;
+    }
+  }
+
+  /* convert to rgb and adjust strength */
+  return xyz_to_rgb(kg, xyz) * 120000.0f;
+}
+
 ccl_device void svm_node_tex_sky(
     KernelGlobals *kg, ShaderData *sd, float *stack, uint4 node, int *offset)
 {
-  /* Define variables */
-  float sunphi, suntheta, radiance_x, radiance_y, radiance_z;
-  float config_x[9], config_y[9], config_z[9];
-
   /* Load data */
   uint dir_offset = node.y;
   uint out_offset = node.z;
   int sky_model = node.w;
 
-  float4 data = read_node_float(kg, offset);
-  sunphi = data.x;
-  suntheta = data.y;
-  radiance_x = data.z;
-  radiance_y = data.w;
-
-  data = read_node_float(kg, offset);
-  radiance_z = data.x;
-  config_x[0] = data.y;
-  config_x[1] = data.z;
-  config_x[2] = data.w;
-
-  data = read_node_float(kg, offset);
-  config_x[3] = data.x;
-  config_x[4] = data.y;
-  config_x[5] = data.z;
-  config_x[6] = data.w;
-
-  data = read_node_float(kg, offset);
-  config_x[7] = data.x;
-  config_x[8] = data.y;
-  config_y[0] = data.z;
-  config_y[1] = data.w;
-
-  data = read_node_float(kg, offset);
-  config_y[2] = data.x;
-  config_y[3] = data.y;
-  config_y[4] = data.z;
-  config_y[5] = data.w;
-
-  data = read_node_float(kg, offset);
-  config_y[6] = data.x;
-  config_y[7] = data.y;
-  config_y[8] = data.z;
-  config_z[0] = data.w;
-
-  data = read_node_float(kg, offset);
-  config_z[1] = data.x;
-  config_z[2] = data.y;
-  config_z[3] = data.z;
-  config_z[4] = data.w;
-
-  data = read_node_float(kg, offset);
-  config_z[5] = data.x;
-  config_z[6] = data.y;
-  config_z[7] = data.z;
-  config_z[8] = data.w;
-
   float3 dir = stack_load_float3(stack, dir_offset);
   float3 f;
 
-  /* Compute Sky */
-  if (sky_model == 0) {
-    f = sky_radiance_old(kg,
-                         dir,
-                         sunphi,
-                         suntheta,
-                         radiance_x,
-                         radiance_y,
-                         radiance_z,
-                         config_x,
-                         config_y,
-                         config_z);
+  /* Preetham and Hosek share the same data */
+  if (sky_model == 0 || sky_model == 1) {
+    /* Define variables */
+    float sunphi, suntheta, radiance_x, radiance_y, radiance_z;
+    float config_x[9], config_y[9], config_z[9];
+
+    float4 data = read_node_float(kg, offset);
+    sunphi = data.x;
+    suntheta = data.y;
+    radiance_x = data.z;
+    radiance_y = data.w;
+
+    data = read_node_float(kg, offset);
+    radiance_z = data.x;
+    config_x[0] = data.y;
+    config_x[1] = data.z;
+    config_x[2] = data.w;
+
+    data = read_node_float(kg, offset);
+    config_x[3] = data.x;
+    config_x[4] = data.y;
+    config_x[5] = data.z;
+    config_x[6] = data.w;
+
+    data = read_node_float(kg, offset);
+    config_x[7] = data.x;
+    config_x[8] = data.y;
+    config_y[0] = data.z;
+    config_y[1] = data.w;
+
+    data = read_node_float(kg, offset);
+    config_y[2] = data.x;
+    config_y[3] = data.y;
+    config_y[4] = data.z;
+    config_y[5] = data.w;
+
+    data = read_node_float(kg, offset);
+    config_y[6] = data.x;
+    config_y[7] = data.y;
+    config_y[8] = data.z;
+    config_z[0] = data.w;
+
+    data = read_node_float(kg, offset);
+    config_z[1] = data.x;
+    config_z[2] = data.y;
+    config_z[3] = data.z;
+    config_z[4] = data.w;
+
+    data = read_node_float(kg, offset);
+    config_z[5] = data.x;
+    config_z[6] = data.y;
+    config_z[7] = data.z;
+    config_z[8] = data.w;
+
+    /* Compute Sky */
+    if (sky_model == 0) {
+      f = sky_radiance_preetham(kg,
+                                dir,
+                                sunphi,
+                                suntheta,
+                                radiance_x,
+                                radiance_y,
+                                radiance_z,
+                                config_x,
+                                config_y,
+                                config_z);
+    }
+    else {
+      f = sky_radiance_hosek(kg,
+                             dir,
+                             sunphi,
+                             suntheta,
+                             radiance_x,
+                             radiance_y,
+                             radiance_z,
+                             config_x,
+                             config_y,
+                             config_z);
+    }
   }
+  /* Nishita */
   else {
-    f = sky_radiance_new(kg,
-                         dir,
-                         sunphi,
-                         suntheta,
-                         radiance_x,
-                         radiance_y,
-                         radiance_z,
-                         config_x,
-                         config_y,
-                         config_z);
+    /* Define variables */
+    float nishita_data[9];
+
+    float4 data = read_node_float(kg, offset);
+    nishita_data[0] = data.x;
+    nishita_data[1] = data.y;
+    nishita_data[2] = data.z;
+    nishita_data[3] = data.w;
+
+    data = read_node_float(kg, offset);
+    nishita_data[4] = data.x;
+    nishita_data[5] = data.y;
+    nishita_data[6] = data.z;
+    nishita_data[7] = data.w;
+
+    data = read_node_float(kg, offset);
+    nishita_data[8] = data.x;
+    uint texture_id = __float_as_uint(data.y);
+
+    /* Compute Sky */
+    f = sky_radiance_nishita(kg, dir, nishita_data, texture_id);
   }
 
   stack_store_float3(stack, out_offset, f);
diff --git a/intern/cycles/kernel/svm/svm_types.h b/intern/cycles/kernel/svm/svm_types.h
index e913d9e0489..f1ebb37e23e 100644
--- a/intern/cycles/kernel/svm/svm_types.h
+++ b/intern/cycles/kernel/svm/svm_types.h
@@ -414,7 +414,7 @@ typedef enum NodeWaveProfile {
   NODE_WAVE_PROFILE_TRI,
 } NodeWaveProfile;
 
-typedef enum NodeSkyType { NODE_SKY_OLD, NODE_SKY_NEW } NodeSkyType;
+typedef enum NodeSkyType { NODE_SKY_PREETHAM, NODE_SKY_HOSEK, NODE_SKY_NISHITA } NodeSkyType;
 
 typedef enum NodeGradientType {
   NODE_BLEND_LINEAR,