Modularize the libdimension codebase.

8 files changed, 1582 insertions, 0 deletions

diff --git a/libdimension/model/objects/cone.c b/libdimension/model/objects/cone.c
new file mode 100644
index 0000000..26e59ca
--- /dev/null
+++ b/libdimension/model/objects/cone.c
@@ -0,0 +1,207 @@
+/*************************************************************************
+ * Copyright (C) 2009-2014 Tavian Barnes <tavianator@tavianator.com>     *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Cones/cylinders.
+ */
+
+#include "internal/polynomial.h"
+#include "dimension/model.h"
+#include <math.h>
+
+/// Cone type.
+typedef struct dmnsn_cone {
+  dmnsn_object object;
+  double r1, r2;
+} dmnsn_cone;
+
+/// Intersection callback for a cone.
+static bool
+dmnsn_cone_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
+                           dmnsn_intersection *intersection)
+{
+  const dmnsn_cone *cone = (const dmnsn_cone *)object;
+  double r1 = cone->r1, r2 = cone->r2;
+
+  // Solve (x0 + nx*t)^2 + (z0 + nz*t)^2 == (((r2 - r1)*(y0 + ny*t) + r1 + r2)/2)^2
+  double poly[3], x[2];
+  poly[2] = l.n.x*l.n.x + l.n.z*l.n.z - l.n.y*l.n.y*(r2 - r1)*(r2 - r1)/4.0;
+  poly[1] = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z)
+            - l.n.y*(r2 - r1)*(l.x0.y*(r2 - r1) + r2 + r1)/2.0;
+  poly[0] = l.x0.x*l.x0.x + l.x0.z*l.x0.z
+            - (l.x0.y*(r2 - r1) + r2 + r1)*(l.x0.y*(r2 - r1) + r2 + r1)/4.0;
+
+  size_t n = dmnsn_polynomial_solve(poly, 2, x);
+
+  if (n > 0) {
+    double t = x[0];
+    dmnsn_vector p;
+    if (n == 2) {
+      t = dmnsn_min(t, x[1]);
+      p = dmnsn_ray_point(l, t);
+
+      if (p.y <= -1.0 || p.y >= 1.0) {
+        t = dmnsn_max(x[0], x[1]);
+        p = dmnsn_ray_point(l, t);
+      }
+    } else {
+      p = dmnsn_ray_point(l, t);
+    }
+
+    if (t >= 0.0 && p.y >= -1.0 && p.y <= 1.0) {
+      double r = ((r2 - r1)*p.y + r1 + r2)/2.0;
+      dmnsn_vector norm = dmnsn_new_vector(p.x, -r*(r2 - r1)/2.0, p.z);
+      intersection->t      = t;
+      intersection->normal = norm;
+      return true;
+    }
+  }
+
+  return false;
+}
+
+/// Inside callback for a cone.
+static bool
+dmnsn_cone_inside_fn(const dmnsn_object *object, dmnsn_vector point)
+{
+  const dmnsn_cone *cone = (const dmnsn_cone *)object;
+  double r1 = cone->r1, r2 = cone->r2;
+  double r = (point.y*(r2 - r1) + r1 + r2)/2.0;
+  return point.x*point.x + point.z*point.z < r*r
+         && point.y > -1.0 && point.y < 1.0;
+}
+
+/// Cone bounding callback.
+static dmnsn_aabb
+dmnsn_cone_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  const dmnsn_cone *cone = (const dmnsn_cone *)object;
+
+  double rmax = dmnsn_max(cone->r1, cone->r2);
+  dmnsn_aabb box = dmnsn_symmetric_aabb(dmnsn_new_vector(rmax, 1.0, rmax));
+  return dmnsn_transform_aabb(trans, box);
+}
+
+/// Cone vtable.
+static const dmnsn_object_vtable dmnsn_cone_vtable = {
+  .intersection_fn = dmnsn_cone_intersection_fn,
+  .inside_fn = dmnsn_cone_inside_fn,
+  .bounding_fn = dmnsn_cone_bounding_fn,
+};
+
+/// Cone cap type.
+typedef struct dmnsn_cone_cap {
+  dmnsn_object object;
+  double r;
+} dmnsn_cone_cap;
+
+/// Cone cap intersection function.
+static bool
+dmnsn_cone_cap_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
+                               dmnsn_intersection *intersection)
+{
+  if (l.n.y != 0.0) {
+    const dmnsn_cone_cap *cap = (const dmnsn_cone_cap *)object;
+    double r = cap->r;
+    double t = -l.x0.y/l.n.y;
+    dmnsn_vector p = dmnsn_ray_point(l, t);
+    if (t >= 0.0 && p.x*p.x + p.z*p.z <= r*r) {
+      intersection->t      = t;
+      intersection->normal = dmnsn_new_vector(0.0, -1.0, 0.0);
+      return true;
+    }
+  }
+
+  return false;
+}
+
+/// Inside callback for a cone cap.
+static bool
+dmnsn_cone_cap_inside_fn(const dmnsn_object *object, dmnsn_vector point)
+{
+  return false;
+}
+
+/// Cone cap bounding callback.
+static dmnsn_aabb
+dmnsn_cone_cap_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  const dmnsn_cone_cap *cap = (const dmnsn_cone_cap *)object;
+  dmnsn_aabb box = dmnsn_symmetric_aabb(dmnsn_new_vector(cap->r, 0.0, cap->r));
+  return dmnsn_transform_aabb(trans, box);
+}
+
+/// Cone cap vtable.
+static const dmnsn_object_vtable dmnsn_cone_cap_vtable = {
+  .intersection_fn = dmnsn_cone_cap_intersection_fn,
+  .inside_fn = dmnsn_cone_cap_inside_fn,
+  .bounding_fn = dmnsn_cone_cap_bounding_fn,
+};
+
+/// Allocate a new cone cap.
+dmnsn_object *
+dmnsn_new_cone_cap(dmnsn_pool *pool, double r)
+{
+  dmnsn_cone_cap *cap = DMNSN_PALLOC(pool, dmnsn_cone_cap);
+  cap->r = r;
+
+  dmnsn_object *object = &cap->object;
+  dmnsn_init_object(object);
+  object->vtable = &dmnsn_cone_cap_vtable;
+  return object;
+}
+
+// Allocate a new cone object
+dmnsn_object *
+dmnsn_new_cone(dmnsn_pool *pool, double r1, double r2, bool open)
+{
+  dmnsn_cone *cone = DMNSN_PALLOC(pool, dmnsn_cone);
+  cone->r1 = r1;
+  cone->r2 = r2;
+
+  dmnsn_object *object = &cone->object;
+  dmnsn_init_object(object);
+  object->vtable = &dmnsn_cone_vtable;
+
+  if (open) {
+    return object;
+  }
+
+  // Implement closed cones as a union with the caps
+  dmnsn_object *cap1 = dmnsn_new_cone_cap(pool, r1);
+  dmnsn_object *cap2 = dmnsn_new_cone_cap(pool, r2);
+  cap1->intrinsic_trans = dmnsn_translation_matrix(
+    dmnsn_new_vector(0.0, -1.0, 0.0)
+  );
+  cap2->intrinsic_trans = dmnsn_translation_matrix(
+    dmnsn_new_vector(0.0, +1.0, 0.0)
+  );
+  // Flip the normal around for the top cap
+  cap2->intrinsic_trans.n[1][1] = -1.0;
+
+  dmnsn_array *withcaps = DMNSN_PALLOC_ARRAY(pool, dmnsn_object *);
+  dmnsn_array_push(withcaps, &cone);
+  dmnsn_array_push(withcaps, &cap1);
+  dmnsn_array_push(withcaps, &cap2);
+  dmnsn_object *cone_cap_union = dmnsn_new_csg_union(pool, withcaps);
+
+  return cone_cap_union;
+}
diff --git a/libdimension/model/objects/csg.c b/libdimension/model/objects/csg.c
new file mode 100644
index 0000000..15008c0
--- /dev/null
+++ b/libdimension/model/objects/csg.c
@@ -0,0 +1,334 @@
+/*************************************************************************
+ * Copyright (C) 2010-2014 Tavian Barnes <tavianator@tavianator.com>     *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Constructive solid geometry.
+ */
+
+#include "internal.h"
+#include "internal/bvh.h"
+#include "dimension/model.h"
+#include <stdlib.h>
+
+////////////
+// Unions //
+////////////
+
+typedef struct {
+  dmnsn_object object;
+  dmnsn_bvh *bvh;
+} dmnsn_csg_union;
+
+/// CSG union intersection callback.
+static bool
+dmnsn_csg_union_intersection_fn(const dmnsn_object *object,
+                                dmnsn_ray ray,
+                                dmnsn_intersection *intersection)
+{
+  const dmnsn_csg_union *csg = (const dmnsn_csg_union *)object;
+  return dmnsn_bvh_intersection(csg->bvh, ray, intersection, true);
+}
+
+/// CSG union inside callback.
+static bool
+dmnsn_csg_union_inside_fn(const dmnsn_object *object, dmnsn_vector point)
+{
+  const dmnsn_csg_union *csg = (const dmnsn_csg_union *)object;
+  return dmnsn_bvh_inside(csg->bvh, point);
+}
+
+/// CSG union precomputation callback.
+static void
+dmnsn_csg_union_precompute_fn(dmnsn_object *object)
+{
+  dmnsn_csg_union *csg = (dmnsn_csg_union *)object;
+  csg->object.trans_inv = dmnsn_identity_matrix();
+
+  dmnsn_bvh *bvh = dmnsn_new_bvh(csg->object.children, DMNSN_BVH_PRTREE);
+  csg->bvh = bvh;
+  csg->object.aabb = dmnsn_bvh_aabb(bvh);
+}
+
+/// CSG union vtable.
+static const dmnsn_object_vtable dmnsn_csg_union_vtable = {
+  .intersection_fn = dmnsn_csg_union_intersection_fn,
+  .inside_fn = dmnsn_csg_union_inside_fn,
+  .precompute_fn = dmnsn_csg_union_precompute_fn,
+};
+
+/// CSG union destruction callback.
+static void
+dmnsn_csg_union_cleanup(void *ptr)
+{
+  dmnsn_csg_union *csg = ptr;
+  dmnsn_delete_bvh(csg->bvh);
+}
+
+// Bulk-load a union
+dmnsn_object *
+dmnsn_new_csg_union(dmnsn_pool *pool, dmnsn_array *objects)
+{
+  dmnsn_csg_union *csg = DMNSN_PALLOC_TIDY(pool, dmnsn_csg_union, dmnsn_csg_union_cleanup);
+  csg->bvh = NULL;
+
+  dmnsn_object *object = &csg->object;
+  dmnsn_init_object(object);
+
+  object->vtable = &dmnsn_csg_union_vtable;
+  object->children = objects;
+  object->split_children = true;
+
+  return object;
+}
+
+/**
+ * Generic CSG intersection callback.
+ * @param[in]  csg           The CSG object.
+ * @param[in]  ray           The intersection ray.
+ * @param[out] intersection  The intersection data.
+ * @param[in]  inside1       Whether the first object is allowed inside the
+ *                           second object.
+ * @param[in]  inside2       Whether the second object is allowed inside the
+ *                           first object.
+ * @return Whether \p ray intersected \p csg.
+ */
+static bool
+dmnsn_csg_intersection_fn(const dmnsn_object *csg, dmnsn_ray ray,
+                          dmnsn_intersection *intersection,
+                          bool inside1, bool inside2)
+{
+  const dmnsn_object *A = *(dmnsn_object **)dmnsn_array_first(csg->children);
+  const dmnsn_object *B = *(dmnsn_object **)dmnsn_array_last(csg->children);
+
+  dmnsn_intersection i1, i2;
+  bool is_i1 = dmnsn_object_intersection(A, ray, &i1);
+  bool is_i2 = dmnsn_object_intersection(B, ray, &i2);
+
+  double oldt = 0.0;
+  while (is_i1) {
+    i1.ray = ray;
+    i1.t += oldt;
+    oldt = i1.t + dmnsn_epsilon;
+
+    dmnsn_vector point = dmnsn_ray_point(i1.ray, i1.t);
+    if (inside2 ^ dmnsn_object_inside(B, point)) {
+      dmnsn_ray newray = ray;
+      newray.x0 = dmnsn_ray_point(ray, i1.t);
+      newray    = dmnsn_ray_add_epsilon(newray);
+      is_i1 = dmnsn_object_intersection(A, newray, &i1);
+    } else {
+      break;
+    }
+  }
+
+  oldt = 0.0;
+  while (is_i2) {
+    i2.ray = ray;
+    i2.t += oldt;
+    oldt = i2.t + dmnsn_epsilon;
+
+    dmnsn_vector point = dmnsn_ray_point(i2.ray, i2.t);
+    if (inside1 ^ dmnsn_object_inside(A, point)) {
+      dmnsn_ray newray = ray;
+      newray.x0 = dmnsn_ray_point(ray, i2.t);
+      newray    = dmnsn_ray_add_epsilon(newray);
+      is_i2 = dmnsn_object_intersection(B, newray, &i2);
+    } else {
+      break;
+    }
+  }
+
+  if (is_i1 && is_i2) {
+    if (i1.t < i2.t) {
+      *intersection = i1;
+    } else {
+      *intersection = i2;
+    }
+  } else if (is_i1) {
+    *intersection = i1;
+  } else if (is_i2) {
+    *intersection = i2;
+  } else {
+    return false;
+  }
+
+  return true;
+}
+
+///////////////////
+// Intersections //
+///////////////////
+
+/// CSG intersection intersection callback.
+static bool
+dmnsn_csg_intersection_intersection_fn(const dmnsn_object *csg,
+                                       dmnsn_ray ray,
+                                       dmnsn_intersection *intersection)
+{
+  return dmnsn_csg_intersection_fn(csg, ray, intersection, true, true);
+}
+
+/// CSG intersection inside callback.
+static bool
+dmnsn_csg_intersection_inside_fn(const dmnsn_object *csg, dmnsn_vector point)
+{
+  const dmnsn_object *A = *(dmnsn_object **)dmnsn_array_first(csg->children);
+  const dmnsn_object *B = *(dmnsn_object **)dmnsn_array_last(csg->children);
+  return dmnsn_object_inside(A, point) && dmnsn_object_inside(B, point);
+}
+
+/// CSG intersection precomputation callback.
+static void
+dmnsn_csg_intersection_precompute_fn(dmnsn_object *csg)
+{
+  dmnsn_object *A = *(dmnsn_object **)dmnsn_array_first(csg->children);
+  dmnsn_object *B = *(dmnsn_object **)dmnsn_array_last(csg->children);
+
+  csg->trans_inv = dmnsn_identity_matrix();
+  csg->aabb.min = dmnsn_vector_max(A->aabb.min, B->aabb.min);
+  csg->aabb.max = dmnsn_vector_min(A->aabb.max, B->aabb.max);
+}
+
+/// CSG intersection vtable.
+static const dmnsn_object_vtable dmnsn_csg_intersection_vtable = {
+  .intersection_fn = dmnsn_csg_intersection_intersection_fn,
+  .inside_fn = dmnsn_csg_intersection_inside_fn,
+  .precompute_fn = dmnsn_csg_intersection_precompute_fn,
+};
+
+dmnsn_object *
+dmnsn_new_csg_intersection(dmnsn_pool *pool, dmnsn_object *A, dmnsn_object *B)
+{
+  dmnsn_object *csg = dmnsn_new_object(pool);
+  csg->vtable = &dmnsn_csg_intersection_vtable;
+
+  csg->children = DMNSN_PALLOC_ARRAY(pool, dmnsn_object *);
+  dmnsn_array_push(csg->children, &A);
+  dmnsn_array_push(csg->children, &B);
+
+  return csg;
+}
+
+/////////////////
+// Differences //
+/////////////////
+
+/// CSG difference intersection callback.
+static bool
+dmnsn_csg_difference_intersection_fn(const dmnsn_object *csg,
+                                     dmnsn_ray ray,
+                                     dmnsn_intersection *intersection)
+{
+  return dmnsn_csg_intersection_fn(csg, ray, intersection, true, false);
+}
+
+/// CSG difference inside callback.
+static bool
+dmnsn_csg_difference_inside_fn(const dmnsn_object *csg, dmnsn_vector point)
+{
+  const dmnsn_object *A = *(dmnsn_object **)dmnsn_array_first(csg->children);
+  const dmnsn_object *B = *(dmnsn_object **)dmnsn_array_last(csg->children);
+  return dmnsn_object_inside(A, point)  && !dmnsn_object_inside(B, point);
+}
+
+/// CSG difference precomputation callback.
+static void
+dmnsn_csg_difference_precompute_fn(dmnsn_object *csg)
+{
+  dmnsn_object *A = *(dmnsn_object **)dmnsn_array_first(csg->children);
+
+  csg->trans_inv = dmnsn_identity_matrix();
+  csg->aabb = A->aabb;
+}
+
+/// CSG difference vtable.
+static const dmnsn_object_vtable dmnsn_csg_difference_vtable = {
+  .intersection_fn = dmnsn_csg_difference_intersection_fn,
+  .inside_fn = dmnsn_csg_difference_inside_fn,
+  .precompute_fn = dmnsn_csg_difference_precompute_fn,
+};
+
+dmnsn_object *
+dmnsn_new_csg_difference(dmnsn_pool *pool, dmnsn_object *A, dmnsn_object *B)
+{
+  dmnsn_object *csg = dmnsn_new_object(pool);
+  csg->vtable = &dmnsn_csg_difference_vtable;
+
+  csg->children = DMNSN_PALLOC_ARRAY(pool, dmnsn_object *);
+  dmnsn_array_push(csg->children, &A);
+  dmnsn_array_push(csg->children, &B);
+
+  return csg;
+}
+
+////////////
+// Merges //
+////////////
+
+/// CSG merge intersection callback.
+static bool
+dmnsn_csg_merge_intersection_fn(const dmnsn_object *csg,
+                                dmnsn_ray ray,
+                                dmnsn_intersection *intersection)
+{
+  return dmnsn_csg_intersection_fn(csg, ray, intersection, false, false);
+}
+
+/// CSG merge inside callback.
+static bool
+dmnsn_csg_merge_inside_fn(const dmnsn_object *csg, dmnsn_vector point)
+{
+  const dmnsn_object *A = *(dmnsn_object **)dmnsn_array_first(csg->children);
+  const dmnsn_object *B = *(dmnsn_object **)dmnsn_array_last(csg->children);
+  return dmnsn_object_inside(A, point) || dmnsn_object_inside(B, point);
+}
+
+/// CSG merge precomputation callback.
+static void
+dmnsn_csg_merge_precompute_fn(dmnsn_object *csg)
+{
+  dmnsn_object *A = *(dmnsn_object **)dmnsn_array_first(csg->children);
+  dmnsn_object *B = *(dmnsn_object **)dmnsn_array_last(csg->children);
+
+  csg->trans_inv = dmnsn_identity_matrix();
+  csg->aabb.min = dmnsn_vector_min(A->aabb.min, B->aabb.min);
+  csg->aabb.max = dmnsn_vector_max(A->aabb.max, B->aabb.max);
+}
+
+/// CSG merge vtable.
+static const dmnsn_object_vtable dmnsn_csg_merge_vtable = {
+  .intersection_fn = dmnsn_csg_merge_intersection_fn,
+  .inside_fn = dmnsn_csg_merge_inside_fn,
+  .precompute_fn = dmnsn_csg_merge_precompute_fn,
+};
+
+dmnsn_object *
+dmnsn_new_csg_merge(dmnsn_pool *pool, dmnsn_object *A, dmnsn_object *B)
+{
+  dmnsn_object *csg = dmnsn_new_object(pool);
+  csg->vtable = &dmnsn_csg_merge_vtable;
+
+  csg->children = DMNSN_PALLOC_ARRAY(pool, dmnsn_object *);
+  dmnsn_array_push(csg->children, &A);
+  dmnsn_array_push(csg->children, &B);
+
+  return csg;
+}
diff --git a/libdimension/model/objects/cube.c b/libdimension/model/objects/cube.c
new file mode 100644
index 0000000..7d6fe0f
--- /dev/null
+++ b/libdimension/model/objects/cube.c
@@ -0,0 +1,154 @@
+/*************************************************************************
+ * Copyright (C) 2009-2014 Tavian Barnes <tavianator@tavianator.com>     *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Cubes.
+ */
+
+#include "dimension/model.h"
+#include <math.h>
+
+/// Intersection callback for a cube.
+static bool
+dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray,
+                           dmnsn_intersection *intersection)
+{
+  // Clip the given ray against the X, Y, and Z slabs
+
+  dmnsn_vector nmin, nmax;
+  double tmin, tmax;
+
+  double tx1 = (-1.0 - ray.x0.x)/ray.n.x;
+  double tx2 = (+1.0 - ray.x0.x)/ray.n.x;
+
+  if (tx1 < tx2) {
+    tmin = tx1;
+    tmax = tx2;
+    nmin = dmnsn_new_vector(-1.0, 0.0, 0.0);
+    nmax = dmnsn_new_vector(+1.0, 0.0, 0.0);
+  } else {
+    tmin = tx2;
+    tmax = tx1;
+    nmin = dmnsn_new_vector(+1.0, 0.0, 0.0);
+    nmax = dmnsn_new_vector(-1.0, 0.0, 0.0);
+  }
+
+  if (tmin > tmax)
+    return false;
+
+  double ty1 = (-1.0 - ray.x0.y)/ray.n.y;
+  double ty2 = (+1.0 - ray.x0.y)/ray.n.y;
+
+  if (ty1 < ty2) {
+    if (ty1 > tmin) {
+      tmin = ty1;
+      nmin = dmnsn_new_vector(0.0, -1.0, 0.0);
+    }
+    if (ty2 < tmax) {
+      tmax = ty2;
+      nmax = dmnsn_new_vector(0.0, +1.0, 0.0);
+    }
+  } else {
+    if (ty2 > tmin) {
+      tmin = ty2;
+      nmin = dmnsn_new_vector(0.0, +1.0, 0.0);
+    }
+    if (ty1 < tmax) {
+      tmax = ty1;
+      nmax = dmnsn_new_vector(0.0, -1.0, 0.0);
+    }
+  }
+
+  if (tmin > tmax)
+    return false;
+
+  double tz1 = (-1.0 - ray.x0.z)/ray.n.z;
+  double tz2 = (+1.0 - ray.x0.z)/ray.n.z;
+
+  if (tz1 < tz2) {
+    if (tz1 > tmin) {
+      tmin = tz1;
+      nmin = dmnsn_new_vector(0.0, 0.0, -1.0);
+    }
+    if (tz2 < tmax) {
+      tmax = tz2;
+      nmax = dmnsn_new_vector(0.0, 0.0, +1.0);
+    }
+  } else {
+    if (tz2 > tmin) {
+      tmin = tz2;
+      nmin = dmnsn_new_vector(0.0, 0.0, +1.0);
+    }
+    if (tz1 < tmax) {
+      tmax = tz1;
+      nmax = dmnsn_new_vector(0.0, 0.0, -1.0);
+    }
+  }
+
+  if (tmin > tmax)
+    return false;
+
+  if (tmin < 0.0) {
+    tmin = tmax;
+    nmin = nmax;
+  }
+
+  if (tmin >= 0.0) {
+    intersection->t      = tmin;
+    intersection->normal = nmin;
+    return true;
+  } else {
+    return false;
+  }
+}
+
+/// Inside callback for a cube.
+static bool
+dmnsn_cube_inside_fn(const dmnsn_object *cube, dmnsn_vector point)
+{
+  return point.x > -1.0 && point.x < 1.0
+      && point.y > -1.0 && point.y < 1.0
+      && point.z > -1.0 && point.z < 1.0;
+}
+
+/// Boundary callback for a cube.
+static dmnsn_aabb
+dmnsn_cube_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  dmnsn_aabb box = dmnsn_symmetric_aabb(dmnsn_new_vector(1.0, 1.0, 1.0));
+  return dmnsn_transform_aabb(trans, box);
+}
+
+/// Cube vtable.
+static const dmnsn_object_vtable dmnsn_cube_vtable = {
+  .intersection_fn = dmnsn_cube_intersection_fn,
+  .inside_fn = dmnsn_cube_inside_fn,
+  .bounding_fn = dmnsn_cube_bounding_fn,
+};
+
+// Allocate a new cube object
+dmnsn_object *
+dmnsn_new_cube(dmnsn_pool *pool)
+{
+  dmnsn_object *cube = dmnsn_new_object(pool);
+  cube->vtable = &dmnsn_cube_vtable;
+  return cube;
+}
diff --git a/libdimension/model/objects/plane.c b/libdimension/model/objects/plane.c
new file mode 100644
index 0000000..b34d8aa
--- /dev/null
+++ b/libdimension/model/objects/plane.c
@@ -0,0 +1,88 @@
+/*************************************************************************
+ * Copyright (C) 2010-2014 Tavian Barnes <tavianator@tavianator.com>     *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Planes.
+ */
+
+#include "dimension/model.h"
+#include <math.h>
+#include <stdlib.h>
+
+/// Plane type.
+typedef struct {
+  dmnsn_object object;
+  dmnsn_vector normal;
+} dmnsn_plane;
+
+/// Returns the closest intersection of `ray' with `plane'.
+static bool
+dmnsn_plane_intersection_fn(const dmnsn_object *object, dmnsn_ray ray,
+                            dmnsn_intersection *intersection)
+{
+  const dmnsn_plane *plane = (const dmnsn_plane *)object;
+  dmnsn_vector normal = plane->normal;
+
+  double den = dmnsn_vector_dot(ray.n, normal);
+  if (den != 0.0) {
+    double t = -dmnsn_vector_dot(ray.x0, normal)/den;
+    if (t >= 0.0) {
+      intersection->t      = t;
+      intersection->normal = normal;
+      return true;
+    }
+  }
+  return false;
+}
+
+/// Return whether a point is inside a plane.
+static bool
+dmnsn_plane_inside_fn(const dmnsn_object *object, dmnsn_vector point)
+{
+  const dmnsn_plane *plane = (const dmnsn_plane *)object;
+  return dmnsn_vector_dot(point, plane->normal) < 0.0;
+}
+
+/// Plane bounding callback.
+static dmnsn_aabb
+dmnsn_plane_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  return dmnsn_infinite_aabb();
+}
+
+/// Plane vtable.
+static const dmnsn_object_vtable dmnsn_plane_vtable = {
+  .intersection_fn = dmnsn_plane_intersection_fn,
+  .inside_fn = dmnsn_plane_inside_fn,
+  .bounding_fn = dmnsn_plane_bounding_fn,
+};
+
+dmnsn_object *
+dmnsn_new_plane(dmnsn_pool *pool, dmnsn_vector normal)
+{
+  dmnsn_plane *plane = DMNSN_PALLOC(pool, dmnsn_plane);
+  plane->normal = normal;
+
+  dmnsn_object *object = &plane->object;
+  dmnsn_init_object(object);
+  object->vtable = &dmnsn_plane_vtable;
+  return object;
+}
diff --git a/libdimension/model/objects/sphere.c b/libdimension/model/objects/sphere.c
new file mode 100644
index 0000000..e1ca784
--- /dev/null
+++ b/libdimension/model/objects/sphere.c
@@ -0,0 +1,116 @@
+/*************************************************************************
+ * Copyright (C) 2009-2014 Tavian Barnes <tavianator@tavianator.com>     *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Spheres.
+ */
+
+#include "internal.h"
+#include "internal/polynomial.h"
+#include "dimension/model.h"
+
+/// Sphere intersection callback.
+static bool
+dmnsn_sphere_intersection_fn(const dmnsn_object *sphere, dmnsn_ray l,
+                             dmnsn_intersection *intersection)
+{
+  // Solve (x0 + nx*t)^2 + (y0 + ny*t)^2 + (z0 + nz*t)^2 == 1
+  double poly[3], x[2];
+  poly[2] = dmnsn_vector_dot(l.n, l.n);
+  poly[1] = 2.0*dmnsn_vector_dot(l.n, l.x0);
+  poly[0] = dmnsn_vector_dot(l.x0, l.x0) - 1.0;
+
+  size_t n = dmnsn_polynomial_solve(poly, 2, x);
+  if (n == 0) {
+    return false;
+  }
+
+  double t = x[0];
+  // Optimize for the case where we're outside the sphere
+  if (dmnsn_likely(n == 2)) {
+    t = dmnsn_min(t, x[1]);
+  }
+
+  intersection->t = t;
+  intersection->normal = dmnsn_ray_point(l, t);
+  return true;
+}
+
+/// Sphere inside callback.
+static bool
+dmnsn_sphere_inside_fn(const dmnsn_object *sphere, dmnsn_vector point)
+{
+  return point.x*point.x + point.y*point.y + point.z*point.z < 1.0;
+}
+
+/// Helper for sphere bounding box calculation.
+static inline double
+dmnsn_implicit_dot(double row[4])
+{
+  double ret = 0.0;
+  for (int i = 0; i < 3; ++i) {
+    ret += row[i]*row[i];
+  }
+  return ret;
+}
+
+/// Sphere bounding callback.
+static dmnsn_aabb
+dmnsn_sphere_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  // Get a tight bound using the quadric representation of a sphere.  For
+  // details, see
+  // http://tavianator.com/2014/06/exact-bounding-boxes-for-spheres-ellipsoids
+
+  dmnsn_aabb box;
+
+  double cx = trans.n[0][3];
+  double dx = sqrt(dmnsn_implicit_dot(trans.n[0]));
+  box.min.x = cx - dx;
+  box.max.x = cx + dx;
+
+  double cy = trans.n[1][3];
+  double dy = sqrt(dmnsn_implicit_dot(trans.n[1]));
+  box.min.y = cy - dy;
+  box.max.y = cy + dy;
+
+  double cz = trans.n[2][3];
+  double dz = sqrt(dmnsn_implicit_dot(trans.n[2]));
+  box.min.z = cz - dz;
+  box.max.z = cz + dz;
+
+  return box;
+}
+
+/// Sphere vtable.
+static const dmnsn_object_vtable dmnsn_sphere_vtable = {
+  .intersection_fn = dmnsn_sphere_intersection_fn,
+  .inside_fn = dmnsn_sphere_inside_fn,
+  .bounding_fn = dmnsn_sphere_bounding_fn,
+};
+
+dmnsn_object *
+dmnsn_new_sphere(dmnsn_pool *pool)
+{
+  dmnsn_object *sphere = dmnsn_new_object(pool);
+  sphere->vtable = &dmnsn_sphere_vtable;
+  return sphere;
+}
diff --git a/libdimension/model/objects/torus.c b/libdimension/model/objects/torus.c
new file mode 100644
index 0000000..b4baebd
--- /dev/null
+++ b/libdimension/model/objects/torus.c
@@ -0,0 +1,173 @@
+/*************************************************************************
+ * Copyright (C) 2010-2014 Tavian Barnes <tavianator@tavianator.com>     *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Torii.  A special case of a quartic.
+ */
+
+#include "internal/polynomial.h"
+#include "dimension/model.h"
+
+/// Torus type.
+typedef struct {
+  dmnsn_object object;
+  double major, minor;
+} dmnsn_torus;
+
+/// Bound the torus in a cylindrical shell.
+static inline bool
+dmnsn_torus_bound_intersection(const dmnsn_torus *torus, dmnsn_ray l)
+{
+  double R = torus->major, r = torus->minor;
+  double rmax = R + r, rmin = R - r;
+  double rmax2 = rmax*rmax, rmin2 = rmin*rmin;
+
+  // Try the caps first
+  double tlower = (-r - l.x0.y)/l.n.y;
+  double tupper = (+r - l.x0.y)/l.n.y;
+  dmnsn_vector lower = dmnsn_ray_point(l, tlower);
+  dmnsn_vector upper = dmnsn_ray_point(l, tupper);
+  double ldist2 = lower.x*lower.x + lower.z*lower.z;
+  double udist2 = upper.x*upper.x + upper.z*upper.z;
+  if ((ldist2 < rmin2 || ldist2 > rmax2) && (udist2 < rmin2 || udist2 > rmax2)) {
+    // No valid intersection with the caps, try the cylinder walls
+    double dist2 = l.x0.x*l.x0.x + l.x0.z*l.x0.z;
+    double bigcyl[3], smallcyl[3];
+    bigcyl[2]   = smallcyl[2] = l.n.x*l.n.x + l.n.z*l.n.z;
+    bigcyl[1]   = smallcyl[1] = 2.0*(l.n.x*l.x0.x + l.n.z*l.x0.z);
+    bigcyl[0]   = dist2 - rmax2;
+    smallcyl[0] = dist2 - rmin2;
+
+    double x[4];
+    size_t n = dmnsn_polynomial_solve(bigcyl, 2, x);
+    n += dmnsn_polynomial_solve(smallcyl, 2, x + n);
+
+    size_t i;
+    for (i = 0; i < n; ++i) {
+      dmnsn_vector p = dmnsn_ray_point(l, x[i]);
+      if (p.y >= -r && p.y <= r)
+        break;
+    }
+
+    if (i == n) {
+      // No valid intersection found
+      return false;
+    }
+  }
+
+  return true;
+}
+
+/// Torus intersection callback.
+static bool
+dmnsn_torus_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
+                            dmnsn_intersection *intersection)
+{
+  const dmnsn_torus *torus = (const dmnsn_torus *)object;
+  double R = torus->major, r = torus->minor;
+  double RR = R*R, rr = r*r;
+
+  if (!dmnsn_torus_bound_intersection(torus, l)) {
+    return false;
+  }
+
+  // This bit of algebra here is correct
+  dmnsn_vector x0mod = dmnsn_new_vector(l.x0.x, -l.x0.y, l.x0.z);
+  dmnsn_vector nmod  = dmnsn_new_vector(l.n.x,  -l.n.y,  l.n.z);
+  double nn      = dmnsn_vector_dot(l.n, l.n);
+  double nx0     = dmnsn_vector_dot(l.n, l.x0);
+  double x0x0    = dmnsn_vector_dot(l.x0, l.x0);
+  double x0x0mod = dmnsn_vector_dot(l.x0, x0mod);
+  double nx0mod  = dmnsn_vector_dot(l.n, x0mod);
+  double nnmod   = dmnsn_vector_dot(l.n, nmod);
+
+  double poly[5];
+  poly[4] = nn*nn;
+  poly[3] = 4*nn*nx0;
+  poly[2] = 2.0*(nn*(x0x0 - rr) + 2.0*nx0*nx0 - RR*nnmod);
+  poly[1] = 4.0*(nx0*(x0x0 - rr) - RR*nx0mod);
+  poly[0] = x0x0*x0x0 + RR*(RR - 2.0*x0x0mod) - rr*(2.0*(RR + x0x0) - rr);
+
+  double x[4];
+  size_t n = dmnsn_polynomial_solve(poly, 4, x);
+  if (n == 0)
+    return false;
+
+  double t = x[0];
+  for (size_t i = 1; i < n; ++i) {
+    t = dmnsn_min(t, x[i]);
+  }
+
+  if (t < 0.0) {
+    return false;
+  }
+
+  dmnsn_vector p = dmnsn_ray_point(l, t);
+  dmnsn_vector center = dmnsn_vector_mul(
+    R,
+    dmnsn_vector_normalized(dmnsn_new_vector(p.x, 0.0, p.z))
+  );
+  dmnsn_vector normal = dmnsn_vector_sub(p, center);
+
+  intersection->t      = t;
+  intersection->normal = normal;
+  return true;
+}
+
+/// Torus inside callback.
+static bool
+dmnsn_torus_inside_fn(const dmnsn_object *object, dmnsn_vector point)
+{
+  const dmnsn_torus *torus = (const dmnsn_torus *)object;
+  double dmajor = torus->major - sqrt(point.x*point.x + point.z*point.z);
+  return dmajor*dmajor + point.y*point.y < torus->minor*torus->minor;
+}
+
+/// Torus bounding callback.
+static dmnsn_aabb
+dmnsn_torus_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  const dmnsn_torus *torus = (const dmnsn_torus *)object;
+
+  double extent = torus->major + torus->minor;
+  dmnsn_aabb box = dmnsn_symmetric_aabb(dmnsn_new_vector(extent, torus->minor, extent));
+  return dmnsn_transform_aabb(trans, box);
+}
+
+/// Torus vtable.
+static const dmnsn_object_vtable dmnsn_torus_vtable = {
+  .intersection_fn = dmnsn_torus_intersection_fn,
+  .inside_fn = dmnsn_torus_inside_fn,
+  .bounding_fn = dmnsn_torus_bounding_fn,
+};
+
+dmnsn_object *
+dmnsn_new_torus(dmnsn_pool *pool, double major, double minor)
+{
+  dmnsn_torus *torus = DMNSN_PALLOC(pool, dmnsn_torus);
+  torus->major = major;
+  torus->minor = minor;
+
+  dmnsn_object *object = &torus->object;
+  dmnsn_init_object(object);
+  object->vtable = &dmnsn_torus_vtable;
+  return object;
+}
diff --git a/libdimension/model/objects/triangle.c b/libdimension/model/objects/triangle.c
new file mode 100644
index 0000000..5af3301
--- /dev/null
+++ b/libdimension/model/objects/triangle.c
@@ -0,0 +1,163 @@
+/*************************************************************************
+ * Copyright (C) 2009-2014 Tavian Barnes <tavianator@tavianator.com>     *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Triangles.  See
+ * http://tavianator.com/2014/05/a-beautiful-raytriangle-intersection-method/
+ * for a description of the intersection algorithm.
+ */
+
+#include "internal.h"
+#include "dimension/model.h"
+
+/// Optimized ray/triangle intersection test.
+static inline bool
+dmnsn_ray_triangle_intersection(dmnsn_ray l, double *t, double *u, double *v)
+{
+  // See the change of basis in dmnsn_triangle_basis()
+  *t = -l.x0.z/l.n.z;
+  *u = l.x0.x + (*t)*l.n.x;
+  *v = l.x0.y + (*t)*l.n.y;
+  return *t >= 0.0 && *u >= 0.0 && *v >= 0.0 && *u + *v <= 1.0;
+}
+
+/// Triangle intersection callback.
+DMNSN_HOT static bool
+dmnsn_triangle_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
+                               dmnsn_intersection *intersection)
+{
+  double t, u, v;
+  if (dmnsn_ray_triangle_intersection(l, &t, &u, &v)) {
+    intersection->t = t;
+    intersection->normal = dmnsn_z;
+    return true;
+  }
+
+  return false;
+}
+
+/// Triangle inside callback.
+static bool
+dmnsn_triangle_inside_fn(const dmnsn_object *object, dmnsn_vector point)
+{
+  return false;
+}
+
+/// Triangle bounding callback.
+static dmnsn_aabb
+dmnsn_triangle_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  dmnsn_vector a = dmnsn_transform_point(trans, dmnsn_zero);
+  dmnsn_vector b = dmnsn_transform_point(trans, dmnsn_x);
+  dmnsn_vector c = dmnsn_transform_point(trans, dmnsn_y);
+
+  dmnsn_aabb box = dmnsn_new_aabb(a, a);
+  box = dmnsn_aabb_swallow(box, b);
+  box = dmnsn_aabb_swallow(box, c);
+  return box;
+}
+
+/// Triangle vtable.
+static const dmnsn_object_vtable dmnsn_triangle_vtable = {
+  .intersection_fn = dmnsn_triangle_intersection_fn,
+  .inside_fn = dmnsn_triangle_inside_fn,
+  .bounding_fn = dmnsn_triangle_bounding_fn,
+};
+
+/// Smooth triangle type.
+typedef struct {
+  dmnsn_object object;
+  dmnsn_vector na, nab, nac;
+} dmnsn_smooth_triangle;
+
+/// Smooth triangle intersection callback.
+DMNSN_HOT static bool
+dmnsn_smooth_triangle_intersection_fn(const dmnsn_object *object, dmnsn_ray l,
+                                      dmnsn_intersection *intersection)
+{
+  const dmnsn_smooth_triangle *triangle = (const dmnsn_smooth_triangle *)object;
+
+  double t, u, v;
+  if (dmnsn_ray_triangle_intersection(l, &t, &u, &v)) {
+    intersection->t = t;
+    intersection->normal = dmnsn_vector_add(
+      triangle->na,
+      dmnsn_vector_add(
+        dmnsn_vector_mul(u, triangle->nab),
+        dmnsn_vector_mul(v, triangle->nac)
+      )
+    );
+    return true;
+  }
+
+  return false;
+}
+
+/// Smooth triangle vtable.
+static const dmnsn_object_vtable dmnsn_smooth_triangle_vtable = {
+  .intersection_fn = dmnsn_smooth_triangle_intersection_fn,
+  .inside_fn = dmnsn_triangle_inside_fn,
+  .bounding_fn = dmnsn_triangle_bounding_fn,
+};
+
+/// Make a change-of-basis matrix.
+static inline dmnsn_matrix
+dmnsn_triangle_basis(dmnsn_vector vertices[3])
+{
+  // The new vector space has corners at <1, 0, 0>, <0, 1, 0>, and 0,
+  // corresponding to the basis (ab, ac, ab X ac).
+  dmnsn_vector ab = dmnsn_vector_sub(vertices[1], vertices[0]);
+  dmnsn_vector ac = dmnsn_vector_sub(vertices[2], vertices[0]);
+  dmnsn_vector normal = dmnsn_vector_cross(ab, ac);
+  return dmnsn_new_matrix4(ab, ac, normal, vertices[0]);
+}
+
+dmnsn_object *
+dmnsn_new_triangle(dmnsn_pool *pool, dmnsn_vector vertices[3])
+{
+  dmnsn_object *object = dmnsn_new_object(pool);
+  object->vtable = &dmnsn_triangle_vtable;
+  object->intrinsic_trans = dmnsn_triangle_basis(vertices);
+  return object;
+}
+
+dmnsn_object *
+dmnsn_new_smooth_triangle(dmnsn_pool *pool, dmnsn_vector vertices[3], dmnsn_vector normals[3])
+{
+  dmnsn_matrix P = dmnsn_triangle_basis(vertices);
+
+  // Transform the given normals.
+  dmnsn_vector na = dmnsn_vector_normalized(dmnsn_transform_normal(P, normals[0]));
+  dmnsn_vector nb = dmnsn_vector_normalized(dmnsn_transform_normal(P, normals[1]));
+  dmnsn_vector nc = dmnsn_vector_normalized(dmnsn_transform_normal(P, normals[2]));
+
+  dmnsn_smooth_triangle *triangle = DMNSN_PALLOC(pool, dmnsn_smooth_triangle);
+  triangle->na  = na;
+  triangle->nab = dmnsn_vector_sub(nb, na);
+  triangle->nac = dmnsn_vector_sub(nc, na);
+
+  dmnsn_object *object = &triangle->object;
+  dmnsn_init_object(object);
+  object->vtable = &dmnsn_smooth_triangle_vtable;
+  object->intrinsic_trans = P;
+
+  return object;
+}
diff --git a/libdimension/model/objects/triangle_fan.c b/libdimension/model/objects/triangle_fan.c
new file mode 100644
index 0000000..93768a9
--- /dev/null
+++ b/libdimension/model/objects/triangle_fan.c
@@ -0,0 +1,347 @@
+/*************************************************************************
+ * Copyright (C) 2014 Tavian Barnes <tavianator@tavianator.com>     *
+ *                                                                       *
+ * This file is part of The Dimension Library.                           *
+ *                                                                       *
+ * The Dimension Library is free software; you can redistribute it and/  *
+ * or modify it under the terms of the GNU Lesser General Public License *
+ * as published by the Free Software Foundation; either version 3 of the *
+ * License, or (at your option) any later version.                       *
+ *                                                                       *
+ * The Dimension Library is distributed in the hope that it will be      *
+ * useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  *
+ * Lesser General Public License for more details.                       *
+ *                                                                       *
+ * You should have received a copy of the GNU Lesser General Public      *
+ * License along with this program.  If not, see                         *
+ * <http://www.gnu.org/licenses/>.                                       *
+ *************************************************************************/
+
+/**
+ * @file
+ * Triangle fans.  See
+ * http://tavianator.com/2014/05/a-beautiful-raymesh-intersection-algorithm/
+ * for a description of the intersection algorithm.
+ */
+
+#include "internal.h"
+#include "dimension/model.h"
+
+/// Triangle fan type.
+typedef struct {
+  dmnsn_object object;
+  size_t ncoeffs;
+  double coeffs[][6];
+} dmnsn_triangle_fan;
+
+/// Change basis from one triangle to the next.
+static inline dmnsn_vector
+dmnsn_change_basis(const double coeffs[6], dmnsn_vector v)
+{
+  return dmnsn_new_vector(
+    coeffs[0]*v.x + coeffs[1]*v.z + v.y,
+    coeffs[2]*v.x + coeffs[3]*v.z,
+    coeffs[4]*v.x + coeffs[5]*v.z
+  );
+}
+
+/// Change basis from one triangle to the next for a normal vector.
+static inline dmnsn_vector
+dmnsn_change_normal_basis(const double coeffs[6], dmnsn_vector n)
+{
+  return dmnsn_new_vector(
+    coeffs[0]*n.x + coeffs[2]*n.y + coeffs[4]*n.z,
+    n.x,
+    coeffs[1]*n.x + coeffs[3]*n.y + coeffs[5]*n.z
+  );
+}
+
+/// Change basis from one triangle to the next for a ray
+static inline dmnsn_ray
+dmnsn_change_ray_basis(const double coeffs[6], dmnsn_ray l)
+{
+  return dmnsn_new_ray(dmnsn_change_basis(coeffs, l.x0), dmnsn_change_basis(coeffs, l.n));
+}
+
+/// Store the compressed incremental matrix.
+static inline void
+dmnsn_compress_coeffs(double coeffs[6], dmnsn_matrix incremental)
+{
+  coeffs[0] = incremental.n[0][0];
+  coeffs[1] = incremental.n[0][2];
+  coeffs[2] = incremental.n[1][0];
+  coeffs[3] = incremental.n[1][2];
+  coeffs[4] = incremental.n[2][0];
+  coeffs[5] = incremental.n[2][2];
+}
+
+/// Decompress the incremental matrix.
+static inline dmnsn_matrix
+dmnsn_decompress_coeffs(const double coeffs[6])
+{
+  dmnsn_matrix incremental = dmnsn_new_matrix(
+    coeffs[0], 1.0, coeffs[1], 0.0,
+    coeffs[2], 0.0, coeffs[3], 0.0,
+    coeffs[4], 0.0, coeffs[5], 0.0
+  );
+  return incremental;
+}
+
+/// Make a change-of-basis matrix for a triangle.
+static inline dmnsn_matrix
+dmnsn_triangle_basis(dmnsn_vector a, dmnsn_vector ab, dmnsn_vector ac)
+{
+  dmnsn_vector normal = dmnsn_vector_cross(ab, ac);
+  return dmnsn_new_matrix4(ab, ac, normal, a);
+}
+
+/// Optimized ray/triangle intersection test.
+static inline bool
+dmnsn_ray_triangle_intersection(dmnsn_ray l, double *t, double *u, double *v)
+{
+  *t = -l.x0.z/l.n.z;
+  *u = l.x0.x + (*t)*l.n.x;
+  *v = l.x0.y + (*t)*l.n.y;
+  return *t >= 0.0 && *u >= 0.0 && *v >= 0.0 && *u + *v <= 1.0;
+}
+
+/// Triangle fan intersection callback.
+DMNSN_HOT static bool
+dmnsn_triangle_fan_intersection_fn(const dmnsn_object *object, dmnsn_ray l, dmnsn_intersection *intersection)
+{
+  const dmnsn_triangle_fan *fan = (const dmnsn_triangle_fan *)object;
+
+  double t, u, v;
+
+  double best_t = INFINITY;
+  if (dmnsn_ray_triangle_intersection(l, &t, &u, &v)) {
+    best_t = t;
+  }
+
+  dmnsn_vector normal = dmnsn_z;
+  dmnsn_vector best_normal = normal;
+
+  for (size_t i = 0; i < fan->ncoeffs; ++i) {
+    const double *coeffs = fan->coeffs[i];
+    l = dmnsn_change_ray_basis(coeffs, l);
+    normal = dmnsn_change_normal_basis(coeffs, normal);
+
+    if (dmnsn_ray_triangle_intersection(l, &t, &u, &v) && t < best_t) {
+      best_t = t;
+      best_normal = normal;
+    }
+  }
+
+  if (!isinf(best_t)) {
+    intersection->t = t;
+    intersection->normal = best_normal;
+    return true;
+  }
+
+  return false;
+}
+
+/// Triangle fan inside callback.
+static bool
+dmnsn_triangle_fan_inside_fn(const dmnsn_object *object, dmnsn_vector point)
+{
+  return false;
+}
+
+/// Computes the bounding box for the first triangle
+static inline dmnsn_aabb
+dmnsn_bound_first_triangle(dmnsn_matrix trans)
+{
+  dmnsn_vector a = dmnsn_transform_point(trans, dmnsn_zero);
+  dmnsn_vector b = dmnsn_transform_point(trans, dmnsn_x);
+  dmnsn_vector c = dmnsn_transform_point(trans, dmnsn_y);
+
+  dmnsn_aabb box = dmnsn_new_aabb(a, a);
+  box = dmnsn_aabb_swallow(box, b);
+  box = dmnsn_aabb_swallow(box, c);
+
+  return box;
+}
+
+/// Triangle fan bounding callback.
+static dmnsn_aabb
+dmnsn_triangle_fan_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  const dmnsn_triangle_fan *fan = (const dmnsn_triangle_fan *)object;
+
+  dmnsn_aabb box = dmnsn_bound_first_triangle(trans);
+
+  for (size_t i = 0; i < fan->ncoeffs; ++i) {
+    dmnsn_matrix incremental = dmnsn_decompress_coeffs(fan->coeffs[i]);
+    trans = dmnsn_matrix_mul(trans, dmnsn_matrix_inverse(incremental));
+    dmnsn_vector vertex = dmnsn_transform_point(trans, dmnsn_y);
+    box = dmnsn_aabb_swallow(box, vertex);
+  }
+
+  return box;
+}
+
+/// Triangle fan vtable.
+static dmnsn_object_vtable dmnsn_triangle_fan_vtable = {
+  .intersection_fn = dmnsn_triangle_fan_intersection_fn,
+  .inside_fn = dmnsn_triangle_fan_inside_fn,
+  .bounding_fn = dmnsn_triangle_fan_bounding_fn,
+};
+
+/// Smooth triangle fan type.
+typedef struct dmnsn_smooth_triangle_fan {
+  dmnsn_object object;
+  dmnsn_vector na, nab, nac;
+  size_t ncoeffs;
+  double coeffs[][9]; ///< 0-6 is same as dmnsn_triangle_fan, 6-9 is the normal
+} dmnsn_smooth_triangle_fan;
+
+/// Smooth triangle fan intersection callback.
+DMNSN_HOT static bool
+dmnsn_smooth_triangle_fan_intersection_fn(const dmnsn_object *object, dmnsn_ray l, dmnsn_intersection *intersection)
+{
+  const dmnsn_smooth_triangle_fan *fan = (const dmnsn_smooth_triangle_fan *)object;
+
+  dmnsn_vector nab = fan->nab;
+  dmnsn_vector nac = fan->nac;
+
+  double t, u, v;
+
+  double best_t = INFINITY;
+  dmnsn_vector best_normal;
+  if (dmnsn_ray_triangle_intersection(l, &t, &u, &v)) {
+    best_t = t;
+    best_normal = dmnsn_vector_add(dmnsn_vector_mul(u, nab), dmnsn_vector_mul(v, nac));
+  }
+
+  for (size_t i = 0; i < fan->ncoeffs; ++i) {
+    const double *coeffs = fan->coeffs[i];
+    l = dmnsn_change_ray_basis(coeffs, l);
+    nab = nac;
+    nac = dmnsn_new_vector(coeffs[6], coeffs[7], coeffs[8]);
+
+    if (dmnsn_ray_triangle_intersection(l, &t, &u, &v) && t < best_t) {
+      best_t = t;
+      best_normal = dmnsn_vector_add(dmnsn_vector_mul(u, nab), dmnsn_vector_mul(v, nac));
+    }
+  }
+
+  if (!isinf(best_t)) {
+    intersection->t = t;
+    intersection->normal = dmnsn_vector_add(fan->na, best_normal);
+    return true;
+  }
+
+  return false;
+}
+
+/// Smooth triangle fan bounding callback.
+static dmnsn_aabb
+dmnsn_smooth_triangle_fan_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans)
+{
+  const dmnsn_smooth_triangle_fan *fan = (const dmnsn_smooth_triangle_fan *)object;
+
+  dmnsn_aabb box = dmnsn_bound_first_triangle(trans);
+
+  for (size_t i = 0; i < fan->ncoeffs; ++i) {
+    dmnsn_matrix incremental = dmnsn_decompress_coeffs(fan->coeffs[i]);
+    trans = dmnsn_matrix_mul(trans, dmnsn_matrix_inverse(incremental));
+    dmnsn_vector vertex = dmnsn_transform_point(trans, dmnsn_y);
+    box = dmnsn_aabb_swallow(box, vertex);
+  }
+
+  return box;
+}
+
+/// Smooth triangle fan vtable.
+static dmnsn_object_vtable dmnsn_smooth_triangle_fan_vtable = {
+  .intersection_fn = dmnsn_smooth_triangle_fan_intersection_fn,
+  .inside_fn = dmnsn_triangle_fan_inside_fn,
+  .bounding_fn = dmnsn_smooth_triangle_fan_bounding_fn,
+};
+
+dmnsn_object *
+dmnsn_new_triangle_fan(dmnsn_pool *pool, dmnsn_vector vertices[], size_t nvertices)
+{
+  dmnsn_assert(nvertices >= 3, "Not enough vertices for one triangle");
+
+  size_t ncoeffs = nvertices - 3;
+  dmnsn_triangle_fan *fan = dmnsn_palloc(pool, sizeof(dmnsn_triangle_fan) + ncoeffs*sizeof(double[6]));
+  fan->ncoeffs = ncoeffs;
+
+  dmnsn_object *object = &fan->object;
+  dmnsn_init_object(object);
+  object->vtable = &dmnsn_triangle_fan_vtable;
+
+  // Compute the initial matrix and the coefficients
+  dmnsn_vector a = vertices[0];
+  dmnsn_vector ab = dmnsn_vector_sub(vertices[1], a);
+  dmnsn_vector ac = dmnsn_vector_sub(vertices[2], a);
+  dmnsn_matrix P = dmnsn_triangle_basis(a, ab, ac);
+  object->intrinsic_trans = P;
+
+  for (size_t i = 0; i < ncoeffs; ++i) {
+    ab = ac;
+    ac = dmnsn_vector_sub(vertices[i + 3], a);
+
+    dmnsn_matrix newP = dmnsn_triangle_basis(a, ab, ac);
+    dmnsn_matrix incremental = dmnsn_matrix_mul(dmnsn_matrix_inverse(newP), P);
+    dmnsn_compress_coeffs(fan->coeffs[i], incremental);
+
+    P = newP;
+  }
+
+  return object;
+}
+
+dmnsn_object *
+dmnsn_new_smooth_triangle_fan(dmnsn_pool *pool, dmnsn_vector vertices[], dmnsn_vector normals[], size_t nvertices)
+{
+  dmnsn_assert(nvertices >= 3, "Not enough vertices for one triangle");
+
+  size_t ncoeffs = nvertices - 3;
+  dmnsn_smooth_triangle_fan *fan = dmnsn_palloc(pool, sizeof(dmnsn_smooth_triangle_fan) + ncoeffs*sizeof(double[9]));
+  fan->ncoeffs = ncoeffs;
+
+  dmnsn_object *object = &fan->object;
+  dmnsn_init_object(object);
+  object->vtable = &dmnsn_smooth_triangle_fan_vtable;
+
+  // Compute the initial matrix
+  dmnsn_vector a = vertices[0];
+  dmnsn_vector ab = dmnsn_vector_sub(vertices[1], a);
+  dmnsn_vector ac = dmnsn_vector_sub(vertices[2], a);
+  dmnsn_matrix P = dmnsn_triangle_basis(a, ab, ac);
+  dmnsn_matrix Pabc = P;
+  object->intrinsic_trans = P;
+
+  // Transform the first three normals
+  dmnsn_vector na = dmnsn_vector_normalized(dmnsn_transform_normal(P, normals[0]));
+  dmnsn_vector nb = dmnsn_vector_normalized(dmnsn_transform_normal(P, normals[1]));
+  dmnsn_vector nc = dmnsn_vector_normalized(dmnsn_transform_normal(P, normals[2]));
+  fan->na = na;
+  fan->nab = dmnsn_vector_sub(nb, na);
+  fan->nac = dmnsn_vector_sub(nc, na);
+
+  // Compute the coefficients
+  for (size_t i = 0; i < ncoeffs; ++i) {
+    ab = ac;
+    ac = dmnsn_vector_sub(vertices[i + 3], a);
+
+    dmnsn_matrix newP = dmnsn_triangle_basis(a, ab, ac);
+    dmnsn_matrix incremental = dmnsn_matrix_mul(dmnsn_matrix_inverse(newP), P);
+    double *coeffs = fan->coeffs[i];
+    dmnsn_compress_coeffs(coeffs, incremental);
+
+    nc = dmnsn_vector_normalized(dmnsn_transform_normal(Pabc, normals[i + 3]));
+    dmnsn_vector nac = dmnsn_vector_sub(nc, na);
+    coeffs[6] = nac.x;
+    coeffs[7] = nac.y;
+    coeffs[8] = nac.z;
+
+    P = newP;
+  }
+
+  return object;
+}