diff --git a/examples/cloud.cpp b/examples/cloud.cpp
index 274d5b3c..adec7a4f 100644
--- a/examples/cloud.cpp
+++ b/examples/cloud.cpp
@@ -223,7 +223,7 @@ struct cloud_app
 			random::generator rng(random::device{});
 			random::uniform_sphere_vector_distribution<float, 3> d;
 
-			std::vector<util::array<math::vector<float, 3>, 3>> grad(4);
+			std::vector<util::ndarray<math::vector<float, 3>, 3>> grad(4);
 			std::vector<float> weights(grad.size());
 			float weight_sum = 0.f;
 
@@ -247,7 +247,7 @@ struct cloud_app
 			// https://gcc.gnu.org/bugzilla/show_bug.cgi?id=89062
 			math::gradient<float> g(std::make_pair(0.2f, 0.f), math::easing_type::quadratic_out, std::pair{0.3f, max_density});
 
-			util::array<std::uint8_t, 3> cloud_data({cloud_size[0], cloud_size[1], cloud_size[2]});
+			util::ndarray<std::uint8_t, 3> cloud_data({cloud_size[0], cloud_size[1], cloud_size[2]});
 
 			for (std::size_t z = 0; z < cloud_data.depth(); ++z)
 			{
@@ -328,7 +328,7 @@ struct cloud_app
 				return math::lerp(v0, v1, t[2]) / 255.f * max_density;
 			};
 
-			util::array<math::vector<float, 4>, 3> cloud_shadow_f(cloud_data.dims());
+			util::ndarray<math::vector<float, 4>, 3> cloud_shadow_f(cloud_data.dims());
 
 			dirs.resize(32);
 
@@ -448,7 +448,7 @@ struct cloud_app
 				}
 			}
 
-			util::array<math::vector<std::uint8_t, 4>, 3> cloud_shadow(cloud_data.dims());
+			util::ndarray<math::vector<std::uint8_t, 4>, 3> cloud_shadow(cloud_data.dims());
 
 			for (auto const & idx : cloud_shadow.indices())
 			{
diff --git a/examples/fibonacci_music_box.cpp b/examples/fibonacci_music_box.cpp
index 176bff5a..bb9f638e 100644
--- a/examples/fibonacci_music_box.cpp
+++ b/examples/fibonacci_music_box.cpp
@@ -22,7 +22,7 @@
 #include <psemek/app/application_base.hpp>
 #include <psemek/app/default_application_factory.hpp>
 #include <psemek/log/log.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 
 #include <unordered_map>
 
@@ -32,7 +32,7 @@ using namespace psemek;
 
 std::vector<std::vector<math::point<int, 2>>> fibonacci_cycles(int n)
 {
-	util::array<int, 2> cycle({n, n}, -1);
+	util::ndarray<int, 2> cycle({n, n}, -1);
 
 	std::vector<std::vector<math::point<int, 2>>> result;
 
diff --git a/examples/fire.cpp b/examples/fire.cpp
index db1b350b..62a2da16 100644
--- a/examples/fire.cpp
+++ b/examples/fire.cpp
@@ -172,7 +172,7 @@ candle_renderer::candle_renderer()
 		random::generator rng;
 		random::uniform_sphere_vector_distribution<float, 2> d;
 
-		util::array<math::vector<float, 2>, 2> grad({16, 16});
+		util::ndarray<math::vector<float, 2>, 2> grad({16, 16});
 
 		for (auto & v : grad) v = d(rng);
 		pcg::perlin<float, 2> perlinx(grad, pcg::seamless);
diff --git a/examples/particle_life_2d.cpp b/examples/particle_life_2d.cpp
index 2b27fa02..85f74e4a 100644
--- a/examples/particle_life_2d.cpp
+++ b/examples/particle_life_2d.cpp
@@ -10,7 +10,7 @@
 #include <psemek/random/uniform.hpp>
 #include <psemek/random/uniform_box.hpp>
 #include <psemek/util/clock.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/log/log.hpp>
 
 #include <format>
@@ -418,9 +418,9 @@ private:
 
 	int types_;
 	std::vector<gfx::color_rgba> colors_;
-	util::array<float, 2> force_constants_;
-	util::array<float, 2> force_distance_;
-	util::array<float, 2> collision_distance_;
+	util::ndarray<float, 2> force_constants_;
+	util::ndarray<float, 2> force_distance_;
+	util::ndarray<float, 2> collision_distance_;
 	std::vector<float> max_force_distance_;
 
 	std::vector<particle> particles_;
diff --git a/examples/physics.cpp b/examples/physics.cpp
index 6a80550d..918e90e0 100644
--- a/examples/physics.cpp
+++ b/examples/physics.cpp
@@ -507,7 +507,7 @@ void physics_demo_app::update()
 
 		std::vector<util::fixed_vector<std::array<int, 2>, 8>> ball_cells(model.points.size());
 
-		util::array<util::fixed_vector<std::uint32_t, 16>, 2> cells({std::ceil(view_region[0].length() / cell_size), std::ceil(view_region[1].length() / cell_size)});
+		util::ndarray<util::fixed_vector<std::uint32_t, 16>, 2> cells({std::ceil(view_region[0].length() / cell_size), std::ceil(view_region[1].length() / cell_size)});
 
 		float dx = view_region[0].length() / cells.width();
 		float dy = view_region[1].length() / cells.height();
diff --git a/examples/srtm.cpp b/examples/srtm.cpp
index cec65db4..1adc1bc6 100644
--- a/examples/srtm.cpp
+++ b/examples/srtm.cpp
@@ -648,7 +648,7 @@ srtm_app::srtm_app(options const &, context const & context)
 	selected_mesh.setup<math::point<float, 3>>();
 
 	{
-		util::array<gfx::color_rgb, 1> colors({16});
+		util::ndarray<gfx::color_rgb, 1> colors({16});
 
 		auto * c = colors.data();
 
@@ -674,7 +674,7 @@ srtm_app::srtm_app(options const &, context const & context)
 	}
 
 	{
-		util::array<gfx::color_rgb, 1> colors({5});
+		util::ndarray<gfx::color_rgb, 1> colors({5});
 
 		colors(0) = {0, 63, 127};
 		colors(1) = {0, 0, 127};
diff --git a/examples/water_2d.cpp b/examples/water_2d.cpp
index 9ace87fb..f9081d00 100644
--- a/examples/water_2d.cpp
+++ b/examples/water_2d.cpp
@@ -4,7 +4,7 @@
 #include <psemek/math/orthographic.hpp>
 #include <psemek/math/camera.hpp>
 #include <psemek/log/log.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/random/device.hpp>
 #include <psemek/random/generator.hpp>
 #include <psemek/random/uniform.hpp>
@@ -36,13 +36,13 @@ struct water_2d_app
 
 	float time = 0.f;
 
-	util::array<float, 2> bed;
-	util::array<float, 2> water;
-	util::array<float, 2> flowx;
-	util::array<float, 2> flowy;
-	util::array<math::vector<float, 2>, 2> velocity;
-	util::array<float, 2> sediment;
-	util::array<float, 2> new_sediment;
+	util::ndarray<float, 2> bed;
+	util::ndarray<float, 2> water;
+	util::ndarray<float, 2> flowx;
+	util::ndarray<float, 2> flowy;
+	util::ndarray<math::vector<float, 2>, 2> velocity;
+	util::ndarray<float, 2> sediment;
+	util::ndarray<float, 2> new_sediment;
 
 	gfx::painter painter;
 
@@ -83,7 +83,7 @@ water_2d_app::water_2d_app(options const &, context const & ctx)
 	// water(N / 2, N / 2) = 100.f;
 
 	random::uniform_sphere_vector_distribution<float, 2> d;
-	util::array<math::vector<float, 2>, 2> grads({9, 9});
+	util::ndarray<math::vector<float, 2>, 2> grads({9, 9});
 	for (auto & v : grads)
 		v = d(rng);
 
diff --git a/examples/water_2d_hex.cpp b/examples/water_2d_hex.cpp
index 80badaf5..08f43d40 100644
--- a/examples/water_2d_hex.cpp
+++ b/examples/water_2d_hex.cpp
@@ -5,7 +5,7 @@
 #include <psemek/math/box.hpp>
 #include <psemek/math/camera.hpp>
 #include <psemek/math/gauss.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/random/generator.hpp>
 #include <psemek/random/device.hpp>
 #include <psemek/random/uniform_sphere.hpp>
@@ -66,11 +66,11 @@ private:
 
 	bool show_velocity_ = false;
 
-	util::array<float, 2> bed_;
-	util::array<float, 2> water_;
-	util::array<float, 2> flow_x_; // flow_x(x, y) is (x-1, y) => (x, y)
-	util::array<float, 2> flow_y_; // flow_y(x, y) is (x, y-1) => (x, y)
-	util::array<float, 2> flow_z_; // flow_z(x, y) is (x, y-1) => (x-1, y)
+	util::ndarray<float, 2> bed_;
+	util::ndarray<float, 2> water_;
+	util::ndarray<float, 2> flow_x_; // flow_x(x, y) is (x-1, y) => (x, y)
+	util::ndarray<float, 2> flow_y_; // flow_y(x, y) is (x, y-1) => (x, y)
+	util::ndarray<float, 2> flow_z_; // flow_z(x, y) is (x, y-1) => (x-1, y)
 };
 
 water_2d_hex_app::water_2d_hex_app(options const &, context const & ctx)
@@ -86,7 +86,7 @@ water_2d_hex_app::water_2d_hex_app(options const &, context const & ctx)
 
 	random::uniform_sphere_vector_distribution<float, 2> d_grad;
 
-	util::array<math::vector<float, 2>, 2> perlin_grad({17, 17});
+	util::ndarray<math::vector<float, 2>, 2> perlin_grad({17, 17});
 	for (auto & v : perlin_grad)
 		v = d_grad(rng_);
 
diff --git a/libs/fonts/source/freetype.cpp b/libs/fonts/source/freetype.cpp
index 30795f96..38063bd3 100644
--- a/libs/fonts/source/freetype.cpp
+++ b/libs/fonts/source/freetype.cpp
@@ -145,7 +145,7 @@ namespace psemek::fonts
 
 			struct page
 			{
-				util::array<std::uint8_t, 2> pixmap;
+				util::ndarray<std::uint8_t, 2> pixmap;
 
 				int current_row_start = 0;
 				int current_row_height = 0;
diff --git a/libs/gfx/include/psemek/gfx/pixmap.hpp b/libs/gfx/include/psemek/gfx/pixmap.hpp
index 00566687..524d28bc 100644
--- a/libs/gfx/include/psemek/gfx/pixmap.hpp
+++ b/libs/gfx/include/psemek/gfx/pixmap.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/gfx/color.hpp>
 #include <psemek/io/stream.hpp>
 
@@ -10,7 +10,7 @@ namespace psemek::gfx
 {
 
 	template <typename Pixel>
-	using basic_pixmap = util::array<Pixel, 2>;
+	using basic_pixmap = util::ndarray<Pixel, 2>;
 
 	using pixmap_monochrome = basic_pixmap<std::uint8_t>;
 	using pixmap_rgb = basic_pixmap<color_rgb>;
@@ -48,7 +48,7 @@ namespace psemek::gfx
 	// Utilities
 
 	template <typename Pixel, std::size_t N>
-	auto to_srgb(util::array<Pixel, N> pm, float g = 1.f / 2.2f)
+	auto to_srgb(util::ndarray<Pixel, N> pm, float g = 1.f / 2.2f)
 	{
 		for (auto & c : pm)
 			c = to_srgb(c, g);
@@ -56,7 +56,7 @@ namespace psemek::gfx
 	}
 
 	template <typename Pixel, std::size_t N>
-	auto to_linear(util::array<Pixel, N> pm, float g = 1.f / 2.2f)
+	auto to_linear(util::ndarray<Pixel, N> pm, float g = 1.f / 2.2f)
 	{
 		return to_srgb(std::move(pm), 1.f / g);
 	}
diff --git a/libs/gfx/include/psemek/gfx/texture.hpp b/libs/gfx/include/psemek/gfx/texture.hpp
index 88cb4f09..de20fa8f 100644
--- a/libs/gfx/include/psemek/gfx/texture.hpp
+++ b/libs/gfx/include/psemek/gfx/texture.hpp
@@ -54,10 +54,10 @@ namespace psemek::gfx
 		void load_srgb(math::vector<std::size_t, D> const & size, Pixel const * data = nullptr);
 
 		template <typename Pixel>
-		void load(util::array<Pixel, D> const & p);
+		void load(util::ndarray<Pixel, D> const & p);
 
 		template <typename Pixel>
-		void load_srgb(util::array<Pixel, D> const & p);
+		void load_srgb(util::ndarray<Pixel, D> const & p);
 
 #ifndef PSEMEK_GLES
 		void pixels(GLenum format, GLenum type, void * data, int layer = 0) const;
@@ -142,7 +142,7 @@ namespace psemek::gfx
 		}
 
 		template <typename Pixel>
-		void load(int f, util::array<Pixel, 2> const & p)
+		void load(int f, util::ndarray<Pixel, 2> const & p)
 		{
 			math::vector<std::size_t, 2> size;
 			for (std::size_t i = 0; i < 2; ++i) size[i] = p.dim(i);
@@ -392,7 +392,7 @@ namespace psemek::gfx
 
 	template <std::size_t D, GLenum Target>
 	template <typename Pixel>
-	void basic_texture<D, Target>::load(util::array<Pixel, D> const & p)
+	void basic_texture<D, Target>::load(util::ndarray<Pixel, D> const & p)
 	{
 		math::vector<std::size_t, D> size;
 		for (std::size_t i = 0; i < D; ++i) size[i] = p.dim(i);
@@ -401,7 +401,7 @@ namespace psemek::gfx
 
 	template <std::size_t D, GLenum Target>
 	template <typename Pixel>
-	void basic_texture<D, Target>::load_srgb(util::array<Pixel, D> const & p)
+	void basic_texture<D, Target>::load_srgb(util::ndarray<Pixel, D> const & p)
 	{
 		math::vector<std::size_t, D> size;
 		for (std::size_t i = 0; i < D; ++i) size[i] = p.dim(i);
diff --git a/libs/gfx/include/psemek/gfx/texture_atlas.hpp b/libs/gfx/include/psemek/gfx/texture_atlas.hpp
index b46d81b4..0060ba42 100644
--- a/libs/gfx/include/psemek/gfx/texture_atlas.hpp
+++ b/libs/gfx/include/psemek/gfx/texture_atlas.hpp
@@ -19,9 +19,9 @@ namespace psemek::gfx
 		texture_2d & texture() { return texture_; }
 		texture_2d const & texture() const { return texture_; }
 
-		util::array<Pixel, 2> const & pixmap() const { return atlas_.array(); }
+		util::ndarray<Pixel, 2> const & pixmap() const { return atlas_.array(); }
 
-		std::pair<texture_view, bool> insert(Key const & key, util::array<Pixel, 2> const & data);
+		std::pair<texture_view, bool> insert(Key const & key, util::ndarray<Pixel, 2> const & data);
 		texture_view find(Key const & key) const;
 		texture_view at(Key const & key) const;
 
@@ -43,7 +43,7 @@ namespace psemek::gfx
 	{}
 
 	template <typename Pixel, typename Key, typename Compare>
-	std::pair<texture_view_2d, bool> texture_atlas_2d<Pixel, Key, Compare>::insert(Key const & key, util::array<Pixel, 2> const & data)
+	std::pair<texture_view_2d, bool> texture_atlas_2d<Pixel, Key, Compare>::insert(Key const & key, util::ndarray<Pixel, 2> const & data)
 	{
 		auto result = atlas_.insert(key, data);
 		if (result.second)
diff --git a/libs/gfx/source/renderer/deferred.cpp b/libs/gfx/source/renderer/deferred.cpp
index 03d6be49..f96a9810 100644
--- a/libs/gfx/source/renderer/deferred.cpp
+++ b/libs/gfx/source/renderer/deferred.cpp
@@ -1069,7 +1069,7 @@ void main(){}
 			camera_distance[i] = dist;
 		}
 
-		util::array<bin, 3> bins({opts.grid_size[0], opts.grid_size[1], opts.grid_size[2]});
+		util::ndarray<bin, 3> bins({opts.grid_size[0], opts.grid_size[1], opts.grid_size[2]});
 
 		math::box<float, 3> lit_bbox;
 		math::box<float, 3> casts_shadow_bbox;
diff --git a/libs/pcg/include/psemek/pcg/chunked_map.hpp b/libs/pcg/include/psemek/pcg/chunked_map.hpp
index a20a6e16..96ef5c17 100644
--- a/libs/pcg/include/psemek/pcg/chunked_map.hpp
+++ b/libs/pcg/include/psemek/pcg/chunked_map.hpp
@@ -1,7 +1,7 @@
 #pragma once
 
 #include <psemek/util/function.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/util/hash_table.hpp>
 #include <psemek/math/box.hpp>
 
@@ -13,7 +13,7 @@ namespace psemek::pcg
 	template <typename T, std::size_t N>
 	struct chunked_map
 	{
-		using generator_func = util::function<util::array<T, N>(math::box<int, N> const &)>;
+		using generator_func = util::function<util::ndarray<T, N>(math::box<int, N> const &)>;
 
 		chunked_map(std::size_t chunk_size, generator_func generator)
 			: chunk_size_(chunk_size)
@@ -68,7 +68,7 @@ namespace psemek::pcg
 	private:
 		int const chunk_size_;
 		generator_func generator_;
-		mutable util::hash_map<math::vector<int, N>, util::array<T, N>> chunks_;
+		mutable util::hash_map<math::vector<int, N>, util::ndarray<T, N>> chunks_;
 
 	};
 
diff --git a/libs/pcg/include/psemek/pcg/lazy_perlin.hpp b/libs/pcg/include/psemek/pcg/lazy_perlin.hpp
index ec97b83a..9abbcbea 100644
--- a/libs/pcg/include/psemek/pcg/lazy_perlin.hpp
+++ b/libs/pcg/include/psemek/pcg/lazy_perlin.hpp
@@ -1,7 +1,7 @@
 #pragma once
 
 #include <psemek/util/function.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/util/hash_table.hpp>
 
 #include <psemek/math/vector.hpp>
@@ -68,7 +68,7 @@ namespace psemek::pcg
 
 		lazy_perlin_view() = default;
 
-		lazy_perlin_view(std::size_t grid_size, math::vector<int, N> const & view_origin, util::array<math::vector<T, N>, N> subgrid)
+		lazy_perlin_view(std::size_t grid_size, math::vector<int, N> const & view_origin, util::ndarray<math::vector<T, N>, N> subgrid)
 			: grid_size_(grid_size)
 			, origin_(view_origin)
 			, subgrid_(std::move(subgrid))
@@ -79,7 +79,7 @@ namespace psemek::pcg
 	private:
 		int grid_size_;
 		math::vector<int, N> origin_;
-		util::array<math::vector<T, N>, N> subgrid_;
+		util::ndarray<math::vector<T, N>, N> subgrid_;
 	};
 
 	template <typename T, std::size_t N>
@@ -143,7 +143,7 @@ namespace psemek::pcg
 		int const grid_size_;
 		math::vector<int, N> origin_;
 		math::vector<int, N> corner_;
-		mutable util::array<std::optional<math::vector<T, N>>, N> grid_;
+		mutable util::ndarray<std::optional<math::vector<T, N>>, N> grid_;
 		generator_func gen_;
 
 		math::vector<T, N> grid_at(math::vector<int, N> const & c) const
@@ -191,7 +191,7 @@ namespace psemek::pcg
 			subgrid_dimensions[i] = subgrid_max - subgrid_origin[i] + 1;
 		}
 
-		util::array<math::vector<T, N>, N> subgrid(subgrid_dimensions);
+		util::ndarray<math::vector<T, N>, N> subgrid(subgrid_dimensions);
 		for (auto const idx : subgrid.indices())
 		{
 			math::vector<int, N> id;
diff --git a/libs/pcg/include/psemek/pcg/perlin.hpp b/libs/pcg/include/psemek/pcg/perlin.hpp
index 374b6bd5..9414b8d5 100644
--- a/libs/pcg/include/psemek/pcg/perlin.hpp
+++ b/libs/pcg/include/psemek/pcg/perlin.hpp
@@ -16,8 +16,8 @@ namespace psemek::pcg
 		static constexpr std::size_t dimension = N;
 
 		perlin() = default;
-		perlin(util::array<math::vector<T, N>, N> grad_map);
-		perlin(util::array<math::vector<T, N>, N> const & grad_map, seamless_tag);
+		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;
@@ -51,16 +51,16 @@ namespace psemek::pcg
 		auto const & grad() const { return grad_map_; }
 
 	private:
-		util::array<math::vector<T, N>, N> grad_map_;
+		util::ndarray<math::vector<T, N>, N> grad_map_;
 	};
 
 	template <typename T, std::size_t N>
-	perlin<T, N>::perlin(util::array<math::vector<T, N>, N> grad_map)
+	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::array<math::vector<T, N>, N> const & grad_map, seamless_tag)
+	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)
diff --git a/libs/pcg/include/psemek/pcg/sample.hpp b/libs/pcg/include/psemek/pcg/sample.hpp
index 659db515..95f27750 100644
--- a/libs/pcg/include/psemek/pcg/sample.hpp
+++ b/libs/pcg/include/psemek/pcg/sample.hpp
@@ -21,7 +21,7 @@ namespace psemek::pcg
 	}
 
 	template <typename Gen, typename R>
-	void sample(Gen && gen, util::array<R, 3> & result)
+	void sample(Gen && gen, util::ndarray<R, 3> & result)
 	{
 		for (std::size_t k = 0; k < result.height(); ++k)
 		{
@@ -55,7 +55,7 @@ namespace psemek::pcg
 	{
 		using R = decltype(gen(0.f, 0.f, 0.f));
 
-		util::array<R, 3> result({width, height, depth});
+		util::ndarray<R, 3> result({width, height, depth});
 
 		sample(std::forward<Gen>(gen), result);
 		return result;
diff --git a/libs/phys/include/psemek/phys/shallow_water.hpp b/libs/phys/include/psemek/phys/shallow_water.hpp
index ff457891..a7b1d6ea 100644
--- a/libs/phys/include/psemek/phys/shallow_water.hpp
+++ b/libs/phys/include/psemek/phys/shallow_water.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/math/math.hpp>
 #include <psemek/math/dual.hpp>
 
@@ -8,16 +8,16 @@ namespace psemek::phys
 {
 
 	template <typename T>
-	std::pair<std::size_t, T> shallow_water_steady_state_single_step(util::array<T, 2> const & depth, util::array<T, 2> & water_height,
-		util::array<T, 2> & x_velocity, util::array<T, 2> & y_velocity,
+	std::pair<std::size_t, T> shallow_water_steady_state_single_step(util::ndarray<T, 2> const & depth, util::ndarray<T, 2> & water_height,
+		util::ndarray<T, 2> & x_velocity, util::ndarray<T, 2> & y_velocity,
 		T gravity, std::size_t max_iterations, T max_error)
 	{
 		std::size_t const n = depth.width() - 2;
 		std::size_t const m = depth.height() - 2;
 
-		util::array<T, 2> delta_height({n + 2, m + 2}, 0.f);
-		util::array<T, 2> delta_x_velocity({n + 1, m + 2}, 0.f);
-		util::array<T, 2> delta_y_velocity({n + 2, m + 1}, 0.f);
+		util::ndarray<T, 2> delta_height({n + 2, m + 2}, 0.f);
+		util::ndarray<T, 2> delta_x_velocity({n + 1, m + 2}, 0.f);
+		util::ndarray<T, 2> delta_y_velocity({n + 2, m + 1}, 0.f);
 
 		auto h_dual = [&](auto i, auto j) { return math::dual<float, 2>{water_height(i, j), {delta_height(i, j), 0.f}}; };
 		auto ux_dual = [&](auto i, auto j) { return math::dual<float, 2>{x_velocity(i, j + 1), {delta_x_velocity(i, j + 1), 0.f}}; };
diff --git a/libs/sir/include/psemek/sir/array.hpp b/libs/sir/include/psemek/sir/array.hpp
index 42c7a797..4a49e409 100644
--- a/libs/sir/include/psemek/sir/array.hpp
+++ b/libs/sir/include/psemek/sir/array.hpp
@@ -2,14 +2,14 @@
 
 #include <psemek/sir/trivial.hpp>
 #include <psemek/sir/container.hpp>
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 
 namespace psemek::sir
 {
 
 	template <typename Stream, typename T, std::size_t N>
-	requires (is_ostream_v<Stream> && !is_custom_v<util::array<T, N>>)
-	void write(Stream & s, util::array<T, N> const & x)
+	requires (is_ostream_v<Stream> && !is_custom_v<util::ndarray<T, N>>)
+	void write(Stream & s, util::ndarray<T, N> const & x)
 	{
 		write(s, x.dims());
 		if constexpr (is_trivial_v<T>)
@@ -23,8 +23,8 @@ namespace psemek::sir
 	}
 
 	template <typename Stream, typename T, std::size_t N>
-	requires (is_istream_v<Stream> && !is_custom_v<util::array<T, N>>)
-	void read(Stream & s, util::array<T, N> & x)
+	requires (is_istream_v<Stream> && !is_custom_v<util::ndarray<T, N>>)
+	void read(Stream & s, util::ndarray<T, N> & x)
 	{
 		std::array<std::size_t, N> dims;
 		read(s, dims);
diff --git a/libs/ui_legacy/include/psemek/ui/grid_layout.hpp b/libs/ui_legacy/include/psemek/ui/grid_layout.hpp
index f7b7c4b7..fc7233e0 100644
--- a/libs/ui_legacy/include/psemek/ui/grid_layout.hpp
+++ b/libs/ui_legacy/include/psemek/ui/grid_layout.hpp
@@ -4,7 +4,7 @@
 #include <psemek/ui/container.hpp>
 #include <psemek/ui/box_shape.hpp>
 
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 
 namespace psemek::ui
 {
@@ -60,7 +60,7 @@ namespace psemek::ui
 		std::vector<math::interval<float>> column_shape_;
 
 	private:
-		util::array<std::shared_ptr<element>, 2> children_;
+		util::ndarray<std::shared_ptr<element>, 2> children_;
 		std::vector<element *> children_range_;
 
 		std::vector<float> row_weight_;
diff --git a/libs/ui_legacy/source/grid_layout.cpp b/libs/ui_legacy/source/grid_layout.cpp
index b35959c7..fb325530 100644
--- a/libs/ui_legacy/source/grid_layout.cpp
+++ b/libs/ui_legacy/source/grid_layout.cpp
@@ -159,7 +159,7 @@ namespace psemek::ui
 	}
 
 	template <int Dimension>
-	static std::vector<math::interval<float>> cell_sizes(util::array<std::shared_ptr<element>, 2> const & children)
+	static std::vector<math::interval<float>> cell_sizes(util::ndarray<std::shared_ptr<element>, 2> const & children)
 	{
 		std::vector<math::interval<float>> result(children.dim(Dimension), {0.f, inf});
 
@@ -175,7 +175,7 @@ namespace psemek::ui
 	}
 
 	template <int Dimension>
-	static std::vector<math::interval<float>> cell_sizes_constrained(util::array<std::shared_ptr<element>, 2> const & children, std::vector<math::interval<float>> const & other_dimension_shapes)
+	static std::vector<math::interval<float>> cell_sizes_constrained(util::ndarray<std::shared_ptr<element>, 2> const & children, std::vector<math::interval<float>> const & other_dimension_shapes)
 	{
 		std::vector<math::interval<float>> result(children.dim(Dimension), {0.f, inf});
 
diff --git a/libs/util/include/psemek/util/atlas.hpp b/libs/util/include/psemek/util/atlas.hpp
index 5c9c0076..c5fe9765 100644
--- a/libs/util/include/psemek/util/atlas.hpp
+++ b/libs/util/include/psemek/util/atlas.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 
 #include <map>
 #include <optional>
@@ -33,12 +33,12 @@ namespace psemek::util
 			, data_(std::move(compare))
 		{}
 
-		util::array<T, 2> const & array() const
+		util::ndarray<T, 2> const & array() const
 		{
 			return array_;
 		}
 
-		std::pair<atlas_part, bool> insert(Key const & key, util::array<T, 2> const & data)
+		std::pair<atlas_part, bool> insert(Key const & key, util::ndarray<T, 2> const & data)
 		{
 			auto it = data_.find(key);
 			if (it != data_.end())
@@ -126,7 +126,7 @@ namespace psemek::util
 		std::size_t const padding_;
 		padding_mode const padding_mode_;
 		std::map<Key, atlas_part, Compare> data_;
-		util::array<T, 2> array_;
+		util::ndarray<T, 2> array_;
 		std::size_t free_start_ = 0;
 	};
 
diff --git a/libs/util/include/psemek/util/array.hpp b/libs/util/include/psemek/util/ndarray.hpp
similarity index 81%
rename from libs/util/include/psemek/util/array.hpp
rename to libs/util/include/psemek/util/ndarray.hpp
index 2a804506..6e629d48 100644
--- a/libs/util/include/psemek/util/array.hpp
+++ b/libs/util/include/psemek/util/ndarray.hpp
@@ -43,7 +43,7 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	struct array
+	struct ndarray
 	{
 		static_assert(N >= 1);
 
@@ -52,17 +52,17 @@ namespace psemek::util
 
 		using dims_type = std::array<std::size_t, N>;
 
-		array();
-		array(dims_type const & dims);
-		array(dims_type const & dims, T const & value);
-		array(dims_type const & dims, std::unique_ptr<T[]> data);
-		array(array &&);
+		ndarray();
+		ndarray(dims_type const & dims);
+		ndarray(dims_type const & dims, T const & value);
+		ndarray(dims_type const & dims, std::unique_ptr<T[]> data);
+		ndarray(ndarray &&);
 
-		array(array const &) = delete;
+		ndarray(ndarray const &) = delete;
 
-		array & operator = (array &&);
+		ndarray & operator = (ndarray &&);
 
-		array & operator = (array const &) = delete;
+		ndarray & operator = (ndarray const &) = delete;
 
 		dims_type const & dims() const
 		{
@@ -117,8 +117,8 @@ namespace psemek::util
 		template <typename I>
 		T const & operator()(std::initializer_list<I> const & index) const;
 
-		array copy() const;
-		array subarray(dims_type const & start, dims_type const & end) const;
+		ndarray copy() const;
+		ndarray subarray(dims_type const & start, dims_type const & end) const;
 
 		void resize(dims_type const & dims);
 
@@ -255,52 +255,52 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	auto array<T, N>::index_begin() const
+	auto ndarray<T, N>::index_begin() const
 	{
 		return indices().begin();
 	}
 
 	template <typename T, std::size_t N>
-	auto array<T, N>::index_end() const
+	auto ndarray<T, N>::index_end() const
 	{
 		return indices().end();
 	}
 
 	template <typename T, std::size_t N>
-	auto array<T, N>::indices() const
+	auto ndarray<T, N>::indices() const
 	{
 		return detail::array_index_range<N>{dims_};
 	}
 
 	template <typename T, std::size_t N>
-	array<T, N>::array()
+	ndarray<T, N>::ndarray()
 	{
 		dims_.fill(0);
 	}
 
 	template <typename T, std::size_t N>
-	array<T, N>::array(dims_type const & dims)
+	ndarray<T, N>::ndarray(dims_type const & dims)
 		: dims_{dims}
 	{
 		data_.reset(new T[size()]);
 	}
 
 	template <typename T, std::size_t N>
-	array<T, N>::array(dims_type const & dims, T const & value)
-		: array(dims)
+	ndarray<T, N>::ndarray(dims_type const & dims, T const & value)
+		: ndarray(dims)
 	{
 		fill(value);
 	}
 
 	template <typename T, std::size_t N>
-	array<T, N>::array(dims_type const & dims, std::unique_ptr<T[]> data)
+	ndarray<T, N>::ndarray(dims_type const & dims, std::unique_ptr<T[]> data)
 		: dims_{dims}
 	{
 		data_ = std::move(data);
 	}
 
 	template <typename T, std::size_t N>
-	array<T, N>::array(array && other)
+	ndarray<T, N>::ndarray(ndarray && other)
 		: data_{std::move(other.data_)}
 		, dims_{other.dims_}
 	{
@@ -308,7 +308,7 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	array<T, N> & array<T, N>::operator = (array && other)
+	ndarray<T, N> & ndarray<T, N>::operator = (ndarray && other)
 	{
 		if (this == &other)
 			return *this;
@@ -321,7 +321,7 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	std::size_t array<T, N>::size() const
+	std::size_t ndarray<T, N>::size() const
 	{
 		return detail::product(dims_);
 	}
@@ -329,7 +329,7 @@ namespace psemek::util
 	template <typename T, std::size_t N>
 	template <typename Index>
 	std::enable_if_t<detail::is_array_v<Index>, T &>
-		array<T, N>::operator()(Index const & index)
+		ndarray<T, N>::operator()(Index const & index)
 	{
 		return data_[detail::index(index, dims_)];
 	}
@@ -337,7 +337,7 @@ namespace psemek::util
 	template <typename T, std::size_t N>
 	template <typename Index>
 	std::enable_if_t<detail::is_array_v<Index>, T const &>
-		array<T, N>::operator()(Index const & index) const
+		ndarray<T, N>::operator()(Index const & index) const
 	{
 		return data_[detail::index(index, dims_)];
 	}
@@ -345,7 +345,7 @@ namespace psemek::util
 	template <typename T, std::size_t N>
 	template <typename ... Ixs>
 	std::enable_if_t<detail::is_index_v<Ixs...>, T &>
-		array<T, N>::operator()(Ixs ... ixs)
+		ndarray<T, N>::operator()(Ixs ... ixs)
 	{
 		static_assert(sizeof...(Ixs) == N);
 		dims_type dims{static_cast<std::size_t>(ixs)...};
@@ -355,7 +355,7 @@ namespace psemek::util
 	template <typename T, std::size_t N>
 	template <typename ... Ixs>
 	std::enable_if_t<detail::is_index_v<Ixs...>, T const &>
-		array<T, N>::operator()(Ixs ... ixs) const
+		ndarray<T, N>::operator()(Ixs ... ixs) const
 	{
 		static_assert(sizeof...(Ixs) == N);
 		dims_type dims{static_cast<std::size_t>(ixs)...};
@@ -364,34 +364,34 @@ namespace psemek::util
 
 	template <typename T, std::size_t N>
 	template <typename I>
-	T & array<T, N>::operator()(std::initializer_list<I> const & index)
+	T & ndarray<T, N>::operator()(std::initializer_list<I> const & index)
 	{
 		return data_[detail::index(std::data(index), dims_)];
 	}
 
 	template <typename T, std::size_t N>
 	template <typename I>
-	T const & array<T, N>::operator()(std::initializer_list<I> const & index) const
+	T const & ndarray<T, N>::operator()(std::initializer_list<I> const & index) const
 	{
 		return data_[detail::index(std::data(index), dims_)];
 	}
 
 	template <typename T, std::size_t N>
-	array<T, N> array<T, N>::copy() const
+	ndarray<T, N> ndarray<T, N>::copy() const
 	{
 		std::unique_ptr<T[]> data(new T[size()]);
 		std::copy(begin(), end(), data.get());
-		return array{dims_, std::move(data)};
+		return ndarray{dims_, std::move(data)};
 	}
 
 	template <typename T, std::size_t N>
-	array<T, N> array<T, N>::subarray(dims_type const & start, dims_type const & end) const
+	ndarray<T, N> ndarray<T, N>::subarray(dims_type const & start, dims_type const & end) const
 	{
 		dims_type size;
 		for (std::size_t i = 0; i < N; ++i)
 			size[i] = end[i] - start[i];
 
-		array<T, N> result(size);
+		ndarray<T, N> result(size);
 		for (auto const & idx : result.indices())
 		{
 			auto jdx = idx;
@@ -404,7 +404,7 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	void array<T, N>::resize(dims_type const & dims)
+	void ndarray<T, N>::resize(dims_type const & dims)
 	{
 		if (dims == dims_)
 			return;
@@ -414,7 +414,7 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	void array<T, N>::resize(dims_type const & dims, T const & value)
+	void ndarray<T, N>::resize(dims_type const & dims, T const & value)
 	{
 		if (dims == dims_)
 			return;
@@ -426,7 +426,7 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	void array<T, N>::assign(dims_type const & dims, T const & value)
+	void ndarray<T, N>::assign(dims_type const & dims, T const & value)
 	{
 		auto const size = detail::product(dims);
 		std::unique_ptr<T[]> data(new T[size]);
@@ -436,33 +436,33 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	std::unique_ptr<T[]> array<T, N>::release()
+	std::unique_ptr<T[]> ndarray<T, N>::release()
 	{
 		dims_.fill(0);
 		return std::move(data_);
 	}
 
 	template <typename T, std::size_t N>
-	bool array<T, N>::empty() const
+	bool ndarray<T, N>::empty() const
 	{
 		return detail::empty(dims_);
 	}
 
 	template <typename T, std::size_t N>
-	void array<T, N>::clear()
+	void ndarray<T, N>::clear()
 	{
 		data_.reset();
 		dims_.fill(0);
 	}
 
 	template <typename T, std::size_t N>
-	void array<T, N>::fill(T const & value)
+	void ndarray<T, N>::fill(T const & value)
 	{
 		std::fill(data_.get(), data_.get() + size(), value);
 	}
 
 	template <typename T, std::size_t N>
-	void array<T, N>::resize_impl(std::unique_ptr<T[]> data, dims_type const & dims)
+	void ndarray<T, N>::resize_impl(std::unique_ptr<T[]> data, dims_type const & dims)
 	{
 		std::array<std::size_t, N> min_dim;
 		for (std::size_t d = 0; d < N; ++d)
@@ -483,7 +483,7 @@ namespace psemek::util
 	}
 
 	template <typename T, std::size_t N>
-	void mirror(array<T, N> & a, std::size_t i)
+	void mirror(ndarray<T, N> & a, std::size_t i)
 	{
 		assert(i < N);
 
@@ -498,11 +498,11 @@ namespace psemek::util
 	}
 
 	template <typename F, typename T, std::size_t N>
-	auto map(F && f, array<T, N> const & a)
+	auto map(F && f, ndarray<T, N> const & a)
 	{
 		using R = std::decay_t<decltype(f(T{}))>;
 
-		array<R, N> r(a.dims());
+		ndarray<R, N> r(a.dims());
 
 		auto begin = a.begin();
 		auto end = a.end();
diff --git a/libs/util/include/psemek/util/spatial_array.hpp b/libs/util/include/psemek/util/spatial_array.hpp
index 671aa19a..d97ac953 100644
--- a/libs/util/include/psemek/util/spatial_array.hpp
+++ b/libs/util/include/psemek/util/spatial_array.hpp
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 #include <psemek/util/exception.hpp>
 #include <psemek/util/to_string.hpp>
 
@@ -49,7 +49,7 @@ namespace psemek::util
 
 	private:
 		Index origin_[N];
-		array<T, N> array_;
+		ndarray<T, N> array_;
 	};
 
 	template <typename T, std::size_t N, typename Index>
@@ -176,7 +176,7 @@ namespace psemek::util
 			}
 		}
 
-		array<T, N> new_array(new_size);
+		ndarray<T, N> new_array(new_size);
 		for (auto & v : new_array)
 			v = T();
 
diff --git a/libs/util/tests/array.cpp b/libs/util/tests/ndarray.cpp
similarity index 72%
rename from libs/util/tests/array.cpp
rename to libs/util/tests/ndarray.cpp
index 2f93addb..6fd2313b 100644
--- a/libs/util/tests/array.cpp
+++ b/libs/util/tests/ndarray.cpp
@@ -1,38 +1,38 @@
 #include <psemek/test/test.hpp>
 
-#include <psemek/util/array.hpp>
+#include <psemek/util/ndarray.hpp>
 
 using namespace psemek::util;
 
-test_case(util_array_empty)
+test_case(util_ndarray_empty)
 {
-	array<int, 1> a;
+	ndarray<int, 1> a;
 	expect_equal(a.size(), 0);
 	expect(a.empty());
 
-	array<int, 2> b;
+	ndarray<int, 2> b;
 	expect_equal(b.size(), 0);
 	expect(b.empty());
 
-	array<int, 3> c;
+	ndarray<int, 3> c;
 	expect_equal(c.size(), 0);
 	expect(c.empty());
 }
 
-test_case(util_array_size)
+test_case(util_ndarray_size)
 {
-	array<int, 1> a({16});
+	ndarray<int, 1> a({16});
 	expect(!a.empty());
 	expect_equal(a.size(), 16);
 	expect_equal(a.width(), 16);
 
-	array<int, 2> b({8, 16});
+	ndarray<int, 2> b({8, 16});
 	expect(!b.empty());
 	expect_equal(b.size(), 8 * 16);
 	expect_equal(b.width(), 8);
 	expect_equal(b.height(), 16);
 
-	array<int, 3> c({4, 8, 16});
+	ndarray<int, 3> c({4, 8, 16});
 	expect(!c.empty());
 	expect_equal(c.size(), 4 * 8 * 16);
 	expect_equal(c.width(), 4);
@@ -40,20 +40,20 @@ test_case(util_array_size)
 	expect_equal(c.depth(), 16);
 }
 
-test_case(util_array_init)
+test_case(util_ndarray_init)
 {
-	array<int, 2> a({16, 16}, 42);
+	ndarray<int, 2> a({16, 16}, 42);
 	for (std::size_t i = 0; i < a.width(); ++i)
 		for (std::size_t j = 0; j < a.height(); ++j)
 			expect_equal(a(i, j), 42);
 }
 
-test_case(util_array_release)
+test_case(util_ndarray_release)
 {
 	std::unique_ptr<int[]> p(new int[256]);
 	auto q = p.get();
 
-	array<int, 2> a({16, 16}, std::move(p));
+	ndarray<int, 2> a({16, 16}, std::move(p));
 	expect_equal(a.size(), 256);
 	expect_equal(a.width(), 16);
 	expect_equal(a.height(), 16);
@@ -66,20 +66,20 @@ test_case(util_array_release)
 	expect_equal_ptr(p.get(), q);
 }
 
-test_case(util_array_move)
+test_case(util_ndarray_move)
 {
-	array<int, 2> a({16, 16}, 42);
+	ndarray<int, 2> a({16, 16}, 42);
 
 	auto p = a.data();
 
-	array<int, 2> b(std::move(a));
+	ndarray<int, 2> b(std::move(a));
 	expect(a.empty());
 	expect_equal(a.size(), 0);
 	expect(!b.empty());
 	expect_equal(b.size(), 256);
 	expect_equal_ptr(b.data(), p);
 
-	array<int, 2> c;
+	ndarray<int, 2> c;
 	c = std::move(b);
 	expect(b.empty());
 	expect_equal(b.size(), 0);
@@ -88,11 +88,11 @@ test_case(util_array_move)
 	expect_equal_ptr(c.data(), p);
 }
 
-test_case(util_array_copy)
+test_case(util_ndarray_copy)
 {
-	array<int, 2> a({16, 16}, 42);
+	ndarray<int, 2> a({16, 16}, 42);
 
-	array<int, 2> b(a.copy());
+	ndarray<int, 2> b(a.copy());
 	expect(!a.empty());
 	expect(!b.empty());
 	expect_equal(a.size(), b.size());
@@ -100,12 +100,12 @@ test_case(util_array_copy)
 	expect_different_ptr(a.data(), b.data());
 }
 
-test_case(util_array_lifetime)
+test_case(util_ndarray_lifetime)
 {
 	std::weak_ptr<int> p;
 
 	{
-		array<std::shared_ptr<int>, 2> a({16, 16}, std::make_shared<int>(42));
+		ndarray<std::shared_ptr<int>, 2> a({16, 16}, std::make_shared<int>(42));
 		p = a(0, 0);
 		for (std::size_t i = 0; i < a.width(); ++i)
 		{
diff --git a/tools/noise-generator/generator.cpp b/tools/noise-generator/generator.cpp
index 9845a371..667af1c4 100644
--- a/tools/noise-generator/generator.cpp
+++ b/tools/noise-generator/generator.cpp
@@ -33,7 +33,7 @@ auto moore_neighbourhood()
 }
 
 template <int D>
-util::array<std::uint8_t, D> generate_perlin(std::array<std::size_t, D> const & size, int octaves, int minoctave, float power,
+util::ndarray<std::uint8_t, D> generate_perlin(std::array<std::size_t, D> const & size, int octaves, int minoctave, float power,
 	bool invert, bool tile, float gamma, bool remap, std::uint64_t seed)
 {
 	std::cout << "Generating " << D << "D perlin noise\n";
@@ -68,7 +68,7 @@ util::array<std::uint8_t, D> generate_perlin(std::array<std::size_t, D> const &
 		for (int d = 0; d < D; ++d)
 			octave_size[d] = (1 << o) + (tile ? 0 : 1);
 
-		util::array<math::vector<float, D>, D> gradients(octave_size);
+		util::ndarray<math::vector<float, D>, D> gradients(octave_size);
 		for (auto & g : gradients)
 			g = random_gradient(rng);
 
@@ -91,7 +91,7 @@ util::array<std::uint8_t, D> generate_perlin(std::array<std::size_t, D> const &
 	for (int d = 0; d < D; ++d)
 		largest_size = std::max<int>(largest_size, size[d]);
 
-	util::array<float, D> result_float(size);
+	util::ndarray<float, D> result_float(size);
 	math::interval<float> value_range;
 
 	for (auto idx : result_float.indices())
@@ -110,7 +110,7 @@ util::array<std::uint8_t, D> generate_perlin(std::array<std::size_t, D> const &
 		value_range |= v;
 	}
 
-	util::array<std::uint8_t, D> result(size);
+	util::ndarray<std::uint8_t, D> result(size);
 
 	for (auto idx : result.indices())
 	{
@@ -127,7 +127,7 @@ util::array<std::uint8_t, D> generate_perlin(std::array<std::size_t, D> const &
 }
 
 template <int D>
-util::array<std::uint8_t, D> generate_worley(std::array<std::size_t, D> const & size, int octaves, int minoctave, float power,
+util::ndarray<std::uint8_t, D> generate_worley(std::array<std::size_t, D> const & size, int octaves, int minoctave, float power,
 	bool invert, bool tile, float gamma, bool remap, std::uint64_t seed)
 {
 	std::cout << "Generating " << D << "D worley noise\n";
@@ -150,7 +150,7 @@ util::array<std::uint8_t, D> generate_worley(std::array<std::size_t, D> const &
 
 	random::generator rng{seed, 0};
 
-	std::vector<util::array<math::point<float, D>, D>> octave_points;
+	std::vector<util::ndarray<math::point<float, D>, D>> octave_points;
 	std::vector<float> octave_weights;
 	float sum_octave_weights = 0.f;
 
@@ -160,7 +160,7 @@ util::array<std::uint8_t, D> generate_worley(std::array<std::size_t, D> const &
 		for (int d = 0; d < D; ++d)
 			octave_size[d] = (1 << o);
 
-		util::array<math::point<float, D>, D> octave(octave_size);
+		util::ndarray<math::point<float, D>, D> octave(octave_size);
 
 		for (auto idx : octave.indices())
 		{
@@ -179,7 +179,7 @@ util::array<std::uint8_t, D> generate_worley(std::array<std::size_t, D> const &
 	for (auto & w : octave_weights)
 		w /= sum_octave_weights;
 
-	util::array<float, D> result_float(size);
+	util::ndarray<float, D> result_float(size);
 	math::interval<float> value_range;
 
 	for (auto idx : result_float.indices())
@@ -249,7 +249,7 @@ util::array<std::uint8_t, D> generate_worley(std::array<std::size_t, D> const &
 		value_range |= v;
 	}
 
-	util::array<std::uint8_t, D> result(size);
+	util::ndarray<std::uint8_t, D> result(size);
 
 	for (auto idx : result.indices())
 	{
@@ -507,7 +507,7 @@ int main(int argc, char * argv[])
 
 	if (*dim == 2)
 	{
-		util::array<std::uint8_t, 2> result;
+		util::ndarray<std::uint8_t, 2> result;
 		if (type == "perlin")
 			result = generate_perlin<2>({*sizex, *sizey}, *octaves, *minoctave, *power, *invert, *tile, *gamma, *remap, *seed);
 		else if (type == "worley")
@@ -539,7 +539,7 @@ int main(int argc, char * argv[])
 
 	if (*dim == 3)
 	{
-		util::array<std::uint8_t, 3> result;
+		util::ndarray<std::uint8_t, 3> result;
 		if (type == "perlin")
 			generate_perlin<3>({*sizex, *sizey, *sizez}, *octaves, *minoctave, *power, *invert, *tile, *gamma, *remap, *seed);
 		else if (type == "worley")