Integrate previously separate unpublished VMQ math library directly into Antkeeper source

CMakeLists.txt

@@ -5,7 +5,6 @@ option(VERSION_STRING "Project version string" "0.0.0")
 project(antkeeper VERSION ${VERSION_STRING} LANGUAGES CXX)
 
 # Find dependency packages
-find_package(vmq REQUIRED CONFIG)
 find_package(dr_wav REQUIRED CONFIG)
 find_package(stb REQUIRED CONFIG)
 find_package(glad REQUIRED CONFIG)
@@ -15,10 +14,8 @@ find_package(SDL2 REQUIRED COMPONENTS SDL2::SDL2-static SDL2::SDL2main CONFIG)
 find_package(OpenAL REQUIRED CONFIG)
 find_library(physfs REQUIRED NAMES physfs-static PATHS "${CMAKE_PREFIX_PATH}/lib")
 
-
 # Determine dependencies
 set(STATIC_LIBS
-	vmq
 	dr_wav
 	stb
 	glad

src/animation/ease.hpp

@@ -51,9 +51,8 @@
 #ifndef ANTKEEPER_EASE_HPP
 #define ANTKEEPER_EASE_HPP
 
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 #include <cmath>
-#include <type_traits>
 
 /**
  * Container for templated easing functions.
@@ -62,7 +61,7 @@
  *
  *     value_type operator+(const value_type&, const value_type&);
  *     value_type operator-(const value_type&, const value_type&);
- *     value_type operator*(const value_type&, scalar_type) const;
+ *     value_type operator*(const value_type&, scalar_type);
  *
  * @tparam T Value type.
  * @tparam S Scalar type.
@@ -72,12 +71,6 @@ struct ease
 {
 	typedef T value_type;
 	typedef S scalar_type;
-
-	/// Linearly interpolates between @p x and @p y.
-	static T linear(const T& x, const T& y, S a);
-	
-	/// Logarithmically interpolates between @p x and @p y.
-	static T log(const T& x, const T& y, S a);
 	
 	static T in_sine(const T& x, const T& y, S a);
 	static T out_sine(const T& x, const T& y, S a);
@@ -120,101 +113,88 @@ struct ease
 	static T in_out_bounce(const T& x, const T& y, S a);
 };
 
-template <typename T, typename S>
-inline T ease<T, S>::linear(const T& x, const T& y, S a)
-{
-	return (y - x) * a + x;
-}
-
-template <typename T, typename S>
-inline T ease<T, S>::log(const T& x, const T& y, S a)
-{
-	//return std::exp(linear(std::log(x), std::log(y), a));
-	return x * std::pow(y / x, a);
-}
-
 template <typename T, typename S>
 T ease<T, S>::in_sine(const T& x, const T& y, S a)
 {
-	return linear(y, x, std::cos(a * vmq::half_pi<S>));
+	return math::lerp(y, x, std::cos(a * math::half_pi<S>));
 }
 
 template <typename T, typename S>
 T ease<T, S>::out_sine(const T& x, const T& y, S a)
 {
-	return linear(x, y, std::sin(a * vmq::half_pi<S>));
+	return math::lerp(x, y, std::sin(a * math::half_pi<S>));
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_out_sine(const T& x, const T& y, S a)
 {
-	return linear(x, y, -(std::cos(a * vmq::pi<S>) - S(1)) * S(0.5));
+	return math::lerp(x, y, -(std::cos(a * math::pi<S>) - S(1)) * S(0.5));
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_quad(const T& x, const T& y, S a)
 {
-	return linear(x, y, a * a);
+	return math::lerp(x, y, a * a);
 }
 
 template <typename T, typename S>
 T ease<T, S>::out_quad(const T& x, const T& y, S a)
 {
-	return linear(x, y, (S(2) - a) * a);
+	return math::lerp(x, y, (S(2) - a) * a);
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_out_quad(const T& x, const T& y, S a)
 {
-	return linear(x, y, (a < S(0.5)) ? S(2) * a * a : -(S(2) * a * a - S(4) * a + S(1)));
+	return math::lerp(x, y, (a < S(0.5)) ? S(2) * a * a : -(S(2) * a * a - S(4) * a + S(1)));
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_cubic(const T& x, const T& y, S a)
 {
-	return linear(x, y, a * a * a);
+	return math::lerp(x, y, a * a * a);
 }
 
 template <typename T, typename S>
 T ease<T, S>::out_cubic(const T& x, const T& y, S a)
 {
-	return linear(x, y, a * ((a - S(3)) * a + S(3)));	
+	return math::lerp(x, y, a * ((a - S(3)) * a + S(3)));	
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_out_cubic(const T& x, const T& y, S a)
 {
-	return linear(x, y, (a < S(0.5)) ? S(4) * a * a * a : S(4) * a * a * a - S(12) * a * a + S(12) * a - 3);
+	return math::lerp(x, y, (a < S(0.5)) ? S(4) * a * a * a : S(4) * a * a * a - S(12) * a * a + S(12) * a - 3);
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_quart(const T& x, const T& y, S a)
 {
-	return linear(x, y, a * a * a * a);
+	return math::lerp(x, y, a * a * a * a);
 }
 
 template <typename T, typename S>
 T ease<T, S>::out_quart(const T& x, const T& y, S a)
 {
-	return linear(x, y, a * (a * ((S(4) - a) * a - S(6)) + S(4)));
+	return math::lerp(x, y, a * (a * ((S(4) - a) * a - S(6)) + S(4)));
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_out_quart(const T& x, const T& y, S a)
 {
-	return linear(x, y, (a < S(0.5)) ? S(8) * a * a * a * a : a * (a * ((S(32) - S(8) * a) * a - S(48)) + S(32)) - S(7));
+	return math::lerp(x, y, (a < S(0.5)) ? S(8) * a * a * a * a : a * (a * ((S(32) - S(8) * a) * a - S(48)) + S(32)) - S(7));
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_quint(const T& x, const T& y, S a)
 {
-	return linear(x, y, a * a * a * a * a);
+	return math::lerp(x, y, a * a * a * a * a);
 }
 
 template <typename T, typename S>
 T ease<T, S>::out_quint(const T& x, const T& y, S a)
 {
-	return linear(x, y, a * (a * (a * ((a - S(5)) * a + S(10)) - S(10)) + S(5)));
+	return math::lerp(x, y, a * (a * (a * ((a - S(5)) * a + S(10)) - S(10)) + S(5)));
 }
 
 template <typename T, typename S>
@@ -222,25 +202,25 @@ T ease::in_out_quint(const T& x, const T& y, S a)
 {
 	if (a < S(0.5))
 	{
-		return linear(x, y, S(16) * a * a * a * a * a);
+		return math::lerp(x, y, S(16) * a * a * a * a * a);
 	}
 	else
 	{
 		a = S(2) * (S(1) - a);
-		return linear(x, y, S(0.5) * (S(2) - a * a * a * a * a));
+		return math::lerp(x, y, S(0.5) * (S(2) - a * a * a * a * a));
 	}
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_expo(const T& x, const T& y, S a)
 {
-	return (a == S(0)) ? x : linear(x, y, std::pow(S(1024), a - S(1)));
+	return (a == S(0)) ? x : math::lerp(x, y, std::pow(S(1024), a - S(1)));
 }
 
 template <typename T, typename S>
 T ease<T, S>::out_expo(const T& x, const T& y, S a)
 {
-	return (a == S(1)) ? y : linear(y, x, std::pow(S(2), S(-10) * a));
+	return (a == S(1)) ? y : math::lerp(y, x, std::pow(S(2), S(-10) * a));
 }
 
 template <typename T, typename S>
@@ -255,19 +235,19 @@ T ease::in_out_expo(const T& x, const T& y, S a)
 		return y;
 	}
 	
-	return linear(x, y, (a < S(0.5)) ? std::pow(S(2), S(20) * a - S(11)) : S(1) - std::pow(S(2), S(9) - S(20) * a));
+	return math::lerp(x, y, (a < S(0.5)) ? std::pow(S(2), S(20) * a - S(11)) : S(1) - std::pow(S(2), S(9) - S(20) * a));
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_circ(const T& x, const T& y, S a)
 {
-	return linear(y, x, std::sqrt(S(1) - a * a));
+	return math::lerp(y, x, std::sqrt(S(1) - a * a));
 }
 
 template <typename T, typename S>
 T ease<T, S>::out_circ(const T& x, const T& y, S a)
 {
-	return linear(x, y, std::sqrt(-(a - S(2)) * a));
+	return math::lerp(x, y, std::sqrt(-(a - S(2)) * a));
 }
 
 template <typename T, typename S>
@@ -275,11 +255,11 @@ T ease::in_out_circ(const T& x, const T& y, S a)
 {
 	if (a < S(0.5))
 	{
-		return linear(x, y, S(0.5) - S(0.5) * std::sqrt(S(1) - S(4) * a * a));
+		return math::lerp(x, y, S(0.5) - S(0.5) * std::sqrt(S(1) - S(4) * a * a));
 	}
 	else
 	{
-		return linear(x, y, S(0.5) * (std::sqrt(S(-4) * (a - S(2)) * a - S(3)) + S(1)));
+		return math::lerp(x, y, S(0.5) * (std::sqrt(S(-4) * (a - S(2)) * a - S(3)) + S(1)));
 	}
 }
 
@@ -287,7 +267,7 @@ template
 T ease<T, S>::in_back(const T& x, const T& y, S a)
 {
 	const S c = S(1.70158);
-	return linear(x, y, a * a * (a * c + a - c));
+	return math::lerp(x, y, a * a * (a * c + a - c));
 }
 
 template <typename T, typename S>
@@ -295,7 +275,7 @@ T ease::out_back(const T& x, const T& y, S a)
 {
 	const S c = S(1.70158);
 	a -= S(1);
-	return linear(x, y, a * a * (a * c + a + c) + S(1));
+	return math::lerp(x, y, a * a * (a * c + a + c) + S(1));
 }
 
 template <typename T, typename S>
@@ -305,12 +285,12 @@ T ease::in_out_back(const T& x, const T& y, S a)
 	
 	if (a < S(0.5))
 	{
-		return linear(x, y, a * a * (a * (S(4) * c + S(4)) - S(2) * c));
+		return math::lerp(x, y, a * a * (a * (S(4) * c + S(4)) - S(2) * c));
 	}
 	else
 	{
 		S b = S(1) - S(2) * a;
-		return linear(x, y, b * b * (a * c + a - c * S(0.5) - S(1)) + S(1));
+		return math::lerp(x, y, b * b * (a * c + a - c * S(0.5) - S(1)) + S(1));
 	}
 }
 
@@ -326,7 +306,7 @@ T ease::in_elastic(const T& x, const T& y, S a)
 		return y;
 	}
 	
-	return linear(x, y, -std::pow(S(1024), a - S(1)) * std::sin(S(20.944) * (a - S(1.075))));
+	return math::lerp(x, y, -std::pow(S(1024), a - S(1)) * std::sin(S(20.944) * (a - S(1.075))));
 }
 
 template <typename T, typename S>
@@ -341,7 +321,7 @@ T ease::out_elastic(const T& x, const T& y, S a)
 		return y;
 	}
 	
-	return linear(x, y, std::pow(S(2), S(-10) * a) * std::sin(S(20.944) * (a - S(0.075))) + S(1));
+	return math::lerp(x, y, std::pow(S(2), S(-10) * a) * std::sin(S(20.944) * (a - S(0.075))) + S(1));
 }
 
 template <typename T, typename S>
@@ -358,19 +338,19 @@ T ease::in_out_elastic(const T& x, const T& y, S a)
 	
 	if (a < S(0.5))
 	{
-		return linear(x, y, std::pow(S(2), S(20) * a - S(11)) * std::sin(S(15.5334) - S(27.5293) * a));
+		return math::lerp(x, y, std::pow(S(2), S(20) * a - S(11)) * std::sin(S(15.5334) - S(27.5293) * a));
 		
 	}
 	else
 	{
-		return linear(y, x, std::pow(2, S(9) - S(20) * a) * std::sin(S(15.5334) - S(27.5293) * a));
+		return math::lerp(y, x, std::pow(2, S(9) - S(20) * a) * std::sin(S(15.5334) - S(27.5293) * a));
 	}
 }
 
 template <typename T, typename S>
 T ease<T, S>::in_bounce(const T& x, const T& y, S a)
 {
-	return linear(x, y, S(1) - ease<S, S>::out_bounce(S(0), S(1), S(1) - a));
+	return math::lerp(x, y, S(1) - ease<S, S>::out_bounce(S(0), S(1), S(1) - a));
 }
 
 template <typename T, typename S>
@@ -399,7 +379,7 @@ T ease::out_bounce(const T& x, const T& y, S a)
 		a = n * a * a + S(0.984375);
 	}
 	
-	return linear(x, y, a);
+	return math::lerp(x, y, a);
 }
 
 template <typename T, typename S>
@@ -407,11 +387,11 @@ T ease::in_out_bounce(const T& x, const T& y, S a)
 {
 	if (a < S(0.5))
 	{
-		return linear(x, y, (S(1) - ease<S, S>::out_bounce(S(0), S(1), S(1) - S(2) * a)) * S(0.5));
+		return math::lerp(x, y, (S(1) - ease<S, S>::out_bounce(S(0), S(1), S(1) - S(2) * a)) * S(0.5));
 	}
 	else
 	{
-		return linear(x, y, (S(1) + ease<S, S>::out_bounce(S(0), S(1), S(2) * a - S(1))) * S(0.5));
+		return math::lerp(x, y, (S(1) + ease<S, S>::out_bounce(S(0), S(1), S(2) * a - S(1))) * S(0.5));
 	}
 }
 

src/application.cpp

@@ -20,7 +20,6 @@
 #include "application.hpp"
 #include "configuration.hpp"
 #include "state/application-states.hpp"
-#include "math.hpp"
 
 // STL
 #include <cstdlib>
@@ -42,6 +41,9 @@
 #include "debug/ansi-codes.hpp"
 #include "debug/console-commands.hpp"
 
+// Math
+#include "math/math.hpp"
+
 // Resources
 #include "resources/resource-manager.hpp"
 
@@ -109,8 +111,7 @@
 // Utilities
 #include "utility/paths.hpp"
 #include "utility/timestamp.hpp"
-
-using namespace vmq::operators;
+#include "utility/fundamental-types.hpp"
 
 application::application(int argc, char** argv):
 	closed(false),
@@ -177,7 +178,7 @@ application::application(int argc, char** argv):
 	// Load configuration
 	//load_config();
 	fullscreen = true;
-	vsync = false;
+	vsync = true;
 	
 	// Parse command line options
 	parse_options(argc, argv);
@@ -518,7 +519,7 @@ application::application(int argc, char** argv):
 	overworld_compositor.add_pass(final_pass);
 	
 	// Setup overworld camera
-	overworld_camera.set_perspective(45.0f * vmq::pi<float> / 180.0f, (float)window_width / (float)window_height, 0.1f, 1000.0f);
+	overworld_camera.set_perspective(math::radians<float>(45.0f), (float)window_width / (float)window_height, 0.1f, 1000.0f);
 	overworld_camera.set_compositor(&overworld_compositor);
 	overworld_camera.set_composite_index(0);
 	overworld_camera.set_active(true);
@@ -548,7 +549,7 @@ application::application(int argc, char** argv):
 	underworld_compositor.add_pass(underworld_final_pass);
 	
 	// Setup underworld camera
-	underworld_camera.set_perspective(45.0f * vmq::pi<float> / 180.0f, (float)window_width / (float)window_height, 0.1f, 1000.0f);
+	underworld_camera.set_perspective(math::radians<float>(45.0f), (float)window_width / (float)window_height, 0.1f, 1000.0f);
 	underworld_camera.look_at({0, 50, 0}, {0, 0, 0}, {0, 0, -1});
 	underworld_camera.set_compositor(&underworld_compositor);
 	underworld_camera.set_composite_index(0);
@@ -751,11 +752,11 @@ application::application(int argc, char** argv):
 			ecs::cavity_component cavity;
 			cavity.position =
 			{
-				frand(-r, r),
-				frand(-r * 2, r),
-				frand(-r, r)
+				math::random(-r, r),
+				math::random(-r * 2, r),
+				math::random(-r, r)
 			};
-			cavity.radius = frand(0.75f, 1.25f);
+			cavity.radius = math::random(0.75f, 1.25f);
 
 			ecs_registry.assign<ecs::cavity_component>(ecs_registry.create(), cavity);
 		});
@@ -800,7 +801,7 @@ application::application(int argc, char** argv):
 	spotlight.set_color({1, 1, 1});
 	spotlight.set_intensity(1.0f);
 	spotlight.set_attenuation({1.0f, 0.09f, 0.032f});
-	spotlight.set_cutoff({vmq::radians(15.0f), vmq::radians(30.0f)});
+	spotlight.set_cutoff({math::radians(15.0f), math::radians(30.0f)});
 	spotlight.update_tweens();
 	spotlight.set_active(false);
 	
@@ -875,12 +876,9 @@ application::application(int argc, char** argv):
 	//model_instance* larva = new model_instance(resource_manager->load<model>("larva.obj"));
 	//underworld_scene.add_object(larva);
 	
-	quaternion<float> flashlight_rotation = vmq::angle_axis(vmq::half_pi<float>, {0.0f, 1.0f, 0.0f});
 	model_instance* flashlight = new model_instance(resource_manager->load<model>("flashlight.obj"));
-	flashlight->set_rotation(flashlight_rotation);
 	underworld_scene.add_object(flashlight);
 	model_instance* flashlight_light_cone = new model_instance(resource_manager->load<model>("flashlight-light-cone.obj"));
-	flashlight_light_cone->set_rotation(flashlight_rotation);
 	underworld_scene.add_object(flashlight_light_cone);
 	
 	control_system->set_flashlight(flashlight, flashlight_light_cone);
@@ -921,7 +919,7 @@ application::application(int argc, char** argv):
 	animator->add_animation(radial_transition_outer->get_animation());
 	
 	// Setup tweens
-	focal_point_tween.set_interpolator(ease<float3>::linear);
+	focal_point_tween.set_interpolator(math::lerp<float3, float>);
 	
 	// Register CLI commands
 	cli.register_command("echo", cc::echo);
@@ -974,7 +972,7 @@ int application::execute()
 
 	// Setup time tween
 	time[0] = time[1] = 0.0;
-	time.set_interpolator(ease<double>::linear);
+	time.set_interpolator(math::lerp<double, double>);
 
 	// Schedule frames until closed
 	while (!closed)
@@ -1033,19 +1031,24 @@ void application::setup_fsm()
 	initial_state = &splash_state;
 }
 
+void application::load_config()
+{
+	
+}
+
 void application::parse_options(int argc, char** argv)
 {
 	cxxopts::Options options("Antkeeper", "Ant colony simulation game");
 	
 	options.add_options()
-		("q,quick-start", "Skips splash screen")
-		("c,continue", "Continues from last save")
+		("c,continue", "Continues from the last save")
+		("d,data", "Sets the data package path", cxxopts::value<std::string>())
+		("f,fullscreen", "Starts in fullscreen mode")
 		("n,new-game", "Starts a new game")
+		("q,quick-start", "Skips to the main menu")
 		("r,reset", "Restores all settings to default")
-		("f,fullscreen", "Starts in fullscreen mode")
-		("w,windowed", "Starts in windowed mode")
 		("v,vsync", "Enables or disables v-sync", cxxopts::value<int>())
-		("d,data", "Specifies the data package path", cxxopts::value<std::string>())
+		("w,windowed", "Starts in windowed mode")
 		;
 	
 	logger.push_task("Parsing command line options");
@@ -1054,29 +1057,23 @@ void application::parse_options(int argc, char** argv)
 	{
 		auto result = options.parse(argc, argv);
 		
-		// --quick-start
-		if (result.count("quick-start"))
-		{
-			logger.log("Skipping splash screen");
-			initial_state = &play_state;
-		}
-		
 		// --continue
 		if (result.count("continue"))
 		{
 			logger.log("Continuing from last save");
 		}
 		
-		// --new-game
-		if (result.count("new-game"))
-		{
-			logger.log("Starting a new game");
-		}
-		
-		// --reset
-		if (result.count("reset"))
+		// --data
+		if (result.count("data"))
 		{
-			logger.log("Restoring all settings to default");
+			data_package_path = result["data"].as<std::string>();
+			
+			if (std::filesystem::path(data_package_path).is_relative())
+			{
+				data_package_path = data_path + data_package_path;
+			}
+			
+			logger.log("Set alternative data package path \"" + data_package_path + "\"");
 		}
 		
 		// --fullscreen
@@ -1086,11 +1083,23 @@ void application::parse_options(int argc, char** argv)
 			fullscreen = true;
 		}
 		
-		// --windowed
-		if (result.count("windowed"))
+		// --new-game
+		if (result.count("new-game"))
 		{
-			logger.log("Starting in windowed mode");
-			fullscreen = false;
+			logger.log("Starting a new game");
+		}
+		
+		// --quick-start
+		if (result.count("quick-start"))
+		{
+			logger.log("Skipping splash screen");
+			initial_state = &play_state;
+		}
+		
+		// --reset
+		if (result.count("reset"))
+		{
+			logger.log("Restoring all settings to default");
 		}
 		
 		// --vsync
@@ -1108,17 +1117,11 @@ void application::parse_options(int argc, char** argv)
 			}
 		}
 		
-		// --data
-		if (result.count("data"))
+		// --windowed
+		if (result.count("windowed"))
 		{
-			data_package_path = result["data"].as<std::string>();
-			
-			if (std::filesystem::path(data_package_path).is_relative())
-			{
-				data_package_path = data_path + data_package_path;
-			}
-			
-			logger.log("Set alternative data package path \"" + data_package_path + "\"");
+			logger.log("Starting in windowed mode");
+			fullscreen = false;
 		}
 	}
 	catch (const std::exception& e)

src/application.hpp

@@ -188,6 +188,7 @@ public:
 
 private:
 	void setup_fsm();
+	void load_config();
 	void parse_options(int argc, char** argv);
 	
 

src/camera-rig.cpp

@@ -19,16 +19,15 @@
 
 #include "orbit-cam.hpp"
 #include "scene/camera.hpp"
+#include "math/constants.hpp"
 #include "configuration.hpp"
 #include <algorithm>
 #include <cmath>
 
-using namespace vmq::operators;
-
 camera_rig::camera_rig():
 	camera(nullptr),
 	translation{0.0f, 0.0f, 0.0f},
-	rotation(vmq::identity_quaternion<float>),
+	rotation(math::identity_quaternion<float>),
 	forward(global_forward),
 	right(global_right),
 	up(global_up)
@@ -53,7 +52,7 @@ void camera_rig::set_translation(const float3& translation)
 	this->translation = translation;
 }
 
-void camera_rig::set_rotation(const vmq::quaternion<float>& rotation)
+void camera_rig::set_rotation(const quaternion_type& rotation)
 {
 	this->rotation = rotation;
 	
@@ -62,4 +61,3 @@ void camera_rig::set_rotation(const vmq::quaternion& rotation)
 	up = rotation * global_up;
 	right = rotation * global_right;
 }
-

src/camera-rig.hpp

@@ -20,9 +20,9 @@
 #ifndef ANTKEEPER_CAMERA_RIG_HPP
 #define ANTKEEPER_CAMERA_RIG_HPP
 
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "math/quaternion-type.hpp"
+#include "math/transform-type.hpp"
+#include "utility/fundamental-types.hpp"
 
 class camera;
 
@@ -32,6 +32,9 @@ class camera;
 class camera_rig
 {
 public:
+	typedef math::quaternion<float> quaternion_type;
+	typedef math::transform<float> transform_type;
+	
 	camera_rig();
 	
 	/**
@@ -61,7 +64,7 @@ public:
 	/**
 	 * Sets the rotation of the camera rig.
 	 */
-	void set_rotation(const vmq::quaternion<float>& rotation);
+	void set_rotation(const quaternion_type& rotation);
 	
 	/**
 	 * Returns the attached camera.
@@ -72,7 +75,7 @@ public:
 	::camera* get_camera();
 	
 	const float3& get_translation() const;
-	const vmq::quaternion<float>& get_rotation() const;
+	const quaternion_type& get_rotation() const;
 	const float3& get_forward() const;
 	const float3& get_right() const;
 	const float3& get_up() const;
@@ -80,7 +83,7 @@ public:
 private:
 	camera* camera;
 	float3 translation;
-	vmq::quaternion<float> rotation;
+	quaternion_type rotation;
 	float3 forward;
 	float3 right;
 	float3 up;
@@ -101,7 +104,7 @@ inline const float3& camera_rig::get_translation() const
 	return translation;
 }
 
-inline const vmq::quaternion<float>& camera_rig::get_rotation() const
+inline const typename camera_rig::quaternion_type& camera_rig::get_rotation() const
 {
 	return rotation;
 }

src/configuration.hpp.in

@@ -20,18 +20,16 @@
 #ifndef ANTKEEPER_CONFIGURATION_HPP
 #define ANTKEEPER_CONFIGURATION_HPP
 
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "math/vector-type.hpp"
 
 #define ANTKEEPER_VERSION_MAJOR @PROJECT_VERSION_MAJOR@
 #define ANTKEEPER_VERSION_MINOR @PROJECT_VERSION_MINOR@
 #define ANTKEEPER_VERSION_PATCH @PROJECT_VERSION_PATCH@
 #define ANTKEEPER_VERSION_STRING "@PROJECT_VERSION@"
 
-constexpr float3 global_forward = {0.0f, 0.0f, -1.0f};
-constexpr float3 global_up = {0.0f, 1.0f, 0.0f};
-constexpr float3 global_right = {1.0f, 0.0f, 0.0f};
+constexpr math::vector<float, 3> global_forward = {0.0f, 0.0f, -1.0f};
+constexpr math::vector<float, 3> global_up = {0.0f, 1.0f, 0.0f};
+constexpr math::vector<float, 3> global_right = {1.0f, 0.0f, 0.0f};
 
 #define MATERIAL_PASS_MAX_AMBIENT_LIGHT_COUNT 1
 #define MATERIAL_PASS_MAX_POINT_LIGHT_COUNT 1

src/entity/components/cavity-component.hpp

@@ -20,8 +20,7 @@
 #ifndef ANTKEEPER_ECS_CAVITY_COMPONENT_HPP
 #define ANTKEEPER_ECS_CAVITY_COMPONENT_HPP
 
-#include <vmq/vmq.hpp>
-using namespace vmq::types;
+#include "math/math.hpp"
 
 namespace ecs {
 

src/entity/components/locomotion-component.hpp

@@ -21,8 +21,7 @@
 #define ANTKEEPER_ECS_PLACEMENT_COMPONENT_HPP
 
 #include "geometry/mesh.hpp"
-#include <vmq/vmq.hpp>
-using namespace vmq::types;
+#include "utility/fundamental-types.hpp"
 
 namespace ecs {
 

src/entity/components/samara-component.hpp

@@ -20,8 +20,7 @@
 #ifndef ANTKEEPER_ECS_SAMARA_COMPONENT_HPP
 #define ANTKEEPER_ECS_SAMARA_COMPONENT_HPP
 
-#include <vmq/vmq.hpp>
-using namespace vmq::types;
+#include "utility/fundamental-types.hpp"
 
 namespace ecs {
 

src/entity/components/transform-component.hpp

@@ -20,13 +20,13 @@
 #ifndef ANTKEEPER_ECS_TRANSFORM_COMPONENT_HPP
 #define ANTKEEPER_ECS_TRANSFORM_COMPONENT_HPP
 
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 
 namespace ecs {
 
 struct transform_component
 {
-	vmq::transform<float> transform;
+	math::transform<float> transform;
 	bool warp;
 };
 

src/entity/entity-commands.hpp

@@ -20,18 +20,14 @@
 #ifndef ANTKEEPER_ECS_ENTITY_COMMANDS_HPP
 #define ANTKEEPER_ECS_ENTITY_COMMANDS_HPP
 
+#include "utility/fundamental-types.hpp"
 #include <entt/entt.hpp>
-#include <vmq/vmq.hpp>
-
-using namespace vmq;
 
 namespace ecs {
 	
 void move_to(entt::registry& registry, entt::entity entity, const float3& position);
 void warp_to(entt::registry& registry, entt::entity entity, const float3& position);
-
 void assign_render_layers(entt::registry& registry, entt::entity entity, unsigned int layers);
-
 void bind_transform(entt::registry& registry, entt::entity entity, entt::entity target);
 
 } // namespace ecs

src/geometry/aabb.hpp

@@ -22,20 +22,18 @@
 
 #include "bounding-volume.hpp"
 #include "sphere.hpp"
-#include <vmq/vmq.hpp>
-#include <array>
+#include "math/math.hpp"
 #include <limits>
 
-using vmq::vector;
-using vmq::matrix;
-using vmq::transform;
-using namespace vmq::operators;
-
 template <class T>
 struct aabb: public bounding_volume<T>
 {
-	vector<T, 3> min_point;
-	vector<T, 3> max_point;
+	typedef math::vector<T, 3> vector_type;
+	typedef math::matrix<T, 4, 4> matrix_type;
+	typedef math::transform<T> transform_type;
+	
+	vector_type min_point;
+	vector_type max_point;
 	
 	/**
 	 * Transforms an AABB.
@@ -44,7 +42,7 @@ struct aabb: public bounding_volume
 	 * @param t Transform by which the AABB should be transformed.
 	 * @return Transformed AABB.
 	 */
-	static aabb transform(const aabb& a, const ::transform<T>& t);
+	static aabb transform(const aabb& a, const transform_type& t);
 	
 	/**
 	 * Transforms an AABB.
@@ -53,9 +51,9 @@ struct aabb: public bounding_volume
 	 * @param m Matrix by which the AABB should be transformed.
 	 * @return Transformed AABB.
 	 */
-	static aabb transform(const aabb& a, const matrix<T, 4, 4>& m);
+	static aabb transform(const aabb& a, const matrix_type& m);
 	
-	aabb(const vector<T, 3>& min_point, const vector<T, 3>& max_point);
+	aabb(const vector_type& min_point, const vector_type& max_point);
 	aabb();
 	
 	virtual bounding_volume_type get_bounding_volume_type() const;
@@ -63,7 +61,7 @@ struct aabb: public bounding_volume
 	virtual bool intersects(const aabb<T>& aabb) const;
 	virtual bool contains(const sphere<T>& sphere) const;
 	virtual bool contains(const aabb<T>& aabb) const;
-	virtual bool contains(const vector<T, 3>& point) const;
+	virtual bool contains(const vector_type& point) const;
 
 	/**
 	 * Returns the position of the specified corner.
@@ -71,18 +69,18 @@ struct aabb: public bounding_volume
 	 * @param index Index of a corner.
 	 * @return Position of the specified corner.
 	 */
-	vector<T, 3> corner(int index) const noexcept;
+	vector_type corner(int index) const noexcept;
 };
 
 template <class T>
-aabb<T> aabb<T>::transform(const aabb& a, const ::transform<T>& t)
+aabb<T> aabb<T>::transform(const aabb& a, const transform_type& t)
 {
-	vector<T, 3> min_point = {std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity()};
-	vector<T, 3> max_point = {-std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()};
+	vector_type min_point = {std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity()};
+	vector_type max_point = {-std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()};
 	
 	for (std::size_t i = 0; i < 8; ++i)
 	{
-		vector<T, 3> transformed_corner = vmq::mul(t, a.corner(i));
+		vector_type transformed_corner = math::mul(t, a.corner(i));
 
 		for (std::size_t j = 0; j < 3; ++j)
 		{
@@ -95,15 +93,15 @@ aabb aabb::transform(const aabb& a, const ::transform& t)
 }
 
 template <class T>
-aabb<T> aabb<T>::transform(const aabb& a, const matrix<T, 4, 4>& m)
+aabb<T> aabb<T>::transform(const aabb& a, const matrix_type& m)
 {
-	vector<T, 3> min_point = {std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity()};
-	vector<T, 3> max_point = {-std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()};
+	vector_type min_point = {std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity()};
+	vector_type max_point = {-std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()};
 	
 	for (std::size_t i = 0; i < 8; ++i)
 	{
-		vector<T, 3> corner = a.corner(i);
-		vector<T, 4> transformed_corner = vmq::mul(m, vector<T, 4>{corner.x, corner.y, corner.z, T(1)});
+		vector_type corner = a.corner(i);
+		vector<T, 4> transformed_corner = math::mul(m, vector<T, 4>{corner.x, corner.y, corner.z, T(1)});
 
 		for (std::size_t j = 0; j < 3; ++j)
 		{
@@ -116,7 +114,7 @@ aabb aabb::transform(const aabb& a, const matrix& m)
 }
 
 template <class T>
-aabb<T>::aabb(const vector<T, 3>& min_point, const vector<T, 3>& max_point):
+aabb<T>::aabb(const vector_type& min_point, const vector_type& max_point):
 	min_point(min_point),
 	max_point(max_point)
 {}
@@ -134,7 +132,7 @@ inline bounding_volume_type aabb::get_bounding_volume_type() const
 template <class T>
 bool aabb<T>::intersects(const sphere<T>& sphere) const
 {
-	const vector<T, 3> radius_vector = {sphere.radius, sphere.radius, sphere.radius};
+	const vector_type radius_vector = {sphere.radius, sphere.radius, sphere.radius};
 	return aabb<T>(min_point - radius_vector, max_point + radius_vector).contains(sphere.center);
 }
 
@@ -175,7 +173,7 @@ bool aabb::contains(const aabb& aabb) const
 }
 
 template <class T>
-bool aabb<T>::contains(const vector<T, 3>& point) const
+bool aabb<T>::contains(const vector_type& point) const
 {
 	if (point.x < min_point.x || point.x > max_point.x)
 		return false;
@@ -187,7 +185,7 @@ bool aabb::contains(const vector& point) const
 }
 
 template <class T>
-vector<T, 3> aabb<T>::corner(int index) const noexcept
+typename aabb<T>::vector_type aabb<T>::corner(int index) const noexcept
 {
 	return
 		{

src/geometry/bounding-volume.hpp

@@ -20,10 +20,10 @@
 #ifndef ANTKEEPER_BOUNDING_VOLUME_HPP
 #define ANTKEEPER_BOUNDING_VOLUME_HPP
 
+#include "math/math.hpp"
 #include <stdexcept>
-#include <vmq/vmq.hpp>
 
-using vmq::vector;
+using math::vector;
 
 template <class T>
 struct sphere;
@@ -91,4 +91,3 @@ bool bounding_volume::intersects(const bounding_volume& volume) const
 }
 
 #endif // ANTKEEPER_BOUNDING_VOLUME_HPP
-

src/geometry/convex-hull.hpp

@@ -20,12 +20,12 @@
 #ifndef ANTKEEPER_CONVEX_HULL_HPP
 #define ANTKEEPER_CONVEX_HULL_HPP
 
+#include "bounding-volume.hpp"
+#include "geometry/plane.hpp"
+#include "geometry/sphere.hpp"
+#include "geometry/aabb.hpp"
 #include <cstdlib>
 #include <vector>
-#include "bounding-volume.hpp"
-#include "plane.hpp"
-#include "sphere.hpp"
-#include "aabb.hpp"
 
 /**
  * A plane-bounded convex hull.
@@ -73,7 +73,7 @@ template
 bool convex_hull<T>::intersects(const sphere<T>& sphere) const
 {
 	for (const plane<T>& plane: planes)
-		if (signed_distance(plane, sphere.center) < -sphere.radius)
+		if (plane.signed_distance(sphere.center) < -sphere.radius)
 			return false;
 	return true;
 }
@@ -87,7 +87,7 @@ bool convex_hull::intersects(const aabb& aabb) const
 		pv.x = (plane.normal.x > T(0)) ? aabb.max_point.x : aabb.min_point.x;
 		pv.y = (plane.normal.y > T(0)) ? aabb.max_point.y : aabb.min_point.y;
 		pv.z = (plane.normal.z > T(0)) ? aabb.max_point.z : aabb.min_point.z;
-		if (signed_distance(plane, pv) < T(0))
+		if (plane.signed_distance(pv) < T(0))
 			return false;
 	}
 	
@@ -98,7 +98,7 @@ template
 bool convex_hull<T>::contains(const sphere<T>& sphere) const
 {
 	for (const plane<T>& plane: planes)
-		if (signed_distance(plane, sphere.center) < sphere.radius)
+		if (plane.signed_distance(sphere.center) < sphere.radius)
 			return false;
 	return true;
 }
@@ -118,7 +118,7 @@ bool convex_hull::contains(const aabb& aabb) const
 		nv.y = (plane.normal.y < T(0)) ? aabb.max_point.y : aabb.min_point.y;
 		nv.z = (plane.normal.z < T(0)) ? aabb.max_point.z : aabb.min_point.z;
 		
-		if (signed_distance(plane, pv) < T(0) || signed_distance(plane, nv) < T(0))
+		if (plane.signed_distance(pv) < T(0) || plane.signed_distance(nv) < T(0))
 			return false;
 	}
 	
@@ -129,7 +129,7 @@ template
 bool convex_hull<T>::contains(const vector<T, 3>& point) const
 {
 	for (const plane<T>& plane: planes)
-		if (signed_distance(plane, point) < T(0))
+		if (plane.signed_distance(point) < T(0))
 			return false;
 
 	return true;

src/geometry/csg.hpp

@@ -20,15 +20,11 @@
 #ifndef ANTKEEPER_CSG_HPP
 #define ANTKEEPER_CSG_HPP
 
+#include "utility/fundamental-types.hpp"
 #include <list>
-#include <vmq/vmq.hpp>
-
 
 namespace csg {
 
-using vmq::float3;
-using vmq::float4;
-
 struct plane
 {
 	float3 normal;

src/geometry/intersection.cpp

@@ -22,10 +22,10 @@
 
 std::tuple<bool, float> ray_plane_intersection(const ray<float>& ray, const plane<float>& plane)
 {
-	float denom = vmq::dot(ray.direction, plane.normal);
+	float denom = math::dot(ray.direction, plane.normal);
 	if (denom != 0.0f)
 	{
-		float t = -(vmq::dot(ray.origin, plane.normal) + plane.distance) / denom;
+		float t = -(math::dot(ray.origin, plane.normal) + plane.distance) / denom;
 		
 		if (t >= 0.0f)
 		{
@@ -43,8 +43,8 @@ std::tuple ray_triangle_intersection(const ray
 	float3 edge20 = c - a;
 
 	// Calculate determinant
-	float3 pv = vmq::cross(ray.direction, edge20);
-	float det = vmq::dot(edge10, pv);
+	float3 pv = math::cross(ray.direction, edge20);
+	float det = math::dot(edge10, pv);
 
 	if (!det)
 	{
@@ -55,7 +55,7 @@ std::tuple ray_triangle_intersection(const ray
 
 	// Calculate u
 	float3 tv = ray.origin - a;
-	float u = vmq::dot(tv, pv) * inverse_det;
+	float u = math::dot(tv, pv) * inverse_det;
 	
 	if (u < 0.0f || u > 1.0f)
 	{
@@ -63,8 +63,8 @@ std::tuple ray_triangle_intersection(const ray
 	}
 
 	// Calculate v
-	float3 qv = vmq::cross(tv, edge10);
-	float v = vmq::dot(ray.direction, qv) * inverse_det;
+	float3 qv = math::cross(tv, edge10);
+	float v = math::dot(ray.direction, qv) * inverse_det;
 
 	if (v < 0.0f || u + v > 1.0f)
 	{
@@ -72,7 +72,7 @@ std::tuple ray_triangle_intersection(const ray
 	}
 
 	// Calculate t
-	float t = vmq::dot(edge20, qv) * inverse_det;
+	float t = math::dot(edge20, qv) * inverse_det;
 
 	if (t > 0.0f)
 	{

src/geometry/intersection.hpp

@@ -20,16 +20,12 @@
 #ifndef ANTKEEPER_INTERSECTION_HPP
 #define ANTKEEPER_INTERSECTION_HPP
 
+#include "geometry/aabb.hpp"
+#include "geometry/mesh.hpp"
+#include "geometry/plane.hpp"
+#include "geometry/ray.hpp"
+#include "utility/fundamental-types.hpp"
 #include <tuple>
-#include <vmq/vmq.hpp>
-#include "aabb.hpp"
-#include "mesh.hpp"
-#include "plane.hpp"
-#include "ray.hpp"
-
-using namespace vmq::types;
-using namespace vmq::operators;
-
 
 /**
  * Tests for intersection between a ray and a plane.

src/geometry/mesh-accelerator.hpp

@@ -20,16 +20,14 @@
 #ifndef ANTKEEPER_MESH_ACCELERATOR_HPP
 #define ANTKEEPER_MESH_ACCELERATOR_HPP
 
-#include "mesh.hpp"
+#include "geometry/mesh.hpp"
 #include "geometry/octree.hpp"
 #include "geometry/aabb.hpp"
-#include "intersection.hpp"
+#include "geometry/intersection.hpp"
+#include "utility/fundamental-types.hpp"
 #include <list>
-#include <unordered_map>
 #include <optional>
-#include <vmq/vmq.hpp>
-using namespace vmq::types;
-using namespace vmq::operators;
+#include <unordered_map>
 
 /**
  * Acceleration structure for querying mesh geometry.

src/geometry/mesh-functions.cpp

@@ -18,13 +18,8 @@
  */
 
 #include "mesh-functions.hpp"
-#include <vmq/vmq.hpp>
-#include <stdexcept>
+#include "math/math.hpp"
 #include <unordered_map>
-#include <map>
-#include <iostream>
-
-using namespace vmq::operators;
 
 struct edge_hasher
 {
@@ -60,10 +55,6 @@ void create_triangle_mesh(mesh& mesh, const std::vector& vertices, const
 
 			if (auto it = edge_map.find({start->index, end->index}); it != edge_map.end())
 			{
-				/*
-				if (it->second->face)
-					std::cout << "THIS EDGE ALREADY HAS A FACE!\n" << std::endl;
-				*/
 				loop[j] = it->second;
 			}
 			else
@@ -92,7 +83,7 @@ void calculate_face_normals(float* normals, const mesh& mesh)
 		const float3& b = reinterpret_cast<const float3&>(face.edge->next->vertex->position);
 		const float3& c = reinterpret_cast<const float3&>(face.edge->previous->vertex->position);
 
-		normal = vmq::normalize(vmq::cross(b - a, c - a));
+		normal = math::normalize(math::cross(b - a, c - a));
 	}
 }
 
@@ -118,5 +109,3 @@ aabb calculate_bounds(const mesh& mesh)
 
 	return aabb<float>{bounds_min, bounds_max};
 }
-
-

src/geometry/mesh-functions.hpp

@@ -20,13 +20,11 @@
 #ifndef ANTKEEPER_MESH_FUNCTIONS_HPP
 #define ANTKEEPER_MESH_FUNCTIONS_HPP
 
-#include "mesh.hpp"
-#include <vmq/vmq.hpp>
+#include "geometry/aabb.hpp"
+#include "geometry/mesh.hpp"
+#include "utility/fundamental-types.hpp"
 #include <array>
 #include <vector>
-#include "geometry/aabb.hpp"
-
-using namespace vmq::types;
 
 /**
  * Creates a triangle mesh from a list of vertices and indices.

src/geometry/plane.hpp

@@ -20,93 +20,97 @@
 #ifndef ANTKEEPER_PLANE_HPP
 #define ANTKEEPER_PLANE_HPP
 
-#include <vmq/vmq.hpp>
-
-using vmq::vector;
-using namespace vmq::operators;
+#include "math/math.hpp"
 
 template <class T>
 struct plane
 {
-	vector<T, 3> normal;
+	typedef math::vector<T, 3> vector_type;
+	
+	/// Plane normal vector.
+	vector_type normal;
+	
+	/// Plane distance.
 	T distance;
 	
 	/**
 	 * Creates a plane given a normal vector and distance.
 	 */
-	plane(const vector<T, 3>& normal, T distance);
+	plane(const vector_type& normal, T distance);
 	
 	/**
 	 * Creates a plane given a normal vector and offset vector.
 	 */
-	plane(const vector<T, 3>& normal, const vector<T, 3>& offset);
+	plane(const vector_type& normal, const vector_type& offset);
 	
 	/**
 	 * Creates a plane given three points.
 	 */
-	plane(const vector<T, 3>& a, const vector<T, 3>& b, const vector<T, 3>& c);
+	plane(const vector_type& a, const vector_type& b, const vector_type& c);
 	
 	/**
 	 * Creates a plane given its coefficients.
 	 *
 	 * @param coefficients Vector containing the plane coefficients, A, B, C and D, as x, y, z, and w, respectively.
 	 */
-	plane(const vector<T, 4>& coefficients);
+	plane(const math::vector<T, 4>& coefficients);
 	
 	/// Creates an uninitialized plane.
 	plane() = default;
+	
+	/**
+	 * Calculates the signed distance between a plane and a vector.
+	 *
+	 * @param v Vector.
+	 * @return Signed distance between the plane and vector.
+	 */
+	T signed_distance(const vector_type& v) const;
+	
+	/**
+	 * Calculates the point of intersection between three planes.
+	 */
+	static vector_type intersection(const plane& p0, const plane& p1, const plane& p2);
 };
 
 template <class T>
-inline plane<T>::plane(const vector<T, 3>& normal, T distance):
+inline plane<T>::plane(const vector_type& normal, T distance):
 	normal(normal),
 	distance(distance)
 {}
 
 template <class T>
-plane<T>::plane(const vector<T, 3>& normal, const vector<T, 3>& offset):
+plane<T>::plane(const vector_type& normal, const vector_type& offset):
 	normal(normal),
-	distance(-vmq::dot(normal, offset))
+	distance(-math::dot(normal, offset))
 {}
 
 template <class T>
-plane<T>::plane(const vector<T, 3>& a, const vector<T, 3>& b, const vector<T, 3>& c)
+plane<T>::plane(const vector_type& a, const vector_type& b, const vector_type& c)
 {
-	normal = vmq::normalize(vmq::cross(c - b, a - b));
-	distance = -(vmq::dot(normal, b));
+	normal = math::normalize(math::cross(c - b, a - b));
+	distance = -(math::dot(normal, b));
 }
 
 template <class T>
-plane<T>::plane(const vector<T, 4>& coefficients)
+plane<T>::plane(const math::vector<T, 4>& coefficients)
 {
-	const vector<T, 3> abc = vmq::resize<3>(coefficients);
-	const float inverse_length = T(1) / vmq::length(abc);
+	const vector_type abc = math::resize<3>(coefficients);
+	const float inverse_length = T(1) / math::length(abc);
 	
 	normal = abc * inverse_length;
 	distance = coefficients[3] * inverse_length;
 }
 
-/**
- * Calculates the signed distance between a plane and a vector.
- *
- * @param p Plane.
- * @param v Vector.
- * @return Signed distance between the plane and a vector.
- */
 template <class T>
-inline T signed_distance(const plane<T>& p, const vector<T, 3>& v)
+inline T plane<T>::signed_distance(const vector_type& v) const
 {
-	return p.distance + vmq::dot(p.normal, v);
+	return distance + math::dot(normal, v);
 }
 
-/**
- * Calculates the point of intersection between three planes.
- */
 template <class T>
-vector<T, 3> intersection(const plane<T>& p0, const plane<T>& p1, const plane<T>& p2)
+typename plane<T>::vector_type plane<T>::intersection(const plane& p0, const plane& p1, const plane& p2)
 {
-	return -(p0.distance * vmq::cross(p1.normal, p2.normal) + p1.distance * vmq::cross(p2.normal, p0.normal) + p2.distance * vmq::cross(p0.normal, p1.normal)) / vmq::dot(p0.normal, vmq::cross(p1.normal, p2.normal));
+	return -(p0.distance * math::cross(p1.normal, p2.normal) + p1.distance * math::cross(p2.normal, p0.normal) + p2.distance * math::cross(p0.normal, p1.normal)) / math::dot(p0.normal, math::cross(p1.normal, p2.normal));
 }
 
 #endif // ANTKEEPER_PLANE_HPP
-

src/geometry/projection.hpp

@@ -0,0 +1,31 @@
+/*
+ * Copyright (C) 2020  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_PROJECTION_HPP
+#define ANTKEEPER_PROJECTION_HPP
+
+#include "math/math.hpp"
+
+template <class T>
+math::vector<T, 3> project_on_plane(const math::vector<T, 3>& v, const math::vector<T, 3>& p, const math::vector<T, 3>& n)
+{
+	return v - n * math::dot(v - p, n);
+}
+
+#endif // ANTKEEPER_PROJECTION_HPP

src/geometry/ray.hpp

@@ -20,22 +20,30 @@
 #ifndef ANTKEEPER_RAY_HPP
 #define ANTKEEPER_RAY_HPP
 
-#include <vmq/vmq.hpp>
-
-using vmq::vector;
-using namespace vmq::operators;
+#include "math/math.hpp"
 
 template <class T>
 struct ray
 {
-	vector<T, 3> origin;
-	vector<T, 3> direction;
-
-	vector<T, 3> extrapolate(T distance) const;
+	typedef math::vector<T, 3> vector_type;
+	
+	/// Origin of the ray.
+	vector_type origin;
+	
+	/// Normalized direction vector of the ray.
+	vector_type direction;
+
+	/**
+	 * Extrapolates from the ray origin along the ray direction vector.
+	 *
+	 * @param distance Distance to extrapolate.
+	 * @return Extrapolated coordinates.
+	 */
+	vector_type extrapolate(T distance) const;
 };
 
 template <class T>
-inline vector<T, 3> ray<T>::extrapolate(T distance) const
+inline typename ray<T>::vector_type ray<T>::extrapolate(T distance) const
 {
 	return origin + direction * distance;
 }

src/geometry/sdf.hpp

@@ -20,9 +20,8 @@
 #ifndef ANTKEEPER_SDF_HPP
 #define ANTKEEPER_SDF_HPP
 
-#include <vmq/vmq.hpp>
-using namespace vmq::types;
-using namespace vmq::operators;
+#include "utility/fundamental-types.hpp"
+#include <algorithm>
 
 /**
  * Contains signed distance functions.
@@ -36,14 +35,14 @@ inline float3 translate(const float3& sample, const float3& offset)
 
 inline float sphere(const float3& p, float r)
 {
-	return vmq::length(p) - r;
+	return math::length(p) - r;
 }
 
 inline float cylinder(const float3& p, float r, float h)
 {
-	float dx = std::abs(vmq::length(vmq::swizzle<0, 2>(p))) - r;
+	float dx = std::abs(math::length(math::swizzle<0, 2>(p))) - r;
 	float dy = std::abs(p[1]) - h;
-	return std::min<float>(std::max<float>(dx, dy), 0.0f) + vmq::length(float2{std::max<float>(dx, 0.0f), std::max<float>(dy, 0.0f)});
+	return std::min<float>(std::max<float>(dx, dy), 0.0f) + math::length(float2{std::max<float>(dx, 0.0f), std::max<float>(dy, 0.0f)});
 }
 
 inline float op_union(float a, float b)

src/geometry/sphere.hpp

@@ -20,21 +20,20 @@
 #ifndef ANTKEEPER_SPHERE_HPP
 #define ANTKEEPER_SPHERE_HPP
 
-#include "bounding-volume.hpp"
-#include "aabb.hpp"
-#include <vmq/vmq.hpp>
+#include "geometry/bounding-volume.hpp"
+#include "geometry/aabb.hpp"
+#include "math/math.hpp"
 #include <algorithm>
 
-using vmq::vector;
-using namespace vmq::operators;
-
 template <class T>
 struct sphere: public bounding_volume<T>
 {
-	vector<T, 3> center;
+	typedef math::vector<T, 3> vector_type;
+	
+	vector_type center;
 	T radius;
 	
-	sphere(const vector<T, 3>& center, T radius);
+	sphere(const vector_type& center, T radius);
 	sphere();
 	
 	virtual bounding_volume_type get_bounding_volume_type() const;
@@ -42,11 +41,11 @@ struct sphere: public bounding_volume
 	virtual bool intersects(const aabb<T>& aabb) const;
 	virtual bool contains(const sphere<T>& sphere) const;
 	virtual bool contains(const aabb<T>& aabb) const;
-	virtual bool contains(const vector<T, 3>& point) const;
+	virtual bool contains(const vector_type& point) const;
 };
 
 template <class T>
-sphere<T>::sphere(const vector<T, 3>& center, T radius):
+sphere<T>::sphere(const vector_type& center, T radius):
 	center(center),
 	radius(radius)
 {}
@@ -64,9 +63,9 @@ inline bounding_volume_type sphere::get_bounding_volume_type() const
 template <class T>
 bool sphere<T>::intersects(const sphere<T>& sphere) const
 {
-	vector<T, 3> d = center - sphere.center;
+	vector_type d = center - sphere.center;
 	float r = radius + sphere.radius;
-	return (vmq::dot(d, d) <= r * r);
+	return (math::dot(d, d) <= r * r);
 }
 
 template <class T>
@@ -82,8 +81,8 @@ bool sphere::contains(const sphere& sphere) const
 	if (containment_radius < T(0))
 		return false;
 	
-	vector<T, 3> d = center - sphere.center;
-	return (vmq::dot(d, d) <= containment_radius * containment_radius);
+	vector_type d = center - sphere.center;
+	return (math::dot(d, d) <= containment_radius * containment_radius);
 }
 
 template <class T>
@@ -91,8 +90,8 @@ bool sphere::contains(const aabb& aabb) const
 {
 	T distance = T(0);
 	
-	vector<T, 3> a = center - aabb.min_point;
-	vector<T, 3> b = center - aabb.max_point;
+	vector_type a = center - aabb.min_point;
+	vector_type b = center - aabb.max_point;
 	
 	distance += std::max(a.x * a.x, b.x * b.x);
 	distance += std::max(a.y * a.y, b.y * b.y);
@@ -102,10 +101,10 @@ bool sphere::contains(const aabb& aabb) const
 }
 
 template <class T>
-bool sphere<T>::contains(const vector<T, 3>& point) const
+bool sphere<T>::contains(const vector_type& point) const
 {
-	vector<T, 3> d = center - point;
-	return (vmq::dot(d, d) <= radius * radius);
+	vector_type d = center - point;
+	return (math::dot(d, d) <= radius * radius);
 }
 
 #endif // ANTKEEPER_SPHERE_HPP

src/geometry/view-frustum.hpp

@@ -20,24 +20,24 @@
 #ifndef ANTKEEPER_VIEW_FRUSTUM_HPP
 #define ANTKEEPER_VIEW_FRUSTUM_HPP
 
-#include "convex-hull.hpp"
+#include "geometry/convex-hull.hpp"
+#include "math/math.hpp"
 #include <array>
-#include <vmq/vmq.hpp>
-
-using vmq::vector;
-using vmq::matrix;
-using namespace vmq::operators;
 
 template <class T>
 class view_frustum
 {
 public:
+	typedef math::vector<T, 3> vector_type;
+	typedef math::matrix<T, 4, 4> matrix_type;
+	typedef plane<T> plane_type;
+	
 	/**
 	 * Creates a view frustum from a view-projection matrix.
 	 *
 	 * @param view_projection View-projection matrix.
 	 */
-	view_frustum(const matrix<T, 4, 4>& view_projection);
+	view_frustum(const matrix_type& view_projection);
 	
 	/// Creates an uninitialized view frustum.
 	view_frustum();
@@ -47,46 +47,46 @@ public:
 	 *
 	 * @param view_projection View-projection matrix.
 	 */
-	void set_matrix(const matrix<T, 4, 4>& view_projection);
+	void set_matrix(const matrix_type& view_projection);
 	
 	/// Returns a convex hull which describes the bounds of the view frustum.
 	const convex_hull<T>& get_bounds() const;
 	
 	/// Returns the left clipping plane.
-	const plane<T>& get_left() const;
+	const plane_type& get_left() const;
 	
 	/// Returns the right clipping plane.
-	const plane<T>& get_right() const;
+	const plane_type& get_right() const;
 	
 	/// Returns the bottom clipping plane.
-	const plane<T>& get_bottom() const;
+	const plane_type& get_bottom() const;
 	
 	/// Returns the top clipping plane.
-	const plane<T>& get_top() const;
+	const plane_type& get_top() const;
 	
 	/// Returns the near clipping plane.
-	const plane<T>& get_near() const;
+	const plane_type& get_near() const;
 	
 	/// Returns the far clipping plane.
-	const plane<T>& get_far() const;
+	const plane_type& get_far() const;
 	
 	/**
 	 * Returns an array containing the corners of the view frustum bounds.
 	 *
 	 * @return Array containing the corners of the view frustum bounds. Corners are stored in the following order: NTL, NTR, NBL, NBR, FTL, FTR, FBL, FBR; where N is near, F is far, T is top, B is bottom, L is left, and R is right, therefore NTL refers to the corner shared between the near, top, and left clipping planes.
 	 */
-	const std::array<vector<T, 3>, 8>& get_corners() const;
+	const std::array<vector_type, 8>& get_corners() const;
 
 private:
-	void recalculate_planes(const matrix<T, 4, 4>& view_projection);
+	void recalculate_planes(const matrix_type& view_projection);
 	void recalculate_corners();
 	
 	convex_hull<T> bounds;
-	std::array<vector<T, 3>, 8> corners;
+	std::array<vector_type, 8> corners;
 };
 
 template <class T>
-view_frustum<T>::view_frustum(const matrix<T, 4, 4>& view_projection):
+view_frustum<T>::view_frustum(const matrix_type& view_projection):
 	bounds(6)
 {
 	set_matrix(view_projection);
@@ -94,11 +94,11 @@ view_frustum::view_frustum(const matrix& view_projection):
 
 template <class T>
 view_frustum<T>::view_frustum():
-	view_frustum(vmq::identity4x4<T>)
+	view_frustum(math::identity4x4<T>)
 {}
 
 template <class T>
-void view_frustum<T>::set_matrix(const matrix<T, 4, 4>& view_projection)
+void view_frustum<T>::set_matrix(const matrix_type& view_projection)
 {
 	recalculate_planes(view_projection);
 	recalculate_corners();
@@ -111,70 +111,70 @@ inline const convex_hull& view_frustum::get_bounds() const
 }
 
 template <class T>
-inline const plane<T>& view_frustum<T>::get_left() const
+inline const typename view_frustum<T>::plane_type& view_frustum<T>::get_left() const
 {
 	return bounds.planes[0];
 }
 
 template <class T>
-inline const plane<T>& view_frustum<T>::get_right() const
+inline const typename view_frustum<T>::plane_type& view_frustum<T>::get_right() const
 {
 	return bounds.planes[1];
 }
 
 template <class T>
-inline const plane<T>& view_frustum<T>::get_bottom() const
+inline const typename view_frustum<T>::plane_type& view_frustum<T>::get_bottom() const
 {
 	return bounds.planes[2];
 }
 
 template <class T>
-inline const plane<T>& view_frustum<T>::get_top() const
+inline const typename view_frustum<T>::plane_type& view_frustum<T>::get_top() const
 {
 	return bounds.planes[3];
 }
 
 template <class T>
-inline const plane<T>& view_frustum<T>::get_near() const
+inline const typename view_frustum<T>::plane_type& view_frustum<T>::get_near() const
 {
 	return bounds.planes[4];
 }
 
 template <class T>
-inline const plane<T>& view_frustum<T>::get_far() const
+inline const typename view_frustum<T>::plane_type& view_frustum<T>::get_far() const
 {
 	return bounds.planes[5];
 }
 
 template <class T>
-inline const std::array<vector<T, 3>, 8>& view_frustum<T>::get_corners() const
+inline const std::array<typename view_frustum<T>::vector_type, 8>& view_frustum<T>::get_corners() const
 {
 	return corners;
 }
 
 template <class T>
-void view_frustum<T>::recalculate_planes(const matrix<T, 4, 4>& view_projection)
+void view_frustum<T>::recalculate_planes(const matrix_type& view_projection)
 {
-	matrix<T, 4, 4> transpose = vmq::transpose(view_projection);
-	bounds.planes[0] = plane<T>(transpose[3] + transpose[0]);
-	bounds.planes[1] = plane<T>(transpose[3] - transpose[0]);
-	bounds.planes[2] = plane<T>(transpose[3] + transpose[1]);
-	bounds.planes[3] = plane<T>(transpose[3] - transpose[1]);
-	bounds.planes[4] = plane<T>(transpose[3] + transpose[2]);
-	bounds.planes[5] = plane<T>(transpose[3] - transpose[2]);
+	matrix_type transpose = math::transpose(view_projection);
+	bounds.planes[0] = plane_type(transpose[3] + transpose[0]);
+	bounds.planes[1] = plane_type(transpose[3] - transpose[0]);
+	bounds.planes[2] = plane_type(transpose[3] + transpose[1]);
+	bounds.planes[3] = plane_type(transpose[3] - transpose[1]);
+	bounds.planes[4] = plane_type(transpose[3] + transpose[2]);
+	bounds.planes[5] = plane_type(transpose[3] - transpose[2]);
 }
 
 template <class T>
 void view_frustum<T>::recalculate_corners()
 {
-	corners[0] = intersection(get_near(), get_top(), get_left());
-	corners[1] = intersection(get_near(), get_top(), get_right());
-	corners[2] = intersection(get_near(), get_bottom(), get_left());
-	corners[3] = intersection(get_near(), get_bottom(), get_right());
-	corners[4] = intersection(get_far(), get_top(), get_left());
-	corners[5] = intersection(get_far(), get_top(), get_right());
-	corners[6] = intersection(get_far(), get_bottom(), get_left());
-	corners[7] = intersection(get_far(), get_bottom(), get_right());
+	corners[0] = plane_type::intersection(get_near(), get_top(), get_left());
+	corners[1] = plane_type::intersection(get_near(), get_top(), get_right());
+	corners[2] = plane_type::intersection(get_near(), get_bottom(), get_left());
+	corners[3] = plane_type::intersection(get_near(), get_bottom(), get_right());
+	corners[4] = plane_type::intersection(get_far(), get_top(), get_left());
+	corners[5] = plane_type::intersection(get_far(), get_top(), get_right());
+	corners[6] = plane_type::intersection(get_far(), get_bottom(), get_left());
+	corners[7] = plane_type::intersection(get_far(), get_bottom(), get_right());
 }
 
 #endif // ANTKEEPER_VIEW_FRUSTUM_HPP

src/math.hpp

@@ -1,75 +0,0 @@
-/*
- * Copyright (C) 2020  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_HPP
-#define ANTKEEPER_MATH_HPP
-
-#include <cstdlib>
-#include <type_traits>
-#include <vmq/vmq.hpp>
-
-using namespace vmq::operators;
-
-inline float frand(float start, float end)
-{
-    float f = (float)std::rand() / RAND_MAX;
-    return start + f * (end - start);
-}
-
-namespace math {
-
-/**
- * Reinterprets data as an `N`-dimensional vector of type `T`.
- *
- * @tparam N Number of vector dimensions.
- * @tparam T Scalar type.
- * @param data Data to reinterpret.
- */
-template <std::size_t N, typename T>
-inline vmq::vector<T, N>& as_vector(T& data)
-{
-	static_assert(std::is_pod<vmq::vector<T, N>>::value);
-	return reinterpret_cast<vmq::vector<T, N>&>(data);
-}
-
-/**
- * Reinterprets data as an `NxM` matrix of type `T`.
- *
- * @tparam N Number of columns.
- * @tparam M Number of rows.
- * @tparam T Element type.
- * @param data Data to reinterpret.
- */
-template <std::size_t N, std::size_t M, typename T>
-inline vmq::matrix<T, N, M>& as_matrix(T& data)
-{
-	static_assert(std::is_pod<vmq::matrix<T, N, M>>::value);
-	return reinterpret_cast<vmq::matrix<T, N, M>&>(data);
-}
-
-template <class T>
-vmq::vector<T, 3> project_on_plane(const vmq::vector<T, 3>& v, const vmq::vector<T, 3>& p, const vmq::vector<T, 3>& n)
-{
-	return v - n * vmq::dot(v - p, n);
-}
-
-} // namespace math
-
-#endif // ANTKEEPER_MATH_HPP
-

src/math/angles.hpp

@@ -0,0 +1,64 @@
+/*
+ * Copyright (C) 2020  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_ANGLES_HPP
+#define ANTKEEPER_MATH_ANGLES_HPP
+
+#include "math/constants.hpp"
+
+namespace math {
+
+/// @addtogroup utility
+/// @{
+
+/**
+ * Converts an angle given in radians to degrees.
+ *
+ * @param radians Angle in radians.
+ * @return Angle in degrees.
+ */
+template <class T>
+T degrees(T radians);
+
+/**
+ * Converts an angle given in degrees to radians.
+ *
+ * @param radians Angle in radians.
+ * @return Angle in degrees.
+ */
+template <class T>
+T radians(T degrees);
+
+template <class T>
+inline T degrees(T radians)
+{
+	return radians * T(180) / pi<T>;
+}
+
+template <class T>
+inline T radians(T degrees)
+{
+	return degrees * pi<T> / T(180);
+}
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_ANGLES_HPP

src/math/constants.hpp

@@ -0,0 +1,104 @@
+/*
+ * Copyright (C) 2020  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_CONSTANTS_HPP
+#define ANTKEEPER_MATH_CONSTANTS_HPP
+
+#include "math/matrix-type.hpp"
+#include "math/quaternion-type.hpp"
+#include "math/transform-type.hpp"
+
+namespace math {
+
+/// @addtogroup utility
+/// @{
+
+/**
+ * Pi.
+ */
+template <class T>
+constexpr T pi = T(3.14159265358979323846);
+
+/**
+ * Pi / 2.
+ */
+template <class T>
+constexpr T half_pi = pi<T> * T(0.5);
+
+/**
+ * Pi * 2.
+ */
+template <class T>
+constexpr T two_pi = pi<T> * T(2);
+
+/**
+ * 2x2 identity matrix.
+ */
+template <class T>
+constexpr matrix<T, 2, 2> identity2x2 =
+{{
+	{1, 0},
+	{0, 1}
+}};
+
+/**
+ * 3x3 identity matrix.
+ */
+template <class T>
+constexpr matrix<T, 3, 3> identity3x3 =
+{{
+	{1, 0, 0},
+	{0, 1, 0},
+	{0, 0, 1}
+}};
+
+/**
+ * 4x4 identity matrix.
+ */
+template <class T>
+constexpr matrix<T, 4, 4> identity4x4 =
+{{
+	{1, 0, 0, 0},
+	{0, 1, 0, 0},
+	{0, 0, 1, 0},
+	{0, 0, 0, 1}
+}};
+
+/**
+ * Unit quaternion.
+ */
+template <class T>
+constexpr quaternion<T> identity_quaternion = {T(1), T(0), T(0), T(0)};
+
+/**
+ * Identity transform.
+ */
+template <class T>
+constexpr transform<T> identity_transform =
+{
+	{T(0), T(0), T(0)},
+	{T(1), T(0), T(0), T(0)},
+	{T(1), T(1), T(1)}
+};
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_CONSTANTS_HPP

src/math/interpolation.hpp

@@ -0,0 +1,77 @@
+/*
+ * Copyright (C) 2020  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_INTERPOLATION_HPP
+#define ANTKEEPER_MATH_INTERPOLATION_HPP
+
+#include "math/matrix-type.hpp"
+#include "math/quaternion-type.hpp"
+#include "math/transform-type.hpp"
+#include <cmath>
+#include <type_traits>
+
+namespace math {
+
+/// @addtogroup utility
+/// @{
+
+/**
+ * Linearly interpolates between @p x and @p y.
+ *
+ * Requires the following operators to be defined:
+ *
+ *     T operator+(const T&, const T&);
+ *     T operator-(const T&, const T&);
+ *     T operator*(const T&, S);
+ *
+ * @tparam T Value type.
+ * @tparam S Scalar type.
+ */
+template <typename T, typename S = float>
+T lerp(const T& x, const T& y, S a);
+
+/**
+ * Logarithmically interpolates between @p x and @p y.
+ *
+ * @warning Undefined behavior when @p x is zero.
+ *
+ * @tparam T Value type.
+ * @tparam S Scalar type.
+ */
+template <typename T, typename S = float>
+T log_lerp(const T& x, const T& y, S a);
+
+template <typename T, typename S>
+inline T lerp(const T& x, const T& y, S a)
+{
+	return (y - x) * a + x;
+}
+
+template <typename T, typename S>
+inline T log_lerp(const T& x, const T& y, S a)
+{
+	//return std::exp(linear(std::log(x), std::log(y), a));
+	return x * std::pow(y / x, a);
+}
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_INTERPOLATION_HPP

src/math/math.hpp

@@ -0,0 +1,79 @@
+/*
+ * Copyright (C) 2020  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_HPP
+#define ANTKEEPER_MATH_HPP
+
+/// Contains mathematical functions and data types.
+namespace math {}
+
+/**
+ * @defgroup vector Vector
+ *
+ * Vector type, functions, and operators.
+ */
+#include "math/vector-type.hpp"
+#include "math/vector-functions.hpp"
+#include "math/vector-operators.hpp"
+
+/**
+ * @defgroup matrix Matrix
+ *
+ * Matrix type, functions, and operators.
+ */
+#include "math/matrix-type.hpp"
+#include "math/matrix-functions.hpp"
+#include "math/matrix-operators.hpp"
+
+/**
+ * @defgroup quaternion Quaternion
+ *
+ * Quaternion type, functions, and operators.
+ */
+#include "math/quaternion-type.hpp"
+#include "math/quaternion-functions.hpp"
+#include "math/quaternion-operators.hpp"
+
+/**
+ * @defgroup transform Transform
+ *
+ * TRS transform type, functions, and operators.
+ */
+#include "math/transform-type.hpp"
+#include "math/transform-functions.hpp"
+#include "math/transform-operators.hpp"
+
+/**
+ * @defgroup io I/O
+ *
+ * Functions and operators that read/write vectors, matrices, or quaternions from/to streams.
+ */
+#include "math/stream-operators.hpp"
+
+/**
+ * @defgroup utility Utility constants and functions
+ *
+ * Commonly used utilities.
+ */
+#include "math/angles.hpp"
+#include "math/constants.hpp"
+#include "math/interpolation.hpp"
+#include "math/random.hpp"
+
+#endif // ANTKEEPER_MATH_HPP

src/math/matrix-functions.hpp

@@ -0,0 +1,795 @@
+/*
+ * Copyright (C) 2020  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_MATRIX_FUNCTIONS_HPP
+#define ANTKEEPER_MATH_MATRIX_FUNCTIONS_HPP
+
+#include "math/matrix-type.hpp"
+#include "math/vector-type.hpp"
+#include "math/vector-functions.hpp"
+#include <type_traits>
+
+namespace math {
+
+/// @addtogroup matrix
+/// @{
+
+/**
+ * Adds two matrices.
+ *
+ * @param x First matrix.
+ * @param y Second matrix.
+ * @return Sum of the two matrices.
+ */
+template <class T>
+matrix<T, 2, 2> add(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y);
+
+/// @copydoc add(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 3, 3> add(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y);
+
+/// @copydoc add(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 4, 4> add(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y);
+
+/**
+ * Reinterprets data as an `NxM` matrix of type `T`.
+ *
+ * @tparam N Number of columns.
+ * @tparam M Number of rows.
+ * @tparam T Element type.
+ * @param data Data to reinterpret.
+ */
+template <std::size_t N, std::size_t M, typename T>
+matrix<T, N, M>& as_matrix(T& data);
+
+/**
+ * Calculates the determinant of a matrix.
+ *
+ * @param m Matrix of which to take the determinant. 
+ */
+template <class T>
+T determinant(const matrix<T, 2, 2>& m);
+
+/// @copydoc determinant(const matrix<T, 2, 2>&)
+template <class T>
+T determinant(const matrix<T, 3, 3>& m);
+
+/// @copydoc determinant(const matrix<T, 2, 2>&)
+template <class T>
+T determinant(const matrix<T, 4, 4>& m);
+
+/**
+ * Calculates the inverse of a matrix.
+ *
+ * @param m Matrix of which to take the inverse.
+ */
+template <class T>
+matrix<T, 2, 2> inverse(const matrix<T, 2, 2>& m);
+
+/// @copydoc inverse(const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 3, 3> inverse(const matrix<T, 3, 3>& m);
+
+/// @copydoc inverse(const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 4, 4> inverse(const matrix<T, 4, 4>& m);
+
+/**
+ * Performs a component-wise multiplication of two matrices.
+ *
+ * @param x First matrix multiplicand.
+ * @param y Second matrix multiplicand.
+ */
+template <class T>
+matrix<T, 2, 2> componentwise_mul(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y);
+
+/// @copydoc componentwise_mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 3, 3> componentwise_mul(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y);
+
+/// @copydoc componentwise_mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 4, 4> componentwise_mul(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y);
+
+/**
+ * Creates a viewing transformation matrix.
+ *
+ * @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>
+matrix<T, 4, 4> look_at(const vector<T, 3>& position, const vector<T, 3>& target, vector<T, 3> up);
+
+/**
+ * Multiplies two matrices.
+ *
+ * @param x First matrix.
+ * @param y Second matrix.
+ * @return Product of the two matrices.
+ */
+template <class T>
+matrix<T, 2, 2> mul(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y);
+
+/// @copydoc mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&);
+template <class T>
+matrix<T, 3, 3> mul(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y);
+
+/// @copydoc mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&);
+template <class T>
+matrix<T, 4, 4> mul(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y);
+
+/**
+ * Multiplies a matrix by a scalar.
+ *
+ * @param m Matrix.
+ * @param s Scalar.
+ * @return Product of the matrix and the scalar..
+ */
+template <class T, std::size_t N, std::size_t M>
+matrix<T, N, M> mul(const matrix<T, N, M>& m, T s);
+
+/**
+ * Transforms a vector by a matrix.
+ *
+ * @param m Matrix.
+ * @param v Vector.
+ * @return Transformed vector.
+ */
+template <class T>
+vector<T, 2> mul(const matrix<T, 2, 2>& m, const vector<T, 2>& v);
+
+/// @copydoc mul(const matrix<T, 2, 2>&, const vector<T, 2>&)
+template <class T>
+vector<T, 3> mul(const matrix<T, 3, 3>& m, const vector<T, 3>& v);
+
+/// @copydoc mul(const matrix<T, 2, 2>&, const vector<T, 2>&)
+template <class T>
+vector<T, 4> mul(const matrix<T, 4, 4>& m, const vector<T, 4>& v);
+
+/**
+ * Creates an orthographic projection matrix.
+ *
+ * @param left Signed distance to the left clipping plane.
+ * @param right Signed distance to the right clipping plane.
+ * @param bottom Signed distance to the bottom clipping plane.
+ * @param top Signed distance to the top clipping plane.
+ * @param z_near Signed distance to the near clipping plane.
+ * @param z_far Signed distance to the far clipping plane.
+ * @return Orthographic projection matrix.
+ */
+template <class T>
+matrix<T, 4, 4> ortho(T left, T right, T bottom, T top, T z_near, T z_far);
+
+/**
+ * Calculates the outer product of a pair of vectors.
+ *
+ * @param c Parameter to be treated as a column vector.
+ * @param r Parameter to be treated as a row vector.
+ */
+template <class T>
+matrix<T, 2, 2> outer_product(const vector<T, 2>& c, const vector<T, 2>& r);
+
+/// @copydoc outer_product(const vector<T, 2>&, const vector<T, 2>&)
+template <class T>
+matrix<T, 3, 3> outer_product(const vector<T, 3>& c, const vector<T, 3>& r);
+
+/// @copydoc outer_product(const vector<T, 2>&, const vector<T, 2>&)
+template <class T>
+matrix<T, 4, 4> outer_product(const vector<T, 4>& c, const vector<T, 4> r);
+
+/**
+ * Creates a perspective projection matrix.
+ *
+ * @param vertical_fov Vertical field of view angle, in radians.
+ * @param aspect_ratio Aspect ratio which determines the horizontal field of view.
+ * @param z_near Distance to the near clipping plane.
+ * @param z_far Distance to the far clipping plane.
+ * @return Perspective projection matrix.
+ */
+template <class T>
+matrix<T, 4, 4> perspective(T vertical_fov, T aspect_ratio, T z_near, T z_far);
+
+/**
+ * Resizes a matrix. Any new elements will be set to `1` if in the diagonal, and `0` otherwise.
+ *
+ * @param m Matrix to resize.
+ * @return Resized matrix.
+ */
+template <std::size_t N1, std::size_t M1, class T, std::size_t N0, std::size_t M0>
+matrix<T, N1, M1> resize(const matrix<T, N0, M0>& m);
+
+/**
+ * Rotates a matrix.
+ *
+ * @param m Matrix to rotate.
+ * @param angle Angle of rotation (in radians).
+ * @param axis Axis of rotation
+ * @return Rotated matrix.
+ */
+template <class T>
+matrix<T, 4, 4> rotate(const matrix<T, 4, 4>& m, T angle, const vector<T, 3>& axis);
+
+/**
+ * Scales a matrix.
+ *
+ * @param m Matrix to scale.
+ * @param v Scale vector.
+ * @return Scaled matrix.
+ */
+template <class T>
+matrix<T, 4, 4> scale(const matrix<T, 4, 4>& m, const vector<T, 3>& v);
+
+/**
+ * Subtracts a matrix from another matrix.
+ *
+ * @param x First matrix.
+ * @param y Second matrix.
+ * @return Difference between the two matrices.
+ */
+template <class T>
+matrix<T, 2, 2> sub(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y);
+
+/// @copydoc sub(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 3, 3> sub(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y);
+
+/// @copydoc sub(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 4, 4> sub(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y);
+
+/**
+ * Translates a matrix.
+ *
+ * @param m Matrix to translate.
+ * @param v Translation vector.
+ * @return Translated matrix.
+ */
+template <class T>
+matrix<T, 4, 4> translate(const matrix<T, 4, 4>& m, const vector<T, 3>& v);
+
+/**
+ * Calculates the transpose of a matrix.
+ *
+ * @param m Matrix of which to take the transpose.
+ */
+template <class T>
+matrix<T, 2, 2> transpose(const matrix<T, 2, 2>& m);
+
+/// @copydoc transpose(const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 3, 3> transpose(const matrix<T, 3, 3>& m);
+
+/// @copydoc transpose(const matrix<T, 2, 2>&)
+template <class T>
+matrix<T, 4, 4> transpose(const matrix<T, 4, 4>& m);
+
+template <class T>
+matrix<T, 2, 2> add(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y)
+{
+	return
+		{{
+			x[0] + y[0],
+			x[1] + y[1]
+		}};
+}
+
+template <class T>
+matrix<T, 3, 3> add(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y)
+{
+	return
+		{{
+			x[0] + y[0],
+			x[1] + y[1],
+			x[2] + y[2]
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> add(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y)
+{
+	return
+		{{
+			x[0] + y[0],
+			x[1] + y[1],
+			x[2] + y[2],
+			x[3] + y[3]
+		}};
+}
+
+template <std::size_t N, std::size_t M, typename T>
+inline matrix<T, N, M>& as_matrix(T& data)
+{
+	static_assert(std::is_pod<matrix<T, N, M>>::value);
+	return reinterpret_cast<matrix<T, N, M>&>(data);
+}
+
+template <class T>
+T determinant(const matrix<T, 2, 2>& m)
+{
+	return m[0][0] * m[1][1] - m[0][1] * m[1][0];
+}
+
+template <class T>
+T determinant(const matrix<T, 3, 3>& m)
+{
+	return m[0][0] * m [1][1] * m[2][2] +
+		m[0][1] * m[1][2] * m[2][0] +
+		m[0][2] * m[1][0] * m[2][1] -
+		m[0][0] * m[1][2] * m[2][1] -
+		m[0][1] * m[1][0] * m[2][2] -
+		m[0][2] * m[1][1] * m[2][0];
+}
+
+template <class T>
+T determinant(const matrix<T, 4, 4>& m)
+{
+	return m[0][3] * m[1][2] * m[2][1] * m[3][0] - m[0][2] * m[1][3] * m[2][1] * m[3][0] -
+		m[0][3] * m[1][1] * m[2][2] * m[3][0] + m[0][1] * m[1][3] * m[2][2] * m[3][0] +
+		m[0][2] * m[1][1] * m[2][3] * m[3][0] - m[0][1] * m[1][2] * m[2][3] * m[3][0] -
+		m[0][3] * m[1][2] * m[2][0] * m[3][1] + m[0][2] * m[1][3] * m[2][0] * m[3][1] +
+		m[0][3] * m[1][0] * m[2][2] * m[3][1] - m[0][0] * m[1][3] * m[2][2] * m[3][1] -
+		m[0][2] * m[1][0] * m[2][3] * m[3][1] + m[0][0] * m[1][2] * m[2][3] * m[3][1] +
+		m[0][3] * m[1][1] * m[2][0] * m[3][2] - m[0][1] * m[1][3] * m[2][0] * m[3][2] -
+		m[0][3] * m[1][0] * m[2][1] * m[3][2] + m[0][0] * m[1][3] * m[2][1] * m[3][2] +
+		m[0][1] * m[1][0] * m[2][3] * m[3][2] - m[0][0] * m[1][1] * m[2][3] * m[3][2] -
+		m[0][2] * m[1][1] * m[2][0] * m[3][3] + m[0][1] * m[1][2] * m[2][0] * m[3][3] +
+		m[0][2] * m[1][0] * m[2][1] * m[3][3] - m[0][0] * m[1][2] * m[2][1] * m[3][3] -
+		m[0][1] * m[1][0] * m[2][2] * m[3][3] + m[0][0] * m[1][1] * m[2][2] * m[3][3];
+}
+
+template <class T>
+matrix<T, 2, 2> inverse(const matrix<T, 2, 2>& m)
+{
+	static_assert(std::is_floating_point<T>::value);
+
+	const T rd(T(1) / determinant(m));
+	return
+		{{
+			{ m[1][1] * rd, -m[0][1] * rd},
+			{-m[1][0] * rd,  m[0][0] * rd}
+		}};
+}
+
+template <class T>
+matrix<T, 3, 3> inverse(const matrix<T, 3, 3>& m)
+{
+	static_assert(std::is_floating_point<T>::value);
+
+	const T rd(T(1) / determinant(m));
+
+	return
+		{{
+			{
+				(m[1][1] * m[2][2] - m[1][2] * m[2][1]) * rd,
+				(m[0][2] * m[2][1] - m[0][1] * m[2][2]) * rd,
+				(m[0][1] * m[1][2] - m[0][2] * m[1][1]) * rd
+			},
+
+			{
+				(m[1][2] * m[2][0] - m[1][0] * m[2][2]) * rd,
+				(m[0][0] * m[2][2] - m[0][2] * m[2][0]) * rd,
+				(m[0][2] * m[1][0] - m[0][0] * m[1][2]) * rd
+			},
+
+			{
+				(m[1][0] * m[2][1] - m[1][1] * m[2][0]) * rd,
+				(m[0][1] * m[2][0] - m[0][0] * m[2][1]) * rd,
+				(m[0][0] * m[1][1] - m[0][1] * m[1][0]) * rd
+			}
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> inverse(const matrix<T, 4, 4>& m)
+{
+	static_assert(std::is_floating_point<T>::value);
+
+	const T rd(T(1) / determinant(m));
+
+	return
+		{{
+			{
+				(m[1][2] * m[2][3] * m[3][1] - m[1][3] * m[2][2] * m[3][1] + m[1][3] * m[2][1] * m[3][2] - m[1][1] * m[2][3] * m[3][2] - m[1][2] * m[2][1] * m[3][3] + m[1][1] * m[2][2] * m[3][3]) * rd,
+				(m[0][3] * m[2][2] * m[3][1] - m[0][2] * m[2][3] * m[3][1] - m[0][3] * m[2][1] * m[3][2] + m[0][1] * m[2][3] * m[3][2] + m[0][2] * m[2][1] * m[3][3] - m[0][1] * m[2][2] * m[3][3]) * rd,
+				(m[0][2] * m[1][3] * m[3][1] - m[0][3] * m[1][2] * m[3][1] + m[0][3] * m[1][1] * m[3][2] - m[0][1] * m[1][3] * m[3][2] - m[0][2] * m[1][1] * m[3][3] + m[0][1] * m[1][2] * m[3][3]) * rd,
+				(m[0][3] * m[1][2] * m[2][1] - m[0][2] * m[1][3] * m[2][1] - m[0][3] * m[1][1] * m[2][2] + m[0][1] * m[1][3] * m[2][2] + m[0][2] * m[1][1] * m[2][3] - m[0][1] * m[1][2] * m[2][3]) * rd
+			},
+
+			{
+				(m[1][3] * m[2][2] * m[3][0] - m[1][2] * m[2][3] * m[3][0] - m[1][3] * m[2][0] * m[3][2] + m[1][0] * m[2][3] * m[3][2] + m[1][2] * m[2][0] * m[3][3] - m[1][0] * m[2][2] * m[3][3]) * rd,
+				(m[0][2] * m[2][3] * m[3][0] - m[0][3] * m[2][2] * m[3][0] + m[0][3] * m[2][0] * m[3][2] - m[0][0] * m[2][3] * m[3][2] - m[0][2] * m[2][0] * m[3][3] + m[0][0] * m[2][2] * m[3][3]) * rd,
+				(m[0][3] * m[1][2] * m[3][0] - m[0][2] * m[1][3] * m[3][0] - m[0][3] * m[1][0] * m[3][2] + m[0][0] * m[1][3] * m[3][2] + m[0][2] * m[1][0] * m[3][3] - m[0][0] * m[1][2] * m[3][3]) * rd,
+				(m[0][2] * m[1][3] * m[2][0] - m[0][3] * m[1][2] * m[2][0] + m[0][3] * m[1][0] * m[2][2] - m[0][0] * m[1][3] * m[2][2] - m[0][2] * m[1][0] * m[2][3] + m[0][0] * m[1][2] * m[2][3]) * rd
+			},
+
+			{
+				(m[1][1] * m[2][3] * m[3][0] - m[1][3] * m[2][1] * m[3][0] + m[1][3] * m[2][0] * m[3][1] - m[1][0] * m[2][3] * m[3][1] - m[1][1] * m[2][0] * m[3][3] + m[1][0] * m[2][1] * m[3][3]) * rd,
+				(m[0][3] * m[2][1] * m[3][0] - m[0][1] * m[2][3] * m[3][0] - m[0][3] * m[2][0] * m[3][1] + m[0][0] * m[2][3] * m[3][1] + m[0][1] * m[2][0] * m[3][3] - m[0][0] * m[2][1] * m[3][3]) * rd,
+				(m[0][1] * m[1][3] * m[3][0] - m[0][3] * m[1][1] * m[3][0] + m[0][3] * m[1][0] * m[3][1] - m[0][0] * m[1][3] * m[3][1] - m[0][1] * m[1][0] * m[3][3] + m[0][0] * m[1][1] * m[3][3]) * rd,
+				(m[0][3] * m[1][1] * m[2][0] - m[0][1] * m[1][3] * m[2][0] - m[0][3] * m[1][0] * m[2][1] + m[0][0] * m[1][3] * m[2][1] + m[0][1] * m[1][0] * m[2][3] - m[0][0] * m[1][1] * m[2][3]) * rd
+			},
+
+			{
+				(m[1][2] * m[2][1] * m[3][0] - m[1][1] * m[2][2] * m[3][0] - m[1][2] * m[2][0] * m[3][1] + m[1][0] * m[2][2] * m[3][1] + m[1][1] * m[2][0] * m[3][2] - m[1][0] * m[2][1] * m[3][2]) * rd,
+				(m[0][1] * m[2][2] * m[3][0] - m[0][2] * m[2][1] * m[3][0] + m[0][2] * m[2][0] * m[3][1] - m[0][0] * m[2][2] * m[3][1] - m[0][1] * m[2][0] * m[3][2] + m[0][0] * m[2][1] * m[3][2]) * rd,
+				(m[0][2] * m[1][1] * m[3][0] - m[0][1] * m[1][2] * m[3][0] - m[0][2] * m[1][0] * m[3][1] + m[0][0] * m[1][2] * m[3][1] + m[0][1] * m[1][0] * m[3][2] - m[0][0] * m[1][1] * m[3][2]) * rd,
+				(m[0][1] * m[1][2] * m[2][0] - m[0][2] * m[1][1] * m[2][0] + m[0][2] * m[1][0] * m[2][1] - m[0][0] * m[1][2] * m[2][1] - m[0][1] * m[1][0] * m[2][2] + m[0][0] * m[1][1] * m[2][2]) * rd
+			}
+		}};
+}
+
+template <class T>
+matrix<T, 2, 2> componentwise_mul(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y)
+{
+	return
+		{{
+			 {x[0][0] * y[0][0], x[0][1] * y[0][1]},
+			 {x[1][0] * y[1][0], x[1][1] * y[1][1]}
+		}};
+}
+
+template <class T>
+matrix<T, 3, 3> componentwise_mul(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y)
+{
+	return
+		{{
+			 {x[0][0] * y[0][0], x[0][1] * y[0][1], x[0][2] * y[0][2]},
+			 {x[1][0] * y[1][0], x[1][1] * y[1][1], x[1][2] * y[1][2]},
+			 {x[2][0] * y[2][0], x[2][1] * y[2][1], x[2][2] * y[2][2]}
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> componentwise_mul(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y)
+{
+	return
+		{{
+			 {x[0][0] * y[0][0], x[0][1] * y[0][1], x[0][2] * y[0][2], x[0][3] * y[0][3]},
+			 {x[1][0] * y[1][0], x[1][1] * y[1][1], x[1][2] * y[1][2], x[1][3] * y[1][3]},
+			 {x[2][0] * y[2][0], x[2][1] * y[2][1], x[2][2] * y[2][2], x[2][3] * y[2][3]},
+			 {x[3][0] * y[3][0], x[3][1] * y[3][1], x[3][2] * y[3][2], x[3][3] * y[3][3]}
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> look_at(const vector<T, 3>& position, const vector<T, 3>& target, vector<T, 3> up)
+{
+	vector<T, 3> forward = normalize(sub(target, position));
+	vector<T, 3> right = normalize(cross(forward, up));
+	up = cross(right, forward);
+
+	matrix<T, 4, 4> m = 
+		{{
+			{right[0], up[0], -forward[0], T(0)},
+			{right[1], up[1], -forward[1], T(0)},
+			{right[2], up[2], -forward[2], T(0)},
+			{T(0), T(0), T(0), T(1)}
+		}};
+	
+	return translate(m, negate(position));
+}
+
+template <class T>
+matrix<T, 2, 2> mul(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y)
+{
+	return
+		{{
+			x[0] * y[0][0] + x[1] * y[0][1],
+			x[0] * y[1][0] + x[1] * y[1][1]
+		}};
+}
+
+template <class T>
+matrix<T, 3, 3> mul(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y)
+{
+	return
+		{{
+			x[0] * y[0][0] + x[1] * y[0][1] + x[2] * y[0][2],
+			x[0] * y[1][0] + x[1] * y[1][1] + x[2] * y[1][2],
+			x[0] * y[2][0] + x[1] * y[2][1] + x[2] * y[2][2]
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> mul(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y)
+{
+	return
+		{{
+			x[0] * y[0][0] + x[1] * y[0][1] + x[2] * y[0][2] + x[3] * y[0][3],
+			x[0] * y[1][0] + x[1] * y[1][1] + x[2] * y[1][2] + x[3] * y[1][3],
+			x[0] * y[2][0] + x[1] * y[2][1] + x[2] * y[2][2] + x[3] * y[2][3],
+			x[0] * y[3][0] + x[1] * y[3][1] + x[2] * y[3][2] + x[3] * y[3][3]
+		}};
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+inline matrix<T, N, M> mul(const matrix<T, N, M>& m, T s, std::index_sequence<I...>)
+{
+	return {{(m[I] * s)...}};
+}
+
+template <class T, std::size_t N, std::size_t M>
+inline matrix<T, N, M> mul(const matrix<T, N, M>& m, T s)
+{
+	return mul(m, s, std::make_index_sequence<N>{});
+}
+
+template <class T>
+vector<T, 2> mul(const matrix<T, 2, 2>& m, const vector<T, 2>& v)
+{
+	return
+		{
+			m[0][0] * v[0] + m[1][0] * v[1],
+			m[0][1] * v[0] + m[1][1] * v[1]
+		};
+}
+
+template <class T>
+vector<T, 3> mul(const matrix<T, 3, 3>& m, const vector<T, 3>& v)
+{
+	return
+		{
+			m[0][0] * v[0] + m[1][0] * v[1] + m[2][0] * v[2],
+			m[0][1] * v[0] + m[1][1] * v[1] + m[2][1] * v[2],
+			m[0][2] * v[0] + m[1][2] * v[1] + m[2][2] * v[2]
+		};
+}
+
+template <class T>
+vector<T, 4> mul(const matrix<T, 4, 4>& m, const vector<T, 4>& v)
+{
+	return
+		{
+			m[0][0] * v[0] + m[1][0] * v[1] + m[2][0] * v[2] + m[3][0] * v[3],
+			m[0][1] * v[0] + m[1][1] * v[1] + m[2][1] * v[2] + m[3][1] * v[3],
+			m[0][2] * v[0] + m[1][2] * v[1] + m[2][2] * v[2] + m[3][2] * v[3],
+			m[0][3] * v[0] + m[1][3] * v[1] + m[2][3] * v[2] + m[3][3] * v[3]
+		};
+}
+
+template <class T>
+matrix<T, 4, 4> ortho(T left, T right, T bottom, T top, T z_near, T z_far)
+{
+	return
+		{{
+			 {T(2) / (right - left), T(0), T(0), T(0)},
+			 {T(0), T(2) / (top - bottom), T(0), T(0)},
+			 {T(0), T(0), T(-2) / (z_far - z_near), T(0)},
+			 {-((right + left) / (right - left)), -((top + bottom) / (top - bottom)), -((z_far + z_near) / (z_far - z_near)), T(1)}
+		}};
+}
+
+template <class T>
+matrix<T, 2, 2> outer_product(const vector<T, 2>& c, const vector<T, 2>& r)
+{
+	return
+		{{
+			 {c[0] * r[0], c[1] * r[0]},
+			 {c[0] * r[1], c[1] * r[1]}
+		}};
+}
+
+template <class T>
+matrix<T, 3, 3> outer_product(const vector<T, 3>& c, const vector<T, 3>& r)
+{
+	return
+		{{
+			 {c[0] * r[0], c[1] * r[0], c[2] * r[0]},
+			 {c[0] * r[1], c[1] * r[1], c[2] * r[1]},
+			 {c[0] * r[2], c[1] * r[2], c[2] * r[2]}
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> outer_product(const vector<T, 4>& c, const vector<T, 4> r)
+{
+	return
+		{{
+			 {c[0] * r[0], c[1] * r[0], c[2] * r[0], c[3] * r[0]},
+			 {c[0] * r[1], c[1] * r[1], c[2] * r[1], c[3] * r[1]},
+			 {c[0] * r[2], c[1] * r[2], c[2] * r[2], c[3] * r[2]},
+			 {c[0] * r[3], c[1] * r[3], c[2] * r[3], c[3] * r[3]}
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> perspective(T vertical_fov, T aspect_ratio, T z_near, T z_far)
+{
+	T half_fov = vertical_fov * T(0.5);
+	T f = std::cos(half_fov) / std::sin(half_fov);
+
+	return
+		{{
+			 {f / aspect_ratio, T(0), T(0), T(0)},
+			 {T(0), f, T(0), T(0)},
+			 {T(0), T(0), (z_far + z_near) / (z_near - z_far), T(-1)},
+			 {T(0), T(0), (T(2) * z_near * z_far) / (z_near - z_far), T(0)}
+		}};
+}
+
+template <std::size_t N1, std::size_t M1, class T, std::size_t N0, std::size_t M0>
+matrix<T, N1, M1> resize(const matrix<T, N0, M0>& m)
+{
+	matrix<T, N1, M1> resized;
+
+	for (std::size_t i = 0; i < N1; ++i)
+	{
+		for (std::size_t j = 0; j < M1; ++j)
+		{
+			resized[i][j] = (i < N0 && j < M0) ? m[i][j] : ((i == j) ? T(1) : T(0));
+		}
+	}
+
+	return resized;
+}
+
+template <class T>
+matrix<T, 4, 4> rotate(const matrix<T, 4, 4>& m, T angle, const vector<T, 3>& axis)
+{
+	const T c = std::cos(angle);
+	const T s = std::sin(angle);
+	const vector<T, 3> temp = mul(axis, T(1) - c);
+
+	matrix<T, 4, 4> rotation;
+	rotation[0][0] = axis[0] * temp[0] + c
+	rotation[0][1] = axis[1] * temp[0] + axis[2] * s;
+	rotation[0][2] = axis[2] * temp[0] - axis[1] * s;
+	rotation[0][3] = T(0);
+	rotation[1][0] = axis[0] * temp[1] - axis[2] * s
+	rotation[1][1] = axis[1] * temp[1] + c;
+	rotation[1][2] = axis[2] * temp[1] + axis[0] * s;
+	rotation[1][3] = T(0);
+	rotation[2][0] = axis[0] * temp[2] + axis[1] * s;
+	rotation[2][1] = axis[1] * temp[2] - axis[0] * s;
+	rotation[2][2] = axis[2] * temp[2] + c
+	rotation[2][3] = T(0);
+	rotation[3][0] = T(0);
+	rotation[3][1] = T(0);
+	rotation[3][2] = T(0);
+	rotation[3][3] = T(1);
+
+	return mul(m, rotation);
+}
+
+template <class T>
+matrix<T, 4, 4> scale(const matrix<T, 4, 4>& m, const vector<T, 3>& v)
+{
+	return mul(m, matrix<T, 4, 4>
+		{{
+			{v[0], T(0), T(0), T(0)},
+			{T(0), v[1], T(0), T(0)},
+			{T(0), T(0), v[2], T(0)},
+			{T(0), T(0), T(0), T(1)}
+		}});
+}
+
+template <class T>
+matrix<T, 2, 2> sub(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y)
+{
+	return
+		{{
+			x[0] - y[0],
+			x[1] - y[1]
+		}};
+}
+
+template <class T>
+matrix<T, 3, 3> sub(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y)
+{
+	return
+		{{
+			x[0] - y[0],
+			x[1] - y[1],
+			x[2] - y[2]
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> sub(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y)
+{
+	return
+		{{
+			x[0] - y[0],
+			x[1] - y[1],
+			x[2] - y[2],
+			x[3] - y[3]
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> translate(const matrix<T, 4, 4>& m, const vector<T, 3>& v)
+{
+	return mul(m, matrix<T, 4, 4>
+		{{
+			{T(1), T(0), T(0), T(0)},
+			{T(0), T(1), T(0), T(0)},
+			{T(0), T(0), T(1), T(0)},
+			{v[0], v[1], v[2], T(1)}
+		}});
+}
+
+template <class T>
+matrix<T, 2, 2> transpose(const matrix<T, 2, 2>& m)
+{
+
+	return
+		{{
+			{
+				m[0][0], m[1][0]
+			},
+
+			{
+				m[0][1], m[1][1]
+			}
+		}};
+}
+
+template <class T>
+matrix<T, 3, 3> transpose(const matrix<T, 3, 3>& m)
+{
+
+	return
+		{{
+			{
+				m[0][0], m[1][0], m[2][0]
+			},
+
+			{
+				m[0][1], m[1][1], m[2][1]
+			},
+
+			{
+				m[0][2], m[1][2], m[2][2]
+			}
+		}};
+}
+
+template <class T>
+matrix<T, 4, 4> transpose(const matrix<T, 4, 4>& m)
+{
+
+	return
+		{{
+			{
+				m[0][0], m[1][0], m[2][0], m[3][0]
+			},
+
+			{
+				m[0][1], m[1][1], m[2][1], m[3][1]
+			},
+
+			{
+				m[0][2], m[1][2], m[2][2], m[3][2]
+			},
+
+			{
+				m[0][3], m[1][3], m[2][3], m[3][3]
+			}
+		}};
+}
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_MATRIX_FUNCTIONS_HPP

src/math/matrix-operators.hpp

@@ -0,0 +1,100 @@
+/*
+ * Copyright (C) 2020  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_MATRIX_OPERATORS_HPP
+#define ANTKEEPER_MATH_MATRIX_OPERATORS_HPP
+
+#include "math/matrix-type.hpp"
+#include "math/matrix-functions.hpp"
+
+namespace math {
+namespace matrix_operators {
+
+/// @addtogroup matrix
+/// @{
+
+/// @copydoc add(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T, std::size_t N, std::size_t M>
+matrix<T, N, M> operator+(const matrix<T, N, M>& x, const matrix<T, N, M>& y);
+
+/// @copydoc mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T, std::size_t N, std::size_t M>
+matrix<T, N, M> operator*(const matrix<T, N, M>& x, const matrix<T, N, M>& y);
+
+/// @copydoc mul(const matrix<T, N, M>&, T)
+template <class T, std::size_t N, std::size_t M>
+matrix<T, N, M> operator*(const matrix<T, N, M>& m, T s);
+
+/// @copydoc mul(const matrix<T, N, M>&, T)
+template <class T, std::size_t N, std::size_t M>
+matrix<T, N, M> operator*(T s, const matrix<T, N, M>& m);
+
+/// @copydoc mul(const matrix<T, 2, 2>&, const vector<T, 2>&)
+template <class T, std::size_t N>
+vector<T, N> operator*(const matrix<T, N, N>& m, const vector<T, N>& v);
+
+/// @copydoc sub(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
+template <class T, std::size_t N, std::size_t M>
+matrix<T, N, M> operator-(const matrix<T, N, M>& x, const matrix<T, N, M>& y);
+
+template <class T, std::size_t N, std::size_t M>
+inline matrix<T, N, M> operator+(const matrix<T, N, M>& x, const matrix<T, N, M>& y)
+{
+	return add(x, y);
+}
+
+template <class T, std::size_t N, std::size_t M>
+inline matrix<T, N, M> operator*(const matrix<T, N, M>& x, const matrix<T, N, M>& y)
+{
+	return mul(x, y);
+}
+
+template <class T, std::size_t N, std::size_t M>
+inline matrix<T, N, M> operator*(const matrix<T, N, M>& m, T s)
+{
+	return mul(m, s);
+}
+
+template <class T, std::size_t N, std::size_t M>
+inline matrix<T, N, M> operator*(T s, const matrix<T, N, M>& m)
+{
+	return mul(m, s);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> operator*(const matrix<T, N, N>& m, const vector<T, N>& v)
+{
+	return mul(m, v);
+}
+
+template <class T, std::size_t N, std::size_t M>
+inline matrix<T, N, M> operator-(const matrix<T, N, M>& x, const matrix<T, N, M>& y)
+{
+	return sub(x, y);
+}
+
+/// @}
+
+} // namespace matrix_operators
+} // namespace math
+
+// Bring matrix operators into global namespace
+using namespace math::matrix_operators;
+
+#endif // ANTKEEPER_MATH_MATRIX_OPERATORS_HPP

src/math/matrix-type.hpp

@@ -0,0 +1,54 @@
+/*
+ * Copyright (C) 2020  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_MATRIX_TYPE_HPP
+#define ANTKEEPER_MATH_MATRIX_TYPE_HPP
+
+#include "math/vector-type.hpp"
+#include <cstddef>
+
+namespace math {
+
+/// @addtogroup matrix
+/// @{
+
+/**
+ * An NxM matrix.
+ *
+ * @tparam T Matrix element type.
+ * @tparam N Number of columns.
+ * @tparam M Number of rows.
+ */
+template <typename T, std::size_t N, std::size_t M>
+struct matrix
+{
+	typedef T element_type;
+	typedef vector<element_type, M> row_type;
+	row_type columns[N];
+
+	inline constexpr row_type& operator[](std::size_t i) noexcept { return columns[i]; }
+	inline constexpr const row_type& operator[](std::size_t i) const noexcept { return columns[i]; }
+};
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_MATRIX_TYPE_HPP
+

src/math/operators.hpp

@@ -0,0 +1,29 @@
+/*
+ * Copyright (C) 2020  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_OPERATORS_HPP
+#define ANTKEEPER_MATH_OPERATORS_HPP
+
+#include "math/vector-operators.hpp"
+#include "math/matrix-operators.hpp"
+#include "math/quaternion-operators.hpp"
+#include "math/transform-operators.hpp"
+#include "math/stream-operators.hpp"
+
+#endif // ANTKEEPER_MATH_OPERATORS_HPP

src/math/quaternion-functions.hpp

@@ -0,0 +1,443 @@
+/*
+ * Copyright (C) 2020  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/matrix-type.hpp"
+#include "math/quaternion-type.hpp"
+#include "math/vector-type.hpp"
+#include "math/vector-functions.hpp"
+#include <cmath>
+
+namespace math {
+
+/// @addtogroup quaternion
+/// @{
+
+/**
+ * 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);
+
+/**
+ * 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);
+
+/**
+ * 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);
+
+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
+		{
+			x.w * (T(1) - a) + y.w * a,
+			x.x * (T(1) - a) + y.x * a,
+			x.y * (T(1) - a) + y.y * a,
+			x.z * (T(1) - a) + y.z * a
+		};
+}
+
+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[0], up[0], -forward[0]},
+			{right[1], up[1], -forward[1]},
+			{right[2], up[2], -forward[2]}
+		}};
+
+	// Convert to quaternion
+	return normalize(quaternion_cast(m));
+}
+
+template <class T>
+matrix<T, 3, 3> matrix_cast(const quaternion<T>& q)
+{
+	T wx = q.w * q.x;
+	T wy = q.w * q.y;
+	T wz = q.w * q.z;
+	T xx = q.x * q.x;
+	T xy = q.x * q.y;
+	T xz = q.x * q.z;
+	T yy = q.y * q.y;
+	T yz = q.y * q.z;
+	T zz = q.z * q.z;
+
+	return
+		{{
+			{T(1) - (yy + zz) * T(2), (xy + wz) * T(2), (xz - wy) * T(2)},
+			{(xy - wz) * T(2), T(1) - (xx + zz) * T(2), (yz + wx) * T(2)},
+			{(xz + wy) * T(2), (yz - wx) * 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 T r = q.w;             // Real part
+	const vector<T, 3>& i = reinterpret_cast<const vector<T, 3>&>(q.x); // Imaginary part
+	return i * dot(i, v) * T(2) + v * (r * r - dot(i, i)) + cross(i, v) * r * 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>
+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>
+quaternion<T> quaternion_cast(const matrix<T, 3, 3>& m)
+{
+	T r;
+	vector<T, 3> i;
+
+	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));
+		r = T(0.25) / s;
+		i =
+			{
+				(m[2][1] - m[1][2]) * s,
+				(m[0][2] - m[2][0]) * s,
+				(m[1][0] - m[0][1]) * 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]);
+			r = (m[2][1] - m[1][2]) / s;
+			i =
+				{
+					T(0.25) * s,
+					(m[0][1] + m[1][0]) / s,
+					(m[0][2] + m[2][0]) / 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]);
+			r = (m[0][2] - m[2][0]) / s;
+			i =
+				{
+					(m[0][1] + m[1][0]) / s,
+					T(0.25) * s,
+					(m[1][2] + m[2][1]) / s
+				};
+		}
+		else
+		{
+			T s = T(2) * std::sqrt(T(1) + m[2][2] - m[0][0] - m[1][1]);
+			r = (m[1][0] - m[0][1]) / s;
+			i = 
+				{
+					(m[0][2] + m[2][0]) / s,
+					(m[1][2] + m[2][1]) / s,
+					T(0.25) * s
+				};
+		}
+	}
+
+	return {r, i.x, i.y, i.z};
+}
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_QUATERNION_FUNCTIONS_HPP
+

src/math/quaternion-operators.hpp

@@ -0,0 +1,111 @@
+/*
+ * Copyright (C) 2020  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"
+
+namespace math {
+namespace quaternion_operators {
+
+/// @addtogroup quaternion
+/// @{
+
+/// @copydoc add(const quaternion<T>&, const quaternion<T>&)
+template <class T>
+quaternion<T> operator+(const quaternion<T>& x, const quaternion<T>& y);
+
+/// @copydoc div(const quaternion<T>&, T)
+template <class T>
+quaternion<T> operator/(const quaternion<T>& q, T s);
+
+/// @copydoc mul(const quaternion<T>&, const quaternion<T>&)
+template <class T>
+quaternion<T> operator*(const quaternion<T>& x, const quaternion<T>& y);
+
+/// @copydoc mul(const quaternion<T>&, T)
+template <class T>
+quaternion<T> operator*(const quaternion<T>& q, T s);
+
+/// @copydoc mul(const quaternion<T>&, const vector<T, 3>&)
+template <class T>
+vector<T, 3> operator*(const quaternion<T>& q, const vector<T, 3>& v);
+
+/// @copydoc sub(const quaternion<T>&, const quaternion<T>&)
+template <class T>
+quaternion<T> operator-(const quaternion<T>& x, const quaternion<T>& y);
+
+/// @copydoc negate(const quaternion<T>&)
+template <class T, std::size_t N>
+quaternion<T> operator-(const quaternion<T>& x);
+
+template <class T>
+inline quaternion<T> operator+(const quaternion<T>& x, const quaternion<T>& y)
+{
+	return add(x, y);
+}
+
+template <class T>
+inline quaternion<T> operator/(const quaternion<T>& q, T s)
+{
+	return div(q, s);
+}
+
+template <class T>
+inline quaternion<T> operator*(const quaternion<T>& x, const quaternion<T>& y)
+{
+	return mul(x, y);
+}
+
+template <class T>
+inline quaternion<T> operator*(const quaternion<T>& q, T s)
+{
+	return mul(q, s);
+}
+
+template <class T>
+inline vector<T, 3> operator*(const quaternion<T>& q, const vector<T, 3>& v)
+{
+	return mul(q, v);
+}
+
+template <class T>
+inline quaternion<T> operator-(const quaternion<T>& x, const quaternion<T>& y)
+{
+	return sub(x, y);
+}
+
+template <class T, std::size_t N>
+inline quaternion<T> operator-(const quaternion<T>& x)
+{
+	return negate(x);
+}
+
+/// @}
+
+} // namespace quaternion_operators
+} // namespace math
+
+/// Bring quaternion operators into global namespace
+using namespace math::quaternion_operators;
+
+#endif // ANTKEEPER_MATH_QUATERNION_OPERATORS_HPP
+

src/math/quaternion-type.hpp

@@ -0,0 +1,48 @@
+/*
+ * Copyright (C) 2020  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
+
+namespace math {
+
+/// @addtogroup quaternion
+/// @{
+
+/**
+ * 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;
+};
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_QUATERNION_TYPE_HPP
+

src/math/random.hpp

@@ -0,0 +1,55 @@
+/*
+ * Copyright (C) 2020  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_RANDOM_HPP
+#define ANTKEEPER_MATH_RANDOM_HPP
+
+#include <cstdlib>
+#include <type_traits>
+
+namespace math {
+
+/// @addtogroup utility
+/// @{
+
+/**
+ * Generates a pseudo-random floating point number on `[start, end)` using std::rand().
+ *
+ * @warning Don't forget to seed with std::srand() before using!
+ *
+ * @param start Start of the range (inclusive).
+ * @param end End of the range (exclusive).
+ * @return Pseudo-random floating point number.
+ */
+template <typename T = float>
+T random(T start, T end);
+
+template <typename T>
+inline T random(T start, T end)
+{
+	static_assert(std::is_floating_point<T>::value);
+	constexpr T rand_max_inverse = T(1) / static_cast<T>(RAND_MAX);
+    return static_cast<T>(std::rand()) * rand_max_inverse * (end - start) + start;
+}
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_RANDOM_HPP

src/math/stream-operators.hpp

@@ -0,0 +1,114 @@
+/*
+ * Copyright (C) 2020  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/vector-type.hpp"
+#include "math/matrix-type.hpp"
+#include "math/quaternion-type.hpp"
+#include <ostream>
+
+namespace math {
+namespace stream_operators {
+
+/// @addtogroup io
+/// @{
+
+/**
+ * Writes the elements of a vector to an output stream, with each element delimeted by a space.
+ *
+ * @param os Output stream.
+ * @param v Vector.
+ * @return Output stream.
+ */
+template <class T, std::size_t N>
+std::ostream& operator<<(std::ostream& os, const vector<T, N>& v);
+
+/**
+ * Writes the elements of a matrix to an output stream, with each element delimeted by a space.
+ *
+ * @param os Output stream.
+ * @param m Matrix.
+ * @return Output stream.
+ */
+template <class T, std::size_t N, std::size_t M>
+std::ostream& operator<<(std::ostream& os, const matrix<T, N, M>& m);
+
+/**
+ * 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 quaternion<T>& q);
+
+template <class T, std::size_t N>
+std::ostream& operator<<(std::ostream& os, const vector<T, N>& v)
+{
+	for (std::size_t i = 0; i < N; ++i)
+	{
+		os << v[i];
+
+		if (i < N - 1)
+		{
+			os << ' ';
+		}
+	}
+
+	return os;
+}
+
+template <class T, std::size_t N, std::size_t M>
+std::ostream& operator<<(std::ostream& os, const matrix<T, N, M>& m)
+{
+	for (std::size_t i = 0; i < N; ++i)
+	{
+		for (std::size_t j = 0; j < M; ++j)
+		{
+			os << m[i][j];
+
+			if (i < N - 1 || j < M - 1)
+			{
+				os << ' ';
+			}
+		}
+	}
+
+	return os;
+}
+
+template <class T>
+std::ostream& operator<<(std::ostream& os, const quaternion<T>& q)
+{
+	os << std::get<0>(q) << ' ' << std::get<1>(q)[0] << ' ' << std::get<1>(q)[1] << ' ' << std::get<1>(q)[2];
+	return os;
+}
+
+/// @}
+
+} // namespace stream_operators
+} // namespace math
+
+// Bring stream operators into global namespace
+using namespace math::stream_operators;
+
+#endif // ANTKEEPER_MATH_STREAM_OPERATORS_HPP

src/math/transform-functions.hpp

@@ -0,0 +1,110 @@
+/*
+ * Copyright (C) 2020  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_TRANSFORM_FUNCTIONS_HPP
+#define ANTKEEPER_MATH_TRANSFORM_FUNCTIONS_HPP
+
+#include "math/transform-type.hpp"
+#include "math/vector-functions.hpp"
+#include "math/matrix-functions.hpp"
+#include "math/quaternion-functions.hpp"
+
+namespace math {
+
+/// @addtogroup transform
+/// @{
+
+/**
+ * Calculates the inverse of a transform.
+ *
+ * @param t Transform of which to take the inverse.
+ */
+template <class T>
+transform<T> inverse(const transform<T>& t);
+
+/**
+ * Converts a transform to a transformation matrix.
+ *
+ * @param t Transform.
+ * @return Matrix representing the transformation described by `t`.
+ */
+template <class T>
+matrix<T, 4, 4> matrix_cast(const transform<T>& t);
+
+/**
+ * Multiplies two transforms.
+ *
+ * @param x First transform.
+ * @param y Second transform.
+ * @return Product of the two transforms.
+ */
+template <class T>
+transform<T> mul(const transform<T>& x, const transform<T>& y);
+
+/**
+ * Multiplies a vector by a transform.
+ *
+ * @param t Transform.
+ * @param v Vector.
+ * @return Product of the transform and vector.
+ */
+template <class T>
+vector<T, 3> mul(const transform<T>& t, const vector<T, 3>& v);
+
+template <class T>
+transform<T> inverse(const transform<T>& t)
+{
+	transform<T> inverse_t;
+	inverse_t.scale = {T(1) / t.scale.x, T(1) / t.scale.y, T(1) / t.scale.z};
+	inverse_t.rotation = conjugate(t.rotation);
+	inverse_t.translation = negate(mul(inverse_t.rotation, t.translation));
+	return inverse_t;
+}
+
+template <class T>
+matrix<T, 4, 4> matrix_cast(const transform<T>& t)
+{
+	matrix<T, 4, 4> transformation = resize<4, 4>(matrix_cast(t.rotation));
+	transformation[3] = {t.translation[0], t.translation[1], t.translation[2], T(1)};
+	return scale(transformation, t.scale);
+}
+
+template <class T>
+transform<T> mul(const transform<T>& x, const transform<T>& y)
+{
+	return
+		{
+			mul(x, y.translation),
+			normalize(mul(x.rotation, y.rotation)),
+			x.scale * y.scale
+		};
+}
+
+template <class T>
+vector<T, 3> mul(const transform<T>& t, const vector<T, 3>& v)
+{
+	return t.translation + (t.rotation * (v * t.scale));
+}
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_TRANSFORM_FUNCTIONS_HPP
+

src/math/transform-operators.hpp

@@ -0,0 +1,77 @@
+/*
+ * Copyright (C) 2020  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_TRANSFORM_OPERATORS_HPP
+#define ANTKEEPER_MATH_TRANSFORM_OPERATORS_HPP
+
+#include "math/transform-type.hpp"
+#include "math/transform-functions.hpp"
+
+namespace math {
+namespace transform_operators {
+
+/// @addtogroup transform
+/// @{
+
+/// @copydoc mul(const transform<T>&, const transform<T>&)
+template <class T>
+transform<T> operator*(const transform<T>& x, const transform<T>& y);
+
+/// @copydoc mul(const transform<T>&, const vector<T, 3>&)
+template <class T>
+vector<T, 3> operator*(const transform<T>& t, const vector<T, 3>& v);
+
+/**
+ * Multiplies two transforms and stores the result in the first transform.
+ *
+ * @param x First transform.
+ * @param y Second transform.
+ * @return Reference to the first transform.
+ */
+template <class T>
+transform<T>& operator*=(transform<T>& x, const transform<T>& y);
+
+template <class T>
+inline transform<T> operator*(const transform<T>& x, const transform<T>& y)
+{
+	return mul(x, y);
+}
+
+template <class T>
+inline vector<T, 3> operator*(const transform<T>& t, const vector<T, 3>& v)
+{
+	return mul(t, v);
+}
+
+template <class T>
+inline transform<T>& operator*=(transform<T>& x, const vector<T, 3>& y)
+{
+	return (x = x * y);
+}
+
+/// @}
+
+} // namespace transform_operators
+} // namespace math
+
+// Bring transform operators into global namespace
+using namespace math::transform_operators;
+
+#endif // ANTKEEPER_MATH_TRANSFORM_OPERATORS_HPP
+

src/math/transform-type.hpp

@@ -0,0 +1,54 @@
+/*
+ * Copyright (C) 2020  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_TRANSFORM_TYPE_HPP
+#define ANTKEEPER_MATH_TRANSFORM_TYPE_HPP
+
+#include "math/vector-type.hpp"
+#include "math/quaternion-type.hpp"
+
+namespace math {
+
+/// @addtogroup transform
+/// @{
+
+/**
+ * Represents 3D TRS transformation.
+ *
+ * @tparam T Scalar type.
+ */
+template <class T>
+struct transform
+{
+	/// Translation vector
+	vector<T, 3> translation;
+
+	/// Rotation quaternion
+	quaternion<T> rotation;
+
+	/// Scale vector
+	vector<T, 3> scale;
+};
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_TRANSFORM_TYPE_HPP
+

src/math/vector-functions.hpp

@@ -0,0 +1,627 @@
+/*
+ * Copyright (C) 2020  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_VECTOR_FUNCTIONS_HPP
+#define ANTKEEPER_MATH_VECTOR_FUNCTIONS_HPP
+
+#include "math/vector-type.hpp"
+#include <algorithm>
+#include <cmath>
+#include <type_traits>
+#include <utility>
+
+namespace math {
+
+/// @addtogroup vector
+/// @{
+
+/**
+ * Adds two vectors.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Sum of the two vectors.
+ */
+template <class T, std::size_t N>
+vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * 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>
+bool all(const vector<bool, N>& x);
+
+/**
+ * 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>
+bool any(const vector<bool, N>& x);
+
+/**
+ * Reinterprets data as an `N`-dimensional vector of type `T`.
+ *
+ * @tparam N Number of vector dimensions.
+ * @tparam T Scalar type.
+ * @param data Data to reinterpret.
+ */
+template <std::size_t N, typename T>
+vector<T, N>& as_vector(T& data);
+
+/**
+ * Clamps the values of a vector's elements.
+ *
+ * @param x Vector to clamp.
+ * @param min_value Minimum element value.
+ * @param max_value Maximum element value.
+ * @return Clamped vector.
+ */
+template <class T, std::size_t N>
+vector<T, N> clamp(const vector<T, N>& x, T min_value, T max_value);
+
+/**
+ * Clamps the length of a vector.
+ *
+ * @param x Vector to clamp.
+ * @param max_length Maximum length.
+ * @return Length-clamped vector.
+ */
+template <class T, std::size_t N>
+vector<T, N> clamp_length(const vector<T, N>& x, T max_length);
+
+/**
+ * Calculate the cross product of two vectors.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Cross product of the two vectors.
+ */
+template <class T>
+vector<T, 3> cross(const vector<T, 3>& x, const vector<T, 3>& y);
+
+/**
+ * Calculates the distance between two points.
+ *
+ * @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>
+vector<T, N> distance(const vector<T, N>& p0, const vector<T, N>& p1);
+
+/**
+ * Calculates the squared distance between two points. The squared distance can be calculated faster than the distance because a call to `std::sqrt` is saved.
+ *
+ * @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>
+vector<T, N> distance_squared(const vector<T, N>& p0, const vector<T, N>& p1);
+
+/**
+ * Divides a vector by another vector.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Result of the division.
+ */
+template <class T, std::size_t N>
+vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Divides a vector by a scalar.
+ *
+ * @param v Vector.
+ * @param s Scalar.
+ * @return Result of the division.
+ */
+template <class T, std::size_t N>
+vector<T, N> div(const vector<T, N>& v, T s);
+
+/**
+ * Calculates the dot product of two vectors.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Dot product of the two vectors.
+ */
+template <class T, std::size_t N>
+T dot(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Compares two vectors for equality
+ *
+ * @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>
+vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Performs a component-wise greater-than comparison of two vectors.
+ *
+ * @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>
+vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Performs a component-wise greater-than or equal-to comparison of two vectors.
+ *
+ * @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>
+vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * 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>
+T length(const vector<T, N>& 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>
+T length_squared(const vector<T, N>& x);
+
+/**
+ * Performs a component-wise less-than comparison of two vectors.
+ *
+ * @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>
+vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Performs a component-wise less-than or equal-to comparison of two vectors.
+ *
+ * @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>
+vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Multiplies two vectors.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Product of the two vectors.
+ */
+template <class T, std::size_t N>
+vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Multiplies a vector by a scalar.
+ *
+ * @param v Vector.
+ * @param s Scalar.
+ * @return Product of the vector and scalar.
+ */
+template <class T, std::size_t N>
+vector<T, N> mul(const vector<T, N>& v, T s);
+
+/**
+ * Negates a vector.
+ *
+ * @param x Vector to negate.
+ * @return Negated vector.
+ */
+template <class T, std::size_t N>
+vector<T, N> negate(const vector<T, N>& x);
+
+/**
+ * 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>
+vector<T, N> normalize(const vector<T, N>& x);
+
+/**
+ * Logically inverts a boolean vector.
+ *
+ * @param x Vector to be inverted.
+ * @return Logically inverted vector.
+ */
+template <class T, std::size_t N>
+vector<bool, N> not(const vector<T, N>& x);
+
+/**
+ * Compares two vectors for inequality
+ *
+ * @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>
+vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Resizes a vector. Any new elements will be set to `0`.
+ *
+ * @param v Vector to resize.
+ * @return Resized vector.
+ */
+template <std::size_t N1, class T, std::size_t N0>
+vector<T, N1> resize(const vector<T, N0>& v);
+
+/**
+ * Subtracts a vector from another vector.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Difference between the two vectors.
+ */
+template <class T, std::size_t N>
+vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Makes an m-dimensional vector by rearranging and/or duplicating elements of an n-dimensional vector.
+ *
+ * @tparam Indices List of indices of elements in the vector `v`.
+ * @tparam T Vector component type.
+ * @tparam N Number of dimensions in vector `v`.
+ * @return Vector containing elements from vector `v` in the order specified by `Indices`. The size of the returned vector is equivalent to the number of indices in `Indices`.
+ */
+template <std::size_t... Indices, class T, std::size_t N>
+vector<T, sizeof...(Indices)> swizzle(const vector<T, N>& v);
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] + y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return add(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <std::size_t N, std::size_t... I>
+inline bool all(const vector<bool, N>& x, std::index_sequence<I...>)
+{
+	return (x[I] && ...);
+}
+
+template <std::size_t N>
+inline bool all(const vector<bool, N>& x)
+{
+	return all(x, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <std::size_t N, std::size_t... I>
+inline bool any(const vector<bool, N>& x, std::index_sequence<I...>)
+{
+	return (x[I] || ...);
+}
+
+template <std::size_t N>
+inline bool any(const vector<bool, N>& x)
+{
+	return any(x, std::make_index_sequence<N>{});
+}
+
+template <std::size_t N, typename T>
+inline vector<T, N>& as_vector(T& data)
+{
+	static_assert(std::is_pod<vector<T, N>>::value);
+	return reinterpret_cast<vector<T, N>&>(data);
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<T, N> clamp(const vector<T, N>& x, T min_value, T max_value, std::index_sequence<I...>)
+{
+	return {std::min<T>(max_value, std::max<T>(min_value, x[I]))...};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> clamp(const vector<T, N>& x, T min_value, T max_value)
+{
+	return clamp(x, min_value, max_value, std::make_index_sequence<N>{});
+}
+
+template <class T, std::size_t N>
+vector<T, N> clamp_length(const vector<T, N>& x, T max_length)
+{
+	T length2 = length_squared(x);
+	return (length2 > max_length * max_length) ? (x * max_length / std::sqrt(length2)) : x;
+}
+
+template <class T>
+inline vector<T, 3> cross(const vector<T, 3>& x, const vector<T, 3>& y)
+{
+	return
+		{
+			x[1] * y[2] - y[1] * x[2],
+			x[2] * y[0] - y[2] * x[0],
+			x[0] * y[1] - y[0] * x[1]
+		};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> distance(const vector<T, N>& p0, const vector<T, N>& p1)
+{
+	static_assert(std::is_floating_point<T>::value);
+	return length(sub(p0, p1));
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> distance_squared(const vector<T, N>& p0, const vector<T, N>& p1)
+{
+	static_assert(std::is_floating_point<T>::value);
+	return length_squared(sub(p0, p1));
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] / y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return div(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<T, N> div(const vector<T, N>& v, T s, std::index_sequence<I...>)
+{
+	return {(v[I] / s)...};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> div(const vector<T, N>& v, T s)
+{
+	return div(v, s, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline T dot(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return ((x[I] * y[I]) + ...);
+}
+
+template <class T, std::size_t N>
+inline T dot(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return dot(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] == y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return equal(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] > y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return greater_than(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] >= y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return greater_than_equal(x, y, std::make_index_sequence<N>{});
+}
+
+template <class T, std::size_t N>
+inline T length(const vector<T, N>& x)
+{
+	static_assert(std::is_floating_point<T>::value);
+	return std::sqrt(dot(x, x));
+}
+
+template <class T, std::size_t N>
+inline T length_squared(const vector<T, N>& x)
+{
+	static_assert(std::is_floating_point<T>::value);
+	return dot(x, x);
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] < y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return less_than(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] <= y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return less_than_equal(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] * y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return mul(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<T, N> mul(const vector<T, N>& v, T s, std::index_sequence<I...>)
+{
+	return {(v[I] * s)...};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> mul(const vector<T, N>& v, T s)
+{
+	return mul(v, s, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<T, N> negate(const vector<T, N>& x, std::index_sequence<I...>)
+{
+	return {(-x[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> negate(const vector<T, N>& x)
+{
+	return negate(x, std::make_index_sequence<N>{});
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> normalize(const vector<T, N>& x)
+{
+	static_assert(std::is_floating_point<T>::value);
+	return mul(x, T(1) / length(x));
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<bool, N> not(const vector<T, N>& x, std::index_sequence<I...>)
+{
+	return {!x[I]...};
+}
+
+template <class T, std::size_t N>
+inline vector<bool, N> not(const vector<T, N>& x)
+{
+	return not(x, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] != y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return not_equal(x, y, std::make_index_sequence<N>{});
+}
+
+template <std::size_t N1, class T, std::size_t N0>
+vector<T, N1> resize(const vector<T, N0>& v)
+{
+	vector<T, N1> resized;
+
+	for (std::size_t i = 0; i < N1; ++i)
+	{
+		resized[i] = (i < N0) ? v[i] : T(0);
+	}
+
+	return resized;
+}
+
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+inline vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+{
+	return {(x[I] - y[I])...};
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return sub(x, y, std::make_index_sequence<N>{});
+}
+
+template <std::size_t... Indices, class T, std::size_t N>
+inline vector<T, sizeof...(Indices)> swizzle(const vector<T, N>& v)
+{
+	return { v[Indices]... };
+}
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_VECTOR_FUNCTIONS_HPP
+

src/math/vector-operators.hpp

@@ -0,0 +1,217 @@
+/*
+ * Copyright (C) 2020  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_VECTOR_OPERATORS_HPP
+#define ANTKEEPER_MATH_VECTOR_OPERATORS_HPP
+
+#include "math/vector-type.hpp"
+#include "math/vector-functions.hpp"
+
+namespace math {
+namespace vector_operators {
+
+/// @addtogroup vector
+/// @{
+
+/// @copydoc add(const vector<T, N>&, const vector<T, N>&)
+template <class T, std::size_t N>
+vector<T, N> operator+(const vector<T, N>& x, const vector<T, N>& y);
+
+/// @copydoc div(const vector<T, N>&, const vector<T, N>&)
+template <class T, std::size_t N>
+vector<T, N> operator/(const vector<T, N>& x, const vector<T, N>& y);
+
+/// @copydoc div(const vector<T, N>&, T)
+template <class T, std::size_t N>
+vector<T, N> operator/(const vector<T, N>& v, T s);
+
+/// @copydoc mul(const vector<T, N>&, const vector<T, N>&)
+template <class T, std::size_t N>
+vector<T, N> operator*(const vector<T, N>& x, const vector<T, N>& y);
+
+/// @copydoc mul(const vector<T, N>&, T)
+template <class T, std::size_t N>
+vector<T, N> operator*(const vector<T, N>& v, T s);
+
+/// @copydoc mul(const vector<T, N>&, T)
+template <class T, std::size_t N>
+vector<T, N> operator*(T s, const vector<T, N>& v);
+
+/// @copydoc negate(const vector<T, N>&)
+template <class T, std::size_t N>
+vector<T, N> operator-(const vector<T, N>& x);
+
+/// @copydoc sub(const vector<T, N>&, const vector<T, N>&)
+template <class T, std::size_t N>
+vector<T, N> operator-(const vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Adds two vectors and stores the result in the first vector.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Reference to the first vector.
+ */
+template <class T, std::size_t N>
+vector<T, N>& operator+=(vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Subtracts two vectors and stores the result in the first vector.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Reference to the first vector.
+ */
+template <class T, std::size_t N>
+vector<T, N>& operator-=(vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Multiplies two vectors and stores the result in the first vector.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Reference to the first vector.
+ */
+template <class T, std::size_t N>
+vector<T, N>& operator*=(vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Multiplies a vector and a scalar and stores the result in the vector.
+ *
+ * @param v Vector.
+ * @param s Scalar.
+ * @return Reference to the vector.
+ */
+template <class T, std::size_t N>
+vector<T, N>& operator*=(vector<T, N>& v, T s);
+
+/**
+ * Divides the first vector by the second vector the result in the first vector.
+ *
+ * @param x First vector.
+ * @param y Second vector.
+ * @return Reference to the first vector.
+ */
+template <class T, std::size_t N>
+vector<T, N>& operator/=(vector<T, N>& x, const vector<T, N>& y);
+
+/**
+ * Divides a vector by a scalar and stores the result in the vector.
+ *
+ * @param v Vector.
+ * @param s Scalar.
+ * @return Reference to the vector.
+ */
+template <class T, std::size_t N>
+vector<T, N>& operator/=(vector<T, N>& v, T s);
+
+template <class T, std::size_t N>
+inline vector<T, N> operator+(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return add(x, y);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> operator-(const vector<T, N>& x)
+{
+	return negate(x);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> operator-(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return sub(x, y);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> operator*(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return mul(x, y);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> operator*(const vector<T, N>& v, T s)
+{
+	return mul(v, s);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> operator*(T s, const vector<T, N>& v)
+{
+	return mul(v, s);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> operator/(const vector<T, N>& x, const vector<T, N>& y)
+{
+	return div(x, y);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N> operator/(const vector<T, N>& v, T s)
+{
+	return div(v, s);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N>& operator+=(vector<T, N>& x, const vector<T, N>& y)
+{
+	return (x = x + y);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N>& operator-=(vector<T, N>& x, const vector<T, N>& y)
+{
+	return (x = x - y);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N>& operator*=(vector<T, N>& x, const vector<T, N>& y)
+{
+	return (x = x * y);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N>& operator*=(vector<T, N>& v, T s)
+{
+	return (v = v * s);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N>& operator/=(vector<T, N>& x, const vector<T, N>& y)
+{
+	return (x = x * y);
+}
+
+template <class T, std::size_t N>
+inline vector<T, N>& operator/=(vector<T, N>& v, T s)
+{
+	return (v = v / s);
+}
+
+/// @}
+
+} // namespace vector_operators
+} // namespace math
+
+// Bring vector operators into global namespace
+using namespace math::vector_operators;
+
+#endif // ANTKEEPER_MATH_VECTOR_OPERATORS_HPP
+

src/math/vector-type.hpp

@@ -0,0 +1,131 @@
+/*
+ * Copyright (C) 2020  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_VECTOR_TYPE_HPP
+#define ANTKEEPER_MATH_VECTOR_TYPE_HPP
+
+#include <array>
+#include <cstddef>
+
+namespace math {
+
+/// @addtogroup vector
+/// @{
+
+/**
+ * An `N`-dimensional Euclidean vector.
+ *
+ * @tparam T Vector component type.
+ * @tparam N Number of dimensions.
+ */
+template <typename T, std::size_t N>
+struct vector
+{
+	typedef T scalar_type;
+	typedef std::array<T, N> array_type;
+	scalar_type components[N];
+
+	inline constexpr scalar_type& operator[](std::size_t i) noexcept { return components[i]; }
+	inline constexpr const scalar_type& operator[](std::size_t i) const noexcept { return components[i]; }
+	inline constexpr operator array_type&() noexcept { return reinterpret_cast<array_type&>(components[0]); }
+	inline constexpr operator const array_type&() const noexcept { return reinterpret_cast<const array_type&>(components[0]); }
+	inline constexpr const scalar_type* data() const noexcept { return components; };
+	inline constexpr scalar_type* data() noexcept { return components; };
+	inline constexpr std::size_t size() const noexcept { return N; };
+};
+
+template <typename T>
+struct vector<T, 1>
+{
+	typedef T scalar_type;
+	typedef std::array<T, 1> array_type;
+	scalar_type x;
+
+	inline constexpr scalar_type& operator[](std::size_t i) noexcept { return *((&x) + i); }
+	inline constexpr const scalar_type& operator[](std::size_t i) const noexcept { return *((&x) + i); }
+	inline constexpr operator array_type&() noexcept { return reinterpret_cast<array_type&>(x); }
+	inline constexpr operator const array_type&() const noexcept { return reinterpret_cast<const array_type&>(x); }
+	inline constexpr const scalar_type* data() const noexcept { return &x; };
+	inline constexpr scalar_type* data() noexcept { return &x; };
+	inline constexpr std::size_t size() const noexcept { return 1; };
+};
+
+template <typename T>
+struct vector<T, 2>
+{
+	typedef T scalar_type;
+	typedef std::array<T, 2> array_type;
+	scalar_type x;
+	scalar_type y;
+
+	inline constexpr scalar_type& operator[](std::size_t i) noexcept { return *((&x) + i); }
+	inline constexpr const scalar_type& operator[](std::size_t i) const noexcept { return *((&x) + i); }
+	inline constexpr operator array_type&() noexcept { return reinterpret_cast<array_type&>(x); }
+	inline constexpr operator const array_type&() const noexcept { return reinterpret_cast<const array_type&>(x); }
+	inline constexpr const scalar_type* data() const noexcept { return &x; };
+	inline constexpr scalar_type* data() noexcept { return &x; };
+	inline constexpr std::size_t size() const noexcept { return 2; };
+};
+
+template <typename T>
+struct vector<T, 3>
+{
+	typedef T scalar_type;
+	typedef std::array<T, 3> array_type;
+	scalar_type x;
+	scalar_type y;
+	scalar_type z;
+
+	inline constexpr scalar_type& operator[](std::size_t i) noexcept { return *((&x) + i); }
+	inline constexpr const scalar_type& operator[](std::size_t i) const noexcept { return *((&x) + i); }
+	inline constexpr operator array_type&() noexcept { return reinterpret_cast<array_type&>(x); }
+	inline constexpr operator const array_type&() const noexcept { return reinterpret_cast<const array_type&>(x); }
+	inline constexpr const scalar_type* data() const noexcept { return &x; };
+	inline constexpr scalar_type* data() noexcept { return &x; };
+	inline constexpr std::size_t size() const noexcept { return 3; };
+
+};
+
+template <typename T>
+struct vector<T, 4>
+{
+	typedef T scalar_type;
+	typedef std::array<T, 4> array_type;
+	scalar_type x;
+	scalar_type y;
+	scalar_type z;
+	scalar_type w;
+
+	inline constexpr scalar_type& operator[](std::size_t i) noexcept { return *((&x) + i); }
+	inline constexpr const scalar_type& operator[](std::size_t i) const noexcept { return *((&x) + i); }
+	inline constexpr operator array_type&() noexcept { return reinterpret_cast<array_type&>(x); }
+	inline constexpr operator const array_type&() const noexcept { return reinterpret_cast<const array_type&>(x); }
+	inline constexpr const scalar_type* data() const noexcept { return &x; };
+	inline constexpr scalar_type* data() noexcept { return &x; };
+	inline constexpr std::size_t size() const noexcept { return 4; };
+};
+
+static_assert(std::is_trivial<vector<float, 4>>::value);
+
+/// @}
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_VECTOR_TYPE_HPP
+

src/nest.cpp

@@ -18,8 +18,7 @@
  */
 
 #include "nest.hpp"
-#include "animation/ease.hpp"
-#include "math.hpp"
+#include "math/math.hpp"
 
 nest::nest()
 {
@@ -30,7 +29,7 @@ float3 nest::extend_shaft(shaft& shaft)
 {
 	float3 dig_position = get_shaft_position(shaft, shaft.current_depth);
 
-	float dr = frand(dig_radius * 0.75f, dig_radius * 1.25f);
+	float dr = math::random(dig_radius * 0.75f, dig_radius * 1.25f);
 
 	shaft.current_depth += dr * 0.1f;
 
@@ -39,15 +38,15 @@ float3 nest::extend_shaft(shaft& shaft)
 
 float3 nest::expand_chamber(chamber& chamber)
 {
-	float dig_angle = frand(0.0f, vmq::two_pi<float>);
-	float2 dig_direction = vmq::normalize(float2{std::cos(dig_angle), std::sin(dig_angle)});
+	float dig_angle = math::random(0.0f, math::two_pi<float>);
+	float2 dig_direction = math::normalize(float2{std::cos(dig_angle), std::sin(dig_angle)});
 
 	float3 chamber_center = get_shaft_position(*chamber.shaft, chamber.depth);
 	float3 dig_position = chamber_center;
 
-	float dr = frand(dig_radius * 0.75f, dig_radius * 1.25f);
+	float dr = math::random(dig_radius * 0.75f, dig_radius * 1.25f);
 
-	float t = frand(0.0f, 1.0f);
+	float t = math::random(0.0f, 1.0f);
 	dig_position.x += dig_direction.x * (chamber.outer_radius - dr) * t;
 	dig_position.z += dig_direction.y * (chamber.outer_radius - dr) * t;
 
@@ -63,8 +62,8 @@ float nest::get_shaft_angle(const shaft& shaft, float depth) const
 {
 	float shaft_length = shaft.depth[1] - shaft.depth[0];
 	float depth_factor = (depth - shaft.depth[0]) / shaft_length;
-	float pitch = ease<float>::linear(shaft.pitch[0], shaft.pitch[1], depth_factor);
-	return shaft.rotation + (depth / pitch) * shaft.chirality * vmq::two_pi<float>;
+	float pitch = math::lerp<float>(shaft.pitch[0], shaft.pitch[1], depth_factor);
+	return shaft.rotation + (depth / pitch) * shaft.chirality * math::two_pi<float>;
 }
 
 float nest::get_shaft_depth(const shaft& shaft, float turns) const
@@ -77,11 +76,11 @@ float3 nest::get_shaft_position(const shaft& shaft, float depth) const
 	float shaft_length = shaft.depth[1] - shaft.depth[0];
 	float depth_factor = (depth - shaft.depth[0]) / shaft_length;
 
-	float pitch = ease<float>::linear(shaft.pitch[0], shaft.pitch[1], depth_factor);
-	float radius = ease<float>::out_expo(shaft.radius[0], shaft.radius[1], depth_factor);
-	float translation_x = ease<float>::linear(shaft.translation[0][0], shaft.translation[1][0], depth_factor);
-	float translation_z = ease<float>::linear(shaft.translation[0][1], shaft.translation[1][1], depth_factor);
-	float angle = shaft.rotation + (depth / pitch) * shaft.chirality  * vmq::two_pi<float>;
+	float pitch = math::lerp<float>(shaft.pitch[0], shaft.pitch[1], depth_factor);
+	float radius = math::lerp<float>(shaft.radius[0], shaft.radius[1], depth_factor);
+	float translation_x = math::lerp<float>(shaft.translation[0][0], shaft.translation[1][0], depth_factor);
+	float translation_z = math::lerp<float>(shaft.translation[0][1], shaft.translation[1][1], depth_factor);
+	float angle = shaft.rotation + (depth / pitch) * shaft.chirality  * math::two_pi<float>;
 
 	float3 position;
 	position[0] = std::cos(angle) * radius + translation_x;

src/nest.hpp

@@ -21,13 +21,10 @@
 #define ANTKEEPER_NEST_HPP
 
 #include "geometry/mesh.hpp"
-#include <vmq/vmq.hpp>
+#include "utility/fundamental-types.hpp"
 #include <array>
 #include <vector>
 
-using namespace vmq::types;
-using namespace vmq::operators;
-
 class nest
 {
 public:

src/orbit-cam.cpp

@@ -19,23 +19,16 @@
 
 #include "orbit-cam.hpp"
 #include "scene/camera.hpp"
+#include "math/math.hpp"
 #include <algorithm>
 #include <cmath>
 
-using namespace vmq::operators;
-
-template <class T>
-static inline T lerp(const T& x, const T& y, float a)
-{
-	return x * (1.0f - a) + y * a;
-}
-
 orbit_cam::orbit_cam():
-	elevation_rotation(vmq::identity_quaternion<float>),
-	azimuth_rotation(vmq::identity_quaternion<float>),
-	target_elevation_rotation(vmq::identity_quaternion<float>),
-	target_azimuth_rotation(vmq::identity_quaternion<float>),
-	target_rotation(vmq::identity_quaternion<float>)
+	elevation_rotation(math::identity_quaternion<float>),
+	azimuth_rotation(math::identity_quaternion<float>),
+	target_elevation_rotation(math::identity_quaternion<float>),
+	target_azimuth_rotation(math::identity_quaternion<float>),
+	target_rotation(math::identity_quaternion<float>)
 {}
 
 orbit_cam::~orbit_cam()
@@ -47,7 +40,7 @@ void orbit_cam::update(float dt)
 
 	// Calculate rotation and target rotation quaternions
 	//rotation = azimuth_rotation * elevation_rotation;
-	target_rotation = vmq::normalize(target_azimuth_rotation * target_elevation_rotation);
+	target_rotation = math::normalize(target_azimuth_rotation * target_elevation_rotation);
 	
 	// Calculate target translation
 	target_translation = target_focal_point + target_rotation * float3{0.0f, 0.0f, target_focal_distance};
@@ -56,15 +49,15 @@ void orbit_cam::update(float dt)
 	//rotation = glm::mix(rotation, target_rotation, interpolation_factor);
 	
 	// Interpolate angles
-	set_elevation(lerp(elevation, target_elevation, interpolation_factor));
-	set_azimuth(lerp(azimuth, target_azimuth, interpolation_factor));
+	set_elevation(math::lerp(elevation, target_elevation, interpolation_factor));
+	set_azimuth(math::lerp(azimuth, target_azimuth, interpolation_factor));
 	
 	// Calculate rotation
-	set_rotation(vmq::normalize(azimuth_rotation * elevation_rotation));
+	set_rotation(math::normalize(azimuth_rotation * elevation_rotation));
 	
 	// Interpolate focal point and focal distance
-	focal_point = vmq::lerp(focal_point, target_focal_point, interpolation_factor);
-	focal_distance = lerp(focal_distance, target_focal_distance, interpolation_factor);
+	focal_point = math::lerp(focal_point, target_focal_point, interpolation_factor);
+	focal_distance = math::lerp(focal_distance, target_focal_distance, interpolation_factor);
 	
 	// Caluclate translation
 	set_translation(focal_point + get_rotation() * float3{0.0f, 0.0f, focal_distance});
@@ -87,7 +80,7 @@ void orbit_cam::update(float dt)
 	// Update camera
 	if (get_camera() != nullptr)
 	{
-		transform<float> transform = vmq::identity_transform<float>;
+		transform_type transform = math::identity_transform<float>;
 		transform.translation = get_translation();
 		transform.rotation = get_rotation();
 		get_camera()->set_transform(transform);
@@ -128,13 +121,13 @@ void orbit_cam::set_focal_distance(float distance)
 void orbit_cam::set_elevation(float angle)
 {
 	elevation = angle;
-	elevation_rotation = vmq::angle_axis(elevation, float3{-1.0f, 0.0f, 0.0f});
+	elevation_rotation = math::angle_axis(elevation, float3{-1.0f, 0.0f, 0.0f});
 }
 
 void orbit_cam::set_azimuth(float angle)
 {
 	azimuth = angle;
-	azimuth_rotation = vmq::angle_axis(azimuth, float3{0.0f, 1.0f, 0.0f});
+	azimuth_rotation = math::angle_axis(azimuth, float3{0.0f, 1.0f, 0.0f});
 }
 
 void orbit_cam::set_target_focal_point(const float3& point)
@@ -150,12 +143,12 @@ void orbit_cam::set_target_focal_distance(float distance)
 void orbit_cam::set_target_elevation(float angle)
 {
 	target_elevation = angle;
-	target_elevation_rotation = vmq::angle_axis(target_elevation, float3{-1.0f, 0.0f, 0.0f});
+	target_elevation_rotation = math::angle_axis(target_elevation, float3{-1.0f, 0.0f, 0.0f});
 }
 
 void orbit_cam::set_target_azimuth(float angle)
 {
 	target_azimuth = angle;
-	target_azimuth_rotation = vmq::angle_axis(target_azimuth, float3{0.0f, 1.0f, 0.0f});
+	target_azimuth_rotation = math::angle_axis(target_azimuth, float3{0.0f, 1.0f, 0.0f});
 }
 

src/orbit-cam.hpp

@@ -56,7 +56,7 @@ public:
 	float get_target_elevation() const;
 	float get_target_azimuth() const;
 	const float3& get_target_translation() const;
-	const vmq::quaternion<float>& get_target_rotation() const;
+	const quaternion_type& get_target_rotation() const;
 	
 private:
 	float3 focal_point;
@@ -69,11 +69,11 @@ private:
 	float target_elevation;
 	float target_azimuth;
 	
-	vmq::quaternion<float> elevation_rotation;
-	vmq::quaternion<float> azimuth_rotation;
-	vmq::quaternion<float> target_elevation_rotation;
-	vmq::quaternion<float> target_azimuth_rotation;
-	vmq::quaternion<float> target_rotation;
+	quaternion_type elevation_rotation;
+	quaternion_type azimuth_rotation;
+	quaternion_type target_elevation_rotation;
+	quaternion_type target_azimuth_rotation;
+	quaternion_type target_rotation;
 	float3 target_translation;
 };
 
@@ -122,7 +122,7 @@ inline const float3& orbit_cam::get_target_translation() const
 	return target_translation;
 }
 
-inline const vmq::quaternion<float>& orbit_cam::get_target_rotation() const
+inline const typename camera_rig::quaternion_type& orbit_cam::get_target_rotation() const
 {
 	return target_rotation;
 }

src/pheromone-matrix.cpp

@@ -18,7 +18,7 @@
  */
 
 #include "pheromone-matrix.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 
 void convolve(pheromone_matrix* matrix, const float* kernel, int kernel_size)
 {
@@ -75,9 +75,9 @@ void evaporate(pheromone_matrix* matrix, float factor)
 
 void diffuse(pheromone_matrix* matrix)
 {
-	const vmq::matrix<float, 3, 3> diffusion_kernel =
-		vmq::mul(
-			vmq::matrix<float, 3, 3>
+	const math::matrix<float, 3, 3> diffusion_kernel =
+		math::mul(
+			math::matrix<float, 3, 3>
 			{{
 				{1, 2, 1},
 				{2, 4, 2},

src/rasterizer/shader-input.hpp

@@ -20,14 +20,13 @@
 #ifndef ANTKEEPER_SHADER_INPUT_HPP
 #define ANTKEEPER_SHADER_INPUT_HPP
 
-#include <vmq/vmq.hpp>
+#include "utility/fundamental-types.hpp"
 #include <string>
 
 class shader_program;
 class texture_2d;
 class texture_cube;
 enum class shader_variable_type;
-using namespace vmq::types;
 
 /**
  * Port through which data can be uploaded to shader variables.

src/renderer/material-property.hpp

@@ -23,13 +23,10 @@
 #include "animation/tween.hpp"
 #include "rasterizer/shader-variable-type.hpp"
 #include "rasterizer/shader-input.hpp"
-#include "animation/ease.hpp"
-#include <vmq/vmq.hpp>
+#include "math/interpolation.hpp"
+#include "utility/fundamental-types.hpp"
 #include <cstdlib>
 
-using namespace vmq::types;
-using namespace vmq::operators;
-
 class material;
 class shader_program;
 class texture_2d;
@@ -191,25 +188,25 @@ inline T material_property::default_interpolator(const T& x, const T& y, doub
 template <>
 inline float material_property<float>::default_interpolator(const float& x, const float& y, double a)
 {
-	return ease<float, float>::linear(x, y, static_cast<float>(a));
+	return math::lerp<float, float>(x, y, static_cast<float>(a));
 }
 
 template <>
 inline float2 material_property<float2>::default_interpolator(const float2& x, const float2& y, double a)
 {
-	return ease<float2, float>::linear(x, y, static_cast<float>(a));
+	return math::lerp<float2, float>(x, y, static_cast<float>(a));
 }
 
 template <>
 inline float3 material_property<float3>::default_interpolator(const float3& x, const float3& y, double a)
 {
-	return ease<float3, float>::linear(x, y, static_cast<float>(a));
+	return math::lerp<float3, float>(x, y, static_cast<float>(a));
 }
 
 template <>
 inline float4 material_property<float4>::default_interpolator(const float4& x, const float4& y, double a)
 {
-	return ease<float4, float>::linear(x, y, static_cast<float>(a));
+	return math::lerp<float4, float>(x, y, static_cast<float>(a));
 }
 
 template <class T>

src/renderer/passes/bloom-pass.cpp

@@ -32,12 +32,10 @@
 #include "rasterizer/texture-filter.hpp"
 #include "renderer/vertex-attributes.hpp"
 #include "renderer/render-context.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 #include <cmath>
 #include <glad/glad.h>
 
-using namespace vmq;
-
 bloom_pass::bloom_pass(::rasterizer* rasterizer, const ::framebuffer* framebuffer, resource_manager* resource_manager):
 	render_pass(rasterizer, framebuffer),
 	source_texture(nullptr),

src/renderer/passes/bloom-pass.hpp

@@ -21,9 +21,7 @@
 #define ANTKEEPER_BLOOM_PASS_HPP
 
 #include "renderer/render-pass.hpp"
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "math/math.hpp"
 
 class shader_program;
 class shader_input;

src/renderer/passes/clear-pass.hpp

@@ -21,9 +21,7 @@
 #define ANTKEEPER_CLEAR_PASS_HPP
 
 #include "renderer/render-pass.hpp"
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "utility/fundamental-types.hpp"
 
 /**
  * Clears the color, depth, or stencil buffer of a render target.

src/renderer/passes/final-pass.cpp

@@ -32,12 +32,10 @@
 #include "rasterizer/texture-filter.hpp"
 #include "renderer/vertex-attributes.hpp"
 #include "renderer/render-context.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 #include <cmath>
 #include <glad/glad.h>
 
-using namespace vmq;
-
 final_pass::final_pass(::rasterizer* rasterizer, const ::framebuffer* framebuffer, resource_manager* resource_manager):
 	render_pass(rasterizer, framebuffer),
 	color_texture(nullptr),

src/renderer/passes/final-pass.hpp

@@ -21,9 +21,7 @@
 #define ANTKEEPER_FINAL_PASS_HPP
 
 #include "renderer/render-pass.hpp"
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "math/math.hpp"
 
 class shader_program;
 class shader_input;

src/renderer/passes/material-pass.cpp

@@ -45,15 +45,13 @@
 #include "scene/spotlight.hpp"
 #include "scene/scene.hpp"
 #include "configuration.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 #include <cmath>
 #include <glad/glad.h>
 #include <iostream>
 
 #include "shadow-map-pass.hpp"
 
-using namespace vmq;
-
 static bool operation_compare(const render_operation& a, const render_operation& b);
 
 material_pass::material_pass(::rasterizer* rasterizer, const ::framebuffer* framebuffer, resource_manager* resource_manager):
@@ -174,7 +172,7 @@ void material_pass::render(render_context* context) const
 					
 					// Transform position into view-space
 					float3 position = light->get_transform_tween().interpolate(context->alpha).translation;
-					float3 view_space_position = vmq::resize<3>(view * float4{position.x, position.y, position.z, 1.0f});
+					float3 view_space_position = math::resize<3>(view * float4{position.x, position.y, position.z, 1.0f});
 					point_light_positions[point_light_count] = view_space_position;
 					
 					point_light_attenuations[point_light_count] = static_cast<const point_light*>(light)->get_attenuation_tween().interpolate(context->alpha);
@@ -192,7 +190,7 @@ void material_pass::render(render_context* context) const
 					
 					// Transform direction into view-space
 					float3 direction = static_cast<const directional_light*>(light)->get_direction_tween().interpolate(context->alpha);
-					float3 view_space_direction = vmq::normalize(vmq::resize<3>(view * vmq::resize<4>(-direction)));
+					float3 view_space_direction = math::normalize(math::resize<3>(view * math::resize<4>(-direction)));
 					directional_light_directions[directional_light_count] = view_space_direction;
 					
 					++directional_light_count;
@@ -209,14 +207,14 @@ void material_pass::render(render_context* context) const
 					
 					// Transform position into view-space
 					float3 position = light->get_transform_tween().interpolate(context->alpha).translation;
-					float3 view_space_position = vmq::resize<3>(view * float4{position.x, position.y, position.z, 1.0f});
+					float3 view_space_position = math::resize<3>(view * float4{position.x, position.y, position.z, 1.0f});
 					spotlight_positions[spotlight_count] = view_space_position;
 					
 					const ::spotlight* spotlight = static_cast<const ::spotlight*>(light);
 					
 					// Transform direction into view-space
 					float3 direction = spotlight->get_direction_tween().interpolate(context->alpha);
-					float3 view_space_direction = vmq::normalize(vmq::resize<3>(view * vmq::resize<4>(-direction)));
+					float3 view_space_direction = math::normalize(math::resize<3>(view * math::resize<4>(-direction)));
 					spotlight_directions[spotlight_count] = view_space_direction;
 					
 					spotlight_attenuations[spotlight_count] = spotlight->get_attenuation_tween().interpolate(context->alpha);
@@ -414,7 +412,7 @@ void material_pass::render(render_context* context) const
 		model = operation.transform;
 		model_view_projection = view_projection * model;
 		model_view = view * model;
-		normal_model_view = vmq::transpose(vmq::inverse(vmq::resize<3, 3>(model_view)));
+		normal_model_view = math::transpose(math::inverse(math::resize<3, 3>(model_view)));
 
 		// Upload operation-dependent parameters
 		if (parameters->model)

src/renderer/passes/material-pass.hpp

@@ -23,6 +23,7 @@
 #include "renderer/render-pass.hpp"
 #include "renderer/material.hpp"
 #include "animation/tween.hpp"
+#include "utility/fundamental-types.hpp"
 #include <unordered_map>
 
 class camera;

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

@@ -34,12 +34,9 @@
 #include "geometry/view-frustum.hpp"
 #include "geometry/aabb.hpp"
 #include "configuration.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 #include <cmath>
 #include <glad/glad.h>
-#include <iostream>
-
-using namespace vmq;
 
 static bool operation_compare(const render_operation& a, const render_operation& b);
 
@@ -74,13 +71,13 @@ shadow_map_pass::shadow_map_pass(::rasterizer* rasterizer, const ::framebuffer*
 	skinned_model_view_projection_input = skinned_shader_program->get_input("model_view_projection");
 	
 	// Calculate bias-tile matrices
-	float4x4 bias_matrix = vmq::translate(vmq::identity4x4<float>, float3{0.5f, 0.5f, 0.5f}) * vmq::scale(vmq::identity4x4<float>, float3{0.5f, 0.5f, 0.5f});
-	float4x4 tile_scale = vmq::scale(vmq::identity4x4<float>, float3{0.5f, 0.5f, 1.0f});
+	float4x4 bias_matrix = math::translate(math::identity4x4<float>, float3{0.5f, 0.5f, 0.5f}) * math::scale(math::identity4x4<float>, float3{0.5f, 0.5f, 0.5f});
+	float4x4 tile_scale = math::scale(math::identity4x4<float>, float3{0.5f, 0.5f, 1.0f});
 	for (int i = 0; i < 4; ++i)
 	{
 		float x = static_cast<float>(i % 2) * 0.5f;
 		float y = static_cast<float>(i / 2) * 0.5f;
-		float4x4 tile_matrix = vmq::translate(vmq::identity4x4<float>, float3{x, y, 0.0f}) * tile_scale;
+		float4x4 tile_matrix = math::translate(math::identity4x4<float>, float3{x, y, 0.0f}) * tile_scale;
 		bias_tile_matrices[i] = tile_matrix * bias_matrix;
 	}
 }
@@ -151,11 +148,11 @@ void shadow_map_pass::render(render_context* context) const
 	}
 	
 	// Calculate a view-projection matrix from the directional light's transform
-	transform<float> light_transform = light->get_transform_tween().interpolate(context->alpha);
+	math::transform<float> light_transform = light->get_transform_tween().interpolate(context->alpha);
 	float3 forward = light_transform.rotation * global_forward;
 	float3 up = light_transform.rotation * global_up;
-	float4x4 light_view = vmq::look_at(light_transform.translation, light_transform.translation + forward, up);
-	float4x4 light_projection = vmq::ortho(-1.0f, 1.0f, -1.0f, 1.0f, -1.0f, 1.0f);
+	float4x4 light_view = math::look_at(light_transform.translation, light_transform.translation + forward, up);
+	float4x4 light_projection = math::ortho(-1.0f, 1.0f, -1.0f, 1.0f, -1.0f, 1.0f);
 	float4x4 light_view_projection = light_projection * light_view;
 	
 	// Get the camera's view matrix
@@ -179,7 +176,7 @@ void shadow_map_pass::render(render_context* context) const
 		// Calculate projection matrix for view camera subfrustum
 		const float subfrustum_near = split_distances[i];
 		const float subfrustum_far = split_distances[i + 1];
-		float4x4 subfrustum_projection = vmq::perspective(camera.get_fov(), camera.get_aspect_ratio(), subfrustum_near, subfrustum_far);
+		float4x4 subfrustum_projection = math::perspective(camera.get_fov(), camera.get_aspect_ratio(), subfrustum_near, subfrustum_far);
 		
 		// Calculate view camera subfrustum
 		view_frustum<float> subfrustum(subfrustum_projection * camera_view);
@@ -222,7 +219,7 @@ void shadow_map_pass::render(render_context* context) const
 		offset.y = std::ceil(offset.y * half_shadow_map_resolution) / half_shadow_map_resolution;
 		
 		// Crop the light view-projection matrix
-		crop_matrix = vmq::translate(vmq::identity4x4<float>, offset) * vmq::scale(vmq::identity4x4<float>, scale);
+		crop_matrix = math::translate(math::identity4x4<float>, offset) * math::scale(math::identity4x4<float>, scale);
 		cropped_view_projection = crop_matrix * light_view_projection;
 		
 		// Calculate shadow matrix

src/renderer/passes/shadow-map-pass.hpp

@@ -21,10 +21,7 @@
 #define ANTKEEPER_SHADOW_MAP_PASS_HPP
 
 #include "renderer/render-pass.hpp"
-
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "utility/fundamental-types.hpp"
 
 class shader_program;
 class shader_input;
@@ -86,4 +83,3 @@ inline const float* shadow_map_pass::get_split_distances() const
 }
 
 #endif // ANTKEEPER_SHADOW_MAP_PASS_HPP
-

src/renderer/passes/sky-pass.cpp

@@ -37,12 +37,10 @@
 #include "scene/ambient-light.hpp"
 #include "scene/directional-light.hpp"
 #include "scene/scene.hpp"
-#include <vmq/vmq.hpp>
+#include "utility/fundamental-types.hpp"
 #include <cmath>
 #include <glad/glad.h>
 
-using namespace vmq;
-
 sky_pass::sky_pass(::rasterizer* rasterizer, const ::framebuffer* framebuffer, resource_manager* resource_manager):
 	render_pass(rasterizer, framebuffer)
 {
@@ -113,7 +111,7 @@ void sky_pass::render(render_context* context) const
 	rasterizer->set_viewport(0, 0, std::get<0>(viewport), std::get<1>(viewport));
 
 	float3 sun_direction = {0, 0, -1};
-	float sun_angular_radius = 3.0f * vmq::pi<float> / 180.0f;
+	float sun_angular_radius = math::radians<float>(3.0f);
 	
 	// Find sun direction
 	const std::list<scene_object_base*>* lights = context->scene->get_objects(light::object_type_id);
@@ -129,9 +127,9 @@ void sky_pass::render(render_context* context) const
 	}
 
 	// Calculate matrix
-	float4x4 model_view = vmq::resize<4, 4>(vmq::resize<3, 3>(context->camera->get_view_tween().interpolate(context->alpha)));
-	float4x4 inverse_projection = vmq::inverse(context->camera->get_projection_tween().interpolate(context->alpha));
-	float4x4 matrix = vmq::inverse(model_view) * inverse_projection;
+	float4x4 model_view = math::resize<4, 4>(math::resize<3, 3>(context->camera->get_view_tween().interpolate(context->alpha)));
+	float4x4 inverse_projection = math::inverse(context->camera->get_projection_tween().interpolate(context->alpha));
+	float4x4 matrix = math::inverse(model_view) * inverse_projection;
 
 	// Change shader program
 	rasterizer->use_program(*shader_program);

src/renderer/passes/ui-pass.cpp

@@ -41,12 +41,10 @@
 #include "scene/directional-light.hpp"
 #include "scene/scene.hpp"
 #include "scene/billboard.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 #include <cmath>
 #include <glad/glad.h>
 
-using namespace vmq;
-
 ui_pass::ui_pass(::rasterizer* rasterizer, const ::framebuffer* framebuffer, resource_manager* resource_manager):
 	render_pass(rasterizer, framebuffer),
 	time(0.0f)

src/renderer/render-context.hpp

@@ -23,19 +23,16 @@
 #include "renderer/render-operation.hpp"
 #include "geometry/plane.hpp"
 #include "geometry/bounding-volume.hpp"
-#include <vmq/vmq.hpp>
+#include "utility/fundamental-types.hpp"
 #include <list>
 
 class camera;
 class scene;
 
-using namespace vmq::types;
-using vmq::transform;
-
 struct render_context
 {
 	const camera* camera;
-	transform<float> camera_transform;
+	math::transform<float> camera_transform;
 	float3 camera_forward;
 	float3 camera_up;
 	const bounding_volume<float>* camera_culling_volume;

src/renderer/render-operation.hpp

@@ -20,14 +20,13 @@
 #ifndef ANTKEEPER_RENDER_OPERATION_HPP
 #define ANTKEEPER_RENDER_OPERATION_HPP
 
+#include "utility/fundamental-types.hpp"
 #include <cstdlib>
-#include <vmq/vmq.hpp>
 
 class pose;
 class material;
 class vertex_array;
 enum class drawing_mode;
-using namespace vmq::types;
 
 /**
  * Describes a render operation with a single mesh and single material.

src/renderer/renderer.cpp

@@ -27,13 +27,12 @@
 #include "scene/lod-group.hpp"
 #include "renderer/model.hpp"
 #include "rasterizer/drawing-mode.hpp"
+#include "math/math.hpp"
+#include "geometry/projection.hpp"
 #include "configuration.hpp"
-#include "math.hpp"
 #include <functional>
 #include <set>
 
-using namespace vmq::operators;
-
 renderer::renderer()
 {
 	// Setup billboard render operation
@@ -166,8 +165,8 @@ void renderer::process_model_instance(render_context& context, const ::model_ins
 		operation.drawing_mode = group->get_drawing_mode();
 		operation.start_index = group->get_start_index();
 		operation.index_count = group->get_index_count();
-		operation.transform = vmq::matrix_cast(model_instance->get_transform_tween().interpolate(context.alpha));
-		operation.depth = signed_distance(context.clip_near, vmq::resize<3>(operation.transform[3]));
+		operation.transform = math::matrix_cast(model_instance->get_transform_tween().interpolate(context.alpha));
+		operation.depth = context.clip_near.signed_distance(math::resize<3>(operation.transform[3]));
 		operation.instance_count = model_instance->get_instance_count();
 
 		context.operations.push_back(operation);
@@ -185,26 +184,26 @@ void renderer::process_billboard(render_context& context, const ::billboard* bil
 	if (!context.camera_culling_volume->intersects(*object_culling_volume))
 		return;
 	
-	transform<float> billboard_transform = billboard->get_transform_tween().interpolate(context.alpha);
+	math::transform<float> billboard_transform = billboard->get_transform_tween().interpolate(context.alpha);
 	billboard_op.material = billboard->get_material();
-	billboard_op.depth = signed_distance(context.clip_near, vmq::resize<3>(billboard_transform.translation));
+	billboard_op.depth = context.clip_near.signed_distance(math::resize<3>(billboard_transform.translation));
 	
 	// Align billboard
 	if (billboard->get_billboard_type() == billboard_type::spherical)
 	{
-		billboard_transform.rotation = vmq::normalize(vmq::look_rotation(context.camera_forward, context.camera_up) * billboard_transform.rotation);
+		billboard_transform.rotation = math::normalize(math::look_rotation(context.camera_forward, context.camera_up) * billboard_transform.rotation);
 	}
 	else if (billboard->get_billboard_type() == billboard_type::cylindrical)
 	{
 		const float3& alignment_axis = billboard->get_alignment_axis();
-		float3 look = vmq::normalize(math::project_on_plane(billboard_transform.translation - context.camera_transform.translation, {0.0f, 0.0f, 0.0f}, alignment_axis));
-		float3 right = vmq::normalize(vmq::cross(alignment_axis, look));
-		look = vmq::cross(right, alignment_axis);
-		float3 up = vmq::cross(look, right);
-		billboard_transform.rotation = vmq::normalize(vmq::look_rotation(look, up) * billboard_transform.rotation);
+		float3 look = math::normalize(project_on_plane(billboard_transform.translation - context.camera_transform.translation, {0.0f, 0.0f, 0.0f}, alignment_axis));
+		float3 right = math::normalize(math::cross(alignment_axis, look));
+		look = math::cross(right, alignment_axis);
+		float3 up = math::cross(look, right);
+		billboard_transform.rotation = math::normalize(math::look_rotation(look, up) * billboard_transform.rotation);
 	}
 	
-	billboard_op.transform = vmq::matrix_cast(billboard_transform);
+	billboard_op.transform = math::matrix_cast(billboard_transform);
 	
 	context.operations.push_back(billboard_op);
 }

src/renderer/simple-render-pass.cpp

@@ -33,11 +33,9 @@
 #include "renderer/render-context.hpp"
 #include "renderer/material.hpp"
 #include "renderer/material-property.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 #include <glad/glad.h>
 
-using namespace vmq;
-
 simple_render_pass::simple_render_pass(::rasterizer* rasterizer, const ::framebuffer* framebuffer, ::shader_program* shader_program):
 	render_pass(rasterizer, framebuffer),
 	shader_program(shader_program),

src/renderer/simple-render-pass.hpp

@@ -22,9 +22,7 @@
 
 #include "renderer/render-pass.hpp"
 #include "animation/tween.hpp"
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "utility/fundamental-types.hpp"
 
 class shader_program;
 class shader_input;

src/resources/material-loader.cpp

@@ -24,6 +24,7 @@
 #include "rasterizer/texture-wrapping.hpp"
 #include "rasterizer/texture-filter.hpp"
 #include "rasterizer/texture-2d.hpp"
+#include "utility/fundamental-types.hpp"
 #include "string-table.hpp"
 #include <cctype>
 #include <vector>

src/resources/mesh-loader.cpp

@@ -20,6 +20,7 @@
 #include "resources/resource-loader.hpp"
 #include "geometry/mesh.hpp"
 #include "geometry/mesh-functions.hpp"
+#include "utility/fundamental-types.hpp"
 #include <sstream>
 #include <stdexcept>
 #include <physfs.h>

src/resources/model-loader.cpp

@@ -23,13 +23,11 @@
 #include "renderer/vertex-attributes.hpp"
 #include "rasterizer/vertex-attribute-type.hpp"
 #include "rasterizer/drawing-mode.hpp"
+#include "utility/fundamental-types.hpp"
 #include <sstream>
 #include <stdexcept>
 #include <limits>
 #include <physfs.h>
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
 
 struct material_group
 {

src/scene/billboard.hpp

@@ -22,6 +22,7 @@
 
 #include "scene/scene-object.hpp"
 #include "geometry/aabb.hpp"
+#include "utility/fundamental-types.hpp"
 
 class material;
 

src/scene/camera.cpp

@@ -20,22 +20,22 @@
 #include "scene/camera.hpp"
 #include "configuration.hpp"
 #include "animation/ease.hpp"
-
-using namespace vmq::operators;
+#include "math/constants.hpp"
+#include "math/interpolation.hpp"
 
 static float4x4 interpolate_view(const camera* camera, const float4x4& x, const float4x4& y, float a)
 {
-	transform<float> transform = camera->get_transform_tween().interpolate(a);
+	math::transform<float> transform = camera->get_transform_tween().interpolate(a);
 	float3 forward = transform.rotation * global_forward;
 	float3 up = transform.rotation * global_up;
-	return vmq::look_at(transform.translation, transform.translation + forward, up);
+	return math::look_at(transform.translation, transform.translation + forward, up);
 }
 
 static float4x4 interpolate_projection(const camera* camera, const float4x4& x, const float4x4& y, float a)
 {
 	if (camera->is_orthographic())
 	{
-		return vmq::ortho(
+		return math::ortho(
 			camera->get_clip_left_tween().interpolate(a),
 			camera->get_clip_right_tween().interpolate(a),
 			camera->get_clip_bottom_tween().interpolate(a),
@@ -45,7 +45,7 @@ static float4x4 interpolate_projection(const camera* camera, const float4x4& x,
 	}
 	else
 	{
-		return vmq::perspective(
+		return math::perspective(
 			camera->get_fov_tween().interpolate(a),
 			camera->get_aspect_ratio_tween().interpolate(a),
 			camera->get_clip_near_tween().interpolate(a),
@@ -62,17 +62,17 @@ camera::camera():
 	compositor(nullptr),
 	composite_index(0),
 	orthographic(true),
-	clip_left(-1.0f, ease<float>::linear),
-	clip_right(1.0f, ease<float>::linear),
-	clip_bottom(-1.0f, ease<float>::linear),
-	clip_top(1.0f, ease<float>::linear),
-	clip_near(-1.0f, ease<float>::linear),
-	clip_far(1.0f, ease<float>::linear),
-	fov(vmq::half_pi<float>, ease<float>::linear),
-	aspect_ratio(1.0f, ease<float>::linear),
-	view(vmq::identity4x4<float>, std::bind(&interpolate_view, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)),
-	projection(vmq::identity4x4<float>, std::bind(&interpolate_projection, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)),
-	view_projection(vmq::identity4x4<float>, std::bind(&interpolate_view_projection, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3))
+	clip_left(-1.0f, math::lerp<float, float>),
+	clip_right(1.0f, math::lerp<float, float>),
+	clip_bottom(-1.0f, math::lerp<float, float>),
+	clip_top(1.0f, math::lerp<float, float>),
+	clip_near(-1.0f, math::lerp<float, float>),
+	clip_far(1.0f, math::lerp<float, float>),
+	fov(math::half_pi<float>, math::lerp<float, float>),
+	aspect_ratio(1.0f, math::lerp<float, float>),
+	view(math::identity4x4<float>, std::bind(&interpolate_view, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)),
+	projection(math::identity4x4<float>, std::bind(&interpolate_projection, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3)),
+	view_projection(math::identity4x4<float>, std::bind(&interpolate_view_projection, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3))
 {}
 
 float3 camera::project(const float3& object, const float4& viewport) const
@@ -85,7 +85,7 @@ float3 camera::project(const float3& object, const float4& viewport) const
 	result[0] = result[0] * viewport[2] + viewport[0];
 	result[1] = result[1] * viewport[3] + viewport[1];
 	
-	return vmq::resize<3>(result);
+	return math::resize<3>(result);
 }
 
 float3 camera::unproject(const float3& window, const float4& viewport) const
@@ -96,9 +96,9 @@ float3 camera::unproject(const float3& window, const float4& viewport) const
 	result[2] = window[2] * 2.0f - 1.0f;
 	result[3] = 1.0f;
 	
-	result = vmq::inverse(view_projection[1]) * result;
+	result = math::inverse(view_projection[1]) * result;
 
-	return vmq::resize<3>(result) * (1.0f / result[3]);
+	return math::resize<3>(result) * (1.0f / result[3]);
 }
 
 void camera::set_perspective(float fov, float aspect_ratio, float clip_near, float clip_far)
@@ -110,7 +110,7 @@ void camera::set_perspective(float fov, float aspect_ratio, float clip_near, flo
 	this->clip_near[1] = clip_near;
 	this->clip_far[1] = clip_far;
 
-	projection[1] = vmq::perspective(fov, aspect_ratio, clip_near, clip_far);
+	projection[1] = math::perspective(fov, aspect_ratio, clip_near, clip_far);
 	
 	// Recalculate view-projection matrix
 	view_projection[1] = projection[1] * view[1];
@@ -130,7 +130,7 @@ void camera::set_orthographic(float clip_left, float clip_right, float clip_bott
 	this->clip_near[1] = clip_near;
 	this->clip_far[1] = clip_far;
 
-	projection[1] = vmq::ortho(clip_left, clip_right, clip_bottom, clip_top, clip_near, clip_far);
+	projection[1] = math::ortho(clip_left, clip_right, clip_bottom, clip_top, clip_near, clip_far);
 	
 	// Recalculate view-projection matrix
 	view_projection[1] = projection[1] * view[1];
@@ -170,7 +170,7 @@ void camera::transformed()
 	// Recalculate view and view-projection matrices
 	float3 forward = get_rotation() * global_forward;
 	float3 up = get_rotation() * global_up;
-	view[1] = vmq::look_at(get_translation(), get_translation() + forward, up);
+	view[1] = math::look_at(get_translation(), get_translation() + forward, up);
 	view_projection[1] = projection[1] * view[1];
 	
 	// Recalculate view frustum

src/scene/camera.hpp

@@ -22,6 +22,7 @@
 
 #include "scene/scene-object.hpp"
 #include "geometry/view-frustum.hpp"
+#include "utility/fundamental-types.hpp"
 
 class compositor;
 

src/scene/directional-light.cpp

@@ -19,14 +19,13 @@
 
 #include "directional-light.hpp"
 #include "configuration.hpp"
-
-using namespace vmq::operators;
+#include "math/math.hpp"
 
 static float3 interpolate_direction(const float3& x, const float3& y, float a)
 {
-	quaternion<float> q0 = vmq::rotation(global_forward, x);
-	quaternion<float> q1 = vmq::rotation(global_forward, y);
-	return vmq::normalize(vmq::slerp(q0, q1, a) * global_forward);
+	math::quaternion<float> q0 = math::rotation(global_forward, x);
+	math::quaternion<float> q1 = math::rotation(global_forward, y);
+	return math::normalize(math::slerp(q0, q1, a) * global_forward);
 }
 
 directional_light::directional_light():
@@ -41,6 +40,6 @@ void directional_light::update_tweens()
 
 void directional_light::transformed()
 {
-	direction[1] = vmq::normalize(get_transform().rotation * global_forward);
+	direction[1] = math::normalize(get_transform().rotation * global_forward);
 }
 

src/scene/directional-light.hpp

@@ -21,9 +21,7 @@
 #define ANTKEEPER_DIRECTIONAL_LIGHT_HPP
 
 #include "scene/light.hpp"
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "utility/fundamental-types.hpp"
 
 class directional_light: public light
 {

src/scene/light.cpp

@@ -19,14 +19,13 @@
 
 #include "scene/light.hpp"
 #include "animation/ease.hpp"
-
-using namespace vmq::operators;
+#include "math/interpolation.hpp"
 
 light::light():
 	bounds(get_translation(), 0.0f),
-	color(float3{1.0f, 1.0f, 1.0f}, ease<float3>::linear),
-	intensity(1.0f, ease<float>::linear),
-	scaled_color(float3{1.0f, 1.0f, 1.0f}, ease<float3>::linear)
+	color(float3{1.0f, 1.0f, 1.0f}, math::lerp<float3, float>),
+	intensity(1.0f, math::lerp<float, float>),
+	scaled_color(float3{1.0f, 1.0f, 1.0f}, math::lerp<float3, float>)
 {}
 
 void light::set_color(const float3& color)

src/scene/light.hpp

@@ -22,6 +22,7 @@
 
 #include "scene/scene-object.hpp"
 #include "geometry/sphere.hpp"
+#include "utility/fundamental-types.hpp"
 
 enum class light_type
 {

src/scene/lod-group.cpp

@@ -37,7 +37,7 @@ void lod_group::resize(std::size_t level_count)
 
 std::size_t lod_group::select_lod(const ::camera& camera) const
 {
-	float distance = signed_distance(camera.get_view_frustum().get_near(), get_translation());
+	float distance = camera.get_view_frustum().get_near().signed_distance(get_translation());
 	
 	if (distance < 300.0f)
 		return 0;

src/scene/point-light.cpp

@@ -18,12 +18,10 @@
  */
 
 #include "point-light.hpp"
-#include "animation/ease.hpp"
-
-using namespace vmq::operators;
+#include "math/interpolation.hpp"
 
 point_light::point_light():
-	attenuation(float3{1, 0, 0}, ease<float3>::linear)
+	attenuation(float3{1, 0, 0}, math::lerp<float3, float>)
 {}
 
 void point_light::set_attenuation(const float3& attenuation)
@@ -36,4 +34,3 @@ void point_light::update_tweens()
 	light::update_tweens();
 	attenuation.update();
 }
-

src/scene/point-light.hpp

@@ -21,6 +21,7 @@
 #define ANTKEEPER_POINT_LIGHT_HPP
 
 #include "scene/light.hpp"
+#include "utility/fundamental-types.hpp"
 
 class point_light: public light
 {

src/scene/scene-object.cpp

@@ -18,24 +18,25 @@
  */
 
 #include "scene/scene-object.hpp"
+#include "math/math.hpp"
 
-static transform<float> transform_interpolate(const transform<float>& x, const transform<float>& y, float a)
+typename scene_object_base::transform_type scene_object_base::interpolate_transforms(const transform_type& x, const transform_type& y, float a)
 {
 	return
 		{
-			vmq::lerp(x.translation, y.translation, a),
-			vmq::slerp(x.rotation, y.rotation, a),
-			vmq::lerp(x.scale, y.scale, a),
+			math::lerp(x.translation, y.translation, a),
+			math::slerp(x.rotation, y.rotation, a),
+			math::lerp(x.scale, y.scale, a),
 		};
 }
 
 scene_object_base::scene_object_base():
 	active(true),
-	transform(vmq::identity_transform<float>, transform_interpolate),
+	transform(math::identity_transform<float>, interpolate_transforms),
 	culling_mask(nullptr)
 {}
 
-void scene_object_base::set_culling_mask(const bounding_volume<float>* culling_mask)
+void scene_object_base::set_culling_mask(const bounding_volume_type* culling_mask)
 {
 	this->culling_mask = culling_mask;
 }
@@ -51,10 +52,10 @@ void scene_object_base::update_tweens()
 	transform.update();
 }
 
-void scene_object_base::look_at(const float3& position, const float3& target, const float3& up)
+void scene_object_base::look_at(const vector_type& position, const vector_type& target, const vector_type& up)
 {
 	transform[1].translation = position;
-	transform[1].rotation = vmq::look_rotation(vmq::normalize(vmq::sub(target, position)), up);
+	transform[1].rotation = math::look_rotation(math::normalize(math::sub(target, position)), up);
 	transformed();
 }
 

src/scene/scene-object.hpp

@@ -20,15 +20,13 @@
 #ifndef ANTKEEPER_SCENE_OBJECT_HPP
 #define ANTKEEPER_SCENE_OBJECT_HPP
 
-#include <atomic>
-#include <cstdlib>
-#include <vmq/vmq.hpp>
 #include "animation/tween.hpp"
 #include "geometry/bounding-volume.hpp"
-
-using namespace vmq::types;
-using vmq::quaternion;
-using vmq::transform;
+#include "math/vector-type.hpp"
+#include "math/quaternion-type.hpp"
+#include "math/transform-type.hpp"
+#include <atomic>
+#include <cstdlib>
 
 /**
  * Internal base class for scene objects.
@@ -36,6 +34,11 @@ using vmq::transform;
 class scene_object_base
 {
 public:
+	typedef math::vector<float, 3> vector_type;
+	typedef math::quaternion<float> quaternion_type;
+	typedef math::transform<float> transform_type;
+	typedef bounding_volume<float> bounding_volume_type;
+	
 	/// Returns the type ID for this scene object type.
 	virtual const std::size_t get_object_type_id() const = 0;
 
@@ -62,32 +65,32 @@ public:
 	/**
 	 *
 	 */
-	void look_at(const float3& position, const float3& target, const float3& up);
+	void look_at(const vector_type& position, const vector_type& target, const vector_type& up);
 
 	/**
 	 * Sets the scene object's transform.
 	 */
-	void set_transform(const transform<float>& transform);
+	void set_transform(const transform_type& transform);
 
 	/**
 	 * Sets the scene object's translation.
 	 */
-	void set_translation(const float3& translation);
+	void set_translation(const vector_type& translation);
 	
 	/**
 	 * Sets the scene object's rotation.
 	 */
-	void set_rotation(const quaternion<float>& rotation);
+	void set_rotation(const quaternion_type& rotation);
 
 	/**
 	 * Sets the scene object's scale.
 	 */
-	void set_scale(const float3& scale);
+	void set_scale(const vector_type& scale);
 	
 	/**
 	 * Sets a culling mask for the object, which will be used for view-frustum culling instead of the object's bounds.
 	 */
-	void set_culling_mask(const bounding_volume<float>* culling_mask);
+	void set_culling_mask(const bounding_volume_type* culling_mask);
 	
 	/// Returns whether the scene object is active.
 	bool is_active() const;
@@ -95,51 +98,54 @@ public:
 	/**
 	 * Returns the transform.
 	 */
-	const transform<float>& get_transform() const;
+	const transform_type& get_transform() const;
 
 	/**
 	 * Returns the transform's translation vector.
 	 */
-	const float3& get_translation() const;
+	const vector_type& get_translation() const;
 
 	/**
 	 * Returns the transform's rotation quaternion.
 	 */
-	const quaternion<float>& get_rotation() const;
+	const quaternion_type& get_rotation() const;
 
 	/**
 	 * Returns the transform's scale vector.
 	 */
-	const float3& get_scale() const;
+	const vector_type& get_scale() const;
 
 	/**
 	 * Returns the transform tween.
 	 */
-	const tween<transform<float>>& get_transform_tween() const;
-	tween<transform<float>>& get_transform_tween();
+	const tween<transform_type>& get_transform_tween() const;
+	tween<transform_type>& get_transform_tween();
 
 	/**
 	 * Returns the bounds of the object.
 	 */
-	virtual const bounding_volume<float>& get_bounds() const = 0;
+	virtual const bounding_volume_type& get_bounds() const = 0;
 	
 	/**
 	 * Returns the culling mask of the object.
 	 */
-	const bounding_volume<float>* get_culling_mask() const;
+	const bounding_volume_type* get_culling_mask() const;
 
 protected:
 	static std::size_t next_object_type_id();
 
 private:
+	/// Interpolates between two transforms.
+	static transform_type interpolate_transforms(const transform_type& x, const transform_type& y, float a);
+	
 	/**
 	 * Called every time the scene object's tranform is changed.
 	 */
 	virtual void transformed();
 
 	bool active;
-	tween<transform<float>> transform;
-	const bounding_volume<float>* culling_mask;
+	tween<transform_type> transform;
+	const bounding_volume_type* culling_mask;
 };
 
 inline void scene_object_base::set_active(bool active)
@@ -147,25 +153,25 @@ inline void scene_object_base::set_active(bool active)
 	this->active = active;
 }
 
-inline void scene_object_base::set_transform(const ::transform<float>& transform)
+inline void scene_object_base::set_transform(const transform_type& transform)
 {
 	this->transform[1] = transform;
 	transformed();
 }
 
-inline void scene_object_base::set_translation(const float3& translation)
+inline void scene_object_base::set_translation(const vector_type& translation)
 {
 	transform[1].translation = translation;
 	transformed();
 }
 
-inline void scene_object_base::set_rotation(const quaternion<float>& rotation)
+inline void scene_object_base::set_rotation(const quaternion_type& rotation)
 {
 	transform[1].rotation = rotation;
 	transformed();
 }
 
-inline void scene_object_base::set_scale(const float3& scale)
+inline void scene_object_base::set_scale(const vector_type& scale)
 {
 	transform[1].scale = scale;
 	transformed();
@@ -176,37 +182,37 @@ inline bool scene_object_base::is_active() const
 	return active;
 }
 
-inline const transform<float>& scene_object_base::get_transform() const
+inline const typename scene_object_base::transform_type& scene_object_base::get_transform() const
 {
 	return transform[1];
 }
 
-inline const float3& scene_object_base::get_translation() const
+inline const typename scene_object_base::vector_type& scene_object_base::get_translation() const
 {
 	return get_transform().translation;
 }
 
-inline const quaternion<float>& scene_object_base::get_rotation() const
+inline const typename scene_object_base::quaternion_type& scene_object_base::get_rotation() const
 {
 	return get_transform().rotation;
 }
 
-inline const float3& scene_object_base::get_scale() const
+inline const typename scene_object_base::vector_type& scene_object_base::get_scale() const
 {
 	return get_transform().scale;
 }
 
-inline const tween<transform<float>>& scene_object_base::get_transform_tween() const
+inline const tween<typename scene_object_base::transform_type>& scene_object_base::get_transform_tween() const
 {
 	return transform;
 }
 
-inline tween<transform<float>>& scene_object_base::get_transform_tween()
+inline tween<typename scene_object_base::transform_type>& scene_object_base::get_transform_tween()
 {
 	return transform;
 }
 
-inline const bounding_volume<float>* scene_object_base::get_culling_mask() const
+inline const typename scene_object_base::bounding_volume_type* scene_object_base::get_culling_mask() const
 {
 	return culling_mask;
 }

src/scene/spotlight.cpp

@@ -19,23 +19,21 @@
 
 #include "spotlight.hpp"
 #include "configuration.hpp"
-#include "animation/ease.hpp"
+#include "math/math.hpp"
 #include <cmath>
 
-using namespace vmq::operators;
-
 static float3 interpolate_direction(const float3& x, const float3& y, float a)
 {
-	quaternion<float> q0 = vmq::rotation(global_forward, x);
-	quaternion<float> q1 = vmq::rotation(global_forward, y);
-	return vmq::normalize(vmq::slerp(q0, q1, a) * global_forward);
+	math::quaternion<float> q0 = math::rotation(global_forward, x);
+	math::quaternion<float> q1 = math::rotation(global_forward, y);
+	return math::normalize(math::slerp(q0, q1, a) * global_forward);
 }
 
 spotlight::spotlight():
 	direction(global_forward, interpolate_direction),
-	attenuation(float3{1, 0, 0}, ease<float3>::linear),
-	cutoff(float2{vmq::pi<float>, vmq::pi<float>}, ease<float2>::linear),
-	cosine_cutoff(float2{std::cos(vmq::pi<float>), std::cos(vmq::pi<float>)}, ease<float2>::linear)
+	attenuation(float3{1, 0, 0}, math::lerp<float3, float>),
+	cutoff(float2{math::pi<float>, math::pi<float>}, math::lerp<float2, float>),
+	cosine_cutoff(float2{std::cos(math::pi<float>), std::cos(math::pi<float>)}, math::lerp<float2, float>)
 {}
 
 void spotlight::set_attenuation(const float3& attenuation)
@@ -60,6 +58,6 @@ void spotlight::update_tweens()
 
 void spotlight::transformed()
 {
-	direction[1] = vmq::normalize(get_transform().rotation * global_forward);
+	direction[1] = math::normalize(get_transform().rotation * global_forward);
 }
 

src/scene/spotlight.hpp

@@ -21,9 +21,7 @@
 #define ANTKEEPER_SPOTLIGHT_HPP
 
 #include "scene/light.hpp"
-#include <vmq/vmq.hpp>
-
-using namespace vmq::types;
+#include "utility/fundamental-types.hpp"
 
 class spotlight: public light
 {

src/state/play-state.cpp

@@ -39,13 +39,11 @@
 #include "entity/archetype.hpp"
 #include "entity/entity-commands.hpp"
 #include "nest.hpp"
-#include "math.hpp"
+#include "math/math.hpp"
+#include "utility/fundamental-types.hpp"
 #include "geometry/mesh-accelerator.hpp"
 #include "behavior/ebt.hpp"
 #include "animation/ease.hpp"
-#include <iostream>
-
-using namespace vmq::operators;
 
 void enter_play_state(application* app)
 {
@@ -92,15 +90,15 @@ void enter_play_state(application* app)
 	
 	for (int i = 0; i < pebble_count; ++i)
 	{
-		float x = frand(-pebble_radius, pebble_radius);
-		float z = frand(-pebble_radius, pebble_radius);
+		float x = math::random(-pebble_radius, pebble_radius);
+		float z = math::random(-pebble_radius, pebble_radius);
 		
 		auto pebble_entity = pebble_archetype->create(ecs_registry);
 		
 		auto& transform = ecs_registry.get<ecs::transform_component>(pebble_entity);
-		transform.transform = vmq::identity_transform<float>;
-		transform.transform.rotation = vmq::angle_axis(frand(0.0f, vmq::two_pi<float>), {0, 1, 0});
-		transform.transform.scale = float3{1, 1, 1} * frand(0.75f, 1.25f);
+		transform.transform = math::identity_transform<float>;
+		transform.transform.rotation = math::angle_axis(math::random(0.0f, math::two_pi<float>), {0, 1, 0});
+		transform.transform.scale = float3{1, 1, 1} * math::random(0.75f, 1.25f);
 		
 		placement.ray.origin = {x, 10000, z};
 		ecs_registry.assign<ecs::placement_component>(pebble_entity, placement);
@@ -133,15 +131,15 @@ void enter_play_state(application* app)
 
 		auto& transform = ecs_registry.get<ecs::transform_component>(samara_entity);
 		float zone = 200.0f;
-		transform.transform = vmq::identity_transform<float>;
-		transform.transform.translation.x = frand(-zone, zone);
-		transform.transform.translation.y = frand(50.0f, 150.0f);
-		transform.transform.translation.z = frand(-zone, zone);
+		transform.transform = math::identity_transform<float>;
+		transform.transform.translation.x = math::random(-zone, zone);
+		transform.transform.translation.y = math::random(50.0f, 150.0f);
+		transform.transform.translation.z = math::random(-zone, zone);
 
 		ecs::samara_component samara_component;
-		samara_component.angle = frand(0.0f, vmq::radians(360.0f));
-		samara_component.direction = vmq::normalize(float3{frand(-1, 1), frand(-1, -5), frand(-1, 1)});
-		samara_component.chirality = (frand(0, 1) < 0.5f) ? -1.0f : 1.0f;
+		samara_component.angle = math::random(0.0f, math::radians(360.0f));
+		samara_component.direction = math::normalize(float3{math::random(-1.0f, 1.0f), math::random(-1.0f, -5.0f), math::random(-1.0f, 1.0f)});
+		samara_component.chirality = (math::random(0.0f, 1.0f) < 0.5f) ? -1.0f : 1.0f;
 
 		ecs_registry.assign_or_replace<ecs::samara_component>(samara_entity, samara_component);
 	}
@@ -150,7 +148,7 @@ void enter_play_state(application* app)
 	ecs::archetype* grass_archetype = resource_manager->load<ecs::archetype>("grassland-grass.ent");
 	auto grass_entity_1 = grass_archetype->create(ecs_registry);
 	auto grass_entity_2 = grass_archetype->create(ecs_registry);
-	ecs_registry.get<ecs::transform_component>(grass_entity_2).transform.rotation = vmq::angle_axis(vmq::radians(120.0f), float3{0, 1, 0});
+	ecs_registry.get<ecs::transform_component>(grass_entity_2).transform.rotation = math::angle_axis(math::radians(120.0f), float3{0, 1, 0});
 	*/
 
 	// Setup overworld camera
@@ -159,7 +157,7 @@ void enter_play_state(application* app)
 	orbit_cam->attach(camera);
 	orbit_cam->set_target_focal_point({0, 0, 0});
 	orbit_cam->set_target_focal_distance(15.0f);
-	orbit_cam->set_target_elevation(vmq::radians(25.0f));
+	orbit_cam->set_target_elevation(math::radians(25.0f));
 	orbit_cam->set_target_azimuth(0.0f);
 	orbit_cam->set_focal_point(orbit_cam->get_target_focal_point());
 	orbit_cam->set_focal_distance(orbit_cam->get_target_focal_distance());
@@ -185,7 +183,7 @@ void enter_play_state(application* app)
 	nest->set_tunnel_radius(tunnel_radius);
 	nest::shaft* central_shaft = nest->get_central_shaft();
 	central_shaft->chirality = 1.0f;
-	central_shaft->rotation = vmq::radians(0.0f);
+	central_shaft->rotation = math::radians(0.0f);
 	central_shaft->depth = {0.0f, 200.0f};
 	central_shaft->radius = {15.0f, 15.0f};
 	central_shaft->pitch = {40.0f, 40.0f};
@@ -196,7 +194,7 @@ void enter_play_state(application* app)
 		nest::chamber chamber;
 		chamber.shaft = central_shaft;
 		chamber.depth = (i + 1) * 50.0f;
-		chamber.rotation = vmq::radians(0.0f);
+		chamber.rotation = math::radians(0.0f);
 		chamber.inner_radius = 4.0f;
 		chamber.outer_radius = 10.0f;
 		central_shaft->chambers.push_back(chamber);
@@ -208,8 +206,8 @@ void enter_play_state(application* app)
 	{
 		ecs::cavity_component cavity;
 		cavity.position = nest->extend_shaft(*nest->get_central_shaft());
-		cavity.position += float3{frand(-shift, shift), frand(-shift, shift), frand(-shift, shift)};
-		cavity.radius = tunnel_radius * frand(1.0f, 1.1f);
+		cavity.position += float3{math::random(-shift, shift), math::random(-shift, shift), math::random(-shift, shift)};
+		cavity.radius = tunnel_radius * math::random(1.0f, 1.1f);
 
 		ecs_registry.assign<ecs::cavity_component>(ecs_registry.create(), cavity);
 	}
@@ -222,8 +220,8 @@ void enter_play_state(application* app)
 		{
 			ecs::cavity_component cavity;
 			cavity.position = nest->expand_chamber(central_shaft->chambers[i]);
-			cavity.position += float3{frand(-shift, shift), frand(-shift, shift), frand(-shift, shift)};
-			cavity.radius = tunnel_radius * frand(1.0f, 1.1f);
+			cavity.position += float3{math::random(-shift, shift), math::random(-shift, shift), math::random(-shift, shift)};
+			cavity.radius = tunnel_radius * math::random(1.0f, 1.1f);
 
 			ecs_registry.assign<ecs::cavity_component>(ecs_registry.create(), cavity);
 		}
@@ -236,7 +234,7 @@ void enter_play_state(application* app)
 		ecs::assign_render_layers(ecs_registry, larva, 1);
 		//ecs::warp_to(ecs_registry, larva, {0, -20, 0});
 		//auto& transform = ecs_registry.get<ecs::transform_component>(larva_entity);
-		//transform.transform = vmq::identity_transform<float>;
+		//transform.transform = math::identity_transform<float>;
 		//transform.transform.translation = nest->get_shaft_position(*central_shaft, central_shaft->depth[1]);
 		//transform.transform.translation.y -= 1.0f;
 	}

src/systems/camera-system.cpp

@@ -25,9 +25,8 @@
 #include "orbit-cam.hpp"
 #include "geometry/mesh.hpp"
 #include "geometry/intersection.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 
-using namespace vmq::operators;
 using namespace ecs;
 
 camera_system::camera_system(entt::registry& registry):
@@ -74,9 +73,9 @@ void camera_system::update(double t, double dt)
 	registry.view<transform_component, collision_component>().each(
 		[&](auto entity, auto& transform, auto& collision)
 		{
-			vmq::transform<float> inverse_transform = vmq::inverse(transform.transform);
+			math::transform<float> inverse_transform = math::inverse(transform.transform);
 			float3 origin = inverse_transform * pick_origin;
-			float3 direction = vmq::normalize(vmq::conjugate(transform.transform.rotation) * pick_direction);
+			float3 direction = math::normalize(math::conjugate(transform.transform.rotation) * pick_direction);
 			ray<float> transformed_ray = {origin, direction};
 
 			// Broad phase AABB test

src/systems/camera-system.hpp

@@ -23,8 +23,7 @@
 #include "entity-system.hpp"
 #include "event/event-handler.hpp"
 #include "input/input-events.hpp"
-#include <vmq/vmq.hpp>
-using namespace vmq::types;
+#include "utility/fundamental-types.hpp"
 
 class camera;
 class orbit_cam;

src/systems/constraint-system.cpp

@@ -23,9 +23,9 @@
 #include "entity/components/copy-scale-component.hpp"
 #include "entity/components/copy-transform-component.hpp"
 #include "entity/components/transform-component.hpp"
+#include "utility/fundamental-types.hpp"
 
 using namespace ecs;
-using namespace vmq;
 
 constraint_system::constraint_system(entt::registry& registry):
 	entity_system(registry)

src/systems/control-system.cpp

@@ -24,9 +24,7 @@
 #include "geometry/intersection.hpp"
 #include "animation/ease.hpp"
 #include "nest.hpp"
-#include <vmq/vmq.hpp>
-
-using namespace vmq::operators;
+#include "math/math.hpp"
 
 control_system::control_system():
 	timestep(0.0f),
@@ -60,14 +58,14 @@ control_system::control_system():
 	}
 
 	zoom_speed = 4.0f; //1
-	min_elevation = vmq::radians(-85.0f);
-	max_elevation = vmq::radians(85.0f);
+	min_elevation = math::radians(-85.0f);
+	max_elevation = math::radians(85.0f);
 	near_focal_distance = 2.0f;
 	far_focal_distance = 200.0f;
 	near_movement_speed = 10.0f;
 	far_movement_speed = 80.0f;
-	near_fov = vmq::radians(80.0f);
-	far_fov = vmq::radians(35.0f);
+	near_fov = math::radians(80.0f);
+	far_fov = math::radians(35.0f);
 	near_clip_near = 0.1f;
 	far_clip_near = 5.0f;
 	near_clip_far = 100.0f;
@@ -98,15 +96,15 @@ void control_system::update(float dt)
 		zoom -= zoom_speed * dt * zoom_out_control.get_current_value();
 	zoom = std::max<float>(0.0f, std::min<float>(1.0f, zoom));
 
-	float focal_distance = ease<float>::log(near_focal_distance, far_focal_distance, 1.0f - zoom);
+	float focal_distance = math::log_lerp<float>(near_focal_distance, far_focal_distance, 1.0f - zoom);
 
-	float fov = ease<float>::log(near_fov, far_fov, 1.0f - zoom);
+	float fov = math::log_lerp<float>(near_fov, far_fov, 1.0f - zoom);
 
 	//float elevation_factor = (orbit_cam->get_target_elevation() - min_elevation) / max_elevation;
-	//fov = ease<float>::log(near_fov, far_fov, elevation_factor);
-	float clip_near = ease<float>::log(near_clip_near, far_clip_near, 1.0f - zoom);
-	float clip_far = ease<float>::log(near_clip_far, far_clip_far, 1.0f - zoom);
-	float movement_speed = ease<float>::log(near_movement_speed * dt, far_movement_speed * dt, 1.0f - zoom);
+	//fov = math::log_lerp<float>(near_fov, far_fov, elevation_factor);
+	float clip_near = math::log_lerp<float>(near_clip_near, far_clip_near, 1.0f - zoom);
+	float clip_far = math::log_lerp<float>(near_clip_far, far_clip_far, 1.0f - zoom);
+	float movement_speed = math::log_lerp<float>(near_movement_speed * dt, far_movement_speed * dt, 1.0f - zoom);
 
 	orbit_cam->set_target_focal_distance(focal_distance);
 	orbit_cam->get_camera()->set_perspective(fov, orbit_cam->get_camera()->get_aspect_ratio(), clip_near, clip_far);
@@ -147,13 +145,13 @@ void control_system::update(float dt)
 		movement[1] += move_back_control.get_current_value();
 	
 	const float deadzone = 0.01f;
-	float magnitude_squared = vmq::length_squared(movement);
+	float magnitude_squared = math::length_squared(movement);
 
 	if (magnitude_squared > deadzone)
 	{
 		if (magnitude_squared > 1.0f)
 		{
-			movement = vmq::normalize(movement);
+			movement = math::normalize(movement);
 		}
 
 		orbit_cam->move(movement * movement_speed);
@@ -189,7 +187,7 @@ void control_system::update(float dt)
 
 	float3 pick_near = camera->unproject({mouse_position[0], viewport[3] - mouse_position[1], 0.0f}, viewport);
 	float3 pick_far = camera->unproject({mouse_position[0], viewport[3] - mouse_position[1], 1.0f}, viewport);
-	float3 pick_direction = vmq::normalize(pick_far - pick_near);
+	float3 pick_direction = math::normalize(pick_far - pick_near);
 	ray<float> picking_ray = {pick_near, pick_direction};
 	plane<float> ground_plane = {float3{0, 1, 0}, 0.0f};
 
@@ -210,7 +208,7 @@ void control_system::update(float dt)
 	
 	
 	float2 viewport_center = {(viewport[0] + viewport[2]) * 0.5f, (viewport[1] + viewport[3]) * 0.5f};
-	float2 mouse_direction = vmq::normalize(mouse_position - viewport_center);
+	float2 mouse_direction = math::normalize(mouse_position - viewport_center);
 	old_mouse_angle = mouse_angle;
 	mouse_angle = std::atan2(-mouse_direction.y, mouse_direction.x);
 	
@@ -218,22 +216,22 @@ void control_system::update(float dt)
 	{
 
 		
-		if (mouse_angle - old_mouse_angle <= -vmq::pi<float>)
+		if (mouse_angle - old_mouse_angle <= -math::pi<float>)
 			flashlight_turns_i -= 1;
-		else if (mouse_angle - old_mouse_angle >= vmq::pi<float>)
+		else if (mouse_angle - old_mouse_angle >= math::pi<float>)
 			flashlight_turns_i += 1;
 		
-		flashlight_turns_f = (mouse_angle) / vmq::two_pi<float>;
+		flashlight_turns_f = (mouse_angle) / math::two_pi<float>;
 		flashlight_turns = flashlight_turns_i - flashlight_turns_f;
 		
 		if (flashlight && nest)
 		{
-			transform<float> flashlight_transform = vmq::identity_transform<float>;
+			math::transform<float> flashlight_transform = math::identity_transform<float>;
 			
 			float flashlight_depth = nest->get_shaft_depth(*nest->get_central_shaft(), flashlight_turns);
 			
 			flashlight_transform.translation = {0.0f, -flashlight_depth, 0.0f};
-			flashlight_transform.rotation = vmq::angle_axis(-flashlight_turns * vmq::two_pi<float> + vmq::half_pi<float>, {0, 1, 0});
+			flashlight_transform.rotation = math::angle_axis(-flashlight_turns * math::two_pi<float> + math::half_pi<float>, {0, 1, 0});
 			
 			flashlight->set_transform(flashlight_transform);
 			flashlight_light_cone->set_transform(flashlight_transform);
@@ -298,7 +296,7 @@ void control_system::handle_event(const mouse_moved_event& event)
 			elevation_factor *= -1.0f;
 		}
 
-		float rotation = vmq::radians(22.5f) * rotation_factor * timestep;
+		float rotation = math::radians(22.5f) * rotation_factor * timestep;
 		float elevation = orbit_cam->get_target_elevation() + elevation_factor * 0.25f * timestep;
 		elevation = std::min<float>(max_elevation, std::max<float>(min_elevation, elevation));
 

src/systems/control-system.hpp

@@ -25,6 +25,7 @@
 #include "input/control.hpp"
 #include "input/control-set.hpp"
 #include "scene/model-instance.hpp"
+#include "utility/fundamental-types.hpp"
 
 class orbit_cam;
 class nest;

src/systems/nest-system.cpp

@@ -19,6 +19,7 @@
 
 #include "nest-system.hpp"
 #include "nest.hpp"
+#include "math/math.hpp"
 
 using namespace ecs;
 
@@ -45,7 +46,7 @@ void nest_system::on_nest_construct(entt::registry& registry, entt::entity entit
 	nest->set_tunnel_radius(1.15f);
 	nest::shaft* central_shaft = nest->get_central_shaft();
 	central_shaft->chirality = -1.0f;
-	central_shaft->rotation = vmq::radians(0.0f);
+	central_shaft->rotation = math::radians(0.0f);
 	central_shaft->depth = {0.0f, 100.0f};
 	central_shaft->current_depth = 0.0f;
 	central_shaft->radius = {0.0f, 5.0f};
@@ -56,7 +57,7 @@ void nest_system::on_nest_construct(entt::registry& registry, entt::entity entit
 		nest::chamber chamber;
 		chamber.shaft = central_shaft;
 		chamber.depth = (i + 1) * 23.0f;
-		chamber.rotation = vmq::radians(0.0f);
+		chamber.rotation = math::radians(0.0f);
 		chamber.inner_radius = 4.0f;
 		chamber.outer_radius = 10.0f;
 		central_shaft->chambers.push_back(chamber);

src/systems/placement-system.cpp

@@ -22,6 +22,7 @@
 #include "entity/components/placement-component.hpp"
 #include "entity/components/transform-component.hpp"
 #include "entity/components/terrain-component.hpp"
+#include "utility/fundamental-types.hpp"
 
 using namespace ecs;
 
@@ -42,9 +43,9 @@ void placement_system::update(double t, double dt)
 				[&](auto entity, auto& transform, auto& collision)
 				{
 					// Transform ray into local space of collision component
-					vmq::transform<float> inverse_transform = vmq::inverse(transform.transform);
+					math::transform<float> inverse_transform = math::inverse(transform.transform);
 					float3 origin = inverse_transform * placement.ray.origin;
-					float3 direction = vmq::normalize(vmq::conjugate(transform.transform.rotation) * placement.ray.direction);
+					float3 direction = math::normalize(math::conjugate(transform.transform.rotation) * placement.ray.direction);
 					ray<float> transformed_ray = {origin, direction};
 
 					// Broad phase AABB test

src/systems/samara-system.cpp

@@ -20,9 +20,9 @@
 #include "samara-system.hpp"
 #include "entity/components/transform-component.hpp"
 #include "entity/components/samara-component.hpp"
-#include "math.hpp"
+#include "math/math.hpp"
+#include "utility/fundamental-types.hpp"
 
-using namespace vmq::operators;
 using namespace ecs;
 
 samara_system::samara_system(entt::registry& registry):
@@ -34,23 +34,23 @@ void samara_system::update(double t, double dt)
 	registry.view<samara_component, transform_component>().each(
 		[&](auto entity, auto& samara, auto& transform)
 		{
-			samara.angle += samara.chirality * vmq::radians(360.0f * 6.0f) * dt;
+			samara.angle += samara.chirality * math::radians(360.0f * 6.0f) * dt;
 
 			transform.transform.translation += samara.direction * 20.0f * (float)dt;
 			transform.transform.rotation =
-				vmq::angle_axis(samara.angle, float3{0, 1, 0}) *
-				vmq::angle_axis(vmq::radians(20.0f), float3{1, 0, 0}) *
-				((samara.chirality < 0.0f) ? vmq::angle_axis(vmq::radians(180.0f), float3{0, 0, -1}) : vmq::quaternion<float>{1, 0, 0, 0});
+				math::angle_axis(samara.angle, float3{0, 1, 0}) *
+				math::angle_axis(math::radians(20.0f), float3{1, 0, 0}) *
+				((samara.chirality < 0.0f) ? math::angle_axis(math::radians(180.0f), float3{0, 0, -1}) : math::quaternion<float>{1, 0, 0, 0});
 
 			if (transform.transform.translation.y < 0.0f)
 			{
 				const float zone = 200.0f;
-				transform.transform.translation.x = frand(-zone, zone);
-				transform.transform.translation.y = frand(100.0f, 150.0f);
-				transform.transform.translation.z = frand(-zone, zone);
+				transform.transform.translation.x = math::random(-zone, zone);
+				transform.transform.translation.y = math::random(100.0f, 150.0f);
+				transform.transform.translation.z = math::random(-zone, zone);
 				transform.warp = true;
 
-				samara.chirality = (frand(0, 1) < 0.5f) ? -1.0f : 1.0f;
+				samara.chirality = (math::random(0.0f, 1.0f) < 0.5f) ? -1.0f : 1.0f;
 			}
 		});
 }

src/systems/subterrain-system.cpp

@@ -32,15 +32,13 @@
 #include "geometry/intersection.hpp"
 #include "scene/scene.hpp"
 #include "scene/model-instance.hpp"
-#include "math.hpp"
-
-using namespace vmq::operators;
-using namespace ecs;
-
+#include "utility/fundamental-types.hpp"
 #include <array>
 #include <iostream>
 #include <limits>
 
+using namespace ecs;
+
 /**
  * An octree containing cubes for the marching cubes algorithm.
  */
@@ -422,7 +420,7 @@ void subterrain_system::regenerate_subterrain_model()
 				edge = edge->previous->symmetric;
 			}
 			while (edge != start);
-			n = vmq::normalize(n);
+			n = math::normalize(n);
 
 			//float3 n = reinterpret_cast<const float3&>(face_normals[i * 3]);
 
@@ -474,10 +472,10 @@ void subterrain_system::dig(const float3& position, float radius)
 			for (int i = 0; i < 8; ++i)
 			{
 				// For outside normals (also set node initial distance to +infinity)
-				//float distance = vmq::length(node->corners[i] - position) - radius;
+				//float distance = math::length(node->corners[i] - position) - radius;
 				// if (distance < node->distances[i])
 
-				float distance = radius - vmq::length(node.corners[i] - position);
+				float distance = radius - math::length(node.corners[i] - position);
 				if (distance > node.distances[i])
 					node.distances[i] = distance;
 			}

src/systems/subterrain-system.hpp

@@ -23,6 +23,7 @@
 #include "entity-system.hpp"
 #include "geometry/mesh.hpp"
 #include "geometry/aabb.hpp"
+#include "utility/fundamental-types.hpp"
 #include <unordered_map>
 
 class resource_manager;
@@ -33,8 +34,6 @@ struct cube_tree;
 class scene;
 class model_instance;
 
-using namespace vmq::types;
-
 template <std::int64_t Mantissa, std::int64_t Exponent>
 struct epsilon
 {
@@ -49,7 +48,7 @@ typedef epsilon<1, -5> epsilon_1en5;
 template <class Epsilon, class T, std::size_t N>
 struct vector_hasher
 {
-	typedef vmq::vector<T, N> vector_type;
+	typedef math::vector<T, N> vector_type;
 
 	std::size_t operator()(const vector_type& v) const noexcept
 	{
@@ -69,7 +68,7 @@ struct vector_hasher
 template <class Epsilon, class T, std::size_t N>
 struct vector_equals
 {
-	typedef vmq::vector<T, N> vector_type;
+	typedef math::vector<T, N> vector_type;
 
 	bool operator()(const vector_type& a, const vector_type& b) const noexcept
 	{

src/systems/terrain-system.cpp

@@ -30,9 +30,9 @@
 #include "rasterizer/vertex-buffer.hpp"
 #include "resources/resource-manager.hpp"
 #include "resources/image.hpp"
+#include "utility/fundamental-types.hpp"
 #include <limits>
 
-using namespace vmq::operators;
 using namespace ecs;
 
 terrain_system::terrain_system(entt::registry& registry, ::resource_manager* resource_manager):
@@ -275,7 +275,7 @@ void terrain_system::update_terrain_model(model* terrain_model, mesh* terrain_me
 				edge = edge->previous->symmetric;
 			}
 			while (edge != start);
-			n = vmq::normalize(n);
+			n = math::normalize(n);
 
 			*(v++) = vertex->position[0];
 			*(v++) = vertex->position[1];
@@ -330,7 +330,7 @@ void terrain_system::on_terrain_construct(entt::registry& registry, entt::entity
 
 	// Assign the entity a transform component
 	transform_component transform;
-	transform.transform = vmq::identity_transform<float>;
+	transform.transform = math::identity_transform<float>;
 	transform.transform.translation = float3{(float)component.x * patch_size, 0.0f, (float)component.z * patch_size};
 	transform.warp = true;
 	registry.assign_or_replace<transform_component>(entity, transform);

src/systems/tool-system.cpp

@@ -25,9 +25,8 @@
 #include "orbit-cam.hpp"
 #include "geometry/mesh.hpp"
 #include "geometry/intersection.hpp"
-#include <vmq/vmq.hpp>
+#include "math/math.hpp"
 
-using namespace vmq::operators;
 using namespace ecs;
 
 tool_system::tool_system(entt::registry& registry):
@@ -47,7 +46,7 @@ void tool_system::update(double t, double dt)
 	float3 pick_near = camera->unproject({mouse_position[0], viewport[3] - mouse_position[1], 0.0f}, viewport);
 	float3 pick_far = camera->unproject({mouse_position[0], viewport[3] - mouse_position[1], 1.0f}, viewport);
 	float3 pick_origin = pick_near;
-	float3 pick_direction = vmq::normalize(pick_far - pick_near);
+	float3 pick_direction = math::normalize(pick_far - pick_near);
 	ray<float> picking_ray = {pick_near, pick_direction};
 
 	float a = std::numeric_limits<float>::infinity();
@@ -58,9 +57,9 @@ void tool_system::update(double t, double dt)
 	registry.view<transform_component, collision_component>().each(
 		[&](auto entity, auto& transform, auto& collision)
 		{
-			vmq::transform<float> inverse_transform = vmq::inverse(transform.transform);
+			math::transform<float> inverse_transform = math::inverse(transform.transform);
 			float3 origin = inverse_transform * pick_origin;
-			float3 direction = vmq::normalize(vmq::conjugate(transform.transform.rotation) * pick_direction);
+			float3 direction = math::normalize(math::conjugate(transform.transform.rotation) * pick_direction);
 			ray<float> transformed_ray = {origin, direction};
 
 			// Broad phase AABB test
@@ -88,13 +87,13 @@ void tool_system::update(double t, double dt)
 	float3 camera_planar_position = float3{camera_position.x, 0, camera_position.z};
 
 	float pick_angle = 0.0f;
-	float3 pick_planar_direction = vmq::normalize(pick_planar_position - camera_planar_position);
+	float3 pick_planar_direction = math::normalize(pick_planar_position - camera_planar_position);
 	float3 camera_planar_focal_point = float3{orbit_cam->get_focal_point().x, 0, orbit_cam->get_focal_point().z};
-	float3 camera_planar_direction = vmq::normalize(camera_planar_focal_point - camera_planar_position);
-	if (std::fabs(vmq::length_squared(camera_planar_direction - pick_planar_direction) > 0.0001f))
+	float3 camera_planar_direction = math::normalize(camera_planar_focal_point - camera_planar_position);
+	if (std::fabs(math::length_squared(camera_planar_direction - pick_planar_direction) > 0.0001f))
 	{
-		pick_angle = std::acos(vmq::dot(camera_planar_direction, pick_planar_direction));
-		if (vmq::dot(vmq::cross(camera_planar_direction, pick_planar_direction), float3{0, 1, 0}) < 0.0f)
+		pick_angle = std::acos(math::dot(camera_planar_direction, pick_planar_direction));
+		if (math::dot(math::cross(camera_planar_direction, pick_planar_direction), float3{0, 1, 0}) < 0.0f)
 			pick_angle = -pick_angle;
 	}
 
@@ -110,7 +109,7 @@ void tool_system::update(double t, double dt)
 				transform.transform.translation = pick;
 			}
 
-			vmq::quaternion<float> rotation = vmq::angle_axis(orbit_cam->get_azimuth() + pick_angle, float3{0, 1, 0});
+			math::quaternion<float> rotation = math::angle_axis(orbit_cam->get_azimuth() + pick_angle, float3{0, 1, 0});
 			transform.transform.rotation = rotation;
 		});
 }
@@ -143,4 +142,3 @@ void tool_system::handle_event(const mouse_moved_event& event)
 		mouse_position[1] = event.y;
 	}
 }
-