diff --git a/include/aare/Cluster.hpp b/include/aare/Cluster.hpp index 7eb1a13..889593b 100644 --- a/include/aare/Cluster.hpp +++ b/include/aare/Cluster.hpp @@ -16,94 +16,61 @@ namespace aare { -template -constexpr bool is_valid_cluster = - std::is_arithmetic_v && std::is_integral_v && - (ClusterSizeX > 0) && (ClusterSizeY > 0); - // requires clause c++20 maybe update template >> + typename CoordType = int16_t> struct Cluster { + + static_assert(std::is_arithmetic_v, "T needs to be an arithmetic type"); + static_assert(std::is_integral_v, + "CoordType needs to be an integral type"); + static_assert(ClusterSizeX > 0 && ClusterSizeY > 0, + "Cluster sizes must be bigger than zero"); + CoordType x; CoordType y; - T data[ClusterSizeX * ClusterSizeY]; + std::array data; static constexpr uint8_t cluster_size_x = ClusterSizeX; static constexpr uint8_t cluster_size_y = ClusterSizeY; using value_type = T; using coord_type = CoordType; - T sum() const { - return std::accumulate(data, data + ClusterSizeX * ClusterSizeY, 0); - } + T sum() const { return std::accumulate(data.begin(), data.end(), T{}); } std::pair max_sum_2x2() const { - constexpr size_t num_2x2_subclusters = - (ClusterSizeX - 1) * (ClusterSizeY - 1); + if constexpr (cluster_size_x == 3 && cluster_size_y == 3) { + std::array sum_2x2_subclusters; + sum_2x2_subclusters[0] = data[0] + data[1] + data[3] + data[4]; + sum_2x2_subclusters[1] = data[1] + data[2] + data[4] + data[5]; + sum_2x2_subclusters[2] = data[3] + data[4] + data[6] + data[7]; + sum_2x2_subclusters[3] = data[4] + data[5] + data[7] + data[8]; + int index = std::max_element(sum_2x2_subclusters.begin(), + sum_2x2_subclusters.end()) - + sum_2x2_subclusters.begin(); + return std::make_pair(sum_2x2_subclusters[index], index); + } else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) { + return std::make_pair(data[0] + data[1] + data[2] + data[3], 0); + } else { + 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] = - data[i * ClusterSizeX + j] + - data[i * ClusterSizeX + j + 1] + - data[(i + 1) * ClusterSizeX + j] + - data[(i + 1) * ClusterSizeX + j + 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] = + data[i * ClusterSizeX + j] + + data[i * ClusterSizeX + j + 1] + + data[(i + 1) * ClusterSizeX + j] + + data[(i + 1) * ClusterSizeX + j + 1]; + } + + int index = std::max_element(sum_2x2_subcluster.begin(), + sum_2x2_subcluster.end()) - + sum_2x2_subcluster.begin(); + return std::make_pair(sum_2x2_subcluster[index], index); } - - int index = std::max_element(sum_2x2_subcluster.begin(), - sum_2x2_subcluster.end()) - - sum_2x2_subcluster.begin(); - return std::make_pair(sum_2x2_subcluster[index], index); - } -}; - -// Specialization for 2x2 clusters (only one sum exists) -template struct Cluster { - int16_t x; - int16_t y; - T data[4]; - - static constexpr uint8_t cluster_size_x = 2; - static constexpr uint8_t cluster_size_y = 2; - using value_type = T; - using coord_type = int16_t; - - T sum() const { return std::accumulate(data, data + 4, 0); } - - std::pair max_sum_2x2() const { - return std::make_pair(data[0] + data[1] + data[2] + data[3], - 0); // Only one possible 2x2 sum - } -}; - -// Specialization for 3x3 clusters -template struct Cluster { - int16_t x; - int16_t y; - T data[9]; - static constexpr uint8_t cluster_size_x = 3; - static constexpr uint8_t cluster_size_y = 3; - using value_type = T; - using coord_type = int16_t; - - T sum() const { return std::accumulate(data, data + 9, 0); } - - std::pair max_sum_2x2() const { - std::array sum_2x2_subclusters; - sum_2x2_subclusters[0] = data[0] + data[1] + data[3] + data[4]; - sum_2x2_subclusters[1] = data[1] + data[2] + data[4] + data[5]; - sum_2x2_subclusters[2] = data[3] + data[4] + data[6] + data[7]; - sum_2x2_subclusters[3] = data[4] + data[5] + data[7] + data[8]; - int index = std::max_element(sum_2x2_subclusters.begin(), - sum_2x2_subclusters.end()) - - sum_2x2_subclusters.begin(); - return std::make_pair(sum_2x2_subclusters[index], index); } }; diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index b3538eb..ea11162 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -77,7 +77,6 @@ class ClusterFinder { int has_center_pixel_y = ClusterSizeY % 2; m_clusters.set_frame_number(frame_number); - std::vector cluster_data(ClusterSizeX * ClusterSizeY); for (int iy = 0; iy < frame.shape(0); iy++) { for (int ix = 0; ix < frame.shape(1); ix++) { @@ -124,8 +123,9 @@ class ClusterFinder { // Store cluster if (value == max) { - // Zero out the cluster data - std::fill(cluster_data.begin(), cluster_data.end(), 0); + ClusterType cluster{}; + cluster.x = ix; + cluster.y = iy; // Fill the cluster data since we have a photon to store // It's worth redoing the look since most of the time we @@ -139,20 +139,15 @@ class ClusterFinder { static_cast(frame(iy + ir, ix + ic)) - static_cast( m_pedestal.mean(iy + ir, ix + ic)); - cluster_data[i] = + cluster.data[i] = tmp; // Watch for out of bounds access i++; } } } - ClusterType new_cluster{}; - new_cluster.x = ix; - new_cluster.y = iy; - std::copy(cluster_data.begin(), cluster_data.end(), - new_cluster.data); // Add the cluster to the output ClusterVector - m_clusters.push_back(new_cluster); + m_clusters.push_back(cluster); } } } diff --git a/include/aare/GainMap.hpp b/include/aare/GainMap.hpp index 41acb33..5311916 100644 --- a/include/aare/GainMap.hpp +++ b/include/aare/GainMap.hpp @@ -44,9 +44,8 @@ class GainMap { cl.data[j] = cl.data[j] * static_cast(m_gain_map(y, x)); } } else { - memset(cl.data, 0, - ClusterSizeX * ClusterSizeY * - sizeof(T)); // clear edge clusters + // clear edge clusters + cl.data.fill(0); } } } diff --git a/python/src/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp index f7fa796..ea02487 100644 --- a/python/src/bind_ClusterVector.hpp +++ b/python/src/bind_ClusterVector.hpp @@ -21,16 +21,14 @@ using namespace aare; #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-parameter" - template void define_ClusterVector(py::module &m, const std::string &typestr) { - using ClusterType = - Cluster; + using ClusterType = Cluster; auto class_name = fmt::format("ClusterVector_{}", typestr); py::class_, void>>( + Cluster, void>>( m, class_name.c_str(), py::buffer_protocol()) @@ -41,10 +39,11 @@ void define_ClusterVector(py::module &m, const std::string &typestr) { self.push_back(cluster); }) - .def("sum", [](ClusterVector &self) { - auto *vec = new std::vector(self.sum()); - return return_vector(vec); - }) + .def("sum", + [](ClusterVector &self) { + auto *vec = new std::vector(self.sum()); + return return_vector(vec); + }) .def_property_readonly("size", &ClusterVector::size) .def("item_size", &ClusterVector::item_size) .def_property_readonly("fmt", @@ -72,32 +71,30 @@ void define_ClusterVector(py::module &m, const std::string &typestr) { ); }); - // Free functions using ClusterVector - m.def("hitmap", - [](std::array image_size, ClusterVector &cv) { - - // Create a numpy array to hold the hitmap - // The shape of the array is (image_size[0], image_size[1]) - // note that the python array is passed as [row, col] which - // is the opposite of the clusters [x,y] - py::array_t hitmap(image_size); - auto r = hitmap.mutable_unchecked<2>(); - - // Initialize hitmap to 0 - for (py::ssize_t i = 0; i < r.shape(0); i++) - for (py::ssize_t j = 0; j < r.shape(1); j++) - r(i, j) = 0; - + // Free functions using ClusterVector + m.def("hitmap", + [](std::array image_size, ClusterVector &cv) { + // Create a numpy array to hold the hitmap + // The shape of the array is (image_size[0], image_size[1]) + // note that the python array is passed as [row, col] which + // is the opposite of the clusters [x,y] + py::array_t hitmap(image_size); + auto r = hitmap.mutable_unchecked<2>(); - // Loop over the clusters and increment the hitmap - // Skip out of bound clusters - for (const auto& cluster : cv) { - auto x = cluster.x; - auto y = cluster.y; - if(x>( + py::class_>( m, class_name.c_str(), py::buffer_protocol()) .def(py::init([](uint8_t x, uint8_t y, py::array_t data) { py::buffer_info buf_info = data.request(); - Type *ptr = static_cast(buf_info.ptr); - Cluster cluster; + Cluster cluster; cluster.x = x; cluster.y = y; - std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY, - cluster.data); // Copy array contents + auto r = data.template unchecked<1>(); // no bounds checks + for (py::ssize_t i = 0; i < data.size(); ++i) { + cluster.data[i] = r(i); + } return cluster; })); @@ -64,9 +65,6 @@ void define_cluster(py::module &m, const std::string &typestr) { */ } - - - template void define_cluster_finder_mt_bindings(py::module &m, @@ -206,6 +204,5 @@ void define_cluster_finder_bindings(py::module &m, const std::string &typestr) { return; }, py::arg(), py::arg("frame_number") = 0); - } #pragma GCC diagnostic pop diff --git a/src/Cluster.test.cpp b/src/Cluster.test.cpp index 879a5e7..ba9cda1 100644 --- a/src/Cluster.test.cpp +++ b/src/Cluster.test.cpp @@ -14,19 +14,6 @@ using namespace aare; -TEST_CASE("Correct Instantiation of Cluster and ClusterVector", - "[.cluster][.instantiation]") { - - CHECK(is_valid_cluster); - CHECK(is_valid_cluster); - CHECK(not is_valid_cluster); - CHECK(not is_valid_cluster); - CHECK(not is_valid_cluster); - - CHECK(not is_cluster_v); - CHECK(is_cluster_v>); -} - TEST_CASE("Test sum of Cluster", "[.cluster]") { Cluster cluster{0, 0, {1, 2, 3, 4}}; diff --git a/src/ClusterFile.test.cpp b/src/ClusterFile.test.cpp index 024bed4..3f15332 100644 --- a/src/ClusterFile.test.cpp +++ b/src/ClusterFile.test.cpp @@ -311,19 +311,19 @@ TEST_CASE("Write cluster with potential padding", "[.files][.ClusterFile]") { CHECK(read_cluster_vector.at(0).x == clustervec.at(0).x); CHECK(read_cluster_vector.at(0).y == clustervec.at(0).y); - CHECK(std::equal(clustervec.at(0).data, clustervec.at(0).data + 9, - read_cluster_vector.at(0).data, [](double a, double b) { - return std::abs(a - b) < - std::numeric_limits::epsilon(); - })); + CHECK(std::equal( + clustervec.at(0).data.begin(), clustervec.at(0).data.end(), + read_cluster_vector.at(0).data.begin(), [](double a, double b) { + return std::abs(a - b) < std::numeric_limits::epsilon(); + })); CHECK(read_cluster_vector.at(1).x == clustervec.at(1).x); CHECK(read_cluster_vector.at(1).y == clustervec.at(1).y); - CHECK(std::equal(clustervec.at(1).data, std::end(clustervec.at(1).data), - read_cluster_vector.at(1).data, [](double a, double b) { - return std::abs(a - b) < - std::numeric_limits::epsilon(); - })); + CHECK(std::equal( + clustervec.at(1).data.begin(), clustervec.at(1).data.end(), + read_cluster_vector.at(1).data.begin(), [](double a, double b) { + return std::abs(a - b) < std::numeric_limits::epsilon(); + })); } TEST_CASE("Read frame and modify cluster data", "[.files][.ClusterFile]") {