Source code for torx.analysis.bspline_m

"""General B-spline basis evaluation utilities."""
import numpy as np
from numba import njit



@njit(cache=True)
def bspline_basis_and_deriv1(u, k, t):
    """
    Return (N, dN) of length k+1: B-spline basis and 1st derivatives at u.

    k is the degree, t is the knot vector.
    N[r] = B_{i-k+r}^k(u), dN[r] = dB_{i-k+r}^k/du, where i is the knot span.
    Also returns i so callers can map basis values to spline coefficients.
    """
    n = len(t) - k - 1

    # Find knot span: largest i s.t. t[i] <= u < t[i+1]
    if u >= t[n]:
        i = n - 1
    elif u <= t[k]:
        i = k
    else:
        lo, hi = k, n
        mid = (lo + hi) // 2
        while u < t[mid] or u >= t[mid + 1]:
            if u < t[mid]:
                hi = mid
            else:
                lo = mid
            mid = (lo + hi) // 2
        i = mid

    N = np.zeros(k + 1)
    N[0] = 1.0
    left = np.empty(k + 1)
    right = np.empty(k + 1)
    for j in range(1, k + 1):
        left[j] = u - t[i + 1 - j]
        right[j] = t[i + j] - u
        saved = 0.0
        for r in range(j):
            denom = right[r + 1] + left[j - r]
            temp = N[r] / denom if denom != 0.0 else 0.0
            N[r] = saved + right[r + 1] * temp
            saved = left[j - r] * temp
        N[j] = saved

    # Degree-(k-1) basis for derivative formula
    Nm1 = np.zeros(k)
    Nm1[0] = 1.0
    for j in range(1, k):
        left[j] = u - t[i + 1 - j]
        right[j] = t[i + j] - u
        saved = 0.0
        for r in range(j):
            denom = right[r + 1] + left[j - r]
            temp = Nm1[r] / denom if denom != 0.0 else 0.0
            Nm1[r] = saved + right[r + 1] * temp
            saved = left[j - r] * temp
        Nm1[j] = saved

    dN = np.zeros(k + 1)
    for r in range(k + 1):
        lv = 0.0
        if r > 0:
            d = t[i + r] - t[i - k + r]
            if d != 0.0:
                lv = k * Nm1[r - 1] / d
        rv = 0.0
        if r < k:
            d = t[i + r + 1] - t[i - k + r + 1]
            if d != 0.0:
                rv = k * Nm1[r] / d
        dN[r] = lv - rv

    return N, dN, i