SplineFunction.h

Go to the documentation of this file.

#ifndef _SplineFunction_Spline_h_
#define _SplineFunction_Spline_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 {
  public:

    typedef std::vector<ABasisFunction> BasisVector;

    Function()
    {
      for (dim_t idim = 0; idim < ADimension; ++idim)
        fCacheX[idim] = std::numeric_limits<double>::quiet_NaN();
    }

    void
    Configure(const dim_t idim, const AKnotVector& knots)
    {
      fKnotVectors[idim] = knots;

      const size_t size = fKnotVectors[idim].GetSize();
      fBasisVectors[idim].resize(size+2);
      for (size_t k = 0; k < size+2; ++k)
        fBasisVectors[idim][k] = ABasisFunction(&fKnotVectors[idim], k);

      fStrides[0] = 0; // mark fStrides for re-computing
    }

    double
    operator()(const double* coefs, const double* x)
      const
    {
      if (fStrides[0] == 0)
      {
        // fStrides need recomputing
        fStrides[0] = 1;
        for (dim_t idim = 0; idim < ADimension-1; ++idim)
          fStrides[idim+1] = GetSize(idim)*fStrides[idim];
      }

      size_t i[ADimension];
      for (size_t idim = 0; idim < ADimension; ++idim)
        i[idim] = fKnotVectors[idim].Locate(x[idim]);

      double result = 0;
      for (size_t j = 0; j < (1 << 2*ADimension); ++j)
      {
        size_t k = 0;
        double basis = 1.0;
        for (dim_t idim = 0; idim < ADimension; ++idim)
        {
          const size_t kk = j % (1 << 2*(idim+1)) / (1 << 2*idim);
          if (x[idim] != fCacheX[idim])
            fCacheBasis[idim][kk] = fBasisVectors[idim][i[idim]+kk].Internal(i[idim], x[idim], 0);
          basis *= fCacheBasis[idim][kk];
          k += fStrides[idim]*(i[idim]+kk);
        }
        result += coefs[k]*basis;
      }

      for (size_t idim = 0; idim < ADimension; ++idim)
        fCacheX[idim] = x[idim];

      return result;
    }

    inline
    const AKnotVector&
    GetKnotVector(const dim_t idim)
      const
    {
      return fKnotVectors[idim];
    }

    inline
    const BasisVector&
    GetBasisVector(const dim_t idim)
      const
    {
      return fBasisVectors[idim];
    }

    inline
    size_t
    GetSize(const dim_t idim)
      const
    {
      return fBasisVectors[idim].size();
    }

    inline
    double
    GetStart(const dim_t idim)
      const
    {
      return fKnotVectors[idim].GetStart();
    }

    inline
    double
    GetStop(const dim_t idim)
      const
    {
      return fKnotVectors[idim].GetStop();
    }

  private:
    AKnotVector fKnotVectors[ADimension];
    BasisVector fBasisVectors[ADimension];
    mutable size_t fStrides[ADimension];
    mutable double fCacheX[ADimension];
    mutable double fCacheBasis[ADimension][4];
  };


  // specialization for ADimension = 1
  template<typename AKnotVector, typename ABasisFunction>
  class Function<1, AKnotVector, ABasisFunction>
  {
  public:

    typedef std::vector<ABasisFunction> BasisVector;
  
    Function() {}
  
    Function(const AKnotVector& knots)
    {
      Configure(0, knots);
    }

    void
    Configure(const dim_t /* idim */, const AKnotVector& knots)
    {
      fKnotVector = knots;

      const size_t size = fKnotVector.GetSize();
      fBasisVector.resize(size+2);
      for (size_t k = 0; k < size+2; ++k)
        fBasisVector[k] = ABasisFunction(&fKnotVector, k);
    }

    inline
    double
    operator()(const double* coefs, const double x, const char derivative=0)
      const
    {
      const size_t i = fKnotVector.Locate(x);
      const size_t kBegin = derivative >= 0 ? i : 0;
      const size_t kEnd = derivative >= 0 ? i+4 : GetSize();
      double result = 0;
      for (size_t k = kBegin; k < kEnd; ++k)
        result += coefs[k]*fBasisVector[k].Internal(i,x,derivative);
      return result;
    }

    inline
    const ABasisFunction&
    operator[](const size_t k)
      const
    {
      return fBasisVector[k];
    }

    inline
    const AKnotVector&
    GetKnotVector()
      const
    {
      return fKnotVector;
    }

    inline
    const AKnotVector&
    GetKnotVector(const dim_t)
      const
    {
      return fKnotVector;
    }

    inline
    const BasisVector&
    GetBasisVector()
      const
    {
      return fBasisVector;
    }

    inline
    const BasisVector&
    GetBasisVector(const dim_t)
      const
    {
      return fBasisVector;
    }

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

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

    inline
    double
    GetStart()
      const
    {
      return fKnotVector.GetStart();
    }

    inline
    double
    GetStart(const dim_t)
      const
    {
      return fKnotVector.GetStart();
    }

    inline
    double
    GetStop()
      const
    {
      return fKnotVector.GetStop();
    }

    inline
    double
    GetStop(const dim_t)
      const
    {
      return fKnotVector.GetStop();
    }

  private:
    AKnotVector fKnotVector;
    BasisVector fBasisVector;
  };

} // NS Spline

#endif