Add lazy perlin noise generator

libs/pcg/include/psemek/pcg/lazy_perlin.hpp

@@ -0,0 +1,85 @@
+#pragma once
+
+#include <psemek/util/function.hpp>
+
+#include <psemek/geom/vector.hpp>
+
+#include <unordered_map>
+
+namespace psemek::pcg
+{
+
+	template <typename T, std::size_t N>
+	struct lazy_perlin
+	{
+		using generator_func = util::function<geom::vector<T, N>(geom::vector<int, N> const &)>;
+
+		lazy_perlin (std::size_t grid_size, generator_func generator)
+			: grid_size_(grid_size)
+			, gen_(std::move(generator))
+		{}
+
+		template <typename ... Args>
+		T operator()(Args const & ... args) const
+		{
+			return (*this)(geom::vector<T, N>{args...});
+		}
+
+		T operator()(geom::vector<T, N> p) const;
+
+	private:
+		int const grid_size_;
+		mutable std::unordered_map<geom::vector<int, N>, geom::vector<T, N>> grid_;
+		generator_func gen_;
+
+		geom::vector<T, N> grid_at(geom::vector<int, N> const & c) const
+		{
+			auto it = grid_.find(c);
+			if (it == grid_.end())
+				it = grid_.insert({c, gen_(c)}).first;
+			return it->second;
+		}
+	};
+
+	template <typename T, std::size_t N>
+	T lazy_perlin<T, N>::operator()(geom::vector<T, N> p) const
+	{
+		geom::vector<int, N> ip;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			ip[i] = (p[i] >= 0) ? (p[i] / grid_size_) : -((- p[i] + grid_size_ - 1) / grid_size_);
+		}
+
+		geom::vector<T, N> t = (p - geom::cast<T>(ip) * T(grid_size_)) / T(grid_size_);
+
+		T values[1 << N];
+
+		for (std::size_t mask = 0; mask < (1 << N); ++mask)
+		{
+			geom::vector<T, N> tt;
+			geom::vector<int, 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] = geom::dot(grid_at(ii), tt);
+		}
+
+		auto smootherstep = [](float x0, float x1, float t)
+		{
+			auto const s = t * t * t * (10. + t * (-15. + 6. * 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] = smootherstep(values[mask], values[mask | (1 << i)], t[i]);
+		}
+
+		return 0.5 * (1. + std::sqrt(2.) * values[0]);
+	}
+
+}