1 /// Custom sized 2-dimension Matrices 2 module gfm.math.matrix; 3 4 import std.math, 5 std.typetuple, 6 std.traits, 7 std..string, 8 std.typecons, 9 std.conv; 10 11 import gfm.math.vector, 12 gfm.math.shapes, 13 gfm.math.quaternion; 14 15 /// Generic non-resizeable matrix with R rows and C columns. 16 /// Intended for 3D use (size 3x3 and 4x4). 17 /// Important: <b>Matrices here are in row-major order whereas OpenGL is column-major.</b> 18 /// Params: 19 /// T = type of elements 20 /// R = number of rows 21 /// C = number of columns 22 struct Matrix(T, int R, int C) 23 { 24 public 25 { 26 static assert(R >= 1 && C >= 1); 27 28 alias Vector!(T, C) row_t; 29 alias Vector!(T, R) column_t; 30 31 enum bool isSquare = (R == C); 32 33 // fields definition 34 union 35 { 36 T[C*R] v; // all elements 37 row_t[R] rows; // all rows 38 T[C][R] c; // components 39 } 40 41 @nogc this(U...)(U values) pure nothrow 42 { 43 static if ((U.length == C*R) && allSatisfy!(isTAssignable, U)) 44 { 45 // construct with components 46 foreach(int i, x; values) 47 v[i] = x; 48 } 49 else static if ((U.length == 1) && (isAssignable!(U[0])) && (!is(U[0] : Matrix))) 50 { 51 // construct with assignment 52 opAssign!(U[0])(values[0]); 53 } 54 else static assert(false, "cannot create a matrix from given arguments"); 55 } 56 57 /// Construct a matrix from columns. 58 @nogc static Matrix fromColumns(column_t[] columns) pure nothrow 59 { 60 assert(columns.length == C); 61 Matrix res; 62 for (int i = 0; i < R; ++i) 63 for (int j = 0; j < C; ++j) 64 { 65 res.c[i][j] = columns[j][i]; 66 } 67 return res; 68 } 69 70 /// Construct a matrix from rows. 71 @nogc static Matrix fromRows(row_t[] rows) pure nothrow 72 { 73 assert(rows.length == R); 74 Matrix res; 75 res.rows[] = rows[]; 76 return res; 77 } 78 79 /// Construct matrix with a scalar. 80 @nogc this(U)(T x) pure nothrow 81 { 82 for (int i = 0; i < _N; ++i) 83 v[i] = x; 84 } 85 86 /// Assign with a scalar. 87 @nogc ref Matrix opAssign(U : T)(U x) pure nothrow 88 { 89 for (int i = 0; i < R * C; ++i) 90 v[i] = x; 91 return this; 92 } 93 94 /// Assign with a samey matrice. 95 @nogc ref Matrix opAssign(U : Matrix)(U x) pure nothrow 96 { 97 for (int i = 0; i < R * C; ++i) 98 v[i] = x.v[i]; 99 return this; 100 } 101 102 /// Assign from other small matrices (same size, compatible type). 103 @nogc ref Matrix opAssign(U)(U x) pure nothrow 104 if (isMatrixInstantiation!U 105 && is(U._T : _T) 106 && (!is(U: Matrix)) 107 && (U._R == R) && (U._C == C)) 108 { 109 for (int i = 0; i < R * C; ++i) 110 v[i] = x.v[i]; 111 return this; 112 } 113 114 /// Assign with a static array of size R * C. 115 @nogc ref Matrix opAssign(U)(U x) pure nothrow 116 if ((isStaticArray!U) 117 && is(typeof(x[0]) : T) 118 && (U.length == R * C)) 119 { 120 for (int i = 0; i < R * C; ++i) 121 v[i] = x[i]; 122 return this; 123 } 124 125 /// Assign with a static array of shape (R, C). 126 @nogc ref Matrix opAssign(U)(U x) pure nothrow 127 if ((isStaticArray!U) && isStaticArray!(typeof(x[0])) 128 && is(typeof(x[0][0]) : T) 129 && (U.length == R) 130 && (x[0].length == C)) 131 { 132 foreach (i; 0..R) 133 rows[i] = x[i]; 134 return this; 135 } 136 137 /// Assign with a dynamic array of size R * C. 138 @nogc ref Matrix opAssign(U)(U x) pure nothrow 139 if ((isDynamicArray!U) 140 && is(typeof(x[0]) : T)) 141 { 142 assert(x.length == R * C); 143 for (int i = 0; i < R * C; ++i) 144 v[i] = x[i]; 145 return this; 146 } 147 148 /// Assign with a dynamic array of shape (R, C). 149 @nogc ref Matrix opAssign(U)(U x) pure nothrow 150 if ((isDynamicArray!U) && isDynamicArray!(typeof(x[0])) 151 && is(typeof(x[0][0]) : T)) 152 { 153 assert(x.length == R); 154 foreach (i; 0..R) 155 { 156 assert(x[i].length == C); 157 rows[i] = x[i]; 158 } 159 return this; 160 } 161 162 /// Return a pointer to content. 163 @nogc inout(T)* ptr() pure inout nothrow @property 164 { 165 return v.ptr; 166 } 167 168 /// Returns a column as a vector 169 /// Returns: column j as a vector. 170 @nogc column_t column(int j) pure const nothrow 171 { 172 column_t res = void; 173 for (int i = 0; i < R; ++i) 174 res.v[i] = c[i][j]; 175 return res; 176 } 177 178 /// Returns a row as a vector 179 /// Returns: row i as a vector. 180 @nogc row_t row(int i) pure const nothrow 181 { 182 return rows[i]; 183 } 184 185 /// Covnerts to pretty string. 186 string toString() const nothrow 187 { 188 try 189 return format("%s", v); 190 catch (Exception e) 191 assert(false); // should not happen since format is right 192 } 193 194 /// Matrix * scalar multiplication. 195 @nogc Matrix opBinary(string op)(T factor) pure const nothrow if (op == "*") 196 { 197 Matrix result = void; 198 199 for (int i = 0; i < R; ++i) 200 { 201 for (int j = 0; j < C; ++j) 202 { 203 result.c[i][j] = c[i][j] * factor; 204 } 205 } 206 return result; 207 } 208 209 /// Matrix * vector multiplication. 210 @nogc column_t opBinary(string op)(row_t x) pure const nothrow if (op == "*") 211 { 212 column_t res = void; 213 for (int i = 0; i < R; ++i) 214 { 215 T sum = 0; 216 for (int j = 0; j < C; ++j) 217 { 218 sum += c[i][j] * x.v[j]; 219 } 220 res.v[i] = sum; 221 } 222 return res; 223 } 224 225 /// Matrix * matrix multiplication. 226 @nogc auto opBinary(string op, U)(U x) pure const nothrow 227 if (isMatrixInstantiation!U && (U._R == C) && (op == "*")) 228 { 229 Matrix!(T, R, U._C) result = void; 230 231 for (int i = 0; i < R; ++i) 232 { 233 for (int j = 0; j < U._C; ++j) 234 { 235 T sum = 0; 236 for (int k = 0; k < C; ++k) 237 sum += c[i][k] * x.c[k][j]; 238 result.c[i][j] = sum; 239 } 240 } 241 return result; 242 } 243 244 /// Matrix add and substraction. 245 @nogc Matrix opBinary(string op, U)(U other) pure const nothrow 246 if (is(U : Matrix) && (op == "+" || op == "-")) 247 { 248 Matrix result = void; 249 250 for (int i = 0; i < R; ++i) 251 { 252 for (int j = 0; j < C; ++j) 253 { 254 mixin("result.c[i][j] = c[i][j] " ~ op ~ " other.c[i][j];"); 255 } 256 } 257 return result; 258 } 259 260 // matrix *= scalar 261 @nogc ref Matrix opOpAssign(string op, U : T)(U x) pure nothrow if (op == "*") 262 { 263 for (int i = 0; i < R * C; ++i) 264 v[i] *= x; 265 return this; 266 } 267 268 /// Assignment operator with another samey matrix. 269 @nogc ref Matrix opOpAssign(string op, U)(U operand) pure nothrow 270 if (is(U : Matrix) && (op == "*" || op == "+" || op == "-")) 271 { 272 mixin("Matrix result = this " ~ op ~ " operand;"); 273 return opAssign!Matrix(result); 274 } 275 276 /// Matrix += <something convertible to a Matrix> 277 /// Matrix -= <something convertible to a Matrix> 278 @nogc ref Matrix opOpAssign(string op, U)(U operand) pure nothrow 279 if ((isConvertible!U) && (op == "*" || op == "+" || op == "-")) 280 { 281 Matrix conv = operand; 282 return opOpAssign!op(conv); 283 } 284 285 /// Cast to other matrix types. 286 /// If the size are different, the resulting matrix is truncated 287 /// and/or filled with identity coefficients. 288 @nogc U opCast(U)() pure const nothrow if (isMatrixInstantiation!U) 289 { 290 U res = U.identity(); 291 enum minR = R < U._R ? R : U._R; 292 enum minC = C < U._C ? C : U._C; 293 for (int i = 0; i < minR; ++i) 294 for (int j = 0; j < minC; ++j) 295 { 296 res.c[i][j] = cast(U._T)(c[i][j]); 297 } 298 return res; 299 } 300 301 @nogc bool opEquals(U)(U other) pure const nothrow if (is(U : Matrix)) 302 { 303 for (int i = 0; i < R * C; ++i) 304 if (v[i] != other.v[i]) 305 return false; 306 return true; 307 } 308 309 @nogc bool opEquals(U)(U other) pure const nothrow 310 if ((isAssignable!U) && (!is(U: Matrix))) 311 { 312 Matrix conv = other; 313 return opEquals(conv); 314 } 315 316 // +matrix, -matrix, ~matrix, !matrix 317 @nogc Matrix opUnary(string op)() pure const nothrow if (op == "+" || op == "-" || op == "~" || op == "!") 318 { 319 Matrix res = void; 320 for (int i = 0; i < N; ++i) 321 mixin("res.v[i] = " ~ op ~ "v[i];"); 322 return res; 323 } 324 325 /// Convert 3x3 rotation matrix to quaternion. 326 /// See_also: 3D Math Primer for Graphics and Game Development. 327 @nogc U opCast(U)() pure const nothrow if (isQuaternionInstantiation!U 328 && is(U._T : _T) 329 && (_R == 3) && (_C == 3)) 330 { 331 T fourXSquaredMinus1 = c[0][0] - c[1][1] - c[2][2]; 332 T fourYSquaredMinus1 = c[1][1] - c[0][0] - c[2][2]; 333 T fourZSquaredMinus1 = c[2][2] - c[0][0] - c[1][1]; 334 T fourWSquaredMinus1 = c[0][0] + c[1][1] + c[2][2]; 335 336 int biggestIndex = 0; 337 T fourBiggestSquaredMinus1 = fourWSquaredMinus1; 338 339 if(fourXSquaredMinus1 > fourBiggestSquaredMinus1) 340 { 341 fourBiggestSquaredMinus1 = fourXSquaredMinus1; 342 biggestIndex = 1; 343 } 344 345 if(fourYSquaredMinus1 > fourBiggestSquaredMinus1) 346 { 347 fourBiggestSquaredMinus1 = fourYSquaredMinus1; 348 biggestIndex = 2; 349 } 350 351 if(fourZSquaredMinus1 > fourBiggestSquaredMinus1) 352 { 353 fourBiggestSquaredMinus1 = fourZSquaredMinus1; 354 biggestIndex = 3; 355 } 356 357 T biggestVal = sqrt(fourBiggestSquaredMinus1 + 1) / 2; 358 T mult = 1 / (biggestVal * 4); 359 360 U quat; 361 switch(biggestIndex) 362 { 363 case 1: 364 quat.w = (c[1][2] - c[2][1]) * mult; 365 quat.x = biggestVal; 366 quat.y = (c[0][1] + c[1][0]) * mult; 367 quat.z = (c[2][0] + c[0][2]) * mult; 368 break; 369 370 case 2: 371 quat.w = (c[2][0] - c[0][2]) * mult; 372 quat.x = (c[0][1] + c[1][0]) * mult; 373 quat.y = biggestVal; 374 quat.z = (c[1][2] + c[2][1]) * mult; 375 break; 376 377 case 3: 378 quat.w = (c[0][1] - c[1][0]) * mult; 379 quat.x = (c[2][0] + c[0][2]) * mult; 380 quat.y = (c[1][2] + c[2][1]) * mult; 381 quat.z = biggestVal; 382 break; 383 384 default: // biggestIndex == 0 385 quat.w = biggestVal; 386 quat.x = (c[1][2] - c[2][1]) * mult; 387 quat.y = (c[2][0] - c[0][2]) * mult; 388 quat.z = (c[0][1] - c[1][0]) * mult; 389 break; 390 } 391 392 return quat; 393 } 394 395 /// Converts a 4x4 rotation matrix to quaternion. 396 @nogc U opCast(U)() pure const nothrow if (isQuaternionInstantiation!U 397 && is(U._T : _T) 398 && (_R == 4) && (_C == 4)) 399 { 400 auto m3 = cast(mat3!T)(this); 401 return cast(U)(m3); 402 } 403 404 static if (isSquare && isFloatingPoint!T && R == 1) 405 { 406 /// Returns an inverted copy of this matrix 407 /// Returns: inverse of matrix. 408 /// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices. 409 @nogc Matrix inverse() pure const nothrow 410 { 411 assert(c[0][0] != 0); // Programming error if matrix is not invertible. 412 return Matrix( 1 / c[0][0]); 413 } 414 } 415 416 static if (isSquare && isFloatingPoint!T && R == 2) 417 { 418 /// Returns an inverted copy of this matrix 419 /// Returns: inverse of matrix. 420 /// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices. 421 @nogc Matrix inverse() pure const nothrow 422 { 423 T det = (c[0][0] * c[1][1] - c[0][1] * c[1][0]); 424 assert(det != 0); // Programming error if matrix is not invertible. 425 T invDet = 1 / det; 426 return Matrix( c[1][1] * invDet, -c[0][1] * invDet, 427 -c[1][0] * invDet, c[0][0] * invDet); 428 } 429 } 430 431 static if (isSquare && isFloatingPoint!T && R == 3) 432 { 433 /// Returns an inverted copy of this matrix 434 /// Returns: inverse of matrix. 435 /// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices. 436 @nogc Matrix inverse() pure const nothrow 437 { 438 T det = c[0][0] * (c[1][1] * c[2][2] - c[2][1] * c[1][2]) 439 - c[0][1] * (c[1][0] * c[2][2] - c[1][2] * c[2][0]) 440 + c[0][2] * (c[1][0] * c[2][1] - c[1][1] * c[2][0]); 441 assert(det != 0); // Programming error if matrix is not invertible. 442 T invDet = 1 / det; 443 444 Matrix res = void; 445 res.c[0][0] = (c[1][1] * c[2][2] - c[2][1] * c[1][2]) * invDet; 446 res.c[0][1] = -(c[0][1] * c[2][2] - c[0][2] * c[2][1]) * invDet; 447 res.c[0][2] = (c[0][1] * c[1][2] - c[0][2] * c[1][1]) * invDet; 448 res.c[1][0] = -(c[1][0] * c[2][2] - c[1][2] * c[2][0]) * invDet; 449 res.c[1][1] = (c[0][0] * c[2][2] - c[0][2] * c[2][0]) * invDet; 450 res.c[1][2] = -(c[0][0] * c[1][2] - c[1][0] * c[0][2]) * invDet; 451 res.c[2][0] = (c[1][0] * c[2][1] - c[2][0] * c[1][1]) * invDet; 452 res.c[2][1] = -(c[0][0] * c[2][1] - c[2][0] * c[0][1]) * invDet; 453 res.c[2][2] = (c[0][0] * c[1][1] - c[1][0] * c[0][1]) * invDet; 454 return res; 455 } 456 } 457 458 static if (isSquare && isFloatingPoint!T && R == 4) 459 { 460 /// Returns an inverted copy of this matrix 461 /// Returns: inverse of matrix. 462 /// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices. 463 @nogc Matrix inverse() pure const nothrow 464 { 465 T det2_01_01 = c[0][0] * c[1][1] - c[0][1] * c[1][0]; 466 T det2_01_02 = c[0][0] * c[1][2] - c[0][2] * c[1][0]; 467 T det2_01_03 = c[0][0] * c[1][3] - c[0][3] * c[1][0]; 468 T det2_01_12 = c[0][1] * c[1][2] - c[0][2] * c[1][1]; 469 T det2_01_13 = c[0][1] * c[1][3] - c[0][3] * c[1][1]; 470 T det2_01_23 = c[0][2] * c[1][3] - c[0][3] * c[1][2]; 471 472 T det3_201_012 = c[2][0] * det2_01_12 - c[2][1] * det2_01_02 + c[2][2] * det2_01_01; 473 T det3_201_013 = c[2][0] * det2_01_13 - c[2][1] * det2_01_03 + c[2][3] * det2_01_01; 474 T det3_201_023 = c[2][0] * det2_01_23 - c[2][2] * det2_01_03 + c[2][3] * det2_01_02; 475 T det3_201_123 = c[2][1] * det2_01_23 - c[2][2] * det2_01_13 + c[2][3] * det2_01_12; 476 477 T det = - det3_201_123 * c[3][0] + det3_201_023 * c[3][1] - det3_201_013 * c[3][2] + det3_201_012 * c[3][3]; 478 assert(det != 0); // Programming error if matrix is not invertible. 479 T invDet = 1 / det; 480 481 T det2_03_01 = c[0][0] * c[3][1] - c[0][1] * c[3][0]; 482 T det2_03_02 = c[0][0] * c[3][2] - c[0][2] * c[3][0]; 483 T det2_03_03 = c[0][0] * c[3][3] - c[0][3] * c[3][0]; 484 T det2_03_12 = c[0][1] * c[3][2] - c[0][2] * c[3][1]; 485 T det2_03_13 = c[0][1] * c[3][3] - c[0][3] * c[3][1]; 486 T det2_03_23 = c[0][2] * c[3][3] - c[0][3] * c[3][2]; 487 T det2_13_01 = c[1][0] * c[3][1] - c[1][1] * c[3][0]; 488 T det2_13_02 = c[1][0] * c[3][2] - c[1][2] * c[3][0]; 489 T det2_13_03 = c[1][0] * c[3][3] - c[1][3] * c[3][0]; 490 T det2_13_12 = c[1][1] * c[3][2] - c[1][2] * c[3][1]; 491 T det2_13_13 = c[1][1] * c[3][3] - c[1][3] * c[3][1]; 492 T det2_13_23 = c[1][2] * c[3][3] - c[1][3] * c[3][2]; 493 494 T det3_203_012 = c[2][0] * det2_03_12 - c[2][1] * det2_03_02 + c[2][2] * det2_03_01; 495 T det3_203_013 = c[2][0] * det2_03_13 - c[2][1] * det2_03_03 + c[2][3] * det2_03_01; 496 T det3_203_023 = c[2][0] * det2_03_23 - c[2][2] * det2_03_03 + c[2][3] * det2_03_02; 497 T det3_203_123 = c[2][1] * det2_03_23 - c[2][2] * det2_03_13 + c[2][3] * det2_03_12; 498 499 T det3_213_012 = c[2][0] * det2_13_12 - c[2][1] * det2_13_02 + c[2][2] * det2_13_01; 500 T det3_213_013 = c[2][0] * det2_13_13 - c[2][1] * det2_13_03 + c[2][3] * det2_13_01; 501 T det3_213_023 = c[2][0] * det2_13_23 - c[2][2] * det2_13_03 + c[2][3] * det2_13_02; 502 T det3_213_123 = c[2][1] * det2_13_23 - c[2][2] * det2_13_13 + c[2][3] * det2_13_12; 503 504 T det3_301_012 = c[3][0] * det2_01_12 - c[3][1] * det2_01_02 + c[3][2] * det2_01_01; 505 T det3_301_013 = c[3][0] * det2_01_13 - c[3][1] * det2_01_03 + c[3][3] * det2_01_01; 506 T det3_301_023 = c[3][0] * det2_01_23 - c[3][2] * det2_01_03 + c[3][3] * det2_01_02; 507 T det3_301_123 = c[3][1] * det2_01_23 - c[3][2] * det2_01_13 + c[3][3] * det2_01_12; 508 509 Matrix res = void; 510 res.c[0][0] = - det3_213_123 * invDet; 511 res.c[1][0] = + det3_213_023 * invDet; 512 res.c[2][0] = - det3_213_013 * invDet; 513 res.c[3][0] = + det3_213_012 * invDet; 514 515 res.c[0][1] = + det3_203_123 * invDet; 516 res.c[1][1] = - det3_203_023 * invDet; 517 res.c[2][1] = + det3_203_013 * invDet; 518 res.c[3][1] = - det3_203_012 * invDet; 519 520 res.c[0][2] = + det3_301_123 * invDet; 521 res.c[1][2] = - det3_301_023 * invDet; 522 res.c[2][2] = + det3_301_013 * invDet; 523 res.c[3][2] = - det3_301_012 * invDet; 524 525 res.c[0][3] = - det3_201_123 * invDet; 526 res.c[1][3] = + det3_201_023 * invDet; 527 res.c[2][3] = - det3_201_013 * invDet; 528 res.c[3][3] = + det3_201_012 * invDet; 529 return res; 530 } 531 } 532 533 /// Returns a transposed copy of this matrix 534 /// Returns: transposed matrice. 535 @nogc Matrix!(T, C, R) transposed() pure const nothrow 536 { 537 Matrix!(T, C, R) res; 538 for (int i = 0; i < C; ++i) 539 for (int j = 0; j < R; ++j) 540 res.c[i][j] = c[j][i]; 541 return res; 542 } 543 544 static if (isSquare && R > 1) 545 { 546 /// Makes a diagonal matrix from a vector. 547 @nogc static Matrix diag(Vector!(T, R) v) pure nothrow 548 { 549 Matrix res = void; 550 for (int i = 0; i < R; ++i) 551 for (int j = 0; j < C; ++j) 552 res.c[i][j] = (i == j) ? v.v[i] : 0; 553 return res; 554 } 555 556 /// In-place translate by (v, 1) 557 @nogc void translate(Vector!(T, R-1) v) pure nothrow 558 { 559 for (int i = 0; i < R; ++i) 560 { 561 T dot = 0; 562 for (int j = 0; j + 1 < C; ++j) 563 dot += v.v[j] * c[i][j]; 564 565 c[i][C-1] += dot; 566 } 567 } 568 569 /// Make a translation matrix. 570 @nogc static Matrix translation(Vector!(T, R-1) v) pure nothrow 571 { 572 Matrix res = identity(); 573 for (int i = 0; i + 1 < R; ++i) 574 res.c[i][C-1] += v.v[i]; 575 return res; 576 } 577 578 /// In-place matrix scaling. 579 void scale(Vector!(T, R-1) v) pure nothrow 580 { 581 for (int i = 0; i < R; ++i) 582 for (int j = 0; j + 1 < C; ++j) 583 c[i][j] *= v.v[j]; 584 } 585 586 /// Make a scaling matrix. 587 @nogc static Matrix scaling(Vector!(T, R-1) v) pure nothrow 588 { 589 Matrix res = identity(); 590 for (int i = 0; i + 1 < R; ++i) 591 res.c[i][i] = v.v[i]; 592 return res; 593 } 594 } 595 596 // rotations are implemented for 3x3 and 4x4 matrices. 597 static if (isSquare && (R == 3 || R == 4) && isFloatingPoint!T) 598 { 599 @nogc public static Matrix rotateAxis(int i, int j)(T angle) pure nothrow 600 { 601 Matrix res = identity(); 602 const T cosa = cos(angle); 603 const T sina = sin(angle); 604 res.c[i][i] = cosa; 605 res.c[i][j] = -sina; 606 res.c[j][i] = sina; 607 res.c[j][j] = cosa; 608 return res; 609 } 610 611 /// Rotate along X axis 612 /// Returns: rotation matrix along axis X 613 alias rotateAxis!(1, 2) rotateX; 614 615 /// Rotate along Y axis 616 /// Returns: rotation matrix along axis Y 617 alias rotateAxis!(2, 0) rotateY; 618 619 /// Rotate along Z axis 620 /// Returns: rotation matrix along axis Z 621 alias rotateAxis!(0, 1) rotateZ; 622 623 /// Similar to the glRotate matrix, however the angle is expressed in radians 624 /// See_also: $(LINK http://www.cs.rutgers.edu/~decarlo/428/gl_man/rotate.html) 625 @nogc static Matrix rotation(T angle, vec3!T axis) pure nothrow 626 { 627 Matrix res = identity(); 628 const T c = cos(angle); 629 const oneMinusC = 1 - c; 630 const T s = sin(angle); 631 axis = axis.normalized(); 632 T x = axis.x, 633 y = axis.y, 634 z = axis.z; 635 T xy = x * y, 636 yz = y * z, 637 xz = x * z; 638 639 res.c[0][0] = x * x * oneMinusC + c; 640 res.c[0][1] = x * y * oneMinusC - z * s; 641 res.c[0][2] = x * z * oneMinusC + y * s; 642 res.c[1][0] = y * x * oneMinusC + z * s; 643 res.c[1][1] = y * y * oneMinusC + c; 644 res.c[1][2] = y * z * oneMinusC - x * s; 645 res.c[2][0] = z * x * oneMinusC - y * s; 646 res.c[2][1] = z * y * oneMinusC + x * s; 647 res.c[2][2] = z * z * oneMinusC + c; 648 return res; 649 } 650 } 651 652 // 4x4 specific transformations for 3D usage 653 static if (isSquare && R == 4 && isFloatingPoint!T) 654 { 655 /// Orthographic projection 656 /// Returns: orthographic projection. 657 @nogc static Matrix orthographic(T left, T right, T bottom, T top, T near, T far) pure nothrow 658 { 659 T dx = right - left, 660 dy = top - bottom, 661 dz = far - near; 662 663 T tx = -(right + left) / dx; 664 T ty = -(top + bottom) / dy; 665 T tz = -(far + near) / dz; 666 667 return Matrix(2 / dx, 0, 0, tx, 668 0, 2 / dy, 0, ty, 669 0, 0, -2 / dz, tz, 670 0, 0, 0, 1); 671 } 672 673 /// Perspective projection 674 /// Returns: perspective projection. 675 @nogc static Matrix perspective(T FOVInRadians, T aspect, T zNear, T zFar) pure nothrow 676 { 677 T f = 1 / tan(FOVInRadians / 2); 678 T d = 1 / (zNear - zFar); 679 680 return Matrix(f / aspect, 0, 0, 0, 681 0, f, 0, 0, 682 0, 0, (zFar + zNear) * d, 2 * d * zFar * zNear, 683 0, 0, -1, 0); 684 } 685 686 /// Look At projection 687 /// Returns: "lookAt" projection. 688 /// Thanks to vuaru for corrections. 689 @nogc static Matrix lookAt(vec3!T eye, vec3!T target, vec3!T up) pure nothrow 690 { 691 vec3!T Z = (eye - target).normalized(); 692 vec3!T X = cross(-up, Z).normalized(); 693 vec3!T Y = cross(Z, -X); 694 695 return Matrix(-X.x, -X.y, -X.z, dot(X, eye), 696 Y.x, Y.y, Y.z, -dot(Y, eye), 697 Z.x, Z.y, Z.z, -dot(Z, eye), 698 0, 0, 0, 1); 699 } 700 701 /// Extract frustum from a 4x4 matrice. 702 @nogc Frustum!T frustum() pure const nothrow 703 { 704 auto left = Plane!T(row(3) + row(0)); 705 auto right = Plane!T(row(3) - row(0)); 706 auto top = Plane!T(row(3) - row(1)); 707 auto bottom = Plane!T(row(3) + row(1)); 708 auto near = Plane!T(row(3) + row(2)); 709 auto far = Plane!T(row(3) - row(2)); 710 return Frustum!T(left, right, top, bottom, near, far); 711 } 712 713 } 714 } 715 716 package 717 { 718 alias T _T; 719 enum _R = R; 720 enum _C = C; 721 } 722 723 private 724 { 725 template isAssignable(T) 726 { 727 enum bool isAssignable = std.traits.isAssignable!(Matrix, T); 728 } 729 730 template isConvertible(T) 731 { 732 enum bool isConvertible = (!is(T : Matrix)) && isAssignable!T; 733 } 734 735 template isTAssignable(U) 736 { 737 enum bool isTAssignable = std.traits.isAssignable!(T, U); 738 } 739 740 template isRowConvertible(U) 741 { 742 enum bool isRowConvertible = is(U : row_t); 743 } 744 745 template isColumnConvertible(U) 746 { 747 enum bool isColumnConvertible = is(U : column_t); 748 } 749 } 750 751 public 752 { 753 /// Construct an identity matrix 754 /// Returns: an identity matrix. 755 /// Note: the identity matrix, while only meaningful for square matrices, 756 /// is also defined for non-square ones. 757 @nogc static Matrix identity() pure nothrow 758 { 759 Matrix res = void; 760 for (int i = 0; i < R; ++i) 761 for (int j = 0; j < C; ++j) 762 res.c[i][j] = (i == j) ? 1 : 0; 763 return res; 764 } 765 766 /// Construct an constant matrix 767 /// Returns: a constant matrice. 768 @nogc static Matrix constant(U)(U x) pure nothrow 769 { 770 Matrix res = void; 771 772 for (int i = 0; i < R * C; ++i) 773 res.v[i] = cast(T)x; 774 return res; 775 } 776 } 777 } 778 779 template isMatrixInstantiation(U) 780 { 781 private static void isMatrix(T, int R, int C)(Matrix!(T, R, C) x) 782 { 783 } 784 785 enum bool isMatrixInstantiation = is(typeof(isMatrix(U.init))); 786 } 787 788 // GLSL is a big inspiration here 789 // we defines types with more or less the same names 790 791 /// 792 template mat2x2(T) { alias Matrix!(T, 2, 2) mat2x2; } 793 /// 794 template mat3x3(T) { alias Matrix!(T, 3, 3) mat3x3; } 795 /// 796 template mat4x4(T) { alias Matrix!(T, 4, 4) mat4x4; } 797 798 // WARNING: in GLSL, first number is _columns_, second is rows 799 // It is the opposite here: first number is rows, second is columns 800 // With this convention mat2x3 * mat3x4 -> mat2x4. 801 802 /// 803 template mat2x3(T) { alias Matrix!(T, 2, 3) mat2x3; } 804 /// 805 template mat2x4(T) { alias Matrix!(T, 2, 4) mat2x4; } 806 /// 807 template mat3x2(T) { alias Matrix!(T, 3, 2) mat3x2; } 808 /// 809 template mat3x4(T) { alias Matrix!(T, 3, 4) mat3x4; } 810 /// 811 template mat4x2(T) { alias Matrix!(T, 4, 2) mat4x2; } 812 /// 813 template mat4x3(T) { alias Matrix!(T, 4, 3) mat4x3; } 814 815 // shorter names for most common matrices 816 alias mat2x2 mat2;/// 817 alias mat3x3 mat3;/// 818 alias mat4x4 mat4;/// 819 820 // Define a lot of type names 821 // Most useful are probably mat4f and mat4d 822 823 alias mat2!float mat2f;/// 824 alias mat2!double mat2d;/// 825 826 alias mat3!float mat3f;/// 827 alias mat3!double mat3d;/// 828 829 alias mat4!float mat4f;/// 830 alias mat4!double mat4d;/// 831 832 alias mat2x2!float mat2x2f;/// 833 alias mat2x2!double mat2x2d;/// 834 835 alias mat3x3!float mat3x3f;/// 836 alias mat3x3!double mat3x3d;/// 837 838 alias mat4x4!float mat4x4f;/// 839 alias mat4x4!double mat4x4d;/// 840 841 unittest 842 { 843 alias mat2i = mat2!int; 844 alias mat2x3f = mat2x3!float; 845 alias mat3x4f = mat3x4!float; 846 alias mat2x4f = mat2x4!float; 847 848 mat2i x = mat2i(0, 1, 849 2, 3); 850 assert(x.c[0][0] == 0 && x.c[0][1] == 1 && x.c[1][0] == 2 && x.c[1][1] == 3); 851 852 vec2i[2] cols = [vec2i(0, 2), vec2i(1, 3)]; 853 mat2i y = mat2i.fromColumns(cols[]); 854 assert(y.c[0][0] == 0 && y.c[0][1] == 1 && y.c[1][0] == 2 && y.c[1][1] == 3); 855 y = mat2i.fromRows(cols[]); 856 assert(y.c[0][0] == 0 && y.c[1][0] == 1 && y.c[0][1] == 2 && y.c[1][1] == 3); 857 y = y.transposed(); 858 859 assert(x == y); 860 x = [0, 1, 2, 3]; 861 assert(x == y); 862 863 mat2i z = x * y; 864 assert(z == mat2i([2, 3, 6, 11])); 865 vec2i vz = z * vec2i(2, -1); 866 assert(vz == vec2i(1, 1)); 867 868 mat2f a = z; 869 mat2d ad = a; 870 ad += a; 871 mat2f w = [4, 5, 6, 7]; 872 z = cast(mat2i)w; 873 assert(w == z); 874 875 { 876 mat2x3f A; 877 mat3x4f B; 878 mat2x4f C = A * B; 879 } 880 881 assert(mat2i.diag(vec2i(1, 2)) == mat2i(1, 0, 882 0, 2)); 883 884 // Construct with a single scalar 885 auto D = mat4f(1.0f); 886 assert(D.v[] == [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ]); 887 888 { 889 double[4][3] starray = [ 890 [ 0, 1, 2, 3], 891 [ 4, 5, 6, 7,], 892 [ 8, 9, 10, 11,], 893 ]; 894 895 // starray has the shape 3x4 896 assert(starray.length == 3); 897 assert(starray[0].length == 4); 898 899 auto m = mat3x4!double(starray); 900 assert(m.v[] == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ]); 901 } 902 903 { 904 auto dynarray = [ 905 [ 0, 1, 2, 3], 906 [ 4, 5, 6, 7,], 907 [ 8, 9, 10, 11,], 908 ]; 909 910 // dynarray has the shape 3x4 911 assert(dynarray.length == 3); 912 assert(dynarray[0].length == 4); 913 914 auto m = mat3x4!double(dynarray); 915 assert(m.v[] == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ]); 916 } 917 } 918 919 // Issue #206 (matrix *= scalar) not yielding matrix * scalar but matrix * matrix(scalar) 920 unittest 921 { 922 mat4f mvp = mat4f.identity; 923 mvp *= 2; 924 assert(mvp == mat4f(2, 0, 0, 0, 925 0, 2, 0, 0, 926 0, 0, 2, 0, 927 0, 0, 0, 2)); 928 929 mvp = mat4f.identity * 2; 930 assert(mvp == mat4f(2, 0, 0, 0, 931 0, 2, 0, 0, 932 0, 0, 2, 0, 933 0, 0, 0, 2)); 934 935 936 mvp = mat4f(1) * mat4f(1); 937 assert(mvp == mat4f(4, 4, 4, 4, 938 4, 4, 4, 4, 939 4, 4, 4, 4, 940 4, 4, 4, 4)); 941 942 mvp = mat4f(1); 943 mvp *= mat4f(1); 944 assert(mvp == mat4f(4, 4, 4, 4, 945 4, 4, 4, 4, 946 4, 4, 4, 4, 947 4, 4, 4, 4)); 948 }