sily.matrix source code

1 /// Flexible templated matrix with some math utils
2 module sily.matrix;
3 
4 import std.math;
5 import std.numeric;
6 import std.conv;
7 import std.traits;
8 import std.typetuple;
9 import std.algorithm;
10 import std.stdio;
11 import std.string;
12 import std.format;
13 
14 import sily.math;
15 import sily.array;
16 import sily.vector;
17 import sily.quat;
18 
19 /// GLSL style aliases
20 alias mat(T, size_t H, size_t W) = Matrix!(T, H, W);
21 /// Ditto
22 alias mat(size_t H, size_t W) = Matrix!(float, H, W);
23 /// Ditto
24 alias mat4 = mat!(4, 4);
25 /// Ditto
26 alias mat3 = mat!(3, 3);
27 /// Ditto
28 alias mat2 = mat!(2, 2);
29 
30 /// Ditto
31 alias imat(size_t H, size_t W) = Matrix!(int, H, W);
32 /// Ditto
33 alias imat4 = imat!(4, 4);
34 /// Ditto
35 alias imat3 = imat!(3, 3);
36 /// Ditto
37 alias imat2 = imat!(2, 2);
38 
39 // mat2x3 // mat[row][col] // mat[1][2] = b
40 // data def = data[row][col] aka T[col][row]
41 // [ [a], [b], [c] ] 
42 // [ [d], [e], [f] ] 
43 /++
44 Matrix implementation, row-major right handed
45 Basic operations:
46 ---
47 // Scalar operations, any of *, /, +, -, % with number
48 m1 + 1;
49 m1 * 2;
50 // etc...
51 // Matrix/Matrix scalar op (on each element)
52 m1 + m2;
53 m1 - m2;
54 m1.scalarMultiply(m2);
55 m1.scalarDivide(m2);
56 // Unary operations
57 -m1; // invert all numbers
58 ~m1; // invert matrix, aka m1^-1, where ~m1*m1=m1.identity
59 // Binary operations
60 m1 * m2; // matrix multiplication
61 vec1x2 * m2x1; // vector by matrix
62 m3x4 * vec4x1; // matrix by vector
63 // Casting
64 cast(vec1x3) mat1x3; // cast to column vector
65 cast(vec3x1) mat1x3; // or row vector
66 cast(mat1x3) vec3x1; // or event vector to mat
67 cast(imat2) fmat2; // and type changing casts
68 cast(imat2) imat4; // changing size (keeps top left)
69 cast(imat4) imat2; // fills new cells with zeros
70 // Data access:
71 m1[1][2]; // access to row 1 column 2
72 m1.data[1][2]; // or explicit access to data
73 ---
74 +/
75 struct Matrix(T, size_t H, size_t W) if (isNumeric!T && W > 0 && H > 0) {
76     /// Matrix data
77     T[W][H] data = 0;
78 
79     /// Alias to allow data access
80     alias data this;
81     /// Alias to data type
82     alias dataType = T;
83 
84     /// Is matrix square (W == H)
85     enum bool isSquare = W == H;
86     
87     /++
88     Alias to matrix type. Can be used to construct matrices
89     of same type
90     ---
91     auto rmat7x2 = Matrix!(real, 7, 2);
92     auto rmat7x2_2 = rmat7x2.MatType;
93     +/
94     alias MatType = Matrix!(T, H, W);
95     /// Matrix size (h*w)
96     enum size_t size = H * W;
97     /// Matrix size (h)
98     enum size_t rows = H;
99     /// Matrix size (w)
100     enum size_t columns = W;
101 
102     /++
103     Constructs Matrix from components. If no components present
104     matrix will be filled with 0
105     Example:
106     ---
107     // Matrix can be constructed manually or with aliases
108     // Fills int matrix 2x2 with 1
109     auto m0 = Matrix!(int, 2, 2)(1);
110     // Fills int matrix 2x2 with [[1, 2], [3, 4]]
111     auto m1 = Matrix!(int, 2, 2)(1, 2, 3, 4);
112     // Same
113     auto m2 = Matrix!(int, 2, 2)([[1, 2], [3, 4]]);
114     // Same
115     mat2 m3 = [1, 2, 3, 4];
116     // Same
117     mat2 m4 = [[1, 2], [3, 4]];
118     // Also you can construct matrix by 
119     // assigning values directly
120     mat!(2, 3) m5 = [[1, 2, 3], [4, 5, 6]];
121     +/
122     this(in T val) {
123         foreach (i; 0..H) {
124             foreach (j; 0..W) {
125                 data[i][j] = val; 
126             }
127         }
128     }
129     /// Ditto
130     this(in T[size] vals...) {
131         foreach (i; 0..H) {
132             foreach (j; 0..W) {
133                 data[i][j] = vals[i * W + j];
134             }    
135         }
136     }
137     /// Ditto
138     this(in T[W][H] vals...) {
139         data = vals;
140     }
141 
142     static if (isSquare && W == 4 && isFloatingPoint!T) {
143         /// Construct matrix from quaternion
144         this(in quat q) {
145             T x2 = 2.0 * q.x * q.x;
146             T y2 = 2.0 * q.y * q.y;
147             T z2 = 2.0 * q.z * q.z;
148             T xy = 2.0 * q.x * q.y;
149             T xz = 2.0 * q.x * q.z;
150             T xw = 2.0 * q.x * q.w;
151             T yz = 2.0 * q.y * q.z;
152             T yw = 2.0 * q.y * q.w;
153             T zw = 2.0 * q.z * q.w;
154 
155             data = [
156                 [1.0 - y2 - z2, xy - zw, xz + yw, 0],
157                 [xy + zw, 1.0 - x2 - z2, yz - xw, 0],
158                 [xz - yw, yz + xw, 1.0 - x2 - y2, 0],
159                 [0.0, 0.0, 0.0, 1.0]
160             ];
161         }
162     }
163     
164 
165     /* -------------------------------------------------------------------------- */
166     /*                         UNARY OPERATIONS OVERRIDES                         */
167     /* -------------------------------------------------------------------------- */
168 
169     /// Scalar matrix addition and subtraction
170     MatType opBinary(string op, R)(in Matrix!(R, H, W) b) const if ( op == "+" || op == "-" ) {
171         MatType ret = MatType();
172         foreach (i; 0 .. H) {
173             foreach (j; 0..W) {
174                 mixin( "ret[i][j] = data[i][j] " ~ op ~ " b.data[i][j];" );
175             }
176         }
177         return ret;
178     }
179 
180     /// Ditto
181     MatType opBinaryRight(string op, R)(in Matrix!(R, H, W) b) const if ( op == "+" || op == "-" ) {
182         MatType ret = MatType();
183         foreach (i; 0 .. H) {
184             foreach (j; 0..W) {
185                 mixin( "ret[i][j] = b.data[i][j] " ~ op ~ " data[i][j];" );
186             }
187         }
188         return ret;
189     }
190 
191     /// Matrix multiplication
192     Matrix!(T, H, U) opBinary(string op, R, size_t U)(in Matrix!(R, W, U) b) const if ( op == "*" ) {
193         Matrix!(T, H, U) ret = Matrix!(T, H, U)(0);
194         foreach (i; 0 .. H) {
195             foreach (j; 0..U) {
196                 foreach (k; 0.. W) {
197                     ret[i][j] += data[i][k] * b.data[k][j];
198                     if (abs(ret[i][j]) <= float.epsilon * 2.0f) ret[i][j] = 0;
199                 }
200             }
201         }
202         return ret;
203     }
204 
205     /// Ditto
206     Matrix!(T, U, W) opBinaryRight(string op, R, size_t U)(in Matrix!(R, U, H) b) const if ( op == "*" ) {
207         Matrix!(T, U, W) ret = Matrix!(T, U, W)(0);
208         foreach (i; 0..U) {
209             foreach (j; 0..W) {
210                 foreach (k; 0..H) {
211                     ret[i][j] += b.data[i][k] * data[k][j];
212                     if (abs(ret[i][j]) <= float.epsilon * 2.0f) ret[i][j] = 0;
213                 }
214             }
215         }
216         return ret;
217     }
218 
219     /// Vector transformation
220     Vector!(R, H) opBinary(string op, R)(in Vector!(R, W) b) const if ( op == "*" ) {
221         Vector!(R, H) ret = 0;
222         foreach (i; 0 .. H) {
223             foreach (k; 0.. W) {
224                 ret[i] += data[i][k] * b.data[k];
225                 if (abs(ret[i]) <= float.epsilon * 2.0f) ret[i] = 0;
226             }
227         }
228         return ret;
229     }
230 
231     /// Ditto
232     Vector!(R, W) opBinaryRight(string op, R)(in Vector!(R, H) b) const if ( op == "*" ) {
233         Vector!(R, W) ret = 0;
234         foreach (j; 0..W) {
235             foreach (k; 0..H) {
236                 ret[j] += b.data[k] * data[k][j];
237                 if (abs(ret[j]) <= float.epsilon * 2.0f) ret[j] = 0;
238             }
239         }
240         return ret;
241     }
242 
243     /// Vector3 transformation
244     Vector!(R, 3) opBinary(string op, R)(in Vector!(R, 3) b) const if ( op == "*" && W == 4 ) {
245         Vector!(R, 3) ret = 0;
246         foreach (i; 0 .. 3) {
247             foreach (k; 0.. W) {
248                 if (k == 3) {
249                     ret[i] += data[i][k];
250                 } else {
251                     ret[i] += data[i][k] * b.data[k];
252                 }
253                 if (abs(ret[i]) <= float.epsilon * 2.0f) ret[i] = 0;
254             }
255         }
256         return ret;
257     }
258 
259     /// Ditto
260     Vector!(R, 3) opBinaryRight(string op, R)(in Vector!(R, 3) b) const if ( op == "*" && H == 4 ) {
261         Vector!(R, 3) ret = 0;
262         foreach (j; 0..3) {
263             foreach (k; 0..H) {
264                 if (k == 3) {
265                     ret[j] += data[k][j];
266                 } else {
267                     ret[j] += b.data[k] * data[k][j];
268                 }
269                 if (abs(ret[j]) <= float.epsilon * 2.0f) ret[j] = 0;
270             }
271         }
272         return ret;
273     }
274 
275     /// Quaternion transformation
276     quat opBinary(string op)(in quat b) const if ( op == "*" && W == 4) {
277         quat ret = 0;
278         foreach (i; 0 .. H) {
279             foreach (k; 0.. W) {
280                 ret[i] += data[i][k] * b.data[k];
281                 if (abs(ret[i]) <= float.epsilon * 2.0f) ret[i] = 0;
282             }
283         }
284         return ret;
285     }
286 
287     /// Ditto 
288     quat opBinaryRight(string op)(in quat b) const if ( op == "*" && H == 4) {
289         quat ret = 0;
290         foreach (j; 0..W) {
291             foreach (k; 0..H) {
292                 ret[j] += b.data[k] * data[k][j];
293                 if (abs(ret[j]) <= float.epsilon * 2.0f) ret[j] = 0;
294             }
295         }
296         return ret;
297     }
298 
299     /// Scalar number operations
300     MatType opBinary(string op, R)(in R b) const if ( isNumeric!R ) {
301         MatType ret = MatType();
302         foreach (i; 0 .. H) {
303             foreach (j; 0..W) {
304                 mixin( "ret[i][j] = data[i][j] " ~ op ~ " b;" );
305             }
306         }
307         return ret;
308     }
309 
310     /// Ditto
311     MatType opBinaryRight(string op, R)(in R b) const if ( isNumeric!R ) {
312         MatType ret = MatType();
313         foreach (i; 0 .. H) {
314             foreach (j; 0..W) {
315                 mixin( "ret[i][j] = b " ~ op ~ " data[i][j];" );
316             }
317         }
318         return ret;
319     }
320 
321     /// opEquals x == y
322     bool opEquals(R)(in Matrix!(R, H, W) b) const {
323         bool eq = true;
324         foreach (i; 0 .. H) {
325             foreach (j; 0..W) {
326                 eq = eq && data[i][j] == b.data[i][j];
327                 if (!eq) break;
328             }
329         }
330         return eq;
331     }
332 
333     /// Ditto
334     bool opEquals(R, size_t V, size_t U)(in Matrix!(R, V, U) b) const if (V != H || U != W) {
335         return false;
336     }
337   
338     /// Ditto
339     bool opEquals(R)(in R[][] b) const if ( isNumeric!R ){
340         bool eq = true;
341         foreach (i; 0 .. H) {
342             foreach (j; 0..W) {
343                 eq = eq && data[i][j] == b[i][j];
344                 if (!eq) break;
345             }
346         }
347         return eq;
348     }
349 
350     /// Ditto 
351     bool opEquals(R)(in R[] b) const if ( isNumeric!R ){
352         bool eq = true;
353         foreach (i; 0 .. H) {
354             foreach (j; 0..W) {
355                 eq = eq && data[i][j] == b[i * W + j];
356                 if (!eq) break;
357             }
358         }
359         return eq;
360     }
361 
362     // Cannot compare matrices
363     // opCmp x [< > <= >=] y
364 
365     /// opUnary [-, +, --, ++, *, ~] x 
366     MatType opUnary(string op)() if(op == "-" || op == "+"){
367         MatType ret = MatType();
368         foreach (i; 0 .. H) {
369             foreach (j; 0..W) {
370                 mixin( "ret[i][j] = " ~ op ~ " data[i][j];" );
371             }
372         }
373         return ret;
374     }
375 
376     /// Invert matrix
377     Matrix!(float, W, W) opUnary(string op)() if(op == "~" && isSquare) {
378         T det = determinant;
379         if (det == 0) throw new Error("Cannot invert single matrix.");
380         Matrix!(float, W, W) inv = 0;
381         T[W][W] adj = adjoint();
382         foreach (i; 0..W) {
383             foreach (j; 0..H) {
384                 inv[i][j] = adj[i][j] / cast(float) det;
385                 if (abs(inv[i][j]) <= float.epsilon * 2.0f) inv[i][j] = 0;
386             }
387         }
388         return inv;
389     }
390     
391     /// Scalar addition and subtraction in place 
392     MatType opOpAssign(string op, R)( in Matrix!(R, H, W) b ) if ( op == "+" || op == "-" ) { 
393         foreach (i; 0 .. H) {
394             foreach (j; 0..W) {
395                 mixin( "data[i] " ~ op ~ "= b.data[i][j];" );
396             }
397         }
398         return this;
399     }
400     
401     // no OpAssign for mat * mat coz this size is constant
402 
403     /// Ditto
404     MatType opOpAssign(string op, R)( in R b ) if ( isNumeric!R ) { 
405         foreach (i; 0 .. H) {
406             foreach (j; 0..W) {
407                 mixin( "data[i] " ~ op ~ "= b" );
408             }
409         }
410         return this;
411     }
412 
413     /// Matrix type conversion
414     R opCast(R)() const if (isMatrix!(R, H, W)){
415         R ret;
416         foreach (i; 0 .. H) {
417             foreach (j; 0 ..W) {
418                 ret[i][j] = cast(R.dataType) data[i][j];
419             }
420         }
421         return ret;
422     }
423     
424     /// Matrix resizing
425     R opCast(R)() const if (isMatrix!(R) && (R.rows != H || R.columns != W)) {
426         size_t V = R.rows;
427         size_t U = R.columns;
428         R ret = 0;
429         size_t _h = V < H ? V : H;
430         size_t _w = U < W ? U : W;
431         foreach (i; 0.._h) {
432             foreach (j; 0.._w) {
433                 ret[i][j] = data[i][j];
434             }
435         }
436         return ret;
437     }
438     
439     /// Matrix to vector cast (column)
440     R opCast(R)() const if (isVector!(R, H) && W == 1) {
441         R ret;
442         foreach (i; 0..H) {
443             ret[i] = cast(ret.dataType) data[i][0];
444         }
445         return ret;
446     }
447     
448     /// Matrix to vector cast (row)
449     R opCast(R)() const if (isVector!(R, W) && H == 1) {
450         R ret;
451         foreach (i; 0..W) {
452             ret[i] = cast(ret.dataType) data[0][i];
453         }
454         return ret;
455     }
456 
457     /// Boolean cast
458     bool opCast(T)() const if (is(T == bool)) {
459         foreach (i; 0 .. H) {
460             foreach (j; 0..W) {
461                 if (data[i][j] != 0) return true;
462             }
463         }
464         return false;
465     }
466 
467     /// Returns hash 
468     size_t toHash() const @safe nothrow {
469         return typeid(data).getHash(&data);
470     }
471     
472     /// Returns copy of matrix
473     public MatType copyof() {
474         return MatType(data);
475     }
476 
477     /// Returns string representation of matrix: `[1.00, 1.00,... , 1.00]`
478     public string toString() const {
479         import std.conv : to;
480         string s;
481         s ~= "[";
482         foreach (i; 0..H) {
483             foreach (j; 0..W) {
484                 s ~= isFloatingPoint!T ? format("%.2f", data[i][j]) : format("%d", data[i][j]);
485                 if (!(i == H - 1 && j == W - 1)) s ~= ", ";
486             }
487         }
488         s ~= "]";
489         return s;
490     }
491     
492     /// Returns string representation of matrix: `1.00, 1.00,... |\n|1.00, ... , 1.00|`
493     public string toStringPretty() const {
494         import std.conv : to;
495         size_t biggest = 1;
496         foreach (i; 0..H) {
497             foreach (j; 0..W) {
498                 string s = isFloatingPoint!T ? format("%.2f", data[i][j]) : format("%d", data[i][j]);
499                 if (s.length > biggest) biggest = s.length;
500             }
501         }
502 
503         string s;
504         foreach (i; 0..H) {
505             s ~= "|";
506             foreach (j; 0..W) {
507                 string _s = isFloatingPoint!T ? format("%.2f", data[i][j]) : format("%d", data[i][j]);
508                 s ~= _s ~ format("%-*s", (biggest + 1) - _s.length, " ");
509             }
510             if (i != H - 1) s ~= "|\n";
511         }
512         s ~= "|";
513         return s;
514     }
515 
516 
517     /// Returns pointer to data
518     public T[W]* ptr() return {
519         return data.ptr;
520     }
521 
522     static if (isSquare) {
523         /// Returns determinant of matrix
524         public T determinant() const {
525             return getDeterminant!(W)(data.to!(T[][]));
526         }
527 
528         private T getDeterminant(size_t S)(T[][] _mat) const if (S > 0) {
529             static if (S == 1) { 
530                 return _mat[0][0];
531             } else {
532                 T d = 0;
533                 T[][] temp = new T[][](S, S);
534                 T sign = 1;
535                 foreach (f; 0..S) {
536                     getCofactor!S(_mat, temp, 0, f);
537                     d += sign * _mat[0][f] * getDeterminant!(S - 1)(temp);
538                     sign = -sign;
539                 }
540                 return d;
541             }
542         }
543 
544         private void getCofactor(size_t S)(T[][] _mat, ref T[][] _temp, size_t y, size_t x) const {
545             int i = 0;
546             int j = 0;
547             foreach (r; 0..S) {
548                 foreach (c; 0..S) {
549                     if (r != y && c != x) {
550                         _temp[i][j++] = _mat[r][c];
551                         if (j == S - 1) {
552                             j = 0;
553                             ++i;
554                         }
555                     }
556                 }
557             }
558         }
559         
560         /// Returns ajdoint matrix
561         public MatType adjoint() const {
562             static if (W == 1) {
563                 return data[0][0];
564             } else {
565                 MatType ret = 0;
566                 int sign = 1;
567                 T[][] temp = new T[][](W, W);
568                 T[][] tdata = data.to!(T[][]);
569                 foreach (i; 0..W) {
570                     foreach (j; 0..W) {
571                         getCofactor!W(tdata, temp, i, j);
572                         sign = ((i + j) % 2 == 0) ? 1 : -1;
573                         ret[j][i] = sign * getDeterminant!(W - 1)(temp);
574                     }
575                 }
576                 return ret;
577             }
578         }
579         
580         /// Returns identity matrix for size
581         public static MatType identity() {
582             MatType ret = 0;
583             foreach (i; 0..W) {
584                 ret[i][i] = 1;
585             }
586             return ret;
587         }
588     
589         static if (isFloatingPoint!T) {
590             /// Inverts current matrix, alias to `data = ~matrix`
591             public MatType invert() {
592                 MatType m = copyof;
593                 data = (~m).data;
594                 return this;
595             }
596             /// Ditto
597             alias inverse = invert;
598         }
599         
600         /// Returns inverted, alias to `return ~matrix`
601         public mat!(H, W) inverted() {
602             MatType m = copyof;
603             return (~m);
604         }
605 
606         /// Ditto
607         alias inversed = inverted;
608 
609     }
610     
611     /// Multiplies each element by each element of matrices
612     public MatType scalarMultiply(R)(in Matrix!(R, W, H) b) {
613         MatType ret = 0;
614         foreach (i; 0 .. H) {
615             foreach (j; 0..W) {
616                 ret[i][j] = data[i][j] * b.data[i][j];
617             }
618         }
619         return ret;
620     }
621 
622     /// Divides each element by each element of matrices
623     public MatType scalarDivide(R)(in Matrix!(R, W, H) b) {
624         MatType ret = 0;
625         foreach (i; 0 .. H) {
626             foreach (j; 0..W) {
627                 ret[i][j] = data[i][j] / b.data[i][j];
628             }
629         }
630         return ret;
631     }
632     
633     /// Transpose of matrix (swaps rows and columns)
634     public Matrix!(T, W, H) transpose() {
635         Matrix!(T, W, H) ret = 0;
636         foreach (i; 0 .. H) {
637             foreach (j; 0..W) {
638                 ret[j][i] = data[i][j];
639             }
640         }
641         return ret;
642     }
643     
644     /// Resize matrix (keep left top). Alias to cast(Matrix!(T, V, U)) Matrix!(T, H, W) 
645     public Matrix!(T, V, U) resize(size_t V, size_t U)() if (V > 0 && U > 0) {
646         return cast(Matrix!(T, V, U)) this;
647     }
648 
649     static if (isSquare && W == 2 && isFloatingPoint!T) {
650     } else
651     static if (isSquare && W == 3 && isFloatingPoint!T) {
652         /// Constructs 2d rotation matrix
653         static MatType rotation(T angle) {
654             return MatType(cos(angle), -sin(angle), 0, sin(angle), cos(angle), 0, 0, 0, 1).fixNegativeZero;
655         }
656 
657         /// Constructs 2d scale matrix
658         static MatType scale(T[2] s, ...) {
659             return MatType(s[0], 0, 0, 0, s[1], 0, 0, 0, 1);
660         }
661         
662         /// Constructs 2d shear matrix
663         static MatType shear(T[2] s, ...) {
664             return MatType(1, s[0], 0, s[1], 1, 0, 0, 0, 1);
665         }
666 
667         /// Constructs 2d translation matrix
668         static MatType translation(T[2] s, ...) {
669             return MatType(1, 0, s[0], 0, 1, s[1], 0, 0, 1);
670         }
671     } else
672     static if (isSquare && W == 4 && isFloatingPoint!T) {
673         /// Constructs 3d rotation matrix from axis and angle
674         static MatType rotation(T[3] axis, T angle) {
675             if (!vec3(axis).isNormalized) axis = vec3(axis).normalize();
676             T x = axis[0];
677             T y = axis[1];
678             T z = axis[2];
679             
680             float sina = sin(angle);
681             float cosa = cos(angle);
682             float t = 1.0f - cosa;
683 
684             MatType ret = MatType(
685                 x * x * t + cosa,     x * y * t - z * sina, x * z * t + y * sina, 0,
686                 y * x * t + z * sina, y * y * t + cosa,     y * z * t - x * sina, 0,
687                 z * x * t - y * sina, z * y * t + x * sina, z * z * t + cosa, 0,
688                 0, 0, 0, 1
689             ).fixNegativeZero;
690 
691             return ret;
692         }
693         
694         /// Constructs 3d rotation matrix from angle on X axis
695         static MatType rotationX(T angle) {
696             return MatType(
697                 1, 0, 0, 0,
698                 0, cos(angle), -sin(angle), 0,
699                 0, sin(angle), cos(angle), 0,
700                 0, 0, 0, 1
701             ).fixNegativeZero;
702         }
703 
704         /// Constructs 3d rotation matrix from angle on Y axis
705         static MatType rotationY(T angle) {
706             return MatType(
707                 cos(angle), 0, sin(angle), 0,
708                 0, 1, 0, 0,
709                 -sin(angle), 0, cos(angle), 0,
710                 0, 0, 0, 1
711             ).fixNegativeZero;
712         }
713 
714         /// Constructs 3d rotation matrix from angle on Z axis
715         static MatType rotationZ(T angle) {
716             return MatType(
717                 cos(angle), -sin(angle), 0, 0,
718                 sin(angle), cos(angle), 0, 0,
719                 0, 0, 1, 0,
720                 0, 0, 0, 1
721             ).fixNegativeZero;
722         }
723         
724         /// Constructs 3d scale matrix
725         static MatType scale(T[3] v, ...) {
726             return MatType(
727                 v[0], 0, 0, 0,
728                 0, v[1], 0, 0,
729                 0, 0, v[2], 0,
730                 0, 0, 0, 1
731             );
732         }
733         
734         /// Constructs 3d translation matrix
735         static MatType translation(T[3] v, ...) {
736             return MatType(
737                 1, 0, 0, v[0],
738                 0, 1, 0, v[1],
739                 0, 0, 1, v[2],
740                 0, 0, 0, 1
741             );
742         }
743         
744         /// Constructs frustum matrix
745         static MatType frustum(T left, T right, T bottom, T top, T near, T far) {
746             T rl = right - left;
747             T tb = top - bottom;
748             T fn = far - near;
749 
750             return MatType(
751                 near * 2.0 / rl, 0, (right + left) / rl, 0,
752                 0, near * 2.0 / tb, (top + bottom) / tb, 0,
753                 0, 0, -(far + near) / fn, -(far * near * 2.0f) / fn,
754                 0, 0, -1.0f, 0
755             );
756         }
757         
758         /++
759         Construct perspective matrix
760         Params:
761             fovy = vertical fov in degrees (int) or radians (double)
762             aspect = screen.width / screen.height
763             near = near cutoff plane
764             fat = far cutoff plane
765         +/
766         static MatType perspective(int fovy, T aspect, T near, T far) {
767             return perspective( fovy * deg2rad, aspect, near, far );
768         }
769 
770         /++ Ditto +/ 
771         static MatType perspective(T fovy, T aspect, T near, T far) {
772             T top = near * tan(fovy * 0.5);
773             T bottom = -top;
774             T right = top * aspect;
775             T left = -right;
776 
777             T rl = right - left;
778             T tb = top - bottom;
779             T fn = far - near;
780             
781             return MatType(
782                 near * 2.0 / rl, 0, (right + left) / rl, 0,
783                 0, near * 2.0 / tb, (top + bottom) / tb, 0,
784                 0, 0, -(far + near) / fn, -(far * near * 2.0) / fn,
785                 0, 0, -1.0f, 0
786             );
787         }
788         
789         /// Constructs orthographic matrix
790         static MatType ortho(T left, T right, T bottom, T top, T near, T far) {
791             T rl = right - left;
792             T tb = top - bottom;
793             T fn = far - near;
794             
795             return MatType(
796                 2.0 / rl, 0, 0, -(left + right) / rl,
797                 0, 2.0 / tb, 0, -(top + bottom) / tb,
798                 0, 0, -2.0 / fn, -(far + near) / fn,
799                 0, 0,  0, 1.0f
800             );
801         }
802         
803         /// Constructs lookAt matrix
804         static MatType lookAt(vec!(T, 3) eye, vec!(T, 3) target, vec!(T, 3) up) {
805             vec3 vz = eye - target;
806             vz.normalize;
807 
808             vec3 vx = up.cross(vz);
809             vx.normalize;
810 
811             vec3 vy = vz.cross(vx);
812 
813             return MatType(
814                 vx.x, vx.y, vx.z, -vx.dot(eye), 
815                 vy.x, vy.y, vy.z, -vy.dot(eye),
816                 vz.x, vz.y, vz.z, -vz.dot(eye),
817                 0, 0, 0, 1
818             );
819         }
820     }
821     
822     /// Removes negative sign from numbers less then float.epsilon*2
823     public MatType fixNegativeZero() {
824         foreach (i; 0 .. H) {
825             foreach (j; 0..W) {
826                 if (abs(data[i][j]) <= float.epsilon * 2.0f) data[i][j] = 0;
827             }
828         }
829         return this;
830     }
831 
832     /++
833     Construct row-major one dimentional array from matrix. 
834 
835     **Do not pass to OpenGL, use** `glbuffer()` **instead**
836     +/
837     public T[W*H] buffer() {
838         T[W*H] ret;
839         foreach (i; 0 .. H) {
840             foreach (j; 0..W) {
841                 ret[i * W + j] = data[i][j];
842             }
843         }
844         return ret;
845     }
846     /// Ditto
847     alias toArray = buffer;
848 
849     /++
850     Construct column-major one dimentional array from matrix.
851 
852     **OpenGL is column-major, please use this method for it.**
853     +/
854     public T[W*H] glbuffer() {
855         T[W*H] ret;
856         foreach (j; 0..W) {
857             foreach (i; 0 .. H) {
858                 ret[j * H + i] = data[i][j];
859             }
860         }
861         return ret;
862     }
863     /// Ditto
864     alias toGlArray = glbuffer;
865 }
866 
867 /// Is V a matrix with any size and any type
868 template isMatrix(M) {
869     enum isMatrix = is(Unqual!M == Matrix!(T, H, W), T, size_t H, size_t W);
870 }
871 
872 /// Is V a vector with size H2xW2 and any type
873 template isMatrix(M, size_t H2, size_t W2) {
874     static if(is(Unqual!M == Matrix!(T, H, W), T, size_t H, size_t W)) {
875         enum isMatrix = H2 == H && W2 == W;
876     } else {
877         enum isMatrix = false;
878     }
879 }
880 
881 /// Is V a vector with any size and type T
882 template isMatrix(M, T) {
883     static if(is(Unqual!M == Matrix!(U, H, W), U, size_t H, size_t W)) {
884         enum isMatrix = is(T == U);
885     } else {
886         enum isMatrix = false;
887     }
888 }
889 
890 /// Is V a vector with size H2xW2 and type T
891 template isMatrix(M, T, size_t H, size_t W) {
892     enum isMatrix = is(Unqual!M == Matrix!(T, H, W));
893 }
894