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
56import package::atmosphere::{
57    types::Atmosphere,
58    bindings::atmosphere,
59    functions::{get_local_r, sample_transmittance_lut},
60};
61
62// Bruneton's paper calls the surface and exosphere "bottom radius" and "top radius".
63// Our `Atmosphere` struct uses `inner_radius` and `outer_radius` for the same values so naming stays
64// meaningful when the atmosphere entity has an arbitrary transform.
65
66/// Mapping from view height (r) and zenith cos angle (mu) to UV coordinates in the transmittance LUT
67/// Assuming r between ground and top atmosphere boundary, and mu= cos(zenith_angle)
68/// Chosen to increase precision near the ground and to work around a discontinuity at the horizon
69/// See Bruneton and Neyret 2008, "Precomputed Atmospheric Scattering" section 4
70fn transmittance_lut_r_mu_to_uv(atm: Atmosphere, r: f32, mu: f32) -> vec2<f32> {
71    // Distance along a horizontal ray from the ground to the top atmosphere boundary
72    let H = sqrt(atm.outer_radius * atm.outer_radius - atm.inner_radius * atm.inner_radius);
73
74    // Distance from a point at height r to the horizon
75    // ignore the case where r <= atmosphere.inner_radius
76    let rho = sqrt(max(r * r - atm.inner_radius * atm.inner_radius, 0.0));
77
78    // Distance from a point at height r to the top atmosphere boundary at zenith angle mu
79    let d = distance_to_top_atmosphere_boundary(atm, r, mu);
80
81    // Minimum and maximum distance to the top atmosphere boundary from a point at height r
82    let d_min = atm.outer_radius - r; // length of the ray straight up to the top atmosphere boundary
83    let d_max = rho + H; // length of the ray to the top atmosphere boundary and grazing the horizon
84
85    let u = (d - d_min) / (d_max - d_min);
86    let v = rho / H;
87    return vec2<f32>(u, v);
88}
89
90/// Inverse of the mapping above, mapping from UV coordinates in the transmittance LUT to view height (r) and zenith cos angle (mu)
91fn transmittance_lut_uv_to_r_mu(uv: vec2<f32>) -> vec2<f32> {
92    // Distance to top atmosphere boundary for a horizontal ray at ground level
93    let H = sqrt(atmosphere.outer_radius * atmosphere.outer_radius - atmosphere.inner_radius * atmosphere.inner_radius);
94
95    // Distance to the horizon, from which we can compute r:
96    let rho = H * uv.y;
97    let r = sqrt(rho * rho + atmosphere.inner_radius * atmosphere.inner_radius);
98
99    // Distance to the top atmosphere boundary for the ray (r,mu), and its minimum
100    // and maximum values over all mu- obtained for (r,1) and (r,mu_horizon) -
101    // from which we can recover mu:
102    let d_min = atmosphere.outer_radius - r;
103    let d_max = rho + H;
104    let d = d_min + uv.x * (d_max - d_min);
105
106    var mu: f32;
107    if d == 0.0 {
108        mu = 1.0;
109    } else {
110        mu = (H * H - rho * rho - d * d) / (2.0 * r * d);
111    }
112
113    mu = clamp(mu, -1.0, 1.0);
114
115    return vec2<f32>(r, mu);
116}
117
118/// Simplified ray-sphere intersection
119/// where:
120/// Ray origin, o = [0,0,r] with r <= atmosphere.outer_radius
121/// mu is the cosine of spherical coordinate theta (-1.0 <= mu <= 1.0)
122/// so ray direction in spherical coordinates is [1,acos(mu),0] which needs to be converted to cartesian
123/// Direction of ray, u = [0,sqrt(1-mu*mu),mu]
124/// Center of sphere, c = [0,0,0]
125/// Radius of sphere, r = atmosphere.outer_radius
126/// This function solves the quadratic equation for line-sphere intersection simplified under these assumptions
127fn distance_to_top_atmosphere_boundary(atm: Atmosphere, r: f32, mu: f32) -> f32 {
128    // ignore the case where r > atm.outer_radius
129    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atm.outer_radius * atm.outer_radius, 0.0);
130    return max(-r * mu + sqrt(positive_discriminant), 0.0);
131}
132
133/// Simplified ray-sphere intersection
134/// as above for intersections with the ground
135fn distance_to_bottom_atmosphere_boundary(r: f32, mu: f32) -> f32 {
136    let positive_discriminant = max(r * r * (mu * mu - 1.0) + atmosphere.inner_radius * atmosphere.inner_radius, 0.0);
137    return max(-r * mu - sqrt(positive_discriminant), 0.0);
138}
139
140fn ray_intersects_ground(r: f32, mu: f32) -> bool {
141    return mu < 0.0 && r * r * (mu * mu - 1.0) + atmosphere.inner_radius * atmosphere.inner_radius >= 0.0;
142}
143
144// NOTICE: This function is copyrighted by Eric Bruneton and INRIA, and falls
145// under the license reproduced in bruneton_functions.wesl (variant of MIT license)
146fn sample_transmittance_lut_segment(r: f32, mu: f32, t: f32) -> vec3<f32> {
147    let r_t = get_local_r(r, mu, t);
148    let mu_t = clamp((r * mu + t) / r_t, -1.0, 1.0);
149
150    if ray_intersects_ground(r, mu) {
151        return min(
152            sample_transmittance_lut(r_t, -mu_t) / sample_transmittance_lut(r, -mu),
153            vec3(1.0)
154        );
155    } else {
156        return min(
157            sample_transmittance_lut(r, mu) / sample_transmittance_lut(r_t, mu_t), vec3(1.0)
158        );
159    }
160}