Integrated spline class and simple spline fitting

8 files changed, 895 insertions, 2 deletions

diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt
index 679c96736..e961e72c5 100644
--- a/unsupported/Eigen/CMakeLists.txt
+++ b/unsupported/Eigen/CMakeLists.txt
@@ -1,6 +1,6 @@
 set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials
                   FFT NonLinearOptimization SparseExtra IterativeSolvers
-                  NumericalDiff Skyline MPRealSupport OpenGLSupport KroneckerProduct
+                  NumericalDiff Skyline MPRealSupport OpenGLSupport KroneckerProduct Splines
    )
 
 install(FILES
diff --git a/unsupported/Eigen/Splines b/unsupported/Eigen/Splines
new file mode 100644
index 000000000..362274157
--- /dev/null
+++ b/unsupported/Eigen/Splines
@@ -0,0 +1,52 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SPLINES_MODULE_H
+#define EIGEN_SPLINES_MODULE_H
+
+namespace Eigen 
+{
+/** \ingroup Unsupported_modules
+  * \defgroup Splines_Module Spline and spline fitting module
+  *
+  * This module provides a simple multi-dimensional spline class while
+  * offering most basic functionality to fit a spline to point sets.
+  *
+  * \code
+  * #include <unsupported/Eigen/Splines>
+  * \endcode
+  */
+//@{
+}
+
+#include "src/Splines/SplineFwd.h"
+#include "src/Splines/Spline.h"
+#include "src/Splines/SplineFitting.h"
+
+namespace Eigen 
+{
+//@}
+}
+
+#endif // EIGEN_SPLINES_MODULE_H
diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt
index 6fd98b310..f3180b52b 100644
--- a/unsupported/Eigen/src/CMakeLists.txt
+++ b/unsupported/Eigen/src/CMakeLists.txt
@@ -9,4 +9,5 @@ ADD_SUBDIRECTORY(NumericalDiff)
 ADD_SUBDIRECTORY(Polynomials)
 ADD_SUBDIRECTORY(Skyline)
 ADD_SUBDIRECTORY(SparseExtra)
-ADD_SUBDIRECTORY(KroneckerProduct)
\ No newline at end of file
+ADD_SUBDIRECTORY(KroneckerProduct)
+ADD_SUBDIRECTORY(Splines)
diff --git a/unsupported/Eigen/src/Splines/Spline.h b/unsupported/Eigen/src/Splines/Spline.h
new file mode 100644
index 000000000..bb4b73e5a
--- /dev/null
+++ b/unsupported/Eigen/src/Splines/Spline.h
@@ -0,0 +1,407 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SPLINE_H
+#define EIGEN_SPLINE_H
+
+#include "SplineFwd.h"
+
+namespace Eigen
+{
+  /**
+  * \class Spline class
+  * \brief A class representing N-D spline curves.
+  * \tparam _Scalar The underlying data type (typically float or double)
+  * \tparam _Dim The curve dimension (e.g. 2 or 3)
+  * \tparam _Degree Per default set to Dynamic; could be set to the actual desired
+  *                 degree for optimization purposes (would result in stack allocation
+  *                 of several temporary variables).
+  **/
+  template <typename _Scalar, int _Dim, int _Degree>
+  class Spline
+  {
+  public:
+    typedef _Scalar Scalar;
+    enum { Dimension = _Dim };
+    enum { Degree = _Degree };
+
+    typedef typename SplineTraits<Spline>::PointType PointType;
+    typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;
+    typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
+    typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
+
+    /**
+    * \brief Creates a spline from a knot vector and control points.
+    **/
+    template <typename OtherVectorType, typename OtherArrayType>
+    Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {}
+
+    template <int OtherDegree>
+    Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) : 
+    m_knots(spline.knots()), m_ctrls(spline.ctrls()) {}
+
+    /* Const access methods for knots and control points. */
+    const KnotVectorType& knots() const { return m_knots; }
+    const ControlPointVectorType& ctrls() const { return m_ctrls; }
+
+    /* Spline evaluation. */
+    PointType operator()(Scalar u) const;
+
+    /* Evaluation of spline derivatives of up-to given order. 
+    * The returned matrix has dimensions Dim-by-(Order+1) containing
+    * the 0-th order up-to Order-th order derivatives.
+    */
+    typename SplineTraits<Spline>::DerivativeType
+      derivatives(Scalar u, DenseIndex order) const;
+
+    /**
+    * Evaluation of spline derivatives of up-to given order.
+    * The function performs identically to derivatives(Scalar, int) but
+    * does not require any heap allocations.
+    * \sa derivatives(Scalar, int)
+    **/    
+    template <int DerivativeOrder>
+    typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
+      derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
+
+    /* Non-zero spline basis functions. */
+    typename SplineTraits<Spline>::BasisVectorType
+      basisFunctions(Scalar u) const;
+
+    /* Non-zero spline basis function derivatives up to given order. 
+    * The order is different from the spline order - it is the order
+    * up to which derivatives will be computed.
+    * \sa basisFunctions(Scalar)
+    */
+    typename SplineTraits<Spline>::BasisDerivativeType
+      basisFunctionDerivatives(Scalar u, DenseIndex order) const;
+
+    /**
+    * Computes non-zero basis function derivatives up to the given derivative order.
+    * As opposed to basisFunctionDerivatives(Scalar, int) this function does not perform
+    * any heap allocations.
+    * \sa basisFunctionDerivatives(Scalar, int)
+    **/    
+    template <int DerivativeOrder>
+    typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
+      basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
+
+    /**
+    * \brief The current spline degree. It's a function of knot size and number 
+    * of controls and thus does not require a dedicated member. 
+    */ 
+    DenseIndex degree() const;
+
+    /** Computes the span within the knot vector in which u falls. */
+    DenseIndex span(Scalar u) const;
+
+    static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots);
+    static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
+
+
+  private:
+    KnotVectorType m_knots; /* Knot vector. */
+    ControlPointVectorType  m_ctrls; /* Control points. */
+  };
+
+  template <typename _Scalar, int _Dim, int _Degree>
+  DenseIndex Spline<_Scalar, _Dim, _Degree>::Span(
+    typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::Scalar u,
+    DenseIndex degree,
+    const typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::KnotVectorType& knots)
+  {
+    // Piegl & Tiller, "The NURBS Book", A2.1 (p. 68)
+    if (u <= knots(0)) return degree;
+    const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
+    return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 );
+  }
+
+  template <typename _Scalar, int _Dim, int _Degree>
+  typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType
+    Spline<_Scalar, _Dim, _Degree>::BasisFunctions(
+    typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
+    DenseIndex degree,
+    const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
+  {
+    typedef typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType BasisVectorType;
+
+    const DenseIndex p = degree;
+    const DenseIndex i = Spline::Span(u, degree, knots);
+
+    const KnotVectorType& U = knots;
+
+    BasisVectorType left(p+1); left(0) = Scalar(0);
+    BasisVectorType right(p+1); right(0) = Scalar(0);        
+
+    VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse();
+    VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u;
+
+    BasisVectorType N(1,p+1);
+    N(0) = Scalar(1);
+    for (DenseIndex j=1; j<=p; ++j)
+    {
+      Scalar saved = Scalar(0);
+      for (DenseIndex r=0; r<j; r++)
+      {
+        const Scalar tmp = N(r)/(right(r+1)+left(j-r));
+        N[r] = saved + right(r+1)*tmp;
+        saved = left(j-r)*tmp;
+      }
+      N(j) = saved;
+    }
+    return N;
+  }
+
+  template <typename _Scalar, int _Dim, int _Degree>
+  DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const
+  {
+    if (_Degree == Dynamic)
+      return m_knots.size() - m_ctrls.cols() - 1;
+    else
+      return _Degree;
+  }
+
+  template <typename _Scalar, int _Dim, int _Degree>
+  DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const
+  {
+    return Spline::Span(u, degree(), knots());
+  }
+
+  /**
+  * \brief A functor for the computation of a spline point.
+  * \sa Piegl & Tiller, "The NURBS Book", A4.1 (p. 124)
+  **/
+  template <typename _Scalar, int _Dim, int _Degree>
+  typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const
+  {
+    enum { Order = SplineTraits<Spline>::OrderAtCompileTime };
+
+    const DenseIndex span = this->span(u);
+    const DenseIndex p = degree();
+    const BasisVectorType basis_funcs = basisFunctions(u);
+
+    const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
+    const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
+    return (ctrl_weights * ctrl_pts).rowwise().sum();
+  }
+
+  /* --------------------------------------------------------------------------------------------- */
+
+  template <typename SplineType, typename DerivativeType>
+  void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
+  {    
+    enum { Dimension = SplineTraits<SplineType>::Dimension };
+    enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
+    enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
+
+    typedef typename SplineTraits<SplineType>::Scalar Scalar;
+
+    typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
+    typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
+
+    typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
+    typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;    
+
+    const DenseIndex p = spline.degree();
+    const DenseIndex span = spline.span(u);
+
+    const DenseIndex n = (std::min)(p, order);
+
+    der.resize(Dimension,n+1);
+
+    // Retrieve the basis function derivatives up to the desired order...    
+    const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);
+
+    // ... and perform the linear combinations of the control points.
+    for (DenseIndex der_order=0; der_order<n+1; ++der_order)
+    {
+      const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
+      const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
+      der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
+    }
+  }
+
+  /**
+  * \brief A functor for the computation of a spline point.
+  * \sa Piegl & Tiller, "The NURBS Book", A4.1 (p. 124)
+  **/
+  template <typename _Scalar, int _Dim, int _Degree>
+  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType
+    Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
+  {
+    typename SplineTraits< Spline >::DerivativeType res;
+    derivativesImpl(*this, u, order, res);
+    return res;
+  }
+
+  template <typename _Scalar, int _Dim, int _Degree>
+  template <int DerivativeOrder>
+  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType
+    Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
+  {
+    typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
+    derivativesImpl(*this, u, order, res);
+    return res;
+  }
+
+  /**
+  * \brief A functor for the computation of a spline's non-zero basis functions.
+  * \sa Piegl & Tiller, "The NURBS Book", A2.2 (p. 70)
+  **/
+  template <typename _Scalar, int _Dim, int _Degree>
+  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType
+    Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const
+  {
+    return Spline::BasisFunctions(u, degree(), knots());
+  }
+
+  /* --------------------------------------------------------------------------------------------- */
+
+  template <typename SplineType, typename DerivativeType>
+  void basisFunctionDerivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_)
+  {
+    enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
+
+    typedef typename SplineTraits<SplineType>::Scalar Scalar;
+    typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
+    typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType;
+    typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
+
+    const KnotVectorType& U = spline.knots();
+
+    const DenseIndex p = spline.degree();
+    const DenseIndex span = spline.span(u);
+
+    const DenseIndex n = (std::min)(p, order);
+
+    N_.resize(n+1, p+1);
+
+    BasisVectorType left = BasisVectorType::Zero(p+1);
+    BasisVectorType right = BasisVectorType::Zero(p+1);
+
+    Matrix<Scalar,Order,Order> ndu(p+1,p+1);
+
+    double saved, temp;
+
+    ndu(0,0) = 1.0;
+
+    DenseIndex j;
+    for (j=1; j<=p; ++j)
+    {
+      left[j] = u-U[span+1-j];
+      right[j] = U[span+j]-u;
+      saved = 0.0;
+
+      for (DenseIndex r=0; r<j; ++r)
+      {
+        /* Lower triangle */
+        ndu(j,r) = right[r+1]+left[j-r];
+        temp = ndu(r,j-1)/ndu(j,r);
+        /* Upper triangle */
+        ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp);
+        saved = left[j-r] * temp;
+      }
+
+      ndu(j,j) = static_cast<Scalar>(saved);
+    }
+
+    for (j = p; j>=0; --j) 
+      N_(0,j) = ndu(j,p);
+
+    // Compute the derivatives
+    DerivativeType a(n+1,p+1);
+    DenseIndex r=0;
+    for (; r<=p; ++r)
+    {
+      DenseIndex s1,s2;
+      s1 = 0; s2 = 1; // alternate rows in array a
+      a(0,0) = 1.0;
+
+      // Compute the k-th derivative
+      for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
+      {
+        double d = 0.0;
+        DenseIndex rk,pk,j1,j2;
+        rk = r-k; pk = p-k;
+
+        if (r>=k)
+        {
+          a(s2,0) = a(s1,0)/ndu(pk+1,rk);
+          d = a(s2,0)*ndu(rk,pk);
+        }
+
+        if (rk>=-1) j1 = 1;
+        else        j1 = -rk;
+
+        if (r-1 <= pk) j2 = k-1;
+        else           j2 = p-r;
+
+        for (j=j1; j<=j2; ++j)
+        {
+          a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
+          d += a(s2,j)*ndu(rk+j,pk);
+        }
+
+        if (r<=pk)
+        {
+          a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
+          d += a(s2,k)*ndu(r,pk);
+        }
+
+        N_(k,r) = static_cast<Scalar>(d);
+        j = s1; s1 = s2; s2 = j; // Switch rows
+      }
+    }
+
+    /* Multiply through by the correct factors */
+    /* (Eq. [2.9])                             */
+    r = p;
+    for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
+    {
+      for (DenseIndex j=p; j>=0; --j) N_(k,j) *= r;
+      r *= p-k;
+    }
+  }
+
+  template <typename _Scalar, int _Dim, int _Degree>
+  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
+    Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
+  {
+    typename SplineTraits< Spline >::BasisDerivativeType der;
+    basisFunctionDerivativesImpl(*this, u, order, der);
+    return der;
+  }
+
+  template <typename _Scalar, int _Dim, int _Degree>
+  template <int DerivativeOrder>
+  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
+    Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
+  {
+    typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der;
+    basisFunctionDerivativesImpl(*this, u, order, der);
+    return der;
+  }
+}
+
+#endif // EIGEN_SPLINE_H
diff --git a/unsupported/Eigen/src/Splines/SplineFitting.h b/unsupported/Eigen/src/Splines/SplineFitting.h
new file mode 100644
index 000000000..6de820af4
--- /dev/null
+++ b/unsupported/Eigen/src/Splines/SplineFitting.h
@@ -0,0 +1,115 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SPLINE_FITTING_H
+#define EIGEN_SPLINE_FITTING_H
+
+#include <numeric>
+
+#include "SplineFwd.h"
+
+#include <Eigen/QR>
+
+namespace Eigen
+{
+  template <typename KnotVectorType>
+  void KnotAveraging(const KnotVectorType& parameters, DenseIndex degree, KnotVectorType& knots)
+  {
+    typedef typename KnotVectorType::Scalar Scalar;
+
+    knots.resize(parameters.size()+degree+1);      
+
+    for (DenseIndex j=1; j<parameters.size()-degree; ++j)
+      knots(j+degree) = parameters.segment(j,degree).mean();
+
+    knots.segment(0,degree+1) = KnotVectorType::Zero(degree+1);
+    knots.segment(knots.size()-degree-1,degree+1) = KnotVectorType::Ones(degree+1);
+  }
+
+  template <typename PointArrayType, typename KnotVectorType>
+  void ChordLengths(const PointArrayType& pts, KnotVectorType& chord_lengths)
+  {
+    typedef typename KnotVectorType::Scalar Scalar;
+
+    const DenseIndex n = pts.cols();
+
+    // 1. compute the column-wise norms
+    chord_lengths.resize(pts.cols());
+    chord_lengths[0] = 0;
+    chord_lengths.rightCols(n-1) = (pts.array().leftCols(n-1) - pts.array().rightCols(n-1)).matrix().colwise().norm();
+
+    // 2. compute the partial sums
+    std::partial_sum(chord_lengths.data(), chord_lengths.data()+n, chord_lengths.data());
+
+    // 3. normalize the data
+    chord_lengths /= chord_lengths(n-1);
+    chord_lengths(n-1) = Scalar(1);
+  }
+
+  template <typename SplineType>
+  struct SplineFitting
+  {
+    template <typename PointArrayType>
+    static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree);
+  };
+
+  template <typename SplineType>
+  template <typename PointArrayType>
+  SplineType SplineFitting<SplineType>::Interpolate(const PointArrayType& pts, DenseIndex degree)
+  {
+    typedef typename SplineType::KnotVectorType::Scalar Scalar;      
+    typedef typename SplineType::KnotVectorType KnotVectorType;
+    typedef typename SplineType::BasisVectorType BasisVectorType;
+    typedef typename SplineType::ControlPointVectorType ControlPointVectorType;      
+
+    typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
+
+    KnotVectorType chord_lengths; // knot parameters
+    ChordLengths(pts, chord_lengths);
+
+    KnotVectorType knots;
+    KnotAveraging(chord_lengths, degree, knots);
+
+    DenseIndex n = pts.cols();
+    MatrixType A = MatrixType::Zero(n,n);
+    for (DenseIndex i=1; i<n-1; ++i)
+    {
+      const DenseIndex span = SplineType::Span(chord_lengths[i], degree, knots);
+
+      // The segment call should somehow be told the spline order at compile time.
+      A.row(i).segment(span-degree, degree+1) = SplineType::BasisFunctions(chord_lengths[i], degree, knots);
+    }
+    A(0,0) = 1.0;
+    A(n-1,n-1) = 1.0;
+
+    HouseholderQR<MatrixType> qr(A);
+
+    // Here, we are creating a temporary due to an Eigen issue.
+    ControlPointVectorType ctrls = qr.solve(MatrixType(pts.transpose())).transpose();
+
+    return SplineType(knots, ctrls);
+  }
+}
+
+#endif // EIGEN_SPLINE_FITTING_H
diff --git a/unsupported/Eigen/src/Splines/SplineFwd.h b/unsupported/Eigen/src/Splines/SplineFwd.h
new file mode 100644
index 000000000..da0e75205
--- /dev/null
+++ b/unsupported/Eigen/src/Splines/SplineFwd.h
@@ -0,0 +1,78 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SPLINES_FWD_H
+#define EIGEN_SPLINES_FWD_H
+
+#include <Eigen/Core>
+
+namespace Eigen
+{
+    template <typename Scalar, int Dim, int Degree = Dynamic> class Spline;
+
+    template < typename SplineType, int _DerivativeOrder = Dynamic > struct SplineTraits {};
+
+// hide specializations from doxygen
+#ifndef DOXYGEN_SHOULD_SKIP_THIS
+
+    template <typename _Scalar, int _Dim, int _Degree>
+    struct SplineTraits< Spline<_Scalar, _Dim, _Degree>, Dynamic >
+    {
+      typedef _Scalar Scalar; /* The underlying scalar value. */
+      enum { Dimension = _Dim }; /* The spline curve's dimension. */
+      enum { Degree = _Degree }; /* The spline curve's degree. */
+
+      enum { OrderAtCompileTime = _Degree==Dynamic ? Dynamic : _Degree+1 };
+      enum { NumOfDerivativesAtCompileTime = OrderAtCompileTime };
+
+      typedef Array<Scalar,1,OrderAtCompileTime> BasisVectorType;
+
+      typedef Array<Scalar,Dynamic,Dynamic,RowMajor,NumOfDerivativesAtCompileTime,OrderAtCompileTime> BasisDerivativeType;
+      typedef Array<Scalar,Dimension,Dynamic,ColMajor,Dimension,NumOfDerivativesAtCompileTime> DerivativeType;
+
+      typedef Array<Scalar,Dimension,1> PointType;
+      typedef Array<Scalar,1,Dynamic> KnotVectorType;
+      typedef Array<Scalar,Dimension,Dynamic> ControlPointVectorType;
+    };
+
+    template < typename _Scalar, int _Dim, int _Degree, int _DerivativeOrder >
+    struct SplineTraits< Spline<_Scalar, _Dim, _Degree>, _DerivativeOrder > : public SplineTraits< Spline<_Scalar, _Dim, _Degree> >
+    {
+      enum { OrderAtCompileTime = _Degree==Dynamic ? Dynamic : _Degree+1 };
+      enum { NumOfDerivativesAtCompileTime = _DerivativeOrder==Dynamic ? Dynamic : _DerivativeOrder+1 };
+
+      typedef Array<_Scalar,Dynamic,Dynamic,RowMajor,NumOfDerivativesAtCompileTime,OrderAtCompileTime> BasisDerivativeType;
+      typedef Array<_Scalar,_Dim,Dynamic,ColMajor,_Dim,NumOfDerivativesAtCompileTime> DerivativeType;
+    };
+
+#endif
+
+    typedef Spline<float,2> Spline2f;
+    typedef Spline<float,3> Spline3f;
+
+    typedef Spline<double,2> Spline2d;
+    typedef Spline<double,3> Spline3d;
+}
+
+#endif // EIGEN_SPLINES_FWD_H
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index 83e947665..8580bb061 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -94,3 +94,4 @@ endif(GSL_FOUND)
 ei_add_test(polynomialsolver " " "${GSL_LIBRARIES}" )
 ei_add_test(polynomialutils)
 ei_add_test(kronecker_product)
+ei_add_test(splines)
diff --git a/unsupported/test/splines.cpp b/unsupported/test/splines.cpp
new file mode 100644
index 000000000..02e170c44
--- /dev/null
+++ b/unsupported/test/splines.cpp
@@ -0,0 +1,239 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010-2011 Hauke Heibel <heibel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+#include <unsupported/Eigen/Splines>
+
+// lets do some explicit instantiations and thus
+// force the compilation of all spline functions...
+template class Spline<double, 2, Dynamic>;
+template class Spline<double, 3, Dynamic>;
+
+template class Spline<double, 2, 2>;
+template class Spline<double, 2, 3>;
+template class Spline<double, 2, 4>;
+template class Spline<double, 2, 5>;
+
+template class Spline<float, 2, Dynamic>;
+template class Spline<float, 3, Dynamic>;
+
+template class Spline<float, 3, 2>;
+template class Spline<float, 3, 3>;
+template class Spline<float, 3, 4>;
+template class Spline<float, 3, 5>;
+
+Spline<double, 2, Dynamic> closed_spline2d()
+{
+  RowVectorXd knots(12);
+  knots << 0,
+    0,
+    0,
+    0,
+    0.867193179093898,
+    1.660330955342408,
+    2.605084834823134,
+    3.484154586374428,
+    4.252699478956276,
+    4.252699478956276,
+    4.252699478956276,
+    4.252699478956276;
+
+  MatrixXd ctrls(8,2);
+  ctrls << -0.370967741935484,   0.236842105263158,
+    -0.231401860693277,   0.442245185027632,
+    0.344361228532831,   0.773369994120753,
+    0.828990216203802,   0.106550882647595,
+    0.407270163678382,  -1.043452922172848,
+    -0.488467813584053,  -0.390098582530090,
+    -0.494657189446427,   0.054804824897884,
+    -0.370967741935484,   0.236842105263158;
+  ctrls.transposeInPlace();
+
+  return Spline<double, 2, Dynamic>(knots, ctrls);
+}
+
+/* create a reference spline */
+Spline<double, 3, Dynamic> spline3d()
+{
+  RowVectorXd knots(11);
+  knots << 0,
+    0,
+    0,
+    0.118997681558377,
+    0.162611735194631,
+    0.498364051982143,
+    0.655098003973841,
+    0.679702676853675,
+    1.000000000000000,
+    1.000000000000000,
+    1.000000000000000;
+
+  MatrixXd ctrls(8,3);
+  ctrls <<    0.959743958516081,   0.340385726666133,   0.585267750979777,
+    0.223811939491137,   0.751267059305653,   0.255095115459269,
+    0.505957051665142,   0.699076722656686,   0.890903252535799,
+    0.959291425205444,   0.547215529963803,   0.138624442828679,
+    0.149294005559057,   0.257508254123736,   0.840717255983663,
+    0.254282178971531,   0.814284826068816,   0.243524968724989,
+    0.929263623187228,   0.349983765984809,   0.196595250431208,
+    0.251083857976031,   0.616044676146639,   0.473288848902729;
+  ctrls.transposeInPlace();
+
+  return Spline<double, 3, Dynamic>(knots, ctrls);
+}
+
+/* compares evaluations against known results */
+void eval_spline3d()
+{
+  Spline3d spline = spline3d();
+
+  RowVectorXd u(10);
+  u << 0.351659507062997,
+    0.830828627896291,
+    0.585264091152724,
+    0.549723608291140,
+    0.917193663829810,
+    0.285839018820374,
+    0.757200229110721,
+    0.753729094278495,
+    0.380445846975357,
+    0.567821640725221;
+
+  MatrixXd pts(10,3);
+  pts << 0.707620811535916,   0.510258911240815,   0.417485437023409,
+    0.603422256426978,   0.529498282727551,   0.270351549348981,
+    0.228364197569334,   0.423745615677815,   0.637687289287490,
+    0.275556796335168,   0.350856706427970,   0.684295784598905,
+    0.514519311047655,   0.525077224890754,   0.351628308305896,
+    0.724152914315666,   0.574461155457304,   0.469860285484058,
+    0.529365063753288,   0.613328702656816,   0.237837040141739,
+    0.522469395136878,   0.619099658652895,   0.237139665242069,
+    0.677357023849552,   0.480655768435853,   0.422227610314397,
+    0.247046593173758,   0.380604672404750,   0.670065791405019;
+  pts.transposeInPlace();
+
+  for (int i=0; i<u.size(); ++i)
+  {
+    Vector3d pt = spline(u(i));
+    VERIFY( (pt - pts.col(i)).norm() < 1e-14 );
+  }
+}
+
+/* compares evaluations on corner cases */
+void eval_spline3d_onbrks()
+{
+  Spline3d spline = spline3d();
+
+  RowVectorXd u = spline.knots();
+
+  MatrixXd pts(11,3);
+  pts <<    0.959743958516081,   0.340385726666133,   0.585267750979777,
+    0.959743958516081,   0.340385726666133,   0.585267750979777,
+    0.959743958516081,   0.340385726666133,   0.585267750979777,
+    0.430282980289940,   0.713074680056118,   0.720373307943349,
+    0.558074875553060,   0.681617921034459,   0.804417124839942,
+    0.407076008291750,   0.349707710518163,   0.617275937419545,
+    0.240037008286602,   0.738739390398014,   0.324554153129411,
+    0.302434111480572,   0.781162443963899,   0.240177089094644,
+    0.251083857976031,   0.616044676146639,   0.473288848902729,
+    0.251083857976031,   0.616044676146639,   0.473288848902729,
+    0.251083857976031,   0.616044676146639,   0.473288848902729;
+  pts.transposeInPlace();
+
+  for (int i=0; i<u.size(); ++i)
+  {
+    Vector3d pt = spline(u(i));
+    VERIFY( (pt - pts.col(i)).norm() < 1e-14 );
+  }
+}
+
+void eval_closed_spline2d()
+{
+  Spline2d spline = closed_spline2d();
+
+  RowVectorXd u(12);
+  u << 0,
+    0.332457030395796,
+    0.356467130532952,
+    0.453562180176215,
+    0.648017921874804,
+    0.973770235555003,
+    1.882577647219307,
+    2.289408593930498,
+    3.511951429883045,
+    3.884149321369450,
+    4.236261590369414,
+    4.252699478956276;
+
+  MatrixXd pts(12,2);
+  pts << -0.370967741935484,   0.236842105263158,
+    -0.152576775123250,   0.448975001279334,
+    -0.133417538277668,   0.461615613865667,
+    -0.053199060826740,   0.507630360006299,
+    0.114249591147281,   0.570414135097409,
+    0.377810316891987,   0.560497102875315,
+    0.665052120135908,  -0.157557441109611,
+    0.516006487053228,  -0.559763292174825,
+    -0.379486035348887,  -0.331959640488223,
+    -0.462034726249078,  -0.039105670080824,
+    -0.378730600917982,   0.225127015099919,
+    -0.370967741935484,   0.236842105263158;
+  pts.transposeInPlace();
+
+  for (int i=0; i<u.size(); ++i)
+  {
+    Vector2d pt = spline(u(i));
+    VERIFY( (pt - pts.col(i)).norm() < 1e-14 );
+  }
+}
+
+void check_global_interpolation2d()
+{
+  typedef Spline2d::PointType PointType;
+  typedef Spline2d::KnotVectorType KnotVectorType;
+  typedef Spline2d::ControlPointVectorType ControlPointVectorType;
+
+  ControlPointVectorType points = ControlPointVectorType::Random(2,100);
+  const Spline2d spline = SplineFitting<Spline2d>::Interpolate(points,3);
+
+  KnotVectorType chord_lengths; // knot parameters
+  Eigen::ChordLengths(points, chord_lengths);
+
+  for (Eigen::DenseIndex i=0; i<points.cols(); ++i)
+  {
+    PointType pt = spline( chord_lengths(i) );
+    PointType ref = points.col(i);
+    VERIFY( (pt - ref).matrix().norm() < 1e-14 );
+  }
+}
+
+
+void test_splines()
+{
+  CALL_SUBTEST( eval_spline3d() );
+  CALL_SUBTEST( eval_spline3d_onbrks() );
+  CALL_SUBTEST( eval_closed_spline2d() );
+  CALL_SUBTEST( check_global_interpolation2d() );
+}