bspline.h

Go to the documentation of this file.

 
#ifndef __SABLIB_BSPLINE_H__
#define __SABLIB_BSPLINE_H__
 
#include <algorithm>
#include <type_traits>
#include <vector>
#include <Eigen/Eigen>
 
namespace sablib {

template <typename Scalar>
class BSpline final
{
    static_assert(std::is_floating_point_v<Scalar>, "Scalar must be a floating-point type.");
 
public:
    BSpline() = delete;

    BSpline(const int degree, const Eigen::VectorX<Scalar> & knots);

    BSpline(
        const int degree, const Eigen::VectorX<Scalar> & knots,
        const Eigen::VectorX<Scalar> & x, const Eigen::VectorX<Scalar> & y
    );

    BSpline(
        const int degree, const Eigen::VectorX<Scalar> & knots, const Eigen::VectorX<Scalar> & coefficients
    );

    int BasisSize() const;

    const Eigen::VectorX<Scalar> Knots() const;

    const Eigen::VectorX<Scalar> Coefficients() const;

    const Eigen::SparseMatrix<Scalar> DesignMatrix(const Eigen::VectorX<Scalar> & x) const;

    void Fit(const Eigen::VectorX<Scalar> & x, const Eigen::VectorX<Scalar> & y);

    const Scalar Interpolate(const Scalar x, const Eigen::VectorX<Scalar> & coefficients) const;

    const Scalar Interpolate(const Scalar x) const;
 
private:
    int sp_degree;                          
    Eigen::VectorX<Scalar> sp_knots;        
    int basis;                              
    Eigen::VectorX<Scalar> sp_coefficients; 

    int FindSpan(const Scalar x) const;

    const Eigen::VectorX<Scalar> BasisFunctions(const int span, const Scalar x) const;
};
 
//
// Implementation of constructor
//
template <typename Scalar>
BSpline<Scalar>::BSpline(
    const int degree, const Eigen::VectorX<Scalar> & knots
) : sp_degree(degree), sp_knots(knots), basis(knots.size() - degree - 1)
{
}
 
//
// Implementation of constructor with fitting
//
template <typename Scalar>
BSpline<Scalar>::BSpline(
    const int degree, const Eigen::VectorX<Scalar> & knots,
    const Eigen::VectorX<Scalar> & x, const Eigen::VectorX<Scalar> & y
) : sp_degree(degree), sp_knots(knots), basis(knots.size() - degree - 1)
{
    Fit(x, y);
}
 
//
// Implementation of constructor with coefficients
//
template <typename Scalar>
BSpline<Scalar>::BSpline(
    const int degree, const Eigen::VectorX<Scalar> & knots, const Eigen::VectorX<Scalar> & coefficients
) : sp_degree(degree), sp_knots(knots), basis(knots.size() - degree - 1), sp_coefficients(coefficients)
{
}
 
//
// Implementation of BasisSize() method
//
template <typename Scalar>
int BSpline<Scalar>::BasisSize() const
{
    return basis;
}
 
//
// Implementation of Knots() method
//
template <typename Scalar>
const Eigen::VectorX<Scalar> BSpline<Scalar>::Knots() const
{
    return sp_knots;
}
 
//
// Implementation of Coefficients() method
//
template <typename Scalar>
const Eigen::VectorX<Scalar> BSpline<Scalar>::Coefficients() const
{
    return sp_coefficients;
}
 
//
// Implementation of FindSpan() method
//
template <typename Scalar>
int BSpline<Scalar>::FindSpan(const Scalar x) const
{
    int n = basis - 1;
 
    if (x >= sp_knots(n + 1)) {
        return n;
    }
 
    if (x <= sp_knots(sp_degree)) {
        return sp_degree;
    }
 
    int low  = sp_degree;
    int high = n + 1;
    int mid  = (low + high) / 2;
 
    while (x < sp_knots(mid) || x >= sp_knots(mid + 1)) {
        if (x < sp_knots(mid)) {
            high = mid;
        }
        else {
            low = mid;
        }
 
        mid = (low + high) / 2;
    }
 
    return mid;
}
 
//
// Implementation of BasisFunctions() method
//
template <typename Scalar>
const Eigen::VectorX<Scalar> BSpline<Scalar>::BasisFunctions(const int span, const Scalar x) const
{
    Eigen::VectorX<Scalar> N(sp_degree + 1);
    Eigen::VectorX<Scalar> left(sp_degree + 1);
    Eigen::VectorX<Scalar> right(sp_degree + 1);
 
    N(0) = 1.0;
 
    for(int j = 1; j <= sp_degree; j++) {
        left(j)  = x - sp_knots(span + 1 - j);
        right(j) = sp_knots(span + j) - x;
 
        Scalar saved = 0.0;
 
        for(int r = 0; r < j; r++) {
            Scalar temp = N(r) / (right(r + 1) + left(j - r));
            N(r) = saved + right(r + 1) * temp;
            saved = left(j - r) * temp;
        }
 
        N(j) = saved;
    }
 
    return N;
}
 
//
// Implementation of DesignMatrix() method
//
template <typename Scalar>
const Eigen::SparseMatrix<Scalar> BSpline<Scalar>::DesignMatrix(const Eigen::VectorX<Scalar> & x) const
{
    int n = x.size();
    std::vector< Eigen::Triplet<Scalar> > triplets;
 
    for(int r = 0; r < n; r++) {
        Scalar xi = x(r);
        int span = FindSpan(xi);
        Eigen::VectorX<Scalar> N = BasisFunctions(span, xi);
 
        for(int j = 0; j <= sp_degree; j++) {
            int col = span - sp_degree + j;
 
            if(col >=0 && col < basis) {
                triplets.emplace_back(r, col, N(j));
            }
        }
    }
 
    Eigen::SparseMatrix<Scalar> B(n , basis);
    B.setFromTriplets(triplets.begin(), triplets.end());
 
    return B;
}
 
//
// Implementation of Fit() method
//
template <typename Scalar>
void BSpline<Scalar>::Fit(const Eigen::VectorX<Scalar> & x, const Eigen::VectorX<Scalar> & y)
{
    auto B = DesignMatrix(x);
    Eigen::SparseMatrix<Scalar> BTB = B.transpose() * B;
    Eigen::SimplicialLDLT< Eigen::SparseMatrix<Scalar> > solver;
 
    solver.compute(BTB);
 
    if(solver.info() != Eigen::Success) {
        throw std::runtime_error("BSpline::Fit(): solver calculation fails.");
    }
 
    sp_coefficients = solver.solve(B.transpose() * y);
}
 
//
// Implementation of Interpolate() method
//
template<typename Scalar>
const Scalar BSpline<Scalar>::Interpolate(const Scalar x, const Eigen::VectorX<Scalar> & coefficients) const
{
    int span = FindSpan(x);
    Eigen::VectorX<Scalar> N = BasisFunctions(span, x);
 
    Scalar y = 0.0;
 
    for(int j = 0; j <= sp_degree; j++) {
        int idx = span - sp_degree + j;
 
        if(idx >=0 && idx < coefficients.size()) {
            y += N(j) * coefficients(idx);
        }
    }
 
    return y;
}
 
//
// Implementation of Interpolate() method with internal coefficients
//
template<typename Scalar>
inline const Scalar BSpline<Scalar>::Interpolate(const Scalar x) const
{
    return Interpolate(x, sp_coefficients);
}
 
} // namespace sablib
 
#endif // __SABLIB_BSPLINE_H__