Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/unsupported/Eigen/src/Splines/CMakeLists.txt b/unsupported/Eigen/src/Splines/CMakeLists.txt
deleted file mode 100644
index 55c6271..0000000
--- a/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/unsupported/Eigen/src/Splines/Spline.h b/unsupported/Eigen/src/Splines/Spline.h
index 771f104..57788c8 100644
--- a/unsupported/Eigen/src/Splines/Spline.h
+++ b/unsupported/Eigen/src/Splines/Spline.h
@@ -44,9 +44,15 @@
/** \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 @@
**/
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 @@
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 @@
**/
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>
@@ -228,8 +249,6 @@
DenseIndex degree,
const typename Spline<_Scalar, _Dim, _Degree>::KnotVectorType& knots)
{
- typedef typename Spline<_Scalar, _Dim, _Degree>::BasisVectorType BasisVectorType;
-
const DenseIndex p = degree;
const DenseIndex i = Spline::Span(u, degree, knots);
@@ -345,20 +364,21 @@
}
/* --------------------------------------------------------------------------------------------- */
-
- 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 +389,7 @@
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 +428,7 @@
// 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 +466,7 @@
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 +475,8 @@
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 +485,21 @@
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/unsupported/Eigen/src/Splines/SplineFitting.h b/unsupported/Eigen/src/Splines/SplineFitting.h
index 0265d53..c761a9b 100644
--- a/unsupported/Eigen/src/Splines/SplineFitting.h
+++ b/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 @@
}
/**
+ * \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 @@
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 @@
**/
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 @@
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/unsupported/Eigen/src/Splines/SplineFwd.h b/unsupported/Eigen/src/Splines/SplineFwd.h
index 49db8d3..0a95fbf 100644
--- a/unsupported/Eigen/src/Splines/SplineFwd.h
+++ b/unsupported/Eigen/src/Splines/SplineFwd.h
@@ -31,6 +31,8 @@
enum { OrderAtCompileTime = _Degree==Dynamic ? Dynamic : _Degree+1 /*!< The spline curve's order at compile-time. */ };
enum { NumOfDerivativesAtCompileTime = OrderAtCompileTime /*!< The number of derivatives defined for the current spline. */ };
+
+ enum { DerivativeMemoryLayout = Dimension==1 ? RowMajor : ColMajor /*!< The derivative type's memory layout. */ };
/** \brief The data type used to store non-zero basis functions. */
typedef Array<Scalar,1,OrderAtCompileTime> BasisVectorType;
@@ -39,13 +41,16 @@
typedef Array<Scalar,Dynamic,Dynamic,RowMajor,NumOfDerivativesAtCompileTime,OrderAtCompileTime> BasisDerivativeType;
/** \brief The data type used to store the spline's derivative values. */
- typedef Array<Scalar,Dimension,Dynamic,ColMajor,Dimension,NumOfDerivativesAtCompileTime> DerivativeType;
+ typedef Array<Scalar,Dimension,Dynamic,DerivativeMemoryLayout,Dimension,NumOfDerivativesAtCompileTime> DerivativeType;
/** \brief The point type the spline is representing. */
typedef Array<Scalar,Dimension,1> PointType;
/** \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;
@@ -62,12 +67,14 @@
{
enum { OrderAtCompileTime = _Degree==Dynamic ? Dynamic : _Degree+1 /*!< The spline curve's order at compile-time. */ };
enum { NumOfDerivativesAtCompileTime = _DerivativeOrder==Dynamic ? Dynamic : _DerivativeOrder+1 /*!< The number of derivatives defined for the current spline. */ };
+
+ enum { DerivativeMemoryLayout = _Dim==1 ? RowMajor : ColMajor /*!< The derivative type's memory layout. */ };
/** \brief The data type used to store the values of the basis function derivatives. */
typedef Array<_Scalar,Dynamic,Dynamic,RowMajor,NumOfDerivativesAtCompileTime,OrderAtCompileTime> BasisDerivativeType;
/** \brief The data type used to store the spline's derivative values. */
- typedef Array<_Scalar,_Dim,Dynamic,ColMajor,_Dim,NumOfDerivativesAtCompileTime> DerivativeType;
+ typedef Array<_Scalar,_Dim,Dynamic,DerivativeMemoryLayout,_Dim,NumOfDerivativesAtCompileTime> DerivativeType;
};
/** \brief 2D float B-spline with dynamic degree. */