diff --git a/include/aare/Interpolator.hpp b/include/aare/Interpolator.hpp index 913e6e2..72e9f9c 100644 --- a/include/aare/Interpolator.hpp +++ b/include/aare/Interpolator.hpp @@ -17,11 +17,27 @@ struct Photon { double energy; }; +struct Coordinate2D { + double x{}; + double y{}; +}; + class Interpolator { - // marginal CDF of eta_x (if rosenblatt applied), conditional - // CDF of eta_x conditioned on eta_y + + /** + * @brief + * marginal CDF of eta_x (if rosenblatt applied), conditional + * CDF of eta_x conditioned on eta_y + * value at (i, j, e): F(eta_x[i] | + *eta_y[j], energy[e]) + */ NDArray m_ietax; - // conditional CDF of eta_y conditioned on eta_x + + /** + * @brief + * conditional CDF of eta_y conditioned on eta_x + * value at (i,j,e): F(eta_y[j] | eta_x[i], energy[e]) + */ NDArray m_ietay; NDArray m_etabinsx; @@ -69,27 +85,26 @@ class Interpolator { * calculate_eta2 * @return interpolated photons (photon positions are given as double but * following row column format e.g. x=0, y=0 means top row and first column - * of frame) + * of frame) (An interpolated photon position of (1.5, 2.5) means that the + * photon hit is estimated to be the pixel center of pixel (1,2)) */ template >> - std::vector interpolate(const ClusterVector &clusters); + std::vector + interpolate(const ClusterVector &clusters) const; + + /** + * @brief transforms the eta values to uniform coordinates based on the CDF + * ieta_x and ieta_y + * @tparam eta Eta to transform + * @return uniform coordinates {x,y} + */ + template + Coordinate2D transform_eta_values(const Eta2 &eta) const; private: /** - * @brief implements underlying interpolation logic based on EtaFunction - * Type - * @tparam EtaFunction Function object that calculates desired eta default: - * @param u: transformed photon position in x between [0,1] - * @param v: transformed photon position in y between [0,1] - * @param c: corner of eta - */ - template - void interpolation_logic(Photon &photon, const double u, const double v, - const corner c = corner::cTopLeft); - - /** - * @brief bilinear interpolation of the transformed eta values + * @brief bilinear interpolation of the transformed eta values * @param ix index of etaX bin * @param iy index of etaY bin * @param ie index of energy bin @@ -98,13 +113,14 @@ class Interpolator { template std::pair bilinear_interpolation(const size_t ix, const size_t iy, const size_t ie, - const Eta2 &eta); + const Eta2 &eta) const; }; template std::pair Interpolator::bilinear_interpolation(const size_t ix, const size_t iy, - const size_t ie, const Eta2 &eta) { + const size_t ie, + const Eta2 &eta) const { auto next_index_y = static_cast(iy + 1) >= m_ietax.shape(1) ? m_ietax.shape(1) - 1 : iy + 1; @@ -144,9 +160,28 @@ Interpolator::bilinear_interpolation(const size_t ix, const size_t iy, return {ietax_interpolated, ietay_interpolated}; } +template +Coordinate2D Interpolator::transform_eta_values(const Eta2 &eta) const { + + // 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, static_cast(eta.sum)); + auto ix = last_smaller(m_etabinsx, eta.x); + auto iy = last_smaller(m_etabinsy, eta.y); + + // TODO: bilinear interpolation only works if all bins have a size > 1 - + // otherwise bilinear interpolation with zero values which skew the + // results + // TODO: maybe trim the bins at the edges with zero values beforehand + // auto [ietax_interpolated, ietay_interpolated] = + // bilinear_interpolation(ix, iy, ie, eta); + + return Coordinate2D{m_ietax(ix, iy, ie), m_ietay(ix, iy, ie)}; +} + template std::vector -Interpolator::interpolate(const ClusterVector &clusters) { +Interpolator::interpolate(const ClusterVector &clusters) const { std::vector photons; photons.reserve(clusters.size()); @@ -159,28 +194,48 @@ Interpolator::interpolate(const ClusterVector &clusters) { photon.y = cluster.y; photon.energy = static_cast(eta.sum); - // std::cout << "eta.x: " << eta.x << " eta.y: " << eta.y << std::endl; + auto uniform_coordinates = transform_eta_values(eta); - // 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); + if (EtaFunction == &calculate_eta2 || + EtaFunction == + &calculate_full_eta2) { + double dX{}, dY{}; - // std::cout << "ix: " << ix << " iy: " << iy << std::endl; - - // TODO: bilinear interpolation only works if all bins have a size > 1 - - // otherwise bilinear interpolation with zero values which skew the - // results - // TODO: maybe trim the bins at the edges with zero values beforehand - // auto [ietax_interpolated, ietay_interpolated] = - // bilinear_interpolation(ix, iy, ie, eta); - - double ietax_interpolated = m_ietax(ix, iy, ie); - double ietay_interpolated = m_ietay(ix, iy, ie); - - interpolation_logic( - photon, ietax_interpolated, ietay_interpolated, eta.c); + // TODO: could also chaneg the sign of the eta calculation + switch (eta.c) { + case corner::cTopLeft: + dX = -1.0; + dY = -1.0; + break; + case corner::cTopRight:; + dX = 0.0; + dY = -1.0; + break; + case corner::cBottomLeft: + dX = -1.0; + dY = 0.0; + break; + case corner::cBottomRight: + dX = 0.0; + dY = 0.0; + break; + } + photon.x = photon.x + 0.5 + uniform_coordinates.x + + dX; // use pixel center + 0.5 + photon.y = + photon.y + 0.5 + uniform_coordinates.y + + dY; // eta2 calculates the ratio between bottom and sum of + // bottom and top shift by 1 add eta value correctly + } else { + photon.x += uniform_coordinates.x; + photon.y += uniform_coordinates.y; + } photons.push_back(photon); } @@ -188,51 +243,4 @@ Interpolator::interpolate(const ClusterVector &clusters) { return photons; } -template -void Interpolator::interpolation_logic(Photon &photon, const double u, - const double v, const corner c) { - - // std::cout << "u: " << u << " v: " << v << std::endl; - - // TODO: try to call this with std::is_same_v and have it constexpr if - // possible - if (EtaFunction == &calculate_eta2 || - EtaFunction == &calculate_full_eta2) { - double dX{}, dY{}; - - // TODO: could also chaneg the sign of the eta calculation - switch (c) { - case corner::cTopLeft: - dX = -1.0; - dY = -1.0; - break; - case corner::cTopRight:; - dX = 0.0; - dY = -1.0; - break; - case corner::cBottomLeft: - dX = -1.0; - dY = 0.0; - break; - case corner::cBottomRight: - dX = 0.0; - dY = 0.0; - break; - } - photon.x = photon.x + 0.5 + u + dX; // use pixel center + 0.5 - photon.y = photon.y + 0.5 + v + - dY; // eta2 calculates the ratio between bottom and sum of - // bottom and top shift by 1 add eta value correctly - } else { - photon.x += u; - photon.y += v; - } -} - } // namespace aare \ No newline at end of file diff --git a/python/src/interpolation.hpp b/python/src/interpolation.hpp deleted file mode 100644 index a9d79c3..0000000 --- a/python/src/interpolation.hpp +++ /dev/null @@ -1,164 +0,0 @@ -// SPDX-License-Identifier: MPL-2.0 -#include "aare/CalculateEta.hpp" -#include "aare/Interpolator.hpp" -#include "aare/NDArray.hpp" -#include "aare/NDView.hpp" -#include "np_helper.hpp" -#include -#include -#include -#include - -namespace py = pybind11; - -#define REGISTER_INTERPOLATOR_ETA2(T, N, M, U) \ - register_interpolate>( \ - interpolator, "_full_eta2", "full eta2"); \ - register_interpolate>( \ - interpolator, "", "eta2"); - -#define REGISTER_INTERPOLATOR_ETA3(T, N, M, U) \ - register_interpolate>( \ - interpolator, "_eta3", "full eta3"); \ - register_interpolate>( \ - interpolator, "_cross_eta3", "cross eta3"); - -template -void register_interpolate(py::class_ &interpolator, - const std::string &typestr = "", - const std::string &doc_string_etatype = "eta2x2") { - - using ClusterType = Cluster; - - const std::string docstring = "interpolation based on " + - doc_string_etatype + - "\n\nReturns:\n interpolated photons"; - - auto function_name = fmt::format("interpolate{}", typestr); - - interpolator.def( - function_name.c_str(), - [](aare::Interpolator &self, - const ClusterVector &clusters) { - auto photons = self.interpolate(clusters); - auto *ptr = new std::vector{photons}; - return return_vector(ptr); - }, - docstring.c_str(), py::arg("cluster_vector")); -} - -void define_interpolation_bindings(py::module &m) { - - PYBIND11_NUMPY_DTYPE(aare::Photon, x, y, energy); - - auto interpolator = - 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)); - }), - R"doc( - Constructor - - Args: - - etacube: - joint distribution of eta_x, eta_y and photon energy (**Note:** for the joint distribution first dimension is eta_x, second: eta_y, third: energy bins.) - xbins: - bin edges of etax - ybins: - bin edges of etay - ebins: - bin edges of photon energy - )doc", - py::arg("etacube"), - py::arg("xbins"), py::arg("ybins"), - py::arg("ebins")) - - .def(py::init( - [](py::array_t xbins, py::array_t ybins, - py::array_t ebins) { - return Interpolator(make_view_1d(xbins), - make_view_1d(ybins), - make_view_1d(ebins)); - }), - R"( - Constructor - - Args: - - xbins: - bin edges of etax - ybins: - bin edges of etay - ebins: - bin edges of photon energy - )", py::arg("xbins"), - py::arg("ybins"), py::arg("ebins")) - .def( - "rosenblatttransform", - [](Interpolator &self, - py::array_t - etacube) { - return self.rosenblatttransform(make_view_3d(etacube)); - }, - R"( - calculated the rosenblatttransform for the given distribution - - etacube: - joint distribution of eta_x, eta_y and photon energy (**Note:** for the joint distribution first dimension is eta_x, second: eta_y, third: energy bins.) - )", - py::arg("etacube")) - .def("get_ietax", - [](Interpolator &self) { - auto *ptr = new NDArray{}; - *ptr = self.get_ietax(); - return return_image_data(ptr); - }, R"(conditional CDF of etax conditioned on etay, marginal CDF of etax (if rosenblatt transform applied))") - .def("get_ietay", [](Interpolator &self) { - auto *ptr = new NDArray{}; - *ptr = self.get_ietay(); - return return_image_data(ptr); - }, R"(conditional CDF of etay conditioned on etax)"); - - REGISTER_INTERPOLATOR_ETA3(int, 3, 3, uint16_t); - REGISTER_INTERPOLATOR_ETA3(float, 3, 3, uint16_t); - REGISTER_INTERPOLATOR_ETA3(double, 3, 3, uint16_t); - - REGISTER_INTERPOLATOR_ETA2(int, 3, 3, uint16_t); - REGISTER_INTERPOLATOR_ETA2(float, 3, 3, uint16_t); - REGISTER_INTERPOLATOR_ETA2(double, 3, 3, uint16_t); - REGISTER_INTERPOLATOR_ETA2(int, 2, 2, uint16_t); - REGISTER_INTERPOLATOR_ETA2(float, 2, 2, uint16_t); - REGISTER_INTERPOLATOR_ETA2(double, 2, 2, uint16_t); - - // 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 984e9f7..f78d2d4 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -10,13 +10,13 @@ #include "bind_ClusterFinderMT.hpp" #include "bind_ClusterVector.hpp" #include "bind_Eta.hpp" +#include "bind_Interpolator.hpp" #include "bind_calibration.hpp" // TODO! migrate the other names #include "ctb_raw_file.hpp" #include "file.hpp" #include "fit.hpp" -#include "interpolation.hpp" #include "jungfrau_data_file.hpp" #include "pedestal.hpp" #include "pixel_map.hpp" diff --git a/python/tests/test_Cluster.py b/python/tests/test_Cluster.py index 944dc3f..cafb118 100644 --- a/python/tests/test_Cluster.py +++ b/python/tests/test_Cluster.py @@ -51,6 +51,11 @@ def test_Interpolator(): cluster = _aare.Cluster3x3i(1,1, np.ones(9, dtype=np.int32)) clustervector.push_back(cluster) + [u,v] = interpolator.transform_eta_values(_aare.Etai()) + + assert u == 0 + assert v == 0 + interpolated_photons = interpolator.interpolate(clustervector) assert interpolated_photons.size == 1 diff --git a/src/Interpolator.cpp b/src/Interpolator.cpp index 1240107..7d95994 100644 --- a/src/Interpolator.cpp +++ b/src/Interpolator.cpp @@ -20,8 +20,8 @@ Interpolator::Interpolator(NDView etacube, NDView xbins, m_ietax = NDArray(etacube); m_ietay = NDArray(etacube); - - // prefix sum - conditional CDF + // TODO: etacube should have different strides energy should come first + // prefix sum - conditional CDF 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++) {