1 /// Affine transforms module
2 module gfx.math.transform;
3 
4 import gfx.math.mat;
5 import gfx.math.vec;
6 import std.meta : allSatisfy;
7 import std.traits : CommonType, isFloatingPoint, isNumeric;
8 
9 pure @safe @nogc nothrow:
10 
11 /// Build a translation matrix.
12 auto translation(X, Y)(in X x, in Y y)
13 {
14     alias ResMat = Mat3x3!(CommonType!(X, Y));
15     return ResMat(
16         1, 0, x,
17         0, 1, y,
18         0, 0, 1,
19     );
20 }
21 
22 /// ditto
23 auto translation(V)(in V v) if (isVec2!V)
24 {
25     return Mat3x3!(V.Component) (
26         1, 0, v.x,
27         0, 1, v.y,
28         0, 0, 1,
29     );
30 }
31 
32 /// ditto
33 auto affineTranslation(X, Y)(in X x, in Y y)
34 {
35     alias ResMat = Mat2x3!(CommonType!(X, Y));
36 
37     return ResMat(
38         1, 0, x,
39         0, 1, y
40     );
41 }
42 
43 /// ditto
44 auto affineTranslation(V)(in V v) if (isVec2!V)
45 {
46     alias ResMat = Mat2x3!(V.Component);
47 
48     return ResMat(
49         1, 0, v.x,
50         0, 1, v.y
51     );
52 }
53 
54 /// ditto
55 auto translation(X, Y, Z)(in X x, in Y y, in Z z)
56 {
57     alias ResMat = Mat4x4!(CommonType!(X, Y, Z));
58     return ResMat(
59         1, 0, 0, x,
60         0, 1, 0, y,
61         0, 0, 1, z,
62         0, 0, 0, 1,
63     );
64 }
65 
66 /// ditto
67 auto translation(V)(in V v) if (isVec3!V)
68 {
69     return Mat4x4!(V.Component) (
70         1, 0, 0, v.x,
71         0, 1, 0, v.y,
72         0, 0, 1, v.z,
73         0, 0, 0, 1,
74     );
75 }
76 
77 unittest
78 {
79     import gfx.math.approx : approxUlp;
80 
81     immutable v2 = dvec(4, 6);
82     assert( approxUlp(translation(2, 7) * dvec(v2, 1), dvec(6, 13, 1)) );
83 
84     immutable v3 = dvec(5, 6, 7);
85     assert( approxUlp(translation(7, 4, 1) * dvec(v3, 1), dvec(12, 10, 8, 1)) );
86 }
87 
88 /// Append a translation transform inferred from arguments to the matrix m.
89 /// This is equivalent to the expression $(D_CODE translation(...) * m)
90 /// but actually save computation by knowing
91 /// where the ones and zeros are in a pure translation matrix.
92 M translate(M, X, Y)(in M m, in X x, in Y y)
93 if (isMat!(3, 3, M) && allSatisfy!(isNumeric, X, Y))
94 {
95     return M (
96         // row 1
97         m[0, 0] + m[2, 0] * x,
98         m[0, 1] + m[2, 1] * x,
99         m[0, 2] + m[2, 2] * x,
100         // row 2
101         m[1, 0] + m[2, 0] * y,
102         m[1, 1] + m[2, 1] * y,
103         m[1, 2] + m[2, 2] * y,
104         // row 3
105         m[2, 0], m[2, 1], m[2, 2]
106     );
107 }
108 
109 /// ditto
110 M translate(M, X, Y)(in M m, in X x, in Y y)
111 if (isMat!(2, 3, M) && allSatisfy!(isNumeric, X, Y))
112 {
113     return M (
114         // row 1
115         m[0, 0],
116         m[0, 1],
117         m[0, 2] + x,
118         // row 2
119         m[1, 0],
120         m[1, 1],
121         m[1, 2] + y,
122     );
123 }
124 
125 /// ditto
126 M translate(M, V)(in M m, in V v)
127 if ((isMat!(2, 3, M) || isMat!(3, 3, M)) && isVec!(2, V))
128 {
129     return translate (m, v[0] ,v[1]);
130 }
131 
132 /// ditto
133 M translate (M, X, Y, Z)(in M m, in X x, in Y y, in Z z)
134 if (isMat!(4, 4, M) && allSatisfy!(isNumeric, X, Y, Z))
135 {
136     return M (
137         // row 1
138         m[0, 0] + m[3, 0] * x,
139         m[0, 1] + m[3, 1] * x,
140         m[0, 2] + m[3, 2] * x,
141         m[0, 3] + m[3, 3] * x,
142         // row 2
143         m[1, 0] + m[3, 0] * y,
144         m[1, 1] + m[3, 1] * y,
145         m[1, 2] + m[3, 2] * y,
146         m[1, 3] + m[3, 3] * y,
147         // row 3
148         m[2, 0] + m[3, 0] * z,
149         m[2, 1] + m[3, 1] * z,
150         m[2, 2] + m[3, 2] * z,
151         m[2, 3] + m[3, 3] * z,
152         // row 4
153         m[3, 0], m[3, 1], m[3, 2], m[3, 3]
154     );
155 }
156 
157 /// ditto
158 M translate (M, X, Y, Z)(in M m, in X x, in Y y, in Z z)
159 if (isMat!(3, 4, M) && allSatisfy!(isNumeric, X, Y, Z))
160 {
161     return M (
162         // row 1
163         m[0, 0] + m[3, 0] * x,
164         m[0, 1] + m[3, 1] * x,
165         m[0, 2] + m[3, 2] * x,
166         m[0, 3] + m[3, 3] * x,
167         // row 2
168         m[1, 0] + m[3, 0] * y,
169         m[1, 1] + m[3, 1] * y,
170         m[1, 2] + m[3, 2] * y,
171         m[1, 3] + m[3, 3] * y,
172         // row 3
173         m[2, 0] + m[3, 0] * z,
174         m[2, 1] + m[3, 1] * z,
175         m[2, 2] + m[3, 2] * z,
176         m[2, 3] + m[3, 3] * z,
177     );
178 }
179 
180 /// ditto
181 M translate(M, V)(in M m, in V v)
182 if ((isMat!(3, 4, M) || isMat!(4, 4, M)) && isVec!(3, V))
183 {
184     return translate (m, v[0] ,v[1], v[2]);
185 }
186 
187 ///
188 unittest
189 {
190     import gfx.math.approx : approxUlp;
191 
192     immutable m = DMat3( 1, 2, 3, 4, 5, 6, 7, 8, 9 );
193 
194     immutable expected = translation(7, 12) * m;  // full multiplication
195     immutable result = translate(m, 7, 12);       // simplified multiplication
196 
197     assert (approxUlp(expected, result));
198 }
199 ///
200 unittest
201 {
202     import gfx.math.approx : approxUlp;
203 
204     immutable m = DMat4( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 );
205 
206     immutable expected = translation(7, 12, 18) * m;  // full multiplication
207     immutable result = translate(m, 7, 12, 18);       // simplified multiplication
208 
209     assert (approxUlp(expected, result));
210 }
211 
212 
213 /// Build a pure 3d rotation matrix with angle in radians
214 auto rotationPure(T, V) (in T angle, in V axis)
215 if (isFloatingPoint!T && isVec!(3, V))
216 {
217     static assert (
218         isFloatingPoint!(V.Component),
219         "rotationPure must be passed a floating point axis"
220     );
221     import std.math : cos, sin;
222     const u = normalize(axis);
223     const c = cos(angle);
224     const s = sin(angle);
225     const c1 = 1 - c;
226     return Mat3x3!(V.Component) (
227         // row 1
228         c1 * u.x * u.x  +  c,
229         c1 * u.x * u.y  -  s * u.z,
230         c1 * u.x * u.z  +  s * u.y,
231         // row 2
232         c1 * u.y * u.x  +  s * u.z,
233         c1 * u.y * u.y  +  c,
234         c1 * u.y * u.z  -  s * u.x,
235         // row 3
236         c1 * u.z * u.x  -  s * u.y,
237         c1 * u.z * u.y  +  s * u.x,
238         c1 * u.z * u.z  +  c,
239     );
240 }
241 
242 /// Build a rotation matrix.
243 /// angle in radians.
244 Mat3x3!T rotation(T) (in T angle) if (isFloatingPoint!T)
245 {
246     import std.math : cos, sin;
247     const c = cast(T) cos(angle);
248     const s = cast(T) sin(angle);
249     return Mat3x3!T (
250         c, -s, 0,
251         s, c, 0,
252         0, 0, 1
253     );
254 }
255 
256 Mat2x3!T affineRotation(T) (in T angle) if (isFloatingPoint!T)
257 {
258     import std.math : cos, sin;
259     const c = cast(T) cos(angle);
260     const s = cast(T) sin(angle);
261     return Mat2x3!T (
262         c, -s, 0,
263         s, c, 0,
264     );
265 }
266 
267 /// ditto
268 auto rotation(T, V) (in T angle, in V axis)
269 if (isVec!(3, V) && isFloatingPoint!T)
270 {
271     static assert (
272         isFloatingPoint!(V.Component),
273         "rotation must be passed a floating point axis"
274     );
275     const m = rotationPure(angle, axis);
276     return mat(
277         vec(m[0], 0),
278         vec(m[1], 0),
279         vec(m[2], 0),
280         vec(0, 0,  0,  1)
281     );
282 }
283 
284 /// ditto
285 auto rotation(T) (in T angle, in T x, in T y, in T z)
286 if (isFloatingPoint!T)
287 {
288     return rotation(angle, vec(x, y, z));
289 }
290 
291 /// Append a rotation transform inferred from arguments to the matrix m.
292 /// This is equivalent to the expression $(D_CODE rotation(...) * m)
293 /// but actually save computation by knowing
294 /// where the ones and zeros are in a pure rotation matrix.
295 M rotate (M, T) (in M m, in T angle)
296 if (isMat!(3, 3, M) && isFloatingPoint!T)
297 {
298     import std.math : cos, sin;
299     immutable c = cos(angle);
300     immutable s = sin(angle);
301     return M (
302         // row 1
303         c * m[0, 0] - s * m[1, 0],
304         c * m[0, 1] - s * m[1, 1],
305         c * m[0, 2] - s * m[1, 2],
306         // row 2
307         s * m[0, 0] + c * m[1, 0],
308         s * m[0, 1] + c * m[1, 1],
309         s * m[0, 2] + c * m[1, 2],
310         // row 3
311         m[2, 0], m[2, 1], m[2, 2]
312     );
313 }
314 
315 /// ditto
316 M rotate (M, T) (in M m, in T angle)
317 if (isMat!(2, 3, M) && isFloatingPoint!T)
318 {
319     import std.math : cos, sin;
320     immutable c = cos(angle);
321     immutable s = sin(angle);
322     return M (
323         // row 1
324         c * m[0, 0] - s * m[1, 0],
325         c * m[0, 1] - s * m[1, 1],
326         c * m[0, 2] - s * m[1, 2],
327         // row 2
328         s * m[0, 0] + c * m[1, 0],
329         s * m[0, 1] + c * m[1, 1],
330         s * m[0, 2] + c * m[1, 2]
331     );
332 }
333 
334 /// ditto
335 M rotate (M, T, V) (in M m, in T angle, in V axis)
336 if (isMat!(4, 4, M) && isFloatingPoint!T && isVec!(3, V))
337 {
338     static assert (
339         isFloatingPoint!(V.Component),
340         "rotate must be passed a floating point axis"
341     );
342     immutable r = rotationPure(angle, axis);
343     return M (
344         // row 1
345         r[0, 0]*m[0, 0] + r[0, 1]*m[1, 0] + r[0, 2]*m[2, 0],
346         r[0, 0]*m[0, 1] + r[0, 1]*m[1, 1] + r[0, 2]*m[2, 1],
347         r[0, 0]*m[0, 2] + r[0, 1]*m[1, 2] + r[0, 2]*m[2, 2],
348         r[0, 0]*m[0, 3] + r[0, 1]*m[1, 3] + r[0, 2]*m[2, 3],
349         // row 2
350         r[1, 0]*m[0, 0] + r[1, 1]*m[1, 0] + r[1, 2]*m[2, 0],
351         r[1, 0]*m[0, 1] + r[1, 1]*m[1, 1] + r[1, 2]*m[2, 1],
352         r[1, 0]*m[0, 2] + r[1, 1]*m[1, 2] + r[1, 2]*m[2, 2],
353         r[1, 0]*m[0, 3] + r[1, 1]*m[1, 3] + r[1, 2]*m[2, 3],
354         // row 3
355         r[2, 0]*m[0, 0] + r[2, 1]*m[1, 0] + r[2, 2]*m[2, 0],
356         r[2, 0]*m[0, 1] + r[2, 1]*m[1, 1] + r[2, 2]*m[2, 1],
357         r[2, 0]*m[0, 2] + r[2, 1]*m[1, 2] + r[2, 2]*m[2, 2],
358         r[2, 0]*m[0, 3] + r[2, 1]*m[1, 3] + r[2, 2]*m[2, 3],
359         // row 4
360         m[3, 0], m[3, 1], m[3, 2], m[3, 3]
361     );
362 }
363 
364 /// ditto
365 M rotate (M, T, V) (in M m, in T angle, in V axis)
366 if (isMat!(3, 4, M) && isVec!(3, V) && isFloatingPoint!T)
367 {
368     static assert (
369         isFloatingPoint!(V.Component),
370         "rotate must be passed a floating point axis"
371     );
372     const r = rotationPure(angle, axis);
373     return M (
374         // row 1
375         r[0, 0]*m[0, 0] + r[0, 1]*m[1, 0] + r[0, 2]*m[2, 0],
376         r[0, 0]*m[0, 1] + r[0, 1]*m[1, 1] + r[0, 2]*m[2, 1],
377         r[0, 0]*m[0, 2] + r[0, 1]*m[1, 2] + r[0, 2]*m[2, 2],
378         r[0, 0]*m[0, 3] + r[0, 1]*m[1, 3] + r[0, 2]*m[2, 3],
379         // row 2
380         r[1, 0]*m[0, 0] + r[1, 1]*m[1, 0] + r[1, 2]*m[2, 0],
381         r[1, 0]*m[0, 1] + r[1, 1]*m[1, 1] + r[1, 2]*m[2, 1],
382         r[1, 0]*m[0, 2] + r[1, 1]*m[1, 2] + r[1, 2]*m[2, 2],
383         r[1, 0]*m[0, 3] + r[1, 1]*m[1, 3] + r[1, 2]*m[2, 3],
384         // row 3
385         r[2, 0]*m[0, 0] + r[2, 1]*m[1, 0] + r[2, 2]*m[2, 0],
386         r[2, 0]*m[0, 1] + r[2, 1]*m[1, 1] + r[2, 2]*m[2, 1],
387         r[2, 0]*m[0, 2] + r[2, 1]*m[1, 2] + r[2, 2]*m[2, 2],
388         r[2, 0]*m[0, 3] + r[2, 1]*m[1, 3] + r[2, 2]*m[2, 3],
389     );
390 }
391 
392 /// ditto
393 M rotate (M, T) (in M m, in T angle, in T x, in T y, in T z)
394 if ((isMat!(3, 4, M) || isMat!(4, 4, M)) && isFloatingPoint!T)
395 {
396     return rotate(m, angle, vec(x, y, z));
397 }
398 
399 ///
400 unittest
401 {
402     import gfx.math.approx : approxUlp;
403     import std.math : PI;
404 
405     immutable m = DMat3( 1, 2, 3, 4, 5, 6, 7, 8, 9 );
406 
407     immutable expected = rotation!double(PI) * m; // full multiplication
408     immutable result = rotate(m, PI);      // simplified multiplication
409 
410     assert (approxUlp(expected, result));
411 }
412 ///
413 unittest
414 {
415     import gfx.math.approx : approxUlp;
416     import std.math : PI;
417 
418     immutable m = DMat4( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 );
419     immutable angle = PI;
420     immutable v = fvec(3, 4, 5);
421 
422     immutable expected = rotation(angle, v) * m; // full multiplication
423     immutable result = rotate(m, angle, v);      // simplified multiplication
424 
425     assert (approxUlp(expected, result));
426 }
427 
428 /// Build a rotation matrix from Euler angles
429 /// The convention taken is Xa, Zb, Xc
430 Mat3!T eulerAngles(T) (in T a, in T b, in T c)
431 if (isFloatingPoint!T)
432 {
433     import std.math : cos, sin;
434 
435     immutable sa = sin(a);
436     immutable sb = sin(b);
437     immutable sc = sin(c);
438 
439     immutable ca = cos(a);
440     immutable cb = cos(b);
441     immutable cc = cos(c);
442 
443     return Mat3!T (
444         cb,     -cc*sb,             sb*sc,
445         ca*sb,  ca*cb*cc - sa*sc,   -cc*sa - ca*cb*sc,
446         sa*sb,  ca*sc + cb*cc*sa,   ca*cc - cb*sa*sc,
447     );
448 }
449 
450 auto eulerAngles(V) (in V angles)
451 if (isVec!(3, V) && isFloatingPoint!(V.Component))
452 {
453     return eulerAngles(angles[0], angles[1], angles[2]);
454 }
455 
456 /// ditto
457 
458 ///
459 unittest {
460     import gfx.math.approx : approxUlpAndAbs;
461     import std.math : PI;
462 
463     const v = fvec(0, 0, 1);
464     // rotate PI/2 around X, then no rotation around Z, then again PI/2 around X
465     const m = eulerAngles(float(PI / 2), 0f, float(PI / 2));
466     const result = m * v;
467     const expected = fvec(0, 0, -1);
468 
469     assert(approxUlpAndAbs(result, expected));
470 }
471 
472 
473 /// Build a scale matrix.
474 Mat3!(CommonType!(X, Y)) scale(X, Y) (in X x, in Y y)
475 if (allSatisfy!(isNumeric, X, Y))
476 {
477     return Mat3!(CommonType!(X, Y))(
478         x, 0, 0,
479         0, y, 0,
480         0, 0, 1,
481     );
482 }
483 
484 /// ditto
485 auto scale(V) (in V v) if (isVec2!V)
486 {
487     return Mat3!(V.Component) (
488         v.x, 0, 0,
489         0, v.y, 0,
490         0, 0, 1,
491     );
492 }
493 
494 /// ditto
495 Mat2x3!(CommonType!(X, Y)) affineScale(X, Y) (in X x, in Y y)
496 if (allSatisfy!(isNumeric, X, Y))
497 {
498     return Mat2x3!(CommonType!(X, Y))(
499         x, 0, 0,
500         0, y, 0,
501     );
502 }
503 
504 /// ditto
505 auto affineScale(V) (in V v) if (isVec2!V)
506 {
507     return Mat2x3!(V.Component) (
508         v.x, 0, 0,
509         0, v.y, 0,
510     );
511 }
512 
513 /// ditto
514 Mat4!(CommonType!(X, Y, Z)) scale (X, Y, Z) (in X x, in Y y, in Z z)
515 if (allSatisfy!(isNumeric, X, Y, Z))
516 {
517     return Mat4!(CommonType!(X, Y, Z))(
518         x, 0, 0, 0,
519         0, y, 0, 0,
520         0, 0, z, 0,
521         0, 0, 0, 1,
522     );
523 }
524 
525 /// ditto
526 auto scale(V) (in V v) if (isVec3!V)
527 {
528     return Mat4!(V.Component) (
529         v.x, 0, 0, 0,
530         0, v.y, 0, 0,
531         0, 0, v.z, 0,
532         0, 0, 0, 1,
533     );
534 }
535 
536 /// Append a scale transform inferred from arguments to the matrix m.
537 /// This is equivalent to the expression $(D_CODE scale(...) * m)
538 /// but actually save computation by knowing
539 /// where the ones and zeros are in a pure scale matrix.
540 M scale (M, X, Y)(in M m, in X x, in Y y)
541 if (isMat!(3, 3, M) && allSatisfy!(isNumeric, X, Y))
542 {
543     return M (
544         // row 1
545         m[0, 0] * x,
546         m[0, 1] * x,
547         m[0, 2] * x,
548         // row 2
549         m[1, 0] * y,
550         m[1, 1] * y,
551         m[1, 2] * y,
552         // row 3
553         m[2, 0], m[2, 1], m[2, 2]
554     );
555 }
556 
557 /// ditto
558 M scale (M, X, Y)(in M m, in X x, in Y y)
559 if (isMat!(2, 3, M) && allSatisfy!(isNumeric, X, Y))
560 {
561     return M (
562         // row 1
563         m[0, 0] * x,
564         m[0, 1] * x,
565         m[0, 2] * x,
566         // row 2
567         m[1, 0] * y,
568         m[1, 1] * y,
569         m[1, 2] * y,
570     );
571 }
572 
573 /// ditto
574 M scale (M, V)(in M m, in V v)
575 if ((isMat!(2, 3, M) || isMat!(3, 3, M)) && isVec!(2, V))
576 {
577     return scale(m, v[0], v[1]);
578 }
579 
580 /// ditto
581 M scale (M, X, Y, Z)(in M m, in X x, in Y y, in Z z)
582 if (isMat!(4, 4, M) && allSatisfy!(isNumeric, X, Y, Z))
583 {
584     return M (
585         // row 1
586         m[0, 0] * x,
587         m[0, 1] * x,
588         m[0, 2] * x,
589         m[0, 3] * x,
590         // row 2
591         m[1, 0] * y,
592         m[1, 1] * y,
593         m[1, 2] * y,
594         m[1, 3] * y,
595         // row 3
596         m[2, 0] * z,
597         m[2, 1] * z,
598         m[2, 2] * z,
599         m[2, 3] * z,
600         // row 4
601         m[3, 0], m[3, 1], m[3, 2], m[3, 3]
602     );
603 }
604 
605 /// ditto
606 M scale (M, X, Y, Z)(in M m, in X x, in Y y, in Z z)
607 if (isMat!(3, 4, M) && allSatisfy!(isNumeric, X, Y, Z))
608 {
609     return M (
610         // row 1
611         m[0, 0] * x,
612         m[0, 1] * x,
613         m[0, 2] * x,
614         m[0, 3] * x,
615         // row 2
616         m[1, 0] * y,
617         m[1, 1] * y,
618         m[1, 2] * y,
619         m[1, 3] * y,
620         // row 3
621         m[2, 0] * z,
622         m[2, 1] * z,
623         m[2, 2] * z,
624         m[2, 3] * z,
625     );
626 }
627 
628 /// ditto
629 M scale (M, V)(in M m, in V v)
630 if ((isMat!(3, 4, M) || isMat!(4, 4, M)) && isVec!(3, V))
631 {
632     return scale(m, v[0], v[1], v[2]);
633 }
634 
635 
636 ///
637 unittest
638 {
639     import gfx.math.approx : approxUlp;
640 
641     immutable m = DMat3( 1, 2, 3, 4, 5, 6, 7, 8, 9 );
642 
643     immutable expected = scale(4, 5) * m; // full multiplication
644     immutable result = scale(m, 4, 5);   // simplified multiplication
645 
646     assert (approxUlp(expected, result));
647 }
648 ///
649 unittest
650 {
651     import gfx.math.approx : approxUlp;
652 
653     immutable m = DMat4( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 );
654 
655     immutable expected = scale(4, 5, 6) * m; // full multiplication
656     immutable result = scale(m, 4, 5, 6);   // simplified multiplication
657 
658     assert (approxUlp(expected, result));
659 }
660 
661 /// Affine matrix multiplication.
662 ///
663 /// Perform a multiplication of two 2x3 or two 3x4 matrices as if their last row
664 /// were [ 0, 0, 1] or [ 0, 0, 0, 1 ].
665 /// Allows manipulation of smaller matrices when only affine transformation are
666 /// required.
667 /// Note: translation, rotation, scaling, shearing and any combination of those
668 /// are affine transforms. Projection is not affine.
669 /// I.e. for 2D, an affine transform is held in 2x3 matrix, 2x2 for rotation and
670 /// scaling and an additional column for translation.
671 /// The same applies with 3D and 3x4 matrices.
672 ///
673 /// This is not implemented as an operator as it is not a mathematical
674 /// operation (ml.columnLength != mr.rowLength).
675 auto affineMult(ML, MR)(in ML ml, in MR mr)
676 if (areMat!(2, 3, ML, MR) || areMat!(3, 4, ML, MR))
677 {
678     alias Comp = CommonType!(ML.Component, MR.Component);
679     enum rowLength = ML.rowLength;
680     enum columnLength = ML.columnLength;
681     alias ResMat = Mat!(Comp, rowLength, columnLength);
682 
683     ResMat res = void;
684     static foreach(r; 0 .. rowLength)
685     {
686         static foreach (c; 0 .. columnLength)
687         {{
688             Comp resComp = 0;
689             static foreach (rc; 0 .. rowLength) // that is columnCount-1
690             {
691                 resComp += ml[r, rc] * mr[rc, c];
692             }
693             static if (c == columnLength-1)
694             {
695                 resComp += ml[r, c]; // that is the last one in the last row
696             }
697             res[r, c] = resComp;
698         }}
699     }
700     return res;
701 }
702 
703 ///
704 unittest
705 {
706     import gfx.math.approx : approxUlp;
707 
708     /// full matrices
709     immutable fm1 = FMat3x3(
710         1, 2, 3,
711         4, 5, 6,
712         0, 0, 1
713     );
714     immutable fm2 = DMat3x3(
715          7,  8,  9,
716         10, 11, 12,
717          0,  0,  1
718     );
719 
720     /// affine matrices
721     immutable am1 = FMat2x3(
722         1, 2, 3,
723         4, 5, 6,
724     );
725     immutable am2 = DMat2x3(
726          7,  8,  9,
727         10, 11, 12,
728     );
729 
730     immutable expected = (fm1 * fm2).slice!(0, 2, 0, 3);
731     immutable result = affineMult(am1, am2);
732     assert( approxUlp(expected, result) );
733 }
734 
735 
736 /// Transform a vector by a matrix in homogenous coordinates.
737 auto transform(V, M)(in V v, in M m)
738 if (isVec!(2, V) && isMat!(3, 3, M))
739 {
740     return (m * vec(v, 1)).xy;
741 }
742 /// ditto
743 auto transform(V, M)(in V v, in M m)
744 if (isVec!(2, V) && isMat!(2, 3, M))
745 {
746     return m * vec(v, 1);
747 }
748 /// ditto
749 auto transform(V, M)(in V v, in M m)
750 if (isVec!(3, V) && isMat!(4, 4, M))
751 {
752     return (m * vec(v, 1)).xyz;
753 }
754 /// ditto
755 auto transform(V, M)(in V v, in M m)
756 if (isVec!(3, V) && isMat!(3, 4, M))
757 {
758     return m * vec(v, 1);
759 }
760 /// ditto
761 auto transform(V, M)(in V v, in M m)
762 if (isVec!(4, V) && isMat!(4, 4, M))
763 {
764     return m * v;
765 }
766 
767 unittest
768 {
769     import gfx.math.approx : approxUlp;
770 
771     // 2x3 matrix can hold affine 2D transforms
772     immutable transl = DMat2x3(
773         1, 0, 3,
774         0, 1, 2,
775     );
776     assert( approxUlp(transform(dvec(3, 5), transl), dvec(6, 7)) );
777 }
778 
779 ///
780 unittest
781 {
782     import gfx.math.approx : approxUlp, approxUlpAndAbs;
783     import std.math : PI;
784 
785     immutable v = dvec(2, 0);
786     auto m = DMat2x3.identity;
787 
788     m = m.rotate(PI/2);
789     assert ( approxUlpAndAbs(transform(v, m), dvec(0, 2)) );
790 
791     m = m.translate(2, 2);
792     assert ( approxUlp(transform(v, m), dvec(2, 4)) );
793 
794     m = m.scale(2, 2);
795     assert ( approxUlp(transform(v, m), dvec(4, 8)) );
796 }
797 
798 ///
799 unittest
800 {
801     import gfx.math.approx : approxUlp;
802 
803     auto st = scale!float(2, 2).translate(3, 1);
804     assert( approxUlp(transform(fvec(0, 0), st), fvec(3, 1)) );
805     assert( approxUlp(transform(fvec(1, 1), st), fvec(5, 3)) );
806 
807     auto ts = translation!float(3, 1).scale(2, 2);
808     assert( approxUlp(transform(fvec(0, 0), ts), fvec(6, 2)) );
809     assert( approxUlp(transform(fvec(1, 1), ts), fvec(8, 4)) );
810 }