diff --git a/src/camera/mod.rs b/src/camera/mod.rs
index 1fba116..1465be7 100644
--- a/src/camera/mod.rs
+++ b/src/camera/mod.rs
@@ -4,11 +4,12 @@ use crate::geometry::HitRecord;
 use crate::geometry::Hittable;
 use crate::interval::Interval;
 use crate::material::Surface;
-use crate::math;
 use crate::ray::Ray;
 use image::RgbImage;
+use itertools::iproduct;
 use nalgebra::Vector3;
 use rand::prelude::*;
+use rand::rngs::ThreadRng;
 use std::f64::consts::PI;
 use std::io::{self, Write};
 use std::sync::Arc;
@@ -22,11 +23,13 @@ pub struct Camera {
     pub image_width: u32,
     pub image_height: u32,
 
+    pub sqrt_spp: u32,
+    pub recip_sqrt_spp: f64,
+
     pub center: Vector3<f64>,
     pub pixel00_loc: Vector3<f64>,
     pub pixel_delta_u: Vector3<f64>,
     pub pixel_delta_v: Vector3<f64>,
-    pub samples: u32,
     pub samples_scale: f64,
     pub max_bounces: u32,
     pub threads: usize,
@@ -38,9 +41,7 @@ pub struct Camera {
     pub background: Vector3<f64>,
 }
 
-fn random_in_unit_disk() -> Vector3<f64> {
-    let mut rng = rand::rng();
-
+fn random_in_unit_disk(rng: &mut ThreadRng) -> Vector3<f64> {
     loop {
         let p = Vector3::new(
             rng.random_range(-1.0f64..1.0f64),
@@ -65,6 +66,21 @@ fn degrees_to_radians(degress: f64) -> f64 {
     degress * PI / 180.0
 }
 
+fn sample_square_stratified(
+    s_i: u32,
+    s_j: u32,
+    recip_sqrt_spp: f64,
+    rng: &mut ThreadRng,
+) -> (f64, f64) {
+    let s_i = s_i as f64;
+    let s_j = s_j as f64;
+
+    (
+        ((s_i + rng.random::<f64>()) * recip_sqrt_spp) - 0.5,
+        ((s_j + rng.random::<f64>()) * recip_sqrt_spp) - 0.5,
+    )
+}
+
 impl Camera {
     pub fn new(options: CameraOptions) -> Self {
         let (image_width, image_height) = options.get_dimensions();
@@ -105,7 +121,9 @@ impl Camera {
         let defocus_disk_v = v * defocus_radius;
 
         let samples = options.samples;
-        let samples_scale = 1.0 / (samples as f64);
+        let sqrt_spp = (samples as f64).sqrt() as u32;
+        let recip_sqrt_spp = 1.0 / (sqrt_spp as f64);
+        let samples_scale = 1.0 / ((sqrt_spp * sqrt_spp) as f64);
 
         let max_bounces = options.max_bounces;
         let threads = options.threads;
@@ -117,7 +135,6 @@ impl Camera {
             pixel00_loc,
             pixel_delta_u,
             pixel_delta_v,
-            samples,
             samples_scale,
             max_bounces,
             threads,
@@ -125,6 +142,8 @@ impl Camera {
             defocus_disk_v,
             defocus_angle,
             background,
+            sqrt_spp,
+            recip_sqrt_spp,
         }
     }
 
@@ -187,8 +206,16 @@ impl Camera {
     }
 
     pub fn get_pixel(&self, world: &Geometry, x: u32, y: u32) -> image::Rgb<u8> {
-        let color = (0..self.samples)
-            .map(|_| self.ray_color(&self.get_ray(x, y), self.max_bounces, world))
+        let mut rng = rand::rng();
+        let color: Vector3<f64> = iproduct!(0..self.sqrt_spp, 0..self.sqrt_spp)
+            .map(|(s_x, s_y)| {
+                self.ray_color(
+                    &self.get_ray(x, y, s_x, s_y, &mut rng),
+                    self.max_bounces,
+                    world,
+                    &mut rng,
+                )
+            })
             .sum::<Vector3<f64>>()
             * self.samples_scale;
 
@@ -200,30 +227,33 @@ impl Camera {
         image::Rgb([r as u8, g as u8, b as u8])
     }
 
-    pub fn get_ray(&self, x: u32, y: u32) -> Ray {
-        let offset = math::random_box(-0.5f64..0.5f64);
+    pub fn get_ray(&self, x: u32, y: u32, s_x: u32, s_y: u32, rng: &mut ThreadRng) -> Ray {
+        let (offset_x, offset_y) = sample_square_stratified(s_x, s_y, self.recip_sqrt_spp, rng);
         let pixel_sample = self.pixel00_loc
-            + (self.pixel_delta_u * (offset.x + x as f64))
-            + (self.pixel_delta_v * (offset.y + y as f64));
+            + (self.pixel_delta_u * (offset_x + x as f64))
+            + (self.pixel_delta_v * (offset_y + y as f64));
 
         let ray_origin = match self.defocus_angle <= 0.0 {
             true => self.center,
-            false => self.defocus_disk_sample(),
+            false => self.defocus_disk_sample(rng),
         };
         let ray_direction = pixel_sample - ray_origin;
 
-        let mut rng = rand::rng();
-        let t = rng.random_range(0.0f64..1.0f64);
-
-        Ray::new(ray_origin, ray_direction, t)
+        Ray::new(ray_origin, ray_direction, rng.random::<f64>())
     }
 
-    pub fn defocus_disk_sample(&self) -> Vector3<f64> {
-        let p = random_in_unit_disk();
+    pub fn defocus_disk_sample(&self, rng: &mut ThreadRng) -> Vector3<f64> {
+        let p = random_in_unit_disk(rng);
         self.center + (p.x * self.defocus_disk_u) + (p.y * self.defocus_disk_v)
     }
 
-    pub fn ray_color(&self, ray: &Ray, depth: u32, world: &Geometry) -> Vector3<f64> {
+    pub fn ray_color(
+        &self,
+        ray: &Ray,
+        depth: u32,
+        world: &Geometry,
+        rng: &mut ThreadRng,
+    ) -> Vector3<f64> {
         if depth == 0 {
             return Vector3::default();
         }
@@ -231,7 +261,7 @@ impl Camera {
         let mut hit_record = HitRecord::default();
         let interval = Interval::new(0.001, f64::INFINITY);
 
-        if !world.hit(ray, &interval, &mut hit_record) {
+        if !world.hit(ray, &interval, &mut hit_record, rng) {
             return self.background;
         }
 
@@ -244,13 +274,13 @@ impl Camera {
 
         if !hit_record
             .material
-            .scatter(ray, &hit_record, &mut attenuation, &mut scattered)
+            .scatter(ray, &hit_record, &mut attenuation, &mut scattered, rng)
         {
             return color_from_emission;
         }
 
         let color_from_scatter =
-            attenuation.component_mul(&self.ray_color(&scattered, depth - 1, world));
+            attenuation.component_mul(&self.ray_color(&scattered, depth - 1, world, rng));
 
         color_from_emission + color_from_scatter
     }
diff --git a/src/geometry/bvh.rs b/src/geometry/bvh.rs
index 1473ba9..9324e67 100644
--- a/src/geometry/bvh.rs
+++ b/src/geometry/bvh.rs
@@ -5,6 +5,7 @@ use crate::geometry::aabb::Aabb;
 use crate::geometry::empty::Empty;
 use crate::interval::Interval;
 use crate::ray::Ray;
+use rand::rngs::ThreadRng;
 
 #[derive(Debug, Clone)]
 pub struct BvhNode {
@@ -56,12 +57,18 @@ impl BvhNode {
 }
 
 impl Hittable for BvhNode {
-    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord) -> bool {
+    fn hit(
+        &self,
+        r: &Ray,
+        interval: &Interval,
+        record: &mut HitRecord,
+        rng: &mut ThreadRng,
+    ) -> bool {
         if !self.bbox.hit(r, interval) {
             return false;
         }
 
-        let left_hit = self.left.hit(r, interval, record);
+        let left_hit = self.left.hit(r, interval, record, rng);
         let new_interval = match left_hit {
             true => &Interval {
                 min: interval.min,
@@ -70,7 +77,7 @@ impl Hittable for BvhNode {
             false => interval,
         };
 
-        let right_hit = self.right.hit(r, new_interval, record);
+        let right_hit = self.right.hit(r, new_interval, record, rng);
 
         left_hit || right_hit
     }
diff --git a/src/geometry/cube.rs b/src/geometry/cube.rs
index 0c2c8ee..2d14504 100644
--- a/src/geometry/cube.rs
+++ b/src/geometry/cube.rs
@@ -8,6 +8,7 @@ use crate::interval::Interval;
 use crate::material::Material;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 
 #[derive(Debug, Clone)]
 pub struct Cube {
@@ -42,8 +43,14 @@ impl Cube {
 }
 
 impl Hittable for Cube {
-    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord) -> bool {
-        self.children.hit(r, interval, record)
+    fn hit(
+        &self,
+        r: &Ray,
+        interval: &Interval,
+        record: &mut HitRecord,
+        rng: &mut ThreadRng,
+    ) -> bool {
+        self.children.hit(r, interval, record, rng)
     }
 
     fn bounding_box(&self) -> Aabb {
diff --git a/src/geometry/empty.rs b/src/geometry/empty.rs
index 5736a49..d070587 100644
--- a/src/geometry/empty.rs
+++ b/src/geometry/empty.rs
@@ -4,6 +4,7 @@ use crate::geometry::Hittable;
 use crate::geometry::aabb::Aabb;
 use crate::interval::Interval;
 use crate::ray::Ray;
+use rand::rngs::ThreadRng;
 
 #[derive(Debug, Clone)]
 pub struct Empty {}
@@ -15,7 +16,7 @@ impl Empty {
 }
 
 impl Hittable for Empty {
-    fn hit(&self, _: &Ray, _: &Interval, _: &mut HitRecord) -> bool {
+    fn hit(&self, _: &Ray, _: &Interval, _: &mut HitRecord, _: &mut ThreadRng) -> bool {
         false
     }
 
diff --git a/src/geometry/mod.rs b/src/geometry/mod.rs
index 723c92b..5468bb5 100644
--- a/src/geometry/mod.rs
+++ b/src/geometry/mod.rs
@@ -24,9 +24,16 @@ use crate::material::lambertian::Lambertian;
 use crate::material::texture::SolidColor;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 
 pub trait Hittable {
-    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord) -> bool;
+    fn hit(
+        &self,
+        r: &Ray,
+        interval: &Interval,
+        record: &mut HitRecord,
+        rng: &mut ThreadRng,
+    ) -> bool;
     fn bounding_box(&self) -> Aabb;
 }
 
@@ -43,16 +50,22 @@ pub enum Geometry {
 }
 
 impl Hittable for Geometry {
-    fn hit(&self, ray: &Ray, interval: &Interval, record: &mut HitRecord) -> bool {
+    fn hit(
+        &self,
+        ray: &Ray,
+        interval: &Interval,
+        record: &mut HitRecord,
+        rng: &mut ThreadRng,
+    ) -> bool {
         match self {
-            Geometry::Empty(geomtry) => geomtry.hit(ray, interval, record),
-            Geometry::Quad(geomtry) => geomtry.hit(ray, interval, record),
-            Geometry::Sphere(geomtry) => geomtry.hit(ray, interval, record),
-            Geometry::BvhNode(geomtry) => geomtry.hit(ray, interval, record),
-            Geometry::Cube(geomtry) => geomtry.hit(ray, interval, record),
-            Geometry::Translate(geomtry) => geomtry.hit(ray, interval, record),
-            Geometry::Rotate(geomtry) => geomtry.hit(ray, interval, record),
-            Geometry::Volume(geomtry) => geomtry.hit(ray, interval, record),
+            Geometry::Empty(geomtry) => geomtry.hit(ray, interval, record, rng),
+            Geometry::Quad(geomtry) => geomtry.hit(ray, interval, record, rng),
+            Geometry::Sphere(geomtry) => geomtry.hit(ray, interval, record, rng),
+            Geometry::BvhNode(geomtry) => geomtry.hit(ray, interval, record, rng),
+            Geometry::Cube(geomtry) => geomtry.hit(ray, interval, record, rng),
+            Geometry::Translate(geomtry) => geomtry.hit(ray, interval, record, rng),
+            Geometry::Rotate(geomtry) => geomtry.hit(ray, interval, record, rng),
+            Geometry::Volume(geomtry) => geomtry.hit(ray, interval, record, rng),
         }
     }
 
diff --git a/src/geometry/quad.rs b/src/geometry/quad.rs
index 61b56f9..1655356 100644
--- a/src/geometry/quad.rs
+++ b/src/geometry/quad.rs
@@ -6,6 +6,7 @@ use crate::interval::Interval;
 use crate::material::Material;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 
 #[derive(Debug, Clone)]
 pub struct Quad {
@@ -63,7 +64,7 @@ fn is_interior(a: f64, b: f64, record: &mut HitRecord) -> bool {
 }
 
 impl Hittable for Quad {
-    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord) -> bool {
+    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord, _: &mut ThreadRng) -> bool {
         let denom = self.normal.dot(&r.direction);
         if denom.abs() < 1e-8 {
             return false;
diff --git a/src/geometry/rotate.rs b/src/geometry/rotate.rs
index ed08299..7b28093 100644
--- a/src/geometry/rotate.rs
+++ b/src/geometry/rotate.rs
@@ -7,6 +7,7 @@ use crate::interval::Interval;
 use crate::ray::Ray;
 use nalgebra::Rotation3;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 
 #[derive(Debug, Clone)]
 pub struct Rotate {
@@ -56,14 +57,20 @@ impl Rotate {
 }
 
 impl Hittable for Rotate {
-    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord) -> bool {
+    fn hit(
+        &self,
+        r: &Ray,
+        interval: &Interval,
+        record: &mut HitRecord,
+        rng: &mut ThreadRng,
+    ) -> bool {
         // Transform the ray from world space to object space.
         let rotated_origin = self.rotation.inverse() * r.origin;
         let rotated_direction = self.rotation.inverse() * r.direction;
         let rotated_ray = Ray::new(rotated_origin, rotated_direction, r.time);
 
         // Determine whether an intersection exists in object space (and if so, where).
-        if !self.geometry.hit(&rotated_ray, interval, record) {
+        if !self.geometry.hit(&rotated_ray, interval, record, rng) {
             return false;
         }
 
diff --git a/src/geometry/sphere.rs b/src/geometry/sphere.rs
index 9cb5988..8fe76e0 100644
--- a/src/geometry/sphere.rs
+++ b/src/geometry/sphere.rs
@@ -6,6 +6,7 @@ use crate::interval::Interval;
 use crate::material::Material;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 use std::f64::consts::PI;
 
 #[derive(Debug, Clone)]
@@ -55,7 +56,7 @@ impl Sphere {
 }
 
 impl Hittable for Sphere {
-    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord) -> bool {
+    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord, _: &mut ThreadRng) -> bool {
         let current_center = self.center.at(r.time);
         let oc = r.origin - current_center;
 
diff --git a/src/geometry/translate.rs b/src/geometry/translate.rs
index 2e236c5..0e35060 100644
--- a/src/geometry/translate.rs
+++ b/src/geometry/translate.rs
@@ -5,6 +5,7 @@ use crate::geometry::aabb::Aabb;
 use crate::interval::Interval;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 
 #[derive(Debug, Clone)]
 pub struct Translate {
@@ -29,10 +30,16 @@ impl Translate {
 }
 
 impl Hittable for Translate {
-    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord) -> bool {
+    fn hit(
+        &self,
+        r: &Ray,
+        interval: &Interval,
+        record: &mut HitRecord,
+        rng: &mut ThreadRng,
+    ) -> bool {
         let offset_ray = Ray::new(r.origin - self.offset, r.direction, r.time);
 
-        if !self.geometry.hit(&offset_ray, interval, record) {
+        if !self.geometry.hit(&offset_ray, interval, record, rng) {
             return false;
         }
 
diff --git a/src/geometry/volume.rs b/src/geometry/volume.rs
index 732714c..1dc5b14 100644
--- a/src/geometry/volume.rs
+++ b/src/geometry/volume.rs
@@ -9,6 +9,7 @@ use crate::material::texture::Texture;
 use crate::ray::Ray;
 use nalgebra::Vector3;
 use rand::prelude::*;
+use rand::rngs::ThreadRng;
 
 #[derive(Debug, Clone)]
 pub struct Volume {
@@ -34,11 +35,20 @@ impl Volume {
 }
 
 impl Hittable for Volume {
-    fn hit(&self, r: &Ray, interval: &Interval, record: &mut HitRecord) -> bool {
+    fn hit(
+        &self,
+        r: &Ray,
+        interval: &Interval,
+        record: &mut HitRecord,
+        rng: &mut ThreadRng,
+    ) -> bool {
         let mut record_a = HitRecord::default();
         let mut record_b = HitRecord::default();
 
-        if !self.boundry.hit(r, &Interval::universe(), &mut record_a) {
+        if !self
+            .boundry
+            .hit(r, &Interval::universe(), &mut record_a, rng)
+        {
             return false;
         }
 
@@ -46,6 +56,7 @@ impl Hittable for Volume {
             r,
             &Interval::new(record_a.t + 0.0001, f64::INFINITY),
             &mut record_b,
+            rng,
         ) {
             return false;
         }
@@ -68,9 +79,7 @@ impl Hittable for Volume {
 
         let ray_length = r.direction.norm();
         let distance_inside_boundary = (record_b.t - record_a.t) * ray_length;
-        let mut rng = rand::rng();
-        let random_double = rng.random_range(0.0f64..1.0f64);
-        let hit_distance = self.neg_inv_density * random_double.ln();
+        let hit_distance = self.neg_inv_density * rng.random::<f64>().ln();
 
         if hit_distance > distance_inside_boundary {
             return false;
diff --git a/src/material/dielectric.rs b/src/material/dielectric.rs
index bed9f46..72d44b5 100644
--- a/src/material/dielectric.rs
+++ b/src/material/dielectric.rs
@@ -7,6 +7,7 @@ use crate::math::refract;
 use crate::ray::Ray;
 use nalgebra::Vector3;
 use rand::prelude::*;
+use rand::rngs::ThreadRng;
 use std::fmt::Debug;
 
 #[derive(Debug, Clone)]
@@ -27,6 +28,7 @@ impl Surface for Dielectric {
         record: &HitRecord,
         attenuation: &mut Vector3<f64>,
         scattered: &mut Ray,
+        rng: &mut ThreadRng,
     ) -> bool {
         attenuation.copy_from(&Vector3::from_element(1.0));
         let r_index = match record.front_face {
@@ -42,13 +44,12 @@ impl Surface for Dielectric {
 
         scattered.origin = record.point;
         scattered.time = r_in.time;
-        let mut rng = rand::rng();
-        let random = rng.random_range(0.0f64..1.0f64);
 
-        scattered.direction = match cannot_refract || (reflectance(cos_theta, r_index) > random) {
-            true => reflect(&normalized_direction, &record.normal),
-            false => refract(&normalized_direction, &record.normal, r_index),
-        };
+        scattered.direction =
+            match cannot_refract || (reflectance(cos_theta, r_index) > rng.random::<f64>()) {
+                true => reflect(&normalized_direction, &record.normal),
+                false => refract(&normalized_direction, &record.normal, r_index),
+            };
 
         true
     }
diff --git a/src/material/isotropic.rs b/src/material/isotropic.rs
index 8208c44..e418312 100644
--- a/src/material/isotropic.rs
+++ b/src/material/isotropic.rs
@@ -6,6 +6,7 @@ use crate::material::texture::Texture;
 use crate::math;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 use std::fmt::Debug;
 
 #[derive(Debug, Clone)]
@@ -29,10 +30,11 @@ impl Surface for Isotropic {
         record: &HitRecord,
         attenuation: &mut Vector3<f64>,
         scattered: &mut Ray,
+        rng: &mut ThreadRng,
     ) -> bool {
         scattered.origin = record.point;
         scattered.time = r_in.time;
-        scattered.direction = math::random_normal();
+        scattered.direction = math::random_normal(rng);
         attenuation.copy_from(&self.texture.sample(record.u, record.v, record.point));
         true
     }
diff --git a/src/material/lambertian.rs b/src/material/lambertian.rs
index 6666f39..c0c479a 100644
--- a/src/material/lambertian.rs
+++ b/src/material/lambertian.rs
@@ -7,6 +7,7 @@ use crate::math;
 use crate::math::near_zero;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 use std::fmt::Debug;
 
 #[derive(Debug, Clone)]
@@ -27,8 +28,9 @@ impl Surface for Lambertian {
         record: &HitRecord,
         attenuation: &mut Vector3<f64>,
         scattered: &mut Ray,
+        rng: &mut ThreadRng,
     ) -> bool {
-        let mut scatter_direction = record.normal + math::random_normal();
+        let mut scatter_direction = record.normal + math::random_normal(rng);
         if near_zero(&scatter_direction) {
             scatter_direction = record.normal;
         }
diff --git a/src/material/light.rs b/src/material/light.rs
index b1f0153..cdd1964 100644
--- a/src/material/light.rs
+++ b/src/material/light.rs
@@ -5,6 +5,7 @@ use crate::material::texture::Sample;
 use crate::material::texture::Texture;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 use std::fmt::Debug;
 
 #[derive(Debug, Clone)]
@@ -19,7 +20,14 @@ impl Light {
 }
 
 impl Surface for Light {
-    fn scatter(&self, _: &Ray, _: &HitRecord, _: &mut Vector3<f64>, _: &mut Ray) -> bool {
+    fn scatter(
+        &self,
+        _: &Ray,
+        _: &HitRecord,
+        _: &mut Vector3<f64>,
+        _: &mut Ray,
+        _: &mut ThreadRng,
+    ) -> bool {
         false
     }
 
diff --git a/src/material/metal.rs b/src/material/metal.rs
index 889f80e..2aef936 100644
--- a/src/material/metal.rs
+++ b/src/material/metal.rs
@@ -5,6 +5,7 @@ use crate::math;
 use crate::math::reflect;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 use std::fmt::Debug;
 
 #[derive(Debug, Clone)]
@@ -26,9 +27,10 @@ impl Surface for Metal {
         record: &HitRecord,
         attenuation: &mut Vector3<f64>,
         scattered: &mut Ray,
+        rng: &mut ThreadRng,
     ) -> bool {
         let mut reflected = reflect(&r_in.direction, &record.normal);
-        reflected = reflected.normalize() + (math::random_normal() * self.roughness);
+        reflected = reflected.normalize() + (math::random_normal(rng) * self.roughness);
 
         scattered.origin = record.point;
         scattered.time = r_in.time;
diff --git a/src/material/mod.rs b/src/material/mod.rs
index e593d33..54ff210 100644
--- a/src/material/mod.rs
+++ b/src/material/mod.rs
@@ -13,6 +13,7 @@ use crate::material::light::Light;
 use crate::material::metal::Metal;
 use crate::ray::Ray;
 use nalgebra::Vector3;
+use rand::rngs::ThreadRng;
 use std::fmt::Debug;
 
 pub trait Surface {
@@ -22,6 +23,7 @@ pub trait Surface {
         record: &HitRecord,
         attenuation: &mut Vector3<f64>,
         scattered: &mut Ray,
+        rng: &mut ThreadRng,
     ) -> bool;
     fn emitted(&self, u: f64, v: f64, p: Vector3<f64>) -> Vector3<f64>;
 }
@@ -42,18 +44,23 @@ impl Surface for Material {
         record: &HitRecord,
         attenuation: &mut Vector3<f64>,
         scattered: &mut Ray,
+        rng: &mut ThreadRng,
     ) -> bool {
         match self {
-            Material::Metal(material) => material.scatter(ray_in, record, attenuation, scattered),
+            Material::Metal(material) => {
+                material.scatter(ray_in, record, attenuation, scattered, rng)
+            }
             Material::Dielectric(material) => {
-                material.scatter(ray_in, record, attenuation, scattered)
+                material.scatter(ray_in, record, attenuation, scattered, rng)
             }
             Material::Lambertian(material) => {
-                material.scatter(ray_in, record, attenuation, scattered)
+                material.scatter(ray_in, record, attenuation, scattered, rng)
+            }
+            Material::Light(material) => {
+                material.scatter(ray_in, record, attenuation, scattered, rng)
             }
-            Material::Light(material) => material.scatter(ray_in, record, attenuation, scattered),
             Material::Isotropic(material) => {
-                material.scatter(ray_in, record, attenuation, scattered)
+                material.scatter(ray_in, record, attenuation, scattered, rng)
             }
         }
     }
diff --git a/src/math/mod.rs b/src/math/mod.rs
index 2b80335..06cb18a 100644
--- a/src/math/mod.rs
+++ b/src/math/mod.rs
@@ -1,28 +1,17 @@
 use nalgebra::Vector3;
 use rand::prelude::*;
-use std::ops::Range;
 
-pub fn random(range: Range<f64>) -> Vector3<f64> {
-    let mut rng = rand::rng();
+pub fn random_vector(rng: &mut ThreadRng) -> Vector3<f64> {
     Vector3::new(
-        rng.random_range(range.clone()),
-        rng.random_range(range.clone()),
-        rng.random_range(range.clone()),
+        rng.random_range(-1.0f64..1.0f64),
+        rng.random_range(-1.0f64..1.0f64),
+        rng.random_range(-1.0f64..1.0f64),
     )
 }
 
-pub fn random_box(range: Range<f64>) -> Vector3<f64> {
-    let mut rng = rand::rng();
-    Vector3::new(
-        rng.random_range(range.clone()),
-        rng.random_range(range.clone()),
-        0.0,
-    )
-}
-
-pub fn random_normal() -> Vector3<f64> {
+pub fn random_normal(rng: &mut ThreadRng) -> Vector3<f64> {
     loop {
-        let position = random(-1.0f64..1.0f64);
+        let position = random_vector(rng);
         let lensq = position.norm_squared();
         if f64::EPSILON < lensq && lensq <= 1.0 {
             return position / lensq.sqrt();
diff --git a/src/noise/mod.rs b/src/noise/mod.rs
index 3cc6c88..d055eec 100644
--- a/src/noise/mod.rs
+++ b/src/noise/mod.rs
@@ -1,6 +1,7 @@
 use itertools::iproduct;
 use nalgebra::Vector3;
 use rand::prelude::*;
+use rand::rngs::ThreadRng;
 
 const PERLIN_POINT_COUNT: usize = 256;
 
@@ -14,13 +15,12 @@ pub struct Perlin {
 
 impl Default for Perlin {
     fn default() -> Self {
-        Perlin::new()
+        Perlin::new(&mut rand::rng())
     }
 }
 
 impl Perlin {
-    pub fn new() -> Self {
-        let mut rng = rand::rng();
+    pub fn new(rng: &mut ThreadRng) -> Self {
         let mut randvec = [Vector3::<f64>::default(); PERLIN_POINT_COUNT];
         for vec in randvec.iter_mut() {
             *vec = Vector3::new(
@@ -32,13 +32,13 @@ impl Perlin {
         }
 
         let mut perm_x = Box::new([0; PERLIN_POINT_COUNT]);
-        Perlin::perlin_generate_perm(&mut perm_x);
+        Perlin::perlin_generate_perm(&mut perm_x, rng);
 
         let mut perm_y = Box::new([0; PERLIN_POINT_COUNT]);
-        Perlin::perlin_generate_perm(&mut perm_y);
+        Perlin::perlin_generate_perm(&mut perm_y, rng);
 
         let mut perm_z = Box::new([0; PERLIN_POINT_COUNT]);
-        Perlin::perlin_generate_perm(&mut perm_z);
+        Perlin::perlin_generate_perm(&mut perm_z, rng);
 
         Perlin {
             randvec: Box::new(randvec),
@@ -69,13 +69,12 @@ impl Perlin {
         Perlin::perlin_interp(&c, u, v, w)
     }
 
-    pub fn perlin_generate_perm(p: &mut [usize; PERLIN_POINT_COUNT]) {
+    pub fn perlin_generate_perm(p: &mut [usize; PERLIN_POINT_COUNT], rng: &mut ThreadRng) {
         p.iter_mut().enumerate().for_each(|(i, p_i)| *p_i = i);
-        Perlin::permute(p, PERLIN_POINT_COUNT);
+        Perlin::permute(p, PERLIN_POINT_COUNT, rng);
     }
 
-    pub fn permute(p: &mut [usize; PERLIN_POINT_COUNT], n: usize) {
-        let mut rng = rand::rng();
+    pub fn permute(p: &mut [usize; PERLIN_POINT_COUNT], n: usize, rng: &mut ThreadRng) {
         (1..n).rev().for_each(|i| p.swap(i, rng.random_range(0..i)));
     }