This document contains a tutorial on Matlab with a principal components analysis for a set of face images as the theme.
I wrote this tutorial while a graduate student in the Artificial Intelligence Laboratory of the Computer Science and Engineering Department at the University of California, San Diego. Now it's here at CSIM-AIT. You might still be able to find the original here.
This work is licensed under a
Creative
Commons Attribution 3.0 License. You are free to use the code here for
any purpose you like, but please acknowledge its source if you do. It's not
required, but if you make improvements I'd appreciate it if you sent me your
updates.
Here's an index to this tutorial:
I used the Ekman and Friesen Pictures of Facial Affect (POFA) for this tutorial. Unfortunately, the faces are copyrighted so if you want them you'll have to pay Paul Ekman some cash to get our online version of the database. Contact Gary Cottrell for details about our versions, or see http://www.paulekman.com/researchproducts.html on getting a CD directly from Ekman.
Don't despair, however. There are a bunch of nice face databases that are free. Check the Face Recognition Home Page for a list of what's available.
I won't go into it in detail here, but the idea is that face images can be economically represented by their projection onto a small number of basis images that are derived by finding the most significant eigenvectors of the pixelwise covariance matrix for a set of training images. A quick Google Search shows that a lot of people like to play with this technique. In my tutorial I simply show how to get some eigenfaces and play with them in Matlab.
If you don't use emacs, the Matlab emacs mode may be one reason to use it. To get Matlab mode working in emacs, put matlab.el somewhere emacs can find it and add something like the following to your .emacs file:
(autoload 'matlab-mode "matlab" "Enter Matlab mode." t)
(setq auto-mode-alist (cons '("\\.m$" . matlab-mode) auto-mode-alist))
(defun my-matlab-mode-hook ()
(setq matlab-indent-function t) ; if you want function bodies indented
(setq fill-column 76) ; where auto-fill should wrap
(matlab-mode-hilit)
(turn-on-auto-fill))
(setq matlab-mode-hook 'my-matlab-mode-hook)
(autoload 'matlab-shell "matlab" "Interactive Matlab mode." t)
(defun my-matlab-shell-mode-hook ()
'())
(setq matlab-mode-hook 'my-matlab-mode-hook)
>> help imread % Online help is useful...
IMREAD Read image from graphics file.
A = IMREAD(FILENAME,FMT) reads the image in FILENAME into
A. If the file contains a grayscale intensity image, A is
a two-dimensional array. If the file contains a truecolor
etc...
>> I = imread('test.png'); % Read image test.png into variable I
>> size(I) % Get the number of rows and columns of I
ans =
320 240
>> I(1:10,1:10) % Display row 1-10 and col 1-10 of I
ans =
26 20 20 22 22 20 19 20 22 22
26 24 19 22 22 19 21 24 24 22
25 24 22 24 21 19 20 22 24 22
25 21 21 24 22 20 21 22 24 20
24 20 22 25 22 19 20 21 25 24
20 19 22 22 24 21 21 24 22 24
20 20 21 25 25 24 24 25 20 20
20 20 21 24 26 25 25 25 20 20
20 19 22 24 24 25 24 24 22 21
21 17 21 25 21 22 24 25 22 21
>> colormap(gray(256)); % Use a 256-value grayscale color map
>> image(I); % Display I as an image
>> daspect([1 1 1]); % Set x:y aspect ratio to be 1:1
If you like the PGM format like I do, here are matlab functions for reading and writing them:
Matlab is really nice for linear algebra stuff and visualization, but sorta sucks when it comes to file I/O. Or rather, it's not much easier than C, although there are functions for reading and writing entire arrays to and from files. The following function, defined in load_images.m, is an example of how to read a bunch of images, make column vectors out of each of them, and return the result.
function [Images,w,h] = load_images(filelist,downscale_f)
%LOAD_IMAGES Load a set of images listed in a file.
%
% [IMGS,W,H] = LOAD_IMAGES(FILELIST) Treat each line of
% the file named FILELIST as the filename of a PGM image,
% and load each image as one column of the return array
% IMGS.
%
% LOAD_IMAGES(FILELIST,DOWNSCALE_F) Do the same thing,
% but downscale each image's width and height by a factor
% of DOWNSCALE_F. Useful if the images are too big to be
% loaded into memory all at once.
% Matthew Dailey 2000
% Check argument consistency
if nargin < 1 | nargin > 2
error('usage: load_images(filelist[,downscale_f]');
end;
if nargin == 1
downscale_f = 1.0;
end;
Images = []; old_w = 0; old_h = 0; w=0; h=0;
% Open input file
numimgs = linecount(filelist);
fid = fopen(filelist,'r');
if fid < 0 | numimgs < 1
error(['Cannot get list of images from file "' filelist, '"']);
end;
% Get the images
for i = 1:numimgs
imgname = fgetl(fid);
if ~isstr(imgname) % EOF is not a string
break; % Exit from loop on EOF
end;
fprintf(1,'loading PGM file %s\n',imgname);
Img = readpgm(imgname); % Read this image as a 2D array
if i==1 % If this is first image, figure things out
old_w = size(Img,2); % - like sizes of the downscaled images
old_h = size(Img,1);
if downscale_f <= 1.0
w = old_w; h = old_h;
else
w = round(old_w/downscale_f); h = round(old_h/downscale_f);
end;
Images = zeros(w*h,numimgs); % - preallocate size of the return matrix
end;
if downscale_f > 1.0
Img = im_resize(Img,w,h); % downscale using bicubic spline interp
end;
Images(1:w*h,i) = reshape(Img',w*h,1); % Make a column vector
end;
fclose(fid); % Close the filelist when done
fprintf(1,'Read %d images.\n',numimgs);
% The function returns the output arguments Images, w, and h here.
>> [Imgs,w,h] = load_images('facelist',4);
Here is how you would convert a column of Imgs back into an image and display it:
>> L = Imgs(:,10); >> L = reshape(L,w,h)'; % Reshapes column vector into a 2D array >> image(L); >> colormap(gray(256)); >> daspect([1 1 1]);
The function pc_evectors uses Turk and Pentland's trick to get the eigenvectors of A*A' from the eigenvectors of A'*A. It uses the function sortem to sort the eigenvectors and eigenvalues by eigenvalue. [Here is a faster version, sortem2, by James Javurek-Humig, which is considerably faster, if you have matlab version 7 or better.] To use pc_evectors, just do:
And to explore the eigenvalue spectrum and how much variance the first n vectors account for, try the following:>> [Vecs,Vals,Psi] = pc_evectors(Imgs,30); % Get top 30 PC evectors of Imgs
>> plot(Vals); % To plot the eigenvalues
>> CVals = zeros(1,length(Vals)); % Allocate a vector same length as Vals
>> CVals(1) = Vals(1);
>> for i = 2:length(Vals) % Accumulate the eigenvalue sum
CVals(i) = CVals(i-1) + Vals(i);
end;
>> CVals = CVals / sum(Vals); % Normalize total sum to 1.0
>> plot(CVals); % Plot the cumulative sum
>> ylim([0 1]); % Set Y-axis limits to be 0-1
One thing that can be fun is to try to figure out what the top few eigenvectors encode. You can make a scatter plot of one component against another as follows:
>> Proj = Vecs(:,1:2)' * Imgs; % Project each image onto first 2 evectors
>> size(Proj)
ans =
2 110
>> plot(Proj(1,:),Proj(2,:),'r.') % To get scatterplot of PC 1 vs PC 2
>> Labels = textread('labels','%s'); % Read labels from file 'labels'
>> text(Proj(1,:),Proj(2,:),Labels); % Add text labels at each plotted point
Finally, you might be interested in determining the reconstruction error involved in representing an image by its projection onto a few eigenvectors. Here is how you would project onto and reconstruct from eigenvectors 1-10.>> pc_ev_1 = Vecs(:,1); % Get PC eigenvector 1 >> pc_ev_1 = reshape(pc_ev_1,w,h)'; % Reshape into 2D array >> imagesc(pc_ev_1); % Display as image scaled 0-255
>> OrigImg = Imgs(:,20); % Grab image 20 >> Projection = Vecs(:,1:10)'*(OrigImg - Psi); % Project onto ev's 1-10 >> Reconstruction = Vecs(:,1:10)*Projection+Psi; % Reconstruct from projection >> Reconstruction = reshape(Reconstruction,w,h)'; >> image(Reconstruction); >> colormap(gray(256)); >> daspect([1 1 1]);