pbrt_r3/samplers/
halton.rs

1use crate::core::lowdiscrepancy::primes::{PRIMES, PRIME_SUMS};
2use crate::core::prelude::*;
3
4use std::cell::Cell;
5use std::sync::Arc;
6use std::sync::OnceLock;
7use std::sync::RwLock;
8
9const K_MAX_RESOLUTION: i32 = 128;
10static RADICAL_INVERSE_PERMUTATIONS: OnceLock<Vec<u16>> = OnceLock::new();
11
12fn get_radical_inverse_permutations() -> &'static Vec<u16> {
13    let v = RADICAL_INVERSE_PERMUTATIONS.get_or_init(|| -> Vec<u16> {
14        let mut v = Vec::new();
15        let mut rng = RNG::new();
16        compute_radical_inverse_permutations(&mut v, &mut rng);
17        return v;
18    });
19    return v;
20}
21
22fn math_mod(a: i64, b: i64) -> u64 {
23    let result = a - (a / b) * b;
24    return if result < 0 {
25        (result + b) as u64
26    } else {
27        result as u64
28    };
29}
30
31fn extended_gcd(a: u64, b: u64) -> (i64, i64) {
32    if b == 0 {
33        return (1, 0);
34    }
35    let d = (a / b) as i64;
36    let (xp, yp) = extended_gcd(b, a % b);
37    return (yp, xp - (d * yp));
38}
39
40fn multiplicative_inverse(a: i64, n: i64) -> u64 {
41    let (x, _y) = extended_gcd(a as u64, n as u64);
42    return math_mod(x, n);
43}
44
45#[derive(Debug, Default, PartialEq, Clone)]
46pub struct HaltonSampler {
47    base: BaseGlobalSampler,
48    base_scales: [i32; 2],
49    base_exponents: [i32; 2],
50    sample_stride: i32,
51    mult_inverse: [i32; 2],
52    pixel_for_offset: Cell<Point2i>,
53    offset_for_current_pixel: Cell<i64>,
54    sample_at_pixel_center: bool,
55}
56
57impl HaltonSampler {
58    pub fn new(samples_per_pixel: u32, sample_bounds: &Bounds2i, sample_at_center: bool) -> Self {
59        let _ = get_radical_inverse_permutations();
60
61        // Find radical inverse base scales and exponents that cover sampling area
62        let res = sample_bounds.max - sample_bounds.min;
63        const BASES: [i32; 2] = [2, 3];
64
65        let mut base_scales = [1, 1];
66        let mut base_exponents = [1, 1];
67        for i in 0..2 {
68            //let base = if i == 0 { 2 } else { 3 };
69            let base = BASES[i];
70            let mut scale = 1;
71            let mut exp = 0;
72            while scale < i32::min(res[i], K_MAX_RESOLUTION) {
73                scale *= base;
74                exp += 1;
75            }
76            base_scales[i] = scale;
77            base_exponents[i] = exp;
78        }
79
80        // Compute stride in samples for visiting each pixel area
81        let sample_stride = base_scales[0] * base_scales[1];
82        let mult_inverse = [
83            multiplicative_inverse(base_scales[1] as i64, base_scales[0] as i64) as i32,
84            multiplicative_inverse(base_scales[0] as i64, base_scales[1] as i64) as i32,
85        ];
86
87        let pixel_for_offset = Cell::new(Point2i::new(i32::MAX, i32::MAX));
88        let offset_for_current_pixel = Cell::new(0);
89        let base = BaseGlobalSampler::new(samples_per_pixel);
90
91        HaltonSampler {
92            base,
93            base_scales,
94            base_exponents,
95            sample_stride,
96            mult_inverse,
97            pixel_for_offset,
98            offset_for_current_pixel,
99            sample_at_pixel_center: sample_at_center,
100        }
101    }
102
103    fn permutation_for_dimension(dim: u32) -> &'static [u16] {
104        let table_size = PRIMES.len();
105        if dim as usize >= table_size {
106            panic!("HaltonSampler can only sample {} dimensions.", table_size);
107        }
108        let start = PRIME_SUMS[dim as usize] as usize;
109        let rip = get_radical_inverse_permutations();
110        return &rip[start..];
111    }
112}
113
114impl GlobalSampler for HaltonSampler {
115    fn get_index_for_sample(&self, sample_num: i64) -> i64 {
116        let sample_stride = self.sample_stride as i64;
117        if self.base.base.current_pixel != self.pixel_for_offset.get() {
118            // Compute Halton sample offset for _currentPixel_
119            self.offset_for_current_pixel.set(0);
120            if sample_stride > 1 {
121                let pm = [
122                    math_mod(
123                        self.base.base.current_pixel[0] as i64,
124                        K_MAX_RESOLUTION as i64,
125                    ),
126                    math_mod(
127                        self.base.base.current_pixel[1] as i64,
128                        K_MAX_RESOLUTION as i64,
129                    ),
130                ];
131                const BASES: [u64; 2] = [2, 3];
132                for i in 0..2 {
133                    let dim_offset =
134                        inverse_radical_inverse(BASES[i], pm[i], self.base_exponents[i] as usize);
135                    let addition = (dim_offset
136                        * ((self.sample_stride / self.base_scales[i]) * self.mult_inverse[i])
137                            as u64) as i64;
138                    self.offset_for_current_pixel
139                        .set(self.offset_for_current_pixel.get() + addition);
140                }
141                let offset_for_current_pixel = self.offset_for_current_pixel.get() % sample_stride;
142                self.offset_for_current_pixel.set(offset_for_current_pixel);
143            }
144            self.pixel_for_offset.set(self.base.base.current_pixel);
145        }
146        return self.offset_for_current_pixel.get() + sample_num * sample_stride;
147    }
148
149    fn sample_dimension(&self, index: i64, dim: u32) -> Float {
150        if self.sample_at_pixel_center && (dim == 0 || dim == 1) {
151            return 0.5;
152        }
153        return match dim {
154            0 => radical_inverse(dim, (index >> self.base_exponents[0]) as u64),
155            1 => radical_inverse(dim, (index / self.base_scales[1] as i64) as u64),
156            _ => scrambled_radical_inverse(
157                dim,
158                index as u64,
159                HaltonSampler::permutation_for_dimension(dim),
160            ),
161        };
162    }
163
164    fn get_base(&mut self) -> &mut BaseGlobalSampler {
165        return &mut self.base;
166    }
167}
168
169impl Sampler for HaltonSampler {
170    fn start_pixel(&mut self, p: &Point2i) {
171        let _p = ProfilePhase::new(Prof::StartPixel);
172
173        self.base.base.start_pixel(p);
174        self.base.dimension = 0;
175        self.base.interval_sample_index = self.get_index_for_sample(0);
176        self.base.array_end_dim = BaseGlobalSampler::array_start_dim
177            + (self.base.base.sample_array1d.len() + 2 * self.base.base.sample_array2d.len())
178                as u32;
179        let samples1d_array_sizes = self.base.base.samples1d_array_sizes.len();
180        // Compute 1D array samples for _GlobalSampler_
181        {
182            let l = samples1d_array_sizes;
183            for i in 0..l {
184                let n_samples =
185                    self.base.base.samples1d_array_sizes[i] * self.base.base.samples_per_pixel;
186                for j in 0..n_samples {
187                    let index = self.get_index_for_sample(j as i64);
188                    self.base.base.sample_array1d[i][j as usize] =
189                        self.sample_dimension(index, self.base.array_end_dim + i as u32);
190                }
191            }
192        }
193        // Compute 2D array samples for _GlobalSampler_
194        {
195            let dim = BaseGlobalSampler::array_start_dim + samples1d_array_sizes as u32;
196            let l = self.base.base.samples2d_array_sizes.len();
197            for i in 0..l {
198                let n_samples =
199                    self.base.base.samples2d_array_sizes[i] * self.base.base.samples_per_pixel;
200                for j in 0..n_samples {
201                    let index = self.get_index_for_sample(j as i64);
202                    self.base.base.sample_array2d[i][j as usize].x =
203                        self.sample_dimension(index, dim);
204                    self.base.base.sample_array2d[i][j as usize].y =
205                        self.sample_dimension(index, dim + 1);
206                }
207            }
208        }
209    }
210
211    fn start_next_sample(&mut self) -> bool {
212        self.base.dimension = 0;
213        self.base.interval_sample_index =
214            self.get_index_for_sample((self.base.base.current_pixel_sample_index + 1) as i64);
215        return self.base.base.start_next_sample();
216    }
217
218    fn set_sample_number(&mut self, sample_num: u32) -> bool {
219        self.base.dimension = 0;
220        self.base.interval_sample_index = self.get_index_for_sample(sample_num as i64);
221        return self.base.base.set_sample_number(sample_num);
222    }
223
224    fn get_1d(&mut self) -> Float {
225        let _p = ProfilePhase::new(Prof::GetSample);
226
227        if self.base.dimension >= BaseGlobalSampler::array_start_dim
228            && self.base.dimension < self.base.array_end_dim
229        {
230            self.base.dimension = self.base.array_end_dim;
231        }
232        let d = self.base.dimension;
233        let x = self.sample_dimension(self.base.interval_sample_index, d);
234        self.base.dimension += 1;
235        return x;
236    }
237
238    fn get_2d(&mut self) -> Vector2f {
239        let _p = ProfilePhase::new(Prof::GetSample);
240
241        if self.base.dimension + 1 >= BaseGlobalSampler::array_start_dim
242            && self.base.dimension < self.base.array_end_dim
243        {
244            self.base.dimension = self.base.array_end_dim;
245        }
246        let d = self.base.dimension;
247        let x = self.sample_dimension(self.base.interval_sample_index, d);
248        let y = self.sample_dimension(self.base.interval_sample_index, d + 1);
249        let p = Point2f::new(x, y);
250        self.base.dimension += 2;
251        return p;
252    }
253
254    fn request_1d_array(&mut self, n: u32) {
255        self.base.request_1d_array(n);
256    }
257    fn request_2d_array(&mut self, n: u32) {
258        self.base.request_2d_array(n);
259    }
260    fn get_1d_array(&mut self, n: u32) -> Option<Vec<Float>> {
261        return self.base.get_1d_array(n);
262    }
263    fn get_2d_array(&mut self, n: u32) -> Option<Vec<Vector2f>> {
264        return self.base.get_2d_array(n);
265    }
266
267    fn clone_with_seed(&self, _seed: u32) -> Arc<RwLock<dyn Sampler>> {
268        return Arc::new(RwLock::new(self.clone()));
269    }
270    fn get_samples_per_pixel(&self) -> u32 {
271        return self.base.base.samples_per_pixel;
272    }
273}
274
275pub fn create_halton_sampler(
276    params: &ParamSet,
277    sample_bounds: &Bounds2i,
278) -> Result<Arc<RwLock<dyn Sampler>>, PbrtError> {
279    let mut nsamp = params.find_one_int("pixelsamples", 16) as u32;
280    let sample_at_center = params.find_one_bool("samplepixelcenter", false);
281
282    // pbrt-r3
283    let _ = get_radical_inverse_permutations();
284    // pbrt-r3
285
286    {
287        let options = PbrtOptions::get();
288        if options.quick_render {
289            nsamp = 1;
290        }
291    }
292
293    return Ok(Arc::new(RwLock::new(HaltonSampler::new(
294        nsamp,
295        sample_bounds,
296        sample_at_center,
297    ))));
298}
pbrt_r3/samplers/halton.rs

pbrt_r3/samplers/
halton.rs
