Fix raytracer bugs: BVH traversal, AABB transforms, root selection, and shading

- BVH: transform AABB using all 8 corners, fix leaf node traversal to check all primitives
- Node: reset AABB from primitive before transform, compute distance in world space
- Primitives: correct quadratic root selection (pick smallest positive), fix t-guards for Circle/RectangleXY, fix Torus AABB orientation
- Ray: fix random_unit_vec to cover all octants, compute reflection outside light loop, add indirect diffuse GI

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

README.md

@@ -1,41 +1,155 @@
-
 # Graphics Project
 
-This is my graphics project that I will be working on for A5
-I will use rlua for interacting with lua files
+# Introduction
 
-## Installation
+![image](img/img2.png)
 
-## Scripting
+This is a project I undertook at the University of Waterloo, where I first started using rust. Because of my inexperience, the code isn't as organised as it would be if I made it today and represents my first steps in computer graphics and the rust language.
 
-V()
-P()
-Scene()
-Scene.addNode()
-Scene.addLight()
-Node()
-Node.translate()
-Node.rotate()
-Node.scale()
-Camera()
-Light()
-Material()
-MaterialRed()
-MaterialBlue()
-MaterialGreen()
-MaterialMagenta()
-MaterialTurquoise()
-Sphere()
-SphereUnit()
-Cube()
-CubeUnit()
-Cone()
-ConeUnit()
-Cyclinder()
-//CylinderUnit()
-Circle()
-CircleUnit()
-Rectangle()
-RectangleUnit()
-Steiner()
-Torus()
+My unique aim was to perform ray intersections on special geometric surfaces, such as the CrossCap surface and the Steiner surface, hence those are among my _primitives_.
+
+# Installation
+
+Clone and run with `cargo run`, however much better performance will be granted with `cargo run --release`.
+
+![example](img/example.png)
+
+# Rhai
+
+Rhai is used as an interactive scripting lang for this project. Examples are found in `rhai/`.
+
+## Full List of rhai commands
+
+```
+/// Basic math types
+
+V(x : float, y : float, z : float) -> Vector3
+    // 3‑dimensional vector, used for directions, colors, etc.
+
+P(x : float, y : float, z : float) -> Position3
+    // 3‑dimensional position vector, used for points in space.
+
+
+/// Scene and graph
+
+Scene() -> Scene
+    // Create an empty scene with no nodes, lights, or camera.
+
+Scene.addNode(node : Node) -> void
+    // Add a node (with geometry or camera) to the scene.
+
+Scene.addLight(light : Light) -> void
+    // Add a light source to the scene.
+
+Scene.setCamera(camera : Camera) -> void
+    // Set the active camera for this scene.
+
+
+/// Nodes and transforms
+
+Node() -> Node
+    // Create an empty node (no mesh / camera / light by default).
+
+Node.translate(x : float, y : float, z : float) -> Node
+    // Apply translation by vector V(x, y, z) in local space, returns self for chaining.
+
+Node.rotate(x : float, y : float, z : float) -> Node
+    // Rotate node by Euler angles (in radians or degrees, implementation‑defined).
+
+Node.scale(x : float, y : float, z : float) -> Node
+    // Non‑uniform scale in local space.
+
+Node.setMaterial(material : Material) -> Node
+    // Set material for this node's mesh (if any).
+
+
+/// Camera
+
+Camera(position : P, target : P, up : V) -> Camera
+    // Create a camera located at `position`, looking at `target`, with `up` as the up direction.
+
+
+/// Lighting
+
+Ambient(color : V) -> AmbientLight
+    // Ambient light contribution with RGB in [0, 1].
+
+Light(position : P, color : V, falloff : V) -> PointLight
+    // Point light at `position` with RGB `color` and falloff parameters (constant, linear, quadratic).
+
+
+/// Materials
+
+Material(kd : V, ks : V, kr : V, shininess : float) -> Material
+    // Phong‑style material:
+    //   kd: diffuse color
+    //   ks: specular color
+    //   kr: reflection / mirror color
+    //   shininess: specular exponent.
+
+MaterialRed() -> Material
+MaterialBlue() -> Material
+MaterialGreen() -> Material
+MaterialMagenta() -> Material
+MaterialTurquoise() -> Material
+    // Convenience materials with predefined colors.
+
+
+/// Primitives
+
+Sphere(pos : P, radius : float) -> Mesh
+    // Sphere centered at `pos` with given radius.
+
+SphereUnit() -> Mesh
+    // Unit sphere at (0, 0, 0) with radius 1.
+
+Cube(pos : P, radius : float, normal : V) -> Mesh
+    // Cube centered at `pos`, edge length = 2 * radius (or radius, implementation‑defined),
+    // `normal` can define an orientation axis.
+
+CubeUnit() -> Mesh
+    // Unit cube at (0, 0, 0).
+
+Cone(radius : float, height : float) -> Mesh
+    // Cone aligned with +Z (for example), base radius and height.
+
+ConeUnit() -> Mesh
+    // Cone with radius 1 and height 1 at the origin.
+
+Cylinder(radius : float, height : float) -> Mesh
+    // Cylinder aligned with +Z, given radius and height.
+
+CylinderUnit() -> Mesh
+    // Cylinder with radius 1 and height 1 at the origin.
+
+Circle(position : P, radius : float, normal : V) -> Mesh
+    // Flat disk at `position` with `normal` orientation and given radius.
+
+CircleUnit() -> Mesh
+    // Unit circle in the XY plane at the origin.
+
+Rectangle(position : P, size : V, normal : V) -> Mesh
+    // Axis‑aligned rectangle centered at `position`, width/height from size.x / size.y, oriented by `normal`.
+
+RectangleUnit() -> Mesh
+    // 1x1 rectangle in the XY plane centered at origin.
+
+
+/// Special / parametric surfaces
+
+Steiner() -> Mesh
+    // A Steiner surface with default parameters and resolution.
+
+Torus(radiusMajor : float, radiusMinor : float) -> Mesh
+    // Torus with major and minor radius, centered at origin.
+
+Roman() -> Mesh
+    // Roman surface with default scale and resolution.
+
+CrossCap() -> Mesh
+    // Cross‑cap surface (Boy's surface variant / projective plane immersion).
+
+Gnonom() -> Mesh
+    // Gnomon‑like parametric surface (implementation‑defined shape).
+
+```

rhai/bvh_test.rhai

@@ -0,0 +1,81 @@
+// BVH Stress Test Scene - Grid of spheres
+let scene = Scene();
+
+// Camera
+let camera = Camera(P(0.0, 2.0, 5.0), P(0.0, 0.0, 0.0), V(0.0, 1.0, 0.0));
+scene.addCamera("main", camera);
+
+// Lights
+let light = Light(P(3.0, 5.0, 3.0), V(0.8, 0.8, 0.8), V(0.05, 0.05, 0.05));
+light.active(true);
+scene.addLight("key", light);
+
+let light2 = Light(P(-3.0, 4.0, -2.0), V(0.4, 0.4, 0.6), V(0.05, 0.05, 0.05));
+light2.active(true);
+scene.addLight("fill", light2);
+
+let ambient = Ambient(V(0.1, 0.1, 0.1));
+ambient.active(true);
+scene.addLight("ambient", ambient);
+
+// Materials
+let kd = V(0.8, 0.2, 0.2);
+let ks = V(0.3, 0.3, 0.3);
+let kr = V(0.0, 0.0, 0.0);
+let red = Material(kd, ks, kr, 10.0);
+scene.addMaterial("red", red);
+
+let kd = V(0.2, 0.2, 0.8);
+let ks = V(0.3, 0.3, 0.3);
+let kr = V(0.0, 0.0, 0.0);
+let blue = Material(kd, ks, kr, 10.0);
+scene.addMaterial("blue", blue);
+
+let kd = V(0.2, 0.8, 0.2);
+let ks = V(0.3, 0.3, 0.3);
+let kr = V(0.0, 0.0, 0.0);
+let green = Material(kd, ks, kr, 10.0);
+scene.addMaterial("green", green);
+
+let kd = V(0.8, 0.8, 0.2);
+let ks = V(0.3, 0.3, 0.3);
+let kr = V(0.0, 0.0, 0.0);
+let yellow = Material(kd, ks, kr, 10.0);
+scene.addMaterial("yellow", yellow);
+
+let kd = V(0.8, 0.4, 0.8);
+let ks = V(0.3, 0.3, 0.3);
+let kr = V(0.1, 0.1, 0.1);
+let purple = Material(kd, ks, kr, 15.0);
+scene.addMaterial("purple", purple);
+
+// Floor
+let floor = RectangleUnit();
+let floor_node = Node(floor, green);
+floor_node.rotate(-90.0, 0.0, 0.0);
+floor_node.translate(0.0, -1.5, 0.0);
+floor_node.scale(5.0, 5.0, 1.0);
+scene.addNode("floor", floor_node);
+
+// Grid of spheres: 5 x 4 x 5 = 100 spheres
+let materials = [red, blue, green, yellow, purple];
+let count = 0;
+
+for x in range(-2, 3) {
+    for y in range(0, 4) {
+        for z in range(-2, 3) {
+            let mat_idx = count % 5;
+            let mat = materials[mat_idx];
+            let sphere = Sphere(P(0.0, 0.0, 0.0), 0.15);
+            let node = Node(sphere, mat);
+            let px = x.to_float() * 0.7;
+            let py = y.to_float() * 0.7 - 1.0;
+            let pz = z.to_float() * 0.7 - 1.0;
+            node.translate(px, py, pz);
+            scene.addNode("s" + count, node);
+            count += 1;
+        }
+    }
+}
+
+scene

rhai/cow.rhai

@@ -1,35 +0,0 @@
-
-let scene = Scene();
-
-let distance = 3.0;
-let falloff = V(0.0,0.0,0.1);
-let camera = Camera( P(0.0,0.0,distance), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("+Z Cam", camera);
-let camera = Camera( P(0.0,distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("+Y Cam", camera);
-let camera = Camera( P(distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("+X Cam", camera);
-let camera = Camera( P(0.0,0.0,-distance), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-Z Cam", camera);
-let camera = Camera( P(0.0,-distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-Y Cam", camera);
-let camera = Camera( P(-distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-X Cam", camera);
-
-let height = 4.0;
-let spacing = 13.0;
-let light = Light(P(0.0,height,spacing), V(0.0,0.3,0.3), falloff);
-scene.addLight("blue", light);
-let light = Light(P(0.0,height,0.0), V(0.0,0.6,0.0), falloff);
-scene.addLight("green", light);
-let light = Light(P(0.0,height,-spacing), V(0.3,0.0,0.0), falloff);
-scene.addLight("red", light);
-
-let material = Material(V(0.2,0.2,0.2), V(0.2, 0.8, 0.8), 10.0);
-scene.addMaterial("bluegreen", material);
-
-let mesh = Mesh("obj/cow.obj", material);
-let mesh_node = Node(mesh);
-scene.addNode("mesh", mesh_node);
-
-scene

rhai/mesh.rhai

@@ -1,66 +0,0 @@
-let scene = Scene();
-
-let material = Material(V(0.2,0.2,0.2), V(0.2, 0.2, 0.2), 10.0);
-scene.addMaterial("material", material);
-let material2 = Material(V(0.2,0.7,0.2), V(0.2, 0.2, 0.2), 10.0);
-scene.addMaterial("mat2", material2);
-
-let camera = Camera(P(0.0,0.0,2.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("Cam", camera);
-
-let falloff = V(0.0,0.0,0.0);
-let light = Light(P(6.0,6.0,6.0), V(0.4,0.4,0.4), falloff);
-light.active(false);
-scene.addLight("white", light);
-let light = Light(P(2.0,0.0,0.0), V(0.0,1.0,0.0), V(0.1, 0.01, 0.001));
-light.active(false);
-scene.addLight("green", light);
-let light = Light(P(-2.0,0.0,0.0), V(1.0,0.0,0.0), V(0.1, 0.01, 0.001));
-light.active(false);
-scene.addLight("red", light);
-
-let light = Ambient(V(0.3,0.3,0.3));
-scene.addLight("ambient", light);
-
-let tri = TriangleUnit();
-let tri_node = Node(tri, material);
-tri_node.active(false);
-scene.addNode("tri", tri_node);
-
-let circle = CircleUnit();
-let circle_node = Node(circle, material);
-circle_node.active(false);
-scene.addNode("circle", circle_node);
-
-let cone = ConeUnit();
-let cone_node = Node(cone, material);
-cone_node.active(false);
-scene.addNode("cone", cone_node);
-
-let torus = Torus(0.5, 1.5);
-let torus_node = Node(torus, material);
-torus_node.active(true);
-scene.addNode("torus", torus_node);
-
-
-let sphere = SphereUnit();
-let sphere_node = Node(sphere, material);
-sphere_node.translate(0.0,0.0,0.0);
-sphere_node.active(false);
-scene.addNode("sphere", sphere_node);
-
-let ground = SphereUnit();
-let ground_node = Node(ground, material2);
-let scale = 2.0;
-ground_node.translate(0.0,-scale*2.0,0.0);
-ground_node.scale(scale,scale,scale);
-ground_node.active(false);
-scene.addNode("ground", ground_node);
-
-let mesh = Mesh("obj/cat.obj");
-let mesh_node = Node(mesh, material);
-mesh_node.active(false);
-scene.addNode("mesh", mesh_node);
-
-
-scene

rhai/scene.rhai

@@ -1,66 +1,174 @@
 let scene = Scene();
 
-let distance = 10.0;
+let distance = 0.99;
 let camera = Camera( P(0.0,0.0,distance), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
 scene.addCamera("+Z Cam", camera);
-let camera = Camera( P(0.0,distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("+Y Cam", camera);
-let camera = Camera( P(distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("+X Cam", camera);
-let camera = Camera( P(0.0,0.0,-distance), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-Z Cam", camera);
-let camera = Camera( P(0.0,-distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-Y Cam", camera);
-let camera = Camera( P(-distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-X Cam", camera);
+// let camera = Camera( P(0.0,distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
+// scene.addCamera("+Y Cam", camera);
+// let camera = Camera( P(distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
+// scene.addCamera("+X Cam", camera);
+// let camera = Camera( P(0.0,0.0,-distance), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
+// scene.addCamera("-Z Cam", camera);
+// let camera = Camera( P(0.0,-distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
+// scene.addCamera("-Y Cam", camera);
+// let camera = Camera( P(-distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
+// scene.addCamera("-X Cam", camera);
 
-let material = Material(V(0.2,0.2,0.2), V(0.2, 0.8, 0.8), 10.0);
-scene.addMaterial("bluegreen", material);
+let falloff = V(0.1, 0.1, 0.15);
 
-let light = Light(P(0.0,7.0,0.0), V(0.0,0.0,1.0), V(0.1, 0.01, 0.001));
-light.active(false);
+let colour = V(0.0,0.5,0.5);
+let pos = P(-0.5,0.9,0.5);
+let light = Light(pos, colour, falloff);
+light.active(true);
 scene.addLight("blue", light);
 
-let light = Light( P(2.0,7.0,0.0), V(0.0,1.0,0.0), V(0.1, 0.01, 0.001));
+let colour = V(0.0,1.0,0.0);
+let pos = P(-0.5,0.9,-0.5);
+let light = Light(pos, colour, falloff);
 light.active(false);
 scene.addLight("green", light);
 
-let light = Light( P(2.0,7.0,2.0), V(1.0,0.0,0.0), V(0.1, 0.01, 0.001));
+let colour = V(1.0,0.0,0.0);
+let light = Light(pos, colour, falloff);
+light.active(false);
 scene.addLight("red", light);
 
+let colour = V(0.7,0.7,0.7);
+let pos = P(0.0,0.9,0.0);
+let light = Light(pos, colour, falloff);
+light.active(true);
+scene.addLight("white", light);
+
 let light = Ambient(V(0.1,0.1,0.1));
+light.active(true);
 scene.addLight("ambient", light);
 
-
-let sphere = Sphere(P(0.0,0.0,0.0), 1.0 );
-let sphere_node = Node(sphere, material);
- scene.addNode("sphere", sphere_node);
-
 //let mesh = Mesh("obj/cow.obj" );
 //let mesh_node = Node(mesh);
 //scene.addNode("mesh", mesh_node);
 
-for i in 0..6 {
-    let sphere = Sphere(P(0.0,0.0,0.0), 2.0 );
-     let sphere_node = Node(sphere, material);
-     sphere_node.translate(4.0*cos(i.to_float()), -4.0, 4.0*sin(i.to_float()));
-     scene.addNode(i.to_string(), sphere_node);
- }
 // let child = sphere_node.child(sphere);
 // child.translate(V(1.0,1.0,1.0));
 //scene.addNode(child);
 
+//let sphere2= SphereUnit();
+//let sphere2_node = Node( sphere2, material2);
+// sphere2_node.rotate(0.1,0.1,0.0);
+// sphere2_node.translate(0.0,1.0,0.0);
+//scene.addNode("sphere2",  sphere2_node);
+
+let kd = V(1.0, 1.0, 1.0); // Diffuse color (white)
+let ks = V(0.0,0.0,0.0); // Specular color (no specular reflection)
+let kr = V(0.0,0.0,0.0); // Reflection color (no reflection)
+let white_wall = Material(kd, ks, kr, 10.0);
+scene.addMaterial("white_wall", white_wall);
+
+let kd = V(1.0, 0.0, 0.0); // Diffuse color (white)
+let ks = V(0.0,0.0,0.0); // Specular color (no specular reflection)
+let kr = V(0.0,0.0,0.0); // Reflection color (no reflection)
+let red_wall = Material(kd, ks, kr, 10.0);
+scene.addMaterial("red_wall", red_wall);
+
+let kd = V(0.0, 1.0, 0.0); // Diffuse color (white)
+let ks = V(0.0,0.0,0.0); // Specular color (no specular reflection)
+let kr = V(0.0,0.0,0.0); // Reflection color (no reflection)
+let green_wall = Material(kd, ks, kr, 10.0);
+scene.addMaterial("green_wall", green_wall);
+
+let kd = V(0.0, 0.0, 1.0); // Diffuse color (white)
+let ks = V(0.0,0.0,0.0); // Specular color (no specular reflection)
+let kr = V(0.0,0.0,0.0); // Reflection color (no reflection)
+let blue_wall = Material(kd, ks, kr, 10.0);
+scene.addMaterial("blue_wall", blue_wall);
+
+
+//Rear wall
+let rectangle1 = RectangleUnit();
+let rectangle_node1 = Node(rectangle1, white_wall);
+rectangle_node1.rotate(0.0, 0.0, 0.0);
+rectangle_node1.translate(0.0, 0.0, -1.0);
+rectangle_node1.active(true);
+scene.addNode("rectangle1", rectangle_node1);
+
+//Behind wall
+// let rectangle6 = RectangleUnit();
+// let rectangle_node6 = Node(rectangle6, white_wall);
+// rectangle_node6.rotate(0.0, 180.0, 0.0);
+// rectangle_node6.translate(0.0, 0.0, 1.0);
+// rectangle_node6.active(true);
+// scene.addNode("rectangle6", rectangle_node6);
+
+//Right wall
+let rectangle2 = RectangleUnit();
+let rectangle_node2 = Node(rectangle2, green_wall);
+rectangle_node2.rotate(0.0, -90.0, 0.0);
+rectangle_node2.translate(1.0, 0.0, 0.0);
+rectangle_node2.active(true);
+scene.addNode("rectangle2", rectangle_node2);
+
+//Floor
+let rectangle3 = RectangleUnit();
+let rectangle_node3 = Node(rectangle3, red_wall);
+rectangle_node3.rotate(0.0, 90.0, 0.0);
+rectangle_node3.translate(-1.0, 0.0, 0.0);
+rectangle_node3.active(true);
+scene.addNode("rectangle3", rectangle_node3);
+
+//Left wall
+let rectangle4 = RectangleUnit();
+let rectangle_node4 = Node(rectangle4, white_wall);
+rectangle_node4.rotate(90.0, 0.0, 0.0);
+rectangle_node4.translate(0.0, 1.0, 0.0);
+rectangle_node4.active(true);
+scene.addNode("rectangle4", rectangle_node4);
+
+//Ceiling
+let rectangle5 = RectangleUnit();
+let rectangle_node5 = Node(rectangle5, white_wall);
+rectangle_node5.rotate(-90.0, 0.0, 0.0);
+
+rectangle_node5.translate(0.0, -1.0, 0.0);
+rectangle_node5.active(true);
+scene.addNode("rectangle5", rectangle_node5);
+
+
+let kd = V(0.0, 0.0, 0.0); // Diffuse color (white)
+let ks = V(0.0,0.0,0.0); // Specular color (no specular reflection)
+let kr = V(1.0,1.0,1.0); // Reflection color (no reflection)
+let reflective = Material(kd, ks, kr, 10.0);
+scene.addMaterial("reflective", reflective);
+
+let sphere = Sphere(P(0.0,0.0,0.0), 0.4 );
+let sphere_node = Node( sphere, reflective);
+sphere_node.translate(0.4, -0.6, 0.0);
+scene.addNode("sphere",sphere_node);
+
+let kd = V(0.3, 0.3, 0.3); // Diffuse color (white)
+let ks = V(0.3,0.3,0.0); // Specular color (no specular reflection)
+let kr = V(0.0,0.0,1.0); // Reflection color (no reflection)
+let shiny = Material(kd, ks, kr, 2.0);
+scene.addMaterial("shiny", shiny);
+
 let cube = CubeUnit();
-let cube_node = Node(cube, material);
-scene.addNode("cube", cube_node);
+let cube_node = Node( cube, shiny);
+cube_node.translate(-0.5,-0.6,0.0);
+cube_node.scale(0.3,0.2,0.2);
+cube_node.rotate(0.0,45.0,30.0);
+scene.addNode("cube",cube_node);
 
-//let gnonom = Gnonom();
-//let gnonom_node = Node(gnonom);
-//scene.addNode("gnonom", gnonom_node);
+let gnonom = Gnonom();
+let gnonom_node = Node(gnonom, shiny);
+gnonom_node.scale(0.2,0.2,0.2);
+gnonom_node.translate(0.0, 0.-0.7, 0.8);
+gnonom_node.rotate(0.0, 45.0, 0.0);
+gnonom_node.active(false);
+scene.addNode("gnonom", gnonom_node);
 
-//let cylinder = Cylinder(2.0,1.0 );
-//let cylinder_node = Node(cylinder);
-//cylinder_node.scale(1.0,1.0,1.0);
-//scene.addNode("cylinder",cylinder_node);
+// let cylinder = Cylinder(2.0, 1.0);
+// let cylinder_node = Node(cylinder, material);
+// cylinder_node.scale(1.0, 1.0, 1.0);
+// scene.addNode("cylinder", cylinder_node);
 
-scene
+//let cone 
+
+scene

rhai/shapes.rhai

@@ -2,7 +2,6 @@
 let scene = Scene();
 
 let distance = 3.0;
-let falloff = V(0.0,0.0,0.1);
 let camera = Camera( P(0.0,0.0,distance), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
 scene.addCamera("+Z Cam", camera);
 let camera = Camera( P(0.0,distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
@@ -16,8 +15,29 @@ scene.addCamera("-Y Cam", camera);
 let camera = Camera( P(-distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
 scene.addCamera("-X Cam", camera);
 
+let kd = V(0.6, 0.3, 0.9); // Diffuse color (white)
+let ks = V(0.0,0.0,0.0); // Specular color (no specular reflection)
+let kr = V(1.0,1.0,1.0); // Reflection color (no reflection)
+let material = Material(kd, ks, kr, 10.0);
+scene.addMaterial("mattee", material);
+
+let kd = V(0.3, 0.3, 0.3); // Diffuse color (white)
+let ks = V(0.7,0.7,0.7); // Specular color (no specular reflection)
+let kr = V(1.0,1.0,1.0); // Reflection color (no reflection)
+let material1 = Material(kd, ks, kr, 10.0);
+scene.addMaterial("mattee", material1);
+
+let kd = V(0.4, 0.0, 0.8); // Diffuse color (white)
+let ks = V(0.0,0.7,0.7); // Specular color (no specular reflection)
+let kr = V(1.0,1.0,1.0); // Reflection color (no reflection)
+let material2 = Material(kd, ks, kr, 10.0);
+scene.addMaterial("reflect", material2);
+
 let height = 4.0;
 let spacing = 4.0;
+
+let falloff = V(0.0,0.0,0.01);
+
 let blue = V(0.0,0.0,0.6);
 let light = Light(P(0.0,height,spacing), blue, falloff);
 scene.addLight("blue", light);
@@ -28,36 +48,43 @@ let red = V(0.6,0.0,0.0);
 let light = Light(P(0.0,height,-spacing), red, falloff);
 scene.addLight("red", light);
 
-let material = Material(V(0.2,0.2,0.2), V(0.2, 0.8, 0.8), 10.0);
-scene.addMaterial("bluegreen", material);
-
 let steiner = Steiner();
-let steiner_node = Node(steiner, material);
+let steiner_node = Node(steiner, material2);
+steiner_node.rotate(90.0,0.0,0.0);
+steiner_node.translate(0.0,0.0,1.0);
  scene.addNode("steiner", steiner_node);
 
 let steiner2 = Steiner2();
 let steiner2_node = Node(steiner2, material2);
+steiner2_node.active(false);
  scene.addNode("steiner2", steiner2_node);
 
 let crosscap = CrossCap();
 let crosscap_node = Node(crosscap, material);
+crosscap_node.active(false);
  scene.addNode("crosscap", crosscap_node);
 
-let p = 1.0;
-let q = 1.0;
+let p = 0.9;
+let q = 0.1;
 let crosscap2 = CrossCap2(p, q);
 let crosscap2_node = Node(crosscap2, material);
+crosscap2_node.active(true);
+crosscap2_node.translate(0.0,0.0,-1.5);
+crosscap2_node.rotate(140.0,0.0,90.0);
  scene.addNode("crosscap2", crosscap2_node);
 
-let k = 0.5;
+let k = 2.0;
 let roman = Roman(k );
 let roman_node = Node(roman, material);
+roman_node.active(false);
  scene.addNode("roman", roman_node);
 
 let inner_rad = 1.0;
-let outer_rad = 0.5;
+let outer_rad = 1.2;
 let torus = Torus(inner_rad, outer_rad );
 let torus_node = Node(torus, material);
+torus_node.scale(0.2,0.2,0.2);
+torus_node.rotate(0.0,70.0,0.0);
  scene.addNode("torus", torus_node);
 
 scene

rhai/space.rhai

@@ -0,0 +1,34 @@
+let scene = Scene();
+
+let falloff = V(0.0,0.0,0.0);
+
+//CAMERAS
+// let camera = Camera( P(100.0,100.0,100.0), P(500.0,500.0,500.0), V(0.0,1.0,0.0));
+// scene.addCamera("Main Camera", camera);
+
+// //Light for the sun
+// let light = Light(P(800.0, 800.0, 250.0), V(1.0, 1.0, 0.929), falloff);
+// scene.addLight("Sun", light);
+
+// //Ball for the sun 
+// let material = Material(V(1.0, 1.0, 0.9), V(0.9, 0.9, 0.9), 10.0);
+// scene.addMaterial("material_sun", material);
+// let sphere = Sphere(P(800.0, 800.0, 200.0), 50.0);
+// let sphere_node = Node(sphere, material);
+//  scene.addNode("sphere", sphere_node);
+
+// //Ball for the planet 
+// let material = Material(V(0.2,0.8,0.2), V(0.2, 0.8, 0.8), 10.0);
+// scene.addMaterial("material_planet", material);
+// let sphere = Sphere(P(500.0, 500.0, 500.0), 50.0);
+// let sphere_node = Node(sphere, material);
+//  scene.addNode("sphere", sphere_node);
+
+let material = 
+
+let steiner = Steiner();
+let steiner_node = Node(steiner, material);
+ scene.addNode("sphere", steiner);
+
+
+scene

rhai/sphere.rhai

@@ -1,34 +0,0 @@
-let scene = Scene();
-
-let distance = 3.0;
-let falloff = V(0.0,0.0,0.1);
-let camera = Camera( P(0.0,0.0,distance), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("+Z Cam", camera);
-let camera = Camera( P(0.0,distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("+Y Cam", camera);
-let camera = Camera( P(distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("+X Cam", camera);
-let camera = Camera( P(0.0,0.0,-distance), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-Z Cam", camera);
-let camera = Camera( P(0.0,-distance,0.1), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-Y Cam", camera);
-let camera = Camera( P(-distance,0.0,0.0), P(0.0,0.0,0.0), V(0.0,1.0,0.0));
-scene.addCamera("-X Cam", camera);
-
-let height = 4.0;
-let spacing = 13.0;
-let light = Light(P(0.0,height,spacing), V(0.0,0.3,0.3), falloff);
-scene.addLight("blue", light);
-let light = Light(P(0.0,height,0.0), V(0.0,0.6,0.0), falloff);
-scene.addLight("green", light);
-let light = Light(P(0.0,height,-spacing), V(0.3,0.0,0.0), falloff);
-scene.addLight("red", light);
-
-let material = Material(V(0.2,0.2,0.2), V(0.2, 0.8, 0.8), 10.0);
-scene.addMaterial("bluegreen", material);
-
-let sphere = Steiner( material);
-let sphere_node = Node(sphere);
- scene.addNode("sphere", sphere_node);
-
-scene

src/bvh.rs

@@ -1,75 +1,85 @@
-use crate::{ray::*, EPSILON};
-use nalgebra::{Point3, Vector3};
+use crate::{node::Node, ray::*, EPSILON};
+use nalgebra::{distance, Matrix4, Point3, Vector3};
+use std::collections::HashMap;
+use std::fmt;
+
+// Debuging statics
+static mut STATIC0: i32 = 0;
+static mut STATIC1: i32 = 0;
+static mut STATIC2: i32 = 0;
+static mut STATIC3: i32 = 0;
+static mut STATIC4: i32 = 0;
 
 // BOUNDING BOX -----------------------------------------------------------------
 #[derive(Clone)]
 pub struct AABB {
     pub bln: Point3<f64>,
     pub trf: Point3<f64>,
+    pub centroid: Point3<f64>,
 }
 
 impl AABB {
     // New box with respective coordinates
     pub fn new(bln: Point3<f64>, trf: Point3<f64>) -> AABB {
-        let bln = bln + Vector3::new(EPSILON, EPSILON, EPSILON);
-        let trf = trf - Vector3::new(EPSILON, EPSILON, EPSILON);
-        AABB { bln, trf }
+        let bln = bln - Vector3::new(EPSILON, EPSILON, EPSILON);
+        let trf = trf + Vector3::new(EPSILON, EPSILON, EPSILON);
+        let centroid = bln + (trf - bln) / 2.0;
+        AABB { bln, trf, centroid }
+    }
+    //Empty box
+    pub fn empty() -> AABB {
+        AABB {
+            bln: Point3::new(f64::MAX, f64::MAX, f64::MAX),
+            trf: Point3::new(f64::MIN, f64::MIN, f64::MIN),
+            centroid: Point3::new(0.0, 0.0, 0.0),
+        }
+    }
+    //Apply a matrix transformation to a box
+    pub fn transform_mut(&mut self, mat: &Matrix4<f64>) {
+        let corners = [
+            Point3::new(self.bln.x, self.bln.y, self.bln.z),
+            Point3::new(self.trf.x, self.bln.y, self.bln.z),
+            Point3::new(self.bln.x, self.trf.y, self.bln.z),
+            Point3::new(self.trf.x, self.trf.y, self.bln.z),
+            Point3::new(self.bln.x, self.bln.y, self.trf.z),
+            Point3::new(self.trf.x, self.bln.y, self.trf.z),
+            Point3::new(self.bln.x, self.trf.y, self.trf.z),
+            Point3::new(self.trf.x, self.trf.y, self.trf.z),
+        ];
+        let mut new_bln = Point3::new(f64::MAX, f64::MAX, f64::MAX);
+        let mut new_trf = Point3::new(f64::MIN, f64::MIN, f64::MIN);
+        for corner in &corners {
+            let t = mat.transform_point(corner);
+            new_bln = new_bln.inf(&t);
+            new_trf = new_trf.sup(&t);
+        }
+        self.bln = new_bln;
+        self.trf = new_trf;
+        self.centroid = self.bln + (self.trf - self.bln) / 2.0;
     }
     // Intersect bounding box exactly
-    pub fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        let bln = &self.bln;
-        let trf = &self.trf;
-        let t1 = (bln - ray.a).component_div(&ray.b);
-        let t2 = (trf - ray.a).component_div(&ray.b);
+    pub fn intersect_ray(&self, ray: &Ray) -> bool {
+        let t1 = (self.bln - ray.a).component_div(&ray.b);
+        let t2 = (self.trf - ray.a).component_div(&ray.b);
 
-        let tmin = t1.inf(&t2).min();
-        let tmax = t1.sup(&t2).max();
+        let tmin = t1.inf(&t2).max();
+        let tmax = t1.sup(&t2).min();
 
-        if tmax >= tmin {
-            let intersect = ray.at_t(tmin);
-
-            // Check if the intersection is inside the box
-            if intersect.x > bln.x
-                || intersect.x < trf.x
-                || intersect.y > bln.y
-                || intersect.y < trf.y
-                || intersect.z > bln.z
-                || intersect.z < trf.z
-            {
-                return true; // Intersection is outside the box
-            }
-        }
-        false
+        tmax >= tmin && tmax > 0.0
     }
-    // Intersect way with some epsilon term
-    pub fn intersect_bounding_box_aprox(&self, ray: &Ray) -> bool {
-        let bln = &self.bln;
-        let trf = &self.trf;
-        let t1 = (bln - ray.a).component_div(&ray.b);
-        let t2 = (trf - ray.a).component_div(&ray.b);
+    // Intersect with some epsilon tolerance
+    pub fn intersect_ray_aprox(&self, ray: &Ray) -> bool {
+        let t1 = (self.bln - ray.a).component_div(&ray.b);
+        let t2 = (self.trf - ray.a).component_div(&ray.b);
 
-        let tmin = t1.inf(&t2).min();
-        let tmax = t1.sup(&t2).max();
+        let tmin = t1.inf(&t2).max();
+        let tmax = t1.sup(&t2).min();
 
-        if tmax >= tmin {
-            let intersect = ray.at_t(tmin);
-
-            // Check if the intersection is inside the box
-            if intersect.x > bln.x - EPSILON
-                || intersect.x < trf.x + EPSILON
-                || intersect.y > bln.y - EPSILON
-                || intersect.y < trf.y + EPSILON
-                || intersect.z > bln.z - EPSILON
-                || intersect.z < trf.z + EPSILON
-            {
-                return true; // Intersection is outside the box
-            }
-        }
-        false
+        tmax >= tmin - EPSILON && tmax > 0.0
     }
     // Get the center of this bounding box
     fn get_centroid(&self) -> Point3<f64> {
-        self.bln + (self.trf - self.bln) / 2.0
+        self.centroid
     }
     // Make a new AABB that contains both
     pub fn join(&self, other: &AABB) -> AABB {
@@ -86,6 +96,20 @@ impl AABB {
             ),
         )
     }
+    //Join mutably
+    pub fn join_mut(&mut self, other: &AABB) {
+        self.bln = Point3::new(
+            self.bln.x.min(other.bln.x),
+            self.bln.y.min(other.bln.y),
+            self.bln.z.min(other.bln.z),
+        );
+        self.trf = Point3::new(
+            self.trf.x.max(other.trf.x),
+            self.trf.y.max(other.trf.y),
+            self.trf.z.max(other.trf.z),
+        );
+        self.centroid = self.bln + (self.trf - self.bln) / 2.0;
+    }
     //Grow the AABB to contain the cover the point
     pub fn grow(&self, other: &Point3<f64>) -> AABB {
         AABB::new(
@@ -101,34 +125,322 @@ impl AABB {
             ),
         )
     }
+    //Grow mutably
+    pub fn grow_mut(&mut self, other: &Point3<f64>) {
+        self.bln = Point3::new(
+            self.bln.x.min(other.x),
+            self.bln.y.min(other.y),
+            self.bln.z.min(other.z),
+        );
+        self.trf = Point3::new(
+            self.trf.x.max(other.x),
+            self.trf.y.max(other.y),
+            self.trf.z.max(other.z),
+        );
+        self.centroid = self.bln + (self.trf - self.bln) / 2.0;
+    }
     // Size of AABB
     pub fn size(&self) -> Vector3<f64> {
         self.trf - self.bln
-    }
-    //Surface area of AABB
+    } //Surface area of AABB
     pub fn surface_area(&self) -> f64 {
         let size = self.size();
         2.0 * (size.x * size.y + size.x * size.z + size.y * size.z)
     }
+    pub fn area(&self) -> f64 {
+        let extent = self.trf - self.bln;
+        return extent.x * extent.y + extent.y * extent.z + extent.z * extent.x;
+    }
     // Volume of the AABB
     pub fn volume(&self) -> f64 {
         let size = self.size();
         size.x * size.y * size.z
     }
 }
-
-pub enum BVHNode<'a> {
-    Leaf {
-        parent: &'a BVHNode<'a>,
-        bounding_box: AABB,
-        depth: u32,
-    },
-    Node {
-        parent: Option<&'a BVHNode<'a>>,
-        child_l: &'a BVHNode<'a>,
-        child_r: &'a BVHNode<'a>,
-        depth: u32,
-    },
+impl fmt::Display for AABB {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        if self.bln[0] == f64::MAX || self.trf[0] == f64::MIN {
+            writeln!(f, "Empty aabb")
+        } else {
+            writeln!(f, "bln: {}\ntrf: {}", self.bln, self.trf)
+        }
+    }
+}
+#[derive(Clone)]
+pub struct BVHNode {
+    aabb: AABB,        //The nodes bounding box
+    l_idx: usize,      //Child node l, the right node is alway l_idx + 1
+    first_prim: usize, //First primitive that the node encapsulates
+    prim_count: usize, //Number of primitives the node encapsulates
 }
 
-impl<'a> BVHNode<'a> {}
+impl BVHNode {
+    pub fn default() -> BVHNode {
+        BVHNode {
+            aabb: AABB::empty(),
+            l_idx: 0,
+            first_prim: 0,
+            prim_count: 0,
+        }
+    }
+}
+
+impl fmt::Display for BVHNode {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        writeln!(f, "l_idx: {}", self.l_idx)?;
+        writeln!(f, "First Prim: {}", self.first_prim)?;
+        writeln!(f, "Prim Count: {}", self.prim_count)?;
+        writeln!(f, "aabb: {}", self.aabb)
+    }
+}
+
+pub struct BVH {
+    bvh_nodes: Vec<BVHNode>, //BVH nodes with AABBs
+    nodes: Vec<Node>,        //Nodes with primitives
+    nodes_used: usize,
+}
+
+impl BVH {
+    //Build a bvh by subdividing recursively
+    pub fn build(in_nodes: &HashMap<String, Node>) -> BVH {
+        /*
+        Make our own vec of nodes so that we can refer to it by index
+        This might be expensive so another method is preferred
+        */
+        let mut nodes = vec![];
+        for (_, node) in in_nodes {
+            nodes.push(node.clone());
+        }
+
+        //A BVH tree will be maximum size of 2*n + 1
+        //Initialise an empty BVHNode with empty AABB
+        let n = nodes.len();
+        let bvh_nodes: Vec<BVHNode> = vec![BVHNode::default(); 2 * n + 1];
+
+        //Begin constructing our BVH tree
+        //One node used to begin with (The root node)
+        let nodes_used = 1;
+        let mut tree = BVH {
+            nodes,
+            bvh_nodes,
+            nodes_used,
+        };
+        // Get the root node at index 0
+        let root = &mut tree.bvh_nodes[0];
+        root.l_idx = 0; //Root node has no left or right child to begin
+        (root.first_prim, root.prim_count) = (0, n); //Make root include all n nodes
+        tree.update_bvh_node_aabb(0); //Create the root nodes AABB on the n primitives
+        tree.subdivide(0); //Sub divide the root node
+        tree
+    }
+    // Will update the node's AABB at bvh_nodes[index]
+    fn update_bvh_node_aabb(&mut self, index: usize) {
+        // We will make his node bound all its primitives
+        let bvh_node = &mut self.bvh_nodes[index]; // Current BVHNode
+        let bvh_node_aabb = &mut bvh_node.aabb; //Current node AABB
+
+        let first_prim = bvh_node.first_prim; //Start index of prim
+        let prim_count = bvh_node.prim_count; //Number of primitives within the nodes aabb
+
+        for i in 0..prim_count {
+            let node = &self.nodes[first_prim + i]; //Get the node from the Vec<Node>
+            bvh_node_aabb.join_mut(&node.aabb); //Join it with the BVH node's AABB
+        }
+
+        // unsafe {
+        //     println!("UPDATE TO AABB ---- {STATIC0}");
+        //     STATIC0 += 1;
+        //     let bvh_node = &mut self.bvh_nodes[index]; //Get the BVHNode we are working on
+        //     println!("{bvh_node}");
+        // }
+    }
+    // Subdivision, will subdivide a split
+    fn subdivide(&mut self, index: usize) {
+        //Get the bvh_node we will be altering
+        // Determine the axis and position of the split plane
+        // Split the group of primitives in two halves using the split plane
+        // Create child nodes for each half
+        // Recurse into each of the child nodes.
+
+        //Leaf node case, we cannot sub-divide any more
+        if self.bvh_nodes[index].prim_count == 1 {
+            return;
+        };
+
+        /* ------------ SUBDIVIDE BY LONGEST AXIS ------------ */
+        //Get information about the node we want to subdivide
+        let (bln, trf) = (
+            self.bvh_nodes[index].aabb.bln,
+            self.bvh_nodes[index].aabb.trf,
+        );
+        let extent = trf - bln;
+        let mut axis = 0; // Assume that x is longest
+        if extent.y > extent.x {
+            axis = 1; // Split y if longest
+        };
+        if extent.z > extent[axis] {
+            axis = 2; // Split z if longest
+        };
+        let split_pos = bln[axis] + extent[axis] * 0.5; // Final split down the middle of AABB
+
+        /* --------- SUBDIVIDE BY Surface Area Heuristic ---------*/
+        // let mut best_axis: Option<usize> = None;
+        // let mut best_pos = 0.0;
+        // let mut best_cost = 1e30;
+        // let first_prim_idx = self.bvh_nodes[index].first_prim;
+        // for axis in 0..3 {
+        //     for i in 0..self.bvh_nodes[index].prim_count {
+        //         let node = &self.nodes[first_prim_idx + i];
+        //         //Get the centroid of the bounding box
+        //         let centroid = node.aabb.get_centroid();
+        //         //Get the candidate position
+        //         let candidate_pos = world_centroid[axis];
+        //         let cost = self.evaluate_sah(&self.bvh_nodes[index], axis, candidate_pos);
+        //         if cost < best_cost {
+        //             best_pos = candidate_pos;
+        //             best_axis = Some(axis);
+        //             best_cost = cost;
+        //         }
+        //     }
+        // }
+        // let axis = match best_axis {
+        //     Some(axis) => axis,
+        //     None => 0,
+        // };
+        // let split_pos = best_pos;
+
+        let left_count;
+        let right_count;
+        let mut i;
+        let mut j;
+        {
+            let bvh_node = &mut self.bvh_nodes[index];
+            i = bvh_node.first_prim; //Start of array
+            j = i + bvh_node.prim_count - 1; //End of array
+            while i <= j {
+                //Perform a quicksort dependent on location
+                let node = &self.nodes[i]; // Node we would like to sort
+                let centroid = node.aabb.get_centroid(); //Centroid of node we would like to sort
+                if centroid[axis] < split_pos {
+                    i += 1; // On Left-Hand-Side
+                } else {
+                    self.nodes.swap(i, j);
+                    j -= 1; // On Right-Hand-Side
+                }
+            }
+            //Now we have two splits
+            //The lhs of the array is in the left split  0..left_count
+            //The rhs of the array is on the right split left_count + 1..n
+            left_count = i - bvh_node.first_prim; //Number of prims on lhs
+            right_count = bvh_node.prim_count - left_count;
+            //println!("SPLIT INTO: {left_count} {right_count}");
+            if left_count == 0 || left_count == bvh_node.prim_count {
+                //Split did nothing
+                return;
+            }
+        }
+        // unsafe {
+        //     println!("SUBDIVIDE: {STATIC1}");
+        //     println!("SPLIT INTO: {left_count} ");
+        //     STATIC1 += 1;
+        // }
+
+        let l_idx = self.nodes_used; //Left child
+        self.bvh_nodes[index].l_idx = l_idx;
+        self.nodes_used = self.nodes_used + 2;
+
+        //Set left node information
+        self.bvh_nodes[l_idx].first_prim = self.bvh_nodes[index].first_prim; //Left split begins at parent split
+        self.bvh_nodes[l_idx].prim_count = left_count; // Left prims
+
+        //Set right node information
+        self.bvh_nodes[l_idx + 1].first_prim = i; // Right split start index
+        self.bvh_nodes[l_idx + 1].prim_count = right_count;
+
+        //Current node is not a leaf node
+        self.bvh_nodes[index].prim_count = 0;
+
+        self.update_bvh_node_aabb(l_idx); //Update AABB for left of split
+        self.update_bvh_node_aabb(l_idx + 1); //Update AABB for right of split
+
+        //Recurse
+        self.subdivide(l_idx); // Subdivide left index
+        self.subdivide(l_idx + 1); // SUbdivide right index
+    }
+    // Traverse the BVH, 0 will be needed to start at root node
+    pub fn traverse(&self, ray: &Ray, idx: usize) -> Option<(&Node, Intersection)> {
+        let bvh_node = &self.bvh_nodes[idx];
+        if !bvh_node.aabb.intersect_ray(&ray) {
+            // No intersection with BVH in world coordinates
+            return None;
+        }
+        if bvh_node.prim_count != 0 {
+            // Leaf node - check all primitives it contains
+            let mut closest: Option<(&Node, Intersection)> = None;
+            let mut closest_dist = f64::MAX;
+            for i in 0..bvh_node.prim_count {
+                let node = &self.nodes[bvh_node.first_prim + i];
+                if !node.active {
+                    continue;
+                }
+                if let Some(intersect) = node.intersect_ray(&ray) {
+                    if intersect.distance >= EPSILON && intersect.distance < closest_dist {
+                        closest_dist = intersect.distance;
+                        closest = Some((node, intersect));
+                    }
+                }
+            }
+            return closest;
+        } else {
+            //Recurse down the BVH
+            //Recurse down the BVH right node
+            let intersect_l = self.traverse(ray, bvh_node.l_idx);
+            let intersect_r = self.traverse(ray, bvh_node.l_idx + 1);
+
+            match (intersect_l, intersect_r) {
+                (None, None) => return None,
+                (Some(intersect), None) => return Some(intersect),
+                (None, Some(intersect)) => return Some(intersect),
+                (Some((node_l, inter_l)), Some((node_r, inter_r))) => {
+                    //Compare intersect distance
+                    let dist_l = distance(&ray.a, &inter_l.point);
+                    let dist_r = distance(&ray.a, &inter_r.point);
+                    if dist_l < dist_r {
+                        return Some((node_l, inter_l));
+                    } else {
+                        return Some((node_r, inter_r));
+                    }
+                }
+            }
+        }
+    }
+    fn evaluate_sah(&self, node: &BVHNode, axis: usize, pos: f64) -> f64 {
+        // determine triangle counts and bounds for this split candidate
+        let mut l_aabb = AABB::empty();
+        let mut r_aabb = AABB::empty();
+        let mut l_count = 0;
+        let mut r_count = 0;
+        for i in 0..node.prim_count {
+            let aabb = self.nodes[node.first_prim + i].get_world_aabb();
+            if aabb.trf[axis] < pos {
+                l_count += 1;
+                l_aabb.join_mut(&aabb);
+            } else {
+                r_count += 1;
+                r_aabb.join_mut(&aabb);
+            }
+        }
+        let cost = l_count as f64 * l_aabb.area() + r_count as f64 * r_aabb.area();
+        if cost > 0.0 { cost } else { 1e30 }
+    }
+}
+
+impl fmt::Display for BVH {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        for (i, node) in self.bvh_nodes.iter().enumerate() {
+            writeln!(f, "Node: {i}")?;
+            writeln!(f, "{node}")?;
+        }
+        write!(f, "")
+    }
+}

src/camera.rs

@@ -51,6 +51,66 @@ impl Camera {
         self.recalculate_matrix();
     }
 
+    /// Get the forward direction vector (from eye toward target)
+    pub fn forward(&self) -> Vector3<f64> {
+        (self.target - self.eye).normalize()
+    }
+
+    /// Get the right direction vector
+    pub fn right(&self) -> Vector3<f64> {
+        self.forward().cross(&self.up).normalize()
+    }
+
+    /// Move the camera forward/backward along its view direction (moves both eye and target)
+    pub fn move_forward(&mut self, amount: f64) {
+        let dir = self.forward() * amount;
+        self.eye += dir;
+        self.target += dir;
+        self.recalculate_matrix();
+    }
+
+    /// Strafe the camera left/right (moves both eye and target)
+    pub fn move_right(&mut self, amount: f64) {
+        let dir = self.right() * amount;
+        self.eye += dir;
+        self.target += dir;
+        self.recalculate_matrix();
+    }
+
+    /// Move the camera up/down along the up vector (moves both eye and target)
+    pub fn move_up(&mut self, amount: f64) {
+        let dir = self.up.normalize() * amount;
+        self.eye += dir;
+        self.target += dir;
+        self.recalculate_matrix();
+    }
+
+    /// Orbit the camera around the target point by yaw (horizontal) and pitch (vertical) angles in radians
+    pub fn orbit(&mut self, yaw: f64, pitch: f64) {
+        let offset = self.eye - self.target;
+        let radius = offset.norm();
+
+        // Current spherical angles
+        let current_pitch = (offset.y / radius).asin();
+        let current_yaw = offset.z.atan2(offset.x);
+
+        let new_yaw = current_yaw + yaw;
+        let new_pitch = (current_pitch + pitch).clamp(
+            -std::f64::consts::FRAC_PI_2 + 0.01,
+            std::f64::consts::FRAC_PI_2 - 0.01,
+        );
+
+        // Convert back to cartesian
+        let new_offset = Vector3::new(
+            radius * new_pitch.cos() * new_yaw.cos(),
+            radius * new_pitch.sin(),
+            radius * new_pitch.cos() * new_yaw.sin(),
+        );
+
+        self.eye = self.target + new_offset;
+        self.recalculate_matrix();
+    }
+
     /// Recalculate the view and inverse view matrices based on the current eye, target, and up vectors
     fn recalculate_matrix(&mut self) {
         self._view = Matrix4::look_at_lh(&self.eye, &self.target, &self.up);

src/gui.rs

@@ -5,41 +5,49 @@ use crate::{
     node::*,
     primitive::*,
     scene::*,
-    state::{INIT_FILE, SAVE_FILE},
+    state::{RaytracingOption, INIT_FILE, SAVE_FILE},
 };
 use imgui::*;
 use nalgebra::{Point3, Vector3};
 use pixels::{wgpu, PixelsContext};
 use rhai::Engine;
-use std::time::Instant;
+use std::time::{Duration, Instant};
 
 //BUFFER CONSTANTS
-const BUFFER_PROPORTION_INIT: f32 = 0.2;
 const BUFFER_PROPORTION_MIN: f32 = 0.1;
 const BUFFER_PROPORTION_MAX: f32 = 1.0;
 
 //RAY CONSTANTS
-const RAYS_INIT: i32 = 100;
-const RAYS_MIN: i32 = 100;
-const RAYS_MAX: i32 = 10000;
+const MIN_THREADS: u32 = 1;
+const MAX_THREADS: u32 = 12;
+const RAYS_MIN: u32 = 100;
+const RAYS_MAX: u32 = 10000;
+const MIN_DEPTH: u8 = 1;
+const MAX_DEPTH: u8 = 10;
+const MIN_SAMPLES: u32 = 1;
+const MAX_SAMPLES: u32 = 10;
+const MIN_RANDOM: f64 = 100.0;
+const MAX_RANDOM: f64 = 1000.0;
+const MIN_EPSILON: f64 = 1e-11;
+const MAX_EPSILON: f64 = 1.0;
+
+//DIFFUSE CONSTANTS
+const MIN_DIFFUSE_RAYS: u8 = 1;
+const MAX_DIFFUSE_RAYS: u8 = 10;
+const MIN_DIFFUSE_COEFFICIENT: f32 = 0.0;
+const MAX_DIFFUSE_COEFFICIENT: f32 = 1.0;
 
 //MATERIAL CONSTANTS
-const MIN_D: f32 = 0.0;
-const MIN_S: f32 = 0.0;
 const MIN_SHINE: f32 = 0.0;
-const MAX_D: f32 = 1.0;
-const MAX_S: f32 = 1.0;
 const MAX_SHINE: f32 = 50.0;
 
 //TRANSFORMATION CONSTANTS
-const MIN_COLOUR: f32 = 0.0;
 const MIN_FALLOFF: f32 = 0.0;
 const MIN_SCALE: f64 = 0.0;
 //const MIN_POSITION: f64 = -10.0;
 const MIN_ROTATION: f64 = -180.0;
 const MIN_TRANSLATE: f64 = -10.0;
 //--
-const MAX_COLOUR: f32 = 1.0;
 const MAX_FALLOFF: f32 = 1.0;
 const MAX_SCALE: f64 = 3.0;
 //const MAX_POSITION: f64 = 10.0;
@@ -47,14 +55,14 @@ const MAX_ROTATION: f64 = 180.0;
 const MAX_TRANSLATE: f64 = 10.0;
 
 // CAMERA CONSTANTS
-const MIN_FOV: f32 = 10.0;
-const MAX_FOV: f32 = 160.0;
+const MIN_FOV: f64 = 10.0;
+const MAX_FOV: f64 = 160.0;
 //const CAMERA_INIT: f32 = 5.0;
 
 /// Manages all state required for rendering Dear ImGui over `Pixels`test.
 pub enum GuiEvent {
-    BufferResize(f32, f32),
-    CameraUpdate(Camera, f32),
+    RaytracerOption(RaytracingOption),
+    CameraUpdate(Camera),
     SceneLoad(Scene),
     SaveImage(String),
 }
@@ -67,17 +75,17 @@ pub struct Gui {
 
     pub event: Option<GuiEvent>,
 
+    render_start: Option<Instant>,
+    render_elapsed: Option<Duration>,
+
     script_filename: String,
     script: String,
     engine: Engine,
     scene: Scene,
 
-    pub ray_num: i32,
-
-    buffer_proportion: f32,
+    raytracing_option: RaytracingOption,
 
     camera: Camera,
-    camera_fov: f32,
 
     image_filename: String,
 }
@@ -122,7 +130,7 @@ impl Gui {
         let renderer = imgui_wgpu::Renderer::new(&mut imgui, device, queue, config);
 
         // Return GUI context
-        Self {
+        let mut gui = Self {
             imgui,
             platform,
             renderer,
@@ -130,18 +138,47 @@ impl Gui {
             last_cursor: None,
             event: None,
 
+            render_start: None,
+            render_elapsed: None,
+
             script_filename: String::from(INIT_FILE),
             script: String::new(),
             engine: init_engine(),
             scene: Scene::empty(),
 
-            ray_num: RAYS_INIT,
-            buffer_proportion: BUFFER_PROPORTION_INIT,
+            raytracing_option: RaytracingOption::default(),
 
             camera: Camera::unit(),
-            camera_fov: 110.0,
 
             image_filename: String::from(SAVE_FILE),
+        };
+
+        // ------------ TESTING CODE (LOAD SCENE ON START) -----------------
+        match std::fs::read_to_string(&mut gui.script_filename) {
+            Ok(script) => {
+                gui.script = script;
+            }
+            Err(e) => println!("{}", e),
+        }
+        match gui.engine.eval(&gui.script) {
+            Ok(scene) => {
+                gui.scene = scene;
+                gui.event = Some(GuiEvent::SceneLoad(gui.scene.clone()));
+            }
+            Err(e) => println!("{e}"),
+        }
+        // ------------ TESTING CODE (LOAD SCENE ON START) -----------------
+        gui
+    }
+
+    pub fn start_render_timer(&mut self) {
+        self.render_start = Some(Instant::now());
+        self.render_elapsed = None;
+    }
+
+    pub fn stop_render_timer(&mut self) {
+        if let Some(start) = self.render_start.take() {
+            self.render_elapsed = Some(start.elapsed());
         }
     }
 
@@ -183,38 +220,131 @@ impl Gui {
 
         //Raytracing options -------------------------------------------
         if CollapsingHeader::new("Raytracer").build(ui) {
-            // Numbers of rays to render
-            ui.slider("# Rays: ", RAYS_MIN, RAYS_MAX, &mut self.ray_num);
-            // Proportion of the window the buffer occupies
             ui.slider(
-                "% Buffer: ",
-                BUFFER_PROPORTION_MIN,
-                BUFFER_PROPORTION_MAX,
-                &mut self.buffer_proportion,
+                "Threads",
+                MIN_THREADS,
+                MAX_THREADS,
+                &mut self.raytracing_option.threads,
             );
+            // Numbers of rays to render per pass
+            Drag::new("Rays Per Pass")
+                .range(RAYS_MIN, RAYS_MAX)
+                .speed(50.0)
+                .build(ui, &mut self.raytracing_option.pixels_per_thread);
+            // Proportion of the window the buffer occupies
+            Drag::new("% Buffer: ")
+                .range(BUFFER_PROPORTION_MIN, BUFFER_PROPORTION_MAX)
+                .speed(0.005)
+                .display_format("%.2f")
+                .build(ui, &mut self.raytracing_option.buffer_proportion);
+            //Clear colour for scene
+            let mut clear_f32 = [
+                self.raytracing_option.clear_color[0] as f32 / 255.0,
+                self.raytracing_option.clear_color[1] as f32 / 255.0,
+                self.raytracing_option.clear_color[2] as f32 / 255.0,
+                self.raytracing_option.clear_color[3] as f32 / 255.0,
+            ];
+            if ui.color_edit4_config("Clear Colour", &mut clear_f32).alpha_bar(true).build() {
+                self.raytracing_option.clear_color = [
+                    (clear_f32[0] * 255.0) as u8, (clear_f32[1] * 255.0) as u8,
+                    (clear_f32[2] * 255.0) as u8, (clear_f32[3] * 255.0) as u8,
+                ];
+            }
+            //Clear colour if no intersect
+            let mut pixel_clear_f32 = [
+                self.raytracing_option.pixel_clear[0] as f32 / 255.0,
+                self.raytracing_option.pixel_clear[1] as f32 / 255.0,
+                self.raytracing_option.pixel_clear[2] as f32 / 255.0,
+                self.raytracing_option.pixel_clear[3] as f32 / 255.0,
+            ];
+            if ui.color_edit4_config("Pixel Clear Colour", &mut pixel_clear_f32).alpha_bar(true).build() {
+                self.raytracing_option.pixel_clear = [
+                    (pixel_clear_f32[0] * 255.0) as u8, (pixel_clear_f32[1] * 255.0) as u8,
+                    (pixel_clear_f32[2] * 255.0) as u8, (pixel_clear_f32[3] * 255.0) as u8,
+                ];
+            }
+            //Ray depth slider
+            ui.slider(
+                "Ray Depth",
+                MIN_DEPTH,
+                MAX_DEPTH,
+                &mut self.raytracing_option.ray_depth,
+            );
+            //Ray samples slider
+            ui.slider(
+                "Ray Samples",
+                MIN_SAMPLES,
+                MAX_SAMPLES,
+                &mut self.raytracing_option.ray_samples,
+            );
+            //Ray randomness
+            Drag::new("Ray Randomness")
+                .range(MIN_RANDOM, MAX_RANDOM)
+                .speed(5.0)
+                .display_format("%.1f")
+                .build(ui, &mut self.raytracing_option.ray_randomness);
+            //Number of diffuse rays
+            ui.slider(
+                "Diffuse Rays",
+                MIN_DIFFUSE_RAYS,
+                MAX_DIFFUSE_RAYS,
+                &mut self.raytracing_option.diffuse_rays,
+            );
+            //Diffuse Coefficient
+            Drag::new("Diffuse Coefficient")
+                .range(MIN_DIFFUSE_COEFFICIENT, MAX_DIFFUSE_COEFFICIENT)
+                .speed(0.005)
+                .display_format("%.3f")
+                .build(ui, &mut self.raytracing_option.diffuse_coefficient);
             // Fov of the buffer
-            ui.slider("fov", MIN_FOV, MAX_FOV, &mut self.camera_fov);
+            ui.slider(
+                "fov",
+                MIN_FOV,
+                MAX_FOV,
+                &mut self.raytracing_option.buffer_fov,
+            );
+            // Enable BVH
+            ui.checkbox("Enable BVH", &mut self.raytracing_option.bvh_active);
+            ui.checkbox("Enable Shadows", &mut self.raytracing_option.shadows);
+            ui.checkbox("Enable Reflections", &mut self.raytracing_option.reflect);
+            ui.checkbox("Enable Specular", &mut self.raytracing_option.specular);
+            ui.checkbox("Enable Diffuse", &mut self.raytracing_option.diffuse);
+            // Render timer display
+            ui.separator();
+            if let Some(start) = &self.render_start {
+                let elapsed = start.elapsed().as_secs_f64();
+                ui.text(format!("Rendering: {:.2}s", elapsed));
+            } else if let Some(elapsed) = &self.render_elapsed {
+                ui.text(format!("Render time: {:.2}s", elapsed.as_secs_f64()));
+            }
+            ui.separator();
             // Apply stored changes
             if ui.button("Apply") {
-                self.event = Some(GuiEvent::BufferResize(
-                    self.buffer_proportion,
-                    self.camera_fov,
-                ));
+                self.event = Some(GuiEvent::RaytracerOption(self.raytracing_option.clone()));
             };
         }
         // CAMERA OPTIONS ----------------------------------------
         if CollapsingHeader::new("Camera").build(ui) {
             // Eye, target and up vector inputs
             ui.text("Camera options:");
-            ui.slider_config("Eye", MIN_TRANSLATE, MAX_TRANSLATE)
-                .build_array(self.camera.eye.coords.as_mut_slice());
-            ui.slider_config("Target", MIN_TRANSLATE, MAX_TRANSLATE)
-                .build_array(self.camera.target.coords.as_mut_slice());
-            ui.slider_config("Up", 0.0, 1.0)
-                .build_array(self.camera.up.as_mut_slice());
+            Drag::new("Eye")
+                .range(MIN_TRANSLATE, MAX_TRANSLATE)
+                .speed(0.05)
+                .display_format("%.2f")
+                .build_array(ui, self.camera.eye.coords.as_mut_slice());
+            Drag::new("Target")
+                .range(MIN_TRANSLATE, MAX_TRANSLATE)
+                .speed(0.05)
+                .display_format("%.2f")
+                .build_array(ui, self.camera.target.coords.as_mut_slice());
+            Drag::new("Up")
+                .range(0.0, 1.0)
+                .speed(0.005)
+                .display_format("%.3f")
+                .build_array(ui, self.camera.up.as_mut_slice());
             if ui.button("Apply Camera") {
                 println!("Camera changed");
-                self.event = Some(GuiEvent::CameraUpdate(self.camera.clone(), self.camera_fov));
+                self.event = Some(GuiEvent::CameraUpdate(self.camera.clone()));
             }
         }
         // SCRIPTING --------------------------------------------
@@ -266,9 +396,7 @@ impl Gui {
         // SCENE --------------------------------------------
         if CollapsingHeader::new("Scene").build(ui) {
             if ui.button("Update Scene") {
-                for (_, node) in &mut self.scene.nodes {
-                    node.compute();
-                }
+                self.scene.compute();
                 self.event = Some(GuiEvent::SceneLoad(self.scene.clone()));
             }
             // Edit transformation of nodes
@@ -277,12 +405,21 @@ impl Gui {
                     ui.checkbox(format!("##active{label}"), &mut node.active);
                     ui.same_line();
                     if let Some(_t) = ui.tree_node(label) {
-                        ui.slider_config("Translation", MIN_TRANSLATE, MAX_TRANSLATE)
-                            .build_array(&mut node.translation);
-                        ui.slider_config("Rotation", MIN_ROTATION, MAX_ROTATION)
-                            .build_array(&mut node.rotation);
-                        ui.slider_config("Scale", MIN_SCALE, MAX_SCALE)
-                            .build_array(&mut node.scale);
+                        Drag::new("Translation")
+                            .range(MIN_TRANSLATE, MAX_TRANSLATE)
+                            .speed(0.05)
+                            .display_format("%.2f")
+                            .build_array(ui, &mut node.translation);
+                        Drag::new("Rotation")
+                            .range(MIN_ROTATION, MAX_ROTATION)
+                            .speed(1.0)
+                            .display_format("%.1f")
+                            .build_array(ui, &mut node.rotation);
+                        Drag::new("Scale")
+                            .range(MIN_SCALE, MAX_SCALE)
+                            .speed(0.01)
+                            .display_format("%.3f")
+                            .build_array(ui, &mut node.scale);
                     }
                 }
             }
@@ -290,11 +427,19 @@ impl Gui {
             if let Some(_t) = ui.tree_node("Materials") {
                 for (label, material) in &mut self.scene.materials {
                     if let Some(_t) = ui.tree_node(label) {
-                        ui.slider_config("ks", MIN_D, MAX_D)
-                            .build_array(material.ks.as_mut_slice());
-                        ui.slider_config("kd", MIN_S, MAX_S)
-                            .build_array(material.kd.as_mut_slice());
-                        ui.slider("shine", MIN_SHINE, MAX_SHINE, &mut material.shininess);
+                        let mut ks_arr: [f32; 3] = material.ks.into();
+                        if ui.color_edit3("ks", &mut ks_arr) {
+                            material.ks = Vector3::from(ks_arr);
+                        }
+                        let mut kd_arr: [f32; 3] = material.kd.into();
+                        if ui.color_edit3("kd", &mut kd_arr) {
+                            material.kd = Vector3::from(kd_arr);
+                        }
+                        Drag::new("shine")
+                            .range(MIN_SHINE, MAX_SHINE)
+                            .speed(0.5)
+                            .display_format("%.1f")
+                            .build(ui, &mut material.shininess);
                     }
                 }
             }
@@ -304,12 +449,20 @@ impl Gui {
                     ui.checkbox(format!("##activelight{label}"), &mut light.active);
                     ui.same_line();
                     if let Some(_t) = ui.tree_node(label) {
-                        ui.slider_config("Colour", MIN_COLOUR, MAX_COLOUR)
-                            .build_array(light.colour.as_mut_slice());
-                        ui.slider_config("Position", MIN_TRANSLATE, MAX_TRANSLATE)
-                            .build_array(light.position.coords.as_mut_slice());
-                        ui.slider_config("Falloff", MIN_FALLOFF, MAX_FALLOFF)
-                            .build_array(light.falloff.as_mut_slice());
+                        let mut colour_arr: [f32; 3] = light.colour.into();
+                        if ui.color_edit3("Colour", &mut colour_arr) {
+                            light.colour = Vector3::from(colour_arr);
+                        }
+                        Drag::new("Position")
+                            .range(MIN_TRANSLATE, MAX_TRANSLATE)
+                            .speed(0.05)
+                            .display_format("%.2f")
+                            .build_array(ui, light.position.coords.as_mut_slice());
+                        Drag::new("Falloff")
+                            .range(MIN_FALLOFF, MAX_FALLOFF)
+                            .speed(0.005)
+                            .display_format("%.3f")
+                            .build_array(ui, light.falloff.as_mut_slice());
                     }
                 }
             }
@@ -318,7 +471,7 @@ impl Gui {
                 for (label, camera) in &self.scene.cameras {
                     if ui.button(label) {
                         self.camera = camera.clone();
-                        self.event = Some(GuiEvent::CameraUpdate(camera.clone(), self.camera_fov));
+                        self.event = Some(GuiEvent::CameraUpdate(camera.clone()));
                     }
                 }
             }
@@ -346,6 +499,11 @@ impl Gui {
         )
     }
 
+    /// Update the GUI's camera to reflect external changes (e.g. from keyboard/mouse movement)
+    pub fn update_camera(&mut self, camera: &Camera) {
+        self.camera = camera.clone();
+    }
+
     /// Handle any outstanding events.
     pub fn handle_event(
         &mut self,
@@ -450,8 +608,8 @@ pub fn init_engine() -> Engine {
         .register_type::<Mesh>()
         .register_fn("Mesh", Mesh::from_file);
     engine
-        .register_type::<Rectangle>()
-        .register_fn("Rectange", Rectangle::new)
-        .register_fn("RectangleUnit", Rectangle::unit);
+        .register_type::<RectangleXY>()
+        .register_fn("Rectange", RectangleXY::new)
+        .register_fn("RectangleUnit", RectangleXY::unit);
     engine
 }

src/main.rs

@@ -1,10 +1,11 @@
 use crate::state::run;
 use error_iter::ErrorIter;
 
-const EPSILON: f64 = 1e-8;
+const EPSILON: f64 = 1e-7;
 const INFINITY: f64 = 1e10;
 
 use log::error;
+//use nalgebra::{Matrix4, RowVector4, Vector3, Vector4};
 use std::env;
 use std::error::Error;
 
@@ -21,13 +22,159 @@ mod state;
 
 fn main() {
     env_logger::init();
-
     env::set_var("RUST_BACKTRACE", "1");
+    //let args: Vec<String> = env::args().collect();
+
+    // let vec = Vector3::new(1.0, 1.0, 1.0);
+    // let translation = Vector3::new(1.0, 1.0, 1.0);
+    // let translation_matrix = Matrix4::new_translation(&translation);
+    // println!(
+    //     "{}, {}",
+    //     translation_matrix,
+    //     translation_matrix.transform_vector(&vec)
+    // );
+    // let mut translation_matrix = translation_matrix.transpose();
+    // translation_matrix.set_row(3, &RowVector4::new(0.0, 0.0, 0.0, 0.0));
+    // println!(
+    //     "{}, {}", //     translation_matrix, //     translation_matrix.transform_vector(&vec)
+    // );
     if let Err(e) = run() {
         println!("Error at runtime: {}", e);
     };
+
+    // if args.len() == 6 {
+    //     let width: usize = args[1].parse().unwrap();
+    //     let height: usize = args[2].parse().unwrap();
+    //     let fovy = args[3].parse::<f64>().unwrap();
+    //     let filename = &args[4];
+    //     let savefile = &args[5];
+    //     headless(
+    //         width,
+    //         height,
+    //         fovy,
+    //         filename.to_string(),
+    //         savefile.to_string(),
+    //     );
+    // } else {
+    //}
 }
 
+// fn headless(width: usize, height: usize, fovy: f64, filename: String, savefile: String) {
+//     let options = Arc::new(RaytracingOption {
+//         threads: 12,
+//         ray_samples: 1,
+//         ray_randomness: 100.0,
+//         clear_color: [0x22, 0x00, 0x11, 0x55],
+//         pixel_clear: [0x55, 0x00, 0x22, 0x55],
+//         pixels_per_thread: 200,
+//         buffer_proportion: 1.0,
+//         buffer_fov: 110.0,
+//         ray_depth: 5,
+//         diffuse_rays: 3,
+//         diffuse_coefficient: 0.8,
+//         bvh_active: false,
+//     });
+//     //Read script from file
+//     let script = match std::fs::read_to_string(&filename) {
+//         Ok(in_script) => in_script,
+//         Err(e) => {
+//             println!("{}", e);
+//             return;
+//         }
+//     };
+//     //Evaluate scene in file
+//     let engine = init_engine();
+//     let scene: Arc<Scene> = match engine.eval(&script) {
+//         Ok(in_scene) => Arc::new(in_scene),
+//         Err(e) => {
+//             println!("{e}");
+//             return;
+//         }
+//     };
+//     //Set the camera
+//     let mut camera = Camera::unit();
+//     for (_, in_camera) in &scene.cameras {
+//         camera = in_camera.clone();
+//     }
+//     //Cast the rays
+//     let rays = Arc::new(Ray::cast_rays(
+//         &camera.eye,
+//         &camera.target,
+//         &camera.up,
+//         fovy,
+//         width as u32,
+//         height as u32,
+//     ));
+//     //Enable bounding volume heirarchy
+//     let bvh;
+//     match options.bvh_active {
+//         true => bvh = Arc::new(Some(BVH::build(&scene.nodes))),
+//         false => bvh = Arc::new(None),
+//     }
+//     //Create our frame and indexer
+//     let size = width * height;
+//     let frame_mutex = Arc::new(Mutex::new(vec![0; size * 4]));
+//     //Multithreading
+//     let mut handles = vec![];
+
+//     for index in 0..size {
+//         for _ in 0..options.threads {
+//             //Get random index from queue
+//             //Create a nre thread for this pixel
+//             let handle = thread::spawn({
+//                 let rays = rays.clone();
+//                 let scene = scene.clone();
+//                 let options = options.clone();
+//                 let bvh = bvh.clone();
+//                 let rays = rays.clone();
+//                 let frame_mutex = frame_mutex.clone();
+//                 move || {
+//                     //Shade colour for selected ray
+//                     let mut colour: Vector3<f32> = Vector3::zeros();
+//                     //Get the ray we want to make
+//                     let shot_ray = &rays[index];
+//                     //Send out ray_samples rays
+//                     for _ in 0..options.ray_samples {
+//                         let point = shot_ray.a;
+//                         let dir = shot_ray.b;
+//                         //Generate a random ray
+//                         let rx = (random::<f64>() - 0.5) / options.ray_randomness;
+//                         let ry = (random::<f64>() - 0.5) / options.ray_randomness;
+//                         let rz = (random::<f64>() - 0.5) / options.ray_randomness;
+//                         let nx = dir.x + rx;
+//                         let ny = dir.y + ry;
+//                         let nz = dir.z + rz;
+//                         let rand_ray = Ray::new(point, Vector3::new(nx, ny, nz));
+
+//                         if let Some(ray_colour) = rand_ray.shade_ray(&scene, 0, &options, &bvh) {
+//                             colour += ray_colour;
+//                         }
+//                     }
+//                     colour = (colour / options.ray_samples as f32) * 255.0;
+//                     let rgba = [colour.x as u8, colour.y as u8, colour.z as u8, 0xff];
+//                     {
+//                         let frame = &mut frame_mutex.lock().unwrap();
+//                         frame[index * 4..(index + 1) * 4].copy_from_slice(&rgba);
+//                     }
+//                 }
+//             });
+//             handles.push(handle);
+//         }
+//         for handle in handles.drain(..) {
+//             handle.join().unwrap();
+//         }
+//     }
+//     use std::path::Path;
+//     image::save_buffer(
+//         Path::new(&savefile),
+//         &frame_mutex.lock().unwrap(),
+//         width as u32,
+//         height as u32,
+//         image::ColorType::Rgba8,
+//     )
+//     .unwrap();
+// }
+
 fn log_error<E: Error + 'static>(method_name: &str, err: E) {
     error!("{method_name}() failed: {err}");
     for source in err.sources().skip(1) {

src/material.rs

@@ -10,10 +10,10 @@ pub struct Material {
 }
 
 impl Material {
-    pub fn new(kd: Vector3<f64>, ks: Vector3<f64>, shininess: f64) -> Material {
+    pub fn new(kd: Vector3<f64>, ks: Vector3<f64>, kr: Vector3<f64>, shininess: f64) -> Material {
         let kd = kd.cast();
         let ks = ks.cast();
-        let kr = ks.cast();
+        let kr = kr.cast();
         let shininess = shininess as f32;
         Material {
             kd,

src/node.rs

@@ -1,12 +1,19 @@
-use crate::{material::Material, primitive::*};
-use nalgebra::{Matrix4, Vector3};
-use std::rc::Rc;
+use crate::{
+    bvh::AABB,
+    material::Material,
+    primitive::*,
+    ray::{Intersection, Ray},
+    EPSILON,
+};
+use nalgebra::{distance, Matrix3, Matrix4, Vector3};
+use std::sync::Arc;
 
 #[derive(Clone)]
 pub struct Node {
     //Primitive
-    pub primitive: Rc<dyn Primitive>,
+    pub primitive: Arc<dyn Primitive>,
     pub material: Material,
+    pub aabb: AABB,
     //Transformations
     pub rotation: [f64; 3],
     pub scale: [f64; 3],
@@ -14,27 +21,30 @@ pub struct Node {
     //Model matricies
     pub model: Matrix4<f64>,
     pub inv_model: Matrix4<f64>,
-
+    pub inv_transpose_model: Matrix3<f64>,
+    //If the node is active
     pub active: bool,
 }
 
 impl Node {
     //New node with no transformations
-    pub fn new(primitive: Rc<dyn Primitive>, material: Material) -> Node {
+    pub fn new(primitive: Arc<dyn Primitive>, material: Material) -> Node {
+        let aabb = primitive.get_aabb();
         Node {
             primitive,
             material,
+            aabb,
             rotation: [0.0, 0.0, 0.0],
             scale: [1.0, 1.0, 1.0],
             translation: [0.0, 0.0, 0.0],
             model: Matrix4::identity(),
             inv_model: Matrix4::identity(),
-
+            inv_transpose_model: Matrix3::identity(),
             active: true,
         }
     }
     //New node with parent transformations
-    pub fn child(self, primitive: Rc<dyn Primitive>) -> Node {
+    pub fn child(self, primitive: Arc<dyn Primitive>) -> Node {
         let mut child = self.clone();
         child.primitive = primitive;
         child
@@ -46,12 +56,6 @@ impl Node {
 
     //Rotate a mesh by adding to its rotation
     pub fn rotate(&mut self, roll: f64, pitch: f64, yaw: f64) {
-        //Convert to radians
-        let roll = roll.to_radians();
-        // Convert pitch and yaw to radians
-        let pitch = pitch.to_radians();
-        let yaw = yaw.to_radians();
-
         // Add the roll, pitch, and yaw to the current rotation
         self.rotation[0] += roll;
         self.rotation[1] += pitch;
@@ -71,9 +75,9 @@ impl Node {
     }
     // Scale a mesh by adding to its current scale
     pub fn scale(&mut self, x: f64, y: f64, z: f64) {
-        self.scale[0] += x;
-        self.scale[1] += y;
-        self.scale[2] += z;
+        self.scale[0] = x;
+        self.scale[1] = y;
+        self.scale[2] = z;
 
         // Recompute the model and inverse model matrices
         self.compute();
@@ -88,10 +92,32 @@ impl Node {
         let scale_matrix = Matrix4::new_nonuniform_scaling(&scale);
         // Rotation matrix
         let (roll, pitch, yaw) = (self.rotation[0], self.rotation[1], self.rotation[2]);
-        let rotation_matrix = Matrix4::from_euler_angles(roll, pitch, yaw);
+        let rotation_matrix =
+            Matrix4::from_euler_angles(roll.to_radians(), pitch.to_radians(), yaw.to_radians());
         // Compute the model matrix by combining the translation, rotation, and scale matrices
         self.model = (translation_matrix * rotation_matrix * scale_matrix).cast();
         // Compute the inverse model matrix by inverting the model matrix
         self.inv_model = self.model.try_inverse().unwrap();
+        self.inv_transpose_model = self.inv_model.transpose().remove_row(3).remove_column(3);
+        // Reset AABB from primitive local space before transforming to world space
+        self.aabb = self.primitive.get_aabb();
+        self.aabb.transform_mut(&self.model);
+    }
+    // Intersection of a ray, will convert to model coords and check
+    pub fn intersect_ray(&self, ray: &Ray) -> Option<Intersection> {
+        let local_ray = ray.transform(&self.inv_model); //Transform from world coordinates
+        if let Some(mut intersect) = self.primitive.intersect_ray(&local_ray) {
+            if intersect.distance < EPSILON {
+                return None;
+            }
+            intersect.transform_mut(&self.model, &self.inv_transpose_model); //Transform to world coords
+            intersect.distance = distance(&intersect.point, &ray.a); // use world-space ray origin
+            return Some(intersect);
+        }
+        return None;
+    }
+    //Gets the bounding box in world coords
+    pub fn get_world_aabb(&self) -> AABB {
+        return self.aabb.clone();
     }
 }

src/primitive.rs

@@ -9,11 +9,11 @@ use nalgebra::{distance, Point3, Vector3};
 use roots::{find_roots_quadratic, find_roots_quartic, Roots};
 use std::fs::File;
 use std::io::{BufRead, BufReader};
-use std::rc::Rc;
+use std::sync::Arc;
 // PRIMITIVE TRAIT -----------------------------------------------------------------
-pub trait Primitive {
+pub trait Primitive: Send + Sync {
     fn intersect_ray(&self, ray: &Ray) -> Option<Intersection>;
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool;
+    fn get_aabb(&self) -> AABB;
 }
 
 // SPHERE -----------------------------------------------------------------
@@ -21,23 +21,14 @@ pub trait Primitive {
 pub struct Sphere {
     position: Point3<f64>,
     radius: f64,
-    bounding_box: AABB,
 }
 
 impl Sphere {
-    pub fn new(position: Point3<f64>, radius: f64) -> Rc<dyn Primitive> {
-        let radius_vec = Vector3::new(radius, radius, radius);
-        let bln = position - radius_vec;
-        let trf = position + radius_vec;
-        let bounding_box = AABB::new(bln, trf);
-        Rc::new(Sphere {
-            position,
-            radius,
-            bounding_box,
-        })
+    pub fn new(position: Point3<f64>, radius: f64) -> Arc<dyn Primitive> {
+        Arc::new(Sphere { position, radius })
     }
 
-    pub fn unit() -> Rc<dyn Primitive> {
+    pub fn unit() -> Arc<dyn Primitive> {
         Sphere::new(Point3::new(0.0, 0.0, 0.0), 1.0)
     }
 }
@@ -56,14 +47,13 @@ impl Primitive for Sphere {
             Roots::No(_) => return None,
             Roots::One([x1]) => x1,
             Roots::Two([x1, x2]) => {
+                // roots are returned in ascending order: x1 <= x2
                 if x1 <= 0.0 && x2 <= 0.0 {
                     return None;
+                } else if x1 <= 0.0 {
+                    x2
                 } else {
-                    if x1.abs() < x2.abs() {
-                        x1
-                    } else {
-                        x2
-                    }
+                    x1
                 }
             }
             _ => return None,
@@ -78,8 +68,12 @@ impl Primitive for Sphere {
         })
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        return self.bounding_box.intersect_bounding_box(ray);
+    fn get_aabb(&self) -> AABB {
+        let radius = self.radius;
+        let radius_vec = Vector3::new(radius, radius, radius);
+        let bln = self.position - radius_vec;
+        let trf = self.position + radius_vec;
+        AABB::new(bln, trf)
     }
 }
 
@@ -90,28 +84,32 @@ pub struct Circle {
     radius: f64,
     normal: Vector3<f64>,
     constant: f64,
-    bounding_box: AABB,
 }
 
 impl Circle {
-    pub fn new(position: Point3<f64>, radius: f64, normal: Vector3<f64>) -> Rc<dyn Primitive> {
-        let radius_vec = Vector3::new(radius, radius, radius);
-        let bln = position - radius_vec;
-        let trf = position + radius_vec;
-        let bounding_box = AABB::new(bln, trf);
-
+    pub fn new(position: Point3<f64>, radius: f64, normal: Vector3<f64>) -> Arc<dyn Primitive> {
         let normal = normal.normalize();
         let constant = normal.dot(&position.coords);
-        Rc::new(Circle {
+        Arc::new(Circle {
             position,
             radius,
             normal,
             constant,
-            bounding_box,
         })
     }
 
-    pub fn unit() -> Rc<dyn Primitive> {
+    pub fn new_unboxed(position: Point3<f64>, radius: f64, normal: Vector3<f64>) -> Circle {
+        let normal = normal.normalize();
+        let constant = normal.dot(&position.coords);
+        Circle {
+            position,
+            radius,
+            normal,
+            constant,
+        }
+    }
+
+    pub fn unit() -> Arc<dyn Primitive> {
         let position = Point3::new(0.0, 0.0, 0.0);
         let normal = Vector3::new(0.0, 0.0, -1.0);
         let radius = 1.0;
@@ -125,9 +123,9 @@ impl Primitive for Circle {
         let n_dot_b = ray.b.dot(&self.normal);
         let t = (self.constant - n_dot_a) / n_dot_b;
 
-        if t > INFINITY {
+        if t <= 0.0 || t > INFINITY {
             return None;
-        };
+        }
 
         let intersect = ray.at_t(t);
         //Distance to center of circle
@@ -144,8 +142,13 @@ impl Primitive for Circle {
         }
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let radius = self.radius;
+        let position = self.position;
+        let radius_vec = Vector3::new(radius, radius, radius);
+        let bln = position - radius_vec;
+        let trf = position + radius_vec;
+        AABB::new(bln, trf)
     }
 }
 
@@ -154,31 +157,27 @@ impl Primitive for Circle {
 pub struct Cylinder {
     radius: f64,
     height: f64,
-    base_circle: Rc<dyn Primitive>,
-    top_circle: Rc<dyn Primitive>,
-    bounding_box: AABB,
+    base_circle: Circle,
+    top_circle: Circle,
 }
 
 impl Cylinder {
-    pub fn new(radius: f64, height: f64) -> Rc<dyn Primitive> {
-        let base_circle = Circle::new(
+    pub fn new(radius: f64, height: f64) -> Arc<dyn Primitive> {
+        let base_circle = Circle::new_unboxed(
             Point3::new(0.0, 0.0, 0.0),
             radius,
             Vector3::new(0.0, -1.0, 0.0),
         );
-        let top_circle = Circle::new(
+        let top_circle = Circle::new_unboxed(
             Point3::new(0.0, height, 0.0),
             radius,
             Vector3::new(0.0, 1.0, 0.0),
         );
-        let bln = Point3::new(-radius, 0.0, -radius);
-        let trf = Point3::new(radius, height, radius);
-        Rc::new(Cylinder {
+        Arc::new(Cylinder {
             radius,
             height,
             base_circle,
             top_circle,
-            bounding_box: AABB { bln, trf },
         })
     }
 }
@@ -197,14 +196,13 @@ impl Primitive for Cylinder {
             Roots::No(_) => return None,
             Roots::One([x1]) => Some(x1),
             Roots::Two([x1, x2]) => {
+                // roots are returned in ascending order: x1 <= x2
                 if x1 <= 0.0 && x2 <= 0.0 {
                     return None;
+                } else if x1 <= 0.0 {
+                    Some(x2)
                 } else {
-                    if x1.abs() < x2.abs() {
-                        Some(x1)
-                    } else {
-                        Some(x2)
-                    }
+                    Some(x1)
                 }
             }
             _ => return None,
@@ -261,8 +259,12 @@ impl Primitive for Cylinder {
         }
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let radius = self.radius;
+        let height = self.height;
+        let bln = Point3::new(-radius, 0.0, -radius);
+        let trf = Point3::new(radius, height, radius);
+        AABB::new(bln, trf)
     }
 }
 
@@ -271,28 +273,24 @@ impl Primitive for Cylinder {
 pub struct Cone {
     height: f64,
     constant: f64,
-    circle: Rc<dyn Primitive>,
-    bounding_box: AABB,
+    circle: Circle,
 }
 
 impl Cone {
-    pub fn new(radius: f64, height: f64) -> Rc<dyn Primitive> {
-        let circle = Circle::new(
+    pub fn new(radius: f64, height: f64) -> Arc<dyn Primitive> {
+        let circle = Circle::new_unboxed(
             Point3::new(0.0, 0.0, 0.0),
             radius,
             Vector3::new(0.0, -1.0, 0.0),
         );
-        let bln = Point3::new(-radius, 0.0, -radius);
-        let trf = Point3::new(radius, height, radius);
         let constant = radius * radius / (height * height);
-        Rc::new(Cone {
+        Arc::new(Cone {
             height,
             constant,
             circle,
-            bounding_box: AABB { bln, trf },
         })
     }
-    pub fn unit() -> Rc<dyn Primitive> {
+    pub fn unit() -> Arc<dyn Primitive> {
         Cone::new(0.5, 1.0)
     }
 
@@ -325,14 +323,13 @@ impl Primitive for Cone {
             Roots::No(_) => None,
             Roots::One([x1]) => Some(x1),
             Roots::Two([x1, x2]) => {
+                // roots are returned in ascending order: x1 <= x2
                 if x1 <= 0.0 && x2 <= 0.0 {
                     None
+                } else if x1 <= 0.0 {
+                    Some(x2)
                 } else {
-                    if x1.abs() < x2.abs() {
-                        Some(x1)
-                    } else {
-                        Some(x2)
-                    }
+                    Some(x1)
                 }
             }
             _ => None,
@@ -359,90 +356,71 @@ impl Primitive for Cone {
             (None, None) => None,
             (Some(cone_intersect), None) => Some(cone_intersect),
             (None, Some(circle_intersect)) => Some(circle_intersect),
-            (Some(cone_intersect), Some(_)) => Some(cone_intersect),
+            (Some(cone_intersect), Some(circle_intersect)) => {
+                let cone_dist = distance(&ray.a, &cone_intersect.point);
+                let circle_dist = distance(&ray.a, &circle_intersect.point);
+                if cone_dist < circle_dist {
+                    Some(cone_intersect)
+                } else {
+                    Some(circle_intersect)
+                }
+            }
         }
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let height = self.height;
+        let radius = (self.constant * height * height).sqrt();
+        let bln = Point3::new(-radius, 0.0, -radius);
+        let trf = Point3::new(radius, height, radius);
+        AABB::new(bln, trf)
     }
 }
 
 // RECTANGLE -----------------------------------------------------------------
+// Normal is (0.0, 0.0, 1.0) always facing towards camera at positive z axis
 #[derive(Clone)]
-pub struct Rectangle {
-    position: Point3<f64>,
-    normal: Vector3<f64>,
-    width_direction: Vector3<f64>,
-    width: f64,
-    height: f64,
-    bounding_box: AABB,
+pub struct RectangleXY {
+    bl: Point3<f64>,
+    tr: Point3<f64>,
 }
 
-impl Rectangle {
-    pub fn new(
-        position: Point3<f64>,
-        normal: Vector3<f64>,
-        width_direction: Vector3<f64>,
-        width: f64,
-        height: f64,
-    ) -> Rc<dyn Primitive> {
-        let normal = normal.normalize();
-        let width_direction = width_direction.normalize();
-        let height_direction = width_direction.cross(&normal);
-        let bln = position - width / 2.0 * width_direction - height / 2.0 * height_direction;
-        let trf = position + width / 2.0 * width_direction + height / 2.0 * height_direction;
-        Rc::new(Rectangle {
-            position,
-            normal: normal.normalize(),
-            width_direction: width_direction.normalize(),
-            width,
-            height,
-            bounding_box: AABB { bln, trf },
-        })
+impl RectangleXY {
+    pub fn new(bl: Point3<f64>, tr: Point3<f64>) -> Arc<dyn Primitive> {
+        Arc::new(RectangleXY { bl, tr })
     }
-    pub fn unit() -> Rc<dyn Primitive> {
-        Rectangle::new(
-            Point3::new(0.0, 0.0, 0.0),
-            Vector3::new(0.0, 1.0, 0.0),
-            Vector3::new(1.0, 0.0, 0.0),
-            2.0,
-            2.0,
-        )
+    pub fn unit() -> Arc<dyn Primitive> {
+        RectangleXY::new(Point3::new(-1.0, -1.0, 0.0), Point3::new(1.0, 1.0, 0.0))
     }
 }
 
-impl Primitive for Rectangle {
+impl Primitive for RectangleXY {
     fn intersect_ray(&self, ray: &Ray) -> Option<Intersection> {
-        let constant = self.position.coords.dot(&self.normal);
-        let denominator = ray.b.dot(&self.normal);
-        let t = (constant - ray.a.coords.dot(&self.normal)) / denominator;
+        let z = self.bl.z;
+        let az = ray.a.z;
+        let bz = ray.b.z;
+        let t = (z - az) / bz;
+        if t <= 0.0 || t > INFINITY {
+            return None;
+        }
+        let intersect = ray.at_t(t);
+        let (ix, iy) = (intersect.x, intersect.y);
 
-        if t > INFINITY {
+        if (ix < self.bl.x) || (ix > self.tr.x) || (iy < self.bl.y) || (iy > self.tr.y) {
             return None;
         }
 
-        let intersect = ray.at_t(t);
-        let height_direction = self.width_direction.cross(&self.normal);
-        let (w2, h2) = (self.width / 2.0, self.height / 2.0);
-        let r1 = w2 * self.width_direction;
-        let r2 = h2 * height_direction;
-        let pi = intersect - self.position;
-        let pi_dot_r1 = pi.dot(&r1);
-        let pi_dot_r2 = pi.dot(&r2);
-
-        if pi_dot_r1 >= -w2 && pi_dot_r1 <= w2 && pi_dot_r2 >= -h2 && pi_dot_r2 <= h2 {
-            return Some(Intersection {
-                point: intersect,
-                normal: self.normal,
-                distance: t,
-            });
-        }
-        None
+        Some(Intersection {
+            point: intersect,
+            normal: Vector3::new(0.0, 0.0, 1.0),
+            distance: t,
+        })
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let bl = self.bl + Vector3::new(0.0, 0.0, -0.1);
+        let tr = self.tr + Vector3::new(0.0, 0.0, 0.1);
+        AABB::new(bl, tr)
     }
 }
 
@@ -451,19 +429,18 @@ impl Primitive for Rectangle {
 pub struct Cube {
     bln: Point3<f64>,
     trf: Point3<f64>,
-    bounding_box: AABB,
 }
 
 impl Cube {
-    pub fn new(bln: Point3<f64>, trf: Point3<f64>) -> Rc<dyn Primitive> {
-        Rc::new(Cube {
-            bln,
-            trf,
-            bounding_box: AABB { bln, trf },
-        })
+    pub fn new(bln: Point3<f64>, trf: Point3<f64>) -> Arc<dyn Primitive> {
+        Arc::new(Cube { bln, trf })
     }
 
-    pub fn unit() -> Rc<dyn Primitive> {
+    pub fn new_unboxed(bln: Point3<f64>, trf: Point3<f64>) -> Cube {
+        Cube { bln, trf }
+    }
+
+    pub fn unit() -> Arc<dyn Primitive> {
         let bln = Point3::new(-1.0, -1.0, -1.0);
         let trf = Point3::new(1.0, 1.0, 1.0);
         Cube::new(bln, trf)
@@ -488,35 +465,42 @@ impl Primitive for Cube {
             let intersect = ray.at_t(tmin);
 
             // Check if the intersection is outside the box
-            if intersect.x < bln.x
-                || intersect.x > trf.x
-                || intersect.y < bln.y
-                || intersect.y > trf.y
-                || intersect.z < bln.z
-                || intersect.z > trf.z
+            if intersect.x < bln.x - EPSILON
+                || intersect.x > trf.x + EPSILON
+                || intersect.y < bln.y - EPSILON
+                || intersect.y > trf.y + EPSILON
+                || intersect.z < bln.z - EPSILON
+                || intersect.z > trf.z + EPSILON
             {
                 return None; // Intersection is outside the box
             }
 
-            //Get normal of intersection point
-            //t1 is bln t2 is trf
-            let normal = if tmin == t1.x {
-                Vector3::new(-1.0, 0.0, 0.0)
-            } else if tmin == t1.y {
-                Vector3::new(0.0, -1.0, 0.0)
-            } else if tmin == t1.z {
-                Vector3::new(0.0, 0.0, -1.0)
-            } else if tmin == t2.x {
-                Vector3::new(1.0, 0.0, 0.0)
-            } else if tmin == t2.y {
-                Vector3::new(0.0, 1.0, 0.0)
-            } else {
-                Vector3::new(0.0, 0.0, 1.0)
-            };
+            // Determine which face was hit by finding the t-value closest to tmin
+            let diffs = [
+                (t1.x - tmin).abs(),
+                (t1.y - tmin).abs(),
+                (t1.z - tmin).abs(),
+                (t2.x - tmin).abs(),
+                (t2.y - tmin).abs(),
+                (t2.z - tmin).abs(),
+            ];
+            let normals = [
+                Vector3::new(-1.0, 0.0, 0.0),
+                Vector3::new(0.0, -1.0, 0.0),
+                Vector3::new(0.0, 0.0, -1.0),
+                Vector3::new(1.0, 0.0, 0.0),
+                Vector3::new(0.0, 1.0, 0.0),
+                Vector3::new(0.0, 0.0, 1.0),
+            ];
+            let min_idx = diffs.iter()
+                .enumerate()
+                .min_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
+                .unwrap().0;
+            let normal = normals[min_idx];
 
             Some(Intersection {
                 point: intersect,
-                normal: normal,
+                normal,
                 distance: tmin,
             })
         } else {
@@ -524,8 +508,8 @@ impl Primitive for Cube {
         }
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        AABB::new(self.bln, self.trf)
     }
 }
 
@@ -539,27 +523,17 @@ pub struct Triangle {
     v: Point3<f64>,
     w: Point3<f64>,
     normal: Vector3<f64>,
-    bounding_box: AABB,
 }
 
 impl Triangle {
-    pub fn new(u: Point3<f64>, v: Point3<f64>, w: Point3<f64>) -> Rc<dyn Primitive> {
+    pub fn new(u: Point3<f64>, v: Point3<f64>, w: Point3<f64>) -> Arc<dyn Primitive> {
         let uv = v - u;
         let uw = w - u;
         let normal = uw.cross(&uv).normalize();
-        let bln = u.inf(&v).inf(&w);
-        let trf = u.sup(&v).sup(&w);
-        let bounding_box = AABB { bln, trf };
-        Rc::new(Triangle {
-            u,
-            v,
-            w,
-            normal,
-            bounding_box,
-        })
+        Arc::new(Triangle { u, v, w, normal })
     }
     #[allow(dead_code)]
-    pub fn unit() -> Rc<dyn Primitive> {
+    pub fn unit() -> Arc<dyn Primitive> {
         let u = Point3::new(-1.0, -1.0, 0.0);
         let v = Point3::new(0.0, 1.0, 0.0);
         let w = Point3::new(1.0, -1.0, 0.0);
@@ -608,8 +582,13 @@ impl Primitive for Triangle {
         None
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let u = self.u;
+        let v = self.v;
+        let w = self.w;
+        let bln = u.inf(&v).inf(&w);
+        let trf = u.sup(&v).sup(&w);
+        AABB::new(bln, trf)
     }
 }
 
@@ -617,17 +596,13 @@ impl Primitive for Triangle {
 #[derive(Clone)]
 pub struct Mesh {
     triangles: Vec<Triangle>,
-    bounding_box: AABB,
 }
 
 impl Mesh {
-    pub fn new(triangles: Vec<Triangle>) -> Rc<dyn Primitive> {
+    pub fn new(triangles: Vec<Triangle>) -> Arc<dyn Primitive> {
         // Calculate the bounding box for the entire mesh based on the bounding boxes of individual triangles
-        let bounding_box = Mesh::compute_bounding_box(&triangles);
-        Rc::new(Mesh {
-            triangles,
-            bounding_box,
-        })
+        let _bounding_box = Mesh::compute_bounding_box(&triangles);
+        Arc::new(Mesh { triangles })
     }
 
     fn compute_bounding_box(triangles: &Vec<Triangle>) -> AABB {
@@ -641,10 +616,10 @@ impl Mesh {
             trf = trf.sup(&triangle.v);
             trf = trf.sup(&triangle.w);
         }
-        AABB { bln, trf }
+        AABB::new(bln, trf)
     }
 
-    pub fn from_file(filename: &str) -> Rc<dyn Primitive> {
+    pub fn from_file(filename: &str) -> Arc<dyn Primitive> {
         let mut triangles: Vec<Triangle> = Vec::new();
         let mut vertices: Vec<Point3<f64>> = Vec::new();
 
@@ -682,19 +657,10 @@ impl Mesh {
                                 let u = vertices[v1 - 1];
                                 let v = vertices[v2 - 1];
                                 let w = vertices[v3 - 1];
-                                let uv = u - v;
-                                let uw = w - v;
-                                let normal = uv.cross(&uw).normalize();
-                                let bln = u.inf(&v).inf(&w);
-                                let trf = u.sup(&v).sup(&w);
-                                let bounding_box = AABB { bln, trf };
-                                triangles.push(Triangle {
-                                    u,
-                                    v,
-                                    w,
-                                    normal,
-                                    bounding_box,
-                                });
+                                let uv = v - u;
+                                let uw = w - u;
+                                let normal = uw.cross(&uv).normalize();
+                                triangles.push(Triangle { u, v, w, normal });
                             }
                         }
                         _ => {}
@@ -727,8 +693,8 @@ impl Primitive for Mesh {
         closest_intersect
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        Mesh::compute_bounding_box(&self.triangles)
     }
 }
 
@@ -737,18 +703,14 @@ impl Primitive for Mesh {
 pub struct Torus {
     inner_rad: f64,
     outer_rad: f64,
-    bounding_box: AABB,
 }
 
 impl Torus {
-    pub fn new(inner_rad: f64, outer_rad: f64) -> Rc<dyn Primitive> {
+    pub fn new(inner_rad: f64, outer_rad: f64) -> Arc<dyn Primitive> {
         // I need to find the bounding box for this shape
-        let trf = Point3::new(1.0, 1.0, 1.0);
-        let bln = Point3::new(-1.0, -1.0, -1.0);
-        Rc::new(Torus {
+        Arc::new(Torus {
             inner_rad,
             outer_rad,
-            bounding_box: AABB { bln, trf },
         })
     }
 }
@@ -854,93 +816,87 @@ impl Primitive for Torus {
         })
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let extent = self.inner_rad + self.outer_rad;
+        let bln = Point3::new(-extent, -extent, -self.outer_rad);
+        let trf = Point3::new(extent, extent, self.outer_rad);
+        AABB::new(bln, trf)
     }
 }
 
 // GNOMON -----------------------------------------------------------------
 #[derive(Clone)]
 pub struct Gnonom {
-    x_cube: Rc<dyn Primitive>,
-    y_cube: Rc<dyn Primitive>,
-    z_cube: Rc<dyn Primitive>,
-    bounding_box: AABB,
+    x_cube: Cube,
+    y_cube: Cube,
+    z_cube: Cube,
 }
 
 impl Gnonom {
     const GNONOM_WIDTH: f64 = 0.1;
     const GNONOM_LENGTH: f64 = 2.0;
-    pub fn new() -> Rc<dyn Primitive> {
-        let x_cube = Cube::new(
+    pub fn new() -> Arc<dyn Primitive> {
+        let x_cube = Cube::new_unboxed(
             Point3::new(0.0, -Self::GNONOM_WIDTH, -Self::GNONOM_WIDTH),
             Point3::new(Self::GNONOM_LENGTH, Self::GNONOM_WIDTH, Self::GNONOM_WIDTH),
         );
-        let y_cube = Cube::new(
+        let y_cube = Cube::new_unboxed(
             Point3::new(-Self::GNONOM_WIDTH, 0.0, -Self::GNONOM_WIDTH),
             Point3::new(Self::GNONOM_WIDTH, Self::GNONOM_LENGTH, Self::GNONOM_WIDTH),
         );
-        let z_cube = Cube::new(
+        let z_cube = Cube::new_unboxed(
             Point3::new(-Self::GNONOM_WIDTH, -Self::GNONOM_WIDTH, 0.0),
             Point3::new(Self::GNONOM_WIDTH, Self::GNONOM_WIDTH, Self::GNONOM_LENGTH),
         );
-        let bounding_box = AABB {
-            bln: Point3::new(
-                -Self::GNONOM_WIDTH,
-                -Self::GNONOM_WIDTH,
-                -Self::GNONOM_WIDTH,
-            ),
-            trf: Point3::new(
-                Self::GNONOM_LENGTH,
-                Self::GNONOM_LENGTH,
-                Self::GNONOM_LENGTH,
-            ),
-        };
-        Rc::new(Gnonom {
+        Arc::new(Gnonom {
             x_cube,
             y_cube,
             z_cube,
-            bounding_box,
         })
     }
 }
 
 impl Primitive for Gnonom {
     fn intersect_ray(&self, ray: &Ray) -> Option<Intersection> {
-        match self.x_cube.intersect_ray(ray) {
-            Some(intersect) => return Some(intersect),
-            None => (),
-        };
-        match self.y_cube.intersect_ray(ray) {
-            Some(intersect) => return Some(intersect),
-            None => (),
-        };
-        match self.z_cube.intersect_ray(ray) {
-            Some(intersect) => return Some(intersect),
-            None => (),
-        };
-        None
+        let mut closest: Option<Intersection> = None;
+        let mut closest_dist = f64::MAX;
+
+        for cube in [&self.x_cube, &self.y_cube, &self.z_cube] {
+            if let Some(intersect) = cube.intersect_ray(ray) {
+                let dist = distance(&ray.a, &intersect.point);
+                if dist < closest_dist {
+                    closest_dist = dist;
+                    closest = Some(intersect);
+                }
+            }
+        }
+        closest
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        AABB::new(
+            Point3::new(
+                -Self::GNONOM_WIDTH,
+                -Self::GNONOM_WIDTH,
+                -Self::GNONOM_WIDTH,
+            ),
+            Point3::new(
+                Self::GNONOM_LENGTH,
+                Self::GNONOM_LENGTH,
+                Self::GNONOM_LENGTH,
+            ),
+        )
     }
 }
 
 // CROSS CAP ---------
 #[derive(Clone)]
-pub struct CrossCap {
-    bounding_box: AABB,
-}
+pub struct CrossCap {}
 
 impl CrossCap {
-    pub fn new() -> Rc<dyn Primitive> {
+    pub fn new() -> Arc<dyn Primitive> {
         // I need to find the bounding box for this shape
-        let trf = Point3::new(1.0, 1.0, 1.0);
-        let bln = Point3::new(-1.0, -1.0, -1.0);
-        Rc::new(CrossCap {
-            bounding_box: AABB { bln, trf },
-        })
+        Arc::new(CrossCap {})
     }
 }
 
@@ -1014,8 +970,10 @@ impl Primitive for CrossCap {
         })
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let trf = Point3::new(1.0, 1.0, 1.0);
+        let bln = Point3::new(-1.0, -1.0, -1.0);
+        AABB::new(bln, trf)
     }
 }
 
@@ -1024,19 +982,12 @@ impl Primitive for CrossCap {
 pub struct CrossCap2 {
     p: f64,
     q: f64,
-    bounding_box: AABB,
 }
 
 impl CrossCap2 {
-    pub fn new(p: f64, q: f64) -> Rc<dyn Primitive> {
+    pub fn new(p: f64, q: f64) -> Arc<dyn Primitive> {
         // I need to find the bounding box for this shape
-        let trf = Point3::new(1.0, 1.0, 1.0);
-        let bln = Point3::new(-1.0, -1.0, -1.0);
-        Rc::new(CrossCap2 {
-            p,
-            q,
-            bounding_box: AABB { bln, trf },
-        })
+        Arc::new(CrossCap2 { p, q })
     }
 }
 
@@ -1135,25 +1086,21 @@ impl Primitive for CrossCap2 {
         })
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let trf = Point3::new(1.0, 1.0, 1.0);
+        let bln = Point3::new(-1.0, -1.0, -1.0);
+        AABB::new(bln, trf)
     }
 }
 
 //  Steiner  ---------
 #[derive(Clone)]
-pub struct Steiner {
-    bounding_box: AABB,
-}
+pub struct Steiner {}
 
 impl Steiner {
-    pub fn new() -> Rc<dyn Primitive> {
+    pub fn new() -> Arc<dyn Primitive> {
         // I need to find the bounding box for this shape
-        let trf = Point3::new(1.0, 1.0, 1.0);
-        let bln = Point3::new(-1.0, -1.0, -1.0);
-        Rc::new(Steiner {
-            bounding_box: AABB { bln, trf },
-        })
+        Arc::new(Steiner {})
     }
 }
 
@@ -1216,25 +1163,21 @@ impl Primitive for Steiner {
         })
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let trf = Point3::new(1.0, 1.0, 1.0);
+        let bln = Point3::new(-1.0, -1.0, -1.0);
+        AABB::new(bln, trf)
     }
 }
 
 //  Steiner 2 ---------
 #[derive(Clone)]
-pub struct Steiner2 {
-    bounding_box: AABB,
-}
+pub struct Steiner2 {}
 
 impl Steiner2 {
-    pub fn new() -> Rc<dyn Primitive> {
+    pub fn new() -> Arc<dyn Primitive> {
         // I need to find the bounding box for this shape
-        let trf = Point3::new(1.0, 1.0, 1.0);
-        let bln = Point3::new(-1.0, -1.0, -1.0);
-        Rc::new(Steiner2 {
-            bounding_box: AABB { bln, trf },
-        })
+        Arc::new(Steiner2 {})
     }
 }
 
@@ -1308,8 +1251,10 @@ impl Primitive for Steiner2 {
         })
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let trf = Point3::new(1.0, 1.0, 1.0);
+        let bln = Point3::new(-1.0, -1.0, -1.0);
+        AABB::new(bln, trf)
     }
 }
 
@@ -1317,18 +1262,12 @@ impl Primitive for Steiner2 {
 #[derive(Clone)]
 pub struct Roman {
     k: f64,
-    bounding_box: AABB,
 }
 
 impl Roman {
-    pub fn new(k: f64) -> Rc<dyn Primitive> {
+    pub fn new(k: f64) -> Arc<dyn Primitive> {
         // I need to find the bounding box for this shape
-        let trf = Point3::new(1.0, 1.0, 1.0);
-        let bln = Point3::new(-1.0, -1.0, -1.0);
-        Rc::new(Roman {
-            k,
-            bounding_box: AABB { bln, trf },
-        })
+        Arc::new(Roman { k })
     }
 }
 
@@ -1419,8 +1358,10 @@ impl Primitive for Roman {
         })
     }
 
-    fn intersect_bounding_box(&self, ray: &Ray) -> bool {
-        self.bounding_box.intersect_bounding_box(ray)
+    fn get_aabb(&self) -> AABB {
+        let trf = Point3::new(1.0, 1.0, 1.0);
+        let bln = Point3::new(-1.0, -1.0, -1.0);
+        AABB::new(bln, trf)
     }
 }
 

src/ray.rs

@@ -1,13 +1,13 @@
-use crate::{node::Node, scene::Scene, EPSILON};
-use nalgebra::{distance, Matrix4, Point3, Vector3};
+use crate::{bvh::BVH, light::Light, node::Node, scene::Scene, state::RaytracingOption, EPSILON};
+use nalgebra::{distance, Matrix3, Matrix4, Point3, Vector3};
 use rand;
 
-const MAX_DEPTH: u8 = 5;
-const DIFFUSE_RAYS: i8 = 5;
-const DIFFUSE_COEFFICIENT: f32 = 0.5;
-
 fn random_vec() -> Vector3<f64> {
-    Vector3::new(rand::random(), rand::random(), rand::random())
+    Vector3::new(
+        rand::random::<f64>() * 2.0 - 1.0,
+        rand::random::<f64>() * 2.0 - 1.0,
+        rand::random::<f64>() * 2.0 - 1.0,
+    )
 }
 fn random_unit_vec() -> Vector3<f64> {
     random_vec().normalize()
@@ -22,15 +22,17 @@ pub struct Intersection {
 }
 //Intersection point including point and normal
 impl Intersection {
-    pub fn transform(&self, trans: &Matrix4<f64>, inv_trans: &Matrix4<f64>) -> Intersection {
-        let point = trans.transform_point(&self.point);
-        let normal = inv_trans.transpose().transform_vector(&self.normal);
+    pub fn transform(&mut self, trans: &Matrix4<f64>, inv_trans: &Matrix4<f64>) -> Intersection {
         Intersection {
-            point,
-            normal,
+            point: trans.transform_point(&self.point),
+            normal: inv_trans.transpose().transform_vector(&self.normal),
             distance: self.distance,
         }
     }
+    pub fn transform_mut(&mut self, trans: &Matrix4<f64>, inv_transpose: &Matrix3<f64>) {
+        self.point = trans.transform_point(&self.point);
+        self.normal = inv_transpose * self.normal;
+    }
 }
 
 // Ray struct represents a ray in 3D space with a starting point 'a' and a direction 'b'
@@ -66,45 +68,45 @@ impl Ray {
             b: trans.transform_vector(&self.b),
         }
     }
+    //Transform mutably
+    pub fn transform_mut(&mut self, trans: &Matrix4<f64>) {
+        self.a = trans.transform_point(&self.a);
+        self.b = trans.transform_vector(&self.b);
+    }
     //This function will determine if the ray hits an object in the scene
-    pub fn hit_scene(&self, scene: &Scene) -> bool {
+    //This is not optimised as it does not include bounding boxes
+    pub fn hit_scene(ray: &Ray, scene: &Scene) -> bool {
         for (_, node) in &scene.nodes {
             if !node.active {
                 continue;
             }
             // Transform ray into local model cordinates
-            let ray = self.transform(&node.inv_model);
-            // Check bounding box intersection
-            if node.primitive.intersect_bounding_box(&ray) {
-                // Check primitive intersection
-                if node.primitive.intersect_ray(&ray).is_some() {
-                    return true;
-                }
+            if node.intersect_ray(&ray).is_some() {
+                return true;
             }
         }
         false
     }
     //This function find the closest intersection point of a ray with an object in the scene
-    pub fn closest_intersect<'a>(&'a self, scene: &'a Scene) -> Option<(&Node, Intersection)> {
+    //Also not optimised, as it does not include bounding boxes
+    pub fn closest_intersect<'a>(
+        ray: &'a Ray,
+        scene: &'a Scene,
+    ) -> Option<(&'a Node, Intersection)> {
         let mut closest_distance = f64::MAX;
         let mut closest_intersect: Option<(&Node, Intersection)> = None;
+        let ray_a = ray.a;
         for (_, node) in &scene.nodes {
+            //position of ray in world coords
             if !node.active {
                 continue;
             }
-            // Transform ray into local model cordinates
-            let ray = self.transform(&node.inv_model);
-            // Check bounding box intersection
-            if node.primitive.intersect_bounding_box(&ray) {
-                // Check primitive intersection
-                if let Some(intersect) = node.primitive.intersect_ray(&ray) {
-                    // Dont intersect with itself
-                    if intersect.distance < EPSILON {
-                        continue;
-                    }
+
+            if node.aabb.intersect_ray(&ray) {
+                //Check node intersection
+                if let Some(intersect) = node.intersect_ray(&ray) {
                     // Check for closest distance by converting to world coords
-                    let intersect = intersect.transform(&node.model, &node.inv_model);
-                    let distance = distance(&ray.a, &intersect.point);
+                    let distance = distance(&ray_a, &intersect.point);
                     if distance < closest_distance {
                         closest_distance = distance;
                         closest_intersect = Some((node, intersect));
@@ -115,17 +117,40 @@ impl Ray {
         closest_intersect
     }
     // This function takes a scene and returns the color of the point where the ray intersects the scene
-    pub fn shade_ray(&self, scene: &Scene, depth: u8) -> Option<Vector3<f32>> {
-        if depth == MAX_DEPTH {
+    pub fn shade_ray(
+        &self,
+        scene: &Scene,
+        depth: u8,
+        options: &RaytracingOption,
+        sbvh: &Option<BVH>,
+    ) -> Option<Vector3<f32>> {
+        //If we have exceeded depth then return
+        if depth == options.ray_depth {
             return None;
         }
-        match self.closest_intersect(scene) {
-            Some((node, intersect)) => {
-                Some(Ray::phong_shade_point(
-                    &scene, &self, &node, &intersect, depth,
-                )) // If there is an intersection, shade it
+        match sbvh {
+            //We have a bvh so use bvh traversal
+            Some(bvh) => {
+                //Intersect the scene with the bvh
+                if let Some((node, intersect)) = bvh.traverse(self, 0) {
+                    return Some(Ray::phong_shade_point(
+                        &scene, &self, &node, &intersect, depth, options, sbvh,
+                    ));
+                }
+                return None;
+            }
+            //We dont have a bvh so use generic algorithm
+            None => {
+                //No BVH given so intersect normally
+                match Ray::closest_intersect(self, scene) {
+                    Some((node, intersect)) => {
+                        Some(Ray::phong_shade_point(
+                            &scene, &self, &node, &intersect, depth, options, sbvh,
+                        )) // If there is an intersection, shade it
+                    }
+                    None => None, // If there is no intersection, return None
+                }
             }
-            None => None, // If there is no intersection, return None
         }
     }
 
@@ -136,16 +161,38 @@ impl Ray {
         node: &Node,
         intersect: &Intersection,
         depth: u8,
+        options: &RaytracingOption,
+        bvh: &Option<BVH>,
     ) -> Vector3<f32> {
         let normal = &intersect.normal;
-        let point = intersect.point;
+        let point = &intersect.point;
         let incidence = &ray.b;
-
         let material = &node.material;
 
-        // Compute the ambient light component and set it as base colour
         let mut colour = Vector3::zeros();
 
+        // Reflection is view-dependent, not light-dependent — compute once
+        let mut reflect = Vector3::zeros();
+        if options.reflect {
+            let reflect_dir = incidence - 2.0 * incidence.dot(&normal) * normal;
+            let reflect_ray = Ray::new(*point, reflect_dir);
+            if let Some(col) = reflect_ray.shade_ray(scene, depth + 1, options, bvh) {
+                reflect += col.component_mul(&material.kr)
+            }
+        }
+
+        // Indirect diffuse (global illumination samples) — compute once
+        let mut indirect = Vector3::zeros();
+        if options.diffuse {
+            for _ in 0..options.diffuse_rays {
+                let diffuse_dir = random_unit_vec();
+                let diffuse_ray = Ray::new(point.clone(), diffuse_dir + normal);
+                if let Some(col) = diffuse_ray.shade_ray(scene, depth + 1, options, bvh) {
+                    indirect += col * options.diffuse_coefficient;
+                }
+            }
+        }
+
         for (_, light) in &scene.lights {
             if !light.active {
                 continue;
@@ -161,67 +208,78 @@ impl Ray {
             let to_light = to_light.normalize();
 
             //Niave Shadows
-            let to_light_ray = Ray::new(point, to_light);
-            if to_light_ray.light_blocked(scene, node) {
-                continue;
+            if options.shadows {
+                let to_light_ray = Ray::new(*point, to_light);
+                if to_light_ray.light_blocked(scene, light, bvh) {
+                    continue;
+                }
             }
 
             let n_dot_l = normal.dot(&to_light).max(0.0) as f32;
 
-            //Reflected component
-            let mut reflect = Vector3::zeros();
-            let reflect_dir = incidence - 2.0 * incidence.dot(&normal) * normal;
-            let reflect_ray = Ray::new(point, reflect_dir);
-            if let Some(col) = reflect_ray.shade_ray(scene, depth + 1) {
-                reflect += col.component_mul(&material.kr)
-            }
-
-            //Diffuse component (Lambertian)
+            //Direct diffuse component (Lambertian)
             let mut diffuse = Vector3::zeros();
-            diffuse += material.kd * n_dot_l;
-            for _ in 0..DIFFUSE_RAYS {
-                let diffuse_dir = random_unit_vec();
-                let diffuse_ray = Ray::new(point.clone(), diffuse_dir + normal);
-                if let Some(col) = diffuse_ray.shade_ray(scene, depth + 1) {
-                    diffuse += col * DIFFUSE_COEFFICIENT;
-                }
+            if options.diffuse {
+                diffuse += material.kd * n_dot_l;
             }
 
             //Specular component
             let mut specular = Vector3::zeros();
-            if n_dot_l < 0.0 {
-                let h = (to_light - incidence).normalize();
-                let n_dot_h = normal.dot(&h).max(0.0) as f32;
-                specular = material.ks * n_dot_h.powf(material.shininess);
+            if options.specular {
+                if n_dot_l > 0.0 {
+                    let h = (to_light - incidence).normalize();
+                    let n_dot_h = normal.dot(&h).max(0.0) as f32;
+                    specular = material.ks * n_dot_h.powf(material.shininess);
+                }
             }
 
             //Falloff
-            // let falloff = 1.0
-            //     / (1.0
-            //         + light.falloff[0]
-            //         + light.falloff[1] * light_distance
-            //         + light.falloff[2] * light_distance * light_distance);
+            let mut falloff = 1.0;
+            if options.falloff {
+                falloff = 1.0
+                    / ((1.0 + light.falloff[0])
+                        + light.falloff[1] * light_distance
+                        + light.falloff[2] * light_distance * light_distance);
+            }
 
-            let intensity = light.colour.component_mul(&(diffuse + reflect + specular));
+            let intensity = light.colour.component_mul(&(diffuse + specular)) * falloff;
             colour += &intensity;
         }
 
+        // Add light-independent terms
+        colour += reflect + indirect;
+
         colour
     }
 
-    pub fn light_blocked(&self, scene: &Scene, _node: &Node) -> bool {
-        for (_, node) in &scene.nodes {
-            if !node.active {
-                continue;
+    pub fn light_blocked(&self, scene: &Scene, light: &Light, bvh: &Option<BVH>) -> bool {
+        let light_distance = distance(&self.a, &light.position);
+        match bvh {
+            Some(bvh) => {
+                //We have a bvh so use bvh traversal
+                if let Some((_, intersect)) = bvh.traverse(self, 0) {
+                    return intersect.distance < light_distance;
+                }
             }
-            let ray = self.transform(&node.inv_model);
-            if node.primitive.intersect_bounding_box(&ray) {
-                if node.primitive.intersect_ray(&ray).is_some() {
-                    return true;
+            None => {
+                for (_, node) in &scene.nodes {
+                    if !node.active {
+                        continue;
+                    }
+                    if node.aabb.intersect_ray(self) {
+                        match node.intersect_ray(self) {
+                            Some(intersect) => {
+                                if intersect.distance < light_distance {
+                                    return true;
+                                }
+                            }
+                            None => continue,
+                        }
+                    }
                 }
             }
         }
-        false
+        return false;
     }
     //Cast a set of rays
     pub fn cast_rays(

src/scene.rs

@@ -35,4 +35,10 @@ impl Scene {
     pub fn add_camera(&mut self, label: String, camera: Camera) {
         self.cameras.insert(label, camera);
     }
+    // Compute all matricies for nodes
+    pub fn compute(&mut self) {
+        for (_, node) in &mut self.nodes {
+            node.compute();
+        }
+    }
 }

src/state.rs

@@ -1,36 +1,87 @@
 //Use linear algebra module
 
+use crate::bvh::BVH;
 use crate::camera::Camera;
 use crate::ray::Ray;
 use crate::{gui::Gui, scene::Scene};
 use crate::{gui::GuiEvent, log_error};
+use std::collections::HashSet;
 use std::path::Path;
+use std::thread;
 
 use nalgebra::Vector3;
 use rand::seq::SliceRandom;
 use rand::{random, thread_rng};
 
 use std::error::Error;
+use std::sync::atomic::{AtomicBool, Ordering};
+use std::sync::{mpsc, Arc, Mutex};
 
 use anyhow::Result;
 use pixels::{Pixels, SurfaceTexture};
 use winit::dpi::{LogicalSize, PhysicalSize};
-use winit::event::{Event, KeyboardInput, MouseButton, VirtualKeyCode, WindowEvent};
+use winit::event::{
+    ElementState, Event, KeyboardInput, MouseButton, VirtualKeyCode, WindowEvent,
+};
 use winit::event_loop::{ControlFlow, EventLoop};
 use winit::window::{Window, WindowBuilder};
 
 const START_WIDTH: i32 = 1200;
 const START_HEIGHT: i32 = 700;
-const RAY_SAMPLES: i8 = 5;
-const RAY_RANDOMNESS: f64 = 100.0;
-const COLOUR_CLEAR: [u8; 4] = [0x22, 0x00, 0x11, 0x55];
-const PIXEL_CLEAR: [u8; 4] = [0x55, 0x00, 0x22, 0x55];
 
 pub const INIT_FILE: &str = "rhai/scene.rhai";
 pub const SAVE_FILE: &str = "img.png";
 
+#[derive(Clone)]
+pub struct RaytracingOption {
+    pub threads: u32,
+    pub ray_samples: u32,
+    pub ray_randomness: f64,
+    pub clear_color: [u8; 4],
+    pub pixel_clear: [u8; 4],
+    pub pixels_per_thread: u32,
+    pub buffer_proportion: f32,
+    pub buffer_fov: f64,
+    pub ray_depth: u8,
+    pub diffuse_rays: u8,
+    pub diffuse_coefficient: f32,
+    pub bvh_active: bool,
+    pub shadows: bool,
+    pub diffuse: bool,
+    pub reflect: bool,
+    pub specular: bool,
+    pub falloff: bool,
+}
+impl RaytracingOption {
+    pub fn default() -> RaytracingOption {
+        RaytracingOption {
+            threads: 12,
+            ray_samples: 1,
+            ray_randomness: 700.0,
+            clear_color: [0x22, 0x00, 0x11, 0x55],
+            pixel_clear: [0x11, 0x00, 0x22, 0x55],
+            pixels_per_thread: 100,
+            buffer_proportion: 1.0,
+            buffer_fov: 70.0,
+            ray_depth: 1,
+            diffuse_rays: 3,
+            diffuse_coefficient: 0.1,
+            bvh_active: false,
+            shadows: true,
+            diffuse: true,
+            reflect: true,
+            specular: true,
+            falloff: true,
+        }
+    }
+}
+
+const CAMERA_MOVE_SPEED: f64 = 0.15;
+const CAMERA_ORBIT_SPEED: f64 = 0.005;
+
 pub struct State {
-    scene: Scene,
+    scene: Arc<Scene>,
+    bvh: Arc<Option<BVH>>,
     camera: Camera,
     window: Window,
 
@@ -40,52 +91,78 @@ pub struct State {
     pixels: Pixels,
     gui: Gui,
 
-    rays: Vec<Ray>,
-    ray_queue: Vec<usize>,
+    rays: Arc<Vec<Ray>>,
+    ray_queue: Arc<Mutex<Vec<usize>>>,
+    raytracing_options: Arc<RaytracingOption>,
+
+    result_rx: mpsc::Receiver<Vec<(usize, [u8; 4])>>,
+    render_active: Arc<AtomicBool>,
+    rendering: bool,
+
+    keys_pressed: HashSet<VirtualKeyCode>,
+    right_mouse_down: bool,
+    last_mouse_pos: Option<(f64, f64)>,
+    camera_dirty: bool,
 }
 
 impl State {
     pub fn new(window: Window, pixels: Pixels, gui: Gui) -> Self {
-        let scene = Scene::empty();
+        let scene = Arc::new(Scene::empty());
         let window_size = window.inner_size();
+        let pixels = pixels;
         let camera = Camera::unit();
-        let rays = Vec::new();
+        let rays = Arc::new(Vec::new());
+        let (_tx, rx) = mpsc::channel();
 
         Self {
             scene,
+            bvh: Arc::new(None),
             camera,
             window,
             buffer_width: window_size.width as u32,
             buffer_height: window_size.height as u32,
-            pixels: pixels,
+            pixels,
             gui,
             rays,
-            ray_queue: Vec::new(),
+            ray_queue: Arc::new(Mutex::new(Vec::new())),
+            raytracing_options: Arc::new(RaytracingOption::default()),
+            result_rx: rx,
+            render_active: Arc::new(AtomicBool::new(false)),
+            rendering: false,
+            keys_pressed: HashSet::new(),
+            right_mouse_down: false,
+            last_mouse_pos: None,
+            camera_dirty: false,
         }
     }
 
     fn update(&mut self) -> Result<(), Box<dyn Error>> {
         if let Some(event) = self.gui.event.take() {
             match event {
-                GuiEvent::BufferResize(proportion, fov) => {
-                    self.resize_buffer(proportion, fov as f64)?
+                GuiEvent::RaytracerOption(options) => {
+                    self.raytracing_options = Arc::new(options);
+                    match self.raytracing_options.bvh_active {
+                        true => self.bvh = Arc::new(Some(BVH::build(&self.scene.nodes))),
+                        false => self.bvh = Arc::new(None),
+                    }
+                    self.resize_buffer()?
                 }
-                GuiEvent::CameraUpdate(camera, fovy) => {
-                    self.rays = Ray::cast_rays(
+                GuiEvent::CameraUpdate(camera) => {
+                    self.rays = Arc::new(Ray::cast_rays(
                         &camera.eye,
                         &camera.target,
                         &camera.up,
-                        fovy as f64,
+                        self.raytracing_options.buffer_fov,
                         self.buffer_width,
                         self.buffer_height,
-                    );
+                    ));
                     self.camera = camera;
-                    self.clear()?;
+                    self.clear_buffer()?;
                     self.reset_queue();
                 }
                 GuiEvent::SceneLoad(scene) => {
-                    self.scene = scene;
-                    self.clear()?;
+                    self.scene = Arc::new(scene);
+                    self.clear_buffer()?;
                     self.reset_queue();
                 }
                 GuiEvent::SaveImage(filename) => {
@@ -103,30 +180,32 @@ impl State {
         Ok(())
     }
 
-    fn resize_buffer(&mut self, proportion: f32, fovy: f64) -> Result<(), Box<dyn Error>> {
+    fn resize_buffer(&mut self) -> Result<(), Box<dyn Error>> {
         // Calculate new buffer dimensions based on proportion
         let size = self.window.inner_size();
+        let proportion = &self.raytracing_options.buffer_proportion;
+        let fovy = self.raytracing_options.buffer_fov;
         self.buffer_width = (size.width as f32 * proportion) as u32;
         self.buffer_height = (size.height as f32 * proportion) as u32;
 
         // Clear the buffer and reset the ray queue
-        self.clear()?;
+        self.clear_buffer()?;
         self.reset_queue();
 
         // Recalculate rays with new buffer dimensions
-        self.rays = Ray::cast_rays(
+        self.rays = Arc::new(Ray::cast_rays(
             &self.camera.eye,
             &self.camera.target,
             &self.camera.up,
             fovy,
             self.buffer_width,
             self.buffer_height,
-        );
+        ));
 
         // Resize buffer and surface
-        let pixels = &mut self.pixels;
-        pixels.resize_surface(size.width, size.height)?;
-        pixels.resize_buffer(self.buffer_width, self.buffer_height)?;
+        self.pixels.resize_surface(size.width, size.height)?;
+        self.pixels
+            .resize_buffer(self.buffer_width, self.buffer_height)?;
 
         Ok(())
     }
@@ -137,77 +216,227 @@ impl State {
     }
 
     fn keyboard_input(&mut self, key: &KeyboardInput) {
-        if let Some(VirtualKeyCode::A) = key.virtual_keycode {
-            // Handle 'A' key event here
-        }
-    }
-
-    fn mouse_input(&mut self, _button: &MouseButton) {
-        // Handle mouse input here
-    }
-
-    fn draw(&mut self) -> Result<(), Box<dyn Error>> {
-        //Draw ray_num in a block
-        let frame = self.pixels.frame_mut();
-        for _ in 0..self.gui.ray_num {
-            //Get random index from queue
-            let index = match self.ray_queue.pop() {
-                Some(index) => index,
-                None => break,
-            };
-            //Shade colour for selected ray
-            let mut colour = Vector3::zeros();
-            for _ in 0..RAY_SAMPLES {
-                let ray = &self.rays[index];
-                let point = ray.a;
-                let dir = ray.b;
-                let rx = (random::<f64>() - 0.5) / RAY_RANDOMNESS;
-                let ry = (random::<f64>() - 0.5) / RAY_RANDOMNESS;
-                let rz = (random::<f64>() - 0.5) / RAY_RANDOMNESS;
-                let nx = dir.x + rx;
-                let ny = dir.y + ry;
-                let nz = dir.z + rz;
-
-                let rand_ray = Ray::new(point, Vector3::new(nx, ny, nz));
-
-                if let Some(ray_colour) = rand_ray.shade_ray(&self.scene, 0) {
-                    colour += ray_colour;
-                };
+        if let Some(keycode) = key.virtual_keycode {
+            match key.state {
+                ElementState::Pressed => {
+                    self.keys_pressed.insert(keycode);
+                }
+                ElementState::Released => {
+                    self.keys_pressed.remove(&keycode);
+                }
             }
-            colour = (colour / RAY_SAMPLES as f32) * 255.0;
-            let rgba = [colour.x as u8, colour.y as u8, colour.z as u8, 0xff];
-            frame[index * 4..(index + 1) * 4].copy_from_slice(&rgba);
         }
-        Ok(())
     }
 
-    fn clear(&mut self) -> Result<(), Box<dyn Error>> {
+    fn mouse_input(&mut self, button: &MouseButton, state: &ElementState) {
+        if *button == MouseButton::Right {
+            self.right_mouse_down = *state == ElementState::Pressed;
+            if !self.right_mouse_down {
+                self.last_mouse_pos = None;
+            }
+        }
+    }
+
+    fn cursor_moved(&mut self, x: f64, y: f64) {
+        if self.right_mouse_down {
+            if let Some((last_x, last_y)) = self.last_mouse_pos {
+                let dx = x - last_x;
+                let dy = y - last_y;
+                self.camera.orbit(
+                    -dx * CAMERA_ORBIT_SPEED,
+                    -dy * CAMERA_ORBIT_SPEED,
+                );
+                self.camera_dirty = true;
+            }
+            self.last_mouse_pos = Some((x, y));
+        }
+    }
+
+    fn process_camera_movement(&mut self) {
+        let speed = CAMERA_MOVE_SPEED;
+
+        if self.keys_pressed.contains(&VirtualKeyCode::W) {
+            self.camera.move_forward(speed);
+            self.camera_dirty = true;
+        }
+        if self.keys_pressed.contains(&VirtualKeyCode::S) {
+            self.camera.move_forward(-speed);
+            self.camera_dirty = true;
+        }
+        if self.keys_pressed.contains(&VirtualKeyCode::A) {
+            self.camera.move_right(-speed);
+            self.camera_dirty = true;
+        }
+        if self.keys_pressed.contains(&VirtualKeyCode::D) {
+            self.camera.move_right(speed);
+            self.camera_dirty = true;
+        }
+        if self.keys_pressed.contains(&VirtualKeyCode::Q) {
+            self.camera.move_up(-speed);
+            self.camera_dirty = true;
+        }
+        if self.keys_pressed.contains(&VirtualKeyCode::E) {
+            self.camera.move_up(speed);
+            self.camera_dirty = true;
+        }
+
+        if self.camera_dirty {
+            self.camera_dirty = false;
+            self.rays = Arc::new(Ray::cast_rays(
+                &self.camera.eye,
+                &self.camera.target,
+                &self.camera.up,
+                self.raytracing_options.buffer_fov,
+                self.buffer_width,
+                self.buffer_height,
+            ));
+            self.gui.update_camera(&self.camera);
+            let _ = self.clear_buffer();
+            self.reset_queue();
+        }
+    }
+
+    fn draw(&mut self) {
+        if !self.rendering {
+            return;
+        }
+
+        // Drain completed results from background workers
+        loop {
+            match self.result_rx.try_recv() {
+                Ok(results) => {
+                    let frame = self.pixels.frame_mut();
+                    for (index, rgba) in results {
+                        frame[index * 4..(index + 1) * 4].copy_from_slice(&rgba);
+                    }
+                }
+                Err(mpsc::TryRecvError::Empty) => break,
+                Err(mpsc::TryRecvError::Disconnected) => {
+                    // All worker threads have finished
+                    self.rendering = false;
+                    self.gui.stop_render_timer();
+                    break;
+                }
+            }
+        }
+    }
+
+    fn clear_buffer(&mut self) -> Result<(), Box<dyn Error>> {
         let frame = self.pixels.frame_mut();
         for pixel in frame.chunks_exact_mut(4) {
-            pixel.copy_from_slice(&COLOUR_CLEAR);
+            pixel.copy_from_slice(&self.raytracing_options.pixel_clear);
         }
         Ok(())
     }
 
     fn reset_queue(&mut self) {
+        // Signal any existing workers to stop
+        self.render_active.store(false, Ordering::Relaxed);
+
+        match self.raytracing_options.bvh_active {
+            true => self.bvh = Arc::new(Some(BVH::build(&self.scene.nodes))),
+            false => self.bvh = Arc::new(None),
+        }
+
+        // Create new shuffled queue
         let size = self.buffer_height as usize * self.buffer_width as usize;
         let mut ray_queue: Vec<usize> = (0..size).collect();
         ray_queue.shuffle(&mut thread_rng());
-        self.ray_queue = ray_queue;
+        self.ray_queue = Arc::new(Mutex::new(ray_queue));
+
+        // Create new channel and active flag
+        let (tx, rx) = mpsc::channel();
+        self.result_rx = rx;
+        let render_active = Arc::new(AtomicBool::new(true));
+        self.render_active = render_active.clone();
+        self.rendering = true;
+
+        // Spawn persistent worker threads
+        let num_threads = self.raytracing_options.threads;
+        let pixels_per_thread = self.raytracing_options.pixels_per_thread;
+
+        for _ in 0..num_threads {
+            let rays = self.rays.clone();
+            let scene = self.scene.clone();
+            let options = self.raytracing_options.clone();
+            let bvh = self.bvh.clone();
+            let queue = self.ray_queue.clone();
+            let tx = tx.clone();
+            let active = render_active.clone();
+
+            thread::spawn(move || {
+                let randomness = options.ray_randomness;
+                let samples = options.ray_samples;
+                let samples_f32 = samples as f32;
+
+                loop {
+                    if !active.load(Ordering::Relaxed) {
+                        break;
+                    }
+
+                    // Pop a batch from the shared queue
+                    let load: Vec<usize> = {
+                        let mut q = queue.lock().unwrap();
+                        let mut batch = Vec::with_capacity(pixels_per_thread as usize);
+                        for _ in 0..pixels_per_thread {
+                            match q.pop() {
+                                Some(index) => batch.push(index),
+                                None => break,
+                            }
+                        }
+                        batch
+                    };
+
+                    if load.is_empty() {
+                        break;
+                    }
+
+                    // Process the batch
+                    let mut results = Vec::with_capacity(load.len());
+                    for index in &load {
+                        let mut colour: Vector3<f32> = Vector3::zeros();
+                        let ray = &rays[*index];
+                        for _ in 0..samples {
+                            let point = ray.a;
+                            let dir = ray.b;
+                            let rx = (random::<f64>() - 0.5) / randomness;
+                            let ry = (random::<f64>() - 0.5) / randomness;
+                            let rz = (random::<f64>() - 0.5) / randomness;
+                            let nx = dir.x + rx;
+                            let ny = dir.y + ry;
+                            let nz = dir.z + rz;
+
+                            let rand_ray = Ray::new(point, Vector3::new(nx, ny, nz));
+
+                            if let Some(ray_colour) =
+                                rand_ray.shade_ray(&scene, 0, &options, &bvh)
+                            {
+                                colour += ray_colour;
+                            }
+                        }
+                        colour = (colour / samples_f32) * 255.0;
+                        let rgba = [colour.x as u8, colour.y as u8, colour.z as u8, 0xff];
+                        results.push((*index, rgba));
+                    }
+
+                    // Send results back to main thread
+                    if tx.send(results).is_err() {
+                        break;
+                    }
+                }
+            });
+        }
+        // Drop our copy of tx so the channel disconnects when all workers finish
+        drop(tx);
+
+        self.gui.start_render_timer();
     }
 
     fn render(&mut self) -> Result<(), Box<dyn Error>> {
         // Update state
         self.update()?;
-        // Draw rays if we have remaining rays in queue
-        if !self.ray_queue.is_empty() {
-            match self.draw() {
-                Err(e) => {
-                    println!("ERROR: {}", e);
-                }
-                _ => {}
-            }
-        }
+        // Collect completed rays from background workers
+        self.draw();
         // Render Gui
         self.gui
             .prepare(&self.window)
@@ -234,7 +463,7 @@ pub fn run() -> Result<(), Box<dyn Error>> {
     let gui = Gui::new(&window, &pixels);
 
     let mut state = State::new(window, pixels, gui);
-    state.resize_buffer(1.0, 90.0)?;
+    state.resize_buffer()?;
 
     event_loop.run(move |event, _, control_flow| {
         state.gui.handle_event(&state.window, &event);
@@ -248,11 +477,17 @@ pub fn run() -> Result<(), Box<dyn Error>> {
                 WindowEvent::CloseRequested => *control_flow = ControlFlow::Exit,
                 WindowEvent::Resized(size) => state.resize(&size).expect("Window Resize Error"),
                 WindowEvent::KeyboardInput { input, .. } => state.keyboard_input(&input),
-                WindowEvent::MouseInput { button, .. } => state.mouse_input(&button),
+                WindowEvent::MouseInput { button, state: elem_state, .. } => {
+                    state.mouse_input(&button, &elem_state)
+                }
+                WindowEvent::CursorMoved { position, .. } => {
+                    state.cursor_moved(position.x, position.y)
+                }
                 _ => {}
             },
 
             Event::RedrawRequested(_) => {
+                state.process_camera_movement();
                 if let Err(_e) = state.render() {
                     *control_flow = ControlFlow::Exit;
                 }