How do I fix the bug in my b-spline implementation: The figure always has a point on 0,0 coordinate

Question

I have implemented Cox deboor's algorithm. However, my b-spline curve always has a point at 0,0 position, making the curve look strange. I am giving below code from scipy with a correct figure, followed by my own implementation which creates the wrong figure. The correct code below is borrowed from https://github.com/kawache/Python-B-spline-examples

    import numpy as np
    from scipy import interpolate

    import matplotlib.pyplot as plt
    Cp=torch.tensor([(3 , 1), (2.5, 4), (0, 1), (-2.5, 4),(-3, 0), (-2.5, -4), (0, -1), (2.5, -4), (3, -1)])


    x=Cp[:,0]
    y=Cp[:,1]
    l=len(x)
    t=np.linspace(0,1,l-2,endpoint=True)
    t=np.append([0,0,0],t)
    t=np.append(t,[1,1,1])
    tck=[t,[x,y],3]
    u3=np.linspace(0,1,(max(l*2,70)),endpoint=True)
    out = interpolate.splev(u3,tck) 
    plt.figure()
    plt.plot(x,y,'k--',label='Control polygon',marker='o',markerfacecolor='red')
    plt.plot(out[0],out[1],'b',linewidth=2.0,label='B-spline curve')
    plt.show()

My code (main recursive function)

    def N_i_p_vec(u,U,i,p):
        def deg_0(u,ui,ui_1):
            u2 = torch.where((ui<=u) &(u<ui_1),1,0)
            return u2
        if p>0:
            exp1 =torch.nan_to_num((u-U[i])/(U[i+p]-U[i])*N_i_p_vec(u,U,i,p-1),0)
            exp2 =torch.nan_to_num (((U[i+p+1]-u)/(U[i+p+1]-U[i+1]))*N_i_p_vec(u,U,i+1,p-1),0)
            return exp1+exp2
        if p==0:   
            return deg_0(u,U[i],U[i+1])

The loop that generates the figure.

  p=3
  m = len(t)-1
  num_basis = m-p-1
  o2 = torch.zeros(2,len(u3))
  for i in range(num_basis):
            N_m= N_i_p_vec(torch.tensor(u3),t,i,p)
            P_ij = Cp[i,:]
            c_tmp =torch.einsum('i,j->ji',N_m,P_ij) 
            o2+=c_tmp

Figure (wrong) itself:

    plt.figure()
    plt.plot(*o2,'b',linewidth=2.0,label='B-spline curve')
    plt.plot(x,y,'k--',label='Control polygon',marker='o',markerfacecolor='red')

UPDATE: Examining the output arrays from both my code and scipy code shows that the last entity in my output array o2 is 0,0, whereas it is 3,-1 for scipy array out

Answer 1

This is happening because the supporting knot interval for basis functions is exclusive at the end. Thus it makes the value of the last basis function at the last point, zero. To overcome this you can check if you are in the last interval, then you should make the interval inclusive at the end:

def in_last_interval(knot, x):
if len(knot) < 2:
    return False
return x >= knot[-2] and x <= knot[-1]

Then your code should be modified as

def N_i_p_vec(u,U,i,p):
    def deg_0(u,ui,ui_1):
        if not in_last_interval(U, u):
            return torch.where((ui<=u) & (u<ui_1),1,0)
        if in_last_interval(U,u):
            return torch.where((ui<=u) & (u<=ui_1),1,0)
    if p>0:
        exp1 =torch.nan_to_num((u-U[i])/(U[i+p]-U[i])*N_i_p_vec(u,U,i,p-1),0)
        exp2 =torch.nan_to_num (((U[i+p+1]-u)/(U[i+p+1]-U[i+1]))*N_i_p_vec(u,U,i+1,p-1),0)
        return exp1+exp2
    if p==0:   
        return deg_0(u,U[i],U[i+1])