1 /** Geometry types and algorithms.
2 
3    Special thanks to:
4    $(UL
5    $(LI Tomasz Stachowiak (h3r3tic): allowed me to use parts of $(LINK2 https://bitbucket.org/h3r3tic/boxen/src/default/src/xf/omg, omg).)
6    $(LI Jakob Øvrum (jA_cOp): improved the code a lot!)
7    $(LI Florian Boesch (___doc__): helps me to understand opengl/complex maths better, see: $(LINK http://codeflow.org/).)
8    $(LI #D on freenode: answered general questions about D.)
9    )
10 
11    Note: All methods marked with pure are weakly pure since, they all access an
12    instance member.  All static methods are strongly pure.
13 
14    TODO Support radian and degree types (from units-d package)
15 
16    TODO Use `sink` as param in `toMathML` and `toLaTeX`
17 
18    TODO Replace toMathML() with fmt argument %M to toString functions
19    TODO Replace toLaTeX() with fmt argument %L to toString functions
20 
21    TODO Optimize using core.simd or std.simd
22    TODO Merge with analyticgeometry
23    TODO Merge with https://github.com/CyberShadow/ae/blob/master/utils/geometry.d
24    TODO Integrate with http://code.dlang.org/packages/blazed2
25    TODO logln, log.warn, log.error, log.info, log.debug
26    TODO Make use of staticReduce etc when they become available in Phobos.
27    TODO Go through all usages of real and use CommonType!(real, E) to make it work when E is a bignum.
28    TODO ead and perhaps make use of http://stackoverflow.com/questions/3098242/fast-vector-struct-that-allows-i-and-xyz-operations-in-d?rq=1
29    TODO Tag member functions in t_geom.d as pure as is done https://github.com/D-Programming-Language/phobos/blob/master/std/bigint.d
30    TODO Remove need to use [] in x[] == y[]
31 
32    See: https://www.google.se/search?q=point+plus+vector
33    See: http://mosra.cz/blog/article.php?a=22-introducing-magnum-a-multiplatform-2d-3d-graphics-engine
34 */
35 module old_geometry;
36 
37 // TODO use import core.simd;
38 import std.math: sqrt, PI, sin, cos, acos;
39 import std.traits: isFloatingPoint, isNumeric, isSigned, isDynamicArray, isAssignable, isArray, CommonType, isInstanceOf;
40 import std..string: format, rightJustify;
41 import std.array: join;
42 import std.algorithm.iteration : map, reduce;
43 import std.algorithm.searching : all, any;
44 import std.algorithm.comparison : min, max;
45 import std.random: uniform;
46 
47 import nxt.mathml;
48 
49 enum isVector(E)     = is(typeof(isVectorImpl(E.init)));
50 enum isPoint(E)      = is(typeof(isPointImpl(E.init)));
51 enum isMatrix(E)     = is(typeof(isMatrixImpl(E.init)));
52 enum isQuaternion(E) = is(typeof(isQuaternionImpl(E.init)));
53 enum isPlane(E)      = is(typeof(isPlaneImpl(E.init)));
54 
55 private void isVectorImpl    (E, uint D)        (Vector    !(E, D)    vec) {}
56 private void isPointImpl     (E, uint D)        (Point     !(E, D)    vec) {}
57 private void isMatrixImpl    (E, uint R, uint C)(Matrix    !(E, R, C) mat) {}
58 private void isQuaternionImpl(E)                (Quaternion!(E)        qu) {}
59 private void isPlaneImpl     (E)                (PlaneT    !(E)         p) {}
60 
61 enum isFixVector(E) = isFix(typeof(isFixVectorImpl(E.init)));
62 enum isFixPoint(E)  = isFix(typeof(isFixPointImpl (E.init)));
63 enum isFixMatrix(E) = isFix(typeof(isFixMatrixImpl(E.init)));
64 
65 private void isFixVectorImpl (E, uint D)        (Vector!(E, D)    vec) {}
66 private void isFixPointImpl  (E, uint D)        (Point !(E, D)    vec) {}
67 private void isFixMatrixImpl (E, uint R, uint C)(Matrix!(E, R, C) mat) {}
68 
69 // See_Also: http://stackoverflow.com/questions/18552454/using-ctfe-to-generate-set-of-struct-aliases/18553026?noredirect=1#18553026
70 version(none)
71 string makeInstanceAliases(in string templateName,
72                            string aliasName = "",
73                            in uint minDimension = 2,
74                            in uint maxDimension = 4,
75                            in string[] elementTypes = defaultElementTypes)
76 in
77 {
78     assert(templateName.length);
79     assert(minDimension <= maxDimension);
80 }
81 do
82 {
83     import std..string : toLower;
84     import std.conv : to;
85     string code;
86     if (!aliasName.length)
87     {
88         aliasName = templateName.toLower;
89     }
90     foreach (immutable n; minDimension .. maxDimension + 1)
91     {
92         foreach (const et; elementTypes) // for each elementtype
93         {
94             immutable prefix = ("alias " ~ templateName ~ "!("~et~", " ~
95                                 to!string(n) ~ ") " ~ aliasName ~ "" ~
96                                 to!string(n));
97             if (et == "float")
98             {
99                 code ~= (prefix ~ ";\n"); // GLSL-style prefix-less single precision
100             }
101             code ~= (prefix ~ et[0] ~ ";\n");
102         }
103     }
104     return code;
105 }
106 
107 version(none)
108 mixin(makeInstanceAliases("Point", "point", 2,3,
109                            ["int", "float", "double", "real"]));
110 
111 /* Copied from https://github.com/CyberShadow/ae/blob/master/utils/geometry.d */
112 auto sqrtx(T)(T x)
113 {
114     static if (is(T : int))
115     {
116         return std.math.sqrt(cast(float)x);
117     }
118     else
119     {
120         return std.math.sqrt(x);
121     }
122 }
123 
124 import std.meta : AliasSeq;
125 
126 version(unittest)
127 {
128     static foreach (T; AliasSeq!(ubyte, int, float, double, real))
129     {
130         static foreach (uint n; 2 .. 4 + 1)
131         {
132         }
133     }
134 
135     alias vec2b = Vector!(byte, 2, false);
136 
137     alias vec2f = Vector!(float, 2, true);
138     alias vec3f = Vector!(float, 3, true);
139 
140     alias vec2d = Vector!(float, 2, true);
141 
142     alias nvec2f = Vector!(float, 2, true);
143 }
144 
145 // mixin(makeInstanceAliases("Vector", "vec", 2,4,
146 //                           ["ubyte", "int", "float", "double", "real"]));
147 
148 ///
149 @safe pure nothrow @nogc unittest
150 {
151     assert(vec2f(2, 3)[] == [2, 3].s);
152     assert(vec2f(2, 3)[0] == 2);
153     assert(vec2f(2) == 2);
154     assert(vec2f(true) == true);
155     assert(vec2b(true) == true);
156     assert(all!"a"(vec2b(true)[]));
157     assert(any!"a"(vec2b(false, true)[]));
158     assert(any!"a"(vec2b(true, false)[]));
159     assert(!any!"a"(vec2b(false, false)[]));
160     assert((vec2f(1, 3)*2.5f)[] == [2.5f, 7.5f].s);
161     nvec2f v = vec2f(3, 4);
162     assert(v[] == nvec2f(0.6, 0.8)[]);
163 }
164 
165 ///
166 @safe unittest
167 {
168     import std.conv : to;
169     auto x = vec2f(2, 3);
170     assert(to!string(vec2f(2, 3)) == `ColumnVector(2,3)`);
171     assert(to!string(transpose(vec2f(11, 22))) == `RowVector(11,22)`);
172     assert(vec2f(11, 22).toLaTeX == `\begin{pmatrix} 11 \\ 22 \end{pmatrix}`);
173     assert(vec2f(11, 22).T.toLaTeX == `\begin{pmatrix} 11 & 22 \end{pmatrix}`);
174 }
175 
176 auto transpose(E, uint D,
177                bool normalizedFlag,
178                Orient orient)(in Vector!(E, D,
179                                          normalizedFlag,
180                                          orient) a)
181 {
182     static if (orient == Orient.row)
183     {
184         return Vector!(E, D, normalizedFlag, Orient.column)(a._vector);
185     }
186     else
187     {
188         return Vector!(E, D, normalizedFlag, Orient.row)(a._vector);
189     }
190 }
191 alias T = transpose; // C++ Armadillo naming convention.
192 
193 auto elementwiseLessThanOrEqual(Ta, Tb,
194                                 uint D,
195                                 bool normalizedFlag,
196                                 Orient orient)(in Vector!(Ta, D, normalizedFlag, orient) a,
197                                                in Vector!(Tb, D, normalizedFlag, orient) b)
198 {
199     Vector!(bool, D) c = void;
200     static foreach (i; 0 .. D)
201     {
202         c[i] = a[i] <= b[i];
203     }
204     return c;
205 }
206 
207 @safe pure nothrow @nogc unittest
208 {
209     assert(elementwiseLessThanOrEqual(vec2f(1, 1),
210                                       vec2f(2, 2)) == vec2b(true, true));
211 }
212 
213 /// Returns: Scalar/Dot-Product of Two Vectors `a` and `b`.
214 T dotProduct(T, U)(in T a, in U b)
215 if (isInstanceOf!(Vector, T) &&
216     isInstanceOf!(Vector, U) &&
217     (T.dimension ==
218      U.dimension))
219 {
220     T c = void;
221     static foreach (i; 0 .. T.dimension)
222     {
223         c[i] = a[i] * b[i];
224     }
225     return c;
226 }
227 alias dot = dotProduct;
228 
229 /// Returns: Outer-Product of Two Vectors `a` and `b`.
230 auto outerProduct(Ta, Tb, uint Da, uint Db)(in Vector!(Ta, Da) a,
231                                             in Vector!(Tb, Db) b)
232 if (Da >= 1 &&
233     Db >= 1)
234 {
235     Matrix!(CommonType!(Ta, Tb), Da, Db) y = void;
236     static foreach (r; 0 .. Da)
237     {
238         static foreach (c; 0 .. Db)
239         {
240             y.at(r,c) = a[r] * b[c];
241         }
242     }
243     return y;
244 }
245 alias outer = outerProduct;
246 
247 /// Returns: Vector/Cross-Product of two 3-Dimensional Vectors.
248 auto crossProduct(T, U)(in T a,
249                         in U b)
250 if (isInstanceOf!(Vector, T) && T.dimension == 3 &&
251     isInstanceOf!(Vector, U) && U.dimension == 3)
252 {
253     return T(a.y * b.z - b.y * a.z,
254              a.z * b.x - b.z * a.x,
255              a.x * b.y - b.x * a.y);
256 }
257 
258 /// Returns: (Euclidean) Distance between `a` and `b`.
259 real distance(T, U)(in T a,
260                     in U b)
261 if ((isInstanceOf!(Vector, T) && // either both vectors
262      isInstanceOf!(Vector, U) &&
263      T.dimension == U.dimension) ||
264     (isPoint!T && // or both points
265      isPoint!U))  // TODO support distance between vector and point
266 {
267     return (a - b).magnitude;
268 }
269 
270 @safe pure nothrow @nogc unittest
271 {
272     auto v1 = vec3f(1, 2, -3);
273     auto v2 = vec3f(1, 3, 2);
274     assert(crossProduct(v1, v2)[] == [13, -5, 1].s);
275     assert(distance(vec2f(0, 0), vec2f(0, 10)) == 10);
276     assert(distance(vec2f(0, 0), vec2d(0, 10)) == 10);
277     assert(dot(v1, v2) == dot(v2, v1)); // commutative
278 }
279 
280 enum Layout { columnMajor, rowMajor }; // Matrix Storage Major Dimension.
281 
282 /// Base template for all matrix-types.
283 /// Params:
284 ///  E = all values get stored as this type
285 ///  rows_ = rows of the matrix
286 ///  cols_ = columns of the matrix
287 ///  layout = matrix layout
288 struct Matrix(E, uint rows_, uint cols_,
289               Layout layout = Layout.rowMajor)
290 if (rows_ >= 1 &&
291     cols_ >= 1)
292 {
293     alias mT = E; /// Internal type of the _matrix
294     static const uint rows = rows_; /// Number of rows
295     static const uint cols = cols_; /// Number of columns
296 
297     /** Matrix $(RED row-major) in memory. */
298     static if (layout == Layout.rowMajor)
299     {
300         E[cols][rows] _matrix; // In C it would be mt[rows][cols], D does it like this: (mt[cols])[rows]
301 
302         ref inout(E) opCall()(uint row, uint col) inout
303         {
304             return _matrix[row][col];
305         }
306 
307         ref inout(E) at()(uint row, uint col) inout
308         {
309             return _matrix[row][col];
310         }
311     }
312     else
313     {
314         E[rows][cols] _matrix; // In C it would be mt[cols][rows], D does it like this: (mt[rows])[cols]
315         ref inout(E) opCall()(uint row, uint col) inout
316         {
317             return _matrix[col][row];
318         }
319 
320         ref inout(E) at()(uint row, uint col) inout
321         {
322             return _matrix[col][row];
323         }
324     }
325     alias _matrix this;
326 
327     /// Returns: The pointer to the stored values as OpenGL requires it.
328     /// Note this will return a pointer to a $(RED row-major) _matrix,
329     /// $(RED this means you've to set the transpose argument to GL_TRUE when passing it to OpenGL).
330     /// Examples:
331     /// ---
332     /// // 3rd argument = GL_TRUE
333     /// glUniformMatrix4fv(programs.main.model, 1, GL_TRUE, mat4.translation(-0.5f, -0.5f, 1.0f).value_ptr);
334     /// ---
335     // @property auto value_ptr() { return _matrix[0].ptr; }
336 
337     /// Returns: The current _matrix formatted as flat string.
338 
339     @property void toString(scope void delegate(scope const(char)[]) @safe sink) const
340     {
341         import std.conv : to;
342         sink(`Matrix(`);
343         sink(to!string(_matrix));
344         sink(`)`);
345     }
346 
347     @property string toLaTeX()() const
348     {
349         string s;
350         static foreach (r; 0 .. rows)
351         {
352             static foreach (c; 0 .. cols)
353             {
354                 s ~= to!string(at(r, c)) ;
355                 if (c != cols - 1) { s ~= ` & `; } // if not last column
356             }
357             if (r != rows - 1) { s ~= ` \\ `; } // if not last row
358         }
359         return `\begin{pmatrix} ` ~ s ~ ` \end{pmatrix}`;
360     }
361 
362     @property string toMathML()() const
363     {
364         // opening
365         string str = `<math><mrow>
366   <mo>❲</mo>
367   <mtable>`;
368 
369         static foreach (r; 0 .. rows)
370         {
371             str ~=  `
372     <mtr>`;
373             static foreach (c; 0 .. cols)
374             {
375                 str ~= `
376       <mtd>
377         <mn>` ~ at(r, c).toMathML ~ `</mn>
378       </mtd>`;
379             }
380             str ~=  `
381     </mtr>`;
382         }
383 
384         // closing
385         str ~= `
386   </mtable>
387   <mo>❳</mo>
388 </mrow></math>
389 `;
390         return str;
391     }
392 
393     /// Returns: The current _matrix as pretty formatted string.
394     @property string toPrettyString()()
395     {
396         string fmtr = "%s";
397 
398         size_t rjust = max(format(fmtr, reduce!(max)(_matrix[])).length,
399                            format(fmtr, reduce!(min)(_matrix[])).length) - 1;
400 
401         string[] outer_parts;
402         foreach (E[] row; _matrix)
403         {
404             string[] inner_parts;
405             foreach (E col; row)
406             {
407                 inner_parts ~= rightJustify(format(fmtr, col), rjust);
408             }
409             outer_parts ~= " [" ~ join(inner_parts, ", ") ~ "]";
410         }
411 
412         return "[" ~ join(outer_parts, "\n")[1..$] ~ "]";
413     }
414 
415     static void isCompatibleMatrixImpl(uint r, uint c)(Matrix!(E, r, c) m) {}
416 
417     enum isCompatibleMatrix(T) = is(typeof(isCompatibleMatrixImpl(T.init)));
418 
419     static void isCompatibleVectorImpl(uint d)(Vector!(E, d) vec) {}
420 
421     enum isCompatibleVector(T) = is(typeof(isCompatibleVectorImpl(T.init)));
422 
423     private void construct(uint i, T, Tail...)(T head, Tail tail)
424     {
425         static if (i >= rows*cols)
426         {
427             static assert(0, "Too many arguments passed to constructor");
428         }
429         else static if (is(T : E))
430         {
431             _matrix[i / cols][i % cols] = head;
432             construct!(i + 1)(tail);
433         }
434         else static if (is(T == Vector!(E, cols)))
435         {
436             static if (i % cols == 0)
437             {
438                 _matrix[i / cols] = head._vector;
439                 construct!(i + T.dimension)(tail);
440             }
441             else
442             {
443                 static assert(0, "Can't convert Vector into the matrix. Maybe it doesn't align to the columns correctly or dimension doesn't fit");
444             }
445         }
446         else
447         {
448             static assert(0, "Matrix constructor argument must be of type " ~ E.stringof ~ " or Vector, not " ~ T.stringof);
449         }
450     }
451 
452     private void construct(uint i)()  // terminate
453     {
454         static assert(i == rows*cols, "Not enough arguments passed to constructor");
455     }
456 
457     /// Constructs the matrix:
458     /// If a single value is passed, the matrix will be cleared with this value (each column in each row will contain this value).
459     /// If a matrix with more rows and columns is passed, the matrix will be the upper left nxm matrix.
460     /// If a matrix with less rows and columns is passed, the passed matrix will be stored in the upper left of an identity matrix.
461     /// It's also allowed to pass vectors and scalars at a time, but the vectors dimension must match the number of columns and align correctly.
462     /// Examples:
463     /// ---
464     /// mat2 m2 = mat2(0.0f); // mat2 m2 = mat2(0.0f, 0.0f, 0.0f, 0.0f);
465     /// mat3 m3 = mat3(m2); // mat3 m3 = mat3(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f);
466     /// mat3 m3_2 = mat3(vec3(1.0f, 2.0f, 3.0f), 4.0f, 5.0f, 6.0f, vec3(7.0f, 8.0f, 9.0f));
467     /// mat4 m4 = mat4.identity; // just an identity matrix
468     /// mat3 m3_3 = mat3(m4); // mat3 m3_3 = mat3.identity
469     /// ---
470     this(Args...)(Args args)
471     {
472         construct!(0)(args);
473     }
474 
475     this(T)(T mat)
476     if (isMatrix!T &&
477         (T.cols >= cols) &&
478         (T.rows >= rows))
479     {
480         _matrix[] = mat._matrix[];
481     }
482 
483     this(T)(T mat)
484     if (isMatrix!T &&
485         (T.cols < cols) &&
486         (T.rows < rows))
487     {
488         makeIdentity();
489         static foreach (r; 0 .. T.rows)
490         {
491             static foreach (c; 0 .. T.cols)
492             {
493                 at(r, c) = mat.at(r, c);
494             }
495         }
496     }
497 
498     this()(E value) { clear(value); }
499 
500     /// Sets all values of the matrix to value (each column in each row will contain this value).
501     void clear()(E value)
502     {
503         static foreach (r; 0 .. rows)
504         {
505             static foreach (c; 0 .. cols)
506             {
507                 at(r,c) = value;
508             }
509         }
510     }
511 
512     static if (rows == cols)
513     {
514         /// Makes the current matrix an identity matrix.
515         void makeIdentity()()
516         {
517             clear(0);
518             static foreach (r; 0 .. rows)
519             {
520                 at(r,r) = 1;
521             }
522         }
523 
524         /// Returns: Identity Matrix.
525         static @property Matrix identity()
526         {
527             Matrix ret;
528             ret.clear(0);
529             static foreach (r; 0 .. rows)
530             {
531                 ret.at(r,r) = 1;
532             }
533 
534             return ret;
535         }
536         alias id = identity;    // shorthand
537 
538         /// Transpose Current Matrix.
539         void transpose()()
540         {
541             _matrix = transposed()._matrix;
542         }
543         alias T = transpose; // C++ Armadillo naming convention.
544 
545         unittest
546         {
547             mat2 m2 = mat2(1.0f);
548             m2.transpose();
549             assert(m2._matrix == mat2(1.0f)._matrix);
550             m2.makeIdentity();
551             assert(m2._matrix == [[1.0f, 0.0f],
552                                   [0.0f, 1.0f]]);
553             m2.transpose();
554             assert(m2._matrix == [[1.0f, 0.0f],
555                                   [0.0f, 1.0f]]);
556             assert(m2._matrix == m2.identity._matrix);
557 
558             mat3 m3 = mat3(1.1f, 1.2f, 1.3f,
559                            2.1f, 2.2f, 2.3f,
560                            3.1f, 3.2f, 3.3f);
561             m3.transpose();
562             assert(m3._matrix == [[1.1f, 2.1f, 3.1f],
563                                   [1.2f, 2.2f, 3.2f],
564                                   [1.3f, 2.3f, 3.3f]]);
565 
566             mat4 m4 = mat4(2.0f);
567             m4.transpose();
568             assert(m4._matrix == mat4(2.0f)._matrix);
569             m4.makeIdentity();
570             assert(m4._matrix == [[1.0f, 0.0f, 0.0f, 0.0f],
571                                   [0.0f, 1.0f, 0.0f, 0.0f],
572                                   [0.0f, 0.0f, 1.0f, 0.0f],
573                                   [0.0f, 0.0f, 0.0f, 1.0f]]);
574             assert(m4._matrix == m4.identity._matrix);
575         }
576 
577     }
578 
579     /// Returns: a transposed copy of the matrix.
580     @property Matrix!(E, cols, rows) transposed()() const
581     {
582         typeof(return) ret;
583         static foreach (r; 0 .. rows)
584         {
585             static foreach (c; 0 .. cols)
586             {
587                 ret.at(c,r) = at(r,c);
588             }
589         }
590         return ret;
591     }
592 
593 }
594 alias mat2i = Matrix!(int, 2, 2);
595 alias mat2 = Matrix!(float, 2, 2);
596 alias mat2d = Matrix!(real, 2, 2);
597 alias mat2r = Matrix!(real, 2, 2);
598 alias mat3 = Matrix!(float, 3, 3);
599 alias mat34 = Matrix!(float, 3, 4);
600 alias mat4 = Matrix!(float, 4, 4);
601 alias mat2_cm = Matrix!(float, 2, 2, Layout.columnMajor);
602 
603 @safe pure nothrow @nogc unittest
604 {
605     auto m = mat2(1, 2,
606                   3, 4);
607     assert(m(0, 0) == 1);
608     assert(m(0, 1) == 2);
609     assert(m(1, 0) == 3);
610     assert(m(1, 1) == 4);
611 }
612 
613 @safe pure nothrow @nogc unittest
614 {
615     auto m = mat2_cm(1, 3,
616                      2, 4);
617     assert(m(0, 0) == 1);
618     assert(m(0, 1) == 2);
619     assert(m(1, 0) == 3);
620     assert(m(1, 1) == 4);
621 }
622 
623 @safe pure nothrow @nogc unittest
624 {
625     alias E = float;
626     immutable a = Vector!(E, 2, false, Orient.column)(1, 2);
627     immutable b = Vector!(E, 3, false, Orient.column)(3, 4, 5);
628     immutable c = outerProduct(a, b);
629     assert(c[] == [[3, 4, 5].s,
630                    [6, 8, 10].s].s);
631 }
632 
633 /// 3-Dimensional Spherical Point with Coordinate Type (Precision) `E`.
634 struct SpherePoint3(E)
635 if (isFloatingPoint!E)
636 {
637     enum D = 3;                 // only in three dimensions
638     alias ElementType = E;
639 
640     /** Construct from Components `args`. */
641     this(T...)(T args)
642     {
643         foreach (immutable ix, arg; args)
644         {
645             _spherePoint[ix] = arg;
646         }
647     }
648     /** Element data. */
649     E[D] _spherePoint;
650     enum dimension = D;
651 
652     @property void toString(scope void delegate(scope const(char)[]) @safe sink) const
653     {
654         sink(`SpherePoint3(`);
655         foreach (const ix, const e ; _spherePoint)
656         {
657             if (ix != 0) { sink(","); }
658             sink(to!string(e));
659         }
660         sink(`)`);
661     }
662 
663     @property string toMathML()() const
664     {
665         // opening
666         string str = `<math><mrow>
667   <mo>(</mo>
668   <mtable>`;
669 
670         static foreach (i; 0 .. D)
671         {
672             str ~= `
673     <mtr>
674       <mtd>
675         <mn>` ~ _spherePoint[i].toMathML ~ `</mn>
676       </mtd>
677     </mtr>`;
678         }
679 
680         // closing
681         str ~= `
682   </mtable>
683   <mo>)</mo>
684 </mrow></math>
685 `;
686         return str;
687     }
688 
689     /** Returns: Area 0 */
690     @property E area() const { return 0; }
691 
692     auto opSlice() { return _spherePoint[]; }
693 }
694 
695 /** Instantiator for `SpherePoint3`. */
696 auto spherePoint(Ts...)(Ts args)
697 if (!is(CommonType!Ts == void))
698 {
699     return SpherePoint3!(CommonType!Ts, args.length)(args);
700 }
701 
702 /** `D`-Dimensional Particle with Cartesian Coordinate `position` and
703     `velocity` of Element (Component) Type `E`.
704 */
705 struct Particle(E, uint D,
706                 bool normalizedVelocityFlag = false)
707 if (D >= 1)
708 {
709     Point!(E, D) position;
710     Vector!(E, D, normalizedVelocityFlag) velocity;
711     E mass;
712 }
713 
714 // mixin(makeInstanceAliases("Particle", "particle", 2,4,
715 //                           ["float", "double", "real"]));
716 
717 /** `D`-Dimensional Particle with Coordinate Position and
718     Direction/Velocity/Force Type (Precision) `E`.
719     F = m*a; where F is force, m is mass, a is acceleration.
720 */
721 struct ForcedParticle(E, uint D,
722                       bool normalizedVelocityFlag = false)
723 if (D >= 1)
724 {
725     Point!(E, D) position;
726     Vector!(E, D, normalizedVelocityFlag) velocity;
727     Vector!(E, D) force;
728     E mass;
729 
730     /// Get acceleration.
731     @property auto acceleration()() const { return force/mass; }
732 }
733 
734 /** `D`-Dimensional Axis-Aligned (Hyper) Cartesian `Box` with Element (Component) Type `E`.
735 
736     Note: We must use inclusive compares betweeen boxes and points in inclusion
737     functions such as inside() and includes() in order for the behaviour of
738     bounding boxes (especially in integer space) to work as desired.
739  */
740 struct Box(E, uint D)
741 if (D >= 1)
742 {
743     this(Vector!(E,D) lh) { min = lh; max = lh; }
744     this(Vector!(E,D) l_,
745          Vector!(E,D) h_) { min = l_; max = h_; }
746 
747     @property void toString(scope void delegate(scope const(char)[]) @safe sink) const
748     {
749         sink(`Box(lower:`);
750         min.toString(sink);
751         sink(`, upper:`);
752         max.toString(sink);
753         sink(`)`);
754     }
755 
756     /// Get Box Center.
757     // TODO @property Vector!(E,D) center() { return (min + max) / 2;}
758 
759     /// Constructs a Box enclosing `points`.
760     static Box fromPoints(in Vector!(E,D)[] points)
761     {
762         Box y;
763         foreach (p; points)
764         {
765             y.expand(p);
766         }
767         return y;
768     }
769 
770     /// Expands the Box, so that $(I v) is part of the Box.
771     ref Box expand(Vector!(E,D) v)
772     {
773         static foreach (i; 0 .. D)
774         {
775             if (min[i] > v[i]) min[i] = v[i];
776             if (max[i] < v[i]) max[i] = v[i];
777         }
778         return this;
779     }
780 
781     /// Expands box by another box `b`.
782     ref Box expand()(Box b)
783     {
784         return this.expand(b.min).expand(b.max);
785     }
786 
787     unittest
788     {
789         immutable auto b = Box(Vector!(E,D)(1),
790                                Vector!(E,D)(3));
791         assert(b.sides == Vector!(E,D)(2));
792         immutable auto c = Box(Vector!(E,D)(0),
793                                Vector!(E,D)(4));
794         assert(c.sides == Vector!(E,D)(4));
795         assert(c.sidesProduct == 4^^D);
796         assert(unite(b, c) == c);
797     }
798 
799     /** Returns: Length of Sides */
800     @property auto sides()() const { return max - min; }
801 
802     /** Returns: Area */
803     @property real sidesProduct()() const
804     {
805         typeof(return) y = 1;
806         foreach (const ref side; sides)
807         {
808             y *= side;
809         }
810         return y;
811     }
812 
813     static if (D == 2)
814     {
815         alias area = sidesProduct;
816     }
817     else static if (D == 3)
818     {
819         alias volume = sidesProduct;
820     }
821     else static if (D >= 4)
822     {
823         alias hyperVolume = sidesProduct;
824     }
825 
826     alias include = expand;
827 
828     Vector!(E,D) min;           /// Low.
829     Vector!(E,D) max;           /// High.
830 
831     /** Either an element in min or max is nan or min <= max. */
832     invariant()
833     {
834         // assert(any!"a==a.nan"(min),
835         //                  all!"a || a == a.nan"(elementwiseLessThanOrEqual(min, max)[]));
836     }
837 }
838 
839 // mixin(makeInstanceAliases("Box","box", 2,4,
840 //                           ["int", "float", "double", "real"]));
841 
842 Box!(E,D) unite(E, uint D)(Box!(E,D) a,
843                            Box!(E,D) b) { return a.expand(b); }
844 Box!(E,D) unite(E, uint D)(Box!(E,D) a,
845                            Vector!(E,D) b) { return a.expand(b); }
846 
847 /** `D`-Dimensional Infinite Cartesian (Hyper)-Plane with Element (Component) Type `E`.
848     See_Also: http://stackoverflow.com/questions/18600328/preferred-representation-of-a-3d-plane-in-c-c
849  */
850 struct Plane(E, uint D)
851 if (D >= 2 &&
852     isFloatingPoint!E)
853 {
854     enum dimension = D;
855 
856     alias ElementType = E;
857 
858     /// Normal type of plane.
859     alias NormalType = Vector!(E, D, true);
860 
861     union
862     {
863         static if (D == 3)
864         {
865             struct
866             {
867                 E a; /// normal.x
868                 E b; /// normal.y
869                 E c; /// normal.z
870             }
871         }
872         NormalType normal;      /// Plane Normal.
873     }
874     E distance;                  /// Plane Constant (Offset from origo).
875 
876     @property void toString(scope void delegate(scope const(char)[]) @safe sink) const
877     {
878         import std.conv : to;
879         sink(`Plane(normal:`);
880         sink(to!string(normal));
881         sink(`, distance:`);
882         sink(to!string(distance));
883         sink(`)`);
884     }
885 
886     /// Constructs the plane, from either four scalars of type $(I E)
887     /// or from a 3-dimensional vector (= normal) and a scalar.
888     static if (D == 2)
889     {
890         this()(E a, E b, E distance)
891         {
892             this.normal.x = a;
893             this.normal.y = b;
894             this.distance = distance;
895         }
896     }
897     else static if (D == 3)
898     {
899         this()(E a, E b, E c, E distance)
900         {
901             this.normal.x = a;
902             this.normal.y = b;
903             this.normal.z = c;
904             this.distance = distance;
905         }
906     }
907 
908     this()(NormalType normal, E distance)
909     {
910         this.normal = normal;
911         this.distance = distance;
912     }
913 
914     /* unittest
915        { */
916     /*     Plane p = Plane(0.0f, 1.0f, 2.0f, 3.0f); */
917     /*     assert(p.normal == N(0.0f, 1.0f, 2.0f)); */
918     /*     assert(p.distance == 3.0f); */
919 
920     /*     p.normal.x = 4.0f; */
921     /*     assert(p.normal == N(4.0f, 1.0f, 2.0f)); */
922     /*     assert(p.x == 4.0f); */
923     /*     assert(p.y == 1.0f); */
924     /*     assert(p.c == 2.0f); */
925     /*     assert(p.distance == 3.0f); */
926     /* } */
927 
928     /// Normalizes the plane inplace.
929     void normalize()()
930     {
931         immutable E det = cast(E)1 / normal.magnitude;
932         normal *= det;
933         distance *= det;
934     }
935 
936     /// Returns: a normalized copy of the plane.
937     /* @property Plane normalized() const { */
938     /*     Plane y = Plane(a, b, c, distance); */
939     /*     y.normalize(); */
940     /*     return y; */
941     /* } */
942 
943 //     unittest {
944 //         Plane p = Plane(0.0f, 1.0f, 2.0f, 3.0f);
945 //         Plane pn = p.normalized();
946 //         assert(pn.normal == N(0.0f, 1.0f, 2.0f).normalized);
947 //         assert(almost_equal(pn.distance, 3.0f / N(0.0f, 1.0f, 2.0f).length));
948 //         p.normalize();
949 //         assert(p == pn);
950 //     }
951 
952     /// Returns: distance from a point to the plane.
953     /// Note: the plane $(RED must) be normalized, the result can be negative.
954     /* E distanceTo(N point) const { */
955     /*     return dot(point, normal) + distance; */
956     /* } */
957 
958     /// Returns: distanceTo from a point to the plane.
959     /// Note: the plane does not have to be normalized, the result can be negative.
960     /* E ndistance(N point) const { */
961     /*     return (dot(point, normal) + distance) / normal.magnitude; */
962     /* } */
963 
964     /* unittest
965        { */
966     /*     Plane p = Plane(-1.0f, 4.0f, 19.0f, -10.0f); */
967     /*     assert(almost_equal(p.ndistance(N(5.0f, -2.0f, 0.0f)), -1.182992)); */
968     /*     assert(almost_equal(p.ndistance(N(5.0f, -2.0f, 0.0f)), */
969     /*                         p.normalized.distanceTo(N(5.0f, -2.0f, 0.0f)))); */
970     /* } */
971 
972     /* bool opEquals(Plane other) const { */
973     /*     return other.normal == normal && other.distance == distance; */
974     /* } */
975 
976 }
977 
978 // mixin(makeInstanceAliases("Plane","plane", 3,4,
979 //                           defaultElementTypes));
980 
981 /** `D`-Dimensional Cartesian (Hyper)-Sphere with Element (Component) Type `E`.
982  */
983 struct Sphere(E, uint D)
984 if (D >= 2 &&
985     isNumeric!E)
986 {
987     alias CenterType = Point!(E, D);
988 
989     CenterType center;
990     E radius;
991 
992     void translate(Vector!(E, D) shift)
993     {
994         center = center + shift; // point + vector => point
995     }
996     alias shift = translate;
997 
998     @property:
999 
1000     E diameter()() const
1001     {
1002         return 2 * radius;
1003     }
1004     static if (D == 2)
1005     {
1006         auto area()() const
1007         {
1008             return PI * radius ^^ 2;
1009         }
1010     }
1011     else static if (D == 3)
1012     {
1013         auto area()() const
1014         {
1015             return 4 * PI * radius ^^ 2;
1016         }
1017         auto volume()() const
1018         {
1019             return 4 * PI * radius ^^ 3 / 3;
1020         }
1021     }
1022     else static if (D >= 4)
1023     {
1024         // See_Also: https://en.wikipedia.org/wiki/Volume_of_an_n-ball
1025         real n = D;
1026         auto volume()() const
1027         {
1028             import std.mathspecial: gamma;
1029             return PI ^^ (n / 2) / gamma(n / 2 + 1) * radius ^^ n;
1030         }
1031     }
1032 }
1033 
1034 auto sphere(C, R)(C center, R radius)
1035 {
1036     return Sphere!(C.ElementType, C.dimension)(center, radius);
1037 }
1038 // TODO Use this instead:
1039 // auto sphere(R, C...)(Point!(CommonType!C, C.length) center, R radius) {
1040 // return Sphere!(CommonType!C, C.length)(center, radius);
1041 // }
1042 
1043 /**
1044    See_Also: http://stackoverflow.com/questions/401847/circle-rectangle-collision-detection-intersect
1045  */
1046 bool intersect(T)(Circle!T circle, Rect!T rect)
1047 {
1048     immutable hw = rect.w/2, hh = rect.h/2;
1049 
1050     immutable dist = Point!T(abs(circle.x - rect.x0 - hw),
1051                              abs(circle.y - rect.y0 - hh));
1052 
1053     if (dist.x > (hw + circle.r)) return false;
1054     if (dist.y > (hh + circle.r)) return false;
1055 
1056     if (dist.x <= hw) return true;
1057     if (dist.y <= hh) return true;
1058 
1059     immutable cornerDistance_sq = ((dist.x - hw)^^2 +
1060                                    (dist.y - hh)^^2);
1061 
1062     return (cornerDistance_sq <= circle.r^^2);
1063 }
1064 
1065 @safe unittest
1066 {
1067     assert(box2f(vec2f(1, 2), vec2f(3, 3)).to!string == `Box(lower:ColumnVector(1,2), upper:ColumnVector(3,3))`);
1068     assert([12, 3, 3].to!string == `[12, 3, 3]`);
1069 
1070     assert(vec2f(2, 3).to!string == `ColumnVector(2,3)`);
1071 
1072     assert(vec2f(2, 3).to!string == `ColumnVector(2,3)`);
1073     assert(vec2f(2, 3).to!string == `ColumnVector(2,3)`);
1074 
1075     assert(vec3f(2, 3, 4).to!string == `ColumnVector(2,3,4)`);
1076 
1077     assert(box2f(vec2f(1, 2),
1078                  vec2f(3, 4)).to!string == `Box(lower:ColumnVector(1,2), upper:ColumnVector(3,4))`);
1079 
1080     assert(vec2i(2, 3).to!string == `ColumnVector(2,3)`);
1081     assert(vec3i(2, 3, 4).to!string == `ColumnVector(2,3,4)`);
1082     assert(vec3i(2, 3, 4).to!string == `ColumnVector(2,3,4)`);
1083 
1084     assert(vec2i(2, 3).toMathML == `<math><mrow>
1085   <mo>⟨</mo>
1086   <mtable>
1087     <mtr>
1088       <mtd>
1089         <mn>2</mn>
1090       </mtd>
1091     </mtr>
1092     <mtr>
1093       <mtd>
1094         <mn>3</mn>
1095       </mtd>
1096     </mtr>
1097   </mtable>
1098   <mo>⟩</mo>
1099 </mrow></math>
1100 `);
1101 
1102     auto m = mat2(1, 2, 3, 4);
1103     assert(m.toLaTeX == `\begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}`);
1104     assert(m.toMathML == `<math><mrow>
1105   <mo>❲</mo>
1106   <mtable>
1107     <mtr>
1108       <mtd>
1109         <mn>1</mn>
1110       </mtd>
1111       <mtd>
1112         <mn>2</mn>
1113       </mtd>
1114     </mtr>
1115     <mtr>
1116       <mtd>
1117         <mn>3</mn>
1118       </mtd>
1119       <mtd>
1120         <mn>4</mn>
1121       </mtd>
1122     </mtr>
1123   </mtable>
1124   <mo>❳</mo>
1125 </mrow></math>
1126 `);
1127 }
1128 
1129 version(unittest)
1130 {
1131     import std.conv : to;
1132     import nxt.array_help : s;
1133 }