Spline.h (18307B)
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 { 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 parameter vectors. */ 49 typedef typename SplineTraits<Spline>::ParameterVectorType ParameterVectorType; 50 51 /** \brief The data type used to store non-zero basis functions. */ 52 typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType; 53 54 /** \brief The data type used to store the values of the basis function derivatives. */ 55 typedef typename SplineTraits<Spline>::BasisDerivativeType BasisDerivativeType; 56 57 /** \brief The data type representing the spline's control points. */ 58 typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType; 59 60 /** 61 * \brief Creates a (constant) zero spline. 62 * For Splines with dynamic degree, the resulting degree will be 0. 63 **/ 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 75 /** 76 * \brief Creates a spline from a knot vector and control points. 77 * \param knots The spline's knot vector. 78 * \param ctrls The spline's control point vector. 79 **/ 80 template <typename OtherVectorType, typename OtherArrayType> 81 Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {} 82 83 /** 84 * \brief Copy constructor for splines. 85 * \param spline The input spline. 86 **/ 87 template <int OtherDegree> 88 Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) : 89 m_knots(spline.knots()), m_ctrls(spline.ctrls()) {} 90 91 /** 92 * \brief Returns the knots of the underlying spline. 93 **/ 94 const KnotVectorType& knots() const { return m_knots; } 95 96 /** 97 * \brief Returns the ctrls of the underlying spline. 98 **/ 99 const ControlPointVectorType& ctrls() const { return m_ctrls; } 100 101 /** 102 * \brief Returns the spline value at a given site \f$u\f$. 103 * 104 * The function returns 105 * \f{align*} 106 * C(u) & = \sum_{i=0}^{n}N_{i,p}P_i 107 * \f} 108 * 109 * \param u Parameter \f$u \in [0;1]\f$ at which the spline is evaluated. 110 * \return The spline value at the given location \f$u\f$. 111 **/ 112 PointType operator()(Scalar u) const; 113 114 /** 115 * \brief Evaluation of spline derivatives of up-to given order. 116 * 117 * The function returns 118 * \f{align*} 119 * \frac{d^i}{du^i}C(u) & = \sum_{i=0}^{n} \frac{d^i}{du^i} N_{i,p}(u)P_i 120 * \f} 121 * for i ranging between 0 and order. 122 * 123 * \param u Parameter \f$u \in [0;1]\f$ at which the spline derivative is evaluated. 124 * \param order The order up to which the derivatives are computed. 125 **/ 126 typename SplineTraits<Spline>::DerivativeType 127 derivatives(Scalar u, DenseIndex order) const; 128 129 /** 130 * \copydoc Spline::derivatives 131 * Using the template version of this function is more efficieent since 132 * temporary objects are allocated on the stack whenever this is possible. 133 **/ 134 template <int DerivativeOrder> 135 typename SplineTraits<Spline,DerivativeOrder>::DerivativeType 136 derivatives(Scalar u, DenseIndex order = DerivativeOrder) const; 137 138 /** 139 * \brief Computes the non-zero basis functions at the given site. 140 * 141 * Splines have local support and a point from their image is defined 142 * by exactly \f$p+1\f$ control points \f$P_i\f$ where \f$p\f$ is the 143 * spline degree. 144 * 145 * This function computes the \f$p+1\f$ non-zero basis function values 146 * for a given parameter value \f$u\f$. It returns 147 * \f{align*}{ 148 * N_{i,p}(u), \hdots, N_{i+p+1,p}(u) 149 * \f} 150 * 151 * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis functions 152 * are computed. 153 **/ 154 typename SplineTraits<Spline>::BasisVectorType 155 basisFunctions(Scalar u) const; 156 157 /** 158 * \brief Computes the non-zero spline basis function derivatives up to given order. 159 * 160 * The function computes 161 * \f{align*}{ 162 * \frac{d^i}{du^i} N_{i,p}(u), \hdots, \frac{d^i}{du^i} N_{i+p+1,p}(u) 163 * \f} 164 * with i ranging from 0 up to the specified order. 165 * 166 * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis function 167 * derivatives are computed. 168 * \param order The order up to which the basis function derivatives are computes. 169 **/ 170 typename SplineTraits<Spline>::BasisDerivativeType 171 basisFunctionDerivatives(Scalar u, DenseIndex order) const; 172 173 /** 174 * \copydoc Spline::basisFunctionDerivatives 175 * Using the template version of this function is more efficieent since 176 * temporary objects are allocated on the stack whenever this is possible. 177 **/ 178 template <int DerivativeOrder> 179 typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType 180 basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const; 181 182 /** 183 * \brief Returns the spline degree. 184 **/ 185 DenseIndex degree() const; 186 187 /** 188 * \brief Returns the span within the knot vector in which u is falling. 189 * \param u The site for which the span is determined. 190 **/ 191 DenseIndex span(Scalar u) const; 192 193 /** 194 * \brief Computes the span within the provided knot vector in which u is falling. 195 **/ 196 static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots); 197 198 /** 199 * \brief Returns the spline's non-zero basis functions. 200 * 201 * The function computes and returns 202 * \f{align*}{ 203 * N_{i,p}(u), \hdots, N_{i+p+1,p}(u) 204 * \f} 205 * 206 * \param u The site at which the basis functions are computed. 207 * \param degree The degree of the underlying spline. 208 * \param knots The underlying spline's knot vector. 209 **/ 210 static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots); 211 212 /** 213 * \copydoc Spline::basisFunctionDerivatives 214 * \param degree The degree of the underlying spline 215 * \param knots The underlying spline's knot vector. 216 **/ 217 static BasisDerivativeType BasisFunctionDerivatives( 218 const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType& knots); 219 220 private: 221 KnotVectorType m_knots; /*!< Knot vector. */ 222 ControlPointVectorType m_ctrls; /*!< Control points. */ 223 224 template <typename DerivativeType> 225 static void BasisFunctionDerivativesImpl( 226 const typename Spline<_Scalar, _Dim, _Degree>::Scalar u, 227 const DenseIndex order, 228 const DenseIndex p, 229 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U, 230 DerivativeType& N_); 231 }; 232 233 template <typename _Scalar, int _Dim, int _Degree> 234 DenseIndex Spline<_Scalar, _Dim, _Degree>::Span( 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> 246 typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType 247 Spline<_Scalar, _Dim, _Degree>::BasisFunctions( 248 typename Spline<_Scalar, _Dim, _Degree>::Scalar u, 249 DenseIndex degree, 250 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots) 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 260 VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse(); 261 VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u; 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> 280 DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const 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> 289 DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const 290 { 291 return Spline::Span(u, degree(), knots()); 292 } 293 294 template <typename _Scalar, int _Dim, int _Degree> 295 typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const 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 342 Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const 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 361 Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const 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> 371 void Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivativesImpl( 372 const typename Spline<_Scalar, _Dim, _Degree>::Scalar u, 373 const DenseIndex order, 374 const DenseIndex p, 375 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U, 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 476 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const 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 495 Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivatives( 496 const typename Spline<_Scalar, _Dim, _Degree>::Scalar u, 497 const DenseIndex order, 498 const DenseIndex degree, 499 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots) 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