Add spline interpolation

CMakeLists.txt

@@ -30,6 +30,7 @@ add_library(vineyard
     src/compressor.cpp
     src/cropper.cpp
     src/filter.cpp
+    src/interp.cpp
     src/log.cpp
     src/mapping.cpp
     src/parse_args.cpp

src/cmd/basis-img.cpp

@@ -2,11 +2,14 @@
 
 #include "algo/otsu.hpp"
 #include "basis.hpp"
+#include "interp.hpp"
 #include "io/hd5.hpp"
 #include "log.hpp"
 #include "parse_args.hpp"
 #include "tensorOps.hpp"
 
+#include <unsupported/Eigen/Splines>
+
 using namespace rl;
 
 int main_basis_img(args::Subparser &parser)
@@ -16,6 +19,9 @@ int main_basis_img(args::Subparser &parser)
 
   args::ValueFlag<Index> ss(parser, "S", "Subsample image", {"subsamp"}, 1);
   args::ValueFlag<Index> spf(parser, "S", "Spokes per frame", {"spf"}, 1);
+  args::ValueFlag<Index> order(parser, "O", "Interpolation order", {"interp-order"}, 3);
+  args::Flag clamp(parser, "C", "Clamp interpolation", {"interp-clamp"});
+  args::ValueFlag<Index> start2(parser, "S", "Second interp start", {"interp2"});
   args::ValueFlag<Index> nBasis(parser, "N", "Number of basis vectors to retain (overrides threshold)", {"nbasis"}, 5);
   args::Flag demean(parser, "C", "Mean-center dynamics", {"demean"});
   args::Flag rotate(parser, "V", "Rotate basis", {"rotate"});
@@ -37,13 +43,40 @@ int main_basis_img(args::Subparser &parser)
   Re3 const toMask = ref.chip<0>(0).abs();
   auto const toMaskMat = CollapseToArray(toMask);
   auto const [thresh, count] = Otsu(toMaskMat);
-  Index const fullN = shape[0] * spf.Get();
-  Eigen::MatrixXf dynamics(fullN, count);
+  Index const f = shape[0];
+  Index const s = shape[0] * spf.Get();
+  Eigen::MatrixXf dynamics(s, count);
   Index col = 0;
-  auto inds = Eigen::ArrayXi::LinSpaced(fullN, 0, shape[0] - 1);
+  Eigen::ArrayXi x1, x2, z1, z2;
+  Index f1, f2, s1, s2;
+  if (start2) {
+    f1 = start2.Get();
+    f2 = f - f1;
+    s1 = f1 * spf.Get();
+    s2 = f2 * spf.Get();
+    x1 = Eigen::ArrayXi::LinSpaced(f1, 0, s1 - 1) + spf.Get() / 2;
+    x2 = Eigen::ArrayXi::LinSpaced(f2, s1, s - 1) + spf.Get() / 2;
+    z1 = Eigen::ArrayXi::LinSpaced(s1, 0, s1 - 1);
+    z2 = Eigen::ArrayXi::LinSpaced(s2, s1, s - 1);
+  } else {
+    x1 = Eigen::ArrayXi::LinSpaced(f, 0, s - 1) + spf.Get() / 2;
+    z1 = Eigen::ArrayXi::LinSpaced(s, 0, s - 1);
+  }
   for (Index ii = 0; ii < realMat.cols(); ii++) {
     if (toMaskMat(ii) >= thresh) {
-      dynamics.col(col) = realMat.col(ii)(inds).transpose();
+      if (spf) {
+        if (start2) {
+          Interpolator interp1(x1, realMat.col(ii).head(f1), order.Get(), clamp);
+          Interpolator interp2(x2, realMat.col(ii).tail(f2), order.Get(), clamp);
+          dynamics.col(col).head(s1) = interp1(z1);
+          dynamics.col(col).tail(s2) = interp2(z2);
+        } else {
+          Interpolator interp(x1, realMat.col(ii), order.Get(), clamp);
+          dynamics.col(col) = interp(z1);
+        }
+      } else {
+        dynamics.col(col) = realMat.col(ii);
+      }
       col++;
     }
   }

src/interp.cpp

@@ -0,0 +1,69 @@
+#include "interp.hpp"
+#include "log.hpp"
+
+#include <set>
+
+namespace rl {
+
+Interpolator::Interpolator()
+  : m_min(0.)
+  , m_width(0.)
+{
+}
+
+// Sort input in ascending order
+
+Interpolator::Interpolator(Eigen::ArrayXi const &x, Eigen::ArrayXf const &y, Eigen::Index const order, bool const c)
+{
+  if (x.size() != y.size()) {
+    Log::Fail("Input vectors to spline must be same size");
+  }
+  if (x.size() == 0) {
+    Log::Fail("Cannot create a spline with no control points");
+  }
+  if (order < 1) {
+    Log::Fail("Interpolation order must be at least 1, was {}", order);
+  }
+  m_min = x[0];
+  m_width = x[x.size() - 1] - m_min;
+  m_clamp = c;
+  Eigen::ArrayXf const scaledx = (x.cast<float>() - m_min) / m_width;
+  m_spline = Eigen::SplineFitting<Spline>::Interpolate(y.transpose(), std::min<int>(x.rows() - 1, order), scaledx.transpose());
+}
+
+auto Interpolator::scale(float const x) const -> float { 
+    auto y = (x - m_min) / m_width;
+    if (m_clamp) {
+        return std::clamp(y, 0.f, 1.f);
+    } else {
+        return y;
+    }
+}
+
+auto Interpolator::scale(Eigen::ArrayXi const &x) const -> Eigen::ArrayXf { 
+    auto y = (x.cast<float>() - m_min) / m_width;
+    if (m_clamp) {
+        return y.cwiseMax(0.f).cwiseMin(1.f);
+    } else {
+        return y;
+    }
+}
+
+auto Interpolator::operator()(float const x) const -> float
+{
+  float const sx = scale(x);
+  float const val = m_spline(sx)[0];
+  return val;
+}
+
+auto Interpolator::operator()(const Eigen::ArrayXi &x) const -> Eigen::ArrayXf
+{
+  auto const sx = scale(x);
+  Eigen::ArrayXf output(sx.rows());
+  for (Eigen::Index i = 0; i < sx.rows(); i++) {
+    output[i] = m_spline(sx[i])[0];
+  }
+  return output;
+}
+
+} // namespace rl

src/interp.hpp

@@ -0,0 +1,33 @@
+#pragma once
+
+#include "Eigen/Core"
+#include <unsupported/Eigen/Splines>
+
+namespace rl {
+
+/*
+ * Eigen spline objects aren't very friendly. Wrap them in a class to do the required
+ * scaling and transposes to get them working.
+ */
+struct Interpolator
+{
+  using Spline = Eigen::Spline<float, 1>;
+
+  Interpolator();
+  Interpolator(
+    Eigen::ArrayXi const &x,
+    Eigen::ArrayXf const &y,
+    Eigen::Index const order = 3,
+    bool const clamp = false);
+  auto operator()(float const x) const -> float;
+  auto operator()(Eigen::ArrayXi const &x) const -> Eigen::ArrayXf;
+
+private:
+  Spline m_spline;
+  float m_min, m_width;
+  bool m_clamp;
+  auto scale(float const x) const -> float;
+  auto scale(Eigen::ArrayXi const &x) const -> Eigen::ArrayXf;
+};
+
+} // namespace rl