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);

% 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.
%
%    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: February 18th, 2004
%      -Revision:    0.6.6
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% variables (parameters) initialisation           %%
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input Defaults              %%   

if (nargin == 0),
    disp(' ');
    disp('Argument(s) ommited. Using defaults instead.');
    disp(' ');
    spline_type = 'multi_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';
         % controls knotpoint 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
    
    save_movie = 0;
         % should movie be saved (boolean)
         % 0: No
         % 1: Yes
    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.
    movie_name = 'undefined';
    format = 'png';
    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_tyype = 'Deafault';
    user_menu = 0;
    number_of_bins=50;
         % 0: parameters need to be set in file
         % 1: menu input
    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;
    show_model_after_optimisation = 0;
    show_specificity_statistics = 0;
         % more booleans which are self-explanatory
    verbose_score = 0;
         % boolean: display score or not after optimisation
    average_smooth_enabled = 0;
    gaussian_smooth_enabled = 0;
    window = 1;
    
    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

    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;
    variation_kept = 0.99;
      % keep 99%
      % number of modes in the model to be accounted for
    shape_weight = 0.06;  
    weighting_normalisation_method = 'edge';
    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';
    save_data = 0;
    load_data = 0;
    save_bumps = 0;
    draw_warp_curve = 0;
    draw_curve_by_images = 0;
    overlap_scores = 0;
    show_target = 0;
        % do not show perfectly aligned images estimate
else
    called_from_gui = 1;
        % indicate that default values were not set and that the function
        % was called from REGISTER
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Optimisation Controllers    %%       

% 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 lenght of the data vector
n_points = 5;
         % an argument that is passed to the model optimisation function.
         % meaning still unknown.

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

warning off;
         % keeps st_pca.m quiet              
error_found = 0;
         % status flag
status = ['                      '];
         % set intiaial status packet
if (run_all_objective_functions == 1),
    colour = 'r';
        % set initial colour for overlaid plots
    style = ':';
        % similarly to style of line
    progress_bar = 2;
end
     % the progress indicators will not work for this routine so disable
     % them


verbose = 'off';
         % be verbose or not during general-purpose optimisation
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
end
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Others (most are obsolete)  %%        
                       
% 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;

frame_counter = 0;
         % counter for video frames. Do not fiddle with this variable.

if (progress_bar == 0),
 progress_indicator_type = 'Hierarchical Text';
elseif (progress_bar == 1),
 progress_indicator_type = 'Progress Bar';
elseif (progress_bar == 2),
 progress_indicator_type = 'Quiet'; 
elseif (progress_bar == 3),
 progress_indicator_type = 'Console Progress Bar';
else
 error('Callback value for progress indicator is not recognised');
end

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

if (user_menu == 1),
  [n_sets, n_iterations, n_images, image_width] = get_parameters_from_user('Default');
end 

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...']);
  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
  end    

  
  if (strcmp(progress_indicator_type,'Progress Bar')),
           h = progressbar( h, 0, 'Fetching Data');
  end
          

  if (load_data == 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
      if (called_from_gui == 1 & save_bumps == 1),
        status = [[status];['(**) Images saved     '];['                      ']];
      end
  else
      output = open([[handle],['.mat']]);
      points = output.points;
      images = output.images;
                  % copy the data structures from the file
  end     
  
  if (save_data == 1)
      save([[handle],['.mat']], 'images','points');
      if (called_from_gui == 1),
        status = [[status];['(**) Data saved       '];['                      ']];
      end
  end
  
  figure;
  [X,Y] = meshgrid(n_images, -image_width/2:1:image_width/2-1);
  Z = zeros(n_images, image_width);
size(X)
size(Y)
size(Z)
size(images)
  for i = 1:n_images,
    Z(:,i) = images(:,i);
  end  
  mesh(X,Y,Z);
  axis([-30 30 -30 30 -1 1]);
  % plot3(1:2,1:2,images(1,:))
  
  
  if (strcmp(progress_indicator_type,'Progress Bar')),
           h = progressbar( h, 0, 'Pre-processing Data');
  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
       
  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

  points = -1 + 2 * (points - 1) / (image_width - 1);
       % Smith: normalise points from -1 to 1
       %% keep = 0.999999; NOTE (RSS Feb. 2004): became obsolete
       % Is this precision required? Perhaps in optimisation.
  ref_points_vec = points(:,1);
  ref_image_vec = images(:,1);
       % set first image generated to be reference
  % 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

       
  if (draw_warp_curve == 1),
       overlap_curve_fig = figure('Name','Overlapping Warp Curves');
  end
  
  if (show_pixels == 1),
    pixels_figure = figure('Name','Pixel Representation');
  end
  
  
  
  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);
    for i=1:n_images,
       subplot(n_images, n_iterations + 1, i * (n_iterations + 1) - n_iterations);
       hold on;
       axis off;
       plot(points(:,i));
       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);
    end
  end
  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));
       % 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
  
  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
           % try all objective functions sequentially
       if strcmp(objective_function, 'all'),
           % meaning that it is in its first iteration
           objective_function = 'model_opt_together';
       elseif strcmp(objective_function, 'model_opt_together'),
           objective_function = 'model_opt_separate';    
       elseif strcmp(objective_function, 'model_opt_separate'),
           objective_function = 'msd_opt_together';
       elseif strcmp(objective_function, 'msd_opt_together'),
           objective_function = 'msd_opt_separate';
       elseif strcmp(objective_function, 'msd_opt_separate'),
           objective_function = 'mi_opt_together';
       elseif strcmp(objective_function, 'mi_opt_together'),
           objective_function = 'mi_opt_separate';
       elseif strcmp(objective_function, 'mi_opt_separate'),
           objective_function = 'nmi_opt_together';
       elseif strcmp(objective_function, 'nmi_opt_together'),
           objective_function = 'nmi_opt_separate';
       elseif strcmp(objective_function, 'nmi_opt_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.
       else
           error('Problem with <Run All Objective Functions> option');
       end
    end   
    
      warped_images = images;
      warped_points = points;
          % set this initial state to plot the examples before change
	
	
	% correctly_warped_points = points;
	% correctly_warped_images = images;
	% for i=1:n_images
	% % now bung in the right deformations to get them back
	% correctly_warped_points(:,i) = linear_warp(ref_points_vec, los(i,white_ctr), his(i,white_ctr), ref_lo, ref_hi);
	% correctly_warped_points_vec = correctly_warped_points(:,i);
	% image_vec = images(:,i);
	% % resample at the warped points - point on reg. grid has value of point at same index in warpy grid
	% % need to interpolate
	% 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);
	% end
	% figure,plot(log(converging_score)),title('Log of score as model converges');
	
      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;
      end   
        % initialise progress bar
      for n = 1:n_iterations,
        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);
                  % create an appearance model...
          show_combined_model(w_c_model,ref_points_vec, 2, 2, [[handle],[': Combined model before iteration number ', num2str(n)]]);
                  % and show the model before each iteration
        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;
           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,['            Estimated remaining time: ', num2str(minutes), ':', num2str(seconds)]);
           else
               h = progressbar( h, step_size/50,['            Estimated remaining time: ', num2str(minutes), ':0', num2str(seconds)]);
           end    
           if ~gui_active,
                    error_found = 1;
                    break; % if aborted was pressed
           end
              % use the third-part 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),
                  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    (strcmp(objective_function,'model_opt_together')),
                 [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);
         elseif(strcmp(objective_function,'void')),        
                 break;
         elseif(strcmp(objective_function,'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);
         elseif(strcmp(objective_function,'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);
         elseif(strcmp(objective_function,'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);
         elseif(strcmp(objective_function,'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);
         elseif(strcmp(objective_function,'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);
         elseif(strcmp(objective_function,'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);
         elseif(strcmp(objective_function,'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);
         elseif(strcmp(objective_function,'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);
                    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);
                    end                    
         elseif(strcmp(objective_function,'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);
                    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);
                    end
         elseif(strcmp(objective_function,'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);         
         elseif(strcmp(objective_function,'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);         
         
         else
                 error(['Unknown objective function: ', objective_function]);
         end
         
         if (draw_warp_curve == 1),
	       figure(overlap_curve_fig);
           hold on;
           plot(warped_points);   
           hold off;
         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));
             axis([-1 image_width -1 1]);
             hold off;
         end 
         
         warped_images(:,i) = warped_image;
         warped_points(:,i) = warped_point;
	
         % 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')),  
                score = find_specificity(c_model, images, spec_iters, ref_points_vec);
            elseif (strcmp(model_score_type, 'Generalisability')),  
                score = find_generalisability(c_model, gen_iters);
            elseif (strcmp(model_score_type, 'Mean MSD')),
                score = measure_model_msd(warped_images);
            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    
         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
      end  % end:iterations
	  if (strcmp(progress_indicator_type,'Progress Bar')),
                h = progressbar( h, 0, 'Finalising');
      end
      if (show_score_vs_iteration),
                    
%           figure('Name', [[handle],[': Score versus iterations for every image warping step - objective function: '],[objective_function]]);
%                  % 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;
%           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);
%           hold off;
          
          if (overlap_scores == 0),
             figure('Name', [[handle],[': Score versus image warping step - objective function: '],[objective_function]]);
                 % create a new figure to show graph of optimisation
          else
              
             figure(101);
          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]]);
              legend('Sequential MSD','Joint MSD','Sequential Model-based','Joint Model-based','Squential MI',...
                     'Joint MI', 'Sequential NMI', 'Joint NMI','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');
                         legend('Specificity', 'Generalisability', 'Mean MSD');
                else 
                         plot (warp_step_score_vector); 
                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
          end   
          hold off;
      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_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 (show_model_after_optimisation),
          show_combined_model(w_c_model,ref_points_vec, 2, 2, [[handle],[': Final combined model built from optimised warps - objective function: '],[objective_function]]);
      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
     
	
	% 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),
        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));
           % title(['Unwarped image #',num2str(i)]);
           axis([0 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 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   
      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 (strcmp(progress_indicator_type,'Progress Bar')),
   progressbar( h,-1 );
      % destroy progress bar
end

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(['| 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_movie == 1)
  movie(M);
  movie2avi(M,movie_name);
  disp(' ');
  disp([movie_name, ' saved']);
          % process the movie and save it
  if (called_from_gui == 1),
      status = [[status];['(**) Video saved      '];['                      ']];
  end
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    

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