Splines: Let KnotAveraging accept parameters outside the unit interval

libeigen/eigen!2842

Closes #764

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

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<Spline>::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<Spline>::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<Spline>::BasisDerivativeType basisFunctionDerivatives(Scalar u, DenseIndex order) const;

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));
 }
 
 /**

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<Spline2d>::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<double>::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());
   }
 }