function build_1d_model(spline_type, placement_type, objective_function)
% BUILD_1D_MODEL: Builds a 1-D model and returns some related statistics.
%                 This is the main routine to be called in the code.
%
% Code written by Katherine Smith, 2003
%
%    GENERAL
%
%      build_1d_model(spline_type, placement_type, objective_function)
%
%    INPUT/S
%
%      -spline_type: 
%           The type of spline to be used. 'single_point' or
%           'multi_point' at present.
%
%      -placement_type:
%           Placement type for 'multi_point':
%                    'grid' or 'edges' or 'random'
%           placement type for 'single_point':
%                    'overlap_grid' or 'edges_and_scale'
%                    or 'random_and_scale'.
%
%      -objective_function:
%           'msd_opt_separate' or 'model_opt_together'.
%           
%    OUTPUT/S
%
%      -None. Figures produced and results displayed are main
%       data to be utilised.
%
%    PENDING WORK
%
%      -Improve readability of file.
%      -Add comments
%
%    KNOWN BUG/S
%
%      -The algorithm is known to be very slow. Some profiling
%       could possibly be used to hasten the process.
%
%    COMMENT/S
%
%      -This function needs some work. It is messy in my opinion.
%      -Smith: the function makes images and points
%      -Many comments added 27/11/03
%      -28/11/03: more comments added.
%
%    RELATED FUNCTION/S
%
%      See other functions in same directory, filenames
%      beginning with "build" in particular. There seem to be much
%      overlap due to cut-and-paste.
%
%    ABOUT
%
%      -Created:     November 23rd, 2003
%      -Last update: November 28th, 2003
%      -Revision:    0.0.7
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% variables (parameters) initialisation

frame_counter = 0;
n_images = 10;
         % how many images to try to register in each set, default was 10
n_iterations = 10;
         % how many iterations of the optimisation to run, default was 50.
n_sets = 1;
         % [used to be 25]
	 % how many sets of images to register
	 % (in order to get better statistics).
image_width = 50;
         % spline type 'single_point' or 'multi_point'
spline_type = 'multi_point';
         % placement type for multipoint: 'grid' or 'edges' or 'random'
         % placement type for singlepoint: 'overlap_grid' or 'edges_and_scale' or 'random_and_scale'
placement_type = 'grid';
objective_function = 'msd_opt_together';
verbose = 'off'; 
         % be verbose or not
n_points = 5;
         % an argument that is passed to the model optimisation function.
         % meaning still unknown.
do_plot = 1; 
         % should model be plotted as it changes? (Boolean)

white_ctr = 0;
         % white counter?
ref_shift = 0.2;
       % shift in reference image??
max_shift = 0.2;
       % maximum allowed shift?
step = 0.1;
       % used below only
los = zeros(n_images,floor((max_shift-ref_shift)/step));
his = zeros(n_images,floor((max_shift-ref_shift)/step));
       % not understood yet
blurred = 0;
       % image blurring enabled???
spec_iters = 25;
       % specificity iterations
gen_iters = 25;
       % generalisability iterations
min_error = 0;
max_error = 1;
       % define error range??
error_step = 0.1;
n_error_steps = (max_error-min_error)/error_step;

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

tic;
       % start counting time
for this_set = 1:n_sets
       % for all sets
   tic;
       % why is this called again??
       % Smith: for each set of n_sets images, build a model
       % calculate obj fn values from these models, plot the mean with error
       % bars given by 1 stddev away
   ['calculating set ', num2str(this_set)]
       % user gets status
   [imagelist images points his los] = make_1d_images(n_images, image_width, 0.2);
       % create a set of images
% window = 9;
% images = average_smooth(images, window); blurred = 1;
% images = gaussian_smooth(images, window); blurred = 1;
       % smooth all images - currently disabled

   points = -1 + 2*(points-1)/(image_width-1);
       % Smith: normalise points from -1 to 1
   keep = 0.999999;
       % is this the precision required???
   ref_points_vec = points(:,1);
   ref_image_vec = images(:,1);
       % set first image generated to be reference
   ref_hi = his(1);
   ref_lo = los(1);
       % and get its upper and lower boundaries

  if(do_plot)
    subplot_fig = figure;
    H=figure(subplot_fig);
    for i=1:n_images,
       subplot(n_images,2,(2*i)-1);
       plot(images(:,i)),title('Unwarped images');
   end
  end
%% NOTE: Change above in traversal (RSS 26/11/03)

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

  warped_images = images;
  warped_points = points;
      % set this initial state to plot the examples before change


% correctly_warped_points = points;
% correctly_warped_images = images;
% for i=1:n_images
% % now bung in the right deformations to get them back
% correctly_warped_points(:,i) = linear_warp(ref_points_vec, los(i,white_ctr), his(i,white_ctr), ref_lo, ref_hi);
% correctly_warped_points_vec = correctly_warped_points(:,i);
% image_vec = images(:,i);
% % resample at the warped points - point on reg. grid has value of point at same index in warpy grid
% % need to interpolate
% correctly_warped_images(:,i) = interp1(ref_points_vec, image_vec, correctly_warped_points_vec);
% correct_model = bufild_model(correctly_warped_images,correctly_warped_points,1,'','edge',0);
% converging_score(i) = measure_model(correct_model.variances,50);
% end
% figure,plot(log(converging_score)),title('Log of score as model converges');


  for n=1:n_iterations
       % iterations aim to get good statistical results by averaging
    disp(['iter ',num2str(n)]);
       % Smith: correct registration
       % attempt to register
    for i=1:n_images
     disp(['Warping image ',num2str(i),' of ',num2str(n_images)]);
     image_vec = warped_images(:,i);
     points_vec = warped_points(:,i);
     if(strcmp(objective_function,'model_opt_together'))
         if(do_plot)
             frame_counter = frame_counter + 1;
             M(frame_counter) = getframe(H);    
             [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*2);
         else
             [param_list, warped_point, warped_image, score] = optimise_all_warps_model([warped_images(:,1:i-1),warped_images(:,i+1:n_images)], [warped_points(:,1:i-1),warped_points(:,i+1:n_images)], image_vec, points_vec,spline_type,placement_type,n_points, verbose, do_plot);
         end
     elseif(strcmp(objective_function,'model_opt_separate'))
         if(do_plot)
             frame_counter = frame_counter + 1;
             M(frame_counter) = getframe(H);    
             [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*2);
         else
             [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);
         end
     elseif(strcmp(objective_function,'msd_opt_together'))
         if(do_plot)
             frame_counter = frame_counter + 1;
             M(frame_counter) = getframe(H);    
             [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*2);
         else
             [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);
         end
     elseif(strcmp(objective_function,'msd_opt_separate'))
         if(do_plot)
             [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*2);
             frame_counter = frame_counter + 1;
             M(frame_counter) = getframe(H);    
         else
             [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);
         end 
     else
         error(['Unknown objective function: ', objective_function]);
    end

    warped_images(:,i) = warped_image;
    warped_points(:,i) = warped_point;

     % figure(subplot_fig),title(['Warped images, after image ',num2str(i),' iteration ',num2str(n)]);
     % dummy = waitforbuttonpress;
     if(n==n_iterations),
        final_score(i, this_set) = score;
     end
   end
        % end:images
   % filename = ['ex1num' num2str(n) '.png']; 
          % set name of file
   %imwrite(warped_images(:,3), filename, 'png');
  end  
        % end:iteration  


% Smith: build model from warped images

  w_c_model = build_model(warped_images, warped_points, keep,'Optimised warp', 'variance');
      % This is where a model is being built for the warped images
% [w_c_model.mean_specificity, w_c_model.std_specificity] = find_specificity(w_c_model, images, spec_iters, ref_points_vec);
% show_combined_model(w_c_model,ref_points_vec, 2, 0.2, 'Combined model built from optimised warps');



  fig_title = 'Automatically aligned: ';
% show_shape_model(w_shape_model, ref_points_vec, ref_image_vec, 2, white_width)
% show_intensity_model(w_intensity_model, 2, white_width, fig_title)
% show_combined_model(w_c_model, ref_points_vec, 2, white_width, fig_title)

  w_intensity_total_vars = w_c_model.intensity_model.total_var;
  w_shape_total_vars = w_c_model.shape_model.total_var;
     % get these two variances to be used later to estimate ''goodness'
     % of the model, possibly using determinant of covariance matrix

  model_score(this_set) = measure_model(w_c_model.variances, 20);
  msd_score(this_set) = measure_model_msd(warped_images);
  shape_modes(this_set) = w_c_model.n_shape_modes;
  intensity_modes(this_set) = size(w_c_model.intensity_model.pcs,2);
  shape_variance(this_set) = sum(w_c_model.shape_model.variances);
  intensity_variance(this_set) = sum(w_c_model.intensity_model.variances);

  t = toc;
  disp(['Time for set ',num2str(this_set),': ',num2str(t)]);
  time(this_set)=t;
  % Get and show time statistics
  
end % end set


% END OF MODEL-BUILDING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now provide some statistics on the simulation/experiment




disp('Mean match score:');
mean(final_score(:))
disp('Std match score:');
std(final_score(:))

disp('Mean model score:');
mean(model_score(:))
disp('Std model score:');
std(model_score(:))

disp('Mean msd score:');
mean(msd_score(:))
disp('Std msd score:');
std(msd_score(:))

disp('Mean shape modes:');
mean(shape_modes(:))
disp('Std shape modes:');
std(shape_modes(:))

disp('Mean intensity modes:');
mean(intensity_modes(:))
disp('Std intensity modes:');
std(intensity_modes(:))

disp('Mean shape variance:');
mean(shape_variance(:))
disp('Std shape variance:');
std(shape_variance(:))

disp('Mean intensity variance:');
mean(intensity_variance(:))
disp('Std intensity variance:');
std(intensity_variance(:))

disp('Mean time: ');
mean(time(:))
disp('Std time: ');
std(time(:))

t = toc;   
       % get time from the counter initiated by <tic>
disp(['Total time: ',num2str(t)]);  
       % display total time
       
movie(M);
movie2avi(M,'set10.avi');

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