Added Sobol sequence and filter

Cargo.toml

@@ -21,6 +21,7 @@ serde_json = "1.0.105"
 png = "0.17.10"
 wavefront_obj = "10.0.0"
 dotenv = "0.15.0"
+sobol_burley = "0.5.0"
 
 [dev-dependencies]
 test-log = "0.2.12"
@@ -30,8 +31,8 @@ criterion = "0.5.1"
 lto = true          # Enable link-time optimization
 codegen-units = 1   # Reduce number of codegen units to increase optimizations
 panic = 'abort'     # Abort on panic
-strip = true        # Strip symbols from binary*
-debug = false
+strip = false        # Strip symbols from binary*
+debug = true
 
 [profile.dev]
 strip = false

benches/default_benchmark.rs

@@ -92,6 +92,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
         samples: 4,
         max_depth: 8,
         sky_color: Color::new(0.5, 0.5, 0.5),
+        pixel_radius: 2.0,
         bvh_split_method: Some(BVHSplitMethod::Mid),
     };
 
@@ -101,6 +102,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
         samples: 4,
         max_depth: 8,
         sky_color: Color::new(0.5, 0.5, 0.5),
+        pixel_radius: 2.0,
         bvh_split_method: Some(BVHSplitMethod::SAH),
     };
 

input/config.json

@@ -1,12 +1,13 @@
 {
   "img_width": 600.0,
   "img_height": 600.0,
-  "samples": 8,
-  "max_depth": 12,
+  "samples": 12,
+  "max_depth": 8,
   "sky_color": {
     "r": 0.3,
     "g": 0.5,
     "b": 0.9
   },
-  "bvh_split_method": "Mid"
+  "pixel_radius": 2.0,
+  "bvh_split_method": "SAH"
 }

input/scene.json

@@ -38,9 +38,9 @@
           "theta_z": 0.0
         },
         "scale" : {
-          "x" : 20.0,
-          "y" : 20.0,
-          "z" : 20.0
+          "x" : 22.0,
+          "y" : 22.0,
+          "z" : 22.0
         }
       }
     },

src/bvh/bvh_node_sah.rs

@@ -60,7 +60,6 @@ impl Buckets {
                 best_split = split;
                 best_cost = plane_split_cost;
             }
-
         }
 
         // set up the vertices for the left and right indices
@@ -84,8 +83,6 @@ impl Buckets {
 /// BVH single bucket struct
 #[derive(Debug, Clone)]
 struct Bucket {
-    /// Interval for the centroids of the elements in this buckt
-    interval: Interval,
 
     /// Indices of the objects that are part of this bucket
     indices: Vec<usize>,
@@ -105,6 +102,7 @@ impl Bucket {
 }
 
 /// BVH Node for the BHV tree (node or leaf), and it consists of references to objects with a hittable trait (elements)
+#[derive(Debug)]
 #[allow(non_camel_case_types)]
 pub enum BVHNode_SAH {
     Leaf {
@@ -140,6 +138,7 @@ pub enum BVHNode_SAH {
 
 /// BVH tree struct
 #[allow(non_camel_case_types)]
+#[derive(Debug)]
 pub struct BVH_SAH<'a> {
     /// The list of nodes of the BVH.
     pub nodes: Vec<BVHNode_SAH>,
@@ -396,35 +395,25 @@ impl BVHNode_SAH {
                         }
                         BVHSplitMethod::SAH => {
                             // define the number of buckets and define the size of one bucket
-                            const NUM_BUCKETS: usize = 12;
+                            const NUM_BUCKETS: usize = 16;
                             let mut bucket_vec: Vec<Bucket> = Vec::with_capacity(NUM_BUCKETS);
-                            let bucket_size: f64 = sort_axis_size / NUM_BUCKETS as f64;
 
                             // create the buckets, with the right interval
-                            for i in 1..NUM_BUCKETS + 1 {
-
+                            for _ in 1..NUM_BUCKETS + 1 {
                                 bucket_vec.push(Bucket {
                                     //set the interval for the bucket, starting at the min of the centroid Aabb
-                                    interval: Interval::new(
-                                        aabb_centroids.min[sort_axis] + 
-                                            ((i - 1) as f64) * bucket_size,
-                                        aabb_centroids.min[sort_axis] + (i as f64 * bucket_size),
-                                    ),
-                                    indices: Vec::new(),
+                                    indices: Vec::with_capacity(indices.len() / NUM_BUCKETS),
                                 });
                             }
 
-                            // walk through each object, and check if its within the interval of each bucket
+                            let k0: f64 = aabb_centroids.min[sort_axis];
+                            let k1: f64 = (NUM_BUCKETS as f64 * (1.0 - 0.0001)) / (aabb_centroids.max[sort_axis] - aabb_centroids.min[sort_axis]);
+
+                            // walk through each object, check the bin number and add to the bucket
                             for index in indices {
-                                for b in &mut bucket_vec {
-                                    // check if this object is in the interval of this bucket
-                                    if b.interval.contains(
-                                        objects[*index].bounding_box().centroid()[sort_axis],
-                                    ) {
-                                        //add the elment to the bucket (this will also expand the bucket bounding box)
-                                        b.indices.push(*index);
-                                    }
-                                }
+                                let bin_id: usize = (k1 * (objects[*index].bounding_box().centroid[sort_axis]-k0)) as usize;
+                                bucket_vec[bin_id].indices.push(*index);
+
                             }
 
                             // add it to the bucket struct, so we can apply the split method

src/linalg/vec3.rs

@@ -53,6 +53,7 @@ impl Vec3 {
     }
 
     //creates a normalized random vector (on unit sphere) on the same hemisphere as given vector
+    #[inline(always)]
     pub fn new_random_on_hemisphere(normal: &Vec3, rng: &mut impl Rng) -> Self {
         let rv = Vec3::new_random_unit_vector(rng);
         if rv.dot(normal) > 0.0 {
@@ -64,6 +65,7 @@ impl Vec3 {
     }
 
     //creates a new random cosine pdf weigthed vector
+    #[inline(always)]
     pub fn random_cosine_direction(rng: &mut impl Rng) -> Self {
         let r1: f64 = rng.gen_range(0.0..1.0);
         let r2: f64 = rng.gen_range(0.0..1.0);
@@ -82,18 +84,22 @@ impl Vec3 {
     }
 
     // Exposing the x, y, z coordinates through functions
+    #[inline(always)]
     pub fn x(&self) -> f64 {
         self.x
     }
 
+    #[inline(always)]
     pub fn y(&self) -> f64 {
         self.y
     }
 
+    #[inline(always)]
     pub fn z(&self) -> f64 {
         self.z
     }
 
+    #[inline(always)]
     pub fn axis(&self, n: usize) -> f64 {
         match n {
             1 => self.y,
@@ -106,6 +112,7 @@ impl Vec3 {
     // calculate the Euclidean distance between two points
     // d(p, w) = sqrt((px - wx)^2 + (py - wy)^2 + (pz - wz)^2)
     // TODO add benchmark and replace with x * x to see if its
+    #[inline(always)]
     pub fn distance(&self, w: &Vec3) -> f64 {
         let dx = (self.x - w.x()).powf(2.0);
         let dy = f64::powf(self.y - w.y(), 2.0);
@@ -115,28 +122,33 @@ impl Vec3 {
 
     // fn length
     // calculate the vector (length) between point p and the origin (0,0,0)
+    #[inline(always)]
     pub fn length(&self) -> f64 {
         self.distance(&Vec3::new(0.0, 0.0, 0.0))
     }
 
     // fn lenght
     // calculates the length, squared
+    #[inline(always)]
     pub fn length_squared(&self) -> f64 {
         self.x * self.x + self.y * self.y + self.z * self.z
     }
 
     // fn abs
     // returns absolute value of vector, which is the same as its length (distance from 0.0)
+    #[inline(always)]
     pub fn abs(&self) -> f64 {
         self.length()
     }
 
     // vector dot multiplication of self with a given vector
+    #[inline(always)]
     pub fn dot(&self, w: &Vec3) -> f64 {
         self.x * w.x + self.y * w.y + self.z * w.z
     }
 
     // vector cross product
+    #[inline(always)]
     pub fn cross(&self, other: &Vec3) -> Vec3 {
         Vec3::new(
             self.y * other.z - self.z * other.y,
@@ -147,12 +159,14 @@ impl Vec3 {
 
     // fn normalized
     // return the normalized vector (= unit vector)
+    #[inline(always)]
     pub fn normalized(&self) -> Self {
         let l = self.length();
         Vec3::new(self.x() / l, self.y() / l, self.z() / l)
     }
 
     // returns a vector that is made from the smallest components of the two vectors
+    #[inline(always)]
     pub fn min(&self, other: Vec3) -> Self {
         Vec3 {
             x: self.x().min(other.x()),
@@ -162,6 +176,7 @@ impl Vec3 {
     }
 
     // returns a vector that is made from the smallest components of the two vectors
+    #[inline(always)]
     pub fn max(&self, other: Vec3) -> Self {
         Vec3 {
             x: self.x().max(other.x()),

src/main.rs

@@ -56,6 +56,7 @@ mod tests {
             samples: 8,
             max_depth: 8,
             sky_color: Color::new(0.5, 0.5, 0.5),
+            pixel_radius: 2.0,
             bvh_split_method: Some(BVHSplitMethod::Mid),
         };
         let json_camera: JSONCamera = JSONCamera {

src/render/camera.rs

@@ -92,22 +92,22 @@ impl Camera {
 
     // Get a randomly-sampled camera ray for the pixel at location x,i, originating from
     // the camera defocus disk.
-    pub fn get_prime_ray(self, x: i64, y: i64, rng: &mut impl Rng) -> Ray {
+    pub fn get_prime_ray(self, x: usize, y: usize, rng: &mut impl Rng, sample_position: Vec<f64>) -> Ray {
         // get a random point in the camera defocus disk to be used as an origin, or the camera center if the defocus angle <= 0 (this is for focus blur)
-        let p: Vec3 = Vec3::new_random_in_unit_disk(rng);
-        let disk_sample: Vec3 =
-            self.camera_center + (self.defocus_disk_u * p.x()) + (self.defocus_disk_v * p.y());
+
         let ray_origin: Vec3 = if self.defocus_angle <= 0.0 {
             self.camera_center
         } else {
+            let p: Vec3 = Vec3::new_random_in_unit_disk(rng);
+            let disk_sample: Vec3 =
+                self.camera_center + (self.defocus_disk_u * p.x()) + (self.defocus_disk_v * p.y());
             disk_sample
         };
-        //println!("origin {:?}", ray_origin);
 
         // establish the the direction of the ray by taking a random sample from the square around the pixel center (this is for anti aliasing)
         let pixel_loc: Vec3 =
             self.pixel00_loc + (self.pixel_delta_u * x as f64) + (self.pixel_delta_v * y as f64);
-        let pixel_sample: Vec3 = pixel_loc + self.pixel_sample_square(rng);
+        let pixel_sample: Vec3 = pixel_loc + Vec3::new(sample_position[0], sample_position[1], 0.);
         let ray_direction: Vec3 = pixel_sample - ray_origin;
 
         Ray::new(ray_origin, ray_direction)
@@ -141,6 +141,7 @@ mod tests {
             samples: 1,
             max_depth: 32,
             sky_color: Color::new(3.0 / 255.0, 165.0 / 255.0, 252.0 / 255.0),
+            pixel_radius: 2.0,
             bvh_split_method: Some(crate::bvh::BVHSplitMethod::Mid),
         };
 

src/render/config.rs

@@ -12,6 +12,7 @@ pub struct Config {
     pub samples: usize,
     pub max_depth: usize,
     pub sky_color: Color,
+    pub pixel_radius: f64,
     pub bvh_split_method: Option<BVHSplitMethod>,
 }
 
@@ -35,6 +36,7 @@ impl Default for Config {
             max_depth: m,
             sky_color,
             bvh_split_method: Some(bvh_split_method),
+            pixel_radius: 2.0,
         }
     }
 }

src/render/filter.rs

@@ -0,0 +1,77 @@
+
+/// Enum for different filter classes
+#[derive(Debug, Clone)]
+pub enum Filter {
+    MitchellNetravali(MitchellNetravali),
+}
+
+impl Filter {
+    pub fn evaluate(&self, sample_x: f64, sample_y: f64) -> f64 {
+        match self {
+            Filter::MitchellNetravali(filter) => filter.evaluate(sample_x, sample_y),
+        }
+    }
+    pub fn get_radius(&self) -> (f64, f64) {
+        match self {
+            Filter::MitchellNetravali(filter) => filter.get_radius(),
+        }
+    }
+}
+
+/// Struct for MitchellNetravali filter
+#[derive(Debug, Default, Copy, Clone)]
+pub struct MitchellNetravali {
+    width: f64,
+    height: f64,
+    inv_width: f64,
+    inv_height: f64,
+    b: f64,
+    c: f64,
+}
+
+impl MitchellNetravali {
+
+    /// Returns a new MN filter
+    pub fn new(w: f64, h: f64, b: f64, c: f64) -> Self {
+        MitchellNetravali {
+            width: w,
+            height: h,
+            inv_width: 1.0 / w,
+            inv_height: 1.0 / h,
+            b,
+            c,
+        }
+    }
+
+    #[inline(always)]
+    fn mitchell_1d(&self, x: f64) -> f64 {
+        let fx = x.abs() * 2.0;
+        if fx < 1.0 {
+            ((12.0 - 9.0 * self.b - 6.0 * self.c) * fx * fx * fx
+                + (-18.0 + 12.0 * self.b + 6.0 * self.c) * fx * fx
+                + (6.0 - 2.0 * self.b))
+                * (1.0 / 6.0)
+        } else if fx < 2.0 {
+            ((-self.b - 6.0 * self.c) * fx * fx * fx
+                + (6.0 * self.b + 30.0 * self.c) * fx * fx
+                + (-12.0 * self.b - 48.0 * self.c) * fx
+                + (8.0 * self.b + 24.0 * self.c))
+                * (1.0 / 6.0)
+        } else {
+            0.0
+        }
+    }
+
+    /// Evaluates the MN filter for a given sample position
+    #[inline(always)]
+    pub fn evaluate(&self, sample_x: f64, sample_y: f64) -> f64 {
+        self.mitchell_1d(sample_x * self.inv_width) * self.mitchell_1d(sample_y * self.inv_height)
+    }
+
+    /// Returns the radius set for the MN filter
+    #[inline(always)]
+    pub fn get_radius(&self) -> (f64, f64) {
+        (self.width, self.height)
+    }
+
+}