math.zig (19617B)
1 const std = @import("std"); 2 3 pub const Vec4 = @Vector(4, f32); 4 pub const Mat4 = [4]Vec4; 5 6 pub fn dot(u: Vec4, v: Vec4) f32 { 7 return @reduce(.Add, u * v); 8 } 9 10 // Adapted from https://geometrian.com/programming/tutorials/cross-product/index.php 11 pub fn vecMul(u: Vec4, v: Vec4) Vec4 { 12 const mask = @Vector(4, i32){ 1, 2, 0, 3 }; 13 const t1 = u * @shuffle(f32, v, undefined, mask); 14 const t2 = v * @shuffle(f32, u, undefined, mask); 15 const t3 = t1 - t2; 16 return @shuffle(f32, t3, undefined, mask); 17 } 18 19 pub fn scale(u: Vec4, s: f32) Vec4 { 20 return u * @as(Vec4, @splat(s)); 21 } 22 23 pub fn norm(u: Vec4) f32 { 24 return @sqrt(dot(u, u)); 25 } 26 27 pub fn normalize(u: Vec4) Vec4 { 28 return u / @as(Vec4, @splat(norm(u))); 29 } 30 31 pub fn matMul(m: Mat4, n: Mat4) Mat4 { 32 var result: Mat4 = undefined; 33 inline for (0..4) |row| { 34 const x = @as(Vec4, @splat(m[row][0])); 35 const y = @as(Vec4, @splat(m[row][1])); 36 const z = @as(Vec4, @splat(m[row][2])); 37 const w = @as(Vec4, @splat(m[row][3])); 38 const xy = @mulAdd(Vec4, x, n[0], y * n[1]); 39 const zw = @mulAdd(Vec4, z, n[2], w * n[3]); 40 result[row] = xy + zw; 41 } 42 return result; 43 } 44 45 pub fn matMulVec(m: Mat4, u: Vec4) Vec4 { 46 return Vec4{ 47 dot(m[0], u), 48 dot(m[1], u), 49 dot(m[2], u), 50 dot(m[3], u), 51 }; 52 } 53 54 test "matMulVec" { 55 const u = matMulVec( 56 .{ 57 .{ 2, 0, 0, 1 }, 58 .{ 0, 3, 0, 1 }, 59 .{ 0, 0, 4, 1 }, 60 .{ 0, 0, 1, 5 }, 61 }, 62 .{ 1, 1, 1, 1 }, 63 ); 64 try std.testing.expectApproxEqAbs(3.0, u[0], 0.000001); 65 try std.testing.expectApproxEqAbs(4.0, u[1], 0.000001); 66 try std.testing.expectApproxEqAbs(5.0, u[2], 0.000001); 67 try std.testing.expectApproxEqAbs(6.0, u[3], 0.000001); 68 } 69 70 pub fn project(m: Mat4, u: Vec4) Vec4 { 71 const v = matMulVec(m, u); 72 return scale(v, 1.0 / v[3]); 73 } 74 75 pub fn lookToLH(eye: Vec4, dir: Vec4, up: Vec4) Mat4 { 76 const f = normalize(dir); 77 const s = normalize(vecMul(up, f)); 78 const u = normalize(vecMul(f, s)); 79 var result = Mat4{ s, u, f, .{ 0, 0, 0, 1 } }; 80 result[0][3] = -dot(s, eye); 81 result[1][3] = -dot(u, eye); 82 result[2][3] = -dot(f, eye); 83 return result; 84 } 85 86 test "lookToLH" { 87 const eye = Vec4{ 5.0, 6.0, 7.0, 1.0 }; 88 const dir = Vec4{ 0.5, 1.0, 0.0, 0.0 }; 89 const up_ref = Vec4{ 0.0, 1.0, 0.0, 0.0 }; 90 const view = lookToLH(eye, dir, up_ref); 91 92 const forward = normalize(dir); 93 const right = normalize(vecMul(up_ref, forward)); 94 const up = normalize(vecMul(forward, right)); 95 96 const eye_cs = matMulVec(view, eye); 97 try std.testing.expectApproxEqAbs(0.0, eye_cs[0], 0.000001); 98 try std.testing.expectApproxEqAbs(0.0, eye_cs[1], 0.000001); 99 try std.testing.expectApproxEqAbs(0.0, eye_cs[2], 0.000001); 100 try std.testing.expectApproxEqAbs(1.0, eye_cs[3], 0.000001); 101 102 const forward_cs = matMulVec(view, eye + forward); 103 try std.testing.expectApproxEqAbs(0.0, forward_cs[0], 0.000001); 104 try std.testing.expectApproxEqAbs(0.0, forward_cs[1], 0.000001); 105 try std.testing.expectApproxEqAbs(1.0, forward_cs[2], 0.000001); 106 try std.testing.expectApproxEqAbs(1.0, forward_cs[3], 0.000001); 107 const backward_cs = matMulVec(view, eye - forward); 108 try std.testing.expectApproxEqAbs(0.0, backward_cs[0], 0.000001); 109 try std.testing.expectApproxEqAbs(0.0, backward_cs[1], 0.000001); 110 try std.testing.expectApproxEqAbs(-1.0, backward_cs[2], 0.000001); 111 try std.testing.expectApproxEqAbs(1.0, backward_cs[3], 0.000001); 112 113 const right_cs = matMulVec(view, eye + right); 114 try std.testing.expectApproxEqAbs(1.0, right_cs[0], 0.000001); 115 try std.testing.expectApproxEqAbs(0.0, right_cs[1], 0.000001); 116 try std.testing.expectApproxEqAbs(0.0, right_cs[2], 0.000001); 117 try std.testing.expectApproxEqAbs(1.0, right_cs[3], 0.000001); 118 const left_cs = matMulVec(view, eye - right); 119 try std.testing.expectApproxEqAbs(-1.0, left_cs[0], 0.000001); 120 try std.testing.expectApproxEqAbs(0.0, left_cs[1], 0.000001); 121 try std.testing.expectApproxEqAbs(0.0, left_cs[2], 0.000001); 122 try std.testing.expectApproxEqAbs(1.0, left_cs[3], 0.000001); 123 124 const up_cs = matMulVec(view, eye + up); 125 try std.testing.expectApproxEqAbs(0.0, up_cs[0], 0.000001); 126 try std.testing.expectApproxEqAbs(1.0, up_cs[1], 0.000001); 127 try std.testing.expectApproxEqAbs(0.0, up_cs[2], 0.000001); 128 try std.testing.expectApproxEqAbs(1.0, up_cs[3], 0.000001); 129 const down_cs = matMulVec(view, eye - up); 130 try std.testing.expectApproxEqAbs(0.0, down_cs[0], 0.000001); 131 try std.testing.expectApproxEqAbs(-1.0, down_cs[1], 0.000001); 132 try std.testing.expectApproxEqAbs(0.0, down_cs[2], 0.000001); 133 try std.testing.expectApproxEqAbs(1.0, down_cs[3], 0.000001); 134 } 135 136 pub fn invLookToLH(eye: Vec4, dir: Vec4, up: Vec4) Mat4 { 137 const f = normalize(dir); 138 const s = normalize(vecMul(up, f)); 139 const u = normalize(vecMul(f, s)); 140 return .{ 141 .{ s[0], u[0], f[0], eye[0] }, 142 .{ s[1], u[1], f[1], eye[1] }, 143 .{ s[2], u[2], f[2], eye[2] }, 144 .{ 0, 0, 0, 1 }, 145 }; 146 } 147 148 test "invLookToLH" { 149 const eye = Vec4{ 5.0, 6.0, 7.0, 1.0 }; 150 const dir = Vec4{ 0.5, 1.0, 0.0, 0.0 }; 151 const up = Vec4{ 0.0, 1.0, 0.0, 0.0 }; 152 const view = lookToLH(eye, dir, up); 153 const inv_view = invLookToLH(eye, dir, up); 154 const result = matMul(view, inv_view); 155 156 try std.testing.expectApproxEqAbs(1.0, result[0][0], 0.000001); 157 try std.testing.expectApproxEqAbs(0.0, result[0][1], 0.000001); 158 try std.testing.expectApproxEqAbs(0.0, result[0][2], 0.000001); 159 try std.testing.expectApproxEqAbs(0.0, result[0][3], 0.000001); 160 161 try std.testing.expectApproxEqAbs(0.0, result[1][0], 0.000001); 162 try std.testing.expectApproxEqAbs(1.0, result[1][1], 0.000001); 163 try std.testing.expectApproxEqAbs(0.0, result[1][2], 0.000001); 164 try std.testing.expectApproxEqAbs(0.0, result[1][3], 0.000001); 165 166 try std.testing.expectApproxEqAbs(0.0, result[2][0], 0.000001); 167 try std.testing.expectApproxEqAbs(0.0, result[2][1], 0.000001); 168 try std.testing.expectApproxEqAbs(1.0, result[2][2], 0.000001); 169 try std.testing.expectApproxEqAbs(0.0, result[2][3], 0.000001); 170 171 try std.testing.expectApproxEqAbs(0.0, result[3][0], 0.000001); 172 try std.testing.expectApproxEqAbs(0.0, result[3][1], 0.000001); 173 try std.testing.expectApproxEqAbs(0.0, result[3][2], 0.000001); 174 try std.testing.expectApproxEqAbs(1.0, result[3][3], 0.000001); 175 } 176 177 pub fn lookAtLH(eye: Vec4, target: Vec4, up: Vec4) Mat4 { 178 return lookToLH(eye, target - eye, up); 179 } 180 181 pub fn perspectiveLH(fovy: f32, aspect: f32, near: f32, far: f32) Mat4 { 182 // we'll split the formation of our perspective matrix into 5 steps: 183 // 1. translate the frustum apex to (0, 0, 0) 184 // 2. map depth to (NDC_Z_MIN, NDC_Z_MAX) or (NDC_Z_MAX, NDC_Z_MIN) 185 // 3. project onto the near plane 186 // 4. map x to (NDC_X_MIN, NDC_X_MAX) and y to (NDC_Y_MIN, NDC_Y_MAX) 187 // 5. combine transformations 188 // 189 // we can skip step 1 since our parameterization can only describe a 190 // frustum with its apex at the origin. 191 // 192 // for step 2, we'll use NDC_Z_MIN=0 and NDC_Z_MAX=1. this works for 193 // metal, directx, and wgpu. for opengl and vulkan, you'll either need to 194 // change your configuration or re-calculate the coefficients used in the 195 // non-linear depth mapping function with NDC_Z_MIN=-1 and NDC_Z_MAX=1. we 196 // will also choose to map the near plane to NDC_Z_MAX and the far plane to 197 // NDC_Z_MIN. this leads to a nicer distribution of values within the 198 // z-buffer if you're dealing with wide ranges of depths. you'll need to 199 // update your configuration or re-calculate the coefficients using a 200 // different mapping. 201 // 202 // for step 3, you could project onto a different plane. i've heard of 203 // people doing this but never done it myself. i'm not aware of what the 204 // advantages or disadvantages would be. 205 // 206 // for step 4, NDC_X/Y_MIN=-1 and NDC_X/Y_MAX=1. i'm not aware of any APIs 207 // that use different values, but if you use one you'll need to recalculate 208 // the depth mapping coefficients. 209 // 210 // in all of the steps, i've assumed a left-handed camera coordinate system 211 // with +z pointing into the screen. to adapt this approach to a 212 // right-handed camera coordinate system, you'll need to be careful about 213 // minus signs anywhere z is used. 214 // 215 // 2. non-linear depth mapping 216 // 217 // d(z) = (a/z) + b 218 // d(z) = (a/z) + b*(z/z) 219 // d(z) = (a + b*z) / z 220 // 221 // as a matrix: 222 // | 1 0 0 0 | 223 // | 0 1 0 0 | 224 // | 0 0 b a | 225 // | 0 0 1 0 | 226 // 227 // this matrix takes advantage of the homogenous divide performed 228 // automatically by most graphics pipelines by setting w=z. 229 // 230 // we want to use a reversed z-buffer so: 231 // d(near) = 1 d(far) = 0 232 // 1 = (a/near)+b 0 = (a/far)+b 233 // 1 = (a/near)-(a/far) 234 // 1 = a(1/near-1/far) 235 // a = 1/(1/near-1/far) 236 // a = near*far/(far-near) 237 // b = -(near*far/(far-near))/far 238 // b = near/(near-far) 239 // 240 // we'll split the formation of our "perspective" matrix into 5 steps: 241 // 1. translate the frustum apex to (0, 0, 0) 242 // 2. map depth to (NDC_Z_MIN, NDC_Z_MAX) or (NDC_Z_MAX, NDC_Z_MIN) 243 // 3. project onto the near plane 244 // 4. map x to (NDC_X_MIN, NDC_X_MAX) and y to (NDC_Y_MIN, NDC_Y_MAX) 245 // 5. combine transformations 246 // 247 // we can skip step 1 since our parameterization can only describe a 248 // frustum with its apex at the origin. 249 // 250 // for step 2, we'll use NDC_Z_MIN=0 and NDC_Z_MAX=1. this works for 251 // metal, directx, and wgpu. for opengl and vulkan, you'll either need to 252 // change your configuration or re-calculate the coefficients used in the 253 // non-linear depth mapping function with NDC_Z_MIN=-1 and NDC_Z_MAX=1. we 254 // will also choose to map the near plane to NDC_Z_MAX and the far plane to 255 // NDC_Z_MIN. this leads to a nicer distribution of values within the 256 // z-buffer if you're dealing with wide ranges of depths. you'll need to 257 // update your configuration or re-calculate the coefficients using a 258 // different mapping. 259 // 260 // for step 3, you could project onto a different plane. i've heard of 261 // people doing this but never done it myself. i'm not aware of what the 262 // advantages or disadvantages would be. 263 // 264 // for step 4, NDC_X/Y_MIN=-1 and NDC_X/Y_MAX=1. i'm not aware of any APIs 265 // that use different values, but if you use one you'll need to recalculate 266 // the depth mapping coefficients. 267 // 268 // in all of the steps, i've assumed a left-handed camera coordinate system 269 // with +z pointing into the screen. to adapt this approach to a 270 // right-handed camera coordinate system, you'll need to be careful about 271 // minus signs anywhere z is used. 272 // 273 // 2. non-linear depth mapping 274 // 275 // d(z) = (a/z) + b 276 // d(z) = (a/z) + b*(z/z) 277 // d(z) = (a + b*z) / z 278 // 279 // as a matrix: 280 // | 1 0 0 0 | 281 // | 0 1 0 0 | 282 // | 0 0 b a | 283 // | 0 0 1 0 | 284 // 285 // this matrix takes advantage of the homogenous divide performed 286 // automatically by most graphics pipelines by setting w=z. 287 // 288 // we want to use a reversed z-buffer so: 289 // d(near) = 1 d(far) = 0 290 // 1 = (a/near)+b 0 = (a/far)+b 291 // 1 = (a/near)-(a/far) 292 // 1 = a(1/near-1/far) 293 // a = 1/(1/near-1/far) 294 // a = near*far/(far-near) 295 // b = -(near*far/(far-near))/far 296 // b = near/(near-far) 297 // 298 // 3. near plane projection: 299 // ^ y 300 // ^ y' | 301 // | | 302 // ------------------------------>| 303 // camera near z 304 // 305 // this projection maps (x, y, z) to (x', y', near). the diagram is 306 // intended to show that the values of y and y' are related by similar 307 // right triangles. this means that 308 // y/z = y'/near 309 // y' = (y*near)/z 310 // 311 // as a matrix: 312 // | near 0 0 0 | 313 // | 0 near 0 0 | 314 // | 0 0 1 0 | 315 // | 0 0 1 0 | 316 // 317 // we utilize the homogenous divide here too, just like we did for our 318 // depth mapping. 319 // 320 // 4. scale x and y 321 // 322 // h = near * tan(fovy/2) 323 // w = near * tan(fovy/2) * aspect 324 // 325 // x' = x / w 326 // y' = y / h 327 // 328 // as a matrix: 329 // | 1/w 0 0 0 | 330 // | 0 1/h 0 0 | 331 // | 0 0 1 0 | 332 // | 0 0 0 1 | 333 // 334 // 5. combine transformations 335 // 336 // combining our previous steps gives us the matrix: 337 // | near/w 0 0 0 | 338 // | 0 near/h 0 0 | 339 // | 0 0 b a | 340 // | 0 0 1 0 | 341 // 342 // and now we can simplify the x and y terms: 343 // near/w = 1 / (tan(fovy/2) * aspect) 344 // near/h = 1 / tan(fovy/2) 345 // 346 // 3. near plane projection: 347 // ^ y 348 // ^ y' | 349 // | | 350 // ------------------------------>| 351 // camera near z 352 // 353 // this projection maps (x, y, z) to (x', y', near). the diagram is 354 // intended to show that the values of y and y' are related by similar 355 // right triangles. this means that 356 // y/z = y'/near 357 // y' = (y*near)/z 358 // 359 // as a matrix: 360 // | near 0 0 0 | 361 // | 0 near 0 0 | 362 // | 0 0 1 0 | 363 // | 0 0 1 0 | 364 // 365 // we utilize the homogenous divide here too, just like we did for our 366 // depth mapping. 367 // 368 // 4. scale x and y 369 // 370 // h = near * tan(fovy/2) 371 // w = near * tan(fovy/2) * aspect 372 // 373 // x' = x / w 374 // y' = y / h 375 // 376 // as a matrix: 377 // | 1/w 0 0 0 | 378 // | 0 1/h 0 0 | 379 // | 0 0 1 0 | 380 // | 0 0 0 1 | 381 // 382 // 5. combine transformations 383 // 384 // combining our previous steps gives us the matrix: 385 // | near/w 0 0 0 | 386 // | 0 near/h 0 0 | 387 // | 0 0 b a | 388 // | 0 0 1 0 | 389 // 390 // and now we can simplify the x and y terms: 391 // near/w = 1 / (tan(fovy/2) * aspect) 392 // near/h = 1 / tan(fovy/2) 393 // 394 // the implementation is pretty simple: 395 const a = near * far / (far - near); 396 const b = near / (near - far); 397 const tan_half_fov = @tan(std.math.degreesToRadians(fovy * 0.5)); 398 const x = 1.0 / (aspect * tan_half_fov); 399 const y = 1.0 / tan_half_fov; 400 return Mat4{ 401 .{ x, 0, 0, 0 }, 402 .{ 0, y, 0, 0 }, 403 .{ 0, 0, b, a }, 404 .{ 0, 0, 1, 0 }, 405 }; 406 } 407 408 test "perspectiveLH" { 409 const fovy = 60.0; 410 const aspect = 16.0 / 9.0; 411 const near = 0.1; 412 const far = 10.0; 413 const proj = perspectiveLH(fovy, aspect, near, far); 414 const tan_half_fov = @tan(std.math.degreesToRadians(fovy * 0.5)); 415 416 const near_h = tan_half_fov * near; 417 const near_w = near_h * aspect; 418 const near_c = project(proj, Vec4{ 0.0, 0.0, near, 1.0 }); 419 try std.testing.expectApproxEqAbs(0.0, near_c[0], 0.000001); 420 try std.testing.expectApproxEqAbs(0.0, near_c[1], 0.000001); 421 try std.testing.expectApproxEqAbs(1.0, near_c[2], 0.000001); 422 try std.testing.expectApproxEqAbs(1.0, near_c[3], 0.000001); 423 const near_tr = project(proj, Vec4{ near_w, near_h, near, 1.0 }); 424 try std.testing.expectApproxEqAbs(1.0, near_tr[0], 0.000001); 425 try std.testing.expectApproxEqAbs(1.0, near_tr[1], 0.000001); 426 try std.testing.expectApproxEqAbs(1.0, near_tr[2], 0.000001); 427 try std.testing.expectApproxEqAbs(1.0, near_tr[3], 0.000001); 428 const near_bl = project(proj, Vec4{ -near_w, -near_h, near, 1.0 }); 429 try std.testing.expectApproxEqAbs(-1.0, near_bl[0], 0.000001); 430 try std.testing.expectApproxEqAbs(-1.0, near_bl[1], 0.000001); 431 try std.testing.expectApproxEqAbs(1.0, near_bl[2], 0.000001); 432 try std.testing.expectApproxEqAbs(1.0, near_bl[3], 0.000001); 433 434 const far_h = tan_half_fov * far; 435 const far_w = far_h * aspect; 436 const far_c = project(proj, Vec4{ 0.0, 0.0, far, 1.0 }); 437 try std.testing.expectApproxEqAbs(0.0, far_c[0], 0.000001); 438 try std.testing.expectApproxEqAbs(0.0, far_c[1], 0.000001); 439 try std.testing.expectApproxEqAbs(0.0, far_c[2], 0.000001); 440 try std.testing.expectApproxEqAbs(1.0, far_c[3], 0.000001); 441 const far_tr = project(proj, Vec4{ far_w, far_h, far, 1.0 }); 442 try std.testing.expectApproxEqAbs(1.0, far_tr[0], 0.000001); 443 try std.testing.expectApproxEqAbs(1.0, far_tr[1], 0.000001); 444 try std.testing.expectApproxEqAbs(0.0, far_tr[2], 0.000001); 445 try std.testing.expectApproxEqAbs(1.0, far_tr[3], 0.000001); 446 const far_bl = project(proj, Vec4{ -far_w, -far_h, far, 1.0 }); 447 try std.testing.expectApproxEqAbs(-1.0, far_bl[0], 0.000001); 448 try std.testing.expectApproxEqAbs(-1.0, far_bl[1], 0.000001); 449 try std.testing.expectApproxEqAbs(0.0, far_tr[2], 0.000001); 450 try std.testing.expectApproxEqAbs(1.0, far_tr[3], 0.000001); 451 } 452 453 pub fn invPerspectiveLH(fovy: f32, aspect: f32, near: f32, far: f32) Mat4 { 454 const inv_a = (far - near) / (near * far); 455 const b = near / (near - far); 456 const tan_half_fov = @tan(std.math.degreesToRadians(fovy * 0.5)); 457 const inv_x = aspect * tan_half_fov; 458 const inv_y = tan_half_fov; 459 return Mat4{ 460 .{ inv_x, 0, 0, 0 }, 461 .{ 0, inv_y, 0, 0 }, 462 .{ 0, 0, 0, 1 }, 463 .{ 0, 0, inv_a, -inv_a * b }, 464 }; 465 } 466 467 test "invPerspectiveLH" { 468 const fovy = 60.0; 469 const aspect = 16.0 / 9.0; 470 const near = 0.1; 471 const far = 10.0; 472 const proj = perspectiveLH(fovy, aspect, near, far); 473 const inv_proj = invPerspectiveLH(fovy, aspect, near, far); 474 const result = matMul(proj, inv_proj); 475 476 try std.testing.expectApproxEqAbs(1.0, result[0][0], 0.000001); 477 try std.testing.expectApproxEqAbs(0.0, result[0][1], 0.000001); 478 try std.testing.expectApproxEqAbs(0.0, result[0][2], 0.000001); 479 try std.testing.expectApproxEqAbs(0.0, result[0][3], 0.000001); 480 481 try std.testing.expectApproxEqAbs(0.0, result[1][0], 0.000001); 482 try std.testing.expectApproxEqAbs(1.0, result[1][1], 0.000001); 483 try std.testing.expectApproxEqAbs(0.0, result[1][2], 0.000001); 484 try std.testing.expectApproxEqAbs(0.0, result[1][3], 0.000001); 485 486 try std.testing.expectApproxEqAbs(0.0, result[2][0], 0.000001); 487 try std.testing.expectApproxEqAbs(0.0, result[2][1], 0.000001); 488 try std.testing.expectApproxEqAbs(1.0, result[2][2], 0.000001); 489 try std.testing.expectApproxEqAbs(0.0, result[2][3], 0.000001); 490 491 try std.testing.expectApproxEqAbs(0.0, result[3][0], 0.000001); 492 try std.testing.expectApproxEqAbs(0.0, result[3][1], 0.000001); 493 try std.testing.expectApproxEqAbs(0.0, result[3][2], 0.000001); 494 try std.testing.expectApproxEqAbs(1.0, result[3][3], 0.000001); 495 }