I have implemented a cubic B-Spline interpolation, not approximation, as follows: import numpy as np import math from geomdl import knotvector def cox_de_boor( d_, t_, k_, knots_): if (d_ == 0): if ( knots_[k_] <= t_ <= knots_[k_+1]): return 1.0 return 0.0 denom_l = (knots_[k_+d_] - knots_[k_]) left = 0.0 if (denom_l != 0.0): left = ((t_ - knots_[k_]) / denom_l) * cox_de_boor(d_-1, t_, k_, knots_) denom_r = (knots_[k_+d_+1] - knots_[k_+1]) right = 0.0 if (denom_r != 0.0): right = ((knots_[k_+d_+1] - t_) / denom_r) * cox_de_boor(d_-1, t_, k_+1, knots_) return left + right def interpolate( d_, P_, n_, ts_, knots_ ): A = np.zeros((n_, n_)) for i in range(n_): for j in range(n_): A[i, j] = cox_de_boor(d_, ts_, j, knots_) control_points = np.linalg.solve(A, P_) return control_points def create_B_spline( d_, P_, t_, knots_): sum = Vector() # just a vector class. for i in range( len(P_) ): sum += P_ * cox_de_boor(d_, t_, i, knots_) return sum def B_spline( points_ ): d = 3 # change to 2 for quadratic. P = np.array( points_ ) n = len( P ) ts = np.linspace( 0.0, 1.0, n ) knots = knotvector.generate( d, n ) # len = n + d + 1 control_points = interpolate( d, P, n, ts, knots) crv_pnts = [] for i in range(10): t = float(i) / 9 crv_pnts.append( create_B_spline(d, control_points, t, knots) ) return crv_pnts control_points = [ [float(i), math.sin(i), 0.0] for i in range(4) ] cps = B_spline( control_points ) Result is OK when interpolating 4 points (control vertices): Result is NOT OK when interpolating 5 points (control vertices): Result is OK when interpolating 6 points (control vertices): and so on... I noticed two things: The spline does not interpolate properly when the number of control vertices is odd. The spline interpolates properly with any number of vertices when the degree becomes quadratic. So, if you change d = 2, in the B_spline function, the curve will interpolate properly for odd and even number of control vertices. The cox de boor function is correct and according to the mathematical expression, but with a small alteration on the 2nd conditional expression t <= t **<=** t[i+1] (see my previous SO question here for more details). Also, I used numpy to solve the linear system, which also works as expected. Other than np.linalg.solve, I have tried np.linalg.lstsq but it returns the same results. I honestly do not know where to attribute this abnormal behaviour. What could cause this issue? Continue reading...