Wrong evaluation of fourth derivative

[mmoelle1]

I just started using TinySpline as a reference implementation in my unit tests. I guess that there is a small error in the evaluation of the fourth derivative. Let the following be the JSON string representing the univariate spline

{
    "degree": 4,
    "dimension": 1,
    "control_points": [
        0,
        0.125,
        0.25,
        0.375,
        0.5,
        0.625,
        0.75,
        0.875,
        1
    ],
    "knots": [
        0,
        0,
        0,
        0,
        0,
        0.20000000000000001,
        0.40000000000000002,
        0.59999999999999998,
        0.80000000000000004,
        1,
        1,
        1,
        1,
        1
    ]
}

Then the line

(tinybspline.derive(4, -1).eval(0.2).result())[0]

evaluates to -729.167 but it should evaluate to -26.0417.

[msteinbeck]

Hi @mmoelle1,

thank you for your report. I'll have a look.

[msteinbeck]

If you derive a spline of degree 4 four times, you get a spline of degree 0 (which is basically a sequence of points; the sequence of points matches exactly the sequence of the control points of the spline). Let's have a look at the control points and knots of your derivative:

Control Points:

[-729.1666666666666, -26.041666666666647, 0.0, 26.04166666666664, 729.1666666666675]

Knots:

[0.0, 0.2, 0.4, 0.6, 0.8, 1.0]

You evaluate the spline at a split point (knot = 0.2)---that is, the spline is discontinuous at this knot. Split points are knots whose multiplicity, s, is equal to the order (degree + 1) of the spline. Because the degree of your derivative is 0, evaluating the spline at any knot with multiplicity 1 (as it is the case with knot 0.2) yields two results. If you access the second result as follows:

(tinybspline.derive(4, -1).eval(0.2).result())[1]

you will get the expected point -26.0417. This is documented in the C header file. Excerpt:

 * There is a special case in which the evaluation of a knot \c u returns two
 * results instead of one. It occurs when the multiplicity of \c u (\c s(u)) is
 * equals to the order of the evaluated spline, indicating that the spline is
 * discontinuous at \c u. This is common practice for B-Splines (and NURBS)
 * consisting of connected Bezier curves where the endpoint of curve \c c_i is
 * equal to the start point of curve \c c_i+1. Yet, the end point of \c c_i and
 * the start point of \c c_i+1 may still be completely different, yielding to
 * visible gaps (if distance of the points is large enough). In such case (\c
 * s(u) == \c order), ::ts_deboornet_points stores only the two resulting
 * points (there is no net to calculate) and ::ts_deboornet_result points to
 * the \e first point in ::ts_deboornet_points. Since having gaps in splines is
 * unusual, both points in ::ts_deboornet_points are generally equal, making it
 * easy to handle this special case by simply calling
 * ::ts_deboornet_result. However, one can access both points if necessary:
 *
 *     ts_deboornet_result(...)[0] ...       // Access the first component of
 *                                           // the first result.
 *
 *     ts_deboornet_result(...)[dim(spline)] // Access the first component of
 *                                           // the second result.
 *
 * As if this wasn't complicated enough, there is an exception for this special
 * case, yielding to exactly one result (just like the regular case) even if \c
 * s(u) == \c order. It occurs when \c u is the lower or upper bound of the
 * domain of the evaluated spline. For instance, if \c b is a spline with
 * domain [0, 1] and \c b is evaluated at \c u = \c 0 or \c u = \c 1, then
 * ::ts_deboornet_result is \e always a single point regardless of the
 * multiplicity of \c u.
 *
 * In summary, there are three different types of evaluation:
 *
 * 1. The regular case, in which all points of the net are returned.
 *
 * 2. A special case, in which two results are returned (required for splines
 * with gaps).
 *
 * 3. The exception of 2., in which exactly one result is returned (even if \c
 * s(u) == \c order).
 *
 * All in all this looks quite complex (and actually it is), but for most
 * applications you do not have to deal with this. Just use
 * ::ts_deboornet_result to access the outcome of De Boor's algorithm.

[msteinbeck]

@mmoelle1: Could I solve your issue? If yes, please close this issue.

[mmoelle1]

Yes, the issue is solved. Thanks!