%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  File: trainEigenEyes.m
%  Matlab script file
%  Date: Oct, 08
%
% Train a PCA model for eye images.  The model of a given
% image is in the form
%        im = sum_k U_k c_k 
% where U_k, k = 1 to nBasis are PCA directions
% The image im has been transformed such that: 
%   im = (I - a * dcIm) * scl
% where dcIm = ones(ny, nx), sum(sum(im)) = 0, sum(sum(im.^2)) = nx*ny
% and size(I) = size(im) = [ny, nx].

% Dependencies, Toolboxes:
%      iseToolbox/
%      utvisToolbox/file/
% Data file: trainSet.mat

% Things still to do:
%   Train multiple PCA patches (multiple means)
%   SPCA

% Author: ADJ
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%% Check Path and Constants  %%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
close all;

if ~exist('rescaleImageVectors','file')
  addpath('./util');
end

%  Check path.  If these aren't on your path, then you will need
%  to add these toolboxes to your path.  See ~jepson/pub/matlab/startup.m
which showIm          % should be iseToolbox\pyrTools\showIm.m
which histo           % should be iseToolbox\pyrTools\MEX\histo.*

%% %%%%%%  Initialize random number generator %%%%%%%%%%%%%%%%%%%%%%%
% Random number generator seed:
seed = round(sum(1000*clock));
rand('state', seed);
seed0 = seed;


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
BATCH = false; % true for no display, but save results in svdEyes.mat
SMOOTH = false; % prefilter images.

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load training data sets.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('trainSet');
nEye = size(eyeIm,2);
clear nonIm;
testImSz = sizeIm;

%% %%%%%%%%%
% Normalize training set.
% %%%%%%%%%%

if SMOOTH % Preblur the eye images
  sigma = 2/pi;
  eyeIm = blurImageVectors(eyeIm, testImSz, sigma);
end

% For each image, remove mean gray-level and scale so that the
% pixel variance is 1.
eyeIm = rescaleImageVectors(eyeIm);

% Check normalization
% Should be close to 0
max(abs(sum(eyeIm,1)))
% Should be 1 
[mx mn] = range2(sqrt(sum(eyeIm.^2,1)/prod(testImSz)))


% %%%%%%%%% Display some of the rescaled eye images %%%%%%%%%%%%%%%%%%%%%%%%
if ~BATCH
  figure(1); clf;
  for k=1:20
    b = 1 + floor(nEye * rand(1,1));
    showIm(reshape(eyeIm(:,b), sizeIm));
    title(sprintf('Eye %d', b));
    fprintf(2, 'Press any key to continue...');
    pause; fprintf(2, 'ok\n');
  end
end

%% %%%%%%%%%
% Compute mean image.
% %%%%%%%%%%

meanIm = sum(eyeIm,2)/nEye;

if ~BATCH % Show the mean image.
  figure(1); clf;
  showIm(reshape(meanIm, sizeIm));
  title('Mean Eye Image');
  fprintf(2, 'Press any key to continue...');
  pause; fprintf(2, 'ok\n');
end  



%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remove meanIm from training set
diffIm = eyeIm - repmat(meanIm, 1, nEye);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the SVD on the brightness normalized, difference from mean images.
% %%%%%%
fprintf(2, 'Calling svd on a %d x %d matrix...', size(eyeIm));
if (size(eyeIm,2) > size(diffIm,1))
  [v s basisIm] = svd(diffIm',0);
else
  [basisIm s v] = svd(diffIm,0);
end
fprintf(2, 'done\n');
sigmaUniEye = diag(s)/sqrt(size(eyeIm,2));
clear v s
sigma = sigmaUniEye;

%% %%%%%%
% Save the results
% %%%%%%%
if BATCH  % Don't overwrite by default.
  save 'svdEyes' sizeIm meanIm sigmaUniEye basisIm;
end

if ~BATCH
  browseEigen(sizeIm, meanIm, basisIm, sigmaUniEye, eyeIm);
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%% End: Train Eigen-Eyes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
