pbrt_r3/core/transform/
interval.rs

1use crate::core::base::*;
2
3use std::ops;
4
5#[derive(Debug, Clone, Copy)]
6pub struct Interval {
7    pub low: Float,
8    pub high: Float,
9}
10
11impl Interval {
12    pub fn new(low: Float, high: Float) -> Self {
13        Self { low, high }
14    }
15
16    pub fn sin(i: Interval) -> Self {
17        assert!(0.0 <= i.low);
18        assert!(i.high <= 2.0001 * PI);
19
20        let mut sin_low = Float::sin(i.low);
21        let mut sin_high = Float::sin(i.high);
22        if sin_low > sin_high {
23            std::mem::swap(&mut sin_low, &mut sin_high);
24        }
25        if i.low < PI / 2.0 && i.high > PI / 2.0 {
26            sin_high = 1.;
27        }
28        if i.low < (3.0 / 2.0) * PI && i.high > (3.0 / 2.0) * PI {
29            sin_low = -1.0;
30        }
31        return Interval::new(sin_low, sin_high);
32    }
33
34    pub fn cos(i: Interval) -> Self {
35        assert!(0.0 <= i.low);
36        assert!(i.high <= 2.0001 * PI);
37
38        let mut cos_low = Float::cos(i.low);
39        let mut cos_high = Float::cos(i.high);
40        if cos_low > cos_high {
41            std::mem::swap(&mut cos_low, &mut cos_high);
42        }
43        if i.low < PI && i.high > PI {
44            cos_low = -1.0;
45        }
46        return Interval::new(cos_low, cos_high);
47    }
48}
49
50impl From<Float> for Interval {
51    fn from(f: Float) -> Self {
52        return Interval::new(f, f);
53    }
54}
55
56impl ops::Add<Interval> for Interval {
57    type Output = Interval;
58    fn add(self, rhs: Interval) -> Self::Output {
59        return Interval::new(self.low + rhs.low, self.high + rhs.high);
60    }
61}
62
63impl ops::Sub<Interval> for Interval {
64    type Output = Interval;
65    fn sub(self, rhs: Interval) -> Self::Output {
66        return Interval::new(self.low - rhs.low, self.high - rhs.high);
67    }
68}
69
70impl ops::Mul<Interval> for Interval {
71    type Output = Interval;
72    fn mul(self, rhs: Interval) -> Self::Output {
73        let albl = self.low * rhs.low;
74        let ahbl = self.high * rhs.low;
75        let albh = self.low * rhs.high;
76        let ahbh = self.high * rhs.high;
77        let min = Float::min(Float::min(albl, ahbl), Float::min(albh, ahbh));
78        let max = Float::max(Float::max(albl, ahbl), Float::max(albh, ahbh));
79        return Interval::new(min, max);
80    }
81}
82
83pub fn interval_find_zeros_(
84    c1: Float,
85    c2: Float,
86    c3: Float,
87    c4: Float,
88    c5: Float,
89    theta: Float,
90    t: Interval,
91    depth: u32,
92    zeros: &mut Vec<Float>,
93) {
94    // Evaluate motion derivative in interval form, return if no zeros
95    let range = Interval::from(c1)
96        + (Interval::from(c2) + Interval::from(c3) * t)
97            * Interval::cos(Interval::from(2.0 * theta) * t)
98        + (Interval::from(c4) + Interval::from(c5) * t)
99            * Interval::sin(Interval::from(2.0 * theta) * t);
100    if range.low > 0.0 || range.high < 0.0 || range.low == range.high {
101        return;
102    }
103    if depth > 0 {
104        // Split the interval and recurse
105        let tm = 0.5 * (t.low + t.high);
106        interval_find_zeros_(
107            c1,
108            c2,
109            c3,
110            c4,
111            c5,
112            theta,
113            Interval::new(t.low, tm),
114            depth - 1,
115            zeros,
116        );
117        interval_find_zeros_(
118            c1,
119            c2,
120            c3,
121            c4,
122            c5,
123            theta,
124            Interval::new(tm, t.high),
125            depth - 1,
126            zeros,
127        );
128    } else {
129        // Use Newton's method to refine zero
130        let mut t_newton = (t.low + t.high) * 0.5;
131        for _i in 0..4 {
132            let f_newton = c1
133                + (c2 + c3 * t_newton) * Float::cos(2.0 * theta * t_newton)
134                + (c4 + c5 * t_newton) * Float::sin(2.0 * theta * t_newton);
135            let f_prime_newton = (c3 + 2.0 * (c4 + c5 * t_newton) * theta)
136                * Float::cos(2.0 * theta * t_newton)
137                + (c5 - 2.0 * (c2 + c3 * t_newton) * theta) * Float::sin(2.0 * theta * t_newton);
138            if f_newton == 0.0 || f_prime_newton == 0.0 {
139                break;
140            }
141            t_newton -= f_newton / f_prime_newton;
142        }
143        if t_newton >= (t.low - 1e-3) && t_newton < (t.high + 1e-3) {
144            zeros.push(t_newton);
145        }
146    }
147}
148
149pub fn interval_find_zeros(
150    c1: Float,
151    c2: Float,
152    c3: Float,
153    c4: Float,
154    c5: Float,
155    theta: Float,
156    t: Interval,
157    depth: u32,
158) -> Vec<Float> {
159    let mut zeros = Vec::new();
160    interval_find_zeros_(c1, c2, c3, c4, c5, theta, t, depth, &mut zeros);
161    return zeros;
162}
pbrt_r3/core/transform/interval.rs

pbrt_r3/core/transform/
interval.rs
