mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-04-16 14:49:39 +08:00
merge with default branch
This commit is contained in:
Gael Guennebaud
commit
3c7686630d
@ -558,6 +558,7 @@ namespace internal {
|
||||
while(m_outerPos<m_end)
|
||||
{
|
||||
m_outerPos++;
|
||||
if(m_outerPos==m_end) break;
|
||||
typename XprType::InnerIterator it(m_block.m_matrix, m_outerPos);
|
||||
// search for the key m_innerIndex in the current outer-vector
|
||||
while(it && it.index() < m_innerIndex) ++it;
|
||||
|
||||
@ -165,7 +165,10 @@ template<typename Lhs, typename Rhs, int InnerSize> struct SparseDenseProductRet
|
||||
|
||||
template<typename Lhs, typename Rhs> struct SparseDenseProductReturnType<Lhs,Rhs,1>
|
||||
{
|
||||
typedef SparseDenseOuterProduct<Lhs,Rhs,false> Type;
|
||||
typedef typename internal::conditional<
|
||||
Lhs::IsRowMajor,
|
||||
SparseDenseOuterProduct<Rhs,Lhs,true>,
|
||||
SparseDenseOuterProduct<Lhs,Rhs,false> >::type Type;
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductReturnType
|
||||
@ -175,7 +178,10 @@ template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductRet
|
||||
|
||||
template<typename Lhs, typename Rhs> struct DenseSparseProductReturnType<Lhs,Rhs,1>
|
||||
{
|
||||
typedef SparseDenseOuterProduct<Rhs,Lhs,true> Type;
|
||||
typedef typename internal::conditional<
|
||||
Rhs::IsRowMajor,
|
||||
SparseDenseOuterProduct<Rhs,Lhs,true>,
|
||||
SparseDenseOuterProduct<Lhs,Rhs,false> >::type Type;
|
||||
};
|
||||
|
||||
namespace internal {
|
||||
@ -260,17 +266,30 @@ class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::InnerIterator : public _LhsNes
|
||||
typedef typename SparseDenseOuterProduct::Index Index;
|
||||
public:
|
||||
EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer)
|
||||
: Base(prod.lhs(), 0), m_outer(outer), m_factor(prod.rhs().coeff(outer))
|
||||
{
|
||||
}
|
||||
: Base(prod.lhs(), 0), m_outer(outer), m_factor(get(prod.rhs(), outer, typename internal::traits<Rhs>::StorageKind() ))
|
||||
{ }
|
||||
|
||||
inline Index outer() const { return m_outer; }
|
||||
inline Index row() const { return Transpose ? Base::row() : m_outer; }
|
||||
inline Index col() const { return Transpose ? m_outer : Base::row(); }
|
||||
inline Index row() const { return Transpose ? m_outer : Base::index(); }
|
||||
inline Index col() const { return Transpose ? Base::index() : m_outer; }
|
||||
|
||||
inline Scalar value() const { return Base::value() * m_factor; }
|
||||
|
||||
protected:
|
||||
static Scalar get(const _RhsNested &rhs, Index outer, Dense = Dense())
|
||||
{
|
||||
return rhs.coeff(outer);
|
||||
}
|
||||
|
||||
static Scalar get(const _RhsNested &rhs, Index outer, Sparse = Sparse())
|
||||
{
|
||||
typename Traits::_RhsNested::InnerIterator it(rhs, outer);
|
||||
if (it && it.index()==0)
|
||||
return it.value();
|
||||
|
||||
return Scalar(0);
|
||||
}
|
||||
|
||||
Index m_outer;
|
||||
Scalar m_factor;
|
||||
};
|
||||
|
||||
@ -105,8 +105,11 @@ template<typename Lhs, typename Rhs> class DenseTimeSparseProduct;
|
||||
template<typename Lhs, typename Rhs, bool Transpose> class SparseDenseOuterProduct;
|
||||
|
||||
template<typename Lhs, typename Rhs> struct SparseSparseProductReturnType;
|
||||
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct DenseSparseProductReturnType;
|
||||
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct SparseDenseProductReturnType;
|
||||
template<typename Lhs, typename Rhs,
|
||||
int InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(internal::traits<Lhs>::ColsAtCompileTime,internal::traits<Rhs>::RowsAtCompileTime)> struct DenseSparseProductReturnType;
|
||||
|
||||
template<typename Lhs, typename Rhs,
|
||||
int InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(internal::traits<Lhs>::ColsAtCompileTime,internal::traits<Rhs>::RowsAtCompileTime)> struct SparseDenseProductReturnType;
|
||||
template<typename MatrixType,int UpLo> class SparseSymmetricPermutationProduct;
|
||||
|
||||
namespace internal {
|
||||
|
||||
@ -9,32 +9,6 @@
|
||||
|
||||
#include "sparse.h"
|
||||
|
||||
template<typename SparseMatrixType, typename DenseMatrix, bool IsRowMajor=SparseMatrixType::IsRowMajor> struct test_outer;
|
||||
|
||||
template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<SparseMatrixType,DenseMatrix,false> {
|
||||
static void run(SparseMatrixType& m2, SparseMatrixType& m4, DenseMatrix& refMat2, DenseMatrix& refMat4) {
|
||||
typedef typename SparseMatrixType::Index Index;
|
||||
Index c = internal::random<Index>(0,m2.cols()-1);
|
||||
Index c1 = internal::random<Index>(0,m2.cols()-1);
|
||||
VERIFY_IS_APPROX(m4=m2.col(c)*refMat2.col(c1).transpose(), refMat4=refMat2.col(c)*refMat2.col(c1).transpose());
|
||||
VERIFY_IS_APPROX(m4=refMat2.col(c1)*m2.col(c).transpose(), refMat4=refMat2.col(c1)*refMat2.col(c).transpose());
|
||||
}
|
||||
};
|
||||
|
||||
template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<SparseMatrixType,DenseMatrix,true> {
|
||||
static void run(SparseMatrixType& m2, SparseMatrixType& m4, DenseMatrix& refMat2, DenseMatrix& refMat4) {
|
||||
typedef typename SparseMatrixType::Index Index;
|
||||
Index r = internal::random<Index>(0,m2.rows()-1);
|
||||
Index c1 = internal::random<Index>(0,m2.cols()-1);
|
||||
VERIFY_IS_APPROX(m4=m2.row(r).transpose()*refMat2.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat2.col(c1).transpose());
|
||||
VERIFY_IS_APPROX(m4=refMat2.col(c1)*m2.row(r), refMat4=refMat2.col(c1)*refMat2.row(r));
|
||||
}
|
||||
};
|
||||
|
||||
// (m2,m4,refMat2,refMat4,dv1);
|
||||
// VERIFY_IS_APPROX(m4=m2.innerVector(c)*dv1.transpose(), refMat4=refMat2.colVector(c)*dv1.transpose());
|
||||
// VERIFY_IS_APPROX(m4=dv1*mcm.col(c).transpose(), refMat4=dv1*refMat2.col(c).transpose());
|
||||
|
||||
template<typename SparseMatrixType> void sparse_product()
|
||||
{
|
||||
typedef typename SparseMatrixType::Index Index;
|
||||
@ -119,7 +93,34 @@ template<typename SparseMatrixType> void sparse_product()
|
||||
VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());
|
||||
|
||||
// sparse * dense and dense * sparse outer product
|
||||
test_outer<SparseMatrixType,DenseMatrix>::run(m2,m4,refMat2,refMat4);
|
||||
{
|
||||
Index c = internal::random<Index>(0,depth-1);
|
||||
Index r = internal::random<Index>(0,rows-1);
|
||||
Index c1 = internal::random<Index>(0,cols-1);
|
||||
Index r1 = internal::random<Index>(0,depth-1);
|
||||
|
||||
VERIFY_IS_APPROX( m4=m2.col(c)*refMat3.col(c1).transpose(), refMat4=refMat2.col(c)*refMat3.col(c1).transpose());
|
||||
VERIFY_IS_APPROX( m4=m2.middleCols(c,1)*refMat3.col(c1).transpose(), refMat4=refMat2.col(c)*refMat3.col(c1).transpose());
|
||||
VERIFY_IS_APPROX(dm4=m2.col(c)*refMat3.col(c1).transpose(), refMat4=refMat2.col(c)*refMat3.col(c1).transpose());
|
||||
|
||||
VERIFY_IS_APPROX(m4=refMat3.col(c1)*m2.col(c).transpose(), refMat4=refMat3.col(c1)*refMat2.col(c).transpose());
|
||||
VERIFY_IS_APPROX(m4=refMat3.col(c1)*m2.middleCols(c,1).transpose(), refMat4=refMat3.col(c1)*refMat2.col(c).transpose());
|
||||
VERIFY_IS_APPROX(dm4=refMat3.col(c1)*m2.col(c).transpose(), refMat4=refMat3.col(c1)*refMat2.col(c).transpose());
|
||||
|
||||
VERIFY_IS_APPROX( m4=refMat3.row(r1).transpose()*m2.col(c).transpose(), refMat4=refMat3.row(r1).transpose()*refMat2.col(c).transpose());
|
||||
VERIFY_IS_APPROX(dm4=refMat3.row(r1).transpose()*m2.col(c).transpose(), refMat4=refMat3.row(r1).transpose()*refMat2.col(c).transpose());
|
||||
|
||||
VERIFY_IS_APPROX( m4=m2.row(r).transpose()*refMat3.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat3.col(c1).transpose());
|
||||
VERIFY_IS_APPROX( m4=m2.middleRows(r,1).transpose()*refMat3.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat3.col(c1).transpose());
|
||||
VERIFY_IS_APPROX(dm4=m2.row(r).transpose()*refMat3.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat3.col(c1).transpose());
|
||||
|
||||
VERIFY_IS_APPROX( m4=refMat3.col(c1)*m2.row(r), refMat4=refMat3.col(c1)*refMat2.row(r));
|
||||
VERIFY_IS_APPROX( m4=refMat3.col(c1)*m2.middleRows(r,1), refMat4=refMat3.col(c1)*refMat2.row(r));
|
||||
VERIFY_IS_APPROX(dm4=refMat3.col(c1)*m2.row(r), refMat4=refMat3.col(c1)*refMat2.row(r));
|
||||
|
||||
VERIFY_IS_APPROX( m4=refMat3.row(r1).transpose()*m2.row(r), refMat4=refMat3.row(r1).transpose()*refMat2.row(r));
|
||||
VERIFY_IS_APPROX(dm4=refMat3.row(r1).transpose()*m2.row(r), refMat4=refMat3.row(r1).transpose()*refMat2.row(r));
|
||||
}
|
||||
|
||||
VERIFY_IS_APPROX(m6=m6*m6, refMat6=refMat6*refMat6);
|
||||
|
||||
|
||||
@ -124,8 +124,9 @@ bool test_redux(const Xpr&, int traversal, int unrolling)
|
||||
|
||||
template<typename Scalar, bool Enable = internal::packet_traits<Scalar>::Vectorizable> struct vectorization_logic
|
||||
{
|
||||
typedef internal::packet_traits<Scalar> PacketTraits;
|
||||
enum {
|
||||
PacketSize = internal::packet_traits<Scalar>::size
|
||||
PacketSize = PacketTraits::size
|
||||
};
|
||||
static void run()
|
||||
{
|
||||
@ -198,7 +199,7 @@ template<typename Scalar, bool Enable = internal::packet_traits<Scalar>::Vectori
|
||||
LinearTraversal,CompleteUnrolling));
|
||||
|
||||
VERIFY(test_assign(Matrix3(),Matrix3().cwiseQuotient(Matrix3()),
|
||||
LinearVectorizedTraversal,CompleteUnrolling));
|
||||
PacketTraits::HasDiv ? LinearVectorizedTraversal : LinearTraversal,CompleteUnrolling));
|
||||
|
||||
VERIFY(test_assign(Matrix<Scalar,17,17>(),Matrix<Scalar,17,17>()+Matrix<Scalar,17,17>(),
|
||||
LinearTraversal,NoUnrolling));
|
||||
@ -256,6 +257,7 @@ void test_vectorization_logic()
|
||||
|
||||
#ifdef EIGEN_VECTORIZE
|
||||
|
||||
CALL_SUBTEST( vectorization_logic<int>::run() );
|
||||
CALL_SUBTEST( vectorization_logic<float>::run() );
|
||||
CALL_SUBTEST( vectorization_logic<double>::run() );
|
||||
CALL_SUBTEST( vectorization_logic<std::complex<float> >::run() );
|
||||
|
||||
@ -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;
|
||||
@ -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);
|
||||
|
||||
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
@ -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
|
||||
@ -49,6 +53,129 @@ namespace Eigen
|
||||
knots.segment(knots.size()-degree-1,degree+1) = KnotVectorType::Ones(degree+1);
|
||||
}
|
||||
|
||||
/**
|
||||
* \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 (unsigned int i = 0; i < numAverageKnots - 1; ++i)
|
||||
{
|
||||
temporaryParameters[0] = averageKnots[i];
|
||||
ParameterVectorType parameterIndices(numParameters);
|
||||
int temporaryParameterIndex = 1;
|
||||
for (int 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 control points (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
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -13,23 +13,23 @@
|
||||
|
||||
namespace Eigen {
|
||||
|
||||
// lets do some explicit instantiations and thus
|
||||
// force the compilation of all spline functions...
|
||||
template class Spline<double, 2, Dynamic>;
|
||||
template class Spline<double, 3, Dynamic>;
|
||||
// lets do some explicit instantiations and thus
|
||||
// force the compilation of all spline functions...
|
||||
template class Spline<double, 2, Dynamic>;
|
||||
template class Spline<double, 3, Dynamic>;
|
||||
|
||||
template class Spline<double, 2, 2>;
|
||||
template class Spline<double, 2, 3>;
|
||||
template class Spline<double, 2, 4>;
|
||||
template class Spline<double, 2, 5>;
|
||||
template class Spline<double, 2, 2>;
|
||||
template class Spline<double, 2, 3>;
|
||||
template class Spline<double, 2, 4>;
|
||||
template class Spline<double, 2, 5>;
|
||||
|
||||
template class Spline<float, 2, Dynamic>;
|
||||
template class Spline<float, 3, Dynamic>;
|
||||
template class Spline<float, 2, Dynamic>;
|
||||
template class Spline<float, 3, Dynamic>;
|
||||
|
||||
template class Spline<float, 3, 2>;
|
||||
template class Spline<float, 3, 3>;
|
||||
template class Spline<float, 3, 4>;
|
||||
template class Spline<float, 3, 5>;
|
||||
template class Spline<float, 3, 2>;
|
||||
template class Spline<float, 3, 3>;
|
||||
template class Spline<float, 3, 4>;
|
||||
template class Spline<float, 3, 5>;
|
||||
|
||||
}
|
||||
|
||||
@ -234,11 +234,48 @@ void check_global_interpolation2d()
|
||||
}
|
||||
}
|
||||
|
||||
void check_global_interpolation_with_derivatives2d()
|
||||
{
|
||||
typedef Spline2d::PointType PointType;
|
||||
typedef Spline2d::KnotVectorType KnotVectorType;
|
||||
|
||||
const unsigned int numPoints = 100;
|
||||
const unsigned int dimension = 2;
|
||||
const unsigned int degree = 3;
|
||||
|
||||
ArrayXXd points = ArrayXXd::Random(dimension, numPoints);
|
||||
|
||||
KnotVectorType knots;
|
||||
Eigen::ChordLengths(points, knots);
|
||||
|
||||
ArrayXXd derivatives = ArrayXXd::Random(dimension, numPoints);
|
||||
VectorXd derivativeIndices(numPoints);
|
||||
|
||||
for (Eigen::DenseIndex i = 0; i < numPoints; ++i)
|
||||
derivativeIndices(i) = static_cast<double>(i);
|
||||
|
||||
const Spline2d spline = SplineFitting<Spline2d>::InterpolateWithDerivatives(
|
||||
points, derivatives, derivativeIndices, degree);
|
||||
|
||||
for (Eigen::DenseIndex i = 0; i < points.cols(); ++i)
|
||||
{
|
||||
PointType point = spline(knots(i));
|
||||
PointType referencePoint = points.col(i);
|
||||
VERIFY_IS_APPROX(point, referencePoint);
|
||||
PointType derivative = spline.derivatives(knots(i), 1).col(1);
|
||||
PointType referenceDerivative = derivatives.col(i);
|
||||
VERIFY_IS_APPROX(derivative, referenceDerivative);
|
||||
}
|
||||
}
|
||||
|
||||
void test_splines()
|
||||
{
|
||||
CALL_SUBTEST( eval_spline3d() );
|
||||
CALL_SUBTEST( eval_spline3d_onbrks() );
|
||||
CALL_SUBTEST( eval_closed_spline2d() );
|
||||
CALL_SUBTEST( check_global_interpolation2d() );
|
||||
for (int i = 0; i < g_repeat; ++i)
|
||||
{
|
||||
CALL_SUBTEST( eval_spline3d() );
|
||||
CALL_SUBTEST( eval_spline3d_onbrks() );
|
||||
CALL_SUBTEST( eval_closed_spline2d() );
|
||||
CALL_SUBTEST( check_global_interpolation2d() );
|
||||
CALL_SUBTEST( check_global_interpolation_with_derivatives2d() );
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user