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 class 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 spline from a knot vector and control points. 56 * \param knots The spline's knot vector. 57 * \param ctrls The spline's control point vector. 58 **/ 59 template <typename OtherVectorType, typename OtherArrayType> 60 Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {} 61 62 /** 63 * \brief Copy constructor for splines. 64 * \param spline The input spline. 65 **/ 66 template <int OtherDegree> 67 Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) : 68 m_knots(spline.knots()), m_ctrls(spline.ctrls()) {} 69 70 /** 71 * \brief Returns the knots of the underlying spline. 72 **/ 73 const KnotVectorType& knots() const { return m_knots; } 74 75 /** 76 * \brief Returns the knots of the underlying spline. 77 **/ 78 const ControlPointVectorType& ctrls() const { return m_ctrls; } 79 80 /** 81 * \brief Returns the spline value at a given site \f$u\f$. 82 * 83 * The function returns 84 * \f{align*} 85 * C(u) & = \sum_{i=0}^{n}N_{i,p}P_i 86 * \f} 87 * 88 * \param u Parameter \f$u \in [0;1]\f$ at which the spline is evaluated. 89 * \return The spline value at the given location \f$u\f$. 90 **/ 91 PointType operator()(Scalar u) const; 92 93 /** 94 * \brief Evaluation of spline derivatives of up-to given order. 95 * 96 * The function returns 97 * \f{align*} 98 * \frac{d^i}{du^i}C(u) & = \sum_{i=0}^{n} \frac{d^i}{du^i} N_{i,p}(u)P_i 99 * \f} 100 * for i ranging between 0 and order. 101 * 102 * \param u Parameter \f$u \in [0;1]\f$ at which the spline derivative is evaluated. 103 * \param order The order up to which the derivatives are computed. 104 **/ 105 typename SplineTraits<Spline>::DerivativeType 106 derivatives(Scalar u, DenseIndex order) const; 107 108 /** 109 * \copydoc Spline::derivatives 110 * Using the template version of this function is more efficieent since 111 * temporary objects are allocated on the stack whenever this is possible. 112 **/ 113 template <int DerivativeOrder> 114 typename SplineTraits<Spline,DerivativeOrder>::DerivativeType 115 derivatives(Scalar u, DenseIndex order = DerivativeOrder) const; 116 117 /** 118 * \brief Computes the non-zero basis functions at the given site. 119 * 120 * Splines have local support and a point from their image is defined 121 * by exactly \f$p+1\f$ control points \f$P_i\f$ where \f$p\f$ is the 122 * spline degree. 123 * 124 * This function computes the \f$p+1\f$ non-zero basis function values 125 * for a given parameter value \f$u\f$. It returns 126 * \f{align*}{ 127 * N_{i,p}(u), \hdots, N_{i+p+1,p}(u) 128 * \f} 129 * 130 * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis functions 131 * are computed. 132 **/ 133 typename SplineTraits<Spline>::BasisVectorType 134 basisFunctions(Scalar u) const; 135 136 /** 137 * \brief Computes the non-zero spline basis function derivatives up to given order. 138 * 139 * The function computes 140 * \f{align*}{ 141 * \frac{d^i}{du^i} N_{i,p}(u), \hdots, \frac{d^i}{du^i} N_{i+p+1,p}(u) 142 * \f} 143 * with i ranging from 0 up to the specified order. 144 * 145 * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis function 146 * derivatives are computed. 147 * \param order The order up to which the basis function derivatives are computes. 148 **/ 149 typename SplineTraits<Spline>::BasisDerivativeType 150 basisFunctionDerivatives(Scalar u, DenseIndex order) const; 151 152 /** 153 * \copydoc Spline::basisFunctionDerivatives 154 * Using the template version of this function is more efficieent since 155 * temporary objects are allocated on the stack whenever this is possible. 156 **/ 157 template <int DerivativeOrder> 158 typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType 159 basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const; 160 161 /** 162 * \brief Returns the spline degree. 163 **/ 164 DenseIndex degree() const; 165 166 /** 167 * \brief Returns the span within the knot vector in which u is falling. 168 * \param u The site for which the span is determined. 169 **/ 170 DenseIndex span(Scalar u) const; 171 172 /** 173 * \brief Computes the spang within the provided knot vector in which u is falling. 174 **/ 175 static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots); 176 177 /** 178 * \brief Returns the spline's non-zero basis functions. 179 * 180 * The function computes and returns 181 * \f{align*}{ 182 * N_{i,p}(u), \hdots, N_{i+p+1,p}(u) 183 * \f} 184 * 185 * \param u The site at which the basis functions are computed. 186 * \param degree The degree of the underlying spline. 187 * \param knots The underlying spline's knot vector. 188 **/ 189 static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots); 190 191 192 private: 193 KnotVectorType m_knots; /*!< Knot vector. */ 194 ControlPointVectorType m_ctrls; /*!< Control points. */ 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 331 Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const 332 { 333 return Spline::BasisFunctions(u, degree(), knots()); 334 } 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 465
