Piecewise cubic Bézier curve interpolation

bezier_cubic_splines.m

clear;

f = @(x) cos(x);
xa = 0;
xb = 2 * pi;
n = 4; % number of cubic polynomials

x = linspace(xa, xb, n + 1); % x_0, ..., x_n
y = f(x); % y_0, ..., y_n

h = zeros(n, 1); % h_0, ..., h_{n - 1}
d = zeros(n, 1); % d_0, ..., d_{n - 1}
for k = 0:n - 1
    h(k + 1) = x(k + 2) - x(k + 1);
    d(k + 1) = (y(k + 2) - y(k + 1)) / h(k + 1);
end

a = h(2:n) / 6; % a_1, ..., a_{n - 1}
b = zeros(n - 1, 1); % b_1, ..., b_{n - 1}
u = zeros(n - 1, 1); % u_1, ..., u_{n - 1}
for i = 1:n - 1
    b(i) = (x(i + 2) - x(i)) / 3;
    u(i) = d(i + 1) - d(i);
end

C = diag(b) + diag(a(1:n - 2), 1) + diag(a(1:n - 2), -1);
m = [0; C\u; 0]; % m_0, ..., m_n (natural spline)

% % Old process: (gives less accurate results)
% NC = zeros(n + 1);
% for k = 2:n
%     NC(k, k) = 2 * (h(k - 1) + h(k));
%     NC(k, k - 1) = h(k - 1);
%     NC(k, k + 1) = h(k);
% end
% % Natural spline conditions
% NC(1, 1) = 1;
% NC(n + 1, n + 1) = 1;
% 
% v = zeros(n + 1, 1);
% for i = 2:n
%     %v(i) = 3 * u(i - 1);
%     v(i) = 3 * ((y(i + 1) - y(i)) / h(i) - (y(i) - y(i - 1)) / h(i - 1));
% end
% m = [0; NC\v; 0]; % m_0, ..., m_n (natural spline)

% Find polynomial coefficients in monomial basis
S = zeros(4, n); % s_{k, 0}, ..., s_{k, 3} for k-th polynomial (k = 0, ..., n - 1)
for k = 0:n - 1
    S(1, k + 1) = y(k + 1);
    S(2, k + 1) = d(k + 1) - (h(k + 1) * (2 * m(k + 1) + m(k + 2))) / 6;
    S(3, k + 1) = m(k + 1) / 2;
    S(4, k + 1) = (m(k + 2) - m(k + 1)) / (6 * h(k + 1));
end

% Spline polynomial basis (1, ..., (x - x_k)^3)
syms t u xk hk;
p = sym('Bas', [4, 1]);
for i = 0:3
    p(i + 1) = (t - xk)^i;
end

% Plot polynomials and compute rough approximation error
fplot(f, [xa, xb]);
hold on;
% err = 0;
for k = 0:n - 1
    interval = [x(k + 1), x(k + 2)];
    pk = S(:, k + 1).' * subs(p, xk, x(k + 1));
    fplot(@(x) double(subs(pk, t, x)), [x(k + 1), x(k + 2)]);
%     err = err + integral(@(x) norm(f(x) - double(subs(pk, t, x))), ...
%                             x(k + 1), x(k + 2), 'ArrayValued', true);
end
hold off;

% Compute Bernstein polynomial basis of degree 3 that, when
% used to express the Bézier curve given by the k-th found
% cubic polynomial, the parameters t=0 and t=1 map the
% x-coordinate to x = x_{k - 1} and x = x_k, respectively.
global pb;
pb = sym('Ber', [4, 1]);
u = (t - xk) / hk;
for i = 0:3
    pb(i + 1) = nchoosek(3, i) * u^i * (1 - u)^(3 - i);
end

Sb = zeros(4, n);
% The x-coordinate of all Bézier curves' control points,
% 4 per polynomial.
cx = zeros(n * 4, 1);
for k = 0:n - 1
    % Make change of basis to Bernstein polynomial basis:
    % Let `basis` be denoted by u and let b be the Bernstein polynomial
    % basis of degree 3. We make use of the following identity:
    % M_{u, b} = M_{e, b} M_{u, e}, where e is the monomial basis. Thus,
    % M_{u, b} = M_{b, e}^-1 M_{u, e}.
    basis = subs(p, xk, x(k + 1));
    Mbe = tobasis(subs(pb, [xk, hk], [x(k + 1), h(k + 1)]));
    Meb = inv(Mbe);
    Mue = tobasis(basis);
    Mub = Meb * Mue;
    Sb(:, k + 1) = Mub * S(:, k + 1);
    
    % Find the x-coordinate of the control points for the
    % Bézier curve representing this cubic polynomial.
    % The y-coordinate of the k-th control point is equal to
    % the {k - 1}-th coefficient of the polynomial in its
    % Bernstein basis.
    cx(4 * k + 1:4 * k + 4) = Meb * [0; 1; 0; 0]; % linear on x
end

cy = Sb(:);
[cx cy]

% To print the expression of the k-th polynomial in
% its Bernstein basis, run:
% Sb(:, k).' * subs(pb, [xk, hk], [x(k + 1), h(k + 1)])

function Mbe = tobasis(basis)
    % Builds a change-of-basis matrix M_{basis, e}, where
    % e is the monomial basis and `basis` is in P^4.
    Mbe = zeros(4);
    for i = 1:length(basis)
        % sym2poly returns coefficients of larger terms first
        coeffs = flip(sym2poly(basis(i)));
        Mbe(1:length(coeffs), i) = coeffs;
    end
end