Basic raytracing

src/main.rs

@@ -1,11 +1,11 @@
 use cgmath::prelude::*;
 
-use cgmath::{vec3, vec4, InnerSpace, Matrix4, Point3, Vector3};
+use cgmath::{vec3, vec4, InnerSpace, Matrix4, Point3, Vector3, Zero};
 
 use image::ImageBuffer;
 use utilities::prelude::*;
 
-use std::f32::{consts::PI as M_PI, MAX};
+use std::f32::{consts::PI as M_PI, MAX, MIN};
 
 struct Triangle {
     points: [Vector3<f32>; 3],
@@ -25,6 +25,92 @@ struct View {
     fov: f32,
 }
 
+struct Ray {
+    origin: Vector3<f32>,
+    direction: Vector3<f32>,
+    inv_direction: Vector3<f32>,
+    signs: Vector3<usize>,
+}
+
+impl Ray {
+    pub fn new(origin: Vector3<f32>, direction: Vector3<f32>) -> Ray {
+        let inv_direction = 1.0 / direction;
+
+        let x = (inv_direction.x < 0.0) as usize;
+        let y = (inv_direction.y < 0.0) as usize;
+        let z = (inv_direction.z < 0.0) as usize;
+
+        Ray {
+            origin,
+            direction,
+            inv_direction,
+            signs: Vector3::new(x, y, z),
+        }
+    }
+}
+
+struct AABB {
+    bounds: [Vector3<f32>; 2],
+    start_index: usize,
+    end_index: usize,
+}
+
+impl AABB {
+    fn new() -> AABB {
+        AABB {
+            bounds: [Vector3::zero(), Vector3::zero()],
+            start_index: 0,
+            end_index: 0,
+        }
+    }
+
+    fn intersect(&self, ray: &Ray) -> bool {
+        let mut tmin = (self.bounds[ray.signs.x].x - ray.origin.x) * ray.inv_direction.x;
+        let mut tmax = (self.bounds[1 - ray.signs.x].x - ray.origin.x) * ray.inv_direction.x;
+        let tymin = (self.bounds[ray.signs.y].y - ray.origin.y) * ray.inv_direction.y;
+        let tymax = (self.bounds[1 - ray.signs.y].y - ray.origin.y) * ray.inv_direction.y;
+
+        if tmin > tymax || tymin > tmax {
+            return false;
+        }
+
+        if tymin > tmin {
+            tmin = tymin;
+        }
+
+        if tymax < tmax {
+            tmax = tymax;
+        }
+
+        let tzmin = (self.bounds[ray.signs.z].z - ray.origin.z) * ray.inv_direction.z;
+        let tzmax = (self.bounds[1 - ray.signs.z].z - ray.origin.z) * ray.inv_direction.z;
+
+        if tmin > tzmax || tzmin > tmax {
+            return false;
+        }
+
+        if tzmin > tmin {
+            tmin = tzmin;
+        }
+
+        if tzmax < tmax {
+            tmax = tzmax;
+        }
+
+        let mut t = tmin;
+
+        if t < 0.0 {
+            t = tmax;
+
+            if t < 0.0 {
+                return false;
+            }
+        }
+
+        true
+    }
+}
+
 fn main() {
     let input_data = [Triangle::new(
         vec3(0.0, 0.0, 0.0),
@@ -38,53 +124,44 @@ fn main() {
         fov: 45.0,
     };
 
-    let (origin, direction) = calculate_ray(640, 360, 1280, 720, &view);
+    debug_raytracer(1280, 720, &view, &input_data);
-
-    let color = pixel_color(origin, direction, &input_data);
-
-    //debug_raytracer(1280, 720, &view, &input_data);
 }
 
 fn f_to_u(color: f32) -> u8 {
     ((color * 255.0) / 1.0) as u8
 }
 
-fn calculate_ray(
+fn calculate_ray(x: u32, y: u32, dim_x: u32, dim_y: u32, view: &View) -> Ray {
-    x: u32,
-    y: u32,
-    dim_x: u32,
-    dim_y: u32,
-    view: &View,
-) -> (Vector3<f32>, Vector3<f32>) {
     let aspect_ratio = dim_x as f32 / dim_y as f32;
 
-    let p_x = (2.0 * ((x as f32 + 0.5) / dim_x as f32) - 1.0)
+    let fov = view.fov / 180.0 * M_PI * 2.0;
-        * (view.fov / 2.0 * M_PI / 180.0).tan()
+
-        * aspect_ratio;
+    let p_x = (1.0 - 2.0 * ((x as f32 + 0.5) / dim_x as f32)) * fov * aspect_ratio;
-    let p_y = (1.0 - 2.0 * (((dim_y - y) as f32 + 0.5) / dim_y as f32))
+    let p_y = (2.0 * (((y) as f32 + 0.5) / dim_y as f32) - 1.0) * fov;
-        * (view.fov / 2.0 * M_PI / 180.0).tan();
 
     let camera_to_world = Matrix4::look_at(
         Point3::from_vec(view.position),
         Point3::from_vec(view.look_at),
         vec3(0.0, 0.0, 1.0),
     )
-    .inverse_transform()
+    .invert()
     .unwrap();
 
     let origin = camera_to_world * vec4(0.0, 0.0, 0.0, 1.0);
     let direction = camera_to_world * vec4(p_x, p_y, -1.0, 1.0);
 
-    (origin.truncate(), direction.truncate().normalize())
+    Ray::new(origin.truncate(), -direction.truncate().normalize())
 }
 
 fn debug_raytracer(dim_x: u32, dim_y: u32, view: &View, data: &[Triangle]) {
     let mut imgbuf = ImageBuffer::new(dim_x, dim_y);
 
     for (x, y, pixel) in imgbuf.enumerate_pixels_mut() {
-        let (origin, direction) = calculate_ray(x, y, dim_x, dim_y, view);
+        let ray = calculate_ray(x, y, dim_x, dim_y, view);
 
-        let color = pixel_color(origin, direction, data);
+        let acceleration_data = create_acceleration_data(data);
+
+        let color = pixel_color(&ray, data, acceleration_data);
 
         *pixel = image::Rgb([f_to_u(color.x), f_to_u(color.y), f_to_u(color.z)]);
     }
@@ -92,21 +169,74 @@ fn debug_raytracer(dim_x: u32, dim_y: u32, view: &View, data: &[Triangle]) {
     imgbuf.save("raytrace.png").unwrap();
 }
 
-fn pixel_color(orig: Vector3<f32>, dir: Vector3<f32>, data: &[Triangle]) -> Vector3<f32> {
+fn create_acceleration_data(input_data: &[Triangle]) -> Vec<AABB> {
+    let mut acceleration_data = Vec::new();
+
+    let accel_count = (input_data.len() as f32 / 24.0).ceil() as usize;
+
+    for i in 0..accel_count {
+        let mut accel = AABB::new();
+
+        let mut bbmin = vec3(MAX, MAX, MAX);
+        let mut bbmax = vec3(MIN, MIN, MIN);
+
+        // 3 vertices per triangles
+        // 8 triangles per acceleration structure
+        let start_index = 8 * i;
+        let mut end_index = 8 * (i + 1);
+
+        // check if end_index exceeds input data length
+        if end_index > input_data.len() {
+            end_index = input_data.len();
+        }
+
+        // set index information
+        accel.start_index = start_index;
+        accel.end_index = end_index;
+
+        // create bounding box
+        for k in start_index..end_index {
+            for vertex in input_data[k].points.iter() {
+                // check minimum values
+                bbmin.x = bbmin.x.min(vertex.x);
+                bbmin.y = bbmin.y.min(vertex.y);
+                bbmin.z = bbmin.z.min(vertex.z);
+
+                // check maximum values
+                bbmax.x = bbmax.x.max(vertex.x);
+                bbmax.y = bbmax.y.max(vertex.y);
+                bbmax.z = bbmax.z.max(vertex.z);
+            }
+        }
+
+        accel.bounds[0] = bbmin;
+        accel.bounds[1] = bbmax;
+
+        acceleration_data.push(accel);
+    }
+
+    acceleration_data
+}
+
+fn pixel_color(ray: &Ray, data: &[Triangle], acceleration_data: Vec<AABB>) -> Vector3<f32> {
     let mut final_color = vec3(0.0, 1.0, 0.0);
     let mut closest_value = MAX;
     let mut closest_index = -1;
 
-    for i in 0..data.len() {
+    for accel in acceleration_data.iter() {
-        let v0 = data[i].points[0];
+        if accel.intersect(ray) {
-        let v1 = data[i].points[1];
+            for i in accel.start_index..accel.end_index {
-        let v2 = data[i].points[2];
+                let v0 = data[i].points[0];
-        let mut t = 0.0;
+                let v1 = data[i].points[1];
+                let v2 = data[i].points[2];
+                let mut t = 0.0;
 
-        if ray_triangle_intersect_naive(orig, dir, v0, v1, v2, &mut t) {
+                if ray_triangle_intersect_mt(ray, v0, v1, v2, &mut t) {
-            if t < closest_value {
+                    if t < closest_value {
-                closest_index = i as i32;
+                        closest_index = i as i32;
-                closest_value = t;
+                        closest_value = t;
+                    }
+                }
             }
         }
     }
@@ -120,8 +250,7 @@ fn pixel_color(orig: Vector3<f32>, dir: Vector3<f32>, data: &[Triangle]) -> Vect
 
 // source: https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/moller-trumbore-ray-triangle-intersection
 fn ray_triangle_intersect_mt(
-    orig: Vector3<f32>,
+    ray: &Ray,
-    dir: Vector3<f32>,
     v0: Vector3<f32>,
     v1: Vector3<f32>,
     v2: Vector3<f32>,
@@ -129,7 +258,7 @@ fn ray_triangle_intersect_mt(
 ) -> bool {
     let v0v1 = v1 - v0;
     let v0v2 = v2 - v0;
-    let pvec = dir.cross(v0v2);
+    let pvec = ray.direction.cross(v0v2);
     let det = v0v1.dot(pvec);
 
     // ray and triangle are parallel if det is close to 0
@@ -139,7 +268,7 @@ fn ray_triangle_intersect_mt(
 
     let inv_det = 1.0 / det;
 
-    let tvec = orig - v0;
+    let tvec = ray.origin - v0;
     let u = tvec.dot(pvec) * inv_det;
 
     if u < 0.0 || u > 1.0 {
@@ -147,7 +276,7 @@ fn ray_triangle_intersect_mt(
     }
 
     let qvec = tvec.cross(v0v1);
-    let v = dir.dot(qvec) * inv_det;
+    let v = ray.direction.dot(qvec) * inv_det;
 
     if v < 0.0 || u + v > 1.0 {
         return false;
@@ -157,75 +286,3 @@ fn ray_triangle_intersect_mt(
 
     true
 }
-
-fn ray_triangle_intersect_naive(
-    orig: Vector3<f32>,
-    dir: Vector3<f32>,
-    v0: Vector3<f32>,
-    v1: Vector3<f32>,
-    v2: Vector3<f32>,
-    t: &mut f32,
-) -> bool {
-    // compute plane's normal
-    let v0v1 = v1 - v0;
-    let v0v2 = v2 - v0;
-    // no need to normalize
-    let N = v0v1.cross(v0v2);
-    let area2 = N.magnitude();
-
-    // Step 1: finding P
-
-    // check if ray and plane are parallel ?
-    let NdotRayDirection = N.dot(dir);
-    if NdotRayDirection.abs() < 0.001
-    // almost 0
-    {
-        return false;
-    } // they are parallel so they don't intersect !
-
-    // compute d parameter using equation 2
-    let d = N.dot(v0);
-
-    // compute t (equation 3)
-    *t = (N.dot(orig) + d) / NdotRayDirection;
-
-    // check if the triangle is in behind the ray
-    if *t < 0.0 {
-        return false;
-    } // the triangle is behind
-
-    // compute the intersection point using equation 1
-    let P = orig + *t * dir;
-
-    // Step 2: inside-outside test
-    let mut C; // vector perpendicular to triangle's plane
-
-    // edge 0
-    let edge0 = v1 - v0;
-    let vp0 = P - v0;
-
-    C = edge0.cross(vp0);
-    if N.dot(C) < 0.0 {
-        return false;
-    } // P is on the right side
-
-    // edge 1
-    let edge1 = v2 - v1;
-    let vp1 = P - v1;
-
-    C = edge1.cross(vp1);
-    if N.dot(C) < 0.0 {
-        return false;
-    } // P is on the right side
-
-    // edge 2
-    let edge2 = v0 - v2;
-    let vp2 = P - v2;
-
-    C = edge2.cross(vp2);
-    if N.dot(C) < 0.0 {
-        return false;
-    } // P is on the right side;
-
-    true // this ray hits the triangle
-}