SplineOperator.h

/*
 * ########################################################################
 * The contents of this file is free and unencumbered software released into the
 * public domain. For more information, please refer to <http://unlicense.org/>
 * ########################################################################
 */
 
#ifndef BSPLINE_OPERATORS_SPLINEOPERATOR_H
#define BSPLINE_OPERATORS_SPLINEOPERATOR_H
 
#include <bspline/Spline.h>
#include <bspline/exceptions/BSplineException.h>
#include <bspline/internal/misc.h>
#include <bspline/operators/GenericOperators.h>
 
namespace bspline::operators {
 
template <typename T, size_t order>
class SplineOperator final : public Operator {
 private:
  Spline<T, order> _s;
 
 public:
  SplineOperator(Spline<T, order> s) : _s(std::move(s)){};
 
  static constexpr size_t outputOrder(size_t inputOrder) {
    return inputOrder + order;
  }
 
  template <size_t size>
  auto transform(const std::array<T, size> &input, const support::Grid<T> &grid,
                 size_t intervalIndex) const {
    static_assert(size >= 1);
    constexpr size_t OUTPUT_SIZE = outputOrder(size - 1) + 1;
 
    if (_s.getSupport().getGrid() != grid) {
      throw exceptions::BSplineException(ErrorCode::DIFFERING_GRIDS);
    }
 
    const auto relativeIndex =
        _s.getSupport().relativeFromAbsolute(intervalIndex);
 
    auto retVal = internal::make_array<T, OUTPUT_SIZE>(static_cast<T>(0));
 
    if (relativeIndex) {
      // The interval is part of the Spline's support.
      const auto &coeffs = _s.getCoefficients()[*relativeIndex];
 
      for (size_t i = 0; i < input.size(); i++) {
        for (size_t j = 0; j < coeffs.size(); j++) {
          retVal[i + j] += input[i] * coeffs[j];
        }
      }
    }
    return retVal;
  }
};
}  // namespace bspline::operators
#endif  // BSPLINE_OPERATORS_SPLINEOPERATOR_H