diff --git a/Shaders/std/sky.glsl b/Shaders/std/sky.glsl index 61f93c84..5a58734f 100644 --- a/Shaders/std/sky.glsl +++ b/Shaders/std/sky.glsl @@ -1,12 +1,13 @@ /* Various sky functions * ===================== * - * Nishita model is based on https://github.com/wwwtyro/glsl-atmosphere(Unlicense License) + * 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 * - Implemented air, dust and ozone density node parameters (see Blender source) + * - Replaced the inner integral calculation with a LUT lookup * * Reference for the sun's limb darkening and ozone calculations: * [Hill] Sebastien Hillaire. Physically Based Sky, Atmosphere and Cloud Rendering in Frostbite @@ -19,15 +20,18 @@ #ifndef _SKY_GLSL_ #define _SKY_GLSL_ -// OpenGl ES doesn't support 1D textures so we use a 1 px height sampler2D here... uniform sampler2D nishitaLUT; -#define PI 3.141592 +#ifndef PI + #define PI 3.141592 +#endif +#ifndef HALF_PI + #define HALF_PI 1.570796 +#endif #define nishita_iSteps 16 -#define nishita_jSteps 8 -// The values here are taken from Cycles code if they +// These values 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 @@ -49,7 +53,17 @@ uniform sampler2D nishitaLUT; // Values from [Hill: 60] #define sun_limb_darkening_col vec3(0.397, 0.503, 0.652) -#define heightToLUT(h) (textureLod(nishitaLUT, vec2(clamp(h * (1 / 60000.0), 0.0, 1.0), 0.0), 0.0).xyz * 10.0) +float random(vec2 coords) { + return fract(sin(dot(coords.xy, vec2(12.9898,78.233))) * 43758.5453); +} + +vec3 nishita_lookupLUT(const float height, const float sunTheta) { + vec2 coords = vec2( + sqrt(height * (1 / nishita_atmo_radius)), + 0.5 + 0.5 * sign(sunTheta - HALF_PI) * sqrt(abs(sunTheta * (1 / HALF_PI) - 1)) + ); + return textureLod(nishitaLUT, coords, 0.0).rgb; +} /* Approximates the density of ozone for a given sample height. Values taken from Cycles code. */ float nishita_density_ozone(const float height) { @@ -112,40 +126,21 @@ vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const floa // Calculate the height of the sample. float iHeight = length(iPos) - rPlanet; - // Calculate the optical depth of the Rayleigh and Mie scattering for this step. - vec3 iLookup = heightToLUT(iHeight); - float odStepRlh = iLookup.x * iStepSize; - float odStepMie = iLookup.y * iStepSize; + // Calculate the optical depth of the Rayleigh and Mie scattering for this step + float odStepRlh = exp(-iHeight / nishita_rayleigh_scale) * density.x * iStepSize; + float odStepMie = exp(-iHeight / nishita_mie_scale) * density.y * iStepSize; // Accumulate optical depth. iOdRlh += odStepRlh; iOdMie += odStepMie; - // Calculate the step size of the secondary ray. - float jStepSize = nishita_rsi(iPos, pSun, nishita_atmo_radius).y / float(nishita_jSteps); + // Idea behind this: "Rotate" everything by iPos (-> iPos is the new zenith) and then all calculations for the + // inner integral only depend on the sample height (iHeight) and sunTheta (angle between sun and new zenith). + float sunTheta = acos(dot(normalize(iPos), normalize(pSun))); + vec3 jODepth = nishita_lookupLUT(iHeight, sunTheta);// * vec3(14000000 / 255, 14000000 / 255, 2000000 / 255); - // Initialize the secondary ray time. - float jTime = 0.0; - - // Initialize optical depth accumulators for the secondary ray. - vec3 jODepth = vec3(0.0); // (Rayleigh, Mie, ozone) - - // Sample the secondary ray. - for (int j = 0; j < nishita_jSteps; j++) { - - // Calculate the secondary ray sample position. - vec3 jPos = iPos + pSun * (jTime + jStepSize * 0.5); - - // Calculate the height of the sample. - float jHeight = length(jPos) - rPlanet; - - // Accumulate the optical depth. - vec3 jLookup = heightToLUT(jHeight); - jODepth += jLookup * jStepSize; - - // Increment the secondary ray time. - jTime += jStepSize; - } + // Apply dithering to reduce visible banding + jODepth += mix(-1000, 1000, random(r.xy)); // Calculate attenuation. vec3 attn = exp(-( diff --git a/Sources/armory/object/Uniforms.hx b/Sources/armory/object/Uniforms.hx index 64d20af2..e492f2fb 100644 --- a/Sources/armory/object/Uniforms.hx +++ b/Sources/armory/object/Uniforms.hx @@ -21,7 +21,7 @@ class Uniforms { public static function textureLink(object: Object, mat: MaterialData, link: String): kha.Image { if (link == "_nishitaLUT") { if (armory.renderpath.Nishita.data == null) armory.renderpath.Nishita.recompute(Scene.active.world); - return armory.renderpath.Nishita.data.optDepthLUT; + return armory.renderpath.Nishita.data.lut; } #if arm_ltc else if (link == "_ltcMat") { diff --git a/Sources/armory/renderpath/Nishita.hx b/Sources/armory/renderpath/Nishita.hx index c968798d..1a4cac92 100644 --- a/Sources/armory/renderpath/Nishita.hx +++ b/Sources/armory/renderpath/Nishita.hx @@ -5,36 +5,77 @@ import kha.graphics4.TextureFormat; import kha.graphics4.Usage; import iron.data.WorldData; +import iron.math.Vec2; +import iron.math.Vec3; import armory.math.Helper; +/** + Utility class to control the Nishita sky model. +**/ class Nishita { public static var data: NishitaData = null; + /** + Recompute the nishita lookup table. Call this function after updating + the sky density settings. + **/ public static function recompute(world: WorldData) { if (world == null || world.raw.sun_direction == null) return; if (data == null) data = new NishitaData(); // TODO - data.recompute(1.0, 1.0, 1.0); + data.computeLUT(new Vec3(1.0, 1.0, 1.0)); } } +/** + This class holds the precalculated result of the inner scattering integral + of the Nishita sky model. The outer integral is calculated in + [`armory/Shaders/std/sky.glsl`](https://github.com/armory3d/armory/blob/master/Shaders/std/sky.glsl). + + @see `armory.renderpath.Nishita` +**/ class NishitaData { - static inline var LUT_WIDTH = 16; - /** Maximum ray height as defined by Cycles **/ - static inline var MAX_HEIGHT = 60000; + public var lut: kha.Image; - static inline var RAYLEIGH_SCALE = 8e3; - static inline var MIE_SCALE = 1.2e3; + /** + The amount of individual sample heights stored in the LUT (and the width + of the LUT image). + **/ + public static var lutHeightSteps = 128; + /** + The amount of individual sun angle steps stored in the LUT (and the + height of the LUT image). + **/ + public static var lutAngleSteps = 128; - public var optDepthLUT: kha.Image; + /** + Amount of steps for calculating the inner scattering integral. Heigher + values are more precise but take longer to compute. + **/ + public static var jSteps = 8; + + /** Radius of the atmosphere in meters. **/ + public static var radiusAtmo = 6420000; + /** + Radius of the planet in meters. The default value is the earth radius as + defined in Cycles. + **/ + public static var radiusPlanet = 6360000; + + /** Rayleigh scattering scale parameter. **/ + public static var rayleighScale = 8e3; + /** Mie scattering scale parameter. **/ + public static var mieScale = 1.2e3; public function new() {} + /** Approximates the density of ozone for a given sample height. **/ function getOzoneDensity(height: FastFloat): FastFloat { + // Values are taken from Cycles code if (height < 10000.0 || height >= 40000.0) { return 0.0; } @@ -45,36 +86,91 @@ class NishitaData { } /** - The RGBA texture layout looks as follows: - R = Rayleigh optical depth at height \in [0, 60000] - G = Mie optical depth at height \in [0, 60000] - B = Ozone optical depth at height \in [0, 60000] - A = Unused + Ray-sphere intersection test that assumes the sphere is centered at the + origin. There is no intersection when result.x > result.y. Otherwise + this function returns the distances to the two intersection points, + which might be equal. **/ - public function recompute(densityFacAir: FastFloat, densityFacDust: FastFloat, densityFacOzone: FastFloat) { - optDepthLUT = kha.Image.create(LUT_WIDTH, 1, TextureFormat.RGBA32, Usage.StaticUsage); + function raySphereIntersection(rayOrigin: Vec3, rayDirection: Vec3, sphereRadius: Int): Vec2 { + // Algorithm is described here: https://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection + var a = rayDirection.dot(rayDirection); + var b = 2.0 * rayDirection.dot(rayOrigin); + var c = rayOrigin.dot(rayOrigin) - (sphereRadius * sphereRadius); + var d = (b * b) - 4.0 * a * c; - var textureData = optDepthLUT.lock(); - for (i in 0...LUT_WIDTH) { - // Get the height for each LUT pixel i (in [-1, 1] range) - var height = (i / LUT_WIDTH) * 2 - 1; + // Ray does not intersect the sphere + if (d < 0.0) return new Vec2(1e5, -1e5); + + return new Vec2( + (-b - Math.sqrt(d)) / (2.0 * a), + (-b + Math.sqrt(d)) / (2.0 * a) + ); + } + + /** + Computes the LUT texture for the given density values. + @param density 3D vector of air density, dust density, ozone density + **/ + public function computeLUT(density: Vec3) { + var imageData = new haxe.io.Float32Array(lutHeightSteps * lutAngleSteps * 4); + + for (x in 0...lutHeightSteps) { + var height = (x / (lutHeightSteps - 1)); // Use quadratic height for better horizon precision - // See https://sebh.github.io/publications/egsr2020.pdf (5.3) - height = 0.5 + 0.5 * Helper.sign(height) * Math.sqrt(Math.abs(height)); - height *= MAX_HEIGHT; // Denormalize + height *= height; + height *= radiusAtmo; // Denormalize - // Make sure we use 32 bit floats - var optDepthRayleigh: FastFloat = Math.exp(-height / RAYLEIGH_SCALE) * densityFacAir; - var optDepthMie: FastFloat = Math.exp(-height / MIE_SCALE) * densityFacDust; - var optDepthOzone: FastFloat = getOzoneDensity(height) * densityFacOzone; + for (y in 0...lutAngleSteps) { + var sunTheta = y / (lutAngleSteps - 1) * 2 - 1; - // 10 is the maximum density, so we divide by it to be able to use normalized values - textureData.set(i * 4 + 0, Std.int(optDepthRayleigh * 255 / 10)); - textureData.set(i * 4 + 1, Std.int(optDepthMie * 255 / 10)); - textureData.set(i * 4 + 2, Std.int(optDepthOzone * 255 / 10)); - textureData.set(i * 4 + 3, 255); // Unused + // Improve horizon precision + // See https://sebh.github.io/publications/egsr2020.pdf (5.3) + sunTheta = Helper.sign(sunTheta) * sunTheta * sunTheta; + sunTheta = sunTheta * Math.PI / 2 + Math.PI / 2; // Denormalize + + var jODepth = sampleSecondaryRay(height, sunTheta, density); + + var pixelIndex = (x + y * lutHeightSteps) * 4; + imageData[pixelIndex + 0] = jODepth.x; + imageData[pixelIndex + 1] = jODepth.y; + imageData[pixelIndex + 2] = jODepth.z; + imageData[pixelIndex + 3] = 1.0; // Unused + } } - optDepthLUT.unlock(); + + lut = kha.Image.fromBytes(imageData.view.buffer, lutHeightSteps, lutAngleSteps, TextureFormat.RGBA128, Usage.StaticUsage); + } + + /** + Calculates the integral for the secondary ray. + **/ + public function sampleSecondaryRay(height: FastFloat, sunTheta: FastFloat, density: Vec3): Vec3 { + // Reconstruct values from the shader + var iPos = new Vec3(0, 0, height + radiusPlanet); + var pSun = new Vec3(0.0, Math.sin(sunTheta), Math.cos(sunTheta)).normalize(); + + var jTime: FastFloat = 0.0; + var jStepSize: FastFloat = raySphereIntersection(iPos, pSun, radiusAtmo).y / jSteps; + + // Optical depth accumulators for the secondary ray (Rayleigh, Mie, ozone) + var jODepth = new Vec3(); + + for (i in 0...jSteps) { + + // Calculate the secondary ray sample position and height + var jPos = iPos.clone().add(pSun.clone().mult(jTime + jStepSize * 0.5)); + var jHeight = jPos.length() - radiusPlanet; + + // Accumulate optical depth + var optDepthRayleigh = Math.exp(-jHeight / rayleighScale) * density.x; + var optDepthMie = Math.exp(-jHeight / mieScale) * density.y; + var optDepthOzone = getOzoneDensity(jHeight) * density.z; + jODepth.addf(optDepthRayleigh, optDepthMie, optDepthOzone); + + jTime += jStepSize; + } + + return jODepth.mult(jStepSize); } } diff --git a/blender/arm/material/cycles_nodes/nodes_texture.py b/blender/arm/material/cycles_nodes/nodes_texture.py index b3b81d1f..c88a8373 100644 --- a/blender/arm/material/cycles_nodes/nodes_texture.py +++ b/blender/arm/material/cycles_nodes/nodes_texture.py @@ -373,6 +373,8 @@ def parse_sky_nishita(node: bpy.types.ShaderNodeTexSky, state: ParserState) -> v curshader = state.curshader curshader.add_include('std/sky.glsl') curshader.add_uniform('vec3 sunDir', link='_sunDirection') + curshader.add_uniform('sampler2D nishitaLUT', link='_nishitaLUT', included=True, + tex_addr_u='clamp', tex_addr_v='clamp') planet_radius = 6360e3 # Earth radius used in Blender ray_origin_z = planet_radius + node.altitude