Generalized atmospheric scattering media

[ecoskey]

Objective

Right now, only a very small set of atmospheres are possible with Bevy's atmospheric scattering system because of the fixed set of scattering terms available. For example, a scene set in a dry, desert environment might want a dense low-lying layer of dust particulate separate from the ambient dust in the atmosphere. This PR introduces a mechanism for generalized scattering media, replacing the fixed scattering terms with a customizable asset.

Solution

#[derive(TypePath, Asset, Clone)]
pub struct ScatteringMedium {
    pub label: Option<Cow<'static, str>>,
    pub falloff_resolution: u32,
    pub phase_resolution: u32,
    pub terms: SmallVec<[ScatteringTerm; 1]>,
}

see other new types

#[derive(Default, Clone)]
pub struct ScatteringTerm {
    pub absorption: Vec3,
    pub scattering: Vec3,
    pub falloff: Falloff,
    pub phase: PhaseFunction,
}

#[derive(Default, Clone)]
pub enum Falloff {
    #[default]
    Linear,
    Exponential {
        scale: f32,
    },
    Tent {
        center: f32,
        width: f32,
    },
    /// A falloff function defined by a custom curve.
    ///
    /// domain: [0, 1),
    /// range: [0, 1],
    Curve(Arc<dyn Curve<f32> + Send + Sync>),
}

#[derive(Clone)]
pub enum PhaseFunction {
    Isotropic,
    Rayleigh,
    Mie {
        /// domain: [-1, 1]
        bias: f32,
    },
    /// A phase function defined by a custom curve.
    ///
    /// domain: [-1, 1]
    /// range: [0, 1]
    Curve(Arc<dyn Curve<f32> + Send + Sync>),
}

ScatteringMedium contains a list of ScatteringTerms, which are processed into a set of two LUTs:

• The "density LUT", a 2D falloff_resolution x 2 LUT which contains the medium's optical density with respect to the atmosphere's "falloff parameter", a linear value which is 1.0 at the planet's surface and 0.0 at the edge of space. Absorption density and scattering density correspond to the first and second rows respectively.
• The "scattering LUT", a 2D falloff_resolution x phase_resolution LUT which contains the medium's scattering density multiplied by the phase function, with the U axis corresponding to the falloff parameter and the V axis corresponding to neg_LdotV * 0.5 + 0.5, where neg_LdotV is the dot product of the light direction and the outgoing view vector.

Testing

• Need to verify output, should be almost exactly the same
• exponential falloff is slightly different now, but verified new parameters against the old in Desmos.

TODOS:

• Docs Docs Docs
• Cleanup
• profile perf
• reduce memory usage/traffic. This approach requires a few extra texture samples in the inner loop of each pass. Each atmosphere LUT is still quite small, but texture samples are expensive and the new LUTs use f32 texels currently.

Showcase

Click to view showcase

fn init_atmosphere(mut commands: Commands, scattering_media: ResMut<Assets<ScatteringMedium>>) {
  let earth_atmosphere = scattering_media.add(
    ScatteringMedium::new(
      256,
      256,
      [
          // rayleigh scattering
          ScatteringTerm {
              absorption: Vec3::ZERO,
              scattering: Vec3::new(5.802e-6, 13.558e-6, 33.100e-6),
              falloff: Falloff::Exponential { scale: 12.5 },
              phase: PhaseFunction::Rayleigh,
          },
          // mie scattering
          ScatteringTerm {
              absorption: Vec3::splat(3.996e-6),
              scattering: Vec3::splat(0.444e-6),
              falloff: Falloff::Exponential { scale: 83.5 },
              phase: PhaseFunction::Mie { bias: 0.8 },
          },
          // ozone
          ScatteringTerm {
              absorption: Vec3::new(0.650e-6, 1.881e-6, 0.085e-6),
              scattering: Vec3::ZERO,
              falloff: Falloff::Tent {
                  center: 0.75,
                  width: 0.3,
              },
              phase: PhaseFunction::Isotropic,
          },
      ],
  ));

  commands.spawn((
    Camera3d::default(), 
    Atmosphere {
      bottom_radius: 6_360_000.0,
      top_radius: 6_460_000.0,
      ground_albedo: Vec3::splat(0.3),
      medium: earth_atmosphere,
    },
  ));
}

[github-actions]

It looks like your PR has been selected for a highlight in the next release blog post, but you didn't provide a release note.

Please review the instructions for writing release notes, then expand or revise the content in the release notes directory to showcase your changes.

[aevyrie]

Consider this an approval for architecture/direction. This is not a review of the implementation, which I won't have time for.

[mate-h]

This is some amazing work! I can't wait to create more types scenes with this new parameterization. I checked out the code locally and started playing around with the new parameterization, and what i noticed is that the sky turns completely black if there is only 1 mie scatterer or no scatterers at all. I posted my example file on this branch:
https://github.com/mate-h/bevy/tree/m/general_scattering_test

Also created this quick chart/ visualization for the exponential falloff: https://observablehq.com/@mateh/exponential-falloff

[mate-h]

Posted a PR to this branch to fix some NaNs in the LUTs in case of single mie scatterer or no scatterers.
ecoskey#8

[mate-h]

The diff is really messed up. Could you sync and merge from main again?

[mate-h]

I think the general architectural decision to go this route makes sense, but I do want to make the user facing API concise as before. we can revise the naming and the user facing API in later PRs. I gave this a thorough test and found that maybe the strength parameter is counterintuitive and I prefer to use the scale_height as described in the paper.
I created a branch as a demo: ecoskey#9
no need to merge this but sort of encapsulates my feedback on the parameterization. the remaining comment I left also also non-blocking.

[mate-h]

Agreed that doing the switch to with_indices is warranted in this case, so let's just stick to it for now until we re-organize the binding in a dedicated PR.

[mate-h]

Since scattering density and the absorption density uses the same distribution, I think it would make more sense to simply store the distribution of the main function (in this case exponential) and then multiplying by the absorption or scattering coefficients by passing them as uniforms to the shader. unless we want to use another compute shader to compute the scattering and absorption coefficients independently which sounds like a niche use case.
After thinking about this some more it would definitely make the code more "messy". So even though baking these values into the lut may not make sense for the use case, the current approach still makes for cleaner code. Also since these luts are super small, wouldn't be saving much VRAM anyway if we reduce the number of rows, and we don't need to worry too much about hitting floating point precision limits. (hopefully f16 texture can store these small values on all devices) .

[ecoskey]

I'm not sure things commute properly to be able to factor that out, at least when there's multiple scattering terms present:

density(p) = d1 * f1(p) + d2 * f2(p) + ... + dN * fN(p) where dN are the optical density coefficients and fN are the falloff distributions, which can be arbitrary functions

Something slightly different that we can do is after calculating all the values, normalize them wrt the greatest density value, and store RGBA16Unorm values in each texel instead

[mate-h]

Oh I see so it's a combination of functions for each scattering term. Okay then never mind it would be hard to factor out any coefficients if the functions can be any density distribution function.

[mate-h]

The diffs here is very ugly. I think there has to be a better to manage bindings in version control and code review

[ecoskey]

going forward I'll try to factor out different chunks of the bindings into functions, after we merge #21625

[mate-h]

Also a dual-lobe Mie scattering function is common for clouds. but can be added later. ref: https://research.nvidia.com/labs/rtr/approximate-mie/publications/approximate-mie.pdf

[mate-h]

Maybe it's better to create two separate functions sample_absorption and sample_density?