diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt
index 699b4c6..f749466 100644
--- a/benchmarks/CMakeLists.txt
+++ b/benchmarks/CMakeLists.txt
@@ -15,7 +15,7 @@ FetchContent_MakeAvailable(benchmark)
 
 add_executable(benchmarks)
 
-target_sources(benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp)
+target_sources(benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp reduce_benchmark.cpp)
 
 # Link Google Benchmark and other necessary libraries
 target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core aare_compiler_flags)
diff --git a/benchmarks/reduce_benchmark.cpp b/benchmarks/reduce_benchmark.cpp
new file mode 100644
index 0000000..2213624
--- /dev/null
+++ b/benchmarks/reduce_benchmark.cpp
@@ -0,0 +1,168 @@
+#include "aare/Cluster.hpp"
+#include <benchmark/benchmark.h>
+
+using namespace aare;
+
+class ClustersForReduceFixture : public benchmark::Fixture {
+  public:
+    Cluster<int, 5, 5> cluster_5x5{};
+    Cluster<int, 3, 3> cluster_3x3{};
+
+  private:
+    using benchmark::Fixture::SetUp;
+
+    void SetUp([[maybe_unused]] const benchmark::State &state) override {
+        int temp_data[25] = {1, 1, 1, 1, 1, 1, 1, 2, 1, 1,
+                             1, 2, 3, 1, 2, 1, 1, 1, 1, 2};
+        std::copy(std::begin(temp_data), std::end(temp_data),
+                  std::begin(cluster_5x5.data));
+
+        cluster_5x5.x = 5;
+        cluster_5x5.y = 5;
+
+        int temp_data2[9] = {1, 1, 1, 2, 3, 1, 2, 2, 1};
+        std::copy(std::begin(temp_data2), std::end(temp_data2),
+                  std::begin(cluster_3x3.data));
+
+        cluster_3x3.x = 5;
+        cluster_3x3.y = 5;
+    }
+
+    // void TearDown(::benchmark::State& state) {
+    // }
+};
+
+template <typename T>
+Cluster<T, 3, 3, int16_t> reduce_to_3x3(const Cluster<T, 5, 5, int16_t> &c) {
+    Cluster<T, 3, 3, int16_t> result;
+
+    // Write out the sums in the hope that the compiler can optimize this
+    std::array<T, 9> sum_3x3_subclusters;
+
+    // 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];
+
+    auto index = std::max_element(sum_3x3_subclusters.begin(),
+                                  sum_3x3_subclusters.end()) -
+                 sum_3x3_subclusters.begin();
+
+    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 result;
+}
+
+BENCHMARK_F(ClustersForReduceFixture, Reduce2x2)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        benchmark::DoNotOptimize(reduce_to_2x2<int, 3, 3, int16_t>(
+            cluster_3x3)); // make sure compiler evaluates the expression
+    }
+}
+
+BENCHMARK_F(ClustersForReduceFixture, SpecificReduce2x2)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        benchmark::DoNotOptimize(reduce_to_2x2<int>(cluster_3x3));
+    }
+}
+
+BENCHMARK_F(ClustersForReduceFixture, Reduce3x3)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        benchmark::DoNotOptimize(
+            reduce_to_3x3<int, 5, 5, int16_t>(cluster_5x5));
+    }
+}
+
+BENCHMARK_F(ClustersForReduceFixture, SpecificReduce3x3)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        benchmark::DoNotOptimize(reduce_to_3x3<int>(cluster_5x5));
+    }
+}
\ No newline at end of file
diff --git a/docs/src/ClusterVector.rst b/docs/src/ClusterVector.rst
index 93d2c7b..92376a0 100644
--- a/docs/src/ClusterVector.rst
+++ b/docs/src/ClusterVector.rst
@@ -12,4 +12,11 @@ ClusterVector
    :members:
    :undoc-members:
    :private-members:
+   
+
+**Free Functions:** 
+
+.. doxygenfunction:: aare::reduce_to_3x3(const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>&)
+   
+.. doxygenfunction:: aare::reduce_to_2x2(const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>&)
    
\ No newline at end of file
diff --git a/docs/src/pyClusterVector.rst b/docs/src/pyClusterVector.rst
index ff115c9..0ca2b25 100644
--- a/docs/src/pyClusterVector.rst
+++ b/docs/src/pyClusterVector.rst
@@ -33,4 +33,17 @@ C++ functions that support the ClusterVector or to view it as a numpy array.
     :members:
     :undoc-members:
     :show-inheritance:
-    :inherited-members:
\ No newline at end of file
+    :inherited-members:
+
+
+**Free Functions:** 
+
+.. autofunction:: reduce_to_3x3
+   :noindex:
+
+   Reduce a single Cluster to 3x3 by taking the 3x3 subcluster with highest photon energy.
+
+.. autofunction:: reduce_to_2x2
+   :noindex:
+
+   Reduce a single Cluster to 2x2 by taking the 2x2 subcluster with highest photon energy.
diff --git a/include/aare/CalculateEta.hpp b/include/aare/CalculateEta.hpp
index db17dad..2b28971 100644
--- a/include/aare/CalculateEta.hpp
+++ b/include/aare/CalculateEta.hpp
@@ -28,7 +28,7 @@ enum class pixel : int {
 template <typename T> struct Eta2 {
     double x;
     double y;
-    int c;
+    int c{0};
     T sum;
 };
 
@@ -70,6 +70,8 @@ calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
     size_t index_bottom_left_max_2x2_subcluster =
         (int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1);
 
+    // calculate direction of gradient
+
     // check that cluster center is in max subcluster
     if (cluster_center_index != index_bottom_left_max_2x2_subcluster &&
         cluster_center_index != index_bottom_left_max_2x2_subcluster + 1 &&
@@ -128,12 +130,15 @@ Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) {
     Eta2<T> eta{};
 
     if ((cl.data[0] + cl.data[1]) != 0)
-        eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
+        eta.x = static_cast<double>(cl.data[1]) /
+                (cl.data[0] + cl.data[1]); // between (0,1) the closer to zero
+                                           // left value probably larger
     if ((cl.data[0] + cl.data[2]) != 0)
-        eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
+        eta.y = static_cast<double>(cl.data[2]) /
+                (cl.data[0] + cl.data[2]); // between (0,1) the closer to zero
+                                           // bottom value probably larger
     eta.sum = cl.sum();
-    eta.c = static_cast<int>(corner::cBottomLeft); // TODO! This is not correct,
-                                                   // but need to put something
+
     return eta;
 }
 
@@ -150,13 +155,11 @@ template <typename T> Eta2<T> calculate_eta3(const Cluster<T, 3, 3> &cl) {
 
     eta.sum = sum;
 
-    eta.c = corner::cBottomLeft;
-
     if ((cl.data[3] + cl.data[4] + cl.data[5]) != 0)
 
         eta.x = static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
 
-                (cl.data[3] + cl.data[4] + cl.data[5]);
+                (cl.data[3] + cl.data[4] + cl.data[5]); // (-1,1)
 
     if ((cl.data[1] + cl.data[4] + cl.data[7]) != 0)
 
diff --git a/include/aare/Cluster.hpp b/include/aare/Cluster.hpp
old mode 100644
new mode 100755
index 889593b..32ab359
--- a/include/aare/Cluster.hpp
+++ b/include/aare/Cluster.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include "logger.hpp"
 #include <algorithm>
 #include <array>
 #include <cstdint>
@@ -74,6 +75,163 @@ 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 reduced cluster
+ */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t>
+Cluster<T, 2, 2, CoordType>
+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();
+
+    int16_t cluster_center_index =
+        (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
+
+    int16_t index_bottom_left_max_2x2_subcluster =
+        (int(index / (ClusterSizeX - 1))) * ClusterSizeX +
+        index % (ClusterSizeX - 1);
+
+    result.x =
+        c.x + (index_bottom_left_max_2x2_subcluster - cluster_center_index) %
+                  ClusterSizeX;
+
+    result.y =
+        c.y - (index_bottom_left_max_2x2_subcluster - cluster_center_index) /
+                  ClusterSizeX;
+    result.data = {
+        c.data[index_bottom_left_max_2x2_subcluster],
+        c.data[index_bottom_left_max_2x2_subcluster + 1],
+        c.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX],
+        c.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX + 1]};
+    return result;
+}
+
+template <typename T>
+Cluster<T, 2, 2, int16_t> reduce_to_2x2(const Cluster<T, 3, 3, int16_t> &c) {
+    Cluster<T, 2, 2, int16_t> result;
+
+    auto [s, i] = c.max_sum_2x2();
+    switch (i) {
+    case 0:
+        result.x = c.x - 1;
+        result.y = c.y + 1;
+        result.data = {c.data[0], c.data[1], c.data[3], c.data[4]};
+        break;
+    case 1:
+        result.x = c.x;
+        result.y = c.y + 1;
+        result.data = {c.data[1], c.data[2], c.data[4], c.data[5]};
+        break;
+    case 2:
+        result.x = c.x - 1;
+        result.y = c.y;
+        result.data = {c.data[3], c.data[4], c.data[6], c.data[7]};
+        break;
+    case 3:
+        result.x = c.x;
+        result.y = c.y;
+        result.data = {c.data[4], c.data[5], c.data[7], c.data[8]};
+        break;
+    }
+
+    return 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) {
+
+    if constexpr (ClusterSizeX == 3 && ClusterSizeY == 3) {
+        return std::make_pair(cluster.sum(), 0);
+    } else {
+
+        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) {
+
+                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;
+                }
+            }
+        }
+
+        return std::make_pair(max_3x3_subcluster_sum, index);
+    }
+}
+
+/**
+ * @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;
+}
+
 // Type Traits for is_cluster_type
 template <typename T>
 struct is_cluster : std::false_type {}; // Default case: Not a Cluster
diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp
index 0c81330..00f4d25 100644
--- a/include/aare/ClusterVector.hpp
+++ b/include/aare/ClusterVector.hpp
@@ -32,8 +32,7 @@ class ClusterVector; // Forward declaration
  */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType>
-class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> 
-{
+class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
 
     std::vector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> m_data{};
     int32_t m_frame_number{0}; // TODO! Check frame number size and type
@@ -173,4 +172,40 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
     }
 };
 
+/**
+ * @brief Reduce a cluster to a 2x2 cluster by selecting the 2x2 block with the
+ * highest sum.
+ * @param cv Clustervector containing clusters to reduce
+ * @return Clustervector with reduced clusters
+ */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
+    const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
+        &cv) {
+    ClusterVector<Cluster<T, 2, 2, CoordType>> result;
+    for (const auto &c : cv) {
+        result.push_back(reduce_to_2x2(c));
+    }
+    return result;
+}
+
+/**
+ * @brief Reduce a cluster to a 3x3 cluster by selecting the 3x3 block with the
+ * highest sum.
+ * @param cv Clustervector containing clusters to reduce
+ * @return Clustervector with reduced clusters
+ */
+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_to_3x3(c));
+    }
+    return result;
+}
+
 } // namespace aare
\ No newline at end of file
diff --git a/python/aare/__init__.py b/python/aare/__init__.py
index 1d34a85..c3b1c5a 100644
--- a/python/aare/__init__.py
+++ b/python/aare/__init__.py
@@ -17,7 +17,7 @@ from .ClusterVector import ClusterVector
 from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
 from ._aare import Interpolator
 from ._aare import calculate_eta2
-
+from ._aare import reduce_to_2x2, reduce_to_3x3
 
 from ._aare import apply_custom_weights
 
diff --git a/python/src/bind_Cluster.hpp b/python/src/bind_Cluster.hpp
index 690d0e8..d483ffd 100644
--- a/python/src/bind_Cluster.hpp
+++ b/python/src/bind_Cluster.hpp
@@ -24,7 +24,8 @@ void define_Cluster(py::module &m, const std::string &typestr) {
     py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
         m, class_name.c_str(), py::buffer_protocol())
 
-        .def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
+        .def(py::init([](uint8_t x, uint8_t y,
+                         py::array_t<Type, py::array::forcecast> data) {
             py::buffer_info buf_info = data.request();
             Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
             cluster.x = x;
@@ -34,31 +35,58 @@ void define_Cluster(py::module &m, const std::string &typestr) {
                 cluster.data[i] = r(i);
             }
             return cluster;
-        }));
+        }))
 
-    /*
-    //TODO! Review if to keep or not
-    .def_property(
-        "data",
-        [](ClusterType &c) -> py::array {
-            return py::array(py::buffer_info(
-                c.data, sizeof(Type),
-                py::format_descriptor<Type>::format(), // Type
-                                                       // format
-                1, // Number of dimensions
-                {static_cast<ssize_t>(ClusterSizeX *
-                                      ClusterSizeY)}, // Shape (flattened)
-                {sizeof(Type)} // Stride (step size between elements)
-                ));
+        // TODO! Review if to keep or not
+        .def_property_readonly(
+            "data",
+            [](Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> &c)
+                -> py::array {
+                return py::array(py::buffer_info(
+                    c.data.data(), sizeof(Type),
+                    py::format_descriptor<Type>::format(), // Type
+                                                           // format
+                    2, // Number of dimensions
+                    {static_cast<ssize_t>(ClusterSizeX),
+                     static_cast<ssize_t>(ClusterSizeY)}, // Shape (flattened)
+                    {sizeof(Type) * ClusterSizeY, sizeof(Type)}
+                    // Stride (step size between elements)
+                    ));
+            })
+
+        .def_readonly("x",
+                      &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::x)
+
+        .def_readonly("y",
+                      &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::y);
+}
+
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t>
+void reduce_to_3x3(py::module &m) {
+
+    m.def(
+        "reduce_to_3x3",
+        [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+            return reduce_to_3x3(cl);
         },
-        [](ClusterType &c, py::array_t<Type> arr) {
-            py::buffer_info buf_info = arr.request();
-            Type *ptr = static_cast<Type *>(buf_info.ptr);
-            std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY,
-                      c.data); // TODO dont iterate over centers!!!
+        py::return_value_policy::move,
+        "Reduce cluster to 3x3 subcluster by taking the 3x3 subcluster with "
+        "the highest photon energy.");
+}
 
-        });
-    */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t>
+void reduce_to_2x2(py::module &m) {
+
+    m.def(
+        "reduce_to_2x2",
+        [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+            return reduce_to_2x2(cl);
+        },
+        py::return_value_policy::move,
+        "Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with "
+        "the highest photon energy.");
 }
 
 #pragma GCC diagnostic pop
\ No newline at end of file
diff --git a/python/src/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp
index 9e9c4ab..597ae5f 100644
--- a/python/src/bind_ClusterVector.hpp
+++ b/python/src/bind_ClusterVector.hpp
@@ -104,4 +104,47 @@ 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_2x2_reduction(py::module &m) {
+    m.def(
+        "reduce_to_2x2",
+        [](const ClusterVector<
+            Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>> &cv) {
+            return new ClusterVector<Cluster<Type, 2, 2, CoordType>>(
+                reduce_to_2x2(cv));
+        },
+        R"(
+        
+        Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with 
+        the highest photon energy."
+        Parameters
+        ----------
+        cv : ClusterVector   
+        )",
+        py::arg("clustervector"));
+}
+
+template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void define_3x3_reduction(py::module &m) {
+
+    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));
+        },
+        R"(
+        
+        Reduce cluster to 3x3 subcluster by taking the 3x3 subcluster with 
+        the highest photon energy."
+        Parameters
+        ----------
+        cv : ClusterVector   
+        )",
+        py::arg("clustervector"));
+}
+
 #pragma GCC diagnostic pop
\ No newline at end of file
diff --git a/python/src/module.cpp b/python/src/module.cpp
index 9d75db6..cd802e8 100644
--- a/python/src/module.cpp
+++ b/python/src/module.cpp
@@ -47,7 +47,9 @@ double, 'f' for float)
     define_ClusterFileSink<T, N, M, U>(m, "Cluster" #N "x" #M #TYPE_CODE);     \
     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);
+    register_calculate_eta<T, N, M, U>(m);                                     \
+    define_2x2_reduction<T, N, M, U>(m);                                       \
+    reduce_to_2x2<T, N, M, U>(m);
 
 PYBIND11_MODULE(_aare, m) {
     define_file_io_bindings(m);
@@ -84,4 +86,30 @@ PYBIND11_MODULE(_aare, m) {
     DEFINE_CLUSTER_BINDINGS(int, 9, 9, uint16_t, i);
     DEFINE_CLUSTER_BINDINGS(double, 9, 9, uint16_t, d);
     DEFINE_CLUSTER_BINDINGS(float, 9, 9, uint16_t, f);
+
+    define_3x3_reduction<int, 3, 3, uint16_t>(m);
+    define_3x3_reduction<double, 3, 3, uint16_t>(m);
+    define_3x3_reduction<float, 3, 3, uint16_t>(m);
+    define_3x3_reduction<int, 5, 5, uint16_t>(m);
+    define_3x3_reduction<double, 5, 5, uint16_t>(m);
+    define_3x3_reduction<float, 5, 5, uint16_t>(m);
+    define_3x3_reduction<int, 7, 7, uint16_t>(m);
+    define_3x3_reduction<double, 7, 7, uint16_t>(m);
+    define_3x3_reduction<float, 7, 7, uint16_t>(m);
+    define_3x3_reduction<int, 9, 9, uint16_t>(m);
+    define_3x3_reduction<double, 9, 9, uint16_t>(m);
+    define_3x3_reduction<float, 9, 9, uint16_t>(m);
+
+    reduce_to_3x3<int, 3, 3, uint16_t>(m);
+    reduce_to_3x3<double, 3, 3, uint16_t>(m);
+    reduce_to_3x3<float, 3, 3, uint16_t>(m);
+    reduce_to_3x3<int, 5, 5, uint16_t>(m);
+    reduce_to_3x3<double, 5, 5, uint16_t>(m);
+    reduce_to_3x3<float, 5, 5, uint16_t>(m);
+    reduce_to_3x3<int, 7, 7, uint16_t>(m);
+    reduce_to_3x3<double, 7, 7, uint16_t>(m);
+    reduce_to_3x3<float, 7, 7, uint16_t>(m);
+    reduce_to_3x3<int, 9, 9, uint16_t>(m);
+    reduce_to_3x3<double, 9, 9, uint16_t>(m);
+    reduce_to_3x3<float, 9, 9, uint16_t>(m);
 }
diff --git a/python/tests/test_Cluster.py b/python/tests/test_Cluster.py
index ddaa6f3..506770b 100644
--- a/python/tests/test_Cluster.py
+++ b/python/tests/test_Cluster.py
@@ -101,6 +101,27 @@ def test_cluster_finder():
     assert clusters.size == 0
 
 
+def test_2x2_reduction(): 
+    """Test 2x2 Reduction"""
+    cluster = _aare.Cluster3x3i(5,5,np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32))
+
+    reduced_cluster = _aare.reduce_to_2x2(cluster) 
+
+    assert reduced_cluster.x == 4
+    assert reduced_cluster.y == 5
+    assert (reduced_cluster.data == np.array([[2, 3], [2, 2]], dtype=np.int32)).all()
+    
+    
+def test_3x3_reduction(): 
+    """Test 3x3 Reduction"""
+    cluster = _aare.Cluster5x5d(5,5,np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0, 2.0, 3.0,
+                                   1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=np.double))
+
+    reduced_cluster = _aare.reduce_to_3x3(cluster) 
+
+    assert reduced_cluster.x == 4
+    assert reduced_cluster.y == 5
+    assert (reduced_cluster.data == np.array([[1.0, 2.0, 1.0], [2.0, 2.0, 3.0], [1.0, 2.0, 1.0]], dtype=np.double)).all()
 
 
 
diff --git a/python/tests/test_ClusterVector.py b/python/tests/test_ClusterVector.py
index b64aeef..fa98e3f 100644
--- a/python/tests/test_ClusterVector.py
+++ b/python/tests/test_ClusterVector.py
@@ -5,7 +5,7 @@ import time
 from pathlib import Path
 import pickle
 
-from aare import ClusterFile
+from aare import ClusterFile, ClusterVector
 from aare import _aare
 from conftest import test_data_path
 
@@ -51,4 +51,36 @@ def test_make_a_hitmap_from_cluster_vector():
     # print(img)
     # print(ref)
     assert (img == ref).all()
-    
\ No newline at end of file
+
+
+def test_2x2_reduction(): 
+    cv = ClusterVector((3,3))
+
+    cv.push_back(_aare.Cluster3x3i(5, 5, np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32)))
+    cv.push_back(_aare.Cluster3x3i(5, 5, np.array([2, 2, 1, 2, 3, 1, 1, 1, 1], dtype=np.int32)))
+
+    reduced_cv = np.array(_aare.reduce_to_2x2(cv), copy=False) 
+
+    assert reduced_cv.size == 2
+    assert reduced_cv[0]["x"] == 4
+    assert reduced_cv[0]["y"] == 5
+    assert (reduced_cv[0]["data"] == np.array([[2, 3], [2, 2]], dtype=np.int32)).all()
+    assert reduced_cv[1]["x"] == 4
+    assert reduced_cv[1]["y"] == 6
+    assert (reduced_cv[1]["data"] == np.array([[2, 2], [2, 3]], dtype=np.int32)).all()
+    
+    
+def test_3x3_reduction(): 
+    cv = _aare.ClusterVector_Cluster5x5d()
+    
+    cv.push_back(_aare.Cluster5x5d(5,5,np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0, 2.0, 3.0,
+                                   1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=np.double)))
+    cv.push_back(_aare.Cluster5x5d(5,5,np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0, 2.0, 3.0,
+                                   1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=np.double)))
+    
+    reduced_cv = np.array(_aare.reduce_to_3x3(cv), copy=False)  
+
+    assert reduced_cv.size == 2
+    assert reduced_cv[0]["x"] == 4
+    assert reduced_cv[0]["y"] == 5
+    assert (reduced_cv[0]["data"] == np.array([[1.0, 2.0, 1.0], [2.0, 2.0, 3.0], [1.0, 2.0, 1.0]], dtype=np.double)).all()
\ No newline at end of file
diff --git a/src/Cluster.test.cpp b/src/Cluster.test.cpp
index ba9cda1..c1e424a 100644
--- a/src/Cluster.test.cpp
+++ b/src/Cluster.test.cpp
@@ -18,4 +18,86 @@ TEST_CASE("Test sum of Cluster", "[.cluster]") {
     Cluster<int, 2, 2> cluster{0, 0, {1, 2, 3, 4}};
 
     CHECK(cluster.sum() == 10);
+}
+
+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}}},
+                        Cluster<int, 2, 2>{4, 6, {1, 2, 3, 4}}),
+        std::make_tuple(
+            ClusterTypes{Cluster<int, 3, 3>{5, 5, {1, 1, 1, 1, 3, 2, 1, 2, 2}}},
+            Cluster<int, 2, 2>{5, 5, {3, 2, 2, 2}}),
+        std::make_tuple(
+            ClusterTypes{Cluster<int, 3, 3>{5, 5, {1, 1, 1, 2, 3, 1, 2, 2, 1}}},
+            Cluster<int, 2, 2>{4, 5, {2, 3, 2, 2}}),
+        std::make_tuple(
+            ClusterTypes{Cluster<int, 3, 3>{5, 5, {2, 2, 1, 2, 3, 1, 1, 1, 1}}},
+            Cluster<int, 2, 2>{4, 6, {2, 2, 2, 3}}),
+        std::make_tuple(
+            ClusterTypes{Cluster<int, 3, 3>{5, 5, {1, 2, 2, 1, 3, 2, 1, 1, 1}}},
+            Cluster<int, 2, 2>{5, 6, {2, 2, 3, 2}}),
+        std::make_tuple(ClusterTypes{Cluster<int, 5, 5>{
+                            5, 5, {1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 3,
+                                   2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}}},
+                        Cluster<int, 2, 2>{5, 6, {2, 2, 3, 2}}),
+        std::make_tuple(ClusterTypes{Cluster<int, 5, 5>{
+                            5, 5, {1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 3,
+                                   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}}},
+                        Cluster<int, 2, 2>{4, 6, {2, 2, 2, 3}}),
+        std::make_tuple(
+            ClusterTypes{Cluster<int, 2, 3>{5, 5, {2, 2, 3, 2, 1, 1}}},
+            Cluster<int, 2, 2>{4, 6, {2, 2, 3, 2}}));
+
+    auto reduced_cluster = std::visit(
+        [](const auto &clustertype) { return reduce_to_2x2(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() + 4,
+                     expected_reduced_cluster.data.begin()));
+}
+
+TEST_CASE("Test reduce to 3x3 Cluster", "[.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