#pragma once

#include <psemek/math/point.hpp>
#include <psemek/math/orientation.hpp>
#include <psemek/math/contains.hpp>
#include <psemek/cg/dcel.hpp>

#include <vector>
#include <set>

namespace psemek::cg
{

	template <typename IndexType = std::size_t, typename RobustTag, typename Iterator>
	auto monotone_triangulation(RobustTag robust_tag, Iterator begin, Iterator end)
	{
		IndexType const count = std::distance(begin, end);

		auto at = [begin](IndexType i){ return *(begin + i); };

		// Initialize dcel

		using dcel_type = dcel<util::empty, util::empty, util::empty, IndexType>;
		using point_handle = typename dcel_type::point_handle;

		dcel_type result;

		result.push_face();
		result.faces[0].edge = count;

		result.points.resize(count);
		result.edges.resize(2 * count);

		for (IndexType i = 0; i < count; ++i)
		{
			result.points[i].edge = i;

			result.edges[i].face = result.null;
			result.edges[i].origin = i;
			result.edges[i].prev = (i + count - 1) % count;
			result.edges[i].next = (i + 1) % count;
			result.edges[i].twin = i + count;

			result.edges[i + count].face = 0;
			result.edges[i + count].origin = (i + 1) % count;
			result.edges[i + count].prev = count + (i + 1) % count;
			result.edges[i + count].next = count + (i + count - 1) % count;
			result.edges[i + count].twin = i;
		}

		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 = math::orientation(RobustTag{}, i0, i1, j0);
				auto const oi1 = math::orientation(RobustTag{}, i0, i1, j1);

				if (oi0 == oi1)
					return oi0 == math::sign_t::positive;

				return math::orientation(RobustTag{}, j0, j1, i0) == math::sign_t::negative;
			}

			bool operator()(IndexType e, point_handle v) const
			{
				return math::orientation(RobustTag{}, at(e), at((e + 1) % count), at(v.index())) == math::sign_t::positive;
			}

			bool operator()(point_handle v, IndexType e) const
			{
				return math::orientation(RobustTag{}, at(e), at((e + 1) % count), at(v.index())) == math::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 = math::orientation(robust_tag, at(p), at(e), at(n));

			if (!lp && !ln && s == math::sign_t::positive)
				return event_type::start;
			else if (lp && ln && s == math::sign_t::positive)
				return event_type::end;
			else if (!lp && !ln && s == math::sign_t::negative)
				return event_type::split;
			else if (lp && ln && s == math::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 = result.edge(i);
			while (true)
			{
				auto f = at(ni.next().origin().index());
				auto p = at(ni.prev().origin().index());
				bool contains;
				if (math::orientation(robust_tag, f, vi, p) == math::sign_t::negative)
					contains = math::orientation(robust_tag, f, vi, vj) == math::sign_t::negative && math::orientation(robust_tag, vj, vi, p) == math::sign_t::negative;
				else
					contains = math::orientation(robust_tag, f, vi, vj) == math::sign_t::negative || math::orientation(robust_tag, vj, vi, p) == math::sign_t::negative;
				if (contains)
					break;
				ni = ni.prev().twin();
			}
			auto pi = ni.prev();

			auto nj = result.edge(j);
			while (true)
			{
				auto f = at(nj.next().origin().index());
				auto p = at(nj.prev().origin().index());
				bool contains;
				if (math::orientation(robust_tag, f, vj, p) == math::sign_t::negative)
					contains = math::orientation(robust_tag, f, vj, vi) == math::sign_t::negative && math::orientation(robust_tag, vi, vj, p) == math::sign_t::negative;
				else
					contains = math::orientation(robust_tag, f, vj, vi) == math::sign_t::negative || math::orientation(robust_tag, vi, vj, p) == math::sign_t::negative;
				if (contains)
					break;
				nj = nj.prev().twin();
			}
			auto pj = nj.prev();

			auto eij = result.push_edge();
			auto eji = result.push_edge();

			eij.origin(result.point(i));
			eij.next(nj); nj.prev(eij);
			eij.prev(pi); pi.next(eij);
			eij.twin(eji);

			eji.origin(result.point(j));
			eji.next(ni); ni.prev(eji);
			eji.prev(pj); pj.next(eji);
			eji.twin(eij);
		};

		// 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(result.point(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(result.point(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(result.point(e)));
					if (auto h = helper[l]; classify(h) == event_type::merge)
						connect(e, h);
					helper[l] = e;
				}
				break;
			}
		}

		fill_missing_faces(result);

		IndexType face_count = result.faces.size();

		// Triangulate monotone parts
		std::vector<IndexType> stack;
		for (IndexType f = 1; f < face_count; ++f)
		{
			auto e = result.face(f).edge();

			// Find leftmost point
			auto le = e;
			for (auto ee = e.next(); ee != e; ee = ee.next())
			{
				if (at(ee.origin().index()) < at(le.origin().index()))
					le = ee;
			}

			bool stack_left = true;

			auto re = le.prev();
			stack.push_back(le.index());
			le = le.next();

			auto emit_triangle = [&](auto e0_id, auto e1_id)
			{
				auto e0 = result.edge(e0_id);
				auto e1 = result.edge(e1_id);

				auto p = e0.prev();
				auto n = e1.next();

				auto face_out = e0.face();
				auto face_in = result.push_face();

				auto e_in = result.push_edge();
				auto e_out = result.push_edge();

				face_in.edge(e_in);
				face_out.edge(e_out);

				e_in.origin(n.origin());
				e_out.origin(e0.origin());

				e_in.twin(e_out);
				e_out.twin(e_in);

				e0.face(face_in);
				e1.face(face_in);
				e_in.face(face_in);
				e_out.face(face_out);

				e0.prev(e_in); e_in.next(e0);
				e1.next(e_in); e_in.prev(e1);

				p.next(e_out); e_out.prev(p);
				n.prev(e_out); e_out.next(n);

				return e_out.index();
			};

			while (true)
			{
				auto ve = le;
				bool left;
				if (at(le.origin().index()) < at(re.origin().index()))
				{
					ve = le;
					le = le.next();
					left = true;
				}
				else
				{
					ve = re;
					re = re.prev();
					left = false;
				}

				auto orientation = [&](auto e0, auto e1, auto e2)
				{
					auto point = [&](auto e){ return at(result.edge(e).origin().index()); };
					return math::orientation(robust_tag, point(e0), point(e1), point(e2));
				};

				if (left && stack_left)
				{
					while (stack.size() >= 2 && orientation(ve.index(), stack.back(), stack[stack.size() - 2]) == math::sign_t::negative)
					{
						auto new_edge = emit_triangle(stack[stack.size() - 2], stack.back());
						stack.pop_back();
						stack.pop_back();
						stack.push_back(new_edge);
					}
					stack.push_back(ve.index());
				}
				else if (!left && !stack_left)
				{
					auto e = ve.index();
					while (stack.size() >= 2 && orientation(ve.index(), stack.back(), stack[stack.size() - 2]) == math::sign_t::positive)
					{
						e = emit_triangle(e, stack.back());
						stack.pop_back();
					}
					stack.push_back(e);
				}
				else if (left && !stack_left)
				{
					auto e = stack.front();
					for (auto it = std::next(stack.begin()); it != stack.end(); ++it)
						e = emit_triangle(*it, e);
					stack.clear();
					stack.push_back(e);
					stack.push_back(ve.index());
				}
				else if (!left && stack_left)
				{
					auto e = ve.index();
					auto last = stack.back();
					for (auto it = stack.begin(); it != std::prev(stack.end()); ++it)
						e = emit_triangle(e, *it);
					stack.clear();
					stack.push_back(last);
					stack.push_back(e);
				}

				stack_left = left;

				if (le == re)
					break;
			}

			auto de = le.index();
			if (stack_left)
			{
				for (IndexType i = 0; i + 2 < stack.size(); ++i)
					de = emit_triangle(de, stack[i]);
			}
			else
			{
				for (IndexType i = stack.size() - 1; i > 1; --i)
					de = emit_triangle(de, stack[i]);
			}

			stack.clear();
		}

		return result;
	}


	template <typename IndexType, typename Iterator>
	auto monotone_triangulation(Iterator begin, Iterator end)
	{
		return monotone_triangulation<IndexType>(math::default_robust_tag, begin, end);
	}

}
