Switch to sobol sampler.

The important thing here is that I figured out how to use the
scrambling parameter properly to decorrelate pixels.  Using the
same approach as with halton (just adding an offset into the sequence)
is very slow with sobol, since moving into the higher samples is
more computationally expensive.  So using the scrambling parameter
instead was important.

src/hash.rs

@@ -3,7 +3,7 @@ pub fn hash_u32(n: u32, seed: u32) -> u32 {
     for _ in 0..3 {
         hash = hash.wrapping_mul(1_936_502_639);
         hash ^= hash.wrapping_shr(16);
-        hash = hash.wrapping_add(seed);
+        hash ^= seed;
     }
 
     hash
@@ -14,7 +14,7 @@ pub fn hash_u64(n: u64, seed: u64) -> u64 {
     for _ in 0..4 {
         hash = hash.wrapping_mul(32_416_190_071 * 314_604_959);
         hash ^= hash.wrapping_shr(32);
-        hash = hash.wrapping_add(seed);
+        hash ^= seed;
     }
 
     hash
@@ -28,7 +28,7 @@ pub fn hash_u32_to_f32(n: u32, seed: u32) -> f32 {
     for _ in 0..3 {
         hash = hash.wrapping_mul(1_936_502_639);
         hash ^= hash.wrapping_shr(16);
-        hash = hash.wrapping_add(seed);
+        hash ^= seed;
     }
     const INV_MAX: f32 = 1.0 / std::u32::MAX as f32;
 

src/renderer.rs

@@ -240,14 +240,12 @@ impl<'a> Renderer<'a> {
             // Generate light paths and initial rays
             for y in bucket.y..(bucket.y + bucket.h) {
                 for x in bucket.x..(bucket.x + bucket.w) {
-                    let scramble = hash_u32(((x as u32) << 16) ^ (y as u32), self.seed);
+                    let pix_id = pixel_id(x, y);
                     for si in 0..self.spp {
                         // Calculate image plane x and y coordinates
                         let (img_x, img_y) = {
-                            let filter_x =
+                            let filter_x = fast_logit(get_sample(4, si as u32, pix_id), 1.5) + 0.5;
-                                fast_logit(get_sample(4, si as u32, scramble), 1.5) + 0.5;
+                            let filter_y = fast_logit(get_sample(5, si as u32, pix_id), 1.5) + 0.5;
-                            let filter_y =
-                                fast_logit(get_sample(5, si as u32, scramble), 1.5) + 0.5;
                             let samp_x = (filter_x + x as f32) * cmpx;
                             let samp_y = (filter_y + y as f32) * cmpy;
                             ((samp_x - 0.5) * x_extent, (0.5 - samp_y) * y_extent)
@@ -257,15 +255,15 @@ impl<'a> Renderer<'a> {
                         let (path, ray) = LightPath::new(
                             &self.scene,
                             (x, y),
+                            pix_id,
                             (img_x, img_y),
                             (
-                                get_sample(0, si as u32, scramble),
+                                get_sample(0, si as u32, pix_id),
-                                get_sample(1, si as u32, scramble),
+                                get_sample(1, si as u32, pix_id),
                             ),
-                            get_sample(2, si as u32, scramble),
+                            get_sample(2, si as u32, pix_id),
-                            map_0_1_to_wavelength(get_sample(3, si as u32, scramble)),
+                            map_0_1_to_wavelength(get_sample(3, si as u32, pix_id)),
                             si as u32,
-                            scramble,
                         );
                         paths.push(path);
                         rays.push(ray, false);
@@ -371,8 +369,8 @@ pub struct LightPath {
     bounce_count: u32,
 
     pixel_co: (u32, u32),
-    lds_offset: u32,
+    pixel_id: u32,
-    lds_scramble: u32,
+    sample_number: u32, // Which sample in the LDS sequence this is.
     dim_offset: Cell<u32>,
     time: f32,
     wavelength: f32,
@@ -391,12 +389,12 @@ impl LightPath {
     fn new(
         scene: &Scene,
         pixel_co: (u32, u32),
+        pixel_id: u32,
         image_plane_co: (f32, f32),
         lens_uv: (f32, f32),
         time: f32,
         wavelength: f32,
-        lds_offset: u32,
+        sample_number: u32,
-        lds_scramble: u32,
     ) -> (LightPath, Ray) {
         (
             LightPath {
@@ -404,8 +402,8 @@ impl LightPath {
                 bounce_count: 0,
 
                 pixel_co: pixel_co,
-                lds_offset: lds_offset,
+                pixel_id: pixel_id,
-                lds_scramble: lds_scramble,
+                sample_number: sample_number,
                 dim_offset: Cell::new(6),
                 time: time,
                 wavelength: wavelength,
@@ -432,7 +430,7 @@ impl LightPath {
     fn next_lds_samp(&self) -> f32 {
         let dimension = self.dim_offset.get();
         self.dim_offset.set(dimension + 1);
-        get_sample(dimension, self.lds_offset, self.lds_scramble)
+        get_sample(dimension, self.sample_number, self.pixel_id)
     }
 
     fn next(
@@ -689,13 +687,26 @@ impl LightPath {
 #[inline(always)]
 fn get_sample(dimension: u32, i: u32, scramble: u32) -> f32 {
     use crate::hash::hash_u32_to_f32;
-    if dimension < halton::MAX_DIMENSION {
+    if dimension < sobol::NUM_DIMENSIONS as u32 {
-        halton::sample(dimension, i + scramble)
+        sobol::sample_with_scramble(dimension, i, hash_u32(dimension, scramble))
     } else {
         hash_u32_to_f32(dimension, i + scramble)
     }
 }
 
+/// Make a unique-ish pixel ID, given its coordinates.
+///
+/// This is a pure, deterministic function.
+#[inline(always)]
+fn pixel_id(x: u32, y: u32) -> u32 {
+    // Pretend the image is 65536 x 65536 resolution,
+    // and do a mapping to a 1d array.  This should produce
+    // unique numbers for any reasonably sized image, and
+    // even for stupid big images will produce sufficiently
+    // unique numbers for our purposes.
+    y * (1 << 16) + x
+}
+
 #[derive(Debug)]
 struct BucketJob {
     x: u32,

sub_crates/sobol/src/matrices.rs

@@ -33,7 +33,7 @@
 pub const NUM_DIMENSIONS: usize = 1024;
 pub const SIZE: usize = 52;
 
-pub const MATRICES: [u32; NUM_DIMENSIONS * SIZE] = [
+pub const MATRICES: &[u32] = &[
     0x80000000u32,
     0x40000000u32,
     0x20000000u32,