Spline.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_SPLINE_H
00011 #define EIGEN_SPLINE_H
00012 
00013 #include "SplineFwd.h"
00014 
00015 namespace Eigen
00016 {
00034   template <typename _Scalar, int _Dim, int _Degree>
00035   class Spline
00036   {
00037   public:
00038     typedef _Scalar Scalar; 
00039     enum { Dimension = _Dim  };
00040     enum { Degree = _Degree  };
00041 
00043     typedef typename SplineTraits<Spline>::PointType PointType;
00044     
00046     typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;
00047     
00049     typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
00050     
00052     typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
00053 
00059     template <typename OtherVectorType, typename OtherArrayType>
00060     Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {}
00061 
00066     template <int OtherDegree>
00067     Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) : 
00068     m_knots(spline.knots()), m_ctrls(spline.ctrls()) {}
00069 
00073     const KnotVectorType& knots() const { return m_knots; }
00074     
00078     const ControlPointVectorType& ctrls() const { return m_ctrls; }
00079 
00091     PointType operator()(Scalar u) const;
00092 
00105     typename SplineTraits<Spline>::DerivativeType
00106       derivatives(Scalar u, DenseIndex order) const;
00107 
00113     template <int DerivativeOrder>
00114     typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
00115       derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
00116 
00133     typename SplineTraits<Spline>::BasisVectorType
00134       basisFunctions(Scalar u) const;
00135 
00149     typename SplineTraits<Spline>::BasisDerivativeType
00150       basisFunctionDerivatives(Scalar u, DenseIndex order) const;
00151 
00157     template <int DerivativeOrder>
00158     typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
00159       basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
00160 
00164     DenseIndex degree() const;
00165 
00170     DenseIndex span(Scalar u) const;
00171 
00175     static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots);
00176     
00189     static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
00190 
00191 
00192   private:
00193     KnotVectorType m_knots; 
00194     ControlPointVectorType  m_ctrls; 
00195   };
00196 
00197   template <typename _Scalar, int _Dim, int _Degree>
00198   DenseIndex Spline<_Scalar, _Dim, _Degree>::Span(
00199     typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::Scalar u,
00200     DenseIndex degree,
00201     const typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::KnotVectorType& knots)
00202   {
00203     // Piegl & Tiller, "The NURBS Book", A2.1 (p. 68)
00204     if (u <= knots(0)) return degree;
00205     const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
00206     return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 );
00207   }
00208 
00209   template <typename _Scalar, int _Dim, int _Degree>
00210   typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType
00211     Spline<_Scalar, _Dim, _Degree>::BasisFunctions(
00212     typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
00213     DenseIndex degree,
00214     const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
00215   {
00216     typedef typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType BasisVectorType;
00217 
00218     const DenseIndex p = degree;
00219     const DenseIndex i = Spline::Span(u, degree, knots);
00220 
00221     const KnotVectorType& U = knots;
00222 
00223     BasisVectorType left(p+1); left(0) = Scalar(0);
00224     BasisVectorType right(p+1); right(0) = Scalar(0);        
00225 
00226     VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse();
00227     VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u;
00228 
00229     BasisVectorType N(1,p+1);
00230     N(0) = Scalar(1);
00231     for (DenseIndex j=1; j<=p; ++j)
00232     {
00233       Scalar saved = Scalar(0);
00234       for (DenseIndex r=0; r<j; r++)
00235       {
00236         const Scalar tmp = N(r)/(right(r+1)+left(j-r));
00237         N[r] = saved + right(r+1)*tmp;
00238         saved = left(j-r)*tmp;
00239       }
00240       N(j) = saved;
00241     }
00242     return N;
00243   }
00244 
00245   template <typename _Scalar, int _Dim, int _Degree>
00246   DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const
00247   {
00248     if (_Degree == Dynamic)
00249       return m_knots.size() - m_ctrls.cols() - 1;
00250     else
00251       return _Degree;
00252   }
00253 
00254   template <typename _Scalar, int _Dim, int _Degree>
00255   DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const
00256   {
00257     return Spline::Span(u, degree(), knots());
00258   }
00259 
00260   template <typename _Scalar, int _Dim, int _Degree>
00261   typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const
00262   {
00263     enum { Order = SplineTraits<Spline>::OrderAtCompileTime };
00264 
00265     const DenseIndex span = this->span(u);
00266     const DenseIndex p = degree();
00267     const BasisVectorType basis_funcs = basisFunctions(u);
00268 
00269     const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
00270     const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
00271     return (ctrl_weights * ctrl_pts).rowwise().sum();
00272   }
00273 
00274   /* --------------------------------------------------------------------------------------------- */
00275 
00276   template <typename SplineType, typename DerivativeType>
00277   void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
00278   {    
00279     enum { Dimension = SplineTraits<SplineType>::Dimension };
00280     enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
00281     enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
00282 
00283     typedef typename SplineTraits<SplineType>::Scalar Scalar;
00284 
00285     typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
00286     typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
00287 
00288     typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
00289     typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;    
00290 
00291     const DenseIndex p = spline.degree();
00292     const DenseIndex span = spline.span(u);
00293 
00294     const DenseIndex n = (std::min)(p, order);
00295 
00296     der.resize(Dimension,n+1);
00297 
00298     // Retrieve the basis function derivatives up to the desired order...    
00299     const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);
00300 
00301     // ... and perform the linear combinations of the control points.
00302     for (DenseIndex der_order=0; der_order<n+1; ++der_order)
00303     {
00304       const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
00305       const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
00306       der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
00307     }
00308   }
00309 
00310   template <typename _Scalar, int _Dim, int _Degree>
00311   typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType
00312     Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
00313   {
00314     typename SplineTraits< Spline >::DerivativeType res;
00315     derivativesImpl(*this, u, order, res);
00316     return res;
00317   }
00318 
00319   template <typename _Scalar, int _Dim, int _Degree>
00320   template <int DerivativeOrder>
00321   typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType
00322     Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
00323   {
00324     typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
00325     derivativesImpl(*this, u, order, res);
00326     return res;
00327   }
00328 
00329   template <typename _Scalar, int _Dim, int _Degree>
00330   typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType
00331     Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const
00332   {
00333     return Spline::BasisFunctions(u, degree(), knots());
00334   }
00335 
00336   /* --------------------------------------------------------------------------------------------- */
00337 
00338   template <typename SplineType, typename DerivativeType>
00339   void basisFunctionDerivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_)
00340   {
00341     enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
00342 
00343     typedef typename SplineTraits<SplineType>::Scalar Scalar;
00344     typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
00345     typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType;
00346     typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
00347 
00348     const KnotVectorType& U = spline.knots();
00349 
00350     const DenseIndex p = spline.degree();
00351     const DenseIndex span = spline.span(u);
00352 
00353     const DenseIndex n = (std::min)(p, order);
00354 
00355     N_.resize(n+1, p+1);
00356 
00357     BasisVectorType left = BasisVectorType::Zero(p+1);
00358     BasisVectorType right = BasisVectorType::Zero(p+1);
00359 
00360     Matrix<Scalar,Order,Order> ndu(p+1,p+1);
00361 
00362     double saved, temp;
00363 
00364     ndu(0,0) = 1.0;
00365 
00366     DenseIndex j;
00367     for (j=1; j<=p; ++j)
00368     {
00369       left[j] = u-U[span+1-j];
00370       right[j] = U[span+j]-u;
00371       saved = 0.0;
00372 
00373       for (DenseIndex r=0; r<j; ++r)
00374       {
00375         /* Lower triangle */
00376         ndu(j,r) = right[r+1]+left[j-r];
00377         temp = ndu(r,j-1)/ndu(j,r);
00378         /* Upper triangle */
00379         ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp);
00380         saved = left[j-r] * temp;
00381       }
00382 
00383       ndu(j,j) = static_cast<Scalar>(saved);
00384     }
00385 
00386     for (j = p; j>=0; --j) 
00387       N_(0,j) = ndu(j,p);
00388 
00389     // Compute the derivatives
00390     DerivativeType a(n+1,p+1);
00391     DenseIndex r=0;
00392     for (; r<=p; ++r)
00393     {
00394       DenseIndex s1,s2;
00395       s1 = 0; s2 = 1; // alternate rows in array a
00396       a(0,0) = 1.0;
00397 
00398       // Compute the k-th derivative
00399       for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
00400       {
00401         double d = 0.0;
00402         DenseIndex rk,pk,j1,j2;
00403         rk = r-k; pk = p-k;
00404 
00405         if (r>=k)
00406         {
00407           a(s2,0) = a(s1,0)/ndu(pk+1,rk);
00408           d = a(s2,0)*ndu(rk,pk);
00409         }
00410 
00411         if (rk>=-1) j1 = 1;
00412         else        j1 = -rk;
00413 
00414         if (r-1 <= pk) j2 = k-1;
00415         else           j2 = p-r;
00416 
00417         for (j=j1; j<=j2; ++j)
00418         {
00419           a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
00420           d += a(s2,j)*ndu(rk+j,pk);
00421         }
00422 
00423         if (r<=pk)
00424         {
00425           a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
00426           d += a(s2,k)*ndu(r,pk);
00427         }
00428 
00429         N_(k,r) = static_cast<Scalar>(d);
00430         j = s1; s1 = s2; s2 = j; // Switch rows
00431       }
00432     }
00433 
00434     /* Multiply through by the correct factors */
00435     /* (Eq. [2.9])                             */
00436     r = p;
00437     for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
00438     {
00439       for (DenseIndex j=p; j>=0; --j) N_(k,j) *= r;
00440       r *= p-k;
00441     }
00442   }
00443 
00444   template <typename _Scalar, int _Dim, int _Degree>
00445   typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
00446     Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
00447   {
00448     typename SplineTraits< Spline >::BasisDerivativeType der;
00449     basisFunctionDerivativesImpl(*this, u, order, der);
00450     return der;
00451   }
00452 
00453   template <typename _Scalar, int _Dim, int _Degree>
00454   template <int DerivativeOrder>
00455   typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
00456     Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
00457   {
00458     typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der;
00459     basisFunctionDerivativesImpl(*this, u, order, der);
00460     return der;
00461   }
00462 }
00463 
00464 #endif // EIGEN_SPLINE_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:26:27