msteinbeck / tinyspline

ANSI C library for NURBS, B-Splines, and Bézier curves with interfaces for C++, C#, D, Go, Java, Javascript, Lua, Octave, PHP, Python, R, and Ruby.
MIT License
1.21k stars 209 forks source link

Artifacts when evaluating seven point spline #239

Closed jay-sridharan closed 8 months ago

jay-sridharan commented 8 months ago

Hello again!

I came across a pretty interesting issue where I am seeing artifacts around the center of a seven point spline. The fifth-degree spline is formulated with seven control points. The details can be seen in the example scripts I've attached, but here is the symptom:

Screenshot from 2024-03-22 12-50-36

The x-axis on the plots are "time" - the spline is evaluated at each successive value every 1ms.

Anyway, the spline is moving along smoothly and then around u = 0.5, there is a discontinuity. This plot was generated using the artifact.c example. Here is the same data in table form:

u x y
0.499531 -127.507957 -189.170307
0.499569 -127.507462 -189.170394
0.499607 -127.506967 -189.170481
0.499644 -127.506484 -189.170566
0.499682 -127.505989 -189.170654
0.499720 -127.505494 -189.170741
0.499758 -127.504999 -189.170829
0.499796 -127.504504 -189.170916
0.499833 -127.504022 -189.171001
0.499871 -127.503526 -189.171088
0.499909 -127.501846 -189.171385
0.499947 -127.501846 -189.171385
0.499985 -127.501846 -189.171385
0.500022 -127.501846 -189.171385
0.500060 -127.501846 -189.171385
0.500098 -127.501846 -189.171385
0.500136 -127.500073 -189.171697
0.500174 -127.499578 -189.171784
0.500211 -127.499096 -189.171869
0.500249 -127.498601 -189.171957
0.500287 -127.498106 -189.172044
0.500325 -127.497610 -189.172131
0.500362 -127.497128 -189.172216
0.500400 -127.496633 -189.172303

Note the repeated values (-189.171385)

I constructed the same spline with python (see artifact.py) and get smooth results

u x y
0.499531 -127.507957 -189.170307
0.499569 -127.507462 -189.170394
0.499607 -127.506967 -189.170481
0.499644 -127.506484 -189.170566
0.499682 -127.505989 -189.170654
0.499720 -127.505494 -189.170741
0.499758 -127.504999 -189.170829
0.499796 -127.504504 -189.170916
0.499833 -127.504022 -189.171001
0.499871 -127.503526 -189.171088
0.499909 -127.503031 -189.171176
0.499947 -127.502536 -189.171263
0.499985 -127.502041 -189.171350
0.500022 -127.501559 -189.171435
0.500060 -127.501064 -189.171523
0.500098 -127.500568 -189.171610
0.500136 -127.500073 -189.171697
0.500174 -127.499578 -189.171784
0.500211 -127.499096 -189.171869
0.500249 -127.498601 -189.171957
0.500287 -127.498106 -189.172044
0.500325 -127.497610 -189.172131
0.500362 -127.497128 -189.172216
0.500400 -127.496633 -189.172303

You can see that the issue is only in the range where u is between 0.499909 and 0.500098, but otherwise all values match!

Here are the two scripts I used to find this issue:

artifact.c

#include "tinyspline.h"

#include <stdlib.h>
#include <stdio.h>

tsBSpline spline;

void setup()
{
    tsReal *ctrlp;

    size_t SPLINE_NUM_CTRLP    = 7;
    size_t SPLINE_DIMENSION    = 3;
    size_t SPLINE_DEGREE       = 5;

    ts_bspline_new(
        SPLINE_NUM_CTRLP,
        SPLINE_DIMENSION,
        SPLINE_DEGREE,
        TS_CLAMPED, /* used to hit first and last control point */
        &spline, /* the spline to setup */
        NULL);

    ts_bspline_control_points(&spline, &ctrlp, NULL);

    ctrlp[ 0] = -133.946950; ctrlp[ 1] = -187.623500; ctrlp[ 2] = 40.000000;
    ctrlp[ 3] = -132.663960; ctrlp[ 4] = -187.967280; ctrlp[ 5] = 40.000000;
    ctrlp[ 6] = -130.097980; ctrlp[ 7] = -188.654840; ctrlp[ 8] = 40.000000;
    ctrlp[ 9] = -127.532000; ctrlp[10] = -189.342400; ctrlp[11] = 40.000000;
    ctrlp[12] = -124.885608; ctrlp[13] = -189.573919; ctrlp[14] = 40.000000;
    ctrlp[15] = -122.239216; ctrlp[16] = -189.805439; ctrlp[17] = 40.000000;
    ctrlp[18] = -120.916020; ctrlp[19] = -189.921198; ctrlp[20] = 40.000000;

    ts_bspline_set_control_points(&spline, ctrlp, NULL);

    free(ctrlp);
}

void tear_down()
{
    ts_bspline_free(&spline);
}

int main(int argc, char** argv)
{
    setup();

    double u[24] = {
        0.499531, 0.499569, 0.499607, 0.499644, 0.499682, 0.499720, 0.499758,
        0.499796, 0.499833, 0.499871, 0.499909, 0.499947, 0.499985,
        0.500022, 0.500060, 0.500098, 0.500136, 0.500174, 0.500211,
        0.500249, 0.500287, 0.500325, 0.500362, 0.500400
    };

    tsReal *result;
    tsDeBoorNet net;

    for(int i = 0; i < 24; i++)
    {
        ts_bspline_eval(&spline, u[i], &net, NULL);
        ts_deboornet_result(&net, &result, NULL);

        printf("u = %f, x = %f, y = %f\n", u[i], result[0], result[1]);

        ts_deboornet_free(&net);
        free(result);
    }

    tear_down();
    return 0; 
}

artifact.py

from scipy import interpolate
import numpy as np

if __name__ == "__main__":

    control_points = np.array([
        np.array([-133.946950, -187.623500, 40.000000]),
        np.array([-132.663960, -187.967280, 40.000000]),
        np.array([-130.097980, -188.654840, 40.000000]),
        np.array([-127.532000, -189.342400, 40.000000]),
        np.array([-124.885608, -189.573919, 40.000000]),
        np.array([-122.239216, -189.805439, 40.000000]),
        np.array([-120.916020, -189.921198, 40.000000])
    ])

    T = [0, 0, 0, 0, 0, 0, 0.5, 1, 1, 1, 1, 1, 1] # clamped knot vector
    k = 5 # 5th degree

    c_x = control_points[:, 0]
    c_y = control_points[:, 1]
    c_z = control_points[:, 2]

    spl_x = interpolate.BSpline(T, c_x, k)
    spl_y = interpolate.BSpline(T, c_y, k)
    spl_z = interpolate.BSpline(T, c_z, k)

    u = [
        0.499531, 0.499569, 0.499607, 0.499644, 0.499682, 0.499720, 0.499758,
        0.499796, 0.499833, 0.499871, 0.499909, 0.499947, 0.499985,
        0.500022, 0.500060, 0.500098, 0.500136, 0.500174, 0.500211,
        0.500249, 0.500287, 0.500325, 0.500362, 0.500400
    ]

    for i in range(0, len(u)):

        x = spl_x(u[i]) #[()]
        y = spl_y(u[i]) #[()]
        z = spl_z(u[i]) #[()]

        print("u = %f, x = %f, y = %f" % (u[i], x, y))

I tried debugging myself, but got a bit lost when it came to the DeBoor algorithm. Hoping you can spot this relatively quickly! Thanks again for the great library!

Jay

msteinbeck commented 8 months ago

@jay-sridharan,

My first guess is that your knot values are too small: https://github.com/msteinbeck/tinyspline/blob/master/src/tinyspline.h#L156 Maybe the knot supplied to the evaluation function is not the knot eventually used for evaluation: https://github.com/msteinbeck/tinyspline/blob/master/src/tinyspline.c#L1222

jay-sridharan commented 8 months ago

Yes, it looks like this was the issue! Thank you for the help - I just reduced the max knots and made the epsilon smaller.