diff --git a/libs/cg/include/psemek/cg/triangulation/monotone.hpp b/libs/cg/include/psemek/cg/triangulation/monotone.hpp
new file mode 100644
index 00000000..10d70c01
--- /dev/null
+++ b/libs/cg/include/psemek/cg/triangulation/monotone.hpp
@@ -0,0 +1,360 @@
+#pragma once
+
+#include <psemek/geom/point.hpp>
+#include <psemek/geom/orientation.hpp>
+#include <psemek/geom/contains.hpp>
+
+#include <vector>
+#include <set>
+
+namespace psemek::cg
+{
+
+	template <typename Iterator, typename OutputIterator, typename IndexType = std::size_t>
+	void monotone_triangulation(Iterator begin, Iterator end, OutputIterator out)
+	{
+		IndexType const count = std::distance(begin, end);
+		IndexType const null = IndexType(-1);
+
+		auto at = [begin](auto i){ return *(begin + i); };
+
+		// Initialize dcel
+
+		std::vector<IndexType> start(count, 0);
+		std::vector<IndexType> next(count, 0);
+		std::vector<IndexType> prev(count, 0);
+		std::vector<IndexType> twin(count, 0);
+
+		for (IndexType i = 0; i < count; ++i)
+		{
+			start[i] = i;
+			next[i] = (i + 1) % count;
+			prev[i] = (i + count - 1) % count;
+			twin[i] = null;
+		}
+
+		auto finish = [&](auto i)
+		{
+			return start[next[i]];
+		};
+
+		struct vertex_id
+		{
+			IndexType index;
+		};
+
+		struct sweepline_compare
+		{
+			struct is_transparent{};
+
+			Iterator begin;
+			IndexType count;
+
+			auto const & at(IndexType i) const
+			{
+				return *(begin + i);
+			}
+
+			bool operator()(IndexType i, IndexType j) const
+			{
+				auto const i0 = at(i);
+				auto const i1 = at((i + 1) % count);
+				auto const j0 = at(j);
+				auto const j1 = at((j + 1) % count);
+
+				auto const oi0 = geom::orientation(i0, i1, j0);
+				auto const oi1 = geom::orientation(i0, i1, j1);
+
+				if (oi0 == oi1)
+					return oi0 == geom::sign_t::positive;
+
+				return geom::orientation(j0, j1, i0) == geom::sign_t::negative;
+			}
+
+			bool operator()(IndexType e, vertex_id const & v) const
+			{
+				return geom::orientation(at(e), at((e + 1) % count), at(v.index)) == geom::sign_t::positive;
+			}
+
+			bool operator()(vertex_id const & v, IndexType e) const
+			{
+				return geom::orientation(at(e), at((e + 1) % count), at(v.index)) == geom::sign_t::negative;
+			}
+		};
+
+		// Edges intersecting sweepline
+		std::set<IndexType, sweepline_compare> sweepline(sweepline_compare{begin, count});
+
+		std::vector<IndexType> events(count, 0);
+		std::vector<IndexType> helper(count, 0);
+
+		for (IndexType i = 0; i < count; ++i)
+			events[i] = i;
+
+		std::sort(events.begin(), events.end(), [&](auto i, auto j){ return at(i) < at(j); });
+
+		enum class event_type
+		{
+			start,
+			end,
+			split,
+			merge,
+			regular,
+		};
+
+		auto classify = [&](auto e) -> event_type
+		{
+			auto p = (e + count - 1) % count;
+			auto n = (e + 1) % count;
+
+			bool lp = at(p) < at(e);
+			bool ln = at(n) < at(e);
+
+			auto s = geom::orientation(at(p), at(e), at(n));
+
+			if (!lp && !ln && s == geom::sign_t::positive)
+				return event_type::start;
+			else if (lp && ln && s == geom::sign_t::positive)
+				return event_type::end;
+			else if (!lp && !ln && s == geom::sign_t::negative)
+				return event_type::split;
+			else if (lp && ln && s == geom::sign_t::negative)
+				return event_type::merge;
+			else
+				return event_type::regular;
+		};
+
+		auto connect = [&](auto i, auto j)
+		{
+			auto const vi = at(i);
+			auto const vj = at(j);
+
+			auto ni = i;
+			while (true)
+			{
+				auto f = at(finish(ni));
+				auto p = at(start[prev[ni]]);
+				bool contains;
+				if (geom::orientation(f, vi, p) == geom::sign_t::negative)
+					contains = geom::orientation(f, vi, vj) == geom::sign_t::negative && geom::orientation(vj, vi, p) == geom::sign_t::negative;
+				else
+					contains = geom::orientation(f, vi, vj) == geom::sign_t::negative || geom::orientation(vj, vi, p) == geom::sign_t::negative;
+				if (contains)
+					break;
+				ni = twin[prev[ni]];
+			}
+			auto pi = prev[ni];
+
+			auto nj = j;
+			while (true)
+			{
+				auto f = at(finish(nj));
+				auto p = at(start[prev[nj]]);
+				bool contains;
+				if (geom::orientation(f, vj, p) == geom::sign_t::negative)
+					contains = geom::orientation(f, vj, vi) == geom::sign_t::negative && geom::orientation(vi, vj, p) == geom::sign_t::negative;
+				else
+					contains = geom::orientation(f, vj, vi) == geom::sign_t::negative || geom::orientation(vi, vj, p) == geom::sign_t::negative;
+				if (contains)
+					break;
+				nj = twin[prev[nj]];
+			}
+			auto pj = prev[nj];
+
+			IndexType eij = start.size();
+			IndexType eji = eij + 1;
+
+			start.emplace_back();
+			start.emplace_back();
+			next.emplace_back();
+			next.emplace_back();
+			prev.emplace_back();
+			prev.emplace_back();
+			twin.emplace_back();
+			twin.emplace_back();
+
+			start[eij] = i;
+			next[eij] = nj;
+			prev[eij] = pi;
+			twin[eij] = eji;
+
+			start[eji] = j;
+			next[eji] = ni;
+			prev[eji] = pj;
+			twin[eji] = eij;
+
+			next[pi] = eij;
+			prev[nj] = eij;
+
+			next[pj] = eji;
+			prev[ni] = eji;
+		};
+
+		// Partition into monotone parts
+		for (auto e : events)
+		{
+			auto p = (e + count - 1) % count;
+
+			switch (classify(e))
+			{
+			case event_type::start:
+				sweepline.insert(e);
+				helper[e] = e;
+				break;
+			case event_type::end:
+				if (auto h = helper[p]; classify(h) == event_type::merge)
+					connect(e, h);
+				sweepline.erase(p);
+				break;
+			case event_type::split:
+				{
+					auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
+					connect(e, helper[l]);
+					helper[l] = e;
+				}
+				sweepline.insert(e);
+				helper[e] = e;
+				break;
+			case event_type::merge:
+				if (auto h = helper[p]; classify(h) == event_type::merge)
+					connect(e, h);
+				sweepline.erase(p);
+				{
+					auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
+					if (auto h = helper[l]; classify(h) == event_type::merge)
+						connect(e, h);
+					helper[l] = e;
+				}
+				break;
+			case event_type::regular:
+				if (at(e) < at((e + 1) % count))
+				{
+					if (auto h = helper[p]; classify(h) == event_type::merge)
+						connect(e, h);
+					sweepline.erase(p);
+					sweepline.insert(e);
+					helper[e] = e;
+				}
+				else
+				{
+					auto l = *std::prev(sweepline.upper_bound(vertex_id{e}));
+					if (auto h = helper[l]; classify(h) == event_type::merge)
+						connect(e, h);
+					helper[l] = e;
+				}
+				break;
+			}
+		}
+
+		// Triangulate monotone parts
+		std::vector<bool> used(start.size(), false);
+		std::vector<IndexType> stack;
+		for (IndexType e = 0; e < start.size(); ++e)
+		{
+			if (used[e]) continue;
+
+			auto le = e;
+			used[e] = true;
+			// Mark the whole part as used and find leftmost point
+			for (auto ee = next[e]; ee != e; ee = next[ee])
+			{
+				used[ee] = true;
+				if (at(start[ee]) < at(start[le]))
+					le = ee;
+			}
+
+			bool stack_left = true;
+
+			auto re = prev[le];
+			stack.push_back(start[le]);
+			le = next[le];
+
+			while (true)
+			{
+				bool last = (le == re);
+
+				IndexType v;
+				bool left;
+				if (at(start[le]) < at(start[re]))
+				{
+					v = start[le];
+					le = next[le];
+					left = true;
+				}
+				else
+				{
+					v = start[re];
+					re = prev[re];
+					left = false;
+				}
+
+				if (left && stack_left)
+				{
+					while (stack.size() >= 2 && geom::orientation(at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::negative)
+					{
+						*out++ = v;
+						*out++ = stack[stack.size() - 2];
+						*out++ = stack.back();
+						stack.pop_back();
+					}
+					stack.push_back(v);
+				}
+				else if (!left && !stack_left)
+				{
+					while (stack.size() >= 2 && geom::orientation(at(v), at(stack.back()), at(stack[stack.size() - 2])) == geom::sign_t::positive)
+					{
+						*out++ = v;
+						*out++ = stack.back();
+						*out++ = stack[stack.size() - 2];
+						stack.pop_back();
+					}
+					stack.push_back(v);
+				}
+				else if (left && !stack_left)
+				{
+					auto last = stack.back();
+					while (stack.size() >= 2)
+					{
+						*out++ = v;
+						*out++ = stack.back();
+						*out++ = stack[stack.size() - 2];
+						stack.pop_back();
+					}
+					stack.pop_back();
+					stack.push_back(last);
+					stack.push_back(v);
+				}
+				else if (!left && stack_left)
+				{
+					auto last = stack.back();
+					while (stack.size() >= 2)
+					{
+						*out++ = v;
+						*out++ = stack[stack.size() - 2];
+						*out++ = stack.back();
+						stack.pop_back();
+					}
+					stack.pop_back();
+					stack.push_back(last);
+					stack.push_back(v);
+				}
+
+				stack_left = left;
+
+				if (last)
+					break;
+			}
+
+			stack.clear();
+		}
+	}
+
+	template <typename IndexType, typename Iterator>
+	std::vector<IndexType> monotone_triangulation(Iterator begin, Iterator end)
+	{
+		std::vector<IndexType> result;
+		monotone_triangulation<Iterator, std::back_insert_iterator<std::vector<IndexType>>, IndexType>(begin, end, std::back_inserter(result));
+		return result;
+	}
+
+}