#ifndef _utl_SplineInterpolator_h_
#define _utl_SplineInterpolator_h_

#include <utl/UniformBSpline.h>
#include <utl/SplineFunction.h>
#include <utl/SplineSolver.h>

#include <vector>
#include <boost/multi_array.hpp>

#include <utl/AugerException.h>

namespace utl {

  namespace Spline {

    /**
      \brief Perform interpolation with B-splines on a regular grid.
      These are thin wrappers over the actual spline library.

      \author Hans Dembinski
      \date Nov 2010
      \ingroup math
    */

    using namespace ::Spline::Uniform;
    using namespace ::Spline;

    class Interpolator1D
    {
    public:

      Interpolator1D() {}

      Interpolator1D(const double xStart,
                     const double xStop, /* inclusive! */
                     const std::vector<double>& values,
                     const BoundaryCondition& bcX = BoundaryCondition());

      double
      operator()(const double x)
        const;

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

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

    protected:
      std::vector<double> fData;
      Function<1, KnotVector, BasisFunction> fSpline;
    };


    class Interpolator2D
    {
    public:

      Interpolator2D() {}

      Interpolator2D(const double xStart,
                     const double xStop, /* inclusive! */
                     const double yStart,
                     const double yStop, /* inclusive! */
                     const boost::multi_array<double,2>& values,
                     const BoundaryCondition bcX = BoundaryCondition(),
                     const BoundaryCondition bcY = BoundaryCondition());

      double
      operator()(const double x, const double y)
        const;

      inline
      double
      GetStart(const unsigned char dimension)
        const
      {
        return fSpline.GetStart(dimension);
      }

      inline
      double
      GetStop(const unsigned char dimension)
        const
      {
        return fSpline.GetStop(dimension);
      }

    protected:
      mutable double fX[2]; // proxy
      std::vector<double> fData;
      Function<2, KnotVector, BasisFunction> fSpline;
    };


    class Interpolator3D
    {
    public:

      Interpolator3D() {}

      Interpolator3D(const double xStart,
                     const double xStop, /* inclusive! */
                     const double yStart,
                     const double yStop, /* inclusive! */
                     const double zStart,
                     const double zStop, /* inclusive! */
                     const boost::multi_array<double,3>& values,
                     const BoundaryCondition bcX = BoundaryCondition(),
                     const BoundaryCondition bcY = BoundaryCondition(),
                     const BoundaryCondition bcZ = BoundaryCondition());

      double
      operator()(const double x, const double y, const double z)
        const;

      inline
      double
      GetStart(const unsigned char dimension)
        const
      {
        return fSpline.GetStart(dimension);
      }

      inline
      double
      GetStop(const unsigned char dimension)
        const
      {
        return fSpline.GetStop(dimension);
      }

    protected:
      mutable double fX[3]; // proxy
      std::vector<double> fData;
      Function<3, KnotVector, BasisFunction> fSpline;
    };


    class VectorInterpolator1D
    {
    public:

      VectorInterpolator1D() {}

      VectorInterpolator1D(const double xStart,
                           const double xStop, /* inclusive! */
                           const boost::multi_array<double,2>& fValue, // first index is function vector index!
                           const BoundaryCondition& bcX = BoundaryCondition());

      void
      operator()(std::vector<double>& result, const double x)
        const;

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

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

    protected:
      std::vector< std::vector<double> > fData;
      Function<1, KnotVector, BasisFunction> fSpline;
    };


    class VectorInterpolator2D
    {
    public:

      VectorInterpolator2D() {}

      VectorInterpolator2D(const double xStart,
                           const double xStop, /* inclusive! */
                           const double yStart,
                           const double yStop, /* inclusive! */
                           const boost::multi_array<double,3>& values, // first index is function vector index!
                           const BoundaryCondition& bcX = BoundaryCondition(),
                           const BoundaryCondition& bcY = BoundaryCondition());

      void
      operator()(std::vector<double>& result, const double x, const double y)
        const;

      inline
      double
      GetStart(const unsigned char dimension)
        const
      {
        return fSpline.GetStart(dimension);
      }

      inline
      double
      GetStop(const unsigned char dimension)
        const
      {
        return fSpline.GetStop(dimension);
      }

    protected:
      mutable double fX[2]; // proxy
      std::vector< std::vector<double> > fData;
      Function<2, KnotVector, BasisFunction> fSpline;
    };


    class VectorInterpolator3D
    {
    public:

      VectorInterpolator3D() {}

      VectorInterpolator3D(const double xStart,
                           const double xStop, /* inclusive! */
                           const double yStart,
                           const double yStop, /* inclusive! */
                           const double zStart,
                           const double zStop, /* inclusive! */
                           const boost::multi_array<double,4>& values, // first index is function vector index!
                           const BoundaryCondition& bcX = BoundaryCondition(),
                           const BoundaryCondition& bcY = BoundaryCondition(),
                           const BoundaryCondition& bcZ = BoundaryCondition());

      void
      operator()(std::vector<double>& result, const double x, const double y, const double z)
        const;

      inline
      double
      GetStart(const unsigned char dimension)
        const
      {
        return fSpline.GetStart(dimension);
      }

      inline
      double
      GetStop(const unsigned char dimension)
        const
      {
        return fSpline.GetStop(dimension);
      }

    protected:
      mutable double fX[3]; // proxy
      std::vector< std::vector<double> > fData;
      Function<3, KnotVector, BasisFunction> fSpline;
    };

  } // ns Spline

} // utl

#endif