ray: introduce a rudimentary ray tracer

My first ray tracer, it only has spheres, planes, and point light sources.

No texture mapping, no soft shadows, no global illumination.

This is all very basic right now, the camera movement is simple and boring, but
sufficient for further development and optimization.

I made some effort to support multiple CPUs, it should detect the number of
CPUs in the system and use enough pthreads to keep them busy.

Jacco Bikker's tutorial on flipcode was the original impetus to do this, and
definitely served as a guide early on.

21 files changed, 1279 insertions, 0 deletions

diff --git a/modules/ray/ray.c b/modules/ray/ray.c
new file mode 100644
index 0000000..60d08cf
--- /dev/null
+++ b/modules/ray/ray.c
@@ -0,0 +1,161 @@
+#include <stdint.h>
+#include <inttypes.h>
+#include <math.h>
+
+#include "fb.h"
+#include "rototiller.h"
+#include "util.h"
+
+#include "ray_camera.h"
+#include "ray_object.h"
+#include "ray_scene.h"
+#include "ray_threads.h"
+
+/* Copyright (C) 2016 Vito Caputo <vcaputo@pengaru.com> */
+
+/* ray trace a simple scene into the fragment */
+static void ray(fb_fragment_t *fragment)
+{
+	static ray_object_t	objects[] = {
+		{
+			.plane = {
+				.type = RAY_OBJECT_TYPE_PLANE,
+				.surface = {
+					.color = { .x = 0.4, .y = 0.2, .z = 0.5 },
+					.diffuse = 1.0f,
+					.specular = 0.2f,
+				},
+				.normal = { .x = 0.0, .y = 1.0, .z = 0.0 },
+				.distance = -3.2f,
+			}
+		}, {
+			.sphere = {
+				.type = RAY_OBJECT_TYPE_SPHERE,
+				.surface = {
+					.color = { .x = 1.0, .y = 0.0, .z = 0.0 },
+					.diffuse = 1.0f,
+					.specular = 0.05f,
+				},
+				.center = { .x = 0.5, .y = 1.0, .z = 0.0 },
+				.radius = 1.2f,
+			}
+		}, {
+			.sphere = {
+				.type = RAY_OBJECT_TYPE_SPHERE,
+				.surface = {
+					.color = { .x = 0.0, .y = 0.0, .z = 1.0 },
+					.diffuse = 1.0f,
+					.specular = 1.0f,
+				},
+				.center = { .x = -2.0, .y = 1.0, .z = 0.0 },
+				.radius = 0.9f,
+			}
+		}, {
+			.sphere = {
+				.type = RAY_OBJECT_TYPE_SPHERE,
+				.surface = {
+					.color = { .x = 0.0, .y = 1.0, .z = 1.0 },
+					.diffuse = 1.0f,
+					.specular = 1.0f,
+				},
+				.center = { .x = 2.0, .y = -1.0, .z = 0.0 },
+				.radius = 1.0f,
+			}
+		}, {
+			.sphere = {
+				.type = RAY_OBJECT_TYPE_SPHERE,
+				.surface = {
+					.color = { .x = 0.0, .y = 1.0, .z = 0.0 },
+					.diffuse = 1.0f,
+					.specular = 1.0f,
+				},
+				.center = { .x = 0.2, .y = -1.25, .z = 0.0 },
+				.radius = 0.6f,
+			}
+		}, {
+			.light = {
+				.type = RAY_OBJECT_TYPE_LIGHT,
+				.brightness = 1.0,
+				.emitter = {
+					.point.type = RAY_LIGHT_EMITTER_TYPE_POINT,
+					.point.center = { .x = 3.0f, .y = 3.0f, .z = 3.0f },
+					.point.surface = {
+						.color = { .x = 1.0f, .y = 1.0f, .z = 1.0f },
+					},
+				}
+			}
+		}
+	};
+
+	ray_camera_t		camera = {
+					.position = { .x = 0.0, .y = 0.0, .z = 6.0 },
+					.orientation = {
+						.yaw = RAY_EULER_DEGREES(0.0f),
+						.pitch = RAY_EULER_DEGREES(0.0f),
+						.roll = RAY_EULER_DEGREES(180.0f),
+					},
+					.focal_length = 700.0f,
+					.width = fragment->width,
+					.height = fragment->height,
+				};
+
+	static ray_scene_t	scene = {
+					.objects = objects,
+					.n_objects = nelems(objects),
+					.lights = &objects[5],
+					.n_lights = 1,
+					.ambient_color = { .x = 1.0f, .y = 1.0f, .z = 1.0f },
+					.ambient_brightness = .04f,
+				};
+	static int		initialized;
+	static ray_threads_t	*threads;
+	static fb_fragment_t	*fragments;
+	static unsigned		ncpus;
+#if 1
+	/* animated point light source */
+	static double		r;
+
+	r += .02;
+
+	scene.lights[0].light.emitter.point.center.x = cosf(r) * 3.5f;
+	scene.lights[0].light.emitter.point.center.z = sinf(r) * 3.5f;
+	camera.orientation.yaw = sinf(r) / 4;
+	camera.orientation.pitch = sinf(r * 10) / 100;
+	camera.orientation.roll = RAY_EULER_DEGREES(180.0f) + cosf(r) / 10;
+	camera.position.x = cosf(r) / 10;
+	camera.position.z = 4.0f + sinf(r) / 10;
+#endif
+
+	if (!initialized) {
+		initialized = 1;
+		ncpus = get_ncpus();
+
+		if (ncpus > 1) {
+			threads = ray_threads_create(ncpus - 1);
+			fragments = malloc(sizeof(fb_fragment_t) * ncpus);
+		}
+	}
+
+	if (ncpus > 1) {
+		/* Always recompute the fragments[] geometry.
+		 * This way the fragment geometry can change at any moment and things will
+		 * continue functioning, which may prove important later on.
+		 * (imagine things like a preview window, or perhaps composite modules
+		 * which call on other modules supplying virtual fragments of varying dimensions..)
+		 */
+		fb_fragment_divide(fragment, ncpus, fragments);
+	} else {
+		fragments = fragment;
+	}
+
+	ray_scene_render_fragments(&scene, &camera, threads, fragments);
+}
+
+
+rototiller_renderer_t	ray_renderer = {
+	.render = ray,
+	.name = "ray",
+	.description = "Multi-threaded ray tracer",
+	.author = "Vito Caputo <vcaputo@pengaru.com>",
+	.license = "GPLv2",
+};
diff --git a/modules/ray/ray.h b/modules/ray/ray.h
new file mode 100644
index 0000000..d33f96a
--- /dev/null
+++ b/modules/ray/ray.h
@@ -0,0 +1,8 @@
+#ifndef _RAY_RAY_H
+#define _RAY_RAY_H
+
+#include "fb.h"
+
+void ray(fb_fragment_t *fragment);
+
+#endif
diff --git a/modules/ray/ray_3f.h b/modules/ray/ray_3f.h
new file mode 100644
index 0000000..8408abb
--- /dev/null
+++ b/modules/ray/ray_3f.h
@@ -0,0 +1,161 @@
+#ifndef _RAY_3F_H
+#define _RAY_3F_H
+
+#include <math.h>
+
+typedef struct ray_3f_t {
+	float	x, y, z;
+} ray_3f_t;
+
+
+/* return the result of (a + b) */
+static inline ray_3f_t ray_3f_add(ray_3f_t *a, ray_3f_t *b)
+{
+	ray_3f_t	res = {
+				.x = a->x + b->x,
+				.y = a->y + b->y,
+				.z = a->z + b->z,
+			};
+
+	return res;
+}
+
+
+/* return the result of (a - b) */
+static inline ray_3f_t ray_3f_sub(ray_3f_t *a, ray_3f_t *b)
+{
+	ray_3f_t	res = {
+				.x = a->x - b->x,
+				.y = a->y - b->y,
+				.z = a->z - b->z,
+			};
+
+	return res;
+}
+
+
+/* return the result of (-v) */
+static inline ray_3f_t ray_3f_negate(ray_3f_t *v)
+{
+	ray_3f_t	res = {
+				.x = -v->x,
+				.y = -v->y,
+				.z = -v->z,
+			};
+
+	return res;
+}
+
+
+/* return the result of (a * b) */
+static inline ray_3f_t ray_3f_mult(ray_3f_t *a, ray_3f_t *b)
+{
+	ray_3f_t	res = {
+				.x = a->x * b->x,
+				.y = a->y * b->y,
+				.z = a->z * b->z,
+			};
+
+	return res;
+}
+
+
+/* return the result of (v * scalar) */
+static inline ray_3f_t ray_3f_mult_scalar(ray_3f_t *v, float scalar)
+{
+	ray_3f_t	res = {
+				.x = v->x * scalar,
+				.y = v->y * scalar,
+				.z = v->z * scalar,
+			};
+
+	return res;
+}
+
+
+/* return the result of (uv / scalar) */
+static inline ray_3f_t ray_3f_div_scalar(ray_3f_t *v, float scalar)
+{
+	ray_3f_t	res = {
+				.x = v->x / scalar,
+				.y = v->y / scalar,
+				.z = v->z / scalar,
+			};
+
+	return res;
+}
+
+
+/* return the result of (a . b) */
+static inline float ray_3f_dot(ray_3f_t *a, ray_3f_t *b)
+{
+	return a->x * b->x + a->y * b->y + a->z * b->z;
+}
+
+
+/* return the length of the supplied vector */
+static inline float ray_3f_length(ray_3f_t *v)
+{
+	return sqrtf(ray_3f_dot(v, v));
+}
+
+
+/* return the normalized form of the supplied vector */
+static inline ray_3f_t ray_3f_normalize(ray_3f_t *v)
+{
+	ray_3f_t	nv;
+	float		f;
+
+	f = 1.0f / ray_3f_length(v);
+
+	nv.x = f * v->x;
+	nv.y = f * v->y;
+	nv.z = f * v->z;
+
+	return nv;
+}
+
+
+/* return the distance between two arbitrary points */
+static inline float ray_3f_distance(ray_3f_t *a, ray_3f_t *b)
+{
+	return sqrtf(powf(a->x - b->x, 2) + powf(a->y - b->y, 2) + powf(a->z - b->z, 2));
+}
+
+
+/* return the cross product of two unit vectors */
+static inline ray_3f_t ray_3f_cross(ray_3f_t *a, ray_3f_t *b)
+{
+	ray_3f_t	product;
+
+	product.x = a->y * b->z - a->z * b->y;
+	product.y = a->z * b->x - a->x * b->z;
+	product.z = a->x * b->y - a->y * b->x;
+
+	return product;
+}
+
+
+/* return the linearly interpolated vector between the two vectors at point alpha (0-1.0) */
+static inline ray_3f_t ray_3f_lerp(ray_3f_t *a, ray_3f_t *b, float alpha)
+{
+	ray_3f_t	lerp_a, lerp_b;
+
+	lerp_a = ray_3f_mult_scalar(a, 1.0f - alpha);
+	lerp_b = ray_3f_mult_scalar(b, alpha);
+
+	return ray_3f_add(&lerp_a, &lerp_b);
+}
+
+
+/* return the normalized linearly interpolated vector between the two vectors at point alpha (0-1.0) */
+static inline ray_3f_t ray_3f_nlerp(ray_3f_t *a, ray_3f_t *b, float alpha)
+{
+	ray_3f_t	lerp;
+
+	lerp = ray_3f_lerp(a, b, alpha);
+
+	return ray_3f_normalize(&lerp);
+}
+
+#endif
diff --git a/modules/ray/ray_camera.c b/modules/ray/ray_camera.c
new file mode 100644
index 0000000..0703c2e
--- /dev/null
+++ b/modules/ray/ray_camera.c
@@ -0,0 +1,85 @@
+#include "fb.h"
+
+#include "ray_camera.h"
+#include "ray_euler.h"
+
+
+/* Produce a vector from the provided orientation vectors and proportions. */
+static ray_3f_t project_corner(ray_3f_t *forward, ray_3f_t *left, ray_3f_t *up, float focal_length, float horiz, float vert)
+{
+	ray_3f_t	tmp;
+	ray_3f_t	corner;
+
+	corner = ray_3f_mult_scalar(forward, focal_length);
+	tmp = ray_3f_mult_scalar(left, horiz);
+	corner = ray_3f_add(&corner, &tmp);
+	tmp = ray_3f_mult_scalar(up, vert);
+	corner = ray_3f_add(&corner, &tmp);
+
+	return ray_3f_normalize(&corner);
+}
+
+
+/* Produce vectors for the corners of the entire camera frame, used for interpolation. */
+static void project_corners(ray_camera_t *camera, ray_camera_frame_t *frame)
+{
+	ray_3f_t	forward, left, up, right, down;
+	float		half_horiz = (float)camera->width / 2.0f;
+	float		half_vert = (float)camera->height / 2.0f;
+
+	ray_euler_basis(&camera->orientation, &forward, &up, &left);
+	right = ray_3f_negate(&left);
+	down = ray_3f_negate(&up);
+
+	frame->nw = project_corner(&forward, &left, &up, camera->focal_length, half_horiz,  half_vert);
+	frame->ne = project_corner(&forward, &right, &up, camera->focal_length, half_horiz,  half_vert);
+	frame->se = project_corner(&forward, &right, &down, camera->focal_length, half_horiz,  half_vert);
+	frame->sw = project_corner(&forward, &left, &down, camera->focal_length, half_horiz,  half_vert);
+}
+
+
+/* Begin a frame for the fragment of camera projection, initializing frame and ray. */
+void ray_camera_frame_begin(ray_camera_t *camera, fb_fragment_t *fragment, ray_ray_t *ray, ray_camera_frame_t *frame)
+{
+	/* References are kept to the camera, fragment, and ray to be traced.
+	 * The ray is maintained as we step through the frame, that is the
+	 * purpose of this api.
+	 *
+	 * Since the ray direction should be a normalized vector, the obvious
+	 * implementation is a bit costly.  The camera frame api hides this
+	 * detail so we can explore interpolation techniques to potentially
+	 * lessen the per-pixel cost.
+	 */
+	frame->camera = camera;
+	frame->fragment = fragment;
+	frame->ray = ray;
+
+	frame->x = frame->y = 0;
+
+	/* From camera->orientation and camera->focal_length compute the vectors
+	 * through the viewport's corners, and place these normalized vectors
+	 * in frame->(nw,ne,sw,se).
+	 *
+	 * These can than be interpolated between to produce the ray vectors
+	 * throughout the frame's fragment.  The efficient option of linear
+	 * interpolation will not maintain the unit vector length, so to
+	 * produce normalized interpolated directions will require the costly
+	 * normalize function.
+	 *
+	 * I'm hoping a simple length correction table can be used to fixup the
+	 * linearly interpolated vectors to make them unit vectors with just
+	 * scalar multiplication instead of the sqrt of normalize.
+	 */
+	project_corners(camera, frame);
+
+	frame->x_delta = 1.0f / (float)camera->width;
+	frame->y_delta = 1.0f / (float)camera->height;
+	frame->x_alpha = frame->x_delta * (float)fragment->x;
+	frame->y_alpha = frame->y_delta * (float)fragment->y;
+
+	frame->cur_w = ray_3f_nlerp(&frame->nw, &frame->sw, frame->y_alpha);
+	frame->cur_e = ray_3f_nlerp(&frame->ne, &frame->se, frame->y_alpha);
+
+	ray->origin = camera->position;
+	ray->direction = frame->cur_w;
+}
diff --git a/modules/ray/ray_camera.h b/modules/ray/ray_camera.h
new file mode 100644
index 0000000..387f8c5
--- /dev/null
+++ b/modules/ray/ray_camera.h
@@ -0,0 +1,77 @@
+#ifndef _RAY_CAMERA_H
+#define _RAY_CAMERA_H
+
+#include <math.h>
+
+#include "fb.h"
+
+#include "ray_3f.h"
+#include "ray_euler.h"
+#include "ray_ray.h"
+
+
+typedef struct ray_camera_t {
+	ray_3f_t	position;		/* position of camera, the origin of all its rays */
+	ray_euler_t	orientation;		/* orientation of the camera */
+	float		focal_length;		/* controls the field of view */
+	unsigned	width;			/* width of camera viewport in pixels */
+	unsigned	height;			/* height of camera viewport in pixels */
+} ray_camera_t;
+
+
+typedef struct ray_camera_frame_t {
+	ray_camera_t	*camera;		/* the camera supplied to frame_begin() */
+	fb_fragment_t	*fragment;		/* the fragment supplied to frame_begin() */
+	ray_ray_t	*ray;			/* the ray supplied to frame_begin(), which gets updated as we step through the frame. */
+
+	ray_3f_t	nw, ne, sw, se;		/* directions pointing through the corners of the frame fragment */
+	ray_3f_t	cur_w, cur_e;		/* current row's west and east ends */
+	float		x_alpha, y_alpha;	/* interpolation position along the x and y axis */
+	float		x_delta, y_delta;	/* interpolation step delta along the x and y axis */
+	unsigned	x, y;			/* integral position within frame fragment */
+} ray_camera_frame_t;
+
+
+void ray_camera_frame_begin(ray_camera_t *camera, fb_fragment_t *fragment, ray_ray_t *ray, ray_camera_frame_t *frame);
+
+
+/* Step the ray through the frame on the x axis, returns 1 when rays remain on this axis, 0 at the end. */
+/* When 1 is returned, frame->ray is left pointing through the new coordinate. */
+static inline int ray_camera_frame_x_step(ray_camera_frame_t *frame)
+{
+	frame->x++;
+
+	if (frame->x >= frame->fragment->width) {
+		frame->x = 0;
+		frame->x_alpha = frame->x_delta * (float)frame->fragment->x;
+		return 0;
+	}
+
+	frame->x_alpha += frame->x_delta;
+	frame->ray->direction = ray_3f_nlerp(&frame->cur_w, &frame->cur_e, frame->x_alpha);
+
+	return 1;
+}
+
+
+/* Step the ray through the frame on the y axis, returns 1 when rays remain on this axis, 0 at the end. */
+/* When 1 is returned, frame->ray is left pointing through the new coordinate. */
+static inline int ray_camera_frame_y_step(ray_camera_frame_t *frame)
+{
+	frame->y++;
+
+	if (frame->y >= frame->fragment->height) {
+		frame->y = 0;
+		frame->y_alpha = frame->y_delta * (float)frame->fragment->y;
+		return 0;
+	}
+
+	frame->y_alpha += frame->y_delta;
+	frame->cur_w = ray_3f_nlerp(&frame->nw, &frame->sw, frame->y_alpha);
+	frame->cur_e = ray_3f_nlerp(&frame->ne, &frame->se, frame->y_alpha);
+	frame->ray->direction = frame->cur_w;
+
+	return 1;
+}
+
+#endif
diff --git a/modules/ray/ray_color.h b/modules/ray/ray_color.h
new file mode 100644
index 0000000..9fe62c1
--- /dev/null
+++ b/modules/ray/ray_color.h
@@ -0,0 +1,29 @@
+#ifndef _RAY_COLOR_H
+#define _RAY_COLOR_H
+
+#include <stdint.h>
+
+#include "ray_3f.h"
+
+typedef ray_3f_t ray_color_t;
+
+/* convert a vector into a packed, 32-bit rgb pixel value */
+static inline uint32_t ray_color_to_uint32_rgb(ray_color_t color) {
+	uint32_t	pixel;
+
+	/* doing this all per-pixel, ugh. */
+
+	if (color.x > 1.0f) color.x = 1.0f;
+	if (color.y > 1.0f) color.y = 1.0f;
+	if (color.z > 1.0f) color.z = 1.0f;
+
+	pixel = (uint32_t)(color.x * 255.0f);
+	pixel <<= 8;
+	pixel |= (uint32_t)(color.y * 255.0f);
+	pixel <<= 8;
+	pixel |= (uint32_t)(color.z * 255.0f);
+
+	return pixel;
+}
+
+#endif
diff --git a/modules/ray/ray_euler.h b/modules/ray/ray_euler.h
new file mode 100644
index 0000000..86f5221
--- /dev/null
+++ b/modules/ray/ray_euler.h
@@ -0,0 +1,45 @@
+#ifndef _RAY_EULER_H
+#define _RAY_EULER_H
+
+#include <math.h>
+
+#include "ray_3f.h"
+
+
+/* euler angles are convenient for describing orientation */
+typedef struct ray_euler_t {
+	float	pitch;	/* pitch in radiasn */
+	float	yaw;	/* yaw in radians */
+	float	roll;	/* roll in radians */
+} ray_euler_t;
+
+
+/* convenience macro for converting degrees to radians */
+#define RAY_EULER_DEGREES(_deg) \
+	(_deg * (2 * M_PI / 360.0f))
+
+
+/* produce basis vectors from euler angles */
+static inline void ray_euler_basis(ray_euler_t *e, ray_3f_t *forward, ray_3f_t *up, ray_3f_t *left)
+{
+	float	cos_yaw = cosf(e->yaw);
+	float	sin_yaw = sinf(e->yaw);
+	float	cos_roll = cosf(e->roll);
+	float	sin_roll = sinf(e->roll);
+	float	cos_pitch = cosf(e->pitch);
+	float	sin_pitch = sinf(e->pitch);
+
+	forward->x = sin_yaw;
+	forward->y = -sin_pitch * cos_yaw;
+	forward->z = cos_pitch * cos_yaw;
+
+	up->x = -cos_yaw * sin_roll;
+	up->y = -sin_pitch * sin_yaw * sin_roll + cos_pitch * cos_roll;
+	up->z = cos_pitch * sin_yaw * sin_roll + sin_pitch * cos_roll;
+
+	left->x = cos_yaw * cos_roll;
+	left->y = sin_pitch * sin_yaw * cos_roll + cos_pitch * sin_roll;
+	left->z = -cos_pitch * sin_yaw * cos_roll + sin_pitch * sin_roll;
+}
+
+#endif
diff --git a/modules/ray/ray_light_emitter.h b/modules/ray/ray_light_emitter.h
new file mode 100644
index 0000000..3b5509e
--- /dev/null
+++ b/modules/ray/ray_light_emitter.h
@@ -0,0 +1,18 @@
+#ifndef _RAY_LIGHT_EMITTER_H
+#define _RAY_LIGHT_EMITTER_H
+
+#include "ray_object_point.h"
+#include "ray_object_sphere.h"
+
+typedef enum ray_light_emitter_type_t {
+	RAY_LIGHT_EMITTER_TYPE_SPHERE,
+	RAY_LIGHT_EMITTER_TYPE_POINT,
+} ray_light_emitter_type_t;
+
+typedef union ray_light_emitter_t {
+	ray_light_emitter_type_t	type;
+	ray_object_sphere_t		sphere;
+	ray_object_point_t		point;
+} ray_light_emitter_t;
+
+#endif
diff --git a/modules/ray/ray_object.c b/modules/ray/ray_object.c
new file mode 100644
index 0000000..8b324e5
--- /dev/null
+++ b/modules/ray/ray_object.c
@@ -0,0 +1,67 @@
+#include "ray_object.h"
+#include "ray_object_light.h"
+#include "ray_object_plane.h"
+#include "ray_object_point.h"
+#include "ray_object_sphere.h"
+#include "ray_ray.h"
+#include "ray_scene.h"
+#include "ray_surface.h"
+
+
+/* Determine if a ray intersects object.
+ * If the object is intersected, store where along the ray the intersection occurs in res_distance.
+ */
+int ray_object_intersects_ray(ray_object_t *object, ray_ray_t *ray, float *res_distance)
+{
+	switch (object->type) {
+	case RAY_OBJECT_TYPE_SPHERE:
+		return ray_object_sphere_intersects_ray(&object->sphere, ray, res_distance);
+
+	case RAY_OBJECT_TYPE_POINT:
+		return ray_object_point_intersects_ray(&object->point, ray, res_distance);
+
+	case RAY_OBJECT_TYPE_PLANE:
+		return ray_object_plane_intersects_ray(&object->plane, ray, res_distance);
+
+	case RAY_OBJECT_TYPE_LIGHT:
+		return ray_object_light_intersects_ray(&object->light, ray, res_distance);
+	}
+}
+
+
+/* Return the surface normal of object @ point */
+ray_3f_t ray_object_normal(ray_object_t *object, ray_3f_t *point)
+{
+	switch (object->type) {
+	case RAY_OBJECT_TYPE_SPHERE:
+		return ray_object_sphere_normal(&object->sphere, point);
+
+	case RAY_OBJECT_TYPE_POINT:
+		return ray_object_point_normal(&object->point, point);
+
+	case RAY_OBJECT_TYPE_PLANE:
+		return ray_object_plane_normal(&object->plane, point);
+
+	case RAY_OBJECT_TYPE_LIGHT:
+		return ray_object_light_normal(&object->light, point);
+	}
+}
+
+
+/* Return the surface of object @ point */
+ray_surface_t ray_object_surface(ray_object_t *object, ray_3f_t *point)
+{
+	switch (object->type) {
+	case RAY_OBJECT_TYPE_SPHERE:
+		return ray_object_sphere_surface(&object->sphere, point);
+
+	case RAY_OBJECT_TYPE_POINT:
+		return ray_object_point_surface(&object->point, point);
+
+	case RAY_OBJECT_TYPE_PLANE:
+		return ray_object_plane_surface(&object->plane, point);
+
+	case RAY_OBJECT_TYPE_LIGHT:
+		return ray_object_light_surface(&object->light, point);
+	}
+}
diff --git a/modules/ray/ray_object.h b/modules/ray/ray_object.h
new file mode 100644
index 0000000..3dc27d1
--- /dev/null
+++ b/modules/ray/ray_object.h
@@ -0,0 +1,25 @@
+#ifndef _RAY_OBJECT_H
+#define _RAY_OBJECT_H
+
+#include "ray_object_light.h"
+#include "ray_object_plane.h"
+#include "ray_object_point.h"
+#include "ray_object_sphere.h"
+#include "ray_object_type.h"
+#include "ray_ray.h"
+#include "ray_scene.h"
+#include "ray_surface.h"
+
+typedef union ray_object_t {
+	ray_object_type_t	type;
+	ray_object_sphere_t	sphere;
+	ray_object_point_t	point;
+	ray_object_plane_t	plane;
+	ray_object_light_t	light;
+} ray_object_t;
+
+int ray_object_intersects_ray(ray_object_t *object, ray_ray_t *ray, float *res_distance);
+ray_3f_t ray_object_normal(ray_object_t *object, ray_3f_t *point);
+ray_surface_t ray_object_surface(ray_object_t *object, ray_3f_t *point);
+
+#endif
diff --git a/modules/ray/ray_object_light.h b/modules/ray/ray_object_light.h
new file mode 100644
index 0000000..73cf917
--- /dev/null
+++ b/modules/ray/ray_object_light.h
@@ -0,0 +1,60 @@
+#ifndef _RAY_OBJECT_LIGHT_H
+#define _RAY_OBJECT_LIGHT_H
+
+#include "ray_light_emitter.h"
+#include "ray_object_light.h"
+#include "ray_object_point.h"
+#include "ray_object_sphere.h"
+#include "ray_object_type.h"
+#include "ray_ray.h"
+#include "ray_scene.h"
+#include "ray_surface.h"
+
+
+typedef struct ray_object_light_t {
+	ray_object_type_t	type;
+	float			brightness;
+	ray_light_emitter_t	emitter;
+} ray_object_light_t;
+
+
+/* TODO: point is really the only one I've implemented... */
+static inline int ray_object_light_intersects_ray(ray_object_light_t *light, ray_ray_t *ray, float *res_distance)
+{
+	switch (light->emitter.type) {
+	case RAY_LIGHT_EMITTER_TYPE_POINT:
+		return ray_object_point_intersects_ray(&light->emitter.point, ray, res_distance);
+
+	case RAY_LIGHT_EMITTER_TYPE_SPHERE:
+		return ray_object_sphere_intersects_ray(&light->emitter.sphere, ray, res_distance);
+	}
+}
+
+
+static inline ray_3f_t ray_object_light_normal(ray_object_light_t *light, ray_3f_t *point)
+{
+	ray_3f_t	normal;
+
+	/* TODO */
+	switch (light->emitter.type) {
+	case RAY_LIGHT_EMITTER_TYPE_SPHERE:
+		return normal;
+
+	case RAY_LIGHT_EMITTER_TYPE_POINT:
+		return normal;
+	}
+}
+
+
+static inline ray_surface_t ray_object_light_surface(ray_object_light_t *light, ray_3f_t *point)
+{
+	switch (light->emitter.type) {
+	case RAY_LIGHT_EMITTER_TYPE_SPHERE:
+		return ray_object_sphere_surface(&light->emitter.sphere, point);
+
+	case RAY_LIGHT_EMITTER_TYPE_POINT:
+		return ray_object_point_surface(&light->emitter.point, point);
+	}
+}
+
+#endif
diff --git a/modules/ray/ray_object_plane.h b/modules/ray/ray_object_plane.h
new file mode 100644
index 0000000..96f204e
--- /dev/null
+++ b/modules/ray/ray_object_plane.h
@@ -0,0 +1,47 @@
+#ifndef _RAY_OBJECT_PLANE_H
+#define _RAY_OBJECT_PLANE_H
+
+#include "ray_object_type.h"
+#include "ray_ray.h"
+#include "ray_scene.h"
+#include "ray_surface.h"
+
+
+typedef struct ray_object_plane_t {
+	ray_object_type_t	type;
+	ray_surface_t		surface;
+	ray_3f_t		normal;
+	float			distance;
+} ray_object_plane_t;
+
+
+static inline int ray_object_plane_intersects_ray(ray_object_plane_t *plane, ray_ray_t *ray, float *res_distance)
+{
+	float	d = ray_3f_dot(&plane->normal, &ray->direction);
+
+	if (d != 0) {
+		float	distance = -(ray_3f_dot(&plane->normal, &ray->origin) + plane->distance) / d;
+
+		if (distance > 0) {
+			*res_distance = distance;
+
+			return 1;
+		}
+	}
+
+	return 0;
+}
+
+
+static inline ray_3f_t ray_object_plane_normal(ray_object_plane_t *plane, ray_3f_t *point)
+{
+	return plane->normal;
+}
+
+
+static inline ray_surface_t ray_object_plane_surface(ray_object_plane_t *plane, ray_3f_t *point)
+{
+	return plane->surface;
+}
+
+#endif
diff --git a/modules/ray/ray_object_point.h b/modules/ray/ray_object_point.h
new file mode 100644
index 0000000..5c83a68
--- /dev/null
+++ b/modules/ray/ray_object_point.h
@@ -0,0 +1,38 @@
+#ifndef _RAY_OBJECT_POINT_H
+#define _RAY_OBJECT_POINT_H
+
+#include "ray_3f.h"
+#include "ray_object_type.h"
+#include "ray_ray.h"
+#include "ray_scene.h"
+#include "ray_surface.h"
+
+
+typedef struct ray_object_point_t {
+	ray_object_type_t	type;
+	ray_surface_t		surface;
+	ray_3f_t		center;
+} ray_object_point_t;
+
+
+static inline int ray_object_point_intersects_ray(ray_object_point_t *point, ray_ray_t *ray, float *res_distance)
+{
+	/* TODO: determine a ray:point intersection */
+	return 0;
+}
+
+
+static inline ray_3f_t ray_object_point_normal(ray_object_point_t *point, ray_3f_t *_point)
+{
+	ray_3f_t	normal;
+
+	return normal;
+}
+
+
+static inline ray_surface_t ray_object_point_surface(ray_object_point_t *point, ray_3f_t *_point)
+{
+	return point->surface;
+}
+
+#endif
diff --git a/modules/ray/ray_object_sphere.h b/modules/ray/ray_object_sphere.h
new file mode 100644
index 0000000..6786bb7
--- /dev/null
+++ b/modules/ray/ray_object_sphere.h
@@ -0,0 +1,66 @@
+#ifndef _RAY_OBJECT_SPHERE_H
+#define _RAY_OBJECT_SPHERE_H
+
+#include <math.h>
+#include <stdio.h>
+
+#include "ray_3f.h"
+#include "ray_color.h"
+#include "ray_object_type.h"
+#include "ray_ray.h"
+#include "ray_scene.h"
+#include "ray_surface.h"
+
+
+typedef struct ray_object_sphere_t {
+	ray_object_type_t	type;
+	ray_surface_t		surface;
+	ray_3f_t		center;
+	float			radius;
+} ray_object_sphere_t;
+
+
+static inline int ray_object_sphere_intersects_ray(ray_object_sphere_t *sphere, ray_ray_t *ray, float *res_distance)
+{
+	ray_3f_t	v = ray_3f_sub(&ray->origin, &sphere->center);
+	float		b = ray_3f_dot(&v, &ray->direction);
+	float		disc = (sphere->radius * sphere->radius) - ray_3f_dot(&v, &v) + (b * b);
+
+	if (disc > 0) {
+		float	i1, i2;
+
+		disc = sqrtf(disc);
+
+		i1 = b - disc;
+		i2 = b + disc;
+
+		if (i2 > 0 && i1 > 0) {
+			*res_distance = i1;
+			return 1;
+		}
+	}
+
+	return 0;
+}
+
+
+/* return the normal of the surface at the specified point */
+static inline ray_3f_t ray_object_sphere_normal(ray_object_sphere_t *sphere, ray_3f_t *point)
+{
+	ray_3f_t	normal;
+	
+	normal = ray_3f_sub(point, &sphere->center);
+	normal = ray_3f_div_scalar(&normal, sphere->radius);	/* normalize without the sqrt() */
+
+	return normal;
+}
+
+
+/* return the surface of the sphere @ point */
+static inline ray_surface_t ray_object_sphere_surface(ray_object_sphere_t *sphere, ray_3f_t *point)
+{
+	/* uniform solids for now... */
+	return sphere->surface;
+}
+
+#endif
diff --git a/modules/ray/ray_object_type.h b/modules/ray/ray_object_type.h
new file mode 100644
index 0000000..6ce20f5
--- /dev/null
+++ b/modules/ray/ray_object_type.h
@@ -0,0 +1,11 @@
+#ifndef _RAY_OBJECT_TYPE_H
+#define _RAY_OBJECT_TYPE_H
+
+typedef enum ray_object_type_t {
+	RAY_OBJECT_TYPE_SPHERE,
+	RAY_OBJECT_TYPE_POINT,
+	RAY_OBJECT_TYPE_PLANE,
+	RAY_OBJECT_TYPE_LIGHT,
+} ray_object_type_t;
+
+#endif
diff --git a/modules/ray/ray_ray.h b/modules/ray/ray_ray.h
new file mode 100644
index 0000000..91469a2
--- /dev/null
+++ b/modules/ray/ray_ray.h
@@ -0,0 +1,11 @@
+#ifndef _RAY_RAY_H
+#define _RAY_RAY_H
+
+#include "ray_3f.h"
+
+typedef struct ray_ray_t {
+	ray_3f_t	origin;
+	ray_3f_t	direction;
+} ray_ray_t;
+
+#endif
diff --git a/modules/ray/ray_scene.c b/modules/ray/ray_scene.c
new file mode 100644
index 0000000..e44990b
--- /dev/null
+++ b/modules/ray/ray_scene.c
@@ -0,0 +1,188 @@
+#include <stdlib.h>
+#include <math.h>
+
+#include "fb.h"
+
+#include "ray_camera.h"
+#include "ray_color.h"
+#include "ray_object.h"
+#include "ray_ray.h"
+#include "ray_scene.h"
+#include "ray_threads.h"
+
+#define MAX_RECURSION_DEPTH	5
+
+
+static ray_color_t trace_ray(ray_scene_t *scene, ray_ray_t *ray, unsigned depth);
+
+
+/* Determine if the ray is obstructed by an object within the supplied distance, for shadows */
+static inline int ray_is_obstructed(ray_scene_t *scene, ray_ray_t *ray, float distance)
+{
+	unsigned	i;
+
+	for (i = 0; i < scene->n_objects; i++) {
+		float	ood;
+
+		if (scene->objects[i].type == RAY_OBJECT_TYPE_LIGHT)
+			continue;
+
+		if (ray_object_intersects_ray(&scene->objects[i], ray, &ood) &&
+		    ood < distance) {
+			return 1;
+		}
+	}
+
+	return 0;
+}
+
+
+/* Determine the color @ distance on ray on object viewed from origin */
+static inline ray_color_t shade_ray(ray_scene_t *scene, ray_ray_t *ray, ray_object_t *object, float distance, unsigned depth)
+{
+	ray_surface_t	surface;
+	ray_color_t	color;
+	ray_3f_t	rvec = ray_3f_mult_scalar(&ray->direction, distance);
+	ray_3f_t	intersection = ray_3f_sub(&ray->origin, &rvec);
+	ray_3f_t	normal = ray_object_normal(object, &intersection);
+	unsigned	i;
+
+	surface = ray_object_surface(object, &intersection);
+	color = ray_3f_mult_scalar(&scene->ambient_color, scene->ambient_brightness);
+	color = ray_3f_mult(&surface.color, &color);
+
+	/* visit lights for shadows and illumination */
+	for (i = 0; i < scene->n_lights; i++) {
+		ray_3f_t	lvec = ray_3f_sub(&scene->lights[i].light.emitter.point.center, &intersection);
+		float		ldist = ray_3f_length(&lvec);
+		float		lvec_normal_dot;
+		ray_ray_t	shadow_ray;
+
+		lvec = ray_3f_mult_scalar(&lvec, (1.0f / ldist)); /* normalize lvec */
+#if 1
+		/* skip this light if it's obstructed,
+		 * we must shift the origin slightly towards the light to prevent
+		 * spurious self-obstruction at the ray:object intersection */
+		/* negate the light vector so it's pointed at the light rather than from it */
+		shadow_ray.direction = ray_3f_negate(&lvec);
+		shadow_ray.origin = ray_3f_mult_scalar(&shadow_ray.direction, 0.00001f);
+		shadow_ray.origin = ray_3f_add(&shadow_ray.origin, &intersection);
+
+		if (ray_is_obstructed(scene, &shadow_ray, ldist))
+			continue;
+#endif
+		lvec_normal_dot = ray_3f_dot(&normal, &lvec);
+
+		if (lvec_normal_dot > 0) {
+#if 1
+			float		rvec_lvec_dot = ray_3f_dot(&ray->direction, &lvec);
+			ray_color_t	diffuse;
+			ray_color_t	specular;
+
+			diffuse = ray_3f_mult_scalar(&surface.color, lvec_normal_dot);
+			diffuse = ray_3f_mult_scalar(&diffuse, surface.diffuse);
+			color = ray_3f_add(&color, &diffuse);
+
+			/* FIXME: assumes light is a point for its color, and 20 is a constant "Phong exponent",
+			 * which should really be object/surface-specific
+			 */
+			specular = ray_3f_mult_scalar(&scene->lights[i].light.emitter.point.surface.color, powf(rvec_lvec_dot, 20));
+			specular = ray_3f_mult_scalar(&specular, surface.specular);
+			color = ray_3f_add(&color, &specular);
+#else
+			ray_color_t	diffuse;
+
+			diffuse = ray_3f_mult_scalar(&surface.color, lvec_normal_dot);
+			color = ray_3f_add(&color, &diffuse);
+#endif
+		}
+	}
+
+	/* generate a reflection ray */
+#if 1
+	float		dot = ray_3f_dot(&ray->direction, &normal);
+	ray_ray_t	reflected_ray = { .direction = ray_3f_mult_scalar(&normal, dot * 2.0f) };
+	ray_3f_t	reflection;
+
+	reflected_ray.origin = intersection;
+	reflected_ray.direction = ray_3f_sub(&ray->direction, &reflected_ray.direction);
+
+	reflection = trace_ray(scene, &reflected_ray, depth + 1);
+	reflection = ray_3f_mult_scalar(&reflection, surface.specular);
+	color = ray_3f_add(&color, &reflection);
+#endif
+
+	/* TODO: generate a refraction ray */
+
+	return color;
+}
+
+
+static ray_color_t trace_ray(ray_scene_t *scene, ray_ray_t *ray, unsigned depth)
+{
+	ray_object_t	*nearest_object = NULL;
+	float		nearest_object_distance = INFINITY;
+	ray_color_t	color = { .x = 0.0, .y = 0.0, .z = 0.0 };
+	unsigned	i;
+
+	depth++;
+	if (depth > MAX_RECURSION_DEPTH)
+		return color;
+
+	for (i = 0; i < scene->n_objects; i++) {
+		float	distance;
+
+		/* Does this ray intersect object? */
+		if (ray_object_intersects_ray(&scene->objects[i], ray, &distance)) {
+
+			/* Is it the nearest intersection? */
+			if (!nearest_object ||
+			    distance < nearest_object_distance) {
+				nearest_object = &scene->objects[i];
+				nearest_object_distance = distance;
+			}
+		}
+	}
+
+	if (nearest_object)
+		color = shade_ray(scene, ray, nearest_object, nearest_object_distance, depth);
+
+	depth--;
+
+	return color;
+}
+
+
+void ray_scene_render_fragment(ray_scene_t *scene, ray_camera_t *camera, fb_fragment_t *fragment)
+{
+	ray_camera_frame_t	frame;
+	ray_ray_t		ray;
+	uint32_t		*buf = fragment->buf;
+	unsigned		stride = fragment->stride / 4;
+
+	ray_camera_frame_begin(camera, fragment, &ray, &frame);
+	do {
+		do {
+			*buf = ray_color_to_uint32_rgb(trace_ray(scene, &ray, 0));
+			buf++;
+		} while (ray_camera_frame_x_step(&frame));
+
+		buf += stride;
+	} while (ray_camera_frame_y_step(&frame));
+}
+
+/* we expect fragments[threads->n_threads + 1], or fragments[1] when threads == NULL */
+void ray_scene_render_fragments(ray_scene_t *scene, ray_camera_t *camera, ray_threads_t *threads, fb_fragment_t *fragments)
+{
+	unsigned	n_threads = threads ? threads->n_threads + 1 : 1;
+	unsigned	i;
+
+	for (i = 1; i < n_threads; i++)
+		ray_thread_fragment_submit(&threads->threads[i - 1], scene, camera, &fragments[i]);
+
+	/* always render the zero fragment in-line */
+	ray_scene_render_fragment(scene, camera, &fragments[0]);
+
+	for (i = 1; i < n_threads; i++)
+		ray_thread_wait_idle(&threads->threads[i - 1]);
+}
diff --git a/modules/ray/ray_scene.h b/modules/ray/ray_scene.h
new file mode 100644
index 0000000..e9781c6
--- /dev/null
+++ b/modules/ray/ray_scene.h
@@ -0,0 +1,27 @@
+#ifndef _RAY_SCENE_H
+#define _RAY_SCENE_H
+
+#include "fb.h"
+
+#include "ray_camera.h"
+#include "ray_color.h"
+#include "ray_ray.h"
+#include "ray_threads.h"
+
+typedef union ray_object_t ray_object_t;
+
+typedef struct ray_scene_t {
+	ray_object_t	*objects;
+	unsigned	n_objects;
+
+	ray_object_t	*lights;
+	unsigned	n_lights;
+
+	ray_color_t	ambient_color;
+	float		ambient_brightness;
+} ray_scene_t;
+
+void ray_scene_render_fragments(ray_scene_t *scene, ray_camera_t *camera, ray_threads_t *threads, fb_fragment_t *fragments);
+void ray_scene_render_fragment(ray_scene_t *scene, ray_camera_t *camera, fb_fragment_t *fragment);
+
+#endif
diff --git a/modules/ray/ray_surface.h b/modules/ray/ray_surface.h
new file mode 100644
index 0000000..b3e3c68
--- /dev/null
+++ b/modules/ray/ray_surface.h
@@ -0,0 +1,14 @@
+#ifndef _RAY_MATERIAL_H
+#define _RAY_MATERIAL_H
+
+#include "ray_3f.h"
+#include "ray_color.h"
+
+/* Surface properties we expect every object to be able to introspect */
+typedef struct ray_surface_t {
+	ray_color_t	color;
+	float		specular;
+	float		diffuse;
+} ray_surface_t;
+
+#endif
diff --git a/modules/ray/ray_threads.c b/modules/ray/ray_threads.c
new file mode 100644
index 0000000..2369687
--- /dev/null
+++ b/modules/ray/ray_threads.c
@@ -0,0 +1,111 @@
+#include <pthread.h>
+#include <stdlib.h>
+
+#include "fb.h"
+
+#include "ray_scene.h"
+#include "ray_threads.h"
+
+#define BUSY_WAIT_NUM	1000000000	/* How much to spin before sleeping in pthread_cond_wait() */
+
+/* for now assuming x86 */
+#define cpu_relax() \
+	__asm__ __volatile__ ( "pause\n" : : : "memory")
+
+/* This is a very simple/naive implementation, there's certainly room for improvement.
+ *
+ * Without the BUSY_WAIT_NUM spinning this approach seems to leave a fairly
+ * substantial proportion of CPU idle while waiting for the render thread to
+ * complete on my core 2 duo.
+ *
+ * It's probably just latency in getting the render thread woken when the work
+ * is submitted, and since the fragments are split equally the main thread gets
+ * a head start and has to wait when it finishes first.  The spinning is just
+ * an attempt to avoid going to sleep while the render threads finish, there
+ * still needs to be improvement in how the work is submitted.
+ *
+ * I haven't spent much time on optimizing the raytracer yet.
+ */
+
+static void * ray_thread_func(void *_thread)
+{
+	ray_thread_t	*thread = _thread;
+
+	for (;;) {
+		pthread_mutex_lock(&thread->mutex);
+		while (thread->fragment == NULL)
+			pthread_cond_wait(&thread->cond, &thread->mutex);
+
+		ray_scene_render_fragment(thread->scene, thread->camera, thread->fragment);
+		thread->fragment = NULL;
+		pthread_mutex_unlock(&thread->mutex);
+		pthread_cond_signal(&thread->cond);
+	}
+
+	return NULL;
+}
+
+
+void ray_thread_fragment_submit(ray_thread_t *thread, ray_scene_t *scene, ray_camera_t *camera, fb_fragment_t *fragment)
+{
+	pthread_mutex_lock(&thread->mutex);
+	while (thread->fragment != NULL)	/* XXX: never true due to ray_thread_wait_idle() */
+		pthread_cond_wait(&thread->cond, &thread->mutex);
+
+	thread->fragment = fragment;
+	thread->scene = scene;
+	thread->camera = camera;
+
+	pthread_mutex_unlock(&thread->mutex);
+	pthread_cond_signal(&thread->cond);
+}
+
+
+void ray_thread_wait_idle(ray_thread_t *thread)
+{
+	unsigned	n;
+
+	/* Spin before going to sleep, the other thread should not take substantially longer. */
+	for (n = 0; thread->fragment != NULL && n < BUSY_WAIT_NUM; n++)
+		cpu_relax();
+
+	pthread_mutex_lock(&thread->mutex);
+	while (thread->fragment != NULL)
+		pthread_cond_wait(&thread->cond, &thread->mutex);
+	pthread_mutex_unlock(&thread->mutex);
+}
+
+
+ray_threads_t * ray_threads_create(unsigned num)
+{
+	ray_threads_t	*threads;
+	unsigned	i;
+
+	threads = malloc(sizeof(ray_threads_t) + sizeof(ray_thread_t) * num);
+	if (!threads)
+		return NULL;
+
+	for (i = 0; i < num; i++) {
+		pthread_mutex_init(&threads->threads[i].mutex, NULL);
+		pthread_cond_init(&threads->threads[i].cond, NULL);
+		threads->threads[i].fragment = NULL;
+		pthread_create(&threads->threads[i].thread, NULL, ray_thread_func, &threads->threads[i]);
+	}
+	threads->n_threads = num;
+
+	return threads;
+}
+
+
+void ray_threads_destroy(ray_threads_t *threads)
+{
+	unsigned	i;
+
+	for (i = 0; i < threads->n_threads; i++)
+		pthread_cancel(threads->threads[i].thread);
+
+	for (i = 0; i < threads->n_threads; i++)
+		pthread_join(threads->threads[i].thread, NULL);
+
+	free(threads);	
+}
diff --git a/modules/ray/ray_threads.h b/modules/ray/ray_threads.h
new file mode 100644
index 0000000..b4c601d
--- /dev/null
+++ b/modules/ray/ray_threads.h
@@ -0,0 +1,30 @@
+#ifndef _RAY_THREADS_H
+#define _RAY_THREADS_H
+
+#include <pthread.h>
+
+typedef struct ray_scene_t ray_scene_t;
+typedef struct ray_camera_t ray_camera_t;
+typedef struct fb_fragment_t fb_fragment_t;
+
+typedef struct ray_thread_t {
+	pthread_t	thread;
+	pthread_mutex_t	mutex;
+	pthread_cond_t	cond;
+	ray_scene_t	*scene;
+	ray_camera_t	*camera;
+	fb_fragment_t	*fragment;
+} ray_thread_t;
+
+typedef struct ray_threads_t {
+	unsigned	n_threads;
+	ray_thread_t	threads[];
+} ray_threads_t;
+
+
+ray_threads_t * ray_threads_create(unsigned num);
+void ray_threads_destroy(ray_threads_t *threads);
+
+void ray_thread_fragment_submit(ray_thread_t *thread, ray_scene_t *scene, ray_camera_t *camera, fb_fragment_t *fragment);
+void ray_thread_wait_idle(ray_thread_t *thread);
+#endif