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;
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 		is(T == Vector!(_), _) &&
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 (!is(T == Vector!(_), _) &&
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 	pure nothrow @safe @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 pure nothrow @safe @nogc unittest {
669 	assert(point(1, 2) + vector(1, 2) == point(2, 4));
670 	assert(point(1, 2) - vector(1, 2) == point(0, 0));
671 }
672 
673 version (unittest)
674 {
675 	import std.conv : to;
676 	import nxt.array_help : s;
677 }