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// Mapping from view height (r) and zenith cos angle (mu) to UV coordinates in the transmittance LUT
64// Assuming r between ground and top atmosphere boundary, and mu= cos(zenith_angle)
65// Chosen to increase precision near the ground and to work around a discontinuity at the horizon
66// See Bruneton and Neyret 2008, "Precomputed Atmospheric Scattering" section 4
67fn transmittance_lut_r_mu_to_uv(atm: Atmosphere, r: f32, mu: f32) -> vec2<f32> {
68  // Distance along a horizontal ray from the ground to the top atmosphere boundary
69    let H = sqrt(atm.top_radius * atm.top_radius - atm.bottom_radius * atm.bottom_radius);
70
71  // Distance from a point at height r to the horizon
72  // ignore the case where r <= atmosphere.bottom_radius
73    let rho = sqrt(max(r * r - atm.bottom_radius * atm.bottom_radius, 0.0));
74
75  // Distance from a point at height r to the top atmosphere boundary at zenith angle mu
76    let d = distance_to_top_atmosphere_boundary(atm, r, mu);
77
78  // Minimum and maximum distance to the top atmosphere boundary from a point at height r
79    let d_min = atm.top_radius - r; // length of the ray straight up to the top atmosphere boundary
80    let d_max = rho + H; // length of the ray to the top atmosphere boundary and grazing the horizon
81
82    let u = (d - d_min) / (d_max - d_min);
83    let v = rho / H;
84    return vec2<f32>(u, v);
85}
86
87// Inverse of the mapping above, mapping from UV coordinates in the transmittance LUT to view height (r) and zenith cos angle (mu)
88fn transmittance_lut_uv_to_r_mu(uv: vec2<f32>) -> vec2<f32> {
89  // Distance to top atmosphere boundary for a horizontal ray at ground level
90    let H = sqrt(atmosphere.top_radius * atmosphere.top_radius - atmosphere.bottom_radius * atmosphere.bottom_radius);
91
92  // Distance to the horizon, from which we can compute r:
93    let rho = H * uv.y;
94    let r = sqrt(rho * rho + atmosphere.bottom_radius * atmosphere.bottom_radius);
95
96  // Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
97  // and maximum values over all mu- obtained for (r,1) and (r,mu_horizon) -
98  // from which we can recover mu:
99    let d_min = atmosphere.top_radius - r;
100    let d_max = rho + H;
101    let d = d_min + uv.x * (d_max - d_min);
102
103    var mu: f32;
104    if d == 0.0 {
105        mu = 1.0;
106    } else {
107        mu = (H * H - rho * rho - d * d) / (2.0 * r * d);
108    }
109
110    mu = clamp(mu, -1.0, 1.0);
111
112    return vec2<f32>(r, mu);
113}
114
115/// Simplified ray-sphere intersection
116/// where:
117/// Ray origin, o = [0,0,r] with r <= atmosphere.top_radius
118/// mu is the cosine of spherical coordinate theta (-1.0 <= mu <= 1.0)
119/// so ray direction in spherical coordinates is [1,acos(mu),0] which needs to be converted to cartesian
120/// Direction of ray, u = [0,sqrt(1-mu*mu),mu]
121/// Center of sphere, c = [0,0,0]
122/// Radius of sphere, r = atmosphere.top_radius
123/// This function solves the quadratic equation for line-sphere intersection simplified under these assumptions
124fn distance_to_top_atmosphere_boundary(atm: Atmosphere, r: f32, mu: f32) -> f32 {
125  // ignore the case where r > atm.top_radius
126    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atm.top_radius * atm.top_radius, 0.0);
127    return max(-r * mu + sqrt(positive_discriminant), 0.0);
128}
129
130/// Simplified ray-sphere intersection
131/// as above for intersections with the ground
132fn distance_to_bottom_atmosphere_boundary(r: f32, mu: f32) -> f32 {
133    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atmosphere.bottom_radius * atmosphere.bottom_radius, 0.0);
134    return max(-r * mu - sqrt(positive_discriminant), 0.0);
135}
136
137fn ray_intersects_ground(r: f32, mu: f32) -> bool {
138    return mu < 0.0 && r * r * (mu * mu - 1.0) + atmosphere.bottom_radius * atmosphere.bottom_radius >= 0.0;
139}