From a10d28667d734dfc27b0eaac5444ff944d046b94 Mon Sep 17 00:00:00 2001 From: 3gg <3gg@shellblade.net> Date: Sat, 4 Dec 2021 18:39:11 -0800 Subject: Initial commit. --- include/math/camera.h | 42 +++++ include/math/defs.h | 51 +++++ include/math/float.h | 9 + include/math/fwd.h | 9 + include/math/mat4.h | 488 ++++++++++++++++++++++++++++++++++++++++++++++++ include/math/spatial3.h | 194 +++++++++++++++++++ include/math/vec.h | 16 ++ include/math/vec2.h | 10 + include/math/vec3.h | 143 ++++++++++++++ include/math/vec4.h | 15 ++ 10 files changed, 977 insertions(+) create mode 100644 include/math/camera.h create mode 100644 include/math/defs.h create mode 100644 include/math/float.h create mode 100644 include/math/fwd.h create mode 100644 include/math/mat4.h create mode 100644 include/math/spatial3.h create mode 100644 include/math/vec.h create mode 100644 include/math/vec2.h create mode 100644 include/math/vec3.h create mode 100644 include/math/vec4.h (limited to 'include') diff --git a/include/math/camera.h b/include/math/camera.h new file mode 100644 index 0000000..1621927 --- /dev/null +++ b/include/math/camera.h @@ -0,0 +1,42 @@ +#pragma once + +#include "mat4.h" +#include "spatial3.h" + +typedef struct Camera { + Spatial3 spatial; + mat4 projection; +} Camera; + +/// Create an orthographic camera. +/// +/// The camera is positioned at the origin with canonical right/up/forward +/// vectors. +/// +/// \param left The coordinate for the left vertical clipping plane. +/// \param right The coordinate for the right vertical clipping plane. +/// \param bottom The coordinate for the bottom horizontal clipping plane. +/// \param top The coordinate for the top horizontal clipping plane. +/// \param near The distance to the near clipping plane. +/// \param far The distance to the far clipping plane. +static inline Camera camera_orthographic(R left, R right, R bottom, R top, + R near, R far) { + return (Camera){.spatial = spatial3_make(), + .projection = + mat4_ortho(left, right, bottom, top, near, far)}; +} + +/// Create a perspective camera. +/// +/// The camera is positioned at the origin with canonical right/up/forward +/// vectors. +/// +/// \param fovy The vertical field of view angle in degrees. +/// \param aspect The aspect ratio that determines the field of view in the +/// x-direction. +/// \param near Distance to the near clipping plane. +/// \param far Distance to the far clipping plane. +static inline Camera camera_perspective(R fovy, R aspect, R near, R far) { + return (Camera){.spatial = spatial3_make(), + .projection = mat4_perspective(fovy, aspect, near, far)}; +} diff --git a/include/math/defs.h b/include/math/defs.h new file mode 100644 index 0000000..acfb037 --- /dev/null +++ b/include/math/defs.h @@ -0,0 +1,51 @@ +#pragma once + +// +// Configuration macros: +// +// MATH_USE_FLOAT +// - Use floats instead of doubles for scalar values. + +#include +#include + +#ifdef MATH_USE_DOUBLE +typedef double R; +#define R_MIN DBL_MIN +#define R_MAX DBL_MAX +#else // floats +typedef float R; +#define R_MIN FLT_MIN +#define R_MAX FLT_MAX +#endif + +#define PI 3.14159265359 +#define INV_PI 0.31830988618 + +/// Radians per degree. +#define TO_RAD (PI / 180.0) + +/// Degrees per radian. +#define TO_DEG (180.0 / PI) + +#ifdef MATH_USE_DOUBLE +static inline double min(R a, R b) { return fmin(a, b); } +static inline double max(R a, R b) { return fmax(a, b); } +#else // floats +static inline float min(R a, R b) { return fminf(a, b); } +static inline float max(R a, R b) { return fmaxf(a, b); } +#endif + +static inline R rabs(R x) { return x >= 0.0 ? x : -x; } +static inline R clamp(R x, R low, R high) { return max(low, min(high, x)); } +static inline R sq(R x) { return x * x; } +static inline R sign(R x) { + if (x < 0) { + return -1; + } else if (x > 0) { + return 1; + } else { + return 0; + } +} +static inline R lerp(R a, R b, R t) { return a + (b - a) * t; } diff --git a/include/math/float.h b/include/math/float.h new file mode 100644 index 0000000..9d289b9 --- /dev/null +++ b/include/math/float.h @@ -0,0 +1,9 @@ +#pragma once + +#include +#include + +/// Compare two floats for equality. +static inline bool float_eq(float a, float b, float eps) { + return fabsf(a - b) <= eps; +} diff --git a/include/math/fwd.h b/include/math/fwd.h new file mode 100644 index 0000000..6c72ded --- /dev/null +++ b/include/math/fwd.h @@ -0,0 +1,9 @@ +/// Forward declarations for all math objects. +#pragma once + +typedef struct Camera Camera; +typedef struct mat4 mat4; +typedef struct spatial3 spatial3; +typedef struct vec2 vec2; +typedef struct vec3 vec3; +typedef struct vec4 vec4; diff --git a/include/math/mat4.h b/include/math/mat4.h new file mode 100644 index 0000000..e6c707a --- /dev/null +++ b/include/math/mat4.h @@ -0,0 +1,488 @@ +#pragma once + +#include "defs.h" +#include "float.h" +#include "vec3.h" +#include "vec4.h" + +#include + +/// A 4x4 column-major matrix. +typedef struct mat4 { + R val[4][4]; +} mat4; + +/// Construct a matrix from 16 values. +static inline mat4 mat4_make(R m00, R m10, R m20, R m30, R m01, R m11, R m21, + R m31, R m02, R m12, R m22, R m32, R m03, R m13, + R m23, R m33) { + mat4 m; + m.val[0][0] = m00; + m.val[0][1] = m01; + m.val[0][2] = m02; + m.val[0][3] = m03; + + m.val[1][0] = m10; + m.val[1][1] = m11; + m.val[1][2] = m12; + m.val[1][3] = m13; + + m.val[2][0] = m20; + m.val[2][1] = m21; + m.val[2][2] = m22; + m.val[2][3] = m23; + + m.val[3][0] = m30; + m.val[3][1] = m31; + m.val[3][2] = m32; + m.val[3][3] = m33; + return m; +} + +/// Construct a matrix from a column-major matrix array. +static inline mat4 mat4_from_array(const R M[16]) { + mat4 m; + m.val[0][0] = M[0]; + m.val[0][1] = M[1]; + m.val[0][2] = M[2]; + m.val[0][3] = M[3]; + + m.val[1][0] = M[4]; + m.val[1][1] = M[5]; + m.val[1][2] = M[6]; + m.val[1][3] = M[7]; + + m.val[2][0] = M[8]; + m.val[2][1] = M[9]; + m.val[2][2] = M[10]; + m.val[2][3] = M[11]; + + m.val[3][0] = M[12]; + m.val[3][1] = M[13]; + m.val[3][2] = M[14]; + m.val[3][3] = M[15]; + return m; +} + +/// Construct the identity matrix. +static inline mat4 mat4_id() { + return mat4_make(1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 1.0); +} + +/// Construct a matrix from 4 column vectors. +static inline mat4 mat4_from_vec4(vec4 v0, vec4 v1, vec4 v2, vec4 v3) { + return mat4_make(v0.x, v0.y, v0.z, v0.w, v1.x, v1.y, v1.z, v1.w, v2.x, v2.y, + v2.z, v2.w, v3.x, v3.y, v3.z, v3.w); +} + +/// Construct a transformation matrix from 4 vectors. +static inline mat4 mat4_from_vec3(vec3 right, vec3 up, vec3 forward, + vec3 position) { + return mat4_make(right.x, right.y, right.z, 0.0, up.x, up.y, up.z, 0.0, + forward.x, forward.y, forward.z, 0.0, position.x, position.y, + position.z, 1.0); +} + +/// Return the value at the specified position. +static inline R mat4_at(mat4 m, int row, int col) { return m.val[col][row]; } + +/// Return the matrix's first column. +static inline vec3 mat4_v0(mat4 m) { return *((vec3*)m.val[0]); } + +/// Return the matrix's second column. +static inline vec3 mat4_v1(mat4 m) { return *((vec3*)m.val[1]); } + +/// Return the matrix's third column. +static inline vec3 mat4_v2(mat4 m) { return *((vec3*)m.val[2]); } + +/// Return the matrix's fourth column. +static inline vec3 mat4_v3(mat4 m) { return *((vec3*)m.val[3]); } + +/// Set the matrix's first column. +static inline void mat4_set_v0(mat4* m, vec3 v) { *((vec3*)m->val[0]) = v; } + +/// Set the matrix's second column. +static inline void mat4_set_v1(mat4* m, vec3 v) { *((vec3*)m->val[1]) = v; } + +/// Set the matrix's third column. +static inline void mat4_set_v2(mat4* m, vec3 v) { *((vec3*)m->val[2]) = v; } + +/// Set the matrix's fourth column. +static inline void mat4_set_v3(mat4* m, vec3 v) { *((vec3*)m->val[3]) = v; } + +/// Set the 3x3 portion of the first matrix equal to the 3x3 portion of the +/// second matrix. +static inline void mat4_set_3x3(mat4* m, mat4 n) { + m->val[0][0] = n.val[0][0]; + m->val[0][1] = n.val[0][1]; + m->val[0][2] = n.val[0][2]; + + m->val[1][0] = n.val[1][0]; + m->val[1][1] = n.val[1][1]; + m->val[1][2] = n.val[1][2]; + + m->val[2][0] = n.val[2][0]; + m->val[2][1] = n.val[2][1]; + m->val[2][2] = n.val[2][2]; +} + +/// Multiply two matrices. +/// A * B = AB. +static inline mat4 mat4_mul(mat4 A, mat4 B) { + R m00 = mat4_at(A, 0, 0) * mat4_at(B, 0, 0) + + mat4_at(A, 0, 1) * mat4_at(B, 1, 0) + + mat4_at(A, 0, 2) * mat4_at(B, 2, 0) + + mat4_at(A, 0, 3) * mat4_at(B, 3, 0); + R m01 = mat4_at(A, 0, 0) * mat4_at(B, 0, 1) + + mat4_at(A, 0, 1) * mat4_at(B, 1, 1) + + mat4_at(A, 0, 2) * mat4_at(B, 2, 1) + + mat4_at(A, 0, 3) * mat4_at(B, 3, 1); + R m02 = mat4_at(A, 0, 0) * mat4_at(B, 0, 2) + + mat4_at(A, 0, 1) * mat4_at(B, 1, 2) + + mat4_at(A, 0, 2) * mat4_at(B, 2, 2) + + mat4_at(A, 0, 3) * mat4_at(B, 3, 2); + R m03 = mat4_at(A, 0, 0) * mat4_at(B, 0, 3) + + mat4_at(A, 0, 1) * mat4_at(B, 1, 3) + + mat4_at(A, 0, 2) * mat4_at(B, 2, 3) + + mat4_at(A, 0, 3) * mat4_at(B, 3, 3); + + R m10 = mat4_at(A, 1, 0) * mat4_at(B, 0, 0) + + mat4_at(A, 1, 1) * mat4_at(B, 1, 0) + + mat4_at(A, 1, 2) * mat4_at(B, 2, 0) + + mat4_at(A, 1, 3) * mat4_at(B, 3, 0); + R m11 = mat4_at(A, 1, 0) * mat4_at(B, 0, 1) + + mat4_at(A, 1, 1) * mat4_at(B, 1, 1) + + mat4_at(A, 1, 2) * mat4_at(B, 2, 1) + + mat4_at(A, 1, 3) * mat4_at(B, 3, 1); + R m12 = mat4_at(A, 1, 0) * mat4_at(B, 0, 2) + + mat4_at(A, 1, 1) * mat4_at(B, 1, 2) + + mat4_at(A, 1, 2) * mat4_at(B, 2, 2) + + mat4_at(A, 1, 3) * mat4_at(B, 3, 2); + R m13 = mat4_at(A, 1, 0) * mat4_at(B, 0, 3) + + mat4_at(A, 1, 1) * mat4_at(B, 1, 3) + + mat4_at(A, 1, 2) * mat4_at(B, 2, 3) + + mat4_at(A, 1, 3) * mat4_at(B, 3, 3); + + R m20 = mat4_at(A, 2, 0) * mat4_at(B, 0, 0) + + mat4_at(A, 2, 1) * mat4_at(B, 1, 0) + + mat4_at(A, 2, 2) * mat4_at(B, 2, 0) + + mat4_at(A, 2, 3) * mat4_at(B, 3, 0); + R m21 = mat4_at(A, 2, 0) * mat4_at(B, 0, 1) + + mat4_at(A, 2, 1) * mat4_at(B, 1, 1) + + mat4_at(A, 2, 2) * mat4_at(B, 2, 1) + + mat4_at(A, 2, 3) * mat4_at(B, 3, 1); + R m22 = mat4_at(A, 2, 0) * mat4_at(B, 0, 2) + + mat4_at(A, 2, 1) * mat4_at(B, 1, 2) + + mat4_at(A, 2, 2) * mat4_at(B, 2, 2) + + mat4_at(A, 2, 3) * mat4_at(B, 3, 2); + R m23 = mat4_at(A, 2, 0) * mat4_at(B, 0, 3) + + mat4_at(A, 2, 1) * mat4_at(B, 1, 3) + + mat4_at(A, 2, 2) * mat4_at(B, 2, 3) + + mat4_at(A, 2, 3) * mat4_at(B, 3, 3); + + R m30 = mat4_at(A, 3, 0) * mat4_at(B, 0, 0) + + mat4_at(A, 3, 1) * mat4_at(B, 1, 0) + + mat4_at(A, 3, 2) * mat4_at(B, 2, 0) + + mat4_at(A, 3, 3) * mat4_at(B, 3, 0); + R m31 = mat4_at(A, 3, 0) * mat4_at(B, 0, 1) + + mat4_at(A, 3, 1) * mat4_at(B, 1, 1) + + mat4_at(A, 3, 2) * mat4_at(B, 2, 1) + + mat4_at(A, 3, 3) * mat4_at(B, 3, 1); + R m32 = mat4_at(A, 3, 0) * mat4_at(B, 0, 2) + + mat4_at(A, 3, 1) * mat4_at(B, 1, 2) + + mat4_at(A, 3, 2) * mat4_at(B, 2, 2) + + mat4_at(A, 3, 3) * mat4_at(B, 3, 2); + R m33 = mat4_at(A, 3, 0) * mat4_at(B, 0, 3) + + mat4_at(A, 3, 1) * mat4_at(B, 1, 3) + + mat4_at(A, 3, 2) * mat4_at(B, 2, 3) + + mat4_at(A, 3, 3) * mat4_at(B, 3, 3); + + return mat4_make(m00, m01, m02, m03, m10, m11, m12, m13, m20, m21, m22, m23, + m30, m31, m32, m33); +} + +/// Return the translation component of the matrix. +static inline mat4 mat4_translation(mat4 m) { + return mat4_make(1.0, 0.0, 0.0, mat4_at(m, 0, 3), 0.0, 1.0, 0.0, + mat4_at(m, 1, 3), 0.0, 0.0, 1.0, mat4_at(m, 2, 3), 0.0, 0.0, + 0.0, 1.0); +} + +/// Return the rotation component of the matrix. +static inline mat4 mat4_rotation(mat4 m) { + return mat4_make(mat4_at(m, 0, 0), mat4_at(m, 0, 1), mat4_at(m, 0, 2), 0.0, + mat4_at(m, 1, 0), mat4_at(m, 1, 1), mat4_at(m, 1, 2), 0.0, + mat4_at(m, 2, 0), mat4_at(m, 2, 1), mat4_at(m, 2, 2), 0.0, + 0.0, 0.0, 0.0, 1.0); +} + +/// Create an X-axis rotation matrix. +static inline mat4 mat4_rotx(R angle) { + const R s = sin(angle); + const R c = cos(angle); + return mat4_make(1, 0, 0, 0, 0, c, -s, 0, 0, s, c, 0, 0, 0, 0, 1); +} + +/// Create a Y-axis rotation matrix. +static inline mat4 mat4_roty(R angle) { + const R s = sin(angle); + const R c = cos(angle); + return mat4_make(c, 0, s, 0, 0, 1, 0, 0, -s, 0, c, 0, 0, 0, 0, 1); +} + +/// Create a Z-axis rotation matrix. +static inline mat4 mat4_rotz(R angle) { + const R s = sin(angle); + const R c = cos(angle); + return mat4_make(c, -s, 0, 0, s, c, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1); +} + +/// Create a rotation matrix. +static inline mat4 mat4_rot(vec3 axis, R angle) { + const R s = sin(angle); + const R c = cos(angle); + const R x = axis.x; + const R y = axis.y; + const R z = axis.z; + const R xy = x * y; + const R xz = x * z; + const R yz = y * z; + const R sx = s * x; + const R sy = s * y; + const R sz = s * z; + const R omc = 1.0 - c; + return mat4_make(c + omc * x * x, omc * xy - sz, omc * xz + sy, 0, + omc * xy + sz, c + omc * y * y, omc * yz - sx, 0, + omc * xz - sy, omc * yz + sx, c + omc * z * z, 0, 0, 0, 0, + 1); +} + +/// Create a scaling matrix. +static inline mat4 mat4_scale(vec3 s) { + return mat4_make(s.x, 0, 0, 0, 0, s.y, 0, 0, 0, 0, s.z, 0, 0, 0, 0, 1); +} + +/// Create a translation matrix. +static inline mat4 mat4_translate(vec3 v) { + return mat4_make(1, 0, 0, v.x, 0, 1, 0, v.y, 0, 0, 1, v.z, 0, 0, 0, 1); +} + +/// The X-axis reflection matrix. +static inline mat4 mat4_reflectx() { + return mat4_make(-1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1); +} + +/// The Y-axis reflection matrix. +static inline mat4 mat4_reflecty() { + return mat4_make(1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1); +} + +/// The Z-axis reflection matrix. +static inline mat4 mat4_reflectz() { + return mat4_make(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1); +} + +/// Create a transformation matrix from the given forward vector. +static inline mat4 mat4_from_forward(vec3 forward) { + const vec3 f = vec3_normalize(forward); + const vec3 r = vec3_normalize(vec3_cross(f, up3())); + const vec3 u = vec3_normalize(vec3_cross(r, f)); + return mat4_make(r.x, u.x, -f.x, 0.0, r.y, u.y, -f.y, 0.0, r.z, u.z, -f.z, + 0.0, 0.0, 0.0, 0.0, 1.0); +} + +/// Create a transformation matrix. +static inline mat4 mat4_lookat(vec3 position, vec3 target, vec3 up) { + const vec3 fwd = vec3_normalize(vec3_sub(target, position)); + const vec3 right = vec3_normalize(vec3_cross(fwd, up)); + up = vec3_normalize(vec3_cross(right, fwd)); + return mat4_from_vec3(right, up, fwd, position); +} + +/// Create an orthographic projection matrix. +/// \param left The coordinate for the left vertical clipping plane. +/// \param right The coordinate for the right vertical clipping plane. +/// \param bottom The coordinate for the bottom horizontal clipping plane. +/// \param top The coordinate for the top horizontal clipping plane. +/// \param near The distance to the near clipping plane. +/// \param far The distance to the far clipping plane. +static inline mat4 mat4_ortho(R left, R right, R bottom, R top, R near, R far) { + const R tx = -(right + left) / (right - left); + const R ty = -(top + bottom) / (top - bottom); + const R tz = -(far + near) / (far - near); + return mat4_make(2 / (right - left), 0, 0, tx, 0, 2 / (top - bottom), 0, ty, + 0, 0, -2 / (far - near), tz, 0, 0, 0, 1); +} + +/// Create a perspective projection matrix. +/// \param fovy The vertical field of view angle in degrees. +/// \param aspect The aspect ratio that determines the field of view in the +/// x-direction. +/// \param near Distance to the near clipping plane. +/// \param far Distance to the far clipping plane. +static inline mat4 mat4_perspective(R fovy, R aspect, R near, R far) { + R f = tan(fovy / 2.0); + assert(f > 0.0); + f = 1.0 / f; + const R a = near - far; + return mat4_make(f / aspect, 0, 0, 0, 0, f, 0, 0, 0, 0, (far + near) / a, + (2 * far * near / a), 0, 0, -1, 0); +} + +/// Create the inverse of a perspective projection matrix. +/// \param fovy The vertical field of view angle in degrees. +/// \param aspect The aspect ratio that determines the field of view in the +/// x-direction. +/// \param near Distance to the near clipping plane. +/// \param far Distance to the far clipping plane. +static inline mat4 mat4_perspective_inverse(R fovy, R aspect, R near, R far) { + R f = tan(fovy / 2.0); + assert(f > 0.0); + f = 1.0 / f; + const R a = far * near; + const R P32 = 0.5 * (near - far) / a; + const R P33 = 0.5 * (far + near) / a; + return mat4_make(aspect / f, 0, 0, 0, 0, 1.0f / f, 0, 0, 0, 0, 0, -1, 0, 0, + P32, P33); +} + +/// Return the matrix's determinant. +static inline R mat4_det(mat4 m) { + const R* M = (const R*)(m.val); + const R inv0 = M[5] * M[10] * M[15] - M[5] * M[11] * M[14] - + M[9] * M[6] * M[15] + M[9] * M[7] * M[14] + + M[13] * M[6] * M[11] - M[13] * M[7] * M[10]; + const R inv1 = -M[4] * M[10] * M[15] + M[4] * M[11] * M[14] + + M[8] * M[6] * M[15] - M[8] * M[7] * M[14] - + M[12] * M[6] * M[11] + M[12] * M[7] * M[10]; + const R inv2 = M[4] * M[9] * M[15] - M[4] * M[11] * M[13] - + M[8] * M[5] * M[15] + M[8] * M[7] * M[13] + + M[12] * M[5] * M[11] - M[12] * M[7] * M[9]; + const R inv3 = -M[4] * M[9] * M[14] + M[4] * M[10] * M[13] + + M[8] * M[5] * M[14] - M[8] * M[6] * M[13] - + M[12] * M[5] * M[10] + M[12] * M[6] * M[9]; + return M[0] * inv0 + M[1] * inv1 + M[2] * inv2 + M[3] * inv3; +} + +/// Invert the matrix. +static inline mat4 mat4_inverse(mat4 m) { + const R* M = (const R*)(m.val); + + R inv[16]; + inv[0] = M[5] * M[10] * M[15] - M[5] * M[11] * M[14] - M[9] * M[6] * M[15] + + M[9] * M[7] * M[14] + M[13] * M[6] * M[11] - M[13] * M[7] * M[10]; + inv[4] = -M[4] * M[10] * M[15] + M[4] * M[11] * M[14] + M[8] * M[6] * M[15] - + M[8] * M[7] * M[14] - M[12] * M[6] * M[11] + M[12] * M[7] * M[10]; + inv[8] = M[4] * M[9] * M[15] - M[4] * M[11] * M[13] - M[8] * M[5] * M[15] + + M[8] * M[7] * M[13] + M[12] * M[5] * M[11] - M[12] * M[7] * M[9]; + inv[12] = -M[4] * M[9] * M[14] + M[4] * M[10] * M[13] + M[8] * M[5] * M[14] - + M[8] * M[6] * M[13] - M[12] * M[5] * M[10] + M[12] * M[6] * M[9]; + inv[1] = -M[1] * M[10] * M[15] + M[1] * M[11] * M[14] + M[9] * M[2] * M[15] - + M[9] * M[3] * M[14] - M[13] * M[2] * M[11] + M[13] * M[3] * M[10]; + inv[5] = M[0] * M[10] * M[15] - M[0] * M[11] * M[14] - M[8] * M[2] * M[15] + + M[8] * M[3] * M[14] + M[12] * M[2] * M[11] - M[12] * M[3] * M[10]; + inv[9] = -M[0] * M[9] * M[15] + M[0] * M[11] * M[13] + M[8] * M[1] * M[15] - + M[8] * M[3] * M[13] - M[12] * M[1] * M[11] + M[12] * M[3] * M[9]; + inv[13] = M[0] * M[9] * M[14] - M[0] * M[10] * M[13] - M[8] * M[1] * M[14] + + M[8] * M[2] * M[13] + M[12] * M[1] * M[10] - M[12] * M[2] * M[9]; + inv[2] = M[1] * M[6] * M[15] - M[1] * M[7] * M[14] - M[5] * M[2] * M[15] + + M[5] * M[3] * M[14] + M[13] * M[2] * M[7] - M[13] * M[3] * M[6]; + inv[6] = -M[0] * M[6] * M[15] + M[0] * M[7] * M[14] + M[4] * M[2] * M[15] - + M[4] * M[3] * M[14] - M[12] * M[2] * M[7] + M[12] * M[3] * M[6]; + inv[10] = M[0] * M[5] * M[15] - M[0] * M[7] * M[13] - M[4] * M[1] * M[15] + + M[4] * M[3] * M[13] + M[12] * M[1] * M[7] - M[12] * M[3] * M[5]; + inv[14] = -M[0] * M[5] * M[14] + M[0] * M[6] * M[13] + M[4] * M[1] * M[14] - + M[4] * M[2] * M[13] - M[12] * M[1] * M[6] + M[12] * M[2] * M[5]; + inv[3] = -M[1] * M[6] * M[11] + M[1] * M[7] * M[10] + M[5] * M[2] * M[11] - + M[5] * M[3] * M[10] - M[9] * M[2] * M[7] + M[9] * M[3] * M[6]; + inv[7] = M[0] * M[6] * M[11] - M[0] * M[7] * M[10] - M[4] * M[2] * M[11] + + M[4] * M[3] * M[10] + M[8] * M[2] * M[7] - M[8] * M[3] * M[6]; + inv[11] = -M[0] * M[5] * M[11] + M[0] * M[7] * M[9] + M[4] * M[1] * M[11] - + M[4] * M[3] * M[9] - M[8] * M[1] * M[7] + M[8] * M[3] * M[5]; + inv[15] = M[0] * M[5] * M[10] - M[0] * M[6] * M[9] - M[4] * M[1] * M[10] + + M[4] * M[2] * M[9] + M[8] * M[1] * M[6] - M[8] * M[2] * M[5]; + + R det = M[0] * inv[0] + M[1] * inv[4] + M[2] * inv[8] + M[3] * inv[12]; + assert(det != 0.0); + det = 1.0 / det; + return mat4_make(inv[0] * det, inv[4] * det, inv[8] * det, inv[12] * det, + inv[1] * det, inv[5] * det, inv[9] * det, inv[13] * det, + inv[2] * det, inv[6] * det, inv[10] * det, inv[14] * det, + inv[3] * det, inv[7] * det, inv[11] * det, inv[15] * det); +} + +/// Invert the transformation matrix. +/// This is much faster than the more general inverse() function, but assumes +/// that the matrix is of the form TR, where T is a translation and R a +/// rotation. +static inline mat4 mat4_inverse_transform(mat4 m) { + const vec3 r = mat4_v0(m); + const vec3 u = mat4_v1(m); + const vec3 f = mat4_v2(m); + const vec3 t = mat4_v3(m); + return mat4_make(r.x, r.y, r.z, -vec3_dot(r, t), u.x, u.y, u.z, + -vec3_dot(u, t), f.x, f.y, f.z, -vec3_dot(f, t), 0.0, 0.0, + 0.0, 1.0); +} + +/// Transpose the matrix. +static inline mat4 mat4_transpose(mat4 m) { + return mat4_make( + mat4_at(m, 0, 0), mat4_at(m, 1, 0), mat4_at(m, 2, 0), mat4_at(m, 3, 0), + mat4_at(m, 0, 1), mat4_at(m, 1, 1), mat4_at(m, 2, 1), mat4_at(m, 3, 1), + mat4_at(m, 0, 2), mat4_at(m, 1, 2), mat4_at(m, 2, 2), mat4_at(m, 3, 2), + mat4_at(m, 0, 3), mat4_at(m, 1, 3), mat4_at(m, 2, 3), mat4_at(m, 3, 3)); +} + +/// Transform the vector with the matrix. +static inline vec3 mat4_mul_vec3(mat4 m, vec3 v, R w) { + vec3 u; + u.x = mat4_at(m, 0, 0) * v.x + mat4_at(m, 0, 1) * v.y + + mat4_at(m, 0, 2) * v.z + mat4_at(m, 0, 3) * w; + u.y = mat4_at(m, 1, 0) * v.x + mat4_at(m, 1, 1) * v.y + + mat4_at(m, 1, 2) * v.z + mat4_at(m, 1, 3) * w; + u.z = mat4_at(m, 2, 0) * v.x + mat4_at(m, 2, 1) * v.y + + mat4_at(m, 2, 2) * v.z + mat4_at(m, 2, 3) * w; + return u; +} + +/// Return the vector multiplied by the matrix. +static inline vec4 mat4_mul_vec4(mat4 m, vec4 v) { + vec4 u; + u.x = mat4_at(m, 0, 0) * v.x + mat4_at(m, 0, 1) * v.y + + mat4_at(m, 0, 2) * v.z + mat4_at(m, 0, 3) * v.w; + u.y = mat4_at(m, 1, 0) * v.x + mat4_at(m, 1, 1) * v.y + + mat4_at(m, 1, 2) * v.z + mat4_at(m, 1, 3) * v.w; + u.z = mat4_at(m, 2, 0) * v.x + mat4_at(m, 2, 1) * v.y + + mat4_at(m, 2, 2) * v.z + mat4_at(m, 2, 3) * v.w; + u.w = mat4_at(m, 3, 0) * v.x + mat4_at(m, 3, 1) * v.y + + mat4_at(m, 3, 2) * v.z + mat4_at(m, 3, 3) * v.w; + return u; +} + +/// Compare two matrices for equality. +/// Returns true if the difference between each ij-value across matrices is +/// within |eps|, false if there is at least one ij-value difference that is +/// greater than eps. +static inline bool mat4_eq(mat4 m, mat4 w, float eps) { + return (float_eq(mat4_at(m, 0, 0), mat4_at(w, 0, 0), eps) && + float_eq(mat4_at(m, 0, 1), mat4_at(w, 0, 1), eps) && + float_eq(mat4_at(m, 0, 2), mat4_at(w, 0, 2), eps) && + float_eq(mat4_at(m, 0, 3), mat4_at(w, 0, 3), eps) && + + float_eq(mat4_at(m, 1, 0), mat4_at(w, 1, 0), eps) && + float_eq(mat4_at(m, 1, 1), mat4_at(w, 1, 1), eps) && + float_eq(mat4_at(m, 1, 2), mat4_at(w, 1, 2), eps) && + float_eq(mat4_at(m, 1, 3), mat4_at(w, 1, 3), eps) && + + float_eq(mat4_at(m, 2, 0), mat4_at(w, 2, 0), eps) && + float_eq(mat4_at(m, 2, 1), mat4_at(w, 2, 1), eps) && + float_eq(mat4_at(m, 2, 2), mat4_at(w, 2, 2), eps) && + float_eq(mat4_at(m, 2, 3), mat4_at(w, 2, 3), eps) && + + float_eq(mat4_at(m, 3, 0), mat4_at(w, 3, 0), eps) && + float_eq(mat4_at(m, 3, 1), mat4_at(w, 3, 1), eps) && + float_eq(mat4_at(m, 3, 2), mat4_at(w, 3, 2), eps) && + float_eq(mat4_at(m, 3, 3), mat4_at(w, 3, 3), eps)); +} diff --git a/include/math/spatial3.h b/include/math/spatial3.h new file mode 100644 index 0000000..f8caf5d --- /dev/null +++ b/include/math/spatial3.h @@ -0,0 +1,194 @@ +#pragma once + +#include "mat4.h" +#include "vec3.h" + +/// An object in 3D space. +typedef struct Spatial3 { + vec3 p; // Position. + vec3 r; // Right vector. + vec3 u; // Up vector. + vec3 f; // Forward vector. +} Spatial3; + +/// Construct a spatial with position 0 and canonical direction vectors. +static inline Spatial3 spatial3_make() { + return (Spatial3){.p = zero3(), .r = right3(), .u = up3(), .f = forward3()}; +} + +/// Return the spatial's transformation matrix (from local to world +/// coordinates). +static inline mat4 spatial3_transform(const Spatial3* spatial) { + const vec3 p = spatial->p; + const vec3 r = spatial->r; + const vec3 u = spatial->u; + const vec3 f = spatial->f; + return mat4_make(r.x, u.x, -f.x, p.x, r.y, u.y, -f.y, p.y, r.z, u.z, -f.z, + p.z, 0.0, 0.0, 0.0, 1.0f); +} + +/// Return the spatial's inverse transformation matrix (from world to local +/// coordinates). +static inline mat4 spatial3_inverse_transform(const Spatial3* spatial) { + return mat4_inverse_transform(spatial3_transform(spatial)); +} + +/// Move the spatial by the given vector. +static inline void spatial3_move(Spatial3* spatial, vec3 v) { + spatial->p = vec3_add(spatial->p, v); +} + +/// Move the spatial along its forward vector. +static inline void spatial3_move_forwards(Spatial3* spatial, R speed) { + spatial->p = vec3_add(spatial->p, vec3_scale(spatial->f, speed)); +} + +/// Move the spatial along its backwards vector. +static inline void spatial3_move_backwards(Spatial3* spatial, R speed) { + spatial->p = vec3_add(spatial->p, vec3_scale(spatial->f, -speed)); +} + +/// Move the spatial along its left vector. +static inline void spatial3_move_left(Spatial3* spatial, R speed) { + spatial->p = vec3_add(spatial->p, vec3_scale(spatial->r, -speed)); +} + +/// Move the spatial along its right vector. +static inline void spatial3_move_right(Spatial3* spatial, R speed) { + spatial->p = vec3_add(spatial->p, vec3_scale(spatial->r, speed)); +} + +/// Move the spatial along the global up vector. +static inline void spatial3_move_up(Spatial3* spatial, R speed) { + spatial->p = vec3_add(spatial->p, vec3_scale(up3(), speed)); +} + +/// Move the spatial along the global down vector. +static inline void spatial3_move_down(Spatial3* spatial, R speed) { + spatial->p = vec3_add(spatial->p, vec3_scale(up3(), -speed)); +} + +/// Rotate the spatial about the given axis by the given angle. +static inline void spatial3_rotate(Spatial3* spatial, vec3 axis, R angle) { + mat4 transf = spatial3_transform(spatial); + const vec3 local_axis = + vec3_normalize(mat4_mul_vec3(mat4_inverse_transform(transf), axis, 0.0)); + transf = mat4_mul(transf, mat4_rot(local_axis, angle)); + spatial->r = vec3_normalize(mat4_v0(transf)); + spatial->u = vec3_normalize(mat4_v1(transf)); + spatial->f = vec3_normalize(vec3_neg(mat4_v2(transf))); +} + +/// Rotate the spatial about its local y axis. +static inline void spatial3_yaw(Spatial3* spatial, const R angle) { + const R sa = sin(angle); + const R ca = cos(angle); + spatial->f = vec3_normalize( + vec3_sub(vec3_scale(spatial->f, ca), vec3_scale(spatial->r, sa))); + spatial->r = vec3_normalize(vec3_cross(spatial->f, spatial->u)); +} + +/// Rotate the spatial about its local x axis. +static inline void spatial3_pitch(Spatial3* spatial, const R angle) { + const R sa = sin(angle); + const R ca = cos(angle); + spatial->f = vec3_normalize( + vec3_add(vec3_scale(spatial->f, ca), vec3_scale(spatial->u, sa))); + spatial->u = vec3_normalize(vec3_cross(spatial->r, spatial->f)); +} + +/// Rotate the spatial about its local z axis. +static inline void spatial3_roll(Spatial3* spatial, const R angle) { + const R sa = sin(angle); + const R ca = cos(angle); + spatial->u = vec3_normalize( + vec3_sub(vec3_scale(spatial->u, ca), vec3_scale(spatial->r, sa))); + spatial->r = vec3_normalize(vec3_cross(spatial->f, spatial->u)); +} + +/// Set the spatial's transformation matrix. +static inline void spatial3_set_transform(Spatial3* spatial, mat4 transform) { + spatial->r = mat4_v0(transform); + spatial->u = mat4_v1(transform); + spatial->f = mat4_v2(transform); + spatial->p = mat4_v3(transform); +} + +static inline void spatial3_set_forward(Spatial3* spatial, vec3 forward) { + spatial->f = vec3_normalize(forward); + // Use aux vector to define right vector orthogonal to forward. + if (vec3_eq(vec3_abs(spatial->f), up3())) { + spatial->r = vec3_normalize(vec3_cross(spatial->f, forward3())); + } else { + spatial->r = vec3_normalize(vec3_cross(spatial->f, up3())); + } + spatial->u = vec3_normalize(vec3_cross(spatial->r, spatial->f)); +} + +/// Make the spatial look at the given target. +static inline void spatial3_lookat(Spatial3* spatial, vec3 target) { + spatial3_set_forward(spatial, vec3_sub(target, spatial->p)); +} + +/// Make the spatial look at the given target. +static inline void spatial3_lookat_spatial(Spatial3* spatial, + const Spatial3* target) { + spatial3_set_forward(spatial, vec3_sub(target->p, spatial->p)); +} + +/// Make the spatial orbit around the given target. +/// \param target Target position. +/// \param radius Radial distance. +/// \param azimuth Azimuthal (horizontal) angle. +/// \param zenith Polar (vertical) angle. +static inline void spatial3_orbit(Spatial3* spatial, vec3 target, R radius, + R azimuth, R zenith) { + const R sx = sin(azimuth); + const R sy = sin(zenith); + const R cx = cos(azimuth); + const R cy = cos(zenith); + spatial->p = (vec3){target.x + radius * cy * sx, target.y + radius * sy, + target.z + radius * cx * cy}; +} + +/// Make the spatial orbit around the given target. +/// \param target Target spatial. +/// \param radius Radial distance. +/// \param azimuth Azimuthal (horizontal) angle. +/// \param zenith Polar (vertical) angle. +static inline void spatial3_orbit_spatial(Spatial3* spatial, + const Spatial3* target, R radius, + R azimuth, R zenith) { + spatial3_orbit(spatial, target->p, radius, azimuth, zenith); +} + +// The functions below are provided so that client code can work with a Spatial3 +// with no assumptions on the data type's definition. + +/// Return the spatial's position. +static inline vec3 spatial3_position(const Spatial3* spatial) { + return spatial->p; +} + +/// Return the spatial's right vector. +static inline vec3 spatial3_right(const Spatial3* spatial) { + return spatial->r; +} + +/// Return the spatial's up vector. +static inline vec3 spatial3_up(const Spatial3* spatial) { return spatial->u; } + +/// Return the spatial's forward vector. +static inline vec3 spatial3_forward(const Spatial3* spatial) { + return spatial->f; +} + +static inline void spatial3_set_position(Spatial3* spatial, vec3 p) { + spatial->p = p; +} + +static inline void spatial3_setx(Spatial3* spatial, R x) { spatial->p.x = x; } + +static inline void spatial3_sety(Spatial3* spatial, R y) { spatial->p.y = y; } + +static inline void spatial3_setz(Spatial3* spatial, R z) { spatial->p.z = z; } diff --git a/include/math/vec.h b/include/math/vec.h new file mode 100644 index 0000000..7b870e1 --- /dev/null +++ b/include/math/vec.h @@ -0,0 +1,16 @@ +// Common functions for vectors. +// +// This header file contains functions that operate with different kinds of +// vectors that would result in circular dependencies if defined in any of the +// vec2/3/4 header files. +#pragma once + +#include "vec2.h" +#include "vec3.h" +#include "vec4.h" + +/// Construct a 3D vector from a 2D one with z = 0. +static inline vec3 vec3_from_vec2(vec2 v) { return vec3{v.x, v.y, 0, 0}; } + +/// Project a 4D vector onto w=0. +static inline vec3 vec3_from_vec4(vec4 v) { return vec3{v.x, v.y, v.z, 0}; } diff --git a/include/math/vec2.h b/include/math/vec2.h new file mode 100644 index 0000000..0a3f692 --- /dev/null +++ b/include/math/vec2.h @@ -0,0 +1,10 @@ +#pragma once + +#include "defs.h" + +typedef struct vec2 { + R x, y; +} vec2; + +/// Construct a vector from 2 coordinates. +static inline vec2 vec2_make(R x, R y) { return (vec2){x, y}; } diff --git a/include/math/vec3.h b/include/math/vec3.h new file mode 100644 index 0000000..3c3b053 --- /dev/null +++ b/include/math/vec3.h @@ -0,0 +1,143 @@ +#pragma once + +#include "defs.h" + +#include +#include + +/// A 3D vector. +typedef struct vec3 { + R x, y, z; +} vec3; + +/// Construct a vector from 3 coordinates. +static inline vec3 vec3_make(R x, R y, R z) { return (vec3){x, y, z}; } + +/// Construct a vector from an array. +static inline vec3 vec3_from_array(const R xyz[3]) { + return (vec3){xyz[0], xyz[1], xyz[2]}; +} + +/// Construct a vector from a single scalar value. +/// x = y = z = val. +static inline vec3 vec3_from_scalar(R val) { return (vec3){val, val, val}; } + +/// Normalize the vector. +static inline vec3 vec3_normalize(vec3 v) { + R n = sqrt(v.x * v.x + v.y * v.y + v.z * v.z); + assert(n > 0); + return (vec3){v.x / n, v.y / n, v.z / n}; +} + +/// Return the vector's ith coordinate. +static inline R vec3_ith(vec3 v, int i) { + assert(i >= 0 && i < 3); + return ((const R*)&v)[i]; +} + +/// Negate the given vector. +static inline vec3 vec3_neg(vec3 v) { return (vec3){-v.x, -v.y, -v.z}; } + +/// Add two vectors. +static inline vec3 vec3_add(vec3 a, vec3 b) { + return (vec3){a.x + b.x, a.y + b.y, a.z + b.z}; +} + +/// Subtract two vectors. +static inline vec3 vec3_sub(vec3 a, vec3 b) { + return (vec3){a.x - b.x, a.y - b.y, a.z - b.z}; +} + +/// Modulate two vectors (component-wise multiplication). +static inline vec3 vec3_mul(vec3 a, vec3 b) { + return (vec3){a.x * b.x, a.y * b.y, a.z * b.z}; +} + +/// Divide two vectors component-wise. +static inline vec3 vec3_div(vec3 a, vec3 b) { + return (vec3){a.x / b.x, a.y / b.y, a.z / b.z}; +} + +/// Scale a vector by a scalar value. +static inline vec3 vec3_scale(vec3 v, R s) { + return (vec3){v.x * s, v.y * s, v.z * s}; +} + +/// Compare two vectors for equality. +static inline bool vec3_eq(vec3 a, vec3 b) { + return a.x == b.x && a.y == b.y && a.z == b.z; +} + +/// Return the absolute value of the vector. +static inline vec3 vec3_abs(vec3 v) { + return (vec3){rabs(v.x), rabs(v.y), rabs(v.z)}; +} + +/// Compare two vectors for inequality. +static inline bool vec3_ne(vec3 a, vec3 b) { return !(vec3_eq(a, b)); } + +/// Return the vector's squared magnitude. +static inline R vec3_norm2(vec3 v) { return v.x * v.x + v.y * v.y + v.z * v.z; } + +/// Return the vector's magnitude. +static inline R vec3_norm(vec3 v) { return sqrt(vec3_norm2(v)); } + +/// Return the squared distance between two points. +static inline R vec3_dist2(vec3 a, vec3 b) { + const vec3 v = vec3_sub(b, a); + return vec3_norm2(v); +} + +/// Return the distance between two points. +static inline R vec3_dist(vec3 a, vec3 b) { return sqrt(vec3_dist2(a, b)); } + +/// Return the given vector divided by its magnitude. +static inline vec3 normalize(vec3 v) { + const R n = vec3_norm(v); + assert(n > 0); + return (vec3){v.x / n, v.y / n, v.z / n}; +} + +/// Return the dot product of two vectors. +static inline R vec3_dot(vec3 a, vec3 b) { + return a.x * b.x + a.y * b.y + a.z * b.z; +} + +/// Return the cross product of two vectors. +static inline vec3 vec3_cross(vec3 a, vec3 b) { + return (vec3){a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, + a.x * b.y - a.y * b.x}; +} + +/// Reflect the vector about the normal. +static inline vec3 vec3_reflect(vec3 v, vec3 n) { + // r = v - 2 * dot(v, n) * n + return vec3_sub(v, vec3_scale(n, 2 * vec3_dot(v, n))); +} + +/// Refract the vector about the normal. +static inline vec3 refract(vec3 v, vec3 n, R e) { + // k = 1 - e^2(1 - dot(n,v) * dot(n,v)) + const R k = 1.0 - e * e * (1.0 - vec3_dot(n, v) * vec3_dot(n, v)); + assert(k >= 0); + // r = e*v - (e * dot(n,v) + sqrt(k)) * n + return vec3_sub(vec3_scale(v, e), + vec3_scale(n, e * vec3_dot(n, v) * sqrt(k))); +} + +/// Elevate the vector to a power. +static inline vec3 vec3_pow(vec3 v, R p) { + return (vec3){pow(v.x, p), pow(v.y, p), pow(v.z, p)}; +} + +/// The (1, 0, 0) vector. +static inline vec3 right3() { return (vec3){1.0, 0.0, 0.0}; } + +/// The (0, 1, 0) vector. +static inline vec3 up3() { return (const vec3){0.0, 1.0, 0.0}; } + +/// The (0, 0, -1) vector. +static inline vec3 forward3() { return (const vec3){0.0, 0.0, -1.0}; } + +/// The (0, 0, 0) vector. +static inline vec3 zero3() { return (const vec3){0.0, 0.0, 0.0}; } diff --git a/include/math/vec4.h b/include/math/vec4.h new file mode 100644 index 0000000..4ab843b --- /dev/null +++ b/include/math/vec4.h @@ -0,0 +1,15 @@ +#pragma once + +#include "defs.h" + +typedef struct vec4 { + R x, y, w, z; +} vec4; + +/// Construct a vector from 4 coordinates. +static inline vec4 vec4_make(R x, R y, R z, R w) { return (vec4){x, y, z, w}; } + +/// Construct a vector from an array. +static inline vec4 vec4_from_array(const R xyzw[4]) { + return (vec4){xyzw[0], xyzw[1], xyzw[2], xyzw[3]}; +} -- cgit v1.2.3