diff --git a/leenkx/Shaders/std/sky.glsl b/leenkx/Shaders/std/sky.glsl index cb79bb27..5243b3a3 100644 --- a/leenkx/Shaders/std/sky.glsl +++ b/leenkx/Shaders/std/sky.glsl @@ -185,19 +185,26 @@ vec3 multi_scatter_sample_lut(vec3 dir, float sun_rotation) { } vec3 multi_scatter_sun_disc(vec3 dir, vec3 sun_dir, float angular_diameter, float intensity) { - float dist = distance(dir, sun_dir) / (angular_diameter * 0.5); - if (dist > 1.0) return vec3(0.0); + float half_diameter = max(angular_diameter * 0.5, 0.0005); + float dist = distance(dir, sun_dir) / half_diameter; + float edge = smoothstep(1.0, 0.95, dist); + if (edge <= 0.0) return vec3(0.0); + + float horizon = multiScatterSunBottom.w; + float dir_elev = asin(clamp(dir.z, -1.0, 1.0)); + float horizon_fade = smoothstep(horizon - half_diameter * 0.1, horizon + half_diameter * 0.1, dir_elev); + edge *= horizon_fade; + if (edge <= 0.0) return vec3(0.0); float invDist = 1.0 - dist; float mu = sqrt(invDist * invDist); vec3 limb_darkening = 1.0 - (1.0 - pow(vec3(mu), sun_limb_darkening_col)); float sun_elev = multiScatterParams.x; - float dir_elev = asin(clamp(dir.z, -1.0, 1.0)); - float t = clamp((dir_elev - (sun_elev - angular_diameter * 0.5)) / angular_diameter, 0.0, 1.0); + float t = clamp((dir_elev - (sun_elev - half_diameter)) / (2.0 * half_diameter), 0.0, 1.0); vec3 sun_color = mix(multiScatterSunBottom.rgb, multiScatterSunTop, t) * intensity * limb_darkening; - return xyz_to_rgb(sun_color); + return xyz_to_rgb(sun_color) * edge; } vec3 multi_scatter_atmosphere(vec3 dir) { diff --git a/leenkx/Sources/leenkx/renderpath/Sky.hx b/leenkx/Sources/leenkx/renderpath/Sky.hx index 8ab18659..15b37f4e 100644 --- a/leenkx/Sources/leenkx/renderpath/Sky.hx +++ b/leenkx/Sources/leenkx/renderpath/Sky.hx @@ -6,9 +6,7 @@ import kha.graphics4.TextureFormat; import kha.graphics4.Usage; import iron.data.WorldData; -import iron.math.Vec2; import iron.math.Vec3; -import iron.math.Vec4; import leenkx.math.Helper; @@ -19,20 +17,60 @@ class Sky { public static var singleScatterData: SingleScatteringData = null; public static var multiScatterData: MultipleScatteringData = null; + // Recomputes the single scattering lookup table after the density settings changed. - /** - Recomputes the single scattering lookup table after the density settings changed. - Do not call this method on every frame (it's slow)! - **/ + static var registered = false; + static inline function register() { + registered = true; + iron.App.notifyOnUpdate(update); + iron.App.notifyOnReset(function() { + iron.App.notifyOnUpdate(update); + }); + } + + static var multiPrevData: haxe.io.Float32Array = null; + static var multiNewData: haxe.io.Float32Array = null; + static var multiBlendData: haxe.io.Float32Array = null; + static var multiBlend: Float = 1.0; // 0 = old data, 1 = new data + static var multiBlendStartTime: Float = 0.0; + + static var recomputeStartTime: Float = 0.0; + static var singlePending: { density: Vec3, data: haxe.io.Float32Array, row: Int } = null; + static var multiPending: { + sunElevation: Float, sunSize: Float, altitude: Float, sunDisc: Bool, density: Vec3, + data: haxe.io.Float32Array, row: Int, transRow: Int, phase: Int, + sdx: Float, sdy: Float, sdz: Float, roz: Float + } = null; + + // incremental sky recompute where larger values do less work per frame + public static var recomputTime: Float = 0.3; + + public static function recomputeSingleScatter(world: WorldData) { + if (!registered) register(); if (world == null || world.raw.sky_density == null) return; if (singleScatterData == null) singleScatterData = new SingleScatteringData(); var density = world.raw.sky_density; - singleScatterData.computeLUT(new Vec3(density[0], density[1], density[2])); + var densityVec = new Vec3(density[0], density[1], density[2]); + + if (!singleScatterData.needsRecompute(densityVec)) return; + + if (singlePending != null) return; + + if (singleScatterData.lut == null) { + singleScatterData.computeLUT(densityVec); + } else { + recomputeStartTime = haxe.Timer.stamp(); + singlePending = { + density: densityVec, + data: new haxe.io.Float32Array(SingleScatteringData.lutHeightSteps * SingleScatteringData.lutAngleSteps * 4), + row: 0 + }; + } } - /** Sets the sky's density parameters and calls `recomputeSingleScatter()` afterwards. **/ + // Sets the sky's density parameters and calls `recomputeSingleScatter()` afterwards. public static function setDensity(world: WorldData, densityAir: FastFloat, densityDust: FastFloat, densityOzone: FastFloat) { if (world == null) return; @@ -46,6 +84,7 @@ class Sky { } public static function recomputeMultiScatter(world: WorldData) { + if (!registered) register(); if (world == null) return; var raw = world.raw; if (raw.sky_density == null) return; @@ -53,6 +92,7 @@ class Sky { if (multiScatterData == null) multiScatterData = new MultipleScatteringData(); var density = raw.sky_density; + var densityVec = new Vec3(density[0], density[1], density[2]); var sunElevation = raw.sky_sun_elevation != null ? raw.sky_sun_elevation : 0.0; var sunRotation = raw.sky_sun_rotation != null ? raw.sky_sun_rotation : 0.0; var sunSize = raw.sky_sun_size != null ? raw.sky_sun_size : 0.545; @@ -60,15 +100,32 @@ class Sky { var altitude = raw.sky_altitude != null ? raw.sky_altitude : 0.0; var sunDisc = raw.sky_sun_disc != null ? raw.sky_sun_disc : 1; - multiScatterData.compute( - sunElevation, - sunRotation, - sunSize, - sunIntensity, - altitude, - sunDisc != 0, - new Vec3(density[0], density[1], density[2]) - ); + if (!multiScatterData.needsRecompute(densityVec, sunElevation, altitude, sunSize, sunDisc != 0)) return; + + if (multiPending != null) return; + + if (multiScatterData.lut == null) { + multiScatterData.compute(sunElevation, sunRotation, sunSize, sunIntensity, + altitude, sunDisc != 0, densityVec); + } else { + recomputeStartTime = haxe.Timer.stamp(); + multiScatterData.setDensities(densityVec); + var needTrans = multiScatterData.needsTransmittance(densityVec); + var altKm = Math.max(1.0, Math.min(99999.0, altitude)) / 1000.0; + var sunZenithCosAngle = Math.cos(Math.PI / 2.0 - sunElevation); + var sdy = Math.sqrt(1.0 - sunZenithCosAngle * sunZenithCosAngle); + var sdz = sunZenithCosAngle; + var sdx = 0.0; + var roz = MultipleScatteringData.EARTH_RADIUS + altKm; + + multiPending = { + sunElevation: sunElevation, sunSize: sunSize, altitude: altitude, + sunDisc: sunDisc != 0, density: densityVec, + data: new haxe.io.Float32Array(MultipleScatteringData.lutWidth * MultipleScatteringData.lutHeight * 4), + row: 0, transRow: 0, phase: needTrans ? 0 : 1, + sdx: sdx, sdy: sdy, sdz: sdz, roz: roz + }; + } } public static function setMultiScatterParams(world: WorldData, sunElevation: FastFloat, sunRotation: FastFloat, @@ -91,6 +148,94 @@ class Sky { recomputeMultiScatter(world); } + + // CPU LUT blending for smooth transition during expensive recomputation + + public static function update() { + if (multiBlend < 1.0 && multiPrevData != null && multiNewData != null) { + multiBlend = Math.min(1.0, (haxe.Timer.stamp() - multiBlendStartTime) / recomputTime); + var lutSize = MultipleScatteringData.lutWidth * MultipleScatteringData.lutHeight * 4; + if (multiBlendData == null || multiBlendData.length != lutSize) { + multiBlendData = new haxe.io.Float32Array(lutSize); + } + for (i in 0...lutSize) { + multiBlendData[i] = multiPrevData[i] * (1.0 - multiBlend) + multiNewData[i] * multiBlend; + } + multiScatterData.applyLUTData(multiBlendData); + if (multiBlend >= 1.0) { + multiPrevData = null; + multiNewData = null; + multiBlendData = null; + } + } + + if (singlePending == null && multiPending == null) return; + var elapsed = haxe.Timer.stamp() - recomputeStartTime; + var fraction = Math.min(1.0, elapsed / recomputTime); + + if (singlePending != null) { + var targetRow = Math.ceil(SingleScatteringData.lutHeightSteps * fraction); + while (singlePending.row < targetRow) { + singleScatterData.computeLUTRow(singlePending.density, singlePending.row, singlePending.data); + singlePending.row++; + } + if (singlePending.row >= SingleScatteringData.lutHeightSteps) { + singleScatterData.applyLUTData(singlePending.data); + singleScatterData.setLastDensity(singlePending.density); + singlePending = null; + } + } + + if (multiPending != null) { + if (multiPending.phase == 0) { + var targetTrans = Math.ceil(MultipleScatteringData.transmittanceResY * fraction); + while (multiPending.transRow < targetTrans) { + multiScatterData.precomputeTransmittanceRow(multiPending.transRow); + multiPending.transRow++; + } + if (multiPending.transRow >= MultipleScatteringData.transmittanceResY) { + multiScatterData.finishTransmittance(multiPending.density); + multiPending.phase = 1; + } + } + if (multiPending.phase == 1) { + var targetRow = Math.ceil(MultipleScatteringData.lutHeight * fraction); + while (multiPending.row < targetRow) { + multiScatterData.computeLUTRow( + multiPending.sdx, multiPending.sdy, multiPending.sdz, multiPending.roz, + multiPending.row, multiPending.data + ); + multiPending.row++; + } + if (multiPending.row >= MultipleScatteringData.lutHeight) { + // old LUT data for blending and use current blended data if a blend is in progress + if (multiBlend < 1.0 && multiBlendData != null) { + multiPrevData = multiBlendData; + } else { + multiPrevData = multiScatterData.lastLUTData; + } + multiNewData = multiPending.data; + multiScatterData.lastLUTData = multiPending.data; + multiScatterData.lastElevation = multiPending.sunElevation; + multiScatterData.lastAltitude = multiPending.altitude; + multiScatterData.lastSunSize = multiPending.sunSize; + multiScatterData.lastSunDisc = multiPending.sunDisc; + multiScatterData.setSunData( + multiPending.sunElevation, multiPending.sunSize, + multiPending.altitude, multiPending.sunDisc + ); + // blend if we have old data, otherwise apply directly + if (multiPrevData != null) { + multiBlend = 0.0; + multiBlendStartTime = haxe.Timer.stamp(); + } else { + multiScatterData.applyLUTData(multiPending.data); + } + multiPending = null; + } + } + } + } } /** @@ -104,6 +249,7 @@ class SingleScatteringData { public var lut: kha.Image; + var lastDensity: Vec3 = null; /** The amount of individual sample heights stored in the LUT (and the width of the LUT image). @@ -147,135 +293,126 @@ class SingleScatteringData { // to include all 21 wavelengths gave unrealistic results. public static var ozoneCoeff = new Vec3(1.59051840791988e-6, 0.00000096707041180970, 0.00000007309568762914); + // salar versions of coefficients + static inline var rx = 5.5e-6; + static inline var ry = 13.0e-6; + static inline var rz = 22.4e-6; + static inline var ox = 1.59051840791988e-6; + static inline var oy = 0.00000096707041180970; + static inline var oz = 0.00000007309568762914; + public function new() {} /** Approximates the density of ozone for a given sample height. **/ - function getOzoneDensity(height: FastFloat): FastFloat { + + static inline function getOzoneDensity(height: FastFloat): FastFloat { // Values are taken from Cycles code - if (height < 10000.0 || height >= 40000.0) { - return 0.0; - } - if (height < 25000.0) { - return (height - 10000.0) / 15000.0; - } + if (height < 10000.0 || height >= 40000.0) return 0.0; + if (height < 25000.0) return (height - 10000.0) / 15000.0; return -((height - 40000.0) / 15000.0); } /** - 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. + Computes the LUT data for the given density values (thread-safe). + Uses scalar math — zero allocations in the inner loop. + @param density 3D vector of air density, dust density, ozone density + @return Float32Array of RGBA128 pixel data **/ - function raySphereIntersection(rayOrigin: Vec3, rayDirection: Vec3, sphereRadius: Float): 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; - - // 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) - ); + public function computeLUTData(density: Vec3): haxe.io.Float32Array { + var imageData = new haxe.io.Float32Array(lutHeightSteps * lutAngleSteps * 4); + for (x in 0...lutHeightSteps) { + computeLUTRow(density, x, imageData); + } + return imageData; } - /** - 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); + // computes a single height row of the LUT for incremental recompute + public function computeLUTRow(density: Vec3, x: Int, imageData: haxe.io.Float32Array) { + var radiusPlanetMeters = radiusPlanet * 1000; + var dAir = density.x; + var dDust = density.y; + var dOzone = density.z; + var invRayleigh = 1.0 / rayleighScale; + var invMie = 1.0 / mieScale; + var invJ = 1.0 / jSteps; - for (x in 0...lutHeightSteps) { - var height = (x / (lutHeightSteps - 1)); + // Use quadratic height for better horizon precision + var height = (x / (lutHeightSteps - 1)); + height *= height; + height *= radiusAtmo * 1000; // Denormalize height + var iz = height + radiusPlanetMeters; - // Use quadratic height for better horizon precision - height *= height; - height *= radiusAtmo * 1000; // Denormalize height + for (y in 0...lutAngleSteps) { + var sunTheta = y / (lutAngleSteps - 1) * 2 - 1; - for (y in 0...lutAngleSteps) { - var sunTheta = y / (lutAngleSteps - 1) * 2 - 1; + // 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 - // 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 sy = Math.sin(sunTheta); + var sz = Math.cos(sunTheta); - 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 + var ozKm = iz * 0.001; + var b = 2.0 * sz * ozKm; + var c = ozKm * ozKm - radiusAtmo * radiusAtmo; + var disc = b * b - 4.0 * c; + var jStepSize: FastFloat; + if (disc < 0.0) { + jStepSize = 1e5 * invJ * 1000; + } else { + var sqd = Math.sqrt(disc); + jStepSize = ((-b + sqd) * 0.5) / jSteps * 1000; } - } + var odRay = 0.0; + var odMie = 0.0; + var odOzone = 0.0; + var jTime = 0.0; + + for (i in 0...jSteps) { + var t = jTime + jStepSize * 0.5; + var jy = sy * t; + var jz = iz + sz * t; + var jHeight = Math.sqrt(jy * jy + jz * jz) - radiusPlanetMeters; + + odRay += Math.exp(-jHeight * invRayleigh) * dAir; + odMie += Math.exp(-jHeight * invMie) * dDust; + odOzone += getOzoneDensity(jHeight) * dOzone; + + jTime += jStepSize; + } + + odRay *= jStepSize; + odMie *= jStepSize; + odOzone *= jStepSize; + + var mie = mieCoeff * odMie; + var pixelIndex = (x + y * lutHeightSteps) * 4; + imageData[pixelIndex + 0] = Math.exp(-(rx * odRay + mie + ox * odOzone)); + imageData[pixelIndex + 1] = Math.exp(-(ry * odRay + mie + oy * odOzone)); + imageData[pixelIndex + 2] = Math.exp(-(rz * odRay + mie + oz * odOzone)); + imageData[pixelIndex + 3] = 1.0; // Unused + } + } + + // creates the LUT image from precomputed data called on main + public function applyLUTData(imageData: haxe.io.Float32Array) { 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 { - var radiusPlanetMeters = radiusPlanet * 1000; + public function needsRecompute(density: Vec3): Bool { + return lut == null || lastDensity == null || + lastDensity.x != density.x || lastDensity.y != density.y || lastDensity.z != density.z; + } - // Reconstruct values from the shader - var iPos = new Vec3(0, 0, height + radiusPlanetMeters); - var pSun = new Vec3(0.0, Math.sin(sunTheta), Math.cos(sunTheta)).normalize(); + public function setLastDensity(density: Vec3) { + lastDensity = new Vec3(density.x, density.y, density.z); + } - var jTime: FastFloat = 0.0; - // We compute the ray-sphere intersection in km to allow larger - // atmosphere radii (radius is squared inside raySphereIntersection()) - var jStepSize: FastFloat = raySphereIntersection(iPos.clone().mult(0.001), pSun, radiusAtmo).y / jSteps; - jStepSize *= 1000; // convert back to m - - // 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() - radiusPlanetMeters; - - // 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; - } - - jODepth.mult(jStepSize); - - // Precalculate a part of the secondary attenuation. - // For one variable (e.g. x) in the vector, the formula is as follows: - // - // attn.x = exp(-(coeffX * (firstOpticalDepth.x + secondOpticalDepth.x))) - // - // We can split that up via: - // - // attn.x = exp(-(coeffX * firstOpticalDepth.x + coeffX * secondOpticalDepth.x)) - // = exp(-(coeffX * firstOpticalDepth.x)) * exp(-(coeffX * secondOpticalDepth.x)) - // - // The first factor of the resulting multiplication is calculated in the - // shader, but we can already precalculate the second one. As a side - // effect this keeps the range of the LUT values small because we don't - // store the optical depth but the attenuation. - var jAttenuation = new Vec3(); - var mie = mieCoeff * jODepth.y; - jAttenuation.addf(mie, mie, mie); - jAttenuation.add(rayleighCoeff.clone().mult(jODepth.x)); - jAttenuation.add(ozoneCoeff.clone().mult(jODepth.z)); - jAttenuation.exp(jAttenuation.mult(-1)); - - return jAttenuation; + public function computeLUT(density: Vec3) { + applyLUTData(computeLUTData(density)); + lastDensity = new Vec3(density.x, density.y, density.z); } } @@ -292,18 +429,26 @@ class MultipleScatteringData { public var sunTop: Vec3; public var earthIntersectionAngle: FastFloat; - public static var lutWidth = 256; - public static var lutHeight = 128; + public static var lutWidth = 128; + public static var lutHeight = 64; - static var transmittanceResX = 256; - static var transmittanceResY = 64; + static var transmittanceResX = 128; + public static var transmittanceResY = 32; var transmittanceLUT: haxe.io.Float32Array; + var transmittanceValid = false; + var transmittanceDensity: Vec3 = null; - static var transmittanceSteps = 64; - static var inScatteringSteps = 64; + public var lastElevation: Float = 0; + public var lastAltitude: Float = 0; + public var lastSunSize: Float = 0; + public var lastSunDisc: Bool = false; + public var lastLUTData: haxe.io.Float32Array = null; + + static var transmittanceSteps = 32; + static var inScatteringSteps = 32; // Atmosphere constants sky_multiple_scattering.cpp - static var EARTH_RADIUS = 6371.0; + public static var EARTH_RADIUS = 6371.0; static var ATMOSPHERE_THICKNESS = 100.0; static var ATMOSPHERE_RADIUS = 6471.0; // EARTH_RADIUS + ATMOSPHERE_THICKNESS @@ -336,6 +481,14 @@ class MultipleScatteringData { var aerosolDensity: Float = 1.0; var ozoneDensity: Float = 1.0; + // avoid allocations in hot loops + static var scratch0: Array = [0, 0, 0, 0]; + static var scratch1: Array = [0, 0, 0, 0]; + static var scratch2: Array = [0, 0, 0, 0]; + static var scratch3: Array = [0, 0, 0, 0]; + static var scratch4: Array = [0, 0, 0, 0]; + static var scratch5: Array = [0, 0, 0, 0]; + public function new() { sunBottom = new Vec3(); sunTop = new Vec3(); @@ -343,81 +496,58 @@ class MultipleScatteringData { transmittanceLUT = new haxe.io.Float32Array(transmittanceResY * transmittanceResX * 4); } - static inline function exp4(a: Array): Array { - return [Math.exp(a[0]), Math.exp(a[1]), Math.exp(a[2]), Math.exp(a[3])]; + static inline function add4(a: Array, b: Array, out: Array) { + out[0] = a[0] + b[0]; out[1] = a[1] + b[1]; out[2] = a[2] + b[2]; out[3] = a[3] + b[3]; } - - static inline function scale4(a: Array, s: Float): Array { - return [a[0] * s, a[1] * s, a[2] * s, a[3] * s]; + static inline function scale4(a: Array, s: Float, out: Array) { + out[0] = a[0] * s; out[1] = a[1] * s; out[2] = a[2] * s; out[3] = a[3] * s; } - - static inline function add4(a: Array, b: Array): Array { - return [a[0] + b[0], a[1] + b[1], a[2] + b[2], a[3] + b[3]]; - } - - static inline function mult4(a: Array, b: Array): Array { - return [a[0] * b[0], a[1] * b[1], a[2] * b[2], a[3] * b[3]]; - } - - static inline function div4(a: Array, b: Array): Array { - return [a[0] / b[0], a[1] / b[1], a[2] / b[2], a[3] / b[3]]; - } - - static inline function max4(a: Array, b: Float): Array { - return [Math.max(a[0], b), Math.max(a[1], b), Math.max(a[2], b), Math.max(a[3], b)]; - } - - static inline function mix4(a: Array, b: Array, t: Float): Array { - return [a[0] + t * (b[0] - a[0]), a[1] + t * (b[1] - a[1]), a[2] + t * (b[2] - a[2]), a[3] + t * (b[3] - a[3])]; - } - // Atmospheric functions - function getMolecularScatteringCoefficient(h: Float): Array { - // MOLECULAR_SCATTERING_COEFFICIENT_BASE * exp(-0.07771971 * h^1.16364243) + function getMolecularScatteringCoefficient(h: Float, out: Array) { var f = Math.exp(-0.07771971 * Math.pow(h, 1.16364243)); - return scale4(MOLECULAR_SCATTERING_COEFFICIENT_BASE, f); + out[0] = MOLECULAR_SCATTERING_COEFFICIENT_BASE[0] * f; + out[1] = MOLECULAR_SCATTERING_COEFFICIENT_BASE[1] * f; + out[2] = MOLECULAR_SCATTERING_COEFFICIENT_BASE[2] * f; + out[3] = MOLECULAR_SCATTERING_COEFFICIENT_BASE[3] * f; } - function getMolecularAbsorptionCoefficient(h: Float): Array { + function getMolecularAbsorptionCoefficient(h: Float, out: Array) { var logH = Math.log(Math.max(h, 1e-4)); var density = 3.78547397e20 * Math.exp(-(logH - 3.22261) * (logH - 3.22261) * 5.55555555 - logH); - var result: Array = []; - for (i in 0...4) { - result[i] = OZONE_ABSORPTION_CROSS_SECTION[i] * OZONE_MEAN_DOBSON * density; - } - return result; + out[0] = OZONE_ABSORPTION_CROSS_SECTION[0] * OZONE_MEAN_DOBSON * density; + out[1] = OZONE_ABSORPTION_CROSS_SECTION[1] * OZONE_MEAN_DOBSON * density; + out[2] = OZONE_ABSORPTION_CROSS_SECTION[2] * OZONE_MEAN_DOBSON * density; + out[3] = OZONE_ABSORPTION_CROSS_SECTION[3] * OZONE_MEAN_DOBSON * density; } - function getAerosolDensity(h: Float): Float { - var division = AEROSOL_BACKGROUND_DENSITY / AEROSOL_BASE_DENSITY; - return AEROSOL_BASE_DENSITY * (Math.exp(-h / AEROSOL_HEIGHT_SCALE) + division); + static inline function getAerosolDensity(h: Float): Float { + return AEROSOL_BASE_DENSITY * (Math.exp(-h / AEROSOL_HEIGHT_SCALE) + AEROSOL_BACKGROUND_DENSITY / AEROSOL_BASE_DENSITY); } - function getAtmosphereCollisionCoefficients(h: Float): { - aerosolAbs: Array, aerosolSca: Array, molAbs: Array, molSca: Array - } { + function getAtmosphereCollisionCoefficients(h: Float) { var localAerosolDensity = getAerosolDensity(h) * aerosolDensity; - var aerosolAbs = scale4(AEROSOL_ABSORPTION_CROSS_SECTION, localAerosolDensity); - var aerosolSca = scale4(AEROSOL_SCATTERING_CROSS_SECTION, localAerosolDensity); - var molAbs = scale4(getMolecularAbsorptionCoefficient(h), ozoneDensity); - var molSca = scale4(getMolecularScatteringCoefficient(h), airDensity); - return { aerosolAbs: aerosolAbs, aerosolSca: aerosolSca, molAbs: molAbs, molSca: molSca }; + scale4(AEROSOL_ABSORPTION_CROSS_SECTION, localAerosolDensity, scratch0); + scale4(AEROSOL_SCATTERING_CROSS_SECTION, localAerosolDensity, scratch1); + getMolecularAbsorptionCoefficient(h, scratch2); + scale4(scratch2, ozoneDensity, scratch2); + getMolecularScatteringCoefficient(h, scratch3); + scale4(scratch3, airDensity, scratch3); } - function molecularPhaseFunction(cosTheta: Float): Float { + static inline function molecularPhaseFunction(cosTheta: Float): Float { return RAYLEIGH_PHASE_SCALE * (1.0 + cosTheta * cosTheta); } - function aerosolPhaseFunction(cosTheta: Float): Float { + static inline function aerosolPhaseFunction(cosTheta: Float): Float { var den = 1.0 + SQR_G + 2.0 * G * cosTheta; return PHASE_ISOTROPIC * (1.0 - SQR_G) / (den * Math.sqrt(den)); } - // nearest positive intersection distance, or -1 if no intersection - function raySphereIntersection(pos: Vec3, dir: Vec3, radius: Float): Float { - var b = pos.dot(dir); - var c = pos.dot(pos) - radius * radius; + static inline function raySphereIntersectScalar(px: Float, py: Float, pz: Float, + dx: Float, dy: Float, dz: Float, radius: Float): Float { + var b = px * dx + py * dy + pz * dz; + var c = px * px + py * py + pz * pz - radius * radius; if (c > 0.0 && b > 0.0) return -1.0; var d = b * b - c; if (d < 0.0) return -1.0; @@ -425,60 +555,98 @@ class MultipleScatteringData { return -b - Math.sqrt(d); } - function sunDirection(cosTheta: Float): Vec3 { - // Place sun at azimuth=0 in the atan(dir.x, dir.y) convention - // Blender uses (-sqrt(1-cos²θ), 0, cosθ) but our azimuth convention - // requires (0, sqrt(1-cos²θ), cosθ) for the sun to be at azimuth=0 - return new Vec3(0.0, Math.sqrt(1.0 - cosTheta * cosTheta), cosTheta); + static inline function spectralToXYZScalars(l0: Float, l1: Float, l2: Float, l3: Float): Vec3 { + return new Vec3( + SPECTRAL_XYZ[0][0] * l0 + SPECTRAL_XYZ[1][0] * l1 + SPECTRAL_XYZ[2][0] * l2 + SPECTRAL_XYZ[3][0] * l3, + SPECTRAL_XYZ[0][1] * l0 + SPECTRAL_XYZ[1][1] * l1 + SPECTRAL_XYZ[2][1] * l2 + SPECTRAL_XYZ[3][1] * l3, + SPECTRAL_XYZ[0][2] * l0 + SPECTRAL_XYZ[1][2] * l1 + SPECTRAL_XYZ[2][2] * l2 + SPECTRAL_XYZ[3][2] * l3 + ); } - - function spectralToXYZ(L: Array): Vec3 { - var x = 0.0, y = 0.0, z = 0.0; - for (i in 0...4) { - x += SPECTRAL_XYZ[i][0] * L[i]; - y += SPECTRAL_XYZ[i][1] * L[i]; - z += SPECTRAL_XYZ[i][2] * L[i]; - } - return new Vec3(x, y, z); + // XYZ directly into imageData avoiding Vec3 allocation in loops + static inline function spectralToXYZ(l0: Float, l1: Float, l2: Float, l3: Float, + imageData: haxe.io.Float32Array, pixelIndex: Int) { + imageData[pixelIndex + 0] = SPECTRAL_XYZ[0][0] * l0 + SPECTRAL_XYZ[1][0] * l1 + SPECTRAL_XYZ[2][0] * l2 + SPECTRAL_XYZ[3][0] * l3; + imageData[pixelIndex + 1] = SPECTRAL_XYZ[0][1] * l0 + SPECTRAL_XYZ[1][1] * l1 + SPECTRAL_XYZ[2][1] * l2 + SPECTRAL_XYZ[3][1] * l3; + imageData[pixelIndex + 2] = SPECTRAL_XYZ[0][2] * l0 + SPECTRAL_XYZ[1][2] * l1 + SPECTRAL_XYZ[2][2] * l2 + SPECTRAL_XYZ[3][2] * l3; + imageData[pixelIndex + 3] = 1.0; } - - function getTransmittance(cosTheta: Float, normalizedAltitude: Float): Array { - var sunDir = sunDirection(cosTheta); + function getTransmittance(cosTheta: Float, normalizedAltitude: Float, out: Array) { + var sy = Math.sqrt(1.0 - cosTheta * cosTheta); + var sz = cosTheta; var distToCenter = EARTH_RADIUS + (ATMOSPHERE_RADIUS - EARTH_RADIUS) * normalizedAltitude; - var rayOrigin = new Vec3(0.0, 0.0, distToCenter); - var tD = raySphereIntersection(rayOrigin, sunDir, ATMOSPHERE_RADIUS); + var tD = raySphereIntersectScalar(0, 0, distToCenter, 0, sy, sz, ATMOSPHERE_RADIUS); var tStep = tD / transmittanceSteps; - var result: Array = [0.0, 0.0, 0.0, 0.0]; + // Accumulate in local scalars (out may alias scratch0 which gets clobbered) + var acc0 = 0.0, acc1 = 0.0, acc2 = 0.0, acc3 = 0.0; for (step in 0...transmittanceSteps) { var t = (step + 0.5) * tStep; - var xT = rayOrigin.clone().add(sunDir.clone().mult(t)); - var altitude = Math.max(xT.length() - EARTH_RADIUS, 0.0); - var coeffs = getAtmosphereCollisionCoefficients(altitude); - var extinction = add4(add4(coeffs.aerosolAbs, coeffs.aerosolSca), add4(coeffs.molAbs, coeffs.molSca)); - result = add4(result, scale4(extinction, tStep)); + var xtY = sy * t; + var xtZ = distToCenter + sz * t; + var altitude = Math.max(Math.sqrt(xtY * xtY + xtZ * xtZ) - EARTH_RADIUS, 0.0); + getAtmosphereCollisionCoefficients(altitude); + // extinction = aerosolAbs + aerosolSca + molAbs + molSca + add4(scratch0, scratch1, scratch4); + add4(scratch2, scratch3, scratch5); + add4(scratch4, scratch5, scratch4); + // result += extinction * tStep + acc0 += scratch4[0] * tStep; + acc1 += scratch4[1] * tStep; + acc2 += scratch4[2] * tStep; + acc3 += scratch4[3] * tStep; } - return exp4(scale4(result, -1.0)); + out[0] = Math.exp(-acc0); + out[1] = Math.exp(-acc1); + out[2] = Math.exp(-acc2); + out[3] = Math.exp(-acc3); } function precomputeTransmittanceLUT() { + var d = new Vec3(airDensity, aerosolDensity, ozoneDensity); + if (!needsTransmittance(d)) return; for (y in 0...transmittanceResY) { - for (x in 0...transmittanceResX) { - var u = x / (transmittanceResX - 1); - var v = y / (transmittanceResY - 1); - var t = getTransmittance(u * 2.0 - 1.0, v); - var idx = (y * transmittanceResX + x) * 4; - transmittanceLUT[idx] = t[0]; - transmittanceLUT[idx + 1] = t[1]; - transmittanceLUT[idx + 2] = t[2]; - transmittanceLUT[idx + 3] = t[3]; - } + precomputeTransmittanceRow(y); + } + finishTransmittance(d); + } + + public function needsTransmittance(density: Vec3): Bool { + return !(transmittanceValid && transmittanceDensity != null && + transmittanceDensity.x == density.x && + transmittanceDensity.y == density.y && + transmittanceDensity.z == density.z); + } + + public function setDensities(density: Vec3) { + airDensity = density.x; + aerosolDensity = density.y; + ozoneDensity = density.z; + } + + public function precomputeTransmittanceRow(y: Int) { + for (x in 0...transmittanceResX) { + var u = x / (transmittanceResX - 1); + var v = y / (transmittanceResY - 1); + getTransmittance(u * 2.0 - 1.0, v, scratch0); + var idx = (y * transmittanceResX + x) * 4; + transmittanceLUT[idx] = scratch0[0]; + transmittanceLUT[idx + 1] = scratch0[1]; + transmittanceLUT[idx + 2] = scratch0[2]; + transmittanceLUT[idx + 3] = scratch0[3]; } } - function lookupTransmittance(cosTheta: Float, normalizedAltitude: Float): Array { + public function finishTransmittance(density: Vec3) { + transmittanceValid = true; + if (transmittanceDensity == null) transmittanceDensity = new Vec3(); + transmittanceDensity.x = density.x; + transmittanceDensity.y = density.y; + transmittanceDensity.z = density.z; + } + + function lookupTransmittance(cosTheta: Float, normalizedAltitude: Float, out: Array) { var u = Math.max(0.0, Math.min(1.0, cosTheta * 0.5 + 0.5)); var v = Math.max(0.0, Math.min(1.0, normalizedAltitude)); var fx = u * (transmittanceResX - 1); @@ -490,118 +658,167 @@ class MultipleScatteringData { var dxF = fx - x1; var dyF = fy - y1; - function getPixel(px: Int, py: Int): Array { - var idx = (py * transmittanceResX + px) * 4; - return [transmittanceLUT[idx], transmittanceLUT[idx + 1], transmittanceLUT[idx + 2], transmittanceLUT[idx + 3]]; - } + var idx00 = (y1 * transmittanceResX + x1) * 4; + var idx01 = (y1 * transmittanceResX + x2) * 4; + var idx10 = (y2 * transmittanceResX + x1) * 4; + var idx11 = (y2 * transmittanceResX + x2) * 4; - var bottom = mix4(getPixel(x1, y1), getPixel(x2, y1), dxF); - var top = mix4(getPixel(x1, y2), getPixel(x2, y2), dxF); - return mix4(bottom, top, dyF); + // Bilinear interpolation directly into out — no scratch arrays needed + for (w in 0...4) { + var b = transmittanceLUT[idx00 + w] + dxF * (transmittanceLUT[idx01 + w] - transmittanceLUT[idx00 + w]); + var t = transmittanceLUT[idx10 + w] + dxF * (transmittanceLUT[idx11 + w] - transmittanceLUT[idx10 + w]); + out[w] = b + dyF * (t - b); + } } - function lookupTransmittanceAtGround(cosTheta: Float): Array { + function lookupTransmittanceAtGround(cosTheta: Float, out: Array) { var u = Math.max(0.0, Math.min(1.0, cosTheta * 0.5 + 0.5)); var fx = u * (transmittanceResX - 1); var x1 = Std.int(fx); var x2 = Std.int(Math.min(x1 + 1, transmittanceResX - 1)); var dxF = fx - x1; - var y = 0; - function getPixel(px: Int): Array { - var idx = (y * transmittanceResX + px) * 4; - return [transmittanceLUT[idx], transmittanceLUT[idx + 1], transmittanceLUT[idx + 2], transmittanceLUT[idx + 3]]; + var idx0 = (0 * transmittanceResX + x1) * 4; + var idx1 = (0 * transmittanceResX + x2) * 4; + for (w in 0...4) { + out[w] = transmittanceLUT[idx0 + w] + dxF * (transmittanceLUT[idx1 + w] - transmittanceLUT[idx0 + w]); } - return mix4(getPixel(x1), getPixel(x2), dxF); } - function lookupTransmittanceToSun(normalizedAltitude: Float): Array { + function lookupTransmittanceToSun(normalizedAltitude: Float, out: Array) { var v = Math.max(0.0, Math.min(1.0, normalizedAltitude)); var fy = v * (transmittanceResY - 1); var y1 = Std.int(fy); var y2 = Std.int(Math.min(y1 + 1, transmittanceResY - 1)); var dyF = fy - y1; var x = transmittanceResX - 1; - function getPixel(py: Int): Array { - var idx = (py * transmittanceResX + x) * 4; - return [transmittanceLUT[idx], transmittanceLUT[idx + 1], transmittanceLUT[idx + 2], transmittanceLUT[idx + 3]]; + var idx0 = (y1 * transmittanceResX + x) * 4; + var idx1 = (y2 * transmittanceResX + x) * 4; + for (w in 0...4) { + out[w] = transmittanceLUT[idx0 + w] + dyF * (transmittanceLUT[idx1 + w] - transmittanceLUT[idx0 + w]); } - return mix4(getPixel(y1), getPixel(y2), dyF); } - // approximation - - function lookupMultiscattering(cosTheta: Float, normalizedHeight: Float, d: Float): Array { + function lookupMultiscattering(cosTheta: Float, normalizedHeight: Float, d: Float, out: Array) { var omega = 2.0 * Math.PI * (1.0 - Math.sqrt(Math.max(0.0, 1.0 - (EARTH_RADIUS / d) * (EARTH_RADIUS / d)))); - var tToGround = lookupTransmittanceAtGround(cosTheta); - var tGroundToSample = div4(lookupTransmittanceToSun(0.0), lookupTransmittanceToSun(normalizedHeight)); - // 2nd order scattering from the ground - var lGround = scale4(mult4(mult4(mult4(scale4(GROUND_ALBEDO, 1.0 / Math.PI), tToGround), tGroundToSample), [cosTheta, cosTheta, cosTheta, cosTheta]), PHASE_ISOTROPIC * omega); - // Fit of Earth's multiple scattering from other points in the atmosphere - var msFactor = 1.0 / (1.0 + 5.0 * Math.exp(-17.92 * cosTheta)); - var lMs = scale4([0.217, 0.347, 0.594, 1.0], 0.02 * msFactor); - return add4(lMs, lGround); + lookupTransmittanceAtGround(cosTheta, scratch0); // tToGround + lookupTransmittanceToSun(0.0, scratch1); // tSun0 + lookupTransmittanceToSun(normalizedHeight, scratch2); // tSunH + // tGroundToSample = tSun0 / tSunH + scratch3[0] = scratch1[0] / scratch2[0]; scratch3[1] = scratch1[1] / scratch2[1]; + scratch3[2] = scratch1[2] / scratch2[2]; scratch3[3] = scratch1[3] / scratch2[3]; + // lGround = GROUND_ALBEDO/PI * tToGround * tGroundToSample * cosTheta * PHASE_ISOTROPIC * omega + var groundScale = (1.0 / Math.PI) * cosTheta * PHASE_ISOTROPIC * omega; + for (w in 0...4) { + out[w] = GROUND_ALBEDO[w] * scratch0[w] * scratch3[w] * groundScale; + } + // lMs = [0.217, 0.347, 0.594, 1.0] * 0.02 * msFactor + var msFactor = 0.02 / (1.0 + 5.0 * Math.exp(-17.92 * cosTheta)); + out[0] += 0.217 * msFactor; + out[1] += 0.347 * msFactor; + out[2] += 0.594 * msFactor; + out[3] += 1.0 * msFactor; } - function getInscattering(sunDir: Vec3, rayOrigin: Vec3, rayDir: Vec3, tD: Float): Array { - var cosTheta = rayDir.clone().mult(-1.0).dot(sunDir); + // Scalar inscattering — sunDir and rayDir as scalars, uses scratch arrays + function getInscatteringScalars( + sdx: Float, sdy: Float, sdz: Float, // sun dir + ox: Float, oy: Float, oz: Float, // ray origin + dx: Float, dy: Float, dz: Float, // ray dir (normalized) + tD: Float, out: Array) { + + // cosTheta = -rayDir . sunDir + var cosTheta = -(dx * sdx + dy * sdy + dz * sdz); var molPhase = molecularPhaseFunction(cosTheta); var aerPhase = aerosolPhaseFunction(cosTheta); var dt = tD / inScatteringSteps; - var lInscattering: Array = [0.0, 0.0, 0.0, 0.0]; - var transmittance: Array = [1.0, 1.0, 1.0, 1.0]; + + out[0] = 0; out[1] = 0; out[2] = 0; out[3] = 0; + // transmittance = [1,1,1,1] + var tr0 = 1.0, tr1 = 1.0, tr2 = 1.0, tr3 = 1.0; + // Accumulate inscattering in local scalars (out may alias scratch arrays) + var li0 = 0.0, li1 = 0.0, li2 = 0.0, li3 = 0.0; for (i in 0...inScatteringSteps) { var t = (i + 0.5) * dt; - var xT = rayOrigin.clone().add(rayDir.clone().mult(t)); - var distToCenter = xT.length(); - var zenithDir = xT.clone().mult(1.0 / distToCenter); + // xT = rayOrigin + rayDir * t + var xtx = ox + dx * t; + var xty = oy + dy * t; + var xtz = oz + dz * t; + var distToCenter = Math.sqrt(xtx * xtx + xty * xty + xtz * xtz); var altitude = Math.max(distToCenter - EARTH_RADIUS, 0.0); var normalizedAltitude = altitude / ATMOSPHERE_THICKNESS; - var sampleCosTheta = zenithDir.dot(sunDir); + // zenithDir = xT / distToCenter + var invDist = 1.0 / distToCenter; + var zdx = xtx * invDist, zdy = xty * invDist, zdz = xtz * invDist; + var sampleCosTheta = zdx * sdx + zdy * sdy + zdz * sdz; - var coeffs = getAtmosphereCollisionCoefficients(altitude); - var extinction = add4(add4(coeffs.aerosolAbs, coeffs.aerosolSca), add4(coeffs.molAbs, coeffs.molSca)); - var tToSun = lookupTransmittance(sampleCosTheta, normalizedAltitude); - var ms = lookupMultiscattering(sampleCosTheta, normalizedAltitude, distToCenter); + getAtmosphereCollisionCoefficients(altitude); + // extinction = aerosolAbs + aerosolSca + molAbs + molSca + add4(scratch0, scratch1, scratch4); + add4(scratch2, scratch3, scratch5); + add4(scratch4, scratch5, scratch4); // extinction in scratch4 + + // Save molSca (scratch3) and aerSca (scratch1) before lookups clobber scratch arrays + var molSca0 = scratch3[0], molSca1 = scratch3[1], molSca2 = scratch3[2], molSca3 = scratch3[3]; + var aerSca0 = scratch1[0], aerSca1 = scratch1[1], aerSca2 = scratch1[2], aerSca3 = scratch1[3]; + + lookupTransmittance(sampleCosTheta, normalizedAltitude, scratch5); // tToSun + lookupMultiscattering(sampleCosTheta, normalizedAltitude, distToCenter, scratch1); // ms + + // stepTransmittance = exp(-extinction * dt) + var st0 = Math.exp(-scratch4[0] * dt); + var st1 = Math.exp(-scratch4[1] * dt); + var st2 = Math.exp(-scratch4[2] * dt); + var st3 = Math.exp(-scratch4[3] * dt); // S = SUN_SPECTRAL_IRRADIANCE * (molSca * (molPhase * tToSun + ms) + aerSca * (aerPhase * tToSun + ms)) - var s: Array = [0.0, 0.0, 0.0, 0.0]; - for (w in 0...4) { - var molTerm = coeffs.molSca[w] * (molPhase * tToSun[w] + ms[w]); - var aerTerm = coeffs.aerosolSca[w] * (aerPhase * tToSun[w] + ms[w]); - s[w] = SUN_SPECTRAL_IRRADIANCE[w] * (molTerm + aerTerm); - } + var s0 = SUN_SPECTRAL_IRRADIANCE[0] * (molSca0 * (molPhase * scratch5[0] + scratch1[0]) + aerSca0 * (aerPhase * scratch5[0] + scratch1[0])); + var s1 = SUN_SPECTRAL_IRRADIANCE[1] * (molSca1 * (molPhase * scratch5[1] + scratch1[1]) + aerSca1 * (aerPhase * scratch5[1] + scratch1[1])); + var s2 = SUN_SPECTRAL_IRRADIANCE[2] * (molSca2 * (molPhase * scratch5[2] + scratch1[2]) + aerSca2 * (aerPhase * scratch5[2] + scratch1[2])); + var s3 = SUN_SPECTRAL_IRRADIANCE[3] * (molSca3 * (molPhase * scratch5[3] + scratch1[3]) + aerSca3 * (aerPhase * scratch5[3] + scratch1[3])); - var stepTransmittance = exp4(scale4(extinction, -dt)); - var cutExt = max4(extinction, 1e-7); - var sInt: Array = []; - for (w in 0...4) { - sInt[w] = (s[w] - s[w] * stepTransmittance[w]) / cutExt[w]; - } - lInscattering = add4(lInscattering, mult4(transmittance, sInt)); - transmittance = mult4(transmittance, stepTransmittance); + // sInt[w] = (s[w] - s[w]*st[w]) / max(ext[w], 1e-7) + var ce0 = scratch4[0] < 1e-7 ? 1e-7 : scratch4[0]; + var ce1 = scratch4[1] < 1e-7 ? 1e-7 : scratch4[1]; + var ce2 = scratch4[2] < 1e-7 ? 1e-7 : scratch4[2]; + var ce3 = scratch4[3] < 1e-7 ? 1e-7 : scratch4[3]; + + var si0 = (s0 - s0 * st0) / ce0; + var si1 = (s1 - s1 * st1) / ce1; + var si2 = (s2 - s2 * st2) / ce2; + var si3 = (s3 - s3 * st3) / ce3; + + // lInscattering += transmittance * sInt + li0 += tr0 * si0; + li1 += tr1 * si1; + li2 += tr2 * si2; + li3 += tr3 * si3; + + // transmittance *= stepTransmittance + tr0 *= st0; tr1 *= st1; tr2 *= st2; tr3 *= st3; } - return lInscattering; - } + out[0] = li0; out[1] = li1; out[2] = li2; out[3] = li3; + } function computeEarthAngle(altitude: Float): Float { return Math.PI / 2.0 - Math.asin(EARTH_RADIUS / (EARTH_RADIUS + altitude / 1000.0)); } function precomputeSun(sunElevation: Float, angularDiameter: Float, altitude: Float): { bottom: Vec3, top: Vec3 } { - var halfAngular = angularDiameter / 2.0; + var halfAngular = Math.max(angularDiameter, 0.001) / 2.0; var solidAngle = 2.0 * Math.PI * (1.0 - Math.cos(halfAngular)); var normalizedAltitude = altitude / ATMOSPHERE_THICKNESS; function getSunXYZ(elevation: Float): Vec3 { var sunZenithCosAngle = Math.cos(Math.PI / 2.0 - elevation); - var tToSun = getTransmittance(sunZenithCosAngle, normalizedAltitude); - var spectrum: Array = []; - for (w in 0...4) { - spectrum[w] = SUN_SPECTRAL_IRRADIANCE[w] * tToSun[w] / solidAngle; - } - return spectralToXYZ(spectrum); + getTransmittance(sunZenithCosAngle, normalizedAltitude, scratch0); + return spectralToXYZScalars( + SUN_SPECTRAL_IRRADIANCE[0] * scratch0[0] / solidAngle, + SUN_SPECTRAL_IRRADIANCE[1] * scratch0[1] / solidAngle, + SUN_SPECTRAL_IRRADIANCE[2] * scratch0[2] / solidAngle, + SUN_SPECTRAL_IRRADIANCE[3] * scratch0[3] / solidAngle + ); } var bottom = getSunXYZ(sunElevation - halfAngular); @@ -609,8 +826,53 @@ class MultipleScatteringData { return { bottom: bottom, top: top }; } - public function compute(sunElevation: Float, sunRotation: Float, sunSize: Float, - sunIntensity: Float, altitude: Float, sunDisc: Bool, density: Vec3) { + // single row of the main LUT for incremental recompute + public function computeLUTRow(sdx: Float, sdy: Float, sdz: Float, roz: Float, + y: Int, imageData: haxe.io.Float32Array) { + + var uvY = (y + 0.5) / lutHeight; + var l = uvY * 2.0 - 1.0; + var elevation = (l < 0 ? -1.0 : 1.0) * l * l * Math.PI / 2.0; + var cosEl = Math.cos(elevation); + var sinEl = Math.sin(elevation); + + for (x in 0...lutWidth) { + var uvX = (x + 0.5) / lutWidth; + var azimuth = 2.0 * Math.PI * uvX; + + var dx = Math.sin(azimuth) * cosEl; + var dy = Math.cos(azimuth) * cosEl; + var dz = sinEl; + + var atmosDist = raySphereIntersectScalar(0, 0, roz, dx, dy, dz, ATMOSPHERE_RADIUS); + var groundDist = -1.0; + if (dz < 0.0) { + groundDist = raySphereIntersectScalar(0, 0, roz, dx, dy, dz, EARTH_RADIUS); + } + var tD = (groundDist < 0.0) ? atmosDist : groundDist; + + getInscatteringScalars(sdx, sdy, sdz, 0, 0, roz, dx, dy, dz, tD, scratch0); + var pixelIndex = (x + y * lutWidth) * 4; + spectralToXYZ(scratch0[0], scratch0[1], scratch0[2], scratch0[3], imageData, pixelIndex); + } + } + + // Sets sun disc data without creating a GPU texture. + public function setSunData(sunElevation: Float, sunSize: Float, altitude: Float, sunDisc: Bool) { + var altKm = Math.max(1.0, Math.min(99999.0, altitude)) / 1000.0; + earthIntersectionAngle = computeEarthAngle(altitude); + if (sunDisc) { + var sunData = precomputeSun(sunElevation, sunSize, altKm); + sunBottom = sunData.bottom; + sunTop = sunData.top; + } else { + sunBottom = new Vec3(); + sunTop = new Vec3(); + } + } + + public function computeData(sunElevation: Float, sunRotation: Float, sunSize: Float, + sunIntensity: Float, altitude: Float, sunDisc: Bool, density: Vec3): { data: haxe.io.Float32Array, sunBottom: Vec3, sunTop: Vec3, earthAngle: FastFloat } { airDensity = density.x; aerosolDensity = density.y; @@ -620,55 +882,57 @@ class MultipleScatteringData { var altKm = Math.max(1.0, Math.min(99999.0, altitude)) / 1000.0; var sunZenithCosAngle = Math.cos(Math.PI / 2.0 - sunElevation); - var sunDir = sunDirection(sunZenithCosAngle); - var rayOrigin = new Vec3(0.0, 0.0, EARTH_RADIUS + altKm); + var sdy = Math.sqrt(1.0 - sunZenithCosAngle * sunZenithCosAngle); + var sdz = sunZenithCosAngle; + var sdx = 0.0; + var roz = EARTH_RADIUS + altKm; - earthIntersectionAngle = computeEarthAngle(altitude); + var earthAngle = computeEarthAngle(altitude); var imageData = new haxe.io.Float32Array(lutWidth * lutHeight * 4); for (y in 0...lutHeight) { - var uvY = (y + 0.5) / lutHeight; - var l = uvY * 2.0 - 1.0; - var elevation = (l < 0 ? -1.0 : 1.0) * l * l * Math.PI / 2.0; - var cosEl = Math.cos(elevation); - var sinEl = Math.sin(elevation); - - for (x in 0...lutWidth) { - var uvX = (x + 0.5) / lutWidth; - var azimuth = 2.0 * Math.PI * uvX; - - // atan(dir.x, dir.y) convention - var rayDir = new Vec3( - Math.sin(azimuth) * cosEl, - Math.cos(azimuth) * cosEl, - sinEl - ).normalize(); - - var atmosDist = raySphereIntersection(rayOrigin, rayDir, ATMOSPHERE_RADIUS); - var groundDist = raySphereIntersection(rayOrigin, rayDir, EARTH_RADIUS); - var tD = (groundDist < 0.0) ? atmosDist : groundDist; - - var radiance = getInscattering(sunDir, rayOrigin, rayDir, tD); - var xyz = spectralToXYZ(radiance); - - var pixelIndex = (x + y * lutWidth) * 4; - imageData[pixelIndex + 0] = xyz.x; - imageData[pixelIndex + 1] = xyz.y; - imageData[pixelIndex + 2] = xyz.z; - imageData[pixelIndex + 3] = 1.0; - } + computeLUTRow(sdx, sdy, sdz, roz, y, imageData); } - lut = kha.Image.fromBytes(imageData.view.buffer, lutWidth, lutHeight, TextureFormat.RGBA128, Usage.StaticUsage); - + var bottom = new Vec3(); + var top = new Vec3(); if (sunDisc) { var sunData = precomputeSun(sunElevation, sunSize, altKm); - sunBottom = sunData.bottom; - sunTop = sunData.top; - } else { - sunBottom.set(0, 0, 0); - sunTop.set(0, 0, 0); + bottom = sunData.bottom; + top = sunData.top; } + + return { data: imageData, sunBottom: bottom, sunTop: top, earthAngle: earthAngle }; + } + + public function applyData(result: { data: haxe.io.Float32Array, sunBottom: Vec3, sunTop: Vec3, earthAngle: FastFloat }) { + lut = kha.Image.fromBytes(result.data.view.buffer, lutWidth, lutHeight, TextureFormat.RGBA128, Usage.StaticUsage); + sunBottom = result.sunBottom; + sunTop = result.sunTop; + earthIntersectionAngle = result.earthAngle; + } + + // LUT data for blending + public function applyLUTData(imageData: haxe.io.Float32Array) { + lut = kha.Image.fromBytes(imageData.view.buffer, lutWidth, lutHeight, TextureFormat.RGBA128, Usage.DynamicUsage); + } + + public function needsRecompute(density: Vec3, sunElevation: Float, altitude: Float, sunSize: Float, sunDisc: Bool): Bool { + return lut == null || needsTransmittance(density) || + sunElevation != lastElevation || altitude != lastAltitude || + sunSize != lastSunSize || sunDisc != lastSunDisc; + } + + public function compute(sunElevation: Float, sunRotation: Float, sunSize: Float, + sunIntensity: Float, altitude: Float, sunDisc: Bool, density: Vec3) { + + var result = computeData(sunElevation, sunRotation, sunSize, sunIntensity, altitude, sunDisc, density); + applyData(result); + lastLUTData = result.data; + lastElevation = sunElevation; + lastAltitude = altitude; + lastSunSize = sunSize; + lastSunDisc = sunDisc; } } diff --git a/leenkx/blender/lnx/make_world.py b/leenkx/blender/lnx/make_world.py index 6c9dec74..fd89eaab 100644 --- a/leenkx/blender/lnx/make_world.py +++ b/leenkx/blender/lnx/make_world.py @@ -239,9 +239,10 @@ def build_node_tree(world: bpy.types.World, frag: Shader, vert: Shader, con: Sha if '_EnvCol' in world.world_defs: frag.add_uniform('vec3 backgroundCol', link='_backgroundCol') frag.write('fragColor.rgb = backgroundCol * envmapStrength;') - - elif '_EnvTex' in world.world_defs and '_EnvLDR' in world.world_defs: - frag.write('fragColor.rgb = pow(fragColor.rgb, vec3(2.2));') + else: + if '_EnvTex' in world.world_defs and '_EnvLDR' in world.world_defs: + frag.write('fragColor.rgb = pow(fragColor.rgb, vec3(2.2));') + frag.write('fragColor.rgb *= envmapStrength;') if '_EnvClouds' in world.world_defs: frag.write('if (pos.z > 0.0) fragColor.rgb = mix(fragColor.rgb, traceClouds(fragColor.rgb, pos), clamp(pos.z * 5.0, 0, 1));')