1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_SPLINE_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_SPLINE_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "SplineFwd.h" 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \ingroup Splines_Module 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class Spline 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A class representing multi-dimensional spline curves. 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The class represents B-splines with non-uniform knot vectors. Each control 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * point of the B-spline is associated with a basis function 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f{align*} 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * C(u) & = \sum_{i=0}^{n}N_{i,p}(u)P_i 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f} 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _Scalar The underlying data type (typically float or double) 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _Dim The curve dimension (e.g. 2 or 3) 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _Degree Per default set to Dynamic; could be set to the actual desired 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * degree for optimization purposes (would result in stack allocation 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of several temporary variables). 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath class Spline 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Scalar Scalar; /*!< The spline curve's scalar type. */ 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Dimension = _Dim /*!< The spline curve's dimension. */ }; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Degree = _Degree /*!< The spline curve's degree. */ }; 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief The point type the spline is representing. */ 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<Spline>::PointType PointType; 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief The data type used to store knot vectors. */ 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief The data type used to store non-zero basis functions. */ 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType; 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief The data type representing the spline's control points. */ 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType; 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Creates a (constant) zero spline. 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * For Splines with dynamic degree, the resulting degree will be 0. 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez **/ 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Spline() 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2)) 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez , m_ctrls(ControlPointVectorType::Zero(2,(Degree==Dynamic ? 1 : Degree+1))) 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // in theory this code can go to the initializer list but it will get pretty 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // much unreadable ... 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { MinDegree = (Degree==Dynamic ? 0 : Degree) }; 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_knots.template segment<MinDegree+1>(0) = Array<Scalar,1,MinDegree+1>::Zero(); 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_knots.template segment<MinDegree+1>(MinDegree+1) = Array<Scalar,1,MinDegree+1>::Ones(); 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Creates a spline from a knot vector and control points. 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param knots The spline's knot vector. 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param ctrls The spline's control point vector. 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename OtherVectorType, typename OtherArrayType> 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {} 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Copy constructor for splines. 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param spline The input spline. 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <int OtherDegree> 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) : 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_knots(spline.knots()), m_ctrls(spline.ctrls()) {} 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Returns the knots of the underlying spline. 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const KnotVectorType& knots() const { return m_knots; } 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Returns the knots of the underlying spline. 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const ControlPointVectorType& ctrls() const { return m_ctrls; } 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Returns the spline value at a given site \f$u\f$. 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The function returns 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f{align*} 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * C(u) & = \sum_{i=0}^{n}N_{i,p}P_i 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f} 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param u Parameter \f$u \in [0;1]\f$ at which the spline is evaluated. 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \return The spline value at the given location \f$u\f$. 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PointType operator()(Scalar u) const; 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Evaluation of spline derivatives of up-to given order. 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The function returns 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f{align*} 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \frac{d^i}{du^i}C(u) & = \sum_{i=0}^{n} \frac{d^i}{du^i} N_{i,p}(u)P_i 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f} 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * for i ranging between 0 and order. 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param u Parameter \f$u \in [0;1]\f$ at which the spline derivative is evaluated. 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param order The order up to which the derivatives are computed. 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits<Spline>::DerivativeType 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derivatives(Scalar u, DenseIndex order) const; 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \copydoc Spline::derivatives 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Using the template version of this function is more efficieent since 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * temporary objects are allocated on the stack whenever this is possible. 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <int DerivativeOrder> 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits<Spline,DerivativeOrder>::DerivativeType 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derivatives(Scalar u, DenseIndex order = DerivativeOrder) const; 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Computes the non-zero basis functions at the given site. 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Splines have local support and a point from their image is defined 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * by exactly \f$p+1\f$ control points \f$P_i\f$ where \f$p\f$ is the 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * spline degree. 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function computes the \f$p+1\f$ non-zero basis function values 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * for a given parameter value \f$u\f$. It returns 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f{align*}{ 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * N_{i,p}(u), \hdots, N_{i+p+1,p}(u) 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f} 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis functions 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * are computed. 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits<Spline>::BasisVectorType 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath basisFunctions(Scalar u) const; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Computes the non-zero spline basis function derivatives up to given order. 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The function computes 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f{align*}{ 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \frac{d^i}{du^i} N_{i,p}(u), \hdots, \frac{d^i}{du^i} N_{i+p+1,p}(u) 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f} 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with i ranging from 0 up to the specified order. 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis function 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * derivatives are computed. 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param order The order up to which the basis function derivatives are computes. 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits<Spline>::BasisDerivativeType 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath basisFunctionDerivatives(Scalar u, DenseIndex order) const; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \copydoc Spline::basisFunctionDerivatives 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Using the template version of this function is more efficieent since 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * temporary objects are allocated on the stack whenever this is possible. 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <int DerivativeOrder> 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const; 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Returns the spline degree. 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex degree() const; 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Returns the span within the knot vector in which u is falling. 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param u The site for which the span is determined. 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex span(Scalar u) const; 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Computes the spang within the provided knot vector in which u is falling. 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots); 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Returns the spline's non-zero basis functions. 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The function computes and returns 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f{align*}{ 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * N_{i,p}(u), \hdots, N_{i+p+1,p}(u) 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f} 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param u The site at which the basis functions are computed. 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param degree The degree of the underlying spline. 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param knots The underlying spline's knot vector. 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath **/ 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath KnotVectorType m_knots; /*!< Knot vector. */ 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ControlPointVectorType m_ctrls; /*!< Control points. */ 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex Spline<_Scalar, _Dim, _Degree>::Span( 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::Scalar u, 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex degree, 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::KnotVectorType& knots) 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Piegl & Tiller, "The NURBS Book", A2.1 (p. 68) 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (u <= knots(0)) return degree; 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u); 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 ); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Spline<_Scalar, _Dim, _Degree>::BasisFunctions( 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Spline<_Scalar, _Dim, _Degree>::Scalar u, 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex degree, 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots) 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType BasisVectorType; 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex p = degree; 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex i = Spline::Span(u, degree, knots); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const KnotVectorType& U = knots; 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BasisVectorType left(p+1); left(0) = Scalar(0); 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BasisVectorType right(p+1); right(0) = Scalar(0); 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse(); 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u; 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BasisVectorType N(1,p+1); 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath N(0) = Scalar(1); 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (DenseIndex j=1; j<=p; ++j) 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar saved = Scalar(0); 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (DenseIndex r=0; r<j; r++) 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar tmp = N(r)/(right(r+1)+left(j-r)); 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath N[r] = saved + right(r+1)*tmp; 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath saved = left(j-r)*tmp; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath N(j) = saved; 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return N; 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (_Degree == Dynamic) 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_knots.size() - m_ctrls.cols() - 1; 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return _Degree; 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Spline::Span(u, degree(), knots()); 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Order = SplineTraits<Spline>::OrderAtCompileTime }; 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex span = this->span(u); 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex p = degree(); 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const BasisVectorType basis_funcs = basisFunctions(u); 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs); 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1); 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return (ctrl_weights * ctrl_pts).rowwise().sum(); 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* --------------------------------------------------------------------------------------------- */ 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename SplineType, typename DerivativeType> 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der) 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Dimension = SplineTraits<SplineType>::Dimension }; 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Order = SplineTraits<SplineType>::OrderAtCompileTime }; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { DerivativeOrder = DerivativeType::ColsAtCompileTime }; 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType; 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType; 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr; 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex p = spline.degree(); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex span = spline.span(u); 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex n = (std::min)(p, order); 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath der.resize(Dimension,n+1); 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Retrieve the basis function derivatives up to the desired order... 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // ... and perform the linear combinations of the control points. 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (DenseIndex der_order=0; der_order<n+1; ++der_order) 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) ); 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1); 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum(); 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline >::DerivativeType res; 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derivativesImpl(*this, u, order, res); 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return res; 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <int DerivativeOrder> 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res; 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath derivativesImpl(*this, u, order, res); 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return res; 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Spline::BasisFunctions(u, degree(), knots()); 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* --------------------------------------------------------------------------------------------- */ 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename SplineType, typename DerivativeType> 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void basisFunctionDerivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_) 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Order = SplineTraits<SplineType>::OrderAtCompileTime }; 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<SplineType>::Scalar Scalar; 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType; 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType; 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const KnotVectorType& U = spline.knots(); 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex p = spline.degree(); 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex span = spline.span(u); 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const DenseIndex n = (std::min)(p, order); 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath N_.resize(n+1, p+1); 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BasisVectorType left = BasisVectorType::Zero(p+1); 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BasisVectorType right = BasisVectorType::Zero(p+1); 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,Order,Order> ndu(p+1,p+1); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double saved, temp; 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ndu(0,0) = 1.0; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex j; 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j=1; j<=p; ++j) 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath left[j] = u-U[span+1-j]; 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath right[j] = U[span+j]-u; 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath saved = 0.0; 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (DenseIndex r=0; r<j; ++r) 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* Lower triangle */ 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ndu(j,r) = right[r+1]+left[j-r]; 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = ndu(r,j-1)/ndu(j,r); 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* Upper triangle */ 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp); 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath saved = left[j-r] * temp; 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ndu(j,j) = static_cast<Scalar>(saved); 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = p; j>=0; --j) 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath N_(0,j) = ndu(j,p); 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the derivatives 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DerivativeType a(n+1,p+1); 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex r=0; 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (; r<=p; ++r) 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex s1,s2; 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath s1 = 0; s2 = 1; // alternate rows in array a 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a(0,0) = 1.0; 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the k-th derivative 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k) 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double d = 0.0; 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex rk,pk,j1,j2; 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rk = r-k; pk = p-k; 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (r>=k) 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a(s2,0) = a(s1,0)/ndu(pk+1,rk); 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath d = a(s2,0)*ndu(rk,pk); 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (rk>=-1) j1 = 1; 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else j1 = -rk; 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (r-1 <= pk) j2 = k-1; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else j2 = p-r; 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j=j1; j<=j2; ++j) 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j); 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath d += a(s2,j)*ndu(rk+j,pk); 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (r<=pk) 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a(s2,k) = -a(s1,k-1)/ndu(pk+1,r); 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath d += a(s2,k)*ndu(r,pk); 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath N_(k,r) = static_cast<Scalar>(d); 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath j = s1; s1 = s2; s2 = j; // Switch rows 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* Multiply through by the correct factors */ 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* (Eq. [2.9]) */ 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = p; 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k) 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (DenseIndex j=p; j>=0; --j) N_(k,j) *= r; 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r *= p-k; 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline >::BasisDerivativeType der; 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath basisFunctionDerivativesImpl(*this, u, order, der); 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return der; 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename _Scalar, int _Dim, int _Degree> 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <int DerivativeOrder> 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der; 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath basisFunctionDerivativesImpl(*this, u, order, der); 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return der; 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_SPLINE_H 475