robotics : 5RI2i

B-Splines And NURBS

Mathematics and code for evaluating B-splines and NURBS curves.

Jan 28, 2024
Markdown
Updated Feb 3, 2024

Introduction

B-splines and NURBS curves are two closely related types of curves defined by a linear combination of basis functions. In fact, B-splines are a subset to NURBS. These curves can be used to describe very intricate and complex trajectories while not taking up a lot of memory.

This post is based on the 1995 book “The NURBS Book” by Les Piegl and Wayne Tiller. Almost all of the mathematics, definitions etc are taken from this book. One exception is De Boor’s algorithm which can be found by simply searching for it, this post does not explain how the algorithm works.

As always, the code presented in this post is ANSI C, not all of the sub-routines are explained or shown here. Visit the repository at github.com/fullnitrous/arm, all relevant code is found at /src/core/.

B-splines

A B-spline curve \(C(u)\) of degree \(p\) with \(n\) control points \(P_i\) is defined as the following.

\[ C(u) = \sum_{i=0}^n N_{i,p}(u)P_i, \quad a \leq u \leq b \]

The definitions is reformulated such that the parameter \(u\) is denoted by \(t\) within the interval \([0,1]\). The degree \(p\) is denoted by \(d\).

\[ C(t) = \sum_{i=0}^n N_{i,d}(t)P_i, \quad i \in [0, 1] \]

The term \(N_{i,d}(t)\) is the \(d\)th-degree B-spline basis function defined on the knot vector \(U\) where \(u_i\) denotes elements in the knot vector.

The type bspline_t is used to store all of the information for a B-spline curve. The struct contains an array of control points p, the knot vector u, degree d and number of control points n.

Due to how the evaluation function works, a maximum degree for the B-spline is spesified in a macro.

The helper function eval_bspline_init takes in the arrays and creates the bspline struct for ease of use.

Approach for Evaluation

There are multiple methods of evaluating \(C(t)\) and its derivatives. A naive approach such as evaluating \(C(t)\) by its equation and the recursive formulas for \(N_{i,d}(t)\) would prove to be too slow for a real-time application.

The approach selected here is to evaluate \(C(t)\) using De Boor’s algorithm which is a numerically stable and efficient. De Boor’s algorithm is also modified to calculate the control points \(Q_i\) for the first and second derivative of \(C(t)\). This way the algorithm simultaneously calculates \(C(t)\), \(C''(t)\) and \(C'''(t)\).

The derivative of a B-spline curve \(C(t)\) of degree \(d\) is another B-spline curve \(C'(t)\) of degree \(d-1\) with new control points \(Q_i\).

The new control points are defined as the following.

\[ Q_i = d\frac{P_{i+1}-P_i}{u_{i+d+i}-u_{i+1}} \]

The derivative \(C'(t)\) then becomes.

\[ C'(t) = \sum_{i=0}^{n-1} N_{i,d-1}(t)Q_i, \quad t \in [0,1] \]

The steps for calculating \(C'(t)\) can be repeated recursively for higher order derivatives.

In order to evaluate the curve, the function eval_bspline is used with the bspline_t passed to f and the parameter value t. This function returns an eval_t with the position and the first and second derivative and the magnitude of the second derivative as well. The routine simply just returns the value from __eval_bspline which actually evaluates the values.

As with the splines in the post about interpolation, the knot vector is searched to find the span k. However, the function eval_findinterval does not work with the knot vector since it is similar but not the same as an array with parameter values. Therefore, the function eval_findspan is used instead, this function is also specified in the referenced book.

The function __eval_bspline below is used to do a complete evaluation of the B-spline given the function f, span k and parameter t.

eval_t __eval_bspline(void* f_, int k, double t) {
    bspline_t* f;
    eval_t e;
    double alpha, tmp, *u;
    int i, j, r, d, n, left, right;

    vec3_t p[EVAL_BSPLINE_MAX_DEG+1];
    vec3_t v[EVAL_BSPLINE_MAX_DEG+0];
    vec3_t a[EVAL_BSPLINE_MAX_DEG-1];

    f = (bspline_t*)f_;
    u = f->u;
    d = f->d;
    n = f->n;
    t = u[d] + t * (u[n] - u[d]);


    for(i = 0; i <= d; i++) {
        p[i] = f->p[i + k - d];
    }

    for(i = 0; i <= d-1; i++) {
        tmp    = d / (u[i+k+1] - u[i+k-d+1]);
        v[i].x = tmp*(p[i+1].x - p[i].x);
        v[i].y = tmp*(p[i+1].y - p[i].y);
        v[i].z = tmp*(p[i+1].z - p[i].z);
    }

    for(i = 0; i <= d - 2; i++) {
        tmp    = (d-1) / (u[i+k+1] - u[i+k-d+2]);
        a[i].x = tmp*(v[i+1].x - v[i].x);
        a[i].y = tmp*(v[i+1].y - v[i].y);
        a[i].z = tmp*(v[i+1].z - v[i].z);
    }

    for(r = 1; r <= d; r++) {
        for(j = d; j >= r; j--) {
            right = j+1+k-r;
            left  = j+k-d;
            alpha = (t-u[left])/(u[right]-u[left]);
            left++;

            p[j].x = (1.0-alpha)*p[j-1].x + alpha*p[j].x;
            p[j].y = (1.0-alpha)*p[j-1].y + alpha*p[j].y;
            p[j].z = (1.0-alpha)*p[j-1].z + alpha*p[j].z;

            if(r <= d-1 && j <= d-1) {
                alpha  = (t-u[left])/(u[right]-u[left]);
                left++;
                v[j].x = (1.0-alpha)*v[j-1].x + alpha*v[j].x;
                v[j].y = (1.0-alpha)*v[j-1].y + alpha*v[j].y;
                v[j].z = (1.0-alpha)*v[j-1].z + alpha*v[j].z;
            }

            if(r <= d-2 && j <= d-2) {
                alpha  = (t-u[left])/(u[right]-u[left]);
                a[j].x = (1.0-alpha)*a[j-1].x + alpha*a[j].x;
                a[j].y = (1.0-alpha)*a[j-1].y + alpha*a[j].y;
                a[j].z = (1.0-alpha)*a[j-1].z + alpha*a[j].z;
            }
        }
    }

    e.p = p[d];

    if(d >= 1) { e.v = v[d-1];     }
    else       { e.v = veczeros(); }

    if(d >= 2) { e.a = a[d-2];     }
    else       { e.a = veczeros(); }

    e.m = vecmag(e.v);

    return e;
}

However, in the case of wanting to only evaluate the position, the function __bspline_fast does exactly that.

B-spline curve, n=6, d=3, Control polygon shown in dashed line. Control points P_i and knot vector U listed in appendix.
B-spline curve, \(n=6\), \(d=3\), Control polygon shown in dashed line. Control points \(P_i\) and knot vector \(U\) listed in appendix.

NURBS

A Non-Uniform Rational B-spline (NURBS) curve \(C(u)\) of degree \(p\) with control points \(P_i\) and corresponding weights \(w_i\) is defined as the following.

\[ C(u) = \frac{\sum_{i=0}^n N_{i,p}(u)w_iP_i}{\sum_{i=0}^nN_{i,p}(u)w_i}, \quad a \leq u \leq b \]

Once again, this definition is reformulated such that the parameter \(u\) is denoted by \(t\) within the interval \([0, 1]\). The degree \(p\) is denoted by \(d\).

\[ C(t) = \frac{\sum_{i=0}^n N_{i,p}(t)w_iP_i}{\sum_{i=0}^nN_{i,p}(t)w_i}, \quad t \in [0,1] \]

The term \(N_{i,d}(t)\) and the knot vector \(U\) is the same as for B-splines.

The type nurbs_t is used to store all of the information for a NURBS curve and is identical to bspline_t besides that the array of points p is an array of four-dimensional vectors.

As with the B-spline curve, the maximum degree is defined.

And a intialization helper function for the nurbs_t type.

Approach for Evaluation

Similarly to B-splines, the approach selected to evaluate NURBS curves is to use De Boor’s algorith. Although in order to evaluate a NURBS curve using De Boor’s algorith, the three-dimensional control points \(P_i\) and weights \(w_i\) have to be converted to homogeneous coordinates such that De Boor’s algorithm is applied on a B-spline curve with four-dimensional control points. The homogeneous coordinates / four-dimensional control points are defined as.

\[ \begin{aligned} H_{i,x} &= w_iP_{i,x} \\ H_{i,y} &= w_iP_{i,y} \\ H_{i,z} &= w_iP_{i,z} \\ H_{i,w} &= w_i \end{aligned} \]

The NURBS curve is then evaluated using \(H\) as the control points for a B-spline curve in four dimensions. After De Boor’s algorithm has been applied to the four-dimensional B-spline curve, the result must be projected back to three dimensions.

\[ \begin{aligned} x_\text{final} &= x_\text{result} / w_\text{result} \\ y_\text{final} &= y_\text{result} / w_\text{result} \\ z_\text{final} &= z_\text{result} / w_\text{result} \\ \end{aligned} \]

The \(k\)th derivative of a NURBS curve is defined as the following

\[ C^{(k)}(t) = \frac{A^{(k)}(t) - \sum_{i=1}^k \begin{pmatrix}k\\i\end{pmatrix}w^{(i)}(t)C^{(k-i)}(t)}{w(t)} \]

Where the \(\begin{pmatrix}k\\i\end{pmatrix}\) is the binomial coefficient. The term \(A^{(k)}(t)\) is a vector from the first three components after evaluating the curve, not having projected it back to three dimensions. The term \(w^{(i)}(t)\) is the same but only for the weight component. Essentially they are this.

\[ \begin{aligned} A^{(k)}(t) &= \begin{bmatrix}x_\text{result} \\ y_\text{result} \\ z_\text{result} \end{bmatrix} \\ w^{(k)}(t) &= w_\text{result} \end{aligned} \]

Where \(x_\text{result}\), \(y_\text{result}\), \(z_\text{result}\) and \(w_\text{result}\), depending on \(k\), are the result of evaluating the four-dimensional B-spline curve or its derivatives. Using the aforementioned equation for the \(k\)th derivative. The first and second derivatives become.

\[ \begin{aligned} C'(t) &= \frac{A'(t)-w'(t)C(t)}{w(t)} \\ C''(t) &= \frac{A''(t)-2w'(t)C''(t)-w''(t)C(t)}{w(t)} \end{aligned} \]

The family of functions for evaluating a NURBS curve are pretty much the same as for a B-spline curve. The function eval_nurbs does a full evaluation using __eval_nurbs and __eval_nurbs_fast is used to only calculate the position of the curve.

eval_t __eval_nurbs(void* f_, int k, double t) {
    nurbs_t* f;
    eval_t e;
    double alpha, tmp, *u;
    int i, j, r, d, n, left, right;

    vec4_t p[EVAL_NURBS_MAX_DEG+1];
    vec4_t v[EVAL_NURBS_MAX_DEG+0];
    vec4_t a[EVAL_NURBS_MAX_DEG-1];

    f = (nurbs_t*)f_;
    u = f->u;
    d = f->d;
    n = f->n;
    t = u[d] + t * (u[n] - u[d]);

    for(i = 0; i <= d; i++) {
        vec4_t v = f->p[i+k-d];
        p[i].x   = v.x * v.w;
        p[i].y   = v.y * v.w;
        p[i].z   = v.z * v.w;
        p[i].w   = v.w;
    }

    for(i = 0; i <= d-1; i++) {
        tmp    = d / (u[i+k+1] - u[i+k-d+1]);
        v[i].x = tmp*(p[i+1].x - p[i].x);
        v[i].y = tmp*(p[i+1].y - p[i].y);
        v[i].z = tmp*(p[i+1].z - p[i].z);
        v[i].w = tmp*(p[i+1].w - p[i].w);
    }

    for(i = 0; i <= d - 2; i++) {
        tmp    = (d-1) / (u[i+k+1] - u[i+k-d+2]);
        a[i].x = tmp*(v[i+1].x - v[i].x);
        a[i].y = tmp*(v[i+1].y - v[i].y);
        a[i].z = tmp*(v[i+1].z - v[i].z);
        a[i].w = tmp*(v[i+1].w - v[i].w);
    }

    for(r = 1; r <= d; r++) {
        for(j = d; j >= r; j--) {
            right = j+1+k-r;
            left  = j+k-d;
            alpha = (t-u[left])/(u[right]-u[left]);
            left++;

            p[j].x = (1.0-alpha)*p[j-1].x + alpha*p[j].x;
            p[j].y = (1.0-alpha)*p[j-1].y + alpha*p[j].y;
            p[j].z = (1.0-alpha)*p[j-1].z + alpha*p[j].z;
            p[j].w = (1.0-alpha)*p[j-1].w + alpha*p[j].w;

            if(r <= d-1 && j <= d-1) {
                alpha  = (t-u[left])/(u[right]-u[left]);
                left++;
                v[j].x = (1.0-alpha)*v[j-1].x + alpha*v[j].x;
                v[j].y = (1.0-alpha)*v[j-1].y + alpha*v[j].y;
                v[j].z = (1.0-alpha)*v[j-1].z + alpha*v[j].z;
                v[j].w = (1.0-alpha)*v[j-1].w + alpha*v[j].w;
            }

            if(r <= d-2 && j <= d-2) {
                alpha  = (t-u[left])/(u[right]-u[left]);
                a[j].x = (1.0-alpha)*a[j-1].x + alpha*a[j].x;
                a[j].y = (1.0-alpha)*a[j-1].y + alpha*a[j].y;
                a[j].z = (1.0-alpha)*a[j-1].z + alpha*a[j].z;
                a[j].w = (1.0-alpha)*a[j-1].w + alpha*a[j].w;
            }
        }
    }

    p[d].x /= p[d].w;
    p[d].y /= p[d].w;
    p[d].z /= p[d].w;

    e.p.x = p[d].x;
    e.p.y = p[d].y;
    e.p.z = p[d].z;

    if(d >= 1) {
        v[d-1].x = (v[d-1].x - v[d-1].w*p[d].x) / p[d].w;
        v[d-1].y = (v[d-1].y - v[d-1].w*p[d].y) / p[d].w;
        v[d-1].z = (v[d-1].z - v[d-1].w*p[d].z) / p[d].w;
        e.v.x    = v[d-1].x;
        e.v.y    = v[d-1].y;
        e.v.z    = v[d-1].z;
    } else {
        e.v = veczeros();
    }

    if(d >= 2) {
        e.a.x = (a[d-2].x - 2.0*v[d-1].w*v[d-1].x - a[d-2].w*p[d].x) / p[d].w;
        e.a.y = (a[d-2].y - 2.0*v[d-1].w*v[d-1].y - a[d-2].w*p[d].y) / p[d].w;
        e.a.z = (a[d-2].z - 2.0*v[d-1].w*v[d-1].z - a[d-2].w*p[d].z) / p[d].w;
    } else {
        e.v = veczeros();
    }

    e.m = vecmag(e.v);

    return e;
}
NURBS curve, n=51, d=3, Control polygon hidden. Control points P_i and weights w_i represented as four-dimensional vectors with knot vector U in appendix.
NURBS curve, \(n=51\), \(d=3\), Control polygon hidden. Control points \(P_i\) and weights \(w_i\) represented as four-dimensional vectors with knot vector \(U\) in appendix.

Appendix

B-spline

The control points in the figure are the following.

Where the following is the knot vector.

NURBS

The control points in the figure are the following.

54.493000, 52.139000, 0.000000, 1.000000
55.507000, 52.139000, 0.000000, 1.000000
56.082000, 49.615000, 0.000000, 1.000000
56.780000, 44.971000, 0.000000, 1.200000
69.575000, 51.358000, 0.000000, 1.000000
77.786000, 58.573000, 0.000000, 1.000000
90.526000, 67.081000, 0.000000, 1.000000
105.973000, 63.801000, 0.000000, 1.000000
100.400000, 47.326000, 0.000000, 1.000000
94.567000, 39.913000, 0.000000, 1.000000
92.369000, 30.485000, 0.000000, 1.000000
83.440000, 33.757000, 0.000000, 2.000000
91.892000, 28.509000, 0.000000, 1.000000
89.444000, 20.393000, 0.000000, 1.000000
83.218000, 15.446000, 0.000000, 5.000000
87.621000, 4.830000, 0.000000, 3.000000
80.945000, 9.267000, 0.000000, 1.000000
79.834000, 14.535000, 0.000000, 1.100000
76.074000, 8.522000, 0.000000, 1.000000
70.183000, 12.550000, 0.000000, 1.000000
64.171000, 16.865000, 0.000000, 1.000000
59.993000, 22.122000, 0.000000, 1.000000
55.680000, 36.359000, 0.000000, 1.000000
56.925000, 24.995000, 0.000000, 1.000000
59.765000, 19.828000, 0.000000, 1.000000
54.493000, 14.940000, 0.000000, 1.000000
49.220000, 19.828000, 0.000000, 1.000000
52.060000, 24.994000, 0.000000, 1.000000
53.305000, 36.359000, 0.000000, 1.000000
48.992000, 22.122000, 0.000000, 1.000000
44.814000, 16.865000, 0.000000, 1.000000
38.802000, 12.551000, 0.000000, 1.000000
32.911000, 8.521000, 0.000000, 1.000000
29.152000, 14.535000, 0.000000, 1.100000
28.040000, 9.267000, 0.000000, 1.000000
21.364000, 4.830000, 0.000000, 3.000000
25.768000, 15.447000, 0.000000, 5.000000
19.539000, 20.391000, 0.000000, 1.000000
17.097000, 28.512000, 0.000000, 1.000000
25.537000, 33.750000, 0.000000, 2.000000
16.602000, 30.496000, 0.000000, 1.000000
14.199000, 39.803000, 0.000000, 1.000000
8.668000, 47.408000, 0.000000, 1.000000
3.000000, 63.794000, 0.000000, 1.000000
18.465000, 67.084000, 0.000000, 1.000000
31.197000, 58.572000, 0.000000, 1.000000
39.411000, 51.358000, 0.000000, 1.000000
52.204000, 44.971000, 0.000000, 1.200000
52.904000, 49.614000, 0.000000, 1.000000
53.478000, 52.139000, 0.000000, 1.000000
54.492000, 52.139000, 0.000000, 1.000000

Where the knot vector is the following.

The knot vector and control points are taken from a research paper I read but I do not remember the name of the paper.