Analitic BSpline basis

sympy_nurbs.py

import sympy as sp
import numpy as np
from matplotlib import pyplot as plt


class Knotvector:

    @classmethod
    def uniform(cls, degree, npts):
        ninter = npts - degree
        one = sp.Rational(1)
        start = [0*one for _ in range(degree)]
        end = [one for _ in range(degree)]
        middle = [i*one/ninter for i in range(ninter+1)]
        vector = start + middle + end
        return cls(vector, degree)

    def __init__(self, vector, degree=None):
        vector = tuple(v if not isinstance(v, int) else sp.Rational(v)
                       for v in vector)
        if degree is None:
            degree = 0
            while vector[degree] == vector[degree+1]:
                degree += 1
        npts = len(vector) - degree - 1

        knots = tuple(sorted(set(vector[degree:npts+1])))
        spans = [0] * len(knots)
        for i, knot in enumerate(knots):
            for j in range(degree, npts):
                if vector[j] == knot:
                    spans[i] = j
                    break
        spans[-1] = spans[-2]

        self.vector = vector
        self.degree = degree
        self.npts = npts
        self.knots = knots
        self.spans = spans

    def __getitem__(self, index):
        return self.vector[index]

    def __str__(self):
        return str(self.vector)


def spline_eval_matrix(knotvector, reqdegree=None):
    """
    Given a knotvector, it has properties like
        - number of points: npts
        - polynomial degree: degree
        - knots: A list of non-repeted knots
        - spans: The span of each knot
    This function returns a matrix of size
        (m) x (j+1) x (j+1)
    which
        - m is the number of segments: len(knots)-1
        - j is the requested degree
    """
    if reqdegree is None:
        reqdegree = knotvector.degree
    j = reqdegree
    knots = knotvector.knots
    spans = knotvector.spans
    niter = len(knots) - 1
    matrix = np.zeros((niter, j+1, j+1), dtype="object")
    if j == 0:
        val = knots[-1] - knots[0]
        matrix.fill(val/val)
        return matrix
    matrix_less1 = spline_eval_matrix(knotvector, j-1)
    for y in range(j):
        for z, sz in enumerate(spans[:-1]):
            i = y + sz - j + 1
            denom = knotvector[i + j] - knotvector[i]
            for k in range(j):
                matrix_less1[z, y, k] /= denom

            a0 = knots[z] - knotvector[i]
            a1 = knots[z + 1] - knots[z]
            b0 = knotvector[i + j] - knots[z]
            b1 = knots[z] - knots[z + 1]
            for k in range(j):
                matrix[z, y, k] += b0 * matrix_less1[z, y, k]
                matrix[z, y, k + 1] += b1 * matrix_less1[z, y, k]
                matrix[z, y + 1, k] += a0 * matrix_less1[z, y, k]
                matrix[z, y + 1, k + 1] += a1 * matrix_less1[z, y, k]
    return matrix


def compute_basis(knotvector: Knotvector, degree: int):
    npts = knotvector.npts
    shape = (npts, npts - degree)
    basis = np.zeros(shape, dtype="object")
    t = sp.symbols("t", real=True)

    matrix3d = spline_eval_matrix(knotvector)
    knots = knotvector.knots
    spans = knotvector.spans

    for j, span in enumerate(spans[:-1]):
        u = (t - knots[j])
        u /= knots[j + 1] - knots[j]
        for y in range(degree + 1):
            i = y + span - degree
            coefs = matrix3d[j, y]
            value = sum(ci * u**i for i, ci in enumerate(coefs))
            basis[i, j] = value

    for i in range(npts):
        for j in range(npts - degree):
            expr = sp.sympify(basis[i, j])
            expr = sp.expand(expr)
            expr = sp.simplify(expr)
            expr = sp.factor(expr)
            basis[i, j] = expr
    return basis


def print_basis(knotvector, basis):
    degree, npts = knotvector.degree, knotvector.npts
    print(f"Basis Function of degree {degree} and {npts} npts")
    print(f"Knotvector: {knotvector}")
    print(f"Expressions:")
    knots = knotvector.knots
    for j, (a, b) in enumerate(zip(knots, knots[1:])):
        print(f"    For interval {j}: [{a}, {b}]")
        for i in range(npts):
            print(f"        {i}: {basis[i, j]}")


def plot_basis(knotvector, basis):
    npts = knotvector.npts
    prop_cycle = plt.rcParams['axes.prop_cycle']
    default_colors = prop_cycle.by_key()['color']
    colors = []
    for i, color in enumerate(default_colors):
        if i == npts:
            break
        colors.append(color)

    tvar = sp.symbols("t", real=True)
    tvals = []
    allfvals = [[] for _ in range(npts)]
    knots = tuple(map(float, knotvector.knots))
    for j, (a, b) in enumerate(zip(knots, knots[1:])):
        tspace = np.linspace(a, b, 129)
        tvals += list(tspace)
        for i in range(npts):
            expr = basis[i, j]
            new_fvals = [expr.subs(tvar, ti) for ti in tspace]
            allfvals[i] += new_fvals

    for i, fvals in enumerate(allfvals):
        plt.plot(tvals, fvals, color=colors[i], label=str(i+1))
    plt.legend()
    plt.show()


def main():
    if True:
        degree, npts = 3, 6
        knotvector = Knotvector.uniform(degree, npts)
    else:
        degree = 2
        knotvector = [-3, -2, -1, 0, 1, 2, 3, 4, 5]
        knotvector = Knotvector(knotvector, degree)
    basis = compute_basis(knotvector, degree=degree)
    print_basis(knotvector, basis)
    plot_basis(knotvector, basis)
    

if __name__ == "__main__":
    main()

This gist translates the knotvector into a set of analitic basis functions, as shown bellow in the output.
Switching the condition at line 163 gives you different outputs, as you can see below:

Basis Function of degree 3 and 6 npts
Knotvector: (0, 0, 0, 0, 1/3, 2/3, 1, 1, 1, 1)
Expressions:
    For interval 0: [0, 1/3]
        0: -(3*t - 1)**3
        1: 9*t*(21*t**2 - 18*t + 4)/4
        2: -9*t**2*(11*t - 6)/4
        3: 9*t**3/2
        4: 0
        5: 0
    For interval 1: [1/3, 2/3]
        0: 0
        1: -(3*t - 2)**3/4
        2: 3*(21*t**3 - 36*t**2 + 18*t - 2)/4
        3: -3*(21*t**3 - 27*t**2 + 9*t - 1)/4
        4: (3*t - 1)**3/4
        5: 0
    For interval 2: [2/3, 1]
        0: 0
        1: 0
        2: -9*(t - 1)**3/2
        3: 9*(t - 1)**2*(11*t - 5)/4
        4: -9*(t - 1)*(21*t**2 - 24*t + 7)/4
        5: (3*t - 2)**3

Other condition:

Basis Function of degree 3 and 6 npts
Knotvector: (-5, -4, -3, -2, -1, 1, 2, 3, 4, 5)
Expressions:
    For interval 0: [-2, -1]
        0: -(t + 1)**3/6
        1: (9*t**3 + 33*t**2 + 27*t + 11)/24
        2: -(7*t**3 + 33*t**2 + 39*t - 1)/24
        3: (t + 2)**3/12
        4: 0
        5: 0
    For interval 1: [-1, 1]
        0: 0
        1: -(t - 1)**3/24
        2: (3*t**3 - 3*t**2 - 9*t + 11)/24
        3: -(3*t**3 + 3*t**2 - 9*t - 11)/24
        4: (t + 1)**3/24
        5: 0
    For interval 2: [1, 2]
        0: 0
        1: 0
        2: -(t - 2)**3/12
        3: (7*t**3 - 33*t**2 + 39*t + 1)/24
        4: -(9*t**3 - 33*t**2 + 27*t - 11)/24
        5: (t - 1)**3/6