% B-spline Tutorial
% by Ashish Myles (marcianx@gmail.com)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

knots = 1:14;
cpts = abs([1:10]-5);
cpts = [0 1 2 3 4 0 1 2 3 4];
cpts = [3 3 3 3 3 0 0 0 0 0];

nk = length(knots);
nc = length(cpts);
deg = nk - nc - 1;

grev = greville(knots, deg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Interpolating samples
    sp = spapi(knots, grev, cpts);

    figure(1);
    plot_spline(sp);
    plot_knots_and_cpts(knots, cpts, 1); % Display interpolated points

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Using samples as control points
    sp = spmak(knots, cpts);

    figure(2);
    plot_spline(sp);
    plot_knots_and_cpts(knots, cpts, 1);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Variation diminution
    cpts = [3 4 0 3 0 2 -1 1 0 2];
    sp = spmak(knots, cpts);

    figure(3);
    plot_spline(sp);
    plot_knots_and_cpts(knots, cpts, 3);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Finite support
    cpts = [0 0 0 0 1 -1 0 0 0 0];
    sp = spmak(knots, cpts);

    figure(4);
    plot_spline(sp);
    plot_knots_and_cpts(knots, cpts, 0);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Periodic
    cpts = [0 1 0 1 0 1 0 1 0 1];
    sp = spmak(knots, cpts);

    figure(5);
    plot_spline(sp);
    plot_knots_and_cpts(knots, cpts, 0);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 1

    % Basis function
    % Uniform
    for d = 2:2
        cpts = zeros(1,2*d+3);
        cpts(d+2) = 1;
        knots = 0:(length(cpts)+d);
        knots(4)=3.5;
        sp = spmak(knots, cpts);

        figure(d+1);
        plot_spline(sp);
        plot_knots_and_cpts(knots, cpts, 4);
        title(sprintf('knots: [ %s] (uniform)', sprintf('%g ', knots)));
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Curve as linear combination of uniform Basis functions
    c = ['b', 'g', 'm', 'c'];
    for d = 2:3
        cpts = [1,3,2,-1];
        nc = length(cpts);
        knots = 1:(d+nc+1);
        sp = spmak(knots, cpts);
        figure(d);
        plot_spline(sp, 2);
        if 1 % display bases?
            for i = 1:nc
                bas_cpts = zeros(1,nc);
                bas_cpts(i) = cpts(i);
                sp = spmak(knots, bas_cpts);
                plot_spline(sp, 1, [c(mod(i-1,length(c))+1)]);
            end
        end
        plot_knots_and_cpts(knots, cpts, 7);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Non-uniform
    c = ['b', 'g', 'm', 'c'];
    d = 3;
    cpts = [1,3,2,-1,2,1,2,1,2];
    nc = length(cpts);
    knots = 1:(d+nc+1);
    knots(5) = 5.8;
    knots(6) = 6.0;
    knots(7) = 6.2;
    knots(8) = 6.4;
    knots(9) = 6.6;
    sp = spmak(knots, cpts);
    for i = 1:2
        figure(i);
        plot_spline(sp, 2);
        plot_knots_and_cpts(knots, cpts, 7);
        if i == 2
            for i = 1:nc
                bas_cpts = zeros(1,nc);
                bas_cpts(i) = cpts(i);
                sp = spmak(knots, bas_cpts);
                plot_spline(sp, 1, [c(mod(i-1,length(c))+1)]);
            end
        end
        plot_knots_and_cpts(knots, cpts, 7);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Another non-uniform
    c = ['b', 'g', 'm', 'c'];
    d = 3;
    cpts = [1,1,1,1,-1,-1,-1,-1];
    nc = length(cpts);
    knots = 1:(d+nc+1);
    knots(5) = 4.70;
    knots(end-2:end) = [4.90, 5.10, 5.30];
%   knots(5) = 5;
%   knots(end-2:end) = [5, 5, 5];
    knots = sort(knots);
    sp = spmak(knots, cpts);
    for i = 1:2
        figure(i);
        plot_spline(sp, 2);
        if i == 2
            for i = 1:nc
                bas_cpts = zeros(1,nc);
                bas_cpts(i) = cpts(i);
                sp = spmak(knots, bas_cpts);
                plot_spline(sp, 1, [c(mod(i-1,length(c))+1)]);
            end
        end
        plot_knots_and_cpts(knots, cpts, 5);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Evaluation by knot insertion
    d = 3;
    cpts = [1,2,2,1];
    nc = length(cpts);
    knots = 1:(d+nc+1);
    kn = 4.5;
    sp = spmak(knots, cpts);
    figure(1);
    plot_spline(sp, 1);
    plot_knots_and_cpts(knots, cpts, 7);
    plot_range_x = [2.5, 6.5];
    plot_range_y = [1.7, 2.1];
    axis([ plot_range_x, plot_range_y ]);
    for i = 1:d
        figure(i+1);
        sp = fnrfn(sp, kn);
        plot_spline(sp, 1);
        plot_knots_and_cpts(sp.knots, sp.coefs, 7);
        plot_knots_and_cpts(knots, cpts, 0);
        axis([ plot_range_x, plot_range_y ]);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Simple spline example
    d = 3;
    cpts = [0,2,3,1,0,2];
    nc = length(cpts);

    % Uniform
    knots = 1:(d+nc+1);
    sp = spmak(knots, cpts);
    figure(1);
    plot_spline(sp,2);
    plot_knots_and_cpts(knots, cpts, 3);
    xlim([3,8]);

    % End point interpolating
    knots = [ones(1,d), 1:(nc-d+1), ones(1,d)*(nc-d+1)];
    sp = spmak(knots, cpts);
    figure(2);
    plot_spline(sp,2);
    plot_knots_and_cpts(knots, cpts, 3);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Knot multiplicity
    d = 3;
    cpts = [1,1,1,2,2,2,1,2,1,1,1];
    nc = length(cpts);
    for i = 1:d+1
        knots = sort([ 1:(d+nc+1-(i-1)), ones(1,i-1)*7 ]);
        sp = spmak(knots, cpts);
        figure(i);
        plot_spline(sp);
        plot_knots_and_cpts(sp.knots, sp.coefs, 5);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Uniform knot insertion
    d = 3;
    cpts = [1,2,2,1];
    nc = length(cpts);
    knots = 1:(d+nc+1);
    sp = spmak(knots, cpts);
    figure(1);
    plot_spline(sp, 1);
    plot_knots_and_cpts(knots, cpts, 7);
    plot_range_x = [2.5, 6.5];
    plot_range_y = [1.7, 2.1];
    axis([ plot_range_x, plot_range_y ]);
    for i = 1:d
        figure(i+1);
        sp = fnrfn(sp, 0.5 * (sp.knots(d:end-d) + sp.knots(d+1:end-d+1)));
        plot_spline(sp, 1);
        plot_knots_and_cpts(sp.knots, sp.coefs, 7);
        plot_knots_and_cpts(knots, cpts, 0);
        axis([ plot_range_x, plot_range_y ]);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Uniform knot insertion
    d = 3;
    cpts = [1,2,2,1];
    nc = length(cpts);
    knots = 1:(d+nc+1);
    sp = spmak(knots, cpts);
    figure(1);
    plot_spline(sp, 1);
    plot_knots_and_cpts(knots, cpts, 7);
    plot_range_x = [2.5, 6.5];
    plot_range_y = [1.7, 2.1];
    axis([ plot_range_x, plot_range_y ]);
    for i = 1:d
        figure(i+1);
        sp = fnrfn(sp, 0.5 * (sp.knots(d:end-d) + sp.knots(d+1:end-d+1)));
        plot_spline(sp, 1);
        plot_knots_and_cpts(sp.knots, sp.coefs, 7);
        plot_knots_and_cpts(knots, cpts, 0);
        axis([ plot_range_x, plot_range_y ]);
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0

    % Differentiation
    d = 3;
    cpts = [1,1,1,2,2,2,1,2,1,1,1];
    nc = length(cpts);
    knots = 1:(d+nc+1);
    sp = spmak(knots, cpts);

    figure(1);
    plot_spline(sp);
    plot_knots_and_cpts(knots, cpts, 7);

    dsp = fnder(sp);
    dsp = spmak(dsp.knots(2:end-1), dsp.coefs(2:end-1));
    figure(2);
    plot_spline(dsp);
    plot_knots_and_cpts(dsp.knots, dsp.coefs, 7);
    plot_knots_and_cpts(knots, dsp.coefs, 0); % Set the limits appropriately

    isp = fnint(sp);
%   figure(3);
%   plot_spline(isp);
%   plot_knots_and_cpts(isp.knots, isp.coefs, 7);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%