detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-02-19 20:48:40 +01:00
Code Issues Actions Packages Releases Activity

Dev/rosenblatttransform (#241)

Browse Source 
- added rosenblatttransform 
- added 3x3 eta methods 
- interpolation can be used with various eta functions
- added documentation for interpolation, eta calculation 
- exposed full eta struct in python 
- disable ClusterFinder for 2x2 clusters 
- factory function for ClusterVector

---------

Co-authored-by: Dhanya Thattil <dhanya.thattil@psi.ch>
Co-authored-by: Erik Fröjdh <erik.frojdh@psi.ch>
This commit is contained in:
mazzol_a
2025-11-21 14:48:46 +01:00
committed by GitHub
parent 7fb500c44c
commit 267ca87ab0
49 changed files with 3253 additions and 1172 deletions
Download Patch File Download Diff File Expand all files Collapse all files

426
include/aare/CalculateEta.hpp
View File

@@ -19,37 +19,124 @@ enum class pixel : int {
pTopRight = 8
};
// TODO: better to have sum after x,y
/**
* eta struct
*/
template <typename T> struct Eta2 {
double x;
double y;
/// @brief eta in x direction
double x{};
/// @brief eta in y direction
double y{};
/// @brief index of subcluster given as corner relative to cluster center
corner c{0};
T sum;
/// @brief photon energy (cluster sum)
T sum{};
};
/**
* @brief Calculate the eta2 values for all clusters in a Clustervector
* @brief Calculate the eta2 values for all clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
NDArray<double, 2> calculate_eta2(const ClusterVector<ClusterType> &clusters) {
NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
std::vector<Eta2<typename ClusterType::value_type>>
calculate_eta2(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters[i]);
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate the full eta2 values for all clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_full_eta2(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_full_eta2(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate eta3 for all 3x3 clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_eta3(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta3(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate cross eta3 for all 3x3 clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_cross_eta3(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_cross_eta3(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief helper function to calculate eta2 x and y values
* @param eta reference to the Eta2 object to update
* @param left_x value of the left pixel
* @param right_x value of the right pixel
* @param bottom_y value of the bottom pixel
* @param top_y value of the top pixel
*/
template <typename T>
inline void calculate_eta2(Eta2<T> &eta, const T left_x, const T right_x,
const T bottom_y, const T top_y) {
if ((right_x + left_x) != 0)
eta.x = static_cast<double>(right_x) /
static_cast<double>(right_x + left_x); // between (0,1) the
// closer to zero left
// value probably larger
if ((top_y + bottom_y) != 0)
eta.y = static_cast<double>(top_y) /
static_cast<double>(top_y + bottom_y); // between (0,1) the
// closer to zero bottom
// value probably larger
}
/**
* @brief Calculate the eta2 values for a generic sized cluster and return them
* in a Eta2 struct containing etay, etax and the index of the respective 2x2
* subcluster.
* in a Eta2 struct containing etay, etax and the index (as corner) of the
* respective 2x2 subcluster relative to the cluster center.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
typename CoordType = uint16_t>
Eta2<T>
calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
@@ -66,67 +153,36 @@ calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
// subcluster top right from center
switch (c) {
case (corner::cTopLeft):
if ((cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]) != 0)
eta.x = static_cast<double>(cl.data[cluster_center_index - 1]) /
static_cast<double>(cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]);
if ((cl.data[cluster_center_index - ClusterSizeX] +
cl.data[cluster_center_index]) != 0)
eta.y = static_cast<double>(
cl.data[cluster_center_index - ClusterSizeX]) /
static_cast<double>(
cl.data[cluster_center_index - ClusterSizeX] +
cl.data[cluster_center_index]);
// dx = 0
// dy = 0
calculate_eta2(eta, cl.data[cluster_center_index - 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index - ClusterSizeX],
cl.data[cluster_center_index]);
// dx = -1
// dy = -1
break;
case (corner::cTopRight):
if (cl.data[cluster_center_index] + cl.data[cluster_center_index + 1] !=
0)
eta.x = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>(cl.data[cluster_center_index] +
cl.data[cluster_center_index + 1]);
if ((cl.data[cluster_center_index - ClusterSizeX] +
cl.data[cluster_center_index]) != 0)
eta.y = static_cast<double>(
cl.data[cluster_center_index - ClusterSizeX]) /
static_cast<double>(
cl.data[cluster_center_index - ClusterSizeX] +
cl.data[cluster_center_index]);
// dx = 1
// dy = 0
calculate_eta2(eta, cl.data[cluster_center_index],
cl.data[cluster_center_index + 1],
cl.data[cluster_center_index - ClusterSizeX],
cl.data[cluster_center_index]);
// dx = 0
// dy = -1
break;
case (corner::cBottomLeft):
if ((cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]) != 0)
eta.x = static_cast<double>(cl.data[cluster_center_index - 1]) /
static_cast<double>(cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]);
if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]) != 0)
eta.y = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>(
cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]);
// dx = 0
// dy = 1
calculate_eta2(eta, cl.data[cluster_center_index - 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index],
cl.data[cluster_center_index + ClusterSizeX]);
// dx = -1
// dy = 0
break;
case (corner::cBottomRight):
if (cl.data[cluster_center_index] + cl.data[cluster_center_index + 1] !=
0)
eta.x = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>(cl.data[cluster_center_index] +
cl.data[cluster_center_index + 1]);
if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]) != 0)
eta.y = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>(
cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]);
// dx = 1
// dy = 1
calculate_eta2(eta, cl.data[cluster_center_index],
cl.data[cluster_center_index + 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index + ClusterSizeX]);
// dx = 0
// dy = 0
break;
}
@@ -135,69 +191,255 @@ calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
return eta;
}
// TODO! Look up eta2 calculation - photon center should be bottom right corner
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) {
/**
* @brief Calculate the eta2 values for a generic sized cluster and return them
* in a Eta2 struct containing etay, etax and the index (as corner) of the
* respective 2x2 subcluster relative to the cluster center.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
Eta2<T> calculate_full_eta2(
const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 1 && ClusterSizeY > 1);
Eta2<T> eta{};
if ((cl.data[0] + cl.data[1]) != 0)
eta.x = static_cast<double>(cl.data[2]) /
(cl.data[2] + cl.data[3]); // 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[1]) /
(cl.data[1] + cl.data[3]); // between (0,1) the closer to zero
// bottom value probably larger
constexpr size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
auto max_sum = cl.max_sum_2x2();
eta.sum = max_sum.sum;
corner c = max_sum.index;
// subcluster top right from center
switch (c) {
case (corner::cTopLeft):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index] +
cl.data[cluster_center_index - ClusterSizeX]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]) /
static_cast<double>(eta.sum);
}
// dx = -1
// dy = -1
break;
case (corner::cTopRight):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index - ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[cluster_center_index] +
cl.data[cluster_center_index + 1]) /
static_cast<double>(eta.sum);
}
// dx = 0
// dy = -1
break;
case (corner::cBottomLeft):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(
cl.data[cluster_center_index + ClusterSizeX] +
cl.data[cluster_center_index + ClusterSizeX - 1]) /
static_cast<double>(eta.sum);
}
// dx = -1
// dy = 0
break;
case (corner::cBottomRight):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index + ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(
cl.data[cluster_center_index + ClusterSizeX] +
cl.data[cluster_center_index + ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
}
// dx = 0
// dy = 0
break;
}
eta.c = c;
return eta;
}
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
Eta2<T> eta{};
// TODO: maybe have as member function of cluster
const uint8_t photon_hit_index =
std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
eta.c = static_cast<corner>(3 - photon_hit_index);
switch (eta.c) {
case corner::cTopLeft:
calculate_eta2(eta, cl.data[2], cl.data[3], cl.data[1], cl.data[3]);
break;
case corner::cTopRight:
calculate_eta2(eta, cl.data[2], cl.data[3], cl.data[0], cl.data[2]);
break;
case corner::cBottomLeft:
calculate_eta2(eta, cl.data[0], cl.data[1], cl.data[1], cl.data[3]);
break;
case corner::cBottomRight:
calculate_eta2(eta, cl.data[0], cl.data[1], cl.data[0], cl.data[2]);
break;
}
eta.sum = cl.sum();
return eta;
}
template <typename T>
Eta2<T> calculate_full_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
Eta2<T> eta{};
eta.sum = cl.sum();
const uint8_t photon_hit_index =
std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
eta.c = static_cast<corner>(3 - photon_hit_index);
if (eta.sum != 0) {
eta.x = static_cast<double>(cl.data[1] + cl.data[3]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[2] + cl.data[3]) /
static_cast<double>(eta.sum);
}
return eta;
}
// TODO generalize
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 1, 2, int16_t> &cl) {
Eta2<T> calculate_eta2(const Cluster<T, 1, 2, uint16_t> &cl) {
Eta2<T> eta{};
eta.x = 0;
eta.y = static_cast<double>(cl.data[0]) / cl.data[1];
eta.y = static_cast<double>(cl.data[1]) / cl.data[0];
eta.sum = cl.sum();
}
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 1, int16_t> &cl) {
Eta2<T> calculate_eta2(const Cluster<T, 2, 1, uint16_t> &cl) {
Eta2<T> eta{};
eta.x = static_cast<double>(cl.data[0]) / cl.data[1];
eta.x = static_cast<double>(cl.data[1]) / cl.data[0];
eta.y = 0;
eta.sum = cl.sum();
}
// calculates Eta3 for 3x3 cluster based on code from analyze_cluster
// TODO only supported for 3x3 Clusters
template <typename T> Eta2<T> calculate_eta3(const Cluster<T, 3, 3> &cl) {
/**
* @brief calculates cross Eta3 for 3x3 cluster
* cross Eta3 calculates the eta by taking into account only the cross pixels
* {top, bottom, left, right, center}
*/
template <typename T, typename CoordType = uint16_t>
Eta2<T> calculate_cross_eta3(const Cluster<T, 3, 3, CoordType> &cl) {
Eta2<T> eta{};
T sum = 0;
T photon_energy = cl.sum();
std::for_each(std::begin(cl.data), std::end(cl.data),
[&sum](T x) { sum += x; });
eta.sum = sum;
eta.sum = photon_energy;
if ((cl.data[3] + cl.data[4] + cl.data[5]) != 0)
eta.x = static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
eta.x =
static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
(cl.data[3] + cl.data[4] + cl.data[5]); // (-1,1)
static_cast<double>(cl.data[3] + cl.data[4] + cl.data[5]); // (-1,1)
if ((cl.data[1] + cl.data[4] + cl.data[7]) != 0)
eta.y = static_cast<double>(-cl.data[1] + cl.data[2 * 3 + 1]) /
(cl.data[1] + cl.data[4] + cl.data[7]);
static_cast<double>(cl.data[1] + cl.data[4] + cl.data[7]);
return eta;
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
Eta2<T> calculate_cross_eta3(
const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 2 && ClusterSizeY > 2,
"calculate_eta3 only defined for clusters larger than 2x2");
if constexpr (ClusterSizeX != 3 || ClusterSizeY != 3) {
auto reduced_cluster = reduce_cluster_to_3x3(cl);
return calculate_cross_eta3(reduced_cluster);
} else {
return calculate_cross_eta3(cl);
}
}
/**
* @brief calculates Eta3 for 3x3 cluster
* It calculates the eta by taking into account all pixels in the 3x3 cluster
*/
template <typename T, typename CoordType = uint16_t>
Eta2<T> calculate_eta3(const Cluster<T, 3, 3, CoordType> &cl) {
Eta2<T> eta{};
T photon_energy = cl.sum();
eta.sum = photon_energy;
// TODO: how do we handle potential arithmetic overflows? - T could be
// uint16
if (photon_energy != 0) {
std::array<T, 2> column_sums{
static_cast<T>(cl.data[0] + cl.data[3] + cl.data[6]),
static_cast<T>(cl.data[2] + cl.data[5] + cl.data[8])};
eta.x = static_cast<double>(-column_sums[0] + column_sums[1]) /
static_cast<double>(photon_energy);
std::array<T, 2> row_sums{
static_cast<T>(cl.data[0] + cl.data[1] + cl.data[2]),
static_cast<T>(cl.data[6] + cl.data[7] + cl.data[8])};
eta.y = static_cast<double>(-row_sums[0] + row_sums[1]) /
static_cast<double>(photon_energy);
}
return eta;
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
Eta2<T>
calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 2 && ClusterSizeY > 2,
"calculate_eta3 only defined for clusters larger than 2x2");
if constexpr (ClusterSizeX != 3 || ClusterSizeY != 3) {
auto reduced_cluster = reduce_cluster_to_3x3(cl);
return calculate_eta3(reduced_cluster);
} else {
return calculate_eta3(cl);
}
}
} // namespace aare

147
include/aare/Cluster.hpp
View File

@@ -18,6 +18,10 @@
namespace aare {
// requires clause c++20 maybe update
/**
* @brief Cluster struct
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
struct Cluster {
@@ -28,8 +32,11 @@ struct Cluster {
static_assert(ClusterSizeX > 0 && ClusterSizeY > 0,
"Cluster sizes must be bigger than zero");
/// @brief Cluster center x coordinate (in pixel coordinates)
CoordType x;
/// @brief Cluster center y coordinate (in pixel coordinates)
CoordType y;
/// @brief Cluster data stored in row-major order starting from top-left
std::array<T, ClusterSizeX * ClusterSizeY> data;
static constexpr uint8_t cluster_size_x = ClusterSizeX;
@@ -37,10 +44,12 @@ struct Cluster {
using value_type = T;
using coord_type = CoordType;
/**
* @brief Sum of all elements in the cluster
*/
T sum() const { return std::accumulate(data.begin(), data.end(), T{}); }
// TODO: handle 1 dimensional clusters
// TODO: change int to corner
/**
* @brief sum of 2x2 subcluster with highest energy
* @return photon energy of subcluster, 2x2 subcluster index relative to
@@ -111,66 +120,71 @@ struct Cluster {
* highest sum.
* @param c Cluster to reduce
* @return reduced cluster
* @note The cluster is filled using row major ordering starting at the top-left
* (thus for a max subcluster in the top left cornern the photon hit is at
* the fourth position)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
typename CoordType = uint16_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;
Cluster<T, 2, 2, CoordType> result{};
auto [sum, index] = c.max_sum_2x2();
int16_t cluster_center_index =
constexpr int16_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
int16_t index_bottom_left_max_2x2_subcluster =
(int(static_cast<int>(index) / (ClusterSizeX - 1))) * ClusterSizeX +
static_cast<int>(index) % (ClusterSizeX - 1);
int16_t index_top_left_max_2x2_subcluster = cluster_center_index;
switch (index) {
case corner::cTopLeft:
index_top_left_max_2x2_subcluster -= (ClusterSizeX + 1);
break;
case corner::cTopRight:
index_top_left_max_2x2_subcluster -= ClusterSizeX;
break;
case corner::cBottomLeft:
index_top_left_max_2x2_subcluster -= 1;
break;
case corner::cBottomRight:
// no change needed
break;
}
result.x =
c.x + (index_bottom_left_max_2x2_subcluster - cluster_center_index) %
ClusterSizeX;
result.x = c.x;
result.y = c.y;
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]};
c.data[index_top_left_max_2x2_subcluster],
c.data[index_top_left_max_2x2_subcluster + 1],
c.data[index_top_left_max_2x2_subcluster + ClusterSizeX],
c.data[index_top_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;
Cluster<T, 2, 2, uint16_t> reduce_to_2x2(const Cluster<T, 3, 3, uint16_t> &c) {
Cluster<T, 2, 2, uint16_t> result{};
auto [s, i] = c.max_sum_2x2();
result.x = c.x;
result.y = c.y;
switch (i) {
case corner::cTopLeft:
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 corner::cTopRight:
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 corner::cBottomLeft:
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 corner::cBottomRight:
result.x = c.x;
result.y = c.y;
result.data = {c.data[4], c.data[5], c.data[7], c.data[8]};
break;
}
@@ -178,43 +192,8 @@ Cluster<T, 2, 2, int16_t> reduce_to_2x2(const Cluster<T, 3, 3, int16_t> &c) {
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.
* @brief Reduce a cluster to a 3x3 cluster
* @param c Cluster to reduce
* @return reduced cluster
*/
@@ -226,40 +205,24 @@ 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);
Cluster<T, 3, 3, CoordType> result{};
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;
result.x = c.x;
result.y = c.y;
int16_t index_3x3_subcluster_cluster_center =
int((cluster_center_index - 1 - ClusterSizeX) / ClusterSizeX) *
(ClusterSizeX - 2) +
(cluster_center_index - 1 - ClusterSizeX) % ClusterSizeX;
result.data = {c.data[cluster_center_index - ClusterSizeX - 1],
c.data[cluster_center_index - ClusterSizeX],
c.data[cluster_center_index - ClusterSizeX + 1],
c.data[cluster_center_index - 1],
c.data[cluster_center_index],
c.data[cluster_center_index + 1],
c.data[cluster_center_index + ClusterSizeX - 1],
c.data[cluster_center_index + ClusterSizeX],
c.data[cluster_center_index + ClusterSizeX + 1]};
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;
}

3
include/aare/ClusterFileSink.hpp
View File

@@ -10,7 +10,8 @@
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
typename = std::enable_if_t<is_cluster_v<ClusterType>>,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFileSink {
ProducerConsumerQueue<ClusterVector<ClusterType>> *m_source;
std::atomic<bool> m_stop_requested{false};

10
include/aare/ClusterFinder.hpp
View File

@@ -10,8 +10,16 @@
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
struct no_2x2_cluster {
constexpr static bool value =
ClusterType::cluster_size_x > 2 && ClusterType::cluster_size_y > 2;
};
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFinder {
Shape<2> m_image_size;
const PEDESTAL_TYPE m_nSigma;

3
include/aare/ClusterFinderMT.hpp
View File

@@ -32,7 +32,8 @@ struct FrameWrapper {
* @tparam CT type of the cluster data
*/
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFinderMT {
protected:

12
include/aare/ClusterVector.hpp
View File

@@ -28,7 +28,7 @@ class ClusterVector; // Forward declaration
* needed.
* @tparam T data type of the pixels in the cluster
* @tparam CoordType data type of the x and y coordinates of the cluster
* (normally int16_t)
* (normally uint16_t)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
@@ -176,9 +176,12 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
* highest sum.
* @param cv Clustervector containing clusters to reduce
* @return Clustervector with reduced clusters
* @note The cluster is filled using row major ordering starting at the top-left
* (thus for a max subcluster in the top left cornern the photon hit is at
* the fourth position)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
typename CoordType>
ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
&cv) {
@@ -190,13 +193,12 @@ ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
}
/**
* @brief Reduce a cluster to a 3x3 cluster by selecting the 3x3 block with the
* highest sum.
* @brief Reduce a cluster to a 3x3 cluster
* @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>
typename CoordType>
ClusterVector<Cluster<T, 3, 3, CoordType>> reduce_to_3x3(
const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
&cv) {

257
include/aare/Interpolator.hpp
View File

@@ -17,7 +17,10 @@ struct Photon {
};
class Interpolator {
// marginal CDF of eta_x (if rosenblatt applied), conditional
// CDF of eta_x conditioned on eta_y
NDArray<double, 3> m_ietax;
// conditional CDF of eta_y conditioned on eta_x
NDArray<double, 3> m_ietay;
NDArray<double, 1> m_etabinsx;
@@ -25,108 +28,210 @@ class Interpolator {
NDArray<double, 1> m_energy_bins;
public:
/**
* @brief Constructor for the Interpolator class
* @param etacube joint distribution of etaX, etaY and photon energy
* @param xbins bin edges for etaX
* @param ybins bin edges for etaY
* @param ebins bin edges for photon energy
* @note note first dimension is etaX, second etaY, third photon energy
*/
Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
NDView<double, 1> ybins, NDView<double, 1> ebins);
/**
* @brief Constructor for the Interpolator class
* @param xbins bin edges for etaX
* @param ybins bin edges for etaY
* @param ebins bin edges for photon energy
*/
Interpolator(NDView<double, 1> xbins, NDView<double, 1> ybins,
NDView<double, 1> ebins);
/**
* @brief transforms the joint eta distribution of etaX and etaY to the two
* independant uniform distributions based on the Roseblatt transform for
* each energy level
* @param etacube joint distribution of etaX, etaY and photon energy
* @note note first dimension is etaX, second etaY, third photon energy
*/
void rosenblatttransform(NDView<double, 3> etacube);
NDArray<double, 3> get_ietax() { return m_ietax; }
NDArray<double, 3> get_ietay() { return m_ietay; }
template <typename ClusterType,
/**
* @brief interpolates the cluster centers for all clusters to a better
* precision
* @tparam ClusterType Type of Clusters to interpolate
* @tparam Etafunction Function object that calculates desired eta default:
* calculate_eta2
* @return interpolated photons (photon positions are given as double but
* following row column format e.g. x=0, y=0 means top row and first column
* of frame)
*/
template <auto EtaFunction = calculate_eta2, typename ClusterType,
typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
private:
/**
* @brief implements underlying interpolation logic based on EtaFunction
* Type
* @tparam EtaFunction Function object that calculates desired eta default:
* @param u: transformed photon position in x between [0,1]
* @param v: transformed photon position in y between [0,1]
* @param c: corner of eta
*/
template <auto EtaFunction, typename ClusterType>
void interpolation_logic(Photon &photon, const double u, const double v,
const corner c = corner::cTopLeft);
/**
* @brief bilinear interpolation of the transformed eta values
* @param ix index of etaX bin
* @param iy index of etaY bin
* @param ie index of energy bin
* @return pair of interpolated transformed eta values (ietax, ietay)
*/
template <typename T>
std::pair<double, double>
bilinear_interpolation(const size_t ix, const size_t iy, const size_t ie,
const Eta2<T> &eta);
};
// TODO: generalize to support any clustertype!!! otherwise add std::enable_if_t
// to only take Cluster2x2 and Cluster3x3
template <typename ClusterType, typename Enable>
template <typename T>
std::pair<double, double>
Interpolator::bilinear_interpolation(const size_t ix, const size_t iy,
const size_t ie, const Eta2<T> &eta) {
auto next_index_y = static_cast<ssize_t>(iy + 1) >= m_ietax.shape(1)
? m_ietax.shape(1) - 1
: iy + 1;
auto next_index_x = static_cast<ssize_t>(ix + 1) >= m_ietax.shape(0)
? m_ietax.shape(0) - 1
: ix + 1;
// bilinear interpolation
double ietax_interp_left = linear_interpolation(
{m_etabinsy(iy), m_etabinsy(iy + 1)},
{m_ietax(ix, iy, ie), m_ietax(ix, next_index_y, ie)}, eta.y);
double ietax_interp_right =
linear_interpolation({m_etabinsy(iy), m_etabinsy(iy + 1)},
{m_ietax(next_index_x, iy, ie),
m_ietax(next_index_x, next_index_y, ie)},
eta.y);
// transformed photon position x between [0,1]
double ietax_interpolated =
linear_interpolation({m_etabinsx(ix), m_etabinsx(ix + 1)},
{ietax_interp_left, ietax_interp_right}, eta.x);
double ietay_interp_left = linear_interpolation(
{m_etabinsx(ix), m_etabinsx(ix + 1)},
{m_ietay(ix, iy, ie), m_ietay(next_index_x, iy, ie)}, eta.x);
double ietay_interp_right =
linear_interpolation({m_etabinsx(ix), m_etabinsx(ix + 1)},
{m_ietay(ix, next_index_y, ie),
m_ietay(next_index_x, next_index_y, ie)},
eta.x);
// transformed photon position y between [0,1]
double ietay_interpolated =
linear_interpolation({m_etabinsy(iy), m_etabinsy(iy + 1)},
{ietay_interp_left, ietay_interp_right}, eta.y);
return {ietax_interpolated, ietay_interpolated};
}
template <auto EtaFunction, typename ClusterType, typename Enable>
std::vector<Photon>
Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
std::vector<Photon> photons;
photons.reserve(clusters.size());
if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
for (const ClusterType &cluster : clusters) {
for (const ClusterType &cluster : clusters) {
auto eta = calculate_eta2(cluster);
auto eta = EtaFunction(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
// auto ie = nearest_index(m_energy_bins, photon.energy)-1;
// auto ix = nearest_index(m_etabinsx, eta.x)-1;
// auto iy = nearest_index(m_etabinsy, eta.y)-1;
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
// std::cout << "eta.x: " << eta.x << " eta.y: " << eta.y << std::endl;
// fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
double dX, dY;
// cBottomLeft = 0,
// cBottomRight = 1,
// cTopLeft = 2,
// cTopRight = 3
// TODO: could also chaneg the sign of the eta calculation
switch (static_cast<corner>(eta.c)) {
case corner::cTopLeft:
dX = 0.0;
dY = 0.0;
break;
case corner::cTopRight:;
dX = 1.0;
dY = 0.0;
break;
case corner::cBottomLeft:
dX = 0.0;
dY = 1.0;
break;
case corner::cBottomRight:
dX = 1.0;
dY = 1.0;
break;
}
photon.x -= m_ietax(ix, iy, ie) - dX;
photon.y -= m_ietay(ix, iy, ie) - dY;
photons.push_back(photon);
}
} else if (clusters.cluster_size_x() == 2 ||
clusters.cluster_size_y() == 2) {
for (const ClusterType &cluster : clusters) {
auto eta = calculate_eta2(cluster);
// std::cout << "ix: " << ix << " iy: " << iy << std::endl;
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
// TODO: bilinear interpolation only works if all bins have a size > 1 -
// otherwise bilinear interpolation with zero values which skew the
// results
// TODO: maybe trim the bins at the edges with zero values beforehand
// auto [ietax_interpolated, ietay_interpolated] =
// bilinear_interpolation(ix, iy, ie, eta);
// Now do some actual interpolation.
// Find which energy bin the cluster is in
// auto ie = nearest_index(m_energy_bins, photon.energy)-1;
// auto ix = nearest_index(m_etabinsx, eta.x)-1;
// auto iy = nearest_index(m_etabinsy, eta.y)-1;
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
double ietax_interpolated = m_ietax(ix, iy, ie);
double ietay_interpolated = m_ietay(ix, iy, ie);
// TODO: why 2?
photon.x -=
m_ietax(ix, iy, ie); // eta goes between 0 and 1 but we could
// move the hit anywhere in the 2x2
photon.y -= m_ietay(ix, iy, ie);
photons.push_back(photon);
}
interpolation_logic<EtaFunction, ClusterType>(
photon, ietax_interpolated, ietay_interpolated, eta.c);
} else {
throw std::runtime_error(
"Only 3x3 and 2x2 clusters are supported for interpolation");
photons.push_back(photon);
}
return photons;
}
template <auto EtaFunction, typename ClusterType>
void Interpolator::interpolation_logic(Photon &photon, const double u,
const double v, const corner c) {
// std::cout << "u: " << u << " v: " << v << std::endl;
// TODO: try to call this with std::is_same_v and have it constexpr if
// possible
if (EtaFunction == &calculate_eta2<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type> ||
EtaFunction == &calculate_full_eta2<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type>) {
double dX{}, dY{};
// TODO: could also chaneg the sign of the eta calculation
switch (c) {
case corner::cTopLeft:
dX = -1.0;
dY = -1.0;
break;
case corner::cTopRight:;
dX = 0.0;
dY = -1.0;
break;
case corner::cBottomLeft:
dX = -1.0;
dY = 0.0;
break;
case corner::cBottomRight:
dX = 0.0;
dY = 0.0;
break;
}
photon.x = photon.x + 0.5 + u + dX; // use pixel center + 0.5
photon.y = photon.y + 0.5 + v +
dY; // eta2 calculates the ratio between bottom and sum of
// bottom and top shift by 1 add eta value correctly
} else {
photon.x += u;
photon.y += v;
}
}
} // namespace aare

15
include/aare/algorithm.hpp
View File

@@ -109,4 +109,19 @@ template <typename Container> bool all_equal(const Container &c) {
return false;
}
/**
* linear interpolation
* @param bin_edge left and right bin edges
* @param bin_values function values at bin edges
* @param coord coordinate to interpolate at
* @return interpolated value at coord
*/
inline double linear_interpolation(const std::pair<double, double> &bin_edge,
const std::pair<double, double> &bin_values,
const double coord) {
const double bin_width = bin_edge.second - bin_edge.first;
return bin_values.first * (1 - (coord - bin_edge.first) / bin_width) +
bin_values.second * (coord - bin_edge.first) / bin_width;
}
} // namespace aare