diff options
Diffstat (limited to 'libdimension/triangle_fan.c')
-rw-r--r-- | libdimension/triangle_fan.c | 186 |
1 files changed, 186 insertions, 0 deletions
diff --git a/libdimension/triangle_fan.c b/libdimension/triangle_fan.c new file mode 100644 index 0000000..48a64ff --- /dev/null +++ b/libdimension/triangle_fan.c @@ -0,0 +1,186 @@ +/************************************************************************* + * 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 "dimension.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 + ); +} + +/** Triangle intersection callback. */ +static bool +dmnsn_triangle_fan_intersection_fn(const dmnsn_object *object, dmnsn_line l, dmnsn_intersection *intersection) +{ + const dmnsn_triangle_fan *fan = (const dmnsn_triangle_fan *)object; + + double t = -l.x0.z/l.n.z; + double u = l.x0.x + t*l.n.x; + double v = l.x0.y + t*l.n.y; + + double best_t = INFINITY; + if (t >= 0.0 && u >= 0.0 && v >= 0.0 && u + v <= 1.0) { + 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_new_line(dmnsn_change_basis(coeffs, l.x0), dmnsn_change_basis(coeffs, l.n)); + normal = dmnsn_change_normal_basis(coeffs, normal); + + t = -l.x0.z/l.n.z; + u = l.x0.x + t*l.n.x; + v = l.x0.y + t*l.n.y; + + if (t >= 0.0 && t < best_t && u >= 0.0 && v >= 0.0 && u + v <= 1.0) { + 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; +} + +/** Triangle fan bounding callback. */ +static dmnsn_bounding_box +dmnsn_triangle_fan_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans) +{ + const dmnsn_triangle_fan *fan = (const dmnsn_triangle_fan *)object; + + 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_bounding_box box = dmnsn_new_bounding_box(a, a); + box = dmnsn_bounding_box_swallow(box, b); + box = dmnsn_bounding_box_swallow(box, c); + + dmnsn_matrix P = trans; + + for (size_t i = 0; i < fan->ncoeffs; ++i) { + const double *coeffs = fan->coeffs[i]; + 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); + P = dmnsn_matrix_mul(P, dmnsn_matrix_inverse(incremental)); + dmnsn_vector d = dmnsn_transform_point(P, dmnsn_y); + box = dmnsn_bounding_box_swallow(box, d); + } + + 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, +}; + +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; + object->bounding_box.min = dmnsn_zero; + object->bounding_box.max = dmnsn_new_vector(1.0, 1.0, 0.0); + + /* 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_vector normal = dmnsn_vector_cross(ab, ac); + dmnsn_matrix P = dmnsn_new_matrix4(ab, ac, normal, a); + object->intrinsic_trans = P; + + for (size_t i = 0; i < ncoeffs; ++i) { + ab = ac; + ac = dmnsn_vector_sub(vertices[i + 3], a); + normal = dmnsn_vector_cross(ab, ac); + + dmnsn_matrix newP = dmnsn_new_matrix4(ab, ac, normal, a); + dmnsn_matrix incremental = dmnsn_matrix_mul(dmnsn_matrix_inverse(newP), P); + double *coeffs = fan->coeffs[i]; + 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]; + + P = newP; + } + + return object; +} |