function figure_handle = show_shape_model(shape_model, ref_points_vec, ref_images_vec, n_sds, image_width)
% SHOW_SHAPE_MODEL: Shows the shape model for each one of the modes.
%
% Code written by Katherine Smith, 2003
%
%    GENERAL
%
%      show_shape_model(shape_model, ref_points_vec, ref_images_vec, n_sds)
%
%    INPUT/S
%
%      -shape_model:
%           The shape model to be shown.
%
%      -ref_points_vec:
%           Vector of reference points.
%
%      -ref_images_vec: 
%           Vector of reference images.
%
%      -n_sds:
%           Number of standard deviations to be shown?
%
%           
%    OUTPUT/S
%
%      -None.
%
%    PENDING WORK
%
%      -
%
%    KNOWN BUG/S
%
%      -None.
%
%    COMMENT/S
%
%      -Number of standard variation input seems to contol the number of
%       figures shown. For 3, it will generate it for -3,-2,-1,-,1,2,3.
%
%    RELATED FUNCTION/S
%
%      SHOW_INTENSITY_MODEL, SHOW_COMBINED_MODEL
%
%    ABOUT
%
%      -Created:     November 23rd, 2003
%      -Last update: Novermber 26th, 2003
%      -Revision:    0.0.2
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================

n_modes = size(shape_model.pcs,2); 
                % get the number of modes
fig = figure('Name',(['Shape Model']));
figure_handle = figure(fig);

for m=1:n_modes              
  for i=-n_sds:n_sds
    % Smith: what are the parameters corresponding to this sd
    % create example with these parameters
    % RSS: it seems like it shows in a figure the deviations from either side
    % in the range of parameters
    example_params = repmat(i*shape_model.sd(1,m),size(shape_model.pcs,1),1);
            % get the parameters by extracting the mode of variation from
            % the model.
    example_points = example_params.*shape_model.pcs(:,m) + shape_model.mean';
            % get the points with the modes of variation applies to the mean
    example = interp1(ref_points_vec, ref_images_vec, example_points);
            % interpolate in 1-D
    figure(fig);
    % subplot(1,n_sds*2+1,i+n_sds+1);
    subplot( n_modes, n_sds * 2 + 1, (m - 1) * (n_sds * 2 + 1) + i + n_sds + 1);
    plot(example);
    axis([0 image_width -0.3 1]);
    if (m == 1 & i ==1),
      text(1,1,'Shape Model');
    end
            % plot figures for that mode. For m modes of varation and n standard
            % deviations, (2*n+1) * m figures will be shown in total??
  end
end