clc;
close all;
clear all;

% Control points
P = [0,0; 0,2; 2,2; 2,0; 4,0; 5.5,2; 6,2; 6,0];

% Cubic Bernstein Polynomials
t = 0:.05:1;
Bezier = [(1-t).^3; 3*t.*(1-t).^2; 3*t.^2.*(1-t); t.^3];

% Bézier extraction matrices (deduced using Blossoming)
% In order to keep the matrices nice and clear, define some constants
a = 1/6;
b = 2/6; % = 1/3
c = 4/6; % = 2/3
d = 7/12;
e = 1/4;
f = 1/2;

figure()

% State here which subfigure you'd like to plot.
for Subfigure = ['a', 'b', 'c', 'd']
    
    switch Subfigure
        case 'a'
            % This is the matrix for Figure 7.19a
            C = [a c a 0 0 0 0 0; ...
                0 c b 0 0 0 0 0; ...
                0 b c 0 0 0 0 0; ...
                0 a c a 0 0 0 0; ...
                0 0 c b 0 0 0 0; ...
                0 0 b c 0 0 0 0; ...
                0 0 a c a 0 0 0; ...
                0 0 0 c b 0 0 0; ...
                0 0 0 b c 0 0 0; ...
                0 0 0 a c a 0 0; ...
                0 0 0 0 c b 0 0; ...
                0 0 0 0 b c 0 0; ...
                0 0 0 0 a c a 0; ...
                0 0 0 0 0 c b 0; ...
                0 0 0 0 0 b c 0; ...
                0 0 0 0 0 a c a];
            
        case 'b'
            % This is the matrix for Figure 7.19b
            C = [a d e 0 0 0 0 0; ...
                0 f f 0 0 0 0 0; ...
                0 0 1 0 0 0 0 0; ...
                0 0 f f 0 0 0 0; ...
                0 0 0 1 0 0 0 0; ...
                0 0 0 f f 0 0 0; ...
                0 0 0 e d a 0 0; ...
                0 0 0 0 c b 0 0; ...
                0 0 0 0 b c 0 0; ...
                0 0 0 0 a c a 0; ...
                0 0 0 0 0 c b 0; ...
                0 0 0 0 0 b c 0; ...
                0 0 0 0 0 a c a];
            
        case 'c'
            % This is the matrix for Figure 7.19c
            C = [a d e 0 0 0 0 0; ...
                0 f f 0 0 0 0 0; ...
                0 0 1 0 0 0 0 0; ...
                0 0 0 1 0 0 0 0; ...
                0 0 0 0 1 0 0 0; ...
                0 0 0 0 f f 0 0; ...
                0 0 0 0 e d a 0; ...
                0 0 0 0 0 c b 0; ...
                0 0 0 0 0 b c 0; ...
                0 0 0 0 0 a c a];
            
        case 'd'
            % This is the matrix for Figure 7.19d
            C = [a d e 0 0 0 0 0; ...
                0 f f 0 0 0 0 0; ...
                0 0 1 0 0 0 0 0; ...
                0 0 0 1 0 0 0 0; ...
                0 0 0 0 1 0 0 0; ...
                0 0 0 0 0 1 0 0; ...
                0 0 0 0 0 f f 0; ...
                0 0 0 0 0 e d a];
    end
    
    % double('a') = 97
    subplot(2,2, double(Subfigure-96) )    
    
    % Plot control polygon and Bézier points
    plot(P(:,1), P(:,2), 'ko--', 'LineWidth', 1)
    axis([-.5, 6.5, -1, 3])
    axis equal
    hold on
    Pnew = C*P;
    plot(Pnew(:,1), Pnew(:,2), 'b*')
    
    % for Figure 7.19a use: for i=[1 4 7 10 13]
    % for Figure 7.19b use: for i=[1 4 7 10]
    % for Figure 7.19c use: for i=[1 4 7]
    % for Figure 7.19d use: for i=[1 5]
    
    if Subfigure == 'd'
        Loop = [1 5];
    else
        Loop = [1:3:size(C,1)-1];
    end
    
    for i = Loop
        Curve11 = Pnew(i,1)*Bezier(1,:) + Pnew(i+1,1)*Bezier(2,:) + Pnew(i+2,1)*Bezier(3,:) + Pnew(i+3,1)*Bezier(4,:);
        Curve12 = Pnew(i,2)*Bezier(1,:) + Pnew(i+1,2)*Bezier(2,:) + Pnew(i+2,2)*Bezier(3,:) + Pnew(i+3,2)*Bezier(4,:);
        
        plot(Curve11,Curve12, '-', 'Color', [1 .7 0], 'LineWidth', 2)
        title(['Figure 7.19' Subfigure])
    end
    
end