Eigen-unsupported  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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 "./InternalHeaderCheck.h"
14
15#include "SplineFwd.h"
16
17namespace Eigen
18{
36 template <typename Scalar_, int Dim_, int _Degree>
37 class Spline
38 {
39 public:
40 typedef Scalar_ Scalar;
41 enum { Dimension = Dim_ };
42 enum { Degree = _Degree };
43
45 typedef typename SplineTraits<Spline>::PointType PointType;
46
48 typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;
49
51 typedef typename SplineTraits<Spline>::ParameterVectorType ParameterVectorType;
52
54 typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
55
57 typedef typename SplineTraits<Spline>::BasisDerivativeType BasisDerivativeType;
58
60 typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
61
67 : m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2))
68 , m_ctrls(ControlPointVectorType::Zero(Dimension,(Degree==Dynamic ? 1 : Degree+1)))
69 {
70 // in theory this code can go to the initializer list but it will get pretty
71 // much unreadable ...
72 enum { MinDegree = (Degree==Dynamic ? 0 : Degree) };
73 m_knots.template segment<MinDegree+1>(0) = Array<Scalar,1,MinDegree+1>::Zero();
74 m_knots.template segment<MinDegree+1>(MinDegree+1) = Array<Scalar,1,MinDegree+1>::Ones();
75 }
76
82 template <typename OtherVectorType, typename OtherArrayType>
83 Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {}
84
89 template <int OtherDegree>
91 m_knots(spline.knots()), m_ctrls(spline.ctrls()) {}
92
96 const KnotVectorType& knots() const { return m_knots; }
97
101 const ControlPointVectorType& ctrls() const { return m_ctrls; }
102
114 PointType operator()(Scalar u) const;
115
128 typename SplineTraits<Spline>::DerivativeType
129 derivatives(Scalar u, DenseIndex order) const;
130
136 template <int DerivativeOrder>
137 typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
138 derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
139
156 typename SplineTraits<Spline>::BasisVectorType
157 basisFunctions(Scalar u) const;
158
172 typename SplineTraits<Spline>::BasisDerivativeType
173 basisFunctionDerivatives(Scalar u, DenseIndex order) const;
174
180 template <int DerivativeOrder>
181 typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
182 basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
183
187 DenseIndex degree() const;
188
193 DenseIndex span(Scalar u) const;
194
198 static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots);
199
212 static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
213
220 const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType& knots);
221
222 private:
223 KnotVectorType m_knots;
224 ControlPointVectorType m_ctrls;
226 template <typename DerivativeType>
227 static void BasisFunctionDerivativesImpl(
229 const DenseIndex order,
230 const DenseIndex p,
232 DerivativeType& N_);
233 };
234
235 template <typename Scalar_, int Dim_, int _Degree>
237 typename SplineTraits< Spline<Scalar_, Dim_, _Degree> >::Scalar u,
238 DenseIndex degree,
239 const typename SplineTraits< Spline<Scalar_, Dim_, _Degree> >::KnotVectorType& knots)
240 {
241 // Piegl & Tiller, "The NURBS Book", A2.1 (p. 68)
242 if (u <= knots(0)) return degree;
243 const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
244 return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 );
245 }
246
247 template <typename Scalar_, int Dim_, int _Degree>
251 DenseIndex degree,
253 {
254 const DenseIndex p = degree;
255 const DenseIndex i = Spline::Span(u, degree, knots);
256
257 const KnotVectorType& U = knots;
258
259 BasisVectorType left(p+1); left(0) = Scalar(0);
260 BasisVectorType right(p+1); right(0) = Scalar(0);
261
264
265 BasisVectorType N(1,p+1);
266 N(0) = Scalar(1);
267 for (DenseIndex j=1; j<=p; ++j)
268 {
269 Scalar saved = Scalar(0);
270 for (DenseIndex r=0; r<j; r++)
271 {
272 const Scalar tmp = N(r)/(right(r+1)+left(j-r));
273 N[r] = saved + right(r+1)*tmp;
274 saved = left(j-r)*tmp;
275 }
276 N(j) = saved;
277 }
278 return N;
279 }
280
281 template <typename Scalar_, int Dim_, int _Degree>
283 {
284 if (_Degree == Dynamic)
285 return m_knots.size() - m_ctrls.cols() - 1;
286 else
287 return _Degree;
288 }
289
290 template <typename Scalar_, int Dim_, int _Degree>
292 {
293 return Spline::Span(u, degree(), knots());
294 }
295
296 template <typename Scalar_, int Dim_, int _Degree>
298 {
299 enum { Order = SplineTraits<Spline>::OrderAtCompileTime };
300
301 const DenseIndex span = this->span(u);
302 const DenseIndex p = degree();
303 const BasisVectorType basis_funcs = basisFunctions(u);
304
305 const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
306 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
307 return (ctrl_weights * ctrl_pts).rowwise().sum();
308 }
309
310 /* --------------------------------------------------------------------------------------------- */
311
312 template <typename SplineType, typename DerivativeType>
313 void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
314 {
315 enum { Dimension = SplineTraits<SplineType>::Dimension };
316 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
317 enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
318
319 typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
320 typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
321 typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;
322
323 const DenseIndex p = spline.degree();
324 const DenseIndex span = spline.span(u);
325
326 const DenseIndex n = (std::min)(p, order);
327
328 der.resize(Dimension,n+1);
329
330 // Retrieve the basis function derivatives up to the desired order...
331 const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);
332
333 // ... and perform the linear combinations of the control points.
334 for (DenseIndex der_order=0; der_order<n+1; ++der_order)
335 {
336 const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
337 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
338 der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
339 }
340 }
341
342 template <typename Scalar_, int Dim_, int _Degree>
343 typename SplineTraits< Spline<Scalar_, Dim_, _Degree> >::DerivativeType
345 {
346 typename SplineTraits< Spline >::DerivativeType res;
347 derivativesImpl(*this, u, order, res);
348 return res;
349 }
350
351 template <typename Scalar_, int Dim_, int _Degree>
352 template <int DerivativeOrder>
353 typename SplineTraits< Spline<Scalar_, Dim_, _Degree>, DerivativeOrder >::DerivativeType
354 Spline<Scalar_, Dim_, _Degree>::derivatives(Scalar u, DenseIndex order) const
355 {
356 typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
357 derivativesImpl(*this, u, order, res);
358 return res;
359 }
360
361 template <typename Scalar_, int Dim_, int _Degree>
362 typename SplineTraits< Spline<Scalar_, Dim_, _Degree> >::BasisVectorType
364 {
365 return Spline::BasisFunctions(u, degree(), knots());
366 }
367
368 /* --------------------------------------------------------------------------------------------- */
369
370
371 template <typename Scalar_, int Dim_, int _Degree>
372 template <typename DerivativeType>
375 const DenseIndex order,
376 const DenseIndex p,
378 DerivativeType& N_)
379 {
380 typedef Spline<Scalar_, Dim_, _Degree> SplineType;
381 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
382
383 const DenseIndex span = SplineType::Span(u, p, U);
384
385 const DenseIndex n = (std::min)(p, order);
386
387 N_.resize(n+1, p+1);
388
389 BasisVectorType left = BasisVectorType::Zero(p+1);
390 BasisVectorType right = BasisVectorType::Zero(p+1);
391
392 Matrix<Scalar,Order,Order> ndu(p+1,p+1);
393
394 Scalar saved, temp; // FIXME These were double instead of Scalar. Was there a reason for that?
395
396 ndu(0,0) = 1.0;
397
398 DenseIndex j;
399 for (j=1; j<=p; ++j)
400 {
401 left[j] = u-U[span+1-j];
402 right[j] = U[span+j]-u;
403 saved = 0.0;
404
405 for (DenseIndex r=0; r<j; ++r)
406 {
407 /* Lower triangle */
408 ndu(j,r) = right[r+1]+left[j-r];
409 temp = ndu(r,j-1)/ndu(j,r);
410 /* Upper triangle */
411 ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp);
412 saved = left[j-r] * temp;
413 }
414
415 ndu(j,j) = static_cast<Scalar>(saved);
416 }
417
418 for (j = p; j>=0; --j)
419 N_(0,j) = ndu(j,p);
420
421 // Compute the derivatives
422 DerivativeType a(n+1,p+1);
423 DenseIndex r=0;
424 for (; r<=p; ++r)
425 {
426 DenseIndex s1,s2;
427 s1 = 0; s2 = 1; // alternate rows in array a
428 a(0,0) = 1.0;
429
430 // Compute the k-th derivative
431 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
432 {
433 Scalar d = 0.0;
434 DenseIndex rk,pk,j1,j2;
435 rk = r-k; pk = p-k;
436
437 if (r>=k)
438 {
439 a(s2,0) = a(s1,0)/ndu(pk+1,rk);
440 d = a(s2,0)*ndu(rk,pk);
441 }
442
443 if (rk>=-1) j1 = 1;
444 else j1 = -rk;
445
446 if (r-1 <= pk) j2 = k-1;
447 else j2 = p-r;
448
449 for (j=j1; j<=j2; ++j)
450 {
451 a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
452 d += a(s2,j)*ndu(rk+j,pk);
453 }
454
455 if (r<=pk)
456 {
457 a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
458 d += a(s2,k)*ndu(r,pk);
459 }
460
461 N_(k,r) = static_cast<Scalar>(d);
462 j = s1; s1 = s2; s2 = j; // Switch rows
463 }
464 }
465
466 /* Multiply through by the correct factors */
467 /* (Eq. [2.9]) */
468 r = p;
469 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
470 {
471 for (j=p; j>=0; --j) N_(k,j) *= r;
472 r *= p-k;
473 }
474 }
475
476 template <typename Scalar_, int Dim_, int _Degree>
477 typename SplineTraits< Spline<Scalar_, Dim_, _Degree> >::BasisDerivativeType
479 {
480 typename SplineTraits<Spline<Scalar_, Dim_, _Degree> >::BasisDerivativeType der;
481 BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
482 return der;
483 }
484
485 template <typename Scalar_, int Dim_, int _Degree>
486 template <int DerivativeOrder>
487 typename SplineTraits< Spline<Scalar_, Dim_, _Degree>, DerivativeOrder >::BasisDerivativeType
488 Spline<Scalar_, Dim_, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
489 {
490 typename SplineTraits< Spline<Scalar_, Dim_, _Degree>, DerivativeOrder >::BasisDerivativeType der;
491 BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
492 return der;
493 }
494
495 template <typename Scalar_, int Dim_, int _Degree>
496 typename SplineTraits<Spline<Scalar_, Dim_, _Degree> >::BasisDerivativeType
499 const DenseIndex order,
500 const DenseIndex degree,
502 {
503 typename SplineTraits<Spline>::BasisDerivativeType der;
504 BasisFunctionDerivativesImpl(u, order, degree, knots, der);
505 return der;
506 }
507}
508
509#endif // EIGEN_SPLINE_H
A class representing multi-dimensional spline curves.
Definition: Spline.h:38
SplineTraits< Spline >::BasisDerivativeType basisFunctionDerivatives(Scalar u, DenseIndex order) const
Computes the non-zero spline basis function derivatives up to given order.
Definition: Spline.h:478
SplineTraits< Spline >::BasisDerivativeType BasisDerivativeType
The data type used to store the values of the basis function derivatives.
Definition: Spline.h:57
const KnotVectorType & knots() const
Returns the knots of the underlying spline.
Definition: Spline.h:96
Scalar_ Scalar
Definition: Spline.h:40
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:236
Spline(const Spline< Scalar, Dimension, OtherDegree > &spline)
Copy constructor for splines.
Definition: Spline.h:90
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:497
SplineTraits< Spline >::ControlPointVectorType ControlPointVectorType
The data type representing the spline's control points.
Definition: Spline.h:60
SplineTraits< Spline >::DerivativeType derivatives(Scalar u, DenseIndex order) const
Evaluation of spline derivatives of up-to given order.
Definition: Spline.h:344
static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType &knots)
Returns the spline's non-zero basis functions.
Definition: Spline.h:249
SplineTraits< Spline >::KnotVectorType KnotVectorType
The data type used to store knot vectors.
Definition: Spline.h:48
Spline()
Creates a (constant) zero spline. For Splines with dynamic degree, the resulting degree will be 0.
Definition: Spline.h:66
PointType operator()(Scalar u) const
Returns the spline value at a given site .
Definition: Spline.h:297
SplineTraits< Spline >::BasisVectorType basisFunctions(Scalar u) const
Computes the non-zero basis functions at the given site.
Definition: Spline.h:363
DenseIndex degree() const
Returns the spline degree.
Definition: Spline.h:282
SplineTraits< Spline >::PointType PointType
The point type the spline is representing.
Definition: Spline.h:45
const ControlPointVectorType & ctrls() const
Returns the ctrls of the underlying spline.
Definition: Spline.h:101
@ Degree
Definition: Spline.h:42
SplineTraits< Spline, DerivativeOrder >::DerivativeType derivatives(Scalar u, DenseIndex order=DerivativeOrder) const
Evaluation of spline derivatives of up-to given order.
SplineTraits< Spline >::ParameterVectorType ParameterVectorType
The data type used to store parameter vectors.
Definition: Spline.h:51
SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType basisFunctionDerivatives(Scalar u, DenseIndex order=DerivativeOrder) const
Computes the non-zero spline basis function derivatives up to given order.
DenseIndex span(Scalar u) const
Returns the span within the knot vector in which u is falling.
Definition: Spline.h:291
SplineTraits< Spline >::BasisVectorType BasisVectorType
The data type used to store non-zero basis functions.
Definition: Spline.h:54
Spline(const OtherVectorType &knots, const OtherArrayType &ctrls)
Creates a spline from a knot vector and control points.
Definition: Spline.h:83
@ Dimension
Definition: Spline.h:41
Namespace containing all symbols from the Eigen library.
const int Dynamic