lib/linear/src/matrix.zig
daab053ee43316e1809a84551d573ddd1e5bf3d2
1 //! Column-major 3x3 and 4x4 matrices over f32.
2 //!
3 //! `cols[j]` holds column j, so a matrix uploads unchanged as a GLSL, Metal,
4 //! or Vulkan `mat3` or `mat4` array. A matrix multiplies column vectors on
5 //! its right. View and projection matrices are right-handed: the camera looks
6 //! down -z with y up, and clip depth runs from 0 at the near plane to 1 at the
7 //! far plane.
8 const std = @import("std");
9 const vector = @import("vector.zig");
10 const Quat = @import("quaternion.zig").Quat;
11
12 const assert = std.debug.assert;
13 const Vec3 = vector.Vec3;
14 const Vec4 = vector.Vec4;
15
16 pub const Mat3 = extern struct {
17 cols: [3]Vec3 = .{
18 .{ .x = 1 },
19 .{ .y = 1 },
20 .{ .z = 1 },
21 },
22
23 pub const identity: Mat3 = .{};
24
25 pub fn fromCols(c0: Vec3, c1: Vec3, c2: Vec3) Mat3 {
26 return .{ .cols = .{ c0, c1, c2 } };
27 }
28
29 pub fn fromScale(s: Vec3) Mat3 {
30 return fromCols(.{ .x = s.x }, .{ .y = s.y }, .{ .z = s.z });
31 }
32
33 /// The rotation matrix of the unit quaternion `q`.
34 pub fn fromQuat(q: Quat) Mat3 {
35 const xx = q.x * q.x;
36 const yy = q.y * q.y;
37 const zz = q.z * q.z;
38 const xy = q.x * q.y;
39 const xz = q.x * q.z;
40 const yz = q.y * q.z;
41 const wx = q.w * q.x;
42 const wy = q.w * q.y;
43 const wz = q.w * q.z;
44 return fromCols(
45 .{ .x = 1 - 2 * (yy + zz), .y = 2 * (xy + wz), .z = 2 * (xz - wy) },
46 .{ .x = 2 * (xy - wz), .y = 1 - 2 * (xx + zz), .z = 2 * (yz + wx) },
47 .{ .x = 2 * (xz + wy), .y = 2 * (yz - wx), .z = 1 - 2 * (xx + yy) },
48 );
49 }
50
51 pub fn row(m: Mat3, index: usize) Vec3 {
52 assert(index < 3);
53 const lanes = [3][3]f32{ m.cols[0].toArray(), m.cols[1].toArray(), m.cols[2].toArray() };
54 return .{ .x = lanes[0][index], .y = lanes[1][index], .z = lanes[2][index] };
55 }
56
57 pub fn transpose(m: Mat3) Mat3 {
58 return fromCols(m.row(0), m.row(1), m.row(2));
59 }
60
61 /// `m v`, summed as `c0 v.x + c1 v.y + c2 v.z`.
62 pub fn mulVec(m: Mat3, v: Vec3) Vec3 {
63 return m.cols[0].scale(v.x).add(m.cols[1].scale(v.y)).add(m.cols[2].scale(v.z));
64 }
65
66 /// `a b`: applying `b`, then `a`.
67 pub fn mul(a: Mat3, b: Mat3) Mat3 {
68 return fromCols(a.mulVec(b.cols[0]), a.mulVec(b.cols[1]), a.mulVec(b.cols[2]));
69 }
70
71 pub fn scale(m: Mat3, s: f32) Mat3 {
72 return fromCols(m.cols[0].scale(s), m.cols[1].scale(s), m.cols[2].scale(s));
73 }
74
75 pub fn determinant(m: Mat3) f32 {
76 return m.cols[0].dot(m.cols[1].cross(m.cols[2]));
77 }
78
79 /// The inverse, or null when the determinant is zero, NaN, or too small
80 /// for its reciprocal to be a finite nonzero f32.
81 pub fn inverse(m: Mat3) ?Mat3 {
82 const inverse_transposed = m.inverseTranspose() orelse return null;
83 return inverse_transposed.transpose();
84 }
85
86 /// The transpose of the inverse, which carries surface normals through
87 /// `m`. Null under the same condition as `inverse`.
88 pub fn inverseTranspose(m: Mat3) ?Mat3 {
89 const c0 = m.cols[1].cross(m.cols[2]);
90 const c1 = m.cols[2].cross(m.cols[0]);
91 const c2 = m.cols[0].cross(m.cols[1]);
92 const inverse_determinant = reciprocal(m.cols[0].dot(c0)) orelse return null;
93 return fromCols(c0, c1, c2).scale(inverse_determinant);
94 }
95 };
96
97 pub const Mat4 = extern struct {
98 cols: [4]Vec4 = .{
99 .{ .x = 1 },
100 .{ .y = 1 },
101 .{ .z = 1 },
102 .{ .w = 1 },
103 },
104
105 pub const identity: Mat4 = .{};
106
107 pub fn fromCols(c0: Vec4, c1: Vec4, c2: Vec4, c3: Vec4) Mat4 {
108 return .{ .cols = .{ c0, c1, c2, c3 } };
109 }
110
111 /// The homogeneous matrix of a linear part and a translation.
112 pub fn fromLinearTranslation(linear: Mat3, translation: Vec3) Mat4 {
113 return fromCols(
114 Vec4.fromVec3(linear.cols[0], 0),
115 Vec4.fromVec3(linear.cols[1], 0),
116 Vec4.fromVec3(linear.cols[2], 0),
117 Vec4.fromVec3(translation, 1),
118 );
119 }
120
121 pub fn row(m: Mat4, index: usize) Vec4 {
122 assert(index < 4);
123 const lanes = [4][4]f32{
124 m.cols[0].toArray(),
125 m.cols[1].toArray(),
126 m.cols[2].toArray(),
127 m.cols[3].toArray(),
128 };
129 return .{ .x = lanes[0][index], .y = lanes[1][index], .z = lanes[2][index], .w = lanes[3][index] };
130 }
131
132 pub fn transpose(m: Mat4) Mat4 {
133 return fromCols(m.row(0), m.row(1), m.row(2), m.row(3));
134 }
135
136 /// `m v`, summed as `c0 v.x + c1 v.y + c2 v.z + c3 v.w`.
137 pub fn mulVec(m: Mat4, v: Vec4) Vec4 {
138 return m.cols[0].scale(v.x)
139 .add(m.cols[1].scale(v.y))
140 .add(m.cols[2].scale(v.z))
141 .add(m.cols[3].scale(v.w));
142 }
143
144 /// `a b`: applying `b`, then `a`.
145 pub fn mul(a: Mat4, b: Mat4) Mat4 {
146 return fromCols(
147 a.mulVec(b.cols[0]),
148 a.mulVec(b.cols[1]),
149 a.mulVec(b.cols[2]),
150 a.mulVec(b.cols[3]),
151 );
152 }
153
154 pub fn scale(m: Mat4, s: f32) Mat4 {
155 return fromCols(m.cols[0].scale(s), m.cols[1].scale(s), m.cols[2].scale(s), m.cols[3].scale(s));
156 }
157
158 pub fn determinant(m: Mat4) f32 {
159 const minors = Minors.of(m);
160 return minors.determinant();
161 }
162
163 /// The inverse, or null when the determinant is zero, NaN, or too small
164 /// for its reciprocal to be a finite nonzero f32.
165 pub fn inverse(m: Mat4) ?Mat4 {
166 const minors = Minors.of(m);
167 const inverse_determinant = reciprocal(minors.determinant()) orelse return null;
168 const a = m.cols;
169 const adjugate = fromCols(
170 .{
171 .x = a[1].y * minors.c5 - a[1].z * minors.c4 + a[1].w * minors.c3,
172 .y = -a[0].y * minors.c5 + a[0].z * minors.c4 - a[0].w * minors.c3,
173 .z = a[3].y * minors.s5 - a[3].z * minors.s4 + a[3].w * minors.s3,
174 .w = -a[2].y * minors.s5 + a[2].z * minors.s4 - a[2].w * minors.s3,
175 },
176 .{
177 .x = -a[1].x * minors.c5 + a[1].z * minors.c2 - a[1].w * minors.c1,
178 .y = a[0].x * minors.c5 - a[0].z * minors.c2 + a[0].w * minors.c1,
179 .z = -a[3].x * minors.s5 + a[3].z * minors.s2 - a[3].w * minors.s1,
180 .w = a[2].x * minors.s5 - a[2].z * minors.s2 + a[2].w * minors.s1,
181 },
182 .{
183 .x = a[1].x * minors.c4 - a[1].y * minors.c2 + a[1].w * minors.c0,
184 .y = -a[0].x * minors.c4 + a[0].y * minors.c2 - a[0].w * minors.c0,
185 .z = a[3].x * minors.s4 - a[3].y * minors.s2 + a[3].w * minors.s0,
186 .w = -a[2].x * minors.s4 + a[2].y * minors.s2 - a[2].w * minors.s0,
187 },
188 .{
189 .x = -a[1].x * minors.c3 + a[1].y * minors.c1 - a[1].z * minors.c0,
190 .y = a[0].x * minors.c3 - a[0].y * minors.c1 + a[0].z * minors.c0,
191 .z = -a[3].x * minors.s3 + a[3].y * minors.s1 - a[3].z * minors.s0,
192 .w = a[2].x * minors.s3 - a[2].y * minors.s1 + a[2].z * minors.s0,
193 },
194 );
195 return adjugate.scale(inverse_determinant);
196 }
197
198 /// The view matrix of a camera at `eye` looking at `target`, or null when
199 /// the two coincide or the view direction is parallel to `up`.
200 pub fn lookAt(eye: Vec3, target: Vec3, up: Vec3) ?Mat4 {
201 const forward = target.sub(eye).normalized(0) orelse return null;
202 const side = forward.cross(up).normalized(0) orelse return null;
203 const camera_up = side.cross(forward);
204 return fromCols(
205 .{ .x = side.x, .y = camera_up.x, .z = -forward.x },
206 .{ .x = side.y, .y = camera_up.y, .z = -forward.y },
207 .{ .x = side.z, .y = camera_up.z, .z = -forward.z },
208 .{ .x = -side.dot(eye), .y = -camera_up.dot(eye), .z = forward.dot(eye), .w = 1 },
209 );
210 }
211
212 /// A perspective projection with vertical field of view `fov_y` radians
213 /// and `aspect` width over height, mapping view depth `-near` to clip
214 /// depth 0 and `-far` to 1.
215 pub fn perspective(fov_y: f32, aspect: f32, near: f32, far: f32) Mat4 {
216 assert(fov_y > 0);
217 assert(fov_y < std.math.pi);
218 assert(aspect > 0);
219 assert(near > 0);
220 assert(far > near);
221 const focal = 1 / @tan(fov_y * 0.5);
222 const depth = 1 / (near - far);
223 return fromCols(
224 .{ .x = focal / aspect },
225 .{ .y = focal },
226 .{ .z = far * depth, .w = -1 },
227 .{ .z = near * far * depth },
228 );
229 }
230
231 /// An orthographic projection of the view-space box between the given
232 /// planes, mapping view depth `-near` to clip depth 0 and `-far` to 1.
233 pub fn orthographic(left: f32, right: f32, bottom: f32, top: f32, near: f32, far: f32) Mat4 {
234 assert(right != left);
235 assert(top != bottom);
236 assert(far != near);
237 const width = 1 / (right - left);
238 const height = 1 / (top - bottom);
239 const depth = 1 / (near - far);
240 return fromCols(
241 .{ .x = 2 * width },
242 .{ .y = 2 * height },
243 .{ .z = depth },
244 .{ .x = -(right + left) * width, .y = -(top + bottom) * height, .z = near * depth, .w = 1 },
245 );
246 }
247 };
248
249 /// The twelve 2x2 minors that the 4x4 determinant and adjugate share: `s`
250 /// from columns 0 and 1, `c` from columns 2 and 3.
251 const Minors = struct {
252 s0: f32,
253 s1: f32,
254 s2: f32,
255 s3: f32,
256 s4: f32,
257 s5: f32,
258 c0: f32,
259 c1: f32,
260 c2: f32,
261 c3: f32,
262 c4: f32,
263 c5: f32,
264
265 fn of(m: Mat4) Minors {
266 const a = m.cols;
267 return .{
268 .s0 = a[0].x * a[1].y - a[1].x * a[0].y,
269 .s1 = a[0].x * a[1].z - a[1].x * a[0].z,
270 .s2 = a[0].x * a[1].w - a[1].x * a[0].w,
271 .s3 = a[0].y * a[1].z - a[1].y * a[0].z,
272 .s4 = a[0].y * a[1].w - a[1].y * a[0].w,
273 .s5 = a[0].z * a[1].w - a[1].z * a[0].w,
274 .c5 = a[2].z * a[3].w - a[3].z * a[2].w,
275 .c4 = a[2].y * a[3].w - a[3].y * a[2].w,
276 .c3 = a[2].y * a[3].z - a[3].y * a[2].z,
277 .c2 = a[2].x * a[3].w - a[3].x * a[2].w,
278 .c1 = a[2].x * a[3].z - a[3].x * a[2].z,
279 .c0 = a[2].x * a[3].y - a[3].x * a[2].y,
280 };
281 }
282
283 fn determinant(minors: Minors) f32 {
284 return minors.s0 * minors.c5 - minors.s1 * minors.c4 + minors.s2 * minors.c3 +
285 minors.s3 * minors.c2 - minors.s4 * minors.c1 + minors.s5 * minors.c0;
286 }
287 };
288
289 /// `1 / value` when it is a finite nonzero f32, else null.
290 fn reciprocal(value: f32) ?f32 {
291 const result = 1 / value;
292 if (!std.math.isFinite(result) or result == 0) return null;
293 return result;
294 }
295
296 comptime {
297 assert(@sizeOf(Mat3) == @sizeOf([9]f32));
298 assert(@alignOf(Mat3) == @alignOf(f32));
299 assert(@sizeOf(Mat4) == @sizeOf([16]f32));
300 assert(@alignOf(Mat4) == @alignOf(f32));
301 }
302
303 const testing = std.testing;
304 const tolerance: f32 = 1e-5;
305
306 fn expectApproxVec3(expected: Vec3, actual: Vec3) !void {
307 try testing.expectApproxEqAbs(expected.x, actual.x, tolerance);
308 try testing.expectApproxEqAbs(expected.y, actual.y, tolerance);
309 try testing.expectApproxEqAbs(expected.z, actual.z, tolerance);
310 }
311
312 fn expectApproxVec4(expected: Vec4, actual: Vec4) !void {
313 try expectApproxVec3(expected.xyz(), actual.xyz());
314 try testing.expectApproxEqAbs(expected.w, actual.w, tolerance);
315 }
316
317 fn expectApproxMat4(expected: Mat4, actual: Mat4) !void {
318 for (expected.cols, actual.cols) |e, a| try expectApproxVec4(e, a);
319 }
320
321 test "identity matrices are the default" {
322 try testing.expectEqual(Mat3.fromScale(Vec3.splat(1)), Mat3{});
323 const v = Vec4.init(1, 2, 3, 4);
324 try testing.expectEqual(v, Mat4.identity.mulVec(v));
325 }
326
327 test "matrices multiply column vectors on the right" {
328 const m = Mat3.fromCols(Vec3.init(1, 2, 3), Vec3.init(4, 5, 6), Vec3.init(7, 8, 10));
329 try testing.expectEqual(Vec3.init(1, 2, 3), m.mulVec(Vec3.init(1, 0, 0)));
330 try testing.expectEqual(Vec3.init(12, 15, 19), m.mulVec(Vec3.init(1, 1, 1)));
331 try testing.expectEqual(Vec3.init(1, 4, 7), m.row(0));
332 try testing.expectEqual(m, m.transpose().transpose());
333 try testing.expectEqual(@as(f32, -3), m.determinant());
334 }
335
336 test "matrix product composes right to left" {
337 const scale = Mat3.fromScale(Vec3.init(2, 3, 4));
338 const turn = Mat3.fromQuat(Quat.fromAxisAngle(Vec3.init(0, 0, 1), std.math.pi / 2.0));
339 const v = Vec3.init(1, 0, 0);
340 try expectApproxVec3(turn.mulVec(scale.mulVec(v)), turn.mul(scale).mulVec(v));
341 try expectApproxVec3(Vec3.init(0, 2, 0), turn.mul(scale).mulVec(v));
342 }
343
344 test "quaternion matrix agrees with quaternion rotation" {
345 const q = Quat.fromAxisAngle(Vec3.init(0.48, 0.6, 0.64), 0.9);
346 const v = Vec3.init(-0.3, 1.2, 0.7);
347 try expectApproxVec3(q.rotate(v), Mat3.fromQuat(q).mulVec(v));
348 }
349
350 test "inverse returns null for singular input" {
351 const flat = Mat3.fromScale(Vec3.init(1, 0, 1));
352 try testing.expectEqual(@as(?Mat3, null), flat.inverse());
353 try testing.expectEqual(@as(?Mat3, null), Mat3.fromScale(Vec3.splat(1e-20)).inverse());
354 try testing.expectEqual(@as(?Mat3, null), Mat3.fromScale(Vec3.splat(std.math.nan(f32))).inverse());
355 try testing.expectEqual(@as(?Mat4, null), Mat4.identity.scale(0).inverse());
356 const repeated = Mat4.fromCols(Vec4.init(1, 2, 3, 4), Vec4.init(1, 2, 3, 4), Vec4.init(0, 0, 1, 0), Vec4.init(0, 0, 0, 1));
357 try testing.expectEqual(@as(?Mat4, null), repeated.inverse());
358 }
359
360 test "inverse undoes an invertible matrix" {
361 const m = Mat3.fromCols(Vec3.init(1, 2, 3), Vec3.init(4, 5, 6), Vec3.init(7, 8, 10));
362 const product = m.mul(m.inverse().?);
363 for (product.cols, Mat3.identity.cols) |a, e| try expectApproxVec3(e, a);
364
365 const n = Mat4.fromCols(
366 Vec4.init(2, 0, 1, 0),
367 Vec4.init(1, 3, 0, 0),
368 Vec4.init(0, 1, 4, 0),
369 Vec4.init(5, -2, 7, 1),
370 );
371 try expectApproxMat4(Mat4.identity, n.mul(n.inverse().?));
372 try expectApproxMat4(Mat4.identity, n.inverse().?.mul(n));
373 try testing.expectApproxEqAbs(@as(f32, 25), n.determinant(), tolerance);
374 }
375
376 test "look-at places the eye at the origin facing down negative z" {
377 const eye = Vec3.init(3, 4, 5);
378 const target = Vec3.init(3, 4, -5);
379 const view = Mat4.lookAt(eye, target, Vec3.init(0, 1, 0)).?;
380 try expectApproxVec4(Vec4.init(0, 0, 0, 1), view.mulVec(Vec4.fromVec3(eye, 1)));
381 try expectApproxVec4(Vec4.init(0, 0, -10, 1), view.mulVec(Vec4.fromVec3(target, 1)));
382 try expectApproxVec4(Vec4.init(0, 1, 0, 0), view.mulVec(Vec4.init(0, 1, 0, 0)));
383 try expectApproxVec4(Vec4.init(1, 0, 0, 0), view.mulVec(Vec4.init(1, 0, 0, 0)));
384 try testing.expectEqual(@as(?Mat4, null), Mat4.lookAt(eye, eye, Vec3.init(0, 1, 0)));
385 try testing.expectEqual(@as(?Mat4, null), Mat4.lookAt(eye, Vec3.init(3, 9, 5), Vec3.init(0, 1, 0)));
386 }
387
388 test "perspective maps near and far planes to depth zero and one" {
389 const projection = Mat4.perspective(std.math.pi / 2.0, 2, 0.5, 100);
390 const near = projection.mulVec(Vec4.init(0, 0, -0.5, 1));
391 const far = projection.mulVec(Vec4.init(0, 0, -100, 1));
392 try testing.expectApproxEqAbs(@as(f32, 0), near.z / near.w, tolerance);
393 try testing.expectApproxEqAbs(@as(f32, 1), far.z / far.w, tolerance);
394 const corner = projection.mulVec(Vec4.init(2, 1, -1, 1));
395 try testing.expectApproxEqAbs(@as(f32, 1), corner.x / corner.w, tolerance);
396 try testing.expectApproxEqAbs(@as(f32, 1), corner.y / corner.w, tolerance);
397 }
398
399 test "orthographic maps its box onto the clip volume" {
400 const projection = Mat4.orthographic(-4, 4, -2, 2, 1, 11);
401 try expectApproxVec4(Vec4.init(-1, -1, 0, 1), projection.mulVec(Vec4.init(-4, -2, -1, 1)));
402 try expectApproxVec4(Vec4.init(1, 1, 1, 1), projection.mulVec(Vec4.init(4, 2, -11, 1)));
403 }