Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
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
15namespace 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>::ParameterVectorType ParameterVectorType;
50
52 typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
53
55 typedef typename SplineTraits<Spline>::BasisDerivativeType BasisDerivativeType;
56
58 typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
59
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
126 typename SplineTraits<Spline>::DerivativeType
127 derivatives(Scalar u, DenseIndex order) const;
128
134 template <int DerivativeOrder>
135 typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
136 derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
137
154 typename SplineTraits<Spline>::BasisVectorType
155 basisFunctions(Scalar u) const;
156
170 typename SplineTraits<Spline>::BasisDerivativeType
171 basisFunctionDerivatives(Scalar u, DenseIndex order) const;
172
178 template <int DerivativeOrder>
179 typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
180 basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
181
185 DenseIndex degree() const;
186
191 DenseIndex span(Scalar u) const;
192
196 static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots);
197
210 static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
211
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>
235 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::Scalar u,
236 DenseIndex degree,
237 const typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::KnotVectorType& knots)
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>
249 DenseIndex 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 {
297 enum { Order = SplineTraits<Spline>::OrderAtCompileTime };
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);
304 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
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 {
313 enum { Dimension = SplineTraits<SplineType>::Dimension };
314 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
315 enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
316
317 typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
318 typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
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 {
344 typename SplineTraits< Spline >::DerivativeType res;
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
352 Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
353 {
354 typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
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;
379 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
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 {
478 typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType der;
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
486 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
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>
494 typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
497 const DenseIndex order,
498 const DenseIndex degree,
500 {
501 typename SplineTraits<Spline>::BasisDerivativeType der;
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:36
SplineTraits< Spline >::ParameterVectorType ParameterVectorType
The data type used to store parameter vectors.
Definition: Spline.h:49
SplineTraits< Spline >::KnotVectorType KnotVectorType
The data type used to store knot vectors.
Definition: Spline.h:46
DenseIndex degree() const
Returns the spline degree.
Definition: Spline.h:280
Spline(const Spline< Scalar, Dimension, OtherDegree > &spline)
Copy constructor for splines.
Definition: Spline.h:88
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
SplineTraits< Spline >::DerivativeType derivatives(Scalar u, DenseIndex order) const
Evaluation of spline derivatives of up-to given order.
Definition: Spline.h:342
SplineTraits< Spline >::BasisVectorType BasisVectorType
The data type used to store non-zero basis functions.
Definition: Spline.h:52
const ControlPointVectorType & ctrls() const
Returns the ctrls of the underlying spline.
Definition: Spline.h:99
Spline()
Creates a (constant) zero spline. For Splines with dynamic degree, the resulting degree will be 0.
Definition: Spline.h:64
PointType operator()(Scalar u) const
Returns the spline value at a given site .
Definition: Spline.h:295
_Scalar Scalar
Definition: Spline.h:38
@ Degree
Definition: Spline.h:40
SplineTraits< Spline >::PointType PointType
The point type the spline is representing.
Definition: Spline.h:43
SplineTraits< Spline >::BasisDerivativeType BasisDerivativeType
The data type used to store the values of the basis function derivatives.
Definition: Spline.h:55
SplineTraits< Spline, DerivativeOrder >::DerivativeType derivatives(Scalar u, DenseIndex order=DerivativeOrder) const
Evaluation of spline derivatives of up-to given order.
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
static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType &knots)
Returns the spline's non-zero basis functions.
Definition: Spline.h:247
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
DenseIndex span(Scalar u) const
Returns the span within the knot vector in which u is falling.
Definition: Spline.h:289
SplineTraits< Spline >::ControlPointVectorType ControlPointVectorType
The data type representing the spline's control points.
Definition: Spline.h:58
SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType basisFunctionDerivatives(Scalar u, DenseIndex order=DerivativeOrder) const
Computes the non-zero spline basis function derivatives up to given order.
Spline(const OtherVectorType &knots, const OtherArrayType &ctrls)
Creates a spline from a knot vector and control points.
Definition: Spline.h:81
@ Dimension
Definition: Spline.h:39
const KnotVectorType & knots() const
Returns the knots of the underlying spline.
Definition: Spline.h:94
SplineTraits< Spline >::BasisVectorType basisFunctions(Scalar u) const
Computes the non-zero basis functions at the given site.
Definition: Spline.h:361
Namespace containing all symbols from the Eigen library.
const int Dynamic