From aba177c536e4c70035a1e280db44de8f398da624 Mon Sep 17 00:00:00 2001
From: Nathan Vegdahl <cessen@cessen.com>
Date: Tue, 16 Aug 2022 23:38:52 -0700
Subject: [PATCH] Use unscrambled Van der Corput for sampling time dimension.

This gives notably better results at powers of two.

Also experimenting with Halton bases 3 and 4 for DoF.
---
 src/renderer.rs | 105 ++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 98 insertions(+), 7 deletions(-)

diff --git a/src/renderer.rs b/src/renderer.rs
index 6fc18d4..ffbcb7b 100644
--- a/src/renderer.rs
+++ b/src/renderer.rs
@@ -232,16 +232,22 @@ impl<'a> Renderer<'a> {
                     // to work well with golden-ratio sampling.  Since we only use sobol
                     // and golden ratio sampling, we do this up-front here rather than in
                     // the samplers themselves.
-                    let si_offset =
-                        owen4_fast(morton::encode(x, y), self.seed).wrapping_mul(self.spp as u32);
+                    let pixel_offset = owen4_fast(morton::encode(x, y), self.seed);
+                    let si_offset = pixel_offset.wrapping_mul(self.spp as u32);
 
                     for si in 0..(self.spp as u32) {
-                        let mut sample_gen = SampleGen::new(si_offset + si, self.seed);
+                        let sii = si_offset.wrapping_add(si);
+                        let mut sample_gen = SampleGen::new(sii, self.seed);
 
                         // Raw sample numbers.
-                        let d0 = golden_ratio_sample(si.wrapping_add(si_offset));
-                        let [d1, d2, d3, _] = sample_gen.next_dims();
-                        sample_gen.inc_padding();
+                        let d0 = golden_ratio_sample(
+                            sii.wrapping_add(rrand::mix32_seed(self.seed, 1111)),
+                        );
+                        let d1 = van_der_corput(sii, rrand::mix32_seed(self.seed, 2222));
+                        let d2 = halton_3(sii, rrand::mix32_seed(self.seed, 3333));
+                        let d3 = halton_5(sii, rrand::mix32_seed(self.seed, 4444));
+                        // let [d2, d3, _, _] = sample_gen.next_dims();
+                        // sample_gen.inc_padding();
                         let [d4, d5, _, _] = sample_gen.next_dims();
                         sample_gen.inc_padding();
 
@@ -584,7 +590,92 @@ impl SampleGen {
     }
 }
 
-/// Golden ratio sampler.
+/// The Van der Corput sequence.
+///
+/// NOTE: use this for motion blur, because the jittering introduced by
+/// Owen Scrambling in the main sampler actually causes banding artifacts
+/// at the sample cell borders, and actually *increases* noise.
+fn van_der_corput(i: u32, seed: u32) -> f32 {
+    use sobol_burley::parts::u32_to_f32_norm;
+    u32_to_f32_norm(i.reverse_bits().wrapping_add(seed))
+}
+
+/// Halton base 3.
+fn halton_3(mut i: u32, seed: u32) -> f32 {
+    const BASE: u32 = 3;
+    const POW: u32 = 13;
+    const TOP_BASE: u32 = {
+        let mut b = 1;
+        let mut i = 0;
+        while i < POW {
+            b *= BASE;
+            i += 1;
+        }
+        b
+    };
+
+    let mut result = 0;
+    for _ in 0..POW {
+        result *= BASE;
+        let next = i / BASE;
+        result += i - (next * BASE);
+        i = next;
+    }
+
+    result = ((result as u64 + seed as u64) % TOP_BASE as u64) as u32;
+
+    result as f32 * (1.0 / TOP_BASE as f32)
+}
+
+/// Halton base 5.
+fn halton_5(mut i: u32, seed: u32) -> f32 {
+    const BASE: u32 = 5;
+    const POW: u32 = 13;
+    const TOP_BASE: u32 = {
+        let mut b = 1;
+        let mut i = 0;
+        while i < POW {
+            b *= BASE;
+            i += 1;
+        }
+        b
+    };
+    const ORDER: [u32; 5] = [0, 3, 1, 4, 2];
+
+    let mut result = 0;
+    for _ in 0..POW {
+        result *= BASE;
+        let next = i / BASE;
+        result += ORDER[(i - (next * BASE)) as usize];
+        i = next;
+    }
+
+    result = ((result as u64 + seed as u64) % TOP_BASE as u64) as u32;
+
+    result as f32 * (1.0 / TOP_BASE as f32)
+}
+
+// fn next_power_of_three(mut n: u64) -> u64 {
+//     n = n.saturating_sub(1);
+//     let mut i = 1;
+//     while n > 0 {
+//         n /= 3;
+//         i *= 3;
+//     }
+//     i
+// }
+
+// fn next_power_of_five(mut n: u64) -> u64 {
+//     n = n.saturating_sub(1);
+//     let mut i = 1;
+//     while n > 0 {
+//         n /= 5;
+//         i *= 5;
+//     }
+//     i
+// }
+
+/// The golden ratio sequence.
 ///
 // NOTE: use this for the wavelength dimension, because
 // due to the nature of hero wavelength sampling this ends up