diff --git a/.clangd b/.clangd
new file mode 100644
index 0000000..6418b4e
--- /dev/null
+++ b/.clangd
@@ -0,0 +1,2 @@
+CompileFlags:
+  Add: [-std=c++23]
diff --git a/CMakePresets.json b/CMakePresets.json
index 28135bc..f2607eb 100644
--- a/CMakePresets.json
+++ b/CMakePresets.json
@@ -12,7 +12,8 @@
             "displayName": "PartStacker default",
             "cacheVariables": {
                 "CMAKE_CXX_STANDARD": "23",
-                "CMAKE_CXX_STANDARD_REQUIRED": "ON"
+                "CMAKE_CXX_STANDARD_REQUIRED": "ON",
+                "CMAKE_EXPORT_COMPILE_COMMANDS": "ON"
             }
         },
         {
diff --git a/src/pstack/calc/CMakeLists.txt b/src/pstack/calc/CMakeLists.txt
index 05cd93d..34ca33f 100644
--- a/src/pstack/calc/CMakeLists.txt
+++ b/src/pstack/calc/CMakeLists.txt
@@ -1,4 +1,5 @@
 add_library(pstack_calc STATIC
+    binpack.cpp
     mesh.cpp
     rotations.cpp
     sinterbox.cpp
@@ -6,6 +7,8 @@ add_library(pstack_calc STATIC
     voxelize.cpp
 )
 target_sources(pstack_calc PUBLIC FILE_SET headers TYPE HEADERS FILES
+    binpack.hpp
+    binpack_types.hpp
     bool.hpp
     mesh.hpp
     part.hpp
diff --git a/src/pstack/calc/binpack.cpp b/src/pstack/calc/binpack.cpp
new file mode 100644
index 0000000..fa05f12
--- /dev/null
+++ b/src/pstack/calc/binpack.cpp
@@ -0,0 +1,105 @@
+#include "pstack/calc/binpack.hpp"
+#include "pstack/calc/binpack_cache.hpp"
+
+namespace pstack::calc::binpack {
+
+PackingResult combine(const std::array<PackingResult, 3> &results,
+                      const ItemPlacement &placement) {
+  PackingResult combined;
+  combined.total_score = results[0].total_score + results[1].total_score +
+                         results[2].total_score + placement.item.score;
+
+  size_t total_placements = results[0].placements.size() +
+                            results[1].placements.size() +
+                            results[2].placements.size() + 1;
+  combined.placements.reserve(total_placements);
+
+  combined.placements.push_back(placement); // Add the current item first
+
+  for (const auto &result : results) {
+    std::copy(result.placements.begin(), result.placements.end(),
+              std::back_inserter(combined.placements));
+  }
+
+  return combined;
+}
+
+PackingResult pack_items_impl(const std::vector<Item> &items,
+                              const PackingBox &into_box,
+                              BinPackCache<int, PackingResult> &cache,
+                              const geo::vector3<int> &min_volume) {
+  auto box_size = into_box.size;
+
+  if (!fits(sort(box_size), min_volume)) {
+    return {};
+  }
+
+  if (auto cached = cache.try_get(box_size); cached) {
+    return *cached;
+  }
+
+  PackingResult best;
+
+  for (auto &item : items) {
+    auto [sX, sY, sZ] = item.size;
+    if (sX > box_size.x || sY > box_size.y || sZ > box_size.z) {
+      continue;
+    }
+    // We test all 6 possible ways of splitting the box
+    std::array<std::array<PackingBox, 3>, 6> boxes;
+    auto offsetX = PackingBox::Offset<geo::Axis::X>{sX};
+    auto offsetY = PackingBox::Offset<geo::Axis::Y>{sY};
+    auto offsetZ = PackingBox::Offset<geo::Axis::Z>{sZ};
+
+    boxes[0] = into_box.split(offsetX, offsetY, offsetZ);
+    boxes[1] = into_box.split(offsetX, offsetZ, offsetY);
+    boxes[2] = into_box.split(offsetY, offsetX, offsetZ);
+    boxes[3] = into_box.split(offsetY, offsetZ, offsetX);
+    boxes[4] = into_box.split(offsetZ, offsetX, offsetY);
+    boxes[5] = into_box.split(offsetZ, offsetY, offsetX);
+
+    for (const auto &box : boxes) {
+      std::array<PackingResult, 3> results;
+      int score = item.score;
+      for (int j = 0; j < 3; j++) {
+        results[j] = pack_items_impl(items, box[j], cache, min_volume);
+        score += results[j].total_score;
+      }
+      if (score > best.total_score) {
+        for (int j = 0; j < 3; j++) {
+          results[j].translate(box[j].origin);
+        }
+        best = combine(results, {item, into_box.origin});
+      }
+    }
+  }
+
+  best.translate(-into_box.origin);
+  geo::vector3<int> aabb = {0, 0, 0};
+  for (auto &p : best.placements) {
+    aabb = geo::component_max(aabb, p.position + p.item.size);
+  }
+  cache.add(aabb, box_size, best);
+  return best;
+}
+
+PackingResult pack_items(const std::vector<Item> &items,
+                         const PackingBox &initial_box) {
+  if (items.empty()) {
+    return {};
+  }
+
+  BinPackCache<int, PackingResult> cache(1000);
+
+  geo::vector3<int> min_volume{std::numeric_limits<int>::max(),
+                               std::numeric_limits<int>::max(),
+                               std::numeric_limits<int>::max()};
+
+  for (const auto &item : items) {
+    min_volume = geo::component_min(min_volume, geo::sort(item.size));
+  }
+
+  return pack_items_impl(items, initial_box, cache, min_volume);
+}
+
+} // namespace pstack::calc::binpack
diff --git a/src/pstack/calc/binpack.hpp b/src/pstack/calc/binpack.hpp
new file mode 100644
index 0000000..b47ac4d
--- /dev/null
+++ b/src/pstack/calc/binpack.hpp
@@ -0,0 +1,12 @@
+#ifndef PSTACK_CALC_BINPACK_HPP
+#define PSTACK_CALC_BINPACK_HPP
+
+#include "binpack_types.hpp"
+
+namespace pstack::calc::binpack {
+
+PackingResult pack_items(const std::vector<Item>& items, const PackingBox& initial_box);
+
+} // namespace pstack::calc::binpack
+
+#endif // PSTACK_CALC_BINPACK_HPP
\ No newline at end of file
diff --git a/src/pstack/calc/binpack_cache.hpp b/src/pstack/calc/binpack_cache.hpp
new file mode 100644
index 0000000..7904ca4
--- /dev/null
+++ b/src/pstack/calc/binpack_cache.hpp
@@ -0,0 +1,70 @@
+#ifndef PSTACK_CALC_SPATIAL_HASH_HPP
+#define PSTACK_CALC_SPATIAL_HASH_HPP
+
+#include "pstack/geo/vector3.hpp"
+#include <unordered_map>
+
+namespace pstack::calc {
+
+template <typename T> struct SpatialKey {
+  T x, y, z;
+  auto operator<=>(const SpatialKey &other) const = default;
+};
+
+template <typename CoordT, typename ValueType> struct BinPackCache {
+  BinPackCache(CoordT cell_size) : cell_size(cell_size) {}
+
+  SpatialKey<CoordT> key_for_vector(const geo::vector3<CoordT> &vector) const {
+    return {vector.x / cell_size, vector.y / cell_size, vector.z / cell_size};
+  }
+
+  void add(const geo::vector3<CoordT> &items_aabb,
+           const geo::vector3<CoordT> &box_size, const ValueType &value) {
+    auto key1 = key_for_vector(items_aabb);
+    auto key2 = key_for_vector(box_size);
+    auto sp_value = CachedValue{value, items_aabb, box_size};
+    for (auto x = key1.x; x <= key2.x; ++x) {
+      for (auto y = key1.y; y <= key2.y; ++y) {
+        for (auto z = key1.z; z <= key2.z; ++z) {
+          cells.insert({{x, y, z}, sp_value});
+        }
+      }
+    }
+  }
+
+  std::optional<ValueType> try_get(const geo::vector3<CoordT> &box_size) const {
+    auto key = key_for_vector(box_size);
+    auto range = cells.equal_range(key);
+    for (auto it = range.first; it != range.second; ++it) {
+      const auto &cached = it->second;
+      if (cached.items_aabb <= box_size && box_size <= cached.box_size) {
+        return cached.value;
+      }
+    }
+    return std::nullopt;
+  }
+
+  struct CachedValue {
+    ValueType value;
+    geo::vector3<CoordT> items_aabb;
+    geo::vector3<CoordT> box_size;
+  };
+
+  std::unordered_multimap<SpatialKey<CoordT>, CachedValue> cells;
+  CoordT cell_size;
+};
+
+} // namespace pstack::calc
+
+namespace std {
+template <typename T> struct hash<pstack::calc::SpatialKey<T>> {
+  std::size_t operator()(const pstack::calc::SpatialKey<T> &k) const {
+    std::size_t h1 = std::hash<T>{}(k.x);
+    std::size_t h2 = std::hash<T>{}(k.y);
+    std::size_t h3 = std::hash<T>{}(k.z);
+    return h1 ^ (h2 << 1) ^ (h3 << 2);
+  }
+};
+} // namespace std
+
+#endif // PSTACK_CALC_SPATIAL_HASH_HPP
diff --git a/src/pstack/calc/binpack_types.hpp b/src/pstack/calc/binpack_types.hpp
new file mode 100644
index 0000000..2aadcd8
--- /dev/null
+++ b/src/pstack/calc/binpack_types.hpp
@@ -0,0 +1,70 @@
+#ifndef PSTACK_CALC_BINPACK_TYPES_HPP
+#define PSTACK_CALC_BINPACK_TYPES_HPP
+
+#include <array>
+#include <cstdint>
+#include <vector>
+
+#include "pstack/geo/vector3.hpp"
+
+namespace pstack::calc::binpack {
+
+struct Item {
+  uint32_t id{};
+  uint32_t score{};
+  geo::vector3<int> size{};
+};
+
+struct PackingBox {
+  geo::vector3<int> origin{};
+  geo::vector3<int> size{};
+
+  template <geo::Axis axis> struct Offset {
+    int value;
+  };
+
+  // Split the box into 3 smaller boxes along 3 given axes at 3 given offsets
+  template <geo::Axis axis1, geo::Axis axis2, geo::Axis axis3>
+  std::array<PackingBox, 3> split(Offset<axis1> offset1, Offset<axis2> offset2,
+                                  Offset<axis3> offset3) const {
+    auto [reduced_1, remaining_1] = this->split(offset1);
+    auto [reduced_2, remaining_2] = reduced_1.split(offset2);
+    auto [_, remaining_3] = reduced_2.split(offset3);
+
+    return {remaining_1, remaining_2, remaining_3};
+  }
+
+private:
+  // Split the box into 2 smaller boxes along a given axis at a given offset
+  template <geo::Axis axis>
+  std::array<PackingBox, 2> split(Offset<axis> offset) const {
+    PackingBox piece1{*this};
+    PackingBox piece2{*this};
+
+    piece1.size[axis] = offset;
+    piece2.size[axis] -= offset;
+    piece2.origin[axis] += offset;
+
+    return {piece1, piece2};
+  }
+};
+
+struct ItemPlacement {
+  Item item;
+  geo::vector3<int> position{};
+};
+
+struct PackingResult {
+  uint32_t total_score{};
+  std::vector<ItemPlacement> placements;
+
+  void translate(const geo::vector3<int> &offset) {
+    for (auto &placement : placements) {
+      placement.position += offset;
+    }
+  }
+};
+
+} // namespace pstack::calc::binpack
+
+#endif // PSTACK_CALC_BINPACK_TYPES_HPP
diff --git a/src/pstack/geo/vector3.hpp b/src/pstack/geo/vector3.hpp
index 3aedbed..a0fa28b 100644
--- a/src/pstack/geo/vector3.hpp
+++ b/src/pstack/geo/vector3.hpp
@@ -2,15 +2,32 @@
 #define PSTACK_GEO_VECTOR3_HPP
 
 #include "pstack/geo/functions.hpp"
+#include <algorithm>
 #include <type_traits>
 
 namespace pstack::geo {
 
+enum class Axis {
+    X = 0,
+    Y = 1,
+    Z = 2
+};
+
 template <class T>
 struct vector3 {
     T x;
     T y;
     T z;
+
+    constexpr T& operator[](const Axis axis) {
+        auto index = static_cast<std::underlying_type_t<Axis>>(axis);
+        return *(&x + index);
+    }
+
+    constexpr const T& operator[](const Axis axis) const {
+        auto index = static_cast<std::underlying_type_t<Axis>>(axis);
+        return *(&x + index);
+    }
 };
 
 template <class T>
@@ -22,6 +39,41 @@ inline constexpr vector3<T> unit_y = { 0, 1, 0 };
 template <class T>
 inline constexpr vector3<T> unit_z = { 0, 0, 1 };
 
+template <class T>
+constexpr auto operator<=>(const vector3<T>& lhs, const vector3<T>& rhs) {
+    if (auto c = lhs.x <=> rhs.x; c != 0) {
+        return c;
+    }
+    if (auto c = lhs.y <=> rhs.y; c != 0) {
+        return c;
+    }
+    return lhs.z <=> rhs.z;
+}
+
+template<class T>
+constexpr vector3<T> component_min(const vector3<T>& lhs, const vector3<T>& rhs) {
+    return { std::min(lhs.x, rhs.x), std::min(lhs.y, rhs.y), std::min(lhs.z, rhs.z) };
+}
+
+template<class T>
+constexpr vector3<T> component_max(const vector3<T>& lhs, const vector3<T>& rhs) {
+    return { std::max(lhs.x, rhs.x), std::max(lhs.y, rhs.y), std::max(lhs.z, rhs.z) };
+}
+
+template <class T>
+constexpr vector3<T> sort(const vector3<T>& v) {
+    vector3<T> result = v;
+    if (result.x > result.z) std::swap(result.x, result.z);
+    if (result.x > result.y) std::swap(result.x, result.y);
+    if (result.y > result.z) std::swap(result.y, result.z);
+    return result;
+}
+
+template <class T>
+constexpr bool fits(const vector3<T>& a, const vector3<T>& b) {
+    return a.x >= b.x && a.y >= b.y && a.z >= b.z;
+}
+
 template <class T>
 constexpr vector3<T> operator+(const vector3<T>& lhs, const vector3<T>& rhs) {
     return { lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z };