diff --git a/CMakeLists.txt b/CMakeLists.txt index 51ed7f5..efb3a0c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 "$" PRIVATE ${lmfit_SOURCE_DIR}/lib ) - - target_link_libraries( aare_core PUBLIC diff --git a/benchmarks/calculateeta_benchmark.cpp b/benchmarks/calculateeta_benchmark.cpp index dc7cd91..609ce89 100644 --- a/benchmarks/calculateeta_benchmark.cpp +++ b/benchmarks/calculateeta_benchmark.cpp @@ -1,3 +1,4 @@ +#include "aare/CalculateEta.hpp" #include "aare/ClusterFile.hpp" #include diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index e9530f6..289647d 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -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 &cl, int32_t *t2, int32_t *t3, - char *quad, double *eta2x, double *eta2y, double *eta3x, - double *eta3y); - -// template >> -// NDArray calculate_eta2(ClusterVector &clusters); - -// TODO: do we need rquire clauses? -// template Eta2 calculate_eta2(const Cluster &cl); - -// template Eta2 calculate_eta2(const Cluster &cl); - -// template >> -// Eta2 calculate_eta2(const ClusterType &cl); - -// template -// Eta2 calculate_eta2( -// const Cluster &cl); - template ClusterFile::ClusterFile( const std::filesystem::path &fname, size_t chunk_size, @@ -374,160 +317,4 @@ ClusterVector ClusterFile::read_frame() { return clusters; } -template >> -NDArray calculate_eta2(const ClusterVector &clusters) { - // TOTO! make work with 2x2 clusters - NDArray eta2({static_cast(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 -Eta2 calculate_eta2( - const Cluster &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 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( - 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( - 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 Eta2 calculate_eta2(const Cluster &cl) { - Eta2 eta{}; - - std::array 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(cl.data[4]) / (cl.data[3] + cl.data[4]); - if ((cl.data[1] + cl.data[4]) != 0) - eta.y = static_cast(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(cl.data[5]) / (cl.data[4] + cl.data[5]); - if ((cl.data[1] + cl.data[4]) != 0) - eta.y = static_cast(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(cl.data[4]) / (cl.data[3] + cl.data[4]); - if ((cl.data[7] + cl.data[4]) != 0) - eta.y = static_cast(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(cl.data[5]) / (cl.data[5] + cl.data[4]); - if ((cl.data[7] + cl.data[4]) != 0) - eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); - eta.c = cTopRight; - break; - } - return eta; -} - -template Eta2 calculate_eta2(const Cluster &cl) { - Eta2 eta{}; - - eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); - eta.y = static_cast(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 Eta2 calculate_eta3(const Cluster &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(-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(-cl.data[1] + cl.data[2 * 3 + 1]) / - - (cl.data[1] + cl.data[4] + cl.data[7]); - - return eta; -} - } // namespace aare diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index ec0fa40..0e47b51 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -30,8 +30,7 @@ template class ClusterVector> { 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> { * @throws std::runtime_error if the cluster size is not 3x3 * @warning Only 3x3 clusters are supported for the 2x2 sum. */ - std::vector sum_2x2() { - std::vector sums(m_size); - const size_t stride = item_size(); + /* only needed to calculate eta + std::vector sum_2x2() { + std::vector 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 total; + auto T_ptr = reinterpret_cast(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 total; - auto T_ptr = reinterpret_cast(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 diff --git a/src/Interpolator.cpp b/src/Interpolator.cpp index d95405a..1c4a385 100644 --- a/src/Interpolator.cpp +++ b/src/Interpolator.cpp @@ -1,4 +1,5 @@ #include "aare/Interpolator.hpp" +#include "aare/CalculateEta.hpp" #include "aare/algorithm.hpp" namespace aare {