diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d81e92..bedb54e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -388,6 +388,7 @@ set(PUBLICHEADERS include/aare/FileInterface.hpp include/aare/FilePtr.hpp include/aare/Frame.hpp + include/aare/PixelHistogram.hpp include/aare/GainMap.hpp include/aare/ROIGeometry.hpp include/aare/DetectorGeometry.hpp @@ -412,6 +413,7 @@ set(SourceFiles ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ROIGeometry.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/PixelHistogram.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/FilePtr.cpp diff --git a/include/aare/PixelHistogram.hpp b/include/aare/PixelHistogram.hpp new file mode 100644 index 0000000..e6a9ea6 --- /dev/null +++ b/include/aare/PixelHistogram.hpp @@ -0,0 +1,28 @@ +#pragma once +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" + +//Lets see if we need to hide it behind a pimpl +#include +#include +namespace bh = boost::histogram; + +namespace aare { +class PixelHistogram{ + using Axes = std::tuple< + bh::axis::regular, + bh::axis::integer, + bh::axis::integer + >; + using Hist = bh::histogram>; + Hist hist_; + + public: + PixelHistogram(int rows, int cols, int n_bins, double xmin, double xmax); + void fill(const NDView& image); + NDArray hdata() const; + NDArray bin_centers() const; + NDArray bin_edges() const; +}; + +} // namespace aare diff --git a/python/src/bind_PixelHistogram.hpp b/python/src/bind_PixelHistogram.hpp new file mode 100644 index 0000000..e3e85aa --- /dev/null +++ b/python/src/bind_PixelHistogram.hpp @@ -0,0 +1,78 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/PixelHistogram.hpp" +#include "np_helper.hpp" + +#include +#include +#include +#include + +namespace py = pybind11; +using namespace ::aare; + +void define_pixel_histogram_bindings(py::module &m) { + py::class_(m, "PixelHistogram", + "A histogram for pixel-wise statistics") + .def(py::init(), + R"( + Initialize a PixelHistogram. + + Args: + rows: Number of rows in the detector + cols: Number of columns in the detector + n_bins: Number of histogram bins + xmin: Minimum value for histogram range + xmax: Maximum value for histogram range + )", + py::arg("rows"), py::arg("cols"), py::arg("n_bins"), + py::arg("xmin"), py::arg("xmax")) + + .def("fill", + [](PixelHistogram &self, + py::array_t image) { + auto view = make_view_2d(image); + self.fill(view); + }, + R"( + Fill the histogram with image data. + + Args: + image: A 2D numpy array of pixel values (dtype: float64) + )", + py::arg("image")) + + .def("hdata", + [](const PixelHistogram &self) { + auto ptr = new NDArray(self.hdata()); + return return_image_data(ptr); + }, + R"( + Get the histogram data as a numpy array. + + Returns: + A 3D numpy array containing the histogram bins for each pixel + )") + + .def("bin_centers", + [](const PixelHistogram &self) { + auto ptr = new NDArray(self.bin_centers()); + return return_image_data(ptr); + }, + R"( + Get the bin centers along the value axis. + + Returns: + A 1D numpy array containing the center values for each histogram bin + )") + .def("bin_edges", + [](const PixelHistogram &self) { + auto ptr = new NDArray(self.bin_edges()); + return return_image_data(ptr); + }, + R"( + Get the bin edges along the value axis. + + Returns: + A 1D numpy array containing the edge values for the histogram bins + )"); +} diff --git a/python/src/module.cpp b/python/src/module.cpp index 19bbe29..02f4e6a 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -12,6 +12,7 @@ #include "bind_Defs.hpp" #include "bind_Eta.hpp" #include "bind_Interpolator.hpp" +#include "bind_PixelHistogram.hpp" #include "bind_PixelMap.hpp" #include "bind_RawFile.hpp" #include "bind_calibration.hpp" @@ -64,6 +65,7 @@ PYBIND11_MODULE(_aare, m) { define_raw_master_file_bindings(m); define_var_cluster_finder_bindings(m); define_pixel_map_bindings(m); + define_pixel_histogram_bindings(m); define_pedestal_bindings(m, "Pedestal_d"); define_pedestal_bindings(m, "Pedestal_f"); define_fit_bindings(m); diff --git a/src/PixelHistogram.cpp b/src/PixelHistogram.cpp new file mode 100644 index 0000000..bbc3de9 --- /dev/null +++ b/src/PixelHistogram.cpp @@ -0,0 +1,70 @@ +#include "aare/PixelHistogram.hpp" + +#include +#include +#include +#include + +namespace aare { + + +PixelHistogram::PixelHistogram(int rows, int cols, int n_bins, double xmin, double xmax): hist_( + bh::axis::regular(n_bins, xmin, xmax, "value"), + bh::axis::integer(0, cols, "y"), + bh::axis::integer(0, rows, "x") + ) { + // Constructor implementation +} + +NDArray PixelHistogram::hdata() const { + const auto bins = static_cast(hist_.axis(0).size()); + const auto cols = static_cast(hist_.axis(1).size()); + const auto rows = static_cast(hist_.axis(2).size()); + + NDArray data({rows, cols, bins}); + + const auto &storage = bh::unsafe_access::storage(hist_); + std::memcpy(data.data(), storage.data(), data.total_bytes()); + + return data; +} + +void PixelHistogram::fill(const NDView &image) { + for (ssize_t row = 0; row < image.shape(0); ++row) { + for (ssize_t col = 0; col < image.shape(1); ++col) { + hist_(image(row, col), col, row); + } + } +} + +NDArray PixelHistogram::bin_centers() const { + const auto& value_axis = hist_.axis(0); + const auto n_bins = static_cast(value_axis.size()); + + NDArray centers({n_bins}); + + for (ssize_t i = 0; i < n_bins; ++i) { + // Get the left and right edges of the bin and compute the center + double left = value_axis.value(i); + double right = value_axis.value(i + 1); + centers(i) = (left + right) / 2.0; + } + + return centers; +} + +NDArray PixelHistogram::bin_edges() const { + const auto& value_axis = hist_.axis(0); + const auto n_bins = static_cast(value_axis.size()); + + NDArray edges({n_bins + 1}); + + for (ssize_t i = 0; i <= n_bins; ++i) { + edges(i) = value_axis.value(i); + } + + return edges; +} + +} // namespace aare