1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
use crate::math::{Dir3, Vec3};
use lazy_static::lazy_static;
use rand::{rngs::ThreadRng, Rng};
use std::f64::consts::{FRAC_PI_2, PI};
lazy_static! {
static ref GOLDEN_RATIO: f64 = (1.0 + 5.0_f64.sqrt()) / 2.0;
}
#[inline]
#[must_use]
pub fn sample_henyey_greenstein(rng: &mut ThreadRng, asym: f64) -> f64 {
debug_assert!(asym.abs() <= 1.0);
if asym.abs() < 1.0e-6 {
return rng.gen_range(-1.0_f64, 1.0).acos();
}
((1.0 + asym.powi(2)
- ((1.0 - asym.powi(2)) / asym.mul_add(rng.gen_range(-1.0, 1.0), 1.0)).powi(2))
/ (2.0 * asym))
.acos()
}
#[inline]
#[must_use]
pub fn sample_normal(rng: &mut ThreadRng) -> f64 {
let a = (-2.0 * rng.gen_range(0.0_f64, 1.0).ln()).sqrt();
let theta = rng.gen_range(0.0, 2.0 * PI);
a * theta.cos()
}
#[inline]
#[must_use]
pub fn sample_gaussian(rng: &mut ThreadRng, mu: f64, sigma: f64) -> f64 {
debug_assert!(sigma > 0.0);
sample_normal(rng).mul_add(sigma, mu)
}
#[inline]
#[must_use]
pub fn rand_isotropic_dir(rng: &mut ThreadRng) -> Dir3 {
let theta = rng.gen_range(0.0, 2.0 * PI);
let z: f64 = rng.gen_range(-1.0, 1.0);
let v = (1.0 - z.powi(2)).sqrt();
let x = v * theta.cos();
let y = v * theta.sin();
Dir3::new_normalize(Vec3::new(x, y, z))
}
#[inline]
#[must_use]
pub fn rand_circle_point(n: i32, max: i32) -> (f64, f64) {
debug_assert!(n >= 0);
debug_assert!(n < max);
let r = f64::from(n) / f64::from(max - 1);
let theta = f64::from(n) * *GOLDEN_RATIO;
(r, theta)
}
#[inline]
#[must_use]
pub fn rand_sphere_point(n: i32, max: i32) -> (f64, f64) {
debug_assert!(n >= 0);
debug_assert!(n < max);
let d = f64::from(1 - max).mul_add(0.5, f64::from(n));
let phi = ((2.0 * d) / f64::from(max)).asin() + FRAC_PI_2;
let theta = ((2.0 * PI) / *GOLDEN_RATIO) * (d % *GOLDEN_RATIO);
(phi, theta)
}
#[inline]
#[must_use]
pub fn rand_hemisphere_point(n: i32, max: i32) -> (f64, f64) {
debug_assert!(n >= 0);
debug_assert!(n < max);
rand_sphere_point(n, max * 2)
}