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