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>::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>
67 Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) :
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>
198 DenseIndex Spline<_Scalar, _Dim, _Degree>::Span(
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>
210 typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType
211 Spline<_Scalar, _Dim, _Degree>::BasisFunctions(
212 typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
213 DenseIndex degree,
214 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
215 {
216 typedef typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType BasisVectorType;
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
226 VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse();
227 VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u;
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>
246 DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const
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>
255 DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const
256 {
257 return Spline::Span(u, degree(), knots());
258 }
259
260 template <typename _Scalar, int _Dim, int _Degree>
261 typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const
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
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
446 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
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
Eigen::Spline::derivatives
SplineTraits< Spline, DerivativeOrder >::DerivativeType derivatives(Scalar u, DenseIndex order=DerivativeOrder) const
Evaluation of spline derivatives of up-to given order.
Definition Spline.h:322
Eigen::Spline::Spline
Spline(const Spline< Scalar, Dimension, OtherDegree > &spline)
Copy constructor for splines.
Definition Spline.h:67
Eigen::Spline::PointType
SplineTraits< Spline >::PointType PointType
The point type the spline is representing.
Definition Spline.h:43
Eigen::Spline::BasisFunctions
static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType &knots)
Returns the spline's non-zero basis functions.
Definition Spline.h:211
Eigen::Spline< float, 2 >::knots
const KnotVectorType & knots() const
Definition Spline.h:73
Eigen::Spline::BasisVectorType
SplineTraits< Spline >::BasisVectorType BasisVectorType
The data type used to store non-zero basis functions.
Definition Spline.h:49
Eigen::Spline< float, 2 >::Dimension
@ Dimension
Definition Spline.h:39
Eigen::Spline::derivatives
SplineTraits< Spline >::DerivativeType derivatives(Scalar u, DenseIndex order) const
Evaluation of spline derivatives of up-to given order.
Definition Spline.h:312
Eigen::Spline::basisFunctionDerivatives
SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType basisFunctionDerivatives(Scalar u, DenseIndex order=DerivativeOrder) const
Computes the non-zero spline basis function derivatives up to given order.
Definition Spline.h:456
Eigen::Spline::KnotVectorType
SplineTraits< Spline >::KnotVectorType KnotVectorType
The data type used to store knot vectors.
Definition Spline.h:46
Eigen::Spline::basisFunctionDerivatives
SplineTraits< Spline >::BasisDerivativeType basisFunctionDerivatives(Scalar u, DenseIndex order) const
Computes the non-zero spline basis function derivatives up to given order.
Definition Spline.h:446
Eigen::Spline::span
DenseIndex span(Scalar u) const
Returns the span within the knot vector in which u is falling.
Definition Spline.h:255
Eigen::Spline::degree
DenseIndex degree() const
Returns the spline degree.
Definition Spline.h:246
Eigen::Spline::ControlPointVectorType
SplineTraits< Spline >::ControlPointVectorType ControlPointVectorType
The data type representing the spline's control points.
Definition Spline.h:52
Eigen::Spline< float, 2 >::ctrls
const ControlPointVectorType & ctrls() const
Definition Spline.h:78
Eigen::Spline< float, 2 >::Degree
@ Degree
Definition Spline.h:40
Eigen::Spline::Span
static DenseIndex Span(typename SplineTraits< Spline >::Scalar u, DenseIndex degree, const typename SplineTraits< Spline >::KnotVectorType &knots)
Computes the spang within the provided knot vector in which u is falling.
Definition Spline.h:198
Eigen::Spline::operator()
PointType operator()(Scalar u) const
Returns the spline value at a given site .
Definition Spline.h:261
Eigen::Spline::Spline
Spline(const OtherVectorType &knots, const OtherArrayType &ctrls)
Creates a spline from a knot vector and control points.
Definition Spline.h:60
Eigen::Spline::basisFunctions
SplineTraits< Spline >::BasisVectorType basisFunctions(Scalar u) const
Computes the non-zero basis functions at the given site.
Definition Spline.h:331
Eigen::Spline::Scalar
_Scalar Scalar
Definition Spline.h:38
Generated by doxygen 1.13.2