✨ Add interpolation classes to remove GSL.

src/pyinterp/core/include/pyinterp/detail/interpolation/akima.hpp

@@ -0,0 +1,135 @@
+// Copyright (c) 2024 CNES
+//
+// All rights reserved. Use of this source code is governed by a
+// BSD-style license that can be found in the LICENSE file.
+#pragma once
+#include "pyinterp/detail/interpolation/interpolator_1d.hpp"
+
+namespace pyinterp::detail::interpolation {
+
+/// Akima interpolation
+template <typename T>
+class Akima : public Interpolator1D<T> {
+ public:
+  using Interpolator1D<T>::Interpolator1D;
+  using Interpolator1D<T>::operator();
+  using Interpolator1D<T>::derivative;
+
+  /// Returns the minimum number of points required for the interpolation.
+  auto min_size() const -> Eigen::Index override { return 5; }
+
+ private:
+  Vector<T> m_{};
+  Vector<T> s_{};
+
+  /// Compute the boundary conditions.
+  virtual auto boundary_condition(T* m, const size_t size) -> void {
+    m[-2] = 3 * m[0] - 2 * m[1];
+    m[-1] = 2 * m[0] - m[1];
+    m[size - 1] = 2 * m[size - 2] - m[size - 3];
+    m[size] = 3 * m[size - 2] - 2 * m[size - 3];
+  }
+
+  /// @brief Compute the coefficients of the interpolation
+  /// @param xa X-coordinates of the data points.
+  /// @param ya Y-coordinates of the data points.
+  auto compute_coefficients(const Eigen::Ref<const Vector<T>>& xa,
+                            const Eigen::Ref<const Vector<T>>& ya)
+      -> void override;
+
+  /// Interpolation
+  /// @param xa X-coordinates of the data points.
+  /// @param ya Y-coordinates of the data points.
+  /// @param x The point where the interpolation must be calculated.
+  /// @param index The index of the last point found in the search.
+  /// @return The interpolated value at the point x.
+  auto operator()(const Eigen::Ref<const Vector<T>>& xa,
+                  const Eigen::Ref<const Vector<T>>& ya, const T& x,
+                  Eigen::Index* index) const -> T override;
+
+  /// @brief Returns the derivative of the interpolation function at the point
+  ///   x.
+  /// @param xa X-coordinates of the data points.
+  /// @param ya Y-coordinates of the data points.
+  /// @param x The point where the derivative must be calculated.
+  /// @param index The index of the last point found in the search.
+  /// @return The derivative of the interpolation function at the point x.
+  auto derivative(const Eigen::Ref<const Vector<T>>& xa,
+                  const Eigen::Ref<const Vector<T>>& ya, const T& x,
+                  Eigen::Index* index) const -> T override;
+};
+
+template <typename T>
+auto Akima<T>::compute_coefficients(const Eigen::Ref<const Vector<T>>& xa,
+                                    const Eigen::Ref<const Vector<T>>& ya)
+    -> void {
+  Interpolator1D<T>::compute_coefficients(xa, ya);
+  auto size = xa.size();
+  if (m_.size() < size + 4) {
+    m_.resize(size + 4);
+    s_.resize(size);
+  }
+
+  // m contains the slopes of the lines between the points. Two extra points
+  // are added at the beginning and end to handle the boundary conditions.
+  auto* m = m_.data() + 2;
+  for (Eigen::Index ix = 0; ix < size - 1; ++ix) {
+    m[ix] = (ya[ix + 1] - ya[ix]) / (xa[ix + 1] - xa[ix]);
+  }
+
+  boundary_condition(m, size);
+
+  // Compute the spline slopes of the lines between the points.
+  for (Eigen::Index ix = 2; ix < size - 2; ++ix) {
+    auto denominator =
+        std::abs(m[ix + 1] - m[ix]) + std::abs(m[ix - 1] - m[ix - 2]);
+    if (denominator != 0) {
+      s_(ix) = std::abs(m[ix + 1] - m[ix]) * m[ix - 1] +
+               std::abs(m[ix - 1] - m[ix - 2]) * m[ix];
+      s_(ix) /= denominator;
+    } else {
+      s_(ix) = (m[ix - 1] + m[ix]) * 0.5;
+    }
+  }
+  s_(0) = m[0];
+  s_(1) = (m[0] + m[2]) * 0.5;
+  s_(size - 2) = (m[size - 3] + m[size - 1]) * 0.5;
+  s_(size - 1) = m[size - 1];
+}
+
+template <typename T>
+auto Akima<T>::operator()(const Eigen::Ref<const Vector<T>>& xa,
+                          const Eigen::Ref<const Vector<T>>& ya, const T& x,
+                          Eigen::Index* index) const -> T {
+  auto search = this->search(xa, x, index);
+  if (!search) {
+    throw std::numeric_limits<T>::quiet_NaN();
+  }
+  auto [i0, i1] = *search;
+  const auto dx = xa(i1) - xa(i0);
+  const auto h = x - xa(i0);
+  const auto ai = ya(i0);
+  const auto bi = s_[i0];
+  const auto ci = (3 * m_[i0] - 2 * s_[i0] - s_[i1]) / dx;
+  const auto di = (s_[i0] + s_[i1] - 2 * m_[i0]) / (dx * dx);
+  return ai + h * (bi + h * (ci + h * di));
+}
+
+template <typename T>
+auto Akima<T>::derivative(const Eigen::Ref<const Vector<T>>& xa,
+                          const Eigen::Ref<const Vector<T>>& ya, const T& x,
+                          Eigen::Index* index) const -> T {
+  auto search = this->search(xa, x, index);
+  if (!search) {
+    throw std::numeric_limits<T>::quiet_NaN();
+  }
+  auto [i0, i1] = *search;
+  const auto dx = xa(i1) - xa(i0);
+  const auto h = x - xa(i0);
+  const auto bi = s_[i0];
+  const auto ci = (3 * m_[i0] - 2 * s_[i0] - s_[i1]) / dx;
+  const auto di = (s_[i0] + s_[i1] - 2 * m_[i0]) / (dx * dx);
+  return bi + h * (2 * ci + h * 3 * di);
+}
+
+}  // namespace pyinterp::detail::interpolation

src/pyinterp/core/include/pyinterp/detail/interpolation/akima_periodic.hpp

@@ -0,0 +1,26 @@
+// Copyright (c) 2024 CNES
+//
+// All rights reserved. Use of this source code is governed by a
+// BSD-style license that can be found in the LICENSE file.
+#pragma once
+#include "pyinterp/detail/interpolation/akima.hpp"
+
+namespace pyinterp::detail::interpolation {
+
+/// Akima periodic interpolation
+template <typename T>
+class AkimaPeriodic : public Akima<T> {
+ public:
+  using Akima<T>::Akima;
+
+ private:
+  /// Compute the boundary conditions.
+  auto boundary_condition(T* m, const size_t size) -> void override {
+    m[-2] = m[size - 3];
+    m[-1] = m[size - 2];
+    m[size - 1] = m[0];
+    m[size] = m[1];
+  }
+};
+
+}  // namespace pyinterp::detail::interpolation

src/pyinterp/core/include/pyinterp/detail/interpolation/bicubic.hpp

@@ -0,0 +1,181 @@
+// Copyright (c) 2024 CNES
+//
+// All rights reserved. Use of this source code is governed by a
+// BSD-style license that can be found in the LICENSE file.
+#pragma once
+
+#include "pyinterp/detail/interpolation/cspline.hpp"
+#include "pyinterp/detail/interpolation/interpolator_2d.hpp"
+
+namespace pyinterp::detail::interpolation {
+
+/// Bicubic interpolation
+/// @tparam T type of the data
+template <typename T>
+class Bicubic : public Interpolator2D<T> {
+ public:
+  using Interpolator2D<T>::Interpolator2D;
+  using Interpolator2D<T>::operator();
+
+  /// Returns the minimum number of points required for the interpolation.
+  auto min_size() const -> Eigen::Index override { return 2; }
+
+ private:
+  /// Interpolation
+  /// @param xa X-coordinates of the data points.
+  /// @param ya Y-coordinates of the data points.
+  /// @param za Z-values of the data points.
+  /// @param x The point where the interpolation must be calculated.
+  /// @param y The point where the interpolation must be calculated.
+  /// @param ix The index of the last point found in the search.
+  /// @param jx The index of the last point found in the search.
+  /// @return The interpolated value at the coordinates x, y.
+  auto operator()(const Eigen::Ref<const Vector<T>> &xa,
+                  const Eigen::Ref<const Vector<T>> &ya,
+                  const Eigen::Ref<const Matrix<T>> &za, const T &x, const T &y,
+                  Eigen::Index *ix, Eigen::Index *jx) const -> T override;
+
+  /// Compute the coefficients of the bicubic interpolation
+  /// @param xa X-coordinates of the data points.
+  /// @param ya Y-coordinates of the data points.
+  /// @param za Z-values of the data points.
+  auto compute_coefficients(const Eigen::Ref<const Vector<T>> &xa,
+                            const Eigen::Ref<const Vector<T>> &ya,
+                            const Eigen::Ref<const Matrix<T>> &za)
+      -> void override;
+
+  Matrix<T> zx_{};
+  Matrix<T> zy_{};
+  Matrix<T> zxy_{};
+};
+
+template <typename T>
+auto Bicubic<T>::compute_coefficients(const Eigen::Ref<const Vector<T>> &xa,
+                                      const Eigen::Ref<const Vector<T>> &ya,
+                                      const Eigen::Ref<const Matrix<T>> &za)
+    -> void {
+  Interpolator2D<T>::compute_coefficients(xa, ya, za);
+  auto xsize = xa.size();
+  auto ysize = ya.size();
+
+  if (zx_.rows() != xsize || zx_.cols() != ysize) {
+    zx_ = Matrix<T>(xsize, ysize);
+    zy_ = Matrix<T>(xsize, ysize);
+    zxy_ = Matrix<T>(xsize, ysize);
+  }
+
+  auto spline = CSpline<T>();
+  for (Eigen::Index j = 0; j < ysize; ++j) {
+    zx_.col(j) = spline.derivative(xa, za.col(j), xa);
+  }
+  for (Eigen::Index i = 0; i < xsize; ++i) {
+    zy_.row(i) = spline.derivative(ya, za.row(i), ya);
+  }
+  for (Eigen::Index j = 0; j < ysize; ++j) {
+    zxy_.col(j) = spline.derivative(xa, zy_.col(j), xa);
+  }
+}
+
+template <typename T>
+auto Bicubic<T>::operator()(const Eigen::Ref<const Vector<T>> &xa,
+                            const Eigen::Ref<const Vector<T>> &ya,
+                            const Eigen::Ref<const Matrix<T>> &za, const T &x,
+                            const T &y, Eigen::Index *ix,
+                            Eigen::Index *jx) const -> T {
+  auto search_x = this->search(xa, x, ix);
+  auto search_y = this->search(ya, y, jx);
+  if (!search_x || !search_y) {
+    throw std::numeric_limits<T>::quiet_NaN();
+  }
+  const auto [i0, i1] = *search_x;
+  const auto [j0, j1] = *search_y;
+  const auto x0 = xa(i0);
+  const auto x1 = xa(i1);
+  const auto y0 = ya(j0);
+  const auto y1 = ya(j1);
+  const auto z00 = za(i0, j0);
+  const auto z01 = za(i0, j1);
+  const auto z10 = za(i1, j0);
+  const auto z11 = za(i1, j1);
+  const auto dx = x1 - x0;
+  const auto dy = y1 - y0;
+  const auto t = (x - x0) / dx;
+  const auto u = (y - y0) / dy;
+  const auto zx00 = zx_(i0, j0) * dx;
+  const auto zx01 = zx_(i0, j1) * dx;
+  const auto zx10 = zx_(i1, j0) * dx;
+  const auto zx11 = zx_(i1, j1) * dx;
+  const auto zy00 = zy_(i0, j0) * dy;
+  const auto zy01 = zy_(i0, j1) * dy;
+  const auto zy10 = zy_(i1, j0) * dy;
+  const auto zy11 = zy_(i1, j1) * dy;
+  const auto zxy00 = zxy_(i0, j0) * (dx * dy);
+  const auto zxy01 = zxy_(i0, j1) * (dx * dy);
+  const auto zxy10 = zxy_(i1, j0) * (dx * dy);
+  const auto zxy11 = zxy_(i1, j1) * (dx * dy);
+  const auto t0 = 1;
+  const auto t1 = t;
+  const auto t2 = t * t;
+  const auto t3 = t * t2;
+  const auto u0 = 1;
+  const auto u1 = u;
+  const auto u2 = u * u;
+  const auto u3 = u * u2;
+
+  auto v = z00;
+  auto z = v * t0 * u0;
+
+  v = zy00;
+  z += v * t0 * u1;
+
+  v = 3 * (-z00 + z01) - 2 * zy00 - zy01;
+  z += v * t0 * u2;
+
+  v = 2 * (z00 - z01) + zy00 + zy01;
+  z += v * t0 * u3;
+
+  v = zx00;
+  z += v * t1 * u0;
+
+  v = zxy00;
+  z += v * t1 * u1;
+
+  v = 3 * (-zx00 + zx01) - 2 * zxy00 - zxy01;
+  z += v * t1 * u2;
+
+  v = 2 * (zx00 - zx01) + zxy00 + zxy01;
+  z += v * t1 * u3;
+
+  v = 3 * (-z00 + z10) - 2 * zx00 - zx10;
+  z += v * t2 * u0;
+
+  v = 3 * (-zy00 + zy10) - 2 * zxy00 - zxy10;
+  z += v * t2 * u1;
+
+  v = 9 * (z00 - z10 + z11 - z01) + 6 * (zx00 - zx01 + zy00 - zy10) +
+      3 * (zx10 - zx11 - zy11 + zy01) + 4 * zxy00 + 2 * (zxy10 + zxy01) + zxy11;
+  z += v * t2 * u2;
+
+  v = 6 * (-z00 + z10 - z11 + z01) + 4 * (-zx00 + zx01) +
+      3 * (-zy00 + zy10 + zy11 - zy01) + 2 * (-zx10 + zx11 - zxy00 - zxy01) -
+      zxy10 - zxy11;
+  z += v * t2 * u3;
+
+  v = 2 * (z00 - z10) + zx00 + zx10;
+  z += v * t3 * u0;
+
+  v = 2 * (zy00 - zy10) + zxy00 + zxy10;
+  z += v * t3 * u1;
+
+  v = 6 * (-z00 + z10 - z11 + z01) + 3 * (-zx00 - zx10 + zx11 + zx01) +
+      4 * (-zy00 + zy10) + 2 * (zy11 - zy01 - zxy00 - zxy10) - zxy11 - zxy01;
+  z += v * t3 * u2;
+
+  v = 4 * (z00 - z10 + z11 - z01) +
+      2 * (zx00 + zx10 - zx11 - zx01 + zy00 - zy10 - zy11 + zy01) + zxy00 +
+      zxy10 + zxy11 + zxy01;
+  z += v * t3 * u3;
+  return z;
+}
+
+}  // namespace pyinterp::detail::interpolation