pbrt_r3/core/base/
functions.rs

1use crate::core::base::Float;
2
3pub const MACHINE_EPSILON: Float = Float::EPSILON * 0.5;
4
5#[inline]
6pub fn gamma(n: Float) -> Float {
7    //return (n * MachineEpsilon) / (1 - n * MachineEpsilon);
8    let e = n * MACHINE_EPSILON;
9    return e / (1.0 - e);
10}
11
12#[inline]
13pub fn gamma_correct(value: Float) -> Float {
14    if value <= 0.0031308 {
15        return 12.92 * value;
16    } else {
17        return 1.055 * Float::powf(value, 1.0 / 2.4) - 0.055;
18    }
19}
20
21#[inline]
22pub fn inverse_gamma_correct(value: Float) -> Float {
23    if value <= 0.04045 {
24        return value * 1.0 / 12.92;
25    } else {
26        return Float::powf((value + 0.055) * 1.0 / 1.055, 2.4);
27    }
28}
29
30#[inline]
31pub fn lerp(t: Float, v1: Float, v2: Float) -> Float {
32    return (1.0 - t) * v1 + t * v2;
33}
34
35#[inline]
36pub fn quadratic(a: Float, b: Float, c: Float) -> Option<(Float, Float)> {
37    let a = a as f64;
38    let b = b as f64;
39    let c = c as f64;
40    // Find quadratic diqscriminant
41    let discrim: f64 = b * b - 4.0 * a * c;
42    if discrim < 0.0 {
43        return None;
44    }
45    let root_discrim = f64::sqrt(discrim);
46    // Compute quadratic _t_ values
47    let q = if b < 0.0 {
48        -0.5 * (b - root_discrim)
49    } else {
50        -0.5 * (b + root_discrim)
51    };
52    let t0 = q / a;
53    let t1 = c / q;
54    if t0 > t1 {
55        std::mem::swap(&mut &t0, &mut &t1);
56    }
57
58    let t0 = t0 as Float;
59    let t1 = t1 as Float;
60    return Some((t0, t1));
61}
62
63#[inline]
64pub fn is_power_of_2(v: u32) -> bool {
65    return (v != 0) && ((v & (v - 1)) == 0);
66}
67
68#[inline]
69pub fn round_up_pow2(v: u32) -> u32 {
70    let mut v = v;
71    v -= 1;
72    v |= v >> 1;
73    v |= v >> 2;
74    v |= v >> 4;
75    v |= v >> 8;
76    v |= v >> 16;
77    return v + 1;
78}
79
80#[inline]
81pub fn is_power_of_2_64(v: u64) -> bool {
82    return (v != 0) && ((v & (v - 1)) == 0);
83}
84
85#[inline]
86pub fn round_up_pow2_64(v: u64) -> u64 {
87    let mut v = v;
88    v -= 1;
89    v |= v >> 1;
90    v |= v >> 2;
91    v |= v >> 4;
92    v |= v >> 8;
93    v |= v >> 16;
94    v |= v >> 32;
95    return v + 1;
96}
97
98#[inline]
99pub fn count_trailing_zeros(v: u32) -> u32 {
100    return v.trailing_zeros();
101}
102
103//, const Predicate &pred
104#[inline]
105pub fn find_interval(v: &[Float], pred: &dyn Fn(&[Float], usize) -> bool) -> usize {
106    let mut first = 0;
107    let mut len = v.len() as i64;
108    while len > 0 {
109        let half = len.wrapping_shr(1);
110        let middle = first + half;
111        // Bisect range based on value of _pred_ at _middle_
112        if pred(v, middle as usize) {
113            first = middle + 1;
114            len -= half + 1;
115        } else {
116            len = half;
117        }
118    }
119    return i64::clamp(first - 1, 0, v.len() as i64 - 2) as usize;
120}
121
122#[inline]
123pub fn find_interval_range(
124    range: (usize, usize),
125    v: &[Float],
126    pred: &dyn Fn(&[Float], usize) -> bool,
127) -> usize {
128    let mut first = range.0;
129    let mut len = range.1;
130    while len > 0 {
131        let half = len.wrapping_shr(1);
132        let middle = first + half;
133        // Bisect range based on value of _pred_ at _middle_
134        if pred(v, middle) {
135            first = middle + 1;
136            len -= half + 1;
137        } else {
138            len = half;
139        }
140    }
141    return usize::clamp(first - 1, 0, range.1 - 2);
142}
143
144#[inline]
145pub fn log2int(v: u32) -> u32 {
146    return 31 - v.leading_zeros();
147}
148
149#[inline]
150pub fn log2_int_64(v: u64) -> u32 {
151    return 63 - v.leading_zeros();
152}
153
154#[inline]
155pub fn erf_inv(x: Float) -> Float {
156    let x = Float::clamp(x, -0.99999, 0.99999);
157    let mut w = -Float::ln((1.0 - x) * (1.0 + x));
158    let mut p;
159    if w < 5.0 {
160        w -= 2.5;
161        p = 2.81022636e-08;
162        p = 3.43273939e-07 + p * w;
163        p = -3.5233877e-06 + p * w;
164        p = -4.39150654e-06 + p * w;
165        p = 0.00021858087 + p * w;
166        p = -0.00125372503 + p * w;
167        p = -0.00417768164 + p * w;
168        p = 0.246640727 + p * w;
169        p = 1.50140941 + p * w;
170    } else {
171        w = Float::sqrt(w) - 3.0;
172        p = -0.000200214257;
173        p = 0.000100950558 + p * w;
174        p = 0.00134934322 + p * w;
175        p = -0.00367342844 + p * w;
176        p = 0.00573950773 + p * w;
177        p = -0.0076224613 + p * w;
178        p = 0.00943887047 + p * w;
179        p = 1.00167406 + p * w;
180        p = 2.83297682 + p * w;
181    }
182    return p * x;
183}
184
185#[inline]
186pub fn erf(x: Float) -> Float {
187    // constants
188    const A1: Float = 0.254829592;
189    const A2: Float = -0.284496736;
190    const A3: Float = 1.421413741;
191    const A4: Float = -1.453152027;
192    const A5: Float = 1.061405429;
193    const P: Float = 0.3275911;
194
195    // Save the sign of x
196    let sign = if x >= 0.0 { 1.0 } else { -1.0 };
197    let x = Float::abs(x);
198
199    // A&S formula 7.1.26
200    let t = 1.0 / (1.0 + P * x);
201    let y = 1.0 - (((((A5 * t + A4) * t) + A3) * t + A2) * t + A1) * t * Float::exp(-x * x);
202
203    return sign * y;
204}
205
206#[cfg(test)]
207mod tests {
208    use super::*;
209
210    #[test]
211    pub fn test_001() {
212        let r1 = round_up_pow2(1);
213        let r2 = round_up_pow2(2);
214        let r3 = round_up_pow2(3);
215        let r4 = round_up_pow2(4);
216        let r5 = round_up_pow2(5);
217        assert_eq!(r1, 1);
218        assert_eq!(r2, 2);
219        assert_eq!(r3, 4);
220        assert_eq!(r4, 4);
221        assert_eq!(r5, 8);
222    }
223
224    #[test]
225    pub fn test_002() {
226        let r2 = log2int(2);
227        let r3 = log2int(3);
228        let r4 = log2int(4);
229        assert_eq!(r2, 1);
230        assert_eq!(r3, 1);
231        assert_eq!(r4, 2);
232    }
233}
pbrt_r3/core/base/functions.rs

pbrt_r3/core/base/
functions.rs
