Basic raytracer written in Rust

Question

I come from a fairly strong background of C and thought that this project would be a good way to get a handle on Rust. Right now, I have everything in one file because I wasn't sure the best way to organize the code (at least not yet). Let me know if I should post anything else (the .toml for the project, things like that).

extern crate image;
extern crate nalgebra;

use image::{Rgb, RgbImage};
use nalgebra::base::{Unit, Vector3};
use nalgebra::geometry::Point3;

trait Hittable {
    fn intersects(&self, ray: &Ray, tmin: f64, tmax: f64, record: &mut HitRecord) -> bool;
}

struct HitRecord {
    t: f64, // time of hit along ray
    n: Unit<Vector3<f64>>, // normal of surface at point
    p: Point3<f64>, // point of intersection
}

impl HitRecord {
    fn new(t: f64, n: Unit<Vector3<f64>>, p: Point3<f64>) -> Self { Self {t, n, p } }
}

struct Ray {
    origin: Point3<f64>,
    dir: Vector3<f64>,
}

impl Ray {
    fn new(origin: Point3<f64>, dir: Vector3<f64>) -> Self { Self { origin, dir } }

    fn at(&self, t: f64) -> Point3<f64> {
        self.origin + self.dir.scale(t)
    }
}

struct Sphere {
    center: Point3<f64>,
    r: f64,
}

impl Sphere {
    fn new(center: Point3<f64>, r: f64) -> Self { Self { center, r } }
}

impl Hittable for Sphere {
    
    fn intersects(&self, ray: &Ray, tmin: f64, tmax: f64, hit_record: &mut HitRecord) -> bool {
        let diff: Vector3<f64> = ray.origin - self.center;
        // get quadratic equation, calculate discriminant
        let a = ray.dir.dot(&ray.dir);
        let b = diff.dot(&ray.dir);
        let c = diff.dot(&diff) - self.r * self.r;
        let disc = b * b - a * c;
        if disc < 0.0 {
            return false; // no need to fill data
        }
        let root = disc.sqrt();
        let ans = (-b - root) / a; // try first solution to equation
        if ans < tmax && ans > tmin {
            hit_record.t = ans;
            hit_record.p = ray.at(ans);
            hit_record.n = Unit::new_normalize(self.center - hit_record.p);
            return true;
        } else {
            // is putting this in an else block necessary? I tried without the else
            // and the compiler said 'if may be missing an else clause', and I'm
            // still not completely sure why that is.
            let ans = (-b + root) / a;
            if ans < tmax && ans > tmin {
                hit_record.t = ans;
                hit_record.p = ray.at(ans);
                hit_record.n = Unit::new_normalize(self.center - hit_record.p);
                return true;
            } else {
                return false;
            }
        }
    }
}


fn main() {
    let image_width: u32 = 512;
    let aspect_ratio = 3.0 / 2.0;
    let image_height: u32 = ((image_width as f64) / aspect_ratio).round() as u32;
    let viewport_height = 2.0;
    let viewport_width = viewport_height * aspect_ratio;
    let focal_length = 1.0;

    let origin: Point3<f64> = Point3::origin();
    let horizontal_offset: Vector3<f64> = Vector3::new(viewport_width, 0.0, 0.0);
    let vertical_offset: Vector3<f64> = Vector3::new(0.0, viewport_height, 0.0);
    // this is the point in world space that represents the bottom left corner of the plane that is being projected onto
    let bottom_left_corner: Point3<f64> = origin - horizontal_offset.scale(0.5) - vertical_offset.scale(0.5) - Vector3::new(0.0, 0.0, focal_length);
    let mut img = RgbImage::new(image_width, image_height);

    let sphere: Sphere = Sphere::new(Point3::new(0.0, 0.0, -1.0), 0.5);

    let sphere_array = [sphere];
    let light: Point3<f64> = Point3::new(0.0, 0.0, 0.0);

    for i in 0u32..image_width {
        for j in 0u32..image_height {
            let u: f64 = (i as f64) / ((image_width - 1) as f64);
            let v: f64 = (j as f64) / ((image_height - 1) as f64);
            let to: Point3<f64> = bottom_left_corner + horizontal_offset.scale(u) + vertical_offset.scale(v);
            let dir: Vector3<f64> = to - origin;
            let ray = Ray::new(origin, dir);
            let color: Rgb<u8> = cast_ray(&ray, &sphere_array, &light);
            img.put_pixel(i, j, color);
        }
    }

    img.save("test.png").unwrap();

}

fn cast_ray(ray: &Ray, array: &[Sphere], light: &Point3<f64>) -> Rgb<u8> {
    // start time at -1 to see if it changes later
    let mut hit_record: HitRecord = HitRecord::new(-1.0, Unit::new_unchecked(Vector3::new(1.0, 0.0, 0.0)), Point3::new(0.0, 0.0, 0.0));
    for sphere in array.iter() {
        if sphere.intersects(ray, 0.0, 10.0, &mut hit_record) {
            break; // this won't work for multiple spheres (yet), need to find closest
        }
    }
    if hit_record.t < 0.0 { // miss
        return Rgb([55, 155, 255]);
    } else {
        let hit: Point3<f64> = hit_record.p;
        let normal: Unit<Vector3<f64>> = hit_record.n;
        let light_dir: Unit<Vector3<f64>> = Unit::new_normalize(hit - light);
        let brightness: f64 = light_dir.dot(normal.as_ref()).max(0.0);
        return Rgb([((255 as f64) * brightness) as u8, 0, 0]);
    }
    
}

A lot of this code feels kind of clunky. For example, in the cast_ray function at the bottom, is there a better way to initialize the hit_record variable that will later be overwritten? Normally in C, I would make those null, but I obviously can't do that and don't want to make them Optional.

In general, I just would like to know if this is 'good' Rust code and follows the generally accepted practices. For example, I assume having an output variable isn't considered best practice, but what would be the alternative? Returning a tuple? Because then I worry that it wouldn't be as fast if I made a new hit_record inside each intersects method and return that.

Finally, in case it's helpful, the libraries I'm using are image to write to image files, and nalgebra to handle a lot of the math for me.

Answer 1

I'm not particularly knowledgeable on Rust, so this will be a light review.

I would suggest leaning on the Rust tooling quite a bit, especially when new to a language. I've heard nothing but praise for the format tool, compiler warnings, and linters.

Running the code through rustfmt only moves newlines and comments around, so good job on having a clean starting point.

---

Clippy (the linter) points out that there are unecessary explicit returns. You can remove the return keyword (and trailing semicolon) from lines at the end of functions like return RGB(...);. For example

if ans < tmax && ans > tmin {
    hit_record.t = ans;
    hit_record.p = ray.at(ans);
    hit_record.n = Unit::new_normalize(self.center - hit_record.p);
    true
} else {
    // is putting this in an else block necessary? I tried without the else
    // and the compiler said 'if may be missing an else clause', and I'm
    // still not completely sure why that is.
    let ans = (-b + root) / a;
    if ans < tmax && ans > tmin {
        hit_record.t = ans;
        hit_record.p = ray.at(ans);
        hit_record.n = Unit::new_normalize(self.center - hit_record.p);
        true
    } else {
        false
    }
}

If you don't like the implicit returns, you can write code with guard clauses instead of if/else (BTW I couldn't reproduce the compiler error you get removing the else block).

if ans < tmax && ans > tmin {
    hit_record.t = ans;
    hit_record.p = ray.at(ans);
    hit_record.n = Unit::new_normalize(self.center - hit_record.p);
    return true;
}

if ans < tmax && ans > tmin {
    hit_record.t = ans;
    hit_record.p = ray.at(ans);
    hit_record.n = Unit::new_normalize(self.center - hit_record.p);
    return true;
}
false

The final warning clippy gives is on the cast of a literal, it suggests changing 255 as f64 to 255_f64, including the type directly in the literal.

---

// get quadratic equation, calculate discriminant
let a = ray.dir.dot(&ray.dir);
let b = diff.dot(&ray.dir);
let c = diff.dot(&diff) - self.r * self.r;
let disc = b * b - a * c;

Do you have a link/diagram/short explanation for this? I don't know the derivation, so I don't know why the discriminant is not the usual disc = b * b - 4 * a * c. Where did the 4 go?

---

let ans = (-b - root) / a;
if ans < tmax && ans > tmin {
    // Block that involves ans
}

let ans = (-b + root) / a;
if ans < tmax && ans > tmin {
    // Block that involves ans
}

I can't tell if this suggestion is actually an improvement without profiling, but if most of the time these conditions fail, you could trade two divisions for two multiplications by scaling tmin and tmax, and only dividing if it passes the checks.

let a_tmin = a * tmin;
let a_tmax = a * tmax;
let scaled_ans = -b - root;
if a_tmin < scaled_ans && scaled_ans < a_tmax {
    let ans = scaled_ans / a;
    // Block that involves ans
}

let scaled_ans = -b + root;
if a_tmin < scaled_ans && scaled_ans < a_tmax {
    let ans = scaled_ans / a;
    // Block that involves ans
}

---

I think changing the intersects function to return an optional hit_record would solve two problems. It can save you from initialising a dummy hit_record, and also removes the output variable. I'm not sure if this code is exactly correct, but something along the lines of

trait Hittable {
    fn intersects(&self, ray: &Ray, tmin: f64, tmax: f64) -> Option<HitRecord>;
}

...

fn intersects(&self, ray: &Ray, tmin: f64, tmax: f64) -> Option<HitRecord> {
    let diff: Vector3<f64> = ray.origin - self.center;
    // get quadratic equation, calculate discriminant
    let a = ray.dir.dot(&ray.dir);
    let b = diff.dot(&ray.dir);
    let c = diff.dot(&diff) - self.r * self.r;
    let disc = b * b - a * c;
    if disc < 0.0 {
        return None; // no need to fill data
    }
    let root = disc.sqrt();
    let ans = (-b - root) / a; // try first solution to equation
    if ans < tmax && ans > tmin {
        let p = ray.at(ans);
        return Some(HitRecord {
            t: ans,
            p,
            n: Unit::new_normalize(self.center - p),
        });
    }

    let ans = (-b + root) / a;
    if ans < tmax && ans > tmin {
        let p = ray.at(ans);
        return Some(HitRecord {
            t: ans,
            p,
            n: Unit::new_normalize(self.center - p),
        });
    }
    None
}


...


    for sphere in array.iter() {
    let result = sphere.intersects(ray, 0.0, 10.0);
    match result {
        Some(record) => {
            hit_record = record;
            break;
        }
        None => continue,
    }

This does highlight that there is repeated code for both potential roots, maybe a quick function to deduplicate would be appropriate here.

This then means the pattern in cast_rays is now "loop over an iterator until the first value is found", so we can turn to the standard library if we would like to. Usually deferring to the standard library is a good idea, the implementation will be well tested, and will require less maintenance on our behalf. I think find_map is the algorithm we want in this case. This is where my Rust knowledge fails me, I think the code would be something like

let hit_record = array
    .iter()
    .find_map(|sphere| sphere.intersects(ray, 0.0, 10.0))
    .unwrap();  // What should happen in there is no intersecting spheres?

But at this point I'm out of my depth.