3x3 reduction for general cluszter sizes

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;
 }
 

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;
 }