2D B-Spline approximation in Matlab

I am dealing with interpolation problems currently. I read about B-Spline approximation. So I tried to implement a matlab script for a better understanding of the B-Spline's mathematics.

In the first case I used my script to approximate a trapezoid B-Splines, which worked very well using four control points and the degree m = 3.

After that I added two more control points to check if the script runs correctly. But the approximation looks weird.

I am using the mahematics from: http://www.tm-mathe.de/Themen/html/funbezierbspline.html

So I think there is a mistake in my script. But I could not figure out where it is.

Main script:

%% Console
close all;
clear all;
format long;
clc;

%% Parameters

% Contron points
C = [0.0 1.0 4.0 5.0 6.0 7.0;  % x components
     0.0 1.0 1.0 0.5 0.5 0.0]; % y components

 C = [0.0 1.0 4.0 5.0 ;  % x components
      0.0 1.0 1.0 0.0]; % y components

% B-Spline Degree
m = 3;

% B-Spline
s = f_Bspline(C, m, 0.001);

%% Plot

figure(1);
plot(C(1,:), C(2,:),'o');
hold on;
plot(C(1,:), C(2,:),'--');
hold on;
plot(s(1,:), s(2,:));
grid on;
grid minor;
xlabel('x');
ylabel('y');
legend({'Control points','Polygon','B-Spline-Approx.'}, 'Location', 'southeast');
ylim([min(C(2,:)) max(C(2,:))*1.25]);

f_BSpline:

function s = f_Bspline(C, m, step)
% Calculates a B-Spline of degree m using the control points C.
% 
% C: 2-dimensional control points (x_0, ... , x_n; y_0, ... , y_n)
% m: B-Spline degree
% step: Step size of t.
%
% s: B-Spline. s(1,:) -> x component, s(2,:) -> y component

%% Parameters

% Control point's x and y components
x = C(1,:);
y = C(2,:);

% Number of control point - 1
n = size(x,2) - 1;

% Knot vector
T = f_BSpline_KnotVector(m,n);

% B-Spline intervall
t = 0:step:(n-m+1);

%% Calculate B-Spline

for z=1:1:size(t,2)
    ti = t(z);
    s(1,z) = 0;
    s(2,z) = 0;
    for i=0:1:n
        % Base B-Spline
        B = f_BSpline_BaseSpline(i, m, ti, T);
        % x component
        s(1,z) = s(1,z) + x(i+1) * B;
        % y component
        s(2,z) = s(2,z) + y(i+1) * B;
    end
end

end

f_BSpline_KnotVector:

function T = f_BSpline_KnotVector(m,n)
% Calculate knot sequence.
% m: Degree of B-Spline
% n: Number of control points - 1
%
% T = [t0, ... , t_n+m+1] : Knot sequence / vector

T = zeros(1, (n+m+2));
for j=0:1:(n+m+1)
    if j <= m
        Ti = 0;
    elseif j >= (m+1) && j <= n
        Ti = j - m;
    elseif j > n
        Ti = n - m + 1;
    end
    T(j+1) = Ti;
end

end

f_Bspline_BaseSpline:

function B = f_BSpline_BaseSpline(i, k, t, T)
% Calculate Base B-Spline B_i,k
% k: B-Spline degree
% T: Knot sequence
% t: Current t parameter
%
% B: Base B-Spline at t.

% Index shift
j = i + 1;

if k == 0
    % End of recusrion
    if t >= T(j) && t <= T(j+1)
        B = 1;
    else
        B = 0;
    end
else
    % Check dividing by zero
    if T(j+k) ~= T(j)
        A = (t -T(j))/(T(j+k) - T(j));
    else
        A = 0;
    end
    if T(j+k+1) ~= T(j+1)
        B = (T(j+k+1) - t) / (T(j+k+1) - T(j+1));
    else
        B = 0;
    end
    % Calculate base B-Spline
    B1 = f_BSpline_BaseSpline(i,   k-1, t, T);
    B2 = f_BSpline_BaseSpline(i+1, k-1, t, T);
    B = A * B1 + B * B2;
end


end

The problem is with your evaluation of degree 0 basis functions at the internal knots. As there are no internal knots in your first example, the problem only becomes visible in the second example.

In f_BSpline_BaseSpline you write

    % End of recusrion
    if t >= T(j) && t <= T(j+1)
        B = 1;
    else
        B = 0;
    end

However, the website you linked suggests comparing with < instead of <=:

    % End of recursion
    if t >= T(j) && t < T(j+1)
        B = 1;
    else
        B = 0;
    end

This fixes the spikes The problem was that at the interior knots your basis functions summed to more than one; note how the spikes point away from the origin.

However, if you look closer at the horizontal axis, you can notice another problem: the curve is now closed! The last degree 0 basis function does have to have a closed support, i.e., in that particular case there should be the <= that we have removed from your code. It is certainly a pity that the website did not explain that. So this is what the final result looks like:

And here is the source code of the corrected f_BSpline_BaseSpline

function B = f_BSpline_BaseSpline(i, k, t, T)
% Calculate Base B-Spline B_i,k
% k: B-Spline degree
% T: Knot sequence
% t: Current t parameter
%
% B: Base B-Spline at t.

% Index shift
j = i + 1;

if k == 0
    % Corrected end of recursion
    if (t >= T(j) && t < T(j+1))
        B = 1;
    % Exception: the last basis function is equal to 1 also at the last knot.
    % There might be a more elegant way of writing it; I am no Matlab expert.
    elseif ( T(j+1) == T(end) && t == T(end))
        B = 1;
    else
        B = 0;
    endif
else
    % Check dividing by zero
    if T(j+k) ~= T(j)
        A = (t -T(j))/(T(j+k) - T(j));
    else
        A = 0;
    end
    if T(j+k+1) ~= T(j+1)
        B = (T(j+k+1) - t) / (T(j+k+1) - T(j+1));
    else
        B = 0;
    end
    % Calculate base B-Spline
    B1 = f_BSpline_BaseSpline(i,   k-1, t, T);
    B2 = f_BSpline_BaseSpline(i+1, k-1, t, T);
    B = A * B1 + B * B2;
end


end