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