25 #ifndef EIGEN_SPLINE_H
26 #define EIGEN_SPLINE_H
28 #include "SplineFwd.h"
49 template <
typename _Scalar,
int _Dim,
int _Degree>
58 typedef typename SplineTraits<Spline>::PointType
PointType;
74 template <
typename OtherVectorType,
typename OtherArrayType>
75 Spline(
const OtherVectorType&
knots,
const OtherArrayType&
ctrls) : m_knots(knots), m_ctrls(ctrls) {}
81 template <
int OtherDegree>
83 m_knots(spline.
knots()), m_ctrls(spline.
ctrls()) {}
120 typename SplineTraits<Spline>::DerivativeType
128 template <
int DerivativeOrder>
129 typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
148 typename SplineTraits<Spline>::BasisVectorType
164 typename SplineTraits<Spline>::BasisDerivativeType
172 template <
int DerivativeOrder>
173 typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
179 DenseIndex
degree()
const;
190 static DenseIndex
Span(
typename SplineTraits<Spline>::Scalar u, DenseIndex
degree,
const typename SplineTraits<Spline>::KnotVectorType&
knots);
212 template <
typename _Scalar,
int _Dim,
int _Degree>
219 if (u <= knots(0))
return degree;
220 const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
221 return static_cast<DenseIndex
>( std::distance(knots.data(), pos) - 1 );
224 template <
typename _Scalar,
int _Dim,
int _Degree>
233 const DenseIndex p = degree;
238 BasisVectorType left(p+1); left(0) =
Scalar(0);
239 BasisVectorType right(p+1); right(0) =
Scalar(0);
241 VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse();
242 VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u;
244 BasisVectorType N(1,p+1);
246 for (DenseIndex j=1; j<=p; ++j)
249 for (DenseIndex r=0; r<j; r++)
251 const Scalar tmp = N(r)/(right(r+1)+left(j-r));
252 N[r] = saved + right(r+1)*tmp;
253 saved = left(j-r)*tmp;
260 template <
typename _Scalar,
int _Dim,
int _Degree>
264 return m_knots.size() - m_ctrls.cols() - 1;
269 template <
typename _Scalar,
int _Dim,
int _Degree>
275 template <
typename _Scalar,
int _Dim,
int _Degree>
278 enum { Order = SplineTraits<Spline>::OrderAtCompileTime };
280 const DenseIndex span = this->span(u);
281 const DenseIndex p = degree();
284 const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
285 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
286 return (ctrl_weights * ctrl_pts).rowwise().sum();
291 template <
typename SplineType,
typename DerivativeType>
292 void derivativesImpl(
const SplineType& spline,
typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
294 enum { Dimension = SplineTraits<SplineType>::Dimension };
295 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
296 enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
298 typedef typename SplineTraits<SplineType>::Scalar Scalar;
300 typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
301 typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
303 typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
304 typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;
306 const DenseIndex p = spline.degree();
307 const DenseIndex span = spline.span(u);
309 const DenseIndex n = (std::min)(p, order);
311 der.resize(Dimension,n+1);
314 const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);
317 for (DenseIndex der_order=0; der_order<n+1; ++der_order)
319 const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
320 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
321 der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
325 template <
typename _Scalar,
int _Dim,
int _Degree>
326 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType
329 typename SplineTraits< Spline >::DerivativeType res;
330 derivativesImpl(*
this, u, order, res);
334 template <
typename _Scalar,
int _Dim,
int _Degree>
335 template <
int DerivativeOrder>
336 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType
339 typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
340 derivativesImpl(*
this, u, order, res);
344 template <
typename _Scalar,
int _Dim,
int _Degree>
345 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType
353 template <
typename SplineType,
typename DerivativeType>
354 void basisFunctionDerivativesImpl(
const SplineType& spline,
typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_)
356 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
358 typedef typename SplineTraits<SplineType>::Scalar Scalar;
359 typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
360 typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType;
361 typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
363 const KnotVectorType& U = spline.knots();
365 const DenseIndex p = spline.degree();
366 const DenseIndex span = spline.span(u);
368 const DenseIndex n = (std::min)(p, order);
372 BasisVectorType left = BasisVectorType::Zero(p+1);
373 BasisVectorType right = BasisVectorType::Zero(p+1);
375 Matrix<Scalar,Order,Order> ndu(p+1,p+1);
384 left[j] = u-U[span+1-j];
385 right[j] = U[span+j]-u;
388 for (DenseIndex r=0; r<j; ++r)
391 ndu(j,r) = right[r+1]+left[j-r];
392 temp = ndu(r,j-1)/ndu(j,r);
394 ndu(r,j) =
static_cast<Scalar
>(saved+right[r+1] * temp);
395 saved = left[j-r] * temp;
398 ndu(j,j) =
static_cast<Scalar
>(saved);
401 for (j = p; j>=0; --j)
405 DerivativeType a(n+1,p+1);
414 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
417 DenseIndex rk,pk,j1,j2;
422 a(s2,0) = a(s1,0)/ndu(pk+1,rk);
423 d = a(s2,0)*ndu(rk,pk);
429 if (r-1 <= pk) j2 = k-1;
432 for (j=j1; j<=j2; ++j)
434 a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
435 d += a(s2,j)*ndu(rk+j,pk);
440 a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
441 d += a(s2,k)*ndu(r,pk);
444 N_(k,r) =
static_cast<Scalar
>(d);
445 j = s1; s1 = s2; s2 = j;
452 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
454 for (DenseIndex j=p; j>=0; --j) N_(k,j) *= r;
459 template <
typename _Scalar,
int _Dim,
int _Degree>
460 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
463 typename SplineTraits< Spline >::BasisDerivativeType der;
464 basisFunctionDerivativesImpl(*
this, u, order, der);
468 template <
typename _Scalar,
int _Dim,
int _Degree>
469 template <
int DerivativeOrder>
470 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
473 typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der;
474 basisFunctionDerivativesImpl(*
this, u, order, der);
479 #endif // EIGEN_SPLINE_H