jc/math/mat.jai
2025-09-06 00:52:48 -06:00

483 lines
13 KiB
Text

Mat2 :: #type,distinct Vec(2 * 2, float);
Mat4 :: #type,distinct Vec(4 * 4, float);
m2 :: (_xx: float, xy: float, yx: float, yy: float) -> Mat2 #expand {
m: Mat2 = ---;
m._00 = _xx;
m._01 = xy;
m._10 = yx;
m._11 = yy;
return m;
}
m2 :: (r0: Vec2, r1: Vec2) -> Mat2 #expand {
return .{
components = .[
r0.x, r0.y,
r1.x, r1.y,
]
};
}
m4 :: (r0: Vec4, r1: Vec4, r2: Vec4, r3: Vec4) -> Mat4 #expand {
return .{
components = .[
r0.x, r0.y, r0.z, r0.w,
r1.x, r1.y, r1.z, r1.w,
r2.x, r2.y, r2.z, r2.w,
r3.x, r3.y, r3.z, r3.w,
],
};
}
m2_identity :: Mat2.{_00=1, _11=1};
m4_identity :: Mat4.{_00=1, _11=1, _22=1, _33=1};
m2d :: (diag: float) -> Mat2 #expand {
r: Mat2;
r._00 = diag;
r._11 = diag;
return r;
}
m4d :: (diag: float) -> Mat4 #expand {
res: Mat4;
res._00 = diag;
res._11 = diag;
res._22 = diag;
res._33 = diag;
return res;
}
operator [] :: inline (m: Mat2, idx: int) -> Vec2 #expand {
bounds_check_index(idx * 2, m.N);
return (m.components[idx * 2]).(*Vec4).*;
}
operator *[] :: inline (m: *Mat2, $$idx: int) -> *Vec2 #no_abc {
N :: Mat2.{}.N; // @note(judah): dumb workaround for not being able to access pointer type field constants
bounds_check_index(idx * 2, N);
return (*m.components[idx * 2]).(*Vec4);
}
operator []= :: inline (m: *Mat2, $$idx: int, value: Vec2) #no_abc {
N :: Mat2.{}.N;
bounds_check_index(idx * 4, N);
ptr := (*m.components[idx * 4]).(*Vec2);
ptr.*= value;
}
operator* :: inline (a: Mat2, b: Mat2) -> Mat2 {
r: Mat2 = ---;
r._00 = a._00*b._00 + a._01*b._10;
r._01 = a._00*b._01 + a._01*b._11;
r._10 = a._10*b._00 + a._11*b._10;
r._11 = a._10*b._01 + a._11*b._11;
return r;
}
operator* :: inline (a: Mat2, b: Vec2) -> Vec2 {
v: Vec2 = ---;
v.x = a._00*b.x + a._01*b.y;
v.y = a._10*b.x + a._11*b.y;
return v;
}
// Note(Jesse): Probably never used except for testing
operator== :: (a: Mat2, b: Mat2) -> bool #no_abc {
for a.components if it != b.components[it_index] return false;
return true;
}
operator [] :: inline (m: Mat4, idx: int) -> Vec4 #expand {
meta.check_bounds(idx * 4, m.N);
return (m.components[idx * 4]).(*Vec4).*;
}
operator *[] :: inline (m: *Mat4, $$idx: int) -> *Vec4 #no_abc {
N :: Mat4.{}.N; // @note(judah): dumb workaround for not being able to access pointer type field constants
meta.check_bounds(idx * 4, N);
return (*m.components[idx * 4]).(*Vec4);
}
operator []= :: inline (m: *Mat4, $$idx: int, value: Vec4) #no_abc {
N :: Mat4.{}.N;
meta.check_bounds(idx * 4, N);
ptr := (*m.components[idx * 4]).(*Vec4);
ptr.*= value;
}
// Note(Jesse): Probably never used except for testing
operator== :: (a: Mat4, b: Mat4) -> bool #no_abc {
for a.components if it != b.components[it_index] return false;
return true;
}
translate :: (from: Mat4, with: Vec3) -> Mat4 {
r := from;
r._03 += with.x*r._00 + with.y*r._01 + with.z*r._02;
r._13 += with.x*r._10 + with.y*r._11 + with.z*r._12;
r._23 += with.x*r._20 + with.y*r._21 + with.z*r._22;
return r;
}
translate :: (with: Vec3) -> Mat4 {
r: Mat4 = m4_identity;
r._03 = with.x;
r._13 = with.y;
r._23 = with.z;
return r;
}
make_look_at :: (camera: Vec3, at: Vec3, up_vector: Vec3) -> Mat4
{
// -Z forward, Y up, X right
forward := normalize(at - camera);
right := normalize(cross(up_vector, forward)); // up_vector
up := normalize(cross(forward, right));
look_at: Mat4;
look_at._00 = right.x;
look_at._01 = right.y;
look_at._02 = right.z;
look_at._10 = up.x;
look_at._11 = up.y;
look_at._12 = up.z;
look_at._20 = -forward.x;
look_at._21 = -forward.y;
look_at._22 = -forward.z;
look_at._33 = 1;
look_at = inline translate(look_at, -camera);
return look_at;
}
operator* :: (a: Mat4, b: Mat4) -> Mat4 {
r: Mat4 = ---;
r._00 = a._00*b._00 + a._01*b._10 + a._02*b._20 + a._03*b._30;
r._01 = a._00*b._01 + a._01*b._11 + a._02*b._21 + a._03*b._31;
r._02 = a._00*b._02 + a._01*b._12 + a._02*b._22 + a._03*b._32;
r._03 = a._00*b._03 + a._01*b._13 + a._02*b._23 + a._03*b._33;
r._10 = a._10*b._00 + a._11*b._10 + a._12*b._20 + a._13*b._30;
r._11 = a._10*b._01 + a._11*b._11 + a._12*b._21 + a._13*b._31;
r._12 = a._10*b._02 + a._11*b._12 + a._12*b._22 + a._13*b._32;
r._13 = a._10*b._03 + a._11*b._13 + a._12*b._23 + a._13*b._33;
r._20 = a._20*b._00 + a._21*b._10 + a._22*b._20 + a._23*b._30;
r._21 = a._20*b._01 + a._21*b._11 + a._22*b._21 + a._23*b._31;
r._22 = a._20*b._02 + a._21*b._12 + a._22*b._22 + a._23*b._32;
r._23 = a._20*b._03 + a._21*b._13 + a._22*b._23 + a._23*b._33;
r._30 = a._30*b._00 + a._31*b._10 + a._32*b._20 + a._33*b._30;
r._31 = a._30*b._01 + a._31*b._11 + a._32*b._21 + a._33*b._31;
r._32 = a._30*b._02 + a._31*b._12 + a._32*b._22 + a._33*b._32;
r._33 = a._30*b._03 + a._31*b._13 + a._32*b._23 + a._33*b._33;
return r;
}
// Note(Jesse): If you want to make it symmetric go ahead, I usually don't
operator* :: (a: Mat4, v: Vec4) -> Vec4 {
r: Vec4 = ---;
r.x = a._00*v.x + a._01*v.y + a._02*v.z + a._03*v.w;
r.y = a._10*v.x + a._11*v.y + a._12*v.z + a._13*v.w;
r.z = a._20*v.x + a._21*v.y + a._22*v.z + a._23*v.w;
r.w = a._30*v.x + a._31*v.y + a._32*v.z + a._33*v.w;
return r;
}
transpose :: (m: Mat2) -> Mat2 {
r: Mat2 = ---;
// Todo(Jesse): This can be one call to simd shuffle
r._00 = m._00;
r._01 = m._10;
r._10 = m._01;
r._11 = m._11;
return r;
}
transpose :: (m: Mat4) -> Mat4 {
r: Mat4 = ---;
r._00 = m._00;
r._01 = m._10;
r._02 = m._20;
r._03 = m._30;
r._10 = m._01;
r._11 = m._11;
r._12 = m._21;
r._13 = m._31;
r._20 = m._02;
r._21 = m._12;
r._22 = m._22;
r._23 = m._32;
r._30 = m._03;
r._31 = m._13;
r._32 = m._23;
r._33 = m._33;
return r;
}
// aspect is width/height
// NDC Depth depth_01
// DirectX: 0-1 true
// Vulkan: 0-1 true
// Opengl: -1-1 false
perspective_projection :: (aspect: float, fov: float, near: float, far: float, $$depth_01: bool) -> Mat4 {
r: Mat4 = ---;
// Y up, X right, Z- forward
fov_part := 1.0/tan(fov*0.5);
denom := (far - near);
r._00 = fov_part/aspect;
r._11 = fov_part;
r._22 = -(far + near)/denom;
r._23 = -2*far*near/denom;
r._32 = -1;
r._33 = 1;
if depth_01 {
// Makes depth range from 0 to 1
// To map -1,1 depth range to 0,1 we transform z as follows: z' = z * 0.5 + 0.5
r._22 = r._22 * 0.5 + r._32 * 0.5;
r._23 = r._23 * 0.5 + r._33 * 0.5;
}
return r;
}
perspective_projection_left_handed :: (aspect: float, fov: float, near: float, far: float, $$depth_01: bool) -> Mat4 {
r: Mat4 = ---;
// Y up, X right, Z- forward
fov_part := 1.0/tan(fov*0.5);
denom := (far - near);
r._00 = fov_part/aspect;
r._11 = fov_part;
r._22 = (far + near)/denom;
r._23 = 2*far*near/denom;
r._32 = -1;
r._33 = 1;
if depth_01 {
// Makes depth range from 0 to 1
// To map -1,1 depth range to 0,1 we transform z as follows: z' = z * 0.5 + 0.5
r._22 = r._22 * 0.5 + r._32 * 0.5;
r._23 = r._23 * 0.5 + r._33 * 0.5;
}
return r;
}
// Note(Jesse): Typically left/right/top/bottom are pixel coords
// NDC Depth depth_01
// DirectX: 0-1 true
// Vulkan: 0-1 true
// Opengl: -1-1 false
orthographic_projection :: (left: float, right: float, top: float, bottom: float, near: float, far: float, $$depth_01: bool) -> Mat4 {
r: Mat4 = m4_identity;
r._00 = 2.0/(right - left);
r._11 = 2.0/(top - bottom);
r._22 = -1.0/(far - near);
if depth_01 == false then r._22 *= 2;
r._30 = -(right + left)/(right - left);
r._31 = -(top + bottom)/(top - bottom);
r._32 = -(far + near)/(far - near);
r._33 = 1;
return r;
}
rotation_mat2 :: (angle: float) -> Mat2 {
m: Mat2 = ---;
m._00 = cos(angle);
m._01 = -sin(angle);
m._10 = sin(angle);
m._11 = cos(angle);
return m;
}
rotation_mat4 :: (axis: Vec3, angle: float) -> Mat4 {
m := m4_identity;
a := normalize(axis);
sin_theta := sin(angle);
cos_theta := cos(angle);
inv_cos := 1.0 - cos_theta;
// Contributors to vector's X
m._00 = (a.x*a.x*inv_cos) + cos_theta;
m._01 = (a.x*a.y*inv_cos) + (a.z*sin_theta);
m._02 = (a.x*a.z*inv_cos) - (a.y*sin_theta);
// Contributors to vector's Y
m._10 = (a.y*a.x*inv_cos) - (a.z*sin_theta);
m._11 = (a.y*a.y*inv_cos) + cos_theta;
m._12 = (a.y*a.z*inv_cos) + (a.x*sin_theta);
// Contributors to vector's Z
m._20 = (a.z*a.x*inv_cos) + (a.y*sin_theta);
m._21 = (a.z*a.y*inv_cos) - (a.x*sin_theta);
m._22 = (a.z*a.z*inv_cos) + cos_theta;
return m;
}
rotation_mat4 :: (q: Quat) -> Mat4 {
m := m4_identity;
xs := 2*q.x;
ys := 2*q.y;
zs := 2*q.z;
wx := q.w*xs;
wy := q.w*ys;
wz := q.w*zs;
_xx := q.x*xs;
xy := q.x*ys;
xz := q.x*zs;
yy := q.y*ys;
yz := q.y*zs;
zz := q.z*zs;
m._00 = 1.0 - (yy + zz);
m._01 = xy - wz;
m._02 = xz + wy;
m._10 = xy + wz;
m._11 = 1.0 - (_xx + zz);
m._12 = yz - wx;
m._20 = xz - wy;
m._21 = yz + wx;
m._22 = 1.0 - (_xx + yy);
return m;
}
determinate :: (m: Mat2) -> float {
return m._00*m._11 - m._01*m._10;
}
determinate :: (m: Mat4) -> float {
c01 := cross(m.v4[0].xyz, m.v4[1].xyz);
c23 := cross(m.v4[2].xyz, m.v4[3].xyz);
b10 := m.v4[0].xyz*m.v4[1].w - m.v4[1].xyz*m.v4[0].w;
b32 := m.v4[2].xyz*m.v4[3].w - m.v4[3].xyz*m.v4[2].w;
return dot(c01, b32) + dot(c23, b10);
}
inverse :: (m: Mat2) -> Mat2 {
r: Mat2 = ---;
invdet := 1.0/determinate(m);
r._00 = invdet* m._11;
r._01 = invdet*-m._01;
r._10 = invdet*-m._10;
r._11 = invdet* m._00;
return r;
}
inverse :: (m: Mat4) -> Mat4 {
c01 := cross(m.v4[0].xyz, m.v4[1].xyz);
c23 := cross(m.v4[2].xyz, m.v4[3].xyz);
b10 := m.v4[0].xyz*m.v4[1].w - m.v4[1].xyz*m.v4[0].w;
b32 := m.v4[2].xyz*m.v4[3].w - m.v4[3].xyz*m.v4[2].w;
det := dot(c01, b32) + dot(c23, b10);
if abs(det) < 0.0001
return m4_identity;
inv_det := 1.0/det;
c01 *= inv_det;
c23 *= inv_det;
b10 *= inv_det;
b32 *= inv_det;
r: Mat4 = ---;
r.v4[0] = v4f(cross(m.v4[1].xyz, b32) + c23*m.v4[1].w, -dot(m.v4[1].xyz, c23));
r.v4[1] = v4f(cross(b32, m.v4[0].xyz) + c23*m.v4[0].w, dot(m.v4[0].xyz, c23));
r.v4[2] = v4f(cross(m.v4[3].xyz, b10) + c01*m.v4[3].w, -dot(m.v4[3].xyz, c01));
r.v4[3] = v4f(cross(b10, m.v4[2].xyz) + c01*m.v4[2].w, dot(m.v4[2].xyz, c01));
return transpose(r);
}
#scope_file;
#if RUN_TESTS #run {
Test(basic.tprint("%: Mat2", UNITS), t => {
{
identity := m2_identity;
a: Mat2 = Mat2.{.[1, 2, 3, 4]};
Expect(a*identity == a);
}
{
rotator := rotation_mat2(from_rad(PI));
v2 := v2f(1.0, 0.5);
v2 = rotator*v2;
Expect(v2_eq(v2, v2f(-1.0, -0.5)));
v2 = rotator*v2;
Expect(v2_eq(v2, v2f(1.0, 0.5)));
}
{
m := m2(1.0, 3.0, 2.0, 2.0);
det := determinate(m);
Expect(det == -4);
}
{
rotator := rotation_mat2(from_rad(PI));
inv_rotator := inverse(rotator);
v2 := v2f(0.1, 1.0);
Expect(inv_rotator*rotator == m2_identity);
Expect(inv_rotator*(rotator*v2) == v2);
v2 = rotator*v2;
v2 = inv_rotator*v2;
}
});
Test(basic.tprint("%: Mat4", UNITS), t => {
{
identity := m4_identity;
a: Mat4 = Mat4.{.[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]};
Expect(a*identity == a);
}
{
UP_VECTOR :: Vec3.{.[0.0, 1.0, 0.0]};
camera := v3f(1.0, 0.0, 1.0);
looking_at := v3f(0.0, 0.0, 0.0);
look_at := make_look_at(camera, looking_at, UP_VECTOR);
}
{
translator := translate(v3f(10.0, 5.0, 2.0));
v3 := v3f(1.0, 2.0, 1.0);
v4 := v4f(v3, 1.0);
Expect(v4 == v4f(1.0, 2.0, 1.0, 1.0));
Expect(translator*v4 == v4f(11.0, 7.0, 3.0, 1.0));
}
{
rotator := rotation_mat4(v3f(0.0, 1.0, 0.0), from_rad(PI));
v4 := v4f(1.0, 0.5, 0.1, 1.0);
v4 = rotator*v4;
Expect(v4_eq(v4, v4f(-1.0, 0.5, -0.1, 1.0)));
v4 = rotator*v4;
Expect(v4_eq(v4, v4f(1.0, 0.5, 0.1, 1.0)));
}
{
rotator_x := rotation_mat4(v3f(1.0, 0.0, 0.0), from_rad(0.4*PI));
rotator_y := rotation_mat4(v3f(0.0, 1.0, 0.0), from_rad(PI));
camera_rotator := rotator_x*rotator_y;
det := determinate(camera_rotator);
}
{
rotator := rotation_mat4(v3f(0.0, 0.0, 1.0), from_rad(PI));
inv_rotator := inverse(rotator);
v4 := v4f(0.1, 1.0, 0.0, 1.0);
Expect(inv_rotator*rotator == m4_identity);
v4 = rotator*v4;
v4 = inv_rotator*v4;
}
});
}