Add util, geom, cg, pcg & gfx libs

CMakeLists.txt

@@ -0,0 +1,13 @@
+cmake_minimum_required(VERSION 3.16)
+project(psemek)
+
+set(CMAKE_CXX_STANDARD 17)
+list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake/modules")
+
+string(TOUPPER ${CMAKE_BUILD_TYPE} PSEMEK_BUILD_TYPE)
+
+if(PSEMEK_BUILD_TYPE STREQUAL "DEBUG")
+	add_definitions("-DPSEMEK_DEBUG=1")
+endif()
+
+add_subdirectory(libs)

cmake/modules/FindGMP.cmake

@@ -0,0 +1,19 @@
+if(GMP_FOUND)
+	set(GMP_FIND_QUIETLY TRUE)
+endif()
+
+find_path(GMP_INCLUDE_DIRS NAMES "gmp.h" PATHS "${GMP_ROOT}/include")
+find_library(GMP_LIBRARIES NAMES "gmp" PATHS "${GMP_ROOT}/lib")
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(GMP DEFAULT_MSG GMP_INCLUDE_DIRS GMP_LIBRARIES)
+
+if(GMP_FOUND AND NOT TARGET GMP)
+	add_library(GMP SHARED IMPORTED)
+	set_target_properties(GMP PROPERTIES
+		IMPORTED_LOCATION "${GMP_LIBRARIES}"
+		INTERFACE_INCLUDE_DIRECTORIES "${GMP_INCLUDE_DIRS}"
+	)
+endif()
+
+mark_as_advanced(GMP_INCLUDE_DIRS GMP_LIBRARIES)

cmake/modules/FindSDL2.cmake

@@ -0,0 +1,19 @@
+if(SDL2_FOUND)
+	set(SDL2_FIND_QUIETLY TRUE)
+endif()
+
+find_path(SDL2_INCLUDE_DIRS NAMES "SDL2/SDL.h" PATHS "${SDL2_ROOT}/include")
+find_library(SDL2_LIBRARIES NAMES "SDL2" PATHS "${SDL2_ROOT}/lib")
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(SDL2 DEFAULT_MSG SDL2_INCLUDE_DIRS SDL2_LIBRARIES)
+
+if(SDL2_FOUND AND NOT TARGET sdl2)
+	add_library(sdl2 SHARED IMPORTED)
+	set_target_properties(sdl2 PROPERTIES
+		IMPORTED_LOCATION "${SDL2_LIBRARIES}"
+		INTERFACE_INCLUDE_DIRECTORIES "${SDL2_INCLUDE_DIRS}"
+	)
+endif()
+
+mark_as_advanced(SDL2_INCLUDE_DIRS SDL2_LIBRARIES)

libs/CMakeLists.txt

@@ -0,0 +1,5 @@
+add_subdirectory(util)
+add_subdirectory(geom)
+add_subdirectory(cg)
+add_subdirectory(pcg)
+add_subdirectory(gfx)

libs/cg/CMakeLists.txt

@@ -0,0 +1,6 @@
+file(GLOB_RECURSE PSEMEK_CG_HEADERS "include/*.hpp")
+
+add_library(cg STATIC ${PSEMEK_CG_HEADERS})
+target_include_directories(cg PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
+target_link_libraries(cg PUBLIC geom)
+set_target_properties(cg PROPERTIES LINKER_LANGUAGE CXX)

libs/cg/include/psemek/cg/bbox.hpp

@@ -0,0 +1,19 @@
+#pragma once
+
+#include <psemek/geom/box.hpp>
+
+namespace cg
+{
+
+	template <typename Iterator>
+	auto bbox (Iterator begin, Iterator end)
+	{
+		using point_type = std::decay_t<decltype(*begin)>;
+
+		box<typename point_type::scalar_type, point_type::dimension> result;
+		for (; begin != end; ++begin)
+			result |= *begin;
+		return result;
+	}
+
+}

libs/cg/include/psemek/cg/clip.hpp

@@ -0,0 +1,174 @@
+#pragma once
+
+#include <psemek/cg/dcel.hpp>
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+
+namespace cg
+{
+
+	// TODO: this was written for Voronoi and is known to fail in some cases
+	template <typename T, typename Edge, typename Face, typename Index>
+	void clip (dcel<point<T, 2>, Edge, Face, Index> & dcel, cg::vector<T, 3> const & eq)
+	{
+		auto outer_face = dcel.face(0);
+
+		auto inside = [eq](auto const & p){
+			return p[0] * eq[0] + p[1] * eq[1] + eq[2] >= 0.f;
+		};
+
+		auto cross = [&inside](auto edge){
+			return !inside(edge.origin().data()) && inside(edge.next().origin().data());
+		};
+
+		auto split_edge = [&dcel, eq](auto edge)
+		{
+			auto q0 = edge.origin().data();
+			auto q1 = edge.next().origin().data();
+
+			float t = - (eq[0] * q0[0] + eq[1] * q0[1] + eq[2]) / (eq[0] * (q1[0] - q0[0]) + eq[1] * (q1[1] - q0[1]));
+			auto q = q0 + t * (q1 - q0);
+
+			auto e = dcel.push_edge();
+			auto twin = dcel.push_edge();
+			auto p = dcel.push_point(q);
+
+			auto ctwin = edge.twin();
+
+			p.edge(e);
+
+			e.origin(p);
+			e.next(edge.next());
+			e.twin(ctwin);
+			e.face(edge.face());
+
+			twin.origin(p);
+			twin.next(ctwin.next());
+			twin.twin(edge);
+			twin.face(ctwin.face());
+
+			edge.next(e);
+			edge.twin(twin);
+
+			ctwin.next(twin);
+			ctwin.twin(e);
+		};
+
+		auto prev_edge = [](auto edge)
+		{
+			auto e = edge;
+			while (true)
+			{
+				e = e.twin();
+				if (e.next() == edge) return e;
+				e = e.next();
+			}
+		};
+
+		auto merge_edge = [&](auto edge)
+		{
+			auto dcel_copy = dcel;
+
+			auto next = edge.next();
+			auto twin = edge.twin();
+			auto ntwin = next.twin();
+
+			assert(ntwin.next() == twin);
+			assert(next.origin() == twin.origin());
+
+			edge.next(next.next());
+			edge.twin(ntwin);
+
+			ntwin.next(twin.next());
+			ntwin.twin(edge);
+
+			edge.face().edge(edge);
+			ntwin.face().edge(ntwin);
+
+			auto face = edge.face();
+			auto p = next.origin();
+			if (next > twin)
+			{
+				dcel.remove_edge(next);
+				dcel.remove_edge(twin);
+			}
+			else
+			{
+				dcel.remove_edge(twin);
+				dcel.remove_edge(next);
+			}
+			dcel.remove_point(p);
+
+			return face.edge();
+		};
+
+		auto current_edge = outer_face.edge();
+
+		bool no_intersection = true;
+
+		while (true)
+		{
+			if (cross(current_edge)) {
+				no_intersection = false;
+				break;
+			}
+
+			current_edge = current_edge.next();
+			if (current_edge == outer_face.edge())
+				break;
+		}
+
+		if (no_intersection)
+			return;
+
+		split_edge(current_edge);
+
+		while (true)
+		{
+			bool end = false;
+			while (true)
+			{
+				auto prev = current_edge.twin().next().twin();
+
+				if (prev.face() != outer_face) break;
+
+				if (cross(prev.twin()))
+				{
+					end = true;
+					split_edge(prev.twin());
+					current_edge = merge_edge(prev.next());
+
+					break;
+				}
+
+				current_edge = merge_edge(prev);
+			}
+
+			if (end) break;
+
+			for (bool finished = false; !finished;)
+			{
+				auto n = current_edge.twin().next();
+
+				if (cross(n))
+				{
+					finished = true;
+					split_edge(n);
+				}
+
+				auto prev = prev_edge(current_edge);
+
+				prev.next(n);
+				n.face().edge(n.next());
+				n.face(outer_face);
+				current_edge.origin(n.next().origin());
+				current_edge.twin().next(n.next());
+				n.next(current_edge);
+				n.origin().edge(n);
+			}
+
+			current_edge = prev_edge(current_edge);
+		}
+	}
+
+}

libs/cg/include/psemek/cg/convex_hull_2d/andrew.hpp

@@ -0,0 +1,53 @@
+#pragma once
+
+#include <psemek/geom/orientation.hpp>
+
+#include <algorithm>
+#include <numeric>
+#include <vector>
+
+namespace cg
+{
+
+	template <typename InIterator, typename OutIterator>
+	OutIterator andrew_convex_hull (InIterator begin, InIterator end, OutIterator out)
+	{
+		// need to store iterators to sort them
+		std::vector<InIterator> its(end - begin);
+		std::iota(its.begin(), its.end(), begin);
+
+		// sort lexicographically
+		std::sort(its.begin(), its.end(), [](auto i1, auto i2){ return *i1 < *i2; });
+
+		std::vector<InIterator> hull;
+		hull.push_back(its.front());
+
+		// find lower half of the hull
+		for (auto it = its.begin() + 1; it != its.end(); ++it)
+		{
+			while (hull.size() >= 2 && orientation(**(hull.end() - 2), **(hull.end() - 1), **it) != sign_t::positive)
+				hull.pop_back();
+
+			hull.push_back(*it);
+		}
+
+		// find upper half of the hull
+		for (auto it = its.rbegin() + 1; it != its.rend(); ++it)
+		{
+			if (*it == *(hull.end() - 2))
+				continue;
+
+			while (hull.size() >= 2 && orientation(**(hull.end() - 2), **(hull.end() - 1), **it) != sign_t::positive)
+				hull.pop_back();
+
+			hull.push_back(*it);
+		}
+
+		// the last upper part point is the first lower part point - remove it
+		hull.pop_back();
+
+		// copy result to output
+		return std::copy(hull.begin(), hull.end(), out);
+	}
+
+}

libs/cg/include/psemek/cg/convex_hull_2d/graham.hpp

@@ -0,0 +1,54 @@
+#pragma once
+
+#include <psemek/geom/orientation.hpp>
+
+#include <algorithm>
+#include <numeric>
+#include <vector>
+
+namespace cg
+{
+
+	template <typename InIterator, typename OutIterator>
+	OutIterator graham_convex_hull (InIterator begin, InIterator end, OutIterator out)
+	{
+		// need to store iterators to sort them
+		std::vector<InIterator> its(end - begin);
+		std::iota(its.begin(), its.end(), begin);
+
+		// find leftmost point & move its iterator to beginning
+		{
+			auto leftmost = std::min_element(its.begin(), its.end(), [](auto i1, auto i2){ return *i1 < *i2; });
+			std::iter_swap(leftmost, its.begin());
+		}
+
+		// sort with respect to angle to leftmost point
+		std::sort(its.begin() + 1, its.end(), [&](auto i1, auto i2){
+			auto o = orientation(*its.front(), *i1, *i2);
+
+			// carefully deal with parallel points
+			if (o == sign_t::positive)
+				return true;
+			else if (o == sign_t::negative)
+				return false;
+
+			return *i1 < *i2;
+		});
+
+		// perform gift wrapping
+
+		auto hull_end = its.begin() + 1;
+
+		for (auto it = its.begin() + 1; it != its.end(); ++it)
+		{
+			while (hull_end - its.begin() >= 2 && orientation(**(hull_end - 2), **(hull_end - 1), **it) != sign_t::positive)
+				--hull_end;
+
+			*hull_end++ = *it;
+		}
+
+		// copy result to output
+		return std::copy(its.begin(), hull_end, out);
+	}
+
+}

libs/cg/include/psemek/cg/convex_hull_2d/jarvis.hpp

@@ -0,0 +1,51 @@
+#pragma once
+
+#include <psemek/geom/orientation.hpp>
+
+namespace cg
+{
+
+	template <typename InIterator, typename OutIterator>
+	OutIterator jarvis_convex_hull (InIterator begin, InIterator end, OutIterator out)
+	{
+		auto first_hull_point = std::min_element(begin, end);
+		auto last_hull_point = first_hull_point;
+
+		*out++ = last_hull_point;
+
+		while (true)
+		{
+			auto it = begin;
+			for (; it != end; ++it)
+			{
+				if (it == last_hull_point) continue;
+
+				bool is_hull_edge = true;
+				for (auto jt = begin; jt != end; ++jt)
+				{
+					if (jt == it) continue;
+					if (jt == last_hull_point) continue;
+
+					if (orientation(*last_hull_point, *it, *jt) != cg::sign_t::positive)
+					{
+						is_hull_edge = false;
+						break;
+					}
+				}
+
+				if (is_hull_edge)
+				{
+					break;
+				}
+			}
+
+			if (it == first_hull_point) break;
+
+			*out++ = it;
+			last_hull_point = it;
+		}
+
+		return out;
+	}
+
+}

libs/cg/include/psemek/cg/convex_hull_2d/quickhull.hpp

@@ -0,0 +1,67 @@
+#pragma once
+
+#include <psemek/geom/orientation.hpp>
+
+#include <algorithm>
+#include <numeric>
+#include <vector>
+
+namespace cg
+{
+
+	namespace detail
+	{
+
+		template <typename InIterator, typename ItIterator, typename OutIterator>
+		OutIterator quick_hull_recursive_helper (InIterator p1, InIterator p2, ItIterator begin, ItIterator end, OutIterator out)
+		{
+			if (begin == end)
+				return *out++ = p1;
+
+			// find point in [begin, end) furthest from segment (p1,p2)
+			auto mid = *std::max_element(begin, end, [&](auto it1, auto it2){
+				// TODO: design a robust predicate
+				return orientation(*it2, (*it2) + ((*p2) - (*p1)), *it1) == sign_t::positive;
+			});
+
+			auto end1 = std::partition(begin, end, [&](auto it){ return orientation(*p1, *it, *mid) == sign_t::positive; });
+			auto end2 = std::partition(end1, end, [&](auto it){ return orientation(*mid, *it, *p2) == sign_t::positive; });
+
+			out = quick_hull_recursive_helper(p1, mid, begin, end1, out);
+			out = quick_hull_recursive_helper(mid, p2, end1, end2, out);
+
+			return out;
+		}
+
+	}
+
+	template <typename InIterator, typename OutIterator>
+	OutIterator quick_hull (InIterator begin, InIterator end, OutIterator out)
+	{
+		// need to store iterators to move sets of points around
+		std::vector<InIterator> its(end - begin);
+		std::iota(its.begin(), its.end(), begin);
+
+		// find leftmost point and move it to the beginning
+		{
+			auto p = std::min_element(its.begin(), its.end(), [](auto it1, auto it2){ return *it1 < *it2; });
+			std::iter_swap(its.begin(), p);
+		}
+
+		// find the next hull point in clockwise order and move it to second place
+		{
+			auto p = std::find_if(its.begin() + 1, its.end(), [&](auto it){
+				return std::all_of(its.begin() + 1, its.end(), [&](auto jt){
+					return it == jt || orientation(*its.front(), *it, *jt) == sign_t::negative;
+				});
+			});
+			std::iter_swap(its.begin() + 1, p);
+		}
+
+		// everything is set up, do the recursion
+		out = detail::quick_hull_recursive_helper(its[0], its[1], its.begin() + 2, its.end(), out);
+
+		return *out++ = *(its.begin() + 1);
+	}
+
+}

libs/cg/include/psemek/cg/dcel.hpp

@@ -0,0 +1,477 @@
+#pragma once
+
+#include <psemek/cg/utility.hpp>
+
+#include <cstddef>
+#include <vector>
+#include <cassert>
+
+namespace cg
+{
+
+	namespace detail
+	{
+
+		template <typename Point, typename Index>
+		struct point_rec : util::ebo_helper<Point>
+		{
+			Index edge;
+		};
+
+		template <typename Edge, typename Index>
+		struct edge_rec : util::ebo_helper<Edge>
+		{
+			Index origin;
+			Index next;
+			Index twin;
+			Index face;
+		};
+
+		template <typename Face, typename Index>
+		struct face_rec : util::ebo_helper<Face>
+		{
+			Index edge;
+		};
+
+		template <typename DCEL, typename Index, typename Tag>
+		struct handle_base
+		{
+			handle_base ( )
+				: owner_(nullptr)
+				, i_(static_cast<Index>(-1))
+			{ }
+
+			handle_base (DCEL * owner, Index i)
+				: owner_(owner)
+				, i_(i)
+			{ }
+
+			DCEL & owner ( ) const { return *owner_; }
+
+			Index index ( ) const { return i_; }
+
+			explicit operator bool ( ) const { return i_ != static_cast<Index>(-1); }
+
+			friend bool operator == (handle_base const & h1, handle_base const & h2)
+			{
+				return h1.i_ == h2.i_;
+			}
+
+			friend bool operator != (handle_base const & h1, handle_base const & h2)
+			{
+				return !(h1 == h2);
+			}
+
+			friend bool operator < (handle_base const & h1, handle_base const & h2)
+			{
+				return h1.i_ < h2.i_;
+			}
+
+			friend bool operator > (handle_base const & h1, handle_base const & h2)
+			{
+				return h2 < h1;
+			}
+
+			friend bool operator <= (handle_base const & h1, handle_base const & h2)
+			{
+				return !(h2 < h1);
+			}
+
+			friend bool operator >= (handle_base const & h1, handle_base const & h2)
+			{
+				return !(h1 < h2);
+			}
+
+		protected:
+			DCEL * owner_;
+			Index i_;
+		};
+
+		struct point_tag { };
+		struct edge_tag { };
+		struct face_tag { };
+
+	}
+
+	template <typename Point, typename Edge = util::empty, typename Face = util::empty, typename Index = std::size_t>
+	struct dcel
+	{
+		static constexpr Index null = static_cast<Index>(-1);
+
+		using point_rec = detail::point_rec<Point, Index>;
+		using edge_rec = detail::edge_rec<Edge, Index>;
+		using face_rec = detail::face_rec<Face, Index>;
+
+		std::vector<point_rec> points;
+		std::vector<edge_rec> edges;
+		std::vector<face_rec> faces;
+
+		// TODO: const handles
+
+		// TODO: iteration over handles
+
+		struct point_handle;
+		struct edge_handle;
+		struct face_handle;
+
+		struct point_handle : detail::handle_base<dcel, Index, detail::point_tag>
+		{
+			using detail::handle_base<dcel, Index, detail::point_tag>::handle_base;
+
+			Point & data ( ) const;
+
+			edge_handle edge ( ) const;
+
+			void edge (edge_handle h);
+
+		protected:
+			point_rec & get ( ) const
+			{
+				assert(this->owner_ != nullptr);
+				assert(this->i_ < this->owner_->points.size());
+				return this->owner_->points[this->i_];
+			}
+		};
+
+		struct edge_handle : detail::handle_base<dcel, Index, detail::edge_tag>
+		{
+			using detail::handle_base<dcel, Index, detail::edge_tag>::handle_base;
+
+			Edge & data ( ) const;
+
+			point_handle origin ( ) const;
+			edge_handle next ( ) const;
+			edge_handle twin ( ) const;
+			face_handle face ( ) const;
+
+			void origin (point_handle h);
+			void next (edge_handle h);
+			void twin (edge_handle h);
+			void face (face_handle h);
+
+		protected:
+			edge_rec & get ( ) const
+			{
+				assert(this->owner_ != nullptr);
+				assert(this->i_ < this->owner_->edges.size());
+				return this->owner_->edges[this->i_];
+			}
+		};
+
+		struct face_handle : detail::handle_base<dcel, Index, detail::face_tag>
+		{
+			using detail::handle_base<dcel, Index, detail::face_tag>::handle_base;
+
+			Face & data ( ) const;
+
+			edge_handle edge ( ) const;
+
+			void edge (edge_handle h);
+
+		protected:
+			face_rec & get ( ) const
+			{
+				assert(this->owner_ != nullptr);
+				assert(this->i_ < this->owner_->faces.size());
+				return this->owner_->faces[this->i_];
+			}
+		};
+
+		point_handle point (Index i)
+		{
+			return {this, i};
+		}
+
+		edge_handle edge (Index i)
+		{
+			return {this, i};
+		}
+
+		face_handle face (Index i)
+		{
+			return {this, i};
+		}
+
+		point_handle push_point (Point data = {})
+		{
+			auto i = static_cast<Index>(points.size());
+			points.push_back({{std::move(data)}, null});
+			return point(i);
+		}
+
+		edge_handle push_edge (Edge data = {})
+		{
+			auto i = static_cast<Index>(edges.size());
+			edges.push_back({{std::move(data)}, null, null, null, null});
+			return edge(i);
+		}
+
+		face_handle push_face (Face data = {})
+		{
+			auto i = static_cast<Index>(faces.size());
+			faces.push_back({{std::move(data)}, null});
+			return face(i);
+		}
+
+		point_handle insert_point (Index i, Point data = {})
+		{
+			points.insert(points.begin() + i, {{std::move(data)}, null});
+			for (auto & e : edges)
+			{
+				if (e.origin >= i)
+					++e.origin;
+			}
+			return point(i);
+		}
+
+		edge_handle insert_edge (Index i, Edge data = {})
+		{
+			edges.insert(edges.begin() + i, {{std::move(data)}, null, null, null, null});
+			for (auto & p : points)
+			{
+				if (p.edge >= i)
+					++p.edge;
+			}
+			for (auto & e : edges)
+			{
+				if (e.next >= i)
+					++e.next;
+				if (e.twin >= i)
+					++e.twin;
+			}
+			for (auto & f : faces)
+			{
+				if (f.edge >= i)
+					++f.edge;
+			}
+			return edge(i);
+		}
+
+		face_handle insert_face (Index i, Face data = {})
+		{
+			faces.insert(faces.begin() + i, {{std::move(data)}, null});
+			for (auto & e : edges)
+			{
+				if (e.face >= i)
+					++e.face;
+			}
+			return face(i);
+		}
+
+		void remove_point (point_handle h)
+		{
+			points.erase(points.begin() + h.index());
+			for (auto & e : edges)
+			{
+				if (e.origin > h.index())
+					--e.origin;
+			}
+		}
+
+		void remove_edge (edge_handle h)
+		{
+			edges.erase(edges.begin() + h.index());
+			for (auto & p : points)
+			{
+				if (p.edge > h.index())
+					--p.edge;
+			}
+			for (auto & e : edges)
+			{
+				if (e.next > h.index())
+					--e.next;
+				if (e.twin > h.index())
+					--e.twin;
+			}
+			for (auto & f : faces)
+			{
+				if (f.edge > h.index())
+					--f.edge;
+			}
+		}
+
+		void remove_face (face_handle h)
+		{
+			faces.erase(faces.begin() + h.index());
+			for (auto & e : edges)
+			{
+				if (e.face > h.index())
+					--e.face;
+			}
+		}
+	};
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	Point & dcel<Point, Edge, Face, Index>::point_handle::data() const
+	{
+		return get().data();
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	typename dcel<Point, Edge, Face, Index>::edge_handle dcel<Point, Edge, Face, Index>::point_handle::edge() const
+	{
+		return edge_handle{this->owner_, get().edge};
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	void dcel<Point, Edge, Face, Index>::point_handle::edge(edge_handle h)
+	{
+		get().edge = h.index();
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	Edge & dcel<Point, Edge, Face, Index>::edge_handle::data() const
+	{
+		return get().data();
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	typename dcel<Point, Edge, Face, Index>::point_handle dcel<Point, Edge, Face, Index>::edge_handle::origin() const
+	{
+		return point_handle{this->owner_, get().origin};
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	typename dcel<Point, Edge, Face, Index>::edge_handle dcel<Point, Edge, Face, Index>::edge_handle::next() const
+	{
+		return edge_handle{this->owner_, get().next};
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	typename dcel<Point, Edge, Face, Index>::edge_handle dcel<Point, Edge, Face, Index>::edge_handle::twin() const
+	{
+		return edge_handle{this->owner_, get().twin};
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	typename dcel<Point, Edge, Face, Index>::face_handle dcel<Point, Edge, Face, Index>::edge_handle::face() const
+	{
+		return face_handle{this->owner_, get().face};
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	void dcel<Point, Edge, Face, Index>::edge_handle::origin(point_handle h)
+	{
+		get().origin = h.index();
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	void dcel<Point, Edge, Face, Index>::edge_handle::next(edge_handle h)
+	{
+		get().next = h.index();
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	void dcel<Point, Edge, Face, Index>::edge_handle::twin(edge_handle h)
+	{
+		get().twin = h.index();
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	void dcel<Point, Edge, Face, Index>::edge_handle::face(face_handle h)
+	{
+		get().face = h.index();
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	Face & dcel<Point, Edge, Face, Index>::face_handle::data() const
+	{
+		return get().data();
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	typename dcel<Point, Edge, Face, Index>::edge_handle dcel<Point, Edge, Face, Index>::face_handle::edge() const
+	{
+		return edge_handle{this->owner_, get().edge};
+	}
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	void dcel<Point, Edge, Face, Index>::face_handle::edge(edge_handle h)
+	{
+		get().edge = h.index();
+	}
+
+	template <
+		typename Point2, typename Edge2 = util::empty, typename Face2 = util::empty, typename Index2 = std::size_t,
+		typename Point, typename Edge, typename Face, typename Index,
+		typename PointFn, typename EdgeFn, typename FaceFn>
+	dcel<Point2, Edge2, Face2, Index2> map(dcel<Point, Edge, Face, Index> const & d, PointFn && point_fn, EdgeFn && edge_fn, FaceFn && face_fn)
+	{
+		using result_type = dcel<Point2, Edge2, Face2, Index2>;
+		result_type result;
+
+		result.points.reserve(d.points.size());
+		result.edges.reserve(d.edges.size());
+		result.faces.reserve(d.faces.size());
+
+		for (auto const & p : d.points)
+		{
+			typename result_type::point_rec rec;
+			rec.edge = static_cast<Index2>(p.edge);
+			rec.data() = point_fn(p.data());
+			result.points.push_back(std::move(rec));
+		}
+
+		for (auto const & e : d.edges)
+		{
+			typename result_type::edge_rec rec;
+			rec.origin = static_cast<Index2>(e.origin);
+			rec.next = static_cast<Index2>(e.next);
+			rec.twin = static_cast<Index2>(e.twin);
+			rec.face = static_cast<Index2>(e.face);
+			rec.data() = edge_fn(e.data());
+			result.edges.push_back(std::move(rec));
+		}
+
+		for (auto const & f : d.faces)
+		{
+			typename result_type::face_rec rec;
+			rec.edge = static_cast<Index2>(f.edge);
+			rec.data() = face_fn(f.data());
+			result.faces.push_back(std::move(rec));
+		}
+
+		return result;
+	}
+
+	template <typename Index = std::size_t, typename Iterator>
+	auto polygon_dcel (Iterator begin, Iterator end)
+	{
+		using point_type = std::decay_t<decltype(*begin)>;
+
+		dcel<point_type, util::empty, util::empty, Index> result;
+
+		auto outer_face = result.push_face();
+		auto inner_face = result.push_face();
+
+		outer_face.edge(result.edge(0));
+		inner_face.edge(result.edge(1));
+
+		for (; begin != end; ++begin)
+		{
+			result.push_point(*begin);
+			result.push_edge();
+			result.push_edge();
+		}
+
+		std::size_t const N = result.points.size();
+
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result.points[i].edge = 2 * i;
+
+			result.edges[2 * i].face = 0;
+			result.edges[2 * i].next = (2 * i + 2) % (2 * N);
+			result.edges[2 * i].twin = 2 * i + 1;
+			result.edges[2 * i].origin = i;
+
+			result.edges[2 * i + 1].face = 1;
+			result.edges[2 * i + 1].next = (2 * i + 2 * N - 1) % (2 * N);
+			result.edges[2 * i + 1].twin = 2 * i;
+			result.edges[2 * i + 1].origin = (i + 1) % N;
+		}
+
+		return result;
+	}
+
+}

libs/cg/include/psemek/cg/dual.hpp

@@ -0,0 +1,50 @@
+#pragma once
+
+#include <psemek/cg/dcel.hpp>
+
+namespace cg
+{
+
+	template <typename Point, typename Edge, typename Face, typename Index>
+	dcel<Face, Edge, Point, Index> dual (dcel<Point, Edge, Face, Index> const & d)
+	{
+		dcel<Face, Edge, Point, Index> result;
+
+		result.edges.reserve(d.edges.size());
+		for (auto const & e : d.edges)
+		{
+			result.push_edge(e.data());
+			auto & n = result.edges.back();
+			n.twin = e.twin;
+			n.face = e.origin;
+			n.origin = d.edges[e.twin].face;
+		}
+
+		for (Index e = 0; e < result.edges.size(); ++e)
+		{
+			auto next = [&](auto i){ return d.edges[d.edges[i].twin].next; };
+
+			Index i = e, n = next(i);
+			for (; n != e; i = n, n = next(n));
+
+			result.edges[e].next = i;
+		}
+
+		result.points.reserve(d.faces.size());
+		for (auto const & f : d.faces)
+		{
+			result.push_point(f.data());
+			result.points.back().edge = d.edges[f.edge].twin;
+		}
+
+		result.faces.reserve(d.points.size());
+		for (auto const & p : d.points)
+		{
+			result.push_face(p.data());
+			result.faces.back().edge = p.edge;
+		}
+
+		return result;
+	}
+
+}

libs/cg/include/psemek/cg/segment_intersection.hpp

@@ -0,0 +1,147 @@
+#pragma once
+
+#include <psemek/geom/intersection.hpp>
+#include <psemek/cg/utility.hpp>
+
+#include <cstddef>
+#include <algorithm>
+#include <variant>
+#include <numeric>
+#include <queue>
+#include <deque>
+#include <set>
+
+namespace cg
+{
+
+	template <typename InIterator, typename OutIterator>
+	OutIterator segment_intersection (InIterator begin, InIterator end, OutIterator out)
+	{
+		using segment = std::decay_t<decltype(*begin)>;
+		using point = typename segment::point_type;
+		using scalar = typename point::scalar_type;
+
+		enum event_type
+		{
+			BEGIN,
+			END,
+			CROSS,
+		};
+
+		// TODO: robust event comparison predicate
+
+		struct event
+		{
+			InIterator it1, it2;
+			point p;
+			event_type type;
+		};
+
+		struct event_comparator
+		{
+			bool operator() (event const & e1, event const & e2) const
+			{
+				return e1.p > e2.p;
+			}
+		};
+
+		std::priority_queue<event, std::vector<event>, event_comparator> queue;
+
+		for (auto it = begin; it != end; ++it)
+		{
+			event e_begin{it, end, (*it)[0], BEGIN};
+			event e_end  {it, end, (*it)[1], END  };
+
+			if (e_begin.p > e_end.p) std::swap(e_begin.p, e_end.p);
+			queue.push(e_begin);
+			queue.push(e_end);
+		}
+
+		auto push_cross_event = [&queue](InIterator it, InIterator jt){
+			auto ion = intersection(*it, *jt);
+			if (auto p = std::get_if<point>(&ion))
+			{
+				queue.push({it, jt, *p, CROSS});
+			}
+		};
+
+		scalar sweep_line = -std::numeric_limits<scalar>::infinity();
+
+		struct status_comparator
+		{
+			scalar & sweep_line;
+
+			bool operator() (InIterator it1, InIterator it2) const
+			{
+				scalar const eps = 1e-4;
+				scalar const y1 = (*it1)[0][1] + ((*it1)[1][1] -(*it1)[0][1]) * (sweep_line + eps - (*it1)[0][0]) / ((*it1)[1][0] - (*it1)[0][0]);
+				scalar const y2 = (*it2)[0][1] + ((*it2)[1][1] -(*it2)[0][1]) * (sweep_line + eps - (*it2)[0][0]) / ((*it2)[1][0] - (*it2)[0][0]);
+
+				return y1 < y2;
+			}
+		};
+
+		std::set<InIterator, status_comparator> status{status_comparator{sweep_line}};
+
+		while (!queue.empty())
+		{
+			event e = queue.top();
+			queue.pop();
+
+			switch (e.type)
+			{
+			case BEGIN: {
+				auto it = status.insert(e.it1).first;
+				if (it != status.begin()) push_cross_event(*std::prev(it), *it);
+				if (std::next(it) != status.end()) push_cross_event(*it, *std::next(it));
+				sweep_line = e.p[0];
+				break;
+			}
+			case END: {
+				auto it = status.find(e.it1);
+				if (it == status.end())
+				{
+					int fuck = 42;
+				}
+				else
+				{
+					if (it != status.begin() && std::next(it) != status.end()) push_cross_event(*std::prev(it), *std::next(it));
+					status.erase(it);
+				}
+				sweep_line = e.p[0];
+				break;
+			}
+			case CROSS: {
+				*out++ = {e.it1, e.it2};
+
+				status.erase(e.it1);
+				status.erase(e.it2);
+				sweep_line = e.p[0];
+				auto it1 = status.insert(e.it1).first;
+				auto it2 = status.insert(e.it2).first;
+				if (it1 == status.end())
+				{
+					int fuck_really = 42;
+				}
+				if (it2 == status.end())
+				{
+					int fuck_really = 42;
+				}
+
+				if (std::next(it2) != it1)
+				{
+					bool check = std::next(it1) == it2;
+
+					int fuck = 42;
+				}
+				if (it2 != status.begin()) push_cross_event(*std::prev(it2), *it2);
+				if (std::next(it1) != status.end()) push_cross_event(*it1, *std::next(it1));
+				break;
+			}
+			}
+		}
+
+		return out;
+	}
+
+}

libs/cg/include/psemek/cg/triangulation/delaunay.hpp

@@ -0,0 +1,103 @@
+#pragma once
+
+#include <psemek/geom/incircle.hpp>
+#include <psemek/cg/triangulation/triangulation.hpp>
+
+#include <queue>
+
+namespace cg
+{
+
+	template <typename Index = std::size_t, typename InputIterator>
+	auto delaunay (InputIterator begin, InputIterator end)
+	{
+		std::vector<Index> edge_queue;
+
+		auto callback = [&](auto & dcel, auto p)
+		{
+			auto outer_face = dcel.face(0);
+
+			// grab interior edges outgoing from p
+			for (auto e = p.edge();;)
+			{
+				auto t = e.twin();
+
+				if (e.face() != outer_face)
+				{
+					if (t.face() != outer_face)
+					{
+						edge_queue.push_back(e.index());
+					}
+
+					auto n = e.next();
+					if (n.twin().face() != outer_face)
+					{
+						edge_queue.push_back(n.index());
+					}
+				}
+
+				e = t.next();
+
+				if (e == p.edge()) break;
+			}
+
+			while (!edge_queue.empty())
+			{
+				auto e = dcel.edge(edge_queue.back());
+				edge_queue.pop_back();
+
+				auto next = e.next();
+				auto prev = next.next();
+				auto twin = e.twin();
+				auto tnext = twin.next();
+				auto tprev = tnext.next();
+
+				auto p0 = prev.origin();
+				auto p1 = e.origin();
+				auto p2 = next.origin();
+				auto p3 = e.twin().next().next().origin();
+
+				// decide if a flip is needed
+				if (in_circle(*p0.data(), *p1.data(), *p2.data(), *p3.data()) != sign_t::positive) continue;
+
+				auto f0 = e.face();
+				auto f1 = twin.face();
+
+				e.origin(p0);
+				e.next(tprev);
+
+				next.next(e);
+
+				prev.face(f1);
+				prev.next(tnext);
+
+				twin.origin(p3);
+				twin.next(prev);
+
+				tnext.next(twin);
+
+				tprev.face(f0);
+				tprev.next(next);
+
+				p1.edge(tnext);
+				p2.edge(next);
+
+				f0.edge(e);
+				f1.edge(twin);
+
+				auto push = [&](auto e)
+				{
+					if (e.twin().face() != outer_face)
+						edge_queue.push_back(e.index());
+				};
+
+				push(next);
+				push(prev);
+				push(tnext);
+				push(tprev);
+			}
+		};
+		return detail::triangulate(begin, end, callback);
+	}
+
+}

libs/cg/include/psemek/cg/triangulation/triangulation.hpp

@@ -0,0 +1,262 @@
+#pragma once
+
+#include <psemek/geom/orientation.hpp>
+#include <psemek/cg/dcel.hpp>
+#include <psemek/cg/utility.hpp>
+
+#include <algorithm>
+#include <numeric>
+
+namespace cg
+{
+
+	namespace detail
+	{
+
+		template <typename Index = std::size_t, typename Iterator, typename Callback>
+		auto triangulate (Iterator begin, Iterator end, Callback && callback)
+		{
+			using point_type = std::decay_t<decltype(*begin)>;
+			static_assert(point_type::dimension == 2);
+
+			using dcel_type = dcel<Iterator, util::empty, util::empty, Index>;
+
+			dcel_type result;
+
+			std::vector<Iterator> sweepline_events;
+			sweepline_events.resize(std::distance(begin, end));
+			std::iota(sweepline_events.begin(), sweepline_events.end(), begin);
+			std::sort(sweepline_events.begin(), sweepline_events.end(), [](auto it, auto jt){ return *it < *jt; });
+
+			auto it = sweepline_events.begin();
+
+			auto outer_face = result.push_face();
+
+			// handle degenerate cases
+			if (sweepline_events.empty()) return result;
+
+			if (sweepline_events.size() == 1)
+			{
+				result.push_point(*it++);
+				return result;
+			}
+
+			{
+				auto const N = sweepline_events.size();
+				result.points.reserve(N);
+				result.edges.reserve(3*(N-1));
+				result.faces.reserve(2*(N-1));
+			}
+
+			typename dcel_type::edge_handle hull_start;
+
+			{
+				auto p0 = result.push_point(*it++);
+				auto p1 = result.push_point(*it++);
+
+				auto e01 = result.push_edge();
+				auto e10 = result.push_edge();
+
+				p0.edge(e01);
+				p1.edge(e10);
+
+				e01.origin(p0);
+				e01.next(e10);
+				e01.face(outer_face);
+				e01.twin(e10);
+				e10.origin(p1);
+				e10.next(e01);
+				e10.face(outer_face);
+				e10.twin(e01);
+				outer_face.edge(e01);
+
+				hull_start = e01;
+			}
+
+			for (; it != sweepline_events.end(); ++it)
+			{
+				// walk along the hull clockwise
+
+				auto cur_hull_edge = hull_start;
+				auto next_hull_edge = cur_hull_edge.next();
+
+				auto hp0 = cur_hull_edge.origin();
+				auto hp1 = next_hull_edge.origin();
+
+				auto move_hull_edge = [&]{
+					cur_hull_edge = next_hull_edge;
+					next_hull_edge = next_hull_edge.next();
+					hp0 = hp1;
+					hp1 = next_hull_edge.origin();
+				};
+
+				bool degenerate = false;
+
+				// find first hull edge visible from new point
+				while (orientation(**it, *hp0.data(), *hp1.data()) != sign_t::positive)
+				{
+					move_hull_edge();
+					if (cur_hull_edge == hull_start)
+					{
+						degenerate = true;
+						break;
+					}
+				}
+
+				if (degenerate)
+				{
+					// the whole hull is just a sequence of points on a line
+					// find the closest & connect via a single edge
+
+					auto q = result.point(0);
+
+					bool vertical = false;
+
+					// find rightmost point
+					for (Index i = 1; i < result.points.size(); ++i)
+					{
+						auto const x = (*result.points[i].data())[0];
+						auto const qx = (*q.data())[0];
+
+						if (x == qx)
+						{
+							vertical = true;
+							break;
+						}
+
+						if (x > qx)
+						{
+							q = result.point(i);
+						}
+					}
+
+					// all points lie on a vertical line - find topmost
+					if (vertical)
+					{
+						q = result.point(0);
+
+						for (Index i = 1; i < result.points.size(); ++i)
+						{
+							auto const y = (*result.points[i].data())[1];
+							auto const qy = (*q.data())[1];
+
+							if (y > qy)
+							{
+								q = result.point(i);
+							}
+						}
+					}
+
+					auto p = result.push_point(*it);
+
+					auto qnext = q.edge();
+					auto qprev = qnext.twin();
+
+					auto pnext = result.push_edge();
+					auto pprev = result.push_edge();
+
+					qprev.next(pprev);
+
+					pprev.next(pnext);
+					pprev.origin(q);
+					pprev.twin(pnext);
+					pprev.face(outer_face);
+
+					pnext.next(qnext);
+					pnext.origin(p);
+					pnext.twin(pprev);
+					pnext.face(outer_face);
+
+					p.edge(pnext);
+
+					continue;
+				}
+
+				auto p = result.push_point(*it);
+
+				auto mid_edge = result.push_edge();
+				mid_edge.origin(hp0);
+				mid_edge.face(outer_face);
+
+				auto first_mid_edge = mid_edge;
+
+				{
+					// fix prev_hull_edge.next
+					auto prev_hull_edge = cur_hull_edge.twin();
+
+					while (prev_hull_edge.face() != outer_face)
+					{
+						prev_hull_edge = prev_hull_edge.next().twin();
+					}
+					prev_hull_edge.next(first_mid_edge);
+				}
+
+				// until current edge is visible
+				while (orientation(**it, *hp0.data(), *hp1.data()) == sign_t::positive)
+				{
+					// fill new triangle
+
+					auto ep0 = result.push_edge();
+					auto e1p = result.push_edge();
+
+					auto f = result.push_face();
+
+					if (!p.edge())
+					{
+						p.edge(ep0);
+					}
+
+					mid_edge.twin(ep0);
+
+					ep0.origin(p);
+					ep0.next(cur_hull_edge);
+					ep0.twin(mid_edge);
+					ep0.face(f);
+
+					e1p.origin(hp1);
+					e1p.next(ep0);
+					e1p.face(f);
+
+					cur_hull_edge.next(e1p);
+					cur_hull_edge.face(f);
+
+					f.edge(ep0);
+
+					mid_edge = e1p;
+
+					move_hull_edge();
+				}
+
+				auto last_mid_edge = result.push_edge();
+				last_mid_edge.origin(p);
+				last_mid_edge.face(outer_face);
+				last_mid_edge.twin(mid_edge);
+				last_mid_edge.next(cur_hull_edge);
+
+				mid_edge.twin(last_mid_edge);
+
+				first_mid_edge.next(last_mid_edge);
+
+				// update first hull edge
+
+				if (hull_start.face() != outer_face)
+				{
+					hull_start = hull_start.next().next().twin();
+					outer_face.edge(hull_start);
+				}
+
+				callback(result, p);
+			}
+
+			return result;
+		}
+
+	}
+
+	template <typename Index = std::size_t, typename InputIterator>
+	auto triangulate (InputIterator begin, InputIterator end)
+	{
+		return detail::triangulate(begin, end, util::nop);
+	}
+
+}

libs/cg/include/psemek/cg/utility.hpp

@@ -0,0 +1,68 @@
+#pragma once
+
+#include <utility>
+#include <variant>
+
+namespace cg::util
+{
+
+	struct empty{};
+
+	template <typename ... Fs>
+	struct overload_impl
+		: Fs ...
+	{
+		using Fs::operator() ...;
+	};
+
+	template <typename ... Fs>
+	auto overload (Fs && ... fs)
+	{
+		return overload_impl<Fs...>{ std::forward<Fs>(fs) ... };
+	}
+
+	template <typename Variant, typename ... Fs>
+	decltype(auto) match (Variant && v, Fs && ... fs)
+	{
+		return std::visit(overload(std::forward<Fs>(fs)...), std::forward<Variant>(v));
+	}
+
+	constexpr auto id = [](auto const & x){ return x; };
+
+	constexpr auto nop = [](auto const & ...){};
+
+	namespace detail
+	{
+
+		template <typename T, bool CanInherit>
+		struct ebo_helper_impl : T
+		{
+			T & data() { return *this; }
+			T const & data() const { return *this; }
+
+			template <typename ... Args>
+			ebo_helper_impl (Args && ... args)
+				: T(std::forward<Args>(args)...)
+			{ }
+		};
+
+		template <typename T>
+		struct ebo_helper_impl<T, false>
+		{
+			T value;
+
+			T & data() { return value; }
+			T const & data() const { return value; }
+
+			template <typename ... Args>
+			ebo_helper_impl (Args && ... args)
+				: value(std::forward<Args>(args)...)
+			{ }
+		};
+
+	}
+
+	template <typename T>
+	using ebo_helper = detail::ebo_helper_impl<T, std::is_class_v<T> && !(std::is_final_v<T>)>;
+
+}

libs/cg/include/psemek/cg/voronoi.hpp

@@ -0,0 +1,194 @@
+#pragma once
+
+#include <psemek/cg/triangulation/delaunay.hpp>
+#include <psemek/cg/dual.hpp>
+#include <psemek/geom/gauss.hpp>
+#include <psemek/cg/bbox.hpp>
+
+namespace cg
+{
+
+	// Point #0 of the resulting dcel corresponds to the outer face,
+	// i.e. the point at infinity, and has unspecified coordinates.
+	template <typename Index = std::size_t, typename Iterator>
+	auto voronoi (Iterator begin, Iterator end)
+	{
+		using point_type = std::decay_t<decltype(*begin)>;
+
+		auto del = dual(delaunay<Index>(begin, end));
+
+		dcel<point_type, util::empty, Iterator, Index> result;
+
+		result.points.reserve(del.points.size());
+		for (std::size_t i = 0; i < del.points.size(); ++i)
+		{
+			auto p = del.point(i);
+
+			point_type q[3];
+
+			auto e = p.edge();
+			q[0] = *e.face().data();
+			e = e.twin().next();
+			q[1] = *e.face().data();
+			e = e.twin().next();
+			q[2] = *e.face().data();
+
+			using T = typename point_type::scalar_type;
+
+			matrix<T, 3, 3> m;
+			vector<T, 3> b;
+
+			for (std::size_t i = 0; i < 3; ++i)
+			{
+				m[i][0] = q[i][0];
+				m[i][1] = q[i][1];
+				m[i][2] = T{1};
+
+				b[i] = q[i][0] * q[i][0] + q[i][1] * q[i][1];
+			}
+
+			gauss(m, b);
+
+			auto newp = result.push_point(point_type{b[0] / 2, b[1] / 2});
+			newp.edge(result.edge(p.edge().index()));
+		}
+
+		result.edges = std::move(del.edges);
+		result.faces = std::move(del.faces);
+
+		return result;
+	}
+
+	// Replaces the #0 infinite point of a dcel with a proper outer face
+	template <typename Point, typename Face, typename Edge, typename Index>
+	void remove_infinite_point (dcel<Point, Face, Edge, Index> & dcel)
+	{
+		// insert outer face at position #0
+		auto outer_face = dcel.insert_face(0);
+
+		std::vector<std::size_t> delete_edges;
+		std::vector<std::size_t> delete_faces;
+
+		auto infinite_point = dcel.point(0);
+
+		// iterate edges of infinite vertex #0 and fix hull
+		for (auto e = infinite_point.edge();;)
+		{
+			for (auto n = e.next(); n.twin().origin() != infinite_point; n = n.next())
+			{
+				n.face(outer_face);
+				n.origin().edge(n);
+				if (!outer_face.edge())
+					outer_face.edge(n);
+			}
+
+			auto n = e;
+			for (; n.next() != e.twin(); n = n.next().twin());
+
+			if (n.origin() != infinite_point)
+			{
+				auto t = e.next();
+				for (; t.twin().origin() == infinite_point; t = t.twin().next());
+
+				n.next(t);
+			}
+
+			delete_faces.push_back(e.face().index());
+			delete_edges.push_back(e.index());
+			e = e.twin();
+			delete_edges.push_back(e.index());
+			e = e.next();
+
+			if (e == infinite_point.edge()) break;
+		}
+
+		// remove edges incident to infinite vertex #0
+		std::sort(delete_edges.begin(), delete_edges.end(), std::greater<>{});
+		for (auto e : delete_edges)
+			dcel.remove_edge(dcel.edge(e));
+
+		// remove vertex #0
+		dcel.remove_point(infinite_point);
+
+		std::sort(delete_faces.begin(), delete_faces.end(), std::greater<>{});
+		for (auto f : delete_faces)
+			dcel.remove_face(dcel.face(f));
+	}
+
+	// iteratively removes degenerate edges
+	template <typename Point, typename Face, typename Edge, typename Index>
+	void remove_degenerate_edges (dcel<Point, Face, Edge, Index> & dcel)
+	{
+		std::vector<std::size_t> delete_points;
+		std::vector<std::size_t> delete_edges;
+
+		for (std::size_t i = 0; i < dcel.points.size(); ++i)
+		{
+			auto p = dcel.point(i);
+			while (p.edge().twin().next() == p.edge())
+			{
+				delete_points.push_back(p.index());
+				delete_edges.push_back(p.edge().index());
+				delete_edges.push_back(p.edge().twin().index());
+
+				p.edge().next().origin().edge(p.edge().next());
+				p.edge().face().edge(p.edge().next());
+
+				auto n = p.edge();
+				while (true)
+				{
+					auto m = n.next().twin();
+					if (m == p.edge()) break;
+					n = m;
+				}
+
+				n.next(p.edge().next());
+
+				p = p.edge().next().origin();
+			}
+		}
+
+		std::sort(delete_points.begin(), delete_points.end(), std::greater<>{});
+		for (auto p : delete_points)
+			dcel.remove_point(dcel.point(p));
+
+		std::sort(delete_edges.begin(), delete_edges.end(), std::greater<>{});
+		for (auto e : delete_edges)
+			dcel.remove_edge(dcel.edge(e));
+	}
+
+	// Outputs 4 points such that the finite faces of the voronoi diagram of (anything
+	// inside the convex hull of input + the output points) covers the convex hull of input
+	template <typename InputIterator, typename OutputIterator>
+	OutputIterator bounded_voronoi_extra_points (InputIterator begin, InputIterator end, OutputIterator out)
+	{
+		using point_type = std::decay_t<decltype(*begin)>;
+
+		static_assert(point_type::dimension == 2);
+
+		auto box = bbox(begin, end);
+
+		point_type const center {(box[0].min + box[0].max) / 2, (box[1].min + box[1].max) / 2};
+
+		auto const R = std::sqrt(box[0].length() * box[0].length() + box[1].length() * box[1].length());
+
+		*out++ = point_type{center[0] + R, center[1]};
+		*out++ = point_type{center[0], center[1] + R};
+		*out++ = point_type{center[0] - R, center[1]};
+		*out++ = point_type{center[0], center[1] - R};
+
+		return out;
+	}
+
+	// Compute a clipped voronoi tessellation of (begin, end - 4)
+	// Last 4 points will be used for temporary data
+	template <typename InputIterator>
+	auto bounded_voronoi (InputIterator begin, InputIterator end)
+	{
+		bounded_voronoi_extra_points(begin, end - 4, end - 4);
+		auto dcel = cg::voronoi(begin, end);
+		cg::remove_infinite_point(dcel);
+		return dcel;
+	}
+
+}

libs/geom/CMakeLists.txt

@@ -0,0 +1,9 @@
+find_package(Boost REQUIRED)
+find_package(GMP REQUIRED)
+
+file(GLOB_RECURSE PSEMEK_GEOM_HEADERS "include/*.hpp")
+file(GLOB_RECURSE PSEMEK_GEOM_SOURCES "source/*.cpp")
+
+add_library(geom STATIC ${PSEMEK_GEOM_HEADERS} ${PSEMEK_GEOM_SOURCES})
+target_include_directories(geom PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
+target_link_libraries(geom PUBLIC Boost::boost GMP)

libs/geom/include/psemek/geom/bezier.hpp

@@ -0,0 +1,41 @@
+#pragma once
+
+#include <psemek/geom/point.hpp>
+
+#include <vector>
+#include <stdexcept>
+
+namespace psemek::geom
+{
+
+	template <typename Point>
+	struct bezier
+	{
+		bezier (std::vector<Point> points)
+			: points_(std::move(points))
+		{
+			if (points_.empty())
+				throw std::runtime_error("Points array should be non-empty");
+			temp_.resize(points_.size());
+		}
+
+		template <typename T>
+		auto operator() (T const & t) const
+		{
+			// In-place de Casteljau's algorithm
+			temp_ = points_;
+
+			for (std::size_t k = 0; k + 1 < points_.size(); ++k)
+			{
+				for (std::size_t i = 0; i + k + 1 < points_.size(); ++i)
+					temp_[i] = lerp(temp_[i], temp_[i + 1], t);
+			}
+			return temp_.front();
+		}
+
+	private:
+		std::vector<Point> const points_;
+		std::vector<Point> mutable temp_;
+	};
+
+}

libs/geom/include/psemek/geom/box.hpp

@@ -0,0 +1,208 @@
+#pragma once
+
+#include <psemek/geom/interval.hpp>
+#include <psemek/geom/point.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t N>
+	struct box
+	{
+		interval<T> axes[N];
+
+		using point_type = point<T, N>;
+		using vector_type = vector<T, N>;
+
+		interval<T> & operator [] (std::size_t i)
+		{
+			return axes[i];
+		}
+
+		interval<T> const & operator [] (std::size_t i) const
+		{
+			return axes[i];
+		}
+
+		// singleton box
+		static box singleton (point_type const & p)
+		{
+			box b;
+			for (std::size_t i = 0; i < N; ++i)
+				b.axes[i] = interval<T>::singleton(p[i]);
+			return b;
+		}
+
+		bool empty ( ) const
+		{
+			for (auto const & i : axes)
+				if (i.empty())
+					return true;
+			return false;
+		}
+
+		T size ( ) const
+		{
+			T result = T{1};
+			for (auto const & i : axes)
+				result *= i.size();
+			return result;
+		}
+
+		box & operator += (vector_type const & delta);
+		box & operator -= (vector_type const & delta);
+
+		box & operator &= (point_type const & p);
+		box & operator |= (point_type const & p);
+
+		box & operator &= (box const & b);
+		box & operator |= (box const & b);
+	};
+
+	template <typename T, std::size_t N>
+	box<T, N> operator + (box<T, N> const & b, vector<T, N> const & delta)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = b[i] + delta[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> operator + (vector<T, N> const & delta, box<T, N> const & b)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = delta[i] + b[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> operator - (box<T, N> const & b, vector<T, N> const & delta)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = b[i] - delta[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> operator & (box<T, N> const & b, point<T, N> const & p)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = b[i] & p[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> operator & (point<T, N> const & p, box<T, N> const & b)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = p[i] & b[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> operator | (box<T, N> const & b, point<T, N> const & p)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = b[i] | p[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> operator | (point<T, N> const & p, box<T, N> const & b)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = p[i] | b[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> operator & (box<T, N> const & b1, box<T, N> const & b2)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = b1[i] & b2[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> operator | (box<T, N> const & b1, box<T, N> const & b2)
+	{
+		box<T, N> result;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			result[i] = b1[i] | b2[i];
+		}
+		return result;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> & box<T, N>::operator += (vector<T, N> const & delta)
+	{
+		return *this = *this + delta;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> & box<T, N>::operator -= (vector<T, N> const & delta)
+	{
+		return *this = *this - delta;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> & box<T, N>::operator &= (point<T, N> const & p)
+	{
+		return *this = *this & p;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> & box<T, N>::operator |= (point<T, N> const & p)
+	{
+		return *this = *this | p;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> & box<T, N>::operator &= (box<T, N> const & b)
+	{
+		return *this = *this & b;
+	}
+
+	template <typename T, std::size_t N>
+	box<T, N> & box<T, N>::operator |= (box<T, N> const & b)
+	{
+		return *this = *this | b;
+	}
+
+	template <typename Stream, typename T, std::size_t N>
+	Stream & operator << (Stream & os, box<T, N> const & b)
+	{
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			if (i != 0) os << 'x';
+			os << b[i];
+		}
+		return os;
+	}
+
+}

libs/geom/include/psemek/geom/camera.hpp

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/matrix.hpp>
+
+namespace psemek::geom
+{
+
+	struct camera
+	{
+		virtual matrix<float, 4, 4> projection() const = 0;
+		virtual matrix<float, 4, 4> view() const = 0;
+
+		virtual matrix<float, 4, 4> transform() const;
+		virtual point<float, 3> position() const;
+
+		// View direction is -axis_z(), or equivalently direction(0,0)
+		// NB: these vectors are not necessarily normalized
+		virtual vector<float, 3> axis_x() const;
+		virtual vector<float, 3> axis_y() const;
+		virtual vector<float, 3> axis_z() const;
+
+		// x, y are in clip-space [-1, 1]
+		// NB: this vector is not necessarily normalized
+		virtual vector<float, 3> direction(float x, float y) const;
+
+		// view frustum is the intersection of half-spaces defined by dot(clip_plane, x) >= 0
+		virtual std::array<vector<float, 4>, 6> clip_planes() const;
+
+		virtual ~camera() = default;
+	};
+
+	struct perspective_camera
+		: camera
+	{
+		float fov_x;
+		float fov_y;
+
+		float near_clip;
+		float far_clip;
+
+		matrix<float, 4, 4> projection() const override;
+	};
+
+	struct spherical_camera
+		: perspective_camera
+	{
+		// assumes up is +Y
+		//
+		float azimuthal_angle;
+		float elevation_angle;
+		float distance;
+		point<float, 3> target;
+
+		matrix<float, 4, 4> view() const override;
+	};
+
+}

libs/geom/include/psemek/geom/contains.hpp

@@ -0,0 +1,37 @@
+#pragma once
+
+#include <psemek/geom/interval.hpp>
+#include <psemek/geom/box.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/simplex.hpp>
+#include <psemek/geom/orientation.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T>
+	bool contains (interval<T> const & i, T const & x)
+	{
+		return (i.min <= x) && (x <= i.max);
+	}
+
+	template <typename T, std::size_t N>
+	bool contains (box<T, N> const & b, point<T, N> const & p)
+	{
+		for (std::size_t i = 0; i < N; ++i)
+			if (!contains(b[i], p[i]))
+				return false;
+		return true;
+	}
+
+	template <typename T>
+	bool contains (triangle<point<T, 2>> const & t, point<T, 2> const & p)
+	{
+		return true
+			&& orientation(t[0], t[1], p) != sign_t::negative
+			&& orientation(t[1], t[2], p) != sign_t::negative
+			&& orientation(t[2], t[0], p) != sign_t::negative
+			;
+	}
+
+}

libs/geom/include/psemek/geom/gauss.hpp

@@ -0,0 +1,238 @@
+#pragma once
+
+#include <psemek/geom/matrix.hpp>
+
+#include <cmath>
+
+namespace psemek::geom
+{
+
+	namespace detail
+	{
+
+		// A helper to treat vectors & matrices uniformly
+
+		template <typename X>
+		struct gauss_helper;
+
+		template <typename X>
+		gauss_helper (X &) -> gauss_helper<X>;
+
+		template <typename T, std::size_t N>
+		struct gauss_helper<vector<T, N>>
+		{
+			vector<T, N> & v;
+
+			static constexpr std::size_t columns ( ) { return 1; }
+			T & at (std::size_t row, std::size_t) { return v[row]; }
+		};
+
+		template <typename T, std::size_t N, std::size_t M>
+		struct gauss_helper<matrix<T, N, M>>
+		{
+			matrix<T, N, M> & m;
+
+			static constexpr std::size_t columns ( ) { return M; }
+			T & at (std::size_t row, std::size_t col) { return m[row][col]; }
+		};
+
+		template <typename F, typename ... Args>
+		void for_each (F && f, Args & ... args)
+		{
+			(f(args), ...);
+		}
+
+	}
+
+	template <typename T, std::size_t N>
+	T det (matrix<T, N, N> m)
+	{
+		using std::abs;
+
+		std::size_t sign = 0;
+
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			// find maximal modulus along i-th column
+			auto M = abs(m[i][i]);
+			std::size_t k = i;
+			for (std::size_t j = i + 1; j < N; ++j)
+			{
+				auto const n = abs(m[j][i]);
+				if (n > M)
+				{
+					M = n;
+					k = j;
+				}
+			}
+
+			// avoid dividing by zero
+			if (M == T{})
+				return T{};
+
+			// swap rows i & k
+			if (i != k)
+			{
+				sign += 1;
+
+				for (std::size_t j = i; j < N; ++j)
+					std::swap(m[i][j], m[k][j]);
+			}
+
+			// zero out the column under i,i
+			for (std::size_t k = i + 1; k < N; ++k)
+			{
+				auto r = m[k][i] / m[i][i];
+				for (std::size_t j = i; j < N; ++j)
+				{
+					m[k][j] -= m[i][j] * r;
+				}
+			}
+		}
+
+		// determinant is product of diagonal entries + sign
+		T res = ((sign % 2) == 0) ? T{1} : -T{1};
+		for (std::size_t i = 0; i < N; ++i)
+			res *= m[i][i];
+		return res;
+	}
+
+	// each element of the rhs can be a vector<N> or a matrix<N, K>
+	// the set of equations (m * x = rhs)... is solved simultaneously column-wise
+	// returns false is the matrix is degenerate
+	template <typename T, std::size_t N, typename ... RHS>
+	bool gauss (matrix<T, N, N> m, RHS & ... rhs)
+	{
+		using std::abs;
+
+		// forward elimination
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			// find maximal modulus along i-th column
+			auto M = abs(m[i][i]);
+			std::size_t k = i;
+			for (std::size_t j = i + 1; j < N; ++j)
+			{
+				auto const n = abs(m[j][i]);
+				if (n > M)
+				{
+					M = n;
+					k = j;
+				}
+			}
+
+			// avoid dividing by zero
+			if (M == T{})
+				return false;
+
+			// swap rows i & k
+			if (i != k)
+			{
+				for (std::size_t j = i; j < N; ++j)
+					std::swap(m[i][j], m[k][j]);
+
+				detail::for_each([i, k](auto & rhs){
+					detail::gauss_helper h{rhs};
+					for (std::size_t c = 0; c < h.columns(); ++c)
+					{
+						std::swap(h.at(i, c), h.at(k, c));
+					}
+				}, rhs...);
+			}
+
+			// make i,i equal 1
+			{
+				auto r = m[i][i];
+				for (std::size_t j = i; j < N; ++j)
+				{
+					m[i][j] /= r;
+				}
+
+				detail::for_each([i, r](auto & rhs){
+					detail::gauss_helper h{rhs};
+					for (std::size_t c = 0; c < h.columns(); ++c)
+					{
+						h.at(i, c) /= r;
+					}
+				}, rhs...);
+			}
+
+			// zero out the column under i,i
+			for (std::size_t k = i + 1; k < N; ++k)
+			{
+				auto r = m[k][i];
+				for (std::size_t j = i; j < N; ++j)
+				{
+					m[k][j] -= m[i][j] * r;
+				}
+
+				detail::for_each([i, k, r](auto & rhs){
+					detail::gauss_helper h{rhs};
+					for (std::size_t c = 0; c < h.columns(); ++c)
+					{
+						h.at(k, c) -= h.at(i, c) * r;
+					}
+				}, rhs...);
+			}
+		}
+
+		// backward elimination
+		for (std::size_t i = N; i --> 0;)
+		{
+			for (std::size_t j = 0; j < i; ++j)
+			{
+				auto r = m[j][i];
+				m[j][i] = T{};
+
+				detail::for_each([i, j, r](auto & rhs){
+					detail::gauss_helper h{rhs};
+					for (std::size_t c = 0; c < h.columns(); ++c)
+					{
+						h.at(j, c) -= r * h.at(i, c);
+					}
+				}, rhs...);
+			}
+		}
+
+		return true;
+	}
+
+	template <typename T, std::size_t N>
+	std::optional<matrix<T, N, N>> inverse (matrix<T, N, N> m)
+	{
+		matrix<T, N, N> r = matrix<T, N, N>::identity();
+		if (!gauss(m, r))
+			return std::nullopt;
+		return r;
+	}
+
+	// Least-squares solution of m*x=b with full-rank m
+	template <typename T, std::size_t N, std::size_t M>
+	std::optional<vector<T, M>> least_squares (matrix<T, N, M> const & m, vector<T, N> const & b)
+	{
+		if constexpr (N == M)
+		{
+			auto rhs = b;
+			if (gauss(m, rhs))
+				return rhs;
+			return {};
+		}
+		else if constexpr (N < M)
+		{
+			auto const mt = transpose(m);
+			auto rhs = b;
+			if (gauss(m * mt, rhs))
+				return mt * rhs;
+			return {};
+		}
+		else // if constexpr (N > M)
+		{
+			auto const mt = transpose(m);
+			auto rhs = mt * b;
+			if (gauss(mt * m, rhs))
+				return rhs;
+			return {};
+		}
+	}
+
+}

libs/geom/include/psemek/geom/gram_schmidt.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+
+namespace psemek::geom
+{
+
+	void gram_schmidt()
+	{}
+
+	template <typename V1, typename ... Vs>
+	void gram_schmidt(V1 & v1, Vs & ... vs)
+	{
+		v1 = normalized(v1);
+
+		((vs -= dot(v1, vs) * v1), ...);
+
+		gram_schmidt(vs...);
+	}
+
+}

libs/geom/include/psemek/geom/homogeneous.hpp

@@ -0,0 +1,45 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/matrix.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t D>
+	vector<T, D+1> homogeneous(vector<T, D> const & v)
+	{
+		vector<T, D+1> result;
+		for (std::size_t i = 0; i < D; ++i)
+			result[i] = v[i];
+		result[D] = T(0);
+		return result;
+	}
+
+	template <typename T, std::size_t D>
+	vector<T, D+1> homogeneous(point<T, D> const & v)
+	{
+		vector<T, D+1> result;
+		for (std::size_t i = 0; i < D; ++i)
+			result[i] = v[i];
+		result[D] = T(1);
+		return result;
+	}
+
+	template <typename T, std::size_t D>
+	matrix<T, D+1, D+1> homogeneous(matrix<T, D, D> const & m)
+	{
+		matrix<T, D+1, D+1> result = matrix<T, D+1, D+1>::zero();
+		for (std::size_t i = 0; i < D; ++i)
+		{
+			for (std::size_t j = 0; j < D; ++j)
+			{
+				result[i][j] = m[i][j];
+			}
+		}
+		result[D][D] = T(1);
+		return result;
+	}
+
+}

libs/geom/include/psemek/geom/incircle.hpp

@@ -0,0 +1,70 @@
+#pragma once
+
+#include <psemek/geom/sign.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/orientation.hpp>
+
+#include <boost/multiprecision/gmp.hpp>
+
+#include <type_traits>
+
+namespace psemek::geom
+{
+
+	template <typename T>
+	std::enable_if_t<!std::is_floating_point_v<T>, sign_t>
+		in_circle (point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2, point<T, 2> const & p3)
+	{
+		auto proj = [](point<T, 2> const & p){
+			auto const x = p[0];
+			auto const y = p[1];
+			return point<T, 3>{ x, y, x*x + y*y };
+		};
+
+		return orientation(proj(p0), proj(p1), proj(p2), proj(p3));
+	}
+
+	template <typename T>
+	std::enable_if_t<std::is_floating_point_v<T>, sign_t>
+		in_circle (point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2, point<T, 2> const & p3)
+	{
+		constexpr T error = std::numeric_limits<T>::epsilon() * T(29) / T(2);
+
+		T const m01 = (p0[0] - p3[0]) * (p1[1] - p3[1]);
+		T const m02 = (p0[0] - p3[0]) * (p2[1] - p3[1]);
+		T const m10 = (p1[0] - p3[0]) * (p0[1] - p3[1]);
+		T const m12 = (p1[0] - p3[0]) * (p2[1] - p3[1]);
+		T const m20 = (p2[0] - p3[0]) * (p0[1] - p3[1]);
+		T const m21 = (p2[0] - p3[0]) * (p1[1] - p3[1]);
+
+		T const d = T(0)
+			+ m01 * p2[0] * p2[0] + m01 * p2[1] * p2[1] - m01 * p3[0] * p3[0] - m01 * p3[1] * p3[1]
+			- m02 * p1[0] * p1[0] - m02 * p1[1] * p1[1] + m02 * p3[0] * p3[0] + m02 * p3[1] * p3[1]
+			- m10 * p2[0] * p2[0] - m10 * p2[1] * p2[1] + m10 * p3[0] * p3[0] + m10 * p3[1] * p3[1]
+			+ m12 * p0[0] * p0[0] + m12 * p0[1] * p0[1] - m12 * p3[0] * p3[0] - m12 * p3[1] * p3[1]
+			+ m20 * p1[0] * p1[0] + m20 * p1[1] * p1[1] - m20 * p3[0] * p3[0] - m20 * p3[1] * p3[1]
+			- m21 * p0[0] * p0[0] - m21 * p0[1] * p0[1] + m21 * p3[0] * p3[0] + m21 * p3[1] * p3[1]
+			;
+
+		T const t = T(0)
+			+ std::abs(m01 * p2[0] * p2[0]) + std::abs(m01 * p2[1] * p2[1]) + std::abs(m01 * p3[0] * p3[0]) + std::abs(m01 * p3[1] * p3[1])
+			+ std::abs(m02 * p1[0] * p1[0]) + std::abs(m02 * p1[1] * p1[1]) + std::abs(m02 * p3[0] * p3[0]) + std::abs(m02 * p3[1] * p3[1])
+			+ std::abs(m10 * p2[0] * p2[0]) + std::abs(m10 * p2[1] * p2[1]) + std::abs(m10 * p3[0] * p3[0]) + std::abs(m10 * p3[1] * p3[1])
+			+ std::abs(m12 * p0[0] * p0[0]) + std::abs(m12 * p0[1] * p0[1]) + std::abs(m12 * p3[0] * p3[0]) + std::abs(m12 * p3[1] * p3[1])
+			+ std::abs(m20 * p1[0] * p1[0]) + std::abs(m20 * p1[1] * p1[1]) + std::abs(m20 * p3[0] * p3[0]) + std::abs(m20 * p3[1] * p3[1])
+			+ std::abs(m21 * p0[0] * p0[0]) + std::abs(m21 * p0[1] * p0[1]) + std::abs(m21 * p3[0] * p3[0]) + std::abs(m21 * p3[1] * p3[1])
+			;
+
+		if (d > t * error)
+			return sign_t::positive;
+		else if (d < - t * error)
+			return sign_t::negative;
+		else
+		{
+			using exact_type = boost::multiprecision::mpq_rational;
+
+			return in_circle(cast<exact_type>(p0), cast<exact_type>(p1), cast<exact_type>(p2), cast<exact_type>(p3));
+		}
+	}
+
+}

libs/geom/include/psemek/geom/intersection.hpp

@@ -0,0 +1,112 @@
+#pragma once
+
+#include <psemek/geom/simplex.hpp>
+#include <psemek/geom/orientation.hpp>
+#include <psemek/geom/contains.hpp>
+
+#include <variant>
+#include <type_traits>
+
+namespace psemek::geom
+{
+
+	// denotes empty intersection
+	struct empty { };
+
+	template <typename Stream>
+	Stream & operator << (Stream & os, empty)
+	{
+		return os << "empty";
+	}
+
+	template <typename Point, typename = std::enable_if_t<Point::dimension == 2>>
+	bool intersect (segment<Point> const & s0, segment<Point> const & s1)
+	{
+		auto const o00 = orientation(s0[0], s0[1], s1[0]);
+		auto const o01 = orientation(s0[0], s0[1], s1[1]);
+		auto const o10 = orientation(s1[0], s1[1], s0[0]);
+		auto const o11 = orientation(s1[0], s1[1], s0[1]);
+
+		return ((o00 != o01) || (o00 == sign_t::zero)) && ((o10 != o11) || (o10 == sign_t::zero));
+	}
+
+	// TODO: robust implementation
+	template <typename Point, typename = std::enable_if_t<Point::dimension == 2>>
+	std::variant<empty, Point, segment<Point>> intersection (segment<Point> s0, segment<Point> s1)
+	{
+		auto const a0 = -det(s1[0] - s0[0], s1[1] - s1[0]);
+		auto const a1 =  det(s0[1] - s0[0], s1[0] - s0[0]);
+		auto const b  = -det(s0[1] - s0[0], s1[1] - s1[0]);
+
+		if (b != 0)
+		{
+			// general case
+
+			auto const t0 = a0 / b;
+			auto const t1 = a1 / b;
+
+			if (t0 < 0 || t0 > 1 || t1 < 0 || t1 > 1)
+				return empty{};
+
+			return s0[0] + t0 * (s0[1] - s0[0]);
+		}
+		else
+		{
+			// collinear segments
+
+			if (a0 != 0)
+			{
+				// segments do not lie on the same line: no intersection
+				return empty{};
+			}
+
+			// if segments are not Y-axis aligned, safe to use X-coordinates to sort them (k = 0)
+			// otherwise use Y-coordinates (k = 1)
+			std::size_t const k = (s0[0][0] != s0[1][0]) ? 0 : 1;
+
+			if (s0[0][k] > s0[1][k])
+				std::swap(s0[0], s0[1]);
+
+			if (s1[0][k] > s1[1][k])
+				std::swap(s1[0], s1[1]);
+
+			auto const r0 = std::max(s0[0][k], s1[0][k]);
+			auto const r1 = std::min(s0[1][k], s1[1][k]);
+
+			if (r0 > r1)
+				return empty{};
+
+			bool const s0_is_first = s0[0][k] < s1[0][k];
+
+			if (r0 == r1)
+			{
+				Point p;
+				p[k] = r0;
+				p[1 - k] = s0_is_first ? s0[1][1 - k] : s1[1][1 - k];
+				return p;
+			}
+			else if (s0_is_first)
+				return segment{ s1[0], s0[1] };
+			else
+				return segment{ s0[0], s1[1] };
+		}
+	}
+
+	template <typename T>
+	bool intersect(triangle<point<T, 2>> const & t0, triangle<point<T, 2>> const & t1)
+	{
+		if (contains(t0, t1[0]) || contains(t0, t1[1]) || contains(t0, t1[2])) return true;
+
+		if (contains(t1, t0[0]) || contains(t1, t0[1]) || contains(t1, t0[2])) return true;
+
+		if (intersect(segment{t0[0], t0[1]}, segment{t1[0], t1[1]})) return true;
+		if (intersect(segment{t0[0], t0[1]}, segment{t1[1], t1[2]})) return true;
+		if (intersect(segment{t0[1], t0[2]}, segment{t1[0], t1[1]})) return true;
+		if (intersect(segment{t0[1], t0[2]}, segment{t1[1], t1[2]})) return true;
+		if (intersect(segment{t0[2], t0[0]}, segment{t1[0], t1[1]})) return true;
+		if (intersect(segment{t0[2], t0[0]}, segment{t1[1], t1[2]})) return true;
+
+		return false;
+	}
+
+}

libs/geom/include/psemek/geom/interval.hpp

@@ -0,0 +1,222 @@
+#pragma once
+
+#include <limits>
+#include <type_traits>
+#include <cmath>
+
+namespace psemek::geom
+{
+
+	// Can be specialized in client code
+	template <typename T>
+	struct limits
+	{
+		static constexpr T min ( )
+		{
+			if constexpr (std::is_floating_point_v<T>)
+			{
+				return -std::numeric_limits<T>::infinity();
+			}
+			else
+			{
+				return std::numeric_limits<T>::min();
+			}
+		}
+
+		static constexpr T max ( )
+		{
+			if constexpr (std::is_floating_point_v<T>)
+			{
+				return std::numeric_limits<T>::infinity();
+			}
+			else
+			{
+				return std::numeric_limits<T>::max();
+			}
+		}
+	};
+
+
+	template <typename T>
+	struct interval_iterator
+	{
+		T value;
+
+		T operator* () const { return value; }
+
+		interval_iterator & operator++ ()
+		{
+			++value;
+			return *this;
+		}
+
+		interval_iterator operator++ (int)
+		{
+			auto copy = *this;
+			++value;
+			return copy;
+		}
+
+		friend bool operator == (interval_iterator const & i1, interval_iterator const & i2)
+		{
+			return i1.value == i2.value;
+		}
+
+		friend bool operator != (interval_iterator const & i1, interval_iterator const & i2)
+		{
+			return !(i1 == i2);
+		}
+	};
+
+	template <typename T>
+	struct interval
+	{
+		T min = limits<T>::max();
+		T max = limits<T>::min();
+
+		static interval singleton (T const & value)
+		{
+			return {value, value};
+		}
+
+		bool empty ( ) const
+		{
+			return min > max;
+		}
+
+		T length ( ) const
+		{
+			return empty() ? T{} : max - min;
+		}
+
+		using iterator = interval_iterator<T>;
+		using const_iterator = iterator;
+
+		iterator begin() const { return {min}; }
+		iterator end() const { return {max}; }
+
+		interval & operator += (T const & delta);
+		interval & operator -= (T const & delta);
+
+		interval & operator &= (T const & a);
+		interval & operator |= (T const & a);
+
+		interval & operator &= (interval const & i);
+		interval & operator |= (interval const & i);
+	};
+
+	template <typename T>
+	interval<T> operator + (interval<T> const & i, T const & delta)
+	{
+		return {i.min + delta, i.max + delta};
+	}
+
+	template <typename T>
+	interval<T> operator + (T const & delta, interval<T> const & i)
+	{
+		return {delta + i.min, delta + i.max};
+	}
+
+	template <typename T>
+	interval<T> operator - (interval<T> const & i, T const & delta)
+	{
+		return {i.min - delta, i.max - delta};
+	}
+
+	template <typename T>
+	interval<T> operator & (interval<T> const & i, T const & a)
+	{
+		using std::min;
+		using std::max;
+		return {max(i.min, a), min(i.max, a)};
+	}
+
+	template <typename T>
+	interval<T> operator & (T const & a, interval<T> const & i)
+	{
+		using std::min;
+		using std::max;
+		return {max(a, i.min), min(a, i.max)};
+	}
+
+	template <typename T>
+	interval<T> operator | (interval<T> const & i, T const & a)
+	{
+		using std::min;
+		using std::max;
+		return {min(i.min, a), max(i.max, a)};
+	}
+
+	template <typename T>
+	interval<T> operator | (T const & a, interval<T> const & i)
+	{
+		using std::min;
+		using std::max;
+		return {min(a, i.min), max(a, i.max)};
+	}
+
+	template <typename T>
+	interval<T> operator & (interval<T> const & i1, interval<T> const & i2)
+	{
+		using std::min;
+		using std::max;
+		return {max(i1.min, i2.min), min(i1.max, i2.max)};
+	}
+
+	template <typename T>
+	interval<T> operator | (interval<T> const & i1, interval<T> const & i2)
+	{
+		using std::min;
+		using std::max;
+		return {min(i1.min, i2.min), max(i1.max, i2.max)};
+	}
+
+	template <typename T>
+	interval<T> & interval<T>::operator += (T const & delta)
+	{
+		return *this = *this + delta;
+	}
+
+	template <typename T>
+	interval<T> & interval<T>::operator -= (T const & delta)
+	{
+		return *this = *this - delta;
+	}
+
+	template <typename T>
+	interval<T> & interval<T>::operator &= (T const & a)
+	{
+		return *this = *this & a;
+	}
+
+	template <typename T>
+	interval<T> & interval<T>::operator |= (T const & a)
+	{
+		return *this = *this | a;
+	}
+
+	template <typename T>
+	interval<T> & interval<T>::operator &= (interval<T> const & i)
+	{
+		return *this = *this & i;
+	}
+
+	template <typename T>
+	interval<T> & interval<T>::operator |= (interval<T> const & i)
+	{
+		return *this = *this | i;
+	}
+
+	template <typename T>
+	T clamp (T x, interval<T> const & i)
+	{
+		return std::max(i.min, std::min(i.max, x));
+	}
+
+	template <typename Stream, typename T>
+	Stream & operator << (Stream & os, interval<T> const & i)
+	{
+		return os << '[' << i.min << " .. " << i.max << ']';
+	}
+
+}

libs/geom/include/psemek/geom/matrix.hpp

@@ -0,0 +1,238 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+
+#include <algorithm>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t R, std::size_t C>
+	struct matrix
+	{
+		static constexpr std::size_t rows = R;
+		static constexpr std::size_t columns = C;
+
+		using scalar_type = T;
+
+		T coords[R * C];
+
+		T * operator [] (std::size_t i)
+		{
+			return coords + C * i;
+		}
+
+		T const * operator [] (std::size_t i) const
+		{
+			return coords + C * i;
+		}
+
+		matrix & operator *= (T const & s);
+		matrix & operator /= (T const & s);
+
+		matrix & operator += (matrix const & v);
+		matrix & operator -= (matrix const & v);
+
+		static matrix zero ( );
+		static matrix identity ( );
+		static matrix scalar (T const & s);
+	};
+
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> matrix<T, R, C>::zero ( )
+	{
+		matrix<T, R, C> m;
+		for (std::size_t i = 0; i < R * C; ++i)
+			m.coords[i] = 0;
+		return m;
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> matrix<T, R, C>::identity ( )
+	{
+		return scalar(T{1});
+	}
+
+	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> m = zero();
+		for (std::size_t i = 0; i < std::min(R, C); ++i)
+			for (std::size_t j = 0; j < C; ++j)
+				m[i][i] = s;
+		return m;
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	bool operator == (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
+	{
+		return std::equal(m1.coords, m1.coords + R * C, m2.coords);
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	bool operator != (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
+	{
+		return !(m1 == m2);
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	bool operator < (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
+	{
+		return std::lexicographical_compare(m1.coords, m1.coords + R * C, m2.coords, m2.coords + R * C);
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	bool operator > (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
+	{
+		return m2 < m1;
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	bool operator <= (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
+	{
+		return !(m2 < m1);
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	bool operator >= (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
+	{
+		return !(m1 < m2);
+	}
+
+	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]);
+		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;
+		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];
+		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;
+		return r;
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> & matrix<T, R, C>::operator *= (T const & s)
+	{
+		*this = *this * s;
+		return *this;
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> & matrix<T, R, C>::operator /= (T const & s)
+	{
+		*this = *this / s;
+		return *this;
+	}
+
+	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];
+		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];
+		return r;
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> & matrix<T, R, C>::operator += (matrix<T, R, C> const & m)
+	{
+		*this = *this + m;
+		return *this;
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> & matrix<T, R, C>::operator -= (matrix<T, R, C> const & m)
+	{
+		*this = *this - m;
+		return *this;
+	}
+
+	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)
+		{
+			r[i] = T{};
+			for (std::size_t j = 0; j < C; ++j)
+				r[i] += m[i][j] * v[j];
+		}
+		return r;
+	}
+
+	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)
+		{
+			r[j] = T{};
+			for (std::size_t i = 0; i < R; ++i)
+				r[i] += v[j] * m[i][j];
+		}
+		return r;
+	}
+
+	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)
+		{
+			for (std::size_t j = 0; j < C; ++j)
+			{
+				r[i][j] = T{};
+
+				for (std::size_t k = 0; k < K; ++k)
+					r[i][j] += m1[i][k] * m2[k][j];
+			}
+		}
+		return r;
+	}
+
+	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)
+				r[j][i] = m[i][j];
+		return r;
+	}
+
+}

libs/geom/include/psemek/geom/orientation.hpp

@@ -0,0 +1,67 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/sign.hpp>
+
+#include <boost/multiprecision/gmp.hpp>
+
+#include <limits>
+#include <type_traits>
+
+namespace psemek::geom
+{
+
+	// TODO: generic implementation
+	template <typename T>
+	std::enable_if_t<!std::is_floating_point_v<T>, sign_t>
+		orientation (point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2)
+	{
+		T const d = det(p1 - p0, p2 - p0);
+
+		if (d > T{})
+			return sign_t::positive;
+		else if (d < T{})
+			return sign_t::negative;
+		else
+			return sign_t::zero;
+	}
+
+	template <typename T>
+	std::enable_if_t<std::is_floating_point_v<T>, sign_t>
+		orientation (point<T, 2> const & p0, point<T, 2> const & p1, point<T, 2> const & p2)
+	{
+		constexpr T error = std::numeric_limits<T>::epsilon() * T(5) / T(2);
+
+		T const d = (p1[0] - p0[0]) * (p2[1] - p0[1])
+			- (p1[1] - p0[1]) * (p2[0] - p0[0]);
+
+		T const t = std::abs((p1[0] - p0[0]) * (p2[1] - p0[1]))
+				+ std::abs((p1[1] - p0[1]) * (p2[0] - p0[0]));
+
+		if (d > t * error)
+			return sign_t::positive;
+		else if (d < - t * error)
+			return sign_t::negative;
+		else
+		{
+			using exact_type = boost::multiprecision::mpq_rational;
+
+			return orientation(cast<exact_type>(p0), cast<exact_type>(p1), cast<exact_type>(p2));
+		}
+	}
+
+	template <typename T>
+	sign_t orientation (point<T, 3> const & p0, point<T, 3> const & p1, point<T, 3> const & p2, point<T, 3> const & p3)
+	{
+		T const d = det(p0 - p3, p1 - p3, p2 - p3);
+
+		if (d > T{})
+			return sign_t::positive;
+		else if (d < T{})
+			return sign_t::negative;
+		else
+			return sign_t::zero;
+	}
+
+}

libs/geom/include/psemek/geom/orthographic.hpp

@@ -0,0 +1,92 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/box.hpp>
+#include <psemek/geom/matrix.hpp>
+#include <psemek/geom/homogeneous.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t D>
+	struct orthographic
+	{
+		using vector_type = vector<T, D>;
+		using point_type = point<T, D>;
+		using homogeneous_matrix_type = matrix<T, D + 1, D + 1>;
+		using box_type = box<T, D>;
+
+		orthographic();
+		orthographic(box_type r);
+
+		box_type box() const;
+		box_type box(box_type r);
+
+		vector_type operator() (vector_type v) const;
+		point_type operator() (point_type p) const;
+
+		homogeneous_matrix_type homogeneous_matrix() const;
+
+		orthographic inverse() const;
+
+	private:
+		box_type r_;
+	};
+
+	template <typename T, std::size_t D>
+	orthographic<T, D>::orthographic()
+	{
+		for (std::size_t i = 0; i < D; ++i)
+		{
+			r_[i] = {T(-1), T(1)};
+		}
+	}
+
+	template <typename T, std::size_t D>
+	orthographic<T, D>::orthographic(box_type r)
+		: r_(r)
+	{ }
+
+	template <typename T, std::size_t D>
+	box<T, D> orthographic<T, D>::box() const
+	{
+		return r_;
+	}
+
+	template <typename T, std::size_t D>
+	box<T, D> orthographic<T, D>::box(box_type r)
+	{
+		std::swap(r, r_);
+		return r;
+	}
+
+	template <typename T, std::size_t D>
+	vector<T, D> orthographic<T, D>::operator() (vector<T, D> v) const
+	{
+		return homogeneous_matrix() * homogeneous(v);
+	}
+
+	template <typename T, std::size_t D>
+	point<T, D> orthographic<T, D>::operator() (point<T, D> p) const
+	{
+		return homogeneous_matrix() * homogeneous(p);
+	}
+
+	template <typename T, std::size_t D>
+	matrix<T, D + 1, D + 1> orthographic<T, D>::homogeneous_matrix() const
+	{
+		auto m = matrix<T, D + 1, D + 1>::zero();
+
+		for (std::size_t d = 0; d < D; ++d)
+			m[d][d] = T(2) / (r_[d].max - r_[d].min);
+
+		for (std::size_t d = 0; d < D; ++d)
+			m[d][D] = - (r_[d].max + r_[d].min) / (r_[d].max - r_[d].min);
+
+		m[D][D] = T(1);
+
+		return m;
+	}
+
+}

libs/geom/include/psemek/geom/permutation.hpp

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/matrix.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t D>
+	struct swap
+	{
+		using vector_type = vector<T, D>;
+		using point_type = point<T, D>;
+		using matrix_type = matrix<T, D, D>;
+
+		swap(std::size_t i, std::size_t j);
+
+		vector_type operator() (vector_type v) const;
+
+		matrix_type matrix() const;
+
+		swap inverse() const;
+
+	private:
+		std::size_t i_, j_;
+	};
+
+	template <typename T, std::size_t D>
+	swap<T, D>::swap (std::size_t i, std::size_t j)
+		: i_(i)
+		, j_(j)
+	{ }
+
+	template <typename T, std::size_t D>
+	vector<T, D> swap<T, D>::operator() (vector_type v) const
+	{
+		std::swap(v[i_], v[j_]);
+		return v;
+	}
+
+	template <typename T, std::size_t D>
+	matrix<T, D, D> swap<T, D>::matrix() const
+	{
+		matrix_type m = matrix_type::identity();
+		m[i_][i_] = 0;
+		m[j_][j_] = 0;
+		m[i_][j_] = 1;
+		m[j_][i_] = 1;
+		return m;
+	}
+
+	template <typename T, std::size_t D>
+	swap<T, D> swap<T, D>::inverse() const
+	{
+		return *this;
+	}
+
+}

libs/geom/include/psemek/geom/perspective.hpp

@@ -0,0 +1,61 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/matrix.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t D>
+	struct perspective;
+
+	template <typename T>
+	struct perspective<T, 3>
+	{
+		using scalar_type = T;
+		using homogeneous_matrix_type = matrix<T, 4, 4>;
+
+		// fov are in radians
+		perspective (scalar_type fov_x, scalar_type fov_y, scalar_type near, scalar_type far);
+
+		perspective (scalar_type left, scalar_type right, scalar_type bottom, scalar_type top, scalar_type near, scalar_type far);
+
+		homogeneous_matrix_type homogeneous_matrix() const;
+
+	private:
+		homogeneous_matrix_type homogeneous_matrix_;
+	};
+
+	template <typename T>
+	perspective<T, 3>::perspective(scalar_type fov_x, scalar_type fov_y, scalar_type near, scalar_type far)
+		: homogeneous_matrix_(homogeneous_matrix_type::zero())
+	{
+		scalar_type dx = 1 / std::tan(fov_x / 2);
+		scalar_type dy = 1 / std::tan(fov_y / 2);
+		homogeneous_matrix_[0][0] = dx;
+		homogeneous_matrix_[1][1] = dy;
+		homogeneous_matrix_[2][2] = (near + far) / (near - far);
+		homogeneous_matrix_[2][3] = 2 * near * far / (near - far);
+		homogeneous_matrix_[3][2] = -1;
+	}
+
+	template <typename T>
+	perspective<T, 3>::perspective(scalar_type left, scalar_type right, scalar_type bottom, scalar_type top, scalar_type near, scalar_type far)
+		: homogeneous_matrix_(homogeneous_matrix_type::zero())
+	{
+		homogeneous_matrix_[0][0] = 2 * near / (right - left);
+		homogeneous_matrix_[1][1] = 2 * near / (top - bottom);
+		homogeneous_matrix_[0][2] = (right + left) / (right - left);
+		homogeneous_matrix_[1][2] = (top + bottom) / (top - bottom);
+		homogeneous_matrix_[2][2] = - (far + near) / (far - near);
+		homogeneous_matrix_[2][3] = - 2 * far * near / (far - near);
+		homogeneous_matrix_[3][2] = -1;
+	}
+
+	template <typename T>
+	matrix<T, 4, 4> perspective<T, 3>::homogeneous_matrix() const
+	{
+		return homogeneous_matrix_;
+	}
+
+}

libs/geom/include/psemek/geom/point.hpp

@@ -0,0 +1,173 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t N>
+	struct point
+	{
+		static constexpr std::size_t dimension = N;
+
+		using scalar_type = T;
+
+		T coords[N];
+
+		point ( ) = default;
+		point (point const &) = default;
+		point (point &) = default;
+		point (point &&) = default;
+		point & operator = (point const &) = default;
+		point & operator = (point &) = default;
+		point & operator = (point &&) = default;
+
+		template <typename ... Args>
+		point (Args && ... args)
+			: coords{ std::forward<Args>(args)... }
+		{
+			static_assert(sizeof...(Args) == N);
+		}
+
+		T & operator[] (std::size_t i)
+		{
+			return coords[i];
+		}
+
+		T const & operator[] (std::size_t i) const
+		{
+			return coords[i];
+		}
+
+		point & operator += (vector<T, N> const & v);
+		point & operator -= (vector<T, N> const & v);
+	};
+
+	template <typename ... Args>
+	point (Args && ...) -> point<std::common_type_t<Args...>, sizeof...(Args)>;
+
+	template <typename T, std::size_t N>
+	bool operator == (point<T, N> const & p1, point<T, N> const & p2)
+	{
+		for (std::size_t i = 0; i < N; ++i)
+			if (p1[i] != p2[i])
+				return false;
+		return true;
+	}
+
+	template <typename T, std::size_t N>
+	bool operator != (point<T, N> const & p1, point<T, N> const & p2)
+	{
+		return !(p1 == p2);
+	}
+
+	template <typename T, std::size_t N>
+	bool operator < (point<T, N> const & p1, point<T, N> const & p2)
+	{
+		return std::lexicographical_compare(p1.coords, p1.coords + N, p2.coords, p2.coords + N);
+	}
+
+	template <typename T, std::size_t N>
+	bool operator > (point<T, N> const & p1, point<T, N> const & p2)
+	{
+		return p2 < p1;
+	}
+
+	template <typename T, std::size_t N>
+	bool operator <= (point<T, N> const & p1, point<T, N> const & p2)
+	{
+		return !(p2 < p1);
+	}
+
+	template <typename T, std::size_t N>
+	bool operator >= (point<T, N> const & p1, point<T, N> const & p2)
+	{
+		return !(p1 < p2);
+	}
+
+	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)
+			r[i] = T1(p[i]);
+		return r;
+	}
+
+	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;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = p[i] + v[i];
+		return r;
+	}
+
+	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;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = v[i] + p[i];
+		return r;
+	}
+
+	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;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = p[i] - v[i];
+		return r;
+	}
+
+	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;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = p1[i] - p2[i];
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	point<T, N> & point<T, N>::operator += (vector<T, N> const & v)
+	{
+		return (*this) = (*this) + v;
+	}
+
+	template <typename T, std::size_t N>
+	point<T, N> & point<T, N>::operator -= (vector<T, N> const & v)
+	{
+		return (*this) = (*this) - v;
+	}
+
+	template <typename T, std::size_t N>
+	point<T, N> lerp (point<T, N> const & p0, point<T, N> const & p1, T const & t)
+	{
+		return p0 + t * (p1 - p0);
+	}
+
+	template <typename T, std::size_t N>
+	T distance_sqr (point<T, N> const & p1, point<T, N> const & p2)
+	{
+		return length_sqr(p2 - p1);
+	}
+
+	template <typename T, std::size_t N>
+	T distance (point<T, N> const & p1, point<T, N> const & p2)
+	{
+		return length(p2 - p1);
+	}
+
+	template <typename Stream, typename T, std::size_t N>
+	Stream & operator << (Stream & os, point<T, N> const & p)
+	{
+		os << '(' << p[0];
+		for (std::size_t i = 1; i < N; ++i)
+			os << ", " << p[i];
+		os << ')';
+		return os;
+	}
+
+}

libs/geom/include/psemek/geom/rotation.hpp

@@ -0,0 +1,184 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/matrix.hpp>
+#include <psemek/geom/homogeneous.hpp>
+
+namespace psemek::geom
+{
+
+	// Rotation in oriented plane (i,j)
+	template <typename T, std::size_t D>
+	struct plane_rotation
+	{
+		using scalar_type = T;
+		using vector_type = geom::vector<T, D>;
+		using point_type = geom::point<T, D>;
+		using matrix_type = geom::matrix<T, D, D>;
+
+		plane_rotation(std::size_t i, std::size_t j, T angle = T(0));
+
+		scalar_type angle() const;
+		scalar_type angle(scalar_type a);
+
+		vector_type operator() (vector_type v) const;
+
+		matrix_type matrix() const;
+
+		plane_rotation inverse() const;
+
+	private:
+		std::size_t const i_;
+		std::size_t const j_;
+		T angle_;
+	};
+
+	// 3D-rotation around an axis
+	template <typename T>
+	struct axis_rotation
+	{
+		using scalar_type = T;
+		using vector_type = geom::vector<T, 3>;
+		using point_type = geom::point<T, 3>;
+		using matrix_type = geom::matrix<T, 3, 3>;
+
+		axis_rotation();
+		axis_rotation(vector_type axis, scalar_type angle);
+
+		vector_type axis() const;
+		vector_type axis(vector_type a);
+
+		scalar_type angle() const;
+		scalar_type angle(scalar_type a);
+
+		vector_type operator() (vector_type v) const;
+
+		matrix_type matrix() const;
+
+		axis_rotation inverse() const;
+
+	private:
+		vector_type axis_;
+		scalar_type angle_;
+	};
+
+	template <typename T, std::size_t D>
+	plane_rotation<T, D>::plane_rotation(std::size_t i, std::size_t j, T angle)
+		: i_(i), j_(j), angle_(angle)
+	{ }
+
+	template <typename T, std::size_t D>
+	T plane_rotation<T, D>::angle() const
+	{
+		return angle_;
+	}
+
+	template <typename T, std::size_t D>
+	T plane_rotation<T, D>::angle(T a)
+	{
+		T t = angle_;
+		angle_ = a;
+		return t;
+	}
+
+	template <typename T, std::size_t D>
+	vector<T, D> plane_rotation<T, D>::operator() (vector_type v) const
+	{
+		T vi = v[i_] * std::cos(angle_) - v[j_] * std::sin(angle_);
+		T vj = v[i_] * std::sin(angle_) + v[j_] * std::cos(angle_);
+		v[i_] = vi;
+		v[j_] = vj;
+		return v;
+	}
+
+	template <typename T, std::size_t D>
+	matrix<T, D, D> plane_rotation<T, D>::matrix() const
+	{
+		matrix_type m = matrix_type::identity();
+		m[i_][i_] =  std::cos(angle_);
+		m[i_][j_] = -std::sin(angle_);
+		m[j_][i_] =  std::sin(angle_);
+		m[j_][j_] =  std::cos(angle_);
+		return m;
+	}
+
+	template <typename T, std::size_t D>
+	plane_rotation<T, D> plane_rotation<T, D>::inverse() const
+	{
+		return {i_, j_, -angle_};
+	}
+
+	template <typename T>
+	axis_rotation<T>::axis_rotation()
+		: axis_rotation(vector_type{T(0), T(0), T(1)}, T(0))
+	{ }
+
+	template <typename T>
+	axis_rotation<T>::axis_rotation(vector_type axis, scalar_type angle)
+		: axis_(axis)
+		, angle_(angle)
+	{ }
+
+	template <typename T>
+	vector<T, 3> axis_rotation<T>::axis() const
+	{
+		return axis_;
+	}
+
+	template <typename T>
+	vector<T, 3> axis_rotation<T>::axis(vector_type a)
+	{
+		auto t = axis_;
+		axis_ = a;
+		return t;
+	}
+
+	template <typename T>
+	T axis_rotation<T>::angle ( ) const
+	{
+		return angle_;
+	}
+
+	template <typename T>
+	T axis_rotation<T>::angle(scalar_type a)
+	{
+		auto t = angle_;
+		angle_ = a;
+		return t;
+	}
+
+	template <typename T>
+	vector<T, 3> axis_rotation<T>::operator() (vector_type v) const
+	{
+		return matrix() * v;
+	}
+
+	template <typename T>
+	matrix<T, 3, 3> axis_rotation<T>::matrix() const
+	{
+		matrix_type m = matrix_type::identity();
+		T c = std::cos(angle_);
+		T s = std::sin(angle_);
+		T x = axis_[0];
+		T y = axis_[1];
+		T z = axis_[2];
+		m[0][0] = c + x * x * (1 - c);
+		m[0][1] = x * y * (1 - c) - z * s;
+		m[0][2] = x * z * (1 - c) + y * s;
+		m[1][0] = y * x * (1 - c) + z * s;
+		m[1][1] = c + y * y * (1 - c);
+		m[1][2] = y * z * (1 - c) - x * s;
+		m[2][0] = z * x * (1 - c) - y * s;
+		m[2][1] = z * y * (1 - c) + x * s;
+		m[2][2] = c + z * z * (1 - c);
+		return m;
+	}
+
+	template <typename T>
+	axis_rotation<T> axis_rotation<T>::inverse() const
+	{
+		return {axis_, -angle_};
+	}
+
+}

libs/geom/include/psemek/geom/scale.hpp

@@ -0,0 +1,93 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/matrix.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t D>
+	struct scale
+	{
+		using vector_type = vector<T, D>;
+		using point_type = point<T, D>;
+		using matrix_type = matrix<T, D, D>;
+
+		scale();
+		scale(T v);
+		scale(vector_type v);
+
+		vector_type vector() const;
+		vector_type vector(vector_type v);
+
+		vector_type operator() (vector_type v) const;
+
+		matrix_type matrix() const;
+
+		scale inverse() const;
+
+	private:
+		vector_type v_;
+	};
+
+	template <typename T, std::size_t D>
+	scale<T, D>::scale()
+	{
+		for (std::size_t i = 0; i < D; ++i)
+			v_[i] = T(1);
+	}
+
+	template <typename T, std::size_t D>
+	scale<T, D>::scale(T v)
+	{
+		for (std::size_t i = 0; i < D; ++i)
+			v_[i] = v;
+	}
+
+	template <typename T, std::size_t D>
+	scale<T, D>::scale(vector_type v)
+		: v_(v)
+	{ }
+
+	template <typename T, std::size_t D>
+	vector<T, D> scale<T, D>::vector() const
+	{
+		return v_;
+	}
+
+	template <typename T, std::size_t D>
+	vector<T, D> scale<T, D>::vector(vector_type v)
+	{
+		auto t = v_;
+		v_ = v;
+		return t;
+	}
+
+	template <typename T, std::size_t D>
+	vector<T, D> scale<T, D>::operator() (vector_type v) const
+	{
+		for (std::size_t i = 0; i < D; ++i)
+			v[i] *= v_[i];
+		return v;
+	}
+
+	template <typename T, std::size_t D>
+	matrix<T, D, D> scale<T, D>::matrix() const
+	{
+		matrix_type m = matrix_type::identity();
+		for (std::size_t i = 0; i < D; ++i)
+			m[i][i] = v_[i];
+		return m;
+	}
+
+	template <typename T, std::size_t D>
+	scale<T, D> scale<T, D>::inverse() const
+	{
+		vector_type v;
+		for (std::size_t i = 0; i < D; ++i)
+			v[i] = T(1) / v_[i];
+		return {v};
+	}
+
+}

libs/geom/include/psemek/geom/sign.hpp

@@ -0,0 +1,26 @@
+#pragma once
+
+namespace psemek::geom
+{
+
+	enum class sign_t : int
+	{
+		positive = 1,
+		zero = 0,
+		negative = -1,
+	};
+
+	template <typename OStream>
+	OStream & operator << (OStream & o, sign_t s)
+	{
+		switch (s)
+		{
+		case sign_t::positive: return o << "positive";
+		case sign_t::zero: return o << "zero";
+		case sign_t::negative: return o << "negative";
+		}
+
+		return o;
+	}
+
+}

libs/geom/include/psemek/geom/simplex.hpp

@@ -0,0 +1,56 @@
+#pragma once
+
+#include <psemek/geom/point.hpp>
+
+#include <type_traits>
+
+namespace psemek::geom
+{
+
+	template <typename Point, std::size_t K>
+	struct simplex
+	{
+		using point_type = Point;
+
+		point_type points[K + 1];
+
+		point_type & operator[] (std::size_t i)
+		{
+			return points[i];
+		}
+
+		point_type const & operator[] (std::size_t i) const
+		{
+			return points[i];
+		}
+	};
+
+	template <typename ... Args>
+	simplex (Args ...) -> simplex<std::common_type_t<Args...>, sizeof...(Args) - 1>;
+
+	template <typename Point>
+	struct segment
+		: simplex<Point, 1>
+	{ };
+
+	template <typename Point>
+	segment (Point, Point) -> segment<Point>;
+
+	template <typename Point>
+	struct triangle
+		: simplex<Point, 2>
+	{ };
+
+	template <typename Point>
+	triangle (Point, Point, Point) -> triangle<Point>;
+
+	template <typename Stream, typename Point, std::size_t K>
+	Stream & operator << (Stream & os, simplex<Point, K> const & s)
+	{
+		os << '(' << s[0];
+		for (std::size_t i = 1; i <= K; ++i)
+			os << ", " << s[i];
+		return os << ')';
+	}
+
+}

libs/geom/include/psemek/geom/translation.hpp

@@ -0,0 +1,85 @@
+#pragma once
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/matrix.hpp>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t D>
+	struct translation
+	{
+		using vector_type = geom::vector<T, D>;
+		using point_type = geom::point<T, D>;
+		using homogeneous_matrix_type = geom::matrix<T, D + 1, D + 1>;
+
+		translation();
+		translation(vector_type v);
+
+		vector_type vector() const;
+		vector_type vector(vector_type v);
+
+		vector_type operator() (vector_type v) const;
+		point_type operator() (point_type p) const;
+
+		homogeneous_matrix_type homogeneous_matrix() const;
+
+		translation inverse() const;
+
+	private:
+		vector_type v_;
+	};
+
+	template <typename T, std::size_t D>
+	translation<T, D>::translation()
+		: translation(vector_type::zero())
+	{ }
+
+	template <typename T, std::size_t D>
+	translation<T, D>::translation(vector_type v)
+		: v_(v)
+	{ }
+
+	template <typename T, std::size_t D>
+	vector<T, D> translation<T, D>::vector() const
+	{
+		return v_;
+	}
+
+	template <typename T, std::size_t D>
+	vector<T, D> translation<T, D>::vector(vector_type v)
+	{
+		auto t = v_;
+		v_ = v;
+		return t;
+	}
+
+	template <typename T, std::size_t D>
+	vector<T, D> translation<T, D>::operator() (vector_type v) const
+	{
+		return v;
+	}
+
+	template <typename T, std::size_t D>
+	point<T, D> translation<T, D>::operator() (point_type p) const
+	{
+		return p + v_;
+	}
+
+	template <typename T, std::size_t D>
+	matrix<T, D + 1, D + 1> translation<T, D>::homogeneous_matrix() const
+	{
+		homogeneous_matrix_type m = homogeneous_matrix_type::identity();
+		for (std::size_t i = 0; i < D; ++i)
+			m[i][D] = v_[i];
+		return m;
+	}
+
+	template <typename T, std::size_t D>
+	translation<T, D> translation<T, D>::inverse() const
+	{
+		return {-v_};
+	}
+
+}

libs/geom/include/psemek/geom/vector.hpp

@@ -0,0 +1,268 @@
+#pragma once
+
+#include <cstddef>
+#include <utility>
+#include <type_traits>
+#include <cmath>
+
+namespace psemek::geom
+{
+
+	template <typename T, std::size_t N>
+	struct vector
+	{
+		static constexpr std::size_t dimension = N;
+
+		using scalar_type = T;
+
+		T coords[N];
+
+		vector ( ) = default;
+		vector (vector const &) = default;
+		vector (vector &) = default;
+		vector (vector &&) = default;
+		vector & operator = (vector const &) = default;
+		vector & operator = (vector &) = default;
+		vector & operator = (vector &&) = default;
+
+		template <typename ... Args>
+		vector (Args && ... args)
+			: coords{ std::forward<Args>(args)... }
+		{
+			static_assert(sizeof...(Args) == N);
+		}
+
+		T & operator[] (std::size_t i)
+		{
+			return coords[i];
+		}
+
+		T const & operator[] (std::size_t i) const
+		{
+			return coords[i];
+		}
+
+		vector & operator *= (T const & s);
+		vector & operator /= (T const & s);
+
+		vector & operator += (vector const & v);
+		vector & operator -= (vector const & v);
+
+		static vector zero();
+	};
+
+	template <typename ... Args>
+	vector (Args && ...) -> vector<std::common_type_t<Args...>, sizeof...(Args)>;
+
+	template <typename T, std::size_t N>
+	bool operator == (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		return std::equal(v1.coords, v1.coords + N, v2.coords);
+	}
+
+	template <typename T, std::size_t N>
+	bool operator != (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		return !(v1 == v2);
+	}
+
+	template <typename T, std::size_t N>
+	bool operator < (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		return std::lexicographical_compare(v1.coords, v1.coords + N, v2.coords, v2.coords + N);
+	}
+
+	template <typename T, std::size_t N>
+	bool operator > (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		return v2 < v1;
+	}
+
+	template <typename T, std::size_t N>
+	bool operator <= (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		return !(v2 < v1);
+	}
+
+	template <typename T, std::size_t N>
+	bool operator >= (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		return !(v1 < v2);
+	}
+
+	template <typename T1, typename T, std::size_t N>
+	vector<T1, N> cast (vector<T, N> const & v)
+	{
+		vector<T1, N> r;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = static_cast<T1>(v[i]);
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> operator * (vector<T, N> const & v, T const & s)
+	{
+		vector<T, N> r;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = v[i] * s;
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> operator * (T const & s, vector<T, N> const & v)
+	{
+		vector<T, N> r;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = s * v[i];
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> operator / (vector<T, N> const & v, T const & s)
+	{
+		vector<T, N> r;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = v[i] / s;
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> operator - (vector<T, N> const & v)
+	{
+		vector<T, N> r;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = -v[i];
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> operator + (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		vector<T, N> r;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = v1[i] + v2[i];
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> operator - (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		vector<T, N> r;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = v1[i] - v2[i];
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> & vector<T, N>::operator *= (T const & s)
+	{
+		return (*this) = (*this) * s;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> & vector<T, N>::operator /= (T const & s)
+	{
+		return (*this) = (*this) / s;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> & vector<T, N>::operator += (vector<T, N> const & v)
+	{
+		return (*this) = (*this) + v;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> & vector<T, N>::operator -= (vector<T, N> const & v)
+	{
+		return (*this) = (*this) - v;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> vector<T, N>::zero()
+	{
+		vector<T, N> r;
+		for (std::size_t i = 0; i < N; ++i)
+			r[i] = T(0);
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	T dot (vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		T r{};
+		for (std::size_t i = 0; i < N; ++i)
+			r += v1[i] * v2[i];
+		return r;
+	}
+
+	template <typename T, std::size_t N>
+	T length_sqr (vector<T, N> const & v)
+	{
+		return dot(v, v);
+	}
+
+	template <typename T, std::size_t N>
+	T length (vector<T, N> const & v)
+	{
+		return std::sqrt(length_sqr(v));
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> normalized (vector<T, N> const & v)
+	{
+		return v / length(v);
+	}
+
+	// TODO: generic implementation
+	template <typename T>
+	T det (vector<T, 2> const & v0, vector<T, 2> const & v1)
+	{
+		return v0[0] * v1[1] - v0[1] * v1[0];
+	}
+
+	template <typename T>
+	T det (vector<T, 3> const & v0, vector<T, 3> const & v1, vector<T, 3> const & v2)
+	{
+		return
+			+ v0[0] * v1[1] * v2[2]
+			- v0[0] * v1[2] * v2[1]
+			- v0[1] * v1[0] * v2[2]
+			+ v0[1] * v1[2] * v2[0]
+			+ v0[2] * v1[0] * v2[1]
+			- v0[2] * v1[1] * v2[0]
+			;
+	}
+
+	template <typename T>
+	vector<T, 3> cross (vector<T, 3> const & v0, vector<T, 3> const & v1)
+	{
+		return vector<T, 3>{
+			v0[1] * v1[2] - v0[2] * v1[1],
+			v0[2] * v1[0] - v0[0] * v1[2],
+			v0[0] * v1[1] - v0[1] * v1[0],
+		};
+	}
+
+	template <typename T>
+	T lerp (T const & x0, T const & x1, T const & t)
+	{
+		return x0 * (T(1) - t) + x1 * t;
+	}
+
+	template <typename T, std::size_t N>
+	vector<T, N> lerp (vector<T, N> const & v0, vector<T, N> const & v1, T const & t)
+	{
+		return v0 * (T(1) - t) + v1 * t;
+	}
+
+	template <typename Stream, typename T, std::size_t N>
+	Stream & operator << (Stream & os, vector<T, N> const & v)
+	{
+		os << '(' << v[0];
+		for (std::size_t i = 1; i < N; ++i)
+			os << ", " << v[i];
+		os << ')';
+		return os;
+	}
+
+}

libs/geom/source/camera.cpp

@@ -0,0 +1,86 @@
+#include <psemek/geom/camera.hpp>
+#include <psemek/geom/gauss.hpp>
+#include <psemek/geom/perspective.hpp>
+#include <psemek/geom/translation.hpp>
+#include <psemek/geom/rotation.hpp>
+#include <psemek/geom/homogeneous.hpp>
+
+#include <cmath>
+
+namespace psemek::geom
+{
+
+	matrix<float, 4, 4> camera::transform() const
+	{
+		return projection() * view();
+	}
+
+	point<float, 3> camera::position() const
+	{
+		vector<float, 4> b{ 0.f, 0.f, 0.f, 1.f };
+		gauss(view(), b);
+		return { b[0], b[1], b[2] };
+	}
+
+	vector<float, 3> camera::axis_x() const
+	{
+		vector<float, 4> b{ 1.f, 0.f, 0.f, 0.f };
+		gauss(view(), b);
+		return { b[0], b[1], b[2] };
+	}
+
+	vector<float, 3> camera::axis_y() const
+	{
+		vector<float, 4> b{ 0.f, 1.f, 0.f, 0.f };
+		gauss(view(), b);
+		return { b[0], b[1], b[2] };
+	}
+
+	vector<float, 3> camera::axis_z() const
+	{
+		vector<float, 4> b{ 0.f, 0.f, 1.f, 0.f };
+		gauss(view(), b);
+		return { b[0], b[1], b[2] };
+	}
+
+	vector<float, 3> camera::direction(float x, float y) const
+	{
+		vector<float, 4> b{ x, y, -1.f, 1.f};
+		gauss(transform(), b);
+
+		point<float, 3> p{ b[0] / b[3], b[1] / b[3], b[2] / b[3] };
+		return p - position();
+	}
+
+	std::array<vector<float, 4>, 6> camera::clip_planes() const
+	{
+		auto const m = transpose(transform());
+
+		std::array<vector<float, 4>, 6> p;
+
+		p[0] = m * vector{1.f, 0.f, 0.f, 1.f};
+		p[1] = m * vector{-1.f, 0.f, 0.f, 1.f};
+		p[2] = m * vector{0.f, 1.f, 0.f, 1.f};
+		p[3] = m * vector{0.f, -1.f, 0.f, 1.f};
+		p[4] = m * vector{0.f, 0.f, 1.f, 1.f};
+		p[5] = m * vector{0.f, 0.f, -1.f, 1.f};
+
+		return p;
+	}
+
+	matrix<float, 4, 4> perspective_camera::projection() const
+	{
+		return perspective<float, 3>(fov_x, fov_y, near_clip, far_clip).homogeneous_matrix();
+	}
+
+	matrix<float, 4, 4> spherical_camera::view() const
+	{
+		return
+			  translation<float, 3>({0.f, 0.f, -distance}).homogeneous_matrix()
+			* homogeneous(plane_rotation<float, 3>(1, 2, elevation_angle).matrix())
+			* homogeneous(plane_rotation<float, 3>(2, 0, azimuthal_angle).matrix())
+			* translation<float, 3>({ -target[0], -target[1], -target[2] }).homogeneous_matrix()
+			;
+	}
+
+}

libs/gfx/CMakeLists.txt

@@ -0,0 +1,9 @@
+set(OpenGL_GL_PREFERENCE GLVND)
+find_package(OpenGL REQUIRED)
+
+file(GLOB_RECURSE PSEMEK_GFX_HEADERS RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "include/*.hpp")
+file(GLOB_RECURSE PSEMEK_GFX_SOURCES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "source/*.cpp")
+
+add_library(gfx ${PSEMEK_GFX_HEADERS} ${PSEMEK_GFX_SOURCES})
+target_include_directories(gfx PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
+target_link_libraries(gfx PUBLIC util geom OpenGL::GL)

libs/gfx/include/psemek/gfx/buffer.hpp

@@ -0,0 +1,47 @@
+#pragma once
+
+#include <psemek/gfx/gl.hpp>
+
+#include <vector>
+
+namespace psemek::gfx
+{
+
+	struct buffer
+	{
+		buffer();
+		buffer(buffer const &) = delete;
+		buffer(buffer &&);
+
+		buffer & operator = (buffer const &) = delete;
+		buffer & operator = (buffer &&);
+
+		~buffer();
+
+		static buffer null();
+
+		GLuint id() const { return id_; }
+
+		void bind() const;
+
+		void load(void const * data, std::size_t size, GLenum usage = gl::STREAM_DRAW);
+
+		template <typename T>
+		void load(T const * data, std::size_t size, GLenum usage = gl::STREAM_DRAW)
+		{
+			load(static_cast<void const *>(data), size * sizeof(T), usage);
+		}
+
+		template <typename T>
+		void load(std::vector<T> const & data, GLenum usage = gl::STREAM_DRAW)
+		{
+			load(data.data(), data.size(), usage);
+		}
+
+	private:
+		GLuint id_;
+
+		explicit buffer(GLuint id);
+	};
+
+}

libs/gfx/include/psemek/gfx/framebuffer.hpp

@@ -0,0 +1,41 @@
+#pragma once
+
+#include <psemek/gfx/gl.hpp>
+#include <psemek/gfx/texture.hpp>
+#include <psemek/gfx/renderbuffer.hpp>
+
+namespace psemek::gfx
+{
+
+	struct framebuffer
+	{
+		framebuffer();
+		framebuffer(framebuffer const &) = delete;
+		framebuffer(framebuffer &&);
+
+		framebuffer & operator = (framebuffer const &) = delete;
+		framebuffer & operator = (framebuffer &&);
+
+		~framebuffer();
+
+		static framebuffer null();
+
+		GLuint id() const { return id_; }
+
+		void bind() const;
+
+		void color(texture_2d const & tex, int attachment = 0);
+		void color(renderbuffer const & rb, int attachment = 0);
+		void depth(renderbuffer const & rb);
+
+		GLenum status() const;
+		bool complete() const;
+		void assert_complete() const;
+
+	private:
+		GLuint id_;
+
+		explicit framebuffer(GLuint id);
+	};
+
+}

libs/gfx/include/psemek/gfx/fullscreen.hpp

@@ -0,0 +1,10 @@
+#pragma once
+
+#include <psemek/gfx/mesh.hpp>
+
+namespace psemek::gfx
+{
+
+	indexed_mesh const & fullscreen_quad();
+
+}

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

@@ -0,0 +1,384 @@
+#pragma once
+
+#include <psemek/gfx/gl.hpp>
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+
+#include <cstddef>
+#include <tuple>
+#include <type_traits>
+#include <initializer_list>
+#include <algorithm>
+#include <vector>
+
+namespace psemek::gfx
+{
+
+	template <typename T>
+	struct normalized
+	{};
+
+	template <typename Attrib>
+	struct attrib_traits;
+
+	template <>
+	struct attrib_traits<std::uint8_t>
+	{
+		using attrib_type = std::uint8_t;
+
+		static constexpr GLint size = 1;
+		static constexpr GLenum type = gl::UNSIGNED_BYTE;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <>
+	struct attrib_traits<std::int8_t>
+	{
+		using attrib_type = std::int8_t;
+
+		static constexpr GLint size = 1;
+		static constexpr GLenum type = gl::BYTE;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <>
+	struct attrib_traits<std::uint16_t>
+	{
+		using attrib_type = std::uint16_t;
+
+		static constexpr GLint size = 1;
+		static constexpr GLenum type = gl::UNSIGNED_SHORT;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <>
+	struct attrib_traits<std::int16_t>
+	{
+		using attrib_type = std::int16_t;
+
+		static constexpr GLint size = 1;
+		static constexpr GLenum type = gl::SHORT;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <>
+	struct attrib_traits<std::uint32_t>
+	{
+		using attrib_type = std::uint32_t;
+
+		static constexpr GLint size = 1;
+		static constexpr GLenum type = gl::UNSIGNED_INT;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <>
+	struct attrib_traits<std::int32_t>
+	{
+		using attrib_type = std::int32_t;
+
+		static constexpr GLint size = 1;
+		static constexpr GLenum type = gl::INT;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <>
+	struct attrib_traits<float>
+	{
+		using attrib_type = float;
+
+		static constexpr GLint size = 1;
+		static constexpr GLenum type = gl::FLOAT;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <>
+	struct attrib_traits<double>
+	{
+		using attrib_type = double;
+
+		static constexpr GLint size = 1;
+		static constexpr GLenum type = gl::DOUBLE;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <typename T, std::size_t N>
+	struct attrib_traits<geom::vector<T, N>>
+	{
+		using attrib_type = geom::vector<T, N>;
+
+		static constexpr GLint size = N;
+		static constexpr GLenum type = attrib_traits<T>::type;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <typename T, std::size_t N>
+	struct attrib_traits<geom::point<T, N>>
+	{
+		using attrib_type = geom::point<T, N>;
+
+		static constexpr GLint size = N;
+		static constexpr GLenum type = attrib_traits<T>::type;
+		static constexpr GLboolean normalized = gl::FALSE_;
+	};
+
+	template <typename T>
+	struct attrib_traits<normalized<T>>
+	{
+		using attrib_type = T;
+
+		static constexpr GLint size = attrib_traits<T>::size;
+		static constexpr GLenum type = attrib_traits<T>::type;
+		static constexpr GLboolean normalized = gl::TRUE_;
+	};
+
+	template <typename ... Attribs>
+	struct vertex
+	{
+		static constexpr std::size_t size = (0 + ... + sizeof(typename attrib_traits<Attribs>::attrib_type));
+
+		using attribs = std::tuple<Attribs...>;
+		using attrib_types = std::tuple<typename attrib_traits<Attribs>::attrib_type...>;
+
+		template <std::size_t I>
+		static constexpr std::size_t offset()
+		{
+			return offset_impl(std::make_index_sequence<I>{});
+		}
+
+		template <std::size_t I>
+		auto & get()
+		{
+			using T = std::tuple_element_t<I, attrib_types>;
+			return *reinterpret_cast<T *>(reinterpret_cast<char *>(&storage) + offset<I>());
+		}
+
+		template <std::size_t I>
+		auto const & get() const
+		{
+			using T = std::tuple_element_t<I, attrib_types>;
+			return *reinterpret_cast<T const *>(reinterpret_cast<char const *>(&storage) + offset<I>());
+		}
+
+		vertex() = default;
+
+		vertex(vertex const & other)
+		{
+			copy_impl(std::make_index_sequence<sizeof...(Attribs)>{}, other);
+		}
+
+		vertex (typename attrib_traits<Attribs>::attrib_type const & ... args)
+		{
+			construct_impl(std::make_index_sequence<sizeof...(Attribs)>{}, args...);
+		}
+
+		vertex & operator = (vertex const & other)
+		{
+			if (this != &other) copy_impl(std::make_index_sequence<sizeof...(Attribs)>{}, other);
+			return *this;
+		}
+
+		std::aligned_storage_t<size, std::max(std::initializer_list<std::size_t>{alignof(Attribs)...})> storage;
+
+	private:
+
+		template <std::size_t ... Is>
+		static constexpr std::size_t offset_impl(std::index_sequence<Is...>)
+		{
+			return (0 + ... + sizeof(std::tuple_element_t<Is, attrib_types>));
+		}
+
+		template <std::size_t ... Is>
+		void construct_impl(std::index_sequence<Is...>, typename attrib_traits<Attribs>::attrib_type const & ... args)
+		{
+			((get<Is>() = args), ...);
+		}
+
+		template <std::size_t ... Is>
+		void copy_impl(std::index_sequence<Is...>, vertex const & other)
+		{
+			((get<Is>() = other.get<Is>()), ...);
+		}
+	};
+
+	namespace detail
+	{
+
+		template <std::size_t ... Is>
+		void enable_attribs(std::index_sequence<Is...>)
+		{
+			int x[] { ((void) gl::EnableVertexAttribArray(Is), 0) ... };
+			(void) x;
+		}
+
+		template <typename Vertex, std::size_t I>
+		void attrib_pointer()
+		{
+			using traits = attrib_traits<std::tuple_element_t<I, typename Vertex::attribs>>;
+			gl::VertexAttribPointer(I, traits::size, traits::type, traits::normalized, sizeof(Vertex), reinterpret_cast<char const *>(Vertex::template offset<I>()));
+		}
+
+		template <typename Vertex, std::size_t ... Is>
+		void attrib_pointers(std::index_sequence<Is...>)
+		{
+			int x[] { ((void) attrib_pointer<Vertex, Is>(), 0) ... };
+			(void) x;
+		}
+
+		template <typename Vertex>
+		struct mesh_setup;
+
+		template <typename ... Attribs>
+		struct mesh_setup<vertex<Attribs...>>
+		{
+			static void setup()
+			{
+				using indices = std::make_index_sequence<sizeof...(Attribs)>;
+
+				enable_attribs(indices{});
+				attrib_pointers<vertex<Attribs...>>(indices{});
+			}
+		};
+
+		template <typename T>
+		struct gl_type;
+
+		template <>
+		struct gl_type<std::uint8_t>
+		{
+			static constexpr GLenum value = gl::UNSIGNED_BYTE;
+		};
+
+		template <>
+		struct gl_type<std::uint16_t>
+		{
+			static constexpr GLenum value = gl::UNSIGNED_SHORT;
+		};
+
+		template <>
+		struct gl_type<std::uint32_t>
+		{
+			static constexpr GLenum value = gl::UNSIGNED_INT;
+		};
+
+		template <typename T>
+		static constexpr GLenum gl_type_v = gl_type<T>::value;
+
+	}
+
+	struct mesh
+	{
+		static mesh null();
+
+		mesh();
+		mesh(mesh &&);
+		mesh(mesh const &) = delete;
+
+		mesh & operator =(mesh &&);
+		mesh & operator =(mesh const &) = delete;
+
+		~mesh();
+
+		template <typename Vertex>
+		void setup()
+		{
+			static_assert(sizeof(Vertex) == Vertex::size);
+
+			gl::BindVertexArray(array_);
+			gl::BindBuffer(gl::ARRAY_BUFFER, buffer_);
+
+			detail::mesh_setup<Vertex>::setup();
+		}
+
+		template <typename Vertex>
+		void load(std::vector<Vertex> const & vertices, GLenum usage = gl::STREAM_DRAW)
+		{
+			gl::BindBuffer(gl::ARRAY_BUFFER, buffer_);
+			gl::BufferData(gl::ARRAY_BUFFER, vertices.size() * sizeof(Vertex), vertices.data(), usage);
+
+			count_ = vertices.size();
+		}
+
+		void draw(GLenum mode) const
+		{
+			if (count_ == 0) return;
+
+			gl::BindVertexArray(array_);
+			gl::DrawArrays(mode, 0, count_);
+		}
+
+		GLsizei count() const
+		{
+			return count_;
+		}
+
+	private:
+		GLuint array_;
+		GLuint buffer_;
+
+		GLsizei count_ = 0;
+
+		mesh(int);
+	};
+
+	struct indexed_mesh
+	{
+		static indexed_mesh null();
+
+		indexed_mesh();
+		indexed_mesh(indexed_mesh &&);
+		indexed_mesh(indexed_mesh const &) = delete;
+
+		indexed_mesh & operator =(indexed_mesh &&);
+		indexed_mesh & operator =(indexed_mesh const &) = delete;
+
+		~indexed_mesh();
+
+		template <typename Vertex>
+		void setup()
+		{
+			gl::BindVertexArray(array_);
+			gl::BindBuffer(gl::ARRAY_BUFFER, buffer_);
+
+			detail::mesh_setup<Vertex>::setup();
+
+			gl::BindBuffer(gl::ELEMENT_ARRAY_BUFFER, index_buffer_);
+		}
+
+		template <typename Vertex, typename Index>
+		void load(std::vector<Vertex> const & vertices, std::vector<Index> const & indices, GLenum usage = gl::STREAM_DRAW)
+		{
+			gl::BindBuffer(gl::ARRAY_BUFFER, buffer_);
+			gl::BufferData(gl::ARRAY_BUFFER, vertices.size() * sizeof(Vertex), vertices.data(), usage);
+			gl::BindBuffer(gl::ARRAY_BUFFER, index_buffer_);
+			gl::BufferData(gl::ARRAY_BUFFER, indices.size() * sizeof(Index), indices.data(), usage);
+
+			count_ = indices.size();
+			index_type_ = detail::gl_type_v<Index>;
+		}
+
+		void draw(GLenum mode) const
+		{
+			if (count_ == 0) return;
+
+			gl::BindVertexArray(array_);
+			gl::DrawElements(mode, count_, index_type_, nullptr);
+		}
+
+		GLsizei count() const
+		{
+			return count_;
+		}
+
+	private:
+		GLuint array_;
+		GLuint buffer_;
+		GLuint index_buffer_;
+
+		GLsizei count_ = 0;
+		GLenum index_type_ = 0;
+
+		indexed_mesh(int);
+	};
+
+}

libs/gfx/include/psemek/gfx/program.hpp

@@ -0,0 +1,94 @@
+#pragma once
+
+#include <psemek/gfx/gl.hpp>
+
+#include <psemek/geom/vector.hpp>
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/matrix.hpp>
+#include <psemek/geom/interval.hpp>
+
+#include <unordered_map>
+#include <string_view>
+
+namespace psemek::gfx
+{
+
+	struct program
+	{
+		program(std::string_view vertex_source, std::string_view fragment_source);
+		program(std::string_view vertex_source, std::string_view geometry_source, std::string_view fragment_source);
+
+		GLuint id() const;
+
+		void bind() const;
+
+		GLint location(char const * name) const;
+
+		struct uniform_proxy
+		{
+			uniform_proxy(GLint location)
+				: location_(location)
+			{}
+
+			GLint location() const { return location_; }
+
+			void operator = (int i);
+			void operator = (unsigned int i);
+			void operator = (float f);
+
+			void operator = (geom::vector<int, 1> const & v);
+			void operator = (geom::vector<int, 2> const & v);
+			void operator = (geom::vector<int, 3> const & v);
+			void operator = (geom::vector<int, 4> const & v);
+
+			void operator = (geom::vector<unsigned int, 1> const & v);
+			void operator = (geom::vector<unsigned int, 2> const & v);
+			void operator = (geom::vector<unsigned int, 3> const & v);
+			void operator = (geom::vector<unsigned int, 4> const & v);
+
+			void operator = (geom::vector<float, 1> const & v);
+			void operator = (geom::vector<float, 2> const & v);
+			void operator = (geom::vector<float, 3> const & v);
+			void operator = (geom::vector<float, 4> const & v);
+
+			void operator = (geom::point<int, 1> const & v);
+			void operator = (geom::point<int, 2> const & v);
+			void operator = (geom::point<int, 3> const & v);
+			void operator = (geom::point<int, 4> const & v);
+
+			void operator = (geom::point<unsigned int, 1> const & v);
+			void operator = (geom::point<unsigned int, 2> const & v);
+			void operator = (geom::point<unsigned int, 3> const & v);
+			void operator = (geom::point<unsigned int, 4> const & v);
+
+			void operator = (geom::point<float, 1> const & v);
+			void operator = (geom::point<float, 2> const & v);
+			void operator = (geom::point<float, 3> const & v);
+			void operator = (geom::point<float, 4> const & v);
+
+			void operator = (geom::matrix<float, 2, 2> const & m);
+			void operator = (geom::matrix<float, 2, 3> const & m);
+			void operator = (geom::matrix<float, 2, 4> const & m);
+			void operator = (geom::matrix<float, 3, 2> const & m);
+			void operator = (geom::matrix<float, 3, 3> const & m);
+			void operator = (geom::matrix<float, 3, 4> const & m);
+			void operator = (geom::matrix<float, 4, 2> const & m);
+			void operator = (geom::matrix<float, 4, 3> const & m);
+			void operator = (geom::matrix<float, 4, 4> const & m);
+
+			void operator = (geom::interval<int> const & i);
+			void operator = (geom::interval<unsigned int> const & i);
+			void operator = (geom::interval<float> const & i);
+
+		private:
+			GLint location_;
+		};
+
+		uniform_proxy operator[] (char const * name) const;
+
+	private:
+		GLuint program_;
+		mutable std::unordered_map<std::string, GLint> uniforms_;
+	};
+
+}

libs/gfx/include/psemek/gfx/render_to_texture.hpp

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <psemek/gfx/texture.hpp>
+#include <psemek/gfx/framebuffer.hpp>
+
+namespace psemek::gfx
+{
+
+	struct render_to_texture
+	{
+		render_to_texture()
+		{
+			texture_.bind();
+			texture_.unbind();
+			framebuffer_.color(texture_);
+		}
+
+		struct framebuffer & framebuffer() { return framebuffer_; }
+		texture_2d & texture() { return texture_; }
+
+		void bind() const
+		{
+			framebuffer_.bind();
+		}
+
+		void unbind() const
+		{
+			framebuffer_.unbind();
+		}
+
+	private:
+		texture_2d texture_;
+		struct framebuffer framebuffer_;
+	};
+
+}

libs/gfx/include/psemek/gfx/renderbuffer.hpp

@@ -0,0 +1,33 @@
+#pragma once
+
+#include <psemek/gfx/gl.hpp>
+
+namespace psemek::gfx
+{
+
+	struct renderbuffer
+	{
+		renderbuffer();
+		renderbuffer(renderbuffer const &) = delete;
+		renderbuffer(renderbuffer &&);
+
+		renderbuffer & operator = (renderbuffer const &) = delete;
+		renderbuffer & operator = (renderbuffer &&);
+
+		~renderbuffer();
+
+		static renderbuffer null();
+
+		GLuint id() const { return id_; }
+
+		void bind() const;
+
+		void storage(GLenum internal_format, GLsizei width, GLsizei height);
+
+	private:
+		GLuint id_;
+
+		explicit renderbuffer(GLuint id);
+	};
+
+}

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

@@ -0,0 +1,115 @@
+#pragma once
+
+#include <psemek/gfx/gl.hpp>
+#include <psemek/util/pixmap.hpp>
+#include <psemek/geom/vector.hpp>
+
+namespace psemek::gfx
+{
+
+	template <typename Pixel>
+	struct pixel_traits;
+
+	template <>
+	struct pixel_traits<std::uint8_t>
+	{
+		static constexpr GLenum internal_format = gl::R8;
+		static constexpr GLenum format = gl::RED;
+		static constexpr GLenum type = gl::UNSIGNED_BYTE;
+	};
+
+	template <>
+	struct pixel_traits<std::array<std::uint8_t, 3>>
+	{
+		static constexpr GLenum internal_format = gl::RGB8;
+		static constexpr GLenum format = gl::RGB;
+		static constexpr GLenum type = gl::UNSIGNED_BYTE;
+	};
+
+	template <>
+	struct pixel_traits<std::array<std::uint8_t, 4>>
+	{
+		static constexpr GLenum internal_format = gl::RGBA8;
+		static constexpr GLenum format = gl::RGBA;
+		static constexpr GLenum type = gl::UNSIGNED_BYTE;
+	};
+
+	template <>
+	struct pixel_traits<float>
+	{
+		static constexpr GLenum internal_format = gl::R32F;
+		static constexpr GLenum format = gl::RED;
+		static constexpr GLenum type = gl::FLOAT;
+	};
+
+	struct texture_2d
+	{
+		texture_2d();
+		texture_2d(texture_2d const &) = delete;
+		texture_2d(texture_2d &&);
+
+		texture_2d & operator = (texture_2d const &) = delete;
+		texture_2d & operator = (texture_2d &&);
+
+		~texture_2d();
+
+		static texture_2d null();
+
+		GLuint id() const { return id_; }
+
+		void bind() const;
+
+		explicit operator bool () const { return id_ != 0; }
+
+		int width() const { return width_; }
+		int height() const { return height_; }
+
+		geom::vector<int, 2> size() const { return {width_, height_}; }
+
+		void load(GLint internal_format, std::size_t width, std::size_t height, GLenum format, GLenum type, const void * data);
+
+		template <typename Pixel>
+		void load(std::size_t width, std::size_t height, Pixel const * data = nullptr)
+		{
+			using traits = pixel_traits<Pixel>;
+			load(traits::internal_format, width, height, traits::format, traits::type, data);
+		}
+
+		template <typename Pixel>
+		void load(util::basic_pixmap<Pixel> const & p)
+		{
+			load(p.width(), p.height(), p.data());
+		}
+
+		void pixels(GLenum format, GLenum type, void * data) const;
+
+		template <typename Pixmap>
+		Pixmap pixels() const
+		{
+			using traits = pixel_traits<typename Pixmap::pixel_type>;
+			Pixmap p(width_, height_);
+			pixels(traits::format, traits::type, p.data());
+			return p;
+		}
+
+		static texture_2d from_data(GLint internal_format, std::size_t width, std::size_t height, GLenum format, GLenum type, const void * data);
+
+		template <typename Pixmap>
+		static texture_2d from_pixmap(Pixmap const & p)
+		{
+			texture_2d t;
+			t.load(p);
+			return t;
+		}
+
+		void generate_mipmap();
+
+	private:
+		GLuint id_;
+		int width_ = 0;
+		int height_ = 0;
+
+		texture_2d(GLuint id);
+	};
+
+}

libs/gfx/source/gfx/buffer.cpp

@@ -0,0 +1,52 @@
+#include <psemek/gfx/buffer.hpp>
+
+namespace psemek::gfx
+{
+
+	buffer::buffer()
+	{
+		gl::GenBuffers(1, &id_);
+	}
+
+	buffer::buffer(buffer && other)
+		: id_(other.id_)
+	{
+		other.id_ = 0;
+	}
+
+	buffer & buffer::operator = (buffer && other)
+	{
+		if (this == &other) return *this;
+
+		gl::DeleteBuffers(1, &id_);
+		id_ = other.id_;
+		other.id_ = 0;
+		return *this;
+	}
+
+	buffer::~buffer()
+	{
+		gl::DeleteBuffers(1, &id_);
+	}
+
+	buffer buffer::null()
+	{
+		return buffer(0);
+	}
+
+	void buffer::bind() const
+	{
+		gl::BindBuffer(gl::ARRAY_BUFFER, id_);
+	}
+
+	void buffer::load(void const * data, std::size_t size, GLenum usage)
+	{
+		bind();
+		gl::BufferData(gl::ARRAY_BUFFER, size, data, usage);
+	}
+
+	buffer::buffer(GLuint id)
+		: id_(id)
+	{}
+
+}

libs/gfx/source/gfx/framebuffer.cpp

@@ -0,0 +1,95 @@
+#include <psemek/gfx/framebuffer.hpp>
+
+namespace psemek::gfx
+{
+
+	static std::string framebuffer_status_string(GLenum status)
+	{
+		switch (status)
+		{
+		case gl::FRAMEBUFFER_UNDEFINED: return "GL_FRAMEBUFFER_UNDEFINED";
+		case gl::FRAMEBUFFER_INCOMPLETE_ATTACHMENT: return "GL_FRAMEBUFFER_INCOMPLETE_ATTACHMENT";
+		case gl::FRAMEBUFFER_INCOMPLETE_MISSING_ATTACHMENT: return "GL_FRAMEBUFFER_INCOMPLETE_MISSING_ATTACHMENT";
+		case gl::FRAMEBUFFER_UNSUPPORTED: return "GL_FRAMEBUFFER_UNSUPPORTED";
+		case gl::FRAMEBUFFER_INCOMPLETE_MULTISAMPLE: return "GL_FRAMEBUFFER_INCOMPLETE_MULTISAMPLE";
+		}
+
+		return "(unknown)";
+	}
+
+	framebuffer::framebuffer()
+	{
+		gl::GenFramebuffers(1, &id_);
+	}
+
+	framebuffer::framebuffer(GLuint id)
+		: id_(id)
+	{}
+
+	framebuffer::framebuffer(framebuffer && other)
+		: id_(other.id_)
+	{
+		other.id_ = 0;
+	}
+
+	framebuffer & framebuffer::operator = (framebuffer && other)
+	{
+		if (this == &other) return *this;
+
+		gl::DeleteFramebuffers(1, &id_);
+		id_ = other.id_;
+		other.id_ = 0;
+		return *this;
+	}
+
+	framebuffer::~framebuffer()
+	{
+		gl::DeleteFramebuffers(1, &id_);
+	}
+
+	framebuffer framebuffer::null()
+	{
+		return framebuffer{0};
+	}
+
+	void framebuffer::bind() const
+	{
+		gl::BindFramebuffer(gl::DRAW_FRAMEBUFFER, id_);
+	}
+
+	void framebuffer::color(texture_2d const & tex, int attachment)
+	{
+		bind();
+		gl::FramebufferTexture2D(gl::DRAW_FRAMEBUFFER, gl::COLOR_ATTACHMENT0 + attachment, gl::TEXTURE_2D, tex.id(), 0);
+	}
+
+	void framebuffer::color(renderbuffer const & rb, int attachment)
+	{
+		bind();
+		gl::FramebufferRenderbuffer(gl::DRAW_FRAMEBUFFER, gl::COLOR_ATTACHMENT0 + attachment, gl::RENDERBUFFER, rb.id());
+	}
+
+	void framebuffer::depth(renderbuffer const & rb)
+	{
+		bind();
+		gl::FramebufferRenderbuffer(gl::DRAW_FRAMEBUFFER, gl::DEPTH_STENCIL_ATTACHMENT, gl::RENDERBUFFER, rb.id());
+	}
+
+	GLenum framebuffer::status() const
+	{
+		bind();
+		return gl::CheckFramebufferStatus(gl::DRAW_FRAMEBUFFER);
+	}
+
+	bool framebuffer::complete() const
+	{
+		return status() == gl::FRAMEBUFFER_COMPLETE;
+	}
+
+	void framebuffer::assert_complete() const
+	{
+		if (auto s = status(); s != gl::FRAMEBUFFER_COMPLETE)
+			throw std::runtime_error("Framebuffer incomplete: " + framebuffer_status_string(s));
+	}
+
+}

libs/gfx/source/gfx/fullscreen.cpp

@@ -0,0 +1,39 @@
+#include <psemek/gfx/fullscreen.hpp>
+
+namespace psemek::gfx
+{
+
+	static indexed_mesh create_fullscreen_quad()
+	{
+		indexed_mesh m;
+
+		using vertex = gfx::vertex<geom::point<float, 2>>;
+
+		std::vector<vertex> vertices(4);
+		std::vector<std::uint8_t> indices(6);
+
+		vertices[0].get<0>() = {-1.f, -1.f};
+		vertices[1].get<0>() = { 1.f, -1.f};
+		vertices[2].get<0>() = {-1.f,  1.f};
+		vertices[3].get<0>() = { 1.f,  1.f};
+
+		indices[0] = 0;
+		indices[1] = 1;
+		indices[2] = 2;
+		indices[3] = 2;
+		indices[4] = 1;
+		indices[5] = 3;
+
+		m.setup<vertex>();
+		m.load(vertices, indices, gl::STATIC_DRAW);
+
+		return m;
+	}
+
+	indexed_mesh const & fullscreen_quad()
+	{
+		static indexed_mesh m = create_fullscreen_quad();
+		return m;
+	}
+
+}

libs/gfx/source/gfx/mesh.cpp

@@ -0,0 +1,129 @@
+#include <psemek/gfx/mesh.hpp>
+
+namespace psemek::gfx
+{
+
+	mesh mesh::null()
+	{
+		return mesh(0);
+	}
+
+	mesh::mesh()
+	{
+		gl::GenVertexArrays(1, &array_);
+		gl::GenBuffers(1, &buffer_);
+	}
+
+	mesh::mesh(mesh && other)
+	{
+		array_ = other.array_;
+		buffer_ = other.buffer_;
+		count_ = other.count_;
+
+		other.array_ = 0;
+		other.buffer_ = 0;
+		other.count_ = 0;
+	}
+
+	mesh & mesh::operator =(mesh && other)
+	{
+		if (this == &other) return *this;
+
+		if (array_)
+		{
+			gl::DeleteVertexArrays(1, &array_);
+			gl::DeleteBuffers(1, &buffer_);
+		}
+
+		array_ = other.array_;
+		buffer_ = other.buffer_;
+		count_ = other.count_;
+
+		other.array_ = 0;
+		other.buffer_ = 0;
+		other.count_ = 0;
+
+		return *this;
+	}
+
+	mesh::~mesh()
+	{
+		gl::DeleteVertexArrays(1, &array_);
+		gl::DeleteBuffers(1, &buffer_);
+	}
+
+	mesh::mesh(int)
+	{
+		array_ = 0;
+		buffer_ = 0;
+	}
+
+	indexed_mesh indexed_mesh::null()
+	{
+		return indexed_mesh(0);
+	}
+
+	indexed_mesh::indexed_mesh()
+	{
+		gl::GenVertexArrays(1, &array_);
+		gl::GenBuffers(1, &buffer_);
+		gl::GenBuffers(1, &index_buffer_);
+	}
+
+	indexed_mesh::indexed_mesh(indexed_mesh && other)
+	{
+		array_ = other.array_;
+		buffer_ = other.buffer_;
+		index_buffer_ = other.index_buffer_;
+		count_ = other.count_;
+		index_type_ = other.index_type_;
+
+		other.array_ = 0;
+		other.buffer_ = 0;
+		other.index_buffer_ = 0;
+		other.count_ = 0;
+		other.index_type_ = 0;
+	}
+
+	indexed_mesh & indexed_mesh::operator =(indexed_mesh && other)
+	{
+		if (this == &other) return *this;
+
+		if (array_)
+		{
+			gl::DeleteVertexArrays(1, &array_);
+			gl::DeleteBuffers(1, &buffer_);
+			gl::DeleteBuffers(1, &index_buffer_);
+		}
+
+		array_ = other.array_;
+		buffer_ = other.buffer_;
+		index_buffer_ = other.index_buffer_;
+		count_ = other.count_;
+		index_type_ = other.index_type_;
+
+		other.array_ = 0;
+		other.buffer_ = 0;
+		other.index_buffer_ = 0;
+		other.count_ = 0;
+		other.index_type_ = 0;
+
+		return *this;
+	}
+
+	indexed_mesh::~indexed_mesh()
+	{
+		gl::DeleteVertexArrays(1, &array_);
+		gl::DeleteBuffers(1, &buffer_);
+		gl::DeleteBuffers(1, &index_buffer_);
+	}
+
+	indexed_mesh::indexed_mesh(int)
+	{
+		array_ = 0;
+		buffer_ = 0;
+		index_buffer_ = 0;
+	}
+
+
+}

libs/gfx/source/gfx/program.cpp

@@ -0,0 +1,300 @@
+#include <psemek/gfx/program.hpp>
+
+#include <psemek/util/to_string.hpp>
+
+#include <memory>
+
+namespace psemek::gfx
+{
+
+	void program::uniform_proxy::operator = (int i)
+	{
+		gl::Uniform1i(location_, i);
+	}
+
+	void program::uniform_proxy::operator = (unsigned int u)
+	{
+		gl::Uniform1ui(location_, u);
+	}
+
+	void program::uniform_proxy::operator = (float f)
+	{
+		gl::Uniform1f(location_, f);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<int, 1> const & v)
+	{
+		gl::Uniform1i(location_, v[0]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<int, 2> const & v)
+	{
+		gl::Uniform2i(location_, v[0], v[1]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<int, 3> const & v)
+	{
+		gl::Uniform3i(location_, v[0], v[1], v[2]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<int, 4> const & v)
+	{
+		gl::Uniform4i(location_, v[0], v[1], v[2], v[3]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<unsigned int, 1> const & v)
+	{
+		gl::Uniform1ui(location_, v[0]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<unsigned int, 2> const & v)
+	{
+		gl::Uniform2ui(location_, v[0], v[1]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<unsigned int, 3> const & v)
+	{
+		gl::Uniform3ui(location_, v[0], v[1], v[2]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<unsigned int, 4> const & v)
+	{
+		gl::Uniform4ui(location_, v[0], v[1], v[2], v[3]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<float, 1> const & v)
+	{
+		gl::Uniform1f(location_, v[0]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<float, 2> const & v)
+	{
+		gl::Uniform2f(location_, v[0], v[1]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<float, 3> const & v)
+	{
+		gl::Uniform3f(location_, v[0], v[1], v[2]);
+	}
+
+	void program::uniform_proxy::operator = (geom::vector<float, 4> const & v)
+	{
+		gl::Uniform4f(location_, v[0], v[1], v[2], v[3]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<int, 1> const & v)
+	{
+		gl::Uniform1i(location_, v[0]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<int, 2> const & v)
+	{
+		gl::Uniform2i(location_, v[0], v[1]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<int, 3> const & v)
+	{
+		gl::Uniform3i(location_, v[0], v[1], v[2]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<int, 4> const & v)
+	{
+		gl::Uniform4i(location_, v[0], v[1], v[2], v[3]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<unsigned int, 1> const & v)
+	{
+		gl::Uniform1ui(location_, v[0]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<unsigned int, 2> const & v)
+	{
+		gl::Uniform2ui(location_, v[0], v[1]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<unsigned int, 3> const & v)
+	{
+		gl::Uniform3ui(location_, v[0], v[1], v[2]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<unsigned int, 4> const & v)
+	{
+		gl::Uniform4ui(location_, v[0], v[1], v[2], v[3]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<float, 1> const & v)
+	{
+		gl::Uniform1f(location_, v[0]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<float, 2> const & v)
+	{
+		gl::Uniform2f(location_, v[0], v[1]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<float, 3> const & v)
+	{
+		gl::Uniform3f(location_, v[0], v[1], v[2]);
+	}
+
+	void program::uniform_proxy::operator = (geom::point<float, 4> const & v)
+	{
+		gl::Uniform4f(location_, v[0], v[1], v[2], v[3]);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 2, 2> const & m)
+	{
+		gl::UniformMatrix2fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 2, 3> const & m)
+	{
+		gl::UniformMatrix2x3fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 2, 4> const & m)
+	{
+		gl::UniformMatrix2x4fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 3, 2> const & m)
+	{
+		gl::UniformMatrix3x2fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 3, 3> const & m)
+	{
+		gl::UniformMatrix3fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 3, 4> const & m)
+	{
+		gl::UniformMatrix3x4fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 4, 2> const & m)
+	{
+		gl::UniformMatrix4x2fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 4, 3> const & m)
+	{
+		gl::UniformMatrix4x3fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::matrix<float, 4, 4> const & m)
+	{
+		gl::UniformMatrix4fv(location_, 1, gl::TRUE_, m.coords);
+	}
+
+	void program::uniform_proxy::operator = (geom::interval<int> const & i)
+	{
+		gl::Uniform2i(location_, i.min, i.max);
+	}
+
+	void program::uniform_proxy::operator = (geom::interval<unsigned int> const & i)
+	{
+		gl::Uniform2ui(location_, i.min, i.max);
+	}
+
+	void program::uniform_proxy::operator = (geom::interval<float> const & i)
+	{
+		gl::Uniform2f(location_, i.min, i.max);
+	}
+
+	static void load_shader(GLuint shader, std::string_view source)
+	{
+		char const * vert_sources[1] { source.data() };
+		GLint vert_sources_len[1] { static_cast<GLint>(source.size()) };
+		gl::ShaderSource(shader, 1, vert_sources, vert_sources_len);
+		gl::CompileShader(shader);
+
+		GLint status;
+		gl::GetShaderiv(shader, gl::COMPILE_STATUS, &status);
+		if (status != gl::TRUE_)
+		{
+			GLint log_len;
+			gl::GetShaderiv(shader, gl::INFO_LOG_LENGTH, &log_len);
+			std::unique_ptr<char[]> log(new char [log_len]);
+			gl::GetShaderInfoLog(shader, log_len, nullptr, log.get());
+			throw std::runtime_error(util::to_string("Shader compilation failed: ", log.get()));
+		}
+	}
+
+	static void load_program(GLuint program, std::vector<GLuint> const & shaders)
+	{
+		for (auto s : shaders)
+			gl::AttachShader(program, s);
+		gl::LinkProgram(program);
+
+		GLint status;
+		gl::GetProgramiv(program, gl::LINK_STATUS, &status);
+		if (status != gl::TRUE_)
+		{
+			GLint log_len;
+			gl::GetProgramiv(program, gl::INFO_LOG_LENGTH, &log_len);
+			std::unique_ptr<char[]> log(new char [log_len]);
+			gl::GetProgramInfoLog(program, log_len, nullptr, log.get());
+			throw std::runtime_error(util::to_string("Program link failed: ", log.get()));
+		}
+
+		for (auto s : shaders)
+			gl::DetachShader(program, s);
+	}
+
+	static GLuint create_program(std::vector<std::pair<GLenum, std::string_view>> const & sources)
+	{
+		std::vector<GLuint> shaders;
+		for (auto const & p : sources)
+		{
+			GLuint sh = gl::CreateShader(p.first);
+			load_shader(sh, p.second);
+		}
+
+		GLuint program = gl::CreateProgram();
+		load_program(program, shaders);
+
+		for (auto s : shaders)
+			gl::DeleteShader(s);
+
+		return program;
+	}
+
+	program::program(std::string_view vertex_source, std::string_view fragment_source)
+	{
+		program_ = create_program({{gl::VERTEX_SHADER, vertex_source}, {gl::FRAGMENT_SHADER, fragment_source}});
+	}
+
+	program::program(std::string_view vertex_source, std::string_view geometry_source, std::string_view fragment_source)
+	{
+		program_ = create_program({{gl::VERTEX_SHADER, vertex_source}, {gl::GEOMETRY_SHADER, geometry_source}, {gl::FRAGMENT_SHADER, fragment_source}});
+	}
+
+	GLuint program::id() const
+	{
+		return program_;
+	}
+
+	void program::bind() const
+	{
+		gl::UseProgram(program_);
+	}
+
+	GLint program::location(char const * name) const
+	{
+		auto it = uniforms_.find(name);
+		if (it == uniforms_.end())
+		{
+			auto l = gl::GetUniformLocation(program_, name);
+			uniforms_[name] = l;
+			return l;
+		}
+		return it->second;
+	}
+
+	program::uniform_proxy program::operator[] (char const * name) const
+	{
+		return {location(name)};
+	}
+
+}

libs/gfx/source/gfx/renderbuffer.cpp

@@ -0,0 +1,52 @@
+#include <psemek/gfx/renderbuffer.hpp>
+
+namespace psemek::gfx
+{
+
+	renderbuffer::renderbuffer()
+	{
+		gl::GenRenderbuffers(1, &id_);
+	}
+
+	renderbuffer::renderbuffer(renderbuffer && other)
+	{
+		id_ = other.id_;
+		other.id_ = 0;
+	}
+
+	renderbuffer & renderbuffer::operator = (renderbuffer && other)
+	{
+		if (this == &other) return *this;
+
+		gl::DeleteRenderbuffers(1, &id_);
+		id_ = other.id_;
+		other.id_ = 0;
+		return *this;
+	}
+
+	renderbuffer::~renderbuffer()
+	{
+		gl::DeleteRenderbuffers(1, &id_);
+	}
+
+	renderbuffer renderbuffer::null()
+	{
+		return renderbuffer(0);
+	}
+
+	void renderbuffer::bind() const
+	{
+		gl::BindRenderbuffer(gl::RENDERBUFFER, id_);
+	}
+
+	void renderbuffer::storage(GLenum internal_format, GLsizei width, GLsizei height)
+	{
+		bind();
+		gl::RenderbufferStorage(gl::RENDERBUFFER, internal_format, width, height);
+	}
+
+	renderbuffer::renderbuffer(GLuint id)
+		: id_(id)
+	{}
+
+}

libs/gfx/source/gfx/texture.cpp

@@ -0,0 +1,84 @@
+#include <psemek/gfx/texture.hpp>
+
+namespace psemek::gfx
+{
+
+	texture_2d::texture_2d()
+	{
+		gl::GenTextures(1, &id_);
+	}
+
+	texture_2d::texture_2d(GLuint id)
+		: id_(id)
+	{}
+
+	texture_2d texture_2d::null()
+	{
+		return texture_2d(0);
+	}
+
+	void texture_2d::bind() const
+	{
+		gl::BindTexture(gl::TEXTURE_2D, id_);
+	}
+
+	texture_2d::texture_2d(texture_2d && other)
+		: id_(other.id_)
+		, width_(other.width_)
+		, height_(other.height_)
+	{
+		other.id_ = 0;
+		other.width_ = 0;
+		other.height_ = 0;
+	}
+
+	texture_2d & texture_2d::operator = (texture_2d && other)
+	{
+		if (this == &other) return *this;
+
+		gl::DeleteTextures(1, &id_);
+		id_ = other.id_;
+		width_ = other.width_;
+		height_ = other.height_;
+		other.id_ = 0;
+		other.width_ = 0;
+		other.height_ = 0;
+
+		return *this;
+	}
+
+	texture_2d::~texture_2d()
+	{
+		gl::DeleteTextures(1, &id_);
+	}
+
+	void texture_2d::load(GLint internal_format, std::size_t width, std::size_t height, GLenum format, GLenum type, const void * data)
+	{
+		bind();
+		gl::TexImage2D(gl::TEXTURE_2D, 0, internal_format, width, height, 0, format, type, data);
+		gl::GenerateMipmap(gl::TEXTURE_2D);
+
+		width_ = width;
+		height_ = height;
+	}
+
+	void texture_2d::pixels(GLenum format, GLenum type, void * data) const
+	{
+		bind();
+		gl::GetTexImage(gl::TEXTURE_2D, 0, format, type, data);
+	}
+
+	texture_2d texture_2d::from_data(GLint internal_format, std::size_t width, std::size_t height, GLenum format, GLenum type, const void * data)
+	{
+		texture_2d tex;
+		tex.load(internal_format, width, height, format, type, data);
+		return tex;
+	}
+
+	void texture_2d::generate_mipmap()
+	{
+		bind();
+		gl::GenerateMipmap(gl::TEXTURE_2D);
+	}
+
+}

libs/pcg/CMakeLists.txt

@@ -0,0 +1,6 @@
+file(GLOB_RECURSE PSEMEK_PCG_HEADERS RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "include/*.hpp")
+file(GLOB_RECURSE PSEMEK_PCG_SOURCES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "source/*.cpp")
+
+add_library(pcg ${PSEMEK_PCG_HEADERS} ${PSEMEK_PCG_SOURCES})
+target_include_directories(pcg PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
+target_link_libraries(pcg PUBLIC util geom)

libs/pcg/include/psemek/pcg/perlin.hpp

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <psemek/util/pixmap.hpp>
+#include <psemek/geom/vector.hpp>
+
+namespace psemek::pcg
+{
+
+	struct perlin
+	{
+		perlin() = default;
+		perlin(util::basic_pixmap<geom::vector<float, 2>> grad_map);
+		perlin(perlin &&) = default;
+
+		perlin & operator= (perlin &&) = default;
+
+		std::size_t width() const
+		{
+			return grad_map_.width() - 1;
+		}
+
+		std::size_t height() const
+		{
+			return grad_map_.height() - 1;
+		}
+
+		// x \in [0.0 .. 1.0]
+		// y \in [0.0 .. 1.0]
+		float operator() (float x, float y) const;
+
+	private:
+		util::basic_pixmap<geom::vector<float, 2>> grad_map_;
+	};
+
+}

libs/pcg/source/perlin.cpp

@@ -0,0 +1,50 @@
+#include <psemek/pcg/perlin.hpp>
+
+#include <psemek/geom/interval.hpp>
+
+#include <psemek/util/assert.hpp>
+
+namespace psemek::pcg
+{
+
+	perlin::perlin(util::basic_pixmap<geom::vector<float, 2>> grad_map)
+		: grad_map_(std::move(grad_map))
+	{}
+
+	static float step(float x0, float x1, float t)
+	{
+		return x0 * (1.f - t) + x1 * t;
+	}
+
+	static float smoothstep(float x0, float x1, float t)
+	{
+		float const s = t * t * (3.f - 2.f * t);
+		return step(x0, x1, s);
+	}
+
+	float perlin::operator() (float x, float y) const
+	{
+		assert(x >= 0.f);
+		assert(y >= 0.f);
+
+		assert(x <= 1.f);
+		assert(y <= 1.f);
+
+		x *= width();
+		y *= height();
+
+		int const ix = geom::clamp<int>(std::floor(x), {0, static_cast<int>(width()) - 1});
+		int const iy = geom::clamp<int>(std::floor(y), {0, static_cast<int>(height()) - 1});
+
+		float const tx = x - ix;
+		float const ty = y - iy;
+
+		float const d00 = tx * grad_map_(ix, iy)[0] + ty * grad_map_(ix, iy)[1];
+		float const d10 = (tx-1.f) * grad_map_(ix+1, iy)[0] + ty * grad_map_(ix+1, iy)[1];
+		float const d01 = tx * grad_map_(ix, iy+1)[0] + (ty-1.f) * grad_map_(ix, iy+1)[1];
+		float const d11 = (tx-1.f) * grad_map_(ix+1, iy+1)[0] + (ty-1.f) * grad_map_(ix+1, iy+1)[1];
+
+		return smoothstep(smoothstep(d00, d10, tx), smoothstep(d01, d11, tx), ty);
+	}
+
+}

libs/util/CMakeLists.txt

@@ -0,0 +1,8 @@
+find_package(Threads)
+
+file(GLOB_RECURSE PSEMEK_UTIL_HEADERS RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "include/*.hpp")
+file(GLOB_RECURSE PSEMEK_UTIL_SOURCES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "source/*.cpp")
+
+add_library(util ${PSEMEK_UTIL_HEADERS} ${PSEMEK_UTIL_SOURCES})
+target_include_directories(util PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include")
+target_link_libraries(util PUBLIC ${CMAKE_THREAD_LIBS_INIT})

libs/util/include/psemek/util/assert.hpp

@@ -0,0 +1,18 @@
+#pragma once
+
+#include <stdexcept>
+
+#undef assert
+
+#ifdef PSEMEK_DEBUG
+#define assert(x) ((void)(!(x) && ::psemek::util::assert_handler(#x, __FILE__, __LINE__)))
+#else
+#define assert(x) ((void)sizeof(x))
+#endif
+
+namespace psemek::util
+{
+
+	[[noreturn]] bool assert_handler (char const * expression, char const * file, int line);
+
+}

libs/util/include/psemek/util/at_scope_exit.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+namespace psemek::util
+{
+
+	template <typename F>
+	struct at_scope_exit
+	{
+		F f;
+
+		at_scope_exit(F f)
+			: f(f)
+		{}
+
+		~at_scope_exit()
+		{
+			f();
+		}
+	};
+
+}

libs/util/include/psemek/util/autoname.hpp

@@ -0,0 +1,5 @@
+#pragma once
+
+#define _UTIL_CAT2(x,y) x##y
+#define _UTIL_CAT(x,y) _UTIL_CAT2(x,y)
+#define autoname _UTIL_CAT(autonamed_,__COUNTER__)

libs/util/include/psemek/util/blob.hpp

@@ -0,0 +1,121 @@
+#pragma once
+
+#include <memory>
+#include <string>
+#include <string_view>
+
+namespace psemek::util
+{
+
+	struct blob
+	{
+		blob() = default;
+		blob(blob const & other);
+		blob(blob && other);
+
+		blob(std::size_t size);
+		blob(std::size_t size, std::unique_ptr<char[]> data);
+		blob(std::size_t size, char * data);
+
+		blob & operator = (blob const & other);
+		blob & operator = (blob && other);
+
+		~ blob() = default;
+
+		char * data() { return data_.get(); }
+		char const * data() const { return data_.get(); }
+
+		std::size_t size() const { return size_; }
+
+		void reset();
+		std::unique_ptr<char[]> release();
+
+		using iterator = char *;
+		using const_iterator = char const *;
+
+		char * begin() { return data(); }
+		char const * begin() const { return data(); }
+
+		char * end() { return data() + size(); }
+		char const * end() const { return data() + size(); }
+
+		char & operator[] (std::size_t i) { return data()[i]; }
+		char const & operator[] (std::size_t i) const { return data()[i]; }
+
+		std::string string() const;
+		std::string_view string_view() const;
+
+	private:
+		std::unique_ptr<char[]> data_;
+		std::size_t size_ = 0;
+	};
+
+	inline blob::blob(blob const & other)
+	{
+		*this = other;
+	}
+
+	inline blob::blob(blob && other)
+	{
+		*this = std::move(other);
+	}
+
+	inline blob::blob(std::size_t size)
+		: data_{new char [size]}
+		, size_{size}
+	{}
+
+	inline blob::blob(std::size_t size, std::unique_ptr<char[]> data)
+		: data_{std::move(data)}
+		, size_{size}
+	{}
+
+	inline blob::blob(std::size_t size, char * data)
+		: data_{data}
+		, size_{size}
+	{}
+
+	inline blob & blob::operator=(blob const & other)
+	{
+		if (this != &other)
+		{
+			data_.reset(new char [other.size()]);
+			std::copy(other.begin(), other.end(), begin());
+			size_ = other.size();
+		}
+		return *this;
+	}
+
+	inline blob & blob::operator=(blob && other)
+	{
+		if (this != &other)
+		{
+			data_ = std::move(other.data_);
+			size_ = other.size();
+			other.size_ = 0;
+		}
+		return *this;
+	}
+
+	inline void blob::reset()
+	{
+		data_.reset();
+		size_ = 0;
+	}
+
+	inline std::unique_ptr<char[]> blob::release()
+	{
+		size_ = 0;
+		return std::move(data_);
+	}
+
+	inline std::string blob::string() const
+	{
+		return std::string(data_.get(), data_.get() + size_);
+	}
+
+	inline std::string_view blob::string_view() const
+	{
+		return std::string_view(data_.get(), size_);
+	}
+}

libs/util/include/psemek/util/clock.hpp

@@ -0,0 +1,49 @@
+#pragma once
+
+#include <chrono>
+
+namespace psemek::util
+{
+
+	template <typename Duration = std::chrono::duration<double>, typename Clock = std::chrono::system_clock>
+	struct clock
+	{
+		typedef Duration duration_type;
+		typedef typename duration_type::rep rep_type;
+
+		typedef Clock clock_type;
+		typedef typename clock_type::time_point time_point_type;
+
+		clock ( )
+		{
+			restart();
+		}
+
+		time_point_type now ( ) const
+		{
+			return clock_type::now();
+		}
+
+		duration_type restart ( )
+		{
+			auto t = now();
+			auto delta = t - start_;
+			start_ = t;
+			return std::chrono::duration_cast<duration_type>(delta);
+		}
+
+		duration_type duration ( ) const
+		{
+			return std::chrono::duration_cast<duration_type>(now() - start_);
+		}
+
+		rep_type count ( ) const
+		{
+			return duration().count();
+		}
+
+	private:
+		time_point_type start_;
+	};
+
+}

libs/util/include/psemek/util/flat_list.hpp

@@ -0,0 +1,77 @@
+#pragma once
+
+#include <vector>
+#include <type_traits>
+#include <initializer_list>
+
+namespace psemek::util
+{
+
+	template <typename T, typename IndexType = std::size_t>
+	struct flat_list
+	{
+		flat_list ( );
+		explicit flat_list (std::size_t count);
+		flat_list (std::size_t count, T const & value);
+		flat_list (flat_list const &);
+		flat_list (flat_list &&);
+		flat_list (std::initializer_list<T> init);
+
+		~ flat_list ( );
+
+		std::size_t size ( ) const;
+
+	private:
+		using item = std::aligned_storage_t<std::max(sizeof(T), sizeof(IndexType))>;
+
+		static constexpr std::size_t nil = static_cast<std::size_t>(-1);
+
+		// Invariant: data_.size() >= size_
+		std::vector<item> data_;
+		std::size_t size_ = 0;
+		std::size_t first_free_ = nil;
+	};
+
+	template <typename T, typename I>
+	flat_list<T, I>::flat_list ( ) = default;
+
+	template <typename T, typename I>
+	flat_list<T, I>::flat_list (std::size_t count)
+		: data_(count)
+		, size_(count)
+	{ }
+
+	template <typename T, typename I>
+	flat_list<T, I>::flat_list (std::size_t count, T const & value)
+		: data_(count, value)
+		, size_(count)
+	{ }
+
+	template <typename T, typename I>
+	flat_list<T, I>::flat_list (flat_list const & other) = default;
+
+	template <typename T, typename I>
+	flat_list<T, I>::flat_list (flat_list && other)
+		: data_(std::move(other.data_))
+		, size_(other.size_)
+		, first_free_(other.first_free_)
+	{
+		other.size_ = 0;
+		other.first_free_ = nil;
+	}
+
+	template <typename T, typename I>
+	flat_list<T, I>::flat_list (std::initializer_list<T> init)
+		: data_(std::move(init))
+		, size_(data_.size())
+	{ }
+
+	template <typename T, typename I>
+	flat_list<T, I>::~flat_list ( ) = default;
+
+	template <typename T, typename I>
+	std::size_t flat_list<T, I>::size ( ) const
+	{
+		return size_;
+	}
+}

libs/util/include/psemek/util/fmap.hpp

@@ -0,0 +1,38 @@
+#pragma once
+
+#include <psemek/util/range.hpp>
+
+#include <optional>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		template <typename F>
+		struct fmap
+		{
+			F f;
+
+			template <typename T>
+			auto operator() (std::optional<T> && x)
+			{
+				using R = decltype(f(*x));
+
+				if (x)
+					return std::optional<R>(f(*x));
+				else
+					return std::optional<R>();
+			}
+		};
+
+	}
+
+	template <typename F>
+	auto fmap (F f)
+	{
+		return detail::fmap<F>{std::move(f)};
+	}
+
+}

libs/util/include/psemek/util/functional.hpp

@@ -0,0 +1,24 @@
+#pragma once
+
+namespace psemek::util
+{
+
+	constexpr auto nop = [](auto const & ...){};
+
+	constexpr auto id = [](auto && x) -> decltype(auto) { return x; };
+
+	template <typename T>
+	auto constant (T const & x)
+	{
+		return [x](auto const & ...){ return x; };
+	}
+
+	template <typename F1, typename F2>
+	auto bind_and (F1 && f1, F2 && f2)
+	{
+		return [=](auto const &... args){
+			return f1(args...) && f2(args...);
+		};
+	}
+
+}

libs/util/include/psemek/util/lazy_range.hpp

@@ -0,0 +1,65 @@
+#pragma once
+
+#include <iterator>
+
+namespace psemek::util
+{
+
+	template <typename Gen>
+	struct lazy_range
+	{
+		Gen generator;
+
+		lazy_range (Gen gen)
+			: generator(std::move(gen))
+		{ }
+
+		using value_type = decltype(generator());
+
+		struct iterator
+		{
+			using iterator_category = std::input_iterator_tag;
+
+			lazy_range & range;
+
+			value_type value;
+
+			value_type const & operator * ()
+			{
+				return value;
+			}
+
+			iterator & operator ++ ()
+			{
+				value = range.generator();
+				return *this;
+			}
+		};
+
+		iterator begin ( )
+		{
+			return { *this, generator() };
+		}
+
+		struct sentinel
+		{
+			using iterator_category = std::input_iterator_tag;
+		};
+
+		sentinel end ( )
+		{
+			return {};
+		}
+
+		friend bool operator == (iterator const &, sentinel)
+		{
+			return false;
+		}
+
+		friend bool operator != (iterator const &, sentinel)
+		{
+			return true;
+		}
+	};
+
+}

libs/util/include/psemek/util/movable_function.hpp

@@ -0,0 +1,84 @@
+#pragma once
+
+#include <memory>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		template <typename Signature>
+		struct movable_function_node_base;
+
+		template <typename R, typename ... Args>
+		struct movable_function_node_base<R(Args...)>
+		{
+			virtual R call (Args ...) = 0;
+
+			virtual ~ movable_function_node_base ( ) = default;
+		};
+
+		template <typename Signature, typename F>
+		struct movable_function_node;
+
+		template <typename R, typename ... Args, typename F>
+		struct movable_function_node<R(Args...), F>
+			: movable_function_node_base<R(Args...)>
+		{
+			F f;
+
+			movable_function_node (F && f)
+				: f(std::move(f))
+			{ }
+
+			R call (Args ... args) override
+			{
+				return f(args...);
+			}
+		};
+
+		// Implemented in cpp to prevent dependency on <functional>
+		[[noreturn]] void bad_function_call ( );
+
+	}
+
+	template <typename Signature>
+	struct movable_function;
+
+	template <typename R, typename ... Args>
+	struct movable_function<R(Args...)>
+	{
+		using signature = R(Args...);
+
+		movable_function ( ) = default;
+
+		template <typename F>
+		movable_function (F f)
+			: p { std::make_unique<detail::movable_function_node<signature, std::remove_cv_t<F>>>(std::move(f)) }
+		{ }
+
+		movable_function (movable_function &&) = default;
+		movable_function & operator = (movable_function &&) = default;
+
+		movable_function (movable_function const &) = delete;
+		movable_function & operator = (movable_function const &) = delete;
+
+		explicit operator bool ( ) const
+		{
+			return static_cast<bool>(p);
+		}
+
+		R operator() (Args ... args) const
+		{
+			if (!p)
+				detail::bad_function_call();
+
+			return p->call(args...);
+		}
+
+	private:
+		std::unique_ptr<detail::movable_function_node_base<signature>> p;
+	};
+
+}

libs/util/include/psemek/util/moving_average.hpp

@@ -0,0 +1,62 @@
+#pragma once
+
+#include <vector>
+
+namespace psemek::util
+{
+
+	template <typename T>
+	struct moving_average
+	{
+		moving_average (std::size_t max)
+			: data_(max, T())
+			, begin_(0)
+			, size_(0)
+			, sum_(T())
+		{ }
+
+		void clear ( )
+		{
+			begin_ = 0;
+			size_ = 0;
+			sum_ = T(0);
+		}
+
+		void push (T x)
+		{
+			if (size_ >= data_.size())
+			{
+				sum_ -= data_[begin_];
+			}
+			else
+				size_++;
+
+			data_[begin_++] = x;
+			sum_ += x;
+
+			if (begin_ == data_.size())
+				begin_ = 0;
+		}
+
+		T sum ( ) const
+		{
+			return sum_;
+		}
+
+		std::size_t count ( ) const
+		{
+			return size_;
+		}
+
+		T average ( ) const
+		{
+			return sum() / count();
+		}
+
+	private:
+		std::vector<T> data_;
+		std::size_t begin_, size_;
+		T sum_;
+	};
+
+}

libs/util/include/psemek/util/multidimensional_array.hpp

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <array>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+		template <typename Type, std::size_t ... Sizes>
+		struct multidimentional_array_impl;
+
+		template <typename Type, std::size_t FirstSize, std::size_t ... Sizes>
+		struct multidimentional_array_impl<Type, FirstSize, Sizes...>
+		{
+			typedef std::array<typename multidimentional_array_impl<Type, Sizes...>::type, FirstSize> type;
+		};
+
+		template <typename Type>
+		struct multidimentional_array_impl<Type>
+		{
+			typedef Type type;
+		};
+	}
+
+	template <typename Type, std::size_t ... Sizes>
+	using multidimentional_array = typename detail::multidimentional_array_impl<Type, Sizes...>::type;
+
+}

libs/util/include/psemek/util/noncopyable.hpp

@@ -0,0 +1,16 @@
+#pragma once
+
+namespace psemek::util
+{
+
+	struct noncopyable
+	{
+		noncopyable (noncopyable const &)              = delete;
+		noncopyable & operator = (noncopyable const &) = delete;
+
+		noncopyable ( )                           = default;
+		noncopyable (noncopyable &&)              = default;
+		noncopyable & operator = (noncopyable &&) = default;
+	};
+
+}

libs/util/include/psemek/util/not_implemented.hpp

@@ -0,0 +1,8 @@
+#pragma once
+
+namespace psemek::util
+{
+
+	void not_implemented ( );
+
+}

libs/util/include/psemek/util/overload.hpp

@@ -0,0 +1,30 @@
+#pragma once
+
+#include <utility>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		template <typename ... Fs>
+		struct overload_impl
+			: Fs...
+		{
+			overload_impl (Fs ... fs)
+				: Fs(fs)...
+			{ }
+
+			using Fs::operator()...;
+		};
+
+	}
+
+	template <typename ... Fs>
+	auto overload (Fs ... fs)
+	{
+		return detail::overload_impl<Fs...>{std::move(fs)...};
+	}
+
+}

libs/util/include/psemek/util/pimpl.hpp

@@ -0,0 +1,68 @@
+#pragma once
+
+#include <type_traits>
+#include <memory>
+
+namespace psemek::util::pimpl
+{
+
+	template <typename Parent>
+	struct impl;
+
+	template <typename Parent, std::size_t Size>
+	struct in_place
+	{
+	protected:
+
+		using impl_type = ::util::pimpl::impl<Parent>;
+
+		impl_type * pimpl ( ) { return reinterpret_cast<impl_type *>(std::addressof(storage_)); }
+		impl_type & impl  ( ) { return *pimpl(); }
+
+		impl_type const * pimpl ( ) const { return reinterpret_cast<impl_type const *>(std::addressof(storage_)); }
+		impl_type const & impl  ( ) const { return *pimpl(); }
+
+		template <typename ... Args>
+		in_place (Args && ... args)
+		{
+			static_assert(sizeof(impl_type) <= Size, "impl storage size too small");
+			new (pimpl()) impl_type (std::forward<Args>(args)...);
+		}
+
+		~ in_place ( )
+		{
+			impl().~impl_type();
+		}
+
+	private:
+
+		typename std::aligned_storage<Size>::type storage_;
+	};
+
+	template <typename Parent>
+	struct dynamic
+	{
+	protected:
+
+		using impl_type = ::util::pimpl::impl<Parent>;
+
+		impl_type * pimpl ( ) { return pointer_.get(); }
+		impl_type & impl  ( ) { return *pimpl(); }
+
+		impl_type const * pimpl ( ) const { return pointer_.get(); }
+		impl_type const & impl  ( ) const { return *pimpl(); }
+
+		template <typename ... Args>
+		dynamic (Args && ... args)
+		{
+			pointer_.reset(new impl_type(std::forward<Args>(args)...));
+		}
+
+		~ dynamic ( ) = default;
+
+	private:
+
+		std::unique_ptr<impl_type> pointer_;
+	};
+
+}

libs/util/include/psemek/util/pixmap.hpp

@@ -0,0 +1,201 @@
+#pragma once
+
+#include <memory>
+#include <cstdint>
+
+namespace psemek::util
+{
+
+	template <typename Pixel>
+	struct basic_pixmap
+	{
+		using pixel_type = Pixel;
+
+		basic_pixmap() = default;
+		basic_pixmap(std::size_t width, std::size_t height);
+		basic_pixmap(std::size_t width, std::size_t height, Pixel value);
+		basic_pixmap(std::size_t width, std::size_t height, std::unique_ptr<Pixel[]> data);
+		basic_pixmap(basic_pixmap &&);
+
+		basic_pixmap(basic_pixmap const &) = delete;
+
+		basic_pixmap & operator = (basic_pixmap &&);
+		basic_pixmap & operator = (basic_pixmap const &) = delete;
+
+		std::size_t width() const { return width_; }
+		std::size_t height() const { return height_; }
+
+		void resize(std::size_t width, std::size_t height);
+		void resize(std::size_t width, std::size_t height, Pixel value);
+
+		Pixel * operator[] (std::size_t row);
+		Pixel const * operator[] (std::size_t row) const;
+
+		Pixel * data() { return data_.get(); }
+		Pixel const * data() const { return data_.get(); }
+
+		Pixel & operator()(std::size_t column, std::size_t row);
+		Pixel const & operator()(std::size_t column, std::size_t row) const;
+
+		std::unique_ptr<Pixel[]> release();
+		void reset(std::size_t width, std::size_t height, std::unique_ptr<Pixel[]> data);
+
+		bool empty() const;
+
+		void clear();
+
+		void fill(Pixel value);
+
+	private:
+		std::size_t width_ = 0;
+		std::size_t height_ = 0;
+		std::unique_ptr<Pixel[]> data_;
+	};
+
+	using pixmap_monochrome = basic_pixmap<std::uint8_t>;
+	using pixmap_rgb = basic_pixmap<std::array<std::uint8_t, 3>>;
+	using pixmap_rgba = basic_pixmap<std::array<std::uint8_t, 4>>;
+	using pixmap_float = basic_pixmap<float>;
+
+	pixmap_monochrome read_pgm(std::istream & is);
+	pixmap_rgb read_ppm(std::istream & is);
+
+	void write_pgm(pixmap_monochrome const & p, std::ostream & os);
+	void write_ppm(pixmap_rgb const & p, std::ostream & os);
+
+	template <typename Pixel>
+	basic_pixmap<Pixel>::basic_pixmap(std::size_t width, std::size_t height)
+		: width_(width)
+		, height_(height)
+	{
+		data_.reset(new Pixel [width_ * height_]);
+	}
+
+	template <typename Pixel>
+	basic_pixmap<Pixel>::basic_pixmap(std::size_t width, std::size_t height, Pixel value)
+		: width_(width)
+		, height_(height)
+	{
+		data_.reset(new Pixel [width_ * height_]);
+		fill(value);
+	}
+
+	template <typename Pixel>
+	basic_pixmap<Pixel>::basic_pixmap(std::size_t width, std::size_t height, std::unique_ptr<Pixel[]> data)
+		: width_(width)
+		, height_(height)
+		, data_(std::move(data))
+	{}
+
+	template <typename Pixel>
+	basic_pixmap<Pixel>::basic_pixmap(basic_pixmap && other)
+		: width_(other.width_)
+		, height_(other.height_)
+		, data_(other.release())
+	{}
+
+	template <typename Pixel>
+	basic_pixmap<Pixel> & basic_pixmap<Pixel>::operator = (basic_pixmap && other)
+	{
+		if (this == &other) return *this;
+
+		width_ = other.width_;
+		height_ = other.height_;
+		data_ = other.release();
+
+		return *this;
+	}
+
+	template <typename Pixel>
+	void basic_pixmap<Pixel>::resize(std::size_t width, std::size_t height)
+	{
+		data_.reset(new Pixel[width * height]);
+		width_ = width;
+		height_ = height;
+	}
+
+	template <typename Pixel>
+	void basic_pixmap<Pixel>::resize(std::size_t width, std::size_t height, Pixel value)
+	{
+		data_.reset(new Pixel[width * height]);
+		width_ = width;
+		height_ = height;
+		fill(value);
+	}
+
+	template <typename Pixel>
+	Pixel * basic_pixmap<Pixel>::operator[] (std::size_t row)
+	{
+		return data() + row * width();
+	}
+
+	template <typename Pixel>
+	Pixel const * basic_pixmap<Pixel>::operator[] (std::size_t row) const
+	{
+		return data() + row * width();
+	}
+
+	template <typename Pixel>
+	Pixel & basic_pixmap<Pixel>::operator()(std::size_t column, std::size_t row)
+	{
+		return (*this)[row][column];
+	}
+
+	template <typename Pixel>
+	Pixel const & basic_pixmap<Pixel>::operator()(std::size_t column, std::size_t row) const
+	{
+		return (*this)[row][column];
+	}
+
+	template <typename Pixel>
+	std::unique_ptr<Pixel[]> basic_pixmap<Pixel>::release()
+	{
+		auto data = std::move(data_);
+		width_ = 0;
+		height_ = 0;
+		return data;
+	}
+
+	template <typename Pixel>
+	void basic_pixmap<Pixel>::reset(std::size_t width, std::size_t height, std::unique_ptr<Pixel[]> data)
+	{
+		width_ = width;
+		height_ = height;
+		data_ = std::move(data);
+	}
+
+	template <typename Pixel>
+	bool basic_pixmap<Pixel>::empty() const
+	{
+		return width() == 0 || height() == 0;
+	}
+
+	template <typename Pixel>
+	void basic_pixmap<Pixel>::clear()
+	{
+		release();
+	}
+
+	template <typename Pixel>
+	void basic_pixmap<Pixel>::fill(Pixel value)
+	{
+		std::fill(data_.get(), data_.get() + width_ * height_, value);
+	}
+
+	template <typename Pixel>
+	void mirror_x(basic_pixmap<Pixel> & p)
+	{
+		for (std::size_t y = 0; y < p.height(); ++y)
+			for (std::size_t x = 0; x < p.width() / 2; ++x)
+				std::swap(p(x, y), p(p.width() - x - 1, y));
+	}
+
+	template <typename Pixel>
+	void mirror_y(basic_pixmap<Pixel> & p)
+	{
+		for (std::size_t y = 0; y < p.height() / 2; ++y)
+			for (std::size_t x = 0; x < p.width(); ++x)
+				std::swap(p(x, y), p(x, p.height() - y - 1));
+	}
+
+}

libs/util/include/psemek/util/pretty_print.hpp

@@ -0,0 +1,44 @@
+#pragma once
+
+#include <chrono>
+#include <iostream>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		template <typename Duration, typename UpTo>
+		struct pretty_print_time_wrapper
+		{
+			Duration d;
+			UpTo up_to;
+		};
+
+		void pretty_print_time (std::ostream & o, std::int64_t d, std::int64_t up_to);
+
+		template <typename Duration, typename UpTo>
+		std::ostream & operator << (std::ostream & o, pretty_print_time_wrapper<Duration, UpTo> t)
+		{
+			std::int64_t const d = std::chrono::duration_cast<std::chrono::nanoseconds>(t.d).count();
+			std::int64_t const up_to = std::chrono::duration_cast<std::chrono::nanoseconds>(t.up_to).count();
+			pretty_print_time(o, d, up_to);
+			return o;
+		}
+
+	}
+
+	template <typename Rep, typename Period, typename UpTo>
+	auto pretty (std::chrono::duration<Rep, Period> d, UpTo up_to)
+	{
+		return detail::pretty_print_time_wrapper<std::chrono::duration<Rep, Period>, UpTo>{d, up_to};
+	}
+
+	template <typename Rep, typename Period>
+	auto pretty (std::chrono::duration<Rep, Period> d)
+	{
+		return pretty(d, std::chrono::seconds{1});
+	}
+
+}

libs/util/include/psemek/util/profiler.hpp

@@ -0,0 +1,27 @@
+#pragma once
+
+#include <psemek/util/clock.hpp>
+#include <psemek/util/pretty_print.hpp>
+
+#include <string>
+
+namespace psemek::util
+{
+
+	struct profiler
+	{
+		profiler (std::string name)
+			: name_(std::move(name))
+		{ }
+
+		~ profiler ( )
+		{
+			std::cout << name_ << ": " << util::pretty(clock_.duration(), std::chrono::microseconds{1}) << "\n";
+		}
+
+	private:
+		std::string const name_;
+		util::clock<> clock_;
+	};
+
+}

libs/util/include/psemek/util/range.hpp

@@ -0,0 +1,80 @@
+#pragma once
+
+#include <iterator>
+#include <utility>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		template <typename Container>
+		auto begin_helper (Container & x)
+		{
+			using std::begin;
+			return begin(x);
+		}
+
+		template <typename Container>
+		auto end_helper (Container & x)
+		{
+			using std::end;
+			return end(x);
+		}
+
+	}
+
+	template <typename Container>
+	auto begin (Container & x)
+	{
+		return detail::begin_helper(x);
+	}
+
+	template <typename Container>
+	auto end (Container & x)
+	{
+		return detail::end_helper(x);
+	}
+
+	template <typename T>
+	struct range_traits
+	{
+		using iterator = decltype(begin(std::declval<T>()));
+		using iterator_traits = std::iterator_traits<iterator>;
+
+		using iterator_category = typename iterator_traits::iterator_category;
+		using value_type        = typename iterator_traits::value_type       ;
+		using difference        = typename iterator_traits::difference       ;
+		using pointer           = typename iterator_traits::pointer          ;
+		using reference         = typename iterator_traits::reference        ;
+	};
+
+	template <typename Iterator>
+	struct range
+	{
+		Iterator it_begin;
+		Iterator it_end;
+
+		Iterator begin() const
+		{
+			return it_begin;
+		}
+
+		Iterator end() const
+		{
+			return it_end;
+		}
+	};
+
+	template <typename Range>
+	auto reversed (Range const & r)
+	{
+		auto it1 = begin(r);
+		auto it2 = end(r);
+
+		using ReverseIterator = std::reverse_iterator<decltype(it1)>;
+		return range<ReverseIterator>{ std::make_reverse_iterator(it2), std::make_reverse_iterator(it1) };
+	}
+
+}

libs/util/include/psemek/util/recursive.hpp

@@ -0,0 +1,29 @@
+#pragma once
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		template <typename F>
+		struct recursive_impl
+		{
+			F f;
+
+			template <typename ... Args>
+			auto operator() (Args && ... args) -> decltype(auto)
+			{
+				return f(*this, std::forward<Args>(args)...);
+			}
+		};
+
+	}
+
+	template <typename F>
+	auto recursive (F f)
+	{
+		return detail::recursive_impl<F>{std::move(f)};
+	}
+
+}

libs/util/include/psemek/util/shared_blob.hpp

@@ -0,0 +1,129 @@
+#pragma once
+
+#include <psemek/util/blob.hpp>
+
+#include <memory>
+#include <string>
+#include <string_view>
+
+namespace psemek::util
+{
+
+	struct shared_blob
+	{
+		shared_blob() = default;
+		shared_blob(shared_blob const & other);
+		shared_blob(shared_blob && other);
+		shared_blob(blob && other);
+
+		shared_blob(std::size_t size);
+		shared_blob(std::size_t size, std::shared_ptr<char[]> data);
+		shared_blob(std::size_t size, char * data);
+
+		shared_blob & operator = (shared_blob const & other);
+		shared_blob & operator = (shared_blob && other);
+		shared_blob & operator = (blob && other);
+
+		~ shared_blob() = default;
+
+		char * data() { return data_.get(); }
+		char const * data() const { return data_.get(); }
+
+		std::size_t size() const { return size_; }
+
+		void reset();
+
+		using iterator = char *;
+		using const_iterator = char const *;
+
+		char * begin() { return data(); }
+		char const * begin() const { return data(); }
+
+		char * end() { return data() + size(); }
+		char const * end() const { return data() + size(); }
+
+		char & operator[] (std::size_t i) { return data()[i]; }
+		char const & operator[] (std::size_t i) const { return data()[i]; }
+
+		std::string string() const;
+		std::string_view string_view() const;
+
+	private:
+		std::shared_ptr<char[]> data_;
+		std::size_t size_ = 0;
+	};
+
+	inline shared_blob::shared_blob(shared_blob const & other)
+	{
+		*this = other;
+	}
+
+	inline shared_blob::shared_blob(shared_blob && other)
+	{
+		*this = std::move(other);
+	}
+
+	inline shared_blob::shared_blob(blob && other)
+	{
+		*this = std::move(other);
+	}
+
+	inline shared_blob::shared_blob(std::size_t size)
+		: data_{new char [size]}
+		, size_{size}
+	{}
+
+	inline shared_blob::shared_blob(std::size_t size, std::shared_ptr<char[]> data)
+		: data_{std::move(data)}
+		, size_{size}
+	{}
+
+	inline shared_blob::shared_blob(std::size_t size, char * data)
+		: data_{data}
+		, size_{size}
+	{}
+
+	inline shared_blob & shared_blob::operator=(shared_blob const & other)
+	{
+		if (this != &other)
+		{
+			data_ = other.data_;
+			size_ = other.size();
+		}
+		return *this;
+	}
+
+	inline shared_blob & shared_blob::operator=(shared_blob && other)
+	{
+		if (this != &other)
+		{
+			data_ = std::move(other.data_);
+			size_ = other.size();
+			other.size_ = 0;
+		}
+		return *this;
+	}
+
+	inline shared_blob & shared_blob::operator=(blob && other)
+	{
+		size_ = other.size();
+		data_ = other.release();
+		return *this;
+	}
+
+	inline void shared_blob::reset()
+	{
+		data_.reset();
+		size_ = 0;
+	}
+
+	inline std::string shared_blob::string() const
+	{
+		return std::string(data_.get(), data_.get() + size_);
+	}
+
+	inline std::string_view shared_blob::string_view() const
+	{
+		return std::string_view(data_.get(), size_);
+	}
+}

libs/util/include/psemek/util/synchronyzed_queue.hpp

@@ -0,0 +1,147 @@
+#pragma once
+
+#include <deque>
+#include <mutex>
+#include <condition_variable>
+#include <optional>
+#include <thread>
+#include <chrono>
+
+namespace psemek::util
+{
+
+	template <typename T>
+	struct synchronized_queue
+	{
+		synchronized_queue (std::size_t max_size = std::numeric_limits<std::size_t>::max()) noexcept
+			: max_size_(max_size)
+		{ }
+
+		std::size_t max_size ( ) const noexcept
+		{
+			return max_size_;
+		}
+
+		void push (T const & x);
+		void push (T && x);
+		T pop ( );
+
+		bool try_push (T const & x);
+		template <typename Rep, typename Period>
+		bool try_push (T const & x, std::chrono::duration<Rep, Period> const & timeout);
+
+		std::optional<T> try_pop ( );
+		template <typename Rep, typename Period>
+		std::optional<T> try_pop (std::chrono::duration<Rep, Period> const & timeout);
+
+		void clear ( );
+
+		// Wait for the queue to become empty
+		// e.g. when no new items are going to be pushed
+		void wait ( );
+
+	private:
+		std::mutex mutex;
+		std::condition_variable push_cv, pop_cv;
+		std::deque<T> queue;
+		std::size_t const max_size_;
+	};
+
+	template <typename T>
+	void synchronized_queue<T>::push (T const & x)
+	{
+		std::unique_lock lock { mutex };
+		push_cv.wait(lock, [this]{ return queue.size() < max_size(); });
+		queue.push_back(x);
+		pop_cv.notify_one();
+	}
+
+	template <typename T>
+	void synchronized_queue<T>::push (T && x)
+	{
+		std::unique_lock lock { mutex };
+		push_cv.wait(lock, [this]{ return queue.size() <  max_size(); });
+		queue.push_back(std::move(x));
+		pop_cv.notify_one();
+	}
+
+	template <typename T>
+	T synchronized_queue<T>::pop ( )
+	{
+		std::unique_lock lock { mutex };
+		pop_cv.wait(lock, [this]{ return !queue.empty(); });
+		T x = std::move(queue.front());
+		queue.pop_front();
+		push_cv.notify_one();
+		return x;
+	}
+
+	template <typename T>
+	bool synchronized_queue<T>::try_push (T const & x)
+	{
+		std::lock_guard lock { mutex };
+		if (queue.size() >=  max_size())
+			return false;
+
+		queue.push_back(x);
+		pop_cv.notify_one();
+		return true;
+	}
+
+	template <typename T>
+	template <typename Rep, typename Period>
+	bool synchronized_queue<T>::try_push (T const & x, std::chrono::duration<Rep, Period> const & timeout)
+	{
+		std::unique_lock lock { mutex };
+		if (push_cv.wait_for(lock, timeout, [this]{ return queue.size() <  max_size(); }))
+		{
+			queue.push_back(std::move(x));
+			pop_cv.notify_one();
+			return true;
+		}
+		return false;
+	}
+
+	template <typename T>
+	std::optional<T> synchronized_queue<T>::try_pop ( )
+	{
+		std::lock_guard lock { mutex };
+		if (queue.empty())
+			return std::nullopt;
+		T x = std::move(queue.front());
+		queue.pop_front();
+		push_cv.notify_one();
+		return { std::move(x) };
+	}
+
+	template <typename T>
+	template <typename Rep, typename Period>
+	std::optional<T> synchronized_queue<T>::try_pop (std::chrono::duration<Rep, Period> const & timeout)
+	{
+		std::unique_lock lock { mutex };
+		if (pop_cv.wait_for(lock, timeout, [this]{ return !queue.empty(); }))
+		{
+			T x = std::move(queue.front());
+			queue.pop_front();
+			push_cv.notify_one();
+			return { std::move(x) };
+		}
+		return std::nullopt;
+	}
+
+	template <typename T>
+	void synchronized_queue<T>::clear ( )
+	{
+		std::lock_guard lock { mutex };
+		queue.clear();
+		push_cv.notify_all();
+	}
+
+	template <typename T>
+	void synchronized_queue<T>::wait ( )
+	{
+		std::unique_lock lock { mutex };
+		push_cv.wait(lock, [this]{ return queue.empty(); });
+	}
+
+}

libs/util/include/psemek/util/thread.hpp

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <thread>
+
+namespace psemek::util
+{
+
+	struct thread
+		: std::thread
+	{
+		template <typename ... Args>
+		thread (Args && ... args)
+			: std::thread(std::forward<Args>(args)...)
+		{ }
+
+		thread (thread &&) = default;
+
+		~ thread ( )
+		{
+			if (joinable())
+				join();
+		}
+	};
+
+}

libs/util/include/psemek/util/threadpool.hpp

@@ -0,0 +1,57 @@
+#pragma once
+
+#include <psemek/util/thread.hpp>
+#include <psemek/util/synchronyzed_queue.hpp>
+#include <psemek/util/movable_function.hpp>
+
+#include <future>
+#include <vector>
+
+namespace psemek::util
+{
+
+	struct threadpool
+	{
+		threadpool ( )
+			: threadpool(std::max(1u, std::thread::hardware_concurrency()))
+		{ }
+
+		threadpool (std::size_t thread_count)
+		{
+			start(thread_count);
+		}
+
+		~ threadpool ( )
+		{
+			stop();
+		}
+
+		template <typename F>
+		auto dispatch (F && f)
+		{
+			using R = decltype(f());
+
+			std::packaged_task<R()> task { std::forward<F>(f) };
+
+			auto result = task.get_future();
+
+			tasks_queue.push(std::move(task));
+
+			return result;
+		}
+
+		void start (std::size_t thread_count);
+
+		void stop ( );
+
+		void wait ( )
+		{
+			tasks_queue.wait();
+		}
+
+	private:
+		std::vector<util::thread> threads;
+		util::synchronized_queue<movable_function<void()>> tasks_queue;
+	};
+
+}

libs/util/include/psemek/util/timer.hpp

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <psemek/util/clock.hpp>
+
+namespace psemek::util
+{
+
+	template <typename Duration = std::chrono::duration<double>, typename Clock = std::chrono::system_clock>
+	struct timer
+		: clock<Duration, Clock>
+	{
+		timer (Duration duration)
+			: duration_(duration)
+		{ }
+
+		explicit operator bool ( )
+		{
+			if (this->duration() >= duration_)
+			{
+				this->restart();
+				return true;
+			}
+
+			return false;
+		}
+
+	private:
+		Duration duration_;
+	};
+
+}

libs/util/include/psemek/util/to_string.hpp

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <sstream>
+#include <stdexcept>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		template <typename Char, typename Traits = std::char_traits<Char>, typename ... Args>
+		std::basic_string<Char, Traits> to_string (Args const & ... args)
+		{
+			std::basic_ostringstream<Char, Traits> oss;
+
+			((oss << args), ...);
+
+			return oss.str();
+		}
+
+	}
+
+	template <typename ... Args>
+	std::string to_string (Args const & ... args)
+	{
+		return detail::to_string<char>(args...);
+	}
+
+	template <typename ... Args>
+	std::wstring to_wstring (Args const & ... args)
+	{
+		return detail::to_string<wchar_t>(args...);
+	}
+
+	template <typename ... Args>
+	std::u32string to_u32string (Args const & ... args)
+	{
+		return detail::to_string<char32_t>(args...);
+	}
+
+	template <typename T, typename Char, typename Traits>
+	T from_string (std::basic_string<Char, Traits> const & s)
+	{
+		std::basic_istringstream<Char, Traits> iss(s);
+		T x;
+		iss >> x;
+		if (!iss)
+			throw std::invalid_argument("failed to parse from string");
+		return x;
+	}
+
+	template <typename T, typename Char>
+	T from_string (Char const * s)
+	{
+		return from_string<T, Char, std::char_traits<Char>>(s);
+	}
+
+}

libs/util/include/psemek/util/unicode.hpp

@@ -0,0 +1,11 @@
+#pragma once
+
+#include <string>
+
+namespace psemek::util
+{
+
+	std::string to_utf8 (std::u32string const & str);
+	std::u32string from_utf8 (std::string const & str);
+
+}

libs/util/source/assert.cpp

@@ -0,0 +1,14 @@
+#include <psemek/util/assert.hpp>
+#include <psemek/util/to_string.hpp>
+
+#include <stdexcept>
+
+namespace psemek::util
+{
+
+	[[noreturn]] bool assert_handler (char const * expression, char const * file, int line)
+	{
+		throw std::runtime_error(to_string(file, ":", line, " Assertion failed: ", expression));
+	}
+
+}

libs/util/source/movable_function.cpp

@@ -0,0 +1,18 @@
+#include <psemek/util/movable_function.hpp>
+
+#include <functional>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		void bad_function_call ( )
+		{
+			throw std::bad_function_call();
+		}
+
+	}
+
+}

libs/util/source/not_implemented.cpp

@@ -0,0 +1,13 @@
+#include <psemek/util/not_implemented.hpp>
+
+#include <stdexcept>
+
+namespace psemek::util
+{
+
+	void not_implemented ( )
+	{
+		throw std::runtime_error("Not implemented");
+	}
+
+}

libs/util/source/pixmap.cpp

@@ -0,0 +1,110 @@
+#include <psemek/util/pixmap.hpp>
+
+#include <iostream>
+
+namespace psemek::util
+{
+
+	pixmap_monochrome read_pgm(std::istream & is)
+	{
+		auto fail = [](std::string str)
+		{
+			throw std::runtime_error("Error loading PGM image: " + str);
+		};
+
+		std::string line;
+		std::getline(is, line);
+
+		bool binary = true;
+
+		if (line == "P2")
+		{
+			binary = false;
+		}
+		else if (line == "P5")
+		{
+			binary = true;
+		}
+		else
+			fail("unknown format " + line);
+
+		std::size_t width, height;
+		std::size_t max;
+		pixmap_monochrome pixmap;
+
+		is >> width >> height >> max;
+
+		if (max != 255)
+			fail("max value " + std::to_string(max) + " is not supported");
+
+		pixmap.resize(width, height);
+
+		if (binary)
+			is.read(reinterpret_cast<char *>(pixmap.data()), width * height * sizeof(pixmap.data()[0]));
+		else
+			fail("P2 format is not supported");
+
+		if (!is)
+			fail("stream error");
+
+		return pixmap;
+	}
+
+	pixmap_rgb load_ppm(std::istream & is)
+	{
+		auto fail = [](std::string str)
+		{
+			throw std::runtime_error("Error loading PPM image: " + str);
+		};
+
+		std::string line;
+		std::getline(is, line);
+
+		bool binary = true;
+
+		if (line == "P3")
+		{
+			binary = false;
+		}
+		else if (line == "P6")
+		{
+			binary = true;
+		}
+		else
+			fail("unknown format " + line);
+
+		std::size_t width, height;
+		std::size_t max;
+		pixmap_rgb pixmap;
+
+		is >> width >> height >> max;
+
+		if (max != 255)
+			fail("max value " + std::to_string(max) + " is not supported");
+
+		pixmap.resize(width, height);
+
+		if (binary)
+			is.read(reinterpret_cast<char *>(pixmap.data()), width * height * sizeof(pixmap.data()[0]));
+		else
+			fail("P3 format is not supported");
+
+		if (!is)
+			fail("stream error");
+
+		return pixmap;
+	}
+
+	void write_pgm(pixmap_monochrome const & p, std::ostream & os)
+	{
+		os << "P5\n" << p.width() << " " << p.height() << "\n255\n";
+		os.write(reinterpret_cast<char const *>(p.data()), p.width() * p.height() * sizeof(p.data()[0]));
+	}
+
+	void write_ppm(pixmap_rgb const & p, std::ostream & os)
+	{
+		os << "P6\n" << p.width() << " " << p.height() << "\n255\n";
+		os.write(reinterpret_cast<char const *>(p.data()), p.width() * p.height() * sizeof(p.data()[0]));
+	}
+
+}

libs/util/source/pretty_print_time.cpp

@@ -0,0 +1,66 @@
+#include <psemek/util/pretty_print.hpp>
+
+namespace psemek::util
+{
+
+	namespace detail
+	{
+
+		void pretty_print_time (std::ostream & o, std::int64_t d, std::int64_t up_to)
+		{
+			static constexpr const std::int64_t durations[8] = {
+				604800 * 1000000000ull,
+				86400 * 1000000000ull,
+				3600 * 1000000000ull,
+				60 * 1000000000ull,
+				1000000000ull,
+				1000000ull,
+				1000ull,
+				1ull
+			};
+
+			static const std::string_view suffixes[8] = {
+				"w",
+				"d",
+				"h",
+				"m",
+				"s",
+				"ms",
+				"us",
+				"ns"
+			};
+
+			bool first = true;
+
+			for (std::size_t i = 0; i < 8; ++i)
+			{
+				if (up_to > durations[i]) break;
+
+				std::int64_t value = d / durations[i];
+				d %= durations[i];
+
+				if (value == 0) continue;
+
+				if (first)
+					first = false;
+				else
+					o << ' ';
+
+				o << value << suffixes[i];
+			}
+
+			// [ first == true ] means no output was generated
+			// Find first matching suffix and output zero
+			if (first) for (std::size_t i = 0; i < 8; ++i)
+			{
+				if (up_to >= durations[i])
+				{
+					o << 0 << suffixes[i];
+					break;
+				}
+			}
+		}
+
+	}
+
+}

libs/util/source/threadpool.cpp

@@ -0,0 +1,35 @@
+#include <psemek/util/threadpool.hpp>
+
+namespace psemek::util
+{
+
+	void threadpool::start (std::size_t thread_count)
+	{
+		for (std::size_t th = 0; th < thread_count; ++th)
+		{
+			threads.emplace_back([this]
+			{
+				while (true)
+				{
+					auto task = tasks_queue.pop();
+
+					if (!task)
+						break;
+
+					task();
+				}
+			});
+		}
+	}
+
+	void threadpool::stop ( )
+	{
+		tasks_queue.clear();
+		for (auto const & thread: threads)
+		{
+			tasks_queue.push({});
+		}
+		threads.clear();
+	}
+
+}

libs/util/source/unicode.cpp

@@ -0,0 +1,21 @@
+#include <psemek/util/unicode.hpp>
+
+#include <codecvt>
+#include <locale>
+
+namespace psemek::util
+{
+
+	using converter = std::wstring_convert<std::codecvt_utf8<char32_t>, char32_t>;
+
+	std::string to_utf8 (std::u32string const & str)
+	{
+		return converter().to_bytes(str);
+	}
+
+	std::u32string from_utf8 (std::string const & str)
+	{
+		return converter().from_bytes(str);
+	}
+
+}

todo.md

@@ -0,0 +1,10 @@
+* Use the same code style everywhere (mostly about spaces after function names)
+* Make sure program & shaders are deleted properly if program creation fails
+* Remove gfx::vertex, setup mesh using attributes directly
+* Design affine transforms in geom & use them instead of matrices when appropriate
+* Create an 'app' module that simplifies application creation
+* Implement pixmap font rendering
+* Create a simple generic primive painter
+* Design resources system
+* Design logging system
+* Design ui system