1// Copyright (c) 2017 Eric Bruneton
2// All rights reserved.
3//
4// Redistribution and use in source and binary forms, with or without
5// modification, are permitted provided that the following conditions
6// are met:
7// 1. Redistributions of source code must retain the above copyright
8//    notice, this list of conditions and the following disclaimer.
9// 2. Redistributions in binary form must reproduce the above copyright
10//    notice, this list of conditions and the following disclaimer in the
11//    documentation and/or other materials provided with the distribution.
12// 3. Neither the name of the copyright holders nor the names of its
13//    contributors may be used to endorse or promote products derived from
14//    this software without specific prior written permission.
15//
16// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
20// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
26// THE POSSIBILITY OF SUCH DAMAGE.
27//
28// Precomputed Atmospheric Scattering
29// Copyright (c) 2008 INRIA
30// All rights reserved.
31//
32// Redistribution and use in source and binary forms, with or without
33// modification, are permitted provided that the following conditions
34// are met:
35// 1. Redistributions of source code must retain the above copyright
36//    notice, this list of conditions and the following disclaimer.
37// 2. Redistributions in binary form must reproduce the above copyright
38//    notice, this list of conditions and the following disclaimer in the
39//    documentation and/or other materials provided with the distribution.
40// 3. Neither the name of the copyright holders nor the names of its
41//    contributors may be used to endorse or promote products derived from
42//    this software without specific prior written permission.
43//
44// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
45// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
46// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
47// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
48// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
49// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
50// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
51// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
52// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
53// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
54// THE POSSIBILITY OF SUCH DAMAGE.
55
56#define_import_path bevy_pbr::atmosphere::bruneton_functions
57
58#import bevy_pbr::atmosphere::{
59    types::Atmosphere,
60    bindings::atmosphere,
61}
62
63// Bruneton's paper calls the surface and exosphere "bottom radius" and "top radius".
64// Our `Atmosphere` struct uses `inner_radius` and `outer_radius` for the same values so naming stays
65// meaningful when the atmosphere entity has an arbitrary transform.
66
67// Mapping from view height (r) and zenith cos angle (mu) to UV coordinates in the transmittance LUT
68// Assuming r between ground and top atmosphere boundary, and mu= cos(zenith_angle)
69// Chosen to increase precision near the ground and to work around a discontinuity at the horizon
70// See Bruneton and Neyret 2008, "Precomputed Atmospheric Scattering" section 4
71fn transmittance_lut_r_mu_to_uv(atm: Atmosphere, r: f32, mu: f32) -> vec2<f32> {
72    // Distance along a horizontal ray from the ground to the top atmosphere boundary
73    let H = sqrt(atm.outer_radius * atm.outer_radius - atm.inner_radius * atm.inner_radius);
74
75    // Distance from a point at height r to the horizon
76    // ignore the case where r <= atmosphere.inner_radius
77    let rho = sqrt(max(r * r - atm.inner_radius * atm.inner_radius, 0.0));
78
79    // Distance from a point at height r to the top atmosphere boundary at zenith angle mu
80    let d = distance_to_top_atmosphere_boundary(atm, r, mu);
81
82    // Minimum and maximum distance to the top atmosphere boundary from a point at height r
83    let d_min = atm.outer_radius - r; // length of the ray straight up to the top atmosphere boundary
84    let d_max = rho + H; // length of the ray to the top atmosphere boundary and grazing the horizon
85
86    let u = (d - d_min) / (d_max - d_min);
87    let v = rho / H;
88    return vec2<f32>(u, v);
89}
90
91// Inverse of the mapping above, mapping from UV coordinates in the transmittance LUT to view height (r) and zenith cos angle (mu)
92fn transmittance_lut_uv_to_r_mu(uv: vec2<f32>) -> vec2<f32> {
93    // Distance to top atmosphere boundary for a horizontal ray at ground level
94    let H = sqrt(atmosphere.outer_radius * atmosphere.outer_radius - atmosphere.inner_radius * atmosphere.inner_radius);
95
96    // Distance to the horizon, from which we can compute r:
97    let rho = H * uv.y;
98    let r = sqrt(rho * rho + atmosphere.inner_radius * atmosphere.inner_radius);
99
100    // Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
101    // and maximum values over all mu- obtained for (r,1) and (r,mu_horizon) -
102    // from which we can recover mu:
103    let d_min = atmosphere.outer_radius - r;
104    let d_max = rho + H;
105    let d = d_min + uv.x * (d_max - d_min);
106
107    var mu: f32;
108    if d == 0.0 {
109        mu = 1.0;
110    } else {
111        mu = (H * H - rho * rho - d * d) / (2.0 * r * d);
112    }
113
114    mu = clamp(mu, -1.0, 1.0);
115
116    return vec2<f32>(r, mu);
117}
118
119/// Simplified ray-sphere intersection
120/// where:
121/// Ray origin, o = [0,0,r] with r <= atmosphere.outer_radius
122/// mu is the cosine of spherical coordinate theta (-1.0 <= mu <= 1.0)
123/// so ray direction in spherical coordinates is [1,acos(mu),0] which needs to be converted to cartesian
124/// Direction of ray, u = [0,sqrt(1-mu*mu),mu]
125/// Center of sphere, c = [0,0,0]
126/// Radius of sphere, r = atmosphere.outer_radius
127/// This function solves the quadratic equation for line-sphere intersection simplified under these assumptions
128fn distance_to_top_atmosphere_boundary(atm: Atmosphere, r: f32, mu: f32) -> f32 {
129    // ignore the case where r > atm.outer_radius
130    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atm.outer_radius * atm.outer_radius, 0.0);
131    return max(-r * mu + sqrt(positive_discriminant), 0.0);
132}
133
134/// Simplified ray-sphere intersection
135/// as above for intersections with the ground
136fn distance_to_bottom_atmosphere_boundary(r: f32, mu: f32) -> f32 {
137    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atmosphere.inner_radius * atmosphere.inner_radius, 0.0);
138    return max(-r * mu - sqrt(positive_discriminant), 0.0);
139}
140
141fn ray_intersects_ground(r: f32, mu: f32) -> bool {
142    return mu < 0.0 && r * r * (mu * mu - 1.0) + atmosphere.inner_radius * atmosphere.inner_radius >= 0.0;
143}