Spline.h
Go to the documentation of this file.
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 
44 
47 
50 
53 
56 
59 
64  Spline()
65  : m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2))
66  , m_ctrls(ControlPointVectorType::Zero(Dimension,(Degree==Dynamic ? 1 : Degree+1)))
67  {
68  // in theory this code can go to the initializer list but it will get pretty
69  // much unreadable ...
70  enum { MinDegree = (Degree==Dynamic ? 0 : Degree) };
71  m_knots.template segment<MinDegree+1>(0) = Array<Scalar,1,MinDegree+1>::Zero();
72  m_knots.template segment<MinDegree+1>(MinDegree+1) = Array<Scalar,1,MinDegree+1>::Ones();
73  }
74 
80  template <typename OtherVectorType, typename OtherArrayType>
81  Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {}
82 
87  template <int OtherDegree>
89  m_knots(spline.knots()), m_ctrls(spline.ctrls()) {}
90 
94  const KnotVectorType& knots() const { return m_knots; }
95 
99  const ControlPointVectorType& ctrls() const { return m_ctrls; }
100 
112  PointType operator()(Scalar u) const;
113 
127  derivatives(Scalar u, DenseIndex order) const;
128 
134  template <int DerivativeOrder>
136  derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
137 
155  basisFunctions(Scalar u) const;
156 
171  basisFunctionDerivatives(Scalar u, DenseIndex order) const;
172 
178  template <int DerivativeOrder>
180  basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
181 
185  DenseIndex degree() const;
186 
191  DenseIndex span(Scalar u) const;
192 
197 
210  static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
211 
217  static BasisDerivativeType BasisFunctionDerivatives(
218  const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType& knots);
219 
220  private:
221  KnotVectorType m_knots;
222  ControlPointVectorType m_ctrls;
224  template <typename DerivativeType>
225  static void BasisFunctionDerivativesImpl(
227  const DenseIndex order,
228  const DenseIndex p,
230  DerivativeType& N_);
231  };
232 
233  template <typename _Scalar, int _Dim, int _Degree>
238  {
239  // Piegl & Tiller, "The NURBS Book", A2.1 (p. 68)
240  if (u <= knots(0)) return degree;
241  const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
242  return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 );
243  }
244 
245  template <typename _Scalar, int _Dim, int _Degree>
251  {
252  const DenseIndex p = degree;
253  const DenseIndex i = Spline::Span(u, degree, knots);
254 
255  const KnotVectorType& U = knots;
256 
257  BasisVectorType left(p+1); left(0) = Scalar(0);
258  BasisVectorType right(p+1); right(0) = Scalar(0);
259 
262 
263  BasisVectorType N(1,p+1);
264  N(0) = Scalar(1);
265  for (DenseIndex j=1; j<=p; ++j)
266  {
267  Scalar saved = Scalar(0);
268  for (DenseIndex r=0; r<j; r++)
269  {
270  const Scalar tmp = N(r)/(right(r+1)+left(j-r));
271  N[r] = saved + right(r+1)*tmp;
272  saved = left(j-r)*tmp;
273  }
274  N(j) = saved;
275  }
276  return N;
277  }
278 
279  template <typename _Scalar, int _Dim, int _Degree>
281  {
282  if (_Degree == Dynamic)
283  return m_knots.size() - m_ctrls.cols() - 1;
284  else
285  return _Degree;
286  }
287 
288  template <typename _Scalar, int _Dim, int _Degree>
290  {
291  return Spline::Span(u, degree(), knots());
292  }
293 
294  template <typename _Scalar, int _Dim, int _Degree>
296  {
298 
299  const DenseIndex span = this->span(u);
300  const DenseIndex p = degree();
301  const BasisVectorType basis_funcs = basisFunctions(u);
302 
303  const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
305  return (ctrl_weights * ctrl_pts).rowwise().sum();
306  }
307 
308  /* --------------------------------------------------------------------------------------------- */
309 
310  template <typename SplineType, typename DerivativeType>
311  void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
312  {
315  enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
316 
319  typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;
320 
321  const DenseIndex p = spline.degree();
322  const DenseIndex span = spline.span(u);
323 
324  const DenseIndex n = (std::min)(p, order);
325 
326  der.resize(Dimension,n+1);
327 
328  // Retrieve the basis function derivatives up to the desired order...
329  const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);
330 
331  // ... and perform the linear combinations of the control points.
332  for (DenseIndex der_order=0; der_order<n+1; ++der_order)
333  {
334  const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
335  const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
336  der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
337  }
338  }
339 
340  template <typename _Scalar, int _Dim, int _Degree>
341  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType
343  {
345  derivativesImpl(*this, u, order, res);
346  return res;
347  }
348 
349  template <typename _Scalar, int _Dim, int _Degree>
350  template <int DerivativeOrder>
351  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType
353  {
355  derivativesImpl(*this, u, order, res);
356  return res;
357  }
358 
359  template <typename _Scalar, int _Dim, int _Degree>
360  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType
362  {
363  return Spline::BasisFunctions(u, degree(), knots());
364  }
365 
366  /* --------------------------------------------------------------------------------------------- */
367 
368 
369  template <typename _Scalar, int _Dim, int _Degree>
370  template <typename DerivativeType>
373  const DenseIndex order,
374  const DenseIndex p,
376  DerivativeType& N_)
377  {
378  typedef Spline<_Scalar, _Dim, _Degree> SplineType;
380 
381  const DenseIndex span = SplineType::Span(u, p, U);
382 
383  const DenseIndex n = (std::min)(p, order);
384 
385  N_.resize(n+1, p+1);
386 
387  BasisVectorType left = BasisVectorType::Zero(p+1);
388  BasisVectorType right = BasisVectorType::Zero(p+1);
389 
390  Matrix<Scalar,Order,Order> ndu(p+1,p+1);
391 
392  Scalar saved, temp; // FIXME These were double instead of Scalar. Was there a reason for that?
393 
394  ndu(0,0) = 1.0;
395 
396  DenseIndex j;
397  for (j=1; j<=p; ++j)
398  {
399  left[j] = u-U[span+1-j];
400  right[j] = U[span+j]-u;
401  saved = 0.0;
402 
403  for (DenseIndex r=0; r<j; ++r)
404  {
405  /* Lower triangle */
406  ndu(j,r) = right[r+1]+left[j-r];
407  temp = ndu(r,j-1)/ndu(j,r);
408  /* Upper triangle */
409  ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp);
410  saved = left[j-r] * temp;
411  }
412 
413  ndu(j,j) = static_cast<Scalar>(saved);
414  }
415 
416  for (j = p; j>=0; --j)
417  N_(0,j) = ndu(j,p);
418 
419  // Compute the derivatives
420  DerivativeType a(n+1,p+1);
421  DenseIndex r=0;
422  for (; r<=p; ++r)
423  {
424  DenseIndex s1,s2;
425  s1 = 0; s2 = 1; // alternate rows in array a
426  a(0,0) = 1.0;
427 
428  // Compute the k-th derivative
429  for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
430  {
431  Scalar d = 0.0;
432  DenseIndex rk,pk,j1,j2;
433  rk = r-k; pk = p-k;
434 
435  if (r>=k)
436  {
437  a(s2,0) = a(s1,0)/ndu(pk+1,rk);
438  d = a(s2,0)*ndu(rk,pk);
439  }
440 
441  if (rk>=-1) j1 = 1;
442  else j1 = -rk;
443 
444  if (r-1 <= pk) j2 = k-1;
445  else j2 = p-r;
446 
447  for (j=j1; j<=j2; ++j)
448  {
449  a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
450  d += a(s2,j)*ndu(rk+j,pk);
451  }
452 
453  if (r<=pk)
454  {
455  a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
456  d += a(s2,k)*ndu(r,pk);
457  }
458 
459  N_(k,r) = static_cast<Scalar>(d);
460  j = s1; s1 = s2; s2 = j; // Switch rows
461  }
462  }
463 
464  /* Multiply through by the correct factors */
465  /* (Eq. [2.9]) */
466  r = p;
467  for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
468  {
469  for (j=p; j>=0; --j) N_(k,j) *= r;
470  r *= p-k;
471  }
472  }
473 
474  template <typename _Scalar, int _Dim, int _Degree>
475  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
477  {
479  BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
480  return der;
481  }
482 
483  template <typename _Scalar, int _Dim, int _Degree>
484  template <int DerivativeOrder>
485  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
487  {
488  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType der;
489  BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
490  return der;
491  }
492 
493  template <typename _Scalar, int _Dim, int _Degree>
497  const DenseIndex order,
498  const DenseIndex degree,
500  {
502  BasisFunctionDerivativesImpl(u, order, degree, knots, der);
503  return der;
504  }
505 }
506 
507 #endif // EIGEN_SPLINE_H
A class representing multi-dimensional spline curves.
Definition: Spline.h:35
void derivativesImpl(const SplineType &spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType &der)
Definition: Spline.h:311
SCALAR Scalar
Definition: bench_gemm.cpp:46
static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType &knots)
Returns the spline&#39;s non-zero basis functions.
Definition: Spline.h:247
SplineTraits< Spline >::PointType PointType
The point type the spline is representing.
Definition: Spline.h:43
KnotVectorType m_knots
Definition: Spline.h:221
SplineTraits< Spline >::DerivativeType derivatives(Scalar u, DenseIndex order) const
Evaluation of spline derivatives of up-to given order.
Definition: Spline.h:342
#define min(a, b)
Definition: datatypes.h:19
SplineTraits< Spline >::BasisDerivativeType BasisDerivativeType
The data type used to store the values of the basis function derivatives.
Definition: Spline.h:55
int n
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Double_ distance(const OrientedPlane3_ &p)
#define N
Definition: gksort.c:12
SplineTraits< Spline >::BasisDerivativeType basisFunctionDerivatives(Scalar u, DenseIndex order) const
Computes the non-zero spline basis function derivatives up to given order.
Definition: Spline.h:476
const ControlPointVectorType & ctrls() const
Returns the ctrls of the underlying spline.
Definition: Spline.h:99
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
static BasisDerivativeType BasisFunctionDerivatives(const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType &knots)
Computes the non-zero spline basis function derivatives up to given order.
Definition: Spline.h:495
static char left
Expression of a fixed-size or dynamic-size sub-vector.
ControlPointVectorType m_ctrls
Definition: Spline.h:222
PointType operator()(Scalar u) const
Returns the spline value at a given site .
Definition: Spline.h:295
DenseIndex span(Scalar u) const
Returns the span within the knot vector in which u is falling.
Definition: Spline.h:289
DenseIndex degree() const
Returns the spline degree.
Definition: Spline.h:280
Spline()
Creates a (constant) zero spline. For Splines with dynamic degree, the resulting degree will be 0...
Definition: Spline.h:64
Expression of the multiple replication of a matrix or vector.
Definition: Replicate.h:61
SplineTraits< Spline >::ParameterVectorType ParameterVectorType
The data type used to store parameter vectors.
Definition: Spline.h:49
const Block< const Derived, 1, internal::traits< Derived >::ColsAtCompileTime, IsRowMajor > ConstRowXpr
Definition: BlockMethods.h:18
static char right
static DenseIndex Span(typename SplineTraits< Spline >::Scalar u, DenseIndex degree, const typename SplineTraits< Spline >::KnotVectorType &knots)
Computes the span within the provided knot vector in which u is falling.
Definition: Spline.h:234
const KnotVectorType & knots() const
Returns the knots of the underlying spline.
Definition: Spline.h:94
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
Definition: Meta.h:66
SplineTraits< Spline >::ControlPointVectorType ControlPointVectorType
The data type representing the spline&#39;s control points.
Definition: Spline.h:58
Spline(const Spline< Scalar, Dimension, OtherDegree > &spline)
Copy constructor for splines.
Definition: Spline.h:88
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
SplineTraits< Spline >::BasisVectorType basisFunctions(Scalar u) const
Computes the non-zero basis functions at the given site.
Definition: Spline.h:361
Spline(const OtherVectorType &knots, const OtherArrayType &ctrls)
Creates a spline from a knot vector and control points.
Definition: Spline.h:81
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:45
float * p
SplineTraits< Spline >::BasisVectorType BasisVectorType
The data type used to store non-zero basis functions.
Definition: Spline.h:52
SplineTraits< Spline >::KnotVectorType KnotVectorType
The data type used to store knot vectors.
Definition: Spline.h:46
const int Dynamic
Definition: Constants.h:22
The matrix class, also used for vectors and row-vectors.
_Scalar Scalar
Definition: Spline.h:38
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_)
Definition: Spline.h:371
std::ptrdiff_t j


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:36:19