summaryrefslogtreecommitdiffstats
path: root/libdimension/model/objects
diff options
context:
space:
mode:
Diffstat (limited to 'libdimension/model/objects')
-rw-r--r--libdimension/model/objects/cone.c207
-rw-r--r--libdimension/model/objects/csg.c334
-rw-r--r--libdimension/model/objects/cube.c154
-rw-r--r--libdimension/model/objects/plane.c88
-rw-r--r--libdimension/model/objects/sphere.c116
-rw-r--r--libdimension/model/objects/torus.c173
-rw-r--r--libdimension/model/objects/triangle.c163
-rw-r--r--libdimension/model/objects/triangle_fan.c347
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;
+}