diff --git a/unsupported/Eigen/src/Splines/Spline.h b/unsupported/Eigen/src/Splines/Spline.h index c7f8b4b08..74eca8e6b 100644 --- a/unsupported/Eigen/src/Splines/Spline.h +++ b/unsupported/Eigen/src/Splines/Spline.h @@ -28,6 +28,12 @@ namespace Eigen { * C(u) & = \sum_{i=0}^{n}N_{i,p}(u)P_i * \f} * + * The spline is parameterized over the domain \f$[u_p; u_{m-p}]\f$ of its knot + * vector \f$\{u_0,\hdots,u_m\}\f$, where \f$p\f$ is the degree. For the clamped + * knot vectors that SplineFitting produces this is \f$[u_0; u_m]\f$, which is + * the interval spanned by the interpolation parameters and therefore not + * necessarily \f$[0;1]\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 @@ -106,7 +112,7 @@ class Spline { * 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. + * \param u Parameter \f$u\f$ in the spline's knot domain at which the spline is evaluated. * \return The spline value at the given location \f$u\f$. **/ PointType operator()(Scalar u) const; @@ -120,7 +126,7 @@ class Spline { * \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 u Parameter \f$u\f$ in the spline's knot domain at which the spline derivative is evaluated. * \param order The order up to which the derivatives are computed. **/ typename SplineTraits::DerivativeType derivatives(Scalar u, DenseIndex order) const; @@ -147,8 +153,8 @@ class Spline { * 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. + * \param u Parameter \f$u\f$ in the spline's knot domain at which the non-zero + * basis functions are computed. **/ typename SplineTraits::BasisVectorType basisFunctions(Scalar u) const; @@ -161,8 +167,8 @@ class Spline { * \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 u Parameter \f$u\f$ in the spline's knot domain at which the non-zero + * basis function derivatives are computed. * \param order The order up to which the basis function derivatives are computed. **/ typename SplineTraits::BasisDerivativeType basisFunctionDerivatives(Scalar u, DenseIndex order) const; diff --git a/unsupported/Eigen/src/Splines/SplineFitting.h b/unsupported/Eigen/src/Splines/SplineFitting.h index b0d5ad3fd..bb39872ff 100644 --- a/unsupported/Eigen/src/Splines/SplineFitting.h +++ b/unsupported/Eigen/src/Splines/SplineFitting.h @@ -31,12 +31,13 @@ namespace Eigen { * * The knots are computed as * \f{align*} - * u_0 & = \hdots = u_p = 0 \\ - * u_{m-p} & = \hdots = u_{m} = 1 \\ + * u_0 & = \hdots = u_p = \bar{u}_0 \\ + * u_{m-p} & = \hdots = u_{m} = \bar{u}_n \\ * u_{j+p} & = \frac{1}{p}\sum_{i=j}^{j+p-1}\bar{u}_i \quad\quad j=1,\hdots,n-p * \f} - * where \f$p\f$ is the degree and \f$m+1\f$ the number of knots - * of the desired interpolating spline. + * where \f$p\f$ is the degree, \f$m+1\f$ the number of knots of the desired + * interpolating spline, and \f$\bar{u}_0,\hdots,\bar{u}_n\f$ the input + * parameters, which may span any interval, not just \f$[0,1]\f$. * * \param[in] parameters The input parameters. During interpolation one for each data point. * \param[in] degree The spline degree which is used during the interpolation. @@ -50,8 +51,10 @@ void KnotAveraging(const KnotVectorType& parameters, DenseIndex degree, KnotVect for (DenseIndex j = 1; j < parameters.size() - degree; ++j) knots(j + degree) = parameters.segment(j, degree).mean(); - knots.segment(0, degree + 1) = KnotVectorType::Zero(degree + 1); - knots.segment(knots.size() - degree - 1, degree + 1) = KnotVectorType::Ones(degree + 1); + // The boundary knots replicate the first and last parameter so that the + // spline domain matches the parameter range, whatever interval it spans. + knots.head(degree + 1).setConstant(parameters(0)); + knots.tail(degree + 1).setConstant(parameters(placeholders::last)); } /** diff --git a/unsupported/test/splines.cpp b/unsupported/test/splines.cpp index 89674c31a..21fe4239f 100644 --- a/unsupported/test/splines.cpp +++ b/unsupported/test/splines.cpp @@ -165,6 +165,43 @@ void check_global_interpolation2d() { } } +// Regression test for issue #764: KnotAveraging hard-coded 0/1 boundary knots, +// producing a non-monotone knot vector for parameters outside [0, 1]. +void check_global_interpolation2d_arbitrary_range() { + typedef Spline2d::PointType PointType; + typedef Spline2d::KnotVectorType KnotVectorType; + typedef Spline2d::ControlPointVectorType ControlPointVectorType; + + // The issue's example: parameters [2,3,4,5] with degree 1 must yield the + // clamped knot vector [2,2,3,4,5,5], not [0,0,3,4,1,1]. + { + KnotVectorType parameters(4), knots, expected(6); + parameters << 2, 3, 4, 5; + expected << 2, 2, 3, 4, 5, 5; + Eigen::KnotAveraging(parameters, 1, knots); + VERIFY_IS_EQUAL(knots.size(), expected.size()); + VERIFY((knots - expected).matrix().norm() == 0.0); + } + + // Interpolation with knot parameters spanning [2, 5] must pass through the + // data points, exactly like the normalized [0, 1] parameterization. + { + ControlPointVectorType points = ControlPointVectorType::Random(2, 100); + + KnotVectorType chord_lengths; + Eigen::ChordLengths(points, chord_lengths); + KnotVectorType parameters = 2.0 + 3.0 * chord_lengths; + + const Spline2d spline = SplineFitting::Interpolate(points, 3, parameters); + + for (Eigen::DenseIndex i = 0; i < points.cols(); ++i) { + PointType pt = spline(parameters(i)); + PointType ref = points.col(i); + VERIFY((pt - ref).matrix().norm() < 32 * NumTraits::epsilon()); + } + } +} + void check_global_interpolation_with_derivatives2d() { typedef Spline2d::PointType PointType; typedef Spline2d::KnotVectorType KnotVectorType; @@ -202,6 +239,7 @@ EIGEN_DECLARE_TEST(splines) { CALL_SUBTEST(eval_spline3d_onbrks()); CALL_SUBTEST(eval_closed_spline2d()); CALL_SUBTEST(check_global_interpolation2d()); + CALL_SUBTEST(check_global_interpolation2d_arbitrary_range()); CALL_SUBTEST(check_global_interpolation_with_derivatives2d()); } }