diff --git a/libs/pcg/include/psemek/pcg/perlin.hpp b/libs/pcg/include/psemek/pcg/perlin.hpp
index 57c940da..96b3fc09 100644
--- a/libs/pcg/include/psemek/pcg/perlin.hpp
+++ b/libs/pcg/include/psemek/pcg/perlin.hpp
@@ -9,11 +9,15 @@
 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(gfx::basic_pixmap<geom::vector<float, 2>> grad_map);
-		perlin(gfx::basic_pixmap<geom::vector<float, 2>> const & grad_map, seamless_tag);
+		perlin(util::array<geom::vector<T, N>, N> grad_map);
+		perlin(util::array<geom::vector<T, N>, N> const & grad_map, seamless_tag);
 		perlin(perlin &&) = default;
 
 		perlin & operator = (perlin &&) = default;
@@ -28,35 +32,90 @@ namespace psemek::pcg
 			return grad_map_.height() - 1;
 		}
 
-		// x \in [0.0 .. 1.0]
-		// y \in [0.0 .. 1.0]
-		float operator()(float x, float y) const;
+		std::size_t depth() const
+		{
+			return grad_map_.depth() - 1;
+		}
+
+		// Coords \in [0.0 .. 1.0]
+		T operator()(geom::vector<T, N> p) const;
 
 		auto & grad() { return grad_map_; }
 		auto const & grad() const { return grad_map_; }
 
 	private:
-		gfx::basic_pixmap<geom::vector<float, 2>> grad_map_;
+		util::array<geom::vector<T, N>, N> grad_map_;
 	};
 
-	struct fractal_perlin
+	template <typename T, std::size_t N>
+	perlin<T, N>::perlin(util::array<geom::vector<T, N>, N> grad_map)
+		: grad_map_{std::move(grad_map)}
+	{}
+
+	template <typename T, std::size_t N>
+	perlin<T, N>::perlin(util::array<geom::vector<T, N>, N> const & grad_map, seamless_tag)
 	{
-		fractal_perlin() = default;
-		fractal_perlin(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> grad_maps);
-		fractal_perlin(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> grad_maps, std::vector<float> weights);
-		fractal_perlin(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> const & grad_maps, seamless_tag);
-		fractal_perlin(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> const & grad_maps, std::vector<float> weights, seamless_tag);
-		fractal_perlin(fractal_perlin &&) = default;
+		auto dims = grad_map.dims();
+		for (std::size_t i = 0; i < N; ++i)
+			dims[i] += 1;
 
-		fractal_perlin & operator = (fractal_perlin &&) = default;
+		grad_map_.resize(dims);
 
-		// x \in [0.0 .. 1.0]
-		// y \in [0.0 .. 1.0]
-		float operator()(float x, float y) const;
+		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);
+		}
+	}
 
-	private:
-		std::vector<perlin> octaves_;
-		std::vector<float> weights_;
-	};
+	template <typename T, std::size_t N>
+	T perlin<T, N>::operator()(geom::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;
+		}
+
+		geom::vector<std::size_t, N> ip;
+		for (std::size_t i = 0; i < N; ++i)
+		{
+			ip[i] = geom::clamp<std::size_t>(std::floor(p[i]), {0, static_cast<std::size_t>(grad_map_.dim(i)) - 2});
+		}
+
+		geom::vector<float, N> t = p - geom::cast<float>(ip);
+
+		T values[1 << N];
+
+		for (std::size_t mask = 0; mask < (1 << N); ++mask)
+		{
+			geom::vector<float, 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] = geom::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]);
+		}
+
+		return 0.5f * (1.f + std::sqrt(2.f) * values[0]);
+	}
 
 }
diff --git a/libs/pcg/source/perlin.cpp b/libs/pcg/source/perlin.cpp
deleted file mode 100644
index 8aac4cb2..00000000
--- a/libs/pcg/source/perlin.cpp
+++ /dev/null
@@ -1,131 +0,0 @@
-#include <psemek/pcg/perlin.hpp>
-
-#include <psemek/geom/interval.hpp>
-
-#include <psemek/util/assert.hpp>
-
-namespace psemek::pcg
-{
-
-	perlin::perlin(gfx::basic_pixmap<geom::vector<float, 2>> grad_map)
-		: grad_map_(std::move(grad_map))
-	{}
-
-	perlin::perlin(gfx::basic_pixmap<geom::vector<float, 2>> const & grad_map, seamless_tag)
-	{
-		auto const w = grad_map.width();
-		auto const h = grad_map.height();
-
-		grad_map_.resize({w + 1, h + 1});
-
-		for (std::size_t j = 0; j <= h; ++j)
-		{
-			for (std::size_t i = 0; i <= w; ++i)
-			{
-				grad_map_(i, j) = grad_map(i % w, j % h);
-			}
-		}
-	}
-
-	static float step(float x0, float x1, float t)
-	{
-		return x0 * (1.f - t) + x1 * t;
-	}
-
-	static float smoothstep(float x0, float x1, float t)
-	{
-		float const s = t * t * (3.f - 2.f * t);
-		return step(x0, x1, s);
-	}
-
-	float perlin::operator()(float x, float y) const
-	{
-		assert(x >= 0.f);
-		assert(y >= 0.f);
-
-		assert(x <= 1.f);
-		assert(y <= 1.f);
-
-		x *= width();
-		y *= height();
-
-		int const ix = geom::clamp<int>(std::floor(x), {0, static_cast<int>(width()) - 1});
-		int const iy = geom::clamp<int>(std::floor(y), {0, static_cast<int>(height()) - 1});
-
-		float const tx = x - ix;
-		float const ty = y - iy;
-
-		float const d00 = tx * grad_map_(ix, iy)[0] + ty * grad_map_(ix, iy)[1];
-		float const d10 = (tx-1.f) * grad_map_(ix+1, iy)[0] + ty * grad_map_(ix+1, iy)[1];
-		float const d01 = tx * grad_map_(ix, iy+1)[0] + (ty-1.f) * grad_map_(ix, iy+1)[1];
-		float const d11 = (tx-1.f) * grad_map_(ix+1, iy+1)[0] + (ty-1.f) * grad_map_(ix+1, iy+1)[1];
-
-		float v = smoothstep(smoothstep(d00, d10, tx), smoothstep(d01, d11, tx), ty) * std::sqrt(2.f);
-		return 0.5f + v * 0.5f;
-	}
-
-	static std::vector<float> estimate_weights(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> const & grad_maps)
-	{
-		std::vector<float> weights(grad_maps.size());
-
-		float sum = 0.f;
-		for (std::size_t i = 0; i < grad_maps.size(); ++i)
-		{
-			auto const & m = grad_maps[i];
-			weights[i] = 1.f / std::sqrt(m.width() * m.width() + m.height() * m.height());
-			sum += weights[i];
-		}
-
-		for (std::size_t i = 0; i < grad_maps.size(); ++i)
-			weights[i] /= sum;
-
-		return weights;
-	}
-
-	fractal_perlin::fractal_perlin(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> grad_maps)
-	{
-		weights_ = estimate_weights(grad_maps);
-
-		octaves_.resize(grad_maps.size());
-		for (std::size_t i = 0; i < grad_maps.size(); ++i)
-			octaves_[i] = perlin(std::move(grad_maps[i]));
-	}
-
-	fractal_perlin::fractal_perlin(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> grad_maps, std::vector<float> weights)
-		: weights_{std::move(weights)}
-	{
-		assert(weights_.size() == grad_maps.size());
-
-		octaves_.resize(grad_maps.size());
-		for (std::size_t i = 0; i < grad_maps.size(); ++i)
-			octaves_[i] = perlin(std::move(grad_maps[i]));
-	}
-
-	fractal_perlin::fractal_perlin(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> const & grad_maps, seamless_tag)
-	{
-		weights_ = estimate_weights(grad_maps);
-
-		octaves_.resize(grad_maps.size());
-		for (std::size_t i = 0; i < grad_maps.size(); ++i)
-			octaves_[i] = perlin(grad_maps[i], seamless);
-	}
-
-	fractal_perlin::fractal_perlin(std::vector<gfx::basic_pixmap<geom::vector<float, 2>>> const & grad_maps, std::vector<float> weights, seamless_tag)
-		: weights_{std::move(weights)}
-	{
-		assert(weights_.size() == grad_maps.size());
-
-		octaves_.resize(grad_maps.size());
-		for (std::size_t i = 0; i < grad_maps.size(); ++i)
-			octaves_[i] = perlin(grad_maps[i], seamless);
-	}
-
-	float fractal_perlin::operator()(float x, float y) const
-	{
-		float v = 0.f;
-		for (std::size_t i = 0; i < octaves_.size(); ++i)
-			v += weights_[i] * octaves_[i](x, y);
-		return v;
-	}
-
-}