init: create git repo

.gitignore

@@ -0,0 +1,3 @@
+.ccls-cache/
+bin/*
+obj/*

Makefile

@@ -0,0 +1,29 @@
+SHELL = /bin/bash
+
+SRC_DIR = ./src
+INC_DIR = ./include
+OBJ_DIR = ./obj
+
+HEADERS   = $(wildcard $(INC_DIR)/*.hpp)
+SRC_FILES = $(wildcard $(SRC_DIR)/*.cpp)
+OBJ_FILES = $(addprefix $(OBJ_DIR)/,$(notdir $(SRC_FILES:.cpp=.o)))
+
+CXXFLAGS += -Wall -std=c++2a -fopenmp -O3 -Ofast -march=native
+
+all: bin/raytracer
+
+$(OBJ_DIR)/%.o: src/%.cpp
+	$(CXX) $(CXXFLAGS) -c -o $@ $^
+
+bin/raytracer: $(OBJ_FILES)
+	@mkdir -p $(OBJ_DIR)
+	$(CXX) $(CXXFLAGS) -o $@ $^
+
+.PHONY: debug
+debug: CXXFLAGS += -O0 -g
+debug: clean bin/raytracer
+
+.PHONY: clean
+clean:
+	$(RM) -r obj/*.o
+	$(RM) bin/*

README

@@ -0,0 +1,27 @@
+A raytracer that will render to any ANSI terminal with support for 24-bit
+colours, at least 80x24 big.
+
+Build with `make' and run the resulting binary in `/bin'.
+
+The raytracer implements the following:
+ - Planes, spheres, and point lights
+ - Reflection
+ - IOR-based refraction
+ - Cook-Torrance BRDF (a few other models are left in the code)
+ - Primitive object materials
+
+You can move around with WASD and SPC/Z; zoom with I/O, change the exposure with
++/-, and quit with Q.
+
+The image output is a 24-bit coloured 80x24 (by default, tweakable in
+`main.cpp'), 2x supersampled. The output uses Unicode characters and differently
+coloured foreground and background to effectively double the output resolution.
+You can read the implementation in `buffer.cpp'.
+
+The raytracing loop is parallelised using an OpenMP pragma, which should be well
+supported on most systems. You can just comment out the pragma (and remove the
+omp flag from the Makefile) to remove parallelisation.
+
+Although the code isn't particularly documented, most functions have a docstring
+in their respective header files. Most fun things to play around with are
+located in `main.cpp'.

include/buffer.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include "./vec.hpp"
+#include "./raytrace.hpp"
+
+struct Frame {
+	v3* data;
+};
+
+void bg_fg(v3* cols, v3& bg, v3& fg, const char*& c);
+
+inline int colcvt(choice_t x) {
+	if (x < 0.0) return 0;
+	if (x > 1.0) return 255;
+	return int(x * 255.0);
+}
+
+void display_frame(const Screen& screen, const Frame& frame);
+
+void render_frame(const Screen& screen, Camera& camera, Frame& frame,
+		  const std::vector<Renderable*>& scene,
+		  const std::vector<Renderable*>& objects,
+		  const std::vector<Renderable*>& lights);

include/objects.hpp

@@ -0,0 +1,38 @@
+#pragma once
+
+#include "./vec.hpp"
+#include "./raytrace.hpp"
+
+#include <stdio.h>
+
+struct Plane : Renderable {
+    v3 n;
+    choice_t offset;
+
+    int surface_type;
+
+    Plane(const v3& n, choice_t offset, const Material& specs, int surface_type);
+    Plane(const v3& n, choice_t offset, const Material& specs);
+
+    inline v3 normal(const v3& _) { return n; }
+    int intersect(const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal) override;
+};
+
+struct Sphere : Renderable {
+    v3 o;
+    choice_t r;
+
+    int surface_type;
+
+    Sphere(const v3& origin, choice_t r, const Material& specs, int surface_type);
+    Sphere(const v3& origin, choice_t r);
+
+    inline v3 normal(const v3& xyz) const { return (xyz - o).norm(); }
+    int intersect(const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal) override;
+    int transmit(const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal) override;
+};
+
+struct Light : Sphere {
+    Light(const v3& origin, choice_t r, const v3& colour);
+    Light(const v3& origin, choice_t r, const Material& specs);
+};

include/raytrace.hpp

@@ -0,0 +1,82 @@
+#pragma once
+
+#include <vector>
+#include "./vec.hpp"
+
+struct Screen {
+    const int W, H;
+    const choice_t ASPECT_RATIO;
+    constexpr Screen(int W, int H, choice_t ASPECT_RATIO)
+	: W(W), H(H), ASPECT_RATIO(ASPECT_RATIO) {}
+};
+
+struct Camera {
+    v3 pos; v3 rotation; v3 translation{0,0,0};
+    choice_t zoom = 0.5;
+    choice_t exposure = -0.6;
+    choice_t dynamic_range = 1;
+};
+
+struct Material {
+    choice_t specular = 0.1;
+    choice_t roughness = 0.5;
+    choice_t ior = 1.2;
+    choice_t transparency = 0.8;
+    v3 colour;
+};
+
+struct Renderable {
+    Material specs;
+
+    Renderable(const Material& specs);
+
+    virtual int intersect(const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal)=0;
+
+    /**
+       Like #intersect, but for transmission rays inside the object. #eye must
+       be on the surface of the object.
+     */
+    virtual int transmit(const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal) { return -1; }
+};
+
+/**
+   Shoots a ray from #eye in direction #dir, (possibly) hitting the closest object from #objs.
+   Expects #dir to be normalised.
+
+   The hit object's distance, location, normal and type will be written in the
+   references so-aptly named. If nothing was hit, the contents of those variables is UB.
+
+   Returns the pointer to the object it hit, or `nullptr' otherwise.
+
+   Respects the initial value of #hit_dist, et al.
+
+   Optional parameter #self is a pointer to an object considered 'self' that
+   will be skipped during collision checks. The objects are compared by pointer
+   address.
+ */
+Renderable* shoot_ray_into_objects(const v3& eye, const v3& dir, const std::vector<Renderable*>& objs,
+				   choice_t& hit_dist, v3& hit_loc, v3& hit_normal, int& hit_type,
+				   Renderable* self = nullptr);
+
+/**
+   Shoots a ray from #eye in (unit) direction #dir into a #scene and returns the
+   colour of that ray.
+
+   #scene must be #objects u #lights.
+
+   #depth is how much the ray tracer is allowed to recurse.
+ */
+v3 ray_trace(const v3& eye, const v3& dir,
+	     const std::vector<Renderable*>& scene,
+	     const std::vector<Renderable*>& objects,
+	     const std::vector<Renderable*>& lights, int depth,
+	     Renderable* skip = nullptr,
+	     choice_t* dist_to_hit = nullptr);
+
+choice_t ggx_ndf(const v3& dir, const v3& h, const v3& normal, choice_t roughness);
+choice_t blinn_phong_ndf(const v3& dir, const v3& h, const v3& normal, choice_t roughness);
+choice_t cook_torrance_gaf(const v3& dir, const v3& h, const v3& normal, const v3& light);
+choice_t cook_torrance_fresnel(const v3& dir, const v3& h, choice_t ior);
+
+v3 cook_torrance(const v3& dir, const v3& normal, const v3& light, choice_t specular, choice_t roughness,
+		 choice_t ior, const v3& light_colour, const v3& surface_colour);

include/vec.hpp

@@ -0,0 +1,105 @@
+#pragma once
+
+#include <math.h>
+
+typedef float choice_t;
+
+#define do_op(o, r, i)							\
+    inline void operator o##= (const r & rhs) { for(unsigned n=0; n<(i); n++) c[n] o##= rhs.c[n]; } \
+    inline void operator o##= (choice_t d) { for(unsigned n=0; n<(i); n++) c[n] o##= d; } \
+    inline r operator o (const r & rhs) const { r tmp(*this); tmp o##= rhs; return tmp; } \
+    inline r operator o (choice_t d) const { r tmp(*this); tmp o##= d; return tmp; }
+
+#define do_ops(c, i)		\
+    do_op(*, c, i);		\
+    do_op(/, c, i);		\
+    do_op(+, c, i);		\
+    do_op(-, c, i);
+
+
+struct v3 {
+    union {
+	choice_t c[4];
+	struct {
+	    choice_t x, y, z, w;
+	};
+    };
+
+    v3(choice_t x, choice_t y, choice_t z) : c{x,y,z} {}
+    v3(choice_t f) : v3(f,f,f) {}
+    v3() : v3(0,0,0) {}
+
+    do_ops(v3, 3);
+
+    inline v3 operator-(void) const { return {-x, -y, -z}; }
+
+    inline choice_t operator^(const v3& rhs) const
+    { return x * rhs.x + y * rhs.y + z * rhs.z; }
+
+    /* Return normalised (unit) vector. */
+    inline v3 norm() const
+    { return *this / len(); }
+
+    inline choice_t len() const
+    { return sqrtf(x*x + y*y + z*z); }
+
+    inline choice_t lensquared() const
+    { return x*x + y*y + z*z; }
+
+    /* Return vector reflected around `n'. Expects `this' and `n' to be unit. */
+    inline v3 ref(const v3& n) const
+    { return n * ((*this ^ n) * 2) - *this;}
+
+    /* Cross product. */
+    inline v3 cross(const v3& b) const
+    { return { y*b.z - z*b.y,
+	       x*b.z - z*b.x,
+	       x*b.y - y*b.x }; }
+
+    inline v3 refract(const v3& n, choice_t ior_ratio) const
+    {
+	return n * sqrtf(1 - powf(ior_ratio, 2) * (1 - powf(*this ^ n, 2))) + (*this - n * (*this ^ n)) * ior_ratio;
+    }
+
+    inline v3 refract(const v3& n, choice_t ior_incident, choice_t ior_refract) const
+    { return refract(n, ior_incident / ior_refract); }
+
+    inline v3 clamp(const v3& from, const v3& to) const
+    { return { fmin(fmax(x, from.x), to.x),
+	       fmin(fmax(y, from.y), to.y),
+	       fmin(fmax(z, from.z), to.z) }; }
+
+    inline v3 exp() const
+    { return { ::exp(x), ::exp(y), ::exp(z) }; }
+
+    inline v3 htan(choice_t exposure = 0) const
+    {
+	return { powf((tanh(x + exposure) + 1) / 2, 2.2), powf((tanh(y + exposure) + 1) / 2, 2.2), powf((tanh(z + exposure) + 1) / 2, 2.2) };
+	return { tanh(x + exposure), tanh(y + exposure), tanh(z + exposure) };
+    }
+};
+
+#undef do_ops
+#undef do_op
+
+#define V3_FMT "v3(%.2f, %.2f, %.2f)"
+#define V3_ARG(v) (v).x, (v).y, (v).z
+
+struct Matrix
+{
+    v3 m[4];
+    void init_rotate(const v3& angle)
+    {
+        choice_t Cx = cos(angle.c[0]), Cy = cos(angle.c[1]), Cz = cos(angle.c[2]);
+        choice_t Sx = sin(angle.c[0]), Sy = sin(angle.c[1]), Sz = sin(angle.c[2]);
+        choice_t sxsz = Sx*Sz, cxsz = Cx*Sz;
+        choice_t cxcz = Cx*Cz, sxcz = Sx*Cz;
+        Matrix result = {{ { Cy*Cz, Cy*Sz, -Sy },
+                           { sxcz*Sy - cxsz, sxsz*Sy + cxcz, Sx*Cy },
+                           { cxcz*Sy + sxsz, cxsz*Sy - sxcz, Cx*Cy } }};
+        *this = result;
+    }
+
+    v3 transform(const v3& vec) const
+    { return { m[0] ^ vec, m[1] ^ vec, m[2] ^ vec }; }
+};

src/buffer.cpp

@@ -0,0 +1,143 @@
+#include "../include/buffer.hpp"
+#include <string>
+
+void bg_fg(v3* cols, v3& bg, v3& fg, const char*& c)
+{
+    // 0 0 0 0
+    // | | | |
+    // | | | \_ top left
+    // | | \_ top right
+    // | \_ bottom left
+    // \_ bottom right
+
+    static const char8_t* chars[16] = {
+	u8"░", // 0000
+	u8"▝", // 0001
+	u8"▘", // 0010
+	u8"▀", // 0011
+	u8"▖", // 0100
+	u8"▌", // 0101
+	u8"▞", // 0110
+	u8"▛", // 0111
+	u8"▗", // 1000
+	u8"▚", // 1001
+	u8"▐", // 1010
+	u8"▜", // 1011
+	u8"▄", // 1100
+	u8"▙", // 1101
+	u8"▟", // 1110
+	u8"▓", // 1111
+    };
+
+    v3 ava = cols[0];
+    v3 avb = cols[1];
+    {
+	choice_t max = -1.0;
+	for (int i = 1; i < 4; ++i) {
+	    choice_t len = (cols[i] - cols[0]).lensquared();
+	    if (len > max) {
+		max = len;
+		avb = cols[i];
+	    }
+	}
+    }
+    int m = 0b0000;
+
+    for (int i = 0; i < 4; ++i) {
+	v3 nava {};
+	v3 navb {};
+	int na = 0;
+	int nb = 0;
+
+	for (int j = 0; j < 4; ++j) {
+	    float lena = (ava - cols[j]).len();
+	    float lenb = (avb - cols[j]).len();
+
+	    if (lena > lenb) {
+		m |= 1 << j;
+		nava += cols[j];
+		++na;
+	    } else {
+		m &= ~(1 << j);
+		navb += cols[j];
+		++nb;
+	    }
+	}
+
+	ava = nava / na;
+	avb = navb / nb;
+    }
+
+    bg = ava;
+    fg = avb;
+    c = (const char*) chars[m];
+    // c = (const char*) chars[15];
+}
+
+void display_frame(const Screen& screen, const Frame& frame)
+{
+    char chr[64];
+    for (int y = 0; y < screen.H; y += 2) {
+	std::string text;
+	for (int x = 0; x < screen.W; x += 2) {
+	    v3 cols[4] = {
+		frame.data[(y + 0) * screen.W + (x + 0)],
+		frame.data[(y + 0) * screen.W + (x + 1)],
+		frame.data[(y + 1) * screen.W + (x + 0)],
+		frame.data[(y + 1) * screen.W + (x + 1)]
+	    };
+
+	    const char* c;
+	    v3 bg;
+	    v3 fg;
+
+	    bg_fg(cols, bg, fg, c);
+
+	    snprintf(chr, 64, "\033[38;2;%d;%d;%dm\033[48;2;%d;%d;%dm%s",
+		     colcvt(bg.x), colcvt(bg.y), colcvt(bg.z),
+		     colcvt(fg.x), colcvt(fg.y), colcvt(fg.z),
+		     c
+		);
+	    text += chr;
+	}
+	printf("%.*s\033[0m\n", (int)text.size(), text.data());
+    }
+}
+
+
+void render_frame(const Screen& screen, Camera& camera, Frame& frame,
+		  const std::vector<Renderable*>& scene,
+		  const std::vector<Renderable*>& objects,
+		  const std::vector<Renderable*>& lights)
+{
+    Matrix m1; m1.init_rotate({camera.rotation.x,0,0});
+    Matrix m2; m2.init_rotate({0,camera.rotation.y,0});
+
+    // The translation vector tells us where to translate the camera,
+    // in 'camera tangent space'. So, we rotate the translation vector,
+    // add it to camera's position, and then reset it to 0.
+    camera.pos += m2.transform(m1.transform(camera.translation));
+    camera.translation = {};
+
+    #pragma omp parallel for
+    for(int y = 0; y < screen.H; y++) {
+	for (int x = 0; x < screen.W; x++) {
+	    v3 camray = {
+		(choice_t)x / (choice_t)(screen.W) - 0.5f,
+		(choice_t)y / (choice_t)(screen.H) - 0.5f,
+		camera.zoom
+	    };
+	    camray.x *= screen.ASPECT_RATIO; // aspect ratio correction
+	    camray = m2.transform(m1.transform(camray.norm()));
+
+	    v3 pix = ray_trace(camera.pos, camray, scene, objects, lights, 3);
+	    pix = pix.htan(camera.exposure);
+
+	    if (pix.x == 0.0) {
+		frame.data[y * screen.W + x] = v3 {0.0, 0.0, 0.0};
+	    } else {
+		frame.data[y * screen.W + x] = pix;
+	    }
+	}
+    }
+}

src/main.cpp

@@ -0,0 +1,176 @@
+#include <string>
+#include <chrono>
+#include <iostream>
+#include <vector>
+#include <stdio.h>
+#include "../include/objects.hpp"
+#include "../include/buffer.hpp"
+
+#include <termios.h>
+#include <unistd.h>
+#include <algorithm>
+
+static std::vector<Renderable*> lights {
+    new Light(v3{0, 0, 19}, 2, v3{1, 1, 1}),
+    new Light(v3{-7, 15, 0}, 2, v3{1, 1, 1}),
+    new Light(v3{17, -9, 16}, 2, v3{1, 1, 1}),
+
+    // new Light(v3{0, 0, 19}, 2, v3{0.8, 0.8, 0.9}),
+    // new Light(v3{-7, 15, 0}, 2, v3{0.5, 0.8, 0.9}),
+    // new Light(v3{17, -9, 16}, 2, v3{0.1, 0.8, 0.4}),
+};
+
+static const Material PLANE_MAT = {
+    .specular = 0.12,
+    .roughness = 0.55,
+    .ior = 2,
+    .transparency = 0,
+    .colour = {1, 1, 1},
+};
+
+static const Material SPHERE_MAT_TRANS = {
+    .specular = 0.85,
+    .roughness = 0.5,
+    .ior = 1.2,
+    .transparency = 0.95,
+    .colour = {0, 1, 0},
+};
+
+static const Material SPHERE_MAT_METAL = {
+    .specular = 0.87,
+    .roughness = 0.53,
+    .ior = 2.3,
+    .transparency = 0.0,
+    .colour = {0, 0, 0},
+};
+
+
+static std::vector<Renderable*> objects
+{
+    new Sphere(v3{10, 0, 6}, 3.0, SPHERE_MAT_TRANS, 1),
+    new Sphere(v3{2, -2, 6}, 3.0, SPHERE_MAT_METAL, 1),
+
+    new Plane{v3{1,  0,  0}, -30, PLANE_MAT},
+    new Plane{v3{-1, 0,  0}, -30, PLANE_MAT},
+    new Plane{v3{0,  1,  0}, -30, PLANE_MAT},
+    new Plane{v3{0, -1,  0}, -30, PLANE_MAT},
+    new Plane{v3{0,  0,  1}, -30, PLANE_MAT},
+    new Plane{v3{0,  0, -1}, -30, PLANE_MAT},
+};
+
+static std::vector<Renderable*> scene;
+
+// static constexpr int W = 294*2, H = 82*2;
+// static constexpr int W = 2 * 640 * 2.9 * 1.8, H = 2 * 480 * 1.5 * 2.5;
+static constexpr int W = 2*80, H = 2*24;
+static constexpr choice_t ASPECT_RATIO = W/(H * 1.9);
+
+static constexpr Screen screen = {W, H, ASPECT_RATIO};
+static Camera camera = { {15.31, -3.3, 7.72}, {0.56, 1.84, 0}, {} };
+
+static inline void gotoxy(int x, int y)
+{ printf("%c[%d;%dH",0x1B, (y+1), (x+1)); }
+
+char getch()
+{
+    char c;
+    // lseek(0, 0, SEEK_END);
+    read(0, &c, 1);
+    // tcflush(0, TCIFLUSH);
+    return c;
+}
+
+inline void clrscr() { printf("\x1B[2J\x1B[H"); }
+
+int main()
+{
+    termios old, current;
+    tcgetattr(0, &old);
+    current = old;
+    current.c_lflag &= ~ICANON & ~ECHO; /* disable buffered i/o, no echo */
+    current.c_cc[VMIN]  = 1;         // control chars (MIN value) = 1
+    current.c_cc[VTIME] = 0;         // control chars (TIME value) = 0 (No time)
+    tcsetattr(0, TCSANOW, &current);
+
+    std::copy(objects.begin(), objects.end(), std::back_inserter(scene));
+    std::copy(lights.begin(), lights.end(), std::back_inserter(scene));
+
+    Frame frame;
+    frame.data = new v3[W * H];
+
+    printf("\e[?1049h");
+
+    // choice_t times = 0;
+    // int samples = 0;
+
+    char c='\n';
+    do {
+	// clrscr();
+	gotoxy(0,0);
+
+	switch(c) {
+	case 65: // Up
+	    camera.rotation.x -= 0.08 * (0.5 / camera.zoom);
+	    break;
+	case 66: // Down
+	    camera.rotation.x += 0.08 * (0.5 / camera.zoom);
+	    break;
+        case 67: // Right
+	    camera.rotation.y -= 0.08 * (0.5 / camera.zoom);
+	    break;
+	case 68: // Left
+	    camera.rotation.y += 0.08 * (0.5 / camera.zoom);
+	    break;
+	case 'a':
+	    camera.translation.x -= 0.5;
+	    break;
+	case 'd':
+	    camera.translation.x += 0.5;;
+	    break;
+	case 'w':
+	    camera.translation.z += 0.5;;
+	    break;
+	case 's':
+	    camera.translation.z -= 0.5;;
+	    break;
+	case ' ':
+	    camera.translation.y -= 0.5;
+	    break;
+	case 'z':
+	    camera.translation.y += 0.5;
+	    break;
+	case '+':
+	    camera.exposure -= 0.015;
+	    break;
+	case '-':
+	    camera.exposure += 0.015;
+	    break;
+	case 'i':
+	    camera.zoom *= 1.015;
+	    break;
+	case 'o':
+	    camera.zoom /= 1.015;
+	    break;
+	case 'p':
+	    printf("pos: " V3_FMT "\nrot: " V3_FMT, V3_ARG(camera.pos), V3_ARG(camera.rotation));
+	    break;
+	case '\n':
+	case 'q':
+	    break;
+	default:
+	    continue;
+        }
+
+	// const auto before = std::chrono::high_resolution_clock::now();
+	render_frame(screen, camera, frame, scene, objects, lights);
+	// const auto after = std::chrono::high_resolution_clock::now();
+	// times += std::chrono::duration<choice_t, std::milli>(after - before).count();
+	// samples++;
+	display_frame(screen, frame);
+	//std::cout <<  1000/(times/samples) << " FPS" << std::endl;
+    } while((c = getch()) != 'q');
+
+    printf("\e[?1049l");
+    tcsetattr(0, TCSANOW, &old);
+    return 0;
+}

src/objects.cpp

@@ -0,0 +1,70 @@
+#include "../include/objects.hpp"
+
+/* Plane */
+
+Plane::Plane(const v3& n, choice_t offset, const Material& specs, int surface_type)
+    : Renderable(specs), n(n), offset(offset), surface_type(surface_type) {}
+Plane::Plane(const v3& n, choice_t offset, const Material& specs) : Plane(n, offset, specs, 2) {}
+
+int Plane::intersect(const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal)
+{
+    const choice_t d = dir ^ n;
+    if(d <= 0) return 0;
+
+    const choice_t k = -((eye ^ n) + offset) / d;
+
+    hit_loc = eye + dir * k;
+    hit_normal = n;
+    hit_dist = (hit_loc - eye).len();
+    return surface_type;
+}
+
+/* Sphere */
+
+Sphere::Sphere(const v3& origin, choice_t r, const Material& specs, int surface_type)
+    : Renderable(specs), o(origin), r(r), surface_type(surface_type) {}
+Sphere::Sphere(const v3& origin, choice_t r) : Sphere(origin, r, Material(), 1) {}
+
+template<int K = 1>
+static inline int _sphere_intersect(Sphere& self, const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal)
+{
+    const v3 d = self.o - eye;
+    const choice_t dist_above_o = dir ^ d;
+
+    // We're facing directly away from the sphere
+    if(dist_above_o < 0) { return -1; }
+
+    // The ray is pointing directly at the center of the sphere
+    // if(fabs(dist_above_o - d.len()) < 1e-5) {
+    // 	hit_loc = eye + d.norm() * (d.len() - self.r);
+    // 	hit_normal = self.normal(hit_loc);
+    // 	hit_dist = (hit_loc - eye).len();
+    // 	return self.surface_type;
+    // }
+
+    const choice_t rejection = sqrtf(d.len() * d.len() - dist_above_o * dist_above_o);
+    // The ray is tangent/missing the sphere
+    if(rejection >= self.r) return -1;
+
+    choice_t inside = sqrt(self.r*self.r - rejection*rejection) * K;
+
+    hit_dist = dist_above_o - inside;
+    hit_loc = eye + dir * hit_dist;
+    hit_normal = self.normal(hit_loc);
+    return self.surface_type;
+}
+
+int Sphere::intersect(const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal)
+{ return _sphere_intersect(*this, eye, dir, hit_dist, hit_loc, hit_normal); }
+
+int Sphere::transmit(const v3& eye, const v3& dir, choice_t& hit_dist, v3& hit_loc, v3& hit_normal)
+{ return _sphere_intersect<-1>(*this, eye, dir, hit_dist, hit_loc, hit_normal); }
+
+
+/* Light */
+
+Light::Light(const v3& origin, choice_t r, const v3& colour) : Light(origin, r,
+								   Material{ .specular = 0,
+									     .roughness = 1,
+									     .colour = colour }) {}
+Light::Light(const v3& origin, choice_t r, const Material& specs) : Sphere(origin, r, specs, 0) {}

src/raytrace.cpp

@@ -0,0 +1,164 @@
+#include "../include/raytrace.hpp"
+#include "../include/objects.hpp"
+#include <math.h>
+#include <random>
+
+Renderable::Renderable(const Material& specs) : specs(specs) {}
+
+Renderable* shoot_ray_into_objects(const v3& eye, const v3& dir, const std::vector<Renderable*>& objs,
+				   choice_t& hit_dist, v3& hit_loc, v3& hit_normal, int& hit_type,
+				   Renderable* self)
+{
+    hit_type = -1;
+    Renderable* hit_object = nullptr;
+
+    // Iterate across objects to find the closest hit
+    for(auto* r : objs) {
+	if(r == self) continue;
+
+	choice_t dist = 1e10;
+	v3 loc, normal;
+	int type = -1;
+
+	type = r->intersect(eye, dir, dist, loc, normal);
+
+	// Hit is closer than closest hit
+	if(type != -1 && dist < 1e9 && dist < hit_dist) {
+	    hit_dist = dist; hit_loc = loc;
+	    hit_normal = normal; hit_type = type;
+	    hit_object = r;
+	}
+    }
+
+    return hit_object;
+}
+
+v3 ray_trace(const v3& eye, const v3& dir,
+	     const std::vector<Renderable*>& scene,
+	     const std::vector<Renderable*>& objects,
+	     const std::vector<Renderable*>& lights, int depth,
+	     Renderable* skip, choice_t* dist_to_hit)
+{
+    choice_t hit_dist = 1e10;
+    v3 hit_loc, hit_normal;
+    int hit_type;
+    Renderable* hit_obj = shoot_ray_into_objects(eye, dir, scene, hit_dist, hit_loc, hit_normal, hit_type, skip);
+
+    if(dist_to_hit != nullptr) *dist_to_hit = hit_dist;
+
+    v3 colour = v3{0.05, 0.05, 0.05};
+
+    if (hit_type == 0) {
+	Light *l = (Light *)hit_obj;
+
+	const choice_t kDistance = powf((eye - l->o).norm() ^ hit_normal, 3);
+
+	colour += ((Light*)hit_obj)->specs.colour * 8 * kDistance;
+	if(depth > 0) colour += ray_trace(eye, dir, scene, objects, lights, depth - 1, l) * (1 - kDistance);
+    }  else if (hit_type > 0) {
+	// Figure out which lights are visible from the hit location
+	choice_t dist;
+	v3 loc, normal;
+	v3 light_dir;
+	int type;
+
+	Light *l = nullptr;
+
+	for (auto* _l : lights) {
+	    l = (Light *)_l;
+	    const v3 to_light = l->o - hit_loc;
+	    light_dir = to_light.norm();
+	    dist = to_light.len();
+	    shoot_ray_into_objects(hit_loc, light_dir, objects,
+				   dist, loc, normal, type, hit_obj);
+
+	    // We hit the light directly: shade the pixel.
+	    if (type == -1) {
+		const choice_t kDistance = to_light.lensquared() / 1000;
+		colour += cook_torrance(dir, hit_normal, light_dir, hit_obj->specs.specular,
+					hit_obj->specs.roughness, hit_obj->specs.ior,
+					hit_obj->specs.colour, l->specs.colour)
+		    * (1 - hit_obj->specs.transparency) / kDistance;
+	    }
+	}
+
+	// Add the colour from the reflection, if we have enough depth
+	const choice_t kReflectionAtt = hit_obj->specs.specular * (1 - hit_obj->specs.roughness);
+	if(kReflectionAtt > 1e-4 && depth > 0) {
+	    colour += ray_trace(hit_loc, -(dir).ref(hit_normal),
+				scene, objects, lights, depth-1)
+		* kReflectionAtt;
+	}
+
+        if (hit_obj->specs.transparency >= 1e-3 && depth > 0) {
+	    choice_t dist = 1e10; v3 loc, normal, light_dir;
+	    const v3 refract_dir = -dir.refract(hit_normal, 1/hit_obj->specs.ior);
+            if(hit_obj->transmit(hit_loc, refract_dir, dist, loc, normal) != -1) {
+		colour += ray_trace(loc, dir, scene, objects, lights, depth-1, hit_obj)
+		    * hit_obj->specs.transparency * (hit_obj->specs.colour + 0.75).norm();
+            }
+        }
+
+	// // Random sample
+	// std::random_device rd;
+	// std::mt19937 gen(rd());
+	// std::uniform_real_distribution<> dis(0.0, 1.0);
+	// const v3 y = hit_normal.cross(dir);
+	// for(int i = 0; i < 2; ++i) {
+	//     const v3 sample_dir = (hit_normal + dir * dis(gen) + y * dis(gen)).norm();
+	//     choice_t dist = 1;
+	//     const v3 hit_colour = ray_trace(hit_loc, sample_dir, scene, objects, lights, 0, hit_obj, &dist);
+	//     colour += hit_colour / (dist * dist);
+	// }
+    }
+
+    return colour;
+}
+
+static constexpr inline choice_t chi(choice_t x) { return x > 0 ? 1 : 0; }
+
+choice_t ggx_ndf(const v3& dir, const v3& h, const v3& normal, choice_t roughness)
+{
+    const choice_t alpha_sq = powf(roughness, 2);
+    const choice_t mn_sq = powf(-dir ^ normal, 2);
+    // return alpha_sq * chi(h ^ normal) / (M_PI * powf((1 - mn_sq * (1 + alpha_sq)), 2));
+    return alpha_sq * chi(h ^ normal) / (M_PI * powf(mn_sq * (alpha_sq + (1 - mn_sq) / mn_sq), 2));
+}
+
+choice_t blinn_phong_ndf(const v3& dir, const v3& h, const v3& normal, choice_t roughness)
+{
+    const choice_t alpha_sq = powf(roughness, 4);
+    return 1 / (M_PI * alpha_sq) * powf(normal ^ h, 2 / alpha_sq - 2);
+}
+
+choice_t ggx_gaf_partial(const v3& dir, const v3& normal, const v3& light)
+{
+
+}
+
+choice_t ggx_gaf(const v3& dir, const v3& normal, const v3& light)
+{
+
+}
+
+choice_t cook_torrance_gaf(const v3& dir, const v3& h, const v3& normal, const v3& light)
+{
+    const choice_t fac = 2 * (h ^ normal) / (-dir ^ h);
+    return fmin(1, fmin(fac * (-dir ^ normal), fac * (light ^ normal)));
+}
+
+choice_t cook_torrance_fresnel(const v3& dir, const v3& h, choice_t ior)
+{
+    const choice_t f_naught = powf(ior - 1, 2) / powf(ior + 1, 2);
+    return f_naught + (1 - f_naught) * powf(1 - (-dir ^ h), 5);
+}
+
+v3 cook_torrance(const v3& dir, const v3& normal, const v3& light, choice_t specular, choice_t roughness, choice_t ior,
+		 const v3& light_colour, const v3& surface_colour)
+{
+    const v3 h = (light - dir).norm();
+    const choice_t r_spec = ggx_ndf(dir, h, light, roughness) *
+	cook_torrance_gaf(dir, h, normal, light) * cook_torrance_fresnel(dir, h, ior) /
+	(4 * (normal ^ light) * (-dir ^ normal));
+    return (surface_colour * (1 - specular) / M_PI + light_colour * (specular * r_spec)) * fabs(normal ^ light);
+}