1use super::vec::{Vec2, Vec3, Vec4};
7
8#[derive(Debug, Clone, Copy, PartialEq)]
13#[non_exhaustive]
14pub struct Mat3 {
15 pub m: [[f64; 3]; 3],
17}
18
19impl Mat3 {
20 pub const IDENTITY: Self = Self {
21 m: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
22 };
23
24 pub const ZERO: Self = Self { m: [[0.0; 3]; 3] };
25
26 pub const fn new(m: [[f64; 3]; 3]) -> Self {
27 Self { m }
28 }
29
30 pub fn from_cols(c0: Vec3, c1: Vec3, c2: Vec3) -> Self {
31 Self {
32 m: [[c0.x, c1.x, c2.x], [c0.y, c1.y, c2.y], [c0.z, c1.z, c2.z]],
33 }
34 }
35
36 pub fn transform_point(&self, v: Vec2) -> Vec2 {
38 let x = self.m[0][0] * v.x + self.m[0][1] * v.y + self.m[0][2];
39 let y = self.m[1][0] * v.x + self.m[1][1] * v.y + self.m[1][2];
40 let w = self.m[2][0] * v.x + self.m[2][1] * v.y + self.m[2][2];
41 if w.abs() < f64::EPSILON {
42 Vec2::new(x, y)
43 } else {
44 Vec2::new(x / w, y / w)
45 }
46 }
47
48 pub fn transform_dir(&self, v: Vec2) -> Vec2 {
50 let x = self.m[0][0] * v.x + self.m[0][1] * v.y;
51 let y = self.m[1][0] * v.x + self.m[1][1] * v.y;
52 Vec2::new(x, y)
53 }
54
55 pub fn mul(&self, rhs: &Self) -> Self {
57 let mut out = Self::ZERO;
58 for i in 0..3 {
59 for j in 0..3 {
60 for k in 0..3 {
61 out.m[i][j] += self.m[i][k] * rhs.m[k][j];
62 }
63 }
64 }
65 out
66 }
67
68 pub fn transpose(&self) -> Self {
69 Self::new([
70 [self.m[0][0], self.m[1][0], self.m[2][0]],
71 [self.m[0][1], self.m[1][1], self.m[2][1]],
72 [self.m[0][2], self.m[1][2], self.m[2][2]],
73 ])
74 }
75
76 pub fn determinant(&self) -> f64 {
77 self.m[0][0] * (self.m[1][1] * self.m[2][2] - self.m[1][2] * self.m[2][1])
78 - self.m[0][1] * (self.m[1][0] * self.m[2][2] - self.m[1][2] * self.m[2][0])
79 + self.m[0][2] * (self.m[1][0] * self.m[2][1] - self.m[1][1] * self.m[2][0])
80 }
81
82 pub fn inverse(&self) -> Option<Self> {
84 let det = self.determinant();
85 if det.abs() < f64::EPSILON {
86 return None;
87 }
88 let inv_det = 1.0 / det;
89
90 Some(Self::new([
91 [
92 (self.m[1][1] * self.m[2][2] - self.m[1][2] * self.m[2][1]) * inv_det,
93 (self.m[0][2] * self.m[2][1] - self.m[0][1] * self.m[2][2]) * inv_det,
94 (self.m[0][1] * self.m[1][2] - self.m[0][2] * self.m[1][1]) * inv_det,
95 ],
96 [
97 (self.m[1][2] * self.m[2][0] - self.m[1][0] * self.m[2][2]) * inv_det,
98 (self.m[0][0] * self.m[2][2] - self.m[0][2] * self.m[2][0]) * inv_det,
99 (self.m[0][2] * self.m[1][0] - self.m[0][0] * self.m[1][2]) * inv_det,
100 ],
101 [
102 (self.m[1][0] * self.m[2][1] - self.m[1][1] * self.m[2][0]) * inv_det,
103 (self.m[0][1] * self.m[2][0] - self.m[0][0] * self.m[2][1]) * inv_det,
104 (self.m[0][0] * self.m[1][1] - self.m[0][1] * self.m[1][0]) * inv_det,
105 ],
106 ]))
107 }
108}
109
110#[derive(Debug, Clone, Copy, PartialEq)]
115#[non_exhaustive]
116pub struct Mat4 {
117 pub m: [[f64; 4]; 4],
119}
120
121impl Mat4 {
122 pub const IDENTITY: Self = Self {
123 m: [
124 [1.0, 0.0, 0.0, 0.0],
125 [0.0, 1.0, 0.0, 0.0],
126 [0.0, 0.0, 1.0, 0.0],
127 [0.0, 0.0, 0.0, 1.0],
128 ],
129 };
130
131 pub const ZERO: Self = Self { m: [[0.0; 4]; 4] };
132
133 pub const fn new(m: [[f64; 4]; 4]) -> Self {
134 Self { m }
135 }
136
137 pub fn transform(&self, v: Vec4) -> Vec4 {
139 Vec4::new(
140 self.m[0][0] * v.x + self.m[0][1] * v.y + self.m[0][2] * v.z + self.m[0][3] * v.w,
141 self.m[1][0] * v.x + self.m[1][1] * v.y + self.m[1][2] * v.z + self.m[1][3] * v.w,
142 self.m[2][0] * v.x + self.m[2][1] * v.y + self.m[2][2] * v.z + self.m[2][3] * v.w,
143 self.m[3][0] * v.x + self.m[3][1] * v.y + self.m[3][2] * v.z + self.m[3][3] * v.w,
144 )
145 }
146
147 pub fn transform_point(&self, v: Vec3) -> Vec3 {
149 self.transform(Vec4::point(v.x, v.y, v.z)).to_vec3()
150 }
151
152 pub fn transform_dir(&self, v: Vec3) -> Vec3 {
154 self.transform(Vec4::direction(v.x, v.y, v.z)).xyz()
155 }
156
157 pub fn mul(&self, rhs: &Self) -> Self {
159 let mut out = Self::ZERO;
160 for i in 0..4 {
161 for j in 0..4 {
162 for k in 0..4 {
163 out.m[i][j] += self.m[i][k] * rhs.m[k][j];
164 }
165 }
166 }
167 out
168 }
169
170 pub fn transpose(&self) -> Self {
171 let mut out = Self::ZERO;
172 for i in 0..4 {
173 for j in 0..4 {
174 out.m[i][j] = self.m[j][i];
175 }
176 }
177 out
178 }
179
180 pub fn to_mat3(&self) -> Mat3 {
182 Mat3::new([
183 [self.m[0][0], self.m[0][1], self.m[0][2]],
184 [self.m[1][0], self.m[1][1], self.m[1][2]],
185 [self.m[2][0], self.m[2][1], self.m[2][2]],
186 ])
187 }
188}
189
190#[cfg(test)]
195mod tests {
196 use super::*;
197
198 #[test]
199 fn mat3_identity_transform() {
200 let p = Vec2::new(5.0, 7.0);
201 let result = Mat3::IDENTITY.transform_point(p);
202 assert!((result.x - 5.0).abs() < f64::EPSILON);
203 assert!((result.y - 7.0).abs() < f64::EPSILON);
204 }
205
206 #[test]
207 fn mat3_multiply_identity() {
208 let m = Mat3::new([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]);
209 let result = m.mul(&Mat3::IDENTITY);
210 assert_eq!(result, m);
211 }
212
213 #[test]
214 fn mat3_determinant() {
215 let m = Mat3::new([[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]]);
216 let det = m.determinant();
217 assert!((det - 1.0).abs() < 1e-10);
218 }
219
220 #[test]
221 fn mat3_inverse_roundtrip() {
222 let m = Mat3::new([[1.0, 2.0, 3.0], [0.0, 1.0, 4.0], [5.0, 6.0, 0.0]]);
223 if let Some(inv) = m.inverse() {
224 let product = m.mul(&inv);
225 for i in 0..3 {
226 for j in 0..3 {
227 let expected = if i == j { 1.0 } else { 0.0 };
228 assert!(
229 (product.m[i][j] - expected).abs() < 1e-10,
230 "m[{i}][{j}] = {} (expected {expected})",
231 product.m[i][j]
232 );
233 }
234 }
235 }
236 }
237
238 #[test]
239 fn mat3_singular_no_inverse() {
240 let m = Mat3::new([[1.0, 2.0, 3.0], [2.0, 4.0, 6.0], [1.0, 1.0, 1.0]]);
241 assert!(m.inverse().is_none());
242 }
243
244 #[test]
245 fn mat4_identity_transform() {
246 let p = Vec3::new(1.0, 2.0, 3.0);
247 let result = Mat4::IDENTITY.transform_point(p);
248 assert!((result.x - 1.0).abs() < f64::EPSILON);
249 assert!((result.y - 2.0).abs() < f64::EPSILON);
250 assert!((result.z - 3.0).abs() < f64::EPSILON);
251 }
252
253 #[test]
254 fn mat4_multiply_identity() {
255 let result = Mat4::IDENTITY.mul(&Mat4::IDENTITY);
256 assert_eq!(result, Mat4::IDENTITY);
257 }
258
259 #[test]
260 fn mat4_transpose() {
261 let m = Mat4::new([
262 [1.0, 2.0, 3.0, 4.0],
263 [5.0, 6.0, 7.0, 8.0],
264 [9.0, 10.0, 11.0, 12.0],
265 [13.0, 14.0, 15.0, 16.0],
266 ]);
267 let t = m.transpose();
268 assert!((t.m[0][1] - 5.0).abs() < f64::EPSILON);
269 assert!((t.m[1][0] - 2.0).abs() < f64::EPSILON);
270 }
271}