Files
LNXSDK/leenkx/Sources/leenkx/renderpath/Sky.hx
2026-06-25 21:33:51 -07:00

939 lines
36 KiB
Haxe

package leenkx.renderpath;
import kha.FastFloat;
import kha.arrays.Float32Array;
import kha.graphics4.TextureFormat;
import kha.graphics4.Usage;
import iron.data.WorldData;
import iron.math.Vec3;
import leenkx.math.Helper;
/**
Utility class to control the sky models (single and multiple scattering).
**/
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.
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;
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.
public static function setDensity(world: WorldData, densityAir: FastFloat, densityDust: FastFloat, densityOzone: FastFloat) {
if (world == null) return;
if (world.raw.sky_density == null) world.raw.sky_density = new Float32Array(3);
var density = world.raw.sky_density;
density[0] = Helper.clamp(densityAir, 0, 10);
density[1] = Helper.clamp(densityDust, 0, 10);
density[2] = Helper.clamp(densityOzone, 0, 10);
recomputeSingleScatter(world);
}
public static function recomputeMultiScatter(world: WorldData) {
if (!registered) register();
if (world == null) return;
var raw = world.raw;
if (raw.sky_density == null) return;
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;
var sunIntensity = raw.sky_sun_intensity != null ? raw.sky_sun_intensity : 1.0;
var altitude = raw.sky_altitude != null ? raw.sky_altitude : 0.0;
var sunDisc = raw.sky_sun_disc != null ? raw.sky_sun_disc : 1;
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,
sunSize: FastFloat, sunIntensity: FastFloat, altitude: FastFloat, sunDisc: Bool,
densityAir: FastFloat, densityDust: FastFloat, densityOzone: FastFloat) {
if (world == null) return;
if (world.raw.sky_density == null) world.raw.sky_density = new Float32Array(3);
var density = world.raw.sky_density;
density[0] = Helper.clamp(densityAir, 0, 10);
density[1] = Helper.clamp(densityDust, 0, 10);
density[2] = Helper.clamp(densityOzone, 0, 10);
world.raw.sky_sun_elevation = sunElevation;
world.raw.sky_sun_rotation = sunRotation;
world.raw.sky_sun_size = sunSize;
world.raw.sky_sun_intensity = sunIntensity;
world.raw.sky_altitude = altitude;
world.raw.sky_sun_disc = sunDisc ? 1 : 0;
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;
}
}
}
}
}
/**
This class holds the precalculated result of the inner scattering integral
of the single scattering sky model. The outer integral is calculated in
[`leenkx/Shaders/std/sky.glsl`](https://github.com/leenkx3d/leenkx/blob/master/Shaders/std/sky.glsl).
@see `leenkx.renderpath.Sky`
**/
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).
**/
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;
/**
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 kilometers. **/
public static var radiusAtmo = 6420.0;
/**
Radius of the planet in kilometers. The default value is the earth radius as
defined in Cycles.
**/
public static var radiusPlanet = 6360.0;
/** Rayleigh scattering coefficient. **/
public static var rayleighCoeff = new Vec3(5.5e-6, 13.0e-6, 22.4e-6);
/** Rayleigh scattering scale parameter. **/
public static var rayleighScale = 8e3;
/** Mie scattering coefficient. **/
public static var mieCoeff = 2e-5;
/** Mie scattering scale parameter. **/
public static var mieScale = 1.2e3;
/** Ozone scattering coefficient. **/
// The ozone absorption coefficients are taken from Cycles code.
// Because Cycles calculates 21 wavelengths, we use the coefficients
// which are closest to the RGB wavelengths (645nm, 510nm, 440nm).
// Precalculating values by simulating Blender's spec_to_xyz() function
// 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. **/
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;
return -((height - 40000.0) / 15000.0);
}
/**
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
**/
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 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;
// Use quadratic height for better horizon precision
var height = (x / (lutHeightSteps - 1));
height *= height;
height *= radiusAtmo * 1000; // Denormalize height
var iz = height + radiusPlanetMeters;
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
var sy = Math.sin(sunTheta);
var sz = Math.cos(sunTheta);
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);
}
public function needsRecompute(density: Vec3): Bool {
return lut == null || lastDensity == null ||
lastDensity.x != density.x || lastDensity.y != density.y || lastDensity.z != density.z;
}
public function setLastDensity(density: Vec3) {
lastDensity = new Vec3(density.x, density.y, density.z);
}
public function computeLUT(density: Vec3) {
applyLUTData(computeLUTData(density));
lastDensity = new Vec3(density.x, density.y, density.z);
}
}
/**
This class precomputes the full sky radiance LUT
Unlike the single scattering model matching Blenders implementation.
**/
class MultipleScatteringData {
public var lut: kha.Image;
public var sunBottom: Vec3;
public var sunTop: Vec3;
public var earthIntersectionAngle: FastFloat;
public static var lutWidth = 128;
public static var lutHeight = 64;
static var transmittanceResX = 128;
public static var transmittanceResY = 32;
var transmittanceLUT: haxe.io.Float32Array;
var transmittanceValid = false;
var transmittanceDensity: Vec3 = null;
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
public static var EARTH_RADIUS = 6371.0;
static var ATMOSPHERE_THICKNESS = 100.0;
static var ATMOSPHERE_RADIUS = 6471.0; // EARTH_RADIUS + ATMOSPHERE_THICKNESS
static var PHASE_ISOTROPIC = 1.0 / (4.0 * Math.PI);
static var RAYLEIGH_PHASE_SCALE = (3.0 / 16.0) * (1.0 / Math.PI);
static var G = 0.8;
static var SQR_G = 0.64;
static var GROUND_ALBEDO: Array<Float> = [0.3, 0.3, 0.3, 0.3];
// Spectral data sampled at 630, 560, 490, 430 nm
static var SUN_SPECTRAL_IRRADIANCE: Array<Float> = [1.679, 1.828, 1.986, 1.307];
static var MOLECULAR_SCATTERING_COEFFICIENT_BASE: Array<Float> = [6.605e-3, 1.067e-2, 1.842e-2, 3.156e-2];
static var OZONE_ABSORPTION_CROSS_SECTION: Array<Float> = [3.472e-25, 3.914e-25, 1.349e-25, 11.03e-27];
static var OZONE_MEAN_DOBSON = 334.5;
static var AEROSOL_ABSORPTION_CROSS_SECTION: Array<Float> = [2.8722e-24, 4.6168e-24, 7.9706e-24, 1.3578e-23];
static var AEROSOL_SCATTERING_CROSS_SECTION: Array<Float> = [1.5908e-22, 1.7711e-22, 2.0942e-22, 2.4033e-22];
static var AEROSOL_BASE_DENSITY = 1.3681e20;
static var AEROSOL_BACKGROUND_DENSITY = 2e6;
static var AEROSOL_HEIGHT_SCALE = 0.73;
static var SPECTRAL_XYZ: Array<Array<Float>> = [
[53.386917738564668023, 22.981337506691024754, 0.0],
[43.904844466369358263, 71.347795700053393866, 0.102506867965741307],
[1.6137278251608962005, 18.422960591455485011, 31.742921188390805758],
[20.762668673810577145, 2.3614213523314368527, 110.48009643252140334]
];
var airDensity: Float = 1.0;
var aerosolDensity: Float = 1.0;
var ozoneDensity: Float = 1.0;
// avoid allocations in hot loops
static var scratch0: Array<Float> = [0, 0, 0, 0];
static var scratch1: Array<Float> = [0, 0, 0, 0];
static var scratch2: Array<Float> = [0, 0, 0, 0];
static var scratch3: Array<Float> = [0, 0, 0, 0];
static var scratch4: Array<Float> = [0, 0, 0, 0];
static var scratch5: Array<Float> = [0, 0, 0, 0];
public function new() {
sunBottom = new Vec3();
sunTop = new Vec3();
earthIntersectionAngle = 0.0;
transmittanceLUT = new haxe.io.Float32Array(transmittanceResY * transmittanceResX * 4);
}
static inline function add4(a: Array<Float>, b: Array<Float>, out: Array<Float>) {
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<Float>, s: Float, out: Array<Float>) {
out[0] = a[0] * s; out[1] = a[1] * s; out[2] = a[2] * s; out[3] = a[3] * s;
}
// Atmospheric functions
function getMolecularScatteringCoefficient(h: Float, out: Array<Float>) {
var f = Math.exp(-0.07771971 * Math.pow(h, 1.16364243));
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, out: Array<Float>) {
var logH = Math.log(Math.max(h, 1e-4));
var density = 3.78547397e20 * Math.exp(-(logH - 3.22261) * (logH - 3.22261) * 5.55555555 - logH);
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;
}
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) {
var localAerosolDensity = getAerosolDensity(h) * aerosolDensity;
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);
}
static inline function molecularPhaseFunction(cosTheta: Float): Float {
return RAYLEIGH_PHASE_SCALE * (1.0 + cosTheta * cosTheta);
}
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
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;
if (d >= b * b) return -b + Math.sqrt(d);
return -b - Math.sqrt(d);
}
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
);
}
// 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, out: Array<Float>) {
var sy = Math.sqrt(1.0 - cosTheta * cosTheta);
var sz = cosTheta;
var distToCenter = EARTH_RADIUS + (ATMOSPHERE_RADIUS - EARTH_RADIUS) * normalizedAltitude;
var tD = raySphereIntersectScalar(0, 0, distToCenter, 0, sy, sz, ATMOSPHERE_RADIUS);
var tStep = tD / transmittanceSteps;
// 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 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;
}
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) {
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];
}
}
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<Float>) {
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);
var fy = v * (transmittanceResY - 1);
var x1 = Std.int(fx);
var y1 = Std.int(fy);
var x2 = Std.int(Math.min(x1 + 1, transmittanceResX - 1));
var y2 = Std.int(Math.min(y1 + 1, transmittanceResY - 1));
var dxF = fx - x1;
var dyF = fy - y1;
var idx00 = (y1 * transmittanceResX + x1) * 4;
var idx01 = (y1 * transmittanceResX + x2) * 4;
var idx10 = (y2 * transmittanceResX + x1) * 4;
var idx11 = (y2 * transmittanceResX + x2) * 4;
// 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, out: Array<Float>) {
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 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]);
}
}
function lookupTransmittanceToSun(normalizedAltitude: Float, out: Array<Float>) {
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;
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]);
}
}
function lookupMultiscattering(cosTheta: Float, normalizedHeight: Float, d: Float, out: Array<Float>) {
var omega = 2.0 * Math.PI * (1.0 - Math.sqrt(Math.max(0.0, 1.0 - (EARTH_RADIUS / d) * (EARTH_RADIUS / d))));
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;
}
// 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<Float>) {
// cosTheta = -rayDir . sunDir
var cosTheta = -(dx * sdx + dy * sdy + dz * sdz);
var molPhase = molecularPhaseFunction(cosTheta);
var aerPhase = aerosolPhaseFunction(cosTheta);
var dt = tD / inScatteringSteps;
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;
// 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;
// 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;
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 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]));
// 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;
}
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 = 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);
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);
var top = getSunXYZ(sunElevation + halfAngular);
return { bottom: bottom, top: top };
}
// 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;
ozoneDensity = density.z;
precomputeTransmittanceLUT();
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 = EARTH_RADIUS + altKm;
var earthAngle = computeEarthAngle(altitude);
var imageData = new haxe.io.Float32Array(lutWidth * lutHeight * 4);
for (y in 0...lutHeight) {
computeLUTRow(sdx, sdy, sdz, roz, y, imageData);
}
var bottom = new Vec3();
var top = new Vec3();
if (sunDisc) {
var sunData = precomputeSun(sunElevation, sunSize, altKm);
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;
}
}