diff --git a/include/box.hpp b/include/box.hpp
index feea59f..6f9766b 100644
--- a/include/box.hpp
+++ b/include/box.hpp
@@ -1,10 +1,12 @@
 #ifndef BOX_HPP
 #define BOX_HPP
 
+#include <array>
+
+#include "primitives.hpp"
 #include "rectangle.hpp"
-#include "rtweekend.hpp"
-#include "visit.hpp"
 
+namespace raytracer::scene {
 /// This class implements a axis aligned cuboid using 6 rectangles
 class box {
  public:
@@ -17,36 +19,12 @@ class box {
       , box_max { p1 }
       , material_type { mat_type } {
     /// Add six sides of the box based on box_min and box_max to sides
-    sides[0] = xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p1.z(), mat_type);
-    sides[1] = xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p0.z(), mat_type);
-    sides[2] = xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p1.y(), mat_type);
-    sides[3] = xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p0.y(), mat_type);
-    sides[4] = yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p1.x(), mat_type);
-    sides[5] = yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p0.x(), mat_type);
-  }
-
-  /// Compute ray interaction with the box
-  bool hit(auto& ctx, const ray& r, real_t min, real_t max, hit_record& rec,
-           material_t& hit_material_type) const {
-    hit_record temp_rec;
-    material_t temp_material_type;
-    auto hit_anything = false;
-    auto closest_so_far = max;
-    // Checking if the ray hits any of the sides
-    for (const auto& side : sides) {
-      if (dev_visit(
-              [&](auto&& arg) {
-                return arg.hit(ctx, r, min, closest_so_far, temp_rec,
-                               temp_material_type);
-              },
-              side)) {
-        hit_anything = true;
-        closest_so_far = temp_rec.t;
-        rec = temp_rec;
-        hit_material_type = temp_material_type;
-      }
-    }
-    return hit_anything;
+    sides[0] = xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p1.z(), material_type);
+    sides[1] = xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p0.z(), material_type);
+    sides[2] = xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p1.y(), material_type);
+    sides[3] = xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p0.y(), material_type);
+    sides[4] = yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p1.x(), material_type);
+    sides[5] = yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p0.x(), material_type);
   }
 
   point box_min;
@@ -54,5 +32,6 @@ class box {
   material_t material_type;
   std::array<rectangle_t, 6> sides;
 };
+} // namespace raytracer::scene
 
 #endif
diff --git a/include/build_parameters.hpp b/include/build_parameters.hpp
index 5ecc388..f04eb8d 100644
--- a/include/build_parameters.hpp
+++ b/include/build_parameters.hpp
@@ -1,6 +1,6 @@
 #ifndef BUILD_PARAMETERS_HPP
 #define BUILD_PARAMETERS_HPP
-namespace buildparams {
+namespace raytracer::buildparams {
 
 #ifdef USE_SINGLE_TASK
 constexpr bool use_single_task = true;
diff --git a/include/camera.hpp b/include/camera.hpp
index cc15f7f..7d4949f 100644
--- a/include/camera.hpp
+++ b/include/camera.hpp
@@ -3,9 +3,11 @@
 
 #include <cmath>
 
+#include "localrandom.hpp"
+#include "primitives.hpp"
 #include "ray.hpp"
-#include "rtweekend.hpp"
 
+namespace raytracer {
 /** Camera model
 
         This implements:
@@ -90,14 +92,14 @@ class camera {
           viewport local coordinates (s,t) based on viewport
           width, height and focus distance
   */
-  ray get_ray(real_t s, real_t t, LocalPseudoRNG& rng) const {
+  ray get_ray(real_t s, real_t t, random::LocalPseudoRNG& rng) const {
     vec rd = lens_radius * rng.in_unit_disk();
     vec offset = u * rd.x() + v * rd.y();
     return { origin + offset,
              lower_left_corner + s * horizontal + t * vertical - origin -
                  offset,
-             rng.float_t(time0, time1) };
+             rng.real(time0, time1) };
   }
 };
-
+} // namespace raytracer
 #endif
diff --git a/include/constant_medium.hpp b/include/constant_medium.hpp
index c4b26e8..37f64ec 100644
--- a/include/constant_medium.hpp
+++ b/include/constant_medium.hpp
@@ -5,9 +5,10 @@
 #include "material.hpp"
 #include "sphere.hpp"
 #include "texture.hpp"
-#include "visit.hpp"
 
-using hittableVolume_t = std::variant<sphere, box>;
+namespace raytracer::scene {
+
+using hittableVolume_t = std::variant<std::monostate, sphere, box>;
 
 /**
  * A ray going through the volume can either make it all the way through
@@ -25,60 +26,9 @@ class constant_medium {
       , neg_inv_density { -1 / d }
       , phase_function { isotropic_material { a } } {}
 
-  bool hit(auto& ctx, const ray& r, real_t min, real_t max, hit_record& rec,
-           material_t& hit_material_type) const {
-    auto& rng = ctx.rng;
-    hit_material_type = phase_function;
-    material_t temp_material_type;
-    hit_record rec1, rec2;
-    if (!dev_visit(
-            [&](auto&& arg) {
-              return arg.hit(ctx, r, -infinity, infinity, rec1,
-                             temp_material_type);
-            },
-            boundary)) {
-      return false;
-    }
-
-    if (!dev_visit(
-            [&](auto&& arg) {
-              return arg.hit(ctx, r, rec1.t + 0.0001f, infinity, rec2,
-                             temp_material_type);
-            },
-            boundary)) {
-      return false;
-    }
-
-    if (rec1.t < min)
-      rec1.t = min;
-    if (rec2.t > max)
-      rec2.t = max;
-    if (rec1.t >= rec2.t)
-      return false;
-    if (rec1.t < 0)
-      rec1.t = 0;
-
-    const auto ray_length = sycl::length(r.direction());
-    /// Distance between the two hitpoints affect of probability
-    /// of the ray hitting a smoke particle
-    const auto distance_inside_boundary = (rec2.t - rec1.t) * ray_length;
-    const auto hit_distance = neg_inv_density * sycl::log(rng.float_t());
-
-    /// With lower density, hit_distance has higher probabilty
-    /// of being greater than distance_inside_boundary
-    if (hit_distance > distance_inside_boundary)
-      return false;
-
-    rec.t = rec1.t + hit_distance / ray_length;
-    rec.p = r.at(rec.t);
-
-    rec.normal = vec { 1, 0, 0 }; // arbitrary
-    rec.front_face = true;        // also arbitrary
-    return true;
-  }
-
   hittableVolume_t boundary;
   real_t neg_inv_density;
   material_t phase_function;
 };
+} // namespace raytracer::scene
 #endif
diff --git a/include/hit_record.hpp b/include/hit_record.hpp
new file mode 100644
index 0000000..ce7338b
--- /dev/null
+++ b/include/hit_record.hpp
@@ -0,0 +1,25 @@
+#ifndef HIT_RECORD_HPP
+#define HIT_RECORD_HPP
+
+#include "primitives.hpp"
+
+namespace raytracer::visitor {
+class hit_record {
+ public:
+  real_t t;
+  point p;         // hit point
+  vec normal;      // normal at hit point
+  bool front_face; // to check if hit point is on the outer surface
+  /*local coordinates for rectangles
+  and mercator coordinates for spheres */
+  real_t u;
+  real_t v;
+
+  // To set if the hit point is on the front face
+  void set_face_normal(const ray& r, const vec& outward_normal) {
+    front_face = dot(r.direction(), outward_normal) < 0;
+    normal = front_face ? outward_normal : vec {} - outward_normal;
+  }
+};
+} // namespace raytracer::visitor
+#endif
\ No newline at end of file
diff --git a/include/hitable.hpp b/include/hitable.hpp
index d5eb79a..0443b11 100644
--- a/include/hitable.hpp
+++ b/include/hitable.hpp
@@ -1,26 +1,20 @@
 #ifndef HITTABLE_H
 #define HITTABLE_H
 
-#include "ray.hpp"
-#include "rtweekend.hpp"
-#include "vec.hpp"
+#include <variant>
 
-class hit_record {
- public:
-  float t;         //
-  point p;         // hit point
-  vec normal;      // normal at hit point
-  bool front_face; // to check if hit point is on the outer surface
-  /*local coordinates for rectangles
-  and mercator coordintes for spheres */
-  float u;
-  float v;
+#include "hit_record.hpp"
+#include "material.hpp"
 
-  // To set if the hit point is on the front face
-  void set_face_normal(const ray& r, const vec& outward_normal) {
-    front_face = dot(r.direction(), outward_normal) < 0;
-    normal = front_face ? outward_normal : vec {} - outward_normal;
-  }
-};
+#include "box.hpp"
+#include "constant_medium.hpp"
+#include "ray.hpp"
+#include "rectangle.hpp"
+#include "sphere.hpp"
+#include "triangle.hpp"
 
+namespace raytracer::scene {
+using hittable_t = std::variant<std::monostate, sphere, xy_rect, triangle, box,
+                                constant_medium>;
+}
 #endif
diff --git a/include/hitable_visitor.hpp b/include/hitable_visitor.hpp
new file mode 100644
index 0000000..ecaacf9
--- /dev/null
+++ b/include/hitable_visitor.hpp
@@ -0,0 +1,274 @@
+#ifndef HITABLE_VISITOR_HPP
+#define HITABLE_VISITOR_HPP
+
+#include "hitable.hpp"
+#include "sycl.hpp"
+#include "task_context.hpp"
+#include "visit.hpp"
+
+namespace raytracer::visitor {
+inline bool
+badouel_ray_triangle_intersec(const ray& r,
+                              raytracer::scene::_triangle_coord const& tri,
+                              real_t min, real_t max, hit_record& rec) {
+  // Get triangle edge vectors and plane normal
+  auto u = tri.v1 - tri.v0;
+  auto v = tri.v2 - tri.v0;
+  vec outward_normal = sycl::cross(u, v);
+
+  auto w0 = r.origin() - tri.v0;
+  auto a = -sycl::dot(outward_normal, w0);
+  auto b = sycl::dot(outward_normal, r.direction());
+
+  // ray is parallel to triangle plane
+  if (sycl::fabs(b) < 0.000001f)
+    return false;
+
+  // intersection point of ray with triangle
+  real_t length = a / b;
+  if (length < 0)
+    return false;
+  else if (length < min || length > max)
+    return false;
+
+  vec hit_pt = r.at(length);
+  auto uu = sycl::dot(u, u);
+  auto uv = sycl::dot(u, v);
+  auto vv = sycl::dot(v, v);
+  auto w = hit_pt - tri.v0;
+  auto wu = sycl::dot(w, u);
+  auto wv = sycl::dot(w, v);
+  auto D = uv * uv - uu * vv;
+
+  auto s = (uv * wv - vv * wu) / D;
+  auto t = (uv * wu - uu * wv) / D;
+  if (s < 0.0f || s > 1.0f || t < 0.0f || (s + t) > 1.0f)
+    return false;
+
+  rec.set_face_normal(r, outward_normal);
+  rec.t = length;
+  rec.p = hit_pt;
+  return true;
+};
+
+inline bool
+moller_trumbore_triangle_intersec(const ray& r,
+                                  raytracer::scene::_triangle_coord const& tri,
+                                  real_t min, real_t max, hit_record& rec) {
+  constexpr auto epsilon = 0.0000001f;
+
+  // Get triangle edge vectors and plane normal
+  auto edge1 = tri.v1 - tri.v0;
+  auto edge2 = tri.v2 - tri.v0;
+  vec h = sycl::cross(r.direction(), edge2);
+  auto a = sycl::dot(edge1, h);
+  auto a_abs = sycl::fabs(a);
+
+  if (a_abs < epsilon)
+    return false; // This ray is parallel to this triangle
+
+  auto a_pos = a > 0.f;
+
+  auto s = r.origin() - tri.v0;
+  auto u = sycl::dot(s, h);
+  auto u_pos = u > 0.f;
+
+  if ((u_pos xor a_pos) || sycl::fabs(u) > a_abs)
+    return false;
+
+  auto q = sycl::cross(s, edge1);
+  auto v = sycl::dot(r.direction(), q);
+  auto v_pos = v > 0.f;
+  if ((v_pos xor a_pos) || (sycl::fabs(u + v) > a_abs))
+    return false;
+
+  auto length = sycl::dot(edge2, q) / a;
+
+  if (length < min || length > max)
+    return false;
+
+  vec hit_pt = r.at(length);
+
+  rec.set_face_normal(r, sycl::cross(edge1, edge2));
+  rec.t = length;
+  rec.p = hit_pt;
+  return true;
+};
+
+template <auto IntersectionStrategy = moller_trumbore_triangle_intersec>
+struct _hitable_hit {
+ private:
+  using material_t = scene::material_t;
+  task_context& ctx;
+  const ray& r;
+  real_t min;
+  real_t max;
+  hit_record& rec;
+  material_t& hit_material_type;
+
+ private:
+ public:
+  _hitable_hit(task_context& ctx, const ray& r, real_t min, real_t max,
+               hit_record& rec, material_t& _hmt)
+      : ctx { ctx }
+      , r { r }
+      , min { min }
+      , max { max }
+      , rec { rec }
+      , hit_material_type { _hmt } {}
+
+  bool operator()(scene::box const& b) const {
+    hit_record temp_rec;
+    scene::material_t temp_material_type;
+    auto hit_anything = false;
+    auto closest_so_far = max;
+    // Checking if the ray hits any of the sides
+    for (const auto& side : b.sides) {
+      if (dev_visit(_hitable_hit(ctx, r, min, closest_so_far, temp_rec,
+                                 temp_material_type),
+                    side)) {
+        hit_anything = true;
+        closest_so_far = temp_rec.t;
+        rec = temp_rec;
+        hit_material_type = temp_material_type;
+      }
+    }
+    return hit_anything;
+  }
+
+  bool operator()(scene::constant_medium const& cm) const {
+    auto& rng = ctx.rng;
+    hit_material_type = cm.phase_function;
+    material_t temp_material_type;
+    hit_record rec1, rec2;
+    if (!dev_visit(
+            _hitable_hit(ctx, r, -infinity, infinity, rec1, temp_material_type),
+            cm.boundary)) {
+      return false;
+    }
+
+    if (!dev_visit(_hitable_hit(ctx, r, rec1.t + 0.0001f, infinity, rec2,
+                                temp_material_type),
+                   cm.boundary)) {
+      return false;
+    }
+
+    if (rec1.t < min)
+      rec1.t = min;
+    if (rec2.t > max)
+      rec2.t = max;
+    if (rec1.t >= rec2.t)
+      return false;
+    if (rec1.t < 0)
+      rec1.t = 0;
+
+    const auto ray_length = sycl::length(r.direction());
+    /// Distance between the two hitpoints affect of probability
+    /// of the ray hitting a smoke particle
+    const auto distance_inside_boundary = (rec2.t - rec1.t) * ray_length;
+    const auto hit_distance = cm.neg_inv_density * sycl::log(rng.real());
+
+    /// With lower density, hit_distance has higher probabilty
+    /// of being greater than distance_inside_boundary
+    if (hit_distance > distance_inside_boundary)
+      return false;
+
+    rec.t = rec1.t + hit_distance / ray_length;
+    rec.p = r.at(rec.t);
+
+    rec.normal = vec { 1, 0, 0 }; // arbitrary
+    rec.front_face = true;        // also arbitrary
+    return true;
+  }
+
+  template <bool x_plane, bool y_plane, bool z_plane>
+  bool operator()(scene::rect<x_plane, y_plane, z_plane> const& pl) {
+    hit_material_type = pl.material_type;
+    auto origin =
+        scene::rect<x_plane, y_plane, z_plane>::plan_basis(r.origin());
+    auto direction =
+        scene::rect<x_plane, y_plane, z_plane>::plan_basis(r.direction());
+
+    auto t = (pl.k - origin.z()) / direction.z();
+    if (t < min || t > max)
+      return false;
+    auto a = origin.x() + t * direction.x();
+    auto b = origin.y() + t * direction.y();
+    if (a < pl.a0 || a > pl.a1 || b < pl.b0 || b > pl.b1)
+      return false;
+    rec.u = (a - pl.a0) / (pl.a1 - pl.a0);
+    rec.v = (b - pl.b0) / (pl.b1 - pl.b0);
+    rec.t = t;
+    rec.p = r.at(rec.t);
+    vec outward_normal = vec(0, 0, 1);
+    rec.set_face_normal(r, outward_normal);
+    return true;
+  }
+
+  /// Compute ray interaction with sphere
+  bool operator()(scene::sphere const& s) {
+    hit_material_type = s.material_type;
+
+    auto actual_center = s.center(r.time());
+
+    /*(P(t)-C).(P(t)-C)=r^2
+    in the above sphere equation P(t) is the point on sphere hit by the ray
+    (A+tb−C)⋅(A+tb−C)=r^2
+    (t^2)b⋅b + 2tb⋅(A−C) + (A−C)⋅(A−C)−r^2 = 0
+    There can 0 or 1 or 2 real roots*/
+    vec oc = r.origin() - actual_center; // oc = A-C
+    auto a = sycl::dot(r.direction(), r.direction());
+    auto b = sycl::dot(oc, r.direction());
+    auto c = sycl::dot(oc, oc) - s.radius * s.radius;
+    auto discriminant = b * b - a * c;
+    // Real roots if discriminant is positive
+    if (discriminant > 0) {
+      // First root
+      auto temp = (-b - sycl::sqrt(discriminant)) / a;
+      if (temp < max && temp > min) {
+        rec.t = temp;
+        // Ray hits the sphere at p
+        rec.p = r.at(rec.t);
+        vec outward_normal = (rec.p - actual_center) / s.radius;
+        // To set if hit point is on the front face and the outward normal in
+        // rec
+        rec.set_face_normal(r, outward_normal);
+        /* Update u and v values in the hit record. Normal of a
+        point is calculated as above. This vector is used to also
+        used to get the mercator coordinates of the hitpoint.*/
+        std::tie(rec.u, rec.v) = mercator_coordinates(rec.normal);
+        return true;
+      }
+      // Second root
+      temp = (-b + sycl::sqrt(discriminant)) / a;
+      if (temp < max && temp > min) {
+        rec.t = temp;
+        // Ray hits the sphere at p
+        rec.p = r.at(rec.t);
+        vec outward_normal = (rec.p - actual_center) / s.radius;
+        rec.set_face_normal(r, outward_normal);
+        // Update u and v values in the hit record
+        std::tie(rec.u, rec.v) = mercator_coordinates(rec.normal);
+        return true;
+      }
+    }
+    // No real roots
+    return false;
+  }
+
+  bool operator()(scene::triangle const & tri) {
+    hit_material_type = tri.material_type;
+    return IntersectionStrategy(r, tri, min, max, rec);
+  }
+
+  bool operator()(std::monostate) {
+    assert(fase && "unreachable");
+    return false;
+  }
+};
+
+using hitable_hit = _hitable_hit<>;
+
+} // namespace raytracer::visitor
+
+#endif
\ No newline at end of file
diff --git a/include/localrandom.hpp b/include/localrandom.hpp
new file mode 100644
index 0000000..651c165
--- /dev/null
+++ b/include/localrandom.hpp
@@ -0,0 +1,71 @@
+#ifndef LOCALRANDOM_HPP
+#define LOCALRANDOM_HPP
+
+#include <xorshift.hpp>
+
+#include "primitives.hpp"
+
+namespace raytracer::random {
+class LocalPseudoRNG {
+ public:
+  inline LocalPseudoRNG(std::uint32_t init_state = xorshift<>::initial_state)
+      : generator { init_state } {}
+
+  // Returns a random float in 0.f 1.
+  inline real_t real() {
+    constexpr real_t scale = 1.f / (uint64_t { 1 } << 32);
+    return generator() * scale;
+  }
+
+  // Returns a random float in min, max
+  inline real_t real(real_t min, real_t max) {
+    // TODO use FMA ?
+    return min + (max - min) * real();
+  }
+
+  // Returns a random vector with coordinates in 0.f 1.
+  inline vec vec_t() { return { real(), real(), real() }; }
+
+  // Returns a random vec with coordinates in min, max
+  inline vec vec_t(real_t min, real_t max) {
+    auto scale = max - min;
+    return vec_t() * scale + min;
+  }
+
+  // Returns a random unit vector
+  inline vec unit_vec() {
+    auto x = real(-1.f, 1.f);
+    auto maxy = sycl::sqrt(1 - x * x);
+    auto y = real(-maxy, maxy);
+    auto absz = sycl::sqrt(maxy * maxy - y * y);
+    auto z = (real() > 0.5) ? absz : -absz;
+    return vec(x, y, z);
+  }
+
+  // Returns a random vector in the unit ball of usual norm
+  inline vec in_unit_ball() {
+    // Polar coordinates r, theta, phi
+    auto r = real();
+    auto theta = real(0, 2 * pi);
+    auto phi = real(0, pi);
+
+    auto plan_seed = r * sycl::sin(phi);
+    auto z = r * sycl::cos(phi);
+
+    return { plan_seed * sycl::cos(theta), plan_seed * sycl::sin(theta), z };
+  }
+
+  // Return a random vector in the unit disk of usual norm in plane x, y
+  inline vec in_unit_disk() {
+    auto x = real(-1.f, 1.f);
+    auto maxy = sycl::sqrt(1 - x * x);
+    auto y = real(-maxy, maxy);
+    return { x, y, 0.f };
+  }
+
+ private:
+  xorshift<> generator;
+};
+} // namespace raytracer::random
+
+#endif
\ No newline at end of file
diff --git a/include/material.hpp b/include/material.hpp
index d11bffb..0e640bf 100644
--- a/include/material.hpp
+++ b/include/material.hpp
@@ -1,55 +1,29 @@
 #ifndef RT_SYCL_MATERIAL_HPP
 #define RT_SYCL_MATERIAL_HPP
 
-#include <iostream>
+#include <variant>
 
-#include "hitable.hpp"
+#include "primitives.hpp"
 #include "texture.hpp"
-#include "vec.hpp"
-#include "visit.hpp"
 
+namespace raytracer::scene {
 struct lambertian_material {
   lambertian_material() = default;
   lambertian_material(const color& a)
       : albedo { solid_texture { a } } {}
   lambertian_material(const texture_t& a)
       : albedo { a } {}
-
-  bool scatter(auto& ctx, const ray& r_in, const hit_record& rec,
-               color& attenuation, ray& scattered) const {
-    auto& rng = ctx.rng;
-    vec scatter_direction = rec.normal + rng.unit_vec();
-    scattered = ray(rec.p, scatter_direction, r_in.time());
-    // Attenuation of the ray hitting the object is modified based on the color
-    // at hit point
-    attenuation *=
-        dev_visit([&](auto&& arg) { return arg.value(ctx, rec); }, albedo);
-    return true;
-  }
-  color emitted(auto&, const hit_record& rec) { return color(0, 0, 0); }
   texture_t albedo;
 };
 
 struct metal_material {
   metal_material() = default;
-  metal_material(const color& a, float f)
+  metal_material(const color& a, real_t f)
       : albedo { a }
       , fuzz { std::clamp(f, 0.0f, 1.0f) } {}
 
-  bool scatter(auto& ctx, const ray& r_in, const hit_record& rec,
-               color& attenuation, ray& scattered) const {
-    auto& rng = ctx.rng;
-    vec reflected = reflect(unit_vector(r_in.direction()), rec.normal);
-    scattered = ray(rec.p, reflected + fuzz * rng.in_unit_ball(), r_in.time());
-    // Attenuation of the ray hitting the object is modified based on the color
-    // at hit point
-    attenuation *= albedo;
-    return (dot(scattered.direction(), rec.normal) > 0);
-  }
-
-  color emitted(auto&, const hit_record& rec) { return color(0, 0, 0); }
   color albedo;
-  float fuzz;
+  real_t fuzz;
 };
 
 struct dielectric_material {
@@ -65,29 +39,6 @@ struct dielectric_material {
     return r0 + (1 - r0) * sycl::pow((1 - cosine), 5.0f);
   }
 
-  bool scatter(auto& ctx, const ray& r_in, const hit_record& rec,
-               color& attenuation, ray& scattered) const {
-    // Attenuation of the ray hitting the object is modified based on the color
-    // at hit point
-    auto& rng = ctx.rng;
-    attenuation *= albedo;
-    float refraction_ratio = rec.front_face ? (1.0f / ref_idx) : ref_idx;
-    vec unit_direction = unit_vector(r_in.direction());
-    float cos_theta = sycl::fmin(-sycl::dot(unit_direction, rec.normal), 1.0f);
-    float sin_theta = sycl::sqrt(1.0f - cos_theta * cos_theta);
-    bool cannot_refract = refraction_ratio * sin_theta > 1.0f;
-    vec direction;
-    if (cannot_refract ||
-        reflectance(cos_theta, refraction_ratio) > rng.float_t())
-      direction = reflect(unit_direction, rec.normal);
-    else
-      direction = refract(unit_direction, rec.normal, refraction_ratio);
-
-    scattered = ray(rec.p, direction, r_in.time());
-    return true;
-  }
-
-  color emitted(auto&, const hit_record& rec) { return color(0, 0, 0); }
   // Refractive index of the glass
   real_t ref_idx;
   // Color of the glass
@@ -101,12 +52,6 @@ struct lightsource_material {
   lightsource_material(const color& a)
       : emit { solid_texture { a } } {}
 
-  template <typename... T> bool scatter(T&...) const { return false; }
-
-  color emitted(auto& ctx, const hit_record& rec) {
-    return dev_visit([&](auto&& arg) { return arg.value(ctx, rec); }, emit);
-  }
-
   texture_t emit;
 };
 
@@ -115,23 +60,12 @@ struct isotropic_material {
       : albedo { solid_texture { a } } {}
   isotropic_material(texture_t& a)
       : albedo { a } {}
-
-  bool scatter(auto& ctx, const ray& r_in, const hit_record& rec,
-               color& attenuation, ray& scattered) const {
-    auto& rng = ctx.rng;
-    scattered = ray(rec.p, rng.in_unit_ball(), r_in.time());
-    attenuation *=
-        dev_visit([&](auto&& arg) { return arg.value(ctx, rec); }, albedo);
-    return true;
-  }
-
-  color emitted(auto&, const hit_record& rec) { return color(0, 0, 0); }
-
   texture_t albedo;
 };
 
 using material_t =
-    std::variant<lambertian_material, metal_material, dielectric_material,
-                 lightsource_material, isotropic_material>;
+    std::variant<std::monostate, lambertian_material, metal_material,
+                 dielectric_material, lightsource_material, isotropic_material>;
+} // namespace raytracer::scene
 
 #endif
diff --git a/include/material_visitor.hpp b/include/material_visitor.hpp
new file mode 100644
index 0000000..8eca982
--- /dev/null
+++ b/include/material_visitor.hpp
@@ -0,0 +1,105 @@
+#ifndef MATERIAL_VISITOR_HPP
+#define MATERIAL_VISITOR_HPP
+
+#include "material.hpp"
+#include "primitives.hpp"
+#include "texture_visitor.hpp"
+#include "visit.hpp"
+
+namespace raytracer::visitor {
+struct material_emitted {
+ private:
+  task_context& ctx;
+  const hit_record& rec;
+
+ public:
+  material_emitted(task_context& ctx, hit_record& rec)
+      : ctx { ctx }
+      , rec { rec } {}
+
+  template <typename M> color operator()(M&&) {
+    return { 0.f, 0.f, 0.f };
+  }
+
+  color operator()(scene::lightsource_material& mat) {
+      return dev_visit(texture_value(ctx, rec), mat.emit);
+  }
+
+  color operator()(std::monostate) {
+    assert(false && "unreachable");
+    return { 0.f, 0.f, 0.f };
+  }
+};
+
+struct material_scatter {
+ private:
+  task_context& ctx;
+  const ray& r_in;
+  const hit_record& rec;
+  color& attenuation;
+  ray& scattered;
+
+ public:
+  material_scatter(auto& ctx, const ray& r_in, const hit_record& rec,
+                           color& attenuation, ray& scattered)
+      : ctx { ctx }
+      , r_in { r_in }
+      , rec { rec }
+      , attenuation { attenuation }
+      , scattered { scattered } {}
+
+  bool operator()(scene::lambertian_material& mat) const {
+    vec scatter_direction = rec.normal + ctx.rng.unit_vec();
+    scattered = ray(rec.p, scatter_direction, r_in.time());
+    // Attenuation of the ray hitting the object is modified based on the color
+    // at hit point
+    attenuation *= dev_visit(texture_value(ctx, rec), mat.albedo);
+    return true;
+  }
+
+  bool operator()(scene::metal_material& mat) {
+    vec reflected = reflect(unit_vector(r_in.direction()), rec.normal);
+    scattered =
+        ray(rec.p, reflected + mat.fuzz * ctx.rng.in_unit_ball(), r_in.time());
+    // Attenuation of the ray hitting the object is modified based on the color
+    // at hit point
+    attenuation *= mat.albedo;
+    return (dot(scattered.direction(), rec.normal) > 0);
+  }
+
+  bool operator()(scene::dielectric_material& mat) {
+    // Attenuation of the ray hitting the object is modified based on the color
+    // at hit point
+    attenuation *= mat.albedo;
+    real_t refraction_ratio =
+        rec.front_face ? (1.0f / mat.ref_idx) : mat.ref_idx;
+    vec unit_direction = unit_vector(r_in.direction());
+    real_t cos_theta = sycl::fmin(-sycl::dot(unit_direction, rec.normal), 1.0f);
+    real_t sin_theta = sycl::sqrt(1.0f - cos_theta * cos_theta);
+    bool cannot_refract = refraction_ratio * sin_theta > 1.0f;
+    vec direction;
+    if (cannot_refract ||
+        mat.reflectance(cos_theta, refraction_ratio) > ctx.rng.real())
+      direction = reflect(unit_direction, rec.normal);
+    else
+      direction = refract(unit_direction, rec.normal, refraction_ratio);
+
+    scattered = ray(rec.p, direction, r_in.time());
+    return true;
+  }
+
+  bool operator()(scene::lightsource_material&) {return false;}
+
+  bool operator()(scene::isotropic_material& mat) {
+    scattered = ray(rec.p, ctx.rng.in_unit_ball(), r_in.time());
+    attenuation *= dev_visit(texture_value(ctx, rec), mat.albedo);
+    return true;
+  }
+
+  bool operator()(std::monostate) {
+    assert(false && "unreachable");
+    return false;
+  }
+};
+} // namespace raytracer::visitor
+#endif
\ No newline at end of file
diff --git a/include/primitives.hpp b/include/primitives.hpp
new file mode 100644
index 0000000..cd78b53
--- /dev/null
+++ b/include/primitives.hpp
@@ -0,0 +1,65 @@
+#ifndef PRIMITIVES_HPP
+#define PRIMITIVES_HPP
+
+#include "sycl.hpp"
+
+namespace raytracer {
+using real_t = float;
+
+// Constants
+
+constexpr real_t infinity = std::numeric_limits<real_t>::infinity();
+constexpr real_t pi = 3.1415926535897932385f;
+
+// type aliases for float3 - vec, point and color
+using real_vec = sycl::float3;
+
+using point = real_vec;
+using color = real_vec;
+using vec = real_vec;
+
+// Utility Functions
+inline real_t degrees_to_radians(real_t degrees) { return degrees * pi / 180.0f; }
+
+// vec Utility Functions
+inline real_t length_squared(const vec& v) {
+  return sycl::fma(v.x(), v.x(), sycl::fma(v.y(), v.y(), v.z() * v.z()));
+}
+
+// Missing operator from the SYCL specification for now
+vec operator-(const vec& u) { return vec(-u.x(), -u.y(), -u.z()); }
+
+// Compute a unit vector from a non-null vector
+inline vec unit_vector(const vec& v) { return v / sycl::length(v); }
+
+// Compute reflected ray's direction
+vec reflect(const vec& v, const vec& n) { return v - 2 * sycl::dot(v, n) * n; }
+
+// Computes refracted ray's direction based on refractive index
+vec refract(const vec& uv, const vec& n, real_t etai_over_etat) {
+  auto cos_theta = sycl::fmin(-sycl::dot(uv, n), 1.0f);
+  vec r_out_perp = etai_over_etat * (uv + cos_theta * n);
+  vec r_out_parallel =
+      -sycl::sqrt(sycl::fabs(1.0f - length_squared(r_out_perp))) * n;
+  return r_out_perp + r_out_parallel;
+}
+
+/* Computes normalised values of theta and phi. The input vector p
+corresponds to a vector passing through the centre of the a sphere
+and the hipoint on the surface of the sphere */
+std::pair<real_t, real_t> mercator_coordinates(const vec& p) {
+  // phi is the angle around the axis
+  auto phi = sycl::atan2(p.z(), p.x());
+  // theta is the angle down from the pole
+  auto theta = sycl::asin(p.y());
+  // theta and phi together constitute the spherical coordinates
+  // phi is between -pi and pi , u is between 0 and 1
+  auto u = 1 - (phi + pi) / (2 * pi);
+  // theta is between -pi/2 and pi/2 , v is between 0 and 1
+  auto v = (theta + pi / 2) / pi;
+  return { u, v };
+}
+
+} // namespace raytracer
+
+#endif
\ No newline at end of file
diff --git a/include/ray.hpp b/include/ray.hpp
index 902e615..d388178 100644
--- a/include/ray.hpp
+++ b/include/ray.hpp
@@ -1,8 +1,9 @@
 #ifndef RT_SYCL_RAY_HPP
 #define RT_SYCL_RAY_HPP
 
-#include "vec.hpp"
+#include "primitives.hpp"
 
+namespace raytracer {
 class ray {
  public:
   ray() = default;
@@ -18,7 +19,7 @@ class ray {
 
   // returns point along the ray at distance t from ray's origin
   // the ray P(t) = Origin + t*direction
-  point at(float t) const { return orig + t * dir; }
+  point at(real_t t) const { return orig + t * dir; }
 
  public:
   // To store the origin and direction of the ray
@@ -26,5 +27,5 @@ class ray {
   vec dir;
   real_t tm;
 };
-
+} // namespace raytracer
 #endif
diff --git a/include/rectangle.hpp b/include/rectangle.hpp
index 7ac9a08..4fa130f 100644
--- a/include/rectangle.hpp
+++ b/include/rectangle.hpp
@@ -2,131 +2,76 @@
 #define RECT_HPP
 
 #include "material.hpp"
-#include "ray.hpp"
-#include "rtweekend.hpp"
-#include "texture.hpp"
-#include "vec.hpp"
+#include "primitives.hpp"
 
 /** The Following classes implement:
 
         -
    https://raytracing.github.io/books/RayTracingTheNextWeek.html#rectanglesandlights/creatingrectangleobjectsa
 */
+namespace raytracer::scene {
 
-class xy_rect {
- public:
-  xy_rect() = default;
-
-  /// x0 <= x1 and y0 <= y1
-  xy_rect(real_t _x0, real_t _x1, real_t _y0, real_t _y1, real_t _k,
-          const material_t& mat_type)
-      : x0 { _x0 }
-      , x1 { _x1 }
-      , y0 { _y0 }
-      , y1 { _y1 }
-      , k { _k }
-      , material_type { mat_type } {}
+template <bool on_x_plane, bool on_y_plane, bool on_z_plane> class rect {
+  static_assert((on_x_plane && !on_y_plane && !on_z_plane) ||
+                    (!on_x_plane && on_y_plane && !on_z_plane) ||
+                    (!on_x_plane && !on_y_plane && on_z_plane),
+                "Only one parameter should be set to true");
 
-  /// Compute ray interaction with rectangle
-  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
-           material_t& hit_material_type) const {
-    hit_material_type = material_type;
-
-    auto t = (k - r.origin().z()) / r.direction().z();
-    if (t < min || t > max)
-      return false;
-    auto x = r.origin().x() + t * r.direction().x();
-    auto y = r.origin().y() + t * r.direction().y();
-    if (x < x0 || x > x1 || y < y0 || y > y1)
-      return false;
-    rec.u = (x - x0) / (x1 - x0);
-    rec.v = (y - y0) / (y1 - y0);
-    rec.t = t;
-    rec.p = r.at(rec.t);
-    vec outward_normal = vec(0, 0, 1);
-    rec.set_face_normal(r, outward_normal);
-    return true;
-  }
-  real_t x0, x1, y0, y1, k;
-  material_t material_type;
-};
-
-class xz_rect {
  public:
-  xz_rect() = default;
-
-  /// x0 <= x1 and z0 <= z1
-  xz_rect(real_t _x0, real_t _x1, real_t _z0, real_t _z1, real_t _k,
-          const material_t& mat_type)
-      : x0 { _x0 }
-      , x1 { _x1 }
-      , z0 { _z0 }
-      , z1 { _z1 }
+  real_t a0, a1, b0, b1, k;
+  material_t material_type;
+  rect(real_t _a0, real_t _a1, real_t _b0, real_t _b1, real_t _k,
+       const material_t& _mat_type)
+      : a0 { _a0 }
+      , a1 { _a1 }
+      , b0 { _b0 }
+      , b1 { _b1 }
       , k { _k }
-      , material_type { mat_type } {}
-
-  /// Compute ray interaction with rectangle
-  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
-           material_t& hit_material_type) const {
-    hit_material_type = material_type;
-
-    auto t = (k - r.origin().y()) / r.direction().y();
-    if (t < min || t > max)
-      return false;
-    auto x = r.origin().x() + t * r.direction().x();
-    auto z = r.origin().z() + t * r.direction().z();
-    if (x < x0 || x > x1 || z < z0 || z > z1)
-      return false;
-    rec.u = (x - x0) / (x1 - x0);
-    rec.v = (z - z0) / (z1 - z0);
-    rec.t = t;
-    rec.p = r.at(rec.t);
-    vec outward_normal = vec(0, 1, 0);
-    rec.set_face_normal(r, outward_normal);
-    return true;
+      , material_type { _mat_type } {}
+    
+  rect() = default;
+
+  static inline real_t first_planar_coord(vec const& values) {
+    if constexpr (on_x_plane) {
+      return values.y();
+    } else if constexpr (on_y_plane) {
+      return values.z();
+    } else {
+      return values.x();
+    }
   }
-  real_t x0, x1, z0, z1, k;
-  material_t material_type;
-};
-
-class yz_rect {
- public:
-  yz_rect() = default;
 
-  /// y0 <= y1 and z0 <= z1
-  yz_rect(real_t _y0, real_t _y1, real_t _z0, real_t _z1, real_t _k,
-          const material_t& mat_type)
-      : y0 { _y0 }
-      , y1 { _y1 }
-      , z0 { _z0 }
-      , z1 { _z1 }
-      , k { _k }
-      , material_type { mat_type } {}
+  static inline real_t second_planar_coord(vec const& values) {
+    if constexpr (on_x_plane) {
+      return values.z();
+    } else if constexpr (on_y_plane) {
+      return values.x();
+    } else {
+      return values.y();
+    }
+  }
 
-  /// Compute ray interaction with rectangle
-  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
-           material_t& hit_material_type) const {
-    hit_material_type = material_type;
+  static inline real_t normal_coord(vec const& values) {
+    if constexpr (on_x_plane) {
+      return values.x();
+    } else if constexpr (on_y_plane) {
+      return values.y();
+    } else {
+      return values.z();
+    }
+  }
 
-    auto t = (k - r.origin().x()) / r.direction().x();
-    if (t < min || t > max)
-      return false;
-    auto y = r.origin().y() + t * r.direction().y();
-    auto z = r.origin().z() + t * r.direction().z();
-    if (y < y0 || y > y1 || z < z0 || z > z1)
-      return false;
-    rec.u = (y - y0) / (y1 - y0);
-    rec.v = (z - z0) / (z1 - z0);
-    rec.t = t;
-    rec.p = r.at(rec.t);
-    vec outward_normal = vec(1, 0, 0);
-    rec.set_face_normal(r, outward_normal);
-    return true;
+  static inline vec plan_basis(vec const& values) {
+    return { first_planar_coord(values), second_planar_coord(values),
+             normal_coord(values) };
   }
-  real_t y0, y1, z0, z1, k;
-  material_t material_type;
 };
 
+using xy_rect = rect<false, false, true>;
+using xz_rect = rect<false, true, false>;
+using yz_rect = rect<true, false, false>;
+
 using rectangle_t = std::variant<xy_rect, xz_rect, yz_rect>;
+} // namespace raytracer::scene
 
 #endif
diff --git a/include/render.hpp b/include/render.hpp
index 565f037..5a7645f 100644
--- a/include/render.hpp
+++ b/include/render.hpp
@@ -3,44 +3,36 @@
 #include <variant>
 #include <vector>
 
-#include "box.hpp"
 #include "build_parameters.hpp"
 #include "camera.hpp"
-#include "constant_medium.hpp"
 #include "hitable.hpp"
+#include "hitable_visitor.hpp"
 #include "material.hpp"
+#include "material_visitor.hpp"
+#include "primitives.hpp"
 #include "ray.hpp"
-#include "rectangle.hpp"
-#include "rtweekend.hpp"
-#include "sphere.hpp"
 #include "sycl.hpp"
+#include "task_context.hpp"
 #include "texture.hpp"
-#include "triangle.hpp"
-#include "vec.hpp"
 #include "visit.hpp"
 
-using hittable_t =
-    std::variant<sphere, xy_rect, triangle, box, constant_medium>;
-
+namespace raytracer {
 template <int width, int height, int samples, int depth>
 inline auto render_pixel(auto& ctx, int x_coord, int y_coord, camera const& cam,
                          auto& hittable_acc, auto fb_acc) {
   auto& rng = ctx.rng;
   auto get_color = [&](const ray& r) {
-    auto hit_world = [&](const ray& r, hit_record& rec,
-                         material_t& material_type) {
-      hit_record temp_rec;
-      material_t temp_material_type;
+    auto hit_world = [&](const ray& r, visitor::hit_record& rec,
+                         scene::material_t& material_type) {
+      visitor::hit_record temp_rec;
+      scene::material_t temp_material_type;
       auto hit_anything = false;
       auto closest_so_far = infinity;
       // Checking if the ray hits any of the spheres
       for (auto i = 0; i < hittable_acc.get_count(); i++) {
-        if (dev_visit(
-                [&](auto&& arg) {
-                  return arg.hit(ctx, r, 0.001f, closest_so_far, temp_rec,
-                                 temp_material_type);
-                },
-                hittable_acc[i])) {
+        if (dev_visit(visitor::hitable_hit(ctx, r, 0.001f, closest_so_far,
+                                           temp_rec, temp_material_type),
+                      hittable_acc[i])) {
           hit_anything = true;
           closest_so_far = temp_rec.t;
           rec = temp_rec;
@@ -54,18 +46,14 @@ inline auto render_pixel(auto& ctx, int x_coord, int y_coord, camera const& cam,
     color cur_attenuation { 1.0f, 1.0f, 1.0f };
     ray scattered;
     color emitted;
-    material_t material_type;
+    scene::material_t material_type;
     for (auto i = 0; i < depth; i++) {
-      hit_record rec;
+      visitor::hit_record rec;
       if (hit_world(cur_ray, rec, material_type)) {
-        emitted = dev_visit([&](auto&& arg) { return arg.emitted(ctx, rec); },
-                            material_type);
-        if (dev_visit(
-                [&](auto&& arg) {
-                  return arg.scatter(ctx, cur_ray, rec, cur_attenuation,
-                                     scattered);
-                },
-                material_type)) {
+        emitted = dev_visit(visitor::material_emitted(ctx, rec), material_type);
+        if (dev_visit(visitor::material_scatter(ctx, cur_ray, rec,
+                                               cur_attenuation, scattered),
+                      material_type)) {
           // On hitting the object, the ray gets scattered
           cur_ray = scattered;
         } else {
@@ -93,8 +81,8 @@ inline auto render_pixel(auto& ctx, int x_coord, int y_coord, camera const& cam,
 
   color final_color(0.0f, 0.0f, 0.0f);
   for (auto i = 0; i < samples; i++) {
-    const auto u = (x_coord + rng.float_t()) / width;
-    const auto v = (y_coord + rng.float_t()) / height;
+    const auto u = (x_coord + rng.real()) / width;
+    const auto v = (y_coord + rng.real()) / height;
     // u and v are points on the viewport
     ray r = cam.get_ray(u, v, rng);
     final_color += get_color(r);
@@ -112,8 +100,8 @@ inline void executor(sycl::handler& cgh, camera const& cam_ptr,
                      auto& hittable_acc, auto& fb_acc, auto& texture_acc) {
   if constexpr (buildparams::use_single_task) {
     cgh.single_task<PixelRender>([=] {
-      LocalPseudoRNG rng;
-      task_context ctx { rng, texture_acc.get_pointer() };
+      random::LocalPseudoRNG rng;
+      visitor::task_context ctx { rng, texture_acc.get_pointer() };
       for (int x_coord = 0; x_coord != width; ++x_coord)
         for (int y_coord = 0; y_coord != height; ++y_coord) {
           render_pixel<width, height, samples, depth>(
@@ -129,8 +117,8 @@ inline void executor(sycl::handler& cgh, camera const& cam_ptr,
       const auto y_coord = gid[0];
       auto init_generator_state =
           std::hash<std::size_t> {}(item.get_linear_id());
-      LocalPseudoRNG rng(init_generator_state);
-      task_context ctx { rng, texture_acc.get_pointer() };
+      random::LocalPseudoRNG rng(init_generator_state);
+      visitor::task_context ctx { rng, texture_acc.get_pointer() };
       render_pixel<width, height, samples, depth>(
           ctx, x_coord, y_coord, cam_ptr, hittable_acc, fb_acc);
     });
@@ -140,12 +128,12 @@ inline void executor(sycl::handler& cgh, camera const& cam_ptr,
 // Render function to call the render kernel
 template <int width, int height, int samples>
 void render(sycl::queue& queue, sycl::buffer<color, 2>& frame_buf,
-            std::vector<hittable_t>& hittables, camera& cam) {
+            std::vector<scene::hittable_t>& hittables, camera& cam) {
   auto constexpr depth = 50;
   const auto nb_hittable = hittables.size();
-  auto hittables_buf = sycl::buffer<hittable_t, 1>(hittables.data(),
+  auto hittables_buf = sycl::buffer<scene::hittable_t, 1>(hittables.data(),
                                                    sycl::range<1>(nb_hittable));
-  auto texture_buf = image_texture::freeze();
+  auto texture_buf = scene::image_texture::freeze();
 
   // Submit command group on device
   queue.submit([&](sycl::handler& cgh) {
@@ -158,3 +146,4 @@ void render(sycl::queue& queue, sycl::buffer<color, 2>& frame_buf,
                                             texture_acc);
   });
 }
+} // namespace raytracer
\ No newline at end of file
diff --git a/include/rtweekend.hpp b/include/rtweekend.hpp
deleted file mode 100644
index db7618d..0000000
--- a/include/rtweekend.hpp
+++ /dev/null
@@ -1,109 +0,0 @@
-#ifndef RT_SYCL_RTWEEKEND_HPP
-#define RT_SYCL_RTWEEKEND_HPP
-
-#include <algorithm>
-#include <cmath>
-#include <cstdlib>
-#include <limits>
-#include <memory>
-#include <random>
-#include <type_traits>
-#include <variant>
-#include <vector>
-
-#include <sycl.hpp>
-
-#include <build_parameters.hpp>
-#include <xorshift.hpp>
-
-// Constants
-
-constexpr float infinity = std::numeric_limits<float>::infinity();
-constexpr float pi = 3.1415926535897932385f;
-
-// type aliases for float3 - vec, point and color
-using point = sycl::float3;
-using color = sycl::float3;
-using vec = sycl::float3;
-
-// Utility Functions
-
-inline float degrees_to_radians(float degrees) { return degrees * pi / 180.0f; }
-
-class LocalPseudoRNG {
- public:
-  inline LocalPseudoRNG(std::uint32_t init_state = xorshift<>::initial_state)
-      : generator { init_state } {}
-
-  // Returns a random float in 0.f 1.
-  inline float float_t() {
-    constexpr float scale = 1.f / (uint64_t { 1 } << 32);
-    return generator() * scale;
-  }
-
-  // Returns a random float in min, max
-  inline float float_t(float min, float max) {
-    // TODO use FMA ?
-    return min + (max - min) * float_t();
-  }
-
-  // Returns a random vector with coordinates in 0.f 1.
-  inline vec vec_t() { return { float_t(), float_t(), float_t() }; }
-
-  // Returns a random vec with coordinates in min, max
-  inline vec vec_t(float min, float max) {
-    auto scale = max - min;
-    return vec_t() * scale + min;
-  }
-
-  // Returns a random unit vector
-  inline vec unit_vec() {
-    auto x = float_t(-1.f, 1.f);
-    auto maxy = sycl::sqrt(1 - x * x);
-    auto y = float_t(-maxy, maxy);
-    auto absz = sycl::sqrt(maxy * maxy - y * y);
-    auto z = (float_t() > 0.5) ? absz : -absz;
-    return vec(x, y, z);
-  }
-
-  // Returns a random vector in the unit ball of usual norm
-  inline vec in_unit_ball() {
-    // Polar coordinates r, theta, phi
-    auto r = float_t();
-    auto theta = float_t(0, 2 * pi);
-    auto phi = float_t(0, pi);
-
-    auto plan_seed = r * sycl::sin(phi);
-    auto z = r * sycl::cos(phi);
-
-    return { plan_seed * sycl::cos(theta), plan_seed * sycl::sin(theta), z };
-  }
-
-  // Return a random vector in the unit disk of usual norm in plane x, y
-  inline vec in_unit_disk() {
-    auto x = float_t(-1.f, 1.f);
-    auto maxy = sycl::sqrt(1 - x * x);
-    auto y = float_t(-maxy, maxy);
-    return { x, y, 0.f };
-  }
-
- private:
-  xorshift<> generator;
-};
-
-/**
- @brief Used as a poorman's cooperative ersatz of device global variable
-        The task context is (manually) passed through the call stack to all
-        kernel callees
- */
-struct task_context {
-  LocalPseudoRNG rng;
-  // See image_texture in texture.hpp for more details
-  sycl::global_ptr<uint8_t> texture_data;
-};
-
-// Common Headers
-#include "ray.hpp"
-#include "vec.hpp"
-
-#endif
diff --git a/include/sphere.hpp b/include/sphere.hpp
index ccc5b26..ab1155d 100644
--- a/include/sphere.hpp
+++ b/include/sphere.hpp
@@ -2,27 +2,10 @@
 #define SPHERE_H
 
 #include "material.hpp"
-#include "ray.hpp"
-#include "rtweekend.hpp"
+#include "primitives.hpp"
 #include "texture.hpp"
-#include "vec.hpp"
-
-/* Computes normalised values of theta and phi. The input vector p
-corresponds to a vector passing through the centre of the a sphere
-and the hipoint on the surface of the sphere */
-std::pair<float, float> mercator_coordinates(const vec& p) {
-  // phi is the angle around the axis
-  auto phi = sycl::atan2(p.z(), p.x());
-  // theta is the angle down from the pole
-  auto theta = sycl::asin(p.y());
-  // theta and phi together constitute the spherical coordinates
-  // phi is between -pi and pi , u is between 0 and 1
-  auto u = 1 - (phi + pi) / (2 * pi);
-  // theta is between -pi/2 and pi/2 , v is between 0 and 1
-  auto v = (theta + pi / 2) / pi;
-  return { u, v };
-}
 
+namespace raytracer::scene {
 class sphere {
  public:
   sphere() = default;
@@ -55,56 +38,6 @@ class sphere {
       return center0 + ((time - time0) / (time1 - time0)) * (center1 - center0);
   }
 
-  /// Compute ray interaction with sphere
-  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
-           material_t& hit_material_type) const {
-    hit_material_type = material_type;
-
-    /*(P(t)-C).(P(t)-C)=r^2
-    in the above sphere equation P(t) is the point on sphere hit by the ray
-    (A+tb−C)⋅(A+tb−C)=r^2
-    (t^2)b⋅b + 2tb⋅(A−C) + (A−C)⋅(A−C)−r^2 = 0
-    There can 0 or 1 or 2 real roots*/
-    vec oc = r.origin() - center(r.time()); // oc = A-C
-    auto a = sycl::dot(r.direction(), r.direction());
-    auto b = sycl::dot(oc, r.direction());
-    auto c = sycl::dot(oc, oc) - radius * radius;
-    auto discriminant = b * b - a * c;
-    // Real roots if discriminant is positive
-    if (discriminant > 0) {
-      // First root
-      auto temp = (-b - sycl::sqrt(discriminant)) / a;
-      if (temp < max && temp > min) {
-        rec.t = temp;
-        // Ray hits the sphere at p
-        rec.p = r.at(rec.t);
-        vec outward_normal = (rec.p - center(r.time())) / radius;
-        // To set if hit point is on the front face and the outward normal in
-        // rec
-        rec.set_face_normal(r, outward_normal);
-        /* Update u and v values in the hit record. Normal of a
-        point is calculated as above. This vector is used to also
-        used to get the mercator coordinates of the hitpoint.*/
-        std::tie(rec.u, rec.v) = mercator_coordinates(rec.normal);
-        return true;
-      }
-      // Second root
-      temp = (-b + sycl::sqrt(discriminant)) / a;
-      if (temp < max && temp > min) {
-        rec.t = temp;
-        // Ray hits the sphere at p
-        rec.p = r.at(rec.t);
-        vec outward_normal = (rec.p - center(r.time())) / radius;
-        rec.set_face_normal(r, outward_normal);
-        // Update u and v values in the hit record
-        std::tie(rec.u, rec.v) = mercator_coordinates(rec.normal);
-        return true;
-      }
-    }
-    // No real roots
-    return false;
-  }
-
   // Geometry properties
   point center0, center1;
   real_t radius;
@@ -115,5 +48,6 @@ class sphere {
   // Material properties
   material_t material_type;
 };
+} // namespace raytracer::scene
 
 #endif
diff --git a/include/task_context.hpp b/include/task_context.hpp
new file mode 100644
index 0000000..413ec46
--- /dev/null
+++ b/include/task_context.hpp
@@ -0,0 +1,18 @@
+#ifndef TASK_CONTEXT_HPP
+#define TASK_CONTEXT_HPP
+
+#include "localrandom.hpp"
+#include "sycl.hpp"
+
+namespace raytracer::visitor {
+/**
+ @brief Used as a poorman's cooperative ersatz of device global variable
+        The task context is (manually) passed through all the visitors
+ */
+struct task_context {
+  raytracer::random::LocalPseudoRNG rng;
+  // See image_texture in texture.hpp for more details
+  sycl::global_ptr<uint8_t> texture_data;
+};
+} // namespace raytracer::visitor
+#endif
\ No newline at end of file
diff --git a/include/texture.hpp b/include/texture.hpp
index 50002cc..d741cc9 100644
--- a/include/texture.hpp
+++ b/include/texture.hpp
@@ -1,54 +1,42 @@
 #ifndef RT_SYCL_TEXTURE_HPP
 #define RT_SYCL_TEXTURE_HPP
-#include "hitable.hpp"
-#include "rtweekend.hpp"
-#include "vec.hpp"
-#include <array>
-#include <cmath>
-#include <iostream>
-#include <optional>
+#include <tuple>
+#include <variant>
 #include <vector>
 
-#include "sycl.hpp"
-
 #define STB_IMAGE_IMPLEMENTATION
 #include <stb/stb_image.h>
 
+#include "primitives.hpp"
+#include "sycl.hpp"
+
+namespace raytracer::scene {
 // Solid texture consists of a single color
 struct solid_texture {
+ private:
+  color color_value;
+
+ public:
   solid_texture() = default;
   solid_texture(const color& c)
       : color_value { c } {}
-  solid_texture(float red, float green, float blue)
+  solid_texture(real_t red, real_t green, real_t blue)
       : solid_texture { color { red, green, blue } } {}
-  // For solid texture, the color is same throughout the sphere
-  color value(auto&, const hit_record&) const { return color_value; }
-
- private:
-  color color_value;
+  color get_color() const { return color_value; }
 };
 
 // Takes two solid_textures to create checker pattern
 struct checker_texture {
-  checker_texture() = default;
+ private:
+  color odd;
+  color even;
 
-  checker_texture(const solid_texture& x, const solid_texture& y)
-      : odd { x }
-      , even { y } {}
+ public:
+  checker_texture() = default;
   checker_texture(const color& c1, const color& c2)
-      : odd { solid_texture { c1 } }
-      , even { solid_texture { c2 } } {}
-  // Color value is different based on normalised spherical coordinates
-  color value(auto& ctx, const hit_record& rec) const {
-    auto sines = sycl::sin(10 * rec.p.x()) * sycl::sin(10 * rec.p.y()) *
-                 sycl::sin(10 * rec.p.z());
-    if (sines < 0)
-      return odd.value(ctx, rec);
-    else
-      return even.value(ctx, rec);
-  }
-  solid_texture odd;
-  solid_texture even;
+      : odd { c1 }
+      , even { c2 } {}
+  auto get_colors() const { return std::tuple(odd, even); };
 };
 
 /**
@@ -77,10 +65,10 @@ struct image_texture {
   std::size_t offset;
 
   /// The repetition rate of the image
-  float cyclic_frequency { 1.f };
+  real_t cyclic_frequency { 1.f };
 
   image_texture(std::size_t _width, std::size_t _height, std::size_t _offset,
-                float _cyclic_frequency)
+                real_t _cyclic_frequency)
       : width { _width }
       , height { _height }
       , offset { _offset }
@@ -95,7 +83,7 @@ struct image_texture {
          the image in the texture
  */
   static image_texture image_texture_factory(const char* file_name,
-                                             float _cyclic_frequency = 1) {
+                                             real_t _cyclic_frequency = 1) {
     assert(!frozen);
     auto components_per_pixel = bytes_per_pixel;
     uint8_t* _data;
@@ -130,30 +118,37 @@ struct image_texture {
                                       { texture_data.size() / 3, 3 } };
   }
 
-  /// Get the color for the texture at the given place
-  /// \todo rename this value() to color() everywhere?
-  color value(auto& ctx, const hit_record& rec) const {
+  /**
+    @brief Get the color for normalised coordinate
+
+    @param u normalised X coordinate
+    @param v normalised Y coordinate
+    @param base_ptr pointer to the image texture global buffer
+    @return color
+   */
+  color color_at(real_t u, real_t v,
+                 sycl::global_ptr<uint8_t> texture_data) const {
     // If texture data is unavailable, return solid cyan
     // The image is repeated by the repetition factor
 
-    std::size_t i =
-        sycl::fmod(rec.u * cyclic_frequency, (float)1) * (width - 1);
+    std::size_t i = sycl::fmod(u * cyclic_frequency, (real_t)1) * (width - 1);
     // The image frame buffer is going downwards, so flip the y axis
     std::size_t j =
-        (1 - sycl::fmod(rec.v * cyclic_frequency, (float)1)) * (height - 1);
+        (1 - sycl::fmod(v * cyclic_frequency, (real_t)1)) * (height - 1);
     std::size_t local_offset = j * width + i;
     std::size_t pix_idx = local_offset + offset;
     auto scale = 1.f / 255;
-    auto& texture_data = ctx.texture_data;
     return { texture_data[pix_idx * 3] * scale,
              texture_data[pix_idx * 3 + 1] * scale,
              texture_data[pix_idx * 3 + 2] * scale };
   }
 };
 
-using texture_t = std::variant<checker_texture, solid_texture, image_texture>;
+using texture_t =
+    std::variant<std::monostate, checker_texture, solid_texture, image_texture>;
 
 // Start filled with the fallback texture (solid blue) for texture load error
 std::vector<uint8_t> image_texture::texture_data { 0, 0, 1 };
 bool image_texture::frozen = false;
+} // namespace raytracer::scene
 #endif
diff --git a/include/texture_visitor.hpp b/include/texture_visitor.hpp
new file mode 100644
index 0000000..bc199e4
--- /dev/null
+++ b/include/texture_visitor.hpp
@@ -0,0 +1,41 @@
+#ifndef TEXTURE_VISITOR_HPP
+#define TEXTURE_VISITOR_HPP
+#include "hit_record.hpp"
+#include "sycl.hpp"
+#include "primitives.hpp"
+#include "task_context.hpp"
+#include "texture.hpp"
+
+namespace raytracer::visitor {
+struct texture_value {
+ private:
+  task_context& ctx;
+  const hit_record& rec;
+
+ public:
+  texture_value(task_context& ctx, const hit_record& rec)
+      : ctx { ctx }
+      , rec { rec } {}
+
+  // For solid texture, the color is same throughout the object
+  color operator()(scene::solid_texture& t) const { return t.get_color(); }
+
+  // Color value is different based on normalised spherical coordinates
+  color operator()(scene::checker_texture& t) const {
+    auto sines = sycl::sin(10 * rec.p.x()) * sycl::sin(10 * rec.p.y()) *
+                 sycl::sin(10 * rec.p.z());
+    auto const & [odd, even] = t.get_colors();
+    return (sines < 0) ? odd : even;
+  }
+
+  color operator()(scene::image_texture& t) const {
+    return t.color_at(rec.u, rec.v, ctx.texture_data);
+  }
+
+  color operator()(std::monostate) {
+    assert(false && "unreachable");
+    return { 0.f, 0.f, 0.f };
+  }
+};
+} // namespace raytracer::visitor
+#endif
\ No newline at end of file
diff --git a/include/triangle.hpp b/include/triangle.hpp
index b1102c8..2ff65ed 100644
--- a/include/triangle.hpp
+++ b/include/triangle.hpp
@@ -2,122 +2,21 @@
 #define TRIANGLE_HPP
 
 #include "material.hpp"
-#include "ray.hpp"
-#include "rtweekend.hpp"
-#include "texture.hpp"
-#include "vec.hpp"
+#include "primitives.hpp"
 
+namespace raytracer::scene {
 struct _triangle_coord {
   point v0, v1, v2;
 };
 
-inline bool badouel_ray_triangle_intersec(const ray& r,
-                                          _triangle_coord const& tri,
-                                          real_t min, real_t max,
-                                          hit_record& rec) {
-  // Get triangle edge vectors and plane normal
-  auto u = tri.v1 - tri.v0;
-  auto v = tri.v2 - tri.v0;
-  vec outward_normal = sycl::cross(u, v);
-
-  auto w0 = r.origin() - tri.v0;
-  auto a = -sycl::dot(outward_normal, w0);
-  auto b = sycl::dot(outward_normal, r.direction());
-
-  // ray is parallel to triangle plane
-  if (sycl::fabs(b) < 0.000001f)
-    return false;
-
-  // intersection point of ray with triangle
-  real_t length = a / b;
-  if (length < 0)
-    return false;
-  else if (length < min || length > max)
-    return false;
-
-  vec hit_pt = r.at(length);
-  auto uu = sycl::dot(u, u);
-  auto uv = sycl::dot(u, v);
-  auto vv = sycl::dot(v, v);
-  auto w = hit_pt - tri.v0;
-  auto wu = sycl::dot(w, u);
-  auto wv = sycl::dot(w, v);
-  auto D = uv * uv - uu * vv;
-
-  auto s = (uv * wv - vv * wu) / D;
-  auto t = (uv * wu - uu * wv) / D;
-  if (s < 0.0f || s > 1.0f || t < 0.0f || (s + t) > 1.0f)
-    return false;
-
-  rec.set_face_normal(r, outward_normal);
-  rec.t = length;
-  rec.p = hit_pt;
-  return true;
-};
-
-inline bool moller_trumbore_triangle_intersec(const ray& r,
-                                              _triangle_coord const& tri,
-                                              real_t min, real_t max,
-                                              hit_record& rec) {
-  constexpr auto epsilon = 0.0000001f;
-
-  // Get triangle edge vectors and plane normal
-  auto edge1 = tri.v1 - tri.v0;
-  auto edge2 = tri.v2 - tri.v0;
-  vec h = sycl::cross(r.direction(), edge2);
-  auto a = sycl::dot(edge1, h);
-  auto a_abs = sycl::fabs(a);
-
-  if (a_abs < epsilon)
-    return false; // This ray is parallel to this triangle
-
-  auto a_pos = a > 0.f;
-
-  auto s = r.origin() - tri.v0;
-  auto u = sycl::dot(s, h);
-  auto u_pos = u > 0.f;
-
-  if ((u_pos xor a_pos) || sycl::fabs(u) > a_abs)
-    return false;
-
-  auto q = sycl::cross(s, edge1);
-  auto v = sycl::dot(r.direction(), q);
-  auto v_pos = v > 0.f;
-  if ((v_pos xor a_pos) || (sycl::fabs(u + v) > a_abs))
-    return false;
-
-  auto length = sycl::dot(edge2, q) / a;
-
-  if (length < min || length > max)
-    return false;
-
-  vec hit_pt = r.at(length);
-
-  rec.set_face_normal(r, sycl::cross(edge1, edge2));
-  rec.t = length;
-  rec.p = hit_pt;
-  return true;
-};
-
-// A triangle based on 3 points
-template <auto IntersectionStrategy = moller_trumbore_triangle_intersec>
-class _triangle : public _triangle_coord {
+class triangle : public _triangle_coord {
  public:
-  _triangle() = default;
-  _triangle(const point& _v0, const point& _v1, const point& _v2,
+  triangle(const point& _v0, const point& _v1, const point& _v2,
             const material_t& mat_type)
       : _triangle_coord { _v0, _v1, _v2 }
       , material_type { mat_type } {}
 
-  /// Compute ray interaction with triangle
-  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
-           material_t& hit_material_type) const {
-    hit_material_type = material_type;
-    return IntersectionStrategy(r, *this, min, max, rec);
-  }
-
   material_t material_type;
 };
-
-using triangle = _triangle<>;
+} // namespace raytracer::scene
 #endif
diff --git a/include/vec.hpp b/include/vec.hpp
deleted file mode 100644
index f9327ac..0000000
--- a/include/vec.hpp
+++ /dev/null
@@ -1,37 +0,0 @@
-#ifndef RT_SYCL_VEC_HPP
-#define RT_SYCL_VEC_HPP
-
-#include "rtweekend.hpp"
-#include <cmath>
-#include <iostream>
-
-using real_t = float;
-
-// vec Utility Functions
-inline float length_squared(const vec& v) {
-  return sycl::fma(v.x(), v.x(), sycl::fma(v.y(), v.y(), v.z() * v.z()));
-}
-
-inline std::ostream& operator<<(std::ostream& out, const vec& v) {
-  return out << v.x() << ' ' << v.y() << ' ' << v.z();
-}
-
-// Missing operator from the SYCL specification for now
-vec operator-(const vec& u) { return vec(-u.x(), -u.y(), -u.z()); }
-
-// Compute a unit vector from a non-null vector
-inline vec unit_vector(const vec& v) { return v / sycl::length(v); }
-
-// Compute reflected ray's direction
-vec reflect(const vec& v, const vec& n) { return v - 2 * sycl::dot(v, n) * n; }
-
-// Computes refracted ray's direction based on refractive index
-vec refract(const vec& uv, const vec& n, float etai_over_etat) {
-  auto cos_theta = sycl::fmin(-sycl::dot(uv, n), 1.0f);
-  vec r_out_perp = etai_over_etat * (uv + cos_theta * n);
-  vec r_out_parallel =
-      -sycl::sqrt(sycl::fabs(1.0f - length_squared(r_out_perp))) * n;
-  return r_out_perp + r_out_parallel;
-}
-
-#endif
diff --git a/src/main.cpp b/src/main.cpp
index df7d03a..940395d 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -11,8 +11,12 @@
 #define STB_IMAGE_WRITE_IMPLEMENTATION
 #include <stb/stb_image_write.h>
 
+#include "build_parameters.hpp"
 #include "render.hpp"
 
+using namespace raytracer::scene;
+using namespace raytracer;
+
 // Function to save image data in ppm format
 void dump_image_ppm(int width, int height, auto& fb_data) {
   std::cout << "P3\n" << width << " " << height << "\n255\n";
@@ -30,7 +34,7 @@ void dump_image_ppm(int width, int height, auto& fb_data) {
   }
 }
 
-void save_image_png(int width, int height, sycl::buffer<color, 2> &fb) {
+void save_image_png(int width, int height, sycl::buffer<color, 2>& fb) {
   constexpr unsigned num_channels = 3;
   auto fb_data = fb.get_access<sycl::access::mode::read>();
 
@@ -40,7 +44,6 @@ void save_image_png(int width, int height, sycl::buffer<color, 2> &fb) {
   int index = 0;
   for (int j = height - 1; j >= 0; --j) {
     for (int i = 0; i < width; ++i) {
-      auto input_index = j * width + i;
       int r = static_cast<int>(
           256 * std::clamp(sycl::sqrt(fb_data[j][i].x()), 0.0f, 0.999f));
       int g = static_cast<int>(
@@ -73,14 +76,14 @@ int main() {
   hittables.emplace_back(sphere(point { 0, -1000, 0 }, 1000, m));
   t = checker_texture(color { 0.9f, 0.9f, 0.9f }, color { 0.4f, 0.2f, 0.1f });
 
-  LocalPseudoRNG rng;
+  random::LocalPseudoRNG rng;
 
   for (int a = -11; a < 11; a++) {
     for (int b = -11; b < 11; b++) {
       // Based on a random variable , the material type is chosen
-      auto choose_mat = rng.float_t();
+      auto choose_mat = rng.real();
       // Spheres are placed at a point randomly displaced from a,b
-      point center(a + 0.9f * rng.float_t(), 0.2f, b + 0.9f * rng.float_t());
+      point center(a + 0.9f * rng.real(), 0.2f, b + 0.9f * rng.real());
       if (sycl::length((center - point(4, 0.2f, 0))) > 0.9f) {
         if (choose_mat < 0.4f) {
           // Lambertian
@@ -90,13 +93,13 @@ int main() {
         } else if (choose_mat < 0.8f) {
           // Lambertian movig spheres
           auto albedo = rng.vec_t() * rng.vec_t();
-          auto center2 = center + point { 0, rng.float_t(0, 0.25f), 0 };
+          auto center2 = center + point { 0, rng.real(0, 0.25f), 0 };
           hittables.emplace_back(sphere(center, center2, 0.0f, 1.0f, 0.2f,
                                         lambertian_material(albedo)));
         } else if (choose_mat < 0.95f) {
           // metal
           auto albedo = rng.vec_t(0.5f, 1);
-          auto fuzz = rng.float_t(0, 0.5f);
+          auto fuzz = rng.real(0, 0.5f);
           hittables.emplace_back(
               sphere(center, 0.2f, metal_material(albedo, fuzz)));
         } else {