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 <[email protected]>
Co-authored-by: atlv <[email protected]>
Co-authored-by: JMS55 <[email protected]>
Co-authored-by: Emerson Coskey <[email protected]>
Co-authored-by: Máté Homolya <[email protected]>

Author: 5 people

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"

assets/models/terrain/terrain.glb

@@ -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/aerial_view_lut.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/bindings.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;
+}