Dev/interpolation documentation (#255)

- added transform_eta_values for easier debugging more control for the
user
- updated Documentation

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