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

refactored and put calculate_eta function in seperate file

Browse Source
This commit is contained in:
Mazzoleni Alice Francesca 2025-03-31 17:35:39 +02:00
parent 7e5f91c6ec
commit e038bd1646
5 changed files with 28 additions and 240 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists.txt
View File

@ -307,6 +307,7 @@ endif()
set(PUBLICHEADERS
include/aare/ArrayExpr.hpp
include/aare/CalculateEta.hpp
include/aare/Cluster.hpp
#include/aare/ClusterFinder.hpp
include/aare/ClusterFile.hpp
@ -331,7 +332,6 @@ set(PUBLICHEADERS
include/aare/RawSubFile.hpp
#include/aare/VarClusterFinder.hpp
include/aare/utils/task.hpp
)
@ -362,8 +362,6 @@ target_include_directories(aare_core PUBLIC
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>" PRIVATE ${lmfit_SOURCE_DIR}/lib
)
target_link_libraries(
aare_core
PUBLIC

1
benchmarks/calculateeta_benchmark.cpp
View File

@ -1,3 +1,4 @@
#include "aare/CalculateEta.hpp"
#include "aare/ClusterFile.hpp"
#include <benchmark/benchmark.h>

213
include/aare/ClusterFile.hpp
View File

@ -9,40 +9,6 @@
namespace aare {
typedef enum {
cBottomLeft = 0,
cBottomRight = 1,
cTopLeft = 2,
cTopRight = 3
} corner;
typedef enum {
pBottomLeft = 0,
pBottom = 1,
pBottomRight = 2,
pLeft = 3,
pCenter = 4,
pRight = 5,
pTopLeft = 6,
pTop = 7,
pTopRight = 8
} pixel;
// TODO: maybe template this!!!!!! why int32_t????
struct Eta2 {
double x;
double y;
int c;
int32_t sum;
};
struct ClusterAnalysis {
uint32_t c;
int32_t tot;
double etax;
double etay;
};
/*
Binary cluster file. Expects data to be layed out as:
int32_t frame_number
@ -126,29 +92,6 @@ class ClusterFile {
void close();
};
int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);
int analyze_cluster(Cluster<int32_t, 3, 3> &cl, int32_t *t2, int32_t *t3,
char *quad, double *eta2x, double *eta2y, double *eta3x,
double *eta3y);
// template <typename ClusterType,
// typename = std::enable_if_t<is_cluster_v<ClusterType>>>
// NDArray<double, 2> calculate_eta2(ClusterVector<ClusterType> &clusters);
// TODO: do we need rquire clauses?
// template <typename T> Eta2 calculate_eta2(const Cluster<T, 3, 3> &cl);
// template <typename T> Eta2 calculate_eta2(const Cluster<T, 2, 2> &cl);
// template <typename ClusterType, std::enable_if_t<is_cluster_v<ClusterType>>>
// Eta2 calculate_eta2(const ClusterType &cl);
// template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
// typename CoordType>
// Eta2 calculate_eta2(
// const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl);
template <typename ClusterType, typename Enable>
ClusterFile<ClusterType, Enable>::ClusterFile(
const std::filesystem::path &fname, size_t chunk_size,
@ -374,160 +317,4 @@ ClusterVector<ClusterType> ClusterFile<ClusterType, Enable>::read_frame() {
return clusters;
}
template <typename ClusterType, std::enable_if_t<is_cluster_v<ClusterType>>>
NDArray<double, 2> calculate_eta2(const ClusterVector<ClusterType> &clusters) {
// TOTO! make work with 2x2 clusters
NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters.at(i));
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
}
return eta2;
}
/**
* @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.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
Eta2 calculate_eta2(
const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
Eta2 eta{};
// TODO loads of overhead for a 2x2 clsuter maybe keep 2x2 calculation
constexpr size_t num_2x2_subclusters =
(ClusterSizeX - 1) * (ClusterSizeY - 1);
std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
for (size_t j = 0; j < ClusterSizeX - 1; ++j)
sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] =
cl.data[i * ClusterSizeX + j] +
cl.data[i * ClusterSizeX + j + 1] +
cl.data[(i + 1) * ClusterSizeX + j] +
cl.data[(i + 1) * ClusterSizeX + j + 1];
}
auto c =
std::max_element(sum_2x2_subcluster.begin(), sum_2x2_subcluster.end()) -
sum_2x2_subcluster.begin();
eta.sum = sum_2x2_subcluster[c];
size_t index_bottom_left_max_2x2_subcluster =
(int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1);
if ((cl.data[index_bottom_left_max_2x2_subcluster] +
cl.data[index_bottom_left_max_2x2_subcluster + 1]) != 0)
eta.x = static_cast<double>(
cl.data[index_bottom_left_max_2x2_subcluster + 1]) /
(cl.data[index_bottom_left_max_2x2_subcluster] +
cl.data[index_bottom_left_max_2x2_subcluster + 1]);
if ((cl.data[index_bottom_left_max_2x2_subcluster] +
cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]) != 0)
eta.y =
static_cast<double>(
cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]) /
(cl.data[index_bottom_left_max_2x2_subcluster] +
cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]);
eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no
// underyling enum class
return eta;
}
/**
* @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2
* struct containing etay, etax and the corner of the cluster.
*/
template <typename T> Eta2 calculate_eta2(const Cluster<T, 3, 3> &cl) {
Eta2 eta{};
std::array<T, 4> tot2;
tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4];
tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5];
tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7];
tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8];
auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin();
eta.sum = tot2[c];
switch (c) {
case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0)
eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[1] + cl.data[4]) != 0)
eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.c = cBottomLeft;
break;
case cBottomRight:
if ((cl.data[2] + cl.data[5]) != 0)
eta.x = static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
if ((cl.data[1] + cl.data[4]) != 0)
eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.c = cBottomRight;
break;
case cTopLeft:
if ((cl.data[7] + cl.data[4]) != 0)
eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopLeft;
break;
case cTopRight:
if ((cl.data[5] + cl.data[4]) != 0)
eta.x = static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopRight;
break;
}
return eta;
}
template <typename T> Eta2 calculate_eta2(const Cluster<T, 2, 2> &cl) {
Eta2 eta{};
eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
eta.sum = cl.data[0] + cl.data[1] + cl.data[2] + cl.data[3];
eta.c = cBottomLeft; // TODO! This is not correct, but need to put something
return eta;
}
// calculates Eta3 for 3x3 cluster based on code from analyze_cluster
// TODO only supported for 3x3 Clusters
template <typename T> Eta2 calculate_eta3(const Cluster<T, 3, 3> &cl) {
Eta2 eta{};
T sum = 0;
std::for_each(std::begin(cl.data), std::end(cl.data),
[&sum](T x) { sum += x; });
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]);
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]);
return eta;
}
} // namespace aare

49
include/aare/ClusterVector.hpp
View File

@ -30,8 +30,7 @@ template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
using value_type = T;
// size_t m_cluster_size_x;
// size_t m_cluster_size_y;
std::byte *m_data{};
size_t m_size{0};
size_t m_capacity;
@ -155,30 +154,32 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
* @throws std::runtime_error if the cluster size is not 3x3
* @warning Only 3x3 clusters are supported for the 2x2 sum.
*/
std::vector<T> sum_2x2() {
std::vector<T> sums(m_size);
const size_t stride = item_size();
/* only needed to calculate eta
std::vector<T> sum_2x2() {
std::vector<T> sums(m_size);
const size_t stride = item_size();
if (ClusterSizeX != 3 || ClusterSizeY != 3) {
throw std::runtime_error(
"Only 3x3 clusters are supported for the 2x2 sum.");
if (ClusterSizeX != 3 || ClusterSizeY != 3) {
throw std::runtime_error(
"Only 3x3 clusters are supported for the 2x2 sum.");
}
std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
for (size_t i = 0; i < m_size; i++) {
std::array<T, 4> total;
auto T_ptr = reinterpret_cast<T *>(ptr);
total[0] = T_ptr[0] + T_ptr[1] + T_ptr[3] + T_ptr[4];
total[1] = T_ptr[1] + T_ptr[2] + T_ptr[4] + T_ptr[5];
total[2] = T_ptr[3] + T_ptr[4] + T_ptr[6] + T_ptr[7];
total[3] = T_ptr[4] + T_ptr[5] + T_ptr[7] + T_ptr[8];
sums[i] = *std::max_element(total.begin(), total.end());
ptr += stride;
}
return sums;
}
std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
for (size_t i = 0; i < m_size; i++) {
std::array<T, 4> total;
auto T_ptr = reinterpret_cast<T *>(ptr);
total[0] = T_ptr[0] + T_ptr[1] + T_ptr[3] + T_ptr[4];
total[1] = T_ptr[1] + T_ptr[2] + T_ptr[4] + T_ptr[5];
total[2] = T_ptr[3] + T_ptr[4] + T_ptr[6] + T_ptr[7];
total[3] = T_ptr[4] + T_ptr[5] + T_ptr[7] + T_ptr[8];
sums[i] = *std::max_element(total.begin(), total.end());
ptr += stride;
}
return sums;
}
*/
/**
* @brief Return the number of clusters in the vector

1
src/Interpolator.cpp
View File

@ -1,4 +1,5 @@
#include "aare/Interpolator.hpp"
#include "aare/CalculateEta.hpp"
#include "aare/algorithm.hpp"
namespace aare {