File indexing completed on 2025-01-18 09:57:10
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #ifndef EIGEN_SPLINE_H
0011 #define EIGEN_SPLINE_H
0012
0013 #include "SplineFwd.h"
0014
0015 namespace Eigen
0016 {
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034 template <typename _Scalar, int _Dim, int _Degree>
0035 class Spline
0036 {
0037 public:
0038 typedef _Scalar Scalar;
0039 enum { Dimension = _Dim };
0040 enum { Degree = _Degree };
0041
0042
0043 typedef typename SplineTraits<Spline>::PointType PointType;
0044
0045
0046 typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;
0047
0048
0049 typedef typename SplineTraits<Spline>::ParameterVectorType ParameterVectorType;
0050
0051
0052 typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;
0053
0054
0055 typedef typename SplineTraits<Spline>::BasisDerivativeType BasisDerivativeType;
0056
0057
0058 typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;
0059
0060
0061
0062
0063
0064 Spline()
0065 : m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2))
0066 , m_ctrls(ControlPointVectorType::Zero(Dimension,(Degree==Dynamic ? 1 : Degree+1)))
0067 {
0068
0069
0070 enum { MinDegree = (Degree==Dynamic ? 0 : Degree) };
0071 m_knots.template segment<MinDegree+1>(0) = Array<Scalar,1,MinDegree+1>::Zero();
0072 m_knots.template segment<MinDegree+1>(MinDegree+1) = Array<Scalar,1,MinDegree+1>::Ones();
0073 }
0074
0075
0076
0077
0078
0079
0080 template <typename OtherVectorType, typename OtherArrayType>
0081 Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {}
0082
0083
0084
0085
0086
0087 template <int OtherDegree>
0088 Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) :
0089 m_knots(spline.knots()), m_ctrls(spline.ctrls()) {}
0090
0091
0092
0093
0094 const KnotVectorType& knots() const { return m_knots; }
0095
0096
0097
0098
0099 const ControlPointVectorType& ctrls() const { return m_ctrls; }
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112 PointType operator()(Scalar u) const;
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126 typename SplineTraits<Spline>::DerivativeType
0127 derivatives(Scalar u, DenseIndex order) const;
0128
0129
0130
0131
0132
0133
0134 template <int DerivativeOrder>
0135 typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
0136 derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154 typename SplineTraits<Spline>::BasisVectorType
0155 basisFunctions(Scalar u) const;
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170 typename SplineTraits<Spline>::BasisDerivativeType
0171 basisFunctionDerivatives(Scalar u, DenseIndex order) const;
0172
0173
0174
0175
0176
0177
0178 template <int DerivativeOrder>
0179 typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
0180 basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;
0181
0182
0183
0184
0185 DenseIndex degree() const;
0186
0187
0188
0189
0190
0191 DenseIndex span(Scalar u) const;
0192
0193
0194
0195
0196 static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots);
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210 static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);
0211
0212
0213
0214
0215
0216
0217 static BasisDerivativeType BasisFunctionDerivatives(
0218 const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType& knots);
0219
0220 private:
0221 KnotVectorType m_knots;
0222 ControlPointVectorType m_ctrls;
0223
0224 template <typename DerivativeType>
0225 static void BasisFunctionDerivativesImpl(
0226 const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
0227 const DenseIndex order,
0228 const DenseIndex p,
0229 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U,
0230 DerivativeType& N_);
0231 };
0232
0233 template <typename _Scalar, int _Dim, int _Degree>
0234 DenseIndex Spline<_Scalar, _Dim, _Degree>::Span(
0235 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::Scalar u,
0236 DenseIndex degree,
0237 const typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::KnotVectorType& knots)
0238 {
0239
0240 if (u <= knots(0)) return degree;
0241 const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
0242 return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 );
0243 }
0244
0245 template <typename _Scalar, int _Dim, int _Degree>
0246 typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType
0247 Spline<_Scalar, _Dim, _Degree>::BasisFunctions(
0248 typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
0249 DenseIndex degree,
0250 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
0251 {
0252 const DenseIndex p = degree;
0253 const DenseIndex i = Spline::Span(u, degree, knots);
0254
0255 const KnotVectorType& U = knots;
0256
0257 BasisVectorType left(p+1); left(0) = Scalar(0);
0258 BasisVectorType right(p+1); right(0) = Scalar(0);
0259
0260 VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse();
0261 VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u;
0262
0263 BasisVectorType N(1,p+1);
0264 N(0) = Scalar(1);
0265 for (DenseIndex j=1; j<=p; ++j)
0266 {
0267 Scalar saved = Scalar(0);
0268 for (DenseIndex r=0; r<j; r++)
0269 {
0270 const Scalar tmp = N(r)/(right(r+1)+left(j-r));
0271 N[r] = saved + right(r+1)*tmp;
0272 saved = left(j-r)*tmp;
0273 }
0274 N(j) = saved;
0275 }
0276 return N;
0277 }
0278
0279 template <typename _Scalar, int _Dim, int _Degree>
0280 DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const
0281 {
0282 if (_Degree == Dynamic)
0283 return m_knots.size() - m_ctrls.cols() - 1;
0284 else
0285 return _Degree;
0286 }
0287
0288 template <typename _Scalar, int _Dim, int _Degree>
0289 DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const
0290 {
0291 return Spline::Span(u, degree(), knots());
0292 }
0293
0294 template <typename _Scalar, int _Dim, int _Degree>
0295 typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const
0296 {
0297 enum { Order = SplineTraits<Spline>::OrderAtCompileTime };
0298
0299 const DenseIndex span = this->span(u);
0300 const DenseIndex p = degree();
0301 const BasisVectorType basis_funcs = basisFunctions(u);
0302
0303 const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
0304 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
0305 return (ctrl_weights * ctrl_pts).rowwise().sum();
0306 }
0307
0308
0309
0310 template <typename SplineType, typename DerivativeType>
0311 void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
0312 {
0313 enum { Dimension = SplineTraits<SplineType>::Dimension };
0314 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
0315 enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };
0316
0317 typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
0318 typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
0319 typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;
0320
0321 const DenseIndex p = spline.degree();
0322 const DenseIndex span = spline.span(u);
0323
0324 const DenseIndex n = (std::min)(p, order);
0325
0326 der.resize(Dimension,n+1);
0327
0328
0329 const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);
0330
0331
0332 for (DenseIndex der_order=0; der_order<n+1; ++der_order)
0333 {
0334 const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
0335 const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
0336 der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
0337 }
0338 }
0339
0340 template <typename _Scalar, int _Dim, int _Degree>
0341 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType
0342 Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
0343 {
0344 typename SplineTraits< Spline >::DerivativeType res;
0345 derivativesImpl(*this, u, order, res);
0346 return res;
0347 }
0348
0349 template <typename _Scalar, int _Dim, int _Degree>
0350 template <int DerivativeOrder>
0351 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType
0352 Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
0353 {
0354 typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
0355 derivativesImpl(*this, u, order, res);
0356 return res;
0357 }
0358
0359 template <typename _Scalar, int _Dim, int _Degree>
0360 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType
0361 Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const
0362 {
0363 return Spline::BasisFunctions(u, degree(), knots());
0364 }
0365
0366
0367
0368
0369 template <typename _Scalar, int _Dim, int _Degree>
0370 template <typename DerivativeType>
0371 void Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivativesImpl(
0372 const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
0373 const DenseIndex order,
0374 const DenseIndex p,
0375 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U,
0376 DerivativeType& N_)
0377 {
0378 typedef Spline<_Scalar, _Dim, _Degree> SplineType;
0379 enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
0380
0381 const DenseIndex span = SplineType::Span(u, p, U);
0382
0383 const DenseIndex n = (std::min)(p, order);
0384
0385 N_.resize(n+1, p+1);
0386
0387 BasisVectorType left = BasisVectorType::Zero(p+1);
0388 BasisVectorType right = BasisVectorType::Zero(p+1);
0389
0390 Matrix<Scalar,Order,Order> ndu(p+1,p+1);
0391
0392 Scalar saved, temp;
0393
0394 ndu(0,0) = 1.0;
0395
0396 DenseIndex j;
0397 for (j=1; j<=p; ++j)
0398 {
0399 left[j] = u-U[span+1-j];
0400 right[j] = U[span+j]-u;
0401 saved = 0.0;
0402
0403 for (DenseIndex r=0; r<j; ++r)
0404 {
0405
0406 ndu(j,r) = right[r+1]+left[j-r];
0407 temp = ndu(r,j-1)/ndu(j,r);
0408
0409 ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp);
0410 saved = left[j-r] * temp;
0411 }
0412
0413 ndu(j,j) = static_cast<Scalar>(saved);
0414 }
0415
0416 for (j = p; j>=0; --j)
0417 N_(0,j) = ndu(j,p);
0418
0419
0420 DerivativeType a(n+1,p+1);
0421 DenseIndex r=0;
0422 for (; r<=p; ++r)
0423 {
0424 DenseIndex s1,s2;
0425 s1 = 0; s2 = 1;
0426 a(0,0) = 1.0;
0427
0428
0429 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
0430 {
0431 Scalar d = 0.0;
0432 DenseIndex rk,pk,j1,j2;
0433 rk = r-k; pk = p-k;
0434
0435 if (r>=k)
0436 {
0437 a(s2,0) = a(s1,0)/ndu(pk+1,rk);
0438 d = a(s2,0)*ndu(rk,pk);
0439 }
0440
0441 if (rk>=-1) j1 = 1;
0442 else j1 = -rk;
0443
0444 if (r-1 <= pk) j2 = k-1;
0445 else j2 = p-r;
0446
0447 for (j=j1; j<=j2; ++j)
0448 {
0449 a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
0450 d += a(s2,j)*ndu(rk+j,pk);
0451 }
0452
0453 if (r<=pk)
0454 {
0455 a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
0456 d += a(s2,k)*ndu(r,pk);
0457 }
0458
0459 N_(k,r) = static_cast<Scalar>(d);
0460 j = s1; s1 = s2; s2 = j;
0461 }
0462 }
0463
0464
0465
0466 r = p;
0467 for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
0468 {
0469 for (j=p; j>=0; --j) N_(k,j) *= r;
0470 r *= p-k;
0471 }
0472 }
0473
0474 template <typename _Scalar, int _Dim, int _Degree>
0475 typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
0476 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
0477 {
0478 typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType der;
0479 BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
0480 return der;
0481 }
0482
0483 template <typename _Scalar, int _Dim, int _Degree>
0484 template <int DerivativeOrder>
0485 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
0486 Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
0487 {
0488 typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType der;
0489 BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
0490 return der;
0491 }
0492
0493 template <typename _Scalar, int _Dim, int _Degree>
0494 typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
0495 Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivatives(
0496 const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
0497 const DenseIndex order,
0498 const DenseIndex degree,
0499 const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
0500 {
0501 typename SplineTraits<Spline>::BasisDerivativeType der;
0502 BasisFunctionDerivativesImpl(u, order, degree, knots, der);
0503 return der;
0504 }
0505 }
0506
0507 #endif