function [] = evaluate_model_from_synthesis_and_training_set( ...
			training_set_path, ...
			synthesis_set_path, ... 
			type_training, ...
			type_synthesis, ...
			show_images, ... 
			smoothing_filter_size, ... 
			filter_type, ...
			smooth_all_images, ...
			smooth_difference_image, ... 
			scale_brightness, ... 
			calculate_fast_shuffle_distance_circle, ... 
			shuffle_radius, ...
			squared_differences, ... 
			distance_of_power_2_in_difference_image, ... 
			calculate_shuffle_distance_square, ...
			shuffle_window, ... 
			resize_training_set, ...
			resize_synthesis_set, ...
			resize_width, ... 
			resize_height, ... 
			focus, ...
			threshold, ... 
			crop_frame_from_synthesis, ... 
            crop_width, ...
			crop_boundaries_from_training, ...
			crop_width_training, ...
            save_filename, ...
            make_shuffle_square_symmetric, ...
            show_distance_matrix_at_end, ...
            use_histogram, ...
            histogram_window, ...
            histogram_bins, ...
            smooth_histograms, ...
            histogram_smoothing_filter_size)
        
% VALUATE_MODEL_FROM_SYNTHESIS_AND_TRAINING_SET: Evaluate model from some synthesis of examples % derived from it and its training set.
%
%
%    GENERAL
%
%      
%
%    INPUT/S
%
%      -Omitted
%
%    OUTPUT/S
%
%      -Omitted
%
%    PENDING WORK
%
%      -
%
%    KNOWN BUG/S
%
%      -None.
%
%    COMMENT/S
%
%      -
%
%    RELATED FUNCTION/S
%
%      
%
%    ABOUT
%
%      -Created:     February 2005
%      -Last update: March 15th, 2005
%      -Revision:    0.0.1
%      -Author:      R. S. Schestowitz, University of Manchester
% ==============================================================	    
	    
	   
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% Start

warning off;
  % just in case

addpath(genpath('/home/S00/schestr0/Model_Evaluation'));
addpath(genpath('/home/S00/schestr0/AART'));

cd(training_set_path);
listing_training = dir(type_training);

% progress bar variables

progress_bar_position = 0;
time_for_this_iteration = 10;
maximum = size(listing_training,1);

for k=1:size(listing_training,1),
   tic;
       % skip the dots
   if (listing_training(k).isdir == 0),
      cd(training_set_path);
      image1=imread(listing_training(k).name);
      if (crop_boundaries_from_training),
        image2=crop_frame(image1, crop_width_training);
      end
      if (focus),
        image1=focus_brain(image1, threshold); % crop brain
      end
      if (resize_training_set),
        image1=imresize(image1,[resize_height ,resize_width],'bicubic');
          % make images the same size by downsizing the training set
      end  
      if (smoothing_filter_size>1 & smooth_all_images),
         image1 = smoothn(image1, smoothing_filter_size, filter_type);
      end 
   end
   cd(synthesis_set_path);
   listing_synthesis = dir(type_synthesis);
   for l=1:size(listing_synthesis,1),
      image2=imread(listing_synthesis(l).name);
      if (crop_frame_from_synthesis),
        image2=crop_frame(image2, crop_width);
      end
      if (focus),
        image2=focus_brain(image2, threshold); % crop brain
      end
      if (resize_synthesis_set),
        image2=imresize(image2,[resize_height ,resize_width],'bicubic');
      end
      
      if (smoothing_filter_size>1 & smooth_all_images),
         image2 = smoothn(image2, smoothing_filter_size, filter_type);
      end    
      
      

      if (calculate_fast_shuffle_distance_circle),
          % Marsland's implementation:
          % difference_image = shuffle_transform(double(image2), double(image1), shuffle_radius);
          DistanceMatrix(k,l) = iman_shuffle_euclid(double(image1),double(image2),shuffle_radius); % Petrovic
      else % slower and more flexible measures go here, ones that scan pixel by pixel naively
          for i=1:size(image2,1),
            for j=1:size(image2,2), % Assumes both images are of same size at this stage
                if (calculate_shuffle_distance_square),
                     if (shuffle_window==3),
                        if (i == 1 || i==size(image2,1) || j == 1 || j==size(image2,2)),
                           difference_image(i,j) = 0;
                        else % currently shuffle with window of size 3x3, implementation that is general (arbitrary size) is not so trivial to implement
                         difference_image(i,j) = uint8(min( ...
                         [abs(double(image2(i,j)) - double(image1(i,j))), ...
                         abs(double(image2(i,j)) - double(image1(i,j-1))), ...
                         abs(double(image2(i,j)) - double(image1(i,j+1))), ...
                         abs(double(image2(i,j)) - double(image1(i-1,j))), ...
                         abs(double(image2(i,j)) - double(image1(i+1,j))), ...
                         abs(double(image2(i,j)) - double(image1(i-1,j-1))), ...
                         abs(double(image2(i,j)) - double(image1(i+1,j-1))), ...
                         abs(double(image2(i,j)) - double(image1(i-1,j+1))), ...
                         abs(double(image2(i,j)) - double(image1(i+1,j+1)))] ...
                         ));
                        end
                     elseif (shuffle_window == 5),
                        if (i == 1 || i == 2 || i == size(image2,1) || i == size(image2,1) - 1 ...
                            || j == 1 || j == 2 || j == size(image2,2) || j == size(image2,2) - 1), 
                           difference_image(i,j) = 0;
                        else

                         if (make_shuffle_square_symmetric),
                            difference_image(i,j) = mean([ ...
                                 uint8(min( ...
                                 [abs(double(image2(i,j)) - double(image1(i,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j-2))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j-2))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j-2))), ...
                                 abs(double(image2(i,j)) - double(image1(i,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i,j-2))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j-2)))] ...
                                 )), ...
                                 uint8(min( ...
                                 [abs(double(image1(i,j)) - double(image2(i,j))), ...
                                 abs(double(image1(i,j)) - double(image2(i,j-1))), ...
                                 abs(double(image1(i,j)) - double(image2(i,j+1))), ...
                                 abs(double(image1(i,j)) - double(image2(i-1,j))), ...
                                 abs(double(image1(i,j)) - double(image2(i+1,j))), ...
                                 abs(double(image1(i,j)) - double(image2(i-1,j-1))), ...
                                 abs(double(image1(i,j)) - double(image2(i+1,j-1))), ...
                                 abs(double(image1(i,j)) - double(image2(i-1,j+1))), ...
                                 abs(double(image1(i,j)) - double(image2(i+1,j+1))), ...
                                 abs(double(image1(i,j)) - double(image2(i-2,j+2))), ...
                                 abs(double(image1(i,j)) - double(image2(i-2,j+1))), ...
                                 abs(double(image1(i,j)) - double(image2(i-2,j))), ...
                                 abs(double(image1(i,j)) - double(image2(i-2,j-1))), ...
                                 abs(double(image1(i,j)) - double(image2(i-2,j-2))), ...
                                 abs(double(image1(i,j)) - double(image2(i+2,j+2))), ...
                                 abs(double(image1(i,j)) - double(image2(i+2,j+1))), ...
                                 abs(double(image1(i,j)) - double(image2(i+2,j))), ...
                                 abs(double(image1(i,j)) - double(image2(i+2,j-1))), ...
                                 abs(double(image1(i,j)) - double(image2(i+2,j-2))), ...
                                 abs(double(image1(i,j)) - double(image2(i+1,j+2))), ...
                                 abs(double(image1(i,j)) - double(image2(i+1,j-2))), ...
                                 abs(double(image1(i,j)) - double(image2(i,j+2))), ...
                                 abs(double(image1(i,j)) - double(image2(i,j-2))), ...
                                 abs(double(image1(i,j)) - double(image2(i-1,j+2))), ...
                                 abs(double(image1(i,j)) - double(image2(i-1,j-2)))]))] ...            
                                 );
                             
                         else
                           difference_image(i,j) = uint8(min( ...
                                 [abs(double(image2(i,j)) - double(image1(i,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i-2,j-2))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j+1))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j-1))), ...
                                 abs(double(image2(i,j)) - double(image1(i+2,j-2))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i+1,j-2))), ...
                                 abs(double(image2(i,j)) - double(image1(i,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i,j-2))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j+2))), ...
                                 abs(double(image2(i,j)) - double(image1(i-1,j-2)))] ...
                                 ));
                         end
                        end
                     else
                         disp('Invalid window size for naive shuffling implementation');
                     end
                     
                     if (squared_differences),
                       difference_image(i,j) = difference_image(i,j) * difference_image(i,j);
                     end
	
                     
                     
                elseif (use_histogram), % histogram-based (locally-orderless images)
                   
                   if (histogram_window==3),
                        if (i == 1 || i==size(image2,1) || j == 1 || j==size(image2,2)),
                           difference_image(i,j) = 0;
                        else % currently shuffle with window of size 3x3, implementation that is genereal is not so trivial
                         histogram1 = hist([double(image1(i,j)), double(image1(i,j-1)), double(image1(i,j+1)), double(image1(i-1,j)), double(image1(i+1,j)), ...
                         double(image1(i-1,j-1)), double(image1(i+1,j-1)), double(image1(i-1,j+1)), double(image1(i+1,j+1)), 255, 0], histogram_bins);
                 
                         histogram2 = hist([double(image2(i,j)), double(image2(i,j-1)), double(image2(i,j+1)), double(image2(i-1,j)), double(image2(i+1,j)), ...
                         double(image2(i-1,j-1)), double(image2(i+1,j-1)), double(image2(i-1,j+1)), double(image2(i+1,j+1)), 255, 0], histogram_bins);
                      
                 
                         if (smooth_histograms),
%                              histogram1 = smoothn(histogram1, smoothing_filter_size, filter_type);
%                              histogram2 = smoothn(histogram2, smoothing_filter_size, filter_type);
                               histogram1 = average_smooth(histogram1', histogram_smoothing_filter_size);
                               histogram2 = average_smooth(histogram2', histogram_smoothing_filter_size);
                               histogram1 = histogram1';
                               histogram2 = histogram2';
                               
                         end
                         
                         sum_of_differences = 0;
                            for m = 1:size(histogram1),
                             sum_of_differences =  sum_of_differences + abs(histogram1(m) - histogram2(m));
                         end    
                          difference_image(i,j) = sum_of_differences;
                        end
                        
                   elseif (histogram_window==5),   
                        if (i == 1 || i == 2 || i == size(image2,1) || i == size(image2,1) - 1 ...
                            || j == 1 || j == 2 || j == size(image2,2) || j == size(image2,2) - 1), 
                           difference_image(i,j) = 0;
                        else
                      
                            histogram1 = hist([double(image1(i,j)), double(image1(i,j-1)), double(image1(i,j+1)), double(image1(i-1,j)), double(image1(i+1,j)), ...
                            double(image1(i-1,j-1)),double(image1(i+1,j-1)), double(image1(i-1,j+1)), double(image1(i+1,j+1)), double(image1(i-2,j+2)), ...
                            double(image1(i-2,j+1)), double(image1(i-2,j)), double(image1(i-2,j-1)), double(image1(i-2,j-2)), double(image1(i+2,j+2)), ...
                            double(image1(i+2,j+1)), double(image1(i+2,j)), double(image1(i+2,j-1)), double(image1(i+2,j-2)), double(image1(i+1,j+2)), ...
                            double(image1(i+1,j-2)), double(image1(i,j+2)), double(image1(i,j-2)), double(image1(i-1,j+2)), double(image1(i-1,j-2)), 255, 0], histogram_bins);
                            
                            histogram2 = hist([double(image2(i,j)), double(image2(i,j-1)), double(image2(i,j+1)), double(image2(i-1,j)), double(image2(i+1,j)), ...
                            double(image2(i-1,j-1)),double(image2(i+1,j-1)), double(image2(i-1,j+1)), double(image2(i+1,j+1)), double(image2(i-2,j+2)), ...
                            double(image2(i-2,j+1)), double(image2(i-2,j)), double(image2(i-2,j-1)), double(image2(i-2,j-2)), double(image2(i+2,j+2)), ...
                            double(image2(i+2,j+1)), double(image2(i+2,j)), double(image2(i+2,j-1)), double(image2(i+2,j-2)), double(image2(i+1,j+2)), ...
                            double(image2(i+1,j-2)), double(image2(i,j+2)), double(image2(i,j-2)), double(image2(i-1,j+2)), double(image2(i-1,j-2)), 255, 0], histogram_bins);
                          
                            if (smooth_histograms),
%                              histogram1 = smoothn(histogram1, smoothing_filter_size, filter_type);
%                              histogram2 = smoothn(histogram2, smoothing_filter_size, filter_type);
                               histogram1 = average_smooth(histogram1', histogram_smoothing_filter_size);
                               histogram2 = average_smooth(histogram2', histogram_smoothing_filter_size);
                               histogram1 = histogram1';
                               histogram2 = histogram2';
                            end
                            
                            sum_of_differences = 0;
                            for m = 1:size(histogram1,2),
                             sum_of_differences =  sum_of_differences + abs(histogram1(m) - histogram2(m));
                            end    
                            difference_image(i,j) = sum_of_differences;
                         end 
                   end
                   
                   
                   

                   
                   if (squared_differences),
                       difference_image(i,j) = difference_image(i,j) * difference_image(i,j);
                   end
                   
                   
                   
                else  % caculate proper intensity difference (i.e. Euclidean)
                  if (squared_differences),
                    difference = double(image2(i,j)) - double(image1(i,j));
                    difference_image(i,j) = uint8(difference * difference);
                  else    
                    difference_image(i,j) = uint8(abs(double(image2(i,j)) - double(image1(i,j))));
                  end
                end
                
                if (smoothing_filter_size>1 & smooth_difference_image),
                  difference_image = smoothn(difference_image, smoothing_filter_size, filter_type);
                end
                  % Note: how does smoothing the actual 2 images affect
                  % measures? Is the smoothn function buggy? (very slow)
                  % this filter must be applied to 2 images /before/ the
                  % difference is calculated.
            end
          end
          if (distance_of_power_2_in_difference_image), % raise distances to the power of 2, divide by 255 to avoid overflow
            DistanceMatrix(k,l) = mean(mean(double(difference_image) / 255 * double(difference_image) / 255));
          else
            DistanceMatrix(k,l) = mean(mean(difference_image));
          end
      end
      if (show_images && calculate_fast_shuffle_distance_circle == 0),
        if (scale_brightness),
          highest_intensity = ceil(double(max(max(difference_image))));
        else
          highest_intensity = 255;
        end
        imshow(difference_image,[0 highest_intensity]);
      end
   end

   progress_bar_position = progress_bar_position + 1 / maximum;
   clc;
   disp(['|=================================================|']);
   progress_string='|';       
   for counter = 1:floor(progress_bar_position * 100 / 2),
      progress_string = [progress_string, '#'];
   end
   disp(progress_string);
   disp(['|================= ',num2str(floor(progress_bar_position * 100)),'% completed =================|']);
      % display progress per cent
   steps_remaining = maximum - k;
   minutes = floor(time_for_this_iteration * steps_remaining / 60);
   seconds = rem(floor(time_for_this_iteration *  steps_remaining), 60);
    disp(' ');
    if (seconds > 9),
        disp(['            Estimated remaining time: ', num2str(minutes), ':', num2str(seconds)]);
        disp(' ');disp(' ');
        disp(['Handling ', synthesis_set_path]);
            % show time indicators
    else
        disp(['            Estimated remaining time: ', num2str(minutes), ':0', num2str(seconds)]);
        disp(' ');disp(' ');
        disp(['Handling ', synthesis_set_path, '...']);
    end
    time_for_this_iteration = toc;
end

DistanceMatrix;
  % shows the distances (in model space) matrix (if semi-colon is omitted)

% echo some values of interest
  
Specificity = mean(mean(DistanceMatrix'))
  % find out the mean of distances from model sythesis for each instance of
  % the training set. Return he mean of that.

Specificity_lower_bound = min(mean(DistanceMatrix'))
Specificity_upper_bound = max(mean(DistanceMatrix'))

Generalisability = mean(min(DistanceMatrix))
  % For each model synthesis, find the nearest example in the trainign set.
  % Return the average of these.
  
Generalisability_lower_bound = min(min(DistanceMatrix))
Generalisability_upper_bound = max(min(DistanceMatrix))
  
  
if (scale_brightness),
  highest_intensity = ceil(double(max(max(DistanceMatrix))));
else
  highest_intensity = 255;
end

if (show_distance_matrix_at_end)
  imshow(DistanceMatrix,[0 highest_intensity]);
end
  % show the distance matrix as an image, using the current figure handle

cd(synthesis_set_path);

save(save_filename);
  % Overwrites existing files of matrices for the same data with different
  % algorithm

cd ~/Model_Evaluation;
  % Return to main directory

% End
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%