#pragma once

#include <psemek/gfx/pixmap.hpp>
#include <psemek/math/vector.hpp>
#include <psemek/pcg/seamless.hpp>

#include <vector>

namespace psemek::pcg
{

	template <typename T, std::size_t N>
	struct perlin
	{
		using value_type = T;
		static constexpr std::size_t dimension = N;

		perlin() = default;
		perlin(util::ndarray<math::vector<T, N>, N> grad_map);
		perlin(util::ndarray<math::vector<T, N>, N> const & grad_map, seamless_tag);
		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;
		}

		std::size_t depth() const
		{
			return grad_map_.depth() - 1;
		}

		// Coords \in [0.0 .. 1.0]
		template <typename ... Args>
		T operator()(Args const & ... args) const
		{
			return (*this)(math::vector{args...});
		}

		// Coords \in [0.0 .. 1.0]
		T operator()(math::vector<T, N> p) const;

		auto & grad() { return grad_map_; }
		auto const & grad() const { return grad_map_; }

	private:
		util::ndarray<math::vector<T, N>, N> grad_map_;
	};

	template <typename T, std::size_t N>
	perlin<T, N>::perlin(util::ndarray<math::vector<T, N>, N> grad_map)
		: grad_map_{std::move(grad_map)}
	{}

	template <typename T, std::size_t N>
	perlin<T, N>::perlin(util::ndarray<math::vector<T, N>, N> const & grad_map, seamless_tag)
	{
		auto dims = grad_map.dims();
		for (std::size_t i = 0; i < N; ++i)
			dims[i] += 1;

		grad_map_.resize(dims);

		for (auto const & idx : grad_map_.indices())
		{
			auto idx_orig = idx;
			for (std::size_t i = 0; i < N; ++i)
				idx_orig[i] %= grad_map.dim(i);
			grad_map_(idx) = grad_map(idx_orig);
		}
	}

	template <typename T, std::size_t N>
	T perlin<T, N>::operator()(math::vector<T, N> p) const
	{
		for (std::size_t i = 0; i < N; ++i)
		{
			assert(p[i] >= 0);
			assert(p[i] <= 1);
			p[i] *= grad_map_.dim(i) - 1;
		}

		math::vector<std::size_t, N> ip;
		for (std::size_t i = 0; i < N; ++i)
		{
			ip[i] = math::clamp<std::size_t>(std::floor(p[i]), {0, static_cast<std::size_t>(grad_map_.dim(i)) - 2});
		}

		math::vector<T, N> t = p - math::cast<T>(ip);

		T values[1 << N];

		for (std::size_t mask = 0; mask < (1 << N); ++mask)
		{
			math::vector<T, N> tt;
			std::array<std::size_t, N> ii;
			for (std::size_t i = 0; i < N; ++i)
			{
				bool m = (mask & (1 << i));
				tt[i] = m ? t[i] - 1 : t[i];
				ii[i] = m ? ip[i] + 1 : ip[i];
			}
			values[mask] = math::dot(grad_map_(ii), tt);
		}

		auto smoothstep = [](auto x0, auto x1, auto t)
		{
			auto const s = t * t * (3 - 2 * t);
			return x0 * (1 - s) + x1 * s;
		};

		for (std::size_t i = N; i --> 0;)
		{
			for (std::size_t mask = 0; mask < (1 << i); ++mask)
				values[mask] = smoothstep(values[mask], values[mask | (1 << i)], t[i]);
		}

		auto const M = std::sqrt(T(N) / T(4));

		return 0.5 * (1. + values[0] / M);
	}

}
