Description of build_1d_bump

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 % ==============================================================
0098 % BUILD_1D_BUMP_MODEL: Builds a 1-D model for some bump data
0099 %                      and returns some related statistics.
0100 %
0101 % Code originally written by Katherine Smith, 2003
0102 %
0103 %    GENERAL
0104 %
0105 %      Omitted. Much of the explanation is included in the section where
0106 %      default arguments are being assigned.
0107 %
0108 %    INPUT/S
0109 %
0110 %      -- No arguments will lead to defaults
0111 %
0112 %      -spline_type:
0113 %           The type of spline to be used. 'single_point' or
0114 %           'multi_point' at present.
0115 %
0116 %      -placement_type:
0117 %           Placement type for 'multi_point':
0118 %                    'grid' or 'edges' or 'random'
0119 %           placement type for 'single_point':
0120 %                    'overlap_grid' or 'edges_and_scale'
0121 %                    or 'random_and_scale'.
0122 %
0123 %      -objective_function:
0124 %           'msd_opt_together' or 'model_opt_separare' or
0125 %           'msd_opt_separate' or 'model_opt_together'.
0126 %
0127 %      .....
0128 %
0129 %    OUTPUT/S
0130 %
0131 %      -status:
0132 %           The status of the operation of this routine (to be returned at
0133 %           the end
0134 %
0135 %      -Figures produced and results displayed are main
0136 %       data to be utilised.
0137 %
0138 %    PENDING WORK
0139 %
0140 %      -
0141 %
0142 %    KNOWN BUG/S
0143 %
0144 %      -
0145 %
0146 %    COMMENT/S
0147 %
0148 %      -Different video outputs are generated but controlling which
0149 %       one is chosen is not trivial. It currently just does the
0150 %       job of experimentation and showing results as pixels
0151 %       or bumps.
0152 %
0153 %    RELATED FUNCTION/S
0154 %
0155 %      See other functions in same directory, filenames
0156 %      beginning with "build" in particular. There seem to be much
0157 %      overlap due to cut-and-paste.
0158 %
0159 %    ABOUT
0160 %
0161 %      -Created:     December 29rd, 2003
0162 %      -Last update: March 30th, 2004
0163 %      -Revision:    0.9.2
0164 %      -Author:      R. S. Schestowitz, University of Manchester
0165 % ==============================================================
0166       
0167 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0168 %% variables (parameters) initialisation           %%
0169   
0170 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0171 %% Input Defaults              %%
0172 
0173 % When the function is not called from the GUI, various arguments must be
0174 % initialised so that the function performs a registration task.
0175 
0176 % A possible extension would allow for _partial_ arguments to be given so
0177 % that a script can run different experiments and then logs are necessary
0178 % too.
0179 
0180 if (nargin == 0),
0181     disp(' ');
0182     disp('Argument(s) omitted. Using defaults instead.');
0183     disp(' ');
0184          % notify user about default values being used.
0185     spline_type = 'single_point';
0186          % SMITH: spline type 'single_point' or 'multi_point'
0187          % SMITH: placement type for multipoint: 'grid' or 'edges' or 'random'
0188          % SMITH: placement type for singlepoint: 'overlap_grid' or 'edges_and_scal or 'random_and_scale'
0189     placement_type = 'random_and_scale';
0190          % controls knot-point selection on the grid
0191     objective_function = 'model_opt_together';
0192          % choice of objective function
0193     apply_intensity_offset = 1;
0194          % 0: Offset should not be applied to data
0195          % 1: Offset enabled
0196     exit_at_end = 0;
0197          % indicates whether AART quits at the end of registration
0198     description = '';
0199          % a name provided for describing registration
0200     cycle = 0;
0201           % cycle size for user display of objective function evaluation
0202     produce_record = 0;
0203           % indicates if results should be displayed to the user in the
0204           % console
0205     save_in_log = 0;
0206           % indicates if an HTML file or network of files will be used to
0207           % record the results
0208     close_at_end = 0;
0209           % indicates if the GUI should be closed when registration is done
0210     save_movie = 0;
0211           % a flag for movie output
0212     force_reference = 0;
0213          % should movie be saved (boolean)
0214          % 0: No
0215          % 1: Yes
0216     frames_per_second = 15;
0217          % number of frames per second for the video generated.
0218     n_standard_deviations_display = 2;
0219          % the number of standard deviation to be displayed when a model is
0220          % presented to the user
0221     n_modes_display = 2;
0222          % as above, for number of modes
0223     perturbation_method = 'CPS warp';
0224          % the type of perturbation applied to the data
0225     offset_extent = 0;
0226          % the extent of perturbation or the level of distortion which
0227          % accompanies the parameter 'perturbation_method'
0228     precision_automation_type = 'default';
0229        % the way to set precision for the optimiser
0230     do_plot = 0;
0231          % should bumps be plotted as it changes? (Boolean)
0232          % (Old comment: this also enables frames being captured which slows down the
0233          % process. Capturing of frames and plotting can potentially be separated to make
0234          % this routine more useful).
0235     movie_type = 1;
0236          % type of movie to save
0237          % 0: bumps as pixels. do_plot may be 0 for a pixel clip to be recorded.
0238          % 1: bump plots. Remember to set do_plot to 1.
0239     show_warps_while_optimising = 0; 
0240          % indicates if the warps applied to the data should be dynamically
0241          % displayed to the user
0242     movie_name = 'undefined';
0243           % the name of the AVI movie to be saved. It will be folloed by
0244           % '.avi' eventually and contain the movie.
0245     format = 'png';
0246          % if data is saved as an image, this paramater defines the format
0247          % as which it shall be saved.
0248     show_warp_targets = 0;
0249          % defines whether or not the target of the warps, i.e. that at
0250          % alignment, should be displayed to the user.
0251     dynamic_precision = 0;
0252        % defines whether or not precision is chosen wisely at each stage
0253     gen_iters = 25;
0254        % generalisability iterations = 25;
0255        % specificity iterations
0256     spec_iters = 25;
0257        % specificity iterations
0258     progress_bar = 1;
0259          % display progress bar or text instead
0260          % 0: Hierarchical Text
0261          % 1: Progress bar
0262          % 2: Quiet
0263          % 3: Console progress bar
0264     model_score_type = 'Default';
0265          % Defined the way in which a model is scored.
0266     force_first_reference = 0;
0267          % depending on its numeric value, this parameter defines the
0268          % reference selection method
0269     user_menu = 0;
0270          % indicates the way in which registration parameters are entered.
0271          % If not via pop-up menu, then probably via top menu in the
0272          % application.
0273     number_of_bins = 50;
0274          % 0: parameters need to be set in file
0275          % 1: menu input
0276     background_red = 1;
0277     background_green = 1;
0278     background_blue = 1;
0279          % the background colourn for the plots
0280     distance_type = 'sum_of_squared_distances';
0281          % The way to meausre distance from mean
0282     show_registration_target = 0;
0283          % indicates if the estimated target of registration should be
0284          % super-imposed in score plots.
0285     produce_statistics = 0;
0286          % should statistics be displayed at the end (boolean)
0287     see_model_before_each_iteration = 0;
0288          % should a combined model be produced and displayed before each
0289          % iteration, this should be set to 1
0290     show_pixels = 0;
0291          % show pixels representation (boolean). Mandatory for pixel movie
0292          % capture
0293     show_score_vs_iteration = 0;
0294          % when enabled, evaluation of the objective function will be
0295          % displayed against the iteration at which it was evaluated.
0296     find_perfect_warp = 1;
0297          % This argument should not be tempered with. It defined the need
0298          % to check for the value and state at conceived registration
0299          % target
0300     show_model_after_optimisation = 0;
0301          % if enabled, then a model produces from the data is displayed
0302          % after registration.
0303     show_specificity_statistics = 0;
0304          % more booleans which are self-explanatory
0305     verbose_score = 0;
0306          % boolean: display score or not after optimisation
0307     show_score_bar = 0;
0308          % The score bar gives a visual indication of the 'goodness' of a
0309          % model construced from the data
0310     average_smooth_enabled = 0;
0311     gaussian_smooth_enabled = 0;
0312          % flags indicating the smoothness to be applied, if any at all
0313     n_points = 5;
0314          % an argument that is passed to the model optimisation function.
0315          % it refers to the number of knotpoint, e.g. for random and scale.
0316     window = 1;
0317          % the size of the window for average smoothing or sigma for
0318          % gaussian smoothing
0319     show_all_score_figure_types = 0;
0320          % defines the need for an abundant number of plots depicting
0321          % quality of model versus iteration
0322     plot_before_and_after = 0;
0323          % will only ever be used if the two above are set
0324     bump_width = 0.2;
0325     bump_width_variation = 0.5;
0326     bump_height = 0.2;
0327     bump_height_variation = 0.3;
0328     bump_position_freedom = 0.9;
0329     initial_diminish_factor = 2;
0330     smoothness_factor = 17;
0331          % set the bump parameters
0332     precision_required = 1e-7;
0333         % the precision required from the optimiser
0334     model_evaluation_method = 1;
0335                       % 1: Determinant
0336                       % 0: MDL
0337                       % 2: Sum of log of Eigenvalues
0338     show_score_using_model = 0;
0339     handle = 'Default';
0340       % give figures a fefault name
0341     n_modes = 20;
0342       % maximal number of modes to be accounted for when evaluating a model
0343     variation_kept = 0.99;
0344       % keep 99%
0345       % number of modes in the model to be accounted for
0346     retain_original_peak = 0;
0347       % determines whether or not to keep the same maximal data height
0348     error_verbosity = 'normal';
0349       % the level of verbosity for errors
0350     initialise_warps_close_to_target = 0;
0351       % determines the requirement to start registration with warps near
0352       % their target (at completion of registration)
0353     shape_weight = 0.06;
0354       % the constant weight Ws for the model constructed. 1 is 100%, 0 is
0355       % 0% and legal range is 0..1
0356     weighting_normalisation_method = 'edge';
0357      % the scheme for Ws weighting selection
0358     run_all_objective_functions = 0;
0359        % do not try running all the objective functions
0360     called_from_gui = 0;
0361        % indicate that GUI was not used in this case
0362     kernel = 'Exponential';
0363        % the kernel type for PDF objective function
0364     save_data = 0;
0365        % Indicated if the data needs to be saved to files after opening or
0366        % generation
0367     peak_type = 'average_position';
0368          % The type of peak to be retained
0369     load_data = 0;
0370        % indicates if data needs to be loaded from file rather than being
0371        % generated.
0372     save_bumps = 0;
0373        % indicated if bumps need to be saved as images (not data)
0374     draw_warp_curve = 0;
0375        % An option for overlapping warps display
0376     draw_curve_by_images = 0;
0377        % The option for warp target maps (images x iterations)
0378     overlap_scores = 0;
0379       % defines the need for scores to be drawn on the same space if
0380       % different objective functions are compared by being run
0381       % sequentially
0382     show_mesh = 0;
0383      % simply indicates if a mesh of the data needs to be displayed
0384     show_target = 0;
0385         % do not show perfectly aligned images estimate
0386     verbose = 'off';
0387          % be verbose or not during general-purpose optimisation
0388     save_stats = 0;
0389        % indicated if statistics need to be saved as HTML
0390     display_shape_model = 0;
0391       % If enabled, shape model is shown
0392     display_intensity_model = 0;
0393       % likewise for intensity model
0394     n_images = 2;
0395              % how many images to try to register in each set
0396     n_iterations = 2;
0397              % how many iterations of the optimisation to run, default was 50.
0398     n_sets = 1;
0399              % how many sets of images to register (in order to get better statistics)
0400     image_width = 50;
0401              % the number of pixels or the lenght of the data vector
0402     figure_handle_number = 1;
0403              % increments applied when figures are saved
0404     plot_alignment = 0; 
0405              % Indicates if a plot for alignment and perturbation should be
0406              % displayed
0407     rotate_surface = 0;
0408              % Indicated if the surface should be rotated horizontally at
0409              % each frame
0410     draw_final_warps = 0;
0411              % indicates if final warps need to be displayed
0412     add_points = 0;
0413              % point stuffing where difference in warp curve is greatest
0414     add_points_cycle = 1;
0415              % the frequency (measured in terms of iteration) for stuffing;
0416 else
0417     called_from_gui = 1;
0418         % indicate that default values were not set and that the function
0419         % was called from REGISTER
0420 end
0421 
0422        
0423 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0424 %% Input / Output              %%
0425 
0426 rand('state',sum(100*clock))
0427           % seeding for random numbers.
0428 
0429 figure_handle_number = 1;
0430           % A number to be incremented every time a figure is saved. It
0431           % manages the number of figures that are generated.
0432 warning off;
0433          % keeps st_pca.m quiet
0434 error_found = 0;
0435          % status flag
0436 current_viewpoint = 45;
0437          % the starting horizontal viewpoint for surfaces
0438          % Optimisation note: initiliase this variables only when surfaces
0439          % are displayed.
0440 status = ['                      '];
0441          % set intitial status message. Messages are aggregated before
0442          % being returned to the GUI and displayed to the user.
0443 if (run_all_objective_functions == 1),
0444          % If all objective functions will run sequentially, then some
0445          % output parameters need to be updated.
0446     colour = 'r';
0447         % set initial colour for overlaid plots
0448     style = ':';
0449         % similarly to style of line in the overlaid plots
0450     progress_bar = 2;
0451 end
0452      % the progress indicators will not work for this routine so disable
0453      % them
0454 
0455 if strcmp(movie_name, 'undefined'),
0456     movie_name = ['movie_', spline_type, '_', placement_type, '_', objective_function,'_width', num2str(image_width), '_images', num2str(n_images), '_iterations', num2str(n_iterations), '.avi'];
0457          % name for the AVI to be saved
0458          % The name will attempt to encompass all construction details
0459          % It will replace the default name given when the function is not
0460          % called from the GUI. If called from the GUI, it usually have a
0461          % user defined name or 'Default.avi'
0462 end
0463 
0464 
0465 frame_counter = 0;
0466          % counter for video frames. Do not fiddle with this variable. It
0467          % keeps track of the current frame number to be captured.
0468          
0469 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0470 %% Others (most are obsolete and are remains of old code.)  %%
0471                        
0472 % white_ctr = 0;
0473          % white counter? Was used for simple bumps in code that is now
0474          % commented out
0475 % ref_shift = 0.2;
0476        % shift in reference image??
0477 % max_shift = 0.2;
0478        % maximum allowed shift?
0479 % step = 0.1;
0480        % used in two lines below only
0481 % los = zeros(n_images,floor((max_shift-ref_shift)/step));
0482 % his = zeros(n_images,floor((max_shift-ref_shift)/step));
0483        % not understood yet
0484 % blurred = 0;
0485        % image blurring has been applied to images (boolean)
0486 % min_error = 0;
0487 % max_error = 1;
0488        % define error range??
0489 % error_step = 0.1;
0490 % n_error_steps = (max_error-min_error) / error_step;
0491 
0492 
0493          
0494 switch model_evaluation_method,
0495     case 0
0496         model_evaluation_method = 'minimum_description_length';
0497     case 1
0498         model_evaluation_method = 'determinant';     
0499     case 2
0500         model_evaluation_method = 'sum_of_log_eigenvalues';
0501     case 3
0502         model_evaluation_method = 'specificity';
0503     case 4
0504         model_evaluation_method = 'generalisability';    
0505     otherwise
0506         error('Unknown model evaluation method. Wrong numeric identifier passed.');
0507 end
0508          % The conversion above can one day be embedded in the GUI to make
0509          % everything more readable. Values are overwritten here too.
0510          
0511 switch progress_bar,
0512     case 0
0513      progress_indicator_type = 'Hierarchical Text';
0514     case 1
0515      progress_indicator_type = 'Progress Bar';
0516     case 2
0517      progress_indicator_type = 'Quiet'; 
0518     case 3
0519      progress_indicator_type = 'Console Progress Bar';
0520     otherwise
0521      error('Callback value for progress indicator is not recognised');
0522 end
0523 
0524          % do type conversions for input to make code below readable
0525          
0526 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0527 %% Begin Set initialisation    %%
0528 
0529 switch user_menu,
0530   case 0
0531     
0532     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0533     %% Optimisation Controllers: Default   %%
0534     
0535     % Original settings from 2003: The values below will be defaults if user menu does not appear
0536     
0537     n_images = 2;
0538              % how many images to try to register in each set
0539     n_iterations = 2;
0540              % how many iterations of the optimisation to run, default was 50.
0541     n_sets = 1;
0542              % how many sets of images to register (in order to get better statistics)
0543     image_width = 50;
0544              % the number of pixels or the length of the data vector
0545              
0546     % A few extra popular choices, added March 2004 to speed up experiments
0547     % construction and customisation
0548     
0549   case 1
0550     [n_sets, n_iterations, n_images, image_width] = get_parameters_from_user('Default');
0551        % Prompts the user for the main registartion parameters using menus
0552 
0553   case 2    
0554     n_images = 10;
0555     n_iterations = 50;
0556     n_sets = 1;
0557     image_width = 50;
0558     
0559   case 3   
0560     n_images = 5;
0561     n_iterations = 5;
0562     n_sets = 1;
0563     image_width = 50; 
0564     
0565   case 4    
0566     n_images = 10;
0567     n_iterations = 10;
0568     n_sets = 1;
0569     image_width = 50;
0570     
0571   case 5   
0572     n_images = 20;
0573     n_iterations = 10;
0574     n_sets = 1;
0575     image_width = 50; 
0576     
0577   case 6
0578     n_images = 20;
0579     n_iterations = 1000;
0580     n_sets = 1;
0581     image_width = 50;
0582     
0583   case 7 
0584     n_images = 10;
0585     n_iterations = 1000;
0586     n_sets = 1;
0587     image_width = 50;
0588     
0589   case 8
0590     n_images = 10;
0591     n_iterations = 10000;
0592     n_sets = 1;
0593     image_width = 50;
0594     
0595   otherwise
0596     error('Unknown parameters option. Internal failure.');
0597     
0598 end 
0599 
0600 % in all the above cases, meaningful names rather than numbers are a
0601 % possible extension which requires time.
0602 
0603 
0604 
0605 
0606 tic;
0607        % start counting time
0608 for this_set = 1:n_sets
0609        % for all sets
0610        % Old: user gets status about set number. Multiple sets generate reliable
0611        %      statistics
0612   tic;
0613        % This is called again (stack of time counters????? Doesn't seem so)
0614        % This will need to be changed (one tic omitted) for the statistics
0615        % to work properly. Only one stopwatch can work at a time.
0616        % Smith: for each set of n_sets images, build a model
0617        % calculate obj fn values from these models, plot the mean with error
0618        % bars given by 1 stddev away
0619   if (strcmp(progress_indicator_type,'Hierarchical Text')),      
0620     disp(['calculating set ', num2str(this_set), ' out of ', num2str(n_sets), ' sets...']);
0621       % display command-line status indicators
0622   elseif (strcmp(progress_indicator_type,'Progress Bar')),
0623       gui_active(1);
0624           % will add an abort button
0625       if (strcmp(getenv('OS'), 'Linux')),
0626          h = progressbar([], 0, 'Preparing for registration', [['Registration Progress on '], [get_host]]);
0627       else
0628          h = progressbar([], 0, 'Preparing for registration', 'Registration Progress');
0629       end
0630         % initiate progress bar if necessary
0631       h = progressbar( h, 0, 'Fetching Data');
0632   end    
0633 
0634   % initiate score bar if necessary
0635 
0636 
0637   if (show_score_bar == 1),
0638      score_bar = waitbar(0, 'Score');
0639   end
0640      
0641 
0642   switch load_data,
0643       
0644    case 0
0645       [imagelist, images, points] = make_1d_bump(n_images, image_width...
0646           , bump_width, bump_width_variation, bump_height, bump_height_variation...
0647           , bump_position_freedom, initial_diminish_factor, smoothness_factor...
0648           , save_bumps, format);
0649        % create a set of images given several parameters
0650        
0651       points = -1 + 2 * (points - 1) / (image_width - 1);
0652         % Smith: normalise points from -1 to 1
0653              
0654       if (called_from_gui == 1 & save_bumps == 1),
0655         status = [[status];['(**) Images saved     '];['                      ']];
0656       end
0657         % notify user using the status bar that data has been saved.
0658  
0659        %% keep = 0.999999; NOTE (RSS Feb. 2004): became obsolete and
0660        %% replaced by <variation kept> which relates to PCA
0661        
0662        % Nov. 2003: Is this precision required? Perhaps in optimisation.
0663        
0664    case 1
0665       output = open([[handle],['.mat']]);
0666       points = output.points;
0667       images = output.images;
0668                   % copy the data structures from the file
0669                   
0670    case 2
0671       output = open('Default.mat');
0672       points = output.points;
0673       images = output.images; 
0674       
0675    case 3
0676       output = open('flat.mat');
0677       points = output.points;
0678       images = output.images; 
0679       
0680    case 4
0681       output = open('plain.mat');
0682       points = output.points;
0683       images = output.images;  
0684       
0685    case 5
0686       output = open('plainran.mat');
0687       points = output.points;
0688       images = output.images;  
0689       
0690    case 6
0691       [imagelist, images, points] = make_elliptic_bumps(n_images, image_width, save_bumps, format, bump_height, bump_position_freedom, bump_height_variation, bump_width_variation);
0692        % create a set of images given several parameters
0693       points = -1 + 2 * (points - 1) / (image_width - 1);
0694              
0695       if (called_from_gui == 1 & save_bumps == 1),
0696         status = [[status];['(**) Images saved     '];['                      ']];
0697       end
0698       
0699    case 7
0700       [imagelist, images, points, dummy1, dummy2] = make_1d_images(n_images, image_width, 0.2);
0701        % 0.2. is the value used by Kate for <white width>
0702       if (save_bumps == 1),
0703           msgbox('Data save is not supported with simple bumps.');
0704       end
0705        % create a set of images given several parameters
0706       points = -1 + 2 * (points - 1) / (image_width - 1);
0707              
0708       if (called_from_gui == 1 & save_bumps == 1),
0709         status = [[status];['(**) Images not saved '];['                      ']];
0710       end    
0711       
0712    otherwise
0713       error('Invalid value passed from open submenu. Internal error.');
0714       
0715   end     
0716   
0717   
0718   if (average_smooth_enabled),
0719       images = average_smooth(images, window);
0720   end  
0721   
0722   if (gaussian_smooth_enabled),
0723       images = gaussian_smooth(images, window); 
0724   end
0725        % smooth all images - used to be disabled and set <blurred> to 1
0726        
0727 
0728   ref_points_vec = points(:,1);
0729   ref_image_vec = images(:,1);
0730        % set first image generated to be reference when calculating a
0731        % perfect warp and also as reference in pair-wise objective
0732        % functions. This has been made more generalised later on.
0733   % ref_hi = his(1);
0734   % ref_lo = los(1);
0735        % and get its upper and lower boundaries. his and lois were
0736        % generated and returned as output for simple bumps. They were
0737        % used to find the perfect warps and defined points of importance in
0738        % the input data.
0739   
0740   % Now choose the image which will not be changed. Note a possible
0741 % optimisation: only if a reference is forced, then the code below should
0742 % be run to choose the reference. The section below that is to do with
0743 % alignment means however that for target to be shown an immutable image
0744 % must always be chosen
0745 
0746       switch force_first_reference,
0747        case 0        
0748            immutable_image = 2;
0749                 % see comments later about the choice of first or second
0750                 % image as an immutable reference. These options will only
0751                 % have an effect if a reference is forced (not as in the
0752                 % default reference in MSD-based algorithms)
0753                 
0754                 % OPTIMISATION NOTE: score can also be set to 0 in all
0755                 % cases without breaking registration. It would then need
0756                 % not be repeated in all the cases below.
0757                 
0758        case 1
0759            immutable_image = 1;
0760            score = 0;
0761                 % the score needs to be initilialised because no previous
0762                 % scoring will be avaliable
0763                 
0764        case 2
0765           % closest data to average (mean)
0766            immutable_image = find_image_closest_to_mean(images, image_width, n_images, 'closest', distance_type);
0767            score = 0;
0768            
0769        case 3
0770               % random reference
0771            immutable_image = ceil(rand * n_images);
0772            score = 0;
0773            
0774        case 4
0775           % closest data to average (mean)
0776            immutable_image = find_image_closest_to_mean(images, image_width, n_images, 'farthest', distance_type);
0777            score = 0;
0778            
0779        case 5
0780            % highest peak
0781            immutable_image = find_image_by_height(images, image_width, n_images, 'highest');
0782            score = 0;
0783            
0784        case 6
0785            % lowest peak
0786            immutable_image = find_image_by_height(images, image_width, n_images, 'lowest');
0787            score = 0;
0788            
0789        otherwise
0790           error('Reference error. Wrong values passed from callback.');
0791           
0792       end  
0793 
0794       % NOTE PROBLEM: If teh refrence chosen is the first then evaluation at
0795       % the start is 0 and the optimisation plot is not suitable
0796        
0797        
0798        
0799        
0800   %%%%% Do warp analysis here
0801   
0802   % The check for perfect warps does not work for corrective offsets,
0803   % non-bump data, and data that is significantly smoothed. In all of these
0804   % cases, the flag below must be set to 0, or else piece wise linear will
0805   % not work.
0806   check_perfect_warps = find_perfect_warp;
0807           % currently always enabled, but computationally expensive and in
0808           % some cases not necessary. A test for the value of this part to
0809           % be implemented in the future.
0810       if (check_perfect_warps == 1),   
0811           if (strcmp(progress_indicator_type,'Progress Bar')),
0812              h = progressbar( h, 0, 'Calculating Correct Warps');
0813           end
0814         % (Old - November 2003 - Smith)
0815           % correctly_warped_points = points;
0816           % correctly_warped_images = images;
0817           for i = 1:n_images,
0818               % Smith: now bung in the right deformations to get them back
0819               % NOTE: THE VARIABLE ref_image_vec BELOW ASSUMED A FIRST
0820               % IMAGE REFERENCE although this is not the case for the <force
0821               % reference> option which makes it the second (implicitly).
0822               % Somthing richer than a second image might be used in
0823               % future.
0824               
0825               %possible optimisation: do not interpolate for the first
0826               %image which is usually a reference and is not affected
0827               %anyhow
0828               %images(:,i)
0829             correctly_warped_points(i,:) = find_target_warp(images(:,immutable_image), images(:,i), image_width);
0830               % correctly_warped_points(i,:)
0831               % get the warp that appears most desirable
0832             correctly_warped_points(i,:) = -1 + 2 * (correctly_warped_points(i,:) - 1) / (image_width - 1);
0833               % normalise points so that they lie in the range x where -1 < x < 1
0834             correctly_warped_points_vec = correctly_warped_points(i,:);              
0835             image_vec = images(:,i);
0836               % RSS: copy original image and copy piece-wise linear warp
0837               % Smith: resample at the warped points - point on reg. grid has value of point at same index in warpy grid
0838               % need to interpolate
0839               % size(ref_points_vec)
0840               % size(image_vec)
0841               % size(correctly_warped_points_vec')
0842             correctly_warped_images(i,:) = interp1(ref_points_vec, image_vec, correctly_warped_points_vec);
0843             % correct_model = build_model(correctly_warped_images,correctly_warped_points,1,'','edge',0);
0844             % converging_score(i) =
0845             % measure_model(correct_model.variances,50);
0846             % correctly_warped_images(i,:)
0847           end
0848           % figure;
0849           % plot(log(converging_score));
0850           % title('Log of score as model converges');
0851       end
0852       % figure;
0853       % (test curve) plot(correctly_warped_points(2,:))
0854       
0855   
0856 
0857       if (show_warp_targets == 1),
0858          warp_handle = figure('Name','On the left: piece-wise linear optimal registration warps; On the right: results of warp');
0859          % title('On the left: piece-wise linear optimal registration warps; On the right: results of warp');
0860          save_handle = figure(warp_handle);
0861          hold on;
0862          subplot(n_images, 2, 1);
0863          for i=1:n_images,
0864              subplot(n_images, 2, i * 2 - 1);
0865              hold on;
0866              axis off;
0867              plot(correctly_warped_points(i,:));
0868              axis([-1 image_width -1 1]);
0869              hold off;
0870          end
0871          for i=1:n_images,
0872              subplot(n_images, 2, i * 2);
0873              hold on;
0874              axis off;
0875              plot(correctly_warped_images(i,:));
0876              axis([0 image_width 0 1]);
0877              hold off;
0878          end
0879          hold off;
0880          if (save_in_log == 1),
0881            saveas(save_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
0882            figure_handle_number = figure_handle_number + 1;
0883          end
0884       end
0885       
0886   
0887     
0888   if (initialise_warps_close_to_target == 1),
0889       % offset_extent = 5;
0890       % perturbation_method = 'CPS warp';
0891       
0892     if (plot_alignment == 1),
0893         % if a plot is required to show this step of perturbation after
0894         % alignment the right-hand side solumn will be occupied by the
0895         % plots
0896       subplot_alignment_fig = figure('Name',[[handle],[': Left: Aligned images ; Right (if existent): Aligned images after perturbation']]);
0897       alignment_handle=figure(subplot_alignment_fig);
0898       for i=1:n_images,
0899        subplot(n_images, 2, 2 * i - 1);
0900        hold on;
0901        plot(correctly_warped_images(i,:));
0902        % title(['Unwarped image #',num2str(i)]);
0903        axis([0 image_width 0 1]);
0904        grid on;
0905        hold off;
0906       end
0907     end
0908       
0909       if (offset_extent == 0),
0910           % if no perturbation is required then the values can be simply
0911           % copied.
0912          for i = 1:n_images,
0913            % copy registration target
0914           points(:,i) = correctly_warped_points(i,:)';
0915            % images(:,i) = correctly_warped_images(i,:)';
0916          end
0917       else
0918           % if some noise is required, then add it in the form of Gaussian
0919           % noise or CPS warps to points and interpolate original images
0920           % accordingly. To be controlled from GUI.
0921            for i = 1:n_images,
0922                 tweaked_points(i,:) = apply_perturbation_to_points(correctly_warped_points(i,:), offset_extent, perturbation_method);
0923         %         figure
0924         %         hold on
0925         %         plot(tweaked_points(i,:));
0926         %         plot(correctly_warped_points(i,:));
0927         %         hold off
0928                    % apply noise to the correct points
0929                 % images(:,i) = interp1(ref_points_vec, images(:,i), tweaked_points(i,:))';
0930                    % interpolation to get current images
0931                 points = tweaked_points';
0932                    % copy the points found to become current
0933                    % inconsistent column and row form is due to initial
0934                    % implementation. Optimally (when our task becomes optimising performance),
0935                    % all transposes should be ridden of.
0936            end
0937            if (plot_alignment == 1),
0938               for i = 1:n_images,
0939                 perturbed_aligned_images(:,i) = interp1(ref_points_vec, images(:,i), points(:,i));
0940               end    
0941                 % inerpolate to show some plots at this earlier stage
0942               alignment_handle = figure(subplot_alignment_fig);
0943               for i = 1:n_images,
0944                    subplot(n_images, 2, 2 * i);
0945                    hold on;
0946                    plot(perturbed_aligned_images(:,i));
0947                    % title(['Unwarped image #',num2str(i)]);
0948                    axis([0 image_width 0 1]);
0949                    grid on;
0950                    hold off;
0951                end
0952              if (save_in_log == 1)
0953                saveas(alignment_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
0954                figure_handle_number = figure_handle_number + 1;
0955              end
0956            end
0957       end
0958       
0959   end
0960   
0961   
0962   if (save_data == 1),
0963       % note that data will be saved with modified warp curves if option is
0964       % set. This is intentional. Data in this cotext is vectors, not
0965       % images.
0966       save([[handle],['.mat']], 'images','points');
0967       if (called_from_gui == 1),
0968         status = [[status];['(**) Data saved       '];['                      ']];
0969       end
0970   elseif (save_data == 2),
0971       % if data needs to be saved as default
0972       save([[handle],['.mat']], 'images','points');
0973       save(['Default.mat'], 'images','points');
0974       if (called_from_gui == 1),
0975         status = [[status];['(**) Data saved (Def) '];['                      ']];
0976       end
0977   end
0978 
0979 
0980   
0981   if (retain_original_peak),  
0982        for current_image = 1:n_images,
0983            maximum_height(current_image) = max(images(:, current_image));
0984               % find out the maximum height of each data instance
0985            
0986            for j=1:image_width,
0987                switch peak_type,
0988                    case 'initial_position'
0989                          if (images(j,current_image) == maximum_height(current_image)),
0990                            maximum_position(current_image) = j;
0991                             % find the position where height is maximal. Currently
0992                             % this search is very inefficient.
0993                          end
0994                          
0995                    case 'average_position'
0996                          if (images(j,current_image) == maximum_height(current_image)),
0997                            maximum_position(current_image) = j;
0998                          end
0999                          % same as above
1000                end
1001            end
1002        end
1003  
1004   
1005        switch peak_type,
1006          case 'average_position'
1007            average_position = floor(mean(maximum_position));
1008            for current_image = 1:n_images,
1009              maximum_position(current_image) = average_position;
1010            end 
1011        end
1012    end
1013            
1014            
1015            
1016            
1017 % plot3(1:2,1:2,images(1,:))
1018 %
1019 %
1020 %       s = surface(L, ...
1021 %         'EdgeColor','none', ...
1022 %         'FaceColor',[0.8 0.2 0.2], ...
1023 %         'FaceLighting','phong', ...
1024 %         'AmbientStrength',0.3, ...
1025 %         'DiffuseStrength',0.6, ...
1026 %         'Clipping','off',...
1027 %         'BackFaceLighting','lit', ...
1028 %         'SpecularStrength',1.1, ...
1029 %         'SpecularColorReflectance',1, ...
1030 %         'SpecularExponent',7);
1031 %     l1 = light('Position',[40 100 20], ...
1032 %         'Style','local', ...
1033 %         'Color',[0 0.7 0.7]);
1034 %     l2 = light('Position',[.5 -1 .4], ...
1035 %         'Color',[1 1 0]);
1036 
1037 % modify the progress bar if required
1038 
1039   if (strcmp(progress_indicator_type,'Progress Bar')),
1040            h = progressbar( h, 0, 'Pre-processing Data');
1041   end
1042   
1043 % depending on whether intensity offset was applied, record the number of
1044 % columns required for the plots
1045        
1046   if (apply_intensity_offset),
1047       vertical_pane_location = 3;
1048                        % Need is for only 3 columns in plot
1049   else
1050       vertical_pane_location = 2;
1051                        % Need is for only 2 columns in plot
1052   end
1053   % Note: intensity offsets are considered a possible way of solving the
1054   % problem, but not natural form for our input data. Hence, saving of data
1055   % precedes this offset.
1056 
1057 
1058 % prepare a figure if a plot of all warp curves needs to be rendered
1059        
1060   if (draw_warp_curve == 1),
1061        overlap_curve_fig = figure('Name','Overlapping Warp Curves');
1062   end
1063   
1064   % if pixel representation is required, prepare the figure
1065   
1066   if (show_pixels == 1),
1067     pixels_figure = figure('Name','Pixel Representation');
1068   end
1069   
1070   % and likewise for meshes/surfaces
1071 
1072   if (show_mesh == 1),
1073     mesh_subfigure = figure('Name','1-D Data Mesh');
1074   end 
1075   
1076   % if a matrix of warps is required, preapre the figure and draw the
1077   % initial warps in its left side
1078   
1079   if (draw_curve_by_images == 1),
1080     curve_subplot_fig = figure('Name',[['Warp Curves for all data - objective function: '],[objective_function]]);
1081     curves = figure(curve_subplot_fig);
1082     % title('Warp Curves for all data');
1083     for i=1:n_images,
1084        subplot(n_images, n_iterations + 1, i * (n_iterations + 1) - n_iterations);
1085        hold on;
1086        axis off;
1087        plot(points(:,i));
1088        plot(correctly_warped_points(i,:),'r-')
1089        if (i == 1 ),
1090          legend('Initial warp', 'Target warp'); 
1091             % only display the legend once at the top
1092        end
1093        axis([-1 image_width -1 1]);
1094        hold off;
1095     end
1096   end
1097   
1098 
1099 %%             END OF INITIALISATION               %%
1100 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1101 
1102 
1103 
1104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1105 %% Warping images, then modelling warped images    %%
1106 %% attempt to register images by warping.          %%
1107 
1108 
1109 
1110   if (plot_before_and_after),
1111     for i=1:n_images,
1112        images_before_offset(:,i) = images(:,i);
1113        % copy these original values in case 3 columns are necessary for plot representation of the data.
1114     end
1115   end
1116   
1117   % draw the original data at the leftmost column
1118   
1119   if (do_plot),
1120     subplot_fig = figure('Name',[[handle],[': Left: Original images ; Right: Images subjected to warping']]);
1121     H=figure(subplot_fig);
1122     for i=1:n_images,
1123        subplot(n_images, vertical_pane_location, (vertical_pane_location * i) - vertical_pane_location + 1);
1124        hold on;
1125        plot(images(:,i));
1126        whitebg(H, [background_red, background_green, background_blue]); 
1127        % title(['Unwarped image #',num2str(i)]);
1128        axis([0 image_width 0 1]);
1129        grid on;
1130        hold off;
1131     end
1132   else
1133     subplot_fig = 0;
1134       % pass a void figure to the optimisation functions
1135   end
1136   
1137   % draw the data after intensity offset has been applied (when applicable)
1138   
1139   if (apply_intensity_offset),
1140       for index = 1:n_images,
1141          peak (index) = max(images(:,index));
1142       end
1143       peak_reference = max (peak);  
1144       for index = 1:n_images,
1145         peak_difference = peak_reference - max(images(:,index));
1146         images(:,index) = images(:,index) + peak_difference;
1147         if (do_plot),
1148             subplot(n_images, vertical_pane_location, (vertical_pane_location * index) - vertical_pane_location + 2);
1149             hold on;
1150             plot(images(:,index));
1151             % title(['Unwarped image #',num2str(i)]);
1152             axis([0 image_width 0 1]);
1153             grid on;
1154             hold off;
1155         end
1156       end        
1157   end 
1158 %% NOTE: Change here in traversal above (RSS, December 2003)
1159   
1160     still_to_go = 1;
1161           % condition for while loop. This will enable the program to
1162           % continually try all objective functions.
1163     
1164     
1165     while (still_to_go == 1),
1166     
1167     if (run_all_objective_functions == 0)
1168         still_to_go = 0;
1169            % break out of while loop after just one iteration
1170     else
1171            % If new objective functions are added, an entry must be added
1172            % here to include it in the more comprehensive cycle.
1173            
1174            % try all objective functions sequentially
1175        switch objective_function,    
1176            case 'all'
1177            % meaning that it is in its first iteration
1178              objective_function = 'model_opt_together';
1179            case 'model_opt_together'
1180              objective_function = 'model_opt_separate';    
1181            case 'model_opt_separate'
1182              objective_function = 'msd_opt_together';
1183            case 'msd_opt_together'
1184              objective_function = 'msd_opt_separate';
1185            case 'msd_opt_separate'
1186              objective_function = 'mi_opt_together';
1187            case 'mi_opt_together'
1188              objective_function = 'mi_opt_separate';
1189            case 'mi_opt_separate'
1190              objective_function = 'nmi_opt_together';
1191            case 'nmi_opt_together'
1192              objective_function = 'nmi_opt_separate';
1193            case 'nmi_opt_separate'          
1194              objective_function = 'mixed_msd_model';
1195            case 'mixed_msd_model'
1196              objective_function = 'model_post_msd';
1197            case 'model_post_msd'        
1198              objective_function = 'tfc_eccv';
1199            case 'tfc_eccv'
1200              objective_function = 'wavelet_complexity';
1201            case 'wavelet_complexity'
1202              objective_function = 'model_pairs_together';
1203            case 'model_pairs_together'
1204              objective_function = 'model_pairs_separate';
1205            case 'model_pairs_separate'
1206              objective_function = 'void'; 
1207              still_to_go = 0;
1208               % signal to break out of the while loop because the test of
1209               % all objective functions is finished.
1210            otherwise
1211              error('Problem with <Run All Objective Functions> option');
1212        end
1213     end   
1214       % warped_images = images
1215       % size(images)
1216       for i=1:n_images,
1217         warped_images(:,i) = interp1(ref_points_vec, images(:,i), points(:,i)')';
1218       end
1219       % size(warped_images)
1220       warped_points = points;
1221           % set this initial state to plot the examples before change
1222   
1223           
1224        
1225       
1226       
1227   
1228       if ((strcmp(progress_indicator_type,'Console Progress Bar')) | (strcmp(progress_indicator_type,'Progress Bar')))
1229          progress_bar_position = 0;
1230          time_for_this_iteration = 10;
1231          total_images = n_images * n_iterations;
1232          step_size = 50 / total_images;
1233             % these variables are used for some progress statistics
1234       end   
1235         % initialise progress bar
1236         
1237 
1238 
1239         
1240       for n = 1:n_iterations,
1241         
1242         if (dynamic_precision == 1),
1243            iterations_ratio = n_iterations / n;
1244              % get iterations ratio. In future, looking at the optimisation
1245              % state and behaviour, for example whether it is stuck,
1246              % should give better results.
1247            precision_required = set_precision(iterations_ratio, precision_automation_type);
1248         end
1249         
1250         if (show_score_bar == 1),
1251             % subplot('position',[0,0,1,0.1]);
1252             Epsilon = 0.0001;
1253             c_model = build_model(warped_images, warped_points, variation_kept, '', weighting_normalisation_method, shape_weight);
1254             % waitbar((measure_model(c_model.variances, n_modes, 2, c_model) - log(Epsilon.^n_modes)) / log(Epsilon.^n_modes), score_bar);
1255             waitbar((measure_model(c_model.variances, n_modes, 'sum_of_log_eigenvalues', c_model) ...
1256                / log(Epsilon.^n_modes) - 0.8) * 5, score_bar);
1257             % (measure_model(c_model.variances, n_modes, 2, c_model) - log(Epsilon.^n_modes)) / log(Epsilon.^n_modes) * 100
1258             % axis off;
1259             % subplot(1,1,1);
1260         end            
1261           
1262           
1263           
1264         if (show_mesh == 1),
1265           meshview = figure(mesh_subfigure);
1266           whitebg(meshview, [background_red, background_green, background_blue]);
1267                % set viewpoint
1268           [X,Y] = meshgrid(1:n_images * 5 + 2, -image_width / 2 - 1:1:image_width / 2);
1269                % two zero buffers and 5 rows of each bump
1270           for i = 1:n_images,
1271             Z(2:image_width + 1,1 + i * 5) = warped_images(1:image_width,i);
1272             Z(2:image_width + 1,1 + i * 5 - 1) = warped_images(1:image_width,i);
1273             Z(2:image_width + 1,1 + i * 5 - 2) = warped_images(1:image_width,i);
1274             Z(2:image_width + 1,1 + i * 5 - 3) = warped_images(1:image_width,i);
1275             Z(2:image_width + 1,1 + i * 5 - 4) = warped_images(1:image_width,i);
1276             % copy the data values
1277           end
1278           Z(:,1) = zeros(1:image_width, 1);
1279                 % set bottom buffer
1280           Z(:, n_images * 5 + 2) = zeros(1:image_width, 1);
1281                 % and upper buffer to close mesh
1282           Z(1,:) = zeros(1, 1:n_images * 5);
1283                 % set left buffer
1284           Z(image_width+2,:) = zeros(1, 1:n_images * 5);
1285                 % and right buffer to close mesh
1286           switch mesh_type,  
1287            case 'Surface Red'  
1288             surf(X,Y,Z,'FaceColor', 'red', 'EdgeColor', 'none');
1289             camlight left;
1290             lighting phong;
1291             axis([0 n_images*5 -image_width/2 image_width/2 0 1]);     
1292            case 'Surface Blue'     
1293             surf(X,Y,Z,'FaceColor', 'blue', 'EdgeColor', 'none');
1294             camlight left;
1295             lighting phong;
1296             axis([0 n_images*5 -image_width/2 image_width/2 0 1]);     
1297            case 'Polygons'   
1298             mesh(X,Y,Z);
1299             axis([0 n_images*5 -image_width/2 image_width/2 0 1]);     
1300            case 'Illuminated Surface'
1301             if (rotate_surface),
1302               current_viewpoint = current_viewpoint + 1;
1303                % shift viewpoint slightly
1304             end
1305             view(current_viewpoint, 37.5); 
1306             surface(X,Y,Z, ...
1307               'EdgeColor','none', ...
1308               'FaceColor',[0.5 0.5 0.5], ...
1309               'FaceLighting','phong', ...
1310               'AmbientStrength',0.3, ...
1311               'DiffuseStrength',0.6, ... 
1312               'Clipping','off',...
1313               'BackFaceLighting','lit', ...
1314               'SpecularStrength',1.1, ...
1315               'SpecularColorReflectance',1, ...
1316               'SpecularExponent',7);
1317             if (n == 1), % set lights only at the start
1318                 l1 = light('Position',[40 100 20], ...
1319                   'Style','local', ...
1320                   'Color',[0 0.5 0.7]);
1321                 l2 = light('Position',[.5 -1 .4], ...
1322                   'Color',[1 1 0]);
1323             end
1324            otherwise
1325             error ('Mesh type not recognised. Internal error.');
1326           end
1327           
1328           axis off; 
1329           if (movie_type == 2 & save_movie == 1),
1330                   frame_counter = frame_counter + 1;
1331                   M(frame_counter) = getframe(meshview);
1332           end
1333           if (save_in_log == 1 & n == n_iterations), 
1334              saveas(meshview,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1335              figure_handle_number = figure_handle_number + 1;
1336           end
1337          
1338         end  
1339           
1340         if (strcmp(progress_indicator_type,'Hierarchical Text'))  
1341            % iterations aim to get good statistical results by averaging
1342           disp(['    +-Running iteration number ',num2str(n), ' out of ',num2str(n_iterations), ' iterations...']);
1343         end
1344         
1345            % Smith: correct registration
1346            % attempt to register
1347                
1348         if (see_model_before_each_iteration),  
1349           w_c_model = build_model(warped_images, warped_points, variation_kept,'Optimised warp', weighting_normalisation_method, shape_weight);
1350                   % create an appearance model...
1351           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)]]);
1352                   % and show the model before each iteration while saving
1353                   % the figure handle
1354           if (movie_type == 3 & save_movie == 1),
1355                 % if moview needs to be saved and it is a model
1356                   frame_counter = frame_counter + 1;
1357                   M(frame_counter) = getframe(model_view);
1358           end        
1359         end          
1360         for i = 1:n_images,
1361            % start iteration time counter
1362          if (strcmp(progress_indicator_type,'Hierarchical Text')),
1363            disp(['        +-Warping image ',num2str(i),' out of ',num2str(n_images), ' images in total...']);
1364          end
1365          
1366          this_step = ((n-1) * n_images) + i;
1367                     
1368          if (strcmp(progress_indicator_type,'Progress Bar')),
1369            tic;
1370            progress_bar_position = progress_bar_position + step_size;
1371            steps_remaining = (n_images * n_iterations) - this_step;
1372            percentage_progress = ceil(this_step / (n_images * n_iterations) * 100);
1373            minutes = floor(time_for_this_iteration * steps_remaining / 60);
1374            seconds = rem(floor(time_for_this_iteration *  steps_remaining), 60);
1375            max_count   = 50;
1376                % The console-based GUI is 50 units wide
1377            if (seconds > 9),
1378                h = progressbar( h, step_size/50,[num2str(percentage_progress),'% Completed ---- Estimated remaining time: ', num2str(minutes), ':', num2str(seconds)]);
1379            else
1380                h = progressbar( h, step_size/50,[num2str(percentage_progress),'% Completed ---- Estimated remaining time: ', num2str(minutes), ':0', num2str(seconds)]);
1381            end    
1382            if ~gui_active,
1383                     error_found = 1;
1384                     break; 
1385                % If aborted was pressed
1386            end
1387               % use the third-party GUI-like progress bar (added February
1388               % 2004)
1389          end
1390          if (strcmp(progress_indicator_type,'Console Progress Bar')),
1391                                % if console progress bar
1392            tic;
1393            progress_bar_position = progress_bar_position + step_size;
1394            clc;
1395            disp(['|=================================================|']);
1396            progress_string='|';       
1397            for counter = 1:floor(progress_bar_position),
1398                progress_string = [progress_string, '#'];
1399            end
1400            disp(progress_string);
1401            disp(['|================= ',num2str(floor(progress_bar_position * 2)),'% completed =================|']);
1402                           % display progress per cent
1403            steps_remaining = (n_images * n_iterations) - this_step;
1404            minutes = floor(time_for_this_iteration * steps_remaining / 60);
1405            seconds = rem(floor(time_for_this_iteration *  steps_remaining), 60);
1406            disp(' ');
1407            if (seconds > 9),
1408              disp(['            Estimated remaining time: ', num2str(minutes), ':', num2str(seconds)]);
1409                           % show time indicators
1410            else
1411              disp(['            Estimated remaining time: ', num2str(minutes), ':0', num2str(seconds)]);
1412            end
1413          end
1414          
1415          image_vec = warped_images(:,i);
1416          points_vec = warped_points(:,i);     
1417          
1418          if (show_pixels == 1),   
1419              appended_vector=[];
1420              for current_image=1:n_images,
1421                appended_vector=[[appended_vector],[warped_images(:,current_image),[warped_images(:,current_image)]],[warped_images(:,current_image)]];
1422                    % append bump data 3 times for width which is visible
1423              end
1424            
1425                  %% appended_vector=appended_vector * 256
1426                  %% use if scale of 0..255 used
1427     
1428              G = figure(pixels_figure);
1429              hold on;
1430               % title('Bumps as pixels');
1431              imshow(appended_vector',[0,1]);
1432                        % show as an image
1433              if (movie_type == 0 & save_movie == 1),
1434                                    % Note: already nested in a show_pixels
1435                                    % 'if' statement
1436                   frame_counter = frame_counter + 1;
1437                   M(frame_counter) = getframe(G);
1438              end
1439              hold off;
1440          end
1441          
1442          if (movie_type == 1 & save_movie == 1 & do_plot == 1),
1443                    frame_counter = frame_counter + 1;
1444                    M(frame_counter) = getframe(H);
1445          end
1446                    % capture pixels frame
1447          
1448          
1449          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
1450                        % Somthing richer than a second image might be used
1451                        % in future e.g. mean image
1452              switch objective_function,     
1453                  case 'model_opt_together'
1454                      
1455 %                      if (i==1 & n==1),
1456 %                          score=0;
1457 %                      end
1458 %
1459 %                      last_score = score;
1460 %
1461 %                      % disp(score);
1462 %                   if (i==3)
1463 %                       figure(101);
1464 %                       plot(image_vec);
1465 %                   end
1466                      [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);
1467                      
1468                      %if (n==1),
1469                      %  [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);
1470                      %else
1471                      %   [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);
1472                      % end
1473                      
1474                      
1475 %                      if (last_score < score),
1476 %                          disp(num2str((n - 1) * i + i))
1477 %                      end
1478 
1479                  case 'model_pairs_together'    
1480                      [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);    
1481                  case 'model_pairs_separate'  
1482                      [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);
1483                  case 'model_opt_separate'
1484                      [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);
1485                  case 'msd_opt_together'         
1486                      [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);
1487                  case 'msd_opt_separate'
1488                      [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);
1489                  case 'mi_opt_together'       
1490                      [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);
1491                  case 'mi_opt_separate'
1492                      [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);
1493                  case 'nmi_opt_together'        
1494                      [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);
1495                  case 'nmi_opt_separate'
1496                      [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);
1497                  case 'mixed_msd_model'
1498                         if (((n / n_iterations < 2 / 10) & (n / n_iterations >= 1 / 10)) |...
1499                             ((n / n_iterations < 4 / 10) & (n / n_iterations >= 3 / 10)) |...
1500                             ((n / n_iterations < 6 / 10) & (n / n_iterations >= 5 / 10)) |...
1501                             ((n / n_iterations < 8 / 10) & (n / n_iterations >= 7 / 10)) |...
1502                             ((n / n_iterations < 10 / 10) & (n / n_iterations >= 9 / 10))),
1503                           [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);
1504                         else
1505                           [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);
1506                         end                    
1507                  case 'model_post_msd'
1508                         if (n / n_iterations < 1 / 2),
1509                             % for the first half
1510                           [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);
1511                         else
1512                             % and for the second half
1513                           [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);
1514                         end
1515                  case 'tfc_eccv'
1516                      % msgbox('TFC ECCV 2004: Under Construction.');
1517                      pdf_evalution_method = 'Default';
1518                      [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);         
1519                  case 'wavelet_complexity'
1520                      % msgbox('Wavelets: Under Construction.');
1521                      wavelet_evalution_method = 'Default';
1522                      [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);                      
1523                  case 'void'
1524                     % This is used for repeated experiments with overlaid
1525                     % plots. It is not well tested and score may need to be
1526                     % assigned to a value here
1527                      break;
1528                      % It is not obsolutely clear id the above is necessary
1529                      % or works
1530                  otherwise
1531                      error(['Unknown objective function: ', objective_function]);
1532              end
1533          else
1534             % score = 0;
1535             warped_image = warped_images(:,i);
1536             warped_point = warped_points(:,i); 
1537                   % copy some values to make up for the skip
1538             if (do_plot == 1),
1539               figure(subplot_fig);
1540               subplot(n_images, vertical_pane_location , vertical_pane_location * i, 'replace');
1541               plot(warped_image);
1542               axis([0 image_width 0 1]);
1543               drawnow;
1544             end
1545                   % draw the plot
1546          end
1547 
1548 
1549          warped_images(:,i) = warped_image;
1550            
1551 
1552          if (retain_original_peak), 
1553            warped_images(maximum_position(i), i) = maximum_height(i);
1554               % always retain the maximum height
1555          end     
1556               
1557               
1558               
1559               
1560          warped_points(:,i) = warped_point;
1561          
1562          if (draw_warp_curve == 1),
1563            figure(overlap_curve_fig);
1564            hold on;
1565            plot(warped_points(:,i));
1566            hold off;
1567            if (save_in_log == 1 & n == n_iterations & i == n_images), 
1568                saveas(overlap_curve_fig,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1569                figure_handle_number = figure_handle_number + 1;
1570            end
1571          end     
1572          
1573          if (draw_curve_by_images == 1),
1574              figure(curves);
1575              subplot(n_images, n_iterations + 1, i * (n_iterations + 1) - (n_iterations - n));
1576              hold on;
1577              axis off;
1578              plot(warped_points(:,i));
1579              axis([-1 image_width -1 1]);
1580              hold off;
1581              if (save_in_log == 1 & n == n_iterations & i == n_images), 
1582                saveas(curves,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1583                figure_handle_number = figure_handle_number + 1;
1584              end
1585          end
1586          if (show_pixels == 1 & save_in_log == 1 & n == n_iterations & i == n_images), 
1587                saveas(G,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1588                figure_handle_number = figure_handle_number + 1;
1589          end
1590          if (do_plot == 1 & save_in_log == 1 & n == n_iterations & i == n_images), 
1591                saveas(H,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1592                figure_handle_number = figure_handle_number + 1;
1593          end
1594          
1595          % figure(subplot_fig),title(['Warped images, after image ',num2str(i),' iteration ',num2str(n)]);
1596          % dummy = waitforbuttonpress;
1597          if (show_score_using_model == 1),           
1598              % if evaluation is through model then always overwrite score
1599             c_model = build_model(warped_images, warped_points, variation_kept, '', weighting_normalisation_method, shape_weight);
1600                 % first build a model
1601             if (strcmp(model_score_type, 'Default')),
1602                 score = measure_model(c_model.variances, n_modes, model_evaluation_method, c_model);
1603             elseif (strcmp(model_score_type, 'Specificity')),
1604                 if (mod(n_images * (n - 1) + i, cycle) == 0),
1605                     % disp('ok')
1606                   [score, std_spec] = find_specificity(c_model, images, spec_iters, ref_points_vec);
1607                   standard_deviation(n_images * (n - 1) + i) = std_spec;
1608                 else
1609                     % set to -1 (or copy from predecessor - the original
1610                     % intension
1611                   score = -1; 
1612                 end
1613             elseif (strcmp(model_score_type, 'Generalisability')), 
1614                if (mod(n_images * (n - 1) + i, cycle) == 0),
1615                   [score, std_gen] = find_generalisability(c_model, gen_iters);
1616                   standard_deviation(n_images * (n - 1) + i) = std_gen;
1617                else   
1618                   score = -1;
1619                end  
1620             elseif (strcmp(model_score_type, 'Mean MSD')),
1621                 if (mod(n_images * (n - 1) + i, cycle) == 0),
1622                     score = measure_model_msd(warped_images);
1623                 else
1624                     score = -1; 
1625                 end
1626             elseif (strcmp(model_score_type, 'All')),
1627               % score = measure_model(c_model.variances, n_modes, model_evaluation_method, c_model);
1628                             % use deafult as score that is primary
1629                 all_scores = [[find_specificity(c_model, images, spec_iters, ref_points_vec)], [find_generalisability(c_model, gen_iters)], [measure_model_msd(warped_images)]];    
1630                            % get all scores above and store them in a
1631                            % matrix
1632             else
1633                error(['Wrong scoring method (', [model_score_type], ') passed. Internal error.']);
1634             end
1635             
1636          end    
1637     
1638     
1639          if(n == n_iterations),
1640             final_score(i, this_set) = score;
1641                 % give the final score for the set if end is reached
1642          end
1643          
1644          score_for_iteration(i,n) = score;
1645                % record the score for this iteration and this image to be later
1646                % plotted for analysis
1647          if (strcmp(model_score_type, 'All')),
1648              all_scores_for_iteration(i,n,1) = all_scores(1);
1649              all_scores_for_iteration(i,n,2) = all_scores(2);
1650              all_scores_for_iteration(i,n,3) = all_scores(3);
1651          end    
1652          if ((strcmp(progress_indicator_type,'Progress Bar')) | (strcmp(progress_indicator_type,'Console Progress Bar'))),
1653            times(this_step) = toc;
1654            time_for_this_iteration = mean (times);
1655          end
1656         end  % end:images
1657         
1658         % The next section deals with point stuffing
1659         
1660         %% Note: cycle scale option here... perform this stage every
1661         %% <add_points_cycle> iterations
1662         if (mod(n, add_points_cycle) == 0)
1663           if (add_points),
1664             for i = 1:n_images,  
1665               for j = 1:image_width - 1, 
1666                 % -Find region in <points> where height difference is high
1667                   difference(j, i) = warped_points(j + 1, i) - warped_points(j, i);
1668               end
1669               
1670               maximum_difference(i) = max(difference(:, i));
1671               for current_position_number = 1:size(difference(:, i), 1),
1672                 if (difference(current_position_number, i) == maximum_difference(i)),
1673                      maximum_difference_location(i) = current_position_number;
1674                         % maximal diffference in position n means that the
1675                         % internal between points(n) and points(n + 1) is
1676                         % maximal
1677                 end
1678               end
1679             end
1680             
1681             
1682               
1683             % -Add points in between that range
1684             % -Set them to some values
1685             % in between.
1686             % -Update <image_width> and all data including
1687             % references
1688             for i = 1:n_images,
1689                 temp_image = warped_images(:, i);
1690                 temp_points = warped_points(:, i);
1691                 for j = maximum_difference_location(i) + 1:image_width,
1692                   warped_points(j + 1,i) = temp_points(j);
1693                   warped_images(j + 1,i) = temp_image(j);
1694               end
1695               
1696                  % shift points one position onwards where stuffing occurs
1697                 warped_points(maximum_difference_location(i) + 1, i) = warped_points(maximum_difference_location(i), i) + maximum_difference(i) / 2;
1698                  % NOTE: THE DIFFERENCE SHOULD NOT BE UNIFORM FOR ALL DATA,
1699                  % It was calculated specifically for each data instance.
1700                  % set new point to lie at the average of its two neighbouring points
1701                 warped_images(maximum_difference_location(i) + 1, i) = mean(warped_images(maximum_difference_location(i) + 1, :)); 
1702                  % -Set new <image> points with similar structures and set their
1703                  % value to one which is projected from the mean.
1704                  % NOTE: REFERENCES ARE BEING CHANGED TOO
1705             end
1706               image_width = image_width + 1;
1707               
1708                 % modify the size of images
1709 %               size(warped_points)
1710 %               size(warped_images)
1711 %               warped_points
1712 %               warped_images
1713             
1714           end % end: point stuffing if
1715         end % end: clcle if
1716         
1717       end  % end:iterations
1718       if (strcmp(progress_indicator_type,'Progress Bar')),
1719                 h = progressbar( h, 0, 'Finalising');
1720       end
1721       if (show_score_vs_iteration),
1722           if (show_all_score_figure_types == 1)            
1723               score_figure = figure('Name', [[handle],[': Score versus iterations for every image warping step - objective function: '],[objective_function]]);
1724                      % create a new figure to show graph of optimisation for
1725                      % all steps and all images individually
1726               hold on;
1727               title('Score versus iterations for every image warping step');
1728               xlabel('Iteration number');
1729               ylabel('Score');
1730               for i=1:n_images,
1731                     plot(score_for_iteration(i,:));
1732               end
1733               hold off;
1734               
1735               figure('Name', [[handle],[': Score versus iterations - objective function: '],[objective_function]]);
1736                      % create a new figure to show graph of optimisation
1737               hold on;
1738               title('Score versus iterations');
1739               xlabel('Iteration number');
1740               ylabel('Mean score');
1741               for i=1:n_iterations,
1742                  mean_score(i) = mean(score_for_iteration(:,i));
1743               end   
1744               plot(mean_score);
1745               hold off;
1746           end   
1747           
1748           if (overlap_scores == 0),
1749              score_figure = figure('Name', [[handle],[': Score versus image warping step - objective function: '],[objective_function]]);
1750                  % create a new figure to show graph of optimisation
1751           else
1752               
1753              score_figure = figure(101);
1754           end
1755           
1756           hold on;
1757           
1758           if (strcmp(model_score_type, 'All')),
1759             title('Proportional score versus image warping step');
1760             xlabel('Warping Step');
1761             ylabel('% Change');  
1762             position_counter = 0;  
1763             for i=1:n_iterations,
1764               for j=1:n_images,
1765                 position_counter = position_counter + 1;
1766                 all_warp_step_score_vector(position_counter,1) = all_scores_for_iteration(j,i,1);
1767                 all_warp_step_score_vector(position_counter,2) = all_scores_for_iteration(j,i,2);
1768                 all_warp_step_score_vector(position_counter,3) = all_scores_for_iteration(j,i,3);
1769               end
1770             end
1771             all_warp_step_percent_vector(1,1) = 0;
1772                                % no change
1773             all_warp_step_percent_vector(2,1) = 0;
1774             all_warp_step_percent_vector(3,1) = 0;
1775                 
1776             for i = 2:n_iterations * n_images,
1777                     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;
1778                     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;
1779                     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;              
1780             end          
1781           else
1782               title('Score versus image warping step');
1783               xlabel('Warping Step');
1784               ylabel('Score');
1785               position_counter = 0;
1786               for i=1:n_iterations,
1787                   for j=1:n_images,
1788                     position_counter = position_counter + 1;
1789                     warp_step_score_vector(position_counter) = score_for_iteration(j,i);
1790                   end
1791               end
1792           end
1793           if (run_all_objective_functions == 1),
1794               if strcmp(colour,'r'),
1795                 colour = 'g';
1796               elseif strcmp(colour,'g')
1797                 colour = 'b';
1798               else
1799                 colour = 'r';
1800               end 
1801               
1802               if strcmp(style,':'),
1803                 style = '-';
1804               elseif strcmp(style,'-')
1805                 style = '--';
1806               elseif strcmp(style,'--')
1807                 style = '-.'; 
1808               else
1809                 style = ':';
1810               end 
1811               plot(warp_step_score_vector,[[colour],[style]]);
1812               % The next line needs to be updated when new objective
1813               % functions are added
1814               legend('Joint MSD','Sequential MSD','Joint Model-based','Sequential Model-based','Joint MI',...
1815                      'Sequential MI', 'Joint NMI', 'Sequential NMI', 'PDF', 'Wavelets', 'Joint Model-based Pair-wise', 'Sequential Model-based Pair-wise', 'Void');
1816           else
1817                 % if a normal single plot will suffice
1818                 if (strcmp(model_score_type, 'All')),
1819                          plot (all_warp_step_percent_vector(1,:), 'r');
1820                          plot (all_warp_step_percent_vector(2,:), 'g');
1821                          plot (all_warp_step_percent_vector(3,:), 'b');
1822                          legend('Specificity', 'Generalisability', 'Mean MSD');
1823                 else 
1824                           if (strcmp(progress_indicator_type,'Progress Bar')),
1825                               h = progressbar( h, 0, 'Generating Scores Plot');
1826                           end    
1827                          % find any -1 entries and "destroy" them to get
1828                          % correct cycle evaluation
1829                          % optimisation tip: use cycle to pick out the
1830                          % right entries rather than the naive -1 flag
1831                          counter = 1;
1832                          for index=1:n_images*n_iterations,
1833                              if (warp_step_score_vector(index) ~= -1),
1834                                  if (strcmp(model_score_type, 'Generalisability') | strcmp(model_score_type, 'Specificity')),
1835                                    standard_deviation_truncated(counter) = standard_deviation (index);
1836                                  end   
1837                                  warp_step_score_vector_truncated(counter) =  warp_step_score_vector(index);
1838                                  counter = counter + 1;
1839                              end
1840                          end    
1841                          if (strcmp(model_score_type, 'Generalisability') | strcmp(model_score_type, 'Specificity')),
1842                              errorbar(warp_step_score_vector_truncated, standard_deviation_truncated);
1843                              title([model_score_type, ' every ', num2str(cycle), ' steps']);
1844                          elseif (strcmp(model_score_type, 'Mean MSD'))
1845                               % for mean MSD
1846                              plot (warp_step_score_vector_truncated); 
1847                              title(['Mean MSD every ', num2str(cycle), ' steps']);
1848                          else % model defaults
1849                              plot (warp_step_score_vector); 
1850                              title(['Registration score versus iterations']);
1851                          end
1852                          
1853                 end
1854                 if (show_target == 1),
1855                     Epsilon = 0.0001;
1856                     plot (repmat([log(Epsilon.^n_modes)], 1,n_iterations * n_images),'--');
1857                     legend('Current Evaluation','Log eigenvalues of perfectly registered images');
1858                 end
1859                 if (show_registration_target == 1),
1860                     w_c_model = build_model(correctly_warped_images', correctly_warped_points', variation_kept, 'Optimised warp', weighting_normalisation_method, shape_weight);      
1861                     model_score =  measure_model(w_c_model.variances, n_modes, model_evaluation_method, w_c_model);
1862                     plot(repmat(model_score, 1,n_iterations * n_images),'r-.');
1863                     if (show_target == 1),
1864                       legend('Current Evaluation','Log eigenvalues of perfectly registered images','Log eigenvalues of registration target');
1865                     else 
1866                       legend('Current Evaluation','Log eigenvalues of registration target');
1867                     end  
1868                 end    
1869           end   
1870           
1871           hold off;
1872           if (save_in_log == 1),
1873                score_handle = figure(score_figure);
1874                saveas(score_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1875                figure_handle_number = figure_handle_number + 1;
1876           end
1877       end
1878       
1879       w_c_model = build_model(warped_images, warped_points, variation_kept, 'Optimised warp', weighting_normalisation_method, shape_weight);      
1880          % Smith: build model from warped images
1881       w_intensity_total_vars = w_c_model.intensity_model.total_var;
1882       w_shape_total_vars = w_c_model.shape_model.total_var;
1883          % get these two values to be used later
1884       model_score(this_set) = measure_model(w_c_model.variances, n_modes, model_evaluation_method, w_c_model);
1885       msd_score(this_set) = measure_model_msd(warped_images);
1886       shape_modes(this_set) = w_c_model.n_shape_modes;
1887       intensity_modes(this_set) = size(w_c_model.intensity_model.pcs,2);
1888       shape_variance(this_set) = sum(w_c_model.shape_model.variances);
1889       intensity_variance(this_set) = sum(w_c_model.intensity_model.variances);
1890          % copy some values to be displayed for statistics
1891       
1892       if (show_model_after_optimisation == 1),
1893           model_handle = show_combined_model(w_c_model, ref_points_vec, n_standard_deviations_display, n_modes_display, [[handle],[': Final combined model built from optimised warps - objective function: '],[objective_function]], image_width);
1894           if (save_in_log == 1),
1895                saveas(model_handle,[[handle],'-', num2str(figure_handle_number)], 'jpg');
1896                figure_handle_number = figure_handle_number + 1;
1897           end
1898       end
1899       
1900       
1901       if (display_intensity_model == 1),
1902         [intensity_model.pcs,intensity_model.variances,intensity_model.params, intensity_model.mean, intensity_model.var_r, intensity_model.total_var, intensity_model.impdata] = st_pca(warped_images', variation_kept);
1903           % perform PCA
1904         intensity_model.sd = std(intensity_model.params);
1905               % get standard deviation from the models (of shape and intensity)
1906         intensity_total_vars = intensity_model.total_var;
1907         total_intensity_variance = sum(intensity_model.variances);
1908         fig_title = 'After registration: ';
1909         intensity_model_handle = show_intensity_model(intensity_model, 2, 2, fig_title, image_width);
1910         if (save_in_log == 1),
1911                saveas(intensity_model_handle, [[handle], '-', num2str(figure_handle_number)], 'jpg');
1912                figure_handle_number = figure_handle_number + 1;
1913         end
1914       end  
1915       
1916       if (display_shape_model == 1),  
1917         [shape_model.pcs,shape_model.variances,shape_model.params, shape_model.mean, shape_model.var_r, shape_model.total_var, shape_model.impdata] = st_pca(warped_points', variation_kept);
1918         shape_model.sd = std(shape_model.params);
1919         shape_total_vars = shape_model.total_var;
1920         total_shape_variance = sum(shape_model.variances);
1921         shape_model_handle = show_shape_model(shape_model, ref_points_vec, ref_image_vec, 2, image_width);
1922         if (save_in_log == 1),
1923                saveas(shape_model_handle, [[handle], '-', num2str(figure_handle_number)], 'jpg');
1924                figure_handle_number = figure_handle_number + 1;
1925         end        
1926       end  
1927           
1928       if (show_specificity_statistics),
1929         [w_c_model.mean_specificity, w_c_model.std_specificity] = find_specificity(w_c_model, images, spec_iters, ref_points_vec);
1930         disp(' ');
1931         disp([['Objective function: '], [objective_function]]);
1932         disp(['Mean of specificity:               ', num2str(w_c_model.mean_specificity)]);
1933         disp(['Standard deviation of specificity: ', num2str(w_c_model.std_specificity)]);
1934       end
1935       
1936       if (draw_final_warps),
1937           final_warps = figure('Name',[['Warp Curves - objective function: '], [objective_function]]);
1938           final_warps_handle = figure(final_warps);
1939           % title('Warp Curves at end);
1940           for i=1:n_images,
1941             subplot(n_images, 1, i);
1942             hold on;
1943             axis off;
1944             plot(warped_points(:,i));
1945             plot(correctly_warped_points(i,:),'r-')
1946             if (i == 1 ),
1947               legend('Final warp', 'Piece-wise linear target warp'); 
1948               % only display the legend once at the top
1949             end
1950             axis([-1 image_width -1 1]);
1951           end
1952           hold off;       
1953           if (save_in_log == 1),
1954                saveas(final_warps_handle, [[handle], '-', num2str(figure_handle_number)], 'jpg');
1955                figure_handle_number = figure_handle_number + 1;
1956           end  
1957        end
1958      
1959     
1960     % fig_title = 'Automatically aligned: ';
1961     % show_shape_model(w_shape_model, ref_points_vec, ref_image_vec, 2, white_width)
1962     % show_intensity_model(w_intensity_model, 2, white_width, fig_title)
1963     % show_combined_model(w_c_model, ref_points_vec, 2, white_width, fig_title)
1964     
1965     
1966     
1967     
1968     %%             END OF MODEL-BUILDING                    %%
1969     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1970     
1971     
1972     
1973     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1974     %% Now provide some statistics on the                   %%
1975     %% simulation/experiment                                %%
1976     
1977       
1978       %% Warping images, then modelling warped images    %%
1979     %% attempt to register images by warping.          %%
1980     
1981       if (plot_before_and_after == 1),
1982         subplot_fig = figure('Name', [[handle],[': Left: Before registration; Right: After registration - Objective function: '], [objective_function]]);
1983         H=figure(subplot_fig);
1984         for i=1:n_images,
1985            subplot(n_images, vertical_pane_location, (vertical_pane_location * i) - vertical_pane_location + 1);
1986            hold on;
1987            plot(images_before_offset(:,i));
1988            % title(['Unwarped image #',num2str(i)]);
1989            axis([0 image_width 0 1]);
1990            grid on;
1991            hold off;
1992         end
1993         if (apply_intensity_offset),
1994           for index = 1:n_images,
1995              peak (index) = max(images(:,index));
1996           end
1997           peak_reference = max (peak);  
1998           for index = 1:n_images,
1999             peak_difference = peak_reference - max(images(:,index));
2000             images(:,index) = images(:,index) + peak_difference;
2001             subplot(n_images, vertical_pane_location, (vertical_pane_location * index) - vertical_pane_location + 2);
2002             hold on;
2003             plot(images(:,index));
2004             % title(['Unwarped image #',num2str(i)]);
2005             axis([0 image_width 0 1]);
2006             grid on;
2007             hold off;
2008           end        
2009         end    
2010     
2011           % pass a void figure to the optimisation functions
2012         for i=1:n_images,
2013            subplot(n_images, vertical_pane_location, i * vertical_pane_location);
2014            hold on;
2015            plot(warped_images(:,i));
2016            % title(['Unwarped image #',num2str(i)]);
2017            axis([0 image_width 0 1]);
2018            grid on;
2019            hold off;
2020         end  
2021         if (save_in_log == 1), 
2022                saveas(subplot_fig,[[handle],'-', num2str(figure_handle_number)], 'jpg');
2023                figure_handle_number = figure_handle_number + 1;
2024         end
2025       else
2026         subplot_fig = 0; 
2027       end
2028     
2029       
2030  end % while
2031  
2032  t = toc;
2033       
2034  if (strcmp(progress_indicator_type,'Hierarchical Text')),
2035          disp(['Time for set ',num2str(this_set),': ',num2str(t)]);
2036  end
2037  time(this_set) = t;
2038  
2039 end % end set
2040 
2041 
2042 
2043 if (produce_statistics == 1),
2044     disp(' ');
2045     disp(['======================================================']); 
2046     disp(['================= SET STATISTICS =====================']);
2047     disp(['| Objective function:                       ', objective_function]);
2048     disp(['| Mean of match score:                      ', num2str(mean(final_score(:)))]);    
2049     disp(['| Standard deviation of match score:        ', num2str(std(final_score(:)))]);
2050     disp(['| Mean of model score:                      ', num2str(mean(model_score(:)))]);
2051     disp(['| Log of mean of model score:               ', num2str(log(mean(model_score(:))))]);
2052     disp(['| Standard deviation of model score:        ', num2str(std(model_score(:)))]);
2053     disp(['| Mean of MSD score:                        ', num2str(mean(msd_score(:)))]);
2054     disp(['| Standard deviation of MSD score:          ', num2str(std(msd_score(:)))]);
2055     disp(['| Mean of shape modes:                      ', num2str(mean(shape_modes(:)))]);
2056     disp(['| Standard deviation of shape modes:        ', num2str(std(shape_modes(:)))]);
2057     disp(['| Mean of intensity modes:                  ', num2str(mean(intensity_modes(:)))]);
2058     disp(['| Standard deviation of intensity modes:    ', num2str(std(intensity_modes(:)))]);
2059     disp(['| Mean of shape variance:                   ', num2str(mean(shape_variance(:)))]);
2060     disp(['| Standard deviation of shape variance:     ', num2str(std(shape_variance(:)))]);        
2061     disp(['| Mean of intensity variance:               ', num2str(mean(intensity_variance(:)))]);    
2062     disp(['| Standard deviation of intensity variance: ', num2str(std(intensity_variance(:)))]);    
2063     disp(['| Mean time:                                ', num2str(mean(time(:)))]);    
2064     disp(['| Standard deviation of time:               ', num2str(std(time(:)))]);    
2065     t = toc;
2066          % get time from the counter, initiated by <tic>
2067     disp(['| Total time:                               ', num2str(t)]);
2068          % display total time
2069     disp(['======================================================']); 
2070 end
2071 
2072 if (save_stats == 1),
2073     fid = fopen([[handle], ['.htm']],'w');
2074     add_html_headers(fid);
2075     fprintf(fid, ['\n<CENTER><TABLE BORDER=BOX CELLPADDING=10 CELLSPACING=3 BACKGROUND="../../scbg.gif"><TR>']);
2076     fprintf(fid, ['\n<TD><H2>Ststistics for ', [handle], '</H2>']);
2077     fprintf(fid, ['\n<BR><U>Objective function:</U> ', objective_function]);
2078     fprintf(fid, ['\n<BR><U>Mean of match score:</U> ', num2str(mean(final_score(:)))]);
2079     fprintf(fid, ['\n<BR><U>Standard deviation of match score:</U> ', num2str(std(final_score(:)))]);
2080     fprintf(fid, ['\n<BR><U>Mean of model score:</U> ', num2str(mean(model_score(:)))]);
2081     fprintf(fid, ['\n<BR><U>Log of mean of model score:</U> ', num2str(log(mean(model_score(:))))]);
2082     fprintf(fid, ['\n<BR><U>Standard deviation of model score:</U> ', num2str(std(model_score(:)))]);
2083     fprintf(fid, ['\n<BR><U>Mean of MSD score:</U> ', num2str(mean(msd_score(:)))]);
2084     fprintf(fid, ['\n<BR><U>Standard deviation of MSD score:</U> ', num2str(std(msd_score(:)))]);
2085     fprintf(fid, ['\n<BR><U>Mean of shape modes:</U> ', num2str(mean(shape_modes(:)))]);
2086     fprintf(fid, ['\n<BR><U>Standard deviation of shape modes:</U> ', num2str(std(shape_modes(:)))]);
2087     fprintf(fid, ['\n<BR><U>Mean of intensity modes:</U> ', num2str(mean(intensity_modes(:)))]);
2088     fprintf(fid, ['\n<BR><U>Standard deviation of intensity modes:</U> ', num2str(std(intensity_modes(:)))]);
2089     fprintf(fid, ['\n<BR><U>Mean of shape variance:</U> ', num2str(mean(shape_variance(:)))]);
2090     fprintf(fid, ['\n<BR><U>Standard deviation of shape variance:</U> ', num2str(std(shape_variance(:)))]);
2091     fprintf(fid, ['\n<BR><U>Mean of intensity variance:</U> ', num2str(mean(intensity_variance(:)))]);   
2092     fprintf(fid, ['\n<BR><U>Standard deviation of intensity variance:</U> ', num2str(std(intensity_variance(:)))]);
2093     fprintf(fid, ['\n<BR><U>Mean time:</U> ', num2str(mean(time(:)))]);
2094     fprintf(fid, ['\n<BR><U>Standard deviation of time:</U> ', num2str(std(time(:)))]);    
2095     fprintf(fid, ['\n<BR><U>Total time:</U> ', num2str(t)]);
2096     fprintf(fid, ['\n</TD></TR>']);     
2097     fprintf(fid, '\n</TABLE></CENTER><BR><BR>');
2098     fclose(fid);
2099 end
2100 
2101 
2102 
2103 if (produce_record == 1),
2104      % In the GUI it has been renamed to "Display Record"
2105     clc;
2106     disp(' ');
2107     disp(' ');
2108     disp('===============================================');
2109     disp(['| Simulation Handle: ', [handle]]);
2110     disp(['| Description: ', [description]]);
2111     disp('|==============================================');
2112     disp(['| AART - ', [get_month], ' ', [get_year], ' - Version ', [get_version], ' - ', [getenv('OS')]]);
2113     disp(['| Run on machine: ', [get_host], '| Domain: ', [get_domain], '| Image width: ', num2str(image_width)]);
2114     disp('|==============================================');
2115     disp(' ');
2116     disp('|== GENERAL ===================================');
2117     disp(['| Objective Function: ', objective_function]);
2118     disp(['| Data type identifier: ', num2str(load_data)]); 
2119     disp(['| Iterations: ', num2str(n_iterations)]);
2120     disp(['| Images: ', num2str(n_images)]);
2121     disp(['| Sets: ', num2str(n_sets)]);
2122     disp(['| Corrective Offset (boolean): ', num2str(apply_intensity_offset)]);
2123     disp('|==============================================');
2124     disp(' ');
2125     disp('|== SMOOTHING =================================');
2126     disp(['| Smoothing window: ', num2str(window)]);
2127     disp(['| Gaussian smoothing (boolean): ', num2str(gaussian_smooth_enabled)]);  
2128     disp(['| Average smoothing (boolean): ', num2str(average_smooth_enabled)]);
2129     disp('|==============================================');
2130     disp(' ');
2131     disp('|== OBJECTIVE FUNCTION-SPECIFIC ===============');
2132      disp(['| Precision Required (All): ', num2str(precision_required)]);
2133     disp(['| Number of bins (MI, NMI): ', num2str(number_of_bins)]);
2134     disp(['| Number of model modes: ', num2str(n_modes)]);
2135     disp(['| Variation Kept (PCA): ', num2str(variation_kept)]);
2136     disp(['| Weighting normalisation method: ', weighting_normalisation_method]);
2137     disp(['| Contant weighting value: ', num2str(shape_weight)]);
2138     disp(['| PDF type: ', kernel]);
2139     disp(['| Wavelet filter: ', filter]);
2140     disp(['| Compression level: ', num2str(filter_level) ]); 
2141     disp(['| Model evaluation method: ',  num2str(model_evaluation_method)]);
2142     disp(['| Specificity iterations: ', num2str(spec_iters)  ]);
2143     disp(['| Generalisability iterations: ', num2str(gen_iters)  ]);
2144     disp(['| Evaluation cycle: ', num2str(cycle)  ]);    
2145     disp('|==============================================');
2146     disp(' ');
2147     disp('|== TRANSFORMATION ============================');
2148     disp(['| CPS Spline type: ', spline_type  ]);
2149     disp(['| Knot-point placement method: ', placement_type]);  
2150     disp(['| Number of knot-points: ', num2str(n_points)]);
2151     disp(['| Automatic precision (boolean): ', num2str(dynamic_precision)]);
2152     disp(['| Initialise near target (boolean): ', num2str(initialise_warps_close_to_target)]);
2153     disp(['| Offset extent: ', num2str(offset_extent)]);
2154     disp(['| Perturbation method: ', [perturbation_method]]);
2155     disp(['| Force reference (boolean): ', num2str(force_reference)]);
2156     disp(['| Reference type: ', num2str(force_first_reference)]);
2157     disp(['| Distance type: ', num2str(distance_type)]);    
2158     disp(['| Retain peak (boolean): ', num2str(retain_original_peak)]);    
2159     disp(['| Peak type: ', [peak_type]]);    
2160     disp(['| Stuff points (boolean): ', num2str(add_points)]);    
2161     disp(['| Point stuffing cycle: ', num2str(add_points_cycle)]);    
2162     disp('|==============================================');
2163     disp(' ');
2164     disp('|== DATA ======================================');
2165     disp(['| Base height variation: ', num2str(bump_height)  ]);
2166     disp(['| Height variation: ', num2str(bump_height_variation) ]);
2167     disp(['| Base width variation: ', num2str(bump_width)  ]);
2168     disp(['| Width variation: ', num2str(bump_width_variation) ]);
2169     disp(['| Position Freedom: ', num2str(bump_position_freedom)  ]);
2170     disp(['| Smoothness: ', num2str(smoothness_factor)]);
2171     disp('|==============================================');
2172     disp(' ');
2173     disp('|==============================================');
2174     disp('| Results and outputs to be located by handle |');
2175     disp('===============================================');    
2176        
2177 end  
2178 
2179 if (save_in_log == 1),
2180      % if the parameters need to be included in the log file
2181     
2182 %     fid_test = fopen('aartlo.htm','r');
2183 %     f = fread(fid_test,1);
2184 %     c = char (f')
2185 %     a=2
2186 %     fclose(fid_test);
2187     if (save_movie == 1),
2188       movie_name = [[handle], '-1'];
2189     end 
2190     % modify the name of the movie to fit the HTML link
2191     % if multiple movies are added in the future, then the HTML structure
2192     % needs to fit the extra links (slots)
2193 
2194     if (strcmp(which('aartlog.htm'), '')),
2195         add_headers = 1;
2196     else
2197         add_headers = 0;
2198     end
2199 
2200     
2201     
2202     fid = fopen('aartlog.htm','a');
2203     if (add_headers == 1),
2204       add_html_headers(fid);
2205     end
2206     
2207     fprintf(fid, ['\n<BR><CENTER><A HREF="index.htm"><IMG src="../../back.jpg" border="0" ALT="Back to sub-index"></A></CENTER><BR><BR><A NAME="', [handle],'"><HR SIZE=5></A><H1>Simulation Handle: ', [handle], '</H1>']);
2208     fprintf(fid, ['\n<BR>', '<CENTER><TABLE BORDER=BOX CELLPADDING=10 CELLSPACING=3 BACKGROUND="../../scbg.gif"><TR><TD>']);
2209     fprintf(fid, ['\n<H2>Main</H2>']);
2210     
2211     if (strcmp(getenv('OS'), 'Linux')), 
2212       fprintf(fid, ['<BR>AART - ', [get_month], ' ', [get_year], ' - Version ', [get_version], ' - ', [getenv('OS')]]);
2213       fprintf(fid, ['\n<BR><U>Time</U>: ', [get_date]]); 
2214       fprintf(fid, ['\n<BR><U>Run on machine:</U> ', [get_host], '<BR><U>Domain:</U> ', [get_domain], '<BR><U>Image width:</U> ', num2str(image_width)]);
2215     end
2216     
2217     fprintf(fid, ['\n<BR><U>Description:</U> ', [description]], '</TD>');
2218     fprintf(fid, ['\n<TD><H2>General</H2>']);
2219     fprintf(fid, ['\n<BR><U>Objective Function:</U> ', objective_function]);
2220     fprintf(fid, ['\n<BR><U>Data type identifier:</U> ', num2str(load_data)]);    
2221     fprintf(fid, ['\n<BR><U>Iterations:</U> ', num2str(n_iterations)]);
2222     fprintf(fid, ['\n<BR><U>Images:</U> ', num2str(n_images)]);
2223     fprintf(fid, ['\n<BR><U>Sets:</U> ', num2str(n_sets)]);
2224     fprintf(fid, ['\n<BR><U>Corrective Offset (boolean):</U> ', num2str(apply_intensity_offset), '</TD></TR><TR><TD>']);
2225     
2226     fprintf(fid, ['\n<BR><BR><H2>Smoothing</H2>']);
2227     fprintf(fid, ['\n<BR><U>Smoothing window:</U> ', num2str(window)]);
2228     fprintf(fid, ['\n<BR><U>Gaussian smoothing (boolean):</U> ', num2str(gaussian_smooth_enabled)]);
2229     fprintf(fid, ['\n<BR><U>Average smoothing (boolean):</U> ', num2str(average_smooth_enabled)]);
2230     
2231     fprintf(fid, ['\n<BR><BR><H2>Objective Function-Specific</H2>']);
2232     fprintf(fid, ['\n<BR><U>Precision Required (All):</U> ', num2str(precision_required)]);
2233     fprintf(fid, ['\n<BR><U>Number of bins (MI, NMI):</U> ', num2str(number_of_bins)]);
2234     fprintf(fid, ['\n<BR><U>Number of model modes:</U> ', num2str(n_modes)]);
2235     fprintf(fid, ['\n<BR><U>Variation Kept (PCA):</U> ', num2str(variation_kept)]);
2236     fprintf(fid, ['\n<BR><U>Weighting normalisation method:</U> ', weighting_normalisation_method]);
2237     fprintf(fid, ['\n<BR><U>Contant weighting value:</U> ', num2str(shape_weight)]);
2238     fprintf(fid, ['\n<BR><U>PDF type:</U> ', kernel]);
2239     fprintf(fid, ['\n<BR><U>Wavelet filter:</U> ', filter]);
2240     fprintf(fid, ['\n<BR><U>Compression level:</U> ', num2str(filter_level)]);
2241     fprintf(fid, ['\n<BR><U>Model evaluation method:</U> ',  num2str(model_evaluation_method)]);
2242     fprintf(fid, ['\n<BR><U>Specificity iterations:</U> ', num2str(spec_iters)]);
2243     fprintf(fid, ['\n<BR><U>Generalisability iterations:</U> ', num2str(gen_iters), '</TD><TD>']);
2244     
2245     fprintf(fid, ['\n<H2>Transformation</H2>']);
2246     fprintf(fid, ['\n<BR><U>Evaluation cycle:</U> ', num2str(cycle)]);  
2247     fprintf(fid, ['\n<BR><U>CPS Spline type:</U> ', spline_type]);
2248     fprintf(fid, ['\n<BR><U>Knot-point placement method:</U> ', placement_type]);
2249     fprintf(fid, ['\n<BR><U>Number of knot-points:</U> ', num2str(n_points)]);
2250     fprintf(fid, ['\n<BR><U>Automatic precision (boolean):</U> ', num2str(dynamic_precision)]);
2251     fprintf(fid, ['\n<BR><U>Initialise near target (boolean):</U> ', num2str(initialise_warps_close_to_target)]);
2252     fprintf(fid, ['\n<BR><U>Offset extent:</U> ', num2str(offset_extent)]);
2253     fprintf(fid, ['\n<BR><U>Perturbation method:</U> ', [perturbation_method]]);
2254     fprintf(fid, ['\n<BR><U>Force reference (boolean):</U> ', num2str(force_reference)]);
2255     fprintf(fid, ['\n<BR><U>Reference type:</U> ', num2str(force_first_reference)]);
2256     fprintf(fid, ['\n<BR><U>Distance type:</U> ', num2str(distance_type)]);    
2257     fprintf(fid, ['\n<BR><U>Retain peak (boolean):</U> ', num2str(retain_original_peak)]);    
2258     fprintf(fid, ['\n<BR><U>Peak type:</U> ', [peak_type]]);
2259     fprintf(fid, ['\n<BR><U>Stuff points (boolean):</U> ', num2str(add_points)]);    
2260     fprintf(fid, ['\n<BR><U>Point stuffing cycle:</U> ', num2str(add_points_cycle)]);  
2261     
2262     fprintf(fid, ['\n<BR><BR><H2>Data</H2>']);    
2263     fprintf(fid, ['\n<BR><U>Base height variation:</U> ', num2str(bump_height)]);
2264     fprintf(fid, ['\n<BR><U>Height variation:</U> ', num2str(bump_height_variation)]);
2265     fprintf(fid, ['\n<BR><U>Base width variation:</U> ', num2str(bump_width)]);
2266     fprintf(fid, ['\n<BR><U>Width variation:</U> ', num2str(bump_width_variation)]);    
2267     fprintf(fid, ['\n<BR><U>Position Freedom:</U> ', num2str(bump_position_freedom)]);
2268     fprintf(fid, ['\n<BR><U>Smoothness:</U> ', num2str(smoothness_factor)]);
2269     fprintf(fid, ['\n</TD></TR>']);
2270     
2271     fprintf(fid, '\n<TR><TD><H2>Image Preview 1</H2><BR><A HREF="');
2272     fprintf(fid, [[handle],'-1.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
2273     fprintf(fid, [[handle],'-1.jpg"></A></TD><TD><H2>Image Preview 2</H2><BR><A HREF="']);
2274     fprintf(fid, [[handle],'-2.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
2275     fprintf(fid, [[handle], '-2.jpg"></A></TD></TR>']);
2276     fprintf(fid, '\n<TR><TD><H2>Image Preview 3</H2><BR><A HREF="');
2277     fprintf(fid, [[handle],'-3.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
2278     fprintf(fid, [[handle],'-3.jpg"></A></TD><TD><H2>Image Preview 4</H2><BR><A HREF="']);
2279     fprintf(fid, [[handle],'-4.jpg"><IMG HEIGHT=200 WIDTH=300 BORDER=0 SRC="']);
2280     fprintf(fid, [[handle], '-4.jpg"></A></TD></TR>']);    
2281     fprintf(fid, '\n<TR><TD><A HREF="');
2282     fprintf(fid, [[handle],'-1.avi">Video #1</A></TD><TD><A HREF="']);
2283     fprintf(fid, [[handle],'-2.avi">Video #2</A></TD></TR>']); 
2284     fprintf(fid, '\n<TR><TD><A HREF="');
2285     fprintf(fid, [[handle],'.mat">Input Data</A></TD><TD><A HREF="']);
2286     fprintf(fid, [[handle],'.htm">Statistics</A></TD></TR>']);     
2287     fprintf(fid, '\n</TABLE></CENTER><BR><BR>');
2288     fclose(fid);
2289     
2290       % save in index file
2291 
2292     fid = fopen('index.htm','a');
2293     if (add_headers == 1),
2294       add_html_headers(fid);
2295       add_table_index(fid);
2296     end
2297     if (strcmp(getenv('OS'), 'Linux')),
2298       fprintf(fid, ['\n<P><CENTER><TABLE WIDTH=100%% BORDER=BOX CELLPADDING=10 CELLSPACING=3 BACKGROUND="../../scbg.gif"><TD ALIGN=CENTER WIDTH=33%%><A HREF="aartlog.htm#', [handle], '">', [handle], '</A></TD><TD ALIGN=CENTER WIDTH=33%%>', [get_date], '</TD><TD ALIGN=CENTER WIDTH=33%%>']);
2299       fprintf(fid, ['\n', [description], '</P></TD></TABLE></CENTER>']);
2300     else
2301       fprintf(fid, ['\n<P><CENTER><TABLE WIDTH=100%% BORDER=BOX CELLPADDING=10 CELLSPACING=3 BACKGROUND="../../scbg.gif"><TD ALIGN=CENTER WIDTH=33%%><A HREF="aartlog.htm#', [handle], '">', [handle], '</A> ', '</TD><TD ALIGN=CENTER WIDTH=33%%>']);
2302       fprintf(fid, ['\n', [description], '</P></TD></TABLE></CENTER>']);
2303     end
2304     fclose(fid);
2305 end
2306 
2307 % a=figure('Name','test');
2308 % b=figure(a);
2309 % hold on;
2310 % plot([-1,2,3]);
2311 % hold off;
2312 % saveas(b,[[handle],'-', num2str(figure_handle_number)], 'jpg');
2313 % figure_handle_number = figure_handle_number + 1;
2314 
2315 
2316 %% END OF STATISTICS                                  %%
2317 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2318 
2319 
2320 if (save_movie == 1),
2321   if (strcmp(progress_indicator_type,'Progress Bar')),
2322            h = progressbar( h, 0, 'Saving Video');
2323   end  
2324   movie(M);
2325   movie2avi(M, movie_name, 'FPS', frames_per_second);
2326   if (called_from_gui == 1),
2327       % send this to the GUI status bar
2328       status = [[status];['(**) Video saved      '];['                      ']];
2329   else
2330       % indicate to user that a video has been save
2331         disp(' ');
2332         disp([movie_name, ' saved']);
2333           % process the movie and save it
2334   end
2335 end
2336 
2337 if (strcmp(progress_indicator_type,'Progress Bar')),
2338    progressbar( h,-1 );
2339       % destroy progress bar
2340 end
2341 
2342 if (error_found == 0),
2343               % if the error flag is not set
2344     if (called_from_gui == 1),
2345       status = [[status];['Registration Completed'];
2346              ['                      '];       
2347              ['Status: OK            ']];
2348          % send message
2349     else
2350       status = 'Status: OK'; 
2351          % indicate that everything was okay
2352     end
2353 elseif (error_found == 1),
2354     if (called_from_gui == 1),
2355       status = [[status];['Registration Completed'];
2356              ['                      '];       
2357              ['Status: Error         '];
2358              ['                      '];
2359              ['Operation Aborted     ']];
2360          % send message about error
2361     else
2362       status = 'Status: Operation Aborted'; 
2363          % indicate that something had gone wrong
2364     end    
2365 else
2366     if (called_from_gui == 1),
2367       status = [[status];['Registration Completed'];
2368              ['                      '];       
2369              ['Status: Error         '];
2370              ['                      '];
2371              ['Exit code unknown     ']];
2372          % send message about error
2373     else
2374       status = 'Status: Error'; 
2375          % indicate that something had gone wrong
2376     end
2377 end    
2378 
2379 if (exit_at_end == 1),
2380     exit;
2381 end  
2382 
2383 if (close_at_end == 1),
2384     close;
2385 end
build_1d_bump_model

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
