diff options
Diffstat (limited to 'libdimension')
| -rw-r--r-- | libdimension/dimension/objects.h | 13 | ||||
| -rw-r--r-- | libdimension/triangle.c | 6 | ||||
| -rw-r--r-- | libdimension/triangle_fan.c | 234 |
3 files changed, 212 insertions, 41 deletions
diff --git a/libdimension/dimension/objects.h b/libdimension/dimension/objects.h index b328025..e5a39c5 100644 --- a/libdimension/dimension/objects.h +++ b/libdimension/dimension/objects.h @@ -49,8 +49,17 @@ dmnsn_object *dmnsn_new_smooth_triangle(dmnsn_pool *pool, dmnsn_vector vertices[ * @param[in] nvertices The number of vertices. * @return A triangle fan. */ -dmnsn_object * -dmnsn_new_triangle_fan(dmnsn_pool *pool, dmnsn_vector vertices[], size_t nvertices); +dmnsn_object *dmnsn_new_triangle_fan(dmnsn_pool *pool, dmnsn_vector vertices[], size_t nvertices); + +/** + * A smooth triangle fan. + * @param[in] pool The memory pool to allocate from. + * @param[in] vertices The vertices of the fan, starting in the center. + * @param[in] vertices The normal vector for each vertex. + * @param[in] nvertices The number of vertices. + * @return A triangle fan. + */ +dmnsn_object *dmnsn_new_smooth_triangle_fan(dmnsn_pool *pool, dmnsn_vector vertices[], dmnsn_vector normals[], size_t nvertices); /** * A plane. diff --git a/libdimension/triangle.c b/libdimension/triangle.c index 4a70c59..fdd2b96 100644 --- a/libdimension/triangle.c +++ b/libdimension/triangle.c @@ -25,7 +25,7 @@ * for a description of the intersection algorithm. */ -#include "dimension.h" +#include "dimension-internal.h" /// Optimized ray/triangle intersection test. static inline bool @@ -39,7 +39,7 @@ dmnsn_ray_triangle_intersection(dmnsn_line l, double *t, double *u, double *v) } /// Triangle intersection callback. -static bool +DMNSN_HOT static bool dmnsn_triangle_intersection_fn(const dmnsn_object *object, dmnsn_line l, dmnsn_intersection *intersection) { @@ -88,7 +88,7 @@ typedef struct { } dmnsn_smooth_triangle; /// Smooth triangle intersection callback. -static bool +DMNSN_HOT static bool dmnsn_smooth_triangle_intersection_fn(const dmnsn_object *object, dmnsn_line l, dmnsn_intersection *intersection) { diff --git a/libdimension/triangle_fan.c b/libdimension/triangle_fan.c index 0165a51..9940614 100644 --- a/libdimension/triangle_fan.c +++ b/libdimension/triangle_fan.c @@ -25,7 +25,7 @@ * for a description of the intersection algorithm. */ -#include "dimension.h" +#include "dimension-internal.h" /// Triangle fan type. typedef struct { @@ -56,18 +56,65 @@ dmnsn_change_normal_basis(const double coeffs[6], dmnsn_vector n) ); } -/// Triangle intersection callback. -static bool +/// Change basis from one triangle to the next for a line +static inline dmnsn_line +dmnsn_change_line_basis(const double coeffs[6], dmnsn_line l) +{ + return dmnsn_new_line(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_line 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_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 t, u, v; double best_t = INFINITY; - if (t >= 0.0 && u >= 0.0 && v >= 0.0 && u + v <= 1.0) { + if (dmnsn_ray_triangle_intersection(l, &t, &u, &v)) { best_t = t; } @@ -76,14 +123,10 @@ dmnsn_triangle_fan_intersection_fn(const dmnsn_object *object, dmnsn_line l, dmn 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)); + l = dmnsn_change_line_basis(coeffs, l); 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) { + if (dmnsn_ray_triangle_intersection(l, &t, &u, &v) && t < best_t) { best_t = t; best_normal = normal; } @@ -105,12 +148,10 @@ 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) +/// Computes the bounding box for the first triangle +static inline dmnsn_bounding_box +dmnsn_bound_first_triangle(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); @@ -119,16 +160,22 @@ dmnsn_triangle_fan_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans) box = dmnsn_bounding_box_swallow(box, b); box = dmnsn_bounding_box_swallow(box, c); - dmnsn_matrix P = trans; + return box; +} + +/// 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_bounding_box box = dmnsn_bound_first_triangle(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); + 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_bounding_box_swallow(box, vertex); } return box; @@ -141,6 +188,78 @@ static dmnsn_object_vtable dmnsn_triangle_fan_vtable = { .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_line 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_line_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_bounding_box +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_bounding_box 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_bounding_box_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) { @@ -158,24 +277,67 @@ dmnsn_new_triangle_fan(dmnsn_pool *pool, dmnsn_vector vertices[], size_t nvertic 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); + 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); - normal = dmnsn_vector_cross(ab, ac); - dmnsn_matrix newP = dmnsn_new_matrix4(ab, ac, normal, 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]; - 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]; + 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; } |