| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #ifndef EIGEN_SPLINE_H |
| #define EIGEN_SPLINE_H |
| |
| #include "./InternalHeaderCheck.h" |
| |
| #include "SplineFwd.h" |
| |
| namespace Eigen |
| { |
| /** |
| * \ingroup Splines_Module |
| * \class Spline |
| * \brief A class representing multi-dimensional spline curves. |
| * |
| * The class represents B-splines with non-uniform knot vectors. Each control |
| * point of the B-spline is associated with a basis function |
| * \f{align*} |
| * C(u) & = \sum_{i=0}^{n}N_{i,p}(u)P_i |
| * \f} |
| * |
| * \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; /*!< The spline curve's scalar type. */ |
| enum { Dimension = Dim_ /*!< The spline curve's dimension. */ }; |
| enum { Degree = Degree_ /*!< The spline curve's degree. */ }; |
| |
| /** \brief The point type the spline is representing. */ |
| typedef typename SplineTraits<Spline>::PointType PointType; |
| |
| /** \brief The data type used to store knot vectors. */ |
| typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType; |
| |
| /** \brief The data type used to store parameter vectors. */ |
| typedef typename SplineTraits<Spline>::ParameterVectorType ParameterVectorType; |
| |
| /** \brief The data type used to store non-zero basis functions. */ |
| typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType; |
| |
| /** \brief The data type used to store the values of the basis function derivatives. */ |
| typedef typename SplineTraits<Spline>::BasisDerivativeType BasisDerivativeType; |
| |
| /** \brief The data type representing the spline's control points. */ |
| typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType; |
| |
| /** |
| * \brief Creates a (constant) zero spline. |
| * For Splines with dynamic degree, the resulting degree will be 0. |
| **/ |
| Spline() |
| : m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2)) |
| , m_ctrls(ControlPointVectorType::Zero(Dimension,(Degree==Dynamic ? 1 : Degree+1))) |
| { |
| // in theory this code can go to the initializer list but it will get pretty |
| // much unreadable ... |
| enum { MinDegree = (Degree==Dynamic ? 0 : Degree) }; |
| m_knots.template segment<MinDegree+1>(0) = Array<Scalar,1,MinDegree+1>::Zero(); |
| m_knots.template segment<MinDegree+1>(MinDegree+1) = Array<Scalar,1,MinDegree+1>::Ones(); |
| } |
| |
| /** |
| * \brief Creates a spline from a knot vector and control points. |
| * \param knots The spline's knot vector. |
| * \param ctrls The spline's control point vector. |
| **/ |
| template <typename OtherVectorType, typename OtherArrayType> |
| Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {} |
| |
| /** |
| * \brief Copy constructor for splines. |
| * \param spline The input spline. |
| **/ |
| template <int OtherDegree> |
| Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) : |
| m_knots(spline.knots()), m_ctrls(spline.ctrls()) {} |
| |
| /** |
| * \brief Returns the knots of the underlying spline. |
| **/ |
| const KnotVectorType& knots() const { return m_knots; } |
| |
| /** |
| * \brief Returns the ctrls of the underlying spline. |
| **/ |
| const ControlPointVectorType& ctrls() const { return m_ctrls; } |
| |
| /** |
| * \brief Returns the spline value at a given site \f$u\f$. |
| * |
| * The function returns |
| * \f{align*} |
| * C(u) & = \sum_{i=0}^{n}N_{i,p}P_i |
| * \f} |
| * |
| * \param u Parameter \f$u \in [0;1]\f$ at which the spline is evaluated. |
| * \return The spline value at the given location \f$u\f$. |
| **/ |
| PointType operator()(Scalar u) const; |
| |
| /** |
| * \brief Evaluation of spline derivatives of up-to given order. |
| * |
| * The function returns |
| * \f{align*} |
| * \frac{d^i}{du^i}C(u) & = \sum_{i=0}^{n} \frac{d^i}{du^i} N_{i,p}(u)P_i |
| * \f} |
| * for i ranging between 0 and order. |
| * |
| * \param u Parameter \f$u \in [0;1]\f$ at which the spline derivative is evaluated. |
| * \param order The order up to which the derivatives are computed. |
| **/ |
| typename SplineTraits<Spline>::DerivativeType |
| derivatives(Scalar u, DenseIndex order) const; |
| |
| /** |
| * \copydoc Spline::derivatives |
| * Using the template version of this function is more efficieent since |
| * temporary objects are allocated on the stack whenever this is possible. |
| **/ |
| template <int DerivativeOrder> |
| typename SplineTraits<Spline,DerivativeOrder>::DerivativeType |
| derivatives(Scalar u, DenseIndex order = DerivativeOrder) const; |
| |
| /** |
| * \brief Computes the non-zero basis functions at the given site. |
| * |
| * Splines have local support and a point from their image is defined |
| * by exactly \f$p+1\f$ control points \f$P_i\f$ where \f$p\f$ is the |
| * spline degree. |
| * |
| * This function computes the \f$p+1\f$ non-zero basis function values |
| * for a given parameter value \f$u\f$. It returns |
| * \f{align*}{ |
| * N_{i,p}(u), \hdots, N_{i+p+1,p}(u) |
| * \f} |
| * |
| * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis functions |
| * are computed. |
| **/ |
| typename SplineTraits<Spline>::BasisVectorType |
| basisFunctions(Scalar u) const; |
| |
| /** |
| * \brief Computes the non-zero spline basis function derivatives up to given order. |
| * |
| * The function computes |
| * \f{align*}{ |
| * \frac{d^i}{du^i} N_{i,p}(u), \hdots, \frac{d^i}{du^i} N_{i+p+1,p}(u) |
| * \f} |
| * with i ranging from 0 up to the specified order. |
| * |
| * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis function |
| * derivatives are computed. |
| * \param order The order up to which the basis function derivatives are computes. |
| **/ |
| typename SplineTraits<Spline>::BasisDerivativeType |
| basisFunctionDerivatives(Scalar u, DenseIndex order) const; |
| |
| /** |
| * \copydoc Spline::basisFunctionDerivatives |
| * Using the template version of this function is more efficieent since |
| * temporary objects are allocated on the stack whenever this is possible. |
| **/ |
| template <int DerivativeOrder> |
| typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType |
| basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const; |
| |
| /** |
| * \brief Returns the spline degree. |
| **/ |
| DenseIndex degree() const; |
| |
| /** |
| * \brief Returns the span within the knot vector in which u is falling. |
| * \param u The site for which the span is determined. |
| **/ |
| DenseIndex span(Scalar u) const; |
| |
| /** |
| * \brief Computes the span within the provided knot vector in which u is falling. |
| **/ |
| static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots); |
| |
| /** |
| * \brief Returns the spline's non-zero basis functions. |
| * |
| * The function computes and returns |
| * \f{align*}{ |
| * N_{i,p}(u), \hdots, N_{i+p+1,p}(u) |
| * \f} |
| * |
| * \param u The site at which the basis functions are computed. |
| * \param degree The degree of the underlying spline. |
| * \param knots The underlying spline's knot vector. |
| **/ |
| static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots); |
| |
| /** |
| * \copydoc Spline::basisFunctionDerivatives |
| * \param degree The degree of the underlying spline |
| * \param knots The underlying spline's knot vector. |
| **/ |
| static BasisDerivativeType BasisFunctionDerivatives( |
| const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType& knots); |
| |
| private: |
| KnotVectorType m_knots; /*!< Knot vector. */ |
| ControlPointVectorType m_ctrls; /*!< Control points. */ |
| |
| template <typename DerivativeType> |
| static void BasisFunctionDerivativesImpl( |
| const typename Spline<Scalar_, Dim_, Degree_>::Scalar u, |
| const DenseIndex order, |
| const DenseIndex p, |
| const typename Spline<Scalar_, Dim_, Degree_>::KnotVectorType& U, |
| DerivativeType& N_); |
| }; |
| |
| 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) |
| { |
| 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()); |
| } |
| |
| 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>::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(); |
| } |
| } |
| |
| 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; |
| } |
| |
| 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 Scalar_, int Dim_, int Degree_> |
| template <typename DerivativeType> |
| void Spline<Scalar_, Dim_, Degree_>::BasisFunctionDerivativesImpl( |
| const typename Spline<Scalar_, Dim_, Degree_>::Scalar u, |
| const DenseIndex order, |
| const DenseIndex p, |
| const typename Spline<Scalar_, Dim_, Degree_>::KnotVectorType& U, |
| DerivativeType& N_) |
| { |
| typedef Spline<Scalar_, Dim_, Degree_> SplineType; |
| enum { Order = SplineTraits<SplineType>::OrderAtCompileTime }; |
| |
| const DenseIndex span = SplineType::Span(u, p, 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); |
| |
| Scalar saved, temp; // FIXME These were double instead of Scalar. Was there a reason for that? |
| |
| 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) |
| { |
| Scalar 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 (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<Scalar_, Dim_, Degree_> >::BasisDerivativeType der; |
| BasisFunctionDerivativesImpl(u, order, degree(), knots(), 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<Scalar_, Dim_, Degree_>, DerivativeOrder >::BasisDerivativeType der; |
| BasisFunctionDerivativesImpl(u, order, degree(), knots(), der); |
| return der; |
| } |
| |
| template <typename Scalar_, int Dim_, int Degree_> |
| typename SplineTraits<Spline<Scalar_, Dim_, Degree_> >::BasisDerivativeType |
| Spline<Scalar_, Dim_, Degree_>::BasisFunctionDerivatives( |
| const typename Spline<Scalar_, Dim_, Degree_>::Scalar u, |
| const DenseIndex order, |
| const DenseIndex degree, |
| const typename Spline<Scalar_, Dim_, Degree_>::KnotVectorType& knots) |
| { |
| typename SplineTraits<Spline>::BasisDerivativeType der; |
| BasisFunctionDerivativesImpl(u, order, degree, knots, der); |
| return der; |
| } |
| } |
| |
| #endif // EIGEN_SPLINE_H |
