diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index fe27c6a..c3089b7 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -107,33 +107,51 @@ class ClusterFinder { */ template void process_pixel(const NDView &frame, - const int iy, const int ix) { + const int iy, const int ix) { constexpr int dy = ClusterSizeY / 2; constexpr int dx = ClusterSizeX / 2; constexpr int has_center_pixel_x = ClusterSizeX % 2; constexpr int has_center_pixel_y = ClusterSizeY % 2; - PEDESTAL_TYPE max = std::numeric_limits::min(); + PEDESTAL_TYPE max = std::numeric_limits::lowest(); PEDESTAL_TYPE total = 0; - const PEDESTAL_TYPE threshold = m_threshold(iy, ix); - const PEDESTAL_TYPE value = m_pd_corrected_frame(iy, ix); + const int cols = static_cast(frame.shape(1)); + const int rows = static_cast(frame.shape(0)); + const auto center = (static_cast(iy) * + static_cast(cols)) + + static_cast(ix); + const auto *corrected = m_pd_corrected_frame.data(); + const PEDESTAL_TYPE threshold = m_threshold.data()[center]; + const PEDESTAL_TYPE value = corrected[center]; if (value < -threshold) return; // NEGATIVE_PEDESTAL, nothing to do for this pixel // TODO! No pedestal update??? - for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { - for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) { - if constexpr (CheckBounds) { - if (ix + ic < 0 || ix + ic >= frame.shape(1) || - iy + ir < 0 || iy + ir >= frame.shape(0)) + if constexpr (CheckBounds) { + for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { + for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) { + const int x = ix + ic; + const int y = iy + ir; + if (x < 0 || x >= cols || y < 0 || y >= rows) continue; + const PEDESTAL_TYPE val = + corrected[(static_cast(y) * cols) + x]; + total += val; + max = std::max(max, val); + } + } + } else { + for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { + const auto *pixel = + corrected + static_cast(iy + ir) * cols + + (ix - dx); + for (int k = 0; k < ClusterSizeX; k++) { + const PEDESTAL_TYPE val = pixel[k]; + total += val; + max = std::max(max, val); } - const PEDESTAL_TYPE val = - m_pd_corrected_frame(iy + ir, ix + ic); - total += val; - max = std::max(max, val); } } @@ -144,7 +162,7 @@ class ClusterFinder { // pass, store the cluster below } else { // m_pedestal.push(iy, ix, frame(iy, ix)); // Safe option - m_pedestal.push(iy, ix, frame(iy, ix)); + m_pedestal.push(iy, ix, frame.data()[center]); return; // It was a pedestal value nothing to store } @@ -154,31 +172,50 @@ class ClusterFinder { cluster.x = ix; cluster.y = iy; - int i = 0; - for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { - for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) { - bool in_bounds = true; - if constexpr (CheckBounds) { - in_bounds = ix + ic >= 0 && ix + ic < frame.shape(1) && - iy + ir >= 0 && iy + ir < frame.shape(0); + if constexpr (CheckBounds) { + int i = 0; + for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { + for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) { + const int x = ix + ic; + const int y = iy + ir; + if (x >= 0 && x < cols && y >= 0 && y < rows) { + const PEDESTAL_TYPE corrected_value = + corrected[(static_cast(y) * cols) + + x]; + if constexpr ( + std::is_integral_v && + std::is_floating_point_v) { + cluster.data[i] = static_cast( + std::lround(corrected_value)); + } else { + cluster.data[i] = + static_cast(corrected_value); + } + } + i++; } - if (in_bounds) { + } + } else { + int i = 0; + for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { + const auto *pixel = + corrected + static_cast(iy + ir) * cols + + (ix - dx); + for (int k = 0; k < ClusterSizeX; k++, i++) { // If the cluster type is an integral type, and the // pedestal is a floating point type then we need to // round the value before storing it if constexpr (std::is_integral_v && std::is_floating_point_v) { - cluster.data[i] = static_cast(std::lround( - m_pd_corrected_frame(iy + ir, ix + ic))); + cluster.data[i] = + static_cast(std::lround(pixel[k])); } // On the other hand if both are floating point or both // are integral then we can just static cast directly else { - cluster.data[i] = static_cast( - m_pd_corrected_frame(iy + ir, ix + ic)); + cluster.data[i] = static_cast(pixel[k]); } } - i++; } } @@ -204,11 +241,23 @@ class ClusterFinder { m_clusters.set_frame_number(frame_number); - m_pd_corrected_frame = frame - m_pedestal.view(); + const int rows = static_cast(frame.shape(0)); const int cols = static_cast(frame.shape(1)); + // TODO! See if we can get the same performace using the operator- + // m_pd_corrected_frame = frame - m_pedestal.view(); + + //here we should be able to safely assume that the frame and corrected frame have the same size + auto n_pixels = frame.size(); + auto pd = m_pedestal.view().data(); + auto corrected = m_pd_corrected_frame.data(); + auto frame_data = frame.data(); + for (size_t i = 0; i < n_pixels; i++) { + corrected[i] = static_cast(frame_data[i]) - pd[i]; + } + // Interior pixels can skip the per-neighbour bounds checks; pixels // within dx/dy of an edge take the bounds-checked path. Iteration order // (row-major, increasing ix) is preserved so results are identical. diff --git a/python/aare/__init__.py b/python/aare/__init__.py index 36e1220..92c3652 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -8,8 +8,10 @@ from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile from ._aare import ( FastPedestal_d, FastPedestal_f, + FastPedestal_i16, Pedestal_d, Pedestal_f, + Pedestal_i16, ClusterFinder_Cluster3x3i, VarClusterFinder, ) diff --git a/python/src/module.cpp b/python/src/module.cpp index 26f9fc6..b668ba8 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -32,6 +32,7 @@ #include "var_cluster.hpp" // Pybind stuff +#include #include #include @@ -73,8 +74,10 @@ PYBIND11_MODULE(_aare, m) { define_pedestal_tracking_pixel_histogram_bindings(m); define_pedestal_bindings(m, "Pedestal_d"); define_pedestal_bindings(m, "Pedestal_f"); + define_pedestal_bindings(m, "Pedestal_i16"); define_fast_pedestal_bindings(m, "FastPedestal_d"); define_fast_pedestal_bindings(m, "FastPedestal_f"); + define_fast_pedestal_bindings(m, "FastPedestal_i16"); define_fit_bindings(m); define_interpolation_bindings(m); define_jungfrau_data_file_io_bindings(m); @@ -106,6 +109,7 @@ PYBIND11_MODULE(_aare, m) { DEFINE_BINDINGS_CLUSTERFINDER(int, 3, 3, uint16_t, i); DEFINE_BINDINGS_CLUSTERFINDER(double, 3, 3, uint16_t, d); DEFINE_BINDINGS_CLUSTERFINDER(float, 3, 3, uint16_t, f); + DEFINE_BINDINGS_CLUSTERFINDER(int16_t, 3, 3, uint16_t, i16); DEFINE_BINDINGS_CLUSTERFINDER(int, 5, 5, uint16_t, i); DEFINE_BINDINGS_CLUSTERFINDER(double, 5, 5, uint16_t, d); diff --git a/python/src/module_config.hpp b/python/src/module_config.hpp index 369140c..f166de5 100644 --- a/python/src/module_config.hpp +++ b/python/src/module_config.hpp @@ -1,4 +1,5 @@ #pragma once +#include //Configure module wide pedestal type for cluster finding using pd_type = double; \ No newline at end of file