Spline.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SPLINE_H
11 #define EIGEN_SPLINE_H
12 
13 #include "SplineFwd.h"
14 
15 namespace Eigen
16 {
34  template <typename _Scalar, int _Dim, int _Degree>
35  class Spline
36  {
37  public:
38  typedef _Scalar Scalar;
39  enum { Dimension = _Dim };
40  enum { Degree = _Degree };
41 
43  typedef typename SplineTraits<Spline>::PointType PointType;
44 
46  typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;
47 
49  typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
50 
52  typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
53 
59  template <typename OtherVectorType, typename OtherArrayType>
60  Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {}
61 
66  template <int OtherDegree>
68  m_knots(spline.knots()), m_ctrls(spline.ctrls()) {}
69 
73  const KnotVectorType& knots() const { return m_knots; }
74 
78  const ControlPointVectorType& ctrls() const { return m_ctrls; }
79 
91  PointType operator()(Scalar u) const;
92 
105  typename SplineTraits<Spline>::DerivativeType
106  derivatives(Scalar u, DenseIndex order) const;
107 
113  template <int DerivativeOrder>
114  typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
115  derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
116 
133  typename SplineTraits<Spline>::BasisVectorType
134  basisFunctions(Scalar u) const;
135 
149  typename SplineTraits<Spline>::BasisDerivativeType
150  basisFunctionDerivatives(Scalar u, DenseIndex order) const;
151 
157  template <int DerivativeOrder>
158  typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
159  basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
160 
164  DenseIndex degree() const;
165 
170  DenseIndex span(Scalar u) const;
171 
175  static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots);
176 
189  static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
190 
191 
192  private:
193  KnotVectorType m_knots;
194  ControlPointVectorType m_ctrls;
195  };
196 
197  template <typename _Scalar, int _Dim, int _Degree>
199  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::Scalar u,
200  DenseIndex degree,
201  const typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::KnotVectorType& knots)
202  {
203  // Piegl & Tiller, "The NURBS Book", A2.1 (p. 68)
204  if (u <= knots(0)) return degree;
205  const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
206  return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 );
207  }
208 
209  template <typename _Scalar, int _Dim, int _Degree>
213  DenseIndex degree,
215  {
217 
218  const DenseIndex p = degree;
219  const DenseIndex i = Spline::Span(u, degree, knots);
220 
221  const KnotVectorType& U = knots;
222 
223  BasisVectorType left(p+1); left(0) = Scalar(0);
224  BasisVectorType right(p+1); right(0) = Scalar(0);
225 
228 
229  BasisVectorType N(1,p+1);
230  N(0) = Scalar(1);
231  for (DenseIndex j=1; j<=p; ++j)
232  {
233  Scalar saved = Scalar(0);
234  for (DenseIndex r=0; r<j; r++)
235  {
236  const Scalar tmp = N(r)/(right(r+1)+left(j-r));
237  N[r] = saved + right(r+1)*tmp;
238  saved = left(j-r)*tmp;
239  }
240  N(j) = saved;
241  }
242  return N;
243  }
244 
245  template <typename _Scalar, int _Dim, int _Degree>
247  {
248  if (_Degree == Dynamic)
249  return m_knots.size() - m_ctrls.cols() - 1;
250  else
251  return _Degree;
252  }
253 
254  template <typename _Scalar, int _Dim, int _Degree>
256  {
257  return Spline::Span(u, degree(), knots());
258  }
259 
260  template <typename _Scalar, int _Dim, int _Degree>
262  {
263  enum { Order = SplineTraits<Spline>::OrderAtCompileTime };
264 
265  const DenseIndex span = this->span(u);
266  const DenseIndex p = degree();
267  const BasisVectorType basis_funcs = basisFunctions(u);
268 
269  const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
270  const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
271  return (ctrl_weights * ctrl_pts).rowwise().sum();
272  }
273 
274  /* --------------------------------------------------------------------------------------------- */
275 
276  template <typename SplineType, typename DerivativeType>
277  void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
278  {
279  enum { Dimension = SplineTraits<SplineType>::Dimension };
280  enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
281  enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
282 
283  typedef typename SplineTraits<SplineType>::Scalar Scalar;
284 
285  typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
286  typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
287 
288  typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
289  typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;
290 
291  const DenseIndex p = spline.degree();
292  const DenseIndex span = spline.span(u);
293 
294  const DenseIndex n = (std::min)(p, order);
295 
296  der.resize(Dimension,n+1);
297 
298  // Retrieve the basis function derivatives up to the desired order...
299  const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);
300 
301  // ... and perform the linear combinations of the control points.
302  for (DenseIndex der_order=0; der_order<n+1; ++der_order)
303  {
304  const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
305  const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
306  der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
307  }
308  }
309 
310  template <typename _Scalar, int _Dim, int _Degree>
311  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType
312  Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
313  {
314  typename SplineTraits< Spline >::DerivativeType res;
315  derivativesImpl(*this, u, order, res);
316  return res;
317  }
318 
319  template <typename _Scalar, int _Dim, int _Degree>
320  template <int DerivativeOrder>
321  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType
322  Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
323  {
324  typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
325  derivativesImpl(*this, u, order, res);
326  return res;
327  }
328 
329  template <typename _Scalar, int _Dim, int _Degree>
330  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType
332  {
333  return Spline::BasisFunctions(u, degree(), knots());
334  }
335 
336  /* --------------------------------------------------------------------------------------------- */
337 
338  template <typename SplineType, typename DerivativeType>
339  void basisFunctionDerivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_)
340  {
341  enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
342 
343  typedef typename SplineTraits<SplineType>::Scalar Scalar;
344  typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
345  typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType;
346  typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
347 
348  const KnotVectorType& U = spline.knots();
349 
350  const DenseIndex p = spline.degree();
351  const DenseIndex span = spline.span(u);
352 
353  const DenseIndex n = (std::min)(p, order);
354 
355  N_.resize(n+1, p+1);
356 
357  BasisVectorType left = BasisVectorType::Zero(p+1);
358  BasisVectorType right = BasisVectorType::Zero(p+1);
359 
360  Matrix<Scalar,Order,Order> ndu(p+1,p+1);
361 
362  double saved, temp;
363 
364  ndu(0,0) = 1.0;
365 
366  DenseIndex j;
367  for (j=1; j<=p; ++j)
368  {
369  left[j] = u-U[span+1-j];
370  right[j] = U[span+j]-u;
371  saved = 0.0;
372 
373  for (DenseIndex r=0; r<j; ++r)
374  {
375  /* Lower triangle */
376  ndu(j,r) = right[r+1]+left[j-r];
377  temp = ndu(r,j-1)/ndu(j,r);
378  /* Upper triangle */
379  ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp);
380  saved = left[j-r] * temp;
381  }
382 
383  ndu(j,j) = static_cast<Scalar>(saved);
384  }
385 
386  for (j = p; j>=0; --j)
387  N_(0,j) = ndu(j,p);
388 
389  // Compute the derivatives
390  DerivativeType a(n+1,p+1);
391  DenseIndex r=0;
392  for (; r<=p; ++r)
393  {
394  DenseIndex s1,s2;
395  s1 = 0; s2 = 1; // alternate rows in array a
396  a(0,0) = 1.0;
397 
398  // Compute the k-th derivative
399  for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
400  {
401  double d = 0.0;
402  DenseIndex rk,pk,j1,j2;
403  rk = r-k; pk = p-k;
404 
405  if (r>=k)
406  {
407  a(s2,0) = a(s1,0)/ndu(pk+1,rk);
408  d = a(s2,0)*ndu(rk,pk);
409  }
410 
411  if (rk>=-1) j1 = 1;
412  else j1 = -rk;
413 
414  if (r-1 <= pk) j2 = k-1;
415  else j2 = p-r;
416 
417  for (j=j1; j<=j2; ++j)
418  {
419  a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
420  d += a(s2,j)*ndu(rk+j,pk);
421  }
422 
423  if (r<=pk)
424  {
425  a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
426  d += a(s2,k)*ndu(r,pk);
427  }
428 
429  N_(k,r) = static_cast<Scalar>(d);
430  j = s1; s1 = s2; s2 = j; // Switch rows
431  }
432  }
433 
434  /* Multiply through by the correct factors */
435  /* (Eq. [2.9]) */
436  r = p;
437  for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
438  {
439  for (DenseIndex j=p; j>=0; --j) N_(k,j) *= r;
440  r *= p-k;
441  }
442  }
443 
444  template <typename _Scalar, int _Dim, int _Degree>
445  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
447  {
448  typename SplineTraits< Spline >::BasisDerivativeType der;
449  basisFunctionDerivativesImpl(*this, u, order, der);
450  return der;
451  }
452 
453  template <typename _Scalar, int _Dim, int _Degree>
454  template <int DerivativeOrder>
455  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
456  Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
457  {
458  typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der;
459  basisFunctionDerivativesImpl(*this, u, order, der);
460  return der;
461  }
462 }
463 
464 #endif // EIGEN_SPLINE_H