function [output] = pcafaces(directory, load_FRGC, remove_bg_outliers, round_up_values, show_data_as_surfaces, show_data_as_mesh, show_data_as_images, apply_smoothing, show_cloud_as_raw_mesh, show_histogram, skip_size, frame_size, cropping_method, nose_range, show_data_again_after_normalisation, surface_lighting, show_xy_images, cutoff_x_top, cutoff_x_bottom, cutoff_y_right, cutoff_y_left, plot_xy_diffs, use_icp2, select_subregion, perform_pca, open_expressions, handle1, handle2, handle3, handle4, show_translation_mask, nose_finding_method, pseudo_sphere_radius, matching_iterator, measures_factor, show_shape_disrepency, residue_filtering_method, apply_hole_remover, perform_just_matching, matching_mode, circle_ratio, expressions_fileName, neutrals_fileName)
% PCAFACES: opens a set of face images and performs a PCA analysis with
% Eigen decomposition
%
%    Licensing
%
%      GNU GPL version 3
%
%    GENERAL
%
%      Script
%
%    INPUT/S
%
%      See function PCAFACES_WRAPPER for inputs documentation
%
%    OUTPUT/S
%
%      -output is a placeholder (may be of use in the future)
%
%    PENDING WORK
%
%      -
%
%    KNOWN BUG/S
%
%      -
%
%    COMMENT/S
%
%      -
%
%    RELATED FUNCTION/S
%
%       PCAFACES_WRAPPER
%
%    ABOUT
%
%      -Created:     February 2011
%      -Last update: April 2011
%      -Revision:    0.0.1
%      -Author:      Dr. R. S. Schestowitz
% ==============================================================

program_path='~'; % extensible in case it moves around the home directory

% ==============================================================
% Data handling
% ==============================================================

rand('seed',sum(100*clock))
% rng('default')
% starting_image_number=0+ceil(rand*50)*2-1; % first image, N+1 (choose odd numbers); this determines whether to skip the first N images

starting_image_number=15*2-1;
% starting_image_number=601;
image_number_offset_size=1; % trick for making more pseudo-pairs from a series of predefined pairs

take_evens=0; % limit to avoid loading PPMs in some datasets

% ==============================================================
% Texas DB
% ==============================================================

texas_mode=1; % whether Texas-style range images should be assumed for the pipeline (hack to make it inclusive of
	      % another database, could use redsign)

% for texas_mode:

  y_half_texas_DB=376; % constants, make code more elegant though
  x_half_texas_DB=251;

  y_half_texas_DB_left=300;
  x_half_texas_DB_left=201;

  y_half_texas_DB_right=300;
  x_half_texas_DB_right=301;

  %  texas_boundaries_i=305;
  %  texas_boundaries_j=170;

  texas_boundaries_i=100;
  texas_boundaries_j=70;
  texas_boundaries_bottom=310;

% ==============================================================
% Perturbation
% ==============================================================

perturb_points=1; % separate points a little
perturb_points_factor=100000;

% ==============================================================
% ICP-related
% ==============================================================

disable_ICP_translation=1; % use translation based on nose finder only if enabled
GIP_ICP=1;
rotate=0; % apply rotation in accordance with ICP
show_rotation_mask=0; % show image after rotation, if any is applied at all

% ==============================================================
% Display options
% ==============================================================

show_subregion=0; % show the subregion in question (3-D)
show_noses=0; % show just the noses for the sake of comparison or debugging

% ==============================================================
% Model assessment
% ==============================================================

% model_comparator='median_abs';
model_comparator='determinant';
n_eigenvalues=80; % number of eigenvalues to use for 'determinant' similarity measure, make 'n_eigenvalues' 0 to just choose all eigenvalues
delta=0.001;

% ==============================================================
% Smoothing
% ==============================================================

two_d_smooth=1;
smooth_1d_method='moving'; % if 1-D
	    % also available:
	    %              'moving'   - Moving average (default)
	    %              'lowess'   - Lowess (linear fit)
	    %              'loess'    - Loess (quadratic fit)
	    %              'sgolay'   - Savitzky-Golay
	    %              'rlowess'  - Robust Lowess (linear fit)
	    %              'rloess'   - Robust Loess (quadratic fit)
	    % Note: argument syntax may vary for some of the above

smooth_param_1=13; % usually the level/window size of smoothing

% ==============================================================
% Nose-finding
% ==============================================================

nose_crop_crop_from_top=4/10; % ranges for scanning for nose in given images
nose_crop_crop_from_top_fallback=5.5/10; % fallback if conditions not met
nose_crop_crop_from_bottom=8.5/10;
nose_steepness=31;

% ==============================================================
% Difference, MDS, and GMDS
% ==============================================================

% difference_type='plain_pixels';
difference_type='mds'; % use 'mds' (also applicable for GMDS), 'gradient' or others

% MDS/GMDS options

apply_gmds=1; % difference_type should also be set to 'mds'
visualise_gmds=1; % whether or not to show GMDS results (assuming it's run)
gmds_N =400; % number of points in GMDS
gmds_levels=15;

mds_skip = 5; % spacing between points to sample for GMDS

scale_index_gmds=1; % used for changing ratio along dimensions

show_mds_images=0; % also applies to visualisation of GMDS (before points are placed)
show_mds_iterations=0; % whether to animate MDS process after its calculations

triangulation_method='delaunay'; % how to turn the set of points into triangles

show_local_stress=0; % for a stress map to be dipslayed as a image

mds_min=-0.25; % MDS-only option
mds_max=0.2;
mds_leap=0.05; % step size between samples

size_factor_mds=0.01; % rescaling MDS based on the data, ad hoc

apply_geodesic_mask=1; % boolean - whether to cull out items outside a particular region 
% geodesic_mask_distance_nose=10.8264e-05+1; % use zeroes to just not take those points into account
% geodesic_mask_distance_eyes=10.8264e-05+1; % area around left and right eyes

geodesic_mask_distance_nose=105; % use zeroes to just not take those points into account
geodesic_mask_distance_eyes=105; 
support_buffer=0;

average_GMDS_sampling=0; % take average in 5x5 window

GMDS_remesh=1; % whether or not to remesh to simplify GMDS

  % remesh_limit=600; % how many vertices to keep for GMDS
  % remesh_limit_second=600; 

remesh_limit=4000; % how many vertices to keep for GMDS
remesh_limit_second=4000; 

% ==============================================================
% FMM-based
% ==============================================================

measure_box_geodesics=1; % whether to measure and echo some dimensions estimates
% penalty_FMM_error=5^5; % the penalty score to use for inability to detect geodesic ring regions
penalty_FMM_error=5; % the penalty score to use for inability to detect geodesic ring regions
error_counter=0; % keeps track of inability to detect geodesic ring regions
rotate_fmm_measure_step=1; % how many degrees to skip at every one step
fmm_include_eyes=1; % whether to cacluate this for both eyes
fmm_include_nose=1; % ...and for nose.
fmm_measure_all_parts=0; % nose and also eyes
show_debug_fmm=0;
nose_weighting=2;

    dilation_start=1;
%   dilation_step=0.0264e-05;
%   max_dilation=2.0264e-05+1;

    dilation_step=5;
    max_dilation=41;

% ==============================================================
% Diffusion
% ==============================================================

n_bins_diffusion=250;
test_diffusion_scores=0;

% ==============================================================
% Custom selection
% ==============================================================

reduce_face_region=1; % set to zero for full face (no cropping beneath eyes and nose)
smooth_residual=0; % smoothing the difference to annul noise
use_second_feature=0; % if enabled, the objective function is dual and PCA is a mixture/hybrid of parameter sources

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gather and echo some information about machine status

system('w');
system('uptime');
% system('top');
system('pwd');
% system('find | grep .mat');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Overriden settings

show_shape_disrepency=0;
% residue_filtering_method='threshold'
matching_mode=0;
select_subregion=0;
residue_filtering_method='none';
apply_hole_remover=0;
% use_icp2;


% nose_finding_method='nearest';
nose_finding_method='fixed_centre';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% start of algorithm

if (perform_just_matching)
  display('Run this function in matching mode')
else % train model first, based on CLI/GUI input
  if (matching_mode)
    display('Opening model data... this can take a while.');
    cd(strcat(program_path,'/MAT'));
    load pcaexperiment
    display('Loaded the following data:');
    matching_mode=1; % override the saved data from the training runtime (could also be set accordinglu before saving in that previous runtime)
    match_score_mean=0;
    match_score_concat=0; % override some counters
    match_score_concat_score=0;
    match_score_concat_latent=0;
    n_images_to_open=943;
    who
    vectorised_residuals_original=vectorised_residuals;
  end

  if (load_FRGC) % for the real experiements, otherwise debugging with reduced set
    [images_list]=load_images('FRGC');
    % size(images_list) % should be 4950 for NIST/FRGC-2.0
  else % deal with pairs


    % [images_list]=load_images('pairs'); % mere samples (as below) are an old, deprecated option



    % [images_list]=load_images('samples'); % see descriptions in loader file

    % [images_list]=load_images('training');

    % [images_list]=load_images('targets_expression_pairs');
    % [images_list]=load_images('unseen_targets_expression_pairs');

    % [images_list]=load_images('targets');



    % [images_list]=load_images('same_person_training_fall');

    % [images_list]=load_images('targets_fall');



    % [images_list]=load_images('same_person_training_fall_and_more');
    
    % [images_list]=load_images('targets_fall_and_more');



    %[images_list]=load_images('Texas3DFRDatabase_training');

% [images_list]=load_images('Texas3DFRDatabase_file_false_1_to_100');
[images_list]=load_images('Texas3DFRDatabase_file_true_95_to_97');
 
% [images_list]=load_images('GIP_December_2011');

    % [images_list]=load_images('Texas3DFRD_synthetic');

    if (take_evens)
      images_list=images_list(2:2:size(images_list,1)); % remove PPMs/metadata, optionally filter for some sets where this omission is necessary
    end

%       % randomising:

%      [images_raw]=load_images('targets');
%      images_raw=images_raw(2:2:size(images_raw,1));
%      images_list=images_raw;
%      for index=1:size(images_raw,1)
%        images_list{index}=images_raw{floor(rand(1)*2000)+1};
%      end
%  
%      directory=strcat(program_path,'/Facial-Expressions-Recognition/Anger.v3r');
%      [images]=load_expressions(directory);

  end


  % images_list
  % size(images_list)

  if (open_expressions) % this overrides the options above for the special case of expressions dataset


    nose_range=nose_range/measures_factor; % some scaling due to granularity differences applied to elevant parametes
    cutoff_x_top=cutoff_x_top/measures_factor;
    cutoff_x_bottom=cutoff_x_bottom/measures_factor;
    cutoff_y_right=cutoff_y_right/measures_factor;
    cutoff_y_left=cutoff_y_left/measures_factor;

    n_images_to_open=load_expressions(expressions_fileName); % get the number of availabke frames
    n_neutral_images=load_expressions(neutrals_fileName); % get the number of neutral images

    if (n_images_to_open>n_neutral_images) % implementation-dependent restriction
      display('Warning: not enough neutral images for pairing with expressions set. Ensure there are enough pairs for training. Will attempt training anyway.');
    end
  else
    n_images_to_open=size(images_list, 1); % just get the size of the data structure
  end
  display(' ');
  display(strcat('Number of images to open: ', num2str(n_images_to_open)));

  for image_n=starting_image_number:1:n_images_to_open,
%  n_images_to_open*2
%  %  images_list{image_n}

    if (open_expressions)
      if (rem(image_n,2))
	if (matching_mode)
	  [header, x, y, z, Textures, FrameProperties] = readGipFile(expressions_fileName, matching_iterator); % open the same example for matching all along
	  offset_size=rem((matching_iterator)*40, size(x,1));   % correct offset of 40 per frame
	else
	  [header, x, y, z, Textures, FrameProperties] = readGipFile(expressions_fileName, (image_n-1)/2); % starting at 0, so offset with minus one
	  offset_size=rem(((image_n-1)/2)*40, size(x,1));
	end
      else  % second image in a pair, handling pairs at a time
	[header, x, y, z, Textures, FrameProperties] = readGipFile(neutrals_fileName, image_n/2-1);
	% use different neutral images; for single reference neutral image just load the same data instance
	 offset_size=rem((image_n/2-1)*40, size(x,1));
      end
      
      xx = x; % start to clip invalid data (reused from sample code of the source of data)
      yy = y;
      zz = z;

  % horizontal shift

  %      for i=1:offset_size,
  %          x(:,i)=xx(:,size(xx,2)-offset_size+i);
  %  	y(:,i)=yy(:,size(xx,2)-offset_size+i);
  %          z(:,i)=zz(:,size(xx,2)-offset_size+i);
  %      end
  %  
  %      for i=offset_size+1:size(xx,2),
  %          x(:,i)=xx(:,i-offset_size);
  %          y(:,i)=yy(:,i-offset_size);
  %          z(:,i)=zz(:,i-offset_size);
  %      end

  % vertical shift - a hack said to be needed due to layout of sequence

      loop_span=size(xx,1);
  %  figure;
  %  imshow(z)
      for i=1:offset_size,
  % i
	  x(i,:)=xx(loop_span-offset_size+i,:);
	  y(i,:)=yy(loop_span-offset_size+i,:);
	  z(i,:)=zz(loop_span-offset_size+i,:);
      end

      for i=offset_size+1:loop_span,
  % i
	  x(i,:)=xx(i-offset_size,:);
	  y(i,:)=yy(i-offset_size,:);
	  z(i,:)=zz(i-offset_size,:);
      end

      carving_threshold=0.3;

      minX = -1.5;
      maxX = 1.5;
      minY = -1.5;
      maxY = 1.5;
      minZ = -2.0;
      maxZ = 1.0;

      x(z<minZ) = 0;
      x(z>maxZ) = 0;
      y(z<minZ) = 0;
      y(z>maxZ) = 0;
      z(z<minZ) = 0;
      z(z>maxZ) = 0;

      TH = 5.0 / 100; %%% mm
      m = [3, 3];
      H = fspecial('average', m);
      XXm = imfilter(x, H);
      YYm = imfilter(y, H);
      ZZm = imfilter(z, H);

      x(abs(x - XXm) > TH) = -1;
      y(abs(y - YYm) > TH) = -1;
      z(abs(z - ZZm) > TH) = -1;

    else

      if (texas_mode)
	if (rem(image_n,2))
	  z=imread(images_list{image_n});
	else
	  if (image_number_offset_size) % non-zero, i.e. offset exists
	    z=imread(images_list{image_n+image_number_offset_size});
	  else
	    z=imread(images_list{image_n});
	  end
	end
	z=z(1:750,1:500);
	x=z;
	y=z;
      else
	images_list{image_n}; % echo for debugging; semi-colon to suppres
      %      [x, y, z, fl]=absload(images_list{image_n}); % get data from file/s
	    %  X, Y, and Z are the 3-D co-ordinates the points, whereas FL specifies if a point is valid

	carving_threshold=30;
      end

%  size(z)
%  figure
%  imshow(z)


% for ccde consistency (inefficient)


      % % for debugging
      % fl
      % mean(x)
      % mean(y)
      % mean(z)
      % size(fl)
      % mean(fl)
      % imshow(fl)
    end

    % % edge analysis

    % figure(500)
    % plot((y(300,1:400))-(y(300,2:401)))
    % figure(600)
    % plot((x(300,1:400))-(x(300,2:401)))

    % [x, y, z, fl]=absload(strcat(program_path,'/NIST/FRGC-2.0-dist/nd1/Spring2004range/'));

    % reorganise the data

    if (round_up_values)
      for i=1:1:size(x,1)
	for j=1:1:size(x,2)
	  combined_data(i,j,:)=[floor(x(i,j)),floor(y(i,j)),floor(z(i,j))];
	  % combined_data(i,j,3)=interp2(combined_data(i,:,3),combined_data(:,j,3),...,'cubic')
	  % combined_data(i,j,:)=[rand,rand,rand];
	end
      end
    else % keep values as they were
      for i=1:1:size(x,1)
	for j=1:1:size(x,2)
	  combined_data(i,j,:)=[x(i,j),y(i,j),z(i,j)];
	end
      end
    end

    % histmax=-2060; % speculative initial estimate

    if (remove_bg_outliers)
      if (open_expressions) % just run a median filter with 11x11 kernel to begin with
	for i=6:1:size(x,1)-6
	  for j=6:1:size(x,2)-6
	      % if (i>5 && i<size(x,1)-5 && j>5 & j<size(x,2)-5)
  % median([[combined_data([i-5:i+5],j,3)]])
  % [combined_data([i-5:i+5],j,3)]
		combined_data(i,j,3)=(median([combined_data([i-5:i+5],j,3)])+median([combined_data(i,[j-5:j+5],3)]))/2;
	      % end
	  end
	end
      else % run the steps applying to the GC dataset
	for i=1:1:size(x,1)
	  for j=1:1:size(x,2)
	    if (combined_data(i,j,1)==-999999) combined_data(i,j,1)=NaN; end
	    if (combined_data(i,j,2)==-999999) combined_data(i,j,2)=NaN; end
	    if (combined_data(i,j,3)==-999999)
		  % detect non-background outliers
		  % size(combined_data(:,:,3))
		  % combined_data(i,j-1:1:j+1,3)=removeoutliers(combined_data(i,j-1:1:j+1,3)); % remove outliers from data using the Thompson Tau method
		  % combined_data(i-1:1:i+1,j,3)=removeoutliers(combined_data(i-1:1:i+1,j,3));
		  if (i>frame_size && j>frame_size && i<size(x,1)-frame_size && j<size(x,2)-frame_size)
		    x_mean=mean(combined_data(i-frame_size:1:i+frame_size,j,3));
		    y_mean=mean(combined_data(i,j-frame_size:1:j+frame_size,3));
		    if (mean([x_mean, y_mean]) > -300000) % old test, not needed if image cropped to face
		    % if (1)
		      combined_data(i,j,3)=median([x_mean, y_mean]);
		      % combined_data(i,j,3)
		    else
		      combined_data(i,j,3)=NaN;
		    end
		  else
		      combined_data(i,j,3)=NaN;
		  end
	    end

	    if (i>frame_size && j>frame_size && i<size(x,1)-frame_size && j<size(x,2)-frame_size) % find face fringes (for the most part)
	      if ((abs(combined_data(i,j,3)-mean(combined_data(i-frame_size:1:i+frame_size,j,3)))>100) && (abs(combined_data(i,j,3)-mean(combined_data(i,j-frame_size:1:j+frame_size,3)))>100))
		% i
		% j
		% combined_data(i-10:1:i+10,j,3)
		without_nas=removeoutliers(combined_data(i-frame_size:1:i+frame_size,j,3));
		data_start=0;
		for k=1:1:size(without_nas,1)
		    if (without_nas(k)>-9999)
		      if (data_start==0)
			data_start=k;
		      end
		    end
		end
		x_mean=mean(without_nas(data_start:1:size(without_nas,1)));

		without_nas=removeoutliers(combined_data(i,j-frame_size:1:j+frame_size,3));
		data_start=0;
		for k=1:1:size(without_nas,2)
		    if (without_nas(k)>-9999)
		      if (data_start==0)
			data_start=k;
		      end
		    end
		end
		% without_nas(data_start:1:size(without_nas,2))
		y_mean=mean(without_nas(data_start:1:size(without_nas,2)));
		mean_without_nas=median([x_mean, y_mean]);
		combined_data(i,j,3)=mean_without_nas;
	      end
	    end
	  end
	end
      end
    end 

    if (apply_smoothing) % smooth data (Gaussian filter, using proprietary MATLAB toolbox)
      if (two_d_smooth)
	for i=1:1:size(x,1)
	  for j=1:1:size(x,2)	
	    if (~(combined_data(i,j,3)>-30000))
	      combined_data(i,j,3)=-10;
	    end
	  end
	end
	combined_data(:,:,3)=smooth2(combined_data(:,:,3),smooth_param_1,smooth_param_1);
	for i=1:1:size(x,1)
	  for j=1:1:size(x,2)	
	    if (combined_data(i,j,3)==-10)
	      combined_data(i,j,3)=NaN;
	    end
	  end
	end
      else
%  figure(60)
%  imshow(combined_data(:,:,3),[-50,50])
	for i=1:1:size(x,1) % first vertical
	  combined_data(i,:,3)=smooth(double(combined_data(i,:,3)), smooth_param_1, smooth_1d_method);
	  % combined_data(i,:,3)=smooth(combined_data(i,:,3),'sgolay' );	
	end
	for j=1:1:size(x,2) % now horizontal
	  combined_data(:,j,3)=smooth(double(combined_data(:,j,3)), smooth_param_1, smooth_1d_method);
	  % combined_data(:,j,3)=smooth(combined_data(:,j,3), 'sgolay');
	end
%  figure(61)
%  imshow(combined_data(:,:,3),[-50,50])
      end
    end

    if (apply_hole_remover)
      combined_data_z=combined_data(:,:,3);
      % size(combined_data_z)
      load filter;

      if (open_expressions==0) 
	Wd=Wd'; % fit dimensions to data
      end

      Wd = Wd>0.3; %99;


      Ws = 1-Wd; % imshow(Wd,[]);
      combined_data_z(isnan(combined_data_z)) = 0;
      Z0 = combined_data_z;

      W = 1-min(bwdist(Ws),3)/3;
      Rval = rand(size(combined_data_z));
      combined_data_z(W>0.3)= Rval(W>0.3);

      iter=20;
  %      combined_data_z=LOD_hole_filling(W,combined_data_z,Z0,55,iter);
  %      combined_data_z=LOD_hole_filling(W,combined_data_z,Z0,5,iter);
  %      combined_data_z=LOD_hole_filling(W,combined_data_z,Z0,1.5,20);


    combined_data_z=hole_filling(Wd,Ws,combined_data_z,Z0);

    combined_data(:,:,3)=combined_data_z;

      % combined_data=LOD_hole_filling(combined_data);
      % size(combined_data_z)
    end


  %    % run a second step, median if needed
  %    
  %    for i=6:1:size(x,1)-6  % run median filter (if not already run, as in the expressions pre-processing step (reused code, edit in tandem if possible)
  %      for j=6:1:size(x,2)-6
  %        combined_data(i,j,3)=(median([combined_data([i-5:i+5],j,3)])+median([combined_data(i,[j-5:j+5],3)]))/2;
  %      end
  %    end

    if strcmp(nose_finding_method, 'ViolaJones') 
	ConvertHaarcasadeXMLOpenCV(strcat(program_path,'/pcafaces/Imported/ViolaJones_version0b/HaarCascades/haarcascade_frontalface_alt.xml'));

	nose_x=100;
	nose_y=100;
    end 

    if strcmp(nose_finding_method, 'icp') % use ICP to template, so open and process templates
      [template_header, template_x, template_y, template_z, template_Textures, template_FrameProperties] = readGipFile(strcat(program_path,'/Facial-Expressions-Recognition/Disguse.v3r', 0));
      
      template_x(template_z<minZ) = 0;
      template_x(template_z>maxZ) = 0;
      template_y(template_z<minZ) = 0;
      template_y(template_z>maxZ) = 0;
      template_z(template_z<minZ) = 0;
      template_z(template_z>maxZ) = 0;
      
      % template=[[template_x]; [template_y]; [template_z]];
      for i=100:1:size(template_x,1)-100
	for j=100:1:size(template_x,2)-100
	  template(i,j,:)=[template_x(i,j),template_y(i,j),template_z(i,j)];
	end
      end

      accumulator=1;
      for i=1:skip_size:size(x,1)
	for j=1:skip_size:size(x,2)
	  if ((combined_data(i,j,1)>-30000) && (combined_data(i,j,2)>-30000) && (combined_data(i,j,3)>-30000))
	    combined_data_x_template(accumulator)=combined_data(i,j,1);
	    combined_data_y_template(accumulator)=combined_data(i,j,2);
	    combined_data_z_template(accumulator)=combined_data(i,j,3);
	    accumulator=accumulator+1;
	  end
	end
      end

      accumulator=1;
      for i=100:skip_size:size(template_x,1)-100
	for j=100:skip_size:size(template_x,2)-100
	  if ((template(i,j,1)>-30000) && (template(i,j,2)>-30000) && (template(i,j,3)>-30000))
	    data_x_template(accumulator)=template(i,j,1);
	    data_y_template(accumulator)=template(i,j,2);
	    data_z_template(accumulator)=template(i,j,3);
	    accumulator=accumulator+1;
	  end
	end
      end

      nose_x=100;
      nose_y=100;
      % size(template)
    end



    if (show_histogram) % visualising for debugging or study purposes
      % figure(1) 
      figure;
      hist(combined_data(:,:,3))
    end

%        figure(1) 
%        imshow(combined_data(:,:,3))

    if (texas_mode)

      histmin=0;
      histmax=0;

    else
      [data,hist_numbers]=hist(combined_data(:,:,3)); % get histogram data
      histmin=hist_numbers(1); % get range
      histmax=hist_numbers(10);

    end

    % now crop the images consistently in the following simple sentence

    % first, find the tip of the nose
    if strcmp(nose_finding_method, 'nearest') % rectangle
      if (texas_mode)
	% size(combined_data)
	mat_temp=combined_data(:,:,3);
	[r c]=find(mat_temp==max(mat_temp(:)))
	      max(max(mat_temp(:)))
	top_candidate_position(1)=floor(mean(r))
	top_candidate_position(2)=floor(mean(c))
      else

	% size(x,1)
	% size(x,1)*4/10
	% (size(x,1)*4/10):size(x,1)
	% figure;
	% imshow(combined_data((size(x,1)*nose_crop_crop_from_top):size(x,1),1:size(x,2)*nose_crop_crop_from_top,3), [-3000,-2000])
	
	
	top_candidate=max(max(combined_data((size(x,1)*nose_crop_crop_from_top):(size(x,1)*nose_crop_crop_from_bottom),1:size(x,2),3))); % get points closest to camera (likely nose, assuming no noise/tilt/etc.)

      % size(x,1)*nose_crop_crop_from_top
      % size(x,1)
      % size(x,2)
      
	top_candidate_position=findit(combined_data((size(x,1)*nose_crop_crop_from_top):(size(x,1)*nose_crop_crop_from_bottom),1:size(x,2),3), 1); % get position of the said point
	top_candidate_position(1)=top_candidate_position(1)+size(x,1)*nose_crop_crop_from_top-1;

	tip_distances_temp=combined_data(top_candidate_position(1),top_candidate_position(2),3)-combined_data(top_candidate_position(1),top_candidate_position(2)-pseudo_sphere_radius,3)+(combined_data(top_candidate_position(1),top_candidate_position(2),3)-combined_data(top_candidate_position(1),top_candidate_position(2)+pseudo_sphere_radius,3));


	if (tip_distances_temp<nose_steepness || ~(tip_distances_temp>-30000))  % maybe forehead of hair
	  top_candidate=max(max(combined_data((size(x,1)*nose_crop_crop_from_top_fallback):(size(x,1)*nose_crop_crop_from_bottom),1:size(x,2),3)));
	  top_candidate_position=findit(combined_data((size(x,1)*nose_crop_crop_from_top_fallback):(size(x,1)*nose_crop_crop_from_bottom),1:size(x,2),3), 1); 
	  top_candidate_position(1)=top_candidate_position(1)+size(x,1)*nose_crop_crop_from_top_fallback-1;
	end
      end


    elseif strcmp(nose_finding_method, 'fixed_centre')
      if (texas_mode)
	top_candidate_position(1)=y_half_texas_DB;
	top_candidate_position(2)=x_half_texas_DB;
      end
    elseif strcmp(nose_finding_method, 'sphere') % distance between points on sphere or ellipsoid
      for i=ceil(size(x,1)/4):1:size(x,1)-ceil(size(x,1)/4)
	for j=ceil(size(x,2)/4):1:size(x,2)-ceil(size(x,2)/4)
	  tip_distances(i,j)=(combined_data(i,j,3)-combined_data(i,j-pseudo_sphere_radius,3))+(combined_data(i,j,3)-combined_data(i,j+pseudo_sphere_radius,3));
	end
	top_candidate_slices(i)=max(tip_distances(i,:));
      end
      top_candidate=max(top_candidate_slices);
      top_candidate_position=findit(top_candidate_slices', 1); % get position of the said point
      findoutput=findit(tip_distances(top_candidate_position(1),:)',1);
      top_candidate_position(2)=findoutput(1)+6;
      % top_candidate_position=[300,300]
    elseif strcmp(nose_finding_method, 'icp') % use ICP to template
      if (use_icp2)
	    % R is the rotation matrix and t the translation
	[Ro, t, corr_out, error_out, data_out] = icp2([combined_data_x_template;combined_data_y_template;combined_data_z_template], [data_x_template;data_y_template;data_z_template],3) % the latter is the reference image
      else
	[Ro, t, corr_out, data_out] = icp_ajmal([combined_data_x_template;combined_data_y_template;combined_data_z_template], [data_x_template;data_y_template;data_z_template],3);
      end

  %  mean(combined_data_x_template)
  %  mean(data_x_template)
  %  mean(combined_data_y_template)
  %  mean(data_y_template)
  %  mean(combined_data_z_template)
  %  mean(data_z_template)

      %  figure
      %  imshow(combined_data(:,:,3), [histmin,histmax])
      %  figure
      %  imshow(combined_data(:,:,2), [histmin,histmax])
      %  figure
      %  imshow(combined_data(:,:,1), [histmin,histmax])
      %        [R, t, corr, error, data2] = icp2(combined_data, template,3) % the latter is the template image
      %      else
      %        [R, t, corr, data2] = icp_ajmal(combined_data, template,3);
      %  combined_data(1,1,1)-combined_data(1,2,1)
      Rotation=Ro;
      translation=t;

      top_candidate_position(1)=nose_x-ceil(t(1));
      top_candidate_position(2)=nose_y-ceil(t(2));
      top_candidate=combined_data(top_candidate_position(1),top_candidate_position(2),3);

    elseif strcmp(nose_finding_method, 'ViolaJones') 
      % imshow((combined_data(1:size(x,1),1:size(x,2),3)+3)*50)
      % mean((combined_data(1:size(x,1),1:size(x,2),3)))
      Options.Resize=false; 
      % figure_to_save=imshow(combined_data(1:size(x,1),1:size(x,2),3),[histmin,histmax]) % clim
      % imsave(figure_to_save)
      % imwrite(imshow(combined_data(1:size(x,1),1:size(x,2),3),[histmin,histmax]));
      % ObjectDetection(strcat(program_path,'/3.jpg'),strcat(program_path,'/pcafaces/Imported/ViolaJones_version0b/HaarCascades/haarcascade_frontalface_alt.mat'),Options) 
  %      figure;
  %      imshow((combined_data(1:size(x,1),1:size(x,2),3)+3)*50, [0,255])
      ObjectDetection((combined_data(1:size(x,1),1:size(x,2),3)+0)/1,strcat(program_path,'/pcafaces/Imported/ViolaJones_version0b/HaarCascades/haarcascade_frontalface_alt.mat'),Options)

      top_candidate_position(1)=nose_x
      top_candidate_position(2)=nose_y
      top_candidate=combined_data(top_candidate_position(1),top_candidate_position(2),3)


    elseif strcmp(nose_finding_method, 'nearcentre') 
      % Find a point near the centre of the image and see if it is the nose (ad hoc)


      % scan the region near the centre, applying simply tests to it.
      
  %      x_start=150;
  %      y_start=150;
  %      x_end_gap=250;
  %      y_end_gap=250;
      if (open_expressions) 
	x_start=270;
	y_start=140;
	x_end_gap=250;
	y_end_gap=250;
      else
	x_start=100;
	y_start=100;
	x_end_gap=100;
	y_end_gap=100;
      end

      for i=x_start:skip_size:size(x,1)-x_end_gap
	for j=y_start:skip_size:size(x,2)-y_end_gap % this is *very* ad hoc
	  tip_distances(i,j)=(combined_data(i,j,3)-combined_data(i,j-pseudo_sphere_radius,3))+(combined_data(i,j,3)-combined_data(i,j+pseudo_sphere_radius,3));
	end
      end

      top_candidate=max(max(tip_distances));
      top_candidate_position=findit(tip_distances, 1); % get position of the said point

    else
      display('Invalid/wrong method for nose extraction/identification');
    end

    % size(combined_data)
    % combined_data(2,34)
    % combined_data(222,334)
    % combined_data(2,4)

    unique(combined_data);

    if (show_data_as_mesh)
      % size(combined_data)
      % mesh(combined_data)
  %      figure(11)

      axes(handle3);
      axis off;
      mesh(combined_data(1:skip_size:size(x,1),1:skip_size:size(x,2),3))
      colorbar
    end

    if (show_data_as_surfaces)
      %    figure(21)
      %    surf(combined_data(:,:,1))
      %  
      %    figure(22)
      %    surf(combined_data(:,:,2))
      if (surface_lighting)
	figure(23)
	axis off;
	surf(combined_data(1:skip_size:size(x,1),1:skip_size:size(x,2),3),'FaceColor', 'red', 'EdgeColor', 'none')
	camlight left;
	lighting phong;
      else
	% figure(23)
	axes(handle3);
	axis off;
	surf(combined_data(1:skip_size:size(x,1),1:skip_size:size(x,2),3))
	colorbar
      end
    end

    if (show_data_as_images)
      %    figure(31)
      %    imshow(-combined_data(:,:,1), [0,255])
      %  
      %    figure(32)
      %    imshow(-combined_data(:,:,2), [0,255])

      % figure(33)

  %      
  %      axes(toFocus.handles.axes16);
  %      axis off;
  %      'Parent', handles.axes16
  %    toClose = findobj(get(0, 'Children'), 'flat', 'Tag', 'myUniqueIdentifyer_');
  %    delete(toClose)
      axes(handle1);
      axis off;
      imshow(combined_data(:,:,3), [histmin,histmax])
    end

    % also use contour() for other debugging purposes

    if (show_cloud_as_raw_mesh) % experimental
      % pointCloud2rawMesh(combined_data, 0.6, [0 0 1])
      % pointCloud2mesh([x,y,z], 0.6, [0 0 1])
      % pointCloud2mesh(combined_data, 0.6, [0 0 1])
    end

    if (texas_mode)
      average_spacing=0.3;
      % now put the nose at the front where z==0, along with the face
      top_candidate=255-combined_data(y_half_texas_DB,x_half_texas_DB,3);
% combined_data(y_half_texas_DB,x_half_texas_DB,3)
      combined_data(:,:,3)=combined_data(:,:,3)+top_candidate;
    else
      % put the nose at the front where z==0, along with the face
      %  top_candidate
      combined_data(:,:,3)=combined_data(:,:,3)-top_candidate;
      % combined_data(:,:,2)=combined_data(:,:,2)-top_candidate;
      % combined_data(:,:,1)=combined_data(:,:,1)-top_candidate;

      % center the face based on the tip of the nose
      %  top_candidate_position
      offset(1)=top_candidate_position(1)-size(x,1)/2;
      offset(2)=top_candidate_position(2)-size(x,2)/2; % distances from centre, separated X & Y just for clarity
      % offset(1)
      % offset
      noncentered_data=combined_data;

  %  
  %      if (offset(1)>0)
  %        if (offset(2)>0)
	  for i=abs(offset(1))+1:1:(size(x,1)-abs(offset(1)+1)) % avoid stepping out of image bounds
	    for j=abs(offset(2))+1:1:(size(x,2)-abs(offset(2)+1))
	      combined_data(i,j,3)=noncentered_data(i+offset(1),j+offset(2),3); % start cleanly on a similar image with the points shifted
	      combined_data(i,j,2)=noncentered_data(i+offset(1),j+offset(2),2);
	      combined_data(i,j,1)=noncentered_data(i+offset(1),j+offset(2),1);
	    end 
	  end
  %        else
  %          for i=abs(offset(1))+1:1:(size(x,1)-abs(offset(1)+1)) % avoid stepping out of image bounds
  %  	  for j=abs(offset(2))-1:1:(size(x,2)-abs(offset(2)-1))
  %  i 
  %  j
  %  
  %  	    combined_data(i,j,3)=noncentered_data(i+offset(1),j+offset(2),3); % start cleanly on a similar image with the points shifted
  %  	    combined_data(i,j,2)=noncentered_data(i+offset(1),j+offset(2),2);
  %  	    combined_data(i,j,1)=noncentered_data(i+offset(1),j+offset(2),1);
  %  	  end 
  %  	end
  %        end
  %      else
  %        if (offset(2)>0)
  %  	for i=abs(offset(1))-1:1:(size(x,1)-abs(offset(1)-1)) % avoid stepping out of image bounds
  %  	  for j=abs(offset(2))+1:1:(size(x,2)-abs(offset(2)+1))
  %  	    combined_data(i,j,3)=noncentered_data(i+offset(1),j+offset(2),3); % start cleanly on a similar image with the points shifted
  %  	    combined_data(i,j,2)=noncentered_data(i+offset(1),j+offset(2),2);
  %  	    combined_data(i,j,1)=noncentered_data(i+offset(1),j+offset(2),1);
  %  	  end 
  %  	end
  %        else
  %          for i=abs(offset(1))-1:1:(size(x,1)-abs(offset(1)-1)) % avoid stepping out of image bounds
  %  	  for j=abs(offset(2))-1:1:(size(x,2)-abs(offset(2)-1))
  %  	    combined_data(i,j,3)=noncentered_data(i+offset(1),j+offset(2),3); % start cleanly on a similar image with the points shifted
  %  	    combined_data(i,j,2)=noncentered_data(i+offset(1),j+offset(2),2);
  %  	    combined_data(i,j,1)=noncentered_data(i+offset(1),j+offset(2),1);
  %  	  end 
  %  	end
  %        end
  %      end

      top_candidate_position=top_candidate_position-offset;  % amend assumed position of nose tip, should be the same as image size cut in half, e.g. 240, 320 for FRGC


      % cut remaining objects now
      if (offset(1)>0)
	combined_data(size(x,1):-1:size(x,1)-abs(offset(1)),:,3)=NaN;
	combined_data(size(x,1):-1:size(x,1)-abs(offset(1)),:,2)=NaN;
	combined_data(size(x,1):-1:size(x,1)-abs(offset(1)),:,1)=NaN;
      else
	combined_data(1:offset(1),:,3)=NaN;
	combined_data(1:offset(1),:,2)=NaN;
	combined_data(1:offset(1),:,1)=NaN;
      end

      if (offset(2)>0)
	combined_data(:,size(x,2):-1:size(x,2)-abs(offset(2)),3)=NaN;
	combined_data(:,size(x,2):-1:size(x,2)-abs(offset(2)),2)=NaN;
	combined_data(:,size(x,2):-1:size(x,2)-abs(offset(2)),1)=NaN;
	% 1:offset(2)
      else
	combined_data(:,1:offset(2),3)=NaN;
	combined_data(:,1:offset(2),2)=NaN;
	combined_data(:,1:offset(2),1)=NaN;
      end

      % calculating distances

      if (plot_xy_diffs)
	figure(500);
	plot(combined_data(size(x,1)/2,1:639,1)-(combined_data(size(x,1)/2,2:size(x,2),1)));
	figure(600);
	plot(combined_data(1:479,size(x,2)/2,2)-(combined_data(2:size(x,1),size(x,2)/2,2)));
      end


  %  size(x,1)
  %  combined_data(size(x,1)/2)
  %  (1+size(x,2)/2-5)
  %  (size(x,2)/2+5)
  %  combined_data(size(x,1)/2)
  %  2+size(x,2)/2-5
  %  size(x,2)/2+5
  %  combined_data(size(x,1)/2-5+1:size(x,1)/2+5+1,size(x,2)/2,2)-(combined_data(size(x,1)/2-5:(size(x,1)/2+5),size(x,2)/2,2))

      y_samples=combined_data(size(x,1)/2,(1+size(x,2)/2-5):(size(x,2)/2+5),1)-(combined_data(size(x,1)/2,(2+size(x,2)/2-5):(size(x,2)/2+5)+1,1));
      x_samples=combined_data(size(x,1)/2-5+1:size(x,1)/2+5+1,size(x,2)/2,2)-(combined_data(size(x,1)/2-5:(size(x,1)/2+5),size(x,2)/2,2));
      x_samples=x_samples';  

      for g=1:size(x_samples,2)
	if ~(x_samples(1,g)>-1000)
	  x_samples(1,g) = -0.5;
	end
      end

  % x_samples

      for g=1:size(y_samples,2)
	if ~(y_samples(1,g)>-1000)
	  y_samples(1,g) = -0.5;
	end
      end

  % y_samples

      spacing_along_y=abs(mean(y_samples)); % sample differences near the centre

      spacing_along_x=abs(mean(x_samples));

      % having sampled around X and Y for real physical distances, assuming correct aspect ratio, get average for later calculation 

      average_spacing=mean([spacing_along_x, spacing_along_y]);

      % remove unwanted points (cropping of face)

      %    for theta=1:360
      %        theta=theta*pi/180;
      %     	 x=cos(top_candidate_position(1))+round((10+i)*(cos(teta)));
      %      	 y=sin(top_candidate_position(2))+round((10+i)*(sin(teta)));
      %      	 combined_data(x,y,3)=1;
      %    end
    end

    if strcmp(cropping_method, 'rectangle') % rectangle
      for i=1:1:size(x,1) 
	for j=1:1:size(x,2)
	  if (i>(top_candidate_position(1)+nose_range) || (i<(top_candidate_position(1)-nose_range)))
	    if (j>(top_candidate_position(2)+nose_range) || (j<(top_candidate_position(2)-nose_range)))L
	      combined_data(i,j,3)=NaN;
	    end 
	  end
	end
      end
    elseif strcmp(cropping_method, 'ellipsoid') % distance between points on sphere or ellipsoid
      for i=1:1:size(x,1) 
	for j=1:1:size(x,2)
	  if (distance([i,j],top_candidate_position)>nose_range)
	    combined_data(i,j,3)=NaN;
	    combined_data(i,j,2)=NaN;
	    combined_data(i,j,1)=NaN;
	  end
	end
      end
    else % circle
      for i=1:1:size(x,1) 
	for j=1:1:size(x,2)
	  if (sqrt(abs(i-top_candidate_position(1))^2+(abs(j-top_candidate_position(2))^2)/circle_ratio)>nose_range/average_spacing)
% +i*circle_ratio*average_spacing
	    combined_data(i,j,3)=NaN;
	    combined_data(i,j,2)=NaN;
	    combined_data(i,j,1)=NaN;
	  end
	end
      end
    end

    image_before_masking=combined_data;
    % size(combined_data)
    % size(image_before_masking)

    % select just the rigid areas of the face (nose and forehead), as a sort of masking

    if (open_expressions) % annul division from earlier
  %      bottom_cutoff_y_right=top_candidate_position(2)+floor(cutoff_y_right*2*measures_factor/average_spacing)
  %      bottom_cutoff_y_left=top_candidate_position(2)-floor(cutoff_y_left*2*measures_factor/average_spacing)
	bottom_cutoff_y_right=top_candidate_position(2)+floor(cutoff_y_right/average_spacing);
	bottom_cutoff_y_left=top_candidate_position(2)-floor(cutoff_y_left/average_spacing);
    else
      bottom_cutoff_y_right=top_candidate_position(2)+floor(cutoff_y_right/average_spacing);
      bottom_cutoff_y_left=top_candidate_position(2)-floor(cutoff_y_left/average_spacing);
    end

    bottom_cutoff_x_top=top_candidate_position(1)-floor(cutoff_x_top/average_spacing);
  %  average_spacing
  %  cutoff_x_bottom
  %  top_candidate_position(1)
    bottom_cutoff_x_bottom=top_candidate_position(1)+floor(cutoff_x_bottom/average_spacing);
  %  size(x,2)
  %  size(bottom_cutoff_x_bottom)
  %  size(combined_data)
  %  mesh(combined_data(:,:,3))
  %  bottom_cutoff_x_bottom

   % bottom_cutoff_x_bottom
   % size(x,1)
    if (reduce_face_region) 
      combined_data(bottom_cutoff_x_bottom:size(x,1),:,3)=NaN;
      combined_data(bottom_cutoff_x_bottom:size(x,1),:,2)=NaN;
      combined_data(bottom_cutoff_x_bottom:size(x,1),:,1)=NaN;

    % size(x,1)
    % size(x,2)

      for i=bottom_cutoff_x_top:bottom_cutoff_x_bottom
	for j=1:1:bottom_cutoff_y_left
	  combined_data(i,j,3)=NaN;
	  combined_data(i,j,2)=NaN;
	  combined_data(i,j,1)=NaN;
	end
	for j=bottom_cutoff_y_right:1:size(x,2)
	  combined_data(i,j,3)=NaN;
	  combined_data(i,j,2)=NaN;
	  combined_data(i,j,1)=NaN;
	end
      end

%  % remove above nose
%  	  combined_data(i, ,3)=NaN;
%  	  combined_data(i,1:300,2)=NaN;
%  	  combined_data(i,1:300,1)=NaN;


    end
    % for debugging purposes, consider more image/video feedback

    if (show_data_again_after_normalisation) % whether or not to visualise the data again after some warping

      if (texas_mode)
	histmin=0;
	histmax=255;
        skip_size=10; % speeding ICP up (it is not necessary)
      else
	[data,hist_numbers]=hist(combined_data(:,:,3)); % get histogram data
	histmin=hist_numbers(1); % get range
	histmax=hist_numbers(10);
      end
      unique(combined_data);

      if (show_data_as_mesh)
  %        figure(111)
	axes(handle3);
	axis off;
	mesh(combined_data(1:skip_size:size(x,1),1:skip_size:size(x,2),3))
	colorbar
      end

      if (show_data_as_surfaces)
	if (surface_lighting)
	  figure(123)
	  axis off;
	  surf(combined_data(1:skip_size:size(x,1),1:skip_size:size(x,2),3),'FaceColor', 'red', 'EdgeColor', 'none')
	  camlight left;
	  lighting phong;
	else
	  % figure(123)
	  axes(handle4);
	  axis off;
	  surf(combined_data(1:skip_size:size(x,1),1:skip_size:size(x,2),3))
	  colorbar
	end
      end

      if (show_data_as_images)
	% figure(133)
	axes(handle2);
	imshow(combined_data(:,:,3), [histmin,histmax])
      end
    end

    if (show_xy_images)
      figure(333)
      imshow(combined_data(:,:,1), [histmin,histmax])
      figure(233)
      imshow(combined_data(:,:,2), [histmin,histmax])
    end

    % ensure X and Y are also cropped like 

    for i=1:1:size(x,1) 
      for j=1:1:size(x,2)
	if (combined_data(i,j,3)==NaN)
	  combined_data(i,j,1)=NaN;
	  combined_data(i,j,2)=NaN
	end
      end
    end

    if (~rem(image_n,2)) % second image in a pair
      display(strcat('* Preparing for ICP applied to image #',num2str(image_n)));
  %        size(x,1)
  %        size(x,2)


      % register the images into a common frame and use ICP
      % REFERENCE image to align the data to common frame of reference is set to be neutral from the same subject as the non-neutral

      % remove unavailable/invalid points


%  figure;
%  surf(combined_data(:,:,1),combined_data(:,:,2), combined_data(:,:,3))
%  figure;
%  surf(combined_data_previous(:,:,1),combined_data_previous(:,:,2), combined_data_previous(:,:,3))

%        if (show_noses)
%  	figure
%  	plot3(combined_data_x, combined_data_y,combined_data_z)
%  	
%  combined_data(:,:,1), combined_data(:,:,2),combined_data(:,:,3)
%  
%        end

      accumulator=1;
      for i=1:skip_size:size(x,1)
	for j=1:skip_size:size(x,2)
	  if ((combined_data(i,j,1)>-30000) && (combined_data(i,j,2)>-30000) && (combined_data(i,j,3)>-30000))
%  	    combined_data_x(accumulator)=combined_data(i,j,1);
%  	    combined_data_y(accumulator)=combined_data(i,j,2);
	    combined_data_x(accumulator)=i;
	    combined_data_y(accumulator)=j;
	    combined_data_z(accumulator)=combined_data(i,j,3);
	    accumulator=accumulator+1;
	  end
	end
      end

  %      size(combined_data_x)
  %      size(combined_data_y)
  %      size(combined_data_z)

      display('	- Sorting data types');
      accumulator=1; % vectorise
      for i=1:skip_size:size(x,1)
	for j=1:skip_size:size(x,2)
	  if ((combined_data_previous(i,j,1)>-30000) && (combined_data_previous(i,j,2)>-30000) && (combined_data_previous(i,j,3)>-30000))
%  	    combined_data_x_previous(accumulator)=combined_data_previous(i,j,1);
%  	    combined_data_y_previous(accumulator)=combined_data_previous(i,j,2);
	    combined_data_x_previous(accumulator)=i;
	    combined_data_y_previous(accumulator)=j;
	    combined_data_z_previous(accumulator)=combined_data_previous(i,j,3);
	    accumulator=accumulator+1;
	  end
	end
      end

  %      size(combined_data_x_previous)
  %      size(combined_data_y_previous)
  %      size(combined_data_z_previous)
      display(strcat(['	- Affine registration of: '], images_list{image_n-1}, [' and '], images_list{image_n+image_number_offset_size}));
      % transform the original image according to ICP output (data2 is just the masked part), consider both ICP and ICP2



  % consistency check:

  %  size(combined_data_x)
  %  size(combined_data_y)
  %  size(combined_data_z)
  %  
  %  size(combined_data_x_previous)
  %  size(combined_data_y_previous)
  %  size(combined_data_z_previous)



      if (show_noses)
%  	figure
%  	plot(combined_data_y)
%  	figure
%  	plot(combined_data_x)
%  	figure
%  	plot(combined_data_z)
	figure
	plot3(combined_data_x, combined_data_y,combined_data_z)
	figure
	plot3(combined_data_x_previous,combined_data_y_previous,combined_data_z_previous)
      end


      if (~disable_ICP_translation) | (rotate)
	if (use_icp2)



	  if (GIP_ICP)
	    %options.method = @build_Q_qi_P2P;

	    options.method = @build_Q_qi_P2Plane;
    %        options.method_coefficient=0;
	    options.number_of_iterations=15;
    %        options.method2 = @compute_photometric_Q;
    %        options.method2_coefficient=0;
	    options.need_geometric_info=1;
	    options.verbose=false;





	    %  ICP_FUNC = @ICP_quad;
	    %  options.method = @build_Q_qi_P2P;
	    %  options.need_geometric_info=0;
	    %  options.verbose=1;
	    %  options.number_of_iterations=50;





	    %  ICP_FUNC = @ICP_quad;
	    %  options.method = @build_Q_qi_P2Plane;
	    %  options.need_geometric_info=0;
	    %  options.verbose=1;
	    %  options.number_of_iterations=50;


	    data1.X(:,:)=combined_data(:,:,1)';
	    data2.X(:,:)=combined_data_previous(:,:,1)';

	    data1.Y(:,:)=combined_data(:,:,2)';
	    data2.Y(:,:)=combined_data_previous(:,:,2)';

	    data1.Z(:,:)=combined_data(:,:,3)';
	    data2.Z(:,:)=combined_data_previous(:,:,3)';

    %        size(data1.X)
    %        size(data2.X)
    %  
    %        size(data1.Y)
    %        size(data2.Y)
    %  
    %        size(data1.Z)
    %        size(data2.Z)



    %  myMesh=load('myMesh.mat')
    %  myMesh=load('face1.mat')
    %  myMesh=load('test_surface_data.mat')
    %  myMesh=myMesh.myMesh



    %  figure; imshow(combined_data_previous(:,:,1),[-40,40]);
    %  figure; imshow(combined_data_previous(:,:,2),[-40,40]);
    %  figure; imshow(combined_data_previous(:,:,3),[-40,40]);


	    for i=1:skip_size:size(x,1)
	      for j=1:skip_size:size(x,2)
		if (~(combined_data_previous(i,j,1)>-30000))
		  combined_data_previous(i,j,1)=0;
		end
		if (~(combined_data_previous(i,j,2)>-30000))
		  combined_data_previous(i,j,2)=0;
		end
		if (~(combined_data_previous(i,j,3)>-30000))
		  combined_data_previous(i,j,3)=0;
		end
	      end
	    end


	    %    size(data1.combined_data(i,j,1))
	    %    size(data2.combined_data_previous(i,j,1))
	    %  
	    %    size(data1.combined_data(i,j,2))
	    %    size(data2.combined_data_previous(i,j,2))
	    %  
	    %    size(data1.combined_data(i,j,3))
	    %    size(data2.combined_data_previous(i,j,3))

	    % load myMesh;

	    myMesh1.X=combined_data_x';
	    myMesh2.X=combined_data_x_previous';

	    myMesh1.Y=combined_data_y';
	    myMesh2.Y=combined_data_y_previous';

	    myMesh1.Z=combined_data_z';
	    myMesh2.Z=combined_data_z_previous';

    %  	runme_photometric1;
    %  	runme_photometric2;
    %  	runme_photometric3;



    %	[R, t, data_out, stats] =  ICP_quad(data1,data2,options)

    % 	[R, t, data_out, stats] =  ICP_quad(myMesh1,myMesh2,options)
	    [R t myMesh3, stats]=ICP_geo_quad(myMesh1,myMesh2);
	  else

	    % R is the rotation matrix and t the translation
	    [R, t, corr_out, error_out, data_out] = icp2([combined_data_x;combined_data_y;combined_data_z], [combined_data_x_previous;combined_data_y_previous;combined_data_z_previous],3); % the latter is the reference image

	    % size(combined_data_x)


	  end
  %  figure
  %  plot3(myMesh1.X, myMesh1.Y, myMesh1.Z,'.r')
  %  figure
  %  plot3(myMesh2.X, myMesh2.Y, myMesh2.Z,'.r')
  %  figure
  %  plot3(myMesh3.X, myMesh3.Y, myMesh3.Z,'.r')

	else

	  if (GIP_ICP)
	    load myMesh;

	    points_max=20000;
	    points_max_previous=20000;
	    points_max=max(size(combined_data_x_previous)); % if exceeded by max
	    points_max_previous=max(size(combined_data_x_previous));


  combined_data_x
  max(combined_data_x)
  % points_max=max(combined_data_x);

	    myMesh1.X=combined_data_x(1:points_max)';
	    myMesh2.X=combined_data_x_previous(1:points_max_previous)';

	    myMesh1.Y=combined_data_y(1:points_max)';
	    myMesh2.Y=combined_data_y_previous(1:points_max_previous)';

	    myMesh1.Z=combined_data_z(1:points_max)';
	    myMesh2.Z=combined_data_z_previous(1:points_max_previous)';

  %    myMesh1
  %    myMesh2

	    % [R, t] = icp_2008(myMesh1,myMesh2,'quat');
	    [R, t] = icp(myMesh1,myMesh2,'pca');
	  else

	    % scale current data example by dividing
	    %	[R, t, corr_out, data_out] = icp_ajmal([combined_data_x/2;combined_data_y/2;combined_data_z/2], [combined_data_x_previous;combined_data_y_previous;combined_data_z_previous],3);

	    [R, t, corr_out, data_out] = icp_ajmal([combined_data_x;combined_data_y;combined_data_z], [combined_data_x_previous;combined_data_y_previous;combined_data_z_previous],3);


	  end
	end
      else
	R=0;
	t=0;
      end

      % also use GIP's implementation of ICP


      % show the suggested transformations

      Rotation=R
      translation=t

      if (disable_ICP_translation)
	t=[0;0;0]
	translation=[0;0;0]
      end

      % image_before_masking_previous = [image_before_masking_previous(1,:,3)+; image_before_masking_previous(2,:,3)+t(2); image_before_masking_previous(3,:,3)+t(3)];    

  %  	  t(1)
  %  	  floor(t(1))

      if (show_translation_mask)
	  before_translation=image_before_masking_previous(:,:,3);
	  % % debugging for translation mask
	  % figure;
	  % imshow(before_translation, [histmin,histmax])
      end

      unaligned_data=image_before_masking_previous;
  %  	  figure
  %  	  imshow(image_before_masking_previous(:,:,3), [histmin,histmax])

      for i=abs(floor(t(1)))+1:1:(size(x,1)-abs(floor(t(1))+1)) % avoid stepping out if image bounds
	for j=abs(floor(t(2)))+1:1:(size(x,2)-abs(floor(t(2))+1))
	    image_before_masking_previous(i,j,3)=unaligned_data(i+floor(t(1)),j+floor(t(2)),3); % start cleanly on a similar image with the points shifted

	    % also shift X and Y
	    image_before_masking_previous(i,j,1)=unaligned_data(i+floor(t(1)),j+floor(t(2)),1); % start cleanly on a similar image with the points shifted

	    image_before_masking_previous(i,j,2)=unaligned_data(i+floor(t(1)),j+floor(t(2)),2); % start cleanly on a similar image with the points shifted

	end 
      end

      % cut remaining objects now
      if (t(1)>0)
	image_before_masking_previous(size(x,1):-1:size(x,1)-abs(floor(t(1))),:,3)=NaN;
	image_before_masking_previous(size(x,1):-1:size(x,1)-abs(floor(t(1))),:,2)=NaN;
	image_before_masking_previous(size(x,1):-1:size(x,1)-abs(floor(t(1))),:,1)=NaN;
      else
	image_before_masking_previous(1:floor(t(1)),:,3)=NaN;
	image_before_masking_previous(1:floor(t(1)),:,2)=NaN;
	image_before_masking_previous(1:floor(t(1)),:,1)=NaN;
      end

      if (t(2)>0)
	image_before_masking_previous(:,size(x,1):-1:size(x,2)-abs(floor(t(2))),3)=NaN;
	image_before_masking_previous(:,size(x,1):-1:size(x,2)-abs(floor(t(2))),2)=NaN;
	image_before_masking_previous(:,size(x,1):-1:size(x,2)-abs(floor(t(2))),1)=NaN;
      else
	image_before_masking_previous(:,1:floor(t(2)),3)=NaN;
	image_before_masking_previous(:,1:floor(t(2)),2)=NaN;
	image_before_masking_previous(:,1:floor(t(2)),1)=NaN;
      end

      image_before_masking_previous(:,:,3)=image_before_masking_previous(:,:,3)+floor(t(3));

      if (show_translation_mask) % debugging for translation mask
	figure;
	imshow(image_before_masking_previous(:,:,3)-before_translation, [1,10])
	% sum(sum(image_before_masking_previous(:,:,3)-before_translation))
      end
  %  	  figure;
  %  	  imshow(image_before_masking_previous(:,:,3)-unaligned_data(:,:,3), [histmin,histmax])


      if (show_rotation_mask) % debugging for rotation mask
	histmin;
	histmax;
	figure;
	imshow(image_before_masking_previous(:,:,1), [-100,100]);
	figure;
	imshow(image_before_masking_previous(:,:,2), [-100,100]);
	figure;
	imshow(image_before_masking_previous(:,:,3), [-100,100]);
      end

% size(myMesh3.X)
      if (rotate)
	% image_before_masking_previous=image_before_masking_previous*R;

%  	 image_before_masking_previous(:,:,1)=myMesh3.X;
%  	 image_before_masking_previous(:,:,2)=myMesh3.Y;
%  	 image_before_masking_previous(:,:,3)=myMesh3.Z;

% R 

	for k = 1 : size(x,1)
	  for j = 1 : size(x,2)

% [image_before_masking_previous(k,j,1) image_before_masking_previous(k,j,2) image_before_masking_previous(k,j,3)]'

	      [XYZ] = R * [image_before_masking_previous(k,j,1) image_before_masking_previous(k,j,2) image_before_masking_previous(k,j,3)]';

% scale current data example by dividing
%	      [XYZ] = R * [image_before_masking_previous(k,j,1)/2 image_before_masking_previous(k,j,2)/2 image_before_masking_previous(k,j,3)/2]';

	      image_before_masking_previous(k,j,1)=XYZ(1);
	      image_before_masking_previous(k,j,2)=XYZ(2);
	      image_before_masking_previous(k,j,3)=XYZ(3);
	  end
	end

	% remove resultant holes

	for i=1:size(x,1)
	  for j=1:size(x,2)
	    if (~(image_before_masking_previous(i,j,1)>-30000))
	      image_before_masking_previous(i,j,1)=0;
	    end
	    if (~(image_before_masking_previous(i,j,2)>-30000))
	      image_before_masking_previous(i,j,2)=0;
	    end
	    if (~(image_before_masking_previous(i,j,3)>-30000))
	      image_before_masking_previous(i,j,3)=0;
	    end
	  end
	end

	if (show_rotation_mask) % debugging for rotation mask
	  figure;
	  imshow(image_before_masking_previous(:,:,3), [-100,100]);
	end

      end

      % sph2cart(rotationMatrix2eulerAngles(R))
      % image_before_masking_previous=rotate(image_before_masking_previous, rotationMatrix2eulerAngles(R));
      
      % function maketform can also help here

      % size(image_before_masking_previous)
      % size(image_before_masking)

      if strcmp(difference_type, 'plain_pixels') % surface-to-surface volumetric measure
%  	figure
%  	imshow(image_before_masking_previous(:,:,3))
%  	figure
%  	imshow(image_before_masking(:,:,3))
%  
%    	figure;
%    	imshow((image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)-image_before_masking_previous(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)).^2)

	scores_diff(image_n/2)=min([mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)-image_before_masking_previous(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)).^2)), ...

	mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,151:1:end-149,3)-image_before_masking_previous(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)).^2)), ...

	mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,152:1:end-148,3)-image_before_masking_previous(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)).^2)), ...

	mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,149:1:end-151,3)-image_before_masking_previous(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)).^2)), ...

	mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,148:1:end-152,3)-image_before_masking_previous(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)).^2)), ...

	mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)-image_before_masking_previous(texas_boundaries_i-2:1:end-texas_boundaries_bottom-2,150:1:end-150,3)).^2)), ...

	mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)-image_before_masking_previous(texas_boundaries_i-1:1:end-texas_boundaries_bottom-1,150:1:end-150,3)).^2)), ...

	mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)-image_before_masking_previous(texas_boundaries_i+1:1:end-texas_boundaries_bottom+1,150:1:end-150,3)).^2)), ...

	mean(mean(abs(image_before_masking(texas_boundaries_i:1:end-texas_boundaries_bottom,150:1:end-150,3)-image_before_masking_previous(texas_boundaries_i+2:1:end-texas_boundaries_bottom+2,150:1:end-150,3)).^2))])



	image_difference=image_before_masking_previous(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:)-image_before_masking(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:); % get point-to-point differences, or shape residuals, <assume images may be of different size

% mean(mean(image_difference))
  %  figure
  %  imshow(image_before_masking_previous(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:), [-100,100]);
  %  
  %  figure
  %  imshow(image_before_masking(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:), [-100,100]);



      elseif strcmp(difference_type, 'mds') 

%  	figure(3)
%  	imshow(image_before_masking(:,:,3),[-200,200])
%  	figure(4)
%  	imshow(image_before_masking_previous(:,:,3),[-200,200])

	% TRIV

%  	for ii=1:10:size(image_before_masking(:,:,3),1)
%  	  for jj=1:10:size(image_before_masking(:,:,3),2)
%  	    if (image_before_masking(ii,jj,1)<9)
%  	        combined_data(ii,jj,1)=0.1;
%  	    end
%  	    if (image_before_masking(ii,jj,2)<9)
%  	        combined_data(ii,jj,2)=0.1;
%  	    end
%  	    if (image_before_masking(ii,jj,3)<9)
%  	         combined_data(ii,jj,3)=0.1;
%  	    end
%  	  end
%  	end

	TRIV=0;
	TRIV_previous=0;
	vertices_x=0;
	vertices_y=0;
	vertices_z=0;
	vertices_x_previous=0;
	vertices_y_previous=0;
	vertices_z_previous=0;


	if (texas_mode) % experimental
	  accumulator=1;
	  accumulator_previous=1;
	  for ii=texas_boundaries_i:mds_skip:size(image_before_masking(:,:,3),1)-texas_boundaries_i
	    for jj=texas_boundaries_j:mds_skip:size(image_before_masking(:,:,3),2)-texas_boundaries_j
	      if (combined_data(ii,jj,3)>10)
		vertices_x(accumulator)=ii/scale_index_gmds;
		vertices_y(accumulator)=jj/scale_index_gmds;
		if (average_GMDS_sampling)
		  vertices_z(accumulator)=mean(mean(image_before_masking(ii-2:ii:ii+2,jj-2:jj:jj+2,3)));
		accumulator=accumulator+1;
		else 
		  vertices_z(accumulator)=image_before_masking(ii,jj,3);
		end
		accumulator=accumulator+1;
	      end
	      if (combined_data_previous(ii,jj,3)>10)
		vertices_x_previous(accumulator_previous)=ii/scale_index_gmds;
		vertices_y_previous(accumulator_previous)=jj/scale_index_gmds;
		if (average_GMDS_sampling)
		  vertices_z_previous(accumulator_previous)=mean(mean(image_before_masking_previous(ii-2:ii:ii+2,jj-2:jj:jj+2,3)));
		else 
		  vertices_z_previous(accumulator_previous)=image_before_masking_previous(ii,jj,3);
		end
		accumulator_previous=accumulator_previous+1;
	      end
	    end
	  end

	  if (apply_geodesic_mask)
%   	  size(vertices_x)
%    	  size(vertices_y)
%  	  length(vertices_x)

	    pos1=length(vertices_x)+1;
	    pos2=length(vertices_x)+2;
	    pos3=length(vertices_x)+3;

	    vertices_x(pos1)=y_half_texas_DB/scale_index_gmds; % make nose centre the first point for geodesic cutoff later
	    vertices_y(pos1)=x_half_texas_DB/scale_index_gmds;
	    vertices_z(pos1)=image_before_masking(y_half_texas_DB,x_half_texas_DB,3);
	    vertices_x_previous(pos1)=y_half_texas_DB/scale_index_gmds;
	    vertices_y_previous(pos1)=x_half_texas_DB/scale_index_gmds;
	    vertices_z_previous(pos1)=image_before_masking_previous(y_half_texas_DB,x_half_texas_DB,3);

	    vertices_x(pos2)=y_half_texas_DB_left/scale_index_gmds;
	    vertices_y(pos2)=x_half_texas_DB_left/scale_index_gmds;
	    vertices_z(pos2)=image_before_masking(y_half_texas_DB_left,x_half_texas_DB_left,3);
	    vertices_x_previous(pos2)=y_half_texas_DB_left/scale_index_gmds;
	    vertices_y_previous(pos2)=x_half_texas_DB_left/scale_index_gmds;
	    vertices_z_previous(pos2)=image_before_masking_previous(y_half_texas_DB_left,x_half_texas_DB_left,3);

	    vertices_x(pos3)=y_half_texas_DB_right/scale_index_gmds;
	    vertices_y(pos3)=x_half_texas_DB_right/scale_index_gmds;
	    vertices_z(pos3)=image_before_masking(y_half_texas_DB_right,x_half_texas_DB_right,3);
	    vertices_x_previous(pos3)=y_half_texas_DB_right/scale_index_gmds;
	    vertices_y_previous(pos3)=x_half_texas_DB_right/scale_index_gmds;
	    vertices_z_previous(pos3)=image_before_masking_previous(y_half_texas_DB_right,x_half_texas_DB_right,3);
	  end

	else % MSD for FRGC
	  accumulator=1;
	  accumulator_previous=1;
	  for ii=130:mds_skip:size(image_before_masking(:,:,3),1)-50
	    for jj=130:mds_skip:size(image_before_masking(:,:,3),2)-50
	      if ((image_before_masking(ii,jj,1)>-30000) && (image_before_masking(ii,jj,2)>-30000) && (combined_data(ii,jj,3)>-30000))
		vertices_x(accumulator)=image_before_masking(ii,jj,1);
		vertices_y(accumulator)=image_before_masking(ii,jj,2);
		vertices_z(accumulator)=image_before_masking(ii,jj,3);
		accumulator=accumulator+1;
	      end
	      if ((image_before_masking_previous(ii,jj,1)>-30000) && (image_before_masking_previous(ii,jj,2)>-30000) && (combined_data_previous(ii,jj,3)>-30000))
		vertices_x_previous(accumulator_previous)=image_before_masking_previous(ii,jj,1);
		vertices_y_previous(accumulator_previous)=image_before_masking_previous(ii,jj,2);
		vertices_z_previous(accumulator_previous)=image_before_masking_previous(ii,jj,3);
		accumulator_previous=accumulator_previous+1;
	      end
	    end
	  end
	end

	% vertices_y

	if strcmp(triangulation_method, 'random') 
	  TRIV=ceil(60*rand(300,3)+1);
	  TRIV_previous=ceil(60*rand(300,3)+1);

    %  size(vertices_x)
    %  [1:(size(vertices_x)-10)]'
    %  
	elseif strcmp(triangulation_method, 'adjacent') 
	  TRIV(:,1)=[1:(size(vertices_x,2)-10)]'
	  TRIV(:,2)=[2:(size(vertices_x,2)-9)]'
	  TRIV(:,3)=[11:size(vertices_x,2)]'
	  TRIV_previous(:,1)=[1:(size(vertices_x_previous,2)-10)]'
	  TRIV_previous(:,2)=[2:(size(vertices_x_previous,2)-9)]'
	  TRIV_previous(:,3)=[11:size(vertices_x_previous,2)]'

	elseif strcmp(triangulation_method, 'delaunay') 

%  	  size(vertices_x)
%  	  size(vertices_y)
	  TRIV=delaunay(vertices_x,vertices_y);
	  TRIV_previous=delaunay(vertices_x_previous,vertices_y_previous);
	else % delaunay by default
	  TRIV=delaunay(vertices_x,vertices_y);
	  TRIV_previous=delaunay(vertices_x_previous,vertices_y_previous);
	end

	if (apply_gmds)
	  tosca1.TRIV=TRIV_previous;
	  tosca1.X=vertices_x_previous;
	  tosca1.Y=vertices_y_previous;
	  tosca1.Z=vertices_z_previous;

	  tosca2.TRIV=TRIV;
	  tosca2.X=vertices_x;
	  tosca2.Y=vertices_y;
	  tosca2.Z=vertices_z;

%  tosca1.X(1)
%  tosca1.Y(1)
%  tosca1.Z(1)
%  tosca2.X(1)
%  tosca2.Y(1)
%  tosca2.Z(1)

    % cd(strcat(program_path,'/MAT'));
    % save tosca tosca1


%   figure;
%   trisurf(tosca1.TRIV,tosca1.X,tosca1.Y,tosca1.Z,'FaceColor', 'green', 'EdgeColor', 'none');

	  if (GMDS_remesh)
        %  size(tosca1.X)
	    tosca_x = remesh(tosca1, set_options('vertices',remesh_limit));
	    tosca_y = remesh(tosca2, set_options('vertices',remesh_limit_second));
%  size(tosca_y)
%  tosca_y
        %  size(tosca1.X)
%  figure;
%  trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z,'FaceColor', 'green', 'EdgeColor', 'none');
	    pos1=length(tosca_y.X)+1;
	    pos2=length(tosca_y.X)+2;
	    pos3=length(tosca_y.X)+3;
	  else
	    tosca_x = tosca1;
	    tosca_y = tosca2;
	  end


	  if (apply_geodesic_mask)

  %	    added_index_x=size(tosca_x.X)+1;
	    added_index_x=pos1;
	    added_index_x_left=pos2;
	    added_index_x_right=pos3;

	    tosca_x.X(added_index_x)=vertices_x_previous(pos1); % restore anchor point at the nose tip
	    tosca_x.Y(added_index_x)=vertices_y_previous(pos1);
	    tosca_x.Z(added_index_x)=vertices_z_previous(pos1);

	    tosca_x.X(added_index_x_left)=vertices_x_previous(pos2); % eye regions, left and right
	    tosca_x.Y(added_index_x_left)=vertices_y_previous(pos2);
	    tosca_x.Z(added_index_x_left)=vertices_z_previous(pos2);

	    tosca_x.X(added_index_x_right)=vertices_x_previous(pos3); 
	    tosca_x.Y(added_index_x_right)=vertices_y_previous(pos3);
	    tosca_x.Z(added_index_x_right)=vertices_z_previous(pos3);

%  	      vertices_x_previous(pos1)
%  	      vertices_y_previous(pos1)
%  	      vertices_z_previous(pos1)
%  
%  	      vertices_x_previous(pos2) 
%  	      vertices_y_previous(pos2)
%  	      vertices_z_previous(pos2)
%  
%  	      vertices_x_previous(pos3) 
%  	      vertices_y_previous(pos3)
%  	      vertices_z_previous(pos3)


  %  figure;
  %  trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z,'FaceColor', 'green', 'EdgeColor', 'none');


  %         added_index_y=size(tosca_y.X)+1;
	    added_index_y=pos1;
	    added_index_y_left=pos2;
	    added_index_y_right=pos3;

	    tosca_y.X(added_index_y)=vertices_x(pos1);
	    tosca_y.Y(added_index_y)=vertices_y(pos1);
	    tosca_y.Z(added_index_y)=vertices_z(pos1);

	    tosca_y.X(added_index_y_left)=vertices_x(pos2);
	    tosca_y.Y(added_index_y_left)=vertices_y(pos2);
	    tosca_y.Z(added_index_y_left)=vertices_z(pos2);

	    tosca_y.X(added_index_y_right)=vertices_x(pos3);
	    tosca_y.Y(added_index_y_right)=vertices_y(pos3);
	    tosca_y.Z(added_index_y_right)=vertices_z(pos3);

%         size(tosca_x.X)
%         size(tosca_x.Y)
%         size(tosca_y.X)
%         size(tosca_y.Y) 
warning off
        if (perturb_points) % to avoid points being dropped, including key points
            tosca_x.X=tosca_x.X+randn(length(tosca_x.X),1)/perturb_points_factor;
            tosca_x.Y=tosca_x.Y+randn(length(tosca_x.Y),1)/perturb_points_factor;
        end
	    tosca_x.TRIV=delaunay(tosca_x.X,tosca_x.Y);
        if (perturb_points)
            tosca_y.X=tosca_y.X+randn(length(tosca_y.X),1)/perturb_points_factor;
            tosca_y.Y=tosca_y.Y+randn(length(tosca_y.Y),1)/perturb_points_factor;
        end
	    tosca_y.TRIV=delaunay(tosca_y.X,tosca_y.Y);
warning on
%         size(tosca_x.TRIV)
%         size(tosca_y.TRIV)
%         size(tosca_x.X)
%         size(tosca_x.Y)
%         size(tosca_y.X)
%         size(tosca_y.Y)
% size(tosca_y.Z)
%         for lll=1:size(tosca_y.X,1)
%         if (tosca_y.X(lll,1)==0)
%             a=lll
%         end
%         end


%  	  if (GMDS_remesh)
%  	    tosca_x = remesh(tosca_x, set_options('vertices',remesh_limit/2));
%  	    tosca_y = remesh(tosca_y, set_options('vertices',remesh_limit_second/2));
%  %  figure;
%  %  trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z,'FaceColor', 'green', 'EdgeColor', 'none');
%  	    
%  	  end


%    tosca_x 
%    figure
%    trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z);

  % tosca_x 
  % tosca_y
%    size(tosca_x.X)
%    size(tosca_x.Y)
%    size(tosca_x.Z)
%    size(tosca_x.TRIV)
%    size(tosca_y.X)

%      figure
%      trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z);
%  
%      figure
%      trisurf(tosca_y.TRIV,tosca_y.X,tosca_y.Y,tosca_y.Z);

% a=0
	    tosca_x_unchanged=tosca_x;

	    [tosca_x,Dx]= apply_mask(tosca_x, [added_index_x,added_index_x_left,added_index_x_right], [geodesic_mask_distance_nose,geodesic_mask_distance_eyes,geodesic_mask_distance_eyes]);
	    tosca_x_post=tosca_x;



	    tosca_y_unchanged=tosca_y;
% a=2
	    [tosca_y,Dy] = apply_mask(tosca_y, [added_index_y,added_index_y_left,added_index_y_right], [geodesic_mask_distance_nose+support_buffer,geodesic_mask_distance_eyes+support_buffer,geodesic_mask_distance_eyes+support_buffer]);
	    tosca_y_post=tosca_y;




%  size(tosca_y)
%  tosca_y
%  size(Dy)


% diffusion_distances_pre(image_n/2)=sum(abs(Dx(end-2,end)-Dy(end-2,end)).^2)+sum(abs(Dx(end-2,end-1)-Dy(end-2,end-1)).^2)+sum(abs(Dx(end-1,end)-Dy(end-1,end)).^2)






% a=6
   %trisurf(TRIV, X, Y, Z, D(:,end)); axis image;
   % drawnow;

%      figure
%      trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z);
%  
%      figure
%      trisurf(tosca_y.TRIV,tosca_y.X,tosca_y.Y,tosca_y.Z);

%  	    tosca_x.TRIV=delaunay(tosca_x.X,tosca_x.Y);
%  	    tosca_y.TRIV=delaunay(tosca_y.X,tosca_y.Y);
triangle_difference(image_n/2)=0;
diffusion_on_geodesic(image_n/2)=0;


gmds_rings_min(image_n/2)=10000000;
gmds_rings(image_n/2)=0;


diffusion_on_geodesic_edges_max(image_n/2)=0;
l_int=0;
	if (measure_box_geodesics)
          geodesic_range_distances(image_n/2)=0
          if (fmm_measure_all_parts)
	      geodesic_range_distances(image_n/2)= ...
              abs((max(tosca_x.X)-min(tosca_x.X))-(max(tosca_y.X)-min(tosca_y.X))) ...
              +abs((max(tosca_x.Y)-min(tosca_x.Y))-(max(tosca_y.Y)-min(tosca_y.Y))) ...
              +abs((max(tosca_x.Z)-min(tosca_x.Z))-(max(tosca_y.Z)-min(tosca_y.Z)));
          end

%  dilation_start
%  dilation_step
%  max_dilation

	      for l=dilation_start:dilation_step:max_dilation,
		l_int=l_int+1; % for non-interger cases, e.g. spectral


%  		tosca_x = apply_mask(tosca_x_unchanged, [added_index_x,added_index_x_left,added_index_x_right], [geodesic_mask_distance_nose-l,geodesic_mask_distance_eyes-l,geodesic_mask_distance_eyes-l]);
%  		tosca_y = apply_mask(tosca_y_unchanged, [added_index_y,added_index_y_left,added_index_y_right], [geodesic_mask_distance_nose+support_buffer-l,geodesic_mask_distance_eyes+support_buffer-l,geodesic_mask_distance_eyes+support_buffer-l]);
%  		
%  		geodesic_range_distances(image_n/2)=geodesic_range_distances(image_n/2)+ ...
%  		abs((max(tosca_x.X)-min(tosca_x.X))-(max(tosca_y.X)-min(tosca_y.X))) ...
%  		+abs((max(tosca_x.Y)-min(tosca_x.Y))-(max(tosca_y.Y)-min(tosca_y.Y))) ...
%  		+abs((max(tosca_x.Z)-min(tosca_x.Z))-(max(tosca_y.Z)-min(tosca_y.Z)))

	        % l % echo number






	    [tosca_x,Dx]= apply_mask(tosca_x_unchanged, [added_index_x,added_index_x_left,added_index_x_right], [geodesic_mask_distance_nose-l,geodesic_mask_distance_eyes-l,geodesic_mask_distance_eyes-l]);
	    tosca_x_post=tosca_x;



	    [tosca_y,Dy] = apply_mask(tosca_y_unchanged, [added_index_y,added_index_y_left,added_index_y_right], [geodesic_mask_distance_nose+support_buffer-l,geodesic_mask_distance_eyes+support_buffer-l,geodesic_mask_distance_eyes+support_buffer-l]);
	    tosca_y_post=tosca_y;




  	  [tosca_x] = init_surface (tosca_x);
	  [tosca_y] = init_surface (tosca_y);

	  [tx, ux, ty, uy, f, rmsdist, maxdist, local_stress] = gmds (tosca_x, tosca_y, gmds_N, gmds_levels);
  
	  if (visualise_gmds)
		  figure(200+image_n);
		  subplot(1,2,1), visualize_shape(tosca_x, tx, ux);
		  subplot(1,2,2), visualize_shape(tosca_y, ty, uy);
	  end   

	  if (show_local_stress)
		    figure(3); imshow(local_stress);
	  end

	  % local_stress
	  local_stress_sum(image_n/2)=sum(sum(local_stress))
	  % local_stress_sum_sqr(image_n/2)=sum(sum(local_stress.^2))

	  % sum_ux(image_n/2)=sum(sum(ux))
	  % sum_uy(image_n/2)=sum(sum(uy))
	  f_acc(image_n/2)=f
	  % rmsdist_acc(image_n/2)=rmsdist
	  % maxdist_acc(image_n/2)=maxdist
	  % size(local_stress)
	  image_difference(:,:,3)=local_stress;

	  gmds_rings_min(image_n/2)=min(gmds_rings_min(image_n/2),f_acc(image_n/2))
	  gmds_rings(image_n/2)=gmds_rings(image_n/2)+f_acc(image_n/2)









		if (fmm_include_nose)
		  [tosca_x,Dx] = apply_mask(tosca_x_unchanged, added_index_x, geodesic_mask_distance_nose-l);
		  [tosca_y,Dy] = apply_mask(tosca_y_unchanged, added_index_y, geodesic_mask_distance_nose+support_buffer-l);

%  length(tosca_x.X)

%   tosca_x 
%      figure
%      trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z,'FaceColor', 'green', 'EdgeColor', 'none');
%  tosca_y
%  length(tosca_y.X)
		  triangle_difference(image_n/2)=abs(length(tosca_x.X)-length(tosca_y.X))+triangle_difference(image_n/2)

% size(Dx)
% size(Dy)
% mean(mean(Dx-Dy))
% mean(mean(Dx(:,-Dy))


% diffusion_distances(image_n/2)=mean(abs(Dx(:,end-2)-Dy(:,end-2)).^2)

  %  max(tosca_x.X)-min(tosca_y.X)
  %  max(tosca_y.X)-min(tosca_y.X)
  %  max(tosca_x.Y)-min(tosca_y.Y)
  %  max(tosca_y.Y)-min(tosca_y.Y)
  %  max(tosca_x.Z)-min(tosca_y.Z)
  %  max(tosca_y.Z)-min(tosca_y.Z)


		  if length(tosca_x.X)>1 && length(tosca_y.X)>1
		    for rotation_angle=0:rotate_fmm_measure_step:180
		      rotAngle = deg2rad(rotation_angle);
		      tosca_x_rotated.X     = tosca_x.X*cos(rotAngle) - tosca_x.Y*sin(rotAngle);
		      tosca_x_rotated.Y     = tosca_x.X*sin(rotAngle) + tosca_x.Y*cos(rotAngle);

		      tosca_y_rotated.X     = tosca_y.X*cos(rotAngle) - tosca_y.Y*sin(rotAngle);
		      tosca_y_rotated.Y     = tosca_y.X*sin(rotAngle) + tosca_y.Y*cos(rotAngle);

      nosepoint(rotation_angle+1, l_int, image_n/2)=nose_weighting*(abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2);

%  		      nosepoint(rotation_angle+1, l_int, image_n/2)=nose_weighting*(abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2 ...
%  		 +abs((max(tosca_x.Z)-min(tosca_x.Z))-(max(tosca_y.Z)-min(tosca_y.Z)))^2);


%        nosepoint(rotation_angle+1, l, image_n/2)=abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2 ...
%  		    +abs((max(tosca_x_rotated.Y)-min(tosca_x_rotated.Y))-(max(tosca_y_rotated.Y)-min(tosca_y_rotated.Y)))^2 ...
%  		    +abs((max(tosca_x.Z)-min(tosca_x.Z))-(max(tosca_y.Z)-min(tosca_y.Z)))^2;

		      geodesic_range_distances(image_n/2)=geodesic_range_distances(image_n/2)+ ...
		  nosepoint(rotation_angle+1, l_int, image_n/2);

		    end



		    if (test_diffusion_scores)
		      D_y=diffusion_distance(tosca_y.TRIV,tosca_y.X,tosca_y.Y,tosca_y.Z);
		      D_x=diffusion_distance(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z);

		      %  figure
		      %  trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z,D_x(1)*1000)
		      %  figure;
		      [hist_bin_x]=sum(hist(D_x,n_bins_diffusion)');

		      %  figure
		      %  trisurf(tosca_y.TRIV,tosca_y.X,tosca_y.Y,tosca_y.Z,D_y(1)*1000)
		      %  figure;
		      [hist_bin_y]=sum(hist(D_y,n_bins_diffusion)');

		      sum(abs(hist_bin_x-hist_bin_y));
		      diffusion_on_geodesic(image_n/2)=diffusion_on_geodesic(image_n/2)+sum(abs(hist_bin_x-hist_bin_y))

		      diffusion_on_geodesic_edges_max(image_n/2)=diffusion_on_geodesic(image_n/2)+abs(max(max(D_y))-max(max(D_x)))
		    end



		  else
		    error_counter=error_counter+1
		    geodesic_range_distances(image_n/2)=geodesic_range_distances(image_n/2)+(penalty_FMM_error*(180/rotate_fmm_measure_step))
		    nosepoint(1:181, l_int, image_n/2)=0;
		  end
		end
		if (l==max_dilation && show_debug_fmm)
		  nosepoint(:, :, image_n/2);
		  size(nosepoint(:, :,image_n/2));
		  figure
		  imshow(nosepoint(:, :,image_n/2)'/255)
		  % nosepoint=rand(181,98,1)*255;
		  size(nosepoint);
		  lines=max_dilation-2;
		  degrees=180/rotate_fmm_measure_step;
		  t_map = linspace(0,lines*pi*2,lines*degrees*2);
		  x_map = t_map.*cos(t_map);
		  y_map = t_map.*sin(t_map);
		  % imshow(x_map,y_map=rand(size(x_map))
		  flat=zeros(1000,1000);
		  % for ll=(size(x_map,2)/2):-1:1
		  for ll=(size(x_map,2)/2):-1:1
		    shade=nosepoint(mod(ll,degrees)+1,floor(ll/degrees)+1 ,image_n/2)'/255;
		    flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+500)=shade;

		    flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+500)=shade;
		    flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+500)=shade;
		    flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+501)=shade;
		    flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+501)=shade;
		    flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+501)=shade;
		    flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+499)=shade;
		    flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+499)=shade;
		    flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+499)=shade;
		  end
		  figure
		  imshow(flat)
		end

%  tosca_x_unchanged
%  tosca_y_unchanged

		if (fmm_include_eyes)
		  tosca_x = apply_mask(tosca_x_unchanged, added_index_x_left, geodesic_mask_distance_eyes-l);
		  tosca_y = apply_mask(tosca_y_unchanged, added_index_y_left, geodesic_mask_distance_eyes+support_buffer-l);
  %  
  %  tosca_x 
  %  tosca_y 
  %  max(tosca_x.X)-min(tosca_y.X)
  %  max(tosca_y.X)-min(tosca_y.X)
  %  max(tosca_x.Y)-min(tosca_y.Y)
  %  max(tosca_y.Y)-min(tosca_y.Y)
  %  max(tosca_x.Z)-min(tosca_y.Z)
  %  max(tosca_y.Z)-min(tosca_y.Z)




		  if length(tosca_x.X)>1 && length(tosca_y.X)>1
		    for rotation_angle=0:rotate_fmm_measure_step:180

  % rotate by 30 clockwise around (0,0)
  % http://en.wikipedia.org/wiki/Rotation_matrix
		      rotAngle = deg2rad(rotation_angle);
		      tosca_x_rotated.X     = tosca_x.X*cos(rotAngle) - tosca_x.Y*sin(rotAngle);
		      tosca_x_rotated.Y     = tosca_x.X*sin(rotAngle) + tosca_x.Y*cos(rotAngle);

		      tosca_y_rotated.X     = tosca_y.X*cos(rotAngle) - tosca_y.Y*sin(rotAngle);
		      tosca_y_rotated.Y     = tosca_y.X*sin(rotAngle) + tosca_y.Y*cos(rotAngle);
 

lefteye(rotation_angle+1, l, image_n/2)=abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2;;

%  		      lefteye(rotation_angle+1, l, image_n/2)=abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2 ...
%  		  +abs((max(tosca_x.Z)-min(tosca_x.Z))-(max(tosca_y.Z)-min(tosca_y.Z)))^2;

%     lefteye(rotation_angle+1, l, image_n/2)=abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2 ...
		  % +abs((max(tosca_x_rotated.Y)-min(tosca_x_rotated.Y))-(max(tosca_y_rotated.Y)-min(tosca_y_rotated.Y)))^2 ...
%  		  +abs((max(tosca_x.Z)-min(tosca_x.Z))-(max(tosca_y.Z)-min(tosca_y.Z)))^2;

		      geodesic_range_distances(image_n/2)=geodesic_range_distances(image_n/2)+ ...
		      lefteye(rotation_angle+1, l, image_n/2); 

  %  figure;
  %  trisurf(tosca_y.TRIV,tosca_y_rotated.X,tosca_y_rotated.Y,tosca_y.Z);
		    end
		  else
		  error_counter=error_counter+1
		    geodesic_range_distances(image_n/2)=geodesic_range_distances(image_n/2)+(penalty_FMM_error*(180/rotate_fmm_measure_step))
		    lefteye(1:181, l, image_n/2)=0;
		  end


		  if (l==max_dilation && show_debug_fmm)
		    flat=zeros(1000,1000);
		    % for ll=(size(x_map,2)/2):-1:1
		    for ll=(size(x_map,2)/2):-1:1
		      shade=lefteye(mod(ll,degrees)+1,floor(ll/degrees)+1 ,image_n/2)'/255;
		      flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+500)=shade;

		      flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+500)=shade;
		      flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+500)=shade;
		      flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+501)=shade;
		      flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+501)=shade;
		      flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+501)=shade;
		      flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+499)=shade;
		      flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+499)=shade;
		      flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+499)=shade;
		    end
		    figure
		    imshow(flat)
		  end

		  tosca_x = apply_mask(tosca_x_unchanged, added_index_x_right, geodesic_mask_distance_eyes-l);
		  tosca_y = apply_mask(tosca_y_unchanged, added_index_y_right, geodesic_mask_distance_eyes+support_buffer-l);



		  if length(tosca_x.X)>1 && length(tosca_y.X)>1
		    for rotation_angle=0:rotate_fmm_measure_step:180

		      rotAngle = deg2rad(rotation_angle);
		      tosca_x_rotated.X     = tosca_x.X*cos(rotAngle) - tosca_x.Y*sin(rotAngle);
		      tosca_x_rotated.Y     = tosca_x.X*sin(rotAngle) + tosca_x.Y*cos(rotAngle);

		      tosca_y_rotated.X     = tosca_y.X*cos(rotAngle) - tosca_y.Y*sin(rotAngle);
		      tosca_y_rotated.Y     = tosca_y.X*sin(rotAngle) + tosca_y.Y*cos(rotAngle);

		      righteye(rotation_angle+1, l, image_n/2)=abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2;

%  		      righteye(rotation_angle+1, l, image_n/2)=abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2 ...
%  		  +abs((max(tosca_x.Z)-min(tosca_x.Z))-(max(tosca_y.Z)-min(tosca_y.Z)))^2;

%      righteye(rotation_angle+1, l, image_n/2)=abs((max(tosca_x_rotated.X)-min(tosca_x_rotated.X))-(max(tosca_y_rotated.X)-min(tosca_y_rotated.X)))^2 ...
		  % +abs((max(tosca_x_rotated.Y)-min(tosca_x_rotated.Y))-(max(tosca_y_rotated.Y)-min(tosca_y_rotated.Y)))^2 ...
%  		  +abs((max(tosca_x.Z)-min(tosca_x.Z))-(max(tosca_y.Z)-min(tosca_y.Z)))^2;

		      geodesic_range_distances(image_n/2)=geodesic_range_distances(image_n/2)+ ...
		  righteye(rotation_angle+1, l, image_n/2);
		     

		    end
		  else
		    error_counter=error_counter+1
		    geodesic_range_distances(image_n/2)=geodesic_range_distances(image_n/2)+(penalty_FMM_error*(180/rotate_fmm_measure_step))
		    righteye(1:181, l, image_n/2)=0;
		  end


		  if (l==max_dilation && show_debug_fmm)
		    flat=zeros(1000,1000);
		    % for ll=(size(x_map,2)/2):-1:1
		    for ll=(size(x_map,2)/2):-1:1
		      shade=righteye(mod(ll,degrees)+1,floor(ll/degrees)+1 ,image_n/2)'/255;
		      flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+500)=shade;

		      flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+500)=shade;
		      flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+500)=shade;
		      flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+501)=shade;
		      flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+501)=shade;
		      flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+501)=shade;
		      flat(floor(x_map(ll*2))'+501,floor(y_map(ll*2))'+499)=shade;
		      flat(floor(x_map(ll*2))'+500,floor(y_map(ll*2))'+499)=shade;
		      flat(floor(x_map(ll*2))'+499,floor(y_map(ll*2))'+499)=shade;
		    end
		    figure
		    imshow(flat)
		  end

	      end
	      end
	    end

	      error_counter
	      geodesic_range_distances
	    if (show_debug_fmm)
	      figure
	      imshow(images_list{image_n-1})
	      figure
	      imshow(images_list{image_n+image_number_offset_size})
	    end

% a=0
	    tosca_x = tosca_x_post;

% a=2
	    tosca_y = tosca_y_post;
% a=6




% combined_data(y_half_texas_DB,x_half_texas_DB,3)
% combined_data_previous(y_half_texas_DB,x_half_texas_DB,3)



%  max(tosca_x.X)-min(tosca_y.X)-max(tosca_y.X)-min(tosca_y.X)
%  max(tosca_x.X)
%  min(tosca_y.X)
%  max(tosca_y.X)
%  min(tosca_y.X)
	  end

	  if (show_mds_images)
	    figure(1),
    % trisurf(TRIV,100*reshape(combined_data(1:10:end,1:10:end,1),size(combined_data,1)*size(combined_data,2)/100,1),100*reshape(combined_data(1:10:end,1:10:end,2),size(combined_data,1)*size(combined_data,2)/100,1),reshape(combined_data(1:10:end,1:10:end,3),size(combined_data,1)*size(combined_data,2)/100,1),'FaceColor', 'red', 'EdgeColor', 'none'); 


	    trisurf(tosca_x.TRIV,tosca_x.X,tosca_x.Y,tosca_x.Z,'FaceColor', 'green', 'EdgeColor', 'none');
	    % trisurf(TRIV,vertices_x,vertices_y,vertices_z,'FaceColor', 'red', 'EdgeColor', 'none');
	    axis off;

	    camlight left;
	    lighting phong;

	    title(strcat('Face as vertices - image #', num2str(image_n)) );
	    drawnow

	    figure(2),
	    
	    trisurf(tosca_y.TRIV,tosca_y.X,tosca_y.Y,tosca_y.Z,'FaceColor', 'green', 'EdgeColor', 'none');
	    % trisurf(TRIV_previous,vertices_x_previous,vertices_y_previous,vertices_z_previous,'FaceColor', 'red', 'EdgeColor', 'none'); 
	    axis off;

	    camlight left;
	    lighting phong;

	    title(strcat('Face as vertices - image #', num2str(image_n)) );
	    drawnow
	  end

  	  [tosca_x] = init_surface (tosca_x);
	  [tosca_y] = init_surface (tosca_y);

	  [tx, ux, ty, uy, f, rmsdist, maxdist, local_stress] = gmds (tosca_x, tosca_y, gmds_N, gmds_levels);
  
	  if (visualise_gmds)
		  figure(200+image_n);
		  subplot(1,2,1), visualize_shape(tosca_x, tx, ux);
		  subplot(1,2,2), visualize_shape(tosca_y, ty, uy);
	  end   

	  if (show_local_stress)
		    figure(3); imshow(local_stress);
	  end

	  % local_stress
	  local_stress_sum(image_n/2)=sum(sum(local_stress))
	  % local_stress_sum_sqr(image_n/2)=sum(sum(local_stress.^2))

	  % sum_ux(image_n/2)=sum(sum(ux))
	  % sum_uy(image_n/2)=sum(sum(uy))
	  f_acc(image_n/2)=f
	  % rmsdist_acc(image_n/2)=rmsdist
	  % maxdist_acc(image_n/2)=maxdist
	  % size(local_stress)
	  image_difference(:,:,3)=local_stress;
	else % just MDS


  % reshape(combined_data(1:10:end,1:10:end,1),size(combined_data,1)*size(combined_data,2)/100,1)

	  if (show_mds_images)
	    figure(1),
    % trisurf(TRIV,100*reshape(combined_data(1:10:end,1:10:end,1),size(combined_data,1)*size(combined_data,2)/100,1),100*reshape(combined_data(1:10:end,1:10:end,2),size(combined_data,1)*size(combined_data,2)/100,1),reshape(combined_data(1:10:end,1:10:end,3),size(combined_data,1)*size(combined_data,2)/100,1),'FaceColor', 'red', 'EdgeColor', 'none'); 

	    trisurf(TRIV,vertices_x,vertices_y,vertices_z,'FaceColor', 'red', 'EdgeColor', 'none');  axis off;

	    camlight left;
	    lighting phong;

	    title(strcat('Face as vertices ', num2str(image_n)) );
	    drawnow


	    figure(2),
	    trisurf(TRIV_previous,vertices_x_previous,vertices_y_previous,vertices_z_previous,'FaceColor', 'red', 'EdgeColor', 'none');  axis off;


	    camlight left;
	    lighting phong;

	    title(strcat('Face as vertices ', num2str(image_n)) );
	    drawnow
	  end

	  % embed using SMACOF
	  disp('Embedding the intrinsic geometry of the face into R^3 using SMACOF (no acceleration)...')
	  options.X0 = [vertices_x',vertices_y',vertices_z'];
	  options.method = 'smacof';
	  options.xhistory = 'on';

	  % size(options.X0)
	  % size(ones(size(vertices_x,2),size(vertices_x,2)))
	  [X_smacof,hist_smacof] = mds(ones(size(vertices_x,2),size(vertices_x,2))*size_factor_mds,options);

	  % embed using SMACOF again
	  disp('Embedding the intrinsic geometry of the face into R^3 using SMACOF (no acceleration)...')
	  options_previous.X0 = [vertices_x_previous',vertices_y_previous',vertices_z_previous'];
	  options_previous.method = 'smacof';
	  options_previous.xhistory = 'on';

	  % size(options_previous.X0)
	  % size(ones(size(vertices_x_previous,2),size(vertices_x_previous,2)))
	  [X_smacof_previous,hist_smacof_previous] = mds(ones(size(vertices_x_previous,2),size(vertices_x_previous,2))*size_factor_mds,options_previous);

	  if (show_mds_iterations)
	    for k = 1:length(hist_smacof.time),

	      %  figure(2), trisurf(rand(6000,3),hist_smacof.X{k}(:,1),hist_smacof.X{k}(:,2),hist_smacof.X{k}(:,3),'FaceColor', 'red', 'EdgeColor', 'none'); axis off;
	      figure(2), trisurf(TRIV,hist_smacof.X{k}(:,1),hist_smacof.X{k}(:,2),hist_smacof.X{k}(:,3),'FaceColor', 'red', 'EdgeColor', 'none'); axis off;
	      camlight left;
	      lighting phong;
	      title(sprintf('Iteration %d, Stress=%g',k,hist_smacof.s(k)));
	      drawnow 
	      pause(0.25);


	      % plot convergence
	      figure(3), semilogy(cumsum(hist_smacof.time),hist_smacof.s,'k');
	      xlabel('CPU time (sec)'); ylabel('stress');
	      title('Convergence speed');
	    end
	  end

	  %  size(TRIV_previous)
	  %  size(TRIV)
	  %  
	  %  size(vertices_x_previous)
	  %  size(vertices_x)

	  if (show_mds_images)
	    figure(3), trisurf(TRIV,hist_smacof.X{length(hist_smacof.time)}(:,1),hist_smacof.X{length(hist_smacof.time)}(:,2),hist_smacof.X{length(hist_smacof.time)}(:,3),'FaceColor', 'red', 'EdgeColor', 'none'); axis off;

	    camlight left;
	    lighting phong;
	    
	    figure(4), trisurf(TRIV_previous,hist_smacof_previous.X{length(hist_smacof_previous.time)}(:,1),hist_smacof_previous.X{length(hist_smacof_previous.time)}(:,2),hist_smacof_previous.X{length(hist_smacof_previous.time)}(:,3),'FaceColor', 'red', 'EdgeColor', 'none'); axis off;

	    camlight left;
	    lighting phong;
	  end

	  gx0=0;
	  for gx=mds_min:mds_leap:mds_max
	    gx0=gx0+1;
	    gy0=0;
	    for gy=mds_min:mds_leap:mds_max
	    gy0=gy0+1;
	      distances(gx0,gy0)=10000;
	      for h=1:size(vertices_x,2)
		if (distances(gx0,gy0)>sqrt(abs(hist_smacof.X{length(hist_smacof.time)}(h,1)-gx)^2+abs(hist_smacof.X{length(hist_smacof.time)}(h,2)-gy)^2))
		  distances(gx0,gy0)=sqrt(abs(hist_smacof.X{length(hist_smacof.time)}(h,1)-gx)^2+abs(hist_smacof.X{length(hist_smacof.time)}(h,2)-gy)^2);
		  z_values(gx0,gy0)=hist_smacof.X{length(hist_smacof.time)}(h,3);
		end
	      end
	    end
	  end

  %  size(vertices_x,2)
  %  sqrt(abs(hist_smacof.X{length(hist_smacof.time)}(h,1)-gx)^2+abs(hist_smacof.X{length(hist_smacof.time)}(h,2)-gy)^2)
  %  distances(gx0,gy0)
  %  hist_smacof.X{length(hist_smacof.time)}(h,3)
  %  
  %  hist_smacof.X{length(hist_smacof.time)}(:,3)
  %  
  %  	distances
  %  z_values




	  gx0=0;
	  for gx=mds_min:mds_leap:mds_max
	    gx0=gx0+1;
	    gy0=0;
	    for gy=mds_min:mds_leap:mds_max
	    gy0=gy0+1;
	      distances_previous(gx0,gy0)=10000;
	      for h=1:size(vertices_x_previous,2)
		if (distances_previous(gx0,gy0)>sqrt(abs(hist_smacof_previous.X{length(hist_smacof_previous.time)}(h,1)-gx)^2+abs(hist_smacof_previous.X{length(hist_smacof_previous.time)}(h,2)-gy)^2))
		  distances_previous(gx0,gy0)=sqrt(abs(hist_smacof_previous.X{length(hist_smacof_previous.time)}(h,1)-gx)^2+abs(hist_smacof_previous.X{length(hist_smacof_previous.time)}(h,2)-gy)^2);
		  z_values_previous(gx0,gy0)=hist_smacof_previous.X{length(hist_smacof_previous.time)}(h,3);
		end
	      end
	    end
	  end


	  % distances_previous
	  % z_values_previous


  %  	for h=1:length(hist_smacof_previous.time)
  %  
  %  	end
  %  
  %  	hist_smacof_previous.X{length(hist_smacof_previous.time)}(h,1)
  %  	hist_smacof_previous.X{length(hist_smacof_previous.time)}(h,2)
  % size(vertices_x_previous)


	  image_difference(:,:,3)=z_values-z_values_previous;
	end



	% figure(2), surf(combined_data(1:10:end,1:10:end,1),combined_data(1:10:end,1:10:end,2),combined_data(1:10:end,1:10:end,3),'FaceColor', 'red', 'EdgeColor', 'none'); 
	% camlight left;
	% lighting phong;
 
      elseif strcmp(difference_type, 'gradient') 
	image_before_masking_previous_temp=image_before_masking_previous(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:);

	image_before_masking_temp=image_before_masking(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:);

	%  figure
	%  imshow(image_before_masking_previous_temp(:,:,3), [-100,100])
	%  figure
	%  imshow(image_before_masking_temp(:,:,3), [-100,100])



	%  size(image_before_masking_previous_temp(:,:,3))
	%  size(diff(image_before_masking_previous_temp(:,:,3)))

	image_before_masking_previous_temp1(:,:,3)=[diff(image_before_masking_previous_temp(:,:,3));zeros(1,size(image_before_masking_previous,2))];

	image_before_masking_temp1(:,:,3)=[diff(image_before_masking_temp(:,:,3));zeros(1,size(image_before_masking,2))];



%    	figure
%    	imshow(image_before_masking_previous_temp1(:,:,3), [-10,10])
%  	figure
%  	imshow(image_before_masking_temp1(:,:,3), [-100,100])



	image_difference=image_before_masking_previous_temp1(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:)-image_before_masking_temp1(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:);



	%  figure
	%  imshow(image_difference(:,:,3), [-10,10])
	if (use_second_feature)

	  image_before_masking_previous_temp2(:,:,3)=[diff(image_before_masking_previous_temp(:,:,3),1,2),zeros(size(image_before_masking_previous,1),1)];

	  image_before_masking_temp2(:,:,3)=[diff(image_before_masking_temp(:,:,3),1,2),zeros(size(image_before_masking,1),1)];

%  
%    		figure
%    		imshow(image_before_masking_previous_temp2(:,:,3), [-10,10])
%  		figure
%  		imshow(image_before_masking_temp2(:,:,3), [-100,100])




	  image_difference2=image_before_masking_previous_temp2(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:)-image_before_masking_temp2(1:min([size(image_before_masking,1), size(image_before_masking_previous,1)]),1:min([size(image_before_masking,2), size(image_before_masking_previous,2)]),:);

%  
%  		figure
%  		imshow(image_difference2(:,:,3), [-10,10])

		end
	end









      % before PCA - show a surface depcting the shape disrepency/residual
      
      % size(image_difference)

      if (show_shape_disrepency)
	figure;
	surf(image_difference(:,:,3),'FaceColor', 'red', 'EdgeColor', 'none')
	camlight left;
	lighting phong;
      end 

    % image_difference(:,:,3)

    % apply binary masks and other filters to remove outliers and border-caused artifacts

      if strcmp(residue_filtering_method, 'none') 

      elseif strcmp(residue_filtering_method, 'threshold') % simple thresholding

	baseline=image_difference(size(image_difference(:,:,3),1)/2,size(image_difference(:,:,3),2)/2,3);

	% this can be vetorised
	for k=1:size(image_difference,1)
	  for l=1:size(image_difference,2)
	  %  if (image_difference(k,l,3)>0.3 | image_difference(k,l,3)<-0.3) image_difference(k,l,3)=0; end
	    if (abs(image_difference(k,l,3)-baseline)>carving_threshold) image_difference(k,l,3)=NaN; end % or set to value "baseline"

	    if (use_second_feature)
		if (abs(image_difference2(k,l,3)-baseline)>carving_threshold) image_difference2(k,l,3)=NaN; end
	    end
	  end
	end


      elseif strcmp(residue_filtering_method, 'binary_mask') % remove images near the border of the sides, but not the other areas
	
      if (open_expressions)
	cd (strcat(program_path,'/pcafaces/Binary_masks'));
	% mask=imread('mask1.png');
	mask=imread('mask2.png');
      else
	cd (strcat(program_path,'/pcafaces/Binary_masks'));
	mask=imread('mask-gc2.png');
      end


	      % several masks are available
	for k=1:size(image_difference,1)
	  for l=1:size(image_difference,2)
	    if (mask(k,l)<20)
	      image_difference(k,l,3)=NaN;
	      if (use_second_feature)
		mage_difference2(k,l,3)=NaN;
	      end
	    end
	  end
	end
      else % use all the above

	baseline=image_difference(size(image_difference(:,:,3),1)/2,size(image_difference(:,:,3),2)/2,3);

	% this can be vetorised
	for k=1:size(image_difference,1)
	  for l=1:size(image_difference,2)
	  %  if (image_difference(k,l,3)>0.3 | image_difference(k,l,3)<-0.3) image_difference(k,l,3)=0; end
	    if (abs(image_difference(k,l,3)-baseline)>carving_threshold) image_difference(k,l,3)=NaN; end % or set to value "baseline"
	    if (use_second_feature)
		if (abs(image_difference2(k,l,3)-baseline)>carving_threshold) image_difference2(k,l,3)=NaN; end 
	    end
	  end
	end

	cd (strcat(program_path,'/pcafaces/Binary_masks'));
	if (open_expressions)
	  cd (strcat(program_path,'/pcafaces/Binary_masks'));
	  % mask=imread('mask1.png');
	  mask=imread('mask2.png');
	else
	  cd (strcat(program_path,'/pcafaces/Binary_masks'));
	  mask=imread('mask-gc2.png');
	end
	      % several masks are available
	for k=1:size(image_difference,1)
	  for l=1:size(image_difference,2)
	    if (mask(k,l)<5)
	      image_difference(k,l,3)=NaN;
	      if (use_second_feature)
		mage_difference2(k,l,3)=NaN;
	      end
	    end
	  end
	end
      end

      if (show_shape_disrepency)
	% figure;
	% imshow(mask)
	figure;
	surf(image_difference(:,:,3),'FaceColor', 'red', 'EdgeColor', 'none')
	camlight left;
	lighting phong;
      end 

      if (select_subregion) % if non-empty, then select just a smaller sample of the points cloud

%  figure;
%  imshow(image_difference, [-40 40])
	subregion=image_difference(select_subregion(1):select_subregion(2),select_subregion(3):select_subregion(4), 3);

	if (show_subregion)
	  figure;
	  imshow(subregion, [-100 100])
	end

	for o=1:size(subregion,1)
	  for p=1:size(subregion,2)
	    if ~(subregion(o,p)>-30000)
	      subregion(o,p)=0;
	    end
	  end
	end
	

	if (smooth_residual)
	  subregion=smooth2(subregion,6);
	end

% subregion
% size(subregion)
% figure;
% imshow(subregion, [-100 100])

	subregion=subregion(1:skip_size:end, 1:skip_size:end);
% size(subregion)
% size(image_difference)

	if (use_second_feature)
	  subregion2=image_difference2(select_subregion(1):select_subregion(2),select_subregion(3):select_subregion(4), 3);

	  if (show_subregion)
	    figure;
	    imshow(subregion2, [-100 100])
	  end

	  for o=1:size(subregion2,1)
	    for p=1:size(subregion2,2)
	      if ~(subregion2(o,p)>-30000)
		subregion2(o,p)=0;
	      end
	    end
	  end
	  

	  if (smooth_residual)
	    subregion2=smooth2(subregion2,6);
	  end


	  subregion2=subregion2(1:skip_size:end, 1:skip_size:end);

	end

	if (matching_mode==0)
	  if (use_second_feature)
	    vectorised_residuals(image_n/2,:)=[reshape(subregion, 1, size(subregion,1)*size(subregion,2)), reshape(subregion2, 1, size(subregion2,1)*size(subregion2,2))];
	  else
	    vectorised_residuals(image_n/2,:)=reshape(subregion, 1, size(subregion,1)*size(subregion,2));
	  end
	end
      else
	if (matching_mode==0)
	  if (use_second_feature)
	    vectorised_residuals(image_n/2,:)=[reshape(image_difference(:,:,3), 1, size(image_difference,1)*size(image_difference,2)), reshape(image_difference2(:,:,3), 1, size(image_difference2,1)*size(image_difference2,2))];
%  size(vectorised_residual)
%  size(image_difference2)
	  else

	    %  image_difference(:,:,3)
	    %  size(image_difference,2)
	    %  size(image_difference,1)
	    %  image_n/2
	    vectorised_residuals(image_n/2,:)=reshape(image_difference(:,:,3), 1, size(image_difference,1)*size(image_difference,2));
	  end
	end
      end

    %      for k=1:size(subregion,1)
    %        for l=1:size(subregion,2)
    %  	image_difference(k,l,3)=subregion(k,l);
    %        end
    %      end

	% size(image_difference)

	% vectorise the residuals

	% size(image_difference)
	% size(reshape(image_difference(:,:,3), 1, size(image_difference,1)*size(image_difference,2)))
	
      if (matching_mode)
	% rather than open one image from the other set, open a whole gallery of them and find the best match
	
%  size(vectorised_residuals,1)
%  size(vectorised_residuals)
%  size(reshape(image_difference(:,:,3), 1, size(image_difference,1)*size(image_difference,2)))

	% see the effect of adding the new residue to the PCA model

	% pc(1,:)
	% pc(:,1)

%  size(pc)
%  size(vectorised_residuals,1)
%  size(vectorised_residuals,2)
%  size(subregion)
%  size(score)
%  size(latent)



%  	shape_mode=reshape(pc(:,10), size(subregion,1), size(subregion,2));
%  figure;
%  	surf(shape_mode)
%  	shape_mode=reshape(pc(10,:), size(subregion,1), size(subregion,2));
%  figure;
%  	surf(shape_mode)
%  
%  	shape_mode=reshape(score(10,:), size(subregion,1), size(subregion,2));
%  figure;
%  	surf(shape_mode)
%  figure;
%  plot(latent)


% figure;
% surf(image_difference(:,:,3))

%  figure;
%  shape_mode=reshape(vectorised_residuals(1,:), size(subregion,1), size(subregion,2));
%  surf(shape_mode)

%  figure;
%  vectorised_residuals(10,:)=vectorised_residuals(10,:).*pc(1,:);
%  shape_mode=reshape(vectorised_residuals(10,:), size(subregion,1), size(subregion,2));
%  surf(shape_mode)

	vectorised_residuals=vectorised_residuals_original; % restore to original model

	if (select_subregion) 
	  if (use_second_feature)
	    vectorised_residuals(size(vectorised_residuals,1)+1,:)=[reshape(subregion, 1, size(subregion,1)*size(subregion,2)), reshape(subregion2, 1, 	size(subregion2,1)*size(subregion2,2))];
	  else
	    vectorised_residuals(size(vectorised_residuals,1)+1,:)=[reshape(subregion, 1, size(subregion,1)*size(subregion,2))];
	  end
	else
	  if (use_second_feature)
	    vectorised_residuals(size(vectorised_residuals,1)+1,:)=[reshape(image_difference(:,:,3), 1, size(image_difference,1)*size(image_difference,2)), reshape(image_difference2(:,:,3), 1, size(image_difference2,1)*size(image_difference2,2))];
	  else
	    vectorised_residuals(size(vectorised_residuals,1)+1,:)=[reshape(image_difference(:,:,3), 1, size(image_difference,1)*size(image_difference,2))];
	  end
	end

        if (perform_pca)
	  % vectorised_residuals
	  [pc_prime,score_prime,latent_prime,tsquare_prime] = princomp(vectorised_residuals);
	  display('New model built with probe residual');

	  if strcmp(model_comparator, 'mean_abs')
	    match_score_concat(image_n/2)=mean(mean(abs(pc_prime-pc)))
	  elseif strcmp(model_comparator, 'median_sqr') 
	    match_score_concat(image_n/2)=median(median(abs((pc_prime-pc)^2)))
	  elseif strcmp(model_comparator, 'mean_sqr') 
	    match_score_concat(image_n/2)=mean(mean(abs((pc_prime-pc)^2)))
	  elseif strcmp(model_comparator, 'median_abs') 
	    match_score_concat(image_n/2)=median(median(abs(pc_prime-pc)))

% size(score_prime)
% size(latent_prime)
	  elseif strcmp(model_comparator, 'scores') 
	    match_score_concat_score(image_n/2)=median(median(abs(score_prime(86,:))))
	  elseif strcmp(model_comparator, 'determinant')
	    if (n_eigenvalues)
	      match_score_concat_latent(image_n/2)=prod(latent_prime(1:n_eigenvalues)+delta)
	    else
	      match_score_concat_latent(image_n/2)=prod(latent_prime+delta)
	    end

%  	    pc_diff=pc_prime-pc;
%  	    pc_diff=p_diff+delta;
%  	    match_score_concat(image_n/2)=pc_diff*pc_diff;
	  else % default
            match_score_concat(image_n/2)=median(median(abs(pc_prime-pc)))
	  end

	  min_value=findit(match_score_concat', 2);

	  display(strcat(['Best concatenated match so far is image #', num2str(min_value(1))]));
	end



	cd(strcat(program_path,'/MAT'));		
	save lastscores match_score_concat match_score_mean


	% tsquare
	% size(tsquare)
%  	size(pc)
%  	pc(1:30,1:30)	
%  	size(score)
%   	size(pc)



%  	if (select_subregion) % use basic formulation from the IJCV paper
%  
%  pc*((pc'*pc).^-1)*pc'
%  figure;
%  imshow(pc*((pc'*pc).^-1)*pc');
%  
%  	  match_score_expression_space=vectorised_residual-pc*((pc'*pc).^-1)*pc'*subregion;
%  	else
%  	  match_score_expression_space=vectorised_residuals-pc*(pc'*pc).^-1*pc'*vectorised_residuals;
%  	end



	% size(reshape(vectorised_residuals(image_n/2,:), size(image_difference,1), size(image_difference,2)))
	% difference_mask=reshape(vectorised_residuals(image_n/2,:), size(image_difference,1), size(image_difference,2));
	% match_score(n)=mean(mean(difference_mask))

	if (select_subregion) 
	  % match_score_mean(image_n/2)=abs(mean(mean(subregion))) % get a score of interest for debugging purposes
	  % match_score_mean(image_n/2)=abs(median(median(subregion.^2))) 

	  match_score_mean(image_n/2)=abs(mean(mean(subregion.^2))) 


%    	  figure;
%    	  imshow(subregion, [-40 40]);
	else
	  match_score_mean(image_n/2)=abs(mean(mean(image_difference(:,:,3))))
	end

	min_value=findit(match_score_mean', 2);
	display(strcat(['Best mean match so far is image #', num2str(min_value(1))]));

	% take the discrepency/residuals, reduce distances using search around the model (but avoid over-fitting), then measure cloudpoint differences

      else % just register and concatenate for PCA later
	  if (open_expressions)
	    display(' ');
	    display(strcat(['	- Registration finished for file:'],[neutrals_fileName],['::'],[num2str(image_n)], [' of '], [' '], [num2str(n_images_to_open*2)]));
	  end
      end

    else % save image for the pairing
      combined_data_previous=combined_data; % remember masked part (for the following interation)
      image_before_masking_previous=image_before_masking;
    end

  end
  if (matching_mode)
    % display best match and zero the results for iterations over probes in later experiments
    display(strcat(['Best man match overall is image #', min_value(1)]));
  end

  if (matching_mode==0)
    % deal with NaN issues by converting them all to zeroes (inelegant workaround)

    for k=1:size(vectorised_residuals,1)
      for l=1:size(vectorised_residuals,2)
	if ~(vectorised_residuals(k,l)>-30000)
	  vectorised_residuals(k,l)=0;
	end
      end
    end


    % use the data acquired from several difference shape residuals to start building a model based on PCA

    if (perform_pca)
      % vectorised_residuals
      % [pc,score,latent,tsquare] = princomp(vectorised_residuals(2,100));
      [pc,score,latent,tsquare] = princomp(vectorised_residuals);
      % pc % principal components
      latent(1:10) % bebugging - top 10 eigenvectors
    end
    cumulative_sum_numbers=cumsum(latent)./sum(latent);
    percent_explained = 100*latent(1:10)/sum(latent(1:10))

    figure;
    pareto(percent_explained)
    xlabel('Principal Component')
    ylabel('Variance Explained (%)')

      % save to disk the files dealt with or the model, for subsequent analysis
    cd(strcat(program_path,'/MAT'));
    save pcaexperiment pc score latent tsquare vectorised_residuals n_images_to_open
  end
end % ended training

  % start using the model matching
  display('Run this function in matching mode')

output='Operation succeeded'; % no error found (i.e. completion)
