EEVEE: Implement the missing Sky texture

I'm not sure if the Sky was deliberately left out or was just waiting for a
better moment, but so many I was disappointed that Sky in EEVEE is
completely white.

There are already 2 implementations (osl and gpu) so this is the third one.
Looking at other cases it seems that we are not supposed to share sources
between cycles and the rest? So the new util_sky_model files are just
copies of what is already in cycles, except that the data file uses the RGB
variant of the Hosek/Wilkie model, because we output RGB anyway (but can be
easily changed to XYZ if desired - the results are nearly identical).
I am not sure if it is okay to pass 3*9 float values as 3 mat4 uniforms (I
wanted to use mat3 but it does not work).
Also, should I cache the sky model data between renders if the parameters
do not change?

Reviewed By: fclem, brecht

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

17 files changed, 557 insertions, 73 deletions

diff --git a/intern/CMakeLists.txt b/intern/CMakeLists.txt
index fa18f4d793a..0758567bb78 100644
--- a/intern/CMakeLists.txt
+++ b/intern/CMakeLists.txt
@@ -30,6 +30,7 @@ add_subdirectory(opensubdiv)
 add_subdirectory(mikktspace)
 add_subdirectory(glew-mx)
 add_subdirectory(eigen)
+add_subdirectory(sky)
 
 if(WITH_AUDASPACE)
   add_subdirectory(audaspace)
diff --git a/intern/cycles/app/CMakeLists.txt b/intern/cycles/app/CMakeLists.txt
index ef374f91a65..a2b0ed03925 100644
--- a/intern/cycles/app/CMakeLists.txt
+++ b/intern/cycles/app/CMakeLists.txt
@@ -35,7 +35,7 @@ if(WITH_CYCLES_OSL)
 endif()
 
 if(NOT CYCLES_STANDALONE_REPOSITORY)
-  list(APPEND LIBRARIES bf_intern_glew_mx bf_intern_guardedalloc bf_intern_numaapi)
+  list(APPEND LIBRARIES bf_intern_glew_mx bf_intern_guardedalloc bf_intern_numaapi bf_intern_sky)
 endif()
 
 if(WITH_CYCLES_LOGGING)
diff --git a/intern/cycles/render/CMakeLists.txt b/intern/cycles/render/CMakeLists.txt
index e37a0407976..6a1335dc5dd 100644
--- a/intern/cycles/render/CMakeLists.txt
+++ b/intern/cycles/render/CMakeLists.txt
@@ -2,6 +2,7 @@
 set(INC
   ..
   ../../glew-mx
+  ../../sky/include
 )
 
 set(INC_SYS
@@ -92,6 +93,7 @@ set(LIB
   cycles_device
   cycles_subd
   cycles_util
+  bf_intern_sky
 )
 
 if(WITH_CYCLES_OSL)
diff --git a/intern/cycles/render/image_sky.cpp b/intern/cycles/render/image_sky.cpp
index 442e1d7941f..24d4834c2fa 100644
--- a/intern/cycles/render/image_sky.cpp
+++ b/intern/cycles/render/image_sky.cpp
@@ -16,10 +16,11 @@
 
 #include "render/image_sky.h"
 
+#include "sky_model.h"
+
 #include "util/util_image.h"
 #include "util/util_logging.h"
 #include "util/util_path.h"
-#include "util/util_sky_model.h"
 #include "util/util_task.h"
 
 CCL_NAMESPACE_BEGIN
@@ -62,17 +63,17 @@ bool SkyLoader::load_pixels(const ImageMetaData &metadata,
   const int rows_per_task = divide_up(1024, width);
   parallel_for(blocked_range<size_t>(0, height, rows_per_task),
                [&](const blocked_range<size_t> &r) {
-                 nishita_skymodel_precompute_texture(pixel_data,
-                                                     metadata.channels,
-                                                     r.begin(),
-                                                     r.end(),
-                                                     width,
-                                                     height,
-                                                     sun_elevation,
-                                                     altitude_f,
-                                                     air_density,
-                                                     dust_density,
-                                                     ozone_density);
+                 SKY_nishita_skymodel_precompute_texture(pixel_data,
+                                                         metadata.channels,
+                                                         r.begin(),
+                                                         r.end(),
+                                                         width,
+                                                         height,
+                                                         sun_elevation,
+                                                         altitude_f,
+                                                         air_density,
+                                                         dust_density,
+                                                         ozone_density);
                });
 
   return true;
diff --git a/intern/cycles/render/nodes.cpp b/intern/cycles/render/nodes.cpp
index ab392839e52..1e09b71b340 100644
--- a/intern/cycles/render/nodes.cpp
+++ b/intern/cycles/render/nodes.cpp
@@ -27,9 +27,10 @@
 #include "render/scene.h"
 #include "render/svm.h"
 
+#include "sky_model.h"
+
 #include "util/util_foreach.h"
 #include "util/util_logging.h"
-#include "util/util_sky_model.h"
 #include "util/util_transform.h"
 
 #include "kernel/svm/svm_color_util.h"
@@ -726,8 +727,8 @@ static void sky_texture_precompute_hosek(SunSky *sunsky,
   float solarElevation = M_PI_2_F - theta;
 
   /* Initialize Sky Model */
-  ArHosekSkyModelState *sky_state;
-  sky_state = arhosek_xyz_skymodelstate_alloc_init(
+  SKY_ArHosekSkyModelState *sky_state;
+  sky_state = SKY_arhosek_xyz_skymodelstate_alloc_init(
       (double)turbidity, (double)ground_albedo, (double)solarElevation);
 
   /* Copy values from sky_state to SunSky */
@@ -741,7 +742,7 @@ static void sky_texture_precompute_hosek(SunSky *sunsky,
   sunsky->radiance_z = (float)sky_state->radiances[2];
 
   /* Free sky_state */
-  arhosekskymodelstate_free(sky_state);
+  SKY_arhosekskymodelstate_free(sky_state);
 }
 
 /* Nishita improved */
@@ -758,7 +759,7 @@ static void sky_texture_precompute_nishita(SunSky *sunsky,
   float pixel_bottom[3];
   float pixel_top[3];
   float altitude_f = (float)altitude;
-  nishita_skymodel_precompute_sun(
+  SKY_nishita_skymodel_precompute_sun(
       sun_elevation, sun_size, altitude_f, air_density, dust_density, pixel_bottom, pixel_top);
   /* send data to svm_sky */
   sunsky->nishita_data[0] = pixel_bottom[0];
diff --git a/intern/cycles/util/CMakeLists.txt b/intern/cycles/util/CMakeLists.txt
index ad4ea9c86e0..23a47e064e2 100644
--- a/intern/cycles/util/CMakeLists.txt
+++ b/intern/cycles/util/CMakeLists.txt
@@ -98,10 +98,6 @@ set(SRC_HEADERS
   util_rect.h
   util_set.h
   util_simd.h
-  util_sky_model.cpp
-  util_sky_model.h
-  util_sky_model_data.h
-  util_sky_nishita.cpp
   util_avxf.h
   util_avxb.h
   util_semaphore.h
diff --git a/intern/sky/CMakeLists.txt b/intern/sky/CMakeLists.txt
new file mode 100644
index 00000000000..d46880367dc
--- /dev/null
+++ b/intern/sky/CMakeLists.txt
@@ -0,0 +1,35 @@
+# ***** BEGIN GPL LICENSE BLOCK *****
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+# ***** END GPL LICENSE BLOCK *****
+
+set(INC
+  include
+)
+
+set(INC_SYS
+
+)
+
+set(SRC
+  source/sky_model.cpp
+  source/sky_nishita.cpp
+)
+
+set(LIB
+)
+
+blender_add_lib(bf_intern_sky "${SRC}" "${INC}" "${INC_SYS}" "${LIB}")
diff --git a/intern/cycles/util/util_sky_model.h b/intern/sky/include/sky_model.h
index 36f1079a16d..9e7700f3042 100644
--- a/intern/cycles/util/util_sky_model.h
+++ b/intern/sky/include/sky_model.h
@@ -298,14 +298,14 @@ HINT #1:   if you want to model the sky of an earth-like planet that orbits
            previous paragraph.
 */
 
-#include "util/util_types.h"
+#ifndef __SKY_MODEL_H__
+#define __SKY_MODEL_H__
 
-CCL_NAMESPACE_BEGIN
+#ifdef __cplusplus
+extern "C" {
+#endif
 
-#ifndef _SKY_MODEL_H_
-#  define _SKY_MODEL_H_
-
-typedef double ArHosekSkyModelConfiguration[9];
+typedef double SKY_ArHosekSkyModelConfiguration[9];
 
 //   Spectral version of the model
 
@@ -335,8 +335,8 @@ typedef double ArHosekSkyModelConfiguration[9];
 
 ---------------------------------------------------------------------------- */
 
-typedef struct ArHosekSkyModelState {
-  ArHosekSkyModelConfiguration configs[11];
+typedef struct SKY_ArHosekSkyModelState {
+  SKY_ArHosekSkyModelConfiguration configs[11];
   double radiances[11];
   double turbidity;
   double solar_radius;
@@ -344,7 +344,7 @@ typedef struct ArHosekSkyModelState {
   double emission_correction_factor_sun[11];
   double albedo;
   double elevation;
-} ArHosekSkyModelState;
+} SKY_ArHosekSkyModelState;
 
 /* ----------------------------------------------------------------------------
 
@@ -355,7 +355,7 @@ typedef struct ArHosekSkyModelState {
 
 ---------------------------------------------------------------------------- */
 
-ArHosekSkyModelState *arhosekskymodelstate_alloc_init(const double solar_elevation,
+SKY_ArHosekSkyModelState *SKY_arhosekskymodelstate_alloc_init(const double solar_elevation,
                                                       const double atmospheric_turbidity,
                                                       const double ground_albedo);
 
@@ -388,31 +388,31 @@ ArHosekSkyModelState *arhosekskymodelstate_alloc_init(const double solar_elevati
 
 ---------------------------------------------------------------------------- */
 
-ArHosekSkyModelState *arhosekskymodelstate_alienworld_alloc_init(
+SKY_ArHosekSkyModelState *SKY_arhosekskymodelstate_alienworld_alloc_init(
     const double solar_elevation,
     const double solar_intensity,
     const double solar_surface_temperature_kelvin,
     const double atmospheric_turbidity,
     const double ground_albedo);
 
-void arhosekskymodelstate_free(ArHosekSkyModelState *state);
+void SKY_arhosekskymodelstate_free(SKY_ArHosekSkyModelState *state);
 
-double arhosekskymodel_radiance(ArHosekSkyModelState *state,
+double SKY_arhosekskymodel_radiance(SKY_ArHosekSkyModelState *state,
                                 double theta,
                                 double gamma,
                                 double wavelength);
 
 // CIE XYZ and RGB versions
 
-ArHosekSkyModelState *arhosek_xyz_skymodelstate_alloc_init(const double turbidity,
+SKY_ArHosekSkyModelState *SKY_arhosek_xyz_skymodelstate_alloc_init(const double turbidity,
                                                            const double albedo,
                                                            const double elevation);
 
-ArHosekSkyModelState *arhosek_rgb_skymodelstate_alloc_init(const double turbidity,
+SKY_ArHosekSkyModelState *SKY_arhosek_rgb_skymodelstate_alloc_init(const double turbidity,
                                                            const double albedo,
                                                            const double elevation);
 
-double arhosek_tristim_skymodel_radiance(ArHosekSkyModelState *state,
+double SKY_arhosek_tristim_skymodel_radiance(SKY_ArHosekSkyModelState *state,
                                          double theta,
                                          double gamma,
                                          int channel);
@@ -421,16 +421,15 @@ double arhosek_tristim_skymodel_radiance(ArHosekSkyModelState *state,
 //   Please read the above description before using this - there are several
 //   caveats!
 
-double arhosekskymodel_solar_radiance(ArHosekSkyModelState *state,
+double SKY_arhosekskymodel_solar_radiance(SKY_ArHosekSkyModelState *state,
                                       double theta,
                                       double gamma,
                                       double wavelength);
 
-#endif  // _SKY_MODEL_H_
 
 /* Nishita improved sky model */
 
-void nishita_skymodel_precompute_texture(float *pixels,
+void SKY_nishita_skymodel_precompute_texture(float *pixels,
                                          int stride,
                                          int start_y,
                                          int end_y,
@@ -442,7 +441,7 @@ void nishita_skymodel_precompute_texture(float *pixels,
                                          float dust_density,
                                          float ozone_density);
 
-void nishita_skymodel_precompute_sun(float sun_elevation,
+void SKY_nishita_skymodel_precompute_sun(float sun_elevation,
                                      float angular_diameter,
                                      float altitude,
                                      float air_density,
@@ -450,4 +449,8 @@ void nishita_skymodel_precompute_sun(float sun_elevation,
                                      float *pixel_bottom,
                                      float *pixel_top);
 
-CCL_NAMESPACE_END
+#ifdef __cplusplus
+}
+#endif
+
+#endif  // __SKY_MODEL_H__
diff --git a/intern/sky/source/sky_float3.h b/intern/sky/source/sky_float3.h
new file mode 100644
index 00000000000..2a9b9c89623
--- /dev/null
+++ b/intern/sky/source/sky_float3.h
@@ -0,0 +1,157 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#ifndef __SKY_FLOAT3_H__
+#define __SKY_FLOAT3_H__
+
+// minimal float3 + util_math.h implementation for nishita sky model
+
+#include <math.h>
+
+#ifndef M_PI_F
+#  define M_PI_F (3.1415926535897932f) /* pi */
+#endif
+#ifndef M_PI_2_F
+#  define M_PI_2_F (1.5707963267948966f) /* pi/2 */
+#endif
+#ifndef M_2PI_F
+#  define M_2PI_F (6.2831853071795864f) /* 2*pi */
+#endif
+
+struct float3 {
+  float x, y, z;
+
+  float3() = default;
+
+  float3(const float *ptr) : x{ptr[0]}, y{ptr[1]}, z{ptr[2]}
+  {
+  }
+
+  float3(const float (*ptr)[3]) : float3((const float *)ptr)
+  {
+  }
+
+  explicit float3(float value) : x(value), y(value), z(value)
+  {
+  }
+
+  explicit float3(int value) : x(value), y(value), z(value)
+  {
+  }
+
+  float3(float x, float y, float z) : x{x}, y{y}, z{z}
+  {
+  }
+
+  operator const float *() const
+  {
+    return &x;
+  }
+
+  operator float *()
+  {
+    return &x;
+  }
+
+  friend float3 operator*(const float3 &a, float b)
+  {
+    return {a.x * b, a.y * b, a.z * b};
+  }
+
+  friend float3 operator*(float b, const float3 &a)
+  {
+    return {a.x * b, a.y * b, a.z * b};
+  }
+
+  friend float3 operator-(const float3 &a, const float3 &b)
+  {
+    return {a.x - b.x, a.y - b.y, a.z - b.z};
+  }
+
+  friend float3 operator-(const float3 &a)
+  {
+    return {-a.x, -a.y, -a.z};
+  }
+
+  float length_squared() const
+  {
+    return x * x + y * y + z * z;
+  }
+
+  float length() const
+  {
+    return sqrt(length_squared());
+  }
+
+  static float distance(const float3 &a, const float3 &b)
+  {
+    return (a - b).length();
+  }
+
+  friend float3 operator+(const float3 &a, const float3 &b)
+  {
+    return {a.x + b.x, a.y + b.y, a.z + b.z};
+  }
+
+  void operator+=(const float3 &b)
+  {
+    this->x += b.x;
+    this->y += b.y;
+    this->z += b.z;
+  }
+
+  friend float3 operator*(const float3 &a, const float3 &b)
+  {
+    return {a.x * b.x, a.y * b.y, a.z * b.z};
+  }
+};
+
+inline float sqr(float a)
+{
+  return a * a;
+}
+
+inline float3 make_float3(float x, float y, float z)
+{
+  return float3(x, y, z);
+}
+
+inline float dot(const float3 &a, const float3 &b)
+{
+  return a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+inline float distance(const float3 &a, const float3 &b)
+{
+  return float3::distance(a, b);
+}
+
+inline float len_squared(float3 f)
+{
+  return f.length_squared();
+}
+
+inline float len(float3 f)
+{
+  return f.length();
+}
+
+inline float reduce_add(float3 f)
+{
+  return f.x + f.y + f.z;
+}
+
+#endif /* __SKY_FLOAT3_H__ */
diff --git a/intern/cycles/util/util_sky_model.cpp b/intern/sky/source/sky_model.cpp
index 8cdad8a90a4..3eca68e076c 100644
--- a/intern/cycles/util/util_sky_model.cpp
+++ b/intern/sky/source/sky_model.cpp
@@ -97,16 +97,14 @@ All instructions on how to use this code are in the accompanying header file.
 
 */
 
-#include "util/util_sky_model.h"
-#include "util/util_sky_model_data.h"
+#include "sky_model.h"
+#include "sky_model_data.h"
 
 #include <assert.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 
-CCL_NAMESPACE_BEGIN
-
 //   Some macro definitions that occur elsewhere in ART, and that have to be
 //   replicated to make this a stand-alone module.
 
@@ -138,7 +136,7 @@ typedef const double *ArHosekSkyModel_Radiance_Dataset;
 // internal functions
 
 static void ArHosekSkyModel_CookConfiguration(ArHosekSkyModel_Dataset dataset,
-                                              ArHosekSkyModelConfiguration config,
+                                              SKY_ArHosekSkyModelConfiguration config,
                                               double turbidity,
                                               double albedo,
                                               double solar_elevation)
@@ -272,7 +270,7 @@ static double ArHosekSkyModel_CookRadianceConfiguration(ArHosekSkyModel_Radiance
   return res;
 }
 
-static double ArHosekSkyModel_GetRadianceInternal(ArHosekSkyModelConfiguration configuration,
+static double ArHosekSkyModel_GetRadianceInternal(SKY_ArHosekSkyModelConfiguration configuration,
                                                   double theta,
                                                   double gamma)
 {
@@ -288,12 +286,12 @@ static double ArHosekSkyModel_GetRadianceInternal(ArHosekSkyModelConfiguration c
           configuration[6] * mieM + configuration[7] * zenith);
 }
 
-void arhosekskymodelstate_free(ArHosekSkyModelState *state)
+void SKY_arhosekskymodelstate_free(SKY_ArHosekSkyModelState *state)
 {
   free(state);
 }
 
-double arhosekskymodel_radiance(ArHosekSkyModelState *state,
+double SKY_arhosekskymodel_radiance(SKY_ArHosekSkyModelState *state,
                                 double theta,
                                 double gamma,
                                 double wavelength)
@@ -324,11 +322,11 @@ double arhosekskymodel_radiance(ArHosekSkyModelState *state,
 
 // xyz and rgb versions
 
-ArHosekSkyModelState *arhosek_xyz_skymodelstate_alloc_init(const double turbidity,
+SKY_ArHosekSkyModelState *SKY_arhosek_xyz_skymodelstate_alloc_init(const double turbidity,
                                                            const double albedo,
                                                            const double elevation)
 {
-  ArHosekSkyModelState *state = ALLOC(ArHosekSkyModelState);
+  SKY_ArHosekSkyModelState *state = ALLOC(SKY_ArHosekSkyModelState);
 
   state->solar_radius = TERRESTRIAL_SOLAR_RADIUS;
   state->turbidity = turbidity;
@@ -345,5 +343,3 @@ ArHosekSkyModelState *arhosek_xyz_skymodelstate_alloc_init(const double turbidit
 
   return state;
 }
-
-CCL_NAMESPACE_END
diff --git a/intern/cycles/util/util_sky_model_data.h b/intern/sky/source/sky_model_data.h
index a2a3935eb84..8d98f84cdae 100644
--- a/intern/cycles/util/util_sky_model_data.h
+++ b/intern/sky/source/sky_model_data.h
@@ -91,8 +91,6 @@ an updated version of this code has been published!
 
 ============================================================================ */
 
-CCL_NAMESPACE_BEGIN
-
 /*
 
 This file contains the coefficient data for the XYZ colour space version of
@@ -3843,5 +3841,3 @@ static const double datasetXYZRad3[] = {
 static const double *datasetsXYZ[] = {datasetXYZ1, datasetXYZ2, datasetXYZ3};
 
 static const double *datasetsXYZRad[] = {datasetXYZRad1, datasetXYZRad2, datasetXYZRad3};
-
-CCL_NAMESPACE_END
diff --git a/intern/cycles/util/util_sky_nishita.cpp b/intern/sky/source/sky_nishita.cpp
index 92397804d43..27286ddecac 100644
--- a/intern/cycles/util/util_sky_nishita.cpp
+++ b/intern/sky/source/sky_nishita.cpp
@@ -14,10 +14,8 @@
  * limitations under the License.
  */
 
-#include "util/util_math.h"
-#include "util/util_sky_model.h"
-
-CCL_NAMESPACE_BEGIN
+#include "sky_model.h"
+#include "sky_float3.h"
 
 /* Constants */
 static const float rayleigh_scale = 8000.0f;        // Rayleigh scale height (m)
@@ -269,7 +267,7 @@ static void single_scattering(float3 ray_dir,
 }
 
 /* calculate texture array */
-void nishita_skymodel_precompute_texture(float *pixels,
+void SKY_nishita_skymodel_precompute_texture(float *pixels,
                                          int stride,
                                          int start_y,
                                          int end_y,
@@ -332,7 +330,7 @@ static void sun_radiation(float3 cam_dir,
   }
 }
 
-void nishita_skymodel_precompute_sun(float sun_elevation,
+void SKY_nishita_skymodel_precompute_sun(float sun_elevation,
                                      float angular_diameter,
                                      float altitude,
                                      float air_density,
@@ -367,5 +365,3 @@ void nishita_skymodel_precompute_sun(float sun_elevation,
   pixel_top[1] = pix_top.y;
   pixel_top[2] = pix_top.z;
 }
-
-CCL_NAMESPACE_END
diff --git a/source/blender/gpu/shaders/material/gpu_shader_material_tex_sky.glsl b/source/blender/gpu/shaders/material/gpu_shader_material_tex_sky.glsl
index 981d17b4283..b2f6147b629 100644
--- a/source/blender/gpu/shaders/material/gpu_shader_material_tex_sky.glsl
+++ b/source/blender/gpu/shaders/material/gpu_shader_material_tex_sky.glsl
@@ -1,4 +1,150 @@
-void node_tex_sky(vec3 co, out vec4 color)
+float sky_angle_between(float thetav, float phiv, float theta, float phi)
+{
+  float cospsi = sin(thetav) * sin(theta) * cos(phi - phiv) + cos(thetav) * cos(theta);
+
+  if (cospsi > 1.0) {
+    return 0.0;
+  }
+  if (cospsi < -1.0) {
+    return M_PI;
+  }
+
+  return acos(cospsi);
+}
+
+vec3 sky_spherical_coordinates(vec3 dir)
+{
+  return vec3(acos(dir.z), atan(dir.x, dir.y), 0);
+}
+
+/* Preetham */
+/* lam03+lam4: 5 floats passed as vec4+float */
+float sky_perez_function(vec4 lam03, float lam4, float theta, float gamma)
+{
+  float ctheta = cos(theta);
+  float cgamma = cos(gamma);
+
+  return (1.0 + lam03[0] * exp(lam03[1] / ctheta)) *
+         (1.0 + lam03[2] * exp(lam03[3] * gamma) + lam4 * cgamma * cgamma);
+}
+
+vec3 xyY_to_xyz(float x, float y, float Y)
+{
+  float X, Z;
+
+  if (y != 0.0) {
+    X = (x / y) * Y;
+  }
+  else {
+    X = 0.0;
+  }
+
+  if (y != 0.0 && Y != 0.0) {
+    Z = ((1.0 - x - y) / y) * Y;
+  }
+  else {
+    Z = 0.0;
+  }
+
+  return vec3(X, Y, Z);
+}
+
+void node_tex_sky_preetham(vec3 co,
+                           vec4 config_Y03,
+                           float config_Y4,
+                           vec4 config_x03,
+                           float config_x4,
+                           vec4 config_y03,
+                           float config_y4,
+                           vec2 sun_angles,
+                           vec3 radiance,
+                           vec3 xyz_to_r,
+                           vec3 xyz_to_g,
+                           vec3 xyz_to_b,
+                           out vec4 color)
+{
+  /* convert vector to spherical coordinates */
+  vec3 spherical = sky_spherical_coordinates(co);
+  float theta = spherical[0];
+  float phi = spherical[1];
+
+  float suntheta = sun_angles[0];
+  float sunphi = sun_angles[1];
+
+  /* angle between sun direction and dir */
+  float gamma = sky_angle_between(theta, phi, suntheta, sunphi);
+
+  /* clamp theta to horizon */
+  theta = min(theta, M_PI_2 - 0.001);
+
+  /* compute xyY color space values */
+  float Y = radiance[0] * sky_perez_function(config_Y03, config_Y4, theta, gamma);
+  float x = radiance[1] * sky_perez_function(config_x03, config_x4, theta, gamma);
+  float y = radiance[2] * sky_perez_function(config_y03, config_y4, theta, gamma);
+
+  /* convert to RGB */
+  vec3 xyz = xyY_to_xyz(x, y, Y);
+  color = vec4(dot(xyz_to_r, xyz), dot(xyz_to_g, xyz), dot(xyz_to_b, xyz), 1);
+}
+
+/* Hosek / Wilkie */
+float sky_radiance_hosekwilkie(
+    vec4 config03, vec4 config47, float config8, float theta, float gamma)
+{
+  float ctheta = cos(theta);
+  float cgamma = cos(gamma);
+
+  float expM = exp(config47[0] * gamma);
+  float rayM = cgamma * cgamma;
+  float mieM = (1.0 + rayM) / pow((1.0 + config8 * config8 - 2.0 * config8 * cgamma), 1.5);
+  float zenith = sqrt(ctheta);
+
+  return (1.0 + config03[0] * exp(config03[1] / (ctheta + 0.01))) *
+         (config03[2] + config03[3] * expM + config47[1] * rayM + config47[2] * mieM +
+          config47[3] * zenith);
+}
+
+void node_tex_sky_hosekwilkie(vec3 co,
+                              vec4 config_x03,
+                              vec4 config_x47,
+                              vec4 config_y03,
+                              vec4 config_y47,
+                              vec4 config_z03,
+                              vec4 config_z47,
+                              vec3 config_xyz8,
+                              vec2 sun_angles,
+                              vec3 radiance,
+                              vec3 xyz_to_r,
+                              vec3 xyz_to_g,
+                              vec3 xyz_to_b,
+                              out vec4 color)
+{
+  /* convert vector to spherical coordinates */
+  vec3 spherical = sky_spherical_coordinates(co);
+  float theta = spherical[0];
+  float phi = spherical[1];
+
+  float suntheta = sun_angles[0];
+  float sunphi = sun_angles[1];
+
+  /* angle between sun direction and dir */
+  float gamma = sky_angle_between(theta, phi, suntheta, sunphi);
+
+  /* clamp theta to horizon */
+  theta = min(theta, M_PI_2 - 0.001);
+
+  vec3 xyz;
+  xyz.x = sky_radiance_hosekwilkie(config_x03, config_x47, config_xyz8[0], theta, gamma) *
+          radiance.x;
+  xyz.y = sky_radiance_hosekwilkie(config_y03, config_y47, config_xyz8[1], theta, gamma) *
+          radiance.y;
+  xyz.z = sky_radiance_hosekwilkie(config_z03, config_z47, config_xyz8[2], theta, gamma) *
+          radiance.z;
+
+  color = vec4(dot(xyz_to_r, xyz), dot(xyz_to_g, xyz), dot(xyz_to_b, xyz), 1);
+}
+
+void node_tex_sky_nishita(vec3 co, out vec4 color)
 {
   color = vec4(1.0);
 }
diff --git a/source/blender/imbuf/IMB_colormanagement.h b/source/blender/imbuf/IMB_colormanagement.h
index 3158e3419b0..4530f6c9fc0 100644
--- a/source/blender/imbuf/IMB_colormanagement.h
+++ b/source/blender/imbuf/IMB_colormanagement.h
@@ -72,6 +72,7 @@ BLI_INLINE float IMB_colormanagement_get_luminance(const float rgb[3]);
 BLI_INLINE unsigned char IMB_colormanagement_get_luminance_byte(const unsigned char[3]);
 BLI_INLINE void IMB_colormangement_xyz_to_rgb(float rgb[3], const float xyz[3]);
 BLI_INLINE void IMB_colormangement_rgb_to_xyz(float xyz[3], const float rgb[3]);
+const float *IMB_colormangement_get_xyz_to_rgb(void);
 
 /* ** Color space transformation functions ** */
 void IMB_colormanagement_transform(float *buffer,
diff --git a/source/blender/imbuf/intern/colormanagement.c b/source/blender/imbuf/intern/colormanagement.c
index 3f5a0f25cc5..6341706bede 100644
--- a/source/blender/imbuf/intern/colormanagement.c
+++ b/source/blender/imbuf/intern/colormanagement.c
@@ -1456,6 +1456,11 @@ bool IMB_colormanagement_space_name_is_data(const char *name)
   return (colorspace && colorspace->is_data);
 }
 
+const float *IMB_colormangement_get_xyz_to_rgb()
+{
+  return &imbuf_xyz_to_rgb[0][0];
+}
+
 /*********************** Threaded display buffer transform routines *************************/
 
 typedef struct DisplayBufferThread {
diff --git a/source/blender/nodes/CMakeLists.txt b/source/blender/nodes/CMakeLists.txt
index 55c5009ea56..2c4ca087efa 100644
--- a/source/blender/nodes/CMakeLists.txt
+++ b/source/blender/nodes/CMakeLists.txt
@@ -38,6 +38,7 @@ set(INC
   ../render/extern/include
   ../../../intern/glew-mx
   ../../../intern/guardedalloc
+  ../../../intern/sky/include
 )
 
 set(INC_SYS
@@ -298,6 +299,7 @@ set(SRC
 
 set(LIB
   bf_functions
+  bf_intern_sky
 )
 
 if(WITH_PYTHON)
diff --git a/source/blender/nodes/shader/nodes/node_shader_tex_sky.c b/source/blender/nodes/shader/nodes/node_shader_tex_sky.c
index 0daa948c139..70f3f80eb4b 100644
--- a/source/blender/nodes/shader/nodes/node_shader_tex_sky.c
+++ b/source/blender/nodes/shader/nodes/node_shader_tex_sky.c
@@ -18,6 +18,7 @@
  */
 
 #include "../node_shader_util.h"
+#include "sky_model.h"
 
 /* **************** OUTPUT ******************** */
 
@@ -53,19 +54,164 @@ static void node_shader_init_tex_sky(bNodeTree *UNUSED(ntree), bNode *node)
   node->storage = tex;
 }
 
+typedef struct SkyModelPreetham {
+  float config_Y[5], config_x[5], config_y[5]; /* named after xyY color space */
+  float radiance[3];
+} SkyModelPreetham;
+
+typedef struct XYZ_to_RGB /* transposed imbuf_xyz_to_rgb, passed as 3x vec3 */
+{
+  float r[3], g[3], b[3];
+} XYZ_to_RGB;
+
+static float sky_perez_function(float *lam, float theta, float gamma)
+{
+  float ctheta = cosf(theta);
+  float cgamma = cosf(gamma);
+
+  return (1.0 + lam[0] * expf(lam[1] / ctheta)) *
+         (1.0 + lam[2] * expf(lam[3] * gamma) + lam[4] * cgamma * cgamma);
+}
+
+static void sky_precompute_old(SkyModelPreetham *sunsky, float sun_angles[], float turbidity)
+{
+  float theta = sun_angles[0];
+  float theta2 = theta * theta;
+  float theta3 = theta2 * theta;
+  float T = turbidity;
+  float T2 = T * T;
+  float chi = (4.0f / 9.0f - T / 120.0f) * (M_PI - 2.0f * theta);
+
+  sunsky->radiance[0] = (4.0453f * T - 4.9710f) * tanf(chi) - 0.2155f * T + 2.4192f;
+  sunsky->radiance[0] *= 0.06f;
+
+  sunsky->radiance[1] = (0.00166f * theta3 - 0.00375f * theta2 + 0.00209f * theta) * T2 +
+                        (-0.02903f * theta3 + 0.06377f * theta2 - 0.03202f * theta + 0.00394f) *
+                            T +
+                        (0.11693f * theta3 - 0.21196f * theta2 + 0.06052f * theta + 0.25886f);
+
+  sunsky->radiance[2] = (0.00275f * theta3 - 0.00610f * theta2 + 0.00317f * theta) * T2 +
+                        (-0.04214f * theta3 + 0.08970f * theta2 - 0.04153f * theta + 0.00516f) *
+                            T +
+                        (0.15346f * theta3 - 0.26756f * theta2 + 0.06670f * theta + 0.26688f);
+
+  sunsky->config_Y[0] = (0.1787f * T - 1.4630f);
+  sunsky->config_Y[1] = (-0.3554f * T + 0.4275f);
+  sunsky->config_Y[2] = (-0.0227f * T + 5.3251f);
+  sunsky->config_Y[3] = (0.1206f * T - 2.5771f);
+  sunsky->config_Y[4] = (-0.0670f * T + 0.3703f);
+
+  sunsky->config_x[0] = (-0.0193f * T - 0.2592f);
+  sunsky->config_x[1] = (-0.0665f * T + 0.0008f);
+  sunsky->config_x[2] = (-0.0004f * T + 0.2125f);
+  sunsky->config_x[3] = (-0.0641f * T - 0.8989f);
+  sunsky->config_x[4] = (-0.0033f * T + 0.0452f);
+
+  sunsky->config_y[0] = (-0.0167f * T - 0.2608f);
+  sunsky->config_y[1] = (-0.0950f * T + 0.0092f);
+  sunsky->config_y[2] = (-0.0079f * T + 0.2102f);
+  sunsky->config_y[3] = (-0.0441f * T - 1.6537f);
+  sunsky->config_y[4] = (-0.0109f * T + 0.0529f);
+
+  sunsky->radiance[0] /= sky_perez_function(sunsky->config_Y, 0, theta);
+  sunsky->radiance[1] /= sky_perez_function(sunsky->config_x, 0, theta);
+  sunsky->radiance[2] /= sky_perez_function(sunsky->config_y, 0, theta);
+}
+
+static void get_XYZ_to_RGB_for_gpu(XYZ_to_RGB *data)
+{
+  const float *xyz_to_rgb = IMB_colormangement_get_xyz_to_rgb();
+  data->r[0] = xyz_to_rgb[0];
+  data->r[1] = xyz_to_rgb[3];
+  data->r[2] = xyz_to_rgb[6];
+  data->g[0] = xyz_to_rgb[1];
+  data->g[1] = xyz_to_rgb[4];
+  data->g[2] = xyz_to_rgb[7];
+  data->b[0] = xyz_to_rgb[2];
+  data->b[1] = xyz_to_rgb[5];
+  data->b[2] = xyz_to_rgb[8];
+}
+
 static int node_shader_gpu_tex_sky(GPUMaterial *mat,
                                    bNode *node,
                                    bNodeExecData *UNUSED(execdata),
                                    GPUNodeStack *in,
                                    GPUNodeStack *out)
 {
-  if (!in[0].link) {
-    in[0].link = GPU_attribute(mat, CD_ORCO, "");
-  }
-
+  node_shader_gpu_default_tex_coord(mat, node, &in[0].link);
   node_shader_gpu_tex_mapping(mat, node, in, out);
+  NodeTexSky *tex = (NodeTexSky *)node->storage;
+  float sun_angles[2]; /* [0]=theta=zenith angle  [1]=phi=azimuth */
+  sun_angles[0] = acosf(tex->sun_direction[2]);
+  sun_angles[1] = atan2f(tex->sun_direction[0], tex->sun_direction[1]);
 
-  return GPU_stack_link(mat, node, "node_tex_sky", in, out);
+  if (tex->sky_model == 0) {
+    /* Preetham */
+    SkyModelPreetham sunsky;
+    sky_precompute_old(&sunsky, sun_angles, tex->turbidity);
+    XYZ_to_RGB xyz_to_rgb;
+    get_XYZ_to_RGB_for_gpu(&xyz_to_rgb);
+    return GPU_stack_link(mat,
+                          node,
+                          "node_tex_sky_preetham",
+                          in,
+                          out,
+                          /* Pass config_Y/x/y as 3x(vec4+float) */
+                          GPU_uniform(&sunsky.config_Y[0]),
+                          GPU_uniform(&sunsky.config_Y[4]),
+                          GPU_uniform(&sunsky.config_x[0]),
+                          GPU_uniform(&sunsky.config_x[4]),
+                          GPU_uniform(&sunsky.config_y[0]),
+                          GPU_uniform(&sunsky.config_y[4]),
+                          GPU_uniform(sun_angles),
+                          GPU_uniform(sunsky.radiance),
+                          GPU_uniform(xyz_to_rgb.r),
+                          GPU_uniform(xyz_to_rgb.g),
+                          GPU_uniform(xyz_to_rgb.b));
+  }
+  else if (tex->sky_model == 1) {
+    /* Hosek / Wilkie */
+    sun_angles[0] = fmin(M_PI_2, sun_angles[0]); /* clamp to horizon */
+    SKY_ArHosekSkyModelState *sky_state = SKY_arhosek_xyz_skymodelstate_alloc_init(
+        tex->turbidity, tex->ground_albedo, fmax(0.0, M_PI_2 - sun_angles[0]));
+    /* Pass sky_state->configs[3][9] as 3*(vec4+vec4)+vec3 */
+    float config_x07[8], config_y07[8], config_z07[8], config_xyz8[3];
+    for (int i = 0; i < 8; ++i) {
+      config_x07[i] = (float)sky_state->configs[0][i];
+      config_y07[i] = (float)sky_state->configs[1][i];
+      config_z07[i] = (float)sky_state->configs[2][i];
+    }
+    for (int i = 0; i < 3; ++i) {
+      config_xyz8[i] = (float)sky_state->configs[i][8];
+    }
+    float radiance[3];
+    for (int i = 0; i < 3; i++) {
+      radiance[i] = sky_state->radiances[i] * (2 * M_PI / 683);
+    }
+    SKY_arhosekskymodelstate_free(sky_state);
+    XYZ_to_RGB xyz_to_rgb;
+    get_XYZ_to_RGB_for_gpu(&xyz_to_rgb);
+    return GPU_stack_link(mat,
+                          node,
+                          "node_tex_sky_hosekwilkie",
+                          in,
+                          out,
+                          GPU_uniform(&config_x07[0]),
+                          GPU_uniform(&config_x07[4]),
+                          GPU_uniform(&config_y07[0]),
+                          GPU_uniform(&config_y07[4]),
+                          GPU_uniform(&config_z07[0]),
+                          GPU_uniform(&config_z07[4]),
+                          GPU_uniform(config_xyz8),
+                          GPU_uniform(sun_angles),
+                          GPU_uniform(radiance),
+                          GPU_uniform(xyz_to_rgb.r),
+                          GPU_uniform(xyz_to_rgb.g),
+                          GPU_uniform(xyz_to_rgb.b));
+  }
+  else {
+    return GPU_stack_link(mat, node, "node_tex_sky_nishita", in, out);
+  }
 }
 
 /* node type definition */