Skip to content

Instantly share code, notes, and snippets.

@LaylBongers

LaylBongers/matrix.h

Created November 27, 2018 14:30
Show Gist options
  • Save LaylBongers/7a861b54e8c6f0b6993632729eb2084c to your computer and use it in GitHub Desktop.
Save LaylBongers/7a861b54e8c6f0b6993632729eb2084c to your computer and use it in GitHub Desktop.
Download ZIP
Raw
matrix.h
#pragma once
#include <stdexcept>
namespace bmath {
template<typename S>
struct Mat4 {
// The matrix defines the direction of the axes of a local grid.
// This matrix is laid out like this:
// [x-axis-x y-axis-x z-axis-x w-axis-x]
// [x-axis-y y-axis-y z-axis-y w-axis-y]
// [x-axis-z y-axis-z z-axis-z w-axis-z]
// [x-axis-w y-axis-w z-axis-w w-axis-w]
// The w column is used for sheering in the 4th dimension, using this the matrix can represent translation.
// All vectors that are at 1 depth in the w dimension are translated, while vectors at 0 depth in the w
// dimension are not. Use the w component on a Vec4 to define if they should be translated or not. Positions
// should be translated, distances should not.
// In calls like the constructor, the separate lines of parameters are columns, not rows!
// Representations of the value
union {
S i[4][4]; // i[column][row]
S v[16]; // c0r0 c0r1 ... c3r2 c3r3
};
// Properties of this matrix type
enum {
width = 4,
height = 4
};
typedef S Scalar;
Mat4<S>() :
i{{1, 0, 0, 0},
{0, 1, 0, 0},
{0, 0, 1, 0},
{0, 0, 0, 1}} {
}
Mat4<S>(S c0r0, S c0r1, S c0r2, S c0r3,
S c1r0, S c1r1, S c1r2, S c1r3,
S c2r0, S c2r1, S c2r2, S c2r3,
S c3r0, S c3r1, S c3r2, S c3r3) :
i{{c0r0, c0r1, c0r2, c0r3},
{c1r0, c1r1, c1r2, c1r3},
{c2r0, c2r1, c2r2, c2r3},
{c3r0, c3r1, c3r2, c3r3}} {
}
};
// Matrix operations
namespace mat {
template<typename S>
Vec4 <S> row(const Mat4<S> &matrix, i32 index);
template<typename S>
Vec4 <S> col(const Mat4<S> &matrix, i32 index);
template<typename S>
Mat4<S> create_identity() {
return Mat4<S>();
}
template<typename S>
Mat4<S> create_perspective(f32 fovy, f32 aspect, f32 nearZ, f32 farZ) {
float_t f = fovy / 2.0f;
f = cosf(f) / sinf(f);
return Mat4<S>(
f / aspect, 0, 0, 0,
0, f, 0, 0,
0, 0, (farZ + nearZ) / (nearZ - farZ), -1.0f,
0, 0, (2.0f * farZ * nearZ) / (nearZ - farZ), 0
);
}
template<typename S>
Mat4<S> create_translation(const Vec3 <S> &amount) {
const auto zero = static_cast<S>(0);
const auto one = static_cast<S>(1);
return Mat4<S>(
one, zero, zero, zero,
zero, one, zero, zero,
zero, zero, one, zero,
amount.x, amount.y, amount.z, one
);
}
template<typename S>
Mat4<S> create_rotate_x(f32 angle) {
return Mat4<S>(
1.0f, 0.0f, 0.0f, 0.0f,
0.0f, cosf(angle), -sinf(angle), 0.0f,
0.0f, sinf(angle), cosf(angle), 0.0f,
0.0f, 0.0f, 0.0f, 1.0f
);
}
template<typename S>
Mat4<S> create_rotate_y(f32 angle) {
return Mat4<S>(
cosf(angle), 0.0f, sinf(angle), 0.0f,
0.0f, 1.0f, 0.0f, 0.0f,
-sinf(angle), 0.0f, cosf(angle), 0.0f,
0.0f, 0.0f, 0.0f, 1.0f
);
}
template<typename S>
Mat4<S> create_rotate_z(f32 angle) {
return Mat4<S>(
cosf(angle), -sinf(angle), 0.0f, 0.0f,
sinf(angle), cosf(angle), 0.0f, 0.0f,
0.0f, 0.0f, 1.0f, 0.0f,
0.0f, 0.0f, 0.0f, 1.0f
);
}
template<typename S>
Mat4<S> create_scale(f32 x, f32 y, f32 z) {
return Mat4<S>(
x, 0.0f, 0.0f, 0.0f,
0.0f, y, 0.0f, 0.0f,
0.0f, 0.0f, z, 0.0f,
0.0f, 0.0f, 0.0f, 1.0f
);
}
template<typename S>
Mat4<S> create_look_at(const Vec3<f32> &eye, const Vec3<f32> &center, const Vec3<f32> &up) {
auto f = vec::norm(vec::sub(center, eye));
auto s = vec::norm(vec::cross(f, up));
auto u = vec::cross(s, f);
return Mat4<S>(
s.x, u.x, -f.x, 0,
s.y, u.y, -f.y, 0,
s.z, u.z, -f.z, 0,
-vec::dot(eye, s), -vec::dot(eye, u), vec::dot(eye, f), 1
);
}
/// Transform a vector using a matrix
template<typename S>
Vec4 <S> transform(const Mat4<S> &matrix, const Vec4 <S> &vector) {
return Vec4<S>(
vec::dot(vector, mat::row(matrix, 0)),
vec::dot(vector, mat::row(matrix, 1)),
vec::dot(vector, mat::row(matrix, 2)),
vec::dot(vector, mat::row(matrix, 3))
);
}
/// Multiply a matrix with another matrix
template<typename S>
Mat4<S> mul(const Mat4<S> &left, const Mat4<S> &right) {
Mat4<f32> out;
for (i32 row_i = 0; row_i < 4; row_i++) {
auto row_v = row(left, row_i);
for (i32 col_i = 0; col_i < 4; col_i++) {
out.i[col_i][row_i] = vec::dot(row_v, col(right, col_i));
}
}
return out;
}
/// Get a row of a matrix
template<typename S>
Vec4 <S> row(const Mat4<S> &matrix, i32 index) {
return Vec4<S>(
matrix.i[0][index],
matrix.i[1][index],
matrix.i[2][index],
matrix.i[3][index]
);
}
/// Get a column of a matrix
template<typename S>
Vec4 <S> col(const Mat4<S> &matrix, i32 index) {
return Vec4<S>(
matrix.i[index][0],
matrix.i[index][1],
matrix.i[index][2],
matrix.i[index][3]
);
}
template<typename S>
Mat4 <S> invert(const Mat4<S> &matrix) {
// Adapted from the MESA implementation
auto inverse = Mat4<S>();
const S *m = matrix.v;
S *inv = inverse.v;
double det;
int i;
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];
det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
if (det == 0) {
throw std::logic_error("This matrix does not have an inverse");
}
det = 1.0 / det;
for (i = 0; i < 16; i++) {
inv[i] = inv[i] * det;
}
return inverse;
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment