diff options
Diffstat (limited to 'libdimension')
23 files changed, 386 insertions, 374 deletions
diff --git a/libdimension/bench/geometry.c b/libdimension/bench/geometry.c index 59d27e3..ae30d17 100644 --- a/libdimension/bench/geometry.c +++ b/libdimension/bench/geometry.c @@ -107,12 +107,6 @@ main(void) }); printf("dmnsn_vector_mul(): %ld\n", sandglass.grains); - // dmnsn_vector_div() - sandglass_bench_fine(&sandglass, { - vector = dmnsn_vector_div(vector, 2.0); - }); - printf("dmnsn_vector_div(): %ld\n", sandglass.grains); - // dmnsn_vector_cross() sandglass_bench_fine(&sandglass, { vector = dmnsn_vector_cross(vector, vector2); diff --git a/libdimension/bench/prtree.c b/libdimension/bench/prtree.c index 6b5e7c1..aa29317 100644 --- a/libdimension/bench/prtree.c +++ b/libdimension/bench/prtree.c @@ -29,7 +29,7 @@ static bool dmnsn_fake_intersection_fn(const dmnsn_object *object, dmnsn_ray ray, dmnsn_intersection *intersection) { - intersection->t = (object->aabb.min.z - ray.x0.z)/ray.n.z; + intersection->t = (object->aabb.min.Z - ray.x0.Z)/ray.n.Z; intersection->normal = dmnsn_x; return true; } @@ -45,13 +45,10 @@ dmnsn_fake_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans) { dmnsn_vector a, b; - a.x = 2.0*((double)rand())/RAND_MAX - 1.0; - a.y = 2.0*((double)rand())/RAND_MAX - 1.0; - a.z = 2.0*((double)rand())/RAND_MAX - 1.0; - - b.x = 2.0*((double)rand())/RAND_MAX - 1.0; - b.y = 2.0*((double)rand())/RAND_MAX - 1.0; - b.z = 2.0*((double)rand())/RAND_MAX - 1.0; + for (unsigned int i = 0; i < 3; ++i) { + a.n[i] = 2.0*((double)rand())/RAND_MAX - 1.0; + b.n[i] = 2.0*((double)rand())/RAND_MAX - 1.0; + } return dmnsn_new_aabb(dmnsn_vector_min(a, b), dmnsn_vector_max(a, b)); } diff --git a/libdimension/bvh/bvh.c b/libdimension/bvh/bvh.c index eab2c28..9f460a2 100644 --- a/libdimension/bvh/bvh.c +++ b/libdimension/bvh/bvh.c @@ -186,7 +186,7 @@ dmnsn_optimize_ray(dmnsn_ray ray) { dmnsn_optimized_ray optray = { .x0 = ray.x0, - .n_inv = dmnsn_new_vector(1.0/ray.n.x, 1.0/ray.n.y, 1.0/ray.n.z) + .n_inv = dmnsn_new_vector(1.0/ray.n.X, 1.0/ray.n.Y, 1.0/ray.n.Z) }; return optray; } @@ -202,20 +202,20 @@ dmnsn_ray_box_intersection(dmnsn_optimized_ray optray, dmnsn_aabb box, double t) // or tmax == -inf, while rays inside the box will have tmin and tmax // unchanged. - double tx1 = (box.min.x - optray.x0.x)*optray.n_inv.x; - double tx2 = (box.max.x - optray.x0.x)*optray.n_inv.x; + double tx1 = (box.min.X - optray.x0.X)*optray.n_inv.X; + double tx2 = (box.max.X - optray.x0.X)*optray.n_inv.X; double tmin = dmnsn_min(tx1, tx2); double tmax = dmnsn_max(tx1, tx2); - double ty1 = (box.min.y - optray.x0.y)*optray.n_inv.y; - double ty2 = (box.max.y - optray.x0.y)*optray.n_inv.y; + double ty1 = (box.min.Y - optray.x0.Y)*optray.n_inv.Y; + double ty2 = (box.max.Y - optray.x0.Y)*optray.n_inv.Y; tmin = dmnsn_max(tmin, dmnsn_min(ty1, ty2)); tmax = dmnsn_min(tmax, dmnsn_max(ty1, ty2)); - double tz1 = (box.min.z - optray.x0.z)*optray.n_inv.z; - double tz2 = (box.max.z - optray.x0.z)*optray.n_inv.z; + double tz1 = (box.min.Z - optray.x0.Z)*optray.n_inv.Z; + double tz2 = (box.max.Z - optray.x0.Z)*optray.n_inv.Z; tmin = dmnsn_max(tmin, dmnsn_min(tz1, tz2)); tmax = dmnsn_min(tmax, dmnsn_max(tz1, tz2)); diff --git a/libdimension/bvh/prtree.c b/libdimension/bvh/prtree.c index c8e4e54..84fca58 100644 --- a/libdimension/bvh/prtree.c +++ b/libdimension/bvh/prtree.c @@ -74,48 +74,48 @@ enum { static int dmnsn_xmin_comp(const void *l, const void *r) { - double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.x; - double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.x; + double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.X; + double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.X; return (lval > rval) - (lval < rval); } static int dmnsn_ymin_comp(const void *l, const void *r) { - double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.y; - double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.y; + double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.Y; + double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.Y; return (lval > rval) - (lval < rval); } static int dmnsn_zmin_comp(const void *l, const void *r) { - double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.z; - double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.z; + double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.min.Z; + double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.min.Z; return (lval > rval) - (lval < rval); } static int dmnsn_xmax_comp(const void *l, const void *r) { - double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.x; - double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.x; + double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.X; + double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.X; return (lval < rval) - (lval > rval); } static int dmnsn_ymax_comp(const void *l, const void *r) { - double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.y; - double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.y; + double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.Y; + double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.Y; return (lval < rval) - (lval > rval); } static int dmnsn_zmax_comp(const void *l, const void *r) { - double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.z; - double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.z; + double lval = (*(const dmnsn_colored_prnode **)l)->node->aabb.max.Z; + double rval = (*(const dmnsn_colored_prnode **)r)->node->aabb.max.Z; return (lval < rval) - (lval > rval); } diff --git a/libdimension/dimension.h b/libdimension/dimension.h index db67a48..2489669 100644 --- a/libdimension/dimension.h +++ b/libdimension/dimension.h @@ -34,9 +34,6 @@ * all rendering-related tasks for Dimension. */ -#ifndef DMNSN_H -#define DMNSN_H - /* Include all modules. */ #include <dimension/base.h> #include <dimension/platform.h> @@ -47,5 +44,3 @@ #include <dimension/pattern.h> #include <dimension/model.h> #include <dimension/render.h> - -#endif /* DMNSN_H */ diff --git a/libdimension/dimension/math.h b/libdimension/dimension/math.h index 603373f..e110ffb 100644 --- a/libdimension/dimension/math.h +++ b/libdimension/dimension/math.h @@ -43,3 +43,14 @@ extern "C" { #endif #endif /* DMNSN_MATH_H */ + +#if defined(DMNSN_SHORT_NAMES) && !defined(DMNSN_SHORT_NAMES_DEFINED) + #define DMNSN_SHORT_NAMES_DEFINED + + /** Short name for the \a x component of a \c dmnsn_vector. */ + #define X n[0] + /** Short name for the \a y component of a \c dmnsn_vector. */ + #define Y n[1] + /** Short name for the \a z component of a \c dmnsn_vector. */ + #define Z n[2] +#endif diff --git a/libdimension/dimension/math/aabb.h b/libdimension/dimension/math/aabb.h index 14cc575..4eb7870 100644 --- a/libdimension/dimension/math/aabb.h +++ b/libdimension/dimension/math/aabb.h @@ -57,8 +57,8 @@ DMNSN_INLINE dmnsn_aabb dmnsn_zero_aabb(void) { dmnsn_aabb box = { - { DMNSN_INFINITY, DMNSN_INFINITY, DMNSN_INFINITY }, - { -DMNSN_INFINITY, -DMNSN_INFINITY, -DMNSN_INFINITY } + { { DMNSN_INFINITY, DMNSN_INFINITY, DMNSN_INFINITY } }, + { { -DMNSN_INFINITY, -DMNSN_INFINITY, -DMNSN_INFINITY } } }; return box; } @@ -68,8 +68,8 @@ DMNSN_INLINE dmnsn_aabb dmnsn_infinite_aabb(void) { dmnsn_aabb box = { - { -DMNSN_INFINITY, -DMNSN_INFINITY, -DMNSN_INFINITY }, - { DMNSN_INFINITY, DMNSN_INFINITY, DMNSN_INFINITY } + { { -DMNSN_INFINITY, -DMNSN_INFINITY, -DMNSN_INFINITY } }, + { { DMNSN_INFINITY, DMNSN_INFINITY, DMNSN_INFINITY } } }; return box; } @@ -94,15 +94,26 @@ dmnsn_symmetric_aabb(dmnsn_vector r) DMNSN_INLINE bool dmnsn_aabb_contains(dmnsn_aabb box, dmnsn_vector p) { - return (p.x >= box.min.x && p.y >= box.min.y && p.z >= box.min.z) - && (p.x <= box.max.x && p.y <= box.max.y && p.z <= box.max.z); + for (unsigned int i = 0; i < 3; ++i) { + if (p.n[i] < box.min.n[i]) { + return false; + } + } + + for (unsigned int i = 0; i < 3; ++i) { + if (p.n[i] > box.max.n[i]) { + return false; + } + } + + return true; } /** Return whether a bounding box is infinite. */ DMNSN_INLINE bool dmnsn_aabb_is_infinite(dmnsn_aabb box) { - return box.min.x == -DMNSN_INFINITY; + return box.min.n[0] == -DMNSN_INFINITY; } /** diff --git a/libdimension/dimension/math/matrix.h b/libdimension/dimension/math/matrix.h index 7471bf5..e67121e 100644 --- a/libdimension/dimension/math/matrix.h +++ b/libdimension/dimension/math/matrix.h @@ -39,152 +39,218 @@ typedef struct dmnsn_matrix { "[%g\t%g\t%g\t%g]\n" \ "[%g\t%g\t%g\t%g]" /** The appropriate arguements to printf() a matrix. */ -#define DMNSN_MATRIX_PRINTF(m) \ - (m).n[0][0], (m).n[0][1], (m).n[0][2], (m).n[0][3], \ - (m).n[1][0], (m).n[1][1], (m).n[1][2], (m).n[1][3], \ - (m).n[2][0], (m).n[2][1], (m).n[2][2], (m).n[2][3], \ +#define DMNSN_MATRIX_PRINTF(M) \ + (M).n[0][0], (M).n[0][1], (M).n[0][2], (M).n[0][3], \ + (M).n[1][0], (M).n[1][1], (M).n[1][2], (M).n[1][3], \ + (M).n[2][0], (M).n[2][1], (M).n[2][2], (M).n[2][3], \ 0.0, 0.0, 0.0, 1.0 -/** Construct a new transformation matrix. */ +/** Create a transformation matrix. */ DMNSN_INLINE dmnsn_matrix -dmnsn_new_matrix(double a0, double a1, double a2, double a3, - double b0, double b1, double b2, double b3, - double c0, double c1, double c2, double c3) +dmnsn_new_matrix(double a0, double b0, double c0, double d0, + double a1, double b1, double c1, double d1, + double a2, double b2, double c2, double d2) { - dmnsn_matrix m = { { { a0, a1, a2, a3 }, - { b0, b1, b2, b3 }, - { c0, c1, c2, c3 } } }; - return m; + dmnsn_matrix M = { + { + { a0, b0, c0, d0 }, + { a1, b1, c1, d1 }, + { a2, b2, c2, d2 } + } + }; + return M; } -/** Construct a new transformation matrix from column vectors. */ +/** Create a transformation matrix from column vectors. */ DMNSN_INLINE dmnsn_matrix -dmnsn_new_matrix4(dmnsn_vector a, dmnsn_vector b, dmnsn_vector c, - dmnsn_vector d) +dmnsn_new_matrix4(dmnsn_vector a, dmnsn_vector b, dmnsn_vector c, dmnsn_vector d) { - dmnsn_matrix m = { { { a.x, b.x, c.x, d.x }, - { a.y, b.y, c.y, d.y }, - { a.z, b.z, c.z, d.z } } }; - return m; + dmnsn_matrix M; + + unsigned int i; + for (i = 0; i < 3; ++i) { + M.n[i][0] = a.n[i]; + } + for (i = 0; i < 3; ++i) { + M.n[i][1] = b.n[i]; + } + for (i = 0; i < 3; ++i) { + M.n[i][2] = c.n[i]; + } + for (i = 0; i < 3; ++i) { + M.n[i][3] = d.n[i]; + } + + return M; } /** Extract column vectors from a matrix. */ DMNSN_INLINE dmnsn_vector dmnsn_matrix_column(dmnsn_matrix M, unsigned int i) { - return dmnsn_new_vector(M.n[0][i], M.n[1][i], M.n[2][i]); + dmnsn_vector v; + unsigned int j; + for (j = 0; j < 3; ++j) { + v.n[j] = M.n[j][i]; + } + return v; } /** Return the identity matrix. */ -dmnsn_matrix dmnsn_identity_matrix(void); +DMNSN_INLINE dmnsn_matrix +dmnsn_identity_matrix(void) +{ + return dmnsn_new_matrix( + 1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0 + ); +} /** - * A scale transformation. - * @param[in] s A vector with components representing the scaling factor in - * each axis. - * @return The transformation matrix. + * Return a scale transformation. + * @param[in] s A vector with components representing the scaling factor in + * each axis. */ -dmnsn_matrix dmnsn_scale_matrix(dmnsn_vector s); +DMNSN_INLINE dmnsn_matrix +dmnsn_scale_matrix(dmnsn_vector s) +{ + return dmnsn_new_matrix( + s.n[0], 0.0, 0.0, 0.0, + 0.0, s.n[1], 0.0, 0.0, + 0.0, 0.0, s.n[2], 0.0 + ); +} + /** - * A translation. - * @param[in] d The vector to translate by. - * @return The transformation matrix. + * Set \p M to a translation matrix. + * @param[in] d The vector to translate by. */ -dmnsn_matrix dmnsn_translation_matrix(dmnsn_vector d); +DMNSN_INLINE dmnsn_matrix +dmnsn_translation_matrix(dmnsn_vector d) +{ + return dmnsn_new_matrix( + 1.0, 0.0, 0.0, d.n[0], + 0.0, 1.0, 0.0, d.n[1], + 0.0, 0.0, 1.0, d.n[2] + ); +} + /** - * A left-handed rotation. - * @param[in] theta A vector representing an axis and angle. - * @f$ axis = \vec{\theta}/|\vec{\theta}| @f$, - * @f$ angle = |\vec{\theta}| @f$ - * @return The transformation matrix. + * Return a rotation matrix. + * @param[in] theta A vector representing an axis and angle. + * @f$ axis = \vec{\theta}/|\vec{\theta}| @f$, + * @f$ angle = |\vec{\theta}| @f$ */ dmnsn_matrix dmnsn_rotation_matrix(dmnsn_vector theta); /** - * An alignment matrix. + * Return an alignment matrix. * @param[in] from The initial vector. * @param[in] to The desired direction. * @param[in] axis1 The first axis about which to rotate. * @param[in] axis2 The second axis about which to rotate. - * @return A transformation matrix that will rotate \p from to \p to. */ -dmnsn_matrix dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to, - dmnsn_vector axis1, dmnsn_vector axis2); +dmnsn_matrix dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis1, dmnsn_vector axis2); /** Invert a matrix. */ -dmnsn_matrix dmnsn_matrix_inverse(dmnsn_matrix A); +dmnsn_matrix dmnsn_matrix_inverse(dmnsn_matrix M); /** Multiply two matricies. */ -dmnsn_matrix dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs); +DMNSN_INLINE dmnsn_matrix dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs) +{ + dmnsn_matrix M; + + unsigned int i, j, k; + for (i = 0; i < 3; ++i) { + for (j = 0; j < 3; ++j) { + M.n[i][j] = 0.0; + for (k = 0; k < 3; ++k) { + M.n[i][j] += lhs.n[i][k]*rhs.n[k][j]; + } + } + + M.n[i][3] = lhs.n[i][3]; + for (k = 0; k < 3; ++k) { + M.n[i][3] += lhs.n[i][k]*rhs.n[k][3]; + } + } + + return M; +} /** Transform a point by a matrix. */ DMNSN_INLINE dmnsn_vector -dmnsn_transform_point(dmnsn_matrix T, dmnsn_vector v) +dmnsn_transform_point(dmnsn_matrix M, dmnsn_vector v) { - /* 9 multiplications, 9 additions */ dmnsn_vector r; - r.x = T.n[0][0]*v.x + T.n[0][1]*v.y + T.n[0][2]*v.z + T.n[0][3]; - r.y = T.n[1][0]*v.x + T.n[1][1]*v.y + T.n[1][2]*v.z + T.n[1][3]; - r.z = T.n[2][0]*v.x + T.n[2][1]*v.y + T.n[2][2]*v.z + T.n[2][3]; + unsigned int i, j; + for (i = 0; i < 3; ++i) { + r.n[i] = M.n[i][3]; + for (j = 0; j < 3; ++j) { + r.n[i] += M.n[i][j]*v.n[j]; + } + } return r; } /** Transform a direction by a matrix. */ DMNSN_INLINE dmnsn_vector -dmnsn_transform_direction(dmnsn_matrix T, dmnsn_vector v) +dmnsn_transform_direction(dmnsn_matrix M, dmnsn_vector v) { - /* 9 multiplications, 6 additions */ dmnsn_vector r; - r.x = T.n[0][0]*v.x + T.n[0][1]*v.y + T.n[0][2]*v.z; - r.y = T.n[1][0]*v.x + T.n[1][1]*v.y + T.n[1][2]*v.z; - r.z = T.n[2][0]*v.x + T.n[2][1]*v.y + T.n[2][2]*v.z; + unsigned int i, j; + for (i = 0; i < 3; ++i) { + r.n[i] = 0.0; + for (j = 0; j < 3; ++j) { + r.n[i] += M.n[i][j]*v.n[j]; + } + } return r; } /** * Transform a pseudovector by a matrix. - * @param[in] Tinv The inverse of the transformation matrix. - * @param[in] v The pseudovector to transform + * @param[in] Minv The inverse of the transformation matrix. + * @param[in] v The pseudovector to transform. * @return The transformed pseudovector. */ DMNSN_INLINE dmnsn_vector -dmnsn_transform_normal(dmnsn_matrix Tinv, dmnsn_vector v) +dmnsn_transform_normal(dmnsn_matrix Minv, dmnsn_vector v) { - /* Multiply by the transpose of the inverse - (9 multiplications, 6 additions) */ + /* Multiply by the transpose of the inverse */ dmnsn_vector r; - r.x = Tinv.n[0][0]*v.x + Tinv.n[1][0]*v.y + Tinv.n[2][0]*v.z; - r.y = Tinv.n[0][1]*v.x + Tinv.n[1][1]*v.y + Tinv.n[2][1]*v.z; - r.z = Tinv.n[0][2]*v.x + Tinv.n[1][2]*v.y + Tinv.n[2][2]*v.z; + unsigned int i, j; + for (i = 0; i < 3; ++i) { + r.n[i] = 0.0; + for (j = 0; j < 3; ++j) { + r.n[i] += Minv.n[j][i]*v.n[j]; + } + } return r; } -/** - * Transform a ray by a matrix. - * \f$ n' = T(l.\vec{x_0} + l.\vec{n}) - T(l.\vec{x_0}) \f$, - * \f$ \vec{x_0}' = T(l.\vec{x_0}) \f$ - */ +/** Transform a ray by a matrix. */ DMNSN_INLINE dmnsn_ray -dmnsn_transform_ray(dmnsn_matrix T, dmnsn_ray l) +dmnsn_transform_ray(dmnsn_matrix M, dmnsn_ray l) { - /* 18 multiplications, 15 additions */ dmnsn_ray ret; - ret.x0 = dmnsn_transform_point(T, l.x0); - ret.n = dmnsn_transform_direction(T, l.n); + ret.x0 = dmnsn_transform_point(M, l.x0); + ret.n = dmnsn_transform_direction(M, l.n); return ret; } /** Transform a bounding box by a matrix. */ -dmnsn_aabb dmnsn_transform_aabb(dmnsn_matrix T, dmnsn_aabb box); +dmnsn_aabb dmnsn_transform_aabb(dmnsn_matrix M, dmnsn_aabb box); /** Return whether a matrix contains any NaN components. */ DMNSN_INLINE bool -dmnsn_matrix_isnan(dmnsn_matrix m) +dmnsn_matrix_isnan(dmnsn_matrix M) { - size_t i, j; + unsigned int i, j; for (i = 0; i < 3; ++i) { for (j = 0; j < 4; ++j) { - if (dmnsn_isnan(m.n[i][j])) { + if (dmnsn_isnan(M.n[i][j])) { return true; } } diff --git a/libdimension/dimension/math/vector.h b/libdimension/dimension/math/vector.h index 8eacee9..90b9a3d 100644 --- a/libdimension/dimension/math/vector.h +++ b/libdimension/dimension/math/vector.h @@ -32,26 +32,24 @@ /** A vector in 3 dimensions. */ typedef struct dmnsn_vector { - double x; /**< The x component. */ - double y; /**< The y component. */ - double z; /**< The z component. */ + double n[3]; /**< The components. */ } dmnsn_vector; /** A standard format string for vectors. */ #define DMNSN_VECTOR_FORMAT "<%g, %g, %g>" /** The appropriate arguements to printf() a vector. */ -#define DMNSN_VECTOR_PRINTF(v) (v).x, (v).y, (v).z +#define DMNSN_VECTOR_PRINTF(v) (v).n[0], (v).n[1], (v).n[2] /* Constants */ /** The zero vector. */ -static const dmnsn_vector dmnsn_zero = { 0.0, 0.0, 0.0 }; +static const dmnsn_vector dmnsn_zero = { { 0.0, 0.0, 0.0 } }; /** The x vector. */ -static const dmnsn_vector dmnsn_x = { 1.0, 0.0, 0.0 }; +static const dmnsn_vector dmnsn_x = { { 1.0, 0.0, 0.0 } }; /** The y vector. */ -static const dmnsn_vector dmnsn_y = { 0.0, 1.0, 0.0 }; +static const dmnsn_vector dmnsn_y = { { 0.0, 1.0, 0.0 } }; /** The z vector. */ -static const dmnsn_vector dmnsn_z = { 0.0, 0.0, 1.0 }; +static const dmnsn_vector dmnsn_z = { { 0.0, 0.0, 1.0 } }; /* Shorthand for vector construction */ @@ -59,7 +57,7 @@ static const dmnsn_vector dmnsn_z = { 0.0, 0.0, 1.0 }; DMNSN_INLINE dmnsn_vector dmnsn_new_vector(double x, double y, double z) { - dmnsn_vector v = { x, y, z }; + dmnsn_vector v = { { x, y, z } }; return v; } @@ -69,8 +67,11 @@ dmnsn_new_vector(double x, double y, double z) DMNSN_INLINE dmnsn_vector dmnsn_vector_negate(dmnsn_vector rhs) { - /* 3 negations */ - dmnsn_vector v = { -rhs.x, -rhs.y, -rhs.z }; + dmnsn_vector v; + unsigned int i; + for (i = 0; i < 3; ++i) { + v.n[i] = -rhs.n[i]; + } return v; } @@ -78,8 +79,11 @@ dmnsn_vector_negate(dmnsn_vector rhs) DMNSN_INLINE dmnsn_vector dmnsn_vector_add(dmnsn_vector lhs, dmnsn_vector rhs) { - /* 3 additions */ - dmnsn_vector v = { lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z }; + dmnsn_vector v; + unsigned int i; + for (i = 0; i < 3; ++i) { + v.n[i] = lhs.n[i] + rhs.n[i]; + } return v; } @@ -87,8 +91,11 @@ dmnsn_vector_add(dmnsn_vector lhs, dmnsn_vector rhs) DMNSN_INLINE dmnsn_vector dmnsn_vector_sub(dmnsn_vector lhs, dmnsn_vector rhs) { - /* 3 additions */ - dmnsn_vector v = { lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z }; + dmnsn_vector v; + unsigned int i; + for (i = 0; i < 3; ++i) { + v.n[i] = lhs.n[i] - rhs.n[i]; + } return v; } @@ -96,17 +103,11 @@ dmnsn_vector_sub(dmnsn_vector lhs, dmnsn_vector rhs) DMNSN_INLINE dmnsn_vector dmnsn_vector_mul(double lhs, dmnsn_vector rhs) { - /* 3 multiplications */ - dmnsn_vector v = { lhs*rhs.x, lhs*rhs.y, lhs*rhs.z }; - return v; -} - -/** Divide a vector by a scalar. */ -DMNSN_INLINE dmnsn_vector -dmnsn_vector_div(dmnsn_vector lhs, double rhs) -{ - /* 3 divisions */ - dmnsn_vector v = { lhs.x/rhs, lhs.y/rhs, lhs.z/rhs }; + dmnsn_vector v; + unsigned int i; + for (i = 0; i < 3; ++i) { + v.n[i] = lhs*rhs.n[i]; + } return v; } @@ -114,18 +115,22 @@ dmnsn_vector_div(dmnsn_vector lhs, double rhs) DMNSN_INLINE double dmnsn_vector_dot(dmnsn_vector lhs, dmnsn_vector rhs) { - /* 3 multiplications, 2 additions */ - return lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z; + double result = 0.0; + unsigned int i; + for (i = 0; i < 3; ++i) { + result += lhs.n[i]*rhs.n[i]; + } + return result; } /** Return the cross product of two vectors. */ DMNSN_INLINE dmnsn_vector dmnsn_vector_cross(dmnsn_vector lhs, dmnsn_vector rhs) { - /* 6 multiplications, 3 additions */ - dmnsn_vector v = { lhs.y*rhs.z - lhs.z*rhs.y, - lhs.z*rhs.x - lhs.x*rhs.z, - lhs.x*rhs.y - lhs.y*rhs.x }; + dmnsn_vector v; + v.n[0] = lhs.n[1]*rhs.n[2] - lhs.n[2]*rhs.n[1]; + v.n[1] = lhs.n[2]*rhs.n[0] - lhs.n[0]*rhs.n[2]; + v.n[2] = lhs.n[0]*rhs.n[1] - lhs.n[1]*rhs.n[0]; return v; } @@ -133,7 +138,6 @@ dmnsn_vector_cross(dmnsn_vector lhs, dmnsn_vector rhs) DMNSN_INLINE dmnsn_vector dmnsn_vector_proj(dmnsn_vector u, dmnsn_vector d) { - /* 1 division, 9 multiplications, 4 additions */ return dmnsn_vector_mul(dmnsn_vector_dot(u, d)/dmnsn_vector_dot(d, d), d); } @@ -141,7 +145,6 @@ dmnsn_vector_proj(dmnsn_vector u, dmnsn_vector d) DMNSN_INLINE double dmnsn_vector_norm(dmnsn_vector n) { - /* 1 sqrt, 3 multiplications, 2 additions */ return sqrt(dmnsn_vector_dot(n, n)); } @@ -149,35 +152,42 @@ dmnsn_vector_norm(dmnsn_vector n) DMNSN_INLINE dmnsn_vector dmnsn_vector_normalized(dmnsn_vector n) { - /* 1 sqrt, 3 divisions, 3 multiplications, 2 additions */ - return dmnsn_vector_div(n, dmnsn_vector_norm(n)); + return dmnsn_vector_mul(1.0/dmnsn_vector_norm(n), n); } /** Return the component-wise minimum of two vectors. */ DMNSN_INLINE dmnsn_vector dmnsn_vector_min(dmnsn_vector a, dmnsn_vector b) { - return dmnsn_new_vector( - dmnsn_min(a.x, b.x), - dmnsn_min(a.y, b.y), - dmnsn_min(a.z, b.z) - ); + dmnsn_vector v; + unsigned int i; + for (i = 0; i < 3; ++i) { + v.n[i] = dmnsn_min(a.n[i], b.n[i]); + } + return v; } /** Return the component-wise maximum of two vectors. */ DMNSN_INLINE dmnsn_vector dmnsn_vector_max(dmnsn_vector a, dmnsn_vector b) { - return dmnsn_new_vector( - dmnsn_max(a.x, b.x), - dmnsn_max(a.y, b.y), - dmnsn_max(a.z, b.z) - ); + dmnsn_vector v; + unsigned int i; + for (i = 0; i < 3; ++i) { + v.n[i] = dmnsn_max(a.n[i], b.n[i]); + } + return v; } /** Return whether a vector contains any NaN components. */ DMNSN_INLINE bool dmnsn_vector_isnan(dmnsn_vector v) { - return dmnsn_isnan(v.x) || dmnsn_isnan(v.y) || dmnsn_isnan(v.z); + unsigned int i; + for (i = 0; i < 3; ++i) { + if (dmnsn_isnan(v.n[i])) { + return true; + } + } + return false; } diff --git a/libdimension/internal.h b/libdimension/internal.h index 3db2612..4b1910f 100644 --- a/libdimension/internal.h +++ b/libdimension/internal.h @@ -27,6 +27,8 @@ #ifndef DMNSN_INTERNAL_H #define DMNSN_INTERNAL_H +#define DMNSN_SHORT_NAMES + #include "dimension/base.h" #include "internal/compiler.h" diff --git a/libdimension/math/matrix.c b/libdimension/math/matrix.c index 25590d8..f0050aa 100644 --- a/libdimension/math/matrix.c +++ b/libdimension/math/matrix.c @@ -27,53 +27,22 @@ #include "dimension/math.h" #include <math.h> -// Identity matrix -dmnsn_matrix -dmnsn_identity_matrix(void) -{ - return dmnsn_new_matrix(1.0, 0.0, 0.0, 0.0, - 0.0, 1.0, 0.0, 0.0, - 0.0, 0.0, 1.0, 0.0); -} - -// Scaling matrix -dmnsn_matrix -dmnsn_scale_matrix(dmnsn_vector s) -{ - return dmnsn_new_matrix(s.x, 0.0, 0.0, 0.0, - 0.0, s.y, 0.0, 0.0, - 0.0, 0.0, s.z, 0.0); -} - -// Translation matrix -dmnsn_matrix -dmnsn_translation_matrix(dmnsn_vector d) -{ - return dmnsn_new_matrix(1.0, 0.0, 0.0, d.x, - 0.0, 1.0, 0.0, d.y, - 0.0, 0.0, 1.0, d.z); -} - // Left-handed rotation matrix; theta/|theta| = axis, |theta| = angle dmnsn_matrix dmnsn_rotation_matrix(dmnsn_vector theta) { - // Two trig calls, 25 multiplications, 13 additions - double angle = dmnsn_vector_norm(theta); if (fabs(angle) < dmnsn_epsilon) { return dmnsn_identity_matrix(); } - dmnsn_vector axis = dmnsn_vector_div(theta, angle); + dmnsn_vector axis = dmnsn_vector_mul(1.0/angle, theta); // Shorthand to make dmnsn_new_matrix() call legible - double s = sin(angle); double t = 1.0 - cos(angle); - - double x = axis.x; - double y = axis.y; - double z = axis.z; + double x = axis.X; + double y = axis.Y; + double z = axis.Z; return dmnsn_new_matrix( 1.0 + t*(x*x - 1.0), -z*s + t*x*y, y*s + t*x*z, 0.0, @@ -87,7 +56,7 @@ static double dmnsn_axis_angle(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis) { from = dmnsn_vector_sub(from, dmnsn_vector_proj(from, axis)); - to = dmnsn_vector_sub(to, dmnsn_vector_proj(to, axis)); + to = dmnsn_vector_sub(to, dmnsn_vector_proj(to, axis)); double fromnorm = dmnsn_vector_norm(from); double tonorm = dmnsn_vector_norm(to); @@ -95,8 +64,8 @@ dmnsn_axis_angle(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis) return 0.0; } - from = dmnsn_vector_div(from, fromnorm); - to = dmnsn_vector_div(to, tonorm); + from = dmnsn_vector_mul(1.0/fromnorm, from); + to = dmnsn_vector_mul(1.0/tonorm, to); double angle = acos(dmnsn_vector_dot(from, to)); @@ -109,8 +78,7 @@ dmnsn_axis_angle(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis) // Alignment matrix dmnsn_matrix -dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to, - dmnsn_vector axis1, dmnsn_vector axis2) +dmnsn_alignment_matrix(dmnsn_vector from, dmnsn_vector to, dmnsn_vector axis1, dmnsn_vector axis2) { double theta1 = dmnsn_axis_angle(from, to, axis1); dmnsn_matrix align1 = dmnsn_rotation_matrix(dmnsn_vector_mul(theta1, axis1)); @@ -141,13 +109,13 @@ static dmnsn_matrix2 dmnsn_matrix2_sub(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs); static dmnsn_matrix2 dmnsn_matrix2_mul(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs); /// Invert a matrix with the slower cofactor algorithm, if partitioning failed. -static dmnsn_matrix dmnsn_matrix_inverse_generic(dmnsn_matrix A); +static dmnsn_matrix dmnsn_matrix_inverse_generic(const dmnsn_matrix M); /// Get the [\p row, \p col] cofactor of A. -static double dmnsn_matrix_cofactor(dmnsn_matrix A, size_t row, size_t col); +static double dmnsn_matrix_cofactor(dmnsn_matrix M, size_t row, size_t col); // Invert a matrix, by partitioning dmnsn_matrix -dmnsn_matrix_inverse(dmnsn_matrix A) +dmnsn_matrix_inverse(dmnsn_matrix M) { // Use partitioning to invert a matrix: // @@ -162,11 +130,8 @@ dmnsn_matrix_inverse(dmnsn_matrix A) // RR = -SS*R*inv(P), and // SS = inv(S - R*inv(P)*Q). - // The algorithm uses 2 inversions, 6 multiplications, and 2 subtractions, - // giving 52 multiplications, 34 additions, and 8 divisions. - dmnsn_matrix2 P, Q, R, S, Pi, RPi, PiQ, RPiQ, PP, QQ, RR, SS; - double Pdet = A.n[0][0]*A.n[1][1] - A.n[0][1]*A.n[1][0]; + double Pdet = M.n[0][0]*M.n[1][1] - M.n[0][1]*M.n[1][0]; if (dmnsn_unlikely(fabs(Pdet) < dmnsn_epsilon)) { // If P is close to singular, try a more generic algorithm; this is very @@ -175,22 +140,24 @@ dmnsn_matrix_inverse(dmnsn_matrix A) // [ 1 1 1 0 ] // [ 0 1 1 0 ] // [ 0 0 0 1 ] - return dmnsn_matrix_inverse_generic(A); + return dmnsn_matrix_inverse_generic(M); } + double Pdet_inv = 1.0/Pdet; + // Partition the matrix - P = dmnsn_new_matrix2(A.n[0][0], A.n[0][1], - A.n[1][0], A.n[1][1]); - Q = dmnsn_new_matrix2(A.n[0][2], A.n[0][3], - A.n[1][2], A.n[1][3]); - R = dmnsn_new_matrix2(A.n[2][0], A.n[2][1], + P = dmnsn_new_matrix2(M.n[0][0], M.n[0][1], + M.n[1][0], M.n[1][1]); + Q = dmnsn_new_matrix2(M.n[0][2], M.n[0][3], + M.n[1][2], M.n[1][3]); + R = dmnsn_new_matrix2(M.n[2][0], M.n[2][1], 0.0, 0.0); - S = dmnsn_new_matrix2(A.n[2][2], A.n[2][3], + S = dmnsn_new_matrix2(M.n[2][2], M.n[2][3], 0.0, 1.0); // Do this inversion ourselves, since we already have the determinant - Pi = dmnsn_new_matrix2( P.n[1][1]/Pdet, -P.n[0][1]/Pdet, - -P.n[1][0]/Pdet, P.n[0][0]/Pdet); + Pi = dmnsn_new_matrix2( P.n[1][1]*Pdet_inv, -P.n[0][1]*Pdet_inv, + -P.n[1][0]*Pdet_inv, P.n[0][0]*Pdet_inv); // Calculate R*inv(P), inv(P)*Q, and R*inv(P)*Q RPi = dmnsn_matrix2_mul(R, Pi); @@ -204,9 +171,11 @@ dmnsn_matrix_inverse(dmnsn_matrix A) PP = dmnsn_matrix2_sub(Pi, dmnsn_matrix2_mul(PiQ, RR)); // Reconstruct the matrix - return dmnsn_new_matrix(PP.n[0][0], PP.n[0][1], QQ.n[0][0], QQ.n[0][1], - PP.n[1][0], PP.n[1][1], QQ.n[1][0], QQ.n[1][1], - RR.n[0][0], RR.n[0][1], SS.n[0][0], SS.n[0][1]); + return dmnsn_new_matrix( + PP.n[0][0], PP.n[0][1], QQ.n[0][0], QQ.n[0][1], + PP.n[1][0], PP.n[1][1], QQ.n[1][0], QQ.n[1][1], + RR.n[0][0], RR.n[0][1], SS.n[0][0], SS.n[0][1] + ); } // For nice shorthand @@ -222,10 +191,9 @@ dmnsn_new_matrix2(double a1, double a2, double b1, double b2) static dmnsn_matrix2 dmnsn_matrix2_inverse(dmnsn_matrix2 A) { - // 4 divisions, 2 multiplications, 1 addition - double det = A.n[0][0]*A.n[1][1] - A.n[0][1]*A.n[1][0]; - return dmnsn_new_matrix2( A.n[1][1]/det, -A.n[0][1]/det, - -A.n[1][0]/det, A.n[0][0]/det); + double inv_det = 1.0/(A.n[0][0]*A.n[1][1] - A.n[0][1]*A.n[1][0]); + return dmnsn_new_matrix2( A.n[1][1]*inv_det, -A.n[0][1]*inv_det, + -A.n[1][0]*inv_det, A.n[0][0]*inv_det); } // Also basically a shorthand @@ -240,7 +208,6 @@ dmnsn_matrix2_negate(dmnsn_matrix2 A) static dmnsn_matrix2 dmnsn_matrix2_sub(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs) { - // 4 additions return dmnsn_new_matrix2( lhs.n[0][0] - rhs.n[0][0], lhs.n[0][1] - rhs.n[0][1], lhs.n[1][0] - rhs.n[1][0], lhs.n[1][1] - rhs.n[1][1] @@ -251,7 +218,6 @@ dmnsn_matrix2_sub(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs) static dmnsn_matrix2 dmnsn_matrix2_mul(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs) { - // 8 multiplications, 4 additions return dmnsn_new_matrix2( lhs.n[0][0]*rhs.n[0][0] + lhs.n[0][1]*rhs.n[1][0], lhs.n[0][0]*rhs.n[0][1] + lhs.n[0][1]*rhs.n[1][1], @@ -262,7 +228,7 @@ dmnsn_matrix2_mul(dmnsn_matrix2 lhs, dmnsn_matrix2 rhs) // Invert a matrix, if partitioning failed (|P| == 0) static dmnsn_matrix -dmnsn_matrix_inverse_generic(dmnsn_matrix A) +dmnsn_matrix_inverse_generic(dmnsn_matrix M) { // For A = [ A' b ] A^-1 = [ A'^-1 -(A'^-1)*b ] // [ 0 ... 0 1 ], [ 0 ... 0 1 ]. @@ -274,20 +240,22 @@ dmnsn_matrix_inverse_generic(dmnsn_matrix A) // Perform a Laplace expansion along the first row to give us the adjugate's // first column and the determinant for (size_t j = 0; j < 3; ++j) { - C = dmnsn_matrix_cofactor(A, 0, j); - det += A.n[0][j]*C; + C = dmnsn_matrix_cofactor(M, 0, j); + det += M.n[0][j]*C; inv.n[j][0] = C; } + double inv_det = 1.0/det; + // Divide the first column by the determinant for (size_t j = 0; j < 3; ++j) { - inv.n[j][0] /= det; + inv.n[j][0] *= inv_det; } // Find the rest of A' for (size_t j = 0; j < 3; ++j) { for (size_t i = 1; i < 3; ++i) { - inv.n[j][i] = dmnsn_matrix_cofactor(A, i, j)/det; + inv.n[j][i] = dmnsn_matrix_cofactor(M, i, j)*inv_det; } inv.n[j][3] = 0.0; } @@ -295,7 +263,7 @@ dmnsn_matrix_inverse_generic(dmnsn_matrix A) // Find the translational component of the inverse for (size_t i = 0; i < 3; ++i) { for (size_t j = 0; j < 3; ++j) { - inv.n[i][3] -= inv.n[i][j]*A.n[j][3]; + inv.n[i][3] -= inv.n[i][j]*M.n[j][3]; } } @@ -306,15 +274,14 @@ dmnsn_matrix_inverse_generic(dmnsn_matrix A) // upper-left 3x3 corner of A by ignoring row `row' and column `col', // times (-1)^(row + col) static double -dmnsn_matrix_cofactor(dmnsn_matrix A, size_t row, size_t col) +dmnsn_matrix_cofactor(dmnsn_matrix M, size_t row, size_t col) { - // 2 multiplications, 1 addition double n[4]; size_t k = 0; for (size_t i = 0; i < 3; ++i) { for (size_t j = 0; j < 3; ++j) { if (i != row && j != col) { - n[k] = A.n[i][j]; + n[k] = M.n[i][j]; ++k; } } @@ -328,61 +295,30 @@ dmnsn_matrix_cofactor(dmnsn_matrix A, size_t row, size_t col) } } -// 4x4 matrix multiplication -dmnsn_matrix -dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs) -{ - // 36 multiplications, 27 additions - dmnsn_matrix r; - - r.n[0][0] = lhs.n[0][0]*rhs.n[0][0] + lhs.n[0][1]*rhs.n[1][0] + lhs.n[0][2]*rhs.n[2][0]; - r.n[0][1] = lhs.n[0][0]*rhs.n[0][1] + lhs.n[0][1]*rhs.n[1][1] + lhs.n[0][2]*rhs.n[2][1]; - r.n[0][2] = lhs.n[0][0]*rhs.n[0][2] + lhs.n[0][1]*rhs.n[1][2] + lhs.n[0][2]*rhs.n[2][2]; - r.n[0][3] = lhs.n[0][0]*rhs.n[0][3] + lhs.n[0][1]*rhs.n[1][3] + lhs.n[0][2]*rhs.n[2][3] + lhs.n[0][3]; - - r.n[1][0] = lhs.n[1][0]*rhs.n[0][0] + lhs.n[1][1]*rhs.n[1][0] + lhs.n[1][2]*rhs.n[2][0]; - r.n[1][1] = lhs.n[1][0]*rhs.n[0][1] + lhs.n[1][1]*rhs.n[1][1] + lhs.n[1][2]*rhs.n[2][1]; - r.n[1][2] = lhs.n[1][0]*rhs.n[0][2] + lhs.n[1][1]*rhs.n[1][2] + lhs.n[1][2]*rhs.n[2][2]; - r.n[1][3] = lhs.n[1][0]*rhs.n[0][3] + lhs.n[1][1]*rhs.n[1][3] + lhs.n[1][2]*rhs.n[2][3] + lhs.n[1][3]; - - r.n[2][0] = lhs.n[2][0]*rhs.n[0][0] + lhs.n[2][1]*rhs.n[1][0] + lhs.n[2][2]*rhs.n[2][0]; - r.n[2][1] = lhs.n[2][0]*rhs.n[0][1] + lhs.n[2][1]*rhs.n[1][1] + lhs.n[2][2]*rhs.n[2][1]; - r.n[2][2] = lhs.n[2][0]*rhs.n[0][2] + lhs.n[2][1]*rhs.n[1][2] + lhs.n[2][2]*rhs.n[2][2]; - r.n[2][3] = lhs.n[2][0]*rhs.n[0][3] + lhs.n[2][1]*rhs.n[1][3] + lhs.n[2][2]*rhs.n[2][3] + lhs.n[2][3]; - - return r; -} - // Give an axis-aligned box that contains the given box transformed by `lhs' dmnsn_aabb -dmnsn_transform_aabb(dmnsn_matrix trans, dmnsn_aabb box) +dmnsn_transform_aabb(dmnsn_matrix M, dmnsn_aabb box) { // Infinite/zero bounding box support - if (isinf(box.min.x)) { + if (isinf(box.min.X)) { return box; } // Taking the "absolute value" of the matrix saves some min/max calculations - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - trans.n[i][j] = fabs(trans.n[i][j]); + dmnsn_vector Mabs[3]; + for (unsigned int i = 0; i < 3; ++i) { + for (unsigned int j = 0; j < 3; ++j) { + Mabs[i].n[j] = fabs(M.n[j][i]); } } - dmnsn_vector Mt = dmnsn_matrix_column(trans, 3); + dmnsn_vector Mt = dmnsn_matrix_column(M, 3); dmnsn_aabb ret = { Mt, Mt }; - dmnsn_vector Mz = dmnsn_matrix_column(trans, 2); - ret.min = dmnsn_vector_add(ret.min, dmnsn_vector_mul(box.min.z, Mz)); - ret.max = dmnsn_vector_add(ret.max, dmnsn_vector_mul(box.max.z, Mz)); - - dmnsn_vector My = dmnsn_matrix_column(trans, 1); - ret.min = dmnsn_vector_add(ret.min, dmnsn_vector_mul(box.min.y, My)); - ret.max = dmnsn_vector_add(ret.max, dmnsn_vector_mul(box.max.y, My)); - - dmnsn_vector Mx = dmnsn_matrix_column(trans, 0); - ret.min = dmnsn_vector_add(ret.min, dmnsn_vector_mul(box.min.x, Mx)); - ret.max = dmnsn_vector_add(ret.max, dmnsn_vector_mul(box.max.x, Mx)); + for (unsigned int i = 0; i < 3; ++i) { + ret.min = dmnsn_vector_add(ret.min, dmnsn_vector_mul(box.min.n[i], Mabs[i])); + ret.max = dmnsn_vector_add(ret.max, dmnsn_vector_mul(box.max.n[i], Mabs[i])); + } return ret; } diff --git a/libdimension/model/objects/cone.c b/libdimension/model/objects/cone.c index 26e59ca..60317bb 100644 --- a/libdimension/model/objects/cone.c +++ b/libdimension/model/objects/cone.c @@ -43,11 +43,11 @@ dmnsn_cone_intersection_fn(const dmnsn_object *object, dmnsn_ray l, // 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; + 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); @@ -58,7 +58,7 @@ dmnsn_cone_intersection_fn(const dmnsn_object *object, dmnsn_ray l, t = dmnsn_min(t, x[1]); p = dmnsn_ray_point(l, t); - if (p.y <= -1.0 || p.y >= 1.0) { + if (p.Y <= -1.0 || p.Y >= 1.0) { t = dmnsn_max(x[0], x[1]); p = dmnsn_ray_point(l, t); } @@ -66,9 +66,9 @@ dmnsn_cone_intersection_fn(const dmnsn_object *object, dmnsn_ray l, 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); + 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; @@ -84,9 +84,9 @@ 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; + 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. @@ -118,12 +118,12 @@ static bool dmnsn_cone_cap_intersection_fn(const dmnsn_object *object, dmnsn_ray l, dmnsn_intersection *intersection) { - if (l.n.y != 0.0) { + 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; + 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) { + 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; @@ -188,12 +188,8 @@ dmnsn_new_cone(dmnsn_pool *pool, double r1, double r2, bool open) // 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) - ); + 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; diff --git a/libdimension/model/objects/cube.c b/libdimension/model/objects/cube.c index 7d6fe0f..0434d89 100644 --- a/libdimension/model/objects/cube.c +++ b/libdimension/model/objects/cube.c @@ -23,6 +23,7 @@ * Cubes. */ +#include "internal.h" #include "dimension/model.h" #include <math.h> @@ -36,8 +37,8 @@ dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray, 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; + 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; @@ -54,8 +55,8 @@ dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray, 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; + 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) { @@ -80,8 +81,8 @@ dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray, 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; + 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) { @@ -124,9 +125,9 @@ dmnsn_cube_intersection_fn(const dmnsn_object *cube, dmnsn_ray ray, 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; + 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. diff --git a/libdimension/model/objects/sphere.c b/libdimension/model/objects/sphere.c index e1ca784..2c8a8ef 100644 --- a/libdimension/model/objects/sphere.c +++ b/libdimension/model/objects/sphere.c @@ -29,8 +29,7 @@ /// Sphere intersection callback. static bool -dmnsn_sphere_intersection_fn(const dmnsn_object *sphere, dmnsn_ray l, - dmnsn_intersection *intersection) +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]; @@ -58,12 +57,12 @@ dmnsn_sphere_intersection_fn(const dmnsn_object *sphere, dmnsn_ray l, 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; + 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]) +dmnsn_implicit_dot(const double row[4]) { double ret = 0.0; for (int i = 0; i < 3; ++i) { @@ -84,18 +83,18 @@ dmnsn_sphere_bounding_fn(const dmnsn_object *object, dmnsn_matrix trans) 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; + 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; + 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; + box.min.Z = cz - dz; + box.max.Z = cz + dz; return box; } diff --git a/libdimension/model/objects/torus.c b/libdimension/model/objects/torus.c index b4baebd..e4894b3 100644 --- a/libdimension/model/objects/torus.c +++ b/libdimension/model/objects/torus.c @@ -41,18 +41,18 @@ dmnsn_torus_bound_intersection(const dmnsn_torus *torus, dmnsn_ray l) 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; + 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; + 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 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[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; @@ -63,7 +63,7 @@ dmnsn_torus_bound_intersection(const dmnsn_torus *torus, dmnsn_ray l) 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) + if (p.Y >= -r && p.Y <= r) break; } @@ -90,8 +90,8 @@ dmnsn_torus_intersection_fn(const dmnsn_object *object, dmnsn_ray l, } // 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); + 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); @@ -123,7 +123,7 @@ dmnsn_torus_intersection_fn(const dmnsn_object *object, dmnsn_ray l, 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_normalized(dmnsn_new_vector(p.X, 0.0, p.Z)) ); dmnsn_vector normal = dmnsn_vector_sub(p, center); @@ -137,8 +137,8 @@ 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; + 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. diff --git a/libdimension/model/objects/triangle.c b/libdimension/model/objects/triangle.c index 5af3301..87372d7 100644 --- a/libdimension/model/objects/triangle.c +++ b/libdimension/model/objects/triangle.c @@ -33,16 +33,15 @@ 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; + *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) +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)) { diff --git a/libdimension/model/objects/triangle_fan.c b/libdimension/model/objects/triangle_fan.c index 93768a9..ddca581 100644 --- a/libdimension/model/objects/triangle_fan.c +++ b/libdimension/model/objects/triangle_fan.c @@ -40,9 +40,9 @@ 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 + 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 ); } @@ -51,9 +51,9 @@ 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 + 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 ); } @@ -80,12 +80,11 @@ dmnsn_compress_coeffs(double coeffs[6], dmnsn_matrix incremental) static inline dmnsn_matrix dmnsn_decompress_coeffs(const double coeffs[6]) { - dmnsn_matrix incremental = dmnsn_new_matrix( + return 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. @@ -100,9 +99,9 @@ dmnsn_triangle_basis(dmnsn_vector a, dmnsn_vector ab, dmnsn_vector ac) 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; + *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; } @@ -149,7 +148,7 @@ dmnsn_triangle_fan_inside_fn(const dmnsn_object *object, dmnsn_vector point) return false; } -/// Computes the bounding box for the first triangle +/// Computes the bounding box for the first triangle. static inline dmnsn_aabb dmnsn_bound_first_triangle(dmnsn_matrix trans) { @@ -336,9 +335,9 @@ dmnsn_new_smooth_triangle_fan(dmnsn_pool *pool, dmnsn_vector vertices[], dmnsn_v 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; + coeffs[6] = nac.X; + coeffs[7] = nac.Y; + coeffs[8] = nac.Z; P = newP; } diff --git a/libdimension/model/pigments/canvas_pigment.c b/libdimension/model/pigments/canvas_pigment.c index bb83b0a..c44ffa6 100644 --- a/libdimension/model/pigments/canvas_pigment.c +++ b/libdimension/model/pigments/canvas_pigment.c @@ -23,6 +23,7 @@ * Image maps. */ +#include "internal.h" #include "dimension/model.h" /// Canvas pigment type. @@ -38,8 +39,8 @@ dmnsn_canvas_pigment_fn(const dmnsn_pigment *pigment, dmnsn_vector v) const dmnsn_canvas_pigment *canvas_pigment = (const dmnsn_canvas_pigment *)pigment; dmnsn_canvas *canvas = canvas_pigment->canvas; - size_t x = llround((fmod(v.x, 1.0) + 1.0)*(canvas->width - 1)); - size_t y = llround((fmod(v.y, 1.0) + 1.0)*(canvas->height - 1)); + size_t x = llround((fmod(v.X, 1.0) + 1.0)*(canvas->width - 1)); + size_t y = llround((fmod(v.Y, 1.0) + 1.0)*(canvas->height - 1)); return dmnsn_canvas_get_pixel(canvas, x%canvas->width, y%canvas->height); } diff --git a/libdimension/model/texture.c b/libdimension/model/texture.c index b7eb7ef..dc25ac9 100644 --- a/libdimension/model/texture.c +++ b/libdimension/model/texture.c @@ -45,8 +45,7 @@ dmnsn_texture_initialize(dmnsn_texture *texture) texture->trans_inv = dmnsn_matrix_inverse(texture->trans); if (!texture->pigment->initialized) { - texture->pigment->trans = dmnsn_matrix_mul(texture->trans, - texture->pigment->trans); + texture->pigment->trans = dmnsn_matrix_mul(texture->trans, texture->pigment->trans); dmnsn_pigment_initialize(texture->pigment); } } diff --git a/libdimension/pattern/checker.c b/libdimension/pattern/checker.c index cce9623..9f378ac 100644 --- a/libdimension/pattern/checker.c +++ b/libdimension/pattern/checker.c @@ -23,15 +23,16 @@ * Checker pattern. */ +#include "internal.h" #include "dimension/pattern.h" /// Checker pattern callback. static double dmnsn_checker_pattern_fn(const dmnsn_pattern *checker, dmnsn_vector v) { - double xmod = fmod(v.x, 2.0); - double ymod = fmod(v.y, 2.0); - double zmod = fmod(v.z, 2.0); + double xmod = fmod(v.X, 2.0); + double ymod = fmod(v.Y, 2.0); + double zmod = fmod(v.Z, 2.0); if (xmod < -dmnsn_epsilon) xmod += 2.0; diff --git a/libdimension/pattern/leopard.c b/libdimension/pattern/leopard.c index 1a7bce0..f17bb19 100644 --- a/libdimension/pattern/leopard.c +++ b/libdimension/pattern/leopard.c @@ -23,6 +23,7 @@ * Leopard pattern. */ +#include "internal.h" #include "dimension/pattern.h" #include <math.h> @@ -30,7 +31,7 @@ static double dmnsn_leopard_pattern_fn(const dmnsn_pattern *leopard, dmnsn_vector v) { - double val = (sin(v.x) + sin(v.y) + sin(v.z))/3.0; + double val = (sin(v.X) + sin(v.Y) + sin(v.Z))/3.0; return val*val; } diff --git a/libdimension/render/render.c b/libdimension/render/render.c index 842b41e..c0878a0 100644 --- a/libdimension/render/render.c +++ b/libdimension/render/render.c @@ -189,10 +189,7 @@ dmnsn_rtstate_initialize(dmnsn_rtstate *state, state->interior = intersection->object->interior; state->r = dmnsn_ray_point(intersection->ray, intersection->t); - state->pigment_r = dmnsn_transform_point( - intersection->object->pigment_trans, - state->r - ); + state->pigment_r = dmnsn_transform_point(intersection->object->pigment_trans, state->r); state->viewer = dmnsn_vector_normalized( dmnsn_vector_negate(intersection->ray.n) ); diff --git a/libdimension/tests/bvh/prtree.c b/libdimension/tests/bvh/prtree.c index fb0cea6..31f1490 100644 --- a/libdimension/tests/bvh/prtree.c +++ b/libdimension/tests/bvh/prtree.c @@ -30,13 +30,13 @@ #include <stdio.h> #include <stdlib.h> -unsigned int calls = 0; +static unsigned int calls = 0; static bool dmnsn_fake_intersection_fn(const dmnsn_object *object, dmnsn_ray ray, dmnsn_intersection *intersection) { - intersection->t = (object->aabb.min.z - ray.x0.z)/ray.n.z; + intersection->t = (object->aabb.min.Z - ray.x0.Z)/ray.n.Z; intersection->normal = dmnsn_x; ++calls; return true; @@ -47,13 +47,10 @@ dmnsn_randomize_aabb(dmnsn_object *object) { dmnsn_vector a, b; - a.x = 2.0*((double)rand())/RAND_MAX - 1.0; - a.y = 2.0*((double)rand())/RAND_MAX - 1.0; - a.z = 2.0*((double)rand())/RAND_MAX - 1.0; - - b.x = 2.0*((double)rand())/RAND_MAX - 1.0; - b.y = 2.0*((double)rand())/RAND_MAX - 1.0; - b.z = 2.0*((double)rand())/RAND_MAX - 1.0; + for (unsigned int i = 0; i < 3; ++i) { + a.n[i] = 2.0*((double)rand())/RAND_MAX - 1.0; + b.n[i] = 2.0*((double)rand())/RAND_MAX - 1.0; + } object->aabb.min = dmnsn_vector_min(a, b); object->aabb.max = dmnsn_vector_max(a, b); |