diff options
Diffstat (limited to 'libdimension/model/objects')
| -rw-r--r-- | libdimension/model/objects/cone.c | 207 | ||||
| -rw-r--r-- | libdimension/model/objects/csg.c | 334 | ||||
| -rw-r--r-- | libdimension/model/objects/cube.c | 154 | ||||
| -rw-r--r-- | libdimension/model/objects/plane.c | 88 | ||||
| -rw-r--r-- | libdimension/model/objects/sphere.c | 116 | ||||
| -rw-r--r-- | libdimension/model/objects/torus.c | 173 | ||||
| -rw-r--r-- | libdimension/model/objects/triangle.c | 163 | ||||
| -rw-r--r-- | libdimension/model/objects/triangle_fan.c | 347 |
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; +} |