CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 3.8)
+cmake_minimum_required(VERSION 3.10)
 project(psemek)
 
 cmake_policy(SET CMP0077 NEW)
@@ -45,8 +45,6 @@ option(PSEMEK_LEGACY_UI "Use legacy UI library" OFF)
 
 message(STATUS "Using backend ${PSEMEK_BACKEND}")
 
-option(PSEMEK_BT_LOG "Add logging for behavior tree actions" OFF)
-
 if(PSEMEK_BACKEND STREQUAL "ANDROID")
 	set(PSEMEK_GL_API gles32)
 	set(PSEMEK_GL_LIBRARIES GLESv3 EGL)
@@ -78,7 +76,7 @@ if(NOT DEFINED PSEMEK_GRAPHICS_API)
 endif()
 
 if(PSEMEK_GRAPHICS_API STREQUAL WEBGPU)
-	find_package(wgpu-native REQUIRED)
+	find_package(wgpu-native 27.0.2.0 REQUIRED)
 endif()
 
 message(STATUS "Using graphics API ${PSEMEK_GRAPHICS_API}")

LICENSE.txt

@@ -0,0 +1,29 @@
+MIT NON-AI NON-NFT License
+
+Copyright (c) 2025, Nikita Lisitsa
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of the software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions.
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+In addition, the following restrictions apply:
+
+1. The Software, any modifications made to it, any artifacts or media generated by it, and any parts thereof, including but not limited to images,
+video, and code snippets, may not be used for the purpose of training or improving machine learning algorithms, or be included in any dataset used
+for training or improving machine learning algorithms, including but not limited to artificial intelligence, natural language processing, large
+language models, or data mining. This condition applies to any derivatives, modifications, or updates based on the Software code. Any usage of the
+Software in an AI-training dataset is considered a breach of this License.
+
+2. The Software, any modifications made to it, any artifacts or media generated by it, and any parts thereof, including but not limited to images,
+video, and code snippets, may not be owned via a blockchain-backed ownership token, or as a means of generating blockchain-backed ownership tokens,
+including but not limited to non-fungible tokens. Any usage of the Software to generate non-fungible tokens is considered a breach of this License.
+
+3. Any person or organization found to be in violation of these restrictions will be subject to legal action and may be held liable for any damages
+resulting from such use.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
+OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

cmake/modules/Findwgpu-native.cmake

@@ -3,16 +3,27 @@ if(wgpu-native_FOUND)
 endif()
 
 # Don't search for include files - these are bundled with psemek-wgpu lib
-find_library(wgpu-native_LIBRARIES NAMES libwgpu_native.a wgpu_native.dll wgpu_native PATHS "${WGPU_NATIVE_ROOT}")
+find_library(wgpu-native_LIBRARY NAMES libwgpu_native.a wgpu_native.dll wgpu_native PATHS "${WGPU_NATIVE_ROOT}" "${WGPU_NATIVE_ROOT}/lib")
+find_file(wgpu-native_VERSION_FILE wgpu-native-git-tag PATHS "${WGPU_NATIVE_ROOT}/wgpu-native-meta")
+
+if(EXISTS "${wgpu-native_VERSION_FILE}")
+	file(READ "${wgpu-native_VERSION_FILE}" wgpu-native_VERSION_NOT_STRIPPED)
+	string(STRIP "${wgpu-native_VERSION_NOT_STRIPPED}" wgpu-native_VERSION_STRIPPED)
+	string(SUBSTRING "${wgpu-native_VERSION_STRIPPED}" 1 -1 wgpu-native_VERSION)
+endif()
 
 include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(wgpu-native DEFAULT_MSG wgpu-native_LIBRARIES)
+find_package_handle_standard_args(wgpu-native
+	REQUIRED_VARS wgpu-native_LIBRARY wgpu-native_VERSION_FILE
+	VERSION_VAR wgpu-native_VERSION
+)
 
 if(wgpu-native_FOUND AND NOT TARGET wgpu-native)
+	set(wgpu-native_LIBRARIES ${wgpu-native_LIBRARY})
 	add_library(wgpu-native STATIC IMPORTED)
 	set_target_properties(wgpu-native PROPERTIES
 		IMPORTED_LOCATION "${wgpu-native_LIBRARIES}"
 	)
 endif()
 
-mark_as_advanced(wgpu-native_LIBRARIES)
+mark_as_advanced(wgpu-native_LIBRARY wgpu-native_VERSION_FILE)

examples/cloud.cpp

@@ -223,7 +223,7 @@ struct cloud_app
 			random::generator rng(random::device{});
 			random::uniform_sphere_vector_distribution<float, 3> d;
 
-			std::vector<util::array<math::vector<float, 3>, 3>> grad(4);
+			std::vector<util::ndarray<math::vector<float, 3>, 3>> grad(4);
 			std::vector<float> weights(grad.size());
 			float weight_sum = 0.f;
 
@@ -247,7 +247,7 @@ struct cloud_app
 			// https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89062
 			math::gradient<float> g(std::make_pair(0.2f, 0.f), math::easing_type::quadratic_out, std::pair{0.3f, max_density});
 
-			util::array<std::uint8_t, 3> cloud_data({cloud_size[0], cloud_size[1], cloud_size[2]});
+			util::ndarray<std::uint8_t, 3> cloud_data({cloud_size[0], cloud_size[1], cloud_size[2]});
 
 			for (std::size_t z = 0; z < cloud_data.depth(); ++z)
 			{
@@ -328,7 +328,7 @@ struct cloud_app
 				return math::lerp(v0, v1, t[2]) / 255.f * max_density;
 			};
 
-			util::array<math::vector<float, 4>, 3> cloud_shadow_f(cloud_data.dims());
+			util::ndarray<math::vector<float, 4>, 3> cloud_shadow_f(cloud_data.dims());
 
 			dirs.resize(32);
 
@@ -448,7 +448,7 @@ struct cloud_app
 				}
 			}
 
-			util::array<math::vector<std::uint8_t, 4>, 3> cloud_shadow(cloud_data.dims());
+			util::ndarray<math::vector<std::uint8_t, 4>, 3> cloud_shadow(cloud_data.dims());
 
 			for (auto const & idx : cloud_shadow.indices())
 			{

examples/fibonacci_music_box.cpp

@@ -22,7 +22,7 @@
 #include <psemek/app/application_base.hpp>
 #include <psemek/app/default_application_factory.hpp>
 #include <psemek/log/log.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 
 #include <unordered_map>
 
@@ -32,7 +32,7 @@ using namespace psemek;
 
 std::vector<std::vector<math::point<int, 2>>> fibonacci_cycles(int n)
 {
-	util::array<int, 2> cycle({n, n}, -1);
+	util::ndarray<int, 2> cycle({n, n}, -1);
 
 	std::vector<std::vector<math::point<int, 2>>> result;
 

examples/fire.cpp

@@ -172,7 +172,7 @@ candle_renderer::candle_renderer()
 		random::generator rng;
 		random::uniform_sphere_vector_distribution<float, 2> d;
 
-		util::array<math::vector<float, 2>, 2> grad({16, 16});
+		util::ndarray<math::vector<float, 2>, 2> grad({16, 16});
 
 		for (auto & v : grad) v = d(rng);
 		pcg::perlin<float, 2> perlinx(grad, pcg::seamless);

examples/particle_life_2d.cpp

@@ -0,0 +1,441 @@
+#include <psemek/app/application_base.hpp>
+#include <psemek/app/default_application_factory.hpp>
+#include <psemek/gfx/painter.hpp>
+#include <psemek/gfx/gl.hpp>
+#include <psemek/math/camera.hpp>
+#include <psemek/math/contains.hpp>
+#include <psemek/cg/kdtree.hpp>
+#include <psemek/random/device.hpp>
+#include <psemek/random/generator.hpp>
+#include <psemek/random/uniform.hpp>
+#include <psemek/random/uniform_box.hpp>
+#include <psemek/util/clock.hpp>
+#include <psemek/util/ndarray.hpp>
+#include <psemek/log/log.hpp>
+
+#include <format>
+
+using namespace psemek;
+
+struct particle
+{
+	math::point<float, 2> position;
+	math::vector<float, 2> velocity;
+	float mass;
+	float radius;
+	int type;
+};
+
+struct particle_life_2d_app
+	: app::application_base
+{
+	particle_life_2d_app(options const &, context const &)
+	{
+		float S = 5.f;
+		float W = 160.f * S;
+		float H =  90.f * S;
+		area_ = {{{-W, W}, {-H, H}}};
+
+		types_ = 10;
+
+		auto seed = random::device{}();
+
+		// Explosions
+		// types_ = 6;
+		// seed = 3940378904;
+
+		// Cells
+		// types_ = 8;
+		// seed = 3071915692;
+
+		// Cells v2
+		// types_ = 6;
+		// seed = 3337663839;
+
+		// Cool cells!
+		types_ = 10;
+		seed = 388834085;
+
+		// One cell, use central gravity
+		// types_ = 10;
+		// seed = 175098945;
+
+		// Mitochondria cell
+		// types_ = 10;
+		// seed = 1763331406;
+
+		// Just fun
+		// types_ = 10;
+		// seed = 1676103531;
+
+		// Another cell
+		// types_ = 10;
+		// seed = 2396049785;
+
+		log::info() << "Seed: " << seed;
+
+		random::generator rng{seed, 0};
+
+		random::uniform_distribution<int> dcolor(63, 255);
+		random::uniform_distribution<int> dtype(0, types_ - 1);
+		random::uniform_box_point_distribution<float, 2> darea(area_);
+
+		for (int i = 0; i < types_; ++i)
+			colors_.push_back({dcolor(rng), dcolor(rng), dcolor(rng), 255});
+
+		for (int i = 0; i < 4 * 1024; ++i)
+		{
+			particles_.push_back({
+				.position = darea(rng),
+				.velocity = {0.f, 0.f},
+				.mass = 1.f,
+				.radius = 1.f,
+				.type = dtype(rng),
+				// .type = random::uniform_from(rng, {4, 6}),
+			});
+		}
+
+		auto dforce = [](auto & rng)
+		{
+			return random::uniform(rng, 20.f, 100.f) * (random::uniform<bool>(rng) ? 1.f : -1.f);
+		};
+
+		random::uniform_distribution<float> ddistance(3.f, 20.f);
+
+		force_constants_.resize({types_, types_});
+		force_distance_.resize({types_, types_});
+		collision_distance_.resize({types_, types_});
+
+		for (int i = 0; i < types_; ++i)
+			for (int j = 0; j < types_; ++j)
+				force_constants_(i, j) = dforce(rng);
+
+		for (int i = 0; i < types_; ++i)
+			for (int j = 0; j < types_; ++j)
+				force_distance_(i, j) = ddistance(rng);
+
+		for (int i = 0; i < types_; ++i)
+			for (int j = 0; j < types_; ++j)
+				collision_distance_(i, j) = random::uniform(rng, 0.125f, 0.75f) * force_distance_(i, j);
+
+		// for (int i = 0; i < types_; ++i)
+		// 	for (int j = 0; j < types_; ++j)
+		// 		force_constants_(i, j) = -100.f * ((i == j) ? 1.f : 2.f);
+
+		max_force_distance_.resize(types_, 2.f);
+		for (int i = 0; i < types_; ++i)
+			for (int j = 0; j < types_; ++j)
+				math::make_max(max_force_distance_[i], force_distance_(i, j));
+	}
+
+	void on_event(app::mouse_wheel_event const & event) override
+	{
+		app::application_base::on_event(event);
+
+		scale_target_ *= std::pow(1.25f, -event.delta);
+	}
+
+	void on_event(app::key_event const & event) override
+	{
+		app::application_base::on_event(event);
+
+		if (event.down && event.key == app::keycode::SPACE)
+			paused_ ^= true;
+	}
+
+	void update() override
+	{
+		float const dt = 0.02f;
+
+		scale_ += (scale_target_ - scale_) * (- std::expm1(- 20.f * dt));
+
+		auto visible_area = area_;
+		visible_area = math::expand(visible_area, 0.5f * (scale_ - 1.f) * visible_area.dimensions());
+
+		float window_aspect_ratio = state().size[0] * 1.f / state().size[1];
+		float area_aspect_ratio = visible_area[0].length() / visible_area[1].length();
+
+		if (area_aspect_ratio < window_aspect_ratio)
+		{
+			view_area_[1] = visible_area[1];
+			view_area_[0] = math::expand(math::interval<float>::singleton(visible_area[0].center()), visible_area[1].length() * window_aspect_ratio * 0.5f);
+		}
+		else
+		{
+			view_area_[0] = visible_area[0];
+			view_area_[1] = math::expand(math::interval<float>::singleton(visible_area[1].center()), visible_area[0].length() / window_aspect_ratio * 0.5f);
+		}
+
+		if (!paused_)
+		for (int step = 0; step < 5; ++step)
+		{
+			cg::kdtree<float, 2, int> kdtree;
+
+			for (int i = 0; i < particles_.size(); ++i)
+				kdtree.insert({particles_[i].position, i});
+
+			float const collision_strength = 200.f;
+
+			mouseover_particle_ = std::nullopt;
+			if (state().mouse_button_down.contains(app::mouse_button::left))
+			{
+				math::point<float, 2> m;
+				m[0] = math::lerp(view_area_[0], state().mouse[0] * 1.f / state().size[0]);
+				m[1] = math::lerp(view_area_[1], 1.f - state().mouse[1] * 1.f / state().size[1]);
+				mouseover_particle_ = kdtree.closest(m).data;
+			}
+
+			#pragma omp parallel for
+			for (int i = 0; i < particles_.size(); ++i)
+			{
+				// ======== kd-tree based algorithm ========
+
+				auto & pi = particles_[i];
+				auto f = math::vector{0.f, 0.f};
+
+				kdtree.closer_than_map(pi.position, max_force_distance_[pi.type], [&](auto const & value){
+					auto j = value.data;
+
+					if (i == j)
+						return;
+
+					auto const & pj = particles_[j];
+
+					auto r = pj.position - pi.position;
+					float l = math::length(r);
+					auto n = r / l;
+
+					// ==== 1/R forces ====
+					// float lg = std::max(0.5f, l);
+					// auto g = n / (lg * lg);
+
+					// if (l < force_distance_(pi.type, pj.type))
+					// 	f += force_constants_(pi.type, pj.type) * g;
+
+					// float collision_distance = pi.radius + pj.radius;
+
+					// if (l < collision_distance)
+					// {
+					// 	float d = std::max(0.25f, l / collision_distance);
+
+					// 	f -= collision_strength * n * (1.f / (d * d) - 1.f);
+					// }
+
+					// ==== Linear forces ====
+					if (l < force_distance_(pi.type, pj.type))
+						f += force_constants_(pi.type, pj.type) * n * (1.f - l / force_distance_(pi.type, pj.type));
+
+					if (l < collision_distance_(pi.type, pj.type))
+						f -= 10.f * std::abs(force_constants_(pi.type, pj.type)) * n * (1.f - l / collision_distance_(pi.type, pj.type));
+
+					(void)collision_strength;
+				});
+
+				pi.velocity += f * dt / pi.mass;
+
+
+
+
+
+
+				// ======== Naive algorithm v2 ========
+
+				// auto & pi = particles_[i];
+				// auto f = math::vector{0.f, 0.f};
+
+				// for (int j = 0; j < particles_.size(); ++j)
+				// {
+				// 	if (i == j)
+				// 		continue;
+
+				// 	auto const & pj = particles_[j];
+
+				// 	auto r = pj.position - pi.position;
+				// 	float l = math::length(r);
+				// 	auto n = r / l;
+				// 	auto g = n / (l * l);
+
+				// 	if (l < force_distance_(pi.type, pj.type))
+				// 		f += force_constants_(pi.type, pj.type) * g;
+
+				// 	if (l < pi.radius + pj.radius)
+				// 	{
+				// 		float d = l / (pi.radius + pj.radius);
+
+				// 		f -= collision_strength * n * (1.f / (d * d) - 1.f);
+				// 	}
+				// }
+
+				// pi.velocity += f * dt / pi.mass;
+
+
+
+
+
+				// ======== Naive algorithm ========
+
+				// for (int j = i + 1; j < particles_.size(); ++j)
+				// {
+				// 	auto & pi = particles_[i];
+				// 	auto & pj = particles_[j];
+
+				// 	auto fi = math::vector{0.f, 0.f};
+				// 	auto fj = math::vector{0.f, 0.f};
+
+				// 	auto r = pj.position - pi.position;
+				// 	float l = math::length(r);
+				// 	auto n = r / l;
+				// 	auto g = n / (l * l);
+
+				// 	if (l < force_distance_(pi.type, pj.type))
+				// 		fi += force_constants_(pi.type, pj.type) * g;
+
+				// 	if (l < force_distance_(pj.type, pi.type))
+				// 		fj -= force_constants_(pj.type, pi.type) * g;
+
+				// 	if (l < pi.radius + pj.radius)
+				// 	{
+				// 		float d = l / (pi.radius + pj.radius);
+
+				// 		auto f = collision_strength * n * (1.f / (d * d) - 1.f);
+				// 		fi -= f;
+				// 		fj += f;
+				// 	}
+
+				// 	pi.velocity += fi * dt / pi.mass;
+				// 	pj.velocity += fj * dt / pj.mass;
+				// }
+			}
+
+			float const friction = 10.f;
+			float const friction_factor = std::exp(- friction * dt);
+			float const wall_distance = 0.f;
+			float const wall_force = 10.f;
+			float const gravity = 0.f;
+			float const center_gravity = 0.f;
+			float const line_gravity = 0.f;
+			float const circle_gravity = 0.f;
+			float const circle_radius = 100.f;
+			bool const circle_wall = false;
+			float const circle_wall_radius = 100.f;
+
+			for (auto & p : particles_)
+			{
+				auto r = p.position - p.position.zero();
+
+				p.velocity[1] -= gravity * dt;
+				p.velocity[1] -= line_gravity * p.position[1] * dt;
+				p.velocity -= center_gravity * r * dt;
+				p.velocity += circle_gravity * r / math::length(r) * (circle_radius - math::length(r)) * dt;
+				p.velocity *= friction_factor;
+				p.position += p.velocity * dt;
+
+				for (int d : {0, 1})
+				{
+					// if (p.position[d] < area_[d].min + p.radius)
+					// {
+					// 	p.position[d] = area_[d].min + p.radius;
+					// 	p.velocity[d] *= -1.f;
+					// }
+
+					// if (p.position[d] > area_[d].max - p.radius)
+					// {
+					// 	p.position[d] = area_[d].max - p.radius;
+					// 	p.velocity[d] *= -1.f;
+					// }
+
+					if (circle_wall)
+					{
+						auto r = p.position - p.position.zero();
+						auto l = math::length(r);
+						if (l > circle_wall_radius)
+						{
+							p.velocity -= wall_force * r / l * (l - circle_wall_radius) * dt;
+						}
+					}
+					else
+					{
+						if (p.position[d] < area_[d].min + wall_distance)
+						{
+							p.velocity[d] += wall_force * (area_[d].min + wall_distance - p.position[d]) * dt;
+						}
+
+						if (p.position[d] > area_[d].max - wall_distance)
+						{
+							p.velocity[d] -= wall_force * (p.position[d] - (area_[d].max - wall_distance)) * dt;
+						}
+					}
+				}
+			}
+		}
+	}
+
+	void present() override
+	{
+		gl::Viewport(0, 0, state().size[0], state().size[1]);
+
+		gl::ClearColor(0.02f, 0.02f, 0.02f, 1.f);
+		gl::Clear(gl::COLOR_BUFFER_BIT);
+
+		for (auto const & p : particles_)
+		{
+			auto c0 = colors_[p.type];
+			c0[3] = 15;
+
+			auto c1 = colors_[p.type];
+			c1[3] = 0;
+
+			painter_.circle(p.position, p.radius * 8.f, c0, c1);
+		}
+
+		for (auto const & p : particles_)
+			painter_.circle(p.position, p.radius, colors_[p.type]);
+
+		painter_.render(math::orthographic_camera{view_area_}.transform());
+
+		if (mouseover_particle_)
+		{
+			int type = particles_[*mouseover_particle_].type;
+			painter_.text({state().size[0] / 2.f, state().size[1] - 20.f}, std::format("Species #{}", type), {
+				.scale = {2.f, 2.f},
+				.y = gfx::painter::y_align::bottom,
+				.c = colors_[type],
+			});
+		}
+
+		painter_.render(math::window_camera{state().size[0], state().size[1]}.transform());
+	}
+
+private:
+	gfx::painter painter_;
+
+	util::clock<std::chrono::duration<float>, std::chrono::high_resolution_clock> clock_;
+
+	math::box<float, 2> area_;
+	math::box<float, 2> view_area_;
+	float scale_ = 1.f;
+	float scale_target_ = 1.f;
+
+	int types_;
+	std::vector<gfx::color_rgba> colors_;
+	util::ndarray<float, 2> force_constants_;
+	util::ndarray<float, 2> force_distance_;
+	util::ndarray<float, 2> collision_distance_;
+	std::vector<float> max_force_distance_;
+
+	std::vector<particle> particles_;
+
+	std::optional<int> mouseover_particle_;
+
+	bool paused_ = true;
+};
+
+namespace psemek::app
+{
+
+	std::unique_ptr<application::factory> make_application_factory()
+	{
+		return default_application_factory<particle_life_2d_app>({.name = "Particle Life 2D"});
+	}
+
+}

examples/physics.cpp

@@ -507,7 +507,7 @@ void physics_demo_app::update()
 
 		std::vector<util::fixed_vector<std::array<int, 2>, 8>> ball_cells(model.points.size());
 
-		util::array<util::fixed_vector<std::uint32_t, 16>, 2> cells({std::ceil(view_region[0].length() / cell_size), std::ceil(view_region[1].length() / cell_size)});
+		util::ndarray<util::fixed_vector<std::uint32_t, 16>, 2> cells({std::ceil(view_region[0].length() / cell_size), std::ceil(view_region[1].length() / cell_size)});
 
 		float dx = view_region[0].length() / cells.width();
 		float dy = view_region[1].length() / cells.height();

examples/soft_creatures_2d.cpp

@@ -751,7 +751,9 @@ private:
 
 			std::vector<creature> evaluated_creatures(population_.size());
 
+#ifdef _OPENMP
 			#pragma omp parallel for
+#endif
 			for (int i = 0; i < population_.size(); ++i)
 			{
 				auto & creature = population_[i];

examples/srtm.cpp

@@ -648,7 +648,7 @@ srtm_app::srtm_app(options const &, context const & context)
 	selected_mesh.setup<math::point<float, 3>>();
 
 	{
-		util::array<gfx::color_rgb, 1> colors({16});
+		util::ndarray<gfx::color_rgb, 1> colors({16});
 
 		auto * c = colors.data();
 
@@ -674,7 +674,7 @@ srtm_app::srtm_app(options const &, context const & context)
 	}
 
 	{
-		util::array<gfx::color_rgb, 1> colors({5});
+		util::ndarray<gfx::color_rgb, 1> colors({5});
 
 		colors(0) = {0, 63, 127};
 		colors(1) = {0, 0, 127};

examples/water_2d.cpp

@@ -2,8 +2,9 @@
 #include <psemek/app/default_application_factory.hpp>
 #include <psemek/gfx/painter.hpp>
 #include <psemek/math/orthographic.hpp>
+#include <psemek/math/camera.hpp>
 #include <psemek/log/log.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/random/device.hpp>
 #include <psemek/random/generator.hpp>
 #include <psemek/random/uniform.hpp>
@@ -31,31 +32,40 @@ struct water_2d_app
 
 	random::generator rng{random::device{}};
 
-	int const N = 128;
+	int const N = 256;
 
 	float time = 0.f;
 
-	util::array<float, 2> bed;
-	util::array<float, 2> water;
-	util::array<float, 2> flowx;
-	util::array<float, 2> flowy;
-	util::array<math::vector<float, 2>, 2> velocity;
-	util::array<float, 2> sediment;
-	util::array<float, 2> new_sediment;
+	util::ndarray<float, 2> bed;
+	util::ndarray<float, 2> water;
+	util::ndarray<float, 2> flowx;
+	util::ndarray<float, 2> flowy;
+	util::ndarray<math::vector<float, 2>, 2> velocity;
+	util::ndarray<float, 2> sediment;
+	util::ndarray<float, 2> new_sediment;
 
 	gfx::painter painter;
 
 	bool paused = true;
 	bool show_water = true;
 	bool show_velocity = false;
+	bool show_particles = false;
 	bool erosion_on = false;
 	bool rain_on = false;
+
+	struct particle
+	{
+		math::point<float, 2> position;
+		int lifetime;
+	};
+
+	std::vector<particle> particles;
 };
 
 water_2d_app::water_2d_app(options const &, context const & ctx)
 {
 	(void)ctx;
-	ctx.vsync(true);
+	// ctx.vsync(true);
 
 	simulation_area[0] = {-1.f, 1.f};
 	simulation_area[1] = {-1.f, 1.f};
@@ -73,7 +83,7 @@ water_2d_app::water_2d_app(options const &, context const & ctx)
 	// water(N / 2, N / 2) = 100.f;
 
 	random::uniform_sphere_vector_distribution<float, 2> d;
-	util::array<math::vector<float, 2>, 2> grads({9, 9});
+	util::ndarray<math::vector<float, 2>, 2> grads({9, 9});
 	for (auto & v : grads)
 		v = d(rng);
 
@@ -86,8 +96,8 @@ water_2d_app::water_2d_app(options const &, context const & ctx)
 			// if (x == N / 2 && (y != N / 4 && y != 3 * N / 4))
 			// 	bed(x, y) = 100.f;
 
-			float tx = (x + 0.5f) / N;
-			float ty = (y + 0.5f) / N;
+			[[maybe_unused]] float tx = (x + 0.5f) / N;
+			[[maybe_unused]] float ty = (y + 0.5f) / N;
 
 			float n = 0.f;
 
@@ -106,7 +116,9 @@ water_2d_app::water_2d_app(options const &, context const & ctx)
 			// n = 4.f * std::abs(2.f * noise(tx, ty) - 1.f) * std::pow(4.f * tx * (1.f - tx), 2.f);
 
 			// River
-			n = std::min(4.f, math::lerp(10.f, 20.f, noise(tx, ty)) * std::abs(ty - math::lerp({0.4f, 0.6f}, noise(tx, 0.f))));
+			// n = std::min(4.f, math::lerp(10.f, 20.f, noise(tx, ty)) * std::abs(ty - math::lerp({0.4f, 0.6f}, noise(tx, 0.f))));
+
+			// n = 0.f;
 
 			// Canyon
 			// n = 10.f * math::clamp(10.f * math::lerp(-1.f, 1.f, noise(tx, ty)), {0.f, 1.f});
@@ -128,13 +140,17 @@ water_2d_app::water_2d_app(options const &, context const & ctx)
 
 			// if (x == 0)
 			// 	water(x, y) = 100.f - bed(x, y);
+
+			water(x, y) = std::max(0.f, 1.f - bed(x, y)) * (0.5f + 0.5f * noise(tx, ty));
+
+			flowx(x, y) = (2.f * ty - 1.f) * water(x, y) * 10.f / N;
 		}
 	}
 
 	// water(N / 2, N / 2) = 10000.f;
 }
 
-float const dt = 0.002f;
+float const dt = 0.001f;
 
 void water_2d_app::update()
 {
@@ -144,11 +160,15 @@ void water_2d_app::update()
 	float const dx = simulation_area[0].length() / N;
 	float const dy = simulation_area[1].length() / N;
 	float const g = 10.f;
-	float const friction = std::pow(0.5f, dt);
+	float const friction = std::pow(1.f, dt);
 	float const viscosity = 0.f;
 	float const sediment_capacity = 0.1f;
 	float const erosion = 1.f;
 	float const deposition = 10.f;
+	float const coriolis = 0.01f;
+	int const max_particles = 16 * 1024;
+	int const particles_per_frame = 16;
+	int const max_particle_lifetime = max_particles / particles_per_frame;
 
 	// Init boundary flows
 	for (int x = 0; x < N; ++x)
@@ -160,13 +180,16 @@ void water_2d_app::update()
 		// if (x < N / 4) flowy(x, N) = -3.f / N;
 
 		// flowy(x, 0) = 5.f * std::sin(10.f * time) / N;
+
+		flowy(x, 0) = 0.f;
+		flowy(x, N) = 0.f;
 	}
 
 	for (int y = 0; y < N; ++y)
 	{
 		// flowx(0, y) =  5.f * std::sin(10.f * time) / N;
 
-		// flowx(0, y) =  1.f / N;
+		// flowx(0, y) =  10.f / N;
 		// flowx(N, y) =  1.f / N;
 
 		// flowx(0, y) = 10.f * ((y + 0.5f) / N) / N;
@@ -175,12 +198,12 @@ void water_2d_app::update()
 		// if (y < N / 4) flowx(0, y) = 3.f / N;
 		// if (y >= 3 * N / 4) flowx(N, y) = -3.f / N;
 
-		if (bed(0, y) < 0.75f)
-			flowx(0, y) = 4.f / N;
+		// if (bed(0, y) < 0.75f)
+		// 	flowx(0, y) = 4.f / N;
 			// flowx(0, y) = 5.f * math::sqr(std::max(0.f, std::sin(15.f * time))) / N;
 
-		if (bed(N - 1, y) < 0.75f)
-			flowx(N, y) = 5.f / N;
+		// if (bed(N - 1, y) < 0.75f)
+		// 	flowx(N, y) = 5.f / N;
 	}
 
 	// Rain :)
@@ -195,9 +218,13 @@ void water_2d_app::update()
 
 	// Update X flows
 	for (int y = 0; y < N; ++y)
+	{
 		for (int x = 1; x < N; ++x)
 			flowx(x, y) = friction * flowx(x, y) + g * dt * (water(x - 1, y) + bed(x - 1, y) - water(x, y) - bed(x, y));
 
+		flowx(0, y) = friction * flowx(0, y) + g * dt * (water(N - 1, y) + bed(N - 1, y) - water(0, y) - bed(0, y));
+	}
+
 	// Update Y flows
 	for (int y = 1; y < N; ++y)
 		for (int x = 0; x < N; ++x)
@@ -236,7 +263,7 @@ void water_2d_app::update()
 		{
 			float outflow = 0.f;
 			outflow += std::max(0.f, - flowx(x    , y    ));
-			outflow += std::max(0.f,   flowx(x + 1, y    ));
+			outflow += std::max(0.f,   flowx((x + 1) % N, y    ));
 			outflow += std::max(0.f, - flowy(x    , y    ));
 			outflow += std::max(0.f,   flowy(x    , y + 1));
 
@@ -247,7 +274,7 @@ void water_2d_app::update()
 				float scale = std::min(1.f, max_outflow / outflow);
 
 				if (flowx(x, y) < 0.f) flowx(x, y) *= scale;
-				if (flowx(x + 1, y) > 0.f) flowx(x + 1, y) *= scale;
+				if (flowx((x + 1) % N, y) > 0.f) flowx((x + 1) % N, y) *= scale;
 
 				if (flowy(x, y) < 0.f) flowy(x, y) *= scale;
 				if (flowy(x, y + 1) > 0.f) flowy(x, y + 1) *= scale;
@@ -262,7 +289,7 @@ void water_2d_app::update()
 		{
 			float w_old = water(x, y);
 
-			water(x, y) += dt / dx / dy * (flowx(x, y) + flowy(x, y) - flowx(x + 1, y) - flowy(x, y + 1));
+			water(x, y) += dt / dx / dy * (flowx(x, y) + flowy(x, y) - flowx((x + 1) % N, y) - flowy(x, y + 1));
 
 			float w_avg = (w_old + water(x, y)) / 2.f;
 
@@ -270,7 +297,7 @@ void water_2d_app::update()
 
 			if (w_avg > 0.f)
 			{
-				v[0] = (flowx(x, y) + flowx(x + 1, y)) / 2.f / dx / w_avg;
+				v[0] = (flowx(x, y) + flowx((x + 1) % N, y)) / 2.f / dx / w_avg;
 				v[1] = (flowy(x, y) + flowy(x, y + 1)) / 2.f / dy / w_avg;
 			}
 
@@ -278,6 +305,20 @@ void water_2d_app::update()
 		}
 	}
 
+	// Add coriolis force
+	for (int y = 0; y < N; ++y)
+	{
+		for (int x = 0; x < N; ++x)
+		{
+			auto n = math::ort(velocity(x, y)) * coriolis * dt / 2.f * std::sin(float(math::pi) / 2.f * ((y + 0.5f) * 2.f / N - 1.f));
+
+			flowx(x, y) += n[0];
+			flowx((x + 1) % N, y) += n[0];
+			flowy(x, y) += n[1];
+			flowy(x, y + 1) += n[1];
+		}
+	}
+
 	if (erosion_on)
 	{
 		// Erosion-deposition
@@ -332,6 +373,36 @@ void water_2d_app::update()
 		std::swap(sediment, new_sediment);
 	}
 
+	// Init particles
+	for (int i = 0; i < particles_per_frame; ++i)
+		if (particles.size() < max_particles)
+		{
+			math::point pos{random::uniform<float>(rng, 0.f, N), random::uniform<float>(rng, 0.f, N)};
+			if (water(std::floor(pos[0]), std::floor(pos[1])) > 1e-3f)
+				particles.push_back({pos, 0});
+		}
+
+	// Update particles
+	for (auto & p : particles)
+	{
+		p.position += velocity(std::min<int>(N - 1, std::floor(p.position[0])), std::min<int>(N - 1, std::floor(p.position[1]))) * (N * dt);
+		p.position[0] = math::fmod(p.position[0], float(N));
+		p.position[1] = math::clamp(p.position[1], {0.f, N});
+		p.lifetime += 1;
+	}
+
+	// Destroy particles
+	for (int i = 0; i < particles.size();)
+	{
+		if (particles[i].lifetime >= max_particle_lifetime)
+		{
+			std::swap(particles[i], particles.back());
+			particles.pop_back();
+		}
+		else
+			++i;
+	}
+
 	time += dt;
 }
 
@@ -371,8 +442,8 @@ void water_2d_app::present()
 			float w = 1.f - std::exp(- 1.f * water(x, y));
 			float s = 1.f - std::exp(- 1.f * sediment(x, y));
 
-			if (water(x, y) > 1e-3f)
-				w = 0.5f + 0.5f * w;
+			// if (water(x, y) > 1e-3f)
+			// 	w = 0.5f + 0.5f * w;
 
 			auto bed_color = gfx::to_coloru8(gfx::color_4f{0.96f, 0.72f, 0.53f, b});
 			// auto bed_color = gfx::to_coloru8(gfx::color_4f{0.3f, 0.5f, 0.1f, b});
@@ -398,20 +469,33 @@ void water_2d_app::present()
 			{
 				auto center = simulation_area.corner((x + 0.5f) * invN, (y + 0.5f) * invN);
 
-				auto v = water(x, y) * velocity(x, y) / 100.f;
+				auto v = water(x, y) * velocity(x, y) / 40.f;
 
-				float max_length = 0.02f;
+				float max_length = 0.05f;
 				auto l = math::length(v);
 
-				auto color = gfx::to_coloru8(gfx::color_4f{1.f, std::exp(- 10.f * l), 0.f, 1.f});
+				// auto color = gfx::to_coloru8(gfx::color_4f{1.f, std::exp(- 10.f * l), 0.f, 1.f});
+				auto color = gfx::color_rgba{255, 255, 255, 255};
 
 				float s = std::min(1.f, l / max_length);
+				// s = 1.f;
 
-				painter.line(center, center + math::normalized(v) * max_length * s, s * 0.004f, s * 0.002f, color, color, false);
+				auto d = math::normalized(v) * max_length * s / 2.f;
+
+				painter.line(center - d, center + d, s * 0.001f, 0.f, color, color, false);
 			}
 		}
 	}
 
+	if (show_particles)
+	{
+		for (auto const & p : particles)
+		{
+			auto pos = math::lerp(simulation_area, math::vector{p.position[0] / N, p.position[1] / N});
+			painter.circle(pos, 0.005f, {255, 255, 255, 127}, 6);
+		}
+	}
+
 	auto mouse = math::cast<float>(state().mouse);
 
 	mouse[0] = math::lerp(view_area[0], mouse[0] / state().size[0]);
@@ -447,6 +531,11 @@ void water_2d_app::present()
 	}
 
 	painter.render(transform);
+
+	{
+		painter.text({20.f, 20.f}, std::format("{} particles", particles.size()), {.scale = {2.f, 2.f}, .x = gfx::painter::x_align::left, .y = gfx::painter::y_align::top,  .c = {0, 0, 0, 255}});
+		painter.render(math::window_camera{state().size[0], state().size[1]}.transform());
+	}
 }
 
 void water_2d_app::on_event(app::resize_event const & event)
@@ -467,6 +556,9 @@ void water_2d_app::on_event(app::key_event const & event)
 	if (event.down && event.key == app::keycode::V)
 		show_velocity ^= true;
 
+	if (event.down && event.key == app::keycode::P)
+		show_particles ^= true;
+
 	if (event.down && event.key == app::keycode::E)
 		erosion_on ^= true;
 

examples/water_2d_hex.cpp

@@ -0,0 +1,349 @@
+#include <psemek/app/application_base.hpp>
+#include <psemek/app/default_application_factory.hpp>
+#include <psemek/gfx/gl.hpp>
+#include <psemek/gfx/painter.hpp>
+#include <psemek/math/box.hpp>
+#include <psemek/math/camera.hpp>
+#include <psemek/math/gauss.hpp>
+#include <psemek/util/ndarray.hpp>
+#include <psemek/random/generator.hpp>
+#include <psemek/random/device.hpp>
+#include <psemek/random/uniform_sphere.hpp>
+#include <psemek/pcg/perlin.hpp>
+#include <psemek/prof/profiler.hpp>
+
+using namespace psemek;
+
+static const math::vector x_axis {1.f, 0.f};
+static const math::vector y_axis {0.5f, std::sqrt(0.75f)};
+static const math::vector z_axis = y_axis - x_axis;
+static const int N = 64;
+
+static const float dt = 0.1f;
+static const float dx = 1.f;
+static float const g = 10.f;
+static float const friction = std::pow(0.875f, dt);
+
+math::point<float, 2> to_world(int x, int y)
+{
+	return math::point{0.f, 0.f} + x_axis * (x - N / 2.f) + y_axis * (y - N / 2.f);
+}
+
+math::point<float, 2> to_grid(math::point<float, 2> const & p)
+{
+	static auto matrix = *math::inverse(math::by_columns(x_axis, y_axis));
+	static auto zero = math::point{0.f, 0.f};
+
+	return zero + matrix * (p - zero) + math::vector{N / 2.f, N / 2.f};
+}
+
+struct water_2d_hex_app
+	: app::application_base
+{
+	water_2d_hex_app(options const &, context const &);
+
+	void update() override;
+	void present() override;
+
+	void on_event(app::resize_event const & event) override;
+	void on_event(app::key_event const & event) override;
+
+	void stop() override;
+
+private:
+	bool vsync_on_ = true;
+	std::function<void(bool)> set_vsync_;
+
+	random::generator rng_{random::device{}};
+
+	float aspect_ratio_ = 1.f;
+	math::vector<int, 2> screen_size_;
+	math::box<float, 2> view_area_;
+	gfx::painter painter_;
+
+	bool paused_ = false;
+	float time_ = 0.f;
+
+	bool show_velocity_ = false;
+
+	util::ndarray<float, 2> bed_;
+	util::ndarray<float, 2> water_;
+	util::ndarray<float, 2> flow_x_; // flow_x(x, y) is (x-1, y) => (x, y)
+	util::ndarray<float, 2> flow_y_; // flow_y(x, y) is (x, y-1) => (x, y)
+	util::ndarray<float, 2> flow_z_; // flow_z(x, y) is (x, y-1) => (x-1, y)
+};
+
+water_2d_hex_app::water_2d_hex_app(options const &, context const & ctx)
+	: set_vsync_(ctx.vsync)
+{
+	set_vsync_(vsync_on_);
+
+	bed_.resize({N + 1, N + 1}, 0.f);
+	water_.resize({N + 1, N + 1}, 0.f);
+	flow_x_.resize({N + 2, N + 1}, 0.f);
+	flow_y_.resize({N + 1, N + 2}, 0.f);
+	flow_z_.resize({N + 2, N + 2}, 0.f);
+
+	random::uniform_sphere_vector_distribution<float, 2> d_grad;
+
+	util::ndarray<math::vector<float, 2>, 2> perlin_grad({17, 17});
+	for (auto & v : perlin_grad)
+		v = d_grad(rng_);
+
+	pcg::perlin<float, 2> noise(std::move(perlin_grad));
+
+	for (int y = 0; y <= N; ++y)
+	{
+		for (int x = 0; x <= N; ++x)
+		{
+			auto q = (to_world(x, y) - math::point{0.f, 0.f}) / (1.f * N);
+			auto p = q + math::vector{0.5f, 0.5f};
+
+			(void)p;
+
+			// Canyon
+			bed_(x, y) = std::max(0.f, 10.f * std::abs(2.f * noise(p) - 1.f) - 1.5f);
+
+			// Islands
+			// bed_(x, y) = 5.f * math::smoothstep(math::clamp(math::unlerp({0.45f, 0.55f}, noise(p) - 1.f * math::length(q)), {0.f, 1.f}));
+
+			water_(x, y) = 0.f;
+		}
+	}
+}
+
+void water_2d_hex_app::update()
+{
+	{
+		float y_extent = (N / 2) * y_axis[1];
+		view_area_[1] = {- y_extent, y_extent};
+		view_area_[0] = {- y_extent * aspect_ratio_, y_extent * aspect_ratio_};
+	}
+
+	if (state().mouse_button_down.contains(app::mouse_button::left))
+	{
+		auto m = math::lerp(view_area_, math::vector{state().mouse[0] * 1.f / screen_size_[0], 1.f - state().mouse[1] * 1.f / screen_size_[1]});
+		auto p = to_grid(m);
+
+		int const R = 4;
+
+		for (int dy = -R; dy <= R; ++dy)
+		{
+			for (int dx = -R; dx <= R; ++dx)
+			{
+				if (dx + dy < -R || dx + dy > R) continue;
+
+				int mx = std::round(p[0]) + dx;
+				int my = std::round(p[1]) + dy;
+
+				if (mx >= 0 && mx <= N && my >= 0 && my <= N)
+				{
+					auto c = to_world(mx, my);
+					water_(mx, my) += 10.f * std::exp(- 0.25f * math::distance_sqr(c, m)) * dt;
+				}
+			}
+		}
+	}
+
+	if (paused_)
+		return;
+
+	prof::profiler prof("update");
+
+	time_ += dt;
+
+	// Init boundary flows
+	for (int i = 0; i <= N / 2; ++i)
+		if (bed_(N / 2 - i, i) < 0.5f)
+			flow_x_(N / 2 - i, i) = 300.f / N;
+			// flow_x_(N / 2 - i, i) = 0*std::pow(std::sin(0.5f * time_) * std::sin(math::pi * (i * 4.f / N)), 5.f) * 300.f / N;
+
+	for (int i = 0; i <= N / 2; ++i)
+		if (bed_(N - i, N / 2 + i) < 0.5f)
+			flow_x_(N - i + 1, N / 2 + i) = 300.f / N;
+
+	// Update X flows
+	for (int y = 0; y <= N; ++y)
+		// for (int x = 1; x <= N; ++x)
+		for (int x = std::max(1, N / 2 + 1 - y); x <= std::min(N, 3 * N / 2 - y); ++x)
+			// if (x + y - 1 >= N / 2 && x + y <= 3 * N / 2)
+				flow_x_(x, y) = friction * flow_x_(x, y) + g * dt * (water_(x - 1, y) + bed_(x - 1, y) - water_(x, y) - bed_(x, y));
+
+	// Update Y flows
+	for (int y = 1; y <= N; ++y)
+		// for (int x = 0; x <= N; ++x)
+		for (int x = std::max(0, N / 2 + 1 - y); x <= std::min(N, 3 * N / 2 - y); ++x)
+			// if (x + y - 1 >= N / 2 && x + y <= 3 * N / 2)
+				flow_y_(x, y) = friction * flow_y_(x, y) + g * dt * (water_(x, y - 1) + bed_(x, y - 1) - water_(x, y) - bed_(x, y));
+
+	// Update Z flows
+	for (int y = 1; y <= N; ++y)
+		// for (int x = 1; x <= N; ++x)
+		for (int x = std::max(1, N / 2 + 1 - y); x <= std::min(N, 3 * N / 2 + 1 - y); ++x)
+			// if (x + y - 1 >= N / 2 && x + y - 1 <= 3 * N / 2)
+				flow_z_(x, y) = friction * flow_z_(x, y) + g * dt * (water_(x, y - 1) + bed_(x, y - 1) - water_(x - 1, y) - bed_(x - 1, y));
+
+	// Scale flows
+	for (int y = 0; y <= N; ++y)
+	{
+		// for (int x = 0; x <= N; ++x)
+		for (int x = std::max(0, N / 2 - y); x <= std::min(N, 3 * N / 2 - y); ++x)
+		{
+			// if (x + y >= N / 2 && x + y <= 3 * N / 2)
+			{
+				float outflow = 0.f;
+
+				float &  fin1 = flow_x_(x    , y    );
+				float &  fin2 = flow_y_(x    , y    );
+				float &  fin3 = flow_z_(x + 1, y    );
+
+				float & fout1 = flow_x_(x + 1, y    );
+				float & fout2 = flow_y_(x    , y + 1);
+				float & fout3 = flow_z_(x    , y + 1);
+
+				outflow += std::max(0.f, -fin1);
+				outflow += std::max(0.f, -fin2);
+				outflow += std::max(0.f, -fin3);
+				outflow += std::max(0.f, fout1);
+				outflow += std::max(0.f, fout2);
+				outflow += std::max(0.f, fout3);
+
+				if (outflow > 0.f)
+				{
+					float max_outflow = water_(x, y) * dx * dx / dt;
+
+					float scale = std::min(1.f, max_outflow / outflow);
+
+					fin1 *= (fin1 < 0.f ? scale : 1.f);
+					fin2 *= (fin2 < 0.f ? scale : 1.f);
+					fin3 *= (fin3 < 0.f ? scale : 1.f);
+
+					fout1 *= (fout1 > 0.f ? scale : 1.f);
+					fout2 *= (fout2 > 0.f ? scale : 1.f);
+					fout3 *= (fout3 > 0.f ? scale : 1.f);
+				}
+			}
+		}
+	}
+
+	// Update water
+	for (int y = 0; y <= N; ++y)
+		// for (int x = 0; x <= N; ++x)
+		for (int x = std::max(0, N / 2 - y); x <= std::min(N, 3 * N / 2 - y); ++x)
+			// if (x + y >= N / 2 && x + y <= 3 * N / 2)
+				water_(x, y) += dt / dx / dx * (flow_x_(x, y) + flow_y_(x, y) + flow_z_(x + 1,y) - flow_x_(x + 1, y) - flow_y_(x, y + 1) - flow_z_(x, y + 1));
+}
+
+void water_2d_hex_app::present()
+{
+	gl::ClearColor(0.f, 0.f, 0.f, 0.f);
+	gl::Clear(gl::COLOR_BUFFER_BIT);
+
+	for (int y = 0; y < N; ++y)
+	{
+		for (int x = 0; x < N; ++x)
+		{
+			auto color = [this](int x, int y)
+			{
+				auto bed = gfx::color_4f{0.9f, 0.7f, 0.5f, - std::expm1(- bed_(x, y))};
+				auto water = gfx::color_4f{0.0f, 0.25f, 1.f, - std::expm1(- water_(x, y))};
+
+				auto color = gfx::color_4f{0.f, 0.f, 0.f, 0.f};
+				color = gfx::blend(color, bed);
+				color = gfx::blend(color, water);
+				return gfx::to_coloru8(color);
+			};
+
+			auto p00 = to_world(x, y);
+			auto p01 = to_world(x + 1, y);
+			auto p10 = to_world(x, y + 1);
+			auto p11 = to_world(x + 1, y + 1);
+
+			auto b00 = color(x, y);
+			auto b01 = color(x + 1, y);
+			auto b10 = color(x, y + 1);
+			auto b11 = color(x + 1, y + 1);
+
+			if (x + y >= N / 2 && x + y < (3 * N) / 2)
+				painter_.triangle(p00, p01, p10, b00, b01, b10);
+
+			if (x + y + 1 >= N / 2 && x + y + 1 < (3 * N) / 2)
+				painter_.triangle(p10, p01, p11, b10, b01, b11);
+		}
+	}
+
+	if (show_velocity_)
+	for (int y = 0; y <= N; ++y)
+	{
+		for (int x = 0; x <= N; ++x)
+		{
+			if (x + y >= N / 2 && x + y <= (3 * N) / 2)
+			{
+				auto p = to_world(x, y);
+
+				auto vx = (flow_x_(x, y) + flow_x_(x + 1, y)) * 0.5f * x_axis;
+				auto vy = (flow_y_(x, y) + flow_y_(x, y + 1)) * 0.5f * y_axis;
+				auto vz = (flow_z_(x + 1, y) + flow_z_(x, y + 1)) * 0.5f * z_axis;
+
+				auto v = (vx + vy + vz);
+
+				auto M = 1.f;
+				auto l = math::length(v);
+				float s = std::min(1.f, l / M);
+				float c = -std::expm1(- 0.1f * l);
+
+				auto color = gfx::to_coloru8(gfx::color_4f{1.f, 1.f - c, 1.f - c, 1.f});
+
+				painter_.line(p, p + v * (s / l), 0.25f * s, 0.f, color, color, true);
+			}
+		}
+	}
+
+	painter_.render(math::orthographic_camera{view_area_}.transform());
+}
+
+void water_2d_hex_app::on_event(app::resize_event const & event)
+{
+	gl::Viewport(0, 0, event.size[0], event.size[1]);
+	screen_size_ = event.size;
+	aspect_ratio_ = (event.size[0] * 1.f) / event.size[1];
+}
+
+void water_2d_hex_app::on_event(app::key_event const & event)
+{
+	if (event.down)
+	{
+		switch (event.key)
+		{
+		case app::keycode::V:
+			vsync_on_ ^= true;
+			set_vsync_(vsync_on_);
+			break;
+		case app::keycode::C:
+			show_velocity_ ^= true;
+			break;
+		case app::keycode::SPACE:
+			paused_ ^= true;
+			break;
+		default:
+			break;
+		}
+	}
+}
+
+
+void water_2d_hex_app::stop()
+{
+	app::application_base::stop();
+	prof::dump();
+}
+
+namespace psemek::app
+{
+
+	std::unique_ptr<application::factory> make_application_factory()
+	{
+		return default_application_factory<water_2d_hex_app>({.name = "Water 2D hex example", .multisampling = 4});
+	}
+
+}

examples/weather.cpp

@@ -0,0 +1,881 @@
+#include <psemek/app/application_base.hpp>
+#include <psemek/app/default_application_factory.hpp>
+#include <psemek/gfx/painter.hpp>
+#include <psemek/gfx/gl.hpp>
+#include <psemek/math/camera.hpp>
+#include <psemek/math/gradient.hpp>
+#include <psemek/random/generator.hpp>
+#include <psemek/random/device.hpp>
+#include <psemek/random/uniform_ball.hpp>
+#include <psemek/pcg/perlin.hpp>
+#include <psemek/pcg/fractal.hpp>
+#include <psemek/util/ndarray.hpp>
+#include <psemek/log/log.hpp>
+#include <psemek/io/file_stream.hpp>
+
+using namespace psemek;
+
+auto make_perlin(random::generator & rng, int min_octave, int max_octave, float power)
+{
+	std::vector<pcg::perlin<float, 2>> octaves;
+	std::vector<float> weights;
+
+	random::uniform_sphere_vector_distribution<float, 2> random_vector{};
+
+	for (int octave = min_octave; octave < max_octave; ++octave)
+	{
+		int size = 1 << octave;
+		util::ndarray<math::vector<float, 2>, 2> gradients({size + 1, size + 1});
+		for (auto & g : gradients)
+			g = random_vector(rng);
+		octaves.emplace_back(std::move(gradients));
+		weights.push_back(std::pow(power, - octave));
+	}
+
+	float weight_sum = 0.f;
+	for (auto w : weights)
+		weight_sum += w;
+	for (auto & w : weights)
+		w /= weight_sum;
+
+	return pcg::fractal<pcg::perlin<float, 2>>(std::move(octaves), std::move(weights));
+}
+
+void make_force_field(random::generator & rng, util::ndarray<math::vector<float, 2>, 2> & result, float scale)
+{
+	auto noise_1 = make_perlin(rng, 4, 6, 2.f);
+	auto noise_2 = make_perlin(rng, 4, 6, 2.f);
+
+	for (int y = 0; y < result.height(); ++y)
+	{
+		for (int x = 0; x < result.width(); ++x)
+		{
+			math::point p{(x + 0.5f) / result.height(), (y + 0.5f) / result.width()};
+			result(x, y) = scale * (math::vector{noise_1(p), noise_2(p)} * 2.f - math::vector{1.f, 1.f});
+		}
+	}
+}
+
+struct weather_app
+	: app::application_base
+{
+	static constexpr int N = 128;
+
+	const bool static_mode = false;
+
+	const float dt = 20.f;
+	const float viscosity = static_mode ? 0.01f : 0.f;
+	const bool temperature_advection = true;
+	const float advection_magnification = 1.f;
+	const float temperature_diffusion = 0.0004f;
+	const float cooling = 0.01f / 300.f;
+	const float cooling_factor = std::exp(- cooling * dt);
+	const float heating = 323.f * (std::exp(cooling * dt) - 1.f) / dt;
+	const float water_heating_factor = 0.9f;
+	const float coriolis = 0.001f;
+	const float coriolis_bands = 6.f;
+	const float band_force = 0.00001f;
+	const float friction = 0.f;
+	const float slope_friction = 1.f;
+	const float slope_force = 0.001f;
+	const float land_force = 0.05f;
+	const float buoyancy_factor = 0.0002f * 0;
+	const float vorticity_confinement = 0.f;
+	const float elevation_temperature_drop = 30.f;
+	const float evaporation = 1.0f;
+	const float max_humidity_factor = 1.f;
+	const float precipitation_factor = 0.0001f;
+	const float force_field_amplitude = 0.00005f;
+	const float random_forces = 0.25f * (static_mode ? 0.f : 1.f);
+	const float force_field_switch_duration = 720.f * 7.5f; // 7.5 days
+	const int force_field_switch_frames = std::round(force_field_switch_duration / dt);
+	// const float friction_factor = 1.f - std::exp(- friction * dt);
+
+	const bool periodic_x = true;
+
+	// random::generator rng{random::device{}};
+	random::generator rng{0, 0};
+
+	gfx::pixmap_rgba biomes_map;
+
+	float expected_temperature_at(int y, bool water) const
+	{
+		// float latitude = (y - N * 0.5f) * 2.f / N;
+		// return std::cos(latitude * float(math::pi));
+
+		return temperature_income_at(y) * (water ? water_heating_factor : 1.f) * dt / (std::exp(cooling * dt) - 1.f);
+	}
+
+	float temperature_income_at(int y) const
+	{
+		// float latitude = (y - N * 0.5f) * 2.f / N;
+		float latitude = y * 1.f / N;
+		return heating * math::lerp(0.75f, 1.f, std::cos(latitude * float(math::pi) / 2.f));
+		// return heating * math::lerp(0.8f, 1.f, 1.f - std::abs(latitude));
+	}
+
+	int wrap(int i) const
+	{
+		return (i + N) % N;
+	}
+
+	weather_app(options const &, context const &)
+	{
+		simulation_box_ = {{{0.f, N}, {0.f, N}}};
+
+		terrain_.resize({N, N}, 0.f);
+		velocity_.resize({N, N});
+		new_velocity_.resize({N, N});
+		pressure_.resize({N, N}, 0.f);
+		temperature_.resize({N, N}, -1.f);
+		new_temperature_.resize({N, N});
+		average_temperature_.resize({N, N}, 0.f);
+		force_field_main_.resize({N, N});
+		force_field_current_.resize({N, N});
+		force_field_next_.resize({N, N});
+		vorticity_.resize({N, N});
+		humidity_.resize({N, N});
+		new_humidity_.resize({N, N});
+		precipitation_.resize({N, N});
+		average_precipitation_.resize({N, N});
+
+
+		auto terrain_noise = make_perlin(rng, 2, 10, 1.6f);
+
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				auto d = math::length(math::vector{x - N / 2.f, y - N / 2.f}) / (N / 2.f);
+				(void)d;
+				float value = terrain_noise((x + 0.5f) / N, (y + 0.5f) / N);
+				value = pow(value, 4.f) - d / 4.f;
+				value = math::lerp(1.f, 16.f, value);
+				terrain_(x, y) = value;
+			}
+		}
+
+		auto heightmap = gfx::read_image<std::uint8_t>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "heightmap_seed_1.png"});
+
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				terrain_(x, y) = ((heightmap(x, y) / 255.f) * 2048.f - 512.f) / 1024.f;
+			}
+		}
+
+		random::uniform_ball_vector_distribution<float, 2> random_velocity{};
+
+		// for (auto & v : velocity_)
+		// {
+		// 	v = random_velocity(rng) * 0.01f;
+		// 	v += std::cos(0.5f * float(math::pi) * latitude * coriolis_bands
+		// }
+
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				float latitude = (N * 0.5f - y) * 2.f / N;
+				velocity_(x, y) = random_velocity(rng) * 0.f + 0.f * math::vector{-std::cos(0.5f * float(math::pi) * latitude * coriolis_bands), 0.f};
+				temperature_(x, y) = expected_temperature_at(y, terrain_(x, y) <= 0.f);
+			}
+		}
+
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				if (terrain_(x, y) > 0.f)
+					continue;
+
+				float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor;
+				humidity_(x, y) = max_humidity + dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f) * (1.f - precipitation_factor * dt) / precipitation_factor / dt;
+				humidity_(x, y) = 0.f;
+			}
+		}
+
+		make_force_field(rng, force_field_main_, 0.5f);
+		make_force_field(rng, force_field_next_, 0.5f);
+
+		biomes_map = gfx::read_image<gfx::color_rgba>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"});
+	}
+
+	void on_event(app::key_event const & event) override
+	{
+		if (event.down && event.key == app::keycode::SPACE)
+			paused_ ^= true;
+
+		if (event.down && event.key == app::keycode::V)
+			show_velocity_ ^= true;
+
+		if (event.down && event.key == app::keycode::T)
+			show_temperature_ ^= true;
+
+		if (event.down && event.key == app::keycode::D)
+			show_temperature_delta_ ^= true;
+
+		if (event.down && event.key == app::keycode::A)
+			show_average_temperature_delta_ ^= true;
+
+		if (event.down && event.key == app::keycode::P)
+			show_pressure_ ^= true;
+
+		if (event.down && event.key == app::keycode::H)
+			show_land_ ^= true;
+
+		if (event.down && event.key == app::keycode::W)
+			show_water_vapor_ ^= true;
+
+		if (event.down && event.key == app::keycode::R)
+			show_precipitation_ ^= true;
+
+		if (event.down && event.key == app::keycode::Q)
+			show_average_precipitation_ ^= true;
+
+		if (event.down && event.key == app::keycode::B)
+			show_biomes_ ^= true;
+	}
+
+	void update() override
+	{
+		if (paused_)
+			return;
+
+		// Update force field
+		if ((frame_ % force_field_switch_frames) == 0)
+		{
+			std::swap(force_field_current_, force_field_next_);
+			make_force_field(rng, force_field_next_, 0.5f);
+		}
+		[[maybe_unused]] float const force_field_t = ((frame_ % force_field_switch_frames) + 0.5f) / force_field_switch_frames;
+
+		int xmin = periodic_x ? 0 : 1;
+		int xmax = periodic_x ? N : N - 1; // exclusive
+
+		// Temperature source
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				bool const is_water = terrain_(x, y) <= 0.f;
+				// temperature_(x, y) = math::lerp(temperature_(x, y), expected_temperature_at(y), 1.f - std::exp(- heating * dt));
+				temperature_(x, y) += dt * temperature_income_at(y) * (is_water ? water_heating_factor : 1.f);
+				temperature_(x, y) *= cooling_factor;
+			}
+		}
+
+		// Evaporation
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				if (terrain_(x, y) <= 0.f)
+					humidity_(x, y) += dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f);
+
+				// float discharge = std::min(humidity_(x, y), precipitation_factor * dt);
+				// float discharge = humidity_(x, y) * precipitation_factor * dt;
+				// float max_humidity = temperature_(x, y) * max_humidity_factor;
+				float max_humidity = temperature_(x, y) * temperature_(x, y) * max_humidity_factor;
+				float discharge = std::max(0.f, humidity_(x, y) - max_humidity) * precipitation_factor * dt;
+				humidity_(x, y) -= discharge;
+				precipitation_(x, y) = discharge / dt;
+			}
+		}
+
+		// Velocity & temperature advection
+		for (int i = 0; i < N; ++i)
+		{
+			new_temperature_(i, 0) = temperature_(i, 0);
+			new_temperature_(i, N - 1) = temperature_(i,  N - 1);
+
+			new_humidity_(i, 0) = humidity_(i, 0);
+			new_humidity_(i, N - 1) = humidity_(i,  N - 1);
+
+			if (!periodic_x)
+			{
+				new_temperature_(0, i) = temperature_(0, i);
+				new_temperature_(N - 1, i) = temperature_(N - 1, i);
+
+				new_humidity_(0, i) = humidity_(0, i);
+				new_humidity_(N - 1, i) = humidity_(N - 1, i);
+			}
+		}
+
+		for (int y = 1; y < N - 1; ++y)
+		{
+			for (int x = xmin; x < xmax; ++x)
+			{
+				auto v = velocity_(x, y);
+				auto p = math::point{x + 0.5f, y + 0.5f} - (advection_magnification * dt) * v;
+				p[0] = p[0] - 0.5f;
+				p[1] = math::clamp(p[1] - 0.5f, {0.f, N - 1.f});
+
+				if (!periodic_x)
+					p[0] = math::clamp(p[0], {0.f, N - 1.f});
+
+				int ix = std::floor(p[0]);
+				int iy = std::min<int>(N - 1, std::floor(p[1]));
+
+				if (!periodic_x)
+					ix = std::min(N - 1, ix);
+
+				float tx = p[0] - ix;
+				float ty = p[1] - iy;
+
+				new_velocity_(x, y) = math::lerp(
+					math::lerp(velocity_(wrap(ix + 0), iy + 0), velocity_(wrap(ix + 1), iy + 0), tx),
+					math::lerp(velocity_(wrap(ix + 0), iy + 1), velocity_(wrap(ix + 1), iy + 1), tx),
+					ty
+				);
+
+				new_temperature_(x, y) = math::lerp(
+					math::lerp(temperature_(wrap(ix + 0), iy + 0), temperature_(wrap(ix + 1), iy + 0), tx),
+					math::lerp(temperature_(wrap(ix + 0), iy + 1), temperature_(wrap(ix + 1), iy + 1), tx),
+					ty
+				);
+
+				new_humidity_(x, y) = math::lerp(
+					math::lerp(humidity_(wrap(ix + 0), iy + 0), humidity_(wrap(ix + 1), iy + 0), tx),
+					math::lerp(humidity_(wrap(ix + 0), iy + 1), humidity_(wrap(ix + 1), iy + 1), tx),
+					ty
+				);
+			}
+		}
+		std::swap(velocity_, new_velocity_);
+		if (temperature_advection) std::swap(temperature_, new_temperature_);
+		std::swap(humidity_, new_humidity_);
+
+		// Apply velocity diffusion
+		for (int y = 0; y < N; ++y)
+			for (int x = 0; x < N; ++x)
+				new_velocity_(x, y) = velocity_(x, y);
+		for (int y = 1; y < N - 1; ++y)
+		{
+			for (int x = xmin; x < xmax; ++x)
+			{
+				// Velocity Laplacian
+				auto laplacian = velocity_(wrap(x + 1), y) + velocity_(wrap(x - 1), y) + velocity_(x, y + 1) + velocity_(x, y - 1) - 4.f * velocity_(x, y);
+				new_velocity_(x, y) = velocity_(x, y) + viscosity * dt * laplacian;
+			}
+		}
+		std::swap(velocity_, new_velocity_);
+
+		// Apply temperature diffusion
+		for (int y = 0; y < N; ++y)
+			for (int x = 0; x < N; ++x)
+				new_temperature_(x, y) = temperature_(x, y);
+		for (int y = 1; y < N - 1; ++y)
+		{
+			for (int x = xmin; x < xmax; ++x)
+			{
+				// Temperature Laplacian
+				auto laplacian = temperature_(wrap(x + 1), y) + temperature_(wrap(x - 1), y) + temperature_(x, y + 1) + temperature_(x, y - 1) - 4.f * temperature_(x, y);
+				new_temperature_(x, y) = temperature_(x, y) + temperature_diffusion * dt * laplacian;
+			}
+		}
+		std::swap(temperature_, new_temperature_);
+
+		// Compute vorticity
+		for (int y = 1; y < N - 1; ++y)
+			for (int x = xmin; x < xmax; ++x)
+				vorticity_(x, y) = (velocity_(x, y + 1)[0] - velocity_(x, y - 1)[0]) / 2.f - (velocity_(wrap(x + 1), y)[1] - velocity_(wrap(x - 1), y)[1]) / 2;
+
+		// Apply forces & friction
+		for (int y = 1; y < N - 1; ++y)
+		{
+			for (int x = xmin; x < xmax; ++x)
+			{
+				// float latitude = (N * 0.5f - y) * 2.f / N;
+				[[maybe_unused]] float latitude = (N - y) * 1.f / N;
+				// velocity_(x, y) += math::ort(velocity_(x, y)) * (coriolis * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands));
+				velocity_(x, y) = math::rotate(velocity_(x, y), coriolis * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands));
+
+				// auto force = force_field_main_(x, y) + random_forces * math::lerp(force_field_current_(x, y), force_field_next_(x, y), force_field_t);
+				// velocity_(x, y) += (dt * force_field_amplitude) * force;
+
+				// velocity_(x, y)[0] += dt * 0.000001f * std::cos(0.5f * float(math::pi) * latitude * 4.f);
+				// velocity_(x, y)[0] += dt * 0.000001f;
+
+				velocity_(x, y)[0] += dt * band_force * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands);
+
+				[[maybe_unused]] math::vector terrain_gradient
+				{
+					(std::max(0.f, terrain_(x + 1, y)) - std::max(0.f, terrain_(x - 1, y))) / 2.f,
+					(std::max(0.f, terrain_(x, y + 1)) - std::max(0.f, terrain_(x, y - 1))) / 2.f,
+				};
+
+				[[maybe_unused]] math::vector temperature_gradient
+				{
+					(temperature_(x + 1, y) - temperature_(x - 1, y)) / 2.f,
+					(temperature_(x, y + 1) - temperature_(x, y - 1)) / 2.f,
+				};
+
+				velocity_(x, y) += temperature_gradient * buoyancy_factor * dt;
+
+				[[maybe_unused]] float slope_factor = std::exp(- dt * slope_friction * math::dot(math::normalized(velocity_(x, y)), terrain_gradient));
+				velocity_(x, y) *= std::min(1.f, slope_factor);
+				// velocity_(x, y) -= terrain_gradient * slope_force * dt;
+
+				// [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * std::pow(math::length(terrain_gradient), 4.f));
+				[[maybe_unused]] float land_factor = std::exp(- dt * land_force * std::pow(std::max(0.f, terrain_(x, y)), 1.f));
+				velocity_(x, y) *= land_factor;
+
+				// Directional external force
+				// velocity_(x, y)[1] += 0.001f * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands);
+				// velocity_(x, y) += math::direction(frame_ * dt * 2.f * float(math::pi) / 10080.f) * 0.001f * dt;
+
+				math::vector vorticity_gradient
+				{
+					(vorticity_(wrap(x + 1), y) - vorticity_(wrap(x - 1), y)) / 2.f,
+					(vorticity_(x, y + 1) - vorticity_(x, y - 1)) / 2.f,
+				};
+				if (auto l = math::length(vorticity_gradient); l > 0.f)
+					vorticity_gradient /= l;
+
+				velocity_(x, y) += vorticity_confinement * dt * vorticity_(x, y) * math::ort(vorticity_gradient);
+
+				float local_friction = friction * terrain_(x, y);
+				// velocity_(x, y) -= local_friction * velocity_(x, y) * math::length(velocity_(x, y));
+				[[maybe_unused]] float local_friction_factor = std::exp(- local_friction * dt);
+				// velocity_(x, y) *= local_friction_factor;
+			}
+		}
+
+		// Solve Poisson equation for pressure
+		for (int iteration = 0; iteration < 16; ++iteration)
+		{
+			int ymin = ((iteration % 2) == 0) ? 1 : N - 2;
+			int ymax = ((iteration % 2) == 0) ? N - 1 : 0;
+			int ystep = ((iteration % 2) == 0) ? 1 : -1;
+
+			for (int y = ymin; y != ymax; y += ystep)
+			{
+				for (int x = xmin; x < xmax; ++x)
+				{
+					// Velocity divergence
+					float divergence = (velocity_(wrap(x + 1), y)[0] - velocity_(wrap(x - 1), y)[0] + velocity_(x, y + 1)[1] - velocity_(x, y - 1)[1]) / 2.f;
+
+					// Gauss-Seidel iteration step
+					pressure_(x, y) = (pressure_(wrap(x - 1), y) + pressure_(wrap(x + 1), y) + pressure_(x, y - 1) + pressure_(x, y + 1) - divergence) / 4.f;
+				}
+			}
+		}
+
+		// Apply boundary conditions for pressure
+		for (int i = 0; i < N; ++i)
+		{
+			if (!periodic_x)
+			{
+				pressure_(0, i) = pressure_(1, i);
+				pressure_(N - 1, i) = pressure_(N - 2, i);
+			}
+			pressure_(i, 0) = pressure_(i, 1);
+			pressure_(i, N - 1) = pressure_(i, N - 2);
+		}
+		if (!periodic_x)
+		{
+			pressure_(0, 0) = (pressure_(0, 1) + pressure_(1, 0)) / 2.f;
+			pressure_(N-1, 0) = (pressure_(N-1, 1) + pressure_(N-2, 0)) / 2.f;
+			pressure_(0, N-1) = (pressure_(0, N-2) + pressure_(1, N-2)) / 2.f;
+			pressure_(N-1, N-1) = (pressure_(N-1, N-2) + pressure_(N-2, N-1)) / 2.f;
+		}
+
+		// Normalize pressure
+		float average_pressure = 0.f;
+		for (auto const & value : pressure_)
+			average_pressure += value;
+		average_pressure /= (1.f * N * N);
+		for (auto & value : pressure_)
+			value -= average_pressure;
+
+		// Project velocity into divergence-free space
+		// by subtracting pressure gradient
+		for (int y = 1; y < N - 1; ++y)
+		{
+			for (int x = xmin; x < xmax; ++x)
+			{
+				// Pressure gradient
+				math::vector gradient{
+					(pressure_(wrap(x + 1), y) - pressure_(wrap(x - 1), y)) / 2.f,
+					(pressure_(x, y + 1) - pressure_(x, y - 1)) / 2.f
+				};
+				velocity_(x, y) -= gradient;
+			}
+		}
+
+		// Apply boundary conditions for velocity
+		for (int i = 0; i < N; ++i)
+		{
+			if (!periodic_x)
+			{
+				float left_boundary_flow = 0.f;//0.01f * std::sin((i * 1.f / N) * float(math::pi) * 4.f);
+				float right_boundary_flow = -left_boundary_flow;
+
+				velocity_(1, i)[0] = left_boundary_flow;
+				velocity_(N-2, i)[0] = right_boundary_flow;
+
+				velocity_(0, i)[0] = - velocity_(1, i)[0];
+				velocity_(N-1, i)[0] = - velocity_(N-2, i)[0];
+			}
+			velocity_(i, 0)[1] = -velocity_(i, 1)[1];
+			velocity_(i, N-2)[1] = -velocity_(i, N-2)[1];
+
+			// velocity_(i, 1)[0] = 0.01f;
+			// velocity_(i, N-2)[0] = 0.01f;
+		}
+
+		// Uncomment to visualize the force field
+		// for (int y = 0; y < N; ++y)
+		// 	for (int x = 0; x < N; ++x)
+		// 		velocity_(x, y) = 100000.f * force_field_(x, y);
+
+		// Uncomment to visualize the terrain gradient field
+		// for (int y = 1; y < N - 1; ++y)
+		// {
+		// 	for (int x = 1; x < N - 1; ++x)
+		// 	{
+		// 		math::vector terrain_gradient
+		// 		{
+		// 			(terrain_(x + 1, y) - terrain_(x - 1, y)) / 2.f,
+		// 			(terrain_(x, y + 1) - terrain_(x, y - 1)) / 2.f,
+		// 		};
+
+		// 		velocity_(x, y) = terrain_gradient;
+		// 	}
+		// }
+
+		// Apply boundary conditions for humidity
+		for (int i = 0; i < N; ++i)
+		{
+			if (!periodic_x)
+			{
+				humidity_(0, i) = 0.f;
+				humidity_(N - 1, i) = 0.f;
+			}
+			humidity_(i, 0) = 0.f;
+			humidity_(i, N - 1) = 0.f;
+		}
+
+		++frame_;
+
+		// Update all-time average temperature & precipitation
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				// float t = 1.f / frame_;
+				// float t = 1.f / std::min(8192, frame_);
+				float t = 1.f / std::min<float>(86400.f / dt, frame_); // year average
+
+				if (static_mode)
+					t = 1.f;
+
+				average_temperature_(x, y) = math::lerp(average_temperature_(x, y), temperature_(x, y), t);
+				average_precipitation_(x, y) = math::lerp(average_precipitation_(x, y), precipitation_(x, y), t);
+			}
+		}
+	}
+
+	void present() override
+	{
+		gl::ClearColor(0.f, 0.f, 0.f, 0.f);
+		gl::Clear(gl::COLOR_BUFFER_BIT);
+
+		float const aspect_ratio = state().size[0] * 1.f / state().size[1];
+		math::box<float, 2> view_box = math::expand(simulation_box_, 1.f);
+		if (view_box[0].length() / view_box[1].length() > aspect_ratio)
+			view_box[1] = math::expand(view_box[1], (view_box[0].length() / aspect_ratio - view_box[1].length()) / 2.f);
+		else
+			view_box[0] = math::expand(view_box[0], (view_box[1].length() * aspect_ratio - view_box[0].length()) / 2.f);
+
+		std::optional<math::vector<int, 2>> mouseover_cell;
+		{
+			auto mouse = math::lerp(view_box, math::vector{state().mouse[0] * 1.f / state().size[0], 1.f - state().mouse[1] * 1.f / state().size[1]});
+			int x = std::floor(mouse[0]);
+			int y = std::floor(mouse[1]);
+
+			if (x >= 0 && x < N && y >= 0 && y < N)
+				mouseover_cell = {x, y};
+		}
+
+		if (mouseover_cell && (state().mouse_button_down.contains(app::mouse_button::left) || state().mouse_button_down.contains(app::mouse_button::right)))
+		{
+			float delta = state().mouse_button_down.contains(app::mouse_button::left) ? 1.f : -1.f;
+
+			int R = 4;
+
+			for (int iy = -R; iy <= R; ++iy)
+			{
+				for (int ix = -R; ix <= R; ++ix)
+				{
+					int x = (*mouseover_cell)[0] + ix;
+					int y = (*mouseover_cell)[1] + iy;
+
+					if (x >= 0 && x < N && y >= 0 && y < N)
+					{
+						auto d = math::vector<float, 2>{ix, iy} / 2.f;
+						terrain_(x, y) += delta * 0.05f * std::exp(- math::dot(d, d));
+					}
+				}
+			}
+		}
+
+		if (mouseover_cell && state().mouse_button_down.contains(app::mouse_button::middle))
+		{
+			int R = 4;
+
+			float average = 0.f;
+			int count = 0;
+
+			for (int iy = -R; iy <= R; ++iy)
+			{
+				for (int ix = -R; ix <= R; ++ix)
+				{
+					int x = (*mouseover_cell)[0] + ix;
+					int y = (*mouseover_cell)[1] + iy;
+
+					if (x >= 0 && x < N && y >= 0 && y < N)
+					{
+						average += terrain_(x, y);
+						count += 1;
+					}
+				}
+			}
+
+			average /= count;
+
+			for (int iy = -R; iy <= R; ++iy)
+			{
+				for (int ix = -R; ix <= R; ++ix)
+				{
+					int x = (*mouseover_cell)[0] + ix;
+					int y = (*mouseover_cell)[1] + iy;
+
+					if (x >= 0 && x < N && y >= 0 && y < N)
+					{
+						auto d = math::vector<float, 2>{ix, iy} / 2.f;
+						terrain_(x, y) += (average - terrain_(x, y)) * 0.05f * std::exp(- math::dot(d, d));
+					}
+				}
+			}
+		}
+
+		[[maybe_unused]] float const pixel_size = view_box[0].length() / state().size[0];
+
+		auto map_color = [](float value, gfx::color_4f const & negative, gfx::color_4f const & positive){
+			return math::lerp(negative, positive, 1.f/ (1.f + std::exp(- value)));
+		};
+
+		auto map_temperature = [&](float value) {
+			return map_color(2.f * std::round(value / 20.f), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f});
+		};
+
+		auto map_biome = [this](float temperature, float precipitation)
+		{
+			auto x = math::clamp<int>(math::unlerp({  -3.f,  5.f}, precipitation) * biomes_map.width() , {0, biomes_map.width()  - 1});
+			auto y = math::clamp<int>(math::unlerp({-10.f, 30.f}, temperature  ) * biomes_map.height(), {0, biomes_map.height() - 1});
+
+			return gfx::to_colorf(biomes_map(x, y));
+		};
+
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				gfx::color_4f color = gfx::color_4f::zero();
+
+				if (show_land_ || show_biomes_)
+				{
+					if (!show_biomes_)
+					{
+						if (terrain_(x, y) <= 0.f)
+							color = {0.5f, 0.5f, 1.f, 1.f};
+						else
+							color = {1.f, 1.f, 1.f, 1.f};
+					}
+					else
+					{
+						if (terrain_(x, y) <= 0.f)
+							color = map_color(8.f * terrain_(x, y), {0.f, 0.f, 0.125f, 1.f}, {0.f, 1.f, 1.5f, 1.f});
+						else
+						{
+							float temperature = average_temperature_(x, y) - 273.f - std::max(0.f, terrain_(x, y)) * elevation_temperature_drop;
+							// float precipitation = (std::log10(std::max(1e-9f, average_precipitation_(x, y))) + 3.f) * 2.f;
+							// float precipitation = std::pow(average_precipitation_(x, y) * 125.f, 2.f) * 8.f;
+							float precipitation = std::log2(std::max(1e-9f, average_precipitation_(x, y)));
+							color = map_biome(temperature, precipitation);
+						}
+					}
+
+					if (show_land_ && x > 0 && x + 1 < N && y > 0 && y + 1 < N)
+					{
+						math::vector terrain_gradient
+						{
+							(terrain_(x + 1, y) - terrain_(x - 1, y)) / 2.f,
+							(terrain_(x, y + 1) - terrain_(x, y - 1)) / 2.f,
+						};
+
+						auto terrain_normal = math::normalized(math::vector{-terrain_gradient[0], -terrain_gradient[1], 0.125f});
+
+						float lightness = 0.5f + 0.5f * math::dot(terrain_normal, math::normalized(math::vector{1.f, 2.f, 3.f}));
+
+						color = gfx::dark(color, 1.f - lightness);
+					}
+				}
+
+				bool const is_water = terrain_(x, y) <= 0.f;
+
+				if (show_temperature_)
+					color = map_temperature(temperature_(x, y) - 273.f);
+
+				if (show_temperature_delta_)
+					color = gfx::blend(color, map_color((temperature_(x, y) - expected_temperature_at(y, is_water)), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f}));
+
+				if (show_average_temperature_delta_)
+					color = gfx::blend(color, map_color((average_temperature_(x, y) - expected_temperature_at(y, is_water)), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f}));
+
+				if (show_pressure_)
+					color = gfx::blend(color, map_color(10000.f * pressure_(x, y), {0.f, 0.f, 1.f, 0.75f}, {1.f, 0.f, 0.f, 0.75f}));
+
+				if (show_water_vapor_)
+					color = gfx::blend(color, map_color(humidity_(x, y) * 0.00001f, {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));
+
+				if (show_precipitation_)
+				{
+					// color = gfx::blend(color, map_color(precipitation_(x, y), {1.f, 1.f, 1.f, -1.f}, {1.f, 1.f, 1.f, 1.f}));
+
+					float alpha = 2.f / (1.f + std::exp(- 0.1f * precipitation_(x, y))) - 1.f;
+
+					gfx::color_4f cloud_color{1.f, 1.f, 1.f, alpha * 0.875f};
+
+					if (y > 0 && y + 1 < N)
+					{
+						if (periodic_x || (x > 0 && x + 1 < N))
+						{
+							math::vector gradient
+							{
+								(precipitation_(wrap(x + 1), y) - precipitation_(wrap(x - 1), y)) / 2.f,
+								(precipitation_(x, y + 1) - precipitation_(x, y - 1)) / 2.f,
+							};
+
+							auto normal = math::normalized(math::vector{-gradient[0], -gradient[1], 2.f});
+							auto lightness = 0.5f + 0.5f * math::dot(normal, math::normalized(math::vector{1.f, 2.f, 3.f}));
+							cloud_color = gfx::dark(cloud_color, 1.f - lightness);
+						}
+					}
+
+					color = gfx::blend(color, cloud_color);
+				}
+
+				if (show_average_precipitation_)
+					color = gfx::blend(color, map_color(average_precipitation_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));
+
+				painter_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color));
+			}
+		}
+
+		if (show_velocity_)
+		{
+			for (int y = 0; y < N; ++y)
+			{
+				for (int x = 0; x < N; ++x)
+				{
+					math::point center{x + 0.5f, y + 0.5f};
+					auto v = velocity_(x, y);
+					auto color = gfx::color_4f::zero();
+					if (auto l = math::length(v); l > 0.f)
+					{
+						float const magnification = 200.f;
+						float const max_length = 1.5f;
+						v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l;
+
+						// color = gfx::lerp(gfx::color_4f{0.5f, 1.f, 0.f, 1.f}, gfx::color_4f{1.f, 0.f, 0.f, 1.f}, 1.f - std::exp(- 0.25f * magnification * l));
+						color = map_color(0.1f * (temperature_(x, y) - 273.f), {0.125f, 0.5f, 1.f, 1.f}, {1.f, 0.5f, 0.125f, 1.f});
+					}
+					auto n = math::ort(v) * 0.3f;
+
+					painter_.triangle(center - v - n, center - v + n, center + v, gfx::to_coloru8(color));
+				}
+			}
+		}
+
+		auto push_text = [&, row = 0](std::string const & text) mutable
+		{
+			painter_.text(view_box.corner(0.f, 1.f) - math::vector{0.f, row * pixel_size * 2.f * 12.f}, text, {.scale = {2.f * pixel_size, - 2.f * pixel_size}, .x = gfx::painter::x_align::left, .y = gfx::painter::y_align::top, .c = {255, 255, 255, 255}});
+			++row;
+		};
+
+		push_text(std::format("Frame {}", frame_));
+		push_text(std::format("Day {:.2f}", frame_ / (720.f / dt)));
+
+		if (mouseover_cell)
+		{
+			int x = (*mouseover_cell)[0];
+			int y = (*mouseover_cell)[1];
+			painter_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, {255, 255, 255, 127});
+
+			push_text(std::format("{} {}", x, y));
+			push_text(std::format("V = {:.3f} {:.3f}", velocity_(x, y)[0] * 1000.f, velocity_(x, y)[1] * 1000.f));
+			push_text(std::format("P = {:.3f}", pressure_(x, y) * 1000.f));
+			push_text(std::format("T = {:.3f}", temperature_(x, y) - 273.f));
+			push_text(std::format("A = {:.3f}", average_temperature_(x, y) - 273.f));
+			push_text(std::format("E = {:.3f}", expected_temperature_at(y, terrain_(x, y) <= 0.f) - 273.f));
+			push_text(std::format("H = {:.3f}", terrain_(x, y)));
+			push_text(std::format("W = {:.3f}", humidity_(x, y)));
+			push_text(std::format("R = {:.3f}", precipitation_(x, y)));
+			push_text(std::format("AR= {:.3f}", average_precipitation_(x, y)));
+		}
+
+		painter_.render(math::orthographic_camera{view_box}.transform());
+	}
+
+private:
+	gfx::painter painter_;
+
+	math::box<float, 2> simulation_box_;
+
+	bool paused_ = true;
+	bool show_velocity_ = true;
+	bool show_temperature_ = false;
+	bool show_temperature_delta_ = false;
+	bool show_average_temperature_delta_ = false;
+	bool show_pressure_ = false;
+	bool show_land_ = true;
+	bool show_biomes_ = true;
+	bool show_water_vapor_ = false;
+	bool show_precipitation_ = false;
+	bool show_average_precipitation_ = false;
+
+	util::ndarray<float, 2> terrain_;
+	util::ndarray<math::vector<float, 2>, 2> velocity_;
+	util::ndarray<math::vector<float, 2>, 2> new_velocity_;
+	util::ndarray<float, 2> pressure_;
+	util::ndarray<float, 2> vorticity_;
+
+	util::ndarray<float, 2> temperature_;
+	util::ndarray<float, 2> new_temperature_;
+	util::ndarray<float, 2> average_temperature_;
+
+	util::ndarray<float, 2> humidity_;
+	util::ndarray<float, 2> new_humidity_;
+	util::ndarray<float, 2> precipitation_;
+	util::ndarray<float, 2> average_precipitation_;
+
+	util::ndarray<math::vector<float, 2>, 2> force_field_main_;
+	util::ndarray<math::vector<float, 2>, 2> force_field_current_;
+	util::ndarray<math::vector<float, 2>, 2> force_field_next_;
+
+	int frame_ = 0;
+};
+
+namespace psemek::app
+{
+
+	std::unique_ptr<application::factory> make_application_factory()
+	{
+		return default_application_factory<weather_app>({.name = "Weather simulation test"});
+	}
+
+}

libs/app/include/psemek/app/application.hpp

@@ -3,6 +3,7 @@
 #include <psemek/app/event_handler.hpp>
 
 #if defined(PSEMEK_GRAPHICS_API_WEBGPU)
+#include <psemek/wgpu/instance.hpp>
 #include <psemek/wgpu/adapter.hpp>
 #include <psemek/wgpu/surface.hpp>
 #include <psemek/wgpu/device.hpp>
@@ -24,9 +25,9 @@ namespace psemek::app
 			int multisampling = 4;
 			bool highdpi = false;
 #if defined(PSEMEK_GRAPHICS_API_WEBGPU)
-			std::vector<wgpu::feature> required_features;
-			std::optional<wgpu::limits> required_limits;
-			std::optional<wgpu::native_limits> required_native_limits;
+			std::vector<wgpu::feature> required_features = {};
+			std::optional<wgpu::limits> required_limits = {};
+			std::optional<wgpu::native_limits> required_native_limits = {};
 #endif
 		};
 
@@ -37,10 +38,12 @@ namespace psemek::app
 			std::function<void(bool)> show_cursor;
 			std::function<void(bool)> relative_mouse_mode;
 			std::function<void(bool)> windowed;
+			std::function<void(bool)> text_input;
 #if defined(PSEMEK_GRAPHICS_API_OPENGL)
 			std::function<void(bool)> vsync;
 #endif
 #if defined(PSEMEK_GRAPHICS_API_WEBGPU)
+			wgpu::instance instance;
 			wgpu::adapter adapter;
 			wgpu::surface surface;
 			wgpu::device device;

libs/app/include/psemek/app/application_base.hpp

@@ -18,6 +18,7 @@ namespace psemek::app
 		void on_event(touch_up_event const &) override;
 		void on_event(touch_move_event const &) override;
 		void on_event(key_event const &) override;
+		void on_event(text_input_event const &) override;
 
 		void stop() override;
 		bool running() const override;

libs/app/include/psemek/app/event_handler.hpp

@@ -16,6 +16,7 @@ namespace psemek::app
 		virtual void on_event(touch_up_event const &) {}
 		virtual void on_event(touch_move_event const &) {}
 		virtual void on_event(key_event const &) {}
+		virtual void on_event(text_input_event const &) {}
 
 		virtual ~event_handler() {}
 	};

libs/app/include/psemek/app/event_state.hpp

@@ -61,4 +61,7 @@ namespace psemek::app
 			state.key_down.erase(event.key);
 	}
 
+	inline void apply(event_state &, text_input_event const &)
+	{}
+
 }

libs/app/include/psemek/app/events.hpp

@@ -2,6 +2,8 @@
 
 #include <psemek/math/point.hpp>
 
+#include <string>
+
 namespace psemek::app
 {
 
@@ -172,4 +174,9 @@ namespace psemek::app
 		bool down;
 	};
 
+	struct text_input_event
+	{
+		std::string input;
+	};
+
 }

libs/app/include/psemek/app/scene_application.hpp

@@ -18,6 +18,7 @@ namespace psemek::app
 		void on_event(touch_up_event const &) override;
 		void on_event(touch_move_event const &) override;
 		void on_event(key_event const &) override;
+		void on_event(text_input_event const &) override;
 
 		void stop() override;
 

libs/app/source/application_base.cpp

@@ -48,6 +48,9 @@ namespace psemek::app
 		apply(state_, event);
 	}
 
+	void application_base::on_event(text_input_event const &)
+	{}
+
 	void application_base::stop()
 	{
 		running_ = false;

libs/app/source/scene_application.cpp

@@ -48,6 +48,11 @@ namespace psemek::app
 		on_event_impl(event);
 	}
 
+	void scene_application::on_event(text_input_event const & event)
+	{
+		on_event_impl(event);
+	}
+
 	void scene_application::stop()
 	{
 		if (auto scene = current_scene())

libs/async/CMakeLists.txt

@@ -3,4 +3,8 @@ file(GLOB_RECURSE PSEMEK_ASYNC_SOURCES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "s
 
 psemek_add_library(psemek-async ${PSEMEK_ASYNC_HEADERS} ${PSEMEK_ASYNC_SOURCES})
 target_include_directories(psemek-async PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
-target_link_libraries(psemek-async PUBLIC psemek-log psemek-util)
+target_link_libraries(psemek-async PUBLIC
+	psemek-log
+	psemek-util
+	psemek-prof
+)

libs/async/include/psemek/async/event_loop.hpp

@@ -17,6 +17,8 @@ namespace psemek::async
 
 		void stop() override;
 
+		void clear() override;
+
 		void wait() override;
 
 		void wait_for(clock::duration period) override;

libs/async/include/psemek/async/executor.hpp

@@ -4,7 +4,6 @@
 #include <psemek/util/function.hpp>
 
 #include <chrono>
-#include <memory>
 
 namespace psemek::async
 {
@@ -30,6 +29,9 @@ namespace psemek::async
 		// NB: the executor must call stop() from destructor.
 		virtual void stop() = 0;
 
+		// Clear the pending tasks queue.
+		virtual void clear() = 0;
+
 		// Wait for all the tasks to be executed.
 		// May take forever.
 		virtual void wait() = 0;

libs/async/include/psemek/async/future.hpp

@@ -40,6 +40,18 @@ namespace psemek::async
 			using type = util::function<void()>;
 		};
 
+		template <typename T>
+		struct future_result_type
+		{
+			using type = T &;
+		};
+
+		template <>
+		struct future_result_type<void>
+		{
+			using type = void;
+		};
+
 		struct cancel_token
 		{};
 
@@ -55,7 +67,8 @@ namespace psemek::async
 			std::weak_ptr<cancel_token> weak_cancel;
 
 			std::mutex then_mutex;
-			typename then_function<T>::type then_func;
+			std::vector<typename then_function<T>::type> then_funcs;
+			std::vector<util::function<void(std::exception_ptr)>> then_exception_funcs;
 		};
 
 	}
@@ -85,14 +98,16 @@ namespace psemek::async
 		using result_type = T;
 
 		future() = default;
-		future(future&&) = default;
+		future(future const &) = default;
+		future(future &&) = default;
 
 		future(std::shared_ptr<detail::task_state<T>> state, std::shared_ptr<detail::cancel_token> cancel)
 			: state_(std::move(state))
 			, cancel_(std::move(cancel))
 		{}
 
-		future & operator = (future&&) = default;
+		future & operator = (future const &) = default;
+		future & operator = (future &&) = default;
 
 		~future()
 		{
@@ -136,7 +151,7 @@ namespace psemek::async
 			return has_value_unsafe();
 		}
 
-		T get()
+		typename detail::future_result_type<T>::type get()
 		{
 			if (!state_)
 				throw empty_future_error{};
@@ -148,7 +163,7 @@ namespace psemek::async
 				if constexpr (std::is_same_v<T, void>)
 					return;
 				else
-					return std::move(*(state_->value));
+					return *(state_->value);
 			}
 			else
 				std::rethrow_exception(state_->exception);
@@ -239,8 +254,10 @@ namespace psemek::async
 				state_->value = value;
 			}
 			state_->value_cv.notify_all();
-			if (state_->then_func)
-				state_->then_func(*state_->value);
+
+			std::lock_guard lock{state_->then_mutex};
+			for (auto & func : state_->then_funcs)
+				func(*state_->value);
 		}
 
 		void set_value(T && value)
@@ -252,8 +269,10 @@ namespace psemek::async
 				state_->value = std::move(value);
 			}
 			state_->value_cv.notify_all();
-			if (state_->then_func)
-				state_->then_func(*state_->value);
+
+			std::lock_guard lock{state_->then_mutex};
+			for (auto & func : state_->then_funcs)
+				func(*state_->value);
 		}
 
 		void set_exception(std::exception_ptr e)
@@ -263,6 +282,10 @@ namespace psemek::async
 				std::lock_guard lock{state_->value_mutex};
 				if (state_->value || state_->exception) throw satisfied_promise_error{};
 				state_->exception = std::move(e);
+
+				std::lock_guard lock_then{state_->then_mutex};
+				for (auto & func : state_->then_exception_funcs)
+					func(state_->exception);
 			}
 			state_->value_cv.notify_all();
 		}
@@ -324,8 +347,10 @@ namespace psemek::async
 				state_->value = true;
 			}
 			state_->value_cv.notify_all();
-			if (state_->then_func)
-				state_->then_func();
+
+			std::lock_guard lock{state_->then_mutex};
+			for (auto & func : state_->then_funcs)
+				func();
 		}
 
 		void set_exception(std::exception_ptr e)
@@ -335,6 +360,10 @@ namespace psemek::async
 				std::lock_guard lock{state_->value_mutex};
 				if (state_->value || state_->exception) throw satisfied_promise_error{};
 				state_->exception = std::move(e);
+
+				std::lock_guard lock_then{state_->then_mutex};
+				for (auto & func : state_->then_exception_funcs)
+					func(state_->exception);
 			}
 			state_->value_cv.notify_all();
 		}
@@ -362,6 +391,14 @@ namespace psemek::async
 		return p.get_future();
 	}
 
+	template <typename T>
+	future<T> make_ready_future(T const & x)
+	{
+		promise<T> p;
+		p.set_value(x);
+		return p.get_future();
+	}
+
 	template <typename Signature>
 	struct packaged_task;
 
@@ -488,7 +525,7 @@ namespace psemek::async
 			auto fut = t.get_future();
 
 			std::lock_guard lock{state_->then_mutex};
-			state_->then_func = std::move(t);
+			state_->then_funcs.push_back(std::move(t));
 			return fut;
 		}
 		else
@@ -498,7 +535,7 @@ namespace psemek::async
 			auto fut = t.get_future();
 
 			std::lock_guard lock{state_->then_mutex};
-			state_->then_func = std::move(t);
+			state_->then_funcs.push_back(std::move(t));
 			return fut;
 		}
 	}

libs/async/include/psemek/async/synchronous_executor.hpp

@@ -14,6 +14,8 @@ namespace psemek::async
 
 		void stop() override;
 
+		void clear() override;
+
 		void wait() override;
 
 		void wait_for(clock::duration period) override;

libs/async/include/psemek/async/threadpool.hpp

@@ -4,7 +4,6 @@
 #include <psemek/util/thread.hpp>
 #include <psemek/util/synchronyzed_queue.hpp>
 
-#include <future>
 #include <vector>
 #include <string>
 #include <mutex>
@@ -32,6 +31,8 @@ namespace psemek::async
 
 		void stop() override;
 
+		void clear() override;
+
 		void wait() override;
 
 		void wait_for(clock::duration period) override;

libs/async/source/event_loop.cpp

@@ -1,7 +1,5 @@
 #include <psemek/async/event_loop.hpp>
 
-#include <algorithm>
-
 namespace psemek::async
 {
 
@@ -20,6 +18,11 @@ namespace psemek::async
 		queue_.clear();
 	}
 
+	void event_loop::clear()
+	{
+		queue_.clear();
+	}
+
 	void event_loop::wait()
 	{
 		while (!queue_.empty())

libs/async/source/synchronous_executor.cpp

@@ -1,8 +1,6 @@
 #include <psemek/async/synchronous_executor.hpp>
 #include <psemek/util/exception.hpp>
 
-#include <stdexcept>
-
 namespace psemek::async
 {
 
@@ -19,6 +17,9 @@ namespace psemek::async
 	void synchronous_executor::stop()
 	{}
 
+	void synchronous_executor::clear()
+	{}
+
 	void synchronous_executor::wait()
 	{}
 

libs/async/source/threadpool.cpp

@@ -3,19 +3,13 @@
 #include <psemek/util/unused.hpp>
 #include <psemek/util/to_string.hpp>
 #include <psemek/util/exception.hpp>
+#include <psemek/prof/profiler.hpp>
 
 #include <psemek/log/log.hpp>
 
 namespace psemek::async
 {
 
-	namespace
-	{
-
-		struct stop_execution{};
-
-	}
-
 	threadpool::threadpool(std::string const & name, std::size_t thread_count)
 		: working_count_{0}
 	{
@@ -25,10 +19,13 @@ namespace psemek::async
 			threads_.emplace_back([this, tname = std::move(tname)]() mutable
 			{
 				log::thread_registrator reg(std::move(tname));
-				for (bool running = true; running;)
+				while (true)
 				{
 					auto task = task_queue_.pop();
 
+					if (!task)
+						break;
+
 					{
 						std::lock_guard lock{working_count_mutex_};
 						++working_count_;
@@ -36,12 +33,9 @@ namespace psemek::async
 
 					try
 					{
+						prof::profiler prof("task");
 						task();
 					}
-					catch (stop_execution const &)
-					{
-						running = false;
-					}
 					catch (util::exception const & e)
 					{
 						log::error() << "Unhandled exception in threadpool executor: " << e;
@@ -86,11 +80,16 @@ namespace psemek::async
 		for (auto const & thread: threads_)
 		{
 			unused(thread);
-			task_queue_.push([]{ throw stop_execution{}; });
+			task_queue_.push(nullptr);
 		}
 		threads_.clear();
 	}
 
+	void threadpool::clear()
+	{
+		task_queue_.clear();
+	}
+
 	void threadpool::wait()
 	{
 		std::unique_lock lock{working_count_mutex_};

libs/audio/include/psemek/audio/channel.hpp

@@ -18,12 +18,12 @@ namespace psemek::audio
 
 		stream_ptr stream() const
 		{
-			return std::atomic_load(&stream_);
+			return stream_.load();
 		}
 
 		stream_ptr stream(stream_ptr new_stream)
 		{
-			return std::atomic_exchange(&stream_, std::move(new_stream));
+			return stream_.exchange(std::move(new_stream));
 		}
 
 		stream_ptr stop()
@@ -37,7 +37,7 @@ namespace psemek::audio
 		}
 
 	private:
-		stream_ptr stream_;
+		std::atomic<stream_ptr> stream_;
 	};
 
 	using channel_ptr = std::shared_ptr<channel>;

libs/bt/CMakeLists.txt

@@ -4,6 +4,3 @@ file(GLOB_RECURSE PSEMEK_BT_SOURCES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "sour
 psemek_add_library(psemek-bt ${PSEMEK_BT_HEADERS} ${PSEMEK_BT_SOURCES})
 target_include_directories(psemek-bt PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
 target_link_libraries(psemek-bt PUBLIC psemek-util psemek-log)
-if(PSEMEK_BT_LOG)
-	target_compile_definitions(psemek-bt PUBLIC -DPSEMEK_BT_LOG=1)
-endif()

libs/bt/include/psemek/bt/conditional.hpp

@@ -39,10 +39,7 @@ namespace psemek::bt
 			{
 				if (condition_result_)
 				{
-					if (*condition_result_)
-						return true_branch_->update(dt, args...);
-					else
-						return false_branch_->update(dt, args...);
+					return current_node()->update(dt, args...);
 				}
 				else
 				{
@@ -51,11 +48,8 @@ namespace psemek::bt
 					if (auto finished = std::get_if<typename node_type::finished>(&result))
 					{
 						condition_result_ = finished->result;
-						if (*condition_result_)
-							true_branch_->start(args...);
-						else
-							false_branch_->start(args...);
-						return running{};
+						current_node()->start(args...);
+						return current_node()->update(Time{}, args...);
 					}
 
 					return result;
@@ -71,10 +65,7 @@ namespace psemek::bt
 		{
 			if (condition_result_)
 			{
-				if (*condition_result_)
-					return true_branch_->event(event, args...);
-				else
-					return false_branch_->event(event, args...);
+				return current_node()->event(event, args...);
 			}
 			else
 				return condition_->event(event, args...);
@@ -85,6 +76,11 @@ namespace psemek::bt
 		node_ptr<tree_type> true_branch_;
 		node_ptr<tree_type> false_branch_;
 		std::optional<bool> condition_result_;
+
+		node<tree_type> * current_node()
+		{
+			return *condition_result_ ? true_branch_.get() : false_branch_.get();
+		}
 	};
 
 	template <typename Tree>

libs/bt/include/psemek/bt/log.hpp

@@ -1,9 +0,0 @@
-#pragma once
-
-#ifdef PSEMEK_BT_LOG
-#include <psemek/log/log.hpp>
-#define bt_log ::psemek::log::debug
-#else
-#include <psemek/util/null_ostream.hpp>
-#define bt_log ::psemek::util::null_ostream
-#endif

libs/bt/include/psemek/bt/selector.hpp

@@ -32,7 +32,7 @@ namespace psemek::bt
 
 		status update(Time dt, Args ... args) override
 		{
-			if (current_index_ < children_.size())
+			while (current_index_ < children_.size())
 			{
 				if (!current_started_)
 				{
@@ -41,13 +41,14 @@ namespace psemek::bt
 				}
 
 				auto result = children_[current_index_]->update(dt, args...);
+				dt = Time{};
 				if (auto f = std::get_if<finished>(&result))
 				{
 					if (f->result)
 						return finished{true};
 					++current_index_;
 					current_started_ = false;
-					return running{};
+					continue;
 				}
 				else
 					return result;

libs/bt/include/psemek/bt/sequence.hpp

@@ -32,7 +32,7 @@ namespace psemek::bt
 
 		status update(Time dt, Args ... args) override
 		{
-			if (current_index_ < children_.size())
+			while (current_index_ < children_.size())
 			{
 				if (!current_started_)
 				{
@@ -41,13 +41,14 @@ namespace psemek::bt
 				}
 
 				auto result = children_[current_index_]->update(dt, args...);
+				dt = Time{};
 				if (auto f = std::get_if<finished>(&result))
 				{
 					if (!f->result)
 						return finished{false};
 					++current_index_;
 					current_started_ = false;
-					return running{};
+					continue;
 				}
 				else
 					return result;

libs/cg/include/psemek/cg/kdtree.hpp

@@ -77,6 +77,9 @@ namespace psemek::cg
 		template <typename Iterator>
 		kdtree(Iterator begin, Iterator end);
 
+		template <typename Iterator>
+		void assign(Iterator begin, Iterator end);
+
 		bool empty() const { return nodes_.empty(); }
 
 		bool insert(value_type && value);
@@ -84,11 +87,19 @@ namespace psemek::cg
 		// TODO: implement
 		bool remove(point_type const & point) const;
 
+		void clear();
+
 		// TODO: implement
 		// TODO: alternative non-const version that allows modifying value.data
 		value_type const * find(point_type const & point) const;
 
-		value_type const & closest(point_type const & target);
+		value_type const & closest(point_type const & target) const;
+
+		template <typename Iterator>
+		Iterator closer_than(point_type const & target, scalar_type max_distance, Iterator out) const;
+
+		template <typename Callback>
+		void closer_than_map(point_type const & target, scalar_type max_distance, Callback && callback) const;
 
 	private:
 
@@ -106,6 +117,7 @@ namespace psemek::cg
 		{
 			value_type value;
 			node_id children[2] {null, null};
+			std::uint32_t size = 0;
 		};
 
 		std::vector<node> nodes_;
@@ -116,6 +128,12 @@ namespace psemek::cg
 		bool insert_impl(value_type && value, node_id id, std::uint32_t split_axis);
 
 		value_type const * closest_impl(point_type const & target, scalar_type & best_distance_sqr, node_id id, std::uint32_t split_axis) const;
+
+		template <typename Iterator>
+		Iterator closer_than_impl(point_type const & target, scalar_type max_distance_sqr, Iterator out, node_id id, std::uint32_t split_axis) const;
+
+		template <typename Callback>
+		void closer_than_map_impl(point_type const & target, scalar_type max_distance_sqr, Callback & callback, node_id id, std::uint32_t split_axis) const;
 	};
 
 	template <typename T, std::size_t N, typename Data>
@@ -125,6 +143,14 @@ namespace psemek::cg
 		build_node_impl(begin, end, 0);
 	}
 
+	template <typename T, std::size_t N, typename Data>
+	template <typename Iterator>
+	void kdtree<T, N, Data>::assign(Iterator begin, Iterator end)
+	{
+		clear();
+		build_node_impl(begin, end, 0);
+	}
+
 	template <typename T, std::size_t N, typename Data>
 	bool kdtree<T, N, Data>::insert(value_type && value)
 	{
@@ -132,7 +158,13 @@ namespace psemek::cg
 	}
 
 	template <typename T, std::size_t N, typename Data>
-	kdtree<T, N, Data>::value_type const & kdtree<T, N, Data>::closest(point_type const & target)
+	void kdtree<T, N, Data>::clear()
+	{
+		nodes_.clear();
+	}
+
+	template <typename T, std::size_t N, typename Data>
+	kdtree<T, N, Data>::value_type const & kdtree<T, N, Data>::closest(point_type const & target) const
 	{
 		if (nodes_.empty())
 			throw util::exception("empty kdtree");
@@ -141,6 +173,26 @@ namespace psemek::cg
 		return *closest_impl(target, best_distance_sqr, 0, 0);
 	}
 
+	template <typename T, std::size_t N, typename Data>
+	template <typename Iterator>
+	Iterator kdtree<T, N, Data>::closer_than(point_type const & target, scalar_type max_distance, Iterator out) const
+	{
+		if (nodes_.empty())
+			return out;
+
+		return closer_than_impl(target, math::sqr(max_distance), out, 0, 0);
+	}
+
+	template <typename T, std::size_t N, typename Data>
+	template <typename Callback>
+	void kdtree<T, N, Data>::closer_than_map(point_type const & target, scalar_type max_distance, Callback && callback) const
+	{
+		if (nodes_.empty())
+			return;
+
+		return closer_than_map_impl(target, math::sqr(max_distance), callback, 0, 0);
+	}
+
 	template <typename T, std::size_t N, typename Data>
 	template <typename Iterator>
 	kdtree<T, N, Data>::node_id kdtree<T, N, Data>::build_node_impl(Iterator begin, Iterator end, std::uint32_t split_axis)
@@ -160,6 +212,7 @@ namespace psemek::cg
 		auto result = static_cast<node_id>(nodes_.size());
 		auto & node = nodes_.emplace_back();
 		node.value = std::move(*middle);
+		node.size = (end - begin);
 
 		nodes_[result].children[0] = build_node_impl(begin, middle, next_axis(split_axis));
 		nodes_[result].children[1] = build_node_impl(std::next(middle), end, next_axis(split_axis));
@@ -174,6 +227,7 @@ namespace psemek::cg
 		{
 			auto & node = nodes_.emplace_back();
 			node.value = std::move(value);
+			node.size = 1;
 			return true;
 		}
 
@@ -184,7 +238,7 @@ namespace psemek::cg
 		if (node_point == point)
 			return false;
 
-		int child = (node.point[split_axis] < point[split_axis]) ? 0 : 1;
+		int child = (point[split_axis] < node_point[split_axis]) ? 0 : 1;
 
 		if (node.children[child] == null)
 		{
@@ -193,10 +247,18 @@ namespace psemek::cg
 
 			auto & child_node = nodes_.emplace_back();
 			child_node.value = std::move(value);
+			child_node.size = 1;
+			nodes_[id].size += 1;
 			return true;
 		}
 
-		return insert_impl(std::move(value), node.children[child], next_axis(split_axis));
+		if (insert_impl(std::move(value), node.children[child], next_axis(split_axis)))
+		{
+			nodes_[id].size += 1;
+			return true;
+		}
+
+		return false;
 	}
 
 	template <typename T, std::size_t N, typename Data>
@@ -224,4 +286,46 @@ namespace psemek::cg
 		return result;
 	}
 
+	template <typename T, std::size_t N, typename Data>
+	template <typename Iterator>
+	Iterator kdtree<T, N, Data>::closer_than_impl(point_type const & target, scalar_type max_distance_sqr, Iterator out, node_id id, std::uint32_t split_axis) const
+	{
+		auto const & node = nodes_[id];
+		auto const & node_point = detail::get_point(node.value);
+
+		if (math::distance_sqr(node_point, target) < max_distance_sqr)
+			*out++ = node.value;
+
+		auto delta = target[split_axis] - node_point[split_axis];
+		auto delta_sqr = math::sqr(delta);
+
+		if (node.children[0] != null && (delta < 0 || delta_sqr < max_distance_sqr))
+			out = closer_than_impl(target, max_distance_sqr, out, node.children[0], next_axis(split_axis));
+
+		if (node.children[1] != null && (delta >= 0 || delta_sqr < max_distance_sqr))
+			out = closer_than_impl(target, max_distance_sqr, out, node.children[1], next_axis(split_axis));
+
+		return out;
+	}
+
+	template <typename T, std::size_t N, typename Data>
+	template <typename Callback>
+	void kdtree<T, N, Data>::closer_than_map_impl(point_type const & target, scalar_type max_distance_sqr, Callback & callback, node_id id, std::uint32_t split_axis) const
+	{
+		auto const & node = nodes_[id];
+		auto const & node_point = detail::get_point(node.value);
+
+		if (math::distance_sqr(node_point, target) < max_distance_sqr)
+			callback(node.value);
+
+		auto delta = target[split_axis] - node_point[split_axis];
+		auto delta_sqr = math::sqr(delta);
+
+		if (node.children[0] != null && (delta < 0 || delta_sqr < max_distance_sqr))
+			closer_than_map_impl(target, max_distance_sqr, callback, node.children[0], next_axis(split_axis));
+
+		if (node.children[1] != null && (delta >= 0 || delta_sqr < max_distance_sqr))
+			closer_than_map_impl(target, max_distance_sqr, callback, node.children[1], next_axis(split_axis));
+	}
+
 }

libs/ecs/include/psemek/ecs/container.hpp

@@ -22,11 +22,8 @@ namespace psemek::ecs
 
 	// TODO:
 	// - Fully document which functions can be called from which callbacks
-	// - Modification callbacks implementation (const-only)
 	// - Tables serialization
 	// - Refactor query caches
-	// - Index API
-	// - Index implementation
 
 	struct container
 	{
@@ -82,6 +79,9 @@ namespace psemek::ecs
 		 */
 		void destroy(handle entity);
 
+		/** Same as destroy(), but via finally() */
+		void destroy_finally(handle entity);
+
 		/** Get an accessor for an entity, which provides access to the entity's components.
 		 *  It is designed to be a single-use object; it cannot be stored as a reference to
 		 *    the entity. Use handle for that instead.
@@ -172,6 +172,10 @@ namespace psemek::ecs
 		template <typename ... Components>
 		void detach(handle entity);
 
+		/** Same as detach(), but via finally() */
+		template <typename ... Components>
+		void detach_finally(handle entity);
+
 		/** Create a query cache that can be used to speed up `apply()` calls.
 		 *
 		 *  The constness of the component types is ignored.
@@ -296,12 +300,7 @@ namespace psemek::ecs
 		 *  The constructor function must have the same signature as a function
 		 *    passed to the `apply<Components...>()` call.
 		 *
-		 *  The constuctor is not considered to be a modification of the entity, i.e. it
-		 *    doesn't trigger modification callbacks.
-		 *
 		 *  @param function A function to be applied to created entity's components
-		 *  @return An ownerwhip token; destroying this token removes
-		 *    the constructor from this container
 		 *  @warning If any two of the passed component types are equal, the call fails with
 		 *    a compilation error
 		 *  @warning If the constructor modifies the entity's archetype (i.e. attaches or
@@ -335,11 +334,7 @@ namespace psemek::ecs
 		 *
 		 *  Note that there is no way to cancel the entity's destruction.
 		 *
-		 *  The destructor is not considered to be a modification of the entity, i.e. it
-		 *    doesn't trigger modification callbacks.
-		 *
 		 *  @param function A function to be applied to a to-be destroyed entity's components
-		 *  @return An ownerwhip token; destroying it removes the destructor from the container
 		 *  @warning If any two of the passed component types are equal, the call fails with
 		 *    a compilation error
 		 *  @warning If the destructor modifies the entity's archetype (i.e. attaches or
@@ -349,45 +344,36 @@ namespace psemek::ecs
 		template <typename ... Components, typename Function>
 		void destructor(Function && function);
 
-		/** Register a component modification callback. Each time an entity that has
-		 *    the specified set of components is modified, and the modification affects
-		 *    one of this callback's component types, the callback is called.
+		/** Call a callback after exiting all currently executing ECS container methods.
 		 *
-		 *  The component types can be const-qualified, in which case the corresponding function must
-		 *    also accept the corresponding components by a const reference.
+		 *  If no ECS container method is currently executed, call the callback immediately instead.
 		 *
-		 *  The component types can be equal to ecs::without<Component>, indicating that entities having
-		 *    this component type will not be watched by this callback. These component types are not
-		 *    included in the called function signature.
+		 *  The function must have one of the following signatures:
+		 *      void()
+		 *      void(container)
 		 *
-		 *  If the modification occurred via an accessor, the callback is called
-		 *    after the accessor is destroyed, allowing for transaction-like modification.
+		 *  This method is meant to be called from inside other callbacks (typically
+		 *    constructors and destructors) and can be used to add extra modifications
+		 *    to entities that don't clash with other callbacks and archetype changes.
 		 *
-		 *  If the modification occurred via an `apply()` or `batch_apply()` call,
-		 *    the callback is called immediately after this call.
+		 *  Generally, if a constructor, destructor, apply or batch apply function
+		 *    changes the archetype of some entities, it should do so using `finally`
+		 *    instead of doing that directly.
 		 *
-		 *  Note that modifications from within a modification callback are not considered
-		 *    as modifications, i.e. they don't recursively invoke modification callbacks.
-		 *
-		 *  The callback function must have the same signature as a function
-		 *    passed to the `apply<Components...>()` call.
-		 *
-		 *  @param function A function to be applied to modified entity's components
-		 *  @return An ownerwhip token; destroying it removes the modification callback from the container
-		 *  @warning If any two of the passed component types are equal, the call fails with
-		 *    a compilation error
-		 *  @warning If the modification callback modifies the entity's archetype (i.e. attaches or
-		 *    detaches components), the behavior is undefined
-		 *
-		 *  TODO: can we allow the callback to modify the archetype?
+		 *  @param function A callback to be called before the topmost currently executing
+		 *    ECS container method returns.
+		 *  @warning If a callback adds new callbacks via `finally`, those will also get executed,
+		 *    which potentially leads to an infinite loop.
 		 */
-		// TODO: implement
-		template <typename ... Components, typename Function>
-		void watch(Function && function);
+		template <typename Function>
+		void finally(Function && function);
 
 		template <typename Index, typename ... Args>
 		Index & index(Args && ... args);
 
+		template <typename Index, typename Factory>
+		Index & index_factory(Factory && factory);
+
 		template <typename ... Components>
 		std::size_t memory_usage();
 
@@ -412,10 +398,18 @@ namespace psemek::ecs
 		util::object_pool<std::vector<util::uuid>> uuid_list_pool_;
 		util::object_pool<util::hash_set<util::uuid>> uuid_set_pool_;
 
+		std::size_t method_recursion_depth_ = 0;
+		std::vector<util::function<void(container &)>> finally_callbacks_;
+
+		#ifdef PSEMEK_DEBUG
+			util::hash_set<ecs::handle> currently_changing_archetype_;
+		#endif
+
 		detail::table * insert_table(std::vector<std::unique_ptr<detail::column>> columns);
 		void do_destroy(handle entity);
 		void remove_row(detail::table & table, std::uint32_t row, util::span<detail::entity_data> entities);
 		void finalize_iteration(detail::table & table);
+		void finalize_method();
 	};
 
 	template <typename Component>
@@ -429,6 +423,8 @@ namespace psemek::ecs
 	{
 		static_assert(detail::all_different_types_v<std::remove_cvref_t<Components>...>, "all component types must be different");
 
+		++method_recursion_depth_;
+
 		(register_component<std::remove_cvref_t<Components>>(), ...);
 
 		detail::component_uuid_helper<std::remove_cvref_t<Components>...> uuids;
@@ -449,6 +445,9 @@ namespace psemek::ecs
 		auto row = table->row_count();
 		auto id = entity_list_.create(table, row);
 		handle handle{id, entity_list_.get_entities()[id].epoch};
+		#ifdef PSEMEK_DEBUG
+			currently_changing_archetype_.insert(handle);
+		#endif
 		[[maybe_unused]] accessor accessor = get(handle);
 
 		table->push_row(handle);
@@ -456,6 +455,14 @@ namespace psemek::ecs
 
 		table->trigger_constructors(*this, row);
 
+		#ifdef PSEMEK_DEBUG
+			currently_changing_archetype_.erase(handle);
+		#endif
+
+		finalize_method();
+
+		--method_recursion_depth_;
+
 		return handle;
 	}
 
@@ -464,6 +471,13 @@ namespace psemek::ecs
 	{
 		static_assert(detail::all_different_types_v<std::remove_cvref_t<Components>...>, "all component types must be different");
 
+		#ifdef PSEMEK_DEBUG
+			assert(!currently_changing_archetype_.contains(entity));
+			currently_changing_archetype_.insert(entity);
+		#endif
+
+		++method_recursion_depth_;
+
 		(register_component<Components>(), ...);
 
 		auto uuids = uuid_list_pool_.get();
@@ -479,6 +493,9 @@ namespace psemek::ecs
 		for (auto const & uuid : attached_uuid_set)
 			uuids.push_back(uuid);
 
+		if (archetype_changed)
+			data->table->trigger_destructors(*this, data->row, attached_uuid_set, {});
+
 		auto table = table_container_.get(uuids);
 
 		if (!table)
@@ -508,13 +525,21 @@ namespace psemek::ecs
 		((accessor.get<std::remove_cvref_t<Components>>() = std::forward<Components>(components)), ...);
 
 		if (archetype_changed)
-			table->trigger_constructors(*this, data->row, attached_uuid_set);
+			table->trigger_constructors(*this, data->row, attached_uuid_set, {});
 
 		attached_uuid_set.clear();
 		uuids.clear();
 
 		uuid_set_pool_.put(std::move(attached_uuid_set));
 		uuid_list_pool_.put(std::move(uuids));
+
+		#ifdef PSEMEK_DEBUG
+			currently_changing_archetype_.erase(entity);
+		#endif
+
+		finalize_method();
+
+		--method_recursion_depth_;
 	}
 
 	template <typename ... Components>
@@ -522,6 +547,13 @@ namespace psemek::ecs
 	{
 		static_assert(detail::all_different_types_v<std::remove_const_t<Components>...>, "all component types must be different");
 
+		#ifdef PSEMEK_DEBUG
+			assert(!currently_changing_archetype_.contains(entity));
+			currently_changing_archetype_.insert(entity);
+		#endif
+
+		++method_recursion_depth_;
+
 		auto detached_uuid_set = uuid_set_pool_.get();
 		(detached_uuid_set.insert(std::remove_const_t<Components>::uuid()), ...);
 
@@ -530,10 +562,9 @@ namespace psemek::ecs
 		auto * data = entity_list_.get_entities().begin() + entity.id;
 		for (auto const & column : data->table->columns())
 		{
-			if (detached_uuid_set.contains(column->uuid()))
-				detached_uuid_set.insert(column->uuid());
-			else
-				uuids.push_back(column->uuid());
+			auto const column_uuid = column->uuid();
+			if (!detached_uuid_set.contains(column_uuid))
+				uuids.push_back(column_uuid);
 		}
 
 		bool const archetype_changed = sizeof...(Components) > 0;
@@ -541,7 +572,7 @@ namespace psemek::ecs
 		auto table = table_container_.get(uuids);
 
 		if (archetype_changed)
-			data->table->trigger_destructors(*this, data->row, detached_uuid_set);
+			data->table->trigger_destructors(*this, data->row, {}, detached_uuid_set);
 
 		// Destructors could lead to reallocation of entity list
 		data = entity_list_.get_entities().begin() + entity.id;
@@ -568,12 +599,32 @@ namespace psemek::ecs
 			data->row = new_row;
 		}
 
+		if (archetype_changed)
+			table->trigger_constructors(*this, data->row, {}, detached_uuid_set);
+
 		detached_uuid_set.clear();
 		uuids.clear();
 		detached_uuid_set.clear();
 
 		uuid_list_pool_.put(std::move(uuids));
 		uuid_set_pool_.put(std::move(detached_uuid_set));
+
+		#ifdef PSEMEK_DEBUG
+			currently_changing_archetype_.erase(entity);
+		#endif
+
+		finalize_method();
+
+		--method_recursion_depth_;
+	}
+
+	template <typename ... Components>
+	void container::detach_finally(handle entity)
+	{
+		finally([entity](container & world){
+			if (world.alive(entity))
+				world.detach<Components...>(entity);
+		});
 	}
 
 	template <typename ... Components>
@@ -590,6 +641,8 @@ namespace psemek::ecs
 	{
 		static_assert(detail::all_different_types_v<std::remove_const_t<Components>...>, "all component types must be different");
 
+		++method_recursion_depth_;
+
 		using invocable_type = typename detail::filter_with<detail::invocable, std::tuple<Components...>, Function>::type;
 
 		static_assert(invocable_type::value, "function is not invocable with these components");
@@ -617,6 +670,10 @@ namespace psemek::ecs
 			finalize_iteration(*entry.table);
 		}
 
+		finalize_method();
+
+		--method_recursion_depth_;
+
 		return cache;
 	}
 
@@ -625,6 +682,8 @@ namespace psemek::ecs
 	{
 		static_assert(detail::all_different_types_v<std::remove_const_t<Components>...>, "all component types must be different");
 
+		++method_recursion_depth_;
+
 		using invocable_type = typename detail::filter_with<detail::batch_invocable, std::tuple<Components...>, Function>::type;
 
 		static_assert(invocable_type::value, "function is not batch-invocable with these components");
@@ -649,6 +708,10 @@ namespace psemek::ecs
 			finalize_iteration(*entry.table);
 		}
 
+		finalize_method();
+
+		--method_recursion_depth_;
+
 		return cache;
 	}
 
@@ -664,9 +727,12 @@ namespace psemek::ecs
 		auto id = next_constructor_id_++;
 
 		auto constructor_factory = [function = std::move(function)](std::vector<std::uint32_t> const & column_indices) -> detail::table_callback {
-			return [function, column_indices](container & container, detail::table & table, std::uint32_t row, util::hash_set<util::uuid> const & attached_components, bool force){
+			return [function, column_indices](container & container, detail::table & table, std::uint32_t row, util::hash_set<util::uuid> const & attached_components, util::hash_set<util::uuid> const & detached_components, bool force){
 
-				bool const invoke = force || (detail::contains_helper<Components>::contains(attached_components) || ...);
+				bool const invoke = force
+					|| (detail::contains_helper<Components>::contains(attached_components) || ...)
+					|| (detail::contains_helper<Components>::contains_without(detached_components) || ...)
+					;
 				if (!invoke)
 					return;
 
@@ -699,9 +765,11 @@ namespace psemek::ecs
 		auto id = next_destructor_id_++;
 
 		auto destructor_factory = [function = std::move(function)](std::vector<std::uint32_t> const & column_indices) -> detail::table_callback {
-			return [function, column_indices](container & container, detail::table & table, std::uint32_t row, util::hash_set<util::uuid> const & detached_components, bool force){
-
-				bool const invoke = force || (detail::contains_helper<Components>::contains(detached_components) || ...);
+			return [function, column_indices](container & container, detail::table & table, std::uint32_t row, util::hash_set<util::uuid> const & attached_components, util::hash_set<util::uuid> const & detached_components, bool force){
+				bool const invoke = force
+					|| (detail::contains_helper<Components>::contains_without(attached_components) || ...)
+					|| (detail::contains_helper<Components>::contains(detached_components) || ...)
+					;
 				if (!invoke)
 					return;
 
@@ -722,12 +790,37 @@ namespace psemek::ecs
 			entry.table->add_destructor({id, destructor_factory(entry.columns_indices)});
 	}
 
+	template <typename Function>
+	void container::finally(Function && function)
+	{
+		util::function<void(container &)> wrapper;
+		if constexpr (std::is_invocable_v<Function, container &>)
+		{
+			wrapper = std::move(function);
+		}
+		else
+		{
+			wrapper = [function = std::move(function)](ecs::container &){ function(); };
+		}
+
+		if (method_recursion_depth_ == 0)
+			wrapper(*this);
+		else
+			finally_callbacks_.push_back(std::move(wrapper));
+	}
+
 	template <typename Index, typename ... Args>
 	Index & container::index(Args && ... args)
 	{
 		return index_container_.get<Index>(*this, std::forward<Args>(args)...);
 	}
 
+	template <typename Index, typename Factory>
+	Index & container::index_factory(Factory && factory)
+	{
+		return index_container_.set<Index>(std::move(factory));
+	}
+
 	template <typename ... Components>
 	std::size_t container::memory_usage()
 	{

libs/ecs/include/psemek/ecs/detail/callback.hpp

@@ -16,7 +16,7 @@ namespace psemek::ecs
 
 		struct table;
 
-		using table_callback = util::function<void(container &, table &, std::uint32_t, util::hash_set<util::uuid> const &, bool)>;
+		using table_callback = util::function<void(container &, table &, std::uint32_t, util::hash_set<util::uuid> const &, util::hash_set<util::uuid> const &, bool)>;
 
 		struct ordered_table_callback
 		{

libs/ecs/include/psemek/ecs/detail/index_container.hpp

@@ -4,6 +4,7 @@
 #include <psemek/util/hash_table.hpp>
 #include <psemek/util/type_name.hpp>
 #include <psemek/util/exception.hpp>
+#include <psemek/util/function.hpp>
 
 #include <memory>
 
@@ -24,19 +25,21 @@ namespace psemek::ecs::detail
 		{
 			auto uuid = Index::uuid();
 			if (auto it = storage_.find(uuid); it != storage_.end())
-				return *reinterpret_cast<Index *>(it->second.get());
+				return *reinterpret_cast<Index *>(it->second());
 
 			if constexpr (std::is_constructible_v<Index, ecs::container &, Args && ...>)
 			{
-				auto ptr = std::make_shared<Index>(container, std::forward<Args>(args)...);
-				storage_.insert({uuid, ptr});
-				return *ptr;
+				auto ptr = std::make_unique<Index>(container, std::forward<Args>(args)...);
+				auto result = ptr.get();
+				storage_.insert({uuid, [ptr = std::move(ptr)]{ return ptr.get(); }});
+				return *result;
 			}
 			else if constexpr (std::is_constructible_v<Index, Args && ...>)
 			{
-				auto ptr = std::make_shared<Index>(std::forward<Args>(args)...);
-				storage_.insert({uuid, ptr});
-				return *ptr;
+				auto ptr = std::make_unique<Index>(std::forward<Args>(args)...);
+				auto result = ptr.get();
+				storage_.insert({uuid, [ptr = std::move(ptr)]{ return ptr.get(); }});
+				return *result;
 			}
 			else
 			{
@@ -44,8 +47,20 @@ namespace psemek::ecs::detail
 			}
 		}
 
+		template <typename Index, typename Factory>
+		Index & set(Factory && factory)
+		{
+			auto uuid = Index::uuid();
+			if (storage_.contains(uuid))
+				throw util::exception("Index " + util::type_name<Index>() + " is already set");
+
+			auto result = factory();
+			storage_[Index::uuid()] = std::move(factory);
+			return *result;
+		}
+
 	private:
-		util::hash_map<util::uuid, std::shared_ptr<void>> storage_;
+		util::hash_map<util::uuid, util::function<void*()>> storage_;
 	};
 
 }

libs/ecs/include/psemek/ecs/detail/table.hpp

@@ -74,10 +74,10 @@ namespace psemek::ecs::detail
 		void add_destructor(ordered_table_callback callback);
 
 		void trigger_constructors(container & container, std::uint32_t row);
-		void trigger_constructors(container & container, std::uint32_t row, util::hash_set<util::uuid> const & attached_components);
+		void trigger_constructors(container & container, std::uint32_t row, util::hash_set<util::uuid> const & attached_components, util::hash_set<util::uuid> const & detached_components);
 
 		void trigger_destructors(container & container, std::uint32_t row);
-		void trigger_destructors(container & container, std::uint32_t row, util::hash_set<util::uuid> const & detached_components);
+		void trigger_destructors(container & container, std::uint32_t row, util::hash_set<util::uuid> const & attached_components, util::hash_set<util::uuid> const & detached_components);
 
 		std::size_t memory_usage() const;
 

libs/ecs/include/psemek/ecs/detail/without.hpp

@@ -87,6 +87,12 @@ namespace psemek::ecs::detail
 		{
 			return container.contains(Component::uuid());
 		}
+
+		template <typename Container>
+		static bool contains_without(Container const &)
+		{
+			return false;
+		}
 	};
 
 	template <typename Component>
@@ -97,6 +103,12 @@ namespace psemek::ecs::detail
 		{
 			return false;
 		}
+
+		template <typename Container>
+		static bool contains_without(Container const & container)
+		{
+			return container.contains(Component::uuid());
+		}
 	};
 
 }

libs/ecs/include/psemek/ecs/dispatcher.hpp

@@ -7,6 +7,19 @@
 namespace psemek::ecs
 {
 
+	namespace detail
+	{
+
+		struct default_context_factory
+		{
+			int operator()() const noexcept
+			{
+				return 0;
+			}
+		};
+
+	}
+
 	struct dispatcher
 	{
 		dispatcher() = default;
@@ -42,12 +55,14 @@ namespace psemek::ecs
 			});
 		}
 
-		template <typename Event, typename ... Components, typename System>
-		void system(System system)
+		template <typename Event, typename ... Components, typename System, typename ContextFactory = detail::default_context_factory>
+		void system(System system, ContextFactory && context_factory = ContextFactory{})
 		{
-			handlers_[Event::uuid()].push_back([system = std::move(system), this, cache = container_->cache<Components...>()](void const * event_ptr) mutable {
+			handlers_[Event::uuid()].push_back([system = std::move(system), this, cache = container_->cache<Components...>(), context_factory = std::move(context_factory)](void const * event_ptr) mutable {
 				auto const & event = *reinterpret_cast<Event const *>(event_ptr);
 
+				[[maybe_unused]] auto context = context_factory();
+
 				container_->apply<Components...>([&]<typename ... FilteredComponents>(ecs::container & container, ecs::handle entity, FilteredComponents & ... filtered_components){
 					if constexpr (std::invocable<System, Event const &, ecs::container &, handle, FilteredComponents & ...>)
 					{

libs/ecs/include/psemek/ecs/handle.hpp

@@ -3,6 +3,7 @@
 #include <psemek/ecs/detail/id.hpp>
 
 #include <iostream>
+#include <format>
 
 namespace psemek::ecs
 {
@@ -46,4 +47,18 @@ namespace std
 		}
 	};
 
+	template <typename Char>
+	struct formatter<::psemek::ecs::handle, Char>
+	{
+		constexpr auto parse(std::format_parse_context & ctx)
+		{
+			return ctx.begin();
+		}
+
+		auto format(::psemek::ecs::handle const & handle, std::format_context & ctx) const
+		{
+			return std::format_to(ctx.out(), "({},{})", handle.id, handle.epoch);
+		}
+	};
+
 }

libs/ecs/source/container.cpp

@@ -1,6 +1,8 @@
 #include <psemek/ecs/container.hpp>
 #include <psemek/util/assert.hpp>
 
+#include <algorithm>
+
 namespace psemek::ecs
 {
 
@@ -11,13 +13,36 @@ namespace psemek::ecs
 
 	void container::destroy(handle entity)
 	{
+		#ifdef PSEMEK_DEBUG
+			assert(!currently_changing_archetype_.contains(entity));
+			currently_changing_archetype_.insert(entity);
+		#endif
+
 		assert(alive(entity));
 
+		++method_recursion_depth_;
+
 		auto const data = entity_list_.get_entities()[entity.id];
 		data.table->trigger_destructors(*this, data.row);
 
 		do_destroy(entity);
 		entity_list_.destroy(entity.id);
+
+		#ifdef PSEMEK_DEBUG
+			currently_changing_archetype_.erase(entity);
+		#endif
+
+		finalize_method();
+
+		--method_recursion_depth_;
+	}
+
+	void container::destroy_finally(handle entity)
+	{
+		finally([entity](container & world){
+			if (world.alive(entity))
+				world.destroy(entity);
+		});
 	}
 
 	accessor container::get(handle entity)
@@ -50,6 +75,8 @@ namespace psemek::ecs
 
 	std::optional<handle> container::try_clone(handle entity)
 	{
+		++method_recursion_depth_;
+
 		auto const data = entity_list_.get_entities()[entity.id];
 		if (!data.table->non_copyable_components().empty())
 			return std::nullopt;
@@ -66,6 +93,10 @@ namespace psemek::ecs
 
 		table->trigger_constructors(*this, row);
 
+		finalize_method();
+
+		--method_recursion_depth_;
+
 		return handle;
 	}
 
@@ -123,6 +154,21 @@ namespace psemek::ecs
 		}
 	}
 
+	void container::finalize_method()
+	{
+		if (method_recursion_depth_ > 1)
+			return;
+
+		while (!finally_callbacks_.empty())
+		{
+			auto callbacks = std::move(finally_callbacks_);
+			finally_callbacks_.clear();
+
+			for (auto & callback : callbacks)
+				callback(*this);
+		}
+	}
+
 	std::size_t container::entity_count()
 	{
 		return entity_list_.size();

libs/ecs/source/detail/table.cpp

@@ -1,6 +1,8 @@
 #include <psemek/ecs/detail/table.hpp>
 #include <psemek/ecs/detail/unordered_component_hash.hpp>
 
+#include <algorithm>
+
 namespace psemek::ecs::detail
 {
 
@@ -159,25 +161,25 @@ namespace psemek::ecs::detail
 	void table::trigger_constructors(container & container, std::uint32_t row)
 	{
 		for (auto const & callback : *constructors_)
-			callback.callback(container, *this, row, column_uuid_set_, true);
+			callback.callback(container, *this, row, column_uuid_set_, {}, true);
 	}
 
-	void table::trigger_constructors(container & container, std::uint32_t row, util::hash_set<util::uuid> const & attached_components)
+	void table::trigger_constructors(container & container, std::uint32_t row, util::hash_set<util::uuid> const & attached_components, util::hash_set<util::uuid> const & detached_components)
 	{
 		for (auto const & callback : *constructors_)
-			callback.callback(container, *this, row, attached_components, false);
+			callback.callback(container, *this, row, attached_components, detached_components, false);
 	}
 
 	void table::trigger_destructors(container & container, std::uint32_t row)
 	{
 		for (auto const & callback : *destructors_)
-			callback.callback(container, *this, row, column_uuid_set_, true);
+			callback.callback(container, *this, row, {}, column_uuid_set_, true);
 	}
 
-	void table::trigger_destructors(container & container, std::uint32_t row, util::hash_set<util::uuid> const & detached_components)
+	void table::trigger_destructors(container & container, std::uint32_t row, util::hash_set<util::uuid> const & attached_components, util::hash_set<util::uuid> const & detached_components)
 	{
 		for (auto const & callback : *destructors_)
-			callback.callback(container, *this, row, detached_components, false);
+			callback.callback(container, *this, row, attached_components, detached_components, false);
 	}
 
 	std::size_t table::memory_usage() const

libs/fonts/include/psemek/fonts/font.hpp

@@ -11,7 +11,7 @@
 #include <vector>
 #include <memory>
 
-namespace psemek::fonts
+namespace psemek::fonts::inline v1
 {
 
 	enum class font_type

libs/fonts/include/psemek/fonts/font_v2.hpp

@@ -92,4 +92,7 @@ namespace psemek::fonts
 
 	std::unique_ptr<font_builder> load_freetype_font(wgpu::device device, std::filesystem::path const & path);
 
+	// Font paths are specified in order
+	std::unique_ptr<font_builder> load_combined_freetype_font(wgpu::device device, std::vector<std::filesystem::path> const & paths);
+
 }

libs/fonts/include/psemek/fonts/kerned_font.hpp

@@ -3,7 +3,7 @@
 #include <psemek/fonts/font.hpp>
 #include <psemek/fonts/bmfont.hpp>
 
-namespace psemek::fonts
+namespace psemek::fonts::inline v1
 {
 
 	struct kerned_font

libs/fonts/include/psemek/fonts/msdf_font.hpp

@@ -3,7 +3,7 @@
 #include <psemek/fonts/font.hpp>
 #include <psemek/fonts/kerned_font.hpp>
 
-namespace psemek::fonts
+namespace psemek::fonts::inline v1
 {
 
 	struct msdf_font

libs/fonts/source/default_fonts.cpp

@@ -13,7 +13,7 @@
 
 #include <sstream>
 
-namespace psemek::fonts
+namespace psemek::fonts::inline v1
 {
 
 	std::unique_ptr<font> make_default_monospace_9x12_font()

libs/fonts/source/freetype.cpp

@@ -72,11 +72,19 @@ namespace psemek::fonts
 		struct face_shared
 		{
 			wgpu::device device;
-			FT_Face face;
+
+			struct face_data
+			{
+				util::blob data;
+				FT_Face face;
+			};
+
+			std::vector<face_data> faces = {};
 
 			~face_shared()
 			{
-				FT_Done_Face(face);
+				for (auto const & face : faces)
+					FT_Done_Face(face.face);
 			}
 		};
 
@@ -88,6 +96,9 @@ namespace psemek::fonts
 				, type_(type)
 				, size_(size)
 			{
+				for (auto const & face : state_->faces)
+					ft_check_result(FT_Set_Char_Size(face.face, 0, size_ << 6, 0, 0), "Failed to select freetype face size: ");
+
 				font::shape("x", {});
 				xheight_ = -shaped_text_[0].position[1].min;
 				shaped_text_.clear();
@@ -100,7 +111,7 @@ namespace psemek::fonts
 
 			std::string_view name() const override
 			{
-				return FT_Get_Postscript_Name(state_->face);
+				return FT_Get_Postscript_Name(state_->faces[0].face);
 			}
 
 			math::vector<int, 2> size() const override
@@ -134,7 +145,7 @@ namespace psemek::fonts
 
 			struct page
 			{
-				util::array<std::uint8_t, 2> pixmap;
+				util::ndarray<std::uint8_t, 2> pixmap;
 
 				int current_row_start = 0;
 				int current_row_height = 0;
@@ -147,6 +158,12 @@ namespace psemek::fonts
 
 			std::vector<page> pages_;
 
+			struct glyph_mapping_data
+			{
+				std::uint32_t font_index;
+				std::uint32_t glyph_index;
+			};
+
 			struct glyph_data
 			{
 				int page;
@@ -155,8 +172,8 @@ namespace psemek::fonts
 				math::vector<int, 2> advance;
 			};
 
-			util::hash_map<char32_t, std::uint32_t> glyph_mapping_;
-			util::hash_map<std::uint32_t, glyph_data> glyphs_;
+			util::hash_map<char32_t, glyph_mapping_data> glyph_mapping_;
+			util::hash_map<std::uint32_t, util::hash_map<std::uint32_t, glyph_data>> glyphs_;
 
 			std::vector<shaped_glyph> shaped_text_;
 
@@ -167,37 +184,53 @@ namespace psemek::fonts
 
 				shaped_text_.clear();
 
-				auto face = state_->face;
-
-				ft_check_result(FT_Set_Char_Size(face, 0, size_ << 6, 0, 0), "Failed to select freetype face size: ");
-
 				bool need_update_pages = false;
 
 				for (char32_t ch : string)
 				{
-					std::uint32_t glyph_id;
+					glyph_mapping_data mapping{0, 0};
 
 					if (auto it = glyph_mapping_.find(ch); it != glyph_mapping_.end())
 					{
-						glyph_id = it->second;
+						mapping = it->second;
 					}
 					else
 					{
-						glyph_id = FT_Get_Char_Index(face, ch);
-						glyph_mapping_[ch] = glyph_id;
+						int glyph_index = 0;
+						for (int font_index = 0; font_index < state_->faces.size(); ++font_index)
+						{
+							glyph_index = FT_Get_Char_Index(state_->faces[font_index].face, ch);
+							if (glyph_index != 0)
+							{
+								mapping = {font_index, glyph_index};
+								break;
+							}
+						}
+
+						if (glyph_index == 0)
+						{
+							mapping.font_index = 0;
+							mapping.glyph_index = FT_Get_Char_Index(state_->faces[0].face, '?');
+						}
+
+						glyph_mapping_[ch] = mapping;
 					}
 
+					auto face = state_->faces[mapping.font_index].face;
+
 					glyph_data * data = nullptr;
 
-					if (auto it = glyphs_.find(glyph_id); it != glyphs_.end())
+					auto & glyphs = glyphs_[mapping.font_index];
+
+					if (auto it = glyphs.find(mapping.glyph_index); it != glyphs.end())
 					{
 						data = &(it->second);
 					}
 					else
 					{
-						data = &glyphs_[glyph_id];
+						data = &glyphs[mapping.glyph_index];
 
-						FT_Load_Glyph(face, glyph_id, type_ == font_type::bitmap ? FT_LOAD_TARGET_LIGHT : FT_LOAD_DEFAULT);
+						FT_Load_Glyph(face, mapping.glyph_index, type_ == font_type::bitmap ? FT_LOAD_TARGET_LIGHT : FT_LOAD_DEFAULT);
 						FT_Render_Glyph(face->glyph, type_ == font_type::bitmap ? FT_RENDER_MODE_NORMAL : FT_RENDER_MODE_SDF);
 
 						if (face->glyph->bitmap.width + 2 * padding > page_size || face->glyph->bitmap.rows + 2 * padding > page_size)
@@ -206,7 +239,7 @@ namespace psemek::fonts
 							data->part = {{{0, 0}, {0, 0}}};
 							data->offset = {0, 0};
 							data->advance = {size_, 0};
-							log::warning() << "Glyph with ID " << glyph_id << " (U+" << std::hex << (int)ch << ") is larger than font atlas page: " << face->glyph->bitmap.width << "x" << face->glyph->bitmap.rows;
+							log::warning() << "Glyph with ID " << mapping.glyph_index << " (U+" << std::hex << (int)ch << ") is larger than font atlas page: " << face->glyph->bitmap.width << "x" << face->glyph->bitmap.rows;
 						}
 						else
 						{
@@ -315,7 +348,7 @@ namespace psemek::fonts
 
 			std::string_view name() const override
 			{
-				return FT_Get_Postscript_Name(state_->face);
+				return FT_Get_Postscript_Name(state_->faces[0].face);
 			}
 
 			std::unique_ptr<font> create(font_type type, int size) override
@@ -333,15 +366,26 @@ namespace psemek::fonts
 
 	std::unique_ptr<font_builder> load_freetype_font(wgpu::device device, std::filesystem::path const & path)
 	{
-		FT_Face face;
-		auto path_str = path.string();
-		ft_check_result(FT_New_Face(ft_library(), path_str.data(), 0, &face), "Failed to load font " + path_str + ": ");
+		return load_combined_freetype_font(device, {path});
+	}
 
-		auto result = std::make_unique<freetype_font_builder>(std::make_shared<face_shared>(device, face));
+	std::unique_ptr<font_builder> load_combined_freetype_font(wgpu::device device, std::vector<std::filesystem::path> const & paths)
+	{
+		auto result = std::make_shared<face_shared>(device);
 
-		log::debug() << "Loaded font " << result->name() << " (" << path << ")";
+		for (auto const & path : paths)
+		{
+			auto font_data = io::read_full(io::file_istream{path});
 
-		return result;
+			FT_Face face = nullptr;
+			ft_check_result(FT_New_Memory_Face(ft_library(), reinterpret_cast<FT_Byte const *>(font_data.data()), font_data.size(), 0, &face), util::to_string("Failed to load font ", path, ": "));
+
+			result->faces.push_back({std::move(font_data), face});
+
+			log::debug() << "Loaded font " << FT_Get_Postscript_Name(face) << " (" << path << ")";
+		}
+
+		return std::make_unique<freetype_font_builder>(std::move(result));
 	}
 
 }

libs/fonts/source/kerned_font.cpp

@@ -3,7 +3,7 @@
 #include <psemek/util/unicode.hpp>
 #include <psemek/util/exception.hpp>
 
-namespace psemek::fonts
+namespace psemek::fonts::inline v1
 {
 
 	kerned_font::kerned_font(bmfont_data data, gfx::texture_2d atlas)

libs/fonts/source/msdf_font.cpp

@@ -1,7 +1,7 @@
 #include <psemek/fonts/msdf_font.hpp>
 #include <psemek/gfx/pixmap.hpp>
 
-namespace psemek::fonts
+namespace psemek::fonts::inline v1
 {
 
 	std::unique_ptr<font> make_msdf_font(io::istream && description, io::istream && texture)

libs/gfx/include/psemek/gfx/color.hpp

@@ -317,6 +317,77 @@ namespace psemek::gfx
 		return generic_color{dark(c.c, darkness)};
 	}
 
+	template <typename T>
+	math::vector<T, 3> rgb_to_hsv(math::vector<T, 3> const & rgb)
+	{
+		T max = std::max({rgb[0], rgb[1], rgb[2]});
+		T min = std::min({rgb[0], rgb[1], rgb[2]});
+		T delta = max - min;
+
+		T hue = T{0};
+
+		if (delta > T{0})
+		{
+			if (rgb[0] > rgb[1] && rgb[0] > rgb[2])
+			{
+				hue = (rgb[1] - rgb[2]) / delta;
+			}
+			else if (rgb[1] > rgb[2])
+			{
+				hue = 2.f + (rgb[2] - rgb[0]) / delta;
+			}
+			else
+			{
+				hue = 4.f + (rgb[0] - rgb[1]) / delta;
+			}
+		}
+
+		hue *= math::rad<T>(60);
+
+		if (hue < 0)
+			hue += math::rad<T>(360);
+
+		math::vector<T, 3> result;
+		result[0] = hue;
+		result[1] = (max > T{0}) ? (delta / max) : T{0};
+		result[2] = max;
+
+		return result;
+	}
+
+	template <typename T>
+	math::vector<T, 3> hsv_to_rgb(math::vector<T, 3> const & hsv)
+	{
+		auto c = hsv[1] * hsv[2];
+		auto x = c * (T{1} - std::abs(std::fmod(hsv[0] / math::rad<T>(60), 2.f) - 1.f));
+		auto m = hsv[2] - c;
+
+		math::vector<T, 3> result;
+
+		if (hsv[0] < math::rad<T>(60))
+			result = {c, x, T{0}};
+		else if (hsv[0] < math::rad<T>(120))
+			result = {x, c, T{0}};
+		else if (hsv[0] < math::rad<T>(180))
+			result = {T{0}, c, x};
+		else if (hsv[0] < math::rad<T>(240))
+			result = {T{0}, x, c};
+		else if (hsv[0] < math::rad<T>(300))
+			result = {x, T{0}, c};
+		else
+			result = {c, T{0}, x};
+
+		result += math::vector{m, m, m};
+
+		return result;
+	}
+
+	template <std::floating_point T, std::size_t N>
+	T luminance(math::vector<T, N> const & c)
+	{
+		return c[0] * T(0.299) + c[1] * T(0.587) + c[2] * T(0.114);
+	}
+
 	std::optional<color_rgba> parse_color(std::string_view const & text);
 
 }

libs/gfx/include/psemek/gfx/detail/stb_image.h

@@ -1,4 +1,4 @@
-/* stb_image - v2.28 - public domain image loader - http://nothings.org/stb
+/* stb_image - v2.30 - public domain image loader - http://nothings.org/stb
                                   no warranty implied; use at your own risk
 
    Do this:
@@ -48,6 +48,8 @@ LICENSE
 
 RECENT REVISION HISTORY:
 
+      2.30  (2024-05-31) avoid erroneous gcc warning
+      2.29  (2023-05-xx) optimizations
       2.28  (2023-01-29) many error fixes, security errors, just tons of stuff
       2.27  (2021-07-11) document stbi_info better, 16-bit PNM support, bug fixes
       2.26  (2020-07-13) many minor fixes
@@ -1072,8 +1074,8 @@ static int stbi__addints_valid(int a, int b)
    return a <= INT_MAX - b;
 }
 
-// returns 1 if the product of two signed shorts is valid, 0 on overflow.
-static int stbi__mul2shorts_valid(short a, short b)
+// returns 1 if the product of two ints fits in a signed short, 0 on overflow.
+static int stbi__mul2shorts_valid(int a, int b)
 {
    if (b == 0 || b == -1) return 1; // multiplication by 0 is always 0; check for -1 so SHRT_MIN/b doesn't overflow
    if ((a >= 0) == (b >= 0)) return a <= SHRT_MAX/b; // product is positive, so similar to mul2sizes_valid
@@ -3384,13 +3386,13 @@ static int stbi__decode_jpeg_header(stbi__jpeg *z, int scan)
    return 1;
 }
 
-static int stbi__skip_jpeg_junk_at_end(stbi__jpeg *j)
+static stbi_uc stbi__skip_jpeg_junk_at_end(stbi__jpeg *j)
 {
    // some JPEGs have junk at end, skip over it but if we find what looks
    // like a valid marker, resume there
    while (!stbi__at_eof(j->s)) {
-      int x = stbi__get8(j->s);
-      while (x == 255) { // might be a marker
+      stbi_uc x = stbi__get8(j->s);
+      while (x == 0xff) { // might be a marker
          if (stbi__at_eof(j->s)) return STBI__MARKER_none;
          x = stbi__get8(j->s);
          if (x != 0x00 && x != 0xff) {
@@ -4176,6 +4178,7 @@ typedef struct
 {
    stbi_uc *zbuffer, *zbuffer_end;
    int num_bits;
+   int hit_zeof_once;
    stbi__uint32 code_buffer;
 
    char *zout;
@@ -4242,9 +4245,20 @@ stbi_inline static int stbi__zhuffman_decode(stbi__zbuf *a, stbi__zhuffman *z)
    int b,s;
    if (a->num_bits < 16) {
       if (stbi__zeof(a)) {
-         return -1;   /* report error for unexpected end of data. */
+         if (!a->hit_zeof_once) {
+            // This is the first time we hit eof, insert 16 extra padding btis
+            // to allow us to keep going; if we actually consume any of them
+            // though, that is invalid data. This is caught later.
+            a->hit_zeof_once = 1;
+            a->num_bits += 16; // add 16 implicit zero bits
+         } else {
+            // We already inserted our extra 16 padding bits and are again
+            // out, this stream is actually prematurely terminated.
+            return -1;
+         }
+      } else {
+         stbi__fill_bits(a);
       }
-      stbi__fill_bits(a);
    }
    b = z->fast[a->code_buffer & STBI__ZFAST_MASK];
    if (b) {
@@ -4309,6 +4323,13 @@ static int stbi__parse_huffman_block(stbi__zbuf *a)
          int len,dist;
          if (z == 256) {
             a->zout = zout;
+            if (a->hit_zeof_once && a->num_bits < 16) {
+               // The first time we hit zeof, we inserted 16 extra zero bits into our bit
+               // buffer so the decoder can just do its speculative decoding. But if we
+               // actually consumed any of those bits (which is the case when num_bits < 16),
+               // the stream actually read past the end so it is malformed.
+               return stbi__err("unexpected end","Corrupt PNG");
+            }
             return 1;
          }
          if (z >= 286) return stbi__err("bad huffman code","Corrupt PNG"); // per DEFLATE, length codes 286 and 287 must not appear in compressed data
@@ -4320,7 +4341,7 @@ static int stbi__parse_huffman_block(stbi__zbuf *a)
          dist = stbi__zdist_base[z];
          if (stbi__zdist_extra[z]) dist += stbi__zreceive(a, stbi__zdist_extra[z]);
          if (zout - a->zout_start < dist) return stbi__err("bad dist","Corrupt PNG");
-         if (zout + len > a->zout_end) {
+         if (len > a->zout_end - zout) {
             if (!stbi__zexpand(a, zout, len)) return 0;
             zout = a->zout;
          }
@@ -4464,6 +4485,7 @@ static int stbi__parse_zlib(stbi__zbuf *a, int parse_header)
       if (!stbi__parse_zlib_header(a)) return 0;
    a->num_bits = 0;
    a->code_buffer = 0;
+   a->hit_zeof_once = 0;
    do {
       final = stbi__zreceive(a,1);
       type = stbi__zreceive(a,2);
@@ -4619,9 +4641,8 @@ enum {
    STBI__F_up=2,
    STBI__F_avg=3,
    STBI__F_paeth=4,
-   // synthetic filters used for first scanline to avoid needing a dummy row of 0s
-   STBI__F_avg_first,
-   STBI__F_paeth_first
+   // synthetic filter used for first scanline to avoid needing a dummy row of 0s
+   STBI__F_avg_first
 };
 
 static stbi_uc first_row_filter[5] =
@@ -4630,29 +4651,56 @@ static stbi_uc first_row_filter[5] =
    STBI__F_sub,
    STBI__F_none,
    STBI__F_avg_first,
-   STBI__F_paeth_first
+   STBI__F_sub // Paeth with b=c=0 turns out to be equivalent to sub
 };
 
 static int stbi__paeth(int a, int b, int c)
 {
-   int p = a + b - c;
-   int pa = abs(p-a);
-   int pb = abs(p-b);
-   int pc = abs(p-c);
-   if (pa <= pb && pa <= pc) return a;
-   if (pb <= pc) return b;
-   return c;
+   // This formulation looks very different from the reference in the PNG spec, but is
+   // actually equivalent and has favorable data dependencies and admits straightforward
+   // generation of branch-free code, which helps performance significantly.
+   int thresh = c*3 - (a + b);
+   int lo = a < b ? a : b;
+   int hi = a < b ? b : a;
+   int t0 = (hi <= thresh) ? lo : c;
+   int t1 = (thresh <= lo) ? hi : t0;
+   return t1;
 }
 
 static const stbi_uc stbi__depth_scale_table[9] = { 0, 0xff, 0x55, 0, 0x11, 0,0,0, 0x01 };
 
+// adds an extra all-255 alpha channel
+// dest == src is legal
+// img_n must be 1 or 3
+static void stbi__create_png_alpha_expand8(stbi_uc *dest, stbi_uc *src, stbi__uint32 x, int img_n)
+{
+   int i;
+   // must process data backwards since we allow dest==src
+   if (img_n == 1) {
+      for (i=x-1; i >= 0; --i) {
+         dest[i*2+1] = 255;
+         dest[i*2+0] = src[i];
+      }
+   } else {
+      STBI_ASSERT(img_n == 3);
+      for (i=x-1; i >= 0; --i) {
+         dest[i*4+3] = 255;
+         dest[i*4+2] = src[i*3+2];
+         dest[i*4+1] = src[i*3+1];
+         dest[i*4+0] = src[i*3+0];
+      }
+   }
+}
+
 // create the png data from post-deflated data
 static int stbi__create_png_image_raw(stbi__png *a, stbi_uc *raw, stbi__uint32 raw_len, int out_n, stbi__uint32 x, stbi__uint32 y, int depth, int color)
 {
-   int bytes = (depth == 16? 2 : 1);
+   int bytes = (depth == 16 ? 2 : 1);
    stbi__context *s = a->s;
    stbi__uint32 i,j,stride = x*out_n*bytes;
    stbi__uint32 img_len, img_width_bytes;
+   stbi_uc *filter_buf;
+   int all_ok = 1;
    int k;
    int img_n = s->img_n; // copy it into a local for later
 
@@ -4664,8 +4712,11 @@ static int stbi__create_png_image_raw(stbi__png *a, stbi_uc *raw, stbi__uint32 r
    a->out = (stbi_uc *) stbi__malloc_mad3(x, y, output_bytes, 0); // extra bytes to write off the end into
    if (!a->out) return stbi__err("outofmem", "Out of memory");
 
+   // note: error exits here don't need to clean up a->out individually,
+   // stbi__do_png always does on error.
    if (!stbi__mad3sizes_valid(img_n, x, depth, 7)) return stbi__err("too large", "Corrupt PNG");
    img_width_bytes = (((img_n * x * depth) + 7) >> 3);
+   if (!stbi__mad2sizes_valid(img_width_bytes, y, img_width_bytes)) return stbi__err("too large", "Corrupt PNG");
    img_len = (img_width_bytes + 1) * y;
 
    // we used to check for exact match between raw_len and img_len on non-interlaced PNGs,
@@ -4673,189 +4724,137 @@ static int stbi__create_png_image_raw(stbi__png *a, stbi_uc *raw, stbi__uint32 r
    // so just check for raw_len < img_len always.
    if (raw_len < img_len) return stbi__err("not enough pixels","Corrupt PNG");
 
+   // Allocate two scan lines worth of filter workspace buffer.
+   filter_buf = (stbi_uc *) stbi__malloc_mad2(img_width_bytes, 2, 0);
+   if (!filter_buf) return stbi__err("outofmem", "Out of memory");
+
+   // Filtering for low-bit-depth images
+   if (depth < 8) {
+      filter_bytes = 1;
+      width = img_width_bytes;
+   }
+
    for (j=0; j < y; ++j) {
-      stbi_uc *cur = a->out + stride*j;
-      stbi_uc *prior;
+      // cur/prior filter buffers alternate
+      stbi_uc *cur = filter_buf + (j & 1)*img_width_bytes;
+      stbi_uc *prior = filter_buf + (~j & 1)*img_width_bytes;
+      stbi_uc *dest = a->out + stride*j;
+      int nk = width * filter_bytes;
       int filter = *raw++;
 
-      if (filter > 4)
-         return stbi__err("invalid filter","Corrupt PNG");
-
-      if (depth < 8) {
-         if (img_width_bytes > x) return stbi__err("invalid width","Corrupt PNG");
-         cur += x*out_n - img_width_bytes; // store output to the rightmost img_len bytes, so we can decode in place
-         filter_bytes = 1;
-         width = img_width_bytes;
+      // check filter type
+      if (filter > 4) {
+         all_ok = stbi__err("invalid filter","Corrupt PNG");
+         break;
       }
-      prior = cur - stride; // bugfix: need to compute this after 'cur +=' computation above
 
       // if first row, use special filter that doesn't sample previous row
       if (j == 0) filter = first_row_filter[filter];
 
-      // handle first byte explicitly
-      for (k=0; k < filter_bytes; ++k) {
-         switch (filter) {
-            case STBI__F_none       : cur[k] = raw[k]; break;
-            case STBI__F_sub        : cur[k] = raw[k]; break;
-            case STBI__F_up         : cur[k] = STBI__BYTECAST(raw[k] + prior[k]); break;
-            case STBI__F_avg        : cur[k] = STBI__BYTECAST(raw[k] + (prior[k]>>1)); break;
-            case STBI__F_paeth      : cur[k] = STBI__BYTECAST(raw[k] + stbi__paeth(0,prior[k],0)); break;
-            case STBI__F_avg_first  : cur[k] = raw[k]; break;
-            case STBI__F_paeth_first: cur[k] = raw[k]; break;
-         }
+      // perform actual filtering
+      switch (filter) {
+      case STBI__F_none:
+         memcpy(cur, raw, nk);
+         break;
+      case STBI__F_sub:
+         memcpy(cur, raw, filter_bytes);
+         for (k = filter_bytes; k < nk; ++k)
+            cur[k] = STBI__BYTECAST(raw[k] + cur[k-filter_bytes]);
+         break;
+      case STBI__F_up:
+         for (k = 0; k < nk; ++k)
+            cur[k] = STBI__BYTECAST(raw[k] + prior[k]);
+         break;
+      case STBI__F_avg:
+         for (k = 0; k < filter_bytes; ++k)
+            cur[k] = STBI__BYTECAST(raw[k] + (prior[k]>>1));
+         for (k = filter_bytes; k < nk; ++k)
+            cur[k] = STBI__BYTECAST(raw[k] + ((prior[k] + cur[k-filter_bytes])>>1));
+         break;
+      case STBI__F_paeth:
+         for (k = 0; k < filter_bytes; ++k)
+            cur[k] = STBI__BYTECAST(raw[k] + prior[k]); // prior[k] == stbi__paeth(0,prior[k],0)
+         for (k = filter_bytes; k < nk; ++k)
+            cur[k] = STBI__BYTECAST(raw[k] + stbi__paeth(cur[k-filter_bytes], prior[k], prior[k-filter_bytes]));
+         break;
+      case STBI__F_avg_first:
+         memcpy(cur, raw, filter_bytes);
+         for (k = filter_bytes; k < nk; ++k)
+            cur[k] = STBI__BYTECAST(raw[k] + (cur[k-filter_bytes] >> 1));
+         break;
       }
 
-      if (depth == 8) {
-         if (img_n != out_n)
-            cur[img_n] = 255; // first pixel
-         raw += img_n;
-         cur += out_n;
-         prior += out_n;
-      } else if (depth == 16) {
-         if (img_n != out_n) {
-            cur[filter_bytes]   = 255; // first pixel top byte
-            cur[filter_bytes+1] = 255; // first pixel bottom byte
-         }
-         raw += filter_bytes;
-         cur += output_bytes;
-         prior += output_bytes;
-      } else {
-         raw += 1;
-         cur += 1;
-         prior += 1;
-      }
+      raw += nk;
 
-      // this is a little gross, so that we don't switch per-pixel or per-component
-      if (depth < 8 || img_n == out_n) {
-         int nk = (width - 1)*filter_bytes;
-         #define STBI__CASE(f) \
-             case f:     \
-                for (k=0; k < nk; ++k)
-         switch (filter) {
-            // "none" filter turns into a memcpy here; make that explicit.
-            case STBI__F_none:         memcpy(cur, raw, nk); break;
-            STBI__CASE(STBI__F_sub)          { cur[k] = STBI__BYTECAST(raw[k] + cur[k-filter_bytes]); } break;
-            STBI__CASE(STBI__F_up)           { cur[k] = STBI__BYTECAST(raw[k] + prior[k]); } break;
-            STBI__CASE(STBI__F_avg)          { cur[k] = STBI__BYTECAST(raw[k] + ((prior[k] + cur[k-filter_bytes])>>1)); } break;
-            STBI__CASE(STBI__F_paeth)        { cur[k] = STBI__BYTECAST(raw[k] + stbi__paeth(cur[k-filter_bytes],prior[k],prior[k-filter_bytes])); } break;
-            STBI__CASE(STBI__F_avg_first)    { cur[k] = STBI__BYTECAST(raw[k] + (cur[k-filter_bytes] >> 1)); } break;
-            STBI__CASE(STBI__F_paeth_first)  { cur[k] = STBI__BYTECAST(raw[k] + stbi__paeth(cur[k-filter_bytes],0,0)); } break;
-         }
-         #undef STBI__CASE
-         raw += nk;
-      } else {
-         STBI_ASSERT(img_n+1 == out_n);
-         #define STBI__CASE(f) \
-             case f:     \
-                for (i=x-1; i >= 1; --i, cur[filter_bytes]=255,raw+=filter_bytes,cur+=output_bytes,prior+=output_bytes) \
-                   for (k=0; k < filter_bytes; ++k)
-         switch (filter) {
-            STBI__CASE(STBI__F_none)         { cur[k] = raw[k]; } break;
-            STBI__CASE(STBI__F_sub)          { cur[k] = STBI__BYTECAST(raw[k] + cur[k- output_bytes]); } break;
-            STBI__CASE(STBI__F_up)           { cur[k] = STBI__BYTECAST(raw[k] + prior[k]); } break;
-            STBI__CASE(STBI__F_avg)          { cur[k] = STBI__BYTECAST(raw[k] + ((prior[k] + cur[k- output_bytes])>>1)); } break;
-            STBI__CASE(STBI__F_paeth)        { cur[k] = STBI__BYTECAST(raw[k] + stbi__paeth(cur[k- output_bytes],prior[k],prior[k- output_bytes])); } break;
-            STBI__CASE(STBI__F_avg_first)    { cur[k] = STBI__BYTECAST(raw[k] + (cur[k- output_bytes] >> 1)); } break;
-            STBI__CASE(STBI__F_paeth_first)  { cur[k] = STBI__BYTECAST(raw[k] + stbi__paeth(cur[k- output_bytes],0,0)); } break;
-         }
-         #undef STBI__CASE
-
-         // the loop above sets the high byte of the pixels' alpha, but for
-         // 16 bit png files we also need the low byte set. we'll do that here.
-         if (depth == 16) {
-            cur = a->out + stride*j; // start at the beginning of the row again
-            for (i=0; i < x; ++i,cur+=output_bytes) {
-               cur[filter_bytes+1] = 255;
-            }
-         }
-      }
-   }
-
-   // we make a separate pass to expand bits to pixels; for performance,
-   // this could run two scanlines behind the above code, so it won't
-   // intefere with filtering but will still be in the cache.
-   if (depth < 8) {
-      for (j=0; j < y; ++j) {
-         stbi_uc *cur = a->out + stride*j;
-         stbi_uc *in  = a->out + stride*j + x*out_n - img_width_bytes;
-         // unpack 1/2/4-bit into a 8-bit buffer. allows us to keep the common 8-bit path optimal at minimal cost for 1/2/4-bit
-         // png guarante byte alignment, if width is not multiple of 8/4/2 we'll decode dummy trailing data that will be skipped in the later loop
+      // expand decoded bits in cur to dest, also adding an extra alpha channel if desired
+      if (depth < 8) {
          stbi_uc scale = (color == 0) ? stbi__depth_scale_table[depth] : 1; // scale grayscale values to 0..255 range
+         stbi_uc *in = cur;
+         stbi_uc *out = dest;
+         stbi_uc inb = 0;
+         stbi__uint32 nsmp = x*img_n;
 
-         // note that the final byte might overshoot and write more data than desired.
-         // we can allocate enough data that this never writes out of memory, but it
-         // could also overwrite the next scanline. can it overwrite non-empty data
-         // on the next scanline? yes, consider 1-pixel-wide scanlines with 1-bit-per-pixel.
-         // so we need to explicitly clamp the final ones
-
+         // expand bits to bytes first
          if (depth == 4) {
-            for (k=x*img_n; k >= 2; k-=2, ++in) {
-               *cur++ = scale * ((*in >> 4)       );
-               *cur++ = scale * ((*in     ) & 0x0f);
+            for (i=0; i < nsmp; ++i) {
+               if ((i & 1) == 0) inb = *in++;
+               *out++ = scale * (inb >> 4);
+               inb <<= 4;
             }
-            if (k > 0) *cur++ = scale * ((*in >> 4)       );
          } else if (depth == 2) {
-            for (k=x*img_n; k >= 4; k-=4, ++in) {
-               *cur++ = scale * ((*in >> 6)       );
-               *cur++ = scale * ((*in >> 4) & 0x03);
-               *cur++ = scale * ((*in >> 2) & 0x03);
-               *cur++ = scale * ((*in     ) & 0x03);
+            for (i=0; i < nsmp; ++i) {
+               if ((i & 3) == 0) inb = *in++;
+               *out++ = scale * (inb >> 6);
+               inb <<= 2;
             }
-            if (k > 0) *cur++ = scale * ((*in >> 6)       );
-            if (k > 1) *cur++ = scale * ((*in >> 4) & 0x03);
-            if (k > 2) *cur++ = scale * ((*in >> 2) & 0x03);
-         } else if (depth == 1) {
-            for (k=x*img_n; k >= 8; k-=8, ++in) {
-               *cur++ = scale * ((*in >> 7)       );
-               *cur++ = scale * ((*in >> 6) & 0x01);
-               *cur++ = scale * ((*in >> 5) & 0x01);
-               *cur++ = scale * ((*in >> 4) & 0x01);
-               *cur++ = scale * ((*in >> 3) & 0x01);
-               *cur++ = scale * ((*in >> 2) & 0x01);
-               *cur++ = scale * ((*in >> 1) & 0x01);
-               *cur++ = scale * ((*in     ) & 0x01);
+         } else {
+            STBI_ASSERT(depth == 1);
+            for (i=0; i < nsmp; ++i) {
+               if ((i & 7) == 0) inb = *in++;
+               *out++ = scale * (inb >> 7);
+               inb <<= 1;
             }
-            if (k > 0) *cur++ = scale * ((*in >> 7)       );
-            if (k > 1) *cur++ = scale * ((*in >> 6) & 0x01);
-            if (k > 2) *cur++ = scale * ((*in >> 5) & 0x01);
-            if (k > 3) *cur++ = scale * ((*in >> 4) & 0x01);
-            if (k > 4) *cur++ = scale * ((*in >> 3) & 0x01);
-            if (k > 5) *cur++ = scale * ((*in >> 2) & 0x01);
-            if (k > 6) *cur++ = scale * ((*in >> 1) & 0x01);
          }
-         if (img_n != out_n) {
-            int q;
-            // insert alpha = 255
-            cur = a->out + stride*j;
+
+         // insert alpha=255 values if desired
+         if (img_n != out_n)
+            stbi__create_png_alpha_expand8(dest, dest, x, img_n);
+      } else if (depth == 8) {
+         if (img_n == out_n)
+            memcpy(dest, cur, x*img_n);
+         else
+            stbi__create_png_alpha_expand8(dest, cur, x, img_n);
+      } else if (depth == 16) {
+         // convert the image data from big-endian to platform-native
+         stbi__uint16 *dest16 = (stbi__uint16*)dest;
+         stbi__uint32 nsmp = x*img_n;
+
+         if (img_n == out_n) {
+            for (i = 0; i < nsmp; ++i, ++dest16, cur += 2)
+               *dest16 = (cur[0] << 8) | cur[1];
+         } else {
+            STBI_ASSERT(img_n+1 == out_n);
             if (img_n == 1) {
-               for (q=x-1; q >= 0; --q) {
-                  cur[q*2+1] = 255;
-                  cur[q*2+0] = cur[q];
+               for (i = 0; i < x; ++i, dest16 += 2, cur += 2) {
+                  dest16[0] = (cur[0] << 8) | cur[1];
+                  dest16[1] = 0xffff;
                }
             } else {
                STBI_ASSERT(img_n == 3);
-               for (q=x-1; q >= 0; --q) {
-                  cur[q*4+3] = 255;
-                  cur[q*4+2] = cur[q*3+2];
-                  cur[q*4+1] = cur[q*3+1];
-                  cur[q*4+0] = cur[q*3+0];
+               for (i = 0; i < x; ++i, dest16 += 4, cur += 6) {
+                  dest16[0] = (cur[0] << 8) | cur[1];
+                  dest16[1] = (cur[2] << 8) | cur[3];
+                  dest16[2] = (cur[4] << 8) | cur[5];
+                  dest16[3] = 0xffff;
                }
             }
          }
       }
-   } else if (depth == 16) {
-      // force the image data from big-endian to platform-native.
-      // this is done in a separate pass due to the decoding relying
-      // on the data being untouched, but could probably be done
-      // per-line during decode if care is taken.
-      stbi_uc *cur = a->out;
-      stbi__uint16 *cur16 = (stbi__uint16*)cur;
-
-      for(i=0; i < x*y*out_n; ++i,cur16++,cur+=2) {
-         *cur16 = (cur[0] << 8) | cur[1];
-      }
    }
 
+   STBI_FREE(filter_buf);
+   if (!all_ok) return 0;
+
    return 1;
 }
 
@@ -5161,9 +5160,11 @@ static int stbi__parse_png_file(stbi__png *z, int scan, int req_comp)
                // non-paletted with tRNS = constant alpha. if header-scanning, we can stop now.
                if (scan == STBI__SCAN_header) { ++s->img_n; return 1; }
                if (z->depth == 16) {
-                  for (k = 0; k < s->img_n; ++k) tc16[k] = (stbi__uint16)stbi__get16be(s); // copy the values as-is
+                  for (k = 0; k < s->img_n && k < 3; ++k) // extra loop test to suppress false GCC warning
+                     tc16[k] = (stbi__uint16)stbi__get16be(s); // copy the values as-is
                } else {
-                  for (k = 0; k < s->img_n; ++k) tc[k] = (stbi_uc)(stbi__get16be(s) & 255) * stbi__depth_scale_table[z->depth]; // non 8-bit images will be larger
+                  for (k = 0; k < s->img_n && k < 3; ++k)
+                     tc[k] = (stbi_uc)(stbi__get16be(s) & 255) * stbi__depth_scale_table[z->depth]; // non 8-bit images will be larger
                }
             }
             break;

libs/gfx/include/psemek/gfx/gltf_accessor_iterator.hpp

@@ -93,14 +93,12 @@ namespace psemek::gfx
 
 			static auto pointer(math::matrix<T, R, C> & value)
 			{
-				return &value[0][0];
+				return value.values().begin();
 			}
 
 			static void finalize(math::matrix<T, R, C> & value)
 			{
-				math::matrix<T, C, R> temp;
-				std::copy(value.coords, value.coords + R * C, temp.coords);
-				value = math::transpose(temp);
+				value = math::transpose(value);
 			}
 		};
 

libs/gfx/include/psemek/gfx/gltf_animation.hpp

@@ -1,9 +1,7 @@
 #pragma once
 
 #include <psemek/gfx/gltf_parser.hpp>
-#include <psemek/math/scale.hpp>
-#include <psemek/math/rotation.hpp>
-#include <psemek/math/translation.hpp>
+#include <psemek/math/trs.hpp>
 #include <psemek/util/exception.hpp>
 
 #include <algorithm>
@@ -147,7 +145,7 @@ namespace psemek::gfx
 
 		math::interval<float> range() const;
 
-		math::affine_transform<float, 3, 3> operator()(float time) const;
+		math::trs<float, 3> operator()(float time) const;
 
 		gltf_scale_animation const & scale() const { return scale_; }
 		gltf_rotation_animation const & rotation() const { return rotation_; }

libs/gfx/include/psemek/gfx/mesh.hpp

@@ -195,11 +195,11 @@ namespace psemek::gfx
 		GLenum index_type() const { return info_.index_type_; }
 
 		buffer & vertex_buffer() { return vertex_buffer_; }
-		buffer & index_buffer() { return index_buffer_; }
+		auto & index_buffer() { return index_buffer_; }
 		buffer & instance_buffer() { return instance_buffer_; }
 
 		buffer const & vertex_buffer() const { return vertex_buffer_; }
-		buffer const & index_buffer() const { return index_buffer_; }
+		auto const & index_buffer() const { return index_buffer_; }
 		buffer const & instance_buffer() const { return instance_buffer_; }
 
 		// Drawing commands
@@ -211,7 +211,7 @@ namespace psemek::gfx
 	private:
 		array array_;
 		buffer vertex_buffer_ = buffer::null();
-		buffer index_buffer_ = buffer::null();
+		basic_buffer<gl::ELEMENT_ARRAY_BUFFER> index_buffer_ = basic_buffer<gl::ELEMENT_ARRAY_BUFFER>::null();
 		buffer instance_buffer_ = buffer::null();
 
 		struct mesh_info
@@ -298,11 +298,8 @@ namespace psemek::gfx
 			assert((index_count % (*n)) == 0);
 
 		if (!index_buffer_)
-		{
-			index_buffer_ = buffer{};
-			array_.bind();
-			gl::BindBuffer(gl::ELEMENT_ARRAY_BUFFER, index_buffer_.id());
-		}
+			index_buffer_ = {};
+		array_.bind();
 		index_buffer_.load(indices, index_count, usage);
 		info_.index_count_ = index_count;
 		info_.indexed_ = true;

libs/gfx/include/psemek/gfx/painter.hpp

@@ -56,6 +56,7 @@ namespace psemek::gfx
 		void triangle(math::point<float, 2> const & p0, math::point<float, 2> const & p1, math::point<float, 2> const & p2, color const & c0, color const & c1, color const & c2);
 		void rect(math::box<float, 2> const & box, color const & c);
 		void circle(math::point<float, 2> const & center, float radius, color const & c, int quality = 24);
+		void circle(math::point<float, 2> const & center, float radius, color const & c0, color const & c1, int quality = 24);
 		void line(math::point<float, 2> const & p0, math::point<float, 2> const & p1, float width, color const & c, bool smooth = true);
 		void line(math::point<float, 2> const & p0, math::point<float, 2> const & p1, float w0, float w1, color const & c0, color const & c1, bool smooth = true);
 		void besier(math::point<float, 2> const & p0, math::point<float, 2> const & p1, math::point<float, 2> const & p2, float width, color const & c, int quality = 8, bool smooth = true);

libs/gfx/include/psemek/gfx/pixmap.hpp

@@ -1,6 +1,6 @@
 #pragma once
 
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/gfx/color.hpp>
 #include <psemek/io/stream.hpp>
 
@@ -10,7 +10,7 @@ namespace psemek::gfx
 {
 
 	template <typename Pixel>
-	using basic_pixmap = util::array<Pixel, 2>;
+	using basic_pixmap = util::ndarray<Pixel, 2>;
 
 	using pixmap_monochrome = basic_pixmap<std::uint8_t>;
 	using pixmap_rgb = basic_pixmap<color_rgb>;
@@ -48,7 +48,7 @@ namespace psemek::gfx
 	// Utilities
 
 	template <typename Pixel, std::size_t N>
-	auto to_srgb(util::array<Pixel, N> pm, float g = 1.f / 2.2f)
+	auto to_srgb(util::ndarray<Pixel, N> pm, float g = 1.f / 2.2f)
 	{
 		for (auto & c : pm)
 			c = to_srgb(c, g);
@@ -56,7 +56,7 @@ namespace psemek::gfx
 	}
 
 	template <typename Pixel, std::size_t N>
-	auto to_linear(util::array<Pixel, N> pm, float g = 1.f / 2.2f)
+	auto to_linear(util::ndarray<Pixel, N> pm, float g = 1.f / 2.2f)
 	{
 		return to_srgb(std::move(pm), 1.f / g);
 	}

libs/gfx/include/psemek/gfx/texture.hpp

@@ -54,10 +54,10 @@ namespace psemek::gfx
 		void load_srgb(math::vector<std::size_t, D> const & size, Pixel const * data = nullptr);
 
 		template <typename Pixel>
-		void load(util::array<Pixel, D> const & p);
+		void load(util::ndarray<Pixel, D> const & p);
 
 		template <typename Pixel>
-		void load_srgb(util::array<Pixel, D> const & p);
+		void load_srgb(util::ndarray<Pixel, D> const & p);
 
 #ifndef PSEMEK_GLES
 		void pixels(GLenum format, GLenum type, void * data, int layer = 0) const;
@@ -142,7 +142,7 @@ namespace psemek::gfx
 		}
 
 		template <typename Pixel>
-		void load(int f, util::array<Pixel, 2> const & p)
+		void load(int f, util::ndarray<Pixel, 2> const & p)
 		{
 			math::vector<std::size_t, 2> size;
 			for (std::size_t i = 0; i < 2; ++i) size[i] = p.dim(i);
@@ -392,7 +392,7 @@ namespace psemek::gfx
 
 	template <std::size_t D, GLenum Target>
 	template <typename Pixel>
-	void basic_texture<D, Target>::load(util::array<Pixel, D> const & p)
+	void basic_texture<D, Target>::load(util::ndarray<Pixel, D> const & p)
 	{
 		math::vector<std::size_t, D> size;
 		for (std::size_t i = 0; i < D; ++i) size[i] = p.dim(i);
@@ -401,7 +401,7 @@ namespace psemek::gfx
 
 	template <std::size_t D, GLenum Target>
 	template <typename Pixel>
-	void basic_texture<D, Target>::load_srgb(util::array<Pixel, D> const & p)
+	void basic_texture<D, Target>::load_srgb(util::ndarray<Pixel, D> const & p)
 	{
 		math::vector<std::size_t, D> size;
 		for (std::size_t i = 0; i < D; ++i) size[i] = p.dim(i);

libs/gfx/include/psemek/gfx/texture_atlas.hpp

@@ -19,9 +19,9 @@ namespace psemek::gfx
 		texture_2d & texture() { return texture_; }
 		texture_2d const & texture() const { return texture_; }
 
-		util::array<Pixel, 2> const & pixmap() const { return atlas_.array(); }
+		util::ndarray<Pixel, 2> const & pixmap() const { return atlas_.array(); }
 
-		std::pair<texture_view, bool> insert(Key const & key, util::array<Pixel, 2> const & data);
+		std::pair<texture_view, bool> insert(Key const & key, util::ndarray<Pixel, 2> const & data);
 		texture_view find(Key const & key) const;
 		texture_view at(Key const & key) const;
 
@@ -43,7 +43,7 @@ namespace psemek::gfx
 	{}
 
 	template <typename Pixel, typename Key, typename Compare>
-	std::pair<texture_view_2d, bool> texture_atlas_2d<Pixel, Key, Compare>::insert(Key const & key, util::array<Pixel, 2> const & data)
+	std::pair<texture_view_2d, bool> texture_atlas_2d<Pixel, Key, Compare>::insert(Key const & key, util::ndarray<Pixel, 2> const & data)
 	{
 		auto result = atlas_.insert(key, data);
 		if (result.second)

libs/gfx/source/gltf_animation.cpp

@@ -90,11 +90,9 @@ namespace psemek::gfx
 		return scale_.range() | rotation_.range() | translation_.range();
 	}
 
-	math::affine_transform<float, 3, 3> gltf_animation::operator()(float time) const
+	math::trs<float, 3> gltf_animation::operator()(float time) const
 	{
-		return math::translation(translation_(time)).transform()
-			* math::quaternion_rotation(rotation_(time)).transform()
-			* math::scale(scale_(time)).transform();
+		return {translation_(time), rotation_(time), scale_(time)};
 	}
 
 }

libs/gfx/source/mesh.cpp

@@ -144,11 +144,8 @@ namespace psemek::gfx
 		assert(vertex_buffer_);
 
 		if (!index_buffer_)
-		{
-			index_buffer_ = buffer{};
-			array_.bind();
-			gl::BindBuffer(gl::ELEMENT_ARRAY_BUFFER, index_buffer_.id());
-		}
+			index_buffer_ = {};
+		array_.bind();
 		vertex_buffer_.load(vertices, vertex_size * vertex_count, usage);
 		index_buffer_.load(indices, index_size(index_type) * index_count, usage);
 		info_.vertex_count_ = vertex_count;

libs/gfx/source/painter.cpp

@@ -216,14 +216,19 @@ namespace psemek::gfx
 	}
 
 	void painter::circle(math::point<float, 2> const & p, float r, color const & c, int quality)
+	{
+		circle(p, r, c, c, quality);
+	}
+
+	void painter::circle(math::point<float, 2> const & p, float r, color const & c0, color const & c1, int quality)
 	{
 		std::uint32_t const base = impl().vertices.size();
 
-		impl().vertices.push_back({{p[0], p[1], 0.f}, c});
+		impl().vertices.push_back({{p[0], p[1], 0.f}, c0});
 		for (int i = 0; i < quality; ++i)
 		{
 			float const a = (math::pi * 2.f * i) / quality;
-			impl().vertices.push_back({{p[0] + r * std::cos(a), p[1] + r * std::sin(a), 0.f}, c});
+			impl().vertices.push_back({{p[0] + r * std::cos(a), p[1] + r * std::sin(a), 0.f}, c1});
 		}
 
 		for (int i = 0; i < quality; ++i)
@@ -301,18 +306,37 @@ namespace psemek::gfx
 	math::vector<float, 2> painter::text_size(std::string_view str, font f)
 	{
 		// TODO: multiline text
-		math::vector<float, 2> s;
+		math::vector<float, 2> font_size;
 		switch (f)
 		{
 		case font::font_9x12:
-			s = {9.f, 12.f};
+			font_size = {9.f, 12.f};
 			break;
 		default:
 			throw util::unknown_enum_value_exception(f);
 		}
 
-		s[0] *= str.size();
-		return s;
+		int max_line_width = 0;
+		int line_count = 0;
+
+		int last_line_start = 0;
+		for (int i = 0; i < str.size(); ++i)
+		{
+			if (str[i] == '\n')
+			{
+				math::make_max(max_line_width, i - last_line_start);
+				last_line_start = i;
+				++line_count;
+			}
+		}
+
+		if (last_line_start + 1 != str.size())
+		{
+			math::make_max(max_line_width, static_cast<int>(str.size()) - last_line_start);
+			++line_count;
+		}
+
+		return {font_size[0] * max_line_width, font_size[1] * line_count};
 	}
 
 	void painter::text(math::point<float, 2> const & p, std::string_view str, text_options const & opts)
@@ -489,17 +513,17 @@ namespace psemek::gfx
 	{
 		auto const size = math::pointwise_mult(text_size(str, opts.f), opts.scale);
 
-		math::vector<float, 3> pen { 0.f, 0.f, 0.f };
+		math::vector<float, 3> origin { 0.f, 0.f, 0.f };
 
 		switch (opts.x)
 		{
 		case x_align::left:
 			break;
 		case x_align::center:
-			pen[0] -= size[0] / 2.f;
+			origin[0] -= size[0] / 2.f;
 			break;
 		case x_align::right:
-			pen[0] -= size[0];
+			origin[0] -= size[0];
 			break;
 		default:
 			throw util::unknown_enum_value_exception(opts.x);
@@ -510,15 +534,17 @@ namespace psemek::gfx
 		case y_align::top:
 			break;
 		case y_align::center:
-			pen[1] -= size[1] / 2.f;
+			origin[1] -= size[1] / 2.f;
 			break;
 		case y_align::bottom:
-			pen[1] -= size[1];
+			origin[1] -= size[1];
 			break;
 		default:
 			throw util::unknown_enum_value_exception(opts.y);
 		}
 
+		auto pen = origin;
+
 		math::vector<float, 3> const sx = {9.f * opts.scale[0], 0.f, 0.f};
 		math::vector<float, 3> const sy = {0.f, 12.f * opts.scale[1], 0.f};
 
@@ -529,6 +555,13 @@ namespace psemek::gfx
 
 		for (char c : str)
 		{
+			if (c == '\n')
+			{
+				pen[0] = origin[0];
+				pen += sy;
+				continue;
+			}
+
 			// Guard against unsigned char
 #ifdef __GNUC__
 #pragma GCC diagnostic push

libs/gfx/source/program.cpp

@@ -153,47 +153,47 @@ namespace psemek::gfx
 
 	void program::uniform_proxy::operator = (math::matrix<float, 2, 2> const & m)
 	{
-		gl::UniformMatrix2fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix2fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 2, 3> const & m)
 	{
-		gl::UniformMatrix3x2fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix3x2fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 2, 4> const & m)
 	{
-		gl::UniformMatrix4x2fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix4x2fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 3, 2> const & m)
 	{
-		gl::UniformMatrix2x3fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix2x3fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 3, 3> const & m)
 	{
-		gl::UniformMatrix3fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix3fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 3, 4> const & m)
 	{
-		gl::UniformMatrix4x3fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix4x3fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 4, 2> const & m)
 	{
-		gl::UniformMatrix2x4fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix2x4fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 4, 3> const & m)
 	{
-		gl::UniformMatrix3x4fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix3x4fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::matrix<float, 4, 4> const & m)
 	{
-		gl::UniformMatrix4fv(location_, 1, gl::TRUE, m.coords);
+		gl::UniformMatrix4fv(location_, 1, gl::TRUE, &m[0][0]);
 	}
 
 	void program::uniform_proxy::operator = (math::interval<int> const & i)

libs/gfx/source/renderer/deferred.cpp

@@ -1069,7 +1069,7 @@ void main(){}
 			camera_distance[i] = dist;
 		}
 
-		util::array<bin, 3> bins({opts.grid_size[0], opts.grid_size[1], opts.grid_size[2]});
+		util::ndarray<bin, 3> bins({opts.grid_size[0], opts.grid_size[1], opts.grid_size[2]});
 
 		math::box<float, 3> lit_bbox;
 		math::box<float, 3> casts_shadow_bbox;

libs/group/include/psemek/group/cyclic.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <psemek/group/permutation.hpp>
 #include <psemek/util/range.hpp>
 
 #include <cstddef>
@@ -55,6 +56,14 @@ namespace psemek::group
 			return repr_;
 		}
 
+		group::permutation<N, Repr> permutation() const
+		{
+			typename group::permutation<N, Repr>::repr_type repr;
+			for (Repr i = 0; i < N; ++i)
+				repr[i] = (i + repr_) % N;
+			return group::permutation<N, Repr>::from_repr(repr);
+		}
+
 		struct value_iterator
 		{
 			using difference_type = Repr;

libs/group/include/psemek/group/dihedral.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <psemek/group/cyclic.hpp>
+#include <psemek/group/permutation.hpp>
 
 #include <cstddef>
 #include <type_traits>
@@ -87,6 +88,25 @@ namespace psemek::group
 			return repr_ >= N;
 		}
 
+		group::permutation<N, Repr> permutation() const
+		{
+			if (repr_ < N)
+			{
+				typename group::permutation<N, Repr>::repr_type repr;
+				for (Repr i = 0; i < N; ++i)
+					repr[i] = (i + repr_) % N;
+				return group::permutation<N, Repr>::from_repr(repr);
+			}
+			else
+			{
+				typename group::permutation<N, Repr>::repr_type repr;
+				for (Repr i = 0; i < N; ++i)
+					repr[i] = (N - i) % N;
+
+				return rotation(repr_ % N).permutation() * group::permutation<N, Repr>::from_repr(repr);
+			}
+		}
+
 		struct value_iterator
 		{
 			using difference_type = Repr;

libs/group/include/psemek/group/permutation.hpp

@@ -0,0 +1,154 @@
+#pragma once
+
+#include <psemek/util/range.hpp>
+
+#include <cstddef>
+#include <type_traits>
+#include <limits>
+#include <array>
+#include <algorithm>
+#include <numeric>
+
+namespace psemek::group
+{
+
+	template <std::size_t N, typename Repr = std::size_t>
+	struct permutation
+	{
+		static_assert(std::is_integral_v<Repr> && std::is_unsigned_v<Repr>);
+		static_assert(std::numeric_limits<Repr>::max() > N);
+
+		using repr_type = std::array<Repr, N>;
+
+		static constexpr std::size_t size()
+		{
+			return N;
+		}
+
+		permutation() = default;
+
+		static permutation identity()
+		{
+			return permutation{};
+		}
+
+		static permutation from_repr(repr_type const & value)
+		{
+			return permutation{value};
+		}
+
+		repr_type const & value() const
+		{
+			return repr_;
+		}
+
+		repr_type operator() (repr_type const & value) const
+		{
+			repr_type result;
+			for (Repr i = 0; i < N; ++i)
+				result[repr_[i]] = value[i];
+			return result;
+		}
+
+		struct value_iterator
+		{
+			using difference_type = Repr;
+			using value_type = permutation<N, Repr>;
+			using pointer = value_type *;
+			using reference	= value_type &;
+			using iterator_category = std::forward_iterator_tag;
+
+			repr_type repr;
+			bool end;
+
+			permutation<N, Repr> operator *() const
+			{
+				return permutation<N, Repr>{repr};
+			}
+
+			value_iterator & operator++()
+			{
+				end = std::next_permutation(repr.begin(), repr.end());
+				return *this;
+			}
+
+			friend bool operator == (value_iterator const & it1, value_iterator const & it2)
+			{
+				return it1.repr == it2.repr && it1.end == it2.end;
+			}
+		};
+
+		static auto values()
+		{
+			repr_type repr;
+			std::iota(repr.begin(), repr.end(), Repr{0});
+			return util::range{value_iterator{repr, false}, value_iterator{repr, true}};
+		}
+
+	private:
+		repr_type repr_{0};
+
+		explicit permutation(repr_type const & repr)
+			: repr_(repr)
+		{}
+	};
+
+	template <std::size_t N, typename Repr>
+	bool operator == (permutation<N, Repr> const & g1, permutation<N, Repr> const & g2)
+	{
+		return g1.value() == g2.value();
+	}
+
+	template <std::size_t N, typename Repr>
+	auto operator <=> (permutation<N, Repr> const & g1, permutation<N, Repr> const & g2)
+	{
+		return g1.value() <=> g2.value();
+	}
+
+	template <std::size_t N, typename Repr>
+	permutation<N, Repr> operator * (permutation<N, Repr> const & g1, permutation<N, Repr> const & g2)
+	{
+		typename permutation<N, Repr>::repr_type repr;
+		for (Repr i = 0; i < N; ++i)
+			repr[i] = g1.value()[g2.value()[i]];
+		return permutation<N, Repr>::from_repr(repr);
+	}
+
+	template <std::size_t N, typename Repr>
+	permutation<N, Repr> inverse(permutation<N, Repr> const & g)
+	{
+		typename permutation<N, Repr>::repr_type repr;
+		for (Repr i = 0; i < N; ++i)
+			repr[g.value()[i]] = i;
+		return permutation<N, Repr>::from_repr(repr);
+	}
+
+	template <typename OStream, std::size_t N, typename Repr>
+	OStream & operator << (OStream & os, permutation<N, Repr> const & g)
+	{
+		os << 'S' << N << '(';
+		for (Repr i = 0; i < N; ++i)
+		{
+			if (i > 0)
+				os << ", ";
+			os << g.value()[i];
+		}
+		os << ')';
+		return os;
+	}
+
+	template <typename OStream, std::size_t N, typename Repr>
+	void write(OStream & out, permutation<N, Repr> const & g)
+	{
+		write(out, g.value());
+	}
+
+	template <typename IStream, std::size_t N, typename Repr>
+	void read(IStream & in, permutation<N, Repr> & g)
+	{
+		typename permutation<N, Repr>::repr_type value;
+		read(in, value);
+		g = permutation<N, Repr>::from_repr(value);
+	}
+
+}

libs/io/source/file_stream.cpp

@@ -1,27 +1,72 @@
 #include <psemek/io/file_stream.hpp>
 #include <psemek/util/system_error.hpp>
+#include <psemek/util/enum.hpp>
 
 #include <cstring>
-#include <codecvt>
 
 namespace psemek::io
 {
 
-	static void throw_fopen [[noreturn]] (std::filesystem::path const & path)
+	namespace
 	{
-		throw util::system_error(std::error_code{errno, std::system_category()}, "Failed to open " + path.string());
-	}
 
-	static FILE * safe_fopen(std::filesystem::path const & path, const char * mode)
-	{
-		std::string path_str = path.string();
-		auto f = std::fopen(path_str.c_str(), mode);
-		if (!f) throw_fopen(path);
-		return f;
+		void throw_fopen [[noreturn]] (std::filesystem::path const & path)
+		{
+			throw util::system_error(std::error_code{errno, std::system_category()}, "Failed to open " + path.string());
+		}
+
+#ifdef _WIN32
+		wchar_t const * fopen_read_mode()
+		{
+			return L"rb";
+		}
+
+		wchar_t const * fopen_write_mode(unsigned flags)
+		{
+			switch (flags)
+			{
+			case 0: return L"wb";
+			case file_ostream::append: return L"ab";
+			}
+
+			throw util::exception("Unknown file_ostream open flags");
+		}
+
+		FILE * safe_fopen(std::filesystem::path const & path, wchar_t const * mode)
+		{
+			auto f = _wfopen(path.c_str(), mode);
+			if (!f) throw_fopen(path);
+			return f;
+		}
+#else
+		char const * fopen_read_mode()
+		{
+			return "rb";
+		}
+
+		char const * fopen_write_mode(unsigned flags)
+		{
+			switch (flags)
+			{
+			case 0: return "wb";
+			case file_ostream::append: return "ab";
+			}
+
+			throw util::exception("Unknown file_ostream open flags");
+		}
+
+		FILE * safe_fopen(std::filesystem::path const & path, char const * mode)
+		{
+			auto f = std::fopen(path.c_str(), mode);
+			if (!f) throw_fopen(path);
+			return f;
+		}
+#endif
+
 	}
 
 	file_istream::file_istream(std::filesystem::path const & path)
-		: file_{safe_fopen(path.c_str(), "rb")}
+		: file_{safe_fopen(path, fopen_read_mode())}
 	{}
 
 	void file_istream::reset()
@@ -43,18 +88,8 @@ namespace psemek::io
 		return std::feof(file_) != 0;
 	}
 
-	static char const * fopen_write_mode(unsigned flags)
-	{
-		switch (flags)
-		{
-		case 0: return "wb";
-		case file_ostream::append: return "ab";
-		default: throw util::exception("Unknown file_ostream open flags");
-		}
-	}
-
 	file_ostream::file_ostream(std::filesystem::path const & path, unsigned flags)
-		: file_{safe_fopen(path.c_str(), fopen_write_mode(flags))}
+		: file_{safe_fopen(path, fopen_write_mode(flags))}
 	{}
 
 	void file_ostream::reset()

libs/journal/CMakeLists.txt

@@ -0,0 +1,6 @@
+file(GLOB_RECURSE PSEMEK_JOURNAL_HEADERS RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "include/*.hpp")
+file(GLOB_RECURSE PSEMEK_JOURNAL_SOURCES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "source/*.cpp")
+
+psemek_add_library(psemek-journal ${PSEMEK_JOURNAL_HEADERS} ${PSEMEK_JOURNAL_SOURCES})
+target_include_directories(psemek-journal PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include" ${SQLite3_INCLUDE_DIRS})
+target_link_libraries(psemek-journal PUBLIC psemek-util psemek-log)

libs/journal/include/psemek/journal/event.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <string>
+#include <vector>
+
+namespace psemek::journal
+{
+
+	struct event_metadata
+	{
+		std::string source_file;
+		int source_line;
+		std::string name;
+		std::vector<std::string> columns;
+	};
+
+	struct event_data
+	{
+		std::string time;
+		std::vector<std::string> values;
+	};
+
+	struct event
+	{
+		event_metadata metadata;
+		event_data data;
+	};
+
+	std::string current_time();
+
+}

libs/journal/include/psemek/journal/journal.hpp

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <psemek/journal/event.hpp>
+#include <psemek/util/pimpl.hpp>
+
+namespace psemek::journal
+{
+
+	struct journal
+	{
+		journal();
+		~journal();
+
+		bool enabled() const;
+		void set_enabled(bool enabled);
+
+		void log_event(event const & event);
+
+		std::vector<std::pair<event_metadata const *, event_data const *>> select();
+
+	private:
+		psemek_declare_pimpl
+	};
+
+}

libs/journal/include/psemek/journal/log_event.hpp

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <psemek/journal/journal.hpp>
+
+#include <boost/preprocessor/seq/for_each.hpp>
+
+#include <format>
+
+#define psemek_journal_log_event_extract_key(key, value) key
+#define psemek_journal_log_event_extract_value(key, value) std::format("{}", value)
+
+#define psemek_journal_log_event_single_attribute_key(r, data, elem) psemek_journal_log_event_extract_key elem,
+#define psemek_journal_log_event_single_attribute_value(r, data, elem) psemek_journal_log_event_extract_value elem,
+
+#define psemek_journal_log_event(JOURNAL, NAME, ATTRIBUTES) \
+	if ((JOURNAL).enabled()) \
+		(JOURNAL).log_event({{ \
+			.source_file = __FILE__, \
+			.source_line = __LINE__, \
+			.name = NAME, \
+			.columns = { BOOST_PP_SEQ_FOR_EACH(psemek_journal_log_event_single_attribute_key, _, ATTRIBUTES) }, \
+		}, { \
+			.time = ::psemek::journal::current_time(), \
+			.values = { BOOST_PP_SEQ_FOR_EACH(psemek_journal_log_event_single_attribute_value, _, ATTRIBUTES) }, \
+		}})

libs/journal/source/event.cpp

@@ -0,0 +1,15 @@
+#include <psemek/journal/event.hpp>
+
+#include <chrono>
+#include <format>
+
+namespace psemek::journal
+{
+
+	std::string current_time()
+	{
+		const auto now = std::chrono::high_resolution_clock::now();
+		return std::format("{:%FT%TZ}", now);
+	}
+
+}

libs/journal/source/journal.cpp

@@ -0,0 +1,100 @@
+#include <psemek/journal/journal.hpp>
+#include <psemek/util/hash_table.hpp>
+
+#include <algorithm>
+
+namespace psemek::journal
+{
+
+	struct journal::impl
+	{
+		struct table
+		{
+			event_metadata metadata;
+			mutable std::vector<event_data> events;
+		};
+
+		struct table_hash
+		{
+			std::uint64_t operator()(std::string const & name) const noexcept
+			{
+				return std::hash<std::string>{}(name);
+			}
+
+			std::uint64_t operator()(table const & table) const noexcept
+			{
+				return std::hash<std::string>{}(table.metadata.name);
+			}
+		};
+
+		struct table_equal
+		{
+			bool operator()(std::string const & name1, table const & table2) const noexcept
+			{
+				return name1 == table2.metadata.name;
+			}
+
+			bool operator()(table const & table1, table const & table2) const noexcept
+			{
+				return table1.metadata.name == table2.metadata.name;
+			}
+		};
+
+		bool enabled = true;
+		util::hash_set<table, table_hash, table_equal> tables;
+
+		void log_event(event const & event)
+		{
+			if (!enabled)
+				return;
+
+			auto table_it = tables.find(event.metadata.name);
+			if (table_it == tables.end())
+				table_it = tables.insert(table{event.metadata, {}}).first;
+			table_it->events.push_back(event.data);
+		}
+
+		std::vector<std::pair<event_metadata const *, event_data const *>> select()
+		{
+			if (!enabled)
+				return {};
+
+			std::vector<std::pair<event_metadata const *, event_data const *>> result;
+
+			for (auto const & table : tables)
+				for (auto const & event : table.events)
+					result.push_back({&table.metadata, &event});
+
+			std::sort(result.begin(), result.end(), [](auto const & e1, auto const & e2){ return e1.second->time < e2.second->time; });
+
+			return result;
+		}
+	};
+
+	journal::journal()
+		: pimpl_(make_impl())
+	{}
+
+	journal::~journal() = default;
+
+	bool journal::enabled() const
+	{
+		return impl().enabled;
+	}
+
+	void journal::set_enabled(bool enabled)
+	{
+		impl().enabled = enabled;
+	}
+
+	void journal::log_event(event const & event)
+	{
+		impl().log_event(event);
+	}
+
+	std::vector<std::pair<event_metadata const *, event_data const *>> journal::select()
+	{
+		return impl().select();
+	}
+
+}

libs/log/source/log.cpp

@@ -24,7 +24,7 @@ namespace psemek::log
 		std::atomic<std::size_t> max_thread_name_length = 3;
 
 		std::mutex thread_names_mutex;
-		util::hash_map<std::thread::id, std::string> thread_names;
+		util::hash_map<std::thread::id, std::unique_ptr<std::string>> thread_names;
 
 		std::mutex sinks_mutex;
 		std::vector<std::unique_ptr<sink>> sinks;
@@ -165,10 +165,14 @@ namespace psemek::log
 			if (it != thread_names.end())
 				throw util::exception("Thread \"" + name + "\" already registered!");
 
-			thread_names[id] = name;
+			thread_names[id] = std::make_unique<std::string>(name);
 		}
 
-		max_thread_name_length = std::max(max_thread_name_length.load(), name.size());
+		{
+			std::size_t current_max = max_thread_name_length.load();
+			while(current_max < name.size() && !max_thread_name_length.compare_exchange_weak(current_max, name.size()))
+			{}
+		}
 
 		put_message(level::info, "Thread \"" + name + "\" registered");
 	}
@@ -187,7 +191,7 @@ namespace psemek::log
 				throw util::exception(os.str());
 			}
 
-			name = std::move(it->second);
+			name = std::move(*(it->second));
 			thread_names.erase(it);
 		}
 
@@ -238,7 +242,7 @@ namespace psemek::log
 			if (it == thread_names.end())
 				thread_name = &unknown_thread_name;
 			else
-				thread_name = &(it->second);
+				thread_name = it->second.get();
 		}
 
 		put_message(l, message, *thread_name);

libs/math/CMakeLists.txt

@@ -1,9 +1,7 @@
 option(PSEMEK_ROBUST_PREDICATES "Use robust geometric predicates" OFF)
 
-if(NOT (CMAKE_SYSTEM_NAME STREQUAL Emscripten))
-	find_package(Boost REQUIRED)
-endif()
 if(PSEMEK_ROBUST_PREDICATES)
+	find_package(Boost REQUIRED CONFIG)
 	find_package(GMP REQUIRED)
 endif()
 
@@ -12,8 +10,9 @@ file(GLOB_RECURSE PSEMEK_GEOM_SOURCES "source/*.cpp")
 
 psemek_add_library(psemek-math ${PSEMEK_GEOM_HEADERS} ${PSEMEK_GEOM_SOURCES})
 target_include_directories(psemek-math PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
-target_link_libraries(psemek-math PUBLIC psemek-util psemek-group Boost::boost)
+target_link_libraries(psemek-math PUBLIC psemek-util psemek-group)
 if(PSEMEK_ROBUST_PREDICATES)
+	target_link_libraries(psemek-math PUBLIC Boost::boost)
 	target_link_libraries(psemek-math PUBLIC gmp)
 endif()
 

libs/math/include/psemek/math/box.hpp

@@ -1,10 +1,10 @@
 #pragma once
 
-#include <psemek/math/detail/array.hpp>
 #include <psemek/math/interval.hpp>
 #include <psemek/math/point.hpp>
 
 #include <iostream>
+#include <format>
 
 namespace psemek::math
 {
@@ -12,7 +12,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	struct box
 	{
-		typename detail::array<interval<T>, N>::type axes;
+		interval<T> axes[N];
 
 		using point_type = point<T, N>;
 		using vector_type = vector<T, N>;
@@ -348,3 +348,26 @@ namespace psemek::math
 	}
 
 }
+
+namespace std
+{
+
+	template <typename T, std::size_t N, typename Char>
+	struct formatter<::psemek::math::box<T, N>, Char>
+		: formatter<::psemek::math::interval<T>, Char>
+	{
+		using formatter<::psemek::math::interval<T>, Char>::parse;
+
+		template <typename FormatContext>
+		auto format(::psemek::math::box<T, N> const & b, FormatContext & ctx) const
+		{
+			for (std::size_t i = 0; i < b.dimension(); ++i)
+			{
+				if (i > 0) ctx.advance_to(std::format_to(ctx.out(), "x"));
+				ctx.advance_to(formatter<::psemek::math::interval<T>, Char>::format(b[i], ctx));
+			}
+			return ctx.out();
+		}
+	};
+
+}

libs/math/include/psemek/math/concat.hpp

@@ -1,18 +1,18 @@
 #pragma once
 
-#include <psemek/math/vector.hpp>
+#include <cstddef>
 
 namespace psemek::math
 {
 
-	template <typename T, std::size_t ... Ns>
-	auto concat(vector<T, Ns> const & ... vs)
+	template <template <typename, std::size_t> typename V, typename T, std::size_t ... Ns>
+	auto concat(V<T, Ns> const & ... vs)
 	{
-		vector<T, (Ns + ...)> result;
+		V<T, (Ns + ...)> result;
 
 		std::size_t i = 0;
 
-		auto apply = [&]<std::size_t N>(vector<T, N> const & v)
+		auto apply = [&]<std::size_t N>(V<T, N> const & v)
 		{
 			for (std::size_t j = 0; j < N;)
 				result[i++] = v[j++];

libs/math/include/psemek/math/detail/array.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <psemek/math/dynamic.hpp>
 #include <psemek/util/exception.hpp>
 
 namespace psemek::math::detail
@@ -9,14 +10,46 @@ namespace psemek::math::detail
 		: util::exception
 	{
 		empty_array_exception(util::stacktrace stacktrace = {})
-			: util::exception("Indexing an empty array", std::move(stacktrace))
+			: util::exception("Indexing into a zero-dimensional array", std::move(stacktrace))
 		{}
 	};
 
+	template <std::size_t N>
+	struct dynamic_dimensions
+	{};
+
+	template <>
+	struct dynamic_dimensions<dynamic>
+	{
+		std::size_t size = 0;
+	};
+
+	template <std::size_t N>
+	void check_dynamic_size(dynamic_dimensions<N> d1, dynamic_dimensions<N> d2)
+	{
+		if constexpr (N == dynamic)
+		{
+			if (d1.size != d2.size)
+				throw dynamic_size_mismatch(d1.size, d2.size);
+		}
+	}
+
 	template <typename T, std::size_t N>
 	struct array
 	{
-		using type = T[N];
+		struct type
+		{
+			static constexpr std::size_t size = N;
+
+			T data[N];
+
+			type(dynamic_dimensions<N>){}
+
+			dynamic_dimensions<N> dimensions() const { return {}; }
+
+			T const & operator[](std::size_t i) const { return data[i]; }
+			T & operator[](std::size_t i) { return data[i]; }
+		};
 	};
 
 	template <typename T>
@@ -24,10 +57,34 @@ namespace psemek::math::detail
 	{
 		struct type
 		{
+			static constexpr std::size_t size = 0;
+
+			type(dynamic_dimensions<0>){}
+
+			dynamic_dimensions<0> dimensions() const { return {}; }
+
 			T const & operator[](std::size_t) const { throw empty_array_exception{}; }
 			T & operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
 
-			type operator + (std::size_t) const { return *this; }
+	template <typename T>
+	struct array<T, dynamic>
+	{
+		struct type
+		{
+			std::size_t size;
+			std::unique_ptr<T[]> data;
+
+			type(dynamic_dimensions<dynamic> dimensions)
+				: size(dimensions.size)
+				, data(std::make_unique_for_overwrite<T[]>(dimensions.size))
+			{}
+
+			dynamic_dimensions<dynamic> dimensions() const { return {.size = size}; }
+
+			T const & operator[](std::size_t i) const { return data[i]; }
+			T & operator[](std::size_t i) { return data[i]; }
 		};
 	};
 

libs/math/include/psemek/math/detail/array_2d.hpp

@@ -0,0 +1,242 @@
+#pragma once
+
+#include <psemek/math/detail/array.hpp>
+
+namespace psemek::math::detail
+{
+
+	template <std::size_t R, std::size_t C>
+	struct dynamic_dimensions_2d
+	{
+		dynamic_dimensions<R> rows = {};
+		dynamic_dimensions<C> columns = {};
+	};
+
+	template <std::size_t R, std::size_t C>
+	void check_dynamic_size(dynamic_dimensions_2d<R, C> d1, dynamic_dimensions_2d<R, C> d2)
+	{
+		check_dynamic_size(d1.rows, d2.rows);
+		check_dynamic_size(d1.columns, d2.columns);
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	struct array_2d
+	{
+		struct type
+		{
+			T data[R * C];
+
+			static constexpr std::size_t rows = R;
+			static constexpr std::size_t columns = C;
+
+			type(dynamic_dimensions_2d<R, C>) {}
+
+			dynamic_dimensions_2d<R, C> dimensions() const { return {}; }
+
+			T const * operator[](std::size_t row) const { return data + C * row; }
+			T * operator[](std::size_t row) { return data + C * row; }
+		};
+	};
+
+	template <typename T, std::size_t R>
+	struct array_2d<T, R, 0>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = R;
+			static constexpr std::size_t columns = 0;
+
+			type(dynamic_dimensions_2d<R, 0>) {}
+
+			dynamic_dimensions_2d<R, 0> dimensions() const { return {}; }
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T>
+	struct array_2d<T, dynamic, 0>
+	{
+		struct type
+		{
+			std::size_t rows;
+			static constexpr std::size_t columns = 0;
+
+			type(dynamic_dimensions_2d<dynamic, 0> dimensions)
+				: rows(dimensions.rows.size)
+			{}
+
+			dynamic_dimensions_2d<dynamic, 0> dimensions() const { return {.rows = rows}; }
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T, std::size_t C>
+	struct array_2d<T, 0, C>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = 0;
+			static constexpr std::size_t columns = C;
+
+			type(dynamic_dimensions_2d<0, C>) {}
+
+			dynamic_dimensions_2d<0, C> dimensions() const { return {}; }
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T>
+	struct array_2d<T, 0, dynamic>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = 0;
+			std::size_t columns;
+
+			type(dynamic_dimensions_2d<0, dynamic> dimensions)
+				: columns(dimensions.columns.size)
+			{}
+
+			dynamic_dimensions_2d<0, dynamic> dimensions() const { return {.columns = columns}; }
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T>
+	struct array_2d<T, 0, 0>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = 0;
+			static constexpr std::size_t columns = 0;
+
+			type(dynamic_dimensions_2d<0, 0>) {}
+
+			dynamic_dimensions_2d<0, 0> dimensions() const { return {}; }
+
+			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
+			T * operator[](std::size_t) { throw empty_array_exception{}; }
+		};
+	};
+
+	template <typename T, std::size_t R>
+	struct array_2d<T, R, dynamic>
+	{
+		struct type
+		{
+			static constexpr std::size_t rows = R;
+			std::size_t columns;
+			std::unique_ptr<T[]> data;
+
+			type(dynamic_dimensions_2d<R, dynamic> dimensions)
+				: columns(dimensions.columns.size)
+				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
+			{}
+
+			dynamic_dimensions_2d<R, dynamic> dimensions() const { return {.columns = columns}; }
+
+			T * operator[] (std::size_t row)
+			{
+				return data.get() + row * columns;
+			}
+
+			T const * operator[] (std::size_t row) const
+			{
+				return data.get() + row * columns;
+			}
+
+			type copy() const
+			{
+				type result;
+				result.columns = columns;
+				result.data = std::make_unique_for_overwrite<T[]>(rows * columns);
+				std::copy(data.get(), data.get() + rows * columns, result.data.get());
+				return result;
+			}
+		};
+	};
+
+	template <typename T, std::size_t C>
+	struct array_2d<T, dynamic, C>
+	{
+		struct type
+		{
+			std::size_t rows;
+			static constexpr std::size_t columns = C;
+			std::unique_ptr<T[]> data;
+
+			type(dynamic_dimensions_2d<dynamic, C> dimensions)
+				: rows(dimensions.rows.size)
+				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
+			{}
+
+			dynamic_dimensions_2d<dynamic, C> dimensions() const { return {.rows = rows}; }
+
+			T * operator[] (std::size_t row)
+			{
+				return data.get() + row * columns;
+			}
+
+			T const * operator[] (std::size_t row) const
+			{
+				return data.get() + row * columns;
+			}
+
+			type copy() const
+			{
+				type result;
+				result.rows = rows;
+				result.data = std::make_unique_for_overwrite<T[]>(rows * columns);
+				std::copy(data.get(), data.get() + rows * columns, result.data.get());
+				return result;
+			}
+		};
+	};
+
+	template <typename T>
+	struct array_2d<T, dynamic, dynamic>
+	{
+		struct type
+		{
+			std::size_t rows;
+			std::size_t columns;
+			std::unique_ptr<T[]> data;
+
+			type(dynamic_dimensions_2d<dynamic, dynamic> dimensions)
+				: rows(dimensions.rows.size)
+				, columns(dimensions.columns.size)
+				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
+			{}
+
+			dynamic_dimensions_2d<dynamic, dynamic> dimensions() const { return {.rows = rows, .columns = columns}; }
+
+			T * operator[] (std::size_t row)
+			{
+				return data.get() + row * columns;
+			}
+
+			T const * operator[] (std::size_t row) const
+			{
+				return data.get() + row * columns;
+			}
+
+			type copy() const
+			{
+				type result;
+				result.rows = rows;
+				result.data = std::make_unique_for_overwrite<T[]>(rows * columns);
+				std::copy(data.get(), data.get() + rows * columns, result.data.get());
+				return result;
+			}
+		};
+	};
+
+}

libs/math/include/psemek/math/dual.hpp

@@ -3,6 +3,7 @@
 #include <psemek/math/vector.hpp>
 
 #include <iostream>
+#include <compare>
 
 namespace psemek::math
 {
@@ -148,6 +149,57 @@ namespace psemek::math
 		return v * inverse(d);
 	}
 
+	// Exists mostly to make max/min/clamp/etc work
+	template <typename T, std::size_t N>
+	std::partial_ordering operator <=> (dual<T, N> const & d1, dual<T, N> const & d2)
+	{
+		return d1.scalar <=> d2.scalar;
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> abs(dual<T, N> const & d)
+	{
+		return d.scalar > T{0} ? d : -d;
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> exp(dual<T, N> const & d)
+	{
+		auto e = std::exp(d.scalar);
+		return {e, e * d.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> log(dual<T, N> const & d)
+	{
+		return {std::log(d.scalar), d.delta / d.scalar};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> sqrt(dual<T, N> const & d)
+	{
+		auto s = std::sqrt(d.scalar);
+		return {s, d.delta / (T{2} * s)};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> sin(dual<T, N> const & d)
+	{
+		return {std::sin(d.scalar), std::cos(d.scalar) * d.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> cos(dual<T, N> const & d)
+	{
+		return {std::cos(d.scalar), - std::sin(d.scalar) * d.delta};
+	}
+
+	template <typename T, std::size_t N>
+	dual<T, N> pow(dual<T, N> const & d1, dual<T, N> const & d2)
+	{
+		return exp(log(d1) * d2);
+	}
+
 	template <typename T, std::size_t N>
 	std::ostream & operator << (std::ostream & os, dual<T, N> const & d)
 	{

libs/math/include/psemek/math/dynamic.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <psemek/util/exception.hpp>
+
+#include <cstddef>
+#include <format>
+
+namespace psemek::math
+{
+
+	constexpr std::size_t dynamic = static_cast<std::size_t>(-1);
+
+	struct dynamic_size_mismatch
+		: util::exception
+	{
+		dynamic_size_mismatch(std::size_t size1, std::size_t size2, util::stacktrace stacktrace = {})
+			: util::exception(std::format("Dynamic array size mismatch: {} != {}", size1, size2), std::move(stacktrace))
+		{}
+	};
+
+}

libs/math/include/psemek/math/interval.hpp

@@ -4,40 +4,73 @@
 #include <limits>
 #include <type_traits>
 #include <cmath>
+#include <format>
 
 namespace psemek::math
 {
 
+	namespace detail
+	{
+
+		template <typename T>
+		concept has_min = requires ()
+		{
+			{ T::min() } -> std::same_as<T>;
+		};
+
+		template <typename T>
+		concept has_max = requires ()
+		{
+			{ T::max() } -> std::same_as<T>;
+		};
+
+	}
+
 	// Can be specialized in client code
 	template <typename T>
 	struct limits
 	{
 		static constexpr T min()
 		{
-			if constexpr (std::is_floating_point_v<T>)
+			if constexpr (detail::has_min<T>)
+			{
+				return T::min();
+			}
+			else if constexpr (std::is_floating_point_v<T>)
 			{
 				return -std::numeric_limits<T>::infinity();
 			}
-			else
+			else if constexpr (std::is_integral_v<T>)
 			{
 				return std::numeric_limits<T>::min();
 			}
+			else
+			{
+				static_assert("unknown type");
+			}
 		}
 
 		static constexpr T max()
 		{
-			if constexpr (std::is_floating_point_v<T>)
+			if constexpr (detail::has_max<T>)
+			{
+				return T::max();
+			}
+			else if constexpr (std::is_floating_point_v<T>)
 			{
 				return std::numeric_limits<T>::infinity();
 			}
-			else
+			else if constexpr (std::is_integral_v<T>)
 			{
 				return std::numeric_limits<T>::max();
 			}
+			else
+			{
+				static_assert("unknown type");
+			}
 		}
 	};
 
-
 	template <typename T>
 	struct interval_iterator
 	{
@@ -97,7 +130,7 @@ namespace psemek::math
 
 		T center() const
 		{
-			return min + (max - min) / 2;
+			return min + (max - min) / T(2);
 		}
 
 		using iterator = interval_iterator<T>;
@@ -302,3 +335,25 @@ namespace psemek::math
 	}
 
 }
+
+namespace std
+{
+
+	template <typename T, typename Char>
+	struct formatter<::psemek::math::interval<T>, Char>
+		: formatter<T, Char>
+	{
+		using formatter<T, Char>::parse;
+
+		template <typename FormatContext>
+		auto format(::psemek::math::interval<T> const & i, FormatContext & ctx) const
+		{
+			ctx.advance_to(std::format_to(ctx.out(), "["));
+			ctx.advance_to(formatter<T, Char>::format(i.min, ctx));
+			ctx.advance_to(std::format_to(ctx.out(), " .. "));
+			ctx.advance_to(formatter<T, Char>::format(i.max, ctx));
+			return std::format_to(ctx.out(), "]");
+		}
+	};
+
+}

libs/math/include/psemek/math/math.hpp

@@ -35,13 +35,13 @@ namespace psemek::math
 	template <typename T>
 	T deg(T x)
 	{
-		return static_cast<T>((180 * x) / pi);
+		return static_cast<T>((180 * x) / static_cast<T>(pi));
 	}
 
 	template <typename T>
 	T rad(T x)
 	{
-		return static_cast<T>((x * pi) / 180);
+		return static_cast<T>((x * static_cast<T>(pi)) / 180);
 	}
 
 	template <typename X, typename T>
@@ -128,19 +128,22 @@ namespace psemek::math
 		return std::pair{x1, x2};
 	}
 
+	// Moves the angle to the range [-pi, pi]
+	template <typename T>
+	T normalize_angle(T a)
+	{
+		while (a < -static_cast<T>(pi))
+			a += static_cast<T>(2 * pi);
+		while (a > static_cast<T>(pi))
+			a -= static_cast<T>(2 * pi);
+		return a;
+	}
+
 	// returns (a1 - a0)
 	template <typename T>
 	T angle_difference(T a0, T a1)
 	{
-		T const x0 = std::cos(a0);
-		T const x1 = std::cos(a1);
-		T const y0 = std::sin(a0);
-		T const y1 = std::sin(a1);
-
-		T const x = x0 * x1 + y0 * y1;
-		T const y = x0 * y1 - y0 * x1;
-
-		return std::atan2(y, x);
+		return normalize_angle(a1 - a0);
 	}
 
 	template <typename T>
@@ -174,6 +177,12 @@ namespace psemek::math
 		return (x / y) + ((x % y) == 0 ? 0 : 1);
 	}
 
+	template <typename T>
+	T fmod(T x, T m)
+	{
+		return (x >= 0) ? std::fmod(x, m) : (m - std::fmod(-x, m));
+	}
+
 	template <typename T>
 	bool make_min(T & target, T const & source)
 	{
@@ -210,13 +219,20 @@ namespace psemek::math
 	template <typename T>
 	T smoothstep(T x)
 	{
-		return x * x * (3 - 2 * x);
+		return x * x * (T{3} - T{2} * x);
 	}
 
 	template <typename T>
 	T smootherstep(T x)
 	{
-		return x * x * x * (10 - x * (15 - 6 * x));
+		auto x2 = x * x;
+		return x * x2 * (T{10} - x * T{15} + x2 * T{6});
+	}
+
+	template <typename T>
+	T inverse_smoothstep(T x)
+	{
+		return T{1} / T{2} - std::sin(std::asin(T{1} - T{2} * x) / T{3});
 	}
 
 }

libs/math/include/psemek/math/matrix.hpp

@@ -1,8 +1,9 @@
 #pragma once
 
-#include <psemek/math/detail/array.hpp>
+#include <psemek/math/detail/array_2d.hpp>
 #include <psemek/math/vector.hpp>
 #include <psemek/math/math.hpp>
+#include <psemek/util/span.hpp>
 
 #include <iostream>
 #include <iomanip>
@@ -18,26 +19,64 @@ namespace psemek::math
 		static constexpr std::size_t static_rows = R;
 		static constexpr std::size_t static_columns = C;
 
-		typename detail::array<T, R * C>::type coords;
+		using dynamic_dimensions_type = detail::dynamic_dimensions_2d<R, C>;
+
+		typename detail::array_2d<T, R, C>::type coords;
+
+		// Internal
+		matrix(dynamic_dimensions_type dimensions);
+
+		matrix();
+
+		template <typename ... Args>
+		matrix(Args && ... args) requires (R != dynamic && C != dynamic && sizeof...(Args) == R * C && detail::all_convertible_to<T, Args...>::value);
+
+		explicit matrix(std::size_t rows) requires(R == dynamic && C != dynamic);
+		explicit matrix(std::size_t columns) requires(R != dynamic && C == dynamic);
+		matrix(std::size_t rows, std::size_t columns) requires(R == dynamic && C == dynamic);
 
 		std::size_t rows() const
 		{
-			return R;
+			return coords.rows;
 		}
 
 		std::size_t columns() const
 		{
-			return C;
+			return coords.columns;
 		}
 
 		auto operator[](std::size_t i)
 		{
-			return coords + C * i;
+			return coords[i];
 		}
 
 		auto operator[](std::size_t i) const
 		{
-			return coords + C * i;
+			return coords[i];
+		}
+
+		util::span<T> values()
+		{
+			return {&coords[0][0], rows() * columns()};
+		}
+
+		util::span<T const> values() const
+		{
+			return {&coords[0][0], rows() * columns()};
+		}
+
+		matrix copy() const
+		{
+			if constexpr (R == dynamic || C == dynamic)
+			{
+				matrix r;
+				r.coords = coords.copy();
+				return r;
+			}
+			else
+			{
+				return *this;
+			}
 		}
 
 		matrix & operator *= (T const & s);
@@ -46,32 +85,72 @@ namespace psemek::math
 		matrix & operator += (matrix const & v);
 		matrix & operator -= (matrix const & v);
 
-		static matrix zero();
-		static matrix identity();
-		static matrix scalar(T const & s);
+		static matrix zero(dynamic_dimensions_type dimensions = {});
+		static matrix identity(dynamic_dimensions_type dimensions = {});
+		static matrix scalar(T const & s, dynamic_dimensions_type dimensions = {});
 		static matrix diagonal(vector<T, std::min(R, C)> const & d);
 	};
 
 	template <typename T, std::size_t R, std::size_t C>
-	matrix<T, R, C> matrix<T, R, C>::zero()
+	void check_dynamic_size(matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		matrix<T, R, C> m;
-		for (std::size_t i = 0; i < R * C; ++i)
-			m.coords[i] = 0;
+		check_dynamic_size(m1.coords.dimensions(), m2.coords.dimensions());
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix(dynamic_dimensions_type dimensions)
+		: coords(dimensions)
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix()
+		: matrix(dynamic_dimensions_type{})
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	template <typename ... Args>
+	matrix<T, R, C>::matrix(Args && ... args) requires (R != dynamic && C != dynamic && sizeof...(Args) == R * C && detail::all_convertible_to<T, Args...>::value)
+		: coords({})
+	{
+		auto out = values().begin();
+		((*out++ = args), ...);
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix(std::size_t rows) requires(R == dynamic && C != dynamic)
+		: coords({.rows = rows})
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix(std::size_t columns) requires(R != dynamic && C == dynamic)
+		: coords({.columns = columns})
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix(std::size_t rows, std::size_t columns) requires(R == dynamic && C == dynamic)
+		: coords({.rows = rows, .columns = columns})
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> matrix<T, R, C>::zero(dynamic_dimensions_type dimensions)
+	{
+		matrix<T, R, C> m(dimensions);
+		for (auto & v : m.values())
+			v = T(0);
 		return m;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
-	matrix<T, R, C> matrix<T, R, C>::identity()
+	matrix<T, R, C> matrix<T, R, C>::identity(dynamic_dimensions_type dimensions)
 	{
-		return scalar(T{1});
+		return scalar(T(1), dimensions);
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
-	matrix<T, R, C> matrix<T, R, C>::scalar(T const & s)
+	matrix<T, R, C> matrix<T, R, C>::scalar(T const & s, dynamic_dimensions_type dimensions)
 	{
-		matrix<T, R, C> m = zero();
-		for (std::size_t i = 0; i < std::min(R, C); ++i)
+		matrix<T, R, C> m = zero(dimensions);
+		for (std::size_t i = 0; i < std::min(m.rows(), m.columns()); ++i)
 			m[i][i] = s;
 		return m;
 	}
@@ -79,7 +158,17 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> matrix<T, R, C>::diagonal(vector<T, std::min(R, C)> const & d)
 	{
-		matrix<T, R, C> m = zero();
+		dynamic_dimensions_type dimensions;
+		if constexpr (R == dynamic)
+		{
+			dimensions.rows = d.dimension();
+		}
+		if constexpr (C == dynamic)
+		{
+			dimensions.columns = d.dimension();
+		}
+
+		matrix<T, R, C> m = zero(dimensions);
 		for (std::size_t i = 0; i < std::min(R, C); ++i)
 			m[i][i] = d[i];
 		return m;
@@ -88,21 +177,33 @@ namespace psemek::math
 	template <typename T1, typename T, std::size_t R, std::size_t C>
 	matrix<T1, R, C> cast(matrix<T, R, C> const & m)
 	{
-		matrix<T1, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = static_cast<T1>(m.coords[i]);
+		matrix<T1, R, C> r(m.coords.dimensions());
+
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = static_cast<T1>(value);
+
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	std::strong_ordering operator <=> (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		for (std::size_t i = 0; i < R * C; ++i)
+		check_dynamic_size(m1, m2);
+		auto begin1 = m1.values().begin();
+		auto begin2 = m2.values().begin();
+		auto end1 = m1.values().end();
+
+		// Intentionally don't use std::lexicographical_compare_three_way to fake strong ordering
+		// (actual ordering might be partial & weak due to NaN's and signed zeros
+		while (begin1 != end1)
 		{
-			if (m1.coords[i] < m2.coords[i])
+			if (*begin1 < *begin2)
 				return std::strong_ordering::less;
-			else if (m1.coords[i] > m2.coords[i])
+			else if (*begin1 > *begin2)
 				return std::strong_ordering::greater;
+			++begin1;
+			++begin2;
 		}
 		return std::strong_ordering::equal;
 	}
@@ -146,27 +247,30 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator * (matrix<T, R, C> const & m, T const & s)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = m.coords[i] * s;
+		matrix<T, R, C> r(m.coords.dimensions());
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = value * s;
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator * (T const & s, matrix<T, R, C> const & m)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = s * m.coords[i];
+		matrix<T, R, C> r(m.coords.dimensions());
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = s * value;
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator / (matrix<T, R, C> const & m, T const & s)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = m.coords[i] / s;
+		matrix<T, R, C> r(m.coords.dimensions());
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = value / s;
 		return r;
 	}
 
@@ -187,27 +291,34 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator - (matrix<T, R, C> const & m)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = -m.coords[i];
+		matrix<T, R, C> r(m.coords.dimensions());
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = -value;
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator + (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = m1.coords[i] + m2.coords[i];
+		check_dynamic_size(m1, m2);
+		matrix<T, R, C> r(m1.coords.dimensions());
+		auto in1 = m1.values().begin();
+		auto in2 = m2.values().begin();
+		for (auto & value : r.values())
+			value = (*in1++) + (*in2++);
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator - (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R * C; ++i)
-			r.coords[i] = m1.coords[i] - m2.coords[i];
+		check_dynamic_size(m1, m2);
+		matrix<T, R, C> r(m1.coords.dimensions());
+		auto in1 = m1.values().begin();
+		auto in2 = m2.values().begin();
+		for (auto & value : r.values())
+			value = (*in1++) - (*in2++);
 		return r;
 	}
 
@@ -228,11 +339,11 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	vector<T, R> operator * (matrix<T, R, C> const & m, vector<T, C> const & v)
 	{
-		vector<T, R> r;
-		for (std::size_t i = 0; i < R; ++i)
+		detail::check_dynamic_size(m.coords.dimensions().columns, v.coords.dimensions());
+		auto r = vector<T, R>::zero(m.coords.dimensions().rows);
+		for (std::size_t i = 0; i < m.rows(); ++i)
 		{
-			r[i] = T{};
-			for (std::size_t j = 0; j < C; ++j)
+			for (std::size_t j = 0; j < m.columns(); ++j)
 				r[i] += m[i][j] * v[j];
 		}
 		return r;
@@ -241,12 +352,12 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	vector<T, C> operator * (vector<T, R> const & v, matrix<T, R, C> const & m)
 	{
-		vector<T, C> r;
-		for (std::size_t j = 0; j < C; ++j)
+		detail::check_dynamic_size(v.coords.dimensions(), m.coords.dimensions().rows);
+		auto r = vector<T, C>::zero(m.coords.dimensions().columns);
+		for (std::size_t i = 0; i < m.rows(); ++i)
 		{
-			r[j] = T{};
-			for (std::size_t i = 0; i < R; ++i)
-				r[i] += v[j] * m[i][j];
+			for (std::size_t j = 0; j < m.columns(); ++j)
+				r[j] += v[i] * m[i][j];
 		}
 		return r;
 	}
@@ -254,14 +365,13 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t K, std::size_t C>
 	matrix<T, R, C> operator * (matrix<T, R, K> const & m1, matrix<T, K, C> const & m2)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
+		detail::check_dynamic_size(m1.coords.dimensions().columns, m2.coords.dimensions().rows);
+		auto r = matrix<T, R, C>::zero();
+		for (std::size_t i = 0; i < m1.rows(); ++i)
 		{
-			for (std::size_t j = 0; j < C; ++j)
+			for (std::size_t j = 0; j < m2.columns(); ++j)
 			{
-				r[i][j] = T{};
-
-				for (std::size_t k = 0; k < K; ++k)
+				for (std::size_t k = 0; k < m1.columns(); ++k)
 					r[i][j] += m1[i][k] * m2[k][j];
 			}
 		}
@@ -271,9 +381,9 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, C, R> transpose(matrix<T, R, C> const & m)
 	{
-		matrix<T, C, R> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
+		matrix<T, C, R> r({.rows = m.coords.dimensions().columns, .columns = m.coords.dimensions().rows});
+		for (std::size_t i = 0; i < m.rows(); ++i)
+			for (std::size_t j = 0; j < m.columns(); ++j)
 				r[j][i] = m[i][j];
 		return r;
 	}
@@ -281,32 +391,35 @@ namespace psemek::math
 	template <typename T, std::size_t N, typename ... Rows>
 	auto by_rows(vector<T, N> const & v, Rows const & ... rows)
 	{
+		(check_dynamic_size(v, rows), ...);
+
 		vector<T, N> m[] = {v, rows...};
 
-		matrix<T, sizeof...(Rows) + 1, N> result;
-		for (std::size_t i = 0; i < result.rows(); ++i)
-			for (std::size_t j = 0; j < result.columns(); ++j)
-				result[i][j] = m[i][j];
-		return result;
+		matrix<T, sizeof...(Rows) + 1, N> r({.columns = v.coords.dimensions()});
+		for (std::size_t i = 0; i < r.rows(); ++i)
+			for (std::size_t j = 0; j < r.columns(); ++j)
+				r[i][j] = m[i][j];
+		return r;
 	}
 
 	template <typename T, std::size_t N, typename ... Columns>
 	auto by_columns(vector<T, N> const & v, Columns const & ... columns)
 	{
+		(check_dynamic_size(v, columns), ...);
 		vector<T, N> m[] = {v, columns...};
 
-		matrix<T, N, sizeof...(Columns) + 1> result;
-		for (std::size_t i = 0; i < result.rows(); ++i)
-			for (std::size_t j = 0; j < result.columns(); ++j)
-				result[i][j] = m[j][i];
-		return result;
+		matrix<T, N, sizeof...(Columns) + 1> r({.rows = v.coords.dimensions()});
+		for (std::size_t i = 0; i < r.rows(); ++i)
+			for (std::size_t j = 0; j < r.columns(); ++j)
+				r[i][j] = m[j][i];
+		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	vector<T, C> row(matrix<T, R, C> const & m, std::size_t i)
 	{
-		vector<T, C> r;
-		for (std::size_t j = 0; j < C; ++j)
+		vector<T, C> r(m.coords.dimensions().columns);
+		for (std::size_t j = 0; j < m.columns(); ++j)
 			r[j] = m[i][j];
 		return r;
 	}
@@ -314,8 +427,8 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	vector<T, R> column(matrix<T, R, C> const & m, std::size_t j)
 	{
-		vector<T, R> r;
-		for (std::size_t i = 0; i < R; ++i)
+		vector<T, R> r(m.coords.dimensions().rows);
+		for (std::size_t i = 0; i < m.rows(); ++i)
 			r[i] = m[i][j];
 		return r;
 	}
@@ -324,28 +437,47 @@ namespace psemek::math
 	T frobenius_norm_sqr(matrix<T, R, C> const & m)
 	{
 		T r{0};
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
-				r += sqr(m[i][j]);
+		for (auto const & value : m.values())
+			r += sqr(value);
+		return r;
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	T linf_norm(matrix<T, R, C> const & m)
+	{
+		using std::abs;
+		using std::max;
+		T r = abs(m[0][0]);
+		for (auto const & value : m.values())
+			make_max(r, abs(value));
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	T frobenius_norm(matrix<T, R, C> const & m)
 	{
-		return std::sqrt(frobenius_norm_sqr(m));
+		T const max = linf_norm(m);
+		if (max == T(0)) return max;
+
+		T sum{};
+		for (auto const & value : m.values())
+			sum += sqr(value / max);
+
+		using std::sqrt;
+		return max * sqrt(sum);
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	T trace(matrix<T, R, C> const & m)
 	{
 		T r{0};
-		for (std::size_t i = 0; i < std::min(R, C); ++i)
+		for (std::size_t i = 0; i < std::min(m.rows(), m.columns()); ++i)
 			r += m[i][i];
 		return r;
 	}
 
 	template <std::size_t R1, std::size_t C1, std::size_t R2, std::size_t C2, typename T, std::size_t R, std::size_t C>
+	requires (R1 != dynamic && C1 != dynamic && R2 != dynamic && C2 != dynamic)
 	matrix<T, R2 - R1, C2 - C1> submatrix(matrix<T, R, C> const & m)
 	{
 		static_assert(R1 < R2);
@@ -354,13 +486,9 @@ namespace psemek::math
 		static_assert(C2 <= C);
 
 		matrix<T, R2 - R1, C2 - C1> result;
-		for (std::size_t i = 0; i < R2 - R1; ++i)
-		{
-			for (std::size_t j = 0; j < C2 - C1; ++j)
-			{
+		for (std::size_t i = 0; i < result.rows(); ++i)
+			for (std::size_t j = 0; j < result.columns(); ++j)
 				result[i][j] = m[i + R1][j + C1];
-			}
-		}
 		return result;
 	}
 
@@ -377,12 +505,12 @@ namespace psemek::math
 		return r;
 	}
 
-	template <typename T, std::size_t N>
-	matrix<T, N, N> outer_product(vector<T, N> const & v1, vector<T, N> const & v2)
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> outer_product(vector<T, R> const & v1, vector<T, C> const & v2)
 	{
-		matrix<T, N, N> r;
-		for (std::size_t i = 0; i < N; ++i)
-			for (std::size_t j = 0; j < N; ++j)
+		matrix<T, R, C> r({.rows = v1.coords.dimensions(), .columns = v2.coords.dimensions()});
+		for (std::size_t i = 0; i < r.rows(); ++i)
+			for (std::size_t j = 0; j < r.columns(); ++j)
 				r[i][j] = v1[i] * v2[j];
 		return r;
 	}
@@ -390,9 +518,9 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	std::ostream & operator << (std::ostream & os, matrix<T, R, C> const & m)
 	{
-		for (std::size_t i = 0; i < R; ++i)
+		for (std::size_t i = 0; i < m.rows(); ++i)
 		{
-			for (std::size_t j = 0; j < C; ++j)
+			for (std::size_t j = 0; j < m.columns(); ++j)
 				os << m[i][j] << ' ';
 			os << '\n';
 		}
@@ -409,9 +537,9 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	std::ostream & operator << (std::ostream & os, setw<T, R, C> const & w)
 	{
-		for (std::size_t i = 0; i < R; ++i)
+		for (std::size_t i = 0; i < w.m.rows(); ++i)
 		{
-			for (std::size_t j = 0; j < C; ++j)
+			for (std::size_t j = 0; j < w.m.columns(); ++j)
 				os << std::fixed << std::right << std::setw(w.w) << w.m[i][j] << ' ';
 			os << '\n';
 		}
@@ -419,3 +547,30 @@ namespace psemek::math
 	}
 
 }
+
+namespace std
+{
+
+	template <typename T, std::size_t R, std::size_t C, typename Char>
+	struct formatter<::psemek::math::matrix<T, R, C>, Char>
+		: formatter<T, Char>
+	{
+		using formatter<T, Char>::parse;
+
+		template <typename FormatContext>
+		auto format(::psemek::math::matrix<T, R, C> const & m, FormatContext & ctx) const
+		{
+			for (std::size_t i = 0; i < m.rows(); ++i)
+			{
+				for (std::size_t j = 0; j < m.columns(); ++j)
+				{
+					ctx.advance_to(formatter<T, Char>::format(m[i][j], ctx));
+					ctx.advance_to(std::format_to(ctx.out(), " "));
+				}
+				ctx.advance_to(std::format_to(ctx.out(), "\n"));
+			}
+			return ctx.out();
+		}
+	};
+
+}

libs/math/include/psemek/math/point.hpp

@@ -6,6 +6,7 @@
 #include <psemek/math/vector.hpp>
 
 #include <iostream>
+#include <format>
 
 namespace psemek::math
 {
@@ -18,17 +19,42 @@ namespace psemek::math
 
 		typename detail::array<T, N>::type coords;
 
-		point() = default;
+		using dynamic_dimensions_type = detail::dynamic_dimensions<N>;
+
+		// Internal
+		point(dynamic_dimensions_type dimensions)
+			: coords(dimensions)
+		{}
+
+		point()
+			: coords({})
+		{}
+
+		explicit point(std::size_t size) requires (N == dynamic)
+			: coords({.size = size})
+		{}
 
 		template <typename ... Args>
-		requires((sizeof...(Args) == N) && detail::all_convertible_to<T, Args...>::value)
+		requires(N != dynamic && sizeof...(Args) == N && detail::all_convertible_to<T, Args...>::value)
 		point(Args && ... args)
-			: coords{ static_cast<T>(std::forward<Args>(args))... }
-		{}
+			: coords({})
+		{
+			auto out = values().begin();
+			((*out++ = args), ...);
+		}
+
+		template <typename ... Args>
+		requires(N == dynamic && detail::all_convertible_to<T, Args...>::value)
+		point(Args && ... args)
+			: coords({.size = sizeof...(Args)})
+		{
+			auto out = values().begin();
+			((*out++ = args), ...);
+		}
 
 		std::size_t dimension() const
 		{
-			return N;
+			return coords.size;
 		}
 
 		T & operator[](std::size_t i)
@@ -41,29 +67,65 @@ namespace psemek::math
 			return coords[i];
 		}
 
+		util::span<T> values()
+		{
+			return {&coords[0], dimension()};
+		}
+
+		util::span<T const> values() const
+		{
+			return {&coords[0], dimension()};
+		}
+
+		point copy() const
+		{
+			if constexpr (N == dynamic)
+			{
+				point result;
+				result.coords = coords.copy();
+				return result;
+			}
+			else
+			{
+				return *this;
+			}
+		}
+
 		point & operator += (vector<T, N> const & v);
 		point & operator -= (vector<T, N> const & v);
 
-		static point zero();
+		static point zero(dynamic_dimensions_type dimensions = {});
 	};
 
 	template <typename ... Args>
 	point(Args && ...) -> point<std::common_type_t<Args...>, sizeof...(Args)>;
 
 	template <typename T, std::size_t N>
-	point<T, N> point<T, N>::zero()
+	void check_dynamic_size(point<T, N> const & p1, point<T, N> const & p2)
 	{
-		point<T, N> p;
-		for (std::size_t i = 0; i < N; ++i)
-			p[i] = 0;
+		check_dynamic_size(p1.coords.dimensions(), p2.coords.dimensions());
+	}
+
+	template <typename T, std::size_t N>
+	void check_dynamic_size(point<T, N> const & p, vector<T, N> const & v)
+	{
+		check_dynamic_size(p.coords.dimensions(), v.coords.dimensions());
+	}
+
+	template <typename T, std::size_t N>
+	point<T, N> point<T, N>::zero(dynamic_dimensions_type dimensions)
+	{
+		point<T, N> p(dimensions);
+		for (std::size_t i = 0; i < p.dimension(); ++i)
+			p[i] = T(0);
 		return p;
 	}
 
 	template <typename T1, typename T, std::size_t N>
 	point<T1, N> cast(point<T, N> const & p)
 	{
-		point<T1, N> r;
-		for (std::size_t i = 0; i < N; ++i)
+		point<T1, N> r(p.coords.dimensions());
+		for (std::size_t i = 0; i < r.dimension(); ++i)
 			r[i] = T1(p[i]);
 		return r;
 	}
@@ -107,7 +169,9 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	std::strong_ordering operator <=> (point<T, N> const & p1, point<T, N> const & p2)
 	{
-		for (std::size_t i = 0; i < N; ++i)
+		check_dynamic_size(p1, p2);
+
+		for (std::size_t i = 0; i < p1.dimension(); ++i)
 		{
 			if (p1[i] < p2[i])
 				return std::strong_ordering::less;
@@ -120,7 +184,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	point<T, N> operator + (point<T, N> const & p, vector<T, N> const & v)
 	{
-		point<T, N> r;
+		check_dynamic_size(p, v);
+		point<T, N> r(p.coords.dimensions());
 		for (std::size_t i = 0; i < N; ++i)
 			r[i] = p[i] + v[i];
 		return r;
@@ -129,7 +194,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	point<T, N> operator + (vector<T, N> const & v, point<T, N> const & p)
 	{
-		point<T, N> r;
+		check_dynamic_size(p, v);
+		point<T, N> r(p.coords.dimensions());
 		for (std::size_t i = 0; i < N; ++i)
 			r[i] = v[i] + p[i];
 		return r;
@@ -138,7 +204,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	point<T, N> operator - (point<T, N> const & p, vector<T, N> const & v)
 	{
-		point<T, N> r;
+		check_dynamic_size(p, v);
+		point<T, N> r(p.coords.dimensions());
 		for (std::size_t i = 0; i < N; ++i)
 			r[i] = p[i] - v[i];
 		return r;
@@ -147,7 +214,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> operator - (point<T, N> const & p1, point<T, N> const & p2)
 	{
-		vector<T, N> r;
+		check_dynamic_size(p1, p2);
+		vector<T, N> r(p1.coords.dimensions());
 		for (std::size_t i = 0; i < N; ++i)
 			r[i] = p1[i] - p2[i];
 		return r;
@@ -198,8 +266,14 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	std::ostream & operator << (std::ostream & os, point<T, N> const & p)
 	{
+		if constexpr (N == 0)
+		{
+			os << "()";
+			return os;
+		}
+
 		os << '(' << p[0];
-		for (std::size_t i = 1; i < N; ++i)
+		for (std::size_t i = 1; i < p.dimension(); ++i)
 			os << ", " << p[i];
 		os << ')';
 		return os;
@@ -208,7 +282,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	bool isfinite(point<T, N> const & p)
 	{
-		for (std::size_t i = 0; i < N; ++i)
+		for (std::size_t i = 0; i < p.dimension(); ++i)
 			if (!std::isfinite(p[i]))
 				return false;
 		return true;
@@ -222,14 +296,40 @@ namespace std
 	template <typename T, std::size_t N>
 	struct hash<::psemek::math::point<T, N>>
 	{
-		std::uint64_t operator()(::psemek::math::point<T, N> const & v) const noexcept
+		std::uint64_t operator()(::psemek::math::point<T, N> const & p) const noexcept
 		{
 			hash<T> h;
 			std::uint64_t r = 0;
 			for (std::size_t i = 0; i < N; ++i)
-				::psemek::util::hash_combine(r, h(v[i]));
+				::psemek::util::hash_combine(r, h(p[i]));
 			return r;
 		}
 	};
 
+	template <typename T, std::size_t N, typename Char>
+	struct formatter<::psemek::math::point<T, N>, Char>
+		: formatter<T, Char>
+	{
+		using formatter<T, Char>::parse;
+
+		template <typename FormatContext>
+		auto format(::psemek::math::point<T, N> const & p, FormatContext & ctx) const
+		{
+			if constexpr (N == 0)
+			{
+				return std::format_to(ctx.out(), "()");
+			}
+			else
+			{
+				ctx.advance_to(std::format_to(ctx.out(), "("));
+				for (std::size_t i = 0; i < p.dimension(); ++i)
+				{
+					if (i > 0) ctx.advance_to(std::format_to(ctx.out(), ", "));
+					ctx.advance_to(formatter<T, Char>::format(p[i], ctx));
+				}
+				return std::format_to(ctx.out(), ")");
+			}
+		}
+	};
+
 }

libs/math/include/psemek/math/quaternion.hpp

@@ -5,6 +5,8 @@
 #include <psemek/math/matrix.hpp>
 #include <psemek/math/interval.hpp>
 
+#include <format>
+
 namespace psemek::math
 {
 
@@ -256,27 +258,24 @@ namespace psemek::math
 	template <typename T>
 	quaternion<T> slerp(quaternion<T> const & q0, quaternion<T> const & q1, T const & t)
 	{
-		// threshold is chosen so that for abs(x) < threshold the second term in
-		// sin(x) Taylor series is less than the minimum value representable by T
-		static auto const threshold = std::pow(6 * std::numeric_limits<T>::min(), T{1}/T{3});
+		using std::sin;
+		using std::acos;
+		using std::abs;
 
-		auto const d = clamp(dot(normalized(q0.coords), normalized(q1.coords)), {T(-1), T(1)});
-		auto const omega = std::acos(std::abs(d));
+		auto const d = dot(q0.coords, q1.coords);
+
+		// prevent division by zero
+		if (d >= T{1})
+			return quaternion<T>{lerp(q0.coords, q1.coords, t)};
+
+		auto const omega = acos(abs(d));
 
 		// NB: the case of omega ~ pi is ambiguous and isn't handled in any special way
 
-		if (std::abs(omega) < threshold)
-		{
-			// prevent division by zero
-			return quaternion<T>{normalized(lerp(q0.coords, q1.coords, t))};
-		}
-		else
-		{
-			auto const s = std::sin(omega);
-			auto const w0 = std::sin((1 - t) * omega) / s;
-			auto const w1 = std::sin(t * omega) / s * ((d > 0) ? 1 : -1);
-			return quaternion<T>{q0.coords * w0 + q1.coords * w1};
-		}
+		auto const s = sin(omega);
+		auto const w0 = sin((1 - t) * omega) / s;
+		auto const w1 = sin(t * omega) / s * ((d > T{0}) ? T{1} : -T{1});
+		return quaternion<T>{q0.coords * w0 + q1.coords * w1};
 	}
 
 	template <typename T>
@@ -338,3 +337,22 @@ namespace psemek::math
 	}
 
 }
+
+namespace std
+{
+
+	template <typename T, typename Char>
+	struct formatter<::psemek::math::quaternion<T>, Char>
+		: formatter<::psemek::math::vector<T, 4>>
+	{
+		using formatter<::psemek::math::vector<T, 4>>::parse;
+
+		template <typename FormatContext>
+		auto format(::psemek::math::quaternion<T> const & q, FormatContext & ctx) const
+		{
+			formatter<::psemek::math::vector<T, 4>>::format(q.coords, ctx);
+			return ctx.out();
+		}
+	};
+
+}

libs/math/include/psemek/math/ray.hpp

@@ -58,6 +58,12 @@ namespace psemek::math
 		return res;
 	}
 
+	template <typename T1, typename T, std::size_t N>
+	ray<T1, N> cast(ray<T, N> const & r)
+	{
+		return {cast<T1>(r.origin), cast<T1>(r.direction)};
+	}
+
 	template <typename T, std::size_t N>
 	std::ostream & operator << (std::ostream & os, ray<T, N> const & r)
 	{

libs/math/include/psemek/math/swizzle.hpp

@@ -24,7 +24,7 @@ namespace psemek::math
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Warray-bounds"
 
-	template <std::size_t ... I>
+	template <int ... I>
 	constexpr auto swizzle = [](auto const & x){
 		static_assert(((I == -1 || I < std::decay_t<decltype(x)>::static_dimension) && ...));
 		using result = typename detail::swizzle_rebind<std::decay_t<decltype(x)>>::template type<sizeof...(I)>;

libs/math/include/psemek/math/trs.hpp

@@ -0,0 +1,88 @@
+#pragma once
+
+#include <psemek/math/scale.hpp>
+#include <psemek/math/rotation.hpp>
+#include <psemek/math/translation.hpp>
+
+namespace psemek::math
+{
+
+	// NB: composition of TRS triplets is not defined as the result
+	//     might not be expressible in TRS form
+
+	template <typename T, std::size_t N>
+	struct trs;
+
+	template <typename T>
+	struct trs<T, 3>
+	{
+		vector<T, 3> translation = {T{0}, T{0}, T{0}};
+		quaternion<T> rotation = quaternion<T>::identity();
+		vector<T, 3> scale = {T{1}, T{1}, T{1}};
+
+		static trs<T, 3> identity()
+		{
+			return {};
+		}
+
+		static trs<T, 3> from(matrix<T, 3, 4> const & m)
+		{
+			trs<T, 3> result;
+
+			result.translation[0] = m[0][3];
+			result.translation[1] = m[1][3];
+			result.translation[2] = m[2][3];
+
+			result.scale[0] = length(vector{m[0][0], m[1][0], m[2][0]});
+			result.scale[1] = length(vector{m[0][1], m[1][1], m[2][1]});
+			result.scale[2] = length(vector{m[0][2], m[1][2], m[2][2]});
+
+			matrix<T, 3, 3> r;
+			for (int i = 0; i < 3; ++i)
+				for (int j = 0; j < 3; ++j)
+					r[i][j] = m[i][j] / result.scale[j];
+
+			result.rotation = quaternion<T>::rotation(r);
+
+			return result;
+		}
+
+		matrix<T, 3, 4> affine_matrix() const
+		{
+			return transform().affine_matrix();
+		}
+
+		matrix<T, 3, 3> linear_matrix() const
+		{
+			return transform().linear_matrix();
+		}
+
+		vector<T, 3> translation_vector() const
+		{
+			return translation;
+		}
+
+		matrix<T, 4, 4> homogeneous_matrix() const
+		{
+			return transform().homogeneous_matrix();
+		}
+
+		vector<T, 3> operator()(vector<T, 3> const & v) const
+		{
+			return rotate(rotation, pointwise_mult(v, scale));
+		}
+
+		point<T, 3> operator()(point<T, 3> const & p) const
+		{
+			return p.zero() + rotate(rotation, pointwise_mult(p - p.zero(), scale)) + translation;
+		}
+
+		affine_transform<T, 3, 3> transform() const
+		{
+			return math::translation(translation).transform()
+				* quaternion_rotation(rotation).transform()
+				* math::scale(scale).transform();
+		}
+	};
+
+}