GeneralBSpline.h

Go to the documentation of this file.

#ifndef _GeneralBSpline_Spline_General_h_
#define _GeneralBSpline_Spline_General_h_

#include <cmath>
#include <vector>
#include <algorithm>
#include <limits>
#include <iostream>
#include <utl/AugerException.h>

namespace Spline {

  typedef unsigned char dim_t;

  template <dim_t ADimension, typename AKnotVector, typename ABasisFunction>
  class Function;

  namespace General {

    class KnotVector {
    public:
      KnotVector() {}

      KnotVector(const std::vector<double>& data)
        : fData(data)
      {
        const size_t n = GetSize();
        if (n<2)
          throw utl::InvalidConfigurationException("size of KnotVector < 2");
        for (size_t i=0;i<n-1;++i)
          if (fData[i]>fData[i+1])
            throw utl::InvalidConfigurationException("knots are not in ascending order");
      }

      inline
      void
      PushBack(const double x)
      {
        fData.push_back(x);
      }

      inline
      double
      operator()(const int i)
        const
      {
        const size_t n = GetSize();

        // virtual knots
        if (i<0)
        {
          const double dx = fData[1]-fData[0];
          return fData[0]+i*dx;
        }

        if (i>int(n-1))
        {
          const double dx = fData[n-1]-fData[n-2];
          return fData[n-1]+(i-(n-1))*dx;
        }

        // real knots
        return fData[i];
      }

      inline
      double
      operator[](const size_t i)
        const
      {
        return fData[i];
      }

      inline
      size_t
      Locate(const double x)
        const
      {
        const size_t n = GetSize();
        if (x >= (*this)[n-1])
          return n-2;
        else
          return std::upper_bound(fData.begin(),fData.end(),x)-fData.begin()-1;
      }

      inline
      size_t
      GetSize()
        const
      {
        return fData.size();
      }

      inline
      double
      GetStart()
        const
      {
        return fData.front();
      }

      inline
      double
      GetStop()
        const
      {
        return fData.back();
      }

    protected:
      std::vector<double> fData;
    };


    class BasisFunction {
    public:
      BasisFunction() { }

      BasisFunction(const KnotVector* const xknot, int j)
        : fKnotPtr(xknot), fIndex(j-3) { }

      BasisFunction(const BasisFunction& other) { *this = other; }

      BasisFunction&
      operator=(const BasisFunction& other)
      {
        fKnotPtr = other.fKnotPtr;
        fIndex = other.fIndex;
        return *this;
      }

      double
      operator()(const double x, const char derivative = 0)
        const
      {
        if (x < GetStart())
          return 0;

        if (x > GetStop())
          return derivative == -1 ? Internal(fKnotPtr->Locate(x), GetStop(), -1) : 0;

        return Internal(fKnotPtr->Locate(x), x, derivative);
      }

      double GetStart() const { return (*fKnotPtr)[std::max(fIndex, 0)]; }

      double GetStop() const { return (*fKnotPtr)[std::min(fIndex+4, int(fKnotPtr->GetSize())-1)]; }

    private:
      double
      Internal(const size_t i, const double x, const char derivative)
        const
      {
        const KnotVector& xknot = *fKnotPtr;

        double t0;
        double t1;
        double t2;
        double t3;
        ComputePolynomialFactors(t0, t1, t2, t3, i - fIndex);

        switch (derivative) {
        case -1: {
          const double a = xknot[i];
          double result =
            x*(t0 + x*(t1/2 + x*(t2/3 + x*t3/4))) -
            a*(t0 + a*(t1/2 + a*(t2/3 + a*t3/4)));
          for (size_t j = std::max(fIndex, 0); j < i; ++j) {
            ComputePolynomialFactors(t0, t1, t2, t3, j - fIndex);
            const double pa = xknot[j];
            const double pb = xknot[j+1];
            result +=
              pb*(t0 + pb*(t1/2 + pb*(t2/3 + pb*t3/4))) -
              pa*(t0 + pa*(t1/2 + pa*(t2/3 + pa*t3/4)));
          }
          return result;
        }
        case 0:
          return t0 + x*(t1 + x*(t2 + x*t3));
        case 1:
          return t1 + x*(2*t2 + x*3*t3);
        case 2:
          return 2*t2 + x*6*t3;
        case 3:
          return 6*t3;
        default:
          return 0;
        }
      }

      void
      ComputePolynomialFactors(double& t0, double& t1, double& t2, double& t3, const size_t i)
        const
      {
        const KnotVector& xknot = *fKnotPtr;

        t0 = 0;
        t1 = 0;
        t2 = 0;
        t3 = 0;

        #define SPLINE_BASIS_POLYNOMAL_FACTORS(ia0, ia1, ia2, ib0, ib1, ib2, ic0, ic1, ic2) \
          t0 += -xknot(fIndex+ia0)*xknot(fIndex+ia1)*xknot(fIndex+ia2)/((xknot(fIndex+ib0)-xknot(fIndex+ic0))*(xknot(fIndex+ib1)-xknot(fIndex+ic1))*(xknot(fIndex+ib2)-xknot(fIndex+ic2))); \
          t1 += (xknot(fIndex+ia0)*xknot(fIndex+ia1)+xknot(fIndex+ia0)*xknot(fIndex+ia2)+xknot(fIndex+ia1)*xknot(fIndex+ia2))/((xknot(fIndex+ib0)-xknot(fIndex+ic0))*(xknot(fIndex+ib1)-xknot(fIndex+ic1))*(xknot(fIndex+ib2)-xknot(fIndex+ic2))); \
          t2 += -(xknot(fIndex+ia0)+xknot(fIndex+ia1)+xknot(fIndex+ia2))/((xknot(fIndex+ib0)-xknot(fIndex+ic0))*(xknot(fIndex+ib1)-xknot(fIndex+ic1))*(xknot(fIndex+ib2)-xknot(fIndex+ic2))); \
          t3 += 1.0/((xknot(fIndex+ib0)-xknot(fIndex+ic0))*(xknot(fIndex+ib1)-xknot(fIndex+ic1))*(xknot(fIndex+ib2)-xknot(fIndex+ic2)));

        switch (i) {
        case 0:
          SPLINE_BASIS_POLYNOMAL_FACTORS(0,0,0,1,2,3,0,0,0);
          break;
        case 1:
          SPLINE_BASIS_POLYNOMAL_FACTORS(0,0,2,2,3,1,0,0,2);
          SPLINE_BASIS_POLYNOMAL_FACTORS(0,1,3,3,2,1,0,1,3);
          SPLINE_BASIS_POLYNOMAL_FACTORS(1,1,4,2,3,1,1,1,4);
          break;
        case 2:
          SPLINE_BASIS_POLYNOMAL_FACTORS(0,3,3,3,1,2,0,3,3);
          SPLINE_BASIS_POLYNOMAL_FACTORS(1,4,3,3,1,2,1,4,3);
          SPLINE_BASIS_POLYNOMAL_FACTORS(2,4,4,3,1,2,2,4,4);
          break;
        case 3:
          SPLINE_BASIS_POLYNOMAL_FACTORS(4,4,4,1,2,3,4,4,4);
          break;
        }

        #undef SPLINE_BASIS_POLYNOMAL_FACTORS
      }

      const KnotVector* fKnotPtr = nullptr;
      int fIndex = 0;

      template <dim_t ADimension, typename AKnotVector, typename ABasisFunction>
      friend class ::Spline::Function;
    };

  }

}


#endif