function [status] = build_1d_bump_model( ...
    spline_type ...
    , placement_type ...
    , objective_function ...
    , apply_intensity_offset ...
    , save_movie ...
    , do_plot ...
    , movie_type ...
    , progress_bar ...
    , user_menu ...
    , produce_statistics ...
    , see_model_before_each_iteration ...
    , show_pixels ...
    , show_score_vs_iteration ...
    , show_model_after_optimisation ...
    , show_specificity_statistics ...
    , verbose_score ...
    , movie_name ...
    , window ...
    , average_smooth_enabled ...
    , gaussian_smooth_enabled ...
    , bump_width ...
    , bump_width_variation ...
    , bump_height ...
    , bump_height_variation ...
    , bump_position_freedom ...
    , initial_diminish_factor ...
    , smoothness_factor ...
    , plot_before_and_after ...
    , number_of_bins ...
    , model_evaluation_method ...
    , show_score_using_model ...
    , handle ...
    , n_modes ...
    , shape_weight ...
    , weighting_normalisation_method ...
    , variation_kept ...
    , run_all_objective_functions ...
    , kernel ...
    , save_data ...
    , load_data ...
    , save_bumps ...
    , draw_warp_curve ...
    , overlap_scores ...
    , show_target ...
    , format ...
    , filter ...
    , filter_level ...
    , model_score_type ...
    , gen_iters ...
    , spec_iters ...
    , draw_curve_by_images ...
    , show_mesh ...
    , mesh_type ...
    , exit_at_end ...
    , show_all_score_figure_types ...
    , close_at_end ...
    , frames_per_second ...
    , show_score_bar ...
    , n_points ...
    , show_warp_targets ...
    , show_registration_target ...
    , produce_record ...
    , initialise_warps_close_to_target ...
    , save_in_log ...
    , perturbation_method ...
    , offset_extent ...
    , force_reference ...
    , cycle ...
    , verbose ...
    , save_stats ...
    , n_modes_display ...
    , n_standard_deviations_display ...
    , show_warps_while_optimising ...
    , description ...
    , find_perfect_warp ...
    , force_first_reference ...
    , display_shape_model ...
    , display_intensity_model ...
    , plot_alignment ...
    , error_verbosity ...
    , precision_required ...
    , background_red ...
    , background_green ...
    , background_blue ...
    , dynamic_precision ...
    , precision_automation_type ...
    , rotate_surface ...
    , draw_final_warps ...
    , distance_type ...
    , retain_original_peak ...
    , peak_type ...
    , add_points ...
    , add_points_cycle ...  
    , timestamp_images ...
    , auto_weighting ...
    , automatic_weighting_type ...
    , show_inverse_projection ...
    , show_model_reconstruction_for_warped_images ...
    , point_insertion_position ...
   )

% BUILD_1D_BUMP_MODEL: Builds a 1-D model for some bump data
%                      and returns some related statistics.
%
% Code originally written by Katherine Smith, 2003
%
%    GENERAL
%
%      Omitted. Much of the explanation is included in the section where
%      default arguments are being assigned.
%
%    INPUT/S
%
%      -- No arguments will lead to defaults
%
%      -spline_type: 
%           The type of spline to be used. 'single_point' or
%           'multi_point' at present.
%
%      -placement_type:
%           Placement type for 'multi_point':
%                    'grid' or 'edges' or 'random'
%           placement type for 'single_point':
%                    'overlap_grid' or 'edges_and_scale'
%                    or 'random_and_scale'.
%
%      -objective_function:
%           'msd_opt_together' or 'model_opt_separare' or
%           'msd_opt_separate' or 'model_opt_together'.
%
%      .....
%           
%    OUTPUT/S
%
%      -status:
%           The status of the operation of this routine (to be returned at
%           the end
%
%      -Figures produced and results displayed are main
%       data to be utilised.
%
%    PENDING WORK
%
%      -
%
%    KNOWN BUG/S
%
%      -
%
%    COMMENT/S
%
%      -Different video outputs are generated but controlling which
%       one is chosen is not trivial. It currently just does the
%       job of experimentation and showing results as pixels
%       or bumps.
%
%    RELATED FUNCTION/S
%
%      See other functions in same directory, filenames
%      beginning with "build" in particular. There seem to be much
%      overlap due to cut-and-paste.
%
%    ABOUT
%
%      -Created:     December 29rd, 2003
%      -Last update: March 30th, 2004
%      -Revision:    0.9.2
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% variables (parameters) initialisation           %%
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input Defaults              %%   

% When the function is not called from the GUI, various arguments must be
% initialised so that the function performs a registration task.

% A possible extension would allow for _partial_ arguments to be given so
% that a script can run different experiments and then logs are necessary
% too.

if (nargin == 0),
    disp(' ');
    disp('Argument(s) omitted. Using defaults instead.');
    disp(' ');
         % notify user about default values being used.
    spline_type = 'single_point';
         % SMITH: spline type 'single_point' or 'multi_point'
         % SMITH: placement type for multipoint: 'grid' or 'edges' or 'random'
         % SMITH: placement type for singlepoint: 'overlap_grid' or 'edges_and_scal or 'random_and_scale'
    placement_type = 'random_and_scale';
         % controls knot-point selection on the grid
    objective_function = 'model_opt_together';
         % choice of objective function
    apply_intensity_offset = 1;
         % 0: Offset should not be applied to data
         % 1: Offset enabled
    exit_at_end = 0;
         % indicates whether AART quits at the end of registration
    description = '';
         % a name provided for describing registration
    cycle = 0;
          % cycle size for user display of objective function evaluation
    produce_record = 0;
          % indicates if results should be displayed to the user in the
          % console
    save_in_log = 0;
          % indicates if an HTML file or network of files will be used to
          % record the results
    close_at_end = 0;
          % indicates if the GUI should be closed when registration is done
    save_movie = 0;
          % a flag for movie output
    force_reference = 0;
         % should movie be saved (boolean)
         % 0: No
         % 1: Yes
    frames_per_second = 15;
         % number of frames per second for the video generated.
    n_standard_deviations_display = 2;
         % the number of standard deviation to be displayed when a model is
         % presented to the user
    n_modes_display = 2;
         % as above, for number of modes
    perturbation_method = 'CPS warp';
         % the type of perturbation applied to the data
    offset_extent = 0;
         % the extent of perturbation or the level of distortion which
         % accompanies the parameter 'perturbation_method'
    precision_automation_type = 'default';
       % the way to set precision for the optimiser
    do_plot = 0;
         % should bumps be plotted as it changes? (Boolean)
         % (Old comment: this also enables frames being captured which slows down the
         % process. Capturing of frames and plotting can potentially be separated to make
         % this routine more useful).
    movie_type = 1;
         % type of movie to save
         % 0: bumps as pixels. do_plot may be 0 for a pixel clip to be recorded.
         % 1: bump plots. Remember to set do_plot to 1.
    show_warps_while_optimising = 0; 
         % indicates if the warps applied to the data should be dynamically
         % displayed to the user
    movie_name = 'undefined';
          % the name of the AVI movie to be saved. It will be folloed by
          % '.avi' eventually and contain the movie.
    format = 'png';
         % if data is saved as an image, this paramater defines the format
         % as which it shall be saved.
    show_warp_targets = 0;
         % defines whether or not the target of the warps, i.e. that at
         % alignment, should be displayed to the user.
    dynamic_precision = 0;
       % defines whether or not precision is chosen wisely at each stage
    gen_iters = 25;
       % generalisability iterations = 25;
       % specificity iterations
    spec_iters = 25;
       % specificity iterations
    progress_bar = 1;
         % display progress bar or text instead
         % 0: Hierarchical Text
         % 1: Progress bar
         % 2: Quiet
         % 3: Console progress bar
    model_score_type = 'Default';
         % Defined the way in which a model is scored.
    force_first_reference = 0;
         % depending on its numeric value, this parameter defines the
         % reference selection method
    user_menu = 0;
         % indicates the way in which registration parameters are entered.
         % If not via pop-up menu, then probably via top menu in the
         % application.
    number_of_bins = 50;
         % 0: parameters need to be set in file
         % 1: menu input
    background_red = 1;
    background_green = 1;
    background_blue = 1;
         % the background colourn for the plots
    distance_type = 'sum_of_squared_distances';
         % The way to meausre distance from mean
    show_registration_target = 0;
         % indicates if the estimated target of registration should be
         % super-imposed in score plots.
    produce_statistics = 0;
         % should statistics be displayed at the end (boolean)
    see_model_before_each_iteration = 0;
         % should a combined model be produced and displayed before each
         % iteration, this should be set to 1
    show_pixels = 0;
         % show pixels representation (boolean). Mandatory for pixel movie
         % capture
    show_score_vs_iteration = 0;
         % when enabled, evaluation of the objective function will be
         % displayed against the iteration at which it was evaluated.
    find_perfect_warp = 1;
         % This argument should not be tempered with. It defined the need
         % to check for the value and state at conceived registration
         % target
    show_model_after_optimisation = 0;
         % if enabled, then a model produces from the data is displayed
         % after registration.
    show_specificity_statistics = 0;
         % more booleans which are self-explanatory
    verbose_score = 0;
         % boolean: display score or not after optimisation
    show_score_bar = 0;
         % The score bar gives a visual indication of the 'goodness' of a
         % model construced from the data
    average_smooth_enabled = 0;
    gaussian_smooth_enabled = 0;
         % flags indicating the smoothness to be applied, if any at all
    n_points = 5;
         % an argument that is passed to the model optimisation function.
         % it refers to the number of knotpoint, e.g. for random and scale.
    window = 1;
         % the size of the window for average smoothing or sigma for
         % gaussian smoothing
    show_all_score_figure_types = 0;
         % defines the need for an abundant number of plots depicting
         % quality of model versus iteration
    plot_before_and_after = 0;
         % will only ever be used if the two above are set
    bump_width = 0.2;
    bump_width_variation = 0.5;
    bump_height = 0.2;
    bump_height_variation = 0.3;
    bump_position_freedom = 0.9;
    initial_diminish_factor = 2;
    smoothness_factor = 17;
         % set the bump parameters
    precision_required = 1e-7;
        % the precision required from the optimiser
    model_evaluation_method = 1;
                      % 1: Determinant
                      % 0: MDL
                      % 2: Sum of log of Eigenvalues
    show_score_using_model = 0;
    handle = 'Default';
      % give figures a fefault name
    n_modes = 20;
      % maximal number of modes to be accounted for when evaluating a model
    variation_kept = 0.99;
      % keep 99%
      % number of modes in the model to be accounted for
    retain_original_peak = 0;
      % determines whether or not to keep the same maximal data height
    error_verbosity = 'normal';
      % the level of verbosity for errors
    initialise_warps_close_to_target = 0;
      % determines the requirement to start registration with warps near
      % their target (at completion of registration)
    shape_weight = 0.06;
      % the constant weight Ws for the model constructed. 1 is 100%, 0 is
      % 0% and legal range is 0..1
    weighting_normalisation_method = 'edge';
     % the scheme for Ws weighting selection
    run_all_objective_functions = 0;
       % do not try running all the objective functions
    called_from_gui = 0;
       % indicate that GUI was not used in this case
    kernel = 'Exponential';
       % the kernel type for PDF objective function
    save_data = 0;
       % Indicated if the data needs to be saved to files after opening or
       % generation
    peak_type = 'average_position';
         % The type of peak to be retained
    load_data = 0;
       % indicates if data needs to be loaded from file rather than being
       % generated.
    save_bumps = 0;
       % indicated if bumps need to be saved as images (not data)
    draw_warp_curve = 0;
       % An option for overlapping warps display
    draw_curve_by_images = 0;
       % The option for warp target maps (images x iterations)
    overlap_scores = 0;
      % defines the need for scores to be drawn on the same space if
      % different objective functions are compared by being run
      % sequentially
    show_mesh = 0;
     % simply indicates if a mesh of the data needs to be displayed
    show_target = 0;
        % do not show perfectly aligned images estimate
    verbose = 'off';
         % be verbose or not during general-purpose optimisation
    save_stats = 0;
       % indicated if statistics need to be saved as HTML
    display_shape_model = 0;
      % If enabled, shape model is shown
    display_intensity_model = 0;
      % likewise for intensity model
    n_images = 2;
             % how many images to try to register in each set
	n_iterations = 2;
             % how many iterations of the optimisation to run, default was 50.
	n_sets = 1;
             % how many sets of images to register (in order to get better statistics)
	image_width = 50;
             % the number of pixels or the lenght of the data vector
    figure_handle_number = 1;
             % increments applied when figures are saved
    plot_alignment = 0; 
             % Indicates if a plot for alignment and perturbation should be
             % displayed
    rotate_surface = 0;
             % Indicated if the surface should be rotated horizontally at
             % each frame
    draw_final_warps = 0;
             % indicates if final warps need to be displayed
    add_points = 0;
             % point stuffing where difference in warp curve is greatest
    timestamp_images = 0;
             % indicates if output images should be stamped with time and
             % date.
    add_points_cycle = 1;
             % the frequency (measured in terms of iteration) for stuffing
    auto_weighting = 'off';
             % enables preset change in model weighting, one which is
             % defined by the user and not variance or edge.
    automatic_weighting_type = 'default';
             % the type of weighting automation
    show_inverse_projection = 'on';
             % whether or not the examples should be reconstructed using
             % the inverse warps
    show_model_reconstruction_for_warped_images = 'off';
             % whether to display model reconstructions or not
    point_insertion_position = 'highest_model_discrepancy';
             % the position to stuff points
             
else
    called_from_gui = 1;
        % indicate that default values were not set and that the function
        % was called from REGISTER
end

       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input / Output              %%

rand('state',sum(100*clock))
          % seeding for random numbers.

figure_handle_number = 1;
          % A number to be incremented every time a figure is saved. It
          % manages the number of figures that are generated.
warning off;
         % keeps st_pca.m quiet
error_found = 0;
         % status flag
current_viewpoint = 45;
         % the starting horizontal viewpoint for surfaces
         % Optimisation note: initiliase this variables only when surfaces
         % are displayed.         
status = ['                      '];
         % set intitial status message. Messages are aggregated before
         % being returned to the GUI and displayed to the user.
if (run_all_objective_functions == 1),
         % If all objective functions will run sequentially, then some
         % output parameters need to be updated.
    colour = 'r';
        % set initial colour for overlaid plots
    style = ':';
        % similarly to style of line in the overlaid plots
    progress_bar = 2;
end
     % the progress indicators will not work for this routine so disable
     % them

if strcmp(movie_name, 'undefined'),
    movie_name = ['movie_', spline_type, '_', placement_type, '_', objective_function,'_width', num2str(image_width), '_images', num2str(n_images), '_iterations', num2str(n_iterations), '.avi'];
         % name for the AVI to be saved
         % The name will attempt to encompass all construction details
         % It will replace the default name given when the function is not
         % called from the GUI. If called from the GUI, it usually have a
         % user defined name or 'Default.avi'
end


frame_counter = 0;
         % counter for video frames. Do not fiddle with this variable. It
         % keeps track of the current frame number to be captured.
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Others (most are obsolete and are remains of old code.)  %%        
                       
% white_ctr = 0;
         % white counter? Was used for simple bumps in code that is now
         % commented out
% ref_shift = 0.2;
       % shift in reference image??
% max_shift = 0.2;
       % maximum allowed shift?
% step = 0.1;
       % used in two lines below only
% los = zeros(n_images,floor((max_shift-ref_shift)/step));
% his = zeros(n_images,floor((max_shift-ref_shift)/step));
       % not understood yet
% blurred = 0;
       % image blurring has been applied to images (boolean)
% min_error = 0;
% max_error = 1;
       % define error range??
% error_step = 0.1;
% n_error_steps = (max_error-min_error) / error_step;


         
switch model_evaluation_method,
	case 0
        model_evaluation_method = 'minimum_description_length';
	case 1
        model_evaluation_method = 'determinant';     
	case 2
        model_evaluation_method = 'sum_of_log_eigenvalues';
	case 3
        model_evaluation_method = 'specificity';
	case 4
        model_evaluation_method = 'generalisability';    
	otherwise
        error('Unknown model evaluation method. Wrong numeric identifier passed.');
end
         % The conversion above can one day be embedded in the GUI to make
         % everything more readable. Values are overwritten here too.
         
switch progress_bar,
	case 0
     progress_indicator_type = 'Hierarchical Text';
	case 1
     progress_indicator_type = 'Progress Bar';
	case 2
     progress_indicator_type = 'Quiet'; 
	case 3
     progress_indicator_type = 'Console Progress Bar';
	otherwise
     error('Callback value for progress indicator is not recognised');
end

         % do type conversions for input to make code below readable
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Begin Set initialisation    %%  

switch user_menu,
  case 0
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Optimisation Controllers: Default   %%       
	
	% Original settings from 2003: The values below will be defaults if user menu does not appear
	
	n_images = 2;
             % how many images to try to register in each set
	n_iterations = 2;
             % how many iterations of the optimisation to run, default was 50.
	n_sets = 1;
             % how many sets of images to register (in order to get better statistics)
	image_width = 50;
             % the number of pixels or the length of the data vector
             
    % A few extra popular choices, added March 2004 to speed up experiments
    % construction and customisation    
    
  case 1
    [n_sets, n_iterations, n_images, image_width] = get_parameters_from_user('Default');
       % Prompts the user for the main registartion parameters using menus

  case 2    
    n_images = 10;
	n_iterations = 50;
	n_sets = 1;
	image_width = 50;
    
  case 3   
    n_images = 5;
	n_iterations = 5;
	n_sets = 1;
	image_width = 50; 
    
  case 4    
    n_images = 10;
	n_iterations = 10;
	n_sets = 1;
	image_width = 50;
    
  case 5   
    n_images = 20;
	n_iterations = 10;
	n_sets = 1;
	image_width = 50; 
    
  case 6
    n_images = 20;
	n_iterations = 1000;
	n_sets = 1;
	image_width = 50;
    
  case 7 
    n_images = 10;
	n_iterations = 1000;
	n_sets = 1;
	image_width = 50;
    
  case 8
    n_images = 10;
	n_iterations = 10000;
	n_sets = 1;
	image_width = 50;
    
  otherwise
    error('Unknown parameters option. Internal failure.');
    
end 

% in all the above cases, meaningful names rather than numbers are a
% possible extension which requires time.


tic;
       % start counting time
for this_set = 1:n_sets
       % for all sets
       % Old: user gets status about set number. Multiple sets generate reliable
       %      statistics
  tic;
       % This is called again (stack of time counters????? Doesn't seem so)
       % This will need to be changed (one tic omitted) for the statistics
       % to work properly. Only one stopwatch can work at a time.
       % Smith: for each set of n_sets images, build a model
       % calculate obj fn values from these models, plot the mean with error
       % bars given by 1 stddev away
  if (strcmp(progress_indicator_type,'Hierarchical Text')),      
    disp(['calculating set ', num2str(this_set), ' out of ', num2str(n_sets), ' sets...']);
      % display command-line status indicators
  elseif (strcmp(progress_indicator_type,'Progress Bar')),
      gui_active(1);
          % will add an abort button
      if (strcmp(getenv('OS'), 'Linux')),
         h = progressbar([], 0, 'Preparing for registration', [['Registration Progress on '], [get_host]]);
      else
         h = progressbar([], 0, 'Preparing for registration', 'Registration Progress');
      end
        % initiate progress bar if necessary
      h = progressbar( h, 0, 'Fetching Data');
  end    

  % initiate score bar if necessary


  if (show_score_bar == 1),
     score_bar = waitbar(0, 'Score');
  end
     

  switch load_data,
      
   case 0
      [imagelist, images, points] = make_1d_bump(n_images, image_width...
          , bump_width, bump_width_variation, bump_height, bump_height_variation...
          , bump_position_freedom, initial_diminish_factor, smoothness_factor...
          , save_bumps, format);
       % create a set of images given several parameters
       
      points = -1 + 2 * (points - 1) / (image_width - 1);
        % Smith: normalise points from -1 to 1
             
      if (called_from_gui == 1 & save_bumps == 1),
        status = [[status];['(**) Images saved     '];['                      ']];
      end
        % notify user using the status bar that data has been saved.
 
       %% keep = 0.999999; NOTE (RSS Feb. 2004): became obsolete and
       %% replaced by <variation kept> which relates to PCA
       
       % Nov. 2003: Is this precision required? Perhaps in optimisation.
       
   case 1
      output = open([[handle],['.mat']]);
      points = output.points;
      images = output.images;
                  % copy the data structures from the file
                  
   case 2
      output = open('Default.mat');
      points = output.points;
      images = output.images; 
      
   case 3
      output = open('flat.mat');
      points = output.points;
      images = output.images; 
      
   case 4
      output = open('plain.mat');
      points = output.points;
      images = output.images;  
      
   case 5
      output = open('plainran.mat');
      points = output.points;
      images = output.images;  
      
   case 6
      [imagelist, images, points] = make_elliptic_bumps(n_images, image_width, save_bumps, format, bump_height, bump_position_freedom, bump_height_variation, bump_width_variation);
       % create a set of images given several parameters
      points = -1 + 2 * (points - 1) / (image_width - 1);
             
      if (called_from_gui == 1 & save_bumps == 1),
        status = [[status];['(**) Images saved     '];['                      ']];
      end
      
   case 7
      [imagelist, images, points, dummy1, dummy2] = make_1d_images(n_images, image_width, 0.2);
       % 0.2. is the value used by Kate for <white width>
      if (save_bumps == 1),
          msgbox('Data save is not supported with simple bumps.');
      end
       % create a set of images given several parameters
      points = -1 + 2 * (points - 1) / (image_width - 1);
             
      if (called_from_gui == 1 & save_bumps == 1),
        status = [[status];['(**) Images not saved '];['                      ']];
      end    
      
   otherwise
      error('Invalid value passed from open submenu. Internal error.');
      
  end     
  
  
  if (average_smooth_enabled),
      images = average_smooth(images, window);
  end  
  
  if (gaussian_smooth_enabled),
      images = gaussian_smooth(images, window); 
  end
       % smooth all images - used to be disabled and set <blurred> to 1
       
  original_image_width = image_width;
    % retain this value in case point insertion is radical

  ref_points_vec = points(:,1);
  ref_image_vec = images(:,1);
  
       % set first image generated to be reference when calculating a
       % perfect warp and also as reference in pair-wise objective
       % functions. This has been made more generalised later on.
  % ref_hi = his(1);
  % ref_lo = los(1);
       % and get its upper and lower boundaries. his and lois were
       % generated and returned as output for simple bumps. They were
       % used to find the perfect warps and defined points of importance in
       % the input data.       
  
  % Now choose the image which will not be changed. Note a possible
% optimisation: only if a reference is forced, then the code below should
% be run to choose the reference. The section below that is to do with
% alignment means however that for target to be shown an immutable image
% must always be chosen

      switch force_first_reference,
       case 0        
           immutable_image = 2;
                % see comments later about the choice of first or second
                % image as an immutable reference. These options will only
                % have an effect if a reference is forced (not as in the
                % default reference in MSD-based algorithms)
                
                % OPTIMISATION NOTE: score can also be set to 0 in all
                % cases without breaking registration. It would then need
                % not be repeated in all the cases below.
                
       case 1
           immutable_image = 1;
           score = 0;
                % the score needs to be initilialised because no previous
                % scoring will be avaliable
                
       case 2
          % closest data to average (mean)
           immutable_image = find_image_closest_to_mean(images, image_width, n_images, 'closest', distance_type);
           score = 0;
           
       case 3
              % random reference
           immutable_image = ceil(rand * n_images);
           score = 0;
           
       case 4
          % closest data to average (mean)
           immutable_image = find_image_closest_to_mean(images, image_width, n_images, 'farthest', distance_type);
           score = 0;
           
       case 5
           % highest peak
           immutable_image = find_image_by_height(images, image_width, n_images, 'highest');
           score = 0;
           
       case 6
           % lowest peak
           immutable_image = find_image_by_height(images, image_width, n_images, 'lowest');
           score = 0;
           
       otherwise
          error('Reference error. Wrong values passed from callback.');
          
      end  

      % NOTE PROBLEM: If teh refrence chosen is the first then evaluation at
      % the start is 0 and the optimisation plot is not suitable
       
       
       
       
  %%%%% Do warp analysis here
  
  % The check for perfect warps does not work for corrective offsets,
  % non-bump data, and data that is significantly smoothed. In all of these
  % cases, the flag below must be set to 0, or else piece wise linear will
  % not work.
  check_perfect_warps = find_perfect_warp;
          % currently always enabled, but computationally expensive and in
          % some cases not necessary. A test for the value of this part to
          % be implemented in the future.
      if (check_perfect_warps == 1),   
          if (strcmp(progress_indicator_type,'Progress Bar')),
             h = progressbar( h, 0, 'Calculating Correct Warps');
          end
		% (Old - November 2003 - Smith)
		  % correctly_warped_points = points;
		  % correctly_warped_images = images;
		  for i = 1:n_images,
		      % Smith: now bung in the right deformations to get them back
              % NOTE: THE VARIABLE ref_image_vec BELOW ASSUMED A FIRST
              % IMAGE REFERENCE although this is not the case for the <force
              % reference> option which makes it the second (implicitly).
              % Somthing richer than a second image might be used in
              % future.
              
              %possible optimisation: do not interpolate for the first
              %image which is usually a reference and is not affected
              %anyhow
              %images(:,i)
		    correctly_warped_points(i,:) = find_target_warp(images(:,immutable_image), images(:,i), image_width);
              % correctly_warped_points(i,:)
              % get the warp that appears most desirable
            correctly_warped_points(i,:) = -1 + 2 * (correctly_warped_points(i,:) - 1) / (image_width - 1);
              % normalise points so that they lie in the range x where -1 < x < 1
		    correctly_warped_points_vec = correctly_warped_points(i,:);              
		    image_vec = images(:,i);
              % RSS: copy original image and copy piece-wise linear warp
		      % Smith: resample at the warped points - point on reg. grid has value of point at same index in warpy grid
		      % need to interpolate
              % size(ref_points_vec)
              % size(image_vec)
              % size(correctly_warped_points_vec')
		    correctly_warped_images(i,:) = interp1(ref_points_vec, image_vec, correctly_warped_points_vec);
		    % correct_model = build_model(correctly_warped_images,correctly_warped_points,1,'','edge',0);
		    % converging_score(i) =
		    % measure_model(correct_model.variances,50);
            % correctly_warped_images(i,:)
          end
		  % figure;
          % plot(log(converging_score));
          % title('Log of score as model converges');
      end
      % figure;
      % (test curve) plot(correctly_warped_points(2,:))
      
  

      if (show_warp_targets == 1),
         warp_handle = figure('Name','On the left: piece-wise linear optimal registration warps; On the right: results of warp');
         % title('On the left: piece-wise linear optimal registration warps; On the right: results of warp');
         save_handle = figure(warp_handle);
         hold on;
         subplot(n_images, 2, 1);
         for i=1:n_images,
             subplot(n_images, 2, i * 2 - 1);
             hold on;
             axis off;
             plot(correctly_warped_points(i,:));
             whitebg(save_handle, [background_red, background_green, background_blue]); 
             axis([-1 image_width -1 1]);
             hold off;
         end
         for i=1:n_images,
             subplot(n_images, 2, i * 2);
             hold on;
             axis off;
             plot(correctly_warped_images(i,:));
             whitebg(save_handle, [background_red, background_green, background_blue]);
             axis([0 image_width 0 1]);
             hold off;
         end
         hold off;
         if (save_in_log == 1),
           if (timestamp_images),
               timestamp_plot;
           end    
           saveas(save_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
           figure_handle_number = figure_handle_number + 1;
         end
      end
      
  
    
  if (initialise_warps_close_to_target == 1),
      % offset_extent = 5;
      % perturbation_method = 'CPS warp';
      
    if (plot_alignment == 1),
        % if a plot is required to show this step of perturbation after
        % alignment the right-hand side solumn will be occupied by the
        % plots
      subplot_alignment_fig = figure('Name',[[handle],[': Left: Aligned images ; Right (if existent): Aligned images after perturbation']]);
      alignment_handle=figure(subplot_alignment_fig);
      for i=1:n_images,
       subplot(n_images, 2, 2 * i - 1);
       hold on;
       plot(correctly_warped_images(i,:));
       whitebg(alignment_handle, [background_red, background_green, background_blue]); 
       % title(['Unwarped image #',num2str(i)]);
       axis([0 image_width 0 1]);
       grid on;
       hold off;
      end
    end
      
      if (offset_extent == 0),
          % if no perturbation is required then the values can be simply
          % copied.
         for i = 1:n_images,
           % copy registration target
          points(:,i) = correctly_warped_points(i,:)';
           % images(:,i) = correctly_warped_images(i,:)';
         end
      else
          % if some noise is required, then add it in the form of Gaussian
          % noise or CPS warps to points and interpolate original images
          % accordingly. To be controlled from GUI.
           for i = 1:n_images,
                tweaked_points(i,:) = apply_perturbation_to_points(correctly_warped_points(i,:), offset_extent, perturbation_method);
		%         figure
		%         hold on
		%         plot(tweaked_points(i,:));
		%         plot(correctly_warped_points(i,:));
		%         hold off
                   % apply noise to the correct points
                % images(:,i) = interp1(ref_points_vec, images(:,i), tweaked_points(i,:))';
                   % interpolation to get current images
                points = tweaked_points';
                   % copy the points found to become current
                   % inconsistent column and row form is due to initial
                   % implementation. Optimally (when our task becomes optimising performance),
                   % all transposes should be ridden of.
           end
           if (plot_alignment == 1),
              for i = 1:n_images,
                perturbed_aligned_images(:,i) = interp1(ref_points_vec, images(:,i), points(:,i));
              end    
                % inerpolate to show some plots at this earlier stage
              alignment_handle = figure(subplot_alignment_fig);
              for i = 1:n_images,
                   subplot(n_images, 2, 2 * i);
                   hold on;
                   plot(perturbed_aligned_images(:,i));
                   whitebg(alignment_handle, [background_red, background_green, background_blue]);
                   % title(['Unwarped image #',num2str(i)]);
                   axis([0 image_width 0 1]);
                   grid on;
                   hold off;
               end
             if (save_in_log == 1),
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(alignment_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
             end
           end
      end
      
  end
  
  
  if (save_data == 1),
      % note that data will be saved with modified warp curves if option is
      % set. This is intentional. Data in this cotext is vectors, not
      % images.
      save([[handle],['.mat']], 'images','points');
      if (called_from_gui == 1),
        status = [[status];['(**) Data saved       '];['                      ']];
      end
  elseif (save_data == 2),
      % if data needs to be saved as default
      save([[handle],['.mat']], 'images','points');
      save(['Default.mat'], 'images','points');
      if (called_from_gui == 1),
        status = [[status];['(**) Data saved (Def) '];['                      ']];
      end
  end


  
  if (retain_original_peak),  
       for current_image = 1:n_images,
           maximum_height(current_image) = max(images(:, current_image));
              % find out the maximum height of each data instance
           
           for j=1:image_width,
               switch peak_type,
                   case 'initial_position'
                         if (images(j,current_image) == maximum_height(current_image)),
                           maximum_position(current_image) = j;
                            % find the position where height is maximal. Currently
                            % this search is very inefficient.  
                         end
                         
                   case 'average_position'
                         if (images(j,current_image) == maximum_height(current_image)),
                           maximum_position(current_image) = j;
                         end
                         % same as above
               end
           end
       end
 
       switch peak_type,
         case 'average_position'
           average_position = floor(mean(maximum_position));
           for current_image = 1:n_images,
             maximum_position(current_image) = average_position;
           end 
       end
   end
           
           
           
           
% plot3(1:2,1:2,images(1,:))
%   
%   
%       s = surface(L, ...
% 	    'EdgeColor','none', ...
% 	    'FaceColor',[0.8 0.2 0.2], ...
% 	    'FaceLighting','phong', ...
% 	    'AmbientStrength',0.3, ...
% 	    'DiffuseStrength',0.6, ... 
% 	    'Clipping','off',...
% 	    'BackFaceLighting','lit', ...
% 	    'SpecularStrength',1.1, ...
% 	    'SpecularColorReflectance',1, ...
% 	    'SpecularExponent',7);
% 	l1 = light('Position',[40 100 20], ...
% 	    'Style','local', ...
% 	    'Color',[0 0.7 0.7]);
% 	l2 = light('Position',[.5 -1 .4], ...
% 	    'Color',[1 1 0]);

% modify the progress bar if required

  if (strcmp(progress_indicator_type,'Progress Bar')),
           h = progressbar( h, 0, 'Pre-processing Data');
  end
  
% depending on whether intensity offset was applied, record the number of
% columns required for the plots
       
  if (apply_intensity_offset),
      vertical_pane_location = 3;
                       % Need is for only 3 columns in plot
  else
      vertical_pane_location = 2;
                       % Need is for only 2 columns in plot
  end
  % Note: intensity offsets are considered a possible way of solving the
  % problem, but not natural form for our input data. Hence, saving of data
  % precedes this offset.


% prepare a figure if a plot of all warp curves needs to be rendered
       
  if (draw_warp_curve == 1),
       overlap_curve_fig = figure('Name','Overlapping Warp Curves');
  end
  
  % if pixel representation is required, prepare the figure
  
  if (show_pixels == 1),
    pixels_figure = figure('Name','Pixel Representation');
  end
  
  % and likewise for meshes/surfaces

  if (show_mesh == 1),
    mesh_subfigure = figure('Name','1-D Data Mesh');
  end 
  
  % if a matrix of warps is required, preapre the figure and draw the
  % initial warps in its left side
  
  if (draw_curve_by_images == 1),
    curve_subplot_fig = figure('Name',[['Warp Curves for all data - objective function: '],[objective_function]]);
    curves = figure(curve_subplot_fig);
    % title('Warp Curves for all data');
    for i=1:n_images,
       subplot(n_images, n_iterations + 1, i * (n_iterations + 1) - n_iterations);
       hold on;
       axis off;
       plot(points(:,i));
       whitebg(curves, [background_red, background_green, background_blue]); 
       plot(correctly_warped_points(i,:),'r-')
       if (i == 1 ),
         legend('Initial warp', 'Target warp'); 
            % only display the legend once at the top
       end
       axis([-1 image_width -1 1]);
       hold off;
    end
  end
  

%%             END OF INITIALISATION               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Warping images, then modelling warped images    %%
%% attempt to register images by warping.          %%



  if (plot_before_and_after),
    for i=1:n_images,
       images_before_offset(:,i) = images(:,i);
       % copy these original values in case 3 columns are necessary for plot representation of the data.   
    end
  end
  
  % draw the original data at the leftmost column
  
  if (do_plot),
    subplot_fig = figure('Name',[[handle],[': Left: Original images ; Right: Images subjected to warping']]);
    H=figure(subplot_fig);
    for i=1:n_images,
       subplot(n_images, vertical_pane_location, (vertical_pane_location * i) - vertical_pane_location + 1);
       hold on;
       plot(images(:,i));
       whitebg(H, [background_red, background_green, background_blue]); 
       % title(['Unwarped image #',num2str(i)]);
       axis([0 image_width 0 1]);
       grid on;
       hold off;
    end
  else
    subplot_fig = 0;
      % pass a void figure to the optimisation functions
  end
  
  % draw the data after intensity offset has been applied (when applicable)
  
  if (apply_intensity_offset),
      for index = 1:n_images,
         peak (index) = max(images(:,index));
      end
      peak_reference = max (peak);  
      for index = 1:n_images,
        peak_difference = peak_reference - max(images(:,index));
        images(:,index) = images(:,index) + peak_difference;
        if (do_plot),
            subplot(n_images, vertical_pane_location, (vertical_pane_location * index) - vertical_pane_location + 2);
            hold on;
            plot(images(:,index));
            % title(['Unwarped image #',num2str(i)]);
            axis([0 image_width 0 1]);
            grid on;
            hold off;
        end
      end        
  end 
%% NOTE: Change here in traversal above (RSS, December 2003)
  
	still_to_go = 1;
          % condition for while loop. This will enable the program to
          % continually try all objective functions.
	
	
	while (still_to_go == 1),
    
    if (run_all_objective_functions == 0)
        still_to_go = 0;
           % break out of while loop after just one iteration
    else
           % If new objective functions are added, an entry must be added
           % here to include it in the more comprehensive cycle.
           
           % try all objective functions sequentially
       switch objective_function,    
           case 'all'
           % meaning that it is in its first iteration
             objective_function = 'model_opt_together';
           case 'model_opt_together'
             objective_function = 'model_opt_separate';    
           case 'model_opt_separate'
             objective_function = 'msd_opt_together';
           case 'msd_opt_together'
             objective_function = 'msd_opt_separate';
           case 'msd_opt_separate'
             objective_function = 'mi_opt_together';
           case 'mi_opt_together'
             objective_function = 'mi_opt_separate';
           case 'mi_opt_separate'
             objective_function = 'nmi_opt_together';
           case 'nmi_opt_together'
             objective_function = 'nmi_opt_separate';
           case 'nmi_opt_separate'          
             objective_function = 'mixed_msd_model';
           case 'mixed_msd_model'
             objective_function = 'model_post_msd';
           case 'model_post_msd'        
             objective_function = 'tfc_eccv';
           case 'tfc_eccv'
             objective_function = 'wavelet_complexity';
           case 'wavelet_complexity'
             objective_function = 'model_pairs_together';
           case 'model_pairs_together'
             objective_function = 'model_pairs_separate';
           case 'model_pairs_separate'
             objective_function = 'void'; 
             still_to_go = 0;
              % signal to break out of the while loop because the test of
              % all objective functions is finished.
           otherwise
             error('Problem with <Run All Objective Functions> option');
       end
    end   
      % warped_images = images
      % size(images)
      for i=1:n_images,
        warped_images(:,i) = interp1(ref_points_vec, images(:,i), points(:,i)')';
      end
      % size(warped_images)
      warped_points = points;
          % set this initial state to plot the examples before change
  
          
       
      
      
  
      if ((strcmp(progress_indicator_type,'Console Progress Bar')) | (strcmp(progress_indicator_type,'Progress Bar')))
         progress_bar_position = 0;
         time_for_this_iteration = 10;
         total_images = n_images * n_iterations;
         step_size = 50 / total_images;
            % these variables are used for some progress statistics
      end   
        % initialise progress bar
        


        
      for n = 1:n_iterations,
        
        if (dynamic_precision == 1),
           iterations_ratio = n_iterations / n;
             % get iterations ratio. In future, looking at the optimisation
             % state and behaviour, for example whether it is stuck,
             % should give better results.
           precision_required = set_precision(iterations_ratio, precision_automation_type);
        end
        
        if (strcmp(auto_weighting, 'on')),
            
           iterations_ratio = n_iterations / n;
           
           shape_weight = set_model_weight(iterations_ratio, automatic_weighting_type);
               % method goes here for setting the weight    
       end
        
        if (show_score_bar == 1),
            % subplot('position',[0,0,1,0.1]);
            Epsilon = 0.0001;
            c_model = build_model(warped_images, warped_points, variation_kept, '', weighting_normalisation_method, shape_weight);
            % waitbar((measure_model(c_model.variances, n_modes, 2, c_model) - log(Epsilon.^n_modes)) / log(Epsilon.^n_modes), score_bar);
            waitbar((measure_model(c_model.variances, n_modes, 'sum_of_log_eigenvalues', c_model) ...
               / log(Epsilon.^n_modes) - 0.8) * 5, score_bar);
            % (measure_model(c_model.variances, n_modes, 2, c_model) - log(Epsilon.^n_modes)) / log(Epsilon.^n_modes) * 100
            % axis off;
            % subplot(1,1,1);
        end            
          
          
          
        if (show_mesh == 1),
          meshview = figure(mesh_subfigure);
          whitebg(meshview, [background_red, background_green, background_blue]);
               % set viewpoint
          [X,Y] = meshgrid(1:n_images * 5 + 2, -image_width / 2 - 1:1:image_width / 2);
               % two zero buffers and 5 rows of each bump
          for i = 1:n_images,
            Z(2:image_width + 1,1 + i * 5) = warped_images(1:image_width,i);
            Z(2:image_width + 1,1 + i * 5 - 1) = warped_images(1:image_width,i);
            Z(2:image_width + 1,1 + i * 5 - 2) = warped_images(1:image_width,i);
            Z(2:image_width + 1,1 + i * 5 - 3) = warped_images(1:image_width,i);
            Z(2:image_width + 1,1 + i * 5 - 4) = warped_images(1:image_width,i);
            % copy the data values
          end
          Z(:,1) = zeros(1:image_width, 1);
                % set bottom buffer
          Z(:, n_images * 5 + 2) = zeros(1:image_width, 1);
                % and upper buffer to close mesh
          Z(1,:) = zeros(1, 1:n_images * 5);
                % set left buffer
          Z(image_width+2,:) = zeros(1, 1:n_images * 5);
                % and right buffer to close mesh      
          switch mesh_type,  
           case 'Surface Red'  
            surf(X,Y,Z,'FaceColor', 'red', 'EdgeColor', 'none');
            camlight left;
            lighting phong;
            axis([0 n_images*5 -image_width/2 image_width/2 0 1]);     
           case 'Surface Blue'     
            surf(X,Y,Z,'FaceColor', 'blue', 'EdgeColor', 'none');
            camlight left;
            lighting phong;
            axis([0 n_images*5 -image_width/2 image_width/2 0 1]);     
           case 'Polygons'   
            mesh(X,Y,Z);
            axis([0 n_images*5 -image_width/2 image_width/2 0 1]);
           case 'Illuminated Surface'
            if (rotate_surface),
              current_viewpoint = current_viewpoint + 1;
               % shift viewpoint slightly
            end
            view(current_viewpoint, 37.5); 
            surface(X,Y,Z, ...
              'EdgeColor','none', ...
              'FaceColor',[0.5 0.5 0.5], ...
              'FaceLighting','phong', ...
              'AmbientStrength',0.3, ...
              'DiffuseStrength',0.6, ... 
              'Clipping','off',...
              'BackFaceLighting','lit', ...
              'SpecularStrength',1.1, ...
              'SpecularColorReflectance',1, ...
              'SpecularExponent',7);
            if (n == 1), % set lights only at the start
                l1 = light('Position',[40 100 20], ...
                  'Style','local', ...
                  'Color',[0 0.5 0.7]);
                l2 = light('Position',[.5 -1 .4], ...
                  'Color',[1 1 0]);
            end
           otherwise
            error ('Mesh type not recognised. Internal error.');
          end
          
          axis off;
          if (movie_type == 2 & save_movie == 1),
                  frame_counter = frame_counter + 1;
                  M(frame_counter) = getframe(meshview);
          end
          if (save_in_log == 1 & n == n_iterations),
             saveas(meshview,[[handle],'-', num2str(figure_handle_number)], 'jpg');
             figure_handle_number = figure_handle_number + 1;
          end
         
        end  
          
        if (strcmp(progress_indicator_type,'Hierarchical Text'))  
           % iterations aim to get good statistical results by averaging
          disp(['    +-Running iteration number ',num2str(n), ' out of ',num2str(n_iterations), ' iterations...']);
        end
        
           % Smith: correct registration
           % attempt to register
               
        if (see_model_before_each_iteration),  
          w_c_model = build_model(warped_images, warped_points, variation_kept,'Optimised warp', weighting_normalisation_method, shape_weight);
                  % create an appearance model...
          model_view = show_combined_model(w_c_model,ref_points_vec, n_standard_deviations_display, n_modes_display, [[handle],[': Combined model before iteration number ', num2str(n)]]);
                  % and show the model before each iteration while saving
                  % the figure handle
          if (movie_type == 3 & save_movie == 1),
                % if movie needs to be saved and it is a model
                  frame_counter = frame_counter + 1;
                  M(frame_counter) = getframe(model_view);
          end        
        end          
        for i = 1:n_images,
           % start iteration time counter   
         if (strcmp(progress_indicator_type,'Hierarchical Text')),
           disp(['        +-Warping image ',num2str(i),' out of ',num2str(n_images), ' images in total...']);
         end
         
         this_step = ((n-1) * n_images) + i;
                    
         if (strcmp(progress_indicator_type,'Progress Bar')),
           tic;
           progress_bar_position = progress_bar_position + step_size;
           steps_remaining = (n_images * n_iterations) - this_step;
           percentage_progress = ceil(this_step / (n_images * n_iterations) * 100);
           minutes = floor(time_for_this_iteration * steps_remaining / 60);
           seconds = rem(floor(time_for_this_iteration *  steps_remaining), 60);
           max_count   = 50;
               % The console-based GUI is 50 units wide
           if (seconds > 9),
               h = progressbar( h, step_size/50,[num2str(percentage_progress),'% Completed ---- Estimated remaining time: ', num2str(minutes), ':', num2str(seconds)]);
           else
               h = progressbar( h, step_size/50,[num2str(percentage_progress),'% Completed ---- Estimated remaining time: ', num2str(minutes), ':0', num2str(seconds)]);
           end    
           if ~gui_active,
                    error_found = 1;
                    break; 
               % If aborted was pressed
           end
              % use the third-party GUI-like progress bar (added February
              % 2004)
         end
         if (strcmp(progress_indicator_type,'Console Progress Bar')),
                               % if console progress bar
           tic;
           progress_bar_position = progress_bar_position + step_size;
           clc;
           disp(['|=================================================|']);
           progress_string='|';       
           for counter = 1:floor(progress_bar_position),
               progress_string = [progress_string, '#'];
           end
           disp(progress_string);
           disp(['|================= ',num2str(floor(progress_bar_position * 2)),'% completed =================|']);
                          % display progress per cent
           steps_remaining = (n_images * n_iterations) - this_step;
           minutes = floor(time_for_this_iteration * steps_remaining / 60);
           seconds = rem(floor(time_for_this_iteration *  steps_remaining), 60);
           disp(' ');
           if (seconds > 9),
             disp(['            Estimated remaining time: ', num2str(minutes), ':', num2str(seconds)]);
                          % show time indicators
           else
             disp(['            Estimated remaining time: ', num2str(minutes), ':0', num2str(seconds)]);
           end
         end
         
         image_vec = warped_images(:,i);
         points_vec = warped_points(:,i);     
         
         if (show_pixels == 1),   
             appended_vector=[];
             for current_image=1:n_images,
               appended_vector=[[appended_vector],[warped_images(:,current_image),[warped_images(:,current_image)]],[warped_images(:,current_image)]];
                   % append bump data 3 times for width which is visible
             end
           
                 %% appended_vector=appended_vector * 256
                 %% use if scale of 0..255 used
	
             G = figure(pixels_figure);
             hold on;
              % title('Bumps as pixels');
             imshow(appended_vector',[0,1]);
                       % show as an image
             if (movie_type == 0 & save_movie == 1),
                                   % Note: already nested in a show_pixels
                                   % 'if' statement
                  frame_counter = frame_counter + 1;
                  M(frame_counter) = getframe(G);
             end
             hold off;
         end
         
         if (movie_type == 1 & save_movie == 1 & do_plot == 1),
                   frame_counter = frame_counter + 1;
                   M(frame_counter) = getframe(H);
         end
                   % capture pixels frame
         
         
         if (i ~= immutable_image | force_reference == 0), % make sure one image (second because score evaluation from a predecessor must exist) remains unchanged if reference is forced        
                       % Somthing richer than a second image might be used
                       % in future e.g. mean image
             switch objective_function,     
                 case 'model_opt_together'
                     
%                   if (i==1 & n==1),
%                          score=0;
%                   end
%                      
%                   last_score = score; 
%                      
%                   % disp(score);
%                   if (i==3) 
%                       figure(101);
%                       plot(image_vec);
%                   end
                    [param_list, warped_point, warped_image, score] = optimise_all_warps_model([warped_images(:,1:i-1), warped_images(:,i+1:n_images)], [warped_points(:,1:i-1), warped_points(:,i+1:n_images)], image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising, precision_required);
                     
                    % if (n==1),
                    %  [param_list, warped_point, warped_image, score] = optimise_all_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, 0);
                    % else
                    %   [param_list, warped_point, warped_image, score] = optimise_all_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, 1, param_list);    
                    % end
                     
                     
%                     if (last_score < score),
%                          disp(num2str((n - 1) * i + i))
%                     end

                 case 'model_pairs_together'    
                     [param_list, warped_point, warped_image, score] = optimise_all_warps_model(warped_images(:,1), warped_points(:,1), image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising, precision_required);    
                 case 'model_pairs_separate'  
                     [param_list, warped_point, warped_image, score] = optimise_warps_model(warped_images(:,1), warped_points(:,1), image_vec, points_vec, spline_type,placement_type,n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising, precision_required);
                 case 'model_opt_separate'
                     [param_list, warped_point, warped_image, score] = optimise_warps_model([warped_images(:,1:i-1), warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type, placement_type,n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising, precision_required);
                 case 'msd_opt_together'         
                     [param_list, warped_point, warped_image, score] = optimise_all_warps_msd(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, precision_required);
                 case 'msd_opt_separate'
                     [param_list, warped_point, warped_image, score] = optimise_warps_msd(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, precision_required);
                 case 'mi_opt_together'       
                     [param_list, warped_point, warped_image, score] = optimise_all_warps_mi(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, number_of_bins, precision_required);
                 case 'mi_opt_separate'
                     [param_list, warped_point, warped_image, score] = optimise_warps_mi(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, number_of_bins, precision_required);
                 case 'nmi_opt_together'        
                     [param_list, warped_point, warped_image, score] = optimise_all_warps_nmi(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, number_of_bins, precision_required);
                 case 'nmi_opt_separate'
                     [param_list, warped_point, warped_image, score] = optimise_warps_nmi(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, number_of_bins, precision_required);
                 case 'mixed_msd_model'
                        if (((n / n_iterations < 2 / 10) & (n / n_iterations >= 1 / 10)) |...
                            ((n / n_iterations < 4 / 10) & (n / n_iterations >= 3 / 10)) |...
                            ((n / n_iterations < 6 / 10) & (n / n_iterations >= 5 / 10)) |...
                            ((n / n_iterations < 8 / 10) & (n / n_iterations >= 7 / 10)) |...
                            ((n / n_iterations < 10 / 10) & (n / n_iterations >= 9 / 10))),
                          [param_list, warped_point, warped_image, score] = optimise_all_warps_msd(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, precision_required);
                        else
                          [param_list, warped_point, warped_image, score] = optimise_all_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising, precision_required);
                        end                    
                 case 'model_post_msd'
                        if (n / n_iterations < 1 / 2),
                            % for the first half
                          [param_list, warped_point, warped_image, score] = optimise_all_warps_msd(ref_image_vec, ref_points_vec, image_vec, points_vec, spline_type, placement_type, n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, precision_required);
                        else
                            % and for the second half
                          [param_list, warped_point, warped_image, score] = optimise_all_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising, precision_required);
                        end
                 case 'tfc_eccv'
                     % msgbox('TFC ECCV 2004: Under Construction.');
                     pdf_evalution_method = 'Default';
                     [param_list, warped_point, warped_image, score] = optimise_groupwise([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, pdf_evalution_method, kernel, precision_required);         
                 case 'wavelet_complexity'
                     % msgbox('Wavelets: Under Construction.');
                     wavelet_evalution_method = 'Default';
                     [param_list, warped_point, warped_image, score] = optimise_wavelet([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot, subplot_fig, i * vertical_pane_location, vertical_pane_location, image_width, n_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, wavelet_evalution_method, filter, filter_level, precision_required);                      
                 case 'void'
                    % This is used for repeated experiments with overlaid
                    % plots. It is not well tested and score may need to be
                    % assigned to a value here
                     break;
                     % It is not obsolutely clear id the above is necessary
                     % or works
                 otherwise
                     error(['Unknown objective function: ', objective_function]);
             end
         else
            % score = 0;
            warped_image = warped_images(:,i);
            warped_point = warped_points(:,i); 
                  % copy some values to make up for the skip
            if (do_plot == 1),
              figure(subplot_fig);
              subplot(n_images, vertical_pane_location , vertical_pane_location * i, 'replace');
              plot(warped_image);
              axis([0 image_width 0 1]);
              drawnow;
            end
                  % draw the plot
         end


         warped_images(:,i) = warped_image;
           

         if (retain_original_peak), 
           warped_images(maximum_position(i), i) = maximum_height(i);
              % always retain the maximum height
         end     
              
              
              
              
         warped_points(:,i) = warped_point;
         
         if (draw_warp_curve == 1),
	       figure(overlap_curve_fig);
           hold on;
           plot(warped_points(:,i));
           whitebg(overlap_curve_fig, [background_red, background_green, background_blue]); 
           hold off;
           if (save_in_log == 1 & n == n_iterations & i == n_images),
               
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(overlap_curve_fig,[[handle],'-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
           end
         end     
         
	     if (draw_curve_by_images == 1),
             figure(curves);
             subplot(n_images, n_iterations + 1, i * (n_iterations + 1) - (n_iterations - n));
             hold on;
             axis off;
             plot(warped_points(:,i));
             whitebg(curves, [background_red, background_green, background_blue]); 
             axis([-1 image_width -1 1]);
             hold off;
             if (save_in_log == 1 & n == n_iterations & i == n_images), 
                 
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(curves,[[handle],'-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
             end
         end
         if (show_pixels == 1 & save_in_log == 1 & n == n_iterations & i == n_images), 
               saveas(G,[[handle],'-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
         end
         if (do_plot == 1 & save_in_log == 1 & n == n_iterations & i == n_images), 
             
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(H,[[handle],'-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
         end
         
         % figure(subplot_fig),title(['Warped images, after image ',num2str(i),' iteration ',num2str(n)]);
         % dummy = waitforbuttonpress;
         if (show_score_using_model == 1),           
             % if evaluation is through model then always overwrite score
            c_model = build_model(warped_images, warped_points, variation_kept, '', weighting_normalisation_method, shape_weight);
                % first build a model
            if (strcmp(model_score_type, 'Default')),
                score = measure_model(c_model.variances, n_modes, model_evaluation_method, c_model);
            elseif (strcmp(model_score_type, 'Specificity')),
                if (mod(n_images * (n - 1) + i, cycle) == 0),
                    % disp('ok')
                  [score, std_spec] = find_specificity(c_model, images, spec_iters, ref_points_vec);
                  standard_deviation(n_images * (n - 1) + i) = std_spec;
                else
                    % set to -1 (or copy from predecessor - the original
                    % intension
                  score = -1; 
                end
            elseif (strcmp(model_score_type, 'Generalisability')), 
               if (mod(n_images * (n - 1) + i, cycle) == 0),
                  [score, std_gen] = find_generalisability(c_model, gen_iters);
                  standard_deviation(n_images * (n - 1) + i) = std_gen;
               else   
                  score = -1;
               end  
            elseif (strcmp(model_score_type, 'Mean MSD')),
                if (mod(n_images * (n - 1) + i, cycle) == 0),
                    score = measure_model_msd(warped_images);
                else
                    score = -1; 
                end
            elseif (strcmp(model_score_type, 'All')),
              % score = measure_model(c_model.variances, n_modes, model_evaluation_method, c_model);
                            % use deafult as score that is primary             
                all_scores = [[find_specificity(c_model, images, spec_iters, ref_points_vec)], [find_generalisability(c_model, gen_iters)], [measure_model_msd(warped_images)]];    
                           % get all scores above and store them in a
                           % matrix
            else
               error(['Wrong scoring method (', [model_score_type], ') passed. Internal error.']);
            end
            
         end    
	
	
         if (n == n_iterations),
            final_score(i, this_set) = score;
                % give the final score for the set if end is reached 
                
            if (strcmp(error_verbosity, 'high')),   
               warped_points(:,i);
               get_inverse_warp(warped_points(:,i)', 'high');
            end
            

            
         end
         
         score_for_iteration(i,n) = score;
               % record the score for this iteration and this image to be later
               % plotted for analysis
         if (strcmp(model_score_type, 'All')),
             all_scores_for_iteration(i,n,1) = all_scores(1);
             all_scores_for_iteration(i,n,2) = all_scores(2);
             all_scores_for_iteration(i,n,3) = all_scores(3);
         end    
         if ((strcmp(progress_indicator_type,'Progress Bar')) | (strcmp(progress_indicator_type,'Console Progress Bar'))),
           times(this_step) = toc;
           time_for_this_iteration = mean (times);
         end
        end  % end:images
        
        
        

        
        
        % The next section deals with point stuffing
        
        %% Note: cycle scale option here... perform this stage every
        %% <add_points_cycle> iterations
        if (mod(n, add_points_cycle) == 0)
          if (add_points),
            if (strcmp(point_insertion_position, 'highest_model_discrepancy')),
                w_c_model = build_model(warped_images, warped_points, variation_kept, 'Optimised warp', weighting_normalisation_method, shape_weight);      
            end    
            for i = 1:n_images,  
              if (strcmp(point_insertion_position, 'highest_model_discrepancy')),
                 mean_intensity = w_c_model.intensity_model.mean;
                 parameters = w_c_model.intensity_model.params(i,:);
                 covariance_matrix = w_c_model.intensity_model.pcs';
                 reconstruction = mean_intensity + parameters * covariance_matrix;
                 reconstructed_image_inverse = interp1(warped_points(:,i), reconstruction, get_inverse_warp(warped_points(:,i)', error_verbosity), 'linear', 0);    
              end    
              for j = 1:image_width - 1, 
                  switch point_insertion_position,
                      case 'steepest_warp'
                         % -Find region in <points> where height difference is high
                         difference(j, i) = warped_points(j + 1, i) - warped_points(j, i);
                         
                      case 'highest_model_discrepancy'
                         difference(j, i) = abs(reconstructed_image_inverse(j) - images(ceil(j / image_width * original_image_width) , i));
                         % size(images)
                         % size(reconstructed_image_inverse)
                         % size(difference)
                  end        
              end
              
              maximum_difference(i) = max(difference(:, i));
              for current_position_number = 1:size(difference(:, i), 1),
                if (difference(current_position_number, i) == maximum_difference(i)),
                     maximum_difference_location(i) = current_position_number;
                        % maximal diffference in position n means that the
                        % internal between points(n) and points(n + 1) is
                        % maximal
                end
              end
            end
            
            
              
            % -Add points in between that range
            % -Set them to some values
            % in between.
            % -Update <image_width> and all data including
            % references
            
            for i = 1:n_images,
                temp_image = warped_images(:, i);
                % temp_original_image = images(:, i);
                temp_points = warped_points(:, i);
                for j = maximum_difference_location(i) + 1:image_width,
                  warped_points(j + 1,i) = temp_points(j);
                  warped_images(j + 1,i) = temp_image(j);
                % images(j + 1,i) = temp_original_image(j);
                end
              
                 % shift points one position onwards where stuffing occurs
                warped_points(maximum_difference_location(i) + 1, i) = warped_points(maximum_difference_location(i), i) + maximum_difference(i) / 2;
                 % NOTE: THE DIFFERENCE SHOULD NOT BE UNIFORM FOR ALL DATA,
                 % It was calculated specifically for each data instance.
                 % set new point to lie at the average of its two neighbouring points 
                warped_images(maximum_difference_location(i) + 1, i) = mean(warped_images(maximum_difference_location(i) + 1, :));
                 % images(maximum_difference_location(i) + 1, i) = mean(warped_images(maximum_difference_location(i) + 1, :));
                 %   % THIS IS DOUBTFUL
                 % -Set new <image> points with similar structures and set their
                 % value to one which is projected from the mean.
                 % NOTE: REFERENCES ARE BEING CHANGED TOO
            end
            
            image_width = image_width + 1;
              
                % modify the size of images
%               size(warped_points)
%               size(warped_images)
%               warped_points
%               warped_images
            
          end % end: point stuffing if
        end % end: cycle if
        
      end  % end:iterations
	  if (strcmp(progress_indicator_type,'Progress Bar')),
                h = progressbar( h, 0, 'Finalising');
      end
      if (show_score_vs_iteration),
          if (show_all_score_figure_types == 1)            
              score_figure = figure('Name', [[handle],[': Score versus iterations for every image warping step - objective function: '],[objective_function]]);
              score_handle = figure(score_figure);
                     % create a new figure to show graph of optimisation for
                     % all steps and all images individually
              hold on;
              
              title('Score versus iterations for every image warping step');
              xlabel('Iteration number');
              ylabel('Score');
              for i=1:n_images,
                    plot(score_for_iteration(i,:));
                end
              whitebg(score_handle, [background_red, background_green, background_blue]);
              hold off;
              
              figure('Name', [[handle],[': Score versus iterations - objective function: '],[objective_function]]);
                     % create a new figure to show graph of optimisation
              hold on;
              title('Score versus iterations');
              xlabel('Iteration number');
              ylabel('Mean score');
              for i=1:n_iterations,
                 mean_score(i) = mean(score_for_iteration(:,i));
              end   
              plot(mean_score);
              whitebg(score_handle, [background_red, background_green, background_blue])
              hold off;
          end   
          
          if (overlap_scores == 0),
             score_figure = figure('Name', [[handle],[': Score versus image warping step - objective function: '],[objective_function]]);
             score_handle = figure(score_figure);
                 % create a new figure to show graph of optimisation
          else
              
             score_figure = figure(101);
             score_handle = figure(score_figure);
          end
          
          hold on;
          
          if (strcmp(model_score_type, 'All')),
            title('Proportional score versus image warping step');
            xlabel('Warping Step');
            ylabel('% Change');  
            position_counter = 0;  
            for i=1:n_iterations,
              for j=1:n_images,
                position_counter = position_counter + 1;
                all_warp_step_score_vector(position_counter,1) = all_scores_for_iteration(j,i,1);
                all_warp_step_score_vector(position_counter,2) = all_scores_for_iteration(j,i,2);
                all_warp_step_score_vector(position_counter,3) = all_scores_for_iteration(j,i,3);
              end
            end
            all_warp_step_percent_vector(1,1) = 0;
                               % no change
            all_warp_step_percent_vector(2,1) = 0;
            all_warp_step_percent_vector(3,1) = 0;
                
            for i = 2:n_iterations * n_images,
                    all_warp_step_percent_vector(1,i) = (all_warp_step_score_vector(i,1) - all_warp_step_score_vector(i-1,1)) / (all_warp_step_score_vector(i-1,1) + 0.001) * 100;
                    all_warp_step_percent_vector(2,i) = (all_warp_step_score_vector(i,2) - all_warp_step_score_vector(i-1,2)) / (all_warp_step_score_vector(i-1,2) + 0.001) * 100;
                    all_warp_step_percent_vector(3,i) = (all_warp_step_score_vector(i,3) - all_warp_step_score_vector(i-1,3)) / (all_warp_step_score_vector(i-1,3) + 0.001) * 100;              
            end          
          else
              title('Score versus image warping step');
              xlabel('Warping Step');
              ylabel('Score');
              position_counter = 0;
              for i=1:n_iterations,
                  for j=1:n_images,
                    position_counter = position_counter + 1;
                    warp_step_score_vector(position_counter) = score_for_iteration(j,i);
                  end
              end
          end
          if (run_all_objective_functions == 1),
              if strcmp(colour,'r'),
                colour = 'g';
              elseif strcmp(colour,'g')
                colour = 'b';
              else
                colour = 'r';
              end 
              
              if strcmp(style,':'),
                style = '-';
              elseif strcmp(style,'-')
                style = '--';
              elseif strcmp(style,'--')
                style = '-.'; 
              else
                style = ':';
              end 
              plot(warp_step_score_vector,[[colour],[style]]);
              
              whitebg(score_handle, [background_red, background_green, background_blue]);
              % The next line needs to be updated when new objective
              % functions are added
              legend('Joint MSD','Sequential MSD','Joint Model-based','Sequential Model-based','Joint MI',...
                     'Sequential MI', 'Joint NMI', 'Sequential NMI', 'PDF', 'Wavelets', 'Joint Model-based Pair-wise', 'Sequential Model-based Pair-wise', 'Void');
          else
                % if a normal single plot will suffice
                if (strcmp(model_score_type, 'All')),
                         plot (all_warp_step_percent_vector(1,:), 'r');
                         plot (all_warp_step_percent_vector(2,:), 'g');
                         plot (all_warp_step_percent_vector(3,:), 'b');
                         whitebg(score_handle, [background_red, background_green, background_blue]);
                         legend('Specificity', 'Generalisability', 'Mean MSD');
                else 
                    	  if (strcmp(progress_indicator_type,'Progress Bar')),
                              h = progressbar( h, 0, 'Generating Scores Plot');
                          end    
                         % find any -1 entries and "destroy" them to get
                         % correct cycle evaluation
                         % optimisation tip: use cycle to pick out the
                         % right entries rather than the naive -1 flag
                         counter = 1;
                         for index=1:n_images*n_iterations,
                             if (warp_step_score_vector(index) ~= -1),
                                 if (strcmp(model_score_type, 'Generalisability') | strcmp(model_score_type, 'Specificity')),
                                   standard_deviation_truncated(counter) = standard_deviation (index);
                                 end   
                                 warp_step_score_vector_truncated(counter) =  warp_step_score_vector(index);
                                 counter = counter + 1;
                             end
                         end    
                         if (strcmp(model_score_type, 'Generalisability') | strcmp(model_score_type, 'Specificity')),
                             errorbar(warp_step_score_vector_truncated, standard_deviation_truncated);
                             whitebg(score_handle, [background_red, background_green, background_blue]);
                             title([model_score_type, ' every ', num2str(cycle), ' steps']);
                         elseif (strcmp(model_score_type, 'Mean MSD'))
                              % for mean MSD
                             plot (warp_step_score_vector_truncated);
                             whitebg(score_handle, [background_red, background_green, background_blue]); 
                             title(['Mean MSD every ', num2str(cycle), ' steps']);
                         else % model defaults
                             plot (warp_step_score_vector);
                             whitebg(score_handle, [background_red, background_green, background_blue]);
                             title(['Registration score versus iterations']);
                         end
                         
                end
                if (show_target == 1),
                    Epsilon = 0.0001;
                    plot (repmat([log(Epsilon.^n_modes)], 1,n_iterations * n_images),'--');
                    legend('Current Evaluation','Log eigenvalues of perfectly registered images');
                end
                if (show_registration_target == 1),
                    w_c_model = build_model(correctly_warped_images', correctly_warped_points', variation_kept, 'Optimised warp', weighting_normalisation_method, shape_weight);      
                    model_score =  measure_model(w_c_model.variances, n_modes, model_evaluation_method, w_c_model);
                    plot(repmat(model_score, 1,n_iterations * n_images),'r-.');
                    if (show_target == 1),
                      legend('Current Evaluation','Log eigenvalues of perfectly registered images','Log eigenvalues of registration target');
                    else 
                      legend('Current Evaluation','Log eigenvalues of registration target');
                    end  
                end    
          end   
          
          hold off;
          if (save_in_log == 1),
              if (timestamp_images),
                 timestamp_plot;
               end
               saveas(score_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
          end
      end
      
      
      
      if (strcmp(show_inverse_projection,'on')),   
              inverse_projection_plot = figure('Name','Warped images projected back to orginal using the inverse warp');   
              for i = 1:n_images,
                warped_image_inverse = interp1(warped_points(:,i), warped_images(:,i), get_inverse_warp(warped_points(:,i)', error_verbosity), 'linear', 0);    
                inverse_projection_plot_handle = figure(inverse_projection_plot);
                hold on;
                subplot(n_images, 1, i);
                plot(warped_image_inverse);
                axis([0 image_width 0 1]);
                whitebg(inverse_projection_plot_handle, [background_red, background_green, background_blue]);
                hold off;
              end
              if (save_in_log == 1),
                        if (timestamp_images),
                             timestamp_plot;
                        end
                        saveas(inverse_projection_plot_handle, [[handle],'-', num2str(figure_handle_number)], 'jpg');
                        figure_handle_number = figure_handle_number + 1;
              end
      end
      
      w_c_model = build_model(warped_images, warped_points, variation_kept, 'Optimised warp', weighting_normalisation_method, shape_weight);      
         % Smith: build model from warped images
      % w_c_model.pcs
      % w_c_model.impdata
      % w_c_model.params
      w_intensity_total_vars = w_c_model.intensity_model.total_var;
      w_shape_total_vars = w_c_model.shape_model.total_var;
         % get these two values to be used later
      model_score(this_set) = measure_model(w_c_model.variances, n_modes, model_evaluation_method, w_c_model);
      msd_score(this_set) = measure_model_msd(warped_images);
      shape_modes(this_set) = w_c_model.n_shape_modes;
      intensity_modes(this_set) = size(w_c_model.intensity_model.pcs,2);
      shape_variance(this_set) = sum(w_c_model.shape_model.variances);
      intensity_variance(this_set) = sum(w_c_model.intensity_model.variances);
         % copy some values to be displayed for statistics
      
      
         
      if (strcmp(show_model_reconstruction_for_warped_images, 'on')),
         reconstruction_figure = figure('Name','Reconstruction of warped images from the model with the inverse warp applied (on the right)');
         for i = 1:n_images,
             mean_intensity = w_c_model.intensity_model.mean;
             parameters = w_c_model.intensity_model.params(i,:);
             covariance_matrix = w_c_model.intensity_model.pcs';
             reconstruction = mean_intensity + parameters * covariance_matrix;
             reconstruction_handle = figure(reconstruction_figure);
             subplot(n_images, 2, i * 2 - 1);
             hold on;
             plot(reconstruction);
             axis([0 image_width 0 1]);
             whitebg(reconstruction_handle, [background_red, background_green, background_blue]);
             hold off;
              % only if the model has been reconstructed, allow for inverse warp
              % applied to it.
             reconstructed_image_inverse = interp1(warped_points(:,i), reconstruction, get_inverse_warp(warped_points(:,i)', error_verbosity), 'linear', 0);    
             subplot(n_images, 2, i * 2);
             hold on;
             plot(reconstructed_image_inverse);
             axis([0 image_width 0 1]);
             whitebg(reconstruction_handle, [background_red, background_green, background_blue]);
             hold off;
         end
         if (save_in_log == 1),
              if (timestamp_images),
                 timestamp_plot;
              end
              saveas(reconstruction_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
              figure_handle_number = figure_handle_number + 1;
         end
      end   
      
      
      if (show_model_after_optimisation == 1),
          model_handle = show_combined_model(w_c_model, ref_points_vec, n_standard_deviations_display, n_modes_display, [[handle],[': Final combined model built from optimised warps - objective function: '],[objective_function]], image_width);
          if (save_in_log == 1),
              
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(model_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
          end
      end
      
      
      if (display_intensity_model == 1),
        [intensity_model.pcs,intensity_model.variances,intensity_model.params, intensity_model.mean, intensity_model.var_r, intensity_model.total_var, intensity_model.impdata] = st_pca(warped_images', variation_kept);
          % perform PCA
        intensity_model.sd = std(intensity_model.params);
              % get standard deviation from the models (of shape and intensity)
        intensity_total_vars = intensity_model.total_var;
        total_intensity_variance = sum(intensity_model.variances);
        fig_title = 'After registration: ';
        intensity_model_handle = show_intensity_model(intensity_model, 2, 2, fig_title, image_width);
        if (save_in_log == 1),
            
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(intensity_model_handle, [[handle], '-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
        end
      end  
      
      if (display_shape_model == 1),  
        [shape_model.pcs,shape_model.variances,shape_model.params, shape_model.mean, shape_model.var_r, shape_model.total_var, shape_model.impdata] = st_pca(warped_points', variation_kept);
        shape_model.sd = std(shape_model.params);
        shape_total_vars = shape_model.total_var;
        total_shape_variance = sum(shape_model.variances);
        shape_model_handle = show_shape_model(shape_model, ref_points_vec, ref_image_vec, 2, image_width);
        if (save_in_log == 1),
            
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(shape_model_handle, [[handle], '-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
        end        
      end  
          
      if (show_specificity_statistics),
        [w_c_model.mean_specificity, w_c_model.std_specificity] = find_specificity(w_c_model, images, spec_iters, ref_points_vec);
        disp(' ');
        disp([['Objective function: '], [objective_function]]);
        disp(['Mean of specificity:               ', num2str(w_c_model.mean_specificity)]);
        disp(['Standard deviation of specificity: ', num2str(w_c_model.std_specificity)]);
      end
      
      if (draw_final_warps),
          final_warps = figure('Name',[['Warp Curves - objective function: '], [objective_function]]);
          final_warps_handle = figure(final_warps);
          % title('Warp Curves at end);
          for i=1:n_images,
            subplot(n_images, 1, i);
            hold on;
            axis off;
            plot(warped_points(:,i));
            plot(correctly_warped_points(i,:),'r-');
            whitebg(final_warps_handle, [background_red, background_green, background_blue])
            if (i == 1 ),
              legend('Final warp', 'Piece-wise linear target warp'); 
              % only display the legend once at the top
            end
            axis([-1 image_width -1 1]);
          end
          hold off;       
          if (save_in_log == 1),
              
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(final_warps_handle, [[handle], '-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
          end  
       end
     
	
	% fig_title = 'Automatically aligned: ';
	% show_shape_model(w_shape_model, ref_points_vec, ref_image_vec, 2, white_width)
	% show_intensity_model(w_intensity_model, 2, white_width, fig_title)
	% show_combined_model(w_c_model, ref_points_vec, 2, white_width, fig_title)
	
	
	
	
	%%             END OF MODEL-BUILDING                    %%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%% Now provide some statistics on the                   %%
	%% simulation/experiment                                %%
	
      
      %% Warping images, then modelling warped images    %%
	%% attempt to register images by warping.          %%
	
      if (plot_before_and_after == 1),
        subplot_fig = figure('Name', [[handle],[': Left: Before registration; Right: After registration - Objective function: '], [objective_function]]);
        H=figure(subplot_fig);
        for i=1:n_images,
           subplot(n_images, vertical_pane_location, (vertical_pane_location * i) - vertical_pane_location + 1);
           hold on;
           plot(images_before_offset(:,i));
           whitebg(H, [background_red, background_green, background_blue]); 
           % title(['Unwarped image #',num2str(i)]);
           axis([0 original_image_width 0 1]);
           grid on;
           hold off;
        end
        if (apply_intensity_offset),
          for index = 1:n_images,
             peak (index) = max(images(:,index));
          end
          peak_reference = max (peak);  
          for index = 1:n_images,
            peak_difference = peak_reference - max(images(:,index));
            images(:,index) = images(:,index) + peak_difference;
            subplot(n_images, vertical_pane_location, (vertical_pane_location * index) - vertical_pane_location + 2);
            hold on;
            plot(images(:,index));
            % title(['Unwarped image #',num2str(i)]);
            axis([0 original_image_width 0 1]);
            grid on;
            hold off;
          end        
        end    
	
          % pass a void figure to the optimisation functions
        for i=1:n_images,
           subplot(n_images, vertical_pane_location, i * vertical_pane_location);
           hold on;
           plot(warped_images(:,i));
           % title(['Unwarped image #',num2str(i)]);
           axis([0 image_width 0 1]);
           grid on;
           hold off;
        end  
        if (save_in_log == 1), 
            
               if (timestamp_images),
                 timestamp_plot;
               end
               saveas(subplot_fig,[[handle],'-', num2str(figure_handle_number)], 'jpg');
               figure_handle_number = figure_handle_number + 1;
        end
      else
        subplot_fig = 0; 
      end
	
      
 end % while 
 
 t = toc;
      
 if (strcmp(progress_indicator_type,'Hierarchical Text')),
         disp(['Time for set ',num2str(this_set),': ',num2str(t)]);
 end
 time(this_set) = t;
 
end % end set



if (produce_statistics == 1),
    disp(' ');
    disp(['======================================================']); 
    disp(['================= SET STATISTICS =====================']);
    disp(['| Objective function:                       ', objective_function]);
	disp(['| Mean of match score:                      ', num2str(mean(final_score(:)))]);	
	disp(['| Standard deviation of match score:        ', num2str(std(final_score(:)))]);
	disp(['| Mean of model score:                      ', num2str(mean(model_score(:)))]);
    disp(['| Log of mean of model score:               ', num2str(log(mean(model_score(:))))]);
	disp(['| Standard deviation of model score:        ', num2str(std(model_score(:)))]);
	disp(['| Mean of MSD score:                        ', num2str(mean(msd_score(:)))]);
	disp(['| Standard deviation of MSD score:          ', num2str(std(msd_score(:)))]);
	disp(['| Mean of shape modes:                      ', num2str(mean(shape_modes(:)))]);
	disp(['| Standard deviation of shape modes:        ', num2str(std(shape_modes(:)))]);
	disp(['| Mean of intensity modes:                  ', num2str(mean(intensity_modes(:)))]);
	disp(['| Standard deviation of intensity modes:    ', num2str(std(intensity_modes(:)))]);
	disp(['| Mean of shape variance:                   ', num2str(mean(shape_variance(:)))]);
    disp(['| Standard deviation of shape variance:     ', num2str(std(shape_variance(:)))]);		
	disp(['| Mean of intensity variance:               ', num2str(mean(intensity_variance(:)))]);	
	disp(['| Standard deviation of intensity variance: ', num2str(std(intensity_variance(:)))]);	
	disp(['| Mean time:                                ', num2str(mean(time(:)))]);	
	disp(['| Standard deviation of time:               ', num2str(std(time(:)))]);	
    t = toc;
         % get time from the counter, initiated by <tic>
    disp(['| Total time:                               ', num2str(t)]);
         % display total time
    disp(['======================================================']); 
end

if (save_stats == 1),
    fid = fopen([[handle], ['.htm']],'w');
    add_html_headers(fid);
    fprintf(fid, ['\n<CENTER><TABLE BORDER=BOX CELLPADDING=10 CELLSPACING=3 BACKGROUND="../../scbg.gif"><TR>']);
    fprintf(fid, ['\n<TD><H2>Ststistics for ', [handle], '</H2>']);
    fprintf(fid, ['\n<BR><U>Objective function:</U> ', objective_function]);
    fprintf(fid, ['\n<BR><U>Mean of match score:</U> ', num2str(mean(final_score(:)))]);
    fprintf(fid, ['\n<BR><U>Standard deviation of match score:</U> ', num2str(std(final_score(:)))]);
    fprintf(fid, ['\n<BR><U>Mean of model score:</U> ', num2str(mean(model_score(:)))]);
    fprintf(fid, ['\n<BR><U>Log of mean of model score:</U> ', num2str(log(mean(model_score(:))))]);
    fprintf(fid, ['\n<BR><U>Standard deviation of model score:</U> ', num2str(std(model_score(:)))]);
    fprintf(fid, ['\n<BR><U>Mean of MSD score:</U> ', num2str(mean(msd_score(:)))]);
    fprintf(fid, ['\n<BR><U>Standard deviation of MSD score:</U> ', num2str(std(msd_score(:)))]);
    fprintf(fid, ['\n<BR><U>Mean of shape modes:</U> ', num2str(mean(shape_modes(:)))]);
    fprintf(fid, ['\n<BR><U>Standard deviation of shape modes:</U> ', num2str(std(shape_modes(:)))]);
    fprintf(fid, ['\n<BR><U>Mean of intensity modes:</U> ', num2str(mean(intensity_modes(:)))]);
    fprintf(fid, ['\n<BR><U>Standard deviation of intensity modes:</U> ', num2str(std(intensity_modes(:)))]);
    fprintf(fid, ['\n<BR><U>Mean of shape variance:</U> ', num2str(mean(shape_variance(:)))]);
    fprintf(fid, ['\n<BR><U>Standard deviation of shape variance:</U> ', num2str(std(shape_variance(:)))]);
    fprintf(fid, ['\n<BR><U>Mean of intensity variance:</U> ', num2str(mean(intensity_variance(:)))]);   
    fprintf(fid, ['\n<BR><U>Standard deviation of intensity variance:</U> ', num2str(std(intensity_variance(:)))]);
    fprintf(fid, ['\n<BR><U>Mean time:</U> ', num2str(mean(time(:)))]);
    fprintf(fid, ['\n<BR><U>Standard deviation of time:</U> ', num2str(std(time(:)))]);    
    fprintf(fid, ['\n<BR><U>Total time:</U> ', num2str(t)]);
    fprintf(fid, ['\n</TD></TR>']);     
    fprintf(fid, '\n</TABLE></CENTER><BR><BR>');
    fclose(fid);
end



if (produce_record == 1),
     % In the GUI it has been renamed to "Display Record"
	clc;
    disp(' ');
    disp(' ');
	disp('===============================================');
	disp(['| Simulation Handle: ', [handle]]);
    disp(['| Description: ', [description]]);
	disp('|==============================================');
	disp(['| AART - ', [get_month], ' ', [get_year], ' - Version ', [get_version], ' - ', [getenv('OS')]]);
	disp(['| Run on machine: ', [get_host], '| Domain: ', [get_domain], '| Image width: ', num2str(image_width)]);
	disp('|==============================================');
	disp(' ');
	disp('|== GENERAL ===================================');
	disp(['| Objective Function: ', objective_function]);
    disp(['| Data type identifier: ', num2str(load_data)]); 
	disp(['| Iterations: ', num2str(n_iterations)]);
	disp(['| Images: ', num2str(n_images)]);
	disp(['| Sets: ', num2str(n_sets)]);
	disp(['| Corrective Offset (boolean): ', num2str(apply_intensity_offset)]);
	disp('|==============================================');
	disp(' ');
	disp('|== SMOOTHING =================================');
	disp(['| Smoothing window: ', num2str(window)]);
	disp(['| Gaussian smoothing (boolean): ', num2str(gaussian_smooth_enabled)]);  
	disp(['| Average smoothing (boolean): ', num2str(average_smooth_enabled)]);
	disp('|==============================================');
	disp(' ');
	disp('|== OBJECTIVE FUNCTION-SPECIFIC ===============');
 	disp(['| Precision Required (All): ', num2str(precision_required)]);
	disp(['| Number of bins (MI, NMI): ', num2str(number_of_bins)]);
	disp(['| Number of model modes: ', num2str(n_modes)]);
	disp(['| Variation Kept (PCA): ', num2str(variation_kept)]);
	disp(['| Weighting normalisation method: ', weighting_normalisation_method]);
	disp(['| Contant weighting value: ', num2str(shape_weight)]);
	disp(['| PDF type: ', kernel]);
	disp(['| Wavelet filter: ', filter]);
	disp(['| Compression level: ', num2str(filter_level) ]); 
	disp(['| Model evaluation method: ',  num2str(model_evaluation_method)]);
	disp(['| Specificity iterations: ', num2str(spec_iters)  ]);
	disp(['| Generalisability iterations: ', num2str(gen_iters)  ]);
	disp(['| Evaluation cycle: ', num2str(cycle)  ]); 
	disp('|==============================================');
	disp(' ');
	disp('|== TRANSFORMATION ============================');
	disp(['| CPS Spline type: ', spline_type  ]);
	disp(['| Knot-point placement method: ', placement_type]);  
	disp(['| Number of knot-points: ', num2str(n_points)]);
    disp(['| Automatic precision (boolean): ', num2str(dynamic_precision)]);
    disp(['| Initialise near target (boolean): ', num2str(initialise_warps_close_to_target)]);
    disp(['| Offset extent: ', num2str(offset_extent)]);
    disp(['| Perturbation method: ', [perturbation_method]]);
    disp(['| Force reference (boolean): ', num2str(force_reference)]);
    disp(['| Reference type: ', num2str(force_first_reference)]);
    disp(['| Distance type: ', num2str(distance_type)]);    
    disp(['| Retain peak (boolean): ', num2str(retain_original_peak)]);    
    disp(['| Peak type: ', [peak_type]]);    
    disp(['| Stuff points (boolean): ', num2str(add_points)]);    
    disp(['| Point stuffing cycle: ', num2str(add_points_cycle)]);
    disp(['| Point stuffing position: ', num2str(point_insertion_position)]);  
    disp(['| Automatic Weighting (Boolean): ', num2str(auto_weighting)]);  
    disp(['| Automatic Weighting Type: ', num2str(automatic_weighting_type)]);
	disp('|==============================================');
	disp(' ');
	disp('|== DATA ======================================');
	disp(['| Base height variation: ', num2str(bump_height)  ]);
	disp(['| Height variation: ', num2str(bump_height_variation) ]);
	disp(['| Base width variation: ', num2str(bump_width)  ]);
	disp(['| Width variation: ', num2str(bump_width_variation) ]);
	disp(['| Position Freedom: ', num2str(bump_position_freedom)  ]);
	disp(['| Smoothness: ', num2str(smoothness_factor)]);
	disp('|==============================================');
	disp(' ');
	disp('|==============================================');
	disp('| Results and outputs to be located by handle |');
	disp('===============================================');    
       
end  

if (save_in_log == 1),
     % if the parameters need to be included in the log file
    
%     fid_test = fopen('aartlo.htm','r');
%     f = fread(fid_test,1);
%     c = char (f')
%     a=2
%     fclose(fid_test);
    if (save_movie == 1),
      movie_name = [[handle], '-1'];
    end 
    % modify the name of the movie to fit the HTML link
    % if multiple movies are added in the future, then the HTML structure
    % needs to fit the extra links (slots)

    if (strcmp(which('aartlog.htm'), '')),
        add_headers = 1;
    else
        add_headers = 0;
    end

    
    
    fid = fopen('aartlog.htm','a');
    if (add_headers == 1),
      add_html_headers(fid);
    end
    
    fprintf(fid, ['\n<BR><CENTER><A HREF="index.htm"><IMG src="../../back.jpg" border="0" ALT="Back to sub-index"></A>     <A HREF="aartlog.htm#top"><IMG src="../../top.jpg" border="0" ALT="Top of page"></A></CENTER>']);
    fprintf(fid, ['\n<BR><BR><A NAME="', [handle],'"><HR SIZE=5></A><H1>Simulation Handle: ', [handle], '</H1>']);
    fprintf(fid, ['\n<BR>', '<CENTER><TABLE BORDER=BOX CELLPADDING=10 CELLSPACING=3 BACKGROUND="../../scbg.gif"><TR><TD>']);
    fprintf(fid, ['\n<H2>Main</H2>']);
    
    if (strcmp(getenv('OS'), 'Linux')), 
      fprintf(fid, ['<BR>AART - ', [get_month], ' ', [get_year], ' - Version ', [get_version], ' - ', [getenv('OS')]]);
      fprintf(fid, ['\n<BR><U>Time</U>: ', [get_date]]); 
      fprintf(fid, ['\n<BR><U>Run on machine:</U> ', [get_host], '<BR><U>Domain:</U> ', [get_domain], '<BR><U>Image width:</U> ', num2str(image_width)]);
    end
    
    fprintf(fid, ['\n<BR><U>Description:</U> ', [description]], '</TD>');
    fprintf(fid, ['\n<TD><H2>General</H2>']);
    fprintf(fid, ['\n<BR><U>Objective Function:</U> ', objective_function]);
    fprintf(fid, ['\n<BR><U>Data type identifier:</U> ', num2str(load_data)]);    
    fprintf(fid, ['\n<BR><U>Iterations:</U> ', num2str(n_iterations)]);
    fprintf(fid, ['\n<BR><U>Images:</U> ', num2str(n_images)]);
    fprintf(fid, ['\n<BR><U>Sets:</U> ', num2str(n_sets)]);
    fprintf(fid, ['\n<BR><U>Corrective Offset (boolean):</U> ', num2str(apply_intensity_offset), '</TD></TR><TR><TD>']);
    
    fprintf(fid, ['\n<BR><BR><H2>Smoothing</H2>']);
    fprintf(fid, ['\n<BR><U>Smoothing window:</U> ', num2str(window)]);
    fprintf(fid, ['\n<BR><U>Gaussian smoothing (boolean):</U> ', num2str(gaussian_smooth_enabled)]);
    fprintf(fid, ['\n<BR><U>Average smoothing (boolean):</U> ', num2str(average_smooth_enabled)]);
    
    fprintf(fid, ['\n<BR><BR><H2>Objective Function-Specific</H2>']);
    fprintf(fid, ['\n<BR><U>Precision Required (All):</U> ', num2str(precision_required)]);
    fprintf(fid, ['\n<BR><U>Number of bins (MI, NMI):</U> ', num2str(number_of_bins)]);
    fprintf(fid, ['\n<BR><U>Number of model modes:</U> ', num2str(n_modes)]);
    fprintf(fid, ['\n<BR><U>Variation Kept (PCA):</U> ', num2str(variation_kept)]);
    fprintf(fid, ['\n<BR><U>Weighting normalisation method:</U> ', weighting_normalisation_method]);
    fprintf(fid, ['\n<BR><U>Contant weighting value:</U> ', num2str(shape_weight)]);
    fprintf(fid, ['\n<BR><U>PDF type:</U> ', kernel]);
    fprintf(fid, ['\n<BR><U>Wavelet filter:</U> ', filter]);
    fprintf(fid, ['\n<BR><U>Compression level:</U> ', num2str(filter_level)]);
    fprintf(fid, ['\n<BR><U>Model evaluation method:</U> ',  num2str(model_evaluation_method)]);
    fprintf(fid, ['\n<BR><U>Specificity iterations:</U> ', num2str(spec_iters)]);
    fprintf(fid, ['\n<BR><U>Generalisability iterations:</U> ', num2str(gen_iters), '</TD><TD>']);
    
    fprintf(fid, ['\n<H2>Transformation</H2>']);
    fprintf(fid, ['\n<BR><U>Evaluation cycle:</U> ', num2str(cycle)]);  
    fprintf(fid, ['\n<BR><U>CPS Spline type:</U> ', spline_type]);
    fprintf(fid, ['\n<BR><U>Knot-point placement method:</U> ', placement_type]);
    fprintf(fid, ['\n<BR><U>Number of knot-points:</U> ', num2str(n_points)]);
    fprintf(fid, ['\n<BR><U>Automatic precision (boolean):</U> ', num2str(dynamic_precision)]);
    fprintf(fid, ['\n<BR><U>Initialise near target (boolean):</U> ', num2str(initialise_warps_close_to_target)]);
    fprintf(fid, ['\n<BR><U>Offset extent:</U> ', num2str(offset_extent)]);
    fprintf(fid, ['\n<BR><U>Perturbation method:</U> ', [perturbation_method]]);
    fprintf(fid, ['\n<BR><U>Force reference (boolean):</U> ', num2str(force_reference)]);
    fprintf(fid, ['\n<BR><U>Reference type:</U> ', num2str(force_first_reference)]);
    fprintf(fid, ['\n<BR><U>Distance type:</U> ', num2str(distance_type)]);    
    fprintf(fid, ['\n<BR><U>Retain peak (boolean):</U> ', num2str(retain_original_peak)]);    
    fprintf(fid, ['\n<BR><U>Peak type:</U> ', [peak_type]]);
    fprintf(fid, ['\n<BR><U>Stuff points (boolean):</U> ', num2str(add_points)]);    
    fprintf(fid, ['\n<BR><U>Point stuffing cycle:</U> ', num2str(add_points_cycle)]); 
    fprintf(fid, ['\n<BR><U>Point stuffing position:</U> ', num2str(point_insertion_position)]);  
    fprintf(fid, ['\n<BR><U>Automatic Weighting (Boolean):</U> ', num2str(auto_weighting)]);  
    fprintf(fid, ['\n<BR><U>Automatic Weighting Type:</U> ', num2str(automatic_weighting_type)]);
    
    fprintf(fid, ['\n<BR><BR><H2>Data</H2>']);    
    fprintf(fid, ['\n<BR><U>Base height variation:</U> ', num2str(bump_height)]);
    fprintf(fid, ['\n<BR><U>Height variation:</U> ', num2str(bump_height_variation)]);
    fprintf(fid, ['\n<BR><U>Base width variation:</U> ', num2str(bump_width)]);
    fprintf(fid, ['\n<BR><U>Width variation:</U> ', num2str(bump_width_variation)]);    
    fprintf(fid, ['\n<BR><U>Position Freedom:</U> ', num2str(bump_position_freedom)]);
    fprintf(fid, ['\n<BR><U>Smoothness:</U> ', num2str(smoothness_factor)]);
    fprintf(fid, ['\n</TD></TR>']);
    
    fprintf(fid, '\n<TR><TD><H2>Image Preview 1</H2><BR><A HREF="');
    fprintf(fid, [[handle],'-1.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
    fprintf(fid, [[handle],'-1.jpg"></A></TD><TD><H2>Image Preview 2</H2><BR><A HREF="']);
    fprintf(fid, [[handle],'-2.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
    fprintf(fid, [[handle], '-2.jpg"></A></TD></TR>']);
    fprintf(fid, '\n<TR><TD><H2>Image Preview 3</H2><BR><A HREF="');
    fprintf(fid, [[handle],'-3.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
    fprintf(fid, [[handle],'-3.jpg"></A></TD><TD><H2>Image Preview 4</H2><BR><A HREF="']);
    fprintf(fid, [[handle],'-4.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
    fprintf(fid, [[handle], '-4.jpg"></A></TD></TR>']);   
    fprintf(fid, '\n<TR><TD><H2>Image Preview 5</H2><BR><A HREF="');
    fprintf(fid, [[handle],'-5.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
    fprintf(fid, [[handle],'-5.jpg"></A></TD><TD><H2>Image Preview 6</H2><BR><A HREF="']);
    fprintf(fid, [[handle],'-6.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
    fprintf(fid, [[handle], '-6.jpg"></A></TD></TR>']);
    fprintf(fid, '\n<TR><TD><A HREF="');
    fprintf(fid, [[handle],'-1.avi">Video #1</A></TD><TD><A HREF="']);
    fprintf(fid, [[handle],'-2.avi">Video #2</A></TD></TR>']); 
    fprintf(fid, '\n<TR><TD><A HREF="');
    fprintf(fid, [[handle],'.mat">Input Data</A></TD><TD><A HREF="']);
    fprintf(fid, [[handle],'.htm">Statistics</A></TD></TR>']);     
    fprintf(fid, '\n</TABLE></CENTER><BR><BR>');
    fclose(fid);
    
      % save in index file

    fid = fopen('index.htm','a');
    
    if (add_headers == 1),
      add_html_headers(fid);
      add_table_index(fid);
    end
    
    if (strcmp(getenv('OS'), 'Linux')),
      fprintf(fid, ['\n<P><CENTER><TABLE WIDTH=100%% BORDER=BOX CELLPADDING=10 CELLSPACING=3 BACKGROUND="../../scbg.gif"><TD ALIGN=CENTER WIDTH=25%%><A HREF="aartlog.htm#', [handle], '">', [handle], '</A></TD><TD ALIGN=CENTER WIDTH=25%%>', [get_date], '</TD><TD ALIGN=CENTER WIDTH=25%%>']);
      fprintf(fid, ['\n', [description], '</TD><TD 25%%>  ']);
    else
      fprintf(fid, ['\n<P><CENTER><TABLE WIDTH=100%% BORDER=BOX CELLPADDING=10 CELLSPACING=3 BACKGROUND="../../scbg.gif"><TD ALIGN=CENTER WIDTH=33%%><A HREF="aartlog.htm#', [handle], '">', [handle], '</A> ', '</TD><TD ALIGN=CENTER WIDTH=33%%>']);      
      fprintf(fid, ['\n', [description], '</TD><TD 33%%>  ']);
    end
    
    
    
    for k = 1:figure_handle_number - 1,
        fprintf(fid, ['<IMG SRC="../../green_square.jpg"> ']); 
    end
    fprintf(fid, ['</TD></TABLE></CENTER></P>']);   
    fclose(fid);
end

% a=figure('Name','test');
% b=figure(a);
% hold on;
% plot([-1,2,3]);
% hold off;
% saveas(b,[[handle],'-', num2str(figure_handle_number)], 'jpg');
% figure_handle_number = figure_handle_number + 1;


%% END OF STATISTICS                                  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if (save_movie == 1),
  if (strcmp(progress_indicator_type,'Progress Bar')),
           h = progressbar( h, 0, 'Saving Video');
  end  
  movie(M);
  movie2avi(M, movie_name, 'FPS', frames_per_second);
  if (called_from_gui == 1),
      % send this to the GUI status bar
      status = [[status];['(**) Video saved      '];['                      ']];
  else
      % indicate to user that a video has been save
        disp(' ');
        disp([movie_name, ' saved']);
          % process the movie and save it
  end
end

if (strcmp(progress_indicator_type,'Progress Bar')),
   progressbar( h,-1 );
      % destroy progress bar
end

if (error_found == 0),
              % if the error flag is not set
	if (called_from_gui == 1),
	  status = [[status];['Registration Completed'];
             ['                      '];       
             ['Status: OK            ']];
         % send message
	else
	  status = 'Status: OK'; 
         % indicate that everything was okay
	end
elseif (error_found == 1),
	if (called_from_gui == 1),
	  status = [[status];['Registration Completed'];
             ['                      '];       
             ['Status: Error         '];
             ['                      '];
             ['Operation Aborted     ']];
         % send message about error
	else
	  status = 'Status: Operation Aborted'; 
         % indicate that something had gone wrong
    end    
else
	if (called_from_gui == 1),
	  status = [[status];['Registration Completed'];
             ['                      '];       
             ['Status: Error         '];
             ['                      '];
             ['Exit code unknown     ']];
         % send message about error
	else
	  status = 'Status: Error'; 
         % indicate that something had gone wrong
    end
end    

if (exit_at_end == 1),
    exit;
end  

if (close_at_end == 1),
    close;
end