diff --git a/include/aare/Cluster.hpp b/include/aare/Cluster.hpp
index 82087df..c4dff2f 100644
--- a/include/aare/Cluster.hpp
+++ b/include/aare/Cluster.hpp
@@ -79,7 +79,7 @@ struct Cluster {
  * @brief Reduce a cluster to a 2x2 cluster by selecting the 2x2 block with the
  * highest sum.
  * @param c Cluster to reduce
- * @return Cluster<T, 2, 2, uint16_t> reduced cluster
+ * @return reduced cluster
  */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = int16_t>
@@ -89,6 +89,7 @@ reduce_to_2x2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
     static_assert(ClusterSizeX >= 2 && ClusterSizeY >= 2,
                   "Cluster sizes must be at least 2x2 for reduction to 2x2");
 
+    // TODO maybe add sanity check and check that center is in max subcluster
     Cluster<T, 2, 2, CoordType> result;
 
     auto [sum, index] = c.max_sum_2x2();
@@ -115,112 +116,88 @@ reduce_to_2x2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
     return result;
 }
 
-template <typename T>
-Cluster<T, 3, 3, uint16_t>
-reduce_5x5_to_3x3(const Cluster<T, 5, 5, uint16_t> &c) {
-    Cluster<T, 3, 3, uint16_t> result;
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t>
+inline std::pair<T, uint16_t>
+max_3x3_sum(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cluster) {
 
-    // Reduce the 5x5 cluster to a 3x3 cluster by selecting the 3x3 block with
-    // the highest sum
-    std::array<T, 9> sum_3x3_subclusters;
+    if constexpr (ClusterSizeX == 3 && ClusterSizeY == 3) {
+        return std::make_pair(cluster.sum(), 0);
+    } else {
 
-    // Write out the sums in the hope that the compiler can optimize this
-    sum_3x3_subclusters[0] = c.data[0] + c.data[1] + c.data[2] + c.data[5] +
-                             c.data[6] + c.data[7] + c.data[10] + c.data[11] +
-                             c.data[12];
-    sum_3x3_subclusters[1] = c.data[1] + c.data[2] + c.data[3] + c.data[6] +
-                             c.data[7] + c.data[8] + c.data[11] + c.data[12] +
-                             c.data[13];
-    sum_3x3_subclusters[2] = c.data[2] + c.data[3] + c.data[4] + c.data[7] +
-                             c.data[8] + c.data[9] + c.data[12] + c.data[13] +
-                             c.data[14];
-    sum_3x3_subclusters[3] = c.data[5] + c.data[6] + c.data[7] + c.data[10] +
-                             c.data[11] + c.data[12] + c.data[15] + c.data[16] +
-                             c.data[17];
-    sum_3x3_subclusters[4] = c.data[6] + c.data[7] + c.data[8] + c.data[11] +
-                             c.data[12] + c.data[13] + c.data[16] + c.data[17] +
-                             c.data[18];
-    sum_3x3_subclusters[5] = c.data[7] + c.data[8] + c.data[9] + c.data[12] +
-                             c.data[13] + c.data[14] + c.data[17] + c.data[18] +
-                             c.data[19];
-    sum_3x3_subclusters[6] = c.data[10] + c.data[11] + c.data[12] + c.data[15] +
-                             c.data[16] + c.data[17] + c.data[20] + c.data[21] +
-                             c.data[22];
-    sum_3x3_subclusters[7] = c.data[11] + c.data[12] + c.data[13] + c.data[16] +
-                             c.data[17] + c.data[18] + c.data[21] + c.data[22] +
-                             c.data[23];
-    sum_3x3_subclusters[8] = c.data[12] + c.data[13] + c.data[14] + c.data[17] +
-                             c.data[18] + c.data[19] + c.data[22] + c.data[23] +
-                             c.data[24];
+        size_t index = 0;
+        T max_3x3_subcluster_sum = 0;
+        for (size_t i = 0; i < ClusterSizeY - 2; ++i) {
+            for (size_t j = 0; j < ClusterSizeX - 2; ++j) {
 
-    auto index = std::max_element(sum_3x3_subclusters.begin(),
-                                  sum_3x3_subclusters.end()) -
-                 sum_3x3_subclusters.begin();
+                T sum = cluster.data[i * ClusterSizeX + j] +
+                        cluster.data[i * ClusterSizeX + j + 1] +
+                        cluster.data[i * ClusterSizeX + j + 2] +
+                        cluster.data[(i + 1) * ClusterSizeX + j] +
+                        cluster.data[(i + 1) * ClusterSizeX + j + 1] +
+                        cluster.data[(i + 1) * ClusterSizeX + j + 2] +
+                        cluster.data[(i + 2) * ClusterSizeX + j] +
+                        cluster.data[(i + 2) * ClusterSizeX + j + 1] +
+                        cluster.data[(i + 2) * ClusterSizeX + j + 2];
+                if (sum > max_3x3_subcluster_sum) {
+                    max_3x3_subcluster_sum = sum;
+                    index = i * (ClusterSizeX - 2) + j;
+                }
+            }
+        }
 
-    switch (index) {
-    case 0:
-        result.x = c.x - 1;
-        result.y = c.y + 1;
-        result.data = {c.data[0], c.data[1],  c.data[2],  c.data[5], c.data[6],
-                       c.data[7], c.data[10], c.data[11], c.data[12]};
-        break;
-    case 1:
-        result.x = c.x;
-        result.y = c.y + 1;
-        result.data = {c.data[1], c.data[2],  c.data[3],  c.data[6], c.data[7],
-                       c.data[8], c.data[11], c.data[12], c.data[13]};
-        break;
-    case 2:
-        result.x = c.x + 1;
-        result.y = c.y + 1;
-        result.data = {c.data[2], c.data[3],  c.data[4],  c.data[7], c.data[8],
-                       c.data[9], c.data[12], c.data[13], c.data[14]};
-        break;
-    case 3:
-        result.x = c.x - 1;
-        result.y = c.y;
-        result.data = {c.data[5],  c.data[6],  c.data[7],
-                       c.data[10], c.data[11], c.data[12],
-                       c.data[15], c.data[16], c.data[17]};
-        break;
-    case 4:
-        result.x = c.x + 1;
-        result.y = c.y;
-        result.data = {c.data[6],  c.data[7],  c.data[8],
-                       c.data[11], c.data[12], c.data[13],
-                       c.data[16], c.data[17], c.data[18]};
-        break;
-    case 5:
-        result.x = c.x + 1;
-        result.y = c.y;
-        result.data = {c.data[7],  c.data[8],  c.data[9],
-                       c.data[12], c.data[13], c.data[14],
-                       c.data[17], c.data[18], c.data[19]};
-        break;
-    case 6:
-        result.x = c.x + 1;
-        result.y = c.y - 1;
-        result.data = {c.data[10], c.data[11], c.data[12],
-                       c.data[15], c.data[16], c.data[17],
-                       c.data[20], c.data[21], c.data[22]};
-        break;
-    case 7:
-        result.x = c.x + 1;
-        result.y = c.y - 1;
-        result.data = {c.data[11], c.data[12], c.data[13],
-                       c.data[16], c.data[17], c.data[18],
-                       c.data[21], c.data[22], c.data[23]};
-        break;
-    case 8:
-        result.x = c.x + 1;
-        result.y = c.y - 1;
-        result.data = {c.data[12], c.data[13], c.data[14],
-                       c.data[17], c.data[18], c.data[19],
-                       c.data[22], c.data[23], c.data[24]};
-        break;
+        return std::make_pair(max_3x3_subcluster_sum, index);
     }
+}
 
-    // do some stuff
+/**
+ * @brief Reduce a cluster to a 3x3 cluster by selecting the 3x3 block with the
+ * highest sum.
+ * @param c Cluster to reduce
+ * @return reduced cluster
+ */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t>
+Cluster<T, 3, 3, CoordType>
+reduce_to_3x3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
+
+    static_assert(ClusterSizeX >= 3 && ClusterSizeY >= 3,
+                  "Cluster sizes must be at least 3x3 for reduction to 3x3");
+
+    Cluster<T, 3, 3, CoordType> result;
+
+    // TODO maybe add sanity check and check that center is in max subcluster
+
+    auto [sum, index] = max_3x3_sum(c);
+
+    int16_t cluster_center_index =
+        (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
+
+    int16_t index_center_max_3x3_subcluster =
+        (int(index / (ClusterSizeX - 2))) * ClusterSizeX + ClusterSizeX +
+        index % (ClusterSizeX - 2) + 1;
+
+    int16_t index_3x3_subcluster_cluster_center =
+        int((cluster_center_index - 1 - ClusterSizeX) / ClusterSizeX) *
+            (ClusterSizeX - 2) +
+        (cluster_center_index - 1 - ClusterSizeX) % ClusterSizeX;
+
+    result.x =
+        c.x + (index % (ClusterSizeX - 2) -
+               (index_3x3_subcluster_cluster_center % (ClusterSizeX - 2)));
+    result.y =
+        c.y - (index / (ClusterSizeX - 2) -
+               (index_3x3_subcluster_cluster_center / (ClusterSizeX - 2)));
+
+    result.data = {c.data[index_center_max_3x3_subcluster - ClusterSizeX - 1],
+                   c.data[index_center_max_3x3_subcluster - ClusterSizeX],
+                   c.data[index_center_max_3x3_subcluster - ClusterSizeX + 1],
+                   c.data[index_center_max_3x3_subcluster - 1],
+                   c.data[index_center_max_3x3_subcluster],
+                   c.data[index_center_max_3x3_subcluster + 1],
+                   c.data[index_center_max_3x3_subcluster + ClusterSizeX - 1],
+                   c.data[index_center_max_3x3_subcluster + ClusterSizeX],
+                   c.data[index_center_max_3x3_subcluster + ClusterSizeX + 1]};
     return result;
 }
 
diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp
index 07e1840..48bbbbd 100644
--- a/include/aare/ClusterVector.hpp
+++ b/include/aare/ClusterVector.hpp
@@ -173,7 +173,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
 };
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
+          typename CoordType = uint16_t>
 ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
     const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
         &cv) {
@@ -184,12 +184,14 @@ ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
     return result;
 }
 
-template <typename T>
-ClusterVector<Cluster<T, 3, 3, uint16_t>>
-reduce_5x5_to_3x3(const ClusterVector<Cluster<T, 5, 5, uint16_t>> &cv) {
-    ClusterVector<Cluster<T, 3, 3, uint16_t>> result;
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+ClusterVector<Cluster<T, 3, 3, CoordType>> reduce_to_3x3(
+    const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
+        &cv) {
+    ClusterVector<Cluster<T, 3, 3, CoordType>> result;
     for (const auto &c : cv) {
-        result.push_back(reduce_5x5_to_3x3(c));
+        result.push_back(reduce_to_3x3(c));
     }
     return result;
 }
diff --git a/python/src/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp
index 5240ad4..5166464 100644
--- a/python/src/bind_ClusterVector.hpp
+++ b/python/src/bind_ClusterVector.hpp
@@ -106,7 +106,7 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
 
 template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = uint16_t>
-void define_reduction(py::module &m) {
+void define_2x2_reduction(py::module &m) {
     m.def("reduce_to_2x2",
           [](const ClusterVector<
               Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>> &cv) {
@@ -115,11 +115,15 @@ void define_reduction(py::module &m) {
           });
 }
 
+template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
 void define_3x3_reduction(py::module &m) {
-    m.def("reduce_5x5_to_3x3",
-          [](const ClusterVector<Cluster<int, 5, 5, uint16_t>> &cv) {
-              return new ClusterVector<Cluster<int, 3, 3, uint16_t>>(
-                  reduce_5x5_to_3x3(cv));
+
+    m.def("reduce_to_3x3",
+          [](const ClusterVector<
+              Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>> &cv) {
+              return new ClusterVector<Cluster<Type, 3, 3, CoordType>>(
+                  reduce_to_3x3(cv));
           });
 }
 
diff --git a/python/src/module.cpp b/python/src/module.cpp
index 43c7c28..feafbd6 100644
--- a/python/src/module.cpp
+++ b/python/src/module.cpp
@@ -48,7 +48,7 @@ double, 'f' for float)
     define_ClusterCollector<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);    \
     define_Cluster<T, N, M, U>(m, #N "x" #M #TYPE_CODE);                       \
     register_calculate_eta<T, N, M, U>(m);                                     \
-    define_reduction<T, N, M, U>(m);
+    define_2x2_reduction<T, N, M, U>(m);
 
 PYBIND11_MODULE(_aare, m) {
     define_file_io_bindings(m);
@@ -86,5 +86,16 @@ PYBIND11_MODULE(_aare, m) {
     DEFINE_CLUSTER_BINDINGS(double, 9, 9, uint16_t, d);
     DEFINE_CLUSTER_BINDINGS(float, 9, 9, uint16_t, f);
 
-    define_3x3_reduction(m);
+    define_3x3_reduction<int, 3, 3>(m);
+    define_3x3_reduction<double, 3, 3>(m);
+    define_3x3_reduction<float, 3, 3>(m);
+    define_3x3_reduction<int, 5, 5>(m);
+    define_3x3_reduction<double, 5, 5>(m);
+    define_3x3_reduction<float, 5, 5>(m);
+    define_3x3_reduction<int, 7, 7>(m);
+    define_3x3_reduction<double, 7, 7>(m);
+    define_3x3_reduction<float, 7, 7>(m);
+    define_3x3_reduction<int, 9, 9>(m);
+    define_3x3_reduction<double, 9, 9>(m);
+    define_3x3_reduction<float, 9, 9>(m);
 }
diff --git a/src/Cluster.test.cpp b/src/Cluster.test.cpp
index ab0b4f9..2d54363 100644
--- a/src/Cluster.test.cpp
+++ b/src/Cluster.test.cpp
@@ -23,6 +23,9 @@ TEST_CASE("Test sum of Cluster", "[.cluster]") {
 using ClusterTypes = std::variant<Cluster<int, 2, 2>, Cluster<int, 3, 3>,
                                   Cluster<int, 5, 5>, Cluster<int, 2, 3>>;
 
+using ClusterTypesLargerThan2x2 =
+    std::variant<Cluster<int, 3, 3>, Cluster<int, 4, 4>, Cluster<int, 5, 5>>;
+
 TEST_CASE("Test reduce to 2x2 Cluster", "[.cluster]") {
     auto [cluster, expected_reduced_cluster] = GENERATE(
         std::make_tuple(ClusterTypes{Cluster<int, 2, 2>{5, 5, {1, 2, 3, 4}}},
@@ -60,4 +63,41 @@ TEST_CASE("Test reduce to 2x2 Cluster", "[.cluster]") {
     CHECK(std::equal(reduced_cluster.data.begin(),
                      reduced_cluster.data.begin() + 4,
                      expected_reduced_cluster.data.begin()));
+}
+
+TEST_CASE("Test reduce to 3x3 Clsuter", "[.cluster]") {
+    auto [cluster, expected_reduced_cluster] = GENERATE(
+        std::make_tuple(ClusterTypesLargerThan2x2{Cluster<int, 3, 3>{
+                            5, 5, {1, 1, 1, 1, 3, 1, 1, 1, 1}}},
+                        Cluster<int, 3, 3>{5, 5, {1, 1, 1, 1, 3, 1, 1, 1, 1}}),
+        std::make_tuple(
+            ClusterTypesLargerThan2x2{Cluster<int, 4, 4>{
+                5, 5, {2, 2, 1, 1, 2, 2, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1}}},
+            Cluster<int, 3, 3>{4, 6, {2, 2, 1, 2, 2, 1, 1, 1, 3}}),
+        std::make_tuple(
+            ClusterTypesLargerThan2x2{Cluster<int, 4, 4>{
+                5, 5, {1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 3, 1, 1, 1, 1, 1}}},
+            Cluster<int, 3, 3>{5, 6, {1, 2, 2, 1, 2, 2, 1, 3, 1}}),
+        std::make_tuple(
+            ClusterTypesLargerThan2x2{Cluster<int, 4, 4>{
+                5, 5, {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 1, 2, 2}}},
+            Cluster<int, 3, 3>{5, 5, {1, 1, 1, 1, 3, 2, 1, 2, 2}}),
+        std::make_tuple(
+            ClusterTypesLargerThan2x2{Cluster<int, 4, 4>{
+                5, 5, {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 1, 2, 2, 1, 1}}},
+            Cluster<int, 3, 3>{4, 5, {1, 1, 1, 2, 2, 3, 2, 2, 1}}),
+        std::make_tuple(ClusterTypesLargerThan2x2{Cluster<int, 5, 5>{
+                            5, 5, {1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 3,
+                                   1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1}}},
+                        Cluster<int, 3, 3>{4, 5, {1, 2, 1, 2, 2, 3, 1, 2, 1}}));
+
+    auto reduced_cluster = std::visit(
+        [](const auto &clustertype) { return reduce_to_3x3(clustertype); },
+        cluster);
+
+    CHECK(reduced_cluster.x == expected_reduced_cluster.x);
+    CHECK(reduced_cluster.y == expected_reduced_cluster.y);
+    CHECK(std::equal(reduced_cluster.data.begin(),
+                     reduced_cluster.data.begin() + 9,
+                     expected_reduced_cluster.data.begin()));
 }
\ No newline at end of file