Dev/expose sum 2x2 to python (#238)
Some checks failed
Build on RHEL8 / build (push) Failing after 3m15s
Build on RHEL9 / build (push) Failing after 3m16s

Saverio requested that max_sum_2x2 exposes index information in  python 
- max_sum_2x2 returns a corner as index
- replaced eta corner with corner enum class
- max_sum_2x2 now returns index as well in python 
- added link to Documenation in README

Note: Some Tests fail in EtaCalculation due to previous PR about
updating Eta 2x2 will fix in other PR
This commit is contained in:
2025-10-27 20:04:16 +01:00
committed by GitHub
14 changed files with 113 additions and 47 deletions

View File

@@ -1,6 +1,10 @@
# aare # aare
Data analysis library for PSI hybrid detectors Data analysis library for PSI hybrid detectors
## Documentation
Detailed documentation including installation can be found in [Documentation](https://slsdetectorgroup.github.io/aare/)
## Build and install ## Build and install

View File

@@ -1,5 +1,13 @@
# Release notes # Release notes
### 2025.10.27
Features:
- max_sum_2x2 including index of subcluster with highest energy is now available from Python API
- eta stores corner as enum class cTopLeft, cTopRight, BottomLeft, cBottomRight indicating 2x2 subcluster with largest energy relative to cluster center
- max_sum_2x2 returns corner as index
### 2025.10.1 ### 2025.10.1
Bugfixes: Bugfixes:

View File

@@ -3,16 +3,10 @@
#include "aare/Cluster.hpp" #include "aare/Cluster.hpp"
#include "aare/ClusterVector.hpp" #include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp" #include "aare/NDArray.hpp"
#include "aare/defs.hpp"
namespace aare { namespace aare {
enum class corner : int {
cTopLeft = 0,
cTopRight = 1,
cBottomLeft = 2,
cBottomRight = 3
};
enum class pixel : int { enum class pixel : int {
pBottomLeft = 0, pBottomLeft = 0,
pBottom = 1, pBottom = 1,
@@ -28,7 +22,7 @@ enum class pixel : int {
template <typename T> struct Eta2 { template <typename T> struct Eta2 {
double x; double x;
double y; double y;
int c{0}; corner c{0};
T sum; T sum;
}; };
@@ -66,11 +60,11 @@ calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX; (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
auto max_sum = cl.max_sum_2x2(); auto max_sum = cl.max_sum_2x2();
eta.sum = max_sum.first; eta.sum = max_sum.sum;
int c = max_sum.second; corner c = max_sum.index;
// subcluster top right from center // subcluster top right from center
switch (static_cast<corner>(c)) { switch (c) {
case (corner::cTopLeft): case (corner::cTopLeft):
if ((cl.data[cluster_center_index - 1] + if ((cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]) != 0) cl.data[cluster_center_index]) != 0)

View File

@@ -8,6 +8,7 @@
#pragma once #pragma once
#include "defs.hpp"
#include <algorithm> #include <algorithm>
#include <array> #include <array>
#include <cstdint> #include <cstdint>
@@ -45,7 +46,7 @@ struct Cluster {
* @return photon energy of subcluster, 2x2 subcluster index relative to * @return photon energy of subcluster, 2x2 subcluster index relative to
* cluster center * cluster center
*/ */
std::pair<T, int> max_sum_2x2() const { Sum_index_pair<T, corner> max_sum_2x2() const {
if constexpr (cluster_size_x == 3 && cluster_size_y == 3) { if constexpr (cluster_size_x == 3 && cluster_size_y == 3) {
std::array<T, 4> sum_2x2_subclusters; std::array<T, 4> sum_2x2_subclusters;
@@ -56,9 +57,11 @@ struct Cluster {
int index = std::max_element(sum_2x2_subclusters.begin(), int index = std::max_element(sum_2x2_subclusters.begin(),
sum_2x2_subclusters.end()) - sum_2x2_subclusters.end()) -
sum_2x2_subclusters.begin(); sum_2x2_subclusters.begin();
return std::make_pair(sum_2x2_subclusters[index], index); return Sum_index_pair<T, corner>{sum_2x2_subclusters[index],
corner{index}};
} else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) { } else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) {
return std::make_pair(data[0] + data[1] + data[2] + data[3], 0); return Sum_index_pair<T, corner>{
data[0] + data[1] + data[2] + data[3], corner{0}};
} else { } else {
constexpr size_t cluster_center_index = constexpr size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX; (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
@@ -97,7 +100,8 @@ struct Cluster {
int index = std::max_element(sum_2x2_subcluster.begin(), int index = std::max_element(sum_2x2_subcluster.begin(),
sum_2x2_subcluster.end()) - sum_2x2_subcluster.end()) -
sum_2x2_subcluster.begin(); sum_2x2_subcluster.begin();
return std::make_pair(sum_2x2_subcluster[index], index); return Sum_index_pair<T, corner>{sum_2x2_subcluster[index],
corner{index}};
} }
} }
}; };
@@ -125,8 +129,8 @@ reduce_to_2x2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX; (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
int16_t index_bottom_left_max_2x2_subcluster = int16_t index_bottom_left_max_2x2_subcluster =
(int(index / (ClusterSizeX - 1))) * ClusterSizeX + (int(static_cast<int>(index) / (ClusterSizeX - 1))) * ClusterSizeX +
index % (ClusterSizeX - 1); static_cast<int>(index) % (ClusterSizeX - 1);
result.x = result.x =
c.x + (index_bottom_left_max_2x2_subcluster - cluster_center_index) % c.x + (index_bottom_left_max_2x2_subcluster - cluster_center_index) %
@@ -149,22 +153,22 @@ Cluster<T, 2, 2, int16_t> reduce_to_2x2(const Cluster<T, 3, 3, int16_t> &c) {
auto [s, i] = c.max_sum_2x2(); auto [s, i] = c.max_sum_2x2();
switch (i) { switch (i) {
case 0: case corner::cTopLeft:
result.x = c.x - 1; result.x = c.x - 1;
result.y = c.y + 1; result.y = c.y + 1;
result.data = {c.data[0], c.data[1], c.data[3], c.data[4]}; result.data = {c.data[0], c.data[1], c.data[3], c.data[4]};
break; break;
case 1: case corner::cTopRight:
result.x = c.x; result.x = c.x;
result.y = c.y + 1; result.y = c.y + 1;
result.data = {c.data[1], c.data[2], c.data[4], c.data[5]}; result.data = {c.data[1], c.data[2], c.data[4], c.data[5]};
break; break;
case 2: case corner::cBottomLeft:
result.x = c.x - 1; result.x = c.x - 1;
result.y = c.y; result.y = c.y;
result.data = {c.data[3], c.data[4], c.data[6], c.data[7]}; result.data = {c.data[3], c.data[4], c.data[6], c.data[7]};
break; break;
case 3: case corner::cBottomRight:
result.x = c.x; result.x = c.x;
result.y = c.y; result.y = c.y;
result.data = {c.data[4], c.data[5], c.data[7], c.data[8]}; result.data = {c.data[4], c.data[5], c.data[7], c.data[8]};

View File

@@ -5,6 +5,7 @@
#include "aare/ClusterFinderMT.hpp" #include "aare/ClusterFinderMT.hpp"
#include "aare/ClusterVector.hpp" #include "aare/ClusterVector.hpp"
#include "aare/ProducerConsumerQueue.hpp" #include "aare/ProducerConsumerQueue.hpp"
#include "aare/defs.hpp"
namespace aare { namespace aare {

View File

@@ -453,7 +453,7 @@ bool ClusterFile<ClusterType, Enable>::is_selected(ClusterType &cl) {
if (m_noise_map) { if (m_noise_map) {
auto sum_1x1 = cl.data[cluster_center_index]; // central pixel auto sum_1x1 = cl.data[cluster_center_index]; // central pixel
auto sum_2x2 = cl.max_sum_2x2().first; // highest sum of 2x2 subclusters auto sum_2x2 = cl.max_sum_2x2().sum; // highest sum of 2x2 subclusters
auto total_sum = cl.sum(); // sum of all pixels auto total_sum = cl.sum(); // sum of all pixels
auto noise = auto noise =

View File

@@ -86,15 +86,14 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
/** /**
* @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in * @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
* each cluster * each cluster
* @return std::vector<T> vector of sums for each cluster * @return vector of sums index pairs for each cluster
*/ */
std::vector<T> sum_2x2() { std::vector<Sum_index_pair<T, corner>> sum_2x2() {
std::vector<T> sums_2x2(m_data.size()); std::vector<Sum_index_pair<T, corner>> sums_2x2(m_data.size());
std::transform(m_data.begin(), m_data.end(), sums_2x2.begin(), std::transform(
[](const ClusterType &cluster) { m_data.begin(), m_data.end(), sums_2x2.begin(),
return cluster.max_sum_2x2().first; [](const ClusterType &cluster) { return cluster.max_sum_2x2(); });
});
return sums_2x2; return sums_2x2;
} }

View File

@@ -331,6 +331,19 @@ enum DACIndex {
SLOW_ADC_TEMP SLOW_ADC_TEMP
}; };
// helper pair class to easily expose in python
template <typename T1, typename T2> struct Sum_index_pair {
T1 sum;
T2 index;
};
enum class corner : int {
cTopLeft = 0,
cTopRight = 1,
cBottomLeft = 2,
cBottomRight = 3
};
enum class TimingMode { Auto, Trigger }; enum class TimingMode { Auto, Trigger };
enum class FrameDiscardPolicy { NoDiscard, Discard, DiscardPartial }; enum class FrameDiscardPolicy { NoDiscard, Discard, DiscardPartial };

View File

@@ -58,7 +58,17 @@ void define_Cluster(py::module &m, const std::string &typestr) {
&Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::x) &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::x)
.def_readonly("y", .def_readonly("y",
&Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::y); &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::y)
.def(
"max_sum_2x2",
[](Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> &self) {
auto max_sum = self.max_sum_2x2();
return py::make_tuple(max_sum.sum,
static_cast<int>(max_sum.index));
},
R"(calculates sum of 2x2 subcluster with highest energy and index relative to cluster center 0: top_left, 1: top_right, 2: bottom_left, 3: bottom_right
)");
} }
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,

View File

@@ -44,11 +44,16 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
auto *vec = new std::vector<Type>(self.sum()); auto *vec = new std::vector<Type>(self.sum());
return return_vector(vec); return return_vector(vec);
}) })
.def("sum_2x2", .def(
"sum_2x2",
[](ClusterVector<ClusterType> &self) { [](ClusterVector<ClusterType> &self) {
auto *vec = new std::vector<Type>(self.sum_2x2()); auto *vec = new std::vector<Sum_index_pair<Type, corner>>(
self.sum_2x2());
return return_vector(vec); return return_vector(vec);
}) },
R"(calculates sum of 2x2 subcluster with highest energy and index relative to cluster center 0: top_left, 1: top_right, 2: bottom_left, 3: bottom_right
)")
.def_property_readonly("size", &ClusterVector<ClusterType>::size) .def_property_readonly("size", &ClusterVector<ClusterType>::size)
.def("item_size", &ClusterVector<ClusterType>::item_size) .def("item_size", &ClusterVector<ClusterType>::item_size)
.def_property_readonly("fmt", .def_property_readonly("fmt",

View File

@@ -114,4 +114,11 @@ PYBIND11_MODULE(_aare, m) {
reduce_to_3x3<int, 9, 9, uint16_t>(m); reduce_to_3x3<int, 9, 9, uint16_t>(m);
reduce_to_3x3<double, 9, 9, uint16_t>(m); reduce_to_3x3<double, 9, 9, uint16_t>(m);
reduce_to_3x3<float, 9, 9, uint16_t>(m); reduce_to_3x3<float, 9, 9, uint16_t>(m);
using Sum_index_pair_d = Sum_index_pair<double, corner>;
PYBIND11_NUMPY_DTYPE(Sum_index_pair_d, sum, index);
using Sum_index_pair_f = Sum_index_pair<float, corner>;
PYBIND11_NUMPY_DTYPE(Sum_index_pair_f, sum, index);
using Sum_index_pair_i = Sum_index_pair<int, corner>;
PYBIND11_NUMPY_DTYPE(Sum_index_pair_i, sum, index);
} }

View File

@@ -86,6 +86,17 @@ def test_calculate_eta():
assert eta2[1,0] == 0.5 assert eta2[1,0] == 0.5
assert eta2[1,1] == 0.4 #2/5 assert eta2[1,1] == 0.4 #2/5
def test_max_sum():
"""Max 2x2 Sum"""
cluster = _aare.Cluster3x3i(5,5,np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32))
max_sum = cluster.max_sum_2x2()
assert max_sum[0] == 9
assert max_sum[1] == 2
def test_cluster_finder(): def test_cluster_finder():
"""Test ClusterFinder""" """Test ClusterFinder"""

View File

@@ -32,6 +32,18 @@ def test_push_back_on_cluster_vector():
assert arr[0]['y'] == 22 assert arr[0]['y'] == 22
def test_max_2x2_sum():
"""max_2x2_sum"""
cv = _aare.ClusterVector_Cluster3x3i()
cv.push_back(_aare.Cluster3x3i(19, 22, np.array([0,1,0,2,3,0,2,1,0], dtype=np.int32)))
cv.push_back(_aare.Cluster3x3i(19, 22, np.ones(9, dtype=np.int32)))
assert cv.size == 2
max_2x2 = cv.sum_2x2()
assert max_2x2.size == 2
assert max_2x2[0]["sum"] == 8
assert max_2x2[0]["index"] == 2
def test_make_a_hitmap_from_cluster_vector(): def test_make_a_hitmap_from_cluster_vector():
cv = _aare.ClusterVector_Cluster3x3i() cv = _aare.ClusterVector_Cluster3x3i()

View File

@@ -21,22 +21,20 @@ using ClusterTypes =
auto get_test_parameters() { auto get_test_parameters() {
return GENERATE( return GENERATE(
std::make_tuple(ClusterTypes{Cluster<int, 2, 2>{0, 0, {1, 2, 3, 1}}}, std::make_tuple(ClusterTypes{Cluster<int, 2, 2>{0, 0, {1, 2, 3, 1}}},
Eta2<int>{2. / 3, 3. / 4, Eta2<int>{2. / 3, 3. / 4, corner::cTopLeft, 7}),
static_cast<int>(corner::cBottomLeft), 7}),
std::make_tuple( std::make_tuple(
ClusterTypes{Cluster<int, 3, 3>{0, 0, {1, 2, 3, 4, 5, 6, 1, 2, 7}}}, ClusterTypes{Cluster<int, 3, 3>{0, 0, {1, 2, 3, 4, 5, 6, 1, 2, 7}}},
Eta2<int>{6. / 11, 2. / 7, static_cast<int>(corner::cTopRight), Eta2<int>{6. / 11, 2. / 7, corner::cBottomRight, 20}),
20}),
std::make_tuple(ClusterTypes{Cluster<int, 5, 5>{ std::make_tuple(ClusterTypes{Cluster<int, 5, 5>{
0, 0, {1, 6, 7, 6, 5, 4, 3, 2, 1, 2, 8, 9, 8, 0, 0, {1, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 9, 8,
1, 4, 5, 6, 7, 8, 4, 1, 1, 1, 1, 1}}}, 1, 4, 1, 6, 7, 8, 1, 1, 1, 1, 1, 1}}},
Eta2<int>{8. / 17, 7. / 15, 9, 30}), Eta2<int>{8. / 17, 7. / 15, corner::cBottomLeft, 30}),
std::make_tuple( std::make_tuple(
ClusterTypes{Cluster<int, 4, 2>{0, 0, {1, 4, 7, 2, 5, 6, 4, 3}}}, ClusterTypes{Cluster<int, 4, 2>{0, 0, {1, 4, 7, 2, 5, 6, 4, 3}}},
Eta2<int>{4. / 10, 4. / 11, 1, 21}), Eta2<int>{4. / 10, 4. / 11, corner::cTopLeft, 21}),
std::make_tuple( std::make_tuple(
ClusterTypes{Cluster<int, 2, 3>{0, 0, {1, 3, 2, 3, 4, 2}}}, ClusterTypes{Cluster<int, 2, 3>{0, 0, {1, 3, 2, 3, 4, 2}}},
Eta2<int>{3. / 5, 2. / 5, 1, 11})); Eta2<int>{3. / 5, 2. / 5, corner::cBottomLeft, 11}));
} }
TEST_CASE("compute_largest_2x2_subcluster", "[eta_calculation]") { TEST_CASE("compute_largest_2x2_subcluster", "[eta_calculation]") {
@@ -91,7 +89,7 @@ TEST_CASE("Calculate eta2 for a 3x3 int32 cluster with the largest 2x2 sum in "
// 30, 23, 5 // 30, 23, 5
auto eta = calculate_eta2(cl); auto eta = calculate_eta2(cl);
CHECK(eta.c == static_cast<int>(corner::cBottomLeft)); CHECK(eta.c == corner::cBottomLeft);
CHECK(eta.x == 50.0 / (20 + 50)); // 4/(3+4) CHECK(eta.x == 50.0 / (20 + 50)); // 4/(3+4)
CHECK(eta.y == 50.0 / (23 + 50)); // 4/(1+4) CHECK(eta.y == 50.0 / (23 + 50)); // 4/(1+4)
CHECK(eta.sum == 30 + 23 + 20 + 50); CHECK(eta.sum == 30 + 23 + 20 + 50);
@@ -120,7 +118,7 @@ TEST_CASE("Calculate eta2 for a 3x3 int32 cluster with the largest 2x2 sum in "
// 8, 12, 5 // 8, 12, 5
auto eta = calculate_eta2(cl); auto eta = calculate_eta2(cl);
CHECK(eta.c == static_cast<int>(corner::cTopLeft)); CHECK(eta.c == corner::cTopLeft);
CHECK(eta.x == 80. / (77 + 80)); // 4/(3+4) CHECK(eta.x == 80. / (77 + 80)); // 4/(3+4)
CHECK(eta.y == 91.0 / (91 + 80)); // 7/(7+4) CHECK(eta.y == 91.0 / (91 + 80)); // 7/(7+4)
CHECK(eta.sum == 77 + 80 + 82 + 91); CHECK(eta.sum == 77 + 80 + 82 + 91);