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 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 475
