Procedural atmospheric scattering

Cargo.toml

@@ -102,41 +102,34 @@ unused_qualifications = "warn"
 
 [features]
 default = [
-  "android-game-activity",
-  "android-game-activity",
-  "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_remote",
+  "bevy_gltf",
   "bevy_render",
-  "bevy_scene",
   "bevy_sprite",
-  "bevy_sprite_picking_backend",
-  "bevy_state",
   "bevy_text",
   "bevy_ui",
-  "bevy_ui_picking_backend",
-  "bevy_winit",
-  "custom_cursor",
-  "default_font",
-  "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
@@ -202,7 +195,6 @@ bevy_sprite = [
   "bevy_render",
   "bevy_core_pipeline",
   "bevy_color",
-  "bevy_sprite_picking_backend",
 ]
 
 # Provides text functionality
@@ -215,7 +207,6 @@ bevy_ui = [
   "bevy_text",
   "bevy_sprite",
   "bevy_color",
-  "bevy_ui_picking_backend",
 ]
 
 # winit window and input backend
@@ -839,6 +830,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"
@@ -1219,7 +1221,7 @@ setup = [
     "curl",
     "-o",
     "assets/models/bunny.meshlet_mesh",
-    "https://raw.githubusercontent.com/JMS55/bevy_meshlet_asset/defbd9b32072624d40d57de7d345c66a9edf5d0b/bunny.meshlet_mesh",
+    "https://raw.githubusercontent.com/JMS55/bevy_meshlet_asset/b6c712cfc87c65de419f856845401aba336a7bcd/bunny.meshlet_mesh",
   ],
 ]
 

crates/bevy_pbr/src/atmosphere/aerial_view_lut.wgsl

@@ -0,0 +1,90 @@
+#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
+        },
+        bruneton_functions::{distance_to_top_atmosphere_boundary, distance_to_bottom_atmosphere_boundary,ray_intersects_ground}
+    }
+}
+
+
+@group(0) @binding(12) 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); //TODO: negate for lighting calcs?
+    let r = view_radius();
+    let mu = ray_dir.y;
+
+    var prev_t = 0.0;
+    var total_inscattering = vec3(0.0);
+    var optical_depth = vec3(0.0);
+    for (var slice_i: i32 = i32(settings.aerial_view_lut_size.z - 1); slice_i >= 0; slice_i--) { //reversed loop to iterate raw depth values near->far 
+        var sum_transmittance = 0.0;
+        for (var step_i: i32 = i32(settings.aerial_view_lut_samples - 1); step_i >= 0; step_i--) { //same here
+            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 step_length = (t_i - prev_t);
+            prev_t = t_i;
+
+            let local_r = get_local_r(r, mu, t_i);
+            if local_r > atmosphere.top_radius { break; }
+            let local_up = get_local_up(r, t_i, ray_dir.xyz);
+
+            let local_atmosphere = sample_atmosphere(local_r);
+            optical_depth += local_atmosphere.extinction * step_length; //TODO: units between step_length and atmosphere
+
+            let transmittance_to_sample = exp(-optical_depth);
+
+            var local_inscattering = sample_local_inscattering(local_atmosphere, transmittance_to_sample, ray_dir.xyz, local_r, local_up);
+            total_inscattering += local_inscattering * step_length;
+            sum_transmittance += transmittance_to_sample.r + transmittance_to_sample.g + transmittance_to_sample.b;
+        }
+        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));
+    }
+}
+
+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);
+}
+
+
+//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; //TODO: remove this step and just use position_ndc_to_world? we didn't need to transform in view space
+
+    // 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); //TODO: correct sign?
+}

crates/bevy_pbr/src/atmosphere/bindings.wgsl

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