hw_skymodel/rgb/
mod.rs

1// An Analytic Model for Full Spectral Sky-Dome Radiance
2// Lukas Hosek & Alexander Wilkie
3// Project page: https://cgg.mff.cuni.cz/projects/SkylightModelling/
4// License file: hosek-wilkie-license.txt
5
6use std::f32::consts::{FRAC_PI_2, PI};
7
8use embed_doc_image::embed_doc_image;
9use thiserror::Error;
10
11/// R, G, or B channel.
12pub enum Channel {
13    R = 0,
14    G = 1,
15    B = 2,
16}
17
18/// Initial parameters for the sky model.
19#[derive(Clone, Copy, PartialEq, Debug)]
20pub struct SkyParams {
21    /// (Solar) elevation is given in radians. Elevation must be between `0..=π/2`.
22    pub elevation: f32,
23    /// Turbidity must be between `1..=10`.
24    pub turbidity: f32,
25    /// (Ground) albedo must be between `0..=1`.
26    pub albedo: [f32; 3],
27}
28
29impl Default for SkyParams {
30    fn default() -> Self {
31        Self {
32            elevation: 0.0,
33            turbidity: 1.0,
34            albedo: [1.0, 1.0, 1.0],
35        }
36    }
37}
38
39#[derive(Error, PartialEq, Debug)]
40pub enum Error {
41    #[error("Elevation must be between 0..=π/2, got {0} instead")]
42    ElevationOutOfRange(f32),
43    #[error("Turbidity must be between 1..=10, got {0} instead")]
44    TurbidityOutOfRange(f32),
45    #[error("Albedo must be between 0..=1, got {0:?} instead")]
46    AlbedoOutOfRange([f32; 3]),
47}
48
49/// The state of the sky model. Ideally, you should keep the state around as
50/// long as any of the `SkyParams` parameters don't change. If any parameters
51/// change, you must re-create the `SkyState` object.
52///
53/// Note: if you are planning to run the model in a shader, you only need to
54/// translate [`SkyState::radiance`] into shader code. The rest can be executed
55/// on the CPU. You grab the raw parameters using [`SkyState::raw`] and upload
56/// the 30 `f32`s to the GPU in a uniform buffer, for example.
57#[derive(Clone, Copy, PartialEq, Debug)]
58pub struct SkyState {
59    params: [f32; 27],
60    radiances: [f32; 3],
61}
62
63impl SkyState {
64    /// Creates `SkyState`.
65    ///
66    /// # Errors
67    /// Can fail if any of the parameters are out of range.
68    pub fn new(sky_params: &SkyParams) -> Result<Self, Error> {
69        // Validate parameters.
70        if !(0.0..=FRAC_PI_2).contains(&sky_params.elevation) {
71            return Err(Error::ElevationOutOfRange(sky_params.elevation));
72        }
73        if !(1.0..=10.0).contains(&sky_params.turbidity) {
74            return Err(Error::TurbidityOutOfRange(sky_params.turbidity));
75        }
76        for albedo in sky_params.albedo {
77            if !(0.0..=1.0).contains(&albedo) {
78                return Err(Error::AlbedoOutOfRange(sky_params.albedo));
79            }
80        }
81
82        // Load datasets.
83        macro_rules! dataset {
84            ($path: literal) => {{
85                // Note: Rust's include_bytes! doesn't guarantee that the byte
86                // slice is aligned. That's why we use include_bytes_aligned
87                // crate to ensure 4-byte alignment so we can cast into
88                // f32-slice.
89                use include_bytes_aligned::include_bytes_aligned;
90                let bytes = include_bytes_aligned!(4, $path);
91                bytemuck::cast_slice::<_, f32>(bytes)
92            }};
93        }
94        let params_r = dataset!("params-r");
95        let params_g = dataset!("params-g");
96        let params_b = dataset!("params-b");
97        let radiances_r = dataset!("radiances-r");
98        let radiances_g = dataset!("radiances-g");
99        let radiances_b = dataset!("radiances-b");
100
101        // Init state.
102        let mut params = [0.0; 27];
103        let mut radiances = [0.0; 3];
104        let elevation = sky_params.elevation;
105        let turbidity = sky_params.turbidity;
106        let albedo = sky_params.albedo;
107        let t = (elevation / (0.5 * PI)).powf(1.0 / 3.0);
108
109        init_params(&mut params[..], params_r, turbidity, albedo[0], t);
110        init_params(&mut params[9..], params_g, turbidity, albedo[1], t);
111        init_params(&mut params[(9 * 2)..], params_b, turbidity, albedo[2], t);
112        init_radiances(&mut radiances[0], radiances_r, turbidity, albedo[0], t);
113        init_radiances(&mut radiances[1], radiances_g, turbidity, albedo[1], t);
114        init_radiances(&mut radiances[2], radiances_b, turbidity, albedo[2], t);
115
116        Ok(Self { params, radiances })
117    }
118
119    /// Evaluates incoming radiance for a given channel. See the figure on how
120    /// `theta` and `gamma` are defined.
121    ///
122    /// ![Coordinate system of the sky model][coordinate-system]
123    ///
124    #[embed_doc_image("coordinate-system", "images/coordinate-system.png")]
125    pub fn radiance(&self, theta: f32, gamma: f32, channel: Channel) -> f32 {
126        let channel = channel as usize;
127        let r = self.radiances[channel];
128        let p = &self.params[(9 * channel)..];
129        let p0 = p[0];
130        let p1 = p[1];
131        let p2 = p[2];
132        let p3 = p[3];
133        let p4 = p[4];
134        let p5 = p[5];
135        let p6 = p[6];
136        let p7 = p[7];
137        let p8 = p[8];
138
139        let cos_gamma = gamma.cos();
140        let cos_gamma2 = cos_gamma * cos_gamma;
141        let cos_theta = theta.cos().abs();
142
143        let exp_m = (p4 * gamma).exp();
144        let ray_m = cos_gamma2;
145        let mie_m_lhs = 1.0 + cos_gamma2;
146        let mie_m_rhs = (1.0 + p8 * p8 - 2.0 * p8 * cos_gamma).powf(1.5);
147        let mie_m = mie_m_lhs / mie_m_rhs;
148        let zenith = cos_theta.sqrt();
149        let radiance_lhs = 1.0 + p0 * (p1 / (cos_theta + 0.01)).exp();
150        let radiance_rhs = p2 + p3 * exp_m + p5 * ray_m + p6 * mie_m + p7 * zenith;
151        let radiance_dist = radiance_lhs * radiance_rhs;
152        r * radiance_dist
153    }
154
155    /// Returns the internal state. Used when [`SkyState::radiance`] is
156    /// implemented and executed externally, for example as a GPU shader.
157    pub fn raw(&self) -> ([f32; 27], [f32; 3]) {
158        (self.params, self.radiances)
159    }
160}
161
162fn init_params(out_params: &mut [f32], dataset: &[f32], turbidity: f32, albedo: f32, t: f32) {
163    let turbidity_int = turbidity.trunc() as usize;
164    let turbidity_rem = turbidity.fract();
165    let turbidity_min = turbidity_int.saturating_sub(1);
166    let turbidity_max = turbidity_int.min(9);
167    let p0 = &dataset[(9 * 6 * turbidity_min)..];
168    let p1 = &dataset[(9 * 6 * turbidity_max)..];
169    let p2 = &dataset[(9 * 6 * 10 + 9 * 6 * turbidity_min)..];
170    let p3 = &dataset[(9 * 6 * 10 + 9 * 6 * turbidity_max)..];
171    let s0 = (1.0 - albedo) * (1.0 - turbidity_rem);
172    let s1 = (1.0 - albedo) * turbidity_rem;
173    let s2 = albedo * (1.0 - turbidity_rem);
174    let s3 = albedo * turbidity_rem;
175
176    for i in 0..9 {
177        out_params[i] += s0 * quintic::<9>(&p0[i..], t);
178        out_params[i] += s1 * quintic::<9>(&p1[i..], t);
179        out_params[i] += s2 * quintic::<9>(&p2[i..], t);
180        out_params[i] += s3 * quintic::<9>(&p3[i..], t);
181    }
182}
183
184fn init_radiances(out_radiance: &mut f32, dataset: &[f32], turbidity: f32, albedo: f32, t: f32) {
185    let turbidity_int = turbidity.trunc() as usize;
186    let turbidity_rem = turbidity.fract();
187    let turbidity_min = turbidity_int.saturating_sub(1);
188    let turbidity_max = turbidity_int.min(9);
189    let p0 = &dataset[(6 * turbidity_min)..];
190    let p1 = &dataset[(6 * turbidity_max)..];
191    let p2 = &dataset[(6 * 10 + 6 * turbidity_min)..];
192    let p3 = &dataset[(6 * 10 + 6 * turbidity_max)..];
193    let s0 = (1.0 - albedo) * (1.0 - turbidity_rem);
194    let s1 = (1.0 - albedo) * turbidity_rem;
195    let s2 = albedo * (1.0 - turbidity_rem);
196    let s3 = albedo * turbidity_rem;
197
198    *out_radiance += s0 * quintic::<1>(p0, t);
199    *out_radiance += s1 * quintic::<1>(p1, t);
200    *out_radiance += s2 * quintic::<1>(p2, t);
201    *out_radiance += s3 * quintic::<1>(p3, t);
202}
203
204fn quintic<const STRIDE: usize>(p: &[f32], t: f32) -> f32 {
205    let t2 = t * t;
206    let t3 = t2 * t;
207    let t4 = t2 * t2;
208    let t5 = t4 * t;
209
210    let inv_t = 1.0 - t;
211    let inv_t2 = inv_t * inv_t;
212    let inv_t3 = inv_t2 * inv_t;
213    let inv_t4 = inv_t2 * inv_t2;
214    let inv_t5 = inv_t4 * inv_t;
215
216    let m0 = p[0] * inv_t5;
217    let m1 = p[STRIDE] * 5.0 * inv_t4 * t;
218    let m2 = p[2 * STRIDE] * 10.0 * inv_t3 * t2;
219    let m3 = p[3 * STRIDE] * 10.0 * inv_t2 * t3;
220    let m4 = p[4 * STRIDE] * 5.0 * inv_t * t4;
221    let m5 = p[5 * STRIDE] * t5;
222
223    m0 + m1 + m2 + m3 + m4 + m5
224}
225
226#[cfg(test)]
227mod tests {
228    use super::*;
229
230    #[test]
231    fn default() {
232        let state = SkyState::new(&SkyParams::default());
233        assert!(state.is_ok());
234    }
235
236    #[test]
237    fn elevation_out_of_range() {
238        let params = SkyParams {
239            elevation: -1.0,
240            ..SkyParams::default()
241        };
242        let state = SkyState::new(&params);
243        assert_eq!(state, Err(Error::ElevationOutOfRange(-1.0)));
244        assert_eq!(
245            "Elevation must be between 0..=π/2, got -1 instead",
246            format!("{}", state.unwrap_err())
247        );
248    }
249
250    #[test]
251    fn turbidity_out_of_range() {
252        let params = SkyParams {
253            turbidity: 0.0,
254            ..SkyParams::default()
255        };
256        let state = SkyState::new(&params);
257        assert_eq!(state, Err(Error::TurbidityOutOfRange(0.0)));
258        assert_eq!(
259            "Turbidity must be between 1..=10, got 0 instead",
260            format!("{}", state.unwrap_err())
261        );
262    }
263
264    #[test]
265    fn albedo_out_of_range() {
266        let params = SkyParams {
267            albedo: [0.0, 2.0, 0.0],
268            ..SkyParams::default()
269        };
270        let state = SkyState::new(&params);
271        assert_eq!(state, Err(Error::AlbedoOutOfRange([0.0, 2.0, 0.0])));
272        assert_eq!(
273            "Albedo must be between 0..=1, got [0.0, 2.0, 0.0] instead",
274            format!("{}", state.unwrap_err())
275        );
276    }
277}
hw_skymodel/rgb/mod.rs

hw_skymodel/rgb/
mod.rs
