nxt.geometry source code

1 /** Geometry and linear algebra primitives.
2 
3    TODO: get rid of instantiators in favor of aliases to reduce template bloat
4    TODO: Merge https://github.com/andrewlalis/dvec
5 */
6 module nxt.geometry;
7 
8 import std.traits: isFloatingPoint, isNumeric, isSigned, isDynamicArray, isAssignable, isArray, CommonType, isInstanceOf;
9 
10 version = use_cconv;            ///< Use cconv for faster compilation.
11 version = unittestAllInstances;
12 
13 version(unittestAllInstances)
14     static immutable defaultElementTypes = ["float", "double", "real"];
15 else
16     static immutable defaultElementTypes = ["double"];
17 
18 enum Orient { column, row } // Vector Orientation.
19 
20 /** `D`-Dimensional Cartesian Point with Coordinate Type (Precision) `E`.
21  */
22 struct Point(E, uint D)
23 if (D >= 1 /* && TODO: extend trait : isNumeric!E */)
24 {
25     alias ElementType = E;
26 
27     this(T...)(T args)
28     {
29         foreach (immutable ix, arg; args)
30             _point[ix] = arg;
31     }
32 
33     /** Element data. */
34     E[D] _point;
35     enum dimension = D;
36 
37     void toString(Sink)(ref scope Sink sink) const
38     {
39         sink(`Point(`);
40         foreach (const ix, const e; _point)
41         {
42             if (ix != 0) { sink(","); }
43             version(use_cconv)
44             {
45                 import nxt.cconv : toStringInSink;
46                 toStringInSink(e, sink);
47             }
48             else
49             {
50                 import std.conv : to;
51                 sink(to!string(e));
52             }
53         }
54         sink(`)`);
55     }
56 
57     @property string toMathML()() const
58     {
59         // opening
60         string str = `<math><mrow>
61   <mo>(</mo>
62   <mtable>`;
63 
64         static foreach (i; 0 .. D)
65         {
66             str ~= `
67     <mtr>
68       <mtd>
69         <mn>` ~ _point[i].toMathML ~ `</mn>
70       </mtd>
71     </mtr>`;
72         }
73 
74         // closing
75         str ~= `
76   </mtable>
77   <mo>)</mo>
78 </mrow></math>
79 `;
80         return str;
81     }
82 
83     /** Returns: Area 0 */
84     @property E area() const => 0;
85 
86     inout(E)[] opSlice() inout => _point[];
87 
88     /** Points +/- Vector => Point */
89     auto opBinary(string op, F)(Vector!(F, D) r) const
90         if ((op == "+") ||
91             (op == "-"))
92     {
93         Point!(CommonType!(E, F), D) y;
94         static foreach (i; 0 .. D)
95             y._point[i] = mixin("_point[i]" ~ op ~ "r._vector[i]");
96         return y;
97     }
98 }
99 
100 /** Instantiator for `Point`. */
101 auto point(Ts...)(Ts args)
102 if (!is(CommonType!Ts == void))
103 	=> Point!(CommonType!Ts, args.length)(args);
104 
105 version(unittest)
106 {
107     alias vec2f = Vector!(float, 2, true);
108     alias vec3f = Vector!(float, 3, true);
109 
110     alias vec2d = Vector!(float, 2, true);
111 
112     alias nvec2f = Vector!(float, 2, true);
113 }
114 
115 /** `D`-Dimensional Vector with Coordinate/Element (Component) Type `E`.
116  *
117  * See_Also: http://physics.stackexchange.com/questions/16850/is-0-0-0-an-undefined-vector
118  */
119 struct Vector(E, uint D,
120               bool normalizedFlag = false, // `true` for unit vectors
121               Orient orient = Orient.column)
122 if (D >= 1 &&
123     (!normalizedFlag ||
124      isFloatingPoint!E)) // only normalize fp for now
125 {
126     // Construct from vector.
127     this(V)(V vec)
128     if (isVector!V &&
129         // TOREVIEW: is(T.E : E) &&
130         (V.dimension >= dimension))
131     {
132         static if (normalizedFlag)
133         {
134             if (vec.normalized)
135             {
136                 immutable vec_norm = vec.magnitude;
137                 static foreach (i; 0 .. D)
138                     _vector[i] = vec._vector[i] / vec_norm;
139                 return;
140             }
141         }
142         static foreach (i; 0 .. D)
143             _vector[i] = vec._vector[i];
144     }
145 
146     /** Construct from Scalar `value`. */
147     this(S)(S scalar)
148     if (isAssignable!(E, S))
149     {
150         static if (normalizedFlag)
151         {
152             import std.math : sqrt;
153             clear(1/sqrt(cast(E)D)); // TODO: costly
154         }
155         else
156             clear(scalar);
157     }
158 
159     /** Construct from combination of arguments. */
160     this(Args...)(Args args) { construct!(0)(args); }
161 
162     enum dimension = D;
163 
164     @property const
165     {
166         string orientationString()()
167 			=> orient == Orient.column ? `Column` : `Row`;
168         string joinString()()
169 			=> orient == Orient.column ? ` \\ ` : ` & `;
170     }
171 
172     void toString(Sink)(ref scope Sink sink) const
173     {
174         sink(orientationString);
175         sink(`Vector(`);
176         foreach (const ix, const e; _vector)
177         {
178             if (ix != 0) { sink(","); }
179             version(use_cconv)
180             {
181                 import nxt.cconv : toStringInSink;
182                 toStringInSink(e, sink);
183             }
184             else
185             {
186                 import std.conv : to;
187                 sink(to!string(e));
188             }
189         }
190         sink(`)`);
191     }
192 
193     /** Returns: LaTeX Encoding of Vector. http://www.thestudentroom.co.uk/wiki/LaTex#Matrices_and_Vectors */
194     @property string toLaTeX()() const
195     	=> `\begin{pmatrix} ` ~ map!(to!string)(_vector[]).join(joinString) ~ ` \end{pmatrix}` ;
196 
197     @property string toMathML()() const
198     {
199         // opening
200         string str = `<math><mrow>
201   <mo>⟨</mo>
202   <mtable>`;
203 
204         if (orient == Orient.row)
205         {
206             str ~=  `
207     <mtr>`;
208         }
209 
210         static foreach (i; 0 .. D)
211         {
212             final switch (orient)
213             {
214             case Orient.column:
215                 str ~= `
216     <mtr>
217       <mtd>
218         <mn>` ~ _vector[i].toMathML ~ `</mn>
219       </mtd>
220     </mtr>`;
221                 break;
222             case Orient.row:
223                 str ~= `
224       <mtd>
225         <mn>` ~ _vector[i].toMathML ~ `</mn>
226       </mtd>`;
227                 break;
228             }
229         }
230 
231         if (orient == Orient.row)
232         {
233             str ~=  `
234     </mtr>`;
235         }
236 
237         // closing
238         str ~= `
239   </mtable>
240   <mo>⟩</mo>
241 </mrow></math>
242 `;
243         return str;
244     }
245 
246     auto randInPlace()(E scaling = 1)
247     {
248         import nxt.random_ex: randInPlace;
249         static if (normalizedFlag &&
250                    isFloatingPoint!E) // cannot use normalized() here (yet)
251         {
252             static if (D == 2)  // randomize on unit circle
253             {
254                 alias P = real; // precision
255                 immutable angle = uniform(0, 2*cast(P)PI);
256                 _vector[0] = scaling*sin(angle);
257                 _vector[1] = scaling*cos(angle);
258             }
259             static if (D == 3)  // randomize on unit sphere: See_Also: http://mathworld.wolfram.com/SpherePointPicking.html
260             {
261                 alias P = real; // precision
262                 immutable u = uniform(0, cast(P)1);
263                 immutable v = uniform(0, cast(P)1);
264                 immutable theta = 2*PI*u;
265                 immutable phi = acos(2*v-1);
266                 _vector[0] = scaling*cos(theta)*sin(phi);
267                 _vector[1] = scaling*sin(theta)*sin(phi);
268                 _vector[2] = scaling*cos(phi);
269             }
270             else
271             {
272                 _vector.randInPlace();
273                 normalize(); // TODO: Turn this into D data restriction instead?
274             }
275         }
276         else
277         {
278             _vector.randInPlace();
279         }
280     }
281 
282     /// Returns: `true` if all values are not `nan` nor `infinite`, otherwise `false`.
283     @property bool ok()() const
284     {
285         static if (isFloatingPoint!E)
286         {
287             foreach (const ref v; _vector)
288                 if (isNaN(v) ||
289                     isInfinity(v))
290                     return false;
291         }
292         return true;
293     }
294     // NOTE: Disabled this because I want same behaviour as MATLAB: bool opCast(T : bool)() const { return ok; }
295     bool opCast(T : bool)() const => all!"a"(_vector[]);
296 
297     /// Returns: Pointer to the coordinates.
298     // @property auto value_ptr() { return _vector.ptr; }
299 
300     /// Sets all values to `value`.
301     void clear(V)(V value)
302     if (isAssignable!(E, V))
303     {
304         static foreach (i; 0 .. D)
305             _vector[i] = value;
306     }
307 
308     /** Returns: Whole Internal Array of E. */
309     auto opSlice() => _vector[];
310     /** Returns: Slice of Internal Array of E. */
311     auto opSlice(uint off, uint len) => _vector[off .. len];
312     /** Returns: Reference to Internal Vector Element. */
313     ref inout(E) opIndex(uint i) inout => _vector[i];
314 
315     bool opEquals(S)(in S scalar) const
316     if (isAssignable!(E, S)) // TODO: is(typeof(E.init != S.init))
317 		=> _vector[] == scalar;
318 
319     bool opEquals(F)(in F vec) const
320     if (isVector!F &&
321         dimension == F.dimension) // TOREVIEW: Use isEquable instead?
322 		=> _vector == vec._vector;
323 
324     bool opEquals(F)(const F[] array) const
325     if (isAssignable!(E, F) &&
326         !isArray!F &&
327         !isVector!F) // TOREVIEW: Use isNotEquable instead?
328     	=> _vector[] == array;
329 
330     static void isCompatibleVectorImpl(uint d)(Vector!(E, d) vec) if (d <= dimension) {}
331     static void isCompatibleMatrixImpl(uint r, uint c)(Matrix!(E, r, c) m) {}
332 
333     enum isCompatibleVector(T) = is(typeof(isCompatibleVectorImpl(T.init)));
334     enum isCompatibleMatrix(T) = is(typeof(isCompatibleMatrixImpl(T.init)));
335 
336     private void construct(uint i)()
337     {
338         static assert(i == D, "Not enough arguments passed to constructor");
339     }
340     private void construct(uint i, T, Tail...)(T head, Tail tail)
341     {
342         static        if (i >= D)
343         {
344             static assert(0, "Too many arguments passed to constructor");
345         }
346         else static if (is(T : E))
347         {
348             _vector[i] = head;
349             construct!(i + 1)(tail);
350         }
351         else static if (isDynamicArray!T)
352         {
353             static assert((Tail.length == 0) && (i == 0), "Dynamic array can not be passed together with other arguments");
354             _vector[] = head[];
355         }
356         else static if (__traits(isStaticArray, T))
357         {
358             _vector[i .. i + T.length] = head[];
359             construct!(i + T.length)(tail);
360         }
361         else static if (isCompatibleVector!T)
362         {
363             _vector[i .. i + T.dimension] = head._vector[];
364             construct!(i + T.dimension)(tail);
365         }
366         else
367             static assert(0, "Vector constructor argument must be of type " ~ E.stringof ~ " or Vector, not " ~ T.stringof);
368     }
369 
370     // private void dispatchImpl(int i, string s, int size)(ref E[size] result) const {
371     //     static if (s.length > 0) {
372     //         result[i] = _vector[coordToIndex!(s[0])];
373     //         dispatchImpl!(i + 1, s[1..$])(result);
374     //     }
375     // }
376 
377     // /// Implements dynamic swizzling.
378     // /// Returns: a Vector
379     // @property Vector!(E, s.length) opDispatch(string s)() const {
380     //     E[s.length] ret;
381     //     dispatchImpl!(0, s)(ret);
382     //     Vector!(E, s.length) ret_vec;
383     //     ret_vec._vector = ret;
384     //     return ret_vec;
385     // }
386 
387     ref inout(Vector) opUnary(string op : "+")() inout => this;
388 
389     Vector opUnary(string op : "-")() const
390     if (isSigned!(E))
391     {
392         Vector y;
393         static foreach (i; 0 .. D)
394             y._vector[i] = - _vector[i];
395         return y;
396     }
397 
398     auto opBinary(string op, T)(in T rhs) const
399     if (op == "+" || op == "-" &&
400         isInstanceOf!(Vector, T) &&
401         dimension && T.dimension)
402     {
403         Vector!(CommonType!(E, F), D) y;
404         static foreach (i; 0 .. D)
405             y._vector[i] = mixin("_vector[i]" ~ op ~ "rhs._vector[i]");
406         return y;
407     }
408 
409     Vector opBinary(string op : "*", F)(in F r) const
410     {
411         Vector!(CommonType!(E, F), D, normalizedFlag) y;
412         static foreach (i; 0 .. dimension)
413             y._vector[i] = _vector[i] * r;
414         return y;
415     }
416 
417     Vector!(CommonType!(E, F), D) opBinary(string op : "*", F)(in Vector!(F, D) r) const
418     {
419         // MATLAB-style Product Behaviour
420         static if (orient == Orient.column &&
421                    r.orient == Orient.row)
422             return outer(this, r);
423         else static if (orient == Orient.row &&
424                         r.orient == Orient.column)
425             return dot(this, r);
426         else
427             static assert(0, "Incompatible vector dimensions.");
428     }
429 
430     /** Multiply with right-hand-side `rhs`. */
431     Vector!(E, T.rows) opBinary(string op : "*", T)(in T rhs) const
432     if (isCompatibleMatrix!T &&
433         (T.cols == dimension))
434     {
435         Vector!(E, T.rows) ret;
436         ret.clear(0);
437         static foreach (c; 0 .. T.cols)
438             static foreach (r; 0 ..  T.rows)
439                 ret._vector[r] += _vector[c] * rhs.at(r,c);
440         return ret;
441     }
442 
443     /** Multiply with left-hand-side `lhs`. */
444     version(none)               // TODO: activate
445     auto opBinaryRight(string op, T)(in T lhs) const
446     if (!isInstanceOf!(Vector, T) &&
447         !isMatrix!T &&
448         !isQuaternion!T)
449 		=> this.opBinary!(op)(lhs);
450 
451     /** TODO: Suitable Restrictions on F. */
452     void opOpAssign(string op, F)(in F r)
453         /* if ((op == "+") || (op == "-") || (op == "*") || (op == "%") || (op == "/") || (op == "^^")) */
454     {
455         static foreach (i; 0 .. dimension)
456             mixin("_vector[i]" ~ op ~ "= r;");
457     }
458     unittest
459     {
460         auto v2 = vec2f(1, 3);
461         v2 *= 5.0f; assert(v2[] == [5, 15].s);
462         v2 ^^= 2; assert(v2[] == [25, 225].s);
463         v2 /= 5; assert(v2[] == [5, 45].s);
464     }
465 
466     void opOpAssign(string op)(in Vector r)
467     if ((op == "+") ||
468         (op == "-"))
469     {
470         static foreach (i; 0 .. dimension)
471             mixin("_vector[i]" ~ op ~ "= r._vector[i];");
472     }
473 
474     /// Returns: Non-Rooted `N` - Norm of `x`.
475     auto nrnNorm(uint N)() const
476     if (isNumeric!E && N >= 1)
477     {
478         static if (isFloatingPoint!E)
479             real y = 0;                 // TOREVIEW: Use maximum precision for now
480         else
481             E y = 0;                // TOREVIEW: Use other precision for now
482         static foreach (i; 0 .. D)
483             y += _vector[i] ^^ N;
484         return y;
485     }
486 
487     /// Returns: Squared Magnitude of x.
488     @property real magnitudeSquared()() const
489     if (isNumeric!E)
490     {
491         static if (normalizedFlag) // cannot use normalized() here (yet)
492             return 1;
493         else
494             return nrnNorm!2;
495     }
496 
497     /// Returns: Magnitude of x.
498     @property real magnitude()() const
499     if (isNumeric!E)
500     {
501         static if (normalizedFlag) // cannot use normalized() here (yet)
502             return 1;
503         else
504         {
505             import std.math : sqrt;
506             return sqrt(magnitudeSquared);
507         }
508     }
509     alias norm = magnitude;
510 
511     static if (isFloatingPoint!E)
512     {
513         /// Normalize `this`.
514         void normalize()()
515         {
516             if (this != 0)         // zero vector have zero magnitude
517             {
518                 immutable m = this.magnitude;
519                 static foreach (i; 0 .. D)
520                     _vector[i] /= m;
521             }
522         }
523 
524         /// Returns: normalizedFlag Copy of this Vector.
525         @property Vector normalized()() const
526         {
527             Vector y = this;
528             y.normalize();
529             return y;
530         }
531 
532         unittest
533         {
534             static if (D == 2 && !normalizedFlag)
535                 assert(Vector(3, 4).magnitude == 5);
536         }
537     }
538 
539     /// Returns: Vector Index at Character Coordinate `coord`.
540     private @property ref inout(E) get_(char coord)() inout
541     {
542         return _vector[coordToIndex!coord];
543     }
544 
545     /// Coordinate Character c to Index
546     template coordToIndex(char c)
547     {
548         static if ((c == 'x'))
549             enum coordToIndex = 0;
550         else static if ((c == 'y'))
551             enum coordToIndex = 1;
552         else static if ((c == 'z'))
553         {
554             static assert(D >= 3, "The " ~ c ~ " property is only available on vectors with a third dimension.");
555             enum coordToIndex = 2;
556         }
557         else static if ((c == 'w'))
558         {
559             static assert(D >= 4, "The " ~ c ~ " property is only available on vectors with a fourth dimension.");
560             enum coordToIndex = 3;
561         }
562         else
563             static assert(0, "Accepted coordinates are x, s, r, u, y, g, t, v, z, p, b, w, q and a not " ~ c ~ ".");
564     }
565 
566     /// Updates the vector with the values from other.
567     void update(in Vector!(E, D) other) { _vector = other._vector; }
568 
569     static if (D == 2)
570     {
571         void set(E x, E y)
572         {
573             _vector[0] = x;
574             _vector[1] = y;
575         }
576     }
577     else static if (D == 3)
578     {
579         void set(E x, E y, E z)
580         {
581             _vector[0] = x;
582             _vector[1] = y;
583             _vector[2] = z;
584         }
585     }
586     else static if (D == 4)
587     {
588         void set(E x, E y, E z, E w)
589         {
590             _vector[0] = x;
591             _vector[1] = y;
592             _vector[2] = z;
593             _vector[3] = w;
594         }
595     }
596 
597     static if (D >= 1) { alias x = get_!'x'; }
598     static if (D >= 2) { alias y = get_!'y'; }
599     static if (D >= 3) { alias z = get_!'z'; }
600     static if (D >= 4) { alias w = get_!'w'; }
601 
602     static if (isNumeric!E)
603     {
604         /* Need these conversions when E is for instance ubyte.
605            See this commit: https://github.com/Dav1dde/gl3n/commit/2504003df4f8a091e58a3d041831dc2522377f95 */
606         static immutable E0 = E(0);
607         static immutable E1 = E(1);
608         static if (dimension == 2)
609         {
610             static immutable Vector e1 = Vector(E1, E0); /// canonical basis for Euclidian space
611             static immutable Vector e2 = Vector(E0, E1); /// ditto
612         }
613         else static if (dimension == 3)
614         {
615             static immutable Vector e1 = Vector(E1, E0, E0); /// canonical basis for Euclidian space
616             static immutable Vector e2 = Vector(E0, E1, E0); /// ditto
617             static immutable Vector e3 = Vector(E0, E0, E1); /// ditto
618         }
619         else static if (dimension == 4)
620         {
621             static immutable Vector e1 = Vector(E1, E0, E0, E0); /// canonical basis for Euclidian space
622             static immutable Vector e2 = Vector(E0, E1, E0, E0); /// ditto
623             static immutable Vector e3 = Vector(E0, E0, E1, E0); /// ditto
624             static immutable Vector e4 = Vector(E0, E0, E0, E1); /// ditto
625         }
626     }
627 
628     version(none)
629     @safe pure nothrow @nogc unittest
630     {
631         static if (isNumeric!E)
632         {
633             assert(vec2.e1[] == [1, 0].s);
634             assert(vec2.e2[] == [0, 1].s);
635 
636             assert(vec3.e1[] == [1, 0, 0].s);
637             assert(vec3.e2[] == [0, 1, 0].s);
638             assert(vec3.e3[] == [0, 0, 1].s);
639 
640             assert(vec4.e1[] == [1, 0, 0, 0].s);
641             assert(vec4.e2[] == [0, 1, 0, 0].s);
642             assert(vec4.e3[] == [0, 0, 1, 0].s);
643             assert(vec4.e4[] == [0, 0, 0, 1].s);
644         }
645     }
646 
647     /** Element data. */
648     E[D] _vector;
649 
650     unittest
651     {
652         // static if (isSigned!(E)) { assert(-Vector!(E,D)(+2),
653         //                                   +Vector!(E,D)(-2)); }
654     }
655 
656 }
657 
658 auto rowVector(Ts...)(Ts args)
659 if (!is(CommonType!Ts == void))
660 	=> Vector!(CommonType!Ts, args.length)(args);
661 alias vector = rowVector; // TODO: Should rowVector or columnVector be default?
662 
663 auto columnVector(Ts...)(Ts args)
664 if (!is(CommonType!Ts == void))
665 	=> Vector!(CommonType!Ts, args.length, false, Orient.column)(args);
666 
667 ///
668 @safe pure nothrow @nogc unittest
669 {
670     assert(point(1, 2) + vector(1, 2) == point(2, 4));
671     assert(point(1, 2) - vector(1, 2) == point(0, 0));
672 }
673 
674 version(unittest)
675 {
676     import std.conv : to;
677     import nxt.array_help : s;
678 }