detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-01-01 18:01:24 +01:00
Code Issues Actions Packages Releases Activity

refactored interpolation function transform_eta_values is public function

Browse Source
This commit is contained in:
mazzol_a
2025-12-11 16:49:10 +01:00
parent 80a2b02345
commit c4d1a8f3f4
5 changed files with 104 additions and 255 deletions
Download Patch File Download Diff File Expand all files Collapse all files

184
include/aare/Interpolator.hpp
View File

@@ -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<double, 3> 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<double, 3> m_ietay;
NDArray<double, 1> 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 <auto EtaFunction = calculate_eta2, typename ClusterType,
typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
std::vector<Photon>
interpolate(const ClusterVector<ClusterType> &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 <typename T>
Coordinate2D transform_eta_values(const Eta2<T> &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 <auto EtaFunction, typename ClusterType>
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 <typename T>
std::pair<double, double>
bilinear_interpolation(const size_t ix, const size_t iy, const size_t ie,
const Eta2<T> &eta);
const Eta2<T> &eta) const;
};
template <typename T>
std::pair<double, double>
Interpolator::bilinear_interpolation(const size_t ix, const size_t iy,
const size_t ie, const Eta2<T> &eta) {
const size_t ie,
const Eta2<T> &eta) const {
auto next_index_y = static_cast<ssize_t>(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 <typename T>
Coordinate2D Interpolator::transform_eta_values(const Eta2<T> &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<double>(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 <auto EtaFunction, typename ClusterType, typename Enable>
std::vector<Photon>
Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) const {
std::vector<Photon> photons;
photons.reserve(clusters.size());
@@ -159,28 +194,48 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(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<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type> ||
EtaFunction ==
&calculate_full_eta2<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type>) {
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<EtaFunction, ClusterType>(
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<ClusterType> &clusters) {
return photons;
}
template <auto EtaFunction, typename ClusterType>
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<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type> ||
EtaFunction == &calculate_full_eta2<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type>) {
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

164
python/src/interpolation.hpp
View File

@@ -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 <cstdint>
#include <filesystem>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
namespace py = pybind11;
#define REGISTER_INTERPOLATOR_ETA2(T, N, M, U) \
register_interpolate<T, N, M, U, aare::calculate_full_eta2<T, N, M, U>>( \
interpolator, "_full_eta2", "full eta2"); \
register_interpolate<T, N, M, U, aare::calculate_eta2<T, N, M, U>>( \
interpolator, "", "eta2");
#define REGISTER_INTERPOLATOR_ETA3(T, N, M, U) \
register_interpolate<T, N, M, U, aare::calculate_eta3<T, N, M, U>>( \
interpolator, "_eta3", "full eta3"); \
register_interpolate<T, N, M, U, aare::calculate_cross_eta3<T, N, M, U>>( \
interpolator, "_cross_eta3", "cross eta3");
template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
typename CoordType = uint16_t, auto EtaFunction>
void register_interpolate(py::class_<aare::Interpolator> &interpolator,
const std::string &typestr = "",
const std::string &doc_string_etatype = "eta2x2") {
using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
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<ClusterType> &clusters) {
auto photons = self.interpolate<EtaFunction, ClusterType>(clusters);
auto *ptr = new std::vector<Photon>{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_<aare::Interpolator>(m, "Interpolator")
.def(py::init(
[](py::array_t<double,
py::array::c_style | py::array::forcecast>
etacube,
py::array_t<double> xbins, py::array_t<double> ybins,
py::array_t<double> 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<double> xbins, py::array_t<double> ybins,
py::array_t<double> 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<double,
py::array::c_style | py::array::forcecast>
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<double, 3>{};
*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<double, 3>{};
*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.
)");
}

2
python/src/module.cpp
View File

@@ -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"

5
python/tests/test_Cluster.py
View File

@@ -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

4
src/Interpolator.cpp
View File

@@ -20,8 +20,8 @@ Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
m_ietax = NDArray<double, 3>(etacube);
m_ietay = NDArray<double, 3>(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++) {