Optimize Nishita sky shader

Shaders/std/sky.glsl

@@ -3,6 +3,7 @@
  * Nishita model is based on https://github.com/wwwtyro/glsl-atmosphere (Unlicense License)
  * Changes to the original implementation:
  *  - r and pSun parameters of nishita_atmosphere() are already normalized
+ *  - Some original parameters of nishita_atmosphere() are replaced with pre-defined values
  */
 
 #ifndef _SKY_GLSL_
@@ -13,6 +14,18 @@
 #define nishita_iSteps 16
 #define nishita_jSteps 8
 
+// The values here are taken from Cycles code if they
+// exist there, otherwise they are taken from the example
+// in the glsl-atmosphere repo
+#define nishita_sun_intensity 22.0
+#define nishita_atmo_radius 6420e3
+#define nishita_rayleigh_scale 8e3
+#define nishita_rayleigh_coeff vec3(5.5e-6, 13.0e-6, 22.4e-6)
+#define nishita_mie_scale 1.2e3
+#define nishita_mie_coeff 2e-5
+#define nishita_mie_dir 0.76 // Aerosols anisotropy ("direction")
+#define nishita_mie_dir_sq 0.5776 // Squared aerosols anisotropy
+
 /* ray-sphere intersection that assumes
  * the sphere is centered at the origin.
  * No intersection when result.x > result.y */
@@ -29,11 +42,15 @@ vec2 nishita_rsi(const vec3 r0, const vec3 rd, const float sr) {
 	);
 }
 
-vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const float iSun, const float rPlanet, const float rAtmos, const vec3 kRlh, const float kMie, const float shRlh, const float shMie, const float g) {
-	// r and pSun must be already normalized!
-
+/*
+ * r: normalized ray direction
+ * r0: ray origin
+ * pSun: normalized sun direction
+ * rPlanet: planet radius
+ */
+vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const float rPlanet) {
 	// Calculate the step size of the primary ray.
-	vec2 p = nishita_rsi(r0, r, rAtmos);
+	vec2 p = nishita_rsi(r0, r, nishita_atmo_radius);
 	if (p.x > p.y) return vec3(0,0,0);
 	p.y = min(p.y, nishita_rsi(r0, r, rPlanet).x);
 	float iStepSize = (p.y - p.x) / float(nishita_iSteps);
@@ -52,9 +69,8 @@ vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const floa
 	// Calculate the Rayleigh and Mie phases.
 	float mu = dot(r, pSun);
 	float mumu = mu * mu;
-	float gg = g * g;
 	float pRlh = 3.0 / (16.0 * PI) * (1.0 + mumu);
-	float pMie = 3.0 / (8.0 * PI) * ((1.0 - gg) * (mumu + 1.0)) / (pow(1.0 + gg - 2.0 * mu * g, 1.5) * (2.0 + gg));
+	float pMie = 3.0 / (8.0 * PI) * ((1.0 - nishita_mie_dir_sq) * (mumu + 1.0)) / (pow(1.0 + nishita_mie_dir_sq - 2.0 * mu * nishita_mie_dir, 1.5) * (2.0 + nishita_mie_dir_sq));
 
 	// Sample the primary ray.
 	for (int i = 0; i < nishita_iSteps; i++) {
@@ -66,15 +82,15 @@ vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const floa
 		float iHeight = length(iPos) - rPlanet;
 
 		// Calculate the optical depth of the Rayleigh and Mie scattering for this step.
-		float odStepRlh = exp(-iHeight / shRlh) * iStepSize;
-		float odStepMie = exp(-iHeight / shMie) * iStepSize;
+		float odStepRlh = exp(-iHeight / nishita_rayleigh_scale) * iStepSize;
+		float odStepMie = exp(-iHeight / nishita_mie_scale) * iStepSize;
 
 		// Accumulate optical depth.
 		iOdRlh += odStepRlh;
 		iOdMie += odStepMie;
 
 		// Calculate the step size of the secondary ray.
-		float jStepSize = nishita_rsi(iPos, pSun, rAtmos).y / float(nishita_jSteps);
+		float jStepSize = nishita_rsi(iPos, pSun, nishita_atmo_radius).y / float(nishita_jSteps);
 
 		// Initialize the secondary ray time.
 		float jTime = 0.0;
@@ -93,15 +109,15 @@ vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const floa
 			float jHeight = length(jPos) - rPlanet;
 
 			// Accumulate the optical depth.
-			jOdRlh += exp(-jHeight / shRlh) * jStepSize;
-			jOdMie += exp(-jHeight / shMie) * jStepSize;
+			jOdRlh += exp(-jHeight / nishita_rayleigh_scale) * jStepSize;
+			jOdMie += exp(-jHeight / nishita_mie_scale) * jStepSize;
 
 			// Increment the secondary ray time.
 			jTime += jStepSize;
 		}
 
 		// Calculate attenuation.
-		vec3 attn = exp(-(kMie * (iOdMie + jOdMie) + kRlh * (iOdRlh + jOdRlh)));
+		vec3 attn = exp(-(nishita_mie_coeff * (iOdMie + jOdMie) + nishita_rayleigh_coeff * (iOdRlh + jOdRlh)));
 
 		// Accumulate scattering.
 		totalRlh += odStepRlh * attn;
@@ -112,7 +128,7 @@ vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const floa
 	}
 
 	// Calculate and return the final color.
-	return iSun * (pRlh * kRlh * totalRlh + pMie * kMie * totalMie);
+	return nishita_sun_intensity * (pRlh * nishita_rayleigh_coeff * totalRlh + pMie * nishita_mie_coeff * totalMie);
 }
 
 #endif

blender/arm/material/cycles_nodes/nodes_texture.py

@@ -373,7 +373,10 @@ def parse_sky_nishita(node: bpy.types.ShaderNodeTexSky, state: ParserState) -> v
     curshader.add_include('std/sky.glsl')
     curshader.add_uniform('vec3 sunDir', link='_sunDirection')
 
-    return 'nishita_atmosphere(n, vec3(0,0,6372e3), sunDir, 22.0, 6371e3, 6471e3, vec3(5.5e-6,13.0e-6,22.4e-6), 21e-6, 8e3, 1.2e3, 0.758)'
+    planet_radius = 6360e3  # Earth radius used in Blender
+    ray_origin_z = planet_radius + node.altitude * 1000
+
+    return f'nishita_atmosphere(n, vec3(0, 0, {ray_origin_z}), sunDir, {planet_radius})'
 
 
 def parse_tex_environment(node: bpy.types.ShaderNodeTexEnvironment, out_socket: bpy.types.NodeSocket, state: ParserState) -> vec3str: