From 33a458c520019e4d09c5a442524fd5eaaf6d90fe Mon Sep 17 00:00:00 2001 From: gingerBill Date: Sat, 28 Dec 2019 23:00:13 +0000 Subject: [PATCH] Update package math/linalg --- core/math/linalg/linalg.odin | 819 ++++++++++++++++++++++++++++++----- 1 file changed, 719 insertions(+), 100 deletions(-) diff --git a/core/math/linalg/linalg.odin b/core/math/linalg/linalg.odin index 1220b952e..2009d478e 100644 --- a/core/math/linalg/linalg.odin +++ b/core/math/linalg/linalg.odin @@ -5,76 +5,292 @@ import "intrinsics" // Generic -dot_vector :: proc(a, b: $T/[$N]$E) -> (c: E) { +@private IS_NUMERIC :: intrinsics.type_is_numeric; +@private IS_QUATERNION :: intrinsics.type_is_quaternion; +@private IS_ARRAY :: intrinsics.type_is_array; + + +vector_dot :: proc(a, b: $T/[$N]$E) -> (c: E) where IS_NUMERIC(E) { for i in 0.. (c: f32) { +quaternion128_dot :: proc(a, b: $T/quaternion128) -> (c: f32) { return real(a)*real(a) + imag(a)*imag(b) + jmag(a)*jmag(b) + kmag(a)*kmag(b); } -dot_quaternion256 :: proc(a, b: $T/quaternion256) -> (c: f64) { +quaternion256_dot :: proc(a, b: $T/quaternion256) -> (c: f64) { return real(a)*real(a) + imag(a)*imag(b) + jmag(a)*jmag(b) + kmag(a)*kmag(b); } -dot :: proc{dot_vector, dot_quaternion128, dot_quaternion256}; +dot :: proc{vector_dot, quaternion128_dot, quaternion256_dot}; -cross2 :: proc(a, b: $T/[2]$E) -> E { +quaternion_inverse :: proc(q: $Q) -> Q where IS_QUATERNION(Q) { + return conj(q) * quaternion(1.0/dot(q, q), 0, 0, 0); +} + + +vector_cross2 :: proc(a, b: $T/[2]$E) -> E where IS_NUMERIC(E) { return a[0]*b[1] - b[0]*a[1]; } -cross3 :: proc(a, b: $T/[3]$E) -> (c: T) { - c[0] = +(a[1]*b[2] - b[1]*a[2]); - c[1] = -(a[2]*b[0] - b[2]*a[0]); - c[2] = +(a[0]*b[1] - b[0]*a[1]); +vector_cross3 :: proc(a, b: $T/[3]$E) -> (c: T) where IS_NUMERIC(E) { + c[0] = a[1]*b[2] - b[1]*a[2]; + c[1] = a[2]*b[0] - b[2]*a[0]; + c[2] = a[0]*b[1] - b[0]*a[1]; return; } -cross :: proc{cross2, cross3}; +vector_cross :: proc{vector_cross2, vector_cross3}; +cross :: vector_cross; - -normalize_vector :: proc(v: $T/[$N]$E) -> T { +vector_normalize :: proc(v: $T/[$N]$E) -> T where IS_NUMERIC(E) { return v / length(v); } -normalize_quaternion128 :: proc(q: $Q/quaternion128) -> Q { +quaternion_normalize :: proc(q: $Q) -> Q where IS_QUATERNION(Q) { return q/abs(q); } -normalize_quaternion256 :: proc(q: $Q/quaternion256) -> Q { - return q/abs(q); -} -normalize :: proc{normalize_vector, normalize_quaternion128, normalize_quaternion256}; +normalize :: proc{vector_normalize, quaternion_normalize}; -normalize0_vector :: proc(v: $T/[$N]$E) -> T { +vector_normalize0 :: proc(v: $T/[$N]$E) -> T where IS_NUMERIC(E) { m := length(v); return m == 0 ? 0 : v/m; } -normalize0_quaternion128 :: proc(q: $Q/quaternion128) -> Q { +quaternion_normalize0 :: proc(q: $Q) -> Q where IS_QUATERNION(Q) { m := abs(q); return m == 0 ? 0 : q/m; } -normalize0_quaternion256 :: proc(q: $Q/quaternion256) -> Q { - m := abs(q); - return m == 0 ? 0 : q/m; -} -normalize0 :: proc{normalize0_vector, normalize0_quaternion128, normalize0_quaternion256}; +normalize0 :: proc{vector_normalize0, quaternion_normalize0}; -length :: proc(v: $T/[$N]$E) -> E { +vector_length :: proc(v: $T/[$N]$E) -> E where IS_NUMERIC(E) { return math.sqrt(dot(v, v)); } -length2 :: proc(v: $T/[$N]$E) -> E { +vector_length2 :: proc(v: $T/[$N]$E) -> E where IS_NUMERIC(E) { return dot(v, v); } +quaternion_length :: proc(q: $Q) -> Q where IS_QUATERNION(Q) { + return abs(q); +} + +quaternion_length2 :: proc(q: $Q) -> Q where IS_QUATERNION(Q) { + return dot(q, q); +} + +length :: proc{vector_length, quaternion_length}; +length2 :: proc{vector_length2, quaternion_length2}; + + + +vector_sin :: proc(angle: $V/[$N]$E) -> V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + s: V; + for i in 0.. V where IS_NUMERIC(E) { + return length(p1 - p0); +} + +vector_reflect :: proc(i, n: $V/[$N]$E) -> V where IS_NUMERIC(E) { + b := n * (2 * dot(n, i)); + return i - b; +} + +vector_refract :: proc(i, n: $V/[$N]$E, eta: E) -> V where IS_NUMERIC(E) { + dv := dot(n, i); + k := 1 - eta*eta - (1 - dv*dv); + a := i * eta; + b := n * eta*dv*math.sqrt(k); + return (a - b) * E(int(k >= 0)); +} + + identity :: proc($T: typeid/[$N][N]$E) -> (m: T) { for i in 0.. (m: [M][N]E) { +transpose :: proc(a: $T/[$N][$M]$E) -> (m: T) { for j in 0.. (m: [M][N]E) { return; } -mul_matrix :: proc(a, b: $M/[$N][N]$E) -> (c: M) - where !intrinsics.type_is_array(E), - intrinsics.type_is_numeric(E) { +matrix_mul :: proc(a, b: $M/[$N][N]$E) -> (c: M) + where !IS_ARRAY(E), + IS_NUMERIC(E) { for i in 0.. (c: M) return; } -mul_matrix_differ :: proc(a: $A/[$J][$I]$E, b: $B/[$K][J]E) -> (c: [K][I]E) - where !intrinsics.type_is_array(E), - intrinsics.type_is_numeric(E), - I != K { +matrix_mul_differ :: proc(a: $A/[$J][$I]$E, b: $B/[$K][J]E) -> (c: [K][I]E) + where !IS_ARRAY(E), + IS_NUMERIC(E), + I != K { for k in 0.. (c: [K][I]E) } -mul_matrix_vector :: proc(a: $A/[$I][$J]$E, b: $B/[I]E) -> (c: B) - where !intrinsics.type_is_array(E), - intrinsics.type_is_numeric(E) { +matrix_mul_vector :: proc(a: $A/[$I][$J]$E, b: $B/[I]E) -> (c: B) + where !IS_ARRAY(E), + IS_NUMERIC(E) { for i in 0.. (c: B) return; } -mul_quaternion128_vector3 :: proc(q: $Q/quaternion128, v: $V/[3]$F/f32) -> V { +quaternion128_mul_vector3 :: proc(q: $Q/quaternion128, v: $V/[3]$F/f32) -> V { Raw_Quaternion :: struct {xyz: [3]f32, r: f32}; q := transmute(Raw_Quaternion)q; @@ -132,7 +348,7 @@ mul_quaternion128_vector3 :: proc(q: $Q/quaternion128, v: $V/[3]$F/f32) -> V { return V(v + q.r*t + cross(q.xyz, t)); } -mul_quaternion256_vector3 :: proc(q: $Q/quaternion256, v: $V/[3]$F/f64) -> V { +quaternion256_mul_vector3 :: proc(q: $Q/quaternion256, v: $V/[3]$F/f64) -> V { Raw_Quaternion :: struct {xyz: [3]f64, r: f64}; q := transmute(Raw_Quaternion)q; @@ -141,16 +357,23 @@ mul_quaternion256_vector3 :: proc(q: $Q/quaternion256, v: $V/[3]$F/f64) -> V { t := cross(2*q.xyz, v); return V(v + q.r*t + cross(q.xyz, t)); } -mul_quaternion_vector3 :: proc{mul_quaternion128_vector3, mul_quaternion256_vector3}; +quaternion_mul_vector3 :: proc{quaternion128_mul_vector3, quaternion256_mul_vector3}; mul :: proc{ - mul_matrix, - mul_matrix_differ, - mul_matrix_vector, - mul_quaternion128_vector3, - mul_quaternion256_vector3, + matrix_mul, + matrix_mul_differ, + matrix_mul_vector, + quaternion128_mul_vector3, + quaternion256_mul_vector3, }; +vector_to_ptr :: proc(v: ^$V/[$N]$E) -> ^E where IS_NUMERIC(E) { + return &v[0]; +} +matrix_to_ptr :: proc(m: ^$A/[$I][$J]$E) -> ^E where IS_NUMERIC(E) { + return &m[0][0]; +} + // Specific @@ -199,6 +422,11 @@ VECTOR3_Y_AXIS :: Vector3{0, 1, 0}; VECTOR3_Z_AXIS :: Vector3{0, 0, 1}; + +vector2_orthogonal :: proc(v: Vector2) -> Vector2 { + return {-v.y, v.x}; +} + vector3_orthogonal :: proc(v: Vector3) -> Vector3 { x := abs(v.x); y := abs(v.y); @@ -206,20 +434,442 @@ vector3_orthogonal :: proc(v: Vector3) -> Vector3 { other: Vector3 = x < y ? (x < z ? {1, 0, 0} : {0, 0, 1}) : (y < z ? {0, 1, 0} : {0, 0, 1}); - return normalize(cross3(v, other)); + return normalize(cross(v, other)); } -vector3_reflect :: proc(i, n: Vector3) -> Vector3 { - b := n * 2 * dot(n, i); - return i - b; + +vector4_srgb_to_linear :: proc(col: Vector4) -> Vector4 { + r := math.pow(col.x, 2.2); + g := math.pow(col.y, 2.2); + b := math.pow(col.z, 2.2); + a := col.w; + return {r, g, b, a}; } -vector3_refract :: proc(i, n: Vector3, eta: Float) -> Vector3 { - dv := dot(n, i); - k := 1 - eta*eta - (1 - dv*dv); - a := i * eta; - b := n * eta*dv*math.sqrt(k); - return (a - b) * Float(int(k >= 0)); +vector4_linear_to_srgb :: proc(col: Vector4) -> Vector4 { + a :: 2.51; + b :: 0.03; + c :: 2.43; + d :: 0.59; + e :: 0.14; + + x := col.x; + y := col.y; + z := col.z; + + x = (x * (a * x + b)) / (x * (c * x + d) + e); + y = (y * (a * y + b)) / (y * (c * y + d) + e); + z = (z * (a * z + b)) / (z * (c * z + d) + e); + + x = math.pow(clamp(x, 0, 1), 1.0 / 2.2); + y = math.pow(clamp(y, 0, 1), 1.0 / 2.2); + z = math.pow(clamp(z, 0, 1), 1.0 / 2.2); + + return {x, y, z, col.w}; +} + +vector4_hsl_to_rgb :: proc(h, s, l: Float, a: Float = 1) -> Vector4 { + hue_to_rgb :: proc(p, q, t0: Float) -> Float { + t := math.mod(t0, 1.0); + switch { + case t < 1.0/6.0: return p + (q - p) * 6.0 * t; + case t < 1.0/2.0: return q; + case t < 2.0/3.0: return p + (q - p) * 6.0 * (2.0/3.0 - t); + } + return p; + } + + r, g, b: Float; + if s == 0 { + r = l; + g = l; + b = l; + } else { + q := l < 0.5 ? l * (1+s) : l+s - l*s; + p := 2*l - q; + r = hue_to_rgb(p, q, h + 1.0/3.0); + g = hue_to_rgb(p, q, h); + b = hue_to_rgb(p, q, h - 1.0/3.0); + } + return {r, g, b, a}; +} + +vector4_rgb_to_hsl :: proc(col: Vector4) -> Vector4 { + r := col.x; + g := col.y; + b := col.z; + a := col.w; + v_min := min(r, g, b); + v_max := max(r, g, b); + h, s, l: Float; + h = 0.0; + s = 0.0; + l = (v_min + v_max) * 0.5; + + if v_max != v_min { + d: = v_max - v_min; + s = l > 0.5 ? d / (2.0 - v_max - v_min) : d / (v_max + v_min); + switch { + case v_max == r: + h = (g - b) / d + (g < b ? 6.0 : 0.0); + case v_max == g: + h = (b - r) / d + 2.0; + case v_max == b: + h = (r - g) / d + 4.0; + } + + h *= 1.0/6.0; + } + + return {h, s, l, a}; +} + + + +quaternion_angle_axis :: proc(angle_radians: Float, axis: Vector3) -> Quaternion { + t := angle_radians*0.5; + w := math.cos(t); + v := normalize(axis) * math.sin(t); + return quaternion(w, v.x, v.y, v.z); +} + +quaternion_from_euler_angles :: proc(pitch, yaw, roll: Float) -> Quaternion { + p := quaternion_angle_axis(pitch, {1, 0, 0}); + y := quaternion_angle_axis(yaw, {0, 1, 0}); + r := quaternion_angle_axis(roll, {0, 0, 1}); + return (y * p) * r; +} + +euler_angles_from_quaternion :: proc(q: Quaternion) -> (roll, pitch, yaw: Float) { + // roll (x-axis rotation) + sinr_cosp: Float = 2 * (real(q)*imag(q) + jmag(q)*kmag(q)); + cosr_cosp: Float = 1 - 2 * (imag(q)*imag(q) + jmag(q)*jmag(q)); + roll = Float(math.atan2(sinr_cosp, cosr_cosp)); + + // pitch (y-axis rotation) + sinp: Float = 2 * (real(q)*kmag(q) - kmag(q)*imag(q)); + if abs(sinp) >= 1 { + pitch = Float(math.copy_sign(math.TAU * 0.25, sinp)); + } else { + pitch = Float(math.asin(sinp)); + } + + // yaw (z-axis rotation) + siny_cosp: Float = 2 * (real(q)*kmag(q) + imag(q)*jmag(q)); + cosy_cosp: Float = 1 - 2 * (jmag(q)*jmag(q) + kmag(q)*kmag(q)); + yaw = Float(math.atan2(siny_cosp, cosy_cosp)); + + return; +} + + +quaternion_nlerp :: proc(a, b: Quaternion, t: Float) -> Quaternion { + c := a + (b-a)*quaternion(t, 0, 0, 0); + return normalize(c); +} + + +quaternion_slerp :: proc(x, y: Quaternion, t: Float) -> Quaternion { + EPSILON :: size_of(Float) == 4 ? 1e-7 : 1e-15; + + a, b := x, y; + cos_angle := dot(a, b); + if cos_angle < 0 { + b = -b; + cos_angle = -cos_angle; + } + if cos_angle > 1 - EPSILON { + return a + (b-a)*quaternion(t, 0, 0, 0); + } + + angle := math.acos(cos_angle); + sin_angle := math.sin(angle); + factor_a, factor_b: Quaternion; + factor_a = quaternion(math.sin((1-t) * angle) / sin_angle, 0, 0, 0); + factor_b = quaternion(math.sin(t * angle) / sin_angle, 0, 0, 0); + + return factor_a * a + factor_b * b; +} + + +quaternion_from_matrix4 :: proc(m: Matrix4) -> Quaternion { + four_x_squared_minus_1, four_y_squared_minus_1, + four_z_squared_minus_1, four_w_squared_minus_1, + four_biggest_squared_minus_1: Float; + + /* xyzw */ + /* 0123 */ + biggest_index := 3; + biggest_value, mult: Float; + + four_x_squared_minus_1 = m[0][0] - m[1][1] - m[2][2]; + four_y_squared_minus_1 = m[1][1] - m[0][0] - m[2][2]; + four_z_squared_minus_1 = m[2][2] - m[0][0] - m[1][1]; + four_w_squared_minus_1 = m[0][0] + m[1][1] + m[2][2]; + + four_biggest_squared_minus_1 = four_w_squared_minus_1; + if four_x_squared_minus_1 > four_biggest_squared_minus_1 { + four_biggest_squared_minus_1 = four_x_squared_minus_1; + biggest_index = 0; + } + if four_y_squared_minus_1 > four_biggest_squared_minus_1 { + four_biggest_squared_minus_1 = four_y_squared_minus_1; + biggest_index = 1; + } + if four_z_squared_minus_1 > four_biggest_squared_minus_1 { + four_biggest_squared_minus_1 = four_z_squared_minus_1; + biggest_index = 2; + } + + biggest_value = math.sqrt(four_biggest_squared_minus_1 + 1) * 0.5; + mult = 0.25 / biggest_value; + + + switch biggest_index { + case 0: + return quaternion( + biggest_value, + (m[0][1] + m[1][0]) * mult, + (m[2][0] + m[0][2]) * mult, + (m[1][2] - m[2][1]) * mult, + ); + case 1: + return quaternion( + (m[0][1] + m[1][0]) * mult, + biggest_value, + (m[1][2] + m[2][1]) * mult, + (m[2][0] - m[0][2]) * mult, + ); + case 2: + return quaternion( + (m[2][0] + m[0][2]) * mult, + (m[1][2] + m[2][1]) * mult, + biggest_value, + (m[0][1] - m[1][0]) * mult, + ); + case 3: + return quaternion( + (m[1][2] - m[2][1]) * mult, + (m[2][0] - m[0][2]) * mult, + (m[0][1] - m[1][0]) * mult, + biggest_value, + ); + } + + return 0; +} + + +quaternion_between_two_vector3 :: proc(from, to: Vector3) -> Quaternion { + EPSILON :: size_of(Float) == 4 ? 1e-7 : 1e-15; + + x := normalize(from); + y := normalize(to); + + cos_theta := dot(x, y); + if abs(cos_theta + 1) < 2*EPSILON { + v := vector3_orthogonal(x); + return quaternion(0, v.x, v.y, v.z); + } + v := cross(x, y); + w := cos_theta + 1; + return Quaternion(normalize(quaternion(w, v.x, v.y, v.z))); +} + + +matrix2_inverse_transpose :: proc(m: Matrix2) -> Matrix2 { + c: Matrix2; + d := m[0][0]*m[1][1] - m[1][0]*m[0][1]; + id := 1.0/d; + c[0][0] = +m[1][1] * id; + c[0][1] = -m[0][1] * id; + c[1][0] = -m[1][0] * id; + c[1][1] = +m[0][0] * id; + return c; +} +matrix2_determinant :: proc(m: Matrix2) -> Float { + return m[0][0]*m[1][1] - m[1][0]*m[0][1]; +} + +matrix2_adjoint :: proc(m: Matrix2) -> Matrix2 { + c: Matrix2; + c[0][0] = +m[1][1]; + c[0][1] = -m[1][0]; + c[1][0] = -m[0][1]; + c[1][1] = +m[0][0]; + return c; +} + + +matrix3_from_quaternion :: proc(q: Quaternion) -> Matrix3 { + xx := imag(q) * imag(q); + xy := imag(q) * jmag(q); + xz := imag(q) * kmag(q); + xw := imag(q) * real(q); + yy := jmag(q) * jmag(q); + yz := jmag(q) * kmag(q); + yw := jmag(q) * real(q); + zz := kmag(q) * kmag(q); + zw := kmag(q) * real(q); + + m: Matrix3; + m[0][0] = 1 - 2 * (yy + zz); + m[1][0] = 2 * (xy - zw); + m[2][0] = 2 * (xz + yw); + + m[0][1] = 2 * (xy + zw); + m[1][1] = 1 - 2 * (xx + zz); + m[2][1] = 2 * (yz - xw); + + m[0][2] = 2 * (xz - yw); + m[1][2] = 2 * (yz + xw); + m[2][2] = 1 - 2 * (xx + yy); + + return m; +} + +matrix3_inverse :: proc(m: Matrix3) -> Matrix3 { + return transpose(matrix3_inverse_transpose(m)); +} + + +matrix3_determinant :: proc(m: Matrix3) -> Float { + a := +m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]); + b := -m[1][0] * (m[0][1] * m[2][2] - m[2][1] * m[0][2]); + c := +m[2][0] * (m[0][1] * m[1][2] - m[1][1] * m[0][2]); + return a + b + c; +} + +matrix3_adjoint :: proc(m: Matrix3) -> Matrix3 { + adjoint: Matrix3; + adjoint[0][0] = +(m[1][1] * m[2][2] - m[1][2] * m[2][1]); + adjoint[1][0] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]); + adjoint[2][0] = +(m[0][1] * m[1][2] - m[0][2] * m[1][1]); + adjoint[0][1] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]); + adjoint[1][1] = +(m[0][0] * m[2][2] - m[0][2] * m[2][0]); + adjoint[2][1] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]); + adjoint[0][2] = +(m[1][0] * m[2][1] - m[1][1] * m[2][0]); + adjoint[1][2] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]); + adjoint[2][2] = +(m[0][0] * m[1][1] - m[0][1] * m[1][0]); + return adjoint; +} + +matrix3_inverse_transpose :: proc(m: Matrix3) -> Matrix3 { + inverse_transpose: Matrix3; + + adjoint := matrix3_adjoint(m); + determinant := matrix3_determinant(m); + inv_determinant := 1.0 / determinant; + for i in 0..<3 { + for j in 0..<3 { + inverse_transpose[i][j] = adjoint[i][j] * inv_determinant; + } + } + return inverse_transpose; +} + + +matrix3_scale :: proc(s: Vector3) -> Matrix3 { + m: Matrix3; + m[0][0] = s[0]; + m[1][1] = s[1]; + m[2][2] = s[2]; + return m; +} + +matrix4_from_quaternion :: proc(q: Quaternion) -> Matrix4 { + m := identity(Matrix4); + + xx := imag(q) * imag(q); + xy := imag(q) * jmag(q); + xz := imag(q) * kmag(q); + xw := imag(q) * real(q); + yy := jmag(q) * jmag(q); + yz := jmag(q) * kmag(q); + yw := jmag(q) * real(q); + zz := kmag(q) * kmag(q); + zw := kmag(q) * real(q); + + m[0][0] = 1 - 2 * (yy + zz); + m[1][0] = 2 * (xy - zw); + m[2][0] = 2 * (xz + yw); + + m[0][1] = 2 * (xy + zw); + m[1][1] = 1 - 2 * (xx + zz); + m[2][1] = 2 * (yz - xw); + + m[0][2] = 2 * (xz - yw); + m[1][2] = 2 * (yz + xw); + m[2][2] = 1 - 2 * (xx + yy); + + return m; +} + +matrix4_from_trs :: proc(t: Vector3, r: Quaternion, s: Vector3) -> Matrix4 { + translation := matrix4_translate(t); + rotation := matrix4_from_quaternion(r); + scale := matrix4_scale(s); + return mul(translation, mul(rotation, scale)); +} + + +matrix4_inverse :: proc(m: Matrix4) -> Matrix4 { + return transpose(matrix4_inverse_transpose(m)); +} + + +matrix4_minor :: proc(m: Matrix4, c, r: int) -> Float { + cut_down: Matrix3; + for i in 0..<3 { + col := i < c ? i : i+1; + for j in 0..<3 { + row := j < r ? j : j+1; + cut_down[i][j] = m[col][row]; + } + } + return matrix3_determinant(cut_down); +} + +matrix4_cofactor :: proc(m: Matrix4, c, r: int) -> Float { + sign := (c + r) % 2 == 0 ? Float(1) : Float(-1); + minor := matrix4_minor(m, c, r); + return sign * minor; +} + +matrix4_adjoint :: proc(m: Matrix4) -> Matrix4 { + adjoint: Matrix4; + for i in 0..<4 { + for j in 0..<4 { + adjoint[i][j] = matrix4_cofactor(m, i, j); + } + } + return adjoint; +} + +matrix4_determinant :: proc(m: Matrix4) -> Float { + adjoint := matrix4_adjoint(m); + determinant: Float = 0; + for i in 0..<4 { + determinant += m[i][0] * adjoint[i][0]; + } + return determinant; + +} + +matrix4_inverse_transpose :: proc(m: Matrix4) -> Matrix4 { + adjoint := matrix4_adjoint(m); + determinant: Float = 0; + for i in 0..<4 { + determinant += m[i][0] * adjoint[i][0]; + } + inv_determinant := 1.0 / determinant; + inverse_transpose: Matrix4; + for i in 0..<4 { + for j in 0..<4 { + inverse_transpose[i][j] = adjoint[i][j] * inv_determinant; + } + } + return inverse_transpose; } @@ -261,16 +911,15 @@ matrix4_rotate :: proc(v: Vector3, angle_radians: Float) -> Matrix4 { return rot; } -scale_matrix4 :: matrix4_scale; -matrix4_scale :: proc(m: Matrix4, v: Vector3) -> Matrix4 { - mm := m; - mm[0][0] *= v[0]; - mm[1][1] *= v[1]; - mm[2][2] *= v[2]; - return mm; +matrix4_scale :: proc(v: Vector3) -> Matrix4 { + m: Matrix4; + m[0][0] = v[0]; + m[1][1] = v[1]; + m[2][2] = v[2]; + m[3][3] = 1; + return m; } -look_at :: matrix4_look_at; matrix4_look_at :: proc(eye, centre, up: Vector3) -> Matrix4 { f := normalize(centre - eye); s := normalize(cross(f, up)); @@ -284,7 +933,6 @@ matrix4_look_at :: proc(eye, centre, up: Vector3) -> Matrix4 { } -perspective :: matrix4_perspective; matrix4_perspective :: proc(fovy, aspect, near, far: Float) -> (m: Matrix4) { tan_half_fovy := math.tan(0.5 * fovy); m[0][0] = 1 / (aspect*tan_half_fovy); @@ -308,41 +956,12 @@ matrix_ortho3d :: proc(left, right, bottom, top, near, far: Float) -> (m: Matrix } -axis_angle :: quaternion_angle_axis; -angle_axis :: quaternion_angle_axis; -quaternion_angle_axis :: proc(angle_radians: Float, axis: Vector3) -> Quaternion { - t := angle_radians*0.5; - w := math.cos(t); - v := normalize(axis) * math.sin(t); - return quaternion(w, v.x, v.y, v.z); -} - -euler_angles :: quaternion_from_euler_angles; -quaternion_from_euler_angles :: proc(pitch, yaw, roll: Float) -> Quaternion { - p := quaternion_angle_axis(pitch, {1, 0, 0}); - y := quaternion_angle_axis(yaw, {0, 1, 0}); - r := quaternion_angle_axis(roll, {0, 0, 1}); - return (y * p) * r; -} - -euler_angles_from_quaternion :: proc(q: Quaternion) -> (roll, pitch, yaw: Float) { - // roll (x-axis rotation) - sinr_cosp: Float = 2 * (real(q)*imag(q) + jmag(q)*kmag(q)); - cosr_cosp: Float = 1 - 2 * (imag(q)*imag(q) + jmag(q)*jmag(q)); - roll = Float(math.atan2(sinr_cosp, cosr_cosp)); - - // pitch (y-axis rotation) - sinp: Float = 2 * (real(q)*kmag(q) - kmag(q)*imag(q)); - if abs(sinp) >= 1 { - pitch = Float(math.copy_sign(math.TAU * 0.25, sinp)); - } else { - pitch = Float(math.asin(sinp)); - } - - // yaw (z-axis rotation) - siny_cosp: Float = 2 * (real(q)*kmag(q) + imag(q)*jmag(q)); - cosy_cosp: Float = 1 - 2 * (jmag(q)*jmag(q) + kmag(q)*kmag(q)); - yaw = Float(math.atan2(siny_cosp, cosy_cosp)); - +matrix4_infinite_perspective :: proc(fovy, aspect, near: Float) -> (m: Matrix4) { + tan_half_fovy := math.tan(0.5 * fovy); + m[0][0] = 1 / (aspect*tan_half_fovy); + m[1][1] = 1 / (tan_half_fovy); + m[2][2] = -1; + m[2][3] = -1; + m[3][2] = -2*near; return; }