% -------------------------------------------------------------------------
% define a cylinder using Bezier curves for the radius function
% and draw the surface of revolution superimposed on the curve
% and the polygon of control points
% -------------------------------------------------------------------------

% setup the plotting parameters
feature('UseGenericOpenGL', 1);
set(gcf, 'Renderer', 'OpenGL');

bcpn = [ % Bezier control polygon for the nose section
  25  19  17  11
   0   2   1   1
   0   0   0   0
   1   1   1   1
];
xmax = max(bcpn(1,:));

bcpt = [ % Bezier control polygon for the tail section
  11   5   4   2   0
   1   1 1.5 1.5   1
   0   0   0   0   0
   1   1   1   1   1
];
xmin = min(bcpt(1,:));

nti = 21; % number of values of t at which the Bezier curve gets evaluated
nxi = 101; % number of gridlines along the length of the cylinder
nyi = 51; % number of gridlines around the circumference of the cylinder
dx = xmax - xmin; % extent of the curve along the x axis
xi = [xmin:dx/(nxi-1):xmax]; % nxi evenly spaced values from xmin to xmax
nxy = nxi*nyi; % total number of grid points

% define the two Bezier curves and concatenate them to make a single curve
P1 = BezierCurve(bcpn, nti); P2 = BezierCurve(bcpt, nti);
r = [P1(1:nti,1:2); P2(2:nti,1:2)]';
% since the curves are evaluated at evenly spaced intervals in 0 <= t <= 1,
% and the cylinder function expects evenly spaced intervals in 0 <= x <= 1,
% we need to resample the curve and rescale by 1/dx.
r = interp1(r(1,:), r(2,:), xi)/dx;
% define the surface of revolution using the cylinder function
[x, y, z] = cylinder(r, nyi-1);
% we now have a surface, ruled with nxi x nyi grid lines, and embedded
% in space such that the grid point p(i,j) = [x(i,j) y(i,j) z(i,j)]'
% so x, y and z are nxi x nyi matrices of coordinates at the grid points
% now, to use our transformation matrices, we need a matrix of vertices,
% so we must make such a matrix out of the x, y and z matrices

% -------------------------------------------------------------------------
% one way to do this is by direct manipulation of the matrices

% vb = [
%   reshape(x, 1, nxy);
%   reshape(y, 1, nxy);
%   reshape(z, 1, nxy);
%   ones(1, nxy)
% ];
% % the cylinder rotates the radius function about the z axis, and our
% % Bezier curve is given by y as a function of x, so we rotate the surface
% % of revolution about y to align its axis of symmetry with the x axis
% R = ConstructRotation('y', pi/2); vb1 = R * vb;
% % recover the matrices of coordinate values in preparation for drawing
% x = reshape(vb1(1,:), nxi, nyi);
% y = reshape(vb1(2,:), nxi, nyi);
% z = reshape(vb1(3,:), nxi, nyi);
% 
% newplot; hold on % display the results
% 
% surf(x, y, z, ... draw the surface of revolution sem-transparent
%   'FaceAlpha', 0.4, ...
%   'FaceLighting', 'phong', 'FaceColor', [0.9 1.0 1.0], ...
%   'EdgeLighting', 'phong', 'EdgeColor', 'none', ...
%   'AmbientStrength', 0.65, 'DiffuseStrength', 0.4, ...
%   'SpecularStrength', 1.0, 'SpecularExponent', 7.0, ...
%   'SpecularColorReflectance', 1.0 ...
% );

% -------------------------------------------------------------------------
% another way to do this is with the surf2patch function

% with surf2patch, we create a structure of the form
%   fvc = struct( ...
%     'faces',[(nxi-1)*(nyi-1) x 4, double], ...
%     'vertices',[(nxi*nyi) x 3, double] ...
%   )
% which can be passed to the patch function for drawing
fvc = surf2patch(x, y, z);
% the cylinder rotates the radius function about the z axis, and our
% Bezier curve is given by y as a function of x, so we rotate the surface
% of revolution about y to align its axis of symmetry with the x axis
R = ConstructRotation('y', pi/2);
% since fvc.vertices stores vertices in rows, we postmultiply by R'
vb = [fvc.vertices ones(size(fvc.vertices,1),1)] * R';
fvc.vertices = vb(:,1:3);

newplot; hold on % display the results

patch(fvc, ... draw the surface of revolution sem-transparent
  'FaceAlpha', 0.4, ...
  'FaceLighting', 'phong', 'FaceColor', [0.9 1.0 1.0], ...
  'EdgeLighting', 'phong', 'EdgeColor', 'none', ...
  'AmbientStrength', 0.65, 'DiffuseStrength', 0.4, ...
  'SpecularStrength', 1.0, 'SpecularExponent', 7.0, ...
  'SpecularColorReflectance', 1.0 ...
);

% -------------------------------------------------------------------------

% overlay the scaled Bezier curve and polygon of control points
PlotFig3(bcpn/dx, [1 1 1 0], [0 0 1]);
PlotFig3(bcpt/dx, [1 1 1 1 0], [0 0 1]);
plot(xi/dx, r, 'color', [1 0 0]);

% lighting and viewing parameters
light( ...
  'Position', [1 0 .5], 'Style', 'infinite', ...
  'Color', [1.0 0.5 0.0] ...
);
view(0, 90); axis([-.1 1.1 -.6 .6 -.6 .6], 'square');
grid on; hold off;
