BSplineGenerator.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/>
 * ########################################################################
 *
 * [1] https://en.wikipedia.org/wiki/B-spline
 */
 
#ifndef BSPLINE_BSPLINEGENERATOR_H
#define BSPLINE_BSPLINEGENERATOR_H
 
#include <bspline/Spline.h>
#include <bspline/exceptions/BSplineException.h>
#include <bspline/operators/CompoundOperators.h>
#include <bspline/operators/Position.h>
#include <bspline/operators/ScalarOperators.h>
 
#include <algorithm>
 
namespace bspline {
using namespace bspline::exceptions;
 
template <typename T>
class BSplineGenerator final {
 private:
  Grid<T> _grid;
 
  std::vector<T> _knots;
 
  Grid<T> generateGrid(std::vector<T> knots) {
    const auto endIterator = std::unique(knots.begin(), knots.end());
    knots.erase(endIterator, knots.end());
    return Grid<T>(std::move(knots));
  }
 
 public:
  explicit BSplineGenerator(std::vector<T> knots)
      : _grid(generateGrid(knots)), _knots(std::move(knots)){};
 
  BSplineGenerator(std::vector<T> knots, const Grid<T> &grid)
      : _grid(grid), _knots(std::move(knots)) {
    // Check, whether the given knots vector and grid are consistent.
    const Grid<T> secondGrid = generateGrid(_knots);
    if (_grid != secondGrid) {
      throw BSplineException(ErrorCode::INCONSISTENT_DATA);
    }
  };
 
  Grid<T> getGrid() const { return _grid; };
 
  template <size_t order>
  std::vector<Spline<T, order>> generateBSplines() const {
    static constexpr size_t k = order + 1;
    if (_knots.size() < k) {
      throw BSplineException(ErrorCode::UNDETERMINED,
                             "The knots vector contains too few elements to "
                             "generate BSplines of the requested order.");
    }
 
    if constexpr (order == 0) {
      return generateZerothOrderSplines();
    } else {
      const auto nextLowerOrderSplines = generateBSplines<order - 1>();
      std::vector<Spline<T, order>> ret;
      ret.reserve(_knots.size() - k);
      for (size_t i = 0; i < _knots.size() - k; i++) {
        ret.push_back(applyRecursionRelation<k>(
            i, nextLowerOrderSplines.at(i), nextLowerOrderSplines.at(i + 1)));
      }
      return ret;
    }
  }
 
 private:
  template <size_t k>
  Spline<T, k - 1> applyRecursionRelation(
      size_t i, const Spline<T, k - 2> &splinei,
      const Spline<T, k - 2> &splineip1) const {
    static_assert(k >= 2, "k has to be at least 2.");
 
    Spline<T, k - 1> ret(_grid);
 
    const T &xi = _knots.at(i);
    const T &xipkm1 = _knots.at(i + k - 1);
    if (xipkm1 > xi) {
      const T prefac = static_cast<T>(1) / (xipkm1 - xi);
      const auto op = prefac * (operators::X<1>{} - xi);
      ret = op * splinei;
    }
 
    const T &xip1 = _knots.at(i + 1);
    const T &xipk = _knots.at(i + k);
    if (xipk > xip1) {
      const T prefac = static_cast<T>(1) / (xipk - xip1);
      const auto op = prefac * (xipk - operators::X<1>{});
      ret += op * splineip1;
    }
 
    return ret;
  }
 
  std::vector<Spline<T, 0>> generateZerothOrderSplines() const {
    std::vector<Spline<T, 0>> ret;
    const size_t numberOfSplines = (_knots.empty()) ? 0 : _knots.size() - 1;
    ret.reserve(numberOfSplines);
 
    for (size_t i = 0; i < numberOfSplines; i++) {
      const T &xi = _knots.at(i);
      const T &xip1 = _knots.at(i + 1);
 
      if (xi > xip1) {
        throw BSplineException(ErrorCode::UNDETERMINED);
      } else if (xi == xip1) {
        ret.push_back(Spline<T, 0>{_grid});
      } else {
        std::vector<std::array<T, 1>> coefficients{{static_cast<T>(1)}};
 
        const size_t gridIndex = _grid.findElement(xi);
 
        ret.push_back(Spline<T, 0>{Support(_grid, gridIndex, gridIndex + 2),
                                   std::move(coefficients)});
      }
    }
    return ret;
  }
};
 
template <size_t order, typename T>
std::vector<Spline<T, order>> generateBSplines(std::vector<T> knots) {
  BSplineGenerator gen(std::move(knots));
  return gen.template generateBSplines<order>();
}
 
}  // namespace bspline
#endif  // BSPLINE_BSPLINEGENERATOR_H