UniformBSpline.h

Go to the documentation of this file.

#ifndef _UniformBSpline_Spline_Uniform_h_
#define _UniformBSpline_Spline_Uniform_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 Uniform {

    class KnotVector {
    public:
      KnotVector() { }

      KnotVector(const size_t n, const double start, const double stop) :
        fN(n),
        fStart(start),
        fDelta(stop-start)
      {
        if (n < 2)
          throw utl::InvalidConfigurationException("KnotVector: size < 2");

        if (stop < start)
          throw utl::InvalidConfigurationException("KnotVector: stop < start");
      }

      void
      PushBack(const double x)
      {
        if (fN == 0)
          fStart = x;
        fDelta = x - fStart;
        ++fN;
      }

      size_t
      Locate(const double x)
        const
      {
        if (x < fStart)
          return 0;

        if (x >= fStart + fDelta)
          return fN-2;

        return size_t((x - fStart)/fDelta*(fN-1));
      }

      double operator[](const size_t i) const
      { return fStart + double(i)/(fN-1)*fDelta; }

      size_t GetSize() const { return fN; }

      double GetStart() const { return fStart; }

      double GetStop() const { return fStart + fDelta; }

    protected:
      size_t fN = 0;
      double fStart = 0;
      double fDelta = 1;
    };


    class BasisFunction {
    public:
      BasisFunction() { }

      BasisFunction(const KnotVector* const xknot, const 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)]; }

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

        const double dx = xknot[1] - xknot[0];
        const double z = (x - xknot[i]) / dx;

        double t0, t1, t2, t3;
        switch (i - fIndex) {
        case 3:
          t0 = 1./6;
          t1 = -0.5;
          t2 = 0.5;
          t3 = -1./6;
          break;
        case 2:
          t0 = 2./3;
          t1 = 0;
          t2 = -1;
          t3 = 0.5;
          break;
        case 1:
          t0 = 1./6;
          t1 = 0.5;
          t2 = 0.5;
          t3 = -0.5;
          break;
        case 0:
          t0 = 0;
          t1 = 0;
          t2 = 0;
          t3 = 1./6;
          break;
        default:
          t0 = 0;
          t1 = 0;
          t2 = 0;
          t3 = 0;
        }

        switch (derivative) {
        case -1: {
          double result = z*(t0 + z*(t1/2 + z*(t2/3 + z*t3/4)))*dx;
          for (size_t j = std::max(fIndex, 0); j < i; ++j) {
            switch (j - fIndex) {
            case 3:
              t0 = 1./6;
              t1 = -0.5;
              t2 = 0.5;
              t3 = -1./6;
              break;
            case 2:
              t0 = 2./3;
              t1 = 0;
              t2 = -1;
              t3 = 0.5;
              break;
            case 1:
              t0 = 1./6;
              t1 = 0.5;
              t2 = 0.5;
              t3 = -0.5;
              break;
            case 0:
              t0 = 0;
              t1 = 0;
              t2 = 0;
              t3 = 1./6;
              break;
            default:
              t0 = 0;
              t1 = 0;
              t2 = 0;
              t3 = 0;
            }
            result += (t0 + (t1/2 + (t2/3 + t3/4)))*dx;
          }
          return result;
        }
        case 0:
          return t0 + z*(t1 + z*(t2 + z*t3));
        case 1:
          return (t1 + z*(2*t2 + z*3*t3))/dx;
        case 2:
          return (2*t2 + z*6*t3)/(dx*dx);
        case 3:
          return 6*t3/(dx*dx*dx);
        default:
          return 0;
        }
      }

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

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

  }

}


#endif