Procedural atmospheric scattering

Cargo.toml

@@ -114,20 +114,18 @@ default = [
   "android_shared_stdcxx",
   "animation",
   "bevy_asset",
+  "bevy_state",
   "bevy_audio",
   "bevy_color",
-  "bevy_core_pipeline",
   "bevy_gilrs",
-  "bevy_gizmos",
   "bevy_gltf",
-  "bevy_mesh_picking_backend",
+  "bevy_scene",
+  "bevy_winit",
+  "bevy_core_pipeline",
   "bevy_pbr",
   "bevy_picking",
   "bevy_render",
-  "bevy_scene",
   "bevy_sprite",
-  "bevy_sprite_picking_backend",
-  "bevy_state",
   "bevy_text",
   "bevy_ui",
   "bevy_ui_picking_backend",
@@ -138,12 +136,16 @@ default = [
   "hdr",
   "multi_threaded",
   "png",
-  "smaa_luts",
-  "sysinfo_plugin",
-  "tonemapping_luts",
+  "hdr",
   "vorbis",
-  "webgl2",
   "x11",
+  "bevy_gizmos",
+  "android_shared_stdcxx",
+  "tonemapping_luts",
+  "smaa_luts",
+  "default_font",
+  "webgl2",
+  "sysinfo_plugin",
 ]
 
 # Provides an implementation for picking meshes
@@ -862,6 +864,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,80 @@
+#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, position_ndc_to_world, depth_ndc_to_view_z,
+            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) //TODO: this approach makes it so closer slices get fewer samples. But we also expect those to have less scattering. So win/win?
+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 = max_atmosphere_distance(r, mu);
+
+    var prev_t = 0.0;
+    var total_inscattering = vec3(0.0);
+    var optical_depth = vec3(0.0);
+    var throughput = vec3(1.0);
+
+    // The aerial view LUT is in NDC space, so it uses bevy's reverse z convention. Since
+    // we write multiple slices from each thread, we need to iterate in order near->far, which 
+    // is why the indices are reversed.
+    for (var slice_i: i32 = i32(settings.aerial_view_lut_size.z - 1); slice_i >= 0; slice_i--) {
+        var sum_transmittance = 0.0;
+        for (var step_i: i32 = i32(settings.aerial_view_lut_samples - 1); step_i >= 0; step_i--) {
+            let sample_depth = depth_at_sample(slice_i, step_i);
+            //view_dir.w is the cosine of the angle between the view vector and the camera forward vector, used to correct the step length.
+            let t_i = -depth_ndc_to_view_z(sample_depth) / ray_dir.w * settings.scene_units_to_km;
+            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);
+            optical_depth += sample_optical_depth;
+
+            // use beer's law to get transmittance from optical density
+            let transmittance_to_sample = exp(-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;
+            }
+
+            sum_transmittance += transmittance_to_sample.r + transmittance_to_sample.g + transmittance_to_sample.b;
+        }
+        //We only have one channel to store transmittance, so we store the mean 
+        let mean_transmittance = sum_transmittance / (f32(settings.aerial_view_lut_samples) * 3.0);
+        textureStore(aerial_view_lut_out, vec3(vec2<i32>(idx.xy), slice_i), vec4(total_inscattering, mean_transmittance));
+    }
+}
+
+// linearly interpolates from 0..1 on the domain of slice_i, using step_i as a substep index
+fn depth_at_sample(slice_i: i32, step_i: i32) -> f32 {
+    return (f32(slice_i) + ((f32(step_i) + 0.5) / f32(settings.aerial_view_lut_samples))) / f32(settings.aerial_view_lut_size.z);
+}

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;
+}