build_1d_bump_model

 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
 ==============================================================

0001 function [status] = build_1d_bump_model( ...
0002     spline_type ...
0003     , placement_type ...
0004     , objective_function ...
0005     , apply_intensity_offset ...
0006     , save_movie ...
0007     , do_plot ...
0008     , movie_type ...
0009     , progress_bar ...
0010     , user_menu ...
0011     , produce_statistics ...
0012     , see_model_before_each_iteration ...
0013     , show_pixels ...
0014     , show_score_vs_iteration ...
0015     , show_model_after_optimisation ...
0016     , show_specificity_statistics ...
0017     , verbose_score ...
0018     , movie_name ...
0019     , window ...
0020     , average_smooth_enabled ...
0021     , gaussian_smooth_enabled ...
0022     , bump_width ...
0023     , bump_width_variation ...
0024     , bump_height ...
0025     , bump_height_variation ...
0026     , bump_position_freedom ...
0027     , initial_diminish_factor ...
0028     , smoothness_factor ...
0029     , plot_before_and_after ...
0030     , number_of_bins ...
0031     , model_evaluation_method ...
0032     , show_score_using_model ...
0033     , handle ...
0034     , n_modes ...
0035     , shape_weight ...
0036     , weighting_normalisation_method ...
0037     , variation_kept ...
0038     , run_all_objective_functions ...
0039     , kernel ...
0040     , save_data ...
0041     , load_data ...
0042     , save_bumps ...
0043     , draw_warp_curve ...
0044     , overlap_scores ...
0045     , show_target ...
0046     , format ...
0047     , filter ...
0048     , filter_level ...
0049     , model_score_type ...
0050     , gen_iters ...
0051     , spec_iters ...
0052     , draw_curve_by_images ...
0053     , show_mesh ...
0054     , mesh_type ...
0055     , exit_at_end ...
0056     , show_all_score_figure_types ...
0057     , close_at_end ...
0058     , frames_per_second ...
0059     , show_score_bar ...
0060     , n_points ...
0061     , show_warp_targets ...
0062     , show_registration_target ...
0063     , produce_record ...
0064     , initialise_warps_close_to_target ...
0065     , save_in_log ...
0066     , perturbation_method ...
0067     , offset_extent ...
0068     , force_reference ...
0069     , cycle ...
0070     , verbose ...
0071     , save_stats ...
0072     , n_modes_display ...
0073     , n_standard_deviations_display ...
0074     , show_warps_while_optimising ...
0075     , description ...
0076     , find_perfect_warp ...
0077     , force_first_reference ...
0078     , display_shape_model ...
0079     , display_intensity_model ...
0080     , plot_alignment ...
0081     , error_verbosity ...
0082     , precision_required ...
0083     , background_red ...
0084     , background_green ...
0085     , background_blue ...
0086     , dynamic_precision ...
0087     , precision_automation_type ...
0088     , rotate_surface ...
0089     , draw_final_warps ...
0090     , distance_type ...
0091     , retain_original_peak ...
0092     , peak_type ...
0093     , add_points ...
0094     , add_points_cycle ...    
0095    )
0096 
0097 % BUILD_1D_BUMP_MODEL: Builds a 1-D model for some bump data
0098 %                      and returns some related statistics.
0099 %
0100 % Code originally written by Katherine Smith, 2003
0101 %
0102 %    GENERAL
0103 %
0104 %      Omitted. Much of the explanation is included in the section where
0105 %      default arguments are being assigned.
0106 %
0107 %    INPUT/S
0108 %
0109 %      -- No arguments will lead to defaults
0110 %
0111 %      -spline_type:
0112 %           The type of spline to be used. 'single_point' or
0113 %           'multi_point' at present.
0114 %
0115 %      -placement_type:
0116 %           Placement type for 'multi_point':
0117 %                    'grid' or 'edges' or 'random'
0118 %           placement type for 'single_point':
0119 %                    'overlap_grid' or 'edges_and_scale'
0120 %                    or 'random_and_scale'.
0121 %
0122 %      -objective_function:
0123 %           'msd_opt_together' or 'model_opt_separare' or
0124 %           'msd_opt_separate' or 'model_opt_together'.
0125 %
0126 %      .....
0127 %
0128 %    OUTPUT/S
0129 %
0130 %      -status:
0131 %           The status of the operation of this routine (to be returned at
0132 %           the end
0133 %
0134 %      -Figures produced and results displayed are main
0135 %       data to be utilised.
0136 %
0137 %    PENDING WORK
0138 %
0139 %      -
0140 %
0141 %    KNOWN BUG/S
0142 %
0143 %      -
0144 %
0145 %    COMMENT/S
0146 %
0147 %      -Different video outputs are generated but controlling which
0148 %       one is chosen is not trivial. It currently just does the
0149 %       job of experimentation and showing results as pixels
0150 %       or bumps.
0151 %
0152 %    RELATED FUNCTION/S
0153 %
0154 %      See other functions in same directory, filenames
0155 %      beginning with "build" in particular. There seem to be much
0156 %      overlap due to cut-and-paste.
0157 %
0158 %    ABOUT
0159 %
0160 %      -Created:     December 29rd, 2003
0161 %      -Last update: March 30th, 2004
0162 %      -Revision:    0.9.2
0163 %      -Author:      R. S. Schestowitz, University of Manchester
0164 % ==============================================================
0165       
0166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0167 %% variables (parameters) initialisation           %%
0168   
0169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0170 %% Input Defaults              %%
0171 
0172 % When the function is not called from the GUI, various arguments must be
0173 % initialised so that the function performs a registration task.
0174 
0175 % A possible extension would allow for _partial_ arguments to be given so
0176 % that a script can run different experiments and then logs are necessary
0177 % too.
0178 
0179 if (nargin == 0),
0180     disp(' ');
0181     disp('Argument(s) omitted. Using defaults instead.');
0182     disp(' ');
0183          % notify user about default values being used.
0184     spline_type = 'single_point';
0185          % SMITH: spline type 'single_point' or 'multi_point'
0186          % SMITH: placement type for multipoint: 'grid' or 'edges' or 'random'
0187          % SMITH: placement type for singlepoint: 'overlap_grid' or 'edges_and_scal or 'random_and_scale'
0188     placement_type = 'random_and_scale';
0189          % controls knot-point selection on the grid
0190     objective_function = 'model_opt_together';
0191          % choice of objective function
0192     apply_intensity_offset = 1;
0193          % 0: Offset should not be applied to data
0194          % 1: Offset enabled
0195     exit_at_end = 0;
0196          % indicates whether AART quits at the end of registration
0197     description = '';
0198          % a name provided for describing registration
0199     cycle = 0;
0200           % cycle size for user display of objective function evaluation
0201     produce_record = 0;
0202           % indicates if results should be displayed to the user in the
0203           % console
0204     save_in_log = 0;
0205           % indicates if an HTML file or network of files will be used to
0206           % record the results
0207     close_at_end = 0;
0208           % indicates if the GUI should be closed when registration is done
0209     save_movie = 0;
0210           % a flag for movie output
0211     force_reference = 0;
0212          % should movie be saved (boolean)
0213          % 0: No
0214          % 1: Yes
0215     frames_per_second = 15;
0216          % number of frames per second for the video generated.
0217     n_standard_deviations_display = 2;
0218          % the number of standard deviation to be displayed when a model is
0219          % presented to the user
0220     n_modes_display = 2;
0221          % as above, for number of modes
0222     perturbation_method = 'CPS warp';
0223          % the type of perturbation applied to the data
0224     offset_extent = 0;
0225          % the extent of perturbation or the level of distortion which
0226          % accompanies the parameter 'perturbation_method'
0227     precision_automation_type = 'default';
0228        % the way to set precision for the optimiser
0229     do_plot = 0;
0230          % should bumps be plotted as it changes? (Boolean)
0231          % (Old comment: this also enables frames being captured which slows down the
0232          % process. Capturing of frames and plotting can potentially be separated to make
0233          % this routine more useful).
0234     movie_type = 1;
0235          % type of movie to save
0236          % 0: bumps as pixels. do_plot may be 0 for a pixel clip to be recorded.
0237          % 1: bump plots. Remember to set do_plot to 1.
0238     show_warps_while_optimising = 0; 
0239          % indicates if the warps applied to the data should be dynamically
0240          % displayed to the user
0241     movie_name = 'undefined';
0242           % the name of the AVI movie to be saved. It will be folloed by
0243           % '.avi' eventually and contain the movie.
0244     format = 'png';
0245          % if data is saved as an image, this paramater defines the format
0246          % as which it shall be saved.
0247     show_warp_targets = 0;
0248          % defines whether or not the target of the warps, i.e. that at
0249          % alignment, should be displayed to the user.
0250     dynamic_precision = 0;
0251        % defines whether or not precision is chosen wisely at each stage
0252     gen_iters = 25;
0253        % generalisability iterations = 25;
0254        % specificity iterations
0255     spec_iters = 25;
0256        % specificity iterations
0257     progress_bar = 1;
0258          % display progress bar or text instead
0259          % 0: Hierarchical Text
0260          % 1: Progress bar
0261          % 2: Quiet
0262          % 3: Console progress bar
0263     model_score_type = 'Default';
0264          % Defined the way in which a model is scored.
0265     force_first_reference = 0;
0266          % depending on its numeric value, this parameter defines the
0267          % reference selection method
0268     user_menu = 0;
0269          % indicates the way in which registration parameters are entered.
0270          % If not via pop-up menu, then probably via top menu in the
0271          % application.
0272     number_of_bins = 50;
0273          % 0: parameters need to be set in file
0274          % 1: menu input
0275     background_red = 1;
0276     background_green = 1;
0277     background_blue = 1;
0278          % the background colourn for the plots
0279     distance_type = 'sum_of_squared_distances';
0280          % The way to meausre distance from mean
0281     show_registration_target = 0;
0282          % indicates if the estimated target of registration should be
0283          % super-imposed in score plots.
0284     produce_statistics = 0;
0285          % should statistics be displayed at the end (boolean)
0286     see_model_before_each_iteration = 0;
0287          % should a combined model be produced and displayed before each
0288          % iteration, this should be set to 1
0289     show_pixels = 0;
0290          % show pixels representation (boolean). Mandatory for pixel movie
0291          % capture
0292     show_score_vs_iteration = 0;
0293          % when enabled, evaluation of the objective function will be
0294          % displayed against the iteration at which it was evaluated.
0295     find_perfect_warp = 1;
0296          % This argument should not be tempered with. It defined the need
0297          % to check for the value and state at conceived registration
0298          % target
0299     show_model_after_optimisation = 0;
0300          % if enabled, then a model produces from the data is displayed
0301          % after registration.
0302     show_specificity_statistics = 0;
0303          % more booleans which are self-explanatory
0304     verbose_score = 0;
0305          % boolean: display score or not after optimisation
0306     show_score_bar = 0;
0307          % The score bar gives a visual indication of the 'goodness' of a
0308          % model construced from the data
0309     average_smooth_enabled = 0;
0310     gaussian_smooth_enabled = 0;
0311          % flags indicating the smoothness to be applied, if any at all
0312     n_points = 5;
0313          % an argument that is passed to the model optimisation function.
0314          % it refers to the number of knotpoint, e.g. for random and scale.
0315     window = 1;
0316          % the size of the window for average smoothing or sigma for
0317          % gaussian smoothing
0318     show_all_score_figure_types = 0;
0319          % defines the need for an abundant number of plots depicting
0320          % quality of model versus iteration
0321     plot_before_and_after = 0;
0322          % will only ever be used if the two above are set
0323     bump_width = 0.2;
0324     bump_width_variation = 0.5;
0325     bump_height = 0.2;
0326     bump_height_variation = 0.3;
0327     bump_position_freedom = 0.9;
0328     initial_diminish_factor = 2;
0329     smoothness_factor = 17;
0330          % set the bump parameters
0331     precision_required = 1e-7;
0332         % the precision required from the optimiser
0333     model_evaluation_method = 1;
0334                       % 1: Determinant
0335                       % 0: MDL
0336                       % 2: Sum of log of Eigenvalues
0337     show_score_using_model = 0;
0338     handle = 'Default';
0339       % give figures a fefault name
0340     n_modes = 20;
0341       % maximal number of modes to be accounted for when evaluating a model
0342     variation_kept = 0.99;
0343       % keep 99%
0344       % number of modes in the model to be accounted for
0345     retain_original_peak = 0;
0346       % determines whether or not to keep the same maximal data height
0347     error_verbosity = 'normal';
0348       % the level of verbosity for errors
0349     initialise_warps_close_to_target = 0;
0350       % determines the requirement to start registration with warps near
0351       % their target (at completion of registration)
0352     shape_weight = 0.06;
0353       % the constant weight Ws for the model constructed. 1 is 100%, 0 is
0354       % 0% and legal range is 0..1
0355     weighting_normalisation_method = 'edge';
0356      % the scheme for Ws weighting selection
0357     run_all_objective_functions = 0;
0358        % do not try running all the objective functions
0359     called_from_gui = 0;
0360        % indicate that GUI was not used in this case
0361     kernel = 'Exponential';
0362        % the kernel type for PDF objective function
0363     save_data = 0;
0364        % Indicated if the data needs to be saved to files after opening or
0365        % generation
0366     peak_type = 'average_position';
0367          % The type of peak to be retained
0368     load_data = 0;
0369        % indicates if data needs to be loaded from file rather than being
0370        % generated.
0371     save_bumps = 0;
0372        % indicated if bumps need to be saved as images (not data)
0373     draw_warp_curve = 0;
0374        % An option for overlapping warps display
0375     draw_curve_by_images = 0;
0376        % The option for warp target maps (images x iterations)
0377     overlap_scores = 0;
0378       % defines the need for scores to be drawn on the same space if
0379       % different objective functions are compared by being run
0380       % sequentially
0381     show_mesh = 0;
0382      % simply indicates if a mesh of the data needs to be displayed
0383     show_target = 0;
0384         % do not show perfectly aligned images estimate
0385     verbose = 'off';
0386          % be verbose or not during general-purpose optimisation
0387     save_stats = 0;
0388        % indicated if statistics need to be saved as HTML
0389     display_shape_model = 0;
0390       % If enabled, shape model is shown
0391     display_intensity_model = 0;
0392       % likewise for intensity model
0393     n_images = 2;
0394              % how many images to try to register in each set
0395     n_iterations = 2;
0396              % how many iterations of the optimisation to run, default was 50.
0397     n_sets = 1;
0398              % how many sets of images to register (in order to get better statistics)
0399     image_width = 50;
0400              % the number of pixels or the lenght of the data vector
0401     figure_handle_number = 1;
0402              % increments applied when figures are saved
0403     plot_alignment = 0; 
0404              % Indicates if a plot for alignment and perturbation should be
0405              % displayed
0406     rotate_surface = 0;
0407              % Indicated if the surface should be rotated horizontally at
0408              % each frame
0409     draw_final_warps = 0;
0410              % indicates if final warps need to be displayed
0411     add_points = 0;
0412              % point stuffing where difference in warp curve is greatest
0413     add_points_cycle = 1;
0414              % the frequency (measured in terms of iteration) for stuffing;
0415 else
0416     called_from_gui = 1;
0417         % indicate that default values were not set and that the function
0418         % was called from REGISTER
0419 end
0420 
0421        
0422 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0423 %% Input / Output              %%
0424 
0425 rand('state',sum(100*clock))
0426           % seeding for random numbers.
0427 
0428 figure_handle_number = 1;
0429           % A number to be incremented every time a figure is saved. It
0430           % manages the number of figures that are generated.
0431 warning off;
0432          % keeps st_pca.m quiet
0433 error_found = 0;
0434          % status flag
0435 current_viewpoint = 45;
0436          % the starting horizontal viewpoint for surfaces
0437          % Optimisation note: initiliase this variables only when surfaces
0438          % are displayed.
0439 status = ['                      '];
0440          % set intitial status message. Messages are aggregated before
0441          % being returned to the GUI and displayed to the user.
0442 if (run_all_objective_functions == 1),
0443          % If all objective functions will run sequentially, then some
0444          % output parameters need to be updated.
0445     colour = 'r';
0446         % set initial colour for overlaid plots
0447     style = ':';
0448         % similarly to style of line in the overlaid plots
0449     progress_bar = 2;
0450 end
0451      % the progress indicators will not work for this routine so disable
0452      % them
0453 
0454 if strcmp(movie_name, 'undefined'),
0455     movie_name = ['movie_', spline_type, '_', placement_type, '_', objective_function,'_width', num2str(image_width), '_images', num2str(n_images), '_iterations', num2str(n_iterations), '.avi'];
0456          % name for the AVI to be saved
0457          % The name will attempt to encompass all construction details
0458          % It will replace the default name given when the function is not
0459          % called from the GUI. If called from the GUI, it usually have a
0460          % user defined name or 'Default.avi'
0461 end
0462 
0463 
0464 frame_counter = 0;
0465          % counter for video frames. Do not fiddle with this variable. It
0466          % keeps track of the current frame number to be captured.
0467          
0468 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0469 %% Others (most are obsolete and are remains of old code.)  %%
0470                        
0471 % white_ctr = 0;
0472          % white counter? Was used for simple bumps in code that is now
0473          % commented out
0474 % ref_shift = 0.2;
0475        % shift in reference image??
0476 % max_shift = 0.2;
0477        % maximum allowed shift?
0478 % step = 0.1;
0479        % used in two lines below only
0480 % los = zeros(n_images,floor((max_shift-ref_shift)/step));
0481 % his = zeros(n_images,floor((max_shift-ref_shift)/step));
0482        % not understood yet
0483 % blurred = 0;
0484        % image blurring has been applied to images (boolean)
0485 % min_error = 0;
0486 % max_error = 1;
0487        % define error range??
0488 % error_step = 0.1;
0489 % n_error_steps = (max_error-min_error) / error_step;
0490 
0491 
0492          
0493 switch model_evaluation_method,
0494     case 0
0495         model_evaluation_method = 'minimum_description_length';
0496     case 1
0497         model_evaluation_method = 'determinant';     
0498     case 2
0499         model_evaluation_method = 'sum_of_log_eigenvalues';
0500     case 3
0501         model_evaluation_method = 'specificity';
0502     case 4
0503         model_evaluation_method = 'generalisability';    
0504     otherwise
0505         error('Unknown model evaluation method. Wrong numeric identifier passed.');
0506 end
0507          % The conversion above can one day be embedded in the GUI to make
0508          % everything more readable. Values are overwritten here too.
0509          
0510 switch progress_bar,
0511     case 0
0512      progress_indicator_type = 'Hierarchical Text';
0513     case 1
0514      progress_indicator_type = 'Progress Bar';
0515     case 2
0516      progress_indicator_type = 'Quiet'; 
0517     case 3
0518      progress_indicator_type = 'Console Progress Bar';
0519     otherwise
0520      error('Callback value for progress indicator is not recognised');
0521 end
0522 
0523          % do type conversions for input to make code below readable
0524          
0525 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0526 %% Begin Set initialisation    %%
0527 
0528 switch user_menu,
0529   case 0
0530     
0531     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0532     %% Optimisation Controllers: Default   %%
0533     
0534     % Original settings from 2003: The values below will be defaults if user menu does not appear
0535     
0536     n_images = 2;
0537              % how many images to try to register in each set
0538     n_iterations = 2;
0539              % how many iterations of the optimisation to run, default was 50.
0540     n_sets = 1;
0541              % how many sets of images to register (in order to get better statistics)
0542     image_width = 50;
0543              % the number of pixels or the length of the data vector
0544              
0545     % A few extra popular choices, added March 2004 to speed up experiments
0546     % construction and customisation
0547     
0548   case 1
0549     [n_sets, n_iterations, n_images, image_width] = get_parameters_from_user('Default');
0550        % Prompts the user for the main registartion parameters using menus
0551 
0552   case 2    
0553     n_images = 10;
0554     n_iterations = 50;
0555     n_sets = 1;
0556     image_width = 50;
0557     
0558   case 3   
0559     n_images = 5;
0560     n_iterations = 5;
0561     n_sets = 1;
0562     image_width = 50; 
0563     
0564   case 4    
0565     n_images = 10;
0566     n_iterations = 10;
0567     n_sets = 1;
0568     image_width = 50;
0569     
0570   case 5   
0571     n_images = 20;
0572     n_iterations = 10;
0573     n_sets = 1;
0574     image_width = 50; 
0575     
0576   case 6
0577     n_images = 20;
0578     n_iterations = 1000;
0579     n_sets = 1;
0580     image_width = 50;
0581     
0582   case 7 
0583     n_images = 10;
0584     n_iterations = 1000;
0585     n_sets = 1;
0586     image_width = 50;
0587     
0588   case 8
0589     n_images = 10;
0590     n_iterations = 10000;
0591     n_sets = 1;
0592     image_width = 50;
0593     
0594   otherwise
0595     error('Unknown parameters option. Internal failure.');
0596     
0597 end 
0598 
0599 % in all the above cases, meaningful names rather than numbers are a
0600 % possible extension which requires time.
0601 
0602 
0603 
0604 
0605 tic;
0606        % start counting time
0607 for this_set = 1:n_sets
0608        % for all sets
0609        % Old: user gets status about set number. Multiple sets generate reliable
0610        %      statistics
0611   tic;
0612        % This is called again (stack of time counters????? Doesn't seem so)
0613        % This will need to be changed (one tic omitted) for the statistics
0614        % to work properly. Only one stopwatch can work at a time.
0615        % Smith: for each set of n_sets images, build a model
0616        % calculate obj fn values from these models, plot the mean with error
0617        % bars given by 1 stddev away
0618   if (strcmp(progress_indicator_type,'Hierarchical Text')),      
0619     disp(['calculating set ', num2str(this_set), ' out of ', num2str(n_sets), ' sets...']);
0620       % display command-line status indicators
0621   elseif (strcmp(progress_indicator_type,'Progress Bar')),
0622       gui_active(1);
0623           % will add an abort button
0624       if (strcmp(getenv('OS'), 'Linux')),
0625          h = progressbar([], 0, 'Preparing for registration', [['Registration Progress on '], [get_host]]);
0626       else
0627          h = progressbar([], 0, 'Preparing for registration', 'Registration Progress');
0628       end
0629         % initiate progress bar if necessary
0630       h = progressbar( h, 0, 'Fetching Data');
0631   end    
0632 
0633   % initiate score bar if necessary
0634 
0635 
0636   if (show_score_bar == 1),
0637      score_bar = waitbar(0, 'Score');
0638   end
0639      
0640 
0641   switch load_data,
0642       
0643    case 0
0644       [imagelist, images, points] = make_1d_bump(n_images, image_width...
0645           , bump_width, bump_width_variation, bump_height, bump_height_variation...
0646           , bump_position_freedom, initial_diminish_factor, smoothness_factor...
0647           , save_bumps, format);
0648        % create a set of images given several parameters
0649        
0650       points = -1 + 2 * (points - 1) / (image_width - 1);
0651         % Smith: normalise points from -1 to 1
0652              
0653       if (called_from_gui == 1 & save_bumps == 1),
0654         status = [[status];['(**) Images saved     '];['                      ']];
0655       end
0656         % notify user using the status bar that data has been saved.
0657  
0658        %% keep = 0.999999; NOTE (RSS Feb. 2004): became obsolete and
0659        %% replaced by <variation kept> which relates to PCA
0660        
0661        % Nov. 2003: Is this precision required? Perhaps in optimisation.
0662        
0663    case 1
0664       output = open([[handle],['.mat']]);
0665       points = output.points;
0666       images = output.images;
0667                   % copy the data structures from the file
0668                   
0669    case 2
0670       output = open('Default.mat');
0671       points = output.points;
0672       images = output.images; 
0673       
0674    case 3
0675       output = open('flat.mat');
0676       points = output.points;
0677       images = output.images; 
0678       
0679    case 4
0680       output = open('plain.mat');
0681       points = output.points;
0682       images = output.images;  
0683       
0684    case 5
0685       output = open('plainran.mat');
0686       points = output.points;
0687       images = output.images;  
0688       
0689    case 6
0690       [imagelist, images, points] = make_elliptic_bumps(n_images, image_width, save_bumps, format, bump_height, bump_position_freedom, bump_height_variation, bump_width_variation);
0691        % create a set of images given several parameters
0692       points = -1 + 2 * (points - 1) / (image_width - 1);
0693              
0694       if (called_from_gui == 1 & save_bumps == 1),
0695         status = [[status];['(**) Images saved     '];['                      ']];
0696       end
0697       
0698    case 7
0699       [imagelist, images, points, dummy1, dummy2] = make_1d_images(n_images, image_width, 0.2);
0700        % 0.2. is the value used by Kate for <white width>
0701       if (save_bumps == 1),
0702           msgbox('Data save is not supported with simple bumps.');
0703       end
0704        % create a set of images given several parameters
0705       points = -1 + 2 * (points - 1) / (image_width - 1);
0706              
0707       if (called_from_gui == 1 & save_bumps == 1),
0708         status = [[status];['(**) Images not saved '];['                      ']];
0709       end    
0710       
0711    otherwise
0712       error('Invalid value passed from open submenu. Internal error.');
0713       
0714   end     
0715   
0716   
0717   if (average_smooth_enabled),
0718       images = average_smooth(images, window);
0719   end  
0720   
0721   if (gaussian_smooth_enabled),
0722       images = gaussian_smooth(images, window); 
0723   end
0724        % smooth all images - used to be disabled and set <blurred> to 1
0725        
0726 
0727   ref_points_vec = points(:,1);
0728   ref_image_vec = images(:,1);
0729        % set first image generated to be reference when calculating a
0730        % perfect warp and also as reference in pair-wise objective
0731        % functions. This has been made more generalised later on.
0732   % ref_hi = his(1);
0733   % ref_lo = los(1);
0734        % and get its upper and lower boundaries. his and lois were
0735        % generated and returned as output for simple bumps. They were
0736        % used to find the perfect warps and defined points of importance in
0737        % the input data.
0738   
0739   % Now choose the image which will not be changed. Note a possible
0740 % optimisation: only if a reference is forced, then the code below should
0741 % be run to choose the reference. The section below that is to do with
0742 % alignment means however that for target to be shown an immutable image
0743 % must always be chosen
0744 
0745       switch force_first_reference,
0746        case 0        
0747            immutable_image = 2;
0748                 % see comments later about the choice of first or second
0749                 % image as an immutable reference. These options will only
0750                 % have an effect if a reference is forced (not as in the
0751                 % default reference in MSD-based algorithms)
0752                 
0753                 % OPTIMISATION NOTE: score can also be set to 0 in all
0754                 % cases without breaking registration. It would then need
0755                 % not be repeated in all the cases below.
0756                 
0757        case 1
0758            immutable_image = 1;
0759            score = 0;
0760                 % the score needs to be initilialised because no previous
0761                 % scoring will be avaliable
0762                 
0763        case 2
0764           % closest data to average (mean)
0765            immutable_image = find_image_closest_to_mean(images, image_width, n_images, 'closest', distance_type);
0766            score = 0;
0767            
0768        case 3
0769               % random reference
0770            immutable_image = ceil(rand * n_images);
0771            score = 0;
0772            
0773        case 4
0774           % closest data to average (mean)
0775            immutable_image = find_image_closest_to_mean(images, image_width, n_images, 'farthest', distance_type);
0776            score = 0;
0777            
0778        case 5
0779            % highest peak
0780            immutable_image = find_image_by_height(images, image_width, n_images, 'highest');
0781            score = 0;
0782            
0783        case 6
0784            % lowest peak
0785            immutable_image = find_image_by_height(images, image_width, n_images, 'lowest');
0786            score = 0;
0787            
0788        otherwise
0789           error('Reference error. Wrong values passed from callback.');
0790           
0791       end  
0792 
0793       % NOTE PROBLEM: If teh refrence chosen is the first then evaluation at
0794       % the start is 0 and the optimisation plot is not suitable
0795        
0796        
0797        
0798        
0799   %%%%% Do warp analysis here
0800   
0801   % The check for perfect warps does not work for corrective offsets,
0802   % non-bump data, and data that is significantly smoothed. In all of these
0803   % cases, the flag below must be set to 0, or else piece wise linear will
0804   % not work.
0805   check_perfect_warps = find_perfect_warp;
0806           % currently always enabled, but computationally expensive and in
0807           % some cases not necessary. A test for the value of this part to
0808           % be implemented in the future.
0809       if (check_perfect_warps == 1),   
0810           if (strcmp(progress_indicator_type,'Progress Bar')),
0811              h = progressbar( h, 0, 'Calculating Correct Warps');
0812           end
0813         % (Old - November 2003 - Smith)
0814           % correctly_warped_points = points;
0815           % correctly_warped_images = images;
0816           for i = 1:n_images,
0817               % Smith: now bung in the right deformations to get them back
0818               % NOTE: THE VARIABLE ref_image_vec BELOW ASSUMED A FIRST
0819               % IMAGE REFERENCE although this is not the case for the <force
0820               % reference> option which makes it the second (implicitly).
0821               % Somthing richer than a second image might be used in
0822               % future.
0823               
0824               %possible optimisation: do not interpolate for the first
0825               %image which is usually a reference and is not affected
0826               %anyhow
0827               %images(:,i)
0828             correctly_warped_points(i,:) = find_target_warp(images(:,immutable_image), images(:,i), image_width);
0829               % correctly_warped_points(i,:)
0830               % get the warp that appears most desirable
0831             correctly_warped_points(i,:) = -1 + 2 * (correctly_warped_points(i,:) - 1) / (image_width - 1);
0832               % normalise points so that they lie in the range x where -1 < x < 1
0833             correctly_warped_points_vec = correctly_warped_points(i,:);              
0834             image_vec = images(:,i);
0835               % RSS: copy original image and copy piece-wise linear warp
0836               % Smith: resample at the warped points - point on reg. grid has value of point at same index in warpy grid
0837               % need to interpolate
0838               % size(ref_points_vec)
0839               % size(image_vec)
0840               % size(correctly_warped_points_vec')
0841             correctly_warped_images(i,:) = interp1(ref_points_vec, image_vec, correctly_warped_points_vec);
0842             % correct_model = build_model(correctly_warped_images,correctly_warped_points,1,'','edge',0);
0843             % converging_score(i) =
0844             % measure_model(correct_model.variances,50);
0845             % correctly_warped_images(i,:)
0846           end
0847           % figure;
0848           % plot(log(converging_score));
0849           % title('Log of score as model converges');
0850       end
0851       % figure;
0852       % (test curve) plot(correctly_warped_points(2,:))
0853       
0854   
0855 
0856       if (show_warp_targets == 1),
0857          warp_handle = figure('Name','On the left: piece-wise linear optimal registration warps; On the right: results of warp');
0858          % title('On the left: piece-wise linear optimal registration warps; On the right: results of warp');
0859          save_handle = figure(warp_handle);
0860          hold on;
0861          subplot(n_images, 2, 1);
0862          for i=1:n_images,
0863              subplot(n_images, 2, i * 2 - 1);
0864              hold on;
0865              axis off;
0866              plot(correctly_warped_points(i,:));
0867              axis([-1 image_width -1 1]);
0868              hold off;
0869          end
0870          for i=1:n_images,
0871              subplot(n_images, 2, i * 2);
0872              hold on;
0873              axis off;
0874              plot(correctly_warped_images(i,:));
0875              axis([0 image_width 0 1]);
0876              hold off;
0877          end
0878          hold off;
0879          if (save_in_log == 1),
0880            saveas(save_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
0881            figure_handle_number = figure_handle_number + 1;
0882          end
0883       end
0884       
0885   
0886     
0887   if (initialise_warps_close_to_target == 1),
0888       % offset_extent = 5;
0889       % perturbation_method = 'CPS warp';
0890       
0891     if (plot_alignment == 1),
0892         % if a plot is required to show this step of perturbation after
0893         % alignment the right-hand side solumn will be occupied by the
0894         % plots
0895       subplot_alignment_fig = figure('Name',[[handle],[': Left: Aligned images ; Right (if existent): Aligned images after perturbation']]);
0896       alignment_handle=figure(subplot_alignment_fig);
0897       for i=1:n_images,
0898        subplot(n_images, 2, 2 * i - 1);
0899        hold on;
0900        plot(correctly_warped_images(i,:));
0901        % title(['Unwarped image #',num2str(i)]);
0902        axis([0 image_width 0 1]);
0903        grid on;
0904        hold off;
0905       end
0906     end
0907       
0908       if (offset_extent == 0),
0909           % if no perturbation is required then the values can be simply
0910           % copied.
0911          for i = 1:n_images,
0912            % copy registration target
0913           points(:,i) = correctly_warped_points(i,:)';
0914            % images(:,i) = correctly_warped_images(i,:)';
0915          end
0916       else
0917           % if some noise is required, then add it in the form of Gaussian
0918           % noise or CPS warps to points and interpolate original images
0919           % accordingly. To be controlled from GUI.
0920            for i = 1:n_images,
0921                 tweaked_points(i,:) = apply_perturbation_to_points(correctly_warped_points(i,:), offset_extent, perturbation_method);
0922         %         figure
0923         %         hold on
0924         %         plot(tweaked_points(i,:));
0925         %         plot(correctly_warped_points(i,:));
0926         %         hold off
0927                    % apply noise to the correct points
0928                 % images(:,i) = interp1(ref_points_vec, images(:,i), tweaked_points(i,:))';
0929                    % interpolation to get current images
0930                 points = tweaked_points';
0931                    % copy the points found to become current
0932                    % inconsistent column and row form is due to initial
0933                    % implementation. Optimally (when our task becomes optimising performance),
0934                    % all transposes should be ridden of.
0935            end
0936            if (plot_alignment == 1),
0937               for i = 1:n_images,
0938                 perturbed_aligned_images(:,i) = interp1(ref_points_vec, images(:,i), points(:,i));
0939               end    
0940                 % inerpolate to show some plots at this earlier stage
0941               alignment_handle = figure(subplot_alignment_fig);
0942               for i = 1:n_images,
0943                    subplot(n_images, 2, 2 * i);
0944                    hold on;
0945                    plot(perturbed_aligned_images(:,i));
0946                    % title(['Unwarped image #',num2str(i)]);
0947                    axis([0 image_width 0 1]);
0948                    grid on;
0949                    hold off;
0950                end
0951              if (save_in_log == 1)
0952                saveas(alignment_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
0953                figure_handle_number = figure_handle_number + 1;
0954              end
0955            end
0956       end
0957       
0958   end
0959   
0960   
0961   if (save_data == 1),
0962       % note that data will be saved with modified warp curves if option is
0963       % set. This is intentional. Data in this cotext is vectors, not
0964       % images.
0965       save([[handle],['.mat']], 'images','points');
0966       if (called_from_gui == 1),
0967         status = [[status];['(**) Data saved       '];['                      ']];
0968       end
0969   elseif (save_data == 2),
0970       % if data needs to be saved as default
0971       save([[handle],['.mat']], 'images','points');
0972       save(['Default.mat'], 'images','points');
0973       if (called_from_gui == 1),
0974         status = [[status];['(**) Data saved (Def) '];['                      ']];
0975       end
0976   end
0977 
0978 
0979   
0980   if (retain_original_peak),  
0981        for current_image = 1:n_images,
0982            maximum_height(current_image) = max(images(:, current_image));
0983               % find out the maximum height of each data instance
0984            
0985            for j=1:image_width,
0986                switch peak_type,
0987                    case 'initial_position'
0988                          if (images(j,current_image) == maximum_height(current_image)),
0989                            maximum_position(current_image) = j;
0990                             % find the position where height is maximal. Currently
0991                             % this search is very inefficient.
0992                          end
0993                          
0994                    case 'average_position'
0995                          if (images(j,current_image) == maximum_height(current_image)),
0996                            maximum_position(current_image) = j;
0997                          end
0998                          % same as above
0999                end
1000            end
1001        end
1002  
1003   
1004        switch peak_type,
1005          case 'average_position'
1006            average_position = floor(mean(maximum_position));
1007            for current_image = 1:n_images,
1008              maximum_position(current_image) = average_position;
1009            end 
1010        end
1011    end
1012            
1013            
1014            
1015            
1016 % plot3(1:2,1:2,images(1,:))
1017 %
1018 %
1019 %       s = surface(L, ...
1020 %         'EdgeColor','none', ...
1021 %         'FaceColor',[0.8 0.2 0.2], ...
1022 %         'FaceLighting','phong', ...
1023 %         'AmbientStrength',0.3, ...
1024 %         'DiffuseStrength',0.6, ...
1025 %         'Clipping','off',...
1026 %         'BackFaceLighting','lit', ...
1027 %         'SpecularStrength',1.1, ...
1028 %         'SpecularColorReflectance',1, ...
1029 %         'SpecularExponent',7);
1030 %     l1 = light('Position',[40 100 20], ...
1031 %         'Style','local', ...
1032 %         'Color',[0 0.7 0.7]);
1033 %     l2 = light('Position',[.5 -1 .4], ...
1034 %         'Color',[1 1 0]);
1035 
1036 % modify the progress bar if required
1037 
1038   if (strcmp(progress_indicator_type,'Progress Bar')),
1039            h = progressbar( h, 0, 'Pre-processing Data');
1040   end
1041   
1042 % depending on whether intensity offset was applied, record the number of
1043 % columns required for the plots
1044        
1045   if (apply_intensity_offset),
1046       vertical_pane_location = 3;
1047                        % Need is for only 3 columns in plot
1048   else
1049       vertical_pane_location = 2;
1050                        % Need is for only 2 columns in plot
1051   end
1052   % Note: intensity offsets are considered a possible way of solving the
1053   % problem, but not natural form for our input data. Hence, saving of data
1054   % precedes this offset.
1055 
1056 
1057 % prepare a figure if a plot of all warp curves needs to be rendered
1058        
1059   if (draw_warp_curve == 1),
1060        overlap_curve_fig = figure('Name','Overlapping Warp Curves');
1061   end
1062   
1063   % if pixel representation is required, prepare the figure
1064   
1065   if (show_pixels == 1),
1066     pixels_figure = figure('Name','Pixel Representation');
1067   end
1068   
1069   % and likewise for meshes/surfaces
1070 
1071   if (show_mesh == 1),
1072     mesh_subfigure = figure('Name','1-D Data Mesh');
1073   end 
1074   
1075   % if a matrix of warps is required, preapre the figure and draw the
1076   % initial warps in its left side
1077   
1078   if (draw_curve_by_images == 1),
1079     curve_subplot_fig = figure('Name',[['Warp Curves for all data - objective function: '],[objective_function]]);
1080     curves = figure(curve_subplot_fig);
1081     % title('Warp Curves for all data');
1082     for i=1:n_images,
1083        subplot(n_images, n_iterations + 1, i * (n_iterations + 1) - n_iterations);
1084        hold on;
1085        axis off;
1086        plot(points(:,i));
1087        plot(correctly_warped_points(i,:),'r-')
1088        if (i == 1 ),
1089          legend('Initial warp', 'Target warp'); 
1090             % only display the legend once at the top
1091        end
1092        axis([-1 image_width -1 1]);
1093        hold off;
1094     end
1095   end
1096   
1097 
1098 %%             END OF INITIALISATION               %%
1099 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1100 
1101 
1102 
1103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1104 %% Warping images, then modelling warped images    %%
1105 %% attempt to register images by warping.          %%
1106 
1107 
1108 
1109   if (plot_before_and_after),
1110     for i=1:n_images,
1111        images_before_offset(:,i) = images(:,i);
1112        % copy these original values in case 3 columns are necessary for plot representation of the data.
1113     end
1114   end
1115   
1116   % draw the original data at the leftmost column
1117   
1118   if (do_plot),
1119     subplot_fig = figure('Name',[[handle],[': Left: Original images ; Right: Images subjected to warping']]);
1120     H=figure(subplot_fig);
1121     for i=1:n_images,
1122        subplot(n_images, vertical_pane_location, (vertical_pane_location * i) - vertical_pane_location + 1);
1123        hold on;
1124        plot(images(:,i));
1125        whitebg(H, [background_red, background_green, background_blue]); 
1126        % title(['Unwarped image #',num2str(i)]);
1127        axis([0 image_width 0 1]);
1128        grid on;
1129        hold off;
1130     end
1131   else
1132     subplot_fig = 0;
1133       % pass a void figure to the optimisation functions
1134   end
1135   
1136   % draw the data after intensity offset has been applied (when applicable)
1137   
1138   if (apply_intensity_offset),
1139       for index = 1:n_images,
1140          peak (index) = max(images(:,index));
1141       end
1142       peak_reference = max (peak);  
1143       for index = 1:n_images,
1144         peak_difference = peak_reference - max(images(:,index));
1145         images(:,index) = images(:,index) + peak_difference;
1146         if (do_plot),
1147             subplot(n_images, vertical_pane_location, (vertical_pane_location * index) - vertical_pane_location + 2);
1148             hold on;
1149             plot(images(:,index));
1150             % title(['Unwarped image #',num2str(i)]);
1151             axis([0 image_width 0 1]);
1152             grid on;
1153             hold off;
1154         end
1155       end        
1156   end 
1157 %% NOTE: Change here in traversal above (RSS, December 2003)
1158   
1159     still_to_go = 1;
1160           % condition for while loop. This will enable the program to
1161           % continually try all objective functions.
1162     
1163     
1164     while (still_to_go == 1),
1165     
1166     if (run_all_objective_functions == 0)
1167         still_to_go = 0;
1168            % break out of while loop after just one iteration
1169     else
1170            % If new objective functions are added, an entry must be added
1171            % here to include it in the more comprehensive cycle.
1172            
1173            % try all objective functions sequentially
1174        switch objective_function,    
1175            case 'all'
1176            % meaning that it is in its first iteration
1177              objective_function = 'model_opt_together';
1178            case 'model_opt_together'
1179              objective_function = 'model_opt_separate';    
1180            case 'model_opt_separate'
1181              objective_function = 'msd_opt_together';
1182            case 'msd_opt_together'
1183              objective_function = 'msd_opt_separate';
1184            case 'msd_opt_separate'
1185              objective_function = 'mi_opt_together';
1186            case 'mi_opt_together'
1187              objective_function = 'mi_opt_separate';
1188            case 'mi_opt_separate'
1189              objective_function = 'nmi_opt_together';
1190            case 'nmi_opt_together'
1191              objective_function = 'nmi_opt_separate';
1192            case 'nmi_opt_separate'          
1193              objective_function = 'mixed_msd_model';
1194            case 'mixed_msd_model'
1195              objective_function = 'model_post_msd';
1196            case 'model_post_msd'        
1197              objective_function = 'tfc_eccv';
1198            case 'tfc_eccv'
1199              objective_function = 'wavelet_complexity';
1200            case 'wavelet_complexity'
1201              objective_function = 'model_pairs_together';
1202            case 'model_pairs_together'
1203              objective_function = 'model_pairs_separate';
1204            case 'model_pairs_separate'
1205              objective_function = 'void'; 
1206              still_to_go = 0;
1207               % signal to break out of the while loop because the test of
1208               % all objective functions is finished.
1209            otherwise
1210              error('Problem with <Run All Objective Functions> option');
1211        end
1212     end   
1213       % warped_images = images
1214       % size(images)
1215       for i=1:n_images,
1216         warped_images(:,i) = interp1(ref_points_vec, images(:,i), points(:,i)')';
1217       end
1218       % size(warped_images)
1219       warped_points = points;
1220           % set this initial state to plot the examples before change
1221   
1222           
1223        
1224       
1225       
1226   
1227       if ((strcmp(progress_indicator_type,'Console Progress Bar')) | (strcmp(progress_indicator_type,'Progress Bar')))
1228          progress_bar_position = 0;
1229          time_for_this_iteration = 10;
1230          total_images = n_images * n_iterations;
1231          step_size = 50 / total_images;
1232             % these variables are used for some progress statistics
1233       end   
1234         % initialise progress bar
1235         
1236 
1237 
1238         
1239       for n = 1:n_iterations,
1240         
1241         if (dynamic_precision == 1),
1242            iterations_ratio = n_iterations / n;
1243              % get iterations ratio. In future, looking at the optimisation
1244              % state and behaviour, for example whether it is stuck,
1245              % should give better results.
1246            precision_required = set_precision(iterations_ratio, precision_automation_type);
1247         end
1248         
1249         if (show_score_bar == 1),
1250             % subplot('position',[0,0,1,0.1]);
1251             Epsilon = 0.0001;
1252             c_model = build_model(warped_images, warped_points, variation_kept, '', weighting_normalisation_method, shape_weight);
1253             % waitbar((measure_model(c_model.variances, n_modes, 2, c_model) - log(Epsilon.^n_modes)) / log(Epsilon.^n_modes), score_bar);
1254             waitbar((measure_model(c_model.variances, n_modes, 'sum_of_log_eigenvalues', c_model) ...
1255                / log(Epsilon.^n_modes) - 0.8) * 5, score_bar);
1256             % (measure_model(c_model.variances, n_modes, 2, c_model) - log(Epsilon.^n_modes)) / log(Epsilon.^n_modes) * 100
1257             % axis off;
1258             % subplot(1,1,1);
1259         end            
1260           
1261           
1262           
1263         if (show_mesh == 1),
1264           meshview = figure(mesh_subfigure);
1265           whitebg(meshview, [background_red, background_green, background_blue]);
1266                % set viewpoint
1267           [X,Y] = meshgrid(1:n_images * 5 + 2, -image_width / 2 - 1:1:image_width / 2);
1268                % two zero buffers and 5 rows of each bump
1269           for i = 1:n_images,
1270             Z(2:image_width + 1,1 + i * 5) = warped_images(1:image_width,i);
1271             Z(2:image_width + 1,1 + i * 5 - 1) = warped_images(1:image_width,i);
1272             Z(2:image_width + 1,1 + i * 5 - 2) = warped_images(1:image_width,i);
1273             Z(2:image_width + 1,1 + i * 5 - 3) = warped_images(1:image_width,i);
1274             Z(2:image_width + 1,1 + i * 5 - 4) = warped_images(1:image_width,i);
1275             % copy the data values
1276           end
1277           Z(:,1) = zeros(1:image_width, 1);
1278                 % set bottom buffer
1279           Z(:, n_images * 5 + 2) = zeros(1:image_width, 1);
1280                 % and upper buffer to close mesh
1281           Z(1,:) = zeros(1, 1:n_images * 5);
1282                 % set left buffer
1283           Z(image_width+2,:) = zeros(1, 1:n_images * 5);
1284                 % and right buffer to close mesh
1285           switch mesh_type,  
1286            case 'Surface Red'  
1287             surf(X,Y,Z,'FaceColor', 'red', 'EdgeColor', 'none');
1288             camlight left;
1289             lighting phong;
1290             axis([0 n_images*5 -image_width/2 image_width/2 0 1]);     
1291            case 'Surface Blue'     
1292             surf(X,Y,Z,'FaceColor', 'blue', 'EdgeColor', 'none');
1293             camlight left;
1294             lighting phong;
1295             axis([0 n_images*5 -image_width/2 image_width/2 0 1]);     
1296            case 'Polygons'   
1297             mesh(X,Y,Z);
1298             axis([0 n_images*5 -image_width/2 image_width/2 0 1]);     
1299            case 'Illuminated Surface'
1300             if (rotate_surface),
1301               current_viewpoint = current_viewpoint + 1;
1302                % shift viewpoint slightly
1303             end
1304             view(current_viewpoint, 37.5); 
1305             surface(X,Y,Z, ...
1306               'EdgeColor','none', ...
1307               'FaceColor',[0.5 0.5 0.5], ...
1308               'FaceLighting','phong', ...
1309               'AmbientStrength',0.3, ...
1310               'DiffuseStrength',0.6, ... 
1311               'Clipping','off',...
1312               'BackFaceLighting','lit', ...
1313               'SpecularStrength',1.1, ...
1314               'SpecularColorReflectance',1, ...
1315               'SpecularExponent',7);
1316             if (n == 1), % set lights only at the start
1317                 l1 = light('Position',[40 100 20], ...
1318                   'Style','local', ...
1319                   'Color',[0 0.5 0.7]);
1320                 l2 = light('Position',[.5 -1 .4], ...
1321                   'Color',[1 1 0]);
1322             end
1323            otherwise
1324             error ('Mesh type not recognised. Internal error.');
1325           end
1326           
1327           axis off; 
1328           if (movie_type == 2 & save_movie == 1),
1329                   frame_counter = frame_counter + 1;
1330                   M(frame_counter) = getframe(meshview);
1331           end
1332           if (save_in_log == 1 & n == n_iterations), 
1333              saveas(meshview,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1334              figure_handle_number = figure_handle_number + 1;
1335           end
1336          
1337         end  
1338           
1339         if (strcmp(progress_indicator_type,'Hierarchical Text'))  
1340            % iterations aim to get good statistical results by averaging
1341           disp(['    +-Running iteration number ',num2str(n), ' out of ',num2str(n_iterations), ' iterations...']);
1342         end
1343         
1344            % Smith: correct registration
1345            % attempt to register
1346                
1347         if (see_model_before_each_iteration),  
1348           w_c_model = build_model(warped_images, warped_points, variation_kept,'Optimised warp', weighting_normalisation_method, shape_weight);
1349                   % create an appearance model...
1350           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)]]);
1351                   % and show the model before each iteration while saving
1352                   % the figure handle
1353           if (movie_type == 3 & save_movie == 1),
1354                 % if moview needs to be saved and it is a model
1355                   frame_counter = frame_counter + 1;
1356                   M(frame_counter) = getframe(model_view);
1357           end        
1358         end          
1359         for i = 1:n_images,
1360            % start iteration time counter
1361          if (strcmp(progress_indicator_type,'Hierarchical Text')),
1362            disp(['        +-Warping image ',num2str(i),' out of ',num2str(n_images), ' images in total...']);
1363          end
1364          
1365          this_step = ((n-1) * n_images) + i;
1366                     
1367          if (strcmp(progress_indicator_type,'Progress Bar')),
1368            tic;
1369            progress_bar_position = progress_bar_position + step_size;
1370            steps_remaining = (n_images * n_iterations) - this_step;
1371            percentage_progress = ceil(this_step / (n_images * n_iterations) * 100);
1372            minutes = floor(time_for_this_iteration * steps_remaining / 60);
1373            seconds = rem(floor(time_for_this_iteration *  steps_remaining), 60);
1374            max_count   = 50;
1375                % The console-based GUI is 50 units wide
1376            if (seconds > 9),
1377                h = progressbar( h, step_size/50,[num2str(percentage_progress),'% Completed ---- Estimated remaining time: ', num2str(minutes), ':', num2str(seconds)]);
1378            else
1379                h = progressbar( h, step_size/50,[num2str(percentage_progress),'% Completed ---- Estimated remaining time: ', num2str(minutes), ':0', num2str(seconds)]);
1380            end    
1381            if ~gui_active,
1382                     error_found = 1;
1383                     break; 
1384                % If aborted was pressed
1385            end
1386               % use the third-party GUI-like progress bar (added February
1387               % 2004)
1388          end
1389          if (strcmp(progress_indicator_type,'Console Progress Bar')),
1390                                % if console progress bar
1391            tic;
1392            progress_bar_position = progress_bar_position + step_size;
1393            clc;
1394            disp(['|=================================================|']);
1395            progress_string='|';       
1396            for counter = 1:floor(progress_bar_position),
1397                progress_string = [progress_string, '#'];
1398            end
1399            disp(progress_string);
1400            disp(['|================= ',num2str(floor(progress_bar_position * 2)),'% completed =================|']);
1401                           % display progress per cent
1402            steps_remaining = (n_images * n_iterations) - this_step;
1403            minutes = floor(time_for_this_iteration * steps_remaining / 60);
1404            seconds = rem(floor(time_for_this_iteration *  steps_remaining), 60);
1405            disp(' ');
1406            if (seconds > 9),
1407              disp(['            Estimated remaining time: ', num2str(minutes), ':', num2str(seconds)]);
1408                           % show time indicators
1409            else
1410              disp(['            Estimated remaining time: ', num2str(minutes), ':0', num2str(seconds)]);
1411            end
1412          end
1413          
1414          image_vec = warped_images(:,i);
1415          points_vec = warped_points(:,i);     
1416          
1417          if (show_pixels == 1),   
1418              appended_vector=[];
1419              for current_image=1:n_images,
1420                appended_vector=[[appended_vector],[warped_images(:,current_image),[warped_images(:,current_image)]],[warped_images(:,current_image)]];
1421                    % append bump data 3 times for width which is visible
1422              end
1423            
1424                  %% appended_vector=appended_vector * 256
1425                  %% use if scale of 0..255 used
1426     
1427              G = figure(pixels_figure);
1428              hold on;
1429               % title('Bumps as pixels');
1430              imshow(appended_vector',[0,1]);
1431                        % show as an image
1432              if (movie_type == 0 & save_movie == 1),
1433                                    % Note: already nested in a show_pixels
1434                                    % 'if' statement
1435                   frame_counter = frame_counter + 1;
1436                   M(frame_counter) = getframe(G);
1437              end
1438              hold off;
1439          end
1440          
1441          if (movie_type == 1 & save_movie == 1 & do_plot == 1),
1442                    frame_counter = frame_counter + 1;
1443                    M(frame_counter) = getframe(H);
1444          end
1445                    % capture pixels frame
1446          
1447          
1448          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
1449                        % Somthing richer than a second image might be used
1450                        % in future e.g. mean image
1451              switch objective_function,     
1452                  case 'model_opt_together'
1453                      
1454 %                      if (i==1 & n==1),
1455 %                          score=0;
1456 %                      end
1457 %
1458 %                      last_score = score;
1459 %
1460 %                      % disp(score);
1461 %                   if (i==3)
1462 %                       figure(101);
1463 %                       plot(image_vec);
1464 %                   end
1465                      [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);
1466                      
1467                      %if (n==1),
1468                      %  [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);
1469                      %else
1470                      %   [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);
1471                      % end
1472                      
1473                      
1474 %                      if (last_score < score),
1475 %                          disp(num2str((n - 1) * i + i))
1476 %                      end
1477 
1478                  case 'model_pairs_together'    
1479                      [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);    
1480                  case 'model_pairs_separate'  
1481                      [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);
1482                  case 'model_opt_separate'
1483                      [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);
1484                  case 'msd_opt_together'         
1485                      [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);
1486                  case 'msd_opt_separate'
1487                      [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);
1488                  case 'mi_opt_together'       
1489                      [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);
1490                  case 'mi_opt_separate'
1491                      [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);
1492                  case 'nmi_opt_together'        
1493                      [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);
1494                  case 'nmi_opt_separate'
1495                      [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);
1496                  case 'mixed_msd_model'
1497                         if (((n / n_iterations < 2 / 10) & (n / n_iterations >= 1 / 10)) |...
1498                             ((n / n_iterations < 4 / 10) & (n / n_iterations >= 3 / 10)) |...
1499                             ((n / n_iterations < 6 / 10) & (n / n_iterations >= 5 / 10)) |...
1500                             ((n / n_iterations < 8 / 10) & (n / n_iterations >= 7 / 10)) |...
1501                             ((n / n_iterations < 10 / 10) & (n / n_iterations >= 9 / 10))),
1502                           [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);
1503                         else
1504                           [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);
1505                         end                    
1506                  case 'model_post_msd'
1507                         if (n / n_iterations < 1 / 2),
1508                             % for the first half
1509                           [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);
1510                         else
1511                             % and for the second half
1512                           [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);
1513                         end
1514                  case 'tfc_eccv'
1515                      % msgbox('TFC ECCV 2004: Under Construction.');
1516                      pdf_evalution_method = 'Default';
1517                      [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);         
1518                  case 'wavelet_complexity'
1519                      % msgbox('Wavelets: Under Construction.');
1520                      wavelet_evalution_method = 'Default';
1521                      [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);                      
1522                  case 'void'
1523                     % This is used for repeated experiments with overlaid
1524                     % plots. It is not well tested and score may need to be
1525                     % assigned to a value here
1526                      break;
1527                      % It is not obsolutely clear id the above is necessary
1528                      % or works
1529                  otherwise
1530                      error(['Unknown objective function: ', objective_function]);
1531              end
1532          else
1533             % score = 0;
1534             warped_image = warped_images(:,i);
1535             warped_point = warped_points(:,i); 
1536                   % copy some values to make up for the skip
1537             if (do_plot == 1),
1538               figure(subplot_fig);
1539               subplot(n_images, vertical_pane_location , vertical_pane_location * i, 'replace');
1540               plot(warped_image);
1541               axis([0 image_width 0 1]);
1542               drawnow;
1543             end
1544                   % draw the plot
1545          end
1546 
1547 
1548          warped_images(:,i) = warped_image;
1549            
1550 
1551          if (retain_original_peak), 
1552            warped_images(maximum_position(i), i) = maximum_height(i);
1553               % always retain the maximum height
1554          end     
1555               
1556               
1557               
1558               
1559          warped_points(:,i) = warped_point;
1560          
1561          if (draw_warp_curve == 1),
1562            figure(overlap_curve_fig);
1563            hold on;
1564            plot(warped_points(:,i));
1565            hold off;
1566            if (save_in_log == 1 & n == n_iterations & i == n_images), 
1567                saveas(overlap_curve_fig,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1568                figure_handle_number = figure_handle_number + 1;
1569            end
1570          end     
1571          
1572          if (draw_curve_by_images == 1),
1573              figure(curves);
1574              subplot(n_images, n_iterations + 1, i * (n_iterations + 1) - (n_iterations - n));
1575              hold on;
1576              axis off;
1577              plot(warped_points(:,i));
1578              axis([-1 image_width -1 1]);
1579              hold off;
1580              if (save_in_log == 1 & n == n_iterations & i == n_images), 
1581                saveas(curves,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1582                figure_handle_number = figure_handle_number + 1;
1583              end
1584          end
1585          if (show_pixels == 1 & save_in_log == 1 & n == n_iterations & i == n_images), 
1586                saveas(G,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1587                figure_handle_number = figure_handle_number + 1;
1588          end
1589          if (do_plot == 1 & save_in_log == 1 & n == n_iterations & i == n_images), 
1590                saveas(H,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1591                figure_handle_number = figure_handle_number + 1;
1592          end
1593          
1594          % figure(subplot_fig),title(['Warped images, after image ',num2str(i),' iteration ',num2str(n)]);
1595          % dummy = waitforbuttonpress;
1596          if (show_score_using_model == 1),           
1597              % if evaluation is through model then always overwrite score
1598             c_model = build_model(warped_images, warped_points, variation_kept, '', weighting_normalisation_method, shape_weight);
1599                 % first build a model
1600             if (strcmp(model_score_type, 'Default')),
1601                 score = measure_model(c_model.variances, n_modes, model_evaluation_method, c_model);
1602             elseif (strcmp(model_score_type, 'Specificity')),
1603                 if (mod(n_images * (n - 1) + i, cycle) == 0),
1604                     % disp('ok')
1605                   [score, std_spec] = find_specificity(c_model, images, spec_iters, ref_points_vec);
1606                   standard_deviation(n_images * (n - 1) + i) = std_spec;
1607                 else
1608                     % set to -1 (or copy from predecessor - the original
1609                     % intension
1610                   score = -1; 
1611                 end
1612             elseif (strcmp(model_score_type, 'Generalisability')), 
1613                if (mod(n_images * (n - 1) + i, cycle) == 0),
1614                   [score, std_gen] = find_generalisability(c_model, gen_iters);
1615                   standard_deviation(n_images * (n - 1) + i) = std_gen;
1616                else   
1617                   score = -1;
1618                end  
1619             elseif (strcmp(model_score_type, 'Mean MSD')),
1620                 if (mod(n_images * (n - 1) + i, cycle) == 0),
1621                     score = measure_model_msd(warped_images);
1622                 else
1623                     score = -1; 
1624                 end
1625             elseif (strcmp(model_score_type, 'All')),
1626               % score = measure_model(c_model.variances, n_modes, model_evaluation_method, c_model);
1627                             % use deafult as score that is primary
1628                 all_scores = [[find_specificity(c_model, images, spec_iters, ref_points_vec)], [find_generalisability(c_model, gen_iters)], [measure_model_msd(warped_images)]];    
1629                            % get all scores above and store them in a
1630                            % matrix
1631             else
1632                error(['Wrong scoring method (', [model_score_type], ') passed. Internal error.']);
1633             end
1634             
1635          end    
1636     
1637     
1638          if(n == n_iterations),
1639             final_score(i, this_set) = score;
1640                 % give the final score for the set if end is reached
1641          end
1642          
1643          score_for_iteration(i,n) = score;
1644                % record the score for this iteration and this image to be later
1645                % plotted for analysis
1646          if (strcmp(model_score_type, 'All')),
1647              all_scores_for_iteration(i,n,1) = all_scores(1);
1648              all_scores_for_iteration(i,n,2) = all_scores(2);
1649              all_scores_for_iteration(i,n,3) = all_scores(3);
1650          end    
1651          if ((strcmp(progress_indicator_type,'Progress Bar')) | (strcmp(progress_indicator_type,'Console Progress Bar'))),
1652            times(this_step) = toc;
1653            time_for_this_iteration = mean (times);
1654          end
1655         end  % end:images
1656         
1657         % The next section deals with point stuffing
1658         
1659         %% Note: cycle scale option here... perform this stage every
1660         %% <add_points_cycle> iterations
1661         if (mod(n, add_points_cycle) == 0)
1662           if (add_points),
1663             for i = 1:n_images,  
1664               for j = 1:image_width - 1, 
1665                 % -Find region in <points> where height difference is high
1666                   difference(j, i) = warped_points(j + 1, i) - warped_points(j, i);
1667               end
1668               
1669               maximum_difference(i) = max(difference(:, i));
1670               for current_position_number = 1:size(difference(:, i), 1),
1671                 if (difference(current_position_number, i) == maximum_difference(i)),
1672                      maximum_difference_location(i) = current_position_number;
1673                         % maximal diffference in position n means that the
1674                         % internal between points(n) and points(n + 1) is
1675                         % maximal
1676                 end
1677               end
1678             end
1679             
1680             
1681               
1682             % -Add points in between that range
1683             % -Set them to some values
1684             % in between.
1685             % -Update <image_width> and all data including
1686             % references
1687             for i = 1:n_images,
1688                 temp_image = warped_images(:, i);
1689                 temp_points = warped_points(:, i);
1690                 for j = maximum_difference_location(i) + 1:image_width,
1691                   warped_points(j + 1,i) = temp_points(j);
1692                   warped_images(j + 1,i) = temp_image(j);
1693               end
1694               
1695                  % shift points one position onwards where stuffing occurs
1696                 warped_points(maximum_difference_location(i) + 1, i) = warped_points(maximum_difference_location(i), i) + maximum_difference(i) / 2;
1697                  % NOTE: THE DIFFERENCE SHOULD NOT BE UNIFORM FOR ALL DATA,
1698                  % It was calculated specifically for each data instance.
1699                  % set new point to lie at the average of its two neighbouring points
1700                 warped_images(maximum_difference_location(i) + 1, i) = mean(warped_images(maximum_difference_location(i) + 1, :)); 
1701                  % -Set new <image> points with similar structures and set their
1702                  % value to one which is projected from the mean.
1703                  % NOTE: REFERENCES ARE BEING CHANGED TOO
1704             end
1705               image_width = image_width + 1;
1706               
1707                 % modify the size of images
1708 %               size(warped_points)
1709 %               size(warped_images)
1710 %               warped_points
1711 %               warped_images
1712             
1713           end % end: point stuffing if
1714         end % end: clcle if
1715         
1716       end  % end:iterations
1717       if (strcmp(progress_indicator_type,'Progress Bar')),
1718                 h = progressbar( h, 0, 'Finalising');
1719       end
1720       if (show_score_vs_iteration),
1721           if (show_all_score_figure_types == 1)            
1722               score_figure = figure('Name', [[handle],[': Score versus iterations for every image warping step - objective function: '],[objective_function]]);
1723                      % create a new figure to show graph of optimisation for
1724                      % all steps and all images individually
1725               hold on;
1726               title('Score versus iterations for every image warping step');
1727               xlabel('Iteration number');
1728               ylabel('Score');
1729               for i=1:n_images,
1730                     plot(score_for_iteration(i,:));
1731               end
1732               hold off;
1733               
1734               figure('Name', [[handle],[': Score versus iterations - objective function: '],[objective_function]]);
1735                      % create a new figure to show graph of optimisation
1736               hold on;
1737               title('Score versus iterations');
1738               xlabel('Iteration number');
1739               ylabel('Mean score');
1740               for i=1:n_iterations,
1741                  mean_score(i) = mean(score_for_iteration(:,i));
1742               end   
1743               plot(mean_score);
1744               hold off;
1745           end   
1746           
1747           if (overlap_scores == 0),
1748              score_figure = figure('Name', [[handle],[': Score versus image warping step - objective function: '],[objective_function]]);
1749                  % create a new figure to show graph of optimisation
1750           else
1751               
1752              score_figure = figure(101);
1753           end
1754           
1755           hold on;
1756           
1757           if (strcmp(model_score_type, 'All')),
1758             title('Proportional score versus image warping step');
1759             xlabel('Warping Step');
1760             ylabel('% Change');  
1761             position_counter = 0;  
1762             for i=1:n_iterations,
1763               for j=1:n_images,
1764                 position_counter = position_counter + 1;
1765                 all_warp_step_score_vector(position_counter,1) = all_scores_for_iteration(j,i,1);
1766                 all_warp_step_score_vector(position_counter,2) = all_scores_for_iteration(j,i,2);
1767                 all_warp_step_score_vector(position_counter,3) = all_scores_for_iteration(j,i,3);
1768               end
1769             end
1770             all_warp_step_percent_vector(1,1) = 0;
1771                                % no change
1772             all_warp_step_percent_vector(2,1) = 0;
1773             all_warp_step_percent_vector(3,1) = 0;
1774                 
1775             for i = 2:n_iterations * n_images,
1776                     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;
1777                     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;
1778                     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;              
1779             end          
1780           else
1781               title('Score versus image warping step');
1782               xlabel('Warping Step');
1783               ylabel('Score');
1784               position_counter = 0;
1785               for i=1:n_iterations,
1786                   for j=1:n_images,
1787                     position_counter = positi