1 module dlangui.core.math3d;
2 
3 import std.math;
4 import std.string : format;
5 
6 /// 2 dimensional vector
7 struct vec2 {
8     union {
9         float[2] vec;
10         struct {
11             float x;
12             float y;
13         }
14     }
15     alias u = x;
16     alias v = y;
17     /// create with all components filled with specified value
18     this(float v) {
19         x = v;
20         y = v;
21     }
22     this(float[2] v) {
23         vec = v;
24     }
25     this(float[] v) {
26         vec = v[0..2];
27     }
28     this(float * v) {
29         vec = v[0..2];
30     }
31     this(const vec2 v) {
32         vec = v.vec;
33     }
34     this(float x, float y) {
35         vec[0] = x;
36         vec[1] = y;
37     }
38     ref vec2 opAssign(float[2] v) return {
39         vec = v;
40         return this;
41     }
42     ref vec2 opAssign(vec2 v) return {
43         vec = v.vec;
44         return this;
45     }
46     /// fill all components of vector with specified value
47     ref vec2 clear(float v) return {
48         vec[0] = vec[1] = v;
49         return this;
50     }
51     /// add value to all components of vector
52     ref vec2 add(float v) return {
53         vec[0] += v;
54         vec[1] += v;
55         return this;
56     }
57     /// multiply all components of vector by value
58     ref vec2 mul(float v) return {
59         vec[0] *= v;
60         vec[1] *= v;
61         return this;
62     }
63     /// subtract value from all components of vector
64     ref vec2 sub(float v) return {
65         vec[0] -= v;
66         vec[1] -= v;
67         return this;
68     }
69     /// divide all components of vector by value
70     ref vec2 div(float v) return {
71         vec[0] /= v;
72         vec[1] /= v;
73         return this;
74     }
75     /// add components of another vector to corresponding components of this vector
76     ref vec2 add(vec2 v) return {
77         vec[0] += v.vec[0];
78         vec[1] += v.vec[1];
79         return this;
80     }
81     /// multiply components of this vector  by corresponding components of another vector
82     ref vec2 mul(vec2 v) return {
83         vec[0] *= v.vec[0];
84         vec[1] *= v.vec[1];
85         return this;
86     }
87     /// subtract components of another vector from corresponding components of this vector
88     ref vec2 sub(vec2 v) return {
89         vec[0] -= v.vec[0];
90         vec[1] -= v.vec[1];
91         return this;
92     }
93     /// divide components of this vector  by corresponding components of another vector
94     ref vec2 div(vec2 v) return {
95         vec[0] /= v.vec[0];
96         vec[1] /= v.vec[1];
97         return this;
98     }
99 
100     /// returns vector rotated 90 degrees counter clockwise
101     vec2 rotated90ccw() const {
102         return vec2(-y, x);
103     }
104 
105     /// returns vector rotated 90 degrees clockwise
106     vec2 rotated90cw() const {
107         return vec2(y, -x);
108     }
109 
110     /// add value to all components of vector
111     vec2 opBinary(string op : "+")(float v) const {
112         vec2 res = this;
113         res.vec[0] += v;
114         res.vec[1] += v;
115         return res;
116     }
117     /// multiply all components of vector by value
118     vec2 opBinary(string op : "*")(float v) const {
119         vec2 res = this;
120         res.vec[0] *= v;
121         res.vec[1] *= v;
122         return res;
123     }
124     /// subtract value from all components of vector
125     vec2 opBinary(string op : "-")(float v) const {
126         vec2 res = this;
127         res.vec[0] -= v;
128         res.vec[1] -= v;
129         return res;
130     }
131     /// divide all components of vector by value
132     vec2 opBinary(string op : "/")(float v) const {
133         vec2 res = this;
134         res.vec[0] /= v;
135         res.vec[1] /= v;
136         return res;
137     }
138 
139 
140     /// add value to all components of vector
141     ref vec2 opOpAssign(string op : "+")(float v) return {
142         vec[0] += v;
143         vec[1] += v;
144         return this;
145     }
146     /// multiply all components of vector by value
147     ref vec2 opOpAssign(string op : "*")(float v) return {
148         vec[0] *= v;
149         vec[1] *= v;
150         return this;
151     }
152     /// subtract value from all components of vector
153     ref vec2 opOpAssign(string op : "-")(float v) return {
154         vec[0] -= v;
155         vec[1] -= v;
156         vec[2] -= v;
157         return this;
158     }
159     /// divide all components of vector by value
160     ref vec2 opOpAssign(string op : "/")(float v) return {
161         vec[0] /= v;
162         vec[1] /= v;
163         vec[2] /= v;
164         return this;
165     }
166 
167     /// by component add values of corresponding components of other vector
168     ref vec2 opOpAssign(string op : "+")(const vec2 v) return {
169         vec[0] += v.vec[0];
170         vec[1] += v.vec[1];
171         return this;
172     }
173     /// by component multiply values of corresponding components of other vector
174     ref vec2 opOpAssign(string op : "*")(const vec2 v) return {
175         vec[0] *= v.vec[0];
176         vec[1] *= v.vec[1];
177         return this;
178     }
179     /// by component subtract values of corresponding components of other vector
180     ref vec2 opOpAssign(string op : "-")(const vec2 v) return {
181         vec[0] -= v.vec[0];
182         vec[1] -= v.vec[1];
183         return this;
184     }
185     /// by component divide values of corresponding components of other vector
186     ref vec2 opOpAssign(string op : "/")(const vec2 v) return {
187         vec[0] /= v.vec[0];
188         vec[1] /= v.vec[1];
189         return this;
190     }
191 
192 
193     /// add value to all components of vector
194     vec2 opBinary(string op : "+")(const vec2 v) const {
195         vec2 res = this;
196         res.vec[0] += v.vec[0];
197         res.vec[1] += v.vec[1];
198         return res;
199     }
200     /// subtract value from all components of vector
201     vec2 opBinary(string op : "-")(const vec2 v) const {
202         vec2 res = this;
203         res.vec[0] -= v.vec[0];
204         res.vec[1] -= v.vec[1];
205         return res;
206     }
207     /// subtract value from all components of vector
208     float opBinary(string op : "*")(const vec3 v) const {
209         return dot(v);
210     }
211     /// dot product (sum of by-component products of vector components)
212     float dot(const vec2 v) const {
213         float res = 0.0f;
214         res += vec[0] * v.vec[0];
215         res += vec[1] * v.vec[1];
216         return res;
217     }
218 
219     /// cross product of 2 vec2 is scalar in Z axis
220     float crossProduct(const vec2 v2) const {
221         return x * v2.y - y * v2.x;
222     }
223 
224     /// returns vector with all components which are negative of components for this vector
225     vec2 opUnary(string op : "-")() const {
226         vec2 ret = this;
227         ret.vec[0] = -vec[0];
228         ret.vec[1] = -vec[1];
229         return ret;
230     }
231 
232 
233     /// sum of squares of all vector components
234     @property float magnitudeSquared() {
235         return vec[0]*vec[0] + vec[1]*vec[1];
236     }
237 
238     /// length of vector
239     @property float magnitude() {
240         return sqrt(magnitudeSquared);
241     }
242 
243     alias length = magnitude;
244 
245     /// normalize vector: make its length == 1
246     void normalize() {
247         div(length);
248     }
249 
250     /// returns normalized copy of this vector
251     @property vec2 normalized() {
252         vec2 res = this;
253         res.normalize();
254         return res;
255     }
256 }
257 
258 /// 3 dimensional vector
259 struct vec3 {
260     union {
261         float[3] vec;
262         struct {
263             float x;
264             float y;
265             float z;
266         }
267     }
268     //@property ref float x() { return vec[0]; }
269     //@property ref float y() { return vec[1]; }
270     //@property ref float z() { return vec[2]; }
271     alias r = x;
272     alias g = y;
273     alias b = z;
274     /// create with all components filled with specified value
275     this(float v) {
276         x = y = z = v;
277     }
278     this(float[3] v) {
279         vec = v;
280     }
281     this(float[] v) {
282         vec = v[0..3];
283     }
284     this(float * v) {
285         vec = v[0..3];
286     }
287     this(const vec3 v) {
288         vec = v.vec;
289     }
290     this(float x, float y, float z) {
291         vec[0] = x;
292         vec[1] = y;
293         vec[2] = z;
294     }
295     ref vec3 opAssign(float[3] v) return {
296         vec = v;
297         return this;
298     }
299     ref vec3 opAssign(vec3 v) return {
300         vec = v.vec;
301         return this;
302     }
303     ref vec3 opAssign(float x, float y, float z) return {
304         vec[0] = x;
305         vec[1] = y;
306         vec[2] = z;
307         return this;
308     }
309     /// fill all components of vector with specified value
310     ref vec3 clear(float v) return {
311         vec[0] = vec[1] = vec[2] = v;
312         return this;
313     }
314     /// add value to all components of vector
315     ref vec3 add(float v) return {
316         vec[0] += v;
317         vec[1] += v;
318         vec[2] += v;
319         return this;
320     }
321     /// multiply all components of vector by value
322     ref vec3 mul(float v) return {
323         vec[0] *= v;
324         vec[1] *= v;
325         vec[2] *= v;
326         return this;
327     }
328     /// subtract value from all components of vector
329     ref vec3 sub(float v) return {
330         vec[0] -= v;
331         vec[1] -= v;
332         vec[2] -= v;
333         return this;
334     }
335     /// divide all components of vector by value
336     ref vec3 div(float v) return {
337         vec[0] /= v;
338         vec[1] /= v;
339         vec[2] /= v;
340         return this;
341     }
342     /// add components of another vector to corresponding components of this vector
343     ref vec3 add(vec3 v) return {
344         vec[0] += v.vec[0];
345         vec[1] += v.vec[1];
346         vec[2] += v.vec[2];
347         return this;
348     }
349     /// multiply components of this vector  by corresponding components of another vector
350     ref vec3 mul(vec3 v) return {
351         vec[0] *= v.vec[0];
352         vec[1] *= v.vec[1];
353         vec[2] *= v.vec[2];
354         return this;
355     }
356     /// subtract components of another vector from corresponding components of this vector
357     ref vec3 sub(vec3 v) return {
358         vec[0] -= v.vec[0];
359         vec[1] -= v.vec[1];
360         vec[2] -= v.vec[2];
361         return this;
362     }
363     /// divide components of this vector  by corresponding components of another vector
364     ref vec3 div(vec3 v) return {
365         vec[0] /= v.vec[0];
366         vec[1] /= v.vec[1];
367         vec[2] /= v.vec[2];
368         return this;
369     }
370 
371     /// add value to all components of vector
372     vec3 opBinary(string op : "+")(float v) const {
373         vec3 res = this;
374         res.vec[0] += v;
375         res.vec[1] += v;
376         res.vec[2] += v;
377         return res;
378     }
379     /// multiply all components of vector by value
380     vec3 opBinary(string op : "*")(float v) const {
381         vec3 res = this;
382         res.vec[0] *= v;
383         res.vec[1] *= v;
384         res.vec[2] *= v;
385         return res;
386     }
387     /// subtract value from all components of vector
388     vec3 opBinary(string op : "-")(float v) const {
389         vec3 res = this;
390         res.vec[0] -= v;
391         res.vec[1] -= v;
392         res.vec[2] -= v;
393         return res;
394     }
395     /// divide all components of vector by value
396     vec3 opBinary(string op : "/")(float v) const {
397         vec3 res = this;
398         res.vec[0] /= v;
399         res.vec[1] /= v;
400         res.vec[2] /= v;
401         return res;
402     }
403 
404 
405     /// add value to all components of vector
406     ref vec3 opOpAssign(string op : "+")(float v) {
407         vec[0] += v;
408         vec[1] += v;
409         vec[2] += v;
410         return this;
411     }
412     /// multiply all components of vector by value
413     ref vec3 opOpAssign(string op : "*")(float v) {
414         vec[0] *= v;
415         vec[1] *= v;
416         vec[2] *= v;
417         return this;
418     }
419     /// subtract value from all components of vector
420     ref vec3 opOpAssign(string op : "-")(float v) {
421         vec[0] -= v;
422         vec[1] -= v;
423         vec[2] -= v;
424         return this;
425     }
426     /// divide all components of vector by value
427     ref vec3 opOpAssign(string op : "/")(float v) {
428         vec[0] /= v;
429         vec[1] /= v;
430         vec[2] /= v;
431         return this;
432     }
433 
434     /// by component add values of corresponding components of other vector
435     ref vec3 opOpAssign(string op : "+")(const vec3 v) {
436         vec[0] += v.vec[0];
437         vec[1] += v.vec[1];
438         vec[2] += v.vec[2];
439         return this;
440     }
441     /// by component multiply values of corresponding components of other vector
442     ref vec3 opOpAssign(string op : "*")(const vec3 v) {
443         vec[0] *= v.vec[0];
444         vec[1] *= v.vec[1];
445         vec[2] *= v.vec[2];
446         return this;
447     }
448     /// by component subtract values of corresponding components of other vector
449     ref vec3 opOpAssign(string op : "-")(const vec3 v) {
450         vec[0] -= v.vec[0];
451         vec[1] -= v.vec[1];
452         vec[2] -= v.vec[2];
453         return this;
454     }
455     /// by component divide values of corresponding components of other vector
456     ref vec3 opOpAssign(string op : "/")(const vec3 v) {
457         vec[0] /= v.vec[0];
458         vec[1] /= v.vec[1];
459         vec[2] /= v.vec[2];
460         return this;
461     }
462 
463 
464     /// add value to all components of vector
465     vec3 opBinary(string op : "+")(const vec3 v) const {
466         vec3 res = this;
467         res.vec[0] += v.vec[0];
468         res.vec[1] += v.vec[1];
469         res.vec[2] += v.vec[2];
470         return res;
471     }
472     /// subtract value from all components of vector
473     vec3 opBinary(string op : "-")(const vec3 v) const {
474         vec3 res = this;
475         res.vec[0] -= v.vec[0];
476         res.vec[1] -= v.vec[1];
477         res.vec[2] -= v.vec[2];
478         return res;
479     }
480     /// subtract value from all components of vector
481     float opBinary(string op : "*")(const vec3 v) const {
482         return dot(v);
483     }
484     /// dot product (sum of by-component products of vector components)
485     float dot(const vec3 v) const {
486         float res = 0.0f;
487         res += vec[0] * v.vec[0];
488         res += vec[1] * v.vec[1];
489         res += vec[2] * v.vec[2];
490         return res;
491     }
492 
493     /// returns vector with all components which are negative of components for this vector
494     vec3 opUnary(string op : "-")() const {
495         vec3 ret = this;
496         ret.vec[0] = -vec[0];
497         ret.vec[1] = -vec[1];
498         ret.vec[2] = -vec[2];
499         return ret;
500     }
501 
502 
503     /// sum of squares of all vector components
504     @property float magnitudeSquared() {
505         return vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2];
506     }
507 
508     /// length of vector
509     @property float magnitude() {
510         return sqrt(magnitudeSquared);
511     }
512 
513     alias length = magnitude;
514 
515     /// normalize vector: make its length == 1
516     void normalize() {
517         div(length);
518     }
519 
520     /// returns normalized copy of this vector
521     @property vec3 normalized() {
522         vec3 res = this;
523         res.normalize();
524         return res;
525     }
526 
527     /// cross product
528     static vec3 crossProduct(const vec3 v1, const vec3 v2) {
529         return vec3(v1.y * v2.z - v1.z * v2.y,
530                     v1.z * v2.x - v1.x * v2.z,
531                     v1.x * v2.y - v1.y * v2.x);
532     }
533 
534     /// multiply vector by matrix
535     vec3 opBinary(string op : "*")(const ref mat4 matrix) const
536     {
537         float xx, yy, zz, ww;
538         xx = x * matrix.m[0*4 + 0] +
539             y * matrix.m[0*4 + 1] +
540             z * matrix.m[0*4 + 2] +
541             matrix.m[0*4 + 3];
542         yy = x * matrix.m[1*4 + 0] +
543             y * matrix.m[1*4 + 1] +
544             z * matrix.m[1*4 + 2] +
545             matrix.m[1*4 + 3];
546         zz = x * matrix.m[2*4 + 0] +
547             y * matrix.m[2*4 + 1] +
548             z * matrix.m[2*4 + 2] +
549             matrix.m[2*4 + 3];
550         ww = x * matrix.m[3*4 + 0] +
551             y * matrix.m[3*4 + 1] +
552             z * matrix.m[3*4 + 2] +
553             matrix.m[3*4 + 3];
554         if (ww == 1.0f)
555             return vec3(xx, yy, zz);
556         else
557             return vec3(xx / ww, yy / ww, zz / ww);
558     }
559 
560     @property string toString() {
561         return "(%f,%f,%f)".format(x, y, z);
562     }
563 }
564 
565 /// 4 component vector
566 struct vec4 {
567     union {
568         float[4] vec;
569         struct {
570             float x;
571             float y;
572             float z;
573             float w;
574         }
575     }
576     alias r = x;
577     alias g = y;
578     alias b = z;
579     alias a = w;
580     /// create with all components filled with specified value
581     this(float v) {
582         x = y = z = w = v;
583     }
584     this(float[4] v) {
585         vec = v;
586     }
587     this(vec4 v) {
588         vec = v.vec;
589     }
590     this(float x, float y, float z, float w) {
591         vec[0] = x;
592         vec[1] = y;
593         vec[2] = z;
594         vec[3] = w;
595     }
596     this(vec3 v) {
597         vec[0] = v.vec[0];
598         vec[1] = v.vec[1];
599         vec[2] = v.vec[2];
600         vec[3] = 1.0f;
601     }
602     ref vec4 opAssign(const float[4] v) return {
603         vec = v;
604         return this;
605     }
606     ref vec4 opAssign(const vec4 v) return {
607         vec = v.vec;
608         return this;
609     }
610     ref vec4 opAssign(float x, float y, float z, float w) return {
611         vec[0] = x;
612         vec[1] = y;
613         vec[2] = z;
614         vec[3] = w;
615         return this;
616     }
617     ref vec4 opAssign(const vec3 v) return {
618         vec[0] = v.vec[0];
619         vec[1] = v.vec[1];
620         vec[2] = v.vec[2];
621         vec[3] = 1.0f;
622         return this;
623     }
624 
625 
626     /// fill all components of vector with specified value
627     ref vec4 clear(float v) return {
628         vec[0] = vec[1] = vec[2] = vec[3] = v;
629         return this;
630     }
631     /// add value to all components of vector
632     ref vec4 add(float v) return {
633         vec[0] += v;
634         vec[1] += v;
635         vec[2] += v;
636         vec[3] += v;
637         return this;
638     }
639     /// multiply all components of vector by value
640     ref vec4 mul(float v) return {
641         vec[0] *= v;
642         vec[1] *= v;
643         vec[2] *= v;
644         vec[3] *= v;
645         return this;
646     }
647     /// subtract value from all components of vector
648     ref vec4 sub(float v) return {
649         vec[0] -= v;
650         vec[1] -= v;
651         vec[2] -= v;
652         vec[3] -= v;
653         return this;
654     }
655     /// divide all components of vector by value
656     ref vec4 div(float v) return {
657         vec[0] /= v;
658         vec[1] /= v;
659         vec[2] /= v;
660         vec[3] /= v;
661         return this;
662     }
663     /// add components of another vector to corresponding components of this vector
664     ref vec4 add(const vec4 v) return {
665         vec[0] += v.vec[0];
666         vec[1] += v.vec[1];
667         vec[2] += v.vec[2];
668         vec[3] += v.vec[3];
669         return this;
670     }
671     /// multiply components of this vector  by corresponding components of another vector
672     ref vec4 mul(vec4 v) return {
673         vec[0] *= v.vec[0];
674         vec[1] *= v.vec[1];
675         vec[2] *= v.vec[2];
676         vec[3] *= v.vec[3];
677         return this;
678     }
679     /// subtract components of another vector from corresponding components of this vector
680     ref vec4 sub(vec4 v) return {
681         vec[0] -= v.vec[0];
682         vec[1] -= v.vec[1];
683         vec[2] -= v.vec[2];
684         vec[3] -= v.vec[3];
685         return this;
686     }
687     /// divide components of this vector  by corresponding components of another vector
688     ref vec4 div(vec4 v) return {
689         vec[0] /= v.vec[0];
690         vec[1] /= v.vec[1];
691         vec[2] /= v.vec[2];
692         vec[3] /= v.vec[3];
693         return this;
694     }
695 
696     /// add value to all components of vector
697     vec4 opBinary(string op : "+")(float v) const return {
698         vec4 res = this;
699         res.vec[0] += v;
700         res.vec[1] += v;
701         res.vec[2] += v;
702         res.vec[3] += v;
703         return res;
704     }
705     /// multiply all components of vector by value
706     vec4 opBinary(string op : "*")(float v) const return {
707         vec4 res = this;
708         res.vec[0] *= v;
709         res.vec[1] *= v;
710         res.vec[2] *= v;
711         res.vec[3] *= v;
712         return res;
713     }
714     /// subtract value from all components of vector
715     vec4 opBinary(string op : "-")(float v) const return {
716         vec4 res = this;
717         res.vec[0] -= v;
718         res.vec[1] -= v;
719         res.vec[2] -= v;
720         res.vec[3] -= v;
721         return res;
722     }
723     /// divide all components of vector by value
724     vec4 opBinary(string op : "/")(float v) const return {
725         vec4 res = this;
726         res.vec[0] /= v;
727         res.vec[1] /= v;
728         res.vec[2] /= v;
729         res.vec[3] /= v;
730         return res;
731     }
732 
733     /// add value to all components of vector
734     ref vec4 opOpAssign(string op : "+")(float v) return {
735         vec[0] += v;
736         vec[1] += v;
737         vec[2] += v;
738         vec[3] += v;
739         return this;
740     }
741     /// multiply all components of vector by value
742     ref vec4 opOpAssign(string op : "*")(float v) return {
743         vec[0] *= v;
744         vec[1] *= v;
745         vec[2] *= v;
746         vec[3] *= v;
747         return this;
748     }
749     /// subtract value from all components of vector
750     ref vec4 opOpAssign(string op : "-")(float v) return {
751         vec[0] -= v;
752         vec[1] -= v;
753         vec[2] -= v;
754         vec[3] -= v;
755         return this;
756     }
757     /// divide all components of vector by value
758     ref vec4 opOpAssign(string op : "/")(float v) return {
759         vec[0] /= v;
760         vec[1] /= v;
761         vec[2] /= v;
762         vec[3] /= v;
763         return this;
764     }
765 
766     /// by component add values of corresponding components of other vector
767     ref vec4 opOpAssign(string op : "+")(const vec4 v) return {
768         vec[0] += v.vec[0];
769         vec[1] += v.vec[1];
770         vec[2] += v.vec[2];
771         vec[3] += v.vec[3];
772         return this;
773     }
774     /// by component multiply values of corresponding components of other vector
775     ref vec4 opOpAssign(string op : "*")(const vec4 v) return {
776         vec[0] *= v.vec[0];
777         vec[1] *= v.vec[1];
778         vec[2] *= v.vec[2];
779         vec[3] *= v.vec[3];
780         return this;
781     }
782     /// by component subtract values of corresponding components of other vector
783     ref vec4 opOpAssign(string op : "-")(const vec4 v) return {
784         vec[0] -= v.vec[0];
785         vec[1] -= v.vec[1];
786         vec[2] -= v.vec[2];
787         vec[3] -= v.vec[3];
788         return this;
789     }
790     /// by component divide values of corresponding components of other vector
791     ref vec4 opOpAssign(string op : "/")(const vec4 v) return {
792         vec[0] /= v.vec[0];
793         vec[1] /= v.vec[1];
794         vec[2] /= v.vec[2];
795         vec[3] /= v.vec[3];
796         return this;
797     }
798 
799 
800 
801     /// add value to all components of vector
802     vec4 opBinary(string op : "+")(const vec4 v) const return {
803         vec4 res = this;
804         res.vec[0] += v.vec[0];
805         res.vec[1] += v.vec[1];
806         res.vec[2] += v.vec[2];
807         res.vec[3] += v.vec[3];
808         return res;
809     }
810     /// subtract value from all components of vector
811     vec4 opBinary(string op : "-")(const vec4 v) const return {
812         vec4 res = this;
813         res.vec[0] -= v.vec[0];
814         res.vec[1] -= v.vec[1];
815         res.vec[2] -= v.vec[2];
816         res.vec[3] -= v.vec[3];
817         return res;
818     }
819     /// subtract value from all components of vector
820     float opBinary(string op : "*")(const vec4 v) const {
821         return dot(v);
822     }
823     /// dot product (sum of by-component products of vector components)
824     float dot(vec4 v) const {
825         float res = 0.0f;
826         res += vec[0] * v.vec[0];
827         res += vec[1] * v.vec[1];
828         res += vec[2] * v.vec[2];
829         res += vec[3] * v.vec[3];
830         return res;
831     }
832 
833     /// returns vector with all components which are negative of components for this vector
834     vec4 opUnary(string op : "-")() const return {
835         vec4 ret = this;
836         ret[0] = -vec[0];
837         ret[1] = -vec[1];
838         ret[2] = -vec[2];
839         ret[3] = -vec[3];
840         return ret;
841     }
842 
843 
844 
845     /// sum of squares of all vector components
846     @property float magnitudeSquared() {
847         return vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2] + vec[3]*vec[3];
848     }
849 
850     /// length of vector
851     @property float magnitude() {
852         return sqrt(magnitudeSquared);
853     }
854 
855     alias length = magnitude;
856 
857     /// normalize vector: make its length == 1
858     void normalize() {
859         div(length);
860     }
861 
862     /// returns normalized copy of this vector
863     @property vec4 normalized() return {
864         vec4 res = this;
865         res.normalize();
866         return res;
867     }
868 
869     /// multiply vector by matrix
870     vec4 opBinary(string op : "*")(const ref mat4 matrix) const
871     {
872         float xx, yy, zz, ww;
873         xx = x * matrix.m[0*4 + 0] +
874              y * matrix.m[0*4 + 1] +
875              z * matrix.m[0*4 + 2] +
876              w * matrix.m[0*4 + 3];
877         yy = x * matrix.m[1*4 + 0] +
878              y * matrix.m[1*4 + 1] +
879              z * matrix.m[1*4 + 2] +
880              w * matrix.m[1*4 + 3];
881         zz = x * matrix.m[2*4 + 0] +
882              y * matrix.m[2*4 + 1] +
883              z * matrix.m[2*4 + 2] +
884              w * matrix.m[2*4 + 3];
885         ww = x * matrix.m[3*4 + 0] +
886              y * matrix.m[3*4 + 1] +
887              z * matrix.m[3*4 + 2] +
888              w * matrix.m[3*4 + 3];
889         return vec4(xx, yy, zz, ww);
890     }
891 
892     @property string toString() {
893         return "(%f,%f,%f,%f)".format(x, y, z, w);
894     }
895 }
896 
897 bool fuzzyNull(float v) {
898     return v < 0.0000001f && v > -0.0000001f;
899 }
900 
901 /// float matrix 4 x 4
902 struct mat4 {
903     float[16] m = [1, 0, 0, 0,  0, 1, 0, 0,  0, 0, 1, 0,  0, 0, 0, 1];
904 
905     @property string dump() const {
906         import std.conv : to;
907         return to!string(m[0..4]) ~ to!string(m[4..8]) ~ to!string(m[8..12]) ~ to!string(m[12..16]);
908     }
909 
910     //alias m this;
911 
912     this(float v) {
913         setDiagonal(v);
914     }
915 
916     this(const ref mat4 v) {
917         m[0..16] = v.m[0..16];
918     }
919     this(const float[16] v) {
920         m[0..16] = v[0..16];
921     }
922 
923     ref mat4 opAssign(const ref mat4 v) return {
924         m[0..16] = v.m[0..16];
925         return this;
926     }
927     ref mat4 opAssign(const  mat4 v) return {
928         m[0..16] = v.m[0..16];
929         return this;
930     }
931     ref mat4 opAssign(const float[16] v) return {
932         m[0..16] = v[0..16];
933         return this;
934     }
935 
936     void setOrtho(float left, float right, float bottom, float top, float nearPlane, float farPlane)
937     {
938         // Bail out if the projection volume is zero-sized.
939         if (left == right || bottom == top || nearPlane == farPlane)
940             return;
941 
942         // Construct the projection.
943         float width = right - left;
944         float invheight = top - bottom;
945         float clip = farPlane - nearPlane;
946         m[0*4 + 0] = 2.0f / width;
947         m[1*4 + 0] = 0.0f;
948         m[2*4 + 0] = 0.0f;
949         m[3*4 + 0] = -(left + right) / width;
950         m[0*4 + 1] = 0.0f;
951         m[1*4 + 1] = 2.0f / invheight;
952         m[2*4 + 1] = 0.0f;
953         m[3*4 + 1] = -(top + bottom) / invheight;
954         m[0*4 + 2] = 0.0f;
955         m[1*4 + 2] = 0.0f;
956         m[2*4 + 2] = -2.0f / clip;
957         m[3*4 + 2] = -(nearPlane + farPlane) / clip;
958         m[0*4 + 3] = 0.0f;
959         m[1*4 + 3] = 0.0f;
960         m[2*4 + 3] = 0.0f;
961         m[3*4 + 3] = 1.0f;
962     }
963 
964     void setPerspective(float angle, float aspect, float nearPlane, float farPlane)
965     {
966         // Bail out if the projection volume is zero-sized.
967         float radians = (angle / 2.0f) * PI / 180.0f;
968         if (nearPlane == farPlane || aspect == 0.0f || radians < 0.0001f)
969             return;
970         float f = 1 / tan(radians);
971         float d = 1 / (nearPlane - farPlane);
972 
973         // Construct the projection.
974         m[0*4 + 0] = f / aspect;
975         m[1*4 + 0] = 0.0f;
976         m[2*4 + 0] = 0.0f;
977         m[3*4 + 0] = 0.0f;
978 
979         m[0*4 + 1] = 0.0f;
980         m[1*4 + 1] = f;
981         m[2*4 + 1] = 0.0f;
982         m[3*4 + 1] = 0.0f;
983 
984         m[0*4 + 2] = 0.0f;
985         m[1*4 + 2] = 0.0f;
986         m[2*4 + 2] = (nearPlane + farPlane) * d;
987         m[3*4 + 2] = 2.0f * nearPlane * farPlane * d;
988 
989         m[0*4 + 3] = 0.0f;
990         m[1*4 + 3] = 0.0f;
991         m[2*4 + 3] = -1.0f;
992         m[3*4 + 3] = 0.0f;
993     }
994 
995     ref mat4 lookAt(const vec3 eye, const vec3 center, const vec3 up) return {
996         vec3 forward = (center - eye).normalized();
997         vec3 side = vec3.crossProduct(forward, up).normalized();
998         vec3 upVector = vec3.crossProduct(side, forward);
999 
1000         mat4 m;
1001         m.setIdentity();
1002         m[0*4 + 0] = side.x;
1003         m[1*4 + 0] = side.y;
1004         m[2*4 + 0] = side.z;
1005         m[3*4 + 0] = 0.0f;
1006         m[0*4 + 1] = upVector.x;
1007         m[1*4 + 1] = upVector.y;
1008         m[2*4 + 1] = upVector.z;
1009         m[3*4 + 1] = 0.0f;
1010         m[0*4 + 2] = -forward.x;
1011         m[1*4 + 2] = -forward.y;
1012         m[2*4 + 2] = -forward.z;
1013         m[3*4 + 2] = 0.0f;
1014         m[0*4 + 3] = 0.0f;
1015         m[1*4 + 3] = 0.0f;
1016         m[2*4 + 3] = 0.0f;
1017         m[3*4 + 3] = 1.0f;
1018 
1019         this *= m;
1020         translate(-eye);
1021         return this;
1022     }
1023 
1024     /// transpose matrix
1025     void transpose() {
1026         float[16] tmp = [
1027             m[0], m[4], m[8], m[12],
1028             m[1], m[5], m[9], m[13],
1029             m[2], m[6], m[10], m[14],
1030             m[3], m[7], m[11], m[15]
1031         ];
1032         m = tmp;
1033     }
1034 
1035     mat4 invert() const
1036     {
1037         float a0 = m[0] * m[5] - m[1] * m[4];
1038         float a1 = m[0] * m[6] - m[2] * m[4];
1039         float a2 = m[0] * m[7] - m[3] * m[4];
1040         float a3 = m[1] * m[6] - m[2] * m[5];
1041         float a4 = m[1] * m[7] - m[3] * m[5];
1042         float a5 = m[2] * m[7] - m[3] * m[6];
1043         float b0 = m[8] * m[13] - m[9] * m[12];
1044         float b1 = m[8] * m[14] - m[10] * m[12];
1045         float b2 = m[8] * m[15] - m[11] * m[12];
1046         float b3 = m[9] * m[14] - m[10] * m[13];
1047         float b4 = m[9] * m[15] - m[11] * m[13];
1048         float b5 = m[10] * m[15] - m[11] * m[14];
1049 
1050         // Calculate the determinant.
1051         float det = a0 * b5 - a1 * b4 + a2 * b3 + a3 * b2 - a4 * b1 + a5 * b0;
1052 
1053         mat4 inverse;
1054 
1055         // Close to zero, can't invert.
1056         if (fabs(det) <= 0.00000001f)
1057             return inverse;
1058 
1059         // Support the case where m == dst.
1060         inverse.m[0]  = m[5] * b5 - m[6] * b4 + m[7] * b3;
1061         inverse.m[1]  = -m[1] * b5 + m[2] * b4 - m[3] * b3;
1062         inverse.m[2]  = m[13] * a5 - m[14] * a4 + m[15] * a3;
1063         inverse.m[3]  = -m[9] * a5 + m[10] * a4 - m[11] * a3;
1064 
1065         inverse.m[4]  = -m[4] * b5 + m[6] * b2 - m[7] * b1;
1066         inverse.m[5]  = m[0] * b5 - m[2] * b2 + m[3] * b1;
1067         inverse.m[6]  = -m[12] * a5 + m[14] * a2 - m[15] * a1;
1068         inverse.m[7]  = m[8] * a5 - m[10] * a2 + m[11] * a1;
1069 
1070         inverse.m[8]  = m[4] * b4 - m[5] * b2 + m[7] * b0;
1071         inverse.m[9]  = -m[0] * b4 + m[1] * b2 - m[3] * b0;
1072         inverse.m[10] = m[12] * a4 - m[13] * a2 + m[15] * a0;
1073         inverse.m[11] = -m[8] * a4 + m[9] * a2 - m[11] * a0;
1074 
1075         inverse.m[12] = -m[4] * b3 + m[5] * b1 - m[6] * b0;
1076         inverse.m[13] = m[0] * b3 - m[1] * b1 + m[2] * b0;
1077         inverse.m[14] = -m[12] * a3 + m[13] * a1 - m[14] * a0;
1078         inverse.m[15] = m[8] * a3 - m[9] * a1 + m[10] * a0;
1079 
1080         float mul = 1.0f / det;
1081         inverse *= mul;
1082         return inverse;
1083     }
1084 
1085     ref mat4 setLookAt(const vec3 eye, const vec3 center, const vec3 up) return {
1086         setIdentity();
1087         lookAt(eye, center, up);
1088         return this;
1089     }
1090 
1091     ref mat4 translate(const vec3 v) return {
1092         m[3*4 + 0] += m[0*4 + 0] * v.x + m[1*4 + 0] * v.y + m[2*4 + 0] * v.z;
1093         m[3*4 + 1] += m[0*4 + 1] * v.x + m[1*4 + 1] * v.y + m[2*4 + 1] * v.z;
1094         m[3*4 + 2] += m[0*4 + 2] * v.x + m[1*4 + 2] * v.y + m[2*4 + 2] * v.z;
1095         m[3*4 + 3] += m[0*4 + 3] * v.x + m[1*4 + 3] * v.y + m[2*4 + 3] * v.z;
1096         return this;
1097     }
1098 
1099     ref mat4 translate(float x, float y, float z) return {
1100         m[3*4 + 0] += m[0*4 + 0] * x + m[1*4 + 0] * y + m[2*4 + 0] * z;
1101         m[3*4 + 1] += m[0*4 + 1] * x + m[1*4 + 1] * y + m[2*4 + 1] * z;
1102         m[3*4 + 2] += m[0*4 + 2] * x + m[1*4 + 2] * y + m[2*4 + 2] * z;
1103         m[3*4 + 3] += m[0*4 + 3] * x + m[1*4 + 3] * y + m[2*4 + 3] * z;
1104         return this;
1105     }
1106 
1107     /// add scalar to all items of matrix
1108     mat4 opBinary(string op : "+")(float v) const {
1109         foreach(ref item; m)
1110             item += v;
1111     }
1112 
1113     /// multiply this matrix by scalar
1114     mat4 opBinary(string op : "-")(float v) const {
1115         foreach(ref item; m)
1116             item -= v;
1117     }
1118 
1119     /// multiply this matrix by scalar
1120     mat4 opBinary(string op : "*")(float v) const {
1121         foreach(ref item; m)
1122             item *= v;
1123     }
1124 
1125     /// multiply this matrix by scalar
1126     mat4 opBinary(string op : "/")(float v) const {
1127         foreach(ref item; m)
1128             item /= v;
1129     }
1130 
1131     /// multiply this matrix by another matrix
1132     mat4 opBinary(string op : "*")(const ref mat4 m2) const {
1133         return mul(this, m2);
1134     }
1135 
1136     /// multiply this matrix by another matrix
1137     void opOpAssign(string op : "*")(const ref mat4 m2) {
1138         this = mul(this, m2);
1139     }
1140 
1141     /// multiply two matrices
1142     static mat4 mul(const ref mat4 m1, const ref mat4 m2) {
1143         mat4 m;
1144         m.m[0*4 + 0] =
1145             m1.m[0*4 + 0] * m2.m[0*4 + 0] +
1146             m1.m[1*4 + 0] * m2.m[0*4 + 1] +
1147             m1.m[2*4 + 0] * m2.m[0*4 + 2] +
1148             m1.m[3*4 + 0] * m2.m[0*4 + 3];
1149         m.m[0*4 + 1] =
1150             m1.m[0*4 + 1] * m2.m[0*4 + 0] +
1151             m1.m[1*4 + 1] * m2.m[0*4 + 1] +
1152             m1.m[2*4 + 1] * m2.m[0*4 + 2] +
1153             m1.m[3*4 + 1] * m2.m[0*4 + 3];
1154         m.m[0*4 + 2] =
1155             m1.m[0*4 + 2] * m2.m[0*4 + 0] +
1156             m1.m[1*4 + 2] * m2.m[0*4 + 1] +
1157             m1.m[2*4 + 2] * m2.m[0*4 + 2] +
1158             m1.m[3*4 + 2] * m2.m[0*4 + 3];
1159         m.m[0*4 + 3] =
1160             m1.m[0*4 + 3] * m2.m[0*4 + 0] +
1161             m1.m[1*4 + 3] * m2.m[0*4 + 1] +
1162             m1.m[2*4 + 3] * m2.m[0*4 + 2] +
1163             m1.m[3*4 + 3] * m2.m[0*4 + 3];
1164         m.m[1*4 + 0] =
1165             m1.m[0*4 + 0] * m2.m[1*4 + 0] +
1166             m1.m[1*4 + 0] * m2.m[1*4 + 1] +
1167             m1.m[2*4 + 0] * m2.m[1*4 + 2] +
1168             m1.m[3*4 + 0] * m2.m[1*4 + 3];
1169         m.m[1*4 + 1] =
1170             m1.m[0*4 + 1] * m2.m[1*4 + 0] +
1171             m1.m[1*4 + 1] * m2.m[1*4 + 1] +
1172             m1.m[2*4 + 1] * m2.m[1*4 + 2] +
1173             m1.m[3*4 + 1] * m2.m[1*4 + 3];
1174         m.m[1*4 + 2] =
1175             m1.m[0*4 + 2] * m2.m[1*4 + 0] +
1176             m1.m[1*4 + 2] * m2.m[1*4 + 1] +
1177             m1.m[2*4 + 2] * m2.m[1*4 + 2] +
1178             m1.m[3*4 + 2] * m2.m[1*4 + 3];
1179         m.m[1*4 + 3] =
1180             m1.m[0*4 + 3] * m2.m[1*4 + 0] +
1181             m1.m[1*4 + 3] * m2.m[1*4 + 1] +
1182             m1.m[2*4 + 3] * m2.m[1*4 + 2] +
1183             m1.m[3*4 + 3] * m2.m[1*4 + 3];
1184         m.m[2*4 + 0] =
1185             m1.m[0*4 + 0] * m2.m[2*4 + 0] +
1186             m1.m[1*4 + 0] * m2.m[2*4 + 1] +
1187             m1.m[2*4 + 0] * m2.m[2*4 + 2] +
1188             m1.m[3*4 + 0] * m2.m[2*4 + 3];
1189         m.m[2*4 + 1] =
1190             m1.m[0*4 + 1] * m2.m[2*4 + 0] +
1191             m1.m[1*4 + 1] * m2.m[2*4 + 1] +
1192             m1.m[2*4 + 1] * m2.m[2*4 + 2] +
1193             m1.m[3*4 + 1] * m2.m[2*4 + 3];
1194         m.m[2*4 + 2] =
1195             m1.m[0*4 + 2] * m2.m[2*4 + 0] +
1196             m1.m[1*4 + 2] * m2.m[2*4 + 1] +
1197             m1.m[2*4 + 2] * m2.m[2*4 + 2] +
1198             m1.m[3*4 + 2] * m2.m[2*4 + 3];
1199         m.m[2*4 + 3] =
1200             m1.m[0*4 + 3] * m2.m[2*4 + 0] +
1201             m1.m[1*4 + 3] * m2.m[2*4 + 1] +
1202             m1.m[2*4 + 3] * m2.m[2*4 + 2] +
1203             m1.m[3*4 + 3] * m2.m[2*4 + 3];
1204         m.m[3*4 + 0] =
1205             m1.m[0*4 + 0] * m2.m[3*4 + 0] +
1206             m1.m[1*4 + 0] * m2.m[3*4 + 1] +
1207             m1.m[2*4 + 0] * m2.m[3*4 + 2] +
1208             m1.m[3*4 + 0] * m2.m[3*4 + 3];
1209         m.m[3*4 + 1] =
1210             m1.m[0*4 + 1] * m2.m[3*4 + 0] +
1211             m1.m[1*4 + 1] * m2.m[3*4 + 1] +
1212             m1.m[2*4 + 1] * m2.m[3*4 + 2] +
1213             m1.m[3*4 + 1] * m2.m[3*4 + 3];
1214         m.m[3*4 + 2] =
1215             m1.m[0*4 + 2] * m2.m[3*4 + 0] +
1216             m1.m[1*4 + 2] * m2.m[3*4 + 1] +
1217             m1.m[2*4 + 2] * m2.m[3*4 + 2] +
1218             m1.m[3*4 + 2] * m2.m[3*4 + 3];
1219         m.m[3*4 + 3] =
1220             m1.m[0*4 + 3] * m2.m[3*4 + 0] +
1221             m1.m[1*4 + 3] * m2.m[3*4 + 1] +
1222             m1.m[2*4 + 3] * m2.m[3*4 + 2] +
1223             m1.m[3*4 + 3] * m2.m[3*4 + 3];
1224         return m;
1225     }
1226 
1227     /// multiply matrix by vec3
1228     vec3 opBinary(string op : "*")(const vec3 vector) const
1229     {
1230         float x, y, z, w;
1231         x = vector.x * m[0*4 + 0] +
1232             vector.y * m[1*4 + 0] +
1233             vector.z * m[2*4 + 0] +
1234             m[3*4 + 0];
1235         y = vector.x * m[0*4 + 1] +
1236             vector.y * m[1*4 + 1] +
1237             vector.z * m[2*4 + 1] +
1238             m[3*4 + 1];
1239         z = vector.x * m[0*4 + 2] +
1240             vector.y * m[1*4 + 2] +
1241             vector.z * m[2*4 + 2] +
1242             m[3*4 + 2];
1243         w = vector.x * m[0*4 + 3] +
1244             vector.y * m[1*4 + 3] +
1245             vector.z * m[2*4 + 3] +
1246             m[3*4 + 3];
1247         if (w == 1.0f)
1248             return vec3(x, y, z);
1249         else
1250             return vec3(x / w, y / w, z / w);
1251     }
1252 
1253     /// multiply matrix by vec4
1254     vec4 opBinary(string op : "*")(const vec4 vector) const
1255     {
1256         float x, y, z, w;
1257         x = vector.x * m[0*4 + 0] +
1258             vector.y * m[1*4 + 0] +
1259             vector.z * m[2*4 + 0] +
1260             vector.w * m[3*4 + 0];
1261         y = vector.x * m[0*4 + 1] +
1262             vector.y * m[1*4 + 1] +
1263             vector.z * m[2*4 + 1] +
1264             vector.w * m[3*4 + 1];
1265         z = vector.x * m[0*4 + 2] +
1266             vector.y * m[1*4 + 2] +
1267             vector.z * m[2*4 + 2] +
1268             vector.w * m[3*4 + 2];
1269         w = vector.x * m[0*4 + 3] +
1270             vector.y * m[1*4 + 3] +
1271             vector.z * m[2*4 + 3] +
1272             vector.w * m[3*4 + 3];
1273         return vec4(x, y, z, w);
1274     }
1275 
1276     /// 2d index by row, col
1277     ref float opIndex(int y, int x) return {
1278         return m[y*4 + x];
1279     }
1280 
1281     /// 2d index by row, col
1282     float opIndex(int y, int x) const return {
1283         return m[y*4 + x];
1284     }
1285 
1286     /// scalar index by rows then (y*4 + x)
1287     ref float opIndex(int index) return {
1288         return m[index];
1289     }
1290 
1291     /// scalar index by rows then (y*4 + x)
1292     float opIndex(int index) const return {
1293         return m[index];
1294     }
1295 
1296     /// set to identity: fill all items of matrix with zero except main diagonal items which will be assigned to 1.0f
1297     ref mat4 setIdentity() return {
1298         return setDiagonal(1.0f);
1299     }
1300     /// set to diagonal: fill all items of matrix with zero except main diagonal items which will be assigned to v
1301     ref mat4 setDiagonal(float v) return {
1302         for (int x = 0; x < 4; x++) {
1303             for (int y = 0; y < 4; y++) {
1304                 if (x == y)
1305                     m[y * 4 + x] = v;
1306                 else
1307                     m[y * 4 + x] = 0.0f;
1308             }
1309         }
1310         return this;
1311     }
1312     /// fill all items of matrix with specified value
1313     ref mat4 fill(float v) return {
1314         foreach(ref f; m)
1315             f = v;
1316         return this;
1317     }
1318     /// fill all items of matrix with zero
1319     ref mat4 setZero() return {
1320         foreach(ref f; m)
1321             f = 0.0f;
1322         return this;
1323     }
1324     /// creates identity matrix
1325     static mat4 identity() {
1326         mat4 res;
1327         return res.setIdentity();
1328     }
1329     /// creates zero matrix
1330     static mat4 zero() {
1331         mat4 res;
1332         return res.setZero();
1333     }
1334 
1335 
1336     /// add value to all components of matrix
1337     void opOpAssign(string op : "+")(float v) {
1338         foreach(ref item; m)
1339             item += v;
1340     }
1341     /// multiply all components of matrix by value
1342     void opOpAssign(string op : "*")(float v) {
1343         foreach(ref item; m)
1344             item *= v;
1345     }
1346     /// subtract value from all components of matrix
1347     void opOpAssign(string op : "-")(float v) {
1348         foreach(ref item; m)
1349             item -= v;
1350     }
1351     /// divide all components of vector by matrix
1352     void opOpAssign(string op : "/")(float v) {
1353         foreach(ref item; m)
1354             item /= v;
1355     }
1356 
1357     /// inplace rotate around Z axis
1358     ref mat4 rotatez(float angle) return {
1359         return rotate(angle, 0, 0, 1);
1360     }
1361 
1362     /// inplace rotate around X axis
1363     ref mat4 rotatex(float angle) return {
1364         return rotate(angle, 1, 0, 0);
1365     }
1366 
1367     /// inplace rotate around Y axis
1368     ref mat4 rotatey(float angle) return {
1369         return rotate(angle, 0, 1, 0);
1370     }
1371 
1372     ref mat4 rotate(float angle, const vec3 axis) return {
1373         return rotate(angle, axis.x, axis.y, axis.z);
1374     }
1375 
1376     ref mat4 rotate(float angle, float x, float y, float z) return {
1377         if (angle == 0.0f)
1378             return this;
1379         mat4 m;
1380         float c, s, ic;
1381         if (angle == 90.0f || angle == -270.0f) {
1382             s = 1.0f;
1383             c = 0.0f;
1384         } else if (angle == -90.0f || angle == 270.0f) {
1385             s = -1.0f;
1386             c = 0.0f;
1387         } else if (angle == 180.0f || angle == -180.0f) {
1388             s = 0.0f;
1389             c = -1.0f;
1390         } else {
1391             float a = angle * PI / 180.0f;
1392             c = cos(a);
1393             s = sin(a);
1394         }
1395         bool quick = false;
1396         if (x == 0.0f) {
1397             if (y == 0.0f) {
1398                 if (z != 0.0f) {
1399                     // Rotate around the Z axis.
1400                     m.setIdentity();
1401                     m.m[0*4 + 0] = c;
1402                     m.m[1*4 + 1] = c;
1403                     if (z < 0.0f) {
1404                         m.m[1*4 + 0] = s;
1405                         m.m[0*4 + 1] = -s;
1406                     } else {
1407                         m.m[1*4 + 0] = -s;
1408                         m.m[0*4 + 1] = s;
1409                     }
1410                     quick = true;
1411                 }
1412             } else if (z == 0.0f) {
1413                 // Rotate around the Y axis.
1414                 m.setIdentity();
1415                 m.m[0*4 + 0] = c;
1416                 m.m[2*4 + 2] = c;
1417                 if (y < 0.0f) {
1418                     m.m[2*4 + 0] = -s;
1419                     m.m[0*4 + 2] = s;
1420                 } else {
1421                     m.m[2*4 + 0] = s;
1422                     m.m[0*4 + 2] = -s;
1423                 }
1424                 quick = true;
1425             }
1426         } else if (y == 0.0f && z == 0.0f) {
1427             // Rotate around the X axis.
1428             m.setIdentity();
1429             m.m[1*4 + 1] = c;
1430             m.m[2*4 + 2] = c;
1431             if (x < 0.0f) {
1432                 m.m[2*4 + 1] = s;
1433                 m.m[1*4 + 2] = -s;
1434             } else {
1435                 m.m[2*4 + 1] = -s;
1436                 m.m[1*4 + 2] = s;
1437             }
1438             quick = true;
1439         }
1440         if (!quick) {
1441             float len = x * x + y * y + z * z;
1442             if (!fuzzyNull(len - 1.0f) && !fuzzyNull(len)) {
1443                 len = sqrt(len);
1444                 x /= len;
1445                 y /= len;
1446                 z /= len;
1447             }
1448             ic = 1.0f - c;
1449             m.m[0*4 + 0] = x * x * ic + c;
1450             m.m[1*4 + 0] = x * y * ic - z * s;
1451             m.m[2*4 + 0] = x * z * ic + y * s;
1452             m.m[3*4 + 0] = 0.0f;
1453             m.m[0*4 + 1] = y * x * ic + z * s;
1454             m.m[1*4 + 1] = y * y * ic + c;
1455             m.m[2*4 + 1] = y * z * ic - x * s;
1456             m.m[3*4 + 1] = 0.0f;
1457             m.m[0*4 + 2] = x * z * ic - y * s;
1458             m.m[1*4 + 2] = y * z * ic + x * s;
1459             m.m[2*4 + 2] = z * z * ic + c;
1460             m.m[3*4 + 2] = 0.0f;
1461             m.m[0*4 + 3] = 0.0f;
1462             m.m[1*4 + 3] = 0.0f;
1463             m.m[2*4 + 3] = 0.0f;
1464             m.m[3*4 + 3] = 1.0f;
1465         }
1466         this *= m;
1467         return this;
1468     }
1469 
1470     ref mat4 rotateX(float angle) return {
1471         return rotate(angle, 1, 0, 0);
1472     }
1473 
1474     ref mat4 rotateY(float angle) return {
1475         return rotate(angle, 0, 1, 0);
1476     }
1477 
1478     ref mat4 rotateZ(float angle) return {
1479         return rotate(angle, 0, 0, 1);
1480     }
1481 
1482     ref mat4 scale(float x, float y, float z) return {
1483         m[0*4 + 0] *= x;
1484         m[0*4 + 1] *= x;
1485         m[0*4 + 2] *= x;
1486         m[0*4 + 3] *= x;
1487         m[1*4 + 0] *= y;
1488         m[1*4 + 1] *= y;
1489         m[1*4 + 2] *= y;
1490         m[1*4 + 3] *= y;
1491         m[2*4 + 0] *= z;
1492         m[2*4 + 1] *= z;
1493         m[2*4 + 2] *= z;
1494         m[2*4 + 3] *= z;
1495         return this;
1496     }
1497 
1498     ref mat4 scale(float v) return {
1499         m[0*4 + 0] *= v;
1500         m[0*4 + 1] *= v;
1501         m[0*4 + 2] *= v;
1502         m[0*4 + 3] *= v;
1503         m[1*4 + 0] *= v;
1504         m[1*4 + 1] *= v;
1505         m[1*4 + 2] *= v;
1506         m[1*4 + 3] *= v;
1507         m[2*4 + 0] *= v;
1508         m[2*4 + 1] *= v;
1509         m[2*4 + 2] *= v;
1510         m[2*4 + 3] *= v;
1511         //m[3*4 + 0] *= v;
1512         //m[3*4 + 1] *= v;
1513         //m[3*4 + 2] *= v;
1514         //m[3*4 + 3] *= v;
1515         return this;
1516     }
1517 
1518     ref mat4 scale(const vec3 v) return {
1519         m[0*4 + 0] *= v.x;
1520         m[0*4 + 1] *= v.x;
1521         m[0*4 + 2] *= v.x;
1522         m[0*4 + 3] *= v.x;
1523         m[1*4 + 0] *= v.y;
1524         m[1*4 + 1] *= v.y;
1525         m[1*4 + 2] *= v.y;
1526         m[1*4 + 3] *= v.y;
1527         m[2*4 + 0] *= v.z;
1528         m[2*4 + 1] *= v.z;
1529         m[2*4 + 2] *= v.z;
1530         m[2*4 + 3] *= v.z;
1531         return this;
1532     }
1533 
1534     static mat4 translation(float x, float y, float z) {
1535         // TODO
1536         mat4 res = 1;
1537         return res;
1538     }
1539 
1540     /**
1541     * Decomposes the scale, rotation and translation components of this matrix.
1542     *
1543     * @param scale The scale.
1544     * @param rotation The rotation.
1545     * @param translation The translation.
1546     */
1547     bool decompose(vec3 * scale, vec4 * rotation, vec3 * translation) const {
1548         if (translation)
1549         {
1550             // Extract the translation.
1551             translation.x = m[12];
1552             translation.y = m[13];
1553             translation.z = m[14];
1554         }
1555 
1556         // Nothing left to do.
1557         if (!scale && !rotation)
1558             return true;
1559 
1560         // Extract the scale.
1561         // This is simply the length of each axis (row/column) in the matrix.
1562         vec3 xaxis = vec3(m[0], m[1], m[2]);
1563         float scaleX = xaxis.length();
1564 
1565         vec3 yaxis = vec3(m[4], m[5], m[6]);
1566         float scaleY = yaxis.length();
1567 
1568         vec3 zaxis = vec3(m[8], m[9], m[10]);
1569         float scaleZ = zaxis.length();
1570 
1571         // Determine if we have a negative scale (true if determinant is less than zero).
1572         // In this case, we simply negate a single axis of the scale.
1573         float det = determinant();
1574         if (det < 0)
1575             scaleZ = -scaleZ;
1576 
1577         if (scale)
1578         {
1579             scale.x = scaleX;
1580             scale.y = scaleY;
1581             scale.z = scaleZ;
1582         }
1583 
1584         // Nothing left to do.
1585         if (!rotation)
1586             return true;
1587 
1588         //// Scale too close to zero, can't decompose rotation.
1589         //if (scaleX < MATH_TOLERANCE || scaleY < MATH_TOLERANCE || fabs(scaleZ) < MATH_TOLERANCE)
1590         //    return false;
1591         // TODO: support rotation
1592         return false;
1593     }
1594 
1595     /**
1596     * Gets the translational component of this matrix in the specified vector.
1597     *
1598     * @param translation A vector to receive the translation.
1599     */
1600     void getTranslation(ref vec3 translation) const
1601     {
1602         decompose(null, null, &translation);
1603     }
1604 
1605     @property float determinant() const
1606     {
1607         float a0 = m[0] * m[5] - m[1] * m[4];
1608         float a1 = m[0] * m[6] - m[2] * m[4];
1609         float a2 = m[0] * m[7] - m[3] * m[4];
1610         float a3 = m[1] * m[6] - m[2] * m[5];
1611         float a4 = m[1] * m[7] - m[3] * m[5];
1612         float a5 = m[2] * m[7] - m[3] * m[6];
1613         float b0 = m[8] * m[13] - m[9] * m[12];
1614         float b1 = m[8] * m[14] - m[10] * m[12];
1615         float b2 = m[8] * m[15] - m[11] * m[12];
1616         float b3 = m[9] * m[14] - m[10] * m[13];
1617         float b4 = m[9] * m[15] - m[11] * m[13];
1618         float b5 = m[10] * m[15] - m[11] * m[14];
1619         // Calculate the determinant.
1620         return (a0 * b5 - a1 * b4 + a2 * b3 + a3 * b2 - a4 * b1 + a5 * b0);
1621     }
1622 
1623     @property vec3 forwardVector() const
1624     {
1625         return vec3(-m[8], -m[9], -m[10]);
1626     }
1627 
1628     @property vec3 backVector() const
1629     {
1630         return vec3(m[8], m[9], m[10]);
1631     }
1632 
1633     void transformVector(ref vec3 v) const {
1634         transformVector(v.x, v.y, v.z, 0, v);
1635     }
1636 
1637     void transformPoint(ref vec3 v) const {
1638         transformVector(v.x, v.y, v.z, 1, v);
1639     }
1640 
1641     void transformVector(float x, float y, float z, float w, ref vec3 dst) const {
1642         dst.x = x * m[0] + y * m[4] + z * m[8] + w * m[12];
1643         dst.y = x * m[1] + y * m[5] + z * m[9] + w * m[13];
1644         dst.z = x * m[2] + y * m[6] + z * m[10] + w * m[14];
1645     }
1646 
1647     static __gshared const mat4 IDENTITY;
1648 }
1649 
1650 
1651 
1652 unittest {
1653     vec3 a, b, c;
1654     a.clear(5);
1655     b.clear(2);
1656     float d = a * b;
1657     auto r1 = a + b;
1658     auto r2 = a - b;
1659     c = a; c += b;
1660     c = a; c -= b;
1661     c = a; c *= b;
1662     c = a; c /= b;
1663     c += 0.3f;
1664     c -= 0.3f;
1665     c *= 0.3f;
1666     c /= 0.3f;
1667     a.x += 0.5f;
1668     a.y += 0.5f;
1669     a.z += 0.5f;
1670     auto v = b.vec;
1671     a = [0.1f, 0.2f, 0.3f];
1672     a.normalize();
1673     c = b.normalized;
1674 }
1675 
1676 unittest {
1677     vec4 a, b, c;
1678     a.clear(5);
1679     b.clear(2);
1680     float d = a * b;
1681     auto r1 = a + b;
1682     auto r2 = a - b;
1683     c = a; c += b;
1684     c = a; c -= b;
1685     c = a; c *= b;
1686     c = a; c /= b;
1687     c += 0.3f;
1688     c -= 0.3f;
1689     c *= 0.3f;
1690     c /= 0.3f;
1691     a.x += 0.5f;
1692     a.y += 0.5f;
1693     a.z += 0.5f;
1694     auto v = b.vec;
1695     a = [0.1f, 0.2f, 0.3f, 0.4f];
1696     a.normalize();
1697     c = b.normalized;
1698 }
1699 
1700 unittest {
1701     mat4 m;
1702     m.setIdentity();
1703     m = [1.0f,2.0f,3.0f,4.0f,5.0f,6.0f,7.0f,8.0f,9.0f,10.0f,11.0f,12.0f,13.0f,14.0f,15.0f,16.0f];
1704     float r;
1705     r = m[1, 3];
1706     m[2, 1] = 0.0f;
1707     m += 1;
1708     m -= 2;
1709     m *= 3;
1710     m /= 3;
1711     m.translate(vec3(2, 3, 4));
1712     m.translate(5, 6, 7);
1713     m.lookAt(vec3(5, 5, 5), vec3(0, 0, 0), vec3(-1, 1, 1));
1714     m.setLookAt(vec3(5, 5, 5), vec3(0, 0, 0), vec3(-1, 1, 1));
1715     m.scale(2,3,4);
1716     m.scale(vec3(2,3,4));
1717 
1718     vec3 vv1 = vec3(1,2,3);
1719     auto p1 = m * vv1;
1720     vec3 vv2 = vec3(3,4,5);
1721     auto p2 = vv2 * m;
1722     auto p3 = vec4(1,2,3,4) * m;
1723     auto p4 = m * vec4(1,2,3,4);
1724 
1725     m.rotate(30, 1, 1, 1);
1726     m.rotateX(10);
1727     m.rotateY(10);
1728     m.rotateZ(10);
1729 }
1730 
1731 /// calculate normal for triangle
1732 vec3 triangleNormal(vec3 p1, vec3 p2, vec3 p3) {
1733     return vec3.crossProduct(p2 - p1, p3 - p2).normalized();
1734 }
1735 
1736 /// calculate normal for triangle
1737 vec3 triangleNormal(float[3] p1, float[3] p2, float[3] p3) {
1738     return vec3.crossProduct(vec3(p2) - vec3(p1), vec3(p3) - vec3(p2)).normalized();
1739 }
1740 
1741 /// Alias for 2d float point
1742 alias PointF = vec2;
1743 
1744 
1745 
1746 
1747 // this form can be used within shaders
1748 /// cubic bezier curve
1749 PointF bezierCubic(const PointF[] cp, float t) pure @nogc @safe
1750     in { assert(cp.length > 3); } do
1751 {
1752     // control points
1753     auto p0 = cp[0];
1754     auto p1 = cp[1];
1755     auto p2 = cp[2];
1756     auto p3 = cp[3];
1757 
1758     float u1 = (1.0 - t);
1759     float u2 = t * t;
1760     // the polynomials
1761     float b3 = u2 * t;
1762     float b2 = 3.0 * u2 * u1;
1763     float b1 = 3.0 * t * u1 * u1;
1764     float b0 = u1 * u1 * u1;
1765     // cubic bezier interpolation
1766     PointF p = p0 * b0 + p1 * b1 + p2 * b2 + p3 * b3;
1767     return p;
1768 }
1769 
1770 /// quadratic bezier curve (not tested)
1771 PointF bezierQuadratic(const PointF[] cp, float t) pure @nogc @safe
1772     in { assert(cp.length > 2); } do
1773 {
1774     auto p0 = cp[0];
1775     auto p1 = cp[1];
1776     auto p2 = cp[2];
1777 
1778     float u1 = (1.0 - t);
1779     float u2 = u1 * u1;
1780 
1781     float b2 = t * t;
1782     float b1 = 2.0 * u1 * t;
1783     float b0 = u2;
1784 
1785     PointF p = p0 * b0 + p1 * b1 + p2 * b2;
1786     return p;
1787 }
1788 
1789 /// cubic bezier (first) derivative
1790 PointF bezierCubicDerivative(const PointF[] cp, float t) pure @nogc @safe 
1791     in { assert(cp.length > 3); } do 
1792 {
1793     auto p0 = cp[0];
1794     auto p1 = cp[1];
1795     auto p2 = cp[2];
1796     auto p3 = cp[3];
1797 
1798     float u1 = (1.0 - t);
1799     float u2 = t * t;
1800     float u3 = 6*(u1)*t;
1801     float d0 = 3 * u1 * u1;
1802     // -3*P0*(1-t)^2 + P1*(3*(1-t)^2 - 6*(1-t)*t) + P2*(6*(1-t)*t - 3*t^2) + 3*P3*t^2
1803     PointF d = p0*(-d0) + p1*(d0 - u3) + p2*(u3 - 3*u2) + (p3*3)*u2;
1804     return d;
1805 }
1806 
1807 /// quadratic bezier (first) derivative
1808 PointF bezierQuadraticDerivative(const PointF[] cp, float t) pure @nogc @safe 
1809     in { assert(cp.length > 2); } do 
1810 {
1811     auto p0 = cp[0];
1812     auto p1 = cp[1];
1813     auto p2 = cp[2];
1814 
1815     float u1 = (1.0 - t);
1816     // -2*(1-t)*(p1-p0) + 2*t*(p2-p1);
1817     PointF d = (p0-p1)*-2*u1 + (p2-p1)*2*t;
1818     return d;
1819 }
1820 
1821 // can't be pure due to normalize & vec2 ctor
1822 /// evaluates cubic bezier direction(tangent) at point t
1823 PointF bezierCubicDirection(const PointF[] cp, float t) {
1824     auto d = bezierCubicDerivative(cp,t);
1825     d.normalize();
1826     return PointF(tan(d.x), tan(d.y));
1827 }
1828 
1829 /// evaluates quadratic bezier direction(tangent) at point t
1830 PointF bezierQuadraticDirection(const PointF[] cp, float t) {
1831     auto d = bezierQuadraticDerivative(cp,t);
1832     d.normalize();
1833     return PointF(tan(d.x), tan(d.y));
1834 }
1835 
1836 /// templated version of bezier flatten curve function, allocates temporary buffer
1837 PointF[] flattenBezier(alias BezierFunc)(const PointF[] cp, int segmentCountInclusive) 
1838     if (is(typeof(BezierFunc)==function)) 
1839 {
1840     if (segmentCountInclusive < 2)
1841         return PointF[].init;
1842     PointF[] coords = new PointF[segmentCountInclusive+1];
1843     flattenBezier!BezierFunc(cp, segmentCountInclusive, coords);
1844     return coords;
1845 }
1846 
1847 /// flatten bezier curve function, writes to provided buffer instead of allocation
1848 void flattenBezier(alias BezierFunc)(const PointF[] cp, int segmentCountInclusive, PointF[] outSegments)
1849     if (is(typeof(BezierFunc)==function)) 
1850 {
1851     if (segmentCountInclusive < 2)
1852         return;
1853     float step = 1f/segmentCountInclusive;
1854     outSegments[0] = BezierFunc(cp, 0);
1855     foreach(i; 1..segmentCountInclusive) {
1856         outSegments[i] = BezierFunc(cp, i*step);
1857     }
1858     outSegments[segmentCountInclusive] = BezierFunc(cp, 1f);
1859 }
1860 
1861 
1862 /// flattens cubic bezier curve, returns PointF[segmentCount+1] array or empty array if <1 segments
1863 PointF[] flattenBezierCubic(const PointF[] cp, int segmentCount) {
1864     return flattenBezier!bezierCubic(cp,segmentCount);
1865 }
1866 
1867 /// flattens quadratic bezier curve, returns PointF[segmentCount+1] array or empty array if <1 segments
1868 PointF[] flattenBezierQuadratic(const PointF[] cp, int segmentCount) {
1869     return flattenBezier!bezierQuadratic(cp, segmentCount);
1870 }
1871 
1872 /// calculates normal vector at point t using direction
1873 PointF bezierCubicNormal(const PointF[] cp, float t) {
1874     auto d = bezierCubicDirection(cp, t);
1875     return d.rotated90ccw;
1876 }
1877 
1878 /// calculates normal vector at point t using direction
1879 PointF bezierQuadraticNormal(const PointF[] cp, float t) {
1880     auto d = bezierQuadraticDerivative(cp, t);
1881     return d.rotated90ccw;
1882 }