#pragma once

#include <psemek/math/box.hpp>
#include <psemek/math/simplex.hpp>
#include <psemek/util/hash_table.hpp>

#include <vector>
#include <cstdint>

namespace psemek::math
{

	template <typename T, std::size_t N>
	std::vector<point<T, N>> vertices(box<T, N> const & b)
	{
		std::vector<point<T, N>> result;

		for (std::size_t mask = 0; mask < (1 << N); ++mask)
		{
			point<T, N> p;
			for (std::size_t i = 0; i < N; ++i)
			{
				p[i] = (mask & (1 << i)) ? b[i].max : b[i].min;
			}
			result.push_back(p);
		}

		return result;
	}

	template <typename Index = std::uint32_t, typename T, std::size_t N>
	std::vector<segment<Index>> edges(box<T, N> const &)
	{
		static_assert(N > 0); // who knows?..

		std::vector<segment<Index>> result;

		for (std::size_t i = 0; i < N; ++i)
		{
			for (std::size_t mask = 0; mask < (1 << (N - 1)); ++mask)
			{
				auto const lo = mask & ((1 << i) - 1);
				auto const hi = mask ^ lo;

				auto const idx0 = lo | (hi << 1);
				auto const idx1 = idx0 | (1 << i);

				result.push_back({static_cast<Index>(idx0), static_cast<Index>(idx1)});
			}
		}

		return result;
	}

	template <typename Index = std::uint32_t, typename T>
	std::vector<triangle<Index>> faces(box<T, 2> const &)
	{
		std::vector<triangle<Index>> result;
		result.push_back({0, 1, 2});
		result.push_back({2, 1, 3});
		return result;
	}

	template <typename Index = std::uint32_t, typename T>
	std::vector<triangle<Index>> faces(box<T, 3> const &)
	{
		std::vector<triangle<Index>> result;
		// -Z
		result.push_back({0, 2, 1});
		result.push_back({1, 2, 3});
		// +Z
		result.push_back({4, 5, 6});
		result.push_back({6, 5, 7});
		// -Y
		result.push_back({0, 1, 4});
		result.push_back({4, 1, 5});
		// +Y
		result.push_back({2, 6, 3});
		result.push_back({6, 7, 3});
		// -X
		result.push_back({0, 4, 2});
		result.push_back({2, 4, 6});
		// +X
		result.push_back({1, 3, 5});
		result.push_back({3, 7, 5});
		return result;
	}

	template <typename Vertex, typename Index, std::size_t N>
	std::vector<simplex<Vertex, N>> deindex(std::vector<Vertex> const & vertices, std::vector<simplex<Index, N>> const & simplices)
	{
		std::vector<simplex<Vertex, N>> result;
		result.reserve(simplices.size());
		for (auto const & s : simplices)
		{
			result.emplace_back();
			for (std::size_t i = 0; i <= N; ++i)
				result.back()[i] = vertices[s[i]];
		}
		return result;
	}

	namespace detail
	{

		template <typename F, typename G, std::size_t ... I>
		auto map_index_sequence(F && f, G && g, std::index_sequence<I...>)
		{
			return f(g(I)...);
		}

	}

	template <typename T, std::size_t N, typename Index>
	std::vector<vector<T, N + 1>> flat_normals(std::vector<point<T, N + 1>> const & vertices, std::vector<simplex<Index, N>> const & simplices)
	{
		std::vector<vector<T, N + 1>> result;
		result.reserve(simplices.size());
		for (auto const & s : simplices)
		{
			auto const n = detail::map_index_sequence([&](auto const & ... vs){ return normal(vs...); }, [&](std::size_t i){ return vertices[s[i]]; }, std::make_index_sequence<N + 1>{});
			result.push_back(n);
		}
		return result;
	}

	template <typename Vertex, typename Index>
	auto smooth_normals(std::vector<Vertex> const & vertices, std::vector<triangle<Index>> const & triangles)
	{
		using vector_type = std::decay_t<decltype(vertices[0] - vertices[0])>;
		std::vector<vector_type> normals(vertices.size(), vector_type::zero());

		for (auto const & t : triangles)
		{
			auto const n = cross(vertices[t[1]] - vertices[t[0]], vertices[t[2]] - vertices[t[0]]);
			normals[t[0]] += n;
			normals[t[1]] += n;
			normals[t[2]] += n;
		}

		for (auto & n : normals)
			n = normalized(n);

		return normals;
	}

	template <typename Index = std::uint32_t, typename Vertex>
	std::vector<triangle<Index>> triangulate_convex(std::vector<Vertex> const & v)
	{
		std::vector<triangle<Index>> result;

		for (std::size_t i = 1; i + 1 < v.size(); ++i)
		{
			result.push_back({static_cast<Index>(0), static_cast<Index>(i), static_cast<Index>(i + 1)});
		}

		return result;
	}

	template <typename Vertex, typename Index>
	void subdivide(std::vector<Vertex> & vertices, std::vector<triangle<Index>> & triangles)
	{
		using edge = std::pair<Index, Index>;

		using scalar = std::decay_t<decltype(vertices[0][0])>;

		auto normalise = [](edge e) -> edge
		{
			if (e.first > e.second)
				std::swap(e.first, e.second);
			return e;
		};

		util::hash_map<edge, Index> midpoints;

		auto add_midpoint = [&](Index i0, Index i1)
		{
			auto ne = normalise({i0, i1});

			if (!midpoints.contains(ne))
			{
				midpoints[ne] = vertices.size();
				vertices.push_back(lerp(vertices[i0], vertices[i1], scalar{0.5}));
			}
		};

		for (auto const & t : triangles)
		{
			add_midpoint(t[0], t[1]);
			add_midpoint(t[1], t[2]);
			add_midpoint(t[2], t[0]);
		}

		auto get_midpoint = [&](Index i0, Index i1)
		{
			return midpoints.at(normalise({i0, i1}));
		};

		std::vector<triangle<Index>> new_triangles;

		for (auto const & t : triangles)
		{
			new_triangles.push_back({t[0], get_midpoint(t[0], t[1]), get_midpoint(t[2], t[0])});
			new_triangles.push_back({t[1], get_midpoint(t[1], t[2]), get_midpoint(t[0], t[1])});
			new_triangles.push_back({t[2], get_midpoint(t[2], t[0]), get_midpoint(t[1], t[2])});
			new_triangles.push_back({get_midpoint(t[0], t[1]), get_midpoint(t[1], t[2]), get_midpoint(t[2], t[0])});
		}

		triangles.swap(new_triangles);
	}

}
