Procedural atmospheric scattering (#16314)

Implement procedural atmospheric scattering from [Sebastien Hillaire's
2020 paper](https://sebh.github.io/publications/egsr2020.pdf). This
approach should scale well even down to mobile hardware, and is
physically accurate.

## Co-author: @mate-h 

He helped massively with getting this over the finish line, ensuring
everything was physically correct, correcting several places where I had
misunderstood or misapplied the paper, and improving the performance in
several places as well. Thanks!

## Credits

@aevyrie: helped find numerous bugs and improve the example to best show
off this feature :)

Built off of @mtsr's original branch, which handled the transmittance
lut (arguably the most important part)

## Showcase: 


![sunset](https://github.com/user-attachments/assets/2eee1f38-f66d-4772-bb72-163e13c719d8)

![twilight](https://github.com/user-attachments/assets/f7d358b6-898d-4df7-becc-188cd753102d)


## For followup

- Integrate with pcwalton's volumetrics code
- refactor/reorganize for better integration with other effects
- have atmosphere transmittance affect directional lights
- add support for generating skybox/environment map

---------

Co-authored-by: Emerson Coskey <56370779+EmersonCoskey@users.noreply.github.com>
Co-authored-by: atlv <email@atlasdostal.com>
Co-authored-by: JMS55 <47158642+JMS55@users.noreply.github.com>
Co-authored-by: Emerson Coskey <coskey@emerlabs.net>
Co-authored-by: Máté Homolya <mate.homolya@gmail.com>

Cargo.toml

@@ -862,6 +862,17 @@ description = "A scene showcasing the atmospheric fog effect"
 category = "3D Rendering"
 wasm = true
 
+[[example]]
+name = "atmosphere"
+path = "examples/3d/atmosphere.rs"
+doc-scrape-examples = true
+
+[package.metadata.example.atmosphere]
+name = "Atmosphere"
+description = "A scene showcasing pbr atmospheric scattering"
+category = "3D Rendering"
+wasm = true
+
 [[example]]
 name = "fog"
 path = "examples/3d/fog.rs"

crates/bevy_pbr/src/atmosphere/aerial_view_lut.wgsl

@@ -0,0 +1,62 @@
+#import bevy_pbr::{
+    mesh_view_types::{Lights, DirectionalLight},
+    atmosphere::{
+        types::{Atmosphere, AtmosphereSettings},
+        bindings::{atmosphere, settings, view, lights, aerial_view_lut_out},
+        functions::{
+            sample_transmittance_lut, sample_atmosphere, rayleigh, henyey_greenstein,
+            sample_multiscattering_lut, AtmosphereSample, sample_local_inscattering,
+            get_local_r, get_local_up, view_radius, uv_to_ndc, max_atmosphere_distance,
+            uv_to_ray_direction
+        },
+    }
+}
+
+
+@group(0) @binding(13) var aerial_view_lut_out: texture_storage_3d<rgba16float, write>;
+
+@compute
+@workgroup_size(16, 16, 1)
+fn main(@builtin(global_invocation_id) idx: vec3<u32>) {
+    if any(idx.xy > settings.aerial_view_lut_size.xy) { return; }
+
+    let uv = (vec2<f32>(idx.xy) + 0.5) / vec2<f32>(settings.aerial_view_lut_size.xy);
+    let ray_dir = uv_to_ray_direction(uv);
+    let r = view_radius();
+    let mu = ray_dir.y;
+    let t_max = settings.aerial_view_lut_max_distance;
+
+    var prev_t = 0.0;
+    var total_inscattering = vec3(0.0);
+    var throughput = vec3(1.0);
+
+    for (var slice_i: u32 = 0; slice_i < settings.aerial_view_lut_size.z; slice_i++) {
+        for (var step_i: u32 = 0; step_i < settings.aerial_view_lut_samples; step_i++) {
+            let t_i = t_max * (f32(slice_i) + ((f32(step_i) + 0.5) / f32(settings.aerial_view_lut_samples))) / f32(settings.aerial_view_lut_size.z);
+            let dt = (t_i - prev_t);
+            prev_t = t_i;
+
+            let local_r = get_local_r(r, mu, t_i);
+            let local_up = get_local_up(r, t_i, ray_dir.xyz);
+
+            let local_atmosphere = sample_atmosphere(local_r);
+            let sample_optical_depth = local_atmosphere.extinction * dt;
+            let sample_transmittance = exp(-sample_optical_depth);
+
+            // evaluate one segment of the integral
+            var inscattering = sample_local_inscattering(local_atmosphere, ray_dir.xyz, local_r, local_up);
+
+            // Analytical integration of the single scattering term in the radiance transfer equation
+            let s_int = (inscattering - inscattering * sample_transmittance) / local_atmosphere.extinction;
+            total_inscattering += throughput * s_int;
+
+            throughput *= sample_transmittance;
+            if all(throughput < vec3(0.001)) {
+                break;
+            }
+        }
+        //We only have one channel to store transmittance, so we store the mean
+        let mean_transmittance = (throughput.r + throughput.g + throughput.b) / 3.0;
+        textureStore(aerial_view_lut_out, vec3(vec2<u32>(idx.xy), slice_i), vec4(total_inscattering, mean_transmittance));
+    }
+}

crates/bevy_pbr/src/atmosphere/bindings.wgsl

@@ -0,0 +1,22 @@
+#define_import_path bevy_pbr::atmosphere::bindings
+
+#import bevy_render::view::View;
+
+#import bevy_pbr::{
+    mesh_view_types::Lights,
+    atmosphere::types::{Atmosphere, AtmosphereSettings, AtmosphereTransforms}
+}
+
+@group(0) @binding(0) var<uniform> atmosphere: Atmosphere;
+@group(0) @binding(1) var<uniform> settings: AtmosphereSettings;
+@group(0) @binding(2) var<uniform> atmosphere_transforms: AtmosphereTransforms;
+@group(0) @binding(3) var<uniform> view: View;
+@group(0) @binding(4) var<uniform> lights: Lights;
+@group(0) @binding(5) var transmittance_lut: texture_2d<f32>;
+@group(0) @binding(6) var transmittance_lut_sampler: sampler;
+@group(0) @binding(7) var multiscattering_lut: texture_2d<f32>;
+@group(0) @binding(8) var multiscattering_lut_sampler: sampler;
+@group(0) @binding(9) var sky_view_lut: texture_2d<f32>;
+@group(0) @binding(10) var sky_view_lut_sampler: sampler;
+@group(0) @binding(11) var aerial_view_lut: texture_3d<f32>;
+@group(0) @binding(12) var aerial_view_lut_sampler: sampler;

crates/bevy_pbr/src/atmosphere/bruneton_functions.wgsl

@@ -0,0 +1,139 @@
+// Copyright (c) 2017 Eric Bruneton
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+// 1. Redistributions of source code must retain the above copyright
+//    notice, this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright
+//    notice, this list of conditions and the following disclaimer in the
+//    documentation and/or other materials provided with the distribution.
+// 3. Neither the name of the copyright holders nor the names of its
+//    contributors may be used to endorse or promote products derived from
+//    this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+// THE POSSIBILITY OF SUCH DAMAGE.
+//
+// Precomputed Atmospheric Scattering
+// Copyright (c) 2008 INRIA
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+// 1. Redistributions of source code must retain the above copyright
+//    notice, this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright
+//    notice, this list of conditions and the following disclaimer in the
+//    documentation and/or other materials provided with the distribution.
+// 3. Neither the name of the copyright holders nor the names of its
+//    contributors may be used to endorse or promote products derived from
+//    this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+// THE POSSIBILITY OF SUCH DAMAGE.
+
+#define_import_path bevy_pbr::atmosphere::bruneton_functions
+
+#import bevy_pbr::atmosphere::{
+    types::Atmosphere,
+    bindings::atmosphere,
+}
+
+// Mapping from view height (r) and zenith cos angle (mu) to UV coordinates in the transmittance LUT
+// Assuming r between ground and top atmosphere boundary, and mu= cos(zenith_angle)
+// Chosen to increase precision near the ground and to work around a discontinuity at the horizon
+// See Bruneton and Neyret 2008, "Precomputed Atmospheric Scattering" section 4
+fn transmittance_lut_r_mu_to_uv(r: f32, mu: f32) -> vec2<f32> {
+  // Distance along a horizontal ray from the ground to the top atmosphere boundary
+    let H = sqrt(atmosphere.top_radius * atmosphere.top_radius - atmosphere.bottom_radius * atmosphere.bottom_radius);
+
+  // Distance from a point at height r to the horizon
+  // ignore the case where r <= atmosphere.bottom_radius
+    let rho = sqrt(max(r * r - atmosphere.bottom_radius * atmosphere.bottom_radius, 0.0));
+
+  // Distance from a point at height r to the top atmosphere boundary at zenith angle mu
+    let d = distance_to_top_atmosphere_boundary(r, mu);
+
+  // Minimum and maximum distance to the top atmosphere boundary from a point at height r
+    let d_min = atmosphere.top_radius - r; // length of the ray straight up to the top atmosphere boundary
+    let d_max = rho + H; // length of the ray to the top atmosphere boundary and grazing the horizon
+
+    let u = (d - d_min) / (d_max - d_min);
+    let v = rho / H;
+    return vec2<f32>(u, v);
+}
+
+// Inverse of the mapping above, mapping from UV coordinates in the transmittance LUT to view height (r) and zenith cos angle (mu)
+fn transmittance_lut_uv_to_r_mu(uv: vec2<f32>) -> vec2<f32> {
+  // Distance to top atmosphere boundary for a horizontal ray at ground level
+    let H = sqrt(atmosphere.top_radius * atmosphere.top_radius - atmosphere.bottom_radius * atmosphere.bottom_radius);
+
+  // Distance to the horizon, from which we can compute r:
+    let rho = H * uv.y;
+    let r = sqrt(rho * rho + atmosphere.bottom_radius * atmosphere.bottom_radius);
+
+  // Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
+  // and maximum values over all mu- obtained for (r,1) and (r,mu_horizon) -
+  // from which we can recover mu:
+    let d_min = atmosphere.top_radius - r;
+    let d_max = rho + H;
+    let d = d_min + uv.x * (d_max - d_min);
+
+    var mu: f32;
+    if d == 0.0 {
+        mu = 1.0;
+    } else {
+        mu = (H * H - rho * rho - d * d) / (2.0 * r * d);
+    }
+
+    mu = clamp(mu, -1.0, 1.0);
+
+    return vec2<f32>(r, mu);
+}
+
+/// Simplified ray-sphere intersection
+/// where:
+/// Ray origin, o = [0,0,r] with r <= atmosphere.top_radius
+/// mu is the cosine of spherical coordinate theta (-1.0 <= mu <= 1.0)
+/// so ray direction in spherical coordinates is [1,acos(mu),0] which needs to be converted to cartesian
+/// Direction of ray, u = [0,sqrt(1-mu*mu),mu]
+/// Center of sphere, c = [0,0,0]
+/// Radius of sphere, r = atmosphere.top_radius
+/// This function solves the quadratic equation for line-sphere intersection simplified under these assumptions
+fn distance_to_top_atmosphere_boundary(r: f32, mu: f32) -> f32 {
+  // ignore the case where r > atmosphere.top_radius
+    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atmosphere.top_radius * atmosphere.top_radius, 0.0);
+    return max(-r * mu + sqrt(positive_discriminant), 0.0);
+}
+
+/// Simplified ray-sphere intersection
+/// as above for intersections with the ground
+fn distance_to_bottom_atmosphere_boundary(r: f32, mu: f32) -> f32 {
+    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atmosphere.bottom_radius * atmosphere.bottom_radius, 0.0);
+    return max(-r * mu - sqrt(positive_discriminant), 0.0);
+}
+
+fn ray_intersects_ground(r: f32, mu: f32) -> bool {
+    return mu < 0.0 && r * r * (mu * mu - 1.0) + atmosphere.bottom_radius * atmosphere.bottom_radius >= 0.0;
+}

crates/bevy_pbr/src/atmosphere/functions.wgsl

@@ -0,0 +1,334 @@
+#define_import_path bevy_pbr::atmosphere::functions
+
+#import bevy_render::maths::{PI, HALF_PI, PI_2, fast_acos, fast_atan2}
+
+#import bevy_pbr::atmosphere::{
+    types::Atmosphere,
+    bindings::{
+        atmosphere, settings, view, lights, transmittance_lut, transmittance_lut_sampler, 
+        multiscattering_lut, multiscattering_lut_sampler, sky_view_lut, sky_view_lut_sampler,
+        aerial_view_lut, aerial_view_lut_sampler, atmosphere_transforms
+    },
+    bruneton_functions::{
+        transmittance_lut_r_mu_to_uv, transmittance_lut_uv_to_r_mu, 
+        ray_intersects_ground, distance_to_top_atmosphere_boundary, 
+        distance_to_bottom_atmosphere_boundary
+    },
+}
+
+// NOTE FOR CONVENTIONS: 
+// r:
+//   radius, or distance from planet center 
+//
+// altitude:
+//   distance from planet **surface**
+//
+// mu:
+//   cosine of the zenith angle of a ray with
+//   respect to the planet normal
+//
+// atmosphere space:
+//   abbreviated as "as" (contrast with vs, cs, ws), this space is similar
+//   to view space, but with the camera positioned horizontally on the planet
+//   surface, so the horizon is a horizontal line centered vertically in the
+//   frame. This enables the non-linear latitude parametrization the paper uses 
+//   to concentrate detail near the horizon 
+
+
+// CONSTANTS
+
+const FRAC_PI: f32 = 0.3183098862; // 1 / π
+const FRAC_2_PI: f32 = 0.15915494309;
+const FRAC_3_16_PI: f32 = 0.0596831036594607509; // 3 / (16π)
+const FRAC_4_PI: f32 = 0.07957747154594767; // 1 / (4π)
+const ROOT_2: f32 = 1.41421356; // √2
+
+// LUT UV PARAMATERIZATIONS
+
+fn unit_to_sub_uvs(val: vec2<f32>, resolution: vec2<f32>) -> vec2<f32> {
+    return (val + 0.5f / resolution) * (resolution / (resolution + 1.0f));
+}
+
+fn sub_uvs_to_unit(val: vec2<f32>, resolution: vec2<f32>) -> vec2<f32> {
+    return (val - 0.5f / resolution) * (resolution / (resolution - 1.0f));
+}
+
+fn multiscattering_lut_r_mu_to_uv(r: f32, mu: f32) -> vec2<f32> {
+    let u = 0.5 + 0.5 * mu;
+    let v = saturate((r - atmosphere.bottom_radius) / (atmosphere.top_radius - atmosphere.bottom_radius)); //TODO
+    return unit_to_sub_uvs(vec2(u, v), vec2<f32>(settings.multiscattering_lut_size));
+}
+
+fn multiscattering_lut_uv_to_r_mu(uv: vec2<f32>) -> vec2<f32> {
+    let adj_uv = sub_uvs_to_unit(uv, vec2<f32>(settings.multiscattering_lut_size));
+    let r = mix(atmosphere.bottom_radius, atmosphere.top_radius, adj_uv.y);
+    let mu = adj_uv.x * 2 - 1;
+    return vec2(r, mu);
+}
+
+fn sky_view_lut_r_mu_azimuth_to_uv(r: f32, mu: f32, azimuth: f32) -> vec2<f32> {
+    let u = (azimuth * FRAC_2_PI) + 0.5;
+
+    let v_horizon = sqrt(r * r - atmosphere.bottom_radius * atmosphere.bottom_radius);
+    let cos_beta = v_horizon / r;
+    let beta = fast_acos(cos_beta);
+    let horizon_zenith = PI - beta;
+    let view_zenith = fast_acos(mu);
+
+    var v: f32;
+    if !ray_intersects_ground(r, mu) {
+        let coord = sqrt(1.0 - view_zenith / horizon_zenith);
+        v = (1.0 - coord) * 0.5;
+    } else {
+        let coord = (view_zenith - horizon_zenith) / beta;
+        v = sqrt(coord) * 0.5 + 0.5;
+    }
+
+    return unit_to_sub_uvs(vec2(u, v), vec2<f32>(settings.sky_view_lut_size));
+}
+
+fn sky_view_lut_uv_to_zenith_azimuth(r: f32, uv: vec2<f32>) -> vec2<f32> {
+    let adj_uv = sub_uvs_to_unit(uv, vec2<f32>(settings.sky_view_lut_size));
+    let azimuth = (adj_uv.x - 0.5) * PI_2;
+
+    let v_horizon = sqrt(r * r - atmosphere.bottom_radius * atmosphere.bottom_radius);
+    let cos_beta = v_horizon / r;
+    let beta = fast_acos(cos_beta);
+    let horizon_zenith = PI - beta;
+
+    var zenith: f32;
+    if adj_uv.y < 0.5 {
+        let coord = 1.0 - 2.0 * adj_uv.y;
+        zenith = horizon_zenith * (1.0 - coord * coord);
+    } else {
+        let coord = 2.0 * adj_uv.y - 1.0;
+        zenith = horizon_zenith + beta * coord * coord;
+    }
+
+    return vec2(zenith, azimuth);
+}
+
+// LUT SAMPLING
+
+fn sample_transmittance_lut(r: f32, mu: f32) -> vec3<f32> {
+    let uv = transmittance_lut_r_mu_to_uv(r, mu);
+    return textureSampleLevel(transmittance_lut, transmittance_lut_sampler, uv, 0.0).rgb;
+}
+
+fn sample_multiscattering_lut(r: f32, mu: f32) -> vec3<f32> {
+    let uv = multiscattering_lut_r_mu_to_uv(r, mu);
+    return textureSampleLevel(multiscattering_lut, multiscattering_lut_sampler, uv, 0.0).rgb;
+}
+
+fn sample_sky_view_lut(r: f32, ray_dir_as: vec3<f32>) -> vec3<f32> {
+    let mu = ray_dir_as.y;
+    let azimuth = fast_atan2(ray_dir_as.x, -ray_dir_as.z);
+    let uv = sky_view_lut_r_mu_azimuth_to_uv(r, mu, azimuth);
+    return textureSampleLevel(sky_view_lut, sky_view_lut_sampler, uv, 0.0).rgb;
+}
+
+//RGB channels: total inscattered light along the camera ray to the current sample.
+//A channel: average transmittance across all wavelengths to the current sample.
+fn sample_aerial_view_lut(uv: vec2<f32>, depth: f32) -> vec4<f32> {
+    let view_pos = view.view_from_clip * vec4(uv_to_ndc(uv), depth, 1.0);
+    let dist = length(view_pos.xyz / view_pos.w) * settings.scene_units_to_m;
+    let uvw = vec3(uv, dist / settings.aerial_view_lut_max_distance);
+    return textureSampleLevel(aerial_view_lut, aerial_view_lut_sampler, uvw, 0.0);
+}
+
+// PHASE FUNCTIONS
+
+// -(L . V) == (L . -V). -V here is our ray direction, which points away from the view 
+// instead of towards it (which would be the *view direction*, V)
+
+// evaluates the rayleigh phase function, which describes the likelihood
+// of a rayleigh scattering event scattering light from the light direction towards the view
+fn rayleigh(neg_LdotV: f32) -> f32 {
+    return FRAC_3_16_PI * (1 + (neg_LdotV * neg_LdotV));
+}
+
+// evaluates the henyey-greenstein phase function, which describes the likelihood
+// of a mie scattering event scattering light from the light direction towards the view
+fn henyey_greenstein(neg_LdotV: f32) -> f32 {
+    let g = atmosphere.mie_asymmetry;
+    let denom = 1.0 + g * g - 2.0 * g * neg_LdotV;
+    return FRAC_4_PI * (1.0 - g * g) / (denom * sqrt(denom));
+}
+
+// ATMOSPHERE SAMPLING
+
+struct AtmosphereSample {
+    /// units: m^-1
+    rayleigh_scattering: vec3<f32>,
+
+    /// units: m^-1
+    mie_scattering: f32,
+
+    /// the sum of scattering and absorption. Since the phase function doesn't
+    /// matter for this, we combine rayleigh and mie extinction to a single 
+    //  value.
+    //
+    /// units: m^-1
+    extinction: vec3<f32>
+}
+
+/// Samples atmosphere optical densities at a given radius
+fn sample_atmosphere(r: f32) -> AtmosphereSample {
+    let altitude = clamp(r, atmosphere.bottom_radius, atmosphere.top_radius) - atmosphere.bottom_radius;
+
+    // atmosphere values at altitude
+    let mie_density = exp(-atmosphere.mie_density_exp_scale * altitude);
+    let rayleigh_density = exp(-atmosphere.rayleigh_density_exp_scale * altitude);
+    var ozone_density: f32 = max(0.0, 1.0 - (abs(altitude - atmosphere.ozone_layer_altitude) / (atmosphere.ozone_layer_width * 0.5)));
+
+    let mie_scattering = mie_density * atmosphere.mie_scattering;
+    let mie_absorption = mie_density * atmosphere.mie_absorption;
+    let mie_extinction = mie_scattering + mie_absorption;
+
+    let rayleigh_scattering = rayleigh_density * atmosphere.rayleigh_scattering;
+    // no rayleigh absorption
+    // rayleigh extinction is the sum of scattering and absorption
+
+    // ozone doesn't contribute to scattering
+    let ozone_absorption = ozone_density * atmosphere.ozone_absorption;
+
+    var sample: AtmosphereSample;
+    sample.rayleigh_scattering = rayleigh_scattering;
+    sample.mie_scattering = mie_scattering;
+    sample.extinction = rayleigh_scattering + mie_extinction + ozone_absorption;
+
+    return sample;
+}
+
+/// evaluates L_scat, equation 3 in the paper, which gives the total single-order scattering towards the view at a single point
+fn sample_local_inscattering(local_atmosphere: AtmosphereSample, ray_dir: vec3<f32>, local_r: f32, local_up: vec3<f32>) -> vec3<f32> {
+    var inscattering = vec3(0.0);
+    for (var light_i: u32 = 0u; light_i < lights.n_directional_lights; light_i++) {
+        let light = &lights.directional_lights[light_i];
+
+        let mu_light = dot((*light).direction_to_light, local_up);
+
+        // -(L . V) == (L . -V). -V here is our ray direction, which points away from the view
+        // instead of towards it (as is the convention for V)
+        let neg_LdotV = dot((*light).direction_to_light, ray_dir);
+
+        // Phase functions give the proportion of light
+        // scattered towards the camera for each scattering type
+        let rayleigh_phase = rayleigh(neg_LdotV);
+        let mie_phase = henyey_greenstein(neg_LdotV);
+        let scattering_coeff = local_atmosphere.rayleigh_scattering * rayleigh_phase + local_atmosphere.mie_scattering * mie_phase;
+
+        let transmittance_to_light = sample_transmittance_lut(local_r, mu_light);
+        let shadow_factor = transmittance_to_light * f32(!ray_intersects_ground(local_r, mu_light));
+
+        // Transmittance from scattering event to light source
+        let scattering_factor = shadow_factor * scattering_coeff;
+
+        // Additive factor from the multiscattering LUT
+        let psi_ms = sample_multiscattering_lut(local_r, mu_light);
+        let multiscattering_factor = psi_ms * (local_atmosphere.rayleigh_scattering + local_atmosphere.mie_scattering);
+
+        inscattering += (*light).color.rgb * (scattering_factor + multiscattering_factor);
+    }
+    return inscattering * view.exposure;
+}
+
+const SUN_ANGULAR_SIZE: f32 = 0.0174533; // angular diameter of sun in radians
+
+fn sample_sun_illuminance(ray_dir_ws: vec3<f32>, transmittance: vec3<f32>) -> vec3<f32> {
+    var sun_illuminance = vec3(0.0);
+    for (var light_i: u32 = 0u; light_i < lights.n_directional_lights; light_i++) {
+        let light = &lights.directional_lights[light_i];
+        let neg_LdotV = dot((*light).direction_to_light, ray_dir_ws);
+        let angle_to_sun = fast_acos(neg_LdotV);
+        let pixel_size = fwidth(angle_to_sun);
+        let factor = smoothstep(0.0, -pixel_size * ROOT_2, angle_to_sun - SUN_ANGULAR_SIZE * 0.5);
+        let sun_solid_angle = (SUN_ANGULAR_SIZE * SUN_ANGULAR_SIZE) * 4.0 * FRAC_PI;
+        sun_illuminance += ((*light).color.rgb / sun_solid_angle) * factor * ray_dir_ws.y;
+    }
+    return sun_illuminance * transmittance * view.exposure;
+}
+
+// TRANSFORM UTILITIES
+
+fn max_atmosphere_distance(r: f32, mu: f32) -> f32 {
+    let t_top = distance_to_top_atmosphere_boundary(r, mu);
+    let t_bottom = distance_to_bottom_atmosphere_boundary(r, mu);
+    let hits = ray_intersects_ground(r, mu);
+    return mix(t_top, t_bottom, f32(hits));
+}
+
+/// Assuming y=0 is the planet ground, returns the view radius in meters
+fn view_radius() -> f32 {
+    return view.world_position.y * settings.scene_units_to_m + atmosphere.bottom_radius;
+}
+
+// We assume the `up` vector at the view position is the y axis, since the world is locally flat/level.
+// t = distance along view ray in atmosphere space
+// NOTE: this means that if your world is actually spherical, this will be wrong.
+fn get_local_up(r: f32, t: f32, ray_dir: vec3<f32>) -> vec3<f32> {
+    return normalize(vec3(0.0, r, 0.0) + t * ray_dir);
+}
+
+// Given a ray starting at radius r, with mu = cos(zenith angle),
+// and a t = distance along the ray, gives the new radius at point t
+fn get_local_r(r: f32, mu: f32, t: f32) -> f32 {
+    return sqrt(t * t + 2.0 * r * mu * t + r * r);
+}
+
+// Convert uv [0.0 .. 1.0] coordinate to ndc space xy [-1.0 .. 1.0]
+fn uv_to_ndc(uv: vec2<f32>) -> vec2<f32> {
+    return uv * vec2(2.0, -2.0) + vec2(-1.0, 1.0);
+}
+
+/// Convert ndc space xy coordinate [-1.0 .. 1.0] to uv [0.0 .. 1.0]
+fn ndc_to_uv(ndc: vec2<f32>) -> vec2<f32> {
+    return ndc * vec2(0.5, -0.5) + vec2(0.5);
+}
+
+/// Converts a direction in world space to atmosphere space
+fn direction_world_to_atmosphere(dir_ws: vec3<f32>) -> vec3<f32> {
+    let dir_as = atmosphere_transforms.atmosphere_from_world * vec4(dir_ws, 0.0);
+    return dir_as.xyz;
+}
+
+/// Converts a direction in atmosphere space to world space
+fn direction_atmosphere_to_world(dir_as: vec3<f32>) -> vec3<f32> {
+    let dir_ws = atmosphere_transforms.world_from_atmosphere * vec4(dir_as, 0.0);
+    return dir_ws.xyz;
+}
+
+// Modified from skybox.wgsl. For this pass we don't need to apply a separate sky transform or consider camera viewport.
+// w component is the cosine of the view direction with the view forward vector, to correct step distance at the edges of the viewport
+fn uv_to_ray_direction(uv: vec2<f32>) -> vec4<f32> {
+    // Using world positions of the fragment and camera to calculate a ray direction
+    // breaks down at large translations. This code only needs to know the ray direction.
+    // The ray direction is along the direction from the camera to the fragment position.
+    // In view space, the camera is at the origin, so the view space ray direction is
+    // along the direction of the fragment position - (0,0,0) which is just the
+    // fragment position.
+    // Use the position on the near clipping plane to avoid -inf world position
+    // because the far plane of an infinite reverse projection is at infinity.
+    let view_position_homogeneous = view.view_from_clip * vec4(
+        uv_to_ndc(uv),
+        1.0,
+        1.0,
+    );
+
+    let view_ray_direction = view_position_homogeneous.xyz / view_position_homogeneous.w;
+    // Transforming the view space ray direction by the inverse view matrix, transforms the
+    // direction to world space. Note that the w element is set to 0.0, as this is a
+    // vector direction, not a position, That causes the matrix multiplication to ignore
+    // the translations from the view matrix.
+    let ray_direction = (view.world_from_view * vec4(view_ray_direction, 0.0)).xyz;
+
+    return vec4(normalize(ray_direction), -view_ray_direction.z);
+}
+
+fn zenith_azimuth_to_ray_dir(zenith: f32, azimuth: f32) -> vec3<f32> {
+    let sin_zenith = sin(zenith);
+    let mu = cos(zenith);
+    let sin_azimuth = sin(azimuth);
+    let cos_azimuth = cos(azimuth);
+    return vec3(sin_azimuth * sin_zenith, mu, -cos_azimuth * sin_zenith);
+}

crates/bevy_pbr/src/atmosphere/mod.rs

@@ -0,0 +1,469 @@
+//! Procedural Atmospheric Scattering.
+//!
+//! This plugin implements [Hillaire's 2020 paper](https://sebh.github.io/publications/egsr2020.pdf)
+//! on real-time atmospheric scattering. While it *will* work simply as a
+//! procedural skybox, it also does much more. It supports dynamic time-of-
+//! -day, multiple directional lights, and since it's applied as a post-processing
+//! effect *on top* of the existing skybox, a starry skybox would automatically
+//! show based on the time of day. Scattering in front of terrain (similar
+//! to distance fog, but more complex) is handled as well, and takes into
+//! account the directional light color and direction.
+//!
+//! Adding the [`Atmosphere`] component to a 3d camera will enable the effect,
+//! which by default is set to look similar to Earth's atmosphere. See the
+//! documentation on the component itself for information regarding its fields.
+//!
+//! Performance-wise, the effect should be fairly cheap since the LUTs (Look
+//! Up Tables) that encode most of the data are small, and take advantage of the
+//! fact that the atmosphere is symmetric. Performance is also proportional to
+//! the number of directional lights in the scene. In order to tune
+//! performance more finely, the [`AtmosphereSettings`] camera component
+//! manages the size of each LUT and the sample count for each ray.
+//!
+//! Given how similar it is to [`crate::volumetric_fog`], it might be expected
+//! that these two modules would work together well. However for now using both
+//! at once is untested, and might not be physically accurate. These may be
+//! integrated into a single module in the future.
+//!
+//! [Shadertoy]: https://www.shadertoy.com/view/slSXRW
+//!
+//! [Unreal Engine Implementation]: https://github.com/sebh/UnrealEngineSkyAtmosphere
+
+mod node;
+pub mod resources;
+
+use bevy_app::{App, Plugin};
+use bevy_asset::load_internal_asset;
+use bevy_core_pipeline::core_3d::graph::Node3d;
+use bevy_ecs::{
+    component::{require, Component},
+    query::{Changed, QueryItem, With},
+    schedule::IntoSystemConfigs,
+    system::{lifetimeless::Read, Query},
+};
+use bevy_math::{UVec2, UVec3, Vec3};
+use bevy_reflect::Reflect;
+use bevy_render::{
+    extract_component::UniformComponentPlugin,
+    render_resource::{DownlevelFlags, ShaderType, SpecializedRenderPipelines},
+};
+use bevy_render::{
+    extract_component::{ExtractComponent, ExtractComponentPlugin},
+    render_graph::{RenderGraphApp, ViewNodeRunner},
+    render_resource::{Shader, TextureFormat, TextureUsages},
+    renderer::RenderAdapter,
+    Render, RenderApp, RenderSet,
+};
+
+use bevy_core_pipeline::core_3d::{graph::Core3d, Camera3d};
+use resources::{
+    prepare_atmosphere_transforms, queue_render_sky_pipelines, AtmosphereTransforms,
+    RenderSkyBindGroupLayouts,
+};
+use tracing::warn;
+
+use self::{
+    node::{AtmosphereLutsNode, AtmosphereNode, RenderSkyNode},
+    resources::{
+        prepare_atmosphere_bind_groups, prepare_atmosphere_textures, AtmosphereBindGroupLayouts,
+        AtmosphereLutPipelines, AtmosphereSamplers,
+    },
+};
+
+mod shaders {
+    use bevy_asset::Handle;
+    use bevy_render::render_resource::Shader;
+
+    pub const TYPES: Handle<Shader> = Handle::weak_from_u128(0xB4CA686B10FA592B508580CCC2F9558C);
+    pub const FUNCTIONS: Handle<Shader> =
+        Handle::weak_from_u128(0xD5524FD88BDC153FBF256B7F2C21906F);
+    pub const BRUNETON_FUNCTIONS: Handle<Shader> =
+        Handle::weak_from_u128(0x7E896F48B707555DD11985F9C1594459);
+    pub const BINDINGS: Handle<Shader> = Handle::weak_from_u128(0x140EFD89B5D4C8490AB895010DFC42FE);
+
+    pub const TRANSMITTANCE_LUT: Handle<Shader> =
+        Handle::weak_from_u128(0xEECBDEDFEED7F4EAFBD401BFAA5E0EFB);
+    pub const MULTISCATTERING_LUT: Handle<Shader> =
+        Handle::weak_from_u128(0x65915B32C44B6287C0CCE1E70AF2936A);
+    pub const SKY_VIEW_LUT: Handle<Shader> =
+        Handle::weak_from_u128(0x54136D7E6FFCD45BE38399A4E5ED7186);
+    pub const AERIAL_VIEW_LUT: Handle<Shader> =
+        Handle::weak_from_u128(0x6FDEC284AD356B78C3A4D8ED4CBA0BC5);
+    pub const RENDER_SKY: Handle<Shader> =
+        Handle::weak_from_u128(0x1951EB87C8A6129F0B541B1E4B3D4962);
+}
+
+#[doc(hidden)]
+pub struct AtmospherePlugin;
+
+impl Plugin for AtmospherePlugin {
+    fn build(&self, app: &mut App) {
+        load_internal_asset!(app, shaders::TYPES, "types.wgsl", Shader::from_wgsl);
+        load_internal_asset!(app, shaders::FUNCTIONS, "functions.wgsl", Shader::from_wgsl);
+        load_internal_asset!(
+            app,
+            shaders::BRUNETON_FUNCTIONS,
+            "bruneton_functions.wgsl",
+            Shader::from_wgsl
+        );
+
+        load_internal_asset!(app, shaders::BINDINGS, "bindings.wgsl", Shader::from_wgsl);
+
+        load_internal_asset!(
+            app,
+            shaders::TRANSMITTANCE_LUT,
+            "transmittance_lut.wgsl",
+            Shader::from_wgsl
+        );
+
+        load_internal_asset!(
+            app,
+            shaders::MULTISCATTERING_LUT,
+            "multiscattering_lut.wgsl",
+            Shader::from_wgsl
+        );
+
+        load_internal_asset!(
+            app,
+            shaders::SKY_VIEW_LUT,
+            "sky_view_lut.wgsl",
+            Shader::from_wgsl
+        );
+
+        load_internal_asset!(
+            app,
+            shaders::AERIAL_VIEW_LUT,
+            "aerial_view_lut.wgsl",
+            Shader::from_wgsl
+        );
+
+        load_internal_asset!(
+            app,
+            shaders::RENDER_SKY,
+            "render_sky.wgsl",
+            Shader::from_wgsl
+        );
+
+        app.register_type::<Atmosphere>()
+            .register_type::<AtmosphereSettings>()
+            .add_plugins((
+                ExtractComponentPlugin::<Atmosphere>::default(),
+                ExtractComponentPlugin::<AtmosphereSettings>::default(),
+                UniformComponentPlugin::<Atmosphere>::default(),
+                UniformComponentPlugin::<AtmosphereSettings>::default(),
+            ));
+    }
+
+    fn finish(&self, app: &mut App) {
+        let Some(render_app) = app.get_sub_app_mut(RenderApp) else {
+            return;
+        };
+
+        let render_adapter = render_app.world().resource::<RenderAdapter>();
+
+        if !render_adapter
+            .get_downlevel_capabilities()
+            .flags
+            .contains(DownlevelFlags::COMPUTE_SHADERS)
+        {
+            warn!("AtmospherePlugin not loaded. GPU lacks support for compute shaders.");
+            return;
+        }
+
+        if !render_adapter
+            .get_texture_format_features(TextureFormat::Rgba16Float)
+            .allowed_usages
+            .contains(TextureUsages::STORAGE_BINDING)
+        {
+            warn!("AtmospherePlugin not loaded. GPU lacks support: TextureFormat::Rgba16Float does not support TextureUsages::STORAGE_BINDING.");
+            return;
+        }
+
+        render_app
+            .init_resource::<AtmosphereBindGroupLayouts>()
+            .init_resource::<RenderSkyBindGroupLayouts>()
+            .init_resource::<AtmosphereSamplers>()
+            .init_resource::<AtmosphereLutPipelines>()
+            .init_resource::<AtmosphereTransforms>()
+            .init_resource::<SpecializedRenderPipelines<RenderSkyBindGroupLayouts>>()
+            .add_systems(
+                Render,
+                (
+                    configure_camera_depth_usages.in_set(RenderSet::ManageViews),
+                    queue_render_sky_pipelines.in_set(RenderSet::Queue),
+                    prepare_atmosphere_textures.in_set(RenderSet::PrepareResources),
+                    prepare_atmosphere_transforms.in_set(RenderSet::PrepareResources),
+                    prepare_atmosphere_bind_groups.in_set(RenderSet::PrepareBindGroups),
+                ),
+            )
+            .add_render_graph_node::<ViewNodeRunner<AtmosphereLutsNode>>(
+                Core3d,
+                AtmosphereNode::RenderLuts,
+            )
+            .add_render_graph_edges(
+                Core3d,
+                (
+                    // END_PRE_PASSES -> RENDER_LUTS -> MAIN_PASS
+                    Node3d::EndPrepasses,
+                    AtmosphereNode::RenderLuts,
+                    Node3d::StartMainPass,
+                ),
+            )
+            .add_render_graph_node::<ViewNodeRunner<RenderSkyNode>>(
+                Core3d,
+                AtmosphereNode::RenderSky,
+            )
+            .add_render_graph_edges(
+                Core3d,
+                (
+                    Node3d::MainOpaquePass,
+                    AtmosphereNode::RenderSky,
+                    Node3d::MainTransparentPass,
+                ),
+            );
+    }
+}
+
+/// This component describes the atmosphere of a planet, and when added to a camera
+/// will enable atmospheric scattering for that camera. This is only compatible with
+/// HDR cameras.
+///
+/// Most atmospheric particles scatter and absorb light in two main ways:
+///
+/// Rayleigh scattering occurs among very small particles, like individual gas
+/// molecules. It's wavelength dependent, and causes colors to separate out as
+/// light travels through the atmosphere. These particles *don't* absorb light.
+///
+/// Mie scattering occurs among slightly larger particles, like dust and sea spray.
+/// These particles *do* absorb light, but Mie scattering and absorption is
+/// *wavelength independent*.
+///
+/// Ozone acts differently from the other two, and is special-cased because
+/// it's very important to the look of Earth's atmosphere. It's wavelength
+/// dependent, but only *absorbs* light. Also, while the density of particles
+/// participating in Rayleigh and Mie scattering falls off roughly exponentially
+/// from the planet's surface, ozone only exists in a band centered at a fairly
+/// high altitude.
+#[derive(Clone, Component, Reflect, ShaderType)]
+#[require(AtmosphereSettings)]
+pub struct Atmosphere {
+    /// Radius of the planet
+    ///
+    /// units: m
+    pub bottom_radius: f32,
+
+    /// Radius at which we consider the atmosphere to 'end' for our
+    /// calculations (from center of planet)
+    ///
+    /// units: m
+    pub top_radius: f32,
+
+    /// An approximation of the average albedo (or color, roughly) of the
+    /// planet's surface. This is used when calculating multiscattering.
+    ///
+    /// units: N/A
+    pub ground_albedo: Vec3,
+
+    /// The rate of falloff of rayleigh particulate with respect to altitude:
+    /// optical density = exp(-rayleigh_density_exp_scale * altitude in meters).
+    ///
+    /// THIS VALUE MUST BE POSITIVE
+    ///
+    /// units: N/A
+    pub rayleigh_density_exp_scale: f32,
+
+    /// The scattering optical density of rayleigh particulate, or how
+    /// much light it scatters per meter
+    ///
+    /// units: m^-1
+    pub rayleigh_scattering: Vec3,
+
+    /// The rate of falloff of mie particulate with respect to altitude:
+    /// optical density = exp(-mie_density_exp_scale * altitude in meters)
+    ///
+    /// THIS VALUE MUST BE POSITIVE
+    ///
+    /// units: N/A
+    pub mie_density_exp_scale: f32,
+
+    /// The scattering optical density of mie particulate, or how much light
+    /// it scatters per meter.
+    ///
+    /// units: m^-1
+    pub mie_scattering: f32,
+
+    /// The absorbing optical density of mie particulate, or how much light
+    /// it absorbs per meter.
+    ///
+    /// units: m^-1
+    pub mie_absorption: f32,
+
+    /// The "asymmetry" of mie scattering, or how much light tends to scatter
+    /// forwards, rather than backwards or to the side.
+    ///
+    /// domain: (-1, 1)
+    /// units: N/A
+    pub mie_asymmetry: f32, //the "asymmetry" value of the phase function, unitless. Domain: (-1, 1)
+
+    /// The altitude at which the ozone layer is centered.
+    ///
+    /// units: m
+    pub ozone_layer_altitude: f32,
+
+    /// The width of the ozone layer
+    ///
+    /// units: m
+    pub ozone_layer_width: f32,
+
+    /// The optical density of ozone, or how much of each wavelength of
+    /// light it absorbs per meter.
+    ///
+    /// units: m^-1
+    pub ozone_absorption: Vec3,
+}
+
+impl Atmosphere {
+    pub const EARTH: Atmosphere = Atmosphere {
+        bottom_radius: 6_360_000.0,
+        top_radius: 6_460_000.0,
+        ground_albedo: Vec3::splat(0.3),
+        rayleigh_density_exp_scale: 1.0 / 8_000.0,
+        rayleigh_scattering: Vec3::new(5.802e-6, 13.558e-6, 33.100e-6),
+        mie_density_exp_scale: 1.0 / 1_200.0,
+        mie_scattering: 3.996e-6,
+        mie_absorption: 0.444e-6,
+        mie_asymmetry: 0.8,
+        ozone_layer_altitude: 25_000.0,
+        ozone_layer_width: 30_000.0,
+        ozone_absorption: Vec3::new(0.650e-6, 1.881e-6, 0.085e-6),
+    };
+
+    pub fn with_density_multiplier(mut self, mult: f32) -> Self {
+        self.rayleigh_scattering *= mult;
+        self.mie_scattering *= mult;
+        self.mie_absorption *= mult;
+        self.ozone_absorption *= mult;
+        self
+    }
+}
+
+impl Default for Atmosphere {
+    fn default() -> Self {
+        Self::EARTH
+    }
+}
+
+impl ExtractComponent for Atmosphere {
+    type QueryData = Read<Atmosphere>;
+
+    type QueryFilter = With<Camera3d>;
+
+    type Out = Atmosphere;
+
+    fn extract_component(item: QueryItem<'_, Self::QueryData>) -> Option<Self::Out> {
+        Some(item.clone())
+    }
+}
+
+/// This component controls the resolution of the atmosphere LUTs, and
+/// how many samples are used when computing them.
+///
+/// The transmittance LUT stores the transmittance from a point in the
+/// atmosphere to the outer edge of the atmosphere in any direction,
+/// parametrized by the point's radius and the cosine of the zenith angle
+/// of the ray.
+///
+/// The multiscattering LUT stores the factor representing luminance scattered
+/// towards the camera with scattering order >2, parametrized by the point's radius
+/// and the cosine of the zenith angle of the sun.
+///
+/// The sky-view lut is essentially the actual skybox, storing the light scattered
+/// towards the camera in every direction with a cubemap.
+///
+/// The aerial-view lut is a 3d LUT fit to the view frustum, which stores the luminance
+/// scattered towards the camera at each point (RGB channels), alongside the average
+/// transmittance to that point (A channel).
+#[derive(Clone, Component, Reflect, ShaderType)]
+pub struct AtmosphereSettings {
+    /// The size of the transmittance LUT
+    pub transmittance_lut_size: UVec2,
+
+    /// The size of the multiscattering LUT
+    pub multiscattering_lut_size: UVec2,
+
+    /// The size of the sky-view LUT.
+    pub sky_view_lut_size: UVec2,
+
+    /// The size of the aerial-view LUT.
+    pub aerial_view_lut_size: UVec3,
+
+    /// The number of points to sample along each ray when
+    /// computing the transmittance LUT
+    pub transmittance_lut_samples: u32,
+
+    /// The number of rays to sample when computing each
+    /// pixel of the multiscattering LUT
+    pub multiscattering_lut_dirs: u32,
+
+    /// The number of points to sample when integrating along each
+    /// multiscattering ray
+    pub multiscattering_lut_samples: u32,
+
+    /// The number of points to sample along each ray when
+    /// computing the sky-view LUT.
+    pub sky_view_lut_samples: u32,
+
+    /// The number of points to sample for each slice along the z-axis
+    /// of the aerial-view LUT.
+    pub aerial_view_lut_samples: u32,
+
+    /// The maximum distance from the camera to evaluate the
+    /// aerial view LUT. The slices along the z-axis of the
+    /// texture will be distributed linearly from the camera
+    /// to this value.
+    ///
+    /// units: m
+    pub aerial_view_lut_max_distance: f32,
+
+    /// A conversion factor between scene units and meters, used to
+    /// ensure correctness at different length scales.
+    pub scene_units_to_m: f32,
+}
+
+impl Default for AtmosphereSettings {
+    fn default() -> Self {
+        Self {
+            transmittance_lut_size: UVec2::new(256, 128),
+            transmittance_lut_samples: 40,
+            multiscattering_lut_size: UVec2::new(32, 32),
+            multiscattering_lut_dirs: 64,
+            multiscattering_lut_samples: 20,
+            sky_view_lut_size: UVec2::new(400, 200),
+            sky_view_lut_samples: 16,
+            aerial_view_lut_size: UVec3::new(32, 32, 32),
+            aerial_view_lut_samples: 10,
+            aerial_view_lut_max_distance: 3.2e4,
+            scene_units_to_m: 1.0,
+        }
+    }
+}
+
+impl ExtractComponent for AtmosphereSettings {
+    type QueryData = Read<AtmosphereSettings>;
+
+    type QueryFilter = (With<Camera3d>, With<Atmosphere>);
+
+    type Out = AtmosphereSettings;
+
+    fn extract_component(item: QueryItem<'_, Self::QueryData>) -> Option<Self::Out> {
+        Some(item.clone())
+    }
+}
+
+fn configure_camera_depth_usages(
+    mut cameras: Query<&mut Camera3d, (Changed<Camera3d>, With<Atmosphere>)>,
+) {
+    for mut camera in &mut cameras {
+        camera.depth_texture_usages.0 |= TextureUsages::TEXTURE_BINDING.bits();
+    }
+}

crates/bevy_pbr/src/atmosphere/multiscattering_lut.wgsl

@@ -0,0 +1,139 @@
+#import bevy_pbr::{
+    mesh_view_types::{Lights, DirectionalLight},
+    atmosphere::{
+        types::{Atmosphere, AtmosphereSettings},
+        bindings::{atmosphere, settings},
+        functions::{
+            multiscattering_lut_uv_to_r_mu, sample_transmittance_lut,
+            get_local_r, get_local_up, sample_atmosphere, FRAC_4_PI,
+            max_atmosphere_distance, rayleigh, henyey_greenstein,
+            zenith_azimuth_to_ray_dir,
+        },
+        bruneton_functions::{
+            distance_to_top_atmosphere_boundary, distance_to_bottom_atmosphere_boundary, ray_intersects_ground
+        }
+    }
+}
+
+#import bevy_render::maths::{PI,PI_2}
+
+const PHI_2: vec2<f32> = vec2(1.3247179572447460259609088, 1.7548776662466927600495087);
+
+@group(0) @binding(13) var multiscattering_lut_out: texture_storage_2d<rgba16float, write>;
+
+fn s2_sequence(n: u32) -> vec2<f32> {
+    return fract(0.5 + f32(n) * PHI_2);
+}
+
+// Lambert equal-area projection. 
+fn uv_to_sphere(uv: vec2<f32>) -> vec3<f32> {
+    let phi = PI_2 * uv.y;
+    let sin_lambda = 2 * uv.x - 1;
+    let cos_lambda = sqrt(1 - sin_lambda * sin_lambda);
+
+    return vec3(cos_lambda * cos(phi), cos_lambda * sin(phi), sin_lambda);
+}
+
+// Shared memory arrays for workgroup communication
+var<workgroup> multi_scat_shared_mem: array<vec3<f32>, 64>;
+var<workgroup> l_shared_mem: array<vec3<f32>, 64>;
+
+@compute 
+@workgroup_size(1, 1, 64)
+fn main(@builtin(global_invocation_id) global_id: vec3<u32>) {
+    var uv = (vec2<f32>(global_id.xy) + 0.5) / vec2<f32>(settings.multiscattering_lut_size);
+
+    let r_mu = multiscattering_lut_uv_to_r_mu(uv);
+    let light_dir = normalize(vec3(0.0, r_mu.y, -1.0));
+
+    let ray_dir = uv_to_sphere(s2_sequence(global_id.z));
+    let ms_sample = sample_multiscattering_dir(r_mu.x, ray_dir, light_dir);
+    
+    // Calculate the contribution for this sample
+    let sphere_solid_angle = 4.0 * PI;
+    let sample_weight = sphere_solid_angle / 64.0;
+    multi_scat_shared_mem[global_id.z] = ms_sample.f_ms * sample_weight;
+    l_shared_mem[global_id.z] = ms_sample.l_2 * sample_weight;
+
+    workgroupBarrier();
+
+    // Parallel reduction bitshift to the right to divide by 2 each step
+    for (var step = 32u; step > 0u; step >>= 1u) {
+        if global_id.z < step {
+            multi_scat_shared_mem[global_id.z] += multi_scat_shared_mem[global_id.z + step];
+            l_shared_mem[global_id.z] += l_shared_mem[global_id.z + step];
+        }
+        workgroupBarrier();
+    }
+
+    if global_id.z > 0u {
+        return;
+    }
+
+    // Apply isotropic phase function
+    let f_ms = multi_scat_shared_mem[0] * FRAC_4_PI;
+    let l_2 = l_shared_mem[0] * FRAC_4_PI;
+    
+    // Equation 10 from the paper: Geometric series for infinite scattering
+    let psi_ms = l_2 / (1.0 - f_ms);
+    textureStore(multiscattering_lut_out, global_id.xy, vec4<f32>(psi_ms, 1.0));
+}
+
+struct MultiscatteringSample {
+    l_2: vec3<f32>,
+    f_ms: vec3<f32>,
+};
+
+fn sample_multiscattering_dir(r: f32, ray_dir: vec3<f32>, light_dir: vec3<f32>) -> MultiscatteringSample {
+    // get the cosine of the zenith angle of the view direction with respect to the light direction
+    let mu_view = ray_dir.y;
+    let t_max = max_atmosphere_distance(r, mu_view);
+
+    let dt = t_max / f32(settings.multiscattering_lut_samples);
+    var optical_depth = vec3<f32>(0.0);
+
+    var l_2 = vec3(0.0);
+    var f_ms = vec3(0.0);
+    var throughput = vec3(1.0);
+    for (var i: u32 = 0u; i < settings.multiscattering_lut_samples; i++) {
+        let t_i = dt * (f32(i) + 0.5);
+        let local_r = get_local_r(r, mu_view, t_i);
+        let local_up = get_local_up(r, t_i, ray_dir);
+
+        let local_atmosphere = sample_atmosphere(local_r);
+        let sample_optical_depth = local_atmosphere.extinction * dt;
+        let sample_transmittance = exp(-sample_optical_depth);
+        optical_depth += sample_optical_depth;
+
+        let mu_light = dot(light_dir, local_up);
+        let scattering_no_phase = local_atmosphere.rayleigh_scattering + local_atmosphere.mie_scattering;
+
+        let ms = scattering_no_phase;
+        let ms_int = (ms - ms * sample_transmittance) / local_atmosphere.extinction;
+        f_ms += throughput * ms_int;
+
+        let transmittance_to_light = sample_transmittance_lut(local_r, mu_light);
+        let shadow_factor = transmittance_to_light * f32(!ray_intersects_ground(local_r, mu_light));
+
+        let s = scattering_no_phase * shadow_factor * FRAC_4_PI;
+        let s_int = (s - s * sample_transmittance) / local_atmosphere.extinction;
+        l_2 += throughput * s_int;
+
+        throughput *= sample_transmittance;
+        if all(throughput < vec3(0.001)) {
+            break;
+        }
+    }
+
+    //include reflected luminance from planet ground 
+    if ray_intersects_ground(r, mu_view) {
+        let transmittance_to_ground = exp(-optical_depth);
+        let local_up = get_local_up(r, t_max, ray_dir);
+        let mu_light = dot(light_dir, local_up);
+        let transmittance_to_light = sample_transmittance_lut(0.0, mu_light);
+        let ground_luminance = transmittance_to_light * transmittance_to_ground * max(mu_light, 0.0) * atmosphere.ground_albedo;
+        l_2 += ground_luminance;
+    }
+
+    return MultiscatteringSample(l_2, f_ms);
+}

crates/bevy_pbr/src/atmosphere/node.rs

@@ -0,0 +1,221 @@
+use bevy_ecs::{query::QueryItem, system::lifetimeless::Read, world::World};
+use bevy_math::{UVec2, Vec3Swizzles};
+use bevy_render::{
+    extract_component::DynamicUniformIndex,
+    render_graph::{NodeRunError, RenderGraphContext, RenderLabel, ViewNode},
+    render_resource::{ComputePass, ComputePassDescriptor, PipelineCache, RenderPassDescriptor},
+    renderer::RenderContext,
+    view::{ViewTarget, ViewUniformOffset},
+};
+
+use crate::ViewLightsUniformOffset;
+
+use super::{
+    resources::{
+        AtmosphereBindGroups, AtmosphereLutPipelines, AtmosphereTransformsOffset,
+        RenderSkyPipelineId,
+    },
+    Atmosphere, AtmosphereSettings,
+};
+
+#[derive(PartialEq, Eq, Debug, Copy, Clone, Hash, RenderLabel)]
+pub enum AtmosphereNode {
+    RenderLuts,
+    RenderSky,
+}
+
+#[derive(Default)]
+pub(super) struct AtmosphereLutsNode {}
+
+impl ViewNode for AtmosphereLutsNode {
+    type ViewQuery = (
+        Read<AtmosphereSettings>,
+        Read<AtmosphereBindGroups>,
+        Read<DynamicUniformIndex<Atmosphere>>,
+        Read<DynamicUniformIndex<AtmosphereSettings>>,
+        Read<AtmosphereTransformsOffset>,
+        Read<ViewUniformOffset>,
+        Read<ViewLightsUniformOffset>,
+    );
+
+    fn run(
+        &self,
+        _graph: &mut RenderGraphContext,
+        render_context: &mut RenderContext,
+        (
+            settings,
+            bind_groups,
+            atmosphere_uniforms_offset,
+            settings_uniforms_offset,
+            atmosphere_transforms_offset,
+            view_uniforms_offset,
+            lights_uniforms_offset,
+        ): QueryItem<Self::ViewQuery>,
+        world: &World,
+    ) -> Result<(), NodeRunError> {
+        let pipelines = world.resource::<AtmosphereLutPipelines>();
+        let pipeline_cache = world.resource::<PipelineCache>();
+        let (
+            Some(transmittance_lut_pipeline),
+            Some(multiscattering_lut_pipeline),
+            Some(sky_view_lut_pipeline),
+            Some(aerial_view_lut_pipeline),
+        ) = (
+            pipeline_cache.get_compute_pipeline(pipelines.transmittance_lut),
+            pipeline_cache.get_compute_pipeline(pipelines.multiscattering_lut),
+            pipeline_cache.get_compute_pipeline(pipelines.sky_view_lut),
+            pipeline_cache.get_compute_pipeline(pipelines.aerial_view_lut),
+        )
+        else {
+            return Ok(());
+        };
+
+        let command_encoder = render_context.command_encoder();
+
+        let mut luts_pass = command_encoder.begin_compute_pass(&ComputePassDescriptor {
+            label: Some("atmosphere_luts_pass"),
+            timestamp_writes: None,
+        });
+
+        fn dispatch_2d(compute_pass: &mut ComputePass, size: UVec2) {
+            const WORKGROUP_SIZE: u32 = 16;
+            let workgroups_x = size.x.div_ceil(WORKGROUP_SIZE);
+            let workgroups_y = size.y.div_ceil(WORKGROUP_SIZE);
+            compute_pass.dispatch_workgroups(workgroups_x, workgroups_y, 1);
+        }
+
+        // Transmittance LUT
+
+        luts_pass.set_pipeline(transmittance_lut_pipeline);
+        luts_pass.set_bind_group(
+            0,
+            &bind_groups.transmittance_lut,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+            ],
+        );
+
+        dispatch_2d(&mut luts_pass, settings.transmittance_lut_size);
+
+        // Multiscattering LUT
+
+        luts_pass.set_pipeline(multiscattering_lut_pipeline);
+        luts_pass.set_bind_group(
+            0,
+            &bind_groups.multiscattering_lut,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+            ],
+        );
+
+        luts_pass.dispatch_workgroups(
+            settings.multiscattering_lut_size.x,
+            settings.multiscattering_lut_size.y,
+            1,
+        );
+
+        // Sky View LUT
+
+        luts_pass.set_pipeline(sky_view_lut_pipeline);
+        luts_pass.set_bind_group(
+            0,
+            &bind_groups.sky_view_lut,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+                atmosphere_transforms_offset.index(),
+                view_uniforms_offset.offset,
+                lights_uniforms_offset.offset,
+            ],
+        );
+
+        dispatch_2d(&mut luts_pass, settings.sky_view_lut_size);
+
+        // Aerial View LUT
+
+        luts_pass.set_pipeline(aerial_view_lut_pipeline);
+        luts_pass.set_bind_group(
+            0,
+            &bind_groups.aerial_view_lut,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+                view_uniforms_offset.offset,
+                lights_uniforms_offset.offset,
+            ],
+        );
+
+        dispatch_2d(&mut luts_pass, settings.aerial_view_lut_size.xy());
+
+        Ok(())
+    }
+}
+
+#[derive(Default)]
+pub(super) struct RenderSkyNode;
+
+impl ViewNode for RenderSkyNode {
+    type ViewQuery = (
+        Read<AtmosphereBindGroups>,
+        Read<ViewTarget>,
+        Read<DynamicUniformIndex<Atmosphere>>,
+        Read<DynamicUniformIndex<AtmosphereSettings>>,
+        Read<AtmosphereTransformsOffset>,
+        Read<ViewUniformOffset>,
+        Read<ViewLightsUniformOffset>,
+        Read<RenderSkyPipelineId>,
+    );
+
+    fn run<'w>(
+        &self,
+        _graph: &mut RenderGraphContext,
+        render_context: &mut RenderContext<'w>,
+        (
+            atmosphere_bind_groups,
+            view_target,
+            atmosphere_uniforms_offset,
+            settings_uniforms_offset,
+            atmosphere_transforms_offset,
+            view_uniforms_offset,
+            lights_uniforms_offset,
+            render_sky_pipeline_id,
+        ): QueryItem<'w, Self::ViewQuery>,
+        world: &'w World,
+    ) -> Result<(), NodeRunError> {
+        let pipeline_cache = world.resource::<PipelineCache>();
+        let Some(render_sky_pipeline) =
+            pipeline_cache.get_render_pipeline(render_sky_pipeline_id.0)
+        else {
+            return Ok(());
+        }; //TODO: warning
+
+        let mut render_sky_pass =
+            render_context
+                .command_encoder()
+                .begin_render_pass(&RenderPassDescriptor {
+                    label: Some("render_sky_pass"),
+                    color_attachments: &[Some(view_target.get_color_attachment())],
+                    depth_stencil_attachment: None,
+                    timestamp_writes: None,
+                    occlusion_query_set: None,
+                });
+
+        render_sky_pass.set_pipeline(render_sky_pipeline);
+        render_sky_pass.set_bind_group(
+            0,
+            &atmosphere_bind_groups.render_sky,
+            &[
+                atmosphere_uniforms_offset.index(),
+                settings_uniforms_offset.index(),
+                atmosphere_transforms_offset.index(),
+                view_uniforms_offset.offset,
+                lights_uniforms_offset.offset,
+            ],
+        );
+        render_sky_pass.draw(0..3, 0..1);
+
+        Ok(())
+    }
+}

crates/bevy_pbr/src/atmosphere/render_sky.wgsl

@@ -0,0 +1,39 @@
+#import bevy_pbr::atmosphere::{
+    types::{Atmosphere, AtmosphereSettings},
+    bindings::{atmosphere, view, atmosphere_transforms},
+    functions::{
+        sample_transmittance_lut, sample_sky_view_lut, 
+        direction_world_to_atmosphere, uv_to_ray_direction,
+        uv_to_ndc, sample_aerial_view_lut, view_radius,
+        sample_sun_illuminance,
+    },
+};
+#import bevy_render::view::View;
+
+#import bevy_core_pipeline::fullscreen_vertex_shader::FullscreenVertexOutput
+
+#ifdef MULTISAMPLED
+@group(0) @binding(13) var depth_texture: texture_depth_multisampled_2d;
+#else
+@group(0) @binding(13) var depth_texture: texture_depth_2d;
+#endif
+
+@fragment
+fn main(in: FullscreenVertexOutput) -> @location(0) vec4<f32> {
+    let depth = textureLoad(depth_texture, vec2<i32>(in.position.xy), 0);
+    if depth == 0.0 {
+        let ray_dir_ws = uv_to_ray_direction(in.uv);
+        let ray_dir_as = direction_world_to_atmosphere(ray_dir_ws.xyz);
+
+        let r = view_radius();
+        let mu = ray_dir_ws.y;
+
+        let transmittance = sample_transmittance_lut(r, mu);
+        let inscattering = sample_sky_view_lut(r, ray_dir_as);
+
+        let sun_illuminance = sample_sun_illuminance(ray_dir_ws.xyz, transmittance);
+        return vec4(inscattering + sun_illuminance, (transmittance.r + transmittance.g + transmittance.b) / 3.0);
+    } else {
+        return sample_aerial_view_lut(in.uv, depth);
+    }
+}

crates/bevy_pbr/src/atmosphere/resources.rs

@@ -0,0 +1,718 @@
+use bevy_core_pipeline::{
+    core_3d::Camera3d, fullscreen_vertex_shader::fullscreen_shader_vertex_state,
+};
+use bevy_ecs::{
+    component::Component,
+    entity::Entity,
+    query::With,
+    resource::Resource,
+    system::{Commands, Query, Res, ResMut},
+    world::{FromWorld, World},
+};
+use bevy_math::{Mat4, Vec3};
+use bevy_render::{
+    camera::Camera,
+    extract_component::ComponentUniforms,
+    render_resource::{binding_types::*, *},
+    renderer::{RenderDevice, RenderQueue},
+    texture::{CachedTexture, TextureCache},
+    view::{ExtractedView, Msaa, ViewDepthTexture, ViewUniform, ViewUniforms},
+};
+
+use crate::{GpuLights, LightMeta};
+
+use super::{shaders, Atmosphere, AtmosphereSettings};
+
+#[derive(Resource)]
+pub(crate) struct AtmosphereBindGroupLayouts {
+    pub transmittance_lut: BindGroupLayout,
+    pub multiscattering_lut: BindGroupLayout,
+    pub sky_view_lut: BindGroupLayout,
+    pub aerial_view_lut: BindGroupLayout,
+}
+
+#[derive(Resource)]
+pub(crate) struct RenderSkyBindGroupLayouts {
+    pub render_sky: BindGroupLayout,
+    pub render_sky_msaa: BindGroupLayout,
+}
+
+impl FromWorld for AtmosphereBindGroupLayouts {
+    fn from_world(world: &mut World) -> Self {
+        let render_device = world.resource::<RenderDevice>();
+        let transmittance_lut = render_device.create_bind_group_layout(
+            "transmittance_lut_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::COMPUTE,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<AtmosphereSettings>(true)),
+                    (
+                        // transmittance lut storage texture
+                        13,
+                        texture_storage_2d(
+                            TextureFormat::Rgba16Float,
+                            StorageTextureAccess::WriteOnly,
+                        ),
+                    ),
+                ),
+            ),
+        );
+
+        let multiscattering_lut = render_device.create_bind_group_layout(
+            "multiscattering_lut_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::COMPUTE,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<AtmosphereSettings>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (
+                        //multiscattering lut storage texture
+                        13,
+                        texture_storage_2d(
+                            TextureFormat::Rgba16Float,
+                            StorageTextureAccess::WriteOnly,
+                        ),
+                    ),
+                ),
+            ),
+        );
+
+        let sky_view_lut = render_device.create_bind_group_layout(
+            "sky_view_lut_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::COMPUTE,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<AtmosphereSettings>(true)),
+                    (2, uniform_buffer::<AtmosphereTransform>(true)),
+                    (3, uniform_buffer::<ViewUniform>(true)),
+                    (4, uniform_buffer::<GpuLights>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (7, texture_2d(TextureSampleType::Float { filterable: true })), //multiscattering lut and sampler
+                    (8, sampler(SamplerBindingType::Filtering)),
+                    (
+                        13,
+                        texture_storage_2d(
+                            TextureFormat::Rgba16Float,
+                            StorageTextureAccess::WriteOnly,
+                        ),
+                    ),
+                ),
+            ),
+        );
+
+        let aerial_view_lut = render_device.create_bind_group_layout(
+            "aerial_view_lut_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::COMPUTE,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<AtmosphereSettings>(true)),
+                    (3, uniform_buffer::<ViewUniform>(true)),
+                    (4, uniform_buffer::<GpuLights>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (7, texture_2d(TextureSampleType::Float { filterable: true })), //multiscattering lut and sampler
+                    (8, sampler(SamplerBindingType::Filtering)),
+                    (
+                        //Aerial view lut storage texture
+                        13,
+                        texture_storage_3d(
+                            TextureFormat::Rgba16Float,
+                            StorageTextureAccess::WriteOnly,
+                        ),
+                    ),
+                ),
+            ),
+        );
+
+        Self {
+            transmittance_lut,
+            multiscattering_lut,
+            sky_view_lut,
+            aerial_view_lut,
+        }
+    }
+}
+
+impl FromWorld for RenderSkyBindGroupLayouts {
+    fn from_world(world: &mut World) -> Self {
+        let render_device = world.resource::<RenderDevice>();
+        let render_sky = render_device.create_bind_group_layout(
+            "render_sky_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::FRAGMENT,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<AtmosphereSettings>(true)),
+                    (2, uniform_buffer::<AtmosphereTransform>(true)),
+                    (3, uniform_buffer::<ViewUniform>(true)),
+                    (4, uniform_buffer::<GpuLights>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (9, texture_2d(TextureSampleType::Float { filterable: true })), //sky view lut and sampler
+                    (10, sampler(SamplerBindingType::Filtering)),
+                    (
+                        // aerial view lut and sampler
+                        11,
+                        texture_3d(TextureSampleType::Float { filterable: true }),
+                    ),
+                    (12, sampler(SamplerBindingType::Filtering)),
+                    (
+                        //view depth texture
+                        13,
+                        texture_2d(TextureSampleType::Depth),
+                    ),
+                ),
+            ),
+        );
+
+        let render_sky_msaa = render_device.create_bind_group_layout(
+            "render_sky_msaa_bind_group_layout",
+            &BindGroupLayoutEntries::with_indices(
+                ShaderStages::FRAGMENT,
+                (
+                    (0, uniform_buffer::<Atmosphere>(true)),
+                    (1, uniform_buffer::<AtmosphereSettings>(true)),
+                    (2, uniform_buffer::<AtmosphereTransform>(true)),
+                    (3, uniform_buffer::<ViewUniform>(true)),
+                    (4, uniform_buffer::<GpuLights>(true)),
+                    (5, texture_2d(TextureSampleType::Float { filterable: true })), //transmittance lut and sampler
+                    (6, sampler(SamplerBindingType::Filtering)),
+                    (9, texture_2d(TextureSampleType::Float { filterable: true })), //sky view lut and sampler
+                    (10, sampler(SamplerBindingType::Filtering)),
+                    (
+                        // aerial view lut and sampler
+                        11,
+                        texture_3d(TextureSampleType::Float { filterable: true }),
+                    ),
+                    (12, sampler(SamplerBindingType::Filtering)),
+                    (
+                        //view depth texture
+                        13,
+                        texture_2d_multisampled(TextureSampleType::Depth),
+                    ),
+                ),
+            ),
+        );
+
+        Self {
+            render_sky,
+            render_sky_msaa,
+        }
+    }
+}
+
+#[derive(Resource)]
+pub struct AtmosphereSamplers {
+    pub transmittance_lut: Sampler,
+    pub multiscattering_lut: Sampler,
+    pub sky_view_lut: Sampler,
+    pub aerial_view_lut: Sampler,
+}
+
+impl FromWorld for AtmosphereSamplers {
+    fn from_world(world: &mut World) -> Self {
+        let render_device = world.resource::<RenderDevice>();
+
+        let base_sampler = SamplerDescriptor {
+            mag_filter: FilterMode::Linear,
+            min_filter: FilterMode::Linear,
+            mipmap_filter: FilterMode::Nearest,
+            ..Default::default()
+        };
+
+        let transmittance_lut = render_device.create_sampler(&SamplerDescriptor {
+            label: Some("transmittance_lut_sampler"),
+            ..base_sampler
+        });
+
+        let multiscattering_lut = render_device.create_sampler(&SamplerDescriptor {
+            label: Some("multiscattering_lut_sampler"),
+            ..base_sampler
+        });
+
+        let sky_view_lut = render_device.create_sampler(&SamplerDescriptor {
+            label: Some("sky_view_lut_sampler"),
+            address_mode_u: AddressMode::Repeat,
+            ..base_sampler
+        });
+
+        let aerial_view_lut = render_device.create_sampler(&SamplerDescriptor {
+            label: Some("aerial_view_lut_sampler"),
+            ..base_sampler
+        });
+
+        Self {
+            transmittance_lut,
+            multiscattering_lut,
+            sky_view_lut,
+            aerial_view_lut,
+        }
+    }
+}
+
+#[derive(Resource)]
+pub(crate) struct AtmosphereLutPipelines {
+    pub transmittance_lut: CachedComputePipelineId,
+    pub multiscattering_lut: CachedComputePipelineId,
+    pub sky_view_lut: CachedComputePipelineId,
+    pub aerial_view_lut: CachedComputePipelineId,
+}
+
+impl FromWorld for AtmosphereLutPipelines {
+    fn from_world(world: &mut World) -> Self {
+        let pipeline_cache = world.resource::<PipelineCache>();
+        let layouts = world.resource::<AtmosphereBindGroupLayouts>();
+
+        let transmittance_lut = pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+            label: Some("transmittance_lut_pipeline".into()),
+            layout: vec![layouts.transmittance_lut.clone()],
+            push_constant_ranges: vec![],
+            shader: shaders::TRANSMITTANCE_LUT,
+            shader_defs: vec![],
+            entry_point: "main".into(),
+            zero_initialize_workgroup_memory: false,
+        });
+
+        let multiscattering_lut =
+            pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+                label: Some("multi_scattering_lut_pipeline".into()),
+                layout: vec![layouts.multiscattering_lut.clone()],
+                push_constant_ranges: vec![],
+                shader: shaders::MULTISCATTERING_LUT,
+                shader_defs: vec![],
+                entry_point: "main".into(),
+                zero_initialize_workgroup_memory: false,
+            });
+
+        let sky_view_lut = pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+            label: Some("sky_view_lut_pipeline".into()),
+            layout: vec![layouts.sky_view_lut.clone()],
+            push_constant_ranges: vec![],
+            shader: shaders::SKY_VIEW_LUT,
+            shader_defs: vec![],
+            entry_point: "main".into(),
+            zero_initialize_workgroup_memory: false,
+        });
+
+        let aerial_view_lut = pipeline_cache.queue_compute_pipeline(ComputePipelineDescriptor {
+            label: Some("aerial_view_lut_pipeline".into()),
+            layout: vec![layouts.aerial_view_lut.clone()],
+            push_constant_ranges: vec![],
+            shader: shaders::AERIAL_VIEW_LUT,
+            shader_defs: vec![],
+            entry_point: "main".into(),
+            zero_initialize_workgroup_memory: false,
+        });
+
+        Self {
+            transmittance_lut,
+            multiscattering_lut,
+            sky_view_lut,
+            aerial_view_lut,
+        }
+    }
+}
+
+#[derive(Component)]
+pub(crate) struct RenderSkyPipelineId(pub CachedRenderPipelineId);
+
+#[derive(Copy, Clone, Hash, PartialEq, Eq)]
+pub(crate) struct RenderSkyPipelineKey {
+    pub msaa_samples: u32,
+    pub hdr: bool,
+}
+
+impl SpecializedRenderPipeline for RenderSkyBindGroupLayouts {
+    type Key = RenderSkyPipelineKey;
+
+    fn specialize(&self, key: Self::Key) -> RenderPipelineDescriptor {
+        let mut shader_defs = Vec::new();
+
+        if key.msaa_samples > 1 {
+            shader_defs.push("MULTISAMPLED".into());
+        }
+        if key.hdr {
+            shader_defs.push("TONEMAP_IN_SHADER".into());
+        }
+
+        RenderPipelineDescriptor {
+            label: Some(format!("render_sky_pipeline_{}", key.msaa_samples).into()),
+            layout: vec![if key.msaa_samples == 1 {
+                self.render_sky.clone()
+            } else {
+                self.render_sky_msaa.clone()
+            }],
+            push_constant_ranges: vec![],
+            vertex: fullscreen_shader_vertex_state(),
+            primitive: PrimitiveState::default(),
+            depth_stencil: None,
+            multisample: MultisampleState {
+                count: key.msaa_samples,
+                mask: !0,
+                alpha_to_coverage_enabled: false,
+            },
+            zero_initialize_workgroup_memory: false,
+            fragment: Some(FragmentState {
+                shader: shaders::RENDER_SKY.clone(),
+                shader_defs,
+                entry_point: "main".into(),
+                targets: vec![Some(ColorTargetState {
+                    format: TextureFormat::Rgba16Float,
+                    blend: Some(BlendState {
+                        color: BlendComponent {
+                            src_factor: BlendFactor::One,
+                            dst_factor: BlendFactor::SrcAlpha,
+                            operation: BlendOperation::Add,
+                        },
+                        alpha: BlendComponent {
+                            src_factor: BlendFactor::Zero,
+                            dst_factor: BlendFactor::One,
+                            operation: BlendOperation::Add,
+                        },
+                    }),
+                    write_mask: ColorWrites::ALL,
+                })],
+            }),
+        }
+    }
+}
+
+pub(super) fn queue_render_sky_pipelines(
+    views: Query<(Entity, &Camera, &Msaa), With<Atmosphere>>,
+    pipeline_cache: Res<PipelineCache>,
+    layouts: Res<RenderSkyBindGroupLayouts>,
+    mut specializer: ResMut<SpecializedRenderPipelines<RenderSkyBindGroupLayouts>>,
+    mut commands: Commands,
+) {
+    for (entity, camera, msaa) in &views {
+        let id = specializer.specialize(
+            &pipeline_cache,
+            &layouts,
+            RenderSkyPipelineKey {
+                msaa_samples: msaa.samples(),
+                hdr: camera.hdr,
+            },
+        );
+        commands.entity(entity).insert(RenderSkyPipelineId(id));
+    }
+}
+
+#[derive(Component)]
+pub struct AtmosphereTextures {
+    pub transmittance_lut: CachedTexture,
+    pub multiscattering_lut: CachedTexture,
+    pub sky_view_lut: CachedTexture,
+    pub aerial_view_lut: CachedTexture,
+}
+
+pub(super) fn prepare_atmosphere_textures(
+    views: Query<(Entity, &AtmosphereSettings), With<Atmosphere>>,
+    render_device: Res<RenderDevice>,
+    mut texture_cache: ResMut<TextureCache>,
+    mut commands: Commands,
+) {
+    for (entity, lut_settings) in &views {
+        let transmittance_lut = texture_cache.get(
+            &render_device,
+            TextureDescriptor {
+                label: Some("transmittance_lut"),
+                size: Extent3d {
+                    width: lut_settings.transmittance_lut_size.x,
+                    height: lut_settings.transmittance_lut_size.y,
+                    depth_or_array_layers: 1,
+                },
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: TextureDimension::D2,
+                format: TextureFormat::Rgba16Float,
+                usage: TextureUsages::STORAGE_BINDING | TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            },
+        );
+
+        let multiscattering_lut = texture_cache.get(
+            &render_device,
+            TextureDescriptor {
+                label: Some("multiscattering_lut"),
+                size: Extent3d {
+                    width: lut_settings.multiscattering_lut_size.x,
+                    height: lut_settings.multiscattering_lut_size.y,
+                    depth_or_array_layers: 1,
+                },
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: TextureDimension::D2,
+                format: TextureFormat::Rgba16Float,
+                usage: TextureUsages::STORAGE_BINDING | TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            },
+        );
+
+        let sky_view_lut = texture_cache.get(
+            &render_device,
+            TextureDescriptor {
+                label: Some("sky_view_lut"),
+                size: Extent3d {
+                    width: lut_settings.sky_view_lut_size.x,
+                    height: lut_settings.sky_view_lut_size.y,
+                    depth_or_array_layers: 1,
+                },
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: TextureDimension::D2,
+                format: TextureFormat::Rgba16Float,
+                usage: TextureUsages::STORAGE_BINDING | TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            },
+        );
+
+        let aerial_view_lut = texture_cache.get(
+            &render_device,
+            TextureDescriptor {
+                label: Some("aerial_view_lut"),
+                size: Extent3d {
+                    width: lut_settings.aerial_view_lut_size.x,
+                    height: lut_settings.aerial_view_lut_size.y,
+                    depth_or_array_layers: lut_settings.aerial_view_lut_size.z,
+                },
+                mip_level_count: 1,
+                sample_count: 1,
+                dimension: TextureDimension::D3,
+                format: TextureFormat::Rgba16Float,
+                usage: TextureUsages::STORAGE_BINDING | TextureUsages::TEXTURE_BINDING,
+                view_formats: &[],
+            },
+        );
+
+        commands.entity(entity).insert({
+            AtmosphereTextures {
+                transmittance_lut,
+                multiscattering_lut,
+                sky_view_lut,
+                aerial_view_lut,
+            }
+        });
+    }
+}
+
+#[derive(Resource, Default)]
+pub struct AtmosphereTransforms {
+    uniforms: DynamicUniformBuffer<AtmosphereTransform>,
+}
+
+impl AtmosphereTransforms {
+    #[inline]
+    pub fn uniforms(&self) -> &DynamicUniformBuffer<AtmosphereTransform> {
+        &self.uniforms
+    }
+}
+
+#[derive(ShaderType)]
+pub struct AtmosphereTransform {
+    world_from_atmosphere: Mat4,
+    atmosphere_from_world: Mat4,
+}
+
+#[derive(Component)]
+pub struct AtmosphereTransformsOffset {
+    index: u32,
+}
+
+impl AtmosphereTransformsOffset {
+    #[inline]
+    pub fn index(&self) -> u32 {
+        self.index
+    }
+}
+
+pub(super) fn prepare_atmosphere_transforms(
+    views: Query<(Entity, &ExtractedView), (With<Atmosphere>, With<Camera3d>)>,
+    render_device: Res<RenderDevice>,
+    render_queue: Res<RenderQueue>,
+    mut atmo_uniforms: ResMut<AtmosphereTransforms>,
+    mut commands: Commands,
+) {
+    let atmo_count = views.iter().len();
+    let Some(mut writer) =
+        atmo_uniforms
+            .uniforms
+            .get_writer(atmo_count, &render_device, &render_queue)
+    else {
+        return;
+    };
+
+    for (entity, view) in &views {
+        let world_from_view = view.world_from_view.compute_matrix();
+        let camera_z = world_from_view.z_axis.truncate();
+        let camera_y = world_from_view.y_axis.truncate();
+        let atmo_z = camera_z
+            .with_y(0.0)
+            .try_normalize()
+            .unwrap_or_else(|| camera_y.with_y(0.0).normalize());
+        let atmo_y = Vec3::Y;
+        let atmo_x = atmo_y.cross(atmo_z).normalize();
+        let world_from_atmosphere = Mat4::from_cols(
+            atmo_x.extend(0.0),
+            atmo_y.extend(0.0),
+            atmo_z.extend(0.0),
+            world_from_view.w_axis,
+        );
+
+        let atmosphere_from_world = world_from_atmosphere.inverse();
+
+        commands.entity(entity).insert(AtmosphereTransformsOffset {
+            index: writer.write(&AtmosphereTransform {
+                world_from_atmosphere,
+                atmosphere_from_world,
+            }),
+        });
+    }
+}
+
+#[derive(Component)]
+pub(crate) struct AtmosphereBindGroups {
+    pub transmittance_lut: BindGroup,
+    pub multiscattering_lut: BindGroup,
+    pub sky_view_lut: BindGroup,
+    pub aerial_view_lut: BindGroup,
+    pub render_sky: BindGroup,
+}
+
+pub(super) fn prepare_atmosphere_bind_groups(
+    views: Query<
+        (Entity, &AtmosphereTextures, &ViewDepthTexture, &Msaa),
+        (With<Camera3d>, With<Atmosphere>),
+    >,
+    render_device: Res<RenderDevice>,
+    layouts: Res<AtmosphereBindGroupLayouts>,
+    render_sky_layouts: Res<RenderSkyBindGroupLayouts>,
+    samplers: Res<AtmosphereSamplers>,
+    view_uniforms: Res<ViewUniforms>,
+    lights_uniforms: Res<LightMeta>,
+    atmosphere_transforms: Res<AtmosphereTransforms>,
+    atmosphere_uniforms: Res<ComponentUniforms<Atmosphere>>,
+    settings_uniforms: Res<ComponentUniforms<AtmosphereSettings>>,
+
+    mut commands: Commands,
+) {
+    if views.iter().len() == 0 {
+        return;
+    }
+
+    let atmosphere_binding = atmosphere_uniforms
+        .binding()
+        .expect("Failed to prepare atmosphere bind groups. Atmosphere uniform buffer missing");
+
+    let transforms_binding = atmosphere_transforms
+        .uniforms()
+        .binding()
+        .expect("Failed to prepare atmosphere bind groups. Atmosphere transforms buffer missing");
+
+    let settings_binding = settings_uniforms.binding().expect(
+        "Failed to prepare atmosphere bind groups. AtmosphereSettings uniform buffer missing",
+    );
+
+    let view_binding = view_uniforms
+        .uniforms
+        .binding()
+        .expect("Failed to prepare atmosphere bind groups. View uniform buffer missing");
+
+    let lights_binding = lights_uniforms
+        .view_gpu_lights
+        .binding()
+        .expect("Failed to prepare atmosphere bind groups. Lights uniform buffer missing");
+
+    for (entity, textures, view_depth_texture, msaa) in &views {
+        let transmittance_lut = render_device.create_bind_group(
+            "transmittance_lut_bind_group",
+            &layouts.transmittance_lut,
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (13, &textures.transmittance_lut.default_view),
+            )),
+        );
+
+        let multiscattering_lut = render_device.create_bind_group(
+            "multiscattering_lut_bind_group",
+            &layouts.multiscattering_lut,
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (5, &textures.transmittance_lut.default_view),
+                (6, &samplers.transmittance_lut),
+                (13, &textures.multiscattering_lut.default_view),
+            )),
+        );
+
+        let sky_view_lut = render_device.create_bind_group(
+            "sky_view_lut_bind_group",
+            &layouts.sky_view_lut,
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (2, transforms_binding.clone()),
+                (3, view_binding.clone()),
+                (4, lights_binding.clone()),
+                (5, &textures.transmittance_lut.default_view),
+                (6, &samplers.transmittance_lut),
+                (7, &textures.multiscattering_lut.default_view),
+                (8, &samplers.multiscattering_lut),
+                (13, &textures.sky_view_lut.default_view),
+            )),
+        );
+
+        let aerial_view_lut = render_device.create_bind_group(
+            "sky_view_lut_bind_group",
+            &layouts.aerial_view_lut,
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (3, view_binding.clone()),
+                (4, lights_binding.clone()),
+                (5, &textures.transmittance_lut.default_view),
+                (6, &samplers.transmittance_lut),
+                (7, &textures.multiscattering_lut.default_view),
+                (8, &samplers.multiscattering_lut),
+                (13, &textures.aerial_view_lut.default_view),
+            )),
+        );
+
+        let render_sky = render_device.create_bind_group(
+            "render_sky_bind_group",
+            if *msaa == Msaa::Off {
+                &render_sky_layouts.render_sky
+            } else {
+                &render_sky_layouts.render_sky_msaa
+            },
+            &BindGroupEntries::with_indices((
+                (0, atmosphere_binding.clone()),
+                (1, settings_binding.clone()),
+                (2, transforms_binding.clone()),
+                (3, view_binding.clone()),
+                (4, lights_binding.clone()),
+                (5, &textures.transmittance_lut.default_view),
+                (6, &samplers.transmittance_lut),
+                (9, &textures.sky_view_lut.default_view),
+                (10, &samplers.sky_view_lut),
+                (11, &textures.aerial_view_lut.default_view),
+                (12, &samplers.aerial_view_lut),
+                (13, view_depth_texture.view()),
+            )),
+        );
+
+        commands.entity(entity).insert(AtmosphereBindGroups {
+            transmittance_lut,
+            multiscattering_lut,
+            sky_view_lut,
+            aerial_view_lut,
+            render_sky,
+        });
+    }
+}

crates/bevy_pbr/src/atmosphere/sky_view_lut.wgsl

@@ -0,0 +1,78 @@
+#import bevy_pbr::{
+    mesh_view_types::Lights,
+    atmosphere::{
+        types::{Atmosphere, AtmosphereSettings},
+        bindings::{atmosphere, view, settings},
+        functions::{
+            sample_atmosphere, get_local_up, AtmosphereSample,
+            sample_local_inscattering, get_local_r, view_radius,
+            max_atmosphere_distance, direction_atmosphere_to_world,
+            sky_view_lut_uv_to_zenith_azimuth, zenith_azimuth_to_ray_dir,
+        },
+    }
+}
+
+#import bevy_render::{
+    view::View,
+    maths::HALF_PI,
+}
+#import bevy_core_pipeline::fullscreen_vertex_shader::FullscreenVertexOutput
+
+@group(0) @binding(13) var sky_view_lut_out: texture_storage_2d<rgba16float, write>;
+
+@compute
+@workgroup_size(16, 16, 1)
+fn main(@builtin(global_invocation_id) idx: vec3<u32>) {
+    let uv = vec2<f32>(idx.xy) / vec2<f32>(settings.sky_view_lut_size);
+
+    let r = view_radius();
+    var zenith_azimuth = sky_view_lut_uv_to_zenith_azimuth(r, uv);
+
+    let ray_dir_as = zenith_azimuth_to_ray_dir(zenith_azimuth.x, zenith_azimuth.y);
+    let ray_dir_ws = direction_atmosphere_to_world(ray_dir_as);
+
+    let mu = ray_dir_ws.y;
+    let t_max = max_atmosphere_distance(r, mu);
+
+    // Raymarch with quadratic distribution
+    let sample_count = mix(1.0, f32(settings.sky_view_lut_samples), clamp(t_max * 0.01, 0.0, 1.0));
+    let sample_count_floor = floor(sample_count);
+    let t_max_floor = t_max * sample_count_floor / sample_count;
+    var total_inscattering = vec3(0.0);
+    var throughput = vec3(1.0);
+    for (var s = 0.0; s < sample_count; s += 1.0) {
+        // Use quadratic distribution like reference
+        var t0 = (s / sample_count_floor);
+        var t1 = ((s + 1.0) / sample_count_floor);
+        t0 = t0 * t0;
+        t1 = t1 * t1;
+        t1 = select(t_max_floor * t1, t_max, t1 > 1.0);
+        let t_i = t_max_floor * t0 + (t1 - t_max_floor * t0) * 0.3;
+        let dt_i = t1 - t_max_floor * t0;
+
+        let local_r = get_local_r(r, mu, t_i);
+        let local_up = get_local_up(r, t_i, ray_dir_ws);
+        let local_atmosphere = sample_atmosphere(local_r);
+
+        let sample_optical_depth = local_atmosphere.extinction * dt_i;
+        let sample_transmittance = exp(-sample_optical_depth);
+
+        let inscattering = sample_local_inscattering(
+            local_atmosphere,
+            ray_dir_ws,
+            local_r,
+            local_up
+        );
+
+        // Analytical integration of the single scattering term in the radiance transfer equation
+        let s_int = (inscattering - inscattering * sample_transmittance) / local_atmosphere.extinction;
+        total_inscattering += throughput * s_int;
+
+        throughput *= sample_transmittance;
+        if all(throughput < vec3(0.001)) {
+            break;
+        }
+    }
+
+    textureStore(sky_view_lut_out, idx.xy, vec4(total_inscattering, 1.0));
+}

crates/bevy_pbr/src/atmosphere/transmittance_lut.wgsl

@@ -0,0 +1,48 @@
+#import bevy_pbr::atmosphere::{
+    types::{Atmosphere, AtmosphereSettings},
+    bindings::{settings, atmosphere},
+    functions::{AtmosphereSample, sample_atmosphere, get_local_r, max_atmosphere_distance},
+    bruneton_functions::{transmittance_lut_uv_to_r_mu, distance_to_bottom_atmosphere_boundary, distance_to_top_atmosphere_boundary},
+}
+
+
+#import bevy_core_pipeline::fullscreen_vertex_shader::FullscreenVertexOutput
+
+@group(0) @binding(13) var transmittance_lut_out: texture_storage_2d<rgba16float, write>;
+
+@compute 
+@workgroup_size(16, 16, 1)
+fn main(@builtin(global_invocation_id) idx: vec3<u32>) {
+    let uv: vec2<f32> = (vec2<f32>(idx.xy) + 0.5) / vec2<f32>(settings.transmittance_lut_size);
+    // map UV coordinates to view height (r) and zenith cos angle (mu)
+    let r_mu = transmittance_lut_uv_to_r_mu(uv);
+
+    // compute the optical depth from view height r to the top atmosphere boundary
+    let optical_depth = ray_optical_depth(r_mu.x, r_mu.y, settings.transmittance_lut_samples);
+    let transmittance = exp(-optical_depth);
+
+    textureStore(transmittance_lut_out, idx.xy, vec4(transmittance, 1.0));
+}
+
+/// Compute the optical depth of the atmosphere from the ground to the top atmosphere boundary
+/// at a given view height (r) and zenith cos angle (mu)
+fn ray_optical_depth(r: f32, mu: f32, sample_count: u32) -> vec3<f32> {
+    let t_max = max_atmosphere_distance(r, mu);
+    var optical_depth = vec3<f32>(0.0f);
+    var prev_t = 0.0f;
+
+    for (var i = 0u; i < sample_count; i++) {
+        let t_i = t_max * (f32(i) + 0.3f) / f32(sample_count);
+        let dt = t_i - prev_t;
+        prev_t = t_i;
+
+        let r_i = get_local_r(r, mu, t_i);
+
+        let atmosphere_sample = sample_atmosphere(r_i);
+        let sample_optical_depth = atmosphere_sample.extinction * dt;
+
+        optical_depth += sample_optical_depth;
+    }
+
+    return optical_depth;
+}

crates/bevy_pbr/src/atmosphere/types.wgsl

@@ -0,0 +1,45 @@
+#define_import_path bevy_pbr::atmosphere::types
+
+struct Atmosphere {
+    // Radius of the planet
+    bottom_radius: f32, // units: m
+
+    // Radius at which we consider the atmosphere to 'end' for out calculations (from center of planet)
+    top_radius: f32, // units: m
+
+    ground_albedo: vec3<f32>,
+
+    rayleigh_density_exp_scale: f32,
+    rayleigh_scattering: vec3<f32>,
+
+    mie_density_exp_scale: f32,
+    mie_scattering: f32, // units: m^-1
+    mie_absorption: f32, // units: m^-1
+    mie_asymmetry: f32, // the "asymmetry" value of the phase function, unitless. Domain: (-1, 1)
+
+    ozone_layer_altitude: f32, // units: m
+    ozone_layer_width: f32, // units: m
+    ozone_absorption: vec3<f32>, // ozone absorption. units: m^-1
+}
+
+struct AtmosphereSettings {
+    transmittance_lut_size: vec2<u32>,
+    multiscattering_lut_size: vec2<u32>,
+    sky_view_lut_size: vec2<u32>,
+    aerial_view_lut_size: vec3<u32>,
+    transmittance_lut_samples: u32,
+    multiscattering_lut_dirs: u32,
+    multiscattering_lut_samples: u32,
+    sky_view_lut_samples: u32,
+    aerial_view_lut_samples: u32,
+    aerial_view_lut_max_distance: f32,
+    scene_units_to_m: f32,
+}
+
+
+// "Atmosphere space" is just the view position with y=0 and oriented horizontally,
+// so the horizon stays a horizontal line in our luts
+struct AtmosphereTransforms {
+    world_from_atmosphere: mat4x4<f32>,
+    atmosphere_from_world: mat4x4<f32>,
+}

crates/bevy_pbr/src/lib.rs

@@ -24,6 +24,7 @@ pub mod experimental {
     }
 }
 
+mod atmosphere;
 mod cluster;
 mod components;
 pub mod decal;
@@ -48,6 +49,7 @@ use crate::material_bind_groups::FallbackBindlessResources;
 
 use bevy_color::{Color, LinearRgba};
 
+pub use atmosphere::*;
 pub use cluster::*;
 pub use components::*;
 pub use extended_material::*;
@@ -342,6 +344,7 @@ impl Plugin for PbrPlugin {
                 SyncComponentPlugin::<SpotLight>::default(),
                 ExtractComponentPlugin::<AmbientLight>::default(),
             ))
+            .add_plugins(AtmospherePlugin)
             .configure_sets(
                 PostUpdate,
                 (

crates/bevy_pbr/src/light/mod.rs

@@ -86,6 +86,8 @@ pub mod light_consts {
         pub const FULL_DAYLIGHT: f32 = 20_000.;
         /// The amount of light (lux) in direct sunlight.
         pub const DIRECT_SUNLIGHT: f32 = 100_000.;
+        /// The amount of light (lux) of raw sunlight, not filtered by the atmosphere.
+        pub const RAW_SUNLIGHT: f32 = 130_000.;
     }
 }
 

crates/bevy_pbr/src/ssao/ssao.wgsl

@@ -10,7 +10,7 @@
 // Source code base on SSRT3 implementation
 // https://github.com/cdrinmatane/SSRT3
 
-#import bevy_pbr::ssao_utils::fast_acos
+#import bevy_render::maths::fast_acos
 
 #import bevy_render::{
     view::View,

crates/bevy_pbr/src/ssao/ssao_utils.wgsl

@@ -11,14 +11,3 @@ fn ssao_multibounce(visibility: f32, base_color: vec3<f32>) -> vec3<f32> {
     let x = vec3<f32>(visibility);
     return max(x, ((x * a + b) * x + c) * x);
 }
-
-fn fast_sqrt(x: f32) -> f32 {
-    return bitcast<f32>(0x1fbd1df5 + (bitcast<i32>(x) >> 1u));
-}
-
-fn fast_acos(in_x: f32) -> f32 {
-    let x = abs(in_x);
-    var res = -0.156583 * x + HALF_PI;
-    res *= fast_sqrt(1.0 - x);
-    return select(PI - res, res, in_x >= 0.0);
-}

crates/bevy_render/src/maths.wgsl

@@ -99,3 +99,56 @@ fn project_onto(lhs: vec3<f32>, rhs: vec3<f32>) -> vec3<f32> {
     let other_len_sq_rcp = 1.0 / dot(rhs, rhs);
     return rhs * dot(lhs, rhs) * other_len_sq_rcp;
 }
+
+// Below are fast approximations of common irrational and trig functions. These
+// are likely most useful when raymarching, for example, where complete numeric
+// accuracy can be sacrificed for greater sample count.
+
+fn fast_sqrt(x: f32) -> f32 {
+    return bitcast<f32>(0x1fbd1df5 + (bitcast<i32>(x) >> 1u));
+}
+
+// Slightly less accurate than fast_acos_4, but much simpler.
+fn fast_acos(in_x: f32) -> f32 {
+    let x = abs(in_x);
+    var res = -0.156583 * x + HALF_PI;
+    res *= fast_sqrt(1.0 - x);
+    return select(PI - res, res, in_x >= 0.0);
+}
+
+// 4th order polynomial approximation
+// 4 VGRP, 16 ALU Full Rate
+// 7 * 10^-5 radians precision
+// Reference : Handbook of Mathematical Functions (chapter : Elementary Transcendental Functions), M. Abramowitz and I.A. Stegun, Ed.
+fn fast_acos_4(x: f32) -> f32 {
+    let x1 = abs(x);
+    let x2 = x1 * x1;
+    let x3 = x2 * x1;
+    var s: f32;
+
+    s = -0.2121144 * x1 + 1.5707288;
+    s = 0.0742610 * x2 + s;
+    s = -0.0187293 * x3 + s;
+    s = fast_sqrt(1.0 - x1) * s;
+
+	// acos function mirroring
+    return select(PI - s, s, x >= 0.0);
+}
+
+fn fast_atan2(y: f32, x: f32) -> f32 {
+    var t0 = max(abs(x), abs(y));
+    var t1 = min(abs(x), abs(y));
+    var t3 = t1 / t0;
+    var t4 = t3 * t3;
+
+    t0 = 0.0872929;
+    t0 = t0 * t4 - 0.301895;
+    t0 = t0 * t4 + 1.0;
+    t3 = t0 * t3;
+
+    t3 = select(t3, (0.5 * PI) - t3, abs(y) > abs(x));
+    t3 = select(t3, PI - t3, x < 0);
+    t3 = select(-t3, t3, y > 0);
+
+    return t3;
+}

crates/bevy_render/src/render_resource/bind_group_layout_entries.rs

@@ -556,4 +556,16 @@ pub mod binding_types {
         }
         .into_bind_group_layout_entry_builder()
     }
+
+    pub fn texture_storage_3d(
+        format: TextureFormat,
+        access: StorageTextureAccess,
+    ) -> BindGroupLayoutEntryBuilder {
+        BindingType::StorageTexture {
+            access,
+            format,
+            view_dimension: TextureViewDimension::D3,
+        }
+        .into_bind_group_layout_entry_builder()
+    }
 }

examples/3d/atmosphere.rs

@@ -0,0 +1,125 @@
+//! This example showcases pbr atmospheric scattering
+
+use std::f32::consts::PI;
+
+use bevy::{
+    core_pipeline::{bloom::Bloom, tonemapping::Tonemapping},
+    pbr::{light_consts::lux, Atmosphere, AtmosphereSettings, CascadeShadowConfigBuilder},
+    prelude::*,
+    render::camera::Exposure,
+};
+
+fn main() {
+    App::new()
+        .add_plugins(DefaultPlugins)
+        .add_systems(Startup, (setup_camera_fog, setup_terrain_scene))
+        .add_systems(Update, dynamic_scene)
+        .run();
+}
+
+fn setup_camera_fog(mut commands: Commands) {
+    commands.spawn((
+        Camera3d::default(),
+        // HDR is required for atmospheric scattering to be properly applied to the scene
+        Camera {
+            hdr: true,
+            ..default()
+        },
+        Transform::from_xyz(-1.2, 0.15, 0.0).looking_at(Vec3::Y * 0.1, Vec3::Y),
+        // This is the component that enables atmospheric scattering for a camera
+        Atmosphere::EARTH,
+        // The scene is in units of 10km, so we need to scale up the
+        // aerial view lut distance and set the scene scale accordingly.
+        // Most usages of this feature will not need to adjust this.
+        AtmosphereSettings {
+            aerial_view_lut_max_distance: 3.2e5,
+            scene_units_to_m: 1e+4,
+            ..Default::default()
+        },
+        // The directional light illuminance  used in this scene
+        // (the one recommended for use with this feature) is
+        // quite bright, so raising the exposure compensation helps
+        // bring the scene to a nicer brightness range.
+        Exposure::SUNLIGHT,
+        // Tonemapper chosen just because it looked good with the scene, any
+        // tonemapper would be fine :)
+        Tonemapping::AcesFitted,
+        // Bloom gives the sun a much more natural look.
+        Bloom::NATURAL,
+    ));
+}
+
+#[derive(Component)]
+struct Terrain;
+
+fn setup_terrain_scene(
+    mut commands: Commands,
+    mut meshes: ResMut<Assets<Mesh>>,
+    mut materials: ResMut<Assets<StandardMaterial>>,
+    asset_server: Res<AssetServer>,
+) {
+    // Configure a properly scaled cascade shadow map for this scene (defaults are too large, mesh units are in km)
+    let cascade_shadow_config = CascadeShadowConfigBuilder {
+        first_cascade_far_bound: 0.3,
+        maximum_distance: 3.0,
+        ..default()
+    }
+    .build();
+
+    // Sun
+    commands.spawn((
+        DirectionalLight {
+            shadows_enabled: true,
+            // lux::RAW_SUNLIGHT is recommended for use with this feature, since
+            // other values approximate sunlight *post-scattering* in various
+            // conditions. RAW_SUNLIGHT in comparison is the illuminance of the
+            // sun unfiltered by the atmosphere, so it is the proper input for
+            // sunlight to be filtered by the atmosphere.
+            illuminance: lux::RAW_SUNLIGHT,
+            ..default()
+        },
+        Transform::from_xyz(1.0, -0.4, 0.0).looking_at(Vec3::ZERO, Vec3::Y),
+        cascade_shadow_config,
+    ));
+
+    let sphere_mesh = meshes.add(Mesh::from(Sphere { radius: 1.0 }));
+
+    // light probe spheres
+    commands.spawn((
+        Mesh3d(sphere_mesh.clone()),
+        MeshMaterial3d(materials.add(StandardMaterial {
+            base_color: Color::WHITE,
+            metallic: 1.0,
+            perceptual_roughness: 0.0,
+            ..default()
+        })),
+        Transform::from_xyz(-0.3, 0.1, -0.1).with_scale(Vec3::splat(0.05)),
+    ));
+
+    commands.spawn((
+        Mesh3d(sphere_mesh.clone()),
+        MeshMaterial3d(materials.add(StandardMaterial {
+            base_color: Color::WHITE,
+            metallic: 0.0,
+            perceptual_roughness: 1.0,
+            ..default()
+        })),
+        Transform::from_xyz(-0.3, 0.1, 0.1).with_scale(Vec3::splat(0.05)),
+    ));
+
+    // Terrain
+    commands.spawn((
+        Terrain,
+        SceneRoot(
+            asset_server.load(GltfAssetLabel::Scene(0).from_asset("models/terrain/terrain.glb")),
+        ),
+        Transform::from_xyz(-1.0, 0.0, -0.5)
+            .with_scale(Vec3::splat(0.5))
+            .with_rotation(Quat::from_rotation_y(PI / 2.0)),
+    ));
+}
+
+fn dynamic_scene(mut suns: Query<&mut Transform, With<DirectionalLight>>, time: Res<Time>) {
+    suns.iter_mut()
+        .for_each(|mut tf| tf.rotate_x(-time.delta_secs() * PI / 10.0));
+}

examples/README.md

@@ -138,6 +138,7 @@ Example | Description
 [Animated Material](../examples/3d/animated_material.rs) | Shows how to animate material properties
 [Anisotropy](../examples/3d/anisotropy.rs) | Displays an example model with anisotropy
 [Anti-aliasing](../examples/3d/anti_aliasing.rs) | Compares different anti-aliasing methods
+[Atmosphere](../examples/3d/atmosphere.rs) | A scene showcasing pbr atmospheric scattering
 [Atmospheric Fog](../examples/3d/atmospheric_fog.rs) | A scene showcasing the atmospheric fog effect
 [Auto Exposure](../examples/3d/auto_exposure.rs) | A scene showcasing auto exposure
 [Blend Modes](../examples/3d/blend_modes.rs) | Showcases different blend modes