
%%% This is a tutorial and exercises for linear visual filters.
%%% This part includes learning filters from an image ensemble
%%% using linear transforms such as PCA and ICA.

%%% Make sure to add paths of all directories associated with
%%% the tutorial. Also, note that some functions used are
%%% similar to built in Matlab functions -- you can use these
%%% or your own versions.

%%% Read the comments and copy each line of code into matlab.
%%% Note that in some places the code is incomplete and you
%%% need to fill in the pieces...
%%% Type 'help <function_name>' in Matlab window for any
%%% function you would like more information on.
%%% Type 'which <function_name>' in Matlab window for the
%%% location of a file.

%%% Odelia Schwartz 2016; previously used in Berkeley Summer Course 2015

%% Add some paths 
path(path,'Images')
path(path,'imageica')
path(path,'imageica/code')
path(path,'imageica/results')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% (1) Optimized filters from images

%%% To optimize from scenes, we will work with a fixed filter
%%% size of 12 by 12. Rather than looking at a full image, we
%%% will consider many samples of 12 by 12 patches from the image.

%%% We will start out by taking many samples of random
%%% patches from an image, and then do a linear transformation
%%% (e.g., PCA, ICA) on the patches.

%%% Read in an image -- try different examples
multFact = 256;
im=vanRead('imk00649.iml');
im = multFact*(im./max(max(im)));
im = im(1:512,1:512);
size(im)
showIm(im)

%%% *To do": Write a Matlab function 
%%% sig = get_patches(im, patch_dim, num_patches)
%%% that takes num_patches random patches from an 
%%% input image and saves each patch as a column
%%% of the output matrix sig.
patch_dim = 12;
num_patches = 10000;
sig = get_patches(im, patch_dim, num_patches);
size(sig)


%%%% *To do*: Look at example random patches in sig by reshaping
%%%% the columns of the matrix to 12 by 12. For example column
%%%% 2 of sig is in Matlab sig(:,2) and you can use Matlab's 
%%%% reshape command.
 

%%% PCA:
%%% *To do *: Consider linear transformations on the image patch ensemble
%%% of the form: newsig = mat * sig        (where * is multiplication)
%%% First compute a PCA transform of the patches sig. 
%%% You can use the function mypca.m, in which you input sig and 
%%% obtain as the output newsig, mat (and also the eigenvectors
%%% and eigenvalues). Do a help on mypca.m


%%% *To do*: Confirm that the transform produced newsig that are
%%% uncorrelated. 


%%% *To do*: Look at the sizes of sig, newsig, and mat.
%%% Plot the eigenvalues.
%%% Look at the principle components (eigenvectors):
%%% The rows of the matrix mat are called projection functions.
%%% They project from the input columns (patches of image)
%%% onto the principle components.
%%% We can view these by reshaping each row (mat(i,:)) to 12 by 12
%%% What do the eigenvectors look like for high eigenvalues? For low ones?


%%% *To do*: We are also interested in the transform:
%%% sig = mat_inverse * newsig
%%% This is a transform from the output of pca back to the
%%% original image patches. The matrix mat_inverse is the
%%% inverse of mat. The columns of the matrix mat_inverse 
%%% are known as basis functions. The input is a linear
%%% combination of the columns of mat_inverse.
%%% To view the basis functions, we need to view columns of mat_inverse
%%% using matlab's inv(mat). Look at different columns of mat_inverse.

mat_inverse = inv(mat);


%%% Since pca is an orthonormal transform mat*mat_inverse is a diagonal
%%% matrix. Note that if we use pca without whitening then the basis
%%% functions are exactly given by the columns of the transpose
%%% of mat since the inverse is equal to the transpose.

%%% We will now look at ICA.
%%% There are various versions of ICA available:
%%% Bell and Sejnowski ICA:
%%% http://cnl.salk.edu/~tony/ica.html
%%% Hoyer and Hyvarinen fast ICA applied to images (imageica):
%%% http://www.cs.helsinki.fi/u/phoyer/software.html
%%% And sparse coding:
%%% Olshausen and Field:
%%% http://redwood.berkeley.edu/bruno/sparsenet/

%%% We will use the fastICA code of Hoyer and Hyvarinen
%%% applied to images
%%% http://www.cs.helsinki.fi/u/phoyer/software.html
%%% We have included it here.

sig = get_patches(im, patch_dim, num_patches);

clear X;
global X;
X = sig;

% Subtract local mean from each patch
% Currently commented out
% X = X-ones(size(X,1),1)*mean(X);

% Pre Whiten
covarianceMatrix = X*X'/size(X,2);
[E, D] = eig(covarianceMatrix);
[dummy,order] = sort(diag(-D));
E = E(:,order);
d = diag(D);
d = real(d.^(-0.5));
D = diag(d(order));
X = D*E'*X;
whiteningMatrix = D*E';
dewhiteningMatrix = E*D^(-1);

p.seed = 1;
p.write = 20;
p.model = 'ica';
p.algorithm = 'fixed-point';
p.components = min(size(X));
numIters = 500;            % run for 500 trials
estimateModif( whiteningMatrix, dewhiteningMatrix, 'tryica.mat', p, numIters );
%%% Original function is estimate.m which runs infinitely; 
%%% estimateModif.m runs for numIters times

%%% Load saved results
load('tryica.mat');

%%% *To do*: The columns of A contain the basis functions; view them...


