1use bytemuck::{Pod, Zeroable};
2#[cfg(test)]
3use pfx_core::daylight::Daylight;
4use pfx_core::daylight::{self, Probe};
5pub use pfx_core::sky::AnalyticSky;
6use pfx_core::sky::{SKY_MODEL_WGSL, SUN_RADIUS};
7use pfx_load::Sky as HdrSky;
8use std::f32::consts::{PI, TAU};
9use std::sync::LazyLock;
10
11#[derive(Clone, Debug)]
12pub enum Environment {
13 Analytic(AnalyticSky),
14 Hdr(HdrSky),
15}
16
17impl From<HdrSky> for Environment {
18 fn from(sky: HdrSky) -> Self {
19 Self::Hdr(sky)
20 }
21}
22
23#[repr(C)]
24#[derive(Clone, Copy, Debug, Pod, Zeroable)]
25pub struct EnvironmentUniform {
26 pub sun: [f32; 4],
27 pub sun_colour_intensity: [f32; 4],
28 pub ambient_turbidity_albedo: [f32; 4],
29 pub albedo_mode: [f32; 4],
30 pub forward: [f32; 4],
31 pub right: [f32; 4],
32 pub up: [f32; 4],
33}
34
35impl Environment {
36 pub fn uniform(&self) -> EnvironmentUniform {
37 let (model, mode) = match self {
38 Self::Analytic(model) => (*model, 0.0),
39 Self::Hdr(_) => (
40 AnalyticSky {
41 sun: [0.0, 1.0, 0.0],
42 sun_colour: [1.0; 3],
43 sun_intensity: 0.0,
44 ambient: 0.0,
45 turbidity: 2.0,
46 ground_albedo: [0.0; 3],
47 },
48 1.0,
49 ),
50 };
51 EnvironmentUniform {
52 sun: [model.sun[0], model.sun[1], model.sun[2], 0.0],
53 sun_colour_intensity: [
54 model.sun_colour[0],
55 model.sun_colour[1],
56 model.sun_colour[2],
57 model.sun_intensity,
58 ],
59 ambient_turbidity_albedo: [
60 model.ambient,
61 model.turbidity,
62 model.ground_albedo[0],
63 model.ground_albedo[1],
64 ],
65 albedo_mode: [model.ground_albedo[2], mode, 0.0, 0.0],
66 forward: [0.0; 4],
67 right: [0.0; 4],
68 up: [0.0; 4],
69 }
70 }
71
72 pub fn from_probes(hour: f64, skies: &[HdrSky], probes: &[Probe]) -> Result<Self, String> {
73 let first = skies.first().ok_or("sky probes need at least one HDR")?;
74 if skies.iter().any(|s| {
75 s.width != first.width
76 || s.height != first.height
77 || s.texels.len() != first.texels.len()
78 }) {
79 return Err("sky probes need matching HDR dimensions".into());
80 }
81 let weights = daylight::probe_weights(hour, probes, skies.len());
82 if weights.iter().sum::<f64>() <= 0.0 {
83 return Err("sky probes have no valid index".into());
84 }
85 let mut texels = vec![[0.0; 4]; first.texels.len()];
86 for (sky, weight) in skies.iter().zip(weights) {
87 for (out, input) in texels.iter_mut().zip(&sky.texels) {
88 for (out_channel, input_channel) in out[..3].iter_mut().zip(&input[..3]) {
89 *out_channel += *input_channel * weight as f32;
90 }
91 }
92 }
93 for pixel in &mut texels {
94 pixel[3] = 1.0;
95 }
96 Ok(Self::Hdr(HdrSky {
97 width: first.width,
98 height: first.height,
99 texels,
100 }))
101 }
102
103 pub fn radiance(&self, direction: [f32; 3]) -> [f32; 3] {
104 match self {
105 Self::Analytic(model) => model.radiance(direction),
106 Self::Hdr(hdr) => {
107 let d = unit(direction);
108 let u = (d[0].atan2(-d[2]) / TAU + 0.5).rem_euclid(1.0);
109 let v = d[1].clamp(-1.0, 1.0).acos() / PI;
110 let x = u * hdr.width as f32 - 0.5;
111 let y = v * hdr.height as f32 - 0.5;
112 let x0 = x.floor();
113 let y0 = y.floor();
114 let tx = x - x0;
115 let ty = y - y0;
116 let at = |ix: f32, iy: f32| {
117 let col = (ix as i32).rem_euclid(hdr.width as i32) as usize;
118 let row = (iy as i32).clamp(0, hdr.height as i32 - 1) as usize;
119 hdr.texels[row * hdr.width as usize + col]
120 };
121 let a = at(x0, y0);
122 let b = at(x0 + 1.0, y0);
123 let c = at(x0, y0 + 1.0);
124 let d = at(x0 + 1.0, y0 + 1.0);
125 [0, 1, 2].map(|i| {
126 (a[i] * (1.0 - tx) + b[i] * tx) * (1.0 - ty)
127 + (c[i] * (1.0 - tx) + d[i] * tx) * ty
128 })
129 }
130 }
131 }
132
133 fn diffuse(&self, direction: [f32; 3]) -> [f32; 3] {
134 match self {
135 Self::Analytic(model) => model.diffuse(direction),
136 Self::Hdr(_) => self.radiance(direction),
137 }
138 }
139}
140
141#[derive(Clone, Debug)]
142pub struct SkyCdf {
143 pub width: u32,
144 pub height: u32,
145 pub rows: Vec<f32>,
146 pub columns: Vec<f32>,
147 pub sun_probability: f32,
148 pub total: f32,
149 sun: Option<[f32; 3]>,
150}
151
152#[repr(C)]
153#[derive(Clone, Copy, Debug, Pod, Zeroable)]
154pub struct SkyCdfUniform {
155 pub width: u32,
156 pub height: u32,
157 pub sun_probability: f32,
158 pub unused: f32,
159}
160
161#[derive(Clone, Copy, Debug)]
162pub struct SkySample {
163 pub direction: [f32; 3],
164 pub radiance: [f32; 3],
165 pub pdf: f32,
166}
167
168impl SkyCdf {
169 pub fn without_sun(mut self) -> Self {
170 self.sun_probability = 0.0;
171 self.sun = None;
172 self
173 }
174
175 pub fn uniform(&self) -> SkyCdfUniform {
176 SkyCdfUniform {
177 width: self.width,
178 height: self.height,
179 sun_probability: self.sun_probability,
180 unused: 0.0,
181 }
182 }
183
184 pub fn build(environment: &Environment, width: u32, height: u32) -> Self {
185 let width = width.max(1);
186 let height = height.max(1);
187 let mut rows = vec![0.0; height as usize];
188 let mut columns = vec![0.0; (width * height) as usize];
189 let mut total = 0.0_f64;
190 for row in 0..height {
191 let theta = PI * (row as f32 + 0.5) / height as f32;
192 let sine = theta.sin();
193 let mut row_total = 0.0_f64;
194 for col in 0..width {
195 let phi = TAU * (col as f32 + 0.5) / width as f32;
196 let direction = [
197 theta.sin() * phi.cos(),
198 theta.cos(),
199 theta.sin() * phi.sin(),
200 ];
201 let color = environment.diffuse(direction);
202 row_total += luminance(color).max(0.0) as f64 * sine as f64;
203 columns[(row * width + col) as usize] = row_total as f32;
204 }
205 if row_total > 0.0 {
206 for col in 0..width {
207 columns[(row * width + col) as usize] /= row_total as f32;
208 }
209 } else {
210 for col in 0..width {
211 columns[(row * width + col) as usize] = (col + 1) as f32 / width as f32;
212 }
213 }
214 total += row_total;
215 rows[row as usize] = total as f32;
216 }
217 if total > 0.0 {
218 for value in &mut rows {
219 *value /= total as f32;
220 }
221 } else {
222 for (index, value) in rows.iter_mut().enumerate() {
223 *value = (index + 1) as f32 / height as f32;
224 }
225 }
226 let (sun_probability, sun) = match environment {
227 Environment::Analytic(model) if model.sun[1] > 0.0 && model.sun_intensity > 0.0 => {
228 let sky_energy = total as f32 * (TAU / width as f32) * (PI / height as f32);
229 let sun_energy = luminance(model.sun_colour) * model.sun_intensity;
230 (
231 sun_energy / (sun_energy + sky_energy).max(1e-6),
232 Some(model.sun),
233 )
234 }
235 _ => (0.0, None),
236 };
237 Self {
238 width,
239 height,
240 rows,
241 columns,
242 sun_probability,
243 total: total as f32,
244 sun,
245 }
246 }
247
248 pub fn sample(&self, environment: &Environment, u: f32, v: f32) -> SkySample {
249 let u = u.clamp(0.0, 1.0 - f32::EPSILON);
250 let v = v.clamp(0.0, 1.0 - f32::EPSILON);
251 let mut cell_pdf = None;
252 let direction = if let Some(sun) = self.sun.filter(|_| u < self.sun_probability) {
253 let cosine = 1.0 - (u / self.sun_probability.max(1e-6)) * (1.0 - SUN_RADIUS.cos());
254 let radial = cosine.acos();
255 let angle = TAU * v;
256 let up = if sun[1].abs() < 0.99 {
257 [0.0, 1.0, 0.0]
258 } else {
259 [1.0, 0.0, 0.0]
260 };
261 let right = unit(cross(up, sun));
262 let up = cross(sun, right);
263 unit(add(
264 add(
265 mul(sun, radial.cos()),
266 mul(right, radial.sin() * angle.cos()),
267 ),
268 mul(up, radial.sin() * angle.sin()),
269 ))
270 } else {
271 let row = self
272 .rows
273 .partition_point(|&value| value <= v)
274 .min(self.height as usize - 1);
275 let row_start = if row == 0 { 0.0 } else { self.rows[row - 1] };
276 let row_end = self.rows[row];
277 let fy = (v - row_start) / (row_end - row_start).max(1e-6);
278 let target =
279 (u - self.sun_probability).max(0.0) / (1.0 - self.sun_probability).max(1e-6);
280 let start = row * self.width as usize;
281 let cdf = &self.columns[start..start + self.width as usize];
282 let col = cdf
283 .partition_point(|&value| value <= target)
284 .min(self.width as usize - 1);
285 let col_start = if col == 0 { 0.0 } else { cdf[col - 1] };
286 let col_end = cdf[col];
287 let fx = (target - col_start) / (col_end - col_start).max(1e-6);
288 let theta = PI * (row as f32 + fy) / self.height as f32;
289 let phi = TAU * (col as f32 + fx) / self.width as f32;
290 let solid_angle =
291 TAU / self.width as f32 * PI / self.height as f32 * theta.sin().max(1e-5);
292 let sun_pdf = match self.sun {
293 Some(sun) if dot(sun, unit_direction(theta, phi)) >= SUN_RADIUS.cos() => {
294 self.sun_probability / (2.0 * PI * (1.0 - SUN_RADIUS.cos()))
295 }
296 _ => 0.0,
297 };
298 cell_pdf = Some(
299 (1.0 - self.sun_probability) * (row_end - row_start) * (col_end - col_start)
300 / solid_angle
301 + sun_pdf,
302 );
303 [
304 theta.sin() * phi.cos(),
305 theta.cos(),
306 theta.sin() * phi.sin(),
307 ]
308 };
309 SkySample {
310 direction,
311 radiance: environment.radiance(direction),
312 pdf: cell_pdf.unwrap_or_else(|| self.pdf(direction)),
313 }
314 }
315
316 pub fn pdf(&self, direction: [f32; 3]) -> f32 {
317 let d = unit(direction);
318 let theta = d[1].clamp(-1.0, 1.0).acos();
319 let phi = d[2].atan2(d[0]).rem_euclid(TAU);
320 let row = ((theta / PI * self.height as f32) as u32).min(self.height - 1) as usize;
321 let col = ((phi / TAU * self.width as f32) as u32).min(self.width - 1) as usize;
322 let start = row * self.width as usize;
323 let row_mass = self.rows[row] - if row == 0 { 0.0 } else { self.rows[row - 1] };
324 let col_mass = self.columns[start + col]
325 - if col == 0 {
326 0.0
327 } else {
328 self.columns[start + col - 1]
329 };
330 let solid_angle = TAU / self.width as f32 * PI / self.height as f32 * theta.sin().max(1e-5);
331 let sky_pdf = (1.0 - self.sun_probability) * row_mass * col_mass / solid_angle;
332 let sun_pdf = if self.sun.is_some_and(|sun| dot(sun, d) >= SUN_RADIUS.cos()) {
333 self.sun_probability / (2.0 * PI * (1.0 - SUN_RADIUS.cos()))
334 } else {
335 0.0
336 };
337 sky_pdf + sun_pdf
338 }
339}
340
341fn unit_direction(theta: f32, phi: f32) -> [f32; 3] {
342 [
343 theta.sin() * phi.cos(),
344 theta.cos(),
345 theta.sin() * phi.sin(),
346 ]
347}
348fn luminance(v: [f32; 3]) -> f32 {
349 0.2126 * v[0] + 0.7152 * v[1] + 0.0722 * v[2]
350}
351fn dot(a: [f32; 3], b: [f32; 3]) -> f32 {
352 a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
353}
354fn add(a: [f32; 3], b: [f32; 3]) -> [f32; 3] {
355 [a[0] + b[0], a[1] + b[1], a[2] + b[2]]
356}
357fn mul(a: [f32; 3], s: f32) -> [f32; 3] {
358 a.map(|v| v * s)
359}
360fn cross(a: [f32; 3], b: [f32; 3]) -> [f32; 3] {
361 [
362 a[1] * b[2] - a[2] * b[1],
363 a[2] * b[0] - a[0] * b[2],
364 a[0] * b[1] - a[1] * b[0],
365 ]
366}
367fn unit(a: [f32; 3]) -> [f32; 3] {
368 mul(a, 1.0 / dot(a, a).sqrt().max(1e-8))
369}
370
371pub static ENVIRONMENT_WGSL: LazyLock<String> =
372 LazyLock::new(|| format!("{}\n{}", SKY_MODEL_WGSL, SKY_SAMPLING_WGSL));
373
374const SKY_SAMPLING_WGSL: &str = r#"
375fn sky_surface_radiance(direction: vec3f) -> vec3f {
376 if sky.albedo_mode.y < 0.5 { return sky_diffuse(direction); }
377 return sky_radiance(direction);
378}
379struct SkyCdfInfo {
380 width: u32,
381 height: u32,
382 sun_probability: f32,
383 unused: f32,
384};
385@group(1) @binding(3) var<storage, read> sky_rows: array<f32>;
386@group(1) @binding(4) var<storage, read> sky_columns: array<f32>;
387@group(1) @binding(5) var<uniform> sky_cdf: SkyCdfInfo;
388struct SkyDirectionSample {
389 direction: vec3f,
390 radiance: vec3f,
391 pdf: f32,
392};
393fn sky_find_row(needle: f32) -> u32 {
394 var low = 0u;
395 var high = sky_cdf.height;
396 loop {
397 if low >= high { break; }
398 let mid = (low + high) / 2u;
399 if sky_rows[mid] <= needle { low = mid + 1u; } else { high = mid; }
400 }
401 return min(low, sky_cdf.height - 1u);
402}
403fn sky_find_column(row: u32, needle: f32) -> u32 {
404 var low = 0u;
405 var high = sky_cdf.width;
406 loop {
407 if low >= high { break; }
408 let mid = (low + high) / 2u;
409 if sky_columns[row * sky_cdf.width + mid] <= needle { low = mid + 1u; } else { high = mid; }
410 }
411 return min(low, sky_cdf.width - 1u);
412}
413fn sky_pdf(direction: vec3f) -> f32 {
414 let d = normalize(direction);
415 let theta = acos(clamp(d.y, -1.0, 1.0));
416 let phi = (atan2(d.z, d.x) + 6.28318531) % 6.28318531;
417 let row = min(u32(theta / 3.14159265 * f32(sky_cdf.height)), sky_cdf.height - 1u);
418 let col = min(u32(phi / 6.28318531 * f32(sky_cdf.width)), sky_cdf.width - 1u);
419 let row_before = select(0.0, sky_rows[max(row, 1u) - 1u], row > 0u);
420 let col_before = select(0.0, sky_columns[row * sky_cdf.width + max(col, 1u) - 1u], col > 0u);
421 let mass = (sky_rows[row] - row_before) * (sky_columns[row * sky_cdf.width + col] - col_before);
422 let area = 6.28318531 / f32(sky_cdf.width) * 3.14159265 / f32(sky_cdf.height) * max(sin(theta), 0.00001);
423 let sun_pdf = select(0.0, sky_cdf.sun_probability / (6.28318531 * (1.0 - cos(0.004675))), sky.sun.y > 0.0 && dot(d, sky.sun.xyz) >= cos(0.004675));
424 return (1.0 - sky_cdf.sun_probability) * mass / area + sun_pdf;
425}
426fn sky_sample(random: vec2f) -> SkyDirectionSample {
427 let u = clamp(random.x, 0.0, 0.99999994);
428 let v = clamp(random.y, 0.0, 0.99999994);
429 var direction: vec3f;
430 if sky_cdf.sun_probability > 0.0 && u < sky_cdf.sun_probability {
431 let cosine = 1.0 - u / sky_cdf.sun_probability * (1.0 - cos(0.004675));
432 let radial = acos(cosine);
433 let angle = 6.28318531 * v;
434 let axis = select(vec3f(1.0, 0.0, 0.0), vec3f(0.0, 1.0, 0.0), abs(sky.sun.y) < 0.99);
435 let right = normalize(cross(axis, sky.sun.xyz));
436 let up = cross(sky.sun.xyz, right);
437 direction = normalize(sky.sun.xyz * cos(radial) + (right * cos(angle) + up * sin(angle)) * sin(radial));
438 } else {
439 let row = sky_find_row(v);
440 let row_before = select(0.0, sky_rows[max(row, 1u) - 1u], row > 0u);
441 let fy = (v - row_before) / max(sky_rows[row] - row_before, 0.000001);
442 let needle = max(u - sky_cdf.sun_probability, 0.0) / max(1.0 - sky_cdf.sun_probability, 0.000001);
443 let col = sky_find_column(row, needle);
444 let col_before = select(0.0, sky_columns[row * sky_cdf.width + max(col, 1u) - 1u], col > 0u);
445 let fx = (needle - col_before) / max(sky_columns[row * sky_cdf.width + col] - col_before, 0.000001);
446 let theta = 3.14159265 * (f32(row) + fy) / f32(sky_cdf.height);
447 let phi = 6.28318531 * (f32(col) + fx) / f32(sky_cdf.width);
448 direction = vec3f(sin(theta) * cos(phi), cos(theta), sin(theta) * sin(phi));
449 let area = 6.28318531 / f32(sky_cdf.width) * 3.14159265 / f32(sky_cdf.height) * max(sin(theta), 0.00001);
450 let mass = (sky_rows[row] - row_before) * (sky_columns[row * sky_cdf.width + col] - col_before);
451 let sun_pdf = select(0.0, sky_cdf.sun_probability / (6.28318531 * (1.0 - cos(0.004675))), sky.sun.y > 0.0 && dot(direction, sky.sun.xyz) >= cos(0.004675));
452 return SkyDirectionSample(direction, sky_radiance(direction), (1.0 - sky_cdf.sun_probability) * mass / area + sun_pdf);
453 }
454 return SkyDirectionSample(direction, sky_radiance(direction), sky_pdf(direction));
455}
456"#;
457
458#[cfg(test)]
459mod tests {
460 use super::*;
461
462 fn model(hour: f64) -> AnalyticSky {
463 AnalyticSky::new(
464 Daylight {
465 hour,
466 day: 172.0,
467 latitude: 45.0,
468 heading: 180.0,
469 },
470 daylight::REFERENCE_HOUR,
471 2.5,
472 [0.2, 0.25, 0.2],
473 )
474 }
475
476 #[test]
477 fn cdf_sums_to_one_and_samples_sun() {
478 let environment = Environment::Analytic(model(12.0));
479 let cdf = SkyCdf::build(&environment, 64, 32);
480 let mut mass = 0.0;
481 for row in 0..cdf.height as usize {
482 let row_mass = cdf.rows[row] - if row == 0 { 0.0 } else { cdf.rows[row - 1] };
483 for col in 0..cdf.width as usize {
484 let start = row * cdf.width as usize;
485 let col_mass = cdf.columns[start + col]
486 - if col == 0 {
487 0.0
488 } else {
489 cdf.columns[start + col - 1]
490 };
491 assert!(col_mass >= -1e-5 && row_mass >= -1e-5);
492 mass += row_mass * col_mass;
493 }
494 }
495 assert!((mass - 1.0).abs() < 1e-4);
496 assert!((cdf.rows.last().unwrap() - 1.0).abs() < 1e-6);
497 assert!(cdf.sun_probability > 0.0);
498 let solar = cdf.sample(&environment, cdf.sun_probability * 0.5, 0.25);
499 assert!(dot(solar.direction, model(12.0).sun) > SUN_RADIUS.cos());
500 assert!(solar.pdf.is_finite() && solar.pdf > 0.0);
501 }
502
503 #[test]
504 fn night_is_darker_and_hdr_probe_blends() {
505 let noon = model(12.0);
506 let night = model(0.0);
507 assert!(noon.radiance([0.0, 1.0, 0.0])[2] > night.radiance([0.0, 1.0, 0.0])[2] * 5.0);
508 let skies = [
509 HdrSky {
510 width: 1,
511 height: 1,
512 texels: vec![[2.0, 0.0, 0.0, 1.0]],
513 },
514 HdrSky {
515 width: 1,
516 height: 1,
517 texels: vec![[0.0, 0.0, 4.0, 1.0]],
518 },
519 ];
520 let probes = [
521 Probe {
522 hour: 6.0,
523 index: 0,
524 },
525 Probe {
526 hour: 18.0,
527 index: 1,
528 },
529 ];
530 let environment = Environment::from_probes(12.0, &skies, &probes).unwrap();
531 assert_eq!(environment.radiance([1.0, 0.0, 0.0]), [1.0, 0.0, 2.0]);
532 let cdf = SkyCdf::build(&environment, 16, 8);
533 assert_eq!(cdf.sun_probability, 0.0);
534 assert!((cdf.rows.last().unwrap() - 1.0).abs() < 1e-6);
535 }
536
537 #[test]
538 fn wgsl_parses() {
539 let module = naga::front::wgsl::parse_str(&ENVIRONMENT_WGSL).unwrap();
540 naga::valid::Validator::new(
541 naga::valid::ValidationFlags::all(),
542 naga::valid::Capabilities::all(),
543 )
544 .validate(&module)
545 .unwrap();
546 }
547
548 #[test]
549 fn sky_sampling_converges_to_analytic_integral() {
550 let model = model(16.0);
551 let environment = Environment::Analytic(model);
552 let full_cdf = SkyCdf::build(&environment, 256, 128);
553 let cdf = full_cdf.clone().without_sun();
554 let mut reference = [0.0_f64; 3];
555 for row in 0..256 {
556 let theta = PI * (row as f32 + 0.5) / 256.0;
557 for col in 0..512 {
558 let phi = TAU * (col as f32 + 0.5) / 512.0;
559 let direction = [
560 theta.sin() * phi.cos(),
561 theta.cos(),
562 theta.sin() * phi.sin(),
563 ];
564 let value = model.diffuse(direction);
565 for channel in 0..3 {
566 reference[channel] += (value[channel] * theta.sin()) as f64;
567 }
568 }
569 }
570 let area = (PI * TAU / (256 * 512) as f32) as f64;
571 let reference = reference.map(|v| v * area);
572 let mut estimate = [0.0_f64; 3];
573 let count = 65536_u32;
574 for i in 0..count {
575 let u = (i as f32 + 0.5) / count as f32;
576 let v = (i.reverse_bits() as f32 + 0.5) * 2.328_306_4e-10;
577 let sample = cdf.sample(&environment, u, v);
578 let value = model.diffuse(sample.direction);
579 for channel in 0..3 {
580 estimate[channel] += (value[channel] / sample.pdf) as f64;
581 }
582 }
583 for channel in 0..3 {
584 let measured = estimate[channel] / count as f64;
585 assert!(
586 (measured / reference[channel] - 1.0).abs() < 0.02,
587 "{measured} vs {}",
588 reference[channel]
589 );
590 }
591 let mut full_estimate = [0.0_f64; 3];
592 for i in 0..count {
593 let u = (i as f32 + 0.5) / count as f32;
594 let v = (i.reverse_bits() as f32 + 0.5) * 2.328_306_4e-10;
595 let sample = full_cdf.sample(&environment, u, v);
596 for (sum, radiance) in full_estimate.iter_mut().zip(sample.radiance) {
597 *sum += (radiance / sample.pdf) as f64;
598 }
599 }
600 for channel in 0..3 {
601 let measured = full_estimate[channel] / count as f64;
602 let expected =
603 reference[channel] + (model.sun_intensity * model.sun_colour[channel]) as f64;
604 assert!(
605 (measured / expected - 1.0).abs() < 0.02,
606 "{measured} vs {expected}"
607 );
608 }
609 }
610}
611
612#[cfg(test)]
613mod density_tests {
614 use super::*;
615
616 fn speckled(width: u32, height: u32) -> HdrSky {
617 let mut state = 0x2545f491_u32;
618 let mut next = move || {
619 state ^= state << 13;
620 state ^= state >> 17;
621 state ^= state << 5;
622 state as f32 / u32::MAX as f32
623 };
624 let texels = (0..width * height)
625 .map(|_| {
626 let value = if next() < 0.11 { 0.0 } else { 20.0 * next() };
627 [value, value * 0.9, value * 0.8, 1.0]
628 })
629 .collect();
630 HdrSky {
631 width,
632 height,
633 texels,
634 }
635 }
636
637 #[test]
638 fn the_pdf_is_the_density_of_what_the_sampler_draws() {
639 let (width, height) = (128_usize, 64_usize);
640 let environment = Environment::Hdr(speckled(width as u32, height as u32));
641 let cdf = SkyCdf::build(&environment, width as u32, height as u32);
642 let mut state = 0x1234567_u32;
643 let mut next = move || {
644 state ^= state << 13;
645 state ^= state >> 17;
646 state ^= state << 5;
647 (state >> 8) as f32 / (1u32 << 24) as f32
648 };
649 let draws = 2_000_000;
650 let mut counts = vec![0u32; width * height];
651 let mut no_pdf = 0;
652 let mut wrong_cell = 0;
653 for _ in 0..draws {
654 let sample = cdf.sample(&environment, next(), next());
655 let d = sample.direction;
656 let theta = d[1].clamp(-1.0, 1.0).acos();
657 let phi = d[2].atan2(d[0]).rem_euclid(TAU);
658 let row = ((theta / PI * height as f32) as usize).min(height - 1);
659 let col = ((phi / TAU * width as f32) as usize).min(width - 1);
660 counts[row * width + col] += 1;
661 if sample.pdf <= 0.0 {
662 no_pdf += 1;
663 continue;
664 }
665 let area = (TAU / width as f32 * PI / height as f32 * theta.sin()) as f64;
666 let row_mass = cdf.rows[row] - if row == 0 { 0.0 } else { cdf.rows[row - 1] };
667 let at = row * width + col;
668 let col_mass = cdf.columns[at] - if col == 0 { 0.0 } else { cdf.columns[at - 1] };
669 let density = row_mass as f64 * col_mass as f64 / area;
670 if (cdf.pdf(d) as f64 / density - 1.0).abs() > 1e-3 {
671 wrong_cell += 1;
672 }
673 assert!(sample.pdf as f64 > 0.0);
674 }
675 println!(
676 "no pdf {no_pdf}, direction lookup in the neighbouring cell {wrong_cell}, of {draws}"
677 );
678 assert_eq!(no_pdf, 0);
679 assert!(wrong_cell < draws / 1_000);
680 for row in 0..height {
681 let row_mass = cdf.rows[row] - if row == 0 { 0.0 } else { cdf.rows[row - 1] };
682 for col in 0..width {
683 let at = row * width + col;
684 let col_mass = cdf.columns[at] - if col == 0 { 0.0 } else { cdf.columns[at - 1] };
685 let expected = row_mass as f64 * col_mass as f64 * draws as f64;
686 if expected < 200.0 {
687 continue;
688 }
689 let got = counts[at] as f64;
690 assert!(
691 (got - expected).abs() < 5.0 * expected.sqrt(),
692 "cell {row},{col}: {got} drawn, {expected} expected"
693 );
694 }
695 }
696 }
697
698 #[test]
699 fn a_sun_adds_its_density_to_every_sample() {
700 let daylight = Daylight {
701 hour: 16.0,
702 day: 172.0,
703 latitude: 45.0,
704 heading: 180.0,
705 };
706 let environment = Environment::Analytic(AnalyticSky::new(
707 daylight,
708 daylight::REFERENCE_HOUR,
709 2.5,
710 [0.2, 0.25, 0.2],
711 ));
712 let cdf = SkyCdf::build(&environment, 128, 64);
713 assert!(cdf.sun_probability > 0.0);
714 let mut state = 0x7654321_u32;
715 let mut next = move || {
716 state ^= state << 13;
717 state ^= state >> 17;
718 state ^= state << 5;
719 (state >> 8) as f32 / (1u32 << 24) as f32
720 };
721 let draws = 300_000;
722 let mut apart = 0;
723 let mut sun_draws = 0;
724 for _ in 0..draws {
725 let sample = cdf.sample(&environment, next(), next());
726 assert!(sample.pdf > 0.0);
727 let looked_up = cdf.pdf(sample.direction);
728 if (sample.pdf / looked_up - 1.0).abs() > 1e-3 {
729 apart += 1;
730 }
731 if dot(sample.direction, cdf.sun.unwrap()) >= SUN_RADIUS.cos() {
732 sun_draws += 1;
733 }
734 }
735 assert!(sun_draws > draws / 100);
736 assert!(apart < draws / 1_000, "{apart}");
737 }
738}