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 }