right-handed DirectionVector

Added missing Vector3::One

CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.18...3.28)
+cmake_minimum_required(VERSION 3.18...3.27)
 PROJECT(J3ML
         VERSION 1.1
         LANGUAGES CXX
@@ -34,7 +34,7 @@ set_target_properties(J3ML PROPERTIES LINKER_LANGUAGE CXX)
 
 CPMAddPackage(
         NAME jtest
-        URL https://git.redacted.cc/josh/jtest/archive/Prerelease-5.zip
+        URL https://git.redacted.cc/josh/jtest/archive/Release-1.4.zip
 )
 
 target_include_directories(J3ML PUBLIC ${jtest_SOURCE_DIR}/include)
@@ -58,4 +58,4 @@ target_link_libraries(MathDemo ${PROJECT_NAME})
 
 if(WIN32)
     #target_compile_options(MathDemo PRIVATE -mwindows)
-endif()
+endif()

Doxyfile

@@ -48,13 +48,13 @@ PROJECT_NAME           = J3ML
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         =
+PROJECT_NUMBER         = 3.1
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
 # quick idea about the purpose of the project. Keep the description short.
 
-PROJECT_BRIEF          = "Linear Algebra, Geometry, and Algorithms in C++"
+PROJECT_BRIEF          = "Linear Algebra, Geometry, and Algorithms in C++ (v3.1)"
 
 # With the PROJECT_LOGO tag one can specify a logo or an icon that is included
 # in the documentation. The maximum height of the logo should not exceed 55

README.md

@@ -4,7 +4,7 @@
 
 Yet Another C++ Math Standard
 
-J3ML is a "Modern C++" C++ library designed to provide comprehensive support for 3D mathematical operations commonly used in computer graphics, game development, physics simulations, and related fields. It offers a wide range of functionalities to simplify the implementation of complex mathematical operations in your projects.
+J3ML is a "Modern C++" library designed to provide comprehensive support for 3D mathematical operations commonly used in computer graphics, game development, physics simulations, and related fields. It offers a wide range of functionalities to simplify the implementation of complex mathematical operations in your projects.
 
 ![Static Badge](https://img.shields.io/badge/Lit-Based-%20)
 
@@ -50,7 +50,7 @@ Documentation is automatically generated from latest commit and is hosted at htt
 
 # Contributing
 
-Contributions to J3ML are welcome! If you find a bug, have a feature request, or would like to contribute code, please submit an issue or pull request to the GitHub repository.
+Contributions to J3ML are welcome! If you find a bug, have a feature request, or would like to contribute code, please submit an issue or pull request to the repository.
 
 # License
 

include/J3ML/Algorithm/Bezier.hpp

@@ -29,17 +29,61 @@ namespace J3ML::Algorithm
     template <typename T>
     inline T Cube(T f) { return f * f * f; }
 
+    /// Four points P0, P1, P2, P3 in the plane space define a cubic Bezier curve.
+    /// The curve can be modeled as a polynomial of third order.
+    /// The curve starts at P0, going toward P1, and arrives at P3 coming from the direction of P2.
+    /// Usually, it will not pass through P1, or P2, these points are only there to provide directional information.
     template <typename T>
     inline T Bezier(float t, const T& p0, const T& p1, const T& p2, const T& p3)
     {
         return Cube(1 - t) * p0 + 3 * Square(1 - t) * t * p1 + 3 * (1 - t) * Square(t) * p2 + Cube(t) * p3;
     }
 
+    /// Computes a point along a 2-dimensional Cubic Bezier curve.
+    /// @param t The normalized distance along the curve to compute, with range of [0, 1].
+    /// @param p0 The start-point of the curve.
+    /// @param p1
+    /// @param p2
+    /// @param p3 The end-point of the curve.
     Vector2 Bezier(float t, const Vector2& p0, const Vector2& p1, const Vector2& p2, const Vector2& p3);
 
-    // Tangent
+    /// Computes a point along the tangent of a 2-dimensional Cubic Bezier Curve.
+    /// @param t The normalized distance along the curve to compute, with range of [0, 1].
+    /// @param p0 The start-point of the curve.
+    /// @param p1
+    /// @param p2
+    /// @param p3 The end-point of the curve.
     Vector2 BezierDerivative(float t, const Vector2& p0, const Vector2& p1, const Vector2& p2, const Vector2& p3);
 
-    // Normal
+    /// Computes a point along the normal of a 2-dimensional Cubic Bezier Curve.
+    /// @param t The normalized distance along the curve to compute, with range of [0, 1].
+    /// @param p0 The start-point of the curve.
+    /// @param p1 The first control point,  which determines the direction at which the curve meets point 0.
+    /// @param p2 The second control point, which determines the direction at which the curve meets point 3.
+    /// @param p3 The end-point of the curve.
     Vector2 BezierNormal(float t, const Vector2& p0, const Vector2& p1, const Vector2& p2, const Vector2& p3);
+
+    /// Computes a point along a 3-dimensional Cubic Bezier curve.
+    /// @param t The normalized distance along the curve to compute, with range of [0, 1].
+    /// @param p0 The start-point of the curve.
+    /// @param p1 The first control point,  which determines the direction at which the curve meets point 0.
+    /// @param p2 The second control point, which determines the direction at which the curve meets point 3.
+    /// @param p3 The end-point of the curve.
+    Vector3 Bezier(float t, const Vector3& p0, const Vector3& p1, const Vector3& p2, const Vector3& p3);
+
+    /// Computes a point along the tangent of a 3-dimensional Cubic Bezier Curve.
+    /// @param t The normalized distance along the curve to compute, with range of [0, 1].
+    /// @param p0 The start-point of the curve.
+    /// @param p1 The first control point,  which determines the direction at which the curve meets point 0.
+    /// @param p2 The second control point, which determines the direction at which the curve meets point 3.
+    /// @param p3 The end-point of the curve.
+    Vector3 BezierDerivative(float t, const Vector3& p0, const Vector3& p1, const Vector3& p2, const Vector3& p3);
+
+    /// Computes a point along the normal of a 3-dimensional Cubic Bezier Curve.
+    /// @param t The normalized distance along the curve to compute, with range of [0, 1].
+    /// @param p0 The start-point of the curve.
+    /// @param p1 The first control point,  which determines the direction at which the curve meets point 0.
+    /// @param p2 The second control point, which determines the direction at which the curve meets point 3.
+    /// @param p3 The end-point of the curve.
+    Vector3 BezierNormal(float t, const Vector3& p0, const Vector3& p1, const Vector3& p2, const Vector3& p3);
 }

include/J3ML/Algorithm/GJK.hpp

@@ -1,4 +1,4 @@
-// @file GJK.hpp
+/// @file GJK.hpp
 /// Implementation of the Gilbert-Johnson-Keerthi (GJK) convex polyhedron intersection test
 
 #pragma once

include/J3ML/Algorithm/Parabola.hpp

@@ -0,0 +1,17 @@
+/// Josh's 3D Math Library
+/// A C++20 Library for 3D Math, Computer Graphics, and Scientific Computing.
+/// Developed and Maintained by Josh O'Leary @ Redacted Software.
+/// Special Thanks to William Tomasine II and Maxine Hayes.
+/// (c) 2024 Redacted Software
+/// This work is dedicated to the public domain.
+
+/// @file Parabola.hpp
+/// @desc Algorithm for calculating a parabolic curve to be used in instantaneous bullet raycasting.
+/// @edit 2024-10-22
+
+#pragma once
+
+namespace J3ML
+{
+
+}

include/J3ML/Algorithm/Triangulate.hpp

@@ -0,0 +1,18 @@
+/// Josh's 3D Math Library
+/// A C++20 Library for 3D Math, Computer Graphics, and Scientific Computing.
+/// Developed and Maintained by Josh O'Leary @ Redacted Software.
+/// Special Thanks to William Tomasine II and Maxine Hayes.
+/// (c) 2024 Redacted Software
+/// This work is dedicated to the public domain.
+
+/// @file Parabola.hpp
+/// @desc Algorithm for calculating a parabolic curve to be used in instantaneous bullet raycasting.
+/// @edit 2024-10-22
+
+#pragma once
+
+
+namespace J3ML
+{
+    /// @see class Polygon::Triangulate for current implementation.
+}

include/J3ML/Geometry.hpp

@@ -14,7 +14,6 @@
 #include <J3ML/Geometry/Polyhedron.hpp>
 #include <J3ML/Geometry/QuadTree.hpp>
 #include <J3ML/Geometry/Ray.hpp>
-#include <J3ML/Geometry/Shape.hpp>
 #include <J3ML/Geometry/Sphere.hpp>
 #include <J3ML/Geometry/Triangle.hpp>
 #include <J3ML/Geometry/Triangle2D.hpp>

include/J3ML/Geometry/AABB.hpp

@@ -14,12 +14,9 @@
 #include <format>
 #include <optional>
 
-#include <J3ML/LinearAlgebra/Common.hpp>
-#include <J3ML/Geometry/Common.hpp>
 #include <J3ML/Geometry/Shape.hpp>
-
-#include "Polygon.hpp"
-#include "Sphere.hpp"
+#include <J3ML/Geometry/Forward.hpp>
+#include <J3ML/LinearAlgebra.hpp>
 #include <J3ML/Algorithm/RNG.hpp>
 
 

include/J3ML/Geometry/AABB2D.hpp

@@ -12,8 +12,18 @@
 
 #pragma once
 
-#include <J3ML/LinearAlgebra/Vector2.hpp>
-#include "Shape.hpp"
+#include <J3ML/LinearAlgebra.hpp>
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
+
+
+/// See Christer Ericson's Real-time Collision Detection, p. 87, or
+/// James Arvo's "Transforming Axis-aligned Bounding Boxes" in Graphics Gems 1, pp. 548-550.
+/// http://www.graphicsgems.org/
+template <typename Matrix>
+void AABB2DTransformAsAABB2D(AABB2D& aabb, Matrix& m);
+
+
 
 namespace J3ML::Geometry
 {
@@ -30,28 +40,30 @@ namespace J3ML::Geometry
         minPoint(min), maxPoint(max)
         {}
 
-        float Width() const;
-        float Height() const;
+        [[nodiscard]] float Width() const;
+        [[nodiscard]] float Height() const;
 
-        float DistanceSq(const Vector2& pt) const;
+        Vector2 Centroid();
+
+        [[nodiscard]] float DistanceSq(const Vector2& pt) const;
 
         void SetNegativeInfinity();
 
         void Enclose(const Vector2& point);
 
-        bool Intersects(const AABB2D& rhs) const;
+        [[nodiscard]] bool Intersects(const AABB2D& rhs) const;
 
-        bool Contains(const AABB2D& rhs) const;
+        [[nodiscard]] bool Contains(const AABB2D& rhs) const;
 
-        bool Contains(const Vector2& pt) const;
+        [[nodiscard]] bool Contains(const Vector2& pt) const;
 
-        bool Contains(int x, int y) const;
+        [[nodiscard]] bool Contains(int x, int y) const;
 
-        bool IsDegenerate() const;
+        [[nodiscard]] bool IsDegenerate() const;
 
-        bool HasNegativeVolume() const;
+        [[nodiscard]] bool HasNegativeVolume() const;
 
-        bool IsFinite() const;
+        [[nodiscard]] bool IsFinite() const;
 
         Vector2 PosInside(const Vector2 &normalizedPos) const;
 
@@ -61,5 +73,30 @@ namespace J3ML::Geometry
 
         AABB2D operator-(const Vector2& pt) const;
 
+        Vector2 CornerPoint(int cornerIndex);
+
+        void TransformAsAABB(const Matrix3x3& transform);
+        void TransformAsAABB(const Matrix4x4& transform);
+
+
     };
+
+    template<typename Matrix>
+    void AABB2DTransformAsAABB2D(AABB2D &aabb, Matrix &m) {
+        float ax = m[0][0] * aabb.minPoint.x;
+        float bx = m[0][0] * aabb.maxPoint.x;
+        float ay = m[0][1] * aabb.minPoint.y;
+        float by = m[0][1] * aabb.maxPoint.y;
+
+        float ax2 = m[1][0] * aabb.minPoint.x;
+        float bx2 = m[1][0] * aabb.maxPoint.x;
+        float ay2 = m[1][1] * aabb.minPoint.y;
+        float by2 = m[1][1] * aabb.maxPoint.y;
+
+        aabb.minPoint.x = J3ML::Math::Min(ax, bx) + J3ML::Math::Min(ay, by) + m[0][3];
+        aabb.maxPoint.x = J3ML::Math::Max(ax, bx) + J3ML::Math::Max(ay, by) + m[0][3];
+
+        aabb.minPoint.y = J3ML::Math::Min(ax2, bx2) + J3ML::Math::Min(ay2, by2) + m[1][3];
+        aabb.maxPoint.y = J3ML::Math::Max(ax2, bx2) + J3ML::Math::Max(ay2, by2) + m[1][3];
+    }
 }

include/J3ML/Geometry/Capsule.hpp

@@ -11,10 +11,10 @@
 
 #pragma once
 
-#include <J3ML/LinearAlgebra/Common.hpp>
-#include <J3ML/Geometry/Common.hpp>
+#include <J3ML/LinearAlgebra.hpp>
 #include <J3ML/Geometry/Shape.hpp>
-#include "Circle.hpp"
+#include <J3ML/Geometry/Forward.hpp>
+#include <J3ML/Geometry/LineSegment.hpp>
 
 namespace J3ML::Geometry
 {
@@ -54,7 +54,7 @@ namespace J3ML::Geometry
 
 
         /// Quickly returns an arbitrary point inside this Capsule. Used in GJK intersection test.
-        [[nodiscard]] inline Vector3 AnyPointFast() const;
+        [[nodiscard]] Vector3 AnyPointFast() const;
 
         /// Generates a point that perhaps lies inside this capsule.
         /** @param height A normalized value between [0,1]. This specifies the point position along the height line of this capsule.

include/J3ML/Geometry/Circle.hpp

@@ -11,6 +11,8 @@
 
 #pragma once
 
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
 #include <J3ML/LinearAlgebra.hpp>
 
 namespace J3ML
@@ -18,7 +20,7 @@ namespace J3ML
 
     /// A two-dimensional circle in 3D space.
     /// This class represents both a hollow circle (only edge) and a solid circle (disc).
-    class Circle
+    class Circle : public Shape
     {
     public:
         /// The center position of this circle.

include/J3ML/Geometry/Forward.hpp

@@ -4,25 +4,25 @@
 // Forward declarations for classes that include each other
 namespace J3ML::Geometry
 {
-    class Shape;
     class AABB2D;
     class AABB;
     class Capsule;
     class Frustum;
     class LineSegment;
+    class LineSegment2D;
     class Line;
     class OBB;
+    class OBB2D;
     class Plane;
     class Polygon;
     class Polyhedron;
     template<typename T> class QuadTree;
     class Ray;
-    class Shape;
     class Sphere;
+    class Circle;
     class Triangle;
     class Triangle2D;
     class TriangleMesh;
-
     template <int N> class PBVolume;
 
 }

include/J3ML/Geometry/Frustum.hpp

@@ -3,9 +3,8 @@
 //
 #pragma once
 
-#include <J3ML/Geometry/Common.hpp>
-#include "Plane.hpp"
-#include "Shape.hpp"
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
 #include <J3ML/LinearAlgebra.hpp>
 
 namespace J3ML::Geometry
@@ -84,7 +83,7 @@ namespace J3ML::Geometry
         /// Specifies whether this frustum is a perspective or an orthographic frustum.
         FrustumType type;
         /// Specifies whether the [-1, 1] or [0, 1] range is used for the post-projective depth range.
-        FrustumProjectiveSpace projectiveSpace;
+        FrustumProjectiveSpace projectiveSpace ;
         /// Specifies the chirality of world and view spaces
         FrustumHandedness handedness;
         /// The eye point of this frustum

include/J3ML/Geometry/Icosahedron.hpp

@@ -0,0 +1,109 @@
+/// Josh's 3D Math Library
+/// A C++20 Library for 3D Math, Computer Graphics, and Scientific Computing.
+/// Developed and Maintained by Josh O'Leary @ Redacted Software.
+/// Special Thanks to William Tomasine II and Maxine Hayes.
+/// (c) 2024 Redacted Software
+/// This work is dedicated to the public domain.
+
+/// @file Icosahedron.hpp
+/// @desc Icosahedron class implementation, borrowed from http://www.songho.ca/opengl/gl_sphere.html#icosphere
+/// @edit 2024-10-22
+
+/** Polyhedron with 12 vertices, 30 edges, and 20 faces (triangles) for OpenGL
+    If the radius is r, then the length of the edge is (r / sin(2pi/5))
+
+    Vertices of icosahedron are constructed with spherical coords by aligning
+    the north pole to (0,0,r) and the south pole to (0,0,-r). Other 10 vertices
+    are computed by rotating 72 degrees along y-axis at the elevation angle
+    +/- 26.565 (=arctan(1/2)).
+
+    The unwrapped (paper model) of icosahedron and texture map look like this:
+    // (S,0)  3S  5S  7S  9S
+    //    /\  /\  /\  /\  /\      : 1st row (5 triangles)       //
+    //   /__\/__\/__\/__\/__\                                   //
+    // T \  /\  /\  /\  /\  /\    : 2nd row (10 triangles)      //
+    //    \/__\/__\/__\/__\/__\                                 //
+    // 2T  \  /\  /\  /\  /\  /   : 3rd row (5 triangles)       //
+    //      \/  \/  \/  \/  \/                                  //
+    //      2S  4S  6S  8S  (10S,3T)
+    // where S = 186/2048 = 0.0908203
+    //       T = 322/1024 = 0.3144531
+    // If a texture size is 2048x1024, S=186, T=322
+
+    AUTHOR: Song Ho Ahn (song.ahn@gmail.com)
+ */
+
+#include <vector>
+#include "Color4.hpp"
+
+#pragma once
+
+namespace J3ML
+{
+
+    class Icosahedron
+    {
+    public:
+        float radius;
+        float edgeLength;
+        Icosahedron(float radius = 1.0f);
+        float Radius() const { return radius; }
+        void Radius(float new_radius) {radius = new_radius;}
+        float EdgeLength() const { return edgeLength;}
+        void EdgeLength(float new_edge_length) { edgeLength = new_edge_length;}
+        // for vertex data
+        unsigned int VertexCount() const;
+        unsigned int NormalCount() const;
+        unsigned int TexCoordCount() const;
+        unsigned int IndexCount() const;
+        unsigned int LineIndexCount() const;
+        unsigned int TriangleCount() const;
+        unsigned int VertexSize() const;
+        unsigned int NormalSize() const;
+        unsigned int TexCoordSize() const;
+        unsigned int IndexSize() const;
+        unsigned int LineIndexSize() const;
+        const float* Vertices() const;
+        const float* Normals() const;
+        const float* TexCoords() const;
+        const unsigned int* Indices() const;
+        const unsigned int* LineIndices() const;
+        // for interleaved vertices: V/N/T
+        unsigned int InterleavedVertexCount() const;
+        unsigned int InterleavedVertexSize() const;
+        int InterleavedStride() const;
+        const float* InterleavedVertices() const;
+
+        // draw in VertexArray mode
+        void draw() const;
+        void drawLines(const Color4& lineColor) const;
+        void drawWithLines(const Color4& lineColor) const;
+
+    protected:
+    private:
+        // static functions
+        static void computeFaceNormal(float v1[3], float v2[3], float v3[3], float n[3]);
+
+        // member functions
+        void updateRadius();
+        std::vector<float> computeVertices();
+        void buildVertices();
+        void buildInterleavedVertices();
+        void addVertices(float v1[3], float v2[3], float v3[3]);
+        void addNormals(float n1[3], float n2[3], float n3[3]);
+        void addTexCoords(float t1[2], float t2[2], float t3[2]);
+        void addIndices(unsigned int i1, unsigned int i2, unsigned int i3);
+        void addLineIndices(unsigned int indexFrom);
+
+        // member vars
+        //float radius;
+        //float edgeLength;
+        std::vector<float> vertices;
+        std::vector<float> normals;
+        std::vector<float> texCoords;
+        std::vector<unsigned int> indices;
+        std::vector<unsigned int> lineIndices;
+
+
+    };
+}

include/J3ML/Geometry/Line.hpp

@@ -1,19 +1,20 @@
 #pragma once
 
 
-#include "LineSegment.hpp"
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/LineSegment.hpp>
 
 namespace J3ML::Geometry
 {
-    class Line
+    class Line : public Shape
     {
     public:
         Vector3 Position;  /// Specifies the origin of this line.
         Vector3 Direction; /// The normalized direction vector of this ray.
-    public:
-        static void
-        ClosestPointLineLine(const Vector3 &v0, const Vector3 &v10, const Vector3 &v2, const Vector3 &v32, float &d,
-                             float &d2);
+
+        Line(const Vector3& position, const Vector3& direction) : Position(position), Direction(direction) {}
+
+        static void ClosestPointLineLine(const Vector3 &v0, const Vector3 &v10, const Vector3 &v2, const Vector3 &v32, float &d, float &d2);
 
         Vector3 ClosestPoint(const J3ML::Geometry::LineSegment &other, float &d, float &d2) const;
 

include/J3ML/Geometry/LineSegment.hpp

@@ -1,16 +1,14 @@
 #pragma once
 
-#include <J3ML/LinearAlgebra/Common.hpp>
-#include <J3ML/Geometry/Common.hpp>
-#include <J3ML/LinearAlgebra/Vector3.hpp>
+#include <J3ML/LinearAlgebra.hpp>
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
 
 namespace J3ML::Geometry
 {
 
-    using LinearAlgebra::Vector3;
-
     /// A line segment in 3D space is a finite line with a start and end point.
-    class LineSegment
+    class LineSegment : public Shape2D
     {
     public:
         Vector3 A; /// The starting point of this line segment.

include/J3ML/Geometry/LineSegment2D.hpp

@@ -0,0 +1,185 @@
+/// Josh's 3D Math Library
+/// A C++20 Library for 3D Math, Computer Graphics, and Scientific Computing.
+/// Developed and Maintained by Josh O'Leary @ Redacted Software.
+/// Special Thanks to William Tomasine II and Maxine Hayes.
+/// (c) 2024 Redacted Software
+/// This work is dedicated to the public domain.
+
+/// @file LineSegment2D.hpp
+/// @desc A 2D representation of a finite line between two points.
+/// @edit 2024-10-22
+
+#pragma once
+
+#include <J3ML/LinearAlgebra/Vector2.hpp>
+#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
+#include <J3ML/LinearAlgebra/Matrix4x4.hpp>
+
+namespace J3ML::Geometry
+{
+    /// A line segment in 2D space is a finite line with a start and end point.
+    class LineSegment2D {
+    public:
+        /// The starting point of this line segment.
+        Vector2 A;
+        /// The end point of this line segment.
+        Vector2 B;
+
+        /// The default constructor does not initialize any members of this class.
+        /** This means that the values of the members a and b are undefined after creating a new LineSegment2D using this
+            default constructor. Remember to assign to them before use.
+            @see A, B. */
+        LineSegment2D() {}
+
+        /// Constructs a line segment through the given end points.
+        LineSegment2D(const Vector2 &a, const Vector2 &b);
+
+        /// Returns a point on the line.
+        /** @param d The normalized distance along the line segment to compute. If a value in the range [0, 1] is passed, then the
+                returned point lies along this line segment. If some other value is specified, the returned point lies on the line
+                defined by this line segment, but *not* inside the interval from a to b.
+            @note The meaning of d here differs from Line2D::GetPoint and Ray2D::GetPoint. For the class LineSegment2D,
+                GetPoint(0) returns a, and getPoint(1) returns b. This means that GetPoint(1) will not generally be exactly one unit
+                away from the starting point of this line segment, as is the case with Line2D and Ray2D.
+            @return (1-d)*a + d*b;
+            @see a, b, Line2D::GetPoint(), Ray2D::GetPoint() */
+        Vector2 GetPoint(float d) const;
+
+        /// Returns the center point of this line segment.
+        /** This function is the same as calling GetPoint(0.5f), but provided here as convenience.
+            @see GetPoint(). */
+        Vector2 CenterPoint() const;
+
+        Vector2 Centroid() const;
+
+        /// Reverses the direction of this line segment.
+        /** This function swaps the start and end points of this line segment so that it runs from b to a.
+            This does not have an effect on the set of points represented by this line segment, but it reverses
+            the direction of the vector returned by Dir().
+            @note This function operates in-place.
+            @see a, b, Dir(). */
+        void Reverse();
+
+        /// Returns the normalized direction vector that points in the direction a->b.
+        /** @note The returned vector is normalized, meaning that its length is 1, not |b-a|.
+            @see a, b. */
+        Vector2 Dir() const;
+
+        /// Quickly returns an arbitrary point inside this LineSegment2D. Used in GJK intersection test.
+        Vector2 AnyPointFast() const;
+
+        /// Computes an extreme point of this LineSegment2D in the given direction.
+        /** An extreme point is a farthest point along this LineSegment2D in the given direction. Given a direction,
+            this point is not necessarily unique.
+            @param direction The direction vector of the direction to find the extreme point. This vector may
+                be un-normalized, but may not be null.
+            @return An extreme point of this LineSegment2D in the given direction. The returned point is always
+                either a or b.
+            @see a, b. **/
+        Vector2 ExtremePoint(const Vector2 &direction) const;
+
+        Vector2 ExtremePoint(const Vector2 &direction, float &projectionDistance) const;
+
+        /// Translates this LineSegment2D in world space.
+        /** @param offset The amount of displacement to apply to this LineSegment2D, in world space coordinates.
+            @see Transform(). */
+        void Translate(const Vector2& offset);
+
+        /// Applies a transformation to this line.
+        /** This function operates in-place.
+            @see Translate(), Transform(), classes Matrix3x3, Matrix4x4, Quaternion*/
+        void Transform(const Matrix3x3&  transform);
+        void Transform(const Matrix4x4&  transform);
+        void Transform(const Quaternion& transform);
+
+        /// Computes the length of this line segment.
+        /** @return |b - a|
+            @see a, b. */
+        float Length() const;
+        /// Computes the squared length of this line segment.
+        /** Calling this function is faster than calling Length(), since this function avoids computing a square root.
+            If you only need to compare lengths to each other and are not interested in the actual length values,
+            you can compare using LengthSq(), instead of Length(), since Sqrt() is an order-preserving,
+            (monotonous and non-decreasing) function. */
+        float LengthSq() const;
+        float LengthSquared() const;
+
+        /// Tests if this line segment is finite
+        /** A line segment is finite if its endpoints a and b do not contain floating-point NaNs for +/- infs
+            in them.
+            @return True if both a and b have finite floating-point values. */
+        bool IsFinite() const;
+
+        /// Tests if this line segment represents the same set of points than the the given line segment.
+        /** @param epsilon Specifies how much distance threshold to allow in the comparison.
+            @return True if a == rhs.a && b == rhs.b, or, a == rhs.b && b == rhs.a, within the given epsilon. */
+        bool Equals(const LineSegment2D& rhs, float epsilon = 1e-3f) const;
+
+        /// Tests if the given point or line segment is contained on this line segment.
+        /** @param epsilon Because a line segment is a one-dimensional object in 3D space, an epsilon value
+                is used as a threshold for this test. This effectively transforms this line segment to a capsule with
+                the radius indicated by this value.
+            @return True if this line segment contains the given point or line segment.
+            @see Intersects, ClosestPoint(), Distance(). */
+        bool Contains(const Vector2& point, float epsilon = 1e-3f) const;
+        bool Contains(const LineSegment2D& rhs, float epsilon = 1e-3f) const;
+
+        /// Computes the closest point on this line segment to the given object.
+        /** @param d If specified, this parameter receives the normalized distance along
+                this line segment which specifies the closest point on this line segment to
+                the specified point.
+            @return The closest point on this line segment to the given object.
+            @see Contains(), Distance(), Intersects(). */
+        Vector2 ClosestPoint(const Vector2& point) const;
+        Vector2 ClosestPoint(const Vector2& point, float& d) const;
+
+
+
+        /** @param d2 [out] If specified, this parameter receives the (normalized, in case of line segment)
+                distance along the other line object which specifies the closest point on that line to
+                this line segment. */
+        Vector2 ClosestPoint(const LineSegment2D& other) const;
+        Vector2 ClosestPoint(const LineSegment2D& other, float& d) const;
+        Vector2 ClosestPoint(const LineSegment2D& other, float& d, float& d2) const;
+
+        /// Computes the distance between this line segment and the given object.
+        /** @param d [out] If specified, this parameter receives the normalized distance along
+                this line segment which specifies the closest point on this line segment to
+                the specified point.
+            @return The distance between this line segment and the given object.
+            @see */
+        float Distance(const Vector2& point) const;
+        float Distance(const Vector2& point, float& d) const;
+
+        /** @param d2 [out] If specified, this parameter receives the (normalized, in case of line segment)
+                distance along the other line object which specifies the closest point on that line to
+                this line segment. */
+        float Distance(const LineSegment2D& other) const;
+        float Distance(const LineSegment2D& other, float& d) const;
+        float Distance(const LineSegment2D& other, float& d, float& d2) const;
+
+        float DistanceSq(const Vector2& point) const;
+        float DistanceSq(const LineSegment2D& other) const;
+
+        /** @param epsilon If testing intersection between two line segments, a distance threshold value is used to account
+                for floating point inaccuracies. */
+        bool Intersects(const LineSegment2D& lineSegment, float epsilon = 1e-3f) const;
+
+
+        /// Projects this LineSegment2D onto the given 1D axis direction vector.
+        /** This function collapses this LineSegment2D onto a 1D axis for the purposes of e.g. separate axis test computations.
+            This function returns a 1D range [outMin, outNax] denoting the interval of the projection.
+            @param direction The 1D axis to project to. This vector may be unnormalized, in which case the output
+                of this function gets scaled by the length of this vector.
+            @param outMin [out] Returns the minimum extent of this vector along the projection axis.
+            @param outMax [out] Returns the maximum extent of this vector along the projection axis. */
+        void ProjectToAxis(const Vector2& direction, float& outMin, float& outMax) const;
+
+    protected:
+    private:
+    };
+
+    LineSegment2D operator * (const Matrix3x3& transform, const LineSegment2D& line);
+    LineSegment2D operator * (const Matrix4x4& transform, const LineSegment2D& line);
+    LineSegment2D operator * (const Quaternion& transform, const LineSegment2D& line);
+}

include/J3ML/Geometry/OBB.hpp

@@ -1,8 +1,8 @@
 #pragma once
 
-#include <J3ML/Geometry/Common.hpp>
-#include <J3ML/Geometry/AABB.hpp>
-#include <J3ML/Geometry/Polyhedron.hpp>
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
+#include <J3ML/LinearAlgebra.hpp>
 
 
 namespace J3ML::Geometry {

include/J3ML/Geometry/PBVolume.hpp

@@ -11,13 +11,11 @@
 
 #pragma once
 
-
 #include <J3ML/Geometry/AABB.hpp>
 #include <J3ML/Geometry/Polyhedron.hpp>
 #include <J3ML/Geometry/Plane.hpp>
+#include <J3ML/Geometry/Sphere.hpp>
 
-#include "Sphere.hpp"
-#include "Sphere.hpp"
 
 namespace J3ML::Geometry
 {

include/J3ML/Geometry/Plane.hpp

@@ -1,14 +1,14 @@
 #pragma once
-#include <J3ML/LinearAlgebra/Vector3.hpp>
-#include "Shape.hpp"
-#include "Ray.hpp"
+
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/LinearAlgebra.hpp>
+#include <J3ML/Geometry/Forward.hpp>
 
 
 namespace J3ML::Geometry
 {
-    using J3ML::LinearAlgebra::Vector3;
 
-    class Plane
+    class Plane : public Shape
     {
     public:
         Vector3 Position;

include/J3ML/Geometry/Polygon.hpp

@@ -11,14 +11,14 @@
 
 #pragma once
 
-#include <J3ML/Geometry/Common.hpp>
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
+#include <J3ML/LinearAlgebra.hpp>
 #include <vector>
-#include "Shape.hpp"
-#include <J3ML/LinearAlgebra/Common.hpp>
 
 namespace J3ML::Geometry {
 
-    class Polygon : public Shape
+    class Polygon : public Shape2D
     {
     public:
         std::vector<Vector3> vertices;

include/J3ML/Geometry/Polyhedron.hpp

@@ -1,17 +1,16 @@
 #pragma once
 
-#include <J3ML/Geometry/Common.hpp>
-#include <J3ML/Geometry/Shape.hpp>
 #include <vector>
-#include <J3ML/LinearAlgebra/Vector3.hpp>
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
+#include <J3ML/LinearAlgebra.hpp>
 
 
 namespace J3ML::Geometry
 {
-    using namespace J3ML::LinearAlgebra;
 
     // Represents a three-dimensional closed geometric solid defined by flat polygonal faces.
-    class Polyhedron
+    class Polyhedron : public Shape
     {
     public:
         // Stores a list of indices of a single face of a Polygon

include/J3ML/Geometry/QuadTree.hpp

@@ -2,14 +2,13 @@
 
 #include <vector>
 #include <cstdint>
-#include <J3ML/LinearAlgebra/Vector2.hpp>
+#include <J3ML/LinearAlgebra.hpp>
+#include <J3ML/Geometry/Forward.hpp>
+
 #include <J3ML/Geometry/AABB2D.hpp>
 
 namespace J3ML::Geometry {
 
-
-    using LinearAlgebra::Vector2;
-
     template<typename T>
     class QuadTree {
         /// A fixed split rule for all QuadTrees: A QuadTree leaf node is only ever split if the leaf contains at least this many objects.
@@ -366,4 +365,4 @@ namespace J3ML::Geometry {
             leaf->objects.pop_back();
         }
     }
-}
+}

include/J3ML/Geometry/Ray.hpp

@@ -10,11 +10,13 @@
 /// @edit 2024-07-06
 #pragma once
 
-#include <J3ML/LinearAlgebra/Vector3.hpp>
+#include <format>
+#include <J3ML/LinearAlgebra.hpp>
 #include <vector>
-#include "TriangleMesh.hpp"
-#include "Frustum.hpp"
-#include "OBB.hpp"
+#include <ostream>
+#include <iosfwd>
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
 
 namespace J3ML::Geometry
 {
@@ -29,12 +31,12 @@ namespace J3ML::Geometry
         Vector3 Intersection;
         Vector3 SurfaceNormal;
         bool Hit;
-        Shape* Target;
-        static RaycastResult NoHit() { return {Vector3::NaN, Vector3::NaN, false, nullptr};}
+        //Shape* Target;
+        static RaycastResult NoHit() { return {Vector3::NaN, Vector3::NaN, false};}
     };
 
     /// A ray in 3D space is a line that starts from an origin point and extends to infinity in one direction
-    class Ray
+    class Ray : public Shape
     {
     public: // Properties
         // The position of this ray.
@@ -151,17 +153,12 @@ namespace J3ML::Geometry
         // Returns a RaycastResult structure containing the first object the ray collides with,
         // the surface intersection point,
         // and the surface normal at the point of intersection.
-        RaycastResult Cast(std::vector<Shape> shapes, float maxDistance = 99999999);
+       // RaycastResult Cast(std::vector<Shape> shapes, float maxDistance = 99999999);
 
         void ProjectToAxis(const Vector3 &direction, float &outMin, float &outMax) const;
 
-        [[nodiscard]] std::string ToString() const {
-            return std::format("Ray(origin:[{}], direction:[{}])", Origin.ToString(), Direction.ToString());
-        }
+        [[nodiscard]] std::string ToString() const;
     };
 
-    inline std::ostream& operator << (std::ostream& o, const Ray& ray) {
-        o << ray.ToString();
-        return o;
-    }
-}
+    std::ostream& operator << (std::ostream& o, const Ray& ray);
+}

include/J3ML/Geometry/Shape.hpp

@@ -1,27 +1,8 @@
 #pragma once
 
-namespace J3ML::Geometry
-{
-    class GeometricPrimitive
-    {
-    public:
-    protected:
-    private:
-    };
+class Shape {
+public:
+    virtual ~Shape() = default;
+};
 
-    class Shape
-    {
-    public:
-        virtual ~Shape() = default; //Polymorphic for dynamic_cast.
-    protected:
-    private:
-    };
-
-    class Shape2D
-    {
-    public:
-        virtual ~Shape2D() = default;
-    protected:
-    private:
-    };
-}
+class Shape2D : public Shape{};

include/J3ML/Geometry/Sphere.hpp

@@ -11,21 +11,16 @@
 
 #pragma once
 
-
-#include <J3ML/Geometry.hpp>
-#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
-#include <J3ML/LinearAlgebra/Matrix4x4.hpp>
-#include <J3ML/Geometry/LineSegment.hpp>
-#include <J3ML/Geometry/TriangleMesh.hpp>
 #include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
+#include <J3ML/LinearAlgebra.hpp>
+
 
 namespace J3ML::Geometry
 {
-    using J3ML::LinearAlgebra::Matrix3x3;
-    using J3ML::LinearAlgebra::Matrix4x4;
 
     // A mathematical representation of a 3-dimensional sphere
-    class Sphere
+    class Sphere : public Shape
     {
     public: // Properties
         Vector3 Position; // The center point of this sphere.
@@ -119,8 +114,10 @@ namespace J3ML::Geometry
         [[nodiscard]] bool Contains(const Vector3& point, float epsilon) const;
 
         [[nodiscard]] bool Contains(const LineSegment& lineseg)  const;
-        TriangleMesh GenerateUVSphere() const {}
-        TriangleMesh GenerateIcososphere() const {}
+        TriangleMesh GenerateUVSphere(int subdivisions = 10.f) const;
+        TriangleMesh GenerateIcososphere() const;
+        TriangleMesh GenerateCubesphere() const;
+
 
         void ProjectToAxis(const Vector3 &direction, float &outMin, float &outMax) const;
     };

include/J3ML/Geometry/Triangle.hpp

@@ -1,13 +1,13 @@
 #pragma once
 
-#include "J3ML/LinearAlgebra/Vector3.hpp"
-#include <J3ML/Geometry/Common.hpp>
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
 #include <J3ML/LinearAlgebra.hpp>
 #include <cfloat>
 
 namespace J3ML::Geometry
 {
-    class Triangle
+    class Triangle : public Shape
     {
     public:
         Vector3 V0;
@@ -137,6 +137,7 @@ namespace J3ML::Geometry
         Vector3 Vertex(int i) const;
 
         LineSegment Edge(int i) const;
+        Triangle(const Vector3& v0, const Vector3& v1, const Vector3& v2) : V0(v0), V1(v1), V2(v2) {}
     };
 
 

include/J3ML/Geometry/Triangle2D.hpp

@@ -1,12 +1,242 @@
 #pragma once
 
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/LinearAlgebra/Forward.hpp>
+#include <J3ML/LinearAlgebra/Vector2.hpp>
+#include <J3ML/Geometry/Forward.hpp>
+#include <J3ML/Algorithm/RNG.hpp>
 
 namespace J3ML::Geometry
 {
-    class Triangle2D {
+    using J3ML::LinearAlgebra::Vector2;
+    using J3ML::Algorithm::RNG;
+    /// Specifies a triangle through three points in 2D space.
+    /// This class stores three member vertices A, B, and C to specify the triangle.
+    /** @note The order in which the vertices are stored in this data structure is important.
+              The triangles (a,b,c) and (a,c,b) are otherwise equivalent, but their plane normals point to the opposite directions.
+        @see PlaneCCW(), PlaneCW() */
+    class Triangle2D : public Shape2D {
+
     public:
-        Vector2 A;
-        Vector2 B;
-        Vector2 C;
+        Vector2 A; /// The first Triangle endpoint.
+        Vector2 B; /// The second Triangle endpoint.
+        Vector2 C; /// The third Triangle endpoint.
+    public:
+        /// The default constructor does not initialize any members of this class.
+        /** This means that the values of the members a, b, and c are undefined after creating a new Triangle2D using this
+            default constructor. Remember to assign to them before use. */
+        Triangle2D() {}
+        /// Constructs a triangle from three given endpoints.
+        /** The normal of the plane will be constructed to point towards the halfspace where
+            the vertices a, b, and c wind in counter-clockwise order. */
+        Triangle2D(const Vector2& a, const Vector2& b, const Vector2& c);
+    public:
+        static int NumFaces() { return 1; }
+        static int NumEdges() { return 3; }
+        static int NumVertices() { return 3; }
+
+        /// Translates this Triangle2D in world space.
+        /** @param offset The amount of displacement to apply to this Triangle2D, in world space coordinates.
+            @see Transform(). */
+        void Translate(const Vector2& offset);
+
+        /// Applies a transformation to this Triangle2D, in-place.
+        /** @see Translate(), classes Matrix3x3, Matrix4x4, Quaternion. */
+        void Transform(const Matrix3x3& transform);
+        void Transform(const Matrix4x4& transform);
+
+        /// Expresses the given point in barycentric (u,v,w) coordinates.
+        /** @note There are two different conventions for representing barycentric coordinates. One uses a
+                (u, v, w) triplet with the equation pt = u*a + v*b + w*c, and the other uses a (u,v) pair
+                with the equation pt = a + u*(b-a) + v*(c-a). These two are equivalent. Use the mappings
+                (u,v) -> (1-u-v, u, v) and (u,v,w)->(v,w) to convert between these two representations.
+            @param point The point of the vector space to express in barycentric coordinates. This point should
+                lie in the plane formed by this triangle.
+            @return The factors (u,v,w) that satisfy the weighted sum equation point == u*a + v*b + w*c.
+            @see BarycentricUV(), BarycentricInsideTriangle(), Point(), http://mathworld.wolfram.com/BarycentricCoordinates.html */
+        Vector3 BarycentricUVW(const Vector2& point) const;
+
+        /// Expresses the given point in barycentric (u, v) coordinates.
+        /** @note There are two different conventions for representing barycentric coordinates. One uses a
+                (u, v, w) triplet with the equation pt = u*a + v*b + w*c, and the other uses a (u,v) pair
+                with the equation pt = a + u*(b-a) + v*(c-a). These two are equivalent. Use the mappings
+                (u,v) -> (1-u-v, u, v) and (u,v,w)->(v,w) to convert between these two representations.
+            @param point The point to express in barycentric coordinates. This point should lie in the plane
+                formed by this triangle.
+            @return The factors (u, v) that satisfy the weighted sum equation point == a + u*(b-a) + v*(c-a)
+            @see BarycentricUVW(), BarycentricInsideTriangle(), Point() */
+        Vector2 BarycentricUV(const Vector2& point) const;
+
+        /// Tests if the given barycentric UVW coordinates lie inside a triangle.
+        /** A barycentric UVW coordinate represents a point inside a triangle if
+                a) 0 <= u,v,w <= 1 and
+                b) u+v+w == 1.
+            @param uvw The barycentric vector containing the barycentric (u,v,w) coordinates.
+            @return True if the given coordinates lie inside a triangle.
+            @see BarycentricUV(), BarycentricUVW(), Point(). */
+        static bool BarycentricInsideTriangle(const Vector3& uvw);
+
+        /// Returns the point at the given barycentric coordinates.
+        /** This function computes the vector space point at the given barycentric coordinates.
+            @param uvw The barycentric UVW coordiante triplet. The condition u+v+w == 1 should hold for the input coordinate.
+                If 0 <= u,v,w <= 1, the returned point lies inside this triangle.
+            @return u*a + v*b + w*c. */
+        Vector2 Point(const Vector3& uvw) const;
+        Vector2 Point(float u, float v, float w) const;
+
+        /** These functions are an alternate form of Point(u,v,w) for the case when the barycentric coordinates are
+            represented as a (u,v) pair and not as a (u,v,w) triplet. This function is provided for convenience
+            and effectively just computes Point(1-u-v, u, v).
+            @param uv The barycentric UV coordinates. If 0 <= u,v <= 1 and u+v <= 1, then the returned point lies inside
+                this triangle.
+            @return a + (b-a)*u + (c-a)*v.
+            @see BarycentricUV(), BarycentricUVW(), BarycentricInsideTriangle(). */
+        Vector2 Point(const Vector2& uv) const;
+        Vector2 Point(float u, float v) const;
+
+        /// Computes the center of mass of this triangle.
+        /** @return The point (a+b+c)/3. */
+        Vector2 Centroid() const;
+        /// Identical to Centroid(), we provide both to enable common signature with AABB2D, and OBB2D to allow them in template algorithms.
+        Vector2 CenterPoint() const { return Centroid();}
+        /// Computes the surface area of this triangle.
+        /** @return The surface area of this triangle.
+            @see Perimeter(), Area2D(), SignedArea(). */
+        float Area() const;
+
+        /// Computes the total edge length of this triangle.
+        /** @return |a-b| + |b-c| + |c-a|.
+            @see Area(), Edge(). */
+        float Perimeter() const;
+
+        /// Returns a pointer to the vertices of this triangle. The array contains three elements.
+        Vector2* VertexArrayPtr() { return &A;}
+        const Vector2* VertexArrayPtr() const { return &A;}
+
+        /// Returns a vertex of this triangle.
+        /** @param i The vertex of the triangle to get: 0, 1, or 2.
+            @return Vertex(0) returns the point A, Vertex(1) returns the point B, and Vertex(2) returns the point C.
+            @note If an index outside [0,2] is passed, an assume() failure occurs, and a NaN vector is returned.
+            @see Edge(). */
+        Vector2 Vertex(int i) const;
+        Vector2 CornerPoint(int i) const { return Vertex(i);} // An alias for Vertex() to be able to mesh Triangle2D into templatized algorithms.
+
+
+        //LineSegment2D Edge(int i) const; // TODO: Implement class LineSegment2D
+
+        /// Quickly returns an arbitrary point inside this Triangle2D. Used in GJK intersection test.
+        inline Vector2 AnyPointFast() const { return A;}
+
+        /// Computes an extreme point of this Triangle2D in the given direction.
+        /** An extreme point is the farthest point of this Triangle2D in the given direction. Given a direction,
+            this point is not necessarily unique.
+            @param direction THe direction vector of the direction to find the extreme point. This vector may
+                be unnormalized, but may not be null.
+            @return An extreme point of this Triangle2D in the given direction. The returned point is always a
+                vertex of this Triangle2D.
+            @see Vertex(). */
+        Vector2 ExtremePoint(const Vector2& direction) const;
+        Vector2 ExtremePoint(const Vector2& direction, float& projectionDistance) const;
+
+        /// Returns the tight AABB2D that encloses this Triangle2D.
+        AABB2D BoundingAABB() const;
+
+        /// Computes the surface area of the given 2D triangle.
+        /** This math library does not have a separate class for for 2D triangles. To compute the area of a 2D triangle,
+            use this Triangle2D class and set z=0 for each coordinate, or use this helper function.
+            @see Area(), SignedArea().*/
+        static float Area2D(const Vector2& p1, const Vector2& p2, const Vector2& p3);
+        /// Relates the barycentric coordinate of the given point to the surface area of triangle abc.
+        /** This function computes the ratio of the signed area of the triangle (point, b, c) to the signed area of
+            the triangle (a,b,c). This is the same as computng the barycentric u-coordinate of the given point
+            on the triangle (a,b,c).
+            @see Area(), Area2D(), BarycentricUVW(). */
+        static float SignedArea(const Vector2& point, const Vector2& a, const Vector2& b, const Vector2& c);
+
+        /// Returns true if this triangle is finite.
+        /** A triangle is <b><i>finite</i></b> if its vertices a,b, and c do not contain floating-point NaNs or +/-infs
+            in them.
+            @return True if each coordinate of each vertex of this triangle has a finite floating-point value.
+            @see A, B, C, IsDegenerate(), ::IsFinite(), IsInf(), IsNan(), IsFinite() */
+        bool IsFinite() const;
+        /// Returns true if this triangle is degenerate.
+        /** A triangle is <b><i>degenerate</i></b> if it is not finite, or if the surface area of the triangle is
+            close to zero.
+            @param epsilon The threshold to test against. If two vertices of this triangle are closer than this, the
+            triangle is considered degenerate.
+            @see a, b, c, IsFinite() */
+        bool IsDegenerate(float epsilon = 1e-3f) const;
+
+        /// Returns true if the triangle defined by the three given points is degenerate.
+        static bool IsDegenerate(const Vector2& p1, const Vector2& p2, const Vector2& p3, float epsilon = 1e-3f);
+
+        /// In some templated algorithms, the input can either be a Triangle2D or a Polygon2D. Provide trivial Polygon2D-specific API
+        /// for compatibility in those template functions.
+        bool IsConvex() const { return true;}
+        bool IsPlanar() const { return true;}
+        bool IsSimple() const { return true;}
+
+        /// Tests if the given object is fully contained inside this triangle.
+        /** @param triangleThickness An epsilon threshold value to use for this test. triangleThicknessSq is the squared version of this parameter.
+                This specifies the maximum distance the given object can lie from the plane defined by this triangle.
+            @see Distance(), Intersects(), ClosestPoint(). */
+        bool Contains(const Vector2& point, float triangleThicknessSq = 1e-5f) const;
+        /// Computes the distance between this triangle and the given object.
+        /** This function finds the nearest pair of points on this and the given object, and computes their distance.
+            If the two objects intersect, or one object is contained inside the other, the returned distance is zero.
+            @see Contains(), Intersects(), ClosestPoint(). */
+        float Distance(const Vector2& point) const;
+        float DistanceSquared(const Vector2& point) const;
+        float DistanceSq(const Vector2& point) const { return DistanceSquared(point); }
+
+        /// A helper function used in line-triangle tests.
+        static float IntersectLineTri(const Vector2& linePos, const Vector2& lineDir,
+                                      const Vector2& v0, const Vector2& v1, const Vector2& v2,
+                                      float& u, float& v);
+
+        /// Projects this triangle onto the given axis.
+        /** This function is used in SAT tests (separate axis theorem) to check the interval this triangle
+            lies in on an 1D line.
+            @param axis THe axis to project onto. This vector can be unnormalized.
+            @param dMin [out] Returns the minimum extent of this triangle on the given axis.
+            @param dMax [out] Returns the maximum extent of this triangle on the given axis. */
+        void ProjectToAxis(const Vector2& axis, float& dMin, float& dMax) const;
+        int UniqueFaceNormals(Vector2* out) const;
+        int UniqueEdgeNormals(Vector2* out) const;
+
+        /// Computes the closest point on this triangle to the given object.
+        /** If the other object intersects this triangle, the function will return an arbitrary point inside
+            the region of intersection.
+            @see Contains(), Distance(), Intersects(), ClosestPointToTriangleEdge(). */
+        Vector2 ClosestPoint(const Vector2& point) const;
+        //Vector2 ClosestPoint(const LineSegment2D& lineSegment, Vector2* otherPt = 0) const;
+
+        /// Generates a random point inside this Triangle2D.
+        /** The points are distrubited uniformly.
+            The implementation of this function is based on Graphics Gems 1, p. 25:
+            "1.5 Generating random points on triangles. Method 2." The Method 1 presented in the book
+            uses a Sqrt() instead of the if().
+            @param rng A pre-seeded random number generator object that is to be used by this function to generate random values.
+            @see class RNG, RandomPointOnEdge(), RandomVertex(), Point() */
+        Vector2 RandomPointInside(RNG& rng) const;
+
+        /// Choose a corner vertex of this triangle2D at random.
+        /** This function returns one of the vertices {a, b, c} uniformly at random.
+            @param rng A pre-seeded random number generator object that is to be used by this function to generate random values.
+            @see class RNG(), RandomPointInside(), RandomPointOnEdge(), Vertex().*/
+        Vector2 RandomVertex(RNG& rng) const;
+        /// Generates a random point on the edge of this Triangle2D.
+        /** The points are distributed uniformly.
+            This function requires that this triangle is not degenerate. This it is, an assume() error will be printed,
+            and the return value will be undefined.
+            @param rng A pre-seeded random number generator object that is to be used by this function to generate random values.
+            @see class RNG, RandomPointInside, RandomVertex(), Edge(), class LineSegment2D, IsDegenerate(). */
+        Vector2 RandomPointOnEdge(RNG& rng) const;
+        std::string ToString() const;
+        bool Equals(const Triangle2D& rhs, float epsilon = 1e-3f) const {
+            return A.Equals(rhs.A, epsilon) && B.Equals(rhs.B, epsilon) && C.Equals(rhs.C, epsilon);
+        }
     };
+
+    std::ostream& operator << (std::ostream& o, const Triangle2D& triangle);
 }

include/J3ML/Geometry/TriangleMesh.hpp

@@ -1,32 +1,31 @@
 #pragma once
 
 #include <vector>
-#include <J3ML/LinearAlgebra/Common.hpp>
-#include <J3ML/Geometry/Common.hpp>
-#include "J3ML/J3ML.hpp"
-
-
-using namespace J3ML::LinearAlgebra;
+#include <J3ML/LinearAlgebra.hpp>
+#include <J3ML/Geometry/Shape.hpp>
+#include <J3ML/Geometry/Forward.hpp>
 
 namespace J3ML::Geometry
 {
 
-    class TriangleMesh
+    class TriangleMesh : public Shape
     {
     public:
         /// Default constructor for a triangle mesh.
         TriangleMesh(int expectedPolygonCount = 1000);
 
     public:
-        //std::vector<Vector3> Vertices;
-        //std::vector<Vector3> Normals;
-        //std::vector<Vector3> UVs;
-        //std::vector<u64> Indices;
+        std::vector<Vector3> Vertices;
+        std::vector<Vector3> Normals;
+        std::vector<Vector3> UVs;
+        std::vector<u64> Indices;
 
         std::vector<float> GenerateVertexList();
         //std::vector<Triangle> GenerateTriangleList();
+    public:
+
     private:
-        std::vector<float> cachedVertexList;
+
         //std::vector<Triangle> cachedTriangleList;
     };
 

include/J3ML/J3ML.hpp

@@ -84,6 +84,9 @@ namespace J3ML::Math::Constants {
     constexpr float RecipSqrt2Pi = 0.3989422804014326779399460599343818684758586311649346576659258296706579258993018385012523339073069364;
     /// pi - https://www.mathsisfun.com/numbers/pi.html
     constexpr float Pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679;
+    /// pi * 0.5.
+    constexpr float HalfPi = 1.5707963267948966192313216916397514420985846996875529104874722961539082031431044993140174126710585;
+
     /// e - https://www.mathsisfun.com/numbers/e-eulers-number.html
     constexpr float EulersNumber = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274;
     /// 2pi - The ratio of a circle's circumferecne to its radius, and the number of radians in one turn.

include/J3ML/LinearAlgebra/Angle2D.hpp

@@ -1,6 +1,5 @@
 #pragma once
-#include "J3ML/J3ML.hpp"
-
+#include <J3ML/J3ML.hpp>
 
 namespace J3ML::LinearAlgebra {
 

include/J3ML/LinearAlgebra/AxisAngle.hpp

@@ -2,28 +2,23 @@
 
 #include <J3ML/LinearAlgebra/EulerAngle.hpp>
 #include <J3ML/LinearAlgebra/Quaternion.hpp>
-#include <J3ML/LinearAlgebra/AxisAngle.hpp>
 #include <J3ML/LinearAlgebra/Vector3.hpp>
 
-namespace J3ML::LinearAlgebra
-{
+namespace J3ML::LinearAlgebra {
+    class AxisAngle;
+}
 
-    /// Transitional datatype, not useful for internal representation of rotation
-    /// But has uses for conversion and manipulation.
-    class AxisAngle {
-    public:
-        Vector3 axis;
-        float angle;
-    public:
-        AxisAngle();
-        explicit AxisAngle(const Quaternion& q);
-        explicit AxisAngle(const EulerAngle& e);
+/// Transitional datatype, not useful for internal representation of rotation
+/// But has uses for conversion and manipulation.
+class J3ML::LinearAlgebra::AxisAngle {
+public:
+    Vector3 axis;
+    // Radians.
+    float angle;
+public:
+    AxisAngle();
+    explicit AxisAngle(const Quaternion& q);
+    explicit AxisAngle(const EulerAngleXYZ& e);
+    AxisAngle(const Vector3& axis, float angle);
 
-        AxisAngle(const Vector3 &axis, float angle);
-
-        EulerAngle ToEulerAngleXYZ() const;
-
-        Quaternion ToQuaternion() const;
-        static AxisAngle FromEulerAngleXYZ(const EulerAngle&);
-    };
-}
+};

include/J3ML/LinearAlgebra/DirectionVector.hpp

@@ -0,0 +1,20 @@
+#pragma once
+#include <J3ML/LinearAlgebra/Vector3.hpp>
+
+namespace J3ML::LinearAlgebra {
+    class DirectionVectorRH;
+}
+
+/// Direction vector of a given Matrix3x3 RotationMatrix in a Right-handed coordinate space.
+class J3ML::LinearAlgebra::DirectionVectorRH : public Vector3 {
+private:
+    // This is purposefully not exposed because these types aren't usually convertable.
+    explicit DirectionVectorRH(const Vector3& rhs);
+public:
+    static DirectionVectorRH Forward(const Matrix3x3& rhs);
+    static DirectionVectorRH Backward(const Matrix3x3& rhs);
+    static DirectionVectorRH Left(const Matrix3x3& rhs);
+    static DirectionVectorRH Right(const Matrix3x3& rhs);
+    static DirectionVectorRH Up(const Matrix3x3& rhs);
+    static DirectionVectorRH Down(const Matrix3x3& rhs);
+};

include/J3ML/LinearAlgebra/EulerAngle.hpp

@@ -5,48 +5,19 @@
 #include <J3ML/LinearAlgebra/AxisAngle.hpp>
 
 namespace J3ML::LinearAlgebra {
+    class EulerAngleXYZ;
+}
 
-    class AxisAngle;
-
-// Essential Reading:
-// http://www.essentialmath.com/GDC2012/GDC2012_JMV_Rotations.pdf
-class EulerAngle {
+class J3ML::LinearAlgebra::EulerAngleXYZ {
 public:
-    EulerAngle();
-    EulerAngle(float pitch, float yaw, float roll);
-    EulerAngle(const Vector3& vec) : pitch(vec.x), yaw(vec.y), roll(vec.z) {}
-
-    AxisAngle ToAxisAngle() const;
-
-    [[nodiscard]] Quaternion ToQuaternion() const;
-
-
-    explicit EulerAngle(const Quaternion& rhs);
-    explicit EulerAngle(const AxisAngle& rhs);
-
-    /// TODO: Implement separate upper and lower bounds
-    /// Preserves internal value of euler angles, normalizes and clamps the output.
-    /// This does not solve gimbal lock!!!
-    float GetPitch(float pitch_limit) const;
-    float GetYaw(float yaw_limit)     const;
-    float GetRoll(float roll_limit)   const;
-
-    bool operator==(const EulerAngle& a) const;
-    void clamp();
-
-    // TODO: Euler Angles do not represent a vector, length doesn't apply, nor is this information meaningful for this data type.
-    // If you need a meaningful representation of length in 3d space, use a vector!!
-    [[nodiscard]] float length() const {
-        return 0;
-    }
-    // TODO: Implement
-    Vector3 unitVector() const;
-
-    EulerAngle movementAngle() const;
 public:
-    float pitch;
-    float yaw;
-    float roll;
+    float roll = 0; // X
+    float pitch = 0; // Y
+    float yaw = 0; // Z
+public:
+    EulerAngleXYZ(float roll, float pitch, float yaw);
+public:
+    explicit EulerAngleXYZ(const Quaternion& rhs);
+    explicit EulerAngleXYZ(const AxisAngle& rhs);
+    explicit EulerAngleXYZ(const Matrix3x3& rhs);
 };
-
-}

include/J3ML/LinearAlgebra/Forward.hpp

@@ -7,7 +7,7 @@ namespace J3ML::LinearAlgebra
     class Vector3; // A type representing a position in a 3-dimensional coordinate space.
     class Vector4; // A type representing a position in a 4-dimensional coordinate space.
     class Angle2D; // Uses x,y components to represent a 2D rotation.
-    class EulerAngle; // Uses pitch,yaw,roll components to represent a 3D orientation.
+    class EulerAngleXYZ; // Uses pitch,yaw,roll components to represent a 3D orientation.
     class AxisAngle; //
     class CoordinateFrame; //
     class Matrix2x2;
@@ -15,6 +15,7 @@ namespace J3ML::LinearAlgebra
     class Matrix4x4;
     class Transform2D;
     class Transform3D;
+    class DirectionVectorRH; // A type representing a direction in 3D space.
     class Quaternion;
 
 

include/J3ML/LinearAlgebra/Matrix2x2.hpp

@@ -14,14 +14,13 @@
 
 
 
-#include <J3ML/LinearAlgebra/Common.hpp>
-
+#include <J3ML/LinearAlgebra.hpp>
 
 namespace J3ML::LinearAlgebra {
     class Matrix2x2 {
     public:
-        enum { Rows = 3 };
-        enum { Cols = 3 };
+        enum { Rows = 2 };
+        enum { Cols = 2 };
         static const Matrix2x2 Zero;
         static const Matrix2x2 Identity;
         static const Matrix2x2 NaN;
@@ -52,6 +51,9 @@ namespace J3ML::LinearAlgebra {
         Vector2 operator * (const Vector2& rhs) const;
         Matrix2x2 operator * (const Matrix2x2 &rhs) const;
 
+        inline float* ptr() { return &elems[0][0];}
+        [[nodiscard]] inline const float* ptr() const {return &elems[0][0];}
+
     protected:
         float elems[2][2];
     };

include/J3ML/LinearAlgebra/Matrix3x3.hpp

@@ -63,7 +63,10 @@ namespace J3ML::LinearAlgebra {
         /// Constructs this matrix3x3 from the given quaternion.
         explicit Matrix3x3(const Quaternion& orientation);
         /// Constructs this matrix3x3 from the given euler angle.
-        explicit Matrix3x3(const EulerAngle& orientation);
+
+        explicit Matrix3x3(const EulerAngleXYZ& orientation);
+
+        explicit Matrix3x3(const AxisAngle& orientation);
 
         /// Constructs this Matrix3x3 from a pointer to an array of floats.
         explicit Matrix3x3(const float *data);
@@ -153,6 +156,7 @@ namespace J3ML::LinearAlgebra {
 
         /// Sets this matrix to perform a rotation about the given axis and angle.
         void SetRotatePart(const Vector3& a, float angle);
+
         void SetRotatePart(const AxisAngle& axisAngle);
         /// Sets this matrix to perform the rotation expressed by the given quaternion.
         void SetRotatePart(const Quaternion& quat);
@@ -236,20 +240,13 @@ namespace J3ML::LinearAlgebra {
         /// Accesses this structure as a float array.
         /// @return A pointer to the upper-left element. The data is contiguous in memory.
         /// ptr[0] gives the element [0][0], ptr[1] is [0][1]. ptr[2] is [0][2]
-        inline float *ptr() { return &elems[0][0];}
-        [[nodiscard]] inline const float *ptr() const {return &elems[0][0];}
+        inline float* ptr() { return &elems[0][0];}
+        [[nodiscard]] inline const float* ptr() const {return &elems[0][0];}
 
-
-        /// Convers this rotation matrix to a quaternion.
-        /// This function assumes that the matrix is orthonormal (no shear or scaling) and does not perform any mirroring (determinant > 0)
-        [[nodiscard]] Quaternion ToQuat() const;
         /// Attempts to convert this matrix to a quaternion. Returns false if the conversion cannot succeed (this matrix was not a rotation
         /// matrix, and there is scaling ,shearing, or mirroring in this matrix)
         bool TryConvertToQuat(Quaternion& q) const;
 
-        /// Converts this rotation matrix to an Euler Angle.
-        [[nodiscard]] EulerAngle ToEulerAngle() const;
-
         /// Returns the main diagonal.
         /// The main diagonal consists of the elements at m[0][0], m[1][1], m[2][2]
         [[nodiscard]] Vector3 Diagonal() const;

include/J3ML/LinearAlgebra/Matrix4x4.hpp

@@ -1,7 +1,7 @@
 #pragma once
 
 
-#include <J3ML/LinearAlgebra/Common.hpp>
+#include <J3ML/LinearAlgebra/Forward.hpp>
 #include <J3ML/Algorithm/RNG.hpp>
 
 #include <algorithm>
@@ -71,8 +71,7 @@ namespace J3ML::LinearAlgebra {
 
         /// Constructs this Matrix4x4 from the given quaternion.
         explicit Matrix4x4(const Quaternion& orientation);
-        /// Constructs this Matrix4x4 from the given Euler Angle.
-        explicit Matrix4x4(const EulerAngle& orientation);
+
 
         /// Constructs this float4x4 from the given quaternion and translation.
         /// Logically, the translation occurs after the rotation has been performed.
@@ -567,8 +566,6 @@ namespace J3ML::LinearAlgebra {
 
         [[nodiscard]] Quaternion ToQuat() const;
 
-        [[nodiscard]] EulerAngle ToEulerAngle() const;
-
         /// Returns true if this Matrix4x4 is equal to the given Matrix4x4, up to given per-element epsilon.
         bool Equals(const Matrix4x4& other, float epsilon = 1e-3f) const;
 

include/J3ML/LinearAlgebra/Quaternion.hpp

@@ -1,261 +1,240 @@
 #pragma once
 
-#include <J3ML/LinearAlgebra/Common.hpp>
+#include <J3ML/LinearAlgebra/Forward.hpp>
 #include <J3ML/Algorithm/RNG.hpp>
 #include <cmath>
 
-namespace J3ML::LinearAlgebra
-{
-    class Quaternion {
-    public:
-        /// The identity quaternion performs no rotation when applied to a vector.
-        static const Quaternion Identity;
-        /// A compile-time constant Quaternion with the value (NAN, NAN, NAN, NAN).
-        /// For this constant, each element has the value of quiet NAN, or Not-A-Number.
-        /// @note Never compare a Quaternion to this value! Due to how IEEE floats work, "nan == nan" returns false!
-        /// That is, nothing is equal to NaN, not even NaN itself!
-        static const Quaternion NaN;
-    public:
-        /// The default constructor does not initialize any member values.
-        Quaternion();
-        /// Copy constructor
-        Quaternion(const Quaternion &rhs) = default;
-        /// Constructs a quaternion from the given data buffer.
-        /// @param data An array of four floats to use for the quaternion, in the order 'x,y,z,w.' (== 'i,j,k,w')
-        explicit Quaternion(const float *data);
+namespace J3ML::LinearAlgebra {
+    class Quaternion;
+}
 
-        explicit Quaternion(const Matrix3x3 &rotationMtrx);
-        explicit Quaternion(const Matrix4x4 &rotationMtrx);
+class J3ML::LinearAlgebra::Quaternion {
+public:
+    float x;
+    float y;
+    float z;
+    float w;
+public:
+    /// The identity quaternion performs no rotation when applied to a vector.
+    static const Quaternion Identity;
+    /// A compile-time constant Quaternion with the value (NAN, NAN, NAN, NAN).
+    /// For this constant, each element has the value of quiet NAN, or Not-A-Number.
+    /// @note Never compare a Quaternion to this value! Due to how IEEE floats work, "nan == nan" returns false!
+    /// That is, nothing is equal to NaN, not even NaN itself!
+    static const Quaternion NaN;
+public:
+    /// The default constructor does not initialize any member values.
+    Quaternion() = default;
+    /// Copy constructor
+    Quaternion(const Quaternion &rhs);
+    /// Quaternion from Matrix3x3
+    explicit Quaternion(const Matrix3x3& ro_mat);
+    /// Quaternion from Matrix4x4 RotatePart.
+    explicit Quaternion(const Matrix4x4& ro_mat);
+    /// Quaternion from EulerAngleXYZ.
+    explicit Quaternion(const EulerAngleXYZ& rhs);
+    /// Quaternion from AxisAngle.
+    explicit Quaternion(const AxisAngle& angle);
+    /// Quaternion from Vector4 (no conversion).
+    explicit Quaternion(const Vector4& vector4);
 
-        /// @param x The factor of i.
-        /// @param y The factor of j.
-        /// @param z The factor of k.
-        /// @param w The scalar factor (or 'w').
-        /// @note The input data is not normalized after construction, this has to be done manually.
-        Quaternion(float X, float Y, float Z, float W);
+    /// @param x The factor of i.
+    /// @param y The factor of j.
+    /// @param z The factor of k.
+    /// @param w The scalar factor (or 'w').
+    /// @note The input data is not normalized after construction, this has to be done manually.
+    Quaternion(float X, float Y, float Z, float W);
 
-        /// Constructs this quaternion by specifying a rotation axis and the amount of rotation to be performed about that axis
-        /// @param rotationAxis The normalized rotation axis to rotate about. If using Vector4 version of the constructor, the w component of this vector must be 0.
-        /// @param rotationAngleRadians The angle to rotate by, in radians. For example, Pi/4.f equals to 45 degrees, Pi/2.f is 90 degrees, etc.
-        /// @see DegToRad()
-        Quaternion(const Vector3 &rotationAxis, float rotationAngleRadians);
-        Quaternion(const Vector4 &rotationAxis, float rotationAngleRadians);
+    /// Creates a LookAt quaternion.
+    /** A LookAt quaternion is a quaternion that orients an object to face towards a specified target direction.
+    @param localForward Specifies the forward direction in the local space of the object. This is the direction
+    the model is facing at in its own local/object space, often +X(1,0,0), +Y(0,1,0), or +Z(0,0,1).
+    The vector to pass in here depends on the conventions you or your modeling software is using, and it is best
+    to pick one convention for all your objects, and be consistent.
+    This input parameter must be a normalized vector.
+    @param targetDirection Specifies the desired world space direction the object should look at. This function
+    will compute a quaternion which will rotate the localForward vector to orient towards this targetDirection
+    vector. This input parameter must be a normalized vector.
+    @param localUp Specifies the up direction in the local space of the object. This is the up direction the model
+    was authored in, often +Y (0,1,0) or +Z (0,0,1). The vector to pass in here depends on the conventions you
+    or your modeling software is using, and it is best to pick one convention for all your objects, and be
+    consistent. This input parameter must be a normalized vector. This vector must be perpendicular to the
+    vector localForward, i.e. localForward.Dot(localUp) == 0.
+    @param worldUp Specifies the global up direction of the scene in world space. Simply rotating one vector to
+    coincide with another (localForward->targetDirection) would cause the up direction of the resulting
+    orientation to drift (e.g. the model could be looking at its target its head slanted sideways). To keep
+    the up direction straight, this function orients the localUp direction of the model to point towards
+    the specified worldUp direction (as closely as possible). The worldUp and targetDirection vectors cannot be
+    collinear, but they do not need to be perpendicular either.
+    @return A quaternion that maps the given local space forward direction vector to point towards the given target
+    direction, and the given local up direction towards the given target world up direction. For the returned
+    quaternion Q it holds that M * localForward = targetDirection, and M * localUp lies in the plane spanned
+    by the vectors targetDirection and worldUp.
+    @see RotateFromTo() */
+    static Quaternion LookAt(const Vector3& localForward, const Vector3& targetDirection, const Vector3& localUp, const Vector3& worldUp);
 
-        explicit Quaternion(const Vector4& vector4);
-        explicit Quaternion(const EulerAngle& angle);
-        explicit Quaternion(const AxisAngle& angle);
+     /// Creates a new quaternion that rotates about the positive X axis by the given rotation.
+     static Quaternion RotateX(float rad);
+     /// Creates a new quaternion that rotates about the positive Y axis by the given rotation.
+     static Quaternion RotateY(float rad);
+     /// Creates a new quaternion that rotates about the positive Z axis by the given rotation.
+     static Quaternion RotateZ(float rad);
 
-        /// Creates a LookAt quaternion.
-        /** A LookAt quaternion is a quaternion that orients an object to face towards a specified target direction.
-            @param localForward Specifies the forward direction in the local space of the object. This is the direction
-                the model is facing at in its own local/object space, often +X(1,0,0), +Y(0,1,0), or +Z(0,0,1).
-                The vector to pass in here depends on the conventions you or your modeling software is using, and it is best
-                to pick one convention for all your objects, and be consistent.
-                This input parameter must be a normalized vector.
-            @param targetDirection Specifies the desired world space direction the object should look at. This function
-                will compute a quaternion which will rotate the localForward vector to orient towards this targetDirection
-                vector. This input parameter must be a normalized vector.
-            @param localUp Specifies the up direction in the local space of the object. This is the up direction the model
-                was authored in, often +Y (0,1,0) or +Z (0,0,1). The vector to pass in here depends on the conventions you
-                or your modeling software is using, and it is best to pick one convention for all your objects, and be
-                consistent. This input parameter must be a normalized vector. This vector must be perpendicular to the
-                vector localForward, i.e. localForward.Dot(localUp) == 0.
-            @param worldUp Specifies the global up direction of the scene in world space. Simply rotating one vector to
-                coincide with another (localForward->targetDirection) would cause the up direction of the resulting
-                orientation to drift (e.g. the model could be looking at its target its head slanted sideways). To keep
-                the up direction straight, this function orients the localUp direction of the model to point towards
-                the specified worldUp direction (as closely as possible). The worldUp and targetDirection vectors cannot be
-                collinear, but they do not need to be perpendicular either.
-            @return A quaternion that maps the given local space forward direction vector to point towards the given target
-                direction, and the given local up direction towards the given target world up direction. For the returned
-                quaternion Q it holds that M * localForward = targetDirection, and M * localUp lies in the plane spanned
-                by the vectors targetDirection and worldUp.
-            @see RotateFromTo() */
-        static Quaternion LookAt(const Vector3& localForward, const Vector3& targetDirection, const Vector3& localUp, const Vector3& worldUp);
+     /// Creates a new quaternion that rotates sourceDirection vector (in world space) to coincide with the
+     /// targetDirection vector (in world space).
+     /// Rotation is performed about the origin.
+     /// The vectors sourceDirection and targetDirection are assumed to be normalized.
+     /// @note There are multiple such rotations - this function returns the rotation that has the shortest angle
+     /// (when decomposed to axis-angle notation).
+     static Quaternion RotateFromTo(const Vector3& sourceDirection, const Vector3& targetDirection);
+     static Quaternion RotateFromTo(const Vector4& sourceDirection, const Vector4& targetDirection);
 
-        /// Creates a new Quaternion that rotates about the given axis by the given angle.
-        static Quaternion RotateAxisAngle(const AxisAngle& axisAngle);
-
-        /// Creates a new quaternion that rotates about the positive X axis by the given rotation.
-        static Quaternion RotateX(float angleRadians);
-        /// Creates a new quaternion that rotates about the positive Y axis by the given rotation.
-        static Quaternion RotateY(float angleRadians);
-        /// Creates a new quaternion that rotates about the positive Z axis by the given rotation.
-        static Quaternion RotateZ(float angleRadians);
-
-        /// Creates a new quaternion that rotates sourceDirection vector (in world space) to coincide with the
-        /// targetDirection vector (in world space).
-        /// Rotation is performed about the origin.
-        /// The vectors sourceDirection and targetDirection are assumed to be normalized.
-        /// @note There are multiple such rotations - this function returns the rotation that has the shortest angle
-        /// (when decomposed to axis-angle notation).
-        static Quaternion RotateFromTo(const Vector3& sourceDirection, const Vector3& targetDirection);
-        static Quaternion RotateFromTo(const Vector4& sourceDirection, const Vector4& targetDirection);
-
-        /// Creates a new quaternion that
-        /// 1. rotates sourceDirection vector to coincide with the targetDirection vector, and then
-        /// 2. rotates sourceDirection2 (which was transformed by 1.) to targetDirection2, but keeping the constraint that
-        /// sourceDirection must look at targetDirection
-        static Quaternion RotateFromTo(const Vector3& sourceDirection, const Vector3& targetDirection, const Vector3& sourceDirection2, const Vector3& targetDirection2);
+     /// Creates a new quaternion that
+     /// 1. rotates sourceDirection vector to coincide with the targetDirection vector, and then
+     /// 2. rotates sourceDirection2 (which was transformed by 1.) to targetDirection2, but keeping the constraint that
+     /// sourceDirection must look at targetDirection
+     static Quaternion RotateFromTo(const Vector3& sourceDirection, const Vector3& targetDirection, const Vector3& sourceDirection2, const Vector3& targetDirection2);
 
 
-        /// Returns a uniformly random unitary quaternion.
-        static Quaternion RandomRotation(RNG &rng);
-    public:
-        void SetFromAxisAngle(const Vector3 &vector3, float between);
+     /// Returns a uniformly random unitary quaternion.
+     static Quaternion RandomRotation(RNG &rng);
+public:
+    /// Inverses this quaternion in-place.
+    /// @note For optimization purposes, this function assumes that the quaternion is unitary, in which
+    /// case the inverse of the quaternion is simply just the same as its conjugate.
+    /// This function does not detect whether the operation succeeded or failed.
+    void Inverse();
 
-        void SetFromAxisAngle(const Vector4 &vector4, float between);
-        void SetFrom(const AxisAngle& angle);
+    /// Returns an inverted copy of this quaternion.
+    [[nodiscard]] Quaternion Inverted() const;
+    /// Computes the conjugate of this quaternion in-place.
+    void Conjugate();
+    /// Returns a conjugated copy of this quaternion.
+    [[nodiscard]] Quaternion Conjugated() const;
 
-        /// Inverses this quaternion in-place.
-        /// @note For optimization purposes, this function assumes that the quaternion is unitary, in which
-        /// case the inverse of the quaternion is simply just the same as its conjugate.
-        /// This function does not detect whether the operation succeeded or failed.
-        void Inverse();
+    /// Inverses this quaternion in-place.
+    /// Call this function when the quaternion is not known beforehand to be normalized.
+    /// This function computes the inverse proper, and normalizes the result.
+    /// @note Because of the normalization, it does not necessarily hold that q * q.InverseAndNormalize() == id.
+    /// @return Returns the old length of this quaternion (not the old length of the inverse quaternion).
+    float InverseAndNormalize();
 
-        /// Returns an inverted copy of this quaternion.
-        [[nodiscard]] Quaternion Inverted() const;
-        /// Computes the conjugate of this quaternion in-place.
-        void Conjugate();
-        /// Returns a conjugated copy of this quaternion.
-        [[nodiscard]] Quaternion Conjugated() const;
+    /// Returns the local +X axis in the post-transformed coordinate space. This is the same as transforming the vector (1,0,0) by this quaternion.
+    [[nodiscard]] Vector3 WorldX() const;
+    /// Returns the local +Y axis in the post-transformed coordinate space. This is the same as transforming the vector (0,1,0) by this quaternion.
+    [[nodiscard]] Vector3 WorldY() const;
+    /// Returns the local +Z axis in the post-transformed coordinate space. This is the same as transforming the vector (0,0,1) by this quaternion.
+    [[nodiscard]] Vector3 WorldZ() const;
+    /// Returns the axis of rotation for this quaternion.
+    [[nodiscard]] Vector3 Axis() const;
 
-        /// Inverses this quaternion in-place.
-        /// Call this function when the quaternion is not known beforehand to be normalized.
-        /// This function computes the inverse proper, and normalizes the result.
-        /// @note Because of the normalization, it does not necessarily hold that q * q.InverseAndNormalize() == id.
-        /// @return Returns the old length of this quaternion (not the old length of the inverse quaternion).
-        float InverseAndNormalize();
+    /// Returns the angle of rotation for this quaternion, in radians.
+    [[nodiscard]] float Angle() const;
 
-        /// Returns the local +X axis in the post-transformed coordinate space. This is the same as transforming the vector (1,0,0) by this quaternion.
-        [[nodiscard]] Vector3 WorldX() const;
-        /// Returns the local +Y axis in the post-transformed coordinate space. This is the same as transforming the vector (0,1,0) by this quaternion.
-        [[nodiscard]] Vector3 WorldY() const;
-        /// Returns the local +Z axis in the post-transformed coordinate space. This is the same as transforming the vector (0,0,1) by this quaternion.
-        [[nodiscard]] Vector3 WorldZ() const;
+    [[nodiscard]] float LengthSquared() const;
+    [[nodiscard]] float Length() const;
 
-        /// Returns the axis of rotation for this quaternion.
-        [[nodiscard]] Vector3 Axis() const;
+    [[nodiscard]] Matrix3x3 ToMatrix3x3() const;
+    [[nodiscard]] Matrix4x4 ToMatrix4x4() const;
 
-        /// Returns the angle of rotation for this quaternion, in radians.
-        [[nodiscard]] float Angle() const;
+    [[nodiscard]] Matrix4x4 ToMatrix4x4(const Vector3 &translation) const;
 
-        [[nodiscard]] float LengthSquared() const;
-        [[nodiscard]] float Length() const;
+    [[nodiscard]] Vector3 Transform(const Vector3& vec) const;
+    [[nodiscard]] Vector3 Transform(float X, float Y, float Z) const;
+    // Note: We only transform the x,y,z components of 4D vectors, w is left untouched
+    [[nodiscard]] Vector4 Transform(const Vector4& vec) const;
+    [[nodiscard]] Vector4 Transform(float X, float Y, float Z, float W) const;
 
-        [[nodiscard]] EulerAngle ToEulerAngle() const;
+    [[nodiscard]] Quaternion Lerp(const Quaternion& b, float t) const;
+    static Quaternion Lerp(const Quaternion &source, const Quaternion& target, float t);
+    [[nodiscard]] Quaternion Slerp(const Quaternion& q2, float t) const;
+    static Quaternion Slerp(const Quaternion &source, const Quaternion& target, float t);
+
+    /// Returns the 'from' vector rotated towards the 'to' vector by the given normalized time parameter.
+    /** This function slerps the given 'form' vector toward the 'to' vector.
+    @param from A normalized direction vector specifying the direction of rotation at t=0
+    @param to A normalized direction vector specifying the direction of rotation at t=1
+    @param t The interpolation time parameter, in the range [0, 1]. Input values outside this range are
+    silently clamped to the [0, 1] interval.
+    @return A spherical linear interpolation of the vector 'from' towards the vector 'to'. */
+    static Vector3 SlerpVector(const Vector3& from, const Vector3& to, float t);
+
+    /// Returns the 'from' vector rotated towards the 'to' vector by the given absolute angle, in radians.
+    /** This function slerps the given 'from' vector towards the 'to' vector.
+    @param from A normalized direction vector specifying the direction of rotation at angleRadians=0.
+    @param to A normalized direction vector specifying the target direction to rotate towards.
+    @param angleRadians The maximum angle to rotate the 'from' vector by, in the range [0, pi]. If the
+    angle between 'from' and 'to' is smaller than this angle, then the vector 'to' is returned.
+    Input values outside this range are silently clamped to the [0, pi] interval.
+    @return A spherical linear interpolation of the vector 'from' towards the vector 'to'. */
+    static Vector3 SlerpVectorAbs(const Vector3 &from, const Vector3& to, float angleRadians);
+
+    /// Normalizes this quaternion in-place.
+    /// @returns false if failure, true if success.
+    [[nodiscard]] bool Normalize();
+    /// Returns a normalized copy of this quaternion.
+    [[nodiscard]] Quaternion Normalized() const;
+
+    /// Returns true if the length of this quaternion is one.
+    [[nodiscard]] bool IsNormalized(float epsilon = 1e-5f) const;
+    [[nodiscard]] bool IsInvertible(float epsilon = 1e-3f) const;
+
+    /// Returns true if the entries of this quaternion are all finite.
+    [[nodiscard]] bool IsFinite() const;
+
+    /// Returns true if this quaternion equals rhs, up to the given epsilon.
+    [[nodiscard]] bool Equals(const Quaternion& rhs, float epsilon = 1e-3f) const;
+
+    /// Compares whether this Quaternion and the given Quaternion are identical bit-by-bit in the underlying representation.
+    /// @note Prefer using this over e.g. memcmp, since there can be SSE-related padding in the structures.
+    bool BitEquals(const Quaternion& rhs) const;
+
+    /// @return A pointer to the first element (x). The data is contiguous in memory.
+    /// ptr[0] gives x, ptr[1] gives y, ptr[2] gives z, ptr[3] gives w.
+    inline float *ptr() { return &x; }
+    [[nodiscard]] inline const float *ptr() const { return &x; }
 
 
-        [[nodiscard]] Matrix3x3 ToMatrix3x3() const;
-        [[nodiscard]] Matrix4x4 ToMatrix4x4() const;
+    // Multiplies two quaternions together.
+    // The product q1 * q2 returns a quaternion that concatenates the two orientation rotations.
+    // The rotation q2 is applied first before q1.
+    Quaternion operator * (const Quaternion& rhs) const;
 
-        [[nodiscard]] Matrix4x4 ToMatrix4x4(const Vector3 &translation) const;
+    // Unsafe
+    Quaternion operator * (float scalar) const;
 
-        [[nodiscard]] Vector3 Transform(const Vector3& vec) const;
-        [[nodiscard]] Vector3 Transform(float X, float Y, float Z) const;
-        // Note: We only transform the x,y,z components of 4D vectors, w is left untouched
-        [[nodiscard]] Vector4 Transform(const Vector4& vec) const;
-        [[nodiscard]] Vector4 Transform(float X, float Y, float Z, float W) const;
+    // Unsafe
+    Quaternion operator / (float scalar) const;
 
-        [[nodiscard]] Quaternion Lerp(const Quaternion& b, float t) const;
-        static Quaternion Lerp(const Quaternion &source, const Quaternion& target, float t);
-        [[nodiscard]] Quaternion Slerp(const Quaternion& q2, float t) const;
-        static Quaternion Slerp(const Quaternion &source, const Quaternion& target, float t);
+    // Transforms the given vector by this Quaternion.
+    Vector3 operator * (const Vector3& rhs) const;
 
-        /// Returns the 'from' vector rotated towards the 'to' vector by the given normalized time parameter.
-        /** This function slerps the given 'form' vector toward the 'to' vector.
-            @param from A normalized direction vector specifying the direction of rotation at t=0
-            @param to A normalized direction vector specifying the direction of rotation at t=1
-            @param t The interpolation time parameter, in the range [0, 1]. Input values outside this range are
-            silently clamped to the [0, 1] interval.
-            @return A spherical linear interpolation of the vector 'from' towards the vector 'to'. */
-        static Vector3 SlerpVector(const Vector3& from, const Vector3& to, float t);
+    Vector4 operator * (const Vector4& rhs) const;
 
-        /// Returns the 'from' vector rotated towards the 'to' vector by the given absolute angle, in radians.
-        /** This function slerps the given 'from' vector towards the 'to' vector.
-            @param from A normalized direction vector specifying the direction of rotation at angleRadians=0.
-            @param to A normalized direction vector specifying the target direction to rotate towards.
-            @param angleRadians The maximum angle to rotate the 'from' vector by, in the range [0, pi]. If the
-                angle between 'from' and 'to' is smaller than this angle, then the vector 'to' is returned.
-                Input values outside this range are silently clamped to the [0, pi] interval.
-            @return A spherical linear interpolation of the vector 'from' towards the vector 'to'. */
-        static Vector3 SlerpVectorAbs(const Vector3 &from, const Vector3& to, float angleRadians);
-
-        /// Normalizes this quaternion in-place.
-        /// Returns the old length of this quaternion, or 0 if normalization failed.
-        float Normalize();
-        /// Returns a normalized copy of this quaternion.
-        [[nodiscard]] Quaternion Normalized() const;
-
-        /// Returns true if the length of this quaternion is one.
-        [[nodiscard]] bool IsNormalized(float epsilon = 1e-5f) const;
-        [[nodiscard]] bool IsInvertible(float epsilon = 1e-3f) const;
-
-        /// Returns true if the entries of this quaternion are all finite.
-        [[nodiscard]] bool IsFinite() const;
-
-        /// Returns true if this quaternion equals rhs, up to the given epsilon.
-        [[nodiscard]] bool Equals(const Quaternion& rhs, float epsilon = 1e-3f) const;
-
-        /// Compares whether this Quaternion and the given Quaternion are identical bit-by-bit in the underlying representation.
-        /// @note Prefer using this over e.g. memcmp, since there can be SSE-related padding in the structures.
-        bool BitEquals(const Quaternion& rhs) const;
-
-        /// @return A pointer to the first element (x). The data is contiguous in memory.
-        /// ptr[0] gives x, ptr[1] gives y, ptr[2] gives z, ptr[3] gives w.
-        inline float *ptr() { return &x; }
-        [[nodiscard]] inline const float *ptr() const { return &x; }
-
-
-        // Multiplies two quaternions together.
-        // The product q1 * q2 returns a quaternion that concatenates the two orientation rotations.
-        // The rotation q2 is applied first before q1.
-        Quaternion operator * (const Quaternion& rhs) const;
-
-        // Unsafe
-        Quaternion operator * (float scalar) const;
-
-        // Unsafe
-        Quaternion operator / (float scalar) const;
-
-        // Transforms the given vector by this Quaternion.
-        Vector3 operator * (const Vector3& rhs) const;
-
-        Vector4 operator * (const Vector4& rhs) const;
-
-        // Divides a quaternion by another. Divison "a / b" results in a quaternion that rotates the orientation b to coincide with orientation of
-        Quaternion operator / (const Quaternion& rhs) const;
-        Quaternion operator + (const Quaternion& rhs) const;
+    // Divides a quaternion by another. Divison "a / b" results in a quaternion that rotates the orientation b to coincide with orientation of
+    Quaternion operator / (const Quaternion& rhs) const;
+    Quaternion operator + (const Quaternion& rhs) const;
 
 
 
-        Quaternion operator + () const;
-        Quaternion operator - () const;
+    Quaternion operator + () const;
+    Quaternion operator - () const;
 
-        /// Computes the dot product of this and the given quaternion.
-        /// Dot product is commutative.
-        [[nodiscard]] float Dot(const Quaternion &quaternion) const;
+    /// Computes the dot product of this and the given quaternion.
+    /// Dot product is commutative.
+    [[nodiscard]] float Dot(const Quaternion &quaternion) const;
 
-        /// Returns the angle between this and the target orientation (the shortest route) in radians.
-        [[nodiscard]] float AngleBetween(const Quaternion& target) const;
-        /// Returns the axis of rotation to get from this orientation to target orientation (the shortest route).
-        [[nodiscard]] Vector3 AxisFromTo(const Quaternion& target) const;
+    /// Returns the angle between this and the target orientation (the shortest route) in radians.
+    [[nodiscard]] float AngleBetween(const Quaternion& target) const;
+    /// Returns the axis of rotation to get from this orientation to target orientation (the shortest route).
+    [[nodiscard]] Vector3 AxisFromTo(const Quaternion& target) const;
 
 
-        [[nodiscard]] AxisAngle ToAxisAngle() const;
-        void SetFromAxisAngle(const AxisAngle& axisAngle);
+    /// Sets this quaternion to represent the same rotation as the given matrix.
+    void Set(const Matrix3x3& matrix);
+    void Set(const Matrix4x4& matrix);
+    void Set(float x, float y, float z, float w);
+    void Set(const Quaternion& q);
+    void Set(const Vector4& v);
 
-        /// Sets this quaternion to represent the same rotation as the given matrix.
-        void Set(const Matrix3x3& matrix);
-        void Set(const Matrix4x4& matrix);
-        void Set(float x, float y, float z, float w);
-        void Set(const Quaternion& q);
-        void Set(const Vector4& v);
-
-    public:
-        float x;
-        float y;
-        float z;
-        float w;
-    };
-}
+};

include/J3ML/LinearAlgebra/Vector2.hpp

@@ -3,6 +3,7 @@
 
 #include <cstddef>
 #include <J3ML/Algorithm/RNG.hpp>
+#include <J3ML/LinearAlgebra/Forward.hpp>
 
 namespace J3ML::LinearAlgebra {
     using namespace J3ML;
@@ -239,8 +240,8 @@ namespace J3ML::LinearAlgebra {
         /// Returns a normalized copy of this vector.
         /** @note If the vector is zero and cannot be normalized, the vector (1, 0) is returned, and an error message is printed.
             If you do not want to generate an error message on failure, but want to handle the failure yourself, use the
-            Normalize() function instead.
-            @see Normalize() */
+            Normalized() function instead.
+            @see Normalized() */
         [[nodiscard]] Vector2 Normalized() const;
 
         /// Linearly interpolates between two points.
@@ -382,7 +383,17 @@ namespace J3ML::LinearAlgebra {
 
         /// Returns a random Vector2 with each entry randomized between the range [minElem, maxElem].
         static Vector2 RandomBox(Algorithm::RNG& rng, float minElem, float maxElem);
+
+        [[nodiscard]] std::string ToString() const;
     };
 
     Vector2 operator*(float lhs, const Vector2 &rhs);
+
+    std::ostream& operator << (std::ostream& out, const Vector2& rhs);
+
+
+    Vector2 Mul2D(const Matrix3x3& transform, const Vector2& v);
+    Vector2 MulPos2D(const Matrix4x4& transform, const Vector2& v);
+    Vector2 MulDir2D(const Matrix4x4& transform, const Vector2& v);
+
 }

include/J3ML/LinearAlgebra/Vector3.hpp

@@ -68,15 +68,15 @@ public:
     /** @note Due to static data initialization order being undefined in C++, do NOT use this
         member to initialize other static data in other compilation units! */
     static const Vector3 NegativeInfinity;
-    /// Specifies a compile-time constant Vector3 with value (1,1,1).
+    /// Specifies a compile-time constant Vector3 with value (1,0,0).
     /** @note Due to static data initialization order being undefined in C++, do NOT use this
         member to initialize other static data in other compilation units! */
     static const Vector3 UnitX;
-    /// Specifies a compile-time constant Vector3 with value (1,1,1).
+    /// Specifies a compile-time constant Vector3 with value (0,1,0).
     /** @note Due to static data initialization order being undefined in C++, do NOT use this
         member to initialize other static data in other compilation units! */
     static const Vector3 UnitY;
-    /// Specifies a compile-time constant Vector3 with value (1,1,1).
+    /// Specifies a compile-time constant Vector3 with value (0,0,1).
     /** @note Due to static data initialization order being undefined in C++, do NOT use this
         member to initialize other static data in other compilation units! */
     static const Vector3 UnitZ;
@@ -178,6 +178,13 @@ public:
         @see ScaledToLength(). */
     float ScaleToLength(float newLength);
 
+    /// Returns a scaled copy of this vector which has its new length as given.
+    /** This function assumes the length of this vector is not zero. In the case of failure, an error message is printed,
+        and the vector (newLength, 0, 0) is returned.
+        @see ScaleToLength(). */
+    [[nodiscard]] Vector3 ScaledToLength(float newLength) const;
+
+
 
     [[nodiscard]] Vector3 ProjectToNorm(const Vector3& direction) const;
 
@@ -225,11 +232,15 @@ public:
         @see Perpendicular(), Cross(). */
     Vector3 AnotherPerpendicular(const Vector3& hint = Vector3(0,1,0), const Vector3& hint2 = Vector3(0,0,1)) const;
 
-    /// Returns a scaled copy of this vector which has its new length as given.
-    /** This function assumes the length of this vector is not zero. In the case of failure, an error message is printed,
-        and the vector (newLength, 0, 0) is returned.
-        @see ScaleToLength(). */
-    Vector3 ScaledToLength(float newLength) const;
+    /// Completes this vector to generate a perpendicular basis.
+    /** This function computes two new vectors b and c which are both orthogonal to this vector and to each other.
+        That is, the set { this, b, c} is an orthogonal set. The vectors b and c that are outputted are also normalized.
+        @param outB [out] Receives vector b.
+        @param outC [out] Receives vector c.
+        @note When calling this function, this vector should not be zero! */
+    void PerpendicularBasis(Vector3& outB, Vector3& outC) const;
+
+
 
     [[nodiscard]] Vector3 Min(const Vector3& min) const;
     static Vector3 Min(const Vector3& a, const Vector3& b, const Vector3& c);
@@ -282,8 +293,8 @@ public:
 
     /// Returns a copy of this vector, resized to have a magnitude of 1, while preserving "direction"
     /// @note If the vector is zero and cannot be normalized, the vector (1, 0, 0) is returned, and an error message is printed.
-    /// If you do not want to generate an error message on failure, but want to handle the failure yourself, use the Normalize() function instead.
-    /// @see Normalize()
+    /// If you do not want to generate an error message on failure, but want to handle the failure yourself, use the Normalized() function instead.
+    /// @see Normalized()
     [[nodiscard]] Vector3 Normalized() const;
     static Vector3 Normalized(const Vector3& targ);
 

include/J3ML/LinearAlgebra/Vector4.hpp

@@ -3,7 +3,7 @@
 #include <cstddef>
 #include <cmath>
 
-#include <J3ML/LinearAlgebra/Common.hpp>
+#include <J3ML/LinearAlgebra/Forward.hpp>
 
 namespace J3ML::LinearAlgebra {
 
@@ -217,7 +217,7 @@ namespace J3ML::LinearAlgebra {
         [[nodiscard]] Vector4 Cross3(const Vector4& rhs) const;
         [[nodiscard]] Vector4 Cross(const Vector4& rhs) const;
 
-        [[nodiscard]] Vector4 Normalize() const;
+        [[nodiscard]] Vector4 Normalized() const;
         [[nodiscard]] Vector4 Lerp(const Vector4& goal, float alpha) const;
 
         [[nodiscard]] float AngleBetween(const Vector4& rhs) const;
@@ -265,6 +265,8 @@ namespace J3ML::LinearAlgebra {
          float y = 0;
          float z = 0;
          float w = 0;
+
+        void Normalize();
     };
 
 }

src/J3ML/Algorithm/Bezier.cpp

@@ -1,11 +1,13 @@
 #include <J3ML/Algorithm/Bezier.hpp>
+#include <J3ML/LinearAlgebra/Vector2.hpp>
+#include <J3ML/LinearAlgebra/Vector3.hpp>
 
 namespace J3ML::Algorithm
 {
     using namespace J3ML::LinearAlgebra;
     Vector2 BezierNormal(float t, const Vector2 &p0, const Vector2 &p1,
                                   const Vector2 &p2, const Vector2 &p3) {
-        auto derived = BezierDerivative(t, p0, p1, p2, p3);
+        Vector2 derived = BezierDerivative(t, p0, p1, p2, p3);
         return derived.Normalized();
     }
 
@@ -16,5 +18,24 @@ namespace J3ML::Algorithm
     Vector2 Bezier(float t, const Vector2 &p0, const Vector2 &p1, const Vector2 &p2, const Vector2 &p3) {
         return {Bezier(t, p0.x, p1.x, p2.x, p3.x), Bezier(t, p0.y, p1.y, p2.y, p3.y)};
     }
+
+    Vector3 BezierDerivative(float t, const Vector3& p0, const Vector3& p1, const Vector3& p2, const Vector3& p3)
+    {
+        return 3 * Square(1 - t) * (p1 - p0) + 6 * (1 - t) * t * (p2 - p1) + 3 * Square(t) * (p3 - p2);
+    }
+
+    Vector3 BezierNormal(float t, const Vector3& p0, const Vector3& p1, const Vector3& p2, const Vector3& p3)
+    {
+        Vector3 derived = BezierDerivative(t, p0, p1, p2, p3);
+        return derived.Normalized();
+    }
+
+    Vector3 Bezier(float t, const Vector3& p0, const Vector3& p1, const Vector3& p2, const Vector3& p3)
+    {
+        return {
+            Bezier(t, p0.x, p1.x, p2.x, p3.x),
+            Bezier(t, p0.y, p1.y, p2.y, p3.y),
+            Bezier(t, p0.z, p1.z, p2.z, p3.z)};
+    }
 }
 

src/J3ML/Algorithm/RNG.cpp

@@ -36,7 +36,7 @@ namespace J3ML::Algorithm {
 
     u32 RNG::Int()
     {
-        assert(modulus != 0);
+        //assert(modulus != 0);
         /// TODO: Convert to using Shrage's method for approximate factorization (Numerical Recipes in C)
 
         // Currently we cast everything to 64-bit to avoid overflow, which is quite dumb.
@@ -73,8 +73,8 @@ namespace J3ML::Algorithm {
     float RNG::Float() {
         u32 i = ((u32)Int() & 0x007FFFFF /* random mantissa */) | 0x3F800000 /* fixed exponent */;
         // TODO: PROPERLY fix this.
-        auto f = (float)i / 1e31; //ReinterpretAs<float, u32>(i); // f is now in range [1, 2[
-        //f -= 1.f; // Map to range [0, 1[
+        auto f = ReinterpretAs<float, u32>(i); // f is now in range [1, 2[
+        f -= 1.f; // Map to range [0, 1[
         assert(f >= 0.f);
         assert(f < 1.f);
         return f;
@@ -82,11 +82,11 @@ namespace J3ML::Algorithm {
 
 
     float RNG::Float01Incl() {
-        for (int i = 0; i < 100; ++i) {
+        for (int i = 0; i < 1000; ++i) {
             u32 val = (u32)Int() & 0x00FFFFFF;
             if (val > 0x800000)
                 continue;
-            else if (val = 0x800000)
+            else if (val == 0x800000)
                 return 1.f;
             else {
                 val |= 0x3F800000;

src/J3ML/Geometry/AABB.cpp

@@ -12,6 +12,7 @@
 #include <J3ML/Geometry/TriangleMesh.hpp>
 #include <J3ML/Geometry/Polyhedron.hpp>
 #include <J3ML/Algorithm/RNG.hpp>
+#include <J3ML/Geometry/OBB.hpp>
 
 namespace J3ML::Geometry {
 
@@ -320,12 +321,12 @@ namespace J3ML::Geometry {
         return result;
     }
 
-    AABB::AABB(const Vector3 &min, const Vector3 &max) : Shape(), minPoint(min), maxPoint(max)
+    AABB::AABB(const Vector3 &min, const Vector3 &max) : minPoint(min), maxPoint(max)
     {
 
     }
 
-    AABB::AABB() : Shape() {}
+    AABB::AABB() {}
 
     float Max(float a, float b)
     {
@@ -528,6 +529,7 @@ namespace J3ML::Geometry {
                 }
             }
         }
+        return mesh;
     }
 
     std::optional<AABB> AABB::Intersection(const AABB& rhs) const {

src/J3ML/Geometry/AABB2D.cpp

@@ -77,4 +77,20 @@ namespace J3ML::Geometry
         a.maxPoint = maxPoint - pt;
         return a;
     }
+
+    Vector2 AABB2D::Centroid() {
+        return (minPoint + (maxPoint / 2.f));
+    }
+
+    Vector2 AABB2D::CornerPoint(int cornerIndex) {
+        assert(0 <= cornerIndex && cornerIndex <= 3);
+        switch(cornerIndex)
+        {
+            default:
+            case 0: return minPoint;
+            case 1: return {minPoint.x, maxPoint.y};
+            case 2: return {maxPoint.x, minPoint.y};
+            case 3: return maxPoint;
+        }
+    }
 }

src/J3ML/Geometry/Capsule.cpp

@@ -189,6 +189,16 @@ namespace J3ML::Geometry
         return 2.f * Math::Pi * r * LineLength() + 4.f * Math::Pi * r * r;
     }
 
+    OBB Capsule::MinimalEnclosingOBB() const {
+        OBB obb;
+        obb.axis[0] = UpDirection();
+        obb.axis[0].PerpendicularBasis(obb.axis[1], obb.axis[2]);
+        obb.pos = Center();
+        obb.r[0] = Height() * 0.5f;
+        obb.r[1] = r;
+        obb.r[2] = r;
+        return obb;
+    }
 
 
 }

src/J3ML/Geometry/Forward.cpp

@@ -1,4 +1,4 @@
-#include <J3ML/Geometry/Common.hpp>
+#include <J3ML/Geometry/Forward.hpp>
 
 namespace J3ML::Geometry {
 

src/J3ML/Geometry/Frustum.cpp

@@ -1,4 +1,3 @@
-#include <J3ML/Geometry/Common.hpp>
 #include <J3ML/Geometry/Frustum.hpp>
 #include <cmath>
 #include <J3ML/Geometry/AABB.hpp>
@@ -370,6 +369,10 @@ namespace J3ML::Geometry
             Vector3 topSide = front + Math::Tan(verticalFov * 0.5f) * up;
             Vector3 right = WorldRight();
             Vector3 topSideNormal = ((handedness == FrustumHandedness::Right) ? Vector3::Cross(right, topSide) : Vector3::Cross(topSide, right)).Normalized();
+            return Plane(pos, topSideNormal);
+        } else
+        {
+            return Plane(NearPlanePos(0.f, 1.f), up);
         }
     }
 
@@ -378,6 +381,7 @@ namespace J3ML::Geometry
             Vector3 bottomSide = front - Math::Tan(verticalFov * 0.5f) * up;
             Vector3 left = -WorldRight();
             Vector3 bottomSideNormal = ((handedness == FrustumHandedness::Right) ? Vector3::Cross(left, bottomSide) : Vector3::Cross(bottomSide, left)).Normalized();
+            return Plane(pos, bottomSideNormal);
         } else {
             return Plane(NearPlanePos(0.f, -1.f), -up);
         }
@@ -391,6 +395,7 @@ namespace J3ML::Geometry
             right.ScaleToLength(Math::Tan(horizontalFov * 0.5f));
             Vector3 rightSide = front + right;
             Vector3 rightSideNormal = ((handedness == FrustumHandedness::Right) ? Vector3::Cross(rightSide, up) : Vector3::Cross(up, rightSide)).Normalized();
+            return Plane(pos, rightSideNormal);
         } else {
             Vector3 right = WorldRight();
             return Plane(NearPlanePos(1.f,0.f), right.Normalized());
@@ -542,7 +547,7 @@ namespace J3ML::Geometry
 
     Matrix4x4 Frustum::ComputeWorldMatrix() const {
         assert(pos.IsFinite());
-        assert(up.IsNormalized(1e-3f));
+        //assert(up.IsNormalized(1e-3f));
         assert(front.IsNormalized(1e-3f));
         assert(up.IsPerpendicular(front));
 
@@ -572,9 +577,9 @@ namespace J3ML::Geometry
         if (type == FrustumType::Invalid || projectiveSpace == FrustumProjectiveSpace::Invalid)
             return Matrix4x4::NaN;
 
-        assert(type == FrustumType::Perspective || type == FrustumType::Orthographic);
-        assert(projectiveSpace == FrustumProjectiveSpace::GL || projectiveSpace == FrustumProjectiveSpace::D3D);
-        assert(handedness == FrustumHandedness::Left || handedness == FrustumHandedness::Right);
+        //assert(type == FrustumType::Perspective || type == FrustumType::Orthographic);
+        //assert(projectiveSpace == FrustumProjectiveSpace::GL || projectiveSpace == FrustumProjectiveSpace::D3D);
+        //assert(handedness == FrustumHandedness::Left || handedness == FrustumHandedness::Right);
 
         if (type == FrustumType::Perspective) {
             if (projectiveSpace == FrustumProjectiveSpace::GL) {
@@ -603,7 +608,7 @@ namespace J3ML::Geometry
                     return Matrix4x4::D3DOrthoProjLH(nearPlaneDistance, farPlaneDistance, orthographicWidth, orthographicHeight);
             }
         }
-        assert(false && "Not all values of Frustum were initialized properly! Please initialize correctly before calling Frustum::ProjectionMatrix()!");
+        //assert(false && "Not all values of Frustum were initialized properly! Please initialize correctly before calling Frustum::ProjectionMatrix()!");
         return Matrix4x4::NaN;
     }
 

src/J3ML/Geometry/Icosahedron.cpp

@@ -0,0 +1,6 @@
+#include <J3ML/Geometry/Icosahedron.hpp>
+
+namespace J3ML
+{
+
+}

src/J3ML/Geometry/LineSegment.cpp

@@ -1,9 +1,9 @@
 #include <J3ML/Geometry/LineSegment.hpp>
 #include "J3ML/Geometry/Capsule.hpp"
 #include <J3ML/Geometry/Line.hpp>
-#include <J3ML/LinearAlgebra/Vector3.hpp>
-#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
+#include <J3ML/LinearAlgebra.hpp>
 #include <J3ML/Geometry/Plane.hpp>
+#include <J3ML/Geometry/Ray.hpp>
 
 namespace J3ML::Geometry {
 
@@ -272,4 +272,4 @@ namespace J3ML::Geometry {
     LineSegment operator*(const Matrix4x4 &transform, const LineSegment &l) {
         return LineSegment(transform.Mul(l.A), transform.Mul(l.B));
     }
-}
+}

src/J3ML/Geometry/LineSegment2D.cpp

@@ -0,0 +1,165 @@
+#include <J3ML/Geometry/LineSegment2D.hpp>
+
+void Line2DClosestPointLineLine(const Vector2& v0, const Vector2& v10, const Vector2& v2, const Vector2& v32, float& d, float& d2) // TODO: Move to Line2D when that exists
+{
+    // assert(!v10.IsZero());
+    // assert(!v32.IsZero());
+    Vector2 v02 = v0 - v2;
+    float d0232 = v02.Dot(v32);
+    float d3210 = v32.Dot(v10);
+    float d3232 = v32.Dot(v32);
+    // assert(d3232 != 0.f); // Don't call with a zero direction vector.
+    float d0210 = v02.Dot(v10);
+    float d1010 = v10.Dot(v10);
+    float denom = d1010*d3232 - d3210*d3210;
+
+    d = (denom !=0) ? (d0232*d3210 - d0210*d3232) / denom : 0;
+
+    d2 = (d0232 + d * d3210) / d3232;
+}
+
+Geometry::LineSegment2D::LineSegment2D(const Vector2 &a, const Vector2 &b)
+        : A(a), B(b) {}
+
+Vector2 Geometry::LineSegment2D::GetPoint(float d) const {
+    return (1.f - d) * A + d * B;
+}
+
+Vector2 Geometry::LineSegment2D::CenterPoint() const {
+    return (A + B) * 0.5f;
+}
+
+Vector2 Geometry::LineSegment2D::Centroid() const {
+    return CenterPoint();
+}
+
+void Geometry::LineSegment2D::Reverse() {
+    Swap(A, B);
+}
+
+Vector2 Geometry::LineSegment2D::Dir() const {
+    return (B - A).Normalized();
+}
+
+Vector2 Geometry::LineSegment2D::AnyPointFast() const { return A; }
+
+Vector2 Geometry::LineSegment2D::ExtremePoint(const Vector2 &direction) const {
+    return Vector2::Dot(direction, B - A) >= 0.f ? B : A;
+}
+
+Vector2 Geometry::LineSegment2D::ExtremePoint(const Vector2 &direction, float &projectionDistance) const {
+    Vector2 extremePoint = ExtremePoint(direction);
+    projectionDistance = extremePoint.Dot(direction);
+    return extremePoint;
+}
+
+void Geometry::LineSegment2D::Translate(const Vector2 &offset) {
+    A += offset;
+    B += offset;
+}
+
+void Geometry::LineSegment2D::Transform(const Matrix3x3 &transform) {
+    A = transform.Mul(A);
+    B = transform.Mul(B);
+}
+
+void Geometry::LineSegment2D::Transform(const Matrix4x4 &transform) {
+    A = transform.Mul(A);
+    B = transform.Mul(B);
+}
+
+void Geometry::LineSegment2D::Transform(const Quaternion &transform) {
+    //A = transform.Mul(A);
+    //B = transform.Mul(B);
+}
+
+float Geometry::LineSegment2D::Length() const { return A.Distance(B); }
+
+float Geometry::LineSegment2D::LengthSq() const { return A.DistanceSq(B); }
+
+float Geometry::LineSegment2D::LengthSquared() const { return LengthSq(); }
+
+bool Geometry::LineSegment2D::IsFinite() const { return A.IsFinite() && B.IsFinite(); }
+
+bool Geometry::LineSegment2D::Equals(const Geometry::LineSegment2D &rhs, float epsilon) const {
+    return (A.Equals(rhs.A, epsilon) && B.Equals(rhs.B, epsilon) ||
+            A.Equals(rhs.B, epsilon) && B.Equals(rhs.A, epsilon));
+}
+
+bool Geometry::LineSegment2D::Contains(const Vector2 &point, float epsilon) const {
+    return ClosestPoint(point).DistanceSq(point) <= epsilon;
+}
+
+bool Geometry::LineSegment2D::Contains(const Geometry::LineSegment2D &rhs, float epsilon) const {
+    return Contains(rhs.A, epsilon) && Contains(rhs.B, epsilon);
+}
+
+Vector2 Geometry::LineSegment2D::ClosestPoint(const Vector2 &point) const {
+    float d;
+    return ClosestPoint(point, d);
+}
+
+Vector2 Geometry::LineSegment2D::ClosestPoint(const Vector2 &point, float &d) const {
+    Vector2 dir = B - A;
+    d = J3ML::Math::Clamp01(Vector2::Dot(point - A, dir) / dir.LengthSquared());
+    return A + d * dir;
+}
+
+Vector2 Geometry::LineSegment2D::ClosestPoint(const Geometry::LineSegment2D &other) const {
+    float d, d2;
+    return ClosestPoint(other, d, d2);
+}
+
+Vector2 Geometry::LineSegment2D::ClosestPoint(const Geometry::LineSegment2D &other, float &d) const {
+    float d2; return ClosestPoint(other, d, d2);
+}
+
+Vector2 Geometry::LineSegment2D::ClosestPoint(const Geometry::LineSegment2D &other, float &d, float &d2) const {
+    Vector2 dir = B - A;
+    Line2DClosestPointLineLine(A, B - A, other.A, other.B - other.A, d, d2);
+
+
+    return Vector2::Zero;
+}
+
+float Geometry::LineSegment2D::Distance(const Vector2 &point) const { float d; return Distance(point, d); }
+
+float Geometry::LineSegment2D::Distance(const Vector2 &point, float &d) const {
+    /// See Christer Ericson's Real-Time Collision Detection, p. 130.
+    Vector2 closestPoint = ClosestPoint(point, d);
+    return closestPoint.Distance(point);
+}
+
+float Geometry::LineSegment2D::Distance(const Geometry::LineSegment2D &other) const { float d, d2; return Distance(other, d, d2);}
+
+float Geometry::LineSegment2D::Distance(const Geometry::LineSegment2D &other, float &d) const { float d2; return Distance(other, d, d2); }
+
+float Geometry::LineSegment2D::Distance(const Geometry::LineSegment2D &other, float &d, float &d2) const {
+    ClosestPoint(other, d, d2);
+    return GetPoint(d).Distance(other.GetPoint(d2));
+}
+
+float Geometry::LineSegment2D::DistanceSq(const Geometry::LineSegment2D &other) const {
+    float d, d2;
+    ClosestPoint(other, d, d2);
+    return GetPoint(d).DistanceSq(other.GetPoint(d2));
+}
+
+float Geometry::LineSegment2D::DistanceSq(const Vector2 &point) const {
+    float d;
+    /// See Christer Ericson's Real-Time Collision Detection, p.130.
+    Vector2 closestPoint = ClosestPoint(point, d);
+    return closestPoint.DistanceSq(point);
+}
+
+bool Geometry::LineSegment2D::Intersects(const Geometry::LineSegment2D &lineSegment, float epsilon) const {
+    return Distance(lineSegment) <= epsilon;
+}
+
+void Geometry::LineSegment2D::ProjectToAxis(const Vector2 &direction, float &outMin, float &outMax) const {
+    outMin = Vector2::Dot(direction, A);
+    outMax = Vector2::Dot(direction, B);
+
+    if (outMax < outMin)
+        Swap(outMin, outMax);
+}

src/J3ML/Geometry/OBB.cpp

@@ -1,8 +1,8 @@
-#include <J3ML/Geometry/Shape.hpp>
 #include <J3ML/Geometry/AABB.hpp>
 #include <J3ML/Geometry/LineSegment.hpp>
 #include <J3ML/Geometry/Polyhedron.hpp>
 #include <J3ML/Geometry/OBB.hpp>
+#include <J3ML/Geometry/Plane.hpp>
 #include <J3ML/Geometry/Sphere.hpp>
 
 namespace J3ML::Geometry {

src/J3ML/Geometry/PBVolume.cpp

@@ -0,0 +1 @@
+#include <J3ML/Geometry/PBVolume.hpp>

src/J3ML/Geometry/Plane.cpp

@@ -2,8 +2,13 @@
 #include <J3ML/Geometry/AABB.hpp>
 #include <J3ML/Geometry/OBB.hpp>
 #include <J3ML/Geometry/Capsule.hpp>
+#include <J3ML/Geometry/Frustum.hpp>
+#include <J3ML/Geometry/Line.hpp>
 #include <J3ML/Geometry/Polygon.hpp>
+#include <J3ML/Geometry/Polyhedron.hpp>
+#include <J3ML/Geometry/Ray.hpp>
 #include <J3ML/Geometry/Sphere.hpp>
+#include <J3ML/Geometry/Triangle.hpp>
 
 namespace J3ML::Geometry
 {
@@ -138,7 +143,7 @@ namespace J3ML::Geometry
 
     Plane::Plane(const Vector3 &pos, const Vector3 &norm) : Position(pos), Normal(norm) {}
 
-    bool Plane::Intersects(J3ML::Geometry::Ray ray, float *dist) const {
+    bool Plane::Intersects(Ray ray, float *dist) const {
         float t;
         bool success = IntersectLinePlane(Normal, this->distance, ray.Origin, ray.Direction, t);
         if (dist)

src/J3ML/Geometry/Polyhedron.cpp

@@ -8,6 +8,10 @@
 #include <J3ML/Geometry/Capsule.hpp>
 #include <set>
 #include <cfloat>
+#include <J3ML/Geometry/Frustum.hpp>
+#include <J3ML/Geometry/OBB.hpp>
+#include <J3ML/Geometry/Plane.hpp>
+#include <J3ML/Geometry/Sphere.hpp>
 
 namespace J3ML::Geometry
 {
@@ -159,7 +163,7 @@ namespace J3ML::Geometry
 
             Vector3 Dir = ((Vector3)v[f[j].v[0]] + (Vector3)v[f[j].v[1]] + (Vector3)v[f[j].v[2]]) * 0.333333333333f - point;
 
-            //if (Dir.Normalize() <= 0.f)
+            //if (Dir.Normalized() <= 0.f)
                 //continue;
             Ray r(Vector3(), Dir);
 
@@ -368,7 +372,7 @@ namespace J3ML::Geometry
             Vector4 b = Vector4(poly.v[face.v[1]], 1.f);
             Vector4 c = Vector4(poly.v[face.v[2]], 1.f);
             Vector4 normal = (b-a).Cross(c-a);
-            normal.Normalize();
+            normal.Normalized();
             return normal;
 //		return ((vec)v[face.v[1]]-(vec)v[face.v[0]]).Cross((vec)v[face.v[2]]-(vec)v[face.v[0]]).Normalized();
         }
@@ -386,7 +390,7 @@ namespace J3ML::Geometry
                 normal.z += (double(poly.v[v0].x) - poly.v[v1].x) * (double(poly.v[v0].y) + poly.v[v1].y); // Project on xy
                 v0 = v1;
             }
-            normal.Normalize();
+            normal.Normalized();
             return normal;
 #if 0
             cv bestNormal;
@@ -397,7 +401,7 @@ namespace J3ML::Geometry
 		{
 			cv c = poly.v[face.v[i]];
 			cv normal = (c-b).Cross(a-b);
-			float len = normal.Normalize();
+			float len = normal.Normalized();
 			if (len > bestLen)
 			{
 				bestLen = len;
@@ -441,7 +445,7 @@ namespace J3ML::Geometry
 			cv edge = cv(poly.v[i]) - cv(poly.v[bestV0]);
 			edge.Normalize();
 			cv normal = bestEdge.Cross(edge);
-			cs len = normal.Normalize();
+			cs len = normal.Normalized();
 			if (len > bestLen)
 			{
 				bestLen = len;

src/J3ML/Geometry/Ray.cpp

@@ -1,6 +1,10 @@
+#include <J3ML/Geometry/AABB.hpp>
+#include <J3ML/Geometry/Line.hpp>
+#include <J3ML/Geometry/LineSegment.hpp>
+#include <J3ML/Geometry/Plane.hpp>
 #include <J3ML/Geometry/Ray.hpp>
 #include <J3ML/Geometry/Sphere.hpp>
-
+#include <ostream>
 
 namespace J3ML::Geometry
 {
@@ -17,6 +21,17 @@ namespace J3ML::Geometry
             outMin = -INFINITY;
     }
 
+    std::string Ray::ToString() const
+    {
+        return std::format("Ray(origin:[{}], direction:[{}])", Origin.ToString(), Direction.ToString());
+    }
+
+    std::ostream& operator<<(std::ostream& o, const Ray& ray)
+    {
+        o << ray.ToString();
+        return o;
+    }
+
     RaycastResult Ray::Cast(const Sphere &target, float maxDistance)
     {
         Vector3 p0 = this->Origin;
@@ -54,7 +69,6 @@ namespace J3ML::Geometry
             intersection,
             normal,
             true,
-            (Shape *) &target
         };
     }
 
@@ -101,7 +115,6 @@ namespace J3ML::Geometry
             intersection,
             normal,
             true,
-            (Shape*)&target
         };
     }
 
@@ -129,7 +142,6 @@ namespace J3ML::Geometry
                         intersection,
                         normal,
                         true,
-                        (Shape*) &target
                     };
 
         return RaycastResult::NoHit();
@@ -199,7 +211,7 @@ namespace J3ML::Geometry
     }
 
     Vector3 Ray::GetPoint(float distance) const {
-        assert(Direction.IsNormalized());
+        //assert(Direction.IsNormalized());
         return Origin + distance * Direction;
     }
 }

src/J3ML/Geometry/Sphere.cpp

@@ -1,4 +1,6 @@
+#include <J3ML/Geometry/LineSegment.hpp>
 #include <J3ML/Geometry/Sphere.hpp>
+#include <J3ML/Geometry/TriangleMesh.hpp>
 
 namespace J3ML::Geometry
 {
@@ -7,6 +9,107 @@ namespace J3ML::Geometry
         return Contains(lineseg.A) && Contains(lineseg.B);
     }
 
+    TriangleMesh Sphere::GenerateUVSphere(int subdivisions) const
+    {
+        // http://www.songho.ca/opengl/gl_sphere.html
+
+        TriangleMesh mesh;
+
+        float x, y, z, xy;                           // Vertex Position
+        float nx, ny, nz, lengthInv = 1.f / Radius;  // Vertex Normal
+        float s, t;                                  // Vertex TexCoord
+
+        int sectorCount = subdivisions;
+        int stackCount = subdivisions;
+
+        float sectorStep = 2.f * Math::Pi / sectorCount;
+        float stackStep = Math::Pi / stackCount;
+        float sectorAngle, stackAngle;
+
+        for (int i = 0; i <= stackCount; ++i)
+        {
+            stackAngle = Math::Pi / 2.f - i * stackStep;   // starting from pi/2 to -pi/2
+            xy = Radius * Math::Cos(stackAngle);           // r * cos(u)
+            z = Radius * Math::Sin(stackAngle);            // r * sin(u)
+
+            // add (sectorCount + 1) vertices per stack
+            // first and last vertices have same position and normal, but different tex coords
+            for (int j = 0; j <= sectorCount; ++j)
+            {
+                sectorAngle = j * sectorStep;   // starting from 0 to 2pi
+
+                // vertex position (x, y, z)
+                x = xy * Math::Cos(sectorAngle);
+                y = xy * Math::Sin(sectorAngle);
+
+                Vector3 vertex = {x, y, z};
+
+                mesh.Vertices.push_back(vertex);
+
+                // normalized vertex normal (nx, ny, nz)
+                nx = x * lengthInv;
+                ny = y * lengthInv;
+                nz = z * lengthInv;
+
+                Vector3 normal = {nx, ny, nz};
+
+                mesh.Normals.push_back(normal);
+
+                // vertex tex coord (s, t) range between [0, 1]
+                s = (float)j / sectorCount;
+                t = (float)i / stackCount;
+
+                Vector2 TexCoords = {s, t};
+
+                mesh.UVs.push_back(normal);
+            }
+        }
+        return mesh;
+    }
+
+    TriangleMesh Sphere::GenerateIcososphere() const
+    {
+
+        // Generate 12 vertices of an icosahedron for a given radius.
+
+        const float h_angle = Math::Pi / 180.f * 72.f; // 72 degree = 360 / 5;
+        const float v_angle = Math::Atan(1.f / 2.f);
+
+        TriangleMesh mesh;
+
+        int i1, i2;
+        float z, xy;
+        float hAngle1 = -Math::Pi / 2.f - h_angle / 2.f;
+        float hAngle2 = -Math::Pi / 2;
+
+        // the first top vertex at (0,0,r)
+        Vector3 top_vertex = {0, 0, Radius};
+
+        // compute 10 vertices at 1st and 2nd rows
+        for (int i = 1; i <= 5; ++i)
+        {
+            i1 = i * 3;       // index for 1st row
+            i2 = (i + 5) * 3; // index for 2nd row
+
+            z  = Radius * Math::Sin(v_angle);   // elevation
+            xy = Radius * Math::Cos(v_angle);   // length on XY plane
+
+            Vector3 vert_0 = {xy * Math::Cos(hAngle1), xy * Math::Sin(hAngle1), z};
+            Vector3 vert_1 = {xy * Math::Cos(hAngle2), xy * Math::Sin(hAngle2), -z};
+
+            // next horizontal angles
+            hAngle1 += h_angle;
+            hAngle2 += h_angle;
+        }
+
+        // the last bottom vertex at (0, 0, -r)
+        i1 = 11 * 3;
+        Vector3 last_vertex = {0,0, -Radius};
+
+        return mesh;
+    }
+
+
     void Sphere::ProjectToAxis(const Vector3 &direction, float &outMin, float &outMax) const
     {
         float d = Vector3::Dot(direction, Position);
@@ -43,11 +146,11 @@ namespace J3ML::Geometry
     }
 
     float Sphere::Volume() const {
-        return 4.f * M_PI * Cube(Radius) / 3.f;
+        return 4.f * Math::Pi * Cube(Radius) / 3.f;
     }
 
     float Sphere::SurfaceArea() const {
-        return 4.f * M_PI * Cube(Radius) / 3.f;
+        return 4.f * Math::Pi * Cube(Radius) / 3.f;
     }
 
     bool Sphere::IsFinite() const {
@@ -108,4 +211,4 @@ namespace J3ML::Geometry
     }
 
 
-}
+}

src/J3ML/Geometry/Triangle.cpp

@@ -4,6 +4,7 @@
 #include <J3ML/Geometry/LineSegment.hpp>
 #include <J3ML/Geometry/Line.hpp>
 #include <J3ML/Geometry/Capsule.hpp>
+#include <J3ML/Geometry/Plane.hpp>
 
 namespace J3ML::Geometry
 {
@@ -44,10 +45,9 @@ namespace J3ML::Geometry
     }
 
     Triangle Triangle::Translated(const Vector3& translation) const {
-        return {
-            V0 + translation,
-            V1 + translation,
-            V2 + translation
+        return {V0 + translation,
+                V1 + translation,
+                V2 + translation
         };
     }
 

src/J3ML/Geometry/Triangle2D.cpp

@@ -0,0 +1,136 @@
+#include <J3ML/Geometry/Triangle2D.hpp>
+#include <J3ML/LinearAlgebra/Vector2.hpp>
+#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
+#include <cfloat>
+
+namespace J3ML::Geometry
+{
+
+    Triangle2D::Triangle2D(const Vector2& a, const Vector2& b, const Vector2& c) : A(a), B(b), C(c) {}
+
+    void Triangle2D::Translate(const Vector2& offset)
+    {
+        A += offset;
+        B += offset;
+        C += offset;
+    }
+
+    void Triangle2D::Transform(const Matrix3x3& transform)
+    {
+        A = Mul2D(transform, A);
+        B = Mul2D(transform, B);
+        C = Mul2D(transform, C);
+    }
+
+    void Triangle2D::Transform(const Matrix4x4& transform)
+    {
+        A = MulPos2D(transform, A);
+        B = MulPos2D(transform, B);
+        C = MulPos2D(transform, C);
+    }
+
+    /// Implementation from Crister Ericson's Real-Time Collision Detection, pp. 51-52.
+    inline float TriArea2D(float x1, float y1, float x2, float y2, float x3, float y3)
+    {
+        return (x1-x2)*(y2-y3) - (x2-x3)*(y1-y2);
+    }
+
+    Vector3 Triangle2D::BarycentricUVW(const Vector2& point) const
+    {
+        // Implementation from Crister Ericson's Real-Time Collision Detection, pp. 51-52, adapted for 2D.
+        float nu, nv;
+
+        nu = TriArea2D(point.x, point.y, B.x, B.y, C.x, C.y);
+        nv = TriArea2D(point.x, point.y, C.x, C.y, A.x, A.y);
+        float u = nu;
+        float v = nv;
+        float w = 1.f - u - v;
+        return Vector3(u,v,w);
+    }
+
+    Vector2 Triangle2D::BarycentricUV(const Vector2& point) const
+    {
+        Vector3 uvw = BarycentricUVW(point);
+        return Vector2(uvw.y, uvw.z);
+    }
+
+    bool Triangle2D::BarycentricInsideTriangle(const Vector3& barycentric)
+    {
+        return barycentric.x >= 0.f && barycentric.y >= 0 && barycentric.z >= 0.f &&
+            Math::EqualAbs(barycentric.x + barycentric.y + barycentric.z, 1.f);
+    }
+
+    Vector2 Triangle2D::Point(float u, float v) const
+    {
+        // In case the triangle is far away from the origin but is small, the elements of 'a' will have large magnitudes,
+        // and the elements of (b-a) and (c-a) will be much smaller quantities. Therefore, be extra careful with the
+        // parenthesis and first sum the small floats together before adding it to the large one.
+        return A + ((B-A) * u + (C-A) * v);
+    }
+
+    Vector2 Triangle2D::Point(float u, float v, float w) const
+    {
+        return u * A + v * B + w * C;
+    }
+
+    Vector2 Triangle2D::Point(const Vector3 &uvw) const
+    {
+        return Point(uvw.x, uvw.y, uvw.z);
+    }
+
+    Vector2 Triangle2D::Point(const Vector2 &uv) const
+    {
+        return Point(uv.x, uv.y);
+    }
+
+    Vector2 Triangle2D::Centroid() const
+    {
+        return (A + B + C) * (1.f/3.f);
+    }
+
+    float Triangle2D::Perimeter() const
+    {
+        return A.Distance(B) + B.Distance(C) + C.Distance(A);
+    }
+
+    Vector2 Triangle2D::Vertex(int i) const
+    {
+        assert(0 <= i);
+        assert(i <= 2);
+        if (i == 0)
+            return A;
+        else if (i == 1)
+            return B;
+        else if (i == 2)
+            return C;
+        else
+            return Vector2::NaN;
+    }
+
+
+    Vector2 Triangle2D::ExtremePoint(const Vector2& direction) const
+    {
+        Vector2 mostExtreme = Vector2::NaN;
+        float mostExtremeDist = -FLT_MAX;
+        for (int i = 0; i < 3; ++i)
+        {
+            Vector2 pt = Vertex(i);
+            float d = Vector2::Dot(direction, pt);
+            if (d > mostExtremeDist)
+            {
+                mostExtremeDist = d;
+                mostExtreme = pt;
+            }
+        }
+        return mostExtreme;
+    }
+
+    //TODO must return on MSVC
+    /*
+    Vector2 Triangle2D::ExtremePoint(const Vector2& direction, float& projectionDistance) const
+    {
+        Vector2 extremePoint = ExtremePoint(direction);
+        projectionDistance = extremePoint.Dot(direction);
+    }
+    */
+}

src/J3ML/J3ML.cpp

@@ -13,6 +13,8 @@
 #include <strstream>
 #include "J3ML/J3ML.hpp"
 
+#include <sstream>
+
 
 float PowUInt(float base, u32 exponent)
 {

src/J3ML/LinearAlgebra/AxisAngle.cpp

@@ -2,60 +2,22 @@
 #include <J3ML/LinearAlgebra/Quaternion.hpp>
 namespace J3ML::LinearAlgebra {
 
-    AxisAngle::AxisAngle() : axis(Vector3::Zero) {}
+    AxisAngle::AxisAngle() : axis(Vector3::Zero), angle(0) {}
 
-    AxisAngle::AxisAngle(const Vector3 &axis, float angle) : axis(axis), angle(angle) {}
+    AxisAngle::AxisAngle(const Vector3& axis, float angle) : axis(axis), angle(angle) {}
 
-    Quaternion AxisAngle::ToQuaternion() const {
-        return {
-                axis.x * std::sin(angle/2),
-                axis.y * std::sin(angle/2),
-                axis.z * std::sin(angle/2),
-                std::cos(angle/2)
-        };
+
+    AxisAngle::AxisAngle(const Quaternion& rhs) {
+        float halfAngle = std::acos(rhs.w);
+        angle = halfAngle * 2.f;
+        float reciprocalSinAngle = 1.f / std::sqrt(1.f - rhs.w*rhs.w);
+
+        axis = { rhs.x*reciprocalSinAngle, rhs.y*reciprocalSinAngle, rhs.z*reciprocalSinAngle };
     }
 
-    AxisAngle::AxisAngle(const Quaternion &q) {
-        auto theta = std::acos(q.w) * 2.f;
-        auto ax = q.x / std::sin(std::acos(theta));
-        auto ay = q.y / std::sin(std::acos(theta));
-        auto az = q.z / std::sin(std::acos(theta));
-    }
-
-    AxisAngle::AxisAngle(const EulerAngle &e) {
-
-        // Assuming the angles are in radians
-
-        float heading = e.pitch;
-        float attitude = e.yaw;
-        float bank = e.roll;
-
-        float c1 = std::cos(heading / 2.f);
-        float s1 = std::sin(heading / 2.f);
-        float c2 = std::cos(attitude / 2.f);
-        float s2 = std::sin(attitude / 2.f);
-        float c3 = std::cos(bank / 2.f);
-        float s3 = std::sin(bank / 2.f);
-
-        float w = c1*c2*c3 - s1*s2*s3;
-        float x = c1*c2*c3 + s1*s2*s3;
-        float y = s1*c2*c3 + c1*s2*s3;
-        float z = c1*s2*c3 - s1*c2*s3;
-
-        angle = 2.f * std::acos(w);
-
-        double norm = x*x + y*y + z*z;
-        if (norm < 0.001) { // when all euler angles are zero angle=0, so
-            // we can set axis to anything to avoid divide by zero
-            x = 1;
-            y = z = 0;
-        } else {
-            norm = std::sqrt(norm);
-            x /= norm;
-            y /= norm;
-            z /= norm;
-        }
-
-        axis = {x, y, z};
+    AxisAngle::AxisAngle(const EulerAngleXYZ& e) {
+        auto a = AxisAngle(Quaternion(e));
+        axis = a.axis;
+        angle = a.angle;
     }
 }

src/J3ML/LinearAlgebra/DirectionVector.cpp

@@ -0,0 +1,38 @@
+#include <J3ML/LinearAlgebra/DirectionVector.hpp>
+#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
+
+DirectionVectorRH::DirectionVectorRH(const Vector3& rhs) {
+    x = rhs.x;
+    y = rhs.y;
+    z = rhs.z;
+}
+
+DirectionVectorRH DirectionVectorRH::Forward(const Matrix3x3& rhs) {
+    return DirectionVectorRH(rhs.Col(2));
+}
+
+DirectionVectorRH DirectionVectorRH::Backward(const Matrix3x3& rhs) {
+    return DirectionVectorRH(-rhs.Col(2));
+}
+
+DirectionVectorRH DirectionVectorRH::Left(const Matrix3x3& rhs) {
+    return DirectionVectorRH(-rhs.Col(0));
+}
+
+DirectionVectorRH DirectionVectorRH::Right(const Matrix3x3& rhs) {
+    return DirectionVectorRH(rhs.Col(0));
+}
+
+DirectionVectorRH DirectionVectorRH::Up(const Matrix3x3 &rhs) {
+    return DirectionVectorRH(rhs.Col(1));
+}
+
+DirectionVectorRH DirectionVectorRH::Down(const Matrix3x3& rhs) {
+    return DirectionVectorRH(-rhs.Col(1));
+}
+
+
+
+
+
+

src/J3ML/LinearAlgebra/EulerAngle.cpp

@@ -1,147 +1,45 @@
 #include <J3ML/LinearAlgebra/EulerAngle.hpp>
+#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
 #include <cmath>
 #include <algorithm>
 
 
 namespace J3ML::LinearAlgebra {
-    EulerAngle::EulerAngle(float pitch, float yaw, float roll): pitch(pitch), yaw(yaw), roll(roll)
-    {}
 
-    float EulerAngle::GetPitch(float pitch_limit) const
-    { return std::clamp( std::remainderf(pitch,360.f), -pitch_limit, pitch_limit); }
-
-    float EulerAngle::GetYaw(float yaw_limit) const
-    { return std::clamp(std::remainderf(yaw, 360.f), -yaw_limit, yaw_limit);   }
-
-    float EulerAngle::GetRoll(float pitch_limit) const
-    { return std::clamp( std::remainderf(pitch,360.f), -pitch_limit, pitch_limit); }
-
-    bool EulerAngle::operator==(const EulerAngle& a) const
-    {
-        return (pitch == a.pitch) && (yaw == a.yaw) && (roll == a.roll);
+    EulerAngleXYZ::EulerAngleXYZ(float roll, float pitch, float yaw) {
+        this->roll = roll;
+        this->pitch = pitch;
+        this->yaw = yaw;
     }
 
-    void EulerAngle::clamp()
-    {
-        if (this->pitch > 89.0f)
-            this->pitch = 89.0f;
-        if (this->pitch <= -89.0f)
-            this->pitch = -89.0f;
-        //TODO: Make this entirely seamless by getting the amount they rotated passed -180 and +180 by.
-        if (this->yaw <= -180.0f)
-            this->yaw = 180.0f;
-        if (this->yaw >= 180.01f)
-            this->yaw = -179.9f;
-        if (this->roll >= 360.0f)
-            this->roll = 0.0;
-        if (this->roll <= -360.0f)
-            this->roll = 0.0;
+    EulerAngleXYZ::EulerAngleXYZ(const AxisAngle& rhs) {
+        *this = EulerAngleXYZ(Quaternion(rhs));
     }
 
-    EulerAngle EulerAngle::movementAngle() const
-    {
-        EulerAngle a;
-        a.pitch = (cos(Math::Radians(yaw)) * cos(Math::Radians(pitch)));
-        a.yaw = -sin(Math::Radians(pitch));
-        a.roll = (sin(Math::Radians(yaw)) * cos(Math::Radians(pitch)));
-        return a;
+    EulerAngleXYZ::EulerAngleXYZ(const Quaternion& q) {
+        float sy = 2 * q.x * q.z + 2 * q.y * q.w;
+        bool gimbal_lock = std::abs(sy) > 0.99999f;
+
+        if (!gimbal_lock)
+            roll = Math::Degrees(std::atan2(-(2 * q.y * q.z - 2 * q.x * q.w),2 * q.w * q.w + 2 * q.z * q.z - 1));
+        else
+            roll = Math::Degrees(std::atan2(2 * q.y * q.z + 2 * q.x * q.w,2 * q.w * q.w + 2 * q.y * q.y - 1));
+
+        pitch = Math::Degrees(std::asin(sy));
+
+        if (!gimbal_lock)
+            yaw = Math::Degrees(std::atan2(-(2 * q.x * q.y - 2 * q.z * q.w),2 * q.w * q.w + 2 * q.x * q.x - 1));
+        else
+            yaw = 0;
     }
 
-    EulerAngle::EulerAngle() : pitch(0), yaw(0), roll(0) {}
+    EulerAngleXYZ::EulerAngleXYZ(const Matrix3x3& rhs) {
+        auto m = rhs.Transposed();
+        auto sy = m.At(0, 2);
+        auto unlocked = std::abs(sy) < 0.99999f;
 
-    EulerAngle::EulerAngle(const AxisAngle &rhs) {
-
-        float x = rhs.axis.x;
-        float y = rhs.axis.y;
-        float z = rhs.axis.z;
-        float angle = rhs.angle;
-
-        double s = std::sin(rhs.angle);
-
-        double c = std::cos(rhs.angle);
-
-        double t = 1-c;
-
-        // if axis is not already normalized then uncomment this
-
-        // double magnitude = std::sqrt(x*x + y*y + z*z);
-        // if (magnitude == 0) throw error;
-        // x /= magnitude;
-        // y /= magnitude;
-        // z /= magnitude;
-
-        if ((x*y*t + z*s) > 0.998) { // North pole singularity detected
-            pitch = 2 * std::atan2(x * std::sin(angle/2.f), std::cos(angle/2.f));
-            yaw = M_PI / 2.f;
-            roll = 0;
-            return;
-        }
-
-        if ((x*y*t + z*s) < -0.998) { // South pole singularity detected
-            pitch = -2 * std::atan2(x * std::sin(angle/2.f), std::cos(angle/2.f));
-            yaw = -M_PI / 2.f;
-            roll = 0;
-            return;
-        }
-
-        pitch = std::atan2(y * s-x * z * t, 1 - (y*y + z*z) * t);
-        yaw = std::asin(x * y * t + z * s);
-        roll = std::atan2(x * s - y * z * t, 1 - (x*x + z*z) * t);
-    }
-
-    AxisAngle EulerAngle::ToAxisAngle() const {
-        auto c1 = std::cos(yaw / 2);
-        auto c2 = std::cos(pitch / 2);
-        auto c3 = std::cos(roll / 2);
-        auto s1 = std::sin(yaw / 2);
-        auto s2 = std::sin(pitch / 2);
-        auto s3 = std::sin(roll / 2);
-
-        auto angle = 2 * std::acos(c1*c2*c3 - s1*s2*s3);
-
-        auto x = s1*s2*c3 + c1*c2*s3;
-        auto y = s1*c2*c3 + c1*s2*s3;
-        auto z = c1*s2*c3 - s1*c2*s3;
-
-        // todo: normalize?
-        // sqrt(x^2 + y^2 + z^2) = sqrt((s1 s2 c3 +c1 c2 s3)^2+(s1 c2 c3 + c1 s2 s3)^2+(c1 s2 c3 - s1 c2 s3)^2)
-
-        return {{x,y,z}, angle};
-    }
-
-    Quaternion EulerAngle::ToQuaternion() const {
-        auto c1 = std::cos(yaw / 2);
-        auto c2 = std::cos(pitch / 2);
-        auto c3 = std::cos(roll / 2);
-        auto s1 = std::sin(yaw / 2);
-        auto s2 = std::sin(pitch / 2);
-        auto s3 = std::sin(roll / 2);
-
-        auto w = c1*c2*c3 - s1*s2*s3;
-        auto x = s1*s2*c3 + c1*c2*s3;
-        auto y = s1*c2*c3 + c1*s2*s3;
-        auto z = c1*s2*c3 - s1*c2*s3;
-
-        return {w,x,y,z};
-    }
-
-    EulerAngle::EulerAngle(const Quaternion &rhs) {
-        double test = rhs.x * rhs.y + rhs.z * rhs.w;
-        if (test > 0.499) { // Singularity at north pole
-            pitch = 2 * std::atan2(rhs.x, rhs.w);
-            yaw = M_PI / 2.f;
-            roll = 0;
-            return;
-        }
-
-        if (test < -0.499) { // Singularity at south pole
-            pitch = -2 * std::atan2(rhs.x, rhs.y);
-            yaw = - M_PI / 2.f;
-            roll = 0;
-            return;
-        }
-        float sqx = rhs.x * rhs.x;
-        float sqy = rhs.y * rhs.y;
-        float sqz = rhs.z * rhs.z;
+        roll = Math::Degrees(unlocked ? std::atan2(-m.At(1, 2), m.At(2, 2)) : std::atan2(m.At(2, 1), m.At(1, 1)));
+        pitch = Math::Degrees(std::asin(sy));
+        yaw = Math::Degrees(unlocked ? std::atan2(-m.At(0, 1), m.At(0, 0)) : 0);
     }
 }

src/J3ML/LinearAlgebra/Matrix3x3.cpp

@@ -109,27 +109,8 @@ namespace J3ML::LinearAlgebra {
         //this->elems[2][2] = r3.z;
     }
 
-    Matrix3x3::Matrix3x3(const Quaternion &orientation) {
-        SetRotatePart(orientation);
-    }
-
-    Matrix3x3::Matrix3x3(const EulerAngle &orientation) {
-        auto sa = std::sin(orientation.pitch);
-        auto ca = std::cos(orientation.pitch);
-        auto sb = std::sin(orientation.roll);
-        auto cb = std::cos(orientation.roll);
-        auto sh = std::sin(orientation.yaw);
-        auto ch = std::cos(orientation.yaw);
-
-        At(0, 0) = ch*ca;
-        At(0, 1) = -ch*sa*cb + sh*sh;
-        At(0, 2) = ch*sa*sb + sh*cb;
-        At(1, 0) = sa;
-        At(1, 1) = ca*cb;
-        At(1, 2) = -ca*cb;
-        At(2, 0) = -sh*ca;
-        At(2, 1) = sh*sa*cb + ch*sb;
-        At(2, 2) = -sh*sa*sb + ch*cb;
+    Matrix3x3::Matrix3x3(const Quaternion& orientation) {
+        *this = Matrix3x3(EulerAngleXYZ(orientation));
     }
 
     float Matrix3x3::Determinant() const {
@@ -207,28 +188,7 @@ namespace J3ML::LinearAlgebra {
         };
     }
 
-    Quaternion Matrix3x3::ToQuat() const {
-        auto m00 = At(0,0);
-        auto m01 = At(0, 1);
-        auto m02 = At(0, 2);
-        auto m10 = At(1,0);
-        auto m11 = At(1, 1);
-        auto m12 = At(1, 2);
-        auto m20 = At(2,0);
-        auto m21 = At(2, 1);
-        auto m22 = At(2, 2);
-
-        auto w = std::sqrt(1.f + m00 + m11 + m22) / 2.f;
-        float w4 = (4.f * w);
-        return {
-                (m21 - m12) / w4,
-                (m02 - m20) / w4,
-                (m10 - m01) / w4,
-                w
-        };
-    }
-
-    void Matrix3x3::SetRotatePart(const Vector3 &a, float angle) {
+    void Matrix3x3::SetRotatePart(const Vector3& a, float angle) {
         float s = std::sin(angle);
         float c = std::cos(angle);
 
@@ -371,9 +331,9 @@ namespace J3ML::LinearAlgebra {
         return m;
     }
 
-    Matrix3x3 Matrix3x3::FromScale(float sx, float sy, float sz) {
+    Matrix3x3 Matrix3x3::FromScale(float sin_roll, float sy, float sz) {
         Matrix3x3 m;
-        m.At(0,0) = sx;
+        m.At(0,0) = sin_roll;
         m.At(1,1) = sy;
         m.At(2,2) = sz;
         return m;
@@ -992,25 +952,12 @@ namespace J3ML::LinearAlgebra {
 
     bool Matrix3x3::TryConvertToQuat(Quaternion &q) const {
         if (IsColOrthogonal() && HasUnitaryScale() && !HasNegativeScale()) {
-            q = ToQuat();
+            q = Quaternion(*this);
             return true;
         }
         return false;
     }
 
-    EulerAngle Matrix3x3::ToEulerAngle() const {
-        auto heading = std::atan2(-At(2, 0), At(0, 0));
-        auto attitude = std::asin(At(1, 0));
-        auto bank = std::atan2(-At(1,2), At(1,1));
-        if (At(1, 0) == 1 || At(1, 0) == -1) // North Pole || South Pole
-        {
-            heading = std::atan2(At(0, 2), At(2,2));
-            bank = 0;
-        }
-
-        return {attitude, heading, bank};
-    }
-
     void Matrix3x3::BatchTransform(Vector3 *pointArray, int numPoints, int stride) const {
         assert(pointArray || numPoints == 0);
         assert(stride >= (int)sizeof(Vector3));
@@ -1119,6 +1066,25 @@ namespace J3ML::LinearAlgebra {
         return m;
     }
 
+    Matrix3x3::Matrix3x3(const EulerAngleXYZ& e) {
+        float cos_roll = std::cos(Math::Radians(e.roll));
+        float sin_roll = std::sin(Math::Radians(e.roll));
+        float cos_pitch = std::cos(Math::Radians(e.pitch));
+        float sin_pitch = std::sin(Math::Radians(e.pitch));
+        float cos_yaw = std::cos(Math::Radians(e.yaw));
+        float sin_yaw = std::sin(Math::Radians(e.yaw));
+
+        Matrix3x3 m;
+        m.SetRow(0, Vector3(cos_pitch * cos_yaw, sin_roll * sin_pitch *cos_yaw + cos_roll * sin_yaw, -cos_roll * sin_pitch *cos_yaw + sin_roll * sin_yaw));
+        m.SetRow(1, Vector3(-cos_pitch * sin_yaw, -sin_roll * sin_pitch * sin_yaw + cos_roll *cos_yaw, cos_roll * sin_pitch * sin_yaw + sin_roll *cos_yaw));
+        m.SetRow(2, Vector3(sin_pitch, -sin_roll * cos_pitch, cos_roll * cos_pitch));
+        *this = m;
+    }
+
+    Matrix3x3::Matrix3x3(const AxisAngle& orientation) {
+        *this = Matrix3x3(Quaternion(orientation));
+    }
+
 
 }
 

src/J3ML/LinearAlgebra/Matrix4x4.cpp

@@ -86,10 +86,6 @@ namespace J3ML::LinearAlgebra {
         Set3x3Part(Matrix3x3(orientation));
     }
 
-    Matrix4x4::Matrix4x4(const EulerAngle &orientation) {
-        Set3x3Part(Matrix3x3(orientation));
-    }
-
     void Matrix4x4::SetTranslatePart(float translateX, float translateY, float translateZ) {
         elems[0][3] = translateX;
         elems[1][3] = translateY;
@@ -676,7 +672,7 @@ namespace J3ML::LinearAlgebra {
 
     Vector3 Matrix4x4::TransformDir(float tx, float ty, float tz) const
     {
-        assert(!this->ContainsProjection()); // This function does not divide by w or output it, so cannot have projection.
+        //assert(!this->ContainsProjection()); // This function does not divide by w or output it, so cannot have projection.
         return Vector3(At(0, 0) * tx + At(0, 1) * ty + At(0, 2) * tz,
                       At(1, 0) * tx + At(1, 1) * ty + At(1, 2) * tz,
                       At(2, 0) * tx + At(2, 1) * ty + At(2, 2) * tz);
@@ -685,7 +681,7 @@ namespace J3ML::LinearAlgebra {
 
     void Matrix4x4::InverseOrthonormal()
     {
-        assert(!ContainsProjection());
+        //assert(!ContainsProjection());
 
         // a) Transposed the top-left 3x3 part in-place to produce R^t.
         Swap(elems[0][1], elems[1][0]);
@@ -777,19 +773,6 @@ namespace J3ML::LinearAlgebra {
         };
     }
 
-    EulerAngle Matrix4x4::ToEulerAngle() const {
-        auto heading = std::atan2(-At(2, 0), At(0, 0));
-        auto attitude = std::asin(At(1, 0));
-        auto bank = std::atan2(-At(1,2), At(1,1));
-        if (At(1, 0) == 1 || At(1, 0) == -1) // North Pole || South Pole
-        {
-            heading = std::atan2(At(0, 2), At(2,2));
-            bank = 0;
-        }
-
-        return {attitude, heading, bank};
-    }
-
     bool Matrix4x4::InverseOrthogonalUniformScale() {
         assert(!ContainsProjection());
         assert(IsColOrthogonal(1e-3f));

src/J3ML/LinearAlgebra/Quaternion.cpp

@@ -1,23 +1,46 @@
-#include <J3ML/LinearAlgebra/Vector3.hpp>
-#include <J3ML/LinearAlgebra/Vector4.hpp>
 #include <J3ML/LinearAlgebra/Matrix3x3.hpp>
 #include <J3ML/LinearAlgebra/Matrix4x4.hpp>
 #include <J3ML/LinearAlgebra/Quaternion.hpp>
 #include <J3ML/LinearAlgebra/AxisAngle.hpp>
+#include <J3ML/LinearAlgebra/EulerAngle.hpp>
+#include <J3ML/LinearAlgebra/Vector4.hpp>
+#include <J3ML/LinearAlgebra/Vector3.hpp>
 
 namespace J3ML::LinearAlgebra {
 
     const Quaternion Quaternion::Identity = Quaternion(0.f, 0.f, 0.f, 1.f);
     const Quaternion Quaternion::NaN      = Quaternion(NAN, NAN, NAN, NAN);
 
-    Quaternion Quaternion::operator-() const
-    {
+    Quaternion Quaternion::operator-() const {
         return {-x, -y, -z, -w};
     }
 
-    Quaternion::Quaternion(const Matrix3x3 &rotationMtrx) {}
+    Quaternion::Quaternion(const Matrix3x3& ro_mat) {
+        auto m = ro_mat.Transposed();
+        auto m00 = m.At(0,0);
+        auto m01 = m.At(0, 1);
+        auto m02 = m.At(0, 2);
+        auto m10 = m.At(1,0);
+        auto m11 = m.At(1, 1);
+        auto m12 = m.At(1, 2);
+        auto m20 = m.At(2,0);
+        auto m21 = m.At(2, 1);
+        auto m22 = m.At(2, 2);
 
-    Quaternion::Quaternion(const Matrix4x4 &rotationMtrx) {}
+        auto field_w = std::sqrt(1.f + m00 + m11 + m22) / 2.f;
+        float w4 = (4.f * field_w);
+
+        x = (m21 - m12) / w4;
+        y = (m02 - m20) / w4;
+        z = (m10 - m01) / w4;
+        w = field_w;
+        Normalize();
+    }
+
+    Quaternion::Quaternion(const Matrix4x4& ro_mat) {
+        auto q = Quaternion(ro_mat.GetRotatePart());
+        x = q.x; y = q.y; z = q.z; w = q.w;
+    }
 
     Vector3 Quaternion::WorldX() const { return Transform(1.f, 0.f, 0.f); }
 
@@ -50,22 +73,6 @@ namespace J3ML::LinearAlgebra {
             return (*this * (t - 1.f) + b * t).Normalized();
     }
 
-    void Quaternion::SetFromAxisAngle(const Vector3 &axis, float angle) {
-        float sinz, cosz;
-        sinz = std::sin(angle*0.5f);
-        cosz = std::cos(angle*0.5f);
-
-        x = axis.x * sinz;
-        y = axis.y * sinz;
-        z = axis.z * sinz;
-        w = cosz;
-    }
-
-    void Quaternion::SetFromAxisAngle(const Vector4 &axis, float angle)
-    {
-        SetFromAxisAngle(Vector3(axis.x, axis.y, axis.z), angle);
-    }
-
     Quaternion Quaternion::operator*(float scalar) const {
         return Quaternion(x * scalar, y * scalar, z * scalar, w * scalar);
     }
@@ -82,8 +89,6 @@ namespace J3ML::LinearAlgebra {
 
     Quaternion Quaternion::operator+() const { return *this; }
 
-    Quaternion::Quaternion() {}
-
     Quaternion::Quaternion(float X, float Y, float Z, float W) : x(X), y(Y), z(Z), w(W) {}
 
     // TODO: implement
@@ -148,20 +153,6 @@ namespace J3ML::LinearAlgebra {
         return (*this * (a * sign) + q2 * b).Normalized();
     }
 
-    AxisAngle Quaternion::ToAxisAngle() const {
-        float halfAngle = std::acos(w);
-        float angle = halfAngle * 2.f;
-        // TODO: Can Implement Fast Inverted Sqrt Here
-        float reciprocalSinAngle = 1.f / std::sqrt(1.f - w*w);
-
-        Vector3 axis = {
-                x*reciprocalSinAngle,
-                y*reciprocalSinAngle,
-                z*reciprocalSinAngle
-        };
-        return AxisAngle(axis, angle);
-    }
-
     float Quaternion::AngleBetween(const Quaternion &target) const {
         Quaternion delta = target / *this;
         return delta.Normalized().Angle();
@@ -212,47 +203,13 @@ namespace J3ML::LinearAlgebra {
         return Vector3(x, y, z) * rcpSinAngle;
     }
 
-    Quaternion::Quaternion(const Vector3 &rotationAxis, float rotationAngleRadians) {
-        SetFromAxisAngle(rotationAxis, rotationAngleRadians);
-    }
-
-    Quaternion::Quaternion(const Vector4 &rotationAxis, float rotationAngleRadians) {
-        SetFromAxisAngle(rotationAxis, rotationAngleRadians);
-    }
-
-    Quaternion::Quaternion(const AxisAngle &angle) {
+    Quaternion::Quaternion(const AxisAngle& angle) {
         double s = std::sin(angle.angle / 2);
         x = angle.axis.x * s;
         y = angle.axis.y * s;
         z = angle.axis.z * s;
         w = std::cos(angle.angle / 2);
-    }
-
-    Quaternion::Quaternion(const EulerAngle &angle) {
-        // Abbreviations for the various angular functions
-        double cr = std::cos(angle.roll * 0.5);
-        double sr = std::sin(angle.roll * 0.5);
-        double cp = std::cos(angle.pitch * 0.5);
-        double sp = std::sin(angle.pitch * 0.5);
-        double cy = std::cos(angle.yaw * 0.5);
-        double sy = std::sin(angle.yaw * 0.5);
-
-        w = cr * cp * cy + sr * sp * sy;
-        x = sr * cp * cy - cr * sp * sy;
-        y = cr * sp * cy + sr * cp * sy;
-        z = cr * cp * sy - sr * sp * cy;
-    }
-
-    void Quaternion::SetFrom(const AxisAngle &angle) {
-        double s = std::sin(angle.angle / 2);
-        x = angle.axis.x * s;
-        y = angle.axis.y * s;
-        z = angle.axis.z * s;
-        w = std::cos(angle.angle / 2);
-    }
-
-    EulerAngle Quaternion::ToEulerAngle() const {
-        return EulerAngle(*this);
+        Normalize();
     }
 
     Quaternion Quaternion::RandomRotation(RNG &rng) {
@@ -271,16 +228,16 @@ namespace J3ML::LinearAlgebra {
         return Quaternion::Identity;
     }
 
-    float Quaternion::Normalize() {
+    bool Quaternion::Normalize() {
         float length = Length();
         if (length < 1e-4f)
-            return 0.f;
+            return false;
         float rcpLength = 1.f / length;
         x *= rcpLength;
         y *= rcpLength;
         z *= rcpLength;
         w *= rcpLength;
-        return length;
+        return true;
     }
 
     bool Quaternion::IsNormalized(float epsilon) const {
@@ -325,9 +282,10 @@ namespace J3ML::LinearAlgebra {
 
     Quaternion Quaternion::LookAt(const Vector3 &localForward, const Vector3 &targetDirection, const Vector3 &localUp,
                                   const Vector3 &worldUp) {
-        return Matrix3x3::LookAt(localForward, targetDirection, localUp, worldUp).ToQuat();
+        return Quaternion(Matrix3x3::LookAt(localForward, targetDirection, localUp, worldUp));
     }
 
+    /*
     Quaternion Quaternion::RotateX(float angleRadians) {
         return {{1,0,0}, angleRadians};
     }
@@ -343,6 +301,7 @@ namespace J3ML::LinearAlgebra {
     Quaternion Quaternion::RotateAxisAngle(const AxisAngle &axisAngle) {
         return {axisAngle.axis, axisAngle.angle};
     }
+    */
 
     Quaternion Quaternion::RotateFromTo(const Vector3 &sourceDirection, const Vector3 &targetDirection) {
         assert(sourceDirection.IsNormalized());
@@ -369,14 +328,6 @@ namespace J3ML::LinearAlgebra {
         return Quaternion::RotateFromTo(sourceDirection.XYZ(), targetDirection.XYZ());
     }
 
-    Quaternion::Quaternion(const float *data) {
-        assert(data);
-        x = data[0];
-        y = data[1];
-        z = data[2];
-        w = data[3];
-    }
-
     Quaternion Quaternion::Lerp(const Quaternion &source, const Quaternion &target, float t) { return source.Lerp(target, t);}
 
     Quaternion Quaternion::Slerp(const Quaternion &source, const Quaternion &target, float t) { return source.Slerp(target, t);}
@@ -401,4 +352,28 @@ namespace J3ML::LinearAlgebra {
     float Quaternion::LengthSquared() const { return x*x + y*y + z*z + w*w;}
 
     float Quaternion::Length() const { return std::sqrt(LengthSquared()); }
+
+    Quaternion::Quaternion(const EulerAngleXYZ& rhs) {
+        float cos_roll = Math::Cos(0.5f * Math::Radians(rhs.roll));
+        float sin_roll = Math::Sin(0.5f * Math::Radians(rhs.roll));
+
+        float cos_pitch = Math::Cos(0.5f * Math::Radians(rhs.pitch));
+        float sin_pitch = Math::Sin(0.5f * Math::Radians(rhs.pitch));
+
+        float cos_yaw = Math::Cos(0.5f * Math::Radians(rhs.yaw));
+        float sin_yaw = Math::Sin(0.5f * Math::Radians(rhs.yaw));
+
+        x = cos_roll * sin_pitch * sin_yaw + sin_roll * cos_pitch * cos_yaw;
+        y = -sin_roll * cos_pitch * sin_yaw + cos_roll * sin_pitch * cos_yaw;
+        z = cos_roll * cos_pitch * sin_yaw + sin_roll * sin_pitch * cos_yaw;
+        w = -sin_roll * sin_pitch * sin_yaw + cos_roll * cos_pitch * cos_yaw;
+        Normalize();
+    }
+
+    Quaternion::Quaternion(const Quaternion& rhs) {
+        x = rhs.x;
+        y = rhs.y;
+        z = rhs.z;
+        w = rhs.w;
+    }
 }

src/J3ML/LinearAlgebra/Vector2.cpp

@@ -3,6 +3,10 @@
 #include <algorithm>
 #include <valarray>
 #include <iostream>
+#include <format>
+#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
+#include <J3ML/LinearAlgebra/Matrix4x4.hpp>
+#include <J3ML/LinearAlgebra/Vector4.hpp>
 
 namespace J3ML::LinearAlgebra {
 
@@ -491,7 +495,27 @@ namespace J3ML::LinearAlgebra {
         }
     }
 
+    std::string Vector2::ToString() const {
+        return std::format("{},{}", x, y);
+    }
+
     Vector2 operator*(float lhs, const Vector2 &rhs) {
         return {lhs * rhs.x, lhs * rhs.y};
     }
+
+    Vector2 Mul2D(const Matrix3x3 &transform, const Vector2 &v) { return transform.Transform(v);}
+
+    std::ostream &operator<<(std::ostream &out, const Vector2 &rhs) {
+        std::string str = rhs.ToString();
+        out << str;
+        return out;
+    }
+
+    Vector2 MulPos2D(const Matrix4x4 &transform, const Vector2 &v) {
+        return transform.Transform(Vector4(v.x, v.y, 0.f, 1.f)).XY();
+    }
+
+    Vector2 MulDir2D(const Matrix4x4 &transform, const Vector2 &v) {
+        return transform.Transform(Vector4(v.x, v.y, 0.f, 0.f)).XY();
+    }
 }

src/J3ML/LinearAlgebra/Vector3.cpp

@@ -2,11 +2,15 @@
 #include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <format>
+#include <ostream>
+#include <J3ML/Geometry/AABB.hpp>
 #include <J3ML/Geometry/Sphere.hpp>
 
 namespace J3ML::LinearAlgebra {
 
     const Vector3 Vector3::Zero     = {0,0,0};
+    const Vector3 Vector3::One      = {1, 1, 1};
     const Vector3 Vector3::Up       = {0, -1, 0};
     const Vector3 Vector3::Down     = {0, 1, 0};
     const Vector3 Vector3::Left     = {-1, 0, 0};
@@ -105,7 +109,7 @@ namespace J3ML::LinearAlgebra {
 
     bool Vector3::operator!=(const Vector3& rhs) const
     {
-        return this->IsWithinMarginOfError(rhs) == false;
+        return !this->IsWithinMarginOfError(rhs);
     }
 
 
@@ -316,11 +320,11 @@ namespace J3ML::LinearAlgebra {
     }
 
     Angle2D Vector3::AngleBetween(const Vector3 &rhs) const {
-        const auto Pi_x_180 = 180.f / M_PI;
+        const auto Pi_x_180 = 180.f / Math::Pi;
         auto dist = this->Distance(rhs);
-        float x = -(asinf((rhs.y - this->y) / dist));
-        float y = (atan2f(rhs.x - this->x,rhs.z - this->z));
-        return {x, y};
+        float dx = -(asinf((rhs.y - this->y) / dist));
+        float dy = (atan2f(rhs.x - this->x,rhs.z - this->z));
+        return {dx, dy};
     }
 
     Angle2D Vector3::AngleBetween(const Vector3 &lhs, const Vector3 &rhs) // TODO: 3D Angle representation?
@@ -329,9 +333,9 @@ namespace J3ML::LinearAlgebra {
     }
 
     Vector3 Vector3::Direction(const Vector3 &rhs) {
-        float x = (cos(Math::Radians(rhs.y)) * cos(Math::Radians(rhs.x)));
-        float y = -sin(Math::Radians(rhs.x));
-        float z = (sin(Math::Radians(rhs.y)) * cos(Math::Radians(rhs.x)));
+        float x = (Math::Cos(Math::Radians(rhs.y)) * Math::Cos(Math::Radians(rhs.x)));
+        float y = -Math::Sin(Math::Radians(rhs.x));
+        float z = (Math::Sin(Math::Radians(rhs.y)) * Math::Cos(Math::Radians(rhs.x)));
         return {x, y, z};
     }
 
@@ -532,7 +536,7 @@ namespace J3ML::LinearAlgebra {
         return std::format("{},{},{}", x, y, z);
     }
 
-    std::ostream & operator<<(std::ostream &o, const Vector3 &vector) {
+    std::ostream& operator<<(std::ostream &o, const Vector3 &vector) {
         o << vector.ToString();
         return o;
     }
@@ -600,5 +604,22 @@ namespace J3ML::LinearAlgebra {
         return *this + Vector3(s, s, s);
     }
 
+    Vector3 Vector3::ScaledToLength(float newLength) const {
+        assert(!IsZero());
+
+        Vector3 v = *this;
+        v.ScaleToLength(newLength);
+        return v;
+    }
+
+    void Vector3::PerpendicularBasis(Vector3 &outB, Vector3 &outC) const {
+        // Pixar orthonormal basis code: https://graphics.pixar.com/library/OrthonormalB/paper.pdf
+        float sign = copysignf(1.0f, z);
+        const float a = -1.0f / (sign + z);
+        const float b = x * y * a;
+        outB = Vector3(1.0f + sign * x * x * a, sign * b,             -sign * x);
+        outC = Vector3(                      b, sign + y * y * a,            -y);
+    }
+
 
 }

src/J3ML/LinearAlgebra/Vector4.cpp

@@ -132,8 +132,8 @@ float Vector4::Magnitude() const
 
 float Vector4::Dot(const Vector4& rhs) const
 {
-    auto a = this->Normalize();
-    auto b = rhs.Normalize();
+    auto a = this->Normalized();
+    auto b = rhs.Normalized();
 
     return  a.x * b.x +
             a.y * b.y +
@@ -155,7 +155,17 @@ Vector4 Vector4::Project(const Vector4& rhs) const
     return rhs * scalar;
 }
 
-Vector4 Vector4::Normalize() const
+void Vector4::Normalize()
+{
+    if (Length() > 0) {
+        x = x / Length();
+        y = y / Length();
+        z = z / Length();
+        w = w / Length();
+    }
+}
+
+Vector4 Vector4::Normalized() const
 {
     if (Length() > 0)
         return {
@@ -274,6 +284,7 @@ Vector4 Vector4::operator-(const Vector4& rhs) const
         if (index == 2) return z;
         if (index == 3) return w;
 
+        return x; // This point should never be reached.
     }
 
     float Vector4::At(int index) const {

tests/Algorithm/RNGTests.hpp

@@ -1,132 +1,149 @@
 #include <J3ML/Algorithm/RNG.hpp>
 #include <jtest/jtest.hpp>
-using J3ML::Algorithm::RNG;
+#include <jtest/Unit.hpp>
 
-int RNGTests()
+jtest::Unit RNGUnit{"RNG"};
+
+
+namespace RNGTests
 {
-    TEST("RNG::IntFast", []{
-        RNG rng;
-        u32 prev = rng.IntFast();
+    inline void Define()
+    {
+        using namespace jtest;
+        using J3ML::Algorithm::RNG;
 
-        for (int i = 0; i < 1000; ++i)
-        {
-            u32 next = rng.IntFast();
-            jtest::check(next != prev);
-            prev = next;
-        }
-    });
+        RNGUnit += Test("IntFast", []{
+            RNG rng;
+            u32 prev = rng.IntFast();
 
-    TEST("RNG::Int", []{
-        RNG rng;
-        assert(rng.lastNumber != 0 || rng.increment != 0);
-        bool allEqual = true;
-        for (int i = 0; i < 1000; ++i)
-        {
-            int prev = rng.Int();
-            int next = rng.Int();
-            jtest::check(prev != 0 || next != 0);
-            if (prev != next)
-                allEqual = false;
-        }
-        jtest::check(!allEqual);
+            for (int i = 0; i < 1000; ++i)
+            {
+                u32 next = rng.IntFast();
+                jtest::check(next != prev);
+                prev = next;
+            }
+        });
 
-    });
+        RNGUnit += Test("Int", []{
+            RNG rng;
+            assert(rng.lastNumber != 0 || rng.increment != 0);
+            bool allEqual = true;
+            for (int i = 0; i < 1000; ++i)
+            {
+                int prev = rng.Int();
+                int next = rng.Int();
+                jtest::check(prev != 0 || next != 0);
+                if (prev != next)
+                    allEqual = false;
+            }
+            jtest::check(!allEqual);
 
-    TEST("Rng::Int_A_B", []{
-        RNG rng;
-        for (int i = 0; i < 1000; ++i)
-        {
-            int a = rng.Int();
-            int b = rng.Int();
-            if (b < a)
-                Swap(a, b);
-            int val = rng.Int(a, b);
-            jtest::check( a <= val);
-            jtest::check(val <= b);
-        }
-    });
+        });
 
-    TEST("Rng::Float", []{
-        RNG rng;
-        bool allEqual = true;
-        for (int i = 0; i < 1000; ++i)
-        {
-            float f = rng.Float();
-            float f2 = rng.Float();
-            jtest::check(f < 1.f);
-            jtest::check(f >= 0.f);
-            jtest::check(f != 0.f || f2 != 0.f);
+        RNGUnit += Test("Int_A_B", []{
+            RNG rng;
+            for (int i = 0; i < 1000; ++i)
+            {
+                int a = rng.Int();
+                int b = rng.Int();
+                if (b < a)
+                    Swap(a, b);
+                int val = rng.Int(a, b);
+                jtest::check( a <= val);
+                jtest::check(val <= b);
+            }
+        });
 
-            if (f != f2)
-                allEqual = false;
-        }
-        jtest::check(!allEqual);
-    });
+        RNGUnit += Test("Float", []{
+            RNG rng;
+            bool allEqual = true;
+            for (int i = 0; i < 1000; ++i)
+            {
+                float f = rng.Float();
+                float f2 = rng.Float();
+                jtest::check(f < 1.f);
+                jtest::check(f >= 0.f);
+                jtest::check(f != 0.f || f2 != 0.f);
 
-    TEST("Rng::Float01Incl", []{
-        RNG rng;
-        bool allEqual = true;
-        for (int i = 0; i < 1000; ++i)
-        {
-            float f = rng.Float01Incl();
-            float f2 = rng.Float01Incl();
+                if (f != f2)
+                    allEqual = false;
+            }
+            jtest::check(!allEqual);
+        });
 
-            jtest::check(f <= 1.f);
-            jtest::check(f >= 0.f);
-            jtest::check(f != 0.f || f2 != 0.f);
-            if (f != f2)
-                allEqual = false;
-        }
-        jtest::check(!allEqual);
-    });
 
-    TEST("Rng::FloatNeg1_1", []{
-        RNG rng;
-        bool allEqual = true;
-        for (int i = 0; i < 1000; ++i)
-        {
-            float f = rng.FloatNeg1_1();
-            float f2 = rng.FloatNeg1_1();
-            jtest::check(f < 1.f);
-            jtest::check(f > -1.f);
-            jtest::check(f != 0.f || f2 != 0.f);
-            if (f != f2)
-                allEqual = false;
-        }
-        jtest::check(!allEqual);
-    });
 
-    TEST("Rng, Float_A_B", []{
-        RNG rng;
-        for (int i = 0; i < 1000; ++i)
-        {
-            float a = rng.Float();
-            float b = rng.Float();
-            if (a == b)
-                continue;
-            if (b < a)
-                Swap(a, b);
+        RNGUnit += Test("Float01Incl", [] {
+            RNG rng;
+            bool allEqual = true;
+            for (int i = 0; i < 1000; ++i)
+            {
+                float f = rng.Float01Incl();
+                float f2 = rng.Float01Incl();
 
-            float f = rng.Float(a, b);
-            jtest::check(a <= f);
-            jtest::check(f < b);
-        }
-    });
+                jtest::check(f <= 1.f);
+                jtest::check(f >= 0.f);
+                jtest::check(f != 0.f || f2 != 0.f);
+                if (f != f2)
+                    allEqual = false;
+            }
+            jtest::check(!allEqual);
+        });
 
-    TEST("Rng::Float_A_B_Incl", []{
-        RNG rng;
-        for (int i = 0; i < 1000; ++i)
-        {
-            float a = rng.Float();
-            float b = rng.Float();
-            if (b > a)
-                Swap(a, b);
+        RNGUnit += Test("FloatNeg1_1", []{
+            RNG rng;
+            bool allEqual = true;
+            for (int i = 0; i < 1000; ++i)
+            {
+                float f = rng.FloatNeg1_1();
+                float f2 = rng.FloatNeg1_1();
+                jtest::check(f < 1.f);
+                jtest::check(f > -1.f);
+                jtest::check(f != 0.f || f2 != 0.f);
+                if (f != f2)
+                    allEqual = false;
+            }
+            jtest::check(!allEqual);
+        });
 
-            float f = rng.FloatIncl(a, b);
-            jtest::check(a <= f);
-            jtest::check(f <= b);
-        }
-    });
+        RNGUnit += Test("Float_A_B", []{
+            RNG rng;
+            for (int i = 0; i < 1000; ++i)
+            {
+                float a = rng.Float();
+                float b = rng.Float();
+                if (a == b)
+                    continue;
+                if (b < a)
+                    Swap(a, b);
 
-    return 0;
+                float f = rng.Float(a, b);
+                jtest::check(a <= f);
+                jtest::check(f < b);
+            }
+        });
+
+
+
+        RNGUnit += Test("Float_A_B_Incl", [] {
+            RNG rng;
+            for (int i = 0; i < 1000; ++i)
+            {
+                float a = rng.Float();
+                float b = rng.Float();
+                if (b < a)
+                    Swap(a, b);
+
+                float f = rng.FloatIncl(a, b);
+                jtest::check(a <= f);
+                jtest::check(f <= b);
+            }
+        });
+
+    }
+
+    inline void Run()
+    {
+        RNGUnit.RunAll();
+    }
 }

tests/Geometry/AABBTests.hpp

@@ -11,81 +11,88 @@
 
 #pragma once
 
-
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
 #include <J3ML/Geometry/AABB.hpp>
+#include <J3ML/Geometry/LineSegment.hpp>
+#include <J3ML/Geometry/Sphere.hpp>
+
+jtest::Unit AABBUnit {"AABB"};
+namespace AABBTests {
+    inline void Define()
+    {
+        using namespace jtest;
+        using namespace J3ML::Geometry;
+
+        AABBUnit += Test("Contains", [] {
+            AABB a ({0,0,0}, {10,10,10});
 
 
-void AABBTests()
-{
+            check(a.Contains({0,0,0}));
+            check(a.Contains({10,10,10}));
+            check(a.Contains({1,2,3}));
+            check(!a.Contains({-1, 2, 3}));
+            check(!a.Contains({1, -2, 3}));
+            check(!a.Contains({1, 2, -3}));
+            check(!a.Contains({11, 2, 3}));
+        });
 
-    using namespace jtest;
-    using namespace J3ML::Geometry;
-
-    TEST("AABB::Contains", [] {
-        AABB a ({0,0,0}, {10,10,10});
+        AABBUnit += Test("ContainsAABB", [] {
+            AABB a ({0,0,0}, {10,10,10});
+            check(a.Contains(a));
+            check(a.Contains(AABB({5, 5, 5}, {6, 6, 6})));
+            check(a.Contains(AABB({5, 5, 5}, {10, 6, 6})));
+            check(!a.Contains(AABB({5,5,5}, {15, 15, 15})));
+            check(!a.Contains(AABB({5,5,5}, {5, 15, 5})));
+            check(!a.Contains(AABB({-5,-5,-5}, {5, 5, 5})));
+            check(!a.Contains(AABB({-5,-5,-5}, {0, 0, 0})));
+        });
 
 
-        check(a.Contains({0,0,0}));
-        check(a.Contains({10,10,10}));
-        check(a.Contains({1,2,3}));
-        check(!a.Contains({-1, 2, 3}));
-        check(!a.Contains({1, -2, 3}));
-        check(!a.Contains({1, 2, -3}));
-        check(!a.Contains({11, 2, 3}));
-    });
+        AABBUnit += Test("ContainsLineSegment", [] {
+            AABB a ({0,0,0}, {10,10,10});
 
-    TEST("AABB::ContainsAABB", [] {
-        AABB a ({0,0,0}, {10,10,10});
-        check(a.Contains(a));
-        check(a.Contains(AABB({5, 5, 5}, {6, 6, 6})));
-        check(a.Contains(AABB({5, 5, 5}, {10, 6, 6})));
-        check(!a.Contains(AABB({5,5,5}, {15, 15, 15})));
-        check(!a.Contains(AABB({5,5,5}, {5, 15, 5})));
-        check(!a.Contains(AABB({-5,-5,-5}, {5, 5, 5})));
-        check(!a.Contains(AABB({-5,-5,-5}, {0, 0, 0})));
-    });
+            check(a.Contains(LineSegment({0,0,0}, {10,10,10})));
+            check(a.Contains(LineSegment({10,10,10}, {0,0,0})));
+            check(a.Contains(LineSegment({5,5,5}, {6,6,6})));
+            check(a.Contains(LineSegment({5,5,5}, {10,6,6})));
+            check(!a.Contains(LineSegment({5,5,5}, {15, 15, 15})));
+            check(!a.Contains(LineSegment({5,5,5}, {5, 15,5})));
+            check(!a.Contains(LineSegment({-5, -5, -5}, {5,5,5})));
+            check(!a.Contains(LineSegment({-5, -5, -5}, {0,0,0})));
+            check(!a.Contains(LineSegment({15, 15, 15}, {0,0,0})));
+        });
+
+        AABBUnit += Test("ContainsSphere", [] {
+            AABB a ({0,0,0}, {10,10,10});
+            check(a.Contains(Sphere({0,0,0}, 0.f)));
+            check(a.Contains(Sphere({5,5,5}, 1.f)));
+            check(!a.Contains(Sphere({5,5,5}, 15.f)));
+            check(!a.Contains(Sphere({9,5,5}, 2.f)));
+            check(!a.Contains(Sphere({1,5,5}, 2.f)));
+            check(!a.Contains(Sphere({-10, -10, -10}, 1000.f)));
+        });
+
+        AABBUnit += Test("IntersectsAABB", [] {
+            AABB a({0,0,0}, {10,10,10});
+            AABB b({5,0,0}, {15,10,10});
+            AABB c({-5,-5,-5}, {0,10,0});
+            AABB d({20,20,20}, {30,30,30});
+            AABB e({1,1,1}, {9,9,9});
+            check(a.Intersects(a));
+            check(a.Intersects(b));
+            check(!a.Intersects(c));
+            check(!a.Intersects(d));
+            check(a.Intersects(e));
+        });
 
 
-    TEST("AABB::ContainsLineSegment", [] {
-        AABB a ({0,0,0}, {10,10,10});
+        AABBUnit += Test("TransformAsAABB", []{ /* TODO: Implement Test Stub */ });
+        AABBUnit += Test("IsDegenerate",    []{ /* TODO: Implement Test Stub */ });
+        AABBUnit += Test("Volume",          []{ /* TODO: Implement Test Stub */ });
+    }
 
-        check(a.Contains(LineSegment({0,0,0}, {10,10,10})));
-        check(a.Contains(LineSegment({10,10,10}, {0,0,0})));
-        check(a.Contains(LineSegment({5,5,5}, {6,6,6})));
-        check(a.Contains(LineSegment({5,5,5}, {10,6,6})));
-        check(!a.Contains(LineSegment({5,5,5}, {15, 15, 15})));
-        check(!a.Contains(LineSegment({5,5,5}, {5, 15,5})));
-        check(!a.Contains(LineSegment({-5, -5, -5}, {5,5,5})));
-        check(!a.Contains(LineSegment({-5, -5, -5}, {0,0,0})));
-        check(!a.Contains(LineSegment({15, 15, 15}, {0,0,0})));
-    });
-
-    TEST("AABB::ContainsSphere", [] {
-        AABB a ({0,0,0}, {10,10,10});
-        check(a.Contains(Sphere({0,0,0}, 0.f)));
-        check(a.Contains(Sphere({5,5,5}, 1.f)));
-        check(!a.Contains(Sphere({5,5,5}, 15.f)));
-        check(!a.Contains(Sphere({9,5,5}, 2.f)));
-        check(!a.Contains(Sphere({1,5,5}, 2.f)));
-        check(!a.Contains(Sphere({-10, -10, -10}, 1000.f)));
-    });
-
-    TEST("AABB::IntersectsAABB", [] {
-        AABB a({0,0,0}, {10,10,10});
-        AABB b({5,0,0}, {15,10,10});
-        AABB c({-5,-5,-5}, {0,10,0});
-        AABB d({20,20,20}, {30,30,30});
-        AABB e({1,1,1}, {9,9,9});
-        check(a.Intersects(a));
-        check(a.Intersects(b));
-        check(!a.Intersects(c));
-        check(!a.Intersects(d));
-        check(a.Intersects(e));
-    });
-
-
-    TEST("AABB::TransformAsAABB", []{ /* TODO: Implement Test Stub */ });
-    TEST("AABB::IsDegenerate",    []{ /* TODO: Implement Test Stub */ });
-    TEST("AABB::Volume",          []{ /* TODO: Implement Test Stub */ });
-
-}
+    inline void Run() {
+        AABBUnit.RunAll();
+    }
+}

tests/Geometry/CommonGeometryTests.hpp

@@ -0,0 +1,67 @@
+#pragma once
+
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
+#include <J3ML/Geometry/Forward.hpp>
+#include <jtest/jtest.hpp>
+
+jtest::Unit CommonGeometryUnit {"CommonGeometry"};
+namespace CommonGeometryTests {
+    inline void Define() {
+        using namespace jtest;
+        using J3ML::Geometry::Interval;
+
+        CommonGeometryUnit += Test("Interval_Intersect", [] {
+            // <- a ->
+            //             <- b ->
+            jtest::check(!(Interval{0, 1}.Intersects({2, 3})));
+
+            //             <- a ->
+            // <- b ->
+            jtest::check(!(Interval{2, 3}.Intersects({0, 1})));
+
+            // <- a ->
+            //    <- b ->
+            jtest::check(!(Interval{2, 4}.Intersects({3, 5})));
+
+            //    <- a ->
+            // <- b ->
+            jtest::check(!(Interval{2, 4}.Intersects({1, 3})));
+
+            // <- a ->
+            //      <- b ->
+            jtest::check(!(Interval{2, 3}.Intersects({3, 5})));
+
+            //      <- a ->
+            // <- b ->
+            jtest::check(!(Interval{3, 5}.Intersects({2, 3})));
+
+            // <- a ->
+            // <- b     ->
+            jtest::check(!(Interval{2, 3}.Intersects({2, 5})));
+
+            // <- a ->
+            // <- b ->
+            jtest::check(!(Interval{2, 3}.Intersects({2, 3})));
+
+            // . a
+            // . b
+            jtest::check(!(Interval{2, 2}.Intersects({2, 2})));
+
+            // <-    a    ->
+            //    <- b ->
+            jtest::check(!(Interval{2, 5}.Intersects({3, 4})));
+
+            //    <- a ->
+            // <-    b    ->
+            jtest::check(!(Interval{3, 4}.Intersects({2, 5})));
+        });
+    }
+
+    inline void Run() {
+        CommonGeometryUnit.RunAll();
+    }
+}
+
+
+

tests/Geometry/FrustumTests.hpp

@@ -12,9 +12,11 @@
 #pragma once
 
 #include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
 
 #include <J3ML/Geometry/Frustum.hpp>
 #include <J3ML/Geometry/PBVolume.hpp>
+#include <J3ML/Geometry/Ray.hpp>
 
 
 Frustum GenIdFrustum(FrustumType t, FrustumHandedness h, FrustumProjectiveSpace p) {
@@ -43,285 +45,295 @@ Frustum GenIdFrustum(FrustumType t, FrustumHandedness h, FrustumProjectiveSpace
         case 7 : f = GenIdFrustum(FrustumType::Orthographic, FrustumHandedness::Right, FrustumProjectiveSpace::D3D); break; \
         } \
 
+jtest::Unit FrustumUnit {"Frustum"};
+namespace FrustumTests
+{
+    inline void Define()
+    {
+        using namespace jtest;
+        using namespace J3ML::Geometry;
+        using namespace J3ML::Math;
 
-void FrustumTests() {
-    using namespace jtest;
-    using namespace J3ML::Geometry;
-    using namespace J3ML::Math;
+        RNG rng;
 
-    RNG rng;
-
-    TEST("Frustum::AspectRatio", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            check(EqualAbs(f.AspectRatio(), 1.f));
-        };
-
-    });
-
-    TEST("Frustum::WorldRight", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            if (f.Handedness() == FrustumHandedness::Right) {
-                check(f.WorldRight().Equals({1, 0, 0}));
-             } else { // In the test func, all cameras lock down to -Z, so left-handed cameras need to point their right towards -X then.
-                 check(f.WorldRight().Equals({-1, 0, 0}));
-             }
-        }
-    });
-
-    TEST("Frustum::Chirality", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            check(f.WorldMatrix().Determinant() > 0.f);
-            check(f.ViewMatrix().Determinant() > 0.f);
-            if (f.Handedness() == FrustumHandedness::Left)
-                check(f.ProjectionMatrix().Determinant4() > 0.f); // left-handed view -> projection space transform does not change handedness.
-            else
-                check(f.ProjectionMatrix().Determinant4() < 0.f); // but right-handed transform should.
-        }
-    });
-
-    TEST("Frustum::Planes", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            check(f.NearPlane().Normal.Equals({0,0,1}));
-            check(Math::EqualAbs(f.NearPlane().distance, -1.f));
-
-            check(f.FarPlane().Normal.Equals({0,0,-1}));
-            check(Math::EqualAbs(f.FarPlane().distance, 100.f));
-
-            for (int i = 0; i < 9; ++i) {
-                Vector3 pt;
-                if (i == 8)
-                    pt = f.CenterPoint();
-                else
-                    pt = f.CornerPoint(i);
-                check(f.NearPlane().SignedDistance(pt) < 1e-3f);
-                check(f.FarPlane().SignedDistance(pt) < 1e-3f);
-                check(f.LeftPlane().SignedDistance(pt) < 1e-3f);
-                check(f.RightPlane().SignedDistance(pt) < 1e-3f);
-                check(f.TopPlane().SignedDistance(pt) < 1e-3f);
-                check(f.BottomPlane().SignedDistance(pt) < 1e-3f);
-                check(f.Contains(pt));
+        FrustumUnit += Test("Frustum::AspectRatio", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                check(EqualAbs(f.AspectRatio(), 1.f));
             }
-        }
-    });
 
+        });
 
-    TEST("Frustum::Corners", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-
-            // Corner points are returned in XYZ order: 0: ---, 1: --+, 2: -+-, 3: -++, 4: +--, 5: +-+, 6: ++-, 7: +++
-            if (f.Type() == FrustumType::Perspective && f.Handedness() == FrustumHandedness::Left) {
-                check(f.CornerPoint(0).Equals({1.f, -1.f, -1.f}));
-                check(f.CornerPoint(1).Equals({100.f, -100.f, -100.f}));
-                check(f.CornerPoint(2).Equals({1.f, 1, -1.f}));
-                check(f.CornerPoint(3).Equals({100.f, 100.f, 100.f}));
-                check(f.CornerPoint(4).Equals({-1.f, -1.f, -1.f}));
-                check(f.CornerPoint(5).Equals({-100.f, -100.f, -100.f}));
-                check(f.CornerPoint(6).Equals({-1.f, 1.f, -1.f}));
-                check(f.CornerPoint(7).Equals({-100.f, 100.f, -100.f}));
-            } else if (f.Type() == FrustumType::Perspective && f.Handedness() == FrustumHandedness::Right) {
-                check(f.CornerPoint(0).Equals({-1.f, -1.f, -1.f}));
-                check(f.CornerPoint(1).Equals({-100.f, -100.f, -100.f}));
-                check(f.CornerPoint(2).Equals({-1.f, 1.f, -1.f}));
-                check(f.CornerPoint(3).Equals({-100.f, 100.f, -100.f}));
-                check(f.CornerPoint(4).Equals({1.f, -1.f, -1.f}));
-                check(f.CornerPoint(5).Equals({100.f, -100.f, -100.f}));
-                check(f.CornerPoint(6).Equals({1.f, 1.f, -1.f}));
-                check(f.CornerPoint(7).Equals({100.f, 100.f, -100.f}));
-            } else if (f.Type() == FrustumType::Orthographic && f.Handedness() == FrustumHandedness::Left) {
-                check(f.CornerPoint(0).Equals({50.f, -50.f, -1.f}));
-                check(f.CornerPoint(1).Equals({50.f, -50.f, -100.f}));
-                check(f.CornerPoint(2).Equals({50.f, 50.f, -1.f}));
-                check(f.CornerPoint(3).Equals({50.f, 50.f, -100.f}));
-                check(f.CornerPoint(4).Equals({-50.f, -50.f, -1.f}));
-                check(f.CornerPoint(5).Equals({-50.f, -50.f, -100.f}));
-                check(f.CornerPoint(6).Equals({-50.f, 50.f, -1.f}));
-                check(f.CornerPoint(7).Equals({-50.f, 50.f, -100.f}));
-            } else if (f.Type() == FrustumType::Orthographic && f.Handedness() == FrustumHandedness::Right) {
-                check(f.CornerPoint(0).Equals({-50.f, -50.f, -1.f}));
-                check(f.CornerPoint(1).Equals({-50.f, -50.f, -100.f}));
-                check(f.CornerPoint(2).Equals({-50.f, 50.f, -1.f}));
-                check(f.CornerPoint(3).Equals({-50.f, 50.f, -100.f}));
-                check(f.CornerPoint(4).Equals({50.f, -50.f, -1.f}));
-                check(f.CornerPoint(5).Equals({50.f, -50.f, -100.f}));
-                check(f.CornerPoint(6).Equals({50.f, 50.f, -1.f}));
-                check(f.CornerPoint(7).Equals({50.f, 50.f, -100.f}));
-            }
-        }
-    });
-
-    TEST("Frustum::ProjectUnprojectSymmetry", [&rng] {
-
-
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            for (int i = 0; i < 10; ++i) {
-                // Orient and locate the Frustum randomly
-                Matrix3x3 rot = Matrix3x3::RandomRotation(rng);
-                f.Transform(rot);
-                f.SetPos(f.Pos() + Vector3::RandomDir(rng, rng.Float(1.f, 100.f)));
-
-                for (int j = 0; j < 100; ++j) {
-                    Vector2 pt = Vector2::RandomBox(rng, -1.f, 1.f);
-                    Vector3 pos = f.NearPlanePos(pt);
-                    Vector3 pt2 = f.Project(pos);
-                    check(pt.Equals(pt2.XY()));
-
-                    pos = f.FarPlanePos(pt);
-                    pt2 = f.Project(pos);
-                    check(pt.Equals(pt2.XY()));
-
-                    pos = f.PointInside(pt.x, pt.y, rng.Float());
-                    pt2 = f.Project(pos);
+        FrustumUnit += Test("Frustum::WorldRight", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                if (f.Handedness() == FrustumHandedness::Right) {
+                    check(f.WorldRight().Equals({1, 0, 0}));
+                } else { // In the test func, all cameras lock down to -Z, so left-handed cameras need to point their right towards -X then.
+                    check(f.WorldRight().Equals({-1, 0, 0}));
                 }
             }
-        }
-    });
+        });
 
-    TEST("Frustum::PlaneNormalsAreCorrect", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            auto pb = f.ToPBVolume();
-            Vector3 corners[8];
-            Plane planes[6];
-            f.GetCornerPoints(corners);
-            f.GetPlanes(planes);
+        FrustumUnit += Test("Frustum::Chirality", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                check(f.WorldMatrix().Determinant() > 0.f);
+                check(f.ViewMatrix().Determinant() > 0.f);
+                if (f.Handedness() == FrustumHandedness::Left)
+                    check(f.ProjectionMatrix().Determinant4() > 0.f); // left-handed view -> projection space transform does not change handedness.
+                else
+                    check(f.ProjectionMatrix().Determinant4() < 0.f); // but right-handed transform should.
+            }
+        });
 
-            for (int i = 0; i < 8; ++i)
-                check(pb.Contains(corners[i]));
+        FrustumUnit += Test("Frustum::Planes", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                check(f.NearPlane().Normal.Equals({0,0,1}));
+                check(Math::EqualAbs(f.NearPlane().distance, -1.f));
 
-            for(int i = 0; i < 6; ++i)
-                for(int j = 0; j < 8; ++j)
-                    check(planes[i].SignedDistance(corners[j]) <= 0.f);
-        }
-    });
+                check(f.FarPlane().Normal.Equals({0,0,-1}));
+                check(Math::EqualAbs(f.FarPlane().distance, 100.f));
 
-    TEST("Frustum::IsConvex", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            Polyhedron p = f.ToPolyhedron();
+                for (int i = 0; i < 9; ++i) {
+                    Vector3 pt;
+                    if (i == 8)
+                        pt = f.CenterPoint();
+                    else
+                        pt = f.CornerPoint(i);
+                    check(f.NearPlane().SignedDistance(pt) < 1e-3f);
+                    check(f.FarPlane().SignedDistance(pt) < 1e-3f);
+                    check(f.LeftPlane().SignedDistance(pt) < 1e-3f);
+                    check(f.RightPlane().SignedDistance(pt) < 1e-3f);
+                    check(f.TopPlane().SignedDistance(pt) < 1e-3f);
+                    check(f.BottomPlane().SignedDistance(pt) < 1e-3f);
+                    check(f.Contains(pt));
+                }
+            }
+        });
 
-            for(int i = 0; i < 6; ++i)
+
+        FrustumUnit += Test("Frustum::Corners", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+
+                // Corner points are returned in XYZ order: 0: ---, 1: --+, 2: -+-, 3: -++, 4: +--, 5: +-+, 6: ++-, 7: +++
+                if (f.Type() == FrustumType::Perspective && f.Handedness() == FrustumHandedness::Left) {
+                    check(f.CornerPoint(0).Equals({1.f, -1.f, -1.f}));
+                    check(f.CornerPoint(1).Equals({100.f, -100.f, -100.f}));
+                    check(f.CornerPoint(2).Equals({1.f, 1, -1.f}));
+                    check(f.CornerPoint(3).Equals({100.f, 100.f, 100.f}));
+                    check(f.CornerPoint(4).Equals({-1.f, -1.f, -1.f}));
+                    check(f.CornerPoint(5).Equals({-100.f, -100.f, -100.f}));
+                    check(f.CornerPoint(6).Equals({-1.f, 1.f, -1.f}));
+                    check(f.CornerPoint(7).Equals({-100.f, 100.f, -100.f}));
+                } else if (f.Type() == FrustumType::Perspective && f.Handedness() == FrustumHandedness::Right) {
+                    check(f.CornerPoint(0).Equals({-1.f, -1.f, -1.f}));
+                    check(f.CornerPoint(1).Equals({-100.f, -100.f, -100.f}));
+                    check(f.CornerPoint(2).Equals({-1.f, 1.f, -1.f}));
+                    check(f.CornerPoint(3).Equals({-100.f, 100.f, -100.f}));
+                    check(f.CornerPoint(4).Equals({1.f, -1.f, -1.f}));
+                    check(f.CornerPoint(5).Equals({100.f, -100.f, -100.f}));
+                    check(f.CornerPoint(6).Equals({1.f, 1.f, -1.f}));
+                    check(f.CornerPoint(7).Equals({100.f, 100.f, -100.f}));
+                } else if (f.Type() == FrustumType::Orthographic && f.Handedness() == FrustumHandedness::Left) {
+                    check(f.CornerPoint(0).Equals({50.f, -50.f, -1.f}));
+                    check(f.CornerPoint(1).Equals({50.f, -50.f, -100.f}));
+                    check(f.CornerPoint(2).Equals({50.f, 50.f, -1.f}));
+                    check(f.CornerPoint(3).Equals({50.f, 50.f, -100.f}));
+                    check(f.CornerPoint(4).Equals({-50.f, -50.f, -1.f}));
+                    check(f.CornerPoint(5).Equals({-50.f, -50.f, -100.f}));
+                    check(f.CornerPoint(6).Equals({-50.f, 50.f, -1.f}));
+                    check(f.CornerPoint(7).Equals({-50.f, 50.f, -100.f}));
+                } else if (f.Type() == FrustumType::Orthographic && f.Handedness() == FrustumHandedness::Right) {
+                    check(f.CornerPoint(0).Equals({-50.f, -50.f, -1.f}));
+                    check(f.CornerPoint(1).Equals({-50.f, -50.f, -100.f}));
+                    check(f.CornerPoint(2).Equals({-50.f, 50.f, -1.f}));
+                    check(f.CornerPoint(3).Equals({-50.f, 50.f, -100.f}));
+                    check(f.CornerPoint(4).Equals({50.f, -50.f, -1.f}));
+                    check(f.CornerPoint(5).Equals({50.f, -50.f, -100.f}));
+                    check(f.CornerPoint(6).Equals({50.f, 50.f, -1.f}));
+                    check(f.CornerPoint(7).Equals({50.f, 50.f, -100.f}));
+                }
+            }
+        });
+
+        FrustumUnit += Test("Frustum::ProjectUnprojectSymmetry", [&rng] {
+
+
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                for (int i = 0; i < 10; ++i) {
+                    // Orient and locate the Frustum randomly
+                    Matrix3x3 rot = Matrix3x3::RandomRotation(rng);
+                    f.Transform(rot);
+                    f.SetPos(f.Pos() + Vector3::RandomDir(rng, rng.Float(1.f, 100.f)));
+
+                    for (int j = 0; j < 100; ++j) {
+                        Vector2 pt = Vector2::RandomBox(rng, -1.f, 1.f);
+                        Vector3 pos = f.NearPlanePos(pt);
+                        Vector3 pt2 = f.Project(pos);
+                        check(pt.Equals(pt2.XY()));
+
+                        pos = f.FarPlanePos(pt);
+                        pt2 = f.Project(pos);
+                        check(pt.Equals(pt2.XY()));
+
+                        pos = f.PointInside(pt.x, pt.y, rng.Float());
+                        pt2 = f.Project(pos);
+                    }
+                }
+            }
+        });
+
+        FrustumUnit += Test("Frustum::PlaneNormalsAreCorrect", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                auto pb = f.ToPBVolume();
+                Vector3 corners[8];
+                Plane planes[6];
+                f.GetCornerPoints(corners);
+                f.GetPlanes(planes);
+
+                for (int i = 0; i < 8; ++i)
+                    check(pb.Contains(corners[i]));
+
+                for(int i = 0; i < 6; ++i)
+                    for(int j = 0; j < 8; ++j)
+                        check(planes[i].SignedDistance(corners[j]) <= 0.f);
+            }
+        });
+
+        FrustumUnit += Test("Frustum::IsConvex", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                Polyhedron p = f.ToPolyhedron();
+
+                for(int i = 0; i < 6; ++i)
+                {
+                    Plane p1 = f.GetPlane(i);
+                    Plane p2 = p.FacePlane(i);
+                    check(p1.Equals(p2));
+                }
+                check(p.EulerFormulaHolds());
+                check(p.IsClosed());
+                check(p.IsConvex());
+                check(!p.IsNull());
+            }
+        });
+
+        FrustumUnit += Test("Plane::ProjectToNegativeHalf", [] {
+            Plane p(Vector3(0,1,0), 50.f);
+
+            Vector3 neg = Vector3(0,-100.f, 0);
+            Vector3 pos = Vector3(0, 100.f, 0);
+            check(neg.Equals(p.ProjectToNegativeHalf(neg)));
+            check(!neg.Equals(p.ProjectToPositiveHalf(neg)));
+
+            check(pos.Equals(p.ProjectToPositiveHalf(pos)));
+            check(!pos.Equals(p.ProjectToNegativeHalf(pos)));
+        });
+
+        FrustumUnit += Test("Frustum::Contains", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                for(int i = 0; i < 8; ++i)
+                {
+                    Vector3 corner = f.CornerPoint(i);
+                    Vector3 closestPoint = f.ClosestPoint(corner);
+                    float distance = f.Distance(corner);
+                    if (!f.Contains(corner) || distance > 1e-4f)
+                        //LOGE("Closest point to %s: %s", corner.ToString().c_str(), closestPoint.ToString().c_str());
+                        check(f.Contains(corner));
+                    check(distance < 10.f);
+                }
+                check(f.Contains(f.CenterPoint()));
+            }
+        });
+
+        FrustumUnit += Test("Frustum::ContainsCorners", [&rng] {
+            constexpr float SCALE = 1e2f;
+
+            Vector3 pt = Vector3::RandomBox(rng, Vector3::FromScalar(-SCALE), Vector3::FromScalar(SCALE));
+            Frustum b;// = RandomFrustumContainingPoint(rng, pt);
+
+            for(int i = 0; i < 9; ++i)
             {
-                Plane p1 = f.GetPlane(i);
-                Plane p2 = p.FacePlane(i);
-                check(p1.Equals(p2));
+                Vector3 point = (i == 8) ? b.CenterPoint() : b.CornerPoint(i);
+
+                check(b.NearPlane().SignedDistance(point) < 1e-3f);
+                check(b.FarPlane().SignedDistance(point) < 1e-3f);
+                check(b.LeftPlane().SignedDistance(point) < 1e-3f);
+                check(b.RightPlane().SignedDistance(point) < 1e-3f);
+                check(b.TopPlane().SignedDistance(point) < 1e-3f);
+                check(b.BottomPlane().SignedDistance(point) < 1e-3f);
+                check(b.Contains(point));
             }
-            check(p.EulerFormulaHolds());
-            check(p.IsClosed());
-            check(p.IsConvex());
-            check(!p.IsNull());
-        }
-    });
+        });
 
-    TEST("Plane::ProjectToNegativeHalf", [] {
-        Plane p(Vector3(0,1,0), 50.f);
+        FrustumUnit += Test("Frustum::Matrices", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                if (f.Handedness() == FrustumHandedness::Right) {
+                    Matrix4x4 wm = f.WorldMatrix();
+                    check(wm.IsIdentity());
 
-        Vector3 neg = Vector3(0,-100.f, 0);
-        Vector3 pos = Vector3(0, 100.f, 0);
-        check(neg.Equals(p.ProjectToNegativeHalf(neg)));
-        check(!neg.Equals(p.ProjectToPositiveHalf(neg)));
+                    Matrix4x4 vm = f.ViewMatrix();
+                    check(vm.IsIdentity());
+                } else {
+                    Matrix4x4 wm = f.WorldMatrix() * Matrix4x4::RotateY(Math::Pi);
+                    check(wm.IsIdentity());
 
-        check(pos.Equals(p.ProjectToPositiveHalf(pos)));
-        check(!pos.Equals(p.ProjectToNegativeHalf(pos)));
-    });
-
-    TEST("Frustum::Contains", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            for(int i = 0; i < 8; ++i)
-            {
-                Vector3 corner = f.CornerPoint(i);
-                Vector3 closestPoint = f.ClosestPoint(corner);
-                float distance = f.Distance(corner);
-                if (!f.Contains(corner) || distance > 1e-4f)
-                    //LOGE("Closest point to %s: %s", corner.ToString().c_str(), closestPoint.ToString().c_str());
-                check(f.Contains(corner));
-                check(distance < 10.f);
+                    Matrix4x4 vm = f.ViewMatrix() * Matrix4x4::RotateY(Math::Pi);
+                    check(vm.IsIdentity());
+                }
             }
-            check(f.Contains(f.CenterPoint()));
-        }
-    });
+        });
 
-    TEST("Frustum::ContainsCorners", [&rng] {
-        constexpr float SCALE = 1e2f;
+        FrustumUnit += Test("Frustum::Projection", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                const float nearD = (f.ProjectiveSpace() == FrustumProjectiveSpace::D3D) ? 0.f : -1.f;
 
-        Vector3 pt = Vector3::RandomBox(rng, Vector3::FromScalar(-SCALE), Vector3::FromScalar(SCALE));
-        Frustum b;// = RandomFrustumContainingPoint(rng, pt);
-
-        for(int i = 0; i < 9; ++i)
-        {
-            Vector3 point = (i == 8) ? b.CenterPoint() : b.CornerPoint(i);
-
-            check(b.NearPlane().SignedDistance(point) < 1e-3f);
-            check(b.FarPlane().SignedDistance(point) < 1e-3f);
-            check(b.LeftPlane().SignedDistance(point) < 1e-3f);
-            check(b.RightPlane().SignedDistance(point) < 1e-3f);
-            check(b.TopPlane().SignedDistance(point) < 1e-3f);
-            check(b.BottomPlane().SignedDistance(point) < 1e-3f);
-            check(b.Contains(point));
-        }
-    });
-
-    TEST("Frustum::Matrices", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            if (f.Handedness() == FrustumHandedness::Right) {
-                Matrix4x4 wm = f.WorldMatrix();
-                check(wm.IsIdentity());
-
-                Matrix4x4 vm = f.ViewMatrix();
-                check(vm.IsIdentity());
-            } else {
-                Matrix4x4 wm = f.WorldMatrix() * Matrix4x4::RotateY(Math::Pi);
-                check(wm.IsIdentity());
-
-                Matrix4x4 vm = f.ViewMatrix() * Matrix4x4::RotateY(Math::Pi);
-                check(vm.IsIdentity());
+                // Corner points are returned in XYZ order: 0: ---, 1: --+, 2: -+-, 3: -++, 4: +--, 5: +-+, 6: ++-, 7: +++
+                check(f.Project(f.CornerPoint(0)).Equals({-1,-1, nearD}));
+                check(f.Project(f.CornerPoint(1)).Equals({-1,-1, 1}));
+                check(f.Project(f.CornerPoint(2)).Equals({-1,1, nearD}));
+                check(f.Project(f.CornerPoint(3)).Equals({-1,1, 1}));
+                check(f.Project(f.CornerPoint(4)).Equals({1,-1, nearD}));
+                check(f.Project(f.CornerPoint(5)).Equals({1,-1, 1}));
+                check(f.Project(f.CornerPoint(6)).Equals({1,1, nearD}));
+                check(f.Project(f.CornerPoint(7)).Equals({1,1, 1}));
             }
-        }
-    });
+        });
 
-    TEST("Frustum::Projection", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            const float nearD = (f.ProjectiveSpace() == FrustumProjectiveSpace::D3D) ? 0.f : -1.f;
-
-            // Corner points are returned in XYZ order: 0: ---, 1: --+, 2: -+-, 3: -++, 4: +--, 5: +-+, 6: ++-, 7: +++
-            check(f.Project(f.CornerPoint(0)).Equals({-1,-1, nearD}));
-            check(f.Project(f.CornerPoint(1)).Equals({-1,-1, 1}));
-            check(f.Project(f.CornerPoint(2)).Equals({-1,1, nearD}));
-            check(f.Project(f.CornerPoint(3)).Equals({-1,1, 1}));
-            check(f.Project(f.CornerPoint(4)).Equals({1,-1, nearD}));
-            check(f.Project(f.CornerPoint(5)).Equals({1,-1, 1}));
-            check(f.Project(f.CornerPoint(6)).Equals({1,1, nearD}));
-            check(f.Project(f.CornerPoint(7)).Equals({1,1, 1}));
-        }
-    });
-
-    TEST("Frustum::UnProject", [] {
-        Frustum f;
-        FOR_EACH_FRUSTUM_CONVENTION(f)
-            if (f.Type() == FrustumType::Perspective) {
-                Ray r = f.UnProject(0, 0);
-                check(r.Origin.Equals(f.Pos()));
-                check(r.Origin.Equals({0,0,0}));
-                check(r.Direction.Equals({0,0,-1}));
+        FrustumUnit += Test("Frustum::UnProject", [] {
+            Frustum f;
+            FOR_EACH_FRUSTUM_CONVENTION(f)
+                if (f.Type() == FrustumType::Perspective) {
+                    Ray r = f.UnProject(0, 0);
+                    check(r.Origin.Equals(f.Pos()));
+                    check(r.Origin.Equals({0,0,0}));
+                    check(r.Direction.Equals({0,0,-1}));
 
 
-                r = f.UnProject(-1, -1);
-                check(r.Origin.Equals(f.Pos()));
-                check(r.Origin.Equals({0,0,0}));
+                    r = f.UnProject(-1, -1);
+                    check(r.Origin.Equals(f.Pos()));
+                    check(r.Origin.Equals({0,0,0}));
+                }
             }
+
+        });
+
+        FrustumUnit += Test("Frustum::UnProjectFromNearPlane", []{});
+
+        FrustumUnit += Test("Frustum::UnProjectFromNearPlane", []{});
     }
 
-    });
+    inline void Run()
+    {
+        //FrustumUnit.RunAll();
+    }
+}
 
-    TEST("Frustum::UnProjectFromNearPlane", []{});
-
-    TEST("Frustum::UnProjectFromNearPlane", []{});
-}

tests/Geometry/Geometry.hpp

@@ -1,58 +0,0 @@
-#pragma once
-
-#include <J3ML/Geometry/Common.hpp>
-#include <jtest/jtest.hpp>
-using J3ML::Geometry::Interval;
-
-int CommonGeometryTests() {
-
-    TEST("Geometry::Interval_Intersect", [] {
-        // <- a ->
-        //             <- b ->
-        jtest::check(!(Interval{0, 1}.Intersects({2, 3})));
-
-        //             <- a ->
-        // <- b ->
-        jtest::check(!(Interval{2, 3}.Intersects({0, 1})));
-
-        // <- a ->
-        //    <- b ->
-        jtest::check(!(Interval{2, 4}.Intersects({3, 5})));
-
-        //    <- a ->
-        // <- b ->
-        jtest::check(!(Interval{2, 4}.Intersects({1, 3})));
-
-        // <- a ->
-        //      <- b ->
-        jtest::check(!(Interval{2, 3}.Intersects({3, 5})));
-
-        //      <- a ->
-        // <- b ->
-        jtest::check(!(Interval{3, 5}.Intersects({2, 3})));
-
-        // <- a ->
-        // <- b     ->
-        jtest::check(!(Interval{2, 3}.Intersects({2, 5})));
-
-        // <- a ->
-        // <- b ->
-        jtest::check(!(Interval{2, 3}.Intersects({2, 3})));
-
-        // . a
-        // . b
-        jtest::check(!(Interval{2, 2}.Intersects({2, 2})));
-
-        // <-    a    ->
-        //    <- b ->
-        jtest::check(!(Interval{2, 5}.Intersects({3, 4})));
-
-        //    <- a ->
-        // <-    b    ->
-        jtest::check(!(Interval{3, 4}.Intersects({2, 5})));
-    });
-    return 0;
-}
-
-
-

tests/Geometry/OBBTests.hpp

@@ -7,4 +7,5 @@
 
 /// @file OBBTests.hpp
 /// @desc Unit tests for the Oriented Bounding Box class.
-/// @edit 2024-07-23
+/// @edit 2024-07-23
+

tests/Geometry/TriangleTests.hpp

@@ -1,94 +1,103 @@
 #include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
+#include <jtest/jtest.hpp>
 #include <J3ML/Geometry/Triangle.hpp>
 
-using J3ML::Geometry::Interval;
-using J3ML::Geometry::Triangle;
-
-int TriangleTests() {
-    TEST("Triangle::FaceNormal", [] {
-        Triangle t(
-             {-1, -1, -1},
-             {0,   1,  0},
-             {1,  -1,  1}
-        );
-
-        jtest::check(t.FaceNormal() == Vector3{4, 0, -4});
-    });
-
-    TEST("Triangle::IntersectTriangle", []
+jtest::Unit TriangleUnit {"Triangle"};
+namespace TriangleTests {
+    inline void Define()
     {
-        Triangle xyTriangle(
-            {0.0f, 0.0f, 0.0f},
-            {1.0f, 1.0f, 0.0f},
-            {2.0f, 0.0f, 0.0f}
-        );
+        using namespace jtest;
+        using J3ML::Geometry::Interval;
+        using J3ML::Geometry::Triangle;
 
-        // Triangle collides with itself
-        jtest::check(Intersects(xyTriangle, xyTriangle));
-        // Translate 1 towards x -- should collide
-        jtest::check(Intersects(xyTriangle, xyTriangle.Translated(Vector3(1.0f, 0.0f, 0.0f))));
-        // Translate 2 towards x -- should collide exactly on V1
-        jtest::check(Intersects(xyTriangle, xyTriangle.Translated(Vector3(2.0f, 0.0f, 0.0f))));
-        // Translate 2 towards negative x -- should collide exactly on V0
-        jtest::check(Intersects(xyTriangle, xyTriangle.Translated(Vector3(-2.0f, 0.0f, 0.0f))));
-        // Translate 3 towards x -- should not collide
-        jtest::check(!Intersects(xyTriangle, xyTriangle.Translated(Vector3(3.0f, 0.0f, 0.0f))));
-        // Translate 3 towards negative x -- should not collide
-        jtest::check(!Intersects(xyTriangle, xyTriangle.Translated(Vector3(-3.0f, 0.0f, 0.0f))));
-        // Translate 1 towards z -- should not collide
-        jtest::check(!Intersects(xyTriangle, xyTriangle.Translated(Vector3(0.0f, 0.0f, 1.0f))));
-        // Triangle collides with contained smaller triangle
-        jtest::check(!Intersects(xyTriangle, xyTriangle.Scaled(Vector3(0.5f, 0.5f, 0.5f)).Translated(Vector3(0.25f, 0.25f, 0.0f))));
+        TriangleUnit += Test("FaceNormal", [] {
+            Triangle t(
+                    {-1, -1, -1},
+                    {0,   1,  0},
+                    {1,  -1,  1}
+            );
 
-        Triangle zxTriangle (
-            {0.0f, 0.0f, 0.0f},
-            {1.0f, 0.0f, 1.0f},
-            {0.0f, 0.0f, 2.0f}
-        );
+            jtest::check(t.FaceNormal() == Vector3{4, 0, -4});
+        });
 
-        // Should collide exactly on V0
-        jtest::check(Intersects(xyTriangle, zxTriangle));
-        // Should collide across xyTriangle's edge and  zxTriangle's face
-        jtest::check(Intersects(xyTriangle, zxTriangle.Translated(Vector3(0.0f, 0.0f, -1.0))));
-        // Should collide exactly on V1
-        jtest::check(Intersects(xyTriangle, zxTriangle.Translated(Vector3(0.0f, 0.0f, -2.0))));
-        // xyTriangle's face should be poked by zxTriangle's V0
-        jtest::check(Intersects(xyTriangle, zxTriangle.Translated(Vector3(1.0f, 1.0f, 0.0f))));
-        // xyTriangle's face should be cut by zxTriangle
-        jtest::check(Intersects(xyTriangle, zxTriangle.Translated(Vector3(1.0f, 1.0f, -0.5f))));
-        // Should not collide
-        jtest::check(!Intersects(xyTriangle, zxTriangle.Translated(Vector3(1.0f, 1.0f, 1.0f))));
-        // Should not collide
-        jtest::check(!Intersects(xyTriangle, zxTriangle.Translated(Vector3(0.0f, 0.0f, -3.0f))));
+        TriangleUnit += Test("IntersectTriangle", []
+        {
+            Triangle xyTriangle(
+                    {0.0f, 0.0f, 0.0f},
+                    {1.0f, 1.0f, 0.0f},
+                    {2.0f, 0.0f, 0.0f}
+            );
 
-        Triangle yxTriangle(
-            {0.0f, 0.0f, 0.0f},
-            {1.0f, 1.0f, 0.0f},
-            {0.0f, 2.0f, 0.0f}
-        );
+            // Triangle collides with itself
+            jtest::check(Intersects(xyTriangle, xyTriangle));
+            // Translate 1 towards x -- should collide
+            jtest::check(Intersects(xyTriangle, xyTriangle.Translated(Vector3(1.0f, 0.0f, 0.0f))));
+            // Translate 2 towards x -- should collide exactly on V1
+            jtest::check(Intersects(xyTriangle, xyTriangle.Translated(Vector3(2.0f, 0.0f, 0.0f))));
+            // Translate 2 towards negative x -- should collide exactly on V0
+            jtest::check(Intersects(xyTriangle, xyTriangle.Translated(Vector3(-2.0f, 0.0f, 0.0f))));
+            // Translate 3 towards x -- should not collide
+            jtest::check(!Intersects(xyTriangle, xyTriangle.Translated(Vector3(3.0f, 0.0f, 0.0f))));
+            // Translate 3 towards negative x -- should not collide
+            jtest::check(!Intersects(xyTriangle, xyTriangle.Translated(Vector3(-3.0f, 0.0f, 0.0f))));
+            // Translate 1 towards z -- should not collide
+            jtest::check(!Intersects(xyTriangle, xyTriangle.Translated(Vector3(0.0f, 0.0f, 1.0f))));
+            // Triangle collides with contained smaller triangle
+            jtest::check(!Intersects(xyTriangle, xyTriangle.Scaled(Vector3(0.5f, 0.5f, 0.5f)).Translated(Vector3(0.25f, 0.25f, 0.0f))));
 
-        // Should collide on V0-V1 edge
-        jtest::check(Intersects(yxTriangle, yxTriangle));
-        // Should not collide
-        jtest::check(!Intersects(xyTriangle, yxTriangle.Translated(Vector3(0.0f, 1.0f, 0.0f))));
-        // Should not collide
-        jtest::check(!Intersects(yxTriangle, yxTriangle.Translated(Vector3(0.0f, 0.0f, 1.0f))));
+            Triangle zxTriangle (
+                    {0.0f, 0.0f, 0.0f},
+                    {1.0f, 0.0f, 1.0f},
+                    {0.0f, 0.0f, 2.0f}
+            );
 
-        Triangle zyInvertedTriangle(
-            {0.0f, 1.0f, -1.0f},
-            {0.0f, 0.0f, 0.0f},
-            {0.0f, 1.0f, 1.0f}
-        );
-        // Should collide exactly on V1
-        jtest::check(Intersects(xyTriangle, zyInvertedTriangle));
-        // Should not collide
-        jtest::check(!Intersects(xyTriangle, zyInvertedTriangle.Translated(Vector3(0.0f, 1.0f, 0.0f))));
-        // Should not collide
-        jtest::check(!Intersects(xyTriangle, zyInvertedTriangle.Translated(Vector3(0.25f, 0.75f, 0.0f))));
-    });
+            // Should collide exactly on V0
+            jtest::check(Intersects(xyTriangle, zxTriangle));
+            // Should collide across xyTriangle's edge and  zxTriangle's face
+            jtest::check(Intersects(xyTriangle, zxTriangle.Translated(Vector3(0.0f, 0.0f, -1.0))));
+            // Should collide exactly on V1
+            jtest::check(Intersects(xyTriangle, zxTriangle.Translated(Vector3(0.0f, 0.0f, -2.0))));
+            // xyTriangle's face should be poked by zxTriangle's V0
+            jtest::check(Intersects(xyTriangle, zxTriangle.Translated(Vector3(1.0f, 1.0f, 0.0f))));
+            // xyTriangle's face should be cut by zxTriangle
+            jtest::check(Intersects(xyTriangle, zxTriangle.Translated(Vector3(1.0f, 1.0f, -0.5f))));
+            // Should not collide
+            jtest::check(!Intersects(xyTriangle, zxTriangle.Translated(Vector3(1.0f, 1.0f, 1.0f))));
+            // Should not collide
+            jtest::check(!Intersects(xyTriangle, zxTriangle.Translated(Vector3(0.0f, 0.0f, -3.0f))));
 
+            Triangle yxTriangle(
+                    {0.0f, 0.0f, 0.0f},
+                    {1.0f, 1.0f, 0.0f},
+                    {0.0f, 2.0f, 0.0f}
+            );
 
-    return 0;
+            // Should collide on V0-V1 edge
+            jtest::check(Intersects(yxTriangle, yxTriangle));
+            // Should not collide
+            jtest::check(!Intersects(xyTriangle, yxTriangle.Translated(Vector3(0.0f, 1.0f, 0.0f))));
+            // Should not collide
+            jtest::check(!Intersects(yxTriangle, yxTriangle.Translated(Vector3(0.0f, 0.0f, 1.0f))));
+
+            Triangle zyInvertedTriangle(
+                    {0.0f, 1.0f, -1.0f},
+                    {0.0f, 0.0f, 0.0f},
+                    {0.0f, 1.0f, 1.0f}
+            );
+            // Should collide exactly on V1
+            jtest::check(Intersects(xyTriangle, zyInvertedTriangle));
+            // Should not collide
+            jtest::check(!Intersects(xyTriangle, zyInvertedTriangle.Translated(Vector3(0.0f, 1.0f, 0.0f))));
+            // Should not collide
+            jtest::check(!Intersects(xyTriangle, zyInvertedTriangle.Translated(Vector3(0.25f, 0.75f, 0.0f))));
+        });
+    }
+
+    inline void Run()
+    {
+        TriangleUnit.RunAll();
+    }
 }
 
 

tests/LinearAlgebra/AxisAngleTests.hpp

@@ -0,0 +1,25 @@
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
+
+jtest::Unit AxisAngleUnit {"AxisAngle"};
+
+namespace AxisAngleTests {
+    inline void Define() {
+        using namespace jtest;
+
+        AxisAngleUnit += Test("From_Quaternion", [] {
+            AxisAngle expected_result({0.3860166, 0.4380138, 0.8118714}, 0.6742209);
+            Quaternion q(0.1276794, 0.1448781, 0.2685358, 0.9437144);
+
+            AxisAngle from_quaternion(q);
+
+            jtest::check(Math::EqualAbs(expected_result.axis.x, from_quaternion.axis.x, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.axis.y, from_quaternion.axis.y, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.axis.z, from_quaternion.axis.z, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.angle, from_quaternion.angle, 1e-6f));
+        });
+    }
+    inline void Run() {
+        AxisAngleUnit.RunAll();
+    }
+}

tests/LinearAlgebra/EulerAngleTests.hpp

@@ -2,6 +2,27 @@
 // Created by josh on 12/26/2023.
 //
 
-int EulerAngleTests() {
-    return 0;
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
+
+jtest::Unit EulerAngleUnit {"EulerAngle_XYZ"};
+
+namespace EulerAngleTests {
+    inline void Define() {
+        using namespace jtest;
+
+        EulerAngleUnit += Test("From_Quaternion", [] {
+            EulerAngleXYZ expected_result(-170, 88, -160);
+            Quaternion q(0.1840604, 0.6952024, 0.1819093, 0.6706149);
+
+            EulerAngleXYZ from_quaternion(q);
+
+            jtest::check(Math::EqualAbs(Math::Radians(expected_result.roll), Math::Radians(from_quaternion.roll), 1e-5f));
+            jtest::check(Math::EqualAbs(Math::Radians(expected_result.pitch), Math::Radians(from_quaternion.pitch), 1e-5f));
+            jtest::check(Math::EqualAbs(Math::Radians(expected_result.yaw), Math::Radians(from_quaternion.yaw), 1e-5f));
+        });
+    }
+    inline void Run() {
+        EulerAngleUnit.RunAll();
+    }
 }

tests/LinearAlgebra/Matrix2x2Tests.hpp

@@ -1,8 +1,19 @@
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
 
-#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
+#include <J3ML/LinearAlgebra/Matrix2x2.hpp>
 
-using namespace J3ML::LinearAlgebra;
+jtest::Unit Matrix2x2Unit {"Matrix2x2"};
+namespace Matrix2x2Tests {
+    inline void Define() {
+        using namespace jtest;
+        using namespace J3ML::LinearAlgebra;
 
-int Matrix2x2Tests() {
-    return 0;
+        Matrix2x2Unit += Test("Not Implemented", [] {
+            throw("Not Implemented");
+        });
+    }
+    inline void Run() {
+        Matrix2x2Unit.RunAll();
+    }
 }

tests/LinearAlgebra/Matrix3x3Tests.hpp

@@ -1,96 +1,139 @@
-#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
 #include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
+#include <J3ML/LinearAlgebra/Matrix3x3.hpp>
+#include <J3ML/LinearAlgebra/Matrix4x4.hpp>
 
-using namespace J3ML::LinearAlgebra;
-
-// TODO: create RNG instance
-
-int Matrix3x3Tests() {
-
-
-    TEST("Mat3x3::Add_Unary", []
+jtest::Unit Matrix3x3Unit {"Matrix3x3"};
+namespace Matrix3x3Tests
+{
+    inline void Define()
     {
-        Matrix3x3 m(1,2,3, 4,5,6, 7,8,9);
-        Matrix3x3 m2 = +m;
-        jtest::check(m.Equals(m2));
-    });
+        using namespace jtest;
+        using namespace J3ML::LinearAlgebra;
 
-    TEST("Mat3x3::Solve_Axb", []{
-        RNG rng;
-        Matrix3x3 A = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
-        bool mayFail = Math::EqualAbs(A.Determinant(), 0.f, 1e-2f);
+        Matrix3x3Unit += Test("AngleTypeRound-TripConversion", [] {
+            EulerAngleXYZ expected_result(8, 60, -27);
 
-        Vector3 b = Vector3::RandomBox(rng, Vector3::FromScalar(-10.f), Vector3::FromScalar(10.f));
+            Matrix3x3 m(expected_result);
+            AxisAngle a(expected_result);
+            Quaternion q(a);
+            Matrix3x3 m2(q);
+            Quaternion q2(m2);
+            AxisAngle a2(q2);
+            EulerAngleXYZ round_trip(a2);
 
-        Vector3 x;
-        bool success = A.SolveAxb(b, x);
-        jtest::check(success || mayFail);
-        if (success)
+            jtest::check(Math::EqualAbs(Math::Radians(expected_result.roll), Math::Radians(round_trip.roll), 1e-6f));
+            jtest::check(Math::EqualAbs(Math::Radians(expected_result.pitch), Math::Radians(round_trip.pitch), 1e-6f));
+            jtest::check(Math::EqualAbs(Math::Radians(expected_result.yaw), Math::Radians(round_trip.yaw), 1e-6f));
+        });
+
+        Matrix3x3Unit += Test("From_EulerAngleXYZ", []{
+            Matrix3x3 expected_result(Vector3(0.4455033, 0.2269952, 0.8660254),
+                                      Vector3(-0.3421816, 0.9370536,  -0.0695866),
+                                      Vector3(-0.8273081, -0.2653369, 0.4951340)
+                                      );
+
+            EulerAngleXYZ e(8, 60, -27);
+            Matrix3x3 from_euler(e);
+
+            jtest::check(Math::EqualAbs(expected_result.At(0, 0), from_euler.At(0, 0), 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.At(0, 1), from_euler.At(0, 1), 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.At(0, 2), from_euler.At(0, 2), 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.At(1, 0), from_euler.At(1, 0), 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.At(1, 1), from_euler.At(1, 1), 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.At(1, 2), from_euler.At(1, 2), 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.At(2, 0), from_euler.At(2, 0), 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.At(2, 1), from_euler.At(2, 1), 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.At(2, 2), from_euler.At(2, 2), 1e-6f));
+        });
+
+        Matrix3x3Unit += Test("Add_Unary", []
         {
-            Vector3 b2 = A*x;
-            jtest::check(b2.Equals(b, 1e-1f));
-        }
-    });
+            Matrix3x3 m(1,2,3, 4,5,6, 7,8,9);
+            Matrix3x3 m2 = +m;
+            jtest::check(m.Equals(m2));
+        });
 
-    TEST("Mat3x3::Inverse_Case", []
-    {
-        Matrix3x3 m(-8.75243664f,6.71196938f,-5.95816374f,6.81996822f,-6.85106039f,2.38949537f,-0.856015682f,3.45762491f,3.311584f);
+        Matrix3x3Unit += Test("Solve_Axb", []{
+            RNG rng;
+            Matrix3x3 A = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
+            bool mayFail = Math::EqualAbs(A.Determinant(), 0.f, 1e-2f);
 
-        bool success = m.Inverse();
-        jtest::check(success);
-    });
+            Vector3 b = Vector3::RandomBox(rng, Vector3::FromScalar(-10.f), Vector3::FromScalar(10.f));
 
-    TEST("Mat3x3::Inverse", []
-    {
-        RNG rng;
-        Matrix3x3 A = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
-        bool mayFail = Math::EqualAbs(A.Determinant(), 0.f, 1e-2f);
+            Vector3 x;
+            bool success = A.SolveAxb(b, x);
+            jtest::check(success || mayFail);
+            if (success)
+            {
+                Vector3 b2 = A*x;
+                jtest::check(b2.Equals(b, 1e-1f));
+            }
+        });
 
-        Matrix3x3 A2 = A;
-        bool success = A2.Inverse();
-        jtest::check(success || mayFail);
-        if (success)
+        Matrix3x3Unit += Test("Inverse_Case", []
         {
-            Matrix3x3 id = A * A2;
-            Matrix3x3 id2 = A2 * A;
-            jtest::check(id.Equals(Matrix3x3::Identity, 0.3f));
-            jtest::check(id2.Equals(Matrix3x3::Identity, 0.3f));
-        }
-    });
+            Matrix3x3 m(-8.75243664f,6.71196938f,-5.95816374f,6.81996822f,-6.85106039f,2.38949537f,-0.856015682f,3.45762491f,3.311584f);
 
+            bool success = m.Inverse();
+            jtest::check(success);
+        });
 
-    TEST("Mat3x3::InverseFast", []
-    {
-        // TODO: Fix implementation of InverseFast
-        /*
-        Matrix3x3 A = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
-        bool mayFail = Math::EqualAbs(A.Determinant(), 0.f, 1e-2f);
-
-        Matrix3x3 A2 = A;
-        bool success = A2.InverseFast();
-        assert(success || mayFail);
-
-        if (success)
+        Matrix3x3Unit += Test("Inverse", []
         {
-            Matrix3x3 id = A * A2;
-            Matrix3x3 id2 = A2 * A;
-            assert(id.Equals(Matrix3x3::Identity, 0.3f));
-            assert(id2.Equals(Matrix3x3::Identity, 0.3f));
-        }
-         */
-    });
+            RNG rng;
+            Matrix3x3 A = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
+            bool mayFail = Math::EqualAbs(A.Determinant(), 0.f, 1e-2f);
 
-    TEST("Mat3x3::MulMat4x4", []
+            Matrix3x3 A2 = A;
+            bool success = A2.Inverse();
+            jtest::check(success || mayFail);
+            if (success)
+            {
+                Matrix3x3 id = A * A2;
+                Matrix3x3 id2 = A2 * A;
+                jtest::check(id.Equals(Matrix3x3::Identity, 0.3f));
+                jtest::check(id2.Equals(Matrix3x3::Identity, 0.3f));
+            }
+        });
+
+
+        Matrix3x3Unit += Test("InverseFast", []
+        {
+            // TODO: Fix implementation of InverseFast
+            /*
+            Matrix3x3 A = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
+            bool mayFail = Math::EqualAbs(A.Determinant(), 0.f, 1e-2f);
+
+            Matrix3x3 A2 = A;
+            bool success = A2.InverseFast();
+            assert(success || mayFail);
+
+            if (success)
+            {
+                Matrix3x3 id = A * A2;
+                Matrix3x3 id2 = A2 * A;
+                assert(id.Equals(Matrix3x3::Identity, 0.3f));
+                assert(id2.Equals(Matrix3x3::Identity, 0.3f));
+            }
+             */
+        });
+
+        Matrix3x3Unit += Test("MulMat4x4", []
+        {
+            RNG rng;
+            Matrix3x3 m = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
+            Matrix4x4 m_ = m;
+            Matrix4x4 m2 = Matrix4x4::RandomGeneral(rng, -10.f, 10.f);
+
+            Matrix4x4 test = m * m2;
+            Matrix4x4 correct = m_ * m2;
+            jtest::check(test.Equals(correct));
+        });
+    }
+
+    inline void Run()
     {
-        RNG rng;
-        Matrix3x3 m = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
-        Matrix4x4 m_ = m;
-        Matrix4x4 m2 = Matrix4x4::RandomGeneral(rng, -10.f, 10.f);
-
-        Matrix4x4 test = m * m2;
-        Matrix4x4 correct = m_ * m2;
-        jtest::check(test.Equals(correct));
-    });
-
-    return 0;
+        Matrix3x3Unit.RunAll();
+    }
 }

tests/LinearAlgebra/Matrix4x4Tests.hpp

@@ -1,129 +1,121 @@
 #include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
 #include <J3ML/LinearAlgebra/Matrix4x4.hpp>
 #include <J3ML/LinearAlgebra/Quaternion.hpp>
+#include <J3ML/Math.hpp>
 
-using namespace J3ML::LinearAlgebra;
+jtest::Unit Matrix4x4Unit {"Matrix4x4"};
+namespace Matrix4x4Tests {
+    inline void Define() {
+        using namespace jtest;
+        using namespace J3ML::LinearAlgebra;
+        using namespace J3ML::Math;
 
-void Matrix4x4_AngleTypeRoundtripConversion()
-{
-    Matrix4x4 matrix;
-    EulerAngle a(Math::Radians(45), Math::Radians(45), Math::Radians(45));
-    Quaternion q(a);
-    matrix.SetRotatePart(q);
-    //matrix.SetRotatePartX(a.pitch);
-    //matrix.SetRotatePartY(a.yaw);
-    //matrix.SetRotatePartZ(a.roll);
-    EulerAngle fromMatrix = matrix.GetRotatePart().ToQuat().ToEulerAngle();
-    jtest::check(a == fromMatrix);
-}
-
-int Matrix4x4Tests() {
-
-    TEST("Mat4x4::Add_Unary", []
-    {
-        Matrix4x4 m(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16);
-        Matrix4x4 m2 = +m;
-        jtest::check(m.Equals(m2));
-    });
-
-    TEST("Mat4x4::Inverse", [] {
-        RNG rng;
-        Matrix4x4 A = Matrix4x4::RandomGeneral(rng, -10.f, 10.f);
-        bool mayFail = Math::EqualAbs(A.Determinant(), 0.f, 1e-2f);
-
-        Matrix4x4 A2 = A;
-        bool success = A2.Inverse();
-
-        jtest::check(success || mayFail);
-
-        if (success)
+        Matrix4x4Unit += Test("Add_Unary", []
         {
-            Matrix4x4 id = A * A2;
-            Matrix4x4 id2 = A2 * A;
+            Matrix4x4 m(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16);
+            Matrix4x4 m2 = +m;
+            jtest::check(m.Equals(m2));
+        });
 
-            jtest::check(id.Equals(Matrix4x4::Identity, 0.3f));
-            jtest::check(id2.Equals(Matrix4x4::Identity, 0.3f));
-        }
-    });
+        Matrix4x4Unit += Test("Inverse", [] {
+            RNG rng;
+            Matrix4x4 A = Matrix4x4::RandomGeneral(rng, -10.f, 10.f);
+            bool mayFail = Math::EqualAbs(A.Determinant(), 0.f, 1e-2f);
 
-    TEST("Mat4x4::Ctor", []{
-        RNG rng;
-        Matrix3x3 m = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
+            Matrix4x4 A2 = A;
+            bool success = A2.Inverse();
 
-        Matrix4x4 m2(m);
+            jtest::check(success || mayFail);
 
-        for (int y = 0; y < 3; ++y)
-            for (int x = 0; x < 3; ++x)
-                assert(Math::EqualAbs(m.At(y, x), m2.At(y, x)));
+            if (success)
+            {
+                Matrix4x4 id = A * A2;
+                Matrix4x4 id2 = A2 * A;
 
-        jtest::check(Math::EqualAbs(m2[0][3], 0.f));
-        jtest::check(Math::EqualAbs(m2[1][3], 0.f));
-        jtest::check(Math::EqualAbs(m2[2][3], 0.f));
+                jtest::check(id.Equals(Matrix4x4::Identity, 0.3f));
+                jtest::check(id2.Equals(Matrix4x4::Identity, 0.3f));
+            }
+        });
 
-        jtest::check(Math::EqualAbs(m2[3][0], 0.f));
-        jtest::check(Math::EqualAbs(m2[3][1], 0.f));
-        jtest::check(Math::EqualAbs(m2[3][2], 0.f));
-        jtest::check(Math::EqualAbs(m2[3][3], 0.f));
-    });
+        Matrix4x4Unit += Test("Ctor", []{
+            RNG rng;
+            Matrix3x3 m = Matrix3x3::RandomGeneral(rng, -10.f, 10.f);
 
-    TEST("Mat4x4::SetRow", []
-    {
-        Matrix4x4 m;
-        m.SetRow(0, 1,2,3,4);
-        m.SetRow(1, Vector4(5,6,7,8));
-        m.SetRow(2, 9,10,11,12);
-        m.SetRow(3, 13, 14, 15, 16);
+            Matrix4x4 m2(m);
+
+            for (int y = 0; y < 3; ++y)
+                for (int x = 0; x < 3; ++x)
+                    assert(Math::EqualAbs(m.At(y, x), m2.At(y, x)));
+
+            jtest::check(Math::EqualAbs(m2[0][3], 0.f));
+            jtest::check(Math::EqualAbs(m2[1][3], 0.f));
+            jtest::check(Math::EqualAbs(m2[2][3], 0.f));
+
+            jtest::check(Math::EqualAbs(m2[3][0], 0.f));
+            jtest::check(Math::EqualAbs(m2[3][1], 0.f));
+            jtest::check(Math::EqualAbs(m2[3][2], 0.f));
+            jtest::check(Math::EqualAbs(m2[3][3], 0.f));
+        });
+
+        Matrix4x4Unit += Test("SetRow", []
+        {
+            Matrix4x4 m;
+            m.SetRow(0, 1,2,3,4);
+            m.SetRow(1, Vector4(5,6,7,8));
+            m.SetRow(2, 9,10,11,12);
+            m.SetRow(3, 13, 14, 15, 16);
 
 
-        Matrix4x4 m3(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16);
-        Matrix4x4 m2;
-        m2.Set(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16);
-        jtest::check(m.Equals(m2));
-        jtest::check(m.Equals(m3));
-    });
+            Matrix4x4 m3(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16);
+            Matrix4x4 m2;
+            m2.Set(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16);
+            jtest::check(m.Equals(m2));
+            jtest::check(m.Equals(m3));
+        });
 
-    TEST("Mat4x4::SwapRows", []
-    {
-        Matrix4x4 m(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16);
-        Matrix4x4 m2(13,14,15,16, 9,10,11,12, 5,6,7,8, 1,2,3,4);
-        m.SwapRows(0,3);
-        m.SwapRows(1,2);
-        jtest::check(m.Equals(m2));
-    });
+        Matrix4x4Unit += Test("SwapRows", []
+        {
+            Matrix4x4 m(1,2,3,4, 5,6,7,8, 9,10,11,12, 13,14,15,16);
+            Matrix4x4 m2(13,14,15,16, 9,10,11,12, 5,6,7,8, 1,2,3,4);
+            m.SwapRows(0,3);
+            m.SwapRows(1,2);
+            jtest::check(m.Equals(m2));
+        });
 
-    TEST("Mat4x4::CtorCols", []
-    {
-        Matrix4x4 m(Vector4(1,2,3,4), Vector4(5,6,7,8), Vector4(9,10,11,12), Vector4(13,14,15,16));
-        Matrix4x4 m2(1,5,9,13, 2,6,10,14, 3,7,11,15, 4,8,12,16);
-        jtest::check(m.Equals(m2));
-    });
+        Matrix4x4Unit += Test("CtorCols", []
+        {
+            Matrix4x4 m(Vector4(1,2,3,4), Vector4(5,6,7,8), Vector4(9,10,11,12), Vector4(13,14,15,16));
+            Matrix4x4 m2(1,5,9,13, 2,6,10,14, 3,7,11,15, 4,8,12,16);
+            jtest::check(m.Equals(m2));
+        });
 
-    TEST("Mat4x4::CtorFromQuat", []
-    {
-        RNG rng;
-        // TODO: Multiple random passes
-        Quaternion q = Quaternion::RandomRotation(rng);
-        Matrix4x4 m(q);
+        Matrix4x4Unit += Test("CtorFromQuat", []
+        {
+            RNG rng;
+            // TODO: Multiple random passes
+            Quaternion q = Quaternion::RandomRotation(rng);
+            Matrix4x4 m(q);
 
-        Vector3 v = Vector3(-1, 5, 20.f);
-        Vector3 v1 = q * v;
-        Vector3 v2 = m.Transform(v); //m.TransformPos(v);
-        jtest::check(v1.Equals(v2));
-    });
+            Vector3 v = Vector3(-1, 5, 20.f);
+            Vector3 v1 = q * v;
+            Vector3 v2 = m.Transform(v); //m.TransformPos(v);
+            jtest::check(v1.Equals(v2));
+        });
 
 
-    TEST("Mat4x4::CtorFromQuatTrans", [] {});
-    TEST("Mat4x4::Translate", [] {});
-    TEST("Mat4x4::Scale", [] {});
-    TEST("Mat4x4::InverseOrthogonalUniformScale", [] {});
-    TEST("Mat4x4::InverseOrthonormal", [] {});
-    TEST("Mat4x4::DeterminantCorrectness", [] { });
-    TEST("Mat4x4::MulMat3x3", [] {});
+        Matrix4x4Unit += Test("CtorFromQuatTrans", [] {});
+        Matrix4x4Unit += Test("Translate", [] {});
+        Matrix4x4Unit += Test("Scale", [] {});
+        Matrix4x4Unit += Test("InverseOrthogonalUniformScale", [] {});
+        Matrix4x4Unit += Test("InverseOrthonormal", [] {});
+        Matrix4x4Unit += Test("DeterminantCorrectness", [] { });
+        Matrix4x4Unit += Test("MulMat3x3", [] {});
+    }
 
-
-    TEST("Mat4x4::AngleTypeRoundtripConversion", Matrix4x4_AngleTypeRoundtripConversion);
-
-    return 0;
+    inline void Run() {
+        Matrix4x4Unit.RunAll();
+    }
 }
 
 
@@ -143,3 +135,4 @@ int Matrix4x4Tests() {
 
 
 
+

tests/LinearAlgebra/QuaternionTests.hpp

@@ -1,77 +1,113 @@
 #pragma once
 
 #include <cmath>
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
 #include <J3ML/LinearAlgebra/Quaternion.hpp>
+#include <J3ML/LinearAlgebra/EulerAngle.hpp>
+#include <J3ML/Algorithm/RNG.hpp>
 
-using namespace J3ML::LinearAlgebra;
+jtest::Unit QuaternionUnit {"Quaternion"};
+namespace QuaternionTests {
 
-Quaternion PreciseSlerp(const Quaternion &a, const Quaternion& b, float t)
-{
-    double angle = a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
-
-    double sign = 1.0;
-
-    if (angle > 0) {
-        angle = -angle;
-        sign = -1.0;
-    }
-    double A;
-    double B;
-
-    if (angle <= 1.0) // perform spherical linear interpolation
+    // This is here to check the accuracy of the Slerp inside the Quaternion class.
+    // Although you don't jtest::check anything :shrug: - Redacted.
+    Quaternion PreciseSlerp(const Quaternion &a, const Quaternion& b, float t)
     {
-        angle = std::acos(angle); // After this, angle is in the range pi / 2 -> 0 as the original angle variable ranged from 0 -> 1.
+        double angle = a.x * b.x + a.y * b.y + a.z * b.z + a.w * b.w;
 
-        double angleT = t * angle;
+        double sign = 1.0;
 
-        double s[3] = { std::sin(angle), std::sin(angle - angleT), std::sin(angleT)};
-        double c = 1.0 / s[0];
+        if (angle > 0) {
+            angle = -angle;
+            sign = -1.0;
+        }
+        double A;
+        double B;
 
-        A = s[1] * c;
-        B = s[2] * c;
-    } else { // If angle is close to taking the denominator to zero, resort to linear interpolation (and normalization).
-        A = 1.0 - t;
-        B = t;
+        if (angle <= 1.0) // perform spherical linear interpolation
+        {
+            angle = std::acos(angle); // After this, angle is in the range pi / 2 -> 0 as the original angle variable ranged from 0 -> 1.
+
+            double angleT = t * angle;
+
+            double s[3] = { std::sin(angle), std::sin(angle - angleT), std::sin(angleT)};
+            double c = 1.0 / s[0];
+
+            A = s[1] * c;
+            B = s[2] * c;
+        } else { // If angle is close to taking the denominator to zero, resort to linear interpolation (and normalization).
+            A = 1.0 - t;
+            B = t;
+        }
+
+        Quaternion C;
+        C.x = (float)(a.x*A*sign + b.x*B);
+        C.y = (float)(a.y*A*sign + b.y*B);
+        C.z = (float)(a.z*A*sign + b.z*B);
+        C.w = (float)(a.w*A*sign + b.w*B);
+        return C.Normalized();
     }
 
-    Quaternion C;
-    C.x = (float)(a.x*A*sign + b.x*B);
-    C.y = (float)(a.y*A*sign + b.y*B);
-    C.z = (float)(a.z*A*sign + b.z*B);
-    C.w = (float)(a.w*A*sign + b.w*B);
-    return C.Normalized();
-}
+    inline void Define()
+    {
+        using namespace jtest;
+        using namespace J3ML::LinearAlgebra;
+        using J3ML::Algorithm::RNG;
 
+        QuaternionUnit += Test("SlerpPrecision", [] {
+            RNG rng;
+            float maxError = 0;
+            float maxLerpError = 0;
+            float magnitudeError = 0;
+            for (int i = 0; i < 10000; ++i)
+            {
+                Quaternion q = Quaternion::RandomRotation(rng);
+                Quaternion q2 = Quaternion::RandomRotation(rng);
+                float t = rng.Float01Incl();
 
-int QuaternionTests()
-{
+                Quaternion correct = PreciseSlerp(q, q2, t);
+                Quaternion fast = q.Slerp(q2, t);
 
-    TEST("Quat::SlerpPrecision", [] {
-        RNG rng;
-        float maxError = 0;
-        float maxLerpError = 0;
-        float magnitudeError = 0;
-        for (int i = 0; i < 10000; ++i)
-        {
-            Quaternion q = Quaternion::RandomRotation(rng);
-            Quaternion q2 = Quaternion::RandomRotation(rng);
-            float t = rng.Float01Incl();
+                magnitudeError = std::max(magnitudeError, std::abs(1.f - fast.LengthSquared()));
+                Quaternion lerp = q.Lerp(q2, t);
+            }
+        });
+        QuaternionUnit += Test("From_EulerAngleXYZ", [] {
+            Quaternion expected_result(0.1819093, 0.6706149, 0.1840604, 0.6952024);
+            EulerAngleXYZ e(10, 88, 20);
+            Quaternion from_euler(e);
 
-            Quaternion correct = PreciseSlerp(q, q2, t);
-            Quaternion fast = q.Slerp(q2, t);
+            jtest::check(Math::EqualAbs(expected_result.x, from_euler.x, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.y, from_euler.y, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.z, from_euler.z, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.w, from_euler.w, 1e-6f));
+        });
 
-            magnitudeError = std::max(magnitudeError, std::abs(1.f - fast.LengthSquared()));
-            Quaternion lerp = q.Lerp(q2, t);
-        }
-    });
-    TEST("Quat::Mat4x4Conversion", [] {});
-    TEST("Quat::MulOpQuat", [] {});
-    TEST("Quat::DivOpQuat", [] {});
-    TEST("Quat::Lerp", [] {});
-    TEST("Quat::RotateFromTo", [] {});
-    TEST("Quat::Transform", [] {});
+        QuaternionUnit += Test("From_AxisAngle", [] {
+            Quaternion expected_result(0.0579133, 0.0782044, 0.1765667, 0.9794664);
+            AxisAngle a({0.2872573, 0.3879036, 0.8757934}, 0.4059981);
+            Quaternion from_axis(a);
 
-    return 0;
+            jtest::check(Math::EqualAbs(expected_result.x, from_axis.x, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.y, from_axis.y, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.z, from_axis.z, 1e-6f));
+            jtest::check(Math::EqualAbs(expected_result.w, from_axis.w, 1e-6f));
+        });
+
+        QuaternionUnit += Test("Mat4x4Conversion", [] { throw("Not Implemented"); });
+        QuaternionUnit += Test("MulOpQuat", [] { throw("Not Implemented"); });
+        QuaternionUnit += Test("DivOpQuat", [] { throw("Not Implemented"); });
+        QuaternionUnit += Test("Lerp", [] { throw("Not Implemented"); });
+        QuaternionUnit += Test("RotateFromTo", [] { throw("Not Implemented"); });
+        QuaternionUnit += Test("Transform", [] { throw("Not Implemented"); });
+
+    }
+
+    inline void Run()
+    {
+        QuaternionUnit.RunAll();
+    }
 }
 
 

tests/LinearAlgebra/Vector2Tests.hpp

@@ -1,191 +1,201 @@
 
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
 #include <J3ML/LinearAlgebra/Vector2.hpp>
 
-using J3ML::LinearAlgebra::Vector2;
-
-int Vector2Tests()
+jtest::Unit Vector2Unit {"Vector2"};
+namespace Vector2Tests
 {
-    TEST("Vector2::Ctor_Default", []
+    inline void Define()
     {
-        // TODO: implement check_eq
-        jtest::check(Vector2() == Vector2::Zero);
-    });
+        using namespace jtest;
+        using J3ML::LinearAlgebra::Vector2;
 
-    TEST("Vector2::Ctor_XY", []
+        Vector2Unit += Test("Ctor_Default", []
+        {
+            // TODO: implement check_eq
+            jtest::check(Vector2() == Vector2::Zero);
+        });
+
+        Vector2Unit += Test("Ctor_XY", []
+        {
+            Vector2 vec (1, 0);
+            // TODO: implement check_eq
+            jtest::check(vec == Vector2::Right);
+        });
+
+        Vector2Unit += Test("Addition_Op", []
+        {
+            Vector2 A (1,1);
+            Vector2 B (2,2);
+
+            Vector2 C (3, 3);
+            // TODO: implement check_eq
+            jtest::check(A+B == C);
+        });
+
+        Vector2Unit += Test("Addition_Method", []
+        {
+            Vector2 A (2,2);
+            Vector2 B (2,2);
+
+            Vector2 C (4, 4);
+            // TODO: implement check_eq
+            jtest::check(A.Add(B) == C);
+        });
+
+        Vector2Unit += Test("Addition_Static", []
+        {
+            Vector2 A (3, 3);
+            Vector2 B (2, 2);
+
+            Vector2 C (5, 5);
+            // TODO: implement check_eq
+            jtest::check(Vector2::Add(A, B) == C);
+        });
+
+        Vector2Unit += Test("Subtract_Op", []
+        {
+            Vector2 A (1,1);
+            Vector2 B (2,2);
+
+            Vector2 C (-1, -1);
+
+            jtest::check(A-B == C);
+        });
+
+        Vector2Unit += Test("Subtract_Method", []
+        {
+            Vector2 A (1,1);
+            Vector2 B (2,2);
+
+            Vector2 C (-1, -1);
+
+            jtest::check(A.Sub(B) == C);
+        });
+
+        Vector2Unit += Test("Subtract_Static", []
+        {
+            Vector2 A (1,1);
+            Vector2 B (2,2);
+
+            Vector2 C (-1, -1);
+
+            jtest::check(Vector2::Sub(A, B) == C);
+        });
+
+        Vector2Unit += Test("Scalar_Multiplication", []
+        {
+            Vector2 A (5, 1);
+            float B = 0.5f;
+            Vector2 C (2.5f, .5f);
+            jtest::check(A*B == C);
+        });
+
+        Vector2Unit += Test("Size", []
+        {
+            jtest::check(sizeof(Vector2) == 8);
+        });
+
+        Vector2Unit += Test("NaN", []
+        {
+
+            jtest::check(Vector2::Zero != Vector2::NaN);
+            jtest::check(Vector2::Up != Vector2::NaN);
+            jtest::check(Vector2::Left != Vector2::NaN);
+            jtest::check(Vector2::Down != Vector2::NaN);
+            jtest::check(Vector2::Right != Vector2::NaN);
+            jtest::check(Vector2::NaN != Vector2::NaN);
+        });
+
+        Vector2Unit += Test("MarginOfError", []
+        {
+            Vector2 A (2,2);
+            Vector2 B (1.85, 1.85);
+
+            jtest::check(A.IsWithinMarginOfError(B, 0.5f));
+        });
+
+        Vector2Unit += Test("Min", []
+        {
+            Vector2 A (2,2);
+            Vector2 B (1.85, 1.85);
+
+            jtest::check( Vector2::Min(A, B) == B);
+        });
+
+        Vector2Unit += Test("Max", []
+        {
+            Vector2 A (2,2);
+            Vector2 B (1.85, 1.85);
+
+            jtest::check( Vector2::Max(A, B) == A);
+        });
+
+        Vector2Unit += Test("Clamp", []
+        {
+            Vector2 Input (0, 20);
+
+            Vector2 Minimum ( 2, 2);
+            Vector2 Maximum (16, 16);
+
+            Vector2 ExpectedResult (2, 16);
+
+            jtest::check(Input.Clamp(Minimum, Maximum) == ExpectedResult);
+        });
+
+        Vector2Unit += Test("DotProduct", []
+        {
+
+            Vector2 A (2, 2);
+            Vector2 B (1, 1);
+            jtest::check(A.Dot(B) == 4.f);
+        });
+
+        Vector2Unit += Test("Project", []
+        {
+            Vector2 Base (4, 4);
+            Vector2 Projected (1, 0);
+
+            Vector2 ExpectedResult (4, 0);
+
+            jtest::check(Base.Project(Projected) == ExpectedResult);
+        });
+
+        Vector2Unit += Test("Normalized", []
+        {
+            Vector2 A(2, 0);
+            Vector2 B(1, 0);
+            jtest::check(A.Normalized() == B);
+        });
+
+        Vector2Unit += Test("Lerp", []
+        {
+            Vector2 A (2,2);
+            Vector2 B (10, 10);
+            Vector2 C (6, 6);
+
+            jtest::check(A.Lerp(B, 0.f) == A);
+            jtest::check(A.Lerp(B, 1.f) == B);
+            jtest::check(A.Lerp(B, 0.5f) == C);
+        });
+
+        Vector2Unit += Test("AngleBetween", []
+        {
+            Vector2 A (0.5f, 0.5);
+            Vector2 B (0.5f, 0.1f);
+
+            A = A.Normalized();
+            B = B.Normalized();
+
+            jtest::check_float_eq(A.AngleBetween(B), 0.58800244);
+        });
+    }
+
+    inline void Run()
     {
-        Vector2 vec (1, 0);
-        // TODO: implement check_eq
-        jtest::check(vec == Vector2::Right);
-    });
-
-    TEST("Vector2::Addition_Op", []
-    {
-        Vector2 A (1,1);
-        Vector2 B (2,2);
-
-        Vector2 C (3, 3);
-        // TODO: implement check_eq
-        jtest::check(A+B == C);
-    });
-
-    TEST("Vector2::Addition_Method", []
-    {
-        Vector2 A (2,2);
-        Vector2 B (2,2);
-
-        Vector2 C (4, 4);
-        // TODO: implement check_eq
-        jtest::check(A.Add(B) == C);
-    });
-
-    TEST("Vector2::Addition_Static", []
-    {
-        Vector2 A (3, 3);
-        Vector2 B (2, 2);
-
-        Vector2 C (5, 5);
-        // TODO: implement check_eq
-        jtest::check(Vector2::Add(A, B) == C);
-    });
-
-    TEST("Vector2::Subtract_Op", []
-    {
-        Vector2 A (1,1);
-        Vector2 B (2,2);
-
-        Vector2 C (-1, -1);
-
-        jtest::check(A-B == C);
-    });
-
-    TEST("Vector2::Subtract_Method", []
-    {
-        Vector2 A (1,1);
-        Vector2 B (2,2);
-
-        Vector2 C (-1, -1);
-
-        jtest::check(A.Sub(B) == C);
-    });
-
-    TEST("Vector2::Subtract_Static", []
-    {
-        Vector2 A (1,1);
-        Vector2 B (2,2);
-
-        Vector2 C (-1, -1);
-
-        jtest::check(Vector2::Sub(A, B) == C);
-    });
-
-    TEST("Vector2::Scalar_Multiplication", []
-    {
-        Vector2 A (5, 1);
-        float B = 0.5f;
-        Vector2 C (2.5f, .5f);
-        jtest::check(A*B == C);
-    });
-
-    TEST("Vector2::Size", []
-    {
-        jtest::check(sizeof(Vector2) == 8);
-    });
-
-    TEST("Vector2::NaN", []
-    {
-
-        jtest::check(Vector2::Zero != Vector2::NaN);
-        jtest::check(Vector2::Up != Vector2::NaN);
-        jtest::check(Vector2::Left != Vector2::NaN);
-        jtest::check(Vector2::Down != Vector2::NaN);
-        jtest::check(Vector2::Right != Vector2::NaN);
-        jtest::check(Vector2::NaN != Vector2::NaN);
-    });
-
-    TEST("Vector2::MarginOfError", []
-    {
-        Vector2 A (2,2);
-        Vector2 B (1.85, 1.85);
-
-        jtest::check(A.IsWithinMarginOfError(B, 0.5f));
-    });
-
-    TEST("Vector2::Min", []
-    {
-        Vector2 A (2,2);
-        Vector2 B (1.85, 1.85);
-
-        jtest::check( Vector2::Min(A, B) == B);
-    });
-
-    TEST("Vector2::Max", []
-    {
-        Vector2 A (2,2);
-        Vector2 B (1.85, 1.85);
-
-        jtest::check( Vector2::Max(A, B) == A);
-    });
-
-    TEST("Vector2::Clamp", []
-    {
-        Vector2 Input (0, 20);
-
-        Vector2 Minimum ( 2, 2);
-        Vector2 Maximum (16, 16);
-
-        Vector2 ExpectedResult (2, 16);
-
-        jtest::check(Input.Clamp(Minimum, Maximum) == ExpectedResult);
-    });
-
-    TEST("Vector2::DotProduct", []
-    {
-
-        Vector2 A (2, 2);
-        Vector2 B (1, 1);
-        jtest::check(A.Dot(B) == 4.f);
-    });
-
-    TEST("Vector2::Project", []
-    {
-        Vector2 Base (4, 4);
-        Vector2 Projected (1, 0);
-
-        Vector2 ExpectedResult (4, 0);
-
-        jtest::check(Base.Project(Projected) == ExpectedResult);
-    });
-
-    TEST("Vector2::Normalize", []
-    {
-        Vector2 A(2, 0);
-        Vector2 B(1, 0);
-        jtest::check(A.Normalized() == B);
-    });
-
-    TEST("Vector2::Lerp", []
-    {
-        Vector2 A (2,2);
-        Vector2 B (10, 10);
-        Vector2 C (6, 6);
-
-        jtest::check(A.Lerp(B, 0.f) == A);
-        jtest::check(A.Lerp(B, 1.f) == B);
-        jtest::check(A.Lerp(B, 0.5f) == C);
-    });
-
-    TEST("Vector2::AngleBetween", []
-    {
-        Vector2 A (0.5f, 0.5);
-        Vector2 B (0.5f, 0.1f);
-
-        A = A.Normalized();
-        B = B.Normalized();
-
-        // TODO: AngleBetween returns not a number
-        // TODO: implement jtest::check_float_eq
-        jtest::check(A.AngleBetween(B) == 0.58800244);
-    });
-
-    return 0;
+        Vector2Unit.RunAll();
+    }
 }
+
+

tests/LinearAlgebra/Vector3Tests.hpp

@@ -1,206 +1,215 @@
-
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
 #include <J3ML/LinearAlgebra/Vector3.hpp>
 
-using J3ML::LinearAlgebra::Vector3;
-
-inline void EXPECT_V3_EQ(const Vector3& lhs, const Vector3& rhs)
-{
-    jtest::check(lhs.x == rhs.x);
-    jtest::check(lhs.y == rhs.y);
-    jtest::check(lhs.z == rhs.z);
-}
-
-int Vector3Tests() {
-    TEST("Vector3::Ctor_Default", []
+jtest::Unit Vector3Unit {"Vector3"};
+namespace Vector3Tests {
+    inline void EXPECT_V3_EQ(const Vector3& lhs, const Vector3& rhs)
     {
-        EXPECT_V3_EQ(Vector3(), Vector3::Zero);
-    });
+        jtest::check(lhs.x == rhs.x);
+        jtest::check(lhs.y == rhs.y);
+        jtest::check(lhs.z == rhs.z);
+    }
 
-    TEST("Vector3::Ctor_XYZ", []
+    inline void Define() {
+        using namespace jtest;
+        using J3ML::LinearAlgebra::Vector3;
+
+        Vector3Unit += Test("Ctor_Default", []
+        {
+            EXPECT_V3_EQ(Vector3(), Vector3::Zero);
+        });
+
+        Vector3Unit += Test("Ctor_XYZ", []
+        {
+            Vector3 Input (0, 1, 0);
+
+            EXPECT_V3_EQ(Input, Vector3::Down);
+        });
+
+        Vector3Unit += Test("Addition_Op", [] {
+            Vector3 A (1,1,1);
+            Vector3 B (2,2,2);
+
+            Vector3 ExpectedResult (3,3,3);
+
+            EXPECT_V3_EQ(A + B, ExpectedResult);
+        });
+
+        Vector3Unit += Test("Addition_Method", [] {
+            Vector3 A (1,1,1);
+            Vector3 B (2,2,2);
+
+            Vector3 ExpectedResult (3,3,3);
+
+            EXPECT_V3_EQ(A.Add(B), ExpectedResult);
+        });
+        Vector3Unit += Test("Addition_Static", [] {
+            Vector3 A (1,1,1);
+            Vector3 B (3,3,3);
+
+            Vector3 ExpectedResult (4,4,4);
+
+            EXPECT_V3_EQ(Vector3::Add(A, B), ExpectedResult);
+        });
+        Vector3Unit += Test("Subtract_Op", [] {
+            Vector3 A (2,2,2);
+            Vector3 B (.5f, .5f, .5f);
+
+            Vector3 ExpectedResult (1.5f, 1.5f, 1.5f);
+
+            EXPECT_V3_EQ(A - B, ExpectedResult);
+        });
+        Vector3Unit += Test("Subtract_Method", [] {
+            Vector3 A (3,3,3);
+            Vector3 B (1,1,1);
+
+            Vector3 ExpectedResult (2,2,2);
+
+            EXPECT_V3_EQ(A.Sub(B), ExpectedResult);
+        });
+        Vector3Unit += Test("V3_Subtract_Static", [] {
+            Vector3 A (4,4,4);
+            Vector3 B (1,1,1);
+
+            Vector3 ExpectedResult (3,3,3);
+
+            EXPECT_V3_EQ(Vector3::Sub(A, B), ExpectedResult);
+        });
+        Vector3Unit += Test("Scalar_Mult_Op", [] {
+            Vector3 A ( 1,1,1);
+            float B = 1.5f;
+
+            Vector3 ExpectedResult (1.5f, 1.5f, 1.5f);
+
+            EXPECT_V3_EQ(A * B, ExpectedResult);
+        });
+        Vector3Unit += Test("Scalar_Mult_Method", [] {
+            Vector3 A (3,3,3);
+            float B = 1.5f;
+
+            Vector3 ExpectedResult (4.5f, 4.5f, 4.5f);
+
+            EXPECT_V3_EQ(A.Mul(B), ExpectedResult);
+        });
+        Vector3Unit += Test("Scalar_Mult_Static", [] {
+            Vector3 A (2,2,2);
+            float B = 1.5f;
+
+            Vector3 ExpectedResult (3.f, 3.f, 3.f);
+
+            EXPECT_V3_EQ(Vector3::Mul(A, B), ExpectedResult);
+        });
+        Vector3Unit += Test("Scalar_Div_Op", [] {
+            Vector3 A (4,4,4);
+            float B = 2.f;
+
+            Vector3 ExpectedResult (2,2,2);
+            EXPECT_V3_EQ(A / B, ExpectedResult);
+        });
+        Vector3Unit += Test("Scalar_Div_Method", [] {
+            Vector3 A (6,6,6);
+            float B = 2.f;
+            Vector3 ExpectedResult ( 3,3,3);
+
+            EXPECT_V3_EQ(A.Div(B), ExpectedResult);
+        });
+        Vector3Unit += Test("Scalar_Div_Static", [] {
+            Vector3 A (3,3,3);
+            float B = 1.5f;
+
+            Vector3 ExpectedResult ( 2.f, 2.f, 2.f);
+
+            EXPECT_V3_EQ(Vector3::Div(A, B), ExpectedResult);
+        });
+        Vector3Unit += Test("Sizeof", [] {
+            jtest::check(sizeof(Vector3) == 12);
+        });
+        Vector3Unit += Test("NaN", [] {
+            jtest::check(Vector3(0, 0, 0) != Vector3::NaN);
+
+        });
+        Vector3Unit += Test("Min", [] {
+            Vector3 Input (2,2,2);
+            Vector3 Minimum (3,3,3);
+            Vector3 ExpectedResult (2,2,2);
+
+            EXPECT_V3_EQ(Input.Min(Minimum), ExpectedResult);
+        });
+        Vector3Unit += Test("Max", [] {
+            Vector3 Input (2,2,2);
+            Vector3 Maximum (3,3,3);
+            Vector3 ExpectedResult (3,3,3);
+
+            EXPECT_V3_EQ(Input.Max(Maximum), ExpectedResult);
+        });
+        Vector3Unit += Test("Clamp", [] {
+            Vector3 Input (5,-1,8);
+            Vector3 Minimum (1,1,1);
+            Vector3 Maximum (5,5,5);
+
+            Vector3 ExpectedResult (5,1,5);
+
+            EXPECT_V3_EQ(Input.Clamp(Minimum, Maximum), ExpectedResult);
+
+        });
+        Vector3Unit += Test("DotProduct", [] {
+            Vector3 A(6,6,6);
+            Vector3 B(1,1,1);
+
+
+            float ExpectedResult = 18;
+
+            // TODO: Add check_float_eq
+            jtest::check(A.Dot(B) == ExpectedResult);
+        });
+        Vector3Unit += Test("CrossProduct", [] {
+            Vector3 A(1,1,1);
+            Vector3 B(2,2,2);
+
+            Vector3 ExpectedResult (0,0,0);
+
+            EXPECT_V3_EQ(A.Cross(B), ExpectedResult);
+
+        });
+        Vector3Unit += Test("Project", [] {
+            Vector3 Base;
+            Vector3 Projection;
+            Vector3 ExpectedResult;
+        });
+        Vector3Unit += Test("Normalized", [] {
+            Vector3 Input (2, 0, 0);
+            Vector3 ExpectedResult (1, 0, 0);
+
+            EXPECT_V3_EQ(Input.Normalized(), ExpectedResult);
+        });
+        Vector3Unit += Test("Lerp", []
+        {
+            Vector3 Start;
+            Vector3 Finish;
+            float Percent = 50;
+            Vector3 ExpectedResult;
+            EXPECT_V3_EQ(Start.Lerp(Finish, Percent), ExpectedResult);
+        });
+        Vector3Unit += Test("AngleBetween", [] {
+            using J3ML::LinearAlgebra::Angle2D;
+            Vector3 A ( .5f, .5f, .5f);
+            Vector3 B (.25f, .75f, .25f);
+            A = A.Normalized();
+            B = B.Normalized();
+            Angle2D ExpectedResult (-0.69791365, -2.3561945);
+            //std::cout << A.AngleBetween(B).x << ", " << A.AngleBetween(B).y << "";
+            auto angle = A.AngleBetween(B);
+
+            jtest::check(angle.x == ExpectedResult.x);
+            jtest::check(angle.y == ExpectedResult.y);
+        });
+    }
+
+    inline void Run()
     {
-        Vector3 Input (0, 1, 0);
-
-        EXPECT_V3_EQ(Input, Vector3::Down);
-    });
-
-    TEST("Vector3::Addition_Op", [] {
-        Vector3 A (1,1,1);
-        Vector3 B (2,2,2);
-
-        Vector3 ExpectedResult (3,3,3);
-
-        EXPECT_V3_EQ(A + B, ExpectedResult);
-    });
-
-    TEST("Vector3::Addition_Method", [] {
-        Vector3 A (1,1,1);
-        Vector3 B (2,2,2);
-
-        Vector3 ExpectedResult (3,3,3);
-
-        EXPECT_V3_EQ(A.Add(B), ExpectedResult);
-    });
-    TEST("Vector3::Addition_Static", [] {
-        Vector3 A (1,1,1);
-        Vector3 B (3,3,3);
-
-        Vector3 ExpectedResult (4,4,4);
-
-        EXPECT_V3_EQ(Vector3::Add(A, B), ExpectedResult);
-    });
-    TEST("Vector3::Subtract_Op", [] {
-        Vector3 A (2,2,2);
-        Vector3 B (.5f, .5f, .5f);
-
-        Vector3 ExpectedResult (1.5f, 1.5f, 1.5f);
-
-        EXPECT_V3_EQ(A - B, ExpectedResult);
-    });
-    TEST("Vector3::Subtract_Method", [] {
-        Vector3 A (3,3,3);
-        Vector3 B (1,1,1);
-
-        Vector3 ExpectedResult (2,2,2);
-
-        EXPECT_V3_EQ(A.Sub(B), ExpectedResult);
-    });
-    TEST("Vector3Test, V3_Subtract_Static", [] {
-        Vector3 A (4,4,4);
-        Vector3 B (1,1,1);
-
-        Vector3 ExpectedResult (3,3,3);
-
-        EXPECT_V3_EQ(Vector3::Sub(A, B), ExpectedResult);
-    });
-    TEST("Vector3::Scalar_Mult_Op", [] {
-        Vector3 A ( 1,1,1);
-        float B = 1.5f;
-
-        Vector3 ExpectedResult (1.5f, 1.5f, 1.5f);
-
-        EXPECT_V3_EQ(A * B, ExpectedResult);
-    });
-    TEST("Vector3::Scalar_Mult_Method", [] {
-        Vector3 A (3,3,3);
-        float B = 1.5f;
-
-        Vector3 ExpectedResult (4.5f, 4.5f, 4.5f);
-
-        EXPECT_V3_EQ(A.Mul(B), ExpectedResult);
-    });
-    TEST("Vector3::Scalar_Mult_Static", [] {
-        Vector3 A (2,2,2);
-        float B = 1.5f;
-
-        Vector3 ExpectedResult (3.f, 3.f, 3.f);
-
-        EXPECT_V3_EQ(Vector3::Mul(A, B), ExpectedResult);
-    });
-    TEST("Vector3::Scalar_Div_Op", [] {
-        Vector3 A (4,4,4);
-        float B = 2.f;
-
-        Vector3 ExpectedResult (2,2,2);
-        EXPECT_V3_EQ(A / B, ExpectedResult);
-    });
-    TEST("Vector3::Scalar_Div_Method", [] {
-        Vector3 A (6,6,6);
-        float B = 2.f;
-        Vector3 ExpectedResult ( 3,3,3);
-
-        EXPECT_V3_EQ(A.Div(B), ExpectedResult);
-    });
-    TEST("Vector3::Scalar_Div_Static", [] {
-        Vector3 A (3,3,3);
-        float B = 1.5f;
-
-        Vector3 ExpectedResult ( 2.f, 2.f, 2.f);
-
-        EXPECT_V3_EQ(Vector3::Div(A, B), ExpectedResult);
-    });
-    TEST("Vector3::Sizeof", [] {
-        jtest::check(sizeof(Vector3) == 12);
-    });
-    TEST("Vector3::NaN", [] {
-        jtest::check(Vector3(0, 0, 0) != Vector3::NaN);
-
-    });
-    TEST("Vector3::Min", [] {
-        Vector3 Input (2,2,2);
-        Vector3 Minimum (3,3,3);
-        Vector3 ExpectedResult (2,2,2);
-
-        EXPECT_V3_EQ(Input.Min(Minimum), ExpectedResult);
-    });
-    TEST("Vector3::Max", [] {
-        Vector3 Input (2,2,2);
-        Vector3 Maximum (3,3,3);
-        Vector3 ExpectedResult (3,3,3);
-
-        EXPECT_V3_EQ(Input.Max(Maximum), ExpectedResult);
-    });
-    TEST("Vector3::Clamp", [] {
-        Vector3 Input (5,-1,8);
-        Vector3 Minimum (1,1,1);
-        Vector3 Maximum (5,5,5);
-
-        Vector3 ExpectedResult (5,1,5);
-
-        EXPECT_V3_EQ(Input.Clamp(Minimum, Maximum), ExpectedResult);
-
-    });
-    TEST("Vector3::DotProduct", [] {
-        Vector3 A(6,6,6);
-        Vector3 B(1,1,1);
-
-
-        float ExpectedResult = 18;
-
-        // TODO: Add check_float_eq
-        jtest::check(A.Dot(B) == ExpectedResult);
-    });
-    TEST("Vector3::CrossProduct", [] {
-        Vector3 A(1,1,1);
-        Vector3 B(2,2,2);
-
-        Vector3 ExpectedResult (0,0,0);
-
-        EXPECT_V3_EQ(A.Cross(B), ExpectedResult);
-
-    });
-    TEST("Vector3::Project", [] {
-        Vector3 Base;
-        Vector3 Projection;
-        Vector3 ExpectedResult;
-    });
-    TEST("Vector3::Normalize", [] {
-        Vector3 Input (2, 0, 0);
-        Vector3 ExpectedResult (1, 0, 0);
-
-        EXPECT_V3_EQ(Input.Normalized(), ExpectedResult);
-    });
-    TEST("Vector3::Lerp", []
-    {
-        Vector3 Start;
-        Vector3 Finish;
-        float Percent = 50;
-        Vector3 ExpectedResult;
-        EXPECT_V3_EQ(Start.Lerp(Finish, Percent), ExpectedResult);
-    });
-    TEST("Vector3::AngleBetween", [] {
-        using J3ML::LinearAlgebra::Angle2D;
-        Vector3 A ( .5f, .5f, .5f);
-        Vector3 B (.25f, .75f, .25f);
-        A = A.Normalized();
-        B = B.Normalized();
-        Angle2D ExpectedResult (-0.69791365, -2.3561945);
-        //std::cout << A.AngleBetween(B).x << ", " << A.AngleBetween(B).y << "";
-        auto angle = A.AngleBetween(B);
-
-        jtest::check(angle.x == ExpectedResult.x);
-        jtest::check(angle.y == ExpectedResult.y);
-    });
-    return 0;
+        Vector3Unit.RunAll();
+    }
 }
 
 

tests/LinearAlgebra/Vector4Tests.hpp

@@ -1,41 +1,55 @@
 
 #include <J3ML/LinearAlgebra/Vector4.hpp>
+#include <jtest/jtest.hpp>
+#include <jtest/Unit.hpp>
 
 using Vector4 = J3ML::LinearAlgebra::Vector4;
 
-void EXPECT_V4_EQ(const Vector4& lhs, const Vector4& rhs)
+namespace jtest
 {
-    jtest::check(lhs.x == rhs.x);
-    jtest::check(lhs.y == rhs.y);
-    jtest::check(lhs.z == rhs.z);
-    jtest::check(lhs.w == rhs.w);
+    void check_v4_eq(const Vector4& lhs, const Vector4& rhs)
+    {
+        check_float_eq(lhs.x, rhs.x);
+        check_float_eq(lhs.y, rhs.y);
+        check_float_eq(lhs.z, rhs.z);
+        check_float_eq(lhs.w, rhs.w);
+    }
 }
 
-int Vector4Tests()
+jtest::Unit Vector4Unit {"Vector4"};
+namespace Vector4Tests
 {
-    TEST("Vector4::Ctor_Default", []
+    inline void Define()
     {
-        EXPECT_V4_EQ(Vector4(), Vector4::Zero);
-    });
+        using namespace jtest;
 
-    TEST("Vector4::Constructor_XYZW", []
+        Vector4Unit += Test("Vector4()", []{
+           jtest::check_v4_eq(Vector4(), Vector4::Zero);
+        });
+
+        Vector4Unit += Test("Vector4(float x, float y, float z, float) w", []{
+            Vector4 Input (0, 1, 0, 1);
+            jtest::check_float_eq(Input.x, 0);
+            jtest::check_float_eq(Input.y, 1);
+            jtest::check_float_eq(Input.z, 0);
+            jtest::check_float_eq(Input.w, 1);
+        });
+
+        Vector4Unit += Test("Vector4::Addition_Op", [] {
+            Vector4 A (1, 1, 1, 1);
+            Vector4 B (2, 2, 2, 2);
+
+            Vector4 ExpectedResult (3, 3, 3, 3);
+
+            jtest::check_v4_eq(A + B, ExpectedResult);
+        });
+    }
+
+    inline void Run()
     {
-        Vector4 Input (0, 1, 0, 1);
-        // TODO: implement jtest::check_float_eq
-        jtest::check(Input.x == 0);
-        jtest::check(Input.y == 1);
-        jtest::check(Input.z == 0);
-        jtest::check(Input.w == 1);
-    });
-
-    TEST("Vector4::Addition_Op", [] {
-        Vector4 A (1, 1, 1, 1);
-        Vector4 B (2, 2, 2, 2);
-
-        Vector4 ExpectedResult (3, 3, 3, 3);
-
-        EXPECT_V4_EQ(A + B, ExpectedResult);
-    });
+        Vector4Unit.RunAll();
+    }
+}
 
     /*TEST(Vector4Test, V4_Addition_Method)
     {
@@ -108,9 +122,4 @@ int Vector4Tests()
     TEST(Vector4Tests, V4_DotProduct) {}
     TEST(Vector4Tests, V4_CrossProduct) {}
     TEST(Vector4Tests, V4_Project) {}
-    TEST(Vector4Test, V4_Normalize) {}*/
-    return 0;
-}
-
-
-
+    TEST(Vector4Test, V4_Normalize) {}*/

tests/MathFuncTests.hpp

@@ -13,71 +13,80 @@
 #pragma once
 
 #include <jtest/jtest.hpp>
-#include "J3ML/J3ML.hpp"
+#include <jtest/Unit.hpp>
+#include <J3ML/J3ML.hpp>
 
-inline void MathFuncTests()
+jtest::Unit MathUnit {"Math"};
+namespace MathTests
 {
-    using namespace jtest;
-    using namespace J3ML::Math;
+    inline void Define()
+    {
+        using namespace jtest;
+        using namespace J3ML::Math;
+        using namespace J3ML;
 
-    TEST("Math::IsFinite", [] {
-        check(IsFinite<u32>(420));
-        check(IsFinite<u64>(25));
-        check(IsFinite(5.f));
-        check(IsFinite(5.0));
-        check(!IsFinite(Infinity));
-        check(!IsFinite(NotANumber));
-    });
+        MathUnit += Test("IsFinite", [] {
+            check(IsFinite<u32>(420));
+            check(IsFinite<u64>(25));
+            check(IsFinite(5.f));
+            check(IsFinite(5.0));
+            check(!IsFinite(Infinity));
+            check(!IsFinite(NotANumber));
+        });
 
-    TEST("Math::IsNotANumber", [] {
-        check(!IsNotANumber(5.f));
-        check(!IsNotANumber(5.0));
-        check(IsNotANumber(NotANumber));
-        check(!IsNotANumber(Infinity));
-    });
+        MathUnit += Test("IsNotANumber", [] {
+            check(!IsNotANumber(5.f));
+            check(!IsNotANumber(5.0));
+            check(IsNotANumber(NotANumber));
+            check(!IsNotANumber(Infinity));
+        });
 
-    TEST("Math::IsInfinite", [] {
-        check(!IsInfinite(5.f));
-        check(!IsInfinite(5.0));
-        check(IsInfinite(Infinity));
-    });
+        MathUnit += Test("IsInfinite", [] {
+            check(!IsInfinite(5.f));
+            check(!IsInfinite(5.0));
+            check(IsInfinite(Infinity));
+        });
 
-    TEST("Math::ReinterpretAsU32", [] {
-        check(ReinterpretAs<u32>(0.f)       == 0x00000000);
-        check(ReinterpretAs<u32>(1.f)       == 0x3F800000);
-        check(ReinterpretAs<u32>(2.f)       == 0x40000000);
-        check(ReinterpretAs<u32>(-1.f)      == 0xBF800000);
-        check(ReinterpretAs<u32>(Infinity)  == 0x7F800000);
-    });
+        MathUnit += Test("ReinterpretAsU32", [] {
+            check(ReinterpretAs<u32>(0.f)       == 0x00000000);
+            check(ReinterpretAs<u32>(1.f)       == 0x3F800000);
+            check(ReinterpretAs<u32>(2.f)       == 0x40000000);
+            check(ReinterpretAs<u32>(-1.f)      == 0xBF800000);
+            check(ReinterpretAs<u32>(Infinity)  == 0x7F800000);
+        });
 
-    TEST("Math::ReinterpretAsFloat", [] {
-        check(ReinterpretAs<float, u32>(0x00000000) == 0.f);
-        check(ReinterpretAs<float, u32>(0x3F800000) == 1.f);
-        check(ReinterpretAs<float, u32>(0x40000000) == 2.f);
-        check(ReinterpretAs<float, u32>(0xBF800000) == -1.f);
-        check(ReinterpretAs<float>(0x7F800000) == Infinity);
-        check(IsNotANumber(ReinterpretAs<float>(0x7F800001)));
-    });
+        MathUnit += Test("ReinterpretAsFloat", [] {
+            check(ReinterpretAs<float, u32>(0x00000000) == 0.f);
+            check(ReinterpretAs<float, u32>(0x3F800000) == 1.f);
+            check(ReinterpretAs<float, u32>(0x40000000) == 2.f);
+            check(ReinterpretAs<float, u32>(0xBF800000) == -1.f);
+            check(ReinterpretAs<float>(0x7F800000) == Infinity);
+            check(IsNotANumber(ReinterpretAs<float>(0x7F800001)));
+        });
 
-    TEST("Math::SqrtVals", [] {
+        MathUnit += Test("SqrtVals", [] {
 
-    });
+        });
 
-    TEST("Math::SqrtPrecision", [] {
+        MathUnit += Test("SqrtPrecision", [] {
 
-    });
+        });
 
-    TEST("Math::SqrtRSqrtPrecision", [] {
+        MathUnit += Test("SqrtRSqrtPrecision", [] {
 
-    });
+        });
 
-    TEST("Math::SqrtRecipPrecision", [] {
+        MathUnit += Test("SqrtRecipPrecision", [] {
 
-    });
+        });
 
-    TEST("Math::Min", [] {});
-    TEST("Math::Max", [] {});
-    TEST("Math::IsPow2", [] {});
+        MathUnit += Test("Min", [] {});
+        MathUnit += Test("Max", [] {});
+        MathUnit += Test("IsPow2", [] {});
+    }
 
-
-}
+    inline void Run()
+    {
+        MathUnit.RunAll();
+    }
+}

tests/tests.cpp

@@ -9,16 +9,14 @@
 /// @desc Includes and runs all project unit tests.
 /// @edit 2024-07-06
 
-
-
 #include "Algorithm/RNGTests.hpp"
-#include "Geometry/Geometry.hpp"
+#include "Geometry/CommonGeometryTests.hpp"
 #include "Geometry/TriangleTests.hpp"
 #include "Geometry/AABBTests.hpp"
 #include "Geometry/FrustumTests.hpp"
 
-
 #include "LinearAlgebra/EulerAngleTests.hpp"
+#include "LinearAlgebra/AxisAngleTests.hpp"
 #include "LinearAlgebra/Matrix2x2Tests.hpp"
 #include "LinearAlgebra/Matrix3x3Tests.hpp"
 #include "LinearAlgebra/Matrix4x4Tests.hpp"
@@ -32,44 +30,110 @@
 // TODO: Fix this once Automatic Test Discovery is implemented in jtest
 // TODO: Add test category blocks as feature to jtest
 
-void MiscTests() {
-    MathFuncTests();
+
+//void GeometryTests() {
+    //CommonGeometryTests();
+    //TriangleTests();
+    //AABBTests();
+    //FrustumTests();
+//}
+
+//void LinearAlgebraTests() {
+    //EulerAngleTests();
+    //Vector2Tests();
+    //Vector3Tests();
+    //Vector4Tests();
+    //QuaternionTests();
+    //Matrix2x2Tests();
+    //Matrix3x3Tests();
+    //Matrix4x4Tests();
+//}
+
+namespace AlgorithmTests
+{
+    void Define()
+    {
+        RNGTests::Define();
+    }
+    void Run()
+    {
+        RNGTests::Run();
+    }
 }
 
-void AlgorithmTests() {
-    RNGTests();
+namespace LinearAlgebraTests
+{
+    void Define()
+    {
+        Vector2Tests::Define();
+        Vector3Tests::Define();
+        Vector4Tests::Define();
+        AxisAngleTests::Define();
+        EulerAngleTests::Define();
+        QuaternionTests::Define();
+        Matrix2x2Tests::Define();
+        Matrix3x3Tests::Define();
+        Matrix4x4Tests::Define();
+
+    }
+    void Run()
+    {
+        Vector2Tests::Run();
+        Vector3Tests::Run();
+        Vector4Tests::Run();
+        AxisAngleTests::Run();
+        EulerAngleTests::Run();
+        QuaternionTests::Run();
+        Matrix2x2Tests::Run();
+        Matrix3x3Tests::Run();
+        Matrix4x4Tests::Run();
+    }
 }
 
-void GeometryTests() {
-    CommonGeometryTests();
-    TriangleTests();
-    AABBTests();
-    FrustumTests();
+namespace GeometryTests
+{
+    void Define()
+    {
+        CommonGeometryTests::Define();
+        TriangleTests::Define();
+        AABBTests::Define();
+        FrustumTests::Define();
+    }
+    void Run()
+    {
+        CommonGeometryTests::Run();
+        TriangleTests::Run();
+        AABBTests::Run();
+        FrustumTests::Run();
+    }
 }
 
-void LinearAlgebraTests() {
-    EulerAngleTests();
-    Vector2Tests();
-    Vector3Tests();
-    Vector4Tests();
-    QuaternionTests();
-    Matrix2x2Tests();
-    Matrix3x3Tests();
-    Matrix4x4Tests();
+
+
+void DefineTests()
+{
+    MathTests::Define();
+    AlgorithmTests::Define();
+    LinearAlgebraTests::Define();
+    GeometryTests::Define();
 }
 
+void RunTests()
+{
+    MathTests::Run();
+    AlgorithmTests::Run();
+    LinearAlgebraTests::Run();
+    GeometryTests::Run();
+}
 
 
 int main(int argc, char** argv)
 {
-    //MiscTests();
-    //AlgorithmTests();
-    //LinearAlgebraTests();
-    GeometryTests();
 
+    DefineTests();
 
+    RunTests();
 
-    jtest::run_tests();
     return 0;
 }