diff options
Diffstat (limited to 'libdimension/dimension/math')
| -rw-r--r-- | libdimension/dimension/math/aabb.h | 25 | ||||
| -rw-r--r-- | libdimension/dimension/math/matrix.h | 210 | ||||
| -rw-r--r-- | libdimension/dimension/math/vector.h | 104 |
3 files changed, 213 insertions, 126 deletions
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; } |