Index: proiecte/PPPP/ica/Arch1.m
===================================================================
--- proiecte/PPPP/ica/Arch1.m	(revision 94)
+++ proiecte/PPPP/ica/Arch1.m	(revision 94)
@@ -0,0 +1,69 @@
+% script Arch1.m
+% Finds ICA representation of train and test images under Architecture I, 
+% described in Bartlett & Sejnowski (1997, 1998), and Bartlett, Movellan & 
+% Sejnowski (2002):  In Architecture I, we load N principal component 
+% eigenvectors into rows of x, and then run ICA on x.
+%
+% Put the aligned training images in the rows of C, one image per row.  
+% In the following examples, there are 500 images of aligned faces of size 
+% 60x60 pixels, so C is 500x3600. 
+%
+% You can use the following matlab code to create C:
+% markFeatures.m collects eye and mouth positions. 
+% align_Faces.m crops, aligns, and scales the face images.
+% loadFaceMat.m loads the images into the rows of C. 
+%
+% This script also calls the matrix of PCA eigenvectors organized in 
+% the columns of V (3600x499), created by [V,R,E] = pcabigFn(C');
+%
+% The ICA representation will be in the rows of F (called B in Bartlett, 
+% Movellan & Sejnowski, 2002): 
+
+[V,R,E] = pcabigFn(C');
+%D = zeroMn(C')'; % D is 500x3600 and D = C-ones(500,1)*mean(C);
+%R = D*V; 	 % R is 500x499 and contains the PCA coefficients;
+
+% We choose to use the first 200 eigenvectors. 
+% (If PCA generalizes better by dropping first few eigenvectors, ICA will too).
+
+%x = V(:,1:200)'; 		% x is 200x3600
+x = V(:,1:3)'; 		% x is 200x3600
+runica 				% calculates wz, w and uu. The matrix x gets 
+				% overwritten by a sphered version of x. 
+%F = R(:,1:200) * inv(w*wz); 	% F is 500x200 and each row contains the 
+F = R(:,1:3) * inv(w*wz); 	% F is 500x200 and each row contains the 
+			    	% ICA1 rep of an image
+
+% Representations of test images under architecture I: 
+% Put original aligned test images in rows of Ctest. 
+
+Dtest = zeroMn(Ctest')'; % For proper testing, subtract the mean of the 
+			 % training images not the test images: 
+			 % Dtest = Ctest-ones(500,1)*mean(C);
+Rtest = Dtest*V;
+%Ftest = Rtest(:,1:200) * inv(w*wz);
+Ftest = Rtest(:,1:3) * inv(w*wz);
+
+% Test nearest neighbor classification using cosine, not euclidean distance, 
+% as similarity measure.
+%
+% First create label vectors. These are column vectors of integers. Lets 
+% say our 500 training examples consisted of 500 different people. Then
+% trainClass = [1:500]'; 
+%
+% We also need the correct class labels of the test examples if we want to 
+% compute percent correct. Lets say the test examples were two images each 
+% of the first 10 individuals. Then 
+% testClass = [1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10]';
+
+trainClass = [1:3]';
+testClass = [1 2 3]';
+
+%We now compute percent correct:
+train_ex = F';
+test_ex = Ftest';
+[pc,rankmat] = nnclassFn(train_ex,test_ex,trainClass,testClass);
+
+%pc is percent correct of first nearest neighbor.
+%rankmat gives the top 30 matches for each test image. 
+
Index: proiecte/PPPP/ica/Arch2.m
===================================================================
--- proiecte/PPPP/ica/Arch2.m	(revision 94)
+++ proiecte/PPPP/ica/Arch2.m	(revision 94)
@@ -0,0 +1,60 @@
+% script Arch2.m
+% Finds ICA representation of train and test images under Architecture II, 
+% described in Bartlett & Sejnowski (1997, 1998), and Bartlett, Movellan & 
+% Sejnowski (2002):  In Architecture II, we load N principal component coefficients
+% into rows of x, and then run ICA on x.
+%
+% Put aligned training images in the rows of C, one image per row.  
+% In the following examples, there are 500 images of aligned faces of size 
+% 60x60 pixels, so C is 500x3600. 
+%
+% You can use the following matlab code to create C:
+% markFeatures.m collects eye and mouth positions. 
+% align_Faces.m crops, aligns, and scales the face images.
+% loadFaceMat.m loads the images into the rows of C. 
+%
+% This script also calls the matrix of PCA eigenvectors organized in 
+% the columns of V (3600x499), created by [V,R,E] = pcabigFn(C');
+%
+% The ICA representation will be in F (called U in Bartlett, Movellan & 
+% Sejnowski, 2002): 
+
+[V,R,E] = pcabigFn(C');
+%D = zeroMn(C')'; % D is 500x3600 and D = C-ones(500,1)*mean(C);
+%R = D*V; 	 % R is 500x499 and contains the PCA coefficients;
+
+x = R(:,1:200)'; 	% x is 200x500;
+runica 			% calculates w, wz, and uu. The matrix x gets overwritten
+			% by a sphered version of x. 
+F = uu'; 		% F is 500x200 and each row contains the ICA2 rep of 1 image. 
+			% F = w * wz * zeroMn(R(:,1:200)')'; is the same thing.
+
+% Representations of test images under architecture II
+% Put original aligned test images in rows of Ctest:
+
+Dtest = zeroMn(Ctest')'; % For proper testing, subtract the mean of the 
+			 % training images not the test images: 
+			 % Dtest = Ctest-ones(500,1)*mean(C);
+Rtest = Dtest*V;
+Ftest = w * wz * zeroMn(Rtest(:,1:200)')';
+
+% Test nearest neighbor classification using cosine, not euclidean distance, 
+% as similarity measure.
+%
+% First create label vectors. These are column vectors of integers. Lets 
+% say our 500 training examples consisted of 500 different people. Then
+% trainClass = [1:500]'; 
+%
+% We also need the correct class labels of the test examples if we want to 
+% compute percent correct. Lets say the test examples were two images each 
+% of the first 10 individuals. Then 
+% testClass = [1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10]';
+
+%We now compute percent correct:
+train_ex = F';
+test_ex = Ftest';
+[pc,rankmat] = nnclassFn(train_ex,test_ex,trainClass,testClass);
+
+%pc is percent correct of first nearest neighbor.
+%rankmat gives the top 30 matches for each test image. 
+
Index: proiecte/PPPP/ica/README
===================================================================
--- proiecte/PPPP/ica/README	(revision 94)
+++ proiecte/PPPP/ica/README	(revision 94)
@@ -0,0 +1,83 @@
+%---------------------------------------------------------
+% 
+% Copyright 1997 Marian Stewart Bartlett
+% This may be copied for personal or academic use.
+% For commercial use, please contact Marian Bartlett 
+% (marni@salk.edu) for a commercial license.
+%
+% Image representations by Marian Bartlett. Revised 7/14/03
+% ICA code by Tony Bell.
+% The ICA method is patented by Bell & Sejnowski, at the Salk Institute. 
+%
+% Please cite Bartlett, M.S. (2001) Face Image Analysis by 
+% Unsupervised Learning. Boston: Kluwer Academic Publishers. 
+%
+% --------------------------------------------------------
+
+
+Based on Bartlett, Movellan, & Sejnowski (2002). Face Recognition by
+Independent Component analysis. IEEE Transactions on Neural Networks
+13(6) p. 1450-1464, and
+
+Bartlett, M.S. (2001) Face Image Analysis by Unsupervised
+Learning. Boston: Kluwer Academic Publishers.
+
+This directory contains 2 matlab scripts for finding the ICA representation 
+of a set of images for recognition:
+
+1. Arch1.m: Gets representation of train and test images under architecture I
+2. Arch2.m: Gets representation of train and test images under architecture II
+
+Read through the comments of these scripts before attempting to run them. 
+
+The above scripts call the following 6 MATLAB files for running infomax ica.
+Written by Tony Bell http://www.cnl.salk.edu/~tony/
+
+1. runica.m, the ica training script which calls the functions below.
+2. sep96.m, the code for one learning pass thru the data
+3. sepout.m, for optional text output
+4. wchange.m, tracks size and direction of weight changes 
+5. spherex.m, spheres the training matrix x.
+6. zeroMn.m: Returns a zero-mean form of the matrix X, where each row has 
+             zero-mean. (This one was added by Marian Bartlett) 
+
+The following variables are used to calculate ica:
+
+sweep:    how many times you've gone thru the data
+P:        how many timepoints in the data
+N:        how many input (mixed) sources there are
+M:        how many outputs you have
+L:        learning rate
+B:        batch-block size (ie: how many presentations per weight update.)
+t:        time index of data
+sources:  NxP matrix of the N sources you read in
+x:        NxP matrix of mixtures
+u:        MxP matrix of hopefully unmixed sources
+a:        NxN mixing matrix
+w:        MxN unmixing matrix (actually w*wz is the full unmixing matrix
+          in this case)
+wz:       zero-phase whitening: a matrix used to remove 
+          correlations from between the mixtures x. Useful as a 
+          preprocessing step.
+noblocks: how many blocks in a sweep;
+oldw:     value of w before the last sweep
+delta:    w-oldw
+olddelta: value of delta before the last sweep
+angle:    angle in degrees between delta and olddelta
+change:   squared length of delta vector 
+Id:       an identity matrix
+permute:  a vector of length P used to scramble the time order of the
+          sources for stationarity during learning.
+
+INITIAL w ADVICE: identity matrix is a good choice, since, for prewhitened
+data, there will be no distracting initial correlations, and the output
+variances will be nicely scaled so <uu^T>=4I, right size to fit the 
+logistic fn (more or less).
+
+LEARNING RATE ADVICE: 
+N=2: L=0.01 works
+N=8-10: L=0.001 is stable. Run this till the 'change' report settles
+down, then anneal a little. L=0.0005,0.0002,0.0001 etc, a few passes
+(= a few 10,000's of data vectors) each pass.
+N>100: L=0.001 works well on sphered image data.
+
Index: proiecte/PPPP/ica/alignFs.asv
===================================================================
--- proiecte/PPPP/ica/alignFs.asv	(revision 94)
+++ proiecte/PPPP/ica/alignFs.asv	(revision 94)
@@ -0,0 +1,99 @@
+%Script align_Faces
+%
+% Aligns the eye positions of a directory of face images. Reads in Labels.mat, 
+% obtained using getLabels.m, writes jpegs to a specified directory.
+% Specify directory paths at the top of the file. 
+%LabelDir is where Labels.mat is.
+%imgDir is where the original images are
+%DestDir is where you want the cropped jpegs to go. 
+
+function [] = alignFs(LabelDir, imgDir, DestDir, labelname);
+
+%LabelDir = 'C:\Users\Andrei\Desktop\yalefaces\'
+%imgDir = 'C:\Users\Andrei\Desktop\yalefaces\normal\'
+%DestDir = 'C:\Users\Andrei\Desktop\yalefaces\'
+  
+homeDir = pwd; 
+
+% CHANGE THESE VARIABLES AS NEEDED
+%XSIZE =  YSIZE = 	%Size of desired cropped image
+%EYES = %Number of pixels desired between the eyes
+%TEETH_EYES = %Desired no. of pixels from teeth to eyes. 
+
+XSIZE = 300; YSIZE = 200;      
+EYES = 50;
+TEETH_EYES = 80; 
+%SEE ALSO PARAMETERS IN CROP ROUTINE. EYES BELOW MIDPOINT.
+ 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Here is where we load the images and do the preprocessing
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      cd (LabelDir)
+      load(labelname)
+      cd (imgDir)
+      r = dir;
+
+for i = 3:(size(r,1))
+   imgName = r(i).name
+   %[X,map] = imread([ t ]);
+      I=imread(imgName);
+      %I=rgb2gray(I);
+
+     %Extract face measurements:
+     %dxeyes =  x distance between eyes in original image                     
+             % Average the locations of inner & outer corners.
+     %dyeyes =  y distance between eyes in original image 
+     %dEeyes = Euclidean norm distance between eyes
+     %dEteeth_eyes = Euclidean norm distance from teeth to midpt b/neyes
+                 
+     [height, width] = size(I);
+     dxeyes = marks(i-2,3) - marks(i-2,1); %Check these. 
+     dyeyes = marks (i-2,4) - marks(i-2,2);
+     dEeyes = sqrt(dxeyes^2 + dyeyes^2);
+     mean_eye_x = mean([marks(i-2,1), marks(i-2,3)]);
+     mean_eye_y = mean([marks(i-2,2), marks(i-2,4)]);
+     dEteeth_eyes = sqrt((marks(i-2,5)-mean_eye_x)^2 + (marks(i-2,6)-mean_eye_y)^2); 
+                                          
+     %scale
+     yscale = TEETH_EYES / dEteeth_eyes; xscale = EYES / dEeyes;
+     height_new = yscale*height; width_new = xscale*width;
+     tmp0=imresize(I,[height_new,width_new],'bicubic');
+   
+      %rotate (Problem: imrotate rotates about the center of the image. 
+      %To avoid losing feature position information, must first center the 
+      %image on the right eye before rotating. 
+      %Then use right eye position to determine cropping.
+      Reye_x = marks(i-2,1);
+      Reye_y = marks(i-2,2);
+
+      %crop a 200x200 window centered on left eye:
+      %Zero-pad to make sure window never falls outside of image. 
+      %W = 100; %Window radius
+      W = 500;  %For bigger images (Gwen's params).
+      padcols = zeros(size(tmp0,1),W); padrows = zeros(W,size(tmp0,2)+W);
+      padcols = uint8(padcols); padrows=uint8(padrows);
+      tmp = [padrows;padcols,tmp0];
+
+      tmpx = xscale*Reye_x - W +W; tmpy = yscale*Reye_y - W +W;     
+      tmp1 = imcrop(tmp,[tmpx,tmpy,2*W,2*W]);
+      %figure(2);imshow(tmp1)
+
+      angle = 180/pi*atan((yscale*dyeyes)/(xscale*dxeyes));
+      tmp2 = imrotate(tmp1,angle,'bicubic','crop');
+      %figure(2); imshow(tmp2);
+
+      %crop
+      % x and y give the upper left corner of cropped image
+      % Reye is centered at (W,W) = (100,100).
+      % For bigger images (W,W) = (500,500)
+      x = W - (XSIZE-EYES)/2;
+      %y = W - YSIZE/2;   %Eyes at midpoint
+      y = W - YSIZE*1/3;  %Face box
+      tmp3=imcrop(tmp2,[x,y,XSIZE,YSIZE]);
+      figure(1); imshow(tmp3);
+
+      %save
+      [imgName, R] = strtok(imgName, '.');
+      fname = [DestDir,imgName, '.pgm'];
+      imwrite(tmp3,fname,'pgm') 
+end
Index: proiecte/PPPP/ica/alignFs.m
===================================================================
--- proiecte/PPPP/ica/alignFs.m	(revision 94)
+++ proiecte/PPPP/ica/alignFs.m	(revision 94)
@@ -0,0 +1,99 @@
+%Script align_Faces
+%
+% Aligns the eye positions of a directory of face images. Reads in Labels.mat, 
+% obtained using getLabels.m, writes jpegs to a specified directory.
+% Specify directory paths at the top of the file. 
+%LabelDir is where Labels.mat is.
+%imgDir is where the original images are
+%DestDir is where you want the cropped jpegs to go. 
+
+function [] = alignFs(LabelDir, imgDir, DestDir, labelname);
+
+%LabelDir = 'C:\Users\Andrei\Desktop\yalefaces\'
+%imgDir = 'C:\Users\Andrei\Desktop\yalefaces\normal\'
+%DestDir = 'C:\Users\Andrei\Desktop\yalefaces\'
+  
+homeDir = pwd; 
+
+% CHANGE THESE VARIABLES AS NEEDED
+%XSIZE =  YSIZE = 	%Size of desired cropped image
+%EYES = %Number of pixels desired between the eyes
+%TEETH_EYES = %Desired no. of pixels from teeth to eyes. 
+
+XSIZE = 300; YSIZE = 200;      
+EYES = 50;
+TEETH_EYES = 80; 
+%SEE ALSO PARAMETERS IN CROP ROUTINE. EYES BELOW MIDPOINT.
+ 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Here is where we load the images and do the preprocessing
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      cd (LabelDir)
+      load(labelname)
+      cd (imgDir)
+      r = dir;
+
+for i = 3:(size(r,1))
+   imgName = r(i).name
+   %[X,map] = imread([ t ]);
+      I=imread(imgName);
+      %I=rgb2gray(I);
+
+     %Extract face measurements:
+     %dxeyes =  x distance between eyes in original image                     
+             % Average the locations of inner & outer corners.
+     %dyeyes =  y distance between eyes in original image 
+     %dEeyes = Euclidean norm distance between eyes
+     %dEteeth_eyes = Euclidean norm distance from teeth to midpt b/neyes
+                 
+     [height, width] = size(I);
+     dxeyes = marks(i-2,3) - marks(i-2,1); %Check these. 
+     dyeyes = marks (i-2,4) - marks(i-2,2);
+     dEeyes = sqrt(dxeyes^2 + dyeyes^2);
+     mean_eye_x = mean([marks(i-2,1), marks(i-2,3)]);
+     mean_eye_y = mean([marks(i-2,2), marks(i-2,4)]);
+     dEteeth_eyes = sqrt((marks(i-2,5)-mean_eye_x)^2 + (marks(i-2,6)-mean_eye_y)^2); 
+                                          
+     %scale
+     yscale = TEETH_EYES / dEteeth_eyes; xscale = EYES / dEeyes;
+     height_new = yscale*height; width_new = xscale*width;
+     tmp0=imresize(I,[height_new,width_new],'bicubic');
+   
+      %rotate (Problem: imrotate rotates about the center of the image. 
+      %To avoid losing feature position information, must first center the 
+      %image on the right eye before rotating. 
+      %Then use right eye position to determine cropping.
+      Reye_x = marks(i-2,1);
+      Reye_y = marks(i-2,2);
+
+      %crop a 200x200 window centered on left eye:
+      %Zero-pad to make sure window never falls outside of image. 
+      %W = 100; %Window radius
+      W = 500;  %For bigger images (Gwen's params).
+      padcols = zeros(size(tmp0,1),W); padrows = zeros(W,size(tmp0,2)+W);
+      padcols = uint8(padcols); padrows=uint8(padrows);
+      tmp = [padrows;padcols,tmp0];
+
+      tmpx = xscale*Reye_x - W +W; tmpy = yscale*Reye_y - W +W;     
+      tmp1 = imcrop(tmp,[tmpx,tmpy,2*W,2*W]);
+      %figure(2);imshow(tmp1)
+
+      angle = 180/pi*atan((yscale*dyeyes)/(xscale*dxeyes));
+      tmp2 = imrotate(tmp1,angle,'bicubic','crop');
+      %figure(2); imshow(tmp2);
+
+      %crop
+      % x and y give the upper left corner of cropped image
+      % Reye is centered at (W,W) = (100,100).
+      % For bigger images (W,W) = (500,500)
+      x = W - (XSIZE-EYES)/2;
+      %y = W - YSIZE/2;   %Eyes at midpoint
+      y = W - YSIZE*1/3;  %Face box
+      tmp3=imcrop(tmp2,[x,y,XSIZE,YSIZE]);
+      figure(1); imshow(tmp3);
+
+      %save
+      [imgName, R] = strtok(imgName, '.');
+      fname = [DestDir,imgName, '.pgm'];
+      imwrite(tmp3,fname,'pgm') 
+end
Index: proiecte/PPPP/ica/align_Faces.m
===================================================================
--- proiecte/PPPP/ica/align_Faces.m	(revision 94)
+++ proiecte/PPPP/ica/align_Faces.m	(revision 94)
@@ -0,0 +1,97 @@
+%Script align_Faces
+%
+% Aligns the eye positions of a directory of face images. Reads in Labels.mat, 
+% obtained using getLabels.m, writes jpegs to a specified directory.
+% Specify directory paths at the top of the file. 
+%LabelDir is where Labels.mat is.
+%imgDir is where the original images are
+%DestDir is where you want the cropped jpegs to go. 
+
+LabelDir = 'C:\Users\Andrei\Desktop\yalefaces\'
+imgDir = 'C:\Users\Andrei\Desktop\yalefaces\normal\'
+DestDir = 'C:\Users\Andrei\Desktop\yalefaces\'
+  
+homeDir = pwd; 
+
+% CHANGE THESE VARIABLES AS NEEDED
+%XSIZE =  YSIZE = 	%Size of desired cropped image
+%EYES = %Number of pixels desired between the eyes
+%TEETH_EYES = %Desired no. of pixels from teeth to eyes. 
+
+XSIZE = 300; YSIZE = 200;      
+EYES = 50;
+TEETH_EYES = 80; 
+%SEE ALSO PARAMETERS IN CROP ROUTINE. EYES BELOW MIDPOINT.
+ 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Here is where we load the images and do the preprocessing
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+      cd (LabelDir)
+      load Labels
+      cd (imgDir)
+      r = dir;
+
+for i = 3:(size(r,1))
+   imgName = r(i).name
+   %[X,map] = imread([ t ]);
+      I=imread(imgName);
+      %I=rgb2gray(I);
+
+     %Extract face measurements:
+     %dxeyes =  x distance between eyes in original image                     
+             % Average the locations of inner & outer corners.
+     %dyeyes =  y distance between eyes in original image 
+     %dEeyes = Euclidean norm distance between eyes
+     %dEteeth_eyes = Euclidean norm distance from teeth to midpt b/neyes
+                 
+     [height, width] = size(I);
+     dxeyes = marks(i-2,3) - marks(i-2,1); %Check these. 
+     dyeyes = marks (i-2,4) - marks(i-2,2);
+     dEeyes = sqrt(dxeyes^2 + dyeyes^2);
+     mean_eye_x = mean([marks(i-2,1), marks(i-2,3)]);
+     mean_eye_y = mean([marks(i-2,2), marks(i-2,4)]);
+     dEteeth_eyes = sqrt((marks(i-2,5)-mean_eye_x)^2 + (marks(i-2,6)-mean_eye_y)^2); 
+                                          
+     %scale
+     yscale = TEETH_EYES / dEteeth_eyes; xscale = EYES / dEeyes;
+     height_new = yscale*height; width_new = xscale*width;
+     tmp0=imresize(I,[height_new,width_new],'bicubic');
+   
+      %rotate (Problem: imrotate rotates about the center of the image. 
+      %To avoid losing feature position information, must first center the 
+      %image on the right eye before rotating. 
+      %Then use right eye position to determine cropping.
+      Reye_x = marks(i-2,1);
+      Reye_y = marks(i-2,2);
+
+      %crop a 200x200 window centered on left eye:
+      %Zero-pad to make sure window never falls outside of image. 
+      %W = 100; %Window radius
+      W = 500;  %For bigger images (Gwen's params).
+      padcols = zeros(size(tmp0,1),W); padrows = zeros(W,size(tmp0,2)+W);
+      padcols = uint8(padcols); padrows=uint8(padrows);
+      tmp = [padrows;padcols,tmp0];
+
+      tmpx = xscale*Reye_x - W +W; tmpy = yscale*Reye_y - W +W;     
+      tmp1 = imcrop(tmp,[tmpx,tmpy,2*W,2*W]);
+      %figure(2);imshow(tmp1)
+
+      angle = 180/pi*atan((yscale*dyeyes)/(xscale*dxeyes));
+      tmp2 = imrotate(tmp1,angle,'bicubic','crop');
+      %figure(2); imshow(tmp2);
+
+      %crop
+      % x and y give the upper left corner of cropped image
+      % Reye is centered at (W,W) = (100,100).
+      % For bigger images (W,W) = (500,500)
+      x = W - (XSIZE-EYES)/2;
+      %y = W - YSIZE/2;   %Eyes at midpoint
+      y = W - YSIZE*1/3;  %Face box
+      tmp3=imcrop(tmp2,[x,y,XSIZE,YSIZE]);
+      figure(1); imshow(tmp3);
+
+      %save
+      [imgName, R] = strtok(imgName, '.');
+      fname = [DestDir,imgName, '.pgm'];
+      imwrite(tmp3,fname,'pgm') 
+end
Index: proiecte/PPPP/ica/cosFn.m
===================================================================
--- proiecte/PPPP/ica/cosFn.m	(revision 94)
+++ proiecte/PPPP/ica/cosFn.m	(revision 94)
@@ -0,0 +1,14 @@
+% function [S] = cosFn(mat1,mat2),
+% Computes the cosine (normalized dot product) between training vectors in 
+% columns of mat1 and test vectors in columns of mat2. Outputs a matrix of 
+% cosines (similarity matrix). 
+%
+% Written by Ian Fasel.
+
+function [S] = cosFn(mat1,mat2),
+  
+  denom = sum(mat1.^2,1)*sum(mat2'.^2,2);
+  denom (find(denom==0)) = 0.00000000000000000000001;
+  numer = mat1'*mat2;
+  
+  S = numer./denom;
Index: proiecte/PPPP/ica/loadFaceMat.m
===================================================================
--- proiecte/PPPP/ica/loadFaceMat.m	(revision 94)
+++ proiecte/PPPP/ica/loadFaceMat.m	(revision 94)
@@ -0,0 +1,18 @@
+%function C = loadFaceMat(imgdir);
+%
+%Loads a directory of images into the rows of C.
+%imgdir is a string with the path to the directory containing the images.
+%For example:
+%imgdir = 'C:\Documents and Settings\Marian Bartlett\My Documents\TestDir\AlignedFaceImages\'
+
+cd (imgDir)
+r = dir;
+
+C = []
+for i = 3:(size(r,1))
+   t = r(i).name;
+        I=imread(t);
+      tmp=mat2gray(double(I));
+      tmp = reshape(tmp,1,size(tmp,1)*size(tmp,2));
+      C = [C;tmp];
+end
Index: proiecte/PPPP/ica/loadTestMat.m
===================================================================
--- proiecte/PPPP/ica/loadTestMat.m	(revision 94)
+++ proiecte/PPPP/ica/loadTestMat.m	(revision 94)
@@ -0,0 +1,18 @@
+%function C = loadFaceMat(imgdir);
+%
+%Loads a directory of images into the rows of C.
+%imgdir is a string with the path to the directory containing the images.
+%For example:
+%imgdir = 'C:\Documents and Settings\Marian Bartlett\My Documents\TestDir\AlignedFaceImages\'
+
+cd ('C:\Users\Andrei\Desktop\yalefaces\test\')
+r = dir;
+
+Ctest = []
+for i = 3:(size(r,1))
+   t = r(i).name;
+        I=imread(t);
+      tmp=mat2gray(double(I));
+      tmp = reshape(tmp,1,size(tmp,1)*size(tmp,2));
+      Ctest = [Ctest;tmp];
+end
Index: proiecte/PPPP/ica/markFeats.m
===================================================================
--- proiecte/PPPP/ica/markFeats.m	(revision 94)
+++ proiecte/PPPP/ica/markFeats.m	(revision 94)
@@ -0,0 +1,34 @@
+%Script markFeatures
+%For marking eye and mouth positions in face images.
+%Writes a file called Labels.mat, in which each row indexes an image, and
+%the columns are [x,y] positions of subject's right eye, [x,y] left eye,
+%and [x,y] of mouth. Specify the image directory and the destination 
+%directory (Where you want the labels saved) at the top of the script.
+
+function [] = markFeats(imgdir, destdir, fname)
+
+cd (imgdir)
+r = dir;
+
+%get marks
+marks = [];
+for i = 3:(size(r,1))
+    t = r(i).name;
+    [X,map] = imread([ t ]);
+
+    figure(1);
+    colormap gray;
+    if isfloat(X)
+        image(gray2ind(mat2gray((X))));
+    else
+        image(X);
+    end
+    title(t);
+    disp 'Click subjects right eye, left eye, then mouth.'
+    [m,n] = ginput(3); pos = round([m,n]);
+    pos = reshape(pos',1,6);
+    marks = [marks; pos];
+end
+
+cd (destdir)
+save(fname, 'marks', 'r')
Index: proiecte/PPPP/ica/markFeatures.m
===================================================================
--- proiecte/PPPP/ica/markFeatures.m	(revision 94)
+++ proiecte/PPPP/ica/markFeatures.m	(revision 94)
@@ -0,0 +1,36 @@
+%Script markFeatures
+%For marking eye and mouth positions in face images.
+%Writes a file called Labels.mat, in which each row indexes an image, and
+%the columns are [x,y] positions of subject's right eye, [x,y] left eye,
+%and [x,y] of mouth. Specify the image directory and the destination 
+%directory (Where you want the labels saved) at the top of the script.
+
+
+imgdir = 'C:\Users\Andrei\Desktop\yalefaces\normal\'
+destdir = 'C:\Users\Andrei\Desktop\yalefaces\'
+
+cd (imgdir)
+r = dir;
+
+%get marks
+marks = [];
+for i = 3:(size(r,1))
+    t = r(i).name;
+    [X,map] = imread([ t ]);
+
+    figure(1);
+    colormap gray;
+    if isfloat(X)
+        image(gray2ind(mat2gray((X))));
+    else
+        image(X);
+    end
+    title(t);
+    disp 'Click subjects right eye, left eye, then mouth.'
+    [m,n] = ginput(3); pos = round([m,n]);
+    pos = reshape(pos',1,6);
+    marks = [marks; pos];
+end
+
+cd (destdir)
+save Labels marks r
Index: proiecte/PPPP/ica/markFeatures_n.asv
===================================================================
--- proiecte/PPPP/ica/markFeatures_n.asv	(revision 94)
+++ proiecte/PPPP/ica/markFeatures_n.asv	(revision 94)
@@ -0,0 +1,34 @@
+%Script markFeatures
+%For marking eye and mouth positions in face images.
+%Writes a file called Labels.mat, in which each row indexes an image, and
+%the columns are [x,y] positions of subject's right eye, [x,y] left eye,
+%and [x,y] of mouth. Specify the image directory and the destination 
+%directory (Where you want the labels saved) at the top of the script.
+
+function [] = markFeats(imgdir, destdir, fname)
+
+cd (imgdir)
+r = dir;
+
+%get marks
+marks = [];
+for i = 3:(size(r,1))
+    t = r(i).name;
+    [X,map] = imread([ t ]);
+
+    figure(1);
+    colormap gray;
+    if isfloat(X)
+        image(gray2ind(mat2gray((X))));
+    else
+        image(X);
+    end
+    title(t);
+    disp 'Click subjects right eye, left eye, then mouth.'
+    [m,n] = ginput(3); pos = round([m,n]);
+    pos = reshape(pos',1,6);
+    marks = [marks; pos];
+end
+
+cd (destdir)
+save fname marks r
Index: proiecte/PPPP/ica/nnClassFn.m
===================================================================
--- proiecte/PPPP/ica/nnClassFn.m	(revision 94)
+++ proiecte/PPPP/ica/nnClassFn.m	(revision 94)
@@ -0,0 +1,60 @@
+%function [testPerf,rankmat,rank] = nnclassFn(train,test,trainClass,answer)
+%
+%Reads in training examples, test examples, class labels of training
+%examples, and correct class of test examples. Data are in columns of train
+%and test, and labels are column vectors.
+% 
+%Note: You will need to create label vectors. TrainClass is a column 
+%vector of integers indicating the identity of the training examples. 
+%e.g. for faces of 3 people with two views each, TrainClass = [1 1 2 2 3 3 ]';
+%Answer contains the correct labels of the test images, which enables
+%us to compute percent correct.  
+%
+%Gets matrix of normalized dot products. Outputs nearest neighbor
+%classification of test examples and percent correct.
+%rankmat gives the top 30 matches for each test image.  rank is a vector
+%containing the percent of times the correct match is in the top N matches.
+
+
+function [testPerf,rankmat,rank] = nnclassFn(train,test,trainClass,answer);
+
+numTest = size(test,2);
+numTrain = size(train,2);
+
+%Get distances to training examples
+%dists = eucDist(test,train); %Outputs a Ntest x Ntrain matrix of Euc dist
+dists=-1 * cosFn(test,train);%Outputs a Ntest x Ntrain matrix of cosines
+
+%sort the rows of dists to find the nearest training example:
+[Sdist,nearest] = sort(dists'); %cols of Sdist are distances in ascend order
+        %1st row of nearest is index of 1st closest training example
+
+%Create vector with nearest example, and vector with class label.
+Nnbr = nearest(1,:);            %First row of nearest contains NN
+%Nnbr = nearest(2,:);
+testClass = trainClass(Nnbr);
+
+correct = find( (testClass - answer == 0));
+testPerf = size(correct,1) / size(answer,1)
+if(size(correct,2)>size(correct,1))
+        testPerf = size(correct,2) / size(answer,2)
+        'check vector orientation'
+end
+
+%get rank = %correct in top N:
+cumtestPerf=0;
+for i = 1:3
+        rankmat(:,i) = trainClass(nearest(i,:)');
+        correcti = find( (rankmat(:,i) - answer == 0));
+        cumtestPerf = cumtestPerf + size(correcti,1) / size(answer,1);
+        rank(i) = cumtestPerf;
+end
+
+%For FERET test, want probeID (answer), then rank, then matched ID no.,
+%then FA flag, then "matching score".  This will be a matrix with: 
+%probe  rank    match                           FAflag          matching score
+%i      1       trainClass(nearest(i,:))        Sdist(:,i)>4.7  1./Sdist(:,i)
+%i      2       OR rankmat(i,:)'
+%i      3
+%i      4
+
Index: proiecte/PPPP/ica/pcaFn.m
===================================================================
--- proiecte/PPPP/ica/pcaFn.m	(revision 94)
+++ proiecte/PPPP/ica/pcaFn.m	(revision 94)
@@ -0,0 +1,34 @@
+%function [U,R,E] = pcaFn(B);
+%
+%Principal components the normal way. Data in columns of B.
+%U is a matrix containing the principal component eigenvectors in 
+%     its columns.
+%R is a matrix containing the principal component representations in its
+%     rows. Each row is the representation for the corresponding column 
+%     of B. 
+%E is the vector of eigenvaules corresponding to the eigenvectors in U.
+
+function [U,R,E] = pcaFn(B);
+
+%Read data into columns of B;
+%B = datamat';
+[N,P] = size(B);
+
+%********subtract mean
+mb=mean(B');
+B=B-(ones(P,1)*mb)';
+
+%********Find eigenvectors vi of BB' (NxN)
+[V,D] = eig (1/(P-1)*(B*B'));
+
+%********Sort eigenvectors
+eigvalvec = max(D);
+[seigvals, index] = sort(eigvalvec); % sort goes low to high
+Vsort = V(:,[fliplr(index)]);
+
+U=Vsort;
+length = sqrt (sum (U.^2));
+U = U ./ (ones(N,1) * length);
+R = B'*U;
+E = fliplr(seigvals);
+
Index: proiecte/PPPP/ica/pcabigFn.m
===================================================================
--- proiecte/PPPP/ica/pcabigFn.m	(revision 94)
+++ proiecte/PPPP/ica/pcabigFn.m	(revision 94)
@@ -0,0 +1,55 @@
+%function [U,R,E] = pcabigFn(B);
+%Compute PCA by calculating smaller covariance matrix and reconstructing
+%eigenvectors of large cov matrix by linear combinations of original data
+%given by the eigenvecs of the smaller cov matrix. 
+%Data in Cols of B. Third version.  
+%
+%***** justification
+%
+%B = N x P data matrix.  N = dim of data  (Data in Cols of B, zero mean)
+%                        P = #examples
+%                        N >> P
+%
+%Want eigenvectors ui of BB' (NxN)
+%Solution:
+%Find eigenvectors vi of B'B (PxP)
+%From def of eigenvector: B'Bvi = di vi ---> BB'Bvi = di Bvi
+%Eigenvecs of BB' are Bvi
+%-------------------------------
+%[V,D] = eig (B'B)
+%Eigenvecs are in cols of V.    (Sorted cols)
+%
+%U = BV;  Cols of U are Bvi (eigenvecs of lg cov mat.) (Gave unit length)
+%R = B'U; Rows of R are pcarep of each observation.
+%E = eigenvalues        (eigenvals of small and large cov mats are equal)
+%*****
+
+function [U,R,E] = pcabigFn(B);
+
+%Read data into columns of B;
+%B = datamat';
+[N,P] = size(B);
+
+%********subtract mean
+mb=mean(B');
+B=B-(ones(P,1)*mb)';
+
+%********Find eigenvectors vi of B'B (PxP)
+[V,D] = eig (1/(P-1)*(B'*B));   %scale factor gives eigvals correct
+                                %magnitude for large cov mat 
+                                %(assuming sample cov)
+                                %(assuming sample cov)
+%********Sort eigenvectors
+eigvalvec = max(D);
+[seigvals, index] = sort(eigvalvec); % sort goes low to high
+Vsort = V(:,[fliplr(index)]);
+
+%********Reconstruct
+U = B*Vsort;  % Cols of U are Bvi. (N-dim Eigvecs)
+
+%********Give eigvecs unit length.  Improves classification.
+length = sqrt (sum (U.^2));
+U = U ./ (ones(N,1) * length);
+
+R = B'*U;  % Rows of R are pcarep of each image.
+E = fliplr(seigvals);
Index: proiecte/PPPP/ica/runica.m
===================================================================
--- proiecte/PPPP/ica/runica.m	(revision 94)
+++ proiecte/PPPP/ica/runica.m	(revision 94)
@@ -0,0 +1,38 @@
+% Copyright 1996 Tony Bell
+% This may be copied for personal or academic use.
+% For commercial use, please contact Tony Bell 
+% (tony@salk.edu) for a commercial license.
+
+
+% Script to run ICA on a matrix of images. Original by Tony Bell. 
+% Modified by Marian Stewart Bartlett.
+
+%Assumes image gravalues are in rows of x. Note x gets overwritten.
+%Will find N independent components, where N is the number of images.
+
+%There must be at least 5 times as many examples (cols of x) as the
+%dimension of the data (rows of x). 
+
+N=size(x,1); P=size(x,2); M=N;	 %M is dimension of the ICA output
+spherex;                          % remove first and second order stats from x
+xx=inv(wz)*x;                     % xx thus holds orig. data, w. mean extracted.
+
+%******** setup various variables
+w=eye(N); count=0; perm=randperm(P); sweep=0; Id=eye(M);
+oldw=w; olddelta=ones(1,N*M); angle=1000; change=1000;
+
+%******** Train. outputs a report every F presentations.
+% Watch "change" get small as it converges. Try annealing learning 
+% rate, L, downwards to 0.0001 towards end.
+% For large numbers of rows in x (e.g. 200), you need to use a low 
+% learning rate (I used 0.0005). Reduce if the output blows 
+% up and becomes NAN. If you have fewer rows, use 0.001 or larger.
+
+B=50; L=0.0005; F=5000; for I=1:1000, sep96; end; 
+B=50; L=0.0003; F=5000; for I=1:200, sep96; end; 
+B=50; L=0.0002; F=5000; for I=1:200, sep96; end;   
+B=50; L=0.0001; F=5000; for I=1:200, sep96; end; 
+
+%******** 
+uu=w*wz*xx;  % make separated output signals.
+cov(uu')     % check output covariance. Should approximate 3.5*I.
Index: proiecte/PPPP/ica/sep96.m
===================================================================
--- proiecte/PPPP/ica/sep96.m	(revision 94)
+++ proiecte/PPPP/ica/sep96.m	(revision 94)
@@ -0,0 +1,30 @@
+% sep96.m implements the learning rule described in Bell \& Sejnowski, Vision
+% Research, in press for 1997, that contained the natural gradient (w'w).
+%
+% Bell & Sejnowski hold the patent for this learning rule. 
+%
+% SEP goes once through the mixed signals, x 
+% (which is of length M), in batch blocks of size B, adjusting weights,
+% w, at the end of each block.
+% sepout is called every F counts.
+%
+% I suggest a learning rate (lrate) of 0.006, and a blocksize (B) of 
+% 300, at least for 2->2 separation.
+% When annealing to the right solution for 10->10, however, lrate of
+% less than 0.0001 and B of 10 were most successful.
+%
+% Copyright 1996 Tony Bell
+% This may be copied for personal or academic use.
+% For commercial use, please contact Tony Bell 
+% (tony@salk.edu) for a commercial license.
+
+x=x(:,perm);
+sweep=sweep+1; t=1;
+noblocks=fix(P/B);
+BI=B*Id;
+for t=t:B:t-1+noblocks*B,
+  count=count+B;
+  u=w*x(:,t:t+B-1); 
+  w=w+L*(BI+(1-2*(1./(1+exp(-u))))*u')*w;
+  if count>F, sepout; count=count-F; end;
+end;
Index: proiecte/PPPP/ica/sepout.m
===================================================================
--- proiecte/PPPP/ica/sepout.m	(revision 94)
+++ proiecte/PPPP/ica/sepout.m	(revision 94)
@@ -0,0 +1,5 @@
+% For reports on training progress
+
+[change,olddelta,angle]=wchange(oldw,w,olddelta); 
+oldw=w;
+fprintf('****sweep=%d, change=%.4f angle=%.1f deg., [N%d,M%d,P%d,B%d,L%.5f]\n',sweep,change,180*angle/pi,N,M,P,B,L);
Index: proiecte/PPPP/ica/spherex.m
===================================================================
--- proiecte/PPPP/ica/spherex.m	(revision 94)
+++ proiecte/PPPP/ica/spherex.m	(revision 94)
@@ -0,0 +1,11 @@
+% SPHEREX - spheres the training vector x.
+%    Requires x, P, to be predefined, and defines mx, c, wz.
+
+fprintf('\nsubtracting mean\n');
+mx=mean(x');
+x=x-(ones(P,1)*mx)';
+fprintf('calculating whitening filter\n');
+c=cov(x');
+wz=2*inv(sqrtm(c));
+fprintf('whitening\n');
+x=wz*x;
Index: proiecte/PPPP/ica/wchange.m
===================================================================
--- proiecte/PPPP/ica/wchange.m	(revision 94)
+++ proiecte/PPPP/ica/wchange.m	(revision 94)
@@ -0,0 +1,7 @@
+% Calculates stats on most recent weight change - magnitude and angle between
+% old and new weight. 
+
+function [change,delta,angle]=wchange(w,oldw,olddelta)
+  [M,N]=size(w); delta=reshape(oldw-w,1,M*N);
+  change=delta*delta';
+  angle=acos((delta*olddelta')/sqrt((delta*delta')*(olddelta*olddelta')));
Index: proiecte/PPPP/ica/zeroMn.m
===================================================================
--- proiecte/PPPP/ica/zeroMn.m	(revision 94)
+++ proiecte/PPPP/ica/zeroMn.m	(revision 94)
@@ -0,0 +1,11 @@
+%function Xzm = zeroMn(X)
+%Returns a zero-mean form of the matrix X. Each row of Xzm will have 
+%zero mean, same as in spherex.m. For PCA, put the observations in cols
+%before doing zeroMn(X).
+
+function Xzm = zeroMn(X);
+
+[N,P] = size(X);
+mx=mean(X');
+Xzm=X-(ones(P,1)*mx)';
+
