From 5d8ad27b21af0d6b22390523263b2d48c442eb9d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Thu, 20 Mar 2025 12:52:04 +0100 Subject: [PATCH] Developer (#138) - Fully functioning variable size cluster finder - Added interpolation - Bit reordering for ADC SAR 05 --------- Co-authored-by: Patrick Co-authored-by: JulianHeymes Co-authored-by: Dhanya Thattil Co-authored-by: xiangyu.xie --- .clang-tidy | 42 +++++++++ .gitea/workflows/cmake_build.yml | 58 ++++++++++++ CMakeLists.txt | 5 ++ conda-recipe/meta.yaml | 3 +- include/aare/ClusterFile.hpp | 10 +++ include/aare/ClusterVector.hpp | 4 + include/aare/Interpolator.hpp | 29 ++++++ include/aare/NDArray.hpp | 15 ++-- include/aare/RawSubFile.hpp | 2 +- include/aare/VarClusterFinder.hpp | 2 +- include/aare/algorithm.hpp | 55 ++++++++++++ include/aare/defs.hpp | 2 + pyproject.toml | 3 +- python/aare/__init__.py | 3 +- python/examples/play.py | 89 +++++++++++------- python/src/cluster.hpp | 10 ++- python/src/cluster_file.hpp | 5 ++ python/src/ctb_raw_file.hpp | 4 +- python/src/file.hpp | 2 + python/src/interpolation.hpp | 58 ++++++++++++ python/src/module.cpp | 2 + python/src/np_helper.hpp | 12 +-- python/src/var_cluster.hpp | 37 +++++++- src/ClusterFile.cpp | 109 ++++++++++++++++++++-- src/Dtype.cpp | 2 +- src/File.cpp | 2 +- src/Interpolator.cpp | 144 ++++++++++++++++++++++++++++++ src/NDArray.test.cpp | 19 ++++ src/RawFile.cpp | 3 +- src/algorithm.test.cpp | 73 +++++++++++++++ 30 files changed, 743 insertions(+), 61 deletions(-) create mode 100644 .clang-tidy create mode 100644 .gitea/workflows/cmake_build.yml create mode 100644 include/aare/Interpolator.hpp create mode 100644 include/aare/algorithm.hpp create mode 100644 python/src/interpolation.hpp create mode 100644 src/Interpolator.cpp create mode 100644 src/algorithm.test.cpp diff --git a/.clang-tidy b/.clang-tidy new file mode 100644 index 0000000..a2ab6c1 --- /dev/null +++ b/.clang-tidy @@ -0,0 +1,42 @@ + +--- +Checks: '*, + -altera-*, + -android-cloexec-fopen, + -cppcoreguidelines-pro-bounds-array-to-pointer-decay, + -cppcoreguidelines-pro-bounds-pointer-arithmetic, + -fuchsia*, + -readability-else-after-return, + -readability-avoid-const-params-in-decls, + -readability-identifier-length, + -cppcoreguidelines-pro-bounds-constant-array-index, + -cppcoreguidelines-pro-type-reinterpret-cast, + -llvm-header-guard, + -modernize-use-nodiscard, + -misc-non-private-member-variables-in-classes, + -readability-static-accessed-through-instance, + -readability-braces-around-statements, + -readability-isolate-declaration, + -readability-implicit-bool-conversion, + -readability-identifier-length, + -readability-identifier-naming, + -hicpp-signed-bitwise, + -hicpp-no-array-decay, + -hicpp-braces-around-statements, + -google-runtime-references, + -google-readability-todo, + -google-readability-braces-around-statements, + -modernize-use-trailing-return-type, + -llvmlibc-*' + +HeaderFilterRegex: \.hpp +FormatStyle: none +CheckOptions: + - { key: readability-identifier-naming.NamespaceCase, value: lower_case } + # - { key: readability-identifier-naming.FunctionCase, value: lower_case } + - { key: readability-identifier-naming.ClassCase, value: CamelCase } + # - { key: readability-identifier-naming.MethodCase, value: CamelCase } + # - { key: readability-identifier-naming.StructCase, value: CamelCase } + # - { key: readability-identifier-naming.VariableCase, value: lower_case } + - { key: readability-identifier-naming.GlobalConstantCase, value: UPPER_CASE } +... diff --git a/.gitea/workflows/cmake_build.yml b/.gitea/workflows/cmake_build.yml new file mode 100644 index 0000000..43a0181 --- /dev/null +++ b/.gitea/workflows/cmake_build.yml @@ -0,0 +1,58 @@ +name: Build the package using cmake then documentation + +on: + workflow_dispatch: + push: + + + +permissions: + contents: read + pages: write + id-token: write + +jobs: + build: + strategy: + fail-fast: false + matrix: + platform: [ubuntu-latest, ] # macos-12, windows-2019] + python-version: ["3.12",] + + runs-on: ${{ matrix.platform }} + + # The setup-miniconda action needs this to activate miniconda + defaults: + run: + shell: "bash -l {0}" + + steps: + - uses: actions/checkout@v4 + + - name: Setup dev env + run: | + sudo apt-get update + sudo apt-get -y install cmake gcc g++ + + - name: Get conda + uses: conda-incubator/setup-miniconda@v3.0.4 + with: + python-version: ${{ matrix.python-version }} + channels: conda-forge + + - name: Prepare + run: conda install doxygen sphinx=7.1.2 breathe pybind11 sphinx_rtd_theme furo nlohmann_json zeromq fmt numpy + + - name: Build library + run: | + mkdir build + cd build + cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_DOCS=ON + make -j 2 + make docs + + + + + + diff --git a/CMakeLists.txt b/CMakeLists.txt index b93b513..4772f0b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -60,6 +60,8 @@ if(AARE_SYSTEM_LIBRARIES) set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE) set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE) set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE) + # Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available + # on conda-forge endif() if(AARE_VERBOSE) @@ -78,6 +80,7 @@ endif() set(CMAKE_EXPORT_COMPILE_COMMANDS ON) if(AARE_FETCH_LMFIT) + #TODO! Should we fetch lmfit from the web or inlcude a tar.gz in the repo? set(lmfit_patch git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/lmfit.patch) FetchContent_Declare( lmfit @@ -343,6 +346,7 @@ set(SourceFiles ${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp @@ -382,6 +386,7 @@ endif() if(AARE_TESTS) set(TestSources + ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index ffa95a7..120854b 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -1,6 +1,7 @@ package: name: aare - version: 2025.2.18 #TODO! how to not duplicate this? + version: 2025.3.18 #TODO! how to not duplicate this? + diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index b796763..5bea342 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -8,11 +8,17 @@ namespace aare { +//TODO! Template this? struct Cluster3x3 { int16_t x; int16_t y; int32_t data[9]; }; +struct Cluster2x2 { + int16_t x; + int16_t y; + int32_t data[4]; +}; typedef enum { cBottomLeft = 0, @@ -37,6 +43,7 @@ struct Eta2 { double x; double y; corner c; + int32_t sum; }; struct ClusterAnalysis { @@ -97,6 +104,8 @@ class ClusterFile { */ ClusterVector read_clusters(size_t n_clusters); + ClusterVector read_clusters(size_t n_clusters, ROI roi); + /** * @brief Read a single frame from the file and return the clusters. The * cluster vector will have the frame number set. @@ -131,5 +140,6 @@ int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad, NDArray calculate_eta2(ClusterVector &clusters); Eta2 calculate_eta2(Cluster3x3 &cl); +Eta2 calculate_eta2(Cluster2x2 &cl); } // namespace aare diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index febf06c..1c15a22 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -231,6 +231,10 @@ template class ClusterVector { return *reinterpret_cast(element_ptr(i)); } + template const V &at(size_t i) const { + return *reinterpret_cast(element_ptr(i)); + } + const std::string_view fmt_base() const { // TODO! how do we match on coord_t? return m_fmt_base; diff --git a/include/aare/Interpolator.hpp b/include/aare/Interpolator.hpp new file mode 100644 index 0000000..4905bce --- /dev/null +++ b/include/aare/Interpolator.hpp @@ -0,0 +1,29 @@ +#pragma once +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include "aare/ClusterVector.hpp" +#include "aare/ClusterFile.hpp" //Cluster_3x3 +namespace aare{ + +struct Photon{ + double x; + double y; + double energy; +}; + +class Interpolator{ + NDArray m_ietax; + NDArray m_ietay; + + NDArray m_etabinsx; + NDArray m_etabinsy; + NDArray m_energy_bins; + public: + Interpolator(NDView etacube, NDView xbins, NDView ybins, NDView ebins); + NDArray get_ietax(){return m_ietax;} + NDArray get_ietay(){return m_ietay;} + + std::vector interpolate(const ClusterVector& clusters); +}; + +} // namespace aare \ No newline at end of file diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index cfa5b5c..45d3a83 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -102,6 +102,9 @@ class NDArray : public ArrayExpr, Ndim> { auto begin() { return data_; } auto end() { return data_ + size_; } + auto begin() const { return data_; } + auto end() const { return data_ + size_; } + using value_type = T; NDArray &operator=(NDArray &&other) noexcept; // Move assign @@ -388,12 +391,12 @@ NDArray NDArray::operator*(const T &value) { result *= value; return result; } -template void NDArray::Print() { - if (shape_[0] < 20 && shape_[1] < 20) - Print_all(); - else - Print_some(); -} +// template void NDArray::Print() { +// if (shape_[0] < 20 && shape_[1] < 20) +// Print_all(); +// else +// Print_some(); +// } template std::ostream &operator<<(std::ostream &os, const NDArray &arr) { diff --git a/include/aare/RawSubFile.hpp b/include/aare/RawSubFile.hpp index 89c278e..1d554e8 100644 --- a/include/aare/RawSubFile.hpp +++ b/include/aare/RawSubFile.hpp @@ -64,7 +64,7 @@ class RawSubFile { size_t bytes_per_frame() const { return m_bytes_per_frame; } size_t pixels_per_frame() const { return m_rows * m_cols; } - size_t bytes_per_pixel() const { return m_bitdepth / 8; } + size_t bytes_per_pixel() const { return m_bitdepth / bits_per_byte; } private: template diff --git a/include/aare/VarClusterFinder.hpp b/include/aare/VarClusterFinder.hpp index d4d51cc..ea62a9d 100644 --- a/include/aare/VarClusterFinder.hpp +++ b/include/aare/VarClusterFinder.hpp @@ -7,7 +7,7 @@ #include "aare/NDArray.hpp" -const int MAX_CLUSTER_SIZE = 200; +const int MAX_CLUSTER_SIZE = 50; namespace aare { template class VarClusterFinder { diff --git a/include/aare/algorithm.hpp b/include/aare/algorithm.hpp new file mode 100644 index 0000000..5d6dc57 --- /dev/null +++ b/include/aare/algorithm.hpp @@ -0,0 +1,55 @@ + +#pragma once +#include +#include +#include +#include + +namespace aare { +/** + * @brief Find the index of the last element smaller than val + * assume a sorted array + */ +template +size_t last_smaller(const T* first, const T* last, T val) { + for (auto iter = first+1; iter != last; ++iter) { + if (*iter > val) { + return std::distance(first, iter-1); + } + } + return std::distance(first, last-1); +} + +template +size_t last_smaller(const NDArray& arr, T val) { + return last_smaller(arr.begin(), arr.end(), val); +} + + +template +size_t nearest_index(const T* first, const T* last, T val) { + auto iter = std::min_element(first, last, + [val](T a, T b) { + return std::abs(a - val) < std::abs(b - val); + }); + return std::distance(first, iter); +} + +template +size_t nearest_index(const NDArray& arr, T val) { + return nearest_index(arr.begin(), arr.end(), val); +} + +template +size_t nearest_index(const std::vector& vec, T val) { + return nearest_index(vec.data(), vec.data()+vec.size(), val); +} + +template +size_t nearest_index(const std::array& arr, T val) { + return nearest_index(arr.data(), arr.data()+arr.size(), val); +} + + + +} // namespace aare \ No newline at end of file diff --git a/include/aare/defs.hpp b/include/aare/defs.hpp index db1a47b..4559882 100644 --- a/include/aare/defs.hpp +++ b/include/aare/defs.hpp @@ -38,6 +38,8 @@ namespace aare { +inline constexpr size_t bits_per_byte = 8; + void assert_failed(const std::string &msg); diff --git a/pyproject.toml b/pyproject.toml index 6dc941e..b9bf7d2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,8 @@ build-backend = "scikit_build_core.build" [project] name = "aare" -version = "2025.2.18" +version = "2025.3.18" + [tool.scikit-build] diff --git a/python/aare/__init__.py b/python/aare/__init__.py index f4c19cc..058d7cf 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -7,11 +7,12 @@ from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder from ._aare import DetectorType from ._aare import ClusterFile from ._aare import hitmap +from ._aare import ROI from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i from ._aare import fit_gaus, fit_pol1 - +from ._aare import Interpolator from .CtbRawFile import CtbRawFile from .RawFile import RawFile from .ScanParameters import ScanParameters diff --git a/python/examples/play.py b/python/examples/play.py index 37754df..da469dc 100644 --- a/python/examples/play.py +++ b/python/examples/play.py @@ -1,50 +1,79 @@ import sys sys.path.append('/home/l_msdetect/erik/aare/build') -#Our normal python imports -from pathlib import Path -import matplotlib.pyplot as plt +from aare._aare import ClusterVector_i, Interpolator + +import pickle import numpy as np +import matplotlib.pyplot as plt import boost_histogram as bh +import torch +import math import time -import aare -data = np.random.normal(10, 1, 1000) +def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2): + """ + Generate a 2D gaussian as position mx, my, with sigma=sigma. + The gaussian is placed on a 2x2 pixel matrix with resolution + res in one dimesion. + """ + x = torch.linspace(0, pixel_size*grid_size, res) + x,y = torch.meshgrid(x,x, indexing="ij") + return 1 / (2*math.pi*sigma**2) * \ + torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2))) -hist = bh.Histogram(bh.axis.Regular(10, 0, 20)) -hist.fill(data) +scale = 1000 #Scale factor when converting to integer +pixel_size = 25 #um +grid = 2 +resolution = 100 +sigma_um = 10 +xa = np.linspace(0,grid*pixel_size,resolution) +ticks = [0, 25, 50] + +hit = np.array((20,20)) +etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl" + +local_resolution = 99 +grid_size = 3 +xaxis = np.linspace(0,grid_size*pixel_size, local_resolution) +t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution) +pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1) +pixels = pixels.numpy() +pixels = (pixels*scale).astype(np.int32) +v = ClusterVector_i(3,3) +v.push_back(1,1, pixels) + +with open(etahist_fname, "rb") as f: + hist = pickle.load(f) +eta = hist.view().copy() +etabinsx = np.array(hist.axes.edges.T[0].flat) +etabinsy = np.array(hist.axes.edges.T[1].flat) +ebins = np.array(hist.axes.edges.T[2].flat) +p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1]) -x = hist.axes[0].centers -y = hist.values() -y_err = np.sqrt(y)+1 -res = aare.fit_gaus(x, y, y_err, chi2 = True) - -t_elapsed = time.perf_counter()-t0 -print(f'Histogram filling took: {t_elapsed:.3f}s {total_clusters/t_elapsed/1e6:.3f}M clusters/s') +#Generate the hit -histogram_data = hist3d.counts() -x = hist3d.axes[2].edges[:-1] -y = histogram_data[100,100,:] -xx = np.linspace(x[0], x[-1]) -# fig, ax = plt.subplots() -# ax.step(x, y, where = 'post') -y_err = np.sqrt(y) -y_err = np.zeros(y.size) -y_err += 1 -# par = fit_gaus2(y,x, y_err) -# ax.plot(xx, gaus(xx,par)) -# print(par) +tmp = p.interpolate(v) +print(f'tmp:{tmp}') +pos = np.array((tmp['x'], tmp['y']))*25 -res = fit_gaus(y,x) -res2 = fit_gaus(y,x, y_err) -print(res) -print(res2) +print(pixels) +fig, ax = plt.subplots(figsize = (7,7)) +ax.pcolormesh(xaxis, xaxis, t) +ax.plot(*pos, 'o') +ax.set_xticks([0,25,50,75]) +ax.set_yticks([0,25,50,75]) +ax.set_xlim(0,75) +ax.set_ylim(0,75) +ax.grid() +print(f'{hit=}') +print(f'{pos=}') \ No newline at end of file diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp index 792b7e6..3db816a 100644 --- a/python/src/cluster.hpp +++ b/python/src/cluster.hpp @@ -20,7 +20,13 @@ template void define_cluster_vector(py::module &m, const std::string &typestr) { auto class_name = fmt::format("ClusterVector_{}", typestr); py::class_>(m, class_name.c_str(), py::buffer_protocol()) - .def(py::init()) + .def(py::init(), + py::arg("cluster_size_x") = 3, py::arg("cluster_size_y") = 3) + .def("push_back", + [](ClusterVector &self, int x, int y, py::array_t data) { + // auto view = make_view_2d(data); + self.push_back(x, y, reinterpret_cast(data.data())); + }) .def_property_readonly("size", &ClusterVector::size) .def("item_size", &ClusterVector::item_size) .def_property_readonly("fmt", @@ -38,6 +44,8 @@ void define_cluster_vector(py::module &m, const std::string &typestr) { auto *vec = new std::vector(self.sum_2x2()); return return_vector(vec); }) + .def_property_readonly("cluster_size_x", &ClusterVector::cluster_size_x) + .def_property_readonly("cluster_size_y", &ClusterVector::cluster_size_y) .def_property_readonly("capacity", &ClusterVector::capacity) .def_property("frame_number", &ClusterVector::frame_number, &ClusterVector::set_frame_number) diff --git a/python/src/cluster_file.hpp b/python/src/cluster_file.hpp index 8a431b5..f587443 100644 --- a/python/src/cluster_file.hpp +++ b/python/src/cluster_file.hpp @@ -31,6 +31,11 @@ void define_cluster_file_io_bindings(py::module &m) { auto v = new ClusterVector(self.read_clusters(n_clusters)); return v; },py::return_value_policy::take_ownership) + .def("read_clusters", + [](ClusterFile &self, size_t n_clusters, ROI roi) { + auto v = new ClusterVector(self.read_clusters(n_clusters, roi)); + return v; + },py::return_value_policy::take_ownership) .def("read_frame", [](ClusterFile &self) { auto v = new ClusterVector(self.read_frame()); diff --git a/python/src/ctb_raw_file.hpp b/python/src/ctb_raw_file.hpp index 9ce656d..56e571b 100644 --- a/python/src/ctb_raw_file.hpp +++ b/python/src/ctb_raw_file.hpp @@ -32,7 +32,7 @@ m.def("adc_sar_05_decode64to16", [](py::array_t input) { } //Create a 2D output array with the same shape as the input - std::vector shape{input.shape(0), input.shape(1)/8}; + std::vector shape{input.shape(0), input.shape(1)/static_cast(bits_per_byte)}; py::array_t output(shape); //Create a view of the input and output arrays @@ -53,7 +53,7 @@ m.def("adc_sar_04_decode64to16", [](py::array_t input) { } //Create a 2D output array with the same shape as the input - std::vector shape{input.shape(0), input.shape(1)/8}; + std::vector shape{input.shape(0), input.shape(1)/static_cast(bits_per_byte)}; py::array_t output(shape); //Create a view of the input and output arrays diff --git a/python/src/file.hpp b/python/src/file.hpp index c3c800c..0d64e16 100644 --- a/python/src/file.hpp +++ b/python/src/file.hpp @@ -195,6 +195,8 @@ void define_file_io_bindings(py::module &m) { py::class_(m, "ROI") .def(py::init<>()) + .def(py::init(), py::arg("xmin"), + py::arg("xmax"), py::arg("ymin"), py::arg("ymax")) .def_readwrite("xmin", &ROI::xmin) .def_readwrite("xmax", &ROI::xmax) .def_readwrite("ymin", &ROI::ymin) diff --git a/python/src/interpolation.hpp b/python/src/interpolation.hpp new file mode 100644 index 0000000..02742e1 --- /dev/null +++ b/python/src/interpolation.hpp @@ -0,0 +1,58 @@ +#include "aare/Interpolator.hpp" +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include "np_helper.hpp" +#include +#include +#include +#include + +namespace py = pybind11; +void define_interpolation_bindings(py::module &m) { + + PYBIND11_NUMPY_DTYPE(aare::Photon, x,y,energy); + + py::class_(m, "Interpolator") + .def(py::init([](py::array_t etacube, py::array_t xbins, + py::array_t ybins, py::array_t ebins) { + return Interpolator(make_view_3d(etacube), make_view_1d(xbins), + make_view_1d(ybins), make_view_1d(ebins)); + })) + .def("get_ietax", [](Interpolator& self){ + auto*ptr = new NDArray{}; + *ptr = self.get_ietax(); + return return_image_data(ptr); + }) + .def("get_ietay", [](Interpolator& self){ + auto*ptr = new NDArray{}; + *ptr = self.get_ietay(); + return return_image_data(ptr); + }) + .def("interpolate", [](Interpolator& self, const ClusterVector& clusters){ + auto photons = self.interpolate(clusters); + auto* ptr = new std::vector{photons}; + return return_vector(ptr); + }); + + // TODO! Evaluate without converting to double + m.def( + "hej", + []() { + // auto boost_histogram = py::module_::import("boost_histogram"); + // py::object axis = + // boost_histogram.attr("axis").attr("Regular")(10, 0.0, 10.0); + // py::object histogram = boost_histogram.attr("Histogram")(axis); + // return histogram; + // return h; + }, + R"( + Evaluate a 1D Gaussian function for all points in x using parameters par. + + Parameters + ---------- + x : array_like + The points at which to evaluate the Gaussian function. + par : array_like + The parameters of the Gaussian function. The first element is the amplitude, the second element is the mean, and the third element is the standard deviation. + )"); +} \ No newline at end of file diff --git a/python/src/module.cpp b/python/src/module.cpp index 70d143f..43f48ba 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -9,6 +9,7 @@ #include "cluster.hpp" #include "cluster_file.hpp" #include "fit.hpp" +#include "interpolation.hpp" //Pybind stuff #include @@ -31,5 +32,6 @@ PYBIND11_MODULE(_aare, m) { define_cluster_collector_bindings(m); define_cluster_file_sink_bindings(m); define_fit_bindings(m); + define_interpolation_bindings(m); } \ No newline at end of file diff --git a/python/src/np_helper.hpp b/python/src/np_helper.hpp index 6e92830..1845196 100644 --- a/python/src/np_helper.hpp +++ b/python/src/np_helper.hpp @@ -40,25 +40,25 @@ template py::array return_vector(std::vector *vec) { } // todo rewrite generic -template auto get_shape_3d(py::array_t arr) { +template auto get_shape_3d(const py::array_t& arr) { return aare::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)}; } -template auto make_view_3d(py::array_t arr) { +template auto make_view_3d(py::array_t& arr) { return aare::NDView(arr.mutable_data(), get_shape_3d(arr)); } -template auto get_shape_2d(py::array_t arr) { +template auto get_shape_2d(const py::array_t& arr) { return aare::Shape<2>{arr.shape(0), arr.shape(1)}; } -template auto get_shape_1d(py::array_t arr) { +template auto get_shape_1d(const py::array_t& arr) { return aare::Shape<1>{arr.shape(0)}; } -template auto make_view_2d(py::array_t arr) { +template auto make_view_2d(py::array_t& arr) { return aare::NDView(arr.mutable_data(), get_shape_2d(arr)); } -template auto make_view_1d(py::array_t arr) { +template auto make_view_1d(py::array_t& arr) { return aare::NDView(arr.mutable_data(), get_shape_1d(arr)); } \ No newline at end of file diff --git a/python/src/var_cluster.hpp b/python/src/var_cluster.hpp index f3a5741..f7b373f 100644 --- a/python/src/var_cluster.hpp +++ b/python/src/var_cluster.hpp @@ -19,15 +19,24 @@ using namespace::aare; void define_var_cluster_finder_bindings(py::module &m) { PYBIND11_NUMPY_DTYPE(VarClusterFinder::Hit, size, row, col, - reserved, energy, max); + reserved, energy, max, rows, cols, enes); py::class_>(m, "VarClusterFinder") .def(py::init, double>()) .def("labeled", [](VarClusterFinder &self) { - auto ptr = new NDArray(self.labeled()); + auto *ptr = new NDArray(self.labeled()); return return_image_data(ptr); }) + .def("set_noiseMap", + [](VarClusterFinder &self, + py::array_t + noise_map) { + auto noise_map_span = make_view_2d(noise_map); + self.set_noiseMap(noise_map_span); + }) + .def("set_peripheralThresholdFactor", + &VarClusterFinder::set_peripheralThresholdFactor) .def("find_clusters", [](VarClusterFinder &self, py::array_t @@ -35,6 +44,30 @@ void define_var_cluster_finder_bindings(py::module &m) { auto view = make_view_2d(img); self.find_clusters(view); }) + .def("find_clusters_X", + [](VarClusterFinder &self, + py::array_t + img) { + auto img_span = make_view_2d(img); + self.find_clusters_X(img_span); + }) + .def("single_pass", + [](VarClusterFinder &self, + py::array_t + img) { + auto img_span = make_view_2d(img); + self.single_pass(img_span); + }) + .def("hits", + [](VarClusterFinder &self) { + auto ptr = new std::vector::Hit>( + self.steal_hits()); + return return_vector(ptr); + }) + .def("clear_hits", + [](VarClusterFinder &self) { + self.clear_hits(); + }) .def("steal_hits", [](VarClusterFinder &self) { auto ptr = new std::vector::Hit>( diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index 2928d26..2e23e09 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -108,6 +108,79 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters) { return clusters; } +ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } + + ClusterVector clusters(3,3); + clusters.reserve(n_clusters); + + int32_t iframe = 0; // frame number needs to be 4 bytes! + size_t nph_read = 0; + uint32_t nn = m_num_left; + uint32_t nph = m_num_left; // number of clusters in frame needs to be 4 + + // auto buf = reinterpret_cast(clusters.data()); + // auto buf = clusters.data(); + + Cluster3x3 tmp; //this would break if the cluster size changes + + // if there are photons left from previous frame read them first + if (nph) { + if (nph > n_clusters) { + // if we have more photons left in the frame then photons to read we + // read directly the requested number + nn = n_clusters; + } else { + nn = nph; + } + //Read one cluster, in the ROI push back + // nph_read += fread((buf + nph_read*clusters.item_size()), + // clusters.item_size(), nn, fp); + for(size_t i = 0; i < nn; i++){ + fread(&tmp, sizeof(tmp), 1, fp); + if(tmp.x >= roi.xmin && tmp.x <= roi.xmax && tmp.y >= roi.ymin && tmp.y <= roi.ymax){ + clusters.push_back(tmp.x, tmp.y, reinterpret_cast(tmp.data)); + nph_read++; + } + } + + m_num_left = nph - nn; // write back the number of photons left + } + + if (nph_read < n_clusters) { + // keep on reading frames and photons until reaching n_clusters + while (fread(&iframe, sizeof(iframe), 1, fp)) { + // read number of clusters in frame + if (fread(&nph, sizeof(nph), 1, fp)) { + if (nph > (n_clusters - nph_read)) + nn = n_clusters - nph_read; + else + nn = nph; + + // nph_read += fread((buf + nph_read*clusters.item_size()), + // clusters.item_size(), nn, fp); + for(size_t i = 0; i < nn; i++){ + fread(&tmp, sizeof(tmp), 1, fp); + if(tmp.x >= roi.xmin && tmp.x <= roi.xmax && tmp.y >= roi.ymin && tmp.y <= roi.ymax){ + clusters.push_back(tmp.x, tmp.y, reinterpret_cast(tmp.data)); + nph_read++; + } + } + m_num_left = nph - nn; + } + if (nph_read >= n_clusters) + break; + } + } + + // Resize the vector to the number of clusters. + // No new allocation, only change bounds. + clusters.resize(nph_read); + return clusters; +} + ClusterVector ClusterFile::read_frame() { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); @@ -268,11 +341,23 @@ ClusterVector ClusterFile::read_frame() { NDArray calculate_eta2(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; + + if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) { + 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; + } + }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 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; + } + }else{ + throw std::runtime_error("Only 3x3 and 2x2 clusters are supported"); } + return eta2; } @@ -290,7 +375,7 @@ Eta2 calculate_eta2(Cluster3x3 &cl) { 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) @@ -333,6 +418,20 @@ Eta2 calculate_eta2(Cluster3x3 &cl) { return eta; } + +Eta2 calculate_eta2(Cluster2x2 &cl) { + Eta2 eta{}; + if ((cl.data[0] + cl.data[1]) != 0) + eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); + if ((cl.data[0] + cl.data[2]) != 0) + 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; +} + + + int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y) { diff --git a/src/Dtype.cpp b/src/Dtype.cpp index 565d509..b818ea3 100644 --- a/src/Dtype.cpp +++ b/src/Dtype.cpp @@ -70,7 +70,7 @@ uint8_t Dtype::bitdepth() const { /** * @brief Get the number of bytes of the data type */ -size_t Dtype::bytes() const { return bitdepth() / 8; } +size_t Dtype::bytes() const { return bitdepth() / bits_per_byte; } /** * @brief Construct a DType object from a TypeIndex diff --git a/src/File.cpp b/src/File.cpp index 1180967..3c68eff 100644 --- a/src/File.cpp +++ b/src/File.cpp @@ -73,7 +73,7 @@ size_t File::tell() const { return file_impl->tell(); } size_t File::rows() const { return file_impl->rows(); } size_t File::cols() const { return file_impl->cols(); } size_t File::bitdepth() const { return file_impl->bitdepth(); } -size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / 8; } +size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / bits_per_byte; } DetectorType File::detector_type() const { return file_impl->detector_type(); } diff --git a/src/Interpolator.cpp b/src/Interpolator.cpp new file mode 100644 index 0000000..7f82533 --- /dev/null +++ b/src/Interpolator.cpp @@ -0,0 +1,144 @@ +#include "aare/Interpolator.hpp" +#include "aare/algorithm.hpp" + +namespace aare { + +Interpolator::Interpolator(NDView etacube, NDView xbins, + NDView ybins, NDView ebins) + : m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins), m_energy_bins(ebins) { + if (etacube.shape(0) != xbins.size() || etacube.shape(1) != ybins.size() || + etacube.shape(2) != ebins.size()) { + throw std::invalid_argument( + "The shape of the etacube does not match the shape of the bins"); + } + + // Cumulative sum in the x direction + for (ssize_t i = 1; i < m_ietax.shape(0); i++) { + for (ssize_t j = 0; j < m_ietax.shape(1); j++) { + for (ssize_t k = 0; k < m_ietax.shape(2); k++) { + m_ietax(i, j, k) += m_ietax(i - 1, j, k); + } + } + } + + // Normalize by the highest row, if norm less than 1 don't do anything + for (ssize_t i = 0; i < m_ietax.shape(0); i++) { + for (ssize_t j = 0; j < m_ietax.shape(1); j++) { + for (ssize_t k = 0; k < m_ietax.shape(2); k++) { + auto val = m_ietax(m_ietax.shape(0) - 1, j, k); + double norm = val < 1 ? 1 : val; + m_ietax(i, j, k) /= norm; + } + } + } + + // Cumulative sum in the y direction + for (ssize_t i = 0; i < m_ietay.shape(0); i++) { + for (ssize_t j = 1; j < m_ietay.shape(1); j++) { + for (ssize_t k = 0; k < m_ietay.shape(2); k++) { + m_ietay(i, j, k) += m_ietay(i, j - 1, k); + } + } + } + + // Normalize by the highest column, if norm less than 1 don't do anything + for (ssize_t i = 0; i < m_ietay.shape(0); i++) { + for (ssize_t j = 0; j < m_ietay.shape(1); j++) { + for (ssize_t k = 0; k < m_ietay.shape(2); k++) { + auto val = m_ietay(i, m_ietay.shape(1) - 1, k); + double norm = val < 1 ? 1 : val; + m_ietay(i, j, k) /= norm; + } + } + } +} + +std::vector Interpolator::interpolate(const ClusterVector& clusters) { + std::vector photons; + photons.reserve(clusters.size()); + + if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) { + for (size_t i = 0; i(i); + Eta2 eta= calculate_eta2(cluster); + + Photon photon; + photon.x = cluster.x; + photon.y = cluster.y; + photon.energy = eta.sum; + + // auto ie = nearest_index(m_energy_bins, photon.energy)-1; + // auto ix = nearest_index(m_etabinsx, eta.x)-1; + // auto iy = nearest_index(m_etabinsy, eta.y)-1; + //Finding the index of the last element that is smaller + //should work fine as long as we have many bins + auto ie = last_smaller(m_energy_bins, photon.energy); + auto ix = last_smaller(m_etabinsx, eta.x); + auto iy = last_smaller(m_etabinsy, eta.y); + + // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy); + + double dX, dY; + int ex, ey; + // cBottomLeft = 0, + // cBottomRight = 1, + // cTopLeft = 2, + // cTopRight = 3 + switch (eta.c) { + case cTopLeft: + dX = -1.; + dY = 0.; + break; + case cTopRight:; + dX = 0.; + dY = 0.; + break; + case cBottomLeft: + dX = -1.; + dY = -1.; + break; + case cBottomRight: + dX = 0.; + dY = -1.; + break; + } + photon.x += m_ietax(ix, iy, ie)*2 + dX; + photon.y += m_ietay(ix, iy, ie)*2 + dY; + photons.push_back(photon); + } + }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){ + for (size_t i = 0; i(i); + Eta2 eta= calculate_eta2(cluster); + + Photon photon; + photon.x = cluster.x; + photon.y = cluster.y; + photon.energy = eta.sum; + + //Now do some actual interpolation. + //Find which energy bin the cluster is in + // auto ie = nearest_index(m_energy_bins, photon.energy)-1; + // auto ix = nearest_index(m_etabinsx, eta.x)-1; + // auto iy = nearest_index(m_etabinsy, eta.y)-1; + //Finding the index of the last element that is smaller + //should work fine as long as we have many bins + auto ie = last_smaller(m_energy_bins, photon.energy); + auto ix = last_smaller(m_etabinsx, eta.x); + auto iy = last_smaller(m_etabinsy, eta.y); + + photon.x += m_ietax(ix, iy, ie)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2 + photon.y += m_ietay(ix, iy, ie)*2; + photons.push_back(photon); + } + + }else{ + throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation"); + } + + + return photons; +} + +} // namespace aare \ No newline at end of file diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp index 942481c..eff3e2c 100644 --- a/src/NDArray.test.cpp +++ b/src/NDArray.test.cpp @@ -2,6 +2,7 @@ #include #include #include +#include using aare::NDArray; using aare::NDView; @@ -34,6 +35,24 @@ TEST_CASE("Construct from an NDView") { } } +TEST_CASE("3D NDArray from NDView"){ + std::vector data(27); + std::iota(data.begin(), data.end(), 0); + NDView view(data.data(), Shape<3>{3, 3, 3}); + NDArray image(view); + REQUIRE(image.shape() == view.shape()); + REQUIRE(image.size() == view.size()); + REQUIRE(image.data() != view.data()); + + for(int64_t i=0; i shape{{20}}; NDArray img(shape, 3); diff --git a/src/RawFile.cpp b/src/RawFile.cpp index e704add..78cb6c5 100644 --- a/src/RawFile.cpp +++ b/src/RawFile.cpp @@ -76,8 +76,7 @@ size_t RawFile::n_mod() const { return n_subfile_parts; } size_t RawFile::bytes_per_frame() { - // return m_rows * m_cols * m_master.bitdepth() / 8; - return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / 8; + return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / bits_per_byte; } size_t RawFile::pixels_per_frame() { // return m_rows * m_cols; diff --git a/src/algorithm.test.cpp b/src/algorithm.test.cpp new file mode 100644 index 0000000..fcfa8d2 --- /dev/null +++ b/src/algorithm.test.cpp @@ -0,0 +1,73 @@ + + +#include +#include + + +TEST_CASE("Find the closed index in a 1D array", "[algorithm]") { + aare::NDArray arr({5}); + for (size_t i = 0; i < arr.size(); i++) { + arr[i] = i; + } + // arr 0, 1, 2, 3, 4 + REQUIRE(aare::nearest_index(arr, 2.3) == 2); + REQUIRE(aare::nearest_index(arr, 2.6) == 3); + REQUIRE(aare::nearest_index(arr, 45.0) == 4); + REQUIRE(aare::nearest_index(arr, 0.0) == 0); + REQUIRE(aare::nearest_index(arr, -1.0) == 0); +} + +TEST_CASE("Passing integers to nearest_index works", "[algorithm]"){ + aare::NDArray arr({5}); + for (size_t i = 0; i < arr.size(); i++) { + arr[i] = i; + } + // arr 0, 1, 2, 3, 4 + REQUIRE(aare::nearest_index(arr, 2) == 2); + REQUIRE(aare::nearest_index(arr, 3) == 3); + REQUIRE(aare::nearest_index(arr, 45) == 4); + REQUIRE(aare::nearest_index(arr, 0) == 0); + REQUIRE(aare::nearest_index(arr, -1) == 0); +} + + +TEST_CASE("nearest_index works with std::vector", "[algorithm]"){ + std::vector vec = {0, 1, 2, 3, 4}; + REQUIRE(aare::nearest_index(vec, 2.123) == 2); + REQUIRE(aare::nearest_index(vec, 2.66) == 3); + REQUIRE(aare::nearest_index(vec, 4555555.0) == 4); + REQUIRE(aare::nearest_index(vec, 0.0) == 0); + REQUIRE(aare::nearest_index(vec, -10.0) == 0); +} + +TEST_CASE("nearest index works with std::array", "[algorithm]"){ + std::array arr = {0, 1, 2, 3, 4}; + REQUIRE(aare::nearest_index(arr, 2.123) == 2); + REQUIRE(aare::nearest_index(arr, 2.501) == 3); + REQUIRE(aare::nearest_index(arr, 4555555.0) == 4); + REQUIRE(aare::nearest_index(arr, 0.0) == 0); + REQUIRE(aare::nearest_index(arr, -10.0) == 0); +} + + +TEST_CASE("last smaller", "[algorithm]"){ + aare::NDArray arr({5}); + for (size_t i = 0; i < arr.size(); i++) { + arr[i] = i; + } + // arr 0, 1, 2, 3, 4 + REQUIRE(aare::last_smaller(arr, -10.0) == 0); + REQUIRE(aare::last_smaller(arr, 0.0) == 0); + REQUIRE(aare::last_smaller(arr, 2.3) == 2); + REQUIRE(aare::last_smaller(arr, 253.) == 4); +} + +TEST_CASE("returns last bin strictly smaller", "[algorithm]"){ + aare::NDArray arr({5}); + for (size_t i = 0; i < arr.size(); i++) { + arr[i] = i; + } + // arr 0, 1, 2, 3, 4 + REQUIRE(aare::last_smaller(arr, 2.0) == 2); + +} \ No newline at end of file