function [param_list, warped_points, warped_image, score] = optimise_all_warps_model(ref_image, ref_points, unwarped_image, unwarped_points, spline_type, placement_method, n_points, verbose, do_plot, subplot_fig, pane, column, image_width, number_of_images, verbose_score, model_evaluation_method, n_modes, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising, precision_required); %use_parameters, parameters)
% OPTIMISE_ALL_WARPS_MDOEL: 
%
% Code written by Katherine Smith, 2003
%
%    GENERAL
%
%      [param_list, warped_points, warped_image, score] = 
%          optimise_all_warps_model(ref_image, ref_points, unwarped_image, 
%          unwarped_points, spline_type, placement_method, n_points, verbose, 
%          do_plot, subplot_fig, pane,
%          image_width, number_of_images, verbose_score)
% 
%    INPUT/S
%
%      -X:
%          X
%           
%    OUTPUT/S
%
%      -X:
%           X
%
%    PENDING WORK
%      
%      
%
%    KNOWN BUG/S
%
%      
%
%    COMMENT/S
%
%      
%
%    RELATED FUNCTION/S
%
%      
%
%    ABOUT
%
%      -Created:     November 23rd, 2003
%      -Last update: January 23rd,  2004
%      -Revision:    0.1.5
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================


% figure(3);
% hold on
% plot(unwarped_image)
% hold off     ** Image is dynamic
[grid, steps] = setup_grid(2, unwarped_image, ref_points, placement_method, n_points);

% POSSIBLE OPTIMISATION: where no randomisation is applied for grid
% selection, e.g. 'overlap grid' and 'grid', the grid needs to be set up ONCE, not every
% iteration. A test for this needs to be put outside this function to avoid
% unnecessary computations.

if (do_plot),
  figure(subplot_fig);
  subplot(number_of_images, column ,pane , 'replace');
  plot(unwarped_image);
  axis([0 image_width 0 1]);
  drawnow;
end

warp_params.green = 'biharmCPS';


% added 16/03/04
% if (~use_parameters),
%   param_list.d = [0,0,0,0,0];
% else
%   for i=1:size(grid),
%     param_list.d(i) = parameters.d(i);
%   end
% end  
% if (use_parameters),
% parameters.d
% end


% 16/03/04:param_list = [];
% also remove 2 last arguments of this function

param_list = [];
start_image = unwarped_image;
start_points = unwarped_points;

%last_score = msd(ref_image, unwarped_image);


if (strcmp(spline_type,'single_point')),

    
      %     steps
      %     grid
      %     size(grid)
  

  last_score = eval_model_cps_warp(zeros(size(grid)), grid, steps, unwarped_image, unwarped_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);  
      % last_score
  
      % if (~use_parameters),
   
  
[params, score] = fminsearch(@eval_model_cps_warp, zeros(size(grid),1), optimset('Display', verbose, 'TolX', precision_required) , grid, steps, start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);     


% 'MaxIter', 1000
% str2num([['1e-'],[num2str(ceil(rand*100))]])
%   score
%   params=0.1;
%   [score1,start_points] = eval_model_cps_warp(params, grid, steps, start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
%   score1
%   params=0.2;
%   [score2,start_points] = eval_model_cps_warp(params, grid, steps, start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
%   score2
%   params=0.2;
%   [score2,start_points] = eval_model_cps_warp(params, grid, steps, start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
%   score2
%   params=0.3;
%   [score3,start_points] = eval_model_cps_warp(params, grid, steps, start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
%   score3
%   params=0.4;
%   [score4,start_points] = eval_model_cps_warp(params, grid, steps, start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
%   score4
%   params=0.5;
%   [score5,start_points] = eval_model_cps_warp(params, grid, steps, start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
%   score5
      % else 
      % [params, score] = fminsearch(@eval_model_cps_warp,parameters.d,optimset('Display',verbose),grid,steps,start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept);
      %  end
      % try randn(size(grid),1)
  param_list(size(param_list,1)+1).d = params;
  [dummy,start_points] = eval_model_cps_warp(params, grid, steps, start_image, start_points, ref_image, ref_points, n_modes, start_points, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
  start_image = interp1(unwarped_points, unwarped_image, start_points, 'linear',0); 
  
  
  
  
elseif (strcmp(spline_type,'multi_point')),
  warped_grid = grid;
  last_score = eval_model_multi_warp(zeros(size(grid)), grid, unwarped_image, unwarped_points, ref_image, ref_points, n_modes, start_points, warp_params, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
  [params, score] = fminsearch(@eval_model_multi_warp, zeros(size(grid)), optimset('Display', verbose, 'TolX', precision_required), grid, unwarped_image, unwarped_points, ref_image, ref_points, n_modes, start_points, warp_params, model_evaluation_method, weighting_normalisation_method, shape_weight, variation_kept, show_warps_while_optimising);
  param_list(size(param_list,1)+1,:).d = params;
  warped_grid = grid + params;
  start_points = nrr_trans_1d(start_points, grid, warped_grid, warp_params,[]);
  start_image = interp1(unwarped_points,unwarped_image,start_points, 'linear',0);
  
  
else
  error('Unknown spline type');  
end




if (verbose_score == 1),
    disp(['Initial score: ',num2str(last_score),' Final score: ', num2str(score), ' Difference: ', num2str(last_score - score)]);
end

% NOTE: It is not yet clear to me if parameters that are recorded would be
% a problem here
if (last_score > score),
  warped_points = start_points;
  warped_image = start_image;
else
  warped_points = unwarped_points;
  warped_image = unwarped_image;    
end    

if (do_plot),
  figure(subplot_fig);
  subplot(number_of_images, column , pane, 'replace');
  plot(warped_image);
  axis([0 image_width 0 1]);
  title(num2str(score));
  drawnow;
end
