Improve and consolidate quaternion struct. Add more quaternion operators

src/ai/steering/agent.hpp

@@ -21,7 +21,7 @@
 #define ANTKEEPER_AI_STEERING_AGENT_HPP
 
 #include "utility/fundamental-types.hpp"
-#include "math/quaternion-type.hpp"
+#include "math/quaternion.hpp"
 
 namespace ai {
 namespace steering {

src/ai/steering/behavior/wander.cpp

@@ -20,8 +20,7 @@
 #include "ai/steering/behavior/wander.hpp"
 #include "ai/steering/behavior/seek.hpp"
 #include "math/random.hpp"
-#include "math/quaternion-functions.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 
 namespace ai {
 namespace steering {

src/entity/commands.cpp

@@ -23,9 +23,7 @@
 #include "game/component/parent.hpp"
 #include "game/component/celestial-body.hpp"
 #include "game/component/terrain.hpp"
-#include "math/quaternion-type.hpp"
-#include "math/quaternion-operators.hpp"
-#include "math/quaternion-functions.hpp"
+#include "math/quaternion.hpp"
 #include <limits>
 
 namespace entity {

src/game/ant/morphogenesis.cpp

@@ -21,7 +21,7 @@
 #include "render/material.hpp"
 #include "render/vertex-attribute.hpp"
 #include "math/transform-operators.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include <unordered_set>
 #include <iostream>
 

src/game/ant/swarm.cpp

@@ -23,7 +23,7 @@
 #include "game/component/model.hpp"
 #include "game/component/picking.hpp"
 #include "resources/resource-manager.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "config.hpp"
 #include <cmath>
 #include <random>
@@ -99,7 +99,7 @@ entity::id create_swarm(game::context& ctx)
 	steering.agent.max_force = 4.0f;
 	steering.agent.max_speed = 5.0f;
 	steering.agent.max_speed_squared = steering.agent.max_speed * steering.agent.max_speed;
-	steering.agent.orientation = math::quaternion<float>::identity;
+	steering.agent.orientation = math::quaternion<float>::identity();
 	steering.agent.forward = steering.agent.orientation * config::global_forward;
 	steering.agent.up = steering.agent.orientation * config::global_up;
 	steering.wander_weight = 1.0f;

src/game/state/splash.cpp

@@ -29,7 +29,7 @@
 #include "debug/logger.hpp"
 #include "resources/resource-manager.hpp"
 #include "render/material-flags.hpp"
-#include "math/matrix.hpp"
+#include "math/linear-algebra.hpp"
 #include <iostream>
 
 namespace game {
@@ -37,7 +37,7 @@ namespace state {
 
 splash::splash(game::context& ctx):
 	game::state::base(ctx)
-{
+{	
 	ctx.logger->push_task("Entering splash state");
 	
 	// Enable color buffer clearing in UI pass

src/game/system/astronomy.cpp

@@ -158,7 +158,7 @@ void astronomy::update(double t, double dt)
 		if (orbit.parent != entt::null)
 		{
 			transform.local.translation = math::normalize(float3(r_eus));
-			transform.local.rotation = math::type_cast<float>(rotation_eus);
+			transform.local.rotation = math::quaternion<float>(rotation_eus);
 			transform.local.scale = {1.0f, 1.0f, 1.0f};
 		}
 	});
@@ -581,7 +581,7 @@ void astronomy::update_icrf_to_eus(const game::component::celestial_body& body,
 			math::transformation::se3<float>
 			{
 				float3(icrf_to_eus.t),
-				math::type_cast<float>(icrf_to_eus.r)
+				math::quaternion<float>(icrf_to_eus.r)
 			}
 		);
 	}

src/game/system/constraint.cpp

@@ -19,7 +19,7 @@
 
 #include "game/system/constraint.hpp"
 #include "game/component/constraint-stack.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "math/transform-operators.hpp"
 
 namespace game {
@@ -311,13 +311,7 @@ void constraint::handle_spring_to_constraint(component::transform& transform, co
 			if (constraint.spring_rotation)
 			{
 				// Update rotation spring target
-				constraint.rotation.x1 =
-				{
-					target_transform->world.rotation.w,
-					target_transform->world.rotation.x,
-					target_transform->world.rotation.y,
-					target_transform->world.rotation.z
-				};
+				constraint.rotation.x1 = float4(target_transform->world.rotation);
 				
 				// Solve rotation spring
 				solve_numeric_spring<float4, float>(constraint.rotation, dt);

src/game/system/steering.cpp

@@ -23,7 +23,7 @@
 #include "entity/id.hpp"
 #include "ai/steering/behavior/wander.hpp"
 #include "ai/steering/behavior/seek.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "config.hpp"
 
 namespace game {

src/game/system/terrain.cpp

@@ -26,7 +26,7 @@
 #include "geom/quadtree.hpp"
 #include "geom/primitive/ray.hpp"
 #include "gl/vertex-attribute.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "render/vertex-attribute.hpp"
 #include "utility/fundamental-types.hpp"
 #include <functional>

src/game/system/terrain.hpp

@@ -23,7 +23,7 @@
 #include "game/system/updatable.hpp"
 #include "game/component/terrain.hpp"
 #include "entity/id.hpp"
-#include "math/quaternion-type.hpp"
+#include "math/quaternion.hpp"
 #include "geom/quadtree.hpp"
 #include "geom/mesh.hpp"
 #include "utility/fundamental-types.hpp"

src/math/linear-algebra.hpp

@@ -0,0 +1,27 @@
+/*
+ * Copyright (C) 2021  Christopher J. Howard
+ *
+ * This file is part of Antkeeper source code.
+ *
+ * Antkeeper source code is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Antkeeper source code is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef ANTKEEPER_MATH_LINEAR_ALGEBRA_HPP
+#define ANTKEEPER_MATH_LINEAR_ALGEBRA_HPP
+
+#include "math/matrix.hpp"
+#include "math/quaternion.hpp"
+#include "math/vector.hpp"
+
+#endif // ANTKEEPER_MATH_LINEAR_ALGEBRA_HPP

src/math/math.hpp

@@ -25,18 +25,13 @@ namespace math {}
 
 #include "math/vector.hpp"
 #include "math/matrix.hpp"
-
-#include "math/quaternion-type.hpp"
-#include "math/quaternion-functions.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 
 #include "math/se3.hpp"
 #include "math/transform-type.hpp"
 #include "math/transform-functions.hpp"
 #include "math/transform-operators.hpp"
 
-#include "math/stream-operators.hpp"
-
 #include "math/angles.hpp"
 #include "math/constants.hpp"
 #include "math/quadrature.hpp"

src/math/matrix.hpp

@@ -365,7 +365,7 @@ constexpr T determinant(const matrix& m) noexcept;
 /**
  * Calculates the inverse of a square matrix.
  *
- * @param m Matrix of which to take the inverse.
+ * @param m Square matrix.
  *
  * @return Inverse of matrix @p m.
  *
@@ -379,6 +379,8 @@ constexpr matrix inverse(const matrix& m) noexcept;
  *
  * @param x First matrix multiplicand.
  * @param y Second matrix multiplicand.
+ *
+ * @return Product of the component-wise multiplcation.
  */
 template <class T, std::size_t N, std::size_t M>
 constexpr matrix<T, N, M> componentwise_mul(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept;
@@ -388,6 +390,7 @@ constexpr matrix componentwise_mul(const matrix& a, const matr
  *
  * @param a First matrix.
  * @param b Second matrix.
+ *
  * @return Result of the division.
  */
 template <class T, std::size_t N, std::size_t M>
@@ -398,6 +401,7 @@ constexpr matrix div(const matrix& a, const matrix& b
  *
  * @param a Matrix.
  * @param b Scalar.
+ *
  * @return Result of the division.
  */
 template <class T, std::size_t N, std::size_t M>
@@ -408,6 +412,7 @@ constexpr matrix div(const matrix& a, T b) noexcept;
  *
  * @param a Scalar.
  * @param b Matrix.
+ *
  * @return Result of the division.
  */
 template <class T, std::size_t N, std::size_t M>
@@ -419,6 +424,7 @@ constexpr matrix div(T a, const matrix& b) noexcept;
  * @param position Position of the view point.
  * @param target Position of the target.
  * @param up Normalized direction of the up vector.
+ *
  * @return Viewing transformation matrix.
  */
 template <class T>
@@ -445,6 +451,7 @@ constexpr matrix mul(const matrix& a, const matrix& b
  *
  * @param a Matrix.
  * @param b Scalar.
+ *
  * @return Product of the matrix and the scalar.
  */
 template <class T, std::size_t N, std::size_t M>
@@ -476,7 +483,8 @@ constexpr typename matrix::row_vector_type mul(const typename matrix
  * Constructs a rotation matrix.
  *
  * @param angle Angle of rotation, in radians.
- * @param axis Axis of rotation
+ * @param axis Axis of rotation.
+ *
  * @return Rotation matrix.
  */
 template <class T>
@@ -486,6 +494,7 @@ matrix rotate(T angle, const vector& axis);
  * Produces a matrix which rotates Cartesian coordinates about the x-axis by a given angle.
  *
  * @param angle Angle of rotation, in radians.
+ *
  * @return Rotation matrix.
  */
 template <class T>
@@ -495,6 +504,7 @@ matrix3 rotate_x(T angle);
  * Produces a matrix which rotates Cartesian coordinates about the y-axis by a given angle.
  *
  * @param angle Angle of rotation, in radians.
+ *
  * @return Rotation matrix.
  */
 template <class T>
@@ -504,6 +514,7 @@ matrix3 rotate_y(T angle);
  * Produces a matrix which rotates Cartesian coordinates about the z-axis by a given angle.
  *
  * @param angle Angle of rotation, in radians.
+ *
  * @return Rotation matrix.
  */
 template <class T>
@@ -514,6 +525,7 @@ matrix3 rotate_z(T angle);
  *
  * @param m Matrix to scale.
  * @param v Scale vector.
+ *
  * @return Scaled matrix.
  */
 template <class T>
@@ -552,11 +564,22 @@ constexpr matrix sub(const matrix& a, T b) noexcept;
 template <class T, std::size_t N, std::size_t M>
 constexpr matrix<T, N, M> sub(T a, const matrix<T, N, M>& b) noexcept;
 
+/**
+ * Calculates the trace of a square matrix.
+ *
+ * @param m Square matrix.
+ *
+ * @return Sum of elements on the main diagonal.
+ */
+template <class T, std::size_t N>
+constexpr T trace(const matrix<T, N, N>& m) noexcept;
+
 /**
  * Translates a matrix.
  *
  * @param m Matrix to translate.
  * @param v Translation vector.
+ *
  * @return Translated matrix.
  */
 template <class T>
@@ -942,6 +965,19 @@ constexpr matrix sub(T a, const matrix& b) noexcept
 	return sub(a, b, std::make_index_sequence<N>{});
 }
 
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+constexpr inline T trace(const matrix<T, N, N>& m, std::index_sequence<I...>) noexcept
+{
+	return ((m[I][I]) + ...);
+}
+
+template <class T, std::size_t N>
+constexpr T trace(const matrix<T, N, N>& m) noexcept
+{
+	return trace(m, std::make_index_sequence<N>{});
+}
+
 template <class T>
 constexpr matrix<T, 4, 4> translate(const matrix<T, 4, 4>& m, const vector<T, 3>& v)
 {
@@ -1079,6 +1115,7 @@ constexpr inline matrix operator-(T a, const matrix& b) noexce
  *
  * @param a First value.
  * @param b Second value.
+ *
  * @return Reference to the first value.
  */
 /// @{
@@ -1099,6 +1136,7 @@ constexpr inline matrix& operator+=(matrix& a, T b) noexcept
  *
  * @param a First value.
  * @param b Second value.
+ *
  * @return Reference to the first value.
  */
 /// @{
@@ -1119,11 +1157,12 @@ constexpr inline matrix& operator-=(matrix& a, T b) noexcept
  *
  * @param a First value.
  * @param b Second value.
+ *
  * @return Reference to the first value.
  */
 /// @{
 template <class T, std::size_t N>
-constexpr inline matrix<T, N, N>& operator*=(matrix<T, N, N>& a, matrix<T, N, N>& b) noexcept
+constexpr inline matrix<T, N, N>& operator*=(matrix<T, N, N>& a, const matrix<T, N, N>& b) noexcept
 {
 	return (a = a * b);
 }
@@ -1139,6 +1178,7 @@ constexpr inline matrix& operator*=(matrix& a, T b) noexcept
  *
  * @param a First value.
  * @param b Second value.
+ *
  * @return Reference to the first value.
  */
 /// @{
@@ -1159,6 +1199,7 @@ constexpr inline matrix& operator/=(matrix& a, T b) noexcept
  *
  * @param os Output stream.
  * @param m Matrix.
+ *
  * @return Output stream.
  */
 template <class T, std::size_t N, std::size_t M>
@@ -1179,6 +1220,7 @@ std::ostream& operator<<(std::ostream& os, const matrix& m)
  *
  * @param is Input stream.
  * @param m Matrix.
+ *
  * @return Input stream.
  */
 template <class T, std::size_t N, std::size_t M>

src/math/operators.hpp

@@ -20,10 +20,6 @@
 #ifndef ANTKEEPER_MATH_OPERATORS_HPP
 #define ANTKEEPER_MATH_OPERATORS_HPP
 
-#include "math/quaternion-operators.hpp"
-#include "math/transform-operators.hpp"
-#include "math/stream-operators.hpp"
-
 namespace math {
 
 /// Mathematical operators.

src/math/polynomial.hpp

@@ -52,39 +52,6 @@ T horner(InputIt first, InputIt last, T x)
  */
 namespace chebyshev {
 
-	/**
-	 * Generates a Chebyshev approximation of a function.
-	 *
-	 * @param[out] first,last Range of Chebyshev polynomial coefficients.
-	 * @param[in] f Unary function to approximate.
-	 * @param[in] min,max Domain of @p f.
-	 */
-	template <class OutputIt, class UnaryOp, class T>
-	void approximate(OutputIt first, OutputIt last, UnaryOp f, T min, T max)
-	{
-		std::size_t n = last - first;
-		const T two_over_n = T(2) / static_cast<T>(n);
-		const T pi_over_n = math::pi<T> / static_cast<T>(n);
-		
-		last = first;
-		for (std::size_t i = 0; i < n; ++i)
-			*(last++) = T(0);
-		
-		for (std::size_t i = 0; i < n; ++i)
-		{
-			const T y = pi_over_n * (static_cast<T>(i) + T(0.5));
-			
-			T x = f(math::map<T>(std::cos(y), T(-1), T(1), min, max)) * two_over_n;
-			
-			*first += x;
-			last = first;
-			for (std::size_t j = 1; j < n; ++j)
-			{
-				*(++last) += x * std::cos(y * static_cast<T>(j));
-			}
-		}
-	}
-
 	/**
 	 * Evaluates a Chebyshev polynomial.
 	 *

src/math/quaternion-functions.hpp

@@ -1,510 +0,0 @@
-/*
- * Copyright (C) 2021  Christopher J. Howard
- *
- * This file is part of Antkeeper source code.
- *
- * Antkeeper source code is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Antkeeper source code is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef ANTKEEPER_MATH_QUATERNION_FUNCTIONS_HPP
-#define ANTKEEPER_MATH_QUATERNION_FUNCTIONS_HPP
-
-#include "math/constants.hpp"
-#include "math/matrix.hpp"
-#include "math/quaternion-type.hpp"
-#include "math/vector.hpp"
-#include <cmath>
-
-namespace math {
-
-/**
- * Adds two quaternions.
- *
- * @param x First quaternion.
- * @param y Second quaternion.
- * @return Sum of the two quaternions.
- */
-template <class T>
-quaternion<T> add(const quaternion<T>& x, const quaternion<T>& y);
-
-/**
- * Calculates the conjugate of a quaternion.
- *
- * @param x Quaternion from which the conjugate will be calculated.
- * @return Conjugate of the quaternion.
- */
-template <class T>
-quaternion<T> conjugate(const quaternion<T>& x);
-
-/**
- * Calculates the dot product of two quaternions.
- *
- * @param x First quaternion.
- * @param y Second quaternion.
- * @return Dot product of the two quaternions.
- */
-template <class T>
-T dot(const quaternion<T>& x, const quaternion<T>& y);
-
-/**
- * Divides a quaternion by a scalar.
- *
- * @param q Quaternion.
- * @param s Scalar.
- * @return Result of the division.
- */
-template <class T>
-quaternion<T> div(const quaternion<T>& q, T s);
-
-/**
- * Calculates the length of a quaternion.
- * 
- * @param x Quaternion to calculate the length of.
- * @return Length of the quaternion.
- */
-template <class T>
-T length(const quaternion<T>& x);
-
-/**
- * Calculates the squared length of a quaternion. The squared length can be calculated faster than the length because a call to `std::sqrt` is saved.
- * 
- * @param x Quaternion to calculate the squared length of.
- * @return Squared length of the quaternion.
- */
-template <class T>
-T length_squared(const quaternion<T>& x);
-
-/**
- * Performs linear interpolation between two quaternions.
- *
- * @param x First quaternion.
- * @param y Second quaternion.
- * @param a Interpolation factor.
- * @return Interpolated quaternion.
- */
-template <class T>
-quaternion<T> lerp(const quaternion<T>& x, const quaternion<T>& y, T a);
-
-/**
- * Creates a unit quaternion rotation using forward and up vectors.
- *
- * @param forward Unit forward vector.
- * @param up Unit up vector.
- * @return Unit rotation quaternion.
- */
-template <class T>
-quaternion<T> look_rotation(const vector<T, 3>& forward, vector<T, 3> up);
-
-/**
- * Converts a quaternion to a rotation matrix.
- *
- * @param q Unit quaternion.
- * @return Matrix representing the rotation described by `q`.
- */
-template <class T>
-matrix<T, 3, 3> matrix_cast(const quaternion<T>& q);
-
-/**
- * Multiplies two quaternions.
- *
- * @param x First quaternion.
- * @param y Second quaternion.
- * @return Product of the two quaternions.
- */
-template <class T>
-quaternion<T> mul(const quaternion<T>& x, const quaternion<T>& y);
-
-/**
- * Multiplies a quaternion by a scalar.
- *
- * @param q Quaternion.
- * @param s Scalar.
- * @return Product of the quaternion and scalar.
- */
-template <class T>
-quaternion<T> mul(const quaternion<T>& q, T s);
-
-/**
- * Rotates a 3-dimensional vector by a quaternion.
- *
- * @param q Unit quaternion.
- * @param v Vector.
- * @return Rotated vector.
- */
-template <class T>
-vector<T, 3> mul(const quaternion<T>& q, const vector<T, 3>& v);
-
-/**
- * Negates a quaternion.
- */
-template <class T>
-quaternion<T> negate(const quaternion<T>& x);
-
-/**
- * Performs normalized linear interpolation between two quaternions.
- *
- * @param x First quaternion.
- * @param y Second quaternion.
- * @param a Interpolation factor.
- * @return Interpolated quaternion.
- */
-template <class T>
-quaternion<T> nlerp(const quaternion<T>& x, const quaternion<T>& y, T a);
-
-/**
- * Normalizes a quaternion.
- */
-template <class T>
-quaternion<T> normalize(const quaternion<T>& x);
-
-/**
- * Creates a rotation from an angle and axis.
- *
- * @param angle Angle of rotation (in radians).
- * @param axis Axis of rotation
- * @return Quaternion representing the rotation.
- */
-template <class T>
-quaternion<T> angle_axis(T angle, const vector<T, 3>& axis);
-
-/**
- * Calculates the minimum rotation between two normalized direction vectors.
- *
- * @param source Normalized source direction.
- * @param destination Normalized destination direction.
- * @return Quaternion representing the minimum rotation between the source and destination vectors.
- */
-template <class T>
-quaternion<T> rotation(const vector<T, 3>& source, const vector<T, 3>& destination);
-
-/**
- * Performs spherical linear interpolation between two quaternions.
- *
- * @param x First quaternion.
- * @param y Second quaternion.
- * @param a Interpolation factor.
- * @return Interpolated quaternion.
- */
-template <class T>
-quaternion<T> slerp(const quaternion<T>& x, const quaternion<T>& y, T a);
-
-/**
- * Subtracts a quaternion from another quaternion.
- *
- * @param x First quaternion.
- * @param y Second quaternion.
- * @return Difference between the quaternions.
- */
-template <class T>
-quaternion<T> sub(const quaternion<T>& x, const quaternion<T>& y);
-
-/**
- * Decomposes a quaternion into swing and twist rotation components.
- *
- * @param[in] q Quaternion to decompose.
- * @param[in] a Axis of twist rotation.
- * @param[out] swing Swing rotation component.
- * @param[out] twist Twist rotation component.
- * @param[in] epsilon Threshold at which a number is considered zero.
- *
- * @see https://www.euclideanspace.com/maths/geometry/rotations/for/decomposition/
- */
-template <class T>
-void swing_twist(const quaternion<T>& q, const vector<T, 3>& a, quaternion<T>& qs, quaternion<T>& qt, T epsilon = T(1e-6));
-
-/**
- * Converts a 3x3 rotation matrix to a quaternion.
- *
- * @param m Rotation matrix.
- * @return Unit quaternion representing the rotation described by `m`.
- */
-template <class T>
-quaternion<T> quaternion_cast(const matrix<T, 3, 3>& m);
-
-/**
- * Types casts each quaternion component and returns a quaternion of the casted type.
- *
- * @tparam T2 Target quaternion component type.
- * @tparam T1 Source quaternion component type.
- * @param q Quaternion to type cast.
- * @return Type-casted quaternion.
- */
-template <class T2, class T1>
-quaternion<T2> type_cast(const quaternion<T1>& q);
-
-template <class T>
-inline quaternion<T> add(const quaternion<T>& x, const quaternion<T>& y)
-{
-	return {x.w + y.w, x.x + y.x, x.y + y.y, x.z + y.z};
-}
-
-template <class T>
-inline quaternion<T> conjugate(const quaternion<T>& x)
-{
-	return {x.w, -x.x, -x.y, -x.z};
-}
-
-template <class T>
-inline T dot(const quaternion<T>& x, const quaternion<T>& y)
-{
-	return {x.w * y.w + x.x * y.x + x.y * y.y + x.z * y.z};
-}
-
-template <class T>
-inline quaternion<T> div(const quaternion<T>& q, T s)
-{
-	return {q.w / s, q.x / s, q.y / s, q.z / s}
-}
-
-template <class T>
-inline T length(const quaternion<T>& x)
-{
-	return std::sqrt(x.w * x.w + x.x * x.x + x.y * x.y + x.z * x.z);
-}
-
-template <class T>
-inline T length_squared(const quaternion<T>& x)
-{
-	return x.w * x.w + x.x * x.x + x.y * x.y + x.z * x.z;
-}
-
-template <class T>
-inline quaternion<T> lerp(const quaternion<T>& x, const quaternion<T>& y, T a)
-{
-	return
-		{
-			(y.w - x.w) * a + x.w,
-			(y.x - x.x) * a + x.x,
-			(y.y - x.y) * a + x.y,
-			(y.z - x.z) * a + x.z
-		};
-}
-
-template <class T>
-quaternion<T> look_rotation(const vector<T, 3>& forward, vector<T, 3> up)
-{
-	vector<T, 3> right = normalize(cross(forward, up));
-	up = cross(right, forward);
-
-	matrix<T, 3, 3> m =
-		{
-			right,
-			up,
-			-forward
-		};
-
-	// Convert to quaternion
-	return normalize(quaternion_cast(m));
-}
-
-template <class T>
-matrix<T, 3, 3> matrix_cast(const quaternion<T>& q)
-{
-	const T xx = q.x * q.x;
-	const T xy = q.x * q.y;
-	const T xz = q.x * q.z;
-	const T xw = q.x * q.w;
-	const T yy = q.y * q.y;
-	const T yz = q.y * q.z;
-	const T yw = q.y * q.w;
-	const T zz = q.z * q.z;
-	const T zw = q.z * q.w;
-
-	return
-		{
-			T(1) - (yy + zz) * T(2), (xy + zw) * T(2), (xz - yw) * T(2),
-			(xy - zw) * T(2), T(1) - (xx + zz) * T(2), (yz + xw) * T(2),
-			(xz + yw) * T(2), (yz - xw) * T(2), T(1) - (xx + yy) * T(2)
-		};
-}
-
-template <class T>
-quaternion<T> mul(const quaternion<T>& x, const quaternion<T>& y)
-{
-	return
-		{
-			-x.x * y.x - x.y * y.y - x.z * y.z + x.w * y.w,
-			 x.x * y.w + x.y * y.z - x.z * y.y + x.w * y.x,
-			-x.x * y.z + x.y * y.w + x.z * y.x + x.w * y.y,
-			 x.x * y.y - x.y * y.x + x.z * y.w + x.w * y.z
-		};
-}
-
-template <class T>
-inline quaternion<T> mul(const quaternion<T>& q, T s)
-{
-	return {q.w * s, q.x * s, q.y * s, q.z * s};
-}
-
-template <class T>
-vector<T, 3> mul(const quaternion<T>& q, const vector<T, 3>& v)
-{
-	const vector<T, 3>& i = reinterpret_cast<const vector<T, 3>&>(q.x);	
-	return add(add(mul(i, dot(i, v) * T(2)), mul(v, (q.w * q.w - dot(i, i)))), mul(cross(i, v), q.w * T(2)));
-}
-
-template <class T>
-inline quaternion<T> negate(const quaternion<T>& x)
-{
-	return {-x.w, -x.x, -x.y, -x.z};
-}
-
-template <class T>
-quaternion<T> nlerp(const quaternion<T>& x, const quaternion<T>& y, T a)
-{
-	return normalize(add(mul(x, T(1) - a), mul(y, a * std::copysign(T(1), dot(x, y)))));
-}
-
-template <class T>
-inline quaternion<T> normalize(const quaternion<T>& x)
-{
-	return mul(x, T(1) / length(x));
-}
-
-template <class T>
-quaternion<T> angle_axis(T angle, const vector<T, 3>& axis)
-{
-	T s = std::sin(angle * T(0.5));
-	return {static_cast<T>(std::cos(angle * T(0.5))), axis.x() * s, axis.y() * s, axis.z() * s};
-}
-
-template <class T>
-quaternion<T> rotation(const vector<T, 3>& source, const vector<T, 3>& destination)
-{
-	quaternion<T> q;
-	q.w = dot(source, destination);
-	reinterpret_cast<vector<T, 3>&>(q.x) = cross(source, destination);
-
-	q.w += length(q);
-	return normalize(q);
-}
-
-template <class T>
-quaternion<T> slerp(const quaternion<T>& x, const quaternion<T>& y, T a)
-{
-	T cos_theta = dot(x, y);
-
-	constexpr T epsilon = T(0.0005);
-	if (cos_theta > T(1) - epsilon)
-	{
-		return normalize(lerp(x, y, a));
-	}
-
-	cos_theta = std::max<T>(T(-1), std::min<T>(T(1), cos_theta));
-	T theta = static_cast<T>(std::acos(cos_theta)) * a;
-
-	quaternion<T> z = normalize(sub(y, mul(x, cos_theta)));
-
-	return add(mul(x, static_cast<T>(std::cos(theta))), mul(z, static_cast<T>(std::sin(theta))));
-}
-
-template <class T>
-inline quaternion<T> sub(const quaternion<T>& x, const quaternion<T>& y)
-{
-	return {x.w - y.w, x.x - y.x, x.y - y.y, x.z - y.z};
-}
-
-template <class T>
-void swing_twist(const quaternion<T>& q, const vector<T, 3>& a, quaternion<T>& qs, quaternion<T>& qt, T epsilon)
-{
-	if (q.x * q.x + q.y * q.y + q.z * q.z > epsilon)
-	{
-		const vector<T, 3> pa = mul(a, (a.x() * q.x + a.y() * q.y + a.z() * q.z));
-		qt = normalize(quaternion<T>{q.w, pa.x(), pa.y(), pa.z()});
-		qs = mul(q, conjugate(qt));
-	}
-	else
-	{
-		qt = angle_axis(pi<T>, a);
-		
-		const vector<T, 3> qa = mul(q, a);
-		const vector<T, 3> sa = cross(a, qa);
-		if (length_squared(sa) > epsilon)
-			qs = angle_axis(std::acos(dot(a, qa)), sa);
-		else
-			qs = quaternion<T>::identity;
-	}
-}
-
-template <class T>
-quaternion<T> quaternion_cast(const matrix<T, 3, 3>& m)
-{
-	T trace = m[0][0] + m[1][1] + m[2][2];
-	
-	if (trace > T(0))
-	{
-		T s = T(0.5) / std::sqrt(trace + T(1));
-		return
-		{
-			T(0.25) / s,
-			(m[1][2] - m[2][1]) * s,
-			(m[2][0] - m[0][2]) * s,
-			(m[0][1] - m[1][0]) * s
-		};
-	}
-	else
-	{
-		if (m[0][0] > m[1][1] && m[0][0] > m[2][2])
-		{
-			T s = T(2) * std::sqrt(T(1) + m[0][0] - m[1][1] - m[2][2]);
-			
-			return
-			{
-				(m[1][2] - m[2][1]) / s,
-				T(0.25) * s,
-				(m[1][0] + m[0][1]) / s,
-				(m[2][0] + m[0][2]) / s
-			};
-		}
-		else if (m[1][1] > m[2][2])
-		{
-			T s = T(2) * std::sqrt(T(1) + m[1][1] - m[0][0] - m[2][2]);
-			return
-			{
-				(m[2][0] - m[0][2]) / s,
-				(m[1][0] + m[0][1]) / s,
-				T(0.25) * s,
-				(m[2][1] + m[1][2]) / s
-			};
-		}
-		else
-		{
-			T s = T(2) * std::sqrt(T(1) + m[2][2] - m[0][0] - m[1][1]);
-			return
-			{
-				(m[0][1] - m[1][0]) / s,
-				(m[2][0] + m[0][2]) / s,
-				(m[2][1] + m[1][2]) / s,
-				T(0.25) * s
-			};
-		}
-	}
-}
-
-template <class T2, class T1>
-inline quaternion<T2> type_cast(const quaternion<T1>& q)
-{
-	return quaternion<T2>
-	{
-		static_cast<T2>(q.w),
-		static_cast<T2>(q.x),
-		static_cast<T2>(q.y),
-		static_cast<T2>(q.z)
-	};
-}
-
-} // namespace math
-
-#endif // ANTKEEPER_MATH_QUATERNION_FUNCTIONS_HPP
-

src/math/quaternion-operators.hpp

@@ -1,97 +0,0 @@
-/*
- * Copyright (C) 2021  Christopher J. Howard
- *
- * This file is part of Antkeeper source code.
- *
- * Antkeeper source code is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Antkeeper source code is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef ANTKEEPER_MATH_QUATERNION_OPERATORS_HPP
-#define ANTKEEPER_MATH_QUATERNION_OPERATORS_HPP
-
-#include "math/quaternion-type.hpp"
-#include "math/quaternion-functions.hpp"
-
-/// @copydoc math::add(const math::quaternion<T>&, const math::quaternion<T>&)
-template <class T>
-math::quaternion<T> operator+(const math::quaternion<T>& x, const math::quaternion<T>& y);
-
-/// @copydoc math::div(const math::quaternion<T>&, T)
-template <class T>
-math::quaternion<T> operator/(const math::quaternion<T>& q, T s);
-
-/// @copydoc math::mul(const math::quaternion<T>&, const math::quaternion<T>&)
-template <class T>
-math::quaternion<T> operator*(const math::quaternion<T>& x, const math::quaternion<T>& y);
-
-/// @copydoc math::mul(const math::quaternion<T>&, T)
-template <class T>
-math::quaternion<T> operator*(const math::quaternion<T>& q, T s);
-
-/// @copydoc math::mul(const math::quaternion<T>&, const math::vector<T, 3>&)
-template <class T>
-math::vector<T, 3> operator*(const math::quaternion<T>& q, const math::vector<T, 3>& v);
-
-/// @copydoc math::sub(const math::quaternion<T>&, const math::quaternion<T>&)
-template <class T>
-math::quaternion<T> operator-(const math::quaternion<T>& x, const math::quaternion<T>& y);
-
-/// @copydoc math::negate(const math::quaternion<T>&)
-template <class T, std::size_t N>
-math::quaternion<T> operator-(const math::quaternion<T>& x);
-
-template <class T>
-inline math::quaternion<T> operator+(const math::quaternion<T>& x, const math::quaternion<T>& y)
-{
-	return add(x, y);
-}
-
-template <class T>
-inline math::quaternion<T> operator/(const math::quaternion<T>& q, T s)
-{
-	return div(q, s);
-}
-
-template <class T>
-inline math::quaternion<T> operator*(const math::quaternion<T>& x, const math::quaternion<T>& y)
-{
-	return mul(x, y);
-}
-
-template <class T>
-inline math::quaternion<T> operator*(const math::quaternion<T>& q, T s)
-{
-	return mul(q, s);
-}
-
-template <class T>
-inline math::vector<T, 3> operator*(const math::quaternion<T>& q, const math::vector<T, 3>& v)
-{
-	return mul(q, v);
-}
-
-template <class T>
-inline math::quaternion<T> operator-(const math::quaternion<T>& x, const math::quaternion<T>& y)
-{
-	return sub(x, y);
-}
-
-template <class T, std::size_t N>
-inline math::quaternion<T> operator-(const math::quaternion<T>& x)
-{
-	return negate(x);
-}
-
-#endif // ANTKEEPER_MATH_QUATERNION_OPERATORS_HPP
-

src/math/quaternion-type.hpp

@@ -1,92 +0,0 @@
-/*
- * Copyright (C) 2021  Christopher J. Howard
- *
- * This file is part of Antkeeper source code.
- *
- * Antkeeper source code is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Antkeeper source code is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef ANTKEEPER_MATH_QUATERNION_TYPE_HPP
-#define ANTKEEPER_MATH_QUATERNION_TYPE_HPP
-
-#include <cmath>
-
-namespace math {
-
-/**
- * A quaternion type is a tuple made of a scalar (real) part and vector (imaginary) part.
- *
- * @tparam T Scalar type.
- */
-template <class T>
-struct quaternion
-{
-	typedef T scalar_type;
-	scalar_type w;
-	scalar_type x;
-	scalar_type y;
-	scalar_type z;
-	
-	/// Rotation identity quaternion.
-	static const quaternion identity;
-	
-	/**
-	 * Returns a quaternion representing a rotation about the x-axis.
-	 *
-	 * @param angle Angle of rotation, in radians.
-	 * @return Quaternion representing an x-axis rotation.
-	 */
-	static quaternion rotate_x(scalar_type angle);
-	
-	/**
-	 * Returns a quaternion representing a rotation about the y-axis.
-	 *
-	 * @param angle Angle of rotation, in radians.
-	 * @return Quaternion representing an y-axis rotation.
-	 */
-	static quaternion rotate_y(scalar_type angle);
-	
-	/**
-	 * Returns a quaternion representing a rotation about the z-axis.
-	 *
-	 * @param angle Angle of rotation, in radians.
-	 * @return Quaternion representing an z-axis rotation.
-	 */
-	static quaternion rotate_z(scalar_type angle);
-};
-
-template <class T>
-constexpr quaternion<T> quaternion<T>::identity = {T(1), T(0), T(0), T(0)};
-
-template <class T>
-quaternion<T> quaternion<T>::rotate_x(scalar_type angle)
-{
-	return {std::cos(angle * T(0.5)), std::sin(angle * T(0.5)), T(0), T(0)};
-}
-
-template <class T>
-quaternion<T> quaternion<T>::rotate_y(scalar_type angle)
-{
-	return {std::cos(angle * T(0.5)), T(0), std::sin(angle * T(0.5)), T(0)};
-}
-
-template <class T>
-quaternion<T> quaternion<T>::rotate_z(scalar_type angle)
-{
-	return {std::cos(angle * T(0.5)), T(0), T(0), std::sin(angle * T(0.5))};
-}
-
-} // namespace math
-
-#endif // ANTKEEPER_MATH_QUATERNION_TYPE_HPP

src/math/quaternion.hpp

@@ -0,0 +1,959 @@
+/*
+ * Copyright (C) 2021  Christopher J. Howard
+ *
+ * This file is part of Antkeeper source code.
+ *
+ * Antkeeper source code is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Antkeeper source code is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef ANTKEEPER_MATH_QUATERNION_HPP
+#define ANTKEEPER_MATH_QUATERNION_HPP
+
+#include "math/constants.hpp"
+#include "math/matrix.hpp"
+#include "math/vector.hpp"
+#include <cmath>
+#include <istream>
+#include <ostream>
+
+namespace math {
+
+/**
+ * A quaternion type is a tuple made of a scalar (real) part and vector (imaginary) part.
+ *
+ * @tparam T Scalar type.
+ */
+template <class T>
+struct quaternion
+{
+	/// Scalar type.
+	typedef T scalar_type;
+	
+	/// Vector type.
+	typedef vector<T, 3> vector_type;
+	
+	/// Rotation matrix type.
+	typedef matrix<T, 3, 3> matrix_type;
+	
+	/// Quaternion real part.
+	scalar_type r;
+	
+	/// Quaternion imaginary part.
+	vector_type i;
+	
+	/// Returns a reference to the quaternion real part.
+	/// @{
+	constexpr inline scalar_type& w() noexcept
+	{
+		return r;
+	}
+	constexpr inline const scalar_type& w() const noexcept
+	{
+		return r;
+	}
+	/// @}
+	
+	/// Returns a reference to the first element of the quaternion imaginary part.
+	/// @{
+	constexpr inline scalar_type& x() noexcept
+	{
+		return i.x();
+	}
+	constexpr inline const scalar_type& x() const noexcept
+	{
+		return i.x();
+	}
+	/// @}
+	
+	/// Returns a reference to the second element of the quaternion imaginary part.
+	/// @{
+	constexpr inline scalar_type& y() noexcept
+	{
+		return i.y();
+	}
+	constexpr inline const scalar_type& y() const noexcept
+	{
+		return i.y();
+	}
+	/// @}
+	
+	/// Returns a reference to the third element of the quaternion imaginary part.
+	/// @{
+	constexpr inline scalar_type& z() noexcept
+	{
+		return i.z();
+	}
+	constexpr inline const scalar_type& z() const noexcept
+	{
+		return i.z();
+	}
+	/// @}
+	
+	/**
+	 * Returns a quaternion representing a rotation about the x-axis.
+	 *
+	 * @param angle Angle of rotation, in radians.
+	 *
+	 * @return Quaternion representing an x-axis rotation.
+	 */
+	static quaternion rotate_x(scalar_type angle)
+	{
+		return {std::cos(angle * T(0.5)), std::sin(angle * T(0.5)), T(0), T(0)};
+	}
+	
+	/**
+	 * Returns a quaternion representing a rotation about the y-axis.
+	 *
+	 * @param angle Angle of rotation, in radians.
+	 *
+	 * @return Quaternion representing an y-axis rotation.
+	 */
+	static quaternion rotate_y(scalar_type angle)
+	{
+		return {std::cos(angle * T(0.5)), T(0), std::sin(angle * T(0.5)), T(0)};
+	}
+	
+	/**
+	 * Returns a quaternion representing a rotation about the z-axis.
+	 *
+	 * @param angle Angle of rotation, in radians.
+	 * @return Quaternion representing an z-axis rotation.
+	 */
+	static quaternion rotate_z(scalar_type angle)
+	{
+		return {std::cos(angle * T(0.5)), T(0), T(0), std::sin(angle * T(0.5))};
+	}
+	
+	/**
+	 * Type-casts the quaternion scalars using `static_cast`.
+	 *
+	 * @tparam U Target scalar type.
+	 *
+	 * @return Type-casted quaternion.
+	 */
+	template <class U>
+	constexpr inline explicit operator quaternion<U>() const noexcept
+	{
+		return {static_cast<U>(r), vector<U, 3>(i)};
+	}
+	
+	/**
+	 * Casts the quaternion to a 4-element vector, with the real part as the first element and the imaginary part as the following three elements.
+	 *
+	 * @return Vector containing the real and imaginary parts of the quaternion.
+	 */
+	constexpr inline explicit operator vector<T, 4>() const noexcept
+	{
+		return {r, i[0], i[1], i[2]};
+	}
+	
+	/// Returns a zero quaternion, where every scalar is equal to zero.
+	static constexpr quaternion zero() noexcept
+	{
+		return {};
+	}
+	
+	/// Returns a rotation identity quaternion.
+	static constexpr quaternion identity() noexcept
+	{
+		return {T{1}, vector_type::zero()};
+	}
+};
+
+/**
+ * Adds two quaternions.
+ *
+ * @param a First quaternion.
+ * @param b Second quaternion.
+ *
+ * @return Sum of the two quaternions.
+ */
+template <class T>
+constexpr quaternion<T> add(const quaternion<T>& a, const quaternion<T>& b) noexcept;
+
+/**
+ * Adds a quaternion and a scalar.
+ *
+ * @param a First value.
+ * @param b Second second value.
+ *
+ * @return Sum of the quaternion and scalar.
+ */
+template <class T>
+constexpr quaternion<T> add(const quaternion<T>& a, T b) noexcept;
+
+/**
+ * Calculates the conjugate of a quaternion.
+ *
+ * @param q Quaternion from which the conjugate will be calculated.
+ *
+ * @return Conjugate of the quaternion.
+ */
+template <class T>
+constexpr quaternion<T> conjugate(const quaternion<T>& q) noexcept;
+
+/**
+ * Calculates the dot product of two quaternions.
+ *
+ * @param a First quaternion.
+ * @param b Second quaternion.
+ *
+ * @return Dot product of the two quaternions.
+ */
+template <class T>
+constexpr T dot(const quaternion<T>& a, const quaternion<T>& b) noexcept;
+
+/**
+ * Divides a quaternion by another quaternion.
+ *
+ * @param a First value.
+ * @param b Second value.
+ *
+ * @return Result of the division.
+ */
+template <class T>
+constexpr quaternion<T> div(const quaternion<T>& a, const quaternion<T>& b) noexcept;
+
+/**
+ * Divides a quaternion by a scalar.
+ *
+ * @param a Quaternion.
+ * @param b Scalar.
+ *
+ * @return Result of the division.
+ */
+template <class T>
+constexpr quaternion<T> div(const quaternion<T>& a, T b) noexcept;
+
+/**
+ * Divides a scalar by a quaternion.
+ *
+ * @param a Scalar.
+ * @param b Quaternion.
+ *
+ * @return Result of the division.
+ */
+template <class T>
+constexpr quaternion<T> div(T a, const quaternion<T>& b) noexcept;
+
+/**
+ * Calculates the length of a quaternion.
+ * 
+ * @param q Quaternion to calculate the length of.
+ *
+ * @return Length of the quaternion.
+ */
+template <class T>
+T length(const quaternion<T>& q);
+
+/**
+ * Calculates the squared length of a quaternion. The squared length can be calculated faster than the length because a call to `std::sqrt` is saved.
+ * 
+ * @param q Quaternion to calculate the squared length of.
+ *
+ * @return Squared length of the quaternion.
+ */
+template <class T>
+constexpr T length_squared(const quaternion<T>& q) noexcept;
+
+/**
+ * Performs linear interpolation between two quaternions.
+ *
+ * @param a First quaternion.
+ * @param b Second quaternion.
+ * @param t Interpolation factor.
+ *
+ * @return Interpolated quaternion.
+ */
+template <class T>
+constexpr quaternion<T> lerp(const quaternion<T>& a, const quaternion<T>& b, T t) noexcept;
+
+/**
+ * Creates a unit quaternion rotation using forward and up vectors.
+ *
+ * @param forward Unit forward vector.
+ * @param up Unit up vector.
+ *
+ * @return Unit rotation quaternion.
+ */
+template <class T>
+quaternion<T> look_rotation(const vector<T, 3>& forward, vector<T, 3> up);
+
+/**
+ * Converts a quaternion to a rotation matrix.
+ *
+ * @param q Unit quaternion.
+ *
+ * @return Matrix representing the rotation described by `q`.
+ */
+template <class T>
+constexpr matrix<T, 3, 3> matrix_cast(const quaternion<T>& q) noexcept;
+
+/**
+ * Multiplies two quaternions.
+ *
+ * @param a First quaternion.
+ * @param b Second quaternion.
+ *
+ * @return Product of the two quaternions.
+ */
+template <class T>
+constexpr quaternion<T> mul(const quaternion<T>& a, const quaternion<T>& b) noexcept;
+
+/**
+ * Multiplies a quaternion by a scalar.
+ *
+ * @param a First value.
+ * @param b Second value.
+ *
+ * @return Product of the quaternion and scalar.
+ */
+template <class T>
+constexpr quaternion<T> mul(const quaternion<T>& a, T b) noexcept;
+
+/**
+ * Calculates the product of a quaternion and a vector.
+ *
+ * @param a First value.
+ * @param b second value.
+ *
+ * @return Product of the quaternion and vector.
+ */
+/// @{
+template <class T>
+constexpr vector<T, 3> mul(const quaternion<T>& a, const vector<T, 3>& b) noexcept;
+template <class T>
+constexpr vector<T, 3> mul(const vector<T, 3>& a, const quaternion<T>& b) noexcept;
+/// @}
+
+/**
+ * Negates a quaternion.
+ *
+ * @param q Quaternion to negate.
+ *
+ * @return Negated quaternion.
+ */
+template <class T>
+constexpr quaternion<T> negate(const quaternion<T>& q) noexcept;
+
+/**
+ * Performs normalized linear interpolation between two quaternions.
+ *
+ * @param a First quaternion.
+ * @param b Second quaternion.
+ * @param t Interpolation factor.
+ *
+ * @return Interpolated quaternion.
+ */
+template <class T>
+quaternion<T> nlerp(const quaternion<T>& a, const quaternion<T>& b, T t);
+
+/**
+ * Normalizes a quaternion.
+ *
+ * @param q Quaternion to normalize.
+ *
+ * @return Normalized quaternion.
+ */
+template <class T>
+quaternion<T> normalize(const quaternion<T>& q);
+
+/**
+ * Creates a rotation from an angle and axis.
+ *
+ * @param angle Angle of rotation (in radians).
+ * @param axis Axis of rotation
+ *
+ * @return Quaternion representing the rotation.
+ */
+template <class T>
+quaternion<T> angle_axis(T angle, const vector<T, 3>& axis);
+
+/**
+ * Calculates the minimum rotation between two normalized direction vectors.
+ *
+ * @param source Normalized source direction.
+ * @param destination Normalized destination direction.
+ *
+ * @return Quaternion representing the minimum rotation between the source and destination vectors.
+ */
+template <class T>
+quaternion<T> rotation(const vector<T, 3>& source, const vector<T, 3>& destination);
+
+/**
+ * Performs spherical linear interpolation between two quaternions.
+ *
+ * @param a First quaternion.
+ * @param b Second quaternion.
+ * @param t Interpolation factor.
+ *
+ * @return Interpolated quaternion.
+ */
+template <class T>
+quaternion<T> slerp(const quaternion<T>& a, const quaternion<T>& b, T t, T error = T{1e-6});
+
+/**
+ * Subtracts a quaternion from another quaternion.
+ *
+ * @param a First quaternion.
+ * @param b Second quaternion.
+ *
+ * @return Difference between the quaternions.
+ */
+template <class T>
+constexpr quaternion<T> sub(const quaternion<T>& a, const quaternion<T>& b) noexcept;
+
+/**
+ * Subtracts a quaternion and a scalar.
+ *
+ * @param a First value.
+ * @param b Second second.
+ *
+ * @return Difference between the quaternion and scalar.
+ */
+/// @{
+template <class T>
+constexpr quaternion<T> sub(const quaternion<T>& a, T b) noexcept;
+template <class T>
+constexpr quaternion<T> sub(T a, const quaternion<T>& b) noexcept;
+/// @}
+
+/**
+ * Decomposes a quaternion into swing and twist rotation components.
+ *
+ * @param[in] q Quaternion to decompose.
+ * @param[in] a Axis of twist rotation.
+ * @param[out] swing Swing rotation component.
+ * @param[out] twist Twist rotation component.
+ * @param[in] error Threshold at which a number is considered zero.
+ *
+ * @see https://www.euclideanspace.com/maths/geometry/rotations/for/decomposition/
+ */
+template <class T>
+void swing_twist(const quaternion<T>& q, const vector<T, 3>& a, quaternion<T>& qs, quaternion<T>& qt, T error = T{1e-6});
+
+/**
+ * Converts a 3x3 rotation matrix to a quaternion.
+ *
+ * @param m Rotation matrix.
+ *
+ * @return Unit quaternion representing the rotation described by @p m.
+ */
+template <class T>
+quaternion<T> quaternion_cast(const matrix<T, 3, 3>& m);
+
+template <class T>
+constexpr inline quaternion<T> add(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return {a.r + b.r, a.i + b.i};
+}
+
+template <class T>
+constexpr inline quaternion<T> add(const quaternion<T>& a, T b) noexcept
+{
+	return {a.r + b, a.i + b};
+}
+
+template <class T>
+constexpr inline quaternion<T> conjugate(const quaternion<T>& q) noexcept
+{
+	return {q.r, -q.i};
+}
+
+template <class T>
+constexpr inline T dot(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return a.r * b.r + dot(a.i, b.i);
+}
+
+template <class T>
+constexpr inline quaternion<T> div(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return {a.r / b.r, a.i / b.i};
+}
+
+template <class T>
+constexpr inline quaternion<T> div(const quaternion<T>& a, T b) noexcept
+{
+	return {a.r / b, a.i / b};
+}
+
+template <class T>
+constexpr inline quaternion<T> div(T a, const quaternion<T>& b) noexcept
+{
+	return {a / b.r, a / b.i};
+}
+
+template <class T>
+inline T length(const quaternion<T>& q)
+{
+	return std::sqrt(length_squared(q));
+}
+
+template <class T>
+constexpr inline T length_squared(const quaternion<T>& q) noexcept
+{
+	return q.r * q.r + length_squared(q.i);
+}
+
+template <class T>
+constexpr inline quaternion<T> lerp(const quaternion<T>& a, const quaternion<T>& b, T t) noexcept
+{
+	return
+	{
+		(b.r - a.r) * t + a.r,
+		(b.i - a.i) * t + a.i
+	};
+}
+
+template <class T>
+quaternion<T> look_rotation(const vector<T, 3>& forward, vector<T, 3> up)
+{
+	vector<T, 3> right = normalize(cross(forward, up));
+	up = cross(right, forward);
+
+	matrix<T, 3, 3> m =
+	{
+		right,
+		up,
+		-forward
+	};
+
+	// Convert to quaternion
+	return normalize(quaternion_cast(m));
+}
+
+template <class T>
+constexpr matrix<T, 3, 3> matrix_cast(const quaternion<T>& q) noexcept
+{
+	const T xx = q.x() * q.x();
+	const T xy = q.x() * q.y();
+	const T xz = q.x() * q.z();
+	const T xw = q.x() * q.w();
+	const T yy = q.y() * q.y();
+	const T yz = q.y() * q.z();
+	const T yw = q.y() * q.w();
+	const T zz = q.z() * q.z();
+	const T zw = q.z() * q.w();
+
+	return
+		{
+			T(1) - (yy + zz) * T(2), (xy + zw) * T(2), (xz - yw) * T(2),
+			(xy - zw) * T(2), T(1) - (xx + zz) * T(2), (yz + xw) * T(2),
+			(xz + yw) * T(2), (yz - xw) * T(2), T(1) - (xx + yy) * T(2)
+		};
+}
+
+template <class T>
+constexpr quaternion<T> mul(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return
+		{
+			-a.x() * b.x() - a.y() * b.y() - a.z() * b.z() + a.w() * b.w(),
+			 a.x() * b.w() + a.y() * b.z() - a.z() * b.y() + a.w() * b.x(),
+			-a.x() * b.z() + a.y() * b.w() + a.z() * b.x() + a.w() * b.y(),
+			 a.x() * b.y() - a.y() * b.x() + a.z() * b.w() + a.w() * b.z()
+		};
+}
+
+template <class T>
+constexpr inline quaternion<T> mul(const quaternion<T>& a, T b) noexcept
+{
+	return {a.r * b, a.i * b};
+}
+
+template <class T>
+constexpr vector<T, 3> mul(const quaternion<T>& a, const vector<T, 3>& b) noexcept
+{
+	return a.i * dot(a.i, b) * T(2) + b * (a.r * a.r - length_squared(a.i)) + cross(a.i, b) * a.r * T(2);
+}
+
+template <class T>
+constexpr inline vector<T, 3> mul(const vector<T, 3>& a, const quaternion<T>& b) noexcept
+{
+	return mul(conjugate(b), a);
+}
+
+template <class T>
+constexpr inline quaternion<T> negate(const quaternion<T>& q) noexcept
+{
+	return {-q.r, -q.i};
+}
+
+template <class T>
+quaternion<T> nlerp(const quaternion<T>& a, const quaternion<T>& b, T t)
+{
+	return normalize(add(mul(a, T(1) - t), mul(b, t * std::copysign(T(1), dot(a, b)))));
+}
+
+template <class T>
+inline quaternion<T> normalize(const quaternion<T>& q)
+{
+	return mul(q, T(1) / length(q));
+}
+
+template <class T>
+quaternion<T> angle_axis(T angle, const vector<T, 3>& axis)
+{
+	angle *= T{0.5};
+	return {std::cos(angle), axis * std::sin(angle)};
+}
+
+template <class T>
+quaternion<T> rotation(const vector<T, 3>& source, const vector<T, 3>& destination)
+{
+	quaternion<T> q = {dot(source, destination), cross(source, destination)};
+	q.w() += length(q);
+	return normalize(q);
+}
+
+template <class T>
+quaternion<T> slerp(const quaternion<T>& a, const quaternion<T>& b, T t, T error)
+{
+	T cos_theta = dot(a, b);
+
+	if (cos_theta > T(1) - error)
+		return normalize(lerp(a, b, t));
+
+	cos_theta = std::max<T>(T(-1), std::min<T>(T(1), cos_theta));
+	const T theta = std::acos(cos_theta) * t;
+	
+	quaternion<T> c = normalize(sub(b, mul(a, cos_theta)));
+	
+	return add(mul(a, std::cos(theta)), mul(c, std::sin(theta)));
+}
+
+template <class T>
+constexpr inline quaternion<T> sub(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return {a.r - b.r, a.i - b.i};
+}
+
+template <class T>
+constexpr inline quaternion<T> sub(const quaternion<T>& a, T b) noexcept
+{
+	return {a.r - b, a.i - b};
+}
+
+template <class T>
+constexpr inline quaternion<T> sub(T a, const quaternion<T>& b) noexcept
+{
+	return {a - b.r, a - b.i};
+}
+
+template <class T>
+void swing_twist(const quaternion<T>& q, const vector<T, 3>& a, quaternion<T>& qs, quaternion<T>& qt, T error)
+{
+	if (length_squared(q.i) > error)
+	{
+		qt = normalize(quaternion<T>{q.w(), a * dot(a, q.i)});
+		qs = mul(q, conjugate(qt));
+	}
+	else
+	{
+		qt = angle_axis(pi<T>, a);
+		
+		const vector<T, 3> qa = mul(q, a);
+		const vector<T, 3> sa = cross(a, qa);
+		if (length_squared(sa) > error)
+			qs = angle_axis(std::acos(dot(a, qa)), sa);
+		else
+			qs = quaternion<T>::identity();
+	}
+}
+
+template <class T>
+quaternion<T> quaternion_cast(const matrix<T, 3, 3>& m)
+{
+	const T t = trace(m);
+	
+	if (t > T(0))
+	{
+		T s = T(0.5) / std::sqrt(t + T(1));
+		return
+		{
+			T(0.25) / s,
+			(m[1][2] - m[2][1]) * s,
+			(m[2][0] - m[0][2]) * s,
+			(m[0][1] - m[1][0]) * s
+		};
+	}
+	else
+	{
+		if (m[0][0] > m[1][1] && m[0][0] > m[2][2])
+		{
+			T s = T(2) * std::sqrt(T(1) + m[0][0] - m[1][1] - m[2][2]);
+			
+			return
+			{
+				(m[1][2] - m[2][1]) / s,
+				T(0.25) * s,
+				(m[1][0] + m[0][1]) / s,
+				(m[2][0] + m[0][2]) / s
+			};
+		}
+		else if (m[1][1] > m[2][2])
+		{
+			T s = T(2) * std::sqrt(T(1) + m[1][1] - m[0][0] - m[2][2]);
+			return
+			{
+				(m[2][0] - m[0][2]) / s,
+				(m[1][0] + m[0][1]) / s,
+				T(0.25) * s,
+				(m[2][1] + m[1][2]) / s
+			};
+		}
+		else
+		{
+			T s = T(2) * std::sqrt(T(1) + m[2][2] - m[0][0] - m[1][1]);
+			return
+			{
+				(m[0][1] - m[1][0]) / s,
+				(m[2][0] + m[0][2]) / s,
+				(m[2][1] + m[1][2]) / s,
+				T(0.25) * s
+			};
+		}
+	}
+}
+
+namespace operators {
+
+/// @copydoc add(const quaternion<T>&, const quaternion<T>&)
+template <class T>
+constexpr inline quaternion<T> operator+(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return add(a, b);
+}
+
+/// @copydoc add(const quaternion<T>&, T)
+/// @{
+template <class T>
+constexpr inline quaternion<T> operator+(const quaternion<T>& a, T b) noexcept
+{
+	return add(a, b);
+}
+template <class T>
+constexpr inline quaternion<T> operator+(T a, const quaternion<T>& b) noexcept
+{
+	return add(b, a);
+}
+/// @}
+
+/// @copydoc div(const quaternion<T>&, const quaternion<T>&)
+template <class T>
+constexpr inline quaternion<T> operator/(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return div(a, b);
+}
+
+/// @copydoc div(const quaternion<T>&, T)
+template <class T>
+constexpr inline quaternion<T> operator/(const quaternion<T>& a, T b) noexcept
+{
+	return div(a, b);
+}
+
+/// @copydoc div(T, const quaternion<T>&)
+template <class T>
+constexpr inline quaternion<T> operator/(T a, const quaternion<T>& b) noexcept
+{
+	return div(a, b);
+}
+
+/// @copydoc mul(const quaternion<T>&, const quaternion<T>&)
+template <class T>
+constexpr inline quaternion<T> operator*(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return mul(a, b);
+}
+
+/// @copydoc mul(const quaternion<T>&, T)
+/// @{
+template <class T>
+constexpr inline quaternion<T> operator*(const quaternion<T>& a, T b) noexcept
+{
+	return mul(a, b);
+}
+template <class T>
+constexpr inline quaternion<T> operator*(T a, const quaternion<T>& b) noexcept
+{
+	return mul(b, a);
+}
+/// @}
+
+/// @copydoc mul(const quaternion<T>&, const vector<T, 3>&)
+template <class T>
+constexpr inline vector<T, 3> operator*(const quaternion<T>& a, const vector<T, 3>& b) noexcept
+{
+	return mul(a, b);
+}
+
+/// @copydoc mul(const vector<T, 3>&, const quaternion<T>&)
+template <class T>
+constexpr inline vector<T, 3> operator*(const vector<T, 3>& a, const quaternion<T>& b) noexcept
+{
+	return mul(a, b);
+}
+
+/// @copydoc sub(const quaternion<T>&, const quaternion<T>&)
+template <class T>
+constexpr inline quaternion<T> operator-(const quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return sub(a, b);
+}
+
+/// @copydoc sub(const quaternion<T>&, T)
+/// @{
+template <class T>
+constexpr inline quaternion<T> operator-(const quaternion<T>& a, T b) noexcept
+{
+	return sub(a, b);
+}
+template <class T>
+constexpr inline quaternion<T> operator-(T a, const quaternion<T>& b) noexcept
+{
+	return sub(a, b);
+}
+/// @}
+
+/// @copydoc negate(const quaternion<T>&)
+template <class T>
+constexpr inline quaternion<T> operator-(const quaternion<T>& q) noexcept
+{
+	return negate(q);
+}
+
+/**
+ * Adds two values and stores the result in the first value.
+ *
+ * @param a First value.
+ * @param b Second value.
+ *
+ * @return Reference to the first value.
+ */
+/// @{
+template <class T>
+constexpr inline quaternion<T>& operator+=(quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return (a = a + b);
+}
+template <class T>
+constexpr inline quaternion<T>& operator+=(quaternion<T>& a, T b) noexcept
+{
+	return (a = a + b);
+}
+/// @}
+
+/**
+ * Subtracts the first value by the second value and stores the result in the first value.
+ *
+ * @param a First value.
+ * @param b Second value.
+ *
+ * @return Reference to the first value.
+ */
+/// @{
+template <class T>
+constexpr inline quaternion<T>& operator-=(quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return (a = a - b);
+}
+template <class T>
+constexpr inline quaternion<T>& operator-=(quaternion<T>& a, T b) noexcept
+{
+	return (a = a - b);
+}
+/// @}
+
+/**
+ * Multiplies two values and stores the result in the first value.
+ *
+ * @param a First value.
+ * @param b Second value.
+ *
+ * @return Reference to the first value.
+ */
+/// @{
+template <class T>
+constexpr inline quaternion<T>& operator*=(quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return (a = a * b);
+}
+template <class T>
+constexpr inline quaternion<T>& operator*=(quaternion<T>& a, T b) noexcept
+{
+	return (a = a * b);
+}
+/// @}
+
+/**
+ * Divides the first value by the second value and stores the result in the first value.
+ *
+ * @param a First value.
+ * @param b Second value.
+ *
+ * @return Reference to the first value.
+ */
+/// @{
+template <class T>
+constexpr inline quaternion<T>& operator/=(quaternion<T>& a, const quaternion<T>& b) noexcept
+{
+	return (a = a / b);
+}
+template <class T>
+constexpr inline quaternion<T>& operator/=(quaternion<T>& a, T b) noexcept
+{
+	return (a = a / b);
+}
+/// @}
+
+/**
+ * Writes the real and imaginary parts of a quaternion to an output stream, with each number delimeted by a space.
+ *
+ * @param os Output stream.
+ * @param q Quaternion.
+ *
+ * @return Output stream.
+ */
+template <class T>
+std::ostream& operator<<(std::ostream& os, const math::quaternion<T>& q)
+{
+	os << q.r << ' ' << q.i;
+	return os;
+}
+
+/**
+ * Reads the real and imaginary parts of a quaternion from an input stream, with each number delimeted by a space.
+ *
+ * @param is Input stream.
+ * @param q Quaternion.
+ *
+ * @return Input stream.
+ */
+template <class T>
+std::istream& operator>>(std::istream& is, const math::quaternion<T>& q)
+{
+	is >> q.r;
+	is >> q.i;
+	return is;
+}
+
+} // namespace operators
+
+} // namespace math
+
+using namespace math::operators;
+
+#endif // ANTKEEPER_MATH_QUATERNION_HPP

src/math/se3.hpp

@@ -21,9 +21,7 @@
 #define ANTKEEPER_MATH_TRANSFORMATION_SE3_HPP
 
 #include "math/vector.hpp"
-#include "math/quaternion-type.hpp"
-#include "math/quaternion-operators.hpp"
-#include "math/quaternion-functions.hpp"
+#include "math/quaternion.hpp"
 
 namespace math {
 namespace transformation {

src/math/stream-operators.hpp

@@ -1,65 +0,0 @@
-/*
- * Copyright (C) 2021  Christopher J. Howard
- *
- * This file is part of Antkeeper source code.
- *
- * Antkeeper source code is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Antkeeper source code is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef ANTKEEPER_MATH_STREAM_OPERATORS_HPP
-#define ANTKEEPER_MATH_STREAM_OPERATORS_HPP
-
-#include "math/quaternion-type.hpp"
-#include <istream>
-#include <ostream>
-
-/**
- * Writes the real and imaginary parts of a quaternion to an output stream, with each number delimeted by a space.
- *
- * @param os Output stream.
- * @param q Quaternion.
- * @return Output stream.
- */
-template <class T>
-std::ostream& operator<<(std::ostream& os, const math::quaternion<T>& q);
-
-/**
- * Reads the real and imaginary parts of a quaternion from an input stream, with each number delimeted by a space.
- *
- * @param is Input stream.
- * @param q Quaternion.
- * @return Input stream.
- */
-template <class T>
-std::istream& operator>>(std::istream& is, const math::quaternion<T>& q);
-
-template <class T>
-std::ostream& operator<<(std::ostream& os, const math::quaternion<T>& q)
-{
-	os << q.w << ' ' << q.x() << ' ' << q.y() << ' ' << q.z();
-	return os;
-}
-
-template <class T>
-std::istream& operator>>(std::istream& is, const math::quaternion<T>& q)
-{
-	is >> q.w;
-	is >> q.x;
-	is >> q.y;
-	is >> q.z;
-	
-	return is;
-}
-
-#endif // ANTKEEPER_MATH_STREAM_OPERATORS_HPP

src/math/transform-functions.hpp

@@ -21,7 +21,7 @@
 #define ANTKEEPER_MATH_TRANSFORM_FUNCTIONS_HPP
 
 #include "math/transform-type.hpp"
-#include "math/quaternion-functions.hpp"
+#include "math/quaternion.hpp"
 
 namespace math {
 

src/math/transform-type.hpp

@@ -21,7 +21,7 @@
 #define ANTKEEPER_MATH_TRANSFORM_TYPE_HPP
 
 #include "math/vector.hpp"
-#include "math/quaternion-type.hpp"
+#include "math/quaternion.hpp"
 
 namespace math {
 
@@ -49,9 +49,9 @@ struct transform
 template <class T>
 constexpr transform<T> transform<T>::identity =
 {
-	{T(0), T(0), T(0)},
-	{T(1), T(0), T(0), T(0)},
-	{T(1), T(1), T(1)}
+	vector<T, 3>::zero(),
+	quaternion<T>::identity(),
+	vector<T, 3>::one()
 };
 
 } // namespace math

src/math/vector.hpp

@@ -293,7 +293,8 @@ using vector4 = vector;
 /**
  * Returns the absolute values of each element.
  *
- * @param x Input vector
+ * @param x Input vector.
+ *
  * @return Absolute values of input vector elements.
  */
 template <class T, std::size_t N>
@@ -304,6 +305,7 @@ constexpr vector abs(const vector& x);
  *
  * @param x First value.
  * @param y Second value.
+ *
  * @return Sum of the two values.
  */
 /// @{
@@ -317,6 +319,7 @@ constexpr vector add(const vector& x, T y) noexcept;
  * Checks if all elements of a boolean vector are `true`.
  *
  * @param x Vector to be tested for truth.
+ *
  * @return `true` if all elements are `true`, `false` otherwise.
  */
 template <std::size_t N>
@@ -326,6 +329,7 @@ constexpr bool all(const vector& x) noexcept;
  * Checks if any elements of a boolean vector are `true`.
  *
  * @param x Vector to be tested for truth.
+ *
  * @return `true` if any elements are `true`, `false` otherwise.
  */
 template <std::size_t N>
@@ -334,7 +338,8 @@ constexpr bool any(const vector& x) noexcept;
 /**
  * Performs a element-wise ceil operation.
  *
- * @param x Input vector
+ * @param x Input vector.
+ *
  * @return Component-wise ceil of input vector.
  */
 template <class T, std::size_t N>
@@ -346,6 +351,7 @@ constexpr vector ceil(const vector& x);
  * @param x Vector to clamp.
  * @param min Minimum value.
  * @param max Maximum value.
+ *
  * @return Clamped vector.
  */
 /// @{
@@ -360,6 +366,7 @@ constexpr vector clamp(const vector& x, T min, T max);
  *
  * @param x Vector to clamp.
  * @param max_length Maximum length.
+ *
  * @return Length-clamped vector.
  */
 template <class T, std::size_t N>
@@ -370,6 +377,7 @@ vector clamp_length(const vector& x, T max_length);
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Cross product of the two vectors.
  */
 template <class T>
@@ -380,6 +388,7 @@ constexpr vector cross(const vector& x, const vector& y) noexc
  *
  * @param p0 First of two points.
  * @param p1 Second of two points.
+ *
  * @return Distance between the two points.
  */
 template <class T, std::size_t N>
@@ -390,6 +399,7 @@ T distance(const vector& p0, const vector& p1);
  *
  * @param p0 First of two points.
  * @param p1 Second of two points.
+ *
  * @return Squared distance between the two points.
  */
 template <class T, std::size_t N>
@@ -400,6 +410,7 @@ constexpr T distance_squared(const vector& p0, const vector& p1) noe
  *
  * @param x First value.
  * @param y Second value.
+ *
  * @return Result of the division.
  */
 /// @{
@@ -416,6 +427,7 @@ constexpr vector div(T x, const vector& y) noexcept;
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Dot product of the two vectors.
  */
 template <class T, std::size_t N>
@@ -426,6 +438,7 @@ constexpr T dot(const vector& x, const vector& y) noexcept;
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
@@ -434,7 +447,8 @@ constexpr vector equal(const vector& x, const vector& y) no
 /**
  * Performs a element-wise floor operation.
  *
- * @param x Input vector
+ * @param x Input vector.
+ *
  * @return Component-wise floor of input vector.
  */
 template <class T, std::size_t N>
@@ -459,7 +473,8 @@ constexpr vector fma(const vector& x, T y, T z);
 /**
  * Returns a vector containing the fractional part of each element.
  *
- * @param x Input vector
+ * @param x Input vector.
+ *
  * @return Fractional parts of input vector.
  */
 template <class T, std::size_t N>
@@ -470,6 +485,7 @@ constexpr vector fract(const vector& x);
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
@@ -480,6 +496,7 @@ constexpr vector greater_than(const vector& x, const vector
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
@@ -489,6 +506,7 @@ constexpr vector greater_than_equal(const vector& x, const vector
  * Calculates the length of a vector.
  *
  * @param x Vector of which to calculate the length.
+ *
  * @return Length of the vector.
  */
 template <class T, std::size_t N>
@@ -498,6 +516,7 @@ T length(const vector& x);
  * Calculates the squared length of a vector. The squared length can be calculated faster than the length because a call to `std::sqrt` is saved.
  *
  * @param x Vector of which to calculate the squared length.
+ *
  * @return Squared length of the vector.
  */
 template <class T, std::size_t N>
@@ -508,6 +527,7 @@ constexpr T length_squared(const vector& x) noexcept;
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
@@ -518,6 +538,7 @@ constexpr vector less_than(const vector& x, const vector& y
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
@@ -570,6 +591,7 @@ constexpr T min(const vector& x);
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Remainders of `x / y`.
  */
 /// @{
@@ -584,6 +606,7 @@ constexpr vector mod(const vector& x, T y);
  *
  * @param x First value.
  * @param y Second value.
+ *
  * @return Product of the two values.
  */
 /// @{
@@ -597,6 +620,7 @@ constexpr vector mul(const vector& x, T y) noexcept;
  * Negates a vector.
  *
  * @param x Vector to negate.
+ *
  * @return Negated vector.
  */
 template <class T, std::size_t N>
@@ -606,6 +630,7 @@ constexpr vector negate(const vector& x) noexcept;
  * Calculates the unit vector in the same direction as the original vector.
  *
  * @param x Vector to normalize.
+ *
  * @return Normalized vector.
  */
 template <class T, std::size_t N>
@@ -615,6 +640,7 @@ vector normalize(const vector& x);
  * Logically inverts a boolean vector.
  *
  * @param x Vector to be inverted.
+ *
  * @return Logically inverted vector.
  */
 template <class T, std::size_t N>
@@ -625,6 +651,7 @@ constexpr vector not(const vector& x) noexcept;
  *
  * @param x First vector.
  * @param y Second vector.
+ *
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
@@ -649,6 +676,7 @@ vector pow(const vector& x, T y);
  * Performs a element-wise round operation.
  *
  * @param x Input vector
+ *
  * @return Component-wise round of input vector.
  */
 template <class T, std::size_t N>
@@ -678,6 +706,7 @@ vector sqrt(const vector& x);
  *
  * @param x First value.
  * @param y Second value.
+ *
  * @return Difference between the two values.
  */
 /// @{
@@ -1422,6 +1451,7 @@ constexpr inline vector operator-(T x, const vector& y) noexcept
  *
  * @param x First value.
  * @param y Second value.
+ *
  * @return Reference to the first value.
  */
 /// @{
@@ -1442,6 +1472,7 @@ constexpr inline vector& operator+=(vector& x, T y) noexcept
  *
  * @param x First value.
  * @param y Second value.
+ *
  * @return Reference to the first value.
  */
 /// @{
@@ -1462,6 +1493,7 @@ constexpr inline vector& operator-=(vector& x, T y) noexcept
  *
  * @param x First value.
  * @param y Second value.
+ *
  * @return Reference to the first value.
  */
 /// @{
@@ -1482,6 +1514,7 @@ constexpr inline vector& operator*=(vector& x, T y) noexcept
  *
  * @param x First value.
  * @param y Second value.
+ *
  * @return Reference to the first value.
  */
 /// @{

src/render/passes/material-pass.cpp

@@ -42,7 +42,7 @@
 #include "scene/point-light.hpp"
 #include "scene/spot-light.hpp"
 #include "config.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "math/projection.hpp"
 #include <cmath>
 #include <glad/glad.h>

src/render/passes/shadow-map-pass.cpp

@@ -34,7 +34,7 @@
 #include "config.hpp"
 #include "math/vector.hpp"
 #include "math/matrix.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "math/projection.hpp"
 #include <cmath>
 #include <glad/glad.h>

src/render/passes/sky-pass.cpp

@@ -70,9 +70,9 @@ sky_pass::sky_pass(gl::rasterizer* rasterizer, const gl::framebuffer* framebuffe
 	sun_luminance_tween(float3{0.0f, 0.0f, 0.0f}, math::lerp<float3, float>),
 	sun_illuminance_tween(float3{0.0f, 0.0f, 0.0f}, math::lerp<float3, float>),
 	icrf_to_eus_translation({0, 0, 0}, math::lerp<float3, float>),
-	icrf_to_eus_rotation(math::quaternion<float>::identity, math::nlerp<float>),
+	icrf_to_eus_rotation(math::quaternion<float>::identity(), math::nlerp<float>),
 	moon_position_tween(float3{0, 0, 0}, math::lerp<float3, float>),
-	moon_rotation_tween(math::quaternion<float>::identity, math::nlerp<float>),
+	moon_rotation_tween(math::quaternion<float>::identity(), math::nlerp<float>),
 	moon_angular_radius_tween(0.0f, math::lerp<float, float>),
 	moon_sunlight_direction_tween(float3{0, 0, 0}, math::lerp<float3, float>),
 	moon_sunlight_illuminance_tween(float3{0, 0, 0}, math::lerp<float3, float>),

src/render/passes/sky-pass.hpp

@@ -25,7 +25,7 @@
 #include "event/event-handler.hpp"
 #include "event/input-events.hpp"
 #include "animation/tween.hpp"
-#include "math/quaternion-type.hpp"
+#include "math/quaternion.hpp"
 #include "gl/shader-program.hpp"
 #include "gl/shader-input.hpp"
 #include "gl/vertex-buffer.hpp"

src/render/renderer.cpp

@@ -31,7 +31,7 @@
 #include "math/matrix.hpp"
 #include "geom/projection.hpp"
 #include "config.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include <functional>
 #include <set>
 

src/resources/entity-archetype-loader.cpp

@@ -275,10 +275,10 @@ static bool load_component_transform(entity::archetype& archetype, const json& e
 	if (element.contains("rotation"))
 	{
 		auto translation = element["rotation"];
-		component.local.rotation.w = translation[0].get<float>();
-		component.local.rotation.x = translation[1].get<float>();
-		component.local.rotation.y = translation[2].get<float>();
-		component.local.rotation.z = translation[3].get<float>();
+		component.local.rotation.w() = translation[0].get<float>();
+		component.local.rotation.x() = translation[1].get<float>();
+		component.local.rotation.y() = translation[2].get<float>();
+		component.local.rotation.z() = translation[3].get<float>();
 	}
 	
 	if (element.contains("scale"))

src/resources/model-loader.cpp

@@ -259,10 +259,10 @@ render::model* resource_loader::load(resource_manager* resource_m
 				{
 					if (rotation_node->size() == 4)
 					{
-						bone_transform.rotation.w = (*rotation_node)[0].get<float>();
-						bone_transform.rotation.x = (*rotation_node)[1].get<float>();
-						bone_transform.rotation.y = (*rotation_node)[2].get<float>();
-						bone_transform.rotation.z = (*rotation_node)[3].get<float>();
+						bone_transform.rotation.w() = (*rotation_node)[0].get<float>();
+						bone_transform.rotation.x() = (*rotation_node)[1].get<float>();
+						bone_transform.rotation.y() = (*rotation_node)[2].get<float>();
+						bone_transform.rotation.z() = (*rotation_node)[3].get<float>();
 					}
 				}
 				

src/scene/camera.cpp

@@ -21,7 +21,7 @@
 #include "config.hpp"
 #include "math/constants.hpp"
 #include "math/interpolation.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "math/projection.hpp"
 
 namespace scene {

src/scene/directional-light.cpp

@@ -19,7 +19,7 @@
 
 #include "directional-light.hpp"
 #include "config.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "math/interpolation.hpp"
 
 namespace scene {

src/scene/object.hpp

@@ -23,7 +23,7 @@
 #include "animation/tween.hpp"
 #include "geom/bounding-volume.hpp"
 #include "math/vector.hpp"
-#include "math/quaternion-type.hpp"
+#include "math/quaternion.hpp"
 #include "math/transform-type.hpp"
 #include "render/context.hpp"
 #include "render/queue.hpp"

src/scene/spot-light.cpp

@@ -19,7 +19,7 @@
 
 #include "spot-light.hpp"
 #include "config.hpp"
-#include "math/quaternion-operators.hpp"
+#include "math/quaternion.hpp"
 #include "math/interpolation.hpp"
 #include <cmath>