diff options
Diffstat (limited to 'eigen/unsupported/Eigen/src/Splines')
| -rw-r--r-- | eigen/unsupported/Eigen/src/Splines/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | eigen/unsupported/Eigen/src/Splines/Spline.h | 74 | ||||
| -rw-r--r-- | eigen/unsupported/Eigen/src/Splines/SplineFitting.h | 274 | ||||
| -rw-r--r-- | eigen/unsupported/Eigen/src/Splines/SplineFwd.h | 3 |
4 files changed, 333 insertions, 24 deletions
diff --git a/eigen/unsupported/Eigen/src/Splines/CMakeLists.txt b/eigen/unsupported/Eigen/src/Splines/CMakeLists.txt deleted file mode 100644 index 55c6271..0000000 --- a/eigen/unsupported/Eigen/src/Splines/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Splines_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Splines_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Splines COMPONENT Devel - ) diff --git a/eigen/unsupported/Eigen/src/Splines/Spline.h b/eigen/unsupported/Eigen/src/Splines/Spline.h index 771f104..627f6e4 100644 --- a/eigen/unsupported/Eigen/src/Splines/Spline.h +++ b/eigen/unsupported/Eigen/src/Splines/Spline.h @@ -44,9 +44,15 @@ namespace Eigen /** \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; @@ -57,7 +63,7 @@ namespace Eigen **/ Spline() : m_knots(1, (Degree==Dynamic ? 2 : 2*Degree+2)) - , m_ctrls(ControlPointVectorType::Zero(2,(Degree==Dynamic ? 1 : Degree+1))) + , 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 ... @@ -88,7 +94,7 @@ namespace Eigen const KnotVectorType& knots() const { return m_knots; } /** - * \brief Returns the knots of the underlying spline. + * \brief Returns the ctrls of the underlying spline. **/ const ControlPointVectorType& ctrls() const { return m_ctrls; } @@ -203,10 +209,25 @@ namespace Eigen **/ 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> @@ -345,20 +366,24 @@ namespace Eigen } /* --------------------------------------------------------------------------------------------- */ - - template <typename SplineType, typename DerivativeType> - void basisFunctionDerivativesImpl(const SplineType& spline, typename SplineType::Scalar u, DenseIndex order, DerivativeType& N_) + + + 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 }; typedef typename SplineTraits<SplineType>::Scalar Scalar; typedef typename SplineTraits<SplineType>::BasisVectorType BasisVectorType; - typedef typename SplineTraits<SplineType>::KnotVectorType KnotVectorType; - - const KnotVectorType& U = spline.knots(); - - const DenseIndex p = spline.degree(); - const DenseIndex span = spline.span(u); + + const DenseIndex span = SplineType::Span(u, p, U); const DenseIndex n = (std::min)(p, order); @@ -369,7 +394,7 @@ namespace Eigen Matrix<Scalar,Order,Order> ndu(p+1,p+1); - double saved, temp; + Scalar saved, temp; // FIXME These were double instead of Scalar. Was there a reason for that? ndu(0,0) = 1.0; @@ -408,7 +433,7 @@ namespace Eigen // Compute the k-th derivative for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k) { - double d = 0.0; + Scalar d = 0.0; DenseIndex rk,pk,j1,j2; rk = r-k; pk = p-k; @@ -446,7 +471,7 @@ namespace Eigen r = p; for (DenseIndex k=1; k<=static_cast<DenseIndex>(n); ++k) { - for (DenseIndex j=p; j>=0; --j) N_(k,j) *= r; + for (j=p; j>=0; --j) N_(k,j) *= r; r *= p-k; } } @@ -455,8 +480,8 @@ namespace Eigen typename SplineTraits< Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const { - typename SplineTraits< Spline >::BasisDerivativeType der; - basisFunctionDerivativesImpl(*this, u, order, der); + typename SplineTraits<Spline<_Scalar, _Dim, _Degree> >::BasisDerivativeType der; + BasisFunctionDerivativesImpl(u, order, degree(), knots(), der); return der; } @@ -465,8 +490,21 @@ namespace Eigen typename SplineTraits< Spline<_Scalar, _Dim, _Degree>, DerivativeOrder >::BasisDerivativeType Spline<_Scalar, _Dim, _Degree>::basisFunctionDerivatives(Scalar u, DenseIndex order) const { - typename SplineTraits< Spline, DerivativeOrder >::BasisDerivativeType der; - basisFunctionDerivativesImpl(*this, u, order, der); + 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; } } diff --git a/eigen/unsupported/Eigen/src/Splines/SplineFitting.h b/eigen/unsupported/Eigen/src/Splines/SplineFitting.h index 0265d53..c761a9b 100644 --- a/eigen/unsupported/Eigen/src/Splines/SplineFitting.h +++ b/eigen/unsupported/Eigen/src/Splines/SplineFitting.h @@ -10,10 +10,14 @@ #ifndef EIGEN_SPLINE_FITTING_H #define EIGEN_SPLINE_FITTING_H +#include <algorithm> +#include <functional> #include <numeric> +#include <vector> #include "SplineFwd.h" +#include <Eigen/LU> #include <Eigen/QR> namespace Eigen @@ -50,6 +54,129 @@ namespace Eigen } /** + * \brief Computes knot averages when derivative constraints are present. + * Note that this is a technical interpretation of the referenced article + * since the algorithm contained therein is incorrect as written. + * \ingroup Splines_Module + * + * \param[in] parameters The parameters at which the interpolation B-Spline + * will intersect the given interpolation points. The parameters + * are assumed to be a non-decreasing sequence. + * \param[in] degree The degree of the interpolating B-Spline. This must be + * greater than zero. + * \param[in] derivativeIndices The indices corresponding to parameters at + * which there are derivative constraints. The indices are assumed + * to be a non-decreasing sequence. + * \param[out] knots The calculated knot vector. These will be returned as a + * non-decreasing sequence + * + * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008. + * Curve interpolation with directional constraints for engineering design. + * Engineering with Computers + **/ + template <typename KnotVectorType, typename ParameterVectorType, typename IndexArray> + void KnotAveragingWithDerivatives(const ParameterVectorType& parameters, + const unsigned int degree, + const IndexArray& derivativeIndices, + KnotVectorType& knots) + { + typedef typename ParameterVectorType::Scalar Scalar; + + DenseIndex numParameters = parameters.size(); + DenseIndex numDerivatives = derivativeIndices.size(); + + if (numDerivatives < 1) + { + KnotAveraging(parameters, degree, knots); + return; + } + + DenseIndex startIndex; + DenseIndex endIndex; + + DenseIndex numInternalDerivatives = numDerivatives; + + if (derivativeIndices[0] == 0) + { + startIndex = 0; + --numInternalDerivatives; + } + else + { + startIndex = 1; + } + if (derivativeIndices[numDerivatives - 1] == numParameters - 1) + { + endIndex = numParameters - degree; + --numInternalDerivatives; + } + else + { + endIndex = numParameters - degree - 1; + } + + // There are (endIndex - startIndex + 1) knots obtained from the averaging + // and 2 for the first and last parameters. + DenseIndex numAverageKnots = endIndex - startIndex + 3; + KnotVectorType averageKnots(numAverageKnots); + averageKnots[0] = parameters[0]; + + int newKnotIndex = 0; + for (DenseIndex i = startIndex; i <= endIndex; ++i) + averageKnots[++newKnotIndex] = parameters.segment(i, degree).mean(); + averageKnots[++newKnotIndex] = parameters[numParameters - 1]; + + newKnotIndex = -1; + + ParameterVectorType temporaryParameters(numParameters + 1); + KnotVectorType derivativeKnots(numInternalDerivatives); + for (DenseIndex i = 0; i < numAverageKnots - 1; ++i) + { + temporaryParameters[0] = averageKnots[i]; + ParameterVectorType parameterIndices(numParameters); + int temporaryParameterIndex = 1; + for (DenseIndex j = 0; j < numParameters; ++j) + { + Scalar parameter = parameters[j]; + if (parameter >= averageKnots[i] && parameter < averageKnots[i + 1]) + { + parameterIndices[temporaryParameterIndex] = j; + temporaryParameters[temporaryParameterIndex++] = parameter; + } + } + temporaryParameters[temporaryParameterIndex] = averageKnots[i + 1]; + + for (int j = 0; j <= temporaryParameterIndex - 2; ++j) + { + for (DenseIndex k = 0; k < derivativeIndices.size(); ++k) + { + if (parameterIndices[j + 1] == derivativeIndices[k] + && parameterIndices[j + 1] != 0 + && parameterIndices[j + 1] != numParameters - 1) + { + derivativeKnots[++newKnotIndex] = temporaryParameters.segment(j, 3).mean(); + break; + } + } + } + } + + KnotVectorType temporaryKnots(averageKnots.size() + derivativeKnots.size()); + + std::merge(averageKnots.data(), averageKnots.data() + averageKnots.size(), + derivativeKnots.data(), derivativeKnots.data() + derivativeKnots.size(), + temporaryKnots.data()); + + // Number of knots (one for each point and derivative) plus spline order. + DenseIndex numKnots = numParameters + numDerivatives + degree + 1; + knots.resize(numKnots); + + knots.head(degree).fill(temporaryKnots[0]); + knots.tail(degree).fill(temporaryKnots.template tail<1>()[0]); + knots.segment(degree, temporaryKnots.size()) = temporaryKnots; + } + + /** * \brief Computes chord length parameters which are required for spline interpolation. * \ingroup Splines_Module * @@ -86,6 +213,7 @@ namespace Eigen struct SplineFitting { typedef typename SplineType::KnotVectorType KnotVectorType; + typedef typename SplineType::ParameterVectorType ParameterVectorType; /** * \brief Fits an interpolating Spline to the given data points. @@ -109,6 +237,52 @@ namespace Eigen **/ template <typename PointArrayType> static SplineType Interpolate(const PointArrayType& pts, DenseIndex degree, const KnotVectorType& knot_parameters); + + /** + * \brief Fits an interpolating spline to the given data points and + * derivatives. + * + * \param points The points for which an interpolating spline will be computed. + * \param derivatives The desired derivatives of the interpolating spline at interpolation + * points. + * \param derivativeIndices An array indicating which point each derivative belongs to. This + * must be the same size as @a derivatives. + * \param degree The degree of the interpolating spline. + * + * \returns A spline interpolating @a points with @a derivatives at those points. + * + * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008. + * Curve interpolation with directional constraints for engineering design. + * Engineering with Computers + **/ + template <typename PointArrayType, typename IndexArray> + static SplineType InterpolateWithDerivatives(const PointArrayType& points, + const PointArrayType& derivatives, + const IndexArray& derivativeIndices, + const unsigned int degree); + + /** + * \brief Fits an interpolating spline to the given data points and derivatives. + * + * \param points The points for which an interpolating spline will be computed. + * \param derivatives The desired derivatives of the interpolating spline at interpolation points. + * \param derivativeIndices An array indicating which point each derivative belongs to. This + * must be the same size as @a derivatives. + * \param degree The degree of the interpolating spline. + * \param parameters The parameters corresponding to the interpolation points. + * + * \returns A spline interpolating @a points with @a derivatives at those points. + * + * \sa Les A. Piegl, Khairan Rajab, Volha Smarodzinana. 2008. + * Curve interpolation with directional constraints for engineering design. + * Engineering with Computers + */ + template <typename PointArrayType, typename IndexArray> + static SplineType InterpolateWithDerivatives(const PointArrayType& points, + const PointArrayType& derivatives, + const IndexArray& derivativeIndices, + const unsigned int degree, + const ParameterVectorType& parameters); }; template <typename SplineType> @@ -151,6 +325,106 @@ namespace Eigen ChordLengths(pts, chord_lengths); return Interpolate(pts, degree, chord_lengths); } + + template <typename SplineType> + template <typename PointArrayType, typename IndexArray> + SplineType + SplineFitting<SplineType>::InterpolateWithDerivatives(const PointArrayType& points, + const PointArrayType& derivatives, + const IndexArray& derivativeIndices, + const unsigned int degree, + const ParameterVectorType& parameters) + { + typedef typename SplineType::KnotVectorType::Scalar Scalar; + typedef typename SplineType::ControlPointVectorType ControlPointVectorType; + + typedef Matrix<Scalar, Dynamic, Dynamic> MatrixType; + + const DenseIndex n = points.cols() + derivatives.cols(); + + KnotVectorType knots; + + KnotAveragingWithDerivatives(parameters, degree, derivativeIndices, knots); + + // fill matrix + MatrixType A = MatrixType::Zero(n, n); + + // Use these dimensions for quicker populating, then transpose for solving. + MatrixType b(points.rows(), n); + + DenseIndex startRow; + DenseIndex derivativeStart; + + // End derivatives. + if (derivativeIndices[0] == 0) + { + A.template block<1, 2>(1, 0) << -1, 1; + + Scalar y = (knots(degree + 1) - knots(0)) / degree; + b.col(1) = y*derivatives.col(0); + + startRow = 2; + derivativeStart = 1; + } + else + { + startRow = 1; + derivativeStart = 0; + } + if (derivativeIndices[derivatives.cols() - 1] == points.cols() - 1) + { + A.template block<1, 2>(n - 2, n - 2) << -1, 1; + + Scalar y = (knots(knots.size() - 1) - knots(knots.size() - (degree + 2))) / degree; + b.col(b.cols() - 2) = y*derivatives.col(derivatives.cols() - 1); + } + + DenseIndex row = startRow; + DenseIndex derivativeIndex = derivativeStart; + for (DenseIndex i = 1; i < parameters.size() - 1; ++i) + { + const DenseIndex span = SplineType::Span(parameters[i], degree, knots); + + if (derivativeIndices[derivativeIndex] == i) + { + A.block(row, span - degree, 2, degree + 1) + = SplineType::BasisFunctionDerivatives(parameters[i], 1, degree, knots); + + b.col(row++) = points.col(i); + b.col(row++) = derivatives.col(derivativeIndex++); + } + else + { + A.row(row++).segment(span - degree, degree + 1) + = SplineType::BasisFunctions(parameters[i], degree, knots); + } + } + b.col(0) = points.col(0); + b.col(b.cols() - 1) = points.col(points.cols() - 1); + A(0,0) = 1; + A(n - 1, n - 1) = 1; + + // Solve + FullPivLU<MatrixType> lu(A); + ControlPointVectorType controlPoints = lu.solve(MatrixType(b.transpose())).transpose(); + + SplineType spline(knots, controlPoints); + + return spline; + } + + template <typename SplineType> + template <typename PointArrayType, typename IndexArray> + SplineType + SplineFitting<SplineType>::InterpolateWithDerivatives(const PointArrayType& points, + const PointArrayType& derivatives, + const IndexArray& derivativeIndices, + const unsigned int degree) + { + ParameterVectorType parameters; + ChordLengths(points, parameters); + return InterpolateWithDerivatives(points, derivatives, derivativeIndices, degree, parameters); + } } #endif // EIGEN_SPLINE_FITTING_H diff --git a/eigen/unsupported/Eigen/src/Splines/SplineFwd.h b/eigen/unsupported/Eigen/src/Splines/SplineFwd.h index 9ea23a9..0a95fbf 100644 --- a/eigen/unsupported/Eigen/src/Splines/SplineFwd.h +++ b/eigen/unsupported/Eigen/src/Splines/SplineFwd.h @@ -48,6 +48,9 @@ namespace Eigen /** \brief The data type used to store knot vectors. */ typedef Array<Scalar,1,Dynamic> KnotVectorType; + + /** \brief The data type used to store parameter vectors. */ + typedef Array<Scalar,1,Dynamic> ParameterVectorType; /** \brief The data type representing the spline's control points. */ typedef Array<Scalar,Dimension,Dynamic> ControlPointVectorType; |