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 (at) 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 spang 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 typedef typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType BasisVectorType; 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 262 VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse(); 263 VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u; 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> 282 DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const 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> 291 DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const 292 { 293 return Spline::Span(u, degree(), knots()); 294 } 295 296 template <typename _Scalar, int _Dim, int _Degree> 297 typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const 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 344 Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const 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 363 Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const 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> 373 void Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivativesImpl( 374 const typename Spline<_Scalar, _Dim, _Degree>::Scalar u, 375 const DenseIndex order, 376 const DenseIndex p, 377 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U, 378 DerivativeType& N_) 379 { 380 typedef Spline<_Scalar, _Dim, _Degree> SplineType; 381 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime }; 382 383 typedef typename SplineTraits<SplineType>::Scalar Scalar; 384 typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType; 385 386 const DenseIndex span = SplineType::Span(u, p, U); 387 388 const DenseIndex n = (std::min)(p, order); 389 390 N_.resize(n+1, p+1); 391 392 BasisVectorType left = BasisVectorType::Zero(p+1); 393 BasisVectorType right = BasisVectorType::Zero(p+1); 394 395 Matrix<Scalar,Order,Order> ndu(p+1,p+1); 396 397 Scalar saved, temp; // FIXME These were double instead of Scalar. Was there a reason for that? 398 399 ndu(0,0) = 1.0; 400 401 DenseIndex j; 402 for (j=1; j<=p; ++j) 403 { 404 left[j] = u-U[span+1-j]; 405 right[j] = U[span+j]-u; 406 saved = 0.0; 407 408 for (DenseIndex r=0; r<j; ++r) 409 { 410 /* Lower triangle */ 411 ndu(j,r) = right[r+1]+left[j-r]; 412 temp = ndu(r,j-1)/ndu(j,r); 413 /* Upper triangle */ 414 ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp); 415 saved = left[j-r] * temp; 416 } 417 418 ndu(j,j) = static_cast<Scalar>(saved); 419 } 420 421 for (j = p; j>=0; --j) 422 N_(0,j) = ndu(j,p); 423 424 // Compute the derivatives 425 DerivativeType a(n+1,p+1); 426 DenseIndex r=0; 427 for (; r<=p; ++r) 428 { 429 DenseIndex s1,s2; 430 s1 = 0; s2 = 1; // alternate rows in array a 431 a(0,0) = 1.0; 432 433 // Compute the k-th derivative 434 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k) 435 { 436 Scalar d = 0.0; 437 DenseIndex rk,pk,j1,j2; 438 rk = r-k; pk = p-k; 439 440 if (r>=k) 441 { 442 a(s2,0) = a(s1,0)/ndu(pk+1,rk); 443 d = a(s2,0)*ndu(rk,pk); 444 } 445 446 if (rk>=-1) j1 = 1; 447 else j1 = -rk; 448 449 if (r-1 <= pk) j2 = k-1; 450 else j2 = p-r; 451 452 for (j=j1; j<=j2; ++j) 453 { 454 a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j); 455 d += a(s2,j)*ndu(rk+j,pk); 456 } 457 458 if (r<=pk) 459 { 460 a(s2,k) = -a(s1,k-1)/ndu(pk+1,r); 461 d += a(s2,k)*ndu(r,pk); 462 } 463 464 N_(k,r) = static_cast<Scalar>(d); 465 j = s1; s1 = s2; s2 = j; // Switch rows 466 } 467 } 468 469 /* Multiply through by the correct factors */ 470 /* (Eq. [2.9]) */ 471 r = p; 472 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k) 473 { 474 for (j=p; j>=0; --j) N_(k,j) *= r; 475 r *= p-k; 476 } 477 } 478 479 template <typename _Scalar, int _Dim, int _Degree> 480 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType 481 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const 482 { 483 typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType der; 484 BasisFunctionDerivativesImpl(u, order, degree(), knots(), der); 485 return der; 486 } 487 488 template <typename _Scalar, int _Dim, int _Degree> 489 template <int DerivativeOrder> 490 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType 491 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const 492 { 493 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType der; 494 BasisFunctionDerivativesImpl(u, order, degree(), knots(), der); 495 return der; 496 } 497 498 template <typename _Scalar, int _Dim, int _Degree> 499 typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType 500 Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivatives( 501 const typename Spline<_Scalar, _Dim, _Degree>::Scalar u, 502 const DenseIndex order, 503 const DenseIndex degree, 504 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots) 505 { 506 typename SplineTraits<Spline>::BasisDerivativeType der; 507 BasisFunctionDerivativesImpl(u, order, degree, knots, der); 508 return der; 509 } 510 } 511 512 #endif // EIGEN_SPLINE_H 513
