renamed interpolator.hpp

RELEASE.md

@@ -6,6 +6,7 @@
 
 - Expanding 24 to 32 bit data
 - Decoding digital data from Mythen 302
+- added ``transform_eta_values``. Function transforms :math:`\eta` to uniform spatial coordinates. Should only be used for easier debugging. 
 
 
 ### 2025.11.21

docs/src/Interpolation.rst

@@ -1,12 +1,43 @@
+.. _Interpolation_C++API:
+
 Interpolation
 ==============
 
-Interpolation class for :math:`\eta` Interpolation. 
+The Interpolation class implements the :math:`\eta`-interpolation method.
+This interpolation technique is based on charge sharing: for detected photon hits (e.g. clusters), it refines the estimated photon hit using information from neighboring pixels.
 
-The Interpolator class provides methods to interpolate the positions of photons based on their :math:`\eta` values. 
+The method relies on the so-called :math:`\eta`-functions, which describe the relationship between the energy measured in the central cluster pixel (the initially estimated photon hit) and the energies measured in its neighboring pixels.
+Depending on how much energy each neighboring pixel receives relative to the central pixel, the estimated photon hit is shifted toward that neighbor by a certain offset to the actual photon hit position in the pixel :math:`(x, y)`.
+
+The mapping between the :math:`\eta` values and the corresponding spatial photon position :math:`(x,y)` can be viewed as an optimal transport problem.
+
+One can readily compute the probability distribution :math:`P_{\eta}` of the :math:`\eta` values by forming a 2D histogram.
+However, the probability distribution :math:`P_{x,y}` of the true photon positions is generally unknown unless the detector is illuminated uniformly (i.e. under flat-field conditions).
+In a flat-field, the photon positions are uniformly distributed.
+
+With this assumption, the problem reduces to determining a transport map :math:`T:(\eta_x,\eta_y) \rightarrow (x,y)`, that pushes forward the distribution of :math:`(\eta_x, \eta_y)` to the known uniform distribution of photon positions of a flatfield.
+
+The map :math:`T` is given by: 
+
+.. math:: 
+    \begin{align*}
+        T_1: & F_{x}^{-1} F_{\eta_x|\eta_y} \\
+        T_2: & F_{y}^{-1} F_{\eta_y|\eta_x}, 
+    \end{align*}
+
+
+where :math:`F_{\eta_x|\eta_y}` and :math:`F_{\eta_y|\eta_x}` are the conditional cumulative distribution functions e.g. :math:`F_{\eta_x|\eta_y}(\eta_x', \eta_y') = P_{\eta_x, \eta_y}(\eta_x \leq \eta_x' | \eta_y = \eta_y')`. 
+And :math:`F_{x}` and :math:`F_{y}` are the cumulative distribution functions of :math:`x` and :math:`y`. Note as :math:`x` and :math:`y` are uniformly distributed :math:`F_{x}` and :math:`F_{y}` are the identity functions. The map :math:`T` thus simplifies to 
+
+.. math:: 
+    \begin{align*}
+        T_1: & F_{\eta_x|\eta_y} \\
+        T_2: & F_{\eta_y|\eta_x}. 
+    \end{align*}
+
+Note that for the implementation :math:`P_{\eta}` is not only a distribution of :math:`\eta_x`, :math:`\eta_y` but also of the estimated photon energy :math:`e`. 
+The energy level correlates slightly with the z-depth. Higher z-depth leads to more charge sharing and a different :math:`\eta` distribution. Thus we create a mapping :math:`T` for each energy level. 
 
-.. warning:: 
-   The interpolation might lead to erroneous photon positions for clusters at the boarders of a frame. Make sure to filter out such cases. 
 
 :math:`\eta`-Functions: 
 ---------------------------
@@ -21,6 +52,8 @@ The Interpolator class provides methods to interpolate the positions of photons
 
 Supported are the following :math:`\eta`-functions: 
 
+:math:`\eta`-Function on 2x2 Clusters: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 .. image:: ../figures/Eta2x2.png
     :target: ../figures/Eta2x2.png
@@ -34,11 +67,18 @@ Supported are the following :math:`\eta`-functions:
    {\color{green}{\eta_y}} = \frac{Q_{1,1}}{Q_{0,1} + Q_{1,1}}
    \end{equation*}
 
+The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change. 
+If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1.
+
+One can apply this :math:`\eta` not only on 2x2 clusters but on clusters with any size. Then the 2x2 subcluster with maximum energy is choosen and the :math:`\eta` function applied on the subcluster.
 
 .. doxygenfunction:: aare::calculate_eta2(const ClusterVector<ClusterType>&)
 
 .. doxygenfunction:: aare::calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&) 
 
+Full :math:`\eta`-Function on 2x2 Clusters: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 .. image:: ../figures/Eta2x2Full.png
     :target: ../figures/Eta2x2Full.png
     :width: 650px
@@ -47,15 +87,20 @@ Supported are the following :math:`\eta`-functions:
 
 .. math:: 
     \begin{equation*}
-        {\color{blue}{\eta_x}}  = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}} \quad \quad
-        {\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}}
+        {\color{blue}{\eta_x}}  = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{1}\sum_j^{1}Q_{i,j}} \quad \quad
+        {\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{1}\sum_j^{1}Q_{i,j}}
     \end{equation*} 
 
+The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change. 
+If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1. 
 
 .. doxygenfunction:: aare::calculate_full_eta2(const ClusterVector<ClusterType>&)
 
 .. doxygenfunction:: aare::calculate_full_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
 
+Full :math:`\eta`-Function on 3x3 Clusters: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 .. image:: ../figures/Eta3x3.png
     :target: ../figures/Eta3x3.png
     :width: 650px
@@ -64,13 +109,18 @@ Supported are the following :math:`\eta`-functions:
 
 .. math::
     \begin{equation*}
-         {\color{blue}{\eta_x}} = \frac{\sum_{i}^{3} Q_{i,2} - \sum_{i}^{3} Q_{i,0}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}} \quad \quad
-        {\color{green}{\eta_y}} = \frac{\sum_{j}^{3} Q_{2,j} - \sum_{j}^{3} Q_{0,j}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}}
+         {\color{blue}{\eta_x}} = \frac{\sum_{i=0}^{2} Q_{i,2} - \sum_{i=0}^{2} Q_{i,0}}{\sum_{i=0}^{2}\sum_{j=0}^{2} Q_{i,j}} \quad \quad
+        {\color{green}{\eta_y}} = \frac{\sum_{j=0}^{2} Q_{2,j} - \sum_{j=0}^{2} Q_{0,j}}{\sum_{i=0}^{2}\sum_{j=0}^{2} Q_{i,j}}
     \end{equation*}
 
-.. doxygenfunction:: aare::calculate_eta3(const ClusterVector<Cluster<T, 3,3, CoordType>>&)
+The :math:`\eta` values can range between -0.5,0.5.
 
-.. doxygenfunction:: aare::calculate_eta3(const Cluster<T, 3, 3, CoordType>&)
+.. doxygenfunction:: aare::calculate_eta3(const ClusterVector<ClusterType>&)
+
+.. doxygenfunction:: aare::calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
+
+Cross :math:`\eta`-Function on 3x3 Clusters: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 .. image:: ../figures/Eta3x3Cross.png
     :target: ../figures/Eta3x3Cross.png
@@ -80,20 +130,28 @@ Supported are the following :math:`\eta`-functions:
 
 .. math:: 
     \begin{equation*}
-       {\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} +  Q_{1,1} + Q_{1,0}} \quad \quad
-        {\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} +  Q_{1,1} + Q_{1,2}}
+       {\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} +  Q_{1,1} + Q_{1,2}} \quad \quad
+        {\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} +  Q_{1,1} + Q_{2,1}}
     \end{equation*}
 
-.. doxygenfunction:: aare::calculate_cross_eta3(const ClusterVector<Cluster<T, 3,3, CoordType>>&)
+The :math:`\eta` values can range between -0.5,0.5. 
 
-.. doxygenfunction:: aare::calculate_cross_eta3(const Cluster<T, 3, 3, CoordType>&)
+.. doxygenfunction:: aare::calculate_cross_eta3(const ClusterVector<ClusterType>&)
+
+.. doxygenfunction:: aare::calculate_cross_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>&)
 
 Interpolation class: 
 ---------------------
 
+.. warning:: 
+   The interpolation might lead to erroneous photon positions for clusters at the borders of a frame. Make sure to filter out such cases. 
+
 .. Warning:: 
     Make sure to use the same :math:`\eta`-function during interpolation as given by the joint :math:`\eta`-distribution passed to the constructor. 
 
+.. Note:: 
+    Make sure to use resonable energy bins, when constructing the joint distribution. If data is too sparse for a given energy the interpolation will lead to erreneous results.
+
 .. doxygenclass:: aare::Interpolator
    :members:
    :undoc-members:

docs/src/python/cluster/pyInterpolation.rst

@@ -1,12 +1,13 @@
 Interpolation 
 ==============
 
-Interpolation class for :math:`\eta` Interpolation. 
+The Interpolation class implements the :math:`\eta`-interpolation method.
+This interpolation technique is based on charge sharing: for detected photon hits (e.g. clusters), it refines the estimated photon hit using information from neighboring pixels.
 
-The Interpolator class provides methods to interpolate the positions of photons based on their :math:`\eta` values. 
+See :ref:`Interpolation_C++API` for a more elaborate documentation and explanation of the method. 
 
-.. warning:: 
-   The interpolation might lead to erroneous photon positions for clusters at the boarders of a frame. Make sure to filter out such cases. 
+:math:`\eta`-Functions: 
+--------------------------
 
 Below is an example of the Eta class of type ``double``. Supported are ``Etaf`` of type ``float`` and ``Etai`` of type ``int``. 
 
@@ -19,6 +20,9 @@ Below is an example of the Eta class of type ``double``. Supported are ``Etaf``
 
 Supported are the following :math:`\eta`-functions: 
 
+:math:`\eta`-Function on 2x2 Clusters: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 .. py:currentmodule:: aare
 
 .. image:: ../../../figures/Eta2x2.png
@@ -33,8 +37,14 @@ Supported are the following :math:`\eta`-functions:
    {\color{green}{\eta_y}} = \frac{Q_{1,1}}{Q_{0,1} + Q_{1,1}}
    \end{equation*}    
 
+The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change. 
+If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1. 
+
 .. autofunction:: calculate_eta2
 
+Full :math:`\eta`-Function on 2x2 Clusters: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 .. image:: ../../../figures/Eta2x2Full.png
     :target: ../../../figures/Eta2x2Full.png
     :width: 650px
@@ -43,12 +53,18 @@ Supported are the following :math:`\eta`-functions:
 
 .. math:: 
     \begin{equation*}
-        {\color{blue}{\eta_x}}  = \frac{Q_{0,1} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}} \quad \quad
-        {\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_i^{2}\sum_j^{2}Q_{i,j}}
+        {\color{blue}{\eta_x}}  = \frac{Q_{0,1} + Q_{1,1}}{\sum_{i=0}^{1}\sum_{j=0}^{1}Q_{i,j}} \quad \quad
+        {\textcolor{green}{\eta_y}} = \frac{Q_{1,0} + Q_{1,1}}{\sum_{i=0}^{1}\sum_{j=0}^{1}Q_{i,j}}
     \end{equation*} 
 
+The :math:`\eta` values can range between 0,1. Note they only range between 0,1 because the position of the center pixel (red) can change. 
+If the center pixel is in the bottom left pixel :math:`\eta_x` will be close to zero. If the center pixel is in the bottom right pixel :math:`\eta_y` will be close to 1. 
+
 .. autofunction:: calculate_full_eta2
 
+Full :math:`\eta`-Function on 3x3 Clusters: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 .. image:: ../../../figures/Eta3x3.png
     :target: ../../../figures/Eta3x3.png
     :width: 650px
@@ -57,12 +73,17 @@ Supported are the following :math:`\eta`-functions:
 
 .. math::
     \begin{equation*}
-         {\color{blue}{\eta_x}} = \frac{\sum_{i}^{3} Q_{i,2} - \sum_{i}^{3} Q_{i,0}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}} \quad \quad
-        {\color{green}{\eta_y}} = \frac{\sum_{j}^{3} Q_{2,j} - \sum_{j}^{3} Q_{0,j}}{\sum_{i}^{3}\sum_{j}^{3} Q_{i,j}}
+         {\color{blue}{\eta_x}} = \frac{\sum_{i=0}^{2} Q_{i,2} - \sum_{i=0}^{2} Q_{i,0}}{\sum_{i=0}^{2}\sum_{j}^{3} Q_{i,j}} \quad \quad
+        {\color{green}{\eta_y}} = \frac{\sum_{j=0}^{2} Q_{2,j} - \sum_{j=0}^{2} Q_{0,j}}{\sum_{i=0}^{2}\sum_{j}^{3} Q_{i,j}}
     \end{equation*}
 
+The :math:`\eta` values can range between -0.5,0.5. 
+
 .. autofunction:: calculate_eta3
 
+Cross :math:`\eta`-Function on 3x3 Clusters: 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 .. image:: ../../../figures/Eta3x3Cross.png
     :target: ../../../figures/Eta3x3Cross.png
     :width: 650px
@@ -71,10 +92,12 @@ Supported are the following :math:`\eta`-functions:
 
 .. math:: 
     \begin{equation*}
-       {\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} +  Q_{1,1} + Q_{1,0}} \quad \quad
-        {\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} +  Q_{1,1} + Q_{1,2}}
+       {\color{blue}{\eta_x}} = \frac{Q_{1,2} - Q_{1,0}}{Q_{1,0} +  Q_{1,1} + Q_{1,2}} \quad \quad
+        {\color{green}{\eta_y}} = \frac{Q_{0,2} - Q_{0,1}}{Q_{0,1} +  Q_{1,1} + Q_{2,1}}
     \end{equation*}
 
+The :math:`\eta` values can range between -0.5,0.5. 
+
 .. autofunction:: calculate_cross_eta3
 
 
@@ -84,6 +107,13 @@ Interpolation class for :math:`\eta`-Interpolation
 .. Warning:: 
     Make sure to use the same :math:`\eta`-function during interpolation as given by the joint :math:`\eta`-distribution passed to the constructor. 
 
+.. Warning:: 
+   The interpolation might lead to erroneous photon positions for clusters at the boarders of a frame. Make sure to filter out such cases.
+
+.. Note:: 
+    Make sure to use resonable energy bins, when constructing the joint distribution. If data is too sparse for a given energy the interpolation will lead to erreneous results.
+
+
 .. py:currentmodule:: aare
 
 .. autoclass:: Interpolator

include/aare/CalculateEta.hpp

@@ -29,7 +29,8 @@ template <typename T> struct Eta2 {
     double x{};
     /// @brief eta in y direction
     double y{};
-    /// @brief index of subcluster given as corner relative to cluster center
+    /// @brief index of subcluster with highest energy value (given as corner
+    /// relative to cluster center)
     corner c{0};
     /// @brief photon energy (cluster sum)
     T sum{};

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<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;
@@ -31,11 +47,11 @@ class Interpolator {
   public:
     /**
      * @brief Constructor for the Interpolator class
-     * @param etacube joint distribution of etaX, etaY and photon energy
+     * @param etacube joint distribution of etaX, etaY and photon energy (note
+     * first dimension is etaX, second etaY, third photon energy)
      * @param xbins bin edges for etaX
      * @param ybins bin edges for etaY
      * @param ebins bin edges for photon energy
-     * @note note first dimension is etaX, second etaY, third photon energy
      */
     Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
                  NDView<double, 1> ybins, NDView<double, 1> ebins);
@@ -53,8 +69,8 @@ class Interpolator {
      * @brief transforms the joint eta distribution of etaX and etaY to the two
      * independant uniform distributions based on the Roseblatt transform for
      * each energy level
-     * @param etacube joint distribution of etaX, etaY and photon energy
-     * @note note first dimension is etaX, second etaY, third photon energy
+     * @param etacube joint distribution of etaX, etaY and photon energy (first
+     * dimension is etaX, second etaY, third photon energy)
      */
     void rosenblatttransform(NDView<double, 3> etacube);
 
@@ -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) corresponds to
+     * an estimated photon hit at 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

python/src/bind_Interpolator.hpp

@@ -48,6 +48,19 @@ void register_interpolate(py::class_<aare::Interpolator> &interpolator,
         docstring.c_str(), py::arg("cluster_vector"));
 }
 
+template <typename Type>
+void register_transform_eta_values(
+    py::class_<aare::Interpolator> &interpolator) {
+    interpolator.def(
+        "transform_eta_values",
+        [](Interpolator &self, const Eta2<Type> &eta) {
+            auto uniform_coord = self.transform_eta_values(eta);
+            return py::make_tuple(uniform_coord.x, uniform_coord.y);
+        },
+        R"(eta.x eta.y transformed to uniform coordinates based on CDF ietax, ietay)",
+        py::arg("Eta"));
+}
+
 void define_interpolation_bindings(py::module &m) {
 
     PYBIND11_NUMPY_DTYPE(aare::Photon, x, y, energy);
@@ -140,6 +153,10 @@ void define_interpolation_bindings(py::module &m) {
     REGISTER_INTERPOLATOR_ETA2(float, 2, 2, uint16_t);
     REGISTER_INTERPOLATOR_ETA2(double, 2, 2, uint16_t);
 
+    register_transform_eta_values<int>(interpolator);
+    register_transform_eta_values<float>(interpolator);
+    register_transform_eta_values<double>(interpolator);
+
     // TODO! Evaluate without converting to double
     m.def(
         "hej",

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"

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

src/Interpolator.cpp

@@ -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++) {