Quelle Spline.h Sprache: C

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 20010-2011 Hauke Heibel <hauke.heibel@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SPLINE_H
#define EIGEN_SPLINE_H

#include "SplineFwd.h"

namespace Eigen
{
    /**
     * \ingroup Splines_Module
     * \class Spline
     * \brief A class representing multi-dimensional spline curves.
     *
     * The class represents B-splines with non-uniform knot vectors. Each control
     * point of the B-spline is associated with a basis function
     * \f{align*}
     *   C(u) & = \sum_{i=0}^{n}N_{i,p}(u)P_i
     * \f}
     *
     * \tparam _Scalar The underlying data type (typically float or double)
     * \tparam _Dim The curve dimension (e.g. 2 or 3)
     * \tparam _Degree Per default set to Dynamic; could be set to the actual desired
     *                degree for optimization purposes (would result in stack allocation
     *                of several temporary variables).
     **/
  template <typename _Scalar, int _Dim, int _Degree>
  class Spline
  {
  public:
    typedef _Scalar Scalar; /*!< The spline curve's scalar type. */
    enum { Dimension = _Dim /*!< The spline curve's dimension. */ };
    enum { Degree = _Degree /*!< The spline curve's degree. */ };

    /** \brief The point type the spline is representing. */
    typedef typename SplineTraits<Spline>::PointType PointType;

    /** \brief The data type used to store knot vectors. */
    typedef typename SplineTraits<Spline>::KnotVectorType KnotVectorType;

    /** \brief The data type used to store parameter vectors. */
    typedef typename SplineTraits<Spline>::ParameterVectorType ParameterVectorType;

    /** \brief The data type used to store non-zero basis functions. */
    typedef typename SplineTraits<Spline>::BasisVectorType BasisVectorType;

    /** \brief The data type used to store the values of the basis function derivatives. */
    typedef typename SplineTraits<Spline>::BasisDerivativeType BasisDerivativeType;

    /** \brief The data type representing the spline's control points. */
    typedef typename SplineTraits<Spline>::ControlPointVectorType ControlPointVectorType;

    /**
    * \brief Creates a (constant) zero spline.
    * For Splines with dynamic degree, the resulting degree will be 0.
    **/
    Spline()
    : m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2))
    , m_ctrls(ControlPointVectorType::Zero(Dimension,(Degree==Dynamic ? 1 : Degree+1)))
    {
      // in theory this code can go to the initializer list but it will get pretty
      // much unreadable ...
      enum { MinDegree = (Degree==Dynamic ? 0 : Degree) };
      m_knots.template segment<MinDegree+1>(0) = Array<Scalar,1,MinDegree+1>::Zero();
      m_knots.template segment<MinDegree+1>(MinDegree+1) = Array<Scalar,1,MinDegree+1>::Ones();
    }

    /**
    * \brief Creates a spline from a knot vector and control points.
    * \param knots The spline's knot vector.
    * \param ctrls The spline's control point vector.
    **/
    template <typename OtherVectorType, typename OtherArrayType>
    Spline(const OtherVectorType& knots, const OtherArrayType& ctrls) : m_knots(knots), m_ctrls(ctrls) {}

    /**
    * \brief Copy constructor for splines.
    * \param spline The input spline.
    **/
    template <int OtherDegree>
    Spline(const Spline<Scalar, Dimension, OtherDegree>& spline) :
    m_knots(spline.knots()), m_ctrls(spline.ctrls()) {}

    /**
     * \brief Returns the knots of the underlying spline.
     **/
    const KnotVectorType& knots() const { return m_knots; }

    /**
     * \brief Returns the ctrls of the underlying spline.
     **/
    const ControlPointVectorType& ctrls() const { return m_ctrls; }

    /**
     * \brief Returns the spline value at a given site \f$u\f$.
     *
     * The function returns
     * \f{align*}
     *   C(u) & = \sum_{i=0}^{n}N_{i,p}P_i
     * \f}
     *
     * \param u Parameter \f$u \in [0;1]\f$ at which the spline is evaluated.
     * \return The spline value at the given location \f$u\f$.
     **/
    PointType operator()(Scalar u) const;

    /**
     * \brief Evaluation of spline derivatives of up-to given order.
     *
     * The function returns
     * \f{align*}
     *   \frac{d^i}{du^i}C(u) & = \sum_{i=0}^{n} \frac{d^i}{du^i} N_{i,p}(u)P_i
     * \f}
     * for i ranging between 0 and order.
     *
     * \param u Parameter \f$u \in [0;1]\f$ at which the spline derivative is evaluated.
     * \param order The order up to which the derivatives are computed.
     **/
    typename SplineTraits<Spline>::DerivativeType
      derivatives(Scalar u, DenseIndex order) const;

    /**
     * \copydoc Spline::derivatives
     * Using the template version of this function is more efficieent since
     * temporary objects are allocated on the stack whenever this is possible.
     **/
    template <int DerivativeOrder>
    typename SplineTraits<Spline,DerivativeOrder>::DerivativeType
      derivatives(Scalar u, DenseIndex order = DerivativeOrder) const;

    /**
     * \brief Computes the non-zero basis functions at the given site.
     *
     * Splines have local support and a point from their image is defined
     * by exactly \f$p+1\f$ control points \f$P_i\f$ where \f$p\f$ is the
     * spline degree.
     *
     * This function computes the \f$p+1\f$ non-zero basis function values
     * for a given parameter value \f$u\f$. It returns
     * \f{align*}{
     *   N_{i,p}(u), \hdots, N_{i+p+1,p}(u)
     * \f}
     *
     * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis functions
     *          are computed.
     **/
    typename SplineTraits<Spline>::BasisVectorType
      basisFunctions(Scalar u) const;

    /**
     * \brief Computes the non-zero spline basis function derivatives up to given order.
     *
     * The function computes
     * \f{align*}{
     *   \frac{d^i}{du^i} N_{i,p}(u), \hdots, \frac{d^i}{du^i} N_{i+p+1,p}(u)
     * \f}
     * with i ranging from 0 up to the specified order.
     *
     * \param u Parameter \f$u \in [0;1]\f$ at which the non-zero basis function
     *          derivatives are computed.
     * \param order The order up to which the basis function derivatives are computes.
     **/
    typename SplineTraits<Spline>::BasisDerivativeType
      basisFunctionDerivatives(Scalar u, DenseIndex order) const;

    /**
     * \copydoc Spline::basisFunctionDerivatives
     * Using the template version of this function is more efficieent since
     * temporary objects are allocated on the stack whenever this is possible.
     **/
    template <int DerivativeOrder>
    typename SplineTraits<Spline,DerivativeOrder>::BasisDerivativeType
      basisFunctionDerivatives(Scalar u, DenseIndex order = DerivativeOrder) const;

    /**
     * \brief Returns the spline degree.
     **/
    DenseIndex degree() const;

    /**
     * \brief Returns the span within the knot vector in which u is falling.
     * \param u The site for which the span is determined.
     **/
    DenseIndex span(Scalar u) const;

    /**
     * \brief Computes the span within the provided knot vector in which u is falling.
     **/
    static DenseIndex Span(typename SplineTraits<Spline>::Scalar u, DenseIndex degree, const typename SplineTraits<Spline>::KnotVectorType& knots);

    /**
     * \brief Returns the spline's non-zero basis functions.
     *
     * The function computes and returns
     * \f{align*}{
     *   N_{i,p}(u), \hdots, N_{i+p+1,p}(u)
     * \f}
     *
     * \param u The site at which the basis functions are computed.
     * \param degree The degree of the underlying spline.
     * \param knots The underlying spline's knot vector.
     **/
    static BasisVectorType BasisFunctions(Scalar u, DenseIndex degree, const KnotVectorType& knots);

    /**
     * \copydoc Spline::basisFunctionDerivatives
     * \param degree The degree of the underlying spline
     * \param knots The underlying spline's knot vector.
     **/
    static BasisDerivativeType BasisFunctionDerivatives(
      const Scalar u, const DenseIndex order, const DenseIndex degree, const KnotVectorType& knots);

  private:
    KnotVectorType m_knots; /*!< Knot vector. */
    ControlPointVectorType  m_ctrls; /*!< Control points. */

    template <typename DerivativeType>
    static void BasisFunctionDerivativesImpl(
      const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
      const DenseIndex order,
      const DenseIndex p,
      const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U,
      DerivativeType& N_);
  };

  template <typename _Scalar, int _Dim, int _Degree>
  DenseIndex Spline<_Scalar, _Dim, _Degree>::Span(
    typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::Scalar u,
    DenseIndex degree,
    const typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::KnotVectorType& knots)
  {
    // Piegl & Tiller, "The NURBS Book", A2.1 (p. 68)
    if (u <= knots(0)) return degree;
    const Scalar* pos = std::upper_bound(knots.data()+degree-1, knots.data()+knots.size()-degree-1, u);
    return static_cast<DenseIndex>( std::distance(knots.data(), pos) - 1 );
  }

  template <typename _Scalar, int _Dim, int _Degree>
  typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType
    Spline<_Scalar, _Dim, _Degree>::BasisFunctions(
    typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
    DenseIndex degree,
    const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
  {
    const DenseIndex p = degree;
    const DenseIndex i = Spline::Span(u, degree, knots);

    const KnotVectorType& U = knots;

    BasisVectorType left(p+1); left(0) = Scalar(0);
    BasisVectorType right(p+1); right(0) = Scalar(0);

    VectorBlock<BasisVectorType,Degree>(left,1,p) = u - VectorBlock<const KnotVectorType,Degree>(U,i+1-p,p).reverse();
    VectorBlock<BasisVectorType,Degree>(right,1,p) = VectorBlock<const KnotVectorType,Degree>(U,i+1,p) - u;

    BasisVectorType N(1,p+1);
    N(0) = Scalar(1);
    for (DenseIndex j=1; j<=p; ++j)
    {
      Scalar saved = Scalar(0);
      for (DenseIndex r=0; r<j; r++)
      {
        const Scalar tmp = N(r)/(right(r+1)+left(j-r));
        N[r] = saved + right(r+1)*tmp;
        saved = left(j-r)*tmp;
      }
      N(j) = saved;
    }
    return N;
  }

  template <typename _Scalar, int _Dim, int _Degree>
  DenseIndex Spline<_Scalar, _Dim, _Degree>::degree() const
  {
    if (_Degree == Dynamic)
      return m_knots.size() - m_ctrls.cols() - 1;
    else
      return _Degree;
  }

  template <typename _Scalar, int _Dim, int _Degree>
  DenseIndex Spline<_Scalar, _Dim, _Degree>::span(Scalar u) const
  {
    return Spline::Span(u, degree(), knots());
  }

  template <typename _Scalar, int _Dim, int _Degree>
  typename Spline<_Scalar, _Dim, _Degree>::PointType Spline<_Scalar, _Dim, _Degree>::operator()(Scalar u) const
  {
    enum { Order = SplineTraits<Spline>::OrderAtCompileTime };

    const DenseIndex span = this->span(u);
    const DenseIndex p = degree();
    const BasisVectorType basis_funcs = basisFunctions(u);

    const Replicate<BasisVectorType,Dimension,1> ctrl_weights(basis_funcs);
    const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(ctrls(),0,span-p,Dimension,p+1);
    return (ctrl_weights * ctrl_pts).rowwise().sum();
  }

  /* --------------------------------------------------------------------------------------------- */

  template <typename SplineType, typename DerivativeType>
  void derivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& der)
  {
    enum { Dimension = SplineTraits<SplineType>::Dimension };
    enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };
    enum { DerivativeOrder = DerivativeType::ColsAtCompileTime };

    typedef typename SplineTraits<SplineType>::ControlPointVectorType ControlPointVectorType;
    typedef typename SplineTraits<SplineType,DerivativeOrder>::BasisDerivativeType BasisDerivativeType;
    typedef typename BasisDerivativeType::ConstRowXpr BasisDerivativeRowXpr;

    const DenseIndex p = spline.degree();
    const DenseIndex span = spline.span(u);

    const DenseIndex n = (std::min)(p, order);

    der.resize(Dimension,n+1);

    // Retrieve the basis function derivatives up to the desired order...
    const BasisDerivativeType basis_func_ders = spline.template basisFunctionDerivatives<DerivativeOrder>(u, n+1);

    // ... and perform the linear combinations of the control points.
    for (DenseIndex der_order=0; der_order<n+1; ++der_order)
    {
      const Replicate<BasisDerivativeRowXpr,Dimension,1> ctrl_weights( basis_func_ders.row(der_order) );
      const Block<const ControlPointVectorType,Dimension,Order> ctrl_pts(spline.ctrls(),0,span-p,Dimension,p+1);
      der.col(der_order) = (ctrl_weights * ctrl_pts).rowwise().sum();
    }
  }

  template <typename _Scalar, int _Dim, int _Degree>
  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::DerivativeType
    Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
  {
    typename SplineTraits< Spline >::DerivativeType res;
    derivativesImpl(*this, u, order, res);
    return res;
  }

  template <typename _Scalar, int _Dim, int _Degree>
  template <int DerivativeOrder>
  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::DerivativeType
    Spline<_Scalar, _Dim, _Degree>::derivatives(Scalar u, DenseIndex order) const
  {
    typename SplineTraits< Spline, DerivativeOrder >::DerivativeType res;
    derivativesImpl(*this, u, order, res);
    return res;
  }

  template <typename _Scalar, int _Dim, int _Degree>
  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisVectorType
    Spline<_Scalar, _Dim, _Degree>::basisFunctions(Scalar u) const
  {
    return Spline::BasisFunctions(u, degree(), knots());
  }

  /* --------------------------------------------------------------------------------------------- */


  template <typename _Scalar, int _Dim, int _Degree>
  template <typename DerivativeType>
  void Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivativesImpl(
    const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
    const DenseIndex order,
    const DenseIndex p,
    const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& U,
    DerivativeType& N_)
  {
    typedef Spline<_Scalar, _Dim, _Degree> SplineType;
    enum { Order = SplineTraits<SplineType>::OrderAtCompileTime };

    const DenseIndex span = SplineType::Span(u, p, U);

    const DenseIndex n = (std::min)(p, order);

    N_.resize(n+1, p+1);

    BasisVectorType left = BasisVectorType::Zero(p+1);
    BasisVectorType right = BasisVectorType::Zero(p+1);

    Matrix<Scalar,Order,Order> ndu(p+1,p+1);

    Scalar saved, temp; // FIXME These were double instead of Scalar. Was there a reason for that?

    ndu(0,0) = 1.0;

    DenseIndex j;
    for (j=1; j<=p; ++j)
    {
      left[j] = u-U[span+1-j];
      right[j] = U[span+j]-u;
      saved = 0.0;

      for (DenseIndex r=0; r<j; ++r)
      {
        /* Lower triangle */
        ndu(j,r) = right[r+1]+left[j-r];
        temp = ndu(r,j-1)/ndu(j,r);
        /* Upper triangle */
        ndu(r,j) = static_cast<Scalar>(saved+right[r+1] * temp);
        saved = left[j-r] * temp;
      }

      ndu(j,j) = static_cast<Scalar>(saved);
    }

    for (j = p; j>=0; --j)
      N_(0,j) = ndu(j,p);

    // Compute the derivatives
    DerivativeType a(n+1,p+1);
    DenseIndex r=0;
    for (; r<=p; ++r)
    {
      DenseIndex s1,s2;
      s1 = 0; s2 = 1; // alternate rows in array a
      a(0,0) = 1.0;

      // Compute the k-th derivative
      for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
      {
        Scalar d = 0.0;
        DenseIndex rk,pk,j1,j2;
        rk = r-k; pk = p-k;

        if (r>=k)
        {
          a(s2,0) = a(s1,0)/ndu(pk+1,rk);
          d = a(s2,0)*ndu(rk,pk);
        }

        if (rk>=-1) j1 = 1;
        else        j1 = -rk;

        if (r-1 <= pk) j2 = k-1;
        else           j2 = p-r;

        for (j=j1; j<=j2; ++j)
        {
          a(s2,j) = (a(s1,j)-a(s1,j-1))/ndu(pk+1,rk+j);
          d += a(s2,j)*ndu(rk+j,pk);
        }

        if (r<=pk)
        {
          a(s2,k) = -a(s1,k-1)/ndu(pk+1,r);
          d += a(s2,k)*ndu(r,pk);
        }

        N_(k,r) = static_cast<Scalar>(d);
        j = s1; s1 = s2; s2 = j; // Switch rows
      }
    }

    /* Multiply through by the correct factors */
    /* (Eq. [2.9])                             */
    r = p;
    for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k)
    {
      for (j=p; j>=0; --j) N_(k,j) *= r;
      r *= p-k;
    }
  }

  template <typename _Scalar, int _Dim, int _Degree>
  typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
    Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
  {
    typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType der;
    BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
    return der;
  }

  template <typename _Scalar, int _Dim, int _Degree>
  template <int DerivativeOrder>
  typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType
    Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const
  {
    typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType der;
    BasisFunctionDerivativesImpl(u, order, degree(), knots(), der);
    return der;
  }

  template <typename _Scalar, int _Dim, int _Degree>
  typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType
  Spline<_Scalar, _Dim, _Degree>::BasisFunctionDerivatives(
    const typename Spline<_Scalar, _Dim, _Degree>::Scalar u,
    const DenseIndex order,
    const DenseIndex degree,
    const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
  {
    typename SplineTraits<Spline>::BasisDerivativeType der;
    BasisFunctionDerivativesImpl(u, order, degree, knots, der);
    return der;
  }
}

#endif // EIGEN_SPLINE_H

quality92%

¤ Dauer der Verarbeitung: 0.17 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.