source: pacpussensors/trunk/Vislab/lib3dv/eigen/unsupported/Eigen/src/Splines/Spline.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 16.4 KB
Line 
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{
17 /**
18 * \ingroup Splines_Module
19 * \class Spline
20 * \brief A class representing multi-dimensional spline curves.
21 *
22 * The class represents B-splines with non-uniform knot vectors. Each control
23 * point of the B-spline is associated with a basis function
24 * \f{align*}
25 * C(u) & = \sum_{i=0}^{n}N_{i,p}(u)P_i
26 * \f}
27 *
28 * \tparam _Scalar The underlying data type (typically float or double)
29 * \tparam _Dim The curve dimension (e.g. 2 or 3)
30 * \tparam _Degree Per default set to Dynamic; could be set to the actual desired
31 * degree for optimization purposes (would result in stack allocation
32 * of several temporary variables).
33 **/
34 template <typename _Scalar, int _Dim, int _Degree>
35 class Spline
36 {
37 public:
38 typedef _Scalar Scalar; /*!< The spline curve's scalar type. */
39 enum { Dimension = _Dim /*!< The spline curve's dimension. */ };
40 enum { Degree = _Degree /*!< The spline curve's degree. */ };
41
42 /** \brief The point type the spline is representing. */
43 typedef typename SplineTraits<Spline>::PointType PointType;
44
45 /** \brief The data type used to store knot vectors. */
46 typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;
47
48 /** \brief The data type used to store non-zero basis functions. */
49 typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
50
51 /** \brief The data type representing the spline's control points. */
52 typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
53
54 /**
55 * \brief Creates a (constant) zero spline.
56 * For Splines with dynamic degree, the resulting degree will be 0.
57 **/
58 Spline()
59 : m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2))
60 , m_ctrls(ControlPointVectorType::Zero(2,(Degree==Dynamic ? 1 : Degree+1)))
61 {
62 // in theory this code can go to the initializer list but it will get pretty
63 // much unreadable ...
64 enum { MinDegree = (Degree==Dynamic ? 0 : Degree) };
65 m_knots.template segment<MinDegree+1>(0) = Array<Scalar,1,MinDegree+1>::Zero();
66 m_knots.template segment<MinDegree+1>(MinDegree+1) = Array<Scalar,1,MinDegree+1>::Ones();
67 }
68
69 /**
70 * \brief Creates a spline from a knot vector and control points.
71 * \param knots The spline's knot vector.
72 * \param ctrls The spline's control point vector.
73 **/
74 template <typename OtherVectorType, typename OtherArrayType>
75 Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {}
76
77 /**
78 * \brief Copy constructor for splines.
79 * \param spline The input spline.
80 **/
81 template <int OtherDegree>
82 Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) :
83 m_knots(spline.knots()), m_ctrls(spline.ctrls()) {}
84
85 /**
86 * \brief Returns the knots of the underlying spline.
87 **/
88 const KnotVectorType& knots() const { return m_knots; }
89
90 /**
91 * \brief Returns the knots of the underlying spline.
92 **/
93 const ControlPointVectorType& ctrls() const { return m_ctrls; }
94
95 /**
96 * \brief Returns the spline value at a given site \f$u\f$.
97 *
98 * The function returns
99 * \f{align*}
100 * C(u) & = \sum_{i=0}^{n}N_{i,p}P_i
101 * \f}
102 *
103 * \param u Parameter \f$u \in [0;1]\f$ at which the spline is evaluated.
104 * \return The spline value at the given location \f$u\f$.
105 **/
106 PointType operator()(Scalar u) const;
107
108 /**
109 * \brief Evaluation of spline derivatives of up-to given order.
110 *
111 * The function returns
112 * \f{align*}
113 * \frac{d^i}{du^i}C(u) & = \sum_{i=0}^{n} \frac{d^i}{du^i} N_{i,p}(u)P_i
114 * \f}
115 * for i ranging between 0 and order.
116 *
117 * \param u Parameter \f$u \in [0;1]\f$ at which the spline derivative is evaluated.
118 * \param order The order up to which the derivatives are computed.
119 **/
120 typename SplineTraits<Spline>::DerivativeType
121 derivatives(Scalar u, DenseIndex order) const;
122
123 /**
124 * \copydoc Spline::derivatives
125 * Using the template version of this function is more efficieent since
126 * temporary objects are allocated on the stack whenever this is possible.
127 **/
128 template <int DerivativeOrder>
129 typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
130 derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
131
132 /**
133 * \brief Computes the non-zero basis functions at the given site.
134 *
135 * Splines have local support and a point from their image is defined
136 * by exactly \f$p+1\f$ control points \f$P_i\f$ where \f$p\f$ is the
137 * spline degree.
138 *
139 * This function computes the \f$p+1\f$ non-zero basis function values
140 * for a given parameter value \f$u\f$. It returns
141 * \f{align*}{
142 * N_{i,p}(u), \hdots, N_{i+p+1,p}(u)
143 * \f}
144 *
145 * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis functions
146 * are computed.
147 **/
148 typename SplineTraits<Spline>::BasisVectorType
149 basisFunctions(Scalar u) const;
150
151 /**
152 * \brief Computes the non-zero spline basis function derivatives up to given order.
153 *
154 * The function computes
155 * \f{align*}{
156 * \frac{d^i}{du^i} N_{i,p}(u), \hdots, \frac{d^i}{du^i} N_{i+p+1,p}(u)
157 * \f}
158 * with i ranging from 0 up to the specified order.
159 *
160 * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis function
161 * derivatives are computed.
162 * \param order The order up to which the basis function derivatives are computes.
163 **/
164 typename SplineTraits<Spline>::BasisDerivativeType
165 basisFunctionDerivatives(Scalar u, DenseIndex order) const;
166
167 /**
168 * \copydoc Spline::basisFunctionDerivatives
169 * Using the template version of this function is more efficieent since
170 * temporary objects are allocated on the stack whenever this is possible.
171 **/
172 template <int DerivativeOrder>
173 typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
174 basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
175
176 /**
177 * \brief Returns the spline degree.
178 **/
179 DenseIndex degree() const;
180
181 /**
182 * \brief Returns the span within the knot vector in which u is falling.
183 * \param u The site for which the span is determined.
184 **/
185 DenseIndex span(Scalar u) const;
186
187 /**
188 * \brief Computes the spang within the provided knot vector in which u is falling.
189 **/
190 static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots);
191
192 /**
193 * \brief Returns the spline's non-zero basis functions.
194 *
195 * The function computes and returns
196 * \f{align*}{
197 * N_{i,p}(u), \hdots, N_{i+p+1,p}(u)
198 * \f}
199 *
200 * \param u The site at which the basis functions are computed.
201 * \param degree The degree of the underlying spline.
202 * \param knots The underlying spline's knot vector.
203 **/
204 static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
205
206
207 private:
208 KnotVectorType m_knots; /*!< Knot vector. */
209 ControlPointVectorType m_ctrls; /*!< Control points. */
210 };
211
212 template <typename _Scalar, int _Dim, int _Degree>
213 DenseIndex Spline<_Scalar, _Dim, _Degree>::Span(
214 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::Scalar u,
215 DenseIndex degree,
216 const typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::KnotVectorType& knots)
217 {
218 // Piegl & Tiller, "The NURBS Book", A2.1 (p. 68)
219 if (u <= knots(0)) return degree;
220 const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
221 return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 );
222 }
223
224 template <typename _Scalar, int _Dim, int _Degree>
225 typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType
226 Spline<_Scalar, _Dim, _Degree>::BasisFunctions(
227 typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
228 DenseIndex degree,
229 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
230 {
231 typedef typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType BasisVectorType;
232
233 const DenseIndex p = degree;
234 const DenseIndex i = Spline::Span(u, degree, knots);
235
236 const KnotVectorType& U = knots;
237
238 BasisVectorType left(p+1); left(0) = Scalar(0);
239 BasisVectorType right(p+1); right(0) = Scalar(0);
240
241 VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse();
242 VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u;
243
244 BasisVectorType N(1,p+1);
245 N(0) = Scalar(1);
246 for (DenseIndex j=1; j<=p; ++j)
247 {
248 Scalar saved = Scalar(0);
249 for (DenseIndex r=0; r<j; r++)
250 {
251 const Scalar tmp = N(r)/(right(r+1)+left(j-r));
252 N[r] = saved + right(r+1)*tmp;
253 saved = left(j-r)*tmp;
254 }
255 N(j) = saved;
256 }
257 return N;
258 }
259
260 template <typename _Scalar, int _Dim, int _Degree>
261 DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const
262 {
263 if (_Degree == Dynamic)
264 return m_knots.size() - m_ctrls.cols() - 1;
265 else
266 return _Degree;
267 }
268
269 template <typename _Scalar, int _Dim, int _Degree>
270 DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const
271 {
272 return Spline::Span(u, degree(), knots());
273 }
274
275 template <typename _Scalar, int _Dim, int _Degree>
276 typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const
277 {
278 enum { Order = SplineTraits<Spline>::OrderAtCompileTime };
279
280 const DenseIndex span = this->span(u);
281 const DenseIndex p = degree();
282 const BasisVectorType basis_funcs = basisFunctions(u);
283
284 const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
285 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
286 return (ctrl_weights * ctrl_pts).rowwise().sum();
287 }
288
289 /* --------------------------------------------------------------------------------------------- */
290
291 template <typename SplineType, typename DerivativeType>
292 void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
293 {
294 enum { Dimension = SplineTraits<SplineType>::Dimension };
295 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
296 enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
297
298 typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
299 typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
300 typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;
301
302 const DenseIndex p = spline.degree();
303 const DenseIndex span = spline.span(u);
304
305 const DenseIndex n = (std::min)(p, order);
306
307 der.resize(Dimension,n+1);
308
309 // Retrieve the basis function derivatives up to the desired order...
310 const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);
311
312 // ... and perform the linear combinations of the control points.
313 for (DenseIndex der_order=0; der_order<n+1; ++der_order)
314 {
315 const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
316 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
317 der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
318 }
319 }
320
321 template <typename _Scalar, int _Dim, int _Degree>
322 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType
323 Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
324 {
325 typename SplineTraits< Spline >::DerivativeType res;
326 derivativesImpl(*this, u, order, res);
327 return res;
328 }
329
330 template <typename _Scalar, int _Dim, int _Degree>
331 template <int DerivativeOrder>
332 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType
333 Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
334 {
335 typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
336 derivativesImpl(*this, u, order, res);
337 return res;
338 }
339
340 template <typename _Scalar, int _Dim, int _Degree>
341 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType
342 Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const
343 {
344 return Spline::BasisFunctions(u, degree(), knots());
345 }
346
347 /* --------------------------------------------------------------------------------------------- */
348
349 template <typename SplineType, typename DerivativeType>
350 void basisFunctionDerivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_)
351 {
352 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
353
354 typedef typename SplineTraits<SplineType>::Scalar Scalar;
355 typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType;
356 typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType;
357
358 const KnotVectorType& U = spline.knots();
359
360 const DenseIndex p = spline.degree();
361 const DenseIndex span = spline.span(u);
362
363 const DenseIndex n = (std::min)(p, order);
364
365 N_.resize(n+1, p+1);
366
367 BasisVectorType left = BasisVectorType::Zero(p+1);
368 BasisVectorType right = BasisVectorType::Zero(p+1);
369
370 Matrix<Scalar,Order,Order> ndu(p+1,p+1);
371
372 double saved, temp;
373
374 ndu(0,0) = 1.0;
375
376 DenseIndex j;
377 for (j=1; j<=p; ++j)
378 {
379 left[j] = u-U[span+1-j];
380 right[j] = U[span+j]-u;
381 saved = 0.0;
382
383 for (DenseIndex r=0; r<j; ++r)
384 {
385 /* Lower triangle */
386 ndu(j,r) = right[r+1]+left[j-r];
387 temp = ndu(r,j-1)/ndu(j,r);
388 /* Upper triangle */
389 ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp);
390 saved = left[j-r] * temp;
391 }
392
393 ndu(j,j) = static_cast<Scalar>(saved);
394 }
395
396 for (j = p; j>=0; --j)
397 N_(0,j) = ndu(j,p);
398
399 // Compute the derivatives
400 DerivativeType a(n+1,p+1);
401 DenseIndex r=0;
402 for (; r<=p; ++r)
403 {
404 DenseIndex s1,s2;
405 s1 = 0; s2 = 1; // alternate rows in array a
406 a(0,0) = 1.0;
407
408 // Compute the k-th derivative
409 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
410 {
411 double d = 0.0;
412 DenseIndex rk,pk,j1,j2;
413 rk = r-k; pk = p-k;
414
415 if (r>=k)
416 {
417 a(s2,0) = a(s1,0)/ndu(pk+1,rk);
418 d = a(s2,0)*ndu(rk,pk);
419 }
420
421 if (rk>=-1) j1 = 1;
422 else j1 = -rk;
423
424 if (r-1 <= pk) j2 = k-1;
425 else j2 = p-r;
426
427 for (j=j1; j<=j2; ++j)
428 {
429 a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
430 d += a(s2,j)*ndu(rk+j,pk);
431 }
432
433 if (r<=pk)
434 {
435 a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
436 d += a(s2,k)*ndu(r,pk);
437 }
438
439 N_(k,r) = static_cast<Scalar>(d);
440 j = s1; s1 = s2; s2 = j; // Switch rows
441 }
442 }
443
444 /* Multiply through by the correct factors */
445 /* (Eq. [2.9]) */
446 r = p;
447 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
448 {
449 for (DenseIndex j=p; j>=0; --j) N_(k,j) *= r;
450 r *= p-k;
451 }
452 }
453
454 template <typename _Scalar, int _Dim, int _Degree>
455 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
456 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
457 {
458 typename SplineTraits< Spline >::BasisDerivativeType der;
459 basisFunctionDerivativesImpl(*this, u, order, der);
460 return der;
461 }
462
463 template <typename _Scalar, int _Dim, int _Degree>
464 template <int DerivativeOrder>
465 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
466 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
467 {
468 typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der;
469 basisFunctionDerivativesImpl(*this, u, order, der);
470 return der;
471 }
472}
473
474#endif // EIGEN_SPLINE_H
Note: See TracBrowser for help on using the repository browser.