diff --git a/dub.json b/dub.json
index b288d74..bcecac5 100644
--- a/dub.json
+++ b/dub.json
@@ -9,8 +9,8 @@
"targetPath": "build",
"sourceFiles-linux": ["build/libvma.a", "build/libstb_image.a", "build/libm3d.a"],
"sourceFiles-windows": [],
- "importPaths": ["src/gears", "src/shared", "src/generated", "external/xxhash", "external/gfm", "external/inteli"],
- "sourcePaths": ["src/gears", "src/shared", "src/generated", "external/xxhash", "external/gfm", "external/inteli"],
+ "importPaths": ["src/gears", "src/shared", "src/generated", "external/xxhash", "external/inteli"],
+ "sourcePaths": ["src/gears", "src/shared", "src/generated", "external/xxhash", "external/inteli"],
"libs-linux": ["xcb", "X11", "X11-xcb", "vulkan", "stdc++"],
"libs-windows": [],
"preGenerateCommands-linux": ["./build.sh", "build/Packer"],
@@ -22,8 +22,8 @@
"targetType": "executable",
"targetPath": "build",
"targetName": "Packer",
- "importPaths": ["src/packer", "src/shared", "src/generated", "external/xxhash", "external/gfm", "external/inteli"],
- "sourcePaths": ["src/packer", "src/shared", "src/generated", "external/xxhash", "external/gfm", "external/inteli"],
+ "importPaths": ["src/packer", "src/shared", "src/generated", "external/xxhash", "external/inteli"],
+ "sourcePaths": ["src/packer", "src/shared", "src/generated", "external/xxhash", "external/inteli"],
"sourceFiles-linux": ["build/libstb_image.a", "build/libm3d.a"],
"preGenerateCommands-linux": ["./build.sh"],
"postGenerateCommands-linux": [],
@@ -35,8 +35,8 @@
"targetType": "executable",
"targetPath": "build",
"targetName": "Codegen",
- "importPaths": ["src/codegen", "src/shared", "external/xxhash", "external/gfm", "external/inteli"],
- "sourcePaths": ["src/codegen", "src/shared", "external/xxhash", "external/gfm", "external/inteli"],
+ "importPaths": ["src/codegen", "src/shared", "external/xxhash", "external/inteli"],
+ "sourcePaths": ["src/codegen", "src/shared", "external/xxhash", "external/inteli"],
"sourceFiles-linux": ["build/libstb_image.a"],
"preGenerateCommands-linux": ["./build.sh"],
"preGenerateCommands-windows": [],
diff --git a/external/gfm/math/box.d b/external/gfm/math/box.d
deleted file mode 100644
index 839e0b4..0000000
--- a/external/gfm/math/box.d
+++ /dev/null
@@ -1,659 +0,0 @@
-/**
- This module implements a generic axis-aligned bounding box (AABB).
- */
-module gfm.math.box;
-
-import std.math,
- std.traits,
- std.conv,
- std.string;
-
-import gfm.math.vector,
- gfm.math.funcs;
-
-/// N-dimensional half-open interval [a, b[.
-struct Box(T, int N)
-{
- static assert(N > 0);
-
- public
- {
- alias bound_t = Vector!(T, N);
-
- bound_t min; // not enforced, the box can have negative volume
- bound_t max;
-
- /// Construct a box which extends between 2 points.
- /// Boundaries: min is inside the box, max is just outside.
- @nogc this(bound_t min_, bound_t max_) pure nothrow
- {
- min = min_;
- max = max_;
- }
-
- static if (N == 1)
- {
- @nogc this(T min_, T max_) pure nothrow
- {
- min.x = min_;
- max.x = max_;
- }
- }
-
- static if (N == 2)
- {
- @nogc this(T min_x, T min_y, T max_x, T max_y) pure nothrow
- {
- min = bound_t(min_x, min_y);
- max = bound_t(max_x, max_y);
- }
- }
-
- static if (N == 3)
- {
- @nogc this(T min_x, T min_y, T min_z, T max_x, T max_y, T max_z) pure nothrow
- {
- min = bound_t(min_x, min_y, min_z);
- max = bound_t(max_x, max_y, max_z);
- }
- }
-
- @property
- {
- /// Returns: Dimensions of the box.
- @nogc bound_t size() pure const nothrow
- {
- return max - min;
- }
-
- /// Sets size of the box assuming min point is the pivot.
- /// Returns: Dimensions of the box.
- @nogc bound_t size(bound_t value) pure nothrow
- {
- max = min + value;
- return value;
- }
-
- /// Returns: Center of the box.
- @nogc bound_t center() pure const nothrow
- {
- return (min + max) / 2;
- }
-
- static if (N >= 1)
- {
- /// Returns: Width of the box, always applicable.
- @nogc T width() pure const nothrow @property
- {
- return max.x - min.x;
- }
-
- /// Sets width of the box assuming min point is the pivot.
- /// Returns: Width of the box, always applicable.
- @nogc T width(T value) pure nothrow @property
- {
- max.x = min.x + value;
- return value;
- }
- }
-
- static if (N >= 2)
- {
- /// Returns: Height of the box, if applicable.
- @nogc T height() pure const nothrow @property
- {
- return max.y - min.y;
- }
-
- /// Sets height of the box assuming min point is the pivot.
- /// Returns: Height of the box, if applicable.
- @nogc T height(T value) pure nothrow @property
- {
- max.y = min.y + value;
- return value;
- }
- }
-
- static if (N >= 3)
- {
- /// Returns: Depth of the box, if applicable.
- @nogc T depth() pure const nothrow @property
- {
- return max.z - min.z;
- }
-
- /// Sets depth of the box assuming min point is the pivot.
- /// Returns: Depth of the box, if applicable.
- @nogc T depth(T value) pure nothrow @property
- {
- max.z = min.z + value;
- return value;
- }
- }
-
- /// Returns: Signed volume of the box.
- @nogc T volume() pure const nothrow
- {
- T res = 1;
- bound_t size = size();
- for(int i = 0; i < N; ++i)
- res *= size[i];
- return res;
- }
-
- /// Returns: true if empty.
- @nogc bool empty() pure const nothrow
- {
- bound_t size = size();
- mixin(generateLoopCode!("if (min[@] == max[@]) return true;", N)());
- return false;
- }
- }
-
- /// Returns: true if it contains point.
- @nogc bool contains(bound_t point) pure const nothrow
- {
- assert(isSorted());
- for(int i = 0; i < N; ++i)
- if ( !(point[i] >= min[i] && point[i] < max[i]) )
- return false;
-
- return true;
- }
-
- static if (N >= 2)
- {
- /// Returns: true if it contains point `x`, `y`.
- @nogc bool contains(T x, T y) pure const nothrow
- {
- assert(isSorted());
- if ( !(x >= min.x && x < max.x) )
- return false;
- if ( !(y >= min.y && y < max.y) )
- return false;
- return true;
- }
- }
-
- static if (N >= 3)
- {
- /// Returns: true if it contains point `x`, `y`, `z`.
- @nogc bool contains(T x, T y, T z) pure const nothrow
- {
- assert(isSorted());
- if ( !(x >= min.x && x < max.x) )
- return false;
- if ( !(y >= min.y && y < max.y) )
- return false;
- if ( !(z >= min.z && z < max.z) )
- return false;
- return true;
- }
- }
-
- /// Returns: true if it contains box other.
- @nogc bool contains(Box other) pure const nothrow
- {
- assert(isSorted());
- assert(other.isSorted());
-
- mixin(generateLoopCode!("if ( (other.min[@] < min[@]) || (other.max[@] > max[@]) ) return false;", N)());
- return true;
- }
-
- /// Euclidean squared distance from a point.
- /// See_also: Numerical Recipes Third Edition (2007)
- @nogc real squaredDistance(bound_t point) pure const nothrow
- {
- assert(isSorted());
- real distanceSquared = 0;
- for (int i = 0; i < N; ++i)
- {
- if (point[i] < min[i])
- distanceSquared += (point[i] - min[i]) ^^ 2;
-
- if (point[i] > max[i])
- distanceSquared += (point[i] - max[i]) ^^ 2;
- }
- return distanceSquared;
- }
-
- /// Euclidean distance from a point.
- /// See_also: squaredDistance.
- @nogc real distance(bound_t point) pure const nothrow
- {
- return sqrt(squaredDistance(point));
- }
-
- /// Euclidean squared distance from another box.
- /// See_also: Numerical Recipes Third Edition (2007)
- @nogc real squaredDistance(Box o) pure const nothrow
- {
- assert(isSorted());
- assert(o.isSorted());
- real distanceSquared = 0;
- for (int i = 0; i < N; ++i)
- {
- if (o.max[i] < min[i])
- distanceSquared += (o.max[i] - min[i]) ^^ 2;
-
- if (o.min[i] > max[i])
- distanceSquared += (o.min[i] - max[i]) ^^ 2;
- }
- return distanceSquared;
- }
-
- /// Euclidean distance from another box.
- /// See_also: squaredDistance.
- @nogc real distance(Box o) pure const nothrow
- {
- return sqrt(squaredDistance(o));
- }
-
- /// Assumes sorted boxes.
- /// This function deals with empty boxes correctly.
- /// Returns: Intersection of two boxes.
- @nogc Box intersection(Box o) pure const nothrow
- {
- assert(isSorted());
- assert(o.isSorted());
-
- // Return an empty box if one of the boxes is empty
- if (empty())
- return this;
-
- if (o.empty())
- return o;
-
- Box result = void;
- for (int i = 0; i < N; ++i)
- {
- T maxOfMins = (min.v[i] > o.min.v[i]) ? min.v[i] : o.min.v[i];
- T minOfMaxs = (max.v[i] < o.max.v[i]) ? max.v[i] : o.max.v[i];
- result.min.v[i] = maxOfMins;
- result.max.v[i] = minOfMaxs >= maxOfMins ? minOfMaxs : maxOfMins;
- }
- return result;
- }
-
- /// Assumes sorted boxes.
- /// This function deals with empty boxes correctly.
- /// Returns: Intersection of two boxes.
- @nogc bool intersects(Box other) pure const nothrow
- {
- Box inter = this.intersection(other);
- return inter.isSorted() && !inter.empty();
- }
-
- /// Extends the area of this Box.
- @nogc Box grow(bound_t space) pure const nothrow
- {
- Box res = this;
- res.min -= space;
- res.max += space;
- return res;
- }
-
- /// Shrink the area of this Box. The box might became unsorted.
- @nogc Box shrink(bound_t space) pure const nothrow
- {
- return grow(-space);
- }
-
- /// Extends the area of this Box.
- @nogc Box grow(T space) pure const nothrow
- {
- return grow(bound_t(space));
- }
-
- /// Translate this Box.
- @nogc Box translate(bound_t offset) pure const nothrow
- {
- return Box(min + offset, max + offset);
- }
-
- static if (N >= 2)
- {
- /// Translate this Box by `x`, `y`.
- @nogc Box translate(T x, T y) pure const nothrow
- {
- Box res = this;
- res.min.x += x;
- res.min.y += y;
- res.max.x += x;
- res.max.y += y;
- return res;
- }
- }
-
- static if (N >= 3)
- {
- /// Translate this Box by `x`, `y`.
- @nogc Box translate(T x, T y, T z) pure const nothrow
- {
- Box res = this;
- res.min.x += x;
- res.min.y += y;
- res.min.z += z;
- res.max.x += x;
- res.max.y += y;
- res.max.z += z;
- return res;
- }
- }
-
- /// Shrinks the area of this Box.
- /// Returns: Shrinked box.
- @nogc Box shrink(T space) pure const nothrow
- {
- return shrink(bound_t(space));
- }
-
- /// Expands the box to include point.
- /// Returns: Expanded box.
- @nogc Box expand(bound_t point) pure const nothrow
- {
- import vector = gfm.math.vector;
- return Box(vector.minByElem(min, point), vector.maxByElem(max, point));
- }
-
- /// Expands the box to include another box.
- /// This function deals with empty boxes correctly.
- /// Returns: Expanded box.
- @nogc Box expand(Box other) pure const nothrow
- {
- assert(isSorted());
- assert(other.isSorted());
-
- // handle empty boxes
- if (empty())
- return other;
- if (other.empty())
- return this;
-
- Box result = void;
- for (int i = 0; i < N; ++i)
- {
- T minOfMins = (min.v[i] < other.min.v[i]) ? min.v[i] : other.min.v[i];
- T maxOfMaxs = (max.v[i] > other.max.v[i]) ? max.v[i] : other.max.v[i];
- result.min.v[i] = minOfMins;
- result.max.v[i] = maxOfMaxs;
- }
- return result;
- }
-
- /// Returns: true if each dimension of the box is >= 0.
- @nogc bool isSorted() pure const nothrow
- {
- for(int i = 0; i < N; ++i)
- {
- if (min[i] > max[i])
- return false;
- }
- return true;
- }
-
- /// Returns: Absolute value of the Box to ensure each dimension of the
- /// box is >= 0.
- @nogc Box abs() pure const nothrow
- {
- Box!(T, N) s = this;
- for (int i = 0; i < N; ++i)
- {
- if (s.min.v[i] > s.max.v[i])
- {
- T tmp = s.min.v[i];
- s.min.v[i] = s.max.v[i];
- s.max.v[i] = tmp;
- }
- }
- return s;
- }
-
- /// Assign with another box.
- @nogc ref Box opAssign(U)(U x) nothrow if (isBox!U)
- {
- static if(is(U.element_t : T))
- {
- static if(U._size == _size)
- {
- min = x.min;
- max = x.max;
- }
- else
- {
- static assert(false, "no conversion between boxes with different dimensions");
- }
- }
- else
- {
- static assert(false, "no conversion from " ~ U.element_t.stringof ~ " to " ~ element_t.stringof);
- }
- return this;
- }
-
- /// Returns: true if comparing equal boxes.
- @nogc bool opEquals(U)(U other) pure const nothrow if (is(U : Box))
- {
- return (min == other.min) && (max == other.max);
- }
-
- /// Cast to other box types.
- @nogc U opCast(U)() pure const nothrow if (isBox!U)
- {
- U b = void;
- for(int i = 0; i < N; ++i)
- {
- b.min[i] = cast(U.element_t)(min[i]);
- b.max[i] = cast(U.element_t)(max[i]);
- }
- return b; // return a box where each element has been casted
- }
-
- static if (N == 2)
- {
- /// Helper function to create rectangle with a given point, width and height.
- static @nogc Box rectangle(T x, T y, T width, T height) pure nothrow
- {
- return Box(x, y, x + width, y + height);
- }
- }
- }
-
- private
- {
- enum _size = N;
- alias T element_t;
- }
-}
-
-/// Instanciate to use a 2D box.
-template box2(T)
-{
- alias Box!(T, 2) box2;
-}
-
-/// Instanciate to use a 3D box.
-template box3(T)
-{
- alias Box!(T, 3) box3;
-}
-
-
-alias box2!int box2i; /// 2D box with integer coordinates.
-alias box3!int box3i; /// 3D box with integer coordinates.
-alias box2!float box2f; /// 2D box with float coordinates.
-alias box3!float box3f; /// 3D box with float coordinates.
-alias box2!double box2d; /// 2D box with double coordinates.
-alias box3!double box3d; /// 3D box with double coordinates.
-
-/// Returns: A 2D rectangle with point `x`,`y`, `width` and `height`.
-box2i rectangle(int x, int y, int width, int height) pure nothrow @nogc
-{
- return box2i(x, y, x + width, y + height);
-}
-
-/// Returns: A 2D rectangle with point `x`,`y`, `width` and `height`.
-box2f rectanglef(float x, float y, float width, float height) pure nothrow @nogc
-{
- return box2f(x, y, x + width, y + height);
-}
-
-/// Returns: A 2D rectangle with point `x`,`y`, `width` and `height`.
-box2d rectangled(double x, double y, double width, double height) pure nothrow @nogc
-{
- return box2d(x, y, x + width, y + height);
-}
-
-
-unittest
-{
- box2i a = box2i(1, 2, 3, 4);
- assert(a.width == 2);
- assert(a.height == 2);
- assert(a.volume == 4);
- box2i b = box2i(vec2i(1, 2), vec2i(3, 4));
- assert(a == b);
-
- box3i q = box3i(-3, -2, -1, 0, 1, 2);
- q.bound_t s = q.bound_t(11, 17, 19);
- q.bound_t q_min = q.min;
- assert((q.size = s) == s);
- assert(q.size == s);
- assert(q.min == q_min);
- assert(q.max == q.min + s);
- assert(q.max - q.min == s);
-
- assert((q.width = s.z) == s.z);
- assert(q.width == s.z);
- assert(q.min.x == q_min.x);
- assert(q.max.x == q.min.x + s.z);
- assert(q.max.x - q.min.x == s.z);
-
- assert((q.height = s.y) == s.y);
- assert(q.height == s.y);
- assert(q.min.y == q_min.y);
- assert(q.max.y == q.min.y + s.y);
- assert(q.max.y - q.min.y == s.y);
-
- assert((q.depth = s.x) == s.x);
- assert(q.depth == s.x);
- assert(q.min.z == q_min.z);
- assert(q.max.z == q.min.z + s.x);
- assert(q.max.z - q.min.z == s.x);
-
- assert(q.size == s.zyx);
-
- box3i n = box3i(2, 1, 0, -1, -2, -3);
- assert(n.abs == box3i(-1, -2, -3, 2, 1, 0));
-
- box2f bf = cast(box2f)b;
- assert(bf == box2f(1.0f, 2.0f, 3.0f, 4.0f));
-
- box3f qf = box3f(-0, 1f, 2.5f, 3.25f, 5.125f, 7.0625f);
- qf.bound_t sf = qf.bound_t(-11.5f, -17.25f, -19.125f);
- qf.bound_t qf_min = qf.min;
- assert((qf.size = sf) == sf);
- assert(qf.size == sf);
- assert(qf.min == qf_min);
- assert(qf.max == qf.min + sf);
- assert(qf.max - qf.min == sf);
-
- assert((qf.width = sf.z) == sf.z);
- assert(qf.width == sf.z);
- assert(qf.min.x == qf_min.x);
- assert(qf.max.x == qf.min.x + sf.z);
- assert(qf.max.x - qf.min.x == sf.z);
-
- assert((qf.height = sf.y) == sf.y);
- assert(qf.height == sf.y);
- assert(qf.min.y == qf_min.y);
- assert(qf.max.y == qf.min.y + sf.y);
- assert(qf.max.y - qf.min.y == sf.y);
-
- assert((qf.depth = sf.x) == sf.x);
- assert(qf.depth == sf.x);
- assert(qf.min.z == qf_min.z);
- assert(qf.max.z == qf.min.z + sf.x);
- assert(qf.max.z - qf.min.z == sf.x);
-
- assert(qf.size == sf.zyx);
-
- box2i c = box2i(0, 0, 1,1);
- assert(c.translate(vec2i(3, 3)) == box2i(3, 3, 4, 4));
- assert(c.translate(3, 3) == box2i(3, 3, 4, 4));
- assert(c.contains(vec2i(0, 0)));
- assert(c.contains(0, 0));
- assert(!c.contains(vec2i(1, 1)));
- assert(!c.contains(1, 1));
- assert(b.contains(b));
- box2i d = c.expand(vec2i(3, 3));
- assert(d.contains(vec2i(2, 2)));
-
- assert(d == d.expand(d));
-
- assert(!box2i(0, 0, 4, 4).contains(box2i(2, 2, 6, 6)));
-
- assert(box2f(0, 0, 0, 0).empty());
- assert(!box2f(0, 2, 1, 1).empty());
- assert(!box2f(0, 0, 1, 1).empty());
-
- assert(box2i(260, 100, 360, 200).intersection(box2i(100, 100, 200, 200)).empty());
-
- // union with empty box is identity
- assert(a.expand(box2i(10, 4, 10, 6)) == a);
-
- // intersection with empty box is empty
- assert(a.intersection(box2i(10, 4, 10, 6)).empty);
-
- assert(box2i.rectangle(1, 2, 3, 4) == box2i(1, 2, 4, 6));
- assert(rectangle(1, 2, 3, 4) == box2i(1, 2, 4, 6));
- assert(rectanglef(1, 2, 3, 4) == box2f(1, 2, 4, 6));
- assert(rectangled(1, 2, 3, 4) == box2d(1, 2, 4, 6));
-}
-
-/// True if `T` is a kind of Box
-enum isBox(T) = is(T : Box!U, U...);
-
-unittest
-{
- static assert( isBox!box2f);
- static assert( isBox!box3d);
- static assert( isBox!(Box!(real, 2)));
- static assert(!isBox!vec2f);
-}
-
-/// Get the numeric type used to measure a box's dimensions.
-alias DimensionType(T : Box!U, U...) = U[0];
-
-///
-unittest
-{
- static assert(is(DimensionType!box2f == float));
- static assert(is(DimensionType!box3d == double));
-}
-
-private
-{
- static string generateLoopCode(string formatString, int N)() pure nothrow
- {
- string result;
- for (int i = 0; i < N; ++i)
- {
- string index = ctIntToString(i);
- // replace all @ by indices
- result ~= formatString.replace("@", index);
- }
- return result;
- }
-
- // Speed-up CTFE conversions
- static string ctIntToString(int n) pure nothrow
- {
- static immutable string[16] table = ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"];
- if (n < 10)
- return table[n];
- else
- return to!string(n);
- }
-}
diff --git a/external/gfm/math/funcs.d b/external/gfm/math/funcs.d
deleted file mode 100644
index fb2afb3..0000000
--- a/external/gfm/math/funcs.d
+++ /dev/null
@@ -1,486 +0,0 @@
-/**
- Useful math functions and range-based statistic computations.
-
- If you need real statistics, consider using the $(WEB github.com/dsimcha/dstats,Dstats) library.
- */
-module gfm.math.funcs;
-
-import std.math,
- std.traits,
- std.range,
- std.math;
-
-import inteli.xmmintrin;
-
-import gfm.math.vector : Vector;
-
-/// Convert from radians to degrees.
-@nogc T degrees(T)(in T x) pure nothrow
-if (isFloatingPoint!T || (is(T : Vector!(U, n), U, int n) && isFloatingPoint!U))
-{
- static if (is(T : Vector!(U, n), U, int n))
- return x * U(180 / PI);
- else
- return x * T(180 / PI);
-}
-
-/// Convert from degrees to radians.
-@nogc T radians(T)(in T x) pure nothrow
-if (isFloatingPoint!T || (is(T : Vector!(U, n), U, int n) && isFloatingPoint!U))
-{
- static if (is(T : Vector!(U, n), U, int n))
- return x * U(PI / 180);
- else
- return x * T(PI / 180);
-}
-
-/// Linear interpolation, akin to GLSL's mix.
-@nogc S lerp(S, T)(S a, S b, T t) pure nothrow
- if (is(typeof(t * b + (1 - t) * a) : S))
-{
- return t * b + (1 - t) * a;
-}
-
-/// Clamp x in [min, max], akin to GLSL's clamp.
-@nogc T clamp(T)(T x, T min, T max) pure nothrow
-{
- if (x < min)
- return min;
- else if (x > max)
- return max;
- else
- return x;
-}
-
-/// Integer truncation.
-@nogc long ltrunc(real x) nothrow // may be pure but trunc isn't pure
-{
- return cast(long)(trunc(x));
-}
-
-/// Integer flooring.
-@nogc long lfloor(real x) nothrow // may be pure but floor isn't pure
-{
- return cast(long)(floor(x));
-}
-
-/// Returns: Fractional part of x.
-@nogc T fract(T)(real x) nothrow
-{
- return x - lfloor(x);
-}
-
-/// Safe asin: input clamped to [-1, 1]
-@nogc T safeAsin(T)(T x) pure nothrow
-{
- return asin(clamp!T(x, -1, 1));
-}
-
-/// Safe acos: input clamped to [-1, 1]
-@nogc T safeAcos(T)(T x) pure nothrow
-{
- return acos(clamp!T(x, -1, 1));
-}
-
-/// Same as GLSL step function.
-/// 0.0 is returned if x < edge, and 1.0 is returned otherwise.
-@nogc T step(T)(T edge, T x) pure nothrow
-{
- return (x < edge) ? 0 : 1;
-}
-
-/// Same as GLSL smoothstep function.
-/// See: http://en.wikipedia.org/wiki/Smoothstep
-@nogc T smoothStep(T)(T a, T b, T t) pure nothrow
-{
- if (t <= a)
- return 0;
- else if (t >= b)
- return 1;
- else
- {
- T x = (t - a) / (b - a);
- return x * x * (3 - 2 * x);
- }
-}
-
-/// Returns: true of i is a power of 2.
-@nogc bool isPowerOf2(T)(T i) pure nothrow if (isIntegral!T)
-{
- assert(i >= 0);
- return (i != 0) && ((i & (i - 1)) == 0);
-}
-
-/// Integer log2
-@nogc int ilog2(T)(T i) nothrow if (isIntegral!T)
-{
- assert(i > 0);
- assert(isPowerOf2(i));
- import core.bitop : bsr;
- return bsr(i);
-}
-
-/// Computes next power of 2.
-@nogc int nextPowerOf2(int i) pure nothrow
-{
- int v = i - 1;
- v |= v >> 1;
- v |= v >> 2;
- v |= v >> 4;
- v |= v >> 8;
- v |= v >> 16;
- v++;
- assert(isPowerOf2(v));
- return v;
-}
-
-/// Computes next power of 2.
-@nogc long nextPowerOf2(long i) pure nothrow
-{
- long v = i - 1;
- v |= v >> 1;
- v |= v >> 2;
- v |= v >> 4;
- v |= v >> 8;
- v |= v >> 16;
- v |= v >> 32;
- v++;
- assert(isPowerOf2(v));
- return v;
-}
-
-/// Computes sin(x)/x accurately.
-/// See_also: $(WEB www.plunk.org/~hatch/rightway.php)
-@nogc T sinOverX(T)(T x) pure nothrow
-{
- if (1 + x * x == 1)
- return 1;
- else
- return sin(x) / x;
-}
-
-
-/// Signed integer modulo a/b where the remainder is guaranteed to be in [0..b[,
-/// even if a is negative. Only support positive dividers.
-@nogc T moduloWrap(T)(T a, T b) pure nothrow if (isSigned!T)
-{
- assert(b > 0);
- if (a >= 0)
- a = a % b;
- else
- {
- auto rem = a % b;
- x = (rem == 0) ? 0 : (-rem + b);
- }
-
- assert(x >= 0 && x < b);
- return x;
-}
-
-unittest
-{
- assert(nextPowerOf2(13) == 16);
-}
-
-/**
- * Find the root of a linear polynomial a + b x = 0
- * Returns: Number of roots.
- */
-@nogc int solveLinear(T)(T a, T b, out T root) pure nothrow if (isFloatingPoint!T)
-{
- if (b == 0)
- {
- return 0;
- }
- else
- {
- root = -a / b;
- return 1;
- }
-}
-
-
-/**
- * Finds the root roots of a quadratic polynomial a + b x + c x^2 = 0
- * Params:
- * a = Coefficient.
- * b = Coefficient.
- * c = Coefficient.
- * outRoots = array of root results, should have room for at least 2 elements.
- * Returns: Number of roots in outRoots.
- */
-@nogc int solveQuadratic(T)(T a, T b, T c, T[] outRoots) pure nothrow if (isFloatingPoint!T)
-{
- assert(outRoots.length >= 2);
- if (c == 0)
- return solveLinear(a, b, outRoots[0]);
-
- T delta = b * b - 4 * a * c;
- if (delta < 0.0 )
- return 0;
-
- delta = sqrt(delta);
- T oneOver2a = 0.5 / a;
-
- outRoots[0] = oneOver2a * (-b - delta);
- outRoots[1] = oneOver2a * (-b + delta);
- return 2;
-}
-
-
-/**
- * Finds the roots of a cubic polynomial a + b x + c x^2 + d x^3 = 0
- * Params:
- * a = Coefficient.
- * b = Coefficient.
- * c = Coefficient.
- * d = Coefficient.
- * outRoots = array of root results, should have room for at least 2 elements.
- * Returns: Number of roots in outRoots.
- * See_also: $(WEB www.codeguru.com/forum/archive/index.php/t-265551.html)
- */
-@nogc int solveCubic(T)(T a, T b, T c, T d, T[] outRoots) pure nothrow if (isFloatingPoint!T)
-{
- assert(outRoots.length >= 3);
- if (d == 0)
- return solveQuadratic(a, b, c, outRoots);
-
- // adjust coefficients
- T a1 = c / d,
- a2 = b / d,
- a3 = a / d;
-
- T Q = (a1 * a1 - 3 * a2) / 9,
- R = (2 * a1 * a1 * a1 - 9 * a1 * a2 + 27 * a3) / 54;
-
- T Qcubed = Q * Q * Q;
- T d2 = Qcubed - R * R;
-
- if (d2 >= 0)
- {
- // 3 real roots
- if (Q < 0.0)
- return 0;
- T P = R / sqrt(Qcubed);
-
- assert(-1 <= P && P <= 1);
- T theta = acos(P);
- T sqrtQ = sqrt(Q);
-
- outRoots[0] = -2 * sqrtQ * cos(theta / 3) - a1 / 3;
- outRoots[1] = -2 * sqrtQ * cos((theta + 2 * T(PI)) / 3) - a1 / 3;
- outRoots[2] = -2 * sqrtQ * cos((theta + 4 * T(PI)) / 3) - a1 / 3;
- return 3;
- }
- else
- {
- // 1 real root
- T e = (sqrt(-d) + abs(R)) ^^ cast(T)(1.0 / 3.0);
- if (R > 0)
- e = -e;
- outRoots[0] = e + Q / e - a1 / 3.0;
- return 1;
- }
-}
-
-/**
- * Returns the roots of a quartic polynomial a + b x + c x^2 + d x^3 + e x^4 = 0
- *
- * Returns number of roots. roots slice should have room for up to 4 elements.
- * Bugs: doesn't pass unit-test!
- * See_also: $(WEB mathworld.wolfram.com/QuarticEquation.html)
- */
-@nogc int solveQuartic(T)(T a, T b, T c, T d, T e, T[] roots) pure nothrow if (isFloatingPoint!T)
-{
- assert(roots.length >= 4);
-
- if (e == 0)
- return solveCubic(a, b, c, d, roots);
-
- // Adjust coefficients
- T a0 = a / e,
- a1 = b / e,
- a2 = c / e,
- a3 = d / e;
-
- // Find a root for the following cubic equation:
- // y^3 - a2 y^2 + (a1 a3 - 4 a0) y + (4 a2 a0 - a1 ^2 - a3^2 a0) = 0
- // aka Resolvent cubic
- T b0 = 4 * a2 * a0 - a1 * a1 - a3 * a3 * a0;
- T b1 = a1 * a3 - 4 * a0;
- T b2 = -a2;
- T[3] resolventCubicRoots;
- int numRoots = solveCubic!T(b0, b1, b2, 1, resolventCubicRoots[]);
- assert(numRoots == 3);
- T y = resolventCubicRoots[0];
- if (y < resolventCubicRoots[1]) y = resolventCubicRoots[1];
- if (y < resolventCubicRoots[2]) y = resolventCubicRoots[2];
-
- // Compute R, D & E
- T R = 0.25f * a3 * a3 - a2 + y;
- if (R < 0.0)
- return 0;
- R = sqrt(R);
-
- T D = void,
- E = void;
- if (R == 0)
- {
- T d1 = 0.75f * a3 * a3 - 2 * a2;
- T d2 = 2 * sqrt(y * y - 4 * a0);
- D = sqrt(d1 + d2) * 0.5f;
- E = sqrt(d1 - d2) * 0.5f;
- }
- else
- {
- T Rsquare = R * R;
- T Rrec = 1 / R;
- T d1 = 0.75f * a3 * a3 - Rsquare - 2 * a2;
- T d2 = 0.25f * Rrec * (4 * a3 * a2 - 8 * a1 - a3 * a3 * a3);
- D = sqrt(d1 + d2) * 0.5f;
- E = sqrt(d1 - d2) * 0.5f;
- }
-
- // Compute the 4 roots
- a3 *= -0.25f;
- R *= 0.5f;
-
- roots[0] = a3 + R + D;
- roots[1] = a3 + R - D;
- roots[2] = a3 - R + E;
- roots[3] = a3 - R - E;
- return 4;
-}
-
-
-unittest
-{
- bool arrayContainsRoot(double[] arr, double root)
- {
- foreach(e; arr)
- if (abs(e - root) < 1e-7)
- return true;
- return false;
- }
-
- // test quadratic
- {
- double[3] roots;
- int numRoots = solveCubic!double(-2, -3 / 2.0, 3 / 4.0, 1 / 4.0, roots[]);
- assert(numRoots == 3);
- assert(arrayContainsRoot(roots[], -4));
- assert(arrayContainsRoot(roots[], -1));
- assert(arrayContainsRoot(roots[], 2));
- }
-
- // test quartic
- {
- double[4] roots;
- int numRoots = solveQuartic!double(0, -2, -1, 2, 1, roots[]);
-
- assert(numRoots == 4);
- assert(arrayContainsRoot(roots[], -2));
- assert(arrayContainsRoot(roots[], -1));
- assert(arrayContainsRoot(roots[], 0));
- assert(arrayContainsRoot(roots[], 1));
- }
-}
-
-/// Arithmetic mean.
-double average(R)(R r) if (isInputRange!R)
-{
- if (r.empty)
- return double.nan;
-
- typeof(r.front()) sum = 0;
- long count = 0;
- foreach(e; r)
- {
- sum += e;
- ++count;
- }
- return sum / count;
-}
-
-/// Minimum of a range.
-double minElement(R)(R r) if (isInputRange!R)
-{
- // do like Javascript for an empty range
- if (r.empty)
- return double.infinity;
-
- return minmax!("<", R)(r);
-}
-
-/// Maximum of a range.
-double maxElement(R)(R r) if (isInputRange!R)
-{
- // do like Javascript for an empty range
- if (r.empty)
- return -double.infinity;
-
- return minmax!(">", R)(r);
-}
-
-/// Variance of a range.
-double variance(R)(R r) if (isForwardRange!R)
-{
- if (r.empty)
- return double.nan;
-
- auto avg = average(r.save); // getting the average
-
- typeof(avg) sum = 0;
- long count = 0;
- foreach(e; r)
- {
- sum += (e - avg) ^^ 2;
- ++count;
- }
- if (count <= 1)
- return 0.0;
- else
- return (sum / (count - 1.0)); // using sample std deviation as estimator
-}
-
-/// Standard deviation of a range.
-double standardDeviation(R)(R r) if (isForwardRange!R)
-{
- return sqrt(variance(r));
-}
-
-private
-{
- typeof(R.front()) minmax(string op, R)(R r) if (isInputRange!R)
- {
- assert(!r.empty);
- auto best = r.front();
- r.popFront();
- foreach(e; r)
- {
- mixin("if (e " ~ op ~ " best) best = e;");
- }
- return best;
- }
-}
-
-/// SSE approximation of reciprocal square root.
-@nogc T inverseSqrt(T)(T x) pure nothrow if (isFloatingPoint!T)
-{
- static if (is(T == float))
- {
- __m128 V = _mm_set_ss(x);
- V = _mm_rsqrt_ss(V);
- return _mm_cvtss_f32(V);
- }
- else
- {
- return 1 / sqrt(x);
- }
-}
-
-unittest
-{
- assert(abs( inverseSqrt!float(1) - 1) < 1e-3 );
- assert(abs( inverseSqrt!double(1) - 1) < 1e-3 );
-}
diff --git a/external/gfm/math/matrix.d b/external/gfm/math/matrix.d
deleted file mode 100644
index d10fb3f..0000000
--- a/external/gfm/math/matrix.d
+++ /dev/null
@@ -1,948 +0,0 @@
-/// Custom sized 2-dimension Matrices
-module gfm.math.matrix;
-
-import std.math,
- std.typetuple,
- std.traits,
- std.string,
- std.typecons,
- std.conv;
-
-import gfm.math.vector,
- gfm.math.shapes,
- gfm.math.quaternion;
-
-/// Generic non-resizeable matrix with R rows and C columns.
-/// Intended for 3D use (size 3x3 and 4x4).
-/// Important: Matrices here are in row-major order whereas OpenGL is column-major.
-/// Params:
-/// T = type of elements
-/// R = number of rows
-/// C = number of columns
-struct Matrix(T, int R, int C)
-{
- public
- {
- static assert(R >= 1 && C >= 1);
-
- alias Vector!(T, C) row_t;
- alias Vector!(T, R) column_t;
-
- enum bool isSquare = (R == C);
-
- // fields definition
- union
- {
- T[C*R] v; // all elements
- row_t[R] rows; // all rows
- T[C][R] c; // components
- }
-
- @nogc this(U...)(U values) pure nothrow
- {
- static if ((U.length == C*R) && allSatisfy!(isTAssignable, U))
- {
- // construct with components
- foreach(int i, x; values)
- v[i] = x;
- }
- else static if ((U.length == 1) && (isAssignable!(U[0])) && (!is(U[0] : Matrix)))
- {
- // construct with assignment
- opAssign!(U[0])(values[0]);
- }
- else static assert(false, "cannot create a matrix from given arguments");
- }
-
- /// Construct a matrix from columns.
- @nogc static Matrix fromColumns(column_t[] columns) pure nothrow
- {
- assert(columns.length == C);
- Matrix res;
- for (int i = 0; i < R; ++i)
- for (int j = 0; j < C; ++j)
- {
- res.c[i][j] = columns[j][i];
- }
- return res;
- }
-
- /// Construct a matrix from rows.
- @nogc static Matrix fromRows(row_t[] rows) pure nothrow
- {
- assert(rows.length == R);
- Matrix res;
- res.rows[] = rows[];
- return res;
- }
-
- /// Construct matrix with a scalar.
- @nogc this(U)(T x) pure nothrow
- {
- for (int i = 0; i < _N; ++i)
- v[i] = x;
- }
-
- /// Assign with a scalar.
- @nogc ref Matrix opAssign(U : T)(U x) pure nothrow
- {
- for (int i = 0; i < R * C; ++i)
- v[i] = x;
- return this;
- }
-
- /// Assign with a samey matrice.
- @nogc ref Matrix opAssign(U : Matrix)(U x) pure nothrow
- {
- for (int i = 0; i < R * C; ++i)
- v[i] = x.v[i];
- return this;
- }
-
- /// Assign from other small matrices (same size, compatible type).
- @nogc ref Matrix opAssign(U)(U x) pure nothrow
- if (isMatrixInstantiation!U
- && is(U._T : _T)
- && (!is(U: Matrix))
- && (U._R == R) && (U._C == C))
- {
- for (int i = 0; i < R * C; ++i)
- v[i] = x.v[i];
- return this;
- }
-
- /// Assign with a static array of size R * C.
- @nogc ref Matrix opAssign(U)(U x) pure nothrow
- if ((isStaticArray!U)
- && is(typeof(x[0]) : T)
- && (U.length == R * C))
- {
- for (int i = 0; i < R * C; ++i)
- v[i] = x[i];
- return this;
- }
-
- /// Assign with a static array of shape (R, C).
- @nogc ref Matrix opAssign(U)(U x) pure nothrow
- if ((isStaticArray!U) && isStaticArray!(typeof(x[0]))
- && is(typeof(x[0][0]) : T)
- && (U.length == R)
- && (x[0].length == C))
- {
- foreach (i; 0..R)
- rows[i] = x[i];
- return this;
- }
-
- /// Assign with a dynamic array of size R * C.
- @nogc ref Matrix opAssign(U)(U x) pure nothrow
- if ((isDynamicArray!U)
- && is(typeof(x[0]) : T))
- {
- assert(x.length == R * C);
- for (int i = 0; i < R * C; ++i)
- v[i] = x[i];
- return this;
- }
-
- /// Assign with a dynamic array of shape (R, C).
- @nogc ref Matrix opAssign(U)(U x) pure nothrow
- if ((isDynamicArray!U) && isDynamicArray!(typeof(x[0]))
- && is(typeof(x[0][0]) : T))
- {
- assert(x.length == R);
- foreach (i; 0..R)
- {
- assert(x[i].length == C);
- rows[i] = x[i];
- }
- return this;
- }
-
- /// Return a pointer to content.
- @nogc inout(T)* ptr() pure inout nothrow @property
- {
- return v.ptr;
- }
-
- /// Returns a column as a vector
- /// Returns: column j as a vector.
- @nogc column_t column(int j) pure const nothrow
- {
- column_t res = void;
- for (int i = 0; i < R; ++i)
- res.v[i] = c[i][j];
- return res;
- }
-
- /// Returns a row as a vector
- /// Returns: row i as a vector.
- @nogc row_t row(int i) pure const nothrow
- {
- return rows[i];
- }
-
- /// Covnerts to pretty string.
- string toString() const nothrow
- {
- try
- return format("%s", v);
- catch (Exception e)
- assert(false); // should not happen since format is right
- }
-
- /// Matrix * scalar multiplication.
- @nogc Matrix opBinary(string op)(T factor) pure const nothrow if (op == "*")
- {
- Matrix result = void;
-
- for (int i = 0; i < R; ++i)
- {
- for (int j = 0; j < C; ++j)
- {
- result.c[i][j] = c[i][j] * factor;
- }
- }
- return result;
- }
-
- /// Matrix * vector multiplication.
- @nogc column_t opBinary(string op)(row_t x) pure const nothrow if (op == "*")
- {
- column_t res = void;
- for (int i = 0; i < R; ++i)
- {
- T sum = 0;
- for (int j = 0; j < C; ++j)
- {
- sum += c[i][j] * x.v[j];
- }
- res.v[i] = sum;
- }
- return res;
- }
-
- /// Matrix * matrix multiplication.
- @nogc auto opBinary(string op, U)(U x) pure const nothrow
- if (isMatrixInstantiation!U && (U._R == C) && (op == "*"))
- {
- Matrix!(T, R, U._C) result = void;
-
- for (int i = 0; i < R; ++i)
- {
- for (int j = 0; j < U._C; ++j)
- {
- T sum = 0;
- for (int k = 0; k < C; ++k)
- sum += c[i][k] * x.c[k][j];
- result.c[i][j] = sum;
- }
- }
- return result;
- }
-
- /// Matrix add and substraction.
- @nogc Matrix opBinary(string op, U)(U other) pure const nothrow
- if (is(U : Matrix) && (op == "+" || op == "-"))
- {
- Matrix result = void;
-
- for (int i = 0; i < R; ++i)
- {
- for (int j = 0; j < C; ++j)
- {
- mixin("result.c[i][j] = c[i][j] " ~ op ~ " other.c[i][j];");
- }
- }
- return result;
- }
-
- // matrix *= scalar
- @nogc ref Matrix opOpAssign(string op, U : T)(U x) pure nothrow if (op == "*")
- {
- for (int i = 0; i < R * C; ++i)
- v[i] *= x;
- return this;
- }
-
- /// Assignment operator with another samey matrix.
- @nogc ref Matrix opOpAssign(string op, U)(U operand) pure nothrow
- if (is(U : Matrix) && (op == "*" || op == "+" || op == "-"))
- {
- mixin("Matrix result = this " ~ op ~ " operand;");
- return opAssign!Matrix(result);
- }
-
- /// Matrix +=
- /// Matrix -=
- @nogc ref Matrix opOpAssign(string op, U)(U operand) pure nothrow
- if ((isConvertible!U) && (op == "*" || op == "+" || op == "-"))
- {
- Matrix conv = operand;
- return opOpAssign!op(conv);
- }
-
- /// Cast to other matrix types.
- /// If the size are different, the resulting matrix is truncated
- /// and/or filled with identity coefficients.
- @nogc U opCast(U)() pure const nothrow if (isMatrixInstantiation!U)
- {
- U res = U.identity();
- enum minR = R < U._R ? R : U._R;
- enum minC = C < U._C ? C : U._C;
- for (int i = 0; i < minR; ++i)
- for (int j = 0; j < minC; ++j)
- {
- res.c[i][j] = cast(U._T)(c[i][j]);
- }
- return res;
- }
-
- @nogc bool opEquals(U)(U other) pure const nothrow if (is(U : Matrix))
- {
- for (int i = 0; i < R * C; ++i)
- if (v[i] != other.v[i])
- return false;
- return true;
- }
-
- @nogc bool opEquals(U)(U other) pure const nothrow
- if ((isAssignable!U) && (!is(U: Matrix)))
- {
- Matrix conv = other;
- return opEquals(conv);
- }
-
- // +matrix, -matrix, ~matrix, !matrix
- @nogc Matrix opUnary(string op)() pure const nothrow if (op == "+" || op == "-" || op == "~" || op == "!")
- {
- Matrix res = void;
- for (int i = 0; i < N; ++i)
- mixin("res.v[i] = " ~ op ~ "v[i];");
- return res;
- }
-
- /// Convert 3x3 rotation matrix to quaternion.
- /// See_also: 3D Math Primer for Graphics and Game Development.
- @nogc U opCast(U)() pure const nothrow if (isQuaternionInstantiation!U
- && is(U._T : _T)
- && (_R == 3) && (_C == 3))
- {
- T fourXSquaredMinus1 = c[0][0] - c[1][1] - c[2][2];
- T fourYSquaredMinus1 = c[1][1] - c[0][0] - c[2][2];
- T fourZSquaredMinus1 = c[2][2] - c[0][0] - c[1][1];
- T fourWSquaredMinus1 = c[0][0] + c[1][1] + c[2][2];
-
- int biggestIndex = 0;
- T fourBiggestSquaredMinus1 = fourWSquaredMinus1;
-
- if(fourXSquaredMinus1 > fourBiggestSquaredMinus1)
- {
- fourBiggestSquaredMinus1 = fourXSquaredMinus1;
- biggestIndex = 1;
- }
-
- if(fourYSquaredMinus1 > fourBiggestSquaredMinus1)
- {
- fourBiggestSquaredMinus1 = fourYSquaredMinus1;
- biggestIndex = 2;
- }
-
- if(fourZSquaredMinus1 > fourBiggestSquaredMinus1)
- {
- fourBiggestSquaredMinus1 = fourZSquaredMinus1;
- biggestIndex = 3;
- }
-
- T biggestVal = sqrt(fourBiggestSquaredMinus1 + 1) / 2;
- T mult = 1 / (biggestVal * 4);
-
- U quat;
- switch(biggestIndex)
- {
- case 1:
- quat.w = (c[1][2] - c[2][1]) * mult;
- quat.x = biggestVal;
- quat.y = (c[0][1] + c[1][0]) * mult;
- quat.z = (c[2][0] + c[0][2]) * mult;
- break;
-
- case 2:
- quat.w = (c[2][0] - c[0][2]) * mult;
- quat.x = (c[0][1] + c[1][0]) * mult;
- quat.y = biggestVal;
- quat.z = (c[1][2] + c[2][1]) * mult;
- break;
-
- case 3:
- quat.w = (c[0][1] - c[1][0]) * mult;
- quat.x = (c[2][0] + c[0][2]) * mult;
- quat.y = (c[1][2] + c[2][1]) * mult;
- quat.z = biggestVal;
- break;
-
- default: // biggestIndex == 0
- quat.w = biggestVal;
- quat.x = (c[1][2] - c[2][1]) * mult;
- quat.y = (c[2][0] - c[0][2]) * mult;
- quat.z = (c[0][1] - c[1][0]) * mult;
- break;
- }
-
- return quat;
- }
-
- /// Converts a 4x4 rotation matrix to quaternion.
- @nogc U opCast(U)() pure const nothrow if (isQuaternionInstantiation!U
- && is(U._T : _T)
- && (_R == 4) && (_C == 4))
- {
- auto m3 = cast(mat3!T)(this);
- return cast(U)(m3);
- }
-
- static if (isSquare && isFloatingPoint!T && R == 1)
- {
- /// Returns an inverted copy of this matrix
- /// Returns: inverse of matrix.
- /// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices.
- @nogc Matrix inverse() pure const nothrow
- {
- assert(c[0][0] != 0); // Programming error if matrix is not invertible.
- return Matrix( 1 / c[0][0]);
- }
- }
-
- static if (isSquare && isFloatingPoint!T && R == 2)
- {
- /// Returns an inverted copy of this matrix
- /// Returns: inverse of matrix.
- /// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices.
- @nogc Matrix inverse() pure const nothrow
- {
- T det = (c[0][0] * c[1][1] - c[0][1] * c[1][0]);
- assert(det != 0); // Programming error if matrix is not invertible.
- T invDet = 1 / det;
- return Matrix( c[1][1] * invDet, -c[0][1] * invDet,
- -c[1][0] * invDet, c[0][0] * invDet);
- }
- }
-
- static if (isSquare && isFloatingPoint!T && R == 3)
- {
- /// Returns an inverted copy of this matrix
- /// Returns: inverse of matrix.
- /// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices.
- @nogc Matrix inverse() pure const nothrow
- {
- T det = c[0][0] * (c[1][1] * c[2][2] - c[2][1] * c[1][2])
- - c[0][1] * (c[1][0] * c[2][2] - c[1][2] * c[2][0])
- + c[0][2] * (c[1][0] * c[2][1] - c[1][1] * c[2][0]);
- assert(det != 0); // Programming error if matrix is not invertible.
- T invDet = 1 / det;
-
- Matrix res = void;
- res.c[0][0] = (c[1][1] * c[2][2] - c[2][1] * c[1][2]) * invDet;
- res.c[0][1] = -(c[0][1] * c[2][2] - c[0][2] * c[2][1]) * invDet;
- res.c[0][2] = (c[0][1] * c[1][2] - c[0][2] * c[1][1]) * invDet;
- res.c[1][0] = -(c[1][0] * c[2][2] - c[1][2] * c[2][0]) * invDet;
- res.c[1][1] = (c[0][0] * c[2][2] - c[0][2] * c[2][0]) * invDet;
- res.c[1][2] = -(c[0][0] * c[1][2] - c[1][0] * c[0][2]) * invDet;
- res.c[2][0] = (c[1][0] * c[2][1] - c[2][0] * c[1][1]) * invDet;
- res.c[2][1] = -(c[0][0] * c[2][1] - c[2][0] * c[0][1]) * invDet;
- res.c[2][2] = (c[0][0] * c[1][1] - c[1][0] * c[0][1]) * invDet;
- return res;
- }
- }
-
- static if (isSquare && isFloatingPoint!T && R == 4)
- {
- /// Returns an inverted copy of this matrix
- /// Returns: inverse of matrix.
- /// Note: Matrix inversion is provided for 1x1, 2x2, 3x3 and 4x4 floating point matrices.
- @nogc Matrix inverse() pure const nothrow
- {
- T det2_01_01 = c[0][0] * c[1][1] - c[0][1] * c[1][0];
- T det2_01_02 = c[0][0] * c[1][2] - c[0][2] * c[1][0];
- T det2_01_03 = c[0][0] * c[1][3] - c[0][3] * c[1][0];
- T det2_01_12 = c[0][1] * c[1][2] - c[0][2] * c[1][1];
- T det2_01_13 = c[0][1] * c[1][3] - c[0][3] * c[1][1];
- T det2_01_23 = c[0][2] * c[1][3] - c[0][3] * c[1][2];
-
- T det3_201_012 = c[2][0] * det2_01_12 - c[2][1] * det2_01_02 + c[2][2] * det2_01_01;
- T det3_201_013 = c[2][0] * det2_01_13 - c[2][1] * det2_01_03 + c[2][3] * det2_01_01;
- T det3_201_023 = c[2][0] * det2_01_23 - c[2][2] * det2_01_03 + c[2][3] * det2_01_02;
- T det3_201_123 = c[2][1] * det2_01_23 - c[2][2] * det2_01_13 + c[2][3] * det2_01_12;
-
- T det = - det3_201_123 * c[3][0] + det3_201_023 * c[3][1] - det3_201_013 * c[3][2] + det3_201_012 * c[3][3];
- assert(det != 0); // Programming error if matrix is not invertible.
- T invDet = 1 / det;
-
- T det2_03_01 = c[0][0] * c[3][1] - c[0][1] * c[3][0];
- T det2_03_02 = c[0][0] * c[3][2] - c[0][2] * c[3][0];
- T det2_03_03 = c[0][0] * c[3][3] - c[0][3] * c[3][0];
- T det2_03_12 = c[0][1] * c[3][2] - c[0][2] * c[3][1];
- T det2_03_13 = c[0][1] * c[3][3] - c[0][3] * c[3][1];
- T det2_03_23 = c[0][2] * c[3][3] - c[0][3] * c[3][2];
- T det2_13_01 = c[1][0] * c[3][1] - c[1][1] * c[3][0];
- T det2_13_02 = c[1][0] * c[3][2] - c[1][2] * c[3][0];
- T det2_13_03 = c[1][0] * c[3][3] - c[1][3] * c[3][0];
- T det2_13_12 = c[1][1] * c[3][2] - c[1][2] * c[3][1];
- T det2_13_13 = c[1][1] * c[3][3] - c[1][3] * c[3][1];
- T det2_13_23 = c[1][2] * c[3][3] - c[1][3] * c[3][2];
-
- T det3_203_012 = c[2][0] * det2_03_12 - c[2][1] * det2_03_02 + c[2][2] * det2_03_01;
- T det3_203_013 = c[2][0] * det2_03_13 - c[2][1] * det2_03_03 + c[2][3] * det2_03_01;
- T det3_203_023 = c[2][0] * det2_03_23 - c[2][2] * det2_03_03 + c[2][3] * det2_03_02;
- T det3_203_123 = c[2][1] * det2_03_23 - c[2][2] * det2_03_13 + c[2][3] * det2_03_12;
-
- T det3_213_012 = c[2][0] * det2_13_12 - c[2][1] * det2_13_02 + c[2][2] * det2_13_01;
- T det3_213_013 = c[2][0] * det2_13_13 - c[2][1] * det2_13_03 + c[2][3] * det2_13_01;
- T det3_213_023 = c[2][0] * det2_13_23 - c[2][2] * det2_13_03 + c[2][3] * det2_13_02;
- T det3_213_123 = c[2][1] * det2_13_23 - c[2][2] * det2_13_13 + c[2][3] * det2_13_12;
-
- T det3_301_012 = c[3][0] * det2_01_12 - c[3][1] * det2_01_02 + c[3][2] * det2_01_01;
- T det3_301_013 = c[3][0] * det2_01_13 - c[3][1] * det2_01_03 + c[3][3] * det2_01_01;
- T det3_301_023 = c[3][0] * det2_01_23 - c[3][2] * det2_01_03 + c[3][3] * det2_01_02;
- T det3_301_123 = c[3][1] * det2_01_23 - c[3][2] * det2_01_13 + c[3][3] * det2_01_12;
-
- Matrix res = void;
- res.c[0][0] = - det3_213_123 * invDet;
- res.c[1][0] = + det3_213_023 * invDet;
- res.c[2][0] = - det3_213_013 * invDet;
- res.c[3][0] = + det3_213_012 * invDet;
-
- res.c[0][1] = + det3_203_123 * invDet;
- res.c[1][1] = - det3_203_023 * invDet;
- res.c[2][1] = + det3_203_013 * invDet;
- res.c[3][1] = - det3_203_012 * invDet;
-
- res.c[0][2] = + det3_301_123 * invDet;
- res.c[1][2] = - det3_301_023 * invDet;
- res.c[2][2] = + det3_301_013 * invDet;
- res.c[3][2] = - det3_301_012 * invDet;
-
- res.c[0][3] = - det3_201_123 * invDet;
- res.c[1][3] = + det3_201_023 * invDet;
- res.c[2][3] = - det3_201_013 * invDet;
- res.c[3][3] = + det3_201_012 * invDet;
- return res;
- }
- }
-
- /// Returns a transposed copy of this matrix
- /// Returns: transposed matrice.
- @nogc Matrix!(T, C, R) transposed() pure const nothrow
- {
- Matrix!(T, C, R) res;
- for (int i = 0; i < C; ++i)
- for (int j = 0; j < R; ++j)
- res.c[i][j] = c[j][i];
- return res;
- }
-
- static if (isSquare && R > 1)
- {
- /// Makes a diagonal matrix from a vector.
- @nogc static Matrix diag(Vector!(T, R) v) pure nothrow
- {
- Matrix res = void;
- for (int i = 0; i < R; ++i)
- for (int j = 0; j < C; ++j)
- res.c[i][j] = (i == j) ? v.v[i] : 0;
- return res;
- }
-
- /// In-place translate by (v, 1)
- @nogc void translate(Vector!(T, R-1) v) pure nothrow
- {
- for (int i = 0; i < R; ++i)
- {
- T dot = 0;
- for (int j = 0; j + 1 < C; ++j)
- dot += v.v[j] * c[i][j];
-
- c[i][C-1] += dot;
- }
- }
-
- /// Make a translation matrix.
- @nogc static Matrix translation(Vector!(T, R-1) v) pure nothrow
- {
- Matrix res = identity();
- for (int i = 0; i + 1 < R; ++i)
- res.c[i][C-1] += v.v[i];
- return res;
- }
-
- /// In-place matrix scaling.
- void scale(Vector!(T, R-1) v) pure nothrow
- {
- for (int i = 0; i < R; ++i)
- for (int j = 0; j + 1 < C; ++j)
- c[i][j] *= v.v[j];
- }
-
- /// Make a scaling matrix.
- @nogc static Matrix scaling(Vector!(T, R-1) v) pure nothrow
- {
- Matrix res = identity();
- for (int i = 0; i + 1 < R; ++i)
- res.c[i][i] = v.v[i];
- return res;
- }
- }
-
- // rotations are implemented for 3x3 and 4x4 matrices.
- static if (isSquare && (R == 3 || R == 4) && isFloatingPoint!T)
- {
- @nogc public static Matrix rotateAxis(int i, int j)(T angle) pure nothrow
- {
- Matrix res = identity();
- const T cosa = cos(angle);
- const T sina = sin(angle);
- res.c[i][i] = cosa;
- res.c[i][j] = -sina;
- res.c[j][i] = sina;
- res.c[j][j] = cosa;
- return res;
- }
-
- /// Rotate along X axis
- /// Returns: rotation matrix along axis X
- alias rotateAxis!(1, 2) rotateX;
-
- /// Rotate along Y axis
- /// Returns: rotation matrix along axis Y
- alias rotateAxis!(2, 0) rotateY;
-
- /// Rotate along Z axis
- /// Returns: rotation matrix along axis Z
- alias rotateAxis!(0, 1) rotateZ;
-
- /// Similar to the glRotate matrix, however the angle is expressed in radians
- /// See_also: $(LINK http://www.cs.rutgers.edu/~decarlo/428/gl_man/rotate.html)
- @nogc static Matrix rotation(T angle, vec3!T axis) pure nothrow
- {
- Matrix res = identity();
- const T c = cos(angle);
- const oneMinusC = 1 - c;
- const T s = sin(angle);
- axis = axis.normalized();
- T x = axis.x,
- y = axis.y,
- z = axis.z;
- T xy = x * y,
- yz = y * z,
- xz = x * z;
-
- res.c[0][0] = x * x * oneMinusC + c;
- res.c[0][1] = x * y * oneMinusC - z * s;
- res.c[0][2] = x * z * oneMinusC + y * s;
- res.c[1][0] = y * x * oneMinusC + z * s;
- res.c[1][1] = y * y * oneMinusC + c;
- res.c[1][2] = y * z * oneMinusC - x * s;
- res.c[2][0] = z * x * oneMinusC - y * s;
- res.c[2][1] = z * y * oneMinusC + x * s;
- res.c[2][2] = z * z * oneMinusC + c;
- return res;
- }
- }
-
- // 4x4 specific transformations for 3D usage
- static if (isSquare && R == 4 && isFloatingPoint!T)
- {
- /// Orthographic projection
- /// Returns: orthographic projection.
- @nogc static Matrix orthographic(T left, T right, T bottom, T top, T near, T far) pure nothrow
- {
- T dx = right - left,
- dy = top - bottom,
- dz = far - near;
-
- T tx = -(right + left) / dx;
- T ty = -(top + bottom) / dy;
- T tz = -(far + near) / dz;
-
- return Matrix(2 / dx, 0, 0, tx,
- 0, 2 / dy, 0, ty,
- 0, 0, -2 / dz, tz,
- 0, 0, 0, 1);
- }
-
- /// Perspective projection
- /// Returns: perspective projection.
- @nogc static Matrix perspective(T FOVInRadians, T aspect, T zNear, T zFar) pure nothrow
- {
- T f = 1 / tan(FOVInRadians / 2);
- T d = 1 / (zNear - zFar);
-
- return Matrix(f / aspect, 0, 0, 0,
- 0, f, 0, 0,
- 0, 0, (zFar + zNear) * d, 2 * d * zFar * zNear,
- 0, 0, -1, 0);
- }
-
- /// Look At projection
- /// Returns: "lookAt" projection.
- /// Thanks to vuaru for corrections.
- @nogc static Matrix lookAt(vec3!T eye, vec3!T target, vec3!T up) pure nothrow
- {
- vec3!T Z = (eye - target).normalized();
- vec3!T X = cross(-up, Z).normalized();
- vec3!T Y = cross(Z, -X);
-
- return Matrix(-X.x, -X.y, -X.z, dot(X, eye),
- Y.x, Y.y, Y.z, -dot(Y, eye),
- Z.x, Z.y, Z.z, -dot(Z, eye),
- 0, 0, 0, 1);
- }
-
- /// Extract frustum from a 4x4 matrice.
- @nogc Frustum!T frustum() pure const nothrow
- {
- auto left = Plane!T(row(3) + row(0));
- auto right = Plane!T(row(3) - row(0));
- auto top = Plane!T(row(3) - row(1));
- auto bottom = Plane!T(row(3) + row(1));
- auto near = Plane!T(row(3) + row(2));
- auto far = Plane!T(row(3) - row(2));
- return Frustum!T(left, right, top, bottom, near, far);
- }
-
- }
- }
-
- package
- {
- alias T _T;
- enum _R = R;
- enum _C = C;
- }
-
- private
- {
- template isAssignable(T)
- {
- enum bool isAssignable = std.traits.isAssignable!(Matrix, T);
- }
-
- template isConvertible(T)
- {
- enum bool isConvertible = (!is(T : Matrix)) && isAssignable!T;
- }
-
- template isTAssignable(U)
- {
- enum bool isTAssignable = std.traits.isAssignable!(T, U);
- }
-
- template isRowConvertible(U)
- {
- enum bool isRowConvertible = is(U : row_t);
- }
-
- template isColumnConvertible(U)
- {
- enum bool isColumnConvertible = is(U : column_t);
- }
- }
-
- public
- {
- /// Construct an identity matrix
- /// Returns: an identity matrix.
- /// Note: the identity matrix, while only meaningful for square matrices,
- /// is also defined for non-square ones.
- @nogc static Matrix identity() pure nothrow
- {
- Matrix res = void;
- for (int i = 0; i < R; ++i)
- for (int j = 0; j < C; ++j)
- res.c[i][j] = (i == j) ? 1 : 0;
- return res;
- }
-
- /// Construct an constant matrix
- /// Returns: a constant matrice.
- @nogc static Matrix constant(U)(U x) pure nothrow
- {
- Matrix res = void;
-
- for (int i = 0; i < R * C; ++i)
- res.v[i] = cast(T)x;
- return res;
- }
- }
-}
-
-template isMatrixInstantiation(U)
-{
- private static void isMatrix(T, int R, int C)(Matrix!(T, R, C) x)
- {
- }
-
- enum bool isMatrixInstantiation = is(typeof(isMatrix(U.init)));
-}
-
-// GLSL is a big inspiration here
-// we defines types with more or less the same names
-
-///
-template mat2x2(T) { alias Matrix!(T, 2, 2) mat2x2; }
-///
-template mat3x3(T) { alias Matrix!(T, 3, 3) mat3x3; }
-///
-template mat4x4(T) { alias Matrix!(T, 4, 4) mat4x4; }
-
-// WARNING: in GLSL, first number is _columns_, second is rows
-// It is the opposite here: first number is rows, second is columns
-// With this convention mat2x3 * mat3x4 -> mat2x4.
-
-///
-template mat2x3(T) { alias Matrix!(T, 2, 3) mat2x3; }
-///
-template mat2x4(T) { alias Matrix!(T, 2, 4) mat2x4; }
-///
-template mat3x2(T) { alias Matrix!(T, 3, 2) mat3x2; }
-///
-template mat3x4(T) { alias Matrix!(T, 3, 4) mat3x4; }
-///
-template mat4x2(T) { alias Matrix!(T, 4, 2) mat4x2; }
-///
-template mat4x3(T) { alias Matrix!(T, 4, 3) mat4x3; }
-
-// shorter names for most common matrices
-alias mat2x2 mat2;///
-alias mat3x3 mat3;///
-alias mat4x4 mat4;///
-
-// Define a lot of type names
-// Most useful are probably mat4f and mat4d
-
-alias mat2!float mat2f;///
-alias mat2!double mat2d;///
-
-alias mat3!float mat3f;///
-alias mat3!double mat3d;///
-
-alias mat4!float mat4f;///
-alias mat4!double mat4d;///
-
-alias mat2x2!float mat2x2f;///
-alias mat2x2!double mat2x2d;///
-
-alias mat3x3!float mat3x3f;///
-alias mat3x3!double mat3x3d;///
-
-alias mat4x4!float mat4x4f;///
-alias mat4x4!double mat4x4d;///
-
-unittest
-{
- alias mat2i = mat2!int;
- alias mat2x3f = mat2x3!float;
- alias mat3x4f = mat3x4!float;
- alias mat2x4f = mat2x4!float;
-
- mat2i x = mat2i(0, 1,
- 2, 3);
- assert(x.c[0][0] == 0 && x.c[0][1] == 1 && x.c[1][0] == 2 && x.c[1][1] == 3);
-
- vec2i[2] cols = [vec2i(0, 2), vec2i(1, 3)];
- mat2i y = mat2i.fromColumns(cols[]);
- assert(y.c[0][0] == 0 && y.c[0][1] == 1 && y.c[1][0] == 2 && y.c[1][1] == 3);
- y = mat2i.fromRows(cols[]);
- assert(y.c[0][0] == 0 && y.c[1][0] == 1 && y.c[0][1] == 2 && y.c[1][1] == 3);
- y = y.transposed();
-
- assert(x == y);
- x = [0, 1, 2, 3];
- assert(x == y);
-
- mat2i z = x * y;
- assert(z == mat2i([2, 3, 6, 11]));
- vec2i vz = z * vec2i(2, -1);
- assert(vz == vec2i(1, 1));
-
- mat2f a = z;
- mat2d ad = a;
- ad += a;
- mat2f w = [4, 5, 6, 7];
- z = cast(mat2i)w;
- assert(w == z);
-
- {
- mat2x3f A;
- mat3x4f B;
- mat2x4f C = A * B;
- }
-
- assert(mat2i.diag(vec2i(1, 2)) == mat2i(1, 0,
- 0, 2));
-
- // Construct with a single scalar
- auto D = mat4f(1.0f);
- assert(D.v[] == [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ]);
-
- {
- double[4][3] starray = [
- [ 0, 1, 2, 3],
- [ 4, 5, 6, 7,],
- [ 8, 9, 10, 11,],
- ];
-
- // starray has the shape 3x4
- assert(starray.length == 3);
- assert(starray[0].length == 4);
-
- auto m = mat3x4!double(starray);
- assert(m.v[] == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ]);
- }
-
- {
- auto dynarray = [
- [ 0, 1, 2, 3],
- [ 4, 5, 6, 7,],
- [ 8, 9, 10, 11,],
- ];
-
- // dynarray has the shape 3x4
- assert(dynarray.length == 3);
- assert(dynarray[0].length == 4);
-
- auto m = mat3x4!double(dynarray);
- assert(m.v[] == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ]);
- }
-}
-
-// Issue #206 (matrix *= scalar) not yielding matrix * scalar but matrix * matrix(scalar)
-unittest
-{
- mat4f mvp = mat4f.identity;
- mvp *= 2;
- assert(mvp == mat4f(2, 0, 0, 0,
- 0, 2, 0, 0,
- 0, 0, 2, 0,
- 0, 0, 0, 2));
-
- mvp = mat4f.identity * 2;
- assert(mvp == mat4f(2, 0, 0, 0,
- 0, 2, 0, 0,
- 0, 0, 2, 0,
- 0, 0, 0, 2));
-
-
- mvp = mat4f(1) * mat4f(1);
- assert(mvp == mat4f(4, 4, 4, 4,
- 4, 4, 4, 4,
- 4, 4, 4, 4,
- 4, 4, 4, 4));
-
- mvp = mat4f(1);
- mvp *= mat4f(1);
- assert(mvp == mat4f(4, 4, 4, 4,
- 4, 4, 4, 4,
- 4, 4, 4, 4,
- 4, 4, 4, 4));
-}
diff --git a/external/gfm/math/package.d b/external/gfm/math/package.d
deleted file mode 100644
index 14b46b7..0000000
--- a/external/gfm/math/package.d
+++ /dev/null
@@ -1,9 +0,0 @@
-module gfm.math;
-
-public import gfm.math.funcs,
- gfm.math.vector,
- gfm.math.box,
- gfm.math.matrix,
- gfm.math.quaternion,
- gfm.math.shapes,
- gfm.math.simplerng;
diff --git a/external/gfm/math/quaternion.d b/external/gfm/math/quaternion.d
deleted file mode 100644
index 5f5c944..0000000
--- a/external/gfm/math/quaternion.d
+++ /dev/null
@@ -1,339 +0,0 @@
-///
-module gfm.math.quaternion;
-
-import std.math,
- std.string;
-
-import gfm.math.vector,
- gfm.math.matrix,
- funcs = gfm.math.funcs;
-
-/// Quaternion implementation.
-/// Holds a rotation + angle in a proper but wild space.
-struct Quaternion(T)
-{
- public
- {
- union
- {
- Vector!(T, 4u) v;
- struct
- {
- T x, y, z, w;
- }
- }
-
- /// Construct a Quaternion from a value.
- @nogc this(U)(U x) pure nothrow if (isAssignable!U)
- {
- opAssign!U(x);
- }
-
- /// Constructs a Quaternion from coordinates.
- /// Warning: order of coordinates is different from storage.
- @nogc this(T qw, T qx, T qy, T qz) pure nothrow
- {
- x = qx;
- y = qy;
- z = qz;
- w = qw;
- }
-
- /// Constructs a Quaternion from axis + angle.
- @nogc static Quaternion fromAxis(Vector!(T, 3) axis, T angle) pure nothrow
- {
- Quaternion q = void;
- axis.normalize();
- T cos_a = cos(angle / 2);
- T sin_a = sin(angle / 2);
- q.x = sin_a * axis.x;
- q.y = sin_a * axis.y;
- q.z = sin_a * axis.z;
- q.w = cos_a;
- return q; // should be normalized
- }
-
- /// Constructs a Quaternion from Euler angles.
- /// All paramers given in radians.
- /// Roll->X axis, Pitch->Y axis, Yaw->Z axis
- /// See_also: $(LINK https://www.cs.princeton.edu/~gewang/projects/darth/stuff/quat_faq.html)
- @nogc static Quaternion fromEulerAngles(T roll, T pitch, T yaw) pure nothrow
- {
- Quaternion q = void;
- T sinPitch = sin(pitch / 2);
- T cosPitch = cos(pitch / 2);
- T sinYaw = sin(yaw / 2);
- T cosYaw = cos(yaw / 2);
- T sinRoll = sin(roll / 2);
- T cosRoll = cos(roll / 2);
- T cosPitchCosYaw = cosPitch * cosYaw;
- T sinPitchSinYaw = sinPitch * sinYaw;
- q.x = sinRoll * cosPitchCosYaw - cosRoll * sinPitchSinYaw;
- q.y = cosRoll * sinPitch * cosYaw + sinRoll * cosPitch * sinYaw;
- q.z = cosRoll * cosPitch * sinYaw - sinRoll * sinPitch * cosYaw;
- q.w = cosRoll * cosPitchCosYaw + sinRoll * sinPitchSinYaw;
- return q;
- }
-
- /// Converts a quaternion to Euler angles.
- /// TODO: adds a EulerAngles type.
- /// Returns: A vector which contains roll-pitch-yaw as x-y-z.
- @nogc vec3!T toEulerAngles() pure const nothrow
- {
- mat3!T m = cast(mat3!T)(this);
-
- T pitch, yaw, roll;
- T s = sqrt(m.c[0][0] * m.c[0][0] + m.c[1][0] * m.c[1][0]);
- if (s > T.epsilon)
- {
- pitch = - atan2(m.c[2][0], s);
- yaw = atan2(m.c[1][0], m.c[0][0]);
- roll = atan2(m.c[2][1], m.c[2][2]);
- }
- else
- {
- pitch = m.c[2][0] < 0.0f ? T(PI) /2 : -T(PI) / 2;
- yaw = -atan2(m.c[0][1], m.c[1][1]);
- roll = 0.0f;
- }
- return vec3!T(roll, pitch, yaw);
- }
-
- /// Assign from another Quaternion.
- @nogc ref Quaternion opAssign(U)(U u) pure nothrow if (isQuaternionInstantiation!U && is(U._T : T))
- {
- v = u.v;
- return this;
- }
-
- /// Assign from a vector of 4 elements.
- @nogc ref Quaternion opAssign(U)(U u) pure nothrow if (is(U : Vector!(T, 4u)))
- {
- v = u;
- return this;
- }
-
- /// Converts to a pretty string.
- string toString() const nothrow
- {
- try
- return format("%s", v);
- catch (Exception e)
- assert(false); // should not happen since format is right
- }
-
- /// Normalizes a quaternion.
- @nogc void normalize() pure nothrow
- {
- v.normalize();
- }
-
- /// Returns: Normalized quaternion.
- @nogc Quaternion normalized() pure const nothrow
- {
- Quaternion res = void;
- res.v = v.normalized();
- return res;
- }
-
- /// Inverses a quaternion in-place.
- @nogc void inverse() pure nothrow
- {
- x = -x;
- y = -y;
- z = -z;
- }
-
- /// Returns: Inverse of quaternion.
- @nogc Quaternion inversed() pure const nothrow
- {
- Quaternion res = void;
- res.v = v;
- res.inverse();
- return res;
- }
-
- @nogc ref Quaternion opOpAssign(string op, U)(U q) pure nothrow
- if (is(U : Quaternion) && (op == "*"))
- {
- T nx = w * q.x + x * q.w + y * q.z - z * q.y,
- ny = w * q.y + y * q.w + z * q.x - x * q.z,
- nz = w * q.z + z * q.w + x * q.y - y * q.x,
- nw = w * q.w - x * q.x - y * q.y - z * q.z;
- x = nx;
- y = ny;
- z = nz;
- w = nw;
- return this;
- }
-
- @nogc ref Quaternion opOpAssign(string op, U)(U operand) pure nothrow if (isConvertible!U)
- {
- Quaternion conv = operand;
- return opOpAssign!op(conv);
- }
-
- @nogc Quaternion opBinary(string op, U)(U operand) pure const nothrow
- if (is(U: Quaternion) || (isConvertible!U))
- {
- Quaternion temp = this;
- return temp.opOpAssign!op(operand);
- }
-
- /// Compare two Quaternions.
- bool opEquals(U)(U other) pure const if (is(U : Quaternion))
- {
- return v == other.v;
- }
-
- /// Compare Quaternion and other types.
- bool opEquals(U)(U other) pure const nothrow if (isConvertible!U)
- {
- Quaternion conv = other;
- return opEquals(conv);
- }
-
- /// Convert to a 3x3 rotation matrix.
- /// TODO: check out why we can't do is(Unqual!U == mat3!T)
- @nogc U opCast(U)() pure const nothrow if (isMatrixInstantiation!U
- && is(U._T : _T)
- && (U._R == 3) && (U._C == 3))
- {
- // do not assume that quaternion is normalized
- T norm = x*x + y*y + z*z + w*w;
- T s = (norm > 0) ? 2 / norm : 0;
-
- T xx = x * x * s, xy = x * y * s, xz = x * z * s, xw = x * w * s,
- yy = y * y * s, yz = y * z * s, yw = y * w * s,
- zz = z * z * s, zw = z * w * s;
- return mat3!(U._T)
- (
- 1 - (yy + zz) , (xy - zw) , (xz + yw) ,
- (xy + zw) , 1 - (xx + zz) , (yz - xw) ,
- (xz - yw) , (yz + xw) , 1 - (xx + yy)
- );
- }
-
- /// Converts a to a 4x4 rotation matrix.
- /// Bugs: check why we can't do is(Unqual!U == mat4!T)
- @nogc U opCast(U)() pure const nothrow if (isMatrixInstantiation!U
- && is(U._T : _T)
- && (U._R == 4) && (U._C == 4))
- {
- auto m3 = cast(mat3!(U._T))(this);
- return cast(U)(m3);
- }
-
- /// Workaround Vector not being constructable through CTFE
- @nogc static Quaternion identity() pure nothrow @property
- {
- Quaternion q;
- q.x = q.y = q.z = 0;
- q.w = 1;
- return q;
- }
- }
-
- private
- {
- alias T _T;
-
- template isAssignable(T)
- {
- enum bool isAssignable =
- is(typeof(
- {
- T x;
- Quaternion v = x;
- }()));
- }
-
- // define types that can be converted to Quaternion, but are not Quaternion
- template isConvertible(T)
- {
- enum bool isConvertible = (!is(T : Quaternion)) && isAssignable!T;
- }
- }
-}
-
-template isQuaternionInstantiation(U)
-{
- private static void isQuaternion(T)(Quaternion!T x)
- {
- }
-
- enum bool isQuaternionInstantiation = is(typeof(isQuaternion(U.init)));
-}
-
-alias Quaternion!float quatf;///
-alias Quaternion!double quatd;///
-
-/// Linear interpolation, for quaternions.
-@nogc Quaternion!T lerp(T)(Quaternion!T a, Quaternion!T b, float t) pure nothrow
-{
- Quaternion!T res = void;
- res.v = funcs.lerp(a.v, b.v, t);
- return res;
-}
-
-/// Nlerp of quaternions
-/// Returns: Nlerp of quaternions.
-/// See_also: $(WEB keithmaggio.wordpress.com/2011/02/15/math-magician-lerp-slerp-and-nlerp/, Math Magician – Lerp, Slerp, and Nlerp)
-@nogc Quaternion!T Nlerp(T)(Quaternion!T a, Quaternion!T b, float t) pure nothrow
-{
- assert(t >= 0 && t <= 1); // else probably doesn't make sense
- Quaternion!T res = void;
- res.v = funcs.lerp(a.v, b.v, t);
- res.v.normalize();
- return res;
-}
-
-/// Slerp of quaternions
-/// Returns: Slerp of quaternions. Slerp is more expensive than Nlerp.
-/// See_also: "Understanding Slerp, Then Not Using It"
-@nogc Quaternion!T slerp(T)(Quaternion!T a, Quaternion!T b, T t) pure nothrow
-{
- assert(t >= 0 && t <= 1); // else probably doesn't make sense
-
- Quaternion!T res = void;
-
- T dotProduct = dot(a.v, b.v);
-
- // spherical lerp always has 2 potential paths
- // here we always take the shortest
- if (dotProduct < 0)
- {
- b.v *= -1;
- dotProduct = dot(a.v, b.v);
- }
-
- immutable T threshold = 10 * T.epsilon; // idMath uses 1e-6f for 32-bits float precision
- if ((1 - dotProduct) > threshold) // if small difference, use lerp
- return lerp(a, b, t);
-
- T theta_0 = funcs.safeAcos(dotProduct); // angle between this and other
- T theta = theta_0 * t; // angle between this and result
-
- vec3!T v2 = dot(b.v, a.v * dotProduct);
- v2.normalize();
-
- res.v = dot(b.v, a.v * dotProduct);
- res.v.normalize();
-
- res.v = a.v * cos(theta) + res.v * sin(theta);
- return res;
-}
-
-unittest
-{
- quatf a = quatf.fromAxis(vec3f(1, 0, 0), 1);
- quatf b = quatf.fromAxis(vec3f(0, 1, 0), 0);
- a = a * b;
-
- quatf c = lerp(a, b, 0.5f);
- quatf d = Nlerp(a, b, 0.1f);
- quatf e = slerp(a, b, 0.0f);
- quatd f = quatd(1.0, 4, 5.0, 6.0);
- quatf g = quatf.fromEulerAngles(-0.1f, 1.2f, -0.3f);
- vec3f ga = g.toEulerAngles();
-}
diff --git a/external/gfm/math/shapes.d b/external/gfm/math/shapes.d
deleted file mode 100644
index d29c246..0000000
--- a/external/gfm/math/shapes.d
+++ /dev/null
@@ -1,645 +0,0 @@
-/**
- This module implements some abstract geometric shapes:
- $(UL
- $(LI Line segments.)
- $(LI Triangle.)
- $(LI Circles/spheres.)
- $(LI Rays)
- $(LI Planes)
- $(LI Frustum)
- )
- */
-module gfm.math.shapes;
-
-import std.math,
- std.traits;
-
-import gfm.math.vector,
- gfm.math.box;
-
-/// A Segment is 2 points.
-/// When considered like a vector, it represents the arrow from a to b.
-struct Segment(T, int N)
-{
- public
- {
- alias Vector!(T, N) point_t;
- point_t a, b;
-
- static if (N == 3 && isFloatingPoint!T)
- {
- /// Segment vs plane intersection.
- /// See_also: "Geometry for Computer Graphics: Formulae, Examples
- /// and Proofs", Vince (2005), p. 62
- /// Returns: $(D true) if the segment crosses the plane exactly
- /// once, or $(D false) if the segment is parallel to or does
- /// not reach the plane
- /// Params:
- /// plane = The plane to intersect.
- /// intersection = Point of intersection.
- /// progress = Set to the point's progress between endpoints
- /// $(D a) at 0.0 and $(D b) at 1.0, or T.infinity
- /// if parallel to the plane
- @nogc bool intersect(Plane!T plane, out point_t intersection, out T progress) pure const nothrow
- {
- // direction vector from a to b
- point_t dir = b-a;
-
- // dot product will be 0 if angle to plane is 0
- T dp = dot(plane.n, dir);
- if (abs(dp) < T.epsilon)
- {
- progress = T.infinity;
- return false; // parallel to plane
- }
-
- // relative distance along segment where intersection happens
- progress = -(dot(plane.n, a) - plane.d) / dp;
-
- intersection = progress*dir + a;
- return progress >= 0 && progress <= 1;
- }
- }
- }
-}
-
-alias Segment!(float, 2) seg2f; /// 2D float segment.
-alias Segment!(float, 3) seg3f; /// 3D float segment.
-alias Segment!(double, 2) seg2d; /// 2D double segment.
-alias Segment!(double, 3) seg3d; /// 3D double segment.
-alias Segment!(int, 2) seg2i; /// 2D integer segment.
-alias Segment!(int, 3) seg3i; /// 3D integer segment.
-
-/// A Triangle is 3 points.
-struct Triangle(T, int N)
-{
- public
- {
- alias Vector!(T, N) point_t;
- point_t a, b, c;
-
- static if (N == 2)
- {
- /// Returns: Area of a 2D triangle.
- @nogc T area() pure const nothrow
- {
- return abs(signedArea());
- }
-
- /// Returns: Signed area of a 2D triangle.
- @nogc T signedArea() pure const nothrow
- {
- return ((b.x * a.y - a.x * b.y)
- + (c.x * b.y - b.x * c.y)
- + (a.x * c.y - c.x * a.y)) / 2;
- }
- }
-
- static if (N == 3 && isFloatingPoint!T)
- {
- /// Returns: Triangle normal.
- @nogc Vector!(T, 3) computeNormal() pure const nothrow
- {
- return cross(b - a, c - a).normalized();
- }
- }
- }
-}
-
-alias Triangle!(float, 2) triangle2f; /// 2D float triangle.
-alias Triangle!(float, 3) triangle3f; /// 3D float triangle.
-alias Triangle!(double, 2) triangle2d; /// 2D double triangle.
-alias Triangle!(double, 3) triangle3d; /// 3D double triangle.
-
-/// A Sphere is a point + a radius.
-struct Sphere(T, int N)
-{
- public nothrow
- {
- alias Vector!(T, N) point_t;
- point_t center;
- T radius;
-
- /// Creates a sphere from a point and a radius.
- @nogc this(in point_t center_, T radius_) pure nothrow
- {
- center = center_;
- radius = radius_;
- }
-
- /// Sphere contains point test.
- /// Returns: true if the point is inside the sphere.
- @nogc bool contains(in Sphere s) pure const nothrow
- {
- if (s.radius > radius)
- return false;
-
- T innerRadius = radius - s.radius;
- return squaredDistanceTo(s.center) < innerRadius * innerRadius;
- }
-
- /// Sphere vs point Euclidean distance squared.
- @nogc T squaredDistanceTo(point_t p) pure const nothrow
- {
- return center.squaredDistanceTo(p);
- }
-
- /// Sphere vs sphere intersection.
- /// Returns: true if the spheres intersect.
- @nogc bool intersects(Sphere s) pure const nothrow
- {
- T outerRadius = radius + s.radius;
- return squaredDistanceTo(s.center) < outerRadius * outerRadius;
- }
-
- static if (isFloatingPoint!T)
- {
- /// Sphere vs point Euclidean distance.
- @nogc T distanceTo(point_t p) pure const nothrow
- {
- return center.distanceTo(p);
- }
-
- static if(N == 2)
- {
- /// Returns: Circle area.
- @nogc T area() pure const nothrow
- {
- return T(PI) * (radius * radius);
- }
- }
- }
- }
-}
-
-alias Sphere!(float, 2) sphere2f; /// 2D float sphere (ie. a circle).
-alias Sphere!(float, 3) sphere3f; /// 3D float sphere.
-alias Sphere!(double, 2) sphere2d; /// 2D double sphere (ie. a circle).
-alias Sphere!(double, 3) sphere3d; /// 3D double sphere (ie. a circle).
-
-
-/// A Ray ir a point + a direction.
-struct Ray(T, int N)
-{
-nothrow:
- public
- {
- alias Vector!(T, N) point_t;
- point_t orig;
- point_t dir;
-
- /// Returns: A point further along the ray direction.
- @nogc point_t progress(T t) pure const nothrow
- {
- return orig + dir * t;
- }
-
- static if (N == 3 && isFloatingPoint!T)
- {
- /// Ray vs triangle intersection.
- /// See_also: "Fast, Minimum Storage Ray/Triangle intersection", Mommer & Trumbore (1997)
- /// Returns: Barycentric coordinates, the intersection point is at $(D (1 - u - v) * A + u * B + v * C).
- @nogc bool intersect(Triangle!(T, 3) triangle, out T t, out T u, out T v) pure const nothrow
- {
- point_t edge1 = triangle.b - triangle.a;
- point_t edge2 = triangle.c - triangle.a;
- point_t pvec = cross(dir, edge2);
- T det = dot(edge1, pvec);
- if (abs(det) < T.epsilon)
- return false; // no intersection
- T invDet = 1 / det;
-
- // calculate distance from triangle.a to ray origin
- point_t tvec = orig - triangle.a;
-
- // calculate U parameter and test bounds
- u = dot(tvec, pvec) * invDet;
- if (u < 0 || u > 1)
- return false;
-
- // prepare to test V parameter
- point_t qvec = cross(tvec, edge1);
-
- // calculate V parameter and test bounds
- v = dot(dir, qvec) * invDet;
- if (v < 0.0 || u + v > 1.0)
- return false;
-
- // calculate t, ray intersects triangle
- t = dot(edge2, qvec) * invDet;
- return true;
- }
-
- /// Ray vs plane intersection.
- /// See_also: "Geometry for Computer Graphics: Formulae, Examples
- /// and Proofs", Vince (2005), p. 62
- /// Returns: $(D true) if the ray crosses the plane exactly once,
- /// or $(D false) if the ray is parallel to or points away from
- /// the plane
- /// Params:
- /// plane = Plane to intersect.
- /// intersection = Point of intersection.
- /// distance = set to the point's distance along the ray, or
- /// T.infinity if parallel to the plane
- @nogc bool intersect(Plane!T plane, out point_t intersection, out T distance) pure const nothrow
- {
- // dot product will be 0 if angle to plane is 0
- T dp = dot(plane.n, dir);
- if (abs(dp) < T.epsilon)
- {
- distance = T.infinity;
- return false; // parallel to plane
- }
-
- // distance along ray where intersection happens
- distance = -(dot(plane.n, orig) - plane.d) / dp;
-
- intersection = distance*dir + orig;
- return distance >= 0;
- }
- }
- }
-}
-
-alias Ray!(float, 2) ray2f; /// 2D float ray.
-alias Ray!(float, 3) ray3f; /// 3D float ray.
-alias Ray!(double, 2) ray2d; /// 2D double ray.
-alias Ray!(double, 3) ray3d; /// 3D double ray.
-
-
-/// 3D plane.
-/// See_also: Flipcode article by Nate Miller $(WEB www.flipcode.com/archives/Plane_Class.shtml).
-struct Plane(T) if (isFloatingPoint!T)
-{
- public
- {
- vec3!T n; /// Normal (always stored normalized).
- T d;
-
- /// Create from four coordinates.
- @nogc this(vec4!T abcd) pure nothrow
- {
- n = vec3!T(abcd.x, abcd.y, abcd.z).normalized();
- d = abcd.w;
- }
-
- /// Create from a point and a normal.
- @nogc this(vec3!T origin, vec3!T normal) pure nothrow
- {
- n = normal.normalized();
- d = -dot(origin, n);
- }
-
- /// Create from 3 non-aligned points.
- @nogc this(vec3!T A, vec3!T B, vec3!T C) pure nothrow
- {
- this(C, cross(B - A, C - A));
- }
-
- /// Assign a plane with another plane.
- @nogc ref Plane opAssign(Plane other) pure nothrow
- {
- n = other.n;
- d = other.d;
- return this;
- }
-
- /// Returns: signed distance between a point and the plane.
- @nogc T signedDistanceTo(vec3!T point) pure const nothrow
- {
- return dot(n, point) + d;
- }
-
- /// Returns: absolute distance between a point and the plane.
- @nogc T distanceTo(vec3!T point) pure const nothrow
- {
- return abs(signedDistanceTo(point));
- }
-
- /// Returns: true if the point is in front of the plane.
- @nogc bool isFront(vec3!T point) pure const nothrow
- {
- return signedDistanceTo(point) >= 0;
- }
-
- /// Returns: true if the point is in the back of the plane.
- @nogc bool isBack(vec3!T point) pure const nothrow
- {
- return signedDistanceTo(point) < 0;
- }
-
- /// Returns: true if the point is on the plane, with a given epsilon.
- @nogc bool isOn(vec3!T point, T epsilon) pure const nothrow
- {
- T sd = signedDistanceTo(point);
- return (-epsilon < sd) && (sd < epsilon);
- }
- }
-}
-
-alias Plane!float planef; /// 3D float plane.
-alias Plane!double planed; /// 3D double plane.
-
-unittest
-{
- auto p = planed(vec4d(1.0, 2.0, 3.0, 4.0));
- auto p2 = planed(vec3d(1.0, 0.0, 0.0),
- vec3d(0.0, 1.0, 0.0),
- vec3d(0.0, 0.0, 1.0));
-
- assert(p2.isOn(vec3d(1.0, 0.0, 0.0), 1e-7));
- assert(p2.isFront(vec3d(1.0, 1.0, 1.0)));
- assert(p2.isBack(vec3d(0.0, 0.0, 0.0)));
-}
-
-/// Plane intersection tests
-@nogc pure nothrow unittest
-{
- void testR(planed p, ray3d r, bool shouldIntersect, double expectedDistance, vec3d expectedPoint = vec3d.init) pure nothrow @nogc
- {
- vec3d point;
- double distance;
- assert(r.intersect(p, point, distance) == shouldIntersect);
- assert(approxEqual(distance, expectedDistance));
- if (shouldIntersect)
- assert(approxEqual(point.v[], expectedPoint.v[]));
- }
- // ray facing plane
- testR(planed(vec4d(1.0, 0.0, 0.0, 1.0)), ray3d(vec3d(2.0, 3.0, 4.0), vec3d(-1.0, 0.0, 0.0)),
- true, 1.0, vec3d(1.0, 3.0, 4.0));
- testR(planed(vec4d(1.0, 1.0, 1.0, -sqrt(3.0))), ray3d(vec3d(1.0, 1.0, 1.0), vec3d(-1.0, -1.0, -1.0).normalized()),
- true, 2.0*sqrt(3.0), vec3d(-1.0, -1.0, -1.0));
- // ray facing away
- testR(planed(vec4d(1.0, 0.0, 0.0, 1.0)), ray3d(vec3d(2.0, 3.0, 4.0), vec3d(1.0, 0.0, 0.0)),
- false, -1.0);
- testR(planed(vec4d(1.0, 0.0, 0.0, 5.0)), ray3d(vec3d(2.0, 3.0, 4.0), vec3d(-1.0, 0.0, 0.0)),
- false, -3.0);
- // ray parallel
- testR(planed(vec4d(0.0, 0.0, 1.0, 1.0)), ray3d(vec3d(1.0, 2.0, 3.0), vec3d(1.0, 0.0, 0.0)),
- false, double.infinity);
-
- void testS(planed p, seg3d s, bool shouldIntersect, double expectedProgress, vec3d expectedPoint = vec3d.init) pure nothrow @nogc
- {
- vec3d point;
- double progress;
- assert(s.intersect(p, point, progress) == shouldIntersect);
- assert(approxEqual(progress, expectedProgress));
- if (shouldIntersect)
- assert(approxEqual(point.v[], expectedPoint.v[]));
- }
- // segment crossing plane
- testS(planed(vec4d(1.0, 0.0, 0.0, 2.0)), seg3d(vec3d(1.0, 2.0, 3.0), vec3d(3.0, 4.0, 5.0)),
- true, 0.5, vec3d(2.0, 3.0, 4.0));
- // segment too short
- testS(planed(vec4d(1.0, 0.0, 0.0, 0.0)), seg3d(vec3d(1.0, 2.0, 3.0), vec3d(3.0, 4.0, 5.0)),
- false, -0.5);
-}
-
-
-/// 3D frustum.
-/// See_also: Flipcode article by Dion Picco $(WEB www.flipcode.com/archives/Frustum_Culling.shtml).
-/// Bugs: verify proper signedness of half-spaces
-struct Frustum(T) if (isFloatingPoint!T)
-{
- public
- {
- enum int LEFT = 0,
- RIGHT = 1,
- TOP = 2,
- BOTTOM = 3,
- NEAR = 4,
- FAR = 5;
-
- Plane!T[6] planes;
-
- /// Create a frustum from 6 planes.
- @nogc this(Plane!T left, Plane!T right, Plane!T top, Plane!T bottom, Plane!T near, Plane!T far) pure nothrow
- {
- planes[LEFT] = left;
- planes[RIGHT] = right;
- planes[TOP] = top;
- planes[BOTTOM] = bottom;
- planes[NEAR] = near;
- planes[FAR] = far;
- }
-
- enum : int
- {
- OUTSIDE, /// object is outside the frustum
- INTERSECT, /// object intersects with the frustum
- INSIDE /// object is inside the frustum
- }
-
- /// Point vs frustum intersection.
- @nogc bool contains(vec3!T point) pure const nothrow
- {
- for(int i = 0; i < 6; ++i)
- {
- T distance = planes[i].signedDistanceTo(point);
-
- if(distance < 0)
- return false;
- }
- return true;
- }
-
- /// Sphere vs frustum intersection.
- /// Returns: Frustum.OUTSIDE, Frustum.INTERSECT or Frustum.INSIDE.
- @nogc int contains(Sphere!(T, 3) sphere) pure const nothrow
- {
- // calculate our distances to each of the planes
- for(int i = 0; i < 6; ++i)
- {
- // find the distance to this plane
- T distance = planes[i].signedDistanceTo(sphere.center);
-
- if(distance < -sphere.radius)
- return OUTSIDE;
-
- else if (distance < sphere.radius)
- return INTERSECT;
- }
-
- // otherwise we are fully in view
- return INSIDE;
- }
-
- /// AABB vs frustum intersection.
- /// Returns: Frustum.OUTSIDE, Frustum.INTERSECT or Frustum.INSIDE.
- @nogc int contains(box3!T box) pure const nothrow
- {
- vec3!T[8] corners;
- int totalIn = 0;
-
- for (int i = 0; i < 2; ++i)
- for (int j = 0; j < 2; ++j)
- for (int k = 0; k < 2; ++k)
- {
- auto x = i == 0 ? box.min.x : box.max.x;
- auto y = j == 0 ? box.min.y : box.max.y;
- auto z = k == 0 ? box.min.z : box.max.z;
- corners[i*4 + j*2 + k] = vec3!T(x, y, z);
- }
-
- // test all 8 corners against the 6 sides
- // if all points are behind 1 specific plane, we are out
- // if we are in with all points, then we are fully in
- for(int p = 0; p < 6; ++p)
- {
- int inCount = 8;
- int ptIn = 1;
-
- for(int i = 0; i < 8; ++i)
- {
- // test this point against the planes
- if (planes[p].isBack(corners[i]))
- {
- ptIn = 0;
- --inCount;
- }
- }
-
- // were all the points outside of plane p?
- if (inCount == 0)
- return OUTSIDE;
-
- // check if they were all on the right side of the plane
- totalIn += ptIn;
- }
-
- // so if totalIn is 6, then all are inside the view
- if(totalIn == 6)
- return INSIDE;
-
- // we must be partly in then otherwise
- return INTERSECT;
- }
-
- }
-}
-
-unittest
-{
- seg2f se;
- triangle3f tr;
- Frustum!double frust;
- planed pl;
-}
-
-/// True if `T` is a kind of Segment
-enum isSegment(T) = is(T : Segment!U, U...);
-
-/// True if `T` is a kind of Triangle
-enum isTriangle(T) = is(T : Triangle!U, U...);
-
-/// True if `T` is a kind of Sphere
-enum isSphere(T) = is(T : Sphere!U, U...);
-
-/// True if `T` is a kind of Ray
-enum isRay(T) = is(T : Ray!U, U...);
-
-/// True if `T` is a kind of Plane
-enum isPlane(T) = is(T : Plane!U, U);
-
-/// True if `T` is a kind of Frustum
-enum isFrustum(T) = is(T : Frustum!U, U);
-
-/// True if `T` is a kind of 2 dimensional Segment
-enum isSegment2D(T) = is(T : Segment!(U, 2), U);
-
-/// True if `T` is a kind of 2 dimensional Triangle
-enum isTriangle2D(T) = is(T : Triangle!(U, 2), U);
-
-/// True if `T` is a kind of 2 dimensional Sphere
-enum isSphere2D(T) = is(T : Sphere!(U, 2), U);
-
-/// True if `T` is a kind of 2 dimensional Ray
-enum isRay2D(T) = is(T : Ray!(U, 2), U);
-
-/// True if `T` is a kind of 3 dimensional Segment
-enum isSegment3D(T) = is(T : Segment!(U, 3), U);
-
-/// True if `T` is a kind of 3 dimensional Triangle
-enum isTriangle3D(T) = is(T : Triangle!(U, 3), U);
-
-/// True if `T` is a kind of 3 dimensional Sphere
-enum isSphere3D(T) = is(T : Sphere!(U, 3), U);
-
-/// True if `T` is a kind of 3 dimensional Ray
-enum isRay3D(T) = is(T : Ray!(U, 3), U);
-
-unittest
-{
- enum ShapeType
- {
- segment,
- triangle,
- sphere,
- ray,
- plane,
- frustum
- }
-
- void test(T, ShapeType type, int dim)()
- {
- with (ShapeType)
- {
- static assert(isSegment!T == (type == segment ));
- static assert(isTriangle!T == (type == triangle));
- static assert(isSphere!T == (type == sphere ));
- static assert(isRay!T == (type == ray ));
- static assert(isPlane!T == (type == plane ));
- static assert(isFrustum!T == (type == frustum ));
-
- static assert(isSegment2D!T == (type == segment && dim == 2));
- static assert(isTriangle2D!T == (type == triangle && dim == 2));
- static assert(isSphere2D!T == (type == sphere && dim == 2));
- static assert(isRay2D!T == (type == ray && dim == 2));
-
- static assert(isSegment3D!T == (type == segment && dim == 3));
- static assert(isTriangle3D!T == (type == triangle && dim == 3));
- static assert(isSphere3D!T == (type == sphere && dim == 3));
- static assert(isRay3D!T == (type == ray && dim == 3));
- }
- }
-
- with (ShapeType)
- {
- // test case type #dimensions
- test!(seg2f , segment , 2);
- test!(seg3d , segment , 3);
- test!(triangle2d , triangle, 2);
- test!(triangle3f , triangle, 3);
- test!(sphere2d , sphere , 2);
- test!(Sphere!(uint, 3), sphere , 3);
- test!(ray2f , ray , 2);
- test!(Ray!(real, 3) , ray , 3);
- test!(planed , plane , 0); // ignore dimension (always 3D)
- test!(Plane!float , plane , 0);
- test!(Frustum!double , frustum , 0);
- }
-}
-
-/// Get the numeric type used to measure a shape's dimensions.
-alias DimensionType(T : Segment!U, U...) = U[0];
-/// ditto
-alias DimensionType(T : Triangle!U, U...) = U[0];
-/// ditto
-alias DimensionType(T : Sphere!U, U...) = U[0];
-/// ditto
-alias DimensionType(T : Ray!U, U...) = U[0];
-/// ditto
-alias DimensionType(T : Plane!U, U) = U;
-/// ditto
-alias DimensionType(T : Frustum!U, U) = U;
-
-///
-unittest
-{
- static assert(is(DimensionType!seg2i == int));
- static assert(is(DimensionType!triangle3d == double));
- static assert(is(DimensionType!sphere2d == double));
- static assert(is(DimensionType!ray3f == float));
- static assert(is(DimensionType!planed == double));
- static assert(is(DimensionType!(Frustum!real) == real));
-}
diff --git a/external/gfm/math/simplerng.d b/external/gfm/math/simplerng.d
deleted file mode 100644
index c2dbb6b..0000000
--- a/external/gfm/math/simplerng.d
+++ /dev/null
@@ -1,485 +0,0 @@
-/**
- Translation of SimpleRNG.
- Removed the builtin RNG to use std.random, but kept the distribution functions.
- John D. Cook confirmed this code as public domain.
-
- Authors: John D. Cook.
- See_also: $(WEB www.johndcook.com/cpp_random_number_generation.html)
- */
-module gfm.math.simplerng;
-
-public import std.random;
-import std.math;
-
-/// Returns: Normal (Gaussian) random sample.
-/// See_also: Box-Muller algorithm.
-double randNormal(RNG)(ref RNG rng, double mean = 0.0, double standardDeviation = 1.0)
-{
- assert(standardDeviation > 0);
- double u1;
-
- do
- {
- u1 = uniform01(rng);
- } while (u1 == 0); // u1 must not be zero
- double u2 = uniform01(rng);
- double r = sqrt(-2.0 * log(u1));
- double theta = 2.0 * double(PI) * u2;
- return mean + standardDeviation * r * sin(theta);
-}
-
-/// Returns: Exponential random sample with specified mean.
-double randExponential(RNG)(ref RNG rng, double mean = 1.0)
-{
- assert(mean > 0);
- return -mean*log(uniform01(rng));
-}
-
-/// Returns: Gamma random sample.
-/// See_also: "A Simple Method for Generating Gamma Variables"
-/// by George Marsaglia and Wai Wan Tsang. ACM Transactions on Mathematical Software
-/// Vol 26, No 3, September 2000, pages 363-372.
-double randGamma(RNG)(ref RNG rng, double shape, double scale)
-{
- double d, c, x, xsquared, v, u;
-
- if (shape >= 1.0)
- {
- d = shape - 1.0/3.0;
- c = 1.0/sqrt(9.0*d);
- for (;;)
- {
- do
- {
- x = randNormal(rng);
- v = 1.0 + c*x;
- }
- while (v <= 0.0);
- v = v*v*v;
- u = uniform01(rng);
- xsquared = x*x;
- if (u < 1.0 -.0331*xsquared*xsquared || log(u) < 0.5*xsquared + d*(1.0 - v + log(v)))
- return scale*d*v;
- }
- }
- else
- {
- assert(shape > 0);
- double g = randGamma(rng, shape+1.0, 1.0);
- double w = uniform01(rng);
- return scale*g*pow(w, 1.0/shape);
- }
-}
-
-/// Returns: Chi-square sample.
-double randChiSquare(RNG)(ref RNG rng, double degreesOfFreedom)
-{
- // A chi squared distribution with n degrees of freedom
- // is a gamma distribution with shape n/2 and scale 2.
- return randGamma(rng, 0.5 * degreesOfFreedom, 2.0);
-}
-
-/// Returns: Inverse-gamma sample.
-double randInverseGamma(RNG)(ref RNG rng, double shape, double scale)
-{
- // If X is gamma(shape, scale) then
- // 1/Y is inverse gamma(shape, 1/scale)
- return 1.0 / randGamma(rng, shape, 1.0 / scale);
-}
-
-/// Returns: Weibull sample.
-double randWeibull(RNG)(ref RNG rng, double shape, double scale)
-{
- assert(shape > 0 && scale > 0);
- return scale * pow(-log(uniform01(rng)), 1.0 / shape);
-}
-
-/// Returns: Cauchy sample.
-double randCauchy(RNG)(ref RNG rng, double median, double scale)
-{
- assert(scale > 0);
- double p = uniform01(rng);
-
- // Apply inverse of the Cauchy distribution function to a uniform
- return median + scale*tan(double(PI)*(p - 0.5));
-}
-
-/// Returns: Student-t sample.
-/// See_also: Seminumerical Algorithms by Knuth.
-double randStudentT(RNG)(ref RNG rng, double degreesOfFreedom)
-{
- assert(degreesOfFreedom > 0);
-
-
- double y1 = getNormal(rng);
- double y2 = getChiSquare(rng, degreesOfFreedom);
- return y1 / sqrt(y2 / degreesOfFreedom);
-}
-
-/// Returns: Laplace distribution random sample (also known as the double exponential distribution).
-double randLaplace(RNG)(ref RNG rng, double mean, double scale)
-{
- double u = uniform01(rng);
- return (u < 0.5) ? (mean + scale*log(2.0*u))
- : (mean - scale*log(2*(1-u)));
-}
-
-/// Returns: Log-normal sample.
-double randLogNormal(RNG)(ref RNG rng, double mu, double sigma)
-{
- return exp(getNormal(rng, mu, sigma));
-}
-
-/// Returns: Beta sample.
-double randBeta(RNG)(ref RNG rng, double a, double b)
-{
- assert(a > 0 && b > 0);
-
- // There are more efficient methods for generating beta samples.
- // However such methods are a little more efficient and much more complicated.
- // For an explanation of why the following method works, see
- // http://www.johndcook.com/distribution_chart.html#gamma_beta
-
- double u = getGamma(rng, a, 1.0);
- double v = getGamma(rng, b, 1.0);
- return u / (u + v);
-}
-
-/// Returns: Poisson sample.
-int randPoisson(RNG)(ref RNG rng, double lambda)
-{
- return (lambda < 30.0) ? poissonSmall(rng, lambda) : poissonLarge(rng, lambda);
-}
-
-private
-{
- int poissonSmall(RNG)(ref RNG rng, double lambda)
- {
- // Algorithm due to Donald Knuth, 1969.
- double p = 1.0, L = exp(-lambda);
- int k = 0;
- do
- {
- k++;
- p *= uniform01(rng);
- }
- while (p > L);
- return k - 1;
- }
-
- int poissonLarge(RNG)(ref RNG rng, double lambda)
- {
- // "Rejection method PA" from "The Computer Generation of Poisson Random Variables" by A. C. Atkinson
- // Journal of the Royal Statistical Society Series C (Applied Statistics) Vol. 28, No. 1. (1979)
- // The article is on pages 29-35. The algorithm given here is on page 32.
-
- double c = 0.767 - 3.36/lambda;
- double beta = double(PI)/sqrt(3.0*lambda);
- double alpha = beta*lambda;
- double k = log(c) - lambda - log(beta);
-
- for(;;)
- {
- double u = uniform01(rng);
- double x = (alpha - log((1.0 - u)/u))/beta;
- int n = cast(int)(floor(x + 0.5));
- if (n < 0)
- continue;
- double v = uniform01(rng);
- double y = alpha - beta*x;
- double temp = 1.0 + exp(y);
- double lhs = y + log(v/(temp*temp));
- double rhs = k + n*log(lambda) - logFactorial(n);
- if (lhs <= rhs)
- return n;
- }
- }
-
- double logFactorial(int n) nothrow
- {
- assert(n >= 0);
- if (n > 254)
- {
- double x = n + 1;
- return (x - 0.5)*log(x) - x + 0.5*log(2*double(PI)) + 1.0/(12.0*x);
- }
- else
- {
- return LOG_FACTORIAL[n];
- }
- }
-}
-
-private static immutable double[255] LOG_FACTORIAL =
-[
- 0.000000000000000,
- 0.000000000000000,
- 0.693147180559945,
- 1.791759469228055,
- 3.178053830347946,
- 4.787491742782046,
- 6.579251212010101,
- 8.525161361065415,
- 10.604602902745251,
- 12.801827480081469,
- 15.104412573075516,
- 17.502307845873887,
- 19.987214495661885,
- 22.552163853123421,
- 25.191221182738683,
- 27.899271383840894,
- 30.671860106080675,
- 33.505073450136891,
- 36.395445208033053,
- 39.339884187199495,
- 42.335616460753485,
- 45.380138898476908,
- 48.471181351835227,
- 51.606675567764377,
- 54.784729398112319,
- 58.003605222980518,
- 61.261701761002001,
- 64.557538627006323,
- 67.889743137181526,
- 71.257038967168000,
- 74.658236348830158,
- 78.092223553315307,
- 81.557959456115029,
- 85.054467017581516,
- 88.580827542197682,
- 92.136175603687079,
- 95.719694542143202,
- 99.330612454787428,
- 102.968198614513810,
- 106.631760260643450,
- 110.320639714757390,
- 114.034211781461690,
- 117.771881399745060,
- 121.533081515438640,
- 125.317271149356880,
- 129.123933639127240,
- 132.952575035616290,
- 136.802722637326350,
- 140.673923648234250,
- 144.565743946344900,
- 148.477766951773020,
- 152.409592584497350,
- 156.360836303078800,
- 160.331128216630930,
- 164.320112263195170,
- 168.327445448427650,
- 172.352797139162820,
- 176.395848406997370,
- 180.456291417543780,
- 184.533828861449510,
- 188.628173423671600,
- 192.739047287844900,
- 196.866181672889980,
- 201.009316399281570,
- 205.168199482641200,
- 209.342586752536820,
- 213.532241494563270,
- 217.736934113954250,
- 221.956441819130360,
- 226.190548323727570,
- 230.439043565776930,
- 234.701723442818260,
- 238.978389561834350,
- 243.268849002982730,
- 247.572914096186910,
- 251.890402209723190,
- 256.221135550009480,
- 260.564940971863220,
- 264.921649798552780,
- 269.291097651019810,
- 273.673124285693690,
- 278.067573440366120,
- 282.474292687630400,
- 286.893133295426990,
- 291.323950094270290,
- 295.766601350760600,
- 300.220948647014100,
- 304.686856765668720,
- 309.164193580146900,
- 313.652829949878990,
- 318.152639620209300,
- 322.663499126726210,
- 327.185287703775200,
- 331.717887196928470,
- 336.261181979198450,
- 340.815058870798960,
- 345.379407062266860,
- 349.954118040770250,
- 354.539085519440790,
- 359.134205369575340,
- 363.739375555563470,
- 368.354496072404690,
- 372.979468885689020,
- 377.614197873918670,
- 382.258588773060010,
- 386.912549123217560,
- 391.575988217329610,
- 396.248817051791490,
- 400.930948278915760,
- 405.622296161144900,
- 410.322776526937280,
- 415.032306728249580,
- 419.750805599544780,
- 424.478193418257090,
- 429.214391866651570,
- 433.959323995014870,
- 438.712914186121170,
- 443.475088120918940,
- 448.245772745384610,
- 453.024896238496130,
- 457.812387981278110,
- 462.608178526874890,
- 467.412199571608080,
- 472.224383926980520,
- 477.044665492585580,
- 481.872979229887900,
- 486.709261136839360,
- 491.553448223298010,
- 496.405478487217580,
- 501.265290891579240,
- 506.132825342034830,
- 511.008022665236070,
- 515.890824587822520,
- 520.781173716044240,
- 525.679013515995050,
- 530.584288294433580,
- 535.496943180169520,
- 540.416924105997740,
- 545.344177791154950,
- 550.278651724285620,
- 555.220294146894960,
- 560.169054037273100,
- 565.124881094874350,
- 570.087725725134190,
- 575.057539024710200,
- 580.034272767130800,
- 585.017879388839220,
- 590.008311975617860,
- 595.005524249382010,
- 600.009470555327430,
- 605.020105849423770,
- 610.037385686238740,
- 615.061266207084940,
- 620.091704128477430,
- 625.128656730891070,
- 630.172081847810200,
- 635.221937855059760,
- 640.278183660408100,
- 645.340778693435030,
- 650.409682895655240,
- 655.484856710889060,
- 660.566261075873510,
- 665.653857411105950,
- 670.747607611912710,
- 675.847474039736880,
- 680.953419513637530,
- 686.065407301994010,
- 691.183401114410800,
- 696.307365093814040,
- 701.437263808737160,
- 706.573062245787470,
- 711.714725802289990,
- 716.862220279103440,
- 722.015511873601330,
- 727.174567172815840,
- 732.339353146739310,
- 737.509837141777440,
- 742.685986874351220,
- 747.867770424643370,
- 753.055156230484160,
- 758.248113081374300,
- 763.446610112640200,
- 768.650616799717000,
- 773.860102952558460,
- 779.075038710167410,
- 784.295394535245690,
- 789.521141208958970,
- 794.752249825813460,
- 799.988691788643450,
- 805.230438803703120,
- 810.477462875863580,
- 815.729736303910160,
- 820.987231675937890,
- 826.249921864842800,
- 831.517780023906310,
- 836.790779582469900,
- 842.068894241700490,
- 847.352097970438420,
- 852.640365001133090,
- 857.933669825857460,
- 863.231987192405430,
- 868.535292100464630,
- 873.843559797865740,
- 879.156765776907600,
- 884.474885770751830,
- 889.797895749890240,
- 895.125771918679900,
- 900.458490711945270,
- 905.796028791646340,
- 911.138363043611210,
- 916.485470574328820,
- 921.837328707804890,
- 927.193914982476710,
- 932.555207148186240,
- 937.921183163208070,
- 943.291821191335660,
- 948.667099599019820,
- 954.046996952560450,
- 959.431492015349480,
- 964.820563745165940,
- 970.214191291518320,
- 975.612353993036210,
- 981.015031374908400,
- 986.422203146368590,
- 991.833849198223450,
- 997.249949600427840,
- 1002.670484599700300,
- 1008.095434617181700,
- 1013.524780246136200,
- 1018.958502249690200,
- 1024.396581558613400,
- 1029.838999269135500,
- 1035.285736640801600,
- 1040.736775094367400,
- 1046.192096209724900,
- 1051.651681723869200,
- 1057.115513528895000,
- 1062.583573670030100,
- 1068.055844343701400,
- 1073.532307895632800,
- 1079.012946818975000,
- 1084.497743752465600,
- 1089.986681478622400,
- 1095.479742921962700,
- 1100.976911147256000,
- 1106.478169357800900,
- 1111.983500893733000,
- 1117.492889230361000,
- 1123.006317976526100,
- 1128.523770872990800,
- 1134.045231790853000,
- 1139.570684729984800,
- 1145.100113817496100,
- 1150.633503306223700,
- 1156.170837573242400,
-];
-
-unittest
-{
- Xorshift32 rng;
- rng.seed(unpredictableSeed());
-
- double x = randNormal!Xorshift32(rng, 0.0, 1.0);
- x = randExponential!Xorshift32(rng);
- x = randGamma!Xorshift32(rng, 1.2, 1.0);
- x = randGamma!Xorshift32(rng, 0.8, 2.0);
- x = randChiSquare!Xorshift32(rng, 2.0);
- x = randInverseGamma!Xorshift32(rng, 1.1, 0.7);
- x = randWeibull!Xorshift32(rng, 3.0, 0.7);
- x = randCauchy!Xorshift32(rng, 5.0, 1.4);
-}
diff --git a/external/gfm/math/vector.d b/external/gfm/math/vector.d
deleted file mode 100644
index 4578db7..0000000
--- a/external/gfm/math/vector.d
+++ /dev/null
@@ -1,770 +0,0 @@
-/// N-dimension vector mathematical object
-module gfm.math.vector;
-
-import std.traits,
- std.math,
- std.conv,
- std.array,
- std.string;
-
-import gfm.math.funcs;
-
-/**
- * Generic 1D small vector.
- * Params:
- * N = number of elements
- * T = type of elements
- */
-struct Vector(T, int N)
-{
-nothrow:
- public
- {
- static assert(N >= 1);
-
- // fields definition
- union
- {
- T[N] v;
- struct
- {
- static if (N >= 1)
- {
- T x;
- alias x r;
- }
- static if (N >= 2)
- {
- T y;
- alias y g;
- }
- static if (N >= 3)
- {
- T z;
- alias z b;
- }
- static if (N >= 4)
- {
- T w;
- alias w a;
- }
- }
- }
-
- /// Construct a Vector with a `T[]` or the values as arguments
- @nogc this(Args...)(Args args) pure nothrow
- {
- static if (args.length == 1)
- {
- // Construct a Vector from a single value.
- opAssign!(Args[0])(args[0]);
- }
- else
- {
- // validate the total argument count across scalars and vectors
- template argCount(T...) {
- static if(T.length == 0)
- enum argCount = 0; // done recursing
- else static if(isVector!(T[0]))
- enum argCount = T[0]._N + argCount!(T[1..$]);
- else
- enum argCount = 1 + argCount!(T[1..$]);
- }
-
- static assert(argCount!Args <= N, "Too many arguments in vector constructor");
-
- int index = 0;
- foreach(arg; args)
- {
- static if (isAssignable!(T, typeof(arg)))
- {
- v[index] = arg;
- index++; // has to be on its own line (DMD 2.068)
- }
- else static if (isVector!(typeof(arg)) && isAssignable!(T, arg._T))
- {
- mixin(generateLoopCode!("v[index + @] = arg[@];", arg._N)());
- index += arg._N;
- }
- else
- static assert(false, "Unrecognized argument in Vector constructor");
- }
- assert(index == N, "Bad arguments in Vector constructor");
- }
- }
-
- size_t toHash() const nothrow @safe
- {
- size_t hash = 0;
- foreach (elem; v) {
- hash = elem.hashOf(hash);
- }
- return hash;
- }
-
- /// Assign a Vector from a compatible type.
- @nogc ref Vector opAssign(U)(U x) pure nothrow if (isAssignable!(T, U))
- {
- mixin(generateLoopCode!("v[@] = x;", N)()); // copy to each component
- return this;
- }
-
- /// Assign a Vector with a static array type.
- @nogc ref Vector opAssign(U)(U arr) pure nothrow if ((isStaticArray!(U) && isAssignable!(T, typeof(arr[0])) && (arr.length == N)))
- {
- mixin(generateLoopCode!("v[@] = arr[@];", N)());
- return this;
- }
-
- /// Assign with a dynamic array.
- /// Size is checked in debug-mode.
- @nogc ref Vector opAssign(U)(U arr) pure nothrow if (isDynamicArray!(U) && isAssignable!(T, typeof(arr[0])))
- {
- assert(arr.length == N);
- mixin(generateLoopCode!("v[@] = arr[@];", N)());
- return this;
- }
-
- /// Assign from a samey Vector.
- @nogc ref Vector opAssign(U)(U u) pure nothrow if (is(U : Vector))
- {
- v[] = u.v[];
- return this;
- }
-
- /// Assign from other vectors types (same size, compatible type).
- @nogc ref Vector opAssign(U)(U x) pure nothrow if (isVector!U
- && isAssignable!(T, U._T)
- && (!is(U: Vector))
- && (U._N == _N))
- {
- mixin(generateLoopCode!("v[@] = x.v[@];", N)());
- return this;
- }
-
- /// Returns: a pointer to content.
- @nogc inout(T)* ptr() pure inout nothrow @property
- {
- return v.ptr;
- }
-
- /// Converts to a pretty string.
- string toString() const nothrow
- {
- try
- return format("%s", v);
- catch (Exception e)
- assert(false); // should not happen since format is right
- }
-
- @nogc bool opEquals(U)(U other) pure const nothrow
- if (is(U : Vector))
- {
- for (int i = 0; i < N; ++i)
- {
- if (v[i] != other.v[i])
- {
- return false;
- }
- }
- return true;
- }
-
- @nogc bool opEquals(U)(U other) pure const nothrow
- if (isConvertible!U)
- {
- Vector conv = other;
- return opEquals(conv);
- }
-
- @nogc Vector opUnary(string op)() pure const nothrow
- if (op == "+" || op == "-" || op == "~" || op == "!")
- {
- Vector res = void;
- mixin(generateLoopCode!("res.v[@] = " ~ op ~ " v[@];", N)());
- return res;
- }
-
- @nogc ref Vector opOpAssign(string op, U)(U operand) pure nothrow
- if (is(U : Vector))
- {
- mixin(generateLoopCode!("v[@] " ~ op ~ "= operand.v[@];", N)());
- return this;
- }
-
- @nogc ref Vector opOpAssign(string op, U)(U operand) pure nothrow if (isConvertible!U)
- {
- Vector conv = operand;
- return opOpAssign!op(conv);
- }
-
- @nogc Vector opBinary(string op, U)(U operand) pure const nothrow
- if (is(U: Vector) || (isConvertible!U))
- {
- Vector result = void;
- static if (is(U: T))
- mixin(generateLoopCode!("result.v[@] = cast(T)(v[@] " ~ op ~ " operand);", N)());
- else
- {
- Vector other = operand;
- mixin(generateLoopCode!("result.v[@] = cast(T)(v[@] " ~ op ~ " other.v[@]);", N)());
- }
- return result;
- }
-
- @nogc Vector opBinaryRight(string op, U)(U operand) pure const nothrow if (isConvertible!U)
- {
- Vector result = void;
- static if (is(U: T))
- mixin(generateLoopCode!("result.v[@] = cast(T)(operand " ~ op ~ " v[@]);", N)());
- else
- {
- Vector other = operand;
- mixin(generateLoopCode!("result.v[@] = cast(T)(other.v[@] " ~ op ~ " v[@]);", N)());
- }
- return result;
- }
-
- @nogc ref T opIndex(size_t i) pure nothrow
- {
- return v[i];
- }
-
- @nogc ref const(T) opIndex(size_t i) pure const nothrow
- {
- return v[i];
- }
-
- @nogc T opIndexAssign(U : T)(U x, size_t i) pure nothrow
- {
- return v[i] = x;
- }
-
-
- /// Implements swizzling.
- ///
- /// Example:
- /// ---
- /// vec4i vi = [4, 1, 83, 10];
- /// assert(vi.zxxyw == [83, 4, 4, 1, 10]);
- /// ---
- @nogc @property auto opDispatch(string op, U = void)() pure const nothrow if (isValidSwizzle!(op))
- {
- alias Vector!(T, op.length) returnType;
- returnType res = void;
- enum indexTuple = swizzleTuple!op;
- foreach(i, index; indexTuple)
- res.v[i] = v[index];
- return res;
- }
-
- /// Support swizzling assignment like in shader languages.
- ///
- /// Example:
- /// ---
- /// vec3f v = [0, 1, 2];
- /// v.yz = v.zx;
- /// assert(v == [0, 2, 0]);
- /// ---
- @nogc @property void opDispatch(string op, U)(U x) pure
- if ((op.length >= 2)
- && (isValidSwizzleUnique!op) // v.xyy will be rejected
- && is(typeof(Vector!(T, op.length)(x)))) // can be converted to a small vector of the right size
- {
- Vector!(T, op.length) conv = x;
- enum indexTuple = swizzleTuple!op;
- foreach(i, index; indexTuple)
- v[index] = conv[i];
- }
-
- /// Casting to small vectors of the same size.
- /// Example:
- /// ---
- /// vec4f vf;
- /// vec4d vd = cast!(vec4d)vf;
- /// ---
- @nogc U opCast(U)() pure const nothrow if (isVector!U && (U._N == _N))
- {
- U res = void;
- mixin(generateLoopCode!("res.v[@] = cast(U._T)v[@];", N)());
- return res;
- }
-
- /// Implement slices operator overloading.
- /// Allows to go back to slice world.
- /// Returns: length.
- @nogc int opDollar() pure const nothrow
- {
- return N;
- }
-
- /// Slice containing vector values
- /// Returns: a slice which covers the whole Vector.
- @nogc T[] opSlice() pure nothrow
- {
- return v[];
- }
-
- /// vec[a..b]
- @nogc T[] opSlice(int a, int b) pure nothrow
- {
- return v[a..b];
- }
-
- /// Squared Euclidean length of the Vector
- /// Returns: squared length.
- @nogc T squaredMagnitude() pure const nothrow
- {
- T sumSquares = 0;
- mixin(generateLoopCode!("sumSquares += v[@] * v[@];", N)());
- return sumSquares;
- }
-
- /// Squared Euclidean distance between this vector and another one
- /// Returns: squared Euclidean distance.
- @nogc T squaredDistanceTo(Vector v) pure const nothrow
- {
- return (v - this).squaredMagnitude();
- }
-
- static if (isFloatingPoint!T)
- {
- /// Euclidean length of the vector
- /// Returns: Euclidean length
- @nogc T magnitude() pure const nothrow
- {
- return sqrt(squaredMagnitude());
- }
-
- /// Inverse Euclidean length of the vector
- /// Returns: Inverse of Euclidean length.
- @nogc T inverseMagnitude() pure const nothrow
- {
- return 1 / sqrt(squaredMagnitude());
- }
-
- alias fastInverseLength = fastInverseMagnitude;
- /// Faster but less accurate inverse of Euclidean length.
- /// Returns: Inverse of Euclidean length.
- @nogc T fastInverseMagnitude() pure const nothrow
- {
- return inverseSqrt(squaredMagnitude());
- }
-
- /// Euclidean distance between this vector and another one
- /// Returns: Euclidean distance between this and other.
- @nogc T distanceTo(Vector other) pure const nothrow
- {
- return (other - this).magnitude();
- }
-
- /// In-place normalization.
- @nogc void normalize() pure nothrow
- {
- auto invMag = inverseMagnitude();
- mixin(generateLoopCode!("v[@] *= invMag;", N)());
- }
-
- /// Returns a normalized copy of this Vector
- /// Returns: Normalized vector.
- @nogc Vector normalized() pure const nothrow
- {
- Vector res = this;
- res.normalize();
- return res;
- }
-
- /// Faster but less accurate in-place normalization.
- @nogc void fastNormalize() pure nothrow
- {
- auto invLength = fastInverseMagnitude();
- mixin(generateLoopCode!("v[@] *= invLength;", N)());
- }
-
- /// Faster but less accurate vector normalization.
- /// Returns: Normalized vector.
- @nogc Vector fastNormalized() pure const nothrow
- {
- Vector res = this;
- res.fastNormalize();
- return res;
- }
-
- static if (N == 3)
- {
- /// Gets an orthogonal vector from a 3-dimensional vector.
- /// Doesn’t normalize the output.
- /// Authors: Sam Hocevar
- /// See_also: Source at $(WEB lolengine.net/blog/2013/09/21/picking-orthogonal-vector-combing-coconuts).
- @nogc Vector getOrthogonalVector() pure const nothrow
- {
- return abs(x) > abs(z) ? Vector(-y, x, 0.0) : Vector(0.0, -z, y);
- }
- }
- }
- }
-
- private
- {
- enum _N = N;
- alias T _T;
-
- // define types that can be converted to this, but are not the same type
- template isConvertible(T)
- {
- enum bool isConvertible = (!is(T : Vector))
- && is(typeof(
- {
- T x;
- Vector v = x;
- }()));
- }
-
- // define types that can't be converted to this
- template isForeign(T)
- {
- enum bool isForeign = (!isConvertible!T) && (!is(T: Vector));
- }
-
- template isValidSwizzle(string op, int lastSwizzleClass = -1)
- {
- static if (op.length == 0)
- enum bool isValidSwizzle = true;
- else
- {
- enum len = op.length;
- enum int swizzleClass = swizzleClassify!(op[0]);
- enum bool swizzleClassValid = (lastSwizzleClass == -1 || (swizzleClass == lastSwizzleClass));
- enum bool isValidSwizzle = (swizzleIndex!(op[0]) != -1)
- && swizzleClassValid
- && isValidSwizzle!(op[1..len], swizzleClass);
- }
- }
-
- template searchElement(char c, string s)
- {
- static if (s.length == 0)
- {
- enum bool result = false;
- }
- else
- {
- enum string tail = s[1..s.length];
- enum bool result = (s[0] == c) || searchElement!(c, tail).result;
- }
- }
-
- template hasNoDuplicates(string s)
- {
- static if (s.length == 1)
- {
- enum bool result = true;
- }
- else
- {
- enum tail = s[1..s.length];
- enum bool result = !(searchElement!(s[0], tail).result) && hasNoDuplicates!(tail).result;
- }
- }
-
- // true if the swizzle has at the maximum one time each letter
- template isValidSwizzleUnique(string op)
- {
- static if (isValidSwizzle!op)
- enum isValidSwizzleUnique = hasNoDuplicates!op.result;
- else
- enum bool isValidSwizzleUnique = false;
- }
-
- template swizzleIndex(char c)
- {
- static if((c == 'x' || c == 'r') && N >= 1)
- enum swizzleIndex = 0;
- else static if((c == 'y' || c == 'g') && N >= 2)
- enum swizzleIndex = 1;
- else static if((c == 'z' || c == 'b') && N >= 3)
- enum swizzleIndex = 2;
- else static if ((c == 'w' || c == 'a') && N >= 4)
- enum swizzleIndex = 3;
- else
- enum swizzleIndex = -1;
- }
-
- template swizzleClassify(char c)
- {
- static if(c == 'x' || c == 'y' || c == 'z' || c == 'w')
- enum swizzleClassify = 0;
- else static if(c == 'r' || c == 'g' || c == 'b' || c == 'a')
- enum swizzleClassify = 1;
- else
- enum swizzleClassify = -1;
- }
-
- template swizzleTuple(string op)
- {
- enum opLength = op.length;
- static if (op.length == 0)
- enum swizzleTuple = [];
- else
- enum swizzleTuple = [ swizzleIndex!(op[0]) ] ~ swizzleTuple!(op[1..op.length]);
- }
- }
-}
-
-/// True if `T` is some kind of `Vector`
-enum isVector(T) = is(T : Vector!U, U...);
-
-///
-unittest
-{
- static assert(isVector!vec2f);
- static assert(isVector!vec3d);
- static assert(isVector!(vec4!real));
- static assert(!isVector!float);
-}
-
-/// Get the numeric type used to measure a vectors's coordinates.
-alias DimensionType(T : Vector!U, U...) = U[0];
-
-///
-unittest
-{
- static assert(is(DimensionType!vec2f == float));
- static assert(is(DimensionType!vec3d == double));
-}
-
-///
-template vec2(T) { alias Vector!(T, 2) vec2; }
-///
-template vec3(T) { alias Vector!(T, 3) vec3; }
-///
-template vec4(T) { alias Vector!(T, 4) vec4; }
-
-alias vec2!int vec2i; ///
-alias vec2!float vec2f; ///
-alias vec2!double vec2d; ///
-
-alias vec3!int vec3i; ///
-alias vec3!float vec3f; ///
-alias vec3!double vec3d; ///
-
-alias vec4!int vec4i; ///
-alias vec4!float vec4f; ///
-alias vec4!double vec4d; ///
-
-private
-{
- static string generateLoopCode(string formatString, int N)() pure nothrow
- {
- string result;
- for (int i = 0; i < N; ++i)
- {
- string index = ctIntToString(i);
- // replace all @ by indices
- result ~= formatString.replace("@", index);
- }
- return result;
- }
-
- // Speed-up CTFE conversions
- static string ctIntToString(int n) pure nothrow
- {
- static immutable string[16] table = ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"];
- if (n < 10)
- return table[n];
- else
- return to!string(n);
- }
-}
-
-
-/// Element-wise minimum.
-@nogc Vector!(T, N) minByElem(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
-{
- import std.algorithm: min;
- Vector!(T, N) res = void;
- mixin(generateLoopCode!("res.v[@] = min(a.v[@], b.v[@]);", N)());
- return res;
-}
-
-/// Element-wise maximum.
-@nogc Vector!(T, N) maxByElem(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
-{
- import std.algorithm: max;
- Vector!(T, N) res = void;
- mixin(generateLoopCode!("res.v[@] = max(a.v[@], b.v[@]);", N)());
- return res;
-}
-
-/// Element-wise absolute value.
-@nogc Vector!(T, N) absByElem(T, int N)(const Vector!(T, N) a) pure nothrow
-{
- Vector!(T, N) res = void;
- mixin(generateLoopCode!("res.v[@] = abs(a.v[@]);", N)());
- return res;
-}
-
-/// Dot product of two vectors
-/// Returns: Dot product.
-@nogc T dot(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
-{
- T sum = 0;
- mixin(generateLoopCode!("sum += a.v[@] * b.v[@];", N)());
- return sum;
-}
-
-/// Cross product of two 3D vectors
-/// Returns: 3D cross product.
-/// Thanks to vuaru for corrections.
-@nogc Vector!(T, 3) cross(T)(const Vector!(T, 3) a, const Vector!(T, 3) b) pure nothrow
-{
- return Vector!(T, 3)(a.y * b.z - a.z * b.y,
- a.z * b.x - a.x * b.z,
- a.x * b.y - a.y * b.x);
-}
-
-/// 3D reflect, like the GLSL function.
-/// Returns: a reflected by normal b.
-@nogc Vector!(T, N) reflect(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
-{
- return a - (2 * dot(b, a)) * b;
-}
-
-///
-@nogc unittest
-{
- // reflect a 2D vector across the x axis (the normal points along the y axis)
- assert(vec2f(1,1).reflect(vec2f(0,1)) == vec2f(1,-1));
- assert(vec2f(1,1).reflect(vec2f(0,-1)) == vec2f(1,-1));
-
- // note that the normal must be, well, normalized:
- assert(vec2f(1,1).reflect(vec2f(0,20)) != vec2f(1,-1));
-
- // think of this like a ball hitting a flat floor at an angle.
- // the x and y components remain unchanged, and the z inverts
- assert(vec3f(2,3,-0.5).reflect(vec3f(0,0,1)) == vec3f(2,3,0.5));
-}
-
-/// Angle between two vectors
-/// Returns: angle between vectors.
-/// See_also: "The Right Way to Calculate Stuff" at $(WEB www.plunk.org/~hatch/rightway.php)
-@nogc T angleBetween(T, int N)(const Vector!(T, N) a, const Vector!(T, N) b) pure nothrow
-{
- auto aN = a.normalized();
- auto bN = b.normalized();
- auto dp = dot(aN, bN);
-
- if (dp < 0)
- return T(PI) - 2 * asin((-bN-aN).magnitude / 2);
- else
- return 2 * asin((bN-aN).magnitude / 2);
-}
-
-static assert(vec2f.sizeof == 8);
-static assert(vec3d.sizeof == 24);
-static assert(vec4i.sizeof == 16);
-
-unittest
-{
- static assert(vec2i.isValidSwizzle!"xyx");
- static assert(!vec2i.isValidSwizzle!"xyz");
- static assert(vec4i.isValidSwizzle!"brra");
- static assert(!vec4i.isValidSwizzle!"rgyz");
- static assert(vec2i.isValidSwizzleUnique!"xy");
- static assert(vec2i.isValidSwizzleUnique!"yx");
- static assert(!vec2i.isValidSwizzleUnique!"xx");
-
- alias vec2l = vec2!long;
- alias vec3ui = vec3!uint;
- alias vec4ub = vec4!ubyte;
-
- assert(vec2l(0, 1) == vec2i(0, 1));
-
- int[2] arr = [0, 1];
- int[] arr2 = new int[2];
- arr2[] = arr[];
- vec2i a = vec2i([0, 1]);
- vec2i a2 = vec2i(0, 1);
- immutable vec2i b = vec2i(0);
- assert(b[0] == 0 && b[1] == 0);
- vec2i c = arr;
- vec2l d = arr2;
- assert(a == a2);
- assert(a == c);
- assert(vec2l(a) == vec2l(a));
- assert(vec2l(a) == d);
-
- int[vec2i] hashMap;
- hashMap[a] = (c - a).squaredMagnitude;
- assert(hashMap[a] == (c - a).squaredMagnitude);
-
- vec4i x = [4, 5, 6, 7];
- assert(x == x);
- --x[0];
- assert(x[0] == 3);
- ++x[0];
- assert(x[0] == 4);
- x[1] &= 1;
- x[2] = 77 + x[2];
- x[3] += 3;
- assert(x == [4, 1, 83, 10]);
- assert(x.xxywz == [4, 4, 1, 10, 83]);
- assert(x.xxxxxxx == [4, 4, 4, 4, 4, 4, 4]);
- assert(x.abgr == [10, 83, 1, 4]);
- assert(a != b);
- x = vec4i(x.xyz, 166);
- assert(x == [4, 1, 83, 166]);
-
- vec2l e = a;
- vec2l f = a + b;
- assert(f == vec2l(a));
-
- vec3ui g = vec3i(78,9,4);
- g ^= vec3i(78,9,4);
- assert(g == vec3ui(0));
- //g[0..2] = 1u;
- //assert(g == [2, 1, 0]);
-
- assert(vec2i(4, 5) + 1 == vec2i(5,6));
- assert(vec2i(4, 5) - 1 == vec2i(3,4));
- assert(1 + vec2i(4, 5) == vec2i(5,6));
- assert(vec3f(1,1,1) * 0 == 0);
- assert(1.0 * vec3d(4,5,6) == vec3f(4,5.0f,6.0));
-
- auto dx = vec2i(1,2);
- auto dy = vec2i(4,5);
- auto dp = dot(dx, dy);
- assert(dp == 14 );
-
- vec3i h = cast(vec3i)(vec3d(0.5, 1.1, -2.2));
- assert(h == [0, 1, -2]);
- assert(h[] == [0, 1, -2]);
- assert(h[1..3] == [1, -2]);
- assert(h.zyx == [-2, 1, 0]);
-
- h.yx = vec2i(5, 2); // swizzle assignment
-
- assert(h.xy == [2, 5]);
- assert(-h[1] == -5);
- assert(++h[0] == 3);
-
- //assert(h == [-2, 1, 0]);
- assert(!__traits(compiles, h.xx = h.yy));
- vec4ub j;
-
- assert(lerp(vec2f(-10, -1), vec2f(10, 1), 0.5) == vec2f(0, 0));
-
- // larger vectors
- alias Vector!(float, 5) vec5f;
- vec5f l = vec5f(1, 2.0f, 3.0, 4u, 5.0L);
- l = vec5f(l.xyz, vec2i(1, 2));
-
- // the ctor should not compile if given too many arguments
- static assert(!is(typeof(vec2f(1, 2, 3))));
- static assert(!is(typeof(vec2f(vec2f(1, 2), 3))));
- static assert( is(typeof(vec3f(vec2f(1, 2), 3))));
- static assert( is(typeof(vec3f(1, 2, 3))));
-
- assert(absByElem(vec3i(-1, 0, 2)) == vec3i(1, 0, 2));
-}
-
diff --git a/src/shared/math.d b/src/shared/math.d
index 9765496..f678709 100644
--- a/src/shared/math.d
+++ b/src/shared/math.d
@@ -542,6 +542,8 @@ struct Matrix(T, int D)
auto l = &this;
auto r = &x;
+
+ // TODO: Test assembly on windows, may need to push registers onto the stack, decrement by 192 and move registers 1 by 1
asm pure @trusted @nogc nothrow
{
mov R8, l;