GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-25 07:43:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Integrated spline class and simple spline fitting

Browse Source
This commit is contained in:
Hauke Heibel 2011-11-25 14:53:40 +01:00
parent 49d652c600
commit b00a33bc70
8 changed files with 895 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
unsupported/Eigen/CMakeLists.txt
View File

@ -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

52
unsupported/Eigen/Splines Normal file
View File

@ -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

1
unsupported/Eigen/src/CMakeLists.txt
View File

@ -10,3 +10,4 @@ ADD_SUBDIRECTORY(Polynomials)
ADD_SUBDIRECTORY(Skyline)
ADD_SUBDIRECTORY(SparseExtra)
ADD_SUBDIRECTORY(KroneckerProduct)
ADD_SUBDIRECTORY(Splines)

407
unsupported/Eigen/src/Splines/Spline.h Normal file
View File

@ -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

115
unsupported/Eigen/src/Splines/SplineFitting.h Normal file
View File

@ -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

78
unsupported/Eigen/src/Splines/SplineFwd.h Normal file
View File

@ -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

1
unsupported/test/CMakeLists.txt
View File

@ -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)

239
unsupported/test/splines.cpp Normal file
View File

@ -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() );
}